Skip to content

EB state redistribution: Fix Itracker abort with vfrac fallback #5630

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
WeiqunZhang merged 9 commits into from
Aug 27, 2026
Merged

EB state redistribution: Fix Itracker abort with vfrac fallback #5630

Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
78b386b
Set qt=0.0 outside domain.
skang67 Aug 15, 2026
8fd3d6a
bug fix qt.
skang67 Aug 15, 2026
d727e03
bugfix qt.
skang67 Aug 15, 2026
eceb5ac
first draft of neighbor selection based on vfrac.
skang67 Aug 16, 2026
0388af8
Clean up.
skang67 Aug 16, 2026
6131291
Remove whitespaces.
skang67 Aug 20, 2026
62926f9
Move variables for 2D.
skang67 Aug 20, 2026
5de0379
Remove whitespaces.
skang67 Aug 20, 2026
9d79917
Adjust indentation.
skang67 Aug 20, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
328 changes: 238 additions & 90 deletions Src/EB/AMReX_EB_StateRedistItracker.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -117,9 +117,53 @@ MakeITracker ( Box const& bx,
int ioff = imap[itracker(i,j,k,1)];
int joff = jmap[itracker(i,j,k,1)];

// Sanity check
if (vfrac(i+ioff,j+joff,k) == 0.) {
amrex::Abort(" Trying to merge with covered cell");
// Variables for vfrac fallback
bool select_by_vfrac = false;
Real second_vfrac = Real(-1.0);
int second_idx = -1;

// Sanity check - if selected neighbor is covered, use vfrac fallback
if (vfrac(i+ioff,j+joff,k) == 0.)
{
Real vfrac_x = (xdir_pls_ok) ? vfrac(i+1,j,k) : Real(-1.0);
int idx_x = 5;
if (xdir_mns_ok && vfrac(i-1,j,k) > vfrac_x) {
vfrac_x = vfrac(i-1,j,k);
idx_x = 4;
}

Real vfrac_y = (ydir_pls_ok) ? vfrac(i,j+1,k) : Real(-1.0);
int idx_y = 7;
if (ydir_mns_ok && vfrac(i,j-1,k) > vfrac_y) {
vfrac_y = vfrac(i,j-1,k);
idx_y = 2;
}

// Select largest and second-largest vfrac
int first_idx;
Real first_vfrac;
if (vfrac_x >= vfrac_y) {
first_idx = idx_x;
first_vfrac = vfrac_x;
second_idx = idx_y;
second_vfrac = vfrac_y;
} else {
first_idx = idx_y;
first_vfrac = vfrac_y;
second_idx = idx_x;
second_vfrac = vfrac_x;
}

if (first_vfrac <= 0.0) {
amrex::Abort("Trying to merge with covered cell - all cardinal neighbors are also covered or outside domain");
}

itracker(i,j,k,1) = first_idx;

ioff = imap[first_idx];
joff = jmap[first_idx];

select_by_vfrac = true;
}

Real sum_vol = vfrac(i,j,k) + vfrac(i+ioff,j+joff,k);
Expand All @@ -135,50 +179,63 @@ MakeITracker ( Box const& bx,
// If the merged cell isn't large enough, we try to merge in the other direction
if (sum_vol < target_volfrac || nx_eq_ny)
{
// Original offset was in y-direction, so we will add to the x-direction
// Note that if we can't because it would go outside the domain, we don't
if (ioff == 0) {
if (nx >= 0 && xdir_pls_ok)
{
itracker(i,j,k,2) = 5;
itracker(i,j,k,0) += 1;
}
else if (nx <= 0 && xdir_mns_ok)
{
itracker(i,j,k,2) = 4;
itracker(i,j,k,0) += 1;
}
// If first neighbor was selected by vfrac, use second-best vfrac candidate for second neighbor
if (select_by_vfrac && second_vfrac > 0.0)
{
itracker(i,j,k,2) = second_idx;
itracker(i,j,k,0) += 1;
int ioff2 = imap[second_idx];
int joff2 = jmap[second_idx];
sum_vol += vfrac(i+ioff2,j+joff2,k);
}
else if (!select_by_vfrac)
{
// Original normal-based perpendicular selection logic
// Original offset was in y-direction, so we will add to the x-direction
// Note that if we can't because it would go outside the domain, we don't
if (ioff == 0) {
if (nx >= 0 && xdir_pls_ok)
{
itracker(i,j,k,2) = 5;
itracker(i,j,k,0) += 1;
}
else if (nx <= 0 && xdir_mns_ok)
{
itracker(i,j,k,2) = 4;
itracker(i,j,k,0) += 1;
}

// Original offset was in x-direction, so we will add to the y-direction
// Note that if we can't because it would go outside the domain, we don't
} else {
if (ny >= 0 && ydir_pls_ok)
{
itracker(i,j,k,2) = 7;
itracker(i,j,k,0) += 1;
}
else if (ny <= 0 && ydir_mns_ok)
{
itracker(i,j,k,2) = 2;
itracker(i,j,k,0) += 1;
// Original offset was in x-direction, so we will add to the y-direction
// Note that if we can't because it would go outside the domain, we don't
} else {
if (ny >= 0 && ydir_pls_ok)
{
itracker(i,j,k,2) = 7;
itracker(i,j,k,0) += 1;
}
else if (ny <= 0 && ydir_mns_ok)
{
itracker(i,j,k,2) = 2;
itracker(i,j,k,0) += 1;
}
}
}

if (itracker(i,j,k,0) > 1)
{
// (i+ioff2,j+joff2) is in the nbhd of (i,j)
int ioff2 = imap[itracker(i,j,k,2)];
int joff2 = jmap[itracker(i,j,k,2)];
if (itracker(i,j,k,0) > 1)
{
// (i+ioff2,j+joff2) is in the nbhd of (i,j)
int ioff2 = imap[itracker(i,j,k,2)];
int joff2 = jmap[itracker(i,j,k,2)];

sum_vol += vfrac(i+ioff2,j+joff2,k);
sum_vol += vfrac(i+ioff2,j+joff2,k);
#if 0
if (debug_verbose > 0)
amrex::Print() << "Cell " << IntVect(i,j) << " with volfrac " << vfrac(i,j,k) <<
" trying to ALSO merge with " << IntVect(i+ioff2,j+joff2) <<
" with volfrac " << vfrac(i+ioff2,j+joff2,k) <<
" to get new sum_vol " << sum_vol << '\n';
if (debug_verbose > 0)
amrex::Print() << "Cell " << IntVect(i,j) << " with volfrac " << vfrac(i,j,k) <<
" trying to ALSO merge with " << IntVect(i+ioff2,j+joff2) <<
" with volfrac " << vfrac(i+ioff2,j+joff2,k) <<
" to get new sum_vol " << sum_vol << '\n';
#endif
}
}
} // end else if (!select_by_vfrac)
}

// Now we merge in the corner direction if we have already claimed two
Expand Down Expand Up @@ -381,11 +438,86 @@ MakeITracker ( Box const& bx,
int joff = jmap[itracker(i,j,k,1)];
int koff = kmap[itracker(i,j,k,1)];

// Sanity check
// Variables for vfrac fallback
bool select_by_vfrac = false;
Real second_vfrac = Real(-1.0);
int second_idx = -1;

// Sanity check - if selected neighbor is covered, use vfrac fallback
if (vfrac(i+ioff,j+joff,k+koff) == 0.)
{
// amrex::Print() << "Cell " << IntVect(i,j,k) << " is trying to merge with cell " << IntVect(i+ioff,j+joff,k+koff) << '\n';
amrex::Abort(" Trying to merge with covered cell");
Real vfrac_x = (xdir_pls_ok) ? vfrac(i+1,j,k) : Real(-1.0);
int idx_x = 5;
if (xdir_mns_ok && vfrac(i-1,j,k) > vfrac_x) { // Strict > means +x wins ties
vfrac_x = vfrac(i-1,j,k);
idx_x = 4;
}

Real vfrac_y = (ydir_pls_ok) ? vfrac(i,j+1,k) : Real(-1.0);
int idx_y = 7;
if (ydir_mns_ok && vfrac(i,j-1,k) > vfrac_y) {
vfrac_y = vfrac(i,j-1,k);
idx_y = 2;
}

Real vfrac_z = (zdir_pls_ok) ? vfrac(i,j,k+1) : Real(-1.0);
int idx_z = 22;
if (zdir_mns_ok && vfrac(i,j,k-1) > vfrac_z) {
vfrac_z = vfrac(i,j,k-1);
idx_z = 13;
}

// Find largest and second-largest vfrac
int first_idx;
Real first_vfrac;

if (vfrac_x >= vfrac_y && vfrac_x >= vfrac_z) {
// X is largest
first_idx = idx_x;
first_vfrac = vfrac_x;
if (vfrac_y >= vfrac_z) {
second_idx = idx_y;
second_vfrac = vfrac_y;
} else {
second_idx = idx_z;
second_vfrac = vfrac_z;
}
} else if (vfrac_y >= vfrac_z) {
// Y is largest
first_idx = idx_y;
first_vfrac = vfrac_y;
if (vfrac_x >= vfrac_z) {
second_idx = idx_x;
second_vfrac = vfrac_x;
} else {
second_idx = idx_z;
second_vfrac = vfrac_z;
}
} else {
// Z is largest
first_idx = idx_z;
first_vfrac = vfrac_z;
if (vfrac_x >= vfrac_y) {
second_idx = idx_x;
second_vfrac = vfrac_x;
} else {
second_idx = idx_y;
second_vfrac = vfrac_y;
}
}

if (first_vfrac <= 0.0) {
amrex::Abort("Trying to merge with covered cell - all cardinal neighbors are also covered or outside domain");
}

// Update itracker and offsets
itracker(i,j,k,1) = first_idx;

ioff = imap[first_idx];
joff = jmap[first_idx];
koff = kmap[first_idx];

select_by_vfrac = true;
}

Real sum_vol = vfrac(i,j,k) + vfrac(i+ioff,j+joff,k+koff);
Expand All @@ -408,62 +540,78 @@ MakeITracker ( Box const& bx,

if ( (sum_vol < target_volfrac) || just_broke_symmetry )
{
// Original offset was in x-direction
if (joff == 0 && koff == 0)
// If first neighbor was selected by vfrac, use second-best vfrac candidate for second neighbor
if (select_by_vfrac && second_vfrac > 0.0)
{
if (nx_eq_ny) {
itracker(i,j,k,2) = (ny > 0) ? 7 : 2;
} else if (nx_eq_nz) {
itracker(i,j,k,2) = (nz > 0) ? 22 : 13;
} else if ( (std::abs(ny) > std::abs(nz)) ) {
itracker(i,j,k,2) = (ny > 0) ? 7 : 2;
} else {
itracker(i,j,k,2) = (nz > 0) ? 22 : 13;
}
itracker(i,j,k,2) = second_idx;
itracker(i,j,k,0) += 1;

// Original offset was in y-direction
} else if (ioff == 0 && koff == 0)
{
if (nx_eq_ny) {
itracker(i,j,k,2) = (nx > 0) ? 5 : 4;
} else if (ny_eq_nz) {
itracker(i,j,k,2) = (nz > 0) ? 22 : 13;
} else if ( (std::abs(nx) > std::abs(nz)) ) {
itracker(i,j,k,2) = (nx > 0) ? 5 : 4;
} else {
itracker(i,j,k,2) = (nz > 0) ? 22 : 13;
}
int ioff2 = imap[second_idx];
int joff2 = jmap[second_idx];
int koff2 = kmap[second_idx];
sum_vol += vfrac(i+ioff2,j+joff2,k+koff2);

// Original offset was in z-direction
} else if (ioff == 0 && joff == 0)
}
else if (!select_by_vfrac)
{
if (nx_eq_nz) {
itracker(i,j,k,2) = (nx > 0) ? 5 : 4;
} else if (ny_eq_nz) {
itracker(i,j,k,2) = (ny > 0) ? 7 : 2;
} else if ( (std::abs(nx) > std::abs(ny)) ) {
itracker(i,j,k,2) = (nx > 0) ? 5 : 4;
} else {
itracker(i,j,k,2) = (ny > 0) ? 7 : 2;
// Original normal-based perpendicular selection logic
// Original offset was in x-direction
if (joff == 0 && koff == 0)
{
if (nx_eq_ny) {
itracker(i,j,k,2) = (ny > 0) ? 7 : 2;
} else if (nx_eq_nz) {
itracker(i,j,k,2) = (nz > 0) ? 22 : 13;
} else if ( (std::abs(ny) > std::abs(nz)) ) {
itracker(i,j,k,2) = (ny > 0) ? 7 : 2;
} else {
itracker(i,j,k,2) = (nz > 0) ? 22 : 13;
}

// Original offset was in y-direction
} else if (ioff == 0 && koff == 0)
{
if (nx_eq_ny) {
itracker(i,j,k,2) = (nx > 0) ? 5 : 4;
} else if (ny_eq_nz) {
itracker(i,j,k,2) = (nz > 0) ? 22 : 13;
} else if ( (std::abs(nx) > std::abs(nz)) ) {
itracker(i,j,k,2) = (nx > 0) ? 5 : 4;
} else {
itracker(i,j,k,2) = (nz > 0) ? 22 : 13;
}

// Original offset was in z-direction
} else if (ioff == 0 && joff == 0)
{
if (nx_eq_nz) {
itracker(i,j,k,2) = (nx > 0) ? 5 : 4;
} else if (ny_eq_nz) {
itracker(i,j,k,2) = (ny > 0) ? 7 : 2;
} else if ( (std::abs(nx) > std::abs(ny)) ) {
itracker(i,j,k,2) = (nx > 0) ? 5 : 4;
} else {
itracker(i,j,k,2) = (ny > 0) ? 7 : 2;
}
}
}

// (i,j,k) merges with at least two cells now
itracker(i,j,k,0) += 1;
// (i,j,k) merges with at least two cells now
itracker(i,j,k,0) += 1;

// (i+ioff2,j+joff2,k+koff2) is in the nbhd of (i,j,k)
int ioff2 = imap[itracker(i,j,k,2)];
int joff2 = jmap[itracker(i,j,k,2)];
int koff2 = kmap[itracker(i,j,k,2)];
// (i+ioff2,j+joff2,k+koff2) is in the nbhd of (i,j,k)
int ioff2 = imap[itracker(i,j,k,2)];
int joff2 = jmap[itracker(i,j,k,2)];
int koff2 = kmap[itracker(i,j,k,2)];

sum_vol += vfrac(i+ioff2,j+joff2,k+koff2);
sum_vol += vfrac(i+ioff2,j+joff2,k+koff2);
#if 0
if (debug_print)
amrex::Print() << "Cell " << IntVect(i,j,k) << " with volfrac " << vfrac(i,j,k) <<
" trying to ALSO merge with " << IntVect(i+ioff2,j+joff2,k+koff2) <<
" with volfrac " << vfrac(i+ioff2,j+joff2,k+koff2) <<
" to get new sum_vol " << sum_vol << '\n';
if (debug_print)
amrex::Print() << "Cell " << IntVect(i,j,k) << " with volfrac " << vfrac(i,j,k) <<
" trying to ALSO merge with " << IntVect(i+ioff2,j+joff2,k+koff2) <<
" with volfrac " << vfrac(i+ioff2,j+joff2,k+koff2) <<
" to get new sum_vol " << sum_vol << '\n';
#endif
} // end else if (!select_by_vfrac)
}

// If the merged cell has merged in two directions, we now merge in the corner direction within the current plane
Expand Down
Loading