erf_grid_utils Namespace Reference

Classes
struct	UniformGridMetadata
struct	GeometricStretch
	A geometrically stretched vertical grid: dz_k = dz0 * ratio^k. More...
struct	InterpolationStencil

Functions
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real	comparison_tolerance (amrex::Real lhs, amrex::Real rhs) noexcept
bool	nearly_equal (amrex::Real lhs, amrex::Real rhs)
std::string	validate_uniform_axis (const amrex::Vector< amrex::Real > &coordinates, int expected_size, const std::string &axis_name, const std::string &field_description, amrex::Real &origin, amrex::Real &spacing)
std::string	validate_matching_grid (const UniformGridMetadata &reference, const UniformGridMetadata &candidate, const std::string &reference_description, const std::string &candidate_description)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE GeometricStretch	geometric_stretch (amrex::Real dz0, amrex::Real z_top, int nz, amrex::Real tol) noexcept
	Solve for the geometric stretch factor that fills a domain of depth z_top with nz layers whose first is dz0 thick. More...
GeometricStretch	build_geometric_zlevels (amrex::Vector< amrex::Real > &zlevels, amrex::Real dz0, amrex::Real z_top, amrex::Real tol)
	Fill a staggered z-level array with a geometrically stretched grid. More...
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE InterpolationStencil	uniform_interpolation_stencil (amrex::Real coordinate, amrex::Real origin, amrex::Real spacing, int point_count) noexcept

Function Documentation

build_geometric_zlevels()

GeometricStretch erf_grid_utils::build_geometric_zlevels ( amrex::Vector< amrex::Real > &  zlevels, amrex::Real  dz0, amrex::Real  z_top, amrex::Real  tol  )

inline

Fill a staggered z-level array with a geometrically stretched grid.

zlevels must already be sized to nz+1; its length is what sets nz. The top level is assigned z_top exactly rather than by accumulation, so the array closes on the domain top to the last bit even though the sum that got there is only accurate to tol.

Parameters

[in,out]	zlevels	staggered z-levels, sized nz+1 on entry and filled on exit
[in]	dz0	requested thickness of the first layer
[in]	z_top	depth of the column to be filled
[in]	tol	convergence tolerance passed to geometric_stretch

Returns

the stretch that was used, so the caller can report or check it

{
    const int nz = static_cast<int>(zlevels.size()) - 1;
    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(nz >= 1,
        "build_geometric_zlevels needs a zlevels array of at least two entries");
 
    const GeometricStretch gs = geometric_stretch(dz0, z_top, nz, tol);
 
    amrex::Real dz = gs.dz0;
    zlevels[0] = amrex::Real(0.0);
    for (int k(1); k<=nz; ++k) {
        zlevels[k] = zlevels[k-1] + dz;
        dz *= gs.ratio;
    }
    zlevels[nz] = z_top;
 
    return gs;
}

Referenced by init_zlevels().

Here is the call graph for this function:

Here is the caller graph for this function:

comparison_tolerance()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real erf_grid_utils::comparison_tolerance ( amrex::Real  lhs, amrex::Real  rhs  )

noexcept

{
    amrex::Real scale = amrex::Real(1.0);
    const amrex::Real lhs_abs = amrex::Math::abs(lhs);
    const amrex::Real rhs_abs = amrex::Math::abs(rhs);
    if (lhs_abs > scale) scale = lhs_abs;
    if (rhs_abs > scale) scale = rhs_abs;
    // NetCDF coordinate variables are commonly stored as float32 even when
    // ERF is built with double precision.  A few float32 ulps at the
    // coordinate magnitude are therefore source-file quantization, not grid
    // nonuniformity.  The Real tolerance is retained for native Real data;
    // the float tolerance is deliberately a small ulp-based allowance rather
    // than an arbitrary absolute threshold.
    const amrex::Real real_ulp =
        amrex::Real(64.0) * std::numeric_limits<amrex::Real>::epsilon() * scale;
    const amrex::Real float_ulp =
        amrex::Real(8.0) * static_cast<amrex::Real>(std::numeric_limits<float>::epsilon()) * scale;
    return (real_ulp > float_ulp) ? real_ulp : float_ulp;
}

Referenced by nearly_equal(), uniform_interpolation_stencil(), and validate_uniform_axis().

Here is the caller graph for this function:

geometric_stretch()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE GeometricStretch erf_grid_utils::geometric_stretch ( amrex::Real  dz0, amrex::Real  z_top, int  nz, amrex::Real  tol  )

noexcept

Solve for the geometric stretch factor that fills a domain of depth z_top with nz layers whose first is dz0 thick.

The layer thicknesses are dz0 * ratio^k, so the column depth is the geometric sum

S(ratio) = dz0 * (ratio^nz - 1) / (ratio - 1),

and the factor wanted is the root of S(ratio) = z_top, found by Newton iteration from ratio = 1.03. The derivative is analytic, so the iteration converges in a handful of passes for any physically sensible (dz0, z_top, nz).

This is the natural construction for an idealized case, where dz0 and z_top are what the user knows a priori and the vertical distribution between them is free. It is not the right construction when a file already carries a vertical distribution worth keeping – see wrf_mean_layer_thickness in ERF_InitFromWRFInput.cpp – because one constant factor cannot represent a stretch ratio that varies with height.

A dz0 at or above the uniform spacing z_top/nz cannot be stretched to fit, since the thicknesses only grow: that case returns the uniform grid with uniform = true, rather than a factor below one that would thin the layers aloft. A single layer is the same situation for a different reason – the sum of one layer is dz0 whatever the factor is, so S has no root and its derivative is identically zero – and is answered the same way.

Parameters

[in]	dz0	requested thickness of the first layer
[in]	z_top	depth of the column to be filled
[in]	nz	number of layers
[in]	tol	convergence tolerance on the column depth, in the same units as z_top

Returns

the stretch factor, the first-layer thickness used, and whether it converged

{
    GeometricStretch gs;
 
    const amrex::Real Nz = static_cast<amrex::Real>(nz);
    if (nz < 1 || !(z_top > amrex::Real(0.0)) || !(dz0 > amrex::Real(0.0))) {
        return gs;
    }
 
    // Default to the uniform grid, which is also the answer whenever the requested first
    // layer is already at least as thick as the uniform spacing, and the only grid there
    // is when the column holds a single layer.
    if (nz == 1 || dz0 >= z_top/Nz) {
        gs.dz0       = z_top/Nz;
        gs.ratio     = amrex::Real(1.0);
        gs.residual  = amrex::Real(0.0);
        gs.uniform   = true;
        gs.converged = true;
        return gs;
    }
 
    amrex::Real SFact = amrex::Real(1.03);
    amrex::Real F     = dz0 * ( (std::pow(SFact,Nz) - amrex::Real(1.0)) / (SFact - amrex::Real(1.0)) ) - z_top;
 
    const int max_iter = 50;
    int iter = 0;
    while (amrex::Math::abs(F) > tol && iter < max_iter) {
        const amrex::Real dFdSF = dz0 * ( Nz * std::pow(SFact,Nz-amrex::Real(1.0)) * (SFact - amrex::Real(1.0))
                                        - std::pow(SFact,Nz) + amrex::Real(1.0) ) /
                                        std::pow(SFact-amrex::Real(1.0),amrex::Real(2.0));
        SFact -= F/dFdSF;
        const amrex::Real SFact_min = amrex::Real(1.0) + tol;
        SFact  = (SFact > SFact_min) ? SFact : SFact_min;
        F      = dz0 * ( (std::pow(SFact,Nz) - amrex::Real(1.0)) / (SFact - amrex::Real(1.0)) ) - z_top;
        ++iter;
    }
 
    gs.dz0       = dz0;
    gs.ratio     = SFact;
    gs.residual  = F;
    gs.uniform   = false;
    gs.converged = (amrex::Math::abs(F) <= tol);
    return gs;
}

Referenced by build_geometric_zlevels().

Here is the caller graph for this function:

nearly_equal()

bool erf_grid_utils::nearly_equal ( amrex::Real  lhs, amrex::Real  rhs  )

inline

{
    return std::abs(lhs - rhs) <= comparison_tolerance(lhs, rhs);
}

Referenced by validate_matching_grid().

Here is the call graph for this function:

Here is the caller graph for this function:

uniform_interpolation_stencil()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE InterpolationStencil erf_grid_utils::uniform_interpolation_stencil ( amrex::Real  coordinate, amrex::Real  origin, amrex::Real  spacing, int  point_count  )

noexcept

{
    InterpolationStencil stencil;
    if (point_count < 2 || spacing <= amrex::Real(0.0)) {
        return stencil;
    }
 
    const amrex::Real upper = origin + static_cast<amrex::Real>(point_count - 1) * spacing;
    const amrex::Real lower_tolerance = comparison_tolerance(coordinate, origin);
    const amrex::Real upper_tolerance = comparison_tolerance(coordinate, upper);
 
    if (coordinate < origin - lower_tolerance || coordinate > upper + upper_tolerance) {
        return stencil;
    }
 
    const amrex::Real bounded_coordinate =
        (coordinate <= origin + lower_tolerance) ? origin :
        ((coordinate >= upper - upper_tolerance) ? upper : coordinate);
    const amrex::Real fractional_index = (bounded_coordinate - origin) / spacing;
    int lower = static_cast<int>(amrex::Math::floor(fractional_index));
    lower = (lower < point_count - 1) ? lower : point_count - 2;
 
    amrex::Real weight = fractional_index - static_cast<amrex::Real>(lower);
    weight = (weight < amrex::Real(0.0)) ? amrex::Real(0.0) :
             ((weight > amrex::Real(1.0)) ? amrex::Real(1.0) : weight);
 
    stencil.lower = lower;
    stencil.weight = weight;
    stencil.inside = true;
    return stencil;
}

Referenced by ForestDrag::define_drag_field().

Here is the call graph for this function:

Here is the caller graph for this function:

validate_matching_grid()

std::string erf_grid_utils::validate_matching_grid ( const UniformGridMetadata &  reference, const UniformGridMetadata &  candidate, const std::string &  reference_description, const std::string &  candidate_description  )

inline

{
    if (candidate.nx != reference.nx || candidate.ny != reference.ny) {
        return candidate_description + " dimensions (" + std::to_string(candidate.nx) +
               " x " + std::to_string(candidate.ny) + ") do not match " +
               reference_description + " dimensions (" + std::to_string(reference.nx) +
               " x " + std::to_string(reference.ny) + ")";
    }
 
    struct MetadataValue {
        const char* name;
        amrex::Real reference_value;
        amrex::Real candidate_value;
    };
    const MetadataValue values[] = {
        {"x spacing", reference.dx, candidate.dx},
        {"y spacing", reference.dy, candidate.dy},
        {"x origin", reference.xmin, candidate.xmin},
        {"y origin", reference.ymin, candidate.ymin}
    };
    for (const auto& value : values) {
        if (!nearly_equal(value.reference_value, value.candidate_value)) {
            std::ostringstream message;
            message << candidate_description << " " << value.name << " ("
                    << value.candidate_value << ") does not match "
                    << reference_description << " (" << value.reference_value << ")";
            return message.str();
        }
    }
 
    return {};
}

Referenced by ForestDrag::ForestDrag().

Here is the call graph for this function:

Here is the caller graph for this function:

validate_uniform_axis()

std::string erf_grid_utils::validate_uniform_axis ( const amrex::Vector< amrex::Real > &  coordinates, int  expected_size, const std::string &  axis_name, const std::string &  field_description, amrex::Real &  origin, amrex::Real &  spacing  )

inline

{
    if (expected_size < 2) {
        return field_description + " requires at least two " + axis_name +
               " coordinates for bilinear interpolation";
    }
    if (coordinates.size() != static_cast<amrex::Long>(expected_size)) {
        return field_description + " has " + axis_name + " coordinate length " +
               std::to_string(coordinates.size()) + ", expected " +
               std::to_string(expected_size);
    }
 
    origin = coordinates.front();
    // Derive the spacing from the two endpoints rather than from the first
    // interval alone.  The returned origin/spacing pair is what
    // uniform_interpolation_stencil uses to reconstruct the axis as
    // origin + index*spacing, so a first interval that is rounded slightly
    // short or long would otherwise accumulate that error over the whole axis
    // -- for a float32 axis with many points the reconstructed far end can
    // drift by (point_count-1) ulps, which is large enough to shift the
    // stencil by a cell near the upper edge.  Averaging over the full extent
    // pins both endpoints instead.  A non-finite last coordinate propagates
    // into spacing and is caught by the isfinite check below.
    spacing = (coordinates[expected_size-1] - coordinates.front()) /
              static_cast<amrex::Real>(expected_size - 1);
    if (!std::isfinite(origin) || !std::isfinite(spacing) || spacing <= amrex::Real(0.0)) {
        return field_description + " must have finite, strictly increasing " +
               axis_name + " coordinates with nonzero spacing";
    }
 
    for (int index = 1; index < expected_size; ++index) {
        const amrex::Real coordinate = coordinates[index];
        if (!std::isfinite(coordinate) || coordinate <= coordinates[index-1]) {
            return field_description + " must have finite, strictly increasing " +
                   axis_name + " coordinates";
        }
        // Compare each interval against the mean spacing instead of
        // accumulating the first spacing over the whole axis.  This avoids
        // rejecting a valid float32 axis because each stored coordinate was
        // rounded independently, while still rejecting materially nonuniform
        // grids.
        //
        // The tolerance is the *larger* of the two allowances, not the
        // smaller: the ulp term covers source-file quantization, which scales
        // with the coordinate magnitude, while the fraction-of-spacing term
        // covers grids whose coordinates are small enough that a few ulps are
        // negligible.  Taking the minimum cancels the ulp allowance in exactly
        // the case it exists for -- a fine grid at a large offset.  Concretely,
        // at a UTM easting near 5e5 the float32 ulp is 0.03125 m, so a 0.4 m
        // raster stores intervals that deviate from the nominal spacing by up
        // to 0.025 m: more than 5% of the spacing, yet far below the ulp
        // allowance.  Integer-metre spacings are unaffected either way, since
        // integers are exactly representable in float32 at UTM magnitudes.
        const amrex::Real local_spacing = coordinate - coordinates[index-1];
        const amrex::Real spacing_tolerance = std::max(
            amrex::Real(0.05) * std::abs(spacing),
            comparison_tolerance(coordinate, coordinates[index-1]));
        if (std::abs(local_spacing - spacing) > spacing_tolerance) {
            return field_description + " has nonuniform " + axis_name +
                   " coordinates, which are not supported";
        }
    }
 
    return {};
}

Here is the call graph for this function: