ERF_NodalReconstruction.H File Reference

#include <algorithm>
#include <cmath>
#include <limits>
#include <string>
#include <utility>
#include <vector>
#include "AMReX.H"
#include "AMReX_Array4.H"
#include "AMReX_BLassert.H"
#include "AMReX_Box.H"
#include "AMReX_FArrayBox.H"
#include "AMReX_GpuDevice.H"
#include "AMReX_Geometry.H"
#include "AMReX_MultiFab.H"
#include "AMReX_ParallelDescriptor.H"
#include "AMReX_Print.H"
#include "AMReX_REAL.H"
#include "ERF_NumericalConstants.H"

Include dependency graph for ERF_NodalReconstruction.H:

This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Classes
struct	SolveInfo
	Convergence and diagnostic information for the reconstruction solve. More...
class	NodalReconstruction

Enumerations
enum class	VariationOperator { FirstDeriv , Laplacian }

Functions
void	deaverage_line (amrex::Real const *t, int n, amrex::Real *s)
std::pair< amrex::FArrayBox, SolveInfo >	deaverage_nodal_slice (amrex::Box const &w_face_box, amrex::FArrayBox const &z_w_face, amrex::FArrayBox const &z_ref)
amrex::FArrayBox	reconstruct_nodal_height_slice (amrex::Box const &w_face_slice_box, amrex::Geometry const &geom, amrex::FArrayBox const &z_w_face_slice, int klev, std::string const &src_name)
void	fill_nodal_level_from_slice (amrex::MultiFab &z_phys_nd, int klev, amrex::FArrayBox const &z_nd_slice)

Enumeration Type Documentation

VariationOperator

enum VariationOperator

strong

Enumerator
FirstDeriv
Laplacian

{
    FirstDeriv,
    Laplacian
};

Function Documentation

deaverage_line()

void deaverage_line ( amrex::Real const *  t, int  n, amrex::Real *  s  )

inline

Exactly de-average one line.

Given cell values t(0..n-1), produce node values s(0..n) satisfying

1/2 ( s(k) + s(k+1) ) = t(k),    k = 0 .. n-1

exactly. The system is underdetermined by exactly one degree of freedom: the recursion s(k+1) = 2 t(k) - s(k) solves it for any seed, and changing the seed adds the sawtooth (-1)^k, which the averaging operator annihilates. So the seed is free, and it is chosen here in closed form to minimize the first-difference energy sum_k ( s(k+1) - s(k) )^2 along the line, i.e. to leave as little of that sawtooth in the answer as the constraint permits.

Parameters

[in]	t	cell values, length n
[in]	n	number of cells
[out]	s	node values, length n+1

{
    using amrex::Real;
 
    // Particular solution, seeded at zero.
    s[0] = zero;
    for (int k(0); k<n; ++k) { s[k+1] = two*t[k] - s[k]; }
 
    // With s(k) = (-1)^k s0 + p(k) we have s(k+1) - s(k) = -2 (-1)^k s0 + dp(k),
    //    and minimizing the sum of squares of that over s0 gives the closed form
    //    s0 = sum_k (-1)^k dp(k) / (2n).
    Real acc(zero), sgn(one);
    for (int k(0); k<n; ++k) { acc += sgn * (s[k+1] - s[k]); sgn = -sgn; }
    Real const s0 = acc / (two * static_cast<Real>(n));
 
    sgn = one;
    for (int k(0); k<=n; ++k) { s[k] += sgn * s0; sgn = -sgn; }
}

Referenced by deaverage_nodal_slice().

Here is the caller graph for this function:

deaverage_nodal_slice()

std::pair<amrex::FArrayBox,SolveInfo> deaverage_nodal_slice ( amrex::Box const &  w_face_box, amrex::FArrayBox const &  z_w_face, amrex::FArrayBox const &  z_ref  )

inline

Exactly de-average one horizontal slice: find nodal S with avg4(S) = T to round-off, rather than the closest such field a smoothness prior will allow.

The four-node averaging operator is a tensor product of two one-dimensional two-point averages, A = Ax (x) Ay. Each of those is surjective – the line recursion above always solves – so A is surjective and an exact solution always exists; the system is underdetermined, not inconsistent, with (nx+1)(ny+1) unknowns against nx*ny equations. Separability means it can be solved in two sweeps: first Ay U = T for U (cell-centered in x, nodal in y), then Ax S = U. Both are O(N) and non-iterative.

What cannot be removed is grid-scale structure. The null space is spanned by (-1)^i a(j) + (-1)^j b(i), in which a does not vary with i, so the amplitude of the sawtooth cannot be varied along a line; rough terrain anywhere on a line leaves a signature everywhere on it. Minimizing the first-difference energy line by line, as deaverage_line does, exhausts that freedom – a globally coupled optimization over the whole null space was measured to give no further improvement. The caller is therefore responsible for deciding whether the result is usable as terrain, and for falling back if it is not.

Parameters

[in]	w_face_box	box of the w-face slice
[in]	z_w_face	w-face heights on w_face_box
[in]	z_ref	direct interpolation of z_w_face to the nodes, for diagnostics only

Returns

the nodal heights and the diagnostics describing them

{
    using amrex::Real;
 
    amrex::Box nd_box = amrex::surroundingNodes(w_face_box);
    nd_box.setSmall(2,0);
    nd_box.setBig(2,0);
 
    amrex::FArrayBox S(nd_box, 1, amrex::The_Managed_Arena());
 
    int const ilo = w_face_box.smallEnd(0), ihi = w_face_box.bigEnd(0);
    int const jlo = w_face_box.smallEnd(1), jhi = w_face_box.bigEnd(1);
    int const nx  = ihi - ilo + 1;
    int const ny  = jhi - jlo + 1;
 
    auto const T_arr = z_w_face.const_array();
    auto const R_arr = z_ref.const_array();
    auto const S_arr = S.array();
 
    // U is the x-average of S, so it is cell-centered in x and nodal in y.
    std::vector<Real> U(static_cast<std::size_t>(nx) * static_cast<std::size_t>(ny+1));
    std::vector<Real> tline(static_cast<std::size_t>(std::max(nx,ny)));
    std::vector<Real> sline(static_cast<std::size_t>(std::max(nx,ny)+1));
 
    // Sweep 1: Ay U = T, one line in y per x-cell.
    for (int i(0); i<nx; ++i) {
        for (int j(0); j<ny; ++j) { tline[static_cast<std::size_t>(j)] = T_arr(ilo+i,jlo+j,0); }
        deaverage_line(tline.data(), ny, sline.data());
        for (int j(0); j<=ny; ++j) {
            U[static_cast<std::size_t>(i)*static_cast<std::size_t>(ny+1) +
              static_cast<std::size_t>(j)] = sline[static_cast<std::size_t>(j)];
        }
    }
 
    // Sweep 2: Ax S = U, one line in x per y-node.
    for (int j(0); j<=ny; ++j) {
        for (int i(0); i<nx; ++i) {
            tline[static_cast<std::size_t>(i)] =
                U[static_cast<std::size_t>(i)*static_cast<std::size_t>(ny+1) +
                  static_cast<std::size_t>(j)];
        }
        deaverage_line(tline.data(), nx, sline.data());
        for (int i(0); i<=nx; ++i) { S_arr(ilo+i,jlo+j,0) = sline[static_cast<std::size_t>(i)]; }
    }
 
    SolveInfo info;
    info.converged      = true;
    info.iterations     = 0;
    info.refinements    = 0;
    info.regularization = zero;
    info.deviation_cap  = zero;
    info.min_value      =  std::numeric_limits<Real>::max();
    info.max_value      = -std::numeric_limits<Real>::max();
 
    for (int j(jlo); j<=jhi+1; ++j) {
        for (int i(ilo); i<=ihi+1; ++i) {
            Real const v = S_arr(i,j,0);
            info.min_value = std::min(info.min_value, v);
            info.max_value = std::max(info.max_value, v);
            info.deviation = std::max(info.deviation, std::abs(v - R_arr(i,j,0)));
        }
    }
    for (int j(jlo); j<=jhi; ++j) {
        for (int i(ilo); i<=ihi; ++i) {
            Real const a4 = fourth * ( S_arr(i,j,0) + S_arr(i+1,j,0)
                                     + S_arr(i,j+1,0) + S_arr(i+1,j+1,0) );
            Real const r4 = fourth * ( R_arr(i,j,0) + R_arr(i+1,j,0)
                                     + R_arr(i,j+1,0) + R_arr(i+1,j+1,0) );
            info.max_average_error    = std::max(info.max_average_error,
                                                 std::abs(a4 - T_arr(i,j,0)));
            info.interp_average_error = std::max(info.interp_average_error,
                                                 std::abs(r4 - T_arr(i,j,0)));
        }
    }
    info.final_residual = info.max_average_error;
 
    return {std::move(S), info};
}

Referenced by reconstruct_nodal_height_slice().

Here is the call graph for this function:

Here is the caller graph for this function:

fill_nodal_level_from_slice()

void fill_nodal_level_from_slice ( amrex::MultiFab &  z_phys_nd, int  klev, amrex::FArrayBox const &  z_nd_slice  )

inline

Copy one reconstructed nodal slice into level klev of a nodal MultiFab, filling the lateral ghost nodes by clamping to the nearest valid node.

Parameters

[in,out]	z_phys_nd	Nodal MultiFab of physical heights.
[in]	klev	Vertical node index to fill.
[in]	z_nd_slice	Nodal slice, indexed at k = 0, spanning the whole domain in (x,y).

{
    amrex::Box const& nd_box = z_nd_slice.box();
    int const ilo = nd_box.smallEnd(0); int const ihi = nd_box.bigEnd(0);
    int const jlo = nd_box.smallEnd(1); int const jhi = nd_box.bigEnd(1);
 
    auto const z_slice_arr = z_nd_slice.const_array();
 
    for ( amrex::MFIter mfi(z_phys_nd); mfi.isValid(); ++mfi ) {
        amrex::Box gbx = mfi.growntilebox();
 
        // Boxes that do not reach this level get their surface data from the
        // ParallelCopy in make_terrain_fitted_coords
        if (klev < gbx.smallEnd(2) || klev > gbx.bigEnd(2)) { continue; }
 
        amrex::Box sbx = amrex::makeSlab(gbx, 2, klev);
        auto const z_arr = z_phys_nd.array(mfi);
        amrex::ParallelFor(sbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
        {
            int const ii = amrex::max(amrex::min(i,ihi),ilo);
            int const jj = amrex::max(amrex::min(j,jhi),jlo);
            z_arr(i,j,k) = z_slice_arr(ii,jj,0);
        });
    }
 
    // The kernels above read z_nd_slice, which every caller destroys as soon as
    // this returns.  Freeing it while the launch is still queued would hand its
    // bytes back to the arena for reuse under a running kernel.
    amrex::Gpu::streamSynchronize();
}

Here is the call graph for this function:

reconstruct_nodal_height_slice()

amrex::FArrayBox reconstruct_nodal_height_slice ( amrex::Box const &  w_face_slice_box, amrex::Geometry const &  geom, amrex::FArrayBox const &  z_w_face_slice, int  klev, std::string const &  src_name  )

inline

Reconstruct one horizontal slice of nodal heights from heights that live at w-faces, reporting diagnostics and aborting if the result is unusable as terrain.

Both NetCDF initialization paths need exactly this operation and differ only in how the w-face slice is assembled (geopotential from PH + PHB for wrfinput, HGT_M for metgrid) and in what they do with the answer, so the solve and its checks live here and are called from both.

The solve is serial and covers the whole domain, so it runs on the I/O rank only and the resulting nodal slice is broadcast. Running it redundantly on every rank would cost a full-domain CG solve and several full-domain managed allocations per rank, i.e. a startup stall and a memory spike that grow with the rank count rather than with the decomposition. z_w_face_slice must be host-accessible and must hold the global slice on the I/O rank.

This is collective: every rank must call it.

Parameters

[in]	w_face_slice_box	Box covering the domain at the w-face staggering, a single k = 0 slab.
[in]	geom	Geometry of the level, used for the grid spacing.
[in]	z_w_face_slice	Heights on w_face_slice_box.
[in]	klev	Vertical index of the slice; used only in messages.
[in]	src_name	Name of the source data; used only in messages.

Returns

The nodal heights on surroundingNodes(w_face_slice_box), indexed at k = 0.

{
    // Fail here, with a message naming the data, rather than deep inside the
    // solver: an empty box means the source heights were never read, and the
    // degenerate loop bounds it produces would otherwise read out of bounds.
    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(w_face_slice_box.ok() && w_face_slice_box.length(2) == 1,
                                     "reconstruct_nodal_height_slice: the input slice is empty "
                                     "or spans more than one z level -- check that the source "
                                     "heights were read successfully.");
 
    // The slice is typically filled by a device kernel, and everything below
    // reads it from the host.
    amrex::Gpu::streamSynchronize();
 
    amrex::Box nd_box = amrex::surroundingNodes(w_face_slice_box);
    nd_box.setSmall(2, 0);
    nd_box.setBig(2, 0);
 
    // Broadcast target; also the return value, so it must be device-readable.
    amrex::FArrayBox z_nd_slice(nd_box, 1, amrex::The_Managed_Arena());
 
    int const ioproc = amrex::ParallelDescriptor::IOProcessorNumber();
 
    // Non-zero if the reconstruction is unusable; decided on the I/O rank and
    // broadcast so that every rank aborts together rather than one rank calling
    // MPI_Abort out from under the others.
    int bad_reconstruction = 0;
 
    if (amrex::ParallelDescriptor::IOProcessor())
    {
        // Range of the source data, used to decide whether a reconstruction is
        // usable as terrain.  Nodes may legitimately over/undershoot the w-face
        // values where the terrain is under-resolved, but only by a fraction of the
        // relief of the layer itself.
        amrex::Real src_min =  std::numeric_limits<amrex::Real>::max();
        amrex::Real src_max = -std::numeric_limits<amrex::Real>::max();
        {
            auto const src_arr = z_w_face_slice.const_array();
            amrex::LoopOnCpu(w_face_slice_box, [=,&src_min,&src_max] (int i, int j, int /*k*/) noexcept
            {
                src_min = std::min(src_min, src_arr(i,j,0));
                src_max = std::max(src_max, src_arr(i,j,0));
            });
        }
        amrex::Real const relief = std::max(src_max - src_min, amrex::Real(1.0));
        amrex::Real const slack  = std::max(amrex::Real(0.5) * relief, amrex::Real(10.0));
 
        auto in_range = [=] (SolveInfo const& si) noexcept -> bool
        {
            return std::isfinite(si.min_value) && std::isfinite(si.max_value) &&
                   (si.min_value >= src_min - slack) && (si.max_value <= src_max + slack);
        };
 
        NodalReconstruction NR_solver(w_face_slice_box, geom);
        amrex::FArrayBox z_ref = NR_solver.makeReference(z_w_face_slice);
 
        // Invert the four-node averaging operator exactly.  The operator is
        // surjective, so this is always possible; what it costs is grid-scale
        // structure in the nodal field, because the sawtooth it cannot control is
        // the operator's own null space.  See deaverage_nodal_slice.
        //
        // Exactness is what the callers actually need.  Each w-face height is the
        // average of the four nodes around it, and each cell-centered height is in
        // turn the average of the two w-faces bracketing it, so avg4(S) = T makes
        // both of them -- and hence what the vertical remap of the atmospheric state
        // searches against, and what the lateral pressure-gradient metric
        // differences -- agree with the source to round-off rather than to within a
        // regularization-dependent error.
        std::pair<amrex::FArrayBox,SolveInfo> result = deaverage_nodal_slice(w_face_slice_box,
                                                                            z_w_face_slice, z_ref);
        SolveInfo info = result.second;
        bool const exact_usable = in_range(info);
 
        if (exact_usable) {
            amrex::Print() << "Nodal reconstruction at k = " << klev << ": exact de-averaging"
                           << "\n    nodal heights in [" << info.min_value << ", " << info.max_value
                           << "] m vs " << src_name << " in [" << src_min << ", " << src_max << "] m"
                           << "\n    max |avg4(nodal) - " << src_name << "| = " << info.max_average_error
                           << " m (direct interpolation gives " << info.interp_average_error << " m)"
                           << "\n    max deviation from direct interpolation = " << info.deviation
                           << " m" << std::endl;
        } else {
            // Rough enough terrain can force an exact solution whose nodal heights are
            // nowhere near the source.  Fall back to the regularized least-squares fit,
            // which stays bounded but no longer reproduces the source exactly, and say so.
            amrex::Print() << "WARNING: exact de-averaging at k = " << klev
                           << " gives nodal heights in [" << info.min_value << ", " << info.max_value
                           << "] m, outside the range of the " << src_name << " ["
                           << src_min << ", " << src_max << "] m.\n"
                           << "         Falling back to the regularized least-squares fit, which "
                              "does NOT reproduce the " << src_name << " exactly; the w-face "
                              "heights will differ from the source by the error reported below.\n";
 
            amrex::Real const tol = amrex::Real(1.e-10);
            result = NR_solver.solve(z_w_face_slice, z_ref, VariationOperator::FirstDeriv, tol);
            info   = result.second;
 
            if (!info.converged) {
                amrex::Print() << "WARNING: Nodal reconstruction did not converge at k = " << klev
                               << "; residual is: " << info.final_residual
                               << " and requested tolerance was: " << tol << "\n";
            }
 
            amrex::Print() << "Nodal reconstruction at k = " << klev << ": "
                           << info.iterations << " CG iterations, " << info.refinements
                           << " refinements, regularization " << info.regularization
                           << "\n    nodal heights in [" << info.min_value << ", " << info.max_value
                           << "] m vs " << src_name << " in [" << src_min << ", " << src_max << "] m"
                           << "\n    max |avg4(nodal) - " << src_name << "| = " << info.max_average_error
                           << " m (direct interpolation gives " << info.interp_average_error << " m)"
                           << "\n    max deviation from direct interpolation = " << info.deviation
                           << " m (cap " << info.deviation_cap << " m)" << std::endl;
 
            if (!in_range(info)) { bad_reconstruction = 1; }
        }
 
        z_nd_slice.copy<amrex::RunOn::Host>(result.first, nd_box, 0, nd_box, 0, 1);
    }
 
    amrex::ParallelDescriptor::Bcast(&bad_reconstruction, 1, ioproc);
    if (bad_reconstruction) {
        amrex::Abort("Nodal reconstruction produced heights far outside the range of the "
                     + src_name + "; the reconstruction is not usable as terrain.");
    }
 
    amrex::ParallelDescriptor::Bcast(z_nd_slice.dataPtr(), z_nd_slice.size(), ioproc);
 
    return z_nd_slice;
}

Here is the call graph for this function: