Address comments from Rene 11/11/2024

Co-authored-by: Rene Gassmoeller <[email protected]>

.vscode/launch.json

@@ -2,13 +2,37 @@
     "version": "0.2.0",
     "configurations": [
         {
-            "name": "C++ launch",
+            "name": "(gdb) Launch",
             "type": "cppdbg",
             "request": "launch",
-            "program": "${workspaceFolder}/build_discrete_phase_funcion_2/aspect-debug",
-            "args": ["/mnt/lochy1/ASPECT_DATA/EntropySubduction/discrete_phase_function/phase_diagram_discrete/visualizing_phase_diagram.prm"],
-            "environment": [{"name": "config", "value": "Debug"}],
-            "cwd": "${workspaceFolder}"
+            // Resolved by CMake Tools:
+            "program": "${command:cmake.launchTargetPath}",
+	    // *** replace --test by the name of the .prm file you want to run: ***
+            "args": ["--test"],
+            "stopAtEntry": true,
+	    // *** append to the path here if you want to run inside some directory like /cookbooks/ ***
+            "cwd": "${workspaceFolder}",
+            "environment": [
+                {
+                    // add the directory where our target was built to the PATHs
+                    // it gets resolved by CMake Tools:
+                    "name": "PATH",
+                    "value": "${env:PATH}:${command:cmake.getLaunchTargetDirectory}"
+                },
+                {
+                    "name": "OTHER_VALUE",
+                    "value": "Something something"
+                }
+            ],
+            "externalConsole": false,
+            "MIMode": "gdb",
+            "setupCommands": [
+                {
+                    "description": "Enable pretty-printing for gdb",
+                    "text": "-enable-pretty-printing",
+                    "ignoreFailures": true
+                }
+            ]
         }
     ]
 }

data/material-model/entropy-table/pyrtable/material_table_temperature_pressure_small.txt

@@ -1,7 +1,7 @@
 # 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  maph,1
+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                

doc/modules/changes/20241107_lhy11009

@@ -1,4 +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.
+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.
 <br>
 (Haoyuan Li, 2024/11/07)

include/aspect/material_model/utilities.h

@@ -551,7 +551,7 @@ namespace aspect
 
       /**
        * A class that bundles functionality to look up the dominant phase in
-       * tables for each respective commposition and export the values
+       * 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
@@ -648,17 +648,44 @@ namespace aspect
 
         private:
           /**
-          * Information about the location of data files.
+          * Directory path where data files are stored.
           */
           std::string data_directory;
+
+          /**
+           * List of file names containing material data for each composition.
+          */
           std::vector<std::string> material_file_names;
+
+          /**
+          * Minimum temperature values for each composition in the P-T table.
+          */
           std::vector<double> minimum_temperature;
+
+          /**
+          * Maximum temperature values for each composition in the P-T table.
+          */
           std::vector<double> maximum_temperature;
+
+          /**
+          * Temperature intervals used for each composition in the P-T table.
+          */
           std::vector<double> interval_temperature;
+
+          /**
+          * Minimum pressure values for each composition in the P-T table.
+          */
           std::vector<double> minimum_pressure;
+
+          /**
+          * Maximum pressure values for each composition in the P-T table.
+          */
           std::vector<double> maximum_pressure;
-          std::vector<double> interval_pressure;
 
+          /**
+          * Pressure intervals used for each composition in the P-T table.
+          */
+          std::vector<double> interval_pressure;
 
           /**
            * List of pointers to objects that read and process data we get from

include/aspect/material_model/visco_plastic.h

@@ -183,7 +183,7 @@ namespace aspect
       public:
         /**
          * Initialization function. Loads the material data and sets up
-         * pointers if requried
+         * pointers if called again.
          */
         void
         initialize () override;

include/aspect/structured_data.h

