ERF_NodalReconstruction.H File Reference

#include <algorithm>
#include <cmath>
#include <limits>
#include <string>
#include <utility>
#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_Constants.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
amrex::FArrayBox	reconstruct_nodal_height_slice (amrex::Box const &cc_slice_box, amrex::Geometry const &geom, amrex::FArrayBox const &z_cc_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

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 &  cc_slice_box, amrex::Geometry const &  geom, amrex::FArrayBox const &  z_cc_slice, int  klev, std::string const &  src_name  )

inline

Reconstruct one horizontal slice of nodal heights from heights that live at cell centers in (x,y), 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 cell-centered 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_cc_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]	cc_slice_box	Box covering the domain, cell-centered in (x,y) and a single k = 0 slab.
[in]	geom	Geometry of the level, used for the grid spacing.
[in]	z_cc_slice	Heights on cc_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(cc_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(cc_slice_box.ok() && cc_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(cc_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())
    {
        // Solve for node values that reproduce the cell-centered values as closely
        // as a bounded, smooth nodal field can.  We deliberately do *not* invert the
        // four-node averaging operator exactly: its symbol vanishes at the grid
        // Nyquist mode, so exact de-averaging amplifies grid-scale content of the
        // terrain without bound and returns nodal heights that are kilometers away
        // from the input terrain.  See the class documentation above for the
        // regularized least-squares formulation used instead.
        amrex::Real const tol = amrex::Real(1.e-10);
        NodalReconstruction NR_solver(cc_slice_box, geom);
        amrex::FArrayBox z_ref = NR_solver.makeReference(z_cc_slice);
        std::pair<amrex::FArrayBox,SolveInfo> result = NR_solver.solve(z_cc_slice, z_ref,
                                                                      VariationOperator::FirstDeriv, tol);
        SolveInfo const& 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";
        }
 
        // Range check on the reconstructed heights.  This is the check that has
        // teeth: max |avg4(nodal) - source| is close to satisfied by construction
        // and cannot detect a blown-up reconstruction.
        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_cc_slice.const_array();
            amrex::LoopOnCpu(cc_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::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;
 
        // Nodes may legitimately over/undershoot the cell values where the terrain
        // is under-resolved, but only by a fraction of the relief of the layer itself.
        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));
 
        if ( !std::isfinite(info.min_value) || !std::isfinite(info.max_value) ||
             (info.min_value < src_min - slack) || (info.max_value > src_max + slack) )
        {
            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: