Skip to content

Test getEBFluxes in CellEB2 #5347

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

Open
eebasso wants to merge 12 commits into
base: development
Choose a base branch
from
+138 −4
Open

Test getEBFluxes in CellEB2 #5347

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
12 commits
Select commit Hold shift + click to select a range
d79d49c
Add regression test for getEBFluxes
eebasso Apr 17, 2026
ea06409
Use barea in writePlotfile
eebasso Apr 17, 2026
d794dd8
Print norms of ebflux error when verbose. Use barea for scaling
eebasso Apr 18, 2026
a6ba496
Clean up mytest_set_ebflux
eebasso Apr 18, 2026
a97d1b6
Adjustments in MyTest::solve
eebasso Apr 18, 2026
0eaaad8
Update comments and gradphi calculation in mytest_set_fluxeb
eebasso Apr 18, 2026
d23c772
Delete accidental code, update comments on grad phi derivation
eebasso Apr 18, 2026
2a024a6
Use apx, apy, apz instead of bndrynorm to compute boundary area norma…
eebasso Apr 18, 2026
3c145d9
Merge branch 'development' of https://github.com/AMReX-Codes/amrex in…
eebasso Apr 22, 2026
c8b062f
Set exact fluxeb to 0 for regular fabtype boxes
eebasso Apr 22, 2026
b6fa994
Rename function to mytest_set_fluxeb_exact
eebasso Apr 22, 2026
fc79fd3
Use const for ebdata values
eebasso Apr 22, 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
2 changes: 2 additions & 0 deletions Tests/LinearSolvers/CellEB2/MyTest.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -65,6 +65,8 @@ private:
amrex::Vector<amrex::MultiFab> acoef;
amrex::Vector<amrex::Array<amrex::MultiFab,AMREX_SPACEDIM> > bcoef;
amrex::Vector<amrex::MultiFab> bcoef_eb;
amrex::Vector<amrex::MultiFab> fluxeb_sol;
amrex::Vector<amrex::MultiFab> fluxeb_exact;
};

#endif
65 changes: 61 additions & 4 deletions Tests/LinearSolvers/CellEB2/MyTest.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -113,6 +113,9 @@ MyTest::solve ()
const Real tol_rel = reltol;
const Real tol_abs = 0.0;
mlmg.solve(amrex::GetVecOfPtrs(phi), amrex::GetVecOfConstPtrs(rhs), tol_rel, tol_abs);

// Test getEBFluxes
mleb.getEBFluxes(amrex::GetVecOfPtrs(fluxeb_sol), amrex::GetVecOfPtrs(phi));
}
else
{
Expand Down Expand Up @@ -155,6 +158,9 @@ MyTest::solve ()
const Real tol_rel = reltol;
const Real tol_abs = 0.0;
mlmg.solve({&phi[ilev]}, {&rhs[ilev]}, tol_rel, tol_abs);

// Test getEBFluxes
mleb.getEBFluxes({&fluxeb_sol[ilev]}, {&phi[ilev]});
}
}

Expand All @@ -173,6 +179,23 @@ MyTest::solve ()
Real norm1 = mf.norm1()*AMREX_D_TERM((1.0/n_cell), *(1.0/n_cell), *(1.0/n_cell));
amrex::Print() << "Level " << ilev << ": weighted max and 1 norms " << norminf << ", " << norm1 << '\n';
}
for (int ilev = 0; ilev <= max_level; ++ilev)
{
const MultiFab& barea = factory[ilev]->getBndryArea().ToMultiFab(0.0, 0.0);
const Real nxyzi = AMREX_D_TERM((1.0/n_cell), *(1.0/n_cell), *(1.0/n_cell));

MultiFab mf_feb_err(fluxeb_sol[ilev].boxArray(), fluxeb_sol[ilev].DistributionMap(), 1, 0);
MultiFab::Copy(mf_feb_err, fluxeb_sol[ilev], 0, 0, 1, 0);
MultiFab::Subtract(mf_feb_err, fluxeb_exact[ilev], 0, 0, 1, 0);

amrex::Print() << "Level " << ilev << ": EB flux error norms "
<< mf_feb_err.norm0() << ", " << mf_feb_err.norm1()*nxyzi << '\n';

MultiFab::Multiply(mf_feb_err, barea, 0, 0, 1, 0);

amrex::Print() << "Level " << ilev << ": EB flux error weighted norms "
<< mf_feb_err.norm0() << ", " << mf_feb_err.norm1()*nxyzi << '\n';
}
}
}

