[WIP] Enable operator split reactions at end of the timestep as well

benchmarks/free_surface_tractions/viscoelastic/free_surface_VE_cylinder_2D_loading.prm

@@ -47,6 +47,7 @@ subsection Solver parameters
     set Reaction solver type                   = fixed step
     set Reaction time step                     = 5000 # larger than maximum Stokes time step
     set Reaction time steps per advection step = 1
+    set Reaction solve strategy                = after nonlinear solver
   end
 end
 

benchmarks/free_surface_tractions/viscoelastic/free_surface_VE_cylinder_3D_loading.prm

@@ -48,6 +48,7 @@ subsection Solver parameters
     set Reaction solver type                   = fixed step
     set Reaction time step                     = 5000 # larger than maximum Stokes time step
     set Reaction time steps per advection step = 1
+    set Reaction solve strategy                = after nonlinear solver
   end
 end
 

benchmarks/infill_density/infill_ascii.prm

@@ -38,6 +38,7 @@ subsection Solver parameters
     set Reaction solver type                   = fixed step
     set Reaction time step                     = 5000 # larger than maximum Stokes time step
     set Reaction time steps per advection step = 1
+    set Reaction solve strategy                = after nonlinear solver
   end
 end
 

benchmarks/viscoelastic_bending_beam/viscoelastic_bending_beam.prm

@@ -49,6 +49,7 @@ subsection Solver parameters
     set Reaction solver type                   = fixed step
     set Reaction time step                     = 5000 # larger than maximum Stokes time step
     set Reaction time steps per advection step = 1
+    set Reaction solve strategy                = after nonlinear solver
   end
 end
 

benchmarks/viscoelastic_plastic_simple_shear/viscoelastic_plastic_simple_shear.prm

@@ -44,6 +44,7 @@ subsection Solver parameters
     set Reaction solver type                   = fixed step
     set Reaction time step                     = 5000 # larger than maximum Stokes time step
     set Reaction time steps per advection step = 1
+    set Reaction solve strategy                = after nonlinear solver
   end
 end
 

benchmarks/viscoelastic_plate_flexure/viscoelastic_plate_flexure.prm

@@ -40,6 +40,7 @@ subsection Solver parameters
     set Reaction solver type                   = fixed step
     set Reaction time step                     = 5000 # larger than maximum Stokes time step
     set Reaction time steps per advection step = 1
+    set Reaction solve strategy                = after nonlinear solver
   end
 end
 

benchmarks/viscoelastic_relaxation/viscoelastic_relaxation.prm

@@ -28,6 +28,7 @@ subsection Solver parameters
     set Reaction solver type                   = fixed step
     set Reaction time step                     = 5000 # larger than maximum Stokes time step
     set Reaction time steps per advection step = 1
+    set Reaction solve strategy                = after nonlinear solver
   end
 end
 

benchmarks/viscoelastic_sheared_torsion/viscoelastic_sheared_torsion.prm

@@ -62,6 +62,7 @@ subsection Solver parameters
     set Reaction solver type                   = fixed step
     set Reaction time step                     = 5000 # larger than maximum Stokes time step
     set Reaction time steps per advection step = 1
+    set Reaction solve strategy                = after nonlinear solver
   end
 end
 

benchmarks/viscoelastic_stress_build_up/viscoelastic_stress_build_up.prm

@@ -56,6 +56,7 @@ subsection Solver parameters
     set Reaction solver type                   = fixed step
     set Reaction time step                     = 5000 # larger than maximum Stokes time step
     set Reaction time steps per advection step = 1
+    set Reaction solve strategy                = after nonlinear solver
   end
 end
 

doc/sphinx/parameters/Solver_20parameters.md

@@ -199,6 +199,15 @@ Once derivatives are used in a Newton method, ASPECT always uses the Eisenstat W
 
 (parameters:Solver_20parameters/Operator_20splitting_20parameters)=
 ## **Subsection:** Solver parameters / Operator splitting parameters
+::::{dropdown} __Parameter:__ {ref}`Reaction solve strategy<parameters:Solver_20parameters/Operator_20splitting_20parameters/Reaction_20solve_20strategy>`
+:name: parameters:Solver_20parameters/Operator_20splitting_20parameters/Reaction_20solve_20strategy
+**Default value:** before nonlinear solver
+
+**Pattern:** [Selection before nonlinear solver|after nonlinear solver ]
+
+**Documentation:** Whether the reaction solve should be done before the nonlinear solver (default) or after the nonlinear solver. The latter strategy should be used when &ldquo;Enable elasticity&rdquo; is set to true. The strategy parameter is only considered if the parameter &ldquo;Use operator splitting&rdquo; is set to true. Units: none.
+::::
+
 ::::{dropdown} __Parameter:__ {ref}`Reaction solver relative tolerance<parameters:Solver_20parameters/Operator_20splitting_20parameters/Reaction_20solver_20relative_20tolerance>`
 :name: parameters:Solver_20parameters/Operator_20splitting_20parameters/Reaction_20solver_20relative_20tolerance
 **Default value:** 1e-6
