Hi
Thank you for the solution. This indeed looks like a nifty way to solve the
issue at hand.
Being a novice in FiPy, and Finite Volume methods in general, I am a bit fuzzy
about a couple of things.
1. The backward Euler method is used in formulating the equation,
"n.grad(phi) = k * (phi_P + n.grad(phi) * dx / 2) " . But Backward Euler is a
first order method, isn't it ? I am a bit confused about the second order
accuracy statement.
2. In, the equation, "eq = (fp.TransientTerm() ==
fp.DiffusionTerm(diffCoeff) - \
fp.ImplicitSourceTerm((mesh.faceNormals * implicitCoeff *
mesh.facesRight).divergence)) " ,
from what I understand the 'ImplicitCoeff' is zeroed out in all other
locations other than the right boundary by the boolean method 'mesh.facesRight'
, right ?
If this is being implemented, do we have the restriction to use fixed dx ?
Although dx appears in the Implicit term, it can just be set to the value of
'dx' of the last node (i.e. length of the last segment in a variable-mesh sized
1D geometry file). Am I missing something here ?
3. Also, I am a bit confused about the negative signs used in the equation
above. The implicitCoeff is given a negative sign, and the equation's
ImplicitSourceTerm is also given a negative sign. From what I undertand,
this is what we intend to do ?
* Assign the diffusion coefficient D throughout the domain, using the
faceVariable method.
* Manually set it to zero at the right boundary
* Then **add** the Diffusion achieved by the div.(D grad( phi)) by using
the ImplicitSourceTerm method, and ensure that this effect gets added only at
the right boundary.
If this is the workflow, then perhaps we just have to declare the
ImplicitSourceTerm to be positive, and use the addition sign instead of
subtraction sign for the implicitSourceTerm in the equation ?
4. Yes, the k needs to be negative for the analytical solution to be stable.
My apologies. I should have been more clear in my earlier posts to the list.
5. Taking one massive step to get to the answer. Here is a big question
(due to my poor knowledge in the subject)
I understand that the implicitness enables us to take these large time-steps
due to the affine and A-stable properties of Implicit Methods. But are we
required to do this ? My problem is that, I am solving a coupled system of
PDEs and my other equations, require me to step very very slowly, due to the
inherent stiffness of my system. The PDEs describe processes of varying
time-constants. So, does taking the timestep to be small inhibit us in
any way ?
Once again, thanks a lot for your solution. Based on the discusssions with Ray
and from this particular method from your posting to the list, I think I
shall be able to solve the problem in 1D. The issues you pointed out about
using this approach in 2D went right over my head. I have to think about this,
when I implement them perhaps in a week or so.
Thanks once again,
Krishna
________________________________
From: [email protected] <[email protected]> on behalf of Daniel Wheeler
<[email protected]>
Sent: Friday, June 10, 2016 6:17 PM
To: Multiple recipients of list
Subject: 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
<[email protected]<mailto:[email protected]>> 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
<[email protected]<mailto:[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
[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[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
([cid:[email protected]])is
[cid:[email protected]]
As you can see, the linearised Boundary condition, is cast in terms of the
field variable, [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,
[cid:[email protected]]
Leading to huge loss of accuracy.
Is there any hope at all in this situation ? :) . Cheers and thanks for your
help thus far.
Krishna
From: [email protected]<mailto:[email protected]>
[mailto:[email protected]<mailto:[email protected]>] On Behalf Of
Raymond Smith
Sent: 09 June 2016 16:06
To: [email protected]<mailto:[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]<mailto:[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.
[cid:[email protected]]
I think we need to solve for the PDE, keeping this implicit BC, rather than
immediately evaluating the term [cid:[email protected]] , since
[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,
[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]>
[mailto:[email protected]<mailto:[email protected]>] On Behalf Of
Raymond Smith
Sent: 09 June 2016 14:23
To: [email protected]<mailto:[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]<mailto:[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]>
[mailto:[email protected]<mailto:[email protected]>] On Behalf Of
Gopalakrishnan, Krishnakumar
Sent: 08 June 2016 23:42
To: [email protected]<mailto:[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]>
[mailto:[email protected]] On Behalf Of Raymond Smith
Sent: 08 June 2016 23:36
To: [email protected]<mailto:[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]<mailto:[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]<mailto:[email protected]>
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
_______________________________________________
fipy mailing list
[email protected]<mailto:[email protected]>
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
_______________________________________________
fipy mailing list
[email protected]<mailto:[email protected]>
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
_______________________________________________
fipy mailing list
[email protected]<mailto:[email protected]>
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
_______________________________________________
fipy mailing list
[email protected]<mailto:[email protected]>
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
--
Daniel Wheeler
_______________________________________________
fipy mailing list
[email protected]
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
