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_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:39

◆ ~eb_aux_()

eb_aux_::~eb_aux_ ( )
10 {
11 }

Member Function Documentation

◆ define()

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

Referenced by redistribute_term().

Here is the caller graph for this function:

◆ getBndryArea()

const MultiFab & eb_aux_::getBndryArea ( ) const
1066 {
1067  AMREX_ASSERT(m_bndryarea != nullptr);
1068  return *m_bndryarea;
1069 }

◆ getBndryCent()

const MultiFab & eb_aux_::getBndryCent ( ) const
1073 {
1074  AMREX_ASSERT(m_bndrycent != nullptr);
1075  return *m_bndrycent;
1076 }

◆ getBndryNorm()

const MultiFab & eb_aux_::getBndryNorm ( ) const
1080 {
1081  AMREX_ASSERT(m_bndrynorm != nullptr);
1082  return *m_bndrynorm;
1083 }

◆ getCentroid()

const MultiFab & eb_aux_::getCentroid ( ) const
1059 {
1060  AMREX_ASSERT(m_volcent != nullptr);
1061  return *m_volcent;
1062 }

Referenced by redistribute_term().

Here is the caller graph for this function:

◆ getFaceCent()

Array< const MultiFab *, AMREX_SPACEDIM > eb_aux_::getFaceCent ( ) const
1094 {
1095  AMREX_ASSERT(m_facecent[0] != nullptr);
1096  return {AMREX_D_DECL(m_facecent[0], m_facecent[1], m_facecent[2])};
1097 }

Referenced by redistribute_term().

Here is the caller graph for this function:

◆ getMultiEBCellFlagFab()

const FabArray< EBCellFlagFab > & eb_aux_::getMultiEBCellFlagFab ( ) const
1045 {
1046  AMREX_ASSERT(m_cellflags != nullptr);
1047  return *m_cellflags;
1048 }

Referenced by redistribute_term().

Here is the caller graph for this function:

◆ getVolFrac()

const MultiFab & eb_aux_::getVolFrac ( ) const
1052 {
1053  AMREX_ASSERT(m_volfrac != nullptr);
1054  return *m_volfrac;
1055 }

Referenced by redistribute_term().

Here is the caller graph for this function:

◆ set_verbose()

void eb_aux_::set_verbose ( )
inline
25 { 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: