Implicit Field Solve Preconditioner based on Curl-Curl Operator

Regression/Checksum/benchmarks_json/test_3d_open_bc_poisson_solver.json

@@ -17,4 +17,4 @@
     "particle_position_z": 10066.329600000008,
     "particle_weight": 19969036501.910976
   }
-}
+}

Source/FieldSolver/ImplicitSolvers/CMakeLists.txt

@@ -2,6 +2,7 @@ foreach(D IN LISTS WarpX_DIMS)
     warpx_set_suffix_dims(SD ${D})
     target_sources(lib_${SD}
       PRIVATE
+        ImplicitSolver.cpp
         SemiImplicitEM.cpp
         ThetaImplicitEM.cpp
         WarpXImplicitOps.cpp

Source/FieldSolver/ImplicitSolvers/ImplicitSolver.H

@@ -1,4 +1,4 @@
-/* Copyright 2024 Justin Angus
+/* Copyright 2024 Justin Angus, Debojyoti Ghosh
  *
  * This file is part of WarpX.
  *
@@ -7,11 +7,13 @@
 #ifndef Implicit_Solver_H_
 #define Implicit_Solver_H_
 
-#include "FieldSolver/ImplicitSolvers/WarpXSolverVec.H"
-#include "NonlinearSolvers/NonlinearSolverLibrary.H"
-
 #include <AMReX_Array.H>
 #include <AMReX_REAL.H>
+#include <AMReX_LO_BCTYPES.H>
+
+#include "FieldSolver/ImplicitSolvers/WarpXSolverVec.H"
+#include "NonlinearSolvers/NonlinearSolverLibrary.H"
+#include "Utils/WarpXAlgorithmSelection.H"
 
 /**
  * \brief Base class for implicit time solvers. The base functions are those
@@ -85,6 +87,16 @@ public:
                               int              a_nl_iter,
                               bool             a_from_jacobian ) = 0;
 
+    [[nodiscard]] virtual amrex::Real theta () const { return 1.0; }
+
+    [[nodiscard]] int numAMRLevels () const { return m_num_amr_levels; }
+
+    [[nodiscard]] const amrex::Geometry& GetGeometry (int) const;
+    [[nodiscard]] const amrex::Array<FieldBoundaryType,AMREX_SPACEDIM>& GetFieldBoundaryLo () const;
+    [[nodiscard]] const amrex::Array<FieldBoundaryType,AMREX_SPACEDIM>& GetFieldBoundaryHi () const;
+    [[nodiscard]] amrex::Array<amrex::LinOpBCType,AMREX_SPACEDIM> GetLinOpBCLo () const;
+    [[nodiscard]] amrex::Array<amrex::LinOpBCType,AMREX_SPACEDIM> GetLinOpBCHi () const;
+
 protected:
 
     /**
@@ -94,6 +106,11 @@ protected:
 
     bool m_is_defined = false;
 
+    /**
+     * \brief Number of AMR levels
+     */
+    int m_num_amr_levels = 1;
+
     /**
      * \brief Nonlinear solver type and object
      */
@@ -140,6 +157,11 @@ protected:
 
     }
 
+    /**
+     * \brief Convert from WarpX FieldBoundaryType to amrex::LinOpBCType
+     */
+    [[nodiscard]] amrex::Array<amrex::LinOpBCType,AMREX_SPACEDIM> convertFieldBCToLinOpBC ( const amrex::Array<FieldBoundaryType,AMREX_SPACEDIM>& ) const;
+
 };
 
 #endif

Source/FieldSolver/ImplicitSolvers/ImplicitSolver.cpp

