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 5 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
63 changes: 59 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,13 +419,15 @@ 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
Expand All @@ -394,6 +447,8 @@ MyTest::initData ()
AMREX_D_DECL(bx_arr,by_arr,bz_arr),
beb_arr,flag_arr,cent_arr,bcent_arr,
dx, lprob_type, bx);
// feb_ex_arr(i,j,k) = 0.0;
mytest_set_fluxeb(i,j,k,feb_ex_arr,flag_arr,ebdata,dx,lprob_type,bx);
});
}

Expand Down
55 changes: 55 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,61 @@ 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 (
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)) {
amrex::EBCellFlag const flag = flag_arr(i,j,k);

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

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;

// Compute grad phi Cartesian components
// 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 = 4 r^3 cos(3 theta) rhat + 3 r^3 sin(3 theta) thhat
// rhat = (x/r) xhat + (y/r) yhat
// thhat = (-y/r) xhat + (x/r) yhat
amrex::Real gradphi_x = r2 * (x * 4.0 * std::cos(3.*theta) - y * 3.0 * std::sin(3.*theta));
amrex::Real gradphi_y = r2 * (y * 4.0 * std::cos(3.*theta) + x * 3.0 * std::sin(3.*theta));

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