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 ]
