Thank you, that’s helpful.
I introduced another auxiliary variable so that every equation only has 2nd
order terms (eq1 was 4th order in both psi and phi).
I used the scipy solvers, which at least solve. I don’t know why the PETSc
solvers are complaining about the diagonals; I’ll investigate that, separately.
While the scipy solvers get solutions to the coupled equations, they still
don’t make a lot of sense. To try to get a handle on things, I slowed the
problem down and got rid of the exponentially increasing time steps. When I do
that, there seems to be some checkerboard instability, particularly on the
auxiliary variables.
I don’t know what the problem is, but have some suspicions on things to examine:
- I suspect the mesh is under-resolved. If you don’t analytically know what
sort of interface thickness to expect, you can at least experiment with finer
meshes to see if things stabilize. When I do that, I can get solutions that
don’t immediately go to zero, and the auxiliary variables retain some sort of
reasonable structure without checker-boarding.
- The checker-boarding makes me wonder if there’s a sign error for some of the
higher order terms. I’m not familiar with this phase field curvature model, but
I’ve put together a notebook that implements a basic phase field crystal model
from [Provatas &
Elder](https://www.physics.mcgill.ca/~provatas/papers/Phase_Field_Methods_text.pdf)
that may give some insights into how to do the separation of variables of a
high-order PDE in a way that FiPy likes:
https://github.com/guyer/phase_field_crystal/blob/main/Figure_8_6.ipynb
I’ll see about doing one of the higher-order conserved models from Provatas &
Elder.
I’ll continue to explore why PETSc has a problem building the matrix.
On Feb 9, 2022, at 1:48 AM, Matan Mussel
<[email protected]<mailto:[email protected]>> wrote:
Hi Jonathan,
Thank you for getting back to me. You are correct, these are the same equations
I previously asked about. This version is almost the original form (which is
attached to this message). Indeed, when using
eq1.solve(dt=dt)
eq2.solve(dt=dt)
I do not get an error message. Unfortunately, I don't get anything reasonable
either. For simplicity, I set sigma to zero and ignore the third equation just
to see if something reasonable occurs for phi (which should maintain its
boundary). Unfortunately, already after the first iteration, the phi field
becomes almost zero everywhere (see attached two images at t=0 and t=0.0025.
May I ask what is the difference in calculation between using:
eq1.solve(dt=dt)
eq2.solve(dt=dt)
Solving these equations separately means that each equation builds a solution
matrix that is implicit in a single variable and any dependence of eq2 on var1
and of eq1 on var2 can only be established by sweeping. All terms involving
other variables are handled explicitly.
Here, you should specify which variable each equation applies to (you might get
lucky, but I wouldn’t count on it):
eq1.solve(var=var1, dt=dt)
eq2.solve(var=var2, dt=dt)
and
eq = eq1 & eq2
eq.solve(dt=dt)
Coupling the equations means that FiPy builds a block matrix, where the block
columns correspond to the individual variables and the block rows correspond
the equations. The entire block matrix is inverted to simultaneously obtain an
implicit solution to all variables. Sweeping may still be necessary to handle
non-linearity.
In this case, do not specify the variable to solve; FiPy handles it.
When encountering difficulties, solving separately is advisable. Some sets of
equations won’t converge at all, but usually, convergence is just harder than
coupled. On the other hand, coupled equations can be fussier until you get them
“right”.
- Jon
?
Thank you and it is really great that you guys at NIST provide this lively
support.
Matan
On Tue, Feb 8, 2022 at 9:25 PM 'Guyer, Jonathan E. Dr. (Fed)' via fipy
<[email protected]<mailto:[email protected]>> wrote:
Are these the same equations you were asking about in
https://github.com/usnistgov/fipy/issues/835? They seem closely related, but
not identical, AFAICT.
Regardless, as explained there, you cannot mix higher-order diffusion terms
with coupled equations.
Can you show the equations you started with, before you started trying to put
them in a form for FiPy?
On Feb 8, 2022, at 12:44 PM, Matan Mussel
<[email protected]<mailto:[email protected]>> wrote:
Hello everyone,
I am new to FiPy and interested in solving a certain curvature model. I thought
using FiPy will make my life easier, but unfortunately I have been struggling
with this for quite some time trying different versions of introducing the
equations into FiPy.
I decided to try this community, hoping to get some insight on what I am doing
wrong. The model includes three variables (two dynamic, one auxiliary) and
three coupled equations running on a 2D grid. Please see a summary of the
model in the attached PDF file. I am also attaching a minimalistic version of
the code. For this version the error message is: "[0] Matrix is missing
diagonal entry 2500" (with 2500 being nx*ny).
Many thanks in advance,
Matan
--
To unsubscribe from this group, send email to
[email protected]<mailto:[email protected]>
View this message at https://list.nist.gov/fipy
---
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected]<mailto:[email protected]>.
<CurvatureModelEquations.pdf><MinimalCurvatureModel.py>
--
To unsubscribe from this group, send email to
[email protected]<mailto:fipy%[email protected]>
View this message at https://list.nist.gov/fipy
---
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected]<mailto:[email protected]>.
<phaseFieldCurvatureModel.pdf><CH_0.0025.png><CH_0.0000.png>
--
To unsubscribe from this group, send email to [email protected]
View this message at https://list.nist.gov/fipy
---
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected].
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
from fipy import CellVariable, Grid2D, GaussianNoiseVariable, DiffusionTerm, TransientTerm, ImplicitSourceTerm, Viewer, MultiViewer
from fipy.tools import numerix
# create a 2D solution mesh
nx = ny = 50
dx = dy = 0.001
mesh = Grid2D(nx=nx, ny=ny, dx=dx, dy=dy)
phi0 = -1.
psi0 = phi0**3 - phi0
sigma0 = 0.001
xi0 = 0.
# create variables
phi = CellVariable(name=r"$\phi$", mesh=mesh, value=phi0, hasOld=True)
psi = CellVariable(name=r"$\psi$", mesh=mesh, value=psi0, hasOld=True)
sigma = CellVariable(name=r"$\sigma$", mesh=mesh, value=sigma0) #hasOld=True
xi = CellVariable(name=r"$\xi$", mesh=mesh, value=xi0)
# Initial condition: initiate phi as a circle of 1. in a field of -1.
Lx = nx * dx
Ly = ny * dy
r0 = np.sqrt(Lx**2 + Ly**2)/4 #radius
x0 = Lx/2 #position of center of circle
y0 = Ly/2
phi.setValue(-1)
phi.setValue(1., where=(mesh.x - x0)**2 + (mesh.y - y0)**2 < r0**2)
# Parameters
kappa = 0.01
epsilon = 0.2
a0 = 10
# Equations
eq1 = (TransientTerm(var=phi) == DiffusionTerm(coeff=kappa, var=xi))
eq2 = (ImplicitSourceTerm(coeff=1., var=psi)
== ImplicitSourceTerm(coeff=(phi**2 - 1), var=phi)
- DiffusionTerm(coeff=epsilon**2, var=phi))
eq3 = (TransientTerm(var=sigma)
== phi.grad.dot(phi.grad)/2 - a0)
eq4 = (ImplicitSourceTerm(coeff=1., var=xi)
== ImplicitSourceTerm(3 * phi**2 - 1, var=psi)
- DiffusionTerm(coeff=epsilon**2, var=psi)
+ DiffusionTerm(coeff=epsilon**2 * sigma, var=phi))
eq = eq1 & eq2 & eq3 & eq4
# Viewers
viewers = [Viewer(vars=var, axes=plt.subplot(2,2,i+1))
for i, var in enumerate((phi, psi, sigma, xi))]
viewer = MultiViewer(viewers=viewers)
viewer.plot()
# Use exponentially increasing time to accelarate evolution as system relaxes towards the equilibrium state
dexp = -9
elapsed = 0.
duration = 500
# solving the equation
while elapsed < duration:
dt = min(100, numerix.exp(dexp))
elapsed += dt
# dexp += 0.01
for sweep in range(3):
eq.sweep(dt=dt)
viewer.plot()
