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 <cmau...@mit.edu> 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.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))
>
>
>
>
>
>
>
>
>
>
> _______________________________________________
> 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 ]
