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 <
k.gopalakrishna...@imperial.ac.uk> 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:* fipy-boun...@nist.gov [mailto:fipy-boun...@nist.gov] *On Behalf
> Of *Raymond Smith
> *Sent:* 09 June 2016 17:03
> *To:* fipy@nist.gov
> *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 <
> 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 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:image001.png@01D1C290.47669780]
>
>
>
> 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:image002.png@01D1C290.47669780] 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:image003.png@01D1C290.47669780])is
>
> [image: cid:image004.png@01D1C290.47669780]
>
>
>
> As you can see, the linearised Boundary condition, is cast in terms of the
> field variable, [image: cid:image003.png@01D1C290.47669780]. 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:image005.png@01D1C290.47669780]
>
>
>
> 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:* fipy-boun...@nist.gov [mailto:fipy-boun...@nist.gov] *On Behalf
> Of *Raymond Smith
> *Sent:* 09 June 2016 16:06
>
>
> *To:* fipy@nist.gov
> *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 <
> k.gopalakrishna...@imperial.ac.uk> 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:image006.png@01D1C290.47669780]
>
>
>
> I think we need to solve for the PDE, keeping this implicit BC, rather
> than immediately evaluating the term [image:
> cid:image007.png@01D1C290.47669780], since [image:
> cid:image008.png@01D1C290.47669780]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:image009.png@01D1C290.47669780] 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:* fipy-boun...@nist.gov [mailto:fipy-boun...@nist.gov] *On Behalf
> Of *Raymond Smith
> *Sent:* 09 June 2016 14:23
>
>
> *To:* fipy@nist.gov
> *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 <
> k.gopalakrishna...@imperial.ac.uk> 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:* fipy-boun...@nist.gov [mailto:fipy-boun...@nist.gov] *On Behalf
> Of *Gopalakrishnan, Krishnakumar
> *Sent:* 08 June 2016 23:42
> *To:* fipy@nist.gov
> *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:* fipy-boun...@nist.gov [mailto:fipy-boun...@nist.gov
> <fipy-boun...@nist.gov>] *On Behalf Of *Raymond Smith
> *Sent:* 08 June 2016 23:36
> *To:* fipy@nist.gov
> *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 <
> k.gopalakrishna...@imperial.ac.uk> 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
> fipy@nist.gov
> http://www.ctcms.nist.gov/fipy
> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
>
>
>
>
> _______________________________________________
> fipy mailing list
> fipy@nist.gov
> http://www.ctcms.nist.gov/fipy
> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
>
>
>
>
> _______________________________________________
> fipy mailing list
> fipy@nist.gov
> http://www.ctcms.nist.gov/fipy
> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
>
>
>
>
> _______________________________________________
> fipy mailing list
> fipy@nist.gov
> http://www.ctcms.nist.gov/fipy
> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
>
>
>
> _______________________________________________
> fipy mailing list
> fipy@nist.gov
> http://www.ctcms.nist.gov/fipy
> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
>
>
_______________________________________________
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
