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

Referenced by redistribute_term().

Here is the caller graph for this function:

◆ getBndryArea()

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

◆ getBndryCent()

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

◆ getBndryNorm()

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

◆ getCentroid()

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

Referenced by redistribute_term().

Here is the caller graph for this function:

◆ getFaceCent()

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

Referenced by redistribute_term().

Here is the caller graph for this function:

◆ getMultiEBCellFlagFab()

const FabArray< EBCellFlagFab > & eb_aux_::getMultiEBCellFlagFab ( ) const
1110 {
1111  AMREX_ASSERT(m_cellflags != nullptr);
1112  return *m_cellflags;
1113 }

Referenced by redistribute_term().

Here is the caller graph for this function:

◆ getVolFrac()

const MultiFab & eb_aux_::getVolFrac ( ) const
1117 {
1118  AMREX_ASSERT(m_volfrac != nullptr);
1119  return *m_volfrac;
1120 }

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: