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

Fix use real bcs #1381

Merged
merged 2 commits into from
Jan 14, 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.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion Source/Advection/AdvectionSrcForState.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -219,7 +219,7 @@ AdvectionSrcForScalars (const Box& bx, const int icomp, const int ncomp,
}
}

amrex::ParallelFor(bx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
ParallelFor(bx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
{
Real invdetJ = (use_terrain) ? 1. / detJ(i,j,k) : 1.;

Expand Down
5 changes: 1 addition & 4 deletions Source/BoundaryConditions/BoundaryConditions_metgrid.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,7 @@ void
ERF::fill_from_metgrid (const Vector<MultiFab*>& mfs,
const Real time,
bool cons_only,
int icomp_cons,
int /*icomp_cons*/,
int ncomp_cons)
{
int lev = 0;
Expand Down Expand Up @@ -64,9 +64,6 @@ ERF::fill_from_metgrid (const Vector<MultiFab*>& mfs,
const auto& dom_lo = amrex::lbound(domain);
const auto& dom_hi = amrex::ubound(domain);

// Offset only applys to cons (we may fill a subset of these vars)
int offset = (var_idx == Vars::cons) ? icomp_cons : 0;

// Loop over each component
for (int comp_idx(0); comp_idx < comp_var[var_idx]; ++comp_idx)
{
Expand Down
23 changes: 6 additions & 17 deletions Source/BoundaryConditions/ERF_FillPatch.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -100,13 +100,8 @@ ERF::FillPatch (int lev, Real time, const Vector<MultiFab*>& mfs, bool fillset)
if (m_r2d) fill_from_bndryregs(mfs,time);

// We call this even if init_type == real because this routine will fill the vertical bcs
(*physbcs[lev])(mfs,icomp_cons,ncomp_cons,ngvect_cons,ngvect_vels,init_type,cons_only,BCVars::cons_bc,time);

/*
// Update vars in the micro model
if (solverChoice.moisture_type != MoistureType::None)
micro.Update_Micro_Vars_Lev(lev, vars_new[lev][Vars::cons]);
*/
(*physbcs[lev])(mfs,icomp_cons,ncomp_cons,ngvect_cons,ngvect_vels,
use_real_bcs,cons_only,BCVars::cons_bc,time);
}

/*
Expand Down Expand Up @@ -134,8 +129,8 @@ ERF::FillPatchMoistVars (int lev, MultiFab& mf)
IntVect ngvect_cons = mf.nGrowVect();
IntVect ngvect_vels = {0,0,0};

if ((init_type != "real") && (init_type != "metgrid")) {
(*physbcs[lev])({&mf},icomp_cons,ncomp_cons,ngvect_cons,ngvect_vels,init_type,cons_only,bccomp_cons);
if (!use_real_bcs) {
(*physbcs[lev])({&mf},icomp_cons,ncomp_cons,ngvect_cons,ngvect_vels,use_real_bcs,cons_only,bccomp_cons);
}

mf.FillBoundary(geom[lev].periodicity());
Expand Down Expand Up @@ -259,18 +254,12 @@ ERF::FillIntermediatePatch (int lev, Real time,
if (m_r2d) fill_from_bndryregs(mfs,time);

// We call this even if init_type == real because this routine will fill the vertical bcs
(*physbcs[lev])(mfs,icomp_cons,ncomp_cons,ngvect_cons,ngvect_vels,init_type,cons_only,BCVars::cons_bc,time);
(*physbcs[lev])(mfs,icomp_cons,ncomp_cons,ngvect_cons,ngvect_vels,use_real_bcs,cons_only,BCVars::cons_bc,time);
// ***************************************************************************

// MOST boundary conditions
if (!(cons_only && ncomp_cons == 1) && m_most && allow_most_bcs)
m_most->impose_most_bcs(lev,mfs,eddyDiffs_lev[lev].get(),z_phys_nd[lev].get());

/*
// Update vars in the micro model
if (solverChoice.moisture_type != MoistureType::None && (!cons_only && ncomp_cons > 2))
micro.Update_Micro_Vars_Lev(lev, *mfs[Vars::cons]);
*/
}

/*
Expand Down Expand Up @@ -331,7 +320,7 @@ ERF::FillCoarsePatch (int lev, Real time, const Vector<MultiFab*>& mfs)
IntVect ngvect_vels = mfs[Vars::xvel]->nGrowVect();
bool cons_only = false;

(*physbcs[lev])(mfs,0,mfs[Vars::cons]->nComp(),ngvect_cons,ngvect_vels,init_type,cons_only,BCVars::cons_bc,time);
(*physbcs[lev])(mfs,0,mfs[Vars::cons]->nComp(),ngvect_cons,ngvect_vels,use_real_bcs,cons_only,BCVars::cons_bc,time);

// ***************************************************************************
// Since lev > 0 here we don't worry about m_r2d or wrfbdy data
Expand Down
2 changes: 1 addition & 1 deletion Source/BoundaryConditions/ERF_FillPatcher.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -147,7 +147,7 @@ void ERFFillPatcher::BuildMask (BoxArray const& fba,

for (auto const& b : com_bl) {
Box com_bx = vbx & b;
amrex::ParallelFor(com_bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
ParallelFor(com_bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
mask_arr(i,j,k) = mask_val;
});
Expand Down
18 changes: 9 additions & 9 deletions Source/BoundaryConditions/ERF_PhysBCFunct.H
Original file line number Diff line number Diff line change
Expand Up @@ -45,18 +45,18 @@ public:
/*
* Impose physical boundary conditions at domain boundaries
*
* @param[out] mfs Vector of MultiFabs to be filled containing, in order: cons, xvel, yvel, and zvel data
* @param[in] icomp_cons starting component for conserved variables
* @param[in] ncomp_cons number of components for conserved variables
* @param[in] nghost_cons number of ghost cells to be filled for conserved variables
* @param[in] nghost_vels number of ghost cells to be filled for velocity components
* @param[in] time time at which the data should be filled
* @param[in] init_type if "real" then we fill boundary conditions for interior locations
* @param[in] cons_only if 1 then only fill conserved variables
* @param[out] mfs Vector of MultiFabs to be filled containing, in order: cons, xvel, yvel, and zvel data
* @param[in] icomp_cons starting component for conserved variables
* @param[in] ncomp_cons number of components for conserved variables
* @param[in] nghost_cons number of ghost cells to be filled for conserved variables
* @param[in] nghost_vels number of ghost cells to be filled for velocity components
* @param[in] time time at which the data should be filled
* @param[in] use_real_bcs if true then we fill boundary conditions for interior locations
* @param[in] cons_only if 1 then only fill conserved variables
*/
void operator() (const amrex::Vector<amrex::MultiFab*>& mfs, int icomp, int ncomp,
amrex::IntVect const& nghost_cons, amrex::IntVect const& nghost_vels,
std::string& init_type, bool cons_only, int bccomp_cons, const amrex::Real time = 0.0);
bool use_real_bcs, bool cons_only, int bccomp_cons, const amrex::Real time = 0.0);

void impose_lateral_xvel_bcs (const amrex::Array4<amrex::Real>& dest_arr,
const amrex::Box& bx, const amrex::Box& domain,
Expand Down
14 changes: 7 additions & 7 deletions Source/BoundaryConditions/ERF_PhysBCFunct.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -12,18 +12,18 @@ using namespace amrex;
* @param[in] ncomp_cons number of components for conserved variables
* @param[in] nghost_cons number of ghost cells to be filled for conserved variables
* @param[in] nghost_vels number of ghost cells to be filled for velocity components
* @param[in] init_type if "real" then we fill boundary conditions for interior locations
* @param[in] cons_only if 1 then only fill conserved variables
* @param[in] use_real_bcs we fill boundary conditions for interior locations
* @param[in] cons_only if 1 then only fill conserved variables
*/

void ERFPhysBCFunct::operator() (const Vector<MultiFab*>& mfs, int icomp_cons, int ncomp_cons,
IntVect const& nghost_cons, IntVect const& nghost_vels,
std::string& init_type, bool cons_only, int bccomp_cons, const Real time)
bool use_real_bcs, bool cons_only, int bccomp_cons, const Real time)
{
BL_PROFILE("ERFPhysBCFunct::()");

#ifndef ERF_USE_NETCDF
amrex::ignore_unused(init_type);
amrex::ignore_unused(use_real_bcs);
#endif

if (m_geom.isAllPeriodic()) return;
Expand Down Expand Up @@ -82,7 +82,7 @@ void ERFPhysBCFunct::operator() (const Vector<MultiFab*>& mfs, int icomp_cons, i
z_nd_arr = m_z_phys_nd->const_array(mfi);
}

if ((init_type != "real") && (init_type != "metgrid"))
if (!use_real_bcs)
{

//! If there are cells not in the valid + periodic grown box
Expand All @@ -109,10 +109,10 @@ void ERFPhysBCFunct::operator() (const Vector<MultiFab*>& mfs, int icomp_cons, i

impose_lateral_zvel_bcs(velz_arr,velx_arr,vely_arr,zbx,domain,z_nd_arr,dxInv,BCVars::zvel_bc);
} // !cons_only
} // init_type != "real" && init_type != "metgrid"
} // use_real_bcs

// Every grid touches the bottom and top domain boundary so we call the vertical bcs
// for every box -- and we need to call these even if init_type == real
// for every box -- and we need to call these even if use_real_bcs
{
impose_vertical_cons_bcs(cons_arr,cbx,domain,z_nd_arr,dxInv,
icomp_cons,ncomp_cons,bccomp_cons);
Expand Down
6 changes: 3 additions & 3 deletions Source/Derive.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,7 @@ void erf_derrhodivide (const amrex::Box& bx,
auto const dat = datfab.array();
auto primitive = derfab.array();

amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
{
const amrex::Real rho = dat(i, j, k, Rho_comp);
const amrex::Real conserved = dat(i, j, k, scalar_index);
Expand Down Expand Up @@ -69,7 +69,7 @@ erf_dersoundspeed (const amrex::Box& bx,
// NOTE: we compute the soundspeed of dry air -- we do not account for any moisture effects here
amrex::Real qv = 0.;

amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
{
const amrex::Real rhotheta = dat(i, j, k, RhoTheta_comp);
const amrex::Real rho = dat(i, j, k, Rho_comp);
Expand Down Expand Up @@ -99,7 +99,7 @@ erf_dertemp (const amrex::Box& bx,
auto const dat = datfab.array();
auto tfab = derfab.array();

amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
{
const amrex::Real rho = dat(i, j, k, Rho_comp);
const amrex::Real rhotheta = dat(i, j, k, RhoTheta_comp);
Expand Down
28 changes: 14 additions & 14 deletions Source/Diffusion/ComputeStrain_N.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -70,7 +70,7 @@ ComputeStrain_N (Box bxcc, Box tbxxy, Box tbxxz, Box tbxyz,
if (xl_v_dir) {
Box planexy = tbxxy; planexy.setBig(0, planexy.smallEnd(0) );
tbxxy.growLo(0,-1);
amrex::ParallelFor(planexy,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
ParallelFor(planexy,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
tau12(i,j,k) = 0.5 * ( (u(i, j, k)/mf_u(i,j,0) - u(i, j-1, k)/mf_u(i,j-1,0))*dxInv[1] +
(-(8./3.) * v(i-1,j,k)/mf_v(i-1,j,0) + 3. * v(i,j,k)/mf_v(i,j,0) - (1./3.) * v(i+1,j,k)/mf_v(i+1,j,0))*dxInv[0] ) * mf_u(i,j,0)*mf_u(i,j,0);
});
Expand All @@ -79,7 +79,7 @@ ComputeStrain_N (Box bxcc, Box tbxxy, Box tbxxz, Box tbxyz,
// note: tilebox xy should be nodal, so i|i-1|i-2 at the bigEnd is analogous to i-1|i|i+1 at the smallEnd
Box planexy = tbxxy; planexy.setSmall(0, planexy.bigEnd(0) );
tbxxy.growHi(0,-1);
amrex::ParallelFor(planexy,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
ParallelFor(planexy,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
tau12(i,j,k) = 0.5 * ( (u(i, j, k)/mf_u(i,j,0) - u(i, j-1, k)/mf_u(i,j-1,0))*dxInv[1] +
-(-(8./3.) * v(i,j,k)/mf_v(i,j,0) + 3. * v(i-1,j,k)/mf_v(i-1,j,0) - (1./3.) * v(i-2,j,k)/mf_v(i-2,j,0))*dxInv[0] ) * mf_u(i,j,0)*mf_u(i,j,0);
});
Expand All @@ -88,7 +88,7 @@ ComputeStrain_N (Box bxcc, Box tbxxy, Box tbxxz, Box tbxyz,
if (xl_w_dir) {
Box planexz = tbxxz; planexz.setBig(0, planexz.smallEnd(0) );
tbxxz.growLo(0,-1);
amrex::ParallelFor(planexz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
ParallelFor(planexz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
tau13(i,j,k) = 0.5 * ( (u(i, j, k) - u(i, j, k-1))*dxInv[2] +
(-(8./3.) * w(i-1,j,k) + 3. * w(i,j,k) - (1./3.) * w(i+1,j,k))*dxInv[0]*mf_u(i,j,0) );
});
Expand All @@ -97,7 +97,7 @@ ComputeStrain_N (Box bxcc, Box tbxxy, Box tbxxz, Box tbxyz,
// note: tilebox xz should be nodal, so i|i-1|i-2 at the bigEnd is analogous to i-1|i|i+1 at the smallEnd
Box planexz = tbxxz; planexz.setSmall(0, planexz.bigEnd(0) );
tbxxz.growHi(0,-1);
amrex::ParallelFor(planexz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
ParallelFor(planexz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
tau13(i,j,k) = 0.5 * ( (u(i, j, k) - u(i, j, k-1))*dxInv[2] +
-(-(8./3.) * w(i,j,k) + 3. * w(i-1,j,k) - (1./3.) * w(i-2,j,k))*dxInv[0]*mf_u(i,j,0) );
});
Expand All @@ -108,7 +108,7 @@ ComputeStrain_N (Box bxcc, Box tbxxy, Box tbxxz, Box tbxyz,
if (yl_u_dir) {
Box planexy = tbxxy; planexy.setBig(1, planexy.smallEnd(1) );
tbxxy.growLo(1,-1);
amrex::ParallelFor(planexy,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
ParallelFor(planexy,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
tau12(i,j,k) = 0.5 * ( (-(8./3.) * u(i,j-1,k)/mf_u(i,j-1,0) + 3. * u(i,j,k)/mf_u(i,j,0) - (1./3.) * u(i,j+1,k)/mf_m(i,j+1,0))*dxInv[1] +
(v(i, j, k)/mf_v(i,j,0) - v(i-1, j, k)/mf_v(i-1,j,0))*dxInv[0] ) * mf_u(i,j,0)*mf_u(i,j,0);
});
Expand All @@ -117,7 +117,7 @@ ComputeStrain_N (Box bxcc, Box tbxxy, Box tbxxz, Box tbxyz,
// note: tilebox xy should be nodal, so j|j-1|j-2 at the bigEnd is analogous to j-1|j|j+1 at the smallEnd
Box planexy = tbxxy; planexy.setSmall(1, planexy.bigEnd(1) );
tbxxy.growHi(1,-1);
amrex::ParallelFor(planexy,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
ParallelFor(planexy,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
tau12(i,j,k) = 0.5 * ( -(-(8./3.) * u(i,j,k)/mf_u(i,j,0) + 3. * u(i,j-1,k)/mf_u(i,j-1,0) - (1./3.) * u(i,j-2,k)/mf_u(i,j-2,0))*dxInv[1] +
(v(i, j, k)/mf_v(i,j,0) - v(i-1, j, k)/mf_v(i-1,j,0))*dxInv[0] ) * mf_u(i,j,0)*mf_u(i,j,0);
});
Expand All @@ -126,7 +126,7 @@ ComputeStrain_N (Box bxcc, Box tbxxy, Box tbxxz, Box tbxyz,
if (yl_w_dir) {
Box planeyz = tbxyz; planeyz.setBig(1, planeyz.smallEnd(1) );
tbxyz.growLo(1,-1);
amrex::ParallelFor(planeyz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
ParallelFor(planeyz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
tau23(i,j,k) = 0.5 * ( (v(i, j, k) - v(i, j, k-1))*dxInv[2] +
(-(8./3.) * w(i,j-1,k) + 3. * w(i,j ,k) - (1./3.) * w(i,j+1,k))*dxInv[1]*mf_v(i,j,0)*mf_v(i,j,0) );
});
Expand All @@ -135,7 +135,7 @@ ComputeStrain_N (Box bxcc, Box tbxxy, Box tbxxz, Box tbxyz,
// note: tilebox yz should be nodal, so j|j-1|j-2 at the bigEnd is analogous to j-1|j|j+1 at the smallEnd
Box planeyz = tbxyz; planeyz.setSmall(1, planeyz.bigEnd(1) );
tbxyz.growHi(1,-1);
amrex::ParallelFor(planeyz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
ParallelFor(planeyz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
tau23(i,j,k) = 0.5 * ( (v(i, j, k) - v(i, j, k-1))*dxInv[2] +
-(-(8./3.) * w(i,j ,k) + 3. * w(i,j-1,k) - (1./3.) * w(i,j-2,k))*dxInv[1]*mf_v(i,j,0)*mf_v(i,j,0) );
});
Expand All @@ -146,7 +146,7 @@ ComputeStrain_N (Box bxcc, Box tbxxy, Box tbxxz, Box tbxyz,
if (zl_u_dir) {
Box planexz = tbxxz; planexz.setBig(2, planexz.smallEnd(2) );
tbxxz.growLo(2,-1);
amrex::ParallelFor(planexz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
ParallelFor(planexz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
tau13(i,j,k) = 0.5 * ( (-(8./3.) * u(i,j,k-1) + 3. * u(i,j,k) - (1./3.) * u(i,j,k+1))*dxInv[2] +
(w(i, j, k) - w(i-1, j, k))*dxInv[0]*mf_u(i,j,0) );
});
Expand All @@ -155,7 +155,7 @@ ComputeStrain_N (Box bxcc, Box tbxxy, Box tbxxz, Box tbxyz,
// note: tilebox xz should be nodal, so k|k-1|k-2 at the bigEnd is analogous to k-1|k|k+1 at the smallEnd
Box planexz = tbxxz; planexz.setSmall(2, planexz.bigEnd(2) );
tbxxz.growHi(2,-1);
amrex::ParallelFor(planexz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
ParallelFor(planexz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
tau13(i,j,k) = 0.5 * ( -(-(8./3.) * u(i,j,k) + 3. * u(i,j,k-1) - (1./3.) * u(i,j,k-2))*dxInv[2] +
(w(i, j, k) - w(i-1, j, k))*dxInv[0]*mf_u(i,j,0) );
});
Expand All @@ -164,7 +164,7 @@ ComputeStrain_N (Box bxcc, Box tbxxy, Box tbxxz, Box tbxyz,
if (zl_v_dir) {
Box planeyz = tbxyz; planeyz.setBig(2, planeyz.smallEnd(2) );
tbxyz.growLo(2,-1);
amrex::ParallelFor(planeyz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
ParallelFor(planeyz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
tau23(i,j,k) = 0.5 * ( (-(8./3.) * v(i,j,k-1) + 3. * v(i,j,k ) - (1./3.) * v(i,j,k+1))*dxInv[2] +
(w(i, j, k) - w(i, j-1, k))*dxInv[1]*mf_v(i,j,0) );
});
Expand All @@ -173,7 +173,7 @@ ComputeStrain_N (Box bxcc, Box tbxxy, Box tbxxz, Box tbxyz,
// note: tilebox yz should be nodal, so k|k-1|k-2 at the bigEnd is analogous to k-1|k|k+1 at the smallEnd
Box planeyz = tbxyz; planeyz.setSmall(2, planeyz.bigEnd(2) );
tbxyz.growHi(2,-1);
amrex::ParallelFor(planeyz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
ParallelFor(planeyz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
tau23(i,j,k) = 0.5 * ( -(-(8./3.) * v(i,j,k ) + 3. * v(i,j,k-1) - (1./3.) * v(i,j,k-2))*dxInv[2] +
(w(i, j, k) - w(i, j-1, k))*dxInv[1]*mf_v(i,j,0) );
});
Expand All @@ -182,14 +182,14 @@ ComputeStrain_N (Box bxcc, Box tbxxy, Box tbxxz, Box tbxyz,
// Fill the remaining cells
//***********************************************************************************
// Cell centered strains
amrex::ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
tau11(i,j,k) = (u(i+1, j , k )/mf_u(i+1,j,0) - u(i, j, k)/mf_u(i,j,0))*dxInv[0]*mf_u(i,j,0)*mf_u(i,j,0);
tau22(i,j,k) = (v(i , j+1, k )/mf_v(i,j+1,0) - v(i, j, k)/mf_v(i,j,0))*dxInv[1]*mf_v(i,j,0)*mf_v(i,j,0);
tau33(i,j,k) = (w(i , j , k+1) - w(i, j, k))*dxInv[2];
});

// Off-diagonal strains
amrex::ParallelFor(tbxxy,tbxxz,tbxyz,
ParallelFor(tbxxy,tbxxz,tbxyz,
[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
tau12(i,j,k) = 0.5 * ( (u(i, j, k)/mf_u(i,j,0) - u(i, j-1, k)/mf_u(i,j-1,0))*dxInv[1] +
(v(i, j, k)/mf_v(i,j,0) - v(i-1, j, k)/mf_v(i-1,j,0))*dxInv[0] ) * mf_u(i,j,0)*mf_u(i,j,0);
Expand Down
Loading
Loading