Hi Jonathan,
First of all, many thanks once again for your timely and helpful response.
> Can you post a minimal code that exhibits the divergence?
I include the minimal code at the end of this message.
>> I'm know it has a solution, because it's been
>> found before by other authors.
>
> Maybe they're all wrong. 8^)
Hmm... skepticism, love it!
> I have a feeling that this is the problem. The FixedFlux boundary
> condition really needs an overhaul. It works reliably for zero flux, but
> doesn't always do the right thing for non-zero values with both convection
> and diffusion.
I can't tell you how relieved is was to hear that.
> You should be able to achieve the same effect by using a source.
> Daniel recently wrote a FAQ for this in the documentation for trunk,
Ok, didn't notice that, I'll have a closer look.
> but in FiPy 1.2 syntax, I *think* you would write:
>
> boundarySource = VectorFaceVariable(mesh=...)
> boundarySource[mesh.getExteriorFaces()] = -15.0 eq1 =
> TransientTerm(coeff=xi**2.0*c) \
> - ExponentialConvectionTerm(coeff = unitVector*xi_f**3.0/R*dR_dt*c_f,
> diffusionTerm=diffTerm) \ + 3.0*c*x_p*xi**2.0*dR_dt/R \
> == diffTerm \
> + boundarySource.getDivergence()
I'll give it a try, but should I adopt what you suggest here and carry it
on to 2.0 or use more traditional boundary conditions?
> There's no need to pass any boundaryConditions to sweep() if you do it
> this way.
Aha, just out of curiosity, I am I right when I say that - also in this
case - when I don't pass any BCs in sweep(), fipy just uses its default
Neumann BC equal to zero and the boundary source compensates for that?
> It converges for me with both 1.2 and trunk, but I could easily not be
> in the same parameter space as you.
I double tested the script I am posting in an initially empty parameter
space (or namespace as IPython calls it) and it exhibited exactly the
behavior I was talking about in my original post. I hope you can spot some
faulty lines.
> We are trying very hard to push a new release out, hopefully in
> another day or two (just waiting on a clean build of the documentation). I
> would suggest you wait for that. At that point, we will be encouraging
> everyone to migrate to the new version.
You guys are on fire! Understood, I'll wait for the new version and
meanwhile, I'll experiment with what you sent me.
Best regards!
Etienne
******************
Begin minimal code
******************
import fipy
import sys
from fipy import *
from numpy import pi
from pdb import set_trace as bp
from pylab import close, ion, plot, figure
ion()
# Mesh
nx = 40
dx = 1.0/40.0
dt = 0.045/500.0
mesh = Grid1D(nx=nx, dx=dx)
# Variable
x_p = CellVariable(mesh=mesh, value=0.3386, hasOld=1)
# Coefficients
c = CellVariable(mesh=mesh, value=7326.0)
c_f = c.getArithmeticFaceValue()
R = 6.6e-5
dR_dt = -3067.8126
D_p = 3.92e-9
unitVector = VectorFaceVariable(mesh=mesh, value=1.0)
xi = CellVariable(mesh=mesh, value=mesh.getCellCenters()[...,0])
xi_f = xi.getArithmeticFaceValue()
# Boundary conditions
bc_xp = (
FixedFlux(faces=mesh.getFacesRight(), value=-15.0),
FixedFlux(faces=mesh.getFacesLeft(), value=0.0)
)
# Equation
diffTerm = ImplicitDiffusionTerm(coeff=xi**2.0*c*D_p/R**2.0)
eq1 = TransientTerm(coeff=xi**2.0*c) \
- ExponentialConvectionTerm(coeff = unitVector*xi_f**3.0/R*dR_dt*c_f,
diffusionTerm=diffTerm) \
+ 3.0*c*x_p*xi**2.0*dR_dt/R \
== diffTerm
close('all')
bp()
res = 4000000000.0
figure()
while res>100000.0:
plot(xi.value, x_p.value)
res = eq1.sweep(
var = x_p,
boundaryConditions=bc_xp,
solver=LinearLUSolver(tolerance=1e-11),
dt=dt)
print res
****************
End minimal code
****************
