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

◆ ~eb_aux_()

eb_aux_::~eb_aux_ ( )
10 {
11 }

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

Referenced by redistribute_term().

Here is the caller graph for this function:

◆ getBndryArea()

const MultiFab & eb_aux_::getBndryArea ( ) const
1101 {
1102  AMREX_ASSERT(m_bndryarea != nullptr);
1103  return *m_bndryarea;
1104 }

◆ getBndryCent()

const MultiFab & eb_aux_::getBndryCent ( ) const
1108 {
1109  AMREX_ASSERT(m_bndrycent != nullptr);
1110  return *m_bndrycent;
1111 }

◆ getBndryNorm()

const MultiFab & eb_aux_::getBndryNorm ( ) const
1115 {
1116  AMREX_ASSERT(m_bndrynorm != nullptr);
1117  return *m_bndrynorm;
1118 }

◆ getCentroid()

const MultiFab & eb_aux_::getCentroid ( ) const
1094 {
1095  AMREX_ASSERT(m_volcent != nullptr);
1096  return *m_volcent;
1097 }

Referenced by redistribute_term().

Here is the caller graph for this function:

◆ getFaceCent()

Array< const MultiFab *, AMREX_SPACEDIM > eb_aux_::getFaceCent ( ) const
1129 {
1130  AMREX_ASSERT(m_facecent[0] != nullptr);
1131  return {AMREX_D_DECL(m_facecent[0], m_facecent[1], m_facecent[2])};
1132 }

Referenced by redistribute_term().

Here is the caller graph for this function:

◆ getMultiEBCellFlagFab()

const FabArray< EBCellFlagFab > & eb_aux_::getMultiEBCellFlagFab ( ) const
1080 {
1081  AMREX_ASSERT(m_cellflags != nullptr);
1082  return *m_cellflags;
1083 }

Referenced by redistribute_term().

Here is the caller graph for this function:

◆ getVolFrac()

const MultiFab & eb_aux_::getVolFrac ( ) const
1087 {
1088  AMREX_ASSERT(m_volfrac != nullptr);
1089  return *m_volfrac;
1090 }

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: