Re: casting implicit Boundary Conditions in FiPy
Hi Krishina and Ray, Thanks for the interesting discussion. I'm not 100% sure about everything that Krishina is asking for in the latter part of the discussion so I'm just going to address the code that Ray has developed below (my code is below). I think there is a way to handle the right hand side boundary condition both implicitly and with second order accuracy using an ImplicitSourceTerm boundary condition trick. Sort of similar to what is here, http://www.ctcms.nist.gov/fipy/documentation/USAGE.html#applying-fixed-flux-boundary-conditions . Like I said, I think the code below is both second order accurate (for fixed dx) and implicit. Extending this to 2D might raise a few issues. The grid spacing needs to be in the coefficient so, obviously, dx needs to be fixed in both x and y directions. Also, I'm not sure if I'm missing a factor of dx in the source term in 1D as I'm using both a divergence and an ImplicitSourceTerm so there is some question about volumes and face areas in 2D as well. I'm also confused about the signs, I had to flip the sign in front of the source to make it work. It seems to look right though in 1D and you can just take one massive time step to get to the answer. This will only work if k is negative otherwise it's unstable, right? The way I came up with the source was the following n.grad(phi) = k * (phi_P + n.grad(phi) * dx / 2) and then solve for n.grad(phi) which gives n.grad(phi) = k * phi_P / (1 - dx * k / 2) where phi_P is the value of phi at the cell next to the boundary with the implicit boundary condition. Then we fake the outbound flux to be the expression on the right of the equation. Just as a general note it would be great in FiPy if we could come up with a nice way to write boundary conditions in scripts that did all these tricks implicitly without having to know all these background details about FV and how FiPy works. import fipy as fp nx = 50 dx = 1. mesh = fp.Grid1D(nx=nx, dx=dx) phi = fp.CellVariable(name="field variable", mesh=mesh, value=1.0) D = 1. k = -1. diffCoeff = fp.FaceVariable(mesh=mesh, value=D) diffCoeff.constrain(0., mesh.facesRight) valueLeft = 0.0 phi.constrain(valueLeft, mesh.facesLeft) #phi.faceGrad.constrain([phi], mesh.facesRight) # This is the problematic BC #phi.faceGrad.constrain(phi.harmonicFaceValue, mesh.facesRight) # This is the problematic BC #phi.faceGrad.constrain([k * phi.harmonicFaceValue], mesh.facesRight) # This is the problematic BC implicitCoeff = -D * k / (1. - k * dx / 2.) eq = (fp.TransientTerm() == fp.DiffusionTerm(diffCoeff) - \ fp.ImplicitSourceTerm((mesh.faceNormals * implicitCoeff * mesh.facesRight).divergence)) timeStep = 0.9 * dx**2 / (2 * D) timeStep = 10.0 steps = 800 viewer = fp.Viewer(vars=phi, datamax=1., datamin=0.) for step in range(steps): eq.solve(var=phi, dt=timeStep) viewer.plot() On Thu, Jun 9, 2016 at 12:02 PM, Raymond Smith wrote: > Hi, Krishna. > > Perhaps I'm misunderstanding something, but I'm still not convinced the > second version you suggested -- c.faceGrad.constrain([-(j_at_c_star + > partial_j_at_op_point*(c.faceValue - c_star))], mesh.facesTop) -- isn't > working like you want. Could you look at the example I suggested to see if > that behaves differently than you expect? > > Here's the code I used. To me it looks very similar in form to c > constraint above and at first glance it seems to behave exactly like we > want -- that is, throughout the time stepping, n*grad(phi) is constrained > to the value, -phi at the surface. Correct me if I'm wrong, but my > impression is that this is the behavior you desire. > > from fipy import * > > nx = 50 > dx = 1. > mesh = Grid1D(nx=nx, dx=dx) > > phi = CellVariable(name="field variable", mesh=mesh, value=1.0) > D = 1. > > valueLeft = 0.0 > phi.constrain(valueLeft, mesh.facesLeft) > #phi.faceGrad.constrain([phi], mesh.facesRight) # This is the problematic > BC > #phi.faceGrad.constrain(phi.harmonicFaceValue, mesh.facesRight) # This is > the problematic BC > phi.faceGrad.constrain([-phi.harmonicFaceValue], mesh.facesRight) # This > is the problematic BC > > eq = TransientTerm() == DiffusionTerm(coeff=D) > > timeStep = 0.9 * dx**2 / (2 * D) > steps = 800 > > viewer = Viewer(vars=phi, datamax=1., datamin=0.) > > for step in range(steps): > eq.solve(var=phi, dt=timeStep) > viewer.plot() > > Cheers, > Ray > > On Thu, Jun 9, 2016 at 11:44 AM, Gopalakrishnan, Krishnakumar < > k.gopalakrishna...@imperial.ac.uk> wrote: > >> Hi Ray, >> >> >> >> Yes. You make a good point. I see that the analytical solution to the >> particular problem I have posted is also zero. >> >> >> >> The reason I posted is because I wanted to present an (oversimplified) >> analogous problem when posting to the group, retaining the generality, >> since many other subject experts might have faced similar situations. >> >> >> >> The actual problem I am solving is the solid diffusion PDE (only 1 >> equation) in a Li-ion battery. I am
Re: casting implicit Boundary Conditions in FiPy
Okay, I see your dilemma more clearly now; thanks for clarifying and yes, I was missing a few things, so this makes more sense now. Anyway, at this point, we're stepping beyond of my FiPy-fu, so I'll be curious to learn as well. Best, Ray On Fri, Jun 10, 2016 at 7:59 AM, Gopalakrishnan, Krishnakumar < k.gopalakrishna...@imperial.ac.uk> wrote: > Hi Ray, > > > > Thanks a lot for your further inputs. > > > > I realise that perhaps I have not been articulating the problem correctly, > and there is a miscommunication. > > > > In your code, you have the k3 term as constant. More importantly, the > boundary constraint code is set outside the loop , i.e. before beginning > time-stepping. > > > > In this case, a single time-step difference will not make a difference, > i.e . an explicit use of ‘j’ from previous time-step will not be too bad > (we’ll have to think a bit more about stability, due to the time-delay > introduced) > > > > However, my problem is quite different. The linearization of my > (non-linear) RHS of the boundary condition is performed around the > operating point c*. The operating point changes at every time-step. Thus, > the line of code implementing the BC constraint must be inside the loop. > The problem now becomes the following: > > > > 1. The code print c.faceValue return a numerical result. That means > that this expression is evaluated numerically, rather than being left as > part of an implicit exprerssion. The only way a numerical solution can be > obtained is to simply return the existing solution, i.e. c.faceValue > simply uses solution values from the previous time-step. > > 2. After computing the solution, the c* point also needs to be > changed to the new operating point, i.e. c* = c. (You always linearise > around your latest operating point – principle of small-signal perturbation > and linearization) > > 3. Thus, in the next iteration and all subsequent time-steps, > c.faceValue – c* term evaluates to zero again. Thus, we are solving the > PDE only with the DC component of the Taylor series expansion (i.e. > effectively ignoring the first-order derivatives too !). > > > > Note that this sharply contrasts with the problem with the implicit PDE > set-up, i.e. one could have an implicit source term, which is a linear > function of solution variable, and fiPy has a built-in method to handle > this in terms of ImplicitSourceTerm(coeff=…). You simply obtain the > linearised mathematical expression of the source term in paper, and feed it > to the coeff term in the ImplicitSourceTerm method. c.faceValue must not > get evaluated numerically, and must be used implicitly as part of obtaining > the solution at this time-step. > > > > I guess describing the problem in words are falling a bit short, since the > idea is complex enough. I have drawn a simple flow-chart of the workflow > for solving this Implicit BC problem. Can you perhaps kindly take a look > at it ? > https://onedrive.live.com/redir?resid=618ADC32B6C2B152!19186&authkey=!AHrlNkAcvENEZ0M&v=3&ithint=photo%2cjpg > > > > Meanwhile, may I also ask if other Fipy users or developers had to deal > with non-linear Implicit Neumann boundary conditions in their problems? > > > > > > Krishna > > > > *From:* fipy-boun...@nist.gov [mailto:fipy-boun...@nist.gov] *On Behalf > Of *Raymond Smith > *Sent:* 10 June 2016 00:16 > > *To:* fipy@nist.gov > *Subject:* Re: casting implicit Boundary Conditions in FiPy > > > > Hi, Krisnha. > > I don't know the details about how FiPy handles the boundary conditions, > so I can't confirm whether or not it's solving it with the current > (implicit) or previous (explicit) time step's value when constraining the > gradient. I think that will have to come from one of the devs or someone > more involved than I. However, I did note that in the simple example code > below, I put it in the same form as you have with the extra terms, > n*grad(c) = k1 + k2*(c - k3) > where k1, k2, and k3 are all constants. > It also seems to work fine. It may be treated explicitly, but it does at > least seem to update in time correctly. Is there a significant problem with > it being explicit (using previous time step value) anyway? The overall > order of accuracy should be the same as implicit (current time step), > although I guess stability could be worse? > > In other words, for the first time step, the (c - cstar) term _should_ be > very close to zero (exactly zero when treated explicitly), because c is > initially equal to cstar at the surface. If your time step is large enough > that c at the surface deviates "significantly" from cstar within a single > time step, then your time steps are probably too large anyway. Of course, > the (c - cstar) should become non-zero at least by the second time step as > long as k1 is finite. > > from fipy import * > > nx = 50 > dx = 1. > mesh = Grid1D(nx=nx, dx=dx) > > phi = CellVariable(name="field variable", mesh=mesh, value=1.0) > D = 1. > k1 = -0.1 > k
