From 50d2c5ebaa1c7424393e17862fb4313fa1f29f47 Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Mon, 6 May 2019 17:20:57 -0700 Subject: BilinearFilter derives from class Filter --- Source/Filter/Filter.cpp | 177 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 177 insertions(+) create mode 100644 Source/Filter/Filter.cpp (limited to 'Source/Filter/Filter.cpp') diff --git a/Source/Filter/Filter.cpp b/Source/Filter/Filter.cpp new file mode 100644 index 000000000..7e13f1cc9 --- /dev/null +++ b/Source/Filter/Filter.cpp @@ -0,0 +1,177 @@ +#include +#include + +#ifdef _OPENMP +#include +#endif + +using namespace amrex; + +#ifdef AMREX_USE_CUDA + +void +Filter::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 + DoFilter(tbx, tmp, dst, 0, dcomp, ncomp); + } +} + +void Filter::DoFilter (const Box& tbx, + Array4 const& tmp, + Array4 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 +Filter::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 + DoFilter(tbx, tmpfab.array(), dstfab.array(), 0, dcomp, ncomp); + } + } +} + +void Filter::DoFilter (const Box& tbx, + Array4 const& tmp, + Array4 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 // #ifdef AMREX_USE_CUDA -- cgit v1.2.3 From bd9e80c56baf4bb6d8760befb3a42896f03fe912 Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Wed, 8 May 2019 11:04:43 -0700 Subject: Apply new C++ NCI filter. OK in 2D --- Source/Filter/Filter.H | 5 ++ Source/Filter/Filter.cpp | 69 +++++++++++++++++++++++++- Source/Filter/NCIGodfreyFilter.cpp | 2 + Source/Particles/PhysicalParticleContainer.cpp | 50 ++++++++++++++++--- 4 files changed, 119 insertions(+), 7 deletions(-) (limited to 'Source/Filter/Filter.cpp') diff --git a/Source/Filter/Filter.H b/Source/Filter/Filter.H index ef322d8e7..86218c913 100644 --- a/Source/Filter/Filter.H +++ b/Source/Filter/Filter.H @@ -14,6 +14,11 @@ public: const amrex::MultiFab& srcmf, int scomp=0, int dcomp=0, int ncomp=10000); + void ApplyStencil (amrex::FArrayBox& dstfab, + const amrex::FArrayBox& srcfab, const amrex::Box& tbx, + const amrex::Box& srcbx, + int scomp=0, int dcomp=0, int ncomp=10000); + // public for cuda void DoFilter(const amrex::Box& tbx, amrex::Array4 const& tmp, diff --git a/Source/Filter/Filter.cpp b/Source/Filter/Filter.cpp index 7e13f1cc9..ab1a6f18e 100644 --- a/Source/Filter/Filter.cpp +++ b/Source/Filter/Filter.cpp @@ -12,7 +12,7 @@ using namespace amrex; void Filter::ApplyStencil (MultiFab& dstmf, const MultiFab& srcmf, int scomp, int dcomp, int ncomp) { - BL_PROFILE("BilinearFilter::ApplyStencil()"); + BL_PROFILE("BilinearFilter::ApplyStencil(MultiFab)"); ncomp = std::min(ncomp, srcmf.nComp()); for (MFIter mfi(dstmf); mfi.isValid(); ++mfi) @@ -43,6 +43,43 @@ Filter::ApplyStencil (MultiFab& dstmf, const MultiFab& srcmf, int scomp, int dco } } +void +Filter::ApplyStencil (FArrayBox& dstfab, const FArrayBox& srcfab, const Box& tbx, + const Box& srcbx, int scomp, int dcomp, int ncomp) +{ + BL_PROFILE("BilinearFilter::ApplyStencil(FArrayBox)"); + ncomp = std::min(ncomp, srcfab.nComp()); + + // for (MFIter mfi(dstmf); mfi.isValid(); ++mfi) + //{ + const auto& src = srcfab.array(mfi); + const auto& dst = dstfab.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(); + const Box& ibx = gbx & srcbx; + + 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 + DoFilter(tbx, tmp, dst, 0, dcomp, ncomp); + //} +} + void Filter::DoFilter (const Box& tbx, Array4 const& tmp, Array4 const& dst, @@ -116,6 +153,36 @@ Filter::ApplyStencil (MultiFab& dstmf, const MultiFab& srcmf, int scomp, int dco } } +void +Filter::ApplyStencil (FArrayBox& dstfab, const FArrayBox& srcfab, const Box& tbx, const Box& srcbx, int scomp, int dcomp, int ncomp) +{ + BL_PROFILE("BilinearFilter::ApplyStencil(FArrayBox)"); + ncomp = std::min(ncomp, srcfab.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 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(); + const Box& ibx = gbx & srcbx; + tmpfab.copy(srcfab, ibx, scomp, ibx, 0, ncomp); + // Apply filter + DoFilter(tbx, tmpfab.array(), dstfab.array(), 0, dcomp, ncomp); + //} + //} +} + void Filter::DoFilter (const Box& tbx, Array4 const& tmp, Array4 const& dst, diff --git a/Source/Filter/NCIGodfreyFilter.cpp b/Source/Filter/NCIGodfreyFilter.cpp index e58cc5ac8..095cb1d8b 100644 --- a/Source/Filter/NCIGodfreyFilter.cpp +++ b/Source/Filter/NCIGodfreyFilter.cpp @@ -13,8 +13,10 @@ NCIGodfreyFilter::NCIGodfreyFilter(godfrey_coeff_set coeff_set_, amrex::Real cdt cdtodz = cdtodz_; l_lower_order_in_v = l_lower_order_in_v_; #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 } diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 4a60fc630..771a41a67 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1179,6 +1179,7 @@ PhysicalParticleContainer::Evolve (int lev, FArrayBox const* byfab = &(By[pti]); FArrayBox const* bzfab = &(Bz[pti]); + Elixir exeli, eyeli, ezeli, bxeli, byeli, bzeli; if (WarpX::use_fdtd_nci_corr) { #if (AMREX_SPACEDIM == 2) @@ -1194,7 +1195,9 @@ PhysicalParticleContainer::Evolve (int lev, // bilinear_filter.ApplyStencil(*j_fp[lev][idim], *current_fp[lev][idim]); filtered_Ex.resize(amrex::convert(tbox,WarpX::Ex_nodal_flag)); - nci_godfrey_filter_exeybz[lev]->ApplyStencil(Ex[pti], filtered_Ex) + exeli = filtered_Ex.elixir(); + //nci_godfrey_filter_exeybz[lev]->ApplyStencil(Ex[pti], filtered_Ex) + nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Ex, Ex[pti], tbox, box); /* WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Ex), BL_TO_FORTRAN_ANYD(filtered_Ex), @@ -1205,7 +1208,9 @@ PhysicalParticleContainer::Evolve (int lev, exfab = &filtered_Ex; filtered_Ez.resize(amrex::convert(tbox,WarpX::Ez_nodal_flag)); - nci_godfrey_filter_bxbyez[lev]->ApplyStencil(Ez[pti], filtered_Ez) + ezeli = filtered_Ez.elixir(); + //nci_godfrey_filter_bxbyez[lev]->ApplyStencil(Ez[pti], filtered_Ez) + nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_Ez, Ez[pti], tbox, box); /* WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Ez), BL_TO_FORTRAN_ANYD(filtered_Ez), @@ -1216,7 +1221,9 @@ PhysicalParticleContainer::Evolve (int lev, ezfab = &filtered_Ez; filtered_By.resize(amrex::convert(tbox,WarpX::By_nodal_flag)); - nci_godfrey_filter_bxbyez[lev]->ApplyStencil(By[pti], filtered_By) + byeli = filtered_By.elixir(); + //nci_godfrey_filter_bxbyez[lev]->ApplyStencil(By[pti], filtered_By) + nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_By, By[pti], tbox, box); /* WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_By), BL_TO_FORTRAN_ANYD(filtered_By), @@ -1228,7 +1235,9 @@ PhysicalParticleContainer::Evolve (int lev, #if (AMREX_SPACEDIM == 3) filtered_Ey.resize(amrex::convert(tbox,WarpX::Ey_nodal_flag)); - nci_godfrey_filter_exeybz[lev]->ApplyStencil(Ey[pti], filtered_Ey) + eyeli = filtered_Ey.elixir(); + //nci_godfrey_filter_exeybz[lev]->ApplyStencil(Ey[pti], filtered_Ey) + nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Ey, Ey[pti], tbox, box); /* WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Ey), BL_TO_FORTRAN_ANYD(filtered_Ey), @@ -1239,7 +1248,9 @@ PhysicalParticleContainer::Evolve (int lev, eyfab = &filtered_Ey; filtered_Bx.resize(amrex::convert(tbox,WarpX::Bx_nodal_flag)); - nci_godfrey_filter_bxbyez[lev]->ApplyStencil(Bx[pti], filtered_Bx) + bxeli = filtered_Bx.elixir(); + //nci_godfrey_filter_bxbyez[lev]->ApplyStencil(Bx[pti], filtered_Bx) + nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_Bx, Bx[pti], tbox, box); /* WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Bx), BL_TO_FORTRAN_ANYD(filtered_Bx), @@ -1250,7 +1261,9 @@ PhysicalParticleContainer::Evolve (int lev, bxfab = &filtered_Bx; filtered_Bz.resize(amrex::convert(tbox,WarpX::Bz_nodal_flag)); - nci_godfrey_filter_exeybz[lev]->ApplyStencil(Bz[pti], filtered_Bz) + bzeli = filtered_Bz.elixir(); + //nci_godfrey_filter_exeybz[lev]->ApplyStencil(Bz[pti], filtered_Bz) + nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Bz, Bz[pti], tbox, box); /* WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Bz), BL_TO_FORTRAN_ANYD(filtered_Bz), @@ -1422,6 +1435,7 @@ PhysicalParticleContainer::Evolve (int lev, const FArrayBox* cbyfab = &(*cBy)[pti]; const FArrayBox* cbzfab = &(*cBz)[pti]; + Elixir exeli, eyeli, ezeli, bxeli, byeli, bzeli; if (WarpX::use_fdtd_nci_corr) { #if (AMREX_SPACEDIM == 2) @@ -1435,51 +1449,75 @@ PhysicalParticleContainer::Evolve (int lev, // both 2d and 3d filtered_Ex.resize(amrex::convert(tbox,WarpX::Ex_nodal_flag)); + exeli = filtered_Ex.elixir(); + nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Ex, Ex[pti], tbox, cbox); + /* WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Ex), BL_TO_FORTRAN_ANYD(filtered_Ex), BL_TO_FORTRAN_ANYD((*cEx)[pti]), mypc.fdtd_nci_stencilz_ex[lev-1].data(), &nstencilz_fdtd_nci_corr); + */ cexfab = &filtered_Ex; filtered_Ez.resize(amrex::convert(tbox,WarpX::Ez_nodal_flag)); + ezeli = filtered_Ez.elixir(); + nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_Ez, Ez[pti], tbox, cbox); + /* WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Ez), BL_TO_FORTRAN_ANYD(filtered_Ez), BL_TO_FORTRAN_ANYD((*cEz)[pti]), mypc.fdtd_nci_stencilz_by[lev-1].data(), &nstencilz_fdtd_nci_corr); cezfab = &filtered_Ez; + */ filtered_By.resize(amrex::convert(tbox,WarpX::By_nodal_flag)); + byeli = filtered_By.elixir(); + nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_By, By[pti], tbox, cbox); + /* WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_By), BL_TO_FORTRAN_ANYD(filtered_By), BL_TO_FORTRAN_ANYD((*cBy)[pti]), mypc.fdtd_nci_stencilz_by[lev-1].data(), &nstencilz_fdtd_nci_corr); + */ cbyfab = &filtered_By; #if (AMREX_SPACEDIM == 3) filtered_Ey.resize(amrex::convert(tbox,WarpX::Ey_nodal_flag)); + eyeli = filtered_Ey.elixir(); + nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Ey, Ey[pti], tbox, cbox); + /* WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Ey), BL_TO_FORTRAN_ANYD(filtered_Ey), BL_TO_FORTRAN_ANYD((*cEy)[pti]), mypc.fdtd_nci_stencilz_ex[lev-1].data(), &nstencilz_fdtd_nci_corr); + */ ceyfab = &filtered_Ey; filtered_Bx.resize(amrex::convert(tbox,WarpX::Bx_nodal_flag)); + bxeli = filtered_Bx.elixir(); + nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_Bx, Bx[pti], tbox, cbox); + /* WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Bx), BL_TO_FORTRAN_ANYD(filtered_Bx), BL_TO_FORTRAN_ANYD((*cBx)[pti]), mypc.fdtd_nci_stencilz_by[lev-1].data(), &nstencilz_fdtd_nci_corr); + */ cbxfab = &filtered_Bx; filtered_Bz.resize(amrex::convert(tbox,WarpX::Bz_nodal_flag)); + bzeli = filtered_Bz.elixir(); + nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Bz, Bz[pti], tbox, cbox); + /* WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Bz), BL_TO_FORTRAN_ANYD(filtered_Bz), BL_TO_FORTRAN_ANYD((*cBz)[pti]), mypc.fdtd_nci_stencilz_ex[lev-1].data(), &nstencilz_fdtd_nci_corr); + */ cbzfab = &filtered_Bz; #endif } -- cgit v1.2.3 From e41e5eac1bccaa5c3732154aa052aee6c32f6518 Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Wed, 8 May 2019 15:24:22 -0700 Subject: fixes, and add calls to the filter --- Source/Filter/Filter.H | 4 +- Source/Filter/Filter.cpp | 70 +++++------- Source/Filter/NCIGodfreyFilter.H | 11 +- Source/Filter/NCIGodfreyFilter.cpp | 11 +- Source/Particles/PhysicalParticleContainer.cpp | 141 +++---------------------- Source/Utils/NCIGodfreyTables.H | 4 +- Source/WarpX.H | 2 +- 7 files changed, 57 insertions(+), 186 deletions(-) (limited to 'Source/Filter/Filter.cpp') diff --git a/Source/Filter/Filter.H b/Source/Filter/Filter.H index 86218c913..d728a16e6 100644 --- a/Source/Filter/Filter.H +++ b/Source/Filter/Filter.H @@ -10,13 +10,13 @@ public: Filter () = default; virtual void ComputeStencils() = 0; + void ApplyStencil(amrex::MultiFab& dstmf, const amrex::MultiFab& srcmf, int scomp=0, int dcomp=0, int ncomp=10000); void ApplyStencil (amrex::FArrayBox& dstfab, const amrex::FArrayBox& srcfab, const amrex::Box& tbx, - const amrex::Box& srcbx, int scomp=0, int dcomp=0, int ncomp=10000); // public for cuda @@ -28,8 +28,8 @@ public: amrex::IntVect stencil_length_each_dir; protected: - amrex::Dim3 slen; amrex::Gpu::ManagedVector stencil_x, stencil_y, stencil_z; + amrex::Dim3 slen; private: diff --git a/Source/Filter/Filter.cpp b/Source/Filter/Filter.cpp index ab1a6f18e..9bbd230a1 100644 --- a/Source/Filter/Filter.cpp +++ b/Source/Filter/Filter.cpp @@ -44,29 +44,22 @@ Filter::ApplyStencil (MultiFab& dstmf, const MultiFab& srcmf, int scomp, int dco } void -Filter::ApplyStencil (FArrayBox& dstfab, const FArrayBox& srcfab, const Box& tbx, - const Box& srcbx, int scomp, int dcomp, int ncomp) +Filter::ApplyStencil (FArrayBox& dstfab, const FArrayBox& srcfab, + const Box& tbx, int scomp, int dcomp, int ncomp) { BL_PROFILE("BilinearFilter::ApplyStencil(FArrayBox)"); ncomp = std::min(ncomp, srcfab.nComp()); + const auto& src = srcfab.array(mfi); + const auto& dst = dstfab.array(mfi); + const Box& gbx = amrex::grow(tbx,stencil_length_each_dir-1); - // for (MFIter mfi(dstmf); mfi.isValid(); ++mfi) - //{ - const auto& src = srcfab.array(mfi); - const auto& dst = dstfab.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(); + // 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(); - const Box& ibx = gbx & srcbx; - - AMREX_PARALLEL_FOR_4D ( gbx, ncomp, i, j, k, n, + // Copy values in srcfab into tmpfab + 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); @@ -75,9 +68,8 @@ Filter::ApplyStencil (FArrayBox& dstfab, const FArrayBox& srcfab, const Box& tbx } }); - // Apply filter - DoFilter(tbx, tmp, dst, 0, dcomp, ncomp); - //} + // Apply filter + DoFilter(tbx, tmp, dst, 0, dcomp, ncomp); } void Filter::DoFilter (const Box& tbx, @@ -154,33 +146,21 @@ Filter::ApplyStencil (MultiFab& dstmf, const MultiFab& srcmf, int scomp, int dco } void -Filter::ApplyStencil (FArrayBox& dstfab, const FArrayBox& srcfab, const Box& tbx, const Box& srcbx, int scomp, int dcomp, int ncomp) +Filter::ApplyStencil (FArrayBox& dstfab, const FArrayBox& srcfab, + const Box& tbx, int scomp, int dcomp, int ncomp) { BL_PROFILE("BilinearFilter::ApplyStencil(FArrayBox)"); ncomp = std::min(ncomp, srcfab.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 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(); - const Box& ibx = gbx & srcbx; - tmpfab.copy(srcfab, ibx, scomp, ibx, 0, ncomp); - // Apply filter - DoFilter(tbx, tmpfab.array(), dstfab.array(), 0, dcomp, ncomp); - //} - //} + FArrayBox tmpfab; + 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; + tmpfab.copy(srcfab, ibx, scomp, ibx, 0, ncomp); + // Apply filter + DoFilter(tbx, tmpfab.array(), dstfab.array(), 0, dcomp, ncomp); } void Filter::DoFilter (const Box& tbx, diff --git a/Source/Filter/NCIGodfreyFilter.H b/Source/Filter/NCIGodfreyFilter.H index 2e6e65772..21e16ff74 100644 --- a/Source/Filter/NCIGodfreyFilter.H +++ b/Source/Filter/NCIGodfreyFilter.H @@ -9,15 +9,9 @@ class NCIGodfreyFilter : public Filter { public: + NCIGodfreyFilter () = default; + NCIGodfreyFilter(godfrey_coeff_set coeff_set_, amrex::Real cdtodz_, amrex::Real l_lower_order_in_v); - /* - { - : coeff_set( coeff_set_ ), - cdtodz( cdtodz_ ), - l_lower_order_in_v( l_lower_order_in_v_ ), - slen - }; - */ virtual void ComputeStencils() override final; @@ -28,7 +22,6 @@ private: godfrey_coeff_set coeff_set; amrex::Real cdtodz; int l_lower_order_in_v; - // amrex::Real* coeff_table; }; diff --git a/Source/Filter/NCIGodfreyFilter.cpp b/Source/Filter/NCIGodfreyFilter.cpp index 095cb1d8b..bab1cc33a 100644 --- a/Source/Filter/NCIGodfreyFilter.cpp +++ b/Source/Filter/NCIGodfreyFilter.cpp @@ -26,6 +26,7 @@ void NCIGodfreyFilter::ComputeStencils(){ Print()<<"slen.x "<(WarpX::noz)}); #endif - // both 2d and 3d - // bilinear_filter.ApplyStencil(*j_fp[lev][idim], *current_fp[lev][idim]); - filtered_Ex.resize(amrex::convert(tbox,WarpX::Ex_nodal_flag)); exeli = filtered_Ex.elixir(); - //nci_godfrey_filter_exeybz[lev]->ApplyStencil(Ex[pti], filtered_Ex) - nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Ex, Ex[pti], tbox, box); - /* - WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Ex), - BL_TO_FORTRAN_ANYD(filtered_Ex), - BL_TO_FORTRAN_ANYD(Ex[pti]), - mypc.fdtd_nci_stencilz_ex[lev].data(), - &nstencilz_fdtd_nci_corr); - */ + nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Ex, Ex[pti], tbox); exfab = &filtered_Ex; filtered_Ez.resize(amrex::convert(tbox,WarpX::Ez_nodal_flag)); ezeli = filtered_Ez.elixir(); - //nci_godfrey_filter_bxbyez[lev]->ApplyStencil(Ez[pti], filtered_Ez) - nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_Ez, Ez[pti], tbox, box); - /* - WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Ez), - BL_TO_FORTRAN_ANYD(filtered_Ez), - BL_TO_FORTRAN_ANYD(Ez[pti]), - mypc.fdtd_nci_stencilz_by[lev].data(), - &nstencilz_fdtd_nci_corr); - */ + nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_Ez, Ez[pti], tbox); ezfab = &filtered_Ez; filtered_By.resize(amrex::convert(tbox,WarpX::By_nodal_flag)); byeli = filtered_By.elixir(); - //nci_godfrey_filter_bxbyez[lev]->ApplyStencil(By[pti], filtered_By) - nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_By, By[pti], tbox, box); - /* - WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_By), - BL_TO_FORTRAN_ANYD(filtered_By), - BL_TO_FORTRAN_ANYD(By[pti]), - mypc.fdtd_nci_stencilz_by[lev].data(), - &nstencilz_fdtd_nci_corr); - */ + nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_By, By[pti], tbox); byfab = &filtered_By; - #if (AMREX_SPACEDIM == 3) filtered_Ey.resize(amrex::convert(tbox,WarpX::Ey_nodal_flag)); eyeli = filtered_Ey.elixir(); - //nci_godfrey_filter_exeybz[lev]->ApplyStencil(Ey[pti], filtered_Ey) - nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Ey, Ey[pti], tbox, box); - /* - WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Ey), - BL_TO_FORTRAN_ANYD(filtered_Ey), - BL_TO_FORTRAN_ANYD(Ey[pti]), - mypc.fdtd_nci_stencilz_ex[lev].data(), - &nstencilz_fdtd_nci_corr); - */ + nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Ey, Ey[pti], tbox); eyfab = &filtered_Ey; filtered_Bx.resize(amrex::convert(tbox,WarpX::Bx_nodal_flag)); bxeli = filtered_Bx.elixir(); - //nci_godfrey_filter_bxbyez[lev]->ApplyStencil(Bx[pti], filtered_Bx) - nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_Bx, Bx[pti], tbox, box); - /* - WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Bx), - BL_TO_FORTRAN_ANYD(filtered_Bx), - BL_TO_FORTRAN_ANYD(Bx[pti]), - mypc.fdtd_nci_stencilz_by[lev].data(), - &nstencilz_fdtd_nci_corr); - */ + nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_Bx, Bx[pti], tbox); bxfab = &filtered_Bx; filtered_Bz.resize(amrex::convert(tbox,WarpX::Bz_nodal_flag)); bzeli = filtered_Bz.elixir(); - //nci_godfrey_filter_exeybz[lev]->ApplyStencil(Bz[pti], filtered_Bz) - nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Bz, Bz[pti], tbox, box); - /* - WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Bz), - BL_TO_FORTRAN_ANYD(filtered_Bz), - BL_TO_FORTRAN_ANYD(Bz[pti]), - mypc.fdtd_nci_stencilz_ex[lev].data(), - &nstencilz_fdtd_nci_corr); - */ + nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Bz, Bz[pti], tbox); bzfab = &filtered_Bz; #endif } @@ -1451,74 +1382,32 @@ PhysicalParticleContainer::Evolve (int lev, // both 2d and 3d filtered_Ex.resize(amrex::convert(tbox,WarpX::Ex_nodal_flag)); exeli = filtered_Ex.elixir(); - nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Ex, Ex[pti], tbox, cbox); - /* - WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Ex), - BL_TO_FORTRAN_ANYD(filtered_Ex), - BL_TO_FORTRAN_ANYD((*cEx)[pti]), - mypc.fdtd_nci_stencilz_ex[lev-1].data(), - &nstencilz_fdtd_nci_corr); - */ + nci_godfrey_filter_exeybz[lev-1]->ApplyStencil(filtered_Ex, (*cEx)[pti], tbox); cexfab = &filtered_Ex; filtered_Ez.resize(amrex::convert(tbox,WarpX::Ez_nodal_flag)); ezeli = filtered_Ez.elixir(); - nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_Ez, Ez[pti], tbox, cbox); - /* - WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Ez), - BL_TO_FORTRAN_ANYD(filtered_Ez), - BL_TO_FORTRAN_ANYD((*cEz)[pti]), - mypc.fdtd_nci_stencilz_by[lev-1].data(), - &nstencilz_fdtd_nci_corr); + nci_godfrey_filter_bxbyez[lev-1]->ApplyStencil(filtered_Ez, (*cEz)[pti], tbox); cezfab = &filtered_Ez; - */ + filtered_By.resize(amrex::convert(tbox,WarpX::By_nodal_flag)); byeli = filtered_By.elixir(); - nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_By, By[pti], tbox, cbox); - /* - WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_By), - BL_TO_FORTRAN_ANYD(filtered_By), - BL_TO_FORTRAN_ANYD((*cBy)[pti]), - mypc.fdtd_nci_stencilz_by[lev-1].data(), - &nstencilz_fdtd_nci_corr); - */ + nci_godfrey_filter_bxbyez[lev-1]->ApplyStencil(filtered_By, (*cBy)[pti], tbox); cbyfab = &filtered_By; - #if (AMREX_SPACEDIM == 3) filtered_Ey.resize(amrex::convert(tbox,WarpX::Ey_nodal_flag)); eyeli = filtered_Ey.elixir(); - nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Ey, Ey[pti], tbox, cbox); - /* - WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Ey), - BL_TO_FORTRAN_ANYD(filtered_Ey), - BL_TO_FORTRAN_ANYD((*cEy)[pti]), - mypc.fdtd_nci_stencilz_ex[lev-1].data(), - &nstencilz_fdtd_nci_corr); - */ + nci_godfrey_filter_exeybz[lev-1]->ApplyStencil(filtered_Ey, (*cEy)[pti], tbox); ceyfab = &filtered_Ey; filtered_Bx.resize(amrex::convert(tbox,WarpX::Bx_nodal_flag)); bxeli = filtered_Bx.elixir(); - nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_Bx, Bx[pti], tbox, cbox); - /* - WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Bx), - BL_TO_FORTRAN_ANYD(filtered_Bx), - BL_TO_FORTRAN_ANYD((*cBx)[pti]), - mypc.fdtd_nci_stencilz_by[lev-1].data(), - &nstencilz_fdtd_nci_corr); - */ + nci_godfrey_filter_bxbyez[lev-1]->ApplyStencil(filtered_Bx, (*cBx)[pti], tbox); cbxfab = &filtered_Bx; filtered_Bz.resize(amrex::convert(tbox,WarpX::Bz_nodal_flag)); bzeli = filtered_Bz.elixir(); - nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Bz, Bz[pti], tbox, cbox); - /* - WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Bz), - BL_TO_FORTRAN_ANYD(filtered_Bz), - BL_TO_FORTRAN_ANYD((*cBz)[pti]), - mypc.fdtd_nci_stencilz_ex[lev-1].data(), - &nstencilz_fdtd_nci_corr); - */ + nci_godfrey_filter_exeybz[lev-1]->ApplyStencil(filtered_Bz, (*cBz)[pti], tbox); cbzfab = &filtered_Bz; #endif } diff --git a/Source/Utils/NCIGodfreyTables.H b/Source/Utils/NCIGodfreyTables.H index c181aead9..b84f2cedb 100644 --- a/Source/Utils/NCIGodfreyTables.H +++ b/Source/Utils/NCIGodfreyTables.H @@ -4,9 +4,10 @@ #define WARPX_GODFREY_COEFF_TABLE_H_ const int tab_width = 4; -const int tab_length = 100; +const int tab_length = 101; const amrex::Real table_nci_godfrey_Ex_Ey_Bz[tab_length][tab_width]{ + -2.47536,2.04288,-0.598163,0.0314711, -2.47536,2.04288,-0.598163,0.0314711, -2.47545,2.04309,-0.598307,0.0315029, -2.4756,2.04342,-0.598549,0.0315558, @@ -110,6 +111,7 @@ const amrex::Real table_nci_godfrey_Ex_Ey_Bz[tab_length][tab_width]{ }; const amrex::Real table_nci_godfrey_Bx_By_Ez[tab_length][tab_width]{ + -2.80862,2.80104,-1.14615,0.154077, -2.80862,2.80104,-1.14615,0.154077, -2.80851,2.80078,-1.14595,0.154027, -2.80832,2.80034,-1.14561,0.153945, diff --git a/Source/WarpX.H b/Source/WarpX.H index 77bc47a7d..9ac64c6c6 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -146,7 +146,7 @@ public: BilinearFilter bilinear_filter; amrex::Vector< std::unique_ptr > nci_godfrey_filter_exeybz; amrex::Vector< std::unique_ptr > nci_godfrey_filter_bxbyez; - + static int num_mirrors; amrex::Vector mirror_z; amrex::Vector mirror_z_width; -- cgit v1.2.3 From 9d25c2b7c4aa473aa088a80f636086ed249b3a7a Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Wed, 8 May 2019 16:14:49 -0700 Subject: Add comments --- .../laser_acceleration/inputs.2d.boost | 14 +++++----- Source/Filter/Filter.H | 9 ++++++ Source/Filter/Filter.cpp | 32 ++++++++++++++++++++++ Source/Filter/NCIGodfreyFilter.cpp | 25 +++++++++++------ Source/Initialization/WarpXInitData.cpp | 3 +- Source/Particles/PhysicalParticleContainer.cpp | 22 +++++++++++++-- Source/Utils/NCIGodfreyTables.H | 7 +++++ Source/WarpX.cpp | 2 ++ 8 files changed, 96 insertions(+), 18 deletions(-) (limited to 'Source/Filter/Filter.cpp') diff --git a/Examples/Physics_applications/laser_acceleration/inputs.2d.boost b/Examples/Physics_applications/laser_acceleration/inputs.2d.boost index b1fb28126..b8f0fff57 100644 --- a/Examples/Physics_applications/laser_acceleration/inputs.2d.boost +++ b/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 @@ -61,11 +61,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 @@ -76,11 +76,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 diff --git a/Source/Filter/Filter.H b/Source/Filter/Filter.H index d728a16e6..cf0848ee9 100644 --- a/Source/Filter/Filter.H +++ b/Source/Filter/Filter.H @@ -9,12 +9,16 @@ class Filter public: Filter () = default; + // This function has to be overriden by derived classes. virtual void ComputeStencils() = 0; + // 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); @@ -25,10 +29,15 @@ public: amrex::Array4 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 stencil_x, stencil_y, stencil_z; + // Length of each stencil. + // In 2D, slen = {length(stencil_x), length(stencil_z), 1} amrex::Dim3 slen; private: diff --git a/Source/Filter/Filter.cpp b/Source/Filter/Filter.cpp index 9bbd230a1..25d935ef6 100644 --- a/Source/Filter/Filter.cpp +++ b/Source/Filter/Filter.cpp @@ -9,6 +9,13 @@ using namespace amrex; #ifdef AMREX_USE_CUDA +/* \brief Apply stencil on MultiFab (GPU version, 2D/3D). + * \param dstmf Destination MultiFab + * \param srcmf source MultiFab + * \param scomp first component of srcmf on which the filter is applied + * \param dcomp first component of dstmf on which the filter is applied + * \param ncomp Number of components on which the filter is applied. + */ void Filter::ApplyStencil (MultiFab& dstmf, const MultiFab& srcmf, int scomp, int dcomp, int ncomp) { @@ -43,6 +50,14 @@ Filter::ApplyStencil (MultiFab& dstmf, const MultiFab& srcmf, int scomp, int dco } } +/* \brief Apply stencil on FArrayBox (GPU version, 2D/3D). + * \param dstfab Destination FArrayBox + * \param srcmf source FArrayBox + * \param tbx Grown box on which srcfab is defined. + * \param scomp first component of srcfab on which the filter is applied + * \param dcomp first component of dstfab on which the filter is applied + * \param ncomp Number of components on which the filter is applied. + */ void Filter::ApplyStencil (FArrayBox& dstfab, const FArrayBox& srcfab, const Box& tbx, int scomp, int dcomp, int ncomp) @@ -72,6 +87,8 @@ Filter::ApplyStencil (FArrayBox& dstfab, const FArrayBox& srcfab, DoFilter(tbx, tmp, dst, 0, dcomp, ncomp); } +/* \brief Apply stencil (2D/3D, CPU/GPU) + */ void Filter::DoFilter (const Box& tbx, Array4 const& tmp, Array4 const& dst, @@ -118,6 +135,13 @@ void Filter::DoFilter (const Box& tbx, #else +/* \brief Apply stencil on MultiFab (CPU version, 2D/3D). + * \param dstmf Destination MultiFab + * \param srcmf source MultiFab + * \param scomp first component of srcmf on which the filter is applied + * \param dcomp first component of dstmf on which the filter is applied + * \param ncomp Number of components on which the filter is applied. + */ void Filter::ApplyStencil (MultiFab& dstmf, const MultiFab& srcmf, int scomp, int dcomp, int ncomp) { @@ -145,6 +169,14 @@ Filter::ApplyStencil (MultiFab& dstmf, const MultiFab& srcmf, int scomp, int dco } } +/* \brief Apply stencil on FArrayBox (CPU version, 2D/3D). + * \param dstfab Destination FArrayBox + * \param srcmf source FArrayBox + * \param tbx Grown box on which srcfab is defined. + * \param scomp first component of srcfab on which the filter is applied + * \param dcomp first component of dstfab on which the filter is applied + * \param ncomp Number of components on which the filter is applied. + */ void Filter::ApplyStencil (FArrayBox& dstfab, const FArrayBox& srcfab, const Box& tbx, int scomp, int dcomp, int ncomp) diff --git a/Source/Filter/NCIGodfreyFilter.cpp b/Source/Filter/NCIGodfreyFilter.cpp index bab1cc33a..34fca7604 100644 --- a/Source/Filter/NCIGodfreyFilter.cpp +++ b/Source/Filter/NCIGodfreyFilter.cpp @@ -9,9 +9,11 @@ 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}; @@ -22,12 +24,6 @@ NCIGodfreyFilter::NCIGodfreyFilter(godfrey_coeff_set coeff_set_, amrex::Real cdt } void NCIGodfreyFilter::ComputeStencils(){ - Print()<<"slen "<ComputeStencils(); nci_godfrey_filter_bxbyez[lev]->ComputeStencils(); diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 0e003927d..1e0d68800 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1105,6 +1105,7 @@ PhysicalParticleContainer::Evolve (int lev, const auto& mypc = WarpX::GetInstance().GetPartContainer(); const int nstencilz_fdtd_nci_corr = mypc.nstencilz_fdtd_nci_corr; + // Get instances of NCI Godfrey filters const auto& nci_godfrey_filter_exeybz = WarpX::GetInstance().nci_godfrey_filter_exeybz; const auto& nci_godfrey_filter_bxbyez = WarpX::GetInstance().nci_godfrey_filter_bxbyez; @@ -1175,31 +1176,40 @@ PhysicalParticleContainer::Evolve (int lev, static_cast(WarpX::noz)}); #endif + // Filter Ex (Both 2D and 3D) filtered_Ex.resize(amrex::convert(tbox,WarpX::Ex_nodal_flag)); + // Safeguard for GPU exeli = filtered_Ex.elixir(); + // Apply filter on Ex, result stored in filtered_Ex nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Ex, Ex[pti], tbox); + // Update exfab reference exfab = &filtered_Ex; + // Filter Ez filtered_Ez.resize(amrex::convert(tbox,WarpX::Ez_nodal_flag)); ezeli = filtered_Ez.elixir(); nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_Ez, Ez[pti], tbox); ezfab = &filtered_Ez; + // Filter By filtered_By.resize(amrex::convert(tbox,WarpX::By_nodal_flag)); byeli = filtered_By.elixir(); nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_By, By[pti], tbox); byfab = &filtered_By; #if (AMREX_SPACEDIM == 3) + // Filter Ey filtered_Ey.resize(amrex::convert(tbox,WarpX::Ey_nodal_flag)); eyeli = filtered_Ey.elixir(); nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Ey, Ey[pti], tbox); eyfab = &filtered_Ey; + // Filter Bx filtered_Bx.resize(amrex::convert(tbox,WarpX::Bx_nodal_flag)); bxeli = filtered_Bx.elixir(); nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_Bx, Bx[pti], tbox); bxfab = &filtered_Bx; + // Filter Bz filtered_Bz.resize(amrex::convert(tbox,WarpX::Bz_nodal_flag)); bzeli = filtered_Bz.elixir(); nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Bz, Bz[pti], tbox); @@ -1378,33 +1388,41 @@ PhysicalParticleContainer::Evolve (int lev, static_cast(WarpX::noy), static_cast(WarpX::noz)}); #endif - - // both 2d and 3d + + // Filter Ex (both 2D and 3D) filtered_Ex.resize(amrex::convert(tbox,WarpX::Ex_nodal_flag)); + // Safeguard for GPU exeli = filtered_Ex.elixir(); + // Apply filter on Ex, result stored in filtered_Ex nci_godfrey_filter_exeybz[lev-1]->ApplyStencil(filtered_Ex, (*cEx)[pti], tbox); + // Update exfab reference cexfab = &filtered_Ex; + // Filter Ez filtered_Ez.resize(amrex::convert(tbox,WarpX::Ez_nodal_flag)); ezeli = filtered_Ez.elixir(); nci_godfrey_filter_bxbyez[lev-1]->ApplyStencil(filtered_Ez, (*cEz)[pti], tbox); cezfab = &filtered_Ez; + // Filter By filtered_By.resize(amrex::convert(tbox,WarpX::By_nodal_flag)); byeli = filtered_By.elixir(); nci_godfrey_filter_bxbyez[lev-1]->ApplyStencil(filtered_By, (*cBy)[pti], tbox); cbyfab = &filtered_By; #if (AMREX_SPACEDIM == 3) + // Filter Ey filtered_Ey.resize(amrex::convert(tbox,WarpX::Ey_nodal_flag)); eyeli = filtered_Ey.elixir(); nci_godfrey_filter_exeybz[lev-1]->ApplyStencil(filtered_Ey, (*cEy)[pti], tbox); ceyfab = &filtered_Ey; + // Filter Bx filtered_Bx.resize(amrex::convert(tbox,WarpX::Bx_nodal_flag)); bxeli = filtered_Bx.elixir(); nci_godfrey_filter_bxbyez[lev-1]->ApplyStencil(filtered_Bx, (*cBx)[pti], tbox); cbxfab = &filtered_Bx; + // Filter Bz filtered_Bz.resize(amrex::convert(tbox,WarpX::Bz_nodal_flag)); bzeli = filtered_Bz.elixir(); nci_godfrey_filter_exeybz[lev-1]->ApplyStencil(filtered_Bz, (*cBz)[pti], tbox); diff --git a/Source/Utils/NCIGodfreyTables.H b/Source/Utils/NCIGodfreyTables.H index b84f2cedb..8cb105aa0 100644 --- a/Source/Utils/NCIGodfreyTables.H +++ b/Source/Utils/NCIGodfreyTables.H @@ -3,9 +3,14 @@ #ifndef WARPX_GODFREY_COEFF_TABLE_H_ #define WARPX_GODFREY_COEFF_TABLE_H_ +// Table width. This is related to the stencil length const int tab_width = 4; +// table length. Each line correspond to 1 value of cdtodz +// (here 101 values). const int tab_length = 101; +// Table of coefficient for Ex, Ey abd Bz +// We typically interpolate between two lines const amrex::Real table_nci_godfrey_Ex_Ey_Bz[tab_length][tab_width]{ -2.47536,2.04288,-0.598163,0.0314711, -2.47536,2.04288,-0.598163,0.0314711, @@ -110,6 +115,8 @@ const amrex::Real table_nci_godfrey_Ex_Ey_Bz[tab_length][tab_width]{ -2.98769,3.19374,-1.41463,0.208607 }; +// Table of coefficient for Bx, By and Ez +// We typically interpolate between two lines const amrex::Real table_nci_godfrey_Bx_By_Ez[tab_length][tab_width]{ -2.80862,2.80104,-1.14615,0.154077, -2.80862,2.80104,-1.14615,0.154077, diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 36dea14b1..5f63e2ad6 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -225,6 +225,8 @@ WarpX::WarpX () insitu_bridge = nullptr; #endif + // NCI Godfrey filters can have different stencils + // at different levels (the stencil depends on c*dt/dz) nci_godfrey_filter_exeybz.resize(nlevs_max); nci_godfrey_filter_bxbyez.resize(nlevs_max); } -- cgit v1.2.3 From 7789d99e979bf27e146e4b9baee263847c8dc7a1 Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Fri, 17 May 2019 17:40:34 -0700 Subject: Fixes suggested by @WeiqunZhang --- Source/Filter/BilinearFilter.H | 4 ++-- Source/Filter/Filter.H | 3 --- Source/Filter/Filter.cpp | 3 ++- Source/Filter/NCIGodfreyFilter.H | 2 +- Source/Filter/NCIGodfreyFilter.cpp | 22 +++++++--------------- 5 files changed, 12 insertions(+), 22 deletions(-) (limited to 'Source/Filter/Filter.cpp') diff --git a/Source/Filter/BilinearFilter.H b/Source/Filter/BilinearFilter.H index ceee56e6a..4fcfec14b 100644 --- a/Source/Filter/BilinearFilter.H +++ b/Source/Filter/BilinearFilter.H @@ -8,8 +8,8 @@ class BilinearFilter : public Filter public: BilinearFilter() = default; - - virtual void ComputeStencils() override final; + + void ComputeStencils(); amrex::IntVect npass_each_dir; diff --git a/Source/Filter/Filter.H b/Source/Filter/Filter.H index cf0848ee9..eaa0498c9 100644 --- a/Source/Filter/Filter.H +++ b/Source/Filter/Filter.H @@ -9,9 +9,6 @@ class Filter public: Filter () = default; - // This function has to be overriden by derived classes. - virtual void ComputeStencils() = 0; - // Apply stencil on MultiFab. // Guard cells are handled inside this function void ApplyStencil(amrex::MultiFab& dstmf, diff --git a/Source/Filter/Filter.cpp b/Source/Filter/Filter.cpp index 25d935ef6..0c1a98443 100644 --- a/Source/Filter/Filter.cpp +++ b/Source/Filter/Filter.cpp @@ -74,6 +74,7 @@ Filter::ApplyStencil (FArrayBox& dstfab, const FArrayBox& srcfab, 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)))) { @@ -189,7 +190,7 @@ Filter::ApplyStencil (FArrayBox& dstfab, const FArrayBox& srcfab, tmpfab.resize(gbx,ncomp); tmpfab.setVal(0.0, gbx, 0, ncomp); // Copy values in srcfab into tmpfab - const Box& ibx = gbx; + const Box& ibx = gbx & srcfab.box(); tmpfab.copy(srcfab, ibx, scomp, ibx, 0, ncomp); // Apply filter DoFilter(tbx, tmpfab.array(), dstfab.array(), 0, dcomp, ncomp); diff --git a/Source/Filter/NCIGodfreyFilter.H b/Source/Filter/NCIGodfreyFilter.H index d75cb6778..b2d86edac 100644 --- a/Source/Filter/NCIGodfreyFilter.H +++ b/Source/Filter/NCIGodfreyFilter.H @@ -13,7 +13,7 @@ public: NCIGodfreyFilter(godfrey_coeff_set coeff_set_, amrex::Real cdtodz_, amrex::Real l_lower_order_in_v_); - virtual void ComputeStencils() override final; + void ComputeStencils(); void getGodfreyCoeffs(godfrey_coeff_set coeff_set_in); diff --git a/Source/Filter/NCIGodfreyFilter.cpp b/Source/Filter/NCIGodfreyFilter.cpp index 28073725a..3f589260a 100644 --- a/Source/Filter/NCIGodfreyFilter.cpp +++ b/Source/Filter/NCIGodfreyFilter.cpp @@ -33,20 +33,6 @@ void NCIGodfreyFilter::ComputeStencils(){ "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"); - - // Godfrey filter uses a table to compute the filter coefficients. - // Get the appropriate table (one for Ex, Ey and Bz, the other - // one for Bx, By, Ez). - Real table_nci[tab_length][tab_width]; - if (coeff_set == godfrey_coeff_set::Ex_Ey_Bz) { - std::copy(&table_nci_godfrey_Ex_Ey_Bz[0][0], - &table_nci_godfrey_Ex_Ey_Bz[0][0]+tab_width*tab_length, - &table_nci[0][0]); - } else if (coeff_set == godfrey_coeff_set::Bx_By_Ez) { - std::copy(&table_nci_godfrey_Bx_By_Ez[0][0], - &table_nci_godfrey_Bx_By_Ez[0][0]+tab_width*tab_length, - &table_nci[0][0]); - } // Interpolate coefficients from the table, and store into prestencil. int index = tab_length*cdtodz; @@ -55,7 +41,13 @@ void NCIGodfreyFilter::ComputeStencils(){ Real weight_right = cdtodz - index/tab_length; Real prestencil[4]; for(int i=0; i Date: Thu, 30 May 2019 16:23:25 -0400 Subject: bug fixes for some refactored GPU filter code --- Source/Filter/Filter.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) (limited to 'Source/Filter/Filter.cpp') diff --git a/Source/Filter/Filter.cpp b/Source/Filter/Filter.cpp index 0c1a98443..f8a2dd6c2 100644 --- a/Source/Filter/Filter.cpp +++ b/Source/Filter/Filter.cpp @@ -64,8 +64,8 @@ Filter::ApplyStencil (FArrayBox& dstfab, const FArrayBox& srcfab, { BL_PROFILE("BilinearFilter::ApplyStencil(FArrayBox)"); ncomp = std::min(ncomp, srcfab.nComp()); - const auto& src = srcfab.array(mfi); - const auto& dst = dstfab.array(mfi); + const auto& src = srcfab.array(); + const auto& dst = dstfab.array(); const Box& gbx = amrex::grow(tbx,stencil_length_each_dir-1); // tmpfab has enough ghost cells for the stencil @@ -74,7 +74,7 @@ Filter::ApplyStencil (FArrayBox& dstfab, const FArrayBox& srcfab, auto const& tmp = tmp_fab.array(); // Copy values in srcfab into tmpfab - const Box& ibx = gbx & srcmf[mfi].box(); + const Box& ibx = gbx & srcfab.box(); AMREX_PARALLEL_FOR_4D ( gbx, ncomp, i, j, k, n, { if (ibx.contains(IntVect(AMREX_D_DECL(i,j,k)))) { -- cgit v1.2.3