ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
Public Member Functions | Private Attributes | List of all members
eb_aux_ Class Reference

#include <ERF_EBAux.H>

Collaboration diagram for eb_aux_:
Collaboration graph
[legend]

Public Member Functions

 eb_aux_ ()
 
 ~eb_aux_ ()
 
void define (int const &a_level, int const &a_idim, amrex::Geometry const &a_geom, amrex::BoxArray const &a_grids, amrex::DistributionMapping const &a_dmap, amrex::Vector< int > const &a_ngrow, amrex::EBFArrayBoxFactory const *a_factory)
 
void set_verbose ()
 
const amrex::FabArray< amrex::EBCellFlagFab > & getMultiEBCellFlagFab () const
 
const amrex::MultiFab & getVolFrac () const
 
const amrex::MultiFab & getCentroid () const
 
const amrex::MultiFab & getBndryArea () const
 
const amrex::MultiFab & getBndryCent () const
 
const amrex::MultiFab & getBndryNorm () const
 
amrex::Array< const amrex::MultiFab *, AMREX_SPACEDIM > getAreaFrac () const
 
amrex::Array< const amrex::MultiFab *, AMREX_SPACEDIM > getFaceCent () const
 

Private Attributes

int m_verbose
 
amrex::FabArray< amrex::EBCellFlagFab > * m_cellflags = nullptr
 
amrex::MultiFab * m_volfrac = nullptr
 
amrex::MultiFab * m_volcent = nullptr
 
amrex::MultiFab * m_bndryarea = nullptr
 
amrex::MultiFab * m_bndrycent = nullptr
 
amrex::MultiFab * m_bndrynorm = nullptr
 
amrex::Array< amrex::MultiFab *, AMREX_SPACEDIM > m_areafrac {{AMREX_D_DECL(nullptr, nullptr, nullptr)}}
 
amrex::Array< amrex::MultiFab *, AMREX_SPACEDIM > m_facecent {{AMREX_D_DECL(nullptr, nullptr, nullptr)}}
 

Constructor & Destructor Documentation

◆ eb_aux_()

eb_aux_::eb_aux_ ( )
16  : m_verbose(0)
17 // ,m_defined(0)
18 {}
eb_aux_::m_verbose
int m_verbose
Definition: ERF_EBAux.H:40

◆ ~eb_aux_()

eb_aux_::~eb_aux_ ( )
11 {
12 }

Member Function Documentation

◆ define()

