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_ ( )
19  : m_verbose(0)
20 // ,m_defined(0)
21 {}
eb_aux_::m_verbose
int m_verbose
Definition: ERF_EBAux.H:40

◆ ~eb_aux_()

eb_aux_::~eb_aux_ ( )
14 {
15 }

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

Referenced by redistribute_term().

Here is the caller graph for this function:

◆ getBndryArea()

const MultiFab & eb_aux_::getBndryArea ( ) const
1173 {
1174  AMREX_ASSERT(m_bndryarea != nullptr);
1175  return *m_bndryarea;
1176 }

◆ getBndryCent()

const MultiFab & eb_aux_::getBndryCent ( ) const
1180 {
1181  AMREX_ASSERT(m_bndrycent != nullptr);
1182  return *m_bndrycent;
1183 }

◆ getBndryNorm()

const MultiFab & eb_aux_::getBndryNorm ( ) const
1187 {
1188  AMREX_ASSERT(m_bndrynorm != nullptr);
1189  return *m_bndrynorm;
1190 }

◆ getCentroid()

const MultiFab & eb_aux_::getCentroid ( ) const
1166 {
1167  AMREX_ASSERT(m_volcent != nullptr);
1168  return *m_volcent;
1169 }

Referenced by redistribute_term().

Here is the caller graph for this function:

◆ getFaceCent()

Array< const MultiFab *, AMREX_SPACEDIM > eb_aux_::getFaceCent ( ) const
1201 {
1202  AMREX_ASSERT(m_facecent[0] != nullptr);
1203  return {AMREX_D_DECL(m_facecent[0], m_facecent[1], m_facecent[2])};
1204 }

Referenced by redistribute_term().

Here is the caller graph for this function:

◆ getMultiEBCellFlagFab()

const FabArray< EBCellFlagFab > & eb_aux_::getMultiEBCellFlagFab ( ) const
1152 {
1153  AMREX_ASSERT(m_cellflags != nullptr);
1154  return *m_cellflags;
1155 }

Referenced by redistribute_term().

Here is the caller graph for this function:

◆ getVolFrac()

const MultiFab & eb_aux_::getVolFrac ( ) const
1159 {
1160  AMREX_ASSERT(m_volfrac != nullptr);
1161  return *m_volfrac;
1162 }

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: