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

Functions

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real compute_if_most_target_vel (const amrex::Real u1_2r, const amrex::Real u2_2r, const amrex::Real delta, const amrex::Real z0, const amrex::Real t_blank, const amrex::Real theta_face, const amrex::Real theta_surf, const amrex::Real tflux_in, const amrex::Real Olen_in, const bool stability_correction, const amrex::Real wind_floor, const amrex::Real psi_cap_factor)
 
void ImmersedForcingTerrain_Xmom (const Box &tbx, const Array4< const Real > &u, const Array4< const Real > &v, const Array4< const Real > &w, const Array4< const Real > &cell_data, const Array4< const Real > &t_blank_arr, const Array4< const Real > &t_blank_xface_arr, const Array4< const Real > &z_cc_arr, const Array4< Real > &xmom_src_arr, const Geometry &geom, const SolverChoice &solverChoice, const Real fac)
 
void ImmersedForcingTerrain_Ymom (const Box &tby, const Array4< const Real > &u, const Array4< const Real > &v, const Array4< const Real > &w, const Array4< const Real > &cell_data, const Array4< const Real > &t_blank_arr, const Array4< const Real > &t_blank_yface_arr, const Array4< const Real > &z_cc_arr, const Array4< Real > &ymom_src_arr, const Geometry &geom, const SolverChoice &solverChoice, const Real fac)
 
static ib_wall::WallCond make_wall_cond (const SolverChoice &sc, const Real time)
 
template<int dir>
void ImmersedForcingTerrain_HorizMom_FractionStress (const Box &tb, 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_face_arr, const Array4< Real > &mom_src_arr, const Geometry &geom, const SolverChoice &solverChoice, const Real dt, const Real time)
 
void ImmersedForcingTerrain_Xmom_FractionStress (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< Real > &xmom_src_arr, const Geometry &geom, const SolverChoice &solverChoice, const Real fac, const Real time)
 
void ImmersedForcingTerrain_Ymom_FractionStress (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< Real > &ymom_src_arr, const Geometry &geom, const SolverChoice &solverChoice, const Real fac, const Real time)
 
void ImmersedForcingTerrain_Scalar_FractionStress (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< Real > &cell_src, const Array4< Real > &wall_hfx, const Array4< Real > &wall_ktarget, const Geometry &geom, const SolverChoice &solverChoice, const Table1D< Real > &r_avg, const Table1D< Real > &t_avg, const Real dt, const Real time, const bool wall_tke, const Real Cmu0)
 
void ImmersedForcingTerrain_Zmom (const Box &tbz, const Array4< const Real > &u, const Array4< const Real > &v, const Array4< const Real > &w, const Array4< const Real > &cell_data, const Array4< const Real > &t_blank_arr, const Array4< const Real > &t_blank_zface_arr, const Array4< const Real > &z_cc_arr, const Array4< Real > &zmom_src_arr, const Geometry &geom, const SolverChoice &solverChoice, const Real fac)
 
void ImmersedForcingBuildings_Xmom (const Box &tbx, const Array4< const Real > &u, const Array4< const Real > &v, const Array4< const Real > &w, const Array4< const Real > &cell_data, const Array4< const Real > &t_blank_arr, const Array4< const Real > &t_blank_xface_arr, const Array4< const Real > &z_cc_arr, const Array4< Real > &xmom_src_arr, const Geometry &geom, const SolverChoice &solverChoice, const Real fac)
 
void ImmersedForcingBuildings_Ymom (const Box &tby, const Array4< const Real > &u, const Array4< const Real > &v, const Array4< const Real > &w, const Array4< const Real > &cell_data, const Array4< const Real > &t_blank_arr, const Array4< const Real > &t_blank_yface_arr, const Array4< const Real > &z_cc_arr, const Array4< Real > &ymom_src_arr, const Geometry &geom, const SolverChoice &solverChoice, const Real fac)
 
void ImmersedForcingBuildings_Zmom (const Box &tbz, const Array4< const Real > &u, const Array4< const Real > &v, const Array4< const Real > &w, const Array4< const Real > &cell_data, const Array4< const Real > &t_blank_arr, const Array4< const Real > &t_blank_zface_arr, const Array4< const Real > &z_cc_arr, const Array4< Real > &zmom_src_arr, const Geometry &geom, const SolverChoice &solverChoice, const Real fac)
 
void ImmersedForcingTerrain_Scalar (const Box &bx, const Array4< const Real > &u, const Array4< const Real > &v, const Array4< const Real > &cell_data, const Array4< const Real > &t_blank_arr, const Array4< const Real > &z_cc_arr, const Array4< Real > &cell_src, const Geometry &geom, const SolverChoice &solverChoice, const Table1D< Real > &r_avg, const Table1D< Real > &t_avg, const Real time)
 
void ImmersedForcingBuildings_Scalar (const Box &bx, const Array4< const Real > &u, const Array4< const Real > &v, const Array4< const Real > &w, const Array4< const Real > &cell_data, const Array4< const Real > &t_blank_arr, const Array4< const Real > &z_cc_arr, const Array4< Real > &cell_src, const Geometry &geom, const SolverChoice &solverChoice, const Table1D< Real > &r_avg, const Table1D< Real > &t_avg, const Real time)
 

Function Documentation

◆ compute_if_most_target_vel()

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

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

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

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)

