Thanks for the advice Mike. None of the solvers have solved this problem. If I limit the number of iterations in each sweep, I would be left with an ever increasing residual as the problem evolved.
~Dan On Tue, Oct 27, 2015 at 12:07 PM, Michael Waters <[email protected]> wrote: > 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 [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 ]
