1 files changed, 78 insertions, 0 deletions
diff --git a/Source/Filter/NCIGodfreyFilter.cpp b/Source/Filter/NCIGodfreyFilter.cpp
new file mode 100644
index 000000000..3f589260a
--- /dev/null
+++ b/Source/Filter/NCIGodfreyFilter.cpp
@@ -0,0 +1,78 @@
+#include <WarpX.H>
+#include <NCIGodfreyFilter.H>
+#include <NCIGodfreyTables.H>
+
+#ifdef _OPENMP
+#include <omp.h>
+#endif
+
+using namespace amrex;
+
+NCIGodfreyFilter::NCIGodfreyFilter(godfrey_coeff_set coeff_set_, amrex::Real cdtodz_, amrex::Real l_lower_order_in_v_){
+    // Store parameters into class data members
+    coeff_set = coeff_set_;
+    cdtodz = cdtodz_;
+    l_lower_order_in_v = l_lower_order_in_v_;
+    // NCI Godfrey filter has fixed size, and is applied along z only.
+#if (AMREX_SPACEDIM == 3)
+    stencil_length_each_dir = {1,1,5};
+    slen = {1,1,5};
+#else
+    stencil_length_each_dir = {1,5};
+    slen = {1,5,1};
+#endif    
+}
+
+void NCIGodfreyFilter::ComputeStencils(){
+    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
+#if ( AMREX_SPACEDIM == 3 )
+        slen.z==5,
+#else
+        slen.y==5, 
+#endif
+        "ERROR: NCI filter requires 5 points in z");
+    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(l_lower_order_in_v==1, 
+        "ERROR: NCI corrector requires l_lower_order_in_v=1, i.e., Galerkin scheme");
+    
+    // Interpolate coefficients from the table, and store into prestencil.
+    int index = tab_length*cdtodz;
+    index = min(index, tab_length-2);
+    index = max(index, 0);
+    Real weight_right = cdtodz - index/tab_length;
+    Real prestencil[4];
+    for(int i=0; i<tab_width; i++){
+        if        (coeff_set == godfrey_coeff_set::Ex_Ey_Bz) {
+            prestencil[i] = (1-weight_right)*table_nci_godfrey_Ex_Ey_Bz[index  ][i] + 
+                                weight_right*table_nci_godfrey_Ex_Ey_Bz[index+1][i];
+        } else if (coeff_set == godfrey_coeff_set::Bx_By_Ez) {
+            prestencil[i] = (1-weight_right)*table_nci_godfrey_Bx_By_Ez[index  ][i] + 
+                                weight_right*table_nci_godfrey_Bx_By_Ez[index+1][i];
+        }
+    }
+
+    // Compute stencil_z
+    stencil_z.resize( 5 );
+    stencil_z[0] =  (256 + 128*prestencil[0] + 96*prestencil[1] + 80*prestencil[2] + 70*prestencil[3]) / 256;
+    stencil_z[1] = -(       64*prestencil[0] + 64*prestencil[1] + 60*prestencil[2] + 56*prestencil[3]) / 256;
+    stencil_z[2] =  (                          16*prestencil[1] + 24*prestencil[2] + 28*prestencil[3]) / 256;
+    stencil_z[3] = -(                                              4*prestencil[2] +  8*prestencil[3]) / 256;
+    stencil_z[4] =  (                                                                 1*prestencil[3]) / 256;
+
+    // Compute stencil_x and stencil_y (no filter in these directions, 
+    // so only 1 coeff, equal to 1)
+    stencil_x.resize(1);
+    stencil_x[0] = 1.;
+#if (AMREX_SPACEDIM == 3)
+    stencil_y.resize(1);
+    stencil_y[0] = 1.;
+#endif
+
+    // Due to the way Filter::DoFilter() is written, 
+    // coefficient 0 has to be /2
+    stencil_x[0] /= 2.;
+#if (AMREX_SPACEDIM == 3)
+    stencil_y[0] /= 2.;
+#endif
+    stencil_z[0] /= 2.;
+    
+}