Skip to content

Commit

Permalink
Fix relaxation zone code. (#1148)
Browse files Browse the repository at this point in the history 
Co-authored-by: Ann Almgren <[email protected]>
  • Loading branch information
@AMLattanzi @asalmgren
AMLattanzi and asalmgren authored Jul 6, 2023
1 parent 3c1b866 commit e03ef78
Show file tree
Hide file tree
Showing 3 changed files with 77 additions and 56 deletions.
50 changes: 25 additions & 25 deletions Source/ERF.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1525,33 +1525,33 @@ ERF::AverageDownTo (int crse_lev) // NOLINT
void
ERF::define_grids_to_evolve (int lev) // NOLINT
{
int width = wrfbdy_width - 1;
Box domain(geom[lev].Domain());
if (lev == 0 && ( init_type == "real" || init_type == "metgrid" ) )
{
Box shrunk_domain(domain);
shrunk_domain.grow(0,-width);
shrunk_domain.grow(1,-width);
grids_to_evolve[lev] = amrex::intersect(grids[lev],shrunk_domain);
} else if (lev == 1) {
Box shrunk_domain(boxes_at_level[lev][0]);
shrunk_domain.grow(0,-width);
shrunk_domain.grow(1,-width);
grids_to_evolve[lev] = amrex::intersect(grids[lev],shrunk_domain);
int width = wrfbdy_set_width;
Box domain(geom[lev].Domain());
if (lev == 0 && ( init_type == "real" || init_type == "metgrid" ) )
{
Box shrunk_domain(domain);
shrunk_domain.grow(0,-width);
shrunk_domain.grow(1,-width);
grids_to_evolve[lev] = amrex::intersect(grids[lev],shrunk_domain);
} else if (lev == 1) {
Box shrunk_domain(boxes_at_level[lev][0]);
shrunk_domain.grow(0,-width);
shrunk_domain.grow(1,-width);
grids_to_evolve[lev] = amrex::intersect(grids[lev],shrunk_domain);
#if 0
if (num_boxes_at_level[lev] > 1) {
for (int i = 1; i < num_boxes_at_level[lev]; i++) {
Box shrunk_domain(boxes_at_level[lev][i]);
shrunk_domain.grow(0,-width);
shrunk_domain.grow(1,-width);
grids_to_evolve[lev] = amrex::intersect(grids_to_evolve[lev],shrunk_domain);
}
}
if (num_boxes_at_level[lev] > 1) {
for (int i = 1; i < num_boxes_at_level[lev]; i++) {
Box shrunk_domain(boxes_at_level[lev][i]);
shrunk_domain.grow(0,-width);
shrunk_domain.grow(1,-width);
grids_to_evolve[lev] = amrex::intersect(grids_to_evolve[lev],shrunk_domain);
}
}
#endif
} else {
// Just copy grids...
grids_to_evolve[lev] = grids[lev];
}
} else {
// Just copy grids...
grids_to_evolve[lev] = grids[lev];
}
}

#ifdef ERF_USE_MULTIBLOCK
Expand Down
2 changes: 1 addition & 1 deletion Source/TimeIntegration/TI_slow_rhs_fun.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -194,7 +194,7 @@

#ifdef ERF_USE_NETCDF
// Populate RHS for relaxation zones
if (init_type == "real" && level == 0)
if (init_type=="real" && level==0 && nrk==0)
compute_interior_ghost_RHS(bdy_time_interval, old_stage_time, slow_dt,
wrfbdy_width-1, wrfbdy_set_width, fine_geom,
S_rhs, S_data, bdy_data_xlo, bdy_data_xhi,
Expand Down
81 changes: 51 additions & 30 deletions Source/Utils/InteriorGhostCells.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -420,7 +420,7 @@ compute_interior_ghost_RHS(const Real& bdy_time_interval,
// 4 halo boxes w/ no ghost cells for CC vars
Box domain = geom.Domain();
Box bx_xlo, bx_xhi, bx_ylo, bx_yhi;
compute_interior_ghost_bxs_xy(tbx, domain, width, 0,
compute_interior_ghost_bxs_xy(tbx, domain, set_width, 0,
bx_xlo, bx_xhi,
bx_ylo, bx_yhi);
S_rhs[IntVar::cons][mfi].template setVal<RunOn::Device>(0.,bx_xlo,Rho_comp,NVAR);
Expand All @@ -437,7 +437,7 @@ compute_interior_ghost_RHS(const Real& bdy_time_interval,
Box domain = geom.Domain();
domain.convert(S_rhs[IntVar::zmom].boxArray().ixType());
Box bx_xlo, bx_xhi, bx_ylo, bx_yhi;
compute_interior_ghost_bxs_xy(tbx, domain, width, 0,
compute_interior_ghost_bxs_xy(tbx, domain, set_width, 0,
bx_xlo, bx_xhi,
bx_ylo, bx_yhi);
S_rhs[IntVar::zmom][mfi].template setVal<RunOn::Device>(0.,bx_xlo);
Expand All @@ -457,6 +457,7 @@ compute_interior_ghost_RHS(const Real& bdy_time_interval,
Box domain = geom.Domain();
domain.convert(S_data[var_idx].boxArray().ixType());
const auto& dom_hi = ubound(domain);
const auto& dom_lo = lbound(domain);

#ifdef _OPENMP
#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
Expand Down Expand Up @@ -517,60 +518,75 @@ compute_interior_ghost_RHS(const Real& bdy_time_interval,
[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
// Corners with x boxes
int j_lo = std::min(j,width-1);
int j_lo = std::min(j-dom_lo.y,width-1);
int j_hi = std::min(dom_hi.y-j,width-1);
int jj = std::min(j_lo,j_hi);
int n = std::min(i,jj) + 1;
int n = std::min(i-dom_lo.x,jj) + 1;
if (n <= Spec_z) {
rhs_arr(i,j,k,icomp) = 0.0;
} else {
Real Factor = (num - Real(n))/denom;
Real delta = arr_xlo(i ,j ,k) - data_arr(i ,j ,k,icomp);
Real delta_xp = arr_xlo(i+1,j ,k) - data_arr(i+1,j ,k,icomp);
Real delta_xm = arr_xlo(i-1,j ,k) - data_arr(i-1,j ,k,icomp);
Real delta_yp = arr_xlo(i ,j+1,k) - data_arr(i ,j+1,k,icomp);
Real delta_ym = arr_xlo(i ,j-1,k) - data_arr(i ,j-1,k,icomp);
Real d = data_arr(i ,j ,k,icomp) + delta_t*rhs_arr(i , j , k,icomp);
Real d_ip1 = data_arr(i+1,j ,k,icomp) + delta_t*rhs_arr(i+1, j , k,icomp);
Real d_im1 = data_arr(i-1,j ,k,icomp) + delta_t*rhs_arr(i-1, j , k,icomp);
Real d_jp1 = data_arr(i ,j+1,k,icomp) + delta_t*rhs_arr(i , j+1, k,icomp);
Real d_jm1 = data_arr(i ,j-1,k,icomp) + delta_t*rhs_arr(i , j-1, k,icomp);
Real delta = arr_xlo(i ,j ,k) - d;
Real delta_xp = arr_xlo(i+1,j ,k) - d_ip1;
Real delta_xm = arr_xlo(i-1,j ,k) - d_im1;
Real delta_yp = arr_xlo(i ,j+1,k) - d_jp1;
Real delta_ym = arr_xlo(i ,j-1,k) - d_jm1;
Real Laplacian = delta_xp + delta_xm + delta_yp + delta_ym - 4.0*delta;
rhs_arr(i,j,k,icomp) = (F1*delta - F2*Laplacian) * Factor;
rhs_arr(i,j,k,icomp) += (F1*delta - F2*Laplacian) * Factor;
}
},
[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
// Corners with x boxes
int j_lo = std::min(j,width-1);
int j_lo = std::min(j-dom_lo.y,width-1);
int j_hi = std::min(dom_hi.y-j,width-1);
int jj = std::min(j_lo,j_hi);
int n = std::min(dom_hi.x-i,jj) + 1;
if (n <= Spec_z) {
rhs_arr(i,j,k,icomp) = 0.0;
} else {
Real Factor = (num - Real(n))/denom;
Real delta = arr_xhi(i ,j ,k) - data_arr(i ,j ,k,icomp);
Real delta_xp = arr_xhi(i+1,j ,k) - data_arr(i+1,j ,k,icomp);
Real delta_xm = arr_xhi(i-1,j ,k) - data_arr(i-1,j ,k,icomp);
Real delta_yp = arr_xhi(i ,j+1,k) - data_arr(i ,j+1,k,icomp);
Real delta_ym = arr_xhi(i ,j-1,k) - data_arr(i ,j-1,k,icomp);
Real d = data_arr(i ,j ,k,icomp) + delta_t*rhs_arr(i , j , k,icomp);
Real d_ip1 = data_arr(i+1,j ,k,icomp) + delta_t*rhs_arr(i+1, j , k,icomp);
Real d_im1 = data_arr(i-1,j ,k,icomp) + delta_t*rhs_arr(i-1, j , k,icomp);
Real d_jp1 = data_arr(i ,j+1,k,icomp) + delta_t*rhs_arr(i , j+1, k,icomp);
Real d_jm1 = data_arr(i ,j-1,k,icomp) + delta_t*rhs_arr(i , j-1, k,icomp);
Real delta = arr_xlo(i ,j ,k) - d;
Real delta_xp = arr_xlo(i+1,j ,k) - d_ip1;
Real delta_xm = arr_xlo(i-1,j ,k) - d_im1;
Real delta_yp = arr_xlo(i ,j+1,k) - d_jp1;
Real delta_ym = arr_xlo(i ,j-1,k) - d_jm1;
Real Laplacian = delta_xp + delta_xm + delta_yp + delta_ym - 4.0*delta;
rhs_arr(i,j,k,icomp) = (F1*delta - F2*Laplacian) * Factor;
rhs_arr(i,j,k,icomp) += (F1*delta - F2*Laplacian) * Factor;
}
});

amrex::ParallelFor(tbx_ylo, tbx_yhi,
[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
// No corners for y boxes
int n = j + 1;
int n = j -dom_lo.y + 1;
if (n <= Spec_z) {
rhs_arr(i,j,k,icomp) = 0.0;
} else {
Real Factor = (num - Real(n))/denom;
Real delta = arr_ylo(i ,j ,k) - data_arr(i ,j ,k,icomp);
Real delta_xp = arr_ylo(i+1,j ,k) - data_arr(i+1,j ,k,icomp);
Real delta_xm = arr_ylo(i-1,j ,k) - data_arr(i-1,j ,k,icomp);
Real delta_yp = arr_ylo(i ,j+1,k) - data_arr(i ,j+1,k,icomp);
Real delta_ym = arr_ylo(i ,j-1,k) - data_arr(i ,j-1,k,icomp);
Real d = data_arr(i ,j ,k,icomp) + delta_t*rhs_arr(i , j , k,icomp);
Real d_ip1 = data_arr(i+1,j ,k,icomp) + delta_t*rhs_arr(i+1, j , k,icomp);
Real d_im1 = data_arr(i-1,j ,k,icomp) + delta_t*rhs_arr(i-1, j , k,icomp);
Real d_jp1 = data_arr(i ,j+1,k,icomp) + delta_t*rhs_arr(i , j+1, k,icomp);
Real d_jm1 = data_arr(i ,j-1,k,icomp) + delta_t*rhs_arr(i , j-1, k,icomp);
Real delta = arr_xlo(i ,j ,k) - d;
Real delta_xp = arr_xlo(i+1,j ,k) - d_ip1;
Real delta_xm = arr_xlo(i-1,j ,k) - d_im1;
Real delta_yp = arr_xlo(i ,j+1,k) - d_jp1;
Real delta_ym = arr_xlo(i ,j-1,k) - d_jm1;
Real Laplacian = delta_xp + delta_xm + delta_yp + delta_ym - 4.0*delta;
rhs_arr(i,j,k,icomp) = (F1*delta - F2*Laplacian) * Factor;
rhs_arr(i,j,k,icomp) += (F1*delta - F2*Laplacian) * Factor;
}
},
[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
Expand All @@ -581,13 +597,18 @@ compute_interior_ghost_RHS(const Real& bdy_time_interval,
rhs_arr(i,j,k,icomp) = 0.0;
} else {
Real Factor = (num - Real(n))/denom;
Real delta = arr_yhi(i ,j ,k) - data_arr(i ,j ,k,icomp);
Real delta_xp = arr_yhi(i+1,j ,k) - data_arr(i+1,j ,k,icomp);
Real delta_xm = arr_yhi(i-1,j ,k) - data_arr(i-1,j ,k,icomp);
Real delta_yp = arr_yhi(i ,j+1,k) - data_arr(i ,j+1,k,icomp);
Real delta_ym = arr_yhi(i ,j-1,k) - data_arr(i ,j-1,k,icomp);
Real d = data_arr(i ,j ,k,icomp) + delta_t*rhs_arr(i , j , k,icomp);
Real d_ip1 = data_arr(i+1,j ,k,icomp) + delta_t*rhs_arr(i+1, j , k,icomp);
Real d_im1 = data_arr(i-1,j ,k,icomp) + delta_t*rhs_arr(i-1, j , k,icomp);
Real d_jp1 = data_arr(i ,j+1,k,icomp) + delta_t*rhs_arr(i , j+1, k,icomp);
Real d_jm1 = data_arr(i ,j-1,k,icomp) + delta_t*rhs_arr(i , j-1, k,icomp);
Real delta = arr_xlo(i ,j ,k) - d;
Real delta_xp = arr_xlo(i+1,j ,k) - d_ip1;
Real delta_xm = arr_xlo(i-1,j ,k) - d_im1;
Real delta_yp = arr_xlo(i ,j+1,k) - d_jp1;
Real delta_ym = arr_xlo(i ,j-1,k) - d_jm1;
Real Laplacian = delta_xp + delta_xm + delta_yp + delta_ym - 4.0*delta;
rhs_arr(i,j,k,icomp) = (F1*delta - F2*Laplacian) * Factor;
rhs_arr(i,j,k,icomp) += (F1*delta - F2*Laplacian) * Factor;
}
});
} // mfi
Expand Down

0 comments on commit e03ef78

Please sign in to comment.