ibseb Namespace Reference

Classes
struct	RayMap

Enumerations
enum	RayHit : int { RAY_SKY = 0 , RAY_GROUND = 1 , RAY_BUILDING = 2 }
	What a ray from a face ends on. More...

Functions
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real	solve_skin_balance (amrex::Real T_start, amrex::Real SW_abs, amrex::Real Q_abs, amrex::Real LW_ext, amrex::Real emis, amrex::Real f_bldg, amrex::Real C_H, amrex::Real exner, amrex::Real theta_air, amrex::Real LE, amrex::Real slab_a, amrex::Real slab_b, amrex::Real T_min, amrex::Real T_max, amrex::Real max_step, amrex::Real tol, int max_iter, int &n_iter, amrex::Real &resid)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real	advance_slab_dirichlet (amrex::Real *T, amrex::Real T_skin, amrex::Real T_int, amrex::Real k, amrex::Real rho_cp, amrex::Real dz, amrex::Real dt, int N)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void	slab_skin_response (const amrex::Real *T, amrex::Real T_int, amrex::Real k, amrex::Real rho_cp, amrex::Real dz, amrex::Real dt, int N, amrex::Real &a, amrex::Real &b)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real	spencer_day_angle (amrex::Real day_of_year)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real	day_of_year_at (amrex::Real day0, amrex::Real t_utc_s)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real	solar_declination (amrex::Real day_of_year)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real	equation_of_time (amrex::Real day_of_year)
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_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real	orbital_hour_angle (double calday, amrex::Real longitude_rad)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real	solar_zenith (amrex::Real latitude_deg, amrex::Real decl, amrex::Real hour_angle)
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)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real	earth_sun_distance_factor (amrex::Real day_of_year)
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)
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)
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)
int	ray_halo_cells (amrex::Real height, amrex::Real dx, amrex::Real sz_min)
amrex::Real	hemisphere_min_sz (int n_az, int n_el)
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 RayMap &m)
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 RayMap &m)
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)

Variables
constexpr amrex::Real	SIGMA_SB = 5.670374419e-8
	Stefan-Boltzmann constant [W/m2/K4]. More...
constexpr int	SLAB_MAX_LAYERS = 32

Enumeration Type Documentation

RayHit

enum ibseb::RayHit : int

What a ray from a face ends on.

Enumerator
RAY_SKY
RAY_GROUND
RAY_BUILDING

: int { RAY_SKY = 0, RAY_GROUND = 1, RAY_BUILDING = 2 };

Function Documentation

advance_slab_dirichlet()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real ibseb::advance_slab_dirichlet	(	amrex::Real *	T,
		amrex::Real	T_skin,
		amrex::Real	T_int,
		amrex::Real	k,
		amrex::Real	rho_cp,
		amrex::Real	dz,
		amrex::Real	dt,
		int	N
	)

Advance one slab column by one step and return the conduction into it.

Solves, for layers l = 0 .. N-1 with implicit Euler, rho_cp dz (T_l^{n+1} - T_l^n) / dt = flux_in - flux_out, the fluxes between neighbouring centres being k (T_{l-1} - T_l) / dz, the flux from the skin 2 k (T_skin - T_0) / dz and the flux to the interior 2 k (T_{N-1} - T_int) / dz. Unconditionally stable.

Parameters

[in,out]	T	Layer temperatures [K], N values, top first.
[in]	T_skin	Skin temperature at the top [K] (Dirichlet).
[in]	T_int	Interior temperature at the bottom [K] (Dirichlet).
[in]	k	Thermal conductivity [W/m/K].
[in]	rho_cp	Volumetric heat capacity [J/m3/K].
[in]	dz	Layer thickness [m].
[in]	dt	Time step [s].
[in]	N	Number of layers, at most SLAB_MAX_LAYERS.

Returns

Conduction into the slab through the skin at the new time, 2 k (T_skin - T_0^{n+1}) / dz [W/m2], positive into the slab.

