Add PMC boundary conditions (#5628)

Docs/source/theory/boundary_conditions.rst

@@ -301,3 +301,23 @@ the right boundary is reflecting.
 
 .. bibliography::
    :keyprefix: bc-
+
+.. _theory-bc-pmc:
+
+Perfect Magnetic Conductor
+----------------------------
+
+This boundary can be used to model a symmetric surface, where charges and current are
+symmetric across the boundary.
+This is equivalent to the Neumann (zero-derivative) boundary condition.
+For the electromagnetic solve, at PMC, the tangential magnetic field and the normal electric
+field are odd across the boundary and set to 0 on the boundary.
+In the guard-cell region, those fields are set equal and
+opposite to the respective field component in the mirror location across the PMC boundary.
+The other components, the normal magnetic field and tangential electric field, are even
+and set equal to the field component in the mirror location in the domain across the PMC boundary.
+
+The PMC boundary condition also impacts the deposition of charge and current density.
+The charge and current densities deposited into the guard cells are reflected back into
+the domain, adding them to the mirror cells in the domain.
+This represents the charge and current from the virtual symmetric particles in the guard cells.

Docs/source/usage/parameters.rst

@@ -533,6 +533,8 @@ Domain Boundary Conditions
 
     * ``pec``: This option can be used to set a Perfect Electric Conductor at the simulation boundary. Please see the :ref:`PEC theory section <theory-bc-pec>` for more details. Note that PEC boundary is invalid at `r=0` for the RZ solver. Please use ``none`` option. This boundary condition does not work with the spectral solver.
 
+    * ``pmc``: This option can be used to set a Perfect Magnetic Conductor at the simulation boundary. Please see the :ref:`PEC theory section <theory-bc-pmc>` for more details. This is equivalent to ``Neumann``. This boundary condition does not work with the spectral solver.
+
     * ``pec_insulator``: This option specifies a mixed perfect electric conductor and insulator boundary, where some part of the
       boundary is PEC and some is insulator. In the insulator portion, the normal fields are extrapolated and the tangential fields
       are either set to the specified value or extrapolated. The region that is insulator is specified using a spatially dependent expression with the insulator being in the area where the value of the expression is greater than zero.

Examples/Tests/pec/CMakeLists.txt

@@ -41,6 +41,16 @@ add_warpx_test(
     OFF  # dependency
 )
 
+add_warpx_test(
+    test_3d_pmc_field  # name
+    3  # dims
+    2  # nprocs
+    inputs_test_3d_pmc_field  # inputs
+    "analysis_pec.py diags/diag1000134"  # analysis
+    "analysis_default_regression.py --path diags/diag1000134"  # checksum
+    OFF  # dependency
+)
+
 add_warpx_test(
     test_2d_pec_field_insulator_implicit  # name
     2  # dims

Examples/Tests/pec/inputs_test_3d_pmc_field

@@ -0,0 +1,54 @@
+# Set-up to test the PMC Boundary condition for the fields
+# Constructive interference between the incident and reflected wave result in a
+# standing wave.
+
+# max step
+max_step = 134
+
+# number of grid points
+amr.n_cell = 32 32 256
+
+# Maximum allowable size of each subdomain
+amr.max_grid_size = 1024
+amr.blocking_factor = 32
+
+amr.max_level = 0
+
+# Geometry
+geometry.dims = 3
+geometry.prob_lo     =  -8.e-6 -8.e-6 -4.e-6
+geometry.prob_hi     =   8.e-6  8.e-6  4.e-6
+
+# Boundary condition
+boundary.field_lo = periodic periodic pmc
+boundary.field_hi = periodic periodic pmc
+
+warpx.serialize_initial_conditions = 1
+
+# Verbosity
+warpx.verbose = 1
+
+# Algorithms
+algo.current_deposition = esirkepov
+# CFL
+warpx.cfl = 0.9
+
+
+my_constants.z1 = -2.e-6
+my_constants.z2 = 2.e-6
+my_constants.wavelength = 1.e-6
+warpx.E_ext_grid_init_style = parse_E_ext_grid_function
+warpx.Ez_external_grid_function(x,y,z) = "0."
+warpx.Ex_external_grid_function(x,y,z) = "0."
+warpx.Ey_external_grid_function(x,y,z) = "((1.e5*sin(2*pi*(z)/wavelength)) * (z<z2) * (z>z1))"
+
+warpx.B_ext_grid_init_style = parse_B_ext_grid_function
+warpx.Bx_external_grid_function(x,y,z)= "(((-1.e5*sin(2*pi*(z)/wavelength))/clight))*(z<z2) * (z>z1) "
+warpx.By_external_grid_function(x,y,z)= "0."
+warpx.Bz_external_grid_function(x,y,z) = "0."
+
+# Diagnostics
+diagnostics.diags_names = diag1
+diag1.intervals = 134
+diag1.diag_type = Full
+diag1.fields_to_plot = Ey Bx

Regression/Checksum/benchmarks_json/test_3d_magnetostatic_eb.json

@@ -1,21 +1,21 @@
 {
   "lev=0": {
-    "Az": 11.358663326449284,
-    "Bx": 111.55929407644248,
-    "By": 111.55929407644244,
-    "Ex": 31257180402.55472,
-    "Ey": 31257180402.55473,
-    "jz": 1034841325.9848926,
-    "phi": 3143521213.0157924,
-    "rho": 3.449203918900721
+    "Az": 11.358663299932457,
+    "Bx": 111.55929615203162,
+    "By": 111.55929615203165,
+    "Ex": 31463410849.74626,
+    "Ey": 31463410849.746258,
+    "jz": 1034841323.6861029,
+    "phi": 3164328318.15416,
+    "rho": 3.4565836983918676
   },
   "beam": {
-    "particle_momentum_x": 1.3604657334742729e-21,
-    "particle_momentum_y": 1.3604657334742772e-21,
-    "particle_momentum_z": 7.150873450281544e-16,
-    "particle_position_x": 11163.99997371537,
-    "particle_position_y": 11163.999973715368,
-    "particle_position_z": 131662.50031035842,
+    "particle_momentum_x": 1.3829464728617761e-21,
+    "particle_momentum_y": 1.3829464728617792e-21,
+    "particle_momentum_z": 7.150871807235339e-16,
+    "particle_position_x": 11163.99997715059,
+    "particle_position_y": 11163.999977150592,
+    "particle_position_z": 131662.5003102683,
     "particle_weight": 20895107655113.465
   }
-}
+}

Regression/Checksum/benchmarks_json/test_3d_magnetostatic_eb_picmi.json

@@ -1,27 +1,27 @@
 {
+  "lev=0": {
+    "Ax": 1.40889223759456e-05,
+    "Ay": 1.4088922375945606e-05,
+    "Az": 11.423480450267745,
+    "Bx": 112.23826705481486,
+    "By": 112.23826705481484,
+    "Bz": 0.00019199345672949735,
+    "Ex": 31557746267.686367,
+    "Ey": 31557746267.686363,
+    "Ez": 3339526660.3539834,
+    "jx": 1980.6549408566577,
+    "jy": 1980.6549408566577,
+    "jz": 1038931605.1197203,
+    "phi": 3171976204.251914,
+    "rho": 3.4840085919357926
+  },
   "beam": {
-    "particle_momentum_x": 1.3878812158350944e-21,
-    "particle_momentum_y": 1.387881215835094e-21,
-    "particle_momentum_z": 7.150872953138685e-16,
-    "particle_position_x": 11163.999973134894,
-    "particle_position_y": 11163.999973134896,
-    "particle_position_z": 131662.5003103311,
+    "particle_momentum_x": 1.4011190163358655e-21,
+    "particle_momentum_y": 1.401119016335865e-21,
+    "particle_momentum_z": 7.15087179293042e-16,
+    "particle_position_x": 11163.99997543546,
+    "particle_position_y": 11163.999975435456,
+    "particle_position_z": 131662.50031026747,
     "particle_weight": 20895107655113.465
-  },
-  "lev=0": {
-    "Ax": 1.408892468360627e-05,
-    "Ay": 1.4088924683606269e-05,
-    "Az": 11.423480469161868,
-    "Bx": 112.23826555908032,
-    "By": 112.2382655590803,
-    "Bz": 0.00019186770330025167,
-    "Ex": 31418238386.183773,
-    "Ey": 31418238386.183773,
-    "Ez": 3461330433.5026026,
-    "jx": 1961.0003914783667,
-    "jy": 1961.0003914783663,
-    "jz": 1038931606.7573991,
-    "phi": 3157908107.1102533,
-    "rho": 3.46977258905983
   }
-}
+}

