EntropyReader Class

benchmarks/entropy_adiabat/plugins/entropy_model.cc

@@ -47,9 +47,9 @@ namespace aspect
                    ExcMessage("The 'entropy model' material model requires the existence of a compositional field "
                               "named 'entropy'. This field does not exist."));
 
-      material_lookup = std::make_unique<Utilities::StructuredDataLookup<2>>(7,1.0);
-      material_lookup->load_file(data_directory+material_file_name,
-                                 this->get_mpi_communicator());
+      entropy_reader = std::make_unique<MaterialUtilities::Lookup::EntropyReader>();
+      entropy_reader->initialize(this->get_mpi_communicator(), data_directory, material_file_name);
+
       lateral_viscosity_prefactor_lookup = std::make_unique<internal::LateralViscosityLookup>(data_directory+lateral_viscosity_file_name,
                                            this->get_mpi_communicator());
     }
@@ -125,20 +125,20 @@ namespace aspect
           // This is a requirement of the projected density approximation for
           // the Stokes equation and not related to the entropy formulation.
           // Also convert pressure from Pa to bar, bar is used in the table.
-          Point<2> entropy_pressure(in.composition[i][entropy_index],
-                                    this->get_adiabatic_conditions().pressure(in.position[i]) / 1.e5);
+          const double entropy = in.composition[i][entropy_index];
+          const double pressure = this->get_adiabatic_conditions().pressure(in.position[i]) / 1.e5;
 
-          out.densities[i]                      = material_lookup->get_data(entropy_pressure,1);
-          out.thermal_expansion_coefficients[i] = material_lookup->get_data(entropy_pressure,2);
-          out.specific_heat[i]                  = material_lookup->get_data(entropy_pressure,3);
+          out.densities[i] = entropy_reader->density(entropy,pressure);
+          out.thermal_expansion_coefficients[i] = entropy_reader->thermal_expansivity(entropy,pressure);
+          out.specific_heat[i] = entropy_reader->specific_heat(entropy,pressure);
 
-          const Tensor<1,2> density_gradient    = material_lookup->get_gradients(entropy_pressure,1);
-          const Tensor<1,2> pressure_unit_vector ({0.0,1.0});
-          out.compressibilities[i]              = (density_gradient * pressure_unit_vector) / out.densities[i];
+          const Tensor<1, 2> density_gradient = entropy_reader->density_gradient(entropy,pressure);
+          const Tensor<1, 2> pressure_unit_vector({0.0, 1.0});
+          out.compressibilities[i] = (density_gradient * pressure_unit_vector) / out.densities[i];
 
 
           // Thermal conductivity can be pressure temperature dependent
-          const double temperature_lookup =  material_lookup->get_data(entropy_pressure,0);
+          const double temperature_lookup =  entropy_reader->temperature(entropy,pressure);
           out.thermal_conductivities[i] = thermal_conductivity(temperature_lookup, in.pressure[i], in.position[i]);
 
           out.entropy_derivative_pressure[i]    = 0.;
@@ -216,8 +216,8 @@ namespace aspect
           // fill seismic velocities outputs if they exist
           if (SeismicAdditionalOutputs<dim> *seismic_out = out.template get_additional_output<SeismicAdditionalOutputs<dim>>())
             {
-              seismic_out->vp[i] = material_lookup->get_data(entropy_pressure,4);
-              seismic_out->vs[i] = material_lookup->get_data(entropy_pressure,5);
+              seismic_out->vp[i] = entropy_reader->seismic_vp(entropy, pressure);
+              seismic_out->vs[i] = entropy_reader->seismic_vs(entropy, pressure);
             }
         }
     }

benchmarks/entropy_adiabat/plugins/entropy_model.h

@@ -28,6 +28,7 @@
 #include <aspect/material_model/rheology/ascii_depth_profile.h>
 #include <aspect/material_model/rheology/drucker_prager.h>
 #include <aspect/material_model/steinberger.h>
