aboutsummaryrefslogtreecommitdiff
path: root/Source
diff options
context:
space:
mode:
Diffstat (limited to 'Source')
-rw-r--r--Source/Diagnostics/ComputeDiagFunctors/CellCenterFunctor.cpp6
-rw-r--r--Source/Diagnostics/ComputeDiagFunctors/DivEFunctor.cpp2
-rw-r--r--Source/Diagnostics/FieldIO.cpp10
-rw-r--r--Source/Diagnostics/WarpXIO.cpp3
-rw-r--r--Source/Utils/Average.H176
-rw-r--r--Source/Utils/Average.cpp109
6 files changed, 213 insertions, 93 deletions
diff --git a/Source/Diagnostics/ComputeDiagFunctors/CellCenterFunctor.cpp b/Source/Diagnostics/ComputeDiagFunctors/CellCenterFunctor.cpp
index f47c9bdcd..f8924a9e6 100644
--- a/Source/Diagnostics/ComputeDiagFunctors/CellCenterFunctor.cpp
+++ b/Source/Diagnostics/ComputeDiagFunctors/CellCenterFunctor.cpp
@@ -27,12 +27,12 @@ CellCenterFunctor::operator()(amrex::MultiFab& mf_dst, int dcomp) const
// All modes > 0
MultiFab::Add(mf_dst_stag, *m_mf_src, ic, 0, 1, m_mf_src->nGrowVect());
}
- Average::ToCellCenter ( mf_dst, mf_dst_stag, dcomp, 0, 0, nComp() );
+ Average::CoarsenAndInterpolate( mf_dst, mf_dst_stag, dcomp, 0, nComp(), 0 );
} else {
- Average::ToCellCenter ( mf_dst, *m_mf_src, dcomp, 0, 0, nComp() );
+ Average::CoarsenAndInterpolate( mf_dst, *m_mf_src, dcomp, 0, nComp(), 0 );
}
#else
// In cartesian geometry, cell-center m_mf_src to mf_dst.
- Average::ToCellCenter ( mf_dst, *m_mf_src, dcomp, 0, 0, nComp() );
+ Average::CoarsenAndInterpolate( mf_dst, *m_mf_src, dcomp, 0, nComp(), 0 );
#endif
}
diff --git a/Source/Diagnostics/ComputeDiagFunctors/DivEFunctor.cpp b/Source/Diagnostics/ComputeDiagFunctors/DivEFunctor.cpp
index 4df506d4a..08e0b6061 100644
--- a/Source/Diagnostics/ComputeDiagFunctors/DivEFunctor.cpp
+++ b/Source/Diagnostics/ComputeDiagFunctors/DivEFunctor.cpp
@@ -16,5 +16,5 @@ DivEFunctor::operator()(amrex::MultiFab& mf_dst, const int dcomp) const
MultiFab divE( ba, warpx.DistributionMap(m_lev), 1, 0 );
warpx.ComputeDivE(divE, m_lev);
- Average::ToCellCenter ( mf_dst, divE, dcomp, 0, 0, 1 );
+ Average::CoarsenAndInterpolate( mf_dst, divE, dcomp, 0, 1, 0 );
}
diff --git a/Source/Diagnostics/FieldIO.cpp b/Source/Diagnostics/FieldIO.cpp
index 4624e5617..770c7cadd 100644
--- a/Source/Diagnostics/FieldIO.cpp
+++ b/Source/Diagnostics/FieldIO.cpp
@@ -273,9 +273,9 @@ AverageAndPackVectorField( MultiFab& mf_avg,
const std::array<std::unique_ptr<MultiFab>,3> &vector_total = vector_field;
#endif
- Average::ToCellCenter( mf_avg, *(vector_total[0]), dcomp , ngrow );
- Average::ToCellCenter( mf_avg, *(vector_total[1]), dcomp+1, ngrow );
- Average::ToCellCenter( mf_avg, *(vector_total[2]), dcomp+2, ngrow );
+ Average::CoarsenAndInterpolate( mf_avg, *(vector_total[0]), dcomp , 0, 1, ngrow );
+ Average::CoarsenAndInterpolate( mf_avg, *(vector_total[1]), dcomp+1, 0, 1, ngrow );
+ Average::CoarsenAndInterpolate( mf_avg, *(vector_total[2]), dcomp+2, 0, 1, ngrow );
}
/** \brief Takes all of the components of the three fields and
@@ -331,7 +331,7 @@ AverageAndPackScalarField (MultiFab& mf_avg,
MultiFab::Copy( mf_avg, *scalar_total, 0, dcomp, 1, ngrow);
} else if ( scalar_total->is_nodal() ){
// - Fully nodal
- Average::ToCellCenter( mf_avg, *scalar_total, dcomp, ngrow, 0, 1 );
+ Average::CoarsenAndInterpolate( mf_avg, *scalar_total, dcomp, 0, 1, ngrow );
} else {
amrex::Abort("Unknown staggering.");
}
@@ -676,7 +676,7 @@ coarsenCellCenteredFields(
BoxArray small_ba = amrex::coarsen(source_mf[lev].boxArray(), coarse_ratio);
coarse_mf.push_back( MultiFab(small_ba, source_mf[lev].DistributionMap(), ncomp, 0) );
- average_down(source_mf[lev], coarse_mf[lev], 0, ncomp, IntVect(coarse_ratio));
+ Average::CoarsenAndInterpolate( coarse_mf[lev], source_mf[lev], 0, 0, ncomp, 0, IntVect(coarse_ratio) );
}
};
diff --git a/Source/Diagnostics/WarpXIO.cpp b/Source/Diagnostics/WarpXIO.cpp
index 4dfe9b3e3..c28fdb713 100644
--- a/Source/Diagnostics/WarpXIO.cpp
+++ b/Source/Diagnostics/WarpXIO.cpp
@@ -10,6 +10,7 @@
#include "WarpX.H"
#include "FieldIO.H"
#include "SliceDiagnostic.H"
+#include "Utils/Average.H"
#ifdef WARPX_USE_OPENPMD
# include "Diagnostics/WarpXOpenPMD.H"
@@ -426,7 +427,7 @@ WarpX::GetCellCenteredData() {
for (int lev = finest_level; lev > 0; --lev)
{
- amrex::average_down(*cc[lev], *cc[lev-1], 0, nc, refRatio(lev-1));
+ Average::CoarsenAndInterpolate( *cc[lev-1], *cc[lev], 0, 0, nc, 0, refRatio(lev-1) );
}
return std::move(cc[0]);
diff --git a/Source/Utils/Average.H b/Source/Utils/Average.H
index bc7d81e66..fb1460413 100644
--- a/Source/Utils/Average.H
+++ b/Source/Utils/Average.H
@@ -8,87 +8,121 @@ namespace Average{
using namespace amrex;
/**
- * \brief Returns the cell-centered average of the floating point data contained
- * in the input
- * <a href="https://amrex-codes.github.io/amrex/doxygen/structamrex_1_1Array4.html">Array4</a>
- * \c mf_in_arr.
+ * \brief Interpolates the floating point data contained in the source Array4
+ * \c mf_fp_arr, extracted from a fine MultiFab, by averaging over either
+ * 1 point or 2 equally distant points.
*
- * \param[in] mf_in_arr floating point data to be averaged
- * \param[in] stag staggering (index type) of the data
- * \param[in] i index along x to access the
- * <a href="https://amrex-codes.github.io/amrex/doxygen/structamrex_1_1Array4.html">Array4</a>
- * \c mf_in_arr containing the data to be averaged
- * \param[in] j index along y to access the
- * <a href="https://amrex-codes.github.io/amrex/doxygen/structamrex_1_1Array4.html">Array4</a>
- * \c mf_in_arr containing the data to be averaged
- * \param[in] k index along z to access the
- * <a href="https://amrex-codes.github.io/amrex/doxygen/structamrex_1_1Array4.html">Array4</a>
- * \c mf_in_arr containing the data to be averaged
- * \param[in] comp index along the fourth component of the
- * <a href="https://amrex-codes.github.io/amrex/doxygen/structamrex_1_1Array4.html">Array4</a>
- * \c mf_in_arr containing the data to be averaged
- * \return averaged field at cell (i,j,k)
+ * \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] np number of points to loop over for interpolation
+ * \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_HOST_DEVICE
AMREX_FORCE_INLINE
- Real ToCellCenter ( Array4<Real const> const& mf_in_arr,
- const IntVect stag,
- const int i,
- const int j,
- const int k,
- const int comp=0 )
+ Real CoarsenAndInterpolateKernel ( Array4<Real const> const& arr_src,
+ int const* const AMREX_RESTRICT sf,
+ int const* const AMREX_RESTRICT sc,
+ int const* const AMREX_RESTRICT cr,
+ int const* const AMREX_RESTRICT np,
+ const int i,
+ const int j,
+ const int k,
+ const int comp )
{
- const int sx = stag[0];
- const int sy = stag[1];
-#if (AMREX_SPACEDIM == 2)
- constexpr int sz = 0;
-#elif (AMREX_SPACEDIM == 3)
- const int sz = stag[2];
-#endif
- return 0.125_rt * ( mf_in_arr(i ,j ,k ,comp)
- + mf_in_arr(i+sx,j ,k ,comp)
- + mf_in_arr(i ,j+sy,k ,comp)
- + mf_in_arr(i ,j ,k+sz,comp)
- + mf_in_arr(i+sx,j+sy,k ,comp)
- + mf_in_arr(i ,j+sy,k+sz,comp)
- + mf_in_arr(i+sx,j ,k+sz,comp)
- + mf_in_arr(i+sx,j+sy,k+sz,comp) );
+ // Indices of source coarse array
+ const int ic[3] = { i, j, k };
+
+ // Index of first point of source Array4 from which interpolation is done
+ int idx_min[3];
+
+ // Compute number of points to loop over (either 1 or 2)
+ // and starting indices of fine array in each direction
+ for ( int l = 0; l < 3; ++l ) {
+ if ( cr[l] == 1 ) idx_min[l] = ic[l]-sc[l]*(1-sf[l]); // no coarsening
+ else idx_min[l] = ic[l]*cr[l]+static_cast<int>(cr[l]/2)*(1-sc[l])-(1-sf[l]);
+ }
+
+ // 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];
+
+ // Interpolate over points computed above
+ 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) {
+ c += arr_src(imin+iref,jmin+jref,kmin+kref,comp);
+ }
+ }
+ }
+ return c / (numx*numy*numz);
};
/**
- * \brief Stores the cell-centered average of the floating point data contained
- * in the input
- * <a href="https://amrex-codes.github.io/amrex/doxygen/classamrex_1_1MultiFab.html">MultiFab</a>
- * \c mf_in into the output
- * <a href="https://amrex-codes.github.io/amrex/doxygen/classamrex_1_1MultiFab.html">MultiFab</a>
- * \c mf_out.
+ * \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] dcomp offset for the fourth component of the coarsened Array4
+ * object, extracted from its MultiFab, where the interpolated
+ * values will be stored
+ * \param[in] scomp offset for the fourth component of the fine Array4
+ * object, extracted from its MultiFab, containing the
+ * 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
+ * \param[in] crse_ratio coarsening ratio between the fine MultiFab \c mf_src
+ * and the coarsened MultiFab \c mf_dst
+ */
+ void CoarsenAndInterpolateLoop ( MultiFab& mf_dst,
+ const MultiFab& mf_src,
+ const int dcomp,
+ const int scomp,
+ const int ncomp,
+ const int ngrow,
+ const IntVect crse_ratio=IntVect(1) );
+
+ /**
+ * \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_out <a href="https://amrex-codes.github.io/amrex/doxygen/classamrex_1_1MultiFab.html">MultiFab</a>
- * containing the floating point data to be filled with cell-centered averages
- * \param[in] mf_in <a href="https://amrex-codes.github.io/amrex/doxygen/classamrex_1_1MultiFab.html">MultiFab</a>
- * containing the floating point data to be averaged
- * \param[in] dcomp offset for the fourth component of the
- * <a href="https://amrex-codes.github.io/amrex/doxygen/structamrex_1_1Array4.html">Array4</a>
- * object, extracted from its
- * <a href="https://amrex-codes.github.io/amrex/doxygen/classamrex_1_1MultiFab.html">MultiFab</a>
- * , where the cell-centered averages will be stored
- * \param[in] scomp optional offset for the fourth component of the
- * <a href="https://amrex-codes.github.io/amrex/doxygen/structamrex_1_1Array4.html">Array4</a>
- * object, extracted from its
- * <a href="https://amrex-codes.github.io/amrex/doxygen/classamrex_1_1MultiFab.html">MultiFab</a>
- * , containing the data to be averaged
- * \param[in] ncomp optional number of components to loop over for both the
- * input and the output
- * <a href="https://amrex-codes.github.io/amrex/doxygen/structamrex_1_1Array4.html">Array4</a>
- * objects, extracted from their respective
- * <a href="https://amrex-codes.github.io/amrex/doxygen/classamrex_1_1MultiFab.html">MultiFab</a>
+ * \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] dcomp offset for the fourth component of the coarsened Array4
+ * object, extracted from its MultiFab, where the interpolated
+ * values will be stored
+ * \param[in] scomp offset for the fourth component of the fine Array4
+ * object, extracted from its MultiFab, containing the
+ * 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
+ * \param[in] crse_ratio coarsening ratio between the fine MultiFab \c mf_src
+ * and the coarsened MultiFab \c mf_dst
*/
- void ToCellCenter ( MultiFab& mf_out,
- const MultiFab& mf_in,
- const int dcomp,
- const int ngrow,
- const int scomp=0,
- const int ncomp=1 );
+ void CoarsenAndInterpolate ( MultiFab& mf_dst,
+ const MultiFab& mf_src,
+ const int dcomp,
+ const int scomp,
+ const int ncomp,
+ const int ngrow,
+ const IntVect crse_ratio=IntVect(1) );
}
#endif // WARPX_AVERAGE_H_
diff --git a/Source/Utils/Average.cpp b/Source/Utils/Average.cpp
index b98f1e948..ba03bf5b8 100644
--- a/Source/Utils/Average.cpp
+++ b/Source/Utils/Average.cpp
@@ -3,27 +3,112 @@
using namespace amrex;
void
-Average::ToCellCenter ( MultiFab& mf_out,
- const MultiFab& mf_in,
- const int dcomp,
- const int ngrow,
- const int scomp,
- const int ncomp )
+Average::CoarsenAndInterpolateLoop ( MultiFab& mf_dst,
+ const MultiFab& mf_src,
+ const int dcomp,
+ const int scomp,
+ const int ncomp,
+ const int 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();
+
+ if ( crse_ratio > IntVect(1) ) AMREX_ALWAYS_ASSERT_WITH_MESSAGE( ngrow == 0,
+ "option of filling guard cells of destination MultiFab with coarsening not supported for this interpolation" );
+
+ AMREX_ALWAYS_ASSERT_WITH_MESSAGE( mf_src.nGrowVect() >= stag_dst-stag_src+IntVect(ngrow),
+ "source fine MultiFab does not have enough guard cells for this interpolation" );
+
+ // Auxiliary integer arrays (always 3D)
+ Gpu::ManagedVector<int> sf_gpuarr, sc_gpuarr, cr_gpuarr, np_gpuarr;
+ sf_gpuarr.resize( 3 ); // staggering of source fine MultiFab
+ sc_gpuarr.resize( 3 ); // staggering of destination coarse MultiFab
+ cr_gpuarr.resize( 3 ); // coarsening ratio
+ np_gpuarr.resize( 3 ); // number of points to loop over for interpolation
+
+ sf_gpuarr[0] = stag_src[0];
+ sf_gpuarr[1] = stag_src[1];
+#if (AMREX_SPACEDIM == 2)
+ sf_gpuarr[2] = 0;
+#elif (AMREX_SPACEDIM == 3)
+ sf_gpuarr[2] = stag_src[2];
+#endif
+
+ sc_gpuarr[0] = stag_dst[0];
+ sc_gpuarr[1] = stag_dst[1];
+#if (AMREX_SPACEDIM == 2)
+ sc_gpuarr[2] = 0;
+#elif (AMREX_SPACEDIM == 3)
+ sc_gpuarr[2] = stag_dst[2];
+#endif
+
+ cr_gpuarr[0] = crse_ratio[0];
+ cr_gpuarr[1] = crse_ratio[1];
+#if (AMREX_SPACEDIM == 2)
+ cr_gpuarr[2] = 1;
+#elif (AMREX_SPACEDIM == 3)
+ cr_gpuarr[2] = crse_ratio[2];
+#endif
+
+ int const* const AMREX_RESTRICT sf = sf_gpuarr.data();
+ int const* const AMREX_RESTRICT sc = sc_gpuarr.data();
+ int const* const AMREX_RESTRICT cr = cr_gpuarr.data();
+ int * const AMREX_RESTRICT np = np_gpuarr.data();
+
+ // Compute number of points to loop over (either 1 or 2) in each direction
+ for ( int l = 0; l < 3; ++l ) {
+ if ( cr[l] == 1 ) np[l] = 1+abs(sf[l]-sc[l]); // no coarsening
+ else np[l] = 2-sf[l];
+ }
+
#ifdef _OPENMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
// Loop over boxes (or tiles if not on GPU)
- for (MFIter mfi( mf_out, TilingIfNotGPU() ); mfi.isValid(); ++mfi)
+ for (MFIter mfi( mf_dst, TilingIfNotGPU() ); mfi.isValid(); ++mfi)
{
- const Box bx = mfi.growntilebox( ngrow );
- Array4<Real> const& mf_out_arr = mf_out.array( mfi );
- Array4<Real const> const& mf_in_arr = mf_in.const_array( mfi );
- const IntVect stag = mf_in.boxArray().ixType().toIntVect();
+ // 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 )
{
- mf_out_arr(i,j,k,n+dcomp) = Average::ToCellCenter( mf_in_arr, stag, i, j, k, n+scomp );
+ arr_dst(i,j,k,n+dcomp) = Average::CoarsenAndInterpolateKernel(
+ arr_src, sf, sc, cr, np, i, j, k, n+scomp );
} );
}
}
+
+void
+Average::CoarsenAndInterpolate ( MultiFab& mf_dst,
+ const MultiFab& mf_src,
+ const int dcomp,
+ const int scomp,
+ const int ncomp,
+ const int ngrow,
+ const IntVect crse_ratio )
+{
+ BL_PROFILE( "Average::CoarsenAndInterpolate" );
+
+ // Convert BoxArray of source MultiFab to staggering of destination MultiFab and coarsen it (if possible)
+ BoxArray ba_tmp = amrex::convert( mf_src.boxArray(), mf_dst.ixType().toIntVect() );
+ AMREX_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() )
+ Average::CoarsenAndInterpolateLoop( mf_dst, mf_src, dcomp, scomp, ncomp, ngrow, 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, 0, MFInfo(), FArrayBoxFactory() );
+ // 2) interpolate from mf_src to mf_tmp (start writing into component 0)
+ Average::CoarsenAndInterpolateLoop( mf_tmp, mf_src, 0, scomp, ncomp, ngrow, crse_ratio );
+ // 3) copy from mf_tmp to mf_dst (with different BoxArray or DistributionMapping)
+ mf_dst.copy( mf_tmp, 0, dcomp, ncomp );
+ }
+}
generated by cgit v1.2.3 (git 2.46.0) at 2025-08-19 06:47:17 +0000