using PECInsulator->IsESet() for matrix BCs.

Author: JustinRayAngus

Source/FieldSolver/ImplicitSolvers/ImplicitSolver.H

@@ -251,7 +251,7 @@ protected:
     /**
      * \brief Convert from WarpX FieldBoundaryType to amrex::LinOpBCType
      */
-    [[nodiscard]] amrex::Array<amrex::LinOpBCType,AMREX_SPACEDIM> convertFieldBCToLinOpBC ( const amrex::Array<FieldBoundaryType,AMREX_SPACEDIM>& ) const;
+    [[nodiscard]] amrex::Array<amrex::LinOpBCType,AMREX_SPACEDIM> convertFieldBCToLinOpBC ( const amrex::Array<FieldBoundaryType,AMREX_SPACEDIM>&, const int bdry_side ) const;
 
 };
 

Source/FieldSolver/ImplicitSolvers/ImplicitSolver.cpp

@@ -52,15 +52,15 @@ const Array<FieldBoundaryType,AMREX_SPACEDIM>& ImplicitSolver::GetFieldBoundaryH
 
 Array<LinOpBCType,AMREX_SPACEDIM> ImplicitSolver::GetLinOpBCLo () const
 {
-    return convertFieldBCToLinOpBC(m_WarpX->GetFieldBoundaryLo());
+    return convertFieldBCToLinOpBC(m_WarpX->GetFieldBoundaryLo(),0);
 }
 
 Array<LinOpBCType,AMREX_SPACEDIM> ImplicitSolver::GetLinOpBCHi () const
 {
-    return convertFieldBCToLinOpBC(m_WarpX->GetFieldBoundaryHi());
+    return convertFieldBCToLinOpBC(m_WarpX->GetFieldBoundaryHi(),1);
 }
 
-Array<LinOpBCType,AMREX_SPACEDIM> ImplicitSolver::convertFieldBCToLinOpBC (const Array<FieldBoundaryType,AMREX_SPACEDIM>& a_fbc) const
+Array<LinOpBCType,AMREX_SPACEDIM> ImplicitSolver::convertFieldBCToLinOpBC (const Array<FieldBoundaryType,AMREX_SPACEDIM>& a_fbc, const int bdry_side) const
 {
     Array<LinOpBCType, AMREX_SPACEDIM> lbc;
     for (auto& bc : lbc) { bc = LinOpBCType::interior; }
@@ -82,10 +82,15 @@ Array<LinOpBCType,AMREX_SPACEDIM> ImplicitSolver::convertFieldBCToLinOpBC (const
             // Also for FieldBoundaryType::PMC
             lbc[i] = LinOpBCType::symmetry;
         } else if (a_fbc[i] == FieldBoundaryType::PECInsulator) {
-            ablastr::warn_manager::WMRecordWarning("Implicit solver",
-                "With PECInsulator, in the Curl-Curl preconditioner Neumann boundary will be used since the full boundary is not yet implemented.",
-                ablastr::warn_manager::WarnPriority::medium);
-            lbc[i] = LinOpBCType::symmetry;
+            const int voltage_driven = m_WarpX->GetPECInsulator_IsESet(i,bdry_side);
+            if (voltage_driven) { // Dirichlet for E
+                lbc[i] = LinOpBCType::Dirichlet;
+            } else { // Dirichlet for B
+                ablastr::warn_manager::WMRecordWarning("Implicit solver with current-driven PECInsulator",
+                    "in the Curl-Curl preconditioner. Neumann boundary will be used since the full boundary is not yet implemented.",
+                    ablastr::warn_manager::WarnPriority::medium);
+                lbc[i] = LinOpBCType::symmetry;
+            }
         } else if (a_fbc[i] == FieldBoundaryType::None) {
             WARPX_ABORT_WITH_MESSAGE("LinOpBCType not set for this FieldBoundaryType");
         } else if (a_fbc[i] == FieldBoundaryType::Open) {

Source/FieldSolver/ImplicitSolvers/ThetaImplicitEM.cpp

@@ -302,8 +302,15 @@ void ThetaImplicitEM::InitializeCurlCurlBCMasks ()
                     val1 = 1.0;
                 }
                 if (bc_type == FieldBoundaryType::PECInsulator) {
-                    val0 = 0.5;
-                    val1 = 1.0;
+                    const int voltage_driven = m_WarpX->GetPECInsulator_IsESet(bdry_dir,bdry_side);
+                    if (voltage_driven) { // Dirichlet boundary for E
+                        val0 = 0.0;
+                        val1 = 0.0;
+                    }
+                    else { // Dirichlet boundary for B
+                        val0 = 0.5;
+                        val1 = 1.0;
+                    }
                 }
 
                 // Set mask values on the boundary cells (same for Ex and Ey for 1D_Z)
@@ -381,9 +388,17 @@ void ThetaImplicitEM::InitializeCurlCurlBCMasks ()
                         val2 = 2.0;
                     }
                     if (bc_type == FieldBoundaryType::PECInsulator) {
-                        val0 = 0.5;
-                        val1 = 1.0;
-                        val2 = 1.0;
+                        const int voltage_driven = m_WarpX->GetPECInsulator_IsESet(bdry_dir,bdry_side);
+                        if (voltage_driven) { // Dirichlet boundary for E
+                            val0 = 0.0;
+                            val1 = 0.0;
+                            val2 = 0.0;
+                        }
+                        else { // Dirichlet boundary for B
+                            val0 = 0.5;
+                            val1 = 1.0;
+                            val2 = 1.0;
+                        }
                     }
 
                     // Set mask values on the boundary cells

Source/WarpX.H

@@ -122,6 +122,9 @@ public:
         return field_boundary_hi;
     }
 
+    [[nodiscard]] int GetPECInsulator_IsESet ( const int  bdry_dir,
+                                               const int  bdry_side ) const;
+
     void InitData ();
 
     void Evolve (int numsteps = -1);

Source/WarpX.cpp

@@ -1914,6 +1914,13 @@ WarpX::ReadParameters ()
     }
 }
 
+int WarpX::GetPECInsulator_IsESet ( const int  bdry_dir,
+                                    const int  bdry_side ) const
+{
+    return pec_insulator_boundary->IsESet(bdry_dir,bdry_side);
+}
+
+
 void
 WarpX::BackwardCompatibility ()
 {