Skip to content

Commit

Permalink
Implement smarter Vsg fix
Browse files Browse the repository at this point in the history
  • Loading branch information
ewquon committed Aug 23, 2024
1 parent 7fcded1 commit b4574fe
Show file tree
Hide file tree
Showing 2 changed files with 10 additions and 10 deletions.
2 changes: 1 addition & 1 deletion Source/BoundaryConditions/MOSTAverage.H
Original file line number Diff line number Diff line change
Expand Up @@ -220,6 +220,6 @@ protected:
bool include_subgrid_vel = false;
amrex::Real m_wstar_beta = 1.2; // Empirical parameter from Beljaars 1995 QJRMS
amrex::Vector<amrex::MultiFab*> m_wstar; // Ptr to 2D mf to hold wstar (maxlev)
amrex::Gpu::DeviceVector<amrex::Real> d_Vsg; // Subgrid velocity scale (Mahrt & Sun 1995 MWR)
amrex::Vector<amrex::Real> m_Vsg; // Subgrid velocity scale (Mahrt & Sun 1995 MWR)
};
#endif
18 changes: 9 additions & 9 deletions Source/BoundaryConditions/MOSTAverage.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -203,21 +203,19 @@ MOSTAverage::MOSTAverage (Vector<Geometry> geom,
// Corrections to the mean surface velocity
pp.query("most.wstar_beta", m_wstar_beta);
pp.query("most.include_subgrid_vel", include_subgrid_vel);
Vector<Real> Vsg(m_maxlev, 0.0);
m_Vsg = Vector<Real>(m_maxlev, 0.0);
if (include_subgrid_vel) {
amrex::Print() << "Subgrid velocity scale correction (by level) : ";
for (int lev(0); lev < m_maxlev; lev++) {
const auto dxArr = m_geom[lev].CellSizeArray();
Real dx = std::sqrt(dxArr[0]*dxArr[1]);
if (dx > 5000.) {
Vsg[lev] = 0.32 * std::pow(dx/5000.-1, 0.33);
m_Vsg[lev] = 0.32 * std::pow(dx/5000.-1, 0.33);
}
amrex::Print() << Vsg[lev];
amrex::Print() << m_Vsg[lev];
}
amrex::Print() << std::endl;
}
d_Vsg.resize(m_maxlev);
Gpu::copy(Gpu::hostToDevice, Vsg.begin(), Vsg.end(), d_Vsg.begin());
}


Expand Down Expand Up @@ -896,6 +894,7 @@ MOSTAverage::compute_plane_averages (int lev)
val_old[iavg] = plane_average[iavg]*d_fact_old;

const Real beta = m_wstar_beta;
const Real Vsg = m_Vsg[lev];

#ifdef _OPENMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
Expand Down Expand Up @@ -931,7 +930,7 @@ MOSTAverage::compute_plane_averages (int lev)
&v_interp, v_mf_arr, z_phys_arr, plo, dxInv, 1);
const Real wst = (wstar_arr) ? beta*wstar_arr(i,j,k) : 0.0;
const Real val = std::sqrt(u_interp*u_interp + v_interp*v_interp
+ wst*wst + d_Vsg[lev]*d_Vsg[lev]);
+ wst*wst + Vsg*Vsg);
Gpu::deviceReduceSum(&plane_avg[iavg], val, handler);
});
} else {
Expand All @@ -948,7 +947,7 @@ MOSTAverage::compute_plane_averages (int lev)
const Real v_val = 0.5 * (v_mf_arr(mi,mj,mk) + v_mf_arr(mi ,mj+1,mk));
const Real wst = (wstar_arr) ? beta*wstar_arr(i,j,k) : 0.0;
const Real val = std::sqrt(u_val*u_val + v_val*v_val
+ wst*wst + d_Vsg[lev]*d_Vsg[lev]);
+ wst*wst + Vsg*Vsg);
Gpu::deviceReduceSum(&plane_avg[iavg], val, handler);
});
}
Expand Down Expand Up @@ -1194,6 +1193,7 @@ MOSTAverage::compute_region_averages (int lev)
int iavg = m_navg - 1;

const Real beta = m_wstar_beta;
const Real Vsg = m_Vsg[lev];

#ifdef _OPENMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
Expand Down Expand Up @@ -1234,7 +1234,7 @@ MOSTAverage::compute_region_averages (int lev)
trilinear_interp_T(xp, yp, zp, &v_interp, v_mf_arr, z_phys_arr, plo, dxInv, 1);
const Real wst = (wstar_arr) ? beta*wstar_arr(i,j,k) : 0.0;
const Real mag = std::sqrt(u_interp*u_interp + v_interp*v_interp
+ wst*wst + d_Vsg[lev]*d_Vsg[lev]);
+ wst*wst + Vsg*Vsg);
Real val = denom * mag * d_fact_new;
ma_arr(i,j,k) += val;
}
Expand All @@ -1259,7 +1259,7 @@ MOSTAverage::compute_region_averages (int lev)
const Real v_val = 0.5 * (v_mf_arr(li,lj,lk) + v_mf_arr(li ,lj+1,lk));
const Real wst = (wstar_arr) ? beta*wstar_arr(i,j,k) : 0.0;
const Real mag = std::sqrt(u_val*u_val + v_val*v_val
+ wst*wst + d_Vsg[lev]*d_Vsg[lev]);
+ wst*wst + Vsg*Vsg);
Real val = denom * mag * d_fact_new;
ma_arr(i,j,k) += val;
}
Expand Down

0 comments on commit b4574fe

Please sign in to comment.