use 4th order finite difference for Ampere's law in RZ

roelof-groenewald · Sep 6, 2023 · 111ec70 · 111ec70
1 parent 5fc56d4
commit 111ec70
Show file tree

Hide file tree

Showing 2 changed files with 66 additions and 4 deletions.
diff --git a/...e/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CylindricalYeeAlgorithm.H b/...e/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CylindricalYeeAlgorithm.H
@@ -104,6 +104,23 @@ struct CylindricalYeeAlgorithm {
         return 1._rt/r * inv_dr*( (r+0.5_rt*dr)*F(i,j,k,comp) - (r-0.5_rt*dr)*F(i-1,j,k,comp) );
     }
 
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    static amrex::Real DownwardDrr_over_r_4thOrder (
+        amrex::Array4<amrex::Real> const& F,
+        amrex::Real const r, amrex::Real const dr,
+        amrex::Real const * const coefs_r, int const n_coefs_r,
+        int const i, int const j, int const k, int const comp ) {
+
+        using namespace amrex;
+        ignore_unused(n_coefs_r);
+
+        Real const inv_dr = coefs_r[0]/12.0;
+        return 1._rt/r * inv_dr*(
+            -(r+1.5_rt*dr)*F(i+1,j,k,comp) + 8.0*(r+0.5_rt*dr)*F(i,j,k,comp)
+            - 8.0*(r-0.5_rt*dr)*F(i-1,j,k,comp) + (r-1.5_rt*dr)*F(i-2,j,k,comp)
+        );
+    }
+
     /**
      * Perform derivative along r on a cell-centered grid, from a nodal field `F` */
     AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
@@ -134,6 +151,24 @@ struct CylindricalYeeAlgorithm {
         return inv_dr*( F(i,j,k,comp) - F(i-1,j,k,comp) );
     }
 
+    /**
+     * Perform derivative along r on a nodal grid, from a cell-centered field `F` */
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    static amrex::Real DownwardDr4thOrder (
+        amrex::Array4<amrex::Real> const& F,
+        amrex::Real const * const coefs_r, int const n_coefs_r,
+        int const i, int const j, int const k, int const comp ) {
+
+        using namespace amrex;
+        ignore_unused(n_coefs_r);
+
+        Real const inv_dr = coefs_r[0]/12.0;
+        return inv_dr*(
+            -F(i+1,j,k,comp) + 8.0*F(i,j,k,comp) - 8.0*F(i-1,j,k,comp)
+            +F(i-2,j,k,comp)
+        );
+    }
+
     /**
      * Perform derivative along z on a cell-centered grid, from a nodal field `F` */
     AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
@@ -162,6 +197,23 @@ struct CylindricalYeeAlgorithm {
         return inv_dz*( F(i,j,k,comp) - F(i,j-1,k,comp) );
     }
 
+    /**
+     * Perform derivative along z on a nodal grid, from a cell-centered field `F` */
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    static amrex::Real DownwardDz4thOrder (
+        amrex::Array4<amrex::Real> const& F,
+        amrex::Real const * const coefs_z, int const n_coefs_z,
+        int const i, int const j, int const k, int const comp ) {
+
+        amrex::ignore_unused(n_coefs_z);
+
+        amrex::Real const inv_dz = coefs_z[0]/12.0;
+        return inv_dz*(
+            -F(i,j+1,k,comp) + 8.0*F(i,j,k,comp) - 8.0*F(i,j-1,k,comp)
+            +F(i,j-2,k,comp)
+        );
+    }
+
 };
 
 #endif // WARPX_FINITE_DIFFERENCE_ALGORITHM_CYLINDRICAL_YEE_H_
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/HybridPICSolveE.cpp b/Source/FieldSolver/FiniteDifferenceSolver/HybridPICSolveE.cpp
@@ -130,8 +130,11 @@ void FiniteDifferenceSolver::CalculateCurrentAmpereCylindrical (
                 if (lr(i, j, 0) <= 0) return;
 #endif
                 // Mode m=0
+                // Jr(i, j, 0, 0) = one_over_mu0 * (
+                //     - T_Algo::DownwardDz(Bt, coefs_z, n_coefs_z, i, j, 0, 0)
+                // );
                 Jr(i, j, 0, 0) = one_over_mu0 * (
-                    - T_Algo::DownwardDz(Bt, coefs_z, n_coefs_z, i, j, 0, 0)
+                    - T_Algo::DownwardDz4thOrder(Bt, coefs_z, n_coefs_z, i, j, 0, 0)
                 );
 
                 // Higher-order modes
@@ -161,9 +164,13 @@ void FiniteDifferenceSolver::CalculateCurrentAmpereCylindrical (
                 // Off-axis, regular curl
                 if (r > 0.5_rt*dr) {
                     // Mode m=0
+                    // Jt(i, j, 0, 0) = one_over_mu0 * (
+                    //     - T_Algo::DownwardDr(Bz, coefs_r, n_coefs_r, i, j, 0, 0)
+                    //     + T_Algo::DownwardDz(Br, coefs_z, n_coefs_z, i, j, 0, 0)
+                    // );
                     Jt(i, j, 0, 0) = one_over_mu0 * (
-                        - T_Algo::DownwardDr(Bz, coefs_r, n_coefs_r, i, j, 0, 0)
-                        + T_Algo::DownwardDz(Br, coefs_z, n_coefs_z, i, j, 0, 0)
+                        - T_Algo::DownwardDr4thOrder(Bz, coefs_r, n_coefs_r, i, j, 0, 0)
+                        + T_Algo::DownwardDz4thOrder(Br, coefs_z, n_coefs_z, i, j, 0, 0)
                     );
 
                     // Higher-order modes
@@ -208,8 +215,11 @@ void FiniteDifferenceSolver::CalculateCurrentAmpereCylindrical (
                 // Off-axis, regular curl
                 if (r > 0.5_rt*dr) {
                     // Mode m=0
+                    // Jz(i, j, 0, 0) = one_over_mu0 * (
+                    //    T_Algo::DownwardDrr_over_r(Bt, r, dr, coefs_r, n_coefs_r, i, j, 0, 0)
+                    // );
                     Jz(i, j, 0, 0) = one_over_mu0 * (
-                       T_Algo::DownwardDrr_over_r(Bt, r, dr, coefs_r, n_coefs_r, i, j, 0, 0)
+                       T_Algo::DownwardDrr_over_r_4thOrder(Bt, r, dr, coefs_r, n_coefs_r, i, j, 0, 0)
                     );
                     // Higher-order modes
                     for (int m=1 ; m<nmodes ; m++) {