ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
Functions
ERF_ImmersedForcing.cpp File Reference
#include "ERF_ImmersedForcing.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.
37 {
38  const Real tiny = std::numeric_limits<amrex::Real>::epsilon();
39  Real psi_m = zero;
40  Real psi_h = zero;
41  Real tang_windspeed2r = std::sqrt(u1_2r * u1_2r + u2_2r * u2_2r);
42 
43  Real ustar = tang_windspeed2r * KAPPA / (std::log(1.5 * delta / z0) - psi_m);
44  Real tflux = (tflux_in != Real(1.e-8)) ? tflux_in : -(theta_face - theta_surf) * ustar * KAPPA / (std::log(1.5 * delta / z0) - psi_h);
45  Real Olen = (Olen_in != Real(1.e-8)) ? Olen_in : -ustar * ustar * ustar * theta_face / (KAPPA * CONST_GRAV * tflux + tiny);
46  Real zeta = 1.5 * delta / Olen;
47 
48  // similarity functions
49  similarity_funs sfuns;
50  if (stability_correction){
51  psi_m = sfuns.calc_psi_m(zeta);
52  psi_h = sfuns.calc_psi_h(zeta);
53  }
54  ustar = tang_windspeed2r * KAPPA / (std::log(1.5 * delta / z0) - psi_m);
55 
56  // prevent some unphysical math
57  if (!(ustar > zero && !std::isnan(ustar))) { ustar = zero; }
58  if (!(ustar < 2.0 && !std::isnan(ustar))) { ustar = 2.0; }
59  if (psi_m > std::log(myhalf * delta / z0)) { psi_m = std::log(myhalf * delta / z0); }
60 
61  Real uTarget = (1 - t_blank) * ustar / KAPPA * (std::log(myhalf * delta / z0) - psi_m);
62  Real u1Target = uTarget * u1_2r / (tiny + tang_windspeed2r);
63 
64  return u1Target;
65 }
KAPPA
constexpr amrex::Real KAPPA
Definition: ERF_Constants.H:63
zero
constexpr amrex::Real zero
Definition: ERF_Constants.H:8
myhalf
constexpr amrex::Real myhalf
Definition: ERF_Constants.H:13
CONST_GRAV
constexpr amrex::Real CONST_GRAV
Definition: ERF_Constants.H:64
z0
Real z0
Definition: ERF_InitCustomPertVels_ScalarAdvDiff.H:8
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_MOSTStress.H:40
similarity_funs::calc_psi_m
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real calc_psi_m(amrex::Real zeta) const
Definition: ERF_MOSTStress.H:105
similarity_funs::calc_psi_h
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real calc_psi_h(amrex::Real zeta) const
Definition: ERF_MOSTStress.H:124

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)

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  const Real tiny = std::numeric_limits<amrex::Real>::epsilon();
993  const Real U_s = one; // unit velocity scale
994 
995  // MOST parameters
996  similarity_funs sfuns;
997  const Real ggg = CONST_GRAV;
998  const Real kappa = KAPPA;
999  const Real z0 = solverChoice.if_z0;
1000  const Real tflux = solverChoice.if_surf_temp_flux;
1001  const Real init_surf_temp = solverChoice.if_init_surf_temp;
1002  const Real surf_heating_rate = solverChoice.if_surf_heating_rate;
1003  const Real Olen_in = solverChoice.if_Olen_in;
1004 
1005  ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
1006  {
1007  const Real t_blank = t_blank_arr(i, j, k);
1008  const Real t_blank_below = t_blank_arr(i, j, k-1);
1009  const Real t_blank_above = t_blank_arr(i, j, k+1);
1010  const Real t_blank_north = t_blank_arr(i , j+1, k);
1011  const Real t_blank_south = t_blank_arr(i , j-1, k);
1012  const Real t_blank_east = t_blank_arr(i+1, j , k);
1013  const Real t_blank_west = t_blank_arr(i-1, j , k);
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  Real drag_coefficient = alpha_h / std::pow(dx_x*dx_y*dx_z, one/three);
1017 
1018  // SURFACE TEMP AND HEATING/COOLING RATE
1019  if (init_surf_temp > zero) {
1020  const Real surf_temp = init_surf_temp + surf_heating_rate*time;
1021  if (t_blank > 0 && (t_blank_above == zero) && (t_blank_below == one)) { // building roof
1022  const Real bc_forcing_rt_srf = -(cell_data(i,j,k,Rho_comp) * surf_temp - cell_data(i,j,k,RhoTheta_comp));
1023  cell_src(i, j, k, RhoTheta_comp) -= drag_coefficient * U_s * bc_forcing_rt_srf;
1024 
1025  } else if (((t_blank > zero && t_blank < t_blank_west && t_blank_east == zero) ||
1026  (t_blank > zero && t_blank < t_blank_east && t_blank_west == zero) ||
1027  (t_blank > zero && t_blank < t_blank_north && t_blank_south == zero) ||
1028  (t_blank > zero && t_blank < t_blank_south && t_blank_north == zero))) {
1029  // this should enter for just building walls
1030  // walls are currently separated to allow for flexibility in the future to heat walls differently
1031 
1032  // south face
1033  if ((t_blank < t_blank_north) && (t_blank_north == one)) {
1034  const Real bc_forcing_rt_srf = -(cell_data(i,j,k,Rho_comp) * surf_temp - cell_data(i,j,k,RhoTheta_comp));
1035  cell_src(i, j, k, RhoTheta_comp) -= drag_coefficient * U_s * bc_forcing_rt_srf;
1036  }
1037 
1038  // north face
1039  if ((t_blank < t_blank_south) && (t_blank_south == one)) {
1040  const Real bc_forcing_rt_srf = -(cell_data(i,j,k,Rho_comp) * surf_temp - cell_data(i,j,k,RhoTheta_comp));
1041  cell_src(i, j, k, RhoTheta_comp) -= drag_coefficient * U_s * bc_forcing_rt_srf;
1042  }
1043 
1044  // west face
1045  if ((t_blank < t_blank_east) && (t_blank_east == one)) {
1046  const Real bc_forcing_rt_srf = -(cell_data(i,j,k,Rho_comp) * surf_temp - cell_data(i,j,k,RhoTheta_comp));
1047  cell_src(i, j, k, RhoTheta_comp) -= drag_coefficient * U_s * bc_forcing_rt_srf;
1048  }
1049 
1050  // east face
1051  if ((t_blank < t_blank_west) && (t_blank_west == one)) {
1052  const Real bc_forcing_rt_srf = -(cell_data(i,j,k,Rho_comp) * surf_temp - cell_data(i,j,k,RhoTheta_comp));
1053  cell_src(i, j, k, RhoTheta_comp) -= drag_coefficient * U_s * bc_forcing_rt_srf;
1054  }
1055 
1056  }
1057  }
1058 
1059  // SURFACE HEAT FLUX
1060  if (tflux != Real(1.e-8)){
1061  const Real ux_cc_2r = myhalf * (u(i ,j ,k+1) + u(i+1,j ,k+1));
1062  const Real uy_cc_2r = myhalf * (v(i ,j ,k+1) + v(i ,j+1,k+1));
1063  const Real h_windspeed2r = std::sqrt(ux_cc_2r * ux_cc_2r + uy_cc_2r * uy_cc_2r);
1064 
1065  const Real theta = cell_data(i,j,k ,RhoTheta_comp) / cell_data(i,j,k ,Rho_comp);
1066  Real theta_neighbor = cell_data(i,j,k+1,RhoTheta_comp) / cell_data(i,j,k+1,Rho_comp);
1067 
1068  if (t_blank > zero && (t_blank_above == zero)) { // building roof
1069  Real psi_m = zero;
1070  Real psi_h = zero;
1071  Real psi_h_neighbor = zero;
1072  Real ustar = h_windspeed2r * kappa / (std::log((1.5) * dx_z / z0) - psi_m);
1073  Real Olen = (Olen_in != Real(1e-8)) ? Olen_in : -ustar * ustar * ustar * theta / (kappa * ggg * tflux + tiny);
1074 
1075  for (int iter = 0; iter < 2; ++iter) {
1076  if (iter > 0) { Olen = -ustar * ustar * ustar * theta / (kappa * ggg * tflux + tiny); }
1077  Real zeta = (myhalf) * dx_z / Olen;
1078  Real zeta_neighbor = (1.5) * dx_z / Olen;
1079 
1080  // similarity functions
1081  psi_m = sfuns.calc_psi_m(zeta);
1082  psi_h = sfuns.calc_psi_h(zeta);
1083  psi_h_neighbor = sfuns.calc_psi_h(zeta_neighbor);
1084  ustar = h_windspeed2r * kappa / (std::log((1.5) * dx_z / z0) - psi_m);
1085  }
1086 
1087  // prevent some unphysical math
1088  if (!(ustar > zero && !std::isnan(ustar))) { ustar = zero; }
1089  if (!(ustar < 2.0 && !std::isnan(ustar))) { ustar = 2.0; }
1090  if (psi_h_neighbor > std::log(1.5 * dx_z / z0)) { psi_h_neighbor = std::log(1.5 * dx_z / z0); }
1091  if (psi_h > std::log(myhalf * dx_z / z0)) { psi_h = std::log(myhalf * dx_z / z0); }
1092 
1093  // We do not know the actual temperature so use cell above
1094  const Real thetastar = theta * ustar * ustar / (kappa * ggg * Olen);
1095  const Real surf_temp = theta_neighbor - thetastar / kappa * (std::log((1.5) * dx_z / z0) - psi_h_neighbor);
1096  const Real tTarget = surf_temp + thetastar / kappa * (std::log((myhalf) * dx_z / z0) - psi_h);
1097 
1098  const Real bc_forcing_rt = -(cell_data(i,j,k,Rho_comp) * tTarget - cell_data(i,j,k,RhoTheta_comp));
1099  cell_src(i, j, k, RhoTheta_comp) -= drag_coefficient * U_s * bc_forcing_rt;
1100 
1101  } else if (((t_blank > zero && t_blank < t_blank_west && t_blank_east == zero) ||
1102  (t_blank > zero && t_blank < t_blank_east && t_blank_west == zero) ||
1103  (t_blank > zero && t_blank < t_blank_north && t_blank_south == zero) ||
1104  (t_blank > zero && t_blank < t_blank_south && t_blank_north == zero))) { // this should enter for just building walls
1105 
1106  Real ux_cellaway = zero;
1107  Real uy_cellaway = zero;
1108  Real uz_cellaway = zero;
1109  Real u1 = zero;
1110  Real u2 = zero;
1111  Real delta = zero;
1112 
1113  // south face
1114  if (t_blank > zero && t_blank < t_blank_north && t_blank_south == zero) {
1115  ux_cellaway = myhalf * (u(i ,j-1,k) + u(i+1,j-1,k ));
1116  uz_cellaway = myhalf * (w(i ,j-1,k) + w(i ,j-1,k+1));
1117  u1 = ux_cellaway;
1118  u2 = uz_cellaway;
1119  delta = dx_y;
1120 
1121  // MOST
1122  theta_neighbor = cell_data(i,j-1,k,RhoTheta_comp) / cell_data(i,j-1,k,Rho_comp);
1123  }
1124 
1125  // north face
1126  if (t_blank > zero && t_blank < t_blank_south && t_blank_north == zero) {
1127  ux_cellaway = myhalf * (u(i ,j+1,k) + u(i+1,j+1,k ));
1128  uz_cellaway = myhalf * (w(i ,j+1,k) + w(i ,j+1,k+1));
1129  u1 = ux_cellaway;
1130  u2 = uz_cellaway;
1131  delta = dx_y;
1132 
1133  // MOST
1134  theta_neighbor = cell_data(i,j+1,k,RhoTheta_comp) / cell_data(i,j+1,k,Rho_comp);
1135  }
1136 
1137  // west face
1138  if (t_blank > zero && t_blank < t_blank_east && t_blank_west == zero) {
1139  uy_cellaway = myhalf * (v(i-1,j ,k) + v(i-1,j+1,k ));
1140  uz_cellaway = myhalf * (w(i-1,j ,k) + w(i-1,j ,k+1));
1141  u1 = uy_cellaway;
1142  u2 = uz_cellaway;
1143  delta = dx_x;
1144 
1145  // MOST
1146  theta_neighbor = cell_data(i-1,j,k,RhoTheta_comp) / cell_data(i-1,j,k,Rho_comp);
1147  }
1148 
1149  // east face
1150  if (t_blank > zero && t_blank < t_blank_west && t_blank_east == zero) {
1151  uy_cellaway = myhalf * (v(i+1,j ,k) + v(i+1,j+1,k ));
1152  uz_cellaway = myhalf * (w(i+1,j ,k) + w(i+1,j ,k+1));
1153  u1 = uy_cellaway;
1154  u2 = uz_cellaway;
1155  delta = dx_x;
1156 
1157  // MOST
1158  theta_neighbor = cell_data(i+1,j,k,RhoTheta_comp) / cell_data(i+1,j,k,Rho_comp);
1159  }
1160 
1161  Real tan_wspd = std::sqrt(u1 * u1 + u2 * u2);
1162 
1163  Real psi_m = zero;
1164  Real psi_h = zero;
1165  Real psi_h_neighbor = zero;
1166  Real ustar = tan_wspd * kappa / (std::log(1.5 * delta / z0) - psi_m);
1167  Real Olen = (Olen_in != Real(1e-8)) ? Olen_in : -ustar * ustar * ustar * theta / (kappa * ggg * tflux + tiny);
1168 
1169  for (int iter = 0; iter < 2; ++iter) {
1170  if (iter > 0) { Olen = -ustar * ustar * ustar * theta / (kappa * ggg * tflux + tiny); }
1171  Real zeta = (myhalf) * delta / Olen;
1172  Real zeta_neighbor = (1.5) * delta / Olen;
1173 
1174  // similarity functions
1175  psi_m = sfuns.calc_psi_m(zeta);
1176  psi_h = sfuns.calc_psi_h(zeta);
1177  psi_h_neighbor = sfuns.calc_psi_h(zeta_neighbor);
1178  ustar = tan_wspd * kappa / (std::log((1.5) * delta / z0) - psi_m);
1179  }
1180 
1181  // prevent some unphysical math
1182  if (!(ustar > zero && !std::isnan(ustar))) { ustar = zero; }
1183  if (!(ustar < 2.0 && !std::isnan(ustar))) { ustar = 2.0; }
1184  if (psi_h_neighbor > std::log(1.5 * delta / z0)) { psi_h_neighbor = std::log(1.5 * delta / z0); }
1185  if (psi_h > std::log(myhalf * delta / z0)) { psi_h = std::log(myhalf * delta / z0); }
1186 
1187  // We do not know the actual temperature so use cell above
1188  const Real thetastar = theta * ustar * ustar / (kappa * ggg * Olen);
1189  const Real surf_temp = theta_neighbor - thetastar / kappa * (std::log((1.5) * delta / z0) - psi_h_neighbor);
1190  const Real tTarget = surf_temp + thetastar / kappa * (std::log((myhalf) * delta / z0) - psi_h);
1191 
1192  const Real bc_forcing_rt = -(cell_data(i,j,k,Rho_comp) * tTarget - cell_data(i,j,k,RhoTheta_comp));
1193  cell_src(i, j, k, RhoTheta_comp) -= drag_coefficient * U_s * bc_forcing_rt;
1194  }
1195  }
1196 
1197  // Force fully immersed cells to planar average rho and theta
1198  if (t_blank == 1.0 && r_avg && t_avg) {
1199  const Real rho_avg = r_avg(k);
1200  const Real theta_avg = t_avg(k) / rho_avg; // Convert from RhoTheta to Theta
1201  const Real rho_cell = cell_data(i,j,k,Rho_comp);
1202  const Real bc_forcing_r = -(rho_avg - rho_cell);
1203  const Real bc_forcing_rt = -(rho_avg * theta_avg - cell_data(i,j,k,RhoTheta_comp));
1204 
1205  cell_src(i, j, k, Rho_comp) -= drag_coefficient * U_s * bc_forcing_r;
1206  cell_src(i, j, k, RhoTheta_comp) -= drag_coefficient * U_s * bc_forcing_rt;
1207  }
1208  });
1209 }
three
constexpr amrex::Real three
Definition: ERF_Constants.H:11
one
constexpr amrex::Real one
Definition: ERF_Constants.H:9
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);})
w
Real w
Definition: ERF_Plotfile2DInterpolator.cpp:22
LsmVar_SLM::theta
@ theta
Definition: ERF_SLM.H:20
SolverChoice::if_Olen_in
amrex::Real if_Olen_in
Input Obukhov length for immersed-forcing MOST [m].
Definition: ERF_DataStruct.H:1933
SolverChoice::if_z0
amrex::Real if_z0
Immersed-forcing roughness length [m].
Definition: ERF_DataStruct.H:1929
SolverChoice::if_Cd_scalar
amrex::Real if_Cd_scalar
Immersed-forcing drag coefficient for scalars.
Definition: ERF_DataStruct.H:1924
SolverChoice::if_init_surf_temp
amrex::Real if_init_surf_temp
Initial immersed-forcing surface temperature [K].
Definition: ERF_DataStruct.H:1931
SolverChoice::if_surf_temp_flux
amrex::Real if_surf_temp_flux
Immersed-forcing surface temperature flux [K m/s].
Definition: ERF_DataStruct.H:1930
SolverChoice::if_surf_heating_rate
amrex::Real if_surf_heating_rate
Immersed-forcing surface heating rate [K/hr].
Definition: ERF_DataStruct.H:1932

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

347 {
348  // geometric properties
349  const Real* dx_arr = geom.CellSize();
350  const Real dx_x = dx_arr[0];
351  const Real dx_y = dx_arr[1];
352  const Real dt = fac;
353 
354  const Real alpha_m = solverChoice.if_Cd_momentum;
355  const Real tiny = std::numeric_limits<amrex::Real>::epsilon();
356  const Real U_s = one; // unit velocity scale
357 
358  // MOST parameters
359  const Real z0 = solverChoice.if_z0;
360  const Real tflux_in = solverChoice.if_surf_temp_flux;
361  const Real Olen_in = solverChoice.if_Olen_in;
362  const bool l_use_most = solverChoice.if_use_most;
363  const bool l_stability_correction = solverChoice.if_stability_correction;
364 
365  // To limit stiffness of drag when using anelastic
366  const Real ws_floor = solverChoice.if_ws_floor;
367  const Real damp_alpha = solverChoice.if_damp_alpha;
368  // Point-implicit alternative to the clamp above; stabilizes both compressible and anelastic
369  const bool l_implicit_drag = solverChoice.if_implicit_drag;
370 
371  const bool is_slow_step = true; // This is determined by calling context
372  const bool use_ImmersedForcing_fast = solverChoice.immersed_forcing_substep;
373  const Real small_volfrac = 0.005;
374 
375  ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
376  {
377  const Real ux = u(i, j, k );
378  const Real uy = fourth * ( v(i, j , k ) + v(i-1, j , k )
379  + v(i, j+1, k ) + v(i-1, j+1, k ) );
380  const Real uz = fourth * ( w(i, j , k ) + w(i-1, j , k )
381  + w(i, j , k+1) + w(i-1, j , k+1) );
382  const amrex::Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz);
383 
384  const Real rho_xface = myhalf * ( cell_data(i,j,k,Rho_comp) + cell_data(i-1,j,k,Rho_comp) );
385  const Real theta_xface = (myhalf * (cell_data(i,j,k,RhoTheta_comp) + cell_data(i-1,j,k, RhoTheta_comp))) / rho_xface;
386 
387  // Use face-centered terrain_blanking if available, otherwise average from cell centers with threshold
388  Real t_blank_raw = (t_blank_xface_arr) ? t_blank_xface_arr(i, j , k ) :
389  myhalf * (t_blank_arr(i, j , k ) + t_blank_arr(i-1, j , k ));
390  const Real t_blank = (t_blank_raw < small_volfrac) ? zero : t_blank_raw;
391 
392  Real t_blank_below_raw = (k == 0) ? zero : (t_blank_xface_arr) ? t_blank_xface_arr(i, j , k-1) :
393  myhalf * (t_blank_arr(i, j , k-1) + t_blank_arr(i-1, j , k-1));
394  const Real t_blank_below = (t_blank_below_raw < small_volfrac) ? zero : t_blank_below_raw;
395 
396  Real t_blank_above_raw = (t_blank_xface_arr) ? t_blank_xface_arr(i, j , k+1) :
397  myhalf * (t_blank_arr(i, j , k+1) + t_blank_arr(i-1, j , k+1));
398  const Real t_blank_above = (t_blank_above_raw < small_volfrac) ? zero : t_blank_above_raw;
399 
400  Real t_blank_north_raw = (t_blank_xface_arr) ? t_blank_xface_arr(i, j+1, k ) :
401  myhalf * (t_blank_arr(i, j+1, k ) + t_blank_arr(i-1, j+1, k ));
402  const Real t_blank_north = (t_blank_north_raw < small_volfrac) ? zero : t_blank_north_raw;
403 
404  Real t_blank_south_raw = (t_blank_xface_arr) ? t_blank_xface_arr(i, j-1, k ) :
405  myhalf * (t_blank_arr(i, j-1, k ) + t_blank_arr(i-1, j-1, k ));
406  const Real t_blank_south = (t_blank_south_raw < small_volfrac) ? zero : t_blank_south_raw;
407 
408  const Real dx_z = (z_cc_arr) ? (z_cc_arr(i,j,k) - z_cc_arr(i,j,k-1)) : dx_arr[2];
409  const Real drag_coefficient = alpha_m / std::pow(dx_x*dx_y*dx_z, one/three);
410  const Real CdM = std::min(drag_coefficient / (windspeed + tiny), drag_coefficient);
411 
412  const Real roof_mask = (t_blank > zero && t_blank < t_blank_below && t_blank_above == zero && l_use_most) ? one : zero; // roof cell
413  const Real south_mask = (t_blank > zero && t_blank <= t_blank_north && t_blank_south == zero && l_use_most) ? one : zero; // south wall cell
414  const Real north_mask = (t_blank > zero && t_blank <= t_blank_south && t_blank_north == zero && l_use_most) ? one : zero; // north wall cell
415  const Real wall_mask = (t_blank > zero && t_blank < one && !l_use_most) ? one : zero; // all walls when NOT using MOST
416  const Real most_mask = roof_mask + south_mask + north_mask; // cells getting MOST treatment
417  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)
418  const Real interior_mask = (t_blank == 1.0) ? one : zero; // interior cell
419 
420  Real drag = zero;
421  Real u1_cellaway = zero;
422  Real u2_cellaway = zero;
423  Real rho_xface_inside = rho_xface;
424  Real theta_surf = theta_xface;
425  Real bc_forcing_x = zero;
426  Real u_target = zero;
427 
428  // roof forcing
429  if (roof_mask == one) {
430  u1_cellaway = u(i, j, k+1) ;
431  u2_cellaway = fourth * ( v(i, j , k+1) + v(i-1, j , k+1)
432  + v(i, j+1, k+1) + v(i-1, j+1, k+1) ) ;
433  rho_xface_inside = myhalf * (cell_data(i,j,k-1,Rho_comp) + cell_data(i-1,j,k-1,Rho_comp));
434  theta_surf = (myhalf * (cell_data(i,j,k-1,RhoTheta_comp) + cell_data(i-1,j,k-1, RhoTheta_comp))) / rho_xface_inside;
435  u_target = compute_if_most_target_vel(u1_cellaway, u2_cellaway, dx_z, z0, t_blank, theta_xface, theta_surf, tflux_in, Olen_in, l_stability_correction);
436  bc_forcing_x = -(u_target - ux); // BC forcing pushes nonrelative velocity toward target velocity
437  drag += bc_forcing_x * roof_mask * rho_xface * CdM * U_s;
438  }
439 
440  // south wall forcing
441  if (south_mask == one) {
442  u1_cellaway = u(i, j-1, k );
443  u2_cellaway = fourth * ( w(i, j-1, k ) + w(i-1, j-1, k )
444  + w(i, j-1, k+1) + w(i-1, j-1, k+1) ) ;
445  rho_xface_inside = myhalf * ( cell_data(i,j+1,k,Rho_comp) + cell_data(i-1,j+1,k,Rho_comp) );
446  theta_surf = (myhalf * (cell_data(i,j+1,k,RhoTheta_comp) + cell_data(i-1,j+1,k, RhoTheta_comp))) / rho_xface_inside;
447  u_target = compute_if_most_target_vel(u1_cellaway, u2_cellaway, dx_y, z0, t_blank, theta_xface, theta_surf, tflux_in, Olen_in, l_stability_correction);
448  bc_forcing_x = -(u_target - ux); // BC forcing pushes nonrelative velocity toward target velocity
449  drag += bc_forcing_x * south_mask * rho_xface * CdM * U_s;
450  }
451 
452  // north wall forcing
453  if (north_mask == one) {
454  u1_cellaway = u(i, j+1, k ) ;
455  u2_cellaway = fourth * ( w(i, j+1, k ) + w(i-1, j+1, k )
456  + w(i, j+1, k+1) + w(i-1, j+1, k+1) ) ;
457  rho_xface_inside = myhalf * ( cell_data(i,j-1,k,Rho_comp) + cell_data(i-1,j-1,k,Rho_comp) );
458  theta_surf = (myhalf * (cell_data(i,j-1,k,RhoTheta_comp) + cell_data(i-1,j-1,k, RhoTheta_comp))) / rho_xface_inside;
459  u_target = compute_if_most_target_vel(u1_cellaway, u2_cellaway, dx_y, z0, t_blank, theta_xface, theta_surf, tflux_in, Olen_in, l_stability_correction);
460  bc_forcing_x = -(u_target - ux); // BC forcing pushes nonrelative velocity toward target velocity
461  drag += bc_forcing_x * north_mask * rho_xface * CdM * U_s;
462  }
463 
464  // wall forcing (if not using most) or east/west walls when using MOST
465  if (wall_mask == one || east_west_mask == one) {
466  drag += (wall_mask + east_west_mask) * t_blank * rho_xface * CdM * ux * windspeed;
467  }
468 
469  // interior cell forcing
470  if (interior_mask == one) {
471  drag += interior_mask * rho_xface * CdM * ux * windspeed;
472  }
473 
474  if (l_implicit_drag) {
475  // point-implicit rescale of the aggregated drag
476  const Real lambda = CdM * ( (roof_mask + south_mask + north_mask) * U_s
477  + (wall_mask + east_west_mask) * t_blank * windspeed
478  + interior_mask * windspeed );
479  xmom_src_arr(i,j,k) -= drag / (one + lambda*dt);
480  } else if (is_slow_step && !use_ImmersedForcing_fast) {
481  // limit drag term for anelastic for numerical stability
482  Real d_drag = dt * -drag; // time step * acceleration like tendency
483  Real wsmax_change = damp_alpha * amrex::max(amrex::Math::abs(ux), ws_floor); // aims to prevent oscillations around 0.
484  if (amrex::Math::abs(ux) < 0.1){ // no damping for smaller velocities
485  wsmax_change =one * amrex::max(amrex::Math::abs(ux), ws_floor);
486  }
487  d_drag = amrex::min(amrex::max(d_drag, -wsmax_change), wsmax_change);
488  xmom_src_arr(i,j,k) += d_drag / dt; // put back as limited tendency
489  } else {
490  xmom_src_arr(i, j, k) -= drag;
491  }
492  });
493 }
fourth
constexpr amrex::Real fourth
Definition: ERF_Constants.H:14
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:25
SolverChoice::if_Cd_momentum
amrex::Real if_Cd_momentum
Immersed-forcing drag coefficient for momentum.
Definition: ERF_DataStruct.H:1923
SolverChoice::if_damp_alpha
amrex::Real if_damp_alpha
Immersed-forcing damping coefficient.
Definition: ERF_DataStruct.H:1937
SolverChoice::if_ws_floor
amrex::Real if_ws_floor
Wind-speed floor for immersed-forcing MOST [m/s].
Definition: ERF_DataStruct.H:1936
SolverChoice::immersed_forcing_substep
bool immersed_forcing_substep
Whether immersed-forcing source terms are applied only during substeps.
Definition: ERF_DataStruct.H:1919
SolverChoice::if_use_most
bool if_use_most
Whether immersed-forcing MOST is enabled.
Definition: ERF_DataStruct.H:1934
SolverChoice::if_implicit_drag
bool if_implicit_drag
Definition: ERF_DataStruct.H:1927
SolverChoice::if_stability_correction
bool if_stability_correction
Whether immersed-forcing stability corrections are enabled.
Definition: ERF_DataStruct.H:1935

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

510 {
511  // geometric properties
512  const Real* dx_arr = geom.CellSize();
513  const Real dx_x = dx_arr[0];
514  const Real dx_y = dx_arr[1];
515  const Real dt = fac;
516 
517  const Real alpha_m = solverChoice.if_Cd_momentum;
518  const Real tiny = std::numeric_limits<amrex::Real>::epsilon();
519  const Real U_s = one; // unit velocity scale
520 
521  // MOST parameters
522  const Real z0 = solverChoice.if_z0;
523  const Real tflux_in = solverChoice.if_surf_temp_flux;
524  const Real Olen_in = solverChoice.if_Olen_in;
525  const bool l_use_most = solverChoice.if_use_most;
526  const bool l_stability_correction = solverChoice.if_stability_correction;
527 
528  // To limit stiffness of drag when using anelastic
529  const Real ws_floor = solverChoice.if_ws_floor;
530  const Real damp_alpha = solverChoice.if_damp_alpha;
531  // Point-implicit alternative to the clamp above; stabilizes both compressible and anelastic
532  const bool l_implicit_drag = solverChoice.if_implicit_drag;
533 
534  const bool is_slow_step = true; // This is determined by calling context
535  const bool use_ImmersedForcing_fast = solverChoice.immersed_forcing_substep;
536  const Real small_volfrac = 0.005;
537 
538  ParallelFor(tby, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
539  {
540  const Real ux = fourth * ( u(i , j , k ) + u(i , j-1, k )
541  + u(i+1, j , k ) + u(i+1, j-1, k ) );
542  const Real uy = v(i, j, k);
543  const Real uz = fourth * ( w(i , j , k ) + w(i , j-1, k )
544  + w(i , j , k+1) + w(i , j-1, k+1) );
545  const amrex::Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz);
546 
547  const Real rho_yface = myhalf * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j-1,k,Rho_comp) );
548  const Real theta_yface = (myhalf * (cell_data(i,j,k ,RhoTheta_comp) + cell_data(i,j-1,k,RhoTheta_comp))) / rho_yface;
549 
550  // Use face-centered terrain_blanking if available, otherwise average from cell centers with threshold
551  Real t_blank_raw = (t_blank_yface_arr) ? t_blank_yface_arr(i , j , k ) :
552  myhalf * (t_blank_arr(i , j , k ) + t_blank_arr(i , j-1, k ));
553  const Real t_blank = (t_blank_raw < small_volfrac) ? zero : t_blank_raw;
554 
555  Real t_blank_below_raw = (k == 0) ? zero : (t_blank_yface_arr) ? t_blank_yface_arr(i , j , k-1) :
556  myhalf * (t_blank_arr(i , j , k-1) + t_blank_arr(i , j-1, k-1));
557  const Real t_blank_below = (t_blank_below_raw < small_volfrac) ? zero : t_blank_below_raw;
558 
559  Real t_blank_above_raw = (t_blank_yface_arr) ? t_blank_yface_arr(i , j , k+1) :
560  myhalf * (t_blank_arr(i , j , k+1) + t_blank_arr(i , j-1, k+1));
561  const Real t_blank_above = (t_blank_above_raw < small_volfrac) ? zero : t_blank_above_raw;
562 
563  Real t_blank_east_raw = (t_blank_yface_arr) ? t_blank_yface_arr(i+1, j , k ) :
564  myhalf * (t_blank_arr(i+1, j , k ) + t_blank_arr(i+1, j-1, k ));
565  const Real t_blank_east = (t_blank_east_raw < small_volfrac) ? zero : t_blank_east_raw;
566 
567  Real t_blank_west_raw = (t_blank_yface_arr) ? t_blank_yface_arr(i-1, j , k ) :
568  myhalf * (t_blank_arr(i-1, j , k ) + t_blank_arr(i-1, j-1, k ));
569  const Real t_blank_west = (t_blank_west_raw < small_volfrac) ? zero : t_blank_west_raw;
570 
571  const Real dx_z = (z_cc_arr) ? (z_cc_arr(i,j,k) - z_cc_arr(i,j,k-1)) : dx_arr[2];
572  const Real drag_coefficient = alpha_m / std::pow(dx_x*dx_y*dx_z, one/three);
573  const Real CdM = std::min(drag_coefficient / (windspeed + tiny), drag_coefficient);
574 
575  const Real roof_mask = (t_blank > zero && t_blank < t_blank_below && t_blank_above == zero && l_use_most) ? one : zero; // roof cell
576  const Real west_mask = (t_blank > zero && t_blank <= t_blank_east && t_blank_west == zero && l_use_most) ? one : zero; // west wall cell
577  const Real east_mask = (t_blank > zero && t_blank <= t_blank_west && t_blank_east == zero && l_use_most) ? one : zero; // east wall cell
578  const Real wall_mask = (t_blank > zero && t_blank < one && !l_use_most) ? one : zero; // all walls when NOT using MOST
579  const Real most_mask = roof_mask + west_mask + east_mask; // cells getting MOST treatment
580  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)
581  const Real interior_mask = (t_blank == 1.0) ? one : zero; // interior cell
582 
583  Real drag = zero;
584  Real u1_cellaway = zero;
585  Real u2_cellaway = zero;
586  Real rho_yface_inside = rho_yface;
587  Real theta_surf = theta_yface;
588  Real bc_forcing_y = zero;
589  Real u_target = zero;
590 
591  // roof forcing
592  if (roof_mask == one) {
593  u1_cellaway = fourth * ( u(i , j , k+1) + u(i , j-1, k+1)
594  + u(i+1, j , k+1) + u(i+1, j-1, k+1) );
595  u2_cellaway = v(i, j, k+1);
596  rho_yface_inside = myhalf * ( cell_data(i,j,k-1,Rho_comp) + cell_data(i,j-1,k-1,Rho_comp) );
597  theta_surf = (myhalf * (cell_data(i,j,k-1,RhoTheta_comp) + cell_data(i,j-1,k-1,RhoTheta_comp))) / rho_yface_inside;
598  u_target = compute_if_most_target_vel(u1_cellaway, u2_cellaway, dx_z, z0, t_blank, theta_yface, theta_surf, tflux_in, Olen_in, l_stability_correction);
599  bc_forcing_y = -(u_target - uy); // BC forcing pushes nonrelative velocity toward target velocity
600  drag += bc_forcing_y * roof_mask * rho_yface * CdM * U_s;
601  }
602 
603  // west wall forcing
604  if (west_mask == one) {
605  u1_cellaway = v(i-1, j , k );
606  u2_cellaway = fourth * ( w(i-1, j , k ) + w(i-1, j-1, k )
607  + w(i-1, j , k+1) + w(i-1, j-1, k+1) );
608  rho_yface_inside = myhalf * ( cell_data(i+1,j,k,Rho_comp) + cell_data(i+1,j-1,k,Rho_comp) );
609  theta_surf = (myhalf * (cell_data(i+1,j,k,RhoTheta_comp) + cell_data(i+1,j-1,k,RhoTheta_comp))) / rho_yface_inside;
610  u_target = compute_if_most_target_vel(u1_cellaway, u2_cellaway, dx_x, z0, t_blank, theta_yface, theta_surf, tflux_in, Olen_in, l_stability_correction);
611  bc_forcing_y = -(u_target - uy); // BC forcing pushes nonrelative velocity toward target velocity
612  drag += bc_forcing_y * west_mask * rho_yface * CdM * U_s;
613  }
614 
615  // east wall forcing
616  if (east_mask == one) {
617  u1_cellaway = v(i+1, j , k );
618  u2_cellaway = fourth * ( w(i+1, j , k ) + w(i+1, j-1, k )
619  + w(i+1, j , k+1) + w(i+1, j-1, k+1) );
620  rho_yface_inside = myhalf * ( cell_data(i-1,j,k,Rho_comp) + cell_data(i-1,j-1,k,Rho_comp) );
621  theta_surf = (myhalf * (cell_data(i-1,j,k,RhoTheta_comp) + cell_data(i-1,j-1,k,RhoTheta_comp))) / rho_yface_inside;
622  u_target = compute_if_most_target_vel(u1_cellaway, u2_cellaway, dx_x, z0, t_blank, theta_yface, theta_surf, tflux_in, Olen_in, l_stability_correction);
623  bc_forcing_y = -(u_target - uy); // BC forcing pushes nonrelative velocity toward target velocity
624  drag += bc_forcing_y * east_mask * rho_yface * CdM * U_s;
625  }
626 
627  // wall forcing (if not using most) or north/south walls when using MOST
628  if (wall_mask == one || north_south_mask == one) {
629  drag += (wall_mask + north_south_mask) * t_blank * rho_yface * CdM * uy * windspeed;
630  }
631 
632  // interior cell forcing
633  if (interior_mask == one) {
634  drag += interior_mask * rho_yface * CdM * uy * windspeed;
635  }
636 
637  if (l_implicit_drag) {
638  // point-implicit rescale of the aggregated drag
639  const Real lambda = CdM * ( (roof_mask + west_mask + east_mask) * U_s
640  + (wall_mask + north_south_mask) * t_blank * windspeed
641  + interior_mask * windspeed );
642  ymom_src_arr(i,j,k) -= drag / (one + lambda*dt);
643  } else if (is_slow_step && !use_ImmersedForcing_fast) {
644  // limit drag term for anelastic for numerical stability
645  Real d_drag = dt * -drag; // time step * acceleration like tendency
646  Real wsmax_change = damp_alpha * amrex::max(amrex::Math::abs(uy), ws_floor); // aims to prevent oscillations around 0.
647  if (amrex::Math::abs(uy) < 0.1){ // no damping for smaller velocities
648  wsmax_change =one * amrex::max(amrex::Math::abs(uy), ws_floor);
649  }
650  d_drag = amrex::min(amrex::max(d_drag, -wsmax_change), wsmax_change);
651  ymom_src_arr(i,j,k) += d_drag / dt; // put back as limited tendency
652  } else {
653  ymom_src_arr(i, j, k) -= drag;
654  }
655  });
656 }

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

673 {
674  // geometric properties
675  const Real* dx_arr = geom.CellSize();
676  const Real dx_x = dx_arr[0];
677  const Real dx_y = dx_arr[1];
678  const Real dt = fac;
679 
680  const Real alpha_m = solverChoice.if_Cd_momentum;
681  const Real tiny = std::numeric_limits<amrex::Real>::epsilon();
682  const Real U_s = one; // unit velocity scale
683 
684  // MOST parameters
685  const Real z0 = solverChoice.if_z0;
686  const Real tflux_in = solverChoice.if_surf_temp_flux;
687  const Real Olen_in = solverChoice.if_Olen_in;
688  const bool l_use_most = solverChoice.if_use_most;
689  const bool l_stability_correction = solverChoice.if_stability_correction;
690 
691  // To limit stiffness of drag when using anelastic
692  const Real ws_floor = solverChoice.if_ws_floor;
693  const Real damp_alpha = solverChoice.if_damp_alpha;
694  // Point-implicit alternative to the clamp above; stabilizes both compressible and anelastic
695  const bool l_implicit_drag = solverChoice.if_implicit_drag;
696 
697  const bool is_slow_step = true; // This is determined by calling context
698  const bool use_ImmersedForcing_fast = solverChoice.immersed_forcing_substep;
699  const Real small_volfrac = 0.005;
700 
701  ParallelFor(tbz, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
702  {
703  const Real ux = fourth * ( u(i , j , k ) + u(i+1, j , k )
704  + u(i , j , k-1) + u(i+1, j , k-1) );
705  const Real uy = fourth * ( v(i, j , k ) + v(i, j+1, k )
706  + v(i, j , k-1) + v(i, j+1, k-1) );
707  const Real uz = w(i, j, k);
708  const amrex::Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz);
709 
710  const Real rho_zface = myhalf * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j,k-1,Rho_comp) );
711  const Real theta_zface = (myhalf * (cell_data(i,j,k,RhoTheta_comp) + cell_data(i,j,k-1,RhoTheta_comp))) / rho_zface;
712 
713  // Use face-centered terrain_blanking if available, otherwise average from cell centers with threshold
714  Real t_blank_raw = (t_blank_zface_arr) ? t_blank_zface_arr(i ,j , k ) :
715  myhalf * (t_blank_arr(i ,j , k) + t_blank_arr(i , j , k-1));
716  const Real t_blank = (t_blank_raw < small_volfrac) ? zero : t_blank_raw;
717 
718  Real t_blank_above_raw = (t_blank_zface_arr) ? t_blank_zface_arr(i ,j , k+1) :
719  myhalf * (t_blank_arr(i ,j , k) + t_blank_arr(i , j , k+1));
720  const Real t_blank_above = (t_blank_above_raw < small_volfrac) ? zero : t_blank_above_raw;
721 
722  Real t_blank_north_raw = (t_blank_zface_arr) ? t_blank_zface_arr(i , j+1, k ) :
723  myhalf * (t_blank_arr(i ,j+1, k) + t_blank_arr(i , j+1, k-1));
724  const Real t_blank_north = (t_blank_north_raw < small_volfrac) ? zero : t_blank_north_raw;
725 
726  Real t_blank_south_raw = (t_blank_zface_arr) ? t_blank_zface_arr(i , j-1, k ) :
727  myhalf * (t_blank_arr(i ,j-1, k) + t_blank_arr(i , j-1, k-1));
728  const Real t_blank_south = (t_blank_south_raw < small_volfrac) ? zero : t_blank_south_raw;
729 
730  Real t_blank_east_raw = (t_blank_zface_arr) ? t_blank_zface_arr(i+1, j , k ) :
731  myhalf * (t_blank_arr(i+1,j , k) + t_blank_arr(i+1, j , k-1));
732  const Real t_blank_east = (t_blank_east_raw < small_volfrac) ? zero : t_blank_east_raw;
733 
734  Real t_blank_west_raw = (t_blank_zface_arr) ? t_blank_zface_arr(i-1, j , k ) :
735  myhalf * (t_blank_arr(i-1,j , k) + t_blank_arr(i-1, j , k-1));
736  const Real t_blank_west = (t_blank_west_raw < small_volfrac) ? zero : t_blank_west_raw;
737 
738  const Real dx_z = (z_cc_arr) ? (z_cc_arr(i,j,k) - z_cc_arr(i,j,k-1)) : dx_arr[2];
739  const Real drag_coefficient = alpha_m / std::pow(dx_x*dx_y*dx_z, one/three);
740  const Real CdM = std::min(drag_coefficient / (windspeed + tiny), drag_coefficient);
741 
742  const Real south_mask = (t_blank > zero && t_blank <= t_blank_north && t_blank_south == zero && l_use_most && k >= 1) ? one : zero; // south wall cell
743  const Real north_mask = (t_blank > zero && t_blank <= t_blank_south && t_blank_north == zero && l_use_most && k >= 1) ? one : zero; // north wall cell
744  const Real west_mask = (t_blank > zero && t_blank <= t_blank_east && t_blank_west == zero && l_use_most && k >= 1) ? one : zero; // west wall cell
745  const Real east_mask = (t_blank > zero && t_blank <= t_blank_west && t_blank_east == zero && l_use_most && k >= 1) ? one : zero; // east wall cell
746  const Real wall_mask = (t_blank > zero && t_blank < one && !l_use_most) ? one : zero; // all walls when NOT using MOST
747  const Real roof_mask = (t_blank > zero && t_blank_above == zero && l_use_most) ? one : zero; // roof cell (horizontal surface) - uses simple drag
748  const Real interior_mask = (t_blank == 1.0) ? one : zero; // interior cell
749 
750  Real drag = zero;
751  Real u1_cellaway = zero;
752  Real u2_cellaway = zero;
753  Real rho_zface_inside = rho_zface;
754  Real theta_surf = theta_zface;
755  Real bc_forcing_z = zero;
756  Real u_target = zero;
757 
758  // south wall forcing
759  if (south_mask == one) {
760  u1_cellaway = fourth * ( u(i , j-1, k ) + u(i+1, j-1, k )
761  + u(i , j-1, k-1) + u(i+1, j-1, k-1) );
762  u2_cellaway = w(i, j-1, k);
763  rho_zface_inside = myhalf * ( cell_data(i,j+1,k,Rho_comp) + cell_data(i,j+1,k-1,Rho_comp) );
764  theta_surf = (myhalf * (cell_data(i,j+1,k,RhoTheta_comp) + cell_data(i,j+1,k-1,RhoTheta_comp))) / rho_zface_inside;
765  u_target = compute_if_most_target_vel(u1_cellaway, u2_cellaway, dx_y, z0, t_blank, theta_zface, theta_surf, tflux_in, Olen_in, l_stability_correction);
766  bc_forcing_z = -(u_target - uz); // BC forcing pushes nonrelative velocity toward target velocity
767  drag += bc_forcing_z * south_mask * rho_zface * CdM * U_s;
768  }
769 
770  // north wall forcing
771  if (north_mask == one) {
772  u1_cellaway = fourth * ( u(i , j+1, k ) + u(i+1, j+1, k )
773  + u(i , j+1, k-1) + u(i+1, j+1, k-1) );
774  u2_cellaway = w(i, j+1, k);
775  rho_zface_inside = myhalf * ( cell_data(i,j-1,k,Rho_comp) + cell_data(i,j-1,k-1,Rho_comp) );
776  theta_surf = (myhalf * (cell_data(i,j-1,k,RhoTheta_comp) + cell_data(i,j-1,k-1,RhoTheta_comp))) / rho_zface_inside;
777  u_target = compute_if_most_target_vel(u1_cellaway, u2_cellaway, dx_y, z0, t_blank, theta_zface, theta_surf, tflux_in, Olen_in, l_stability_correction);
778  bc_forcing_z = -(u_target - uz); // BC forcing pushes nonrelative velocity toward target velocity
779  drag += bc_forcing_z * north_mask * rho_zface * CdM * U_s;
780  }
781 
782  // west wall forcing
783  if (west_mask == one) {
784  u1_cellaway = fourth * ( v(i-1, j , k ) + v(i-1, j+1, k )
785  + v(i-1, j , k-1) + v(i-1, j+1, k-1) );
786  u2_cellaway = w(i-1, j, k);
787  rho_zface_inside = myhalf * ( cell_data(i+1,j,k,Rho_comp) + cell_data(i+1,j,k-1,Rho_comp) );
788  theta_surf = (myhalf * (cell_data(i+1,j,k,RhoTheta_comp) + cell_data(i+1,j,k-1,RhoTheta_comp))) / rho_zface_inside;
789  u_target = compute_if_most_target_vel(u1_cellaway, u2_cellaway, dx_x, z0, t_blank, theta_zface, theta_surf, tflux_in, Olen_in, l_stability_correction);
790  bc_forcing_z = -(u_target - uz); // BC forcing pushes nonrelative velocity toward target velocity
791  drag += bc_forcing_z * west_mask * rho_zface * CdM * U_s;
792  }
793 
794  // east wall forcing
795  if (east_mask == one) {
796  u1_cellaway = fourth * ( v(i+1, j , k ) + v(i+1, j+1, k )
797  + v(i+1, j , k-1) + v(i+1, j+1, k-1) );
798  u2_cellaway = w(i+1, j, k);
799  rho_zface_inside = myhalf * ( cell_data(i-1,j,k,Rho_comp) + cell_data(i-1,j,k-1,Rho_comp) );
800  theta_surf = (myhalf * (cell_data(i-1,j,k,RhoTheta_comp) + cell_data(i-1,j,k-1,RhoTheta_comp))) / rho_zface_inside;
801  u_target = compute_if_most_target_vel(u1_cellaway, u2_cellaway, dx_x, z0, t_blank, theta_zface, theta_surf, tflux_in, Olen_in, l_stability_correction);
802  bc_forcing_z = -(u_target - uz); // BC forcing pushes nonrelative velocity toward target velocity
803  drag += bc_forcing_z * east_mask * rho_zface * CdM * U_s;
804  }
805 
806  // wall forcing (if not using most) or roof when using MOST
807  if (wall_mask == one || roof_mask == one) {
808  drag += (wall_mask + roof_mask) * t_blank * rho_zface * CdM * uz * windspeed;
809  }
810 
811  // interior cell forcing
812  if (interior_mask == one) {
813  drag += interior_mask * rho_zface * CdM * uz * windspeed;
814  }
815 
816  if (l_implicit_drag) {
817  // point-implicit rescale of the aggregated drag
818  const Real lambda = CdM * ( (south_mask + north_mask + west_mask + east_mask) * U_s
819  + (wall_mask + roof_mask) * t_blank * windspeed
820  + interior_mask * windspeed );
821  zmom_src_arr(i,j,k) -= drag / (one + lambda*dt);
822  } else if (is_slow_step && !use_ImmersedForcing_fast) {
823  // limit drag term for anelastic for numerical stability
824  Real d_drag = dt * -drag; // time step * acceleration like tendency
825  Real wsmax_change = damp_alpha * amrex::max(amrex::Math::abs(uz), ws_floor); // aims to prevent oscillations around 0.
826  if (amrex::Math::abs(uz) < 0.1){ // no damping for smaller velocities
827  wsmax_change = one * amrex::max(amrex::Math::abs(uz), ws_floor);
828  }
829  d_drag = amrex::min(amrex::max(d_drag, -wsmax_change), wsmax_change);
830  zmom_src_arr(i,j,k) += d_drag / dt; // put back as limited tendency
831  } else {
832  zmom_src_arr(i, j, k) -= drag;
833  }
834  });
835 }

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)

852 {
853  // geometric properties
854  const Real* dx_arr = geom.CellSize();
855  const Real dx_x = dx_arr[0];
856  const Real dx_y = dx_arr[1];
857 
858  const Real alpha_h = solverChoice.if_Cd_scalar;
859  const Real tiny = std::numeric_limits<amrex::Real>::epsilon();
860  const Real U_s = one; // unit velocity scale
861 
862  // MOST parameters
863  similarity_funs sfuns;
864  const Real ggg = CONST_GRAV;
865  const Real kappa = KAPPA;
866  const Real z0 = solverChoice.if_z0;
867  const Real tflux = solverChoice.if_surf_temp_flux;
868  const Real init_surf_temp = solverChoice.if_init_surf_temp;
869 
870  // Note this has been converted to K / s when it was read in;
871  const Real surf_heating_rate = solverChoice.if_surf_heating_rate;
872 
873  const Real Olen_in = solverChoice.if_Olen_in;
874 
875  ParallelFor(bx, [=]
876  AMREX_GPU_DEVICE(int i, int j, int k) noexcept
877  {
878  const Real dx_z = (z_cc_arr) ? (z_cc_arr(i,j,k) - z_cc_arr(i,j,k-1)) : dx_arr[2];
879  const Real drag_coefficient = alpha_h / std::pow(dx_x*dx_y*dx_z, one/three);
880 
881  const Real t_blank = t_blank_arr(i, j, k);
882  const Real t_blank_above = t_blank_arr(i, j, k+1);
883  const Real ux_cc_2r = myhalf * (u(i ,j ,k+1) + u(i+1,j ,k+1));
884  const Real uy_cc_2r = myhalf * (v(i ,j ,k+1) + v(i ,j+1,k+1));
885  const Real h_windspeed2r = std::sqrt(ux_cc_2r * ux_cc_2r + uy_cc_2r * uy_cc_2r);
886 
887  const Real theta = cell_data(i,j,k ,RhoTheta_comp) / cell_data(i,j,k ,Rho_comp);
888  const Real theta_neighbor = cell_data(i,j,k+1,RhoTheta_comp) / cell_data(i,j,k+1,Rho_comp);
889 
890  // SURFACE TEMP AND HEATING/COOLING RATE
891  if (init_surf_temp > zero) {
892  if (t_blank > 0 && (t_blank_above == zero)) { // force to MOST value
893  const Real surf_temp = init_surf_temp + surf_heating_rate*time;
894  const Real bc_forcing_rt_srf = -(cell_data(i,j,k-1,Rho_comp) * surf_temp - cell_data(i,j,k-1,RhoTheta_comp));
895  cell_src(i, j, k, RhoTheta_comp) -= drag_coefficient * U_s * bc_forcing_rt_srf;
896  }
897  }
898 
899  // SURFACE HEAT FLUX
900  if (tflux != Real(1e-8)){
901  if (t_blank > 0 && (t_blank_above == zero)) { // force to MOST value
902  Real psi_m = zero;
903  Real psi_h = zero;
904  Real psi_h_neighbor = zero;
905  Real ustar = h_windspeed2r * kappa / (std::log((Real(1.5)) * dx_z / z0) - psi_m);
906  const Real Olen = -ustar * ustar * ustar * theta / (kappa * ggg * tflux + tiny);
907  const Real zeta = (myhalf) * dx_z / Olen;
908  const Real zeta_neighbor = (Real(1.5)) * dx_z / Olen;
909 
910  // similarity functions
911  psi_m = sfuns.calc_psi_m(zeta);
912  psi_h = sfuns.calc_psi_h(zeta);
913  psi_h_neighbor = sfuns.calc_psi_h(zeta_neighbor);
914  ustar = h_windspeed2r * kappa / (std::log((Real(1.5)) * dx_z / z0) - psi_m);
915 
916  // prevent some unphysical math
917  if (!(ustar > zero && !std::isnan(ustar))) { ustar = zero; }
918  if (!(ustar < two && !std::isnan(ustar))) { ustar = two; }
919  if (psi_h_neighbor > std::log(Real(1.5) * dx_z / z0)) { psi_h_neighbor = std::log(Real(1.5) * dx_z / z0); }
920  if (psi_h > std::log(myhalf * dx_z / z0)) { psi_h = std::log(myhalf * dx_z / z0); }
921 
922  // We do not know the actual temperature so use cell above
923  const Real thetastar = theta * ustar * ustar / (kappa * ggg * Olen);
924  const Real surf_temp = theta_neighbor - thetastar / kappa * (std::log((Real(1.5)) * dx_z / z0) - psi_h_neighbor);
925  const Real tTarget = surf_temp + thetastar / kappa * (std::log((myhalf) * dx_z / z0) - psi_h);
926 
927  const Real bc_forcing_rt = -(cell_data(i,j,k,Rho_comp) * tTarget - cell_data(i,j,k,RhoTheta_comp));
928  cell_src(i, j, k, RhoTheta_comp) -= drag_coefficient * U_s * bc_forcing_rt;
929  }
930  }
931 
932  // OBUKHOV LENGTH
933  if (Olen_in != Real(1e-8)){
934  if (t_blank > 0 && (t_blank_above == zero)) { // force to MOST value
935  const Real Olen = Olen_in;
936  const Real zeta = (myhalf) * dx_z / Olen;
937  const Real zeta_neighbor = (Real(1.5)) * dx_z / Olen;
938 
939  // similarity functions
940  const Real psi_m = sfuns.calc_psi_m(zeta);
941  const Real psi_h = sfuns.calc_psi_h(zeta);
942  const Real psi_h_neighbor = sfuns.calc_psi_h(zeta_neighbor);
943  const Real ustar = h_windspeed2r * kappa / (std::log((Real(1.5)) * dx_z / z0) - psi_m);
944 
945  // We do not know the actual temperature so use cell above
946  const Real thetastar = theta * ustar * ustar / (kappa * ggg * Olen);
947  const Real surf_temp = theta_neighbor - thetastar / kappa * (std::log((Real(1.5)) * dx_z / z0) - psi_h_neighbor);
948  const Real tTarget = surf_temp + thetastar / kappa * (std::log((myhalf) * dx_z / z0) - psi_h);
949 
950  const Real bc_forcing_rt = -(cell_data(i,j,k,Rho_comp) * tTarget - cell_data(i,j,k,RhoTheta_comp));
951  cell_src(i, j, k, RhoTheta_comp) -= drag_coefficient * U_s * bc_forcing_rt;
952  }
953  }
954 
955  // Force fully immersed cells to planar average rho and theta
956  if (t_blank == one && r_avg && t_avg) {
957  const Real rho_avg = r_avg(k);
958  const Real theta_avg = t_avg(k) / rho_avg; // Convert from RhoTheta to Theta
959  const Real rho_cell = cell_data(i,j,k,Rho_comp);
960  const Real bc_forcing_r = -(rho_avg - rho_cell);
961  const Real bc_forcing_rt = -(rho_avg * theta_avg - cell_data(i,j,k,RhoTheta_comp));
962 
963  cell_src(i, j, k, Rho_comp) -= drag_coefficient * U_s * bc_forcing_r;
964  cell_src(i, j, k, RhoTheta_comp) -= drag_coefficient * U_s * bc_forcing_rt;
965  }
966  });
967 }
two
constexpr amrex::Real two
Definition: ERF_Constants.H:10

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

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

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

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

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

295 {
296  // geometric properties
297  const Real* dx_arr = geom.CellSize();
298  const Real dx_x = dx_arr[0];
299  const Real dx_y = dx_arr[1];
300  const Real dt = fac; // fac is actually dt in the calling code
301 
302  const Real alpha_m = solverChoice.if_Cd_momentum;
303  const Real tiny = std::numeric_limits<amrex::Real>::epsilon();
304  const bool l_implicit_drag = solverChoice.if_implicit_drag;
305 
306  const Real small_volfrac = 0.005;
307 
308  ParallelFor(tbz, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
309  {
310  const Real ux = fourth * ( u(i , j , k ) + u(i+1, j , k )
311  + u(i , j , k-1) + u(i+1, j , k-1) );
312  const Real uy = fourth * ( v(i , j , k ) + v(i , j+1, k )
313  + v(i , j , k-1) + v(i , j+1, k-1) );
314  const Real uz = w(i, j, k);
315  const Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz);
316  // Use face-centered terrain_blanking if available, otherwise average from cell centers
317  Real t_blank_raw = (t_blank_zface_arr) ? t_blank_zface_arr(i, j, k) :
318  myhalf * (t_blank_arr(i, j, k) + t_blank_arr(i, j, k-1));
319  const Real t_blank = (t_blank_raw < small_volfrac) ? zero : t_blank_raw;
320 
321  const Real dx_z = (z_cc_arr) ? (z_cc_arr(i,j,k) - z_cc_arr(i,j,k-1)) : dx_arr[2];
322  const Real drag_coefficient = alpha_m / std::pow(dx_x*dx_y*dx_z, one/three);
323  const Real CdM = std::min(drag_coefficient / (windspeed + tiny), drag_coefficient);
324 
325  const Real rho_zface = myhalf * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j,k-1,Rho_comp) );
326  const Real lambda = t_blank * CdM * windspeed; // linear drag rate [1/s]
327  const Real fac_local = l_implicit_drag ? lambda / (one + lambda*dt) : lambda; // point-implicit rescale (else explicit)
328  zmom_src_arr(i, j, k) -= fac_local * rho_zface * uz;
329  });
330 }

Referenced by make_mom_sources().

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