Merge pull request #6888 from lhy11009/prescribed_solution_initial_temperature

Prescribed solution initial temperature

Author: gassmoeller

doc/modules/changes/20260307_lhy11009

@@ -0,0 +1,3 @@
+Added: a new instance of 'initial temperature' in the prescribed solution plugin system.
+<br>
+(Haoyuan Li, 2025/12/30)

doc/sphinx/parameters/Material_20model.md

@@ -200,7 +200,7 @@ Viscous stress may also be limited by a non-linear stress limiter that has a for
 
  The value for the components of this formula and additional parameters are read from the parameter file in subsection  ’Material model/Visco Plastic’.
 
-‘viscoelastic’: An implementation of a simple linear viscoelastic rheology that only includes the deviatoric components of elasticity. Specifically, the viscoelastic rheology only takes into account the elastic shear strength (e.g., shear modulus), while the tensile and volumetric strength (e.g., Young’s and bulk modulus) are not considered. The model is incompressible and allows specifying an arbitrary number of compositional fields, where each field represents a different rock type or component of the viscoelastic stress tensor. The stress tensor in 2d and 3d, respectively, contains 3 or 6 components. The compositional fields representing these components must be named and listed in a very specific format, which is designed to minimize mislabeling stress tensor components as distinct ’compositional rock types’ (or vice versa). For 2d models, the first six compositional fields must be labeled ’stress\_xx’, ’stress\_yy’ and ’stress\_xy’, ’stress\_xx\_old’, ’stress\_yy\_old’ and ’stress\_xy\_old’, In 3d, the first twelve compositional fields must be labeled ’stress\_xx’, ’stress\_yy’, ’stress\_zz’, ’stress\_xy’, ’stress\_xz’, ’stress\_yz’, ’stress\_xx\_old’, ’stress\_yy\_old’, ’stress\_zz\_old’, ’stress\_xy\_old’, ’stress\_xz\_old’, ’stress\_yz\_old’.
+‘viscoelastic’: An implementation of a simple linear viscoelastic rheology that only includes the deviatoric components of elasticity. Specifically, the viscoelastic rheology only takes into account the elastic shear strength (e.g., shear modulus), while the tensile and volumetric strength (e.g., Young’s and bulk modulus) are not considered. The model is incompressible and allows specifying an arbitrary number of compositional fields, where each field represents a different rock type or component of the viscoelastic stress tensor. The stress tensor in 2d and 3d, respectively, contains 3 or 6 components. The compositional fields representing these components must be named and listed in a very specific format, which is designed to minimize mislabeling stress tensor components as distinct ’compositional rock types’ (or vice versa). For 2d models, the first six compositional fields of type stress must be labeled ’ve\_stress\_xx’, ’ve\_stress\_yy’ and ’ve\_stress\_xy’, ’ve\_stress\_xx\_old’, ’ve\_stress\_yy\_old’ and ’ve\_stress\_xy\_old’, In 3d, the first twelve compositional fields of type stress must be labeled ’ve\_stress\_xx’, ’ve\_stress\_yy’, ’ve\_stress\_zz’, ’ve\_stress\_xy’, ’ve\_stress\_xz’, ’ve\_stress\_yz’, ’ve\_stress\_xx\_old’, ’ve\_stress\_yy\_old’, ’ve\_stress\_zz\_old’,  ’ve\_stress\_xy\_old’, ’ve\_stress\_xz\_old’, ’ve\_stress\_yz\_old’.
 
  Expanding the model to include non-linear viscous flow (e.g., diffusion/dislocation creep) and plasticity would produce a constitutive relationship commonly referred to as partial elastoviscoplastic (e.g., pEVP) in the geodynamics community. While extensively discussed and applied within the geodynamics literature, notable references include: Moresi et al. (2003), J. Comp. Phys., v. 184, p. 476-497. Gerya and Yuen (2007), Phys. Earth. Planet. Inter., v. 163, p. 83-105. Gerya (2010), Introduction to Numerical Geodynamic Modeling. Kaus (2010), Tectonophysics, v. 484, p. 36-47. Choi et al. (2013), J. Geophys. Res., v. 118, p. 2429-2444. Keller et al. (2013), Geophys. J. Int., v. 195, p. 1406-1442.
 

doc/sphinx/parameters/Prescribed_20solution.md

@@ -9,17 +9,63 @@
 :name: parameters:Prescribed_20solution/List_20of_20model_20names
 **Default value:**
 
-**Pattern:** [MultipleSelection temperature function|velocity function ]
+**Pattern:** [MultipleSelection initial temperature|temperature function|velocity function ]
 
 **Documentation:** A comma-separated list of prescribed solution models that will be used to compute the solution in certain regions. These plugins are loaded in the order given, and are combined via the operators listed in ’List of model operators’.
 
 The following prescribed solution models are available:
 
+‘initial temperature’: Prescribe the temperature in a selected region using the active initial temperature model. The selected region is defined through an indicator function. At locations where the indicator value is greater than 0.5, the temperature is constrained to the initial temperature evaluated at that position.
+
 ‘temperature function’: Prescribe the temperature in terms of an explicit formula. The format of these functions follows the syntax understood by the muparser library.
 
 ‘velocity function’: Prescribe the velocity in terms of an explicit formula. The format of these functions follows the syntax understood by the muparser library, see {ref}`sec:run-aspect:parameters-overview:muparser-format`.
 ::::
 
+(parameters:Prescribed_20solution/Initial_20temperature)=
+## **Subsection:** Prescribed solution / Initial temperature
+::::{dropdown} __Parameter:__ {ref}`Coordinate system<parameters:Prescribed_20solution/Initial_20temperature/Coordinate_20system>`
+:name: parameters:Prescribed_20solution/Initial_20temperature/Coordinate_20system
+**Default value:** cartesian
+
+**Pattern:** [Selection cartesian|spherical|depth ]
+
+**Documentation:** A selection that determines the assumed coordinate system for the indicator function variables. Allowed values are &lsquo;cartesian&rsquo;, &lsquo;spherical&rsquo;, and &lsquo;depth&rsquo;.
+::::
+
+(parameters:Prescribed_20solution/Initial_20temperature/Indicator_20function)=
+## **Subsection:** Prescribed solution / Initial temperature / Indicator function
+::::{dropdown} __Parameter:__ {ref}`Function constants<parameters:Prescribed_20solution/Initial_20temperature/Indicator_20function/Function_20constants>`
+:name: parameters:Prescribed_20solution/Initial_20temperature/Indicator_20function/Function_20constants
+**Default value:**
+
+**Pattern:** [Anything]
+
+**Documentation:** Sometimes it is convenient to use symbolic constants in the expression that describes the function, rather than having to use its numeric value everywhere the constant appears. These values can be defined using this parameter, in the form &lsquo;var1=value1, var2=value2, ...&rsquo;.
+
+A typical example would be to set this runtime parameter to &lsquo;pi=3.1415926536&rsquo; and then use &lsquo;pi&rsquo; in the expression of the actual formula. (That said, for convenience this class actually defines both &lsquo;pi&rsquo; and &lsquo;Pi&rsquo; by default, but you get the idea.)
+::::
+
+::::{dropdown} __Parameter:__ {ref}`Function expression<parameters:Prescribed_20solution/Initial_20temperature/Indicator_20function/Function_20expression>`
+:name: parameters:Prescribed_20solution/Initial_20temperature/Indicator_20function/Function_20expression
+**Default value:** 0
+
+**Pattern:** [Anything]
+
+**Documentation:** The formula that denotes the function you want to evaluate for particular values of the independent variables. This expression may contain any of the usual operations such as addition or multiplication, as well as all of the common functions such as &lsquo;sin&rsquo; or &lsquo;cos&rsquo;. In addition, it may contain expressions like &lsquo;if(x>0, 1, -1)&rsquo; where the expression evaluates to the second argument if the first argument is true, and to the third argument otherwise. For a full overview of possible expressions accepted see the documentation of the muparser library at http://muparser.beltoforion.de/.
+
+If the function you are describing represents a vector-valued function with multiple components, then separate the expressions for individual components by a semicolon.
+::::
+
+::::{dropdown} __Parameter:__ {ref}`Variable names<parameters:Prescribed_20solution/Initial_20temperature/Indicator_20function/Variable_20names>`
+:name: parameters:Prescribed_20solution/Initial_20temperature/Indicator_20function/Variable_20names
+**Default value:** x,y,t
+
+**Pattern:** [Anything]
+
+**Documentation:** The names of the variables as they will be used in the function, separated by commas. By default, the names of variables at which the function will be evaluated are &lsquo;x&rsquo; (in 1d), &lsquo;x,y&rsquo; (in 2d) or &lsquo;x,y,z&rsquo; (in 3d) for spatial coordinates and &lsquo;t&rsquo; for time. You can then use these variable names in your function expression and they will be replaced by the values of these variables at which the function is currently evaluated. However, you can also choose a different set of names for the independent variables at which to evaluate your function expression. For example, if you work in spherical coordinates, you may wish to set this input parameter to &lsquo;r,phi,theta,t&rsquo; and then use these variable names in your function expression.
+::::
+
 (parameters:Prescribed_20solution/Temperature_20function)=
 ## **Subsection:** Prescribed solution / Temperature function
 ::::{dropdown} __Parameter:__ {ref}`Coordinate system<parameters:Prescribed_20solution/Temperature_20function/Coordinate_20system>`

include/aspect/prescribed_solution/initial_temperature.h

@@ -0,0 +1,108 @@
+/*
+  Copyright (C) 2011 - 2026 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_prescribed_solution_initial_temperature_h
+#define _aspect_prescribed_solution_initial_temperature_h
+
+#include <aspect/initial_temperature/interface.h>
+#include <aspect/prescribed_solution/interface.h>
+#include <aspect/simulator_access.h>
+#include <deal.II/base/parsed_function.h>
+
+namespace aspect
+{
+  namespace PrescribedSolution
+  {
+    /**
+     * A class that prescribes temperature in a selected region using
+     * the value returned by the active initial temperature model.
+     *
+     * The region is selected through an indicator function. Where the
+     * indicator is greater than 0.5, the temperature is constrained
+     * to the initial temperature value evaluated at that position.
+     *
+     * @ingroup PrescribedSolution
+     */
+    template <int dim>
+    class InitialTemperature
+      : public Interface<dim>,
+        public SimulatorAccess<dim>
+    {
+      public:
+        /**
+         * Constructor.
+         */
+        InitialTemperature ();
+
+        /**
+        * Store a shared pointer to the initial temperature manager so the
+        * plugin can safely access it after initialization.
+        */
+        void initialize () override;
+
+        /**
+         * Update the current time in the indicator function.
+         */
+        void update () 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;
+
+        /**
+         * Decide and assign cell-wise constraints for temperature DoFs.
+         */
+        void constrain_solution (
+          const typename DoFHandler<dim>::active_cell_iterator &cell,
+          const std::vector<Point<dim>> &positions,
+          const std::vector<unsigned int> &component_indices,
+          std::vector<bool> &should_be_constrained,
+          std::vector<double> &solution) override;
+
+      private:
+        /**
+         * Indicator function for selecting where the temperature should
+         * be prescribed.
+         */
+        Functions::ParsedFunction<dim> prescribed_temperature_indicator_function;
+
+        /**
+         * The coordinate representation used to evaluate the indicator
+         * function. Possible choices are cartesian, spherical, and depth.
+         */
+        Utilities::Coordinates::CoordinateSystem coordinate_system;
+
+        /**
+         * Shared pointer to the initial temperature manager. We keep this
+         * alive because the simulator may release its own pointer after
+         * initialization.
+         */
+        std::shared_ptr<const aspect::InitialTemperature::Manager<dim>> initial_temperature_manager;
+    };
+  }
+}
+
+#endif

source/prescribed_solution/initial_temperature.cc

@@ -0,0 +1,170 @@
+/*
+  Copyright (C) 2011 - 2026 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_temperature/interface.h>
+#include <aspect/prescribed_solution/initial_temperature.h>
+#include <aspect/geometry_model/interface.h>
+
+namespace aspect
+{
+  namespace PrescribedSolution
+  {
+    template <int dim>
+    InitialTemperature<dim>::InitialTemperature ()
+      :
+      prescribed_temperature_indicator_function(1)
+    {}
+
+
+
+    template <int dim>
+    void
+    InitialTemperature<dim>::initialize ()
+    {
+      initial_temperature_manager =
+        this->get_initial_temperature_manager_pointer();
+    }
+
+
+
+    template <int dim>
+    void
+    InitialTemperature<dim>::update ()
+    {
+      if (this->convert_output_to_years())
+        prescribed_temperature_indicator_function.set_time(this->get_time() / year_in_seconds);
+      else
+        prescribed_temperature_indicator_function.set_time(this->get_time());
+    }
+
+
+
+    template <int dim>
+    void
+    InitialTemperature<dim>::constrain_solution (const typename DoFHandler<dim>::active_cell_iterator &/*cell*/,
+                                                 const std::vector<Point<dim>> &positions,
+                                                 const std::vector<unsigned int> &component_indices,
+                                                 std::vector<bool> &should_be_constrained,
+                                                 std::vector<double> &solution)
+    {
+
+      const unsigned int temperature_index =
+        this->introspection().component_indices.temperature;
+
+      for (unsigned int q=0; q<positions.size(); ++q)
+        {
+          if (component_indices[q] != temperature_index)
+            continue;
+
+          const auto point =
+            this->get_geometry_model().cartesian_to_other_coordinates(positions[q], coordinate_system);
+
+          const double indicator =
+            prescribed_temperature_indicator_function.value(Utilities::convert_array_to_point<dim>(point.get_coordinates()), 0);
+
+          if (indicator > 0.5)
+            {
+              should_be_constrained[q] = true;
+              solution[q] = initial_temperature_manager->initial_temperature(positions[q]);
+            }
+        }
+    }
+
+
+
+    template <int dim>
+    void
+    InitialTemperature<dim>::declare_parameters (ParameterHandler &prm)
+    {
+      prm.enter_subsection("Prescribed solution");
+      {
+        prm.enter_subsection("Initial temperature");
+        {
+          prm.declare_entry("Coordinate system",
+                            "cartesian",
+                            Patterns::Selection("cartesian|spherical|depth"),
+                            "A selection that determines the assumed coordinate "
+                            "system for the indicator function variables. "
+                            "Allowed values are `cartesian', `spherical', and `depth'.");
+
+          prm.enter_subsection("Indicator function");
+          {
+            Functions::ParsedFunction<dim>::declare_parameters(prm, 1);
+          }
+          prm.leave_subsection();
+        }
+        prm.leave_subsection();
+      }
+      prm.leave_subsection();
+    }
+
+
+
+    template <int dim>
+    void
+    InitialTemperature<dim>::parse_parameters (ParameterHandler &prm)
+    {
+      prm.enter_subsection("Prescribed solution");
+      {
+        prm.enter_subsection("Initial temperature");
+        {
+          coordinate_system =
+            Utilities::Coordinates::string_to_coordinate_system(prm.get("Coordinate system"));
+
+          prm.enter_subsection("Indicator function");
+          {
+            try
+              {
+                prescribed_temperature_indicator_function.parse_parameters(prm);
+              }
+            catch (...)
+              {
+                std::cerr << "ERROR: FunctionParser failed to parse\n"
+                          << "\t'Prescribed solution.Initial temperature.Indicator function'\n"
+                          << "with expression\n"
+                          << "\t'" << prm.get("Function expression") << "'\n";
+                throw;
+              }
+          }
+          prm.leave_subsection();
+        }
+        prm.leave_subsection();
+      }
+      prm.leave_subsection();
+    }
+  }
+}
+
+
+
+// explicit instantiations
+namespace aspect
+{
+  namespace PrescribedSolution
+  {
+    ASPECT_REGISTER_PRESCRIBED_SOLUTION(InitialTemperature,
+                                        "initial temperature",
+                                        "Prescribe the temperature in a selected region using the active "
+                                        "initial temperature model. The selected region is defined through "
+                                        "an indicator function. At locations where the indicator value is greater than 0.5, "
+                                        "the temperature is constrained to the initial temperature evaluated "
+                                        "at that position.")
+  }
+}