Prescribed condition from initial composition

lhy11009 · lhy11009 · commit ec34707012fe · 2026-03-23T07:33:18.000-07:00
diff --git a/doc/modules/changes/20260311_lhy11009 b/doc/modules/changes/20260311_lhy11009
@@ -0,0 +1,5 @@
+Added: a new instance of 'initial composition' in the prescribed solution plugin system.
+This plugin prescribes compositional fields from the initial composition model within a
+region defined by an indicator function.
+<br>
+(Haoyuan Li, 2026/03/11)
diff --git a/doc/sphinx/parameters/Prescribed_20solution.md b/doc/sphinx/parameters/Prescribed_20solution.md
@@ -9,19 +9,65 @@
 :name: parameters:Prescribed_20solution/List_20of_20model_20names
 **Default value:**
 
-**Pattern:** [MultipleSelection initial temperature|temperature function|velocity function ]
+**Pattern:** [MultipleSelection initial composition|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 &rsquo;List of model operators&rsquo;.
 
 The following prescribed solution models are available:
 
+&lsquo;initial composition&rsquo;: Prescribe compositional fields from the initial composition model within a region defined by an indicator function.
+
 &lsquo;initial temperature&rsquo;: 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.
 
 &lsquo;temperature function&rsquo;: Prescribe the temperature in terms of an explicit formula. The format of these functions follows the syntax understood by the muparser library.
 
 &lsquo;velocity function&rsquo;: 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_20composition)=
+## **Subsection:** Prescribed solution / Initial composition
+::::{dropdown} __Parameter:__ {ref}`Coordinate system<parameters:Prescribed_20solution/Initial_20composition/Coordinate_20system>`
+:name: parameters:Prescribed_20solution/Initial_20composition/Coordinate_20system
+**Default value:** cartesian
+
+**Pattern:** [Selection cartesian|spherical|depth ]
+
+**Documentation:** Coordinate system used for evaluating the indicator function.
+::::
+
+(parameters:Prescribed_20solution/Initial_20composition/Indicator_20function)=
+## **Subsection:** Prescribed solution / Initial composition / Indicator function
+::::{dropdown} __Parameter:__ {ref}`Function constants<parameters:Prescribed_20solution/Initial_20composition/Indicator_20function/Function_20constants>`
+:name: parameters:Prescribed_20solution/Initial_20composition/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_20composition/Indicator_20function/Function_20expression>`
+:name: parameters:Prescribed_20solution/Initial_20composition/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_20composition/Indicator_20function/Variable_20names>`
+:name: parameters:Prescribed_20solution/Initial_20composition/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/Initial_20temperature)=
 ## **Subsection:** Prescribed solution / Initial temperature
 ::::{dropdown} __Parameter:__ {ref}`Coordinate system<parameters:Prescribed_20solution/Initial_20temperature/Coordinate_20system>`
diff --git a/include/aspect/prescribed_solution/initial_composition.h b/include/aspect/prescribed_solution/initial_composition.h
@@ -0,0 +1,98 @@
+/*
+ Copyright (C) 2015 - 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/>.
+*/
+
+#ifndef _aspect_prescribed_solution_initial_composition_h
+#define _aspect_prescribed_solution_initial_composition_h
+
+#include <aspect/initial_composition/interface.h>
+#include <aspect/prescribed_solution/interface.h>
+#include <aspect/simulator_access.h>
+
+#include <deal.II/base/parsed_function.h>
+
+namespace aspect
+{
+  namespace PrescribedSolution
+  {
+    /**
+     * Prescribe compositional fields using the initial composition model.
+     * The region where the constraint applies is determined by an
+     * indicator function.
+     */
+    template <int dim>
+    class InitialComposition
+      : public Interface<dim>,
+        public SimulatorAccess<dim>
+    {
+      public:
+
+        InitialComposition ();
+
+        /**
+         * Store shared pointer to initial composition manager.
+         */
+        void initialize () override;
+
+        /**
+         * Update function time.
+         */
+        void update () override;
+
+        /**
+         * Declare parameters.
+         */
+        static void declare_parameters (ParameterHandler &prm);
+
+        /**
+         * Parse parameters.
+         */
+        void parse_parameters (ParameterHandler &prm) override;
+
+        /**
+         * Constrain compositional solution.
+         */
+        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 defining region where composition
+         * is prescribed.
+         */
+        Functions::ParsedFunction<dim> indicator_function;
+
+        /**
+         * Coordinate system used for evaluating indicator.
+         */
+        Utilities::Coordinates::CoordinateSystem coordinate_system;
+
+        /**
+         * Pointer to initial composition manager.
+         */
+        std::shared_ptr<const aspect::InitialComposition::Manager<dim>> initial_composition_manager;
+    };
+  }
+}
+
+#endif
diff --git a/source/prescribed_solution/initial_composition.cc b/source/prescribed_solution/initial_composition.cc
@@ -0,0 +1,175 @@
+/*
+  Copyright (C) 2011 - 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 <aspect/initial_composition/interface.h>
+#include <aspect/prescribed_solution/initial_composition.h>
+
+namespace aspect
+{
+  namespace PrescribedSolution
+  {
+
+    template <int dim>
+    InitialComposition<dim>::InitialComposition ()
+      :
+      indicator_function(1)
+    {}
+
+
+
+    template <int dim>
+    void
+    InitialComposition<dim>::initialize ()
+    {
+      initial_composition_manager =
+        this->get_initial_composition_manager_pointer();
+    }
+
+
+
+    template <int dim>
+    void
+    InitialComposition<dim>::update ()
+    {
+      if (this->convert_output_to_years())
+        indicator_function.set_time(this->get_time() / year_in_seconds);
+      else
+        indicator_function.set_time(this->get_time());
+    }
+
+
+
+    template <int dim>
+    void
+    InitialComposition<dim>::constrain_solution (const typename DoFHandler<dim>::active_cell_iterator &,
+                                                 const std::vector<Point<dim>> &positions,
+                                                 const std::vector<unsigned int> &component_indices,
+                                                 std::vector<bool> &should_be_constrained,
+                                                 std::vector<double> &solution)
+    {
+      // Determine the component range corresponding to compositional fields.
+      // These component indices are originally mapped from local DoFs and
+      // determine which DoFs in the system belong to compositional fields.
+      const unsigned int first_comp =
+        this->introspection().component_indices.compositional_fields[0];
+
+      const unsigned int last_comp =
+        this->introspection().component_indices.compositional_fields[this->introspection().n_compositional_fields-1];
+
+      const auto &geometry_model = this->get_geometry_model();
+
+      for (unsigned int q=0; q<positions.size(); ++q)
+        {
+          const unsigned int component = component_indices[q];
+
+          if (component < first_comp || component > last_comp)
+            continue;
+
+          const auto point = geometry_model.cartesian_to_other_coordinates(
+                               positions[q], coordinate_system);
+
+          const double indicator = indicator_function.value(
+                                     Utilities::convert_array_to_point<dim>(point.get_coordinates()));
+
+          if (indicator > 0.5)
+            {
+              const unsigned int field_index = component - first_comp;
+
+              solution[q] =
+                initial_composition_manager->initial_composition(
+                  positions[q], field_index);
+
+              should_be_constrained[q] = true;
+            }
+        }
+    }
+
+
+
+    template <int dim>
+    void
+    InitialComposition<dim>::declare_parameters (ParameterHandler &prm)
+    {
+      prm.enter_subsection("Prescribed solution");
+      {
+        prm.enter_subsection("Initial composition");
+        {
+
+          prm.declare_entry("Coordinate system",
+                            "cartesian",
+                            Patterns::Selection("cartesian|spherical|depth"),
+                            "Coordinate system used for evaluating the indicator function.");
+
+          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
+    InitialComposition<dim>::parse_parameters (ParameterHandler &prm)
+    {
+      prm.enter_subsection("Prescribed solution");
+      {
+        prm.enter_subsection("Initial composition");
+        {
+
+          coordinate_system =
+            Utilities::Coordinates::string_to_coordinate_system(
+              prm.get("Coordinate system"));
+
+          prm.enter_subsection("Indicator function");
+          {
+            indicator_function.parse_parameters(prm);
+          }
+          prm.leave_subsection();
+
+        }
+        prm.leave_subsection();
+      }
+      prm.leave_subsection();
+    }
+  }
+}
+
+
+
+namespace aspect
+{
+  namespace PrescribedSolution
+  {
+
+    ASPECT_REGISTER_PRESCRIBED_SOLUTION(InitialComposition,
+                                        "initial composition",
+                                        "Prescribe compositional fields from the initial composition model "
+                                        "within a region defined by an indicator function.")
+
+  }
+}
diff --git a/tests/prescribed_solution_composition_initial_condition.prm b/tests/prescribed_solution_composition_initial_condition.prm
@@ -0,0 +1,41 @@
+#########################################################
+# This is modified from the composition_passive.prm
+# parameter file by prescribing the compositional solution
+# from initial compositional fields.
+
+include $ASPECT_SOURCE_DIR/cookbooks/composition_passive/composition_passive.prm
+
+# Modify the original cookbook to end earlier
+set Start time                             = 0
+set End time                               = 1
+set Use years instead of seconds           = false
+
+# Modify the original cookbook to have coarser meshes
+subsection Mesh refinement
+  set Initial adaptive refinement        = 0
+  set Initial global refinement          = 3
+  set Time steps between mesh refinement = 0
+end
+
+# The next section prescribes the composition inside the domain
+# using a spatially and temporally varying function. An indicator
+# function selects the region where the composition is constrained.
+subsection Prescribed solution
+    set List of model names = initial composition
+    subsection Initial composition
+        subsection Indicator function
+            set Variable names = x, y
+            set Function expression = (x<1) ? 1:0
+        end
+    end
+end
+
+subsection Postprocess
+  set List of postprocessors = visualization
+
+  subsection Visualization
+    set Interpolate output = false
+    set Time between graphical output = 0.1
+    set Output format = gnuplot
+  end
+end
diff --git a/tests/prescribed_solution_composition_initial_condition/screen-output b/tests/prescribed_solution_composition_initial_condition/screen-output