Expand All @@ -188,20 +211,25 @@ MyTest::writePlotfile ()
if (gpu_regtest) {
for (int ilev = 0; ilev <= max_level; ++ilev) {
const MultiFab& vfrc = factory[ilev]->getVolFrac();
plotmf[ilev].define(grids[ilev],dmap[ilev],3,0);
const MultiFab& barea = factory[ilev]->getBndryArea().ToMultiFab(0.0, 0.0);
plotmf[ilev].define(grids[ilev],dmap[ilev], 6, 0);
MultiFab::Copy(plotmf[ilev], phi[ilev], 0, 0, 1, 0);
MultiFab::Copy(plotmf[ilev], phiexact[ilev], 0, 1, 1, 0);
MultiFab::Copy(plotmf[ilev], vfrc, 0, 2, 1, 0);
MultiFab::Copy(plotmf[ilev], fluxeb_sol[ilev], 0, 3, 1, 0);
MultiFab::Copy(plotmf[ilev], fluxeb_exact[ilev], 0, 4, 1, 0);
MultiFab::Copy(plotmf[ilev], barea, 0, 5, 1, 0);
}
WriteMultiLevelPlotfile(plot_file_name, max_level+1,
amrex::GetVecOfConstPtrs(plotmf),
{"phi","exact","vfrac"},
{"phi","exact","vfrac","eb flux","eb flux exact","barea"},
geom, 0.0, Vector<int>(max_level+1,0),
Vector<IntVect>(max_level,IntVect{2}));
} else {
for (int ilev = 0; ilev <= max_level; ++ilev) {
const MultiFab& vfrc = factory[ilev]->getVolFrac();
plotmf[ilev].define(grids[ilev],dmap[ilev],5,0);
const MultiFab& barea = factory[ilev]->getBndryArea().ToMultiFab(0.0, 0.0);
plotmf[ilev].define(grids[ilev],dmap[ilev], 10, 0);

MultiFab::Copy(plotmf[ilev], phi[ilev], 0, 0, 1, 0);

Expand All @@ -215,10 +243,27 @@ MyTest::writePlotfile ()
MultiFab::Multiply(plotmf[ilev], vfrc, 0, 3, 1, 0);

MultiFab::Copy(plotmf[ilev], vfrc, 0, 4, 1, 0);

// EB fluxes
MultiFab::Copy(plotmf[ilev], fluxeb_sol[ilev], 0, 5, 1, 0);

MultiFab::Copy(plotmf[ilev], fluxeb_exact[ilev], 0, 6, 1, 0);

MultiFab::Copy( plotmf[ilev], fluxeb_sol[ilev], 0, 7, 1, 0);
MultiFab::Subtract(plotmf[ilev], fluxeb_exact[ilev], 0, 7, 1, 0);

MultiFab::Copy( plotmf[ilev], fluxeb_sol[ilev], 0, 8, 1, 0);
MultiFab::Subtract(plotmf[ilev], fluxeb_exact[ilev], 0, 8, 1, 0);
MultiFab::Multiply(plotmf[ilev], barea, 0, 8, 1, 0);

MultiFab::Copy(plotmf[ilev], barea, 0, 9, 1, 0);
}
WriteMultiLevelPlotfile(plot_file_name, max_level+1,
amrex::GetVecOfConstPtrs(plotmf),
{"phi","exact","error","error*vfrac","vfrac"},
{
"phi", "exact", "error", "error*vfrac", "vfrac",
"eb flux", "eb flux exact", "flux error", "(flux error)*barea", "barea"
},
geom, 0.0, Vector<int>(max_level+1,0),
Vector<IntVect>(max_level,IntVect{2}));
}
Expand Down Expand Up @@ -318,6 +363,8 @@ MyTest::initData ()
acoef.resize(nlevels);
bcoef.resize(nlevels);
bcoef_eb.resize(nlevels);
fluxeb_sol.resize(nlevels);
fluxeb_exact.resize(nlevels);

for (int ilev = 0; ilev < nlevels; ++ilev)
{
Expand All @@ -338,13 +385,16 @@ MyTest::initData ()
}
bcoef_eb[ilev].define(grids[ilev], dmap[ilev], 1, 0, MFInfo(), *factory[ilev]);
bcoef_eb[ilev].setVal(1.0);
fluxeb_sol[ilev].define(grids[ilev], dmap[ilev], 1, 1, MFInfo(), *factory[ilev]);
fluxeb_exact[ilev].define(grids[ilev], dmap[ilev], 1, 0, MFInfo(), *factory[ilev]);

phi[ilev].setVal(0.0);
rhs[ilev].setVal(0.0);
acoef[ilev].setVal(0.0);
for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
bcoef[ilev][idim].setVal(1.0);
}
fluxeb_sol[ilev].setVal(0.0);

const auto dx = geom[ilev].CellSizeArray();
const Box& domainbox = geom[ilev].Domain();
Expand All @@ -359,6 +409,7 @@ MyTest::initData ()
#endif
for (MFIter mfi(phiexact[ilev]); mfi.isValid(); ++mfi)
{
const amrex::EBData& ebdata = factory[ilev]->getEBData(mfi);
const Box& bx = mfi.validbox();
const Box& nbx = amrex::surroundingNodes(bx);
Array4<Real> const& phi_arr = phi[ilev].array(mfi);
Expand All @@ -368,20 +419,25 @@ MyTest::initData ()
AMREX_D_TERM(Array4<Real> const& bx_arr = bcoef[ilev][0].array(mfi);,
Array4<Real> const& by_arr = bcoef[ilev][1].array(mfi);,
Array4<Real> const& bz_arr = bcoef[ilev][2].array(mfi););
Array4<Real> const& feb_ex_arr = fluxeb_exact[ilev].array(mfi);

auto fabtyp = flags[mfi].getType(bx);
if (FabType::covered == fabtyp) {
amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
phi_ex_arr(i,j,k) = 0.0;
phi_eb_arr(i,j,k) = 0.0;
feb_ex_arr(i,j,k) = 0.0;
Copy link
Copy Markdown
Member

@WeiqunZhang WeiqunZhang Apr 18, 2026 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Need to do this for FabType::regular too.

Copy link
Copy Markdown
Member

@WeiqunZhang WeiqunZhang Apr 20, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What I was talking about is this branch, not about regular cells in else branch.

            } else if (FabType::regular == fabtyp) {
                amrex::ParallelFor(nbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                {
                    mytest_set_phi_reg(i,j,k,phi_ex_arr,rhs_arr,
                                       AMREX_D_DECL(bx_arr,by_arr,bz_arr),
                                       dx, lprob_type, bx);
                });

It never touches feb_ex_arr.

Copy link
Copy Markdown
Member

@WeiqunZhang WeiqunZhang Apr 20, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

BTW, you can remove this comment // feb_ex_arr(i,j,k) = 0.0;.

Copy link
Copy Markdown
Contributor Author

@eebasso eebasso Apr 22, 2026 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Okay, should be done now. Would using fluxeb_exact.setVal(0.0) at the same place where fluxeb_sol.setVal(0.0) occurs be appropriate? I noticed that this was not done for phiexact

});
} else if (FabType::regular == fabtyp) {
amrex::ParallelFor(nbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
mytest_set_phi_reg(i,j,k,phi_ex_arr,rhs_arr,
AMREX_D_DECL(bx_arr,by_arr,bz_arr),
dx, lprob_type, bx);
if (bx.contains(IntVect(AMREX_D_DECL(i,j,k)))) {
feb_ex_arr(i,j,k) = 0.0;
}
});
} else {
Array4<Real> const& beb_arr = bcoef_eb[ilev].array(mfi);
Expand All @@ -394,6 +450,7 @@ MyTest::initData ()
AMREX_D_DECL(bx_arr,by_arr,bz_arr),
beb_arr,flag_arr,cent_arr,bcent_arr,
dx, lprob_type, bx);
mytest_set_fluxeb_exact(i,j,k,feb_ex_arr,flag_arr,ebdata,dx,lprob_type,bx);
});
}

Expand Down
75 changes: 75 additions & 0 deletions Tests/LinearSolvers/CellEB2/MyTest_K.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -136,6 +136,81 @@ void mytest_set_phi_eb (int i, int j, int k, amrex::Array4<amrex::Real> const& p
}
}

