keep support for the old method

geodynamics · Jun 14, 2024 · 0999eb5 · 0999eb5
1 parent 314f518
commit 0999eb5
Show file tree

Hide file tree

Showing 15 changed files with 198 additions and 66 deletions.
diff --git a/include/aspect/parameters.h b/include/aspect/parameters.h
@@ -381,6 +381,37 @@ namespace aspect
       }
     };
 
+    /**
+     * This enum represents the different choices for the reaction solver.
+     * See @p reaction_solver_type.
+     */
+    struct ReactionSolverType
+    {
+      enum Kind
+      {
+        ARKode,
+        fixed_step
+      };
+
+      static const std::string pattern()
+      {
+        return "ARKode|fixed step";
+      }
+
+      static Kind
+      parse(const std::string &input)
+      {
+        if (input == "ARKode")
+          return ARKode;
+        else if (input == "fixed step")
+          return fixed_step;
+        else
+          AssertThrow(false, ExcNotImplemented());
+
+        return Kind();
+      }
+    };
+
     /**
      * Use the struct aspect::CompositionalFieldDescription
      */
@@ -512,6 +543,8 @@ namespace aspect
     bool                           AMG_output_details;
 
     // subsection: Operator splitting parameters
+    typename ReactionSolverType::Kind reaction_solver_type;
+    double                         ARKode_relative_tolerance;
     double                         reaction_time_step;
     unsigned int                   reaction_steps_per_advection_step;
 

diff --git a/source/material_model/melt_boukare.cc b/source/material_model/melt_boukare.cc
@@ -1130,8 +1130,10 @@ namespace aspect
           if (this->convert_output_to_years() == true)
             melting_time_scale *= year_in_seconds;
 
-          AssertThrow(this->get_parameters().use_operator_splitting,
-                      ExcMessage("The melt boukare material model has to be used with oprator splitting."));
+          AssertThrow(this->get_parameters().use_operator_splitting &&
+                      this->get_parameters().reaction_solver_type == Parameters<dim>::ReactionSolverType::fixed_step,
+                      ExcMessage("The melt boukare material model has to be used with operator splitting, "
+                                 "and the reaction solver needs to be `fixed step'."));
 
           AssertThrow(melting_time_scale >= this->get_parameters().reaction_time_step,
                       ExcMessage("The reaction time step " + Utilities::to_string(this->get_parameters().reaction_time_step)

diff --git a/source/material_model/melt_global.cc b/source/material_model/melt_global.cc
@@ -488,12 +488,13 @@ namespace aspect
 
           if (this->get_parameters().use_operator_splitting)
             {
-              AssertThrow(melting_time_scale >= this->get_parameters().reaction_time_step,
-                          ExcMessage("The reaction time step " + Utilities::to_string(this->get_parameters().reaction_time_step)
-                                     + " in the operator splitting scheme is too large to compute melting rates! "
-                                     "You have to choose it in such a way that it is smaller than the 'Melting time scale for "
-                                     "operator splitting' chosen in the material model, which is currently "
-                                     + Utilities::to_string(melting_time_scale) + "."));
+              if (this->get_parameters().reaction_solver_type == Parameters<dim>::ReactionSolverType::fixed_step)
+                AssertThrow(melting_time_scale >= this->get_parameters().reaction_time_step,
+                            ExcMessage("The reaction time step " + Utilities::to_string(this->get_parameters().reaction_time_step)
+                                       + " in the operator splitting scheme is too large to compute melting rates! "
+                                       "You have to choose it in such a way that it is smaller than the 'Melting time scale for "
+                                       "operator splitting' chosen in the material model, which is currently "
+                                       + Utilities::to_string(melting_time_scale) + "."));
               AssertThrow(melting_time_scale > 0,
                           ExcMessage("The Melting time scale for operator splitting must be larger than 0!"));
               AssertThrow(this->introspection().compositional_name_exists("porosity"),

diff --git a/source/material_model/reaction_model/katz2003_mantle_melting.cc b/source/material_model/reaction_model/katz2003_mantle_melting.cc
@@ -608,20 +608,25 @@ namespace aspect
             freezing_rate /= year_in_seconds;
           }
 
-        AssertThrow(melting_time_scale >= this->get_parameters().reaction_time_step,
-                    ExcMessage("The reaction time step " + Utilities::to_string(this->get_parameters().reaction_time_step)
-                               + " in the operator splitting scheme is too large to compute melting rates! "
-                               "You have to choose it in such a way that it is smaller than the 'Melting time scale for "
-                               "operator splitting' chosen in the material model, which is currently "
-                               + Utilities::to_string(melting_time_scale) + "."));
         AssertThrow(melting_time_scale > 0,
                     ExcMessage("The Melting time scale for operator splitting must be larger than 0!"));
-        AssertThrow(freezing_rate * this->get_parameters().reaction_time_step <= 1.0,
-                    ExcMessage("The reaction time step " + Utilities::to_string(this->get_parameters().reaction_time_step)
-                               + " in the operator splitting scheme is too large to compute freezing rates! "
-                               "You have to choose it in such a way that it is smaller than the inverse of the "
-                               "'Freezing rate' chosen in the material model, which is currently "
-                               + Utilities::to_string(1.0/freezing_rate) + "."));
+
+        if (this->get_parameters().reaction_solver_type == Parameters<dim>::ReactionSolverType::fixed_step)
+          {
+            AssertThrow(melting_time_scale >= this->get_parameters().reaction_time_step,
+                        ExcMessage("The reaction time step " + Utilities::to_string(this->get_parameters().reaction_time_step)
+                                   + " in the operator splitting scheme is too large to compute melting rates! "
+                                   "You have to choose it in such a way that it is smaller than the 'Melting time scale for "
+                                   "operator splitting' chosen in the material model, which is currently "
+                                   + Utilities::to_string(melting_time_scale) + "."));
+
+            AssertThrow(freezing_rate * this->get_parameters().reaction_time_step <= 1.0,
+                        ExcMessage("The reaction time step " + Utilities::to_string(this->get_parameters().reaction_time_step)
+                                   + " in the operator splitting scheme is too large to compute freezing rates! "
+                                   "You have to choose it in such a way that it is smaller than the inverse of the "
+                                   "'Freezing rate' chosen in the material model, which is currently "
+                                   + Utilities::to_string(1.0/freezing_rate) + "."));
+          }
       }
     }
   }

