Thank you for the note. I found the following to work:
------------------------------------------------------------------------------
#defining solution variable
phi = CellVariable(name="solution variable",
mesh=mesh,
value=0.)
flux = FaceVariable(name='constant flux',
mesh=mesh,
value=1.)
#defining boundary conditions
phi.constrain(0.,mesh.facesRight)
#defining governing equations
termf = (mesh.facesLeft * flux).divergence
eqn = TransientTerm() == DiffusionTerm(coeff=diffCoeff) - D*termf
----------------------------------------------------------------------------------------------
For others reference: trying to extend this to the variable case by defining
something like
flux = (phi**2/(1+phi).faceGrad
does not work. But I found from "examples.convection.robin" an example of a
constraint which is dependent on the potential, which is (contrary to what I
said earlier) done by the constrain function, as in:
phi.faceGrad.constrain(-D*(phi.faceValue)**2/(1+phi.faceValue),mesh.facesLeft)
________________________________________
From: [email protected] [[email protected]] on behalf of Guyer,
Jonathan E. Dr. (Fed) [[email protected]]
Sent: Friday, June 08, 2018 2:56 PM
To: FIPY
Subject: Re: Flux BC Dependent on Potential
David -
The problem is here:
diffCoeff=Variable(value=1.)
diffCoeff.constrain(0., mesh.facesLeft)
It doesn't mean anything to constrain a Variable at particular faces. A
Variable doesn't exist on a mesh; it's just a value.
You can either do
diffCoeff=CellVariable(mesh=mesh, value=1.)
diffCoeff.constrain(0., mesh.facesLeft)
or
diffCoeff=Variable(value=1.)
and leave off the constraint.
Taking the divergence 51 faces to get 50 cells is not a problem.
DiffusionTerm is declared right here on the page you linked:
https://www.ctcms.nist.gov/fipy/fipy/generated/fipy.terms.html#fipy.terms.diffusionTerm.DiffusionTerm
DiffusionTerm is an implicit term, which is what you usually want.
ExplicitDiffusionTerm is generally not what you want.
> On Jun 4, 2018, at 11:13 PM, Ollodart, David Bernd <[email protected]>
> wrote:
>
> Hello,
>
> I have a reaction-diffusion equation which has a flux boundary condition
> dependent on the potential value at the boundary. I have seen in
> examples.diffusion.mesh1d how to set the boundary conditions to a constant or
> function of independent variables, but trying to insert the dependent
> variable, or its face value (evaluated at some face(s)) directly into
> .constrain does not work. I have seen at
> https://na01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fwww.mail-archive.com%2Ffipy%40nist.gov%2Fmsg00389.html&data=02%7C01%7Cjonathan.guyer%40nist.gov%7Cb48df7e9f52d4b681f4608d5ca9277f3%7C2ab5d82fd8fa4797a93e054655c61dec%7C1%7C1%7C636637652788127952&sdata=5f05NC56AzhaVB%2BGKg9gncSQTlDSsUQugUfNQzYNRDs%3D&reserved=0
> a question exactly on this topic but it uses FixedFlux which is now
> deprecated: https://www.ctcms.nist.gov/fipy/documentation/USAGE.html. But I
> have problems trying to implement the new methods suggested under "Applying
> fixed flux boundary conditions". By the way,!
there seems to be a typo; it should read ExplicitDiffusionTerm(diffCoeff)
rather than DiffusionCoeff(diffCoeff), there being no DiffusionCoeff() in the
fipy.terms package:
https://www.ctcms.nist.gov/fipy/fipy/generated/fipy.terms.html.
>
> I am first trying to use a constant external flux in the new method, which
> could be achieved by potential.faceGrad.constrain(1., mesh.facesLeft), then
> extend it to a potential-dependent flux. I have modified
> examples.diffusion.mesh1d as below. I define the exterior flux on the whole
> interior mesh, because the product with mesh.facesLeft makes it only non-zero
> on the left face. The script is:
> ---------------------------------------------------------------------------------------
> from fipy import *
> from scipy.special import erf
>
> #defining mesh of independent variables
> nx = 50
> dx = 1.
> mesh = Grid1D(nx=nx, dx=dx)
> timeStepDuration = 0.9*dx**2/(2*1)
> steps = 10
> x = mesh.cellCenters[0]
> t = timeStepDuration * steps
>
> #setting parameter values
> diffCoeff=Variable(value=1.)
> diffCoeff.constrain(0., mesh.facesLeft)
>
> #defining solution variable and known analytical
> phi = CellVariable(name="solution variable",
> mesh=mesh,
> value=0.)
>
> #defining boundary conditions
> phi.constrain(1.,mesh.facesRight)
> exteriorFlux = FaceVariable(name="fixed flux",
> mesh=mesh,
> value=1.,
> rank=1)
>
> #defining governing equations (might be matrix or vector valued for multiple
> dependent variables)
> eqn = TransientTerm() == ExplicitDiffusionTerm(coeff=diffCoeff) +
> (mesh.facesLeft * exteriorFlux).divergence
>
> #defining the viewer
> viewer = Viewer(vars=phi,
> datamin=0.,datamax=1.)
>
> #iterating through solution
> for step in range(steps):
> eqn.solve(var=phi,dt=timeStepDuration)
> viewer.plot()
>
> raw_input("Diffusion profile. Press <return> to proceed...")
> --------------------------------------------------------------------------------------
> However, this throws an IndexError: too many indices for array. If i print
> the product of exteriorFlux and mesh.facesLeft, it is an array of 51
> elements, the leftmost being 1, the others being 0, which is 1 greater than
> the size of the array found from mesh.cellCenters which is what the
> CellVariable is defined on. This error results if one specifies flux on one
> or both sides of the mesh, the latter being the example provided. Since there
> is 1 more face than there are cells, this result might be expected for
> exteriorFlux, but then which mesh should be chosen? Trying to define a
> different mesh throws an error for the product mesh.facesLeft * exteriorFlux.
>
> Thank you,
>
> David Ollodart
> _______________________________________________
> 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 ]