@@ -250,10 +250,11 @@ namespace aspect
 
         /**
         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 get_table_points(const unsigned int dimension) const;
+        unsigned get_number_of_coordinates(const unsigned int dimension) const;
 
       private:
         /**

source/material_model/utilities.cc

@@ -1163,9 +1163,10 @@ namespace aspect
       PhaseFunctionDiscrete<dim>::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 "
+                     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<double>(material_file_names.size());
         maximum_temperature = std::vector<double>(material_file_names.size());
         interval_temperature = std::vector<double>(material_file_names.size());
@@ -1180,12 +1181,14 @@ namespace aspect
             material_lookup[i]->load_file(data_directory+material_file_names[i],
                                           this->get_mpi_communicator());
 
+            assert(material_lookup[i]->has_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_table_points(0)-1);
+            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_table_points(1)-1);
+            interval_pressure[i] = (maximum_pressure[i]-minimum_pressure[i])/(material_lookup[i]->get_number_of_coordinates(1)-1);
           }
       }
 
@@ -1200,47 +1203,47 @@ namespace aspect
         // lookup the most abundant phase
         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++)
+        for (unsigned 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[i])
+            if (in.phase_index < start_phase_transition_index + n_phase_transitions_per_chemical_composition[n_relevant_fields])
               {
-                n_comp = i;
+                n_comp = n_relevant_fields ;
                 break;
               }
-            start_phase_transition_index += n_phase_transitions_per_chemical_composition[i];
+            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 (in.pressure/1e5 >= minimum_pressure[n_comp] && in.pressure/1e5 < maximum_pressure[n_comp], ExcInternalError());
+        Assert (pressure_in_bar >= minimum_pressure[n_comp] && pressure_in_bar < maximum_pressure[n_comp], ExcInternalError());
 
         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[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]));
+        const unsigned i_p = static_cast<unsigned>(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 maph = unsigned(material_lookup[n_comp]->get_data(temperature_pressure, 7));
+        const unsigned dominant_phase_index = unsigned(material_lookup[n_comp]->get_data(temperature_pressure, 7));
 
-        // match the dominant phase index to one of the transtion indicators
-        bool matched_maph = false;
-        unsigned matched_phase_transition_index = 0;
+        // match the dominant phase index to one of the transition indicators
+        unsigned matched_phase_transition_index = numbers::invalid_unsigned_int;
         for (unsigned i = start_phase_transition_index;
              i < start_phase_transition_index + n_phase_transitions_per_chemical_composition[n_comp];
              i++)
           {
-            if (unsigned(transition_indicators[i]) == maph)
+            if (unsigned(transition_indicators[i]) == dominant_phase_index)
               {
-                matched_maph = true;
                 matched_phase_transition_index = i;
               }
           }
 
         // determine the value of phase function to facilitate the exact transition
-        if (matched_maph && in.phase_index <= matched_phase_transition_index)
+        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;
@@ -1249,16 +1252,15 @@ namespace aspect
       }
 
 
-
       template <int dim>
       double
       PhaseFunctionDiscrete<dim>::compute_derivative (const PhaseFunctionInputs<dim> &in) const
       {
-        return 0;
+        // make sure the phase derivative of this phase function is not asked for.
+        AssertThrow(false, ExcNotImplemented());
       }
 
 
-
       template <int dim>
       unsigned int
       PhaseFunctionDiscrete<dim>::
@@ -1352,12 +1354,7 @@ namespace aspect
                            "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. If there are "
-                           "exactly as many files as compositional fields, the fields are "
-                           "assumed to represent the fractions of different materials "
-                           "and the average property is computed as a sum of "
-                           "the value of the compositional field times the "
-                           "material property of that field.");
+                           "modified by the compositional fields.");
       }
 
 
@@ -1366,24 +1363,26 @@ namespace aspect
       void
       PhaseFunctionDiscrete<dim>::parse_parameters (ParameterHandler &prm)
       {
-        // Establish that a background field is required here
-        const bool has_background_field = true;
-
         // Retrieve the list of composition names
         const std::vector<std::string> 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"));
 
-        n_phase_transitions_per_composition = std::make_unique<std::vector<unsigned int>>();
+        // Make options file for parsing maps to double arrays
+        std::vector<std::string> chemical_field_names = this->introspection().chemical_composition_field_names();
+        chemical_field_names.insert(chemical_field_names.begin(),"background");
 
-        transition_indicators          = Utilities::parse_map_to_double_array (prm.get("Phase transition indicators"),
-                                                                               list_of_composition_names,
-                                                                               has_background_field,
-                                                                               "Phase transition indicators",
-                                                                               true,
-                                                                               n_phase_transitions_per_composition,
-                                                                               true);
+        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<unsigned int>();
+
+        transition_indicators          = Utilities::MapParsing::parse_map_to_double_array (prm.get("Phase transition indicators"),
+                                         options);
+
+        n_phase_transitions_per_composition = std::make_unique<std::vector<unsigned int>>(options.n_values_per_key);
 
         n_phases_total = 0;
         n_phases_per_composition.clear();
@@ -1393,7 +1392,6 @@ namespace aspect
             n_phases_total += n+1;
           }
 
-        Assert(has_background_field == true, ExcInternalError());
         // 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};

source/structured_data.cc

@@ -147,7 +147,7 @@ namespace aspect
 
     template <int dim>
     unsigned
-    StructuredDataLookup<dim>::get_table_points(const unsigned int dimension) const
+    StructuredDataLookup<dim>::get_number_of_coordinates(const unsigned int dimension) const
     {
       return table_points[dimension];
     }