I had something similar happen with an Eigenvalue problem I was working
on. Two things affected the stability of my convergence. The first was
which solver I was using, turns out that the PCG solvers worked best for
me. The second was limiting the number of iterations in each sweep. I am
no expert in nonlinear diffusion but these may help.
Best of luck,
-mike waters
On 10/27/15 9:11 AM, Dan Brunner wrote:
I'm working on a nonlinear diffusion problem (specifically finding how
long a spatially-dependent heat flux can be held on a surface before
it reaches the melting point, 'simple' example at the end) and I am
seeking assistance in setting up the sweeping. I've looked through all
of the FiPy examples to learn what I can. If I set it up to sweep
until a certain residual threshold is reached, each iteration on the
time loop takes longer and longer because the initial residual keeps
increasing and the residual decrease between sweeps keeps decreasing.
This usually goes on until the residual bounces up and down in a sweep
loop. If I try some logic to decrease the sweep dt, the sweep
residuals blow up.
Does anyone have guidance on what to do here and/or can refer me to a
resource where I may learn how to do this better?
Thanks
~Dan
from fipy import *
import numpy as np
q0=5e8
qLambda=1e-2
Lx=qLambda*4.0
dx=qLambda/20.0
Ly=0.02
dy=Ly/2.0
T0=300.0
peakTemp=3695.0
resMin=1e-2
mesh=Gmsh2D('''
Lx=%(Lx)g;
Ly=%(Ly)g;
dx=%(dx)g;
dy=%(dy)g;
Point(1) = { 0, 0, 0, dx};
Point(2) = {Lx, 0, 0, dx};
Point(3) = {Lx, Ly, 0, dy};
Point(4) = { 0, Ly, 0, dy};
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 1};
Line Loop(5) = {1, 2, 3, 4};
Plane Surface(6) = {5};
''' % locals())
X, Y = mesh.faceCenters
time = Variable()
temperature = CellVariable(name = "temperature",
mesh = mesh,
value = T0,
hasOld=1)
# thermal properties of tungsten
rho = 19300.0
cp = (5.33e+01*temperature**0.0+
4.50e-01*temperature**1.0-
8.30e-04*temperature**2.0+
7.07e-07*temperature**3.0-
2.73e-10*temperature**4.0+
3.91e-14*temperature**5.0)
kappa = (2.32e+02*temperature**0.0-
2.67e-01*temperature**1.0+
2.42e-04*temperature**2.0-
1.13e-07*temperature**3.0+
2.50e-11*temperature**4.0-
1.99e-15*temperature**5.0)
eq = TransientTerm(coeff=rho*cp) == DiffusionTerm(coeff=kappa)
dTime = 0.9*dx**2.0/(2.0*5e-5)
if __name__ == '__main__':
viewer = Viewer(vars=temperature, datamin=0., datamax=4000.)
viewer.plot()
while temperature().max() < peakTemp:
time.setValue(time()+dTime)
temperature.updateOld()
facesDecayBottom = (mesh.facesBottom & (X>=(qLambda)))
valueDecayBottom =
-q0*np.exp(-(X-qLambda)/qLambda)/kappa.getFaceValue()
temperature.faceGrad.constrain(0.0, mesh.facesBottom)
temperature.faceGrad.constrain(valueDecayBottom, facesDecayBottom)
temperature.faceGrad.constrain(0.0, mesh.facesTop)
res=1e9
while res>resMin:
res = eq.sweep(var=temperature,dt=dTime)
print res
print time, temperature().max()
if __name__ == '__main__':
viewer.plot()
_______________________________________________
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 ]
