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 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)
 
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 
)

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.
Returns
Target velocity component.
38 {
39  const Real tiny = std::numeric_limits<amrex::Real>::epsilon();
40  Real psi_m = zero;
41  Real psi_h = zero;
42  Real tang_windspeed2r = std::sqrt(u1_2r * u1_2r + u2_2r * u2_2r);
43 
44  Real ustar = tang_windspeed2r * KAPPA / (std::log(1.5 * delta / z0) - psi_m);
45  Real tflux = (tflux_in != Real(1.e-8)) ? tflux_in : -(theta_face - theta_surf) * ustar * KAPPA / (std::log(1.5 * delta / z0) - psi_h);
46  Real Olen = (Olen_in != Real(1.e-8)) ? Olen_in : -ustar * ustar * ustar * theta_face / (KAPPA * CONST_GRAV * tflux + tiny);
47  Real zeta = 1.5 * delta / Olen;
48 
49  // similarity functions
50  similarity_funs sfuns;
51  if (stability_correction){
52  psi_m = sfuns.calc_psi_m(zeta);
53  psi_h = sfuns.calc_psi_h(zeta);
54  }
55  ustar = tang_windspeed2r * KAPPA / (std::log(1.5 * delta / z0) - psi_m);
56 
57  // prevent some unphysical math
58  if (!(ustar > zero && !std::isnan(ustar))) { ustar = zero; }
59  if (!(ustar < 2.0 && !std::isnan(ustar))) { ustar = 2.0; }
60  if (psi_m > std::log(myhalf * delta / z0)) { psi_m = std::log(myhalf * delta / z0); }
61 
62  Real uTarget = (1 - t_blank) * ustar / KAPPA * (std::log(myhalf * delta / z0) - psi_m);
63  Real u1Target = uTarget * u1_2r / (tiny + tang_windspeed2r);
64 
65  return u1Target;
66 }
KAPPA
constexpr amrex::Real KAPPA
Definition: ERF_Constants.H:55
CONST_GRAV
constexpr amrex::Real CONST_GRAV
Definition: ERF_Constants.H:56
z0
Real z0
Definition: ERF_InitCustomPertVels_ScalarAdvDiff.H:8
zero
constexpr amrex::Real zero
Definition: ERF_NumericalConstants.H:29
myhalf
constexpr amrex::Real myhalf
Definition: ERF_NumericalConstants.H:34
Real
amrex::Real Real
Definition: ERF_ShocInterface.H:19
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)

