ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
Functions
ERF_ImmersedForcing.cpp File Reference
#include "ERF_ImmersedForcing.H"
#include "ERF_Constants.H"
#include "ERF_TI_slow_headers.H"
#include "ERF_SrcHeaders.H"
#include "ERF_ImmersedWallStability.H"
Include dependency graph for ERF_ImmersedForcing.cpp:

Functions

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real compute_if_most_target_vel (const amrex::Real u1_2r, const amrex::Real u2_2r, const amrex::Real delta, const amrex::Real z0, const amrex::Real t_blank, const amrex::Real theta_face, const amrex::Real theta_surf, const amrex::Real tflux_in, const amrex::Real Olen_in, const bool stability_correction, const amrex::Real wind_floor, const amrex::Real psi_cap_factor)
 
void ImmersedForcingTerrain_Xmom (const Box &tbx, const Array4< const Real > &u, const Array4< const Real > &v, const Array4< const Real > &w, const Array4< const Real > &cell_data, const Array4< const Real > &t_blank_arr, const Array4< const Real > &t_blank_xface_arr, const Array4< const Real > &z_cc_arr, const Array4< Real > &xmom_src_arr, const Geometry &geom, const SolverChoice &solverChoice, const Real fac)
 
void ImmersedForcingTerrain_Ymom (const Box &tby, const Array4< const Real > &u, const Array4< const Real > &v, const Array4< const Real > &w, const Array4< const Real > &cell_data, const Array4< const Real > &t_blank_arr, const Array4< const Real > &t_blank_yface_arr, const Array4< const Real > &z_cc_arr, const Array4< Real > &ymom_src_arr, const Geometry &geom, const SolverChoice &solverChoice, const Real fac)
 
void ImmersedForcingTerrain_Zmom (const Box &tbz, const Array4< const Real > &u, const Array4< const Real > &v, const Array4< const Real > &w, const Array4< const Real > &cell_data, const Array4< const Real > &t_blank_arr, const Array4< const Real > &t_blank_zface_arr, const Array4< const Real > &z_cc_arr, const Array4< Real > &zmom_src_arr, const Geometry &geom, const SolverChoice &solverChoice, const Real fac)
 
void ImmersedForcingBuildings_Xmom (const Box &tbx, const Array4< const Real > &u, const Array4< const Real > &v, const Array4< const Real > &w, const Array4< const Real > &cell_data, const Array4< const Real > &t_blank_arr, const Array4< const Real > &t_blank_xface_arr, const Array4< const Real > &z_cc_arr, const Array4< Real > &xmom_src_arr, const Geometry &geom, const SolverChoice &solverChoice, const Real fac)
 
void ImmersedForcingBuildings_Ymom (const Box &tby, const Array4< const Real > &u, const Array4< const Real > &v, const Array4< const Real > &w, const Array4< const Real > &cell_data, const Array4< const Real > &t_blank_arr, const Array4< const Real > &t_blank_yface_arr, const Array4< const Real > &z_cc_arr, const Array4< Real > &ymom_src_arr, const Geometry &geom, const SolverChoice &solverChoice, const Real fac)
 
void ImmersedForcingBuildings_Zmom (const Box &tbz, const Array4< const Real > &u, const Array4< const Real > &v, const Array4< const Real > &w, const Array4< const Real > &cell_data, const Array4< const Real > &t_blank_arr, const Array4< const Real > &t_blank_zface_arr, const Array4< const Real > &z_cc_arr, const Array4< Real > &zmom_src_arr, const Geometry &geom, const SolverChoice &solverChoice, const Real fac)
 
void ImmersedForcingTerrain_Scalar (const Box &bx, const Array4< const Real > &u, const Array4< const Real > &v, const Array4< const Real > &cell_data, const Array4< const Real > &t_blank_arr, const Array4< const Real > &z_cc_arr, const Array4< Real > &cell_src, const Geometry &geom, const SolverChoice &solverChoice, const Table1D< Real > &r_avg, const Table1D< Real > &t_avg, const Real time)
 
void ImmersedForcingBuildings_Scalar (const Box &bx, const Array4< const Real > &u, const Array4< const Real > &v, const Array4< const Real > &w, const Array4< const Real > &cell_data, const Array4< const Real > &t_blank_arr, const Array4< const Real > &z_cc_arr, const Array4< Real > &cell_src, const Geometry &geom, const SolverChoice &solverChoice, const Table1D< Real > &r_avg, const Table1D< Real > &t_avg, const Real time)
 

Function Documentation

◆ compute_if_most_target_vel()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real compute_if_most_target_vel ( const amrex::Real  u1_2r,
const amrex::Real  u2_2r,
const amrex::Real  delta,
const amrex::Real  z0,
const amrex::Real  t_blank,
const amrex::Real  theta_face,
const amrex::Real  theta_surf,
const amrex::Real  tflux_in,
const amrex::Real  Olen_in,
const bool  stability_correction,
const amrex::Real  wind_floor,
const amrex::Real  psi_cap_factor 
)

Compute the target velocity using Monin-Obukhov Similarity Theory for immersed forcing.

Parameters
[in]u1_2rFirst tangential velocity component.
[in]u2_2rSecond tangential velocity component.
[in]deltaDistance from the surface.
[in]z0Roughness length.
[in]t_blankVolume fraction.
[in]theta_facePotential temperature at the face.
[in]theta_surfPotential temperature at the surface.
[in]tflux_inSurface heat flux.
[in]Olen_inObukhov length.
[in]stability_correctionWhether to apply stability corrections.
[in]wind_floorWind-speed floor [m/s] behind the friction velocity of the stability estimate (erf.if_stability_wind_floor).
[in]psi_cap_factorCap on psi_m as a fraction of ln(z/z0) (erf.if_psi_cap_factor).
Returns
Target velocity component.
43 {
44  const Real tiny = std::numeric_limits<amrex::Real>::epsilon();
45  Real psi_m = zero;
46  Real psi_h = zero;
47  Real tang_windspeed2r = std::sqrt(u1_2r * u1_2r + u2_2r * u2_2r);
48 
49  // The stability estimate uses the wind floored at wind_floor, a derived Obukhov length bounded
50  // so that |zeta| <= ib_stability::zeta_max() as on flat ground, and a prescribed one as given
51  // with its zeta held there; the target below uses the actual wind
52  Real ustar = ib_stability::floored_wind(tang_windspeed2r, wind_floor) * KAPPA / (std::log(1.5 * delta / z0) - psi_m);
53  Real tflux = (tflux_in != Real(1.e-8)) ? tflux_in : -(theta_face - theta_surf) * ustar * KAPPA / (std::log(1.5 * delta / z0) - psi_h);
54  Real Olen = (Olen_in != Real(1.e-8)) ? Olen_in
55  : ib_stability::bounded_obukhov_length(-ustar * ustar * ustar * theta_face / (KAPPA * CONST_GRAV * tflux + tiny), Real(1.5) * delta);
56  Real zeta = ib_stability::bounded_zeta(Real(1.5) * delta, Olen);
57 
58  // similarity functions
59  similarity_funs sfuns;
60  if (stability_correction){
61  psi_m = sfuns.calc_psi_m(zeta);
62  psi_h = sfuns.calc_psi_h(zeta);
63  }
64  ustar = tang_windspeed2r * KAPPA / (std::log(1.5 * delta / z0) - ib_stability::psi_m_for_ustar(psi_m, Real(1.5) * delta, z0, psi_cap_factor));
65 
66  // prevent some unphysical math
67  ustar = ib_stability::clamped_ustar(ustar);
68  psi_m = ib_stability::capped_psi(psi_m, myhalf * delta, z0, psi_cap_factor);
69 
70  Real uTarget = (1 - t_blank) * ustar / KAPPA * (std::log(myhalf * delta / z0) - psi_m);
71  Real u1Target = uTarget * u1_2r / (tiny + tang_windspeed2r);
72 
73  return u1Target;
74 }
KAPPA
constexpr amrex::Real KAPPA
Definition: ERF_Constants.H:67
CONST_GRAV
constexpr amrex::Real CONST_GRAV
Definition: ERF_Constants.H:68
z0
Real z0
Definition: ERF_InitCustomPertVels_ScalarAdvDiff.H:8
zero
constexpr amrex::Real zero
Definition: ERF_NumericalConstants.H:36
myhalf
constexpr amrex::Real myhalf
Definition: ERF_NumericalConstants.H:41
Real
amrex::Real Real
Definition: ERF_ShocInterface.H:19
LsmVar_SLM::ustar
@ ustar
Definition: ERF_SLM.H:28
ib_stability::bounded_zeta
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real bounded_zeta(amrex::Real z, amrex::Real Olen) noexcept
Definition: ERF_ImmersedWallStability.H:62
ib_stability::floored_wind
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real floored_wind(amrex::Real speed, amrex::Real wind_floor) noexcept
Definition: ERF_ImmersedWallStability.H:38
ib_stability::capped_psi
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real capped_psi(amrex::Real psi, amrex::Real z, amrex::Real z0, amrex::Real factor) noexcept
Definition: ERF_ImmersedWallStability.H:85
ib_stability::bounded_obukhov_length
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real bounded_obukhov_length(amrex::Real Olen, amrex::Real z_ref) noexcept
Definition: ERF_ImmersedWallStability.H:51
ib_stability::psi_m_for_ustar
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real psi_m_for_ustar(amrex::Real psi_m, amrex::Real z, amrex::Real z0, amrex::Real factor) noexcept
Definition: ERF_ImmersedWallStability.H:98
ib_stability::clamped_ustar
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real clamped_ustar(amrex::Real ustar) noexcept
Definition: ERF_ImmersedWallStability.H:109
module_model_constants::epsilon
real(c_double), parameter epsilon
Definition: ERF_module_model_constants.F90:12
similarity_funs
Definition: ERF_MOSTUtils.H:37
similarity_funs::calc_psi_m
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real calc_psi_m(amrex::Real zeta) const
Definition: ERF_MOSTUtils.H:102
similarity_funs::calc_psi_h
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real calc_psi_h(amrex::Real zeta) const
Definition: ERF_MOSTUtils.H:121

Referenced by ImmersedForcingBuildings_Xmom(), ImmersedForcingBuildings_Ymom(), and ImmersedForcingBuildings_Zmom().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ ImmersedForcingBuildings_Scalar()

void ImmersedForcingBuildings_Scalar ( const Box &  bx,
const Array4< const Real > &  u,
const Array4< const Real > &  v,
const Array4< const Real > &  w,
const Array4< const Real > &  cell_data,
const Array4< const Real > &  t_blank_arr,
const Array4< const Real > &  z_cc_arr,
const Array4< Real > &  cell_src,
const Geometry &  geom,
const SolverChoice &  solverChoice,
const Table1D< Real > &  r_avg,
const Table1D< Real > &  t_avg,
const Real  time 
)

Apply buildings immersed forcing to scalars (Rho, RhoTheta)

