From d79d49cfecbe4a0ee90bccd5d4aaed986f9bcbc0 Mon Sep 17 00:00:00 2001 From: Edward Basso Date: Fri, 17 Apr 2026 16:41:13 -0700 Subject: [PATCH 01/11] Add regression test for getEBFluxes --- Tests/LinearSolvers/CellEB2/MyTest.H | 2 + Tests/LinearSolvers/CellEB2/MyTest.cpp | 34 +++++++- Tests/LinearSolvers/CellEB2/MyTest_K.H | 108 +++++++++++++++++++++++++ 3 files changed, 140 insertions(+), 4 deletions(-) diff --git a/Tests/LinearSolvers/CellEB2/MyTest.H b/Tests/LinearSolvers/CellEB2/MyTest.H index c91c3d5cc11..ddd1378d8a7 100644 --- a/Tests/LinearSolvers/CellEB2/MyTest.H +++ b/Tests/LinearSolvers/CellEB2/MyTest.H @@ -65,6 +65,8 @@ private: amrex::Vector acoef; amrex::Vector > bcoef; amrex::Vector bcoef_eb; + amrex::Vector fluxeb_sol; + amrex::Vector fluxeb_exact; }; #endif diff --git a/Tests/LinearSolvers/CellEB2/MyTest.cpp b/Tests/LinearSolvers/CellEB2/MyTest.cpp index 8b9bb19d38b..b295ff2757a 100644 --- a/Tests/LinearSolvers/CellEB2/MyTest.cpp +++ b/Tests/LinearSolvers/CellEB2/MyTest.cpp @@ -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); + + // Get EB fluxes + mleb.getEBFluxes(amrex::GetVecOfPtrs(fluxeb_sol), amrex::GetVecOfPtrs(phi)); } else { @@ -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); + + // Get EB fluxes + mleb.getEBFluxes({&fluxeb_sol[ilev]}, {&phi[ilev]}); } } @@ -188,20 +194,22 @@ 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); + plotmf[ilev].define(grids[ilev],dmap[ilev], 5, 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); } WriteMultiLevelPlotfile(plot_file_name, max_level+1, amrex::GetVecOfConstPtrs(plotmf), - {"phi","exact","vfrac"}, + {"phi","exact","vfrac","EB flux sol","EB flux exact"}, geom, 0.0, Vector(max_level+1,0), Vector(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); + plotmf[ilev].define(grids[ilev],dmap[ilev], 8, 0); MultiFab::Copy(plotmf[ilev], phi[ilev], 0, 0, 1, 0); @@ -215,10 +223,18 @@ 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); } WriteMultiLevelPlotfile(plot_file_name, max_level+1, amrex::GetVecOfConstPtrs(plotmf), - {"phi","exact","error","error*vfrac","vfrac"}, + {"phi","exact","error","error*vfrac","vfrac","EB flux sol","EB flux exact","EB flux error"}, geom, 0.0, Vector(max_level+1,0), Vector(max_level,IntVect{2})); } @@ -318,6 +334,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) { @@ -338,6 +356,8 @@ 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); @@ -345,6 +365,7 @@ MyTest::initData () 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(); @@ -359,6 +380,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 const& phi_arr = phi[ilev].array(mfi); @@ -368,6 +390,7 @@ MyTest::initData () AMREX_D_TERM(Array4 const& bx_arr = bcoef[ilev][0].array(mfi);, Array4 const& by_arr = bcoef[ilev][1].array(mfi);, Array4 const& bz_arr = bcoef[ilev][2].array(mfi);); + Array4 const& feb_ex_arr = fluxeb_exact[ilev].array(mfi); auto fabtyp = flags[mfi].getType(bx); if (FabType::covered == fabtyp) { @@ -375,6 +398,7 @@ MyTest::initData () { phi_ex_arr(i,j,k) = 0.0; phi_eb_arr(i,j,k) = 0.0; + feb_ex_arr(i,j,k) = 0.0; }); } else if (FabType::regular == fabtyp) { amrex::ParallelFor(nbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept @@ -394,6 +418,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); }); } diff --git a/Tests/LinearSolvers/CellEB2/MyTest_K.H b/Tests/LinearSolvers/CellEB2/MyTest_K.H index deafabf9d7d..377ee771a47 100644 --- a/Tests/LinearSolvers/CellEB2/MyTest_K.H +++ b/Tests/LinearSolvers/CellEB2/MyTest_K.H @@ -136,6 +136,114 @@ void mytest_set_phi_eb (int i, int j, int k, amrex::Array4 const& p } } +AMREX_GPU_DEVICE AMREX_FORCE_INLINE +void mytest_set_fluxeb ( + int i, int j, int k, + amrex::Array4 const& fluxeb_exact, + amrex::Array4 const& flag_arr, + amrex::EBData const& ebdata, + amrex::GpuArray 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(i,j,k,0); + // amrex::Real const ny = ebdata.get(i,j,k,1); + amrex::Real const bcentx = ebdata.get(i,j,k,0); + amrex::Real const bcenty = ebdata.get(i,j,k,1); + +#if (AMREX_SPACEDIM != 3) + amrex::Real apxm = ebdata.get(i ,j,k); + amrex::Real apxp = ebdata.get(i+1,j,k); + amrex::Real apym = ebdata.get(i,j ,k); + amrex::Real apyp = ebdata.get(i,j+1,k); + amrex::Real dapx = (apxm-apxp)*dx[1]; + amrex::Real dapy = (apym-apyp)*dx[0]; + amrex::Real anorm = std::hypot(dapx,dapy) + 1.e-30*std::sqrt(dx[0]*dx[1]); +#else + amrex::Real apxm = ebdata.get(i ,j,k); + amrex::Real apxp = ebdata.get(i+1,j,k); + amrex::Real apym = ebdata.get(i,j ,k); + amrex::Real apyp = ebdata.get(i,j+1,k); + amrex::Real apzm = ebdata.get(i,j,k ); + amrex::Real apzp = ebdata.get(i,j,k+1); + amrex::Real dapx = (apxm-apxp)*dx[1]*dx[2]; + amrex::Real dapy = (apym-apyp)*dx[0]*dx[2]; + amrex::Real dapz = (apzm-apzp)*dx[0]*dx[1]; + amrex::Real anorm = std::sqrt(dapx*dapx+dapy*dapy+dapz*dapz); +#endif + amrex::Real anorminv = amrex::Real(1.0)/anorm; + amrex::Real nx = dapx * anorminv; + amrex::Real ny = dapy * anorminv; + + amrex::Real const nx_ebdata = ebdata.get(i,j,k,0); + amrex::Real const ny_ebdata = ebdata.get(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; + amrex::Real r = std::sqrt(r2); + amrex::Real r3 = r*r2; + + // amrex::Real gradphi_x = 0.5 * r3 * (-7.0*std::cos(2.*theta) - std::cos(4.*theta)); + // amrex::Real gradphi_y = 0.5 * r3 * (-7.0*std::sin(2.*theta) + std::sin(4.*theta)); + + 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 th = atan2(y,x); + amrex::Real costh = std::cos(th); + amrex::Real sinth = std::sin(th); + + // theta = pi - th + // phi = r^4 cos(theta) = r^4 cos(3*(pi - th)) + amrex::Real gradphi_r = +4.0 * r3 * std::cos(3.*theta); + amrex::Real gradphi_t = +3.0 * r3 * std::sin(3.*theta); + amrex::Real gradphi_cyl_x = costh*gradphi_r - sinth*gradphi_t; + amrex::Real gradphi_cyl_y = sinth*gradphi_r + costh*gradphi_t; + + if (std::hypot(nx - nx_ebdata, ny - ny_ebdata) >= 1.0e-7) { + amrex::Print() << "\n"; + amrex::Print() << "i j k: " << i << " " << j << " " << k << "\n"; + amrex::Print() << " x = " << x << ", y = " << y << "\n"; + amrex::Print() << " r = " << r << ", th (degrees) = " << 180.0*th/pi << "\n"; + amrex::Print() << " n from ap = " << nx << ", " << ny << "\n"; + amrex::Print() << " n_ebdata = " << nx_ebdata << ", " << ny_ebdata << "\n"; + amrex::Print() << "\n"; + } + + // nx = nx_ebdata; + // ny = ny_ebdata; + // BL_ASSERT(std::abs(nx - nx_ebdata) < 1.0e-7); + // BL_ASSERT(std::abs(ny - ny_ebdata) < 1.0e-7); + + BL_ASSERT(std::abs(gradphi_x - gradphi_cyl_x) < 1.0e-9); + BL_ASSERT(std::abs(gradphi_y - gradphi_cyl_y) < 1.0e-9); + + 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; + } + // fluxeb_exact(i,j,k) = 7.0; + } +} + AMREX_GPU_DEVICE AMREX_FORCE_INLINE void mytest_set_phi_boundary (int i, int j, int k, amrex::Array4 const& phi, amrex::GpuArray const& dx, From ea0640979d037d6d5d063a5802e95587c400f44c Mon Sep 17 00:00:00 2001 From: Edward Basso Date: Fri, 17 Apr 2026 16:52:56 -0700 Subject: [PATCH 02/11] Use barea in writePlotfile --- Tests/LinearSolvers/CellEB2/MyTest.cpp | 20 ++++++++++++++++---- 1 file changed, 16 insertions(+), 4 deletions(-) diff --git a/Tests/LinearSolvers/CellEB2/MyTest.cpp b/Tests/LinearSolvers/CellEB2/MyTest.cpp index b295ff2757a..77ead73887d 100644 --- a/Tests/LinearSolvers/CellEB2/MyTest.cpp +++ b/Tests/LinearSolvers/CellEB2/MyTest.cpp @@ -194,22 +194,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], 5, 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, 2, 1, 0); } WriteMultiLevelPlotfile(plot_file_name, max_level+1, amrex::GetVecOfConstPtrs(plotmf), - {"phi","exact","vfrac","EB flux sol","EB flux exact"}, + {"phi","exact","vfrac","eb flux","eb flux exact","barea"}, geom, 0.0, Vector(max_level+1,0), Vector(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], 8, 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); @@ -231,10 +234,19 @@ MyTest::writePlotfile () 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","EB flux sol","EB flux exact","EB flux error"}, + { + "phi", "exact", "error", "error*vfrac", "vfrac", + "eb flux", "eb flux exact", "flux error", "(flux error)*barea", "barea" + }, geom, 0.0, Vector(max_level+1,0), Vector(max_level,IntVect{2})); } From d794dd816fdb98b6cdc2cb652986be60f48aec27 Mon Sep 17 00:00:00 2001 From: Edward Basso Date: Fri, 17 Apr 2026 17:07:32 -0700 Subject: [PATCH 03/11] Print norms of ebflux error when verbose. Use barea for scaling --- Tests/LinearSolvers/CellEB2/MyTest.cpp | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/Tests/LinearSolvers/CellEB2/MyTest.cpp b/Tests/LinearSolvers/CellEB2/MyTest.cpp index 77ead73887d..85fd6f27eb6 100644 --- a/Tests/LinearSolvers/CellEB2/MyTest.cpp +++ b/Tests/LinearSolvers/CellEB2/MyTest.cpp @@ -179,6 +179,19 @@ 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); + + 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); + MultiFab::Multiply(mf_feb_err, barea, 0, 0, 1, 0); + + Real norminf_feb = mf_feb_err.norm0(); + Real norm1_feb = mf_feb_err.norm1()*AMREX_D_TERM((1.0/n_cell), *(1.0/n_cell), *(1.0/n_cell)); + amrex::Print() << "Level " << ilev << ": EB flux error norms " << norminf_feb << ", " << norm1_feb << '\n'; + } } } From a6ba4961c452c12521143ff6fea445f688e0f159 Mon Sep 17 00:00:00 2001 From: Edward Basso Date: Fri, 17 Apr 2026 17:34:07 -0700 Subject: [PATCH 04/11] Clean up mytest_set_ebflux --- Tests/LinearSolvers/CellEB2/MyTest_K.H | 77 ++++---------------------- 1 file changed, 12 insertions(+), 65 deletions(-) diff --git a/Tests/LinearSolvers/CellEB2/MyTest_K.H b/Tests/LinearSolvers/CellEB2/MyTest_K.H index 377ee771a47..52d7b82c7b2 100644 --- a/Tests/LinearSolvers/CellEB2/MyTest_K.H +++ b/Tests/LinearSolvers/CellEB2/MyTest_K.H @@ -154,93 +154,40 @@ void mytest_set_fluxeb ( if (flag.isCovered()) { fluxeb_exact(i,j,k) = 0.0; } else if (flag.isSingleValued()) { - // amrex::Real const nx = ebdata.get(i,j,k,0); - // amrex::Real const ny = ebdata.get(i,j,k,1); + amrex::Real const nx = ebdata.get(i,j,k,0); + amrex::Real const ny = ebdata.get(i,j,k,1); amrex::Real const bcentx = ebdata.get(i,j,k,0); amrex::Real const bcenty = ebdata.get(i,j,k,1); -#if (AMREX_SPACEDIM != 3) - amrex::Real apxm = ebdata.get(i ,j,k); - amrex::Real apxp = ebdata.get(i+1,j,k); - amrex::Real apym = ebdata.get(i,j ,k); - amrex::Real apyp = ebdata.get(i,j+1,k); - amrex::Real dapx = (apxm-apxp)*dx[1]; - amrex::Real dapy = (apym-apyp)*dx[0]; - amrex::Real anorm = std::hypot(dapx,dapy) + 1.e-30*std::sqrt(dx[0]*dx[1]); -#else - amrex::Real apxm = ebdata.get(i ,j,k); - amrex::Real apxp = ebdata.get(i+1,j,k); - amrex::Real apym = ebdata.get(i,j ,k); - amrex::Real apyp = ebdata.get(i,j+1,k); - amrex::Real apzm = ebdata.get(i,j,k ); - amrex::Real apzp = ebdata.get(i,j,k+1); - amrex::Real dapx = (apxm-apxp)*dx[1]*dx[2]; - amrex::Real dapy = (apym-apyp)*dx[0]*dx[2]; - amrex::Real dapz = (apzm-apzp)*dx[0]*dx[1]; - amrex::Real anorm = std::sqrt(dapx*dapx+dapy*dapy+dapz*dapz); -#endif - amrex::Real anorminv = amrex::Real(1.0)/anorm; - amrex::Real nx = dapx * anorminv; - amrex::Real ny = dapy * anorminv; - - amrex::Real const nx_ebdata = ebdata.get(i,j,k,0); - amrex::Real const ny_ebdata = ebdata.get(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; - amrex::Real r = std::sqrt(r2); - amrex::Real r3 = r*r2; - - // amrex::Real gradphi_x = 0.5 * r3 * (-7.0*std::cos(2.*theta) - std::cos(4.*theta)); - // amrex::Real gradphi_y = 0.5 * r3 * (-7.0*std::sin(2.*theta) + std::sin(4.*theta)); + // 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 th = atan2(y,x); - amrex::Real costh = std::cos(th); - amrex::Real sinth = std::sin(th); - - // theta = pi - th - // phi = r^4 cos(theta) = r^4 cos(3*(pi - th)) - amrex::Real gradphi_r = +4.0 * r3 * std::cos(3.*theta); - amrex::Real gradphi_t = +3.0 * r3 * std::sin(3.*theta); - amrex::Real gradphi_cyl_x = costh*gradphi_r - sinth*gradphi_t; - amrex::Real gradphi_cyl_y = sinth*gradphi_r + costh*gradphi_t; - - if (std::hypot(nx - nx_ebdata, ny - ny_ebdata) >= 1.0e-7) { - amrex::Print() << "\n"; - amrex::Print() << "i j k: " << i << " " << j << " " << k << "\n"; - amrex::Print() << " x = " << x << ", y = " << y << "\n"; - amrex::Print() << " r = " << r << ", th (degrees) = " << 180.0*th/pi << "\n"; - amrex::Print() << " n from ap = " << nx << ", " << ny << "\n"; - amrex::Print() << " n_ebdata = " << nx_ebdata << ", " << ny_ebdata << "\n"; - amrex::Print() << "\n"; - } - - // nx = nx_ebdata; - // ny = ny_ebdata; - // BL_ASSERT(std::abs(nx - nx_ebdata) < 1.0e-7); - // BL_ASSERT(std::abs(ny - ny_ebdata) < 1.0e-7); - - BL_ASSERT(std::abs(gradphi_x - gradphi_cyl_x) < 1.0e-9); - BL_ASSERT(std::abs(gradphi_y - gradphi_cyl_y) < 1.0e-9); - 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) + // 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; } - // fluxeb_exact(i,j,k) = 7.0; } } From a97d1b6d9a5da8a21e2b15db4fc097b12174b09e Mon Sep 17 00:00:00 2001 From: Edward Basso Date: Fri, 17 Apr 2026 17:55:29 -0700 Subject: [PATCH 05/11] Adjustments in MyTest::solve --- Tests/LinearSolvers/CellEB2/MyTest.cpp | 16 ++++++++++------ 1 file changed, 10 insertions(+), 6 deletions(-) diff --git a/Tests/LinearSolvers/CellEB2/MyTest.cpp b/Tests/LinearSolvers/CellEB2/MyTest.cpp index 85fd6f27eb6..223b0db3abf 100644 --- a/Tests/LinearSolvers/CellEB2/MyTest.cpp +++ b/Tests/LinearSolvers/CellEB2/MyTest.cpp @@ -114,7 +114,7 @@ MyTest::solve () const Real tol_abs = 0.0; mlmg.solve(amrex::GetVecOfPtrs(phi), amrex::GetVecOfConstPtrs(rhs), tol_rel, tol_abs); - // Get EB fluxes + // Test getEBFluxes mleb.getEBFluxes(amrex::GetVecOfPtrs(fluxeb_sol), amrex::GetVecOfPtrs(phi)); } else @@ -159,7 +159,7 @@ MyTest::solve () const Real tol_abs = 0.0; mlmg.solve({&phi[ilev]}, {&rhs[ilev]}, tol_rel, tol_abs); - // Get EB fluxes + // Test getEBFluxes mleb.getEBFluxes({&fluxeb_sol[ilev]}, {&phi[ilev]}); } } @@ -182,15 +182,19 @@ MyTest::solve () 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); - Real norminf_feb = mf_feb_err.norm0(); - Real norm1_feb = mf_feb_err.norm1()*AMREX_D_TERM((1.0/n_cell), *(1.0/n_cell), *(1.0/n_cell)); - amrex::Print() << "Level " << ilev << ": EB flux error norms " << norminf_feb << ", " << norm1_feb << '\n'; + amrex::Print() << "Level " << ilev << ": EB flux error weighted norms " + << mf_feb_err.norm0() << ", " << mf_feb_err.norm1()*nxyzi << '\n'; } } } @@ -214,7 +218,7 @@ MyTest::writePlotfile () 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, 2, 1, 0); + MultiFab::Copy(plotmf[ilev], barea, 0, 5, 1, 0); } WriteMultiLevelPlotfile(plot_file_name, max_level+1, amrex::GetVecOfConstPtrs(plotmf), From 0eaaad8959c322f65c56fd7b4927b056b4f05642 Mon Sep 17 00:00:00 2001 From: Edward Basso Date: Sat, 18 Apr 2026 14:24:46 -0700 Subject: [PATCH 06/11] Update comments and gradphi calculation in mytest_set_fluxeb --- Tests/LinearSolvers/CellEB2/MyTest_K.H | 21 +++++++++++++++------ 1 file changed, 15 insertions(+), 6 deletions(-) diff --git a/Tests/LinearSolvers/CellEB2/MyTest_K.H b/Tests/LinearSolvers/CellEB2/MyTest_K.H index 52d7b82c7b2..f7571ec09a9 100644 --- a/Tests/LinearSolvers/CellEB2/MyTest_K.H +++ b/Tests/LinearSolvers/CellEB2/MyTest_K.H @@ -164,17 +164,22 @@ void mytest_set_fluxeb ( amrex::Real theta = std::atan2(x,y) + 0.5*pi; amrex::Real r2 = x*x+y*y; - // Compute grad phi Cartesian components + // 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 = 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)); + // 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; @@ -191,6 +196,10 @@ void mytest_set_fluxeb ( } } +AMREX_GPU_DEVICE AMREX_FORCE_INLINE +void compute_bnorm(int i, int j, int k, ) + + AMREX_GPU_DEVICE AMREX_FORCE_INLINE void mytest_set_phi_boundary (int i, int j, int k, amrex::Array4 const& phi, amrex::GpuArray const& dx, From d23c7726b4ab3f0b119f962c825668a074766241 Mon Sep 17 00:00:00 2001 From: Edward Basso Date: Sat, 18 Apr 2026 15:29:37 -0700 Subject: [PATCH 07/11] Delete accidental code, update comments on grad phi derivation --- Tests/LinearSolvers/CellEB2/MyTest_K.H | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/Tests/LinearSolvers/CellEB2/MyTest_K.H b/Tests/LinearSolvers/CellEB2/MyTest_K.H index f7571ec09a9..74cf1d32b70 100644 --- a/Tests/LinearSolvers/CellEB2/MyTest_K.H +++ b/Tests/LinearSolvers/CellEB2/MyTest_K.H @@ -170,8 +170,8 @@ void mytest_set_fluxeb ( // = 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)) + // = 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)) @@ -196,10 +196,6 @@ void mytest_set_fluxeb ( } } -AMREX_GPU_DEVICE AMREX_FORCE_INLINE -void compute_bnorm(int i, int j, int k, ) - - AMREX_GPU_DEVICE AMREX_FORCE_INLINE void mytest_set_phi_boundary (int i, int j, int k, amrex::Array4 const& phi, amrex::GpuArray const& dx, From 2a024a68e6a325233ebe29bee99f95eaf56e6e96 Mon Sep 17 00:00:00 2001 From: Edward Basso Date: Sat, 18 Apr 2026 15:44:56 -0700 Subject: [PATCH 08/11] Use apx, apy, apz instead of bndrynorm to compute boundary area normal vector --- Tests/LinearSolvers/CellEB2/MyTest_K.H | 19 +++++++++++++++++-- 1 file changed, 17 insertions(+), 2 deletions(-) diff --git a/Tests/LinearSolvers/CellEB2/MyTest_K.H b/Tests/LinearSolvers/CellEB2/MyTest_K.H index 74cf1d32b70..791e6303a75 100644 --- a/Tests/LinearSolvers/CellEB2/MyTest_K.H +++ b/Tests/LinearSolvers/CellEB2/MyTest_K.H @@ -154,11 +154,26 @@ void mytest_set_fluxeb ( if (flag.isCovered()) { fluxeb_exact(i,j,k) = 0.0; } else if (flag.isSingleValued()) { - amrex::Real const nx = ebdata.get(i,j,k,0); - amrex::Real const ny = ebdata.get(i,j,k,1); amrex::Real const bcentx = ebdata.get(i,j,k,0); amrex::Real const bcenty = ebdata.get(i,j,k,1); + // Compute boundary area normal vector + amrex::Real apxm = ebdata.get(i ,j,k); + amrex::Real apxp = ebdata.get(i+1,j,k); + amrex::Real dapx = (apxm-apxp)/dx[0]; + amrex::Real apym = ebdata.get(i,j ,k); + amrex::Real apyp = ebdata.get(i,j+1,k); + amrex::Real dapy = (apym-apyp)/dx[1]; +#if (AMREX_SPACEDIM == 3) + amrex::Real apzm = ebdata.get(i,j,k ); + amrex::Real apzp = ebdata.get(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; From c8b062febff1a54e7cb91ad64f3ac2812e652c1a Mon Sep 17 00:00:00 2001 From: Edward Basso Date: Wed, 22 Apr 2026 14:03:16 -0700 Subject: [PATCH 09/11] Set exact fluxeb to 0 for regular fabtype boxes --- Tests/LinearSolvers/CellEB2/MyTest.cpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/Tests/LinearSolvers/CellEB2/MyTest.cpp b/Tests/LinearSolvers/CellEB2/MyTest.cpp index 223b0db3abf..f77e0c37c0e 100644 --- a/Tests/LinearSolvers/CellEB2/MyTest.cpp +++ b/Tests/LinearSolvers/CellEB2/MyTest.cpp @@ -435,6 +435,9 @@ MyTest::initData () 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 const& beb_arr = bcoef_eb[ilev].array(mfi); From b6fa9949b3ec5b2b33969d9f1f03bae9e47cc807 Mon Sep 17 00:00:00 2001 From: Edward Basso Date: Wed, 22 Apr 2026 14:06:42 -0700 Subject: [PATCH 10/11] Rename function to mytest_set_fluxeb_exact --- Tests/LinearSolvers/CellEB2/MyTest.cpp | 3 +-- Tests/LinearSolvers/CellEB2/MyTest_K.H | 2 +- 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/Tests/LinearSolvers/CellEB2/MyTest.cpp b/Tests/LinearSolvers/CellEB2/MyTest.cpp index f77e0c37c0e..64bf3ff9750 100644 --- a/Tests/LinearSolvers/CellEB2/MyTest.cpp +++ b/Tests/LinearSolvers/CellEB2/MyTest.cpp @@ -450,8 +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); - // feb_ex_arr(i,j,k) = 0.0; - mytest_set_fluxeb(i,j,k,feb_ex_arr,flag_arr,ebdata,dx,lprob_type,bx); + mytest_set_fluxeb_exact(i,j,k,feb_ex_arr,flag_arr,ebdata,dx,lprob_type,bx); }); } diff --git a/Tests/LinearSolvers/CellEB2/MyTest_K.H b/Tests/LinearSolvers/CellEB2/MyTest_K.H index 791e6303a75..367513f4447 100644 --- a/Tests/LinearSolvers/CellEB2/MyTest_K.H +++ b/Tests/LinearSolvers/CellEB2/MyTest_K.H @@ -137,7 +137,7 @@ void mytest_set_phi_eb (int i, int j, int k, amrex::Array4 const& p } AMREX_GPU_DEVICE AMREX_FORCE_INLINE -void mytest_set_fluxeb ( +void mytest_set_fluxeb_exact ( int i, int j, int k, amrex::Array4 const& fluxeb_exact, amrex::Array4 const& flag_arr, From fc79fd392d46e822494370e59c0bf78fe83c22cb Mon Sep 17 00:00:00 2001 From: Edward Basso Date: Wed, 22 Apr 2026 14:13:03 -0700 Subject: [PATCH 11/11] Use const for ebdata values --- Tests/LinearSolvers/CellEB2/MyTest_K.H | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/Tests/LinearSolvers/CellEB2/MyTest_K.H b/Tests/LinearSolvers/CellEB2/MyTest_K.H index 367513f4447..fc7584bfc03 100644 --- a/Tests/LinearSolvers/CellEB2/MyTest_K.H +++ b/Tests/LinearSolvers/CellEB2/MyTest_K.H @@ -149,24 +149,24 @@ void mytest_set_fluxeb_exact ( amrex::IntVect iv(AMREX_D_DECL(i,j,k)); if (vbx.contains(iv)) { - amrex::EBCellFlag const flag = flag_arr(i,j,k); + const amrex::EBCellFlag flag = flag_arr(i,j,k); if (flag.isCovered()) { fluxeb_exact(i,j,k) = 0.0; } else if (flag.isSingleValued()) { - amrex::Real const bcentx = ebdata.get(i,j,k,0); - amrex::Real const bcenty = ebdata.get(i,j,k,1); + const amrex::Real bcentx = ebdata.get(i,j,k,0); + const amrex::Real bcenty = ebdata.get(i,j,k,1); // Compute boundary area normal vector - amrex::Real apxm = ebdata.get(i ,j,k); - amrex::Real apxp = ebdata.get(i+1,j,k); + const amrex::Real apxm = ebdata.get(i ,j,k); + const amrex::Real apxp = ebdata.get(i+1,j,k); amrex::Real dapx = (apxm-apxp)/dx[0]; - amrex::Real apym = ebdata.get(i,j ,k); - amrex::Real apyp = ebdata.get(i,j+1,k); + const amrex::Real apym = ebdata.get(i,j ,k); + const amrex::Real apyp = ebdata.get(i,j+1,k); amrex::Real dapy = (apym-apyp)/dx[1]; #if (AMREX_SPACEDIM == 3) - amrex::Real apzm = ebdata.get(i,j,k ); - amrex::Real apzp = ebdata.get(i,j,k+1); + const amrex::Real apzm = ebdata.get(i,j,k ); + const amrex::Real apzp = ebdata.get(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));