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 < [email protected]> 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:* [email protected] [mailto:[email protected]] *On Behalf > Of *Raymond Smith > *Sent:* 10 June 2016 00:16 > > *To:* [email protected] > *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 > k2 = 0.1 > phistar = 1.0 > > 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 > phi.faceGrad.constrain([k1 + k2*(phistar - phi.harmonicFaceValue)], > mesh.facesRight) > > 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 6:32 PM, Gopalakrishnan, Krishnakumar < > [email protected]> wrote: > > Hi Ray, > > > > Thank you for replying. Yes, for this problem solved by your code, it > doesn’t have any other terms inside the [] in > phi.faceGrad.constrain([-phi.harmonicFaceValue], > mesh.facesRight). > > > > What is really confusing me is that, the variable to solve for is phi, > right ? And, we need the solution for this particular time-step within > the loop. But the phi.harmonicFaceValue method performs a numerical > evaluation , i.e. substituting the value of field variable (solution) from > the previous time-step. > > > > If this is indeed the case, there is a problem, because unlike the minimal > example which works correctly, the real problem has the term (c.faceValue > - c_star) . Upon seeing the c.faceValue term, fiPy automatically > substitutes the numerical value of the field variable, ie. boundary value > from the previous time-step. However, this is the same as c_star, which > I methodically verified. Thus, the two terms gets cancelled out, > effectively eliminating the first-order derivative term from the > taylor-series expansion of the non-linear j, thereby leading to > first-order truncation errors. > > > > Is my understanding of how the whole process works correct ? What I am > saying is, there should be a way to cast the (linearised) solution > variable term in the BC, in an implicit format (without substituting the > numerical value from the previous time-step), alike the > ImplicitSourceTerm(coeff=…) style of casting. > > > > > > The solution variable ‘c’ must be obtained by solving the PDE, the > left-boundary no-flux BC , and this right boundary *Implicit* BC, at this > time-step, ie. without substituting c.facevalue (numerical values from > previous time-step) > > > > If we do this c.faceValue substitution, we are simply left with > > > > > > i.e. evaluation of the j at a given DC operating point, introducing first > order errors in the local time-step. > > > > Am I on the wrong track here ? > > > > > > Krishna > > > > *From:* [email protected] [mailto:[email protected]] *On Behalf > Of *Raymond Smith > *Sent:* 09 June 2016 17:03 > *To:* [email protected] > *Subject:* Re: casting implicit Boundary Conditions in FiPy > > > > 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 < > [email protected]> 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 solving this PDE in a pseudo-2D > domain. i.e. I have defined a Cartesian 2D space, wherein the y-coordinate > corresponds to the radial direction. The bottom face corresponds to > particle centres, and the top face corresponding to surface of each > spherical particle. The x-axis co-ordinate corresponds to particles along > the width (or thickness) of the positive electrode domain. Diffusion of Li > is restricted to be within the solid particle (i.e. y-direction only), by > defining a suitable tensor diffusion coefficient as described in the > Anisotropic diffusion example and FAQ in FiPy. I have normalised my x and > y dimensions to have a length of unity. > > > > Now, the boundary condition along the top face is > > [image: cid:[email protected]] > > > > Now, j is non-linear (Butler-Volmer), and I am using a Taylor-expanded > linear version for this boundary condition. All other field variables[image: > cid:[email protected]] are assumed as constants. The idea > is to set up the infrastructure and solve this problem independently, > before worrying upon the rubrics of setting up the coupled system. In a > similar fashion, I have built up and solved the solid phase potential PDE > (thanks to your help for pointing out about the implicit source term). > Thus, the idea is to build up the coupled P2D Newman model piecemeal. > > > > The linearised version of my BC’s RHS at a given operating point ([image: > cid:[email protected]])is > > [image: cid:[email protected]] > > > > As you can see, the linearised Boundary condition, is cast in terms of the > field variable, [image: cid:[email protected]]. Hence, we > need it in an implicit form corresponding to (pseudocode: > c.faceGrad.constrain([-( (j_at_c_star - partial_j_at_op_point*c_star) + > coeff = partial_j_at_op_point)],mesh.facesTop) , or something of this > form/meaning. (just like the very useful ImplicitSourceterm method) > > > > If I instead apply the c.faceValue method,i.e. using it in setting the BC > as > > > > c.faceGrad.constrain([-(j_at_c_star + partial_j_at_op_point*(c.faceValue - > c_star))], mesh.facesTop), then c.faceValue gets immediately evaluated > at the operating point, c_star, and we are left with 0 multiplying the > first-order derivative. > > > > ie. the Boundary conditions becomes, > > [image: cid:[email protected]] > > > > Leading to huge loss of accuracy. > > > > Is there any hope at all in this situation ? J . Cheers and thanks for > your help thus far. > > > > > > Krishna > > > > *From:* [email protected] [mailto:[email protected]] *On Behalf > Of *Raymond Smith > *Sent:* 09 June 2016 16:06 > > > *To:* [email protected] > *Subject:* Re: casting implicit Boundary Conditions in FiPy > > > > Hi, Krishna. > > Could you give a bit more detail and/or an example about how you know it's > doing the wrong thing throughout the solution process? In the example you > sent, the correct solution is the same (c(x, t) = 0) whether you set > n*grad(phi) to zero or to phi at the boundary, so it's not a good example > for concluding that it's not behaving as you'd expect. It's helpful here to > find a situation in which you know that analytical solution to confirm one > way or the other. For example, you should be able to get the solution to > the following problem using a Fourier series expansion: > > dc/dt = Laplacian(c) > > c(t=0) = 1 > > x=0: c = 0 > > x=1: c - dc/dx = 0 > > Ray > > > > On Thu, Jun 9, 2016 at 10:52 AM, Gopalakrishnan, Krishnakumar < > [email protected]> wrote: > > Hi Ray, > > > > Thanks for your help. > > > > But when I apply phi.harmonicFaceValue , it is immediately evaluated to a > numerical result (a zero vector in this case, since initial value = 0, the > data type is > fipy.variables.harmonicCellToFaceVariable._HarmonicCellToFaceVariable', > i.e. the boundary condition is not remaining implicit. > > > > The same is the case with the examples.convection.robin example. Here, > the phi.faceValue method is used. However, this also results in an > immediate numerical evaluation. > > > > However, what is actually required is that, the BC must remain implicit > (in variable form, without getting numerically evaluated), being cast in > terms of the field variable being solved for. Then the solver needs to > solve the PDE on the domain to yield the solution of the field variable. > > > > [image: cid:[email protected]] > > > > I think we need to solve for the PDE, keeping this implicit BC, rather > than immediately evaluating the term [image: > cid:[email protected]], since [image: > cid:[email protected]]is the field variable to be solved > for, i.e. there ought to be some way to cast the Boundary condition as > implicit. > > > > In FiPy, I have previously set up an implicit source term, [image: > cid:[email protected]] by using the following code > snippet, ImplicitSourceTerm(coeff=k) . Perhaps there might be an > equivalent method in FiPy to set up the implicit BC, I think ? > > > > > > Krishna > > > > > > > > > > *From:* [email protected] [mailto:[email protected]] *On Behalf > Of *Raymond Smith > *Sent:* 09 June 2016 14:23 > > > *To:* [email protected] > *Subject:* Re: casting implicit Boundary Conditions in FiPy > > > > Oh, right, the boundary condition is applied on a face, so you need the > facevalue of phi: > > phi.faceGrad.constrain([phi.harmonicFaceValue]) > > Ray > > > > On Thu, Jun 9, 2016 at 7:28 AM, Gopalakrishnan, Krishnakumar < > [email protected]> wrote: > > Hi ray, > > > > Casting the implicit PDE does not work for my problem. FiPy throws up a > ton of errors. > > I am attaching a minimal example (based off example1.mesh.1D) > > > > *from *fipy *import ** > > nx = 50 > dx = 1. > mesh = Grid1D(nx=nx, dx=dx) > > phi = CellVariable(name=*"field variable"*, mesh=mesh, value=0.0) > D = 1. > > valueLeft = 0.0 > phi.constrain(valueLeft, mesh.facesLeft) > phi.faceGrad.constrain([phi], mesh.facesRight) > > *# This is the problematic BC *eq = TransientTerm() == DiffusionTerm(coeff > =D) > > timeStep = 0.9 * dx**2 / (2 * D) > steps = 100 > > viewer = Viewer(vars=phi) > > *for *step *in *range(steps): > eq.solve(var=phi, dt=timeStep) > viewer.plot() > > > > The errors are as follows: > > > > line 22, in <module> > > eq.solve(var=phi, dt=timeStep) > > \fipy\terms\term.py", line 211, in solve > > solver = self._prepareLinearSystem(var, solver, boundaryConditions, > dt) > > \fipy\terms\term.py", line 170, in _prepareLinearSystem > > buildExplicitIfOther=self._buildExplcitIfOther) > > \fipy\terms\binaryTerm.py", line 68, in _buildAndAddMatrices > > buildExplicitIfOther=buildExplicitIfOther) > > \fipy\terms\unaryTerm.py", line 99, in _buildAndAddMatrices > > diffusionGeomCoeff=diffusionGeomCoeff) > > \fipy\terms\abstractDiffusionTerm.py", line 337, in _buildMatrix > > nthCoeffFaceGrad = coeff[numerix.newaxis] * > var.faceGrad[:,numerix.newaxis] > > \fipy\variables\variable.py", line 1575, in __getitem__ > > unit=self.unit, > > \fipy\variables\variable.py", line 255, in _getUnit > > return self._extractUnit(self.value) > > \fipy\variables\variable.py", line 561, in _getValue > > value[..., mask] = numerix.array(constraint.value)[..., mask] > > IndexError: index 50 is out of bounds for axis 1 with size 50 > > > > I have tried including the implicit BC within the time-stepper loop, but > that does not still help. > > > > > > Best Regards > > > > Krishna > > > > > > > > *From:* [email protected] [mailto:[email protected]] *On Behalf > Of *Gopalakrishnan, Krishnakumar > *Sent:* 08 June 2016 23:42 > *To:* [email protected] > *Subject:* RE: casting implicit Boundary Conditions in FiPy > > > > Hi Raymond, > > > > Sorry, it was a typo. > > > > Yes, It is indeed d (phi)/dx, the spatial derivative BC. I shall try > setting phi.faceGrad.constrain([k*phi], mesh.facesRight), and see if it > will work. > > > > Thanks for pointing this out. > > > > > > Krishna > > > > *From:* [email protected] [mailto:[email protected] > <[email protected]>] *On Behalf Of *Raymond Smith > *Sent:* 08 June 2016 23:36 > *To:* [email protected] > *Subject:* Re: casting implicit Boundary Conditions in FiPy > > > > Hi, Krishna. > > Just to make sure, do you mean that the boundary condition is a derivative > with respect to the spatial variable or with respect to time as-written? If > you mean spatial, such that d\phi/dx = k*phi, have you tried > > phi.faceGrad.constrain(k*phi) and that didn't work? > > If you mean that its value is prescribed by its rate of change, then I'm > not sure the best way to do it. Could you maybe do it explicitly? > - Store the values from the last time step with hasOld set to True in the > creation of the cell variable > - In each time step, calculate the backward-Euler time derivative > manually and then set the value of phi with the phi.constrain method. > > > > Ray > > > > On Wed, Jun 8, 2016 at 6:26 PM, Gopalakrishnan, Krishnakumar < > [email protected]> wrote: > > I am trying to solve the standard fickean diffusion equation on a 1D > uniform mesh in (0,1) > > > > $$\frac{\partial \phi}{\partial t} = \nabla.(D \nabla\phi)$$ > > > > with a suitable initial value for $\phi(x,t)$. > > > > The problem is that, one of my boundary conditions is *implicit*, i.e. is > a function of the field variable being solved for. > > > > $ \frac{\partial\phi}{\partial t} = k \phi $ , at the right boundary edge, > k = constant > > > > The left BC is not a problem, it is just a standard no-flux BC. > > > > How do I cast this *implicit BC* in FiPy ? Any help/pointers will be much > appreciated. > > > > > > Best regards > > > > Krishna > > Imperial College London > > > _______________________________________________ > fipy mailing list > [email protected] > http://www.ctcms.nist.gov/fipy > [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ] > > > > > _______________________________________________ > fipy mailing list > [email protected] > http://www.ctcms.nist.gov/fipy > [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ] > > > > > _______________________________________________ > fipy mailing list > [email protected] > http://www.ctcms.nist.gov/fipy > [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ] > > > > > _______________________________________________ > fipy mailing list > [email protected] > http://www.ctcms.nist.gov/fipy > [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ] > > > > > _______________________________________________ > fipy mailing list > [email protected] > http://www.ctcms.nist.gov/fipy > [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ] > > > > _______________________________________________ > fipy mailing list > [email protected] > http://www.ctcms.nist.gov/fipy > [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ] > >
_______________________________________________ fipy mailing list [email protected] http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
