Dan sent me a script to demonstrate what he was seeing. The steadily increasing solve time was because he was applying new constraints at each time step to achieve boundary conditions that evolved with time and position.
The problem with this is that FiPy does not replace constraints when new ones are applied, rather all constraints that have ever been applied get calculated and applied in succession, leading to the increasing solve time. It's possible to flush out constraints before applying new ones, but the better solution, which I sent to Dan, is to define the constraint once as a Variable expression that has the appropriate dependencies on position and time. On Dec 23, 2015, at 10:55 AM, Dan Brunner <[email protected]> wrote: > 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 ] _______________________________________________ fipy mailing list [email protected] http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
