Skip to content
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.

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

3d anisotropic eb #3907

Merged
merged 21 commits into from Apr 30, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
35 changes: 20 additions & 15 deletions Src/EB/AMReX_EB2_3D_C.cpp
Expand Up @@ -32,18 +32,18 @@ void set_eb_data (const int i, const int j, const int k,
Array4<Real const> const& fcx, Array4<Real const> const& fcy,
Array4<Real const> const& fcz, Array4<Real const> const& m2x,
Array4<Real const> const& m2y, Array4<Real const> const& m2z,
GpuArray<Real,AMREX_SPACEDIM> const& dx,
Array4<Real> const& vfrac, Array4<Real> const& vcent,
Array4<Real> const& barea, Array4<Real> const& bcent,
Array4<Real> const& bnorm, Real small_volfrac,
bool& is_small_cell, bool& is_multicut) noexcept
{
Real axm = apx(i,j,k);
Real axp = apx(i+1,j,k);
Real aym = apy(i,j,k);
Real ayp = apy(i,j+1,k);
Real azm = apz(i,j,k);
Real azp = apz(i,j,k+1);

const Real axm = apx(i,j,k);
const Real axp = apx(i+1,j,k);
const Real aym = apy(i,j,k);
const Real ayp = apy(i,j+1,k);
const Real azm = apz(i,j,k);
const Real azp = apz(i,j,k+1);
// Check for small cell first
if (((axm == 0.0_rt && axp == 0.0_rt) &&
(aym == 0.0_rt && ayp == 0.0_rt) &&
Expand Down Expand Up @@ -76,10 +76,10 @@ void set_eb_data (const int i, const int j, const int k,
return;
}

Real dapx = axm - axp;
Real dapy = aym - ayp;
Real dapz = azm - azp;
Real apnorm = std::sqrt(dapx*dapx+dapy*dapy+dapz*dapz);
Real dapx = (axm - axp)*(dx[1]*dx[2]);
Real dapy = (aym - ayp)*(dx[0]*dx[2]);
Real dapz = (azm - azp)*(dx[0]*dx[1]);
const Real apnorm = std::sqrt(dapx*dapx+dapy*dapy+dapz*dapz) + 1.e-30_rt*std::sqrt(dx[0]*dx[1]*dx[2]);
Copy link
Member

Choose a reason for hiding this comment

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

The second term has a different units than the first. The first term is L^2, whereas the second is L^(3/2).

if (apnorm == 0.0_rt) {
bool maybe_multi_cuts = (axm == 0.0_rt && axp == 0.0_rt) ||
(aym == 0.0_rt && ayp == 0.0_rt) ||
Expand All @@ -96,10 +96,13 @@ void set_eb_data (const int i, const int j, const int k,
Real nx = dapx * apnorminv;
Real ny = dapy * apnorminv;
Real nz = dapz * apnorminv;
const Real bareascaling = std::sqrt( (nx*dx[0])*(nx*dx[0]) +
(ny*dx[1])*(ny*dx[1]) +
(nz*dx[2])*(nz*dx[2]) );
bnorm(i,j,k,0) = nx;
bnorm(i,j,k,1) = ny;
bnorm(i,j,k,2) = nz;
barea(i,j,k) = nx*dapx + ny*dapy + nz*dapz;
barea(i,j,k) = (nx*dapx/(dx[1]*dx[2]) + ny*dapy/(dx[0]*dx[2]) + nz*dapz/(dx[0]*dx[1]));

Real aax = 0.5_rt*(axm+axp);
Real aay = 0.5_rt*(aym+ayp);
Expand All @@ -121,7 +124,7 @@ void set_eb_data (const int i, const int j, const int k,
return;
}

Real bainv = 1.0_rt/barea(i,j,k);
Real bainv = bareascaling*bareascaling/apnorm;
Copy link
Member

Choose a reason for hiding this comment

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

This is the only place bareascaling is used. barescaling is computed with std::sqrt(..). So maybe sqrt is unnecessary.

bcent(i,j,k,0) = bainv * (Bx + nx*vfrac(i,j,k));
bcent(i,j,k,1) = bainv * (By + ny*vfrac(i,j,k));
bcent(i,j,k,2) = bainv * (Bz + nz*vfrac(i,j,k));
Expand Down Expand Up @@ -310,6 +313,7 @@ void set_eb_cell (int i, int j, int k,
Array4<Real const> const& fcx, Array4<Real const> const& fcy,
Array4<Real const> const& fcz, Array4<Real const> const& m2x,
Array4<Real const> const& m2y, Array4<Real const> const& m2z,
GpuArray<Real,AMREX_SPACEDIM> const& dx,
Array4<Real> const& vfrac, Array4<Real> const& vcent,
Array4<Real> const& barea, Array4<Real> const& bcent,
Array4<Real> const& bnorm, Real small_volfrac,
Expand Down Expand Up @@ -341,7 +345,7 @@ void set_eb_cell (int i, int j, int k,
barea(i,j,k) = 0.0_rt;
} else {
set_eb_data(i, j , k, cell, apx, apy, apz, fcx, fcy, fcz, m2x, m2y, m2z,
vfrac, vcent, barea, bcent, bnorm, small_volfrac,
dx, vfrac, vcent, barea, bcent, bnorm, small_volfrac,
is_small_cell, is_multicut);
}
}
Expand Down Expand Up @@ -773,6 +777,7 @@ void build_cells (Box const& bx, Array4<EBCellFlag> const& cell,
Array4<Real const> const& fcx, Array4<Real const> const& fcy,
Array4<Real const> const& fcz, Array4<Real const> const& m2x,
Array4<Real const> const& m2y, Array4<Real const> const& m2z,
GpuArray<Real,AMREX_SPACEDIM> const& dx,
Array4<Real> const& vfrac, Array4<Real> const& vcent,
Array4<Real> const& barea, Array4<Real> const& bcent,
Array4<Real> const& bnorm, Array4<EBCellFlag> const& ctmp,
Expand All @@ -790,7 +795,7 @@ void build_cells (Box const& bx, Array4<EBCellFlag> const& cell,
bool is_small_cell = false;
bool is_multicut = false;
set_eb_cell(i, j, k, cell, apx, apy, apz, fcx, fcy, fcz, m2x, m2y, m2z,
vfrac, vcent, barea, bcent, bnorm, small_volfrac,
dx, vfrac, vcent, barea, bcent, bnorm, small_volfrac,
is_small_cell, is_multicut);
if (is_small_cell) {
Gpu::Atomic::Add(dp, 1);
Expand Down
1 change: 1 addition & 0 deletions Src/EB/AMReX_EB2_C.H
Expand Up @@ -64,6 +64,7 @@ void build_cells (Box const& bx, Array4<EBCellFlag> const& cell,
Array4<Real const> const& fcx, Array4<Real const> const& fcy,
Array4<Real const> const& fcz, Array4<Real const> const& m2x,
Array4<Real const> const& m2y, Array4<Real const> const& m2z,
GpuArray<Real,AMREX_SPACEDIM> const& dx,
Array4<Real> const& vfrac, Array4<Real> const& vcent,
Array4<Real> const& barea, Array4<Real> const& bcent,
Array4<Real> const& bnorm, Array4<EBCellFlag> const& ctmp,
Expand Down
2 changes: 1 addition & 1 deletion Src/EB/AMReX_EB2_Level.H
Expand Up @@ -422,7 +422,7 @@ GShopLevel<G>::define_fine (G const& gshop, const Geometry& geom,
Array4<EBCellFlag> const& cfgtmp = cellflagtmp.array();

build_cells(vbx, cfg, ftx, fty, ftz, apx, apy, apz,
fcx, fcy, fcz, xm2, ym2, zm2, vfr, ctr,
fcx, fcy, fcz, xm2, ym2, zm2, dx, vfr, ctr,
bar, bct, bnm, cfgtmp, lst,
small_volfrac, geom, extend_domain_face, cover_multiple_cuts,
nsm, nmc);
Expand Down