Refactor the Newton line search.

Author: arnoldkk13

include/aspect/simulator.h

@@ -1337,6 +1337,16 @@ namespace aspect
       void interpolate_onto_velocity_system(const TensorFunction<1,dim> &func,
                                             LinearAlgebra::Vector &vec) const;
 
+      /**
+       * Perform a Newton line search to determine the optimal step length
+       * along the search direction.
+       *
+       * This function is implemented in
+       * <code>source/simulator/helper_functions.cc</code>
+       */
+      double perform_line_search(const DefectCorrectionResiduals &dcr,
+                                 const bool use_picard,
+                                 LinearAlgebra::BlockVector &search_direction);
 
       /**
        * Add constraints to the given @p constraints object that are required

source/simulator/helper_functions.cc

@@ -2629,6 +2629,92 @@ namespace aspect
     return initial_newton_residual;
   }
 
+  /**
+   * Performs a Newton line search for the do_one_defect_correction_Stokes_step function
+   */
+  template <int dim>
+  double Simulator<dim>::perform_line_search(const DefectCorrectionResiduals &dcr,
+                                             const bool use_picard,
+                                             LinearAlgebra::BlockVector &search_direction)
+  {
+    LinearAlgebra::BlockVector backup_linearization_point(introspection.index_sets.stokes_partitioning, mpi_communicator);
+    double step_length_factor = 1.0;
+    unsigned int line_search_iteration = 0;
+    // Many parts of the solver depend on the block layout (velocity = 0,
+    // pressure = 1). For example the linearized_stokes_initial_guess vector or the StokesBlock matrix
+    // wrapper. Let us make sure that this holds (and shorten their names):
+    const unsigned int pressure_block_index = (parameters.include_melt_transport) ?
+                                              introspection.variable("fluid pressure").block_index
+                                              : introspection.block_indices.pressure;
+    const unsigned int velocity_block_index = introspection.block_indices.velocities;
+    Assert(velocity_block_index == 0, ExcNotImplemented());
+    Assert(pressure_block_index == 1, ExcNotImplemented());
+    (void) pressure_block_index;
+    // Do the loop for the line search at least once with the full step length.
+    // If line search is disabled we will exit the loop in the first iteration.
+
+    do
+      {
+        if (line_search_iteration == 0)
+          {
+            // backup our starting point for the line search
+            backup_linearization_point.block(pressure_block_index) = current_linearization_point.block(pressure_block_index);
+            backup_linearization_point.block(velocity_block_index) = current_linearization_point.block(velocity_block_index);
+          }
+        else
+          {
+            // undo the last iteration and try again with smaller step length
+            current_linearization_point.block(pressure_block_index) = backup_linearization_point.block(pressure_block_index);
+            current_linearization_point.block(velocity_block_index) = backup_linearization_point.block(velocity_block_index);
+            search_direction.block(pressure_block_index) *= step_length_factor;
+            search_direction.block(velocity_block_index) *= step_length_factor;
+          }
+
+        // Update the current linearization point with the search direction
+        current_linearization_point.block(pressure_block_index) += search_direction.block(pressure_block_index);
+        current_linearization_point.block(velocity_block_index) += search_direction.block(velocity_block_index);
+
+        // Rebuild the rhs to determine the new residual.
+        assemble_newton_stokes_matrix = rebuild_stokes_preconditioner = false;
+        rebuild_stokes_matrix = !boundary_velocity_manager.get_prescribed_boundary_velocity_indicators().empty();
+
+        assemble_stokes_system();
+
+        const double test_velocity_residual = system_rhs.block(velocity_block_index).l2_norm();
+        const double test_pressure_residual = system_rhs.block(pressure_block_index).l2_norm();
+        const double test_residual = std::sqrt(test_velocity_residual * test_velocity_residual
+                                               + test_pressure_residual * test_pressure_residual);
+
+        // Determine if the residual has decreased sufficiently.
+        const double alpha = 1e-4;
+        if (test_residual < (1.0 - alpha * step_length_factor) * dcr.residual
+            ||
+            line_search_iteration >= newton_handler->parameters.max_newton_line_search_iterations
+            ||
+            use_picard)
+          {
+            pcout << "      Newton system information: Norm of the rhs: " << test_residual
+                  << ", Derivative scaling factor: " << newton_handler->parameters.newton_derivative_scaling_factor << std::endl;
+            return test_residual;
+          }
+        else
+          {
+            pcout << "   Line search iteration " << line_search_iteration << ", with norm of the rhs "
+                  << test_residual << " and going to " << (1.0 - alpha * step_length_factor) * dcr.residual
+                  << ", relative residual: " << test_residual/dcr.initial_residual << std::endl;
+
+            // The current search direction has not decreased the residual
+            // enough, so we take a smaller step and try again.
+            // TODO: make a parameter out of this.
+            step_length_factor = 2.0/3.0;
+          }
+
+        ++line_search_iteration;
+        Assert(line_search_iteration <= newton_handler->parameters.max_newton_line_search_iterations,
+               ExcInternalError());
+      }
+    while (true);
+  }
 
 
   template <int dim>
@@ -2779,6 +2865,9 @@ namespace aspect
   template void Simulator<dim>::restore_outflow_boundary_ids(const unsigned int boundary_id_offset); \
   template void Simulator<dim>::check_consistency_of_boundary_conditions() const; \
   template double Simulator<dim>::compute_initial_newton_residual(); \
+  template double Simulator<dim>::perform_line_search(const DefectCorrectionResiduals &dcr, \
+                                                      const bool use_picard, \
+                                                      LinearAlgebra::BlockVector &search_direction); \
   template double Simulator<dim>::compute_Eisenstat_Walker_linear_tolerance(const bool EisenstatWalkerChoiceOne, \
                                                                             const double maximum_linear_stokes_solver_tolerance, \
                                                                             const double linear_stokes_solver_tolerance, \

source/simulator/solver_schemes.cc

@@ -530,6 +530,7 @@ namespace aspect
   }
 
 
+
   template <int dim>
   void Simulator<dim>::do_one_defect_correction_Stokes_step(DefectCorrectionResiduals &dcr,
                                                             const bool use_picard)
@@ -702,7 +703,9 @@ namespace aspect
                                                             dcr.residual,
                                                             dcr.residual_old);
 
-                pcout << "   The linear solver tolerance is set to " << parameters.linear_stokes_solver_tolerance << std::endl;
+                pcout << "   The linear solver tolerance is set to "
+                      << parameters.linear_stokes_solver_tolerance
+                      << ". ";
               }
           }
 
@@ -737,79 +740,11 @@ namespace aspect
         dcr.residual_old = dcr.residual;
 
         // We have computed an update. Prepare for a line search.
-        LinearAlgebra::BlockVector backup_linearization_point(introspection.index_sets.stokes_partitioning, mpi_communicator);
-
-        double step_length_factor = 1.0;
-        unsigned int line_search_iteration = 0;
-
-        // Do the loop for the line search at least once with the full step length.
-        // If line search is disabled we will exit the loop in the first iteration.
-        do
-          {
-            if (line_search_iteration == 0)
-              {
-                // backup our starting point for the line search
-                backup_linearization_point.block(pressure_block_index) = current_linearization_point.block(pressure_block_index);
-                backup_linearization_point.block(velocity_block_index) = current_linearization_point.block(velocity_block_index);
-              }
-            else
-              {
-                // undo the last iteration and try again with smaller step length
-                current_linearization_point.block(pressure_block_index) = backup_linearization_point.block(pressure_block_index);
-                current_linearization_point.block(velocity_block_index) = backup_linearization_point.block(velocity_block_index);
-                search_direction.block(pressure_block_index) *= step_length_factor;
-                search_direction.block(velocity_block_index) *= step_length_factor;
-              }
-
-            // Update the current linearization point with the search direction
-            current_linearization_point.block(pressure_block_index) += search_direction.block(pressure_block_index);
-            current_linearization_point.block(velocity_block_index) += search_direction.block(velocity_block_index);
-
-            // Rebuild the rhs to determine the new residual.
-            assemble_newton_stokes_matrix = rebuild_stokes_preconditioner = false;
-            rebuild_stokes_matrix = !boundary_velocity_manager.get_prescribed_boundary_velocity_indicators().empty();
-
-            assemble_stokes_system();
-
-            const double test_velocity_residual = system_rhs.block(velocity_block_index).l2_norm();
-            const double test_pressure_residual = system_rhs.block(pressure_block_index).l2_norm();
-            const double test_residual = std::sqrt(test_velocity_residual * test_velocity_residual
-                                                   + test_pressure_residual * test_pressure_residual);
-
-            // Determine if the residual has decreased sufficiently.
-            const double alpha = 1e-4;
-            if (test_residual < (1.0 - alpha * step_length_factor) * dcr.residual
-                ||
-                line_search_iteration >= newton_handler->parameters.max_newton_line_search_iterations
-                ||
-                use_picard)
-              {
-                pcout << "      Newton system information: Norm of the rhs: " << test_residual
-                      << ", Derivative scaling factor: " << newton_handler->parameters.newton_derivative_scaling_factor << std::endl;
-                dcr.residual = test_residual;
-                break;
-              }
-            else
-              {
-                pcout << "   Line search iteration " << line_search_iteration << ", with norm of the rhs "
-                      << test_residual << " and going to " << (1.0 - alpha * step_length_factor) * dcr.residual
-                      << ", relative residual: " << test_residual/dcr.initial_residual << std::endl;
-
-                // The current search direction has not decreased the residual
-                // enough, so we take a smaller step and try again.
-                // TODO: make a parameter out of this.
-                step_length_factor = 2.0/3.0;
-              }
-
-            ++line_search_iteration;
-            Assert(line_search_iteration <= newton_handler->parameters.max_newton_line_search_iterations,
-                   ExcInternalError());
-          }
-        while (true);
+        dcr.residual = perform_line_search(dcr, use_picard, search_direction);
       }
 
 
-    if (use_picard == true)
+    if (use_picard)
       {
         // When we are using (defect corrected) Picard, keep the
         // newton_derivative_scaling_factor at zero. The newton_derivative_scaling_factor