Entropy initial temperature lookup

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

@@ -0,0 +1,3 @@
+# H/nR [K], depth [km]
+0 0.0000000e+00
+0 3.0000000e+06

doc/modules/changes/20230913_haoyuan

@@ -0,0 +1,7 @@
+ New: There is now an entropy table lookup model
+ in the initial composition model.
+ This model takes an initial temperature field and
+ converts it to an initial entropy field using
+ a lookup table.
+ <br>
+ (Haoyuan Li, 2023/08/11)

include/aspect/initial_composition/entropy_table_lookup.h

@@ -0,0 +1,111 @@
+/*
+  Copyright (C) 2017 - 2023 by the authors of the ASPECT code.
+
+  This file is part of ASPECT.
+
+  ASPECT is free software; you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation; either version 2, or (at your option)
+  any later version.
+
+  ASPECT is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with ASPECT; see the file LICENSE.  If not see
+  <http://www.gnu.org/licenses/>.
+*/
+
+
+#ifndef _aspect_initial_composition_entropy_table_lookup_h
+#define _aspect_initial_composition_entropy_table_lookup_h
+
+#include <aspect/initial_composition/interface.h>
+#include <aspect/initial_temperature/interface.h>
+#include <aspect/simulator_access.h>
+#include <aspect/utilities.h>
+
+namespace aspect
+{
+  namespace InitialComposition
+  {
+    using namespace dealii;
+
+    /**
+     * A class that implements initial conditions for the entropy field
+     * Note that this plugin only
+     * works if there is a compositional field called 'entropy'.
+     * All compositional fields except entropy are not changed by this plugin.
+     *
+     * @ingroup InitialCompositionModels
+     */
+    template <int dim>
+    class EntropyTableLookUp : public Interface<dim>,
+      public SimulatorAccess<dim>
+    {
+      public:
+        /**
+         * Initialize the plugin.
+         */
+        void initialize () override;
+
+        /**
+         * Return the initial composition as a function of position and number
+         * of compositional field.
+         */
+        double initial_composition (const Point<dim> &position,
+                                    const unsigned int compositional_index) const override;
+
+        /**
+         * Declare the parameters this class takes through input files.
+         */
+        static
+        void
+        declare_parameters (ParameterHandler &prm);
+
+        /**
+         * Read the parameters this class declares from the parameter file.
+         */
+        void
+        parse_parameters (ParameterHandler &prm) override;
+
+      private:
+        /**
+         * Information about the location of data files.
+         */
+        std::string data_directory;
+        std::string material_file_name;
+
+        /**
+         * Index of the entropy in the compositional fields
+         */
+        unsigned entropy_index;
+
+        /**
+         * A shared pointer to the initial temperature object
+         * that ensures that the current object can continue
+         * to access the initial temperature object beyond the
+         * first time step.
+         */
+        std::shared_ptr<const aspect::InitialTemperature::Manager<dim>> initial_temperature_manager;
+
+        /**
+         * A shared pointer to the initial composition object
+         * that ensures that the current object can continue
+         * to access the initial composition object beyond the
+         * first time step.
+         */
+        std::shared_ptr<const aspect::InitialComposition::Manager<dim>> initial_composition_manager;
+
+        /**
+         * Pointer to the StructuredDataLookup object that holds the material data.
+         */
+        std::unique_ptr<Utilities::StructuredDataLookup<2>> material_lookup;
+    };
+  }
+}
+
+
+#endif

source/initial_composition/entropy_table_lookup.cc

@@ -0,0 +1,138 @@
+/*
+  Copyright (C) 2017 - 2023 by the authors of the ASPECT code.
+
+  This file is part of ASPECT.
+
+  ASPECT is free software; you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation; either version 2, or (at your option)
+  any later version.
+
+  ASPECT is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with ASPECT; see the file LICENSE.  If not see
+  <http://www.gnu.org/licenses/>.
+*/
+
+
+#include <aspect/initial_composition/entropy_table_lookup.h>
+#include <aspect/initial_temperature/interface.h>
+#include <aspect/adiabatic_conditions/interface.h>
+
+
+namespace aspect
+{
+  namespace InitialComposition
+  {
+    template <int dim>
+    void
+    EntropyTableLookUp<dim>::initialize()
+    {
+      AssertThrow (this->introspection().compositional_name_exists("entropy"),
+                   ExcMessage("The 'entropy model' material model requires the existence of a compositional field "
+                              "named 'entropy'. This field does not exist."));
+
+      // Make sure we keep track of the initial temperature manager and
+      // that it continues to live beyond the time when the simulator
+      // class releases its pointer to it.
+      initial_temperature_manager = this->get_initial_temperature_manager_pointer();
+
+      // Make sure we keep track of the initial composition manager and
+      // that it continues to live beyond the time when the simulator
+      // class releases its pointer to it.
+      initial_composition_manager = this->get_initial_composition_manager_pointer();
+
+      entropy_index = this->introspection().compositional_index_for_name("entropy");
+
+      material_lookup = std::make_unique<Utilities::StructuredDataLookup<2>>(7,1.0);
+      material_lookup->load_file(data_directory+material_file_name,
+                                 this->get_mpi_communicator());
+    }
+
+
+    template <int dim>
+    double
+    EntropyTableLookUp<dim>::
+    initial_composition (const Point<dim> &position,
+                         const unsigned int compositional_index) const
+    {
+      if (compositional_index == entropy_index)
+        {
+          const double temperature = initial_temperature_manager->initial_temperature(position);
+          const double pressure = this->get_adiabatic_conditions().pressure(position);
+
+          // Convert pressure from Pa to bar, bar is used in the table.
+          Point<2> temperature_pressure(temperature, pressure / 1.e5);
+
+          const double entropy = material_lookup->get_data(temperature_pressure, 0);
+
+          return entropy;
+        }
+      return 0.0;
+    }
+
+    template <int dim>
+    void
+    EntropyTableLookUp<dim>::declare_parameters (ParameterHandler &prm)
+    {
+      prm.enter_subsection("Initial composition model");
+      {
+        prm.enter_subsection("Entropy table lookup");
+        {
+          prm.declare_entry ("Data directory", "$ASPECT_SOURCE_DIR/data/material-model/entropy-table/pyrtable/",
+                             Patterns::DirectoryName (),
+                             "The path to the model data. The path may also include the special "
+                             "text '$ASPECT_SOURCE_DIR' which will be interpreted as the path "
+                             "in which the ASPECT source files were located when ASPECT was "
+                             "compiled. This interpretation allows, for example, to reference "
+                             "files located in the `data/' subdirectory of ASPECT.");
+          prm.declare_entry ("Material file name", "material_table_temperature_pressure.txt",
+                             Patterns::List (Patterns::Anything()),
+                             "The file name of the material data.");
+        }
+        prm.leave_subsection();
+      }
+      prm.leave_subsection();
+    }
+
+
+    template <int dim>
+    void
+    EntropyTableLookUp<dim>::parse_parameters (ParameterHandler &prm)
+    {
+      prm.enter_subsection("Initial composition model");
+      {
+        prm.enter_subsection("Entropy table lookup");
+        {
+          data_directory              = Utilities::expand_ASPECT_SOURCE_DIR(prm.get("Data directory"));
+          material_file_name          = prm.get("Material file name");
+        }
+        prm.leave_subsection();
+      }
+      prm.leave_subsection();
+    }
+
+  }
+}
+
+// explicit instantiations
+namespace aspect
+{
+  namespace InitialComposition
+  {
+    ASPECT_REGISTER_INITIAL_COMPOSITION_MODEL(EntropyTableLookUp,
+                                              "entropy table lookup",
+                                              "A class that implements initial conditions for the entropy field "
+                                              "by converting the initial temperature field through a look up table"
+                                              "Note that this plugin only works if there is a compositional field"
+                                              "called `entropy', and an additional look up table that can convert"
+                                              "pressure and temperature to entropy for each composition in the model."
+                                              "For all compositional fields except entropy this plugin returns 0.0, "
+                                              "and they are therefore not changed as long as the default `add' "
+                                              "operator is selected for this plugin.")
+  }
+}

tests/entropy_initial_lookup.cc

@@ -0,0 +1,21 @@
+/*
+  Copyright (C) 2022 by the authors of the ASPECT code.
+
+  This file is part of ASPECT.
+
+  ASPECT is free software; you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation; either version 2, or (at your option)
+  any later version.
+
+  ASPECT is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with ASPECT; see the file LICENSE.  If not see
+  <http://www.gnu.org/licenses/>.
+*/
+
+#include "../benchmarks/entropy_adiabat/plugins/entropy_advection.cc"