AMREX_GPU_DEVICE AMREX_FORCE_INLINE
void mytest_set_fluxeb_exact (
int i, int j, int k,
amrex::Array4<amrex::Real> const& fluxeb_exact,
amrex::Array4<amrex::EBCellFlag const> const& flag_arr,
amrex::EBData const& ebdata,
amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> const& dx,
int prob_type, amrex::Box const& vbx)
{
constexpr amrex::Real pi = 3.1415926535897932;
amrex::IntVect iv(AMREX_D_DECL(i,j,k));

if (vbx.contains(iv)) {
const amrex::EBCellFlag flag = flag_arr(i,j,k);

if (flag.isCovered()) {
fluxeb_exact(i,j,k) = 0.0;
} else if (flag.isSingleValued()) {
const amrex::Real bcentx = ebdata.get<amrex::EBData_t::bndrycent>(i,j,k,0);
const amrex::Real bcenty = ebdata.get<amrex::EBData_t::bndrycent>(i,j,k,1);

// Compute boundary area normal vector
const amrex::Real apxm = ebdata.get<amrex::EBData_t::apx>(i ,j,k);
const amrex::Real apxp = ebdata.get<amrex::EBData_t::apx>(i+1,j,k);
amrex::Real dapx = (apxm-apxp)/dx[0];
const amrex::Real apym = ebdata.get<amrex::EBData_t::apy>(i,j ,k);
const amrex::Real apyp = ebdata.get<amrex::EBData_t::apy>(i,j+1,k);
amrex::Real dapy = (apym-apyp)/dx[1];
#if (AMREX_SPACEDIM == 3)
const amrex::Real apzm = ebdata.get<amrex::EBData_t::apz>(i,j,k );
const amrex::Real apzp = ebdata.get<amrex::EBData_t::apz>(i,j,k+1);
amrex::Real dapz = (apzm-apzp)/dx[2];
#endif
amrex::Real anorm = std::sqrt(AMREX_D_TERM(dapx*dapx, + dapy*dapy, + dapz*dapz));
amrex::Real anorminv = amrex::Real(1.0)/anorm;
amrex::Real nx = dapx * anorminv;
amrex::Real ny = dapy * anorminv;

amrex::Real x = (i+0.5+bcentx)*dx[0] - 0.5;
amrex::Real y = (j+0.5+bcenty)*dx[1] - 0.5;
amrex::Real theta = std::atan2(x,y) + 0.5*pi;
amrex::Real r2 = x*x+y*y;

// Derivation of grad phi Cartesian components:
// th = atan2(y,x) is the polar angle
// theta = atan2(x,y) + pi/2
// = pi/2 - atan2(y,x) + pi/2
// = pi - th
// phi = r^4 cos(3 theta)
// = r^4 cos(3 (pi - th))
// grad phi = rhat (+4 r^3 cos(3 theta)) + thhat (+3 r^3 sin(3 theta))
// rhat = (+x/r) xhat + (+y/r) yhat
// thhat = (-y/r) xhat + (+x/r) yhat
// grad_x phi = (+x/r)(+4 r^3 cos(3 theta)) + (-y/r)(+3 r^3 sin(3 theta))
// grad_y phi = (-y/r)(+4 r^3 cos(3 theta)) + (+x/r)(+3 r^3 sin(3 theta))
amrex::Real gradphi_r_over_r = 4.0 * r2 * std::cos(3.*theta);
amrex::Real gradphi_t_over_r = 3.0 * r2 * std::sin(3.*theta);
amrex::Real gradphi_x = x * gradphi_r_over_r - y * gradphi_t_over_r;
amrex::Real gradphi_y = y * gradphi_r_over_r + x * gradphi_t_over_r;

amrex::Real bx = 1.0;
amrex::Real by = 1.0;
if (prob_type == 2) {
bx = 1.0 - r2;
by = 1.0 - r2;
}
// Flux on EB surface = - (area normal).(b tensor).(grad phi)
amrex::Real n_b_gradphi = nx * bx * gradphi_x + ny * by * gradphi_y;
fluxeb_exact(i,j,k) = -n_b_gradphi;
} else {
fluxeb_exact(i,j,k) = 0.0;
}
}
}

AMREX_GPU_DEVICE AMREX_FORCE_INLINE
void mytest_set_phi_boundary (int i, int j, int k, amrex::Array4<amrex::Real> const& phi,
amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> const& dx,
Expand Down
Loading