Hi Jon, Thank you so much for taking the time to explain everything!
Clara > On 22 Apr 2019, at 13:31, Guyer, Jonathan E. Dr. (Fed) via fipy > <[email protected]> wrote: > > Also, the natural boundary condition for FiPy is no-flux, so there's no need > to do > > T_var.faceGrad.constrain(0, where=mesh.facesLeft) > >> On Apr 22, 2019, at 1:28 PM, Guyer, Jonathan E. Dr. (Fed) via fipy >> <[email protected]> wrote: >> >> Clara - >> >> - The cause of everything dropping to zero is that T_var.old is initialized >> to zero. This isn't a particularly useful initial value, but that's what it >> is. Usually, T_var.updateOld() is the first thing called in the time step >> loop, so this initial condition doesn't matter, but because of the way >> you've structured your adaptive time steps, the first sweep is done with an >> erroneous initial condition. >> >> Call T_var.updateOld() one time before you start taking time steps. >> >> - Your sweep tolerance is too severe. I adjusted to 1.0e-6 and the adaptive >> stepper works much better. >> >> - The way you define and update dt, you'll get recursion errors eventually. >> Instead, define >> >> dt = 0.9*dr**2/(2.0*Diff[0].value) >> >> If you're curious, this is because: >> >> * Diff is a FaceVariable >> * so, Diff[0] is a UnaryOperator Variable >> * so, dt = 0.9*dr**2/(2.0*Diff[0]) is a BinaryOperator Variable >> * so, dt *= 1.1/dt *= 0.8 is a BinaryOperator Variable that gets nested >> deeper >> and deeper every time step >> >> Diff[0].value short circuits that >> >> - Your problem is completely linear, so sweeps aren't necessary at this point >> >> - This equation is fully implicit, so there's no need to start with such a >> tiny initial time step >> >> - Jon >> >>> On Apr 19, 2019, at 2:42 PM, Clara Maurel <[email protected]> wrote: >>> >>> Hello everyone, >>> >>> I apologize in advance for this very basic question… >>> I am trying to solve a 1D diffusion problem where I have a cooling planet >>> (core+mantle) with an initial temperature profile and a diffusion >>> coefficient that takes different values in the core and in the mantle. I >>> want a fixed boundary condition at the surface (T = 200 K) and a zero >>> gradient at the center of the body, but the temperature at the center (left >>> of the mesh) is not fixed and should cool. >>> When I run the script, the whole profile immediately goes down to zero >>> (except the end right of the mesh) and diffuses from there. I must be doing >>> something wrong with the left boundary condition but I can’t figure out >>> what. >>> >>> Could someone help me to find out what I missed? Any hint would be greatly >>> appreciated! >>> >>> Thank you in advance for your help! >>> Clara >>> >>> Here is my simple script: >>> >>> from datetime import datetime >>> startTime = datetime.now() >>> import numpy as np >>> import scipy >>> from fipy import * >>> from fipy.solvers.pysparse import LinearLUSolver as Solver >>> from sympy import * >>> from scipy import interpolate >>> >>> def Get_closest_id(L,value): >>> return list(L).index(min(L, key=lambda x:abs(x-value))) >>> >>> #------------------------------------# >>> # DIFFUSION EQUATION: INITIALIZATION # >>> #------------------------------------# >>> >>> ## Dimensions (km) >>> R_body = 200.0 >>> R_core = 50.0 >>> >>> ## Temperatures (K) >>> T_core = 1200.0 >>> T_surf = 200.0 >>> >>> ## Mesh >>> nr = 1000 >>> dr = R_body / nr >>> mesh = Grid1D(nx=nr, dx=dr) >>> >>> ## Variable >>> T_var = CellVariable(name="T", mesh=mesh, hasOld=True) >>> >>> ## Initial temperature profile >>> slope = (T_core-T_surf)/(R_core-R_body) >>> intercept = T_core-slope*R_core >>> T_var[:] = np.array([T_core]*int(R_core/dr) + >>> list(slope*np.arange(R_core+dr, R_body+dr, dr)+intercept)) >>> >>> ## Initial conditions (fixed surface value, null gradient at the center of >>> the core) >>> T_var.constrain(T_surf, mesh.facesRight) >>> T_var.faceGrad.constrain(0, where=mesh.facesLeft) >>> >>> ## Diffusion coefficient (different in core and mantle, in km2 My-1) >>> Diff = FaceVariable(mesh=mesh) >>> X = mesh.faceCenters[0] >>> Diff.setValue(150.0, where=(R_core >= X)) >>> Diff.setValue(15.0, where=(R_core < X) & (R_body >= X)) >>> >>> ## Equation >>> eq = TransientTerm(coeff=1.) == DiffusionTerm(Diff) >>> >>> if __name__ == "__main__": >>> viewer = Viewer(vars=(T_var),ymin=0,ymax=1400) >>> viewer.plot() >>> >>> if __name__ == '__main__': >>> raw_input("Press <return> to proceed...") >>> >>> ## Time and integration parameters (time in My; 3.154e13 = 1 My in s) >>> time = 0.0 >>> duration = 1000.0 >>> dt = 0.9*dr**2/(2.0*Diff[0]) >>> >>> total_sweeps = 4 >>> tolerance = 1.0e-10 >>> >>> solver = Solver() >>> >>> #-------------------------------# >>> # DIFFUSION EQUATION: MAIN LOOP # >>> #-------------------------------# >>> >>> while time < duration: >>> >>> res0 = eq.sweep(T_var, dt=dt, solver=solver) >>> >>> for sweeps in range(total_sweeps): >>> res = eq.sweep(T_var, dt=dt, solver=solver) >>> >>> if res < res0 * tolerance: >>> time += dt >>> dt *= 1.1 >>> T_var.updateOld() >>> print dt, res >>> if __name__ == '__main__': >>> viewer.plot() >>> >>> else: >>> dt *= 0.8 >>> T_var[:] = T_var.old >>> print "warning " + str(res) >>> >>> print datetime.now() - startTime >>> >>> if __name__ == '__main__': >>> raw_input("Press <return> to proceed...") >>> >>> >>> _______________________________________________ >>> 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 ]
smime.p7s
Description: S/MIME cryptographic signature
_______________________________________________ fipy mailing list [email protected] http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