1376 {
1377  // geometric properties
1378  const auto dx_arr = geom.CellSizeArray();
1379  const Real dx_x = dx_arr[0];
1380  const Real dx_y = dx_arr[1];
1381 
1382  const Real alpha_h = solverChoice.if_Cd_scalar;
1383  // erf.if_snap_partial_cells: read the cell blanking snapped to solid (1)
1384  // or fluid (0) at half, so a height-map building becomes the same
1385  // staircase of whole cells an exact box is: the thermal conditions sit
1386  // on the boundary solid cells, roofs included, found from the neighbour
1387  // blanking. Off (the default), the raw fractions are used and nothing
1388  // below changes. See SolverChoice::if_snap_partial_cells.
1389  const bool l_snap = solverChoice.if_snap_partial_cells;
1390  auto snapb = [=] AMREX_GPU_DEVICE (amrex::Real b) noexcept -> amrex::Real {
1391  return l_snap ? ((b >= myhalf) ? one : zero) : b;
1392  };
1393  const Real tiny = std::numeric_limits<amrex::Real>::epsilon();
1394  const Real U_s = one; // unit velocity scale
1395 
1396  // MOST parameters
1397  similarity_funs sfuns;
1398  const Real ggg = CONST_GRAV;
1399  const Real kappa = KAPPA;
1400  const Real z0 = solverChoice.if_z0;
1401  const Real tflux = solverChoice.if_surf_temp_flux;
1402  const Real init_surf_temp = solverChoice.if_init_surf_temp;
1403  const Real surf_heating_rate = solverChoice.if_surf_heating_rate;
1404  const Real stab_wind_floor = solverChoice.if_stability_wind_floor;
1405 
1406  ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
1407  {
1408  const Real t_blank = snapb(t_blank_arr(i, j, k));
1409  const Real t_blank_below = snapb(t_blank_arr(i, j, k-1));
1410  const Real t_blank_above = snapb(t_blank_arr(i, j, k+1));
1411  const Real t_blank_north = snapb(t_blank_arr(i , j+1, k));
1412  const Real t_blank_south = snapb(t_blank_arr(i , j-1, k));
1413  const Real t_blank_east = snapb(t_blank_arr(i+1, j , k));
1414  const Real t_blank_west = snapb(t_blank_arr(i-1, j , k));
1415 
1416  const Real dx_z = (z_cc_arr) ? (z_cc_arr(i,j,k) - z_cc_arr(i,j,k-1)) : dx_arr[2];
1417  Real drag_coefficient = alpha_h / std::pow(dx_x*dx_y*dx_z, one/three);
1418 
1419  // Wall cells of the building, named by the face they carry: a
1420  // partially blanked cell with a more solid neighbour behind it and a
1421  // fluid one in front, or, with the snap (every solid cell at 1), a
1422  // solid cell with a solid neighbour behind and a fluid one in front.
1423  const bool south_face = l_snap ? (t_blank == one && t_blank_north == one && t_blank_south == zero)
1424  : (t_blank > zero && t_blank < t_blank_north && t_blank_south == zero);
1425  const bool north_face = l_snap ? (t_blank == one && t_blank_south == one && t_blank_north == zero)
1426  : (t_blank > zero && t_blank < t_blank_south && t_blank_north == zero);
1427  const bool west_face = l_snap ? (t_blank == one && t_blank_east == one && t_blank_west == zero)
1428  : (t_blank > zero && t_blank < t_blank_east && t_blank_west == zero);
1429  const bool east_face = l_snap ? (t_blank == one && t_blank_west == one && t_blank_east == zero)
1430  : (t_blank > zero && t_blank < t_blank_west && t_blank_east == zero);
1431 
1432  // SURFACE TEMP AND HEATING/COOLING RATE
1433  if (init_surf_temp > zero) {
1434  const Real surf_temp = init_surf_temp + surf_heating_rate*time;
1435  if (t_blank > 0 && (t_blank_above == zero) && (t_blank_below == one)) { // building roof
1436  const Real bc_forcing_rt_srf = -(cell_data(i,j,k,Rho_comp) * surf_temp - cell_data(i,j,k,RhoTheta_comp));
1437  cell_src(i, j, k, RhoTheta_comp) -= drag_coefficient * U_s * bc_forcing_rt_srf;
1438 
1439  } else if (east_face || west_face || south_face || north_face) {
1440  // this should enter for just building walls
1441  // walls are currently separated to allow for flexibility in the future to heat walls differently
1442 
1443  // south face
1444  if (l_snap ? south_face : ((t_blank < t_blank_north) && (t_blank_north == one))) {
1445  const Real bc_forcing_rt_srf = -(cell_data(i,j,k,Rho_comp) * surf_temp - cell_data(i,j,k,RhoTheta_comp));
1446  cell_src(i, j, k, RhoTheta_comp) -= drag_coefficient * U_s * bc_forcing_rt_srf;
1447  }
1448 
1449  // north face
1450  if (l_snap ? north_face : ((t_blank < t_blank_south) && (t_blank_south == one))) {
1451  const Real bc_forcing_rt_srf = -(cell_data(i,j,k,Rho_comp) * surf_temp - cell_data(i,j,k,RhoTheta_comp));
1452  cell_src(i, j, k, RhoTheta_comp) -= drag_coefficient * U_s * bc_forcing_rt_srf;
1453  }
1454 
1455  // west face
1456  if (l_snap ? west_face : ((t_blank < t_blank_east) && (t_blank_east == one))) {
1457  const Real bc_forcing_rt_srf = -(cell_data(i,j,k,Rho_comp) * surf_temp - cell_data(i,j,k,RhoTheta_comp));
1458  cell_src(i, j, k, RhoTheta_comp) -= drag_coefficient * U_s * bc_forcing_rt_srf;
1459  }
1460 
1461  // east face
1462  if (l_snap ? east_face : ((t_blank < t_blank_west) && (t_blank_west == one))) {
1463  const Real bc_forcing_rt_srf = -(cell_data(i,j,k,Rho_comp) * surf_temp - cell_data(i,j,k,RhoTheta_comp));
1464  cell_src(i, j, k, RhoTheta_comp) -= drag_coefficient * U_s * bc_forcing_rt_srf;
1465  }
1466 
1467  }
1468  }
1469 
1470  // SURFACE HEAT FLUX
1471  if (tflux != Real(1.e-8)){
1472  const Real ux_cc_2r = myhalf * (u(i ,j ,k+1) + u(i+1,j ,k+1));
1473  const Real uy_cc_2r = myhalf * (v(i ,j ,k+1) + v(i ,j+1,k+1));
1474  // The heat transfer uses the wind floored at erf.if_stability_wind_floor and an Obukhov
1475  // length bounded so that |zeta| <= ib_stability::zeta_max(), as on flat ground
1476  const Real h_windspeed2r = ib_stability::floored_wind(std::sqrt(ux_cc_2r * ux_cc_2r + uy_cc_2r * uy_cc_2r), stab_wind_floor);
1477 
1478  const Real theta = cell_data(i,j,k ,RhoTheta_comp) / cell_data(i,j,k ,Rho_comp);
1479  Real theta_neighbor = cell_data(i,j,k+1,RhoTheta_comp) / cell_data(i,j,k+1,Rho_comp);
1480 
1481  if (t_blank > zero && (t_blank_above == zero)) { // building roof
1482  Real psi_m = zero;
1483  Real psi_h = zero;
1484  Real psi_h_neighbor = zero;
1485  Real ustar = h_windspeed2r * kappa / (std::log((1.5) * dx_z / z0) - psi_m);
1486  // (start-up refuses erf.if_Olen with a surface flux, so L always comes from the flux here)
1487  Real Olen = ib_stability::bounded_obukhov_length(-ustar * ustar * ustar * theta / (kappa * ggg * tflux + tiny), Real(1.5) * dx_z);
1488 
1489  for (int iter = 0; iter < 2; ++iter) {
1490  if (iter > 0) { Olen = ib_stability::bounded_obukhov_length(-ustar * ustar * ustar * theta / (kappa * ggg * tflux + tiny), Real(1.5) * dx_z); }
1491  Real zeta = ib_stability::bounded_zeta(myhalf * dx_z, Olen);
1492  Real zeta_neighbor = ib_stability::bounded_zeta(Real(1.5) * dx_z, Olen);
1493 
1494  // similarity functions
1495  psi_m = sfuns.calc_psi_m(zeta);
1496  psi_h = sfuns.calc_psi_h(zeta);
1497  psi_h_neighbor = sfuns.calc_psi_h(zeta_neighbor);
1498  ustar = h_windspeed2r * kappa / (std::log((1.5) * dx_z / z0) - psi_m);
1499  }
1500 
1501  // prevent some unphysical math
1502  ustar = ib_stability::clamped_ustar(ustar);
1503  // caps of ln(z / z0): erf.if_psi_cap_factor stays off the temperature forcing, whose
1504  // theta* comes from an L not made consistent with u* (the legacy law, kept as it is;
1505  // the fraction-stress law of erf-model/ERF#4134 forms them consistently)
1506  psi_h_neighbor = ib_stability::capped_psi(psi_h_neighbor, Real(1.5) * dx_z, z0, one);
1507  psi_h = ib_stability::capped_psi(psi_h, myhalf * dx_z, z0, one);
1508 
1509  // We do not know the actual temperature so use cell above
1510  const Real thetastar = theta * ustar * ustar / (kappa * ggg * Olen);
1511  const Real surf_temp = theta_neighbor - thetastar / kappa * (std::log((1.5) * dx_z / z0) - psi_h_neighbor);
1512  const Real tTarget = surf_temp + thetastar / kappa * (std::log((myhalf) * dx_z / z0) - psi_h);
1513 
1514  const Real bc_forcing_rt = -(cell_data(i,j,k,Rho_comp) * tTarget - cell_data(i,j,k,RhoTheta_comp));
1515  cell_src(i, j, k, RhoTheta_comp) -= drag_coefficient * U_s * bc_forcing_rt;
1516 
1517  } else if (east_face || west_face || south_face || north_face) { // this should enter for just building walls
1518 
1519  Real ux_cellaway = zero;
1520  Real uy_cellaway = zero;
1521  Real uz_cellaway = zero;
1522  Real u1 = zero;
1523  Real u2 = zero;
1524  Real delta = zero;
1525 
1526  // south face
1527  if (south_face) {
1528  ux_cellaway = myhalf * (u(i ,j-1,k) + u(i+1,j-1,k ));
1529  uz_cellaway = myhalf * (w(i ,j-1,k) + w(i ,j-1,k+1));
1530  u1 = ux_cellaway;
1531  u2 = uz_cellaway;
1532  delta = dx_y;
1533 
1534  // MOST
1535  theta_neighbor = cell_data(i,j-1,k,RhoTheta_comp) / cell_data(i,j-1,k,Rho_comp);
1536  }
1537 
1538  // north face
1539  if (north_face) {
1540  ux_cellaway = myhalf * (u(i ,j+1,k) + u(i+1,j+1,k ));
1541  uz_cellaway = myhalf * (w(i ,j+1,k) + w(i ,j+1,k+1));
1542  u1 = ux_cellaway;
1543  u2 = uz_cellaway;
1544  delta = dx_y;
1545 
1546  // MOST
1547  theta_neighbor = cell_data(i,j+1,k,RhoTheta_comp) / cell_data(i,j+1,k,Rho_comp);
1548  }
1549 
1550  // west face
1551  if (west_face) {
1552  uy_cellaway = myhalf * (v(i-1,j ,k) + v(i-1,j+1,k ));
1553  uz_cellaway = myhalf * (w(i-1,j ,k) + w(i-1,j ,k+1));
1554  u1 = uy_cellaway;
1555  u2 = uz_cellaway;
1556  delta = dx_x;
1557 
1558  // MOST
1559  theta_neighbor = cell_data(i-1,j,k,RhoTheta_comp) / cell_data(i-1,j,k,Rho_comp);
1560  }
1561 
1562  // east face
1563  if (east_face) {
1564  uy_cellaway = myhalf * (v(i+1,j ,k) + v(i+1,j+1,k ));
1565  uz_cellaway = myhalf * (w(i+1,j ,k) + w(i+1,j ,k+1));
1566  u1 = uy_cellaway;
1567  u2 = uz_cellaway;
1568  delta = dx_x;
1569 
1570  // MOST
1571  theta_neighbor = cell_data(i+1,j,k,RhoTheta_comp) / cell_data(i+1,j,k,Rho_comp);
1572  }
1573 
1574  Real tan_wspd = ib_stability::floored_wind(std::sqrt(u1 * u1 + u2 * u2), stab_wind_floor);
1575 
1576  Real psi_m = zero;
1577  Real psi_h = zero;
1578  Real psi_h_neighbor = zero;
1579  Real ustar = tan_wspd * kappa / (std::log(1.5 * delta / z0) - psi_m);
1580  Real Olen = ib_stability::bounded_obukhov_length(-ustar * ustar * ustar * theta / (kappa * ggg * tflux + tiny), Real(1.5) * delta);
1581 
1582  for (int iter = 0; iter < 2; ++iter) {
1583  if (iter > 0) { Olen = ib_stability::bounded_obukhov_length(-ustar * ustar * ustar * theta / (kappa * ggg * tflux + tiny), Real(1.5) * delta); }
1584  Real zeta = ib_stability::bounded_zeta(myhalf * delta, Olen);
1585  Real zeta_neighbor = ib_stability::bounded_zeta(Real(1.5) * delta, Olen);
1586 
1587  // similarity functions
1588  psi_m = sfuns.calc_psi_m(zeta);
1589  psi_h = sfuns.calc_psi_h(zeta);
1590  psi_h_neighbor = sfuns.calc_psi_h(zeta_neighbor);
1591  ustar = tan_wspd * kappa / (std::log((1.5) * delta / z0) - psi_m);
1592  }
1593 
1594  // prevent some unphysical math
1595  ustar = ib_stability::clamped_ustar(ustar);
1596  // caps of ln(z / z0), as above (erf-model/ERF#4134)
1597  psi_h_neighbor = ib_stability::capped_psi(psi_h_neighbor, Real(1.5) * delta, z0, one);
1598  psi_h = ib_stability::capped_psi(psi_h, myhalf * delta, z0, one);
1599 
1600  // We do not know the actual temperature so use cell above
1601  const Real thetastar = theta * ustar * ustar / (kappa * ggg * Olen);
1602  const Real surf_temp = theta_neighbor - thetastar / kappa * (std::log((1.5) * delta / z0) - psi_h_neighbor);
1603  const Real tTarget = surf_temp + thetastar / kappa * (std::log((myhalf) * delta / z0) - psi_h);
1604 
1605  const Real bc_forcing_rt = -(cell_data(i,j,k,Rho_comp) * tTarget - cell_data(i,j,k,RhoTheta_comp));
1606  cell_src(i, j, k, RhoTheta_comp) -= drag_coefficient * U_s * bc_forcing_rt;
1607  }
1608  }
1609 
1610  // Force fully immersed cells to planar average rho and theta
1611  if (t_blank == 1.0 && r_avg && t_avg) {
1612  const Real rho_avg = r_avg(k);
1613  const Real theta_avg = t_avg(k) / rho_avg; // Convert from RhoTheta to Theta
1614  const Real rho_cell = cell_data(i,j,k,Rho_comp);
1615  const Real bc_forcing_r = -(rho_avg - rho_cell);
1616  const Real bc_forcing_rt = -(rho_avg * theta_avg - cell_data(i,j,k,RhoTheta_comp));
1617 
1618  cell_src(i, j, k, Rho_comp) -= drag_coefficient * U_s * bc_forcing_r;
1619  cell_src(i, j, k, RhoTheta_comp) -= drag_coefficient * U_s * bc_forcing_rt;
1620  }
1621  });
1622 }
Rho_comp
#define Rho_comp
Definition: ERF_IndexDefines.H:53
RhoTheta_comp
#define RhoTheta_comp
Definition: ERF_IndexDefines.H:54
ParallelFor
ParallelFor(fab_box, [=] AMREX_GPU_DEVICE(int i, int j, int k) { qrcuten_arr(i, j, k)=Real(0);qscuten_arr(i, j, k)=Real(0);qicuten_arr(i, j, k)=Real(0);})
three
constexpr amrex::Real three
Definition: ERF_NumericalConstants.H:39
one
constexpr amrex::Real one
Definition: ERF_NumericalConstants.H:37
w
Real w
Definition: ERF_Plotfile2DInterpolator.cpp:23
MicVar_Kess::theta
@ theta
Definition: ERF_Kessler.H:26
SolverChoice::if_z0
amrex::Real if_z0
Immersed-forcing roughness length [m].
Definition: ERF_DataStruct.H:2706
SolverChoice::if_snap_partial_cells
bool if_snap_partial_cells
Definition: ERF_DataStruct.H:2712
SolverChoice::if_Cd_scalar
amrex::Real if_Cd_scalar
Immersed-forcing drag coefficient for scalars.
Definition: ERF_DataStruct.H:2700
SolverChoice::if_init_surf_temp
amrex::Real if_init_surf_temp
Initial immersed-forcing surface temperature [K].
Definition: ERF_DataStruct.H:2708
SolverChoice::if_surf_temp_flux
amrex::Real if_surf_temp_flux
Immersed-forcing surface temperature flux [K m/s].
Definition: ERF_DataStruct.H:2707
SolverChoice::if_surf_heating_rate
amrex::Real if_surf_heating_rate
Immersed-forcing surface heating rate [K/hr].
Definition: ERF_DataStruct.H:2709
SolverChoice::if_stability_wind_floor
amrex::Real if_stability_wind_floor
Wind-speed floor [m/s] behind the friction velocity of the wall law's stability estimate (0: none).
Definition: ERF_DataStruct.H:2721

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

