follow

geodynamics · Nov 6, 2024 · d07d104 · d07d104
1 parent 367a1ad
commit d07d104
Show file tree

Hide file tree

Showing 4 changed files with 34 additions and 5,711 deletions.
diff --git a/source/material_model/utilities.cc b/source/material_model/utilities.cc
@@ -1200,22 +1200,31 @@ namespace aspect
         double function_value;
 
         // lookup the most abundant phase
-        const unsigned i = 0;
+        unsigned n_phases_composition = 0;
+        unsigned n_comp = 0;
+        for (auto i : this->introspection().chemical_composition_field_indices())
+          {
+            n_phases_composition += n_phases_per_chemical_composition[i];
+            if (in.phase_index < n_phases_composition){
+              n_comp = i;
+              break;
+            }
+          }
 
-        Assert (in.temperature >= minimum_temperature[i] && in.temperature < maximum_temperature[i], ExcInternalError());
-        Assert (in.pressure/1e5 >= minimum_pressure[i] && in.pressure/1e5 < maximum_pressure[i], ExcInternalError());
+        Assert (in.temperature >= minimum_temperature[n_comp] && in.temperature < maximum_temperature[n_comp], ExcInternalError());
+        Assert (in.pressure/1e5 >= minimum_pressure[n_comp] && in.pressure/1e5 < maximum_pressure[n_comp], ExcInternalError());
 
-        const std::vector<double>& temperature_points = material_lookup[i]->get_interpolation_point_coordinates(0);
-        const std::vector<double>& pressure_points = material_lookup[i]->get_interpolation_point_coordinates(1);
+        const std::vector<double>& temperature_points = material_lookup[n_comp]->get_interpolation_point_coordinates(0);
+        const std::vector<double>& pressure_points = material_lookup[n_comp]->get_interpolation_point_coordinates(1);
 
         // round T and p points to exact values in the table
-        const unsigned i_T = static_cast<unsigned>(std::round((in.temperature - minimum_temperature[i]) / interval_temperature[i]));
-        const unsigned i_p = static_cast<unsigned>(std::round((in.pressure/1e5 - minimum_pressure[i]) / interval_pressure[i]));
+        const unsigned i_T = static_cast<unsigned>(std::round((in.temperature - minimum_temperature[n_comp]) / interval_temperature[n_comp]));
+        const unsigned i_p = static_cast<unsigned>(std::round((in.pressure/1e5 - 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 maph = unsigned(material_lookup[i]->get_data(temperature_pressure, 7));
+        const unsigned maph = unsigned(material_lookup[n_comp]->get_data(temperature_pressure, 7));
 
         // determine the value of phase function
         if (in.phase_index < maph)
@@ -1232,7 +1241,7 @@ namespace aspect
       double
       PhaseFunctionDiscrete<dim>::compute_derivative (const PhaseFunctionInputs<dim> &in) const
       {
-        return 0.0;
+        return 0;
       }
 
 

diff --git a/tests/visco_plastic_phases_discrete.prm b/tests/visco_plastic_phases_discrete.prm
@@ -114,7 +114,7 @@ subsection Boundary velocity model
 end
 
 subsection Postprocess
-  set List of postprocessors = visualization
+  set List of postprocessors = visualization, material statistics
 
   subsection Visualization
     set Output format                 = gnuplot