Move PETSc Initialize and Finalize to amrex::Initialize and Finalize

Close #5030.

Author: WeiqunZhang

Src/Base/AMReX.cpp

@@ -29,6 +29,10 @@
 #include <AMReX_Sundials.H>
 #endif
 
+#ifdef AMREX_USE_PETSC
+#include <petscsys.h>
+#endif
+
 #ifdef AMREX_USE_CUPTI
 #include <AMReX_CuptiTrace.H>
 #endif
@@ -135,6 +139,12 @@ namespace {
     std::vector<std::string> command_arguments;
 }
 
+#ifdef AMREX_USE_PETSC
+namespace {
+    bool petsc_initialized_by_amrex = false;
+}
+#endif
+
 namespace {
     std::streamsize  prev_out_precision;
     std::streamsize  prev_err_precision;
@@ -762,6 +772,20 @@ amrex::Initialize (int& argc, char**& argv, bool build_parm_parse,
     }
 #endif
 
+#ifdef AMREX_USE_PETSC
+    {
+        PetscBool petsc_is_initialized = PETSC_FALSE;
+        PetscInitialized(&petsc_is_initialized);
+        if (!petsc_is_initialized) {
+            PetscErrorCode ierr = PetscInitialize(&argc, &argv, nullptr, nullptr);
+            if (ierr != 0) {
+                amrex::Abort("PetscInitialize failed");
+            }
+            petsc_initialized_by_amrex = true;
+        }
+    }
+#endif
+
 #ifdef AMREX_USE_SUNDIALS
     sundials::Initialize(amrex::OpenMP::get_max_threads());
 #endif
@@ -797,6 +821,13 @@ amrex::Finalize (amrex::AMReX* pamrex)
 
     AMReX::erase(pamrex);
 
+#ifdef AMREX_USE_PETSC
+    if (petsc_initialized_by_amrex) {
+        PetscFinalize();
+        petsc_initialized_by_amrex = false;
+    }
+#endif
+
 #ifdef AMREX_USE_HYPRE
     if (init_hypre) { HYPRE_Finalize(); }
 #endif

Src/Extern/PETSc/AMReX_PETSc.cpp

@@ -87,9 +87,6 @@ PETScABecLap::PETScABecLap (const BoxArray& grids, const DistributionMapping& dm
 
     diaginv.define(grids,dmap,ncomp,0);
 
-    PETSC_COMM_WORLD = comm;
-    PetscInitialize(nullptr, nullptr, nullptr, nullptr);
-
     solver = std::make_unique<amrex_KSP>();
     A = std::make_unique<amrex_Mat>();
     b = std::make_unique<amrex_Vec>();
@@ -106,8 +103,6 @@ PETScABecLap::~PETScABecLap ()
     m_factory = nullptr;
     m_bndry = nullptr;
     m_maxorder = -1;
-
-    PetscFinalize();
 }
 
 void
@@ -181,8 +176,8 @@ void
 PETScABecLap::prepareSolver ()
 {
     int num_procs, myid;
-    MPI_Comm_size(PETSC_COMM_WORLD, &num_procs);
-    MPI_Comm_rank(PETSC_COMM_WORLD, &myid);
+    MPI_Comm_size(comm, &num_procs);
+    MPI_Comm_rank(comm, &myid);
 
     const BoxArray& ba = acoefs.boxArray();
     const DistributionMapping& dm = acoefs.DistributionMap();
@@ -341,7 +336,7 @@ PETScABecLap::prepareSolver ()
     Vector<PetscInt> ncells_allprocs(num_procs);
     MPI_Allgather(&ncells_proc, sizeof(PetscInt), MPI_CHAR,
                   ncells_allprocs.data(), sizeof(PetscInt), MPI_CHAR,
-                  PETSC_COMM_WORLD);
+                  comm);
     PetscInt proc_begin = 0;
     for (int i = 0; i < myid; ++i) {
         proc_begin += ncells_allprocs[i];
@@ -401,7 +396,7 @@ PETScABecLap::prepareSolver ()
     // estimated amount of block off diag elements
     PetscInt o_nz  = d_nz / 2;
     if (A->a) { MatDestroy(&A->a); }
-    MatCreate(PETSC_COMM_WORLD, &A->a);
+    MatCreate(comm, &A->a);
     MatSetType(A->a, MATMPIAIJ);
     MatSetSizes(A->a, ncells_proc, ncells_proc, ncells_world, ncells_world);
     MatMPIAIJSetPreallocation(A->a, d_nz, nullptr, o_nz, nullptr );
@@ -589,7 +584,7 @@ PETScABecLap::prepareSolver ()
     MatAssemblyEnd(A->a, MAT_FINAL_ASSEMBLY);
     // create solver
     if (solver->a) { KSPDestroy(&solver->a); }
-    KSPCreate(PETSC_COMM_WORLD, &solver->a);
+    KSPCreate(comm, &solver->a);
     KSPSetOperators(solver->a, A->a, A->a);
 
     // Set up preconditioner
@@ -607,7 +602,7 @@ PETScABecLap::prepareSolver ()
     // create b & x
     if (x->a) { VecDestroy(&x->a); }
     if (b->a) { VecDestroy(&b->a); }
-    VecCreateMPI(PETSC_COMM_WORLD, ncells_proc, ncells_world, &x->a);
+    VecCreateMPI(comm, ncells_proc, ncells_world, &x->a);
     VecDuplicate(x->a, &b->a);
 }