Skip to content

Commit

Permalink
fix outflow bc for boundary normal velocities. (#1374)
Browse files Browse the repository at this point in the history 
Co-authored-by: Aaron Lattanzi <[email protected]>
  • Loading branch information
@AMLattanzi
AMLattanzi and Aaron Lattanzi authored Jan 11, 2024
1 parent 60668b2 commit 73c2293
Show file tree
Hide file tree
Showing 5 changed files with 49 additions and 13 deletions.
26 changes: 20 additions & 6 deletions Source/BoundaryConditions/BoundaryConditions_xvel.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -62,7 +62,8 @@ void ERFPhysBCFunct::impose_lateral_xvel_bcs (const Array4<Real>& dest_arr,
Box bx_xlo_face(bx); bx_xlo_face.setSmall(0,dom_lo.x ); bx_xlo_face.setBig(0,dom_lo.x );
Box bx_xhi_face(bx); bx_xhi_face.setSmall(0,dom_hi.x+1); bx_xhi_face.setBig(0,dom_hi.x+1);
ParallelFor(
bx_xlo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) {
bx_xlo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n)
{
int iflip = dom_lo.x - i;
if (bc_ptr[n].lo(0) == ERFBCType::ext_dir) {
dest_arr(i,j,k) = l_bc_extdir_vals_d[n][0];
Expand All @@ -72,16 +73,23 @@ void ERFPhysBCFunct::impose_lateral_xvel_bcs (const Array4<Real>& dest_arr,
dest_arr(i,j,k) = dest_arr(iflip,j,k);
} else if (bc_ptr[n].lo(0) == ERFBCType::reflect_odd) {
dest_arr(i,j,k) = -dest_arr(iflip,j,k);
} else if (bc_ptr[n].lo(0) == ERFBCType::neumann_int) {
dest_arr(i,j,k) = (4.0*dest_arr(dom_lo.x+1,j,k) - dest_arr(dom_lo.x+2,j,k))/3.0;
}
},
// We only set the values on the domain faces themselves if EXT_DIR
bx_xlo_face, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) {
if (bc_ptr[n].lo(0) == ERFBCType::ext_dir)
bx_xlo_face, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n)
{
if (bc_ptr[n].lo(0) == ERFBCType::ext_dir) {
dest_arr(i,j,k) = l_bc_extdir_vals_d[n][0];
} else if (bc_ptr[n].lo(0) == ERFBCType::neumann_int) {
dest_arr(i,j,k) = (4.0*dest_arr(dom_lo.x+1,j,k) - dest_arr(dom_lo.x+2,j,k))/3.0;
}
}
);
ParallelFor(
bx_xhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) {
bx_xhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n)
{
int iflip = 2*(dom_hi.x + 1) - i;
if (bc_ptr[n].hi(0) == ERFBCType::ext_dir) {
dest_arr(i,j,k) = l_bc_extdir_vals_d[n][3];
Expand All @@ -91,12 +99,18 @@ void ERFPhysBCFunct::impose_lateral_xvel_bcs (const Array4<Real>& dest_arr,
dest_arr(i,j,k) = dest_arr(iflip,j,k);
} else if (bc_ptr[n].hi(0) == ERFBCType::reflect_odd) {
dest_arr(i,j,k) = -dest_arr(iflip,j,k);
} else if (bc_ptr[n].hi(0) == ERFBCType::neumann_int) {
dest_arr(i,j,k) = (4.0*dest_arr(dom_hi.x,j,k) - dest_arr(dom_hi.x-1,j,k))/3.0;
}
},
// We only set the values on the domain faces themselves if EXT_DIR
bx_xhi_face, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) {
if (bc_ptr[n].hi(0) == ERFBCType::ext_dir)
bx_xhi_face, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n)
{
if (bc_ptr[n].hi(0) == ERFBCType::ext_dir) {
dest_arr(i,j,k) = l_bc_extdir_vals_d[n][3];
} else if (bc_ptr[n].hi(0) == ERFBCType::neumann_int) {
dest_arr(i,j,k) = (4.0*dest_arr(dom_hi.x,j,k) - dest_arr(dom_hi.x-1,j,k))/3.0;
}
}
);
} // not periodic in x
Expand Down
26 changes: 20 additions & 6 deletions Source/BoundaryConditions/BoundaryConditions_yvel.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -94,7 +94,8 @@ void ERFPhysBCFunct::impose_lateral_yvel_bcs (const Array4<Real>& dest_arr,
Box bx_yhi_face(bx); bx_yhi_face.setSmall(1,dom_hi.y+1); bx_yhi_face.setBig(1,dom_hi.y+1);

ParallelFor(
bx_ylo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) {
bx_ylo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n)
{
int jflip = dom_lo.y-j;
if (bc_ptr[n].lo(1) == ERFBCType::ext_dir) {
dest_arr(i,j,k) = l_bc_extdir_vals_d[n][1];
Expand All @@ -104,16 +105,23 @@ void ERFPhysBCFunct::impose_lateral_yvel_bcs (const Array4<Real>& dest_arr,
dest_arr(i,j,k) = dest_arr(i,jflip,k);
} else if (bc_ptr[n].lo(1) == ERFBCType::reflect_odd) {
dest_arr(i,j,k) = -dest_arr(i,jflip,k);
} else if (bc_ptr[n].lo(1) == ERFBCType::neumann_int) {
dest_arr(i,j,k) = (4.0*dest_arr(i,dom_lo.y+1,k) - dest_arr(i,dom_lo.y+2,k))/3.0;
}
},
// We only set the values on the domain faces themselves if EXT_DIR
bx_ylo_face, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) {
if (bc_ptr[n].lo(1) == ERFBCType::ext_dir)
bx_ylo_face, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n)
{
if (bc_ptr[n].lo(1) == ERFBCType::ext_dir) {
dest_arr(i,j,k) = l_bc_extdir_vals_d[n][1];
} else if (bc_ptr[n].lo(1) == ERFBCType::neumann_int) {
dest_arr(i,j,k) = (4.0*dest_arr(i,dom_lo.y+1,k) - dest_arr(i,dom_lo.y+2,k))/3.0;
}
}
);
ParallelFor(
bx_yhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) {
bx_yhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n)
{
int jflip = 2*(dom_hi.y + 1) - j;
if (bc_ptr[n].hi(1) == ERFBCType::ext_dir) {
dest_arr(i,j,k) = l_bc_extdir_vals_d[n][4];
Expand All @@ -123,12 +131,18 @@ void ERFPhysBCFunct::impose_lateral_yvel_bcs (const Array4<Real>& dest_arr,
dest_arr(i,j,k) = dest_arr(i,jflip,k);
} else if (bc_ptr[n].hi(1) == ERFBCType::reflect_odd) {
dest_arr(i,j,k) = -dest_arr(i,jflip,k);
} else if (bc_ptr[n].hi(1) == ERFBCType::neumann_int) {
dest_arr(i,j,k) = (4.0*dest_arr(i,dom_hi.y,k) - dest_arr(i,dom_hi.y-1,k))/3.0;
}
},
// We only set the values on the domain faces themselves if EXT_DIR
bx_yhi_face, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) {
if (bc_ptr[n].lo(4) == ERFBCType::ext_dir)
bx_yhi_face, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n)
{
if (bc_ptr[n].hi(1) == ERFBCType::ext_dir) {
dest_arr(i,j,k) = l_bc_extdir_vals_d[n][4];
} else if (bc_ptr[n].hi(1) == ERFBCType::neumann_int) {
dest_arr(i,j,k) = (4.0*dest_arr(i,dom_hi.y,k) - dest_arr(i,dom_hi.y-1,k))/3.0;
}
}
);
}
Expand Down
7 changes: 6 additions & 1 deletion Source/BoundaryConditions/BoundaryConditions_zvel.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -235,7 +235,7 @@ void ERFPhysBCFunct::impose_vertical_zvel_bcs (const Array4<Real>& dest_arr,
// *******************************************************

AMREX_ALWAYS_ASSERT(bc_ptr_w_h[0].hi(2) == ERFBCType::ext_dir ||
bc_ptr_w_h[0].hi(2) == ERFBCType::foextrap);
bc_ptr_w_h[0].hi(2) == ERFBCType::neumann_int);

// NOTE: if we set SlipWall at top, that generates ERFBCType::ext_dir which sets w=0 here
// NOTE: if we set Outflow at top, that generates ERFBCType::foextrap which doesn't touch w here
Expand All @@ -248,6 +248,11 @@ void ERFPhysBCFunct::impose_vertical_zvel_bcs (const Array4<Real>& dest_arr,
dest_arr(i,j,k) = l_bc_extdir_vals_d[0][5];
}
});
} else if (bc_ptr_w_h[0].hi(2) == ERFBCType::neumann_int) {
ParallelFor(makeSlab(bx,2,dom_hi.z+1), [=] AMREX_GPU_DEVICE (int i, int j, int k)
{
dest_arr(i,j,k) = (4.0*dest_arr(i,j,dom_hi.z) - dest_arr(i,j,dom_hi.z-1))/3.0;
});
}
Gpu::streamSynchronize();
}
1 change: 1 addition & 0 deletions Source/IndexDefines.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -133,6 +133,7 @@ enum mathematicalBndryTypes : int {
MOST = 101,
ext_dir_ingested = 102,
neumann = 103,
neumann_int = 104
};
}

Expand Down
2 changes: 2 additions & 0 deletions Source/Initialization/ERF_init_bcs.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -277,9 +277,11 @@ void ERF::init_bcs ()
if (side == Orientation::low) {
for (int i = 0; i < AMREX_SPACEDIM; i++)
domain_bcs_type[BCVars::xvel_bc+i].setLo(dir, ERFBCType::foextrap);
domain_bcs_type[BCVars::xvel_bc+dir].setLo(dir, ERFBCType::neumann_int);
} else {
for (int i = 0; i < AMREX_SPACEDIM; i++)
domain_bcs_type[BCVars::xvel_bc+i].setHi(dir, ERFBCType::foextrap);
domain_bcs_type[BCVars::xvel_bc+dir].setHi(dir, ERFBCType::neumann_int);
}
}
else if (bct == ERF_BC::inflow)
Expand Down

0 comments on commit 73c2293

Please sign in to comment.