Skip to content

Commit

Permalink
Correct z pos locations for MOST average and include an assert to ens…
Browse files Browse the repository at this point in the history 
…ure the correct k index is found for interpolation. (erf-model#1465)
  • Loading branch information
@AMLattanzi
AMLattanzi authored Feb 28, 2024
1 parent 4efa42a commit 5da3c98
Show file tree
Hide file tree
Showing 7 changed files with 24 additions and 30 deletions.
8 changes: 1 addition & 7 deletions Exec/RegTests/Terrain3d_Hemisphere/inputs_most_test
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,13 +12,6 @@ geometry.prob_hi = 10. 1. 2.
amr.n_cell = 256 8 64 # dx=dy=dz=100 m, Straka et al 1993 / Xue et al 2000

geometry.is_periodic = 1 1 0

#xlo.type = "Inflow"
#xhi.type = "Outflow"
#xlo.velocity = 1. 0. 0.
#xlo.density = 1.16
#xlo.theta = 300.
#xlo.scalar = 0.

zhi.type = "SlipWall"

Expand All @@ -40,6 +33,7 @@ erf.most.zref = 0.075 # QUERY DISTANCE (HEIGHT OR NORM LENGTH)
erf.most.radius = 1 # SPECIFIED REGION RADIUS
#-----------------------------------------------------------------
#erf.most.time_window = 50.0 # WINDOW FOR TIME AVG
#=================================================================


# TIME STEP CONTROL
Expand Down
10 changes: 7 additions & 3 deletions Source/BoundaryConditions/MOSTAverage.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -111,8 +111,9 @@ public:
const int interp_comp)
{
// Search to get z/k
bool found = false;
amrex::Real zval= 0.0;
int kmax = ubound(z_arr).z;
int kmax = ubound(z_arr).z;
int i_new = (int) (xp * dxi[0] - 0.5);
int j_new = (int) (yp * dxi[1] - 0.5);
amrex::Real z_target = zp;
Expand All @@ -122,11 +123,14 @@ public:
amrex::Real z_hi = 0.25 * ( z_arr(i_new,j_new ,lk+1) + z_arr(i_new+1,j_new ,lk+1)
+ z_arr(i_new,j_new+1,lk+1) + z_arr(i_new+1,j_new+1,lk+1) );
if (z_target > z_lo && z_target < z_hi){
found = true;
zval = (amrex::Real) lk + ((z_target - z_lo) / (z_hi - z_lo)) + 0.5;
break;
}
}

AMREX_ASSERT_WITH_MESSAGE(found, "MOSTAverage: Height above terrain not found, try increasing z_ref!");

const amrex::RealVect lx((xp - plo[0])*dxi[0] + 0.5,
(yp - plo[1])*dxi[1] + 0.5,
zval);
Expand All @@ -136,8 +140,8 @@ public:
int i = ijk[0]; int j = ijk[1]; int k = ijk[2];

// Weights
const amrex::RealVect sx_hi = lx - ijk;
const amrex::RealVect sx_lo = 1 - sx_hi;
const amrex::RealVect sx_hi = lx - ijk;
const amrex::RealVect sx_lo = 1.0 - sx_hi;

for (int n = 0; n < interp_comp; n++)
interp_vals[n] = sx_lo[0]*sx_lo[1]*sx_lo[2]*interp_array(i-1, j-1, k-1,n) +
Expand Down
26 changes: 14 additions & 12 deletions Source/BoundaryConditions/MOSTAverage.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -361,7 +361,7 @@ MOSTAverage::set_norm_indices_T()
const auto dxInv = m_geom[lev].InvCellSizeArray();
IntVect ng = m_k_indx[lev]->nGrowVect(); ng[2]=0;
for (MFIter mfi(*m_k_indx[lev], TileNoZ()); mfi.isValid(); ++mfi) {
Box npbx = mfi.tilebox(); npbx.convert({1,1,0});
Box npbx = mfi.tilebox(IntVect(1,1,0),IntVect(1,1,0));
Box gpbx = mfi.growntilebox(ng);
const auto z_phys_arr = m_z_phys_nd[lev]->const_array(mfi);
auto i_arr = m_i_indx[lev]->array(mfi);
Expand Down Expand Up @@ -429,20 +429,22 @@ MOSTAverage::set_z_positions_T()
const auto dx = m_geom[lev].CellSizeArray();
IntVect ng = m_x_pos[lev]->nGrowVect(); ng[2]=0;
for (MFIter mfi(*m_x_pos[lev], TileNoZ()); mfi.isValid(); ++mfi) {
Box npbx = mfi.tilebox(); npbx.convert({1,1,0});
Box npbx = mfi.tilebox(IntVect(1,1,0),IntVect(1,1,0));
Box gpbx = mfi.growntilebox(ng);
RealBox grb{gpbx,dx.data(),base.dataPtr()};

const auto z_phys_arr = m_z_phys_nd[lev]->const_array(mfi);
auto x_pos_arr = m_x_pos[lev]->array(mfi);
auto y_pos_arr = m_y_pos[lev]->array(mfi);
auto z_pos_arr = m_z_pos[lev]->array(mfi);
auto x_pos_arr = m_x_pos[lev]->array(mfi);
auto y_pos_arr = m_y_pos[lev]->array(mfi);
auto z_pos_arr = m_z_pos[lev]->array(mfi);
ParallelFor(npbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
// Final position at end of vector
x_pos_arr(i,j,k) = ((Real) i + 0.5) * dx[0];
y_pos_arr(i,j,k) = ((Real) j + 0.5) * dx[1];
z_pos_arr(i,j,k) = z_phys_arr(i,j,k) + d_zref;
Real z_bot_face = 0.25 * ( z_phys_arr(i ,j ,k) + z_phys_arr(i+1,j ,k)
+ z_phys_arr(i ,j+1,k) + z_phys_arr(i+1,j+1,k) );
z_pos_arr(i,j,k) = z_bot_face + d_zref;

// Destination position must be contained on the current process!
Real pos[] = {x_pos_arr(i,j,k),y_pos_arr(i,j,k),0.5*dx[2]};
Expand Down Expand Up @@ -789,9 +791,9 @@ MOSTAverage::compute_region_averages(int lev)
for (int lj(-d_radius); lj <= (d_radius); ++lj) {
for (int li(-d_radius); li <= (d_radius); ++li) {
Real interp{0};
Real xp = x_pos_arr(i,j,k) + li*dx[0];
Real yp = y_pos_arr(i,j,k) + lj*dx[1];
Real zp = z_pos_arr(i,j,k) + met_h_zeta*lk*dx[2];
Real xp = x_pos_arr(i+li,j+lj,k);
Real yp = y_pos_arr(i+li,j+lj,k);
Real zp = z_pos_arr(i+li,j+lj,k) + met_h_zeta*lk*dx[2];
trilinear_interp_T(xp, yp, zp, &interp, mf_arr, z_phys_arr, plo, dxInv, 1);
Real val = denom * interp * d_fact_new;
ma_arr(i,j,k) += val;
Expand Down Expand Up @@ -861,9 +863,9 @@ MOSTAverage::compute_region_averages(int lev)
for (int li(-d_radius); li <= (d_radius); ++li) {
Real u_interp{0};
Real v_interp{0};
Real xp = x_pos_arr(i,j,k) + li*dx[0];
Real yp = y_pos_arr(i,j,k) + lj*dx[1];
Real zp = z_pos_arr(i,j,k) + met_h_zeta*lk*dx[2];
Real xp = x_pos_arr(i+li,j+lj,k);
Real yp = y_pos_arr(i+li,j+lj,k);
Real zp = z_pos_arr(i+li,j+lj,k) + met_h_zeta*lk*dx[2];
trilinear_interp_T(xp, yp, zp, &u_interp, u_mf_arr, z_phys_arr, plo, dxInv, 1);
trilinear_interp_T(xp, yp, zp, &v_interp, v_mf_arr, z_phys_arr, plo, dxInv, 1);
Real mag = std::sqrt(u_interp*u_interp + v_interp*v_interp);
Expand Down
5 changes: 1 addition & 4 deletions Source/Diffusion/ComputeTurbulentViscosity.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,7 +47,6 @@ void ComputeTurbulentViscosityLES (const amrex::MultiFab& Tau11, const amrex::Mu
const amrex::Geometry& geom,
const amrex::MultiFab& mapfac_u, const amrex::MultiFab& mapfac_v,
const std::unique_ptr<amrex::MultiFab>& z_phys_nd,
const int& qstate_size,
const TurbChoice& turbChoice, const Real const_grav, std::unique_ptr<ABLMost>& most)
{
const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> cellSizeInv = geom.InvCellSizeArray();
Expand Down Expand Up @@ -407,7 +406,6 @@ void ComputeTurbulentViscosity (const amrex::MultiFab& xvel , const amrex::Multi
const TurbChoice& turbChoice, const Real const_grav,
std::unique_ptr<ABLMost>& most,
const amrex::BCRec* bc_ptr,
const int& qstate_size,
bool vert_only)
{
BL_PROFILE_VAR("ComputeTurbulentViscosity()",ComputeTurbulentViscosity);
Expand Down Expand Up @@ -438,8 +436,7 @@ void ComputeTurbulentViscosity (const amrex::MultiFab& xvel , const amrex::Multi
cons_in, eddyViscosity,
Hfx1, Hfx2, Hfx3, Diss,
geom, mapfac_u, mapfac_v,
z_phys_nd, qstate_size,
turbChoice, const_grav, most);
z_phys_nd, turbChoice, const_grav, most);
}

if (turbChoice.pbl_type != PBLType::None) {
Expand Down
1 change: 0 additions & 1 deletion Source/Diffusion/EddyViscosity.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,6 @@ ComputeTurbulentViscosity (const amrex::MultiFab& xvel , const amrex::MultiFab&
const TurbChoice& turbChoice, const amrex::Real const_grav,
std::unique_ptr<ABLMost>& most,
const amrex::BCRec* bc_ptr,
const int& qstate_size,
bool vert_only = false);

AMREX_GPU_DEVICE
Expand Down
2 changes: 0 additions & 2 deletions Source/Initialization/ERF_init_custom.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,8 +27,6 @@ ERF::init_custom (int lev)
{
auto& lev_new = vars_new[lev];

int n_qstate = micro.Get_Qstate_Size();

MultiFab r_hse(base_state[lev], make_alias, 0, 1); // r_0 is first component
MultiFab p_hse(base_state[lev], make_alias, 1, 1); // p_0 is second component

Expand Down
2 changes: 1 addition & 1 deletion Source/TimeIntegration/ERF_advance_dycore.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -188,7 +188,7 @@ void ERF::advance_dycore(int level,
*eddyDiffs, *Hfx1, *Hfx2, *Hfx3, *Diss, // to be updated
fine_geom, *mapfac_u[level], *mapfac_v[level],
z_phys_nd[level], tc, solverChoice.gravity,
m_most, bc_ptr_d, micro.Get_Qstate_Size());
m_most, bc_ptr_d);
}

// ***********************************************************************************************
Expand Down

0 comments on commit 5da3c98

Please sign in to comment.