31 Real small_volfrac = 1.e-14;
33 pp.queryAdd(
"small_volfrac", small_volfrac);
34 const Real small_value = 1.e-15;
36 const IntVect vdim(IntVect::TheDimensionVector(a_idim));
38 const BoxArray& grids = amrex::convert(a_grids, vdim);
40 m_cellflags =
new FabArray<EBCellFlagFab>(grids, a_dmap, 1, a_ngrow[0], MFInfo(),
41 DefaultFabFactory<EBCellFlagFab>());
45 for (MFIter mfi(*
m_cellflags,
false); mfi.isValid(); ++mfi) {
46 auto& fab = (*m_cellflags)[mfi];
47 fab.setType(FabType::singlevalued);
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());
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());
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());
66 for (
int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
74 const auto& FlagFab = a_factory->getMultiEBCellFlagFab();
76 for (MFIter mfi(*
m_cellflags,
false); mfi.isValid(); ++mfi) {
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);
83 GpuArray<Real, AMREX_SPACEDIM> dx = a_geom.CellSizeArray();
84 bool l_periodic = a_geom.isPeriodic(a_idim);
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);
92 if (FlagFab[mfi].getType(bx) == FabType::covered ) {
94 ParallelFor(tbx, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept
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;
104 }
else if (FlagFab[mfi].getType(bx) == FabType::regular ) {
106 ParallelFor(tbx, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept
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;
116 }
else if (FlagFab[mfi].getType(bx) == FabType::singlevalued ) {
121 Array4<EBCellFlag const>
const& flag = FlagFab.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();
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);
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]);
145 const IntVect dom_grown_lo = dom_grown.smallEnd();
146 const IntVect dom_grown_hi = dom_grown.bigEnd();
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
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();
162 #ifndef AMREX_USE_GPU
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
174 aux_flag(i,j,k).setCovered();
175 aux_flag(i,j,k).setDisconnected();
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;
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;
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;
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;
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;
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;
206 aux_barea(i,j,k) = 0.0;
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;
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;
217 IntVect iv_hi(i,j,k);
218 IntVect iv_lo(iv_hi - vdim);
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();
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));
232 if (at_lo_boundary) {
235 lo_isRegular =
false;
236 lo_isSingleValued =
false;
237 }
else if (hi_isRegular) {
238 lo_isCovered =
false;
240 lo_isSingleValued =
false;
241 }
else if (hi_isSingleValued) {
242 if (almostEqual(afrac(i,j,k),0.0)) {
244 lo_isRegular =
false;
245 lo_isSingleValued =
false;
246 }
else if (almostEqual(afrac(i,j,k),1.0)) {
247 lo_isCovered =
false;
249 lo_isSingleValued =
false;
251 lo_isCovered =
false;
252 lo_isRegular =
false;
253 lo_isSingleValued =
true;
261 if (at_hi_boundary) {
264 hi_isRegular =
false;
265 hi_isSingleValued =
false;
266 }
else if (lo_isRegular) {
267 hi_isCovered =
false;
269 hi_isSingleValued =
false;
270 }
else if (lo_isSingleValued) {
271 if (almostEqual(afrac(i,j,k),0.0)) {
273 hi_isRegular =
false;
274 hi_isSingleValued =
false;
275 }
else if (almostEqual(afrac(i,j,k),1.0)) {
276 hi_isCovered =
false;
278 hi_isSingleValued =
false;
280 hi_isCovered =
false;
281 hi_isRegular =
false;
282 hi_isSingleValued =
true;
288 if ( lo_isCovered && hi_isCovered) {
292 }
else if ( lo_isRegular && hi_isRegular) {
294 aux_flag(i,j,k).setRegular();
295 aux_flag(i,j,k).setConnected();
297 aux_vfrac(i,j,k) = 1.0;
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;
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;
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;
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;
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;
325 #ifndef AMREX_USE_GPU
326 if (verbose) { Print() <<
"\ncell: " << amrex::IntVect(i,j,k) <<
"\n"; }
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};
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));
341 if (at_lo_boundary) {
342 lo_point[idim] += 1.0;
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];
350 Real norm = sqrt( bnorm_x*bnorm_x + bnorm_y*bnorm_y + bnorm_z*bnorm_z);
352 RealVect bnorm_isoparam ( bnorm_x / norm, bnorm_y / norm, bnorm_z / norm);
354 lo_normal = bnorm_isoparam;
360 RealBox lo_rbx(lo_arr.data(), hi_arr.data());
362 eb_cut_cell_ lo_eb_cc(flag(iv_lo), lo_rbx, lo_point, lo_normal);
366 AMREX_ASSERT( !lo_isCovered || lo_eb_cc.isCovered() );
367 AMREX_ASSERT( !lo_isRegular || lo_eb_cc.isRegular() );
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));
376 if (at_hi_boundary) {
377 hi_point[idim] += -1.0;
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];
385 Real norm = sqrt( bnorm_x*bnorm_x + bnorm_y*bnorm_y + bnorm_z*bnorm_z);
387 RealVect bnorm_isoparam ( bnorm_x / norm, bnorm_y / norm, bnorm_z / norm);
389 hi_normal = bnorm_isoparam;
395 RealBox hi_rbx(lo_arr.data(), hi_arr.data());
397 eb_cut_cell_ hi_eb_cc(flag(iv_hi), hi_rbx, hi_point, hi_normal);
401 AMREX_ASSERT( !hi_isCovered || hi_eb_cc.isCovered() );
402 AMREX_ASSERT( !hi_isRegular || hi_eb_cc.isRegular() );
405 #if defined(AMREX_DEBUG) || defined(AMREX_TESTING) || 1
415 eb_cut_cell_ hi_hi_eb_cc(flag(iv_hi), lo_rbx, hi_point, hi_normal);
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()
430 AMREX_ALWAYS_ASSERT( !hi_isRegular || hi_hi_eb_cc.isRegular() );
431 AMREX_ALWAYS_ASSERT( !hi_isCovered || hi_hi_eb_cc.isCovered() );
438 if ( hi_isSingleValued ) {
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];
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];
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];
455 Real const apnorm = std::sqrt(adx*adx + ady*ady + adz*adz);
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);
462 #ifndef AMREX_USE_GPU
463 if ( !amrex::almostEqual(dot_normals, 1.0) ) {
464 Print() <<
"\nFail: check-1 dot_normals " << dot_normals
470 }
else if (verbose) {
471 Print() <<
"Pass: dot_normals = 1.0\n";
475 AMREX_ALWAYS_ASSERT( amrex::almostEqual(dot_normals, 1.0) );
480 #ifndef AMREX_USE_GPU
481 Real const abs_err = std::abs( hi_eb_cc.areaHi(idim) - hi_hi_eb_cc.areaLo(idim) );
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)
488 }
else if (verbose) {
489 Print() <<
"Pass: hi_eb_cc.areaHi = hi_hi_eb_cc.areaLo"
490 <<
" abs_err: " << abs_err <<
"\n";
492 AMREX_ALWAYS_ASSERT( abs_err < machine_tol );
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 ) {
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)
508 }
else if (verbose) {
509 Print() <<
"Pass: hi_eb_cc.areaLo = afrac = " << afrac(iv_hi)
510 <<
" abs_err: " << abs_err <<
"\n";
513 AMREX_ALWAYS_ASSERT( abs_err < compare_tol );
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 ) {
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)
532 }
else if (verbose) {
533 Print() <<
"Pass: hi_eb_cc + hi_hi_eb_cc = vfrac = " << vfrac(iv_hi)
534 <<
" abs_err: " << abs_err <<
"\n";
537 AMREX_ALWAYS_ASSERT( abs_err < compare_tol );
547 if (lo_eb_cc.isCovered() && hi_eb_cc.isCovered()) {
551 }
else if (lo_eb_cc.isRegular() && hi_eb_cc.isRegular()) {
553 aux_flag(i,j,k).setRegular();
554 aux_flag(i,j,k).setConnected();
556 aux_vfrac(i,j,k) = 1.0;
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;
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;
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;
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;
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;
579 }
else if ( (lo_eb_cc.isRegular() && hi_eb_cc.isCovered())
580 || (lo_eb_cc.isCovered() && hi_eb_cc.isRegular()) ) {
583 #ifndef AMREX_USE_GPU
584 Print()<<
"eb_aux_ / Check: Regular and Covered cut cells are facing each other." << std::endl;
591 aux_flag(i,j,k).setSingleValued();
595 Real lo_vol {lo_eb_cc.volume()};
596 Real hi_vol {hi_eb_cc.volume()};
598 aux_vfrac(i,j,k) = lo_vol + hi_vol;
610 RealVect lo_vcent {lo_eb_cc.centVol()};
611 RealVect hi_vcent {hi_eb_cc.centVol()};
613 lo_vcent[idim] = lo_vcent[idim] - 0.5;
614 hi_vcent[idim] = hi_vcent[idim] + 0.5;
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);
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)};
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)};
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;
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;
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;
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;
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)};
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)};
669 aux_fcent_x(i,j,k,0) = lo_centLo_x[1];
670 aux_fcent_x(i,j,k,1) = lo_centLo_x[2];
671 aux_fcent_y(i,j,k,0) = (aux_afrac_y(i,j,k) > 0.0)
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)
675 aux_fcent_y(i,j,k,1) = (aux_afrac_y(i,j,k) > 0.0)
676 ? ( lo_areaLo_y * lo_centLo_y[2]
677 + hi_areaLo_y * hi_centLo_y[2] ) / aux_afrac_y(i,j,k)
679 aux_fcent_z(i,j,k,0) = (aux_afrac_z(i,j,k) > 0.0)
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)
683 aux_fcent_z(i,j,k,1) = (aux_afrac_z(i,j,k) > 0.0)
684 ? ( lo_areaLo_z * lo_centLo_z[1]
685 + hi_areaLo_z * hi_centLo_z[1] ) / aux_afrac_z(i,j,k)
687 }
else if (idim == 1) {
688 aux_fcent_x(i,j,k,0) = (aux_afrac_x(i,j,k) > 0.0)
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)
692 aux_fcent_x(i,j,k,1) = (aux_afrac_x(i,j,k) > 0.0)
693 ? ( lo_areaLo_x * lo_centLo_x[2]
694 + hi_areaLo_x * hi_centLo_x[2] ) / aux_afrac_x(i,j,k)
696 aux_fcent_y(i,j,k,0) = lo_centLo_y[0];
697 aux_fcent_y(i,j,k,1) = lo_centLo_y[2];
698 aux_fcent_z(i,j,k,0) = (aux_afrac_z(i,j,k) > 0.0)
699 ? ( lo_areaLo_z * lo_centLo_z[0]
700 + hi_areaLo_z * hi_centLo_z[0] ) / aux_afrac_z(i,j,k)
702 aux_fcent_z(i,j,k,1) = (aux_afrac_z(i,j,k) > 0.0)
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)
706 }
else if (idim == 2) {
707 aux_fcent_x(i,j,k,0) = (aux_afrac_x(i,j,k) > 0.0)
708 ? ( lo_areaLo_x * lo_centLo_x[1]
709 + hi_areaLo_x * hi_centLo_x[1] ) / aux_afrac_x(i,j,k)
711 aux_fcent_x(i,j,k,1) = (aux_afrac_x(i,j,k) > 0.0)
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)
715 aux_fcent_y(i,j,k,0) = (aux_afrac_y(i,j,k) > 0.0)
716 ? ( lo_areaLo_y * lo_centLo_y[0]
717 + hi_areaLo_y * hi_centLo_y[0] ) / aux_afrac_y(i,j,k)
719 aux_fcent_y(i,j,k,1) = (aux_afrac_y(i,j,k) > 0.0)
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)
723 aux_fcent_z(i,j,k,0) = lo_centLo_z[0];
724 aux_fcent_z(i,j,k,1) = lo_centLo_z[1];
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)};
733 aux_fcent_x(i+1,j,k,0) = hi_centHi_x[1];
734 aux_fcent_x(i+1,j,k,1) = hi_centHi_x[2];
735 }
else if (idim == 1) {
736 aux_fcent_x(i+1,j,k,0) = (aux_afrac_x(i+1,j,k) > 0.0)
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)
740 aux_fcent_x(i+1,j,k,1) = (aux_afrac_x(i+1,j,k) > 0.0)
741 ? ( lo_areaHi_x * lo_centHi_x[2]
742 + hi_areaHi_x * hi_centHi_x[2] ) / aux_afrac_x(i+1,j,k)
744 }
else if (idim == 2) {
745 aux_fcent_x(i+1,j,k,0) = (aux_afrac_x(i+1,j,k) > 0.0)
746 ? ( lo_areaHi_x * lo_centHi_x[1]
747 + hi_areaHi_x * hi_centHi_x[1] ) / aux_afrac_x(i+1,j,k)
749 aux_fcent_x(i+1,j,k,1) = (aux_afrac_x(i+1,j,k) > 0.0)
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)
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)};
761 aux_fcent_y(i,j+1,k,0) = (aux_afrac_y(i,j+1,k) > 0.0)
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)
765 aux_fcent_y(i,j+1,k,1) = (aux_afrac_y(i,j+1,k) > 0.0)
766 ? ( lo_areaHi_y * lo_centHi_y[2]
767 + hi_areaHi_y * hi_centHi_y[2] ) / aux_afrac_y(i,j+1,k)
769 }
else if (idim == 1) {
770 aux_fcent_y(i,j+1,k,0) = lo_centHi_y[0];
771 aux_fcent_y(i,j+1,k,1) = lo_centHi_y[2];
772 }
else if (idim == 2) {
773 aux_fcent_y(i,j+1,k,0) = (aux_afrac_y(i,j+1,k) > 0.0)
774 ? ( lo_areaHi_y * lo_centHi_y[0]
775 + hi_areaHi_y * hi_centHi_y[0] ) / aux_afrac_y(i,j+1,k)
777 aux_fcent_y(i,j+1,k,1) = (aux_afrac_y(i,j+1,k) > 0.0)
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)
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)};
789 aux_fcent_z(i,j,k+1,0) = (aux_afrac_z(i,j,k+1) > 0.0)
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)
793 aux_fcent_z(i,j,k+1,1) = (aux_afrac_z(i,j,k+1) > 0.0)
794 ? ( lo_areaHi_z * lo_centHi_z[1]
795 + hi_areaHi_z * hi_centHi_z[1] ) / aux_afrac_z(i,j,k+1)
797 }
else if (idim == 1) {
798 aux_fcent_z(i,j,k+1,0) = (aux_afrac_z(i,j,k+1) > 0.0)
799 ? ( lo_areaHi_z * lo_centHi_z[0]
800 + hi_areaHi_z * hi_centHi_z[0] ) / aux_afrac_z(i,j,k+1)
802 aux_fcent_z(i,j,k+1,1) = (aux_afrac_z(i,j,k+1) > 0.0)
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)
806 }
else if (idim == 2) {
807 aux_fcent_z(i,j,k+1,0) = lo_centHi_z[0];
808 aux_fcent_z(i,j,k+1,1) = lo_centHi_z[1];
814 Real lo_areaBoun {lo_eb_cc.areaBoun()};
815 Real hi_areaBoun {hi_eb_cc.areaBoun()};
817 aux_barea(i,j,k) = lo_areaBoun + hi_areaBoun;
821 RealVect lo_centBoun {lo_eb_cc.centBoun()};
822 RealVect hi_centBoun {hi_eb_cc.centBoun()};
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);
826 aux_bcent(i,j,k,1) = ( lo_areaBoun * lo_centBoun[1] + hi_areaBoun * hi_centBoun[1] ) / aux_barea(i,j,k);
827 aux_bcent(i,j,k,2) = ( lo_areaBoun * lo_centBoun[2] + hi_areaBoun * hi_centBoun[2] ) / aux_barea(i,j,k);
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);
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);
831 aux_bcent(i,j,k,2) = ( lo_areaBoun * lo_centBoun[2] + hi_areaBoun * hi_centBoun[2] ) / aux_barea(i,j,k);
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);
834 aux_bcent(i,j,k,1) = ( lo_areaBoun * lo_centBoun[1] + hi_areaBoun * hi_centBoun[1] ) / aux_barea(i,j,k);
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);
840 RealVect eb_normal = ( lo_areaBoun * lo_normal + hi_areaBoun * hi_normal )/ aux_barea(i,j,k);
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];
854 ParallelFor(bx, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept
856 if (aux_vfrac(i,j,k) < small_volfrac) {
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;
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;
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;
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;
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;
885 aux_barea(i,j,k) = 0.0;
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;
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;
895 aux_flag(i,j,k).setCovered();
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;
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());
924 for (MFIter mfi(*
m_cellflags,
false); mfi.isValid(); ++mfi) {
926 const Box& bx = mfi.validbox();
927 const Box domain = surroundingNodes(a_geom.Domain(), a_idim);
929 if (FlagFab[mfi].getType(bx) == FabType::singlevalued ) {
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);
936 ParallelFor(bx, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept
938 EB2::build_cellflag_from_ap (i, j, k, aux_flag, aux_afrac_x, aux_afrac_y, aux_afrac_z);
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);
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));
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;
960 ParallelFor(bx_face_x_lo, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept
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);
967 ParallelFor(bx_face_x_hi, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept
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);
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));
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;
989 ParallelFor(bx_face_y_lo, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept
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);
996 ParallelFor(bx_face_y_hi, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept
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);
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));
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;
1018 ParallelFor(bx_face_z_lo, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept
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);
1025 ParallelFor(bx_face_z_hi, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept
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);
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real pp(amrex::Real y)
Definition: ERF_MicrophysicsUtils.H:230
amrex::Real Real
Definition: ERF_ShocInterface.H:16
amrex::FabArray< amrex::EBCellFlagFab > * m_cellflags
Definition: ERF_EBAux.H:43
amrex::MultiFab * m_bndrynorm
Definition: ERF_EBAux.H:48
amrex::MultiFab * m_bndryarea
Definition: ERF_EBAux.H:46
amrex::MultiFab * m_volfrac
Definition: ERF_EBAux.H:44
amrex::MultiFab * m_bndrycent
Definition: ERF_EBAux.H:47
amrex::Array< amrex::MultiFab *, AMREX_SPACEDIM > m_facecent
Definition: ERF_EBAux.H:51
amrex::MultiFab * m_volcent
Definition: ERF_EBAux.H:45
amrex::Array< amrex::MultiFab *, AMREX_SPACEDIM > m_areafrac
Definition: ERF_EBAux.H:50
Definition: ERF_EBCutCell.H:14
real(c_double), parameter epsilon
Definition: ERF_module_model_constants.F90:12