diff --git a/data/material-model/entropy-table/pyrtable/material_table_temperature_pressure_small.txt b/data/material-model/entropy-table/pyrtable/material_table_temperature_pressure_small.txt
new file mode 100644
index 00000000000..e2f674c7fc8
--- /dev/null
+++ b/data/material-model/entropy-table/pyrtable/material_table_temperature_pressure_small.txt
@@ -0,0 +1,8 @@
+# This is a data output from HeFESTo
+# Independent variables are T(K) and P(bar)
+# POINTS: 2 2
+T(K) P(bar) s,J/K/kg rho,kg/m3 alpha,1/K cp,J/K/kg vp,km/s vs,km/s h,J/kg dominant_phase_index,1
+250.0 0.0 509.74740059944844 3332.8967443100287 1.9781287327429287e-05 748.5636509347647 8.116376883289359 4.7495273882440685 -13445454.390915 0.0
+4000.0 0.0 3804.799798825367 2879.31713365359 0.00012376961201237974 1865.6003296108431 5.89690417370827 3.068046128481904 -8401781.805825338 0.0
+250.0 400781.25 326.1989261514726 4738.970221094591 9.15389351114028e-06 565.4531513648064 11.974450477467947 6.779618381631337 -3722250.1404551915 1.0
+4000.0 400781.25 3394.028743328079 4373.768154238075 2.6303932404721963e-05 1362.8625515701538 10.986106548987486 5.8231315598976465 873079.6281181354 1.0
diff --git a/doc/modules/changes/20241107_lhy11009 b/doc/modules/changes/20241107_lhy11009
new file mode 100644
index 00000000000..d09500f84e4
--- /dev/null
+++ b/doc/modules/changes/20241107_lhy11009
@@ -0,0 +1,6 @@
+Added: there is now a new class of phase function that handles discrete phase transitions
+by looking up the most dominant phases in a lookup table. This function can be used
+to make the rheology of the visco-plastic material model dependent on the dominant
+mineral phase.
+
+(Haoyuan Li, 2024/11/07)
diff --git a/include/aspect/material_model/utilities.h b/include/aspect/material_model/utilities.h
index 82c0ff75bdb..fb48807a162 100644
--- a/include/aspect/material_model/utilities.h
+++ b/include/aspect/material_model/utilities.h
@@ -549,6 +549,191 @@ namespace aspect
unsigned int phase_index;
};
+ /**
+ * A class that bundles functionality to look up the dominant phase in
+ * tables for each respective composition and export the values
+ * of phase functions. The class can handle arbitrary numbers of
+ * dominant phases for each composition, but the calling side
+ * has to determine how to use the return values of this object
+ * (e.g. in terms of density or viscosity).
+ */
+ template
+ class PhaseFunctionDiscrete: public ::aspect::SimulatorAccess
+ {
+ public:
+
+ /**
+ * The initialization process loads the contents of the material files
+ * for the respective compositions.
+ */
+ void initialize();
+
+ /**
+ * Percentage of material that has already undergone the phase
+ * transition to the higher-pressure material. For this class
+ * this function only returns 1.0 or 0.0, depending on whether
+ * the selected phase transition has been crossed or not.
+ */
+ double compute_value (const PhaseFunctionInputs &in) const;
+
+ /**
+ * No valid implementation exists for this function, as the derivative of a
+ * discrete function is undefined at locations of phase jumps. This function raises an
+ * error to ensure that a phase derivative request is not made for this phase function.
+ */
+ double compute_derivative () const;
+
+ /**
+ * Return the total number of phase transitions.
+ */
+ unsigned int n_phase_transitions () const;
+
+ /**
+ * Return the total number of phases.
+ */
+ unsigned int n_phases () const;
+
+ /**
+ * Return the total number of phases over all chemical compositions.
+ */
+ unsigned int n_phases_over_all_chemical_compositions () const;
+
+ /**
+ * Return how many phase transitions there are for each chemical composition.
+ */
+ const std::vector &
+ n_phase_transitions_for_each_chemical_composition () const;
+
+ /**
+ * Return how many phases there are for each chemical composition.
+ */
+ const std::vector &
+ n_phases_for_each_chemical_composition () const;
+
+ /**
+ * Return how many phase transitions there are for each composition.
+ * Note, that most likely you only need the number of phase transitions
+ * for each chemical composition, so use the function above instead.
+ * This function is only kept for backward compatibility.
+ */
+ const std::vector &
+ n_phase_transitions_for_each_composition () const;
+
+ /**
+ * Return how many phases there are for each composition.
+ * Note, that most likely you only need the number of phase transitions
+ * for each chemical composition, so use the function above instead.
+ * This function is only kept for backward compatibility.
+ */
+ const std::vector &
+ n_phases_for_each_composition () const;
+
+ /**
+ * Declare the parameters this class takes through input files.
+ * Note that this class does not declare its own subsection,
+ * i.e. the parameters will be declared in the subsection that
+ * was active before calling this function.
+ */
+ static
+ void
+ declare_parameters (ParameterHandler &prm);
+
+ /**
+ * Read the parameters this class declares from the parameter file.
+ * Note that this class does not declare its own subsection,
+ * i.e. the parameters will be parsed from the subsection that
+ * was active before calling this function.
+ */
+ void
+ parse_parameters (ParameterHandler &prm);
+
+
+ private:
+ /**
+ * Directory path where data files are stored.
+ */
+ std::string data_directory;
+
+ /**
+ * List of file names containing material data for each composition.
+ */
+ std::vector material_file_names;
+
+ /**
+ * Minimum temperature values for each composition in the P-T table.
+ */
+ std::vector minimum_temperature;
+
+ /**
+ * Maximum temperature values for each composition in the P-T table.
+ */
+ std::vector maximum_temperature;
+
+ /**
+ * Temperature intervals used for each composition in the P-T table.
+ */
+ std::vector interval_temperature;
+
+ /**
+ * Minimum pressure values for each composition in the P-T table.
+ */
+ std::vector minimum_pressure;
+
+ /**
+ * Maximum pressure values for each composition in the P-T table.
+ */
+ std::vector maximum_pressure;
+
+ /**
+ * Pressure intervals used for each composition in the P-T table.
+ */
+ std::vector interval_pressure;
+
+ /**
+ * List of pointers to objects that read and process data we get from
+ * material data files. There is one pointer/object per lookup file.
+ */
+ std::vector>> material_lookup;
+
+ /**
+ * List of phase indicators of the most dominant phases in the material data files
+ * to construct the different phase transitions in this class. For a description of
+ * the use of the phase indicators, please see the documentation of the input parameter
+ * 'Phase transition indicators' in the function declare_parameters().
+ */
+ std::vector transition_indicators;
+
+ /**
+ * A vector that stores how many phase transitions there are for each compositional field.
+ */
+ std::unique_ptr> n_phase_transitions_per_composition;
+
+ /**
+ * A vector that stores how many phases there are for each compositional field.
+ */
+ std::vector n_phases_per_composition;
+
+ /**
+ * A vector that stores how many phase transitions there are for each chemical compositional field.
+ */
+ std::vector n_phase_transitions_per_chemical_composition;
+
+ /**
+ * A vector that stores how many phases there are for each chemical compositional field.
+ */
+ std::vector n_phases_per_chemical_composition;
+
+ /**
+ * Total number of phases over all compositional fields
+ */
+ unsigned int n_phases_total;
+
+ /**
+ * Total number of phases over all compositional fields
+ */
+ unsigned int n_phases_total_chemical_compositions;
+ };
+
/**
* A class that bundles functionality to compute the values and
* derivatives of phase functions. The class can handle arbitrary
diff --git a/include/aspect/material_model/visco_plastic.h b/include/aspect/material_model/visco_plastic.h
index cd9666d5fc2..446ce21f57b 100644
--- a/include/aspect/material_model/visco_plastic.h
+++ b/include/aspect/material_model/visco_plastic.h
@@ -181,6 +181,12 @@ namespace aspect
class ViscoPlastic : public MaterialModel::Interface, public ::aspect::SimulatorAccess
{
public:
+ /**
+ * Initialization function. Loads the material data and sets up
+ * pointers if it is required.
+ */
+ void
+ initialize () override;
void evaluate(const MaterialModel::MaterialModelInputs &in,
MaterialModel::MaterialModelOutputs &out) const override;
@@ -261,6 +267,16 @@ namespace aspect
*/
MaterialUtilities::PhaseFunction phase_function;
+ /**
+ * Determines whether to look up the dominant phases for each composition in its respective lookup table.
+ */
+ bool use_dominant_phase_for_viscosity;
+
+ /**
+ * Object that handles discrete phase transitions for the rheology if requested by the variable use_dominant_phase_for_viscosity.
+ */
+ std::unique_ptr> phase_function_discrete;
+
};
}
diff --git a/include/aspect/structured_data.h b/include/aspect/structured_data.h
index 121e6d821aa..341a871187a 100644
--- a/include/aspect/structured_data.h
+++ b/include/aspect/structured_data.h
@@ -248,6 +248,15 @@ namespace aspect
*/
double get_maximum_component_value(const unsigned int component) const;
+ /**
+ * Retrieve the number of table points for a given dimension.
+ * Equivalent to calling get_interpolation_point_coordinates().size().
+ *
+ * @param dimension The index of the dimension for which to get the number of table points.
+ * @return The number of points along the specified dimension.
+ */
+ unsigned int get_number_of_coordinates(const unsigned int dimension) const;
+
private:
/**
* The number of data components read in (=columns in the data file).
diff --git a/source/material_model/utilities.cc b/source/material_model/utilities.cc
index 5c686673e83..0b8c1c3dfd0 100644
--- a/source/material_model/utilities.cc
+++ b/source/material_model/utilities.cc
@@ -1142,6 +1142,8 @@ namespace aspect
return averaged_parameter;
}
+
+
template
PhaseFunctionInputs::PhaseFunctionInputs(const double temperature_,
const double pressure_,
@@ -1159,6 +1161,265 @@ namespace aspect
+ template
+ void
+ PhaseFunctionDiscrete::initialize()
+ {
+ AssertThrow (this->introspection().get_number_of_fields_of_type(CompositionalFieldDescription::chemical_composition)+1 == material_file_names.size(),
+ ExcMessage(" The 'discrete phase function' requires that the number of compositional fields of type chemical composition plus one (for a background field) "
+ "matches the number of lookup tables."));
+
+
+ minimum_temperature = std::vector(material_file_names.size());
+ maximum_temperature = std::vector(material_file_names.size());
+ interval_temperature = std::vector(material_file_names.size());
+ minimum_pressure = std::vector(material_file_names.size());
+ maximum_pressure = std::vector(material_file_names.size());
+ interval_pressure = std::vector(material_file_names.size());
+ for (unsigned int i = 0; i < material_file_names.size(); ++i)
+ {
+ material_lookup
+ .push_back(std::make_unique>(8,1.0));
+
+ material_lookup[i]->load_file(data_directory+material_file_names[i],
+ this->get_mpi_communicator());
+
+ Assert(material_lookup[i]->has_equidistant_coordinates(), ExcMessage("The loaded lookup tables do not use equidistant coordinates. The class 'PhaseFunctionDiscrete' cannot currently handle non equidistant coordinates."));
+
+ minimum_temperature[i] = material_lookup[i]->get_interpolation_point_coordinates(0).front();
+ maximum_temperature[i] = material_lookup[i]->get_interpolation_point_coordinates(0).back();
+ interval_temperature[i] = (maximum_temperature[i]-minimum_temperature[i])/(material_lookup[i]->get_number_of_coordinates(0)-1);
+ minimum_pressure[i] = material_lookup[i]->get_interpolation_point_coordinates(1).front();
+ maximum_pressure[i] = material_lookup[i]->get_interpolation_point_coordinates(1).back();
+ interval_pressure[i] = (maximum_pressure[i]-minimum_pressure[i])/(material_lookup[i]->get_number_of_coordinates(1)-1);
+ }
+ }
+
+
+ template
+ double
+ PhaseFunctionDiscrete::compute_value (const PhaseFunctionInputs &in) const
+ {
+ // the percentage of material that has undergone the transition
+ double function_value;
+
+ // lookup the most abundant phase
+ unsigned int start_phase_transition_index = 0;
+ unsigned int n_comp = 0;
+ for (unsigned int n_relevant_fields = 0 ; n_relevant_fields < this->introspection().n_chemical_composition_fields() + 1 ; n_relevant_fields ++)
+ {
+ if (in.phase_index < start_phase_transition_index + n_phase_transitions_per_chemical_composition[n_relevant_fields])
+ {
+ n_comp = n_relevant_fields ;
+ break;
+ }
+ start_phase_transition_index += n_phase_transitions_per_chemical_composition[n_relevant_fields];
+ }
+
+ const double pressure_in_bar = in.pressure/1e5;
+
+ Assert (in.temperature >= minimum_temperature[n_comp] && in.temperature < maximum_temperature[n_comp], ExcInternalError());
+ Assert (pressure_in_bar >= minimum_pressure[n_comp] && pressure_in_bar < maximum_pressure[n_comp], ExcInternalError());
+
+ const std::vector &temperature_points = material_lookup[n_comp]->get_interpolation_point_coordinates(0);
+ const std::vector &pressure_points = material_lookup[n_comp]->get_interpolation_point_coordinates(1);
+
+ // round T and p points to exact values in the table
+ const unsigned int i_T = static_cast(std::round((in.temperature - minimum_temperature[n_comp]) / interval_temperature[n_comp]));
+ const unsigned int i_p = static_cast(std::round((pressure_in_bar - minimum_pressure[n_comp]) / interval_pressure[n_comp]));
+
+ Point<2> temperature_pressure(temperature_points[i_T], pressure_points[i_p]);
+
+ // determine the dominant phase index
+ const unsigned int dominant_phase_index = static_cast(std::round(material_lookup[n_comp]->get_data(temperature_pressure, 7)));
+
+ // match the dominant phase index to one of the transition indicators
+ unsigned int matched_phase_transition_index = numbers::invalid_unsigned_int;
+ for (unsigned int i = start_phase_transition_index;
+ i < start_phase_transition_index + n_phase_transitions_per_chemical_composition[n_comp];
+ i++)
+ {
+ if (transition_indicators[i] == dominant_phase_index)
+ {
+ matched_phase_transition_index = i;
+ }
+ }
+
+ // determine the value of phase function to facilitate the exact transition
+ if ((matched_phase_transition_index != numbers::invalid_unsigned_int) && in.phase_index <= matched_phase_transition_index)
+ function_value = 1.0;
+ else
+ function_value = 0.0;
+
+ return function_value;
+ }
+
+
+ template
+ double
+ PhaseFunctionDiscrete::compute_derivative () const
+ {
+ // raises an error to ensure that a phase derivative request is not made for this phase function.
+ AssertThrow(false, ExcNotImplemented());
+ }
+
+
+ template
+ unsigned int
+ PhaseFunctionDiscrete::
+ n_phase_transitions () const
+ {
+ return transition_indicators.size();
+ }
+
+
+
+ template
+ unsigned int
+ PhaseFunctionDiscrete::
+ n_phases () const
+ {
+ return n_phases_total;
+ }
+
+
+
+ template
+ unsigned int
+ PhaseFunctionDiscrete::
+ n_phases_over_all_chemical_compositions () const
+ {
+ return n_phases_total_chemical_compositions;
+ }
+
+
+
+ template
+ const std::vector &
+ PhaseFunctionDiscrete::n_phase_transitions_for_each_composition () const
+ {
+ return *n_phase_transitions_per_composition;
+ }
+
+
+
+ template
+ const std::vector &
+ PhaseFunctionDiscrete::n_phases_for_each_composition () const
+ {
+ return n_phases_per_composition;
+ }
+
+
+
+ template
+ const std::vector &
+ PhaseFunctionDiscrete::n_phase_transitions_for_each_chemical_composition () const
+ {
+ return n_phase_transitions_per_chemical_composition;
+ }
+
+
+
+ template
+ const std::vector &
+ PhaseFunctionDiscrete::n_phases_for_each_chemical_composition () const
+ {
+ return n_phases_per_chemical_composition;
+ }
+
+
+
+ template
+ void
+ PhaseFunctionDiscrete::declare_parameters (ParameterHandler &prm)
+ {
+ prm.declare_entry ("Phase transition indicators", "",
+ Patterns::Anything(),
+ "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. ");
+
+ prm.declare_entry ("Data directory", "$ASPECT_SOURCE_DIR/data/material-model/entropy-table/pyrtable",
+ Patterns::DirectoryName (),
+ "The path to the model data. The path may also include the special "
+ "text '$ASPECT_SOURCE_DIR' which will be interpreted as the path "
+ "in which the ASPECT source files were located when ASPECT was "
+ "compiled. This interpretation allows, for example, to reference "
+ "files located in the `data/' subdirectory of ASPECT. ");
+
+ prm.declare_entry ("Material file names", "material_table_temperature_pressure_small.txt",
+ Patterns::List (Patterns::Anything()),
+ "The file names of the material data (material "
+ "data is assumed to be in order with the ordering "
+ "of the compositional fields). Note that there are "
+ "two options on how many files need to be listed "
+ "here: 1. If only one file is provided, it is used "
+ "for the whole model domain, and compositional fields "
+ "are ignored. 2. If there is one more file name than the "
+ "number of compositional fields, then the first file is "
+ "assumed to define a `background composition' that is "
+ "modified by the compositional fields. These data files need "
+ "to have the same structure as the one necessary for equation "
+ "of state plus a new column for the phase indexes, which amounts "
+ "to 8 columns in total.");
+ }
+
+
+
+ template
+ void
+ PhaseFunctionDiscrete::parse_parameters (ParameterHandler &prm)
+ {
+ // Retrieve the list of composition names
+ const std::vector list_of_composition_names = this->introspection().get_composition_names();
+
+ data_directory = Utilities::expand_ASPECT_SOURCE_DIR(prm.get ("Data directory"));
+ material_file_names = Utilities::split_string_list(prm.get ("Material file names"));
+
+ // Make options file for parsing maps to double arrays
+ std::vector chemical_field_names = this->introspection().chemical_composition_field_names();
+ chemical_field_names.insert(chemical_field_names.begin(),"background");
+
+ Utilities::MapParsing::Options options(chemical_field_names, "Phase transition indicators");
+ options.allow_missing_keys = true;
+ options.allow_multiple_values_per_key = true;
+ options.store_values_per_key = true;
+ options.n_values_per_key = std::vector();
+
+ const std::vector transition_indicators_double = Utilities::MapParsing::parse_map_to_double_array (prm.get("Phase transition indicators"),
+ options);
+
+ transition_indicators.reserve(transition_indicators_double.size());
+ for (const double transition_indicator: transition_indicators_double)
+ transition_indicators.push_back(static_cast(std::round(transition_indicator)));
+ n_phase_transitions_per_composition = std::make_unique>(options.n_values_per_key);
+
+ n_phases_total = 0;
+ n_phases_per_composition.clear();
+ for (unsigned int n : *n_phase_transitions_per_composition)
+ {
+ n_phases_per_composition.push_back(n+1);
+ n_phases_total += n+1;
+ }
+
+ // The background field is always the first composition
+ n_phases_per_chemical_composition = {n_phases_per_composition[0]};
+ n_phase_transitions_per_chemical_composition = {n_phases_per_composition[0] - 1};
+ n_phases_total_chemical_compositions = n_phases_per_composition[0];
+ for (auto i : this->introspection().chemical_composition_field_indices())
+ {
+ n_phases_per_chemical_composition.push_back(n_phases_per_composition[i+1]);
+ n_phase_transitions_per_chemical_composition.push_back(n_phases_per_composition[i+1] - 1);
+ n_phases_total_chemical_compositions += n_phases_per_composition[i+1];
+ }
+ }
+
+
template
double
PhaseFunction::compute_value (const PhaseFunctionInputs &in) const
@@ -1540,7 +1801,8 @@ namespace aspect
const unsigned int, \
MaterialModelOutputs &); \
template struct PhaseFunctionInputs; \
- template class PhaseFunction;
+ template class PhaseFunction; \
+ template class PhaseFunctionDiscrete;
ASPECT_INSTANTIATE(INSTANTIATE)
diff --git a/source/material_model/visco_plastic.cc b/source/material_model/visco_plastic.cc
index f96fb8bc4b4..d76bac94313 100644
--- a/source/material_model/visco_plastic.cc
+++ b/source/material_model/visco_plastic.cc
@@ -30,6 +30,17 @@ namespace aspect
{
namespace MaterialModel
{
+ template
+ void
+ ViscoPlastic::initialize()
+ {
+ if (use_dominant_phase_for_viscosity)
+ {
+ phase_function_discrete->initialize();
+ }
+ }
+
+
template
bool
ViscoPlastic::
@@ -102,6 +113,10 @@ namespace aspect
// While the number of phases is fixed, the value of the phase function is updated for every point
std::vector phase_function_values(phase_function.n_phase_transitions(), 0.0);
+ std::vector phase_function_discrete_values = (use_dominant_phase_for_viscosity?
+ std::vector(phase_function_discrete->n_phase_transitions(), 0.0): std::vector());
+
+
// Loop through all requested points
for (unsigned int i=0; i < in.n_evaluation_points(); ++i)
{
@@ -180,14 +195,31 @@ namespace aspect
// to compute the elastic force term.
bool plastic_yielding = false;
IsostrainViscosities isostrain_viscosities;
+
if (in.requests_property(MaterialProperties::viscosity) || in.requests_property(MaterialProperties::additional_outputs))
{
// Currently, the viscosities for each of the compositional fields are calculated assuming
// isostrain amongst all compositions, allowing calculation of the viscosity ratio.
// TODO: This is only consistent with viscosity averaging if the arithmetic averaging
// scheme is chosen. It would be useful to have a function to calculate isostress viscosities.
- isostrain_viscosities =
- rheology->calculate_isostrain_viscosities(in, i, volume_fractions, phase_function_values, n_phase_transitions_for_each_chemical_composition);
+ // Methods for averaging among phases depend on whether the most dominant phases are looked up
+ // for each composition in its respective lookup table. A separate phase function class
+ // manages the averaging if the user chooses this option.
+ if (use_dominant_phase_for_viscosity)
+ {
+ for (unsigned int j=0; j < phase_function_discrete->n_phase_transitions(); ++j)
+ {
+ phase_inputs.phase_index = j;
+ phase_function_discrete_values[j] = phase_function_discrete->compute_value(phase_inputs);
+ }
+ isostrain_viscosities =
+ rheology->calculate_isostrain_viscosities(in, i, volume_fractions, phase_function_discrete_values, phase_function_discrete->n_phase_transitions_for_each_chemical_composition());
+ }
+ else
+ {
+ isostrain_viscosities =
+ rheology->calculate_isostrain_viscosities(in, i, volume_fractions, phase_function_values, n_phase_transitions_for_each_chemical_composition);
+ }
// The isostrain condition implies that the viscosity averaging should be arithmetic (see above).
// We have given the user freedom to apply alternative bounds, because in diffusion-dominated
@@ -298,8 +330,16 @@ namespace aspect
{
prm.enter_subsection ("Visco Plastic");
{
+ prm.declare_entry ("Use dominant phase for viscosity","false",
+ Patterns::Bool (),
+ "Whether to look up the dominant phase for each composition in its respective "
+ "material data file to calculate viscosity. This allows each phase to have distinct "
+ "rheological parameterizations.");
+
MaterialUtilities::PhaseFunction::declare_parameters(prm);
+ MaterialUtilities::PhaseFunctionDiscrete::declare_parameters(prm);
+
EquationOfState::MulticomponentIncompressible::declare_parameters (prm);
Rheology::ViscoPlastic::declare_parameters(prm);
@@ -375,7 +415,20 @@ namespace aspect
rheology = std::make_unique>();
rheology->initialize_simulator (this->get_simulator());
- rheology->parse_parameters(prm, std::make_unique>(n_phases_for_each_chemical_composition));
+
+ use_dominant_phase_for_viscosity = prm.get_bool ("Use dominant phase for viscosity");
+ if (use_dominant_phase_for_viscosity)
+ {
+ phase_function_discrete = std::make_unique>();
+ phase_function_discrete->initialize_simulator (this->get_simulator());
+ phase_function_discrete->parse_parameters (prm);
+ rheology->parse_parameters(prm, std::make_unique>(phase_function_discrete->n_phases_for_each_chemical_composition()));
+ }
+ else
+ {
+ rheology->parse_parameters(prm, std::make_unique>(n_phases_for_each_chemical_composition));
+ }
+
}
prm.leave_subsection();
}
diff --git a/source/structured_data.cc b/source/structured_data.cc
index 62da14cd9b8..0caefa1d25a 100644
--- a/source/structured_data.cc
+++ b/source/structured_data.cc
@@ -145,6 +145,13 @@ namespace aspect
return maximum_component_value[component];
}
+ template
+ unsigned
+ StructuredDataLookup::get_number_of_coordinates(const unsigned int dimension) const
+ {
+ return table_points[dimension];
+ }
+
namespace
diff --git a/tests/visco_plastic_phases_discrete.prm b/tests/visco_plastic_phases_discrete.prm
new file mode 100644
index 00000000000..208e158a8af
--- /dev/null
+++ b/tests/visco_plastic_phases_discrete.prm
@@ -0,0 +1,113 @@
+# This prm file is used to generate a phase diagram.
+# and test the implementation of discrete phase functions
+# in the visco-plastic material model.
+# In the outputs, the material statistics are checked for
+# correct viscosity values.
+
+set Dimension = 2
+set CFL number = 1.0
+set End time = 0
+set Adiabatic surface temperature = 1673.0
+set Use years in output instead of seconds = true
+set Nonlinear solver scheme = no Advection, no Stokes
+
+# 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 = 800e3
+ 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
+
+subsection Gravity model
+ set Model name = vertical
+
+ subsection Vertical
+ set Magnitude = 10.0
+ end
+end
+
+# Temperature boundary and initial conditions
+# The initial temperature increases linearly from
+# the left to the right boundary.
+subsection Boundary temperature model
+ set Fixed temperature boundary indicators = left, right
+ set List of model names = box
+
+ subsection Box
+ set Left temperature = 273
+ set Right temperature = 2273
+ 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.0, Tr=2273
+ set Function expression = Tl * (x - XMAX)/(-XMAX) + Tr * x / XMAX
+ end
+end
+
+# Value for material model
+# A set of parameters are assigned to the conventional
+# phase function interface. A single phase with index 1
+# is added to the background composition with a smaller
+# viscosity by using the discrete phase transition.
+subsection Material model
+ set Model name = visco plastic
+
+ subsection Visco Plastic
+ set Reference temperature = 273
+ set Maximum viscosity = 1e24
+ set Phase transition depths = background:410e3|520e3|560e3|670e3|670e3|670e3|670e3
+ set Phase transition widths = background:5e3|5e3|5e3|5e3|5e3|5e3|5e3
+ set Phase transition temperatures = background:1662.0|1662.0|1662.0|1662.0|1662.0|1662.0|1662.0
+ set Phase transition Clapeyron slopes = background:4e6|4.1e6|4e6|-2e6|4e6|-3.1e6|1.3e6
+ set Densities = background: 3300.0|3394.4|3442.1|3453.2|3617.6|3691.5|3774.7|3929.1
+ set Heat capacities = 1250.0
+ set Thermal expansivities = 0.0
+ # The options to look up the dominant phase
+ set Use dominant phase for viscosity = true
+ set Data directory = $ASPECT_SOURCE_DIR/data/material-model/entropy-table/pyrtable/
+ set Material file names = material_table_temperature_pressure_small.txt
+ set Viscous flow law = diffusion
+ set Viscosity averaging scheme = harmonic
+ set Phase transition indicators = background:1
+ set Prefactors for diffusion creep = background:5e-21|5e-19
+ set Grain size exponents for diffusion creep = 0.0
+ set Activation energies for diffusion creep = 0.0
+ set Activation volumes for diffusion creep = 0.0
+ end
+end
+
+subsection Boundary velocity model
+ set Tangential velocity boundary indicators = top, bottom, left, right
+end
+
+subsection Postprocess
+ set List of postprocessors = visualization, material statistics
+
+ subsection Visualization
+ set Output format = gnuplot
+ set List of output variables = material properties
+ set Time between graphical output = 0e6
+
+ subsection Material properties
+ set List of material properties = density, thermal expansivity, specific heat, viscosity
+ end
+ end
+end
diff --git a/tests/visco_plastic_phases_discrete/screen-output b/tests/visco_plastic_phases_discrete/screen-output
new file mode 100644
index 00000000000..5a70aa37ce8
--- /dev/null
+++ b/tests/visco_plastic_phases_discrete/screen-output
@@ -0,0 +1,14 @@
+
+Number of active cells: 2,500 (on 1 levels)
+Number of degrees of freedom: 33,204 (20,402+2,601+10,201)
+
+*** Timestep 0: t=0 years, dt=0 years
+
+ Postprocessing:
+ Writing graphical output: output-visco_plastic_phases_discrete/solution/solution-00000
+ Average density / Average viscosity / Total mass: 3477 kg/m^3, 7.569e+19 Pa s, 2.226e+15 kg
+
+Termination requested by criterion: end time
+
+
+
diff --git a/tests/visco_plastic_phases_discrete/statistics b/tests/visco_plastic_phases_discrete/statistics
new file mode 100644
index 00000000000..329b18a90d6
--- /dev/null
+++ b/tests/visco_plastic_phases_discrete/statistics
@@ -0,0 +1,12 @@
+# 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 nonlinear iterations
+# 8: Visualization file name
+# 9: Average density (kg/m^3)
+# 10: Average viscosity (Pa s)
+# 11: Total mass (kg)
+0 0.000000000000e+00 0.000000000000e+00 2500 23003 10201 0 output-visco_plastic_phases_discrete/solution/solution-00000 3.47743073e+03 7.56900000e+19 2.22555567e+15
diff --git a/tests/visco_plastic_phases_discrete_comp.prm b/tests/visco_plastic_phases_discrete_comp.prm
new file mode 100644
index 00000000000..259de2573c9
--- /dev/null
+++ b/tests/visco_plastic_phases_discrete_comp.prm
@@ -0,0 +1,146 @@
+# This prm file is used to generate a phase diagram.
+# and test the implementation of discrete phase functions
+# in the visco-plastic material model. Specifically, this test
+# checks the implementation for multiple compositions.
+# In the outputs, the material statistics are checked for
+# correct viscosity values.
+
+set Dimension = 2
+set CFL number = 1.0
+set End time = 0
+set Adiabatic surface temperature = 1673.0
+set Use years in output instead of seconds = true
+set Nonlinear solver scheme = no Advection, no Stokes
+
+# 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 = 800e3
+ 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
+
+subsection Gravity model
+ set Model name = vertical
+
+ subsection Vertical
+ set Magnitude = 10.0
+ end
+end
+
+# Temperature boundary and initial conditions
+# The initial temperature increases linearly from
+# the left to the right boundary.
+subsection Boundary temperature model
+ set Fixed temperature boundary indicators = left, right
+ set List of model names = box
+
+ subsection Box
+ set Left temperature = 273
+ set Right temperature = 2273
+ 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.0, Tr=2273
+ set Function expression = Tl * (x - XMAX)/(-XMAX) + Tr * x / XMAX
+ end
+end
+
+# Comment the following 3 sections for steinberg model and pyrolitic lookup table
+# Fields of composition
+subsection Compositional fields
+ set Number of fields = 1
+ set Names of fields = spharz
+ set Compositional field methods = field
+end
+
+# Initial composition model
+subsection Initial composition model
+ set List of model names = function
+
+ subsection Function
+ set Coordinate system = cartesian
+ set Function expression = 1.0
+ end
+end
+
+# Boundary composition model
+subsection Boundary composition model
+ set List of model names = box
+ set Fixed composition boundary indicators = top, bottom, left, right
+
+ subsection Box
+ set Bottom composition = 1.0
+ set Left composition = 1.0
+ set Right composition = 1.0
+ set Top composition = 1.0
+ end
+end
+
+# Value for material model
+# A set of parameters are assigned to the conventional
+# phase function interface. A single phase with index 1
+# is added to the background composition with a smaller
+# viscosity by using the discrete phase transition.
+subsection Material model
+ set Model name = visco plastic
+
+ subsection Visco Plastic
+ set Reference temperature = 273
+ set Maximum viscosity = 1e25
+ set Minimum viscosity = 1e13
+ set Phase transition depths = background:410e3|520e3|560e3|670e3|670e3|670e3|670e3, spharz: 410e3|520e3|560e3|670e3|670e3|670e3|670e3
+ set Phase transition widths = background:5e3|5e3|5e3|5e3|5e3|5e3|5e3, spharz: 5e3|5e3|5e3|5e3|5e3|5e3|5e3
+ set Phase transition temperatures = background:1662.0|1662.0|1662.0|1662.0|1662.0|1662.0|1662.0, spharz: 1662.0|1662.0|1662.0|1662.0|1662.0|1662.0|1662.0
+ set Phase transition Clapeyron slopes = background:4e6|4.1e6|4e6|-2e6|4e6|-3.1e6|1.3e6, spharz: 4e6|4.1e6|4e6|-2e6|4e6|-3.1e6|1.3e6
+ set Densities = background: 3300.0|3394.4|3442.1|3453.2|3617.6|3691.5|3774.7|3929.1, spharz: 3235.0|3372.3|3441.7|3441.7|3680.8|3717.8|3759.4|3836.6
+ set Heat capacities = 1250.0
+ set Thermal expansivities = 0.0
+ # The options to look up the dominant phase
+ set Use dominant phase for viscosity = true
+ set Material file names = material_table_temperature_pressure_small.txt, material_table_temperature_pressure_small.txt
+ set Phase transition indicators = background:1, spharz:1
+ set Viscous flow law = diffusion
+ set Viscosity averaging scheme = harmonic
+ set Data directory = $ASPECT_SOURCE_DIR/data/material-model/entropy-table/pyrtable/
+ set Prefactors for diffusion creep = background:5e-21|5e-19, spharz:5e-25|5e-15
+ set Grain size exponents for diffusion creep = 0.0
+ set Activation energies for diffusion creep = 0.0
+ set Activation volumes for diffusion creep = 0.0
+ end
+end
+
+subsection Boundary velocity model
+ set Tangential velocity boundary indicators = top, bottom, left, right
+end
+
+subsection Postprocess
+ set List of postprocessors = visualization, material statistics
+
+ subsection Visualization
+ set Output format = gnuplot
+ set List of output variables = material properties
+ set Time between graphical output = 0e6
+
+ subsection Material properties
+ set List of material properties = density, thermal expansivity, specific heat, viscosity
+ end
+ end
+end
diff --git a/tests/visco_plastic_phases_discrete_comp/screen-output b/tests/visco_plastic_phases_discrete_comp/screen-output
new file mode 100644
index 00000000000..f39316c890f
--- /dev/null
+++ b/tests/visco_plastic_phases_discrete_comp/screen-output
@@ -0,0 +1,14 @@
+
+Number of active cells: 2,500 (on 1 levels)
+Number of degrees of freedom: 43,405 (20,402+2,601+10,201+10,201)
+
+*** Timestep 0: t=0 years, dt=0 years
+
+ Postprocessing:
+ Writing graphical output: output-visco_plastic_phases_discrete_comp/solution/solution-00000
+ Average density / Average viscosity / Total mass: 3428 kg/m^3, 7.6e+23 Pa s, 2.194e+15 kg
+
+Termination requested by criterion: end time
+
+
+
diff --git a/tests/visco_plastic_phases_discrete_comp/statistics b/tests/visco_plastic_phases_discrete_comp/statistics
new file mode 100644
index 00000000000..d5ad1dc97ae
--- /dev/null
+++ b/tests/visco_plastic_phases_discrete_comp/statistics
@@ -0,0 +1,13 @@
+# 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: Average density (kg/m^3)
+# 11: Average viscosity (Pa s)
+# 12: Total mass (kg)
+0 0.000000000000e+00 0.000000000000e+00 2500 23003 10201 10201 0 output-visco_plastic_phases_discrete_comp/solution/solution-00000 3.42821554e+03 7.60000000e+23 2.19405794e+15
diff --git a/tests/visco_plastic_phases_discrete_comp_invert.prm b/tests/visco_plastic_phases_discrete_comp_invert.prm
new file mode 100644
index 00000000000..8c4cb508d45
--- /dev/null
+++ b/tests/visco_plastic_phases_discrete_comp_invert.prm
@@ -0,0 +1,147 @@
+# This prm file is used to generate a phase diagram.
+# and test the implementation of discrete phase functions
+# in the visco-plastic material model. Specifically, this test
+# inverts the indicator for phase transition lookup from 0 to 1
+# and checks the implementation for multiple compositions.
+# In the outputs, the material statistics are checked for
+# correct viscosity values.
+
+set Dimension = 2
+set CFL number = 1.0
+set End time = 0
+set Adiabatic surface temperature = 1673.0
+set Use years in output instead of seconds = true
+set Nonlinear solver scheme = no Advection, no Stokes
+
+# 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 = 800e3
+ 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
+
+subsection Gravity model
+ set Model name = vertical
+
+ subsection Vertical
+ set Magnitude = 10.0
+ end
+end
+
+# Temperature boundary and initial conditions
+# The initial temperature increases linearly from
+# the left to the right boundary.
+subsection Boundary temperature model
+ set Fixed temperature boundary indicators = left, right
+ set List of model names = box
+
+ subsection Box
+ set Left temperature = 273
+ set Right temperature = 2273
+ 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.0, Tr=2273
+ set Function expression = Tl * (x - XMAX)/(-XMAX) + Tr * x / XMAX
+ end
+end
+
+# Comment the following 3 sections for steinberg model and pyrolitic lookup table
+# Fields of composition
+subsection Compositional fields
+ set Number of fields = 1
+ set Names of fields = spharz
+ set Compositional field methods = field
+end
+
+# Initial composition model
+subsection Initial composition model
+ set List of model names = function
+
+ subsection Function
+ set Coordinate system = cartesian
+ set Function expression = 1.0
+ end
+end
+
+# Boundary composition model
+subsection Boundary composition model
+ set List of model names = box
+ set Fixed composition boundary indicators = top, bottom, left, right
+
+ subsection Box
+ set Bottom composition = 1.0
+ set Left composition = 1.0
+ set Right composition = 1.0
+ set Top composition = 1.0
+ end
+end
+
+# Value for material model
+# A set of parameters are assigned to the conventional
+# phase function interface. A single phase with index 1
+# is added to the background composition with a smaller
+# viscosity by using the discrete phase transition.
+subsection Material model
+ set Model name = visco plastic
+
+ subsection Visco Plastic
+ set Reference temperature = 273
+ set Maximum viscosity = 1e25
+ set Minimum viscosity = 1e13
+ set Phase transition depths = background:410e3|520e3|560e3|670e3|670e3|670e3|670e3, spharz: 410e3|520e3|560e3|670e3|670e3|670e3|670e3
+ set Phase transition widths = background:5e3|5e3|5e3|5e3|5e3|5e3|5e3, spharz: 5e3|5e3|5e3|5e3|5e3|5e3|5e3
+ set Phase transition temperatures = background:1662.0|1662.0|1662.0|1662.0|1662.0|1662.0|1662.0, spharz: 1662.0|1662.0|1662.0|1662.0|1662.0|1662.0|1662.0
+ set Phase transition Clapeyron slopes = background:4e6|4.1e6|4e6|-2e6|4e6|-3.1e6|1.3e6, spharz: 4e6|4.1e6|4e6|-2e6|4e6|-3.1e6|1.3e6
+ set Densities = background: 3300.0|3394.4|3442.1|3453.2|3617.6|3691.5|3774.7|3929.1, spharz: 3235.0|3372.3|3441.7|3441.7|3680.8|3717.8|3759.4|3836.6
+ set Heat capacities = 1250.0
+ set Thermal expansivities = 0.0
+ # The options to look up the dominant phase
+ set Use dominant phase for viscosity = true
+ set Material file names = material_table_temperature_pressure_small.txt, material_table_temperature_pressure_small.txt
+ set Phase transition indicators = background:0, spharz:0
+ set Viscous flow law = diffusion
+ set Viscosity averaging scheme = harmonic
+ set Data directory = $ASPECT_SOURCE_DIR/data/material-model/entropy-table/pyrtable/
+ set Prefactors for diffusion creep = background:5e-21|5e-19, spharz:5e-25|5e-15
+ set Grain size exponents for diffusion creep = 0.0
+ set Activation energies for diffusion creep = 0.0
+ set Activation volumes for diffusion creep = 0.0
+ end
+end
+
+subsection Boundary velocity model
+ set Tangential velocity boundary indicators = top, bottom, left, right
+end
+
+subsection Postprocess
+ set List of postprocessors = visualization, material statistics
+
+ subsection Visualization
+ set Output format = gnuplot
+ set List of output variables = material properties
+ set Time between graphical output = 0e6
+
+ subsection Material properties
+ set List of material properties = density, thermal expansivity, specific heat, viscosity
+ end
+ end
+end
diff --git a/tests/visco_plastic_phases_discrete_comp_invert/screen-output b/tests/visco_plastic_phases_discrete_comp_invert/screen-output
new file mode 100644
index 00000000000..00a874a6275
--- /dev/null
+++ b/tests/visco_plastic_phases_discrete_comp_invert/screen-output
@@ -0,0 +1,14 @@
+
+Number of active cells: 2,500 (on 1 levels)
+Number of degrees of freedom: 43,405 (20,402+2,601+10,201+10,201)
+
+*** Timestep 0: t=0 years, dt=0 years
+
+ Postprocessing:
+ Writing graphical output: output-visco_plastic_phases_discrete_comp_invert/solution/solution-00000
+ Average density / Average viscosity / Total mass: 3428 kg/m^3, 2.4e+23 Pa s, 2.194e+15 kg
+
+Termination requested by criterion: end time
+
+
+
diff --git a/tests/visco_plastic_phases_discrete_comp_invert/statistics b/tests/visco_plastic_phases_discrete_comp_invert/statistics
new file mode 100644
index 00000000000..64778aeba4e
--- /dev/null
+++ b/tests/visco_plastic_phases_discrete_comp_invert/statistics
@@ -0,0 +1,13 @@
+# 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: Average density (kg/m^3)
+# 11: Average viscosity (Pa s)
+# 12: Total mass (kg)
+0 0.000000000000e+00 0.000000000000e+00 2500 23003 10201 10201 0 output-visco_plastic_phases_discrete_comp_invert/solution/solution-00000 3.42821554e+03 2.40000000e+23 2.19405794e+15
diff --git a/tests/visco_plastic_phases_discrete_mixed.prm b/tests/visco_plastic_phases_discrete_mixed.prm
new file mode 100644
index 00000000000..d1dc980bea6
--- /dev/null
+++ b/tests/visco_plastic_phases_discrete_mixed.prm
@@ -0,0 +1,147 @@
+# This prm file is used to generate a phase diagram.
+# and test the implementation of discrete phase functions
+# in the visco-plastic material model. Specifically, this test
+# checks the implementation for multiple compositions with
+# the two compositions mixed by half and half.
+# In the outputs, the material statistics are checked for
+# correct viscosity values.
+
+set Dimension = 2
+set CFL number = 1.0
+set End time = 0
+set Adiabatic surface temperature = 1673.0
+set Use years in output instead of seconds = true
+set Nonlinear solver scheme = no Advection, no Stokes
+
+# 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 = 800e3
+ 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
+
+subsection Gravity model
+ set Model name = vertical
+
+ subsection Vertical
+ set Magnitude = 10.0
+ end
+end
+
+# Temperature boundary and initial conditions
+# The initial temperature increases linearly from
+# the left to the right boundary.
+subsection Boundary temperature model
+ set Fixed temperature boundary indicators = left, right
+ set List of model names = box
+
+ subsection Box
+ set Left temperature = 273
+ set Right temperature = 2273
+ 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.0, Tr=2273
+ set Function expression = Tl * (x - XMAX)/(-XMAX) + Tr * x / XMAX
+ end
+end
+
+# Comment the following 3 sections for steinberg model and pyrolitic lookup table
+# Fields of composition
+subsection Compositional fields
+ set Number of fields = 1
+ set Names of fields = spharz
+ set Compositional field methods = field
+end
+
+# Initial composition model
+subsection Initial composition model
+ set List of model names = function
+
+ subsection Function
+ set Coordinate system = cartesian
+ set Function expression = 0.5
+ end
+end
+
+# Boundary composition model
+subsection Boundary composition model
+ set List of model names = box
+ set Fixed composition boundary indicators = top, bottom, left, right
+
+ subsection Box
+ set Bottom composition = 0.5
+ set Left composition = 0.5
+ set Right composition = 0.5
+ set Top composition = 0.5
+ end
+end
+
+# Value for material model
+# A set of parameters are assigned to the conventional
+# phase function interface. A single phase with index 1
+# is added to the background composition with a smaller
+# viscosity by using the discrete phase transition.
+subsection Material model
+ set Model name = visco plastic
+
+ subsection Visco Plastic
+ set Reference temperature = 273
+ set Maximum viscosity = 1e25
+ set Minimum viscosity = 1e13
+ set Phase transition depths = background:410e3|520e3|560e3|670e3|670e3|670e3|670e3, spharz: 410e3|520e3|560e3|670e3|670e3|670e3|670e3
+ set Phase transition widths = background:5e3|5e3|5e3|5e3|5e3|5e3|5e3, spharz: 5e3|5e3|5e3|5e3|5e3|5e3|5e3
+ set Phase transition temperatures = background:1662.0|1662.0|1662.0|1662.0|1662.0|1662.0|1662.0, spharz: 1662.0|1662.0|1662.0|1662.0|1662.0|1662.0|1662.0
+ set Phase transition Clapeyron slopes = background:4e6|4.1e6|4e6|-2e6|4e6|-3.1e6|1.3e6, spharz: 4e6|4.1e6|4e6|-2e6|4e6|-3.1e6|1.3e6
+ set Densities = background: 3300.0|3394.4|3442.1|3453.2|3617.6|3691.5|3774.7|3929.1, spharz: 3235.0|3372.3|3441.7|3441.7|3680.8|3717.8|3759.4|3836.6
+ set Heat capacities = 1250.0
+ set Thermal expansivities = 0.0
+ # The options to look up the dominant phase
+ set Use dominant phase for viscosity = true
+ set Material file names = material_table_temperature_pressure_small.txt, material_table_temperature_pressure_small.txt
+ set Phase transition indicators = background:1, spharz:1
+ set Viscous flow law = diffusion
+ set Viscosity averaging scheme = harmonic
+ set Data directory = $ASPECT_SOURCE_DIR/data/material-model/entropy-table/pyrtable/
+ set Prefactors for diffusion creep = background:5e-21|5e-19, spharz:5e-25|5e-15
+ set Grain size exponents for diffusion creep = 0.0
+ set Activation energies for diffusion creep = 0.0
+ set Activation volumes for diffusion creep = 0.0
+ end
+end
+
+subsection Boundary velocity model
+ set Tangential velocity boundary indicators = top, bottom, left, right
+end
+
+subsection Postprocess
+ set List of postprocessors = visualization, material statistics
+
+ subsection Visualization
+ set Output format = gnuplot
+ set List of output variables = material properties
+ set Time between graphical output = 0e6
+
+ subsection Material properties
+ set List of material properties = density, thermal expansivity, specific heat, viscosity
+ end
+ end
+end
diff --git a/tests/visco_plastic_phases_discrete_mixed/screen-output b/tests/visco_plastic_phases_discrete_mixed/screen-output
new file mode 100644
index 00000000000..8600e529e02
--- /dev/null
+++ b/tests/visco_plastic_phases_discrete_mixed/screen-output
@@ -0,0 +1,14 @@
+
+Number of active cells: 2,500 (on 1 levels)
+Number of degrees of freedom: 43,405 (20,402+2,601+10,201+10,201)
+
+*** Timestep 0: t=0 years, dt=0 years
+
+ Postprocessing:
+ Writing graphical output: output-visco_plastic_phases_discrete_mixed/solution/solution-00000
+ Average density / Average viscosity / Total mass: 3453 kg/m^3, 1.509e+20 Pa s, 2.21e+15 kg
+
+Termination requested by criterion: end time
+
+
+
diff --git a/tests/visco_plastic_phases_discrete_mixed/statistics b/tests/visco_plastic_phases_discrete_mixed/statistics
new file mode 100644
index 00000000000..99827427e07
--- /dev/null
+++ b/tests/visco_plastic_phases_discrete_mixed/statistics
@@ -0,0 +1,13 @@
+# 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: Average density (kg/m^3)
+# 11: Average viscosity (Pa s)
+# 12: Total mass (kg)
+0 0.000000000000e+00 0.000000000000e+00 2500 23003 10201 10201 0 output-visco_plastic_phases_discrete_mixed/solution/solution-00000 3.45280169e+03 1.50873851e+20 2.20979308e+15