follow

geodynamics · Nov 5, 2024 · dad0466 · dad0466
1 parent 4704f45
commit dad0466
Show file tree

Hide file tree

Showing 4 changed files with 42 additions and 5 deletions.
diff --git a/include/aspect/material_model/utilities.h b/include/aspect/material_model/utilities.h
@@ -653,6 +653,13 @@ namespace aspect
           */
           std::string data_directory;
           std::vector<std::string> material_file_names;
+          std::vector<double> minimum_temperature;
+          std::vector<double> maximum_temperature;
+          std::vector<double> interval_temperature;
+          std::vector<double> minimum_pressure;
+          std::vector<double> maximum_pressure;
+          std::vector<double> interval_pressure;
+
 
           /**
            * List of pointers to objects that read and process data we get from

diff --git a/include/aspect/structured_data.h b/include/aspect/structured_data.h
@@ -248,6 +248,9 @@ namespace aspect
          */
         double get_maximum_component_value(const unsigned int component) const;
 
+        //todo_PT
+        unsigned get_table_points(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
@@ -1167,15 +1167,27 @@ namespace aspect
                      ExcMessage(" The 'discrete phase function' requires that the number of compositional fields "
                                 "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());
+        minimum_pressure = std::vector<double>(material_file_names.size());
+        maximum_pressure = std::vector<double>(material_file_names.size());
+        interval_pressure = std::vector<double>(material_file_names.size());
         for (unsigned i = 0; i < material_file_names.size(); ++i)
           {
             material_lookup
             .push_back(std::make_unique<Utilities::StructuredDataLookup<2>>(8,1.0));
 
             material_lookup[i]->load_file(data_directory+material_file_names[i],
                                           this->get_mpi_communicator());
-          }
 
+            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);
+            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);
+          }
       }
 
 
@@ -1189,12 +1201,19 @@ namespace aspect
 
         // lookup the most abundant phase
 
-        // Convert pressure from Pa to bar, bar is used in the table.
-        Point<2> temperature_pressure(in.temperature, in.pressure / 1.e5);
+        const unsigned i = 0;
+
+        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());
+
+        // const unsigned maph = unsigned(material_lookup[0]->get_data(temperature_pressure, 7));
+        const std::vector<double>& maph_points = material_lookup[0]->get_interpolation_point_coordinates(7);
 
-        const unsigned maph = unsigned(material_lookup[0]->get_data(temperature_pressure, 7));
+        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_point = i_p * material_lookup[0]->get_table_points(0) + i_T;
 
-        // Convert pressure from Pa to bar, bar is used in the table.
+        const unsigned maph = unsigned(maph_points[i_point]);
 
         // determine the value
         if (in.phase_index < maph)

diff --git a/source/structured_data.cc b/source/structured_data.cc
@@ -144,6 +144,14 @@ namespace aspect
     {
       return maximum_component_value[component];
     }
+
+    //todo_PT
+    template <int dim>
+    unsigned 
+    StructuredDataLookup<dim>::get_table_points(const unsigned int dimension) const
+    {
+      return table_points[dimension];
+    }