Boundary Condition : PEC (#1767)

* Read boundary and set periodicity, enumerate BC types, added support for periodic

* separate particle and field boudnary structs

* PEC Etangential Bnormal and setting default pml=0

* eol

* add PEC and FieldBoundary routines

* eol

* remove duplicate struct

* update guard cells for PEC

* host device functions

* remove duplicate function

* fix compilation error by adding ;

* temporary logic to set PEC boundary

* eol

* add documentation for PEC and abort if PEC is used for RZ

* documentation for PEC functions

* remove print statements

* Adding field CI test for PEC for 3D with and without MR

* fix eol

* Apply suggestions from code review

* Update Docs/source/usage/parameters.rst

* Update Source/BoundaryConditions/WarpX_PEC.cpp

* particle PEC CI test

* rename analysis scripts for PEC

* use iside == 0 as condition

* ParallelFor over guard cells. Set field values on boundary and guard in the kernel

* ijk_guard not needed anymore

* remove unnecessary nbody generated by yt

* pi and clight are built-in

* cleanup

* eol

* fix bug in computing ig and ijk_valid

* bug fix in guard-cell update in PEC changes benchmark for particles

* fix bug in the ig compute

* reset particle PEC benchmark

* fix logic to set PEC with current interface

* adding a none type for RZ

* reset benchmarks for 2dLaserInjection, 2dgalilean_hybrid and 2dcomoving CI tests due to differences in PEC reflection

* ensure noz shape factor is used for RZ and k=0 for RZ

* fix XZ/RZ logic for Ey,Etheta,By,Btheta. Add ncomponent for RZ, generically.

* fix the field update for corner cases and mixed boundaries. rename ijk_valid to ijk_mirror since ijk_mirror need not be in valid region

* set GuardCell and Boundary cell flags

* fix compilation bug

* suppress unused variable warning

* eol whitespace

* fix doc

* Apply suggestions from code review

Co-authored-by: Andy Nonaka <AJNonaka@lbl.gov>

* RZ abort message

* reset benchmark after fixing guard-cell updates

* fix typo

* RZ doc

* Update Source/Utils/WarpXUtil.cpp

* Update Examples/Tests/PEC/inputs_particle_PEC_3d

* changes suggested from code-review

* set Enormal and Btangential in guard cell across PEC boundary

* add more info about guard cell update in the docs

* fix typo

* ijk mirror only if ig>0

* update benchmarks

* ignore unused k

* temporarily revert guardcell damping for Cartesian

* fix typo

* dont set Etangential and Bnormal to zero if field is not nodal on boundary

* set bindary int as boolean

* fix high to hi

* set default nox to 1

* particle shape for PEC test

* update benchmarks for 2D comoving and galilean

* reset benchmark for RZ test-cases

* Add in benchmarks generated by @EZoni

* Changing transverse boundaries to periodic for psatd galilean cases

* resetting benchmark for comoving_hybrid_2d and galilean_hybrid_2d

* at-least one guard cell must be filled for PEC

* abort if PEC is used for PSATD

* use new boundary interface to set none for rmin and rmax, since pec does not work for psatd

* Update Source/WarpX.cpp

* add doc mentioning PEC does not work for PSATD

* missing semicolon

* Update Source/BoundaryConditions/WarpX_PEC.H

Co-authored-by: Edoardo Zoni <59625522+EZoni@users.noreply.github.com>

* Update Source/BoundaryConditions/WarpX_PEC.H

* lengthy comment for normal and tangential components at domain boundaries

* eol fix

* grammar

* field_hi

* reset benchmarks

Co-authored-by: Andy Nonaka <AJNonaka@lbl.gov>
Co-authored-by: Edoardo Zoni <59625522+EZoni@users.noreply.github.com>

1 files changed, 26 insertions, 6 deletions

diff --git a/Source/Utils/WarpXUtil.cpp b/Source/Utils/WarpXUtil.cpp
index d85cd17a4..e9cb6e34f 100644
--- a/Source/Utils/WarpXUtil.cpp
+++ b/Source/Utils/WarpXUtil.cpp
@@ -446,6 +446,7 @@ void ReadBCParams ()
     ParmParse pp_geometry("geometry");
     ParmParse pp_warpx("warpx");
     ParmParse pp_algo("algo");
+    int maxwell_solver_id = GetAlgorithmInteger(pp_algo, "maxwell_solver");
     if (pp_geometry.queryarr("is_periodic", geom_periodicity)) {
         // set default field and particle boundary appropriately
         for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
@@ -458,15 +459,24 @@ void ReadBCParams ()
             } else {
                 // if non-periodic and do_pml=0, then set default boundary to PEC
                 int pml_input = 1;
+                int silverMueller_input = 0;
                 pp_warpx.query("do_pml", pml_input);
-                if (pml_input == 0) {
-                    WarpX::field_boundary_lo[idim] = FieldBoundaryType::PEC;
-                    WarpX::field_boundary_hi[idim] = FieldBoundaryType::PEC;
+                pp_warpx.query("do_silver_mueller", silverMueller_input);
+                if (pml_input == 0 and silverMueller_input == 0) {
+                    if (maxwell_solver_id == MaxwellSolverAlgo::PSATD) {
+                        WarpX::field_boundary_lo[idim] = FieldBoundaryType::None;
+                        WarpX::field_boundary_hi[idim] = FieldBoundaryType::None;
+                    } else {
+                        WarpX::field_boundary_lo[idim] = FieldBoundaryType::PEC;
+                        WarpX::field_boundary_hi[idim] = FieldBoundaryType::PEC;
+                    }
+#ifdef WARPX_DIM_RZ
+                    if (idim == 0) WarpX::field_boundary_lo[idim] = FieldBoundaryType::None;
+#endif
                 }
             }
         }
-        // Temporarily setting default boundary to Damped until PEC Boundary Type is enabled
-        int maxwell_solver_id = GetAlgorithmInteger(pp_algo, "maxwell_solver");
+        // Temporarily setting default boundary to Damped until new boundary interface is introduced
         if (maxwell_solver_id == MaxwellSolverAlgo::PSATD) {
             ParmParse pp_psatd("psatd");
             int do_moving_window = 0;
@@ -528,9 +538,19 @@ void ReadBCParams ()
                 WarpX::particle_boundary_lo[idim] = ParticleBoundaryType::Periodic;
                 WarpX::particle_boundary_hi[idim] = ParticleBoundaryType::Periodic;
             }
-
+        }
+        if (maxwell_solver_id == MaxwellSolverAlgo::PSATD) {
+            if (WarpX::field_boundary_lo[idim] == FieldBoundaryType::PEC ||
+                WarpX::field_boundary_hi[idim] == FieldBoundaryType::PEC) {
+                amrex::Abort(" PEC boundary not implemented for PSATD, yet!");
+            }
         }
     }
+#ifdef WARPX_DIM_RZ
+    // Ensure code aborts if PEC is specified at r=0 for RZ
+    AMREX_ALWAYS_ASSERT_WITH_MESSAGE( WarpX::field_boundary_lo[0] == FieldBoundaryType::None,
+        "Error : Field boundary at r=0 must be ``none``. \n");
+#endif
 
     pp_geometry.addarr("is_periodic", geom_periodicity);
 }