1#ifndef AMREX_FILLPATCHER_H_
2#define AMREX_FILLPATCHER_H_
3#include <AMReX_Config.H>
76template <
class MF = MultiFab>
129 template <
typename BC,
135 int scomp,
int dcomp,
int ncomp,
136 BC& cbc,
int cbccomp, BC& fbc,
int fbccomp,
138 PreInterpHook
const& pre_interp = {},
139 PostInterpHook
const& post_interp = {});
160 template <
typename BC,
161 typename PreInterpHook=NullInterpHook<MF>,
162 typename PostInterpHook=NullInterpHook<MF> >
166 int scomp,
int dcomp,
int ncomp,
167 BC& cbc,
int cbccomp,
169 PreInterpHook
const& pre_interp = {},
170 PostInterpHook
const& post_interp = {});
181 template <std::
size_t order>
182 requires (order == 3 || order == 4)
202 template <
typename BC>
203 void fillRK (
int stage,
int iteration,
int ncycle, MF& mf,
Real time,
221 std::unique_ptr<MF> m_cf_crse_data_tmp;
222 std::unique_ptr<MF> m_cf_fine_data;
223 Real m_dt_coarse = std::numeric_limits<Real>::lowest();
244 m_eb_index_space(eb_index_space),
245 m_sfine(fba, fdm, 1, nghost,
MFInfo().SetAlloc(false))
248 "FillPatcher<MF>: MF must be FabArray type");
251 for (
int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
258template <
typename BC,
typename PreInterpHook,
typename PostInterpHook>
263 int scomp,
int dcomp,
int ncomp,
264 BC& cbc,
int cbccomp,
265 BC& fbc,
int fbccomp,
267 PreInterpHook
const& pre_interp,
268 PostInterpHook
const& post_interp)
273 !fmf.empty() && fmf.
size() == ft.
size() &&
274 m_fba == fmf[0]->boxArray() && m_fdm == fmf[0]->DistributionMap());
276 fillCoarseFineBoundary(mf, nghost, time, cmf, ct, scomp, dcomp, ncomp,
277 cbc, cbccomp, bcs, bcscomp, pre_interp, post_interp);
280 m_fgeom, fbc, fbccomp);
289 m_fgeom, m_cgeom, m_eb_index_space);
293template <
typename BC,
typename PreInterpHook,
typename PostInterpHook>
298 int scomp,
int dcomp,
int ncomp,
299 BC& cbc,
int cbccomp,
301 PreInterpHook
const& pre_interp,
302 PostInterpHook
const& post_interp)
307 nghost.
allLE(m_nghost) &&
308 m_fba == mf.boxArray() &&
309 m_fdm == mf.DistributionMap() &&
310 m_cba == cmf[0]->boxArray() &&
311 m_cdm == cmf[0]->DistributionMap() &&
313 m_ncomp == cmf[0]->nComp());
315 auto const& fpc = getFPinfo();
317 if ( ! fpc.ba_crse_patch.empty())
319 if (m_cf_fine_data ==
nullptr) {
320 m_cf_fine_data = std::make_unique<MF>
321 (detail::make_mf_fine_patch<MF>(fpc, m_ncomp));
324 int ncmfs = cmf.
size();
325 for (
int icmf = 0; icmf < ncmfs; ++icmf) {
327 auto it = std::find_if(m_cf_crse_data.begin(), m_cf_crse_data.end(),
328 [=] (
auto const&
x) {
329 return amrex::almostEqual(x.first,t,5);
332 if (it == std::end(m_cf_crse_data)) {
333 MF mf_crse_patch = detail::make_mf_crse_patch<MF>(fpc, m_ncomp);
334 mf_crse_patch.ParallelCopy(*cmf[icmf], m_cgeom.periodicity());
336 std::pair<Real,std::unique_ptr<MF>> tmp;
338 tmp.second = std::make_unique<MF>(std::move(mf_crse_patch));
339 m_cf_crse_data.push_back(std::move(tmp));
343 if (m_cf_crse_data_tmp ==
nullptr) {
344 m_cf_crse_data_tmp = std::make_unique<MF>
345 (detail::make_mf_crse_patch<MF>(fpc, m_ncomp));
348 int const ng_space_interp = 8;
349 Box domain = m_cgeom.growPeriodicDomain(ng_space_interp);
353 if (m_cf_crse_data.size() == 1) {
355 }
else if (m_cf_crse_data.size() == 2) {
356 Real const teps = std::abs(m_cf_crse_data[1].first -
357 m_cf_crse_data[0].first) * 1.e-3_rt;
358 if (time > m_cf_crse_data[0].first - teps &&
359 time < m_cf_crse_data[0].first + teps) {
361 }
else if (time > m_cf_crse_data[1].first - teps &&
362 time < m_cf_crse_data[1].first + teps) {
369 if (idata == 0 || idata == 1) {
370 auto const& dst = m_cf_crse_data_tmp->arrays();
371 auto const& src = m_cf_crse_data[idata].second->const_arrays();
376 dst[bi](i,j,k,n) = src[bi](i,j,k,n+scomp);
379 }
else if (idata == 2) {
380 Real t0 = m_cf_crse_data[0].first;
381 Real t1 = m_cf_crse_data[1].first;
382 Real alpha = (t1-time)/(t1-t0);
383 Real beta = (time-t0)/(t1-t0);
384 auto const& a = m_cf_crse_data_tmp->arrays();
385 auto const& a0 = m_cf_crse_data[0].second->const_arrays();
386 auto const& a1 = m_cf_crse_data[1].second->const_arrays();
392 = alpha*a0[bi](i,j,k,scomp+n)
393 + beta*a1[bi](i,j,k,scomp+n);
399 amrex::Abort(
"FillPatcher: High order interpolation in time not supported. Or FillPatcher was not properly deleted.");
404 cbc(*m_cf_crse_data_tmp, 0, ncomp, m_cf_crse_data_tmp->nGrowVect(), time, cbccomp);
406 detail::call_interp_hook(pre_interp, *m_cf_crse_data_tmp, 0, ncomp);
409 ncomp,
IntVect(0), m_cgeom, m_fgeom,
411 mf.ixType()),nghost),
412 m_ratio, m_interp, bcs, bcscomp);
414 detail::call_interp_hook(post_interp, *m_cf_fine_data, scomp, ncomp);
416 mf.ParallelCopy(*m_cf_fine_data, scomp, dcomp, ncomp,
IntVect{0}, nghost);
420template <
typename MF>
421template <std::
size_t order>
422requires (order == 3 || order == 4)
426 BL_PROFILE(
"FillPatcher::storeRKCoarseData()");
428 m_cf_crse_data.resize(order+1);
430 auto const& fpc = getFPinfo();
432 for (
auto& tmf : m_cf_crse_data) {
433 tmf.first = std::numeric_limits<Real>::lowest();
434 tmf.second = std::make_unique<MF>(detail::make_mf_crse_patch<MF>(fpc, m_ncomp));
436 m_cf_crse_data[0].second->ParallelCopy(S_old, m_cgeom.periodicity());
437 for (std::size_t i = 0; i < order; ++i) {
438 m_cf_crse_data[i+1].second->ParallelCopy(RK_k[i], m_cgeom.periodicity());
442template <
typename MF>
443template <
typename BC>
445 MF& mf,
Real time, BC& cbc, BC& fbc,
449 int rk_order = m_cf_crse_data.size()-1;
450 if (rk_order != 3 && rk_order != 4) {
451 amrex::Abort(
"FillPatcher: unsupported RK order "+std::to_string(rk_order));
456 auto const& fpc = getFPinfo();
457 if (m_cf_crse_data_tmp ==
nullptr) {
458 m_cf_crse_data_tmp = std::make_unique<MF>
459 (detail::make_mf_crse_patch<MF>(fpc, m_ncomp));
462 auto const& u = m_cf_crse_data_tmp->arrays();
463 auto const& u0 = m_cf_crse_data[0].second->const_arrays();
464 auto const& k1 = m_cf_crse_data[1].second->const_arrays();
465 auto const& k2 = m_cf_crse_data[2].second->const_arrays();
466 auto const& k3 = m_cf_crse_data[3].second->const_arrays();
468 Real dtc = m_dt_coarse;
472 int const ng_space_interp = 8;
473 Box cdomain = m_cgeom.growPeriodicDomain(ng_space_interp);
474 cdomain.
convert(m_cf_crse_data_tmp->ixType());
478 Real b1 = xsi -
Real(5./6.)*xsi*xsi;
494 Real kk1 = k1[bi](i,j,k,n);
495 Real kk2 = k2[bi](i,j,k,n);
496 Real kk3 = k3[bi](i,j,k,n);
497 Real uu = b1*kk1 + b2*kk2 + b3*kk3;
498 u[bi](i,j,k,n) = u0[bi](i,j,k,n) + dtc*uu;
501 }
else if (stage == 2) {
506 Real kk1 = k1[bi](i,j,k,n);
507 Real kk2 = k2[bi](i,j,k,n);
508 Real kk3 = k3[bi](i,j,k,n);
509 Real uu = b1*kk1 + b2*kk2 + b3*kk3;
510 Real ut = c1*kk1 + c2*kk2 + c3*kk3;
511 u[bi](i,j,k,n) = u0[bi](i,j,k,n) + dtc*(uu + r*ut);
514 }
else if (stage == 3) {
519 Real kk1 = k1[bi](i,j,k,n);
520 Real kk2 = k2[bi](i,j,k,n);
521 Real kk3 = k3[bi](i,j,k,n);
522 Real uu = b1*kk1 + b2*kk2 + b3*kk3;
523 Real ut = c1*kk1 + c2*kk2 + c3*kk3;
524 Real utt = d1*kk1 + d2*kk2 + d3*kk3;
525 u[bi](i,j,k,n) = u0[bi](i,j,k,n) + dtc*
526 (uu + Real(0.5)*r*ut + Real(0.25)*r*r*utt);
530 }
else if (rk_order == 4) {
531 auto const& k4 = m_cf_crse_data[4].second->const_arrays();
533 Real xsi3 = xsi2*xsi;
559 Real kk1 = k1[bi](i,j,k,n);
560 Real kk2 = k2[bi](i,j,k,n);
561 Real kk3 = k3[bi](i,j,k,n);
562 Real kk4 = k4[bi](i,j,k,n);
563 Real uu = b1*kk1 + b2*kk2 + b3*kk3 + b4*kk4;
564 u[bi](i,j,k,n) = u0[bi](i,j,k,n) + dtc*uu;
567 }
else if (stage == 2) {
572 Real kk1 = k1[bi](i,j,k,n);
573 Real kk2 = k2[bi](i,j,k,n);
574 Real kk3 = k3[bi](i,j,k,n);
575 Real kk4 = k4[bi](i,j,k,n);
576 Real uu = b1*kk1 + b2*kk2 + b3*kk3 + b4*kk4;
577 Real ut = c1*kk1 + c2*kk2 + c3*kk3 + c4*kk4;
578 u[bi](i,j,k,n) = u0[bi](i,j,k,n) + dtc*(uu + Real(0.5)*r*ut);
581 }
else if (stage == 3 || stage == 4) {
584 Real at = (stage == 3) ?
Real(0.5)*r : r;
586 Real attt = (stage == 3) ?
Real(0.0625)*r3 :
Real(0.125)*r3;
592 Real kk1 = k1[bi](i,j,k,n);
593 Real kk2 = k2[bi](i,j,k,n);
594 Real kk3 = k3[bi](i,j,k,n);
595 Real kk4 = k4[bi](i,j,k,n);
596 Real uu = b1*kk1 + b2*kk2 + b3*kk3 + b4*kk4;
597 Real ut = c1*kk1 + c2*kk2 + c3*kk3 + c4*kk4;
598 Real utt = d1*kk1 + d2*kk2 + d3*kk3 + d4*kk4;
599 Real uttt = e1*kk1 + e2*kk2 + e3*kk3 + e4*kk4;
600 u[bi](i,j,k,n) = u0[bi](i,j,k,n) + dtc *
601 (uu + at*ut + att*utt + attt*(uttt+akk*(kk3-kk2)));
609 cbc(*m_cf_crse_data_tmp, 0, m_ncomp, m_cf_crse_data_tmp->nGrowVect(), time, 0);
611 if (m_cf_fine_data ==
nullptr) {
612 m_cf_fine_data = std::make_unique<MF>(detail::make_mf_fine_patch<MF>(fpc, m_ncomp));
616 m_ncomp,
IntVect(0), m_cgeom, m_fgeom,
618 mf.ixType()),m_nghost),
619 m_ratio, m_interp, bcs, 0);
622 mf.ParallelCopy(*m_cf_fine_data, 0, 0, m_ncomp,
IntVect(0), m_nghost);
624 mf.FillBoundary(m_fgeom.periodicity());
625 fbc(mf, 0, m_ncomp, m_nghost, time, 0);
#define BL_PROFILE(a)
Definition AMReX_BLProfiler.H:551
#define AMREX_ASSERT(EX)
Definition AMReX_BLassert.H:38
#define AMREX_ALWAYS_ASSERT(EX)
Definition AMReX_BLassert.H:50
High-level FillPatch helpers for AMR coarse-to-fine synchronization.
#define AMREX_GPU_DEVICE
Definition AMReX_GpuQualifiers.H:18
Reference-counted collection of Boxes.
Definition AMReX_BoxArray.H:676
IndexType ixType() const noexcept
Return index type of this BoxArray.
Definition AMReX_BoxArray.H:1252
__host__ __device__ IntVectND< dim > length() const noexcept
Return the length of the BoxND.
Definition AMReX_Box.H:167
__host__ __device__ BoxND & convert(IndexTypeND< dim > typ) noexcept
Convert the BoxND from the current type into the argument type. This may change the BoxND coordinates...
Definition AMReX_Box.H:1008
__host__ __device__ bool contains(const IntVectND< dim > &p) const noexcept
Return true if argument is contained within BoxND.
Definition AMReX_Box.H:233
Calculates the distribution of FABs to MPI processes.
Definition AMReX_DistributionMapping.H:51
Definition AMReX_EB2.H:28
static const FPinfo & TheFPinfo(const FabArrayBase &srcfa, const FabArrayBase &dstfa, const IntVect &dstng, const BoxConverter &coarsener, const Geometry &fgeom, const Geometry &cgeom, const EB2::IndexSpace *)
Definition AMReX_FabArrayBase.cpp:2066
FillPatcher is for filling a fine level MultiFab/FabArray.
Definition AMReX_FillPatcher.H:78
void storeRKCoarseData(Real time, Real dt, MF const &S_old, Array< MF, order > const &RK_k)
Store coarse AMR level data for RK3 and RK4.
Definition AMReX_FillPatcher.H:423
void fillCoarseFineBoundary(MF &mf, IntVect const &nghost, Real time, Vector< MF * > const &cmf, Vector< Real > const &ct, int scomp, int dcomp, int ncomp, BC &cbc, int cbccomp, Vector< BCRec > const &bcs, int bcscomp, PreInterpHook const &pre_interp={}, PostInterpHook const &post_interp={})
Function to fill data at coarse/fine boundary only.
Definition AMReX_FillPatcher.H:295
void fill(MF &mf, IntVect const &nghost, Real time, Vector< MF * > const &cmf, Vector< Real > const &ct, Vector< MF * > const &fmf, Vector< Real > const &ft, int scomp, int dcomp, int ncomp, BC &cbc, int cbccomp, BC &fbc, int fbccomp, Vector< BCRec > const &bcs, int bcscomp, PreInterpHook const &pre_interp={}, PostInterpHook const &post_interp={})
Function to fill data.
Definition AMReX_FillPatcher.H:260
void fillRK(int stage, int iteration, int ncycle, MF &mf, Real time, BC &cbc, BC &fbc, Vector< BCRec > const &bcs)
Fill ghost cells of fine AMR level for RK3 and RK4.
Definition AMReX_FillPatcher.H:444
FillPatcher(BoxArray const &fba, DistributionMapping const &fdm, Geometry const &fgeom, BoxArray const &cba, DistributionMapping const &cdm, Geometry const &cgeom, IntVect const &nghost, int ncomp, InterpBase *interp, EB2::IndexSpace const *eb_index_space=EB2::TopIndexSpaceIfPresent())
Constructor of FillPatcher.
Definition AMReX_FillPatcher.H:229
Rectangular problem domain geometry.
Definition AMReX_Geometry.H:75
const Box & Domain() const noexcept
Returns our rectangular domain.
Definition AMReX_Geometry.H:216
__host__ __device__ bool cellCentered() const noexcept
True if the IndexTypeND is CELL based in all directions.
Definition AMReX_IndexType.H:102
__host__ __device__ bool nodeCentered() const noexcept
True if the IndexTypeND is NODE based in all directions.
Definition AMReX_IndexType.H:108
__host__ __device__ constexpr bool allLE(const IntVectND< dim > &rhs) const noexcept
Returns true if this is less than or equal to argument for all components. NOTE: This is NOT a strict...
Definition AMReX_IntVect.H:492
Definition AMReX_InterpBase.H:34
Definition AMReX_InterpBase.H:20
This class is a thin wrapper around std::vector. Unlike vector, Vector::operator[] provides bound che...
Definition AMReX_Vector.H:29
Long size() const noexcept
Definition AMReX_Vector.H:54
amrex_real Real
Floating Point Type for Fields.
Definition AMReX_REAL.H:79
__host__ __device__ BoxND< dim > convert(const BoxND< dim > &b, const IntVectND< dim > &typ) noexcept
Return a copy of b converted to the nodal flags typ.
Definition AMReX_Box.H:1630
__host__ __device__ BoxND< dim > grow(const BoxND< dim > &b, int i) noexcept
Return a copy of b grown uniformly by i cells in every direction.
Definition AMReX_Box.H:1326
__host__ __device__ BoxND< dim > refine(const BoxND< dim > &b, int ref_ratio) noexcept
Return a copy of b refined by the isotropic ratio ref_ratio.
Definition AMReX_Box.H:1510
std::array< T, N > Array
Definition AMReX_Array.H:31
const IndexSpace * TopIndexSpaceIfPresent() noexcept
Return the top IndexSpace if one has been built (nullptr otherwise).
Definition AMReX_EB2.cpp:93
void streamSynchronize() noexcept
Definition AMReX_GpuDevice.H:310
Definition AMReX_Amr.cpp:50
void ParallelFor(TypeList< CTOs... > ctos, std::array< int, sizeof...(CTOs)> const &runtime_options, T N, F &&f)
Definition AMReX_CTOParallelForImpl.H:202
void FillPatchSingleLevel(MF &mf, IntVect const &nghost, Real time, const Vector< MF * > &smf, const Vector< Real > &stime, int scomp, int dcomp, int ncomp, const Geometry &geom, BC &physbcf, int bcfcomp)
FillPatch with data from the current level.
Definition AMReX_FillPatchUtil_I.H:75
IntVectND< 3 > IntVect
IntVect is an alias for amrex::IntVectND instantiated with AMREX_SPACEDIM.
Definition AMReX_BaseFwd.H:38
void Abort(const std::string &msg)
Print a fatal-error message to stderr and abort execution.
Definition AMReX.cpp:241
void FillPatchInterp(MultiFab &mf_fine_patch, int fcomp, MultiFab const &mf_crse_patch, int ccomp, int ncomp, IntVect const &ng, const Geometry &cgeom, const Geometry &fgeom, Box const &dest_domain, const IntVect &ratio, MFInterpolater *mapper, const Vector< BCRec > &bcs, int bcscomp)
Helper that applies a MFInterpolater to fill a fine patch from a coarse patch.
Definition AMReX_FillPatchUtil.cpp:139
Definition AMReX_FabArrayBase.H:304
Definition AMReX_TypeTraits.H:27
FabArray memory allocation information.
Definition AMReX_FabArray.H:68
Definition AMReX_FillPatchUtil.H:39