Merge pull request #124 from ECP-WarpX/clean_godfrey

Move NCI Godfrey filter to C++

Author: MaxThevenet

Examples/Physics_applications/laser_acceleration/inputs.2d.boost

@@ -1,8 +1,8 @@
 #################################
 ######### BOX PARAMETERS ########
 #################################
-# max_step = 2700
-stop_time = 1.9e-12
+max_step = 1000
+# stop_time = 1.9e-12
 amr.n_cell =  128 1024
 amr.max_grid_size = 64
 amr.blocking_factor = 32
@@ -39,7 +39,7 @@ warpx.serialize_ics = 1
 #################################
 ####### BOOST PARAMETERS ########
 #################################
-warpx.gamma_boost = 10.
+warpx.gamma_boost = 30.
 warpx.boost_direction = z
 warpx.do_boosted_frame_diagnostic = 1
 warpx.num_snapshots_lab = 7
@@ -62,11 +62,11 @@ electrons.momentum_distribution_type = "gaussian"
 electrons.xmin = -120.e-6
 electrons.xmax =  120.e-6
 electrons.zmin = 0.5e-3
-electrons.zmax = .0035
+electrons.zmax = 1.
 electrons.profile                 = "predefined"
 electrons.predefined_profile_name = "parabolic_channel"
 #         predefined_profile_params = z_start   ramp_up   plateau   ramp_down   rc       n0
-electrons.predefined_profile_params = .5e-3     .5e-3     2.e-3     .5e-3       50.e-6   3.5e24
+electrons.predefined_profile_params = .5e-3     .5e-3     2.e-3     .5e-3       50.e-6   3.5e25
 electrons.do_continuous_injection = 1
 
 ions.charge = q_e
@@ -77,11 +77,11 @@ ions.momentum_distribution_type = "gaussian"
 ions.xmin = -120.e-6
 ions.xmax =  120.e-6
 ions.zmin = 0.5e-3
-ions.zmax = .0035
+ions.zmax = 1.
 ions.profile                 = "predefined"
 ions.predefined_profile_name = "parabolic_channel"
 #    predefined_profile_params = z_start   ramp_up   plateau   ramp_down   rc       n0
-ions.predefined_profile_params = .5e-3     .5e-3     2.e-3     .5e-3       50.e-6   3.5e24
+ions.predefined_profile_params = .5e-3     .5e-3     2.e-3     .5e-3       50.e-6   3.5e25
 ions.do_continuous_injection = 1
 
 beam.charge = -q_e

Examples/Tests/SingleParticle/bilinear_filter_analysis.py

@@ -0,0 +1,38 @@
+#! /usr/bin/env python
+
+import sys
+import yt ; yt.funcs.mylog.setLevel(0)
+import numpy as np
+from scipy import signal
+
+# Build Jx without filter (from other simulation)
+my_F_nofilter = np.zeros([16,16])
+my_F_nofilter[8,8] = -1.601068065642412e-11
+my_F_nofilter[8,7] = -1.601068065642412e-11
+
+# Build 2D filter
+filter0 = np.array([.25,.5,.25])
+my_order = [1,5]
+my_filterx = filter0
+my_filtery = filter0
+while my_order[0]>1:
+    my_filterx = np.convolve(my_filterx,filter0)
+    my_order[0] -= 1
+while my_order[1]>1:
+    my_filtery = np.convolve(my_filtery,filter0)
+    my_order[1] -= 1
+my_filter = my_filterx[:,None]*my_filtery
+
+# Apply filter. my_F_filetered is the theoretical value for filtered field
+my_F_filtered = signal.convolve2d(my_F_nofilter, my_filter, boundary='symm', mode='same')
+
+# Get simulation result for F_filtered
+filename = sys.argv[1]
+ds = yt.load( filename )
+sl = yt.SlicePlot(ds, 2, 'jx', aspect=1)
+all_data_level_0 = ds.covering_grid(level=0,left_edge=ds.domain_left_edge, dims=ds.domain_dimensions)
+F_filtered = all_data_level_0['boxlib', 'jx'].v.squeeze()
+
+# Compare theory and PIC for filtered value
+error = np.sum( np.abs(F_filtered - my_F_filtered) ) / np.sum( np.abs(my_F_filtered) )
+assert( error < 1.e-14 )

Examples/Tests/SingleParticle/inputs

@@ -0,0 +1,25 @@
+max_step = 1
+amr.n_cell = 16 16
+amr.max_level = 0
+amr.blocking_factor = 8
+amr.max_grid_size = 8
+amr.plot_int = 1
+geometry.coord_sys   = 0
+geometry.is_periodic = 0 0
+geometry.prob_lo     = -8 -12
+geometry.prob_hi     =  8  12
+warpx.do_pml = 0
+algo.current_deposition = 1
+algo.charge_deposition = 1
+algo.field_gathering = 1
+algo.particle_pusher = 0
+warpx.cfl = 1.0
+
+particles.nspecies = 1
+particles.species_names = electron
+electron.charge = -q_e
+electron.mass = m_e
+electron.injection_style = "SingleParticle"
+electron.single_particle_pos = 0.0   0.0 0.0
+electron.single_particle_vel = 1.e20 0.0 0.0  # gamma*beta
+electron.single_particle_weight = 1.0

Regression/WarpX-tests.ini

@@ -102,6 +102,21 @@ compileTest = 0
 doVis = 0
 analysisRoutine = Examples/Modules/nci_corrector/ncicorr_analysis.py
 
+[bilinear_filter]
+buildDir = .
+inputFile = Examples/Tests/SingleParticle/inputs
+runtime_params = warpx.use_filter=1 warpx.filter_npass_each_dir=1 5
+dim = 2
+addToCompileString =
+restartTest = 0
+useMPI = 1
+numprocs = 2
+useOMP = 1
+numthreads = 2
+compileTest = 0
+doVis = 0
+analysisRoutine = Examples/Tests/SingleParticle/bilinear_filter_analysis.py
+
 [Langmuir_2d]
 buildDir = .
 inputFile = Examples/Tests/Langmuir/inputs.rt

Source/Filter/BilinearFilter.H

@@ -1,30 +1,17 @@
-#include <AMReX_MultiFab.H>
-#include <AMReX_CudaContainers.H>
+#include <Filter.H>
 
-#ifndef WARPX_FILTER_H_
-#define WARPX_FILTER_H_
+#ifndef WARPX_BILIN_FILTER_H_
+#define WARPX_BILIN_FILTER_H_
 
-class BilinearFilter
+class BilinearFilter : public Filter
 {
 public:
-    BilinearFilter () = default;
+
+    BilinearFilter() = default;
 
     void ComputeStencils();
-    void ApplyStencil(amrex::MultiFab& dstmf, const amrex::MultiFab& srcmf, int scomp=0, int dcomp=0, int ncomp=10000);
 
     amrex::IntVect npass_each_dir;
-    amrex::IntVect stencil_length_each_dir;
-
-    // public for cuda
-    void Filter(const amrex::Box& tbx,
-                amrex::Array4<amrex::Real const> const& tmp,
-                amrex::Array4<amrex::Real      > const& dst,
-                int scomp, int dcomp, int ncomp);
-
-private:
-
-    amrex::Gpu::ManagedVector<amrex::Real> stencil_x, stencil_y, stencil_z;
 
-    amrex::Dim3 slen;
 };
-#endif
+#endif // #ifndef WARPX_BILIN_FILTER_H_

Source/Filter/BilinearFilter.cpp

@@ -68,173 +68,3 @@ void BilinearFilter::ComputeStencils(){
     slen.z = 1;
 #endif
 }
-
-
-#ifdef AMREX_USE_CUDA
-
-void
-BilinearFilter::ApplyStencil (MultiFab& dstmf, const MultiFab& srcmf, int scomp, int dcomp, int ncomp)
-{
-    BL_PROFILE("BilinearFilter::ApplyStencil()");
-    ncomp = std::min(ncomp, srcmf.nComp());
-
-    for (MFIter mfi(dstmf); mfi.isValid(); ++mfi)
-    {
-        const auto& src = srcmf.array(mfi);
-        const auto& dst = dstmf.array(mfi);
-        const Box& tbx = mfi.growntilebox();
-        const Box& gbx = amrex::grow(tbx,stencil_length_each_dir-1);
-
-        // tmpfab has enough ghost cells for the stencil
-        FArrayBox tmp_fab(gbx,ncomp);
-        Elixir tmp_eli = tmp_fab.elixir();  // Prevent the tmp data from being deleted too early
-        auto const& tmp = tmp_fab.array();
-
-        // Copy values in srcfab into tmpfab
-        const Box& ibx = gbx & srcmf[mfi].box();
-        AMREX_PARALLEL_FOR_4D ( gbx, ncomp, i, j, k, n,
-        {
-            if (ibx.contains(IntVect(AMREX_D_DECL(i,j,k)))) {
-                tmp(i,j,k,n) = src(i,j,k,n+scomp);
-            } else {
-                tmp(i,j,k,n) = 0.0;
-            }
-        });
-
-        // Apply filter
-        Filter(tbx, tmp, dst, 0, dcomp, ncomp);
-    }
-}
-
-void BilinearFilter::Filter (const Box& tbx,
-                             Array4<Real const> const& tmp,
-                             Array4<Real      > const& dst,
-                             int scomp, int dcomp, int ncomp)
-{
-    amrex::Real const* AMREX_RESTRICT sx = stencil_x.data();
-    amrex::Real const* AMREX_RESTRICT sy = stencil_y.data();
-    amrex::Real const* AMREX_RESTRICT sz = stencil_z.data();
-    Dim3 slen_local = slen;
-    AMREX_PARALLEL_FOR_4D ( tbx, ncomp, i, j, k, n,
-    {
-        Real d = 0.0;
-
-        for         (int iz=0; iz < slen_local.z; ++iz){
-            for     (int iy=0; iy < slen_local.y; ++iy){
-                for (int ix=0; ix < slen_local.x; ++ix){
-#if (AMREX_SPACEDIM == 3)        
-                    Real sss = sx[ix]*sy[iy]*sz[iz];
-#else
-                    Real sss = sx[ix]*sz[iy];
-#endif                        
-#if (AMREX_SPACEDIM == 3)
-                    d += sss*( tmp(i-ix,j-iy,k-iz,scomp+n)
-                              +tmp(i+ix,j-iy,k-iz,scomp+n)
-                              +tmp(i-ix,j+iy,k-iz,scomp+n)
-                              +tmp(i+ix,j+iy,k-iz,scomp+n)
-                              +tmp(i-ix,j-iy,k+iz,scomp+n)
-                              +tmp(i+ix,j-iy,k+iz,scomp+n)
-                              +tmp(i-ix,j+iy,k+iz,scomp+n)
-                              +tmp(i+ix,j+iy,k+iz,scomp+n));
-#else
-                    d += sss*( tmp(i-ix,j-iy,k,scomp+n)
-                              +tmp(i+ix,j-iy,k,scomp+n)
-                              +tmp(i-ix,j+iy,k,scomp+n)
-                              +tmp(i+ix,j+iy,k,scomp+n));
-#endif
-                }
-            }
-        }
-
-        dst(i,j,k,dcomp+n) = d;
-    });
-}
-
-#else
-
-void
-BilinearFilter::ApplyStencil (MultiFab& dstmf, const MultiFab& srcmf, int scomp, int dcomp, int ncomp)
-{
-    BL_PROFILE("BilinearFilter::ApplyStencil()");
-    ncomp = std::min(ncomp, srcmf.nComp());
-#ifdef _OPENMP
-#pragma omp parallel
-#endif
-    {
-        FArrayBox tmpfab;
-        for (MFIter mfi(dstmf,true); mfi.isValid(); ++mfi){
-            const auto& srcfab = srcmf[mfi];
-            auto& dstfab = dstmf[mfi];
-            const Box& tbx = mfi.growntilebox();
-            const Box& gbx = amrex::grow(tbx,stencil_length_each_dir-1);
-            // tmpfab has enough ghost cells for the stencil
-            tmpfab.resize(gbx,ncomp);
-            tmpfab.setVal(0.0, gbx, 0, ncomp);
-            // Copy values in srcfab into tmpfab
-            const Box& ibx = gbx & srcfab.box();
-            tmpfab.copy(srcfab, ibx, scomp, ibx, 0, ncomp);
-            // Apply filter
-            Filter(tbx, tmpfab.array(), dstfab.array(), 0, dcomp, ncomp);
-        }
-    }
-}
-
-void BilinearFilter::Filter (const Box& tbx,
-                             Array4<Real const> const& tmp,
-                             Array4<Real      > const& dst,
-                             int scomp, int dcomp, int ncomp)
-{
-    const auto lo = amrex::lbound(tbx);
-    const auto hi = amrex::ubound(tbx);
-    // tmp and dst are of type Array4 (Fortran ordering)
-    amrex::Real const* AMREX_RESTRICT sx = stencil_x.data();
-    amrex::Real const* AMREX_RESTRICT sy = stencil_y.data();
-    amrex::Real const* AMREX_RESTRICT sz = stencil_z.data();
-    for (int n = 0; n < ncomp; ++n) {
-        // Set dst value to 0.
-        for         (int k = lo.z; k <= hi.z; ++k) {
-            for     (int j = lo.y; j <= hi.y; ++j) {
-                for (int i = lo.x; i <= hi.x; ++i) {
-                    dst(i,j,k,dcomp+n) = 0.0;
-                }
-            }
-        }
-        // 3 nested loop on 3D stencil
-        for         (int iz=0; iz < slen.z; ++iz){
-            for     (int iy=0; iy < slen.y; ++iy){
-                for (int ix=0; ix < slen.x; ++ix){
-#if (AMREX_SPACEDIM == 3)        
-                    Real sss = sx[ix]*sy[iy]*sz[iz];
-#else
-                    Real sss = sx[ix]*sz[iy];
-#endif
-                    // 3 nested loop on 3D array
-                    for         (int k = lo.z; k <= hi.z; ++k) {
-                        for     (int j = lo.y; j <= hi.y; ++j) {
-                            AMREX_PRAGMA_SIMD
-                            for (int i = lo.x; i <= hi.x; ++i) {
-#if (AMREX_SPACEDIM == 3)
-                                dst(i,j,k,dcomp+n) += sss*(tmp(i-ix,j-iy,k-iz,scomp+n)
-                                                          +tmp(i+ix,j-iy,k-iz,scomp+n)
-                                                          +tmp(i-ix,j+iy,k-iz,scomp+n)
-                                                          +tmp(i+ix,j+iy,k-iz,scomp+n)
-                                                          +tmp(i-ix,j-iy,k+iz,scomp+n)
-                                                          +tmp(i+ix,j-iy,k+iz,scomp+n)
-                                                          +tmp(i-ix,j+iy,k+iz,scomp+n)
-                                                          +tmp(i+ix,j+iy,k+iz,scomp+n));
-#else
-                                dst(i,j,k,dcomp+n) += sss*(tmp(i-ix,j-iy,k,scomp+n)
-                                                          +tmp(i+ix,j-iy,k,scomp+n)
-                                                          +tmp(i-ix,j+iy,k,scomp+n)
-                                                          +tmp(i+ix,j+iy,k,scomp+n));
-#endif
-                            }
-                        }
-                    }
-                }
-            }
-        }
-    }
-}
-
-#endif

Source/Filter/Filter.H

@@ -0,0 +1,43 @@
+#include <AMReX_MultiFab.H>
+#include <AMReX_CudaContainers.H>
+
+#ifndef WARPX_FILTER_H_
+#define WARPX_FILTER_H_
+
+class Filter
+{
+public:
+    Filter () = default;
+
+    // Apply stencil on MultiFab.
+    // Guard cells are handled inside this function
+    void ApplyStencil(amrex::MultiFab& dstmf, 
+                              const amrex::MultiFab& srcmf, int scomp=0, 
+                              int dcomp=0, int ncomp=10000);
+
+    // Apply stencil on a FabArray.
+    void ApplyStencil (amrex::FArrayBox& dstfab, 
+                       const amrex::FArrayBox& srcfab, const amrex::Box& tbx, 
+                       int scomp=0, int dcomp=0, int ncomp=10000);
+
+    // public for cuda
+    void DoFilter(const amrex::Box& tbx,
+                          amrex::Array4<amrex::Real const> const& tmp,
+                          amrex::Array4<amrex::Real      > const& dst,
+                          int scomp, int dcomp, int ncomp);
+
+    // In 2D, stencil_length_each_dir = {length(stencil_x), length(stencil_z)}
+    amrex::IntVect stencil_length_each_dir;
+
+protected:
+    // Stencil along each direction.
+    // in 2D, stencil_y is not initialized.
+    amrex::Gpu::ManagedVector<amrex::Real> stencil_x, stencil_y, stencil_z;
+    // Length of each stencil.
+    // In 2D, slen = {length(stencil_x), length(stencil_z), 1}
+    amrex::Dim3 slen;
+
+private:
+
+};
+#endif // #ifndef WARPX_FILTER_H_