1097 {
1098  // geometric properties
1099  const Real* dx_arr = geom.CellSize();
1100  const Real dx_x = dx_arr[0];
1101  const Real dx_y = dx_arr[1];
1102 
1103  const Real alpha_h = solverChoice.if_Cd_scalar;
1104  // erf.if_snap_partial_cells: read the cell blanking snapped to solid (1)
1105  // or fluid (0) at half, so a height-map building becomes the same
1106  // staircase of whole cells an exact box is: the thermal conditions sit
1107  // on the boundary solid cells, roofs included, found from the neighbour
1108  // blanking. Off (the default), the raw fractions are used and nothing
1109  // below changes. See SolverChoice::if_snap_partial_cells.
1110  const bool l_snap = solverChoice.if_snap_partial_cells;
1111  auto snapb = [=] AMREX_GPU_DEVICE (amrex::Real b) noexcept -> amrex::Real {
1112  return l_snap ? ((b >= myhalf) ? one : zero) : b;
1113  };
1114  const Real tiny = std::numeric_limits<amrex::Real>::epsilon();
1115  const Real U_s = one; // unit velocity scale
1116 
1117  // MOST parameters
1118  similarity_funs sfuns;
1119  const Real ggg = CONST_GRAV;
1120  const Real kappa = KAPPA;
1121  const Real z0 = solverChoice.if_z0;
1122  const Real tflux = solverChoice.if_surf_temp_flux;
1123  const Real init_surf_temp = solverChoice.if_init_surf_temp;
1124  const Real surf_heating_rate = solverChoice.if_surf_heating_rate;
1125  const Real Olen_in = solverChoice.if_Olen_in;
1126 
1127  ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
1128  {
1129  const Real t_blank = snapb(t_blank_arr(i, j, k));
1130  const Real t_blank_below = snapb(t_blank_arr(i, j, k-1));
1131  const Real t_blank_above = snapb(t_blank_arr(i, j, k+1));
1132  const Real t_blank_north = snapb(t_blank_arr(i , j+1, k));
1133  const Real t_blank_south = snapb(t_blank_arr(i , j-1, k));
1134  const Real t_blank_east = snapb(t_blank_arr(i+1, j , k));
1135  const Real t_blank_west = snapb(t_blank_arr(i-1, j , k));
1136 
1137  const Real dx_z = (z_cc_arr) ? (z_cc_arr(i,j,k) - z_cc_arr(i,j,k-1)) : dx_arr[2];
1138  Real drag_coefficient = alpha_h / std::pow(dx_x*dx_y*dx_z, one/three);
1139 
1140  // Wall cells of the building, named by the face they carry: a
1141  // partially blanked cell with a more solid neighbour behind it and a
1142  // fluid one in front, or, with the snap (every solid cell at 1), a
1143  // solid cell with a solid neighbour behind and a fluid one in front.
1144  const bool south_face = l_snap ? (t_blank == one && t_blank_north == one && t_blank_south == zero)
1145  : (t_blank > zero && t_blank < t_blank_north && t_blank_south == zero);
1146  const bool north_face = l_snap ? (t_blank == one && t_blank_south == one && t_blank_north == zero)
1147  : (t_blank > zero && t_blank < t_blank_south && t_blank_north == zero);
1148  const bool west_face = l_snap ? (t_blank == one && t_blank_east == one && t_blank_west == zero)
1149  : (t_blank > zero && t_blank < t_blank_east && t_blank_west == zero);
1150  const bool east_face = l_snap ? (t_blank == one && t_blank_west == one && t_blank_east == zero)
1151  : (t_blank > zero && t_blank < t_blank_west && t_blank_east == zero);
1152 
1153  // SURFACE TEMP AND HEATING/COOLING RATE
1154  if (init_surf_temp > zero) {
1155  const Real surf_temp = init_surf_temp + surf_heating_rate*time;
1156  if (t_blank > 0 && (t_blank_above == zero) && (t_blank_below == one)) { // building roof
1157  const Real bc_forcing_rt_srf = -(cell_data(i,j,k,Rho_comp) * surf_temp - cell_data(i,j,k,RhoTheta_comp));
1158  cell_src(i, j, k, RhoTheta_comp) -= drag_coefficient * U_s * bc_forcing_rt_srf;
1159 
1160  } else if (east_face || west_face || south_face || north_face) {
1161  // this should enter for just building walls
1162  // walls are currently separated to allow for flexibility in the future to heat walls differently
1163 
1164  // south face
1165  if (l_snap ? south_face : ((t_blank < t_blank_north) && (t_blank_north == one))) {
1166  const Real bc_forcing_rt_srf = -(cell_data(i,j,k,Rho_comp) * surf_temp - cell_data(i,j,k,RhoTheta_comp));
1167  cell_src(i, j, k, RhoTheta_comp) -= drag_coefficient * U_s * bc_forcing_rt_srf;
1168  }
1169 
1170  // north face
1171  if (l_snap ? north_face : ((t_blank < t_blank_south) && (t_blank_south == one))) {
1172  const Real bc_forcing_rt_srf = -(cell_data(i,j,k,Rho_comp) * surf_temp - cell_data(i,j,k,RhoTheta_comp));
1173  cell_src(i, j, k, RhoTheta_comp) -= drag_coefficient * U_s * bc_forcing_rt_srf;
1174  }
1175 
1176  // west face
1177  if (l_snap ? west_face : ((t_blank < t_blank_east) && (t_blank_east == one))) {
1178  const Real bc_forcing_rt_srf = -(cell_data(i,j,k,Rho_comp) * surf_temp - cell_data(i,j,k,RhoTheta_comp));
1179  cell_src(i, j, k, RhoTheta_comp) -= drag_coefficient * U_s * bc_forcing_rt_srf;
1180  }
1181 
1182  // east face
1183  if (l_snap ? east_face : ((t_blank < t_blank_west) && (t_blank_west == one))) {
1184  const Real bc_forcing_rt_srf = -(cell_data(i,j,k,Rho_comp) * surf_temp - cell_data(i,j,k,RhoTheta_comp));
1185  cell_src(i, j, k, RhoTheta_comp) -= drag_coefficient * U_s * bc_forcing_rt_srf;
1186  }
1187 
1188  }
1189  }
1190 
1191  // SURFACE HEAT FLUX
1192  if (tflux != Real(1.e-8)){
1193  const Real ux_cc_2r = myhalf * (u(i ,j ,k+1) + u(i+1,j ,k+1));
1194  const Real uy_cc_2r = myhalf * (v(i ,j ,k+1) + v(i ,j+1,k+1));
1195  const Real h_windspeed2r = std::sqrt(ux_cc_2r * ux_cc_2r + uy_cc_2r * uy_cc_2r);
1196 
1197  const Real theta = cell_data(i,j,k ,RhoTheta_comp) / cell_data(i,j,k ,Rho_comp);
1198  Real theta_neighbor = cell_data(i,j,k+1,RhoTheta_comp) / cell_data(i,j,k+1,Rho_comp);
1199 
1200  if (t_blank > zero && (t_blank_above == zero)) { // building roof
1201  Real psi_m = zero;
1202  Real psi_h = zero;
1203  Real psi_h_neighbor = zero;
1204  Real ustar = h_windspeed2r * kappa / (std::log((1.5) * dx_z / z0) - psi_m);
1205  Real Olen = (Olen_in != Real(1e-8)) ? Olen_in : -ustar * ustar * ustar * theta / (kappa * ggg * tflux + tiny);
1206 
1207  for (int iter = 0; iter < 2; ++iter) {
1208  if (iter > 0) { Olen = -ustar * ustar * ustar * theta / (kappa * ggg * tflux + tiny); }
1209  Real zeta = (myhalf) * dx_z / Olen;
1210  Real zeta_neighbor = (1.5) * dx_z / Olen;
1211 
1212  // similarity functions
1213  psi_m = sfuns.calc_psi_m(zeta);
1214  psi_h = sfuns.calc_psi_h(zeta);
1215  psi_h_neighbor = sfuns.calc_psi_h(zeta_neighbor);
1216  ustar = h_windspeed2r * kappa / (std::log((1.5) * dx_z / z0) - psi_m);
1217  }
1218 
1219  // prevent some unphysical math
1220  if (!(ustar > zero && !std::isnan(ustar))) { ustar = zero; }
1221  if (!(ustar < 2.0 && !std::isnan(ustar))) { ustar = 2.0; }
1222  if (psi_h_neighbor > std::log(1.5 * dx_z / z0)) { psi_h_neighbor = std::log(1.5 * dx_z / z0); }
1223  if (psi_h > std::log(myhalf * dx_z / z0)) { psi_h = std::log(myhalf * dx_z / z0); }
1224 
1225  // We do not know the actual temperature so use cell above
1226  const Real thetastar = theta * ustar * ustar / (kappa * ggg * Olen);
1227  const Real surf_temp = theta_neighbor - thetastar / kappa * (std::log((1.5) * dx_z / z0) - psi_h_neighbor);
1228  const Real tTarget = surf_temp + thetastar / kappa * (std::log((myhalf) * dx_z / z0) - psi_h);
1229 
1230  const Real bc_forcing_rt = -(cell_data(i,j,k,Rho_comp) * tTarget - cell_data(i,j,k,RhoTheta_comp));
1231  cell_src(i, j, k, RhoTheta_comp) -= drag_coefficient * U_s * bc_forcing_rt;
1232 
1233  } else if (east_face || west_face || south_face || north_face) { // this should enter for just building walls
1234 
1235  Real ux_cellaway = zero;
1236  Real uy_cellaway = zero;
1237  Real uz_cellaway = zero;
1238  Real u1 = zero;
1239  Real u2 = zero;
1240  Real delta = zero;
1241 
1242  // south face
1243  if (south_face) {
1244  ux_cellaway = myhalf * (u(i ,j-1,k) + u(i+1,j-1,k ));
1245  uz_cellaway = myhalf * (w(i ,j-1,k) + w(i ,j-1,k+1));
1246  u1 = ux_cellaway;
1247  u2 = uz_cellaway;
1248  delta = dx_y;
1249 
1250  // MOST
1251  theta_neighbor = cell_data(i,j-1,k,RhoTheta_comp) / cell_data(i,j-1,k,Rho_comp);
1252  }
1253 
1254  // north face
1255  if (north_face) {
1256  ux_cellaway = myhalf * (u(i ,j+1,k) + u(i+1,j+1,k ));
1257  uz_cellaway = myhalf * (w(i ,j+1,k) + w(i ,j+1,k+1));
1258  u1 = ux_cellaway;
1259  u2 = uz_cellaway;
1260  delta = dx_y;
1261 
1262  // MOST
1263  theta_neighbor = cell_data(i,j+1,k,RhoTheta_comp) / cell_data(i,j+1,k,Rho_comp);
1264  }
1265 
1266  // west face
1267  if (west_face) {
1268  uy_cellaway = myhalf * (v(i-1,j ,k) + v(i-1,j+1,k ));
1269  uz_cellaway = myhalf * (w(i-1,j ,k) + w(i-1,j ,k+1));
1270  u1 = uy_cellaway;
1271  u2 = uz_cellaway;
1272  delta = dx_x;
1273 
1274  // MOST
1275  theta_neighbor = cell_data(i-1,j,k,RhoTheta_comp) / cell_data(i-1,j,k,Rho_comp);
1276  }
1277 
1278  // east face
1279  if (east_face) {
1280  uy_cellaway = myhalf * (v(i+1,j ,k) + v(i+1,j+1,k ));
1281  uz_cellaway = myhalf * (w(i+1,j ,k) + w(i+1,j ,k+1));
1282  u1 = uy_cellaway;
1283  u2 = uz_cellaway;
1284  delta = dx_x;
1285 
1286  // MOST
1287  theta_neighbor = cell_data(i+1,j,k,RhoTheta_comp) / cell_data(i+1,j,k,Rho_comp);
1288  }
1289 
1290  Real tan_wspd = std::sqrt(u1 * u1 + u2 * u2);
1291 
1292  Real psi_m = zero;
1293  Real psi_h = zero;
1294  Real psi_h_neighbor = zero;
1295  Real ustar = tan_wspd * kappa / (std::log(1.5 * delta / z0) - psi_m);
1296  Real Olen = (Olen_in != Real(1e-8)) ? Olen_in : -ustar * ustar * ustar * theta / (kappa * ggg * tflux + tiny);
1297 
1298  for (int iter = 0; iter < 2; ++iter) {
1299  if (iter > 0) { Olen = -ustar * ustar * ustar * theta / (kappa * ggg * tflux + tiny); }
1300  Real zeta = (myhalf) * delta / Olen;
1301  Real zeta_neighbor = (1.5) * delta / Olen;
1302 
1303  // similarity functions
1304  psi_m = sfuns.calc_psi_m(zeta);
1305  psi_h = sfuns.calc_psi_h(zeta);
1306  psi_h_neighbor = sfuns.calc_psi_h(zeta_neighbor);
1307  ustar = tan_wspd * kappa / (std::log((1.5) * delta / z0) - psi_m);
1308  }
1309 
1310  // prevent some unphysical math
1311  if (!(ustar > zero && !std::isnan(ustar))) { ustar = zero; }
1312  if (!(ustar < 2.0 && !std::isnan(ustar))) { ustar = 2.0; }
1313  if (psi_h_neighbor > std::log(1.5 * delta / z0)) { psi_h_neighbor = std::log(1.5 * delta / z0); }
1314  if (psi_h > std::log(myhalf * delta / z0)) { psi_h = std::log(myhalf * delta / z0); }
1315 
1316  // We do not know the actual temperature so use cell above
1317  const Real thetastar = theta * ustar * ustar / (kappa * ggg * Olen);
1318  const Real surf_temp = theta_neighbor - thetastar / kappa * (std::log((1.5) * delta / z0) - psi_h_neighbor);
1319  const Real tTarget = surf_temp + thetastar / kappa * (std::log((myhalf) * delta / z0) - psi_h);
1320 
1321  const Real bc_forcing_rt = -(cell_data(i,j,k,Rho_comp) * tTarget - cell_data(i,j,k,RhoTheta_comp));
1322  cell_src(i, j, k, RhoTheta_comp) -= drag_coefficient * U_s * bc_forcing_rt;
1323  }
1324  }
1325 
1326  // Force fully immersed cells to planar average rho and theta
1327  if (t_blank == 1.0 && r_avg && t_avg) {
1328  const Real rho_avg = r_avg(k);
1329  const Real theta_avg = t_avg(k) / rho_avg; // Convert from RhoTheta to Theta
1330  const Real rho_cell = cell_data(i,j,k,Rho_comp);
1331  const Real bc_forcing_r = -(rho_avg - rho_cell);
1332  const Real bc_forcing_rt = -(rho_avg * theta_avg - cell_data(i,j,k,RhoTheta_comp));
1333 
1334  cell_src(i, j, k, Rho_comp) -= drag_coefficient * U_s * bc_forcing_r;
1335  cell_src(i, j, k, RhoTheta_comp) -= drag_coefficient * U_s * bc_forcing_rt;
1336  }
1337  });
1338 }
Rho_comp
#define Rho_comp
Definition: ERF_IndexDefines.H:39
RhoTheta_comp
#define RhoTheta_comp
Definition: ERF_IndexDefines.H:40
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:32
one
constexpr amrex::Real one
Definition: ERF_NumericalConstants.H:30
w
Real w
Definition: ERF_Plotfile2DInterpolator.cpp:22
LsmVar_SLM::theta
@ theta
Definition: ERF_SLM.H:19
SolverChoice::if_Olen_in
amrex::Real if_Olen_in
Input Obukhov length for immersed-forcing MOST [m].
Definition: ERF_DataStruct.H:2046
SolverChoice::if_z0
amrex::Real if_z0
Immersed-forcing roughness length [m].
Definition: ERF_DataStruct.H:2042
SolverChoice::if_snap_partial_cells
bool if_snap_partial_cells
Definition: ERF_DataStruct.H:2048
SolverChoice::if_Cd_scalar
amrex::Real if_Cd_scalar
Immersed-forcing drag coefficient for scalars.
Definition: ERF_DataStruct.H:2037
SolverChoice::if_init_surf_temp
amrex::Real if_init_surf_temp
Initial immersed-forcing surface temperature [K].
Definition: ERF_DataStruct.H:2044
SolverChoice::if_surf_temp_flux
amrex::Real if_surf_temp_flux
Immersed-forcing surface temperature flux [K m/s].
Definition: ERF_DataStruct.H:2043
SolverChoice::if_surf_heating_rate
amrex::Real if_surf_heating_rate
Immersed-forcing surface heating rate [K/hr].
Definition: ERF_DataStruct.H:2045

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

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

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

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

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

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

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)

961 {
962  // geometric properties
963  const Real* dx_arr = geom.CellSize();
964  const Real dx_x = dx_arr[0];
965  const Real dx_y = dx_arr[1];
966 
967  const Real alpha_h = solverChoice.if_Cd_scalar;
968  // The terrain kernels keep the raw fractions: erf.if_snap_partial_cells
969  // applies to the buildings kernels only (the wall law below is weighted
970  // by the fluid fraction of the face, which a snapped face has not).
971  const Real tiny = std::numeric_limits<amrex::Real>::epsilon();
972  const Real U_s = one; // unit velocity scale
973 
974  // MOST parameters
975  similarity_funs sfuns;
976  const Real ggg = CONST_GRAV;
977  const Real kappa = KAPPA;
978  const Real z0 = solverChoice.if_z0;
979  const Real tflux = solverChoice.if_surf_temp_flux;
980  const Real init_surf_temp = solverChoice.if_init_surf_temp;
981 
982  // Note this has been converted to K / s when it was read in;
983  const Real surf_heating_rate = solverChoice.if_surf_heating_rate;
984 
985  const Real Olen_in = solverChoice.if_Olen_in;
986 
987  ParallelFor(bx, [=]
988  AMREX_GPU_DEVICE(int i, int j, int k) noexcept
989  {
990  const Real dx_z = (z_cc_arr) ? (z_cc_arr(i,j,k) - z_cc_arr(i,j,k-1)) : dx_arr[2];
991  const Real drag_coefficient = alpha_h / std::pow(dx_x*dx_y*dx_z, one/three);
992 
993  const Real t_blank = t_blank_arr(i, j, k);
994  const Real t_blank_above = t_blank_arr(i, j, k+1);
995  const Real ux_cc_2r = myhalf * (u(i ,j ,k+1) + u(i+1,j ,k+1));
996  const Real uy_cc_2r = myhalf * (v(i ,j ,k+1) + v(i ,j+1,k+1));
997  const Real h_windspeed2r = std::sqrt(ux_cc_2r * ux_cc_2r + uy_cc_2r * uy_cc_2r);
998 
999  const Real theta = cell_data(i,j,k ,RhoTheta_comp) / cell_data(i,j,k ,Rho_comp);
1000  const Real theta_neighbor = cell_data(i,j,k+1,RhoTheta_comp) / cell_data(i,j,k+1,Rho_comp);
1001 
1002  // SURFACE TEMP AND HEATING/COOLING RATE
1003  if (init_surf_temp > zero) {
1004  if (t_blank > 0 && (t_blank_above == zero)) { // force to MOST value
1005  const Real surf_temp = init_surf_temp + surf_heating_rate*time;
1006  const Real bc_forcing_rt_srf = -(cell_data(i,j,k-1,Rho_comp) * surf_temp - cell_data(i,j,k-1,RhoTheta_comp));
1007  cell_src(i, j, k, RhoTheta_comp) -= drag_coefficient * U_s * bc_forcing_rt_srf;
1008  }
1009  }
1010 
1011  // SURFACE HEAT FLUX
1012  if (tflux != Real(1e-8)){
1013  if (t_blank > 0 && (t_blank_above == zero)) { // force to MOST value
1014  Real psi_m = zero;
1015  Real psi_h = zero;
1016  Real psi_h_neighbor = zero;
1017  Real ustar = h_windspeed2r * kappa / (std::log((Real(1.5)) * dx_z / z0) - psi_m);
1018  const Real Olen = -ustar * ustar * ustar * theta / (kappa * ggg * tflux + tiny);
1019  const Real zeta = (myhalf) * dx_z / Olen;
1020  const Real zeta_neighbor = (Real(1.5)) * dx_z / Olen;
1021 
1022  // similarity functions
1023  psi_m = sfuns.calc_psi_m(zeta);
1024  psi_h = sfuns.calc_psi_h(zeta);
1025  psi_h_neighbor = sfuns.calc_psi_h(zeta_neighbor);
1026  ustar = h_windspeed2r * kappa / (std::log((Real(1.5)) * dx_z / z0) - psi_m);
1027 
1028  // prevent some unphysical math
1029  if (!(ustar > zero && !std::isnan(ustar))) { ustar = zero; }
1030  if (!(ustar < two && !std::isnan(ustar))) { ustar = two; }
1031  if (psi_h_neighbor > std::log(Real(1.5) * dx_z / z0)) { psi_h_neighbor = std::log(Real(1.5) * dx_z / z0); }
1032  if (psi_h > std::log(myhalf * dx_z / z0)) { psi_h = std::log(myhalf * dx_z / z0); }
1033 
1034  // We do not know the actual temperature so use cell above
1035  const Real thetastar = theta * ustar * ustar / (kappa * ggg * Olen);
1036  const Real surf_temp = theta_neighbor - thetastar / kappa * (std::log((Real(1.5)) * dx_z / z0) - psi_h_neighbor);
1037  const Real tTarget = surf_temp + thetastar / kappa * (std::log((myhalf) * dx_z / z0) - psi_h);
1038 
1039  const Real bc_forcing_rt = -(cell_data(i,j,k,Rho_comp) * tTarget - cell_data(i,j,k,RhoTheta_comp));
1040  cell_src(i, j, k, RhoTheta_comp) -= drag_coefficient * U_s * bc_forcing_rt;
1041  }
1042  }
1043 
1044  // OBUKHOV LENGTH
1045  if (Olen_in != Real(1e-8)){
1046  if (t_blank > 0 && (t_blank_above == zero)) { // force to MOST value
1047  const Real Olen = Olen_in;
1048  const Real zeta = (myhalf) * dx_z / Olen;
1049  const Real zeta_neighbor = (Real(1.5)) * dx_z / Olen;
1050 
1051  // similarity functions
1052  const Real psi_m = sfuns.calc_psi_m(zeta);
1053  const Real psi_h = sfuns.calc_psi_h(zeta);
1054  const Real psi_h_neighbor = sfuns.calc_psi_h(zeta_neighbor);
1055  const Real ustar = h_windspeed2r * kappa / (std::log((Real(1.5)) * dx_z / z0) - psi_m);
1056 
1057  // We do not know the actual temperature so use cell above
1058  const Real thetastar = theta * ustar * ustar / (kappa * ggg * Olen);
1059  const Real surf_temp = theta_neighbor - thetastar / kappa * (std::log((Real(1.5)) * dx_z / z0) - psi_h_neighbor);
1060  const Real tTarget = surf_temp + thetastar / kappa * (std::log((myhalf) * dx_z / z0) - psi_h);
1061 
1062  const Real bc_forcing_rt = -(cell_data(i,j,k,Rho_comp) * tTarget - cell_data(i,j,k,RhoTheta_comp));
1063  cell_src(i, j, k, RhoTheta_comp) -= drag_coefficient * U_s * bc_forcing_rt;
1064  }
1065  }
1066 
1067  // Force fully immersed cells to planar average rho and theta
1068  if (t_blank == one && r_avg && t_avg) {
1069  const Real rho_avg = r_avg(k);
1070  const Real theta_avg = t_avg(k) / rho_avg; // Convert from RhoTheta to Theta
1071  const Real rho_cell = cell_data(i,j,k,Rho_comp);
1072  const Real bc_forcing_r = -(rho_avg - rho_cell);
1073  const Real bc_forcing_rt = -(rho_avg * theta_avg - cell_data(i,j,k,RhoTheta_comp));
1074 
1075  cell_src(i, j, k, Rho_comp) -= drag_coefficient * U_s * bc_forcing_r;
1076  cell_src(i, j, k, RhoTheta_comp) -= drag_coefficient * U_s * bc_forcing_rt;
1077  }
1078  });
1079 }
two
constexpr amrex::Real two
Definition: ERF_NumericalConstants.H:31

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

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

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

