[0] Matrix is missing diagonal entry ###

[Matan Mussel]

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

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

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?

--
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.
<CurvatureModelEquations.pdf><MinimalCurvatureModel.py>

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

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 <mata...@gmail.com> 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

<phaseFieldCurvatureModel.pdf><CH_0.0025.png><CH_0.0000.png>

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

Matan -

On Feb 14, 2022, at 3:57 AM, Matan Mussel <mata...@gmail.com> 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).

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

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

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).

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

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). 

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

On Feb 20, 2022, at 10:38 AM, Matan Mussel <mata...@gmail.com> 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 <mata...@gmail.com> 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

 

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

I’m glad I could help.

On Feb 20, 2022, at 2:23 AM, Matan Mussel <mata...@gmail.com> 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.

The chain rule is our friend!

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.  

I haven’t observed either, but maybe I didn’t run long enough. The octagonal symmetry doesn’t surprise me; it’s almost certainly “feeling” the edges of the mesh.

Thank you again for your significant help,

Matan

 

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

On Feb 21, 2022, at 7:48 PM, 'Guyer, Jonathan E. Dr. (Fed)' via fipy <fi...@list.nist.gov> wrote:

On Feb 20, 2022, at 2:23 AM, Matan Mussel <mata...@gmail.com> wrote:

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.  

I haven’t observed either, but maybe I didn’t run long enough. The octagonal symmetry doesn’t surprise me; it’s almost certainly “feeling” the edges of the mesh.

OK, I ran longer and see what you mean. I still think the octagonal symmetry is arising because of the box size. As to the divergence, I’m only guessing, but there are three length scales in your equations:

- one associated with the terms 2nd-order in \phi. This is the \tanh solution that describes the diffuse interface. This length is of order \epsilon.

- one associated with the 4th-order term in \phi. I’m not positive, but I think this one may be order \epsilon, too.

- one associated with the term(s) in sigma and phi. It’s this one that I’m guessing is the issue. As the problem destabilizes, you can see undulations in \xi along the interface as sigma gets large.

My guess is that one or more of these length scales is under-resolved. The paper says that \sigma is small and essentially homogeneous, but they never quantify it.