Clara Maurel
2018-01-10 21:43:29 UTC
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.00001).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))
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.00001).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))