193 {
194  // geometric properties
195  const Real* dx_arr = geom.CellSize();
196  const Real dx_x = dx_arr[0];
197  const Real dx_y = dx_arr[1];
198  const Real dt = fac; // fac is actually dt in the calling code
199 
200  const Real alpha_m = solverChoice.if_Cd_momentum;
201  const Real tiny = std::numeric_limits<amrex::Real>::epsilon();
202  const Real U_s = one; // unit velocity scale
203  const bool l_implicit_drag = solverChoice.if_implicit_drag;
204 
205  // MOST parameters
206  similarity_funs sfuns;
207  const Real ggg = CONST_GRAV;
208  const Real kappa = KAPPA;
209  const Real z0 = solverChoice.if_z0;
210  const Real tflux_in = solverChoice.if_surf_temp_flux;
211  const Real Olen_in = solverChoice.if_Olen_in;
212  const bool l_use_most = solverChoice.if_use_most;
213 
214  const Real small_volfrac = 0.005;
215  // The terrain kernels keep the raw fractions: erf.if_snap_partial_cells
216  // applies to the buildings kernels only (the wall law below is weighted
217  // by the fluid fraction of the face, which a snapped face has not).
218 
219  ParallelFor(tby, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
220  {
221  const Real ux = fourth * ( u(i , j , k ) + u(i , j-1, k )
222  + u(i+1, j , k ) + u(i+1, j-1, k ) );
223  const Real uy = v(i, j, k);
224  const Real uz = fourth * ( w(i , j , k ) + w(i , j-1, k )
225  + w(i , j , k+1) + w(i , j-1, k+1) );
226  const Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz);
227  // Use face-centered terrain_blanking if available, otherwise average from cell centers
228  Real t_blank_raw = (t_blank_yface_arr) ? t_blank_yface_arr(i, j, k) :
229  myhalf * (t_blank_arr(i, j, k) + t_blank_arr(i, j-1, k));
230  const Real t_blank = (t_blank_raw < small_volfrac) ? zero : t_blank_raw;
231 
232  Real t_blank_above_raw = (t_blank_yface_arr) ? t_blank_yface_arr(i, j, k+1) :
233  myhalf * (t_blank_arr(i, j, k+1) + t_blank_arr(i, j-1, k+1));
234  const Real t_blank_above = (t_blank_above_raw < small_volfrac) ? zero : t_blank_above_raw;
235 
236  const Real dx_z = (z_cc_arr) ? (z_cc_arr(i,j,k) - z_cc_arr(i,j,k-1)) : dx_arr[2];
237  const Real drag_coefficient = alpha_m / std::pow(dx_x*dx_y*dx_z, one/three);
238  const Real CdM = std::min(drag_coefficient / (windspeed + tiny), drag_coefficient);
239 
240  const Real rho_yface = myhalf * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j-1,k,Rho_comp) );
241 
242  if ((t_blank > 0 && (t_blank_above == zero)) && l_use_most) { // force to MOST value
243  // calculate tangential velocity one cell above
244  const Real ux2r = fourth * ( u(i , j , k+1) + u(i , j-1, k+1)
245  + u(i+1, j , k+1) + u(i+1, j-1, k+1) );
246  const Real uy2r = v(i, j, k+1) ;
247  const Real h_windspeed2r = std::sqrt(ux2r * ux2r + uy2r * uy2r);
248 
249  // MOST
250  const Real theta_yface = (myhalf * (cell_data(i,j,k ,RhoTheta_comp) + cell_data(i,j-1,k, RhoTheta_comp))) / rho_yface;
251  const Real rho_yface_below = myhalf * ( cell_data(i,j,k-1,Rho_comp) + cell_data(i,j-1,k-1,Rho_comp) );
252  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;
253  const Real theta_surf = theta_yface_below;
254 
255  Real psi_m = zero;
256  Real psi_h = zero;
257  Real ustar = h_windspeed2r * kappa / (std::log(Real(1.5) * dx_z / z0) - psi_m); // calculated from bottom of cell. Maintains flexibility for different Vf values
258  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);
259  Real Olen = (Olen_in != Real(1e-8)) ? Olen_in : -ustar * ustar * ustar * theta_yface / (kappa * ggg * tflux + tiny);
260  Real zeta = Real(1.5) * dx_z / Olen;
261 
262  // similarity functions
263  psi_m = sfuns.calc_psi_m(zeta);
264  psi_h = sfuns.calc_psi_h(zeta);
265  ustar = h_windspeed2r * kappa / (std::log(Real(1.5) * dx_z / z0) - psi_m);
266 
267  // prevent some unphysical math
268  if (!(ustar > zero && !std::isnan(ustar))) { ustar = zero; }
269  if (!(ustar < two && !std::isnan(ustar))) { ustar = two; }
270  if (psi_m > std::log(myhalf * dx_z / z0)) { psi_m = std::log(myhalf * dx_z / z0); }
271 
272  // determine target velocity
273  const Real uTarget = ustar / kappa * (std::log(myhalf * dx_z / z0) - psi_m);
274  Real uyTarget = uTarget * uy2r / (tiny + h_windspeed2r);
275  const Real bc_forcing_y = -(uyTarget - uy); // BC forcing pushes nonrelative velocity toward target velocity
276  const Real lambda = (1 - t_blank) * CdM * U_s; // affine relaxation rate toward MOST target [1/s]
277  const Real fac_local = l_implicit_drag ? lambda / (one + lambda*dt) : lambda; // point-implicit rescale (else explicit)
278  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.
279  } else {
280  const Real lambda = t_blank * CdM * windspeed; // linear drag rate [1/s]
281  const Real fac_local = l_implicit_drag ? lambda / (one + lambda*dt) : lambda; // point-implicit rescale (else explicit)
282  ymom_src_arr(i, j, k) -= fac_local * rho_yface * uy;
283  }
284  });
285 }

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

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

Referenced by make_mom_sources().

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