TopographicShading Class Reference

#include <ERF_TopographicShading.H>

Collaboration diagram for TopographicShading:

Classes
struct	Level

Public Types
enum	RadiationField {   CosZen , Shortwave , DirectVis , DirectNir ,   DiffuseVis , DiffuseNir , Longwave , NRad }
enum	Diagnostic {   Shadow , DirectFactor , Incidence , UnshadedSW ,   TerrainSW , NDiag }
enum	GeometryField {   Height , MapX , MapY , Latitude ,   Longitude , CosAlpha , SinAlpha , NGeom }

Public Member Functions
void	define_level (int lev, const amrex::BoxArray &ba2d, const amrex::DistributionMapping &dm, const amrex::Geometry &geom)
bool	geometry_ready (int lev) const
bool	forcing_ready (int lev) const
void	initialize_geometry (int lev, const TopographicShadingChoice &choice, const amrex::MultiFab &z_nd, const amrex::MultiFab &map_x, const amrex::MultiFab &map_y, const amrex::MultiFab *lat, const amrex::MultiFab *lon)
amrex::Vector< amrex::MultiFab * >	raw_fields (int lev)
amrex::MultiFab *	raw (int lev)
const amrex::MultiFab *	diagnostics (int lev) const
void	set_sun (int lev, double calday, double declin)
void	interpolate_from_coarse (int lev)
void	apply (int lev, const TopographicShadingChoice &choice, const amrex::Vector< amrex::MultiFab * > &destinations)

Private Attributes
amrex::Vector< std::unique_ptr< Level > >	m_levels

Detailed Description

Surface-only correction of the radiation supplied to land models. The atmospheric fluxes and heating remain the column solver's outputs. Raw forcing is retained separately, including across radiation intervals, so a refined level never interpolates an already corrected parent flux.

Member Enumeration Documentation

Diagnostic

enum TopographicShading::Diagnostic

Enumerator
Shadow
DirectFactor
Incidence
UnshadedSW
TerrainSW
NDiag

{ Shadow, DirectFactor, Incidence, UnshadedSW, TerrainSW, NDiag };

GeometryField

enum TopographicShading::GeometryField

Enumerator
Height
MapX
MapY
Latitude
Longitude
CosAlpha
SinAlpha
NGeom

{ Height, MapX, MapY, Latitude, Longitude, CosAlpha, SinAlpha, NGeom };

RadiationField

enum TopographicShading::RadiationField

Enumerator
CosZen
Shortwave
DirectVis
DirectNir
DiffuseVis
DiffuseNir
Longwave
NRad

{ CosZen, Shortwave, DirectVis, DirectNir, DiffuseVis, DiffuseNir, Longwave, NRad };

Member Function Documentation

apply()

void TopographicShading::apply	(	int	lev,
		const TopographicShadingChoice &	choice,
		const amrex::Vector< amrex::MultiFab * > &	destinations
	)

{
    AMREX_ALWAYS_ASSERT(geometry_ready(lev) && forcing_ready(lev));
    auto& l = *m_levels.at(lev);
    const Box dom = l.geom.Domain();
    const int ilo = dom.smallEnd(0), ihi = dom.bigEnd(0);
    const int jlo = dom.smallEnd(1), jhi = dom.bigEnd(1);
    const bool px = l.geom.isPeriodic(0), py = l.geom.isPeriodic(1);
    const auto dx = l.geom.CellSizeArray();
    const double calday = l.calday, declin = l.declin;
    Array4<const Real> terrain;
    for (MFIter mfi(l.patch); mfi.isValid(); ++mfi) { terrain = l.patch.const_array(mfi); }
    for (MFIter mfi(l.raw); mfi.isValid(); ++mfi) {
        auto g = l.ground.const_array(mfi);
        auto r = l.raw.const_array(mfi);
        auto diag = l.diagnostic.array(mfi);
        GpuArray<Array4<Real>, NRad> out;
        for (int n = 0; n < NRad; ++n) {
            out[n] = n < int(destinations.size()) && destinations[n]
                       ? destinations[n]->array(mfi) : Array4<Real>{};
        }
        ParallelFor(mfi.validbox(), [=] AMREX_GPU_DEVICE (int i, int j, int k) {
            auto s = erf_topography::sun_direction(calday, declin, g(i,j,k,Latitude), g(i,j,k,Longitude));
            if (choice.fixed_mu > 0) {
                const Real az = choice.fixed_azimuth*PI/Real(180);
                const Real h = std::sqrt(Real(1)-choice.fixed_mu*choice.fixed_mu);
                s = {h*std::sin(az), h*std::cos(az), choice.fixed_mu};
            }
            const Real east = s.x, north = s.y;
            s.x = g(i,j,k,CosAlpha)*east + g(i,j,k,SinAlpha)*north;
            s.y = -g(i,j,k,SinAlpha)*east + g(i,j,k,CosAlpha)*north;
            const Real dxp = dx[0]/g(i,j,k,MapX), dyp = dx[1]/g(i,j,k,MapY);
            const int im = px ? erf_topography::wrap(i-1,ilo,ihi-ilo+1) : amrex::max(i-1,ilo);
            const int ip = px ? erf_topography::wrap(i+1,ilo,ihi-ilo+1) : amrex::min(i+1,ihi);
            const int jm = py ? erf_topography::wrap(j-1,jlo,jhi-jlo+1) : amrex::max(j-1,jlo);
            const int jp = py ? erf_topography::wrap(j+1,jlo,jhi-jlo+1) : amrex::min(j+1,jhi);
            const Real hx = (terrain(ip,j,0)-terrain(im,j,0))/(dxp*(px ? 2 : amrex::max(1,ip-im)));
            const Real hy = (terrain(i,jp,0)-terrain(i,jm,0))/(dyp*(py ? 2 : amrex::max(1,jp-jm)));
            // RRTMGP's direct flux uses its interval-averaged cosine. The
            // midpoint sun supplies direction only; scaling by its cosine
            // would amplify the averaged flux near sunrise or sunset.
            const Real mu0 = r(i,j,k,CosZen);
            const bool has_direct_sun = mu0 > Real(1.e-4);
            const bool midpoint_daytime = s.z > Real(1.e-4);
            const Real inc = midpoint_daytime ? erf_topography::incidence(hx,hy,s) : Real(0);
            const bool shadow = has_direct_sun && choice.topo_shading &&
                (!midpoint_daytime || erf_topography::shadowed(
                    i,j,terrain,ilo,ihi,jlo,jhi,px,py,dxp,dyp,s,choice.shadlen));
            const Real factor = shadow ? Real(0) :
                ((has_direct_sun && choice.slope_rad) ? inc/mu0 : Real(1));
            const Real direct = r(i,j,k,DirectVis) + r(i,j,k,DirectNir);
            const Real diffuse = r(i,j,k,DiffuseVis) + r(i,j,k,DiffuseNir);
            // Preserve the exact raw broadband flux when no adjustment is needed.
            const Real sw = factor == Real(1) ? r(i,j,k,Shortwave) :
                erf_topography::corrected_shortwave(direct,diffuse,factor);
            diag(i,j,k,Shadow) = shadow ? Real(1) : Real(0);
            diag(i,j,k,DirectFactor) = factor;
            diag(i,j,k,Incidence) = inc;
            diag(i,j,k,UnshadedSW) = r(i,j,k,Shortwave);
            diag(i,j,k,TerrainSW) = sw;
            for (int n = 0; n < NRad; ++n) {
                if (out[n]) {
                    out[n](i,j,0) = n == Shortwave ? sw :
                        ((n == DirectVis || n == DirectNir) ? factor*r(i,j,k,n) : r(i,j,k,n));
                }
            }
        });
    }
}

Here is the call graph for this function:

define_level()

void TopographicShading::define_level	(	int	lev,
		const amrex::BoxArray &	ba2d,
		const amrex::DistributionMapping &	dm,
		const amrex::Geometry &	geom
	)

{
    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(ba.isDisjoint(),
        "Topographic radiation requires horizontally disjoint column boxes; "
        "set amr.no_box_split_dir = 2");
    if (lev >= int(m_levels.size())) { m_levels.resize(lev+1); }
    m_levels[lev] = std::make_unique<Level>();
    auto& l = *m_levels[lev];
    l.geom = geom;
    l.ground.define(ba, dm, NGeom, 0);
    l.raw.define(ba, dm, NRad, 0);
    l.diagnostic.define(ba, dm, NDiag, 0);
    l.raw.setVal(missing);
    l.diagnostic.setVal(-999);
    for (int n = 0; n < NRad; ++n) {
        l.aliases.push_back(std::make_unique<MultiFab>(l.raw, make_alias, n, 1));
    }
    // Descendant caches include this level's terrain outside their own grids.
    for (int f = lev+1; f < int(m_levels.size()); ++f) {
        if (m_levels[f]) { m_levels[f]->geometry_ready = false; }
    }
}

Here is the call graph for this function:

diagnostics()

const MultiFab * TopographicShading::diagnostics

(

int

lev

)

const

{
    return lev < int(m_levels.size()) && m_levels[lev] ? &m_levels[lev]->diagnostic : nullptr;
}

forcing_ready()

bool TopographicShading::forcing_ready

(

int

lev

)

const

{
    return lev < int(m_levels.size()) && m_levels[lev] && m_levels[lev]->forcing_ready;
}

geometry_ready()

bool TopographicShading::geometry_ready

(

int

lev

)

const

{
    return lev < int(m_levels.size()) && m_levels[lev] && m_levels[lev]->geometry_ready;
}

initialize_geometry()

void TopographicShading::initialize_geometry	(	int	lev,
		const TopographicShadingChoice &	choice,
		const amrex::MultiFab &	z_nd,
		const amrex::MultiFab &	map_x,
		const amrex::MultiFab &	map_y,
		const amrex::MultiFab *	lat,
		const amrex::MultiFab *	lon
	)

{
    auto& l = *m_levels.at(lev);
    const Box dom = makeSlab(l.geom.Domain(), 2, 0);
    const auto dx = l.geom.CellSizeArray();
    l.ground.setVal(missing);
    if (lev > 0) {
        AMREX_ALWAYS_ASSERT(geometry_ready(lev-1));
        prolong(l.ground, m_levels[lev-1]->ground,
                horizontal_ratio(m_levels[lev-1]->geom, l.geom), NGeom);
    }
    l.ground.ParallelCopy(map_x, 0, MapX, 1);
    l.ground.ParallelCopy(map_y, 0, MapY, 1);
    l.ground.setVal(choice.latitude, Latitude, 1);
    l.ground.setVal(choice.longitude, Longitude, 1);
    l.ground.setVal(1, CosAlpha, 1);
    l.ground.setVal(0, SinAlpha, 1);
    if (lat && lon) {
        l.ground.ParallelCopy(*lat, 0, Latitude, 1);
        l.ground.ParallelCopy(*lon, 0, Longitude, 1);
        MultiFab c(l.ground, make_alias, CosAlpha, 1);
        MultiFab s(l.ground, make_alias, SinAlpha, 1);
        fill_grid_rotation(c, s, *lat, *lon, l.geom);
    }
    const int kbot = l.geom.Domain().smallEnd(2);
    for (MFIter mfi(l.ground); mfi.isValid(); ++mfi) {
        // An elevated nested box has no ground of its own. Use the ancestor's.
        if (z_nd.boxArray()[mfi.index()].smallEnd(2) != kbot) { continue; }
        const auto z = z_nd.const_array(mfi);
        const auto g = l.ground.array(mfi);
        ParallelFor(mfi.validbox(), [=] AMREX_GPU_DEVICE (int i, int j, int k) {
            g(i,j,k,Height) = Real(0.25)*(z(i,j,kbot) + z(i+1,j,kbot)
                                        + z(i,j+1,kbot) + z(i+1,j+1,kbot));
        });
    }
    AMREX_ALWAYS_ASSERT(l.ground.min(Height) > missing);
    AMREX_ALWAYS_ASSERT(l.ground.min(MapX) > 0 && l.ground.min(MapY) > 0);
 
    // One search patch per owning rank, enlarged to include every possible ray
    // and the slope stencil. ParallelCopy brings remote terrain across MPI
    // boundaries once, at initialization/regrid, rather than at every sun update.
    // This patch bounds all boxes owned by the rank, not just their individual
    // ray neighborhoods. SFC or knapsack placement can assign distant boxes to
    // one rank, making its patch nearly the full 2D domain and increasing cache
    // memory and copy traffic even when shadlen is small. apply() currently reads
    // one contiguous terrain patch per rank; a per-box cache would need a
    // corresponding change there.
    const auto& ba = l.ground.boxArray();
    const auto& dm = l.ground.DistributionMap();
    BoxList boxes;
    Vector<int> owners;
    const Real max_map[] = {l.ground.max(MapX), l.ground.max(MapY)};
    for (int rank = 0; rank < ParallelDescriptor::NProcs(); ++rank) {
        Box patch;
        for (int b = 0; b < int(ba.size()); ++b) {
            if (dm[b] == rank) { patch.minBox(ba[b]); }
        }
        if (!patch.ok()) { continue; }
        for (int d = 0; d < 2; ++d) {
            const Real cells = choice.topo_shading ? choice.shadlen*max_map[d]/dx[d] : 0;
            const int ng = int(std::min(Real(dom.length(d)), std::ceil(cells))) + 2;
            patch.grow(d, ng);
            // Wrapped coordinates must be addressable in this one patch.
            if (l.geom.isPeriodic(d) &&
                (patch.smallEnd(d) < dom.smallEnd(d) || patch.bigEnd(d) > dom.bigEnd(d))) {
                patch.setSmall(d, dom.smallEnd(d));
                patch.setBig(d, dom.bigEnd(d));
            }
        }
        boxes.push_back(patch & dom);
        owners.push_back(rank);
    }
    l.patch.define(BoxArray(std::move(boxes)), DistributionMapping(std::move(owners)), 1, 0);
    l.patch.setVal(missing);
    for (int c = 0; c < lev; ++c) {
        AMREX_ALWAYS_ASSERT(geometry_ready(c));
        prolong(l.patch, m_levels[c]->ground, horizontal_ratio(m_levels[c]->geom, l.geom), 1);
    }
    l.patch.ParallelCopy(l.ground, Height, 0, 1);
    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(l.patch.min(0) > missing,
        "Topographic search patch contains terrain not supplied by any AMR level");
    l.geometry_ready = true;
}

Here is the call graph for this function:

interpolate_from_coarse()

void TopographicShading::interpolate_from_coarse

(

int

lev

)

{
    AMREX_ALWAYS_ASSERT(lev > 0 && forcing_ready(lev-1));
    auto& l = *m_levels.at(lev);
    const auto& c = *m_levels[lev-1];
    prolong(l.raw, c.raw, horizontal_ratio(c.geom, l.geom), NRad);
    set_sun(lev, c.calday, c.declin);
}

Here is the call graph for this function:

raw()

MultiFab * TopographicShading::raw

(

int

lev

)

{ return &m_levels.at(lev)->raw; }

raw_fields()

Vector< MultiFab * > TopographicShading::raw_fields

(

int

lev

)

{
    Vector<MultiFab*> fields;
    for (auto& f : m_levels.at(lev)->aliases) { fields.push_back(f.get()); }
    return fields;
}

set_sun()

void TopographicShading::set_sun	(	int	lev,
		double	calday,
		double	declin
	)

{
    auto& l = *m_levels.at(lev);
    l.calday = calday;
    l.declin = declin;
    l.forcing_ready = true;
}

Member Data Documentation

m_levels

amrex::Vector<std::unique_ptr<Level> > TopographicShading::m_levels

private

---

The documentation for this class was generated from the following files:

• Source/Radiation/ERF_TopographicShading.H
• Source/Radiation/ERF_TopographicShading.cpp