erf_terrain_file Namespace Reference

Classes
struct	HeightTable

Functions
std::string	TableProblem (bool is_usgs, int nx, int ny, const Vector< Real > &x, const Vector< Real > &y, const Vector< Real > &z, const GpuArray< Real, AMREX_SPACEDIM > &prob_lo, const GpuArray< Real, AMREX_SPACEDIM > &prob_hi, bool require_domain_coverage)
std::string	TableSummary (const std::string &fname, bool is_usgs, int nx, int ny, const Vector< Real > &x, const Vector< Real > &y, const Vector< Real > &z)
const HeightTable &	ReadHeightTable (const std::string &fname, bool is_usgs, const GpuArray< Real, AMREX_SPACEDIM > &prob_lo, const GpuArray< Real, AMREX_SPACEDIM > &prob_hi, bool require_domain_coverage)
Array< Real, 4 >	FaceBandSlopes (const MultiFab &z_phys_nd, const MultiFab &mf_mx, const MultiFab &mf_my, const Geometry &geom, int band)
std::string	OpenFaceSlopeWarning (int lev, int band, const Array< Real, 4 > &slopes, const Array< bool, 4 > &open, const Array< std::string, 4 > &bc_names)
const HeightTable &	ReadHeightTable (const std::string &fname, bool is_usgs, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &prob_lo, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &prob_hi, bool require_domain_coverage=true)
std::string	TableProblem (bool is_usgs, int nx, int ny, const amrex::Vector< amrex::Real > &x, const amrex::Vector< amrex::Real > &y, const amrex::Vector< amrex::Real > &z, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &prob_lo, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &prob_hi, bool require_domain_coverage=true)
std::string	TableSummary (const std::string &fname, bool is_usgs, int nx, int ny, const amrex::Vector< amrex::Real > &x, const amrex::Vector< amrex::Real > &y, const amrex::Vector< amrex::Real > &z)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real	CellTerrainSlope (int i, int j, int k, const amrex::Array4< const amrex::Real > &z_nd, amrex::Real mfx, amrex::Real mfy, amrex::Real dxinv, amrex::Real dyinv) noexcept
amrex::Array< amrex::Real, 4 >	FaceBandSlopes (const amrex::MultiFab &z_phys_nd, const amrex::MultiFab &mf_mx, const amrex::MultiFab &mf_my, const amrex::Geometry &geom, int band)
std::string	OpenFaceSlopeWarning (int lev, int band, const amrex::Array< amrex::Real, 4 > &slopes, const amrex::Array< bool, 4 > &open, const amrex::Array< std::string, 4 > &bc_names)

Variables
constexpr amrex::Real	open_face_slope_warn = amrex::Real(0.3)

Function Documentation

CellTerrainSlope()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real erf_terrain_file::CellTerrainSlope ( int  i, int  j, int  k, const amrex::Array4< const amrex::Real > &  z_nd, amrex::Real  mfx, amrex::Real  mfy, amrex::Real  dxinv, amrex::Real  dyinv  )

noexcept

Physical slope |grad h| of the terrain over cell (i,j,k): the terrain drops of ComputeTerrainCellDrops (the same drops as the diffusive limits and the slope factor) times the map factor over the cell size in each direction.

Parameters

[in]	i,j,k	cell indices
[in]	z_nd	nodal heights
[in]	mfx,mfy	map factors of the cell
[in]	dxinv,dyinv	inverse horizontal cell sizes

{
    const TerrainCellDrops d = ComputeTerrainCellDrops(i, j, k, z_nd);
    const amrex::Real hx = mfx * dxinv * d.dx_drop;
    const amrex::Real hy = mfy * dyinv * d.dy_drop;
    return std::sqrt(hx*hx + hy*hy);
}

Referenced by FaceBandSlopes().

Here is the call graph for this function:

Here is the caller graph for this function:

FaceBandSlopes() [1/2]

amrex::Array<amrex::Real,4> erf_terrain_file::FaceBandSlopes	(	const amrex::MultiFab &	z_phys_nd,
		const amrex::MultiFab &	mf_mx,
		const amrex::MultiFab &	mf_my,
		const amrex::Geometry &	geom,
		int	band
	)

Largest CellTerrainSlope, on one level, over the surface cells within band cells of each lateral domain face, in the order xlo, xhi, ylo, yhi; -1 for a face whose band holds no cell of this level. Reduced over all ranks.

Parameters

[in]	z_phys_nd	nodal heights of the level (surface at the domain's lowest k)
[in]	mf_mx	map factor in x at cell centres (ERF's mapfac m_x, on ba2d)
[in]	mf_my	map factor in y at cell centres (m_y)
[in]	geom	geometry of the level
[in]	band	width of the band in cells (>= 1)

FaceBandSlopes() [2/2]

Array<Real,4> erf_terrain_file::FaceBandSlopes	(	const MultiFab &	z_phys_nd,
		const MultiFab &	mf_mx,
		const MultiFab &	mf_my,
		const Geometry &	geom,
		int	band
	)

{
    AMREX_ALWAYS_ASSERT(band >= 1);
    const Box& dom = geom.Domain();
    const int klo = dom.smallEnd(2);
    const Real dxinv = geom.InvCellSize(0);
    const Real dyinv = geom.InvCellSize(1);
 
    // The bands of cells next to each lateral face, at the surface only
    Array<Box,4> bands;
    for (int f = 0; f < 4; ++f) {
        Box b(dom);
        b.setRange(2, klo);
        const int dir = f / 2;
        if (f % 2 == 0) {
            b.setBig(dir, amrex::min(dom.smallEnd(dir) + band - 1, dom.bigEnd(dir)));
        } else {
            b.setSmall(dir, amrex::max(dom.bigEnd(dir) - band + 1, dom.smallEnd(dir)));
        }
        bands[f] = b;
    }
 
    Array<Real,4> smax{Real(-1.), Real(-1.), Real(-1.), Real(-1.)};
    for (int f = 0; f < 4; ++f) {
        ReduceOps<ReduceOpMax> reduce_op;
        ReduceData<Real> reduce_data(reduce_op);
        using ReduceTuple = typename decltype(reduce_data)::Type;
        bool any = false;
        for (MFIter mfi(z_phys_nd); mfi.isValid(); ++mfi) {
            // z_phys_nd is nodal: its valid cells are the enclosed cells of the nodal box
            const Box bx = enclosedCells(mfi.validbox()) & bands[f];
            if (bx.isEmpty()) { continue; }
            any = true;
            const Array4<const Real> z_nd = z_phys_nd.const_array(mfi);
            // the map factors live on the same boxes flattened to one k (ba2d)
            const Array4<const Real> mfx = mf_mx.const_array(mfi);
            const Array4<const Real> mfy = mf_my.const_array(mfi);
            reduce_op.eval(bx, reduce_data, [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple
            {
                return {CellTerrainSlope(i, j, k, z_nd, mfx(i,j,0), mfy(i,j,0), dxinv, dyinv)};
            });
        }
        if (any) { smax[f] = amrex::get<0>(reduce_data.value(reduce_op)); }
    }
    ParallelDescriptor::ReduceRealMax(smax.data(), 4);
    return smax;
}

Referenced by ERF::WarnSteepTerrainAtOpenFaces().

Here is the call graph for this function:

Here is the caller graph for this function:

OpenFaceSlopeWarning() [1/2]

std::string erf_terrain_file::OpenFaceSlopeWarning	(	int	lev,
		int	band,
		const amrex::Array< amrex::Real, 4 > &	slopes,
		const amrex::Array< bool, 4 > &	open,
		const amrex::Array< std::string, 4 > &	bc_names
	)

Warning text for the inflow and outflow faces of one level whose band slope (from FaceBandSlopes) exceeds open_face_slope_warn, or an empty string if there is none. The caller prints it only when no steep-terrain option is on at that level.

Parameters

[in]	lev	level, for the text
[in]	band	band width in cells, for the text
[in]	slopes	band slopes in the order xlo, xhi, ylo, yhi (-1 = no cells)
[in]	open	whether each face is an inflow or outflow face
[in]	bc_names	boundary type of each face, for the text

OpenFaceSlopeWarning() [2/2]

std::string erf_terrain_file::OpenFaceSlopeWarning	(	int	lev,
		int	band,
		const Array< Real, 4 > &	slopes,
		const Array< bool, 4 > &	open,
		const Array< std::string, 4 > &	bc_names
	)

{
    static const char* face_names[4] = {"xlo", "xhi", "ylo", "yhi"};
    std::ostringstream os;
    for (int f = 0; f < 4; ++f) {
        if (open[f] && slopes[f] > open_face_slope_warn) {
            os << "WARNING: at level " << lev << " the terrain slope |grad h| reaches " << slopes[f]
               << " within " << band << " cells of the " << face_names[f] << " face ("
               << bc_names[f] << "), above " << open_face_slope_warn << "\n";
        }
    }
    if (os.str().empty()) { return ""; }
    os << "    Flattening the terrain within " << band << " cells of inflow and outflow faces avoids"
       << " this.  With an eddy-diffusivity closure, erf.smag2d_slope_limiter or"
       << " erf.smag2d_kh_cap (Smagorinsky2D), erf.implicit_terrain_metric (needs"
       << " erf.vert_implicit_fac > 0) or erf.diffusive_dt_limit limit the terrain-metric"
       << " diffusion on steep cells; none of them is on here\n";
    return os.str();
}

Referenced by ERF::WarnSteepTerrainAtOpenFaces().

Here is the caller graph for this function:

ReadHeightTable() [1/2]

const HeightTable& erf_terrain_file::ReadHeightTable	(	const std::string &	fname,
		bool	is_usgs,
		const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &	prob_lo,
		const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &	prob_hi,
		bool	require_domain_coverage = true
	)

The height table in file fname, read, checked (TableProblem) and summarised (TableSummary) on the I/O rank the first time the file is asked for, then broadcast; later calls return the same table without reading the file again. Aborts, naming the file and the line, on a file that cannot be read as the layout says.

Parameters

[in]	fname	file name
[in]	is_usgs	USGS layout
[in]	prob_lo,prob_hi	domain extent, for the coverage check
[in]	require_domain_coverage	false for local building patches; all other checks remain

ReadHeightTable() [2/2]

const HeightTable& erf_terrain_file::ReadHeightTable	(	const std::string &	fname,
		bool	is_usgs,
		const GpuArray< Real, AMREX_SPACEDIM > &	prob_lo,
		const GpuArray< Real, AMREX_SPACEDIM > &	prob_hi,
		bool	require_domain_coverage
	)

{
    // Include the coverage policy so a cached building patch cannot bypass the terrain check
    // if the same file is subsequently requested as terrain.
    static std::map<std::string, HeightTable> tables;
    const std::string key = std::string(require_domain_coverage ? "covered:" : "partial:")
                          + (is_usgs ? "USGS:" : "ERF:") + fname;
    auto it = tables.find(key);
    if (it == tables.end()) {
        it = tables.emplace(key, read_and_broadcast(fname, is_usgs, prob_lo, prob_hi,
                                                    require_domain_coverage)).first;
    }
    return it->second;
}

Referenced by ProblemBase::read_custom_terrain().

Here is the caller graph for this function:

TableProblem() [1/2]

std::string erf_terrain_file::TableProblem	(	bool	is_usgs,
		int	nx,
		int	ny,
		const amrex::Vector< amrex::Real > &	x,
		const amrex::Vector< amrex::Real > &	y,
		const amrex::Vector< amrex::Real > &	z,
		const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &	prob_lo,
		const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &	prob_hi,
		bool	require_domain_coverage = true
	)

The first problem found in a height table, or an empty string if there is none. Problems: fewer points than the layout needs (ERF format: one in each direction, USGS: two; one value in a direction is not checked further), list sizes that do not match nx and ny, a value that is not finite, x or y values that do not increase strictly (in the USGS format: the x values of each row, and the y values from row to row), and a table that does not cover the domain: each end of the x and y lists may stop at most one table spacing inside prob_lo or prob_hi (the layout of a file with one height per cell starting at prob_lo), with an allowance of 1e-3 of that spacing for rounded coordinates. Beyond the table the reader only extends the heights from its edge: linearly before the first point, and (ERF format) by holding the last value after the last point; USGS extends linearly at both ends. Setting require_domain_coverage to false skips only the coverage check, for building files that describe a local patch. Coordinate ordering and all format checks still apply.

Parameters

[in]	is_usgs	USGS layout (x varies fastest, one x per point)
[in]	nx,ny	table size
[in]	x,y,z	table coordinates and heights
[in]	prob_lo,prob_hi	domain extent
[in]	require_domain_coverage	whether the table must cover the domain

TableProblem() [2/2]

std::string erf_terrain_file::TableProblem	(	bool	is_usgs,
		int	nx,
		int	ny,
		const Vector< Real > &	x,
		const Vector< Real > &	y,
		const Vector< Real > &	z,
		const GpuArray< Real, AMREX_SPACEDIM > &	prob_lo,
		const GpuArray< Real, AMREX_SPACEDIM > &	prob_hi,
		bool	require_domain_coverage
	)

{
    std::ostringstream os;
    const int nmin = is_usgs ? 2 : 1;
    if (nx < nmin || ny < nmin) {
        os << "nx = " << nx << " and ny = " << ny << " must both be at least " << nmin
           << (is_usgs ? " in the USGS format" : "");
        return os.str();
    }
    const Long npts = Long(nx) * Long(ny);
    const Long nxv  = is_usgs ? npts : Long(nx);
    if (Long(x.size()) != nxv || Long(y.size()) != Long(ny) || Long(z.size()) != npts) {
        os << "expected " << nxv << " x values, " << ny << " y values and " << npts
           << " heights, found " << x.size() << ", " << y.size() << " and " << z.size();
        return os.str();
    }
    // A file cannot hold nan or inf (the reader stops on them as "not a number"); this is a
    // guard for tables built otherwise
    auto finite = [] (const Vector<Real>& v, const char* name) -> std::string {
        for (Long n = 0; n < Long(v.size()); ++n) {
            if (!std::isfinite(v[n])) {
                std::ostringstream o;
                o << name << " value " << n << " is not finite (" << v[n] << ")";
                return o.str();
            }
        }
        return "";
    };
    for (const auto& s : {finite(x,"x"), finite(y,"y"), finite(z,"height")}) {
        if (!s.empty()) { return s; }
    }
 
    std::string s = list_problem("y", ny, [&] (int m) { return y[m]; }, prob_lo[1], prob_hi[1],
                                 require_domain_coverage);
    if (!s.empty()) { return s; }
    if (is_usgs) {
        for (int j = 0; j < ny; ++j) {
            s = list_problem("x", nx, [&] (int m) { return x[j*nx + m]; }, prob_lo[0], prob_hi[0],
                             require_domain_coverage);
            if (!s.empty()) {
                return s + " (row " + std::to_string(j) + " of the x-fastest USGS layout)";
            }
        }
    } else {
        s = list_problem("x", nx, [&] (int m) { return x[m]; }, prob_lo[0], prob_hi[0],
                         require_domain_coverage);
        if (!s.empty()) { return s; }
    }
    return "";
}

Here is the call graph for this function:

TableSummary() [1/2]

std::string erf_terrain_file::TableSummary	(	const std::string &	fname,
		bool	is_usgs,
		int	nx,
		int	ny,
		const amrex::Vector< amrex::Real > &	x,
		const amrex::Vector< amrex::Real > &	y,
		const amrex::Vector< amrex::Real > &	z
	)

Text describing a valid height table: its size and extent, the order its heights are read in, where the lowest and highest heights sit, and the largest slopes |dh/dx| and |dh/dy| between neighbouring table points with where they sit. Meant to be printed once per file, so that a transposed or mis-ordered file can be spotted at start-up.

Parameters

[in]	fname	file name, for the first line
[in]	is_usgs	USGS layout
[in]	nx,ny	table size
[in]	x,y,z	table coordinates and heights (must pass TableProblem)

TableSummary() [2/2]

std::string erf_terrain_file::TableSummary	(	const std::string &	fname,
		bool	is_usgs,
		int	nx,
		int	ny,
		const Vector< Real > &	x,
		const Vector< Real > &	y,
		const Vector< Real > &	z
	)

{
    int imin = 0, jmin = 0, imax = 0, jmax = 0;
    Real zmin = z_at(is_usgs, nx, ny, z, 0, 0), zmax = zmin;
    Real sx = Real(-1.), sy = Real(-1.);
    Real sx_x = 0, sx_y = 0, sy_x = 0, sy_y = 0;
    for (int j = 0; j < ny; ++j) {
        for (int i = 0; i < nx; ++i) {
            const Real h = z_at(is_usgs, nx, ny, z, i, j);
            if (h < zmin) { zmin = h; imin = i; jmin = j; }
            if (h > zmax) { zmax = h; imax = i; jmax = j; }
            if (i+1 < nx) {
                const Real x0 = x_at(is_usgs, nx, x, i, j), x1 = x_at(is_usgs, nx, x, i+1, j);
                const Real s = std::abs(z_at(is_usgs, nx, ny, z, i+1, j) - h) / (x1 - x0);
                if (s > sx) { sx = s; sx_x = Real(0.5)*(x0 + x1); sx_y = y[j]; }
            }
            if (j+1 < ny) {
                const Real s = std::abs(z_at(is_usgs, nx, ny, z, i, j+1) - h) / (y[j+1] - y[j]);
                if (s > sy) { sy = s; sy_x = x_at(is_usgs, nx, x, i, j); sy_y = Real(0.5)*(y[j] + y[j+1]); }
            }
        }
    }
 
    std::ostringstream os;
    os << "Height file " << fname << ": " << nx << " x " << ny << " heights over x "
       << x_at(is_usgs, nx, x, 0, 0) << " .. " << x_at(is_usgs, nx, x, nx-1, 0)
       << ", y " << y[0] << " .. " << y[ny-1] << "\n";
    if (is_usgs) {
        os << "    read in USGS order: x varies fastest, height n is at (x[n], y[n / nx])\n";
    } else {
        os << "    read in ERF order: y varies fastest, height n is at (x[n / ny], y[n % ny])\n";
    }
    os << "    lowest  " << zmin << " m at (x, y) = (" << x_at(is_usgs, nx, x, imin, jmin)
       << ", " << y[jmin] << ")\n";
    os << "    highest " << zmax << " m at (x, y) = (" << x_at(is_usgs, nx, x, imax, jmax)
       << ", " << y[jmax] << ")\n";
    if (sx < Real(0.)) {
        os << "    one x value: constant in x";
    } else {
        os << "    largest |dh/dx| " << sx << " at (" << sx_x << ", " << sx_y << ")";
    }
    if (sy < Real(0.)) {
        os << ", one y value: constant in y\n";
    } else {
        os << ", largest |dh/dy| " << sy << " at (" << sy_x << ", " << sy_y << ")\n";
    }
    return os.str();
}

Variable Documentation

open_face_slope_warn

constexpr amrex::Real erf_terrain_file::open_face_slope_warn = amrex::Real(0.3)

constexpr

Terrain slope |grad h| above which a band next to an inflow or outflow face is reported.

Referenced by OpenFaceSlopeWarning(), and ERF::WarnSteepTerrainAtOpenFaces().