Hi all, I have a question regarding a zero flux boundary condition in my model. First, I’ll input the model equations and the relevant boundary conditions, as well as the domain:
As you can see, I have Robin boundary conditions at \partial e^-, which I enter into the weak form after applying the divergence theorem. In particular, I have a zero-flux boundary condition. However, when I evaluate \boldsymbol{N}_2|_{\partial e^-}\cdot\hat{\boldsymbol{e}}_1 along boundary \partial e^-, I do not obtain exactly zero; instead the value of flux oscillates between very small (~1e-9) positive and negative values, which increase in size at later time steps (although they remain very small). Why might this be happening? Is it simply due to numerical error? If so, how can I improve my code such that the zero-flux boundary condition is better obeyed? Would a more refined mesh help? Please let me know if you require more information to answer my question. Many thanks.
Can you upload ( not paste into message ) your code or a minimal
working example? I can’t be sure of getting to it but am looking at
drift-diffusion and electrochemistry somewhat.
Hi Mike, I understand if you’re busy. If you have the time, here’s my code: forForum2.edp (21.6 KB)
Note, I have implemented a Newton method to deal with the non-linearity; I put the equations into weak form and then expanded about an approximate solution. So, the code solves for a correction which is then appended to the approximation. Note also that I’ve implemented a varying time step (the solution evolves more rapidly towards the start and so initially requires a very small time step, which is then increased at later times). I started using freefem a few months ago, so my code could no doubt be improved; one thing I’d like to try and do is use the ‘varf’ keyword rather than ‘problem’ as that would allow me to try out some different solvers.
I can upload the rest of the boundary conditions if it would help, just say if you’d like me to.
I was looking at some etching and used IIRC SNESsolve for the initial
solution and had some remaining flux but don’t remember where I left it.
I could upload that but I use several idp files I made and would have to
look for the code at this point.
Okay, only look for it if you’ve got time. At the very least, any gut instinct ideas that you might have about why the flux isn’t completely zero would be very useful. I can only imagine that it’s simply a result of some kind of numerical error. I should mention that the aspect ratio of my adapted mesh is poor (i.e. close to the boundaries, the triangles are extremely long and thin, which is something that I’m going to ask about correcting in a separate post), which could also have something to do with it. What do you think?
It is due to numerical error. Your Neumann BC can only ever be as accurate as your discrete approximation of the derivatives. You can use a finer mesh or higher order elements if this is a problem, but I don’t expect errors of the order of 1e-9 to be an issue…
Is that an absolute number or relative to something? In actual fact their
will be a flux of reactants getting to the surface and typically
engaging in something like A+ + e- → A which then creates a concentation
gradient and net flux away from the boundary.
Hi Mike. Species 2 in the above model is not involved in any reactions at boundary \partial e^-. The flux is absolute and should be equal to zero. But from Chris’ answer, it seems to be due to numerical error, which makes sense.
Hi Chris, this makes a lot of sense. I’ll investigate the effect of increasing the mesh fineness. When you say higher order elements, are you referring to P1, P2 etc.? Apologies, I’m quite rusty with the theory behind FEM.
Yes, that’s exactly right. In general, you can reduce numerical error by h-refinement (finer mesh) or p-refinement (higher order functional approximation). Note that your estimate of how accurately the Neumann condition is satisfied is also dependent on the functional approximation you make via your discretization. An error of 1e-9 is likely just numerical noise, as your actual errors (the difference between the continuous solution and the discrete one) are probably much higher unless your discretization is extremely resolved.
Okay, this is very useful. One final question: do you think the code’s failure to better meet the zero-flux BC is partly due to the poor aspect ratio of my adapted mesh (as seen in my other question that you answered today)?
I still haven’t looked at your code but in general there can be “high” fields
at the interface and the balance between diffusion ( gradient of carriers)
and drift( gradient of voltage) may be an issue of local mesh parameters.
It may be interesting to see if the “noise” correlates with any particular
local mesh features or if a more uniform mesh ( squares with diagonals)
changes the noise distirbution.
But still there is not really much of an “absolute” nu,mber notion.
1e-9 could be a nanoamp per something which may be important .
Probably you do have two large components for a drift term
and a diffusion term but important current amounts may be
due to a “small difference” between these two. Does the noise
effect anything you care about?
Hi Mike, I agree, it would certainly be interesting to see if the noise changes with respect to the mesh structure. It is hard to tell if the noise is affecting the solution significantly; I really only have a small number of contour plots in the literature (and a few voltage-time plots) to go by. I’ll look at trying to find a way to reduce the noise, which will then allow me to see if it has any affect upon the solution.
I just plot drift and diffusion separately. I guess also you can do this just in self consistent
iterations between charge density and voltage and then think about the issues
with consistency between the Lagrangian elements and the exponential dependence
on voltage - this is why I’m looking at other interpolation schemes. If you have a linear
potential and are near equ you may be better off interpolating the carrier density
with an exponential. I pushed that too while trying to get datascope out and
btw I would not discount the utility of 1D models for a starting point.
