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
> <FIPY@nist.gov> 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 <cmau...@mit.edu> 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
>> fipy@nist.gov
>> http://www.ctcms.nist.gov/fipy
>> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
>
>
> _______________________________________________
> fipy mailing list
> fipy@nist.gov
> http://www.ctcms.nist.gov/fipy
> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
_______________________________________________
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
