From 4a84e9009ae2c065cb3e2ff9ab57c125d4e2ca40 Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Fri, 15 Mar 2019 17:15:40 -0700 Subject: implemented stencil application. --- Source/Parallelization/WarpXComm.cpp | 44 ++++++++++++++++++++++++++++++++++-- 1 file changed, 42 insertions(+), 2 deletions(-) (limited to 'Source/Parallelization/WarpXComm.cpp') diff --git a/Source/Parallelization/WarpXComm.cpp b/Source/Parallelization/WarpXComm.cpp index 09767c265..1d2993d8e 100644 --- a/Source/Parallelization/WarpXComm.cpp +++ b/Source/Parallelization/WarpXComm.cpp @@ -1,4 +1,3 @@ - #include #include @@ -378,6 +377,47 @@ WarpX::SyncCurrent () Vector,3> > j_cp(finest_level+1); Vector,3> > j_buf(finest_level+1); + if (WarpX::use_filter) { + for (int lev = 0; lev <= finest_level; ++lev) { + IntVect ng = current_fp[lev][0]->nGrowVect(); + ng += bilinear_filter.stencil_length_each_dir-1; + //ng += 1; + for (int idim = 0; idim < 3; ++idim) { + j_fp[lev][idim].reset(new MultiFab(current_fp[lev][idim]->boxArray(), + current_fp[lev][idim]->DistributionMap(), + 1, ng)); + bilinear_filter.ApplyStencil(*j_fp[lev][idim], *current_fp[lev][idim]); + std::swap(j_fp[lev][idim], current_fp[lev][idim]); + } + } + for (int lev = 1; lev <= finest_level; ++lev) { + IntVect ng = current_cp[lev][0]->nGrowVect(); + ng += 1; + for (int idim = 0; idim < 3; ++idim) { + j_cp[lev][idim].reset(new MultiFab(current_cp[lev][idim]->boxArray(), + current_cp[lev][idim]->DistributionMap(), + 1, ng)); + // applyFilter(*j_cp[lev][idim], *current_cp[lev][idim]); + bilinear_filter.ApplyStencil(*j_cp[lev][idim], *current_cp[lev][idim]); + std::swap(j_cp[lev][idim], current_cp[lev][idim]); + } + } + for (int lev = 1; lev <= finest_level; ++lev) { + if (current_buf[lev][0]) { + IntVect ng = current_buf[lev][0]->nGrowVect(); + ng += 1; + for (int idim = 0; idim < 3; ++idim) { + j_buf[lev][idim].reset(new MultiFab(current_buf[lev][idim]->boxArray(), + current_buf[lev][idim]->DistributionMap(), + 1, ng)); + //applyFilter(*j_buf[lev][idim], *current_buf[lev][idim]); + bilinear_filter.ApplyStencil(*j_buf[lev][idim], *current_buf[lev][idim]); + std::swap(*j_buf[lev][idim], *current_buf[lev][idim]); + } + } + } + } + /* if (WarpX::use_filter) { for (int lev = 0; lev <= finest_level; ++lev) { IntVect ng = current_fp[lev][0]->nGrowVect(); @@ -415,7 +455,7 @@ WarpX::SyncCurrent () } } } - + */ // Sum up fine patch for (int lev = 0; lev <= finest_level; ++lev) { -- cgit v1.2.3 From e6e7c1bab398f87a6bcef0af0e25204dce72e5a2 Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Mon, 18 Mar 2019 15:29:12 -0700 Subject: externalize filter function --- Source/Filter/BilinearFilter.H | 1 + Source/Filter/BilinearFilter.cpp | 49 ++++++++++++++++++++---------------- Source/Parallelization/WarpXComm.cpp | 1 + 3 files changed, 29 insertions(+), 22 deletions(-) (limited to 'Source/Parallelization/WarpXComm.cpp') diff --git a/Source/Filter/BilinearFilter.H b/Source/Filter/BilinearFilter.H index fbc572d64..ecdd81950 100644 --- a/Source/Filter/BilinearFilter.H +++ b/Source/Filter/BilinearFilter.H @@ -10,6 +10,7 @@ public: // void ApplyStencils(); //void ApplyStencils(amrex::MultiFab* dstmf, amrex::MultiFab* srcmf); void ApplyStencil (amrex::MultiFab& dstmf, const amrex::MultiFab& srcmf, int scomp = 0, int dcomp = 0, int ncomp = 10000); + void filter_2d(const amrex::Box& tbx, const amrex::Box& gbx, amrex::FArrayBox &tmpfab, amrex::FArrayBox &dstfab, int ncomp); amrex::IntVect npass_each_dir; amrex::IntVect stencil_length_each_dir; diff --git a/Source/Filter/BilinearFilter.cpp b/Source/Filter/BilinearFilter.cpp index fb89ac72d..e77da41ca 100644 --- a/Source/Filter/BilinearFilter.cpp +++ b/Source/Filter/BilinearFilter.cpp @@ -1,6 +1,7 @@ #include #include #include +#include #ifdef _OPENMP #include @@ -72,10 +73,8 @@ BilinearFilter::ApplyStencil (MultiFab& dstmf, const MultiFab& srcmf, int scomp, { Print()<<"stencil_x: "; for(int i=0; i const& tmparr = tmpfab.array(); - Array4 const& dstarr = dstfab.array(); - - Print()<<"here --- 5\n"; - - for(int i=loVect[0]; i const& tmparr = tmpfab.array(); + Array4 const& dstarr = dstfab.array(); + // filter_2d() + for(int i=loVect[0]; i Date: Mon, 18 Mar 2019 18:17:48 -0700 Subject: works with 0 pass --- Source/Filter/BilinearFilter.H | 4 ++-- Source/Filter/BilinearFilter.cpp | 33 +++++++++++++++++++++++++-------- Source/Parallelization/WarpXComm.cpp | 27 ++++++++++++++------------- 3 files changed, 41 insertions(+), 23 deletions(-) (limited to 'Source/Parallelization/WarpXComm.cpp') diff --git a/Source/Filter/BilinearFilter.H b/Source/Filter/BilinearFilter.H index ecdd81950..75d37c73c 100644 --- a/Source/Filter/BilinearFilter.H +++ b/Source/Filter/BilinearFilter.H @@ -9,8 +9,8 @@ public: void ComputeStencils(); // void ApplyStencils(); //void ApplyStencils(amrex::MultiFab* dstmf, amrex::MultiFab* srcmf); - void ApplyStencil (amrex::MultiFab& dstmf, const amrex::MultiFab& srcmf, int scomp = 0, int dcomp = 0, int ncomp = 10000); - void filter_2d(const amrex::Box& tbx, const amrex::Box& gbx, amrex::FArrayBox &tmpfab, amrex::FArrayBox &dstfab, int ncomp); + void ApplyStencil (amrex::MultiFab& dstmf, const amrex::MultiFab& srcmf, int scomp = 0, int dcomp = 0, int ncomp=10000); + void filter_2d(const amrex::Box& tbx, const amrex::Box& gbx, amrex::FArrayBox &tmpfab, amrex::FArrayBox &dstfab, int ncomp=10000); amrex::IntVect npass_each_dir; amrex::IntVect stencil_length_each_dir; diff --git a/Source/Filter/BilinearFilter.cpp b/Source/Filter/BilinearFilter.cpp index e77da41ca..0b4ab6e2a 100644 --- a/Source/Filter/BilinearFilter.cpp +++ b/Source/Filter/BilinearFilter.cpp @@ -84,6 +84,9 @@ BilinearFilter::ApplyStencil (MultiFab& dstmf, const MultiFab& srcmf, int scomp, const Box& ibx = gbx & srcfab.box(); tmpfab.copy(srcfab, ibx, scomp, ibx, 0, ncomp); + + //#define BL_TO_FORTRAN_N_ANYD(x,n) (x).dataPtr(n), AMREX_ARLIM_ANYD((x).loVect()), AMREX_ARLIM_ANYD((x).hiVect()) + //WRPX_FILTER(BL_TO_FORTRAN_BOX(tbx), // BL_TO_FORTRAN_ANYD(tmpfab), // BL_TO_FORTRAN_N_ANYD(dstfab,dcomp), @@ -95,22 +98,36 @@ BilinearFilter::ApplyStencil (MultiFab& dstmf, const MultiFab& srcmf, int scomp, void BilinearFilter::filter_2d(const Box& tbx, const Box& gbx, FArrayBox &tmpfab, FArrayBox &dstfab, int ncomp) { - const int* loVect = gbx.loVect(); - const int* hiVect = gbx.hiVect(); - Print()<<"loVect: "; for(int i=0; i<2;i++){Print()< const& tmparr = tmpfab.array(); Array4 const& dstarr = dstfab.array(); - // filter_2d() - for(int i=loVect[0]; inGrowVect(); - ng += bilinear_filter.stencil_length_each_dir-1; - //ng += 1; + //ng += bilinear_filter.stencil_length_each_dir-1; + ng += 1; for (int idim = 0; idim < 3; ++idim) { + Print()<<"idim "<boxArray(), current_fp[lev][idim]->DistributionMap(), 1, ng)); @@ -601,7 +602,7 @@ WarpX::SyncRho (amrex::Vector >& rhof, rho_f_g[lev].reset(new MultiFab(rhof[lev]->boxArray(), rhof[lev]->DistributionMap(), ncomp, ng)); - applyFilter(*rho_f_g[lev], *rhof[lev]); + bilinear_filter.ApplyStencil(*rho_f_g[lev], *rhof[lev]); std::swap(rho_f_g[lev], rhof[lev]); } for (int lev = 1; lev <= finest_level; ++lev) { @@ -611,7 +612,7 @@ WarpX::SyncRho (amrex::Vector >& rhof, rho_c_g[lev].reset(new MultiFab(rhoc[lev]->boxArray(), rhoc[lev]->DistributionMap(), ncomp, ng)); - applyFilter(*rho_c_g[lev], *rhoc[lev]); + bilinear_filter.ApplyStencil(*rho_c_g[lev], *rhoc[lev]); std::swap(rho_c_g[lev], rhoc[lev]); } for (int lev = 1; lev <= finest_level; ++lev) { @@ -622,7 +623,7 @@ WarpX::SyncRho (amrex::Vector >& rhof, rho_buf_g[lev].reset(new MultiFab(charge_buf[lev]->boxArray(), charge_buf[lev]->DistributionMap(), ncomp, ng)); - applyFilter(*rho_buf_g[lev], *charge_buf[lev]); + bilinear_filter.ApplyStencil(*rho_buf_g[lev], *charge_buf[lev]); std::swap(*rho_buf_g[lev], *charge_buf[lev]); } } @@ -754,7 +755,7 @@ WarpX::ApplyFilterandSumBoundaryJ (int lev, PatchType patch_type) IntVect ng = j[idim]->nGrowVect(); ng += 1; MultiFab jf(j[idim]->boxArray(), j[idim]->DistributionMap(), 1, ng); - applyFilter(jf, *j[idim]); + bilinear_filter.ApplyStencil(jf, *j[idim]); jf.SumBoundary(period); MultiFab::Copy(*j[idim], jf, 0, 0, 1, 0); } else { @@ -796,12 +797,12 @@ WarpX::AddCurrentFromFineLevelandSumBoundary (int lev) ng += 1; MultiFab jfc(current_cp[lev+1][idim]->boxArray(), current_cp[lev+1][idim]->DistributionMap(), 1, ng); - applyFilter(jfc, *current_cp[lev+1][idim]); + bilinear_filter.ApplyStencil(jfc, *current_cp[lev+1][idim]); // buffer patch of fine level MultiFab jfb(current_buf[lev+1][idim]->boxArray(), current_buf[lev+1][idim]->DistributionMap(), 1, ng); - applyFilter(jfb, *current_buf[lev+1][idim]); + bilinear_filter.ApplyStencil(jfb, *current_buf[lev+1][idim]); MultiFab::Add(jfb, jfc, 0, 0, 1, ng); mf.ParallelAdd(jfb, 0, 0, 1, ng, IntVect::TheZeroVector(), period); @@ -816,7 +817,7 @@ WarpX::AddCurrentFromFineLevelandSumBoundary (int lev) ng += 1; MultiFab jf(current_cp[lev+1][idim]->boxArray(), current_cp[lev+1][idim]->DistributionMap(), 1, ng); - applyFilter(jf, *current_cp[lev+1][idim]); + bilinear_filter.ApplyStencil(jf, *current_cp[lev+1][idim]); mf.ParallelAdd(jf, 0, 0, 1, ng, IntVect::TheZeroVector(), period); jf.SumBoundary(period); MultiFab::Copy(*current_cp[lev+1][idim], jf, 0, 0, 1, 0); @@ -865,7 +866,7 @@ WarpX::ApplyFilterandSumBoundaryRho (int lev, PatchType patch_type, int icomp, i IntVect ng = r->nGrowVect(); ng += 1; MultiFab rf(r->boxArray(), r->DistributionMap(), ncomp, ng); - applyFilter(rf, *r, icomp, 0, ncomp); + bilinear_filter.ApplyStencil(rf, *r, icomp, 0, ncomp); rf.SumBoundary(period); MultiFab::Copy(*r, rf, 0, icomp, ncomp, 0); } else { @@ -902,12 +903,12 @@ WarpX::AddRhoFromFineLevelandSumBoundary(int lev, int icomp, int ncomp) ng += 1; MultiFab rhofc(rho_cp[lev+1]->boxArray(), rho_cp[lev+1]->DistributionMap(), ncomp, ng); - applyFilter(rhofc, *rho_cp[lev+1], icomp, 0, ncomp); + bilinear_filter.ApplyStencil(rhofc, *rho_cp[lev+1], icomp, 0, ncomp); // buffer patch of fine level MultiFab rhofb(charge_buf[lev+1]->boxArray(), charge_buf[lev+1]->DistributionMap(), ncomp, ng); - applyFilter(rhofb, *charge_buf[lev+1], icomp, 0, ncomp); + bilinear_filter.ApplyStencil(rhofb, *charge_buf[lev+1], icomp, 0, ncomp); MultiFab::Add(rhofb, rhofc, 0, 0, ncomp, ng); mf.ParallelAdd(rhofb, 0, 0, ncomp, ng, IntVect::TheZeroVector(), period); @@ -920,7 +921,7 @@ WarpX::AddRhoFromFineLevelandSumBoundary(int lev, int icomp, int ncomp) IntVect ng = rho_cp[lev+1]->nGrowVect(); ng += 1; MultiFab rf(rho_cp[lev+1]->boxArray(), rho_cp[lev+1]->DistributionMap(), ncomp, ng); - applyFilter(rf, *rho_cp[lev+1], icomp, 0, ncomp); + bilinear_filter.ApplyStencil(rf, *rho_cp[lev+1], icomp, 0, ncomp); mf.ParallelAdd(rf, 0, 0, ncomp, ng, IntVect::TheZeroVector(), period); rf.SumBoundary(0, ncomp, period); MultiFab::Copy(*rho_cp[lev+1], rf, 0, icomp, ncomp, 0); -- cgit v1.2.3 From 2e3c6010ea627fc7995308ea0fc56cac212d8dbb Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Tue, 19 Mar 2019 10:08:39 -0700 Subject: cleaning. still works for no pass, but not >=1 pass --- Source/Filter/BilinearFilter.H | 3 +- Source/Filter/BilinearFilter.cpp | 97 +++--------------------------------- Source/Parallelization/WarpXComm.cpp | 1 - Source/WarpX.cpp | 2 +- 4 files changed, 9 insertions(+), 94 deletions(-) (limited to 'Source/Parallelization/WarpXComm.cpp') diff --git a/Source/Filter/BilinearFilter.H b/Source/Filter/BilinearFilter.H index 75d37c73c..e3823bf86 100644 --- a/Source/Filter/BilinearFilter.H +++ b/Source/Filter/BilinearFilter.H @@ -10,7 +10,8 @@ public: // void ApplyStencils(); //void ApplyStencils(amrex::MultiFab* dstmf, amrex::MultiFab* srcmf); void ApplyStencil (amrex::MultiFab& dstmf, const amrex::MultiFab& srcmf, int scomp = 0, int dcomp = 0, int ncomp=10000); - void filter_2d(const amrex::Box& tbx, const amrex::Box& gbx, amrex::FArrayBox &tmpfab, amrex::FArrayBox &dstfab, int ncomp=10000); + void filter_2d(const amrex::Box& tbx, const amrex::Box& gbx, amrex::FArrayBox &tmpfab, amrex::FArrayBox &dstfab, int dcomp=0, int ncomp=10000); + // void filter_2d(const amrex::Box& tbx, const amrex::Box& gbx, amrex::FArrayBox &tmpfab, amrex::FArrayBox &dstfab, int ncomp=10000); amrex::IntVect npass_each_dir; amrex::IntVect stencil_length_each_dir; diff --git a/Source/Filter/BilinearFilter.cpp b/Source/Filter/BilinearFilter.cpp index 0b4ab6e2a..2a423f8fe 100644 --- a/Source/Filter/BilinearFilter.cpp +++ b/Source/Filter/BilinearFilter.cpp @@ -83,110 +83,25 @@ BilinearFilter::ApplyStencil (MultiFab& dstmf, const MultiFab& srcmf, int scomp, tmpfab.setVal(0.0, gbx, 0, ncomp); const Box& ibx = gbx & srcfab.box(); tmpfab.copy(srcfab, ibx, scomp, ibx, 0, ncomp); - - - //#define BL_TO_FORTRAN_N_ANYD(x,n) (x).dataPtr(n), AMREX_ARLIM_ANYD((x).loVect()), AMREX_ARLIM_ANYD((x).hiVect()) - - //WRPX_FILTER(BL_TO_FORTRAN_BOX(tbx), - // BL_TO_FORTRAN_ANYD(tmpfab), - // BL_TO_FORTRAN_N_ANYD(dstfab,dcomp), - // ncomp); - filter_2d(tbx, gbx, tmpfab, dstfab, ncomp); + filter_2d(tbx, gbx, tmpfab, dstfab, dcomp, ncomp); } } } -void BilinearFilter::filter_2d(const Box& tbx, const Box& gbx, FArrayBox &tmpfab, FArrayBox &dstfab, int ncomp) +void BilinearFilter::filter_2d(const Box& tbx, const Box& gbx, FArrayBox &tmpfab, FArrayBox &dstfab, int dcomp, int ncomp) { - // ncomp = std::min(ncomp, srcmf.nComp()); - Print()<<"ncomp "< const& tmparr = tmpfab.array(); Array4 const& dstarr = dstfab.array(); - // for(int c=0; carray(mfi); - auto const& srcmf_fab = srcmf->array(mfi); - - amrex::ParallelFor(tbx, - [=] AMREX_GPU_DEVICE (int j, int k, int l){ - warpx_apply_filter_2d(j,k,l,dstmf_fab, srcmf_fab, - stencil_x,stencil_z, - stencil_length_each_dir[0], - stencil_length_each_dir[1]); - }); - } - } -} - */ diff --git a/Source/Parallelization/WarpXComm.cpp b/Source/Parallelization/WarpXComm.cpp index 678ecc940..ff185eca8 100644 --- a/Source/Parallelization/WarpXComm.cpp +++ b/Source/Parallelization/WarpXComm.cpp @@ -380,7 +380,6 @@ WarpX::SyncCurrent () if (WarpX::use_filter) { for (int lev = 0; lev <= finest_level; ++lev) { IntVect ng = current_fp[lev][0]->nGrowVect(); - //ng += bilinear_filter.stencil_length_each_dir-1; ng += 1; for (int idim = 0; idim < 3; ++idim) { Print()<<"idim "< Date: Tue, 19 Mar 2019 10:13:22 -0700 Subject: further cleaning --- Source/Parallelization/WarpXComm.cpp | 40 +----------------------------------- 1 file changed, 1 insertion(+), 39 deletions(-) (limited to 'Source/Parallelization/WarpXComm.cpp') diff --git a/Source/Parallelization/WarpXComm.cpp b/Source/Parallelization/WarpXComm.cpp index ff185eca8..f80258fa3 100644 --- a/Source/Parallelization/WarpXComm.cpp +++ b/Source/Parallelization/WarpXComm.cpp @@ -417,45 +417,7 @@ WarpX::SyncCurrent () } } } - /* - if (WarpX::use_filter) { - for (int lev = 0; lev <= finest_level; ++lev) { - IntVect ng = current_fp[lev][0]->nGrowVect(); - ng += 1; - for (int idim = 0; idim < 3; ++idim) { - j_fp[lev][idim].reset(new MultiFab(current_fp[lev][idim]->boxArray(), - current_fp[lev][idim]->DistributionMap(), - 1, ng)); - applyFilter(*j_fp[lev][idim], *current_fp[lev][idim]); - std::swap(j_fp[lev][idim], current_fp[lev][idim]); - } - } - for (int lev = 1; lev <= finest_level; ++lev) { - IntVect ng = current_cp[lev][0]->nGrowVect(); - ng += 1; - for (int idim = 0; idim < 3; ++idim) { - j_cp[lev][idim].reset(new MultiFab(current_cp[lev][idim]->boxArray(), - current_cp[lev][idim]->DistributionMap(), - 1, ng)); - applyFilter(*j_cp[lev][idim], *current_cp[lev][idim]); - std::swap(j_cp[lev][idim], current_cp[lev][idim]); - } - } - for (int lev = 1; lev <= finest_level; ++lev) { - if (current_buf[lev][0]) { - IntVect ng = current_buf[lev][0]->nGrowVect(); - ng += 1; - for (int idim = 0; idim < 3; ++idim) { - j_buf[lev][idim].reset(new MultiFab(current_buf[lev][idim]->boxArray(), - current_buf[lev][idim]->DistributionMap(), - 1, ng)); - applyFilter(*j_buf[lev][idim], *current_buf[lev][idim]); - std::swap(*j_buf[lev][idim], *current_buf[lev][idim]); - } - } - } - } - */ + // Sum up fine patch for (int lev = 0; lev <= finest_level; ++lev) { -- cgit v1.2.3 From 1abb26092bbaece522dbe5266118768985b97706 Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Tue, 19 Mar 2019 13:09:02 -0700 Subject: addd npass capability for 2d. ATTENTION number of ghosts cells unchanged --- Source/Filter/BilinearFilter.H | 7 ++----- Source/Filter/BilinearFilter.cpp | 9 ++++----- Source/Parallelization/WarpXComm.cpp | 36 ++++++++++++++++++++++++------------ Source/WarpX.H | 3 --- Source/WarpX.cpp | 27 --------------------------- 5 files changed, 30 insertions(+), 52 deletions(-) (limited to 'Source/Parallelization/WarpXComm.cpp') diff --git a/Source/Filter/BilinearFilter.H b/Source/Filter/BilinearFilter.H index e3823bf86..ada0f9e74 100644 --- a/Source/Filter/BilinearFilter.H +++ b/Source/Filter/BilinearFilter.H @@ -7,11 +7,8 @@ public: BilinearFilter () = default; void ComputeStencils(); - // void ApplyStencils(); - //void ApplyStencils(amrex::MultiFab* dstmf, amrex::MultiFab* srcmf); - void ApplyStencil (amrex::MultiFab& dstmf, const amrex::MultiFab& srcmf, int scomp = 0, int dcomp = 0, int ncomp=10000); - void filter_2d(const amrex::Box& tbx, const amrex::Box& gbx, amrex::FArrayBox &tmpfab, amrex::FArrayBox &dstfab, int dcomp=0, int ncomp=10000); - // void filter_2d(const amrex::Box& tbx, const amrex::Box& gbx, amrex::FArrayBox &tmpfab, amrex::FArrayBox &dstfab, int ncomp=10000); + void ApplyStencil(amrex::MultiFab& dstmf, const amrex::MultiFab& srcmf, int scomp=0, int dcomp=0, int ncomp=10000); + void Filter2d(const amrex::Box& tbx, const amrex::Box& gbx, amrex::FArrayBox &tmpfab, amrex::FArrayBox &dstfab, int dcomp=0, int ncomp=10000); amrex::IntVect npass_each_dir; amrex::IntVect stencil_length_each_dir; diff --git a/Source/Filter/BilinearFilter.cpp b/Source/Filter/BilinearFilter.cpp index c90f56161..d4016a7dc 100644 --- a/Source/Filter/BilinearFilter.cpp +++ b/Source/Filter/BilinearFilter.cpp @@ -70,14 +70,14 @@ BilinearFilter::ApplyStencil (MultiFab& dstmf, const MultiFab& srcmf, int scomp, #pragma omp parallel #endif { - Print()<<"stencil_x: "; for(int i=0; inGrowVect(); - ng += 1; + ng += bilinear_filter.stencil_length_each_dir-1; + //ng += 1; for (int idim = 0; idim < 3; ++idim) { Print()<<"idim "<boxArray(), @@ -392,7 +393,8 @@ WarpX::SyncCurrent () } for (int lev = 1; lev <= finest_level; ++lev) { IntVect ng = current_cp[lev][0]->nGrowVect(); - ng += 1; + ng += bilinear_filter.stencil_length_each_dir-1; + //ng += 1; for (int idim = 0; idim < 3; ++idim) { j_cp[lev][idim].reset(new MultiFab(current_cp[lev][idim]->boxArray(), current_cp[lev][idim]->DistributionMap(), @@ -405,7 +407,8 @@ WarpX::SyncCurrent () for (int lev = 1; lev <= finest_level; ++lev) { if (current_buf[lev][0]) { IntVect ng = current_buf[lev][0]->nGrowVect(); - ng += 1; + ng += bilinear_filter.stencil_length_each_dir-1; + // ng += 1; for (int idim = 0; idim < 3; ++idim) { j_buf[lev][idim].reset(new MultiFab(current_buf[lev][idim]->boxArray(), current_buf[lev][idim]->DistributionMap(), @@ -559,7 +562,8 @@ WarpX::SyncRho (amrex::Vector >& rhof, for (int lev = 0; lev <= finest_level; ++lev) { const int ncomp = rhof[lev]->nComp(); IntVect ng = rhof[lev]->nGrowVect(); - ng += 1; + //ng += 1; + ng += bilinear_filter.stencil_length_each_dir-1; rho_f_g[lev].reset(new MultiFab(rhof[lev]->boxArray(), rhof[lev]->DistributionMap(), ncomp, ng)); @@ -569,7 +573,8 @@ WarpX::SyncRho (amrex::Vector >& rhof, for (int lev = 1; lev <= finest_level; ++lev) { const int ncomp = rhoc[lev]->nComp(); IntVect ng = rhoc[lev]->nGrowVect(); - ng += 1; + //ng += 1; + ng += bilinear_filter.stencil_length_each_dir-1; rho_c_g[lev].reset(new MultiFab(rhoc[lev]->boxArray(), rhoc[lev]->DistributionMap(), ncomp, ng)); @@ -580,7 +585,8 @@ WarpX::SyncRho (amrex::Vector >& rhof, if (charge_buf[lev]) { const int ncomp = charge_buf[lev]->nComp(); IntVect ng = charge_buf[lev]->nGrowVect(); - ng += 1; + //ng += 1; + ng += bilinear_filter.stencil_length_each_dir-1; rho_buf_g[lev].reset(new MultiFab(charge_buf[lev]->boxArray(), charge_buf[lev]->DistributionMap(), ncomp, ng)); @@ -714,7 +720,8 @@ WarpX::ApplyFilterandSumBoundaryJ (int lev, PatchType patch_type) for (int idim = 0; idim < 3; ++idim) { if (use_filter) { IntVect ng = j[idim]->nGrowVect(); - ng += 1; + //ng += 1; + ng += bilinear_filter.stencil_length_each_dir-1; MultiFab jf(j[idim]->boxArray(), j[idim]->DistributionMap(), 1, ng); bilinear_filter.ApplyStencil(jf, *j[idim]); jf.SumBoundary(period); @@ -755,7 +762,8 @@ WarpX::AddCurrentFromFineLevelandSumBoundary (int lev) { // coarse patch of fine level IntVect ng = current_cp[lev+1][idim]->nGrowVect(); - ng += 1; + //ng += 1; + ng += bilinear_filter.stencil_length_each_dir-1; MultiFab jfc(current_cp[lev+1][idim]->boxArray(), current_cp[lev+1][idim]->DistributionMap(), 1, ng); bilinear_filter.ApplyStencil(jfc, *current_cp[lev+1][idim]); @@ -775,7 +783,8 @@ WarpX::AddCurrentFromFineLevelandSumBoundary (int lev) { // coarse patch of fine level IntVect ng = current_cp[lev+1][idim]->nGrowVect(); - ng += 1; + //ng += 1; + ng += bilinear_filter.stencil_length_each_dir-1; MultiFab jf(current_cp[lev+1][idim]->boxArray(), current_cp[lev+1][idim]->DistributionMap(), 1, ng); bilinear_filter.ApplyStencil(jf, *current_cp[lev+1][idim]); @@ -825,7 +834,8 @@ WarpX::ApplyFilterandSumBoundaryRho (int lev, PatchType patch_type, int icomp, i if (r == nullptr) return; if (use_filter) { IntVect ng = r->nGrowVect(); - ng += 1; + //ng += 1; + ng += bilinear_filter.stencil_length_each_dir-1; MultiFab rf(r->boxArray(), r->DistributionMap(), ncomp, ng); bilinear_filter.ApplyStencil(rf, *r, icomp, 0, ncomp); rf.SumBoundary(period); @@ -861,7 +871,8 @@ WarpX::AddRhoFromFineLevelandSumBoundary(int lev, int icomp, int ncomp) { // coarse patch of fine level IntVect ng = rho_cp[lev+1]->nGrowVect(); - ng += 1; + //ng += 1; + ng += bilinear_filter.stencil_length_each_dir-1; MultiFab rhofc(rho_cp[lev+1]->boxArray(), rho_cp[lev+1]->DistributionMap(), ncomp, ng); bilinear_filter.ApplyStencil(rhofc, *rho_cp[lev+1], icomp, 0, ncomp); @@ -880,7 +891,8 @@ WarpX::AddRhoFromFineLevelandSumBoundary(int lev, int icomp, int ncomp) else if (use_filter) // but no buffer { IntVect ng = rho_cp[lev+1]->nGrowVect(); - ng += 1; + //ng += 1; + ng += bilinear_filter.stencil_length_each_dir-1; MultiFab rf(rho_cp[lev+1]->boxArray(), rho_cp[lev+1]->DistributionMap(), ncomp, ng); bilinear_filter.ApplyStencil(rf, *rho_cp[lev+1], icomp, 0, ncomp); mf.ParallelAdd(rf, 0, 0, ncomp, ng, IntVect::TheZeroVector(), period); diff --git a/Source/WarpX.H b/Source/WarpX.H index c25323553..6f0d49e24 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -388,9 +388,6 @@ private: void LoadBalance (); - static void applyFilter (amrex::MultiFab& dstmf, const amrex::MultiFab& srcmf, - int scomp = 0, int dcomp = 0, int ncomp = 10000); - void BuildBufferMasks (); const amrex::iMultiFab* getCurrentBufferMasks (int lev) const { return current_buffer_masks[lev].get(); diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 09f2c875b..685c12902 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -955,33 +955,6 @@ WarpX::ComputeDivE (MultiFab& divE, int dcomp, } } -void -WarpX::applyFilter (MultiFab& dstmf, const MultiFab& srcmf, int scomp, int dcomp, int ncomp) -{ - 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,1); - tmpfab.resize(gbx,ncomp); - tmpfab.setVal(0.0, gbx, 0, ncomp); - const Box& ibx = gbx & srcfab.box(); - tmpfab.copy(srcfab, ibx, scomp, ibx, 0, ncomp); - WRPX_FILTER(BL_TO_FORTRAN_BOX(tbx), - BL_TO_FORTRAN_ANYD(tmpfab), - BL_TO_FORTRAN_N_ANYD(dstfab,dcomp), - ncomp); - } - } - } - void WarpX::BuildBufferMasks () { -- cgit v1.2.3 From 5ba52fe850624317746d49e2a5ca3d847213078b Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Tue, 19 Mar 2019 13:13:48 -0700 Subject: minor cleaning --- Source/Filter/BilinearFilter.cpp | 2 +- Source/Initialization/WarpXInitData.cpp | 1 - Source/Parallelization/WarpXComm.cpp | 1 + Source/WarpX.H | 1 - 4 files changed, 2 insertions(+), 3 deletions(-) (limited to 'Source/Parallelization/WarpXComm.cpp') diff --git a/Source/Filter/BilinearFilter.cpp b/Source/Filter/BilinearFilter.cpp index d4016a7dc..8f0747afb 100644 --- a/Source/Filter/BilinearFilter.cpp +++ b/Source/Filter/BilinearFilter.cpp @@ -82,7 +82,7 @@ BilinearFilter::ApplyStencil (MultiFab& dstmf, const MultiFab& srcmf, int scomp, tmpfab.setVal(0.0, gbx, 0, ncomp); const Box& ibx = gbx & srcfab.box(); tmpfab.copy(srcfab, ibx, scomp, ibx, 0, ncomp); - filter_2d(tbx, gbx, tmpfab, dstfab, dcomp, ncomp); + Filter2d(tbx, gbx, tmpfab, dstfab, dcomp, ncomp); } } } diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp index 3b9fcdb65..55df6db9c 100644 --- a/Source/Initialization/WarpXInitData.cpp +++ b/Source/Initialization/WarpXInitData.cpp @@ -39,7 +39,6 @@ WarpX::InitData () WarpX::InitNCICorrector(); } - Print()<<"HERE I AM\n"; if (WarpX::use_filter) { WarpX::InitFilter(); } diff --git a/Source/Parallelization/WarpXComm.cpp b/Source/Parallelization/WarpXComm.cpp index dd9ab7142..51a5e64de 100644 --- a/Source/Parallelization/WarpXComm.cpp +++ b/Source/Parallelization/WarpXComm.cpp @@ -387,6 +387,7 @@ WarpX::SyncCurrent () j_fp[lev][idim].reset(new MultiFab(current_fp[lev][idim]->boxArray(), current_fp[lev][idim]->DistributionMap(), 1, ng)); + // applyFilter(*j_fp[lev][idim], *current_fp[lev][idim]); bilinear_filter.ApplyStencil(*j_fp[lev][idim], *current_fp[lev][idim]); std::swap(j_fp[lev][idim], current_fp[lev][idim]); } diff --git a/Source/WarpX.H b/Source/WarpX.H index 6f0d49e24..b24d0a44a 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -139,7 +139,6 @@ public: } static amrex::IntVect filter_npass_each_dir; - // static amrex::Vector filter_npass_each_dir; BilinearFilter bilinear_filter; void ComputeDt (); -- cgit v1.2.3 From ae4ecc49ac31e7b168ccf96101f33b019feaf7fe Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Tue, 19 Mar 2019 17:46:50 -0700 Subject: final cleaning, and make sure that default value is single-pass --- Source/Filter/BilinearFilter.cpp | 8 -------- Source/Parallelization/WarpXComm.cpp | 16 ---------------- Source/WarpX.cpp | 6 +----- 3 files changed, 1 insertion(+), 29 deletions(-) (limited to 'Source/Parallelization/WarpXComm.cpp') diff --git a/Source/Filter/BilinearFilter.cpp b/Source/Filter/BilinearFilter.cpp index dd77889c0..a539d9a6e 100644 --- a/Source/Filter/BilinearFilter.cpp +++ b/Source/Filter/BilinearFilter.cpp @@ -41,7 +41,6 @@ void compute_stencil(Vector &stencil, int npass){ } void BilinearFilter::ComputeStencils(){ - Print()<<"npass_each_dir "<0!"); FArrayBox tmpfab; @@ -98,8 +94,6 @@ void BilinearFilter::Filter2d(const Box& tbx, FArrayBox &tmpfab, FArrayBox &dstf { const int* loVector = tbx.loVect(); const int* hiVector = tbx.hiVect(); - // const int* loVectorg = gbx.loVect(); - // const int* hiVectorg = gbx.hiVect(); Array4 const& tmparr = tmpfab.array(); Array4 const& dstarr = dstfab.array(); for(int i=loVector[0]; i<=hiVector[0]; i++){ @@ -118,8 +112,6 @@ void BilinearFilter::Filter3d(const Box& tbx, FArrayBox &tmpfab, FArrayBox &dstf { const int* loVector = tbx.loVect(); const int* hiVector = tbx.hiVect(); - // const int* loVectorg = gbx.loVect(); - // const int* hiVectorg = gbx.hiVect(); Array4 const& tmparr = tmpfab.array(); Array4 const& dstarr = dstfab.array(); for(int i=loVector[0]; i<=hiVector[0]; i++){ diff --git a/Source/Parallelization/WarpXComm.cpp b/Source/Parallelization/WarpXComm.cpp index 51a5e64de..38ab3a4a5 100644 --- a/Source/Parallelization/WarpXComm.cpp +++ b/Source/Parallelization/WarpXComm.cpp @@ -381,13 +381,10 @@ WarpX::SyncCurrent () for (int lev = 0; lev <= finest_level; ++lev) { IntVect ng = current_fp[lev][0]->nGrowVect(); ng += bilinear_filter.stencil_length_each_dir-1; - //ng += 1; for (int idim = 0; idim < 3; ++idim) { - Print()<<"idim "<boxArray(), current_fp[lev][idim]->DistributionMap(), 1, ng)); - // applyFilter(*j_fp[lev][idim], *current_fp[lev][idim]); bilinear_filter.ApplyStencil(*j_fp[lev][idim], *current_fp[lev][idim]); std::swap(j_fp[lev][idim], current_fp[lev][idim]); } @@ -395,12 +392,10 @@ WarpX::SyncCurrent () for (int lev = 1; lev <= finest_level; ++lev) { IntVect ng = current_cp[lev][0]->nGrowVect(); ng += bilinear_filter.stencil_length_each_dir-1; - //ng += 1; for (int idim = 0; idim < 3; ++idim) { j_cp[lev][idim].reset(new MultiFab(current_cp[lev][idim]->boxArray(), current_cp[lev][idim]->DistributionMap(), 1, ng)); - // applyFilter(*j_cp[lev][idim], *current_cp[lev][idim]); bilinear_filter.ApplyStencil(*j_cp[lev][idim], *current_cp[lev][idim]); std::swap(j_cp[lev][idim], current_cp[lev][idim]); } @@ -409,12 +404,10 @@ WarpX::SyncCurrent () if (current_buf[lev][0]) { IntVect ng = current_buf[lev][0]->nGrowVect(); ng += bilinear_filter.stencil_length_each_dir-1; - // ng += 1; for (int idim = 0; idim < 3; ++idim) { j_buf[lev][idim].reset(new MultiFab(current_buf[lev][idim]->boxArray(), current_buf[lev][idim]->DistributionMap(), 1, ng)); - //applyFilter(*j_buf[lev][idim], *current_buf[lev][idim]); bilinear_filter.ApplyStencil(*j_buf[lev][idim], *current_buf[lev][idim]); std::swap(*j_buf[lev][idim], *current_buf[lev][idim]); } @@ -563,7 +556,6 @@ WarpX::SyncRho (amrex::Vector >& rhof, for (int lev = 0; lev <= finest_level; ++lev) { const int ncomp = rhof[lev]->nComp(); IntVect ng = rhof[lev]->nGrowVect(); - //ng += 1; ng += bilinear_filter.stencil_length_each_dir-1; rho_f_g[lev].reset(new MultiFab(rhof[lev]->boxArray(), rhof[lev]->DistributionMap(), @@ -574,7 +566,6 @@ WarpX::SyncRho (amrex::Vector >& rhof, for (int lev = 1; lev <= finest_level; ++lev) { const int ncomp = rhoc[lev]->nComp(); IntVect ng = rhoc[lev]->nGrowVect(); - //ng += 1; ng += bilinear_filter.stencil_length_each_dir-1; rho_c_g[lev].reset(new MultiFab(rhoc[lev]->boxArray(), rhoc[lev]->DistributionMap(), @@ -586,7 +577,6 @@ WarpX::SyncRho (amrex::Vector >& rhof, if (charge_buf[lev]) { const int ncomp = charge_buf[lev]->nComp(); IntVect ng = charge_buf[lev]->nGrowVect(); - //ng += 1; ng += bilinear_filter.stencil_length_each_dir-1; rho_buf_g[lev].reset(new MultiFab(charge_buf[lev]->boxArray(), charge_buf[lev]->DistributionMap(), @@ -721,7 +711,6 @@ WarpX::ApplyFilterandSumBoundaryJ (int lev, PatchType patch_type) for (int idim = 0; idim < 3; ++idim) { if (use_filter) { IntVect ng = j[idim]->nGrowVect(); - //ng += 1; ng += bilinear_filter.stencil_length_each_dir-1; MultiFab jf(j[idim]->boxArray(), j[idim]->DistributionMap(), 1, ng); bilinear_filter.ApplyStencil(jf, *j[idim]); @@ -763,7 +752,6 @@ WarpX::AddCurrentFromFineLevelandSumBoundary (int lev) { // coarse patch of fine level IntVect ng = current_cp[lev+1][idim]->nGrowVect(); - //ng += 1; ng += bilinear_filter.stencil_length_each_dir-1; MultiFab jfc(current_cp[lev+1][idim]->boxArray(), current_cp[lev+1][idim]->DistributionMap(), 1, ng); @@ -784,7 +772,6 @@ WarpX::AddCurrentFromFineLevelandSumBoundary (int lev) { // coarse patch of fine level IntVect ng = current_cp[lev+1][idim]->nGrowVect(); - //ng += 1; ng += bilinear_filter.stencil_length_each_dir-1; MultiFab jf(current_cp[lev+1][idim]->boxArray(), current_cp[lev+1][idim]->DistributionMap(), 1, ng); @@ -835,7 +822,6 @@ WarpX::ApplyFilterandSumBoundaryRho (int lev, PatchType patch_type, int icomp, i if (r == nullptr) return; if (use_filter) { IntVect ng = r->nGrowVect(); - //ng += 1; ng += bilinear_filter.stencil_length_each_dir-1; MultiFab rf(r->boxArray(), r->DistributionMap(), ncomp, ng); bilinear_filter.ApplyStencil(rf, *r, icomp, 0, ncomp); @@ -872,7 +858,6 @@ WarpX::AddRhoFromFineLevelandSumBoundary(int lev, int icomp, int ncomp) { // coarse patch of fine level IntVect ng = rho_cp[lev+1]->nGrowVect(); - //ng += 1; ng += bilinear_filter.stencil_length_each_dir-1; MultiFab rhofc(rho_cp[lev+1]->boxArray(), rho_cp[lev+1]->DistributionMap(), ncomp, ng); @@ -892,7 +877,6 @@ WarpX::AddRhoFromFineLevelandSumBoundary(int lev, int icomp, int ncomp) else if (use_filter) // but no buffer { IntVect ng = rho_cp[lev+1]->nGrowVect(); - //ng += 1; ng += bilinear_filter.stencil_length_each_dir-1; MultiFab rf(rho_cp[lev+1]->boxArray(), rho_cp[lev+1]->DistributionMap(), ncomp, ng); bilinear_filter.ApplyStencil(rf, *rho_cp[lev+1], icomp, 0, ncomp); diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 914649941..6e8375aa3 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -346,7 +346,7 @@ WarpX::ReadParameters () // Read filter and fill IntVect filter_npass_each_dir with // proper size for AMREX_SPACEDIM pp.query("use_filter", use_filter); - Vector parse_filter_npass_each_dir(AMREX_SPACEDIM,0); + Vector parse_filter_npass_each_dir(AMREX_SPACEDIM,1); pp.queryarr("filter_npass_each_dir", parse_filter_npass_each_dir); filter_npass_each_dir[0] = parse_filter_npass_each_dir[0]; filter_npass_each_dir[1] = parse_filter_npass_each_dir[1]; @@ -623,10 +623,6 @@ WarpX::AllocLevelData (int lev, const BoxArray& ba, const DistributionMapping& d ngJz = std::max(ngJz,2); } - //Print()<<"ngz "< Date: Thu, 21 Mar 2019 11:06:51 -0700 Subject: fixes proposed by Weiqun --- Source/Filter/BilinearFilter.H | 3 +++ Source/Filter/BilinearFilter.cpp | 2 ++ Source/Parallelization/WarpXComm.cpp | 1 - 3 files changed, 5 insertions(+), 1 deletion(-) (limited to 'Source/Parallelization/WarpXComm.cpp') diff --git a/Source/Filter/BilinearFilter.H b/Source/Filter/BilinearFilter.H index e13079f61..5a5ee6c87 100644 --- a/Source/Filter/BilinearFilter.H +++ b/Source/Filter/BilinearFilter.H @@ -1,3 +1,6 @@ +#include +#include + #ifndef WARPX_FILTER_H_ #define WARPX_FILTER_H_ diff --git a/Source/Filter/BilinearFilter.cpp b/Source/Filter/BilinearFilter.cpp index 0bed49fab..2c9209d69 100644 --- a/Source/Filter/BilinearFilter.cpp +++ b/Source/Filter/BilinearFilter.cpp @@ -8,6 +8,7 @@ using namespace amrex; +namespace { void compute_stencil(Vector &stencil, int npass){ Vector old_s(1+npass,0.); Vector new_s(1+npass,0.); @@ -39,6 +40,7 @@ void compute_stencil(Vector &stencil, int npass){ // is corrent even when npass = 0 stencil = old_s; } +} void BilinearFilter::ComputeStencils(){ BL_PROFILE("BilinearFilter::ComputeStencils()"); diff --git a/Source/Parallelization/WarpXComm.cpp b/Source/Parallelization/WarpXComm.cpp index 38ab3a4a5..b1740df5b 100644 --- a/Source/Parallelization/WarpXComm.cpp +++ b/Source/Parallelization/WarpXComm.cpp @@ -704,7 +704,6 @@ WarpX::RestrictCurrentFromFineToCoarsePatch (int lev) void WarpX::ApplyFilterandSumBoundaryJ (int lev, PatchType patch_type) { - Print()<<"HERE AND SHOULD NOT BE\n"; const int glev = (patch_type == PatchType::fine) ? lev : lev-1; const auto& period = Geom(glev).periodicity(); auto& j = (patch_type == PatchType::fine) ? current_fp[lev] : current_cp[lev]; -- cgit v1.2.3 From d3dc85b6afc48d543f5279f0e8d80140d6a3b36a Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Thu, 21 Mar 2019 16:24:31 -0700 Subject: Added comments --- Source/Filter/BilinearFilter.cpp | 8 +++++++- Source/Parallelization/WarpXComm.cpp | 16 ++++++++++++++++ 2 files changed, 23 insertions(+), 1 deletion(-) (limited to 'Source/Parallelization/WarpXComm.cpp') diff --git a/Source/Filter/BilinearFilter.cpp b/Source/Filter/BilinearFilter.cpp index d153206ce..17857a416 100644 --- a/Source/Filter/BilinearFilter.cpp +++ b/Source/Filter/BilinearFilter.cpp @@ -84,10 +84,13 @@ BilinearFilter::ApplyStencil (MultiFab& dstmf, const MultiFab& srcmf, int scomp, 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, dstfab, 0, dcomp, ncomp); } } @@ -100,10 +103,12 @@ void BilinearFilter::Filter (const Box& tbx, FArrayBox const& tmpfab, FArrayBox const auto hi = amrex::ubound(tbx); const auto tmp = tmpfab.array(); const auto dst = dstfab.array(); + // 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) { @@ -111,7 +116,7 @@ void BilinearFilter::Filter (const Box& tbx, FArrayBox const& tmpfab, FArrayBox } } } - + // 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){ @@ -120,6 +125,7 @@ void BilinearFilter::Filter (const Box& tbx, FArrayBox const& tmpfab, FArrayBox #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 diff --git a/Source/Parallelization/WarpXComm.cpp b/Source/Parallelization/WarpXComm.cpp index b1740df5b..28971eb0c 100644 --- a/Source/Parallelization/WarpXComm.cpp +++ b/Source/Parallelization/WarpXComm.cpp @@ -382,11 +382,18 @@ WarpX::SyncCurrent () IntVect ng = current_fp[lev][0]->nGrowVect(); ng += bilinear_filter.stencil_length_each_dir-1; for (int idim = 0; idim < 3; ++idim) { + // Create new MultiFab j_jp with enough guard cells for the + // (potentially large) stencil of the multi-pass bilinear filter. j_fp[lev][idim].reset(new MultiFab(current_fp[lev][idim]->boxArray(), current_fp[lev][idim]->DistributionMap(), 1, ng)); + // Apply the filter to current_fp, store the result in j_fp. bilinear_filter.ApplyStencil(*j_fp[lev][idim], *current_fp[lev][idim]); + // Then swap j_fp and current_fp std::swap(j_fp[lev][idim], current_fp[lev][idim]); + // At this point, current_fp may have false values close to the + // edges of each FAB. This will be solved with a SumBoundary later. + // j_fp contains the exact MultiFab current_fp before this step. } } for (int lev = 1; lev <= finest_level; ++lev) { @@ -419,6 +426,9 @@ WarpX::SyncCurrent () for (int lev = 0; lev <= finest_level; ++lev) { const auto& period = Geom(lev).periodicity(); + // When using a bilinear filter with many passes, current_fp has + // temporarily more ghost cells here, so that its value inside + // the domain is correct at the end of this stage. current_fp[lev][0]->SumBoundary(period); current_fp[lev][1]->SumBoundary(period); current_fp[lev][2]->SumBoundary(period); @@ -460,8 +470,14 @@ WarpX::SyncCurrent () for (int lev = 0; lev <= finest_level; ++lev) { for (int idim = 0; idim < 3; ++idim) { + // swap j_fp and current_fp so that j_fp has correct values inside + // the domain and wrong number of ghost cells. + // current_fp has right number of ghost cells. std::swap(j_fp[lev][idim], current_fp[lev][idim]); + // Then copy the interior of j_fp to current_fp. MultiFab::Copy(*current_fp[lev][idim], *j_fp[lev][idim], 0, 0, 1, 0); + // current_fp has right number of ghost cells and + // correct filtered values here. } } for (int lev = 1; lev <= finest_level; ++lev) -- cgit v1.2.3