From 76f3c0cb1ba717de8fa571a06fbb0267b4f83237 Mon Sep 17 00:00:00 2001
From: Remi Lehe <remi.lehe@normalesup.org>
Date: Mon, 22 Apr 2019 09:56:04 -0700
Subject: Add calculation of Fonberg coefficients

---
 .../FieldSolver/SpectralSolver/SpectralKSpace.cpp  | 76 +++++++++++++++-------
 1 file changed, 54 insertions(+), 22 deletions(-)

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

diff --git a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp
index 90c20b01d..77aa7342f 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp
@@ -66,9 +66,38 @@ SpectralKSpace::getKComponent( const DistributionMapping& dm, const int i_dim )
     return k_comp;
 }
 
+SpectralShiftFactor
+SpectralKSpace::getSpectralShiftFactor(
+        const DistributionMapping& dm, const int i_dim, const int shift_type ) const
+{
+    // Initialize an empty ManagedVector in each box
+    SpectralShiftFactor shift_factor = SpectralShiftFactor( spectralspace_ba, dm );
+    // Loop over boxes
+    for ( MFIter mfi(spectralspace_ba, dm); mfi.isValid(); ++mfi ){
+        const ManagedVector<Real>& k = k_vec[i_dim][mfi];
+        ManagedVector<Complex>& shift = shift_factor[mfi];
+
+        // Allocate shift coefficients
+        shift.resize( k.size() );
+
+        // Fill the shift coefficients
+        Real sign = 0;
+        switch (shift_type){
+            case ShiftType::CenteredToNodal: sign = -1.; break;
+            case ShiftType::NodalToCentered: sign = 1.;
+        }
+        constexpr Complex I{0,1};
+        for (int i=0; i<k.size(); i++ ){
+            shift[i] = std::exp( I*sign*k[i]*0.5*dx[i_dim] );
+        }
+    }
+    return shift_factor;
+}
+
 KVectorComponent
 SpectralKSpace::getModifiedKComponent(
-        const DistributionMapping& dm, const int i_dim, const int order ) const
+        const DistributionMapping& dm, const int i_dim,
+        const int n_order, const bool nodal ) const
 {
     // Initialize an empty ManagedVector in each box
     KVectorComponent modified_k_comp = KVectorComponent( spectralspace_ba, dm );
@@ -90,30 +119,33 @@ SpectralKSpace::getModifiedKComponent(
     return modified_k_comp;
 }
 
-SpectralShiftFactor
-SpectralKSpace::getSpectralShiftFactor(
-        const DistributionMapping& dm, const int i_dim, const int shift_type ) const
+/* TODO: Documentation: point to Fonberg paper ; explain recurrence relation
+ */
+Array<Real>
+getFonbergStencilCoefficients( const int n_order, const bool nodal )
 {
-    // Initialize an empty ManagedVector in each box
-    SpectralShiftFactor shift_factor = SpectralShiftFactor( spectralspace_ba, dm );
-    // Loop over boxes
-    for ( MFIter mfi(spectralspace_ba, dm); mfi.isValid(); ++mfi ){
-        const ManagedVector<Real>& k = k_vec[i_dim][mfi];
-        ManagedVector<Complex>& shift = shift_factor[mfi];
+    AMREX_ALWAYS_ASSERT_WITH_MESSAGE( n_order%2 == 0, "n_order should be even.");
+    const int m = n_order/2;
+    Array<Real> stencil_coef;
+    stencil_coef.resize( m+1 );
 
-        // Allocate shift coefficients
-        shift.resize( k.size() );
-
-        // Fill the shift coefficients
-        Real sign = 0;
-        switch (shift_type){
-            case ShiftType::CenteredToNodal: sign = -1.; break;
-            case ShiftType::NodalToCentered: sign = 1.;
+    // Coefficients for nodal (a.k.a. centered) finite-difference
+    if (nodal == true) {
+        stencil_coef[0] = -2.; // First coefficient
+        for (int n=1; n<m+1; n++){ // Get the other coefficients by recurrence
+            stencil_coef[n] = - (m+1-n)*1./(m+n)*stencil_coef[n-1];
         }
-        constexpr Complex I{0,1};
-        for (int i=0; i<k.size(); i++ ){
-            shift[i] = std::exp( I*sign*k[i]*0.5*dx[i_dim] );
+    }
+    // Coefficients for staggered finite-difference
+    else {
+        Real prod = 1.;
+        for (int k=1; k<m+1; k++){
+            prod *= (m+k)*1./(4*k);
+        }
+        stencil_coef[0] = 4*m*prod*prod; // First coefficient
+        for (int n=1; n<m+1; n++){ // Get the other coefficients by recurrence
+            stencil_coef[n] = - ((2*n-3)*(m+1-n))*1./((2*n-1)*(m-1+n))*stencil_coef[n-1];
         }
     }
-    return shift_factor;
+    return stencil_coef;
 }
-- 
cgit v1.2.3