code improvements

Author: gassmoeller

include/aspect/simulator.h

@@ -2056,7 +2056,7 @@ namespace aspect
       bool                                                      rebuild_sparsity_and_matrices;
       bool                                                      rebuild_stokes_matrix;
       bool                                                      rebuild_stokes_preconditioner;
-      bool                                                      assemble_defect_correction_stokes_system;
+      bool                                                      assemble_newton_stokes_system;
 
       /**
        * @}

source/simulator/assembly.cc

@@ -259,24 +259,24 @@ namespace aspect
 
     // Set up the set of assemblers depending on the mode of equations
     if (!parameters.include_melt_transport
-        && !assemble_defect_correction_stokes_system)
+        && !assemble_newton_stokes_system)
       {
         set_advection_assemblers();
         set_stokes_assemblers();
       }
     else if (parameters.include_melt_transport
-             && !assemble_defect_correction_stokes_system)
+             && !assemble_newton_stokes_system)
       {
         melt_handler->set_assemblers(*assemblers);
       }
     else if (!parameters.include_melt_transport
-             && assemble_defect_correction_stokes_system)
+             && assemble_newton_stokes_system)
       {
         set_advection_assemblers();
         newton_handler->set_assemblers(*assemblers);
       }
     else if (parameters.include_melt_transport
-             && assemble_defect_correction_stokes_system)
+             && assemble_newton_stokes_system)
       {
         AssertThrow (false,
                      ExcMessage("The melt implementation does not support the Newton solver at the moment."));
@@ -337,7 +337,7 @@ namespace aspect
     scratch.material_model_inputs.requested_properties
       = MaterialModel::MaterialProperties::equation_of_state_properties |
         MaterialModel::MaterialProperties::viscosity |
-        (parameters.include_melt_transport || assemble_defect_correction_stokes_system
+        (parameters.include_melt_transport || assemble_newton_stokes_system
          ?
          MaterialModel::MaterialProperties::additional_outputs
          :
@@ -621,7 +621,7 @@ namespace aspect
 
     // initialize the material model data on the cell
     const bool need_viscosity =
-      assemble_defect_correction_stokes_system || this->parameters.enable_prescribed_dilation || rebuild_stokes_matrix;
+      assemble_newton_stokes_system || this->parameters.enable_prescribed_dilation || rebuild_stokes_matrix;
 
     scratch.material_model_inputs.reinit  (scratch.finite_element_values,
                                            cell,
@@ -649,7 +649,7 @@ namespace aspect
 
     scratch.finite_element_values[introspection.extractors.velocities].get_function_values(current_linearization_point,
         scratch.velocity_values);
-    if (assemble_defect_correction_stokes_system)
+    if (assemble_newton_stokes_system)
       scratch.finite_element_values[introspection.extractors.velocities].get_function_divergences(current_linearization_point,scratch.velocity_divergence);
     if (parameters.formulation_mass_conservation == Parameters<dim>::Formulation::MassConservation::hydrostatic_compression)
       scratch.finite_element_values[introspection.extractors.temperature].get_function_gradients(current_linearization_point,
@@ -755,7 +755,7 @@ namespace aspect
   {
     std::string timer_section_name = "Assemble Stokes system";
 
-    if (assemble_defect_correction_stokes_system)
+    if (assemble_newton_stokes_system)
       {
         if (!rebuild_stokes_matrix && !stokes_matrix_free)
           timer_section_name += " rhs";

source/simulator/core.cc

@@ -248,11 +248,7 @@ namespace aspect
 
     rebuild_stokes_matrix (true),
     rebuild_stokes_preconditioner (true),
-    // note that assemble_defect_correction_stokes_system is
-    // initialized to false even for defect correction solvers
-    // because the first nonlinear iteration is always
-    // solved as the full system
-    assemble_defect_correction_stokes_system (false)
+    assemble_newton_stokes_system (false)
   {
     wall_timer.start();
 
@@ -438,11 +434,11 @@ namespace aspect
         melt_handler->initialize();
       }
 
-
     // Initialize the Newton handler for all defect correction solver schemes
     // (defect correction Picard and Newton solvers)
     if (Parameters<dim>::is_defect_correction(parameters.nonlinear_solver))
       {
+        assemble_newton_stokes_system = true
         newton_handler->initialize_simulator(*this);
         newton_handler->parameters.parse_parameters(prm);
       }
@@ -647,13 +643,6 @@ namespace aspect
     for (auto &particle_manager : particle_managers)
       particle_manager.backup_particles();
 
-    if (Parameters<dim>::is_defect_correction(parameters.nonlinear_solver))
-      {
-        // Ensure we always begin by solving the normal Stokes system
-        // even for defect correction / Newton solvers
-        assemble_defect_correction_stokes_system = false;
-        set_assemblers();
-      }
 
     // then interpolate the current boundary velocities. copy constraints
     // into current_constraints and then add to current_constraints
@@ -1454,7 +1443,7 @@ namespace aspect
         (introspection.n_components,
          [&] (const dealii::Point<dim> &x) -> Tensor<1,dim>
         {
-          if (!assemble_defect_correction_stokes_system || (assemble_defect_correction_stokes_system && nonlinear_iteration == 0))
+          if (!assemble_newton_stokes_system || (assemble_newton_stokes_system && nonlinear_iteration == 0))
             return boundary_velocity_manager.boundary_velocity(boundary_id, x);
           else
             return Tensor<1,dim>();

source/simulator/helper_functions.cc

@@ -2611,16 +2611,19 @@ namespace aspect
     temp_velocity_linearization_point
     .swap(current_linearization_point.block(introspection.block_indices.velocities));
 
-    // Rebuild the whole system to compute the rhs.
-    const bool assemble_dc_Stokes_backup = assemble_defect_correction_stokes_system;
-    assemble_defect_correction_stokes_system = true;
-    rebuild_stokes_preconditioner = false;
+    // Ensure we assemble the Newton residual, no matter which system we are solving
+    const bool assemble_dc_Stokes_backup = assemble_newton_stokes_system;
+    assemble_newton_stokes_system = true;
+    // Ensure the assemblers are consistent with the system we are assembling
+    if (assemble_dc_Stokes_backup != assemble_newton_stokes_system)
+      set_assemblers();
 
     // Technically we only need the rhs, but we have asserts in place that check if
     // the system is assembled correctly when boundary conditions are prescribed, so we assemble the whole system.
     // TODO: This is a waste of time in the first nonlinear iteration. Check if we can modify the asserts in the
     // assemble_stokes_system() function to only assemble the RHS.
     rebuild_stokes_matrix = boundary_velocity_manager.get_prescribed_boundary_velocity_indicators().size()!=0;
+    rebuild_stokes_preconditioner = false;
 
     compute_current_constraints ();
     assemble_stokes_system();
@@ -2633,8 +2636,13 @@ namespace aspect
     current_linearization_point.block(introspection.block_indices.velocities)
     .swap(temp_velocity_linearization_point);
 
-    // Undo the change to this state variable
-    assemble_defect_correction_stokes_system = assemble_dc_Stokes_backup;
+    // Undo the change to the state of assembly, assemblers, and constraints
+    if (assemble_dc_Stokes_backup != assemble_newton_stokes_system)
+      {
+        assemble_newton_stokes_system = assemble_dc_Stokes_backup;
+        set_assemblers();
+        compute_current_constraints ();
+      }
 
     pcout << "   Initial Newton Stokes residual = " << initial_newton_residual << ", v = " << initial_newton_residual_vel << ", p = " << initial_newton_residual_p << std::endl << std::endl;
     return initial_newton_residual;

source/simulator/solver.cc

@@ -595,15 +595,15 @@ namespace aspect
       {
         outputs = stokes_matrix_free->solve(system_matrix,
                                             system_rhs,
-                                            assemble_defect_correction_stokes_system,
+                                            assemble_newton_stokes_system,
                                             last_pressure_normalization_adjustment,
                                             solution_vector);
       }
     else if (parameters.use_direct_stokes_solver)
       {
         outputs = stokes_direct->solve(system_matrix,
                                        system_rhs,
-                                       assemble_defect_correction_stokes_system,
+                                       assemble_newton_stokes_system,
                                        last_pressure_normalization_adjustment,
                                        solution_vector);
       }
@@ -676,7 +676,7 @@ namespace aspect
         // linearized_stokes_variables has a different
         // layout than current_linearization_point, which also contains all the
         // other solution variables.
-        if (assemble_defect_correction_stokes_system == false)
+        if (assemble_newton_stokes_system == false)
           {
             linearized_stokes_initial_guess.block (velocity_block_index) = current_linearization_point.block (velocity_block_index);
             linearized_stokes_initial_guess.block (pressure_block_index) = current_linearization_point.block (pressure_block_index);
@@ -702,7 +702,7 @@ namespace aspect
         linearized_stokes_initial_guess.block (pressure_block_index) /= pressure_scaling;
 
         double solver_tolerance = 0;
-        if (assemble_defect_correction_stokes_system == false)
+        if (assemble_newton_stokes_system == false)
           {
             // (ab)use the distributed solution vector to temporarily put a residual in
             // (we don't care about the residual vector -- all we care about is the
@@ -941,7 +941,7 @@ namespace aspect
         // do some cleanup now that we have the solution
         remove_nullspace(solution_vector, distributed_stokes_solution);
 
-        if (assemble_defect_correction_stokes_system == false)
+        if (assemble_newton_stokes_system == false)
           outputs.pressure_normalization_adjustment = normalize_pressure(solution_vector);
       }
 

source/simulator/solver_schemes.cc

@@ -611,13 +611,18 @@ namespace aspect
         dcr.switch_initial_residual = dcr.initial_residual;
         dcr.residual_old = dcr.initial_residual;
         dcr.residual = dcr.initial_residual;
+
+        // Ensure we always begin by solving the normal Stokes system
+        // even for defect correction / Newton solvers
+        assemble_newton_stokes_system = false;
+        set_assemblers();
+        compute_current_constraints ();
       }
     else if (nonlinear_iteration == 1)
       {
         // Switch to solving the defect correction system
-        assemble_defect_correction_stokes_system = true;
-
-        // Replace the assemblers and constraints
+        // and use the correct assemblers and constraints
+        assemble_newton_stokes_system = true;
         set_assemblers();
         compute_current_constraints ();
       }