void eb_aux_::define ( int const &  a_level,
int const &  a_idim,
amrex::Geometry const &  a_geom,
amrex::BoxArray const &  a_grids,
amrex::DistributionMapping const &  a_dmap,
amrex::Vector< int > const &  a_ngrow,
amrex::EBFArrayBoxFactory const *  a_factory 
)
29 {
30  // Box dbox(a_geom.Domain());
31 
32  // small_volfrac
33  Real small_volfrac = 1.e-14;
34  ParmParse pp("eb2");
35  pp.queryAdd("small_volfrac", small_volfrac);
36  const Real small_value = 1.e-15;
37 
38  const IntVect vdim(IntVect::TheDimensionVector(a_idim));
39 
40  const BoxArray& grids = amrex::convert(a_grids, vdim);
41 
42  m_cellflags = new FabArray<EBCellFlagFab>(grids, a_dmap, 1, a_ngrow[0], MFInfo(),
43  DefaultFabFactory<EBCellFlagFab>());
44 
45  // Set m_cellflags type to singlevalued
46  m_cellflags->setVal(EBCellFlag::TheDefaultCell());
47  for (MFIter mfi(*m_cellflags, false); mfi.isValid(); ++mfi) {
48  auto& fab = (*m_cellflags)[mfi];
49  fab.setType(FabType::singlevalued);
50  }
51 
52  m_volfrac = new MultiFab(grids, a_dmap, 1, a_ngrow[1], MFInfo(), FArrayBoxFactory());
53  m_volcent = new MultiFab(grids, a_dmap, AMREX_SPACEDIM, a_ngrow[2], MFInfo(), FArrayBoxFactory());
54 
55  for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
56  const BoxArray& faceba = amrex::convert(a_grids, IntVect::TheDimensionVector(idim));
57  if (idim == a_idim) {
58  m_areafrac[idim] = new MultiFab(a_grids, a_dmap, 1, a_ngrow[1]+1, MFInfo(), FArrayBoxFactory());
59  m_facecent[idim] = new MultiFab(a_grids, a_dmap, AMREX_SPACEDIM-1, a_ngrow[2], MFInfo(), FArrayBoxFactory());
60  } else {
61  m_areafrac[idim] = new MultiFab(faceba, a_dmap, 1, a_ngrow[1], MFInfo(), FArrayBoxFactory());
62  m_facecent[idim] = new MultiFab(faceba, a_dmap, AMREX_SPACEDIM-1, a_ngrow[2], MFInfo(), FArrayBoxFactory());
63  }
64  }
65 
66  m_bndryarea = new MultiFab(grids, a_dmap, 1, a_ngrow[2], MFInfo(), FArrayBoxFactory());
67  m_bndrycent = new MultiFab(grids, a_dmap, AMREX_SPACEDIM, a_ngrow[2], MFInfo(), FArrayBoxFactory());
68  m_bndrynorm = new MultiFab(grids, a_dmap, AMREX_SPACEDIM, a_ngrow[2], MFInfo(), FArrayBoxFactory());
69 
70  // Initialize with zeros
71  m_volfrac->setVal(0.0);
72  m_volcent->setVal(0.0);
73 
74  for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
75  m_areafrac[idim]->setVal(0.0);
76  m_facecent[idim]->setVal(0.0);
77  }
78 
79  m_bndryarea->setVal(0.0);
80  m_bndrycent->setVal(0.0);
81  m_bndrynorm->setVal(0.0);
82 
83  const auto& FlagFab = a_factory->getMultiEBCellFlagFab(); // EBFArrayBoxFactory, EBDataCollection
84 
85  for (MFIter mfi(*m_cellflags, false); mfi.isValid(); ++mfi) {
86 
87  const Box& bx = mfi.validbox();
88  const Box& bx_grown = mfi.growntilebox();
89  const Box tbx = mfi.nodaltilebox(a_idim);
90  const Box domain = surroundingNodes(a_geom.Domain(), a_idim);
91 
92  GpuArray<Real, AMREX_SPACEDIM> dx = a_geom.CellSizeArray();
93  bool l_periodic = a_geom.isPeriodic(a_idim);
94 
95  Array4<EBCellFlag> const& aux_flag = m_cellflags->array(mfi);
96  Array4<Real> const& aux_vfrac = m_volfrac->array(mfi);
97  Array4<Real> const& aux_afrac_x = m_areafrac[0]->array(mfi);
98  Array4<Real> const& aux_afrac_y = m_areafrac[1]->array(mfi);
99  Array4<Real> const& aux_afrac_z = m_areafrac[2]->array(mfi);
100 
101  if (FlagFab[mfi].getType(bx) == FabType::covered ) {
102 
103  ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
104  {
105  aux_flag(i,j,k).setCovered();
106  aux_flag(i,j,k).setDisconnected();
107  if (i==bx.bigEnd(0)) {
108  aux_flag(i+1,j,k).setCovered();
109  }
110  if (j==bx.bigEnd(1)) {
111  aux_flag(i,j+1,k).setCovered();
112  }
113  if (k==bx.bigEnd(2)) {
114  aux_flag(i,j,k+1).setCovered();
115  }
116  });
117 
118  } else if (FlagFab[mfi].getType(bx) == FabType::regular ) {
119 
120  ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
121  {
122  aux_flag(i,j,k).setRegular();
123  aux_flag(i,j,k).setDisconnected();
124  aux_vfrac(i,j,k) = 1.0;
125  aux_afrac_x(i,j,k) = 1.0;
126  aux_afrac_y(i,j,k) = 1.0;
127  aux_afrac_z(i,j,k) = 1.0;
128  if (i==bx.bigEnd(0)) {
129  aux_flag(i+1,j,k).setRegular();
130  aux_vfrac(i+1,j,k) = 1.0;
131  aux_afrac_x(i+1,j,k) = 1.0;
132  }
133  if (j==bx.bigEnd(1)) {
134  aux_flag(i,j+1,k).setRegular();
135  aux_vfrac(i,j+1,k) = 1.0;
136  aux_afrac_y(i,j+1,k) = 1.0;
137  }
138  if (k==bx.bigEnd(2)) {
139  aux_flag(i,j,k+1).setRegular();
140  aux_vfrac(i,j,k+1) = 1.0;
141  aux_afrac_z(i,j,k+1) = 1.0;
142  }
143  });
144 
145  } else if (FlagFab[mfi].getType(bx) == FabType::singlevalued ) {
146 
147  // Initialization
148 
149  // CC cell quantities
150  Array4<EBCellFlag const> const& flag = FlagFab.const_array(mfi);
151  Array4<Real const> const& afrac = (a_factory->getAreaFrac()[a_idim])->const_array(mfi);
152  Array4<Real const> const& bnorm = a_factory->getBndryNormal()[mfi].const_array();
153  Array4<Real const> const& bcent = a_factory->getBndryCent()[mfi].const_array();
154 
155  // aux quantities
156  Array4<Real> const& aux_vcent = m_volcent->array(mfi);
157  Array4<Real> const& aux_fcent_x = m_facecent[0]->array(mfi);
158  Array4<Real> const& aux_fcent_y = m_facecent[1]->array(mfi);
159  Array4<Real> const& aux_fcent_z = m_facecent[2]->array(mfi);
160  Array4<Real> const& aux_barea = m_bndryarea->array(mfi);
161  Array4<Real> const& aux_bcent = m_bndrycent->array(mfi);
162  Array4<Real> const& aux_bnorm = m_bndrynorm->array(mfi);
163 
164  // Extended domain in the direction of periodicity
165  Box dom_grown = domain;
166  for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
167  if (a_geom.isPeriodic(idim)) {
168  dom_grown.grow(idim, a_ngrow[0]);
169  }
170  }
171 
172  const IntVect dom_grown_lo = dom_grown.smallEnd();
173  const IntVect dom_grown_hi = dom_grown.bigEnd();
174 
175  BoxList diffList = boxDiff(bx_grown, bx);
176  for (const Box& b : diffList) {
177  ParallelFor(b, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
178  {
179  if ( i < dom_grown_lo[0] || i > dom_grown_hi[0] ||
180  j < dom_grown_lo[1] || j > dom_grown_hi[1] ||
181  k < dom_grown_lo[2] || k > dom_grown_hi[2] ) {
182  aux_flag(i,j,k).setCovered();
183  aux_flag(i,j,k).setDisconnected();
184  }
185  });
186  }
187 
188 #ifndef AMREX_USE_GPU
189  int const verbose=m_verbose;
190 #endif
191 
192  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
193  {
194  // defaults to covered and disconnected.
195  aux_flag(i,j,k).setCovered();
196  aux_flag(i,j,k).setDisconnected();
197 
198  if (i==bx.bigEnd(0)) {
199  aux_flag(i+1,j,k).setCovered();
200  }
201  if (j==bx.bigEnd(1)) {
202  aux_flag(i,j+1,k).setCovered();
203  }
204  if (k==bx.bigEnd(2)) {
205  aux_flag(i,j,k+1).setCovered();
206  }
207 
208  // Index for low and hi cells
209  IntVect iv_hi(i,j,k);
210  IntVect iv_lo(iv_hi - vdim);
211 
212  bool lo_isCovered = flag(iv_lo).isCovered();
213  bool hi_isCovered = flag(iv_hi).isCovered();
214  bool lo_isRegular = flag(iv_lo).isRegular();
215  bool hi_isRegular = flag(iv_hi).isRegular();
216  bool lo_isSingleValued = flag(iv_lo).isSingleValued();
217  bool hi_isSingleValued = flag(iv_hi).isSingleValued();
218 
219  const bool at_lo_boundary = (!l_periodic && iv_hi[a_idim]==domain.smallEnd(a_idim));
220  const bool at_hi_boundary = (!l_periodic && iv_hi[a_idim]==domain.bigEnd(a_idim));
221 
222  // Treatment of lower boundary
223 
224  if (at_lo_boundary) {
225  if (hi_isCovered) {
226  lo_isCovered = true;
227  lo_isRegular = false;
228  lo_isSingleValued = false;
229  } else if (hi_isRegular) {
230  lo_isCovered = false;
231  lo_isRegular = true;
232  lo_isSingleValued = false;
233  } else if (hi_isSingleValued) {
234  if (almostEqual(afrac(i,j,k),0.0)) {
235  lo_isCovered = true;
236  lo_isRegular = false;
237  lo_isSingleValued = false;
238  } else if (almostEqual(afrac(i,j,k),1.0)) {
239  lo_isCovered = false;
240  lo_isRegular = true;
241  lo_isSingleValued = false;
242  } else {
243  lo_isCovered = false;
244  lo_isRegular = false;
245  lo_isSingleValued = true;
246  iv_lo = iv_hi; // At the lower boundary, low cell takes the values of the high cell.
247  }
248  }
249  }
250 
251  // Treatment of upper boundary
252 
253  if (at_hi_boundary) {
254  if (lo_isCovered) { // Covered
255  hi_isCovered = true;
256  hi_isRegular = false;
257  hi_isSingleValued = false;
258  } else if (lo_isRegular) { // Regular
259  hi_isCovered = false;
260  hi_isRegular = true;
261  hi_isSingleValued = false;
262  } else if (lo_isSingleValued) { // SingleValued
263  if (almostEqual(afrac(i,j,k),0.0)) { //Covered
264  hi_isCovered = true;
265  hi_isRegular = false;
266  hi_isSingleValued = false;
267  } else if (almostEqual(afrac(i,j,k),1.0)) { //Regular
268  hi_isCovered = false;
269  hi_isRegular = true;
270  hi_isSingleValued = false;
271  } else { // SingleValued
272  hi_isCovered = false;
273  hi_isRegular = false;
274  hi_isSingleValued = true;
275  iv_hi = iv_lo; // At the upper boundary, hi cell takes the values of the low cell.
276  }
277  }
278  }
279 
280  if ( lo_isCovered && hi_isCovered) {
281 
282  // defaults to covered and disconnected.
283 
284  } else if ( lo_isRegular && hi_isRegular) {
285 
286  aux_flag(i,j,k).setRegular();
287  aux_flag(i,j,k).setConnected();
288 
289  aux_vfrac(i,j,k) = 1.0;
290 
291  aux_afrac_x(i,j,k) = 1.0;
292  aux_afrac_y(i,j,k) = 1.0;
293  aux_afrac_z(i,j,k) = 1.0;
294 
295  if (i==bx.bigEnd(0)) {
296  aux_afrac_x(i+1,j,k) = 1.0;
297  }
298  if (j==bx.bigEnd(1)) {
299  aux_afrac_y(i,j+1,k) = 1.0;
300  }
301  if (k==bx.bigEnd(2)) {
302  aux_afrac_z(i,j,k+1) = 1.0;
303  }
304 
305  } else {
306 
307 #ifndef AMREX_USE_GPU
308  if (verbose) { Print() << "\ncell: " << amrex::IntVect(i,j,k) << "\n"; }
309 #endif
310  Array<Real,AMREX_SPACEDIM> lo_arr = {-0.5,-0.5,-0.5};
311  Array<Real,AMREX_SPACEDIM> hi_arr = { 0.5, 0.5, 0.5};
312 
313  //-----------------------
314  // Low EB cut cell
315  //-----------------------
316 
317  // Map bcent and bnorm to the isoparametric space for anisotropic grids.
318  // (This step is needed because bcent in AMReX is isotropically normalized.)
319 
320  RealVect lo_point (bcent(iv_lo,0), bcent(iv_lo,1), bcent(iv_lo,2));
321  RealVect lo_normal(bnorm(iv_lo,0), bnorm(iv_lo,1), bnorm(iv_lo,2));
322 
323  if (at_lo_boundary) { // At lower boundary
324  lo_point[a_idim] += 1.0; // Move the boundary centroid upward in the a_idim direction.
325  }
326 
327  if (lo_isSingleValued ) {
328  Real bnorm_x = bnorm(iv_lo,0) * dx[0];
329  Real bnorm_y = bnorm(iv_lo,1) * dx[1];
330  Real bnorm_z = bnorm(iv_lo,2) * dx[2];
331 
332  Real norm = sqrt( bnorm_x*bnorm_x + bnorm_y*bnorm_y + bnorm_z*bnorm_z);
333 
334  RealVect bnorm_isoparam ( bnorm_x / norm, bnorm_y / norm, bnorm_z / norm);
335 
336  lo_normal = bnorm_isoparam;
337  }
338 
339  // High side of low cell
340  lo_arr[a_idim] = 0.0;
341  hi_arr[a_idim] = 0.5;
342  RealBox lo_rbx(lo_arr.data(), hi_arr.data());
343 
344  eb_cut_cell_ lo_eb_cc(flag(iv_lo), lo_rbx, lo_point, lo_normal);
345 
346  // cell iv_lo covered (regular) imples lo_eb_cc is covered (regular)
347  // The inverse is not always true.
348  AMREX_ASSERT( !lo_isCovered || lo_eb_cc.isCovered() );
349  AMREX_ASSERT( !lo_isRegular || lo_eb_cc.isRegular() );
350 
351  //-----------------------
352  // High EB cut cell
353  //-----------------------
354 
355  RealVect hi_point (bcent(iv_hi,0), bcent(iv_hi,1), bcent(iv_hi,2));
356  RealVect hi_normal(bnorm(iv_hi,0), bnorm(iv_hi,1), bnorm(iv_hi,2));
357 
358  if (at_hi_boundary) {
359  hi_point[a_idim] += -1.0; // Move the boundary centroid downward in the a_idim direction.
360  }
361 
362  if (hi_isSingleValued ) {
363  Real bnorm_x = bnorm(iv_hi,0) * dx[0];
364  Real bnorm_y = bnorm(iv_hi,1) * dx[1];
365  Real bnorm_z = bnorm(iv_hi,2) * dx[2];
366 
367  Real norm = sqrt( bnorm_x*bnorm_x + bnorm_y*bnorm_y + bnorm_z*bnorm_z);
368 
369  RealVect bnorm_isoparam ( bnorm_x / norm, bnorm_y / norm, bnorm_z / norm);
370 
371  hi_normal = bnorm_isoparam;
372  }
373 
374  // Low side of high cell
375  lo_arr[a_idim] = -0.5;
376  hi_arr[a_idim] = 0.0;
377  RealBox hi_rbx(lo_arr.data(), hi_arr.data());
378 
379  eb_cut_cell_ hi_eb_cc(flag(iv_hi), hi_rbx, hi_point, hi_normal);
380 
381  // cell iv_hi covered (regular) imples hi_eb_cc is covered (regular)
382  // The inverse is not always true.
383  AMREX_ASSERT( !hi_isCovered || hi_eb_cc.isCovered() );
384  AMREX_ASSERT( !hi_isRegular || hi_eb_cc.isRegular() );
385 
386 #if 0
387 #if defined(AMREX_DEBUG) || defined(AMREX_TESTING) || 1
388 
389  { /***************************** SANITY CHECK ***********************\
390  * Perform some basic sanity checks to verify that what we computed *
391  * for cell (i,j,k) compares to what we know to be true. *
392  \******************************************************************/
393 
394  // Compute the cut-cell for the high side of the high cell. This is
395  // only needed for sanity checks.
396 
397  eb_cut_cell_ hi_hi_eb_cc(flag(iv_hi), lo_rbx, hi_point, hi_normal);
398 
399  // cell iv_hi covered (regular) imples hi_hi_eb_cc is covered (regular)
400  // The inverse is not always true.
401 #ifndef AMREX_USE_GPU
402  if ( !(!hi_isRegular || hi_hi_eb_cc.isRegular()) ||
403  !(!hi_isCovered || hi_hi_eb_cc.isCovered()) ) {
404  Print() << "flag(iv_hi) and hi_hi_eb_cc flags do not agree\n"
405  << "\n isRegular() " << hi_isRegular << " " << hi_hi_eb_cc.isRegular()
406  << "\n isCovered() " << hi_isCovered << " " << hi_hi_eb_cc.isCovered()
407  << "\n";
408  }
409 #endif
410  // If cell iv_hi is regular or covered, then hi_hi_eb_cc must also
411  // be regular or covered. The inverse is not true.
412  AMREX_ALWAYS_ASSERT( !hi_isRegular || hi_hi_eb_cc.isRegular() );
413  AMREX_ALWAYS_ASSERT( !hi_isCovered || hi_hi_eb_cc.isCovered() );
414 
415  // The area and volume fractions that are computed for the scalar grid
416  // are slightly different than those we compute from the geometric
417  // reconstruction using the EB point and normal. However, we expect
418  // that the area fractions computed here will give back the same
419  // normal we used to compute them.
420  if ( hi_isSingleValued ) {
421 
422  Real const adx = (a_idim == 0)
423  ? (hi_eb_cc.areaLo(0) - hi_hi_eb_cc.areaHi(0)) * dx[1] * dx[2]
424  : (hi_eb_cc.areaLo(0) + hi_hi_eb_cc.areaLo(0)) * dx[1] * dx[2]
425  - (hi_eb_cc.areaHi(0) + hi_hi_eb_cc.areaHi(0)) * dx[1] * dx[2];
426 
427  Real const ady = (a_idim == 1)
428  ? (hi_eb_cc.areaLo(1) - hi_hi_eb_cc.areaHi(1)) * dx[0] * dx[2]
429  : (hi_eb_cc.areaLo(1) + hi_hi_eb_cc.areaLo(1)) * dx[0] * dx[2]
430  - (hi_eb_cc.areaHi(1) + hi_hi_eb_cc.areaHi(1)) * dx[0] * dx[2];
431 
432  Real const adz = (a_idim == 2)
433  ? (hi_eb_cc.areaLo(2) - hi_hi_eb_cc.areaHi(2)) * dx[0] * dx[1]
434  : (hi_eb_cc.areaLo(2) + hi_hi_eb_cc.areaLo(2)) * dx[0] * dx[1]
435  - (hi_eb_cc.areaHi(2) + hi_hi_eb_cc.areaHi(2)) * dx[0] * dx[1];
436 
437  Real const apnorm = std::sqrt(adx*adx + ady*ady + adz*adz);
438 
439  // EB normal
440  Real const apnorminv = 1. / apnorm;
441  RealVect const normal(adx*apnorminv, ady*apnorminv, adz*apnorminv);
442  Real const dot_normals = normal.dotProduct(hi_normal);
443 
444 #ifndef AMREX_USE_GPU
445  if ( !amrex::almostEqual(dot_normals, 1.0) ) {
446  Print() << "\nFail: check-1 dot_normals " << dot_normals
447  << '\n';
448 
449  hi_eb_cc.debug();
450  hi_hi_eb_cc.debug();
451 
452  } else if (verbose) {
453  Print() << "Pass: dot_normals = 1.0\n";
454 
455  }
456 #endif
457  AMREX_ALWAYS_ASSERT( amrex::almostEqual(dot_normals, 1.0) );
458  }
459 
460  // The a_idim area of hi_eb_cc.areaHi() should equal hi_hi_eb_cc.areaLo()
461  {
462 #ifndef AMREX_USE_GPU
463  Real const abs_err = std::abs( hi_eb_cc.areaHi(a_idim) - hi_hi_eb_cc.areaLo(a_idim) );
464  Real machine_tol = 10.0*std::numeric_limits<amrex::Real>::epsilon();
465  if ( abs_err >= machine_tol ) {
466  Print() << "\nFail: check-2 area abs_err: " << abs_err
467  << "\n hi_eb_cc.areaHi " << hi_eb_cc.areaHi(a_idim)
468  << "\n hi_hi_eb_cc.areaLo " << hi_hi_eb_cc.areaLo(a_idim)
469  << '\n';
470  } else if (verbose) {
471  Print() << "Pass: hi_eb_cc.areaHi = hi_hi_eb_cc.areaLo"
472  << " abs_err: " << abs_err << "\n";
473  }
474  AMREX_ALWAYS_ASSERT( abs_err < machine_tol );
475 #endif
476  }
477 
478  // The low-side area of hi_eb_cc should equal a_idim afrac.
479  { Real const abs_err = amrex::max(std::abs(lo_eb_cc.areaHi(a_idim) - afrac(iv_hi)),
480  std::abs(hi_eb_cc.areaLo(a_idim) - afrac(iv_hi)));
481  Real compare_tol = 5.0e-6;
482 #ifndef AMREX_USE_GPU
483  if ( abs_err >= compare_tol ) {
484  //hi_eb_cc.debug();
485  Print() << "\nFail: check-3 area abs_err " << abs_err
486  << "\n hi_eb_cc.areaLo(" << a_idim << ") = " << hi_eb_cc.areaLo(a_idim)
487  << "\n lo_eb_cc.areaHi(" << a_idim << ") = " << lo_eb_cc.areaHi(a_idim)
488  << "\n afrac" << iv_hi << " = " << afrac(iv_hi)
489  << '\n';
490  } else if (verbose) {
491  Print() << "Pass: hi_eb_cc.areaLo = afrac = " << afrac(iv_hi)
492  << " abs_err: " << abs_err << "\n";
493  }
494 #endif
495  AMREX_ALWAYS_ASSERT( abs_err < compare_tol );
496  }
497 
498  // The combined volumes of hi_eb_cc.areaHi() and hi_hi_eb_cc should
499  // equal vfrac(iv_hi).
500  { Real const vol = hi_eb_cc.volume() + hi_hi_eb_cc.volume();
501  Real const abs_err = amrex::Math::abs(vfrac(iv_hi) - vol);
502  Real compare_tol = 5.0e-6;
503 #ifndef AMREX_USE_GPU
504  if ( abs_err >= compare_tol ) {
505  hi_eb_cc.debug();
506  hi_hi_eb_cc.debug();
507  amrex::Print() << "\nFail: check-4 volume abs_err: " << abs_err
508  << "\n point: " << hi_point
509  << "\n normal: " << hi_normal
510  << "\n hi_eb_cc.volume() " << hi_eb_cc.volume()
511  << "\n hi_hi_eb_cc.volume() " << hi_hi_eb_cc.volume()
512  << "\n vfrac: " << vfrac(iv_hi)
513  << '\n';
514  } else if (verbose) {
515  Print() << "Pass: hi_eb_cc + hi_hi_eb_cc = vfrac = " << vfrac(iv_hi)
516  << " abs_err: " << abs_err << "\n";
517  }
518 #endif
519  AMREX_ALWAYS_ASSERT( abs_err < compare_tol );
520  }
521  } //
522 #endif
523 #endif // 0
524 
525  //-----------------------
526  // Fill out aux_ arrays
527  //-----------------------
528 
529  if (lo_eb_cc.isCovered() && hi_eb_cc.isCovered()) {
530 
531  // defaults to covered and disconnected.
532 
533  } else if (lo_eb_cc.isRegular() && hi_eb_cc.isRegular()) {
534 
535  aux_flag(i,j,k).setRegular();
536  aux_flag(i,j,k).setConnected();
537 
538  aux_vfrac(i,j,k) = 1.0;
539 
540  aux_afrac_x(i,j,k) = 1.0;
541  aux_afrac_y(i,j,k) = 1.0;
542  aux_afrac_z(i,j,k) = 1.0;
543 
544  aux_fcent_x(i,j,k,0) = 0.0; aux_fcent_x(i,j,k,1) = 0.0;
545  aux_fcent_y(i,j,k,0) = 0.0; aux_fcent_y(i,j,k,1) = 0.0;
546  aux_fcent_z(i,j,k,0) = 0.0; aux_fcent_z(i,j,k,1) = 0.0;
547 
548  if (i==bx.bigEnd(0)) {
549  aux_afrac_x(i+1,j,k) = 1.0;
550  aux_fcent_x(i+1,j,k,0) = 0.0; aux_fcent_x(i+1,j,k,1) = 0.0;
551  }
552  if (j==bx.bigEnd(1)) {
553  aux_afrac_y(i,j+1,k) = 1.0;
554  aux_fcent_y(i,j+1,k,0) = 0.0; aux_fcent_y(i,j+1,k,1) = 0.0;
555  }
556  if (k==bx.bigEnd(2)) {
557  aux_afrac_z(i,j,k+1) = 1.0;
558  aux_fcent_z(i,j,k+1,0) = 0.0; aux_fcent_z(i,j,k+1,1) = 0.0;
559  }
560 
561  } else if ( (lo_eb_cc.isRegular() && hi_eb_cc.isCovered())
562  || (lo_eb_cc.isCovered() && hi_eb_cc.isRegular()) ) {
563 
564  // This is a problematic situation.
565 #ifndef AMREX_USE_GPU
566  Print()<< "eb_aux_ / Check: Regular and Covered cut cells are facing each other." << std::endl;
567 #endif
568 
569  } else {
570 
571  // 0. Cell Flag
572 
573  aux_flag(i,j,k).setSingleValued();
574 
575  // 1. Volume Fraction
576 
577  Real lo_vol {lo_eb_cc.volume()};
578  Real hi_vol {hi_eb_cc.volume()};
579 
580  aux_vfrac(i,j,k) = lo_vol + hi_vol;
581 
582  // 2. Volume Centroid
583 
584  /* centVol() returns the coordinates based on m_rbx.
585  The coordinates in the a_idim direction are in [0.0,0.5] for the low cell and in [-0.5,0.0] for the hi cell.
586  Therefore, they need to be mapped to the eb_aux space, by shifting:
587  x' = x - 0.5 (low cell), x + 0.5 (hi cell) if a_idim = 0
588  y' = y - 0.5 (low cell), y + 0.5 (hi cell) if a_idim = 1
589  z' = z - 0.5 (low cell), z + 0.5 (hi cell) if a_idim = 2
590  */
591 
592  RealVect lo_vcent {lo_eb_cc.centVol()};
593  RealVect hi_vcent {hi_eb_cc.centVol()};
594 
595  lo_vcent[a_idim] = lo_vcent[a_idim] - 0.5;
596  hi_vcent[a_idim] = hi_vcent[a_idim] + 0.5;
597 
598  aux_vcent(i,j,k,0) = ( lo_vol * lo_vcent[0] + hi_vol * hi_vcent[0] ) / aux_vfrac(i,j,k);
599  aux_vcent(i,j,k,1) = ( lo_vol * lo_vcent[1] + hi_vol * hi_vcent[1] ) / aux_vfrac(i,j,k);
600  aux_vcent(i,j,k,2) = ( lo_vol * lo_vcent[2] + hi_vol * hi_vcent[2] ) / aux_vfrac(i,j,k);
601 
602  // 3. Area Fraction
603 
604  Real lo_areaLo_x {lo_eb_cc.areaLo(0)};
605  Real lo_areaLo_y {lo_eb_cc.areaLo(1)};
606  Real lo_areaLo_z {lo_eb_cc.areaLo(2)};
607 
608  Real hi_areaLo_x {hi_eb_cc.areaLo(0)};
609  Real hi_areaLo_y {hi_eb_cc.areaLo(1)};
610  Real hi_areaLo_z {hi_eb_cc.areaLo(2)};
611 
612  aux_afrac_x(i,j,k) = (a_idim == 0) ? lo_areaLo_x : lo_areaLo_x + hi_areaLo_x;
613  aux_afrac_y(i,j,k) = (a_idim == 1) ? lo_areaLo_y : lo_areaLo_y + hi_areaLo_y;
614  aux_afrac_z(i,j,k) = (a_idim == 2) ? lo_areaLo_z : lo_areaLo_z + hi_areaLo_z;
615 
616  if (i==bx.bigEnd(0)) {
617  Real lo_areaHi_x {lo_eb_cc.areaHi(0)};
618  Real hi_areaHi_x {hi_eb_cc.areaHi(0)};
619  aux_afrac_x(i+1,j,k) = (a_idim == 0) ? hi_areaHi_x : lo_areaHi_x + hi_areaHi_x;
620  }
621  if (j==bx.bigEnd(1)) {
622  Real lo_areaHi_y {lo_eb_cc.areaHi(1)};
623  Real hi_areaHi_y {hi_eb_cc.areaHi(1)};
624  aux_afrac_y(i,j+1,k) = (a_idim == 1) ? hi_areaHi_y : lo_areaHi_y + hi_areaHi_y;
625  }
626  if (k==bx.bigEnd(2)) {
627  Real lo_areaHi_z {lo_eb_cc.areaHi(2)};
628  Real hi_areaHi_z {hi_eb_cc.areaHi(2)};
629  aux_afrac_z(i,j,k+1) = (a_idim == 2) ? hi_areaHi_z : lo_areaHi_z + hi_areaHi_z;
630  }
631 
632  // 4. Face Centroid
633 
634  /* fcentLo returns the coordinates based on m_rbx.
635  The coordinates in the a_idim direction are in [0.0,0.5] for the low cell and in [-0.5,0.0] for the hi cell.
636  Therefore, they need to be mapped to the eb_aux space, by shifting:
637  x' = x - 0.5 (low cell), x + 0.5 (hi cell) if a_idim = 0
638  y' = y - 0.5 (low cell), y + 0.5 (hi cell) if a_idim = 1
639  z' = z - 0.5 (low cell), z + 0.5 (hi cell) if a_idim = 2
640  */
641 
642  RealVect lo_centLo_x {lo_eb_cc.centLo(0)};
643  RealVect lo_centLo_y {lo_eb_cc.centLo(1)};
644  RealVect lo_centLo_z {lo_eb_cc.centLo(2)};
645 
646  RealVect hi_centLo_x {hi_eb_cc.centLo(0)};
647  RealVect hi_centLo_y {hi_eb_cc.centLo(1)};
648  RealVect hi_centLo_z {hi_eb_cc.centLo(2)};
649 
650  if (a_idim == 0) {
651  aux_fcent_x(i,j,k,0) = lo_centLo_x[1]; // y
652  aux_fcent_x(i,j,k,1) = lo_centLo_x[2]; // z
653  aux_fcent_y(i,j,k,0) = (aux_afrac_y(i,j,k) > 0.0) // x (mapped)
654  ? ( lo_areaLo_y * (lo_centLo_y[0] - 0.5)
655  + hi_areaLo_y * (hi_centLo_y[0] + 0.5) ) / aux_afrac_y(i,j,k)
656  : 0.0;
657  aux_fcent_y(i,j,k,1) = (aux_afrac_y(i,j,k) > 0.0) // z
658  ? ( lo_areaLo_y * lo_centLo_y[2]
659  + hi_areaLo_y * hi_centLo_y[2] ) / aux_afrac_y(i,j,k)
660  : 0.0;
661  aux_fcent_z(i,j,k,0) = (aux_afrac_z(i,j,k) > 0.0) // x (mapped)
662  ? ( lo_areaLo_z * (lo_centLo_z[0] - 0.5)
663  + hi_areaLo_z * (hi_centLo_z[0] + 0.5) ) / aux_afrac_z(i,j,k)
664  : 0.0;
665  aux_fcent_z(i,j,k,1) = (aux_afrac_z(i,j,k) > 0.0) // y
666  ? ( lo_areaLo_z * lo_centLo_z[1]
667  + hi_areaLo_z * hi_centLo_z[1] ) / aux_afrac_z(i,j,k)
668  : 0.0;
669  } else if (a_idim == 1) {
670  aux_fcent_x(i,j,k,0) = (aux_afrac_x(i,j,k) > 0.0) // y (mapped)
671  ? ( lo_areaLo_x * (lo_centLo_x[1] - 0.5)
672  + hi_areaLo_x * (hi_centLo_x[1] + 0.5) ) / aux_afrac_x(i,j,k)
673  : 0.0;
674  aux_fcent_x(i,j,k,1) = (aux_afrac_x(i,j,k) > 0.0) // z
675  ? ( lo_areaLo_x * lo_centLo_x[2]
676  + hi_areaLo_x * hi_centLo_x[2] ) / aux_afrac_x(i,j,k)
677  : 0.0;
678  aux_fcent_y(i,j,k,0) = lo_centLo_y[0]; // x
679  aux_fcent_y(i,j,k,1) = lo_centLo_y[2]; // z
680  aux_fcent_z(i,j,k,0) = (aux_afrac_z(i,j,k) > 0.0) // x
681  ? ( lo_areaLo_z * lo_centLo_z[0]
682  + hi_areaLo_z * hi_centLo_z[0] ) / aux_afrac_z(i,j,k)
683  : 0.0;
684  aux_fcent_z(i,j,k,1) = (aux_afrac_z(i,j,k) > 0.0) // y (mapped)
685  ? ( lo_areaLo_z * (lo_centLo_z[1] - 0.5)
686  + hi_areaLo_z * (hi_centLo_z[1] + 0.5) ) / aux_afrac_z(i,j,k)
687  : 0.0;
688  } else if (a_idim == 2) {
689  aux_fcent_x(i,j,k,0) = (aux_afrac_x(i,j,k) > 0.0) // y
690  ? ( lo_areaLo_x * lo_centLo_x[1]
691  + hi_areaLo_x * hi_centLo_x[1] ) / aux_afrac_x(i,j,k)
692  : 0.0;
693  aux_fcent_x(i,j,k,1) = (aux_afrac_x(i,j,k) > 0.0) // z (mapped)
694  ? ( lo_areaLo_x * (lo_centLo_x[2] - 0.5)
695  + hi_areaLo_x * (hi_centLo_x[2] + 0.5) ) / aux_afrac_x(i,j,k)
696  : 0.0;
697  aux_fcent_y(i,j,k,0) = (aux_afrac_y(i,j,k) > 0.0) // x
698  ? ( lo_areaLo_y * lo_centLo_y[0]
699  + hi_areaLo_y * hi_centLo_y[0] ) / aux_afrac_y(i,j,k)
700  : 0.0;
701  aux_fcent_y(i,j,k,1) = (aux_afrac_y(i,j,k) > 0.0) // z (mapped)
702  ? ( lo_areaLo_y * (lo_centLo_y[2] - 0.5)
703  + hi_areaLo_y * (hi_centLo_y[2] + 0.5) ) / aux_afrac_y(i,j,k)
704  : 0.0;
705  aux_fcent_z(i,j,k,0) = lo_centLo_z[0]; // x
706  aux_fcent_z(i,j,k,1) = lo_centLo_z[1]; // y
707  }
708 
709  if (i==bx.bigEnd(0)) {
710  Real lo_areaHi_x {lo_eb_cc.areaHi(0)};
711  Real hi_areaHi_x {hi_eb_cc.areaHi(0)};
712  RealVect lo_centHi_x {lo_eb_cc.centHi(0)};
713  RealVect hi_centHi_x {hi_eb_cc.centHi(0)};
714  if (a_idim == 0) {
715  aux_fcent_x(i+1,j,k,0) = hi_centHi_x[1]; // y
716  aux_fcent_x(i+1,j,k,1) = hi_centHi_x[2]; // z
717  } else if (a_idim == 1) {
718  aux_fcent_x(i+1,j,k,0) = (aux_afrac_x(i+1,j,k) > 0.0) // y (mapped)
719  ? ( lo_areaHi_x * (lo_centHi_x[1] - 0.5)
720  + hi_areaHi_x * (hi_centHi_x[1] + 0.5) ) / aux_afrac_x(i+1,j,k)
721  : 0.0;
722  aux_fcent_x(i+1,j,k,1) = (aux_afrac_x(i+1,j,k) > 0.0) // z
723  ? ( lo_areaHi_x * lo_centHi_x[2]
724  + hi_areaHi_x * hi_centHi_x[2] ) / aux_afrac_x(i+1,j,k)
725  : 0.0;
726  } else if (a_idim == 2) {
727  aux_fcent_x(i+1,j,k,0) = (aux_afrac_x(i+1,j,k) > 0.0) // y
728  ? ( lo_areaHi_x * lo_centHi_x[1]
729  + hi_areaHi_x * hi_centHi_x[1] ) / aux_afrac_x(i+1,j,k)
730  : 0.0;
731  aux_fcent_x(i+1,j,k,1) = (aux_afrac_x(i+1,j,k) > 0.0) // z (mapped)
732  ? ( lo_areaHi_x * (lo_centHi_x[2] - 0.5)
733  + hi_areaHi_x * (hi_centHi_x[2] + 0.5) ) / aux_afrac_x(i+1,j,k)
734  : 0.0;
735  }
736  }
737  if (j==bx.bigEnd(1)) {
738  Real lo_areaHi_y {lo_eb_cc.areaHi(1)};
739  Real hi_areaHi_y {hi_eb_cc.areaHi(1)};
740  RealVect lo_centHi_y {lo_eb_cc.centHi(1)};
741  RealVect hi_centHi_y {hi_eb_cc.centHi(1)};
742  if (a_idim == 0) {
743  aux_fcent_y(i,j+1,k,0) = (aux_afrac_y(i,j+1,k) > 0.0) // x (mapped)
744  ? ( lo_areaHi_y * (lo_centHi_y[0] - 0.5)
745  + hi_areaHi_y * (hi_centHi_y[0] + 0.5) ) / aux_afrac_y(i,j+1,k)
746  : 0.0;
747  aux_fcent_y(i,j+1,k,1) = (aux_afrac_y(i,j+1,k) > 0.0) // z
748  ? ( lo_areaHi_y * lo_centHi_y[2]
749  + hi_areaHi_y * hi_centHi_y[2] ) / aux_afrac_y(i,j+1,k)
750  : 0.0;
751  } else if (a_idim == 1) {
752  aux_fcent_y(i,j+1,k,0) = lo_centHi_y[0]; // x
753  aux_fcent_y(i,j+1,k,1) = lo_centHi_y[2]; // z
754  } else if (a_idim == 2) {
755  aux_fcent_y(i,j+1,k,0) = (aux_afrac_y(i,j+1,k) > 0.0) // x
756  ? ( lo_areaHi_y * lo_centHi_y[0]
757  + hi_areaHi_y * hi_centHi_y[0] ) / aux_afrac_y(i,j+1,k)
758  : 0.0;
759  aux_fcent_y(i,j+1,k,1) = (aux_afrac_y(i,j+1,k) > 0.0) // z (mapped)
760  ? ( lo_areaHi_y * (lo_centHi_y[2] - 0.5)
761  + hi_areaHi_y * (hi_centHi_y[2] + 0.5) ) / aux_afrac_y(i,j+1,k)
762  : 0.0;
763  }
764  }
765  if (k==bx.bigEnd(2)) {
766  Real lo_areaHi_z {lo_eb_cc.areaHi(2)};
767  Real hi_areaHi_z {hi_eb_cc.areaHi(2)};
768  RealVect lo_centHi_z {lo_eb_cc.centHi(2)};
769  RealVect hi_centHi_z {hi_eb_cc.centHi(2)};
770  if (a_idim == 0) {
771  aux_fcent_z(i,j,k+1,0) = (aux_afrac_z(i,j,k+1) > 0.0) // x (mapped)
772  ? ( lo_areaHi_z * (lo_centHi_z[0] - 0.5)
773  + hi_areaHi_z * (hi_centHi_z[0] + 0.5) ) / aux_afrac_z(i,j,k+1)
774  : 0.0;
775  aux_fcent_z(i,j,k+1,1) = (aux_afrac_z(i,j,k+1) > 0.0) // y
776  ? ( lo_areaHi_z * lo_centHi_z[1]
777  + hi_areaHi_z * hi_centHi_z[1] ) / aux_afrac_z(i,j,k+1)
778  : 0.0;
779  } else if (a_idim == 1) {
780  aux_fcent_z(i,j,k+1,0) = (aux_afrac_z(i,j,k+1) > 0.0) // x
781  ? ( lo_areaHi_z * lo_centHi_z[0]
782  + hi_areaHi_z * hi_centHi_z[0] ) / aux_afrac_z(i,j,k+1)
783  : 0.0;
784  aux_fcent_z(i,j,k+1,1) = (aux_afrac_z(i,j,k+1) > 0.0) // y (mapped)
785  ? ( lo_areaHi_z * (lo_centHi_z[1] - 0.5)
786  + hi_areaHi_z * (hi_centHi_z[1] + 0.5) ) / aux_afrac_z(i,j,k+1)
787  : 0.0;
788  } else if (a_idim == 2) {
789  aux_fcent_z(i,j,k+1,0) = lo_centHi_z[0]; // x
790  aux_fcent_z(i,j,k+1,1) = lo_centHi_z[1]; // y
791  }
792  }
793 
794  // 5. Boundary Area
795 
796  Real lo_areaBoun {lo_eb_cc.areaBoun()};
797  Real hi_areaBoun {hi_eb_cc.areaBoun()};
798 
799  aux_barea(i,j,k) = lo_areaBoun + hi_areaBoun;
800 
801  // 6. Boundary Centroid
802 
803  RealVect lo_centBoun {lo_eb_cc.centBoun()};
804  RealVect hi_centBoun {hi_eb_cc.centBoun()};
805 
806  if (a_idim == 0) {
807  aux_bcent(i,j,k,0) = ( lo_areaBoun * (lo_centBoun[0]-0.5) + hi_areaBoun * (hi_centBoun[0]+0.5) ) / aux_barea(i,j,k); // x (mapped)
808  aux_bcent(i,j,k,1) = ( lo_areaBoun * lo_centBoun[1] + hi_areaBoun * hi_centBoun[1] ) / aux_barea(i,j,k); // y
809  aux_bcent(i,j,k,2) = ( lo_areaBoun * lo_centBoun[2] + hi_areaBoun * hi_centBoun[2] ) / aux_barea(i,j,k); // z
810  } else if (a_idim == 1) {
811  aux_bcent(i,j,k,0) = ( lo_areaBoun * lo_centBoun[0] + hi_areaBoun * hi_centBoun[0] ) / aux_barea(i,j,k); // x
812  aux_bcent(i,j,k,1) = ( lo_areaBoun * (lo_centBoun[1]-0.5) + hi_areaBoun * (hi_centBoun[1]+0.5) ) / aux_barea(i,j,k); // y (mapped)
813  aux_bcent(i,j,k,2) = ( lo_areaBoun * lo_centBoun[2] + hi_areaBoun * hi_centBoun[2] ) / aux_barea(i,j,k); // z
814  } else if (a_idim == 2) {
815  aux_bcent(i,j,k,0) = ( lo_areaBoun * lo_centBoun[0] + hi_areaBoun * hi_centBoun[0] ) / aux_barea(i,j,k); // x
816  aux_bcent(i,j,k,1) = ( lo_areaBoun * lo_centBoun[1] + hi_areaBoun * hi_centBoun[1] ) / aux_barea(i,j,k); // y
817  aux_bcent(i,j,k,2) = ( lo_areaBoun * (lo_centBoun[2]-0.5) + hi_areaBoun * (hi_centBoun[2]+0.5) ) / aux_barea(i,j,k); // z (mapped)
818  }
819 
820  // 7. Boundary Normal
821 
822  RealVect eb_normal = ( lo_areaBoun * lo_normal + hi_areaBoun * hi_normal )/ aux_barea(i,j,k);
823 
824  aux_bnorm(i,j,k,0) = eb_normal[0];
825  aux_bnorm(i,j,k,1) = eb_normal[1];
826  aux_bnorm(i,j,k,2) = eb_normal[2];
827 
828  }
829 
830  } // flag(iv_lo) and flag(iv_hi)
831 
832  });
833 
834  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
835  {
836  if (aux_vfrac(i,j,k) < small_volfrac)
837  {
838  aux_vfrac(i,j,k) = 0.0;
839  }
840  });
841 
842  } // if (FlagFab[mfi].getType(bx) == FabType::singlevalued )
843 
844  } // MFIter
845 
846  // We FillBoundary volfrac here so that we can use tests on volfrac in ghost cells below
847  m_volfrac->FillBoundary(a_geom.periodicity());
848 
849  for (MFIter mfi(*m_cellflags, false); mfi.isValid(); ++mfi) {
850 
851  const Box& bx = mfi.validbox();
852  const Box& bx_grown = mfi.growntilebox();
853 
854  Array4<EBCellFlag> const& aux_flag = m_cellflags->array(mfi);
855  Array4<Real> const& aux_vfrac = m_volfrac->array(mfi);
856  Array4<Real> const& aux_afrac_x = m_areafrac[0]->array(mfi);
857  Array4<Real> const& aux_afrac_y = m_areafrac[1]->array(mfi);
858  Array4<Real> const& aux_afrac_z = m_areafrac[2]->array(mfi);
859 
860  Array4<Real> const& aux_vcent = m_volcent->array(mfi);
861  Array4<Real> const& aux_fcent_x = m_facecent[0]->array(mfi);
862  Array4<Real> const& aux_fcent_y = m_facecent[1]->array(mfi);
863  Array4<Real> const& aux_fcent_z = m_facecent[2]->array(mfi);
864  Array4<Real> const& aux_barea = m_bndryarea->array(mfi);
865  Array4<Real> const& aux_bcent = m_bndrycent->array(mfi);
866  Array4<Real> const& aux_bnorm = m_bndrynorm->array(mfi);
867 
868 
869  if (FlagFab[mfi].getType(bx) == FabType::singlevalued ) {
870 
871  // Corrections for small cells
872  Box my_xbx(bx); if (a_idim == 0) my_xbx.growLo(0,1);
873  ParallelFor(my_xbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
874  {
875  if (aux_vfrac(i,j,k) < small_volfrac) {
876  aux_afrac_x(i ,j ,k ) = 0.0;
877  aux_afrac_x(i+1,j ,k ) = 0.0;
878  }
879  });
880 
881  Box my_ybx(bx); if (a_idim == 1) my_ybx.growLo(1,1);
882  ParallelFor(my_ybx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
883  {
884  if (aux_vfrac(i,j,k) < small_volfrac) {
885  aux_afrac_y(i ,j ,k ) = 0.0;
886  aux_afrac_y(i ,j+1,k ) = 0.0;
887  }
888  });
889 
890  Box my_zbx(bx); if (a_idim == 2) my_zbx.growLo(2,1);
891  ParallelFor(my_zbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
892  {
893  if (aux_vfrac(i,j,k) < small_volfrac) {
894  aux_afrac_z(i ,j ,k+1) = 0.0;
895  aux_afrac_z(i ,j ,k ) = 0.0;
896  }
897  });
898 
899  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
900  {
901  if (aux_vfrac(i,j,k) < small_volfrac)
902  {
903  aux_vcent(i,j,k,0) = 0.0;
904  aux_vcent(i,j,k,1) = 0.0;
905  aux_vcent(i,j,k,2) = 0.0;
906 
907  aux_fcent_x(i ,j ,k ,0) = 0.0;
908  aux_fcent_x(i ,j ,k ,1) = 0.0;
909  aux_fcent_x(i+1,j ,k ,0) = 0.0;
910  aux_fcent_x(i+1,j ,k ,1) = 0.0;
911 
912  aux_fcent_y(i ,j ,k ,0) = 0.0;
913  aux_fcent_y(i ,j ,k ,1) = 0.0;
914  aux_fcent_y(i ,j+1,k ,0) = 0.0;
915  aux_fcent_y(i ,j+1,k ,1) = 0.0;
916 
917  aux_fcent_z(i ,j ,k ,0) = 0.0;
918  aux_fcent_z(i ,j ,k ,1) = 0.0;
919  aux_fcent_z(i ,j ,k+1,0) = 0.0;
920  aux_fcent_z(i ,j ,k+1,1) = 0.0;
921 
922  aux_barea(i,j,k) = 0.0;
923 
924  aux_bcent(i,j,k,0) = 0.0;
925  aux_bcent(i,j,k,1) = 0.0;
926  aux_bcent(i,j,k,2) = 0.0;
927 
928  aux_bnorm(i,j,k,0) = 0.0;
929  aux_bnorm(i,j,k,1) = 0.0;
930  aux_bnorm(i,j,k,2) = 0.0;
931 
932  aux_flag(i,j,k).setCovered();
933  }
934 
935  if (std::abs(aux_vcent(i,j,k,0)) < small_value) aux_vcent(i,j,k,0) = 0.0;
936  if (std::abs(aux_vcent(i,j,k,1)) < small_value) aux_vcent(i,j,k,1) = 0.0;
937  if (std::abs(aux_vcent(i,j,k,2)) < small_value) aux_vcent(i,j,k,2) = 0.0;
938  if (std::abs(aux_bcent(i,j,k,0)) < small_value) aux_bcent(i,j,k,0) = 0.0;
939  if (std::abs(aux_bcent(i,j,k,1)) < small_value) aux_bcent(i,j,k,1) = 0.0;
940  if (std::abs(aux_bcent(i,j,k,2)) < small_value) aux_bcent(i,j,k,2) = 0.0;
941  });
942 
943  // Area fraction MultiFab has one more slice at bigEnd(idim),
944  // and this slice is not filled by fillBoundary(), for higher levels.
945  // (Lower level might be filled by fillBoundary().)
946  // Fill the ghost region for the last slice at bigEnd(idim)
947  // by the value of the nearst point. And let fillBoundary() overwrite it.
948 
949  Box upper_slab = makeSlab(bx_grown, a_idim, bx.bigEnd(a_idim)+1);
950  Box bx_grown_1 = bx; bx_grown_1.grow(a_idim,1);
951  BoxList slab_diffList = boxDiff(upper_slab, bx_grown_1);
952 
953  for (const Box& b : slab_diffList) {
954  ParallelFor(b, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
955  {
956  IntVect iv(AMREX_D_DECL(i,j,k));
957  IntVect iv_nearest = iv;
958  for (int d=0; d<AMREX_SPACEDIM; ++d) {
959  iv_nearest[d] = Clamp(iv[d], bx_grown_1.smallEnd(d), bx_grown_1.bigEnd(d));
960  }
961  aux_afrac_x(iv) = aux_afrac_x(iv_nearest);
962  });
963  }
964 
965  } // if (FlagFab[mfi].getType(bx) == FabType::singlevalued )
966 
967  } // MFIter
968 
969  // Fill Boundary
970 
971  // The FB call for volfrac is done above
972  // m_volfrac->FillBoundary(a_geom.periodicity());
973 
974  m_volcent->FillBoundary(a_geom.periodicity());
975  for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
976  m_areafrac[idim]->FillBoundary(a_geom.periodicity());
977  m_facecent[idim]->FillBoundary(a_geom.periodicity());
978  }
979  m_bndryarea->FillBoundary(a_geom.periodicity());
980  m_bndrycent->FillBoundary(a_geom.periodicity());
981  m_bndrynorm->FillBoundary(a_geom.periodicity());
982 
983  // Set Connectivities
984  for (MFIter mfi(*m_cellflags, false); mfi.isValid(); ++mfi) {
985 
986  const Box& bx = mfi.validbox();
987  const Box domain = surroundingNodes(a_geom.Domain(), a_idim);
988 
989  if (FlagFab[mfi].getType(bx) == FabType::singlevalued ) {
990 
991  Array4<EBCellFlag> const& aux_flag = m_cellflags->array(mfi);
992  Array4<Real> const& aux_afrac_x = m_areafrac[0]->array(mfi);
993  Array4<Real> const& aux_afrac_y = m_areafrac[1]->array(mfi);
994  Array4<Real> const& aux_afrac_z = m_areafrac[2]->array(mfi);
995 
996  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
997  {
998  EB2::build_cellflag_from_ap (i, j, k, aux_flag, aux_afrac_x, aux_afrac_y, aux_afrac_z);
999  });
1000 
1001  // Set disconnected non-periodicfaces
1002 
1003  bool l_periodic_x = a_geom.isPeriodic(0);
1004  bool l_periodic_y = a_geom.isPeriodic(1);
1005  bool l_periodic_z = a_geom.isPeriodic(2);
1006 
1007  if (!l_periodic_x) {
1008  const Box dom_grown = grow(grow(domain,1,1),2,1);
1009  const Box bx_grown = grow(grow( bx,1,1),2,1);
1010  const Box bx_face_x_lo = bx_grown & makeSlab(dom_grown,0,domain.smallEnd(0));
1011  const Box bx_face_x_hi = bx_grown & makeSlab(dom_grown,0,domain.bigEnd(0));
1012 
1013  ParallelFor(bx_face_x_lo, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
1014  {
1015  for(int kk(-1); kk<=1; kk++) {
1016  for(int jj(-1); jj<=1; jj++) {
1017  aux_flag(i,j,k).setDisconnected(-1,jj,kk);
1018  }}
1019  });
1020  ParallelFor(bx_face_x_hi, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
1021  {
1022  for(int kk(-1); kk<=1; kk++) {
1023  for(int jj(-1); jj<=1; jj++) {
1024  aux_flag(i,j,k).setDisconnected( 1,jj,kk);
1025  }}
1026  });
1027  }
1028 
1029  if (!l_periodic_y) {
1030  const Box dom_grown = grow(grow(domain,0,1),2,1);
1031  const Box bx_grown = grow(grow( bx,0,1),2,1);
1032  const Box bx_face_y_lo = bx_grown & makeSlab(dom_grown,1,domain.smallEnd(1));
1033  const Box bx_face_y_hi = bx_grown & makeSlab(dom_grown,1,domain.bigEnd(1));
1034 
1035  ParallelFor(bx_face_y_lo, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
1036  {
1037  for(int kk(-1); kk<=1; kk++) {
1038  for(int ii(-1); ii<=1; ii++) {
1039  aux_flag(i,j,k).setDisconnected(ii,-1,kk);
1040  }}
1041  });
1042  ParallelFor(bx_face_y_hi, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
1043  {
1044  for(int kk(-1); kk<=1; kk++) {
1045  for(int ii(-1); ii<=1; ii++) {
1046  aux_flag(i,j,k).setDisconnected(ii, 1,kk);
1047  }}
1048  });
1049  }
1050 
1051  if (!l_periodic_z) {
1052  const Box dom_grown = grow(grow(domain,0,1),1,1);
1053  const Box bx_grown = grow(grow( bx,0,1),1,1);
1054  const Box bx_face_z_lo = bx_grown & makeSlab(dom_grown,2,domain.smallEnd(2));
1055  const Box bx_face_z_hi = bx_grown & makeSlab(dom_grown,2,domain.bigEnd(2));
1056 
1057  ParallelFor(bx_face_z_lo, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
1058  {
1059  for(int jj(-1); jj<=1; jj++) {
1060  for(int ii(-1); ii<=1; ii++) {
1061  aux_flag(i,j,k).setDisconnected(ii,jj,-1);
1062  }}
1063  });
1064  ParallelFor(bx_face_z_hi, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
1065  {
1066  for(int jj(-1); jj<=1; jj++) {
1067  for(int ii(-1); ii<=1; ii++) {
1068  aux_flag(i,j,k).setDisconnected(ii,jj, 1);
1069  }}
1070  });
1071  }
1072 
1073  } // FabType::singlevalued
1074 
1075  } // MFIter
1076 
1077  // Set disconnected zero-volume-fraction cells
1078  // (equivalent to eb_::set_connection_flags for CC grids)
1079 
1080  for (MFIter mfi(*m_cellflags, false); mfi.isValid(); ++mfi) {
1081 
1082  const Box& bx = mfi.validbox();
1083  const Box gbx = amrex::grow(bx, m_cellflags->nGrow()-1); // Leave one cell layer
1084 
1085  Array4<EBCellFlag> const& aux_flag = m_cellflags->array(mfi);
1086  Array4<Real> const& aux_vfrac = m_volfrac->array(mfi);
1087 
1088  ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
1089  {
1090  for(int kk(-1); kk<=1; kk++) {
1091  for(int jj(-1); jj<=1; jj++) {
1092  for(int ii(-1); ii<=1; ii++)
1093  {
1094  if (aux_vfrac(i+ii,j+jj,k+kk) == 0.0) {
1095  aux_flag(i,j,k).setDisconnected(ii,jj,kk);
1096  }
1097  }}}
1098  });
1099 
1100  ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
1101  {
1102  if (aux_vfrac(i,j,k)==0.0) {
1103  aux_flag(i,j,k).setCovered();
1104  }
1105  });
1106 
1107  } // MFIter
1108 
1109  // Fill Boundary
1110 
1111  m_cellflags->FillBoundary(a_geom.periodicity());
1112 
1113 }
pp
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real pp(amrex::Real y)
Definition: ERF_MicrophysicsUtils.H:233
Real
amrex::Real Real
Definition: ERF_ShocInterface.H:19
eb_aux_::m_cellflags
amrex::FabArray< amrex::EBCellFlagFab > * m_cellflags
Definition: ERF_EBAux.H:44
eb_aux_::m_bndrynorm
amrex::MultiFab * m_bndrynorm
Definition: ERF_EBAux.H:49
eb_aux_::m_bndryarea
amrex::MultiFab * m_bndryarea
Definition: ERF_EBAux.H:47
eb_aux_::m_volfrac
amrex::MultiFab * m_volfrac
Definition: ERF_EBAux.H:45
eb_aux_::m_bndrycent
amrex::MultiFab * m_bndrycent
Definition: ERF_EBAux.H:48
eb_aux_::m_facecent
amrex::Array< amrex::MultiFab *, AMREX_SPACEDIM > m_facecent
Definition: ERF_EBAux.H:52
eb_aux_::m_volcent
amrex::MultiFab * m_volcent
Definition: ERF_EBAux.H:46
eb_aux_::m_areafrac
amrex::Array< amrex::MultiFab *, AMREX_SPACEDIM > m_areafrac
Definition: ERF_EBAux.H:51
eb_cut_cell_
Definition: ERF_EBCutCell.H:14
module_model_constants::epsilon
real(c_double), parameter epsilon
Definition: ERF_module_model_constants.F90:12

