Boundary composition plugin based on mesh deformation output

doc/sphinx/parameters/Boundary_20composition_20model.md

@@ -35,7 +35,7 @@ This parameter only describes which boundaries have a fixed composition, but not
 :name: parameters:Boundary_20composition_20model/List_20of_20model_20names
 **Default value:**
 
-**Pattern:** [MultipleSelection ascii data|box|box with lithosphere boundary indicators|function|initial composition|spherical constant ]
+**Pattern:** [MultipleSelection ascii data|box|box with lithosphere boundary indicators|function|initial composition|mesh deformation|spherical constant ]
 
 **Documentation:** A comma-separated list of boundary composition models that will be used to initialize the composition. These plugins are loaded in the order given, and modify the existing composition field via the operators listed in ’List of model operators’.
 
@@ -57,6 +57,10 @@ The format of these functions follows the syntax understood by the muparser libr
 
 Because this class simply takes what the initial composition had described, this class can not know certain pieces of information such as the minimal and maximal composition on the boundary. For operations that require this, for example in post-processing, this boundary composition model must therefore be told what the minimal and maximal values on the boundary are. This is done using parameters set in section “Boundary composition model/Initial composition”.
 
+‘mesh deformation’: A model in which the composition at the boundary is retrieved from the active mesh deformation plugins.
+
+The active mesh deformation plugins can each return a value for the compositional fields that are prescribed on a boundary; their values are summed per field. If a mesh deformation plugin does not implement the boundary composition function, a default value of zero is returned.
+
 ‘spherical constant’: A model in which the composition is chosen constant on the inner and outer boundaries of a sphere, spherical shell, chunk or ellipsoidal chunk. Parameters are read from subsection ’Spherical constant’.
 ::::
 

include/aspect/boundary_composition/mesh_deformation.h

@@ -0,0 +1,65 @@
+/*
+  Copyright (C) 2013 - 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_boundary_composition_mesh_deformation_h
+#define _aspect_boundary_composition_mesh_deformation_h
+
+#include <aspect/mesh_deformation/interface.h>
+#include <aspect/boundary_composition/interface.h>
+#include <aspect/simulator_access.h>
+
+
+namespace aspect
+{
+  namespace BoundaryComposition
+  {
+    /**
+     * A class that implements a composition boundary condition that
+     * is set by the active mesh deformation plugins. Their returned
+     * values at a given point on a boundary are summed.
+     *
+     * @ingroup BoundaryCompositions
+     */
+    template <int dim>
+    class MeshDeformation : public Interface<dim>, public SimulatorAccess<dim>
+    {
+      public:
+        /**
+         * Initialization function. This function is called once at the
+         * beginning of the program after parse_parameters is run.
+         */
+        void
+        initialize () override;
+
+        /**
+         * This function returns the composition value returned by the
+         * mesh deformation interface, which is the sum of the boundary compositions
+         * provided by the active mesh deformation objects.
+         */
+        double boundary_composition (const types::boundary_id boundary_indicator,
+                                     const Point<dim> &position,
+                                     const unsigned int compositional_field) const override;
+    };
+  }
+}
+
+
+#endif

include/aspect/mesh_deformation/fastscape.h

@@ -71,6 +71,19 @@ namespace aspect
                                                  AffineConstraints<double> &mesh_velocity_constraints,
                                                  const std::set<types::boundary_id> &boundary_ids) const override;
 
+        /**
+         * This function returns the compositional field value at the
+         * requested point on the requested boundary. If the fastscape
+         * plugin cannot provide a boundary condition for a field (e.g.
+         * for a field storing viscoelastic stresses), a value of zero
+         * is returned.
+         *
+         * @copydoc aspect::MeshDeformation::Interface::boundary_composition()
+         */
+        double boundary_composition (const types::boundary_id boundary_indicator,
+                                     const Point<dim> &position,
+                                     const unsigned int compositional_field) const override;
+
         /**
          * Returns whether or not the plugin requires surface stabilization
          */

include/aspect/mesh_deformation/interface.h

@@ -137,6 +137,26 @@ namespace aspect
         compute_velocity_constraints_on_boundary(const DoFHandler<dim> &mesh_deformation_dof_handler,
                                                  AffineConstraints<double> &mesh_velocity_constraints,
                                                  const std::set<types::boundary_id> &boundary_ids) const;
+
+        /**
+         * Return the composition that is to hold at a particular position on
+         * the boundary of the domain.
+         *
+         * @param boundary_indicator The boundary indicator of the part of the
+         * boundary of the domain on which the point is located at which we
+         * are requesting the composition.
+         * @param position The position of the point at which we ask for the
+         * composition.
+         * @param compositional_field The index of the compositional field
+         * between 0 and @p parameters.n_compositional_fields.
+         * @return Boundary value of the compositional field @p
+         * compositional_field at the position @p position.
+         */
+        virtual
+        double
+        boundary_composition (const types::boundary_id boundary_indicator,
+                              const Point<dim> &position,
+                              const unsigned int compositional_field) const;
     };
 
 
@@ -378,6 +398,18 @@ namespace aspect
         const Mapping<dim> &
         get_level_mapping(const unsigned int level) const;
 
+        /**
+         * Loop over all mesh deformation objects that are active on the boundary
+         * with the ID boundary_indicator and sum their contributions to the compositional field
+         * with index compositional_field at the given position.
+         * For example, the plugin might want to set the field representing sediment to 1
+         * if deposition occurs, or a field representing deposition depth to the depth below sea level.
+         */
+        double
+        boundary_composition (const types::boundary_id boundary_indicator,
+                              const Point<dim> &position,
+                              const unsigned int compositional_field) const;
+
         /**
          * For the current plugin subsystem, write a connection graph of all of the
          * plugins we know about, in the format that the

source/boundary_composition/mesh_deformation.cc

@@ -0,0 +1,73 @@
+/*
+  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/boundary_composition/mesh_deformation.h>
+#include <aspect/initial_composition/interface.h>
+
+
+namespace aspect
+{
+  namespace BoundaryComposition
+  {
+// ------------------------------ MeshDeformation -------------------
+
+    template <int dim>
+    void
+    MeshDeformation<dim>::initialize()
+    {
+      // Check that mesh deformation is actually active. The mesh deformation plugins
+      // are not initialized yet, so we can only check the input parameters.
+      AssertThrow(this->get_parameters().mesh_deformation_enabled == true,
+                  ExcMessage ("The boundary composition plugin ``mesh deformation'' can only be used when a mesh deformation plugin is active."));
+    }
+
+
+
+    template <int dim>
+    double
+    MeshDeformation<dim>::
+    boundary_composition (const types::boundary_id boundary_indicator,
+                          const Point<dim> &position,
+                          const unsigned int compositional_field) const
+    {
+      return this->get_mesh_deformation_handler().boundary_composition(boundary_indicator, position, compositional_field);
+    }
+
+  }
+}
+
+// explicit instantiations
+namespace aspect
+{
+  namespace BoundaryComposition
+  {
+    ASPECT_REGISTER_BOUNDARY_COMPOSITION_MODEL(MeshDeformation,
+                                               "mesh deformation",
+                                               "A model in which the composition at the boundary "
+                                               "is retrieved from the active mesh deformation plugins."
+                                               "\n\n"
+                                               "The active mesh deformation plugins can each return "
+                                               "a value for the compositional fields that are prescribed on "
+                                               "a boundary; their values are summed per field. If a mesh "
+                                               "deformation plugin does not implement the boundary composition "
+                                               "function, a default value of zero is returned.")
+  }
+}

source/mesh_deformation/fastscape.cc

@@ -1591,6 +1591,56 @@ namespace aspect
 
 
 
+    template <int dim>
+    double
+    FastScape<dim>::
+    boundary_composition (const types::boundary_id boundary_indicator,
+                          const Point<dim> &position,
+                          const unsigned int compositional_field) const
+    {
+      // FastScape is only applied to the top boundary of the model domain.
+      // If a composition value is requested for any other boundary,
+      // return zero.
+      if (boundary_indicator != this->get_geometry_model().translate_symbolic_boundary_name_to_id ("top"))
+        return 0.0;
+
+      // Two fields often used in conjunction with the FastScape plugin
+      // are sediment_age and deposition_depth. If the fields exist, and
+      // their boundary values are requested, set them here.
+      if ( this->introspection().compositional_name_exists("sediment_age") &&
+           compositional_field == this->introspection().compositional_index_for_name("sediment_age"))
+        {
+          return this->get_parameters().convert_to_years ? this->get_time()/year_in_seconds : this->get_time();
+        }
+      else if ( this->introspection().compositional_name_exists("deposition_depth") &&
+                compositional_field == this->introspection().compositional_index_for_name("deposition_depth"))
+        {
+          // Get the time-dependent sea level if necessary.
+          const double current_sea_level = use_sea_level_function
+                                           ? sea_level_function.value(Point<1>())
+                                           : sea_level_constant_value;
+
+          // FastScape only works on box geometries, so the last component of the position is the height.
+          // The sea level is defined with respect to the original, unperturbed height of the box.
+          // Sediments deposited below sea level will have a positive deposition depth.
+          const GeometryModel::Box<dim> *box_geometry
+            = dynamic_cast<const GeometryModel::Box<dim>*> (&this->get_geometry_model());
+          const GeometryModel::TwoMergedBoxes<dim> *two_merged_boxes_geometry
+            = dynamic_cast<const GeometryModel::TwoMergedBoxes<dim>*> (&this->get_geometry_model());
+          const Point<dim> origin = (box_geometry != nullptr
+                                     ? box_geometry->get_origin()
+                                     : two_merged_boxes_geometry->get_origin());
+          const Point<dim> extents = (box_geometry != nullptr
+                                      ? box_geometry->get_extents()
+                                      : two_merged_boxes_geometry->get_extents());
+          return origin[dim-1] + extents[dim-1] + current_sea_level - position[dim-1];
+        }
+      else
+        return 0.0;
+    }
+
+
+
     template <int dim>
     template <class Archive>
     void FastScape<dim>::serialize (Archive &ar, const unsigned int)

source/mesh_deformation/interface.cc

@@ -217,6 +217,19 @@ namespace aspect
 
 
 
+    template <int dim>
+    double
+    Interface<dim>::
+    boundary_composition (const types::boundary_id /*boundary_indicator*/,
+                          const Point<dim> &/*position*/,
+                          const unsigned int /*compositional_field*/) const
+    {
+      AssertThrow(false, ExcMessage("The boundary_composition function is not implemented for this mesh deformation plugin."));
+      return 0.0;
+    }
+
+
+
     template <int dim>
     MeshDeformationHandler<dim>::MeshDeformationHandler (Simulator<dim> &simulator)
       : sim(simulator),  // reference to the simulator that owns the MeshDeformationHandler
@@ -643,6 +656,26 @@ namespace aspect
 
 
 
+    template <int dim>
+    double
+    MeshDeformationHandler<dim>::boundary_composition (const types::boundary_id boundary_indicator,
+                                                       const Point<dim> &position,
+                                                       const unsigned int compositional_field) const
+    {
+      double composition = 0.0;
+
+      // Loop over all mesh deformation objects that are assigned to this
+      // boundary indicator and sum their contributions
+      for (const auto &deformation_object : mesh_deformation_objects.at(boundary_indicator))
+        {
+          composition += deformation_object->boundary_composition(boundary_indicator,
+                                                                  position,
+                                                                  compositional_field);
+        }
+
+      return composition;
+    }
+
     template <int dim>
     void MeshDeformationHandler<dim>::make_constraints()
     {