From cae3e9ce6769fd5c626e27fad473da17a2b027fb Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Tue, 7 May 2019 17:47:56 -0700 Subject: NCI filter initialized. See how to apply it --- Source/Particles/PhysicalParticleContainer.cpp | 40 ++++++++++++++++++++++++++ 1 file changed, 40 insertions(+) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 17e6d98d9..4a60fc630 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1104,6 +1104,26 @@ PhysicalParticleContainer::Evolve (int lev, const auto& mypc = WarpX::GetInstance().GetPartContainer(); const int nstencilz_fdtd_nci_corr = mypc.nstencilz_fdtd_nci_corr; + + + + + + + + + const auto& nci_godfrey_filter_exeybz = WarpX::GetInstance().nci_godfrey_filter_exeybz; + const auto& nci_godfrey_filter_bxbyez = WarpX::GetInstance().nci_godfrey_filter_bxbyez; + + + + + + + + + + BL_ASSERT(OnSameGrids(lev,jx)); MultiFab* cost = WarpX::getCosts(lev); @@ -1171,53 +1191,73 @@ PhysicalParticleContainer::Evolve (int lev, #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)); + nci_godfrey_filter_exeybz[lev]->ApplyStencil(Ex[pti], filtered_Ex) + /* 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); + */ exfab = &filtered_Ex; filtered_Ez.resize(amrex::convert(tbox,WarpX::Ez_nodal_flag)); + nci_godfrey_filter_bxbyez[lev]->ApplyStencil(Ez[pti], filtered_Ez) + /* 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); + */ ezfab = &filtered_Ez; filtered_By.resize(amrex::convert(tbox,WarpX::By_nodal_flag)); + nci_godfrey_filter_bxbyez[lev]->ApplyStencil(By[pti], filtered_By) + /* 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); + */ byfab = &filtered_By; #if (AMREX_SPACEDIM == 3) filtered_Ey.resize(amrex::convert(tbox,WarpX::Ey_nodal_flag)); + nci_godfrey_filter_exeybz[lev]->ApplyStencil(Ey[pti], filtered_Ey) + /* 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); + */ eyfab = &filtered_Ey; filtered_Bx.resize(amrex::convert(tbox,WarpX::Bx_nodal_flag)); + nci_godfrey_filter_bxbyez[lev]->ApplyStencil(Bx[pti], filtered_Bx) + /* 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); + */ bxfab = &filtered_Bx; filtered_Bz.resize(amrex::convert(tbox,WarpX::Bz_nodal_flag)); + nci_godfrey_filter_exeybz[lev]->ApplyStencil(Bz[pti], filtered_Bz) + /* 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); + */ bzfab = &filtered_Bz; #endif } -- 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/Particles/PhysicalParticleContainer.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/Particles/PhysicalParticleContainer.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/Particles/PhysicalParticleContainer.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 6bca0956123895d7641bf710d69aff04419470af Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Wed, 8 May 2019 17:28:47 -0700 Subject: further cleaning and fix bug when using MR --- Source/Filter/NCIGodfreyFilter.cpp | 1 - Source/Initialization/WarpXInitData.cpp | 2 -- Source/Particles/MultiParticleContainer.H | 11 ---------- Source/Particles/MultiParticleContainer.cpp | 2 -- Source/Particles/PhysicalParticleContainer.cpp | 29 ++++++++++++-------------- Source/WarpX.cpp | 2 +- 6 files changed, 14 insertions(+), 33 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Filter/NCIGodfreyFilter.cpp b/Source/Filter/NCIGodfreyFilter.cpp index 34fca7604..28073725a 100644 --- a/Source/Filter/NCIGodfreyFilter.cpp +++ b/Source/Filter/NCIGodfreyFilter.cpp @@ -55,7 +55,6 @@ void NCIGodfreyFilter::ComputeStencils(){ Real weight_right = cdtodz - index/tab_length; Real prestencil[4]; for(int i=0; ifdtd_nci_stencilz_ex.resize(max_level+1); - mypc->fdtd_nci_stencilz_by.resize(max_level+1); for (int lev = 0; lev <= max_level; ++lev) { const Geometry& gm = Geom(lev); diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H index 0c5e49c04..cc9dc1f59 100644 --- a/Source/Particles/MultiParticleContainer.H +++ b/Source/Particles/MultiParticleContainer.H @@ -178,17 +178,6 @@ public: void UpdateContinuousInjectionPosition(amrex::Real dt) const; int doContinuousInjection() const; - // - // Parameters for the Cherenkov corrector in the FDTD solver. - // Both stencils are calculated ar runtime. - // - // Number of coefficients for the stencil of the NCI corrector. - // The stencil is applied in the z direction only. - static constexpr int nstencilz_fdtd_nci_corr=5; - - amrex::Vector > fdtd_nci_stencilz_ex; - amrex::Vector > fdtd_nci_stencilz_by; - std::vector GetSpeciesNames() const { return species_names; } PhysicalParticleContainer& GetPCtmp () { return *pc_tmp; } diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 440906348..983530569 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -8,8 +8,6 @@ using namespace amrex; -constexpr int MultiParticleContainer::nstencilz_fdtd_nci_corr; - MultiParticleContainer::MultiParticleContainer (AmrCore* amr_core) { ReadParameters(); diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 1e0d68800..d98957c28 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1102,9 +1102,6 @@ PhysicalParticleContainer::Evolve (int lev, const std::array& dx = WarpX::CellSize(lev); const std::array& cdx = WarpX::CellSize(std::max(lev-1,0)); - 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; @@ -1181,38 +1178,38 @@ PhysicalParticleContainer::Evolve (int lev, // 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); + nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Ex, Ex[pti], filtered_Ex.box()); // 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); + nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_Ez, Ez[pti], filtered_Ez.box()); 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); + nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_By, By[pti], filtered_By.box()); 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); + nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Ey, Ey[pti], filtered_Ey.box()); 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); + nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_Bx, Bx[pti], filtered_Bx.box()); 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); + nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Bz, Bz[pti], filtered_Bz.box()); bzfab = &filtered_Bz; #endif } @@ -1388,44 +1385,44 @@ PhysicalParticleContainer::Evolve (int lev, static_cast(WarpX::noy), 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-1]->ApplyStencil(filtered_Ex, (*cEx)[pti], tbox); + nci_godfrey_filter_exeybz[lev-1]->ApplyStencil(filtered_Ex, (*cEx)[pti], filtered_Ex.box()); // 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); + nci_godfrey_filter_bxbyez[lev-1]->ApplyStencil(filtered_Ez, (*cEz)[pti], filtered_Ez.box()); 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); + nci_godfrey_filter_bxbyez[lev-1]->ApplyStencil(filtered_By, (*cBy)[pti], filtered_By.box()); 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); + nci_godfrey_filter_exeybz[lev-1]->ApplyStencil(filtered_Ey, (*cEy)[pti], filtered_Ey.box()); 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); + nci_godfrey_filter_bxbyez[lev-1]->ApplyStencil(filtered_Bx, (*cBx)[pti], filtered_Bx.box()); 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); + nci_godfrey_filter_exeybz[lev-1]->ApplyStencil(filtered_Bz, (*cBz)[pti], filtered_Bz.box()); cbzfab = &filtered_Bz; #endif } diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 5f63e2ad6..b24058a0e 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -651,7 +651,7 @@ WarpX::AllocLevelData (int lev, const BoxArray& ba, const DistributionMapping& d int ngz_nonci = (ngz_tmp % 2) ? ngz_tmp+1 : ngz_tmp; // Always even number int ngz; if (WarpX::use_fdtd_nci_corr) { - int ng = ngz_tmp + (mypc->nstencilz_fdtd_nci_corr-1); + int ng = ngz_tmp + 4; ngz = (ng % 2) ? ng+1 : ng; } else { ngz = ngz_nonci; -- cgit v1.2.3 From 97e66052cc96ed40bef74ed3d62bed3fdb26427a Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Wed, 8 May 2019 18:11:33 -0700 Subject: remove redundant Elixir declaration --- Source/Particles/PhysicalParticleContainer.cpp | 1 - 1 file changed, 1 deletion(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index d98957c28..680086786 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1374,7 +1374,6 @@ 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) -- cgit v1.2.3 From 5308cbefa8fbde30bed8a88943d58ec646731d9e Mon Sep 17 00:00:00 2001 From: Andrew Myers Date: Thu, 23 May 2019 15:27:17 -0700 Subject: fix a bunch of unused variable / parameter shadowing warnings --- Source/Diagnostics/BoostedFrameDiagnostic.H | 2 +- Source/Diagnostics/BoostedFrameDiagnostic.cpp | 5 --- Source/Diagnostics/FieldIO.cpp | 1 - Source/Evolve/WarpXEvolveEM.cpp | 28 +++++++------- Source/FieldSolver/WarpXPushFieldsEM.cpp | 48 ++++++++++++------------ Source/FortranInterface/WarpX_picsar.F90 | 7 +++- Source/Laser/LaserParticleContainer.cpp | 6 +-- Source/Parallelization/WarpXComm.cpp | 52 +++++++++++++------------- Source/Particles/PhysicalParticleContainer.cpp | 3 -- Source/Particles/WarpXParticleContainer.cpp | 44 +++++++++++----------- Source/Utils/WarpXMovingWindow.cpp | 6 +-- Source/WarpX.cpp | 1 - 12 files changed, 96 insertions(+), 107 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Diagnostics/BoostedFrameDiagnostic.H b/Source/Diagnostics/BoostedFrameDiagnostic.H index 55a124c51..3a09259b0 100644 --- a/Source/Diagnostics/BoostedFrameDiagnostic.H +++ b/Source/Diagnostics/BoostedFrameDiagnostic.H @@ -40,8 +40,8 @@ class BoostedFrameDiagnostic { amrex::Real current_z_lab; amrex::Real current_z_boost; - std::vector mesh_field_names_; int ncomp_to_dump_; + std::vector mesh_field_names_; int file_num; int initial_i; diff --git a/Source/Diagnostics/BoostedFrameDiagnostic.cpp b/Source/Diagnostics/BoostedFrameDiagnostic.cpp index f62369091..57872790f 100644 --- a/Source/Diagnostics/BoostedFrameDiagnostic.cpp +++ b/Source/Diagnostics/BoostedFrameDiagnostic.cpp @@ -497,8 +497,6 @@ LorentzTransformZ(MultiFab& data, Real gamma_boost, Real beta_boost, int ncomp) // Transform the transverse E and B fields. Note that ez and bz are not // changed by the tranform. Real e_lab, b_lab, j_lab, r_lab; - int i0 = 0; - int i4 = 4; e_lab = gamma_boost * (arr(i, j, k, 0) + beta_boost*clight*arr(i, j, k, 4)); b_lab = gamma_boost * (arr(i, j, k, 4) + beta_boost*arr(i, j, k, 0)/clight); @@ -718,7 +716,6 @@ writeLabFrameData(const MultiFab* cell_centered_data, if (buff_counter_[i] == 0) { // ... reset fields buffer data_buffer_[i] if (WarpX::do_boosted_frame_fields) { - const int ncomp = cell_centered_data->nComp(); Box buff_box = geom.Domain(); buff_box.setSmall(boost_direction_, i_lab - num_buffer_ + 1); buff_box.setBig(boost_direction_, i_lab); @@ -726,7 +723,6 @@ writeLabFrameData(const MultiFab* cell_centered_data, buff_ba.maxSize(max_box_size_); DistributionMapping buff_dm(buff_ba); data_buffer_[i].reset( new MultiFab(buff_ba, buff_dm, ncomp_to_dump, 0) ); - // data_buffer_[i].reset( new MultiFab(buff_ba, buff_dm, ncomp, 0) ); } // ... reset particle buffer particles_buffer_[i] if (WarpX::do_boosted_frame_particles) @@ -961,7 +957,6 @@ LabSnapShot(Real t_lab_in, Real t_boost, RealBox prob_domain_lab, my_bfd(bfd) { Real zmin_lab = prob_domain_lab_.lo(AMREX_SPACEDIM-1); - Real zmax_lab = prob_domain_lab_.hi(AMREX_SPACEDIM-1); current_z_lab = 0.0; current_z_boost = 0.0; updateCurrentZPositions(t_boost, my_bfd.inv_gamma_boost_, my_bfd.inv_beta_boost_); diff --git a/Source/Diagnostics/FieldIO.cpp b/Source/Diagnostics/FieldIO.cpp index ced1f8ea3..e3d44d1fc 100644 --- a/Source/Diagnostics/FieldIO.cpp +++ b/Source/Diagnostics/FieldIO.cpp @@ -660,7 +660,6 @@ getInterpolatedScalar( interpolated_F->setVal(0.); // Loop through the boxes and interpolate the values from the _cp data - const int use_limiter = 0; #ifdef _OPEMP #pragma omp parallel #endif diff --git a/Source/Evolve/WarpXEvolveEM.cpp b/Source/Evolve/WarpXEvolveEM.cpp index dab58f95b..ad7c7d840 100644 --- a/Source/Evolve/WarpXEvolveEM.cpp +++ b/Source/Evolve/WarpXEvolveEM.cpp @@ -514,7 +514,7 @@ WarpX::ComputeDt () * simulation box passes input parameter zmax_plasma_to_compute_max_step. */ void -WarpX::computeMaxStepBoostAccelerator(amrex::Geometry geom){ +WarpX::computeMaxStepBoostAccelerator(amrex::Geometry a_geom){ // Sanity checks: can use zmax_plasma_to_compute_max_step only if // the moving window and the boost are all in z direction. AMREX_ALWAYS_ASSERT_WITH_MESSAGE( @@ -534,7 +534,7 @@ WarpX::computeMaxStepBoostAccelerator(amrex::Geometry geom){ // Lower end of the simulation domain. All quantities are given in boosted // frame except zmax_plasma_to_compute_max_step. - const Real zmin_domain_boost = geom.ProbLo(AMREX_SPACEDIM-1); + const Real zmin_domain_boost = a_geom.ProbLo(AMREX_SPACEDIM-1); // End of the plasma: Transform input argument // zmax_plasma_to_compute_max_step to boosted frame. const Real len_plasma_boost = zmax_plasma_to_compute_max_step/gamma_boost; @@ -591,19 +591,19 @@ WarpX::applyMirrors(Real time){ NullifyMF(Bz, lev, z_min, z_max); if (lev>0){ // Get coarse patch field MultiFabs - MultiFab& Ex = *Efield_cp[lev][0].get(); - MultiFab& Ey = *Efield_cp[lev][1].get(); - MultiFab& Ez = *Efield_cp[lev][2].get(); - MultiFab& Bx = *Bfield_cp[lev][0].get(); - MultiFab& By = *Bfield_cp[lev][1].get(); - MultiFab& Bz = *Bfield_cp[lev][2].get(); + MultiFab& cEx = *Efield_cp[lev][0].get(); + MultiFab& cEy = *Efield_cp[lev][1].get(); + MultiFab& cEz = *Efield_cp[lev][2].get(); + MultiFab& cBx = *Bfield_cp[lev][0].get(); + MultiFab& cBy = *Bfield_cp[lev][1].get(); + MultiFab& cBz = *Bfield_cp[lev][2].get(); // Set each field to zero between z_min and z_max - NullifyMF(Ex, lev, z_min, z_max); - NullifyMF(Ey, lev, z_min, z_max); - NullifyMF(Ez, lev, z_min, z_max); - NullifyMF(Bx, lev, z_min, z_max); - NullifyMF(By, lev, z_min, z_max); - NullifyMF(Bz, lev, z_min, z_max); + NullifyMF(cEx, lev, z_min, z_max); + NullifyMF(cEy, lev, z_min, z_max); + NullifyMF(cEz, lev, z_min, z_max); + NullifyMF(cBx, lev, z_min, z_max); + NullifyMF(cBy, lev, z_min, z_max); + NullifyMF(cBz, lev, z_min, z_max); } } } diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index a66d14980..c53e13f8f 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -17,30 +17,30 @@ using namespace amrex; void -WarpX::EvolveB (Real dt) +WarpX::EvolveB (Real a_dt) { for (int lev = 0; lev <= finest_level; ++lev) { - EvolveB(lev, dt); + EvolveB(lev, a_dt); } } void -WarpX::EvolveB (int lev, Real dt) +WarpX::EvolveB (int lev, Real a_dt) { BL_PROFILE("WarpX::EvolveB()"); - EvolveB(lev, PatchType::fine, dt); + EvolveB(lev, PatchType::fine, a_dt); if (lev > 0) { - EvolveB(lev, PatchType::coarse, dt); + EvolveB(lev, PatchType::coarse, a_dt); } } void -WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real dt) +WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real a_dt) { const int patch_level = (patch_type == PatchType::fine) ? lev : lev-1; const std::array& dx = WarpX::CellSize(patch_level); - Real dtsdx = dt/dx[0], dtsdy = dt/dx[1], dtsdz = dt/dx[2]; + Real dtsdx = a_dt/dx[0], dtsdy = a_dt/dx[1], dtsdz = a_dt/dx[2]; MultiFab *Ex, *Ey, *Ez, *Bx, *By, *Bz; if (patch_type == PatchType::fine) @@ -164,30 +164,30 @@ WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real dt) } void -WarpX::EvolveE (Real dt) +WarpX::EvolveE (Real a_dt) { for (int lev = 0; lev <= finest_level; ++lev) { - EvolveE(lev, dt); + EvolveE(lev, a_dt); } } void -WarpX::EvolveE (int lev, Real dt) +WarpX::EvolveE (int lev, Real a_dt) { BL_PROFILE("WarpX::EvolveE()"); - EvolveE(lev, PatchType::fine, dt); + EvolveE(lev, PatchType::fine, a_dt); if (lev > 0) { - EvolveE(lev, PatchType::coarse, dt); + EvolveE(lev, PatchType::coarse, a_dt); } } void -WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real dt) +WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real a_dt) { - const Real mu_c2_dt = (PhysConst::mu0*PhysConst::c*PhysConst::c) * dt; - const Real c2dt = (PhysConst::c*PhysConst::c) * dt; + const Real mu_c2_dt = (PhysConst::mu0*PhysConst::c*PhysConst::c) * a_dt; + const Real c2dt = (PhysConst::c*PhysConst::c) * a_dt; int patch_level = (patch_type == PatchType::fine) ? lev : lev-1; const std::array& dx = WarpX::CellSize(patch_level); @@ -358,27 +358,27 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real dt) } void -WarpX::EvolveF (Real dt, DtType dt_type) +WarpX::EvolveF (Real a_dt, DtType a_dt_type) { if (!do_dive_cleaning) return; for (int lev = 0; lev <= finest_level; ++lev) { - EvolveF(lev, dt, dt_type); + EvolveF(lev, a_dt, a_dt_type); } } void -WarpX::EvolveF (int lev, Real dt, DtType dt_type) +WarpX::EvolveF (int lev, Real a_dt, DtType a_dt_type) { if (!do_dive_cleaning) return; - EvolveF(lev, PatchType::fine, dt, dt_type); - if (lev > 0) EvolveF(lev, PatchType::coarse, dt, dt_type); + EvolveF(lev, PatchType::fine, a_dt, a_dt_type); + if (lev > 0) EvolveF(lev, PatchType::coarse, a_dt, a_dt_type); } void -WarpX::EvolveF (int lev, PatchType patch_type, Real dt, DtType dt_type) +WarpX::EvolveF (int lev, PatchType patch_type, Real a_dt, DtType a_dt_type) { if (!do_dive_cleaning) return; @@ -388,7 +388,7 @@ WarpX::EvolveF (int lev, PatchType patch_type, Real dt, DtType dt_type) int patch_level = (patch_type == PatchType::fine) ? lev : lev-1; const auto& dx = WarpX::CellSize(patch_level); - const std::array dtsdx {dt/dx[0], dt/dx[1], dt/dx[2]}; + const std::array dtsdx {a_dt/dx[0], a_dt/dx[1], a_dt/dx[2]}; MultiFab *Ex, *Ey, *Ez, *rho, *F; if (patch_type == PatchType::fine) @@ -408,12 +408,12 @@ WarpX::EvolveF (int lev, PatchType patch_type, Real dt, DtType dt_type) F = F_cp[lev].get(); } - const int rhocomp = (dt_type == DtType::FirstHalf) ? 0 : 1; + const int rhocomp = (a_dt_type == DtType::FirstHalf) ? 0 : 1; MultiFab src(rho->boxArray(), rho->DistributionMap(), 1, 0); ComputeDivE(src, 0, {Ex,Ey,Ez}, dx); MultiFab::Saxpy(src, -mu_c2, *rho, rhocomp, 0, 1, 0); - MultiFab::Saxpy(*F, dt, src, 0, 0, 1, 0); + MultiFab::Saxpy(*F, a_dt, src, 0, 0, 1, 0); if (do_pml && pml[lev]->ok()) { diff --git a/Source/FortranInterface/WarpX_picsar.F90 b/Source/FortranInterface/WarpX_picsar.F90 index a4cc99926..c17e8861b 100644 --- a/Source/FortranInterface/WarpX_picsar.F90 +++ b/Source/FortranInterface/WarpX_picsar.F90 @@ -265,8 +265,10 @@ subroutine warpx_charge_deposition(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,n real(amrex_real), intent(IN OUT):: rho(*) real(amrex_real), intent(IN) :: rmin, dr +#ifdef WARPX_RZ integer(c_long) :: type_rz_depose = 1 - +#endif + ! Compute the number of valid cells and guard cells integer(c_long) :: rho_nvalid(AMREX_SPACEDIM), rho_nguards(AMREX_SPACEDIM) rho_nvalid = rho_ntot - 2*rho_ng @@ -385,8 +387,9 @@ subroutine warpx_charge_deposition(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,n real(amrex_real), intent(IN OUT):: jx(*), jy(*), jz(*) real(amrex_real), intent(IN) :: rmin, dr +#ifdef WARPX_RZ integer(c_long) :: type_rz_depose = 1 - +#endif ! Compute the number of valid cells and guard cells integer(c_long) :: jx_nvalid(AMREX_SPACEDIM), jy_nvalid(AMREX_SPACEDIM), jz_nvalid(AMREX_SPACEDIM), & jx_nguards(AMREX_SPACEDIM), jy_nguards(AMREX_SPACEDIM), jz_nguards(AMREX_SPACEDIM) diff --git a/Source/Laser/LaserParticleContainer.cpp b/Source/Laser/LaserParticleContainer.cpp index 52d506bb8..2f964b6f3 100644 --- a/Source/Laser/LaserParticleContainer.cpp +++ b/Source/Laser/LaserParticleContainer.cpp @@ -78,14 +78,14 @@ LaserParticleContainer::LaserParticleContainer (AmrCore* amr_core, int ispecies, parser.define(field_function); parser.registerVariables({"X","Y","t"}); - ParmParse pp("my_constants"); + ParmParse ppc("my_constants"); std::set symbols = parser.symbols(); symbols.erase("X"); symbols.erase("Y"); symbols.erase("t"); // after removing variables, we are left with constants for (auto it = symbols.begin(); it != symbols.end(); ) { Real v; - if (pp.query(it->c_str(), v)) { + if (ppc.query(it->c_str(), v)) { parser.setConstant(*it, v); it = symbols.erase(it); } else { @@ -429,8 +429,6 @@ LaserParticleContainer::Evolve (int lev, { Real wt = amrex::second(); - const Box& box = pti.validbox(); - auto& attribs = pti.GetAttribs(); auto& wp = attribs[PIdx::w ]; diff --git a/Source/Parallelization/WarpXComm.cpp b/Source/Parallelization/WarpXComm.cpp index 348b31c8b..1cf347420 100644 --- a/Source/Parallelization/WarpXComm.cpp +++ b/Source/Parallelization/WarpXComm.cpp @@ -79,7 +79,7 @@ WarpX::UpdateAuxilaryData () MultiFab::Subtract(dBz, *Bfield_cp[lev][2], 0, 0, 1, ng); const Real* dx = Geom(lev-1).CellSize(); - const int ref_ratio = refRatio(lev-1)[0]; + const int rr = refRatio(lev-1)[0]; #ifdef _OPENMP #pragma omp parallel #endif @@ -89,7 +89,7 @@ WarpX::UpdateAuxilaryData () { Box ccbx = mfi.fabbox(); ccbx.enclosedCells(); - ccbx.coarsen(ref_ratio).refine(ref_ratio); // so that ccbx is coarsenable + ccbx.coarsen(rr).refine(rr); // so that ccbx is coarsenable const FArrayBox& cxfab = dBx[mfi]; const FArrayBox& cyfab = dBy[mfi]; @@ -106,18 +106,18 @@ WarpX::UpdateAuxilaryData () BL_TO_FORTRAN_ANYD(cxfab), BL_TO_FORTRAN_ANYD(cyfab), BL_TO_FORTRAN_ANYD(czfab), - dx, &ref_ratio,&use_limiter); + dx, &rr,&use_limiter); #else amrex_interp_div_free_bfield(ccbx.loVect(), ccbx.hiVect(), BL_TO_FORTRAN_ANYD(bfab[0]), BL_TO_FORTRAN_ANYD(bfab[2]), BL_TO_FORTRAN_ANYD(cxfab), BL_TO_FORTRAN_ANYD(czfab), - dx, &ref_ratio,&use_limiter); + dx, &rr,&use_limiter); amrex_interp_cc_bfield(ccbx.loVect(), ccbx.hiVect(), BL_TO_FORTRAN_ANYD(bfab[1]), BL_TO_FORTRAN_ANYD(cyfab), - &ref_ratio,&use_limiter); + &rr,&use_limiter); #endif for (int idim = 0; idim < 3; ++idim) @@ -153,7 +153,7 @@ WarpX::UpdateAuxilaryData () MultiFab::Subtract(dEy, *Efield_cp[lev][1], 0, 0, 1, ng); MultiFab::Subtract(dEz, *Efield_cp[lev][2], 0, 0, 1, ng); - const int ref_ratio = refRatio(lev-1)[0]; + const int rr = refRatio(lev-1)[0]; #ifdef _OPEMP #pragma omp parallel #endif @@ -163,7 +163,7 @@ WarpX::UpdateAuxilaryData () { Box ccbx = mfi.fabbox(); ccbx.enclosedCells(); - ccbx.coarsen(ref_ratio).refine(ref_ratio); // so that ccbx is coarsenable + ccbx.coarsen(rr).refine(rr); // so that ccbx is coarsenable const FArrayBox& cxfab = dEx[mfi]; const FArrayBox& cyfab = dEy[mfi]; @@ -180,18 +180,18 @@ WarpX::UpdateAuxilaryData () BL_TO_FORTRAN_ANYD(cxfab), BL_TO_FORTRAN_ANYD(cyfab), BL_TO_FORTRAN_ANYD(czfab), - &ref_ratio,&use_limiter); + &rr,&use_limiter); #else amrex_interp_efield(ccbx.loVect(), ccbx.hiVect(), BL_TO_FORTRAN_ANYD(efab[0]), BL_TO_FORTRAN_ANYD(efab[2]), BL_TO_FORTRAN_ANYD(cxfab), BL_TO_FORTRAN_ANYD(czfab), - &ref_ratio,&use_limiter); + &rr,&use_limiter); amrex_interp_nd_efield(ccbx.loVect(), ccbx.hiVect(), BL_TO_FORTRAN_ANYD(efab[1]), BL_TO_FORTRAN_ANYD(cyfab), - &ref_ratio); + &rr); #endif for (int idim = 0; idim < 3; ++idim) @@ -364,7 +364,7 @@ WarpX::SyncCurrent () current_cp[lev][1]->setVal(0.0); current_cp[lev][2]->setVal(0.0); - const IntVect& ref_ratio = refRatio(lev-1); + const IntVect& rr = refRatio(lev-1); std::array fine { current_fp[lev][0].get(), current_fp[lev][1].get(), @@ -372,7 +372,7 @@ WarpX::SyncCurrent () std::array< MultiFab*,3> crse { current_cp[lev][0].get(), current_cp[lev][1].get(), current_cp[lev][2].get() }; - SyncCurrent(fine, crse, ref_ratio[0]); + SyncCurrent(fine, crse, rr[0]); } Vector,3> > j_fp(finest_level+1); @@ -524,10 +524,10 @@ WarpX::SyncCurrent () void WarpX::SyncCurrent (const std::array& fine, const std::array< amrex::MultiFab*,3>& crse, - int ref_ratio) + int rr) { - BL_ASSERT(ref_ratio == 2); - const IntVect& ng = (fine[0]->nGrowVect() + 1) /ref_ratio; + BL_ASSERT(rr == 2); + const IntVect& ng = (fine[0]->nGrowVect() + 1) /rr; #ifdef _OPEMP #pragma omp parallel @@ -539,7 +539,7 @@ WarpX::SyncCurrent (const std::array& fine, for (MFIter mfi(*crse[idim],true); mfi.isValid(); ++mfi) { const Box& bx = mfi.growntilebox(ng); - Box fbx = amrex::grow(amrex::refine(bx,ref_ratio),1); + Box fbx = amrex::grow(amrex::refine(bx,rr),1); ffab.resize(fbx); fbx &= (*fine[idim])[mfi].box(); ffab.setVal(0.0); @@ -564,8 +564,8 @@ WarpX::SyncRho (amrex::Vector >& rhof, for (int lev = 1; lev <= finest_level; ++lev) { rhoc[lev]->setVal(0.0); - const IntVect& ref_ratio = refRatio(lev-1); - SyncRho(*rhof[lev], *rhoc[lev], ref_ratio[0]); + const IntVect& rr = refRatio(lev-1); + SyncRho(*rhof[lev], *rhoc[lev], rr[0]); } Vector > rho_f_g(finest_level+1); @@ -673,10 +673,10 @@ WarpX::SyncRho (amrex::Vector >& rhof, * averaging the values of the charge density of the fine patch (on the same level). */ void -WarpX::SyncRho (const MultiFab& fine, MultiFab& crse, int ref_ratio) +WarpX::SyncRho (const MultiFab& fine, MultiFab& crse, int rr) { - BL_ASSERT(ref_ratio == 2); - const IntVect& ng = (fine.nGrowVect()+1)/ref_ratio; + BL_ASSERT(rr == 2); + const IntVect& ng = (fine.nGrowVect()+1)/rr; const int nc = fine.nComp(); #ifdef _OPEMP @@ -687,7 +687,7 @@ WarpX::SyncRho (const MultiFab& fine, MultiFab& crse, int ref_ratio) for (MFIter mfi(crse,true); mfi.isValid(); ++mfi) { const Box& bx = mfi.growntilebox(ng); - Box fbx = amrex::grow(amrex::refine(bx,ref_ratio),1); + Box fbx = amrex::grow(amrex::refine(bx,rr),1); ffab.resize(fbx, nc); fbx &= fine[mfi].box(); ffab.setVal(0.0); @@ -710,7 +710,7 @@ WarpX::RestrictCurrentFromFineToCoarsePatch (int lev) current_cp[lev][1]->setVal(0.0); current_cp[lev][2]->setVal(0.0); - const IntVect& ref_ratio = refRatio(lev-1); + const IntVect& rr = refRatio(lev-1); std::array fine { current_fp[lev][0].get(), current_fp[lev][1].get(), @@ -718,7 +718,7 @@ WarpX::RestrictCurrentFromFineToCoarsePatch (int lev) std::array< MultiFab*,3> crse { current_cp[lev][0].get(), current_cp[lev][1].get(), current_cp[lev][2].get() }; - SyncCurrent(fine, crse, ref_ratio[0]); + SyncCurrent(fine, crse, rr[0]); } void @@ -824,8 +824,8 @@ WarpX::RestrictRhoFromFineToCoarsePatch (int lev) { if (rho_fp[lev]) { rho_cp[lev]->setVal(0.0); - const IntVect& ref_ratio = refRatio(lev-1); - SyncRho(*rho_fp[lev], *rho_cp[lev], ref_ratio[0]); + const IntVect& rr = refRatio(lev-1); + SyncRho(*rho_fp[lev], *rho_cp[lev], rr[0]); } } diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 37c136a3d..212084e64 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1024,9 +1024,6 @@ PhysicalParticleContainer::FieldGather (int lev, { const std::array& dx = WarpX::CellSize(lev); - // WarpX assumes the same number of guard cells for Ex, Ey, Ez, Bx, By, Bz - long ng = Ex.nGrow(); - BL_ASSERT(OnSameGrids(lev,Ex)); MultiFab* cost = WarpX::getCosts(lev); diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index c7ffe956b..9791eee80 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -237,10 +237,10 @@ WarpXParticleContainer::AddNParticles (int lev, if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { - auto& particle_tile = DefineAndReturnParticleTile(0, 0, 0); - particle_tile.push_back_real(particle_comps["xold"], x[i]); - particle_tile.push_back_real(particle_comps["yold"], y[i]); - particle_tile.push_back_real(particle_comps["zold"], z[i]); + auto& ptile = DefineAndReturnParticleTile(0, 0, 0); + ptile.push_back_real(particle_comps["xold"], x[i]); + ptile.push_back_real(particle_comps["yold"], y[i]); + ptile.push_back_real(particle_comps["zold"], z[i]); } particle_tile.push_back(p); @@ -255,10 +255,10 @@ WarpXParticleContainer::AddNParticles (int lev, if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { - auto& particle_tile = DefineAndReturnParticleTile(0, 0, 0); - particle_tile.push_back_real(particle_comps["uxold"], vx + ibegin, vx + iend); - particle_tile.push_back_real(particle_comps["uyold"], vy + ibegin, vy + iend); - particle_tile.push_back_real(particle_comps["uzold"], vz + ibegin, vz + iend); + auto& ptile = DefineAndReturnParticleTile(0, 0, 0); + ptile.push_back_real(particle_comps["uxold"], vx + ibegin, vx + iend); + ptile.push_back_real(particle_comps["uyold"], vy + ibegin, vy + iend); + ptile.push_back_real(particle_comps["uzold"], vz + ibegin, vz + iend); } for (int comp = PIdx::uz+1; comp < PIdx::nattribs; ++comp) @@ -737,7 +737,6 @@ WarpXParticleContainer::DepositCharge (Vector >& rho, const auto& gm = m_gdb->Geom(lev); const auto& ba = m_gdb->ParticleBoxArray(lev); - const auto& dm = m_gdb->DistributionMap(lev); const Real* dx = gm.CellSize(); const Real* plo = gm.ProbLo(); @@ -807,36 +806,36 @@ WarpXParticleContainer::GetChargeDensity (int lev, bool local) #ifdef _OPENMP #pragma omp parallel -#endif { +#endif Cuda::ManagedDeviceVector xp, yp, zp; - FArrayBox local_rho; +#ifdef _OPENMP + FArrayBox rho_loc; +#endif for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) { - const Box& box = pti.validbox(); - auto& wp = pti.GetAttribs(PIdx::w); const long np = pti.numParticles(); pti.GetPosition(xp, yp, zp); - const std::array& xyzmin_tile = WarpX::LowerCorner(pti.tilebox(), lev); - const std::array& xyzmin_grid = WarpX::LowerCorner(box, lev); - // Data on the grid Real* data_ptr; FArrayBox& rhofab = (*rho)[pti]; #ifdef _OPENMP + const std::array& xyzmin_tile = WarpX::LowerCorner(pti.tilebox(), lev); Box tile_box = convert(pti.tilebox(), IntVect::TheUnitVector()); const std::array& xyzmin = xyzmin_tile; tile_box.grow(ng); - local_rho.resize(tile_box); - local_rho = 0.0; - data_ptr = local_rho.dataPtr(); - auto rholen = local_rho.length(); + rho_loc.resize(tile_box); + rho_loc = 0.0; + data_ptr = rho_loc.dataPtr(); + auto rholen = rho_loc.length(); #else + const Box& box = pti.validbox(); + const std::array& xyzmin_grid = WarpX::LowerCorner(box, lev); const std::array& xyzmin = xyzmin_grid; data_ptr = rhofab.dataPtr(); auto rholen = rhofab.length(); @@ -874,10 +873,9 @@ WarpXParticleContainer::GetChargeDensity (int lev, bool local) #endif #ifdef _OPENMP - rhofab.atomicAdd(local_rho); -#endif + rhofab.atomicAdd(rho_loc); } - +#endif } if (!local) rho->SumBoundary(gm.periodicity()); diff --git a/Source/Utils/WarpXMovingWindow.cpp b/Source/Utils/WarpXMovingWindow.cpp index 5f9e36bf6..a0ab1f26f 100644 --- a/Source/Utils/WarpXMovingWindow.cpp +++ b/Source/Utils/WarpXMovingWindow.cpp @@ -5,7 +5,7 @@ using namespace amrex; void -WarpX::UpdatePlasmaInjectionPosition (Real dt) +WarpX::UpdatePlasmaInjectionPosition (Real a_dt) { int dir = moving_window_dir; // Continuously inject plasma in new cells (by default only on level 0) @@ -14,12 +14,12 @@ WarpX::UpdatePlasmaInjectionPosition (Real dt) // call to this function, and injection position needs to be updated current_injection_position -= WarpX::beta_boost * #if ( AMREX_SPACEDIM == 3 ) - WarpX::boost_direction[dir] * PhysConst::c * dt; + WarpX::boost_direction[dir] * PhysConst::c * a_dt; #elif ( AMREX_SPACEDIM == 2 ) // In 2D, dir=0 corresponds to x and dir=1 corresponds to z // This needs to be converted in order to index `boost_direction` // which has 3 components, for both 2D and 3D simulations. - WarpX::boost_direction[2*dir] * PhysConst::c * dt; + WarpX::boost_direction[2*dir] * PhysConst::c * a_dt; #endif } } diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 6b6752bf1..3d7f7dcc5 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -1001,7 +1001,6 @@ WarpX::ComputeDivE (MultiFab& divE, int dcomp, void WarpX::BuildBufferMasks () { - int ngbuffer = std::max(n_field_gather_buffer, n_current_deposition_buffer); for (int lev = 1; lev <= maxLevel(); ++lev) { for (int ipass = 0; ipass < 2; ++ipass) -- cgit v1.2.3 From 46b9fde64b935703e9bda18f6150a165c14d4032 Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Fri, 24 May 2019 16:19:14 -0700 Subject: option do symmetrize gaussian particle beam --- Source/Initialization/PlasmaInjector.H | 1 + Source/Initialization/PlasmaInjector.cpp | 1 + Source/Particles/PhysicalParticleContainer.H | 5 +- Source/Particles/PhysicalParticleContainer.cpp | 91 +++++++++++++++++--------- 4 files changed, 66 insertions(+), 32 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Initialization/PlasmaInjector.H b/Source/Initialization/PlasmaInjector.H index cd5802a55..f998e217e 100644 --- a/Source/Initialization/PlasmaInjector.H +++ b/Source/Initialization/PlasmaInjector.H @@ -293,6 +293,7 @@ public: amrex::Real z_rms; amrex::Real q_tot; long npart; + int do_symmetrize = 0; bool radially_weighted = true; diff --git a/Source/Initialization/PlasmaInjector.cpp b/Source/Initialization/PlasmaInjector.cpp index 52b5d8b61..f9642d1b6 100644 --- a/Source/Initialization/PlasmaInjector.cpp +++ b/Source/Initialization/PlasmaInjector.cpp @@ -286,6 +286,7 @@ PlasmaInjector::PlasmaInjector(int ispecies, const std::string& name) pp.get("z_rms", z_rms); pp.get("q_tot", q_tot); pp.get("npart", npart); + pp.query("do_symmetrize", do_symmetrize); gaussian_beam = true; parseMomentum(pp); } diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index df1e0296a..bd8cfb47e 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -94,7 +94,10 @@ public: void AddGaussianBeam(amrex::Real x_m, amrex::Real y_m, amrex::Real z_m, amrex::Real x_rms, amrex::Real y_rms, amrex::Real z_rms, - amrex::Real q_tot, long npart); + amrex::Real q_tot, long npart, int do_symmetrize); + + void CheckAndAddParticle(amrex::Real x, amrex::Real y, amrex::Real z, + std::array u, amrex::Real weight); virtual void GetParticleSlice(const int direction, const amrex::Real z_old, const amrex::Real z_new, const amrex::Real t_boost, diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 212084e64..38ea55d6e 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -181,7 +181,8 @@ void PhysicalParticleContainer::MapParticletoBoostedFrame(Real& x, Real& y, Real void PhysicalParticleContainer::AddGaussianBeam(Real x_m, Real y_m, Real z_m, Real x_rms, Real y_rms, Real z_rms, - Real q_tot, long npart) { + Real q_tot, long npart, + int do_symmetrize) { const Geometry& geom = m_gdb->Geom(0); RealBox containing_bx = geom.ProbDomain(); @@ -191,13 +192,15 @@ PhysicalParticleContainer::AddGaussianBeam(Real x_m, Real y_m, Real z_m, std::normal_distribution disty(y_m, y_rms); std::normal_distribution distz(z_m, z_rms); - std::array attribs; - attribs.fill(0.0); - if (ParallelDescriptor::IOProcessor()) { - std::array u; - Real weight; - for (long i = 0; i < npart; ++i) { + std::array u; + Real weight; + // If do_symmetrize, create 4x fewer particles, and + // Replicate each particle 4 times (x,y) (-x,y) (x,-y) (-x,-y) + if (do_symmetrize){ + npart /= 4; + } + for (long i = 0; i < npart; ++i) { #if ( AMREX_SPACEDIM == 3 | WARPX_RZ) weight = q_tot/npart/charge; Real x = distx(mt); @@ -209,35 +212,60 @@ PhysicalParticleContainer::AddGaussianBeam(Real x_m, Real y_m, Real z_m, Real y = 0.; Real z = distz(mt); #endif - if (plasma_injector->insideBounds(x, y, z)) { - plasma_injector->getMomentum(u, x, y, z); - if (WarpX::gamma_boost > 1.) { - MapParticletoBoostedFrame(x, y, z, u); - } - attribs[PIdx::ux] = u[0]; - attribs[PIdx::uy] = u[1]; - attribs[PIdx::uz] = u[2]; - attribs[PIdx::w ] = weight; - - if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) - { - auto& particle_tile = DefineAndReturnParticleTile(0, 0, 0); - particle_tile.push_back_real(particle_comps["xold"], x); - particle_tile.push_back_real(particle_comps["yold"], y); - particle_tile.push_back_real(particle_comps["zold"], z); - - particle_tile.push_back_real(particle_comps["uxold"], u[0]); - particle_tile.push_back_real(particle_comps["uyold"], u[1]); - particle_tile.push_back_real(particle_comps["uzold"], u[2]); - } - - AddOneParticle(0, 0, 0, x, y, z, attribs); + if (plasma_injector->insideBounds(x, y, z)) { + plasma_injector->getMomentum(u, x, y, z); + if (do_symmetrize){ + // Add four particles to the beam: + // (x,ux,y,uy) (-x,-ux,y,uy) (x,ux,-y,-uy) (-x,-ux,-y,-uy) + std::array u_tmp = u; + for (int ix=0; ix<2; ix++){ + for (int iy=0; iy<2; iy++){ + x *= std::pow(-1,ix); + u_tmp[0] *= std::pow(-1,ix); + y *= std::pow(-1,iy); + u_tmp[1] *= std::pow(-1,iy); + CheckAndAddParticle(x, y, z, u_tmp, weight); + } + } + } else { + CheckAndAddParticle(x, y, z, u, weight); + } } } } Redistribute(); } +void +PhysicalParticleContainer::CheckAndAddParticle(Real x, Real y, Real z, + std::array u, + Real weight) +{ + std::array attribs; + attribs.fill(0.0); + + if (WarpX::gamma_boost > 1.) { + MapParticletoBoostedFrame(x, y, z, u); + } + attribs[PIdx::ux] = u[0]; + attribs[PIdx::uy] = u[1]; + attribs[PIdx::uz] = u[2]; + attribs[PIdx::w ] = weight; + + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) + { + auto& particle_tile = DefineAndReturnParticleTile(0, 0, 0); + particle_tile.push_back_real(particle_comps["xold"], x); + particle_tile.push_back_real(particle_comps["yold"], y); + particle_tile.push_back_real(particle_comps["zold"], z); + + particle_tile.push_back_real(particle_comps["uxold"], u[0]); + particle_tile.push_back_real(particle_comps["uyold"], u[1]); + particle_tile.push_back_real(particle_comps["uzold"], u[2]); + } + AddOneParticle(0, 0, 0, x, y, z, attribs); +} + void PhysicalParticleContainer::AddParticles (int lev) { @@ -263,7 +291,8 @@ PhysicalParticleContainer::AddParticles (int lev) plasma_injector->y_rms, plasma_injector->z_rms, plasma_injector->q_tot, - plasma_injector->npart); + plasma_injector->npart, + plasma_injector->do_symmetrize); return; -- cgit v1.2.3 From 5473027f14805cae5f00f21a8eef2ca988d9f8ab Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Fri, 24 May 2019 16:27:09 -0700 Subject: add doc --- Docs/source/running_cpp/parameters.rst | 9 +++++++++ Source/Particles/PhysicalParticleContainer.cpp | 3 +++ 2 files changed, 12 insertions(+) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Docs/source/running_cpp/parameters.rst b/Docs/source/running_cpp/parameters.rst index ac378dc03..2ab071673 100644 --- a/Docs/source/running_cpp/parameters.rst +++ b/Docs/source/running_cpp/parameters.rst @@ -187,6 +187,15 @@ Particle initialization * ``NRandomPerCell``: injection with a fixed number of randomly-distributed particles per cell. This requires the additional parameter ``.num_particles_per_cell``. + * ``gaussian_beam``: Inject particle beam with gaussian distribution in + space in all directions. This requires additional parameters: + ``.q_tot`` (beam charge), + ``.npart`` (number of particles in the beam), + ``.x/y/z_m`` (average position in `x/y/z`), + ``.x/y/z_rms`` (standard deviation in `x/y/z`), + and optional argument ``.do_symmetrize`` (whether to + symmetrize the beam in the x and y directions). + * ``.do_continuous_injection`` (`0` or `1`) Whether to inject particles during the simulation, and not only at initialization. This can be required whith a moving window and/or when diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 38ea55d6e..bf93375fe 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -244,6 +244,7 @@ PhysicalParticleContainer::CheckAndAddParticle(Real x, Real y, Real z, std::array attribs; attribs.fill(0.0); + // update attribs with input arguments if (WarpX::gamma_boost > 1.) { MapParticletoBoostedFrame(x, y, z, u); } @@ -254,6 +255,7 @@ PhysicalParticleContainer::CheckAndAddParticle(Real x, Real y, Real z, if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { + // need to create old values auto& particle_tile = DefineAndReturnParticleTile(0, 0, 0); particle_tile.push_back_real(particle_comps["xold"], x); particle_tile.push_back_real(particle_comps["yold"], y); @@ -263,6 +265,7 @@ PhysicalParticleContainer::CheckAndAddParticle(Real x, Real y, Real z, particle_tile.push_back_real(particle_comps["uyold"], u[1]); particle_tile.push_back_real(particle_comps["uzold"], u[2]); } + // add particle AddOneParticle(0, 0, 0, x, y, z, attribs); } -- cgit v1.2.3 From b7dbfc6b22fb3d94310cc5f5e4b2f34fccfe49ce Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Mon, 27 May 2019 07:33:26 -0700 Subject: typo --- Source/Particles/PhysicalParticleContainer.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index bf93375fe..c762bdbb3 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -224,7 +224,7 @@ PhysicalParticleContainer::AddGaussianBeam(Real x_m, Real y_m, Real z_m, u_tmp[0] *= std::pow(-1,ix); y *= std::pow(-1,iy); u_tmp[1] *= std::pow(-1,iy); - CheckAndAddParticle(x, y, z, u_tmp, weight); + CheckAndAddParticle(x, y, z, u_tmp, weight/4); } } } else { -- cgit v1.2.3 From 2ddfcb10aeb31b496153658269147e203a362657 Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Tue, 28 May 2019 17:02:14 -0700 Subject: fix issue for symmetric beam --- Source/Particles/PhysicalParticleContainer.cpp | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index c762bdbb3..2c501985f 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -215,14 +215,16 @@ PhysicalParticleContainer::AddGaussianBeam(Real x_m, Real y_m, Real z_m, if (plasma_injector->insideBounds(x, y, z)) { plasma_injector->getMomentum(u, x, y, z); if (do_symmetrize){ + std::array u_tmp; + Real x_tmp, y_tmp, z_tmp; // Add four particles to the beam: // (x,ux,y,uy) (-x,-ux,y,uy) (x,ux,-y,-uy) (-x,-ux,-y,-uy) - std::array u_tmp = u; for (int ix=0; ix<2; ix++){ for (int iy=0; iy<2; iy++){ - x *= std::pow(-1,ix); + u_tmp = u; + x_tmp = x*std::pow(-1,ix); u_tmp[0] *= std::pow(-1,ix); - y *= std::pow(-1,iy); + y_tmp = y*std::pow(-1,iy); u_tmp[1] *= std::pow(-1,iy); CheckAndAddParticle(x, y, z, u_tmp, weight/4); } -- cgit v1.2.3 From 6c20c43e0ecc9a03b21ef6c57adb2f13911004e2 Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Tue, 28 May 2019 17:19:53 -0700 Subject: typos --- Source/Particles/PhysicalParticleContainer.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 2c501985f..6faf7a127 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -216,7 +216,7 @@ PhysicalParticleContainer::AddGaussianBeam(Real x_m, Real y_m, Real z_m, plasma_injector->getMomentum(u, x, y, z); if (do_symmetrize){ std::array u_tmp; - Real x_tmp, y_tmp, z_tmp; + Real x_tmp, y_tmp; // Add four particles to the beam: // (x,ux,y,uy) (-x,-ux,y,uy) (x,ux,-y,-uy) (-x,-ux,-y,-uy) for (int ix=0; ix<2; ix++){ @@ -226,7 +226,8 @@ PhysicalParticleContainer::AddGaussianBeam(Real x_m, Real y_m, Real z_m, u_tmp[0] *= std::pow(-1,ix); y_tmp = y*std::pow(-1,iy); u_tmp[1] *= std::pow(-1,iy); - CheckAndAddParticle(x, y, z, u_tmp, weight/4); + CheckAndAddParticle(x_tmp, y_tmp, z, + u_tmp, weight/4); } } } else { -- cgit v1.2.3