1132 {
1133  // geometric properties
1134  const Real* dx_arr = geom.CellSize();
1135  const Real dx_x = dx_arr[0];
1136  const Real dx_y = dx_arr[1];
1137 
1138  const Real alpha_h = solverChoice.if_Cd_scalar;
1139  // erf.if_snap_partial_cells: read the cell blanking snapped to solid (1)
1140  // or fluid (0) at half, so a height-map building becomes the same
1141  // staircase of whole cells an exact box is: the thermal conditions sit
1142  // on the boundary solid cells, roofs included, found from the neighbour
1143  // blanking. Off (the default), the raw fractions are used and nothing
1144  // below changes. See SolverChoice::if_snap_partial_cells.
1145  const bool l_snap = solverChoice.if_snap_partial_cells;
1146  auto snapb = [=] AMREX_GPU_DEVICE (amrex::Real b) noexcept -> amrex::Real {
1147  return l_snap ? ((b >= myhalf) ? one : zero) : b;
1148  };
1149  const Real tiny = std::numeric_limits<amrex::Real>::epsilon();
1150  const Real U_s = one; // unit velocity scale
1151 
1152  // MOST parameters
1153  similarity_funs sfuns;
1154  const Real ggg = CONST_GRAV;
1155  const Real kappa = KAPPA;
1156  const Real z0 = solverChoice.if_z0;
1157  const Real tflux = solverChoice.if_surf_temp_flux;
1158  const Real init_surf_temp = solverChoice.if_init_surf_temp;
1159  const Real surf_heating_rate = solverChoice.if_surf_heating_rate;
1160  const Real stab_wind_floor = solverChoice.if_stability_wind_floor;
1161 
1162  ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
1163  {
1164  const Real t_blank = snapb(t_blank_arr(i, j, k));
1165  const Real t_blank_below = snapb(t_blank_arr(i, j, k-1));
1166  const Real t_blank_above = snapb(t_blank_arr(i, j, k+1));
1167  const Real t_blank_north = snapb(t_blank_arr(i , j+1, k));
1168  const Real t_blank_south = snapb(t_blank_arr(i , j-1, k));
1169  const Real t_blank_east = snapb(t_blank_arr(i+1, j , k));
1170  const Real t_blank_west = snapb(t_blank_arr(i-1, j , k));
1171 
1172  const Real dx_z = (z_cc_arr) ? (z_cc_arr(i,j,k) - z_cc_arr(i,j,k-1)) : dx_arr[2];
1173  Real drag_coefficient = alpha_h / std::pow(dx_x*dx_y*dx_z, one/three);
1174 
1175  // Wall cells of the building, named by the face they carry: a
1176  // partially blanked cell with a more solid neighbour behind it and a
1177  // fluid one in front, or, with the snap (every solid cell at 1), a
1178  // solid cell with a solid neighbour behind and a fluid one in front.
1179  const bool south_face = l_snap ? (t_blank == one && t_blank_north == one && t_blank_south == zero)
1180  : (t_blank > zero && t_blank < t_blank_north && t_blank_south == zero);
1181  const bool north_face = l_snap ? (t_blank == one && t_blank_south == one && t_blank_north == zero)
1182  : (t_blank > zero && t_blank < t_blank_south && t_blank_north == zero);
1183  const bool west_face = l_snap ? (t_blank == one && t_blank_east == one && t_blank_west == zero)
1184  : (t_blank > zero && t_blank < t_blank_east && t_blank_west == zero);
1185  const bool east_face = l_snap ? (t_blank == one && t_blank_west == one && t_blank_east == zero)
1186  : (t_blank > zero && t_blank < t_blank_west && t_blank_east == zero);
1187 
1188  // SURFACE TEMP AND HEATING/COOLING RATE
1189  if (init_surf_temp > zero) {
1190  const Real surf_temp = init_surf_temp + surf_heating_rate*time;
1191  if (t_blank > 0 && (t_blank_above == zero) && (t_blank_below == one)) { // building roof
1192  const Real bc_forcing_rt_srf = -(cell_data(i,j,k,Rho_comp) * surf_temp - cell_data(i,j,k,RhoTheta_comp));
1193  cell_src(i, j, k, RhoTheta_comp) -= drag_coefficient * U_s * bc_forcing_rt_srf;
1194 
1195  } else if (east_face || west_face || south_face || north_face) {
1196  // this should enter for just building walls
1197  // walls are currently separated to allow for flexibility in the future to heat walls differently
1198 
1199  // south face
1200  if (l_snap ? south_face : ((t_blank < t_blank_north) && (t_blank_north == one))) {
1201  const Real bc_forcing_rt_srf = -(cell_data(i,j,k,Rho_comp) * surf_temp - cell_data(i,j,k,RhoTheta_comp));
1202  cell_src(i, j, k, RhoTheta_comp) -= drag_coefficient * U_s * bc_forcing_rt_srf;
1203  }
1204 
1205  // north face
1206  if (l_snap ? north_face : ((t_blank < t_blank_south) && (t_blank_south == one))) {
1207  const Real bc_forcing_rt_srf = -(cell_data(i,j,k,Rho_comp) * surf_temp - cell_data(i,j,k,RhoTheta_comp));
1208  cell_src(i, j, k, RhoTheta_comp) -= drag_coefficient * U_s * bc_forcing_rt_srf;
1209  }
1210 
1211  // west face
1212  if (l_snap ? west_face : ((t_blank < t_blank_east) && (t_blank_east == one))) {
1213  const Real bc_forcing_rt_srf = -(cell_data(i,j,k,Rho_comp) * surf_temp - cell_data(i,j,k,RhoTheta_comp));
1214  cell_src(i, j, k, RhoTheta_comp) -= drag_coefficient * U_s * bc_forcing_rt_srf;
1215  }
1216 
1217  // east face
1218  if (l_snap ? east_face : ((t_blank < t_blank_west) && (t_blank_west == one))) {
1219  const Real bc_forcing_rt_srf = -(cell_data(i,j,k,Rho_comp) * surf_temp - cell_data(i,j,k,RhoTheta_comp));
1220  cell_src(i, j, k, RhoTheta_comp) -= drag_coefficient * U_s * bc_forcing_rt_srf;
1221  }
1222 
1223  }
1224  }
1225 
1226  // SURFACE HEAT FLUX
1227  if (tflux != Real(1.e-8)){
1228  const Real ux_cc_2r = myhalf * (u(i ,j ,k+1) + u(i+1,j ,k+1));
1229  const Real uy_cc_2r = myhalf * (v(i ,j ,k+1) + v(i ,j+1,k+1));
1230  // The heat transfer uses the wind floored at erf.if_stability_wind_floor and an Obukhov
1231  // length bounded so that |zeta| <= ib_stability::zeta_max(), as on flat ground
1232  const Real h_windspeed2r = ib_stability::floored_wind(std::sqrt(ux_cc_2r * ux_cc_2r + uy_cc_2r * uy_cc_2r), stab_wind_floor);
1233 
1234  const Real theta = cell_data(i,j,k ,RhoTheta_comp) / cell_data(i,j,k ,Rho_comp);
1235  Real theta_neighbor = cell_data(i,j,k+1,RhoTheta_comp) / cell_data(i,j,k+1,Rho_comp);
1236 
1237  if (t_blank > zero && (t_blank_above == zero)) { // building roof
1238  Real psi_m = zero;
1239  Real psi_h = zero;
1240  Real psi_h_neighbor = zero;
1241  Real ustar = h_windspeed2r * kappa / (std::log((1.5) * dx_z / z0) - psi_m);
1242  // (start-up refuses erf.if_Olen with a surface flux, so L always comes from the flux here)
1243  Real Olen = ib_stability::bounded_obukhov_length(-ustar * ustar * ustar * theta / (kappa * ggg * tflux + tiny), Real(1.5) * dx_z);
1244 
1245  for (int iter = 0; iter < 2; ++iter) {
1246  if (iter > 0) { Olen = ib_stability::bounded_obukhov_length(-ustar * ustar * ustar * theta / (kappa * ggg * tflux + tiny), Real(1.5) * dx_z); }
1247  Real zeta = ib_stability::bounded_zeta(myhalf * dx_z, Olen);
1248  Real zeta_neighbor = ib_stability::bounded_zeta(Real(1.5) * dx_z, Olen);
1249 
1250  // similarity functions
1251  psi_m = sfuns.calc_psi_m(zeta);
1252  psi_h = sfuns.calc_psi_h(zeta);
1253  psi_h_neighbor = sfuns.calc_psi_h(zeta_neighbor);
1254  ustar = h_windspeed2r * kappa / (std::log((1.5) * dx_z / z0) - psi_m);
1255  }
1256 
1257  // prevent some unphysical math
1258  ustar = ib_stability::clamped_ustar(ustar);
1259  // caps of ln(z / z0): erf.if_psi_cap_factor stays off the temperature forcing, whose
1260  // theta* comes from an L not made consistent with u* (the legacy law, kept as it is;
1261  // the fraction-stress law of erf-model/ERF#4134 forms them consistently)
1262  psi_h_neighbor = ib_stability::capped_psi(psi_h_neighbor, Real(1.5) * dx_z, z0, one);
1263  psi_h = ib_stability::capped_psi(psi_h, myhalf * dx_z, z0, one);
1264 
1265  // We do not know the actual temperature so use cell above
1266  const Real thetastar = theta * ustar * ustar / (kappa * ggg * Olen);
1267  const Real surf_temp = theta_neighbor - thetastar / kappa * (std::log((1.5) * dx_z / z0) - psi_h_neighbor);
1268  const Real tTarget = surf_temp + thetastar / kappa * (std::log((myhalf) * dx_z / z0) - psi_h);
1269 
1270  const Real bc_forcing_rt = -(cell_data(i,j,k,Rho_comp) * tTarget - cell_data(i,j,k,RhoTheta_comp));
1271  cell_src(i, j, k, RhoTheta_comp) -= drag_coefficient * U_s * bc_forcing_rt;
1272 
1273  } else if (east_face || west_face || south_face || north_face) { // this should enter for just building walls
1274 
1275  Real ux_cellaway = zero;
1276  Real uy_cellaway = zero;
1277  Real uz_cellaway = zero;
1278  Real u1 = zero;
1279  Real u2 = zero;
1280  Real delta = zero;
1281 
1282  // south face
1283  if (south_face) {
1284  ux_cellaway = myhalf * (u(i ,j-1,k) + u(i+1,j-1,k ));
1285  uz_cellaway = myhalf * (w(i ,j-1,k) + w(i ,j-1,k+1));
1286  u1 = ux_cellaway;
1287  u2 = uz_cellaway;
1288  delta = dx_y;
1289 
1290  // MOST
1291  theta_neighbor = cell_data(i,j-1,k,RhoTheta_comp) / cell_data(i,j-1,k,Rho_comp);
1292  }
1293 
1294  // north face
1295  if (north_face) {
1296  ux_cellaway = myhalf * (u(i ,j+1,k) + u(i+1,j+1,k ));
1297  uz_cellaway = myhalf * (w(i ,j+1,k) + w(i ,j+1,k+1));
1298  u1 = ux_cellaway;
1299  u2 = uz_cellaway;
1300  delta = dx_y;
1301 
1302  // MOST
1303  theta_neighbor = cell_data(i,j+1,k,RhoTheta_comp) / cell_data(i,j+1,k,Rho_comp);
1304  }
1305 
1306  // west face
1307  if (west_face) {
1308  uy_cellaway = myhalf * (v(i-1,j ,k) + v(i-1,j+1,k ));
1309  uz_cellaway = myhalf * (w(i-1,j ,k) + w(i-1,j ,k+1));
1310  u1 = uy_cellaway;
1311  u2 = uz_cellaway;
1312  delta = dx_x;
1313 
1314  // MOST
1315  theta_neighbor = cell_data(i-1,j,k,RhoTheta_comp) / cell_data(i-1,j,k,Rho_comp);
1316  }
1317 
1318  // east face
1319  if (east_face) {
1320  uy_cellaway = myhalf * (v(i+1,j ,k) + v(i+1,j+1,k ));
1321  uz_cellaway = myhalf * (w(i+1,j ,k) + w(i+1,j ,k+1));
1322  u1 = uy_cellaway;
1323  u2 = uz_cellaway;
1324  delta = dx_x;
1325 
1326  // MOST
1327  theta_neighbor = cell_data(i+1,j,k,RhoTheta_comp) / cell_data(i+1,j,k,Rho_comp);
1328  }
1329 
1330  Real tan_wspd = ib_stability::floored_wind(std::sqrt(u1 * u1 + u2 * u2), stab_wind_floor);
1331 
1332  Real psi_m = zero;
1333  Real psi_h = zero;
1334  Real psi_h_neighbor = zero;
1335  Real ustar = tan_wspd * kappa / (std::log(1.5 * delta / z0) - psi_m);
1336  Real Olen = ib_stability::bounded_obukhov_length(-ustar * ustar * ustar * theta / (kappa * ggg * tflux + tiny), Real(1.5) * delta);
1337 
1338  for (int iter = 0; iter < 2; ++iter) {
1339  if (iter > 0) { Olen = ib_stability::bounded_obukhov_length(-ustar * ustar * ustar * theta / (kappa * ggg * tflux + tiny), Real(1.5) * delta); }
1340  Real zeta = ib_stability::bounded_zeta(myhalf * delta, Olen);
1341  Real zeta_neighbor = ib_stability::bounded_zeta(Real(1.5) * delta, Olen);
1342 
1343  // similarity functions
1344  psi_m = sfuns.calc_psi_m(zeta);
1345  psi_h = sfuns.calc_psi_h(zeta);
1346  psi_h_neighbor = sfuns.calc_psi_h(zeta_neighbor);
1347  ustar = tan_wspd * kappa / (std::log((1.5) * delta / z0) - psi_m);
1348  }
1349 
1350  // prevent some unphysical math
1351  ustar = ib_stability::clamped_ustar(ustar);
1352  // caps of ln(z / z0), as above (erf-model/ERF#4134)
1353  psi_h_neighbor = ib_stability::capped_psi(psi_h_neighbor, Real(1.5) * delta, z0, one);
1354  psi_h = ib_stability::capped_psi(psi_h, myhalf * delta, z0, one);
1355 
1356  // We do not know the actual temperature so use cell above
1357  const Real thetastar = theta * ustar * ustar / (kappa * ggg * Olen);
1358  const Real surf_temp = theta_neighbor - thetastar / kappa * (std::log((1.5) * delta / z0) - psi_h_neighbor);
1359  const Real tTarget = surf_temp + thetastar / kappa * (std::log((myhalf) * delta / z0) - psi_h);
1360 
1361  const Real bc_forcing_rt = -(cell_data(i,j,k,Rho_comp) * tTarget - cell_data(i,j,k,RhoTheta_comp));
1362  cell_src(i, j, k, RhoTheta_comp) -= drag_coefficient * U_s * bc_forcing_rt;
1363  }
1364  }
1365 
1366  // Force fully immersed cells to planar average rho and theta
1367  if (t_blank == 1.0 && r_avg && t_avg) {
1368  const Real rho_avg = r_avg(k);
1369  const Real theta_avg = t_avg(k) / rho_avg; // Convert from RhoTheta to Theta
1370  const Real rho_cell = cell_data(i,j,k,Rho_comp);
1371  const Real bc_forcing_r = -(rho_avg - rho_cell);
1372  const Real bc_forcing_rt = -(rho_avg * theta_avg - cell_data(i,j,k,RhoTheta_comp));
1373 
1374  cell_src(i, j, k, Rho_comp) -= drag_coefficient * U_s * bc_forcing_r;
1375  cell_src(i, j, k, RhoTheta_comp) -= drag_coefficient * U_s * bc_forcing_rt;
1376  }
1377  });
1378 }
Rho_comp
#define Rho_comp
Definition: ERF_IndexDefines.H:42
RhoTheta_comp
#define RhoTheta_comp
Definition: ERF_IndexDefines.H:43
ParallelFor
ParallelFor(fab_box, [=] AMREX_GPU_DEVICE(int i, int j, int k) { qrcuten_arr(i, j, k)=Real(0);qscuten_arr(i, j, k)=Real(0);qicuten_arr(i, j, k)=Real(0);})
three
constexpr amrex::Real three
Definition: ERF_NumericalConstants.H:39
one
constexpr amrex::Real one
Definition: ERF_NumericalConstants.H:37
w
Real w
Definition: ERF_Plotfile2DInterpolator.cpp:22
MicVar_Kess::theta
@ theta
Definition: ERF_Kessler.H:26
SolverChoice::if_z0
amrex::Real if_z0
Immersed-forcing roughness length [m].
Definition: ERF_DataStruct.H:2323
SolverChoice::if_snap_partial_cells
bool if_snap_partial_cells
Definition: ERF_DataStruct.H:2329
SolverChoice::if_Cd_scalar
amrex::Real if_Cd_scalar
Immersed-forcing drag coefficient for scalars.
Definition: ERF_DataStruct.H:2318
SolverChoice::if_init_surf_temp
amrex::Real if_init_surf_temp
Initial immersed-forcing surface temperature [K].
Definition: ERF_DataStruct.H:2325
SolverChoice::if_surf_temp_flux
amrex::Real if_surf_temp_flux
Immersed-forcing surface temperature flux [K m/s].
Definition: ERF_DataStruct.H:2324
SolverChoice::if_surf_heating_rate
amrex::Real if_surf_heating_rate
Immersed-forcing surface heating rate [K/hr].
Definition: ERF_DataStruct.H:2326
SolverChoice::if_stability_wind_floor
amrex::Real if_stability_wind_floor
Wind-speed floor [m/s] behind the friction velocity of the wall law's stability estimate (0: none).
Definition: ERF_DataStruct.H:2333

Referenced by make_sources().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ ImmersedForcingBuildings_Xmom()

void ImmersedForcingBuildings_Xmom ( const Box &  tbx,
const Array4< const Real > &  u,
const Array4< const Real > &  v,
const Array4< const Real > &  w,
const Array4< const Real > &  cell_data,
const Array4< const Real > &  t_blank_arr,
const Array4< const Real > &  t_blank_xface_arr,
const Array4< const Real > &  z_cc_arr,
const Array4< Real > &  xmom_src_arr,
const Geometry &  geom,
const SolverChoice &  solverChoice,
const Real  fac 
)

Apply buildings immersed forcing to X-momentum

375 {
376  // geometric properties
377  const Real* dx_arr = geom.CellSize();
378  const Real dx_x = dx_arr[0];
379  const Real dx_y = dx_arr[1];
380  const Real dt = fac;
381 
382  const Real alpha_m = solverChoice.if_Cd_momentum;
383  const Real tiny = std::numeric_limits<amrex::Real>::epsilon();
384  const Real U_s = one; // unit velocity scale
385 
386  // MOST parameters
387  const Real z0 = solverChoice.if_z0;
388  const Real tflux_in = solverChoice.if_surf_temp_flux;
389  const Real Olen_in = solverChoice.if_Olen_in;
390  const Real stab_wind_floor = solverChoice.if_stability_wind_floor;
391  const Real psi_cap_factor = solverChoice.if_psi_cap_factor;
392  const bool l_use_most = solverChoice.if_use_most;
393  const bool l_stability_correction = solverChoice.if_stability_correction;
394 
395  // To limit stiffness of drag when using anelastic
396  const Real ws_floor = solverChoice.if_ws_floor;
397  const Real damp_alpha = solverChoice.if_damp_alpha;
398  // Point-implicit alternative to the clamp above; stabilizes both compressible and anelastic
399  const bool l_implicit_drag = solverChoice.if_implicit_drag;
400 
401  const bool is_slow_step = true; // This is determined by calling context
402  const bool use_ImmersedForcing_fast = solverChoice.immersed_forcing_substep;
403  const Real small_volfrac = 0.005;
404  // erf.if_snap_partial_cells: read the cell blanking snapped to solid (1)
405  // or fluid (0) at half, so a height-map building becomes the same
406  // staircase of whole cells an exact box is. A face is solid when either
407  // cell it joins is; a face between a solid and a fluid cell is
408  // wall-normal and gets the interior drag toward zero (no penetration); a
409  // face between two solid cells carries the roof or wall law of its row
410  // (the full log-law target, not the partial-cell weighted one) or the
411  // interior drag, never both; the partial-cell branches (wall_mask,
412  // east_west_mask and the like, faces with 0 < t_blank < 1) do not arise
413  // under the snap. Off (the default), the raw fractions
414  // are used and nothing below changes: a boundary face of an exact box
415  // carries the wall law and the interior drag together as before. See
416  // SolverChoice::if_snap_partial_cells.
417  const bool l_snap = solverChoice.if_snap_partial_cells;
418  auto snapb = [=] AMREX_GPU_DEVICE (amrex::Real b) noexcept -> amrex::Real {
419  return l_snap ? ((b >= myhalf) ? one : zero) : b;
420  };
421  // Blanking of an x-face (i, j, k): with the snap, from the two cells it
422  // joins (solid when either is), so both builds give the same staircase;
423  // otherwise the face-centred fraction when the build has it, else the
424  // mean of the two cells.
425  auto fb = [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept -> amrex::Real {
426  if (l_snap) { return amrex::max(snapb(t_blank_arr(i, j, k)), snapb(t_blank_arr(i-1, j, k))); }
427  return (t_blank_xface_arr) ? t_blank_xface_arr(i, j, k)
428  : myhalf * (t_blank_arr(i, j, k) + t_blank_arr(i-1, j, k));
429  };
430 
431 
432  ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
433  {
434  const Real ux = u(i, j, k );
435  const Real uy = fourth * ( v(i, j , k ) + v(i-1, j , k )
436  + v(i, j+1, k ) + v(i-1, j+1, k ) );
437  const Real uz = fourth * ( w(i, j , k ) + w(i-1, j , k )
438  + w(i, j , k+1) + w(i-1, j , k+1) );
439  const amrex::Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz);
440 
441  const Real rho_xface = myhalf * ( cell_data(i,j,k,Rho_comp) + cell_data(i-1,j,k,Rho_comp) );
442  const Real theta_xface = (myhalf * (cell_data(i,j,k,RhoTheta_comp) + cell_data(i-1,j,k, RhoTheta_comp))) / rho_xface;
443 
444  // Use face-centered terrain_blanking if available, otherwise average from cell centers with threshold
445  Real t_blank_raw = fb(i, j, k);
446  const Real t_blank = (t_blank_raw < small_volfrac) ? zero : t_blank_raw;
447  // With the snap the boundary solid face stands for the wall layer and
448  // takes the full log-law target; the partial-cell weight (1 - t_blank)
449  // of compute_if_most_target_vel() would make it zero (a no-slip
450  // staircase), so the snap path passes a zero blanking to it.
451  const Real t_blank_law = l_snap ? zero : t_blank;
452 
453  Real t_blank_below_raw = (k == 0) ? zero : fb(i, j, k-1);
454  const Real t_blank_below = (t_blank_below_raw < small_volfrac) ? zero : t_blank_below_raw;
455 
456  Real t_blank_above_raw = fb(i, j, k+1);
457  const Real t_blank_above = (t_blank_above_raw < small_volfrac) ? zero : t_blank_above_raw;
458 
459  Real t_blank_north_raw = fb(i, j+1, k);
460  const Real t_blank_north = (t_blank_north_raw < small_volfrac) ? zero : t_blank_north_raw;
461 
462  Real t_blank_south_raw = fb(i, j-1, k);
463  const Real t_blank_south = (t_blank_south_raw < small_volfrac) ? zero : t_blank_south_raw;
464 
465  const Real dx_z = (z_cc_arr) ? (z_cc_arr(i,j,k) - z_cc_arr(i,j,k-1)) : dx_arr[2];
466  const Real drag_coefficient = alpha_m / std::pow(dx_x*dx_y*dx_z, one/three);
467  const Real CdM = std::min(drag_coefficient / (windspeed + tiny), drag_coefficient);
468  // With the snap a face joining a solid and a fluid cell is wall-normal:
469  // it gets the interior drag (no penetration), not a wall law, whatever
470  // row it lies in.
471  const bool normal_face = l_snap && (snapb(t_blank_arr(i, j, k)) != snapb(t_blank_arr(i-1, j, k)));
472 
473  // With the snap every solid face has t_blank = 1: a roof face lies in
474  // the top solid row (t_blank <= t_blank_below, the face above fluid),
475  // and a face carrying a wall law is not also an interior face.
476  const Real roof_mask = (!normal_face && t_blank > zero && (l_snap ? t_blank <= t_blank_below : t_blank < t_blank_below) && t_blank_above == zero && l_use_most) ? one : zero; // roof cell
477  const Real south_mask = (!normal_face && t_blank > zero && t_blank <= t_blank_north && t_blank_south == zero && l_use_most) ? one : zero; // south wall cell
478  const Real north_mask = (!normal_face && t_blank > zero && t_blank <= t_blank_south && t_blank_north == zero && l_use_most) ? one : zero; // north wall cell
479  const Real wall_mask = (t_blank > zero && t_blank < one && !l_use_most) ? one : zero; // all walls when NOT using MOST
480  const Real most_mask = roof_mask + south_mask + north_mask; // cells getting MOST treatment
481  const Real east_west_mask = (t_blank > zero && t_blank < one && l_use_most && most_mask == zero) ? one : zero; // partial cells not covered by MOST (east/west walls)
482  const Real interior_mask = (t_blank == 1.0 && !(l_snap && most_mask > zero)) ? one : zero; // interior cell
483 
484  Real drag = zero;
485  Real u1_cellaway = zero;
486  Real u2_cellaway = zero;
487  Real rho_xface_inside = rho_xface;
488  Real theta_surf = theta_xface;
489  Real bc_forcing_x = zero;
490  Real u_target = zero;
491 
492  // roof forcing
493  if (roof_mask == one) {
494  u1_cellaway = u(i, j, k+1) ;
495  u2_cellaway = fourth * ( v(i, j , k+1) + v(i-1, j , k+1)
496  + v(i, j+1, k+1) + v(i-1, j+1, k+1) ) ;
497  rho_xface_inside = myhalf * (cell_data(i,j,k-1,Rho_comp) + cell_data(i-1,j,k-1,Rho_comp));
498  theta_surf = (myhalf * (cell_data(i,j,k-1,RhoTheta_comp) + cell_data(i-1,j,k-1, RhoTheta_comp))) / rho_xface_inside;
499  u_target = compute_if_most_target_vel(u1_cellaway, u2_cellaway, dx_z, z0, t_blank_law, theta_xface, theta_surf, tflux_in, Olen_in, l_stability_correction, stab_wind_floor, psi_cap_factor);
500  bc_forcing_x = -(u_target - ux); // BC forcing pushes nonrelative velocity toward target velocity
501  drag += bc_forcing_x * roof_mask * rho_xface * CdM * U_s;
502  }
503 
504  // south wall forcing
505  if (south_mask == one) {
506  u1_cellaway = u(i, j-1, k );
507  u2_cellaway = fourth * ( w(i, j-1, k ) + w(i-1, j-1, k )
508  + w(i, j-1, k+1) + w(i-1, j-1, k+1) ) ;
509  rho_xface_inside = myhalf * ( cell_data(i,j+1,k,Rho_comp) + cell_data(i-1,j+1,k,Rho_comp) );
510  theta_surf = (myhalf * (cell_data(i,j+1,k,RhoTheta_comp) + cell_data(i-1,j+1,k, RhoTheta_comp))) / rho_xface_inside;
511  u_target = compute_if_most_target_vel(u1_cellaway, u2_cellaway, dx_y, z0, t_blank_law, theta_xface, theta_surf, tflux_in, Olen_in, l_stability_correction, stab_wind_floor, psi_cap_factor);
512  bc_forcing_x = -(u_target - ux); // BC forcing pushes nonrelative velocity toward target velocity
513  drag += bc_forcing_x * south_mask * rho_xface * CdM * U_s;
514  }
515 
516  // north wall forcing
517  if (north_mask == one) {
518  u1_cellaway = u(i, j+1, k ) ;
519  u2_cellaway = fourth * ( w(i, j+1, k ) + w(i-1, j+1, k )
520  + w(i, j+1, k+1) + w(i-1, j+1, k+1) ) ;
521  rho_xface_inside = myhalf * ( cell_data(i,j-1,k,Rho_comp) + cell_data(i-1,j-1,k,Rho_comp) );
522  theta_surf = (myhalf * (cell_data(i,j-1,k,RhoTheta_comp) + cell_data(i-1,j-1,k, RhoTheta_comp))) / rho_xface_inside;
523  u_target = compute_if_most_target_vel(u1_cellaway, u2_cellaway, dx_y, z0, t_blank_law, theta_xface, theta_surf, tflux_in, Olen_in, l_stability_correction, stab_wind_floor, psi_cap_factor);
524  bc_forcing_x = -(u_target - ux); // BC forcing pushes nonrelative velocity toward target velocity
525  drag += bc_forcing_x * north_mask * rho_xface * CdM * U_s;
526  }
527 
528  // wall forcing (if not using most) or east/west walls when using MOST
529  if (wall_mask == one || east_west_mask == one) {
530  drag += (wall_mask + east_west_mask) * t_blank * rho_xface * CdM * ux * windspeed;
531  }
532 
533  // interior cell forcing
534  if (interior_mask == one) {
535  drag += interior_mask * rho_xface * CdM * ux * windspeed;
536  }
537 
538  if (l_implicit_drag) {
539  // point-implicit rescale of the aggregated drag
540  const Real lambda = CdM * ( (roof_mask + south_mask + north_mask) * U_s
541  + (wall_mask + east_west_mask) * t_blank * windspeed
542  + interior_mask * windspeed );
543  xmom_src_arr(i,j,k) -= drag / (one + lambda*dt);
544  } else if (is_slow_step && !use_ImmersedForcing_fast) {
545  // limit drag term for anelastic for numerical stability
546  Real d_drag = dt * -drag; // time step * acceleration like tendency
547  Real wsmax_change = damp_alpha * amrex::max(amrex::Math::abs(ux), ws_floor); // aims to prevent oscillations around 0.
548  if (amrex::Math::abs(ux) < 0.1){ // no damping for smaller velocities
549  wsmax_change =one * amrex::max(amrex::Math::abs(ux), ws_floor);
550  }
551  d_drag = amrex::min(amrex::max(d_drag, -wsmax_change), wsmax_change);
552  xmom_src_arr(i,j,k) += d_drag / dt; // put back as limited tendency
553  } else {
554  xmom_src_arr(i, j, k) -= drag;
555  }
556  });
557 }
compute_if_most_target_vel
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real compute_if_most_target_vel(const amrex::Real u1_2r, const amrex::Real u2_2r, const amrex::Real delta, const amrex::Real z0, const amrex::Real t_blank, const amrex::Real theta_face, const amrex::Real theta_surf, const amrex::Real tflux_in, const amrex::Real Olen_in, const bool stability_correction, const amrex::Real wind_floor, const amrex::Real psi_cap_factor)
Definition: ERF_ImmersedForcing.cpp:29
fourth
constexpr amrex::Real fourth
Definition: ERF_NumericalConstants.H:42
SolverChoice::if_psi_cap_factor
amrex::Real if_psi_cap_factor
Cap on psi_m in the immersed momentum wall law as a fraction of ln(z/z0) (1: none).
Definition: ERF_DataStruct.H:2334
SolverChoice::if_Cd_momentum
amrex::Real if_Cd_momentum
Immersed-forcing drag coefficient for momentum.
Definition: ERF_DataStruct.H:2317
SolverChoice::if_Olen_in
amrex::Real if_Olen_in
Input Obukhov length for immersed-forcing MOST [m].
Definition: ERF_DataStruct.H:2327
SolverChoice::if_damp_alpha
amrex::Real if_damp_alpha
Immersed-forcing damping coefficient.
Definition: ERF_DataStruct.H:2332
SolverChoice::if_ws_floor
amrex::Real if_ws_floor
Wind-speed floor [m/s] in the per-step change limit of the building wall-law forcing.
Definition: ERF_DataStruct.H:2331
SolverChoice::immersed_forcing_substep
bool immersed_forcing_substep
Whether immersed-forcing source terms are applied only during substeps.
Definition: ERF_DataStruct.H:2313
SolverChoice::if_use_most
bool if_use_most
Whether immersed-forcing MOST is enabled.
Definition: ERF_DataStruct.H:2328
SolverChoice::if_implicit_drag
bool if_implicit_drag
Definition: ERF_DataStruct.H:2321
SolverChoice::if_stability_correction
bool if_stability_correction
Whether immersed-forcing stability corrections are enabled.
Definition: ERF_DataStruct.H:2330

Referenced by make_mom_sources().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ ImmersedForcingBuildings_Ymom()

void ImmersedForcingBuildings_Ymom ( const Box &  tby,
const Array4< const Real > &  u,
const Array4< const Real > &  v,
const Array4< const Real > &  w,
const Array4< const Real > &  cell_data,
const Array4< const Real > &  t_blank_arr,
const Array4< const Real > &  t_blank_yface_arr,
const Array4< const Real > &  z_cc_arr,
const Array4< Real > &  ymom_src_arr,
const Geometry &  geom,
const SolverChoice &  solverChoice,
const Real  fac 
)

Apply buildings immersed forcing to Y-momentum

574 {
575  // geometric properties
576  const Real* dx_arr = geom.CellSize();
577  const Real dx_x = dx_arr[0];
578  const Real dx_y = dx_arr[1];
579  const Real dt = fac;
580 
581  const Real alpha_m = solverChoice.if_Cd_momentum;
582  const Real tiny = std::numeric_limits<amrex::Real>::epsilon();
583  const Real U_s = one; // unit velocity scale
584 
585  // MOST parameters
586  const Real z0 = solverChoice.if_z0;
587  const Real tflux_in = solverChoice.if_surf_temp_flux;
588  const Real Olen_in = solverChoice.if_Olen_in;
589  const Real stab_wind_floor = solverChoice.if_stability_wind_floor;
590  const Real psi_cap_factor = solverChoice.if_psi_cap_factor;
591  const bool l_use_most = solverChoice.if_use_most;
592  const bool l_stability_correction = solverChoice.if_stability_correction;
593 
594  // To limit stiffness of drag when using anelastic
595  const Real ws_floor = solverChoice.if_ws_floor;
596  const Real damp_alpha = solverChoice.if_damp_alpha;
597  // Point-implicit alternative to the clamp above; stabilizes both compressible and anelastic
598  const bool l_implicit_drag = solverChoice.if_implicit_drag;
599 
600  const bool is_slow_step = true; // This is determined by calling context
601  const bool use_ImmersedForcing_fast = solverChoice.immersed_forcing_substep;
602  const Real small_volfrac = 0.005;
603  // erf.if_snap_partial_cells: read the cell blanking snapped to solid (1)
604  // or fluid (0) at half, so a height-map building becomes the same
605  // staircase of whole cells an exact box is. A face is solid when either
606  // cell it joins is; a face between a solid and a fluid cell is
607  // wall-normal and gets the interior drag toward zero (no penetration); a
608  // face between two solid cells carries the roof or wall law of its row
609  // (the full log-law target, not the partial-cell weighted one) or the
610  // interior drag, never both; the partial-cell branches (wall_mask,
611  // east_west_mask and the like, faces with 0 < t_blank < 1) do not arise
612  // under the snap. Off (the default), the raw fractions
613  // are used and nothing below changes: a boundary face of an exact box
614  // carries the wall law and the interior drag together as before. See
615  // SolverChoice::if_snap_partial_cells.
616  const bool l_snap = solverChoice.if_snap_partial_cells;
617  auto snapb = [=] AMREX_GPU_DEVICE (amrex::Real b) noexcept -> amrex::Real {
618  return l_snap ? ((b >= myhalf) ? one : zero) : b;
619  };
620  // Blanking of a y-face (i, j, k): with the snap, from the two cells it
621  // joins (solid when either is), so both builds give the same staircase;
622  // otherwise the face-centred fraction when the build has it, else the
623  // mean of the two cells.
624  auto fb = [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept -> amrex::Real {
625  if (l_snap) { return amrex::max(snapb(t_blank_arr(i, j, k)), snapb(t_blank_arr(i, j-1, k))); }
626  return (t_blank_yface_arr) ? t_blank_yface_arr(i, j, k)
627  : myhalf * (t_blank_arr(i, j, k) + t_blank_arr(i, j-1, k));
628  };
629 
630 
631  ParallelFor(tby, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
632  {
633  const Real ux = fourth * ( u(i , j , k ) + u(i , j-1, k )
634  + u(i+1, j , k ) + u(i+1, j-1, k ) );
635  const Real uy = v(i, j, k);
636  const Real uz = fourth * ( w(i , j , k ) + w(i , j-1, k )
637  + w(i , j , k+1) + w(i , j-1, k+1) );
638  const amrex::Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz);
639 
640  const Real rho_yface = myhalf * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j-1,k,Rho_comp) );
641  const Real theta_yface = (myhalf * (cell_data(i,j,k ,RhoTheta_comp) + cell_data(i,j-1,k,RhoTheta_comp))) / rho_yface;
642 
643  // Use face-centered terrain_blanking if available, otherwise average from cell centers with threshold
644  Real t_blank_raw = fb(i, j, k);
645  const Real t_blank = (t_blank_raw < small_volfrac) ? zero : t_blank_raw;
646  // With the snap the boundary solid face stands for the wall layer and
647  // takes the full log-law target; the partial-cell weight (1 - t_blank)
648  // of compute_if_most_target_vel() would make it zero (a no-slip
649  // staircase), so the snap path passes a zero blanking to it.
650  const Real t_blank_law = l_snap ? zero : t_blank;
651 
652  Real t_blank_below_raw = (k == 0) ? zero : fb(i, j, k-1);
653  const Real t_blank_below = (t_blank_below_raw < small_volfrac) ? zero : t_blank_below_raw;
654 
655  Real t_blank_above_raw = fb(i, j, k+1);
656  const Real t_blank_above = (t_blank_above_raw < small_volfrac) ? zero : t_blank_above_raw;
657 
658  Real t_blank_east_raw = fb(i+1, j, k);
659  const Real t_blank_east = (t_blank_east_raw < small_volfrac) ? zero : t_blank_east_raw;
660 
661  Real t_blank_west_raw = fb(i-1, j, k);
662  const Real t_blank_west = (t_blank_west_raw < small_volfrac) ? zero : t_blank_west_raw;
663 
664  const Real dx_z = (z_cc_arr) ? (z_cc_arr(i,j,k) - z_cc_arr(i,j,k-1)) : dx_arr[2];
665  const Real drag_coefficient = alpha_m / std::pow(dx_x*dx_y*dx_z, one/three);
666  const Real CdM = std::min(drag_coefficient / (windspeed + tiny), drag_coefficient);
667  // With the snap a face joining a solid and a fluid cell is wall-normal:
668  // it gets the interior drag (no penetration), not a wall law, whatever
669  // row it lies in.
670  const bool normal_face = l_snap && (snapb(t_blank_arr(i, j, k)) != snapb(t_blank_arr(i, j-1, k)));
671 
672  // As in the x-momentum: with the snap a roof face lies in the top
673  // solid row and a wall-law face is not also an interior face.
674  const Real roof_mask = (!normal_face && t_blank > zero && (l_snap ? t_blank <= t_blank_below : t_blank < t_blank_below) && t_blank_above == zero && l_use_most) ? one : zero; // roof cell
675  const Real west_mask = (!normal_face && t_blank > zero && t_blank <= t_blank_east && t_blank_west == zero && l_use_most) ? one : zero; // west wall cell
676  const Real east_mask = (!normal_face && t_blank > zero && t_blank <= t_blank_west && t_blank_east == zero && l_use_most) ? one : zero; // east wall cell
677  const Real wall_mask = (t_blank > zero && t_blank < one && !l_use_most) ? one : zero; // all walls when NOT using MOST
678  const Real most_mask = roof_mask + west_mask + east_mask; // cells getting MOST treatment
679  const Real north_south_mask = (t_blank > zero && t_blank < one && l_use_most && most_mask == zero) ? one : zero; // partial cells not covered by MOST (north/south walls)
680  const Real interior_mask = (t_blank == 1.0 && !(l_snap && most_mask > zero)) ? one : zero; // interior cell
681 
682  Real drag = zero;
683  Real u1_cellaway = zero;
684  Real u2_cellaway = zero;
685  Real rho_yface_inside = rho_yface;
686  Real theta_surf = theta_yface;
687  Real bc_forcing_y = zero;
688  Real u_target = zero;
689 
690  // roof forcing
691  if (roof_mask == one) {
692  u1_cellaway = fourth * ( u(i , j , k+1) + u(i , j-1, k+1)
693  + u(i+1, j , k+1) + u(i+1, j-1, k+1) );
694  u2_cellaway = v(i, j, k+1);
695  rho_yface_inside = myhalf * ( cell_data(i,j,k-1,Rho_comp) + cell_data(i,j-1,k-1,Rho_comp) );
696  theta_surf = (myhalf * (cell_data(i,j,k-1,RhoTheta_comp) + cell_data(i,j-1,k-1,RhoTheta_comp))) / rho_yface_inside;
697  u_target = compute_if_most_target_vel(u1_cellaway, u2_cellaway, dx_z, z0, t_blank_law, theta_yface, theta_surf, tflux_in, Olen_in, l_stability_correction, stab_wind_floor, psi_cap_factor);
698  bc_forcing_y = -(u_target - uy); // BC forcing pushes nonrelative velocity toward target velocity
699  drag += bc_forcing_y * roof_mask * rho_yface * CdM * U_s;
700  }
701 
702  // west wall forcing
703  if (west_mask == one) {
704  u1_cellaway = v(i-1, j , k );
705  u2_cellaway = fourth * ( w(i-1, j , k ) + w(i-1, j-1, k )
706  + w(i-1, j , k+1) + w(i-1, j-1, k+1) );
707  rho_yface_inside = myhalf * ( cell_data(i+1,j,k,Rho_comp) + cell_data(i+1,j-1,k,Rho_comp) );
708  theta_surf = (myhalf * (cell_data(i+1,j,k,RhoTheta_comp) + cell_data(i+1,j-1,k,RhoTheta_comp))) / rho_yface_inside;
709  u_target = compute_if_most_target_vel(u1_cellaway, u2_cellaway, dx_x, z0, t_blank_law, theta_yface, theta_surf, tflux_in, Olen_in, l_stability_correction, stab_wind_floor, psi_cap_factor);
710  bc_forcing_y = -(u_target - uy); // BC forcing pushes nonrelative velocity toward target velocity
711  drag += bc_forcing_y * west_mask * rho_yface * CdM * U_s;
712  }
713 
714  // east wall forcing
715  if (east_mask == one) {
716  u1_cellaway = v(i+1, j , k );
717  u2_cellaway = fourth * ( w(i+1, j , k ) + w(i+1, j-1, k )
718  + w(i+1, j , k+1) + w(i+1, j-1, k+1) );
719  rho_yface_inside = myhalf * ( cell_data(i-1,j,k,Rho_comp) + cell_data(i-1,j-1,k,Rho_comp) );
720  theta_surf = (myhalf * (cell_data(i-1,j,k,RhoTheta_comp) + cell_data(i-1,j-1,k,RhoTheta_comp))) / rho_yface_inside;
721  u_target = compute_if_most_target_vel(u1_cellaway, u2_cellaway, dx_x, z0, t_blank_law, theta_yface, theta_surf, tflux_in, Olen_in, l_stability_correction, stab_wind_floor, psi_cap_factor);
722  bc_forcing_y = -(u_target - uy); // BC forcing pushes nonrelative velocity toward target velocity
723  drag += bc_forcing_y * east_mask * rho_yface * CdM * U_s;
724  }
725 
726  // wall forcing (if not using most) or north/south walls when using MOST
727  if (wall_mask == one || north_south_mask == one) {
728  drag += (wall_mask + north_south_mask) * t_blank * rho_yface * CdM * uy * windspeed;
729  }
730 
731  // interior cell forcing
732  if (interior_mask == one) {
733  drag += interior_mask * rho_yface * CdM * uy * windspeed;
734  }
735 
736  if (l_implicit_drag) {
737  // point-implicit rescale of the aggregated drag
738  const Real lambda = CdM * ( (roof_mask + west_mask + east_mask) * U_s
739  + (wall_mask + north_south_mask) * t_blank * windspeed
740  + interior_mask * windspeed );
741  ymom_src_arr(i,j,k) -= drag / (one + lambda*dt);
742  } else if (is_slow_step && !use_ImmersedForcing_fast) {
743  // limit drag term for anelastic for numerical stability
744  Real d_drag = dt * -drag; // time step * acceleration like tendency
745  Real wsmax_change = damp_alpha * amrex::max(amrex::Math::abs(uy), ws_floor); // aims to prevent oscillations around 0.
746  if (amrex::Math::abs(uy) < 0.1){ // no damping for smaller velocities
747  wsmax_change =one * amrex::max(amrex::Math::abs(uy), ws_floor);
748  }
749  d_drag = amrex::min(amrex::max(d_drag, -wsmax_change), wsmax_change);
750  ymom_src_arr(i,j,k) += d_drag / dt; // put back as limited tendency
751  } else {
752  ymom_src_arr(i, j, k) -= drag;
753  }
754  });
755 }