Referenced by eb_::make_all_factories().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ getAreaFrac()

Array< const MultiFab *, AMREX_SPACEDIM > eb_aux_::getAreaFrac ( ) const
1159 {
1160  AMREX_ASSERT(m_areafrac[0] != nullptr);
1161  return {AMREX_D_DECL(m_areafrac[0], m_areafrac[1], m_areafrac[2])};
1162 }

Referenced by redistribute_term().

Here is the caller graph for this function:

◆ getBndryArea()

const MultiFab & eb_aux_::getBndryArea ( ) const
1138 {
1139  AMREX_ASSERT(m_bndryarea != nullptr);
1140  return *m_bndryarea;
1141 }

◆ getBndryCent()

const MultiFab & eb_aux_::getBndryCent ( ) const
1145 {
1146  AMREX_ASSERT(m_bndrycent != nullptr);
1147  return *m_bndrycent;
1148 }

◆ getBndryNorm()

const MultiFab & eb_aux_::getBndryNorm ( ) const
1152 {
1153  AMREX_ASSERT(m_bndrynorm != nullptr);
1154  return *m_bndrynorm;
1155 }

◆ getCentroid()

const MultiFab & eb_aux_::getCentroid ( ) const
1131 {
1132  AMREX_ASSERT(m_volcent != nullptr);
1133  return *m_volcent;
1134 }

