Solver problems with a different KWC Phase field implementation

[vitesh shah]

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

[vitesh shah]

I added the .py file as attaching the notebook is somehow not supported

[Guyer, Jonathan E. Dr. (Fed)]

Vitesh -

I get this warning from SciPy and PySparse, but not from PETSc or Trilinos. Evidently, there’s something LU doesn’t like about it. I don’t think it actually matters, though. The solution proceeds anyway. Often, the initial matrix can raise objections that don’t cause the overall evolution of the problem to fail.

A bigger issue is that you have several cases where you are casting away the “Variable” nature of your variables. Everywhere where you write `something.value`, you get the instantaneous value of `something`; it won’t update as the problem involves. Just use `something`.

-    return 1. + beta_r*(1.-np.exp(-beta_s*phi.value*one_over(phi.grad.mag.value)))

+    return 1. + beta_r*(1.-np.exp(-beta_s*phi*one_over(phi.grad.mag)))

-    diffusionCoeff = (g(phiFace,a).value * IGamma.value + 2.0*h(phiFace,d,e))

+    diffusionCoeff = (g(phiFace,a) * IGamma + 2.0*h(phiFace,d,e))

More subtly, appending `()` to a FiPy Variable also returns its instantaneous value (`something()` is equivalent to `something.value`).

 def build_E_st_equation(phi,theta,C_d,E_st):

-    prefactor = -E_st*((phi - phi.old()) > 0.0)

+    prefactor = -E_st*((phi - phi.old) > 0.0)

     prefactor = prefactor*localize_E(theta)

-    return TransientTerm() == prefactor*(phi - phi.old())/dt

+    return TransientTerm() == prefactor*(phi - phi.old)/dt

Finally, you need to be careful with calling functions. You want to be sure that your functions return a Variable expression, not just the instantaneous value. Most of your helper functions are OK, but `one_over` won’t worked as desired because `np.where` and `np.isclose` discard the Variable character and just return numpy arrays. I recommend:

-  return np.where(np.isclose(f,0.0,atol=gamma),1.0/gamma,np.tanh(f/gamma)/f)

+  return 1.0 / (f + gamma)

I suppose you could do 

+  return 1.0 / (f + (abs(f) < gamma) * gamma)

but I doubt it makes much difference.

After that, the solution evolves, but it’s not stable. 

Dropping the time step to 2e-4 or 2e-5 gets a smooth evolution toward something like Fig. 2(b) in Abrivard et al. (I don’t know if it’s supposed to, as I haven’t actually read that paper, but it does).

- Jon

--
To unsubscribe from this group, send email to fipy+uns...@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+uns...@list.nist.gov.
<KWC_with_stored_energy_example.py>