diff --git a/source/material_model/reactive_fluid_transport.cc b/source/material_model/reactive_fluid_transport.cc
@@ -509,12 +509,13 @@ namespace aspect
 
           if (this->get_parameters().use_operator_splitting)
             {
-              AssertThrow(fluid_reaction_time_scale >= this->get_parameters().reaction_time_step,
-                          ExcMessage("The reaction time step " + Utilities::to_string(this->get_parameters().reaction_time_step)
-                                     + " in the operator splitting scheme is too large to compute fluid release rates! "
-                                     "You have to choose it in such a way that it is smaller than the 'Fluid reaction time scale for "
-                                     "operator splitting' chosen in the material model, which is currently "
-                                     + Utilities::to_string(fluid_reaction_time_scale) + "."));
+              if (this->get_parameters().reaction_solver_type == Parameters<dim>::ReactionSolverType::fixed_step)
+                AssertThrow(fluid_reaction_time_scale >= this->get_parameters().reaction_time_step,
+                            ExcMessage("The reaction time step " + Utilities::to_string(this->get_parameters().reaction_time_step)
+                                       + " in the operator splitting scheme is too large to compute fluid release rates! "
+                                       "You have to choose it in such a way that it is smaller than the 'Fluid reaction time scale for "
+                                       "operator splitting' chosen in the material model, which is currently "
+                                       + Utilities::to_string(fluid_reaction_time_scale) + "."));
               AssertThrow(fluid_reaction_time_scale > 0,
                           ExcMessage("The Fluid reaction time scale for operator splitting must be larger than 0!"));
             }

diff --git a/source/simulator/helper_functions.cc b/source/simulator/helper_functions.cc
@@ -1681,6 +1681,17 @@ namespace aspect
     LinearAlgebra::BlockVector distributed_reaction_vector (introspection.index_sets.system_partitioning,
                                                             mpi_communicator);
 
+    // we use a different (potentially smaller) time step than in the advection scheme.
+    // and for the fixed step scheme, we want all of our reaction time steps (within one advection step) to have the same size
+    const unsigned int number_of_reaction_steps = std::max(static_cast<unsigned int>(time_step / parameters.reaction_time_step),
+                                                           std::max(parameters.reaction_steps_per_advection_step,1U));
+
+    const double reaction_time_step_size = time_step / static_cast<double>(number_of_reaction_steps);
+
+    if (parameters.reaction_solver_type == Parameters<dim>::ReactionSolverType::fixed_step)
+      Assert (reaction_time_step_size > 0,
+              ExcMessage("Reaction time step must be greater than 0."));
+
     pcout << "   Solving composition reactions... " << std::flush;
 
 
@@ -1855,41 +1866,70 @@ namespace aspect
           initial_values_C = in.composition;
           initial_values_T = in.temperature;
 
-          ode.explicit_function = [&] (const double /*time*/,
-                                       const VectorType &y,
-                                       VectorType &ydot)
-          {
-            copy_fields_into_material_model_inputs (y, in);
-            material_model->fill_additional_material_model_inputs(in, solution, fe_values, introspection);
-            material_model->evaluate(in, out);
-            heating_model_manager.evaluate(in, out, heating_model_outputs);
-            copy_rates_into_one_vector (reaction_rate_outputs, heating_model_outputs, ydot);
-          };
-
-          // Make the reaction time steps: We have to update the values of compositional fields and the temperature.
-          // We can reuse the same material model inputs and outputs structure for each reaction time step.
-          // We store the computed updates to temperature and composition in a separate (accumulated_reactions) vector,
-          // so that we can later copy it over to the solution vector.
-          const unsigned int iteration_count = ode.solve_ode(fields);
-
-          total_iteration_count += iteration_count;
-          number_of_solves += 1;
-
-          for (unsigned int j=0; j<n_q_points; ++j)
-            for (unsigned int f=0; f<n_fields; ++f)
+          if (parameters.reaction_solver_type == Parameters<dim>::ReactionSolverType::ARKode)
+            {
+
+              ode.explicit_function = [&] (const double /*time*/,
+                                           const VectorType &y,
+                                           VectorType &ydot)
               {
-                if (f==0)
+                copy_fields_into_material_model_inputs (y, in);
+                material_model->fill_additional_material_model_inputs(in, solution, fe_values, introspection);
+                material_model->evaluate(in, out);
+                heating_model_manager.evaluate(in, out, heating_model_outputs);
+                copy_rates_into_one_vector (reaction_rate_outputs, heating_model_outputs, ydot);
+              };
+
+              // Make the reaction time steps: We have to update the values of compositional fields and the temperature.
+              // We can reuse the same material model inputs and outputs structure for each reaction time step.
+              // We store the computed updates to temperature and composition in a separate (accumulated_reactions) vector,
+              // so that we can later copy it over to the solution vector.
+              const unsigned int iteration_count = ode.solve_ode(fields);
+
+              total_iteration_count += iteration_count;
+              number_of_solves += 1;
+
+              for (unsigned int j=0; j<n_q_points; ++j)
+                for (unsigned int f=0; f<n_fields; ++f)
                   {
-                    in.temperature[j]          = fields[j*n_fields+f];
-                    accumulated_reactions_T[j] = in.temperature[j] - initial_values_T[j];
-                  }
-                else
-                  {
-                    in.composition[j][f-1]          = fields[j*n_fields+f];
-                    accumulated_reactions_C[j][f-1] = in.composition[j][f-1] - initial_values_C[j][f-1];
+                    if (f==0)
+                      {
+                        in.temperature[j]          = fields[j*n_fields+f];
+                        accumulated_reactions_T[j] = in.temperature[j] - initial_values_T[j];
+                      }
+                    else
+                      {
+                        in.composition[j][f-1]          = fields[j*n_fields+f];
+                        accumulated_reactions_C[j][f-1] = in.composition[j][f-1] - initial_values_C[j][f-1];
+                      }
                   }
-              }
+            }
+
+          else if (parameters.reaction_solver_type == Parameters<dim>::ReactionSolverType::fixed_step)
+            {
+              for (unsigned int i=0; i<number_of_reaction_steps; ++i)
+                {
+                  // Loop over composition element
+                  material_model->fill_additional_material_model_inputs(in, solution, fe_values, introspection);
+
+                  material_model->evaluate(in, out);
+                  heating_model_manager.evaluate(in, out, heating_model_outputs);
 
+                  for (unsigned int j=0; j<n_q_points; ++j)
+                    {
+                      for (unsigned int c=0; c<introspection.n_compositional_fields; ++c)
+                        {
+                          // simple forward euler
+                          in.composition[j][c] = in.composition[j][c]
+                                                 + reaction_time_step_size * reaction_rate_outputs->reaction_rates[j][c];
+                          accumulated_reactions_C[j][c] += reaction_time_step_size * reaction_rate_outputs->reaction_rates[j][c];
+                        }
+                      in.temperature[j] = in.temperature[j]
+                                          + reaction_time_step_size * heating_model_outputs.rates_of_temperature_change[j];
+                      accumulated_reactions_T[j] += reaction_time_step_size * heating_model_outputs.rates_of_temperature_change[j];
+                    }
+                }
+            }
 
           cell->get_dof_indices (local_dof_indices);
           const unsigned int component_idx_T = introspection.component_indices.temperature;
@@ -1948,11 +1988,14 @@ namespace aspect
 
     initialize_current_linearization_point();
 
-    const double average_iteration_count = number_of_solves > 0
-                                           ?
-                                           total_iteration_count / number_of_solves
-                                           :
-                                           total_iteration_count;
+    double average_iteration_count = number_of_reaction_steps;
+    if (parameters.reaction_solver_type == Parameters<dim>::ReactionSolverType::ARKode)
+      {
+        if (number_of_solves > 0)
+          average_iteration_count = total_iteration_count / number_of_solves;
+        else
+          average_iteration_count = total_iteration_count;
+      }
 
     pcout << "in "
           << average_iteration_count