diff --git a/doc/modules/changes/20260713_lhy11009 b/doc/modules/changes/20260713_lhy11009
new file mode 100644
index 00000000000..8d7a65bafda
--- /dev/null
+++ b/doc/modules/changes/20260713_lhy11009
@@ -0,0 +1,6 @@
+Added: A new composition type of "phase kinetics"
+which is used to modify composition properties in
+visco plastic material models base on phase transition
+kinetics.
+
+(Haoyuan Li, 2026/07/13)
diff --git a/doc/sphinx/parameters/Compositional_20fields.md b/doc/sphinx/parameters/Compositional_20fields.md
index 64cb47d1e64..cbfe52171d6 100644
--- a/doc/sphinx/parameters/Compositional_20fields.md
+++ b/doc/sphinx/parameters/Compositional_20fields.md
@@ -65,7 +65,7 @@ The format of valid entries for this parameter is that of a map given as “
:name: parameters:Compositional_20fields/Types_20of_20fields
**Default value:** unspecified
-**Pattern:** [List of <[Selection chemical composition|stress|strain|grain size|porosity|density|entropy|generic|unspecified ]> of length 0...4294967295 (inclusive)]
+**Pattern:** [List of <[Selection chemical composition|stress|strain|grain size|porosity|density|entropy|phase kinetics|generic|unspecified ]> of length 0...4294967295 (inclusive)]
**Documentation:** A comma separated list denoting a “type” for each of the compositional fields requested. ASPECT uses these types to determine how fields are handled when evaluating the material model and when solving the equations as described below.
@@ -76,7 +76,7 @@ Each entry of the list must be one of several recognized types: * “chemica
* “porosity”: This type of field represents porosity in a two-phase flow or Darcy flow system. Note that setting the type of a compositional field to “porosity” does not automatically enable melt transport, which is done with the parameter “Melt settings/Include melt transport”.
* “density”: This type of field is a finite-element field representation of the density in the model. This field type is not usually used except for the projected density approximation of the compressible Stokes equations, which uses this field type to compute gradients and time-derivatives of the density.
* “entropy”: This type of field represents entropy. If one or more entropy fields are found in a model, they automatically replace temperature as the main thermodynamic state variable in the model. The temperature equation is then automatically changed to a pure diffusion equation, which is coupled to the entropy advection equation as described in the paper {cite}`dannberg:etal:2022`.
-* “generic”: The generic type is intended to be a placeholder type that is not used by any component of ASPECT unless in user-provided source code.
+* “phase kinetics”: This type of field represents the progress of a phase transition controlled by reaction kinetics. It will only be considered in material models that include models for time-dependent phase transition kinetics.* “generic”: The generic type is intended to be a placeholder type that is not used by any component of ASPECT unless in user-provided source code.
* “unspecified”: The unspecified type is intended to tell ASPECT that the user has not explicitly indicated the type of this field. ASPECT will then try to detect the type automatically based on the name, but will default to “chemical composition” if the name does not correspond to a known type.
Note that while ASPECT’s functionality can make use of the field types, not all of the code will make use of it. It is the user’s responsibility to check that the chosen material model and other plugins interpret the compositional fields as intended.
diff --git a/doc/sphinx/parameters/Material_20model.md b/doc/sphinx/parameters/Material_20model.md
index c384bf74150..9b13956c331 100644
--- a/doc/sphinx/parameters/Material_20model.md
+++ b/doc/sphinx/parameters/Material_20model.md
@@ -1539,6 +1539,15 @@ Units: \si{\pascal\second}
**Documentation:** A list of depths where phase transitions occur. Values must monotonically increase. Units: \si{\meter}.
::::
+::::{dropdown} __Parameter:__ {ref}`Phase transition kinetics mapping`
+:name: parameters:Material_20model/Grain_20size_20model/Phase_20transition_20kinetics_20mapping
+**Default value:** -1
+
+**Pattern:** [Anything]
+
+**Documentation:** A list of indices that maps each phase transition to a phase-kinetics compositional field. For example, an entry of 0 indicates that the corresponding phase transition uses the 0th phase-kinetics composition. All following phase transitions will be affected by the metastablity of that transition kinetics. A negative value means the phase transition is assumed to be equilibrium.
+::::
+
::::{dropdown} __Parameter:__ {ref}`Phase transition pressure widths`
:name: parameters:Material_20model/Grain_20size_20model/Phase_20transition_20pressure_20widths
**Default value:**
@@ -1903,6 +1912,15 @@ Units: \si{\pascal\second}
**Documentation:** A list of depths where phase transitions occur. Values must monotonically increase. Units: \si{\meter}.
::::
+::::{dropdown} __Parameter:__ {ref}`Phase transition kinetics mapping`
+:name: parameters:Material_20model/Latent_20heat/Phase_20transition_20kinetics_20mapping
+**Default value:** -1
+
+**Pattern:** [Anything]
+
+**Documentation:** A list of indices that maps each phase transition to a phase-kinetics compositional field. For example, an entry of 0 indicates that the corresponding phase transition uses the 0th phase-kinetics composition. All following phase transitions will be affected by the metastablity of that transition kinetics. A negative value means the phase transition is assumed to be equilibrium.
+::::
+
::::{dropdown} __Parameter:__ {ref}`Phase transition pressure widths`
:name: parameters:Material_20model/Latent_20heat/Phase_20transition_20pressure_20widths
**Default value:**
@@ -3421,6 +3439,15 @@ A typical example would be to set this runtime parameter to ‘pi=3.14159265
**Documentation:** A list of depths where phase transitions occur. Values must monotonically increase. Units: \si{\meter}.
::::
+::::{dropdown} __Parameter:__ {ref}`Phase transition kinetics mapping`
+:name: parameters:Material_20model/Multicomponent_20compressible/Phase_20transition_20kinetics_20mapping
+**Default value:** -1
+
+**Pattern:** [Anything]
+
+**Documentation:** A list of indices that maps each phase transition to a phase-kinetics compositional field. For example, an entry of 0 indicates that the corresponding phase transition uses the 0th phase-kinetics composition. All following phase transitions will be affected by the metastablity of that transition kinetics. A negative value means the phase transition is assumed to be equilibrium.
+::::
+
::::{dropdown} __Parameter:__ {ref}`Phase transition pressure widths`
:name: parameters:Material_20model/Multicomponent_20compressible/Phase_20transition_20pressure_20widths
**Default value:**
@@ -5206,6 +5233,15 @@ Note that melt does not freeze unless the ’Freezing rate’ parameter
**Documentation:** A list of phase indicators in a look-up table for each phase transition. This parameter selectively assign different rheologies to specific phases, rather than having a unique rheology for each phase in the table. For example, if the table has phases 0, 1, and 2, and one only want a distinct rheology for phase 2, then only phase 2 is needed in the list of indicator. And phases 0, 1 will just be assigned the rheology of the base phase.
::::
+::::{dropdown} __Parameter:__ {ref}`Phase transition kinetics mapping`
+:name: parameters:Material_20model/Visco_20Plastic/Phase_20transition_20kinetics_20mapping
+**Default value:** -1
+
+**Pattern:** [Anything]
+
+**Documentation:** A list of indices that maps each phase transition to a phase-kinetics compositional field. For example, an entry of 0 indicates that the corresponding phase transition uses the 0th phase-kinetics composition. All following phase transitions will be affected by the metastablity of that transition kinetics. A negative value means the phase transition is assumed to be equilibrium.
+::::
+
::::{dropdown} __Parameter:__ {ref}`Phase transition pressure widths`
:name: parameters:Material_20model/Visco_20Plastic/Phase_20transition_20pressure_20widths
**Default value:**
diff --git a/include/aspect/introspection.h b/include/aspect/introspection.h
index f6de4afe47b..3de6af3bef0 100644
--- a/include/aspect/introspection.h
+++ b/include/aspect/introspection.h
@@ -63,8 +63,9 @@ namespace aspect
porosity = 4,
density = 5,
entropy = 6,
- generic = 7,
- unspecified = 8
+ phase_kinetics = 7,
+ generic = 8,
+ unspecified = 9
} type;
/**
@@ -94,6 +95,8 @@ namespace aspect
return CompositionalFieldDescription::density;
else if (input == "entropy")
return CompositionalFieldDescription::entropy;
+ else if (input == "phase kinetics")
+ return CompositionalFieldDescription::phase_kinetics;
else if (input == "generic")
return CompositionalFieldDescription::generic;
else if (input == "unspecified")
diff --git a/include/aspect/material_model/equation_of_state/interface.h b/include/aspect/material_model/equation_of_state/interface.h
index dfb7dd55cf3..8fe223967c4 100644
--- a/include/aspect/material_model/equation_of_state/interface.h
+++ b/include/aspect/material_model/equation_of_state/interface.h
@@ -116,6 +116,8 @@ namespace aspect
phase_average_equation_of_state_outputs(const EquationOfStateOutputs &eos_outputs_all_phases,
const std::vector &phase_function_values,
const std::vector &n_phase_transitions_per_composition,
+ const std::vector &phase_kinetics_values,
+ const std::vector &phase_kinetics_mapping,
EquationOfStateOutputs &eos_outputs);
}
}
diff --git a/include/aspect/material_model/utilities.h b/include/aspect/material_model/utilities.h
index 014032af3f3..6cc516ec58b 100644
--- a/include/aspect/material_model/utilities.h
+++ b/include/aspect/material_model/utilities.h
@@ -516,7 +516,8 @@ namespace aspect
* on phase functions and the change of variables on the trajectory of phase boundaries.
* Thus on a single phase boundary, values of variables change gradually from one phase
* to the other. The values of the phase function used to average the properties varies
- * between 0 and 1.
+ * between 0 and 1. In addition, values in @p phase_kinetics_values is used the modify
+ * the results based on phase transition kinetics.
*/
double phase_average_value (const std::vector &phase_function_values,
const std::vector &n_phase_transitions_per_composition,
@@ -524,7 +525,32 @@ namespace aspect
const unsigned int composition_index,
const PhaseUtilities::PhaseAveragingOperation operation = PhaseUtilities::arithmetic);
+ double phase_average_value (const std::vector &phase_function_values,
+ const std::vector &n_phase_transitions_per_composition,
+ const std::vector &phase_kinetics_values,
+ const std::vector &phase_kinetics_mapping,
+ const std::vector ¶meter_values,
+ const unsigned int composition_index,
+ const PhaseUtilities::PhaseAveragingOperation operation = PhaseUtilities::arithmetic);
+ /**
+ * Extract the phase-transition kinetics information associated with a
+ * particular composition. The number of phase transitions for each
+ * composition is given in @p n_phase_transitions_per_composition, while
+ * @p phase_transition_kinetics_mapping contains one entry for every phase
+ * transition across all compositions. A non-negative mapping entry gives
+ * the index of the corresponding value in @p phase_kinetics_values, whereas
+ * a negative entry indicates that the phase transition is not associated
+ * with a kinetics variable.
+ * For the composition specified by @p composition_index, this function
+ * returns a pair containing the relevant phase-kinetics values and their
+ * corresponding local phase-transition indices.
+ */
+ std::pair, std::vector>
+ map_transition_kinetics_to_phases (const std::vector &n_phase_transitions_per_composition,
+ const std::vector &phase_transition_kinetics_mapping,
+ const std::vector &phase_kinetics_values,
+ const unsigned int composition_index);
/**
* A data structure with all inputs for the
@@ -787,6 +813,13 @@ namespace aspect
*/
unsigned int n_phases_over_all_chemical_compositions () const;
+ /**
+ * Return the mapping to the transition kinetics
+ */
+ const std::vector
+ get_transition_kinetics_mapping () const;
+
+
/**
* Return the Clapeyron slope (dp/dT of the transition) for
* phase transition number @p phase_transition_index.
@@ -861,6 +894,7 @@ namespace aspect
std::vector transition_slopes;
std::vector transition_temperature_upper_limits;
std::vector transition_temperature_lower_limits;
+ std::vector transition_kinetics_mapping;
/**
* Whether to define the phase transitions based on depth, or pressure.
diff --git a/source/material_model/equation_of_state/interface.cc b/source/material_model/equation_of_state/interface.cc
index d414d108ca0..83dd09701ed 100644
--- a/source/material_model/equation_of_state/interface.cc
+++ b/source/material_model/equation_of_state/interface.cc
@@ -45,22 +45,36 @@ namespace aspect
phase_average_equation_of_state_outputs(const EquationOfStateOutputs &eos_outputs_all_phases,
const std::vector &phase_function_values,
const std::vector &n_phase_transitions_per_composition,
+ const std::vector &phase_kinetics_values,
+ const std::vector &phase_transition_kinetics_mapping,
EquationOfStateOutputs &eos_outputs)
{
for (unsigned int c=0; c (const EquationOfStateOutputs &, \
const std::vector &phase_function_values, \
const std::vector &n_phase_transitions_per_composition, \
+ const std::vector &phase_kinetics_values, \
+ const std::vector &phase_transition_kinetics_mapping,\
EquationOfStateOutputs &);
ASPECT_INSTANTIATE(INSTANTIATE)
diff --git a/source/material_model/multicomponent_compressible.cc b/source/material_model/multicomponent_compressible.cc
index 2198806329c..32df2b37b76 100644
--- a/source/material_model/multicomponent_compressible.cc
+++ b/source/material_model/multicomponent_compressible.cc
@@ -80,6 +80,8 @@ namespace aspect
phase_average_equation_of_state_outputs(eos_outputs_all_phases,
phase_function_values,
n_phase_transitions_for_each_chemical_composition,
+ {},
+ {},
eos_outputs);
// Calculate volume fractions from mass fractions
diff --git a/source/material_model/utilities.cc b/source/material_model/utilities.cc
index 83a297430d3..3bfb8b6d601 100644
--- a/source/material_model/utilities.cc
+++ b/source/material_model/utilities.cc
@@ -1099,13 +1099,41 @@ namespace aspect
}
+ double
+ phase_average_value (const std::vector &phase_function_values,
+ const std::vector &n_phase_transitions_per_composition,
+ const std::vector ¶meter_values,
+ const unsigned int composition_index,
+ const PhaseUtilities::PhaseAveragingOperation operation)
+ {
+ return phase_average_value(phase_function_values,
+ n_phase_transitions_per_composition,
+ {},
+ {},
+ parameter_values,
+ composition_index,
+ operation);
+ }
+
double phase_average_value (const std::vector &phase_function_values,
const std::vector &n_phase_transitions_per_composition,
+ const std::vector &phase_kinetics_values,
+ const std::vector &phase_kinetics_mapping,
const std::vector ¶meter_values,
const unsigned int composition_index,
const PhaseUtilities::PhaseAveragingOperation operation)
{
+ std::pair, std::vector> composition_kinetics {{}, {}};
+
+ if (phase_kinetics_values.size() > 0)
+ composition_kinetics = map_transition_kinetics_to_phases (n_phase_transitions_per_composition,
+ phase_kinetics_mapping,
+ phase_kinetics_values,
+ composition_index);
+
+ const std::vector &kinetics_values = composition_kinetics.first;
+ const std::vector &kinetics_mapping = composition_kinetics.second;
// Calculate base index and assign base value
unsigned int start_phase_index = 0;
for (unsigned int i=0; i, std::vector>
+ map_transition_kinetics_to_phases (const std::vector &n_phase_transitions_per_composition,
+ const std::vector &phase_transition_kinetics_mapping,
+ const std::vector &phase_kinetics_values,
+ const unsigned int composition_index)
+ {
+ AssertIndexRange(composition_index, n_phase_transitions_per_composition.size());
+
+ unsigned int start_transition_index = 0;
+ for (unsigned int i=0; i>(options.n_values_per_key);
n_phases_total = 0;
diff --git a/source/material_model/visco_plastic.cc b/source/material_model/visco_plastic.cc
index 612842b3bbc..e2be0a21dac 100644
--- a/source/material_model/visco_plastic.cc
+++ b/source/material_model/visco_plastic.cc
@@ -145,10 +145,22 @@ namespace aspect
phase_function_values[j] = phase_function.compute_value(phase_inputs);
}
+ // Include phase transition kinetics
+ const std::vector phase_kinetics_indices =
+ this->introspection().get_indices_for_fields_of_type(CompositionalFieldDescription::phase_kinetics);
+
+ std::vector phase_kinetics_values(phase_kinetics_indices.size(), 0.0);
+ for (unsigned int j=0; j phase_transition_kinetics_mapping = phase_function.get_transition_kinetics_mapping();
+
// Average by value of gamma function to get value of compositions
phase_average_equation_of_state_outputs(eos_outputs_all_phases,
phase_function_values,
n_phase_transitions_for_each_chemical_composition,
+ phase_kinetics_values,
+ phase_transition_kinetics_mapping,
eos_outputs);
const std::vector volume_fractions = MaterialUtilities::compute_only_composition_fractions(in.composition[i],
diff --git a/source/simulator/parameters.cc b/source/simulator/parameters.cc
index cde5893e858..1da15ae95cf 100644
--- a/source/simulator/parameters.cc
+++ b/source/simulator/parameters.cc
@@ -1352,7 +1352,7 @@ namespace aspect
Patterns::List(Patterns::Anything()),
"A user-defined name for each of the compositional fields requested.");
prm.declare_entry ("Types of fields", "unspecified",
- Patterns::List (Patterns::Selection("chemical composition|stress|strain|grain size|porosity|density|entropy|generic|unspecified")),
+ Patterns::List (Patterns::Selection("chemical composition|stress|strain|grain size|porosity|density|entropy|phase kinetics|generic|unspecified")),
"A comma separated list denoting a ``type'' for each of the "
"compositional fields requested. ASPECT uses these types to "
"determine how fields are handled when evaluating the "
@@ -1399,6 +1399,10 @@ namespace aspect
"which is coupled to the entropy advection equation "
"as described in the paper \\cite{dannberg:etal:2022}."
"\n"
+ "* ``phase kinetics'': This type of field represents the progress of a "
+ "phase transition controlled by reaction kinetics. It will only be "
+ "considered in material models that include models for time-dependent "
+ "phase transition kinetics."
"* ``generic'': The generic type is intended to be a placeholder type "
"that is not used by any component of ASPECT unless in user-"
"provided source code."
@@ -2203,6 +2207,8 @@ namespace aspect
x_compositional_field_types[i] = "porosity";
else if (names_of_compositional_fields[i] == "density_field")
x_compositional_field_types[i] = "density";
+ else if (names_of_compositional_fields[i].find("phase_kinetics") != std::string::npos)
+ x_compositional_field_types[i] = "phase kinetics";
else
x_compositional_field_types[i] = "chemical composition";
}
diff --git a/tests/check_compositional_field_names.cc b/tests/check_compositional_field_names.cc
index 497899e8c35..1913512747f 100644
--- a/tests/check_compositional_field_names.cc
+++ b/tests/check_compositional_field_names.cc
@@ -47,6 +47,8 @@ void f(const aspect::SimulatorAccess &simulator_access,
std::cout << c_names[i] << " is of type generic" << std::endl;
if (descriptions[i].type == aspect::CompositionalFieldDescription::stress)
std::cout << c_names[i] << " is of type stress" << std::endl;
+ if (descriptions[i].type == aspect::CompositionalFieldDescription::phase_kinetics)
+ std::cout << c_names[i] << " is of type phase kinetics" << std::endl;
if (descriptions[i].type == aspect::CompositionalFieldDescription::unspecified)
std::cout << c_names[i] << " is of type unspecified" << std::endl;
}
diff --git a/tests/phase_kinetics_composition_metastable.prm b/tests/phase_kinetics_composition_metastable.prm
new file mode 100755
index 00000000000..cf71bab0657
--- /dev/null
+++ b/tests/phase_kinetics_composition_metastable.prm
@@ -0,0 +1,127 @@
+set Dimension = 2
+set Start time = 0
+set End time = 0
+set Pressure normalization = surface
+set Surface pressure = 0
+set Adiabatic surface temperature = 1662
+set Nonlinear solver scheme = no Advection, no Stokes # iterated Advection and Newton Stokes
+set Resume computation = false
+
+subsection Discretization
+ set Composition polynomial degree = 2
+ set Stokes velocity polynomial degree = 2
+ set Temperature polynomial degree = 2
+ set Use discontinuous composition discretization = true
+ subsection Stabilization parameters
+ set Use limiter for discontinuous composition solution = true
+ set Global composition maximum = 1
+ set Global composition minimum = 0
+ end
+end
+
+
+subsection Formulation
+ set Formulation = custom
+ set Mass conservation = incompressible
+ set Temperature equation = reference density profile
+end
+
+# model geometry
+subsection Geometry model
+ set Model name = box
+ subsection Box
+ set X repetitions = 50
+ set Y repetitions = 50
+ set X extent = 800e3
+ set Y extent = 1200e3
+ end
+end
+
+# Mesh refinement specifications
+subsection Mesh refinement
+ set Initial adaptive refinement = 0
+ set Initial global refinement = 0
+ set Time steps between mesh refinement = 0
+end
+
+
+# Temperature boundary and initial conditions maintain
+# a linear increase in temperature along the x axis
+subsection Boundary temperature model
+ set Fixed temperature boundary indicators = left, right
+ set List of model names = box
+ subsection Box
+ set Left temperature = 273.15
+ set Right temperature = 2273.15
+ end
+end
+
+
+subsection Initial temperature model
+ set Model name = function
+ subsection Function
+ set Coordinate system = cartesian
+ set Variable names = x, y
+ set Function constants = XMAX=800e3, Tl=273.15, Tr=2273.15
+ set Function expression = Tl * (x - XMAX)/(-XMAX) + Tr * x / XMAX
+ end
+end
+
+
+# Include crust and barzburgite for chemical compositions
+subsection Compositional fields
+ set Number of fields = 3
+ set Names of fields = crust, harzburgite, phase_kinetics_1
+ set Compositional field methods = field, field, field
+end
+
+
+# Initial composition model
+# Assign the background composition for this test
+subsection Initial composition model
+ set List of model names = function
+ subsection Function
+ set Coordinate system = cartesian
+ set Variable names = x, y
+ set Function expression = 0.0 ; 0.0; ((x < 500e3)? 0.0: 1.0)
+ end
+end
+
+# Phase transition kinetics mapping to map phase transitions and transition kinetics
+subsection Material model
+ set Model name = visco plastic
+ set Material averaging = harmonic average
+ subsection Visco Plastic
+ set Reference strain rate = 1.e-15
+ set Minimum viscosity = 1e18
+ set Maximum viscosity = 1e24
+ set Phase transition temperatures = background:1662.0|1662.0|1662.0|1662.0, crust: 1173.0|1662.0|1662.0, harzburgite: 1662.0|1662.0|1662.0|1662.0
+ set Phase transition depths = background:410e3|520e3|560e3|670e3, crust: 80e3|665e3|720e3, harzburgite: 410e3|520e3|560e3|670e3
+ set Phase transition widths = background:5e3, crust: 5e3, harzburgite: 5e3
+ set Phase transition Clapeyron slopes = background:4e6|4.1e6|4e6|-2e6, crust: 0.0|4e6|1.3e6, harzburgite: 4e6|4.1e6|4e6|-2e6
+ set Thermal diffusivities = 1.0e-6
+ set Heat capacities = 1250.0
+ set Thermal expansivities = 3.1e-5
+ set Densities = background: 3300.0|3394.4|3442.1|3453.2|3617.6,\
+ harzburgite: 3235.0|3372.3|3441.7|3441.7|3680.8,\
+ crust: 3000.0|3540.0|3613.0|3871.7
+ set Phase transition kinetics mapping = background: 0|-1|-1|-1, harzburgite: -1, crust:-1
+ end
+end
+
+subsection Heating model
+ set List of model names = adiabatic heating, latent heat, shear heating
+end
+
+subsection Gravity model
+ set Model name = ascii data
+end
+
+subsection Postprocess
+ set List of postprocessors = visualization, temperature statistics, composition statistics, material statistics
+ subsection Visualization
+ set List of output variables = material properties, heating
+ set Output format = gnuplot
+ set Time between graphical output = 0
+ end
+end
diff --git a/tests/phase_kinetics_composition_metastable/screen-output b/tests/phase_kinetics_composition_metastable/screen-output
new file mode 100644
index 00000000000..9063c02ba3d
--- /dev/null
+++ b/tests/phase_kinetics_composition_metastable/screen-output
@@ -0,0 +1,16 @@
+
+Number of active cells: 2,500 (on 1 levels)
+Number of degrees of freedom: 100,704 (20,402+2,601+10,201+22,500+22,500+22,500)
+
+*** Timestep 0: t=0 years, dt=0 years
+
+ Postprocessing:
+ Writing graphical output: output-phase_kinetics_composition_metastable/solution/solution-00000
+ Temperature min/avg/max: 273.1 K, 1273 K, 2273 K
+ Compositions min/max/mass: 0/0/0 // 0/0/0 // 0/1/3.616e+11
+ Average density / Average viscosity / Total mass: 3430 kg/m^3, 5.453e+23 Pa s, 3.293e+15 kg
+
+Termination requested by criterion: end time
+
+
+
diff --git a/tests/phase_kinetics_composition_metastable/statistics b/tests/phase_kinetics_composition_metastable/statistics
new file mode 100644
index 00000000000..c91ad382223
--- /dev/null
+++ b/tests/phase_kinetics_composition_metastable/statistics
@@ -0,0 +1,26 @@
+# 1: Time step number
+# 2: Time (years)
+# 3: Time step size (years)
+# 4: Number of mesh cells
+# 5: Number of Stokes degrees of freedom
+# 6: Number of temperature degrees of freedom
+# 7: Number of degrees of freedom for all compositions
+# 8: Number of nonlinear iterations
+# 9: Visualization file name
+# 10: Minimal temperature (K)
+# 11: Average temperature (K)
+# 12: Maximal temperature (K)
+# 13: Average nondimensional temperature (K)
+# 14: Minimal value for composition crust
+# 15: Maximal value for composition crust
+# 16: Global mass for composition crust
+# 17: Minimal value for composition harzburgite
+# 18: Maximal value for composition harzburgite
+# 19: Global mass for composition harzburgite
+# 20: Minimal value for composition phase_kinetics_1
+# 21: Maximal value for composition phase_kinetics_1
+# 22: Global mass for composition phase_kinetics_1
+# 23: Average density (kg/m^3)
+# 24: Average viscosity (Pa s)
+# 25: Total mass (kg)
+0 0.000000000000e+00 0.000000000000e+00 2500 23003 10201 67500 0 output-phase_kinetics_composition_metastable/solution/solution-00000 2.73150000e+02 1.27315000e+03 2.27315000e+03 5.00000000e-01 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 3.61600000e+11 3.43044552e+03 5.45342439e+23 3.29322770e+15