Merge pull request #6685 from bangerth/wrap-solvers

When solvers fail, make sure we get an exception message only on process zero.

Author: gassmoeller

source/mesh_deformation/diffusion.cc

@@ -373,7 +373,21 @@ namespace aspect
       this->get_pcout() << "   Solving mesh surface diffusion" << std::endl;
       SolverControl solver_control(5*system_rhs.size(), this->get_parameters().linear_stokes_solver_tolerance*system_rhs.l2_norm());
       SolverCG<LinearAlgebra::Vector> cg(solver_control);
-      cg.solve (matrix, solution, system_rhs, preconditioner_mass);
+      try
+        {
+          cg.solve (matrix, solution, system_rhs, preconditioner_mass);
+        }
+      catch (const std::exception &exc)
+        {
+          // if the solver fails, report the error from processor 0 with some additional
+          // information about its location, and throw a quiet exception on all other
+          // processors
+          Utilities::throw_linear_solver_failure_exception("iterative diffusion surface deformation solver",
+                                                           "MeshDeformation::Diffusion::diffuse_boundary()",
+                                                           std::vector<SolverControl> {solver_control},
+                                                           exc,
+                                                           this->get_mpi_communicator());
+        }
 
       // Distribute constraints on mass matrix
       matrix_constraints.distribute (solution);

source/mesh_deformation/free_surface.cc

@@ -201,7 +201,22 @@ namespace aspect
 
       SolverControl solver_control(5*rhs.size(), this->get_parameters().linear_stokes_solver_tolerance*rhs.l2_norm());
       SolverCG<LinearAlgebra::Vector> cg(solver_control);
-      cg.solve (mass_matrix, dist_solution, rhs, preconditioner_mass);
+
+      try
+        {
+          cg.solve (mass_matrix, dist_solution, rhs, preconditioner_mass);
+        }
+      catch (const std::exception &exc)
+        {
+          // if the solver fails, report the error from processor 0 with some additional
+          // information about its location, and throw a quiet exception on all other
+          // processors
+          Utilities::throw_linear_solver_failure_exception("iterative free surface solver",
+                                                           "MeshDeformation::FreeSurface::project_velocity_onto_boundary()",
+                                                           std::vector<SolverControl> {solver_control},
+                                                           exc,
+                                                           this->get_mpi_communicator());
+        }
 
       mass_matrix_constraints.distribute (dist_solution);
       output = dist_solution;

source/mesh_deformation/interface.cc

@@ -946,8 +946,22 @@ namespace aspect
       SolverControl solver_control(5*rhs.size(), tolerance * rhs.l2_norm());
       SolverCG<LinearAlgebra::Vector> cg(solver_control);
 
-      cg.solve (mesh_matrix, solution, rhs, preconditioner_stiffness);
-      this->get_pcout() << "   Solving mesh displacement system... " << solver_control.last_step() <<" iterations."<< std::endl;
+      try
+        {
+          cg.solve (mesh_matrix, solution, rhs, preconditioner_stiffness);
+          this->get_pcout() << "   Solving mesh displacement system... " << solver_control.last_step() <<" iterations."<< std::endl;
+        }
+      catch (const std::exception &exc)
+        {
+          // if the solver fails, report the error from processor 0 with some additional
+          // information about its location, and throw a quiet exception on all other
+          // processors
+          Utilities::throw_linear_solver_failure_exception("iterative mesh displacement solver",
+                                                           "MeshDeformationHandler::compute_mesh_displacements()",
+                                                           std::vector<SolverControl> {solver_control},
+                                                           exc,
+                                                           this->get_mpi_communicator());
+        }
 
       mesh_velocity_constraints.distribute (solution);
 
@@ -1222,8 +1236,23 @@ namespace aspect
       SolverCG<dealii::LinearAlgebra::distributed::Vector<double>> cg(solver_control_mf);
 
       mesh_velocity_constraints.set_zero(solution);
-      cg.solve(laplace_operator, solution, rhs, preconditioner);
-      this->get_pcout() << "   Solving mesh displacement system... " << solver_control_mf.last_step() <<" iterations."<< std::endl;
+
+      try
+        {
+          cg.solve(laplace_operator, solution, rhs, preconditioner);
+          this->get_pcout() << "   Solving mesh displacement system... " << solver_control_mf.last_step() <<" iterations."<< std::endl;
+        }
+      catch (const std::exception &exc)
+        {
+          // if the solver fails, report the error from processor 0 with some additional
+          // information about its location, and throw a quiet exception on all other
+          // processors
+          Utilities::throw_linear_solver_failure_exception("iterative mesh displacement solver",
+                                                           "MeshDeformationHandler::compute_mesh_displacements_gmg()",
+                                                           std::vector<SolverControl> {solver_control_mf},
+                                                           exc,
+                                                           this->get_mpi_communicator());
+        }
 
       mesh_velocity_constraints.distribute(solution);
       solution.update_ghost_values();

source/simulator/melt.cc

@@ -1262,18 +1262,31 @@ namespace aspect
       Amg_data.aggregation_threshold = 0.02;
       preconditioner.initialize(system_matrix.block(block_idx, block_idx));
 
+      this->get_pcout() << "   Solving fluid velocity system... " << std::flush;
+
       SolverControl solver_control(system_rhs.block(block_idx).size(),
                                    1e-8*system_rhs.block(block_idx).l2_norm());
       SolverCG<LinearAlgebra::Vector> cg(solver_control);
 
-      this->get_pcout() << "   Solving fluid velocity system... " << std::flush;
-
-      cg.solve (system_matrix.block(block_idx, block_idx),
-                distributed_solution.block(block_idx),
-                system_rhs.block(block_idx),
-                preconditioner);
-
-      this->get_pcout() << solver_control.last_step() <<" iterations."<< std::endl;
+      try
+        {
+          cg.solve (system_matrix.block(block_idx, block_idx),
+                    distributed_solution.block(block_idx),
+                    system_rhs.block(block_idx),
+                    preconditioner);
+          this->get_pcout() << solver_control.last_step() <<" iterations."<< std::endl;
+        }
+      catch (const std::exception &exc)
+        {
+          // if the solver fails, report the error from processor 0 with some additional
+          // information about its location, and throw a quiet exception on all other
+          // processors
+          Utilities::throw_linear_solver_failure_exception("iterative melt solver",
+                                                           "MeltHandler::compute_melt_variables",
+                                                           std::vector<SolverControl> {solver_control},
+                                                           exc,
+                                                           this->get_mpi_communicator());
+        }
 
       this->get_current_constraints().distribute (distributed_solution);
       solution.block(block_idx) = distributed_solution.block(block_idx);