Merge pull request #6659 from gassmoeller/iterated_advection_no_stokes

Add Iterated advection, no stokes solver scheme

Author: jdannberg

doc/modules/changes/20250902_gassmoeller

@@ -0,0 +1,5 @@
+New: Added a new nonlinear solver scheme 'iterated Advection, no Stokes'
+that iterates the advection equation and uses a prescribed Stokes
+solution.
+<br>
+(Rene Gassmoeller, 2025/09/02)

include/aspect/parameters.h

@@ -74,6 +74,7 @@ namespace aspect
         single_Advection_iterated_Stokes,
         single_Advection_iterated_defect_correction_Stokes,
         single_Advection_iterated_Newton_Stokes,
+        iterated_Advection_no_Stokes,
         iterated_Advection_and_Stokes,
         iterated_Advection_and_defect_correction_Stokes,
         iterated_Advection_and_Newton_Stokes

include/aspect/simulator.h

@@ -575,6 +575,19 @@ namespace aspect
        */
       void solve_single_advection_iterated_newton_stokes (bool use_newton_iterations);
 
+      /**
+       * This function implements one scheme for the various
+       * steps necessary to assemble and solve the nonlinear problem.
+       *
+       * The `iterated Advection, no Stokes' scheme iterates the temperature and other
+       * advection systems and instead of solving for the Stokes system,
+       * a prescribed velocity and pressure are used."
+       *
+       * This function is implemented in
+       * <code>source/simulator/solver_schemes.cc</code>.
+       */
+      void solve_iterated_advection_no_stokes ();
+
       /**
        * This function implements one scheme for the various
        * steps necessary to assemble and solve the nonlinear problem.

source/material_model/entropy_model.cc

@@ -45,6 +45,7 @@ namespace aspect
         // models splits temperature diffusion from entropy advection.
         switch (parameters.nonlinear_solver)
           {
+            case Parameters<dim>::NonlinearSolver::Kind::iterated_Advection_no_Stokes:
             case Parameters<dim>::NonlinearSolver::Kind::iterated_Advection_and_Stokes:
             case Parameters<dim>::NonlinearSolver::Kind::iterated_Advection_and_defect_correction_Stokes:
             case Parameters<dim>::NonlinearSolver::Kind::iterated_Advection_and_Newton_Stokes:
@@ -76,7 +77,9 @@ namespace aspect
       CitationInfo::add("entropy");
 
       AssertThrow (this->get_parameters().formulation_mass_conservation ==
-                   Parameters<dim>::Formulation::MassConservation::projected_density_field,
+                   Parameters<dim>::Formulation::MassConservation::projected_density_field ||
+                   this->get_parameters().nonlinear_solver == Parameters<dim>::NonlinearSolver::single_Advection_no_Stokes ||
+                   this->get_parameters().nonlinear_solver == Parameters<dim>::NonlinearSolver::iterated_Advection_no_Stokes,
                    ExcMessage("The 'entropy model' material model was only tested with the "
                               "'projected density field' approximation "
                               "for the mass conservation equation, which is not selected."));

source/simulator/core.cc

@@ -378,17 +378,18 @@ namespace aspect
     boundary_velocity_manager.parse_parameters (prm);
 
     // Make sure we only have a prescribed Stokes plugin if needed
-    if (parameters.nonlinear_solver == NonlinearSolver::single_Advection_no_Stokes)
+    if (parameters.nonlinear_solver == NonlinearSolver::single_Advection_no_Stokes ||
+        parameters.nonlinear_solver == NonlinearSolver::iterated_Advection_no_Stokes)
       {
         AssertThrow(prescribed_stokes_solution.get()!=nullptr,
-                    ExcMessage("For the 'single Advection, no Stokes' solver scheme you need to provide a Stokes plugin!")
+                    ExcMessage("For the selected nonlinear solver scheme you need to provide a prescribed Stokes plugin!")
                    );
       }
     else
       {
         AssertThrow(prescribed_stokes_solution.get()==nullptr,
                     ExcMessage("The prescribed stokes plugin you selected only works with the solver "
-                               "scheme 'single Advection, no Stokes'.")
+                               "scheme 'single Advection, no Stokes' or 'iterated Advection, no Stokes'.")
                    );
       }
 
@@ -872,6 +873,7 @@ namespace aspect
           case Parameters<dim>::NonlinearSolver::Kind::single_Advection_iterated_Stokes:
           case Parameters<dim>::NonlinearSolver::Kind::single_Advection_iterated_defect_correction_Stokes:
           case Parameters<dim>::NonlinearSolver::Kind::single_Advection_iterated_Newton_Stokes:
+          case Parameters<dim>::NonlinearSolver::Kind::iterated_Advection_no_Stokes:
           case Parameters<dim>::NonlinearSolver::Kind::iterated_Advection_and_Stokes:
           case Parameters<dim>::NonlinearSolver::Kind::iterated_Advection_and_defect_correction_Stokes:
           case Parameters<dim>::NonlinearSolver::Kind::iterated_Advection_and_Newton_Stokes:
@@ -911,6 +913,7 @@ namespace aspect
 
           case Parameters<dim>::NonlinearSolver::Kind::no_Advection_no_Stokes:
           case Parameters<dim>::NonlinearSolver::Kind::single_Advection_no_Stokes:
+          case Parameters<dim>::NonlinearSolver::Kind::iterated_Advection_no_Stokes:
             return false;
         }
       Assert(false, ExcNotImplemented());
@@ -1981,6 +1984,12 @@ namespace aspect
               break;
             }
 
+            case NonlinearSolver::iterated_Advection_no_Stokes:
+            {
+              solve_iterated_advection_no_stokes();
+              break;
+            }
+
             case NonlinearSolver::iterated_Advection_and_Stokes:
             {
               solve_iterated_advection_and_stokes();

source/simulator/helper_functions.cc

@@ -2577,14 +2577,15 @@ namespace aspect
                                  +
                                  "> is listed for a boundary condition, but is not used by the geometry model."));
 
-    if (parameters.nonlinear_solver == NonlinearSolver::single_Advection_no_Stokes)
+    if (parameters.nonlinear_solver == NonlinearSolver::single_Advection_no_Stokes ||
+        parameters.nonlinear_solver == NonlinearSolver::iterated_Advection_no_Stokes)
       {
         // make sure that there are no listed velocity boundary conditions
         for (unsigned int i=0; i<4; ++i)
           AssertThrow (boundary_indicator_lists[i].empty(),
-                       ExcMessage ("With the solver scheme `single Advection, no Stokes', "
-                                   "one cannot set boundary conditions for velocity or traction, "
-                                   "but a boundary condition has been set."));
+                       ExcMessage ("With the solver schemes `single Advection, no Stokes' or "
+                                   "'iterated Advection, no Stokes' one cannot set boundary conditions "
+                                   "for velocity or traction, but a boundary condition has been set."));
       }
   }
 

source/simulator/parameters.cc

@@ -216,6 +216,7 @@ namespace aspect
                                                "single Advection, iterated Stokes|"
                                                "single Advection, iterated defect correction Stokes|"
                                                "single Advection, iterated Newton Stokes|"
+                                               "iterated Advection, no Stokes|"
                                                "iterated Advection and Stokes|"
                                                "iterated Advection and defect correction Stokes|"
                                                "iterated Advection and Newton Stokes";
@@ -254,6 +255,10 @@ namespace aspect
                        "the temperature and composition equations once at the beginning of each time step and "
                        "then iterates out the solution of the Stokes equation, using Newton iterations for the "
                        "Stokes system.\n"
+                       "The `iterated Advection, no Stokes' scheme iterates the temperature and other advection "
+                       "equations, and instead of solving for the Stokes system, a prescribed "
+                       "velocity and pressure are used. This is useful for kinematic models and advection benchmarks "
+                       "with nonlinear processes in the advection equations.\n"
                        "The `iterated Advection and Stokes' scheme iterates out the nonlinearity "
                        "by alternating the solution of the temperature, composition and Stokes systems.\n"
                        "The `iterated Advection and defect correction Stokes' scheme iterates by alternating the "
@@ -1595,6 +1600,8 @@ namespace aspect
         nonlinear_solver = NonlinearSolver::single_Advection_iterated_defect_correction_Stokes;
       else if (solver_scheme == "single Advection, iterated Newton Stokes")
         nonlinear_solver = NonlinearSolver::single_Advection_iterated_Newton_Stokes;
+      else if (solver_scheme == "iterated Advection, no Stokes")
+        nonlinear_solver = NonlinearSolver::iterated_Advection_no_Stokes;
       else if (solver_scheme == "iterated Advection and Stokes")
         nonlinear_solver = NonlinearSolver::iterated_Advection_and_Stokes;
       else if (solver_scheme == "iterated Advection and defect correction Stokes")

source/simulator/solver_schemes.cc

@@ -1222,6 +1222,120 @@ namespace aspect
 
 
 
+  template <int dim>
+  void Simulator<dim>::solve_iterated_advection_no_stokes ()
+  {
+    double initial_temperature_residual = 0;
+    std::vector<double> initial_composition_residual (introspection.n_compositional_fields,0);
+
+    const unsigned int max_nonlinear_iterations =
+      (pre_refinement_step < parameters.initial_adaptive_refinement)
+      ?
+      std::min(parameters.max_nonlinear_iterations,
+               parameters.max_nonlinear_iterations_in_prerefinement)
+      :
+      parameters.max_nonlinear_iterations;
+
+    SolverControl nonlinear_solver_control(max_nonlinear_iterations,
+                                           parameters.nonlinear_tolerance);
+
+    double relative_residual = std::numeric_limits<double>::max();
+    nonlinear_iteration = 0;
+
+    do
+      {
+        for (auto &particle_manager : particle_managers)
+          {
+            // Restore particles through stored copy of particle handler,
+            // but only if they have already been displaced in a nonlinear
+            // iteration (in the assemble_and_solve_composition call).
+            if (nonlinear_iteration > 0)
+              particle_manager.restore_particles();
+
+            // Apply a particle update if required by the particle properties.
+            // Apply the update even if nonlinear_iteration == 0,
+            // because the signal could be used, for example, to apply operator
+            // splitting on the particles, and this would need to be
+            // applied at the beginning of the timestep and after
+            // every restore_particles().
+            signals.post_restore_particles(particle_manager);
+          }
+
+        const double relative_temperature_residual =
+          assemble_and_solve_temperature(initial_temperature_residual,
+                                         nonlinear_iteration == 0 ? &initial_temperature_residual : nullptr);
+
+        const std::vector<double>  relative_composition_residual =
+          assemble_and_solve_composition(initial_composition_residual,
+                                         nonlinear_iteration == 0 ? &initial_composition_residual : nullptr);
+
+        // write the residual output in the same order as the solutions
+        pcout << "      Relative nonlinear residuals (temperature"
+              << (introspection.n_compositional_fields > 0 ? ", compositional fields" : "")
+              << "): "
+              << relative_temperature_residual;
+        for (unsigned int c=0; c<introspection.n_compositional_fields; ++c)
+          pcout << ", " << relative_composition_residual[c];
+        pcout << std::endl;
+
+        // Find the maximum residual of the individual equations
+        relative_residual = relative_temperature_residual;
+        for (unsigned int c=0; c<introspection.n_compositional_fields; ++c)
+          relative_residual = std::max(relative_composition_residual[c],relative_residual);
+
+        pcout << "      Relative nonlinear residual (total system) after nonlinear iteration " << nonlinear_iteration+1
+              << ": " << relative_residual
+              << std::endl
+              << std::endl;
+
+        if (parameters.run_postprocessors_on_nonlinear_iterations)
+          postprocess ();
+
+        ++nonlinear_iteration;
+      }
+    while (nonlinear_solver_control.check(nonlinear_iteration, relative_residual) == SolverControl::iterate);
+
+    // Assign Stokes solution
+    {
+      TimerOutput::Scope timer (computing_timer, "Interpolate Stokes solution");
+
+      LinearAlgebra::BlockVector distributed_stokes_solution (introspection.index_sets.system_partitioning, mpi_communicator);
+
+      VectorFunctionFromVectorFunctionObject<dim> func(
+        [&](const Point<dim> &p, Vector<double> &result)
+      {
+        prescribed_stokes_solution->stokes_solution(p, result);
+      },
+      0,
+      parameters.include_melt_transport ? 2*dim+3 : dim+1, // velocity and pressure
+      introspection.n_components);
+
+      VectorTools::interpolate (*mapping, dof_handler, func, distributed_stokes_solution);
+
+      // distribute hanging node and other constraints
+      current_constraints.distribute (distributed_stokes_solution);
+
+      solution.block(introspection.block_indices.velocities) =
+        distributed_stokes_solution.block(introspection.block_indices.velocities);
+      solution.block(introspection.block_indices.pressure) =
+        distributed_stokes_solution.block(introspection.block_indices.pressure);
+
+      if (parameters.include_melt_transport)
+        {
+          const unsigned int block_u_f = introspection.variable("fluid velocity").block_index;
+          const unsigned int block_p_f = introspection.variable("fluid pressure").block_index;
+          solution.block(block_u_f) = distributed_stokes_solution.block(block_u_f);
+          solution.block(block_p_f) = distributed_stokes_solution.block(block_p_f);
+        }
+
+    }
+
+    AssertThrow(nonlinear_solver_control.last_check() != SolverControl::failure, ExcNonlinearSolverNoConvergence());
+    signals.post_nonlinear_solver(nonlinear_solver_control);
+  }
+
+
+
   template <int dim>
   void Simulator<dim>::solve_iterated_advection_and_stokes ()
   {
@@ -1490,6 +1604,7 @@ namespace aspect
   template void Simulator<dim>::solve_single_advection_iterated_stokes(); \
   template void Simulator<dim>::solve_single_advection_iterated_defect_correction_stokes(); \
   template void Simulator<dim>::solve_single_advection_iterated_newton_stokes(bool); \
+  template void Simulator<dim>::solve_iterated_advection_no_stokes(); \
   template void Simulator<dim>::solve_iterated_advection_and_stokes(); \
   template void Simulator<dim>::solve_iterated_advection_and_defect_correction_stokes(); \
   template void Simulator<dim>::solve_iterated_advection_and_newton_stokes(bool); \

tests/entropy_multicomponent_11_edge_case.prm

@@ -11,7 +11,7 @@ set Use years instead of seconds           = true
 set End time                               = 4e6
 set Maximum time step                      = 2e5
 set Output directory                       = output-entropy-multi-edge-case-change-adiabatic-p
-set Nonlinear solver scheme                = iterated Advection and Stokes 
+set Nonlinear solver scheme                = iterated Advection, no Stokes
 set Use operator splitting                 = true
 
 subsection Solver parameters
@@ -22,13 +22,6 @@ subsection Solver parameters
   end
 end
 
-set Surface pressure                       = 1.2e9
-set Pressure normalization                 = volume
-
-subsection Formulation
-  set Mass conservation = projected density field
-end
-
 subsection Geometry model
   set Model name = box
 
@@ -41,14 +34,12 @@ end
 
 subsection Adiabatic conditions model
  subsection Compute profile
- 
   set Use surface condition function = true
 
   	subsection Surface condition function
-    set Variable names = x,t 
-    set Function expression = 1.2e9+4e2*t;0
+      set Variable names = x,t
+      set Function expression = 1.2e9+4e2*t;0
     end
-
  end
 end
 
@@ -65,10 +56,6 @@ subsection Compositional fields
   set Compositional field methods = field, field, field, prescribed field
 end
 
-subsection Boundary velocity model
-  set Tangential velocity boundary indicators =  left, right, top, bottom
-end
-
 subsection Gravity model
   set Model name = vertical
 
@@ -81,10 +68,8 @@ subsection Initial composition model
   set List of model names = function
 
   subsection Function
-   
     set Function expression = 0.3; 1610; 2006; 0.0
-
-    end
+  end
 end
 
 subsection Material model
@@ -110,6 +95,19 @@ subsection Mesh refinement
   set Time steps between mesh refinement       = 0
 end
 
+subsection Prescribed Stokes solution
+  set Model name = function
+
+  # The same function as we use for the adiabatic pressure
+  subsection Pressure function
+    set Function expression = 1.2e9+4e2*t
+  end
+
+  subsection Velocity function
+    set Function expression = 0.0; 0.0
+  end
+end
+
 subsection Postprocess
   set List of postprocessors = visualization, velocity statistics, temperature statistics, mass flux statistics, composition statistics