Silver Mueller in RZ (#1804)

* Added RZ Er, Et, Ez and Br, Bt and Bz

* Answered my own question

* Making Yee and cells domain check for RZ too

* Added relevant coefficients

* Derivative in z requires coefsz and n_coefsz

* Bt computed with UpwardDz

* UpwardDz might require CylindricalYeeAlgorithm.H

* Added mode 0 Bz

* Added higher-order modes for Bt

* Added higher-order modes for Bz

* Fix to EOLs

* Fix typo

* Added cylindrical specific parameters

* Fix error of #endif

* rmin also needed in RZ parameters

* T_Algo needed for RZ -> different initialization in .H

* Fix private / public function

* Replacing T_Algo by CylindricalYeeAlgorithm

* Fix typo

* ParallelFor for Br, Bt and Bz separated

* Compiled after removing unecessary Br

* Changes suggested by reviewer

* No need to compute r before if()

* Corrected real and imaginary parts of Bz

* Remove vscode file

Co-authored-by: Remi Lehe <remi.lehe@normalesup.org>

1 files changed, 81 insertions, 4 deletions

diff --git a/Source/FieldSolver/FiniteDifferenceSolver/ApplySilverMuellerBoundary.cpp b/Source/FieldSolver/FiniteDifferenceSolver/ApplySilverMuellerBoundary.cpp
index abceead1f..21abb27e4 100644
--- a/Source/FieldSolver/FiniteDifferenceSolver/ApplySilverMuellerBoundary.cpp
+++ b/Source/FieldSolver/FiniteDifferenceSolver/ApplySilverMuellerBoundary.cpp
@@ -9,6 +9,9 @@
 #include "FiniteDifferenceSolver.H"
 #include "Utils/WarpXConst.H"
 #include <AMReX_Gpu.H>
+#ifdef WARPX_DIM_RZ
+#   include "FiniteDifferenceAlgorithms/CylindricalYeeAlgorithm.H"
+#endif
 
 using namespace amrex;
 
@@ -21,10 +24,6 @@ void FiniteDifferenceSolver::ApplySilverMuellerBoundary (
     amrex::Box domain_box,
     amrex::Real const dt ) {
 
-#ifdef WARPX_DIM_RZ
-    amrex::Abort("The Silver-Mueller boundary conditions cannot be used in RZ geometry.");
-#else
-
     // Ensure that we are using the Yee solver
     if (m_fdtd_algo != MaxwellSolverAlgo::Yee) {
         amrex::Abort("The Silver-Mueller boundary conditions can only be used with the Yee solver.");
@@ -33,6 +32,84 @@ void FiniteDifferenceSolver::ApplySilverMuellerBoundary (
     // Ensure that we are using the cells the domain
     domain_box.enclosedCells();
 
+#ifdef WARPX_DIM_RZ
+    // Calculate relevant coefficients
+    amrex::Real const cdt = PhysConst::c*dt;
+    amrex::Real const cdt_over_dr = cdt*m_stencil_coefs_r[0];
+    amrex::Real const coef1_r = (1._rt - cdt_over_dr)/(1._rt + cdt_over_dr);
+    amrex::Real const coef2_r = 2._rt*cdt_over_dr/(1._rt + cdt_over_dr) / PhysConst::c;
+    amrex::Real const coef3_r = cdt/(1._rt + cdt_over_dr) / PhysConst::c;
+
+    // Extract stencil coefficients
+    Real const * const AMREX_RESTRICT coefs_z = m_stencil_coefs_z.dataPtr();
+    int const n_coefs_z = m_stencil_coefs_z.size();
+
+    // Extract cylindrical specific parameters
+    Real const dr = m_dr;
+    int const nmodes = m_nmodes;
+    Real const rmin = m_rmin;
+
+    // tiling is usually set by TilingIfNotGPU()
+    // but here, we set it to false because of potential race condition,
+    // since we grow the tiles by one guard cell after creating them.
+    for ( MFIter mfi(*Efield[0], false); mfi.isValid(); ++mfi ) {
+        // Extract field data for this grid/tile
+        Array4<Real> const& Er = Efield[0]->array(mfi);
+        Array4<Real> const& Et = Efield[1]->array(mfi);
+        Array4<Real> const& Ez = Efield[2]->array(mfi);
+        Array4<Real> const& Bt = Bfield[1]->array(mfi);
+        Array4<Real> const& Bz = Bfield[2]->array(mfi);
+
+        // Extract tileboxes for which to loop
+        Box tbt  = mfi.tilebox(Bfield[1]->ixType().toIntVect());
+        Box tbz  = mfi.tilebox(Bfield[2]->ixType().toIntVect());
+
+        // We will modify the first (i.e. innermost) guard cell
+        // (if it is outside of the physical domain)
+        // Thus, the tileboxes here are grown by 1 guard cell
+        tbt.grow(1);
+        tbz.grow(1);
+
+        // Loop over the cells
+        amrex::ParallelFor(tbt, tbz,
+            [=] AMREX_GPU_DEVICE (int i, int j, int /*k*/){
+
+                // At the +r boundary (innermost guard cell)
+                if ( i==domain_box.bigEnd(0)+1 ){
+                    // Mode 0
+                    Bt(i,j,0,0) = coef1_r*Bt(i,j,0,0) - coef2_r*Ez(i,j,0,0)
+                        + coef3_r*CylindricalYeeAlgorithm::UpwardDz(Er, coefs_z, n_coefs_z, i, j, 0, 0);
+                    for (int m=1; m<nmodes; m++) { // Higher-order modes
+                        // Real part
+                        Bt(i,j,0,2*m-1) = coef1_r*Bt(i,j,0,2*m-1) - coef2_r*Ez(i,j,0,2*m-1)
+                            + coef3_r*CylindricalYeeAlgorithm::UpwardDz(Er, coefs_z, n_coefs_z, i, j, 0, 2*m-1);
+                        // Imaginary part
+                        Bt(i,j,0,2*m) = coef1_r*Bt(i,j,0,2*m) - coef2_r*Ez(i,j,0,2*m)
+                            + coef3_r*CylindricalYeeAlgorithm::UpwardDz(Er, coefs_z, n_coefs_z, i, j, 0, 2*m);
+                    }
+                }
+            },
+            [=] AMREX_GPU_DEVICE (int i, int j, int /*k*/){
+
+                // At the +r boundary (innermost guard cell)
+                if ( i==domain_box.bigEnd(0)+1 ){
+                    Real const r = rmin + (i + 0.5_rt)*dr; // r on nodal point (Bz is cell-centered in r)
+                    // Mode 0
+                    Bz(i,j,0,0) = coef1_r*Bz(i,j,0,0) + coef2_r*Et(i,j,0,0) - coef3_r*Et(i,j,0,0)/r;
+                    for (int m=1; m<nmodes; m++) { // Higher-order modes
+                        // Real part
+                        Bz(i,j,0,2*m-1) = coef1_r*Bz(i,j,0,2*m-1) + coef2_r*Et(i,j,0,2*m-1)
+                            - coef3_r/r*(Et(i,j,0,2*m-1) - m*Er(i,j,0,2*m));
+                        // Imaginary part
+                        Bz(i,j,0,2*m) = coef1_r*Bz(i,j,0,2*m) + coef2_r*Et(i,j,0,2*m)
+                            - coef3_r/r*(Et(i,j,0,2*m) + m*Er(i,j,0,2*m-1));
+                    }
+                }
+            }
+        );
+    }
+#else
+
     // Calculate relevant coefficients
     amrex::Real const cdt_over_dx = PhysConst::c*dt*m_stencil_coefs_x[0];
     amrex::Real const coef1_x = (1._rt - cdt_over_dx)/(1._rt + cdt_over_dx);