Yes, indeed the residual is converged and I'm not sweeping well past convergence. Even if I do one sweep per time step, the time it takes to do each sweep increases with each time step, independent of the time step size. Could this be a memory issue?
~Dan On Mon, Dec 21, 2015 at 10:42 AM, Guyer, Jonathan E. Dr. < [email protected]> wrote: > 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 ] > >
_______________________________________________ fipy mailing list [email protected] http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
