From 73c22935ab4e75966f63341e276c770d106f350b Mon Sep 17 00:00:00 2001 From: "Aaron M. Lattanzi" <103702284+AMLattanzi@users.noreply.github.com> Date: Wed, 10 Jan 2024 17:26:38 -0800 Subject: [PATCH] fix outflow bc for boundary normal velocities. (#1374) Co-authored-by: Aaron Lattanzi --- .../BoundaryConditions_xvel.cpp | 26 ++++++++++++++----- .../BoundaryConditions_yvel.cpp | 26 ++++++++++++++----- .../BoundaryConditions_zvel.cpp | 7 ++++- Source/IndexDefines.H | 1 + Source/Initialization/ERF_init_bcs.cpp | 2 ++ 5 files changed, 49 insertions(+), 13 deletions(-) diff --git a/Source/BoundaryConditions/BoundaryConditions_xvel.cpp b/Source/BoundaryConditions/BoundaryConditions_xvel.cpp index 041b068bf..216c70675 100644 --- a/Source/BoundaryConditions/BoundaryConditions_xvel.cpp +++ b/Source/BoundaryConditions/BoundaryConditions_xvel.cpp @@ -62,7 +62,8 @@ void ERFPhysBCFunct::impose_lateral_xvel_bcs (const Array4& 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]; @@ -72,16 +73,23 @@ void ERFPhysBCFunct::impose_lateral_xvel_bcs (const Array4& 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]; @@ -91,12 +99,18 @@ void ERFPhysBCFunct::impose_lateral_xvel_bcs (const Array4& 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 diff --git a/Source/BoundaryConditions/BoundaryConditions_yvel.cpp b/Source/BoundaryConditions/BoundaryConditions_yvel.cpp index a7fd08b80..ebfbde9cf 100644 --- a/Source/BoundaryConditions/BoundaryConditions_yvel.cpp +++ b/Source/BoundaryConditions/BoundaryConditions_yvel.cpp @@ -94,7 +94,8 @@ void ERFPhysBCFunct::impose_lateral_yvel_bcs (const Array4& 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]; @@ -104,16 +105,23 @@ void ERFPhysBCFunct::impose_lateral_yvel_bcs (const Array4& 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]; @@ -123,12 +131,18 @@ void ERFPhysBCFunct::impose_lateral_yvel_bcs (const Array4& 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; + } } ); } diff --git a/Source/BoundaryConditions/BoundaryConditions_zvel.cpp b/Source/BoundaryConditions/BoundaryConditions_zvel.cpp index b41f5c279..ee2f3c8ee 100644 --- a/Source/BoundaryConditions/BoundaryConditions_zvel.cpp +++ b/Source/BoundaryConditions/BoundaryConditions_zvel.cpp @@ -235,7 +235,7 @@ void ERFPhysBCFunct::impose_vertical_zvel_bcs (const Array4& 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 @@ -248,6 +248,11 @@ void ERFPhysBCFunct::impose_vertical_zvel_bcs (const Array4& 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(); } diff --git a/Source/IndexDefines.H b/Source/IndexDefines.H index 8ece30b94..10c894d59 100644 --- a/Source/IndexDefines.H +++ b/Source/IndexDefines.H @@ -133,6 +133,7 @@ enum mathematicalBndryTypes : int { MOST = 101, ext_dir_ingested = 102, neumann = 103, + neumann_int = 104 }; } diff --git a/Source/Initialization/ERF_init_bcs.cpp b/Source/Initialization/ERF_init_bcs.cpp index 4f2413ef5..e7249a4e4 100644 --- a/Source/Initialization/ERF_init_bcs.cpp +++ b/Source/Initialization/ERF_init_bcs.cpp @@ -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)