Implemented implicit field solve preconditioner based on the curl-cur…

…l MLMG solver
LLNL · Sep 18, 2024 · af654f8 · af654f8
1 parent 55dd436
commit af654f8
Show file tree

Hide file tree

Showing 13 changed files with 659 additions and 22 deletions.
diff --git a/Source/FieldSolver/ImplicitSolvers/CMakeLists.txt b/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

diff --git a/Source/FieldSolver/ImplicitSolvers/ImplicitSolver.H b/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 (const int) const;
+    [[nodiscard]] const amrex::Array<FieldBoundaryType,AMREX_SPACEDIM>& GetFieldBoundaryLo () const;
+    [[nodiscard]] const amrex::Array<FieldBoundaryType,AMREX_SPACEDIM>& GetFieldBoundaryHi () const;
+    [[nodiscard]] const amrex::Array<amrex::LinOpBCType,AMREX_SPACEDIM> GetLinOpBCLo () const;
+    [[nodiscard]] const 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
+     */
+    amrex::Array<amrex::LinOpBCType,AMREX_SPACEDIM> convertFieldBCToLinOpBC ( const amrex::Array<FieldBoundaryType,AMREX_SPACEDIM>& ) const;
+
 };
 
 #endif
diff --git a/Source/FieldSolver/ImplicitSolvers/ImplicitSolver.cpp b/Source/FieldSolver/ImplicitSolvers/ImplicitSolver.cpp
@@ -0,0 +1,57 @@
+#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();
+}
+
+const Array<LinOpBCType,AMREX_SPACEDIM> ImplicitSolver::GetLinOpBCLo () const
+{
+    return convertFieldBCToLinOpBC(m_WarpX->GetFieldBoundaryLo());
+}
+
+const 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 (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;
+}
diff --git a/Source/FieldSolver/ImplicitSolvers/Make.package b/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

diff --git a/Source/FieldSolver/ImplicitSolvers/ThetaImplicitEM.H b/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:
 

diff --git a/Source/FieldSolver/ImplicitSolvers/ThetaImplicitEM.cpp b/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, FieldType::Efield_fp );
     m_Eold.Define( m_E );
 
     // Define Bold MultiFab
-    const int num_levels = 1;
-    m_Bold.resize(num_levels); // size is number of levels
-    for (int lev = 0; lev < num_levels; ++lev) {
+    m_Bold.resize(m_num_amr_levels); // size is number of levels
+    for (int lev = 0; lev < m_num_amr_levels; ++lev) {
         for (int n=0; n<3; n++) {
             const amrex::MultiFab& Bfp = m_WarpX->getField( FieldType::Bfield_fp,lev,n);
             m_Bold[lev][n] = std::make_unique<amrex::MultiFab>( Bfp.boxArray(),

diff --git a/Source/FieldSolver/ImplicitSolvers/WarpXSolverVec.H b/Source/FieldSolver/ImplicitSolvers/WarpXSolverVec.H
@@ -74,6 +74,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.getArrayVecType(),
                 a_solver_vec.getScalarVecType() );
@@ -291,7 +292,7 @@ private:
     warpx::fields::FieldType m_scalar_type = warpx::fields::FieldType::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;

diff --git a/Source/FieldSolver/ImplicitSolvers/WarpXSolverVec.cpp b/Source/FieldSolver/ImplicitSolvers/WarpXSolverVec.cpp
@@ -23,6 +23,8 @@ void WarpXSolverVec::Define ( WarpX*     a_WarpX,
         m_warpx_ptr_defined = true;
     }
 
+    m_num_amr_levels = 1;
+
     m_array_type = a_array_type;
     m_scalar_type = a_scalar_type;