Reproduce "Melting of PCM" benchmark from danaila2014newton

[agzimmerman]

We've successfully tested the basic capability in issue #14, but it still may be difficult to obtain a solution for the exact parameters from danaila2014newton.

[agzimmerman]

Reducing r to 0.005 (from 0.05) makes the Newton method rapidly diverge.

[agzimmerman]

I was able to converge one time step of the melting octadecane PCM benchmark here: https://github.com/geo-fluid-dynamics/phaseflow-fenics/tree/benchmark/melting_octadecane_pcm

There are some deficiencies:

• It required about 150 Newton iterations, and changing the initial mesh refinement or time step size makes the Newton solver blow up.
• Using a regularization smoothing factor, r, smaller than around 0.05 makes the Newton solver blow up.
• The error is being estimated already below 1.e-4 with the initial grid (refined near the hot wall), and so no adaptive mesh refinement is being performed, and so the error is not being controlled as time goes on. Setting the tolerance to 1.e-5 ends up refining every cell, which then does not converge. Maybe this behavior could be related to the large smoothing factor r.
• A second time step can also converge given enough Newton iterations, but convergence behavior only gets worse (because the mesh is not being adaptively refined).

[agzimmerman]

Maybe relaxing the Newton solver would do the trick, since this has seemed to handle rapid divergence pretty well before. I added issue #154 for this.

[agzimmerman]

So far relaxing the solver isn't helping.

[agzimmerman]

In this lid-driven cavity example, the metric u[0]fenics.dx was not useful, while u[0]u[0]*fenics.dx worked very well. So maybe I just need to get more creative with the metric for this melting case.

[agzimmerman]

I just noticed that when solving the nonlinear variational problem, a Krylov solver is being mentioned, e.g.

Solving nonlinear variational problem. Newton iteration 0: r (abs) = 1.817e+01 (tol = 1.000e-08) r (rel) = 1.000e+00 (tol = 1.000e-04) Solving linear system of size 8244 x 8244 (PETSc LU solver, umfpack). PETSc Krylov solver starting to solve 8244 x 8244 system. Newton iteration 1: r (abs) = 1.182e-02 (tol = 1.000e-08) r (rel) = 6.507e-04 (tol = 1.000e-04) Solving linear system of size 8244 x 8244 (PETSc LU solver, umfpack). PETSc Krylov solver starting to solve 8244 x 8244 system. Newton iteration 2: r (abs) = 1.146e-03 (tol = 1.000e-08) r (rel) = 6.306e-05 (tol = 1.000e-04) Newton solver finished in 2 iterations and 2 linear solver iterations.

This is perplexing. I began exploring this in issue #156 .

Having had no luck pinpointing this problem yet, I think it's about time I try making a minimal script that solves the melting octadecane PCM benchmark (without Phaseflow). If the new script also somehow invokes a Krylov method, then maybe that will get me closer to a MWE that reproduces the issue.

[agzimmerman]

At this point, I think I should verify the coupled problem with MMS (see issue #158) before getting too worried about the benchmark case not working. Maybe my implementation is wrong.

[agzimmerman]

As of PR #172 , we have a pretty good result for the octadecane PCM melting benchmark. The PCI doesn't advance quite far enough. We still need to study the sensitivity of parameters such as the time step size, the initial mesh, the adaptive solver tolerance, etc.