I’ve been able to get your full system of equations to solve:
https://gist.github.com/guyer/0d3dcb36303bfe0a6a9cdfb327ec3ff2
Notable changes I needed to make (in addition to transcribing it to a Jupyter
notebook):
- \epsilon is on the order of the interface thickness in Cahn-Hilliard-like
models. Campelo & Hernández-Machado say in Sec. 4 they set \epsilon "to be
equal to the mesh size”, so I reduced `epsilon` from 0.2 to `dx`. Your mesh was
actually over-resolved; the interface thickness of phi was vastly larger than
your solution domain.
- The solutions then had some structure, but were unstable. I reduced the time
step to 1e-5 (and got rid of the exponentially increasing time steps; you don’t
need any of that).
- At this point, without sigma, the solution was stable with a nice diffuse
interface on phi.
- I tried different forms of both the sigma equation and the sigma-containing
term in the xi equation (the mess inside the Laplacian of the phi equation). I
had to add an arbitrary factor to reduce the strength of the coupling, but I
was finally able to get stable solutions to all four equations. Instead of the
arbitrary factor I introduced, you’ll want to check the derivation of the sigma
equation, which I don’t see in Campelo & Hernández-Machado; there should be
some scaling between sigma and phi.
- Your notes on the sigma equation said
$\frac{\partial\sigma}{\partial t} &=
\frac{1}{2}\lvert\nabla\sigma\rvert^2 - a_0$
but your Python code says
eq3 = (TransientTerm(var=sigma) == phi.grad.dot(phi.grad)/2 - a0)
which is
$\frac{\partial\sigma}{\partial t} &= \frac{1}{2}\lvert\nabla\phi\rvert^2
- a_0$
which makes more sense (otherwise there’s no coupling between sigma and the
rest of the problem).
On Feb 16, 2022, at 6:32 PM, Guyer, Jonathan E. Dr. (Fed)
<[email protected]<mailto:[email protected]>> wrote:
I’ve now had a chance to post a 6th-order, conserved PFC model:
https://github.com/guyer/phase_field_crystal/blob/main/Figure_8_14.ipynb
It doesn’t work with PETSc, but it works with SciPy, Trilinos, and PySparse
(the latter two only with Python 2.7, unfortunately).
On Feb 14, 2022, at 10:39 AM, Guyer, Jonathan E. Dr. (Fed)
<[email protected]<mailto:[email protected]>> wrote:
Matan -
On Feb 14, 2022, at 3:57 AM, Matan Mussel
<[email protected]<mailto:[email protected]>> wrote:
Can you please explain in more detail how to change the solver? In the terminal
window I typed "FIPY_SOLVERS=scipy". However, if I then type
"python3 -c "from fipy import *; print(DefaultSolver)", the output I receive is
<class 'fipy.solvers.petsc.linearGMRESSolver.LinearGMRESSolver'>. This, to the
best of my understanding, indicates that the solver I am using is still PETSc.
Indeed, I still get the error about the diagonals and cannot see the solutions
you describe.
When you put `FIPY_SOLVERS=scipy` on a line by itself, the shell sets the
variable and then immediately forgets it.
You either need to set the variable as part of the same shell command:
FIPY_SOLVERS=scipy python3 -c "from fipy import *; print(DefaultSolver)"
or you need to make the variable part of the environment in the shell session.
For bash:
export FIPY_SOLVERS=scipy
python3 -c "from fipy import *; print(DefaultSolver)"
Alternatively, you can pass a flag:
python3 -c "from fipy import *; print(DefaultSolver)" --scipy
Depending on your platform and shell, you might need to do:
python3 -c "import os; os.environ['FIPY_SOLVERS'] = ‘scipy’; from fipy import
*; print(DefaultSolver)"
I was also thinking that there may be a sign-error but I double checked with
past published works and all use the same equation with similar signs (e.g.,
equation 29 here<https://link.springer.com/article/10.1140/epje/i2005-10079-5>).
Thank you for the reference; I’ll take a look.
I’ve implemented the 6th-order, conserved PFC model in Provatas & Elder (their
Eq. (8.89)) and get the same general behavior with PETSc. Solutions look
reasonable with SciPy, though. Given that, I’d guess that you’re right that you
don’t have sign errors (I couldn’t identify any, either), but that your mesh or
time steps may not be adequately resolved for your choice of parameters. I need
to do a few more diagnostics, but, hopefully today, I’ll post the 6th-order PFC
notebook so you can use it as a guide.
I think I understand why PETSc is failing. When FiPy builds the matrix for
coupled equations, it has an empty diagonal block. This is mathematically OK,
because we are free to swap rows or columns. SciPy doesn’t seem to mind, but my
guess is that PETSc is expecting a block-diagonal matrix. I’ll need to research
whether that constraint can be relaxed (PETSc takes a *lot* of options) or if
FiPy can be forced to build a block-diagonal matrix (my recollection of our
algorithm is that it should be already, so that’s something to figure out,
anyway).
- Jon
Thanks again and best regards,
Matan
On Thu, Feb 10, 2022 at 10:43 PM 'Guyer, Jonathan E. Dr. (Fed)' via fipy
<[email protected]<mailto:[email protected]>> wrote:
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<https://gcc02.safelinks.protection.outlook.com/?url=https%3A%2F%2Fwww.physics.mcgill.ca%2F~provatas%2Fpapers%2FPhase_Field_Methods_text.pdf&data=04%7C01%7Cjonathan.guyer%40nist.gov%7C3fb7ce1509494fa08f9308d9ef985070%7C2ab5d82fd8fa4797a93e054655c61dec%7C1%7C0%7C637804259729967807%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000&sdata=A4CxNkNL8xnDX3JLyAmu0g0s%2FkHFYu7XXdSFZZ9KIgQ%3D&reserved=0>)
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]<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]>.
--
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].
