diff --git a/source/material_model/utilities.cc b/source/material_model/utilities.cc index 3009ff57eee..6d07c8c3049 100644 --- a/source/material_model/utilities.cc +++ b/source/material_model/utilities.cc @@ -1163,7 +1163,7 @@ namespace aspect void PhaseFunctionDiscrete::initialize() { - AssertThrow (this->introspection().get_number_of_fields_of_type(CompositionalFieldDescription::chemical_composition) == material_file_names.size(), + 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 " "matches the number of lookup tables.")); @@ -1200,15 +1200,15 @@ namespace aspect double function_value; // lookup the most abundant phase - unsigned n_phases_composition = 0; - unsigned n_comp = 0; - for (auto i : this->introspection().chemical_composition_field_indices()) + unsigned int start_phase_transition_index = 0; + unsigned int n_comp = 0; + for (unsigned i = 0 ; i < this->introspection().n_chemical_composition_fields() + 1 ; i++) { - n_phases_composition += n_phases_per_chemical_composition[i]; - if (in.phase_index < n_phases_composition){ + if (in.phase_index < start_phase_transition_index + n_phase_transitions_per_chemical_composition[i]){ n_comp = i; break; } + start_phase_transition_index += n_phase_transitions_per_chemical_composition[i]; } Assert (in.temperature >= minimum_temperature[n_comp] && in.temperature < maximum_temperature[n_comp], ExcInternalError()); @@ -1227,7 +1227,7 @@ namespace aspect const unsigned maph = unsigned(material_lookup[n_comp]->get_data(temperature_pressure, 7)); // determine the value of phase function - if (in.phase_index < maph) + if (in.phase_index - start_phase_transition_index < maph) function_value = 1.0; else function_value = 0.0; diff --git a/tests/visco_plastic_phases_discrete.prm b/tests/visco_plastic_phases_discrete.prm index 21709d95f03..4abf1adcfd3 100644 --- a/tests/visco_plastic_phases_discrete.prm +++ b/tests/visco_plastic_phases_discrete.prm @@ -1,6 +1,8 @@ # 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 +# correst viscosity values. set Dimension = 2 set CFL number = 1.0 @@ -95,12 +97,12 @@ subsection Material model 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 Viscous flow law = diffusion - set Viscosity averaging scheme = harmonic - # for discrete phase function 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 phases = background:1 set Prefactors for diffusion creep = background:5e-21|5e-19, spharz:5e-21 set Grain size exponents for diffusion creep = 0.0 diff --git a/tests/visco_plastic_phases_discrete_comp.prm b/tests/visco_plastic_phases_discrete_comp.prm new file mode 100644 index 00000000000..4c574c44b8b --- /dev/null +++ b/tests/visco_plastic_phases_discrete_comp.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. +# In the outputs, the material statistics are checked for +# correst 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 phases = 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_mixed.prm b/tests/visco_plastic_phases_discrete_mixed.prm new file mode 100644 index 00000000000..71e347bac15 --- /dev/null +++ b/tests/visco_plastic_phases_discrete_mixed.prm @@ -0,0 +1,148 @@ +# 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 +# correst 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 phases = 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