Hello,
I am currently using Fipy, to solve a set of phase field equations. These equations are based on the KWC model but with an additional extension made to it, which enables it to account for stored energy of a material to be a driving force. Currently, I have three equations which I attempt to solve using Fipy: 1. Phase equation 2. Theta (misorientation) equation 3. Equation for the stored energy. If I solve without the stored energy, the Fipy solver works fine and I get expected behavior (in 1D at least). But when I also attempt to solve equation for stored energy, I get a runtime warning: “/nethome/v.shah/.local/lib/python3.8/site-packages/fipy/solvers/scipy/linearLUSolver.py:41: RuntimeWarning: invalid value encountered in double_scalars if (numerix.sqrt(numerix.sum(errorVector**2)) / error0) <= self.tolerance:” So I suspect that solver is not handling these equations properly. And this affects the solution that we get. We expect that for 1D case, the phase boundary should move when stored energy is introduced (in normal phase field it should remain stationary). But this behavior is not observed. Also when I compare the analytical rates and the rates that I observe from Fipy and they have quite a different behavior. I am attaching a jupyter notebook here, so that the whole code (with references to the papers) can be looked at. (link: if the attachment doesnt work) It would be great, if you can give me some suggestions or point out some mistakes that I might have made. With regards, Vitesh Shah -- To unsubscribe from this group, send email to fipy+unsubscr...@list.nist.gov View this message at https://list.nist.gov/fipy --- To unsubscribe from this group and stop receiving emails from it, send an email to fipy+unsubscr...@list.nist.gov.
#!/usr/bin/env python # coding: utf-8 # # KWC model with stored energy as extra driving force # We try to implement a phase field model described by Abrivard et al. (2012). # The link to the paper is https://doi.org/10.1080/14786435.2012.713135 # # This notebook uses the example code of for 1D as an example: # https://www.ctcms.nist.gov/fipy/examples/phase/generated/examples.phase.impingement.mesh40x1.html # ### Import necessary modules # In[1]: import numpy as np import os from scipy import constants # In[2]: import matplotlib as mpl mpl.use('Agg') import matplotlib.pyplot as PyPlot import matplotlib.ticker as ticker from matplotlib.font_manager import FontProperties get_ipython().run_line_magic('matplotlib', 'inline') # In[3]: from fipy import CellVariable, ModularVariable, Grid1D, TransientTerm, DiffusionTerm, ExplicitDiffusionTerm, ImplicitSourceTerm, GeneralSolver, Viewer from fipy.tools import numerix from fipy.tools.parser import parse numberOfElements = parse('--numberOfElements', action = 'store', type = 'int', default = 40000) numberOfSteps = parse('--numberOfSteps', action = 'store', type = 'int', default = 100) # ### Build geometry # In[4]: steps = numberOfSteps nx = 200 Lx = 1.0 * nx / 100. dx = Lx / nx mesh = Grid1D(dx=dx, nx=nx) # In[5]: dx,Lx # ### Parameters # In[6]: timeStepDuration = 2.0E-03 #0.02 beta_Theta_0 = 1.126e-15 beta_phi_0 = 1.126e-17 #beta = 1e5 mu = 1e3 gamma= 1e3 epsilon = 0.8 s = 0.01 alpha = 0.04 #0.04 #0.04 omega = 7.5 a = np.array([1.0,0.95,1.35]) # static parameters - 2.1,0.95,1.35 d = 0.45 e = 1.0E-03 # the second paper of Abrivard uses 0.015 originally 1e-03 beta_s = 3.0e-06 beta_r = 1.e6 Q = 2.72e-19 C_d = 20.0 # In[7]: temperature = 573. # In[8]: dt=timeStepDuration # ## Defining functions for calculations # In[9]: def beta(beta_0,Q,T): return beta_0 * np.exp(Q/(constants.k*T)) # In[10]: def tau_phi(beta_phi_0,T,Q,omega): return beta(beta_phi_0,Q,T)*omega**2. # In[11]: def tau_Theta(beta_Theta_0,T,Q,omega,phi,beta_r,beta_s): return beta(beta_Theta_0,Q,T)*omega**2. # In[12]: def P_m(phi): return phi**2. # In[13]: def P_r(phi,beta_r,beta_s): return 1. + beta_r*(1.-np.exp(-beta_s*phi.value*one_over(phi.grad.mag.value))) # In[14]: def g(phi,a): return a[0]*phi + a[1]*phi**2. + a[2]*phi**3. # In[15]: def g_dash(phi,a): return a[0] + 2.*a[1]*phi + 3.*a[2]*phi**2. # In[16]: def h(phi,d,e): return d*phi**2. + e # In[17]: def h_dash(phi,d): return 2.*d*phi # In[18]: def one_over(f,gamma=1.e-12): """ Gamma in the function of Kobayashi is upper limit of 1/0.0. Therefore, the Gamma is defined as 1/gamma here. Parameters: ----------- f : np.ndarray Field gamma : float tolerance to check for zero """ return np.where(np.isclose(f,0.0,atol=gamma),1.0/gamma,np.tanh(f/gamma)/f) # In[19]: def localize_E(theta): return numerix.tanh(theta.old.grad.mag**2.0) # ## Defining phase field variables # ### Phase variable $ \phi $ # In[20]: phi = CellVariable( name='phase field', mesh=mesh, value=1., hasOld=1 ) # the original example doesn#t have the attribute hasOld, but we need it here for the evolution equation of stored energy # ## Orientation variable $ \theta $ # In[21]: theta = ModularVariable( name='theta', mesh=mesh, value=1., hasOld=1 ) # In[22]: theta.setValue(0., where=mesh.cellCenters[0] > Lx / 2.) # ## Stored energy field $ E_{st} $ # In[23]: E_st = CellVariable(name='Stored_energy', mesh=mesh,value = 0.0, hasOld=True) # If the value of $E_{st}$ is yero everywhere, then it behaves like regular KWC implementation. # Is setValue if made 0.0, we get regular KWC behaviour # In[24]: E_st.setValue(0.9, where=mesh.cellCenters[0] > Lx / 2.) # ## Bulding equations # ### Evolution equation for phase # # $ \tau_\phi \frac{\partial \phi}{\partial t}$ = $ \alpha^2 \nabla^2 \phi - f^{'}(\phi) - C E_{st} - g^{'}(\phi) | \nabla \theta | - h^{'}(\phi) | \nabla \theta |^2 $ # # For this case we have assumed, $f(\phi) = \frac{\omega^2}{2} (1 - \phi)^2$ # # The expanded form of the equation split in implicit and source terms is: # # $\tau_\phi \frac{\partial \phi}{\partial t}$ = $ \alpha^2 \nabla^2 \phi + \omega^{2} - \omega^{2} \phi - C E_{st} - a_{0}| \nabla \theta | - 2 a_{1} \phi | \nabla \theta | - 3a_{2} \phi^{2} | \nabla \theta | - 2 d \phi | \nabla \theta |^2 $ # # $\tau_\phi \frac{\partial \phi}{\partial t}$ = $ \alpha^2 \nabla^2 \phi - \bigl( \omega^{2} + \bigl[ 2 a_{1} + 3a_{2} \phi + 2 d | \nabla \theta | \bigr] | \nabla \theta | \bigr) \phi + \omega^{2} - C E_{st} - a_{0}| \nabla \theta |$ # In[25]: def buildPhiEquation(phi, theta, E_st): old_theta_grad = theta.old.grad.mag implicitSource = omega**2.0 implicitSource += (2.0*a[1] + 3.0*a[2]*phi + 2*d*old_theta_grad)*old_theta_grad sourceCoeff = omega**2.0 - a[0]*old_theta_grad - E_st*(omega**2.0) #assuming c/w^2 = 1, therefore c = w^2 return TransientTerm(tau_phi(beta_phi_0,temperature,Q,omega)) == ExplicitDiffusionTerm(alpha**2) - ImplicitSourceTerm(implicitSource) + sourceCoeff # In[26]: phiEq = buildPhiEquation(phi, theta, E_st) # ### Evolution equation for theta # # $ \tau_\theta \frac{\partial \theta}{\partial t}$ = $\nabla \biggl ( 2 h (\phi) + \frac{g(\phi)}{| \nabla \theta | } \biggr) \nabla \theta$ # # Expanding it leads to: # $ \tau_\theta \frac{\partial \theta}{\partial t}$ = $\nabla \biggl ( 2 (d \phi^2 + e) + \frac{a_0 \phi + a_{1}\phi^2 + a_{3}\phi^3}{| \nabla \theta| } \biggr) \nabla \theta $ # # Moreover, $\tau_\theta$ is defined as, # # $\tau_\theta $ = $\beta_{\theta} \mathrm{exp} (\frac{Q}{k_{B} T}) P_{m}(\phi) P_{r}(\phi,\nabla \phi)$ # # $\tau_\theta $ = $\beta_{\theta} \mathrm{exp} (\frac{Q}{k_{B} T}) \phi^2 (1 + \beta_{r} [ 1 - \mathrm{exp}(-\beta_{s} \frac{\phi}{| \nabla \phi |})])$ # In[27]: def buildThetaEquation(phi, theta): phiFace = phi.arithmeticFaceValue phiSq = phiFace * phiFace gradMag = theta.faceGrad.mag eps = 1. / gamma / 10. gradMag += (gradMag < eps) * eps IGamma = (gradMag > 1. / gamma) * (1 / gradMag - gamma) + gamma diffusionCoeff = (g(phiFace,a).value * IGamma.value + 2.0*h(phiFace,d,e)) thetaGradDiff = theta.faceGrad - theta.faceGradNoMod sourceCoeff = (diffusionCoeff * thetaGradDiff).divergence return TransientTerm(tau_Theta(beta_Theta_0,temperature,Q,omega,phi,beta_r,beta_s) * P_m(phi) * P_r(phi,beta_r,beta_s)) == DiffusionTerm(diffusionCoeff) + sourceCoeff # In[28]: thetaEq = buildThetaEquation(phi, theta) # ### Evolution of stored energy # The evolution rate of the stored energy $E_{st}$ is given as, # # $\frac{\partial{E_{st}}}{\partial t}$ = $- E_{st} C_{d} \mathrm{tanh}(|\nabla \theta|^2) \frac{\partial \phi}{\partial t}$ if $\frac{\partial \phi}{\partial t} > 0$ # # $\frac{\partial{E_{st}}}{\partial t}$ = 0 if $\frac{\partial \phi}{\partial t} \leq 0$ # # This equation is solved as described here: # # https://tsapps.nist.gov/publication/get_pdf.cfm?pub_id=854461 # In[29]: def build_E_st_equation(phi,theta,C_d,E_st): prefactor = -E_st*((phi - phi.old()) > 0.0) prefactor = prefactor*localize_E(theta) return TransientTerm() == prefactor*(phi - phi.old())/dt # In[30]: E_st_eq = build_E_st_equation(phi,theta,C_d,E_st) # In[31]: if __name__ == '__main__': phiViewer = Viewer(vars=phi, datamin=0., datamax=1.) thetaViewer = Viewer(vars=theta, datamin=-numerix.pi/3.0, datamax=numerix.pi/3.0) EstViewer = Viewer(vars=E_st, datamin=0., datamax=1.) phiViewer.plot() thetaViewer.plot() EstViewer.plot() # #### Solve per time step # In[32]: from builtins import range for i in range(steps): theta.updateOld() phi.updateOld() E_st.updateOld() thetaEq.solve(theta, dt=timeStepDuration, solver=GeneralSolver(iterations=3000, tolerance=1e-15)) phiEq.solve(phi, dt=timeStepDuration, solver=GeneralSolver(iterations=3000, tolerance=1e-15)) E_st_eq.solve(E_st,dt=timeStepDuration, solver=GeneralSolver(iterations=3000, tolerance=1e-15)) if __name__ == '__main__': if (i+1)%100 == 0: phiViewer.plot() thetaViewer.plot() # ## Visualization # In[33]: phiViewer.plot() # In[34]: thetaViewer.plot() # In[35]: PyPlot.plot((phi-phi.old).value) PyPlot.title('phi_dot') # #### Different rates in the $\phi$ equation plotted # If we sum these indivual contributions, we get the 'expected' evolution rate of $\phi$ # In[36]: PyPlot.figure(figsize=(12,5)) PyPlot.plot((alpha**2.0)*np.gradient(np.gradient(phi.value)),label='divergence_phi') PyPlot.plot((omega**2.0)*(1.0 - phi.value),label='$ \omega^2 (1 -\phi) $') PyPlot.plot(-(g_dash(phi,a)*theta.old.grad.mag).value,label='$ -g(\phi) grad(theta) $') PyPlot.plot(-(h_dash(phi,d)*(theta.old.grad.mag)**2.0).value,label='$-h(\phi) sq(grad(Theta)) $') PyPlot.plot(-E_st.value*(omega**2.0),label='$- E \omega^2$') PyPlot.legend() # In[37]: example_phi_dot = (alpha**2.0)*np.gradient(np.gradient(phi.value)) + (omega**2.0)*(1.0 - phi.value) - (g_dash(phi,a)*theta.old.grad.mag).value - (h_dash(phi,d)*(theta.old.grad.mag)**2.0).value - E_st.value*(omega**2.0) # In[38]: PyPlot.plot(example_phi_dot) PyPlot.title('phi_dot_expected') # In[ ]:
