diff --git a/Exec/RegTests/Terrain3d_Hemisphere/inputs_most_test b/Exec/RegTests/Terrain3d_Hemisphere/inputs_most_test index 69c2d147d..063fbf980 100644 --- a/Exec/RegTests/Terrain3d_Hemisphere/inputs_most_test +++ b/Exec/RegTests/Terrain3d_Hemisphere/inputs_most_test @@ -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" @@ -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 diff --git a/Source/BoundaryConditions/MOSTAverage.H b/Source/BoundaryConditions/MOSTAverage.H index 0d4a67842..cf55c7b0e 100644 --- a/Source/BoundaryConditions/MOSTAverage.H +++ b/Source/BoundaryConditions/MOSTAverage.H @@ -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; @@ -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); @@ -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) + diff --git a/Source/BoundaryConditions/MOSTAverage.cpp b/Source/BoundaryConditions/MOSTAverage.cpp index 73c474d04..6e00c6a89 100644 --- a/Source/BoundaryConditions/MOSTAverage.cpp +++ b/Source/BoundaryConditions/MOSTAverage.cpp @@ -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); @@ -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]}; @@ -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; @@ -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); diff --git a/Source/Diffusion/ComputeTurbulentViscosity.cpp b/Source/Diffusion/ComputeTurbulentViscosity.cpp index 7c26e1ecd..6e281e25e 100644 --- a/Source/Diffusion/ComputeTurbulentViscosity.cpp +++ b/Source/Diffusion/ComputeTurbulentViscosity.cpp @@ -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& z_phys_nd, - const int& qstate_size, const TurbChoice& turbChoice, const Real const_grav, std::unique_ptr& most) { const amrex::GpuArray cellSizeInv = geom.InvCellSizeArray(); @@ -407,7 +406,6 @@ void ComputeTurbulentViscosity (const amrex::MultiFab& xvel , const amrex::Multi const TurbChoice& turbChoice, const Real const_grav, std::unique_ptr& most, const amrex::BCRec* bc_ptr, - const int& qstate_size, bool vert_only) { BL_PROFILE_VAR("ComputeTurbulentViscosity()",ComputeTurbulentViscosity); @@ -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) { diff --git a/Source/Diffusion/EddyViscosity.H b/Source/Diffusion/EddyViscosity.H index c8fca0993..116535a15 100644 --- a/Source/Diffusion/EddyViscosity.H +++ b/Source/Diffusion/EddyViscosity.H @@ -20,7 +20,6 @@ ComputeTurbulentViscosity (const amrex::MultiFab& xvel , const amrex::MultiFab& const TurbChoice& turbChoice, const amrex::Real const_grav, std::unique_ptr& most, const amrex::BCRec* bc_ptr, - const int& qstate_size, bool vert_only = false); AMREX_GPU_DEVICE diff --git a/Source/Initialization/ERF_init_custom.cpp b/Source/Initialization/ERF_init_custom.cpp index 5ea10911a..e598f304e 100644 --- a/Source/Initialization/ERF_init_custom.cpp +++ b/Source/Initialization/ERF_init_custom.cpp @@ -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 diff --git a/Source/TimeIntegration/ERF_advance_dycore.cpp b/Source/TimeIntegration/ERF_advance_dycore.cpp index 4c14bc263..a47da7214 100644 --- a/Source/TimeIntegration/ERF_advance_dycore.cpp +++ b/Source/TimeIntegration/ERF_advance_dycore.cpp @@ -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); } // ***********************************************************************************************