On Feb 20, 2022, at 10:38 AM, Matan Mussel
<[email protected]<mailto:[email protected]>> wrote:
Jonathan,
A few comments regarding your version of the code:
1. Equation 4 is missing one term (see equation29 in Campelo 2006) which is: "+
fp.ImplicitSourceTerm(coeff=epsilon**2 * phi.faceGrad.divergence, var=sigma)"
Compelo 2006 Eq. (29) has:
\ldots + \epsilon^2 \bar{sigma}(\vec{x})\nabla^2\phi +
\epsilon^2\nabla\bar{\sigma}(\vec{x})\cdot\nabla\phi
but these two terms are equivalent to
\ldots + \nabla\cdot\left(\epsilon^2 \bar{sigma}(\vec{x})\nabla\phi \right)
In this divergence form, it is precisely represented by
fp.DiffusionTerm(coeff=epsilon**2 * sigma, var=phi)
or by
fp.ConvectionTerm(coeff=epsilon**2 * phi.faceGrad, var=sigma)
The second form (specifically using a PowerLawConvectionTerm) couples the
equations together more effectively.
2. Instead of using the factor 1e3 in terms containing sigma in equation 4, we
can move it to equation 3 as a time scale in order to slow down the transient
change of sigma.
Yes, this is what I meant by "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.” A rigorous derivation of the sigma equation might show you the
proportionality that should arise, but more likely sigma will have its own
mobility that you have some freedom to pick.
3. I changed the code to cylindrical coordinates and I am now trying to
identify known solutions (e.g., figure 1 in Campelo 2006). Unfortunately, the
parameter values are not well specified. So far, I was unable to obtain
physically reasonable solutions. If I'll be able to identify a regime that
works I will upload the code here.
That would be interesting to see. I agree that the paper is a little vague on
important details. They say that the three free parameters are \epsilon, a_0,
and V_0, but they never define or use V_0 and they never give any value for
a_0. For that matter, they say that \epsilon is set to the mesh size, but I see
no sign that they ever define that, either.
Thank you again for your great help,
Matan
On Sun, Feb 20, 2022 at 9:23 AM Matan Mussel
<[email protected]<mailto:[email protected]>> wrote:
Hi Jonathan,
Thank you very much, the modifications you made to the code are incredibly
helpful! In particular, the way you rewrote the equation for sigma is something
I did not know how to do alone.
I started examining this version of the code. I am still wondering about the
equations since after ~500dt the circle takes an octagonal shape, and after
~1100dt the solution diverges. I'll try to use parameters with a known
solution. Since the solutions I know of were solved in cylindrical coordinates,
I'll try first to convert the code to work with a cylindrical grid.
Thank you again for your significant help,
Matan
On Fri, Feb 18, 2022 at 4:16 AM 'Guyer, Jonathan E. Dr. (Fed)' via fipy
<[email protected]<mailto:[email protected]>> wrote:
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%7Cc07d57c3579d40c16b1a08d9f48726a1%7C2ab5d82fd8fa4797a93e054655c61dec%7C1%7C0%7C637809683582655874%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000&sdata=JdrLEp84TWV%2FSeuiGL8Ej0SzKMr%2Fd%2BirtmFtKY1XcmI%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]<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].
