As you point out, it’s not a simple problem, and no, it’s not enough to use the compressibility from the “multicomponent compressible” model. We discuss the whole problem in more detail here https://doi.org/10.1093/gji/ggaa078, specifically, we also describe an example concerning phase transitions. The parameter files we used can be found in this repo GitHub - gassmoeller/formulations-of-compressible-mantle-convection-data, and the one in the folder slab-dehydration-reaction is the one with the phase transition.
Specifically, the problem is that even if you use a compressible model, the mass conservation equation needs to compute the density gradient somehow. Most implemetations in compressible mantle convection codes assume that the density follows some sort of reference profile (and we also have this option: The Anelastic Liquid Approximation). However, that means that the mass only changes in the correct way if you are on this reference profile (in terms of temperature, composition etc.). So if you whole model is just basalt/eclogite, that might be okay, but if you have other compositions, this will not work. So what we do instead is that we have the option to interpolate the density onto a finite-element field and then compute the density gradient from that, which always takes into account the actual density in the model. To use that option, you need to pick
set Mass conservation = projected density field
set Temperature equation = real density
in your input file. You also need to have a compositional field with the name projected_density that is used for the purpose of interpolating the density, for example:
subsection Compositional fields
set Number of fields = 2
set Names of fields = crust, projected_density
set Compositional field methods = field, prescribed field
But for this to work, the material model needs to actually provide the outputs that can be interpolated onto the finite element field. The visco-plastic model does not do that at the moment (since it was written to be incompressible). It would need to have a function like this:
template <int dim>
NameOfMaterialModel<dim>::create_additional_named_outputs (MaterialModel::MaterialModelOutputs<dim> &out) const
if (out.template get_additional_output<PrescribedFieldOutputs<dim>>() == NULL)
const unsigned int n_points = out.n_evaluation_points();
(new MaterialModel::PrescribedFieldOutputs<dim> (n_points, this->n_compositional_fields())));
that creates the outputs, and then fill them in the evaluate function:
if (PrescribedFieldOutputs<dim> *prescribed_field_out = out.template get_additional_output<PrescribedFieldOutputs<dim>>())
for (unsigned int i=0; i < in.n_evaluation_points(); ++i)
prescribed_field_out->prescribed_field_outputs[i][density_field_index] = out.densities[i];
For an exmaple, you can have a look at benchmarks/compressibility_formulations/plugins/compressibility_formumations.cc, and the steinberger material model has this as well. So you would have to add this to the material model you want to use.
For the latent heat, the material model needs to compute terms that describe how much the energy should change, and I do not think that the viscoplastic/multicomponent compressible model computes these (but I could be wrong, so maybe someone who actually uses these models can correct me). The latent heat material model is one that computes them (it fills the variables out.entropy_derivative_pressure and out.entropy_derivative_temperature), and you can then use the “latent heat” heating model to actually include the terms in the energy equation. But the latent heat material model only allows for 2 compositions, so that may not be exactly what you want either.
If you don’t want to make all of the required changes in the code, one option would also be to use the steinberger material model as the one that computes density/thermal expansivity/specific heat/compressibility when you are compositing material models. It has all of the required terms implemented, but it requires a table of material properties computed with a mineral physics software for this to work.
I hope that helps a bit, and don’t hesitate to ask if something is unclear, this is quite a lot of information and I hope I didn’t forget anything. I have only been using this feature for applications in the sublithospheric mantle, so I’ve only added the functionality to the material models designed for those sort of applications.