diff options
Diffstat (limited to 'Source')
32 files changed, 640 insertions, 546 deletions
diff --git a/Source/Diagnostics/BTDiagnostics.cpp b/Source/Diagnostics/BTDiagnostics.cpp index 2b4efcb5e..84a0501b5 100644 --- a/Source/Diagnostics/BTDiagnostics.cpp +++ b/Source/Diagnostics/BTDiagnostics.cpp @@ -15,12 +15,12 @@ #include "Diagnostics/FlushFormats/FlushFormat.H" #include "ComputeDiagFunctors/BackTransformParticleFunctor.H" #include "Utils/Algorithms/IsIn.H" -#include "Utils/CoarsenIO.H" #include "Utils/Parser/ParserUtils.H" #include "Utils/TextMsg.H" #include "Utils/WarpXConst.H" #include "WarpX.H" +#include <ablastr/coarsen/sample.H> #include <ablastr/utils/Communication.H> #include <ablastr/utils/SignalHandling.H> #include <ablastr/warn_manager/WarnManager.H> @@ -772,8 +772,8 @@ BTDiagnostics::PrepareFieldDataForOutput () // Flattening out MF over levels for (int lev = warpx.finestLevel(); lev > 0; --lev) { - CoarsenIO::Coarsen( *m_cell_centered_data[lev-1], *m_cell_centered_data[lev], 0, 0, - m_cellcenter_varnames.size(), 0, WarpX::RefRatio(lev-1) ); + ablastr::coarsen::sample::Coarsen(*m_cell_centered_data[lev - 1], *m_cell_centered_data[lev], 0, 0, + m_cellcenter_varnames.size(), 0, WarpX::RefRatio(lev-1) ); } int num_BT_functors = 1; diff --git a/Source/Diagnostics/ComputeDiagFunctors/CellCenterFunctor.cpp b/Source/Diagnostics/ComputeDiagFunctors/CellCenterFunctor.cpp index 2dac5fb00..f43714a9c 100644 --- a/Source/Diagnostics/ComputeDiagFunctors/CellCenterFunctor.cpp +++ b/Source/Diagnostics/ComputeDiagFunctors/CellCenterFunctor.cpp @@ -1,11 +1,12 @@ #include "CellCenterFunctor.H" -#include "Utils/CoarsenIO.H" #include "Utils/TextMsg.H" #ifdef WARPX_DIM_RZ # include "WarpX.H" #endif +#include <ablastr/coarsen/sample.H> + #include <AMReX.H> #include <AMReX_IntVect.H> #include <AMReX_MultiFab.H> @@ -35,14 +36,14 @@ CellCenterFunctor::operator()(amrex::MultiFab& mf_dst, int dcomp, const int /*i_ // All modes > 0 amrex::MultiFab::Add(mf_dst_stag, *m_mf_src, ic, 0, 1, m_mf_src->nGrowVect()); } - CoarsenIO::Coarsen( mf_dst, mf_dst_stag, dcomp, 0, nComp(), 0, m_crse_ratio); + ablastr::coarsen::sample::Coarsen( mf_dst, mf_dst_stag, dcomp, 0, nComp(), 0, m_crse_ratio); } else { - CoarsenIO::Coarsen( mf_dst, *m_mf_src, dcomp, 0, nComp(), 0, m_crse_ratio); + ablastr::coarsen::sample::Coarsen( mf_dst, *m_mf_src, dcomp, 0, nComp(), 0, m_crse_ratio); } #else // In cartesian geometry, coarsen and interpolate from simulation MultiFab, m_mf_src, // to output diagnostic MultiFab, mf_dst. - CoarsenIO::Coarsen( mf_dst, *m_mf_src, dcomp, 0, nComp(), mf_dst.nGrowVect(), m_crse_ratio); + ablastr::coarsen::sample::Coarsen(mf_dst, *m_mf_src, dcomp, 0, nComp(), mf_dst.nGrowVect(), m_crse_ratio); amrex::ignore_unused(m_lev, m_convertRZmodes2cartesian); #endif } diff --git a/Source/Diagnostics/ComputeDiagFunctors/DivBFunctor.cpp b/Source/Diagnostics/ComputeDiagFunctors/DivBFunctor.cpp index 2dd401a28..093b4960e 100644 --- a/Source/Diagnostics/ComputeDiagFunctors/DivBFunctor.cpp +++ b/Source/Diagnostics/ComputeDiagFunctors/DivBFunctor.cpp @@ -1,8 +1,9 @@ #include "DivBFunctor.H" -#include "Utils/CoarsenIO.H" #include "WarpX.H" +#include <ablastr/coarsen/sample.H> + #include <AMReX_IntVect.H> #include <AMReX_MultiFab.H> @@ -25,7 +26,7 @@ DivBFunctor::operator()(amrex::MultiFab& mf_dst, int dcomp, const int /*i_buffer amrex::MultiFab divB( warpx.boxArray(m_lev), warpx.DistributionMap(m_lev), warpx.ncomps, ng ); warpx.ComputeDivB(divB, 0, m_arr_mf_src, WarpX::CellSize(m_lev) ); // // Coarsen and Interpolate from divB to coarsened/reduced_domain mf_dst - // CoarsenIO::Coarsen( mf_dst, divB, dcomp, 0, nComp(), 0, m_crse_ratio); + // ablastr::coarsen::sample::Coarsen( mf_dst, divB, dcomp, 0, nComp(), 0, m_crse_ratio); #ifdef WARPX_DIM_RZ if (m_convertRZmodes2cartesian) { // In cylindrical geometry, sum real part of all modes of divE in @@ -41,15 +42,15 @@ DivBFunctor::operator()(amrex::MultiFab& mf_dst, int dcomp, const int /*i_buffer amrex::MultiFab::Add(mf_dst_stag, divB, ic, 0, 1, divB.nGrowVect()); } // TODO check if coarsening is needed, otherwise copy - CoarsenIO::Coarsen( mf_dst, mf_dst_stag, dcomp, 0, nComp(), 0, m_crse_ratio); + ablastr::coarsen::sample::Coarsen( mf_dst, mf_dst_stag, dcomp, 0, nComp(), 0, m_crse_ratio); } else { // TODO check if coarsening is needed, otherwise copy - CoarsenIO::Coarsen( mf_dst, divB, dcomp, 0, nComp(), 0, m_crse_ratio); + ablastr::coarsen::sample::Coarsen( mf_dst, divB, dcomp, 0, nComp(), 0, m_crse_ratio); } #else // In cartesian geometry, coarsen and interpolate from simulation MultiFab, divE, // to output diagnostic MultiFab, mf_dst. - CoarsenIO::Coarsen( mf_dst, divB, dcomp, 0, nComp(), 0, m_crse_ratio); + ablastr::coarsen::sample::Coarsen(mf_dst, divB, dcomp, 0, nComp(), 0, m_crse_ratio); amrex::ignore_unused(m_convertRZmodes2cartesian); #endif } diff --git a/Source/Diagnostics/ComputeDiagFunctors/DivEFunctor.cpp b/Source/Diagnostics/ComputeDiagFunctors/DivEFunctor.cpp index 3859b859c..0fd525d4f 100644 --- a/Source/Diagnostics/ComputeDiagFunctors/DivEFunctor.cpp +++ b/Source/Diagnostics/ComputeDiagFunctors/DivEFunctor.cpp @@ -1,12 +1,13 @@ #include "DivEFunctor.H" -#include "Utils/CoarsenIO.H" #include "Utils/TextMsg.H" #ifdef WARPX_DIM_RZ # include "Utils/WarpXAlgorithmSelection.H" #endif #include "WarpX.H" +#include <ablastr/coarsen/sample.H> + #include <AMReX.H> #include <AMReX_BoxArray.H> #include <AMReX_IntVect.H> @@ -55,14 +56,14 @@ DivEFunctor::operator()(amrex::MultiFab& mf_dst, const int dcomp, const int /*i_ // Real part of all modes > 0 amrex::MultiFab::Add(mf_dst_stag, divE, ic, 0, 1, divE.nGrowVect()); } - CoarsenIO::Coarsen( mf_dst, mf_dst_stag, dcomp, 0, nComp(), 0, m_crse_ratio); + ablastr::coarsen::sample::Coarsen( mf_dst, mf_dst_stag, dcomp, 0, nComp(), 0, m_crse_ratio); } else { - CoarsenIO::Coarsen( mf_dst, divE, dcomp, 0, nComp(), 0, m_crse_ratio); + ablastr::coarsen::sample::Coarsen( mf_dst, divE, dcomp, 0, nComp(), 0, m_crse_ratio); } #else // In cartesian geometry, coarsen and interpolate from simulation MultiFab, divE, // to output diagnostic MultiFab, mf_dst. - CoarsenIO::Coarsen( mf_dst, divE, dcomp, 0, nComp(), 0, m_crse_ratio); + ablastr::coarsen::sample::Coarsen(mf_dst, divE, dcomp, 0, nComp(), 0, m_crse_ratio); amrex::ignore_unused(m_convertRZmodes2cartesian); #endif } diff --git a/Source/Diagnostics/ComputeDiagFunctors/PartPerCellFunctor.cpp b/Source/Diagnostics/ComputeDiagFunctors/PartPerCellFunctor.cpp index dd663cfb3..493c52e2e 100644 --- a/Source/Diagnostics/ComputeDiagFunctors/PartPerCellFunctor.cpp +++ b/Source/Diagnostics/ComputeDiagFunctors/PartPerCellFunctor.cpp @@ -2,9 +2,10 @@ #include "Diagnostics/ComputeDiagFunctors/ComputeDiagFunctor.H" #include "Particles/MultiParticleContainer.H" -#include "Utils/CoarsenIO.H" #include "WarpX.H" +#include <ablastr/coarsen/sample.H> + #include <AMReX_BLassert.H> #include <AMReX_IntVect.H> #include <AMReX_MultiFab.H> @@ -36,5 +37,5 @@ PartPerCellFunctor::operator()(amrex::MultiFab& mf_dst, const int dcomp, const i // Compute ppc which includes a summation over all species. warpx.GetPartContainer().Increment(ppc_mf, m_lev); // Coarsen and interpolate from ppc_mf to the output diagnostic MultiFab, mf_dst. - CoarsenIO::Coarsen(mf_dst, ppc_mf, dcomp, 0, nComp(), 0, m_crse_ratio); + ablastr::coarsen::sample::Coarsen(mf_dst, ppc_mf, dcomp, 0, nComp(), 0, m_crse_ratio); } diff --git a/Source/Diagnostics/ComputeDiagFunctors/PartPerGridFunctor.cpp b/Source/Diagnostics/ComputeDiagFunctors/PartPerGridFunctor.cpp index 45a0ad115..5c637f2c7 100644 --- a/Source/Diagnostics/ComputeDiagFunctors/PartPerGridFunctor.cpp +++ b/Source/Diagnostics/ComputeDiagFunctors/PartPerGridFunctor.cpp @@ -2,9 +2,10 @@ #include "Diagnostics/ComputeDiagFunctors/ComputeDiagFunctor.H" #include "Particles/MultiParticleContainer.H" -#include "Utils/CoarsenIO.H" #include "WarpX.H" +#include <ablastr/coarsen/sample.H> + #include <AMReX_BLassert.H> #include <AMReX_Config.H> #include <AMReX_FArrayBox.H> @@ -48,5 +49,5 @@ PartPerGridFunctor::operator()(amrex::MultiFab& mf_dst, const int dcomp, const i } // Coarsen and interpolate from ppg_mf to the output diagnostic MultiFab, mf_dst. - CoarsenIO::Coarsen(mf_dst, ppg_mf, dcomp, 0, nComp(), 0, m_crse_ratio); + ablastr::coarsen::sample::Coarsen(mf_dst, ppg_mf, dcomp, 0, nComp(), 0, m_crse_ratio); } diff --git a/Source/Diagnostics/ComputeDiagFunctors/ParticleReductionFunctor.cpp b/Source/Diagnostics/ComputeDiagFunctors/ParticleReductionFunctor.cpp index 52952a373..bdc524898 100644 --- a/Source/Diagnostics/ComputeDiagFunctors/ParticleReductionFunctor.cpp +++ b/Source/Diagnostics/ComputeDiagFunctors/ParticleReductionFunctor.cpp @@ -4,9 +4,10 @@ #include "Diagnostics/ComputeDiagFunctors/ComputeDiagFunctor.H" #include "Particles/MultiParticleContainer.H" #include "Particles/WarpXParticleContainer.H" -#include "Utils/CoarsenIO.H" #include "WarpX.H" +#include <ablastr/coarsen/sample.H> + #include <AMReX_Array.H> #include <AMReX_BLassert.H> #include <AMReX_IntVect.H> @@ -153,5 +154,5 @@ ParticleReductionFunctor::operator() (amrex::MultiFab& mf_dst, const int dcomp, } // Coarsen and interpolate from ppc_mf to the output diagnostic MultiFab, mf_dst. - CoarsenIO::Coarsen(mf_dst, red_mf, dcomp, 0, nComp(), 0, m_crse_ratio); + ablastr::coarsen::sample::Coarsen(mf_dst, red_mf, dcomp, 0, nComp(), 0, m_crse_ratio); } diff --git a/Source/Diagnostics/ComputeDiagFunctors/RhoFunctor.cpp b/Source/Diagnostics/ComputeDiagFunctors/RhoFunctor.cpp index de9257722..6f4fcda4d 100644 --- a/Source/Diagnostics/ComputeDiagFunctors/RhoFunctor.cpp +++ b/Source/Diagnostics/ComputeDiagFunctors/RhoFunctor.cpp @@ -8,10 +8,11 @@ #endif #include "Particles/MultiParticleContainer.H" #include "Particles/WarpXParticleContainer.H" -#include "Utils/CoarsenIO.H" #include "Utils/TextMsg.H" #include "WarpX.H" +#include <ablastr/coarsen/sample.H> + #include <AMReX.H> #include <AMReX_IntVect.H> #include <AMReX_MultiFab.H> @@ -83,14 +84,14 @@ RhoFunctor::operator() ( amrex::MultiFab& mf_dst, const int dcomp, const int /*i // Real part of all modes > 0 amrex::MultiFab::Add( mf_dst_stag, *rho, ic, 0, 1, rho->nGrowVect() ); } - CoarsenIO::Coarsen( mf_dst, mf_dst_stag, dcomp, 0, nComp(), 0, m_crse_ratio ); + ablastr::coarsen::sample::Coarsen( mf_dst, mf_dst_stag, dcomp, 0, nComp(), 0, m_crse_ratio ); } else { - CoarsenIO::Coarsen( mf_dst, *rho, dcomp, 0, nComp(), 0, m_crse_ratio ); + ablastr::coarsen::sample::Coarsen( mf_dst, *rho, dcomp, 0, nComp(), 0, m_crse_ratio ); } #else // In Cartesian geometry, coarsen and interpolate from temporary MultiFab rho // to output diagnostic MultiFab mf_dst - CoarsenIO::Coarsen( mf_dst, *rho, dcomp, 0, nComp(), mf_dst.nGrowVect(), m_crse_ratio ); + ablastr::coarsen::sample::Coarsen(mf_dst, *rho, dcomp, 0, nComp(), mf_dst.nGrowVect(), m_crse_ratio ); amrex::ignore_unused(m_convertRZmodes2cartesian); #endif } diff --git a/Source/Diagnostics/FieldIO.cpp b/Source/Diagnostics/FieldIO.cpp index 0f2b86050..0b7cc9356 100644 --- a/Source/Diagnostics/FieldIO.cpp +++ b/Source/Diagnostics/FieldIO.cpp @@ -7,9 +7,10 @@ */ #include "FieldIO.H" -#include "Utils/CoarsenIO.H" #include "Utils/TextMsg.H" +#include <ablastr/coarsen/sample.H> + #include <AMReX.H> #include <AMReX_IntVect.H> #include <AMReX_MultiFab.H> @@ -183,9 +184,9 @@ AverageAndPackVectorField( MultiFab& mf_avg, const std::array<std::unique_ptr<MultiFab>,3> &vector_total = vector_field; #endif - CoarsenIO::Coarsen( mf_avg, *(vector_total[0]), dcomp , 0, 1, ngrow ); - CoarsenIO::Coarsen( mf_avg, *(vector_total[1]), dcomp+1, 0, 1, ngrow ); - CoarsenIO::Coarsen( mf_avg, *(vector_total[2]), dcomp+2, 0, 1, ngrow ); + ablastr::coarsen::sample::Coarsen(mf_avg, *(vector_total[0]), dcomp , 0, 1, ngrow ); + ablastr::coarsen::sample::Coarsen(mf_avg, *(vector_total[1]), dcomp + 1, 0, 1, ngrow ); + ablastr::coarsen::sample::Coarsen(mf_avg, *(vector_total[2]), dcomp + 2, 0, 1, ngrow ); } /** \brief Take a MultiFab `scalar_field` @@ -220,7 +221,7 @@ AverageAndPackScalarField (MultiFab& mf_avg, MultiFab::Copy( mf_avg, *scalar_total, 0, dcomp, 1, ngrow); } else if ( scalar_total->is_nodal() ){ // - Fully nodal - CoarsenIO::Coarsen( mf_avg, *scalar_total, dcomp, 0, 1, ngrow ); + ablastr::coarsen::sample::Coarsen(mf_avg, *scalar_total, dcomp, 0, 1, ngrow ); } else { amrex::Abort(Utils::TextMsg::Err("Unknown staggering.")); } diff --git a/Source/Diagnostics/ReducedDiags/FieldMaximum.cpp b/Source/Diagnostics/ReducedDiags/FieldMaximum.cpp index 804d1641a..a36c349dd 100644 --- a/Source/Diagnostics/ReducedDiags/FieldMaximum.cpp +++ b/Source/Diagnostics/ReducedDiags/FieldMaximum.cpp @@ -7,10 +7,11 @@ #include "FieldMaximum.H" -#include "Utils/CoarsenIO.H" #include "Utils/TextMsg.H" #include "WarpX.H" +#include <ablastr/coarsen/sample.H> + #include <AMReX_Algorithm.H> #include <AMReX_Array.H> #include <AMReX_Array4.H> @@ -192,65 +193,65 @@ void FieldMaximum::ComputeDiags (int step) reduceEx_op.eval(box, reduceEx_data, [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple { - const Real Ex_interp = CoarsenIO::Interp(arrEx, Extype, cellCenteredtype, - reduction_coarsening_ratio, i, j, k, reduction_comp); + const Real Ex_interp = ablastr::coarsen::sample::Interp(arrEx, Extype, cellCenteredtype, + reduction_coarsening_ratio, i, j, k, reduction_comp); return amrex::Math::abs(Ex_interp); }); reduceEy_op.eval(box, reduceEy_data, [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple { - const Real Ey_interp = CoarsenIO::Interp(arrEy, Eytype, cellCenteredtype, - reduction_coarsening_ratio, i, j, k, reduction_comp); + const Real Ey_interp = ablastr::coarsen::sample::Interp(arrEy, Eytype, cellCenteredtype, + reduction_coarsening_ratio, i, j, k, reduction_comp); return amrex::Math::abs(Ey_interp); }); reduceEz_op.eval(box, reduceEz_data, [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple { - const Real Ez_interp = CoarsenIO::Interp(arrEz, Eztype, cellCenteredtype, - reduction_coarsening_ratio, i, j, k, reduction_comp); + const Real Ez_interp = ablastr::coarsen::sample::Interp(arrEz, Eztype, cellCenteredtype, + reduction_coarsening_ratio, i, j, k, reduction_comp); return amrex::Math::abs(Ez_interp); }); reduceBx_op.eval(box, reduceBx_data, [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple { - const Real Bx_interp = CoarsenIO::Interp(arrBx, Bxtype, cellCenteredtype, - reduction_coarsening_ratio, i, j, k, reduction_comp); + const Real Bx_interp = ablastr::coarsen::sample::Interp(arrBx, Bxtype, cellCenteredtype, + reduction_coarsening_ratio, i, j, k, reduction_comp); return amrex::Math::abs(Bx_interp); }); reduceBy_op.eval(box, reduceBy_data, [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple { - const Real By_interp = CoarsenIO::Interp(arrBy, Bytype, cellCenteredtype, - reduction_coarsening_ratio, i, j, k, reduction_comp); + const Real By_interp = ablastr::coarsen::sample::Interp(arrBy, Bytype, cellCenteredtype, + reduction_coarsening_ratio, i, j, k, reduction_comp); return amrex::Math::abs(By_interp); }); reduceBz_op.eval(box, reduceBz_data, [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple { - const Real Bz_interp = CoarsenIO::Interp(arrBz, Bztype, cellCenteredtype, - reduction_coarsening_ratio, i, j, k, reduction_comp); + const Real Bz_interp = ablastr::coarsen::sample::Interp(arrBz, Bztype, cellCenteredtype, + reduction_coarsening_ratio, i, j, k, reduction_comp); return amrex::Math::abs(Bz_interp); }); reduceE_op.eval(box, reduceE_data, [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple { - const Real Ex_interp = CoarsenIO::Interp(arrEx, Extype, cellCenteredtype, - reduction_coarsening_ratio, i, j, k, reduction_comp); - const Real Ey_interp = CoarsenIO::Interp(arrEy, Eytype, cellCenteredtype, - reduction_coarsening_ratio, i, j, k, reduction_comp); - const Real Ez_interp = CoarsenIO::Interp(arrEz, Eztype, cellCenteredtype, - reduction_coarsening_ratio, i, j, k, reduction_comp); + const Real Ex_interp = ablastr::coarsen::sample::Interp(arrEx, Extype, cellCenteredtype, + reduction_coarsening_ratio, i, j, k, reduction_comp); + const Real Ey_interp = ablastr::coarsen::sample::Interp(arrEy, Eytype, cellCenteredtype, + reduction_coarsening_ratio, i, j, k, reduction_comp); + const Real Ez_interp = ablastr::coarsen::sample::Interp(arrEz, Eztype, cellCenteredtype, + reduction_coarsening_ratio, i, j, k, reduction_comp); return Ex_interp*Ex_interp + Ey_interp*Ey_interp + Ez_interp*Ez_interp; }); reduceB_op.eval(box, reduceB_data, [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple { - const Real Bx_interp = CoarsenIO::Interp(arrBx, Bxtype, cellCenteredtype, - reduction_coarsening_ratio, i, j, k, reduction_comp); - const Real By_interp = CoarsenIO::Interp(arrBy, Bytype, cellCenteredtype, - reduction_coarsening_ratio, i, j, k, reduction_comp); - const Real Bz_interp = CoarsenIO::Interp(arrBz, Bztype, cellCenteredtype, - reduction_coarsening_ratio, i, j, k, reduction_comp); + const Real Bx_interp = ablastr::coarsen::sample::Interp(arrBx, Bxtype, cellCenteredtype, + reduction_coarsening_ratio, i, j, k, reduction_comp); + const Real By_interp = ablastr::coarsen::sample::Interp(arrBy, Bytype, cellCenteredtype, + reduction_coarsening_ratio, i, j, k, reduction_comp); + const Real Bz_interp = ablastr::coarsen::sample::Interp(arrBz, Bztype, cellCenteredtype, + reduction_coarsening_ratio, i, j, k, reduction_comp); return Bx_interp*Bx_interp + By_interp*By_interp + Bz_interp*Bz_interp; }); } diff --git a/Source/Diagnostics/ReducedDiags/FieldMomentum.cpp b/Source/Diagnostics/ReducedDiags/FieldMomentum.cpp index 45a5cc6cb..8ae51b0e6 100644 --- a/Source/Diagnostics/ReducedDiags/FieldMomentum.cpp +++ b/Source/Diagnostics/ReducedDiags/FieldMomentum.cpp @@ -7,11 +7,12 @@ #include "FieldMomentum.H" -#include "Utils/CoarsenIO.H" #include "Utils/TextMsg.H" #include "Utils/WarpXConst.H" #include "WarpX.H" +#include <ablastr/coarsen/sample.H> + #include <AMReX_Array.H> #include <AMReX_Array4.H> #include <AMReX_Box.H> @@ -162,13 +163,13 @@ void FieldMomentum::ComputeDiags (int step) reduce_ops.eval(box, reduce_data, [=] AMREX_GPU_DEVICE (int i, int j, int k) -> amrex::GpuTuple<Real, Real, Real> { - const amrex::Real Ex_cc = CoarsenIO::Interp(Ex_arr, Ex_stag, cc, cr, i, j, k, comp); - const amrex::Real Ey_cc = CoarsenIO::Interp(Ey_arr, Ey_stag, cc, cr, i, j, k, comp); - const amrex::Real Ez_cc = CoarsenIO::Interp(Ez_arr, Ez_stag, cc, cr, i, j, k, comp); + const amrex::Real Ex_cc = ablastr::coarsen::sample::Interp(Ex_arr, Ex_stag, cc, cr, i, j, k, comp); + const amrex::Real Ey_cc = ablastr::coarsen::sample::Interp(Ey_arr, Ey_stag, cc, cr, i, j, k, comp); + const amrex::Real Ez_cc = ablastr::coarsen::sample::Interp(Ez_arr, Ez_stag, cc, cr, i, j, k, comp); - const amrex::Real Bx_cc = CoarsenIO::Interp(Bx_arr, Bx_stag, cc, cr, i, j, k, comp); - const amrex::Real By_cc = CoarsenIO::Interp(By_arr, By_stag, cc, cr, i, j, k, comp); - const amrex::Real Bz_cc = CoarsenIO::Interp(Bz_arr, Bz_stag, cc, cr, i, j, k, comp); + const amrex::Real Bx_cc = ablastr::coarsen::sample::Interp(Bx_arr, Bx_stag, cc, cr, i, j, k, comp); + const amrex::Real By_cc = ablastr::coarsen::sample::Interp(By_arr, By_stag, cc, cr, i, j, k, comp); + const amrex::Real Bz_cc = ablastr::coarsen::sample::Interp(Bz_arr, Bz_stag, cc, cr, i, j, k, comp); return {Ey_cc * Bz_cc - Ez_cc * By_cc, Ez_cc * Bx_cc - Ex_cc * Bz_cc, diff --git a/Source/Diagnostics/ReducedDiags/FieldProbeParticleContainer.cpp b/Source/Diagnostics/ReducedDiags/FieldProbeParticleContainer.cpp index 35e25fb4b..54d31ed5a 100644 --- a/Source/Diagnostics/ReducedDiags/FieldProbeParticleContainer.cpp +++ b/Source/Diagnostics/ReducedDiags/FieldProbeParticleContainer.cpp @@ -12,7 +12,6 @@ #include "Particles/Pusher/GetAndSetPosition.H" #include "Particles/Pusher/UpdatePosition.H" #include "Particles/ParticleBoundaries_K.H" -#include "Utils/CoarsenMR.H" #include "Utils/TextMsg.H" #include "Utils/WarpXAlgorithmSelection.H" #include "Utils/WarpXConst.H" diff --git a/Source/Diagnostics/ReducedDiags/FieldReduction.H b/Source/Diagnostics/ReducedDiags/FieldReduction.H index b6b918946..ca1b766d3 100644 --- a/Source/Diagnostics/ReducedDiags/FieldReduction.H +++ b/Source/Diagnostics/ReducedDiags/FieldReduction.H @@ -9,9 +9,10 @@ #define WARPX_DIAGNOSTICS_REDUCEDDIAGS_FIELDREDUCTION_H_ #include "ReducedDiags.H" -#include "Utils/CoarsenIO.H" #include "WarpX.H" +#include <ablastr/coarsen/sample.H> + #include <AMReX_Array.H> #include <AMReX_Box.H> #include <AMReX_Config.H> @@ -157,18 +158,18 @@ public: const amrex::Real y = (j + 0.5_rt)*dx[1] + real_box.lo(1); const amrex::Real z = (k + 0.5_rt)*dx[2] + real_box.lo(2); #endif - const amrex::Real Ex_interp = CoarsenIO::Interp(arrEx, Extype, cellCenteredtype, - reduction_coarsening_ratio, i, j, k, reduction_comp); - const amrex::Real Ey_interp = CoarsenIO::Interp(arrEy, Eytype, cellCenteredtype, - reduction_coarsening_ratio, i, j, k, reduction_comp); - const amrex::Real Ez_interp = CoarsenIO::Interp(arrEz, Eztype, cellCenteredtype, - reduction_coarsening_ratio, i, j, k, reduction_comp); - const amrex::Real Bx_interp = CoarsenIO::Interp(arrBx, Bxtype, cellCenteredtype, - reduction_coarsening_ratio, i, j, k, reduction_comp); - const amrex::Real By_interp = CoarsenIO::Interp(arrBy, Bytype, cellCenteredtype, - reduction_coarsening_ratio, i, j, k, reduction_comp); - const amrex::Real Bz_interp = CoarsenIO::Interp(arrBz, Bztype, cellCenteredtype, - reduction_coarsening_ratio, i, j, k, reduction_comp); + const amrex::Real Ex_interp = ablastr::coarsen::sample::Interp(arrEx, Extype, cellCenteredtype, + reduction_coarsening_ratio, i, j, k, reduction_comp); + const amrex::Real Ey_interp = ablastr::coarsen::sample::Interp(arrEy, Eytype, cellCenteredtype, + reduction_coarsening_ratio, i, j, k, reduction_comp); + const amrex::Real Ez_interp = ablastr::coarsen::sample::Interp(arrEz, Eztype, cellCenteredtype, + reduction_coarsening_ratio, i, j, k, reduction_comp); + const amrex::Real Bx_interp = ablastr::coarsen::sample::Interp(arrBx, Bxtype, cellCenteredtype, + reduction_coarsening_ratio, i, j, k, reduction_comp); + const amrex::Real By_interp = ablastr::coarsen::sample::Interp(arrBy, Bytype, cellCenteredtype, + reduction_coarsening_ratio, i, j, k, reduction_comp); + const amrex::Real Bz_interp = ablastr::coarsen::sample::Interp(arrBz, Bztype, cellCenteredtype, + reduction_coarsening_ratio, i, j, k, reduction_comp); return reduction_function_parser(x, y, z, Ex_interp, Ey_interp, Ez_interp, Bx_interp, By_interp, Bz_interp); }); diff --git a/Source/Diagnostics/WarpXIO.cpp b/Source/Diagnostics/WarpXIO.cpp index 1785a04c1..65d71442f 100644 --- a/Source/Diagnostics/WarpXIO.cpp +++ b/Source/Diagnostics/WarpXIO.cpp @@ -13,7 +13,6 @@ #endif #include "FieldIO.H" #include "Particles/MultiParticleContainer.H" -#include "Utils/CoarsenIO.H" #include "Utils/TextMsg.H" #include "Utils/WarpXProfilerWrapper.H" #include "WarpX.H" diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/FieldAccessorFunctors.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/FieldAccessorFunctors.H index 58c0837a8..d4fdf207e 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/FieldAccessorFunctors.H +++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/FieldAccessorFunctors.H @@ -9,10 +9,11 @@ #define WARPX_FIELD_ACCESSOR_FUNCTORS_H #include "WarpX.H" -#include "Utils/CoarsenIO.H" #include "FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H" + #include <AMReX_Array.H> + /** * \brief Functor that returns the division of the source m_field Array4 value by macroparameter obtained using m_parameter, at the respective (i,j,k). diff --git a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveE.cpp b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveE.cpp index ab2e25014..22a222726 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveE.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveE.cpp @@ -8,11 +8,12 @@ # include "FiniteDifferenceAlgorithms/FieldAccessorFunctors.H" #endif #include "MacroscopicProperties/MacroscopicProperties.H" -#include "Utils/CoarsenIO.H" #include "Utils/TextMsg.H" #include "Utils/WarpXAlgorithmSelection.H" #include "WarpX.H" +#include <ablastr/coarsen/sample.H> + #include <AMReX.H> #include <AMReX_Array4.H> #include <AMReX_Config.H> @@ -112,7 +113,7 @@ void FiniteDifferenceSolver::MacroscopicEvolveECartesian ( amrex::MultiFab& epsilon_mf = macroscopic_properties->getepsilon_mf(); amrex::MultiFab& mu_mf = macroscopic_properties->getmu_mf(); - // Index type required for calling CoarsenIO::Interp to interpolate macroscopic + // Index type required for calling ablastr::coarsen::sample::Interp to interpolate macroscopic // properties from their respective staggering to the Ex, Ey, Ez locations amrex::GpuArray<int, 3> const& sigma_stag = macroscopic_properties->sigma_IndexType; amrex::GpuArray<int, 3> const& epsilon_stag = macroscopic_properties->epsilon_IndexType; @@ -178,11 +179,11 @@ void FiniteDifferenceSolver::MacroscopicEvolveECartesian ( if (lx(i, j, k) <= 0) return; #endif // Interpolate conductivity, sigma, to Ex position on the grid - amrex::Real const sigma_interp = CoarsenIO::Interp( sigma_arr, sigma_stag, - Ex_stag, macro_cr, i, j, k, scomp); + amrex::Real const sigma_interp = ablastr::coarsen::sample::Interp(sigma_arr, sigma_stag, + Ex_stag, macro_cr, i, j, k, scomp); // Interpolated permittivity, epsilon, to Ex position on the grid - amrex::Real const epsilon_interp = CoarsenIO::Interp( eps_arr, epsilon_stag, - Ex_stag, macro_cr, i, j, k, scomp); + amrex::Real const epsilon_interp = ablastr::coarsen::sample::Interp(eps_arr, epsilon_stag, + Ex_stag, macro_cr, i, j, k, scomp); amrex::Real alpha = T_MacroAlgo::alpha( sigma_interp, epsilon_interp, dt); amrex::Real beta = T_MacroAlgo::beta( sigma_interp, epsilon_interp, dt); Ex(i, j, k) = alpha * Ex(i, j, k) @@ -202,11 +203,11 @@ void FiniteDifferenceSolver::MacroscopicEvolveECartesian ( #endif #endif // Interpolate conductivity, sigma, to Ey position on the grid - amrex::Real const sigma_interp = CoarsenIO::Interp( sigma_arr, sigma_stag, - Ey_stag, macro_cr, i, j, k, scomp); + amrex::Real const sigma_interp = ablastr::coarsen::sample::Interp(sigma_arr, sigma_stag, + Ey_stag, macro_cr, i, j, k, scomp); // Interpolated permittivity, epsilon, to Ey position on the grid - amrex::Real const epsilon_interp = CoarsenIO::Interp( eps_arr, epsilon_stag, - Ey_stag, macro_cr, i, j, k, scomp); + amrex::Real const epsilon_interp = ablastr::coarsen::sample::Interp(eps_arr, epsilon_stag, + Ey_stag, macro_cr, i, j, k, scomp); amrex::Real alpha = T_MacroAlgo::alpha( sigma_interp, epsilon_interp, dt); amrex::Real beta = T_MacroAlgo::beta( sigma_interp, epsilon_interp, dt); @@ -222,11 +223,11 @@ void FiniteDifferenceSolver::MacroscopicEvolveECartesian ( if (lz(i,j,k) <= 0) return; #endif // Interpolate conductivity, sigma, to Ez position on the grid - amrex::Real const sigma_interp = CoarsenIO::Interp( sigma_arr, sigma_stag, - Ez_stag, macro_cr, i, j, k, scomp); + amrex::Real const sigma_interp = ablastr::coarsen::sample::Interp(sigma_arr, sigma_stag, + Ez_stag, macro_cr, i, j, k, scomp); // Interpolated permittivity, epsilon, to Ez position on the grid - amrex::Real const epsilon_interp = CoarsenIO::Interp( eps_arr, epsilon_stag, - Ez_stag, macro_cr, i, j, k, scomp); + amrex::Real const epsilon_interp = ablastr::coarsen::sample::Interp(eps_arr, epsilon_stag, + Ez_stag, macro_cr, i, j, k, scomp); amrex::Real alpha = T_MacroAlgo::alpha( sigma_interp, epsilon_interp, dt); amrex::Real beta = T_MacroAlgo::beta( sigma_interp, epsilon_interp, dt); diff --git a/Source/Parallelization/WarpXComm.cpp b/Source/Parallelization/WarpXComm.cpp index fed81f2fc..8585da632 100644 --- a/Source/Parallelization/WarpXComm.cpp +++ b/Source/Parallelization/WarpXComm.cpp @@ -13,13 +13,13 @@ # include "BoundaryConditions/PML_RZ.H" #endif #include "Filter/BilinearFilter.H" -#include "Utils/CoarsenMR.H" #include "Utils/TextMsg.H" #include "Utils/WarpXAlgorithmSelection.H" #include "Utils/WarpXProfilerWrapper.H" #include "WarpXComm_K.H" #include "WarpXSumGuardCells.H" +#include <ablastr/coarsen/average.H> #include <ablastr/utils/Communication.H> #include <AMReX.H> @@ -887,9 +887,9 @@ WarpX::SyncCurrent ( std::array< MultiFab*,3> crse { J_cp[lev][0].get(), J_cp[lev][1].get(), J_cp[lev][2].get() }; - CoarsenMR::Coarsen( *crse[0], *fine[0], refinement_ratio ); - CoarsenMR::Coarsen( *crse[1], *fine[1], refinement_ratio ); - CoarsenMR::Coarsen( *crse[2], *fine[2], refinement_ratio ); + ablastr::coarsen::average::Coarsen(*crse[0], *fine[0], refinement_ratio ); + ablastr::coarsen::average::Coarsen(*crse[1], *fine[1], refinement_ratio ); + ablastr::coarsen::average::Coarsen(*crse[2], *fine[2], refinement_ratio ); } // For each level @@ -915,7 +915,7 @@ WarpX::SyncRho () { rho_cp[lev]->setVal(0.0); const IntVect& refinement_ratio = refRatio(lev-1); - CoarsenMR::Coarsen( *rho_cp[lev], *rho_fp[lev], refinement_ratio ); + ablastr::coarsen::average::Coarsen(*rho_cp[lev], *rho_fp[lev], refinement_ratio ); } // For each level @@ -947,9 +947,9 @@ void WarpX::RestrictCurrentFromFineToCoarsePatch ( std::array< MultiFab*,3> crse { J_cp[lev][0].get(), J_cp[lev][1].get(), J_cp[lev][2].get() }; - CoarsenMR::Coarsen( *crse[0], *fine[0], refinement_ratio ); - CoarsenMR::Coarsen( *crse[1], *fine[1], refinement_ratio ); - CoarsenMR::Coarsen( *crse[2], *fine[2], refinement_ratio ); + ablastr::coarsen::average::Coarsen(*crse[0], *fine[0], refinement_ratio ); + ablastr::coarsen::average::Coarsen(*crse[1], *fine[1], refinement_ratio ); + ablastr::coarsen::average::Coarsen(*crse[2], *fine[2], refinement_ratio ); } void WarpX::ApplyFilterJ ( @@ -1126,7 +1126,7 @@ void WarpX::RestrictRhoFromFineToCoarsePatch ( if (charge_fp[lev]) { charge_cp[lev]->setVal(0.0); const IntVect& refinement_ratio = refRatio(lev-1); - CoarsenMR::Coarsen( *charge_cp[lev], *charge_fp[lev], refinement_ratio ); + ablastr::coarsen::average::Coarsen(*charge_cp[lev], *charge_fp[lev], refinement_ratio ); } } diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 34e25d1a2..e56c94583 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -15,13 +15,13 @@ #include "Pusher/GetAndSetPosition.H" #include "Pusher/UpdatePosition.H" #include "ParticleBoundaries_K.H" -#include "Utils/CoarsenMR.H" #include "Utils/TextMsg.H" #include "Utils/WarpXAlgorithmSelection.H" #include "Utils/WarpXConst.H" #include "Utils/WarpXProfilerWrapper.H" #include "WarpX.H" +#include <ablastr/coarsen/average.H> #include <ablastr/utils/Communication.H> #include <AMReX.H> @@ -714,7 +714,7 @@ WarpXParticleContainer::DepositCharge (amrex::Vector<std::unique_ptr<amrex::Mult MultiFab coarsened_fine_data(coarsened_fine_BA, fine_dm, rho[lev+1]->nComp(), ngrow ); coarsened_fine_data.setVal(0.0); - CoarsenMR::Coarsen( coarsened_fine_data, *rho[lev+1], m_gdb->refRatio(lev) ); + ablastr::coarsen::average::Coarsen(coarsened_fine_data, *rho[lev + 1], m_gdb->refRatio(lev) ); ablastr::utils::communication::ParallelAdd(*rho[lev], coarsened_fine_data, 0, 0, rho[lev]->nComp(), amrex::IntVect::TheZeroVector(), diff --git a/Source/Utils/CMakeLists.txt b/Source/Utils/CMakeLists.txt index 19fe6bc34..0cbd98780 100644 --- a/Source/Utils/CMakeLists.txt +++ b/Source/Utils/CMakeLists.txt @@ -1,7 +1,5 @@ target_sources(WarpX PRIVATE - CoarsenIO.cpp - CoarsenMR.cpp Interpolate.cpp MPIInitHelpers.cpp ParticleUtils.cpp diff --git a/Source/Utils/CoarsenIO.cpp b/Source/Utils/CoarsenIO.cpp deleted file mode 100644 index 7357dc923..000000000 --- a/Source/Utils/CoarsenIO.cpp +++ /dev/null @@ -1,148 +0,0 @@ -#include "CoarsenIO.H" - -#include "Utils/TextMsg.H" - -#include <AMReX_BLProfiler.H> -#include <AMReX_BLassert.H> -#include <AMReX_Box.H> -#include <AMReX_BoxArray.H> -#include <AMReX_Config.H> -#include <AMReX_DistributionMapping.H> -#include <AMReX_FArrayBox.H> -#include <AMReX_FabArray.H> -#include <AMReX_GpuControl.H> -#include <AMReX_GpuLaunch.H> -#include <AMReX_IndexType.H> -#include <AMReX_MFIter.H> -#include <AMReX_MultiFab.H> - -using namespace amrex; - -void -CoarsenIO::Loop ( MultiFab& mf_dst, - const MultiFab& mf_src, - const int dcomp, - const int scomp, - const int ncomp, - const IntVect ngrowvect, - const IntVect crse_ratio ) -{ - // Staggering of source fine MultiFab and destination coarse MultiFab - const IntVect stag_src = mf_src.boxArray().ixType().toIntVect(); - const IntVect stag_dst = mf_dst.boxArray().ixType().toIntVect(); - - if ( crse_ratio > IntVect(1) ) WARPX_ALWAYS_ASSERT_WITH_MESSAGE( ngrowvect == IntVect(0), - "option of filling guard cells of destination MultiFab with coarsening not supported for this interpolation" ); - - WARPX_ALWAYS_ASSERT_WITH_MESSAGE( mf_src.nGrowVect() >= stag_dst-stag_src+ngrowvect, - "source fine MultiFab does not have enough guard cells for this interpolation" ); - - // Auxiliary integer arrays (always 3D) - GpuArray<int,3> sf; // staggering of source fine MultiFab - GpuArray<int,3> sc; // staggering of destination coarse MultiFab - GpuArray<int,3> cr; // coarsening ratio - - sf[0] = stag_src[0]; -#if defined(WARPX_DIM_1D_Z) - sf[1] = 0; -#else - sf[1] = stag_src[1]; -#endif -#if (AMREX_SPACEDIM <= 2) - sf[2] = 0; -#elif defined(WARPX_DIM_3D) - sf[2] = stag_src[2]; -#endif - - sc[0] = stag_dst[0]; -#if defined(WARPX_DIM_1D_Z) - sc[1] = 0; -#else - sc[1] = stag_dst[1]; -#endif -#if (AMREX_SPACEDIM <= 2) - sc[2] = 0; -#elif defined(WARPX_DIM_3D) - sc[2] = stag_dst[2]; -#endif - - cr[0] = crse_ratio[0]; -#if defined(WARPX_DIM_1D_Z) - cr[1] = 1; -#else - cr[1] = crse_ratio[1]; -#endif -#if (AMREX_SPACEDIM <= 2) - cr[2] = 1; -#elif defined(WARPX_DIM_3D) - cr[2] = crse_ratio[2]; -#endif - -#ifdef AMREX_USE_OMP -#pragma omp parallel if (Gpu::notInLaunchRegion()) -#endif - // Loop over boxes (or tiles if not on GPU) - for (MFIter mfi( mf_dst, TilingIfNotGPU() ); mfi.isValid(); ++mfi) - { - // Tiles defined at the coarse level - const Box& bx = mfi.growntilebox( ngrowvect ); - Array4<Real> const& arr_dst = mf_dst.array( mfi ); - Array4<Real const> const& arr_src = mf_src.const_array( mfi ); - ParallelFor( bx, ncomp, - [=] AMREX_GPU_DEVICE( int i, int j, int k, int n ) - { - arr_dst(i,j,k,n+dcomp) = CoarsenIO::Interp( - arr_src, sf, sc, cr, i, j, k, n+scomp ); - } ); - } -} - -void -CoarsenIO::Coarsen ( MultiFab& mf_dst, - const MultiFab& mf_src, - const int dcomp, - const int scomp, - const int ncomp, - const int ngrow, - const IntVect crse_ratio ) -{ - amrex::IntVect ngrowvect(ngrow); - Coarsen(mf_dst, - mf_src, - dcomp, - scomp, - ncomp, - ngrowvect, - crse_ratio); -} - -void -CoarsenIO::Coarsen ( MultiFab& mf_dst, - const MultiFab& mf_src, - const int dcomp, - const int scomp, - const int ncomp, - const IntVect ngrowvect, - const IntVect crse_ratio ) -{ - BL_PROFILE("CoarsenIO::Coarsen()"); - - // Convert BoxArray of source MultiFab to staggering of destination MultiFab and coarsen it - BoxArray ba_tmp = amrex::convert( mf_src.boxArray(), mf_dst.ixType().toIntVect() ); - WARPX_ALWAYS_ASSERT_WITH_MESSAGE( ba_tmp.coarsenable( crse_ratio ), - "source MultiFab converted to staggering of destination MultiFab is not coarsenable" ); - ba_tmp.coarsen( crse_ratio ); - - if ( ba_tmp == mf_dst.boxArray() and mf_src.DistributionMap() == mf_dst.DistributionMap() ) - CoarsenIO::Loop( mf_dst, mf_src, dcomp, scomp, ncomp, ngrowvect, crse_ratio ); - else - { - // Cannot coarsen into MultiFab with different BoxArray or DistributionMapping: - // 1) create temporary MultiFab on coarsened version of source BoxArray with same DistributionMapping - MultiFab mf_tmp( ba_tmp, mf_src.DistributionMap(), ncomp, ngrowvect, MFInfo(), FArrayBoxFactory() ); - // 2) interpolate from mf_src to mf_tmp (start writing into component 0) - CoarsenIO::Loop( mf_tmp, mf_src, 0, scomp, ncomp, ngrowvect, crse_ratio ); - // 3) copy from mf_tmp to mf_dst (with different BoxArray or DistributionMapping) - mf_dst.ParallelCopy( mf_tmp, 0, dcomp, ncomp ); - } -} diff --git a/Source/Utils/CoarsenMR.H b/Source/Utils/CoarsenMR.H deleted file mode 100644 index 1191b0c7e..000000000 --- a/Source/Utils/CoarsenMR.H +++ /dev/null @@ -1,154 +0,0 @@ -#ifndef WARPX_COARSEN_MR_H_ -#define WARPX_COARSEN_MR_H_ - -#include <AMReX_Array.H> -#include <AMReX_Array4.H> -#include <AMReX_Extension.H> -#include <AMReX_GpuQualifiers.H> -#include <AMReX_REAL.H> - -#include <AMReX_BaseFwd.H> - -#include <cstdlib> - -namespace CoarsenMR{ - - using namespace amrex; - - /** - * \brief Interpolates the floating point data contained in the source Array4 - * \c arr_src, extracted from a fine MultiFab, with weights defined in - * such a way that the total charge is preserved. - * - * \param[in] arr_src floating point data to be interpolated - * \param[in] sf staggering of the source fine MultiFab - * \param[in] sc staggering of the destination coarsened MultiFab - * \param[in] cr coarsening ratio along each spatial direction - * \param[in] i index along x of the coarsened Array4 to be filled - * \param[in] j index along y of the coarsened Array4 to be filled - * \param[in] k index along z of the coarsened Array4 to be filled - * \param[in] comp index along the fourth component of the Array4 \c arr_src - * containing the data to be interpolated - * - * \return interpolated field at cell (i,j,k) of a coarsened Array4 - */ - AMREX_GPU_DEVICE - AMREX_FORCE_INLINE - Real Interp ( Array4<Real const> const& arr_src, - GpuArray<int,3> const& sf, - GpuArray<int,3> const& sc, - GpuArray<int,3> const& cr, - const int i, - const int j, - const int k, - const int comp ) - { - // Indices of destination array (coarse) - const int ic[3] = { i, j, k }; - - // Number of points and starting indices of source array (fine) - int np[3], idx_min[3]; - - // Compute number of points - for ( int l = 0; l < 3; ++l ) { - if ( cr[l] == 1 ) np[l] = 1; // no coarsening - else np[l] = cr[l]*(1-sf[l])*(1-sc[l]) // cell-centered - +(2*(cr[l]-1)+1)*sf[l]*sc[l]; // nodal - } - - // Compute starting indices of source array (fine) - for ( int l = 0; l < 3; ++l ) { - if ( cr[l] == 1 ) idx_min[l] = ic[l]; // no coarsening - else idx_min[l] = ic[l]*cr[l]*(1-sf[l])*(1-sc[l]) // cell-centered - +(ic[l]*cr[l]-cr[l]+1)*sf[l]*sc[l]; // nodal - } - - // Auxiliary integer variables - const int numx = np[0]; - const int numy = np[1]; - const int numz = np[2]; - const int imin = idx_min[0]; - const int jmin = idx_min[1]; - const int kmin = idx_min[2]; - const int sfx = sf[0]; - const int sfy = sf[1]; - const int sfz = sf[2]; - const int scx = sc[0]; - const int scy = sc[1]; - const int scz = sc[2]; - const int crx = cr[0]; - const int cry = cr[1]; - const int crz = cr[2]; - int ii, jj, kk; - Real wx, wy, wz; - - // Add neutral elements (=0) beyond guard cells in source array (fine) - auto const arr_src_safe = [arr_src] - AMREX_GPU_DEVICE (int const ix, int const iy, int const iz, int const n) noexcept - { - return arr_src.contains( ix, iy, iz ) ? arr_src(ix,iy,iz,n) : 0.0_rt; - }; - - // Interpolate over points computed above. Weights are computed in order - // to guarantee total charge conservation for both cell-centered data - // (equal weights) and nodal data (weights depend on distance between - // points on fine and coarse grids). Terms multiplied by (1-sf)*(1-sc) - // are ON for cell-centered data and OFF for nodal data, while terms - // multiplied by sf*sc are ON for nodal data and OFF for cell-centered data. - // Python script Source/Utils/check_interp_points_and_weights.py can be - // used to check interpolation points and weights in 1D. - Real c = 0.0_rt; - for (int kref = 0; kref < numz; ++kref) { - for (int jref = 0; jref < numy; ++jref) { - for (int iref = 0; iref < numx; ++iref) { - ii = imin+iref; - jj = jmin+jref; - kk = kmin+kref; - wx = (1.0_rt/static_cast<Real>(numx))*(1-sfx)*(1-scx) // if cell-centered - +((amrex::Math::abs(crx-amrex::Math::abs(ii-i*crx)))/static_cast<Real>(crx*crx))*sfx*scx; // if nodal - wy = (1.0_rt/static_cast<Real>(numy))*(1-sfy)*(1-scy) // if cell-centered - +((amrex::Math::abs(cry-amrex::Math::abs(jj-j*cry)))/static_cast<Real>(cry*cry))*sfy*scy; // if nodal - wz = (1.0_rt/static_cast<Real>(numz))*(1-sfz)*(1-scz) // if cell-centered - +((amrex::Math::abs(crz-amrex::Math::abs(kk-k*crz)))/static_cast<Real>(crz*crz))*sfz*scz; // if nodal - c += wx*wy*wz*arr_src_safe(ii,jj,kk,comp); - } - } - } - return c; - } - - /** - * \brief Loops over the boxes of the coarsened MultiFab \c mf_dst and fills - * them by interpolating the data contained in the fine MultiFab \c mf_src. - * - * \param[in,out] mf_dst coarsened MultiFab containing the floating point data - * to be filled by interpolating the source fine MultiFab - * \param[in] mf_src fine MultiFab containing the floating point data to be interpolated - * \param[in] ncomp number of components to loop over for the coarsened - * Array4 extracted from the coarsened MultiFab \c mf_dst - * \param[in] ngrow number of guard cells to fill along each spatial direction - * \param[in] crse_ratio coarsening ratio between the fine MultiFab \c mf_src - * and the coarsened MultiFab \c mf_dst along each spatial direction - */ - void Loop ( MultiFab& mf_dst, - const MultiFab& mf_src, - const int ncomp, - const IntVect ngrow, - const IntVect crse_ratio ); - - /** - * \brief Stores in the coarsened MultiFab \c mf_dst the values obtained by - * interpolating the data contained in the fine MultiFab \c mf_src. - * - * \param[in,out] mf_dst coarsened MultiFab containing the floating point data - * to be filled by interpolating the fine MultiFab \c mf_src - * \param[in] mf_src fine MultiFab containing the floating point data to be interpolated - * \param[in] crse_ratio coarsening ratio between the fine MultiFab \c mf_src - * and the coarsened MultiFab \c mf_dst along each spatial direction - */ - void Coarsen ( MultiFab& mf_dst, - const MultiFab& mf_src, - const IntVect crse_ratio ); -} - -#endif // WARPX_COARSEN_MR_H_ diff --git a/Source/Utils/CoarsenMR.cpp b/Source/Utils/CoarsenMR.cpp deleted file mode 100644 index 549ff6ea2..000000000 --- a/Source/Utils/CoarsenMR.cpp +++ /dev/null @@ -1,104 +0,0 @@ -#include "CoarsenMR.H" - -#include "Utils/TextMsg.H" - -#include <AMReX_BLProfiler.H> -#include <AMReX_BLassert.H> -#include <AMReX_BoxArray.H> -#include <AMReX_Config.H> -#include <AMReX_GpuControl.H> -#include <AMReX_GpuLaunch.H> -#include <AMReX_IndexType.H> -#include <AMReX_IntVect.H> -#include <AMReX_MFIter.H> -#include <AMReX_MultiFab.H> - -using namespace amrex; - -void -CoarsenMR::Loop ( MultiFab& mf_dst, - const MultiFab& mf_src, - const int ncomp, - const IntVect ngrow, - const IntVect crse_ratio ) -{ - // Staggering of source fine MultiFab and destination coarse MultiFab - const IntVect stag_src = mf_src.boxArray().ixType().toIntVect(); - const IntVect stag_dst = mf_dst.boxArray().ixType().toIntVect(); - - // Auxiliary integer arrays (always 3D) - GpuArray<int,3> sf; // staggering of source fine MultiFab - GpuArray<int,3> sc; // staggering of destination coarse MultiFab - GpuArray<int,3> cr; // coarsening ratio - - sf[0] = stag_src[0]; -#if defined(WARPX_DIM_1D_Z) - sf[1] = 0; -#else - sf[1] = stag_src[1]; -#endif -#if (AMREX_SPACEDIM <= 2) - sf[2] = 0; -#elif defined(WARPX_DIM_3D) - sf[2] = stag_src[2]; -#endif - - sc[0] = stag_dst[0]; -#if defined(WARPX_DIM_1D_Z) - sc[1] = 0; -#else - sc[1] = stag_dst[1]; -#endif -#if (AMREX_SPACEDIM <= 2) - sc[2] = 0; -#elif defined(WARPX_DIM_3D) - sc[2] = stag_dst[2]; -#endif - - cr[0] = crse_ratio[0]; -#if defined(WARPX_DIM_1D_Z) - cr[1] = 1; -#else - cr[1] = crse_ratio[1]; -#endif -#if (AMREX_SPACEDIM <= 2) - cr[2] = 1; -#elif defined(WARPX_DIM_3D) - cr[2] = crse_ratio[2]; -#endif - -#ifdef AMREX_USE_OMP -#pragma omp parallel if (Gpu::notInLaunchRegion()) -#endif - // Loop over boxes (or tiles if not on GPU) - for (MFIter mfi( mf_dst, TilingIfNotGPU() ); mfi.isValid(); ++mfi) - { - // Tiles defined at the coarse level - const Box& bx = mfi.growntilebox( ngrow ); - Array4<Real> const& arr_dst = mf_dst.array( mfi ); - Array4<Real const> const& arr_src = mf_src.const_array( mfi ); - ParallelFor( bx, ncomp, - [=] AMREX_GPU_DEVICE( int i, int j, int k, int n ) - { - arr_dst(i,j,k,n) = CoarsenMR::Interp( - arr_src, sf, sc, cr, i, j, k, n ); - } ); - } -} - -void -CoarsenMR::Coarsen ( MultiFab& mf_dst, - const MultiFab& mf_src, - const IntVect crse_ratio ) -{ - BL_PROFILE("CoarsenMR::Coarsen()"); - - WARPX_ALWAYS_ASSERT_WITH_MESSAGE( mf_src.ixType() == mf_dst.ixType(), - "source MultiFab and destination MultiFab have different IndexType" ); - - // Number of guard cells to fill on coarse patch and number of components - const IntVect ngrow = ( mf_src.nGrowVect() + 1 ) / crse_ratio; - const int ncomp = mf_src.nComp(); - - CoarsenMR::Loop( mf_dst, mf_src, ncomp, ngrow, crse_ratio ); -} diff --git a/Source/Utils/Make.package b/Source/Utils/Make.package index f43953720..3d6da0de9 100644 --- a/Source/Utils/Make.package +++ b/Source/Utils/Make.package @@ -3,8 +3,6 @@ CEXE_sources += WarpXTagging.cpp CEXE_sources += WarpXUtil.cpp CEXE_sources += WarpXVersion.cpp CEXE_sources += WarpXAlgorithmSelection.cpp -CEXE_sources += CoarsenIO.cpp -CEXE_sources += CoarsenMR.cpp CEXE_sources += Interpolate.cpp CEXE_sources += IntervalsParser.cpp CEXE_sources += MPIInitHelpers.cpp diff --git a/Source/Utils/check_interp_points_and_weights.py b/Source/Utils/check_interp_points_and_weights.py index 4d54af01c..a1d17c8dd 100644 --- a/Source/Utils/check_interp_points_and_weights.py +++ b/Source/Utils/check_interp_points_and_weights.py @@ -21,7 +21,8 @@ # For MR applications only the cases sc=sf=0 and sc=sf=1 are considered. Terms # multiplied by (1-sf)*(1-sc) are ON for cell-centered data and OFF for nodal data, # while terms multiplied by sf*sc are ON for nodal data and OFF for cell-centered -# data. C++ implementation in Source/Utils/CoarsenMR.H/.cpp and Source/Utils/CoarsenIO.H/.cpp +# data. C++ implementation in Source/ablastr/coarsen/average.(H/.cpp) and +# Source/ablastr/coarsen/sample.(H/.cpp) #------------------------------------------------------------------------------- import sys diff --git a/Source/ablastr/CMakeLists.txt b/Source/ablastr/CMakeLists.txt index 224145a98..0224f1a8b 100644 --- a/Source/ablastr/CMakeLists.txt +++ b/Source/ablastr/CMakeLists.txt @@ -1,4 +1,5 @@ #add_subdirectory(fields) +add_subdirectory(coarsen) #add_subdirectory(particles) #add_subdirectory(profiler) add_subdirectory(utils) diff --git a/Source/ablastr/Make.package b/Source/ablastr/Make.package index 46b3f1858..b10d9f629 100644 --- a/Source/ablastr/Make.package +++ b/Source/ablastr/Make.package @@ -1,5 +1,6 @@ #CEXE_sources += ParticleBoundaries.cpp +include $(WARPX_HOME)/Source/ablastr/coarsen/Make.package include $(WARPX_HOME)/Source/ablastr/particles/Make.package include $(WARPX_HOME)/Source/ablastr/utils/Make.package include $(WARPX_HOME)/Source/ablastr/warn_manager/Make.package diff --git a/Source/ablastr/coarsen/CMakeLists.txt b/Source/ablastr/coarsen/CMakeLists.txt new file mode 100644 index 000000000..7396c8d8a --- /dev/null +++ b/Source/ablastr/coarsen/CMakeLists.txt @@ -0,0 +1,5 @@ +target_sources(ablastr + PRIVATE + average.cpp + sample.cpp +) diff --git a/Source/ablastr/coarsen/Make.package b/Source/ablastr/coarsen/Make.package new file mode 100644 index 000000000..4fcd0e2ec --- /dev/null +++ b/Source/ablastr/coarsen/Make.package @@ -0,0 +1,4 @@ +CEXE_sources += average.cpp +CEXE_sources += sample.cpp + +VPATH_LOCATIONS += $(WARPX_HOME)/Source/ablastr/coarsen diff --git a/Source/ablastr/coarsen/average.H b/Source/ablastr/coarsen/average.H new file mode 100644 index 000000000..269403f7b --- /dev/null +++ b/Source/ablastr/coarsen/average.H @@ -0,0 +1,191 @@ +/* Copyright 2022 Edoardo Zoni, Remi Lehe, Prabhat Kumar, Axel Huebl + * + * This file is part of ABLASTR. + * + * License: BSD-3-Clause-LBNL + */ +#ifndef ABLASTR_COARSEN_AVERAGE_H_ +#define ABLASTR_COARSEN_AVERAGE_H_ + + +#include <AMReX_Array.H> +#include <AMReX_Array4.H> +#include <AMReX_BLassert.H> +#include <AMReX_Extension.H> +#include <AMReX_GpuQualifiers.H> +#include <AMReX_IntVect.H> +#include <AMReX_Math.H> +#include <AMReX_REAL.H> + +#include <AMReX_BaseFwd.H> + +#include <cstdlib> + + +/** Mesh Coarsening by Averaging + * + * These methods are mostly used for mesh-refinement. + */ +namespace ablastr::coarsen::average +{ + /** + * \brief Interpolates the floating point data contained in the source Array4 + * \c arr_src, extracted from a fine MultiFab, with weights defined in + * such a way that the total charge is preserved. + * + * The input (sf) and output (sc) staggering need to be the same. + * + * \param[in] arr_src floating point data to be interpolated + * \param[in] sf staggering of the source fine MultiFab + * \param[in] sc staggering of the destination coarsened MultiFab + * \param[in] cr coarsening ratio along each spatial direction + * \param[in] i index along x of the coarsened Array4 to be filled + * \param[in] j index along y of the coarsened Array4 to be filled + * \param[in] k index along z of the coarsened Array4 to be filled + * \param[in] comp index along the fourth component of the Array4 \c arr_src + * containing the data to be interpolated + * + * \return interpolated field at cell (i,j,k) of a coarsened Array4 + */ + AMREX_GPU_DEVICE + AMREX_FORCE_INLINE + amrex::Real + Interp ( + amrex::Array4<amrex::Real const> const &arr_src, + amrex::GpuArray<int, 3> const &sf, + amrex::GpuArray<int, 3> const &sc, + amrex::GpuArray<int, 3> const &cr, + int const i, + int const j, + int const k, + int const comp + ) + { + using namespace amrex::literals; + + AMREX_ASSERT_WITH_MESSAGE(sf[0] == sc[0], "Interp: Staggering for component 0 does not match!"); + AMREX_ASSERT_WITH_MESSAGE(sf[1] == sc[1], "Interp: Staggering for component 1 does not match!"); + AMREX_ASSERT_WITH_MESSAGE(sf[2] == sc[2], "Interp: Staggering for component 2 does not match!"); + + // Indices of destination array (coarse) + int const ic[3] = {i, j, k}; + + // Number of points and starting indices of source array (fine) + int np[3], idx_min[3]; + + // Compute number of points + for (int l = 0; l < 3; ++l) { + if (cr[l] == 1) + np[l] = 1; // no coarsening + else + np[l] = cr[l] * (1 - sf[l]) * (1 - sc[l]) // cell-centered + + (2 * (cr[l] - 1) + 1) * sf[l] * sc[l]; // nodal + } + + // Compute starting indices of source array (fine) + for (int l = 0; l < 3; ++l) { + if (cr[l] == 1) + idx_min[l] = ic[l]; // no coarsening + else + idx_min[l] = ic[l] * cr[l] * (1 - sf[l]) * (1 - sc[l]) // cell-centered + + (ic[l] * cr[l] - cr[l] + 1) * sf[l] * sc[l]; // nodal + } + + // Auxiliary integer variables + int const numx = np[0]; + int const numy = np[1]; + int const numz = np[2]; + int const imin = idx_min[0]; + int const jmin = idx_min[1]; + int const kmin = idx_min[2]; + int const sfx = sf[0]; + int const sfy = sf[1]; + int const sfz = sf[2]; + int const scx = sc[0]; + int const scy = sc[1]; + int const scz = sc[2]; + int const crx = cr[0]; + int const cry = cr[1]; + int const crz = cr[2]; + int ii, jj, kk; + amrex::Real wx, wy, wz; + + // Add neutral elements (=0) beyond guard cells in source array (fine) + auto const arr_src_safe = [arr_src] + AMREX_GPU_DEVICE(int const ix, int const iy, int const iz, int const n) noexcept { + return arr_src.contains(ix, iy, iz) ? arr_src(ix, iy, iz, n) : 0.0_rt; + }; + + // Interpolate over points computed above. Weights are computed in order + // to guarantee total charge conservation for both cell-centered data + // (equal weights) and nodal data (weights depend on distance between + // points on fine and coarse grids). Terms multiplied by (1-sf)*(1-sc) + // are ON for cell-centered data and OFF for nodal data, while terms + // multiplied by sf*sc are ON for nodal data and OFF for cell-centered data. + // Python script Source/Utils/check_interp_points_and_weights.py can be + // used to check interpolation points and weights in 1D. + amrex::Real c = 0.0_rt; + for (int kref = 0; kref < numz; ++kref) { + for (int jref = 0; jref < numy; ++jref) { + for (int iref = 0; iref < numx; ++iref) { + ii = imin + iref; + jj = jmin + jref; + kk = kmin + kref; + wx = (1.0_rt / static_cast<amrex::Real>(numx)) * (1 - sfx) * (1 - scx) // if cell-centered + + ((amrex::Math::abs(crx - amrex::Math::abs(ii - i * crx))) / + static_cast<amrex::Real>(crx * crx)) * sfx * scx; // if nodal + wy = (1.0_rt / static_cast<amrex::Real>(numy)) * (1 - sfy) * (1 - scy) // if cell-centered + + ((amrex::Math::abs(cry - amrex::Math::abs(jj - j * cry))) / + static_cast<amrex::Real>(cry * cry)) * sfy * scy; // if nodal + wz = (1.0_rt / static_cast<amrex::Real>(numz)) * (1 - sfz) * (1 - scz) // if cell-centered + + ((amrex::Math::abs(crz - amrex::Math::abs(kk - k * crz))) / + static_cast<amrex::Real>(crz * crz)) * sfz * scz; // if nodal + c += wx * wy * wz * arr_src_safe(ii, jj, kk, comp); + } + } + } + return c; + } + + /** + * \brief Loops over the boxes of the coarsened MultiFab \c mf_dst and fills + * them by interpolating the data contained in the fine MultiFab \c mf_src. + * + * \param[in,out] mf_dst coarsened MultiFab containing the floating point data + * to be filled by interpolating the source fine MultiFab + * \param[in] mf_src fine MultiFab containing the floating point data to be interpolated + * \param[in] ncomp number of components to loop over for the coarsened + * Array4 extracted from the coarsened MultiFab \c mf_dst + * \param[in] ngrow number of guard cells to fill along each spatial direction + * \param[in] crse_ratio coarsening ratio between the fine MultiFab \c mf_src + * and the coarsened MultiFab \c mf_dst along each spatial direction + */ + void + Loop ( + amrex::MultiFab & mf_dst, + amrex::MultiFab const & mf_src, + int const ncomp, + amrex::IntVect const ngrow, + amrex::IntVect const crse_ratio + ); + + /** + * \brief Stores in the coarsened MultiFab \c mf_dst the values obtained by + * interpolating the data contained in the fine MultiFab \c mf_src. + * + * \param[in,out] mf_dst coarsened MultiFab containing the floating point data + * to be filled by interpolating the fine MultiFab \c mf_src + * \param[in] mf_src fine MultiFab containing the floating point data to be interpolated + * \param[in] crse_ratio coarsening ratio between the fine MultiFab \c mf_src + * and the coarsened MultiFab \c mf_dst along each spatial direction + */ + void + Coarsen ( + amrex::MultiFab & mf_dst, + amrex::MultiFab const & mf_src, + amrex::IntVect const crse_ratio + ); + +} // namespace ablastr::coarsen::average + +#endif // ABLASTR_COARSEN_AVERAGE_H_ diff --git a/Source/ablastr/coarsen/average.cpp b/Source/ablastr/coarsen/average.cpp new file mode 100644 index 000000000..2e333867d --- /dev/null +++ b/Source/ablastr/coarsen/average.cpp @@ -0,0 +1,114 @@ +/* Copyright 2022 Edoardo Zoni, Remi Lehe, Prabhat Kumar + * + * This file is part of ABLASTR. + * + * License: BSD-3-Clause-LBNL + */ +#include "average.H" + +#include "ablastr/utils/TextMsg.H" + +#include <AMReX_BLProfiler.H> +#include <AMReX_BLassert.H> +#include <AMReX_BoxArray.H> +#include <AMReX_Config.H> +#include <AMReX_GpuControl.H> +#include <AMReX_GpuLaunch.H> +#include <AMReX_IntVect.H> +#include <AMReX_MFIter.H> +#include <AMReX_MultiFab.H> + + +namespace ablastr::coarsen::average +{ + void + Loop ( + amrex::MultiFab & mf_dst, + amrex::MultiFab const & mf_src, + int const ncomp, + amrex::IntVect const ngrow, + amrex::IntVect const crse_ratio + ) + { + // Staggering of source fine MultiFab and destination coarse MultiFab + amrex::IntVect const stag_src = mf_src.boxArray().ixType().toIntVect(); + amrex::IntVect const stag_dst = mf_dst.boxArray().ixType().toIntVect(); + + // Auxiliary integer arrays (always 3D) + amrex::GpuArray<int, 3> sf; // staggering of source fine MultiFab + amrex::GpuArray<int, 3> sc; // staggering of destination coarse MultiFab + amrex::GpuArray<int, 3> cr; // coarsening ratio + + sf[0] = stag_src[0]; +#if defined(WARPX_DIM_1D_Z) + sf[1] = 0; +#else + sf[1] = stag_src[1]; +#endif +#if (AMREX_SPACEDIM <= 2) + sf[2] = 0; +#elif defined(WARPX_DIM_3D) + sf[2] = stag_src[2]; +#endif + + sc[0] = stag_dst[0]; +#if defined(WARPX_DIM_1D_Z) + sc[1] = 0; +#else + sc[1] = stag_dst[1]; +#endif +#if (AMREX_SPACEDIM <= 2) + sc[2] = 0; +#elif defined(WARPX_DIM_3D) + sc[2] = stag_dst[2]; +#endif + + cr[0] = crse_ratio[0]; +#if defined(WARPX_DIM_1D_Z) + cr[1] = 1; +#else + cr[1] = crse_ratio[1]; +#endif +#if (AMREX_SPACEDIM <= 2) + cr[2] = 1; +#elif defined(WARPX_DIM_3D) + cr[2] = crse_ratio[2]; +#endif + +#ifdef AMREX_USE_OMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + // Loop over boxes (or tiles if not on GPU) + for (amrex::MFIter mfi(mf_dst, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) { + // Tiles defined at the coarse level + amrex::Box const & bx = mfi.growntilebox(ngrow); + amrex::Array4<amrex::Real> const &arr_dst = mf_dst.array(mfi); + amrex::Array4<amrex::Real const> const &arr_src = mf_src.const_array(mfi); + ParallelFor(bx, ncomp, + [=] AMREX_GPU_DEVICE(int i, int j, int k, int n) { + arr_dst(i, j, k, n) = Interp( + arr_src, sf, sc, cr, i, j, k, n); + }); + } + } + + void + Coarsen ( + amrex::MultiFab & mf_dst, + amrex::MultiFab const & mf_src, + amrex::IntVect const crse_ratio + ) + { + BL_PROFILE("ablastr::coarsen::Coarsen()"); + + ABLASTR_ALWAYS_ASSERT_WITH_MESSAGE(mf_src.ixType() == mf_dst.ixType(), + "source MultiFab and destination MultiFab have different IndexType"); + + // Number of guard cells to fill on coarse patch and number of components + const amrex::IntVect ngrow = (mf_src.nGrowVect() + 1) / crse_ratio; + const int ncomp = mf_src.nComp(); + + Loop(mf_dst, mf_src, ncomp, ngrow, crse_ratio); + } + +} // namespace ablastr::coarsen::average diff --git a/Source/Utils/CoarsenIO.H b/Source/ablastr/coarsen/sample.H index 0b53831e6..1390cbebb 100644 --- a/Source/Utils/CoarsenIO.H +++ b/Source/ablastr/coarsen/sample.H @@ -1,21 +1,33 @@ -#ifndef WARPX_COARSEN_IO_H_ -#define WARPX_COARSEN_IO_H_ +/* Copyright 2022 Edoardo Zoni, Remi Lehe, David Grote, Axel Huebl + * + * This file is part of ABLASTR. + * + * License: BSD-3-Clause-LBNL + */ +#ifndef ABLASTR_COARSEN_SAMPLE_H_ +#define ABLASTR_COARSEN_SAMPLE_H_ + #include <AMReX_Array.H> #include <AMReX_Array4.H> +#include <AMReX_BLassert.H> #include <AMReX_Extension.H> #include <AMReX_GpuQualifiers.H> #include <AMReX_IntVect.H> +#include <AMReX_Math.H> #include <AMReX_REAL.H> #include <AMReX_BaseFwd.H> #include <cstdlib> -namespace CoarsenIO{ - - using namespace amrex; +/** Mesh Coarsening by Sampling + * + * These methods are mostly used for I/O. + */ +namespace ablastr::coarsen::sample +{ /** * \brief Interpolates the floating point data contained in the source Array4 * \c arr_src, extracted from a fine MultiFab, by averaging over either @@ -35,15 +47,19 @@ namespace CoarsenIO{ */ AMREX_GPU_DEVICE AMREX_FORCE_INLINE - Real Interp ( Array4<Real const> const& arr_src, - GpuArray<int,3> const& sf, - GpuArray<int,3> const& sc, - GpuArray<int,3> const& cr, - const int i, - const int j, - const int k, - const int comp ) + amrex::Real Interp ( + amrex::Array4<amrex::Real const> const& arr_src, + amrex::GpuArray<int,3> const& sf, + amrex::GpuArray<int,3> const& sc, + amrex::GpuArray<int,3> const& cr, + const int i, + const int j, + const int k, + const int comp + ) { + using namespace amrex::literals; + // Indices of destination array (coarse) const int ic[3] = { i, j, k }; @@ -70,19 +86,18 @@ namespace CoarsenIO{ const int jmin = idx_min[1]; const int kmin = idx_min[2]; int ii, jj, kk; - Real wx, wy, wz; + amrex::Real const wx = 1.0_rt / static_cast<amrex::Real>(numx); + amrex::Real const wy = 1.0_rt / static_cast<amrex::Real>(numy); + amrex::Real const wz = 1.0_rt / static_cast<amrex::Real>(numz); // Interpolate over points computed above - Real c = 0.0_rt; + amrex::Real c = 0.0_rt; for (int kref = 0; kref < numz; ++kref) { for (int jref = 0; jref < numy; ++jref) { for (int iref = 0; iref < numx; ++iref) { ii = imin+iref; jj = jmin+jref; kk = kmin+kref; - wx = 1.0_rt/static_cast<Real>(numx); - wy = 1.0_rt/static_cast<Real>(numy); - wz = 1.0_rt/static_cast<Real>(numz); c += wx*wy*wz*arr_src(ii,jj,kk,comp); } } @@ -109,13 +124,13 @@ namespace CoarsenIO{ * \param[in] crse_ratio coarsening ratio between the fine MultiFab \c mf_src * and the coarsened MultiFab \c mf_dst along each spatial direction */ - void Loop ( MultiFab& mf_dst, - const MultiFab& mf_src, + void Loop ( amrex::MultiFab& mf_dst, + const amrex::MultiFab& mf_src, const int dcomp, const int scomp, const int ncomp, - const IntVect ngrow, - const IntVect crse_ratio=IntVect(1) ); + const amrex::IntVect ngrow, + const amrex::IntVect crse_ratio=amrex::IntVect(1) ); /** * \brief Stores in the coarsened MultiFab \c mf_dst the values obtained by @@ -136,20 +151,21 @@ namespace CoarsenIO{ * \param[in] crse_ratio coarsening ratio between the fine MultiFab \c mf_src * and the coarsened MultiFab \c mf_dst along each spatial direction */ - void Coarsen ( MultiFab& mf_dst, - const MultiFab& mf_src, + void Coarsen ( amrex::MultiFab& mf_dst, + const amrex::MultiFab& mf_src, const int dcomp, const int scomp, const int ncomp, const int ngrow, - const IntVect crse_ratio=IntVect(1) ); - void Coarsen ( MultiFab& mf_dst, - const MultiFab& mf_src, + const amrex::IntVect crse_ratio=amrex::IntVect(1) ); + void Coarsen ( amrex::MultiFab& mf_dst, + const amrex::MultiFab& mf_src, const int dcomp, const int scomp, const int ncomp, - const IntVect ngrowvect, - const IntVect crse_ratio=IntVect(1) ); -} + const amrex::IntVect ngrowvect, + const amrex::IntVect crse_ratio=amrex::IntVect(1) ); + +} // namespace ablastr::coarsen::sample -#endif // WARPX_COARSEN_IO_H_ +#endif // ABLASTR_COARSEN_SAMPLE_H_ diff --git a/Source/ablastr/coarsen/sample.cpp b/Source/ablastr/coarsen/sample.cpp new file mode 100644 index 000000000..e77869017 --- /dev/null +++ b/Source/ablastr/coarsen/sample.cpp @@ -0,0 +1,160 @@ +/* Copyright 2022 Edoardo Zoni, Remi Lehe, David Grote, Axel Huebl + * + * This file is part of ABLASTR. + * + * License: BSD-3-Clause-LBNL + */ +#include "sample.H" + +#include "ablastr/utils/TextMsg.H" + +#include <AMReX_BLProfiler.H> +#include <AMReX_BLassert.H> +#include <AMReX_BoxArray.H> +#include <AMReX_Config.H> +#include <AMReX_GpuControl.H> +#include <AMReX_GpuLaunch.H> +#include <AMReX_IntVect.H> +#include <AMReX_MFIter.H> +#include <AMReX_MultiFab.H> + + +namespace ablastr::coarsen::sample +{ + void + Loop ( + amrex::MultiFab& mf_dst, + const amrex::MultiFab& mf_src, + const int dcomp, + const int scomp, + const int ncomp, + const amrex::IntVect ngrowvect, + const amrex::IntVect crse_ratio + ) + { + // Staggering of source fine MultiFab and destination coarse MultiFab + const amrex::IntVect stag_src = mf_src.boxArray().ixType().toIntVect(); + const amrex::IntVect stag_dst = mf_dst.boxArray().ixType().toIntVect(); + + if ( crse_ratio > amrex::IntVect(1) ) + ABLASTR_ALWAYS_ASSERT_WITH_MESSAGE( ngrowvect == amrex::IntVect(0), + "option of filling guard cells of destination MultiFab with coarsening not supported for this interpolation" ); + + ABLASTR_ALWAYS_ASSERT_WITH_MESSAGE( mf_src.nGrowVect() >= stag_dst-stag_src+ngrowvect, + "source fine MultiFab does not have enough guard cells for this interpolation" ); + + // Auxiliary integer arrays (always 3D) + amrex::GpuArray<int, 3> sf; // staggering of source fine MultiFab + amrex::GpuArray<int, 3> sc; // staggering of destination coarse MultiFab + amrex::GpuArray<int, 3> cr; // coarsening ratio + + sf[0] = stag_src[0]; +#if defined(WARPX_DIM_1D_Z) + sf[1] = 0; +#else + sf[1] = stag_src[1]; +#endif +#if (AMREX_SPACEDIM <= 2) + sf[2] = 0; +#elif defined(WARPX_DIM_3D) + sf[2] = stag_src[2]; +#endif + + sc[0] = stag_dst[0]; +#if defined(WARPX_DIM_1D_Z) + sc[1] = 0; +#else + sc[1] = stag_dst[1]; +#endif +#if (AMREX_SPACEDIM <= 2) + sc[2] = 0; +#elif defined(WARPX_DIM_3D) + sc[2] = stag_dst[2]; +#endif + + cr[0] = crse_ratio[0]; +#if defined(WARPX_DIM_1D_Z) + cr[1] = 1; +#else + cr[1] = crse_ratio[1]; +#endif +#if (AMREX_SPACEDIM <= 2) + cr[2] = 1; +#elif defined(WARPX_DIM_3D) + cr[2] = crse_ratio[2]; +#endif + +#ifdef AMREX_USE_OMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + // Loop over boxes (or tiles if not on GPU) + for (amrex::MFIter mfi( mf_dst, amrex::TilingIfNotGPU() ); mfi.isValid(); ++mfi) + { + // Tiles defined at the coarse level + const amrex::Box& bx = mfi.growntilebox( ngrowvect ); + amrex::Array4<amrex::Real> const& arr_dst = mf_dst.array( mfi ); + amrex::Array4<amrex::Real const> const& arr_src = mf_src.const_array( mfi ); + ParallelFor( bx, ncomp, + [=] AMREX_GPU_DEVICE( int i, int j, int k, int n ) + { + arr_dst(i,j,k,n+dcomp) = Interp( + arr_src, sf, sc, cr, i, j, k, n+scomp ); + } ); + } + } + + void + Coarsen ( + amrex::MultiFab& mf_dst, + const amrex::MultiFab& mf_src, + const int dcomp, + const int scomp, + const int ncomp, + const int ngrow, + const amrex::IntVect crse_ratio + ) + { + amrex::IntVect ngrowvect(ngrow); + Coarsen(mf_dst, + mf_src, + dcomp, + scomp, + ncomp, + ngrowvect, + crse_ratio); + } + + void + Coarsen ( + amrex::MultiFab& mf_dst, + const amrex::MultiFab& mf_src, + const int dcomp, + const int scomp, + const int ncomp, + const amrex::IntVect ngrowvect, + const amrex::IntVect crse_ratio + ) + { + BL_PROFILE("sample::Coarsen()"); + + // Convert BoxArray of source MultiFab to staggering of destination MultiFab and coarsen it + amrex::BoxArray ba_tmp = amrex::convert( mf_src.boxArray(), mf_dst.ixType().toIntVect() ); + ABLASTR_ALWAYS_ASSERT_WITH_MESSAGE( ba_tmp.coarsenable( crse_ratio ), + "source MultiFab converted to staggering of destination MultiFab is not coarsenable" ); + ba_tmp.coarsen( crse_ratio ); + + if ( ba_tmp == mf_dst.boxArray() and mf_src.DistributionMap() == mf_dst.DistributionMap() ) + Loop( mf_dst, mf_src, dcomp, scomp, ncomp, ngrowvect, crse_ratio ); + else + { + // Cannot coarsen into MultiFab with different BoxArray or DistributionMapping: + // 1) create temporary MultiFab on coarsened version of source BoxArray with same DistributionMapping + amrex::MultiFab mf_tmp( ba_tmp, mf_src.DistributionMap(), ncomp, ngrowvect, amrex::MFInfo(), amrex::FArrayBoxFactory() ); + // 2) interpolate from mf_src to mf_tmp (start writing into component 0) + Loop( mf_tmp, mf_src, 0, scomp, ncomp, ngrowvect, crse_ratio ); + // 3) copy from mf_tmp to mf_dst (with different BoxArray or DistributionMapping) + mf_dst.ParallelCopy( mf_tmp, 0, dcomp, ncomp ); + } + } + +} // namespace ablastr::coarsen::sample |