From 2873bf9036f45399305928f677827ced75b6fed8 Mon Sep 17 00:00:00 2001 From: anne-glerum Date: Fri, 20 Mar 2026 14:08:19 +0100 Subject: [PATCH 01/15] Fix typo --- source/simulator/parameters.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/source/simulator/parameters.cc b/source/simulator/parameters.cc index 4b948299999..d6853cc9754 100644 --- a/source/simulator/parameters.cc +++ b/source/simulator/parameters.cc @@ -639,9 +639,9 @@ namespace aspect "`Reaction solver relative tolerance'. " "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 " + "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 `fixed step' scheme, a method that a forward Euler scheme and a " + "to use the `fixed step' scheme, a method that uses a forward Euler scheme and a " "fixed number of steps given by the `Reaction time step' and " "`Reaction time steps per advection step' parameters. "); From 25e82dbb892185d65ad249433b9fbef1f5297367 Mon Sep 17 00:00:00 2001 From: anne-glerum Date: Fri, 20 Mar 2026 14:15:41 +0100 Subject: [PATCH 02/15] Add bool for reaction solve at end of time step --- include/aspect/parameters.h | 1 + source/simulator/parameters.cc | 8 ++++++++ 2 files changed, 9 insertions(+) diff --git a/include/aspect/parameters.h b/include/aspect/parameters.h index 02308c545b1..064a1c6bbd6 100644 --- a/include/aspect/parameters.h +++ b/include/aspect/parameters.h @@ -614,6 +614,7 @@ 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; // subsection: Diffusion solver parameters double diffusion_length_scale; diff --git a/source/simulator/parameters.cc b/source/simulator/parameters.cc index d6853cc9754..9bb1340debe 100644 --- a/source/simulator/parameters.cc +++ b/source/simulator/parameters.cc @@ -678,6 +678,14 @@ namespace aspect "this criterion and the ``Reaction time step'', whichever yields the " "smaller time step. " "Units: none."); + + prm.declare_entry ("Apply reaction solve at end of time step", "false", + Patterns::Bool (), + "Whether the reaction solve should be done at the beginning of a time step " + "(if set to false, default) or at the end of a time step (if true). " + "This is only used if the parameter ``Use operator splitting'' is set " + "to true. " + "Units: none."); } prm.leave_subsection (); prm.enter_subsection ("Diffusion solver parameters"); From 980ad9dfcc07471e5ac1a7d55616faf1c33c6e82 Mon Sep 17 00:00:00 2001 From: anne-glerum Date: Fri, 20 Mar 2026 14:32:19 +0100 Subject: [PATCH 03/15] Add compute_reaction call at end of timestep --- source/simulator/core.cc | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/source/simulator/core.cc b/source/simulator/core.cc index cc19177aa07..43709e680b1 100644 --- a/source/simulator/core.cc +++ b/source/simulator/core.cc @@ -1973,7 +1973,8 @@ namespace aspect } // Compute the reactions of compositional fields and temperature in case of operator splitting. - if (parameters.use_operator_splitting) + // Only if specifically requested, do we compute the reactions at the end of the time step instead. + if (parameters.use_operator_splitting && !parameters.apply_reaction_solve_at_end_of_time_step) compute_reactions (); try @@ -2263,6 +2264,13 @@ namespace aspect const double new_time_step_size = time_stepping_manager.get_next_time_step_size(); + // Compute the reactions of compositional fields and temperature here after + // the time step solve loop instead of at the beginning of the next time step + // if requested by the user. + if (parameters.use_operator_splitting && + parameters.apply_reaction_solve_at_end_of_time_step) + compute_reactions (); + // if we postprocess nonlinear iterations, this function is called within // solve_timestep () in the individual solver schemes if (!time_stepping_manager.should_repeat_time_step() From d96c5ccfd481517df323b44d4ed09f3d51e87842 Mon Sep 17 00:00:00 2001 From: anne-glerum Date: Fri, 20 Mar 2026 14:56:22 +0100 Subject: [PATCH 04/15] Assert reactions at end of time step for elasticity --- source/material_model/rheology/elasticity.cc | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/source/material_model/rheology/elasticity.cc b/source/material_model/rheology/elasticity.cc index 164dd76ddf6..0b23cf36c5c 100644 --- a/source/material_model/rheology/elasticity.cc +++ b/source/material_model/rheology/elasticity.cc @@ -200,6 +200,11 @@ namespace aspect AssertThrow(!this->get_parameters().use_operator_splitting, ExcMessage("If stresses are tracked on particles, the stress update is applied by the particle property 'elastic stress' " "and operator splitting should not be turned on. ")); + // The operator split should happen at the end of the timestep, as it add the dynamic stress update from the current time step. + if (this->get_parameters().use_operator_splitting) + AssertThrow(this->get_parameters().apply_reaction_solve_at_end_of_time_step == true, + ExcMessage("The operator splitting scheme used to update the stresses should be applied at the end of each " + "time step by setting 'Apply reaction solve at end of time step' to true.")); // Check that 3+3 in 2D or 6+6 in 3D stress fields exist. AssertThrow((this->introspection().get_number_of_fields_of_type(CompositionalFieldDescription::stress) == 2*SymmetricTensor<2,dim>::n_independent_components), From a4a1461f8d2748fd50a6c85092581cc8c2877173 Mon Sep 17 00:00:00 2001 From: anne-glerum Date: Fri, 20 Mar 2026 15:25:38 +0100 Subject: [PATCH 05/15] Update comments elasticity reaction_rates --- source/material_model/rheology/elasticity.cc | 47 +++++++++----------- 1 file changed, 22 insertions(+), 25 deletions(-) diff --git a/source/material_model/rheology/elasticity.cc b/source/material_model/rheology/elasticity.cc index 0b23cf36c5c..b95eed083ee 100644 --- a/source/material_model/rheology/elasticity.cc +++ b/source/material_model/rheology/elasticity.cc @@ -602,7 +602,7 @@ namespace aspect } // The following function computes the reaction rates for the operator - // splitting step that at the beginning of the new timestep $t+dtc$ updates the + // splitting step that at the end of the timestep $t$ updates the // stored compositions $tau^{0\mathrm{adv}}$ at time $t$ to $tau^{t}$. // This update consists of the stress change resulting from system evolution, // but does not advect or rotate the stress tensor. Advection is done by @@ -633,11 +633,10 @@ namespace aspect if (this->get_timestep_number() == 0) return; - // At the moment when the reaction rates are required (at the beginning of the timestep), - // the solution vector 'solution' holds the stress from the previous timestep, - // advected into the new position of the previous timestep, so $\tau^{t}_{0adv}$. - // This is the same as the vector 'old_solution' holds. At later moments during the current timestep, - // 'solution' will hold the current_linearization_point instead of the solution of the previous timestep. + // At the moment when the reaction rates are required (at the end of the timestep), + // the solution vector 'solution' holds the stress from the previous timestep + // advected and rotated into the new position of the current timestep, so $\tau^{t}_{0adv}$. + // The vector 'old_solution' holds the full deviatoric stress solution of the previous timestep. // // In case fields are used to track the stresses, MaterialModelInputs are based on 'solution' // when calling the MaterialModel for the reaction rates. When particles are used, MaterialModelInputs @@ -686,18 +685,19 @@ namespace aspect for (unsigned int i = 0; i < in.n_evaluation_points(); ++i) { - // Get $\tau^{0adv}$ of the previous timestep t from the compositional fields. - // This stress includes the rotation and advection of the previous timestep, + // Get $\tau^{0adv}$ from the compositional fields. + // This stress includes the rotation and advection of the current timestep, // i.e., the reaction term (which prescribes the change in stress due to rotation - // over the previous timestep) has already been applied during the previous timestep. + // over the timestep) has already been applied during the current timestep. const SymmetricTensor<2, dim> stress_0_t (Utilities::Tensors::to_symmetric_tensor(&in.composition[i][stress_start_index], &in.composition[i][stress_start_index]+n_independent_components)); - // Get the old stress that is used to interpolate to timestep $t+\Delta t_c$. It is stored on the + // Get the old stress that is used to interpolate to the current computational time step + // if it differs from the elastic time step. It is stored on the // second set of n_independent_components fields, e.g. in 2D on field 3, 4 and 5. - // The old stress was advected into the previous timestep, but not rotated. - // Below we update it to full stress of the previous timestep, so that it can be - // advected into the current timestep. + // The old stress was advected into the current timestep, but not rotated. + // Below we update it to full stress of the current timestep after using it to + // compute the full stress. const SymmetricTensor<2, dim> stress_old (Utilities::Tensors::to_symmetric_tensor(&in.composition[i][stress_start_index+n_independent_components], &in.composition[i][stress_start_index+n_independent_components]+n_independent_components)); @@ -720,16 +720,11 @@ namespace aspect + (1. - timestep_ratio) * (1. - effective_creep_viscosity / elastic_viscosity) * stress_old; // Fill reaction rates. - // During this timestep, the reaction rates will be multiplied + // The reaction rates will be multiplied // with the current timestep size to turn the rate of change into a change. - // However, this update belongs - // to the previous timestep. Therefore we divide by the - // current timestep and multiply with the previous one. - // When multiplied with the current timestep, this will give - // (rate * previous_dt / current_dt) * current_dt = rate * previous_dt = previous_change. - // previous_change = stress_t - stress_0_t. - // To compute the rate we should return to the operator splitting scheme, - // we therefore divide the change in stress by the current timestep current_dt (= dtc). + // The change in stress is the difference between the fully updated stress + // and the advected and rotated stress. Dividing by the current time step + // size (= dtc) gives the rate of change. const double dtc = timestep_ratio * elastic_timestep(); const SymmetricTensor<2, dim> stress_update = (stress_t - stress_0_t) / dtc; @@ -739,9 +734,11 @@ namespace aspect &reaction_rate_out->reaction_rates[i][stress_start_index]+n_independent_components); // Also update the second set of stresses, stress_old, with the newly computed stress, - // which in the rest of the timestep will serve as the old stress advected but not rotated - // into the current timestep. This function fill_reaction_rates is only called at the - // beginning of the timestep, and so this update only happens once. + // which in the next timestep will serve as the old stress advected but not rotated + // into the current timestep. For fields, this function fill_reaction_rates is only + // called at the end of the timestep, and so this update only happens once. + // For particles, this function is called each time they have been restored to their + // original position from the beginning of the time step. const SymmetricTensor<2, dim> stress_old_update = (stress_t - stress_old) / dtc; Utilities::Tensors::unroll_symmetric_tensor_into_array(stress_old_update, From 463cd109f4da0a2674a00b3b65d06be878763eb3 Mon Sep 17 00:00:00 2001 From: anne-glerum Date: Fri, 20 Mar 2026 16:04:10 +0100 Subject: [PATCH 06/15] Actually read in parameter --- source/simulator/parameters.cc | 1 + 1 file changed, 1 insertion(+) diff --git a/source/simulator/parameters.cc b/source/simulator/parameters.cc index 9bb1340debe..bb6720a9446 100644 --- a/source/simulator/parameters.cc +++ b/source/simulator/parameters.cc @@ -1672,6 +1672,7 @@ namespace aspect if (convert_to_years == true) reaction_time_step *= year_in_seconds; reaction_steps_per_advection_step = prm.get_integer ("Reaction time steps per advection step"); + apply_reaction_solve_at_end_of_time_step = prm.get_bool ("Apply reaction solve at end of time step"); } prm.leave_subsection (); prm.enter_subsection ("Diffusion solver parameters"); From d11e40d34911d5495dda6ebf0125edb4bf345f79 Mon Sep 17 00:00:00 2001 From: anne-glerum Date: Fri, 27 Mar 2026 14:55:59 +0100 Subject: [PATCH 07/15] Fix typo --- source/material_model/rheology/elasticity.cc | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/source/material_model/rheology/elasticity.cc b/source/material_model/rheology/elasticity.cc index b95eed083ee..5708c32ede7 100644 --- a/source/material_model/rheology/elasticity.cc +++ b/source/material_model/rheology/elasticity.cc @@ -200,7 +200,8 @@ namespace aspect AssertThrow(!this->get_parameters().use_operator_splitting, ExcMessage("If stresses are tracked on particles, the stress update is applied by the particle property 'elastic stress' " "and operator splitting should not be turned on. ")); - // The operator split should happen at the end of the timestep, as it add the dynamic stress update from the current time step. + // The operator split should happen at the end of the timestep, + // as it adds the dynamic stress update from the current time step. if (this->get_parameters().use_operator_splitting) AssertThrow(this->get_parameters().apply_reaction_solve_at_end_of_time_step == true, ExcMessage("The operator splitting scheme used to update the stresses should be applied at the end of each " From 74a36eb2ae34e30ac1df8031a425991d44a08ca9 Mon Sep 17 00:00:00 2001 From: anne-glerum Date: Fri, 27 Mar 2026 15:32:09 +0100 Subject: [PATCH 08/15] Use struct for reaction solve strategy and move call to after nonlinear solver --- include/aspect/parameters.h | 28 ++++++++++++++++++++++++++++ source/simulator/core.cc | 18 ++++++++---------- source/simulator/parameters.cc | 15 ++++++++------- 3 files changed, 44 insertions(+), 17 deletions(-) diff --git a/include/aspect/parameters.h b/include/aspect/parameters.h index 064a1c6bbd6..4b9f88edc6d 100644 --- a/include/aspect/parameters.h +++ b/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 */ @@ -615,6 +642,7 @@ namespace aspect 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; diff --git a/source/simulator/core.cc b/source/simulator/core.cc index 43709e680b1..6faf4e5ae44 100644 --- a/source/simulator/core.cc +++ b/source/simulator/core.cc @@ -1972,9 +1972,9 @@ namespace aspect signals.post_mesh_deformation(*this); } - // Compute the reactions of compositional fields and temperature in case of operator splitting. - // Only if specifically requested, do we compute the reactions at the end of the time step instead. - if (parameters.use_operator_splitting && !parameters.apply_reaction_solve_at_end_of_time_step) + // Compute the reactions of compositional fields and temperature in case of operator splitting + // before the nonlinear solver instead of after. + if (parameters.use_operator_splitting && parameters.reaction_strategy == Parameters::ReactionStrategy::before_nonlinear_solver) compute_reactions (); try @@ -2117,6 +2117,11 @@ namespace aspect AssertThrow(false, ExcNotImplemented()); } } + + // Compute the reactions of compositional fields and temperature here after + // the nonlinear solver instead of before. + if (parameters.use_operator_splitting && parameters.reaction_strategy == Parameters::ReactionStrategy::after_nonlinear_solver) + compute_reactions (); } @@ -2264,13 +2269,6 @@ namespace aspect const double new_time_step_size = time_stepping_manager.get_next_time_step_size(); - // Compute the reactions of compositional fields and temperature here after - // the time step solve loop instead of at the beginning of the next time step - // if requested by the user. - if (parameters.use_operator_splitting && - parameters.apply_reaction_solve_at_end_of_time_step) - compute_reactions (); - // if we postprocess nonlinear iterations, this function is called within // solve_timestep () in the individual solver schemes if (!time_stepping_manager.should_repeat_time_step() diff --git a/source/simulator/parameters.cc b/source/simulator/parameters.cc index bb6720a9446..f3134327de3 100644 --- a/source/simulator/parameters.cc +++ b/source/simulator/parameters.cc @@ -679,12 +679,13 @@ namespace aspect "smaller time step. " "Units: none."); - prm.declare_entry ("Apply reaction solve at end of time step", "false", - Patterns::Bool (), - "Whether the reaction solve should be done at the beginning of a time step " - "(if set to false, default) or at the end of a time step (if true). " - "This is only used if the parameter ``Use operator splitting'' is set " - "to true. " + prm.declare_entry ("Reaction solve strategy", "before nonlinear solver", + Patterns::Selection("before nonlinear solver|after nonlinear solver"), + "Whether the reaction solve should be done before the nonlinear solver " + "(default) or after the nonlinear solver. The latter strategy should be " + "used when ``Enable elasticity'' is set to true. " + "The strategy parameter is only considered if the parameter " + "``Use operator splitting'' is set to true. " "Units: none."); } prm.leave_subsection (); @@ -1672,7 +1673,7 @@ namespace aspect if (convert_to_years == true) reaction_time_step *= year_in_seconds; reaction_steps_per_advection_step = prm.get_integer ("Reaction time steps per advection step"); - apply_reaction_solve_at_end_of_time_step = prm.get_bool ("Apply reaction solve at end of time step"); + reaction_strategy = ReactionStrategy::parse(prm.get("Reaction solve strategy")); } prm.leave_subsection (); prm.enter_subsection ("Diffusion solver parameters"); From 5eb6c50600cb45c79ebae575f96387c2bf7d8227 Mon Sep 17 00:00:00 2001 From: anne-glerum Date: Fri, 27 Mar 2026 16:33:50 +0100 Subject: [PATCH 09/15] Use new parameter for AssertThrow --- source/material_model/rheology/elasticity.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/source/material_model/rheology/elasticity.cc b/source/material_model/rheology/elasticity.cc index 5708c32ede7..7bc10fc1f77 100644 --- a/source/material_model/rheology/elasticity.cc +++ b/source/material_model/rheology/elasticity.cc @@ -203,9 +203,9 @@ namespace aspect // The operator split should happen at the end of the timestep, // as it adds the dynamic stress update from the current time step. if (this->get_parameters().use_operator_splitting) - AssertThrow(this->get_parameters().apply_reaction_solve_at_end_of_time_step == true, + AssertThrow(this->get_parameters().reaction_strategy == Parameters::ReactionStrategy::after_nonlinear_solver, ExcMessage("The operator splitting scheme used to update the stresses should be applied at the end of each " - "time step by setting 'Apply reaction solve at end of time step' to true.")); + "time step by setting 'Reaction solve strategy' to 'after nonlinear solver'.")); // Check that 3+3 in 2D or 6+6 in 3D stress fields exist. AssertThrow((this->introspection().get_number_of_fields_of_type(CompositionalFieldDescription::stress) == 2*SymmetricTensor<2,dim>::n_independent_components), From 4b200d17eb80f13fb61be75b4ae208f97099ab36 Mon Sep 17 00:00:00 2001 From: anne-glerum Date: Fri, 27 Mar 2026 16:35:00 +0100 Subject: [PATCH 10/15] Use another signal to update the particle stresses with the current solution --- .../aspect/particle/property/elastic_stress.h | 7 ++-- include/aspect/simulator_signals.h | 10 +++++ source/particle/property/elastic_stress.cc | 41 +++++++++---------- source/simulator/core.cc | 10 +++++ source/simulator/solver_schemes.cc | 2 - 5 files changed, 43 insertions(+), 27 deletions(-) diff --git a/include/aspect/particle/property/elastic_stress.h b/include/aspect/particle/property/elastic_stress.h index b88cb969532..f37dfda5373 100644 --- a/include/aspect/particle/property/elastic_stress.h +++ b/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 &particle_manager) const; diff --git a/include/aspect/simulator_signals.h b/include/aspect/simulator_signals.h index 6634947420f..8db0611dc8b 100644 --- a/include/aspect/simulator_signals.h +++ b/include/aspect/simulator_signals.h @@ -331,6 +331,16 @@ namespace aspect */ boost::signals2::signal &)> 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 &)> post_nonlinear_solver_loop; + }; diff --git a/source/particle/property/elastic_stress.cc b/source/particle/property/elastic_stress.cc index 62a228b1171..84bd428c827 100644 --- a/source/particle/property/elastic_stress.cc +++ b/source/particle/property/elastic_stress.cc @@ -75,9 +75,9 @@ namespace aspect stress_field_indices.begin(), stress_field_indices.end(), std::inserter(non_stress_field_indices, non_stress_field_indices.begin())); - // Connect to the signal after particles are restored at the beginning of - // a nonlinear iteration of iterative advection schemes. - this->get_signals().post_restore_particles.connect( + // Connect to the signal after the nonlinear solver to apply + // the stress update of the current time step. + this->get_signals().post_nonlinear_solver_loop.connect( [&](typename Particle::Manager &particle_manager) { this->update_particles(particle_manager); @@ -86,7 +86,6 @@ namespace aspect } - template void ElasticStress::update_particles(typename Particle::Manager &particle_manager) const @@ -167,11 +166,11 @@ namespace aspect const std::shared_ptr> reaction_rate_outputs = material_outputs_cell.template get_additional_output_object>(); - // Collect the values of the old solution restricted to the current cell's DOFs - small_vector old_solution_values(this->get_fe().dofs_per_cell); - cell->get_dof_values(this->get_old_solution(), - old_solution_values.begin(), - old_solution_values.end()); + // Collect the values of the solution restricted to the current cell's DOFs + small_vector solution_values(this->get_fe().dofs_per_cell); + cell->get_dof_values(this->get_solution(), + solution_values.begin(), + solution_values.end()); EvaluationFlags::EvaluationFlags evaluation_flags_union = EvaluationFlags::nothing; for (const auto flag : evaluation_flags) @@ -183,15 +182,15 @@ namespace aspect // Reinitialize and evaluate the requested solution values and gradients evaluators->reinit(cell, {positions.data(), positions.size()}); - evaluators->evaluate({old_solution_values.data(),old_solution_values.size()}, + evaluators->evaluate({solution_values.data(),solution_values.size()}, evaluation_flags); } // To store the old solutions - std::vector> old_solution(n_particles_in_cell,small_vector(evaluators->n_components(), numbers::signaling_nan())); + std::vector> solution(n_particles_in_cell,small_vector(evaluators->n_components(), numbers::signaling_nan())); // To store the old gradients - std::vector,50>> old_gradients(n_particles_in_cell,small_vector,50>(evaluators->n_components(), numbers::signaling_nan>())); + std::vector,50>> gradients(n_particles_in_cell,small_vector,50>(evaluators->n_components(), numbers::signaling_nan>())); // Loop over all particles in the cell auto particle = particles_in_cell.begin(); @@ -199,26 +198,26 @@ namespace aspect { // Evaluate the old solution, but only if it is requested in the update_flags if (evaluation_flags_union & EvaluationFlags::values) - evaluators->get_solution(i, {&old_solution[i][0],old_solution[i].size()}, evaluation_flags); + evaluators->get_solution(i, {&solution[i][0],solution[i].size()}, evaluation_flags); // Evaluate the old gradients, but only if they are requested in the update_flags if (evaluation_flags_union & EvaluationFlags::gradients) - evaluators->get_gradients(i, {&old_gradients[i][0],old_gradients[i].size()}, evaluation_flags); + evaluators->get_gradients(i, {&gradients[i][0],gradients[i].size()}, evaluation_flags); // Fill material model input // Get the real location of the particle material_inputs_cell.position[i] = particle->get_location(); - material_inputs_cell.temperature[i] = old_solution[i][this->introspection().component_indices.temperature]; + material_inputs_cell.temperature[i] = solution[i][this->introspection().component_indices.temperature]; - material_inputs_cell.pressure[i] = old_solution[i][this->introspection().component_indices.pressure]; + material_inputs_cell.pressure[i] = solution[i][this->introspection().component_indices.pressure]; for (unsigned int d = 0; d < dim; ++d) - material_inputs_cell.velocity[i][d] = old_solution[i][this->introspection().component_indices.velocities[d]]; + material_inputs_cell.velocity[i][d] = solution[i][this->introspection().component_indices.velocities[d]]; - // Fill the non-stress composition inputs with the old_solution. + // Fill the non-stress composition inputs with the current solution. for (const unsigned int &n : non_stress_field_indices) - material_inputs_cell.composition[i][n] = old_solution[i][this->introspection().component_indices.compositional_fields[n]]; + material_inputs_cell.composition[i][n] = solution[i][this->introspection().component_indices.compositional_fields[n]]; // Retrieve the ve_stress_* values from the particles and fields and // fill the material model inputs with a weighted average of the two. @@ -227,14 +226,14 @@ namespace aspect for (unsigned int n = 0; n < 2*SymmetricTensor<2,dim>::n_independent_components; ++n) { const double particle_stress_value = particle->get_properties()[data_position + n]; - const double field_stress_value = old_solution[i][this->introspection().component_indices.compositional_fields[stress_field_indices[n]]]; + const double field_stress_value = solution[i][this->introspection().component_indices.compositional_fields[stress_field_indices[n]]]; const double stress_value = particle_weight * particle_stress_value + (1-particle_weight) * field_stress_value; material_inputs_cell.composition[i][stress_field_indices[n]] = stress_value; } Tensor<2,dim> grad_u; for (unsigned int d=0; d::ReactionStrategy::after_nonlinear_solver) diff --git a/source/simulator/solver_schemes.cc b/source/simulator/solver_schemes.cc index e86aa22216d..469aaeb4605 100644 --- a/source/simulator/solver_schemes.cc +++ b/source/simulator/solver_schemes.cc @@ -238,8 +238,6 @@ namespace aspect // Signal that particles have been restored (or we just started a time // step, which means they are in the same state as post restore). - // This signal can be used to apply operator - // splitting on the particle properties. signals.post_restore_particles(particle_manager); // Advect particles, but not if we are in the initial refinement stage From 813b13200fb38d0908fc3617569f15cb2d66956c Mon Sep 17 00:00:00 2001 From: anne-glerum Date: Fri, 27 Mar 2026 17:31:30 +0100 Subject: [PATCH 11/15] Update elasticity comments --- source/material_model/rheology/elasticity.cc | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/source/material_model/rheology/elasticity.cc b/source/material_model/rheology/elasticity.cc index 7bc10fc1f77..fc0fec1ee7f 100644 --- a/source/material_model/rheology/elasticity.cc +++ b/source/material_model/rheology/elasticity.cc @@ -608,7 +608,9 @@ namespace aspect // This update consists of the stress change resulting from system evolution, // but does not advect or rotate the stress tensor. Advection is done by // solving the advection equation and the stress tensor is rotated through - // the source term (reaction_terms) of that same equation. + // the source term (reaction_terms) of that same equation. Or, in case stresses + // are tracked on particles, by advection of the particles and the updating + // of the their properties with the reaction_terms. template void Elasticity::fill_reaction_rates (const MaterialModel::MaterialModelInputs &in, @@ -634,20 +636,16 @@ namespace aspect if (this->get_timestep_number() == 0) return; - // At the moment when the reaction rates are required (at the end of the timestep), + // At the moment when the reaction rates are required (after the nonlinear solver loop), // the solution vector 'solution' holds the stress from the previous timestep // advected and rotated into the new position of the current timestep, so $\tau^{t}_{0adv}$. - // The vector 'old_solution' holds the full deviatoric stress solution of the previous timestep. // // In case fields are used to track the stresses, MaterialModelInputs are based on 'solution' // when calling the MaterialModel for the reaction rates. When particles are used, MaterialModelInputs - // for this function are filled with the old_solution (including for the strain rate), except for the + // for this function are filled with the solution (including for the strain rate), except for the // compositions that represent the stress tensor components, these are taken directly from the // particles in the property plugin by default (although this can be changed from the input file). - // As the particles are restored to their pre-advection location at the beginning of each nonlinear iteration, - // their values and positions correspond to the old solution. - // This means that in both cases we can use 'in' to get to the $\tau^{t}_{0adv}$ and velocity/strain rate of the - // previous timestep. + // This means that in both cases we can use 'in' to get to the $\tau^{t}_{0adv}$ and velocity/strain rate. // TODO The additional outputs include the reaction rates, so we also have to fill the reaction_rates // if additional_outputs are requested. if (in.current_cell.state() == IteratorState::valid && this->get_timestep_number() > 0 && From 6d6a4bd6c5787e78d486b4e6bb352e1513d2f32e Mon Sep 17 00:00:00 2001 From: anne-glerum Date: Fri, 27 Mar 2026 17:32:03 +0100 Subject: [PATCH 12/15] Update parameter documentation --- doc/sphinx/parameters/Solver_20parameters.md | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/doc/sphinx/parameters/Solver_20parameters.md b/doc/sphinx/parameters/Solver_20parameters.md index 11079d227c1..a3fcb59692e 100644 --- a/doc/sphinx/parameters/Solver_20parameters.md +++ b/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` +: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 “Enable elasticity” is set to true. The strategy parameter is only considered if the parameter “Use operator splitting” is set to true. Units: none. +:::: + ::::{dropdown} __Parameter:__ {ref}`Reaction solver relative tolerance` :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 ‘ARKode’, 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 ‘Reaction solver relative tolerance’. 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 ‘fixed step’ scheme, a method that a forward Euler scheme and a fixed number of steps given by the ‘Reaction time step’ and ‘Reaction time steps per advection step’ 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 ‘ARKode’, 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 ‘Reaction solver relative tolerance’. 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 ‘fixed step’ scheme, a method that uses a forward Euler scheme and a fixed number of steps given by the ‘Reaction time step’ and ‘Reaction time steps per advection step’ parameters. :::: ::::{dropdown} __Parameter:__ {ref}`Reaction time step` From 1349f59e29b275d346e48de0e948f4e83f71b973 Mon Sep 17 00:00:00 2001 From: anne-glerum Date: Mon, 30 Mar 2026 13:49:58 +0200 Subject: [PATCH 13/15] Update tests and benchmarks --- benchmarks/infill_density/infill_ascii.prm | 1 + .../viscoelastic_bending_beam/viscoelastic_bending_beam.prm | 1 + .../viscoelastic_plastic_simple_shear.prm | 1 + .../viscoelastic_plate_flexure/viscoelastic_plate_flexure.prm | 1 + benchmarks/viscoelastic_relaxation/viscoelastic_relaxation.prm | 1 + .../viscoelastic_sheared_torsion.prm | 1 + .../viscoelastic_stress_build_up.prm | 1 + tests/boundary_traction_function_cartesian_free_surface.prm | 1 + tests/boundary_traction_function_spherical_free_surface.prm | 1 + tests/free_surface_timestep_repeat.prm | 1 + tests/maximum_horizontal_compressive_stress_case_elasticity.prm | 1 + tests/newton_postprocess_nonlinear.prm | 1 + tests/skip_refinement.prm | 1 + tests/vep_chemical_compositions.prm | 1 + tests/visco_elastic_top_beam.prm | 1 + ...viscoelastoplastic_yield_plastic_viscous_strain_weakening.prm | 1 + tests/visualization_shear_stress_viscoelastic.prm | 1 + tests/visualization_stress_viscoelastic.prm | 1 + 18 files changed, 18 insertions(+) diff --git a/benchmarks/infill_density/infill_ascii.prm b/benchmarks/infill_density/infill_ascii.prm index fa569dafc69..f5691599da4 100644 --- a/benchmarks/infill_density/infill_ascii.prm +++ b/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 diff --git a/benchmarks/viscoelastic_bending_beam/viscoelastic_bending_beam.prm b/benchmarks/viscoelastic_bending_beam/viscoelastic_bending_beam.prm index 554d1d7d077..f171b5e6d7b 100644 --- a/benchmarks/viscoelastic_bending_beam/viscoelastic_bending_beam.prm +++ b/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 diff --git a/benchmarks/viscoelastic_plastic_simple_shear/viscoelastic_plastic_simple_shear.prm b/benchmarks/viscoelastic_plastic_simple_shear/viscoelastic_plastic_simple_shear.prm index 304ff5f495f..17fe9e5b8e9 100644 --- a/benchmarks/viscoelastic_plastic_simple_shear/viscoelastic_plastic_simple_shear.prm +++ b/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 diff --git a/benchmarks/viscoelastic_plate_flexure/viscoelastic_plate_flexure.prm b/benchmarks/viscoelastic_plate_flexure/viscoelastic_plate_flexure.prm index 343085bfbc0..81d2ecb055d 100644 --- a/benchmarks/viscoelastic_plate_flexure/viscoelastic_plate_flexure.prm +++ b/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 diff --git a/benchmarks/viscoelastic_relaxation/viscoelastic_relaxation.prm b/benchmarks/viscoelastic_relaxation/viscoelastic_relaxation.prm index 30190ff61aa..5cfba766d2c 100644 --- a/benchmarks/viscoelastic_relaxation/viscoelastic_relaxation.prm +++ b/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 diff --git a/benchmarks/viscoelastic_sheared_torsion/viscoelastic_sheared_torsion.prm b/benchmarks/viscoelastic_sheared_torsion/viscoelastic_sheared_torsion.prm index dd5ba914dd3..00b867fba94 100644 --- a/benchmarks/viscoelastic_sheared_torsion/viscoelastic_sheared_torsion.prm +++ b/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 diff --git a/benchmarks/viscoelastic_stress_build_up/viscoelastic_stress_build_up.prm b/benchmarks/viscoelastic_stress_build_up/viscoelastic_stress_build_up.prm index 25c99c07d23..ae087e5d0af 100644 --- a/benchmarks/viscoelastic_stress_build_up/viscoelastic_stress_build_up.prm +++ b/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 diff --git a/tests/boundary_traction_function_cartesian_free_surface.prm b/tests/boundary_traction_function_cartesian_free_surface.prm index e559eca4b07..5d043db735d 100644 --- a/tests/boundary_traction_function_cartesian_free_surface.prm +++ b/tests/boundary_traction_function_cartesian_free_surface.prm @@ -21,6 +21,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 diff --git a/tests/boundary_traction_function_spherical_free_surface.prm b/tests/boundary_traction_function_spherical_free_surface.prm index d4272bdad63..8fc095dedbb 100644 --- a/tests/boundary_traction_function_spherical_free_surface.prm +++ b/tests/boundary_traction_function_spherical_free_surface.prm @@ -21,6 +21,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 diff --git a/tests/free_surface_timestep_repeat.prm b/tests/free_surface_timestep_repeat.prm index 2e8172e663c..adeabfcaf62 100644 --- a/tests/free_surface_timestep_repeat.prm +++ b/tests/free_surface_timestep_repeat.prm @@ -26,6 +26,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 diff --git a/tests/maximum_horizontal_compressive_stress_case_elasticity.prm b/tests/maximum_horizontal_compressive_stress_case_elasticity.prm index fe47f76268f..12ef6f28ec4 100644 --- a/tests/maximum_horizontal_compressive_stress_case_elasticity.prm +++ b/tests/maximum_horizontal_compressive_stress_case_elasticity.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 diff --git a/tests/newton_postprocess_nonlinear.prm b/tests/newton_postprocess_nonlinear.prm index 56995821da4..97c3d50b913 100644 --- a/tests/newton_postprocess_nonlinear.prm +++ b/tests/newton_postprocess_nonlinear.prm @@ -23,6 +23,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 diff --git a/tests/skip_refinement.prm b/tests/skip_refinement.prm index a43cb37f3ef..557400d27ea 100644 --- a/tests/skip_refinement.prm +++ b/tests/skip_refinement.prm @@ -23,6 +23,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 diff --git a/tests/vep_chemical_compositions.prm b/tests/vep_chemical_compositions.prm index 7329d28492a..3656fd1cf68 100644 --- a/tests/vep_chemical_compositions.prm +++ b/tests/vep_chemical_compositions.prm @@ -20,6 +20,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 diff --git a/tests/visco_elastic_top_beam.prm b/tests/visco_elastic_top_beam.prm index 6faad1fe082..91686a601d7 100644 --- a/tests/visco_elastic_top_beam.prm +++ b/tests/visco_elastic_top_beam.prm @@ -22,6 +22,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 subsection Stokes solver parameters set Number of cheap Stokes solver steps = 2000 diff --git a/tests/viscoelastoplastic_yield_plastic_viscous_strain_weakening.prm b/tests/viscoelastoplastic_yield_plastic_viscous_strain_weakening.prm index 307ad02e709..fdf3491c0fb 100644 --- a/tests/viscoelastoplastic_yield_plastic_viscous_strain_weakening.prm +++ b/tests/viscoelastoplastic_yield_plastic_viscous_strain_weakening.prm @@ -17,6 +17,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 diff --git a/tests/visualization_shear_stress_viscoelastic.prm b/tests/visualization_shear_stress_viscoelastic.prm index 379b38f534c..48c276fe82d 100644 --- a/tests/visualization_shear_stress_viscoelastic.prm +++ b/tests/visualization_shear_stress_viscoelastic.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 diff --git a/tests/visualization_stress_viscoelastic.prm b/tests/visualization_stress_viscoelastic.prm index dec19272a78..4ef8ae91126 100644 --- a/tests/visualization_stress_viscoelastic.prm +++ b/tests/visualization_stress_viscoelastic.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 From 590769d9658aa0ac0df433ee7effb8dcf571e480 Mon Sep 17 00:00:00 2001 From: anne-glerum Date: Mon, 30 Mar 2026 14:10:45 +0200 Subject: [PATCH 14/15] Update comments rheology --- source/material_model/rheology/elasticity.cc | 27 ++++++++++---------- 1 file changed, 13 insertions(+), 14 deletions(-) diff --git a/source/material_model/rheology/elasticity.cc b/source/material_model/rheology/elasticity.cc index fc0fec1ee7f..062626ac284 100644 --- a/source/material_model/rheology/elasticity.cc +++ b/source/material_model/rheology/elasticity.cc @@ -357,7 +357,7 @@ namespace aspect // Create the ReactionRateOutputs that are necessary for the operator splitting // step (either on the fields or directly on the particles) // that sets both sets of stresses to the total stress of the - // previous timestep. + // current timestep after the nonlinear solver is done. if (out.template has_additional_output_object>() == false && (this->get_parameters().use_operator_splitting || (this->get_parameters().mapped_particle_properties).count(this->introspection().compositional_index_for_name("ve_stress_xx")))) { @@ -425,7 +425,7 @@ namespace aspect // (the force term goes into the rhs of the momentum equation). // This happens after the advection equations have been solved, and hence in.composition // contains the rotated and advected stresses $tau^{0adv}$. - // Only at the beginning of the next timestep do we add the stress update of the + // Only after the nonlinear solver loop do we add the stress update of the // current timestep to the stress stored in the compositional fields, giving // $\tau{t+\Delta t_c}$ with $t+\Delta t_c$ being the current timestep. const SymmetricTensor<2,dim> stress_0_advected (Utilities::Tensors::to_symmetric_tensor(&in.composition[i][stress_start_index], @@ -467,7 +467,7 @@ namespace aspect in.strain_rate[i], stress_0_advected, stress_old, effective_creep_viscosity, average_elastic_shear_moduli[i]); } - // Apply the stress update to get the total stress of timestep t. + // Apply the stress update to get the total stress of timestep $t+\Delta t_c$. const SymmetricTensor<2, dim> stress = 2. * effective_creep_viscosity * deviatoric_strain_rate + viscosity_ratio * stress_0_advected + (1. - timestep_ratio) * (1. - viscosity_ratio) * stress_old; @@ -520,7 +520,7 @@ namespace aspect const double elastic_viscosity = calculate_elastic_viscosity(average_elastic_shear_moduli[i]); - // Apply the stress update to get the total deviatoric stress of timestep t. + // Apply the stress update to get the total deviatoric stress of timestep $t+\Delta t_c$. elastic_additional_out->deviatoric_stress[i] = 2. * eta * deviatoric_strain_rate + eta / elastic_viscosity * stress_0_advected + (1. - timestep_ratio) * (1. - eta / elastic_viscosity) * stress_old; elastic_additional_out->elastic_shear_moduli[i] = average_elastic_shear_moduli[i]; elastic_additional_out->elastic_viscosity[i] = elastic_viscosity; @@ -603,8 +603,8 @@ namespace aspect } // The following function computes the reaction rates for the operator - // splitting step that at the end of the timestep $t$ updates the - // stored compositions $tau^{0\mathrm{adv}}$ at time $t$ to $tau^{t}$. + // splitting step that at the end of the timestep $t+\Delta t_c$ updates the + // stored compositions $tau^{0\mathrm{adv}}$ at time $t+\Delta t_c$ to $tau^{t+\Delta t_c}$. // This update consists of the stress change resulting from system evolution, // but does not advect or rotate the stress tensor. Advection is done by // solving the advection equation and the stress tensor is rotated through @@ -638,14 +638,14 @@ namespace aspect // At the moment when the reaction rates are required (after the nonlinear solver loop), // the solution vector 'solution' holds the stress from the previous timestep - // advected and rotated into the new position of the current timestep, so $\tau^{t}_{0adv}$. + // advected and rotated into the new position of the current timestep, so $\tau^{t+\Delta t_c}_{0adv}$. // // In case fields are used to track the stresses, MaterialModelInputs are based on 'solution' // when calling the MaterialModel for the reaction rates. When particles are used, MaterialModelInputs // for this function are filled with the solution (including for the strain rate), except for the // compositions that represent the stress tensor components, these are taken directly from the // particles in the property plugin by default (although this can be changed from the input file). - // This means that in both cases we can use 'in' to get to the $\tau^{t}_{0adv}$ and velocity/strain rate. + // This means that in both cases we can use 'in' to get to the $\tau^{t+\Delta t_c}_{0adv}$ and velocity/strain rate. // TODO The additional outputs include the reaction rates, so we also have to fill the reaction_rates // if additional_outputs are requested. if (in.current_cell.state() == IteratorState::valid && this->get_timestep_number() > 0 && @@ -700,7 +700,7 @@ namespace aspect const SymmetricTensor<2, dim> stress_old (Utilities::Tensors::to_symmetric_tensor(&in.composition[i][stress_start_index+n_independent_components], &in.composition[i][stress_start_index+n_independent_components]+n_independent_components)); - // $\eta^{t}_{effcreep}$. This viscosity has been calculated with the timestep_ratio dtc/dte. + // $\eta^{t+\Delta t_c}_{effcreep}$. This viscosity has been calculated with the timestep_ratio dtc/dte. const double effective_creep_viscosity = effective_creep_viscosities[i]; // $\eta_{el} = G \Delta t_c$. @@ -712,7 +712,7 @@ namespace aspect // The ratio between the computational and elastic timestep $\frac{\Delta t_c} / {\Delta t_{el}}$. const double timestep_ratio = calculate_timestep_ratio(); - // Compute the total stress at time t. + // Compute the total stress at time $t+\Delta t_c$. const SymmetricTensor<2, dim> stress_t = 2. * effective_creep_viscosity * Utilities::Tensors::consistent_deviator(in.strain_rate[i]) + effective_creep_viscosity / elastic_viscosity * stress_0_t @@ -893,8 +893,8 @@ namespace aspect { // Get the compositional fields from the previous timestep $t$. // The 'old_solution' has been updated to the full stress tensor - // of time $t$ by the operator splitting step at the beginning - // of the current timestep. + // of time $t$ by the operator splitting step after the nonlinear + // solver loop of the previous timestep. std::vector old_solution_values(this->get_fe().dofs_per_cell); in.current_cell->get_dof_values(this->get_old_solution(), old_solution_values.begin(), @@ -936,8 +936,7 @@ namespace aspect // which include the stresses stored on the particles. // The computation of the reaction terms (during which this function is called) happens // after the particles have been restored to their position and values - // at the beginning of the timestep (i.e., the position and values of the previous timestep) and after they - // have been updated to the full stress of the previous timestep by the reaction rates. + // at the beginning of the timestep (i.e., the position and values of the previous timestep). for (unsigned int i = 0; i < in.n_evaluation_points(); ++i) stress_t[i] = Utilities::Tensors::to_symmetric_tensor(&in.composition[i][stress_start_index], &in.composition[i][stress_start_index]+n_independent_components); From e6c9e4df4b93c71a6f1ad011a8ef477ed5ae42e6 Mon Sep 17 00:00:00 2001 From: anne-glerum Date: Mon, 30 Mar 2026 17:37:11 +0200 Subject: [PATCH 15/15] Update two more prms --- .../viscoelastic/free_surface_VE_cylinder_2D_loading.prm | 1 + .../viscoelastic/free_surface_VE_cylinder_3D_loading.prm | 1 + 2 files changed, 2 insertions(+) diff --git a/benchmarks/free_surface_tractions/viscoelastic/free_surface_VE_cylinder_2D_loading.prm b/benchmarks/free_surface_tractions/viscoelastic/free_surface_VE_cylinder_2D_loading.prm index 0fcca6715e6..0fb5f174fae 100644 --- a/benchmarks/free_surface_tractions/viscoelastic/free_surface_VE_cylinder_2D_loading.prm +++ b/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 diff --git a/benchmarks/free_surface_tractions/viscoelastic/free_surface_VE_cylinder_3D_loading.prm b/benchmarks/free_surface_tractions/viscoelastic/free_surface_VE_cylinder_3D_loading.prm index f6add59ac8c..ad96b33f360 100644 --- a/benchmarks/free_surface_tractions/viscoelastic/free_surface_VE_cylinder_3D_loading.prm +++ b/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