642 {
643  // geometric properties
644  const auto dx_arr = geom.CellSizeArray();
645  const Real dx_x = dx_arr[0];
646  const Real dx_y = dx_arr[1];
647  const Real dt = fac;
648 
649  const Real alpha_m = solverChoice.if_Cd_momentum;
650  const Real tiny = std::numeric_limits<amrex::Real>::epsilon();
651  const Real U_s = one; // unit velocity scale
652 
653  // MOST parameters
654  const Real z0 = solverChoice.if_z0;
655  const Real tflux_in = solverChoice.if_surf_temp_flux;
656  const Real Olen_in = solverChoice.if_Olen_in;
657  const Real stab_wind_floor = solverChoice.if_stability_wind_floor;
658  const Real psi_cap_factor = solverChoice.if_psi_cap_factor;
659  const bool l_use_most = solverChoice.if_use_most;
660  const bool l_stability_correction = solverChoice.if_stability_correction;
661 
662  // To limit stiffness of drag when using anelastic
663  const Real ws_floor = solverChoice.if_ws_floor;
664  const Real damp_alpha = solverChoice.if_damp_alpha;
665  // Point-implicit alternative to the clamp above; stabilizes both compressible and anelastic
666  const bool l_implicit_drag = solverChoice.if_implicit_drag;
667  // erf.if_implicit_projection: the linear drag of the solid is applied inside the anelastic
668  // projection (ERF::make_if_projection_sigma), not here
669  const bool l_ip = solverChoice.if_implicit_projection;
670 
671  const bool is_slow_step = true; // This is determined by calling context
672  const bool use_ImmersedForcing_fast = solverChoice.immersed_forcing_substep;
673  const Real small_volfrac = solverChoice.if_small_volfrac;
674  // erf.if_snap_partial_cells: read the cell blanking snapped to solid (1)
675  // or fluid (0) at half, so a height-map building becomes the same
676  // staircase of whole cells an exact box is. A face is solid when either
677  // cell it joins is; a face between a solid and a fluid cell is
678  // wall-normal and gets the interior drag toward zero (no penetration); a
679  // face between two solid cells carries the roof or wall law of its row
680  // (the full log-law target, not the partial-cell weighted one) or the
681  // interior drag, never both; the partial-cell branches (wall_mask,
682  // east_west_mask and the like, faces with 0 < t_blank < 1) do not arise
683  // under the snap. Off (the default), the raw fractions
684  // are used and nothing below changes: a boundary face of an exact box
685  // carries the wall law and the interior drag together as before. See
686  // SolverChoice::if_snap_partial_cells.
687  const bool l_snap = solverChoice.if_snap_partial_cells;
688  // The face blanking and its neighbours are read through ERF_ImmersedBuildingMasks.H, as the
689  // implicit projection reads them
690  const int dom_klo = geom.Domain().smallEnd(2);
691 
692 
693  ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
694  {
695  const Real ux = u(i, j, k );
696  const Real uy = fourth * ( v(i, j , k ) + v(i-1, j , k )
697  + v(i, j+1, k ) + v(i-1, j+1, k ) );
698  const Real uz = fourth * ( w(i, j , k ) + w(i-1, j , k )
699  + w(i, j , k+1) + w(i-1, j , k+1) );
700  const amrex::Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz);
701 
702  const Real rho_xface = myhalf * ( cell_data(i,j,k,Rho_comp) + cell_data(i-1,j,k,Rho_comp) );
703  const Real theta_xface = (myhalf * (cell_data(i,j,k,RhoTheta_comp) + cell_data(i-1,j,k, RhoTheta_comp))) / rho_xface;
704 
705  // the face, its neighbours and the raw blanking the explicit path zeroes on (ERF_ImmersedBuildingMasks.H)
706  const if_bld::HorizFace fc = if_bld::horiz_face<0>(t_blank_arr, t_blank_xface_arr, i, j, k, dom_klo, l_snap, small_volfrac);
707  const Real t_blank = fc.t;
708  // With the snap the boundary solid face stands for the wall layer and
709  // takes the full log-law target; the partial-cell weight (1 - t_blank)
710  // of compute_if_most_target_vel() would make it zero (a no-slip
711  // staircase), so the snap path passes a zero blanking to it.
712  const Real t_blank_law = l_snap ? zero : t_blank;
713  const Real t_blank_below = fc.below;
714  const Real t_blank_above = fc.above;
715  const Real t_blank_south = fc.lo;
716  const Real t_blank_north = fc.hi;
717 
718  const Real dx_z = (z_cc_arr) ? (z_cc_arr(i,j,k) - z_cc_arr(i,j,k-1)) : dx_arr[2];
719  const Real drag_coefficient = alpha_m / std::pow(dx_x*dx_y*dx_z, one/three);
720  const Real CdM = std::min(drag_coefficient / (windspeed + tiny), drag_coefficient);
721  // With the snap a face joining a solid and a fluid cell is wall-normal:
722  // it gets the interior drag (no penetration), not a wall law, whatever
723  // row it lies in.
724  const bool normal_face = fc.normal;
725 
726  // With the snap every solid face has t_blank = 1: a roof face lies in
727  // the top solid row (t_blank <= t_blank_below, the face above fluid),
728  // and a face carrying a wall law is not also an interior face.
729  // roof, south/north wall law, partial-wall and interior drag (ERF_ImmersedBuildingMasks.H)
730  if_bld::HorizMasks mk = if_bld::horiz_masks(t_blank, t_blank_below, t_blank_above, t_blank_south, t_blank_north,
731  normal_face, l_use_most, l_snap);
732  // erf.if_implicit_projection: the linear drag is in the projection, and a face the explicit
733  // path would zero after the projection (solid on both sides) is held there, with no wall law
734  if (l_ip) {
735  mk.wall = zero; mk.side = zero; mk.interior = zero;
736  if (fc.raw == one) { mk.roof = zero; mk.lo = zero; mk.hi = zero; }
737  }
738  const Real roof_mask = mk.roof;
739  const Real south_mask = mk.lo;
740  const Real north_mask = mk.hi;
741  const Real wall_mask = mk.wall;
742  const Real east_west_mask = mk.side;
743  const Real interior_mask = mk.interior;
744 
745  Real drag = zero;
746  Real u1_cellaway = zero;
747  Real u2_cellaway = zero;
748  Real rho_xface_inside = rho_xface;
749  Real theta_surf = theta_xface;
750  Real bc_forcing_x = zero;
751  Real u_target = zero;
752 
753  // roof forcing
754  if (roof_mask == one) {
755  u1_cellaway = u(i, j, k+1) ;
756  u2_cellaway = fourth * ( v(i, j , k+1) + v(i-1, j , k+1)
757  + v(i, j+1, k+1) + v(i-1, j+1, k+1) ) ;
758  rho_xface_inside = myhalf * (cell_data(i,j,k-1,Rho_comp) + cell_data(i-1,j,k-1,Rho_comp));
759  theta_surf = (myhalf * (cell_data(i,j,k-1,RhoTheta_comp) + cell_data(i-1,j,k-1, RhoTheta_comp))) / rho_xface_inside;
760  u_target = compute_if_most_target_vel(u1_cellaway, u2_cellaway, dx_z, z0, t_blank_law, theta_xface, theta_surf, tflux_in, Olen_in, l_stability_correction, stab_wind_floor, psi_cap_factor);
761  bc_forcing_x = -(u_target - ux); // BC forcing pushes nonrelative velocity toward target velocity
762  drag += bc_forcing_x * roof_mask * rho_xface * CdM * U_s;
763  }
764 
765  // south wall forcing
766  if (south_mask == one) {
767  u1_cellaway = u(i, j-1, k );
768  u2_cellaway = fourth * ( w(i, j-1, k ) + w(i-1, j-1, k )
769  + w(i, j-1, k+1) + w(i-1, j-1, k+1) ) ;
770  rho_xface_inside = myhalf * ( cell_data(i,j+1,k,Rho_comp) + cell_data(i-1,j+1,k,Rho_comp) );
771  theta_surf = (myhalf * (cell_data(i,j+1,k,RhoTheta_comp) + cell_data(i-1,j+1,k, RhoTheta_comp))) / rho_xface_inside;
772  u_target = compute_if_most_target_vel(u1_cellaway, u2_cellaway, dx_y, z0, t_blank_law, theta_xface, theta_surf, tflux_in, Olen_in, l_stability_correction, stab_wind_floor, psi_cap_factor);
773  bc_forcing_x = -(u_target - ux); // BC forcing pushes nonrelative velocity toward target velocity
774  drag += bc_forcing_x * south_mask * rho_xface * CdM * U_s;
775  }
776 
777  // north wall forcing
778  if (north_mask == one) {
779  u1_cellaway = u(i, j+1, k ) ;
780  u2_cellaway = fourth * ( w(i, j+1, k ) + w(i-1, j+1, k )
781  + w(i, j+1, k+1) + w(i-1, j+1, k+1) ) ;
782  rho_xface_inside = myhalf * ( cell_data(i,j-1,k,Rho_comp) + cell_data(i-1,j-1,k,Rho_comp) );
783  theta_surf = (myhalf * (cell_data(i,j-1,k,RhoTheta_comp) + cell_data(i-1,j-1,k, RhoTheta_comp))) / rho_xface_inside;
784  u_target = compute_if_most_target_vel(u1_cellaway, u2_cellaway, dx_y, z0, t_blank_law, theta_xface, theta_surf, tflux_in, Olen_in, l_stability_correction, stab_wind_floor, psi_cap_factor);
785  bc_forcing_x = -(u_target - ux); // BC forcing pushes nonrelative velocity toward target velocity
786  drag += bc_forcing_x * north_mask * rho_xface * CdM * U_s;
787  }
788 
789  // wall forcing (if not using most) or east/west walls when using MOST
790  if (wall_mask == one || east_west_mask == one) {
791  drag += (wall_mask + east_west_mask) * t_blank * rho_xface * CdM * ux * windspeed;
792  }
793 
794  // interior cell forcing
795  if (interior_mask == one) {
796  drag += interior_mask * rho_xface * CdM * ux * windspeed;
797  }
798 
799  if (l_implicit_drag) {
800  // point-implicit rescale of the aggregated drag
801  const Real lambda = CdM * ( (roof_mask + south_mask + north_mask) * U_s
802  + (wall_mask + east_west_mask) * t_blank * windspeed
803  + interior_mask * windspeed );
804  xmom_src_arr(i,j,k) -= drag / (one + lambda*dt);
805  } else if (is_slow_step && !use_ImmersedForcing_fast) {
806  // limit drag term for anelastic for numerical stability
807  Real d_drag = dt * -drag; // time step * acceleration like tendency
808  Real wsmax_change = damp_alpha * amrex::max(amrex::Math::abs(ux), ws_floor); // aims to prevent oscillations around 0.
809  if (amrex::Math::abs(ux) < 0.1){ // no damping for smaller velocities
810  wsmax_change =one * amrex::max(amrex::Math::abs(ux), ws_floor);
811  }
812  d_drag = amrex::min(amrex::max(d_drag, -wsmax_change), wsmax_change);
813  xmom_src_arr(i,j,k) += d_drag / dt; // put back as limited tendency
814  } else {
815  xmom_src_arr(i, j, k) -= drag;
816  }
817  });
818 }
compute_if_most_target_vel
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real compute_if_most_target_vel(const amrex::Real u1_2r, const amrex::Real u2_2r, const amrex::Real delta, const amrex::Real z0, const amrex::Real t_blank, const amrex::Real theta_face, const amrex::Real theta_surf, const amrex::Real tflux_in, const amrex::Real Olen_in, const bool stability_correction, const amrex::Real wind_floor, const amrex::Real psi_cap_factor)
Definition: ERF_ImmersedForcing.cpp:31
fourth
constexpr amrex::Real fourth
Definition: ERF_NumericalConstants.H:42
MORRInd::fc
@ fc
Definition: ERF_AdvanceMorrison.cpp:141
if_bld::horiz_masks
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE HorizMasks horiz_masks(amrex::Real t, amrex::Real below, amrex::Real above, amrex::Real lo, amrex::Real hi, bool normal_face, bool use_most, bool snap) noexcept
Definition: ERF_ImmersedBuildingMasks.H:140
SolverChoice::if_psi_cap_factor
amrex::Real if_psi_cap_factor
Cap on psi_m in the immersed momentum wall law as a fraction of ln(z/z0) (1: none).
Definition: ERF_DataStruct.H:2722
SolverChoice::if_Cd_momentum
amrex::Real if_Cd_momentum
Immersed-forcing drag coefficient for momentum.
Definition: ERF_DataStruct.H:2699
SolverChoice::if_Olen_in
amrex::Real if_Olen_in
Input Obukhov length for immersed-forcing MOST [m].
Definition: ERF_DataStruct.H:2710
SolverChoice::if_small_volfrac
amrex::Real if_small_volfrac
Immersed-forcing solid fraction below which a cell or face counts as fluid; read from eb2....
Definition: ERF_DataStruct.H:2723
SolverChoice::if_implicit_projection
bool if_implicit_projection
Definition: ERF_DataStruct.H:2704
SolverChoice::if_damp_alpha
amrex::Real if_damp_alpha
Immersed-forcing damping coefficient.
Definition: ERF_DataStruct.H:2720
SolverChoice::if_ws_floor
amrex::Real if_ws_floor
Wind-speed floor [m/s] in the per-step change limit of the building wall-law forcing.
Definition: ERF_DataStruct.H:2719
SolverChoice::immersed_forcing_substep
bool immersed_forcing_substep
Whether immersed-forcing source terms are applied only during substeps.
Definition: ERF_DataStruct.H:2695
SolverChoice::if_use_most
bool if_use_most
Whether immersed-forcing MOST is enabled.
Definition: ERF_DataStruct.H:2711
SolverChoice::if_implicit_drag
bool if_implicit_drag
Definition: ERF_DataStruct.H:2703
SolverChoice::if_stability_correction
bool if_stability_correction
Whether immersed-forcing stability corrections are enabled.
Definition: ERF_DataStruct.H:2713
if_bld::HorizFace
Definition: ERF_ImmersedBuildingMasks.H:60
if_bld::HorizMasks
Definition: ERF_ImmersedBuildingMasks.H:119
if_bld::HorizMasks::lo
amrex::Real lo
Definition: ERF_ImmersedBuildingMasks.H:120
if_bld::HorizMasks::hi
amrex::Real hi
Definition: ERF_ImmersedBuildingMasks.H:120
if_bld::HorizMasks::roof
amrex::Real roof
Definition: ERF_ImmersedBuildingMasks.H:120
if_bld::HorizMasks::wall
amrex::Real wall
Definition: ERF_ImmersedBuildingMasks.H:120
if_bld::HorizMasks::side
amrex::Real side
Definition: ERF_ImmersedBuildingMasks.H:120
if_bld::HorizMasks::interior
amrex::Real interior
Definition: ERF_ImmersedBuildingMasks.H:120

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

