Fix compilation errors

5 files changed, 80 insertions, 42 deletions

diff --git a/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp b/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp
index c08dee04f..72d3a1135 100644
--- a/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp
+++ b/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp
@@ -1,7 +1,7 @@
 #include "WarpXAlgorithmSelection.H"
 #include "FiniteDifferenceSolver.H"
 #ifdef WARPX_DIM_RZ
-    #include "FiniteDifferenceAlgorithms/YeeAlgorithm.H"
+    #include "FiniteDifferenceAlgorithms/CylindricalYeeAlgorithm.H"
 #else
     #include "FiniteDifferenceAlgorithms/YeeAlgorithm.H"
     #include "FiniteDifferenceAlgorithms/CKCAlgorithm.H"
@@ -14,19 +14,20 @@ void FiniteDifferenceSolver::EvolveB ( VectorField& Bfield,
                                        VectorField const& Efield,
                                        amrex::Real const dt ) {
 
+   // Select algorithm (The choice of algorithm is a runtime option,
+   // but we compile code for each algorithm, using templates)
 #ifdef WARPX_DIM_RZ
-    EvolveBCylindrical( Bfield, Efield, dt );
+    if (m_fdtd_algo == MaxwellSolverAlgo::Yee){
+        EvolveBCylindrical <CylindricalYeeAlgorithm> ( Bfield, Efield, dt );
 #else
-    // Select algorithm (The choice of algorithm is a runtime option,
-    // but we compile code for each algorithm, using templates)
     if (m_fdtd_algo == MaxwellSolverAlgo::Yee){
         EvolveBCartesian <YeeAlgorithm> ( Bfield, Efield, dt );
     } else if (m_fdtd_algo == MaxwellSolverAlgo::CKC) {
         EvolveBCartesian <CKCAlgorithm> ( Bfield, Efield, dt );
+#endif
     } else {
         amrex::Abort("Unknown algorithm");
     }
-#endif
 
 }
 
@@ -110,12 +111,16 @@ void FiniteDifferenceSolver::EvolveBCylindrical ( VectorField& Bfield,
         auto const& Ez = Efield[2]->array(mfi);
 
         // Extract stencil coefficients
-        Real const dr = m_dr;
         Real const* AMREX_RESTRICT coefs_r = stencil_coefs_r.dataPtr();
         int const n_coefs_r = stencil_coefs_r.size();
         Real const* AMREX_RESTRICT coefs_z = stencil_coefs_z.dataPtr();
         int const n_coefs_z = stencil_coefs_z.size();
 
+        // Extract cylindrical specific parameters
+        Real const dr = m_dr;
+        int const nmodes = m_nmodes;
+        Real const rmin = m_rmin;
+
         // Extract tileboxes for which to loop
         const Box& tbr  = mfi.tilebox(Bfield[0]->ixType().ixType());
         const Box& tbt  = mfi.tilebox(Bfield[1]->ixType().ixType());
@@ -128,9 +133,11 @@ void FiniteDifferenceSolver::EvolveBCylindrical ( VectorField& Bfield,
                 Real const r = rmin + i*dr; // r on nodal point (Br is nodal in r)
                 Br(i, j, 0, 0) += dt * T_Algo::UpwardDz(Et, coefs_z, n_coefs_z, i, j, 0, 0); // Mode m=0
                 for (int m=1; m<nmodes; m++) { // Higher-order modes
-                    Br(i, j, 0, 2*m-1) += dt*( T_Algo::UpwardDz(Et, inv_dz, i, j, 0, 2*m-1)
+                    Br(i, j, 0, 2*m-1) += dt*(
+                        T_Algo::UpwardDz(Et, coefs_z, n_coefs_z, i, j, 0, 2*m-1)
                         - m * T_Algo::DivideByR(Ez, r, dr, m, i, j, 0, 2*m  ));  // Real part
-                    Br(i, j, 0, 2*m  ) += dt*( T_Algo::UpwardDz(Et, inv_dz, i, j, 0, 2*m  )
+                    Br(i, j, 0, 2*m  ) += dt*(
+                        T_Algo::UpwardDz(Et, coefs_z, n_coefs_z, i, j, 0, 2*m  )
                         + m * T_Algo::DivideByR(Ez, r, dr, m, i, j, 0, 2*m-1)); // Imaginary part
                 }
                 // Ensure that Br remains 0 on axis (except for m=1)
@@ -144,24 +151,28 @@ void FiniteDifferenceSolver::EvolveBCylindrical ( VectorField& Bfield,
             },
 
             [=] AMREX_GPU_DEVICE (int i, int j, int k){
-                Bt(i, j, 0, 0) += dt*( T_Algo::UpwardDr(Ez, inv_dr, i, j, 0, 0)
-                                     - T_Algo::UpwardDz(Er, inv_dz, i, j, 0, 0)); // Mode m=0
+                Bt(i, j, 0, 0) += dt*(
+                    T_Algo::UpwardDr(Ez, coefs_r, n_coefs_r, i, j, 0, 0)
+                    - T_Algo::UpwardDz(Er, coefs_z, n_coefs_z, i, j, 0, 0)); // Mode m=0
                 for (int m=1 ; m<nmodes ; m++) { // Higher-order modes
-                    Bt(i, j, 0, 2*m-1) += dt*( T_Algo::UpwardDr(Ez, inv_dr, i, j, 0, 2*m-1)
-                                         - T_Algo::UpwardDz(Er, inv_dz, i, j, 0, 2*m-1)); // Real part
-                    Bt(i, j, 0, 2*m  ) += dt*( T_Algo::UpwardDr(Ez, inv_dr, i, j, 0, 2*m  )
-                                         - T_Algo::UpwardDz(Er, inv_dz, i, j, 0, 2*m  )); // Imaginary part
+                    Bt(i, j, 0, 2*m-1) += dt*(
+                        T_Algo::UpwardDr(Ez, coefs_r, n_coefs_r, i, j, 0, 2*m-1)
+                        - T_Algo::UpwardDz(Er, coefs_z, n_coefs_z, i, j, 0, 2*m-1)); // Real part
+                    Bt(i, j, 0, 2*m  ) += dt*(
+                        T_Algo::UpwardDr(Ez, coefs_r, n_coefs_r, i, j, 0, 2*m  )
+                        - T_Algo::UpwardDz(Er, coefs_z, n_coefs_z, i, j, 0, 2*m  )); // Imaginary part
                 }
             },
 
             [=] AMREX_GPU_DEVICE (int i, int j, int k){
                 Real const r = rmin + (i + 0.5)*dr; // r on a cell-centered grid (Bz is cell-centered in r)
-                Bz(i, j, 0, 0) += dt*( - T_Algo::UpwardDrr_over_r(Et, r, dr, i, j, 0, 0));
+                Bz(i, j, 0, 0) += dt*( - T_Algo::UpwardDrr_over_r(Et, r, dr, coefs_r, n_coefs_r, i, j, 0, 0));
                 for (int m=1 ; m<nmodes ; m++) { // Higher-order modes
                     Bz(i, j, 0, 2*m-1) += dt*( m * Er(i, j, 0, 2*m  )/r
-                        - T_Algo::UpwardDrr_over_r(Et, r, dr, i, j, 0, 2*m-1)); // Real part
+                        - T_Algo::UpwardDrr_over_r(Et, r, dr, coefs_r, n_coefs_r, i, j, 0, 2*m-1)); // Real part
                     Bz(i, j, 0, 2*m  ) += dt*(-m * Er(i, j, 0, 2*m-1)/r
-                        - T_Algo::UpwardDrr_over_r(Et, r, dr, i, j, 0, 2*m  )); // Imaginary part
+                        - T_Algo::UpwardDrr_over_r(Et, r, dr, coefs_r, n_coefs_r, i, j, 0, 2*m  )); // Imaginary part
+                }
             }
 
         );
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CylindricalYeeAlgorithm.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CylindricalYeeAlgorithm.H
index 9258fdff4..5d160d865 100644
--- a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CylindricalYeeAlgorithm.H
+++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CylindricalYeeAlgorithm.H
@@ -19,6 +19,21 @@ struct CylindricalYeeAlgorithm {
         stencil_coefs_z[0] = 1./cell_size[2];  // 1./dz
     }
 
+    /** Applies the differential operator `1/r * d(rF)/dr`,
+     * where `F` is on a *nodal* grid in `r`
+     * and the differential operator is evaluated on *cell-centered* grid.
+     * The input parameter `r` is given at the cell-centered position*/
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    static amrex::Real UpwardDrr_over_r(
+        amrex::Array4<amrex::Real> const& F,
+        amrex::Real const r, int const dr,
+        amrex::Real const* coefs_r, int const n_coefs_r,
+        int const i, int const j, int const k, int const comp ) {
+
+        amrex::Real inv_dr = coefs_r[0];
+        return inv_dr*( (r+0.5*dr)*F(i+1,j,k,comp) - (r-0.5*dr)*F(i,j,k,comp) );
+    };
+
     AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
     static amrex::Real UpwardDr(
         amrex::Array4<amrex::Real> const& F,
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H
index deb89b1d6..9ce910e3d 100644
--- a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H
+++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H
@@ -2,7 +2,9 @@
 #define WARPX_FINITE_DIFFERENCE_SOLVER_H_
 
 #include "WarpXAlgorithmSelection.H"
-#ifndef WARPX_DIM_RZ
+#ifdef WARPX_DIM_RZ
+    #include "FiniteDifferenceAlgorithms/CylindricalYeeAlgorithm.H"
+#else
     #include "FiniteDifferenceAlgorithms/YeeAlgorithm.H"
     #include "FiniteDifferenceAlgorithms/CKCAlgorithm.H"
 #endif
@@ -21,29 +23,7 @@ class FiniteDifferenceSolver
 
         // Constructor
         FiniteDifferenceSolver ( int const fdtd_algo,
-            std::array<amrex::Real,3> cell_size ) {
-
-            // Register the type of finite-difference algorithm
-            m_fdtd_algo = fdtd_algo;
-
-            // Calculate coefficients of finite-difference stencil
-#ifdef WARPX_DIM_RZ
-            m_dr = cell_size[0];
-            if (fdtd_algo == MaxwellSolverAlgo::Yee){
-                CylindricalYeeAlgorithm::InitializeStencilCoefficients( cell_size,
-                    stencil_coefs_r, stencil_coefs_z );
-#else
-            if (fdtd_algo == MaxwellSolverAlgo::Yee){
-                YeeAlgorithm::InitializeStencilCoefficients( cell_size,
-                    stencil_coefs_x, stencil_coefs_y, stencil_coefs_z );
-            } else if (fdtd_algo == MaxwellSolverAlgo::CKC) {
-                CKCAlgorithm::InitializeStencilCoefficients( cell_size,
-                    stencil_coefs_x, stencil_coefs_y, stencil_coefs_z );
-#endif
-            } else {
-                amrex::Abort("Unknown algorithm");
-            }
-        };
+            std::array<amrex::Real,3> cell_size );
 
         void EvolveB ( VectorField& Bfield,
                        VectorField const& Efield,
@@ -53,7 +33,8 @@ class FiniteDifferenceSolver
         int m_fdtd_algo;
 
 #ifdef WARPX_DIM_RZ
-        amrex::Real m_dr;
+        amrex::Real m_dr, m_rmin;
+        amrex::Real m_nmodes;
         amrex::Gpu::ManagedVector<amrex::Real> stencil_coefs_r;
         amrex::Gpu::ManagedVector<amrex::Real> stencil_coefs_z;
 #else
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.cpp b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.cpp
new file mode 100644
index 000000000..d89e6e9d3
--- /dev/null
+++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.cpp
@@ -0,0 +1,30 @@
+#include "FiniteDifferenceSolver.H"
+#include "WarpX.H"
+
+// Constructor
+FiniteDifferenceSolver::FiniteDifferenceSolver ( int const fdtd_algo,
+    std::array<amrex::Real,3> cell_size ) {
+
+    // Register the type of finite-difference algorithm
+    m_fdtd_algo = fdtd_algo;
+
+    // Calculate coefficients of finite-difference stencil
+#ifdef WARPX_DIM_RZ
+    m_dr = cell_size[0];
+    m_nmodes = WarpX::GetInstance().n_rz_azimuthal_modes;
+    m_rmin = WarpX::GetInstance().Geom(0).ProbLo(0);
+    if (fdtd_algo == MaxwellSolverAlgo::Yee){
+        CylindricalYeeAlgorithm::InitializeStencilCoefficients( cell_size,
+            stencil_coefs_r, stencil_coefs_z );
+#else
+    if (fdtd_algo == MaxwellSolverAlgo::Yee){
+        YeeAlgorithm::InitializeStencilCoefficients( cell_size,
+            stencil_coefs_x, stencil_coefs_y, stencil_coefs_z );
+    } else if (fdtd_algo == MaxwellSolverAlgo::CKC) {
+        CKCAlgorithm::InitializeStencilCoefficients( cell_size,
+            stencil_coefs_x, stencil_coefs_y, stencil_coefs_z );
+#endif
+    } else {
+        amrex::Abort("Unknown algorithm");
+    }
+};
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/Make.package b/Source/FieldSolver/FiniteDifferenceSolver/Make.package
index 27ac23846..289ed98f2 100644
--- a/Source/FieldSolver/FiniteDifferenceSolver/Make.package
+++ b/Source/FieldSolver/FiniteDifferenceSolver/Make.package
@@ -1,4 +1,5 @@
 CEXE_headers += FiniteDifferenceSolver.H
+CEXE_sources += FiniteDifferenceSolver.cpp
 CEXE_sources += EvolveB.cpp
 
 INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/FieldSolver/FiniteDifferenceSolver