Referenced by redistribute_term().

Here is the caller graph for this function:

◆ getFaceCent()

Array< const MultiFab *, AMREX_SPACEDIM > eb_aux_::getFaceCent ( ) const
1166 {
1167  AMREX_ASSERT(m_facecent[0] != nullptr);
1168  return {AMREX_D_DECL(m_facecent[0], m_facecent[1], m_facecent[2])};
1169 }

Referenced by redistribute_term().

Here is the caller graph for this function:

◆ getMultiEBCellFlagFab()

const FabArray< EBCellFlagFab > & eb_aux_::getMultiEBCellFlagFab ( ) const
1117 {
1118  AMREX_ASSERT(m_cellflags != nullptr);
1119  return *m_cellflags;
1120 }

Referenced by redistribute_term().

Here is the caller graph for this function:

◆ getVolFrac()

const MultiFab & eb_aux_::getVolFrac ( ) const
1124 {
1125  AMREX_ASSERT(m_volfrac != nullptr);
1126  return *m_volfrac;
1127 }

Referenced by redistribute_term().

Here is the caller graph for this function:

◆ set_verbose()

void eb_aux_::set_verbose ( )
inline
26 { m_verbose = 1; }

Member Data Documentation

◆ m_areafrac

amrex::Array<amrex::MultiFab*,AMREX_SPACEDIM> eb_aux_::m_areafrac {{AMREX_D_DECL(nullptr, nullptr, nullptr)}}
private

