Hi, Kyle.
I've pasted my transient script here
<https://gist.github.com/raybsmith/34d49404d0117ffd74ff> and attached the
file. The script shows how I selected a different solver
<http://www.ctcms.nist.gov/fipy/fipy/generated/fipy.solvers.html> algorithm
and how I did sweeps in the transient steps, much like the "uncoupled"
version of the coupled diffusion example
<http://www.ctcms.nist.gov/fipy/examples/diffusion/generated/examples.diffusion.coupled.html#module-examples.diffusion.coupled>.
That said, I wouldn't trust the way I managed the time step size in that
file. It works for the specific case I solved, but the way to adjust the
time step can be very problem specific, and I don't think my time step
halving is even doing the right thing in that example. I haven't updated
that problem much because it was a toy example problem to see if FiPy might
be an appropriate tool for the types of problems I was working on
(unfortunately, the answer to that became to do my own spatial
discretization so I could use a more sophisticated time stepping algorithm).
Hope this helps,
Ray
On Fri, Dec 5, 2014 at 9:26 PM, Kyle Briton Lawlor <[email protected]>
wrote:
> Hi, Ray.
>
> I am wondering if you might be willing to share your code with me. I am
> curious to see how you handled sweeping, with the transient terms of your
> equations for the electrolyte problem. I am also interested to see how you
> switch the GMRESSolver and if this might help me with my problem. I’ll
> leave it to you to decide if you want to share this is the mailing list, if
> it all.
>
> Best,
> Kyle
import fipy as fp
Lx = 1.
nx = 1500
dx = Lx/nx
mesh = fp.Grid1D(nx=nx, dx=dx)
phi = fp.CellVariable(name="elec. potential", mesh=mesh, value=0.,
hasOld=True)
c_m = fp.CellVariable(name="conc.", mesh=mesh, value=1.0, hasOld=True)
V = -2.5e-0
c0 = 1.
c_m.setValue(c0)
D_m = 2.03e-9 # m^2/s
D_p = 1.33e-9 # m^2/s
D_p = D_p/D_m
D_m = 1.
eq1 = (0 == fp.TransientTerm(coeff=1., var=c_m)
+ fp.DiffusionTerm(coeff=-D_m, var=c_m)
+ fp.DiffusionTerm(coeff=D_m*c_m.faceValue, var=phi))
eq2 = (0 == 0.
+ fp.DiffusionTerm(coeff=-(D_p - D_m), var=c_m)
+ fp.DiffusionTerm(coeff=-(D_p + D_m)*c_m.faceValue, var=phi))
phi.constrain(0., mesh.facesLeft)
phi.constrain(V, mesh.facesRight)
c_m.constrain(c0, mesh.facesLeft)
c_m.faceGrad.constrain(c_m.faceValue*phi.faceGrad, mesh.facesRight)
Xcc = mesh.cellCenters[0]
J = 1 - fp.numerix.exp(V)
c_m_analyt = fp.CellVariable(name="analyt. conc.", mesh=mesh)
c_m_analyt.setValue(1 - J*Xcc)
phi_analyt = fp.CellVariable(name="analyt. phi", mesh=mesh)
phi_analyt.setValue(fp.numerix.log(1-J*Xcc))
eq = eq1 & eq2
viewer = fp.Viewer(vars=(c_m, phi))
dt = 5e-4
tf = 8.
t = 0.
success = True
nItMax = 200
while t < tf:
if success:
c_m.updateOld()
phi.updateOld()
res = 1e10
nIt = 0
while res > 5e-2 and nIt < nItMax:
res = eq.sweep(dt=dt,
solver=fp.LinearGMRESSolver(),
)
# Print every 10th sweep
if not nIt % 10:
print "Iteration:", nIt, "Residual:", res
nIt += 1
if nIt < 5:
dt *= 1.01
success == True
elif nIt > 150:
dt *= 0.5
success = False
t += dt
print "\n time: {t}, dt: {dt}\n".format(t=t, dt=dt)
viewer.plot()
anstol = 1e-3
print c_m.allclose(c_m_analyt, rtol=anstol, atol=anstol)
print phi.allclose(phi_analyt, rtol=anstol, atol=anstol)
_______________________________________________
fipy mailing list
[email protected]
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
