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.

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

WIP: interpolate perturbations not full state #1911

Closed
wants to merge 6 commits into from
Closed

WIP: interpolate perturbations not full state #1911

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
4953735
WIP: interpolate perturbations not full state
asalmgren Oct 26, 2024
8603c24
fix order
asalmgren Oct 26, 2024
2bf257c
need to pass old and new base state to FillPatch for lev > 0
asalmgren Oct 26, 2024
a279392
more places to modify FillPatch call
asalmgren Oct 26, 2024
1107e4a
fix typo
asalmgren Oct 26, 2024
f26bcd7
remove extra swap
asalmgren Oct 26, 2024
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
126 changes: 122 additions & 4 deletions Source/BoundaryConditions/ERF_BoundaryConditions_basestate.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,5 +1,6 @@
#include "AMReX_PhysBCFunct.H"
#include <ERF_PhysBCFunct.H>
#include <ERF_Constants.H>

using namespace amrex;

Expand All @@ -15,9 +16,14 @@ void ERFPhysBCFunct_base::impose_lateral_basestate_bcs (const Array4<Real>& dest
int ncomp, const IntVect& nghost)
{
BL_PROFILE_VAR("impose_lateral_base_bcs()",impose_lateral_base_bcs);
//
// Note that the "bx" that comes in here has already been grown in the lateral directions
// but not in the vertical
//

const int* bxlo = bx.loVect();
const int* bxhi = bx.hiVect();

const int* dlo = domain.loVect();
const int* dhi = domain.hiVect();

Expand Down Expand Up @@ -63,6 +69,84 @@ void ERFPhysBCFunct_base::impose_lateral_basestate_bcs (const Array4<Real>& dest
// Populate ghost cells on lo-x and hi-x domain boundaries
Box bx_xlo(bx); bx_xlo.setBig (0,dom_lo.x-nghost[0]);
Box bx_xhi(bx); bx_xhi.setSmall(0,dom_hi.x+nghost[0]);

ParallelFor(
bx_xlo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n)
{
int dest_comp = n;
int l_bc_type = bc_ptr[n].lo(0);
int iflip = dom_lo.x - 1 - i;
if (l_bc_type == ERFBCType::foextrap) {
dest_arr(i,j,k,dest_comp) = dest_arr(dom_lo.x,j,k,dest_comp);
} else if (l_bc_type == ERFBCType::open) {
dest_arr(i,j,k,dest_comp) = dest_arr(dom_lo.x,j,k,dest_comp);
} else if (l_bc_type == ERFBCType::reflect_even) {
dest_arr(i,j,k,dest_comp) = dest_arr(iflip,j,k,dest_comp);
}
},
bx_xhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n)
{
int dest_comp = n;
int h_bc_type = bc_ptr[n].hi(0);
int iflip = 2*dom_hi.x + 1 - i;
if (h_bc_type == ERFBCType::foextrap) {
dest_arr(i,j,k,dest_comp) = dest_arr(dom_hi.x,j,k,dest_comp);
} else if (h_bc_type == ERFBCType::open) {
dest_arr(i,j,k,dest_comp) = dest_arr(dom_hi.x,j,k,dest_comp);
} else if (h_bc_type == ERFBCType::reflect_even) {
dest_arr(i,j,k,dest_comp) = dest_arr(iflip,j,k,dest_comp);
}
}
);
}

if (!is_periodic_in_y)
{
// Populate ghost cells on lo-y and hi-y domain boundaries
Box bx_ylo(bx); bx_ylo.setBig (1,dom_lo.y-nghost[1]);
Box bx_yhi(bx); bx_yhi.setSmall(1,dom_hi.y+nghost[1]);
if (bx_ylo.smallEnd(2) != domain.smallEnd(2)) bx_ylo.growLo(2,nghost[2]);
if (bx_ylo.bigEnd(2) != domain.bigEnd(2)) bx_ylo.growHi(2,nghost[2]);
if (bx_yhi.smallEnd(2) != domain.smallEnd(2)) bx_yhi.growLo(2,nghost[2]);
if (bx_yhi.bigEnd(2) != domain.bigEnd(2)) bx_yhi.growHi(2,nghost[2]);
ParallelFor(
bx_ylo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n)
{
int dest_comp = n;
int l_bc_type = bc_ptr[n].lo(1);
int jflip = dom_lo.y - 1 - j;
if (l_bc_type == ERFBCType::foextrap) {
dest_arr(i,j,k,dest_comp) = dest_arr(i,dom_lo.y,k,dest_comp);
} else if (l_bc_type == ERFBCType::open) {
dest_arr(i,j,k,dest_comp) = dest_arr(i,dom_lo.y,k,dest_comp);
} else if (l_bc_type == ERFBCType::reflect_even) {
dest_arr(i,j,k,dest_comp) = dest_arr(i,jflip,k,dest_comp);
}

},
bx_yhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n)
{
int dest_comp = n;
int h_bc_type = bc_ptr[n].hi(1);
int jflip = 2*dom_hi.y + 1 - j;
if (h_bc_type == ERFBCType::foextrap) {
dest_arr(i,j,k,dest_comp) = dest_arr(i,dom_hi.y,k,dest_comp);
} else if (h_bc_type == ERFBCType::open) {
dest_arr(i,j,k,dest_comp) = dest_arr(i,dom_hi.y,k,dest_comp);
} else if (h_bc_type == ERFBCType::reflect_even) {
dest_arr(i,j,k,dest_comp) = dest_arr(i,jflip,k,dest_comp);
}
}
);
}

// Next do ghost cells in x-direction but not reaching out in y
// The corners we miss here will be covered in the y-loop below or by periodicity
if (!is_periodic_in_x)
{
// Populate ghost cells on lo-x and hi-x domain boundaries
Box bx_xlo(bx); bx_xlo.setBig (0,dom_lo.x-1);
Box bx_xhi(bx); bx_xhi.setSmall(0,dom_hi.x+1);
if (bx_xlo.smallEnd(2) != domain.smallEnd(2)) bx_xlo.growLo(2,nghost[2]);
if (bx_xlo.bigEnd(2) != domain.bigEnd(2)) bx_xlo.growHi(2,nghost[2]);
if (bx_xhi.smallEnd(2) != domain.smallEnd(2)) bx_xhi.growLo(2,nghost[2]);
Expand All @@ -79,6 +163,11 @@ void ERFPhysBCFunct_base::impose_lateral_basestate_bcs (const Array4<Real>& dest
dest_arr(i,j,k,dest_comp) = dest_arr(dom_lo.x,j,k,dest_comp);
} else if (l_bc_type == ERFBCType::reflect_even) {
dest_arr(i,j,k,dest_comp) = dest_arr(iflip,j,k,dest_comp);
} else if (l_bc_type == ERFBCType::reflect_odd) {
dest_arr(i,j,k,dest_comp) = -dest_arr(iflip,j,k,dest_comp);
} else if (l_bc_type == ERFBCType::hoextrapcc) {
Real delta_i = (dom_lo.x - i);
dest_arr(i,j,k,dest_comp) = (1.0 + delta_i)*dest_arr(dom_lo.x,j,k,dest_comp) - delta_i*dest_arr(dom_lo.x+1,j,k,dest_comp) ;
}
},
bx_xhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n)
Expand All @@ -92,6 +181,11 @@ void ERFPhysBCFunct_base::impose_lateral_basestate_bcs (const Array4<Real>& dest
dest_arr(i,j,k,dest_comp) = dest_arr(dom_hi.x,j,k,dest_comp);
} else if (h_bc_type == ERFBCType::reflect_even) {
dest_arr(i,j,k,dest_comp) = dest_arr(iflip,j,k,dest_comp);
} else if (h_bc_type == ERFBCType::reflect_odd) {
dest_arr(i,j,k,dest_comp) = -dest_arr(iflip,j,k,dest_comp);
} else if (h_bc_type == ERFBCType::hoextrapcc) {
Real delta_i = (i - dom_hi.x);
dest_arr(i,j,k,dest_comp) = (1.0 + delta_i)*dest_arr(dom_hi.x,j,k,dest_comp) - delta_i*dest_arr(dom_hi.x-1,j,k,dest_comp) ;
}
}
);
Expand All @@ -100,8 +194,8 @@ void ERFPhysBCFunct_base::impose_lateral_basestate_bcs (const Array4<Real>& dest
if (!is_periodic_in_y)
{
// Populate ghost cells on lo-y and hi-y domain boundaries
Box bx_ylo(bx); bx_ylo.setBig (1,dom_lo.y-nghost[1]);
Box bx_yhi(bx); bx_yhi.setSmall(1,dom_hi.y+nghost[1]);
Box bx_ylo(bx); bx_ylo.setBig (1,dom_lo.y-1);
Box bx_yhi(bx); bx_yhi.setSmall(1,dom_hi.y+1);
if (bx_ylo.smallEnd(2) != domain.smallEnd(2)) bx_ylo.growLo(2,nghost[2]);
if (bx_ylo.bigEnd(2) != domain.bigEnd(2)) bx_ylo.growHi(2,nghost[2]);
if (bx_yhi.smallEnd(2) != domain.smallEnd(2)) bx_yhi.growLo(2,nghost[2]);
Expand All @@ -118,6 +212,11 @@ void ERFPhysBCFunct_base::impose_lateral_basestate_bcs (const Array4<Real>& dest
dest_arr(i,j,k,dest_comp) = dest_arr(i,dom_lo.y,k,dest_comp);
} else if (l_bc_type == ERFBCType::reflect_even) {
dest_arr(i,j,k,dest_comp) = dest_arr(i,jflip,k,dest_comp);
} else if (l_bc_type == ERFBCType::reflect_odd) {
dest_arr(i,j,k,dest_comp) = -dest_arr(i,jflip,k,dest_comp);
} else if (l_bc_type == ERFBCType::hoextrapcc) {
Real delta_j = (dom_lo.y - j);
dest_arr(i,j,k,dest_comp) = (1.0 + delta_j)*dest_arr(i,dom_lo.y,k,dest_comp) - delta_j*dest_arr(i,dom_lo.y+1,k,dest_comp) ;
}

},
Expand All @@ -132,6 +231,11 @@ void ERFPhysBCFunct_base::impose_lateral_basestate_bcs (const Array4<Real>& dest
dest_arr(i,j,k,dest_comp) = dest_arr(i,dom_hi.y,k,dest_comp);
} else if (h_bc_type == ERFBCType::reflect_even) {
dest_arr(i,j,k,dest_comp) = dest_arr(i,jflip,k,dest_comp);
} else if (h_bc_type == ERFBCType::reflect_odd) {
dest_arr(i,j,k,dest_comp) = -dest_arr(i,jflip,k,dest_comp);
} else if (h_bc_type == ERFBCType::hoextrapcc) {
Real delta_j = (j - dom_hi.y);
dest_arr(i,j,k,dest_comp) = (1.0 + delta_j)*dest_arr(i,dom_hi.y,k,dest_comp) - delta_j*dest_arr(i,dom_hi.y-1,k,dest_comp);
}
}
);
Expand All @@ -147,16 +251,30 @@ void ERFPhysBCFunct_base::impose_vertical_basestate_bcs (const Array4<Real>& des
const auto& dom_lo = lbound(domain);
const auto& dom_hi = ubound(domain);

const Real hz = Real(0.5) * m_geom.CellSize(2);

Box bx_zlo1(bx); bx_zlo1.setBig(2,dom_lo.z-1); if (bx_zlo1.ok()) bx_zlo1.setSmall(2,dom_lo.z-1);
ParallelFor(
bx_zlo1, [=] AMREX_GPU_DEVICE (int i, int j, int k)
{
dest_arr(i,j,k,BaseState::r0_comp) = dest_arr(i,j,dom_lo.z,BaseState::r0_comp);
dest_arr(i,j,k,BaseState::p0_comp) = p_0 -
dest_arr(i,j,k,BaseState::r0_comp) * hz * 9.81;
dest_arr(i,j,k,BaseState::pi0_comp) = dest_arr(i,j,dom_lo.z,BaseState::pi0_comp);
dest_arr(i,j,k,BaseState::th0_comp) = dest_arr(i,j,dom_lo.z,BaseState::th0_comp);
}
);

Box bx_zlo(bx); bx_zlo.setBig(2,dom_lo.z-2);
Box bx_zhi(bx); bx_zhi.setSmall(2,dom_hi.z+2);
Box bx_zhi(bx); bx_zhi.setSmall(2,dom_hi.z+1);
ParallelFor(
bx_zlo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n)
{
dest_arr(i,j,k,n) = dest_arr(i,j,dom_lo.z-1,n);
},
bx_zhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n)
{
dest_arr(i,j,k,n) = dest_arr(i,j,dom_hi.z+1,n);
dest_arr(i,j,k,n) = dest_arr(i,j,dom_hi.z,n);
}
);
}
19 changes: 14 additions & 5 deletions Source/BoundaryConditions/ERF_FillCoarsePatch.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -29,11 +29,20 @@ ERF::FillCoarsePatch (int lev, Real time)
//****************************************************************************************************************
//
bool cons_only = false;
FillPatch(lev-1, time, {&vars_new[lev-1][Vars::cons], &vars_new[lev-1][Vars::xvel],
&vars_new[lev-1][Vars::yvel], &vars_new[lev-1][Vars::zvel]},
{&vars_new[lev-1][Vars::cons],
&rU_new[lev-1], &rV_new[lev-1], &rW_new[lev-1]},
false, cons_only);
if (lev == 1) {
FillPatch(lev-1, time, {&vars_new[lev-1][Vars::cons], &vars_new[lev-1][Vars::xvel],
&vars_new[lev-1][Vars::yvel], &vars_new[lev-1][Vars::zvel]},
{&vars_new[lev-1][Vars::cons],
&rU_new[lev-1], &rV_new[lev-1], &rW_new[lev-1]},
cons_only);
} else {
FillPatch(lev-1, time, {&vars_new[lev-1][Vars::cons], &vars_new[lev-1][Vars::xvel],
&vars_new[lev-1][Vars::yvel], &vars_new[lev-1][Vars::zvel]},
{&vars_new[lev-1][Vars::cons],
&rU_new[lev-1], &rV_new[lev-1], &rW_new[lev-1]},
base_state[lev-1], base_state[lev-1],
false, cons_only);
}

//
// ************************************************
Expand Down
51 changes: 40 additions & 11 deletions Source/BoundaryConditions/ERF_FillIntermediatePatch.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -65,14 +65,17 @@ ERF::FillIntermediatePatch (int lev, Real time,
AMREX_ALWAYS_ASSERT(mfs_vel.size() == Vars::NumTypes);

// Enforce no penetration for thin immersed body
if (xflux_imask[lev]) {
ApplyMask(*mfs_mom[IntVars::xmom], *xflux_imask[lev]);
}
if (yflux_imask[lev]) {
ApplyMask(*mfs_mom[IntVars::ymom], *yflux_imask[lev]);
}
if (zflux_imask[lev]) {
ApplyMask(*mfs_mom[IntVars::zmom], *zflux_imask[lev]);
if (!cons_only) {
// Enforce no penetration for thin immersed body
if (xflux_imask[lev]) {
ApplyMask(*mfs_mom[IntVars::xmom], *xflux_imask[lev]);
}
if (yflux_imask[lev]) {
ApplyMask(*mfs_mom[IntVars::ymom], *yflux_imask[lev]);
}
if (zflux_imask[lev]) {
ApplyMask(*mfs_mom[IntVars::zmom], *zflux_imask[lev]);
}
}

// We always come in to this call with updated momenta but we need to create updated velocity
Expand Down Expand Up @@ -104,8 +107,18 @@ ERF::FillIntermediatePatch (int lev, Real time,
Vector<Real> ftime = {time,time};

// Impose physical bc's on coarse data (note time and 0 are not used)
(*physbcs_cons[lev-1])(vars_old[lev-1][Vars::cons],0,ncomp_cons,IntVect{ng_cons},time,BCVars::cons_bc,true);
(*physbcs_cons[lev-1])(vars_new[lev-1][Vars::cons],0,ncomp_cons,IntVect{ng_cons},time,BCVars::cons_bc,true);
(*physbcs_cons[lev-1])(vars_old[lev-1][Vars::cons],icomp_cons,ncomp_cons,IntVect{ng_cons},time,BCVars::cons_bc,true);
(*physbcs_cons[lev-1])(vars_new[lev-1][Vars::cons],icomp_cons,ncomp_cons,IntVect{ng_cons},time,BCVars::cons_bc,true);

// Subtract rho_0 from rho before we interpolate -- note we only subtract
// on valid region of mf since the ghost cells will be filled below
if (icomp_cons == 0) {
MultiFab::Subtract(mf,base_state[lev],BaseState::r0_comp,Rho_comp,1,IntVect{0});
MultiFab::Subtract(vars_old[lev-1][Vars::cons], base_state[lev-1],
BaseState::r0_comp,Rho_comp,1,vars_old[lev-1][Vars::cons].nGrowVect());
MultiFab::Subtract(vars_new[lev-1][Vars::cons], base_state[lev-1],
BaseState::r0_comp,Rho_comp,1,vars_new[lev-1][Vars::cons].nGrowVect());
}

// Call FillPatchTwoLevels which ASSUMES that all ghost cells have already been filled
mapper = &cell_cons_interp;
Expand All @@ -115,6 +128,22 @@ ERF::FillIntermediatePatch (int lev, Real time,
refRatio(lev-1), mapper, domain_bcs_type,
icomp_cons);

if (icomp_cons == 0) {
// Restore the coarse values to what they were
MultiFab::Add(vars_old[lev-1][Vars::cons], base_state[lev-1],
BaseState::r0_comp,Rho_comp,1,vars_old[lev-1][Vars::cons].nGrowVect());
MultiFab::Add(vars_new[lev-1][Vars::cons], base_state[lev-1],
BaseState::r0_comp,Rho_comp,1,vars_new[lev-1][Vars::cons].nGrowVect());

// Set values in the cells outside the domain boundary so that we can do the Add
// without worrying about uninitialized values outside the domain -- these
// will be filled in the physbcs call
mf.setDomainBndry(1.234e20,0,1,geom[lev]);

// Add rho_0 back to rho after we interpolate -- on all the valid + ghost region
MultiFab::Add(mf, base_state[lev],BaseState::r0_comp,Rho_comp,1,IntVect{ng_cons});
}

// *****************************************************************************************

if (!cons_only)
Expand Down Expand Up @@ -201,7 +230,7 @@ ERF::FillIntermediatePatch (int lev, Real time,
#ifdef ERF_USE_NETCDF
// We call this here because it is an ERF routine
if (use_real_bcs && (lev==0)) {
fill_from_realbdy(mfs_vel,time,cons_only,icomp_cons,ncomp_cons,ngvect_cons, ngvect_vels);
fill_from_realbdy(mfs_vel,time,cons_only,icomp_cons,ncomp_cons,ngvect_cons,ngvect_vels);
do_fb = false;
}
#endif
Expand Down
Loading
Loading