Referenced by define(), and getAreaFrac().

◆ m_bndryarea

amrex::MultiFab* eb_aux_::m_bndryarea = nullptr
private

Referenced by define(), and getBndryArea().

◆ m_bndrycent

amrex::MultiFab* eb_aux_::m_bndrycent = nullptr
private

Referenced by define(), and getBndryCent().

◆ m_bndrynorm

amrex::MultiFab* eb_aux_::m_bndrynorm = nullptr
private

Referenced by define(), and getBndryNorm().

◆ m_cellflags

amrex::FabArray<amrex::EBCellFlagFab>* eb_aux_::m_cellflags = nullptr
private

Referenced by define(), and getMultiEBCellFlagFab().

◆ m_facecent

amrex::Array<amrex::MultiFab*,AMREX_SPACEDIM> eb_aux_::m_facecent {{AMREX_D_DECL(nullptr, nullptr, nullptr)}}
private

Referenced by define(), and getFaceCent().

◆ m_verbose

int eb_aux_::m_verbose
private

Referenced by define(), and set_verbose().

◆ m_volcent

amrex::MultiFab* eb_aux_::m_volcent = nullptr
private

Referenced by define(), and getCentroid().

◆ m_volfrac

amrex::MultiFab* eb_aux_::m_volfrac = nullptr
private

Referenced by define(), and getVolFrac().

The documentation for this class was generated from the following files: