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 ]
