1D diffusion problem

[Clara Maurel]

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...")

smime.p7s
Description: S/MIME cryptographic signature

_______________________________________________
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]