Hm, changing the right boundary condition to c + dc/dx = 0 may make a somewhat more reasonable output.
Going back to my original question, when I change your script to solve that problem, the solution at least looks _reasonable_ to me, although again, worth checking against an analytical solution. That is, the solution value at the right decreases and the slope at the right seems to decrease proportionally. On Thu, Jun 9, 2016 at 11:06 AM, Raymond Smith <[email protected]> wrote: > 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. >> >> >> >> >> >> I think we need to solve for the PDE, keeping this implicit BC, rather >> than immediately evaluating the term , since 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, 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 ]