835 {
836  // geometric properties
837  const auto dx_arr = geom.CellSizeArray();
838  const Real dx_x = dx_arr[0];
839  const Real dx_y = dx_arr[1];
840  const Real dt = fac;
841 
842  const Real alpha_m = solverChoice.if_Cd_momentum;
843  const Real tiny = std::numeric_limits<amrex::Real>::epsilon();
844  const Real U_s = one; // unit velocity scale
845 
846  // MOST parameters
847  const Real z0 = solverChoice.if_z0;
848  const Real tflux_in = solverChoice.if_surf_temp_flux;
849  const Real Olen_in = solverChoice.if_Olen_in;
850  const Real stab_wind_floor = solverChoice.if_stability_wind_floor;
851  const Real psi_cap_factor = solverChoice.if_psi_cap_factor;
852  const bool l_use_most = solverChoice.if_use_most;
853  const bool l_stability_correction = solverChoice.if_stability_correction;
854 
855  // To limit stiffness of drag when using anelastic
856  const Real ws_floor = solverChoice.if_ws_floor;
857  const Real damp_alpha = solverChoice.if_damp_alpha;
858  // Point-implicit alternative to the clamp above; stabilizes both compressible and anelastic
859  const bool l_implicit_drag = solverChoice.if_implicit_drag;
860  // erf.if_implicit_projection: the linear drag of the solid is applied inside the anelastic
861  // projection (ERF::make_if_projection_sigma), not here
862  const bool l_ip = solverChoice.if_implicit_projection;
863 
864  const bool is_slow_step = true; // This is determined by calling context
865  const bool use_ImmersedForcing_fast = solverChoice.immersed_forcing_substep;
866  const Real small_volfrac = solverChoice.if_small_volfrac;
867  // erf.if_snap_partial_cells: read the cell blanking snapped to solid (1)
868  // or fluid (0) at half, so a height-map building becomes the same
869  // staircase of whole cells an exact box is. A face is solid when either
870  // cell it joins is; a face between a solid and a fluid cell is
871  // wall-normal and gets the interior drag toward zero (no penetration); a
872  // face between two solid cells carries the roof or wall law of its row
873  // (the full log-law target, not the partial-cell weighted one) or the
874  // interior drag, never both; the partial-cell branches (wall_mask,
875  // east_west_mask and the like, faces with 0 < t_blank < 1) do not arise
876  // under the snap. Off (the default), the raw fractions
877  // are used and nothing below changes: a boundary face of an exact box
878  // carries the wall law and the interior drag together as before. See
879  // SolverChoice::if_snap_partial_cells.
880  const bool l_snap = solverChoice.if_snap_partial_cells;
881  // The face blanking and its neighbours are read through ERF_ImmersedBuildingMasks.H, as the
882  // implicit projection reads them
883  const int dom_klo = geom.Domain().smallEnd(2);
884 
885 
886  ParallelFor(tby, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
887  {
888  const Real ux = fourth * ( u(i , j , k ) + u(i , j-1, k )
889  + u(i+1, j , k ) + u(i+1, j-1, k ) );
890  const Real uy = v(i, j, k);
891  const Real uz = fourth * ( w(i , j , k ) + w(i , j-1, k )
892  + w(i , j , k+1) + w(i , j-1, k+1) );
893  const amrex::Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz);
894 
895  const Real rho_yface = myhalf * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j-1,k,Rho_comp) );
896  const Real theta_yface = (myhalf * (cell_data(i,j,k ,RhoTheta_comp) + cell_data(i,j-1,k,RhoTheta_comp))) / rho_yface;
897 
898  // the face, its neighbours and the raw blanking the explicit path zeroes on (ERF_ImmersedBuildingMasks.H)
899  const if_bld::HorizFace fc = if_bld::horiz_face<1>(t_blank_arr, t_blank_yface_arr, i, j, k, dom_klo, l_snap, small_volfrac);
900  const Real t_blank = fc.t;
901  // With the snap the boundary solid face stands for the wall layer and
902  // takes the full log-law target; the partial-cell weight (1 - t_blank)
903  // of compute_if_most_target_vel() would make it zero (a no-slip
904  // staircase), so the snap path passes a zero blanking to it.
905  const Real t_blank_law = l_snap ? zero : t_blank;
906  const Real t_blank_below = fc.below;
907  const Real t_blank_above = fc.above;
908  const Real t_blank_west = fc.lo;
909  const Real t_blank_east = fc.hi;
910 
911  const Real dx_z = (z_cc_arr) ? (z_cc_arr(i,j,k) - z_cc_arr(i,j,k-1)) : dx_arr[2];
912  const Real drag_coefficient = alpha_m / std::pow(dx_x*dx_y*dx_z, one/three);
913  const Real CdM = std::min(drag_coefficient / (windspeed + tiny), drag_coefficient);
914  // With the snap a face joining a solid and a fluid cell is wall-normal:
915  // it gets the interior drag (no penetration), not a wall law, whatever
916  // row it lies in.
917  const bool normal_face = fc.normal;
918 
919  // As in the x-momentum: with the snap a roof face lies in the top
920  // solid row and a wall-law face is not also an interior face.
921  // roof, west/east wall law, partial-wall and interior drag (ERF_ImmersedBuildingMasks.H)
922  if_bld::HorizMasks mk = if_bld::horiz_masks(t_blank, t_blank_below, t_blank_above, t_blank_west, t_blank_east,
923  normal_face, l_use_most, l_snap);
924  // erf.if_implicit_projection: as in the x-momentum
925  if (l_ip) {
926  mk.wall = zero; mk.side = zero; mk.interior = zero;
927  if (fc.raw == one) { mk.roof = zero; mk.lo = zero; mk.hi = zero; }
928  }
929  const Real roof_mask = mk.roof;
930  const Real west_mask = mk.lo;
931  const Real east_mask = mk.hi;
932  const Real wall_mask = mk.wall;
933  const Real north_south_mask = mk.side;
934  const Real interior_mask = mk.interior;
935 
936  Real drag = zero;
937  Real u1_cellaway = zero;
938  Real u2_cellaway = zero;
939  Real rho_yface_inside = rho_yface;
940  Real theta_surf = theta_yface;
941  Real bc_forcing_y = zero;
942  Real u_target = zero;
943 
944  // roof forcing
945  if (roof_mask == one) {
946  u1_cellaway = fourth * ( u(i , j , k+1) + u(i , j-1, k+1)
947  + u(i+1, j , k+1) + u(i+1, j-1, k+1) );
948  u2_cellaway = v(i, j, k+1);
949  rho_yface_inside = myhalf * ( cell_data(i,j,k-1,Rho_comp) + cell_data(i,j-1,k-1,Rho_comp) );
950  theta_surf = (myhalf * (cell_data(i,j,k-1,RhoTheta_comp) + cell_data(i,j-1,k-1,RhoTheta_comp))) / rho_yface_inside;
951  u_target = compute_if_most_target_vel(u1_cellaway, u2_cellaway, dx_z, z0, t_blank_law, theta_yface, theta_surf, tflux_in, Olen_in, l_stability_correction, stab_wind_floor, psi_cap_factor);
952  bc_forcing_y = -(u_target - uy); // BC forcing pushes nonrelative velocity toward target velocity
953  drag += bc_forcing_y * roof_mask * rho_yface * CdM * U_s;
954  }
955 
956  // west wall forcing
957  if (west_mask == one) {
958  u1_cellaway = v(i-1, j , k );
959  u2_cellaway = fourth * ( w(i-1, j , k ) + w(i-1, j-1, k )
960  + w(i-1, j , k+1) + w(i-1, j-1, k+1) );
961  rho_yface_inside = myhalf * ( cell_data(i+1,j,k,Rho_comp) + cell_data(i+1,j-1,k,Rho_comp) );
962  theta_surf = (myhalf * (cell_data(i+1,j,k,RhoTheta_comp) + cell_data(i+1,j-1,k,RhoTheta_comp))) / rho_yface_inside;
963  u_target = compute_if_most_target_vel(u1_cellaway, u2_cellaway, dx_x, z0, t_blank_law, theta_yface, theta_surf, tflux_in, Olen_in, l_stability_correction, stab_wind_floor, psi_cap_factor);
964  bc_forcing_y = -(u_target - uy); // BC forcing pushes nonrelative velocity toward target velocity
965  drag += bc_forcing_y * west_mask * rho_yface * CdM * U_s;
966  }
967 
968  // east wall forcing
969  if (east_mask == one) {
970  u1_cellaway = v(i+1, j , k );
971  u2_cellaway = fourth * ( w(i+1, j , k ) + w(i+1, j-1, k )
972  + w(i+1, j , k+1) + w(i+1, j-1, k+1) );
973  rho_yface_inside = myhalf * ( cell_data(i-1,j,k,Rho_comp) + cell_data(i-1,j-1,k,Rho_comp) );
974  theta_surf = (myhalf * (cell_data(i-1,j,k,RhoTheta_comp) + cell_data(i-1,j-1,k,RhoTheta_comp))) / rho_yface_inside;
975  u_target = compute_if_most_target_vel(u1_cellaway, u2_cellaway, dx_x, z0, t_blank_law, theta_yface, theta_surf, tflux_in, Olen_in, l_stability_correction, stab_wind_floor, psi_cap_factor);
976  bc_forcing_y = -(u_target - uy); // BC forcing pushes nonrelative velocity toward target velocity
977  drag += bc_forcing_y * east_mask * rho_yface * CdM * U_s;
978  }
979 
980  // wall forcing (if not using most) or north/south walls when using MOST
981  if (wall_mask == one || north_south_mask == one) {
982  drag += (wall_mask + north_south_mask) * t_blank * rho_yface * CdM * uy * windspeed;
983  }
984 
985  // interior cell forcing
986  if (interior_mask == one) {
987  drag += interior_mask * rho_yface * CdM * uy * windspeed;
988  }
989 
990  if (l_implicit_drag) {
991  // point-implicit rescale of the aggregated drag
992  const Real lambda = CdM * ( (roof_mask + west_mask + east_mask) * U_s
993  + (wall_mask + north_south_mask) * t_blank * windspeed
994  + interior_mask * windspeed );
995  ymom_src_arr(i,j,k) -= drag / (one + lambda*dt);
996  } else if (is_slow_step && !use_ImmersedForcing_fast) {
997  // limit drag term for anelastic for numerical stability
998  Real d_drag = dt * -drag; // time step * acceleration like tendency
999  Real wsmax_change = damp_alpha * amrex::max(amrex::Math::abs(uy), ws_floor); // aims to prevent oscillations around 0.
1000  if (amrex::Math::abs(uy) < 0.1){ // no damping for smaller velocities
1001  wsmax_change =one * amrex::max(amrex::Math::abs(uy), ws_floor);
1002  }
1003  d_drag = amrex::min(amrex::max(d_drag, -wsmax_change), wsmax_change);
1004  ymom_src_arr(i,j,k) += d_drag / dt; // put back as limited tendency
1005  } else {
1006  ymom_src_arr(i, j, k) -= drag;
1007  }
1008  });
1009 }

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

