diff --git a/Source/BoundaryConditions/BoundaryConditions_netcdf.cpp b/Source/BoundaryConditions/BoundaryConditions_netcdf.cpp index e294168..db371e9 100644 --- a/Source/BoundaryConditions/BoundaryConditions_netcdf.cpp +++ b/Source/BoundaryConditions/BoundaryConditions_netcdf.cpp @@ -79,7 +79,6 @@ REMORA::fill_from_bdyfiles (MultiFab& mf_to_fill, const MultiFab& mf_mask, const const bool apply_south = domain_bcs_type[bccomp+icomp].lo(1) == REMORABCType::clamped; const bool apply_north = domain_bcs_type[bccomp+icomp].hi(1) == REMORABCType::clamped; - #ifdef AMREX_USE_OMP #pragma omp parallel if (Gpu::notInLaunchRegion()) #endif @@ -139,25 +138,26 @@ REMORA::fill_from_bdyfiles (MultiFab& mf_to_fill, const MultiFab& mf_mask, const }); } - if (!xlo_ylo.isEmpty() && apply_east && apply_south) { + // If we've applied boundary conditions to either side, update the corner + if (!xlo_ylo.isEmpty() && (apply_east || apply_south)) { ParallelFor(xlo_ylo, [=] AMREX_GPU_DEVICE (int i, int j, int k) { dest_arr(i,j,k,icomp+icomp_to_fill) = 0.5 * (dest_arr(i,dom_lo.y+mf_index_type[1],k,icomp+icomp_to_fill) + dest_arr(dom_lo.x+mf_index_type[0],j,k,icomp+icomp_to_fill)); }); } - if (!xlo_yhi.isEmpty() && apply_east && apply_north) { + if (!xlo_yhi.isEmpty() && (apply_east || apply_north)) { ParallelFor(xlo_yhi, [=] AMREX_GPU_DEVICE (int i, int j, int k) { dest_arr(i,j,k,icomp+icomp_to_fill) = 0.5 * (dest_arr(i,dom_hi.y-mf_index_type[1],k,icomp+icomp_to_fill) + dest_arr(dom_lo.x+mf_index_type[0],j,k,icomp+icomp_to_fill)); }); } - if (!xhi_ylo.isEmpty() && apply_west && apply_south) { + if (!xhi_ylo.isEmpty() && (apply_west || apply_south)) { ParallelFor(xhi_ylo, [=] AMREX_GPU_DEVICE (int i, int j, int k) { dest_arr(i,j,k,icomp+icomp_to_fill) = 0.5 * (dest_arr(i,dom_lo.y+mf_index_type[1],k,icomp+icomp_to_fill) + dest_arr(dom_hi.x-mf_index_type[0],j,k,icomp+icomp_to_fill)); }); } - if (!xhi_yhi.isEmpty() && apply_west && apply_north) { + if (!xhi_yhi.isEmpty() && (apply_west || apply_north)) { ParallelFor(xhi_yhi, [=] AMREX_GPU_DEVICE (int i, int j, int k) { dest_arr(i,j,k,icomp+icomp_to_fill) = 0.5 * (dest_arr(i,dom_hi.y-mf_index_type[1],k,icomp+icomp_to_fill) + dest_arr(dom_hi.x-mf_index_type[0],j,k,icomp+icomp_to_fill)); diff --git a/Source/BoundaryConditions/BoundaryConditions_xvel.cpp b/Source/BoundaryConditions/BoundaryConditions_xvel.cpp index d4007cf..831d67a 100644 --- a/Source/BoundaryConditions/BoundaryConditions_xvel.cpp +++ b/Source/BoundaryConditions/BoundaryConditions_xvel.cpp @@ -171,22 +171,36 @@ void REMORAPhysBCFunct::impose_xvel_bcs (const Array4& dest_arr, const Box Box xlo_yhi = xlo & yhi; Box xhi_ylo = xhi & ylo; Box xhi_yhi = xhi & yhi; - ParallelFor(xlo_ylo, [=] AMREX_GPU_DEVICE (int i, int j, int k) - { - dest_arr(i,j,k) = 0.5 * (dest_arr(i,dom_lo.y,k) + dest_arr(dom_lo.x+1,j,k)); - }); - ParallelFor(xlo_yhi, [=] AMREX_GPU_DEVICE (int i, int j, int k) - { - dest_arr(i,j,k) = 0.5 * (dest_arr(i,dom_hi.y,k) + dest_arr(dom_lo.x+1,j,k)); - }); - ParallelFor(xhi_ylo, [=] AMREX_GPU_DEVICE (int i, int j, int k) - { - dest_arr(i,j,k) = 0.5 * (dest_arr(i,dom_lo.y,k) + dest_arr(dom_hi.x,j,k)); - }); - ParallelFor(xhi_yhi, [=] AMREX_GPU_DEVICE (int i, int j, int k) - { - dest_arr(i,j,k) = 0.5 * (dest_arr(i,dom_hi.y,k) + dest_arr(dom_hi.x,j,k)); - }); + + const bool clamp_east = m_domain_bcs_type[bccomp].lo(0) == REMORABCType::clamped; + const bool clamp_west = m_domain_bcs_type[bccomp].hi(0) == REMORABCType::clamped; + const bool clamp_south = m_domain_bcs_type[bccomp].lo(1) == REMORABCType::clamped; + const bool clamp_north = m_domain_bcs_type[bccomp].hi(1) == REMORABCType::clamped; + + if (!clamp_east && !clamp_south) { + ParallelFor(xlo_ylo, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + dest_arr(i,j,k) = 0.5 * (dest_arr(i,dom_lo.y,k) + dest_arr(dom_lo.x+1,j,k)); + }); + } + if (!clamp_east && !clamp_north) { + ParallelFor(xlo_yhi, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + dest_arr(i,j,k) = 0.5 * (dest_arr(i,dom_hi.y,k) + dest_arr(dom_lo.x+1,j,k)); + }); + } + if (!clamp_west && !clamp_south) { + ParallelFor(xhi_ylo, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + dest_arr(i,j,k) = 0.5 * (dest_arr(i,dom_lo.y,k) + dest_arr(dom_hi.x,j,k)); + }); + } + if (!clamp_west && !clamp_north) { + ParallelFor(xhi_yhi, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + dest_arr(i,j,k) = 0.5 * (dest_arr(i,dom_hi.y,k) + dest_arr(dom_hi.x,j,k)); + }); + } } diff --git a/Source/BoundaryConditions/BoundaryConditions_yvel.cpp b/Source/BoundaryConditions/BoundaryConditions_yvel.cpp index c89d50a..a45ff6e 100644 --- a/Source/BoundaryConditions/BoundaryConditions_yvel.cpp +++ b/Source/BoundaryConditions/BoundaryConditions_yvel.cpp @@ -173,22 +173,35 @@ void REMORAPhysBCFunct::impose_yvel_bcs (const Array4& dest_arr, const Box Box xlo_yhi = xlo & yhi; Box xhi_ylo = xhi & ylo; Box xhi_yhi = xhi & yhi; - ParallelFor(xlo_ylo, [=] AMREX_GPU_DEVICE (int i, int j, int k) - { - dest_arr(i,j,k) = 0.5 * (dest_arr(i,dom_lo.y+1,k) + dest_arr(dom_lo.x,j,k)); - }); - ParallelFor(xlo_yhi, [=] AMREX_GPU_DEVICE (int i, int j, int k) - { - dest_arr(i,j,k) = 0.5 * (dest_arr(i,dom_hi.y,k) + dest_arr(dom_lo.x,j,k)); - }); - ParallelFor(xhi_ylo, [=] AMREX_GPU_DEVICE (int i, int j, int k) - { - dest_arr(i,j,k) = 0.5 * (dest_arr(i,dom_lo.y+1,k) + dest_arr(dom_hi.x,j,k)); - }); - ParallelFor(xhi_yhi, [=] AMREX_GPU_DEVICE (int i, int j, int k) - { - dest_arr(i,j,k) = 0.5 * (dest_arr(i,dom_hi.y,k) + dest_arr(dom_hi.x,j,k)); - }); + const bool clamp_east = m_domain_bcs_type[bccomp].lo(0) == REMORABCType::clamped; + const bool clamp_west = m_domain_bcs_type[bccomp].hi(0) == REMORABCType::clamped; + const bool clamp_south = m_domain_bcs_type[bccomp].lo(1) == REMORABCType::clamped; + const bool clamp_north = m_domain_bcs_type[bccomp].hi(1) == REMORABCType::clamped; + + if (!clamp_east && !clamp_south) { + ParallelFor(xlo_ylo, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + dest_arr(i,j,k) = 0.5 * (dest_arr(i,dom_lo.y+1,k) + dest_arr(dom_lo.x,j,k)); + }); + } + if (!clamp_east && !clamp_north) { + ParallelFor(xlo_yhi, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + dest_arr(i,j,k) = 0.5 * (dest_arr(i,dom_hi.y,k) + dest_arr(dom_lo.x,j,k)); + }); + } + if (!clamp_west && !clamp_south) { + ParallelFor(xhi_ylo, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + dest_arr(i,j,k) = 0.5 * (dest_arr(i,dom_lo.y+1,k) + dest_arr(dom_hi.x,j,k)); + }); + } + if (!clamp_west && !clamp_north) { + ParallelFor(xhi_yhi, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + dest_arr(i,j,k) = 0.5 * (dest_arr(i,dom_hi.y,k) + dest_arr(dom_hi.x,j,k)); + }); + } } Gpu::streamSynchronize(); diff --git a/Source/TimeIntegration/REMORA_gls.cpp b/Source/TimeIntegration/REMORA_gls.cpp index a34d769..9bc4536 100644 --- a/Source/TimeIntegration/REMORA_gls.cpp +++ b/Source/TimeIntegration/REMORA_gls.cpp @@ -34,6 +34,10 @@ REMORA::gls_prestep (int lev, MultiFab* mf_gls, MultiFab* mf_tke, const auto dlo = amrex::lbound(domain); const auto dhi = amrex::ubound(domain); + GeometryData const& geomdata = geom[0].data(); + bool is_periodic_in_x = geomdata.isPeriodic(0); + bool is_periodic_in_y = geomdata.isPeriodic(1); + int ncomp = 1; Vector bcrs_x(ncomp); Vector bcrs_y(ncomp); @@ -77,11 +81,11 @@ REMORA::gls_prestep (int lev, MultiFab* mf_gls, MultiFab* mf_tke, // Adjust boundaries // TODO: Make sure indices match with what ROMS does - if (i == dlo.x-1 && (bcr_x.lo(0) == REMORABCType::ext_dir or ic_bc_type==IC_BC_Type::Real)) { + if (i == dlo.x-1 && !is_periodic_in_x) { grad_im1 = tke(i,j,k,nstp) - tke(i-1,j,k,nstp); gradL_im1 = gls(i,j,k,nstp) - gls(i-1,j,k,nstp); } - else if (i == dhi.x+1 && (bcr_x.hi(0) == REMORABCType::ext_dir or ic_bc_type==IC_BC_Type::Real)) { + else if (i == dhi.x+1 && !is_periodic_in_x) { grad_ip1 = tke(i,j,k,nstp) - tke(i-1,j,k,nstp); gradL_ip1 = gls(i,j,k,nstp) - gls(i-1,j,k,nstp); } @@ -104,11 +108,11 @@ REMORA::gls_prestep (int lev, MultiFab* mf_gls, MultiFab* mf_tke, // Adjust boundaries // TODO: Make sure indices match with what ROMS does - if (j == dlo.y-1 && (bcr_y.lo(1) == REMORABCType::ext_dir or ic_bc_type==IC_BC_Type::Real)) { + if (j == dlo.y-1 && !is_periodic_in_y) { grad_jm1 = tke(i,j,k,nstp) - tke(i,j-1,k,nstp); gradL_jm1 = gls(i,j,k,nstp) - gls(i,j-1,k,nstp); } - else if (j == dhi.y+1 && (bcr_y.hi(1) == REMORABCType::ext_dir or ic_bc_type==IC_BC_Type::Real)) { + else if (j == dhi.y+1 && !is_periodic_in_y) { grad_jp1 = tke(i,j,k,nstp) - tke(i,j-1,k,nstp); gradL_jp1 = gls(i,j,k,nstp) - gls(i,j-1,k,nstp); } @@ -371,6 +375,14 @@ REMORA::gls_corrector (int lev, MultiFab* mf_gls, MultiFab* mf_tke, MultiFab mf_BCK(ba,dm,1,IntVect(NGROW,NGROW,0)); MultiFab mf_BCP(ba,dm,1,IntVect(NGROW,NGROW,0)); + const Box& domain = geom[0].Domain(); + const auto dlo = amrex::lbound(domain); + const auto dhi = amrex::ubound(domain); + + GeometryData const& geomdata = geom[0].data(); + bool is_periodic_in_x = geomdata.isPeriodic(0); + bool is_periodic_in_y = geomdata.isPeriodic(1); + for ( MFIter mfi(*mf_gls, TilingIfNotGPU()); mfi.isValid(); ++mfi ) { @@ -391,10 +403,6 @@ REMORA::gls_corrector (int lev, MultiFab* mf_gls, MultiFab* mf_tke, auto CF = mf_CF.array(mfi); auto shear2_cached = mf_shear2_cached.array(mfi); - const Box& domain = geom[0].Domain(); - const auto dlo = amrex::lbound(domain); - const auto dhi = amrex::ubound(domain); - ParallelFor(gbx1D, [=] AMREX_GPU_DEVICE (int i, int j, int ) { CF(i,j,0) = 0.0_rt; @@ -532,12 +540,12 @@ REMORA::gls_corrector (int lev, MultiFab* mf_gls, MultiFab* mf_tke, { Real gradK, gradK_ip1, gradP, gradP_ip1; - if (i == dlo.x-1 && (bcr_x.lo(0) == REMORABCType::ext_dir or ic_bc_type==IC_BC_Type::Real)) { + if (i == dlo.x-1 && !is_periodic_in_x) { gradK_ip1 = tke(i+1,j,k,2)-tke(i ,j,k,2); gradK = gradK_ip1; gradP_ip1 = gls(i+1,j,k,2)-gls(i ,j,k,2); gradP = gradP_ip1; - } else if (i == dhi.x+1 && (bcr_x.hi(0) == REMORABCType::ext_dir or ic_bc_type==IC_BC_Type::Real)) { + } else if (i == dhi.x+1 && !is_periodic_in_x) { gradK = tke(i ,j,k,2)-tke(i-1,j,k,2); gradK_ip1 = gradK; gradP = gls(i ,j,k,2)-gls(i-1,j,k,2); @@ -570,11 +578,11 @@ REMORA::gls_corrector (int lev, MultiFab* mf_gls, MultiFab* mf_tke, Real gradP = (gls(i,j ,k,2)-gls(i,j-1,k,2)) * mskv(i,j ,0); Real gradP_jp1 = (gls(i,j+1,k,2)-gls(i,j ,k,2)) * mskv(i,j+1,0); - if (j == dlo.y-1 && (bcr_y.lo(1) == REMORABCType::ext_dir or ic_bc_type==IC_BC_Type::Real)) { + if (j == dlo.y-1 && !is_periodic_in_y) { gradK = gradK_jp1; gradP = gradP_jp1; } - else if (j == dhi.y+1 && (bcr_y.hi(1) == REMORABCType::ext_dir or ic_bc_type==IC_BC_Type::Real)) { + else if (j == dhi.y+1 && !is_periodic_in_y) { gradK_jp1 = gradK; gradP_jp1 = gradP; } diff --git a/Source/TimeIntegration/REMORA_rhs_uv_2d.cpp b/Source/TimeIntegration/REMORA_rhs_uv_2d.cpp index b5432c7..9c687fb 100644 --- a/Source/TimeIntegration/REMORA_rhs_uv_2d.cpp +++ b/Source/TimeIntegration/REMORA_rhs_uv_2d.cpp @@ -30,6 +30,10 @@ REMORA::rhs_uv_2d (const Box& xbx, const Box& ybx, const auto dlo = amrex::lbound(domain); const auto dhi = amrex::ubound(domain); + GeometryData const& geomdata = geom[0].data(); + bool is_periodic_in_x = geomdata.isPeriodic(0); + bool is_periodic_in_y = geomdata.isPeriodic(1); + int ncomp = 1; Vector bcrs_x(ncomp); Vector bcrs_y(ncomp); @@ -82,11 +86,11 @@ REMORA::rhs_uv_2d (const Box& xbx, const Box& ybx, Real Huxx_i = DUon(i-1,j,0)-2.0_rt*DUon(i ,j,0)+DUon(i+1,j,0); Real Huxx_ip1 = DUon(i ,j,0)-2.0_rt*DUon(i+1,j,0)+DUon(i+2,j,0); - if (i == dlo.x && (bcr_x.lo(0) == REMORABCType::ext_dir or ic_bc_type==IC_BC_Type::Real)) { + if (i == dlo.x && !is_periodic_in_x) { uxx_i = uxx_ip1; Huxx_i = Huxx_ip1; } - else if (i == dhi.x && (bcr_x.hi(0) == REMORABCType::ext_dir or ic_bc_type==IC_BC_Type::Real)) { + else if (i == dhi.x && !is_periodic_in_x) { uxx_ip1 = uxx_i; Huxx_ip1 = Huxx_i; } @@ -110,9 +114,9 @@ REMORA::rhs_uv_2d (const Box& xbx, const Box& ybx, Real Hvxx_i = DVom(i-1,j,0)-2.0_rt*DVom(i ,j,0)+DVom(i+1,j,0); Real Hvxx_im1 = DVom(i-2,j,0)-2.0_rt*DVom(i-1,j,0)+DVom(i ,j,0); - if (j == dlo.y and (bcr_y.lo(1) == REMORABCType::ext_dir or ic_bc_type==IC_BC_Type::Real)) { + if (j == dlo.y and !is_periodic_in_y) { uee_jm1 = uee_j; - } else if (j == dhi.y+1 and (bcr_y.hi(1) == REMORABCType::ext_dir or ic_bc_type==IC_BC_Type::Real)) { + } else if (j == dhi.y+1 and !is_periodic_in_y) { uee_j = uee_jm1; } @@ -174,9 +178,9 @@ REMORA::rhs_uv_2d (const Box& xbx, const Box& ybx, Real Huee_j = DUon(i,j-1,0)-2.0_rt*DUon(i,j ,0)+DUon(i,j+1,0); Real Huee_jm1 = DUon(i,j-2,0)-2.0_rt*DUon(i,j-1,0)+DUon(i,j ,0); - if (i == dlo.x and (bcr_x.lo(0) == REMORABCType::ext_dir or ic_bc_type==IC_BC_Type::Real)) { + if (i == dlo.x and !is_periodic_in_x) { vxx_im1 = vxx_i; - } else if (i == dhi.x + 1 and (bcr_x.hi(0) == REMORABCType::ext_dir or ic_bc_type==IC_BC_Type::Real)) { + } else if (i == dhi.x + 1 and !is_periodic_in_x) { vxx_i = vxx_im1; } @@ -195,11 +199,11 @@ REMORA::rhs_uv_2d (const Box& xbx, const Box& ybx, Real Hvee_j = DVom(i,j-1,0)-2.0_rt*DVom(i,j ,0)+DVom(i,j+1,0); Real Hvee_jp1 = DVom(i,j ,0)-2.0_rt*DVom(i,j+1,0)+DVom(i,j+2,0); - if (j == dlo.y and (bcr_y.lo(1) == REMORABCType::ext_dir or ic_bc_type==IC_BC_Type::Real)) { + if (j == dlo.y and !is_periodic_in_y) { vee_j = vee_jp1; Hvee_j = Hvee_jp1; } - else if (j == dhi.y and (bcr_y.hi(1) == REMORABCType::ext_dir or ic_bc_type==IC_BC_Type::Real)) { + else if (j == dhi.y and !is_periodic_in_y) { vee_jp1 = vee_j; Hvee_jp1 = Hvee_j; } diff --git a/Source/TimeIntegration/REMORA_rhs_uv_3d.cpp b/Source/TimeIntegration/REMORA_rhs_uv_3d.cpp index 0f99ad9..223980e 100644 --- a/Source/TimeIntegration/REMORA_rhs_uv_3d.cpp +++ b/Source/TimeIntegration/REMORA_rhs_uv_3d.cpp @@ -51,6 +51,10 @@ REMORA::rhs_uv_3d (const Box& xbx, const Box& ybx, const auto dlo = amrex::lbound(domain); const auto dhi = amrex::ubound(domain); + GeometryData const& geomdata = geom[0].data(); + bool is_periodic_in_x = geomdata.isPeriodic(0); + bool is_periodic_in_y = geomdata.isPeriodic(1); + int ncomp = 1; Vector bcrs_x(ncomp); Vector bcrs_y(ncomp); @@ -108,11 +112,11 @@ REMORA::rhs_uv_3d (const Box& xbx, const Box& ybx, Real Huxx_i = Huon(i-1,j,k)-2.0_rt*Huon(i ,j,k)+Huon(i+1,j,k); Real Huxx_ip1 = Huon(i ,j,k)-2.0_rt*Huon(i+1,j,k)+Huon(i+2,j,k); - if (i == dlo.x && (bcr_x.lo(0) == REMORABCType::ext_dir or ic_bc_type==IC_BC_Type::Real)) { + if (i == dlo.x && !is_periodic_in_x) { uxx_i = uxx_ip1; Huxx_i = Huxx_ip1; } - else if (i == dhi.x && (bcr_x.hi(0) == REMORABCType::ext_dir or ic_bc_type==IC_BC_Type::Real)) { + else if (i == dhi.x && !is_periodic_in_x) { uxx_ip1 = uxx_i; Huxx_ip1 = Huxx_i; } @@ -135,9 +139,9 @@ REMORA::rhs_uv_3d (const Box& xbx, const Box& ybx, Real uee_jm1 = uold(i,j-2,k,nrhs) - 2.0_rt*uold(i,j-1,k,nrhs) + uold(i ,j,k,nrhs); Real uee_j = uold(i,j-1,k,nrhs) - 2.0_rt*uold(i,j ,k,nrhs) + uold(i,j+1,k,nrhs); - if (j == dlo.y and (bcr_y.lo(1) == REMORABCType::ext_dir or ic_bc_type==IC_BC_Type::Real)) { + if (j == dlo.y and !is_periodic_in_y) { uee_jm1 = uee_j; - } else if (j == dhi.y+1 and (bcr_y.hi(1) == REMORABCType::ext_dir or ic_bc_type==IC_BC_Type::Real)) { + } else if (j == dhi.y+1 and !is_periodic_in_y) { uee_j = uee_jm1; } @@ -263,9 +267,9 @@ REMORA::rhs_uv_3d (const Box& xbx, const Box& ybx, Real vxx_im1 = vold(i-2,j,k,nrhs)-2.0_rt*vold(i-1,j,k,nrhs)+vold(i ,j,k,nrhs); Real vxx_i = vold(i-1,j,k,nrhs)-2.0_rt*vold(i ,j,k,nrhs)+vold(i+1,j,k,nrhs); - if (i == dlo.x and (bcr_x.lo(0) == REMORABCType::ext_dir or ic_bc_type==IC_BC_Type::Real)) { + if (i == dlo.x and !is_periodic_in_x) { vxx_im1 = vxx_i; - } else if (i == dhi.x+1 and (bcr_x.hi(0) == REMORABCType::ext_dir or ic_bc_type==IC_BC_Type::Real)) { + } else if (i == dhi.x+1 and !is_periodic_in_x) { vxx_i = vxx_im1; } @@ -291,11 +295,11 @@ REMORA::rhs_uv_3d (const Box& xbx, const Box& ybx, Real Hvee_j = Hvom(i,j-1,k)-2.0_rt*Hvom(i,j ,k)+Hvom(i,j+1,k); Real Hvee_jp1 = Hvom(i,j ,k)-2.0_rt*Hvom(i,j+1,k)+Hvom(i,j+2,k); - if (j == dlo.y and (bcr_y.lo(1) == REMORABCType::ext_dir or ic_bc_type==IC_BC_Type::Real)) { + if (j == dlo.y and !is_periodic_in_y) { vee_j = vee_jp1; Hvee_j = Hvee_jp1; } - else if (j == dhi.y and (bcr_y.hi(1) == REMORABCType::ext_dir or ic_bc_type==IC_BC_Type::Real)) { + else if (j == dhi.y and !is_periodic_in_y) { vee_jp1 = vee_j; Hvee_jp1 = Hvee_j; } diff --git a/Source/TimeIntegration/REMORA_update_massflux_3d.cpp b/Source/TimeIntegration/REMORA_update_massflux_3d.cpp index 0b7dfe8..98757e8 100644 --- a/Source/TimeIntegration/REMORA_update_massflux_3d.cpp +++ b/Source/TimeIntegration/REMORA_update_massflux_3d.cpp @@ -39,6 +39,10 @@ REMORA::update_massflux_3d (const Box& bx, const auto dlo = amrex::lbound(domain); const auto dhi = amrex::ubound(domain); + GeometryData const& geomdata = geom[0].data(); + bool is_periodic_in_x = geomdata.isPeriodic(0); + bool is_periodic_in_y = geomdata.isPeriodic(1); + auto ic_bc_type = solverChoice.ic_bc_type; int ncomp = 1; @@ -59,10 +63,6 @@ REMORA::update_massflux_3d (const Box& bx, FArrayBox fab_CF(bxD,1,amrex::The_Async_Arena()); fab_CF.template setVal(0.); auto CF=fab_CF.array(); - // auto geomdata = Geom(0).data(); - // bool NSPeriodic = geomdata.isPeriodic(1); - // bool EWPeriodic = geomdata.isPeriodic(0); - //Copied depth of water column calculation from DepthStretchTransform //Compute thicknesses of U-boxes DC(i,j,0:N-1), total depth of the water column DC(i,j,-1) @@ -84,19 +84,22 @@ REMORA::update_massflux_3d (const Box& bx, DC(i,j,-1) = 1.0_rt / DC(i,j,-1); CF(i,j,0) = DC(i,j,-1) * (CF(i,j,0) - Dphi_avg1(i,j,0)); + // In order to agree with ROMS on the boundaries, the corner points shouldn't actually + // be updated with CF for clamped E/W, wall N/S boundaries. This doesn't seem to affect + // the interior valid points, though for (int k=0; k<=N; k++) { - if (i == dlo.x-joff && (bcr_x.lo(0) == REMORABCType::ext_dir or ic_bc_type==IC_BC_Type::Real)) { + if (i == dlo.x-joff && !is_periodic_in_x) { phi(i,j,k,nnew) -= CF(i,j,0); phi(i,j,k,nnew) *= msk(i,j,0); - } else if (i == dhi.x+1 && (bcr_x.hi(0) == REMORABCType::ext_dir or ic_bc_type==IC_BC_Type::Real)) { + } else if (i == dhi.x+1 && !is_periodic_in_x) { phi(i,j,k,nnew) -= CF(i,j,0); phi(i,j,k,nnew) *= msk(i,j,0); } - if (j == dlo.y-ioff && (bcr_y.lo(1) == REMORABCType::ext_dir or ic_bc_type==IC_BC_Type::Real)) { + if (j == dlo.y-ioff && !is_periodic_in_y) { phi(i,j,k,nnew) -= CF(i,j,0); phi(i,j,k,nnew) *= msk(i,j,0); - } else if (j == dhi.y+1 && (bcr_y.hi(1) == REMORABCType::ext_dir or ic_bc_type==IC_BC_Type::Real)) { + } else if (j == dhi.y+1 && !is_periodic_in_y) { phi(i,j,k,nnew) -= CF(i,j,0); phi(i,j,k,nnew) *= msk(i,j,0); }