{
    if (N < 1 || N > SLAB_MAX_LAYERS) return 0.0;
    const amrex::Real Fourier = k * dt / (rho_cp * dz * dz);   // interior Fourier number
    // Coefficients of a l T_{l-1} + b_l T_l + c_l T_{l+1} = d_l (all-plus form).
    amrex::Real alpha[SLAB_MAX_LAYERS], gamma[SLAB_MAX_LAYERS];
    // Row 0: skin above at half spacing (coefficient 2 Fourier), layer 1 below.
    {
        const amrex::Real b0 = 1.0 + 2.0 * Fourier + ((N > 1) ? Fourier : 2.0 * Fourier);
        const amrex::Real c0 = (N > 1) ? -Fourier : 0.0;
        const amrex::Real d0 = T[0] + 2.0 * Fourier * T_skin + ((N > 1) ? 0.0 : 2.0 * Fourier * T_int);
        alpha[0] = c0 / b0;
        gamma[0] = d0 / b0;
    }
    // Interior rows.
    for (int l = 1; l < N - 1; ++l) {
        const amrex::Real a = -Fourier, b = 1.0 + 2.0 * Fourier, c = -Fourier;
        const amrex::Real bm = b - a * alpha[l - 1];
        gamma[l] = (T[l] - a * gamma[l - 1]) / bm;
        alpha[l] = c / bm;
    }
    // Last row (N > 1): interior below at half spacing (coefficient 2 Fourier).
    if (N > 1) {
        const int l = N - 1;
        const amrex::Real a = -Fourier, b = 1.0 + Fourier + 2.0 * Fourier;
        const amrex::Real bm = b - a * alpha[l - 1];
        gamma[l] = (T[l] + 2.0 * Fourier * T_int - a * gamma[l - 1]) / bm;
        alpha[l] = 0.0;
    }
    // Back substitution.
    T[N - 1] = gamma[N - 1];
    for (int l = N - 2; l >= 0; --l) { T[l] = gamma[l] - alpha[l] * T[l + 1]; }
    return 2.0 * k * (T_skin - T[0]) / dz;
}

Referenced by IBFaceSet::compute_ground(), slab_skin_response(), and IBFaceSet::solve_balance().

Here is the caller graph for this function:

clear_sky_diffuse_h()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real ibseb::clear_sky_diffuse_h	(	amrex::Real	cos_zenith,
		amrex::Real	S0,
		amrex::Real	tau,
		amrex::Real	distance_factor,
		amrex::Real	k_d
	)

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

{
    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);
}

Referenced by IBFaceSet::compute_shortwave().

Here is the caller graph for this function:

clear_sky_dni()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real ibseb::clear_sky_dni	(	amrex::Real	cos_zenith,
		amrex::Real	S0,
		amrex::Real	tau,
		amrex::Real	distance_factor
	)

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.

{
    if (cos_zenith <= 1.0e-3) return 0.0;
    return S0 * distance_factor * std::pow(tau, 1.0 / cos_zenith);
}

Referenced by IBFaceSet::compute_shortwave().

Here is the caller graph for this function:

day_of_year_at()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real ibseb::day_of_year_at	(	amrex::Real	day0,
		amrex::Real	t_utc_s
	)

Day of year at a UTC time: the start day plus the whole days elapsed since 00:00 UTC of that day (t_utc_s may exceed 86400 s in a multi-day run). The Spencer series repeat every 365 days, so the result needs no wrapping.

{
    return day0 + std::floor(t_utc_s / 86400.0);
}

Referenced by IBFaceSet::compute_shortwave().

Here is the caller graph for this function:

earth_sun_distance_factor()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real ibseb::earth_sun_distance_factor

(

amrex::Real

day_of_year

)

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.

{
    const amrex::Real G = spencer_day_angle(day_of_year);
    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);
}

Referenced by IBFaceSet::compute_shortwave().

Here is the call graph for this function:

Here is the caller graph for this function:

equation_of_time()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real ibseb::equation_of_time

(

amrex::Real

day_of_year

)

Equation of time [minutes] from the day of year (Spencer 1971).

{
    const amrex::Real gamma = spencer_day_angle(day_of_year);
    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));
}

Referenced by solar_hour_angle().

Here is the call graph for this function:

Here is the caller graph for this function:

hemisphere_direction()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void ibseb::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
	)

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.

{
    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];
}

Referenced by IBFaceSet::compute_view_fractions().

Here is the call graph for this function:

Here is the caller graph for this function:

hemisphere_min_sz()

amrex::Real ibseb::hemisphere_min_sz ( int  n_az, int  n_el  )

inline

Sine of the elevation of the shallowest ray of the hemisphere sampling (hemisphere_direction()) with n_az azimuths and n_el elevations: on a roof the ray furthest from the normal, cos(theta_max); on a wall, where a ray's sine of elevation is sin(theta) sin(phi), the ray nearest the normal at the azimuth nearest the horizontal, sin(theta_min) sin(pi / n_az). With an odd n_az a wall has level rays (phi = pi), which the walk follows for two rounds of the map.

{
    const amrex::Real u_min = 0.5 / n_el, u_max = (n_el - 0.5) / n_el;
    const amrex::Real roof  = std::sqrt(1.0 - u_max);
    const amrex::Real wall  = std::sqrt(u_min) * std::sin(PI / n_az);
    return amrex::min(roof, wall);
}

Referenced by IBFaceSet::build_ray_map().

Here is the caller graph for this function:

orbital_hour_angle()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real ibseb::orbital_hour_angle	(	double	calday,
		amrex::Real	longitude_rad
	)

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.

{
    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));
}

Referenced by ERF::ibseb_set_two_stream_sun().

Here is the caller graph for this function:

ray_blocked()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool ibseb::ray_blocked	(	amrex::Real	x0,
		amrex::Real	y0,
		amrex::Real	z0,
		amrex::Real	sx,
		amrex::Real	sy,
		amrex::Real	sz,
		const RayMap &	m
	)

Whether the ray from a face toward the sun hits a building.

Buildings are the column tops of the rank's ray map (RayMap), 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 or leaves the map (the far field is open ground; in a periodic direction the map spans, the ray wraps round instead).

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.

Returns

true when blocked (the face is in shadow).

{
    if (sz <= 0.0) return true;              // sun at or below the horizon
    if (z0 >= m.z_max) return false;         // nothing stands this high
    const amrex::Real dx = m.dx, dy = m.dy;
    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 = m.ilo + static_cast<int>(std::floor((x - m.x_lo) / dx));
    int cj = m.jlo + static_cast<int>(std::floor((y - m.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 = ((m.x_lo + (ci - m.ilo + 1) * dx) - x) / sx;
    else if (stepx < 0) tmx = ((m.x_lo + (ci - m.ilo) * dx) - x) / sx;
    if (stepy > 0)      tmy = ((m.y_lo + (cj - m.jlo + 1) * dy) - y) / sy;
    else if (stepy < 0) tmy = ((m.y_lo + (cj - m.jlo) * dy) - y) / sy;
    // Columns crossed before the ray passes the tallest one, plus a margin,
    // and never more than the halo's reach (see RayMap).
    const amrex::Real t_max = (m.z_max - z) / sz;
    const amrex::Real steps = amrex::min(t_max * (std::abs(sx) / dx + std::abs(sy) / dy), amrex::Real(m.max_steps));
    const int max_iter = static_cast<int>(amrex::min(amrex::Real(1.0e9), steps)) + 8;
    for (int it = 0; it < max_iter; ++it) {
        // Column (ci, cj) with the ray entering at height z.
        if (!m.locate(ci, cj)) return false;
        if (z + eps < m.column_top(ci, cj)) return true;
        // Step to the next column.
        amrex::Real t;
        if (tmx < tmy) { t = tmx; tmx += tdx; ci += stepx; }
        else           { t = tmy; tmy += tdy; cj += stepy; }
        z = z0 + t0 * sz + t * sz;
        if (z >= m.z_max) return false;
    }
    return false;
}

Referenced by IBFaceSet::compute_shortwave(), and IBFaceSet::compute_shortwave_two_stream().

Here is the call graph for this function:

Here is the caller graph for this function:

ray_halo_cells()

int ibseb::ray_halo_cells ( amrex::Real  height, amrex::Real  dx, amrex::Real  sz_min  )

inline

Halo of a ray map, in cells: how far a ray no shallower than sz_min (the sine of its elevation) can travel horizontally before it rises above height or falls that far, height / tan(elevation), rounded up and with one cell for the column it starts in. At least one cell, so a face always finds the column of its solid neighbour; a level ray would reach for ever, which the caller caps at the domain.

{
    if (height <= 0.0) { return 1; }
    const amrex::Real cz = std::sqrt(amrex::max(amrex::Real(0.0), amrex::Real(1.0) - sz_min * sz_min));
    if (sz_min <= 0.0 || cz <= 0.0) { return std::numeric_limits<int>::max() / 4; }
    const amrex::Real reach = height * cz / sz_min / dx;
    if (reach >= amrex::Real(std::numeric_limits<int>::max() / 8)) { return std::numeric_limits<int>::max() / 4; }
    return amrex::max(1, static_cast<int>(std::ceil(reach)) + 1);
}

Referenced by IBFaceSet::build_ray_map().

Here is the caller graph for this function:

ray_hit()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE int ibseb::ray_hit	(	amrex::Real	x0,
		amrex::Real	y0,
		amrex::Real	z0,
		amrex::Real	sx,
		amrex::Real	sy,
		amrex::Real	sz,
		const RayMap &	m
	)

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 leaves the map, or has crossed the domain twice where the map wraps, 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. Leaving the map counts as sky for a rising ray and as ground for a falling one: the far field is open ground. Where the map spans a periodic direction the ray wraps round instead.

{
    const amrex::Real dx = m.dx, dy = m.dy;
    const amrex::Real eps = 1.0e-6 * amrex::min(dx, dy);
    if (sz > 0.0 && z0 >= m.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 = m.ilo + static_cast<int>(std::floor((x - m.x_lo) / dx));
    int cj = m.jlo + static_cast<int>(std::floor((y - m.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 = ((m.x_lo + (ci - m.ilo + 1) * dx) - x) / sx;
    else if (stepx < 0) tmx = ((m.x_lo + (ci - m.ilo) * dx) - x) / sx;
    if (stepy > 0)      tmy = ((m.y_lo + (cj - m.jlo + 1) * dy) - y) / sy;
    else if (stepy < 0) tmy = ((m.y_lo + (cj - m.jlo) * 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;
    }
    // Columns crossed before the ray passes the tallest column or reaches
    // the ground, plus a margin, and never more than the halo's reach (see
    // RayMap). A ray within 1e-9 of level (phi = pi of an odd n_az gives sz
    // of 1e-16) is a level ray and gets two rounds of the map.
    const bool level = (std::abs(sz) < 1.0e-9);
    const amrex::Real t_max = (sz > 0.0) ? (m.z_max - z) / sz : ((sz < 0.0) ? (z - m.z_ground) / (-sz) : big);
    const amrex::Real rounds = amrex::Real(2 * ((m.top.end[0] - m.top.begin[0]) + (m.top.end[1] - m.top.begin[1])));
    const amrex::Real steps = level ? rounds
                            : amrex::min(t_max * (std::abs(sx) / dx + std::abs(sy) / dy), amrex::Real(m.max_steps));
    const int max_iter = static_cast<int>(amrex::min(amrex::Real(1.0e9), steps)) + 8;
    for (int it = 0; it < max_iter; ++it) {
        if (!m.locate(ci, cj)) return (sz >= 0.0) ? RAY_SKY : RAY_GROUND;
        const amrex::Real top = m.column_top(ci, cj);
        const bool solid = (top > m.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 <= m.z_ground + eps) return RAY_GROUND;
        }
        if (tmx < tmy) { tmx += tdx; ci += stepx; }
        else           { tmy += tdy; cj += stepy; }
        z = z_out;
        if (sz > 0.0 && z >= m.z_max) return RAY_SKY;
    }
    return (sz >= 0.0) ? RAY_SKY : RAY_GROUND;
}

Referenced by IBFaceSet::compute_view_fractions().

Here is the call graph for this function:

Here is the caller graph for this function:

slab_skin_response()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void ibseb::slab_skin_response	(	const amrex::Real *	T,
		amrex::Real	T_int,
		amrex::Real	k,
		amrex::Real	rho_cp,
		amrex::Real	dz,
		amrex::Real	dt,
		int	N,
		amrex::Real &	a,
		amrex::Real &	b
	)

Linear response of the implicit slab step to the skin temperature.

The step of advance_slab_dirichlet() is linear in T_skin, so the conduction it returns is G = a T_skin - b for the slab's current state, interior temperature and time step. The two coefficients come from two trial steps on copies of the column (T_skin = 0 and T_skin = 1), the same arithmetic as the real step, so a balance solved with this response and then advanced with the resulting skin temperature closes to rounding. The column itself is not modified.

Parameters

[in]	T	Layer temperatures [K] at the start of the step.
[out]	a	Slope of G with respect to T_skin [W/m2/K], positive.
[out]	b	Offset: G at T_skin = 0 is -b [W/m2]. Other arguments as in advance_slab_dirichlet().

{
    if (N < 1 || N > SLAB_MAX_LAYERS) { a = 0.0; b = 0.0; return; }
    // Zero-initialized: only the first N entries are ever used, but the
    // compiler cannot see that through the inlined solve and warns otherwise.
    amrex::Real W[SLAB_MAX_LAYERS] = {};
    for (int l = 0; l < N; ++l) { W[l] = T[l]; }
    const amrex::Real G0 = advance_slab_dirichlet(W, 0.0, T_int, k, rho_cp, dz, dt, N);
    for (int l = 0; l < N; ++l) { W[l] = T[l]; }
    const amrex::Real G1 = advance_slab_dirichlet(W, 1.0, T_int, k, rho_cp, dz, dt, N);
    a = G1 - G0;
    b = -G0;
}

Referenced by IBFaceSet::solve_balance().

Here is the call graph for this function:

Here is the caller graph for this function:

solar_azimuth()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real ibseb::solar_azimuth	(	amrex::Real	latitude_deg,
		amrex::Real	decl,
		amrex::Real	hour_angle,
		amrex::Real	zenith
	)

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.

{
    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;
}

Referenced by IBFaceSet::compute_shortwave().

Here is the caller graph for this function:

solar_declination()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real ibseb::solar_declination

(

amrex::Real

day_of_year

)

Solar declination [rad] from the day of year (Spencer 1971).

{
    const amrex::Real gamma = spencer_day_angle(day_of_year);
    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.002697 * std::cos(3.0 * gamma) + 0.00148 * std::sin(3.0 * gamma);
}

Referenced by IBFaceSet::compute_shortwave().

Here is the call graph for this function:

Here is the caller graph for this function:

solar_hour_angle()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real ibseb::solar_hour_angle	(	amrex::Real	time_utc_s,
		amrex::Real	longitude_deg,
		amrex::Real	day_of_year
	)

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::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;
}

Referenced by IBFaceSet::compute_shortwave().

Here is the call graph for this function:

Here is the caller graph for this function:

solar_zenith()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real ibseb::solar_zenith	(	amrex::Real	latitude_deg,
		amrex::Real	decl,
		amrex::Real	hour_angle
	)

Zenith angle [rad] from the latitude, declination and 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);
}

Referenced by IBFaceSet::compute_shortwave().

Here is the caller graph for this function:

solve_skin_balance()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real ibseb::solve_skin_balance	(	amrex::Real	T_start,
		amrex::Real	SW_abs,
		amrex::Real	Q_abs,
		amrex::Real	LW_ext,
		amrex::Real	emis,
		amrex::Real	f_bldg,
		amrex::Real	C_H,
		amrex::Real	exner,
		amrex::Real	theta_air,
		amrex::Real	LE,
		amrex::Real	slab_a,
		amrex::Real	slab_b,
		amrex::Real	T_min,
		amrex::Real	T_max,
		amrex::Real	max_step,
		amrex::Real	tol,
		int	max_iter,
		int &	n_iter,
		amrex::Real &	resid
	)

Solve the face balance for the skin temperature.

Parameters

[in]	T_start	Starting guess, the skin temperature of the previous step [K].
[in]	SW_abs	Absorbed shortwave [W/m2].
[in]	Q_abs	Absorbed external flux eps Q_ext [W/m2].
[in]	LW_ext	Incoming longwave from sky and ground [W/m2].
[in]	emis	Face emissivity.
[in]	f_bldg	Building view fraction of the face.
[in]	C_H	Sensible coefficient rho c_p kappa u* / ln(delta/z0h) [W/m2/K].
[in]	exner	Exner function of the fluid cell, T_a / theta_a.
[in]	theta_air	Potential temperature of the fluid cell [K].
[in]	LE	Latent flux, frozen [W/m2].
[in]	slab_a,slab_b	Slab response G = a T_s - b [W/m2/K], [W/m2].
[in]	T_min,T_max	Bounds on the skin temperature [K].
[in]	max_step	Largest Newton step [K].
[in]	tol	Convergence: stop when a step changes T_s by less than this [K].
[in]	max_iter	Iteration cap.
[out]	n_iter	Iterations taken.
[out]	resid	|F| at the returned temperature [W/m2]; non-zero when a bound was hit or the cap was reached.

Returns

The skin temperature [K], within the bounds.

{
    using amrex::Real;
    const Real e_eff = emis * (Real(1.0) - f_bldg) * SIGMA_SB;   // effective emission coefficient
    auto residual = [&](Real T) {
        const Real T4 = T * T * T * T;
        return SW_abs + Q_abs + emis * LW_ext - e_eff * T4
             - C_H * (T / exner - theta_air) - LE - (slab_a * T - slab_b);
    };
    Real T = amrex::max(amrex::min(T_start, T_max), T_min);
    n_iter = 0;
    Real F = residual(T);
    for (int it = 0; it < max_iter; ++it) {
        ++n_iter;
        const Real Fp = -Real(4.0) * e_eff * T * T * T - C_H / exner - slab_a;
        Real step = -F / Fp;                                     // Fp < 0 always
        step = amrex::max(amrex::min(step, max_step), -max_step);
        const Real T_new = amrex::max(amrex::min(T + step, T_max), T_min);
        const Real dT = std::abs(T_new - T);
        T = T_new;
        F = residual(T);
        if (dT < tol) break;
    }
    resid = std::abs(F);
    return T;
}

Referenced by IBFaceSet::solve_balance().

Here is the caller graph for this function:

spencer_day_angle()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real ibseb::spencer_day_angle

(

amrex::Real

day_of_year

)

Day angle [rad] of Spencer's (1971) series, 2 pi (day - 1) / 365. The series' coefficients were fitted with this day angle; the declination, the equation of time and the Earth-Sun distance below all use it, so day 366 gives the same values as day 1.

{
    return 2.0 * PI * (day_of_year - 1.0) / 365.0;
}

Referenced by earth_sun_distance_factor(), equation_of_time(), and solar_declination().

Here is the caller graph for this function:

sun_vector()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void ibseb::sun_vector	(	amrex::Real	zenith,
		amrex::Real	azimuth,
		amrex::Real &	sx,
		amrex::Real &	sy,
		amrex::Real &	sz
	)

Unit vector toward the sun from zenith and azimuth (east +x, north +y).

{
    const amrex::Real s = std::sin(zenith);
    sx = s * std::sin(azimuth);
    sy = s * std::cos(azimuth);
    sz = std::cos(zenith);
}

Referenced by IBFaceSet::compute_shortwave().

Here is the caller graph for this function:

Variable Documentation

SIGMA_SB

constexpr amrex::Real ibseb::SIGMA_SB = 5.670374419e-8

constexpr

Stefan-Boltzmann constant [W/m2/K4].

Referenced by IBFaceSet::compute_longwave(), IBFaceSet::compute_longwave_two_stream(), IBFaceSet::solve_balance(), and solve_skin_balance().

SLAB_MAX_LAYERS

constexpr int ibseb::SLAB_MAX_LAYERS = 32

constexpr

Largest number of layers the stack arrays allow. The kernels keep three such arrays per thread, so the bound is kept small for GPU occupancy; 32 layers resolve any wall (erf.ibseb.n_slab_layers is checked against it).

Referenced by advance_slab_dirichlet(), IBFaceSet::compute_ground(), slab_skin_response(), and IBFaceSet::solve_balance().