1026 {
1027  // geometric properties
1028  const auto dx_arr = geom.CellSizeArray();
1029  const Real dx_x = dx_arr[0];
1030  const Real dx_y = dx_arr[1];
1031  const Real dt = fac;
1032 
1033  const Real alpha_m = solverChoice.if_Cd_momentum;
1034  const Real tiny = std::numeric_limits<amrex::Real>::epsilon();
1035  const Real U_s = one; // unit velocity scale
1036 
1037  // MOST parameters
1038  const Real z0 = solverChoice.if_z0;
1039  const Real tflux_in = solverChoice.if_surf_temp_flux;
1040  const Real Olen_in = solverChoice.if_Olen_in;
1041  const Real stab_wind_floor = solverChoice.if_stability_wind_floor;
1042  const Real psi_cap_factor = solverChoice.if_psi_cap_factor;
1043  const bool l_use_most = solverChoice.if_use_most;
1044  const bool l_stability_correction = solverChoice.if_stability_correction;
1045 
1046  // To limit stiffness of drag when using anelastic
1047  const Real ws_floor = solverChoice.if_ws_floor;
1048  const Real damp_alpha = solverChoice.if_damp_alpha;
1049  // Point-implicit alternative to the clamp above; stabilizes both compressible and anelastic
1050  const bool l_implicit_drag = solverChoice.if_implicit_drag;
1051  // erf.if_implicit_projection: the linear drag of the solid is applied inside the anelastic
1052  // projection (ERF::make_if_projection_sigma), not here
1053  const bool l_ip = solverChoice.if_implicit_projection;
1054 
1055  const bool is_slow_step = true; // This is determined by calling context
1056  const bool use_ImmersedForcing_fast = solverChoice.immersed_forcing_substep;
1057  const Real small_volfrac = solverChoice.if_small_volfrac;
1058  // erf.if_snap_partial_cells: read the cell blanking snapped to solid (1)
1059  // or fluid (0) at half, so a height-map building becomes the same
1060  // staircase of whole cells an exact box is. A face is solid when either
1061  // cell it joins is; a face between a solid and a fluid cell is
1062  // wall-normal and gets the interior drag toward zero (no penetration); a
1063  // face between two solid cells carries the roof or wall law of its row
1064  // (the full log-law target, not the partial-cell weighted one) or the
1065  // interior drag, never both; the partial-cell branches (wall_mask,
1066  // east_west_mask and the like, faces with 0 < t_blank < 1) do not arise
1067  // under the snap. Off (the default), the raw fractions
1068  // are used and nothing below changes: a boundary face of an exact box
1069  // carries the wall law and the interior drag together as before. See
1070  // SolverChoice::if_snap_partial_cells.
1071  const bool l_snap = solverChoice.if_snap_partial_cells;
1072  // The face blanking and its neighbours are read through ERF_ImmersedBuildingMasks.H, as the
1073  // implicit projection reads them
1074 
1075 
1076  ParallelFor(tbz, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
1077  {
1078  const Real ux = fourth * ( u(i , j , k ) + u(i+1, j , k )
1079  + u(i , j , k-1) + u(i+1, j , k-1) );
1080  const Real uy = fourth * ( v(i, j , k ) + v(i, j+1, k )
1081  + v(i, j , k-1) + v(i, j+1, k-1) );
1082  const Real uz = w(i, j, k);
1083  const amrex::Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz);
1084 
1085  const Real rho_zface = myhalf * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j,k-1,Rho_comp) );
1086  const Real theta_zface = (myhalf * (cell_data(i,j,k,RhoTheta_comp) + cell_data(i,j,k-1,RhoTheta_comp))) / rho_zface;
1087 
1088  // the face, its neighbours and the raw blanking the explicit path zeroes on (ERF_ImmersedBuildingMasks.H)
1089  const if_bld::VertFace fc = if_bld::vert_face(t_blank_arr, t_blank_zface_arr, i, j, k, l_snap, small_volfrac);
1090  const Real t_blank = fc.t;
1091  // With the snap the boundary solid face stands for the wall layer and
1092  // takes the full log-law target; the partial-cell weight (1 - t_blank)
1093  // of compute_if_most_target_vel() would make it zero (a no-slip
1094  // staircase), so the snap path passes a zero blanking to it.
1095  const Real t_blank_law = l_snap ? zero : t_blank;
1096  const Real t_blank_above = fc.above;
1097  const Real t_blank_south = fc.south;
1098  const Real t_blank_north = fc.north;
1099  const Real t_blank_west = fc.west;
1100  const Real t_blank_east = fc.east;
1101 
1102  const Real dx_z = (z_cc_arr) ? (z_cc_arr(i,j,k) - z_cc_arr(i,j,k-1)) : dx_arr[2];
1103  const Real drag_coefficient = alpha_m / std::pow(dx_x*dx_y*dx_z, one/three);
1104  const Real CdM = std::min(drag_coefficient / (windspeed + tiny), drag_coefficient);
1105  // With the snap a face joining a solid and a fluid cell is wall-normal:
1106  // it gets the interior drag (no penetration), not a wall law, whatever
1107  // row it lies in.
1108  const bool normal_face = fc.normal;
1109 
1110  // wall law on the four sides, partial-wall, roof and interior drag (ERF_ImmersedBuildingMasks.H);
1111  // with the snap a face carrying a wall law or the roof drag is not also an interior face
1112  if_bld::VertMasks mk = if_bld::vert_masks(t_blank, t_blank_above, t_blank_south, t_blank_north,
1113  t_blank_west, t_blank_east, k >= 1, normal_face, l_use_most, l_snap);
1114  // erf.if_implicit_projection: as in the x-momentum
1115  if (l_ip) {
1116  mk.wall = zero; mk.roof = zero; mk.interior = zero;
1117  if (fc.raw == one) { mk.south = zero; mk.north = zero; mk.west = zero; mk.east = zero; }
1118  }
1119  const Real south_mask = mk.south;
1120  const Real north_mask = mk.north;
1121  const Real west_mask = mk.west;
1122  const Real east_mask = mk.east;
1123  const Real wall_mask = mk.wall;
1124  const Real roof_mask = mk.roof;
1125  const Real interior_mask = mk.interior;
1126 
1127  Real drag = zero;
1128  Real u1_cellaway = zero;
1129  Real u2_cellaway = zero;
1130  Real rho_zface_inside = rho_zface;
1131  Real theta_surf = theta_zface;
1132  Real bc_forcing_z = zero;
1133  Real u_target = zero;
1134 
1135  // south wall forcing
1136  if (south_mask == one) {
1137  u1_cellaway = fourth * ( u(i , j-1, k ) + u(i+1, j-1, k )
1138  + u(i , j-1, k-1) + u(i+1, j-1, k-1) );
1139  u2_cellaway = w(i, j-1, k);
1140  rho_zface_inside = myhalf * ( cell_data(i,j+1,k,Rho_comp) + cell_data(i,j+1,k-1,Rho_comp) );
1141  theta_surf = (myhalf * (cell_data(i,j+1,k,RhoTheta_comp) + cell_data(i,j+1,k-1,RhoTheta_comp))) / rho_zface_inside;
1142  u_target = compute_if_most_target_vel(u1_cellaway, u2_cellaway, dx_y, z0, t_blank_law, theta_zface, theta_surf, tflux_in, Olen_in, l_stability_correction, stab_wind_floor, psi_cap_factor);
1143  bc_forcing_z = -(u_target - uz); // BC forcing pushes nonrelative velocity toward target velocity
1144  drag += bc_forcing_z * south_mask * rho_zface * CdM * U_s;
1145  }
1146 
1147  // north wall forcing
1148  if (north_mask == one) {
1149  u1_cellaway = fourth * ( u(i , j+1, k ) + u(i+1, j+1, k )
1150  + u(i , j+1, k-1) + u(i+1, j+1, k-1) );
1151  u2_cellaway = w(i, j+1, k);
1152  rho_zface_inside = myhalf * ( cell_data(i,j-1,k,Rho_comp) + cell_data(i,j-1,k-1,Rho_comp) );
1153  theta_surf = (myhalf * (cell_data(i,j-1,k,RhoTheta_comp) + cell_data(i,j-1,k-1,RhoTheta_comp))) / rho_zface_inside;
1154  u_target = compute_if_most_target_vel(u1_cellaway, u2_cellaway, dx_y, z0, t_blank_law, theta_zface, theta_surf, tflux_in, Olen_in, l_stability_correction, stab_wind_floor, psi_cap_factor);
1155  bc_forcing_z = -(u_target - uz); // BC forcing pushes nonrelative velocity toward target velocity
1156  drag += bc_forcing_z * north_mask * rho_zface * CdM * U_s;
1157  }
1158 
1159  // west wall forcing
1160  if (west_mask == one) {
1161  u1_cellaway = fourth * ( v(i-1, j , k ) + v(i-1, j+1, k )
1162  + v(i-1, j , k-1) + v(i-1, j+1, k-1) );
1163  u2_cellaway = w(i-1, j, k);
1164  rho_zface_inside = myhalf * ( cell_data(i+1,j,k,Rho_comp) + cell_data(i+1,j,k-1,Rho_comp) );
1165  theta_surf = (myhalf * (cell_data(i+1,j,k,RhoTheta_comp) + cell_data(i+1,j,k-1,RhoTheta_comp))) / rho_zface_inside;
1166  u_target = compute_if_most_target_vel(u1_cellaway, u2_cellaway, dx_x, z0, t_blank_law, theta_zface, theta_surf, tflux_in, Olen_in, l_stability_correction, stab_wind_floor, psi_cap_factor);
1167  bc_forcing_z = -(u_target - uz); // BC forcing pushes nonrelative velocity toward target velocity
1168  drag += bc_forcing_z * west_mask * rho_zface * CdM * U_s;
1169  }
1170 
1171  // east wall forcing
1172  if (east_mask == one) {
1173  u1_cellaway = fourth * ( v(i+1, j , k ) + v(i+1, j+1, k )
1174  + v(i+1, j , k-1) + v(i+1, j+1, k-1) );
1175  u2_cellaway = w(i+1, j, k);
1176  rho_zface_inside = myhalf * ( cell_data(i-1,j,k,Rho_comp) + cell_data(i-1,j,k-1,Rho_comp) );
1177  theta_surf = (myhalf * (cell_data(i-1,j,k,RhoTheta_comp) + cell_data(i-1,j,k-1,RhoTheta_comp))) / rho_zface_inside;
1178  u_target = compute_if_most_target_vel(u1_cellaway, u2_cellaway, dx_x, z0, t_blank_law, theta_zface, theta_surf, tflux_in, Olen_in, l_stability_correction, stab_wind_floor, psi_cap_factor);
1179  bc_forcing_z = -(u_target - uz); // BC forcing pushes nonrelative velocity toward target velocity
1180  drag += bc_forcing_z * east_mask * rho_zface * CdM * U_s;
1181  }
1182 
1183  // wall forcing (if not using most) or roof when using MOST
1184  if (wall_mask == one || roof_mask == one) {
1185  drag += (wall_mask + roof_mask) * t_blank * rho_zface * CdM * uz * windspeed;
1186  }
1187 
1188  // interior cell forcing
1189  if (interior_mask == one) {
1190  drag += interior_mask * rho_zface * CdM * uz * windspeed;
1191  }
1192 
1193  if (l_implicit_drag) {
1194  // point-implicit rescale of the aggregated drag
1195  const Real lambda = CdM * ( (south_mask + north_mask + west_mask + east_mask) * U_s
1196  + (wall_mask + roof_mask) * t_blank * windspeed
1197  + interior_mask * windspeed );
1198  zmom_src_arr(i,j,k) -= drag / (one + lambda*dt);
1199  } else if (is_slow_step && !use_ImmersedForcing_fast) {
1200  // limit drag term for anelastic for numerical stability
1201  Real d_drag = dt * -drag; // time step * acceleration like tendency
1202  Real wsmax_change = damp_alpha * amrex::max(amrex::Math::abs(uz), ws_floor); // aims to prevent oscillations around 0.
1203  if (amrex::Math::abs(uz) < 0.1){ // no damping for smaller velocities
1204  wsmax_change = one * amrex::max(amrex::Math::abs(uz), ws_floor);
1205  }
1206  d_drag = amrex::min(amrex::max(d_drag, -wsmax_change), wsmax_change);
1207  zmom_src_arr(i,j,k) += d_drag / dt; // put back as limited tendency
1208  } else {
1209  zmom_src_arr(i, j, k) -= drag;
1210  }
1211  });
1212 }
if_bld::vert_masks
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE VertMasks vert_masks(amrex::Real t, amrex::Real above, amrex::Real south, amrex::Real north, amrex::Real west, amrex::Real east, bool above_floor, bool normal_face, bool use_most, bool snap) noexcept
Definition: ERF_ImmersedBuildingMasks.H:185
if_bld::vert_face
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE VertFace vert_face(const amrex::Array4< const amrex::Real > &t, const amrex::Array4< const amrex::Real > &tf, int i, int j, int k, bool snap, amrex::Real small) noexcept
Definition: ERF_ImmersedBuildingMasks.H:100
if_bld::VertFace
Definition: ERF_ImmersedBuildingMasks.H:94
if_bld::VertMasks
Definition: ERF_ImmersedBuildingMasks.H:167
if_bld::VertMasks::interior
amrex::Real interior
Definition: ERF_ImmersedBuildingMasks.H:168
if_bld::VertMasks::east
amrex::Real east
Definition: ERF_ImmersedBuildingMasks.H:168
if_bld::VertMasks::north
amrex::Real north
Definition: ERF_ImmersedBuildingMasks.H:168
if_bld::VertMasks::wall
amrex::Real wall
Definition: ERF_ImmersedBuildingMasks.H:168
if_bld::VertMasks::roof
amrex::Real roof
Definition: ERF_ImmersedBuildingMasks.H:168
if_bld::VertMasks::south
amrex::Real south
Definition: ERF_ImmersedBuildingMasks.H:168
if_bld::VertMasks::west
amrex::Real west
Definition: ERF_ImmersedBuildingMasks.H:168