+#include <aspect/material_model/utilities.h>
 
 namespace aspect
 {
@@ -158,9 +159,10 @@ namespace aspect
         std::string lateral_viscosity_file_name;
 
         /**
-         * Pointer to the StructuredDataLookup object that holds the material data.
+         * Pointer to the EntropyReader that reads in material data for
+         * given entropy and pressure.
          */
-        std::unique_ptr<Utilities::StructuredDataLookup<2>> material_lookup;
+        std::unique_ptr<MaterialUtilities::Lookup::EntropyReader> entropy_reader;
 
         /**
          * Pointer to an object that reads and processes data for the lateral

include/aspect/material_model/utilities.h

@@ -31,7 +31,14 @@
 namespace aspect
 {
   template <int dim> class SimulatorAccess;
+  namespace Utilities
+  {
+    using namespace dealii;
+    using namespace dealii::Utilities;
 
+    template <int dim>
+    class StructuredDataLookup;
+  }
   namespace MaterialModel
   {
     using namespace dealii;
@@ -259,6 +266,78 @@ namespace aspect
                           const bool interpol,
                           const MPI_Comm &comm);
         };
+
+        /**
+        * This class reads in an entropy-pressure material table and looks up material
+        * properties for the given entropy and pressure.
+        */
+        class EntropyReader
+        {
+          public:
+
+            /**
+             * Read the table.
+             */
+            void
+            initialize(const MPI_Comm comm,
+                       const std::string data_directory,
+                       const std::string material_file_name);
+
+            /**
+             * Returns the specific heat for a given entropy and pressure.
+             */
+            double
+            specific_heat(const double entropy,
+                          const double pressure) const;
+
+            /**
+             * Returns the density for a given entropy and pressure.
+             */
+            double
+            density(const double entropy,
+                    const double pressure) const;
+
+            /**
+             * Returns the thermal_expansivity for a given entropy and pressure.
+             */
+            double
+            thermal_expansivity(const double entropy,
+                                const double pressure) const;
+
+            /**
+             * Returns the temperature for a given entropy and pressure.
+             */
+            double
+            temperature(const double entropy,
+                        const double pressure) const;
+
+            /**
+             * Returns the seismic p wave velocity for a given entropy and pressure.
+             */
+            double
+            seismic_vp(const double entropy,
+                       const double pressure) const;
+
+            /**
+             * Returns the seismic s wave velocity for a given entropy and pressure.
+             */
+            double
+            seismic_vs(const double entropy,
+                       const double pressure) const;
+
+            /**
+             * Returns density gradient for a given entropy and pressure.
+             */
+            Tensor<1, 2>
+            density_gradient(const double entropy,
+                             const double pressure) const;
+
+          private:
+            /**
+             * The StucturedDataLookup object that stores the material data.
+             */
+            std::unique_ptr<Utilities::StructuredDataLookup<2>> material_lookup;
+        };
       }
 
       /**

source/material_model/utilities.cc

@@ -21,7 +21,9 @@
 
 #include <aspect/global.h>
 #include <aspect/simulator_access.h>
+#include <aspect/structured_data.h>
 
+#include <aspect/geometry_model/interface.h>
 #include <aspect/adiabatic_conditions/interface.h>
 #include <aspect/gravity_model/interface.h>
 
@@ -786,6 +788,90 @@ namespace aspect
           AssertThrow(i == n_temperature*n_pressure, ExcMessage("Material table size not consistent with header."));
 
         }
+
+
+
+        void
+        EntropyReader::initialize(const MPI_Comm comm,
+                                  const std::string data_directory,
+                                  const std::string material_file_name)
+        {
+          material_lookup = std::make_unique<Utilities::StructuredDataLookup<2>>(7,1.0);
+          material_lookup->load_file(data_directory+material_file_name,
+                                     comm);
+        }
+
+
+
+        double
+        EntropyReader::specific_heat(const double entropy,
+                                     const double pressure) const
+        {
+          const double specific_heat = material_lookup->get_data({entropy,pressure}, 3);
+          return specific_heat;
+        }
+
+
+
+        double
+        EntropyReader::density(const double entropy,
+                               const double pressure) const
+        {
+          const double density = material_lookup->get_data({entropy,pressure}, 1);
+          return density;
+        }
+
+
+
+        double
+        EntropyReader::thermal_expansivity(const double entropy,
+                                           const double pressure) const
+        {
+          const double thermal_expansivity = material_lookup->get_data({entropy,pressure}, 2);
+          return thermal_expansivity;
+        }
+
+
+
+        double
+        EntropyReader::temperature(const double entropy,
+                                   const double pressure) const
+        {
+          const double temperature = material_lookup->get_data({entropy,pressure}, 0);
+          return temperature;
+        }
+
+
+
+        double
+        EntropyReader::seismic_vp(const double entropy,
+                                  const double pressure) const
+        {
+          const double seismic_vp = material_lookup->get_data({entropy,pressure}, 4);
+          return seismic_vp;
+        }
+
+
+
+        double
+        EntropyReader::seismic_vs(const double entropy,
+                                  const double pressure) const
+        {
+          const double seismic_vs = material_lookup->get_data({entropy,pressure}, 5);
+          return seismic_vs;
+        }
+
+
+
+        Tensor<1, 2>
+        EntropyReader::density_gradient(const double entropy,
+                                        const double pressure) const
+        {
+          const Tensor<1, 2> density_gradient= material_lookup->get_gradients({entropy,pressure}, 1);
+          return density_gradient;
+        }
+
+
       }