Dan - I recommend you plot the residual as a function of sweeps for large numbers of sweeps at each time step. Are you converged? Are you continuing to sweep well after you converged? You might find that a smaller time step will allow you to converge with a smaller number of sweeps. While the diffusion equation is unconditionally stable, accuracy deteriorates at larger time steps, especially with your non-linear coefficients.
On Dec 17, 2015, at 9:33 AM, Dan Brunner <[email protected]> wrote: > I'd like to revisit this. Now I've progressed to doing simulations where the > heat flux profile is swept back and forth on the surface. This requires both > small time steps to get the details of the movement of the heat flux profile > as well as a long total time for the simulation. Doing this I noticed that as > the simulation progresses in time, each successive time step takes longer. > And no matter what I cange in the simulation it does not change this fact. > I've narrowed it down to the 'res = eq.sweep(var=temperature,dt=dTime)' part. > Each time this while loop is called, the time it takes to execute increases > by a constant factor. It starts from ~0.5s and increases by ~0.01s with each > sweep, quickly growing to levels too slow to both simulating. Any guidence on > this would be appreciated. > > ~Dan > > On Sat, Oct 31, 2015 at 2:04 PM, Guyer, Jonathan E. Dr. > <[email protected]> wrote: > Glad to hear it. > > On Oct 31, 2015, at 8:45 AM, Dan Brunner <[email protected]> wrote: > > > That did it. Thanks! > > > > ~Dan > > > > On Thu, Oct 29, 2015 at 10:18 AM, Dan Brunner <[email protected]> wrote: > > Thanks for the help Jon. I did not realize that the coeff in the > > TransientTerm was assumed to only have a time dependence. I agree that > > derivations of the heat equation are often not careful about non-linear > > coefficients. I'll try your reformulation of it later this week. We have > > some other tools in-house to benchmark it with. > > > > ~Dan > > > > On Wed, Oct 28, 2015 at 5:08 PM, Guyer, Jonathan E. Dr. > > <[email protected]> wrote: > > Dan - > > > > Interesting problem. What I actually find is that the residual decays or at > > least holds fairly steady until some point at which the residual rises > > rapidly, regardless of the number of sweeps. > > > > I was ultimately able to make some headway by looking at the form of your > > equations. By writing > > > > TransientTerm(coeff=rho*cp) > > > > you are saying that cp (and rho) have time dependence, but not spatial > > dependence. Clearly, in your case, the do have spatial dependence. In > > theory, sweeping should resolve some of this, but is apparently inadequate. > > > > Some derivations lump everything together as > > > > alpha = kappa/(rho*cp) > > > > and solve > > > > TransientTerm(coeff=1.) == DiffusionTerm(coeff=alpha) > > > > but I don't think this is right. > > > > The derivations of the heat equation I've seen that are careful about > > non-uniform coefficients (which is not most derivations of the heat > > equation that I've seen) end up with > > > > \rho c_p \frac{\partial T}{\partial t} = \nabla\cdot(\kappa\nabla T) > > > > This can be rearranged as > > > > \frac{\partial T}{\partial t} = \nabla\cdot(\kappa/(\rho c_p)\nabla T) > > - \nabla(1/(\rho c_p))\cdot\kappa\nabla T > > > > which FiPy can express as > > > > eq = (fp.TransientTerm(coeff=1.) > > == fp.DiffusionTerm(coeff=kappa/(rho*cp)) > > - (1/(rho*cp)).grad.dot(kappa*temperature.grad)) > > > > > > I don't know if it's right, but it solves. > > > > You also should be able to increase the time step considerably. > > > > - Jon > > > > > > > > On Oct 27, 2015, at 9:11 AM, Dan Brunner <[email protected]> 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 ] > > > > > > > > _______________________________________________ > > 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 ] _______________________________________________ fipy mailing list [email protected] http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
