ERF_IBSEBSolar.H

Go to the documentation of this file.

#ifndef ERF_IBSEB_SOLAR_H
#define ERF_IBSEB_SOLAR_H
 
#include <AMReX_REAL.H>
#include <AMReX_GpuQualifiers.H>
#include <AMReX_Algorithm.H>
#include <ERF_NumericalConstants.H>
#include <cmath>
 
/**
 * \file ERF_IBSEBSolar.H
 * \brief Sun position, clear-sky shortwave, and the ray cast that decides
 *        whether a building face is in shadow.
 *
 * The solar geometry (Spencer's declination and equation of time, the hour
 * angle from UTC and the site longitude, the zenith, the clockwise-from-north
 * azimuth, the Earth-Sun distance) is kept here: the radiation models use
 * the CESM orbital formulation of ERF_OrbCosZenith.H, which gives no
 * azimuth and has no equation of time; orbital_hour_angle() turns its
 * calendar day into the hour angle these functions take, for
 * erf.ibseb.sun_mode = two_stream. The clear-sky formulas follow Bird (1984) as used by the SLUCM
 * branch. The ray cast is new here.
 *
 * **Sun vector convention.** ``ibseb_sun_vector()`` returns the unit vector
 * from a surface toward the sun, with azimuth measured clockwise from north:
 * east is +x, north is +y, so ``s = (sin z sin a, sin z cos a, cos z)`` with
 * z the zenith angle and a the azimuth. A face with outward normal n receives
 * direct beam only when ``n . s > 0``.
 *
 * References: Spencer (1971); Duffie and Beckman (1991); Iqbal (1983);
 * Bird and Hulstrom (1981).
 */
namespace ibseb {
 
/** Solar declination [rad] from the day of year (Spencer 1971). */
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
amrex::Real solar_declination (amrex::Real day_of_year)
{
    const amrex::Real gamma = 2.0 * PI * (day_of_year - 1.0) / 365.25;
    return 0.006918
        - 0.399912 * std::cos(gamma) + 0.070257 * std::sin(gamma)
        - 0.006758 * std::cos(2.0 * gamma) + 0.000907 * std::sin(2.0 * gamma)
        - 0.00248 * std::cos(3.0 * gamma) + 0.00031 * std::sin(3.0 * gamma);
}
 
/** Equation of time [minutes] from the day of year (Spencer 1971). */
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
amrex::Real equation_of_time (amrex::Real day_of_year)
{
    const amrex::Real gamma = 2.0 * PI * (day_of_year - 1.0) / 365.25;
    return 229.18 * (0.000075 + 0.001868 * std::cos(gamma) - 0.032077 * std::sin(gamma)
                     - 0.014615 * std::cos(2.0 * gamma) - 0.040849 * std::sin(2.0 * gamma));
}
 
/**
 * Hour angle [rad] from the UTC time of day and the site: solar time is UTC
 * plus one hour per 15 degrees of longitude (east positive) plus the equation
 * of time, and the hour angle is 15 degrees per hour from solar noon,
 * negative in the morning. The time zone plays no part (it would cancel).
 */
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
amrex::Real solar_hour_angle (amrex::Real time_utc_s, amrex::Real longitude_deg, amrex::Real day_of_year)
{
    amrex::Real t_utc_h = std::fmod(time_utc_s / 3600.0, 24.0);
    if (t_utc_h < 0.0) { t_utc_h += 24.0; }
    const amrex::Real t_solar_h = t_utc_h + longitude_deg / 15.0 + equation_of_time(day_of_year) / 60.0;
    return 15.0 * (t_solar_h - 12.0) * PI / 180.0;
}
 
/**
 * Hour angle [rad] of the orbital formula the radiation models use
 * (orbital_cos_zenith_instant() of ERF_OrbCosZenith.H): its
 * ``cos z = sin(lat) sin(decl) - cos(lat) cos(decl) cos(2 pi frac(calday) + lon)``
 * is the usual form with the hour angle ``2 pi frac(calday) + lon - pi``. The
 * fraction of the calendar day is the UTC time of day and there is no equation
 * of time, so local solar time is UTC plus one hour per 15 degrees of
 * longitude. Wrapped to [-pi, pi), negative in the morning.
 */
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
amrex::Real orbital_hour_angle (double calday, amrex::Real longitude_rad)
{
    const double two_pi = 2.0 * double(PI);
    double h = two_pi * (calday - std::floor(calday)) + double(longitude_rad) - double(PI);
    h = std::fmod(h + double(PI), two_pi);
    if (h < 0.0) { h += two_pi; }
    return static_cast<amrex::Real>(h - double(PI));
}
 
/** Zenith angle [rad] from the latitude, declination and hour angle. */
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
amrex::Real solar_zenith (amrex::Real latitude_deg, amrex::Real decl, amrex::Real hour_angle)
{
    const amrex::Real lat = latitude_deg * PI / 180.0;
    amrex::Real cz = std::sin(lat) * std::sin(decl) + std::cos(lat) * std::cos(decl) * std::cos(hour_angle);
    if (cz < -1.0) { cz = -1.0; }
    if (cz >  1.0) { cz =  1.0; }
    return std::acos(cz);
}
 
/**
 * Azimuth [rad] clockwise from north (east is pi/2) from the latitude,
 * declination, hour angle and zenith; north when the sun is at the zenith or
 * the site is at a pole, where the azimuth is undefined.
 */
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
amrex::Real solar_azimuth (amrex::Real latitude_deg, amrex::Real decl,
                           amrex::Real hour_angle, amrex::Real zenith)
{
    const amrex::Real lat = latitude_deg * PI / 180.0;
    const amrex::Real sz = std::sin(zenith);
    const amrex::Real cz = std::cos(zenith);
    if (sz < 1.0e-6) { return 0.0; }
    const amrex::Real cos_lat = std::cos(lat);
    if (std::abs(cos_lat) < 1.0e-6) { return 0.0; }
    // The hour angle is negative before solar noon, when the sun is east of
    // south, hence the sign for the clockwise-from-north convention.
    const amrex::Real sin_az = -std::sin(hour_angle) * std::cos(decl) / sz;
    const amrex::Real cos_az = (std::sin(decl) - cz * std::sin(lat)) / (sz * cos_lat);
    amrex::Real az = std::atan2(sin_az, cos_az);
    if (az < 0.0) { az += 2.0 * PI; }
    return az;
}
 
/**
 * Earth-Sun distance factor (d0/d)^2 that scales the solar constant for the
 * day of year (Spencer 1971): about 1.034 at perihelion in early January,
 * 0.967 at aphelion in early July, annual mean 1.
 */
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
amrex::Real earth_sun_distance_factor (amrex::Real day_of_year)
{
    const amrex::Real G = 2.0 * PI * (day_of_year - 1.0) / 365.0;
    return 1.000110 + 0.034221 * std::cos(G) + 0.001280 * std::sin(G)
         + 0.000719 * std::cos(2.0 * G) + 0.000077 * std::sin(2.0 * G);
}
 
/** Unit vector toward the sun from zenith and azimuth (east +x, north +y). */
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
void sun_vector (amrex::Real zenith, amrex::Real azimuth,
                 amrex::Real& sx, amrex::Real& sy, amrex::Real& sz)
{
    const amrex::Real s = std::sin(zenith);
    sx = s * std::sin(azimuth);
    sy = s * std::cos(azimuth);
    sz = std::cos(zenith);
}
 
/**
 * Clear-sky direct-normal irradiance [W/m2] (Bird form): the solar constant,
 * corrected for the Earth-Sun distance, attenuated by ``tau^(1/cos z)``.
 * Zero when the sun is below the horizon.
 */
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
amrex::Real clear_sky_dni (amrex::Real cos_zenith, amrex::Real S0, amrex::Real tau,
                           amrex::Real distance_factor)
{
    if (cos_zenith <= 1.0e-3) return 0.0;
    return S0 * distance_factor * std::pow(tau, 1.0 / cos_zenith);
}
 
/**
 * Clear-sky diffuse irradiance on a horizontal surface [W/m2]: a fraction
 * ``k_d`` of what the direct beam lost to the atmosphere reaches the ground
 * as diffuse light (Liu and Jordan form).
 */
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
amrex::Real clear_sky_diffuse_h (amrex::Real cos_zenith, amrex::Real S0, amrex::Real tau,
                                 amrex::Real distance_factor, amrex::Real k_d)
{
    if (cos_zenith <= 1.0e-3) return 0.0;
    const amrex::Real tr = std::pow(tau, 1.0 / cos_zenith);
    return k_d * S0 * distance_factor * cos_zenith * (1.0 - tr);
}
 
/**
 * Top height of a column of the level. The column map covers the bounding box
 * of the built columns only (origin ``i0, j0`` in 0-based domain columns,
 * ``bw x bh`` values indexed ``(ci - i0) * bh + (cj - j0)``); every other
 * column is open ground. An entry is the domain-relative k of the highest
 * solid cell of its column, or -1 where the column is fluid, so the top is
 * one cell above that index and a fluid column sits at the ground.
 */
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
amrex::Real column_top (const int* col_top, int ci, int cj,
                        int i0, int j0, int bw, int bh,
                        amrex::Real z_ground, amrex::Real dz)
{
    ci -= i0; cj -= j0;
    if (ci < 0 || ci >= bw || cj < 0 || cj >= bh) { return z_ground; }
    const int kt = col_top[ci * bh + cj];
    return (kt < 0) ? z_ground : z_ground + (kt + 1) * dz;
}
 
/**
 * Whether the ray from a face toward the sun hits a building.
 *
 * Buildings are represented by the top height of the columns of the level
 * (``col_top`` over the built bounding box, see column_top(), replicated on
 * every rank), so a ray is blocked wherever its height on entering a
 * column is below that column's top: buildings stand on the ground, so the
 * lowest point of the ray inside a column is at its entry, and the height
 * only grows along a ray toward a sun above the horizon. The walk over the
 * columns is a 2D digital differential analyser. It stops when the ray
 * rises above the tallest column, leaves the domain in a non-periodic
 * direction, or has travelled ``max_path``.
 *
 * The start point is nudged a little along the ray so a face does not
 * test its own column at its own height: a wall starts in its fluid cell's
 * column, whose top is below the face (no overhangs), and a roof starts a
 * hair above its own top.
 *
 * @return true when blocked (the face is in shadow).
 */
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
bool ray_blocked (amrex::Real x0, amrex::Real y0, amrex::Real z0,
                  amrex::Real sx, amrex::Real sy, amrex::Real sz,
                  const int* col_top, int nx, int ny, int i0, int j0, int bw, int bh,
                  amrex::Real x_lo, amrex::Real y_lo, amrex::Real dx, amrex::Real dy, amrex::Real dz,
                  bool per_x, bool per_y, amrex::Real z_ground, amrex::Real z_max, amrex::Real max_path)
{
    if (sz <= 0.0) return true;              // sun at or below the horizon
    if (z0 >= z_max) return false;           // nothing stands this high
    const amrex::Real eps = 1.0e-6 * amrex::min(dx, dy);
    const amrex::Real t0  = 1.0e-3 * amrex::min(dx, dy);
    amrex::Real x = x0 + t0 * sx, y = y0 + t0 * sy, z = z0 + t0 * sz;
    int ci = static_cast<int>(std::floor((x - x_lo) / dx));
    int cj = static_cast<int>(std::floor((y - y_lo) / dy));
    const int stepx = (sx > 0.0) ? 1 : ((sx < 0.0) ? -1 : 0);
    const int stepy = (sy > 0.0) ? 1 : ((sy < 0.0) ? -1 : 0);
    const amrex::Real big = 1.0e30;
    const amrex::Real tdx = (stepx != 0) ? dx / std::abs(sx) : big;
    const amrex::Real tdy = (stepy != 0) ? dy / std::abs(sy) : big;
    // Distance along the ray to the next column boundary in x and in y.
    amrex::Real tmx = big, tmy = big;
    if (stepx > 0)      tmx = ((x_lo + (ci + 1) * dx) - x) / sx;
    else if (stepx < 0) tmx = ((x_lo + ci * dx) - x) / sx;
    if (stepy > 0)      tmy = ((y_lo + (cj + 1) * dy) - y) / sy;
    else if (stepy < 0) tmy = ((y_lo + cj * dy) - y) / sy;
    amrex::Real t = 0.0;
    const int max_iter = 4 * (nx + ny) + 8;
    for (int it = 0; it < max_iter; ++it) {
        // Column (ci, cj) with the ray entering at height z.
        if (ci < 0 || ci >= nx) {
            if (!per_x) return false;
            ci = (ci % nx + nx) % nx;
        }
        if (cj < 0 || cj >= ny) {
            if (!per_y) return false;
            cj = (cj % ny + ny) % ny;
        }
        if (z + eps < column_top(col_top, ci, cj, i0, j0, bw, bh, z_ground, dz)) return true;
        // Step to the next column.
        if (tmx < tmy) { t = tmx; tmx += tdx; ci += stepx; }
        else           { t = tmy; tmy += tdy; cj += stepy; }
        z = z0 + t0 * sz + t * sz;
        if (z >= z_max) return false;
        if (t > max_path) return false;
    }
    return false;
}
 
/// What a ray from a face ends on.
enum RayHit : int { RAY_SKY = 0, RAY_GROUND = 1, RAY_BUILDING = 2 };
 
/**
 * Where a ray from a face ends: sky, ground or a building, for any
 * direction (the hemisphere sampling of the view fractions). The same column walk as
 * ray_blocked(), with the height test made direction-aware:
 *
 * - rising or level rays (``sz >= 0``) are blocked by a column whose top is
 *   above the entry height; a level ray that crosses ``max_path`` without a
 *   hit ends at the horizon, which counts as sky;
 * - falling rays (``sz < 0``) are blocked by a solid column whose top is
 *   above the exit height (the ray descends inside the column) and reach
 *   the ground where the exit height is at or below ``z_ground``.
 *
 * Solid columns are those with a top above the ground. The domain is
 * wrapped in periodic directions; leaving it in a non-periodic direction
 * counts as sky for a rising ray and as ground for a falling one.
 */
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
int ray_hit (amrex::Real x0, amrex::Real y0, amrex::Real z0,
             amrex::Real sx, amrex::Real sy, amrex::Real sz,
             const int* col_top, int nx, int ny, int i0, int j0, int bw, int bh,
             amrex::Real x_lo, amrex::Real y_lo, amrex::Real dx, amrex::Real dy, amrex::Real dz,
             bool per_x, bool per_y, amrex::Real z_ground, amrex::Real z_max,
             amrex::Real max_path)
{
    const amrex::Real eps = 1.0e-6 * amrex::min(dx, dy);
    if (sz > 0.0 && z0 >= z_max) return RAY_SKY;
    const amrex::Real t0 = 1.0e-3 * amrex::min(dx, dy);
    amrex::Real x = x0 + t0 * sx, y = y0 + t0 * sy, z = z0 + t0 * sz;
    int ci = static_cast<int>(std::floor((x - x_lo) / dx));
    int cj = static_cast<int>(std::floor((y - y_lo) / dy));
    const int stepx = (sx > 0.0) ? 1 : ((sx < 0.0) ? -1 : 0);
    const int stepy = (sy > 0.0) ? 1 : ((sy < 0.0) ? -1 : 0);
    const amrex::Real big = 1.0e30;
    const amrex::Real tdx = (stepx != 0) ? dx / std::abs(sx) : big;
    const amrex::Real tdy = (stepy != 0) ? dy / std::abs(sy) : big;
    amrex::Real tmx = big, tmy = big;
    if (stepx > 0)      tmx = ((x_lo + (ci + 1) * dx) - x) / sx;
    else if (stepx < 0) tmx = ((x_lo + ci * dx) - x) / sx;
    if (stepy > 0)      tmy = ((y_lo + (cj + 1) * dy) - y) / sy;
    else if (stepy < 0) tmy = ((y_lo + cj * dy) - y) / sy;
    if (stepx == 0 && stepy == 0) {          // vertical ray
        if (sz > 0.0) return RAY_SKY;
        // straight down from a wall face never happens (the face is vertical);
        // from a roof the hemisphere points up. Treat as ground for safety.
        return RAY_GROUND;
    }
    amrex::Real t = 0.0;
    const int max_iter = 4 * (nx + ny) + 8;
    for (int it = 0; it < max_iter; ++it) {
        if (ci < 0 || ci >= nx) {
            if (!per_x) return (sz >= 0.0) ? RAY_SKY : RAY_GROUND;
            ci = (ci % nx + nx) % nx;
        }
        if (cj < 0 || cj >= ny) {
            if (!per_y) return (sz >= 0.0) ? RAY_SKY : RAY_GROUND;
            cj = (cj % ny + ny) % ny;
        }
        const amrex::Real top = column_top(col_top, ci, cj, i0, j0, bw, bh, z_ground, dz);
        const bool solid = (top > z_ground + eps);
        const amrex::Real t_next = amrex::min(tmx, tmy);
        const amrex::Real z_out  = z0 + t0 * sz + t_next * sz;
        if (sz >= 0.0) {
            if (solid && z + eps < top) return RAY_BUILDING;
        } else {
            if (solid && z + eps < top) return RAY_BUILDING;     // entered below the top
            if (solid && z_out < top - eps) return RAY_BUILDING; // descends into it
            if (z_out <= z_ground + eps) return RAY_GROUND;
        }
        if (tmx < tmy) { t = tmx; tmx += tdx; ci += stepx; }
        else           { t = tmy; tmy += tdy; cj += stepy; }
        z = z_out;
        if (sz > 0.0 && z >= z_max) return RAY_SKY;
        if (t > max_path) return (sz >= 0.0) ? RAY_SKY : RAY_GROUND;
    }
    return (sz >= 0.0) ? RAY_SKY : RAY_GROUND;
}
 
/**
 * Direction of the ``(ia, ie)`` sample of a cosine-weighted hemisphere
 * around the outward normal of a face with direction ``dir`` and normal
 * sign ``nsign`` (+1 along the axis, -1 against it). Stratified in the
 * azimuth ``phi = 2 pi (ia + 1/2) / n_az`` and in ``u = (ie + 1/2) / n_el``
 * with ``theta = asin(sqrt(u))`` from the normal, so every sample carries
 * the same weight and the counts are view factors.
 */
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
void hemisphere_direction (int dir, int nsign, int ia, int ie, int n_az, int n_el,
                           amrex::Real& sx, amrex::Real& sy, amrex::Real& sz)
{
    const amrex::Real u     = (ie + 0.5) / n_el;
    const amrex::Real theta = std::asin(std::sqrt(u));
    const amrex::Real phi   = 2.0 * PI * (ia + 0.5) / n_az;
    const amrex::Real cn = std::cos(theta), ct = std::sin(theta) * std::cos(phi), st = std::sin(theta) * std::sin(phi);
    // Normal and two tangents for each axis.
    amrex::Real n[3] = {0.0, 0.0, 0.0}, t1[3] = {0.0, 0.0, 0.0}, t2[3] = {0.0, 0.0, 0.0};
    n[dir] = nsign;
    if (dir == 0)      { t1[1] = 1.0; t2[2] = 1.0; }
    else if (dir == 1) { t1[0] = 1.0; t2[2] = 1.0; }
    else               { t1[0] = 1.0; t2[1] = 1.0; }
    sx = cn * n[0] + ct * t1[0] + st * t2[0];
    sy = cn * n[1] + ct * t1[1] + st * t2[1];
    sz = cn * n[2] + ct * t1[2] + st * t2[2];
}
 
} // namespace ibseb
 
#endif