Referenced by make_mom_sources().

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

◆ ImmersedForcingTerrain_HorizMom_FractionStress()

template<int dir>
void ImmersedForcingTerrain_HorizMom_FractionStress ( const Box &  tb,
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_face_arr,
const Array4< Real > &  mom_src_arr,
const Geometry &  geom,
const SolverChoice &  solverChoice,
const Real  dt,
const Real  time 
)

Fraction-stress wall law (erf.if_wall_form = fraction_stress) for a horizontal momentum component, see ERF_ImmersedWallCell.H: the wall stress in the wall cell of each face column, the drag of the terrain kernels elsewhere (solid cells, and partial cells under the wall cell). dir = 0 forces x-momentum, dir = 1 y-momentum.

404 {
405  const Real* dx_arr = geom.CellSize();
406  const Real dx_x = dx_arr[0];
407  const Real dx_y = dx_arr[1];
408  const Real dz = dx_arr[2]; // fraction_stress runs on a constant-dz mesh
409  const Real alpha_m = solverChoice.if_Cd_momentum;
410  const Real z0 = solverChoice.if_z0;
411  const bool l_implicit_drag = solverChoice.if_implicit_drag;
412  // erf.if_implicit_projection: the drag of the solid is applied inside the anelastic projection
413  const bool l_drag = !solverChoice.if_implicit_projection;
414  const Real tiny = std::numeric_limits<Real>::epsilon();
415  const Real small = solverChoice.if_small_volfrac;
416  const ib_wall::WallCond cond = make_wall_cond(solverChoice, time);
417  constexpr int io = (dir == 0) ? 1 : 0;
418  constexpr int jo = (dir == 1) ? 1 : 0;
419  const int klo = geom.Domain().smallEnd(2);
420  const int khi = geom.Domain().bigEnd(2);
421 
422  ParallelFor(tb, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
423  {
424  // solid fraction on this face column, as the terrain kernels read it, with the column
425  // ends clamped to the domain as make_ib_wall_face_masks does
426  const Real b_km1 = if_face_beta(t_blank_arr, t_blank_face_arr, i, j, k-1, io, jo, klo, khi, small);
427  const Real b_k = if_face_beta(t_blank_arr, t_blank_face_arr, i, j, k, io, jo, klo, khi, small);
428  const Real b_kp1 = if_face_beta(t_blank_arr, t_blank_face_arr, i, j, k+1, io, jo, klo, khi, small);
429  const Real rho_face = myhalf * ( cell_data(i,j,k,Rho_comp) + cell_data(i-io,j-jo,k,Rho_comp) );
430  const Real un = (dir == 0) ? u(i, j, k) : v(i, j, k);
431 
432  if (ib_wall::is_wall_cell(b_km1, b_k, b_kp1, small)) {
433  const Real ut = if_face_speed<dir>(u, v, i, j, k);
434  const Real ut_ref = if_face_speed<dir>(u, v, i, j, k+1);
435  const ib_wall::WallState ws = ib_wall::wall_state(ut, ut_ref,
436  if_face_theta(cell_data, i, j, k, io, jo),
437  if_face_theta(cell_data, i, j, k+1, io, jo),
438  ib_wall::wall_beta(b_k, small), dz, z0, cond);
439  mom_src_arr(i, j, k) -= ib_wall::stress_rate(ws.ustar, ut, dz, dt) * rho_face * un;
440  } else if (l_drag) {
441  const Real wz = (dir == 0) ? fourth * ( w(i, j, k) + w(i-1, j, k) + w(i, j, k+1) + w(i-1, j, k+1) )
442  : fourth * ( w(i, j, k) + w(i, j-1, k) + w(i, j, k+1) + w(i, j-1, k+1) );
443  const Real uh = if_face_speed<dir>(u, v, i, j, k);
444  const Real windspeed = std::sqrt(uh * uh + wz * wz);
445  const Real drag_coefficient = alpha_m / std::pow(dx_x*dx_y*dz, one/three);
446  const Real CdM = std::min(drag_coefficient / (windspeed + tiny), drag_coefficient);
447  const Real lambda = b_k * CdM * windspeed; // linear drag rate [1/s]; zero in fluid cells
448  const Real fac_local = l_implicit_drag ? lambda / (one + lambda*dt) : lambda;
449  mom_src_arr(i, j, k) -= fac_local * rho_face * un;
450  }
451  });
452 }
make_wall_cond
static ib_wall::WallCond make_wall_cond(const SolverChoice &sc, const Real time)
Definition: ERF_ImmersedForcing.cpp:328
khi
const int khi
Definition: ERF_InitCustomPert_Bubble.H:25
klo
const int klo
Definition: ERF_InitCustomPert_KE.H:9
WDM6SedCellScratch::dz
@ dz
Definition: ERF_AdvanceWDM6.cpp:273
ib_wall::is_wall_cell
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool is_wall_cell(amrex::Real b_km1, amrex::Real b_k, amrex::Real b_kp1, amrex::Real small) noexcept
Definition: ERF_ImmersedWallCell.H:57
ib_wall::wall_beta
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real wall_beta(amrex::Real b_k, amrex::Real small) noexcept
Definition: ERF_ImmersedWallCell.H:65
ib_wall::wall_state
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE WallState wall_state(amrex::Real ut, amrex::Real ut_ref, amrex::Real theta, amrex::Real theta_ref, amrex::Real beta, amrex::Real dz, amrex::Real z0, const WallCond &c) noexcept
Definition: ERF_ImmersedWallCell.H:161
ib_wall::stress_rate
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real stress_rate(amrex::Real ustar, amrex::Real ut, amrex::Real dz, amrex::Real dt) noexcept
Definition: ERF_ImmersedWallCell.H:234
ib_wall::WallCond
Definition: ERF_ImmersedWallCell.H:121
ib_wall::WallState
Definition: ERF_ImmersedWallCell.H:134
ib_wall::WallState::ustar
amrex::Real ustar
Definition: ERF_ImmersedWallCell.H:135
Here is the call 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)

