Combine interpolating fields in one operation

doc/modules/changes/20250113_gassmoeller

@@ -0,0 +1,5 @@
+Changed: Prescribed compositional fields are now copied from the material
+model output in one combined operation (instead of being copied field by field).
+This change improves the efficiency of the copy process.
+<br>
+(Rene Gassmoeller, 2025/01/13)

include/aspect/simulator.h

@@ -1370,6 +1370,25 @@ namespace aspect
        */
       void apply_limiter_to_dg_solutions (const AdvectionField &advection_field);
 
+      /**
+       * Compute the unique support points for the advection fields @p advection_fields.
+       * The support points are collected by taking the union of the support points
+       * of all given advection fields and filtering out duplicate points. The resulting
+       * set of points is written into @p unique_support_points. @p support_point_index_by_field
+       * is a vector of vectors that contains the indices of the support points for each field.
+       * I.e. support_point_index_by_field[i][j] contains the j-th support point
+       * for the i-th field in @p advection_fields and
+       * unique_support_points[support_point_index_by_field[i][j]] is its location.
+       *
+       * For the common case that all fields have the same support points, each vector in
+       * @p support_point_index_by_field will contain the same indices for all fields.
+       *
+       * Note that existing content of @p unique_support_points and @p support_point_index_by_field
+       * will be overwritten in this function.
+       */
+      void compute_unique_advection_support_points (const std::vector<AdvectionField> &advection_fields,
+                                                    std::vector<Point<dim>> &unique_support_points,
+                                                    std::vector<std::vector<unsigned int>> &support_point_index_by_field) const;
 
       /**
        * Compute the reactions in case of operator splitting:
@@ -1439,7 +1458,7 @@ namespace aspect
        * This function is implemented in
        * <code>source/simulator/helper_functions.cc</code>.
        */
-      void interpolate_material_output_into_advection_field (const AdvectionField &adv_field);
+      void interpolate_material_output_into_advection_field (const std::vector<AdvectionField> &adv_field);
 
 
       /**

source/simulator/solver_schemes.cc

@@ -157,7 +157,7 @@ namespace aspect
             {
               TimerOutput::Scope timer (computing_timer, "Interpolate prescribed temperature");
 
-              interpolate_material_output_into_advection_field(adv_field);
+              interpolate_material_output_into_advection_field({adv_field});
 
               // Also set the old_solution block to the prescribed field. The old
               // solution is the one that is used to assemble the diffusion system in
@@ -180,14 +180,15 @@ namespace aspect
 
           TimerOutput::Scope timer (computing_timer, "Interpolate prescribed temperature");
 
-          interpolate_material_output_into_advection_field(adv_field);
+          interpolate_material_output_into_advection_field({adv_field});
 
           // Call the signal in case the user wants to do something with the variable:
           SolverControl dummy;
           signals.post_advection_solver(*this,
                                         adv_field.is_temperature(),
                                         adv_field.compositional_variable,
                                         dummy);
+
           break;
         }
 
@@ -245,6 +246,7 @@ namespace aspect
       Assert(residual->size() == introspection.n_compositional_fields, ExcInternalError());
 
     std::vector<AdvectionField> fields_advected_by_particles;
+    std::vector<AdvectionField> fields_interpolated_from_material_output;
 
     for (unsigned int c=0; c < introspection.n_compositional_fields; ++c)
       {
@@ -263,7 +265,7 @@ namespace aspect
                 {
                   TimerOutput::Scope timer (computing_timer, "Interpolate prescribed composition");
 
-                  interpolate_material_output_into_advection_field(adv_field);
+                  interpolate_material_output_into_advection_field({adv_field});
 
                   // Also set the old_solution block to the prescribed field. The old
                   // solution is the one that is used to assemble the diffusion system in
@@ -298,16 +300,8 @@ namespace aspect
 
             case Parameters<dim>::AdvectionFieldMethod::prescribed_field:
             {
-              TimerOutput::Scope timer (computing_timer, "Interpolate prescribed composition");
-
-              interpolate_material_output_into_advection_field(adv_field);
-
-              // Call the signal in case the user wants to do something with the variable:
-              SolverControl dummy;
-              signals.post_advection_solver(*this,
-                                            adv_field.is_temperature(),
-                                            adv_field.compositional_variable,
-                                            dummy);
+              // handle all prescribed fields together to increase efficiency
+              fields_interpolated_from_material_output.push_back(adv_field);
               break;
             }
 
@@ -334,6 +328,23 @@ namespace aspect
           }
       }
 
+    if (fields_interpolated_from_material_output.size() > 0)
+      {
+        TimerOutput::Scope timer (computing_timer, "Interpolate prescribed composition");
+
+        interpolate_material_output_into_advection_field(fields_interpolated_from_material_output);
+
+        for (const auto &adv_field: fields_interpolated_from_material_output)
+          {
+            // Call the signal in case the user wants to do something with the variable:
+            SolverControl dummy;
+            signals.post_advection_solver(*this,
+                                          adv_field.is_temperature(),
+                                          adv_field.compositional_variable,
+                                          dummy);
+          }
+      }
+
     if (fields_advected_by_particles.size() > 0)
       interpolate_particle_properties(fields_advected_by_particles);
 

tests/entropy_half_space/statistics

@@ -8,8 +8,8 @@
 # 8: Number of nonlinear iterations
 # 9: Iterations for temperature solver
 # 10: Iterations for composition solver 1
-# 11: Iterations for composition solver 2
-# 12: Iterations for composition solver 3
+# 11: Iterations for composition solver 3
+# 12: Iterations for composition solver 2
 # 13: Iterations for Stokes solver
 # 14: Velocity iterations in Stokes preconditioner
 # 15: Schur complement iterations in Stokes preconditioner

tests/multiple_prescribed_fields.cc

@@ -0,0 +1,95 @@
+/*
+  Copyright (C) 2022 - 2024 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/material_model/simple.h>
+#include <aspect/simulator_access.h>
+
+namespace aspect
+{
+  namespace MaterialModel
+  {
+    using namespace dealii;
+
+    template <int dim>
+    class PrescribedFieldsMaterial : public MaterialModel::Simple<dim>
+    {
+      public:
+
+        void evaluate(const MaterialModel::MaterialModelInputs<dim> &in,
+                      MaterialModel::MaterialModelOutputs<dim> &out) const override;
+
+        void
+        create_additional_named_outputs (MaterialModel::MaterialModelOutputs<dim> &out) const override;
+    };
+
+  }
+}
+
+namespace aspect
+{
+  namespace MaterialModel
+  {
+
+    template <int dim>
+    void
+    PrescribedFieldsMaterial<dim>::
+    evaluate(const MaterialModel::MaterialModelInputs<dim> &in,
+             MaterialModel::MaterialModelOutputs<dim> &out) const
+    {
+      Simple<dim>::evaluate(in, out);
+
+      // set up variable to interpolate prescribed field outputs onto compositional fields
+      PrescribedFieldOutputs<dim> *prescribed_field_out = out.template get_additional_output<PrescribedFieldOutputs<dim>>();
+
+      if (prescribed_field_out != nullptr)
+        for (unsigned int i=0; i < in.n_evaluation_points(); ++i)
+          {
+            const double y = in.position[i](1);
+            for (unsigned int j=0; j<this->n_compositional_fields(); ++j)
+              prescribed_field_out->prescribed_field_outputs[i][j] = y*j;
+          }
+    }
+
+
+    template <int dim>
+    void
+    PrescribedFieldsMaterial<dim>::create_additional_named_outputs (MaterialModel::MaterialModelOutputs<dim> &out) const
+    {
+      if (out.template get_additional_output<PrescribedFieldOutputs<dim>>() == nullptr)
+        {
+          const unsigned int n_points = out.n_evaluation_points();
+          out.additional_outputs.push_back(
+            std::make_unique<MaterialModel::PrescribedFieldOutputs<dim>> (n_points,this->n_compositional_fields()));
+        }
+    }
+  }
+}
+
+// explicit instantiations
+namespace aspect
+{
+  namespace MaterialModel
+  {
+    ASPECT_REGISTER_MATERIAL_MODEL(PrescribedFieldsMaterial,
+                                   "prescribed fields material",
+                                   "A simple material model that is like the "
+                                   "'Simple' model, but creates multiple prescribed field outputs.")
+  }
+}

tests/multiple_prescribed_fields.prm

@@ -0,0 +1,83 @@
+# A testcase that demonstrates the named additional outputs postprocessor for multiple
+# named outputs. The test is based on prescribed_field_with_diffusion only
+# with more than one additional output.
+
+set Dimension = 2
+set End time                               = 0
+set Start time                             = 0
+set Adiabatic surface temperature          = 0
+set Surface pressure                       = 0
+set Use years in output instead of seconds = false
+
+subsection Gravity model
+  set Model name = vertical
+
+  subsection Vertical
+    set Magnitude = 1.0
+  end
+end
+
+subsection Geometry model
+  set Model name = box
+
+  subsection Box
+    set X extent = 1
+    set Y extent = 1
+  end
+end
+
+subsection Compositional fields
+  set Number of fields = 3
+  set Names of fields  = prescribed_1, prescribed_2, prescribed_3
+  set Compositional field methods = prescribed field, prescribed field, prescribed field
+end
+
+subsection Initial temperature model
+  set Model name = function
+
+  subsection Function
+    set Function expression = 0.1
+  end
+end
+
+subsection Initial composition model
+  set Model name = function
+
+  subsection Function
+    set Function expression = 0.0; 0.2; 0.3
+  end
+end
+
+subsection Material model
+  set Model name = prescribed fields material
+
+  subsection Simple model
+    set Reference density             = 1
+    set Reference specific heat       = 1
+    set Reference temperature         = 0
+    set Thermal conductivity          = 1
+    set Thermal expansion coefficient = 1
+    set Viscosity                     = 1
+  end
+end
+
+subsection Mesh refinement
+  set Initial adaptive refinement        = 0
+  set Initial global refinement          = 1
+end
+
+subsection Boundary velocity model
+  set Tangential velocity boundary indicators = 0, 1, 2, 3
+end
+
+subsection Postprocess
+  set List of postprocessors = visualization, composition statistics
+
+  subsection Visualization
+    set Interpolate output = false
+    set List of output variables      = named additional outputs
+    set Number of grouped files       = 0
+    set Output format                 = gnuplot
+    set Time between graphical output = 0
+  end
+end

tests/multiple_prescribed_fields/screen-output

@@ -0,0 +1,19 @@
+
+Loading shared library <./libmultiple_prescribed_fields.debug.so>
+
+Number of active cells: 4 (on 2 levels)
+Number of degrees of freedom: 159 (50+9+25+25+25+25)
+
+*** Timestep 0:  t=0 seconds, dt=0 seconds
+   Solving temperature system... 0 iterations.
+   Copying properties into prescribed compositional fields... done.
+   Solving Stokes system (GMG)... 6+0 iterations.
+
+   Postprocessing:
+     Writing graphical output:  output-multiple_prescribed_fields/solution/solution-00000
+     Compositions min/max/mass: 0/0/0 // 0/1/0.5 // 0/2/1
+
+Termination requested by criterion: end time
+
+
+

tests/multiple_prescribed_fields/statistics

@@ -0,0 +1,25 @@
+# 1: Time step number
+# 2: Time (seconds)
+# 3: Time step size (seconds)
+# 4: Number of mesh cells
+# 5: Number of Stokes degrees of freedom
+# 6: Number of temperature degrees of freedom
+# 7: Number of degrees of freedom for all compositions
+# 8: Iterations for temperature solver
+# 9: Iterations for composition solver 1
+# 10: Iterations for composition solver 2
+# 11: Iterations for composition solver 3
+# 12: Iterations for Stokes solver
+# 13: Velocity iterations in Stokes preconditioner
+# 14: Schur complement iterations in Stokes preconditioner
+# 15: Visualization file name
+# 16: Minimal value for composition prescribed_1
+# 17: Maximal value for composition prescribed_1
+# 18: Global mass for composition prescribed_1
+# 19: Minimal value for composition prescribed_2
+# 20: Maximal value for composition prescribed_2
+# 21: Global mass for composition prescribed_2
+# 22: Minimal value for composition prescribed_3
+# 23: Maximal value for composition prescribed_3
+# 24: Global mass for composition prescribed_3
+0 0.000000000000e+00 0.000000000000e+00 4 59 25 75 0 4294967294 4294967294 4294967294 5 7 7 output-multiple_prescribed_fields/solution/solution-00000 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 5.00000000e-01 0.00000000e+00 2.00000000e+00 1.00000000e+00