SNES and PETSc solvers call MPI_Abort on exceptions

include/bout/boutcomm.hxx

@@ -4,9 +4,9 @@
 * 
 *
 **************************************************************************
-* Copyright 2010 B.D.Dudson, S.Farley, M.V.Umansky, X.Q.Xu
+* Copyright 2010 - 2026 BOUT++ contributors
 *
-* Contact: Ben Dudson, bd512@york.ac.uk
+* Contact: Ben Dudson, dudson2@llnl.gov
 * 
 * This file is part of BOUT++.
 *
@@ -49,6 +49,8 @@ public:
   static int rank(); ///< Rank: my processor number
   static int size(); ///< Size: number of processors
 
+  static void abort(int errorcode); ///< MPI abort
+
   // Setting options
   void setComm(MPI_Comm c);
 

src/solver/impls/petsc/petsc.cxx

@@ -2,7 +2,7 @@
  * Interface to PETSc solver
  *
  **************************************************************************
- * Copyright 2010 - 2025 BOUT++ contributors
+ * Copyright 2010 - 2026 BOUT++ contributors
  *
  * Contact: Ben Dudson, dudson2@llnl.gov
  *
@@ -942,12 +942,10 @@ PetscErrorCode PetscSolver::rhs(BoutReal t, Vec udata, Vec dudata, bool linear)
   } catch (BoutException& e) {
     // Simulation might fail, e.g. negative densities
     // if timestep too large
-    output_warn.write("WARNING: BoutException thrown: {}\n", e.what());
-
-    // Tell SNES that the input was out of domain
-    SNESSetFunctionDomainError(snes);
-    // Note: Returning non-zero error here leaves vectors in locked state
-    return 0;
+    output_error.write("BoutException thrown: {}\n", e.what());
+    // There is no way to recover and synchronise MPI ranks
+    // unless they all threw an exception at the same point.
+    BoutComm::abort(1);
   }
 
   // Save derivatives to PETSc

src/solver/impls/snes/snes.cxx

@@ -1108,7 +1108,10 @@ int SNESSolver::run() {
           run_rhs(simtime);
         } catch (BoutException& e) {
           output_error.write("ERROR: BoutException thrown: {}\n", e.what());
-          return 1;
+          // Abort simulation. There is no way to recover and
+          // synchronise unless all processors throw exceptions
+          // together
+          BoutComm::abort(1);
         }
 
         // Copy derivatives back
@@ -1262,6 +1265,9 @@ int SNESSolver::run() {
       run_rhs(target); // Run RHS to calculate auxilliary variables
     } catch (BoutException& e) {
       output_error.write("ERROR: BoutException thrown: {}\n", e.what());
+      // Abort simulation. There is no way to recover unless
+      // all processors throw an exception at the same point.
+      BoutComm::abort(1);
     }
 
     if (call_monitors(target, s, getNumberOutputSteps()) != 0) {
@@ -1538,7 +1544,9 @@ PetscErrorCode SNESSolver::rhs_function(Vec x, Vec f, bool linear) {
     // Simulation might fail, e.g. negative densities
     // if timestep too large
     output_warn.write("WARNING: BoutException thrown: {}\n", e.what());
-    return 2;
+    // Abort simulation. Unless all processors threw an exception
+    // at the same point, there is no way to synchronise again.
+    BoutComm::abort(2);
   }
 
   // Copy derivatives back

src/sys/boutcomm.cxx

@@ -58,6 +58,8 @@ int BoutComm::size() {
   return NPES;
 }
 
+void BoutComm::abort(int errorcode) { MPI_Abort(get(), errorcode); }
+
 BoutComm* BoutComm::getInstance() {
   if (instance == nullptr) {
     // Create the singleton object