Re: Spline interpolation and fipy variable
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
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 ]