@@ -214,7 +223,7 @@ Once derivatives are used in a Newton method, ASPECT always uses the Eisenstat W
 
 **Pattern:** [Selection ARKode|fixed step ]
 
-**Documentation:** This parameter determines what solver will be used when the reactions are computed within the operator splitting scheme. For reactions where the reaction rate is a known, finite quantity, the appropriate choice is &lsquo;ARKode&rsquo;, which uses an ODE solver from SUNDIALs ARKode (adaptive-step additive Runge Kutta ODE solver methods) to compute the solution. ARKode will pick a reasonable step size based on the reaction rate and the given &lsquo;Reaction solver relative tolerance&rsquo;. However, in some cases we have instantaneous reactions, where we know the new value of a compositional field (and the reaction rate would be infinite), or reaction where we need to know or be able to control the step size we use to compute the reactions. In theses cases, it is appropriate to use the &lsquo;fixed step&rsquo; scheme, a method that a forward Euler scheme and a fixed number of steps given by the &lsquo;Reaction time step&rsquo; and &lsquo;Reaction time steps per advection step&rsquo; parameters.
+**Documentation:** This parameter determines what solver will be used when the reactions are computed within the operator splitting scheme. For reactions where the reaction rate is a known, finite quantity, the appropriate choice is &lsquo;ARKode&rsquo;, which uses an ODE solver from SUNDIALs ARKode (adaptive-step additive Runge Kutta ODE solver methods) to compute the solution. ARKode will pick a reasonable step size based on the reaction rate and the given &lsquo;Reaction solver relative tolerance&rsquo;. However, in some cases we have instantaneous reactions, where we know the new value of a compositional field (and the reaction rate would be infinite), or reactions where we need to know or be able to control the step size we use to compute the reactions. In theses cases, it is appropriate to use the &lsquo;fixed step&rsquo; scheme, a method that uses a forward Euler scheme and a fixed number of steps given by the &lsquo;Reaction time step&rsquo; and &lsquo;Reaction time steps per advection step&rsquo; parameters.
 ::::
 
 ::::{dropdown} __Parameter:__ {ref}`Reaction time step<parameters:Solver_20parameters/Operator_20splitting_20parameters/Reaction_20time_20step>`

include/aspect/parameters.h

@@ -474,6 +474,33 @@ namespace aspect
       }
     };
 
+    /**
+     * This enum represents the different choices for when to
+     * apply operator splitting. See @p reaction_strategy.
+     */
+    struct ReactionStrategy
+    {
+      enum Kind
+      {
+        before_nonlinear_solver,
+        after_nonlinear_solver
+      };
+
+      static
+      Kind
+      parse(const std::string &input)
+      {
+        if (input == "before nonlinear solver")
+          return before_nonlinear_solver;
+        else if (input == "after nonlinear solver")
+          return after_nonlinear_solver;
+        else
+          AssertThrow(false, ExcNotImplemented());
+
+        return Kind();
+      }
+    };
+
     /**
      * Use the struct aspect::CompositionalFieldDescription
      */
@@ -614,6 +641,8 @@ namespace aspect
     double                         ARKode_relative_tolerance;
     double                         reaction_time_step;
     unsigned int                   reaction_steps_per_advection_step;
+    bool                           apply_reaction_solve_at_end_of_time_step;
+    typename ReactionStrategy::Kind reaction_strategy;
 
     // subsection: Diffusion solver parameters
     double                         diffusion_length_scale;

include/aspect/particle/property/elastic_stress.h

@@ -51,10 +51,9 @@ namespace aspect
           void initialize () override;
 
           /**
-           * Function to update particles after they have been restored
-           * to their position and values from the beginning of the timestep.
-           * This restoring happens at the beginning of nonlinear iterations
-           * of iterative Advection solver schemes.
+           * Function to update the stresses on the particles after the
+           * nonlinear solver, similar to the operator split for stresses
+           * carried on fields.
            */
           void
           update_particles (typename Particle::Manager<dim> &particle_manager) const;

include/aspect/simulator_signals.h

@@ -331,6 +331,16 @@ namespace aspect
      */
     boost::signals2::signal<void (Particle::Manager<dim> &)> post_restore_particles;
 
+
+    /**
+     * A signal that is triggered after particles have been restored to their position
+     * and property values from the beginning of the current timestep. This happens
+     * at the beginning of each nonlinear iteration (except for the first iteration of the timestep)
+     * of iterative advection schemes.
+     * Parameters are a reference to a ParticleManager.
+     */
+    boost::signals2::signal<void (Particle::Manager<dim> &)> post_nonlinear_solver_loop;
+
   };