@@ -0,0 +1,58 @@
+#include "ImplicitSolver.H"
+#include "WarpX.H"
+
+using namespace amrex;
+
+const Geometry& ImplicitSolver::GetGeometry (const int a_lvl) const
+{
+    AMREX_ASSERT((a_lvl >= 0) && (a_lvl < m_num_amr_levels));
+    return m_WarpX->GetGeometry(a_lvl);
+}
+
+const Array<FieldBoundaryType,AMREX_SPACEDIM>& ImplicitSolver::GetFieldBoundaryLo () const
+{
+    return m_WarpX->GetFieldBoundaryLo();
+}
+
+const Array<FieldBoundaryType,AMREX_SPACEDIM>& ImplicitSolver::GetFieldBoundaryHi () const
+{
+    return m_WarpX->GetFieldBoundaryHi();
+}
+
+Array<LinOpBCType,AMREX_SPACEDIM> ImplicitSolver::GetLinOpBCLo () const
+{
+    return convertFieldBCToLinOpBC(m_WarpX->GetFieldBoundaryLo());
+}
+
+Array<LinOpBCType,AMREX_SPACEDIM> ImplicitSolver::GetLinOpBCHi () const
+{
+    return convertFieldBCToLinOpBC(m_WarpX->GetFieldBoundaryHi());
+}
+
+Array<LinOpBCType,AMREX_SPACEDIM> ImplicitSolver::convertFieldBCToLinOpBC (const Array<FieldBoundaryType,AMREX_SPACEDIM>& a_fbc) const
+{
+    Array<LinOpBCType, AMREX_SPACEDIM> lbc;
+    for (auto& bc : lbc) { bc = LinOpBCType::interior; }
+    for (int i = 0; i < AMREX_SPACEDIM; i++) {
+        if (a_fbc[i] == FieldBoundaryType::PML) {
+            WARPX_ABORT_WITH_MESSAGE("LinOpBCType not set for this FieldBoundaryType");
+        } else if (a_fbc[i] == FieldBoundaryType::Periodic) {
+            lbc[i] = LinOpBCType::Periodic;
+        } else if (a_fbc[i] == FieldBoundaryType::PEC) {
+            WARPX_ABORT_WITH_MESSAGE("LinOpBCType not set for this FieldBoundaryType");
+        } else if (a_fbc[i] == FieldBoundaryType::PMC) {
+            WARPX_ABORT_WITH_MESSAGE("LinOpBCType not set for this FieldBoundaryType");
+        } else if (a_fbc[i] == FieldBoundaryType::Damped) {
+            WARPX_ABORT_WITH_MESSAGE("LinOpBCType not set for this FieldBoundaryType");
+        } else if (a_fbc[i] == FieldBoundaryType::Absorbing_SilverMueller) {
+            WARPX_ABORT_WITH_MESSAGE("LinOpBCType not set for this FieldBoundaryType");
+        } else if (a_fbc[i] == FieldBoundaryType::Neumann) {
+            lbc[i] = LinOpBCType::Neumann;
+        } else if (a_fbc[i] == FieldBoundaryType::None) {
+            WARPX_ABORT_WITH_MESSAGE("LinOpBCType not set for this FieldBoundaryType");
+        } else if (a_fbc[i] == FieldBoundaryType::Open) {
+            WARPX_ABORT_WITH_MESSAGE("LinOpBCType not set for this FieldBoundaryType");
+        }
+    }
+    return lbc;
+}

Source/FieldSolver/ImplicitSolvers/Make.package

@@ -1,3 +1,4 @@
+CEXE_sources += ImplicitSolver.cpp
 CEXE_sources += SemiImplicitEM.cpp
 CEXE_sources += ThetaImplicitEM.cpp
 CEXE_sources += WarpXImplicitOps.cpp

Source/FieldSolver/ImplicitSolvers/ThetaImplicitEM.H

@@ -79,7 +79,7 @@ public:
                       int              a_nl_iter,
                       bool             a_from_jacobian ) override;
 
-    [[nodiscard]] amrex::Real theta () const { return m_theta; }
+    [[nodiscard]] amrex::Real theta () const override { return m_theta; }
 
 private:
 

Source/FieldSolver/ImplicitSolvers/ThetaImplicitEM.cpp

@@ -19,15 +19,15 @@ void ThetaImplicitEM::Define ( WarpX* const  a_WarpX )
 
     // Retain a pointer back to main WarpX class
     m_WarpX = a_WarpX;
+    m_num_amr_levels = 1;
 
     // Define E and Eold vectors
     m_E.Define( m_WarpX, "Efield_fp" );
     m_Eold.Define( m_E );
 
     // Define B_old MultiFabs
     using ablastr::fields::Direction;
-    const int num_levels = 1;
-    for (int lev = 0; lev < num_levels; ++lev) {
+    for (int lev = 0; lev < m_num_amr_levels; ++lev) {
         const auto& ba_Bx = m_WarpX->m_fields.get(FieldType::Bfield_fp, Direction{0}, lev)->boxArray();
         const auto& ba_By = m_WarpX->m_fields.get(FieldType::Bfield_fp, Direction{1}, lev)->boxArray();
         const auto& ba_Bz = m_WarpX->m_fields.get(FieldType::Bfield_fp, Direction{2}, lev)->boxArray();

Source/FieldSolver/ImplicitSolvers/WarpXSolverVec.H

@@ -75,6 +75,7 @@ public:
     void Define ( const WarpXSolverVec&  a_solver_vec )
     {
         assertIsDefined( a_solver_vec );
+        m_num_amr_levels = a_solver_vec.m_num_amr_levels;
         Define( WarpXSolverVec::m_WarpX,
                 a_solver_vec.getVectorType(),
                 a_solver_vec.getScalarType() );
@@ -300,7 +301,7 @@ private:
     std::string m_scalar_type_name = "none";
 
     static constexpr int m_ncomp = 1;
-    static constexpr int m_num_amr_levels = 1;
+    int m_num_amr_levels = 1;
 
     inline static bool m_warpx_ptr_defined = false;
     inline static WarpX* m_WarpX = nullptr;

Source/FieldSolver/ImplicitSolvers/WarpXSolverVec.cpp

@@ -34,6 +34,8 @@ void WarpXSolverVec::Define ( WarpX*  a_WarpX,
         m_warpx_ptr_defined = true;
     }
 
+    m_num_amr_levels = 1;
+
     m_vector_type_name = a_vector_type_name;
     m_scalar_type_name = a_scalar_type_name;