Referenced by make_mom_sources().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ ImmersedForcingBuildings_Zmom()

void ImmersedForcingBuildings_Zmom ( const Box &  tbz,
const Array4< const Real > &  u,
const Array4< const Real > &  v,
const Array4< const Real > &  w,
const Array4< const Real > &  cell_data,
const Array4< const Real > &  t_blank_arr,
const Array4< const Real > &  t_blank_zface_arr,
const Array4< const Real > &  z_cc_arr,
const Array4< Real > &  zmom_src_arr,
const Geometry &  geom,
const SolverChoice &  solverChoice,
const Real  fac 
)

Apply buildings immersed forcing to Z-momentum

772 {
773  // geometric properties
774  const Real* dx_arr = geom.CellSize();
775  const Real dx_x = dx_arr[0];
776  const Real dx_y = dx_arr[1];
777  const Real dt = fac;
778 
779  const Real alpha_m = solverChoice.if_Cd_momentum;
780  const Real tiny = std::numeric_limits<amrex::Real>::epsilon();
781  const Real U_s = one; // unit velocity scale
782 
783  // MOST parameters
784  const Real z0 = solverChoice.if_z0;
785  const Real tflux_in = solverChoice.if_surf_temp_flux;
786  const Real Olen_in = solverChoice.if_Olen_in;
787  const Real stab_wind_floor = solverChoice.if_stability_wind_floor;
788  const Real psi_cap_factor = solverChoice.if_psi_cap_factor;
789  const bool l_use_most = solverChoice.if_use_most;
790  const bool l_stability_correction = solverChoice.if_stability_correction;
791 
792  // To limit stiffness of drag when using anelastic
793  const Real ws_floor = solverChoice.if_ws_floor;
794  const Real damp_alpha = solverChoice.if_damp_alpha;
795  // Point-implicit alternative to the clamp above; stabilizes both compressible and anelastic
796  const bool l_implicit_drag = solverChoice.if_implicit_drag;
797 
798  const bool is_slow_step = true; // This is determined by calling context
799  const bool use_ImmersedForcing_fast = solverChoice.immersed_forcing_substep;
800  const Real small_volfrac = 0.005;
801  // erf.if_snap_partial_cells: read the cell blanking snapped to solid (1)
802  // or fluid (0) at half, so a height-map building becomes the same
803  // staircase of whole cells an exact box is. A face is solid when either
804  // cell it joins is; a face between a solid and a fluid cell is
805  // wall-normal and gets the interior drag toward zero (no penetration); a
806  // face between two solid cells carries the roof or wall law of its row
807  // (the full log-law target, not the partial-cell weighted one) or the
808  // interior drag, never both; the partial-cell branches (wall_mask,
809  // east_west_mask and the like, faces with 0 < t_blank < 1) do not arise
810  // under the snap. Off (the default), the raw fractions
811  // are used and nothing below changes: a boundary face of an exact box
812  // carries the wall law and the interior drag together as before. See
813  // SolverChoice::if_snap_partial_cells.
814  const bool l_snap = solverChoice.if_snap_partial_cells;
815  auto snapb = [=] AMREX_GPU_DEVICE (amrex::Real b) noexcept -> amrex::Real {
816  return l_snap ? ((b >= myhalf) ? one : zero) : b;
817  };
818  // Blanking of a z-face (i, j, k): with the snap, from the two cells it
819  // joins (solid when either is), so both builds give the same staircase;
820  // otherwise the face-centred fraction when the build has it, else the
821  // mean of the two cells.
822  auto fb = [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept -> amrex::Real {
823  if (l_snap) { return amrex::max(snapb(t_blank_arr(i, j, k)), snapb(t_blank_arr(i, j, k-1))); }
824  return (t_blank_zface_arr) ? t_blank_zface_arr(i, j, k)
825  : myhalf * (t_blank_arr(i, j, k) + t_blank_arr(i, j, k-1));
826  };
827 
828 
829  ParallelFor(tbz, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
830  {
831  const Real ux = fourth * ( u(i , j , k ) + u(i+1, j , k )
832  + u(i , j , k-1) + u(i+1, j , k-1) );
833  const Real uy = fourth * ( v(i, j , k ) + v(i, j+1, k )
834  + v(i, j , k-1) + v(i, j+1, k-1) );
835  const Real uz = w(i, j, k);
836  const amrex::Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz);
837 
838  const Real rho_zface = myhalf * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j,k-1,Rho_comp) );
839  const Real theta_zface = (myhalf * (cell_data(i,j,k,RhoTheta_comp) + cell_data(i,j,k-1,RhoTheta_comp))) / rho_zface;
840 
841  // Use face-centered terrain_blanking if available, otherwise average from cell centers with threshold
842  Real t_blank_raw = fb(i, j, k);
843  const Real t_blank = (t_blank_raw < small_volfrac) ? zero : t_blank_raw;
844  // With the snap the boundary solid face stands for the wall layer and
845  // takes the full log-law target; the partial-cell weight (1 - t_blank)
846  // of compute_if_most_target_vel() would make it zero (a no-slip
847  // staircase), so the snap path passes a zero blanking to it.
848  const Real t_blank_law = l_snap ? zero : t_blank;
849 
850  Real t_blank_above_raw = fb(i, j, k+1);
851  const Real t_blank_above = (t_blank_above_raw < small_volfrac) ? zero : t_blank_above_raw;
852 
853  Real t_blank_north_raw = fb(i, j+1, k);
854  const Real t_blank_north = (t_blank_north_raw < small_volfrac) ? zero : t_blank_north_raw;
855 
856  Real t_blank_south_raw = fb(i, j-1, k);
857  const Real t_blank_south = (t_blank_south_raw < small_volfrac) ? zero : t_blank_south_raw;
858 
859  Real t_blank_east_raw = fb(i+1, j, k);
860  const Real t_blank_east = (t_blank_east_raw < small_volfrac) ? zero : t_blank_east_raw;
861 
862  Real t_blank_west_raw = fb(i-1, j, k);
863  const Real t_blank_west = (t_blank_west_raw < small_volfrac) ? zero : t_blank_west_raw;
864 
865  const Real dx_z = (z_cc_arr) ? (z_cc_arr(i,j,k) - z_cc_arr(i,j,k-1)) : dx_arr[2];
866  const Real drag_coefficient = alpha_m / std::pow(dx_x*dx_y*dx_z, one/three);
867  const Real CdM = std::min(drag_coefficient / (windspeed + tiny), drag_coefficient);
868  // With the snap a face joining a solid and a fluid cell is wall-normal:
869  // it gets the interior drag (no penetration), not a wall law, whatever
870  // row it lies in.
871  const bool normal_face = l_snap && (snapb(t_blank_arr(i, j, k)) != snapb(t_blank_arr(i, j, k-1)));
872 
873  const Real south_mask = (!normal_face && t_blank > zero && t_blank <= t_blank_north && t_blank_south == zero && l_use_most && k >= 1) ? one : zero; // south wall cell
874  const Real north_mask = (!normal_face && t_blank > zero && t_blank <= t_blank_south && t_blank_north == zero && l_use_most && k >= 1) ? one : zero; // north wall cell
875  const Real west_mask = (!normal_face && t_blank > zero && t_blank <= t_blank_east && t_blank_west == zero && l_use_most && k >= 1) ? one : zero; // west wall cell
876  const Real east_mask = (!normal_face && t_blank > zero && t_blank <= t_blank_west && t_blank_east == zero && l_use_most && k >= 1) ? one : zero; // east wall cell
877  const Real wall_mask = (t_blank > zero && t_blank < one && !l_use_most) ? one : zero; // all walls when NOT using MOST
878  const Real roof_mask = (!normal_face && t_blank > zero && t_blank_above == zero && l_use_most) ? one : zero; // roof cell (horizontal surface) - uses simple drag
879  // With the snap a face carrying a wall law or the roof drag is not also an interior face.
880  const Real most_mask = south_mask + north_mask + west_mask + east_mask + roof_mask;
881  const Real interior_mask = (t_blank == 1.0 && !(l_snap && most_mask > zero)) ? one : zero; // interior cell
882 
883  Real drag = zero;
884  Real u1_cellaway = zero;
885  Real u2_cellaway = zero;
886  Real rho_zface_inside = rho_zface;
887  Real theta_surf = theta_zface;
888  Real bc_forcing_z = zero;
889  Real u_target = zero;
890 
891  // south wall forcing
892  if (south_mask == one) {
893  u1_cellaway = fourth * ( u(i , j-1, k ) + u(i+1, j-1, k )
894  + u(i , j-1, k-1) + u(i+1, j-1, k-1) );
895  u2_cellaway = w(i, j-1, k);
896  rho_zface_inside = myhalf * ( cell_data(i,j+1,k,Rho_comp) + cell_data(i,j+1,k-1,Rho_comp) );
897  theta_surf = (myhalf * (cell_data(i,j+1,k,RhoTheta_comp) + cell_data(i,j+1,k-1,RhoTheta_comp))) / rho_zface_inside;
898  u_target = compute_if_most_target_vel(u1_cellaway, u2_cellaway, dx_y, z0, t_blank_law, theta_zface, theta_surf, tflux_in, Olen_in, l_stability_correction, stab_wind_floor, psi_cap_factor);
899  bc_forcing_z = -(u_target - uz); // BC forcing pushes nonrelative velocity toward target velocity
900  drag += bc_forcing_z * south_mask * rho_zface * CdM * U_s;
901  }
902 
903  // north wall forcing
904  if (north_mask == one) {
905  u1_cellaway = fourth * ( u(i , j+1, k ) + u(i+1, j+1, k )
906  + u(i , j+1, k-1) + u(i+1, j+1, k-1) );
907  u2_cellaway = w(i, j+1, k);
908  rho_zface_inside = myhalf * ( cell_data(i,j-1,k,Rho_comp) + cell_data(i,j-1,k-1,Rho_comp) );
909  theta_surf = (myhalf * (cell_data(i,j-1,k,RhoTheta_comp) + cell_data(i,j-1,k-1,RhoTheta_comp))) / rho_zface_inside;
910  u_target = compute_if_most_target_vel(u1_cellaway, u2_cellaway, dx_y, z0, t_blank_law, theta_zface, theta_surf, tflux_in, Olen_in, l_stability_correction, stab_wind_floor, psi_cap_factor);
911  bc_forcing_z = -(u_target - uz); // BC forcing pushes nonrelative velocity toward target velocity
912  drag += bc_forcing_z * north_mask * rho_zface * CdM * U_s;
913  }
914 
915  // west wall forcing
916  if (west_mask == one) {
917  u1_cellaway = fourth * ( v(i-1, j , k ) + v(i-1, j+1, k )
918  + v(i-1, j , k-1) + v(i-1, j+1, k-1) );
919  u2_cellaway = w(i-1, j, k);
920  rho_zface_inside = myhalf * ( cell_data(i+1,j,k,Rho_comp) + cell_data(i+1,j,k-1,Rho_comp) );
921  theta_surf = (myhalf * (cell_data(i+1,j,k,RhoTheta_comp) + cell_data(i+1,j,k-1,RhoTheta_comp))) / rho_zface_inside;
922  u_target = compute_if_most_target_vel(u1_cellaway, u2_cellaway, dx_x, z0, t_blank_law, theta_zface, theta_surf, tflux_in, Olen_in, l_stability_correction, stab_wind_floor, psi_cap_factor);
923  bc_forcing_z = -(u_target - uz); // BC forcing pushes nonrelative velocity toward target velocity
924  drag += bc_forcing_z * west_mask * rho_zface * CdM * U_s;
925  }
926 
927  // east wall forcing
928  if (east_mask == one) {
929  u1_cellaway = fourth * ( v(i+1, j , k ) + v(i+1, j+1, k )
930  + v(i+1, j , k-1) + v(i+1, j+1, k-1) );
931  u2_cellaway = w(i+1, j, k);
932  rho_zface_inside = myhalf * ( cell_data(i-1,j,k,Rho_comp) + cell_data(i-1,j,k-1,Rho_comp) );
933  theta_surf = (myhalf * (cell_data(i-1,j,k,RhoTheta_comp) + cell_data(i-1,j,k-1,RhoTheta_comp))) / rho_zface_inside;
934  u_target = compute_if_most_target_vel(u1_cellaway, u2_cellaway, dx_x, z0, t_blank_law, theta_zface, theta_surf, tflux_in, Olen_in, l_stability_correction, stab_wind_floor, psi_cap_factor);
935  bc_forcing_z = -(u_target - uz); // BC forcing pushes nonrelative velocity toward target velocity
936  drag += bc_forcing_z * east_mask * rho_zface * CdM * U_s;
937  }
938 
939  // wall forcing (if not using most) or roof when using MOST
940  if (wall_mask == one || roof_mask == one) {
941  drag += (wall_mask + roof_mask) * t_blank * rho_zface * CdM * uz * windspeed;
942  }
943 
944  // interior cell forcing
945  if (interior_mask == one) {
946  drag += interior_mask * rho_zface * CdM * uz * windspeed;
947  }
948 
949  if (l_implicit_drag) {
950  // point-implicit rescale of the aggregated drag
951  const Real lambda = CdM * ( (south_mask + north_mask + west_mask + east_mask) * U_s
952  + (wall_mask + roof_mask) * t_blank * windspeed
953  + interior_mask * windspeed );
954  zmom_src_arr(i,j,k) -= drag / (one + lambda*dt);
955  } else if (is_slow_step && !use_ImmersedForcing_fast) {
956  // limit drag term for anelastic for numerical stability
957  Real d_drag = dt * -drag; // time step * acceleration like tendency
958  Real wsmax_change = damp_alpha * amrex::max(amrex::Math::abs(uz), ws_floor); // aims to prevent oscillations around 0.
959  if (amrex::Math::abs(uz) < 0.1){ // no damping for smaller velocities
960  wsmax_change = one * amrex::max(amrex::Math::abs(uz), ws_floor);
961  }
962  d_drag = amrex::min(amrex::max(d_drag, -wsmax_change), wsmax_change);
963  zmom_src_arr(i,j,k) += d_drag / dt; // put back as limited tendency
964  } else {
965  zmom_src_arr(i, j, k) -= drag;
966  }
967  });
968 }

