Re: 1D diffusion problem
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
> 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
>> 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 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 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)
>
Re: 1D diffusion problem
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
> 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 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 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)
Re: 1D diffusion problem
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 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 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 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
fipy@
