[WIP] new melt formulation

[jdannberg]

This is a first try to change our melt transport formulation to the one from https://doi.org/10.1029/2023JB027469
@tjhei Can you have a look? I used the melt_transport_convergence_test to test the implementation, but the results are quite wrong (the velocities a little bit, the pressures by a lot). For now, I also hard-coded e to be 1.e-5 / viscosity_c; in the paper they use 1.e-10 / viscosity_c as far as I can tell, but if I make e any smaller the solver does not converge any more. That is even though this problem never goes to small melt fractions, so we should not even need any regularization. And for the testcase, I thought that even if I made a mistake in the computation of the pressure somehow, at least the solid velocity should still be correct, which it isn't.

EDIT: We now think the solution is correct, based on the fact that we get the correct velocities and fluid pressure for the melt_transport_convergence_simple test. The other pressures are off by a constant factor, probably due to the pressure normalization.

Now the question is: Is this new formulation actually better than what we have?
One issue is that we cannot normalize the fluid pressure any more, because it explicitly shows up in the equations. In addition, the fluid pressure is defined to be zero in regions without melt, which makes it jump at the melt/no melt boundary. Therefore, it might be better to apply the regularization to both p_f and p_t in the second equation. This would make the fluid pressure equal to the total pressure for zero porosity (no jump) and we can again normalize the pressure because only p_t - p_f appears in the equations, so we could change both by the same amount without changing the solution otherwise.

This still means that in the case of zero porosity, the second equation would be eps (p_t - p_f) = 0 with a small eps (in the paper, they choose 1e-10). We need to see if that is a problem for the linear solver.

Other things we still need to do:

• Try with regularization -- we need the additional regularization (whole p_f row), otherwise the fluid pressure values have a really big spike at the interface between melt and not melt regions
• Pressure normalization (apply to both total and fluid pressure: let p_t be zero at the surface and apply the normalization to p_f as well)
• change porosity and other advection equations
• remove melt_cell
• rename solution variables
• Compute inverse compaction pressure rather than compaction pressure in the material models
• Update the cookbooks/tests accordingly: Some material models had minimum porosities for the compaction viscosity; these should now instead be translated into how large the regularization is. For reference: Reactive fluid transport: 1e-8. Katz 2003: 1e-4. Shear bands: 1e-4. Solubility: 1e-8.
• Think about scaling for regularization (it should scale with some reference viscosity, but we cannot use the actual compaction viscosity because that becomes infinity at some point)
• do not average divergence of velocity for melt advection

Test cases to do:

• Convergence test for the Arbogast with both May and our new regularization
• Check if solid pressure is correct
• Solitary wave
• Time-dependent global 2D box
• Update the analytical solutions for all p_c cases + recompute the RHS terms for the manufactured solutions for p_c

Already done:

• think about the preconditioner
• Assemble new terms
• remove K_D scaling (p_c_scale)
• recompute other melt outputs: Melt velocity is already correct, solid pressure has been changed
• coupling
• we checked that the boundary conditions should be okay as they are

For new features/models or changes of existing features:

• I have tested my new feature locally to ensure it is correct.
• I have added a changelog entry in the doc/modules/changes directory that will inform other users of my change.

[tjhei]

we will need to rename scratch.phi_p_c ...

[tjhei]

Is the 1/ correct?

[tjhei]

can you multiply with regularization below and remove the outer 1/ ?

[jdannberg]

Line 230:
replace one_over_eta with one_over_viscosity_c_e

[jdannberg]

Maybe 1./viscosity_c needs to be -1./viscosity_c
Edit: Nevermind, it should stay a plus

[jdannberg]

times (-1) and this should be one_over_viscosity_c

[jdannberg]

times (-1) and this should be one_over_viscosity_c

[jdannberg]

Add the e_squared here.