Referenced by make_mom_sources().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ ImmersedForcingTerrain_Scalar()

void ImmersedForcingTerrain_Scalar ( const Box &  bx,
const Array4< const Real > &  u,
const Array4< const Real > &  v,
const Array4< const Real > &  cell_data,
const Array4< const Real > &  t_blank_arr,
const Array4< const Real > &  z_cc_arr,
const Array4< Real > &  cell_src,
const Geometry &  geom,
const SolverChoice &  solverChoice,
const Table1D< Real > &  r_avg,
const Table1D< Real > &  t_avg,
const Real  time 
)

Apply terrain immersed forcing to scalars (Rho, RhoTheta)

985 {
986  // geometric properties
987  const Real* dx_arr = geom.CellSize();
988  const Real dx_x = dx_arr[0];
989  const Real dx_y = dx_arr[1];
990 
991  const Real alpha_h = solverChoice.if_Cd_scalar;
992  // The terrain kernels keep the raw fractions: erf.if_snap_partial_cells
993  // applies to the buildings kernels only (the wall law below is weighted
994  // by the fluid fraction of the face, which a snapped face has not).
995  const Real tiny = std::numeric_limits<amrex::Real>::epsilon();
996  const Real U_s = one; // unit velocity scale
997 
998  // MOST parameters
999  similarity_funs sfuns;
1000  const Real ggg = CONST_GRAV;
1001  const Real kappa = KAPPA;
1002  const Real z0 = solverChoice.if_z0;
1003  const Real tflux = solverChoice.if_surf_temp_flux;
1004  const Real init_surf_temp = solverChoice.if_init_surf_temp;
1005 
1006  // Note this has been converted to K / s when it was read in;
1007  const Real surf_heating_rate = solverChoice.if_surf_heating_rate;
1008 
1009  const Real Olen_in = solverChoice.if_Olen_in;
1010  const Real stab_wind_floor = solverChoice.if_stability_wind_floor;
1011 
1012  ParallelFor(bx, [=]
1013  AMREX_GPU_DEVICE(int i, int j, int k) noexcept
1014  {
1015  const Real dx_z = (z_cc_arr) ? (z_cc_arr(i,j,k) - z_cc_arr(i,j,k-1)) : dx_arr[2];
1016  const Real drag_coefficient = alpha_h / std::pow(dx_x*dx_y*dx_z, one/three);
1017 
1018  const Real t_blank = t_blank_arr(i, j, k);
1019  const Real t_blank_above = t_blank_arr(i, j, k+1);
1020  const Real ux_cc_2r = myhalf * (u(i ,j ,k+1) + u(i+1,j ,k+1));
1021  const Real uy_cc_2r = myhalf * (v(i ,j ,k+1) + v(i ,j+1,k+1));
1022  // The heat transfer uses the wind floored at erf.if_stability_wind_floor and an Obukhov
1023  // length bounded so that |zeta| <= ib_stability::zeta_max(), as on flat ground
1024  const Real h_windspeed2r = ib_stability::floored_wind(std::sqrt(ux_cc_2r * ux_cc_2r + uy_cc_2r * uy_cc_2r), stab_wind_floor);
1025 
1026  const Real theta = cell_data(i,j,k ,RhoTheta_comp) / cell_data(i,j,k ,Rho_comp);
1027  const Real theta_neighbor = cell_data(i,j,k+1,RhoTheta_comp) / cell_data(i,j,k+1,Rho_comp);
1028 
1029  // SURFACE TEMP AND HEATING/COOLING RATE
1030  if (init_surf_temp > zero) {
1031  if (t_blank > 0 && (t_blank_above == zero)) { // force to MOST value
1032  const Real surf_temp = init_surf_temp + surf_heating_rate*time;
1033  const Real bc_forcing_rt_srf = -(cell_data(i,j,k-1,Rho_comp) * surf_temp - cell_data(i,j,k-1,RhoTheta_comp));
1034  cell_src(i, j, k, RhoTheta_comp) -= drag_coefficient * U_s * bc_forcing_rt_srf;
1035  }
1036  }
1037 
1038  // SURFACE HEAT FLUX
1039  if (tflux != Real(1e-8)){
1040  if (t_blank > 0 && (t_blank_above == zero)) { // force to MOST value
1041  Real psi_m = zero;
1042  Real psi_h = zero;
1043  Real psi_h_neighbor = zero;
1044  Real ustar = h_windspeed2r * kappa / (std::log((Real(1.5)) * dx_z / z0) - psi_m);
1045  const Real Olen = ib_stability::bounded_obukhov_length(-ustar * ustar * ustar * theta / (kappa * ggg * tflux + tiny),
1046  Real(1.5) * dx_z);
1047  const Real zeta = ib_stability::bounded_zeta(myhalf * dx_z, Olen);
1048  const Real zeta_neighbor = ib_stability::bounded_zeta(Real(1.5) * dx_z, Olen);
1049 
1050  // similarity functions
1051  psi_m = sfuns.calc_psi_m(zeta);
1052  psi_h = sfuns.calc_psi_h(zeta);
1053  psi_h_neighbor = sfuns.calc_psi_h(zeta_neighbor);
1054  ustar = h_windspeed2r * kappa / (std::log((Real(1.5)) * dx_z / z0) - psi_m);
1055 
1056  // prevent some unphysical math
1057  ustar = ib_stability::clamped_ustar(ustar);
1058  // caps of ln(z / z0): erf.if_psi_cap_factor stays off the temperature forcing, whose
1059  // theta* comes from an L not made consistent with u* (the legacy law, kept as it is;
1060  // the fraction-stress law of erf-model/ERF#4134 forms them consistently)
1061  psi_h_neighbor = ib_stability::capped_psi(psi_h_neighbor, Real(1.5) * dx_z, z0, one);
1062  psi_h = ib_stability::capped_psi(psi_h, myhalf * dx_z, z0, one);
1063 
1064  // We do not know the actual temperature so use cell above
1065  const Real thetastar = theta * ustar * ustar / (kappa * ggg * Olen);
1066  const Real surf_temp = theta_neighbor - thetastar / kappa * (std::log((Real(1.5)) * dx_z / z0) - psi_h_neighbor);
1067  const Real tTarget = surf_temp + thetastar / kappa * (std::log((myhalf) * dx_z / z0) - psi_h);
1068 
1069  const Real bc_forcing_rt = -(cell_data(i,j,k,Rho_comp) * tTarget - cell_data(i,j,k,RhoTheta_comp));
1070  cell_src(i, j, k, RhoTheta_comp) -= drag_coefficient * U_s * bc_forcing_rt;
1071  }
1072  }
1073 
1074  // OBUKHOV LENGTH
1075  if (Olen_in != Real(1e-8)){
1076  if (t_blank > 0 && (t_blank_above == zero)) { // force to MOST value
1077  // the prescribed length as given, with its zeta held within +-zeta_max()
1078  const Real Olen = Olen_in;
1079  const Real zeta = ib_stability::bounded_zeta(myhalf * dx_z, Olen);
1080  const Real zeta_neighbor = ib_stability::bounded_zeta(Real(1.5) * dx_z, Olen);
1081 
1082  // similarity functions; psi_h capped at ln(z / z0) as on the other temperature branches
1083  const Real psi_m = sfuns.calc_psi_m(zeta);
1084  const Real psi_h = ib_stability::capped_psi(sfuns.calc_psi_h(zeta), myhalf * dx_z, z0, one);
1085  const Real psi_h_neighbor = ib_stability::capped_psi(sfuns.calc_psi_h(zeta_neighbor), Real(1.5) * dx_z, z0, one);
1086  // psi_m at 0.5 dz against ln(1.5 dz / z0) can cancel; u* is held in [0, 2] m/s as on the
1087  // other branches. In unstable air psi_h >= psi_m and psi_h grows with |zeta|, so wherever
1088  // psi_m reaches ln(1.5 dz / z0) both psi_h caps above bind, the two log brackets close
1089  // and the target is the cell above whatever u* is: the clamp only guards the arithmetic.
1090  const Real ustar = ib_stability::clamped_ustar(h_windspeed2r * kappa / (std::log((Real(1.5)) * dx_z / z0) - psi_m));
1091 
1092  // We do not know the actual temperature so use cell above
1093  const Real thetastar = theta * ustar * ustar / (kappa * ggg * Olen);
1094  const Real surf_temp = theta_neighbor - thetastar / kappa * (std::log((Real(1.5)) * dx_z / z0) - psi_h_neighbor);
1095  const Real tTarget = surf_temp + thetastar / kappa * (std::log((myhalf) * dx_z / z0) - psi_h);
1096 
1097  const Real bc_forcing_rt = -(cell_data(i,j,k,Rho_comp) * tTarget - cell_data(i,j,k,RhoTheta_comp));
1098  cell_src(i, j, k, RhoTheta_comp) -= drag_coefficient * U_s * bc_forcing_rt;
1099  }
1100  }
1101 
1102  // Force fully immersed cells to planar average rho and theta
1103  if (t_blank == one && r_avg && t_avg) {
1104  const Real rho_avg = r_avg(k);
1105  const Real theta_avg = t_avg(k) / rho_avg; // Convert from RhoTheta to Theta
1106  const Real rho_cell = cell_data(i,j,k,Rho_comp);
1107  const Real bc_forcing_r = -(rho_avg - rho_cell);
1108  const Real bc_forcing_rt = -(rho_avg * theta_avg - cell_data(i,j,k,RhoTheta_comp));
1109 
1110  cell_src(i, j, k, Rho_comp) -= drag_coefficient * U_s * bc_forcing_r;
1111  cell_src(i, j, k, RhoTheta_comp) -= drag_coefficient * U_s * bc_forcing_rt;
1112  }
1113  });
1114 }

Referenced by make_sources().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ ImmersedForcingTerrain_Xmom()

void ImmersedForcingTerrain_Xmom ( const Box &  tbx,
const Array4< const Real > &  u,
const Array4< const Real > &  v,
const Array4< const Real > &  w,
const Array4< const Real > &  cell_data,
const Array4< const Real > &  t_blank_arr,
const Array4< const Real > &  t_blank_xface_arr,
const Array4< const Real > &  z_cc_arr,
const Array4< Real > &  xmom_src_arr,
const Geometry &  geom,
const SolverChoice &  solverChoice,
const Real  fac 
)

Apply terrain immersed forcing to X-momentum

91 {
92  // geometric properties
93  const Real* dx_arr = geom.CellSize();
94  const Real dx_x = dx_arr[0];
95  const Real dx_y = dx_arr[1];
96  const Real dt = fac; // fac is actually dt in the calling code
97 
98  const Real alpha_m = solverChoice.if_Cd_momentum;
99  const Real tiny = std::numeric_limits<amrex::Real>::epsilon();
100  const Real U_s = one; // unit velocity scale
101  const bool l_implicit_drag = solverChoice.if_implicit_drag;
102 
103  // MOST parameters
104  similarity_funs sfuns;
105  const Real ggg = CONST_GRAV;
106  const Real kappa = KAPPA;
107  const Real z0 = solverChoice.if_z0;
108  const Real tflux_in = solverChoice.if_surf_temp_flux;
109  const Real Olen_in = solverChoice.if_Olen_in;
110  const Real stab_wind_floor = solverChoice.if_stability_wind_floor;
111  const Real psi_cap_factor = solverChoice.if_psi_cap_factor;
112  const bool l_use_most = solverChoice.if_use_most;
113 
114  const Real small_volfrac = 0.005;
115  // The terrain kernels keep the raw fractions: erf.if_snap_partial_cells
116  // applies to the buildings kernels only (the wall law below is weighted
117  // by the fluid fraction of the face, which a snapped face has not).
118 
119  ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
120  {
121  const Real ux = u(i, j, k);
122  const Real uy = fourth * ( v(i, j , k ) + v(i-1, j , k )
123  + v(i, j+1, k ) + v(i-1, j+1, k ) );
124  const Real uz = fourth * ( w(i, j , k ) + w(i-1, j , k )
125  + w(i, j , k+1) + w(i-1, j , k+1) );
126  const Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz);
127  // Use face-centered terrain_blanking if available, otherwise average from cell centers
128  Real t_blank_raw = (t_blank_xface_arr) ? t_blank_xface_arr(i, j, k) :
129  myhalf * (t_blank_arr(i, j, k) + t_blank_arr(i-1, j, k));
130  // Threshold: if averaged value is below small_volfrac, set to zero
131  const Real t_blank = (t_blank_raw < small_volfrac) ? zero : t_blank_raw;
132 
133  Real t_blank_above_raw = (t_blank_xface_arr) ? t_blank_xface_arr(i, j, k+1) :
134  myhalf * (t_blank_arr(i, j, k+1) + t_blank_arr(i-1, j, k+1));
135  const Real t_blank_above = (t_blank_above_raw < small_volfrac) ? zero : t_blank_above_raw;
136 
137  const Real dx_z = (z_cc_arr) ? (z_cc_arr(i,j,k) - z_cc_arr(i,j,k-1)) : dx_arr[2];
138  const Real drag_coefficient = alpha_m / std::pow(dx_x*dx_y*dx_z, one/three);
139  const Real CdM = std::min(drag_coefficient / (windspeed + tiny), drag_coefficient);
140 
141  const Real rho_xface = myhalf * ( cell_data(i,j,k,Rho_comp) + cell_data(i-1,j,k,Rho_comp) );
142 
143  if ((t_blank > 0 && (t_blank_above == zero)) && l_use_most) { // force to MOST value
144  // calculate tangential velocity one cell above
145  const Real ux2r = u(i, j, k+1) ;
146  const Real uy2r = fourth * ( v(i, j , k+1) + v(i-1, j , k+1)
147  + v(i, j+1, k+1) + v(i-1, j+1, k+1) ) ;
148  const Real h_windspeed2r = std::sqrt(ux2r * ux2r + uy2r * uy2r);
149 
150  // MOST
151  const Real theta_xface = (myhalf * (cell_data(i,j,k ,RhoTheta_comp) + cell_data(i-1,j,k, RhoTheta_comp))) / rho_xface;
152  const Real rho_xface_below = myhalf * ( cell_data(i,j,k-1,Rho_comp) + cell_data(i-1,j,k-1,Rho_comp) );
153  const Real theta_xface_below = (myhalf * (cell_data(i,j,k-1,RhoTheta_comp) + cell_data(i-1,j,k-1, RhoTheta_comp))) / rho_xface_below;
154  const Real theta_surf = theta_xface_below;
155 
156  Real psi_m = zero;
157  Real psi_h = zero;
158  // The stability estimate uses the wind floored at erf.if_stability_wind_floor, a derived
159  // Obukhov length bounded so that |zeta| <= ib_stability::zeta_max() as on flat ground,
160  // and a prescribed one as given with its zeta held there; the target uses the actual wind
161  Real ustar = ib_stability::floored_wind(h_windspeed2r, stab_wind_floor) * kappa / (std::log(Real(1.5) * dx_z / z0) - psi_m); // calculated from bottom of cell. Maintains flexibility for different Vf values
162  Real tflux = (tflux_in != Real(1e-8)) ? tflux_in : -(theta_xface - theta_surf) * ustar * kappa / (std::log(Real(1.5) * dx_z / z0) - psi_h);
163  Real Olen = (Olen_in != Real(1e-8)) ? Olen_in
164  : ib_stability::bounded_obukhov_length(-ustar * ustar * ustar * theta_xface / (kappa * ggg * tflux + tiny), Real(1.5) * dx_z);
165  Real zeta = ib_stability::bounded_zeta(Real(1.5) * dx_z, Olen);
166 
167  // similarity functions
168  psi_m = sfuns.calc_psi_m(zeta);
169  psi_h = sfuns.calc_psi_h(zeta);
170  ustar = h_windspeed2r * kappa / (std::log(Real(1.5) * dx_z / z0) - ib_stability::psi_m_for_ustar(psi_m, Real(1.5) * dx_z, z0, psi_cap_factor));
171 
172  // prevent some unphysical math
173  ustar = ib_stability::clamped_ustar(ustar);
174  psi_m = ib_stability::capped_psi(psi_m, myhalf * dx_z, z0, psi_cap_factor);
175 
176  // determine target velocity
177  const Real uTarget = ustar / kappa * (std::log(myhalf * dx_z / z0) - psi_m);
178  Real uxTarget = uTarget * ux2r / (tiny + h_windspeed2r);
179  const Real bc_forcing_x = -(uxTarget - ux); // BC forcing pushes nonrelative velocity toward target velocity
180  const Real lambda = (1-t_blank) * CdM * U_s; // affine relaxation rate toward MOST target [1/s]
181  const Real fac_local = l_implicit_drag ? lambda / (one + lambda*dt) : lambda; // point-implicit rescale (else explicit)
182  xmom_src_arr(i, j, k) -= fac_local * rho_xface * bc_forcing_x; // if Vf low, force more strongly to MOST. If high, less forcing.
183  } else {
184  const Real lambda = t_blank * CdM * windspeed; // linear drag rate [1/s]
185  const Real fac_local = l_implicit_drag ? lambda / (one + lambda*dt) : lambda; // point-implicit rescale (else explicit)
186  xmom_src_arr(i, j, k) -= fac_local * rho_xface * ux;
187  }
188  });
189 }

