Merge pull request #5331 from jdannberg/SUNDIALS_reaction

use SUNDIALS to compute reactions

benchmarks/solubility/solubility.prm

@@ -5,11 +5,14 @@
 # depth and below 60 km depth, and a zero solubility in between. When
 # the upwelling material reaches the boundary at 60 km depth, water is
 # being released and can move relative to the solid as a free fluid
-# phase. Due to its lower density, it moves upwards faster than the
+# phase. Due to its lower density, it moves upwards twice as fast as the
 # solid, until it reaches 30 km depth where it is reabsorbed into the
 # solid. In steady state, the water content should be the same in the
-# top and bottom layer, and in the middle layer it should be lower,
-# proportional to the ratio between solid and melt velocity.
+# top and bottom layer, and in the middle layer it should be zero (at
+# least the water_content field, representing the bound water in the
+# solid). The porosity (free water) should be 0.1, which is half of the
+# bound water content in the upper and lower layer (which follows from
+# mass conservation, since the free fluid moves twice as fast).
 
 set Additional shared libraries            = ./plugin/libsolubility.so
 set Adiabatic surface temperature          = 1600
@@ -27,18 +30,6 @@ set End time                               = 6e6
 # the operator splitting scheme.
 set Use operator splitting                     = true
 
-subsection Solver parameters
-  subsection Operator splitting parameters
-    # We choose the size of the reaction time step as 200 years, small enough
-    # so that it can accurately model melting and freezing.
-    set Reaction time step                     = 1e4
-
-    # Additionally, we always want to do at least 10 operator splitting time
-    # steps in each model time step, to accurately compute the reactions.
-    set Reaction time steps per advection step = 10
-  end
-end
-
 # There are two compositional fields, one that tracks the amount of free water
 # (the porosity) and one that tracks the amount of bound water (the water_content).
 subsection Compositional fields

cookbooks/mid_ocean_ridge/mid_ocean_ridge.prm

@@ -24,16 +24,6 @@ set End time                               = 8e6
 set Use operator splitting                     = true
 
 subsection Solver parameters
-  subsection Operator splitting parameters
-    # We choose the size of the reaction time step as 200 years, small enough
-    # so that it can accurately model melting and freezing.
-    set Reaction time step                     = 2e2
-
-    # Additionally, we always want to do at least 10 operator splitting time
-    # steps in each model time step, to accurately compute the reactions.
-    set Reaction time steps per advection step = 10
-  end
-
   # Because this model includes strong localized viscosity contrasts we
   # increase the robustness of the solver at the cost of memory consumption.
   subsection Stokes solver parameters

cookbooks/tian_parameterization_kinematic_slab/coupled-two-phase-tian-parameterization-kinematic-slab.prm

@@ -55,13 +55,7 @@ subsection Discretization
   end
 end
 
-# Do 10 reaction time steps per advection time step.
 subsection Solver parameters
-  subsection Operator splitting parameters
-    set Reaction time step                     = 100
-    set Reaction time steps per advection step = 10
-  end
-
   # Set a stricter linear solver tolerance to ensure convergence
   subsection Stokes solver parameters
     set Linear solver tolerance             = 1e-10

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;
 

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)

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"),

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) + "."));
+          }
       }
     }
   }

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!"));
             }