Line data Source code
1 : #ifndef AMREX_FillPatchUtil_I_H_
2 : #define AMREX_FillPatchUtil_I_H_
3 : #include <AMReX_Config.H>
4 :
5 : namespace amrex {
6 :
7 : namespace detail {
8 :
9 : template <typename F, typename MF>
10 132 : auto call_interp_hook (F const& f, MF& mf, int icomp, int ncomp)
11 : -> decltype(f(mf[0],Box(),icomp,ncomp))
12 : {
13 : #ifdef AMREX_USE_OMP
14 : #pragma omp parallel if (Gpu::notInLaunchRegion())
15 : #endif
16 564 : for (MFIter mfi(mf); mfi.isValid(); ++mfi) {
17 432 : auto& dfab = mf[mfi];
18 432 : const Box& dbx = dfab.box();
19 432 : f(dfab, dbx, icomp, ncomp);
20 : }
21 132 : }
22 :
23 : template <typename F, typename MF>
24 : auto call_interp_hook (F const& f, MF& mf, int icomp, int ncomp)
25 : -> decltype(f(mf,icomp,ncomp))
26 : {
27 : f(mf, icomp, ncomp);
28 : }
29 :
30 : }
31 :
32 : template <typename Interp>
33 : bool ProperlyNested (const IntVect& ratio, const IntVect& blocking_factor, int ngrow,
34 : const IndexType& boxType, Interp* mapper)
35 : {
36 : int ratio_max = ratio[0];
37 : #if (AMREX_SPACEDIM > 1)
38 : ratio_max = std::max(ratio_max, ratio[1]);
39 : #endif
40 : #if (AMREX_SPACEDIM == 3)
41 : ratio_max = std::max(ratio_max, ratio[2]);
42 : #endif
43 : // There are at least this many coarse cells outside fine grids
44 : // (except at physical boundaries).
45 : const IntVect& nbuf = blocking_factor / ratio_max;
46 :
47 : Box crse_box(IntVect(AMREX_D_DECL(0 ,0 ,0 )), IntVect(AMREX_D_DECL(4*nbuf[0]-1,
48 : 4*nbuf[1]-1,
49 : 4*nbuf[2]-1)));
50 : crse_box.convert(boxType);
51 : Box fine_box(nbuf, IntVect(AMREX_D_DECL(3*nbuf[0]-1,3*nbuf[1]-1,3*nbuf[2]-1)));
52 : fine_box.convert(boxType);
53 : fine_box.refine(ratio_max);
54 : fine_box.grow(ngrow);
55 : const Box& fine_box_coarsened = mapper->CoarseBox(fine_box, ratio_max);
56 : return crse_box.contains(fine_box_coarsened);
57 : }
58 :
59 : template <typename MF, typename BC>
60 : std::enable_if_t<IsFabArray<MF>::value>
61 370 : FillPatchSingleLevel (MF& mf, Real time,
62 : const Vector<MF*>& smf, const Vector<Real>& stime,
63 : int scomp, int dcomp, int ncomp,
64 : const Geometry& geom,
65 : BC& physbcf, int bcfcomp)
66 : {
67 370 : FillPatchSingleLevel(mf, mf.nGrowVect(), time, smf, stime, scomp, dcomp, ncomp,
68 : geom, physbcf, bcfcomp);
69 370 : }
70 :
71 : template <typename MF, typename BC>
72 : std::enable_if_t<IsFabArray<MF>::value>
73 436 : FillPatchSingleLevel (MF& mf, IntVect const& nghost, Real time,
74 : const Vector<MF*>& smf, const Vector<Real>& stime,
75 : int scomp, int dcomp, int ncomp,
76 : const Geometry& geom,
77 : BC& physbcf, int bcfcomp)
78 : {
79 : BL_PROFILE("FillPatchSingleLevel");
80 :
81 : AMREX_ASSERT(scomp+ncomp <= smf[0]->nComp());
82 : AMREX_ASSERT(dcomp+ncomp <= mf.nComp());
83 : AMREX_ASSERT(smf.size() == stime.size());
84 : AMREX_ASSERT(!smf.empty());
85 : AMREX_ASSERT(nghost.allLE(mf.nGrowVect()));
86 :
87 436 : if (smf.size() == 1)
88 : {
89 436 : if (&mf == smf[0] && scomp == dcomp) {
90 360 : mf.FillBoundary(dcomp, ncomp, nghost, geom.periodicity());
91 : } else {
92 76 : mf.ParallelCopy(*smf[0], scomp, dcomp, ncomp, IntVect{0}, nghost, geom.periodicity());
93 : }
94 : }
95 0 : else if (smf.size() == 2)
96 : {
97 : BL_ASSERT(smf[0]->boxArray() == smf[1]->boxArray());
98 0 : MF raii;
99 : MF * dmf;
100 : int destcomp;
101 : bool sameba;
102 0 : if (mf.boxArray() == smf[0]->boxArray() &&
103 0 : mf.DistributionMap() == smf[0]->DistributionMap())
104 : {
105 0 : dmf = &mf;
106 0 : destcomp = dcomp;
107 0 : sameba = true;
108 : } else {
109 0 : raii.define(smf[0]->boxArray(), smf[0]->DistributionMap(), ncomp, 0,
110 0 : MFInfo(), smf[0]->Factory());
111 :
112 0 : dmf = &raii;
113 0 : destcomp = 0;
114 0 : sameba = false;
115 : }
116 :
117 0 : if ((dmf != smf[0] && dmf != smf[1]) || scomp != dcomp)
118 : {
119 : #ifdef AMREX_USE_OMP
120 : #pragma omp parallel if (Gpu::notInLaunchRegion())
121 : #endif
122 0 : for (MFIter mfi(*dmf,TilingIfNotGPU()); mfi.isValid(); ++mfi)
123 : {
124 0 : const Box& bx = mfi.tilebox();
125 0 : const Real t0 = stime[0];
126 0 : const Real t1 = stime[1];
127 0 : auto const sfab0 = smf[0]->array(mfi);
128 0 : auto const sfab1 = smf[1]->array(mfi);
129 0 : auto dfab = dmf->array(mfi);
130 :
131 0 : if (time == t0)
132 : {
133 0 : AMREX_HOST_DEVICE_PARALLEL_FOR_4D ( bx, ncomp, i, j, k, n,
134 : {
135 : dfab(i,j,k,n+destcomp) = sfab0(i,j,k,n+scomp);
136 : });
137 : }
138 0 : else if (time == t1)
139 : {
140 0 : AMREX_HOST_DEVICE_PARALLEL_FOR_4D ( bx, ncomp, i, j, k, n,
141 : {
142 : dfab(i,j,k,n+destcomp) = sfab1(i,j,k,n+scomp);
143 : });
144 : }
145 0 : else if (! amrex::almostEqual(t0,t1))
146 : {
147 0 : Real alpha = (t1-time)/(t1-t0);
148 0 : Real beta = (time-t0)/(t1-t0);
149 0 : AMREX_HOST_DEVICE_PARALLEL_FOR_4D ( bx, ncomp, i, j, k, n,
150 : {
151 : dfab(i,j,k,n+destcomp) = alpha*sfab0(i,j,k,n+scomp)
152 : + beta*sfab1(i,j,k,n+scomp);
153 : });
154 : }
155 : else
156 : {
157 0 : AMREX_HOST_DEVICE_PARALLEL_FOR_4D ( bx, ncomp, i, j, k, n,
158 : {
159 : dfab(i,j,k,n+destcomp) = sfab0(i,j,k,n+scomp);
160 : });
161 : }
162 : }
163 : }
164 :
165 0 : if (sameba)
166 : {
167 : // Note that when sameba is true mf's BoxArray is nonoverlapping.
168 : // So FillBoundary is safe.
169 0 : mf.FillBoundary(dcomp, ncomp, nghost, geom.periodicity());
170 : }
171 : else
172 : {
173 0 : IntVect src_ngrow = IntVect::TheZeroVector();
174 0 : IntVect dst_ngrow = nghost;
175 :
176 0 : mf.ParallelCopy(*dmf, 0, dcomp, ncomp, src_ngrow, dst_ngrow, geom.periodicity());
177 : }
178 : }
179 : else {
180 : amrex::Abort("FillPatchSingleLevel: high-order interpolation in time not implemented yet");
181 : }
182 :
183 436 : physbcf(mf, dcomp, ncomp, nghost, time, bcfcomp);
184 436 : }
185 :
186 : void FillPatchInterp (MultiFab& mf_fine_patch, int fcomp, MultiFab const& mf_crse_patch, int ccomp,
187 : int ncomp, IntVect const& ng, const Geometry& cgeom, const Geometry& fgeom,
188 : Box const& dest_domain, const IntVect& ratio,
189 : MFInterpolater* mapper, const Vector<BCRec>& bcs, int bcscomp);
190 :
191 : template <typename MF, typename Interp>
192 : std::enable_if_t<IsFabArray<MF>::value && !std::is_same_v<Interp,MFInterpolater>>
193 66 : FillPatchInterp (MF& mf_fine_patch, int fcomp, MF const& mf_crse_patch, int ccomp,
194 : int ncomp, IntVect const& ng, const Geometry& cgeom, const Geometry& fgeom,
195 : Box const& dest_domain, const IntVect& ratio,
196 : Interp* mapper, const Vector<BCRec>& bcs, int bcscomp)
197 : {
198 : BL_PROFILE("FillPatchInterp(Fab)");
199 :
200 66 : Box const& cdomain = amrex::convert(cgeom.Domain(), mf_fine_patch.ixType());
201 66 : int idummy=0;
202 : #ifdef AMREX_USE_OMP
203 : #pragma omp parallel if (Gpu::notInLaunchRegion())
204 : #endif
205 : {
206 132 : Vector<BCRec> bcr(ncomp);
207 282 : for (MFIter mfi(mf_fine_patch); mfi.isValid(); ++mfi)
208 : {
209 216 : auto& sfab = mf_crse_patch[mfi];
210 216 : const Box& sbx = sfab.box();
211 :
212 216 : auto& dfab = mf_fine_patch[mfi];
213 432 : Box const& dbx = amrex::grow(mfi.validbox(),ng) & dest_domain;
214 :
215 216 : amrex::setBC(sbx,cdomain,bcscomp,0,ncomp,bcs,bcr);
216 216 : mapper->interp(sfab, ccomp, dfab, fcomp, ncomp, dbx, ratio,
217 : cgeom, fgeom, bcr, idummy, idummy, RunOn::Gpu);
218 : }
219 : }
220 66 : }
221 :
222 : template <typename MF>
223 : std::enable_if_t<IsFabArray<MF>::value>
224 : FillPatchInterp (MF& mf_fine_patch, int fcomp, MF const& mf_crse_patch, int ccomp,
225 : int ncomp, IntVect const& ng, const Geometry& cgeom, const Geometry& fgeom,
226 : Box const& dest_domain, const IntVect& ratio,
227 : InterpBase* mapper, const Vector<BCRec>& bcs, int bcscomp)
228 : {
229 : if (dynamic_cast<MFInterpolater*>(mapper)) {
230 : FillPatchInterp(mf_fine_patch, fcomp, mf_crse_patch, ccomp,
231 : ncomp, ng, cgeom, fgeom, dest_domain, ratio,
232 : static_cast<MFInterpolater*>(mapper), bcs, bcscomp);
233 : } else if (dynamic_cast<Interpolater*>(mapper)) {
234 : FillPatchInterp(mf_fine_patch, fcomp, mf_crse_patch, ccomp,
235 : ncomp, ng, cgeom, fgeom, dest_domain, ratio,
236 : static_cast<Interpolater*>(mapper), bcs, bcscomp);
237 : } else {
238 : amrex::Abort("FillPatchInterp: unknown InterpBase");
239 : }
240 : }
241 :
242 : template <typename MF, typename iMF, typename Interp>
243 : std::enable_if_t<IsFabArray<MF>::value && !std::is_same_v<Interp,MFInterpolater>>
244 0 : InterpFace (Interp *interp,
245 : MF const& mf_crse_patch, int crse_comp,
246 : MF& mf_refined_patch, int fine_comp,
247 : int ncomp, const IntVect& ratio,
248 : const iMF& solve_mask, const Geometry& crse_geom, const Geometry& fine_geom,
249 : int bcscomp, RunOn gpu_or_cpu,
250 : const Vector<BCRec>& bcs)
251 : {
252 0 : Vector<BCRec> bcr(ncomp);
253 0 : Box const& cdomain = amrex::convert(crse_geom.Domain(), mf_crse_patch.ixType());
254 0 : for (MFIter mfi(mf_refined_patch);mfi.isValid(); ++mfi)
255 : {
256 0 : auto& sfab = mf_crse_patch[mfi];
257 0 : const Box& sbx = sfab.box();
258 0 : auto& dfab = mf_refined_patch[mfi];
259 0 : Box const& dbx = dfab.box();
260 0 : auto& ifab = solve_mask[mfi];
261 0 : amrex::setBC(sbx,cdomain,bcscomp,0,ncomp,bcs,bcr);
262 0 : interp->interp_face(sfab,crse_comp,dfab,fine_comp,ncomp,
263 : dbx, ratio, ifab, crse_geom, fine_geom,
264 : bcr, bcscomp, gpu_or_cpu);
265 : }
266 0 : }
267 :
268 : template <typename MF, typename iMF>
269 : std::enable_if_t<IsFabArray<MF>::value>
270 : InterpFace (InterpBase *interp,
271 : MF const& mf_crse_patch, int crse_comp,
272 : MF& mf_refined_patch, int fine_comp,
273 : int ncomp, const IntVect& ratio,
274 : const iMF& solve_mask, const Geometry& crse_geom, const Geometry& fine_geom,
275 : int bccomp, RunOn gpu_or_cpu,
276 : const Vector<BCRec>& bcs)
277 : {
278 : if (dynamic_cast<MFInterpolater*>(interp)){
279 : InterpFace(static_cast<MFInterpolater*>(interp),
280 : mf_crse_patch, crse_comp,mf_refined_patch, fine_comp,
281 : ncomp, ratio, solve_mask, crse_geom, fine_geom, bccomp,
282 : gpu_or_cpu, bcs);
283 : }
284 : else if (dynamic_cast<Interpolater*>(interp)){
285 : InterpFace(static_cast<Interpolater*>(interp),
286 : mf_crse_patch, crse_comp,mf_refined_patch, fine_comp,
287 : ncomp, ratio, solve_mask, crse_geom, fine_geom, bccomp,
288 : gpu_or_cpu, bcs);
289 : }
290 : else {
291 : amrex::Abort("InterpFace: unknown InterpBase");
292 : }
293 : }
294 :
295 :
296 : namespace detail {
297 :
298 : // ======== FArrayBox
299 :
300 : template <typename MF,
301 : std::enable_if_t<std::is_same_v<typename MF::FABType::value_type,
302 : FArrayBox>,
303 : int> = 0>
304 66 : MF make_mf_crse_patch (FabArrayBase::FPinfo const& fpc, int ncomp)
305 : {
306 66 : MF mf_crse_patch(fpc.ba_crse_patch, fpc.dm_patch, ncomp, 0, MFInfo(),
307 66 : *fpc.fact_crse_patch);
308 66 : return mf_crse_patch;
309 : }
310 :
311 : template <typename MF,
312 : std::enable_if_t<std::is_same_v<typename MF::FABType::value_type,
313 : FArrayBox>,
314 : int> = 0>
315 0 : MF make_mf_crse_patch (FabArrayBase::FPinfo const& fpc, int ncomp, IndexType idx_type)
316 : {
317 0 : MF mf_crse_patch(amrex::convert(fpc.ba_crse_patch, idx_type), fpc.dm_patch,
318 0 : ncomp, 0, MFInfo(), *fpc.fact_crse_patch);
319 0 : return mf_crse_patch;
320 : }
321 :
322 : template <typename MF,
323 : std::enable_if_t<std::is_same_v<typename MF::FABType::value_type,
324 : FArrayBox>,
325 : int> = 0>
326 66 : MF make_mf_fine_patch (FabArrayBase::FPinfo const& fpc, int ncomp)
327 : {
328 66 : MF mf_fine_patch(fpc.ba_fine_patch, fpc.dm_patch, ncomp, 0, MFInfo(),
329 66 : *fpc.fact_fine_patch);
330 66 : return mf_fine_patch;
331 : }
332 :
333 : template <typename MF,
334 : std::enable_if_t<std::is_same_v<typename MF::FABType::value_type,
335 : FArrayBox>,
336 : int> = 0>
337 : MF make_mf_fine_patch (FabArrayBase::FPinfo const& fpc, int ncomp, IndexType idx_type)
338 : {
339 : MF mf_fine_patch(amrex::convert(fpc.ba_fine_patch, idx_type), fpc.dm_patch,
340 : ncomp, 0, MFInfo(), *fpc.fact_fine_patch);
341 : return mf_fine_patch;
342 : }
343 :
344 : template <typename MF,
345 : std::enable_if_t<std::is_same_v<typename MF::FABType::value_type,
346 : FArrayBox>,
347 : int> = 0>
348 0 : MF make_mf_refined_patch (FabArrayBase::FPinfo const& fpc, int ncomp, IndexType idx_type, IntVect ratio)
349 : {
350 0 : MF mf_refined_patch(amrex::convert( amrex::refine( amrex::coarsen(fpc.ba_fine_patch, ratio), ratio), idx_type),
351 0 : fpc.dm_patch, ncomp, 0, MFInfo(), *fpc.fact_fine_patch);
352 0 : return mf_refined_patch;
353 : }
354 :
355 : template <typename MF,
356 : std::enable_if_t<std::is_same_v<typename MF::FABType::value_type,
357 : FArrayBox>,
358 : int> = 0>
359 : MF make_mf_crse_mask (FabArrayBase::FPinfo const& fpc, int ncomp, IndexType idx_type, IntVect ratio)
360 : {
361 : MF mf_crse_mask(amrex::convert(amrex::coarsen(fpc.ba_fine_patch, ratio), idx_type),
362 : fpc.dm_patch, ncomp, 0, MFInfo(), *fpc.fact_fine_patch);
363 : return mf_crse_mask;
364 : }
365 :
366 : template <typename MF,
367 : std::enable_if_t<std::is_same_v<typename MF::FABType::value_type,
368 : FArrayBox>,
369 : int> = 0>
370 66 : void mf_set_domain_bndry (MF &mf, Geometry const & geom)
371 : {
372 66 : mf.setDomainBndry(std::numeric_limits<Real>::quiet_NaN(), geom);
373 66 : }
374 :
375 :
376 : // ======== Not FArrayBox
377 :
378 : template <typename MF,
379 : std::enable_if_t<!std::is_same_v<typename MF::FABType::value_type,
380 : FArrayBox>,
381 : int> = 0>
382 0 : MF make_mf_crse_patch (FabArrayBase::FPinfo const& fpc, int ncomp)
383 : {
384 0 : return MF(fpc.ba_crse_patch, fpc.dm_patch, ncomp, 0);
385 : }
386 :
387 : template <typename MF,
388 : std::enable_if_t<!std::is_same_v<typename MF::FABType::value_type,
389 : FArrayBox>,
390 : int> = 0>
391 0 : MF make_mf_crse_patch (FabArrayBase::FPinfo const& fpc, int ncomp, IndexType idx_type)
392 : {
393 0 : return MF(amrex::convert(fpc.ba_crse_patch, idx_type), fpc.dm_patch, ncomp, 0);
394 : }
395 :
396 : template <typename MF,
397 : std::enable_if_t<!std::is_same_v<typename MF::FABType::value_type,
398 : FArrayBox>,
399 : int> = 0>
400 0 : MF make_mf_fine_patch (FabArrayBase::FPinfo const& fpc, int ncomp)
401 : {
402 0 : return MF(fpc.ba_fine_patch, fpc.dm_patch, ncomp, 0);
403 : }
404 :
405 : template <typename MF,
406 : std::enable_if_t<!std::is_same_v<typename MF::FABType::value_type,
407 : FArrayBox>,
408 : int> = 0>
409 : MF make_mf_fine_patch (FabArrayBase::FPinfo const& fpc, int ncomp, IndexType idx_type)
410 : {
411 : return MF(amrex::convert(fpc.ba_fine_patch, idx_type), fpc.dm_patch, ncomp, 0);
412 : }
413 :
414 : template <typename MF,
415 : std::enable_if_t<!std::is_same_v<typename MF::FABType::value_type,
416 : FArrayBox>,
417 : int> = 0>
418 0 : MF make_mf_refined_patch (FabArrayBase::FPinfo const& fpc, int ncomp, IndexType idx_type, IntVect ratio)
419 : {
420 0 : return MF(amrex::convert( amrex::refine( amrex::coarsen(fpc.ba_fine_patch, ratio), ratio), idx_type), fpc.dm_patch, ncomp, 0);
421 : }
422 :
423 : template <typename MF,
424 : std::enable_if_t<!std::is_same_v<typename MF::FABType::value_type,
425 : FArrayBox>,
426 : int> = 0>
427 0 : MF make_mf_crse_mask (FabArrayBase::FPinfo const& fpc, int ncomp, IndexType idx_type, IntVect ratio)
428 : {
429 0 : return MF(amrex::convert(amrex::coarsen(fpc.ba_fine_patch, ratio), idx_type), fpc.dm_patch, ncomp, 0);
430 : }
431 :
432 : template <typename MF,
433 : std::enable_if_t<!std::is_same_v<typename MF::FABType::value_type,
434 : FArrayBox>,
435 : int> = 0>
436 0 : void mf_set_domain_bndry (MF &/*mf*/, Geometry const & /*geom*/)
437 : {
438 : // nothing
439 0 : }
440 :
441 : template <typename MF, typename BC, typename Interp, typename PreInterpHook, typename PostInterpHook>
442 : std::enable_if_t<IsFabArray<MF>::value,int>
443 66 : FillPatchTwoLevels_doit (MF& mf, IntVect const& nghost, Real time,
444 : const Vector<MF*>& cmf, const Vector<Real>& ct,
445 : const Vector<MF*>& fmf, const Vector<Real>& ft,
446 : int scomp, int dcomp, int ncomp,
447 : const Geometry& cgeom, const Geometry& fgeom,
448 : BC& cbc, int cbccomp,
449 : BC& fbc, int fbccomp,
450 : const IntVect& ratio,
451 : Interp* mapper,
452 : const Vector<BCRec>& bcs, int bcscomp,
453 : const PreInterpHook& pre_interp,
454 : const PostInterpHook& post_interp,
455 : EB2::IndexSpace const* index_space,
456 : bool return_error_code = false)
457 : {
458 : BL_PROFILE("FillPatchTwoLevels");
459 :
460 66 : int success_code = return_error_code ? 0 : -1;
461 66 : int failure_code = 1;
462 :
463 66 : if (nghost.max() > 0 || mf.getBDKey() != fmf[0]->getBDKey())
464 : {
465 66 : const InterpolaterBoxCoarsener& coarsener = mapper->BoxCoarsener(ratio);
466 :
467 : // Test for Face-centered data
468 198 : if ( AMREX_D_TERM( mf.ixType().nodeCentered(0),
469 : + mf.ixType().nodeCentered(1),
470 66 : + mf.ixType().nodeCentered(2) ) == 1 )
471 : {
472 0 : if ( !dynamic_cast<Interpolater*>(mapper) ){
473 : amrex::Abort("This interpolater has not yet implemented a version for face-based data");
474 : }
475 :
476 : // Convert to cell-centered MF meta-data for FPInfo.
477 0 : MF mf_cc_dummy( amrex::convert(mf.boxArray(), IntVect::TheZeroVector()),
478 0 : mf.DistributionMap(), ncomp, nghost, MFInfo().SetAlloc(false) );
479 0 : MF fmf_cc_dummy( amrex::convert(fmf[0]->boxArray(), IntVect::TheZeroVector()),
480 0 : fmf[0]->DistributionMap(), ncomp, nghost, MFInfo().SetAlloc(false) );
481 :
482 0 : const FabArrayBase::FPinfo& fpc = FabArrayBase::TheFPinfo(fmf_cc_dummy, mf_cc_dummy,
483 : nghost,
484 : coarsener,
485 : fgeom,
486 : cgeom,
487 : index_space);
488 :
489 0 : if ( ! fpc.ba_crse_patch.empty())
490 : {
491 0 : if (return_error_code) {
492 0 : BoxArray const& cba = amrex::convert(cmf[0]->boxArray(), IntVect(0));
493 0 : if (!cba.contains(fpc.ba_crse_patch,cgeom.periodicity())) {
494 0 : return failure_code;
495 : }
496 : }
497 :
498 0 : MF mf_crse_patch = make_mf_crse_patch<MF> (fpc, ncomp, mf.boxArray().ixType());
499 : // Must make sure fine exists under needed coarse faces.
500 : // It stores values for the final (interior) interpolation,
501 : // which is done from this fine MF that's been partially filled
502 : // (with only faces overlying coarse having valid data).
503 0 : MF mf_refined_patch = make_mf_refined_patch<MF> (fpc, ncomp, mf.boxArray().ixType(), ratio);
504 0 : auto solve_mask = make_mf_crse_mask<iMultiFab>(fpc, ncomp, mf.boxArray().ixType(), ratio);
505 :
506 0 : mf_set_domain_bndry(mf_crse_patch, cgeom);
507 0 : FillPatchSingleLevel(mf_crse_patch, time, cmf, ct, scomp, 0, ncomp,
508 : cgeom, cbc, cbccomp);
509 :
510 0 : mf_set_domain_bndry(mf_refined_patch, fgeom);
511 0 : FillPatchSingleLevel(mf_refined_patch, time, fmf, ft, scomp, 0, ncomp,
512 : fgeom, fbc, fbccomp);
513 :
514 : // Aliased MFs, used to allow CPC caching.
515 0 : MF mf_known( amrex::coarsen(fmf[0]->boxArray(), ratio), fmf[0]->DistributionMap(),
516 0 : ncomp, nghost, MFInfo().SetAlloc(false) );
517 0 : MF mf_solution( amrex::coarsen(mf_refined_patch.boxArray(), ratio), mf_refined_patch.DistributionMap(),
518 0 : ncomp, 0, MFInfo().SetAlloc(false) );
519 :
520 0 : const FabArrayBase::CPC mask_cpc( mf_solution, IntVect::TheZeroVector(),
521 0 : mf_known, IntVect::TheZeroVector(),
522 0 : cgeom.periodicity());
523 :
524 0 : solve_mask.setVal(1); // Values to solve.
525 0 : solve_mask.setVal(0, mask_cpc, 0, 1); // Known values.
526 :
527 0 : detail::call_interp_hook(pre_interp, mf_crse_patch, 0, ncomp);
528 :
529 0 : InterpFace(mapper, mf_crse_patch, 0, mf_refined_patch, 0, ncomp,
530 : ratio, solve_mask, cgeom, fgeom, bcscomp, RunOn::Gpu, bcs);
531 :
532 0 : detail::call_interp_hook(post_interp, mf_refined_patch, 0, ncomp);
533 :
534 0 : bool aliasing = false;
535 0 : for (auto const& fmf_a : fmf) {
536 0 : aliasing = aliasing || (&mf == fmf_a);
537 : }
538 0 : if (aliasing) {
539 0 : mf.ParallelCopyToGhost(mf_refined_patch, 0, dcomp, ncomp,
540 : IntVect{0}, nghost);
541 : } else {
542 0 : mf.ParallelCopy(mf_refined_patch, 0, dcomp, ncomp,
543 : IntVect{0}, nghost);
544 : }
545 : }
546 : }
547 : else
548 : {
549 66 : const FabArrayBase::FPinfo& fpc = FabArrayBase::TheFPinfo(*fmf[0], mf,
550 : nghost,
551 : coarsener,
552 : fgeom,
553 : cgeom,
554 : index_space);
555 :
556 66 : if ( ! fpc.ba_crse_patch.empty())
557 : {
558 66 : if (return_error_code) {
559 0 : BoxArray const& cba = cmf[0]->boxArray();
560 0 : if (!cba.contains(fpc.ba_crse_patch,cgeom.periodicity())) {
561 0 : return failure_code;
562 : }
563 : }
564 :
565 132 : MF mf_crse_patch = make_mf_crse_patch<MF>(fpc, ncomp);
566 66 : mf_set_domain_bndry (mf_crse_patch, cgeom);
567 :
568 66 : FillPatchSingleLevel(mf_crse_patch, time, cmf, ct, scomp, 0, ncomp, cgeom, cbc, cbccomp);
569 :
570 66 : MF mf_fine_patch = make_mf_fine_patch<MF>(fpc, ncomp);
571 :
572 66 : detail::call_interp_hook(pre_interp, mf_crse_patch, 0, ncomp);
573 :
574 66 : FillPatchInterp(mf_fine_patch, 0, mf_crse_patch, 0,
575 : ncomp, IntVect(0), cgeom, fgeom,
576 66 : amrex::grow(amrex::convert(fgeom.Domain(),mf.ixType()),nghost),
577 : ratio, mapper, bcs, bcscomp);
578 :
579 66 : detail::call_interp_hook(post_interp, mf_fine_patch, 0, ncomp);
580 :
581 66 : mf.ParallelCopy(mf_fine_patch, 0, dcomp, ncomp, IntVect{0}, nghost);
582 : }
583 : }
584 : }
585 :
586 66 : FillPatchSingleLevel(mf, nghost, time, fmf, ft, scomp, dcomp, ncomp,
587 : fgeom, fbc, fbccomp);
588 :
589 66 : return success_code;
590 : }
591 :
592 : template <typename MF, typename BC, typename Interp, typename PreInterpHook, typename PostInterpHook>
593 : std::enable_if_t<IsFabArray<MF>::value>
594 : FillPatchTwoLevels_doit (Array<MF*, AMREX_SPACEDIM> const& mf, IntVect const& nghost, Real time,
595 : const Vector<Array<MF*, AMREX_SPACEDIM> >& cmf, const Vector<Real>& ct,
596 : const Vector<Array<MF*, AMREX_SPACEDIM> >& fmf, const Vector<Real>& ft,
597 : int scomp, int dcomp, int ncomp,
598 : const Geometry& cgeom, const Geometry& fgeom,
599 : Array<BC, AMREX_SPACEDIM>& cbc, const Array<int, AMREX_SPACEDIM>& cbccomp,
600 : Array<BC, AMREX_SPACEDIM>& fbc, const Array<int, AMREX_SPACEDIM>& fbccomp,
601 : const IntVect& ratio,
602 : Interp* mapper,
603 : const Array<Vector<BCRec>, AMREX_SPACEDIM>& bcs, const Array<int, AMREX_SPACEDIM>& bcscomp,
604 : const PreInterpHook& pre_interp,
605 : const PostInterpHook& post_interp,
606 : EB2::IndexSpace const* index_space)
607 : {
608 : BL_PROFILE("FillPatchTwoLevels (Array<MF*>)");
609 :
610 : using FAB = typename MF::FABType::value_type;
611 : using iFAB = typename iMultiFab::FABType::value_type;
612 :
613 : AMREX_ASSERT(AMREX_D_TERM(mf[0]->ixType().nodeCentered(0),
614 : && mf[1]->ixType().nodeCentered(1),
615 : && mf[2]->ixType().nodeCentered(2)));
616 :
617 : // These need to be true: (ba[0] == ba[1] == ba[2]) & (dm[0] == dm[1] == dm[2]).
618 : // Debatable whether these are required, or will be enforced elsewhere prior to this func.
619 : AMREX_ASSERT(AMREX_D_TERM(true,
620 : && BoxArray::SameRefs(mf[0]->boxArray(), mf[1]->boxArray()),
621 : && BoxArray::SameRefs(mf[0]->boxArray(), mf[2]->boxArray())));
622 : /*
623 : AMREX_ASSERT(AMREX_D_TERM(true,
624 : && DistributionMapping::SameRefs(mf[0]->DistributionMap(), mf[1]->DistributionMap()),
625 : && DistributionMapping::SameRefs(mf[0]->DistributionMap(), mf[2]->DistributionMap())));
626 : */
627 :
628 :
629 : // Test all of them?
630 : if (nghost.max() > 0 || mf[0]->getBDKey() != fmf[0][0]->getBDKey())
631 : {
632 : const InterpolaterBoxCoarsener& coarsener = mapper->BoxCoarsener(ratio);
633 :
634 : // Convert to cell-centered MF meta-data for FPInfo.
635 : MF mf_cc_dummy( amrex::convert(mf[0]->boxArray(), IntVect::TheZeroVector()),
636 : mf[0]->DistributionMap(), ncomp, nghost, MFInfo().SetAlloc(false) );
637 : MF fmf_cc_dummy( amrex::convert(fmf[0][0]->boxArray(), IntVect::TheZeroVector()),
638 : fmf[0][0]->DistributionMap(), ncomp, nghost, MFInfo().SetAlloc(false) );
639 :
640 : const FabArrayBase::FPinfo& fpc = FabArrayBase::TheFPinfo(fmf_cc_dummy, mf_cc_dummy,
641 : nghost,
642 : coarsener,
643 : fgeom,
644 : cgeom,
645 : index_space);
646 :
647 : if ( !fpc.ba_crse_patch.empty() )
648 : {
649 : Array<MF, AMREX_SPACEDIM> mf_crse_patch;
650 : Array<MF, AMREX_SPACEDIM> mf_refined_patch;
651 : Array<iMultiFab, AMREX_SPACEDIM> solve_mask;
652 :
653 : for (int d=0; d<AMREX_SPACEDIM; ++d)
654 : {
655 : mf_crse_patch[d] = make_mf_crse_patch<MF> (fpc, ncomp, mf[d]->boxArray().ixType());
656 : mf_refined_patch[d] = make_mf_refined_patch<MF> (fpc, ncomp, mf[d]->boxArray().ixType(), ratio);
657 : solve_mask[d] = make_mf_crse_mask<iMultiFab>(fpc, ncomp, mf[d]->boxArray().ixType(), ratio);
658 :
659 : mf_set_domain_bndry(mf_crse_patch[d], cgeom);
660 : Vector<MF*> cmf_time;
661 : for (const auto & mfab : cmf)
662 : { cmf_time.push_back(mfab[d]); }
663 :
664 : FillPatchSingleLevel(mf_crse_patch[d], time, cmf_time, ct, scomp, 0, ncomp,
665 : cgeom, cbc[d], cbccomp[d]);
666 :
667 : mf_set_domain_bndry(mf_refined_patch[d], fgeom);
668 : Vector<MF*> fmf_time;
669 : for (const auto & mfab : fmf)
670 : { fmf_time.push_back(mfab[d]); }
671 :
672 : FillPatchSingleLevel(mf_refined_patch[d], time, fmf_time, ft, scomp, 0, ncomp,
673 : fgeom, fbc[d], fbccomp[d]);
674 :
675 :
676 : // Aliased MFs, used to allow CPC caching.
677 : MF mf_known( amrex::coarsen(fmf[0][d]->boxArray(), ratio), fmf[0][d]->DistributionMap(),
678 : ncomp, nghost, MFInfo().SetAlloc(false) );
679 : MF mf_solution( amrex::coarsen(mf_refined_patch[d].boxArray(), ratio), mf_refined_patch[d].DistributionMap(),
680 : ncomp, 0, MFInfo().SetAlloc(false) );
681 :
682 : const FabArrayBase::CPC mask_cpc( mf_solution, IntVect::TheZeroVector(),
683 : mf_known, IntVect::TheZeroVector(),
684 : cgeom.periodicity() );
685 :
686 : solve_mask[d].setVal(1); // Values to solve.
687 : solve_mask[d].setVal(0, mask_cpc, 0, 1); // Known values.
688 : }
689 :
690 : int idummy=0;
691 : #ifdef AMREX_USE_OMP
692 : // bool cc = fpc.ba_crse_patch.ixType().cellCentered();
693 : bool cc = false; // can anything be done to allow threading, or can the OpenMP just be removed?
694 : #pragma omp parallel if (cc && Gpu::notInLaunchRegion() )
695 : #endif
696 : {
697 : Vector<Array<BCRec, AMREX_SPACEDIM> > bcr(ncomp);
698 : for (MFIter mfi(mf_refined_patch[0]); mfi.isValid(); ++mfi)
699 : {
700 : Array<FAB*, AMREX_SPACEDIM> sfab{ AMREX_D_DECL( &(mf_crse_patch[0][mfi]),
701 : &(mf_crse_patch[1][mfi]),
702 : &(mf_crse_patch[2][mfi]) )};
703 : Array<FAB*, AMREX_SPACEDIM> dfab{ AMREX_D_DECL( &(mf_refined_patch[0][mfi]),
704 : &(mf_refined_patch[1][mfi]),
705 : &(mf_refined_patch[2][mfi]) )};
706 : Array<iFAB*, AMREX_SPACEDIM> mfab{ AMREX_D_DECL( &(solve_mask[0][mfi]),
707 : &(solve_mask[1][mfi]),
708 : &(solve_mask[2][mfi]) )};
709 :
710 : const Box& sbx_cc = amrex::convert(sfab[0]->box(), IntVect::TheZeroVector());
711 : const Box& dbx_cc = amrex::convert(dfab[0]->box(), IntVect::TheZeroVector());
712 :
713 : for (int d=0; d<AMREX_SPACEDIM; ++d)
714 : {
715 : Vector<BCRec> bcr_d(ncomp);
716 :
717 : amrex::setBC(sfab[d]->box(), amrex::convert(cgeom.Domain(), mf[d]->ixType()),
718 : bcscomp[d],0,ncomp,bcs[d],bcr_d);
719 :
720 : for (int n=0; n<ncomp; ++n)
721 : { bcr[n][d] = bcr_d[n]; }
722 : }
723 :
724 : pre_interp(sfab, sbx_cc, 0, ncomp);
725 :
726 : mapper->interp_arr(sfab, 0, dfab, 0, ncomp, dbx_cc, ratio, mfab,
727 : cgeom, fgeom, bcr, idummy, idummy, RunOn::Gpu);
728 :
729 : post_interp(dfab, dbx_cc, 0, ncomp);
730 : }
731 : }
732 :
733 : for (int d=0; d<AMREX_SPACEDIM; ++d)
734 : {
735 : bool aliasing = false;
736 : for (auto const& fmf_a : fmf) {
737 : aliasing = aliasing || (mf[d] == fmf_a[d]);
738 : }
739 : if (aliasing) {
740 : mf[d]->ParallelCopyToGhost(mf_refined_patch[d], 0, dcomp, ncomp,
741 : IntVect{0}, nghost);
742 : } else {
743 : mf[d]->ParallelCopy(mf_refined_patch[d], 0, dcomp, ncomp,
744 : IntVect{0}, nghost);
745 : }
746 : }
747 : }
748 : }
749 :
750 : for (int d=0; d<AMREX_SPACEDIM; ++d)
751 : {
752 : Vector<MF*> fmf_time;
753 : for (auto const& ffab : fmf)
754 : { fmf_time.push_back(ffab[d]); }
755 :
756 : FillPatchSingleLevel(*mf[d], nghost, time, fmf_time, ft, scomp, dcomp, ncomp,
757 : fgeom, fbc[d], fbccomp[d]);
758 : }
759 : }
760 :
761 : } // namespace detail
762 :
763 : template <typename MF, typename BC, typename Interp, typename PreInterpHook, typename PostInterpHook>
764 : std::enable_if_t<IsFabArray<MF>::value>
765 : FillPatchTwoLevels (MF& mf, IntVect const& nghost, Real time,
766 : const Vector<MF*>& cmf, const Vector<Real>& ct,
767 : const Vector<MF*>& fmf, const Vector<Real>& ft,
768 : int scomp, int dcomp, int ncomp,
769 : const Geometry& cgeom, const Geometry& fgeom,
770 : BC& cbc, int cbccomp,
771 : BC& fbc, int fbccomp,
772 : const IntVect& ratio,
773 : Interp* mapper,
774 : const Vector<BCRec>& bcs, int bcscomp,
775 : const PreInterpHook& pre_interp,
776 : const PostInterpHook& post_interp)
777 : {
778 : #ifdef AMREX_USE_EB
779 : EB2::IndexSpace const* index_space = EB2::TopIndexSpaceIfPresent();
780 : #else
781 : EB2::IndexSpace const* index_space = nullptr;
782 : #endif
783 : detail::FillPatchTwoLevels_doit(mf,nghost,time,cmf,ct,fmf,ft,
784 : scomp,dcomp,ncomp,cgeom,fgeom,
785 : cbc,cbccomp,fbc,fbccomp,ratio,mapper,bcs,bcscomp,
786 : pre_interp,post_interp,index_space);
787 : }
788 :
789 : template <typename MF, typename BC, typename Interp, typename PreInterpHook, typename PostInterpHook>
790 : std::enable_if_t<IsFabArray<MF>::value>
791 66 : FillPatchTwoLevels (MF& mf, Real time,
792 : const Vector<MF*>& cmf, const Vector<Real>& ct,
793 : const Vector<MF*>& fmf, const Vector<Real>& ft,
794 : int scomp, int dcomp, int ncomp,
795 : const Geometry& cgeom, const Geometry& fgeom,
796 : BC& cbc, int cbccomp,
797 : BC& fbc, int fbccomp,
798 : const IntVect& ratio,
799 : Interp* mapper,
800 : const Vector<BCRec>& bcs, int bcscomp,
801 : const PreInterpHook& pre_interp,
802 : const PostInterpHook& post_interp)
803 : {
804 : #ifdef AMREX_USE_EB
805 : EB2::IndexSpace const* index_space = EB2::TopIndexSpaceIfPresent();
806 : #else
807 66 : EB2::IndexSpace const* index_space = nullptr;
808 : #endif
809 :
810 66 : detail::FillPatchTwoLevels_doit(mf,mf.nGrowVect(),time,cmf,ct,fmf,ft,
811 : scomp,dcomp,ncomp,cgeom,fgeom,
812 : cbc,cbccomp,fbc,fbccomp,ratio,mapper,bcs,bcscomp,
813 : pre_interp,post_interp,index_space);
814 66 : }
815 :
816 : template <typename MF, typename BC, typename Interp, typename PreInterpHook, typename PostInterpHook>
817 : std::enable_if_t<IsFabArray<MF>::value>
818 : FillPatchTwoLevels (Array<MF*, AMREX_SPACEDIM> const& mf, IntVect const& nghost, Real time,
819 : const Vector<Array<MF*, AMREX_SPACEDIM> >& cmf, const Vector<Real>& ct,
820 : const Vector<Array<MF*, AMREX_SPACEDIM> >& fmf, const Vector<Real>& ft,
821 : int scomp, int dcomp, int ncomp,
822 : const Geometry& cgeom, const Geometry& fgeom,
823 : Array<BC, AMREX_SPACEDIM>& cbc, const Array<int, AMREX_SPACEDIM>& cbccomp,
824 : Array<BC, AMREX_SPACEDIM>& fbc, const Array<int, AMREX_SPACEDIM>& fbccomp,
825 : const IntVect& ratio,
826 : Interp* mapper,
827 : const Array<Vector<BCRec>, AMREX_SPACEDIM>& bcs, const Array<int, AMREX_SPACEDIM>& bcscomp,
828 : const PreInterpHook& pre_interp,
829 : const PostInterpHook& post_interp)
830 : {
831 : #ifdef AMREX_USE_EB
832 : EB2::IndexSpace const* index_space = EB2::TopIndexSpaceIfPresent();
833 : #else
834 : EB2::IndexSpace const* index_space = nullptr;
835 : #endif
836 :
837 : detail::FillPatchTwoLevels_doit(mf,nghost,time,cmf,ct,fmf,ft,
838 : scomp,dcomp,ncomp,cgeom,fgeom,
839 : cbc,cbccomp,fbc,fbccomp,ratio,mapper,bcs,bcscomp,
840 : pre_interp,post_interp,index_space);
841 : }
842 :
843 : template <typename MF, typename BC, typename Interp, typename PreInterpHook, typename PostInterpHook>
844 : std::enable_if_t<IsFabArray<MF>::value>
845 : FillPatchTwoLevels (Array<MF*, AMREX_SPACEDIM> const& mf, IntVect const& nghost, Real time,
846 : const Vector<Array<MF*, AMREX_SPACEDIM> >& cmf, const Vector<Real>& ct,
847 : const Vector<Array<MF*, AMREX_SPACEDIM> >& fmf, const Vector<Real>& ft,
848 : int scomp, int dcomp, int ncomp,
849 : const Geometry& cgeom, const Geometry& fgeom,
850 : Array<BC, AMREX_SPACEDIM>& cbc, int cbccomp,
851 : Array<BC, AMREX_SPACEDIM>& fbc, int fbccomp,
852 : const IntVect& ratio,
853 : Interp* mapper,
854 : const Array<Vector<BCRec>, AMREX_SPACEDIM>& bcs, int bcscomp,
855 : const PreInterpHook& pre_interp,
856 : const PostInterpHook& post_interp)
857 : {
858 : #ifdef AMREX_USE_EB
859 : EB2::IndexSpace const* index_space = EB2::TopIndexSpaceIfPresent();
860 : #else
861 : EB2::IndexSpace const* index_space = nullptr;
862 : #endif
863 :
864 : Array<int, AMREX_SPACEDIM> cbccomp_arr = {AMREX_D_DECL(cbccomp,cbccomp,cbccomp)};
865 : Array<int, AMREX_SPACEDIM> fbccomp_arr = {AMREX_D_DECL(fbccomp,fbccomp,fbccomp)};
866 : Array<int, AMREX_SPACEDIM> bcscomp_arr = {AMREX_D_DECL(bcscomp,bcscomp,bcscomp)};
867 :
868 : detail::FillPatchTwoLevels_doit(mf,nghost,time,cmf,ct,fmf,ft,
869 : scomp,dcomp,ncomp,cgeom,fgeom,
870 : cbc,cbccomp_arr,fbc,fbccomp_arr,ratio,mapper,bcs,bcscomp_arr,
871 : pre_interp,post_interp,index_space);
872 : }
873 :
874 : template <typename MF, typename BC, typename Interp, typename PreInterpHook, typename PostInterpHook>
875 : std::enable_if_t<IsFabArray<MF>::value>
876 : FillPatchTwoLevels (Array<MF*, AMREX_SPACEDIM> const& mf, Real time,
877 : const Vector<Array<MF*, AMREX_SPACEDIM> >& cmf, const Vector<Real>& ct,
878 : const Vector<Array<MF*, AMREX_SPACEDIM> >& fmf, const Vector<Real>& ft,
879 : int scomp, int dcomp, int ncomp,
880 : const Geometry& cgeom, const Geometry& fgeom,
881 : Array<BC, AMREX_SPACEDIM>& cbc, int cbccomp,
882 : Array<BC, AMREX_SPACEDIM>& fbc, int fbccomp,
883 : const IntVect& ratio,
884 : Interp* mapper,
885 : const Array<Vector<BCRec>, AMREX_SPACEDIM>& bcs, int bcscomp,
886 : const PreInterpHook& pre_interp,
887 : const PostInterpHook& post_interp)
888 : {
889 : #ifdef AMREX_USE_EB
890 : EB2::IndexSpace const* index_space = EB2::TopIndexSpaceIfPresent();
891 : #else
892 : EB2::IndexSpace const* index_space = nullptr;
893 : #endif
894 :
895 : Array<int, AMREX_SPACEDIM> cbccomp_arr = {AMREX_D_DECL(cbccomp,cbccomp,cbccomp)};
896 : Array<int, AMREX_SPACEDIM> fbccomp_arr = {AMREX_D_DECL(fbccomp,fbccomp,fbccomp)};
897 : Array<int, AMREX_SPACEDIM> bcscomp_arr = {AMREX_D_DECL(bcscomp,bcscomp,bcscomp)};
898 :
899 : detail::FillPatchTwoLevels_doit(mf,mf[0]->nGrowVect(),time,cmf,ct,fmf,ft,
900 : scomp,dcomp,ncomp,cgeom,fgeom,
901 : cbc,cbccomp_arr,fbc,fbccomp_arr,ratio,mapper,bcs,bcscomp_arr,
902 : pre_interp,post_interp,index_space);
903 : }
904 :
905 : #ifdef AMREX_USE_EB
906 : template <typename MF, typename BC, typename Interp, typename PreInterpHook, typename PostInterpHook>
907 : std::enable_if_t<IsFabArray<MF>::value>
908 : FillPatchTwoLevels (MF& mf, IntVect const& nghost, Real time,
909 : const EB2::IndexSpace& index_space,
910 : const Vector<MF*>& cmf, const Vector<Real>& ct,
911 : const Vector<MF*>& fmf, const Vector<Real>& ft,
912 : int scomp, int dcomp, int ncomp,
913 : const Geometry& cgeom, const Geometry& fgeom,
914 : BC& cbc, int cbccomp,
915 : BC& fbc, int fbccomp,
916 : const IntVect& ratio,
917 : Interp* mapper,
918 : const Vector<BCRec>& bcs, int bcscomp,
919 : const PreInterpHook& pre_interp,
920 : const PostInterpHook& post_interp)
921 : {
922 : detail::FillPatchTwoLevels_doit(mf,nghost,time,cmf,ct,fmf,ft,
923 : scomp,dcomp,ncomp,cgeom,fgeom,
924 : cbc,cbccomp,fbc,fbccomp,ratio,mapper,bcs,bcscomp,
925 : pre_interp,post_interp,&index_space);
926 : }
927 :
928 : template <typename MF, typename BC, typename Interp, typename PreInterpHook, typename PostInterpHook>
929 : std::enable_if_t<IsFabArray<MF>::value>
930 : FillPatchTwoLevels (MF& mf, Real time,
931 : const EB2::IndexSpace& index_space,
932 : const Vector<MF*>& cmf, const Vector<Real>& ct,
933 : const Vector<MF*>& fmf, const Vector<Real>& ft,
934 : int scomp, int dcomp, int ncomp,
935 : const Geometry& cgeom, const Geometry& fgeom,
936 : BC& cbc, int cbccomp,
937 : BC& fbc, int fbccomp,
938 : const IntVect& ratio,
939 : Interp* mapper,
940 : const Vector<BCRec>& bcs, int bcscomp,
941 : const PreInterpHook& pre_interp,
942 : const PostInterpHook& post_interp)
943 : {
944 : detail::FillPatchTwoLevels_doit(mf,mf.nGrowVect(),time,cmf,ct,fmf,ft,
945 : scomp,dcomp,ncomp,cgeom,fgeom,
946 : cbc,cbccomp,fbc,fbccomp,ratio,mapper,bcs,bcscomp,
947 : pre_interp,post_interp,&index_space);
948 : }
949 : #endif
950 :
951 : template <typename MF, typename BC, typename Interp, typename PreInterpHook, typename PostInterpHook>
952 : std::enable_if_t<IsFabArray<MF>::value>
953 0 : InterpFromCoarseLevel (MF& mf, Real time,
954 : const MF& cmf, int scomp, int dcomp, int ncomp,
955 : const Geometry& cgeom, const Geometry& fgeom,
956 : BC& cbc, int cbccomp,
957 : BC& fbc, int fbccomp,
958 : const IntVect& ratio,
959 : Interp* mapper,
960 : const Vector<BCRec>& bcs, int bcscomp,
961 : const PreInterpHook& pre_interp,
962 : const PostInterpHook& post_interp)
963 : {
964 0 : InterpFromCoarseLevel(mf,mf.nGrowVect(),time,cmf,scomp,dcomp,ncomp,cgeom,fgeom,
965 : cbc,cbccomp,fbc,fbccomp,ratio,mapper,bcs,bcscomp,
966 : pre_interp,post_interp);
967 0 : }
968 :
969 : template <typename MF, typename BC, typename Interp, typename PreInterpHook, typename PostInterpHook>
970 : std::enable_if_t<IsFabArray<MF>::value>
971 : InterpFromCoarseLevel (Array<MF*, AMREX_SPACEDIM> const& mf, Real time,
972 : const Array<MF*, AMREX_SPACEDIM>& cmf, int scomp, int dcomp, int ncomp,
973 : const Geometry& cgeom, const Geometry& fgeom,
974 : Array<BC, AMREX_SPACEDIM>& cbc, int cbccomp,
975 : Array<BC, AMREX_SPACEDIM>& fbc, int fbccomp,
976 : const IntVect& ratio,
977 : Interp* mapper,
978 : const Array<Vector<BCRec>, AMREX_SPACEDIM>& bcs, int bcscomp,
979 : const PreInterpHook& pre_interp,
980 : const PostInterpHook& post_interp)
981 : {
982 : InterpFromCoarseLevel(mf,mf[0]->nGrowVect(),time,cmf,scomp,dcomp,ncomp,cgeom,fgeom,
983 : cbc,cbccomp,fbc,fbccomp,ratio,mapper,bcs,bcscomp,
984 : pre_interp,post_interp);
985 : }
986 :
987 : template <typename MF, typename BC, typename Interp, typename PreInterpHook, typename PostInterpHook>
988 : std::enable_if_t<IsFabArray<MF>::value>
989 0 : InterpFromCoarseLevel (MF& mf, IntVect const& nghost, Real time,
990 : const MF& cmf, int scomp, int dcomp, int ncomp,
991 : const Geometry& cgeom, const Geometry& fgeom,
992 : BC& cbc, int cbccomp,
993 : BC& fbc, int fbccomp,
994 : const IntVect& ratio,
995 : Interp* mapper,
996 : const Vector<BCRec>& bcs, int bcscomp,
997 : const PreInterpHook& pre_interp,
998 : const PostInterpHook& post_interp)
999 : {
1000 : using FAB = typename MF::FABType::value_type;
1001 :
1002 0 : const InterpolaterBoxCoarsener& coarsener = mapper->BoxCoarsener(ratio);
1003 :
1004 0 : const BoxArray& ba = mf.boxArray();
1005 0 : const DistributionMapping& dm = mf.DistributionMap();
1006 :
1007 0 : const IndexType& typ = ba.ixType();
1008 :
1009 : BL_ASSERT(typ == cmf.boxArray().ixType());
1010 :
1011 0 : Box fdomain_g( amrex::convert(fgeom.Domain(),mf.ixType()) );
1012 0 : for (int i = 0; i < AMREX_SPACEDIM; ++i) {
1013 0 : if (fgeom.isPeriodic(i)) {
1014 0 : fdomain_g.grow(i,nghost[i]);
1015 : }
1016 : }
1017 :
1018 0 : BoxArray ba_crse_patch(ba.size());
1019 : { // TODO: later we might want to cache this
1020 0 : for (int i = 0, N = ba.size(); i < N; ++i)
1021 : {
1022 0 : Box bx = amrex::convert(amrex::grow(ba[i],nghost), typ);
1023 0 : bx &= fdomain_g;
1024 0 : ba_crse_patch.set(i, coarsener.doit(bx));
1025 : }
1026 : }
1027 :
1028 0 : MF mf_crse_patch;
1029 : #ifdef AMREX_USE_EB
1030 : if (EB2::TopIndexSpaceIfPresent()) {
1031 : auto factory = makeEBFabFactory(cgeom, ba_crse_patch, dm, {0,0,0}, EBSupport::basic);
1032 : mf_crse_patch.define(ba_crse_patch, dm, ncomp, 0, MFInfo(), *factory);
1033 : } else
1034 : #endif
1035 : {
1036 0 : mf_crse_patch.define(ba_crse_patch, dm, ncomp, 0);
1037 : }
1038 0 : detail::mf_set_domain_bndry (mf_crse_patch, cgeom);
1039 :
1040 0 : mf_crse_patch.ParallelCopy(cmf, scomp, 0, ncomp, cgeom.periodicity());
1041 :
1042 0 : cbc(mf_crse_patch, 0, ncomp, mf_crse_patch.nGrowVect(), time, cbccomp);
1043 :
1044 0 : detail::call_interp_hook(pre_interp, mf_crse_patch, 0, ncomp);
1045 :
1046 0 : FillPatchInterp(mf, dcomp, mf_crse_patch, 0, ncomp, nghost, cgeom, fgeom, fdomain_g,
1047 : ratio, mapper, bcs, bcscomp);
1048 :
1049 : #ifdef AMREX_USE_OMP
1050 : #pragma omp parallel if (Gpu::notInLaunchRegion())
1051 : #endif
1052 0 : for (MFIter mfi(mf); mfi.isValid(); ++mfi)
1053 : {
1054 0 : FAB& dfab = mf[mfi];
1055 0 : Box dfab_bx = dfab.box();
1056 0 : dfab_bx.grow(nghost-mf.nGrowVect());
1057 0 : const Box& dbx = dfab_bx & fdomain_g;
1058 :
1059 0 : post_interp(dfab, dbx, dcomp, ncomp);
1060 : }
1061 :
1062 0 : fbc(mf, dcomp, ncomp, nghost, time, fbccomp);
1063 0 : }
1064 :
1065 : template <typename MF, typename BC, typename Interp, typename PreInterpHook, typename PostInterpHook>
1066 : std::enable_if_t<IsFabArray<MF>::value>
1067 : InterpFromCoarseLevel (Array<MF*, AMREX_SPACEDIM> const& mf, IntVect const& nghost, Real time,
1068 : const Array<MF*, AMREX_SPACEDIM>& cmf, int scomp, int dcomp, int ncomp,
1069 : const Geometry& cgeom, const Geometry& fgeom,
1070 : Array<BC, AMREX_SPACEDIM>& cbc, int cbccomp,
1071 : Array<BC, AMREX_SPACEDIM>& fbc, int fbccomp,
1072 : const IntVect& ratio,
1073 : Interp* mapper,
1074 : const Array<Vector<BCRec>, AMREX_SPACEDIM>& bcs, int bcscomp,
1075 : const PreInterpHook& pre_interp,
1076 : const PostInterpHook& post_interp)
1077 : {
1078 : using FAB = typename MF::FABType::value_type;
1079 : using iFAB = typename iMultiFab::FABType::value_type;
1080 :
1081 : const InterpolaterBoxCoarsener& coarsener = mapper->BoxCoarsener(ratio);
1082 : const BoxArray& ba = mf[0]->boxArray();
1083 : const DistributionMapping& dm = mf[0]->DistributionMap();
1084 :
1085 : AMREX_ASSERT(AMREX_D_TERM(mf[0]->ixType().nodeCentered(0),
1086 : && mf[1]->ixType().nodeCentered(1),
1087 : && mf[2]->ixType().nodeCentered(2)));
1088 :
1089 : // These need to be true: (ba[0] == ba[1] == ba[2]) & (dm[0] == dm[1] == dm[2]).
1090 : // Debatable whether these are required, or will be enforced elsewhere prior to this func.
1091 : AMREX_ASSERT(AMREX_D_TERM(true,
1092 : && BoxArray::SameRefs(mf[0]->boxArray(), mf[1]->boxArray()),
1093 : && BoxArray::SameRefs(mf[0]->boxArray(), mf[2]->boxArray())));
1094 : /*
1095 : AMREX_ASSERT(AMREX_D_TERM(true,
1096 : && DistributionMapping::SameRefs(mf[0]->DistributionMap(), mf[1]->DistributionMap()),
1097 : && DistributionMapping::SameRefs(mf[0]->DistributionMap(), mf[2]->DistributionMap())));
1098 : */
1099 :
1100 : // If needed, adjust to fully overlap the coarse cells.
1101 : IntVect nghost_adj = nghost;
1102 : for (int d=0; d<AMREX_SPACEDIM; ++d) {
1103 : if (nghost[d] % ratio[d] != 0) {
1104 : nghost_adj[d] += ratio[d] - (nghost[d] % ratio[d]);
1105 : }
1106 : }
1107 :
1108 : Array<MF*, AMREX_SPACEDIM> mf_local = mf;
1109 : int dcomp_adj = dcomp;
1110 : Array<std::unique_ptr<MF>, AMREX_SPACEDIM> mf_temp;
1111 : if (! nghost.allGE(nghost_adj)) {
1112 : for (int d=0; d<AMREX_SPACEDIM; ++d) {
1113 : mf_temp[d] = std::make_unique<MF>(mf[d]->boxArray(),
1114 : mf[d]->DistributionMap(), ncomp, nghost_adj);
1115 : mf_local[d] = mf_temp[d].get();
1116 : }
1117 : dcomp_adj = 0;
1118 : }
1119 :
1120 : // Create a cell-centered boxArray of the region to interp.
1121 : // Convert this boxArray and domain as needed.
1122 : Box fdomain = amrex::convert(fgeom.Domain(), IntVect::TheZeroVector());
1123 : Box fdomain_g(fdomain);
1124 : for (int d = 0; d < AMREX_SPACEDIM; ++d) {
1125 : if (fgeom.isPeriodic(d)) {
1126 : fdomain_g.grow(d,nghost_adj[d]);
1127 : }
1128 : }
1129 :
1130 : // Build patches, using domain to account for periodic bcs.
1131 : BoxArray ba_crse_patch(ba.size());
1132 : { // TODO: later we might want to cache this
1133 : for (int i = 0, N = ba.size(); i < N; ++i)
1134 : {
1135 : Box bx = amrex::convert(amrex::grow(ba[i], nghost_adj), IntVect::TheZeroVector());
1136 : bx &= fdomain_g;
1137 : ba_crse_patch.set(i, coarsener.doit(bx));
1138 : }
1139 : }
1140 :
1141 : Array<MF, AMREX_SPACEDIM> mf_crse_patch;
1142 : for (int d = 0; d<AMREX_SPACEDIM; ++d)
1143 : {
1144 : IndexType typ = mf[d]->boxArray().ixType();
1145 : BoxArray ba_crse_idxed = amrex::convert(ba_crse_patch, typ);
1146 :
1147 : #ifdef AMREX_USE_EB
1148 : auto crse_factory = makeEBFabFactory(cgeom, ba_crse_idxed, dm, {0,0,0}, EBSupport::basic);
1149 : mf_crse_patch[d].define(ba_crse_idxed, dm, ncomp, 0, MFInfo(), *crse_factory);
1150 : #else
1151 : mf_crse_patch[d].define(ba_crse_idxed, dm, ncomp, 0);
1152 : #endif
1153 : detail::mf_set_domain_bndry(mf_crse_patch[d], cgeom);
1154 :
1155 : mf_crse_patch[d].ParallelCopy(*(cmf[d]), scomp, 0, ncomp, cgeom.periodicity());
1156 : cbc[d](mf_crse_patch[d], 0, ncomp, mf_crse_patch[d].nGrowVect(), time, cbccomp);
1157 : }
1158 :
1159 : int idummy1=0, idummy2=0;
1160 : #ifdef AMREX_USE_OMP
1161 : #pragma omp parallel if (Gpu::notInLaunchRegion())
1162 : #endif
1163 : {
1164 : Vector<Array<BCRec, AMREX_SPACEDIM> > bcr(ncomp);
1165 :
1166 : // Empty containers describing that all points must be solved (no mask).
1167 : Array<iFAB*, AMREX_SPACEDIM> mfab{ AMREX_D_DECL( nullptr, nullptr, nullptr ) };
1168 :
1169 : for (MFIter mfi(mf_crse_patch[0]); mfi.isValid(); ++mfi)
1170 : {
1171 : Array<FAB*, AMREX_SPACEDIM> sfab{ AMREX_D_DECL( &(mf_crse_patch[0][mfi]),
1172 : &(mf_crse_patch[1][mfi]),
1173 : &(mf_crse_patch[2][mfi]) )};
1174 : Array<FAB*, AMREX_SPACEDIM> dfab{ AMREX_D_DECL( &(*mf_local[0])[mfi],
1175 : &(*mf_local[1])[mfi],
1176 : &(*mf_local[2])[mfi] )};
1177 :
1178 : const Box& sbx_cc = amrex::convert(sfab[0]->box(), IntVect::TheZeroVector());
1179 : Box dfab_cc = amrex::convert(dfab[0]->box(), IntVect::TheZeroVector());
1180 : const Box& dbx_cc = dfab_cc & fdomain_g;
1181 :
1182 : for (int d=0; d<AMREX_SPACEDIM; ++d)
1183 : {
1184 : Vector<BCRec> bcr_d(ncomp);
1185 :
1186 : amrex::setBC(sfab[d]->box(),
1187 : amrex::convert(cgeom.Domain(), sfab[d]->box().ixType()),
1188 : bcscomp,0,ncomp,bcs[d],bcr_d);
1189 :
1190 : for (int n=0; n<ncomp; ++n)
1191 : { bcr[n][d] = bcr_d[n]; }
1192 : }
1193 :
1194 : pre_interp(sfab, sbx_cc, 0, ncomp);
1195 :
1196 : mapper->interp_arr(sfab, 0, dfab, 0, ncomp, dbx_cc, ratio, mfab,
1197 : cgeom, fgeom, bcr, idummy1, idummy2, RunOn::Gpu);
1198 :
1199 : post_interp(dfab, dbx_cc, 0, ncomp);
1200 : }
1201 : }
1202 :
1203 : for (int d=0; d<AMREX_SPACEDIM; ++d)
1204 : {
1205 : if (mf[d] != mf_local[d]) {
1206 : amrex::Copy(*mf[d], *mf_local[d], 0, dcomp_adj, ncomp, nghost);
1207 : }
1208 :
1209 : fbc[d](*mf[d], dcomp, ncomp, nghost, time, fbccomp);
1210 : }
1211 : }
1212 :
1213 : template <typename MF, typename BC, typename Interp>
1214 : std::enable_if_t<IsFabArray<MF>::value>
1215 : FillPatchNLevels (MF& mf, int level, const IntVect& nghost, Real time,
1216 : const Vector<Vector<MF*>>& smf, const Vector<Vector<Real>>& st,
1217 : int scomp, int dcomp, int ncomp,
1218 : const Vector<Geometry>& geom,
1219 : Vector<BC>& bc, int bccomp,
1220 : const Vector<IntVect>& ratio,
1221 : Interp* mapper,
1222 : const Vector<BCRec>& bcr, int bcrcomp)
1223 : {
1224 : BL_PROFILE("FillPatchNLevels");
1225 :
1226 : #ifdef AMREX_USE_EB
1227 : EB2::IndexSpace const* index_space = EB2::TopIndexSpaceIfPresent();
1228 : #else
1229 : EB2::IndexSpace const* index_space = nullptr;
1230 : #endif
1231 :
1232 : AMREX_ALWAYS_ASSERT(level < int(geom.size()) &&
1233 : level < int(bc.size()) &&
1234 : level < int(ratio.size()+1));
1235 : if (level == 0) {
1236 : FillPatchSingleLevel(mf, nghost, time, smf[0], st[0], scomp, dcomp, ncomp, geom[0],
1237 : bc[0], bccomp);
1238 : } else if (level >= int(smf.size())) {
1239 : BoxArray const& ba = mf.boxArray();
1240 : auto const& typ = ba.ixType();
1241 : Box domain_g = geom[level].growPeriodicDomain(nghost);
1242 : domain_g.convert(typ);
1243 : BoxArray cba;
1244 : {
1245 : BoxList cbl(typ);
1246 : cbl.reserve(ba.size());
1247 : for (int i = 0, N = int(ba.size()); i < N; ++i) {
1248 : cbl.push_back(mapper->CoarseBox(amrex::grow(ba[i],nghost), ratio[level-1]));
1249 : }
1250 : cba = BoxArray(std::move(cbl));
1251 : }
1252 : MultiFab cmf_tmp;
1253 : #ifdef AMREX_USE_EB
1254 : if (index_space) {
1255 : auto factory = makeEBFabFactory(index_space, geom[level-1], cba,
1256 : mf.DistributionMap(), {0,0,0},
1257 : EBSupport::basic);
1258 : cmf_tmp.define(cba, mf.DistributionMap(), ncomp, 0, MFInfo(), *factory);
1259 : } else
1260 : #endif
1261 : {
1262 : cmf_tmp.define(cba, mf.DistributionMap(), ncomp, 0);
1263 : }
1264 : FillPatchNLevels(cmf_tmp, level-1, IntVect(0), time, smf, st, scomp, 0, ncomp,
1265 : geom, bc, bccomp, ratio, mapper, bcr, bcrcomp);
1266 : FillPatchInterp(mf, dcomp, cmf_tmp, 0, ncomp, nghost, geom[level-1], geom[level],
1267 : domain_g, ratio[level-1], mapper, bcr, bcrcomp);
1268 : } else {
1269 : NullInterpHook<typename MF::FABType::value_type> hook{};
1270 : int error_code = detail::FillPatchTwoLevels_doit(mf, nghost, time,
1271 : smf[level-1], st[level-1],
1272 : smf[level ], st[level ],
1273 : scomp, dcomp, ncomp,
1274 : geom[level-1], geom[level],
1275 : bc[level-1], bccomp,
1276 : bc[level ], bccomp,
1277 : ratio[level-1], mapper, bcr, bcrcomp,
1278 : hook, hook, index_space, true);
1279 : if (error_code == 0) { return; }
1280 :
1281 : BoxArray cba;
1282 : {
1283 : BoxArray const& ba = mf.boxArray();
1284 : BoxList cbl(mf.ixType());
1285 : cbl.reserve(ba.size());
1286 : for (int i = 0; i < int(ba.size()); ++i) {
1287 : cbl.push_back(mapper->CoarseBox(amrex::grow(ba[i], nghost), ratio[level-1]));
1288 : }
1289 : cba = BoxArray(std::move(cbl));
1290 : }
1291 : MultiFab cmf_tmp;
1292 : #ifdef AMREX_USE_EB
1293 : if (index_space) {
1294 : auto factory = makeEBFabFactory(index_space, geom[level-1], cba,
1295 : mf.DistributionMap(), {0,0,0},
1296 : EBSupport::basic);
1297 : cmf_tmp.define(cba, mf.DistributionMap(), ncomp, 0, MFInfo(), *factory);
1298 : } else
1299 : #endif
1300 : {
1301 : cmf_tmp.define(cba, mf.DistributionMap(), ncomp, 0);
1302 : }
1303 : FillPatchNLevels(cmf_tmp, level-1, IntVect(0), time, smf, st, scomp, 0, ncomp,
1304 : geom, bc, bccomp, ratio, mapper, bcr, bcrcomp);
1305 : Vector<MF*> cmf{&cmf_tmp};
1306 : Vector<MF*> fmf = smf[level];
1307 : Vector<MF> fmf_raii;
1308 : if (scomp != 0) {
1309 : for (auto const* p : fmf) {
1310 : fmf_raii.emplace_back(*p, amrex::make_alias, scomp, ncomp);
1311 : }
1312 : }
1313 : detail::FillPatchTwoLevels_doit(mf, nghost, time,
1314 : cmf, {time},
1315 : fmf, st[level],
1316 : 0, dcomp, ncomp,
1317 : geom[level-1], geom[level],
1318 : bc[level-1], bccomp,
1319 : bc[level ], bccomp,
1320 : ratio[level-1], mapper, bcr, bccomp,
1321 : hook, hook, index_space);
1322 : }
1323 : }
1324 :
1325 : }
1326 :
1327 : #endif
|