Referenced by make_mom_sources().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ ImmersedForcingTerrain_Ymom()

void ImmersedForcingTerrain_Ymom ( const Box &  tby,
const Array4< const Real > &  u,
const Array4< const Real > &  v,
const Array4< const Real > &  w,
const Array4< const Real > &  cell_data,
const Array4< const Real > &  t_blank_arr,
const Array4< const Real > &  t_blank_yface_arr,
const Array4< const Real > &  z_cc_arr,
const Array4< Real > &  ymom_src_arr,
const Geometry &  geom,
const SolverChoice &  solverChoice,
const Real  fac 
)

Apply terrain immersed forcing to Y-momentum

206 {
207  // geometric properties
208  const Real* dx_arr = geom.CellSize();
209  const Real dx_x = dx_arr[0];
210  const Real dx_y = dx_arr[1];
211  const Real dt = fac; // fac is actually dt in the calling code
212 
213  const Real alpha_m = solverChoice.if_Cd_momentum;
214  const Real tiny = std::numeric_limits<amrex::Real>::epsilon();
215  const Real U_s = one; // unit velocity scale
216  const bool l_implicit_drag = solverChoice.if_implicit_drag;
217 
218  // MOST parameters
219  similarity_funs sfuns;
220  const Real ggg = CONST_GRAV;
221  const Real kappa = KAPPA;
222  const Real z0 = solverChoice.if_z0;
223  const Real tflux_in = solverChoice.if_surf_temp_flux;
224  const Real Olen_in = solverChoice.if_Olen_in;
225  const Real stab_wind_floor = solverChoice.if_stability_wind_floor;
226  const Real psi_cap_factor = solverChoice.if_psi_cap_factor;
227  const bool l_use_most = solverChoice.if_use_most;
228 
229  const Real small_volfrac = 0.005;
230  // The terrain kernels keep the raw fractions: erf.if_snap_partial_cells
231  // applies to the buildings kernels only (the wall law below is weighted
232  // by the fluid fraction of the face, which a snapped face has not).
233 
234  ParallelFor(tby, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
235  {
236  const Real ux = fourth * ( u(i , j , k ) + u(i , j-1, k )
237  + u(i+1, j , k ) + u(i+1, j-1, k ) );
238  const Real uy = v(i, j, k);
239  const Real uz = fourth * ( w(i , j , k ) + w(i , j-1, k )
240  + w(i , j , k+1) + w(i , j-1, k+1) );
241  const Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz);
242  // Use face-centered terrain_blanking if available, otherwise average from cell centers
243  Real t_blank_raw = (t_blank_yface_arr) ? t_blank_yface_arr(i, j, k) :
244  myhalf * (t_blank_arr(i, j, k) + t_blank_arr(i, j-1, k));
245  const Real t_blank = (t_blank_raw < small_volfrac) ? zero : t_blank_raw;
246 
247  Real t_blank_above_raw = (t_blank_yface_arr) ? t_blank_yface_arr(i, j, k+1) :
248  myhalf * (t_blank_arr(i, j, k+1) + t_blank_arr(i, j-1, k+1));
249  const Real t_blank_above = (t_blank_above_raw < small_volfrac) ? zero : t_blank_above_raw;
250 
251  const Real dx_z = (z_cc_arr) ? (z_cc_arr(i,j,k) - z_cc_arr(i,j,k-1)) : dx_arr[2];
252  const Real drag_coefficient = alpha_m / std::pow(dx_x*dx_y*dx_z, one/three);
253  const Real CdM = std::min(drag_coefficient / (windspeed + tiny), drag_coefficient);
254 
255  const Real rho_yface = myhalf * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j-1,k,Rho_comp) );
256 
257  if ((t_blank > 0 && (t_blank_above == zero)) && l_use_most) { // force to MOST value
258  // calculate tangential velocity one cell above
259  const Real ux2r = fourth * ( u(i , j , k+1) + u(i , j-1, k+1)
260  + u(i+1, j , k+1) + u(i+1, j-1, k+1) );
261  const Real uy2r = v(i, j, k+1) ;
262  const Real h_windspeed2r = std::sqrt(ux2r * ux2r + uy2r * uy2r);
263 
264  // MOST
265  const Real theta_yface = (myhalf * (cell_data(i,j,k ,RhoTheta_comp) + cell_data(i,j-1,k, RhoTheta_comp))) / rho_yface;
266  const Real rho_yface_below = myhalf * ( cell_data(i,j,k-1,Rho_comp) + cell_data(i,j-1,k-1,Rho_comp) );
267  const Real theta_yface_below = (myhalf * (cell_data(i,j,k-1,RhoTheta_comp) + cell_data(i,j-1,k-1, RhoTheta_comp))) / rho_yface_below;
268  const Real theta_surf = theta_yface_below;
269 
270  Real psi_m = zero;
271  Real psi_h = zero;
272  // The stability estimate uses the wind floored at erf.if_stability_wind_floor, a derived
273  // Obukhov length bounded so that |zeta| <= ib_stability::zeta_max() as on flat ground,
274  // and a prescribed one as given with its zeta held there; the target uses the actual wind
275  Real ustar = ib_stability::floored_wind(h_windspeed2r, stab_wind_floor) * kappa / (std::log(Real(1.5) * dx_z / z0) - psi_m); // calculated from bottom of cell. Maintains flexibility for different Vf values
276  Real tflux = (tflux_in != Real(1e-8)) ? tflux_in : -(theta_yface - theta_surf) * ustar * kappa / (std::log(Real(1.5) * dx_z / z0) - psi_h);
277  Real Olen = (Olen_in != Real(1e-8)) ? Olen_in
278  : ib_stability::bounded_obukhov_length(-ustar * ustar * ustar * theta_yface / (kappa * ggg * tflux + tiny), Real(1.5) * dx_z);
279  Real zeta = ib_stability::bounded_zeta(Real(1.5) * dx_z, Olen);
280 
281  // similarity functions
282  psi_m = sfuns.calc_psi_m(zeta);
283  psi_h = sfuns.calc_psi_h(zeta);
284  ustar = h_windspeed2r * kappa / (std::log(Real(1.5) * dx_z / z0) - ib_stability::psi_m_for_ustar(psi_m, Real(1.5) * dx_z, z0, psi_cap_factor));
285 
286  // prevent some unphysical math
287  ustar = ib_stability::clamped_ustar(ustar);
288  psi_m = ib_stability::capped_psi(psi_m, myhalf * dx_z, z0, psi_cap_factor);
289 
290  // determine target velocity
291  const Real uTarget = ustar / kappa * (std::log(myhalf * dx_z / z0) - psi_m);
292  Real uyTarget = uTarget * uy2r / (tiny + h_windspeed2r);
293  const Real bc_forcing_y = -(uyTarget - uy); // BC forcing pushes nonrelative velocity toward target velocity
294  const Real lambda = (1 - t_blank) * CdM * U_s; // affine relaxation rate toward MOST target [1/s]
295  const Real fac_local = l_implicit_drag ? lambda / (one + lambda*dt) : lambda; // point-implicit rescale (else explicit)
296  ymom_src_arr(i, j, k) -= fac_local * rho_yface * bc_forcing_y; // if Vf low, force more strongly to MOST. If high, less forcing.
297  } else {
298  const Real lambda = t_blank * CdM * windspeed; // linear drag rate [1/s]
299  const Real fac_local = l_implicit_drag ? lambda / (one + lambda*dt) : lambda; // point-implicit rescale (else explicit)
300  ymom_src_arr(i, j, k) -= fac_local * rho_yface * uy;
301  }
302  });
303 }

Referenced by make_mom_sources().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ ImmersedForcingTerrain_Zmom()

void ImmersedForcingTerrain_Zmom ( const Box &  tbz,
const Array4< const Real > &  u,
const Array4< const Real > &  v,
const Array4< const Real > &  w,
const Array4< const Real > &  cell_data,
const Array4< const Real > &  t_blank_arr,
const Array4< const Real > &  t_blank_zface_arr,
const Array4< const Real > &  z_cc_arr,
const Array4< Real > &  zmom_src_arr,
const Geometry &  geom,
const SolverChoice &  solverChoice,
const Real  fac 
)

Apply terrain immersed forcing to Z-momentum

320 {
321  // geometric properties
322  const Real* dx_arr = geom.CellSize();
323  const Real dx_x = dx_arr[0];
324  const Real dx_y = dx_arr[1];
325  const Real dt = fac; // fac is actually dt in the calling code
326 
327  const Real alpha_m = solverChoice.if_Cd_momentum;
328  const Real tiny = std::numeric_limits<amrex::Real>::epsilon();
329  const bool l_implicit_drag = solverChoice.if_implicit_drag;
330 
331  const Real small_volfrac = 0.005;
332  // The terrain kernels keep the raw fractions: erf.if_snap_partial_cells
333  // applies to the buildings kernels only (the wall law below is weighted
334  // by the fluid fraction of the face, which a snapped face has not).
335 
336  ParallelFor(tbz, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
337  {
338  const Real ux = fourth * ( u(i , j , k ) + u(i+1, j , k )
339  + u(i , j , k-1) + u(i+1, j , k-1) );
340  const Real uy = fourth * ( v(i , j , k ) + v(i , j+1, k )
341  + v(i , j , k-1) + v(i , j+1, k-1) );
342  const Real uz = w(i, j, k);
343  const Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz);
344  // Use face-centered terrain_blanking if available, otherwise average from cell centers
345  Real t_blank_raw = (t_blank_zface_arr) ? t_blank_zface_arr(i, j, k) :
346  myhalf * (t_blank_arr(i, j, k) + t_blank_arr(i, j, k-1));
347  const Real t_blank = (t_blank_raw < small_volfrac) ? zero : t_blank_raw;
348 
349  const Real dx_z = (z_cc_arr) ? (z_cc_arr(i,j,k) - z_cc_arr(i,j,k-1)) : dx_arr[2];
350  const Real drag_coefficient = alpha_m / std::pow(dx_x*dx_y*dx_z, one/three);
351  const Real CdM = std::min(drag_coefficient / (windspeed + tiny), drag_coefficient);
352 
353  const Real rho_zface = myhalf * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j,k-1,Rho_comp) );
354  const Real lambda = t_blank * CdM * windspeed; // linear drag rate [1/s]
355  const Real fac_local = l_implicit_drag ? lambda / (one + lambda*dt) : lambda; // point-implicit rescale (else explicit)
356  zmom_src_arr(i, j, k) -= fac_local * rho_zface * uz;
357  });
358 }

Referenced by make_mom_sources().

Here is the call graph for this function:
Here is the caller graph for this function: