Hi, Dan. I was able validate that your code correctly implements the Implicit Neumann Boundary Condition (in 1-D).
Here is a link to the plot of the solution after a transient simulation for 0.2 seconds. The plot on the left is generated in Mathematica, and that to the right is that generated by the FiPy’s matplotlib viewer class. The two are identical. https://imperialcollegelondon.box.com/s/lb8v1iqxb3rqhxsphn9cmhwcmh5jzpiz Can we add this as a feature request in the project’s github page, perhaps for a future release? Does Implicit Neumann BCs occur often in other disciplines ? Once again, thanks for your help in providing this solution approach. Krishna From: fipy-boun...@nist.gov [mailto:fipy-boun...@nist.gov] On Behalf Of Daniel Wheeler Sent: 10 June 2016 18:18 To: Multiple recipients of list <fipy@nist.gov> 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 <smit...@mit.edu<mailto:smit...@mit.edu>> 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<mailto: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 [cid:image001.png@01D1C700.5100A820] 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:image002.png@01D1C700.5100A820] 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:image003.png@01D1C700.5100A820])is [cid:image004.png@01D1C700.5100A820] As you can see, the linearised Boundary condition, is cast in terms of the field variable, [cid:image003.png@01D1C700.5100A820] . 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:image005.png@01D1C700.5100A820] 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: fipy-boun...@nist.gov<mailto:fipy-boun...@nist.gov> [mailto:fipy-boun...@nist.gov<mailto:fipy-boun...@nist.gov>] On Behalf Of Raymond Smith Sent: 09 June 2016 16:06 To: fipy@nist.gov<mailto: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<mailto: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. [cid:image006.png@01D1C700.5100A820] I think we need to solve for the PDE, keeping this implicit BC, rather than immediately evaluating the term [cid:image007.png@01D1C700.5100A820] , since [cid:image008.png@01D1C700.5100A820] 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:image009.png@01D1C700.5100A820] 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> [mailto:fipy-boun...@nist.gov<mailto:fipy-boun...@nist.gov>] On Behalf Of Raymond Smith Sent: 09 June 2016 14:23 To: fipy@nist.gov<mailto: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<mailto: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> [mailto:fipy-boun...@nist.gov<mailto:fipy-boun...@nist.gov>] On Behalf Of Gopalakrishnan, Krishnakumar Sent: 08 June 2016 23:42 To: fipy@nist.gov<mailto: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> [mailto:fipy-boun...@nist.gov] On Behalf Of Raymond Smith Sent: 08 June 2016 23:36 To: fipy@nist.gov<mailto: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<mailto: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<mailto: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<mailto: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<mailto: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<mailto: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<mailto:fipy@nist.gov> http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ] -- Daniel Wheeler
_______________________________________________ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
