Re: Spline interpolation and fipy variable

[Daniel Wheeler]

On Fri, Jan 12, 2018 at 12:42 PM, Clara Maurel  wrote:
> Hi Daniel,
>
> Thank you very much! The code runs for me as well: fantastic!
>
> I have to admit that the distinction between cell and face variables remains
> a bit obscure to me…

The domain is divided into cells or boxes in 2D. The faces are where
the boxes meet and the cell centers are located at the center of the
boxes. FiPy calculates/resolves fluxes between the boxes.

> The cell variable X_var is the composition of the
> system I am modelling. Then "G" in d2G is the Gibbs free energy (called f
> here
> https://www.ctcms.nist.gov/fipy/examples/cahnHilliard/generated/examples.cahnHilliard.mesh2DCoupled.html#module-examples.cahnHilliard.mesh2DCoupled),
> which is itself a function of the composition X_var (the spline function).
> I don’t know if these info are more useful to figure out if d2G should be
> cell or face variable? As I understood it, because I want to solve for
> X_var, this should be a cell variable and nothing else.

X_var is a cell variable but the diffusion coefficient can either be a
cell variable or a face variable in FiPy. FiPy converts the diffusion
coefficient to a face variable if it's defined as a cell variable in
the input file. FiPy calculates the flux between cells so the
diffusion coefficient is only used at the faces where the flux is
calculated. d2G has values which are defined over the domain. How do
those values correspond to a location in that domain? Why does d2G
have a shape of (1001,) when there are only 1000 cells?

> By higher dimension do you also mean 2D? I know 2D requires a lot of
> changes, but I thought 2D required only a change of grid. In any case, 2D is
> not my top priority at the moment.

Ok, this is only being used for 1D right now so you're fine.

> Thank you again very much for your help!

Good luck with it.

-- 
Daniel Wheeler

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

Re: Spline interpolation and fipy variable

[Clara Maurel]

Hi Daniel,

Thank you very much! The code runs for me as well: fantastic!

I have to admit that the distinction between cell and face variables remains a 
bit obscure to me… The cell variable X_var is the composition of the system I 
am modelling. Then "G" in d2G is the Gibbs free energy (called f here 
https://www.ctcms.nist.gov/fipy/examples/cahnHilliard/generated/examples.cahnHilliard.mesh2DCoupled.html#module-examples.cahnHilliard.mesh2DCoupled
 
),
 which is itself a function of the composition X_var (the spline function).
I don’t know if these info are more useful to figure out if d2G should be cell 
or face variable? As I understood it, because I want to solve for X_var, this 
should be a cell variable and nothing else.

By higher dimension do you also mean 2D? I know 2D requires a lot of changes, 
but I thought 2D required only a change of grid. In any case, 2D is not my top 
priority at the moment.

Thank you again very much for your help!
Clara


> On 12 Jan 2018, at 12:04, Daniel Wheeler  wrote:
> 
> Hi Clara,
> 
> The following changes made it run for me. I turned d2G into a face
> variable and K into a variable that is updated in the sweep loop. I
> also removed redefining the equation. FiPy's designed so you only need
> to to create the equation once. I'm not sure what the values of d2G
> correspond to. Which spatial location do they correspond to? Faces or
> Cells. I'm assuming faces since there are 1001 of them. Also, this
> might need to be adjusted for higher dimensions.
> 
> Cheers,
> 
> Daniel
> 
> 
> $ git diff
> diff --git a/Spinodal_fipy_FeNi_varT_runs_DEBUG.py
> b/Spinodal_fipy_FeNi_varT_runs_DEBUG.py
> index 715e79b..a82387e 100644
> --- a/Spinodal_fipy_FeNi_varT_runs_DEBUG.py
> +++ b/Spinodal_fipy_FeNi_varT_runs_DEBUG.py
> @@ -289,20 +289,21 @@ def
> Get_spline_Chebnodes(X,N,T0,R,RefFe,RefNi,L0ex,L1ex,L2ex,Tmag,Bmag,p_fcc,A_f
> return d2G
> 
> 
> +
> ## THIS IS WHERE I WOULD LIKE TO HAVE A FIPY VARIABLE TYPE FOR d2G.
> d2G = 
> Get_spline_Chebnodes(Xf,100,T0,R,RefFe,RefNi,L0ex,L1ex,L2ex,Tmag,Bmag,p_fcc,A_fcc,Nv,Na)*1.0e-9
> -
> +d2G_var = FaceVariable(mesh=mesh, value=d2G)
> 
> ## Initial diffusion coefficient
> -Diff = Mob*d2G
> +Diff = Mob*d2G_var
> 
> 
> ## Initial gradient energy (depending on the interfacial energy we want)
> K = 2.0*kappa[id_mean_spino]
> -
> +K_var = Variable(K)
> 
> ## Initial equation
> -eq = TransientTerm(coeff=1.) == DiffusionTerm(Diff) -
> DiffusionTerm(coeff=(Mob, K))
> +eq = TransientTerm(coeff=1.) == DiffusionTerm(Diff) -
> DiffusionTerm(coeff=(Mob, K_var))
> 
> 
> plt.figure(figsize=(8,6))
> @@ -390,16 +391,16 @@ while time < duration:
> Mob = M[step_T][id_mean_T]
> 
> d2G = 
> Get_spline_Chebnodes(Xf,100,T_spino[step_T],R,RefFe,RefNi,L0ex,L1ex,L2ex,Tmag,Bmag,p_fcc,A_fcc,Nv,Na)*1.0e-9
> -
> -Diff = Mob*d2G
> +d2G_var[:] = d2G
> 
> K = 2.0*kappa[step_T]
> +K_var.setValue(K)
> print 'Mob, Diff, K'
> print Mob, Diff, K
> print 'MAX CONCENTRATION'
> print max(X_var)
> 
> -eq = TransientTerm(coeff=1.) == DiffusionTerm(Diff) -
> DiffusionTerm(coeff=(Mob, K))
> +
> 
> 
> print datetime.now() - startTime
> (END)
> 
> 
> On Thu, Jan 11, 2018 at 2:35 PM, Clara Maurel  wrote:
>> Hi Daniel,
>> 
>> Thank you for taking the time to look at this! My main code calls several
>> text files and subroutine. I attach everything that is needed to run the
>> code:
> 
> -- 
> Daniel Wheeler
> ___
> fipy mailing list
> fipy@nist.gov
> 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
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]

Re: Spline interpolation and fipy variable

[Daniel Wheeler]

Hi Clara,

The following changes made it run for me. I turned d2G into a face
variable and K into a variable that is updated in the sweep loop. I
also removed redefining the equation. FiPy's designed so you only need
to to create the equation once. I'm not sure what the values of d2G
correspond to. Which spatial location do they correspond to? Faces or
Cells. I'm assuming faces since there are 1001 of them. Also, this
might need to be adjusted for higher dimensions.

Cheers,

Daniel


$ git diff
diff --git a/Spinodal_fipy_FeNi_varT_runs_DEBUG.py
b/Spinodal_fipy_FeNi_varT_runs_DEBUG.py
index 715e79b..a82387e 100644
--- a/Spinodal_fipy_FeNi_varT_runs_DEBUG.py
+++ b/Spinodal_fipy_FeNi_varT_runs_DEBUG.py
@@ -289,20 +289,21 @@ def
Get_spline_Chebnodes(X,N,T0,R,RefFe,RefNi,L0ex,L1ex,L2ex,Tmag,Bmag,p_fcc,A_f
 return d2G


+
 ## THIS IS WHERE I WOULD LIKE TO HAVE A FIPY VARIABLE TYPE FOR d2G.
 d2G = 
Get_spline_Chebnodes(Xf,100,T0,R,RefFe,RefNi,L0ex,L1ex,L2ex,Tmag,Bmag,p_fcc,A_fcc,Nv,Na)*1.0e-9
-
+d2G_var = FaceVariable(mesh=mesh, value=d2G)

 ## Initial diffusion coefficient
-Diff = Mob*d2G
+Diff = Mob*d2G_var


 ## Initial gradient energy (depending on the interfacial energy we want)
 K = 2.0*kappa[id_mean_spino]
-
+K_var = Variable(K)

 ## Initial equation
-eq = TransientTerm(coeff=1.) == DiffusionTerm(Diff) -
DiffusionTerm(coeff=(Mob, K))
+eq = TransientTerm(coeff=1.) == DiffusionTerm(Diff) -
DiffusionTerm(coeff=(Mob, K_var))


 plt.figure(figsize=(8,6))
@@ -390,16 +391,16 @@ while time < duration:
 Mob = M[step_T][id_mean_T]

 d2G = 
Get_spline_Chebnodes(Xf,100,T_spino[step_T],R,RefFe,RefNi,L0ex,L1ex,L2ex,Tmag,Bmag,p_fcc,A_fcc,Nv,Na)*1.0e-9
-
-Diff = Mob*d2G
+d2G_var[:] = d2G

 K = 2.0*kappa[step_T]
+K_var.setValue(K)
 print 'Mob, Diff, K'
 print Mob, Diff, K
 print 'MAX CONCENTRATION'
 print max(X_var)

-eq = TransientTerm(coeff=1.) == DiffusionTerm(Diff) -
DiffusionTerm(coeff=(Mob, K))
+


 print datetime.now() - startTime
(END)


On Thu, Jan 11, 2018 at 2:35 PM, Clara Maurel  wrote:
> Hi Daniel,
>
> Thank you for taking the time to look at this! My main code calls several
> text files and subroutine. I attach everything that is needed to run the
> code:

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

Re: Spline interpolation and fipy variable

[Clara Maurel]

Hi Daniel,Thank you for taking the time to look at this! My main code calls several text files and subroutine. I attach everything that is needed to run the code:Spinodal_fipy_FeNi_varT_runs_DEBUG.pyfrom datetime import datetime
startTime = datetime.now()
import numpy as np
import scipy
from scipy import integrate
from operator import itemgetter
from collections import OrderedDict
from scipy.interpolate import InterpolatedUnivariateSpline
from fipy import *
from fipy.solvers.pysparse import LinearLUSolver as Solver
import Calculate_free_energy_theory_analytic_function as CalcFreeEnergyAn
import Calculate_diffusion_coefficient_2 as CalcDiff
from sympy import *
from scipy import interpolate
import matplotlib.pyplot as plt


#--#
# USEFUL FUNCTIONS #
#--#

def Get_closest_id(L,value):

return list(L).index(min(L, key=lambda x:abs(x-value)))

def calc_mol_frac(X_wt_EoI,AtomWeight_EoI,X_wt,AtomWeight):

# EoI = element of interest; X_wt and AtomWeight are lists of other compounds.

m_EoI = X_wt_EoI / AtomWeight_EoI
m_others = 0
for k in np.arange(len(X_wt)):
m_others += X_wt[k]/AtomWeight[k]

return m_EoI / float(m_EoI+m_others)

def Tot_energy(X,G,id,K):

return 1.0 * float((G(X)*1.0e-9 + (K/2.0)*(X.grad.mag)**2).cellVolumeAverage) * X.mesh.numberOfCells

def Poly_fit(x,coeffs):
n = len(coeffs)
f = 0
for j in np.arange(n):
f += coeffs[n-j-1]*x**j
return f

def d2Poly_fit(x,coeffs):
n = len(coeffs)
f = 0
for j in np.arange(2,n):
f += j*(j-1)*coeffs[n-j-1]*x**(j-2)
return f

def d2_spline(x,G):
tck = interpolate.splrep(x, G, s=0)
return interpolate.splev(x, tck, der=2)

#---#
# CONSTANTS #
#---#

# Composition
comp = 39.0

# Cooling rate: dT = -Cr dt
Myr_in_s = 3.1536e13 # 1 Myr in s
cool = 10.0

Cr = cool/Myr_in_s # K/s

# Boltzmann constant:
kb = 1.38064852e-5 # nm2 kg s-2 K-1

# Perfect gas constant
R = 8.31

# Avogadro number
Na = 6.02e23

# Atoms per m3 (FCC = 4 atoms per unit cell)
Nv = 4.0/(0.357e-9)**3

#---#
# VARIABLES #
#---#

X = np.arange(0.01,0.8,0.001)

#!!
dT = 0.01
#T = [350.0+273,500.0+273.0]
T = np.arange(210.0+273.0,405.0+dT+273.0,dT)
T = T[::-1]

#--#
# REFERENCE FOR THE FREE ENERGY DATA   #
# in J m-3 ==> kg m-1 s-2 ==> 1e-9 kg nm-1 s-2 #
#--#


ref_gibbs = 'DeKeyzer' # 'FactSage' or 'DeKeyzer'


# Pure elements FCC (Dinsdale 1991): a + b*T + c*T*ln(T) + d*T^2 + e*T^(-1) + f*T^3

aFe = -236.7
bFe = 132.416
cFe = -24.6643
dFe = -0.00375752
eFe = 77358.5
fFe = -5.89269e-8
RefFe = [aFe,bFe,cFe,dFe,eFe,fFe]

aNi = -5179.159
bNi = 117.854
cNi = -22.096
dNi = -4.8407e-3
eNi = 0.0
fNi = 0.0
RefNi = [aNi,bNi,cNi,dNi,eNi,fNi]

# Magnetic contribution (Factsage and De Keyzer et al.)

p_fcc = 0.28
A_fcc = (518.0/1125.0) + (11692.0/15975.0)*(1.0/p_fcc-1.0)

TmagFe = -201.0
TmagNi = 633.0
Tmag0 = 2133.0
Tmag1 = -682.0
Tmag = [TmagFe,TmagNi,Tmag0,Tmag1]

BmagFe = -2.1
BmagNi = 0.52
Bmag0 = 9.55
Bmag1 = 7.23
Bmag2 = 5.93
Bmag3 = 6.18
Bmag = [BmagFe,BmagNi,Bmag0,Bmag1,Bmag2,Bmag3]

# Excess Gibbs energy

if ref_gibbs == 'FactSage':

L00exFeNi = -12054.355
L01exFeNi = 3.27413
L0ex = [L00exFeNi,L01exFeNi]

L10exFeNi = 11082.131
L11exFeNi = -4.450770
L1ex = [L10exFeNi,L11exFeNi]

L20exFeNi = -725.8050
L21exFeNi = 0.0
L2ex = [L20exFeNi,L21exFeNi]

if ref_gibbs == 'DeKeyzer':

L00exFeNi = -14084.0
L01exFeNi = 3.27413
L0ex = [L00exFeNi,L01exFeNi]

L10exFeNi = 14000.0
L11exFeNi = -4.45077
L1ex = [L10exFeNi,L11exFeNi]

L20exFeNi = -3000.0
L21exFeNi = 0.0
L2ex = [L20exFeNi,L21exFeNi]

#-#
# SPINODAL BOUNDARIES #
#-#

T_spino = []
X_spino_low = []
X_spino_high = []

fp = open('Spinodal_FeNi_vol_'+ref_gibbs+'_0.5_fit.txt','r')

for j, line in enumerate(fp):
cols = line.split()
if j%50 == 0:
T_spino.append(float(cols[0])+273.0)
X_spino_low.append(float(cols[1]))
X_spino_high.append(float(cols[2]))
fp.close()

#-#
# DIFFUSION COEFFICIENT from Yang and Goldstein (2004) in m2 s-1 ==> 1e18 nm2 s-1 #
#-#

print "Calculating D..."

Tc, D_tmp, slope, intercept, lat_tmp = CalcDiff.Diff_coef_m(X,T_spino)

D = [[D_tmp[k][j]*1.0e18 for k in np.arange(len(X))] for j in np.arange(len(T_spino))]
Dfit = [[np.exp(slope[k]*(1.0/T_spino[j])+intercept[k])*1.0e18 for k in np.arange(len(X))] for j in np.arange(len(T_spino))]

#plt.semilogy([1.0/T_spino[k] for k in np.arange(len(T_spino))],[D[k][200] for k in np.arange(len(T))])
#plt.semilogy([1.0/T_spino[k] for k in np.arange(len(T_spino))],[Dfit[k][200] for k in np.arange

Re: Spline interpolation and fipy variable

[Daniel Wheeler]

Hi Clara,

Can you post the entire working code so I can try and debug?

Cheers,

Daniel

On Wed, Jan 10, 2018 at 4:43 PM, Clara Maurel  wrote:
> Hello,
>
> I asked a similar questions several months ago and someone nicely answered
> with some ideas, which however did not solve my problem. So I am just trying
> another time!
>
> I want to model the dynamics of a system using the  Cahn Hilliard equation,
> with a diffusion coefficient that depends on the cell variable
> (concentration X) and time.
>
> The diffusion coefficient (Diff) is proportional to the second derivative of
> the free energy (d2G). The latter function is not continuous and I would
> like to interpolate it with a B-spline (scipy.interpolate.UnivariateSpline).
> FYI, I first started with a polynomial interpolation: the fit was not
> satisfactory in terms of the physics I want to capture, but the fipy code
> worked well.
>
> Problems come when I want to use my spline function with the cell variable
> X: this always return an array where fipy would like a fipy variable, and I
> don’t know if there would be a way to do so. Would anyone have an idea on
> how I could do that?
> The suggestion I had before was to set Diff as a cell variable and update
> the value of Diff at each time step. However, doing so the behaviour of the
> system was not physical.
>
> Diff = CellVariable(mesh=mesh)
> Diff.setValue(d2G(Xf))
>
>
> Any other idea would be greatly appreciated! I would be happy to send my
> code if necessary.
>
> Thank you in advance,
> Clara
>
>
>
>
>
> Here is the part of my code in question:
>
> ## Mesh
> L = 200.
> nx = 1000
> dx = L/nx
>
> mesh = PeriodicGrid1D(nx=nx, dx=dx)
>
> ## Variable
> X_var = CellVariable(name=r"$X_{at}$", mesh=mesh, hasOld=True)
>
>
> ## Mean initial concentration
> X_mean = float(comp)/100.0
> X_mean_txt = str(comp)
> id_mean = Get_closest_id(X,X_mean)
> id_mean_spino = Get_closest_id(X_spino_low,X_mean)
>
>
> ## Initial temperature corresponding to X_mean on spinodal boundary
> T0 = T_spino[id_mean_spino]
>
> ## Initial thermal fluctuations
> noise = GaussianNoiseVariable(mesh=mesh, mean=X_mean,
> variance=0.1).value
> X_var[:] =  noise
> X_var.updateOld()
> Xf = X_var.arithmeticFaceValue
>
> ## Initial mobility
> Mob = M[id_mean_spino][id_mean]
>
> ## This interpolate the second derivative of the free energy with a B-spline
> and return Spline(Xf) WHICH IS AN ARRAY...
> d2G =
> Get_spline_Chebnodes(Xf,100,T0,R,RefFe,RefNi,L0ex,L1ex,L2ex,Tmag,Bmag,p_fcc,A_fcc,Nv,Na)
>
> ## Initial diffusion coefficient
> Diff = Mob*d2G
>
> ## Initial gradient energy (depending on the interfacial energy we want)
> K = 2.0*kappa[id_mean_spino]
>
> ## Initial equation
> eq = TransientTerm(coeff=1.) == DiffusionTerm(Diff) -
> DiffusionTerm(coeff=(Mob, K))
>
> ## Time and integration parameters
> time = (T0-T_spino[id_mean_spino])/Cr
> duration = (T0-T_spino[-2])/Cr # in s
> dt = 3.154e7 # 1 year in s
> dt_max = 3.154e13 # 1,000,000 years in s
> total_sweeps = 4
> tolerance = 1.0e-3
> step_T = id_mean_spino
>
> solver = Solver()
>
> ## Start looping
> while time < duration:
>
>
>
> res0 = eq.sweep(X_var, dt=dt, solver=solver)
>
>
>
> next_time_T = (T0-T_spino[step_T+1])/Cr
> print "Next temperature step is at: " + str(next_time_T/Myr_in_s) + "
> Myr"
> print "Temperature: "+str(T_spino[step_T]-273)
>
>
>
> while time < next_time_T:
>
>
>
> for sweeps in range(total_sweeps):
> res = eq.sweep(X_var, dt=dt, solver=solver)
>
> if res < res0 * tolerance:
> time += dt
> dt *= 1.1
> dt = min(dt_max, dt)
> X_var.updateOld()
> Xf = X_var.arithmeticFaceValue
>
>
> else:
> dt *= 0.8
> X_var[:] = X_var.old
> Xf = X_var.arithmeticFaceValue
>
> print "Arrived at a temperature step!"
> step_T += 1
> id_mean_T = Get_closest_id(X,np.mean(X_var))
>
> Mob = M[step_T][id_mean_T]
>
>
>
> d2G =
> Get_spline_Chebnodes(Xf,100,T_spino[step_T],R,RefFe,RefNi,L0ex,L1ex,L2ex,Tmag,Bmag,p_fcc,A_fcc,Nv,Na)
>
>
>
> Diff = Mob*d2G
>
>
>
> K = 2.0*kappa[step_T]
>
> eq = TransientTerm(coeff=1.) == DiffusionTerm(Diff) -
> DiffusionTerm(coeff=(Mob, K))
>
>
>
>
>
>
>
>
>
>
> ___
> fipy mailing list
> fipy@nist.gov
> http://www.ctcms.nist.gov/fipy
>   [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
>



-- 
Daniel Wheeler

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

Spline interpolation and fipy variable

[Clara Maurel]

Hello,

I asked a similar questions several months ago and someone nicely answered with 
some ideas, which however did not solve my problem. So I am just trying another 
time!

I want to model the dynamics of a system using the  Cahn Hilliard equation, 
with a diffusion coefficient that depends on the cell variable (concentration 
X) and time.

The diffusion coefficient (Diff) is proportional to the second derivative of 
the free energy (d2G). The latter function is not continuous and I would like 
to interpolate it with a B-spline (scipy.interpolate.UnivariateSpline). FYI, I 
first started with a polynomial interpolation: the fit was not satisfactory in 
terms of the physics I want to capture, but the fipy code worked well.

Problems come when I want to use my spline function with the cell variable X: 
this always return an array where fipy would like a fipy variable, and I don’t 
know if there would be a way to do so. Would anyone have an idea on how I could 
do that?
The suggestion I had before was to set Diff as a cell variable and update the 
value of Diff at each time step. However, doing so the behaviour of the system 
was not physical.

Diff = CellVariable(mesh=mesh)
Diff.setValue(d2G(Xf))


Any other idea would be greatly appreciated! I would be happy to send my code 
if necessary.

Thank you in advance,
Clara





Here is the part of my code in question:

## Mesh
L = 200.
nx = 1000
dx = L/nx

mesh = PeriodicGrid1D(nx=nx, dx=dx)

## Variable
X_var = CellVariable(name=r"$X_{at}$", mesh=mesh, hasOld=True)


## Mean initial concentration
X_mean = float(comp)/100.0
X_mean_txt = str(comp)
id_mean = Get_closest_id(X,X_mean)
id_mean_spino = Get_closest_id(X_spino_low,X_mean)


## Initial temperature corresponding to X_mean on spinodal boundary
T0 = T_spino[id_mean_spino]

## Initial thermal fluctuations
noise = GaussianNoiseVariable(mesh=mesh, mean=X_mean, variance=0.1).value
X_var[:] =  noise
X_var.updateOld()
Xf = X_var.arithmeticFaceValue

## Initial mobility
Mob = M[id_mean_spino][id_mean]

## This interpolate the second derivative of the free energy with a B-spline 
and return Spline(Xf) WHICH IS AN ARRAY...
d2G = 
Get_spline_Chebnodes(Xf,100,T0,R,RefFe,RefNi,L0ex,L1ex,L2ex,Tmag,Bmag,p_fcc,A_fcc,Nv,Na)

## Initial diffusion coefficient
Diff = Mob*d2G

## Initial gradient energy (depending on the interfacial energy we want)
K = 2.0*kappa[id_mean_spino]

## Initial equation
eq = TransientTerm(coeff=1.) == DiffusionTerm(Diff) - DiffusionTerm(coeff=(Mob, 
K))

## Time and integration parameters
time = (T0-T_spino[id_mean_spino])/Cr
duration = (T0-T_spino[-2])/Cr # in s
dt = 3.154e7 # 1 year in s
dt_max = 3.154e13 # 1,000,000 years in s
total_sweeps = 4
tolerance = 1.0e-3
step_T = id_mean_spino

solver = Solver()

## Start looping
while time < duration:

res0 = eq.sweep(X_var, dt=dt, solver=solver)

next_time_T = (T0-T_spino[step_T+1])/Cr
print "Next temperature step is at: " + str(next_time_T/Myr_in_s) + " Myr"
print "Temperature: "+str(T_spino[step_T]-273)

while time < next_time_T:

for sweeps in range(total_sweeps):
res = eq.sweep(X_var, dt=dt, solver=solver)

if res < res0 * tolerance:
time += dt
dt *= 1.1
dt = min(dt_max, dt)
X_var.updateOld()
Xf = X_var.arithmeticFaceValue

else:
dt *= 0.8
X_var[:] = X_var.old
Xf = X_var.arithmeticFaceValue

print "Arrived at a temperature step!"
step_T += 1
id_mean_T = Get_closest_id(X,np.mean(X_var))

Mob = M[step_T][id_mean_T]

d2G = 
Get_spline_Chebnodes(Xf,100,T_spino[step_T],R,RefFe,RefNi,L0ex,L1ex,L2ex,Tmag,Bmag,p_fcc,A_fcc,Nv,Na)

Diff = Mob*d2G
   
K = 2.0*kappa[step_T]

eq = TransientTerm(coeff=1.) == DiffusionTerm(Diff) - 
DiffusionTerm(coeff=(Mob, K))











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 ]