Regression/Checksum/benchmarks_json/test_3d_pmc_field.json

@@ -0,0 +1,6 @@
+{
+  "lev=0": {
+    "Bx": 4.1354151621557795,
+    "Ey": 8373879983.480644
+  }
+}

Source/BoundaryConditions/WarpXFieldBoundaries.cpp

@@ -56,36 +56,68 @@ void WarpX::ApplyEfieldBoundary(const int lev, PatchType patch_type, amrex::Real
     if (::isAnyBoundary<FieldBoundaryType::PEC>(field_boundary_lo, field_boundary_hi)) {
         if (patch_type == PatchType::fine) {
             PEC::ApplyPECtoEfield(
-                    {m_fields.get(FieldType::Efield_fp, Direction{0}, lev),
-                     m_fields.get(FieldType::Efield_fp, Direction{1}, lev),
-                     m_fields.get(FieldType::Efield_fp, Direction{2}, lev)},
-                    field_boundary_lo, field_boundary_hi,
+                    m_fields.get_alldirs(FieldType::Efield_fp, lev),
+                    field_boundary_lo, field_boundary_hi, FieldBoundaryType::PEC,
                     get_ng_fieldgather(), Geom(lev),
                     lev, patch_type, ref_ratio);
             if (::isAnyBoundary<FieldBoundaryType::PML>(field_boundary_lo, field_boundary_hi)) {
                 // apply pec on split E-fields in PML region
                 const bool split_pml_field = true;
                 PEC::ApplyPECtoEfield(
                     m_fields.get_alldirs(FieldType::pml_E_fp, lev),
-                    field_boundary_lo, field_boundary_hi,
+                    field_boundary_lo, field_boundary_hi, FieldBoundaryType::PEC,
                     get_ng_fieldgather(), Geom(lev),
                     lev, patch_type, ref_ratio,
                     split_pml_field);
             }
         } else {
             PEC::ApplyPECtoEfield(
-                {m_fields.get(FieldType::Efield_cp,Direction{0},lev),
-                 m_fields.get(FieldType::Efield_cp,Direction{1},lev),
-                 m_fields.get(FieldType::Efield_cp,Direction{2},lev)},
-                field_boundary_lo, field_boundary_hi,
-                get_ng_fieldgather(), Geom(lev),
-                lev, patch_type, ref_ratio);
+                    m_fields.get_alldirs(FieldType::Efield_cp, lev),
+                    field_boundary_lo, field_boundary_hi, FieldBoundaryType::PEC,
+                    get_ng_fieldgather(), Geom(lev),
+                    lev, patch_type, ref_ratio);
             if (::isAnyBoundary<FieldBoundaryType::PML>(field_boundary_lo, field_boundary_hi)) {
                 // apply pec on split E-fields in PML region
                 const bool split_pml_field = true;
                 PEC::ApplyPECtoEfield(
                     m_fields.get_alldirs(FieldType::pml_E_cp, lev),
-                    field_boundary_lo, field_boundary_hi,
+                    field_boundary_lo, field_boundary_hi, FieldBoundaryType::PEC,
+                    get_ng_fieldgather(), Geom(lev),
+                    lev, patch_type, ref_ratio,
+                    split_pml_field);
+            }
+        }
+    }
+
+    if (::isAnyBoundary<FieldBoundaryType::PMC>(field_boundary_lo, field_boundary_hi)) {
+        if (patch_type == PatchType::fine) {
+            PEC::ApplyPECtoBfield(
+                    m_fields.get_alldirs(FieldType::Efield_fp, lev),
+                    field_boundary_lo, field_boundary_hi, FieldBoundaryType::PMC,
+                    get_ng_fieldgather(), Geom(lev),
+                    lev, patch_type, ref_ratio);
+            if (::isAnyBoundary<FieldBoundaryType::PML>(field_boundary_lo, field_boundary_hi)) {
+                // apply pec on split E-fields in PML region
+                const bool split_pml_field = true;
+                PEC::ApplyPECtoBfield(
+                    m_fields.get_alldirs(FieldType::pml_E_fp, lev),
+                    field_boundary_lo, field_boundary_hi, FieldBoundaryType::PMC,
+                    get_ng_fieldgather(), Geom(lev),
+                    lev, patch_type, ref_ratio,
+                    split_pml_field);
+            }
+        } else {
+            PEC::ApplyPECtoBfield(
+                    m_fields.get_alldirs(FieldType::Efield_cp, lev),
+                    field_boundary_lo, field_boundary_hi, FieldBoundaryType::PMC,
+                    get_ng_fieldgather(), Geom(lev),
+                    lev, patch_type, ref_ratio);
+            if (::isAnyBoundary<FieldBoundaryType::PML>(field_boundary_lo, field_boundary_hi)) {
+                // apply pec on split E-fields in PML region
+                const bool split_pml_field = true;
+                PEC::ApplyPECtoBfield(
+                    m_fields.get_alldirs(FieldType::pml_E_cp, lev),
+                    field_boundary_lo, field_boundary_hi, FieldBoundaryType::PMC,
                     get_ng_fieldgather(), Geom(lev),
                     lev, patch_type, ref_ratio,
                     split_pml_field);
@@ -152,19 +184,31 @@ void WarpX::ApplyBfieldBoundary (const int lev, PatchType patch_type, DtType a_d
 
     if (::isAnyBoundary<FieldBoundaryType::PEC>(field_boundary_lo, field_boundary_hi)) {
         if (patch_type == PatchType::fine) {
-            PEC::ApplyPECtoBfield( {
-                m_fields.get(FieldType::Bfield_fp,Direction{0},lev),
-                m_fields.get(FieldType::Bfield_fp,Direction{1},lev),
-                m_fields.get(FieldType::Bfield_fp,Direction{2},lev) },
-                field_boundary_lo, field_boundary_hi,
+            PEC::ApplyPECtoBfield(
+                m_fields.get_alldirs(FieldType::Bfield_fp, lev),
+                field_boundary_lo, field_boundary_hi, FieldBoundaryType::PEC,
                 get_ng_fieldgather(), Geom(lev),
                 lev, patch_type, ref_ratio);
         } else {
-            PEC::ApplyPECtoBfield( {
-                m_fields.get(FieldType::Bfield_cp,Direction{0},lev),
-                m_fields.get(FieldType::Bfield_cp,Direction{1},lev),
-                m_fields.get(FieldType::Bfield_cp,Direction{2},lev) },
-                field_boundary_lo, field_boundary_hi,
+            PEC::ApplyPECtoBfield(
+                m_fields.get_alldirs(FieldType::Bfield_cp, lev),
+                field_boundary_lo, field_boundary_hi, FieldBoundaryType::PEC,
+                get_ng_fieldgather(), Geom(lev),
+                lev, patch_type, ref_ratio);
+        }
+    }
+
+    if (::isAnyBoundary<FieldBoundaryType::PMC>(field_boundary_lo, field_boundary_hi)) {
+        if (patch_type == PatchType::fine) {
+            PEC::ApplyPECtoEfield(
+                m_fields.get_alldirs(FieldType::Bfield_fp, lev),
+                field_boundary_lo, field_boundary_hi, FieldBoundaryType::PMC,
+                get_ng_fieldgather(), Geom(lev),
+                lev, patch_type, ref_ratio);
+        } else {
+            PEC::ApplyPECtoEfield(
+                m_fields.get_alldirs(FieldType::Bfield_cp, lev),
+                field_boundary_lo, field_boundary_hi, FieldBoundaryType::PMC,
                 get_ng_fieldgather(), Geom(lev),
                 lev, patch_type, ref_ratio);
         }
@@ -224,7 +268,8 @@ void WarpX::ApplyRhofieldBoundary (const int lev, MultiFab* rho,
 {
     if (::isAnyBoundary<ParticleBoundaryType::Reflecting>(particle_boundary_lo, particle_boundary_hi) ||
         ::isAnyBoundary<ParticleBoundaryType::Thermal>(particle_boundary_lo, particle_boundary_hi) ||
-        ::isAnyBoundary<FieldBoundaryType::PEC>(field_boundary_lo, field_boundary_hi))
+        ::isAnyBoundary<FieldBoundaryType::PEC>(field_boundary_lo, field_boundary_hi) ||
+        ::isAnyBoundary<FieldBoundaryType::PMC>(field_boundary_lo, field_boundary_hi))
     {
         PEC::ApplyReflectiveBoundarytoRhofield(rho,
             field_boundary_lo, field_boundary_hi,
@@ -239,7 +284,8 @@ void WarpX::ApplyJfieldBoundary (const int lev, amrex::MultiFab* Jx,
 {
     if (::isAnyBoundary<ParticleBoundaryType::Reflecting>(particle_boundary_lo, particle_boundary_hi) ||
         ::isAnyBoundary<ParticleBoundaryType::Thermal>(particle_boundary_lo, particle_boundary_hi) ||
-        ::isAnyBoundary<FieldBoundaryType::PEC>(field_boundary_lo, field_boundary_hi))
+        ::isAnyBoundary<FieldBoundaryType::PEC>(field_boundary_lo, field_boundary_hi) ||
+        ::isAnyBoundary<FieldBoundaryType::PMC>(field_boundary_lo, field_boundary_hi))
     {
         PEC::ApplyReflectiveBoundarytoJfield(Jx, Jy, Jz,
             field_boundary_lo, field_boundary_hi,

Source/BoundaryConditions/WarpX_PEC.H

@@ -33,6 +33,7 @@ namespace PEC {
                         std::array<amrex::MultiFab*, 3> Efield,
                         const amrex::Array<FieldBoundaryType,AMREX_SPACEDIM>& field_boundary_lo,
                         const amrex::Array<FieldBoundaryType,AMREX_SPACEDIM>& field_boundary_hi,
+                        FieldBoundaryType bc_type,
                         const amrex::IntVect& ng_fieldgather, const amrex::Geometry& geom,
                         int lev, PatchType patch_type, const amrex::Vector<amrex::IntVect>& ref_ratios,
                         bool split_pml_field = false);
@@ -54,8 +55,10 @@ namespace PEC {
                         std::array<amrex::MultiFab*, 3> Bfield,
                         const amrex::Array<FieldBoundaryType,AMREX_SPACEDIM>& field_boundary_lo,
                         const amrex::Array<FieldBoundaryType,AMREX_SPACEDIM>& field_boundary_hi,
+                        FieldBoundaryType bc_type,
                         const amrex::IntVect& ng_fieldgather, const amrex::Geometry& geom,
-                        int lev, PatchType patch_type, const amrex::Vector<amrex::IntVect>& ref_ratios);
+                        int lev, PatchType patch_type, const amrex::Vector<amrex::IntVect>& ref_ratios,
+                        bool split_pml_field = false);
 
     /**
      * \brief Reflects charge density deposited over the PEC boundary back into