From 3d08a02a6d3401d5e0da2fdac5e271ab2425c51e Mon Sep 17 00:00:00 2001
From: Edoardo Zoni <59625522+EZoni@users.noreply.github.com>
Date: Mon, 21 Mar 2022 16:24:07 -0700
Subject: Hybrid Algorithms: Add FDTD Support for Centering (#2327)

* Hybrid Algorithms: Add FDTD Support for Centering

* Current Centering Incompatible with Esirkepov Deposition

* Remove FDTD v. PSATD Runtime Conditions

* Add CI Test

* Rename coefs as coeffs in getFornbergStencilCoefficients

* amrex::Vector: Use at() Instead of Operator []

* Reset Benchmark

* Reset Benchmark

* Cleaning

* Reduce Size of New CI Test

* Reset Benchmark

* Fix Bug
---
 .../FieldSolver/SpectralSolver/SpectralKSpace.cpp  | 41 ++--------------------
 1 file changed, 2 insertions(+), 39 deletions(-)

(limited to 'Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp')

diff --git a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp
index ca26682a9..b9afc63b4 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp
@@ -7,6 +7,7 @@
  */
 #include "SpectralKSpace.H"
 
+#include "WarpX.H"
 #include "Utils/TextMsg.H"
 #include "Utils/WarpXConst.H"
 
@@ -216,7 +217,7 @@ SpectralKSpace::getModifiedKComponent( const DistributionMapping& dm,
     } else {
 
         // Compute real-space stencil coefficients
-        Vector<Real> h_stencil_coef = getFornbergStencilCoefficients(n_order, nodal);
+        Vector<Real> h_stencil_coef = WarpX::getFornbergStencilCoefficients(n_order, nodal);
         Gpu::DeviceVector<Real> d_stencil_coef(h_stencil_coef.size());
         Gpu::copyAsync(Gpu::hostToDevice, h_stencil_coef.begin(), h_stencil_coef.end(),
                        d_stencil_coef.begin());
@@ -277,41 +278,3 @@ SpectralKSpace::getModifiedKComponent( const DistributionMapping& dm,
     }
     return modified_k_comp;
 }
-
-Vector<Real>
-getFornbergStencilCoefficients(const int n_order, const bool nodal)
-{
-    WARPX_ALWAYS_ASSERT_WITH_MESSAGE(n_order % 2 == 0, "n_order must be even");
-
-    const int m = n_order / 2;
-    Vector<Real> coefs;
-    coefs.resize(m);
-
-    // There are closed-form formula for these coefficients, but they result in
-    // an overflow when evaluated numerically. One way to avoid the overflow is
-    // to calculate the coefficients by recurrence.
-
-    // Coefficients for nodal (that is, centered) finite-difference approximation
-    if (nodal == true) {
-       // First coefficient
-       coefs[0] = m * 2. / (m+1);
-       // Other coefficients by recurrence
-       for (int n = 1; n < m; n++) {
-           coefs[n] = - (m-n) * 1. / (m+n+1) * coefs[n-1];
-       }
-    }
-    // Coefficients for staggered finite-difference approximation
-    else {
-       Real prod = 1.;
-       for (int k = 1; k < m+1; k++) {
-           prod *= (m + k) / (4. * k);
-       }
-       // First coefficient
-       coefs[0] = 4 * m * prod * prod;
-       // Other coefficients by recurrence
-       for (int n = 1; n < m; n++) {
-           coefs[n] = - ((2*n-1) * (m-n)) * 1. / ((2*n+1) * (m+n)) * coefs[n-1];
-       }
-    }
-    return coefs;
-}
-- 
cgit v1.2.3