1229 {
1230  // geometric properties
1231  const auto dx_arr = geom.CellSizeArray();
1232  const Real dx_x = dx_arr[0];
1233  const Real dx_y = dx_arr[1];
1234 
1235  const Real alpha_h = solverChoice.if_Cd_scalar;
1236  // The terrain kernels keep the raw fractions: erf.if_snap_partial_cells
1237  // applies to the buildings kernels only (the wall law below is weighted
1238  // by the fluid fraction of the face, which a snapped face has not).
1239  const Real tiny = std::numeric_limits<amrex::Real>::epsilon();
1240  const Real U_s = one; // unit velocity scale
1241 
1242  // MOST parameters
1243  similarity_funs sfuns;
1244  const Real ggg = CONST_GRAV;
1245  const Real kappa = KAPPA;
1246  const Real z0 = solverChoice.if_z0;
1247  const Real tflux = solverChoice.if_surf_temp_flux;
1248  const Real init_surf_temp = solverChoice.if_init_surf_temp;
1249 
1250  // Note this has been converted to K / s when it was read in;
1251  const Real surf_heating_rate = solverChoice.if_surf_heating_rate;
1252 
1253  const Real Olen_in = solverChoice.if_Olen_in;
1254  const Real stab_wind_floor = solverChoice.if_stability_wind_floor;
1255 
1256  ParallelFor(bx, [=]
1257  AMREX_GPU_DEVICE(int i, int j, int k) noexcept
1258  {
1259  const Real dx_z = (z_cc_arr) ? (z_cc_arr(i,j,k) - z_cc_arr(i,j,k-1)) : dx_arr[2];
1260  const Real drag_coefficient = alpha_h / std::pow(dx_x*dx_y*dx_z, one/three);
1261 
1262  const Real t_blank = t_blank_arr(i, j, k);
1263  const Real t_blank_above = t_blank_arr(i, j, k+1);
1264  const Real ux_cc_2r = myhalf * (u(i ,j ,k+1) + u(i+1,j ,k+1));
1265  const Real uy_cc_2r = myhalf * (v(i ,j ,k+1) + v(i ,j+1,k+1));
1266  // The heat transfer uses the wind floored at erf.if_stability_wind_floor and an Obukhov
1267  // length bounded so that |zeta| <= ib_stability::zeta_max(), as on flat ground
1268  const Real h_windspeed2r = ib_stability::floored_wind(std::sqrt(ux_cc_2r * ux_cc_2r + uy_cc_2r * uy_cc_2r), stab_wind_floor);
1269 
1270  const Real theta = cell_data(i,j,k ,RhoTheta_comp) / cell_data(i,j,k ,Rho_comp);
1271  const Real theta_neighbor = cell_data(i,j,k+1,RhoTheta_comp) / cell_data(i,j,k+1,Rho_comp);
1272 
1273  // SURFACE TEMP AND HEATING/COOLING RATE
1274  if (init_surf_temp > zero) {
1275  if (t_blank > 0 && (t_blank_above == zero)) { // force to MOST value
1276  const Real surf_temp = init_surf_temp + surf_heating_rate*time;
1277  const Real bc_forcing_rt_srf = -(cell_data(i,j,k-1,Rho_comp) * surf_temp - cell_data(i,j,k-1,RhoTheta_comp));
1278  cell_src(i, j, k, RhoTheta_comp) -= drag_coefficient * U_s * bc_forcing_rt_srf;
1279  }
1280  }
1281 
1282  // SURFACE HEAT FLUX
1283  if (tflux != Real(1e-8)){
1284  if (t_blank > 0 && (t_blank_above == zero)) { // force to MOST value
1285  Real psi_m = zero;
1286  Real psi_h = zero;
1287  Real psi_h_neighbor = zero;
1288  Real ustar = h_windspeed2r * kappa / (std::log((Real(1.5)) * dx_z / z0) - psi_m);
1289  const Real Olen = ib_stability::bounded_obukhov_length(-ustar * ustar * ustar * theta / (kappa * ggg * tflux + tiny),
1290  Real(1.5) * dx_z);
1291  const Real zeta = ib_stability::bounded_zeta(myhalf * dx_z, Olen);
1292  const Real zeta_neighbor = ib_stability::bounded_zeta(Real(1.5) * dx_z, Olen);
1293 
1294  // similarity functions
1295  psi_m = sfuns.calc_psi_m(zeta);
1296  psi_h = sfuns.calc_psi_h(zeta);
1297  psi_h_neighbor = sfuns.calc_psi_h(zeta_neighbor);
1298  ustar = h_windspeed2r * kappa / (std::log((Real(1.5)) * dx_z / z0) - psi_m);
1299 
1300  // prevent some unphysical math
1301  ustar = ib_stability::clamped_ustar(ustar);
1302  // caps of ln(z / z0): erf.if_psi_cap_factor stays off the temperature forcing, whose
1303  // theta* comes from an L not made consistent with u* (the legacy law, kept as it is;
1304  // the fraction-stress law of erf-model/ERF#4134 forms them consistently)
1305  psi_h_neighbor = ib_stability::capped_psi(psi_h_neighbor, Real(1.5) * dx_z, z0, one);
1306  psi_h = ib_stability::capped_psi(psi_h, myhalf * dx_z, z0, one);
1307 
1308  // We do not know the actual temperature so use cell above
1309  const Real thetastar = theta * ustar * ustar / (kappa * ggg * Olen);
1310  const Real surf_temp = theta_neighbor - thetastar / kappa * (std::log((Real(1.5)) * dx_z / z0) - psi_h_neighbor);
1311  const Real tTarget = surf_temp + thetastar / kappa * (std::log((myhalf) * dx_z / z0) - psi_h);
1312 
1313  const Real bc_forcing_rt = -(cell_data(i,j,k,Rho_comp) * tTarget - cell_data(i,j,k,RhoTheta_comp));
1314  cell_src(i, j, k, RhoTheta_comp) -= drag_coefficient * U_s * bc_forcing_rt;
1315  }
1316  }
1317 
1318  // OBUKHOV LENGTH
1319  if (Olen_in != Real(1e-8)){
1320  if (t_blank > 0 && (t_blank_above == zero)) { // force to MOST value
1321  // the prescribed length as given, with its zeta held within +-zeta_max()
1322  const Real Olen = Olen_in;
1323  const Real zeta = ib_stability::bounded_zeta(myhalf * dx_z, Olen);
1324  const Real zeta_neighbor = ib_stability::bounded_zeta(Real(1.5) * dx_z, Olen);
1325 
1326  // similarity functions; psi_h capped at ln(z / z0) as on the other temperature branches
1327  const Real psi_m = sfuns.calc_psi_m(zeta);
1328  const Real psi_h = ib_stability::capped_psi(sfuns.calc_psi_h(zeta), myhalf * dx_z, z0, one);
1329  const Real psi_h_neighbor = ib_stability::capped_psi(sfuns.calc_psi_h(zeta_neighbor), Real(1.5) * dx_z, z0, one);
1330  // psi_m at 0.5 dz against ln(1.5 dz / z0) can cancel; u* is held in [0, 2] m/s as on the
1331  // other branches. In unstable air psi_h >= psi_m and psi_h grows with |zeta|, so wherever
1332  // psi_m reaches ln(1.5 dz / z0) both psi_h caps above bind, the two log brackets close
1333  // and the target is the cell above whatever u* is: the clamp only guards the arithmetic.
1334  const Real ustar = ib_stability::clamped_ustar(h_windspeed2r * kappa / (std::log((Real(1.5)) * dx_z / z0) - psi_m));
1335 
1336  // We do not know the actual temperature so use cell above
1337  const Real thetastar = theta * ustar * ustar / (kappa * ggg * Olen);
1338  const Real surf_temp = theta_neighbor - thetastar / kappa * (std::log((Real(1.5)) * dx_z / z0) - psi_h_neighbor);
1339  const Real tTarget = surf_temp + thetastar / kappa * (std::log((myhalf) * dx_z / z0) - psi_h);
1340 
1341  const Real bc_forcing_rt = -(cell_data(i,j,k,Rho_comp) * tTarget - cell_data(i,j,k,RhoTheta_comp));
1342  cell_src(i, j, k, RhoTheta_comp) -= drag_coefficient * U_s * bc_forcing_rt;
1343  }
1344  }
1345 
1346  // Force fully immersed cells to planar average rho and theta
1347  if (t_blank == one && r_avg && t_avg) {
1348  const Real rho_avg = r_avg(k);
1349  const Real theta_avg = t_avg(k) / rho_avg; // Convert from RhoTheta to Theta
1350  const Real rho_cell = cell_data(i,j,k,Rho_comp);
1351  const Real bc_forcing_r = -(rho_avg - rho_cell);
1352  const Real bc_forcing_rt = -(rho_avg * theta_avg - cell_data(i,j,k,RhoTheta_comp));
1353 
1354  cell_src(i, j, k, Rho_comp) -= drag_coefficient * U_s * bc_forcing_r;
1355  cell_src(i, j, k, RhoTheta_comp) -= drag_coefficient * U_s * bc_forcing_rt;
1356  }
1357  });
1358 }

Referenced by make_sources().

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

◆ ImmersedForcingTerrain_Scalar_FractionStress()

void ImmersedForcingTerrain_Scalar_FractionStress ( 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< Real > &  cell_src,
const Array4< Real > &  wall_hfx,
const Array4< Real > &  wall_ktarget,
const Geometry &  geom,
const SolverChoice &  solverChoice,
const Table1D< Real > &  r_avg,
const Table1D< Real > &  t_avg,
const Real  dt,
const Real  time,
const bool  wall_tke,
const Real  Cmu0 
)

Fraction-stress wall law for the scalars (erf.if_wall_form = fraction_stress), see ERF_ImmersedWallCell.H: the wall heat flux -u* theta* in the wall cell of each column, with the u* and theta* of the momentum law, and the relaxation of fully immersed cells to the planar average of the legacy kernel. The kinematic wall flux -u* theta* is stored in wall_hfx at the wall cell (zero elsewhere), for the TKE buoyancy on the wall face; with wall_tke the wall value of the TKE is stored in wall_ktarget at the wall cell (-1 elsewhere), and a time factor of 1 leaves the pin to erf_slow_rhs_post and the implicit TKE solve.

513 {
514  const Real* dx_arr = geom.CellSize();
515  const Real dz = dx_arr[2]; // fraction_stress runs on a constant-dz mesh
516  const Real tke_time_factor = solverChoice.if_wall_tke_time_factor;
517  // a time factor of 1 pins the wall cell at the wall value after every stage (erf_slow_rhs_post
518  // and the implicit TKE solve read wall_ktarget), as erf.dirichlet_k does on flat ground;
519  // larger factors relax toward it here as a source
520  const bool tke_pin = (tke_time_factor <= one);
521  const Real drag_coefficient = solverChoice.if_Cd_scalar / std::pow(dx_arr[0]*dx_arr[1]*dz, one/three);
522  const Real z0 = solverChoice.if_z0;
523  const Real small = solverChoice.if_small_volfrac;
524  const ib_wall::WallCond cond = make_wall_cond(solverChoice, time);
525  const int klo = geom.Domain().smallEnd(2);
526  const int khi = geom.Domain().bigEnd(2);
527  const bool have_avg = (r_avg && t_avg);
528 
529  ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
530  {
531  const Real b_k = if_cell_beta(t_blank_arr, i, j, k, klo, khi, small);
532  if (wall_hfx) { wall_hfx(i, j, k) = zero; }
533  if (wall_ktarget) { wall_ktarget(i, j, k) = -one; } // no wall value: not a wall cell
534 
535  if (ib_wall::is_wall_cell(if_cell_beta(t_blank_arr, i, j, k-1, klo, khi, small), b_k,
536  if_cell_beta(t_blank_arr, i, j, k+1, klo, khi, small), small)) {
537  // horizontal speed at the cell centres of the wall cell and the cell above
538  const Real uc0 = myhalf * (u(i, j, k) + u(i+1, j, k));
539  const Real vc0 = myhalf * (v(i, j, k) + v(i, j+1, k));
540  const Real uc1 = myhalf * (u(i, j, k+1) + u(i+1, j, k+1));
541  const Real vc1 = myhalf * (v(i, j, k+1) + v(i, j+1, k+1));
542  const Real rho = cell_data(i, j, k, Rho_comp);
543  const Real theta = cell_data(i, j, k, RhoTheta_comp) / rho;
544  const Real theta_ref = cell_data(i, j, k+1, RhoTheta_comp) / cell_data(i, j, k+1, Rho_comp);
545  const ib_wall::WallState ws = ib_wall::wall_state(std::sqrt(uc0*uc0 + vc0*vc0), std::sqrt(uc1*uc1 + vc1*vc1),
546  theta, theta_ref,
547  ib_wall::wall_beta(b_k, small), dz, z0, cond);
548  cell_src(i, j, k, RhoTheta_comp) += rho * ib_wall::heat_source(ws, theta, cond, dz, dt);
549  if (wall_hfx) { wall_hfx(i, j, k) = -ws.ustar * ws.thetastar; }
550  if (wall_tke) {
551  const Real k_wall = ib_wall::wall_tke(ws, ib_wall::wall_beta(b_k, small), dz, z0, theta, Cmu0);
552  if (wall_ktarget) { wall_ktarget(i, j, k) = k_wall; }
553  if (!tke_pin) {
554  // relax the wall-cell TKE toward its wall value over tke_time_factor steps
555  const Real k_now = cell_data(i, j, k, RhoKE_comp) / rho;
556  cell_src(i, j, k, RhoKE_comp) += rho * (k_wall - k_now) / (tke_time_factor * dt);
557  }
558  }
559  } else if (b_k == one && have_avg) {
560  // fully immersed cells: relax to the planar average, as the legacy kernel does
561  const Real rho_avg = r_avg(k);
562  const Real theta_avg = t_avg(k) / rho_avg;
563  cell_src(i, j, k, Rho_comp) -= drag_coefficient * (cell_data(i,j,k,Rho_comp) - rho_avg);
564  cell_src(i, j, k, RhoTheta_comp) -= drag_coefficient * (cell_data(i,j,k,RhoTheta_comp) - rho_avg * theta_avg);
565  }
566  });
567 }
RhoKE_comp
#define RhoKE_comp
Definition: ERF_IndexDefines.H:55
MicVar_Kess::rho
@ rho
Definition: ERF_Kessler.H:25
ib_wall::wall_tke
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real wall_tke(const WallState &w, amrex::Real beta, amrex::Real dz, amrex::Real z0, amrex::Real theta, amrex::Real Cmu0) noexcept
Definition: ERF_ImmersedWallCell.H:251
ib_wall::heat_source
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real heat_source(const WallState &w, amrex::Real theta, const WallCond &c, amrex::Real dz, amrex::Real dt) noexcept
Definition: ERF_ImmersedWallCell.H:215
SolverChoice::if_wall_tke_time_factor
amrex::Real if_wall_tke_time_factor
1 (the default) pins the wall cell after every stage and in the implicit TKE solve; a larger value re...
Definition: ERF_DataStruct.H:2717
ib_wall::WallState::thetastar
amrex::Real thetastar
Definition: ERF_ImmersedWallCell.H:136

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

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

Referenced by make_mom_sources().

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

◆ ImmersedForcingTerrain_Xmom_FractionStress()

void ImmersedForcingTerrain_Xmom_FractionStress ( 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< Real > &  xmom_src_arr,
const Geometry &  geom,
const SolverChoice &  solverChoice,
const Real  fac,
const Real  time 
)
466 {
467  ImmersedForcingTerrain_HorizMom_FractionStress<0>(tbx, u, v, w, cell_data, t_blank_arr, t_blank_xface_arr,
468  xmom_src_arr, geom, solverChoice, fac, time);
469 }

Referenced by make_mom_sources().

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

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

Referenced by make_mom_sources().

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

◆ ImmersedForcingTerrain_Ymom_FractionStress()

void ImmersedForcingTerrain_Ymom_FractionStress ( 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< Real > &  ymom_src_arr,
const Geometry &  geom,
const SolverChoice &  solverChoice,
const Real  fac,
const Real  time 
)
483 {
484  ImmersedForcingTerrain_HorizMom_FractionStress<1>(tby, u, v, w, cell_data, t_blank_arr, t_blank_yface_arr,
485  ymom_src_arr, geom, solverChoice, fac, time);
486 }

Referenced by make_mom_sources().

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

585 {
586  // geometric properties
587  const auto dx_arr = geom.CellSizeArray();
588  const Real dx_x = dx_arr[0];
589  const Real dx_y = dx_arr[1];
590  const Real dt = fac; // fac is actually dt in the calling code
591 
592  const Real alpha_m = solverChoice.if_Cd_momentum;
593  const Real tiny = std::numeric_limits<amrex::Real>::epsilon();
594  const bool l_implicit_drag = solverChoice.if_implicit_drag;
595  // erf.if_implicit_projection: the drag of the solid is applied inside the anelastic projection
596  const bool l_drag = !solverChoice.if_implicit_projection;
597 
598  const Real small_volfrac = solverChoice.if_small_volfrac;
599  // The terrain kernels keep the raw fractions: erf.if_snap_partial_cells
600  // applies to the buildings kernels only (the wall law below is weighted
601  // by the fluid fraction of the face, which a snapped face has not).
602 
603  ParallelFor(tbz, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
604  {
605  const Real ux = fourth * ( u(i , j , k ) + u(i+1, j , k )
606  + u(i , j , k-1) + u(i+1, j , k-1) );
607  const Real uy = fourth * ( v(i , j , k ) + v(i , j+1, k )
608  + v(i , j , k-1) + v(i , j+1, k-1) );
609  const Real uz = w(i, j, k);
610  const Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz);
611  // Use face-centered terrain_blanking if available, otherwise average from cell centers
612  Real t_blank_raw = (t_blank_zface_arr) ? t_blank_zface_arr(i, j, k) :
613  myhalf * (t_blank_arr(i, j, k) + t_blank_arr(i, j, k-1));
614  const Real t_blank = (t_blank_raw < small_volfrac) ? zero : t_blank_raw;
615 
616  const Real dx_z = (z_cc_arr) ? (z_cc_arr(i,j,k) - z_cc_arr(i,j,k-1)) : dx_arr[2];
617  const Real drag_coefficient = alpha_m / std::pow(dx_x*dx_y*dx_z, one/three);
618  const Real CdM = std::min(drag_coefficient / (windspeed + tiny), drag_coefficient);
619 
620  const Real rho_zface = myhalf * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j,k-1,Rho_comp) );
621  const Real lambda = t_blank * CdM * windspeed; // linear drag rate [1/s]
622  const Real fac_local = l_implicit_drag ? lambda / (one + lambda*dt) : lambda; // point-implicit rescale (else explicit)
623  if (l_drag) { zmom_src_arr(i, j, k) -= fac_local * rho_zface * uz; }
624  });
625 }

Referenced by make_mom_sources().

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

◆ make_wall_cond()

static ib_wall::WallCond make_wall_cond ( const SolverChoice &  sc,
const Real  time 
)
static

Surface condition of the fraction-stress wall law from the erf.if_* inputs at this time: a given Obukhov length, a given heat flux, a given surface temperature (erf.if_init_surf_temp, changing at erf.if_surf_heating_rate, which is stored in K/s), or none (neutral). The inputs guarantee that at most one is set.

329 {
330  ib_wall::WallCond c;
331  if (sc.if_Olen_in != Real(1e-8)) {
332  c.type = ib_wall::WallCond::obukhov;
333  c.Linv = one / sc.if_Olen_in;
334  } else if (sc.if_surf_temp_flux != Real(1e-8)) {
335  c.type = ib_wall::WallCond::heat_flux;
336  c.q = sc.if_surf_temp_flux;
337  } else if (sc.if_init_surf_temp > zero) {
338  c.type = ib_wall::WallCond::surface_temp;
339  c.theta_s = sc.if_init_surf_temp + sc.if_surf_heating_rate * time;
340  }
341  return c;
342 }
ib_wall::WallCond::heat_flux
static constexpr int heat_flux
given kinematic surface heat flux q [K m/s]
Definition: ERF_ImmersedWallCell.H:123
ib_wall::WallCond::Linv
amrex::Real Linv
Definition: ERF_ImmersedWallCell.H:129
ib_wall::WallCond::obukhov
static constexpr int obukhov
given Obukhov length (inverse Linv) [1/m]
Definition: ERF_ImmersedWallCell.H:125
ib_wall::WallCond::type
int type
Definition: ERF_ImmersedWallCell.H:126
ib_wall::WallCond::surface_temp
static constexpr int surface_temp
given surface potential temperature theta_s [K]
Definition: ERF_ImmersedWallCell.H:124
ib_wall::WallCond::q
amrex::Real q
Definition: ERF_ImmersedWallCell.H:127
ib_wall::WallCond::theta_s
amrex::Real theta_s
Definition: ERF_ImmersedWallCell.H:128

Referenced by ImmersedForcingTerrain_HorizMom_FractionStress(), and ImmersedForcingTerrain_Scalar_FractionStress().

Here is the caller graph for this function: