geo-fluid-dynamics / sapphire Goto Github PK

From my old notes, there were some surprises that may be worth considering the next time I look at this:

• The pressure error magnitude is proportional to $\tau$
• Increasing $\sigma$ does not reduce pressure error, actually observed a small increase
• Decreasing the Grashof number doesn't substantially affect the pressure error

I copied the code from assigning initial values to create the following for plotting manufactured solutions:

def plot_manufactured_solution(
        Model,
        meshsize = 128,
        parameters = {},
        time = None):
        
    MMSVerificationModel = make_mms_verification_model_class(Model)
    
    model = MMSVerificationModel(meshsize = meshsize)

    model.init_manufactured_solution()
    
    if time is not None:
    
        model.time.assign(time)
    
    for u_m, V in zip(
            model.manufactured_solution, model.function_space):

        model.solution.assign(fe.interpolate(u_m, V))

    model.plot()

but this did not have the expected result; which makes me worry about how initial values are being assigned in the MMS code. By nature the MMS code shouldn't be passing if the initial values aren't being set correctly, so why is the behavior different here?

This happens when solving the convection-coupled phase-change model with low smoothing and a poor initial guess. Maybe I can demonstrate this for the 1D enthalpy phase-change model and ask the Firedrake Slack team.

Currently we use a direct solver, which is already slow in 2D and entirely impractical in 3D. This only gets worse when we want to use our model as a forward model for inverse problems.

Patrick Farrell showed me that applying the following patch allows us to solve the nonlinear problem often with much less regularization:

diff --git a/sapphire/simulations/convection_coupled_phasechange.py b/sapphire/simulations/convection_coupled_phasechange.py
index f043774..c31fb6c 100644
--- a/sapphire/simulations/convection_coupled_phasechange.py
+++ b/sapphire/simulations/convection_coupled_phasechange.py
@@ -274,6 +274,9 @@ class Simulation(sapphire.simulation.Simulation):
                 "snes_max_it": 24,
                 "snes_monitor": None,
                 "snes_converged_reason": False,
+                "snes_linesearch_type": "l2",
+                "snes_linesearch_maxstep": 1.0,
+                "snes_linesearch_damping": 0.9,
                 "ksp_type": "preonly", 
                 "pc_type": "lu", 
                 "mat_type": "aij",

This worked fine for the octadecane melting regression test. I want to try it with even smaller $\sigma$, both for octadecane melting and water freezing regression tests, and then for higher fidelity simulations to the simulated times where we compare to benchmarks in our paper.

There is a still a comment in the docstring about the pressure accuracy not being verified, which is no longer true. That was fixed in #89 .

Whoops! Second order time accuracy never would have worked with this.

phil_t = bdf2( phil(T_np1, Cl_np1), phil(T_n, Cl_n), phil(T_np1, Cl_np1))

But don't commit any of the docs to the repository. That creates quite a lot of noise and increases the size of the repo rapidly.

Before starting this, we should verify the binary alloy enthalpy model (issue #27) and the enthalpy-porosity model (issue #5).

Based on this paper: https://www.researchgate.net/publication/226804783_Heat_transfer_during_the_melting_of_ice_around_a_horizontal_isothermal_cylinder

The current enthalpy-porosity model should be enough for a nice qualitative result. A nice quantitative result, and perhaps easier set up in general, would first require phase-dependent heat capacity and thermal conductivity (for which we have issue #37).

This would be our first demonstration of an interesting geometry.

And leave the water freezing benchmark to the Sapphire-Ice portfolio.

• Replace portfolio/benchmarks section with examples section
  - [ ] Scale streamline thickness with max speed from entire simulation
  - [ ] Update octadecane melting gif
• Add water freezing gif
• Add gallium melting gif
• Add freezing salt water from above (with brine plume) once that simulation is in the master branch. Use matplotlib.cm.ocean colormap if showing porosity

https://pandas.pydata.org/

https://docs.python-guide.org/writing/gotchas/#mutable-default-arguments

For the most part I did an okay job of using None as the default argument value and then setting the mutable value in __init__; but there are many cases where I use a dict which is mutable and can lead to all sorts of unintended behavior. For example, the base class Simulation has a dict as the default value for solver_parameters.

Have the plotter class get all of its data from HDF5 outputs, i.e. there should be no need to instantiate the original Model.

e.g. I just copy/pasted the benchmark code in commit fe878b8 .

I'm trying this on the following branch: https://github.com/alexanderzimmerman/fempy/tree/model/heat2

I'm currently running into an "OSError: ... undefined symbol", which has come up in this probably unrelated Firedrake issue: firedrakeproject/firedrake#581

I opened a discussion about using time as a variable on the Firedrake Slack channel. I need to provide them a better MWE showing my requirements for the time variable, e.g. an actual MMS example.

e.g. perhaps we should be using GLS instead of the penalty method.

First, reproduce this result from Phaseflow (See geo-fluid-dynamics/phaseflow-fenics#83):

I have been trying to resolve a sharper interface, and it's making it so that the freezing front is stuck:

Before digging into that, I should make sure we get the first result with the oversmoothed interface. Also note that, for the old Phaseflow result, not as much work was put into setting up the example with physically relevant parameters. While the movement of the interface is qualitatively wrong for the latest FemPy result, the magnitude of the concentration at the solidification front is about right.

@Article{michalek2003simulations,
title={Simulations of the water freezing process--numerical benchmarks},
author={MICHA{\L}EK, TOMASZ and KOWALEWSKI, TOMASZ A},
journal={Task Quarterly},
volume={7},
number={3},
pages={389--408},
year={2003}
}

Our current approach to regularizing the convection-coupled phase-change nonlinear problem works very well for the benchmarks we simulate. In short, we solve for small $sigma$ via first solving a sequence of over-regularized (larger $sigma$) problems, which essentially improve the initial guess for Newton's method. Similar approaches are referred to as "continuation".

Patrick Farrell, who develops UFL/FEniCS/Firedrake, said he has been working on a general regularization method for problems with local nonlinearities/irregularities, and that it might work well for our convection-coupled phase-change problem.

So now I've made this branch and added a test to try his idea, or perhaps to try other regularization methods.

The existing regression test (here) for melting octadecane passes. I ran it with

$ python3 -m pytest -v -s -k "validate and melt_octadecane and (not without_auto_smoothing)" sapphire/tests/ > log__test_melt_octadecane_with_auto_smoothing.txt

which output log__test_melt_octadecane_with_auto_smoothing.txt

I added a new test to the same script which uses a simple default solver

class OctadecaneMeltingSimulationWithoutAutoSmoothing(
        OctadecaneMeltingSimulation):
    
    def solve(self, *args, **kwargs):
    
        return sapphire.simulation.Simulation.solve(
            self,
            *args,
            parameters = {
                "snes_type": "newtonls",
                "snes_max_it": 50,
                "snes_monitor": True,
                "snes_converged_reason": True,
                "ksp_type": "preonly", 
                "pc_type": "lu", 
                "mat_type": "aij",
                "pc_factor_mat_solver_type": "mumps"},
            **kwargs)
            
    def run(self, *args, **kwargs):
        
        return sapphire.benchmarks.melt_octadecane_in_cavity.Simulation.run(
            self,
            *args,
            solve = self.solve,
            **kwargs)

We currently expect the test to fail, which it does. I run it with

$ python3 -m pytest -v -s -k "validate and melt_octadecane and without_auto_smoothing" sapphire/tests/ > log__test_melt_octadecane_without_auto_smoothing.txt

which output log__test_melt_octadecane_without_auto_smoothing.txt

Ideally there would be a set of solver parameters, which we could set in the definition of solve above, for which the test passes. With that, we would want to try the same for the water freezing benchmark, and for even smaller $sigma$, e.g. on the order of 0.001 which has been commonly used for some of my recent results.

e.g. to what extent do we depend on VTK? I personally had trouble getting Python VTK set up on my Ubuntu 18.04 system. For now I should at least document how I got VTK working on my system.

Right now we assume constant thermal conductivity and constant heat capacity. For water-ice, these vary significantly between the solid and liquid phases.

I asked how to do this on Slack and it seems there's a way for us to redefine a SNES Monitor:

"
Julian Andrej [3:54 PM]
"I'd implement a SNES Monitor for that
You can then access the jacobian after every nonlinear iteration
so if you have

solver = NonlinearVariationalSolver(...)
mysnes = solver.snes

def monitor(snes, its, fgnorm):
  myjacobian = snes.jacobian

snes.setMonitor(monitor)

myjacobian then contains a Petsc Mat
"

Right now we assume constant thermal conductivity and constant heat capacity. For water-ice, these vary significantly between the solid and liquid phases.

First one should look at the enthalpy model alone (issue #36).

See existing work on Rayleigh-Bernard convection: https://www.firedrakeproject.org/demos/rayleigh-benard.py.html

After correcting the scaling of energy diffusion with non-unit Reynolds number in #94 , the spatial accuracy MMS verification for unsteady Navier-Stokes-Boussinesq didn't require nearly as large of a grid (when using Re = 20). Maybe there are other tests which can be made less expensive now.

In this case, we'd only use Firedrake for its key capabilities, which include assembling system matrices (with Dirichlet BC's) and symbolically deriving Jacobians. We would try ignoring the PETSc interface component of Firedrake and see how far we can get with scipy.

In tests/verification/test__unsteady_navier_stokes_boussinesq.py::test__verify_first_order_temporal_convergence_via_mms, after generating

Delta_t, errors, temporal_orders
0.5, [None, 0.011107126169912135, 0.025971296457037198], None
0.25, [None, 0.006439718841114489, 0.014279010781708205], [None, 0.7864159791056798, 0.8630219973726072]
0.125, [None, 0.003462654871468789, 0.007475212830499297], [None, 0.8951191025760915, 0.9337094754111044]

the nonlinear solver diverges with Delta_t = 0.06125 after solving at t = 0.8125.

https://scikits-odes.readthedocs.io/en/latest/solvers.html

https://computation.llnl.gov/projects/sundials/ida

So far I've been estimating the time by looking at the earliest and latest modification times in the output files. It would be more convenient to have this in the report.

Some old notes, running binary alloy model:

parameters = {
"temperature_rayleigh_number": 8.,
"concentration_rayleigh_number": -9.,
"prandtl_number": 7.,
"stefan_number": 0.2,
"schmidt_number": 6.,
"pure_liquidus_temperature": 0.,
"liquidus_slope": -0.11,
"phase_interface_smoothing": 1./16.,
"autosmooth_maxcount": 16}

m, Delta_t, smin
4, 1./4., 0.0870
4, 1./8., 0.1519
4, 1./16., 0.0758
4, 1./32., 0.0917
4, 1./64., 0.0640
4, 1./128., 0.0625
4, 1./256., 0.0661
4, 1./512., 0.0635
8, 1./64., 0.0804
16, 1./64., 0.0842
32, 1./64., 0.0867

I also have many results with different parameters from running Phaseflow.

Right now we're only using fully implicit backward difference formulas of first, second, and third order. Thankfully, the time discretization of any of our models is well encapsulated, so we can test mixed explicit-implicit methods such as Crank Nicholson and Runge-Kutta. @tomoleary said that Crank Nicholson and RK4 are needed for plugging our models into inverse problems because they are self-adjoint.

Right now, solutions are only written to VTK files. For post-processing, this is only accurate for CG1 bases. For restarting simulations, this is mostly useless.

This was easy with FEniCS, because they provide a Docker image; but how do we do this with Firedrake?

Unfortunately pathlib isn't so seemless for Python 3.5: https://stackoverflow.com/questions/42694112/when-using-pathlib-getting-error-typeerror-invalid-file-posixpathexample-t

I did some development in Python 3.6 and now that I'm running with Python 3.5, pathlib is having some trouble.

We could consider requiring Python 3.6 for a bunch of reasons; but Firedrake supports Python 3.5+ so we should probably do the same.

If it shows up with help, it should have a docstring.

Right now, solutions are only written to VTK files. For post-processing, this is only accurate for CG1 bases. For restarting simulations, this is mostly useless.

sapphire.Simulation should have a simplified interface to firedrake.DumbCheckpoint for checkpointing and restarting.

https://docs.python-guide.org/writing/logging/