ERF_TerrainMetrics.cpp File Reference

#include <ERF_TerrainMetrics.H>
#include <ERF_Utils.H>
#include <AMReX_ParmParse.H>
#include <ERF_Constants.H>
#include <ERF_Interpolation_1D.H>
#include <algorithm>
#include <array>
#include <cmath>
#include <limits>
#include <map>
#include <numeric>
#include <utility>
#include <vector>

Include dependency graph for ERF_TerrainMetrics.cpp:

Functions
void	validate_flat_terrain (int lev, const MultiFab &z_phys_nd, const Vector< Real > &z_levels_h)
void	validate_flat_top (int lev, const Geometry &geom, const MultiFab &z_phys_nd)
void	init_default_zphys (int, const Geometry &geom, MultiFab &z_phys_nd, MultiFab &z_phys_cc, Real z_offset)
BoxArray	join_boxes_stacked_in_z (const BoxArray &ba)
void	make_terrain_fitted_coords (int lev, const Geometry &geom, MultiFab &z_phys_nd, Vector< Real > const &z_levels_h, GpuArray< ERF_BC, AMREX_SPACEDIM *2 > &phys_bc_type, FineTerrain fine_terrain, MultiFab const *z_phys_interp)
FineTerrain	which_fine_terrain ()
void	init_fine_terrain_grid (int lev, Geometry const &geom, MultiFab &z_phys_nd, Vector< Real > const &z_levels_h, FineTerrain fine_terrain, MultiFab const *z_phys_interp)
void	init_which_terrain_grid (int lev, Geometry const &geom, MultiFab &z_phys_nd, Vector< Real > const &z_levels_h, FineTerrain fine_terrain, MultiFab const *z_phys_interp)
void	make_J (const Geometry &geom, MultiFab &z_phys_nd, MultiFab &detJ_cc)
void	make_areas (const Geometry &geom, MultiFab &z_phys_nd, MultiFab &ax, MultiFab &ay, MultiFab &az)
void	make_zcc (const Geometry &geom, MultiFab &z_phys_nd, MultiFab &z_phys_cc)
Real	get_dzmin_terrain (MultiFab &z_phys_nd)

Function Documentation

get_dzmin_terrain()

Real get_dzmin_terrain

(

MultiFab &

z_phys_nd

)

Computation min dz at cell-center

{
    auto const& ma_z_nd_arr = z_phys_nd.const_arrays();
    GpuTuple<Real> min = ParReduce(TypeList<ReduceOpMin>{},
                                   TypeList<Real>{},
                                   z_phys_nd, IntVect(0),
                                   [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept
                                   -> GpuTuple<Real>
                                   {
                                       amrex::Real dz = Compute_Z_AtWFace(i,j,k+1,ma_z_nd_arr[box_no]) -
                                                        Compute_Z_AtWFace(i,j,k  ,ma_z_nd_arr[box_no]);
                                       return { dz };
                                   });
    Real r = (get<0>(min) + std::numeric_limits<amrex::Real>::epsilon());
    ParallelDescriptor::ReduceRealMin(r);
    return r;
}

Referenced by ERF::init_zphys(), ERF::ReadCheckpointFile(), and ERF::remake_zphys().

Here is the call graph for this function:

Here is the caller graph for this function:

init_default_zphys()

void init_default_zphys	(	int	,
		const Geometry &	geom,
		MultiFab &	z_phys_nd,
		MultiFab &	z_phys_cc,
		Real	z_offset
	)

Define a default z_phys so we have it even if a completely regular grid This will be over-written if we use z_levels, or grid stretching, or terrain-fitted grids

{
    const auto& dx = geom.CellSize();
    Real dz = dx[2];
    Real z_bot = geom.ProbLo(2) - z_offset;
 
    for (MFIter mfi(z_phys_nd,true); mfi.isValid(); ++mfi)
    {
        const Box& bx = mfi.growntilebox();
        const Array4< Real> z_nd_arr = z_phys_nd.array(mfi);
        ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
        {
            z_nd_arr(i,j,k) = z_bot + k * dz;
        });
    }
 
    for (MFIter mfi(z_phys_cc,true); mfi.isValid(); ++mfi)
    {
        const Box& bx = mfi.growntilebox();
        const Array4< Real> z_cc_arr = z_phys_cc.array(mfi);
        ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
        {
            z_cc_arr(i,j,k) = z_bot + (k + myhalf) * dz;
        });
    }
}

Referenced by ERF::init_stuff().

Here is the call graph for this function:

Here is the caller graph for this function:

init_fine_terrain_grid()

void init_fine_terrain_grid	(	int	lev,
		Geometry const &	geom,
		MultiFab &	z_phys_nd,
		Vector< Real > const &	z_levels_h,
		FineTerrain	fine_terrain,
		MultiFab const *	z_phys_interp
	)

Finish the terrain-fitted mesh on a fine level built from the coarse level.

On entry z_phys_nd holds the mesh interpolated from the coarse level in its valid region, with – for Transform – the fine-resolution terrain already written into the k=0 slab. The nodes outside the domain are still bogus_large_value, because neither the interpolation nor make_terrain_fitted_coords()'s setDomainBndry() call gave them a usable value, so we fill them here the same way the BTF branch does.

{
    AMREX_ALWAYS_ASSERT(lev > 0);
    AMREX_ALWAYS_ASSERT(fine_terrain != FineTerrain::None);
 
    const Box& domain = geom.Domain();
    int domlo_x = domain.smallEnd(0); int domhi_x = domain.bigEnd(0) + 1;
    int domlo_y = domain.smallEnd(1); int domhi_y = domain.bigEnd(1) + 1;
    int domlo_z = domain.smallEnd(2);
 
    IntVect ngrowVect = z_phys_nd.nGrowVect();
 
    if (fine_terrain == FineTerrain::Transform) {
        AMREX_ALWAYS_ASSERT(z_phys_interp != nullptr);
 
        const BoxArray& ba = z_phys_nd.boxArray();
 
        //
        // The correction is anchored at the surface, so a box that does not reach the
        // surface would keep the uncorrected interpolated mesh and would not join the
        // box below it.  Say so rather than silently building a tangled mesh.
        //
        for (int i = 0; i < ba.size(); i++) {
            if (ba[i].smallEnd(2) != domlo_z) {
                Abort("erf.amr_terrain_refinement = transform requires that every box on a fine level reach the surface: increase amr.max_grid_size_z or use erf.amr_terrain_refinement = interpolate");
            }
        }
 
        //
        // Decay over the height of the fine grids as a whole and never over the height
        // of an individual box, so that the mesh does not depend on the decomposition.
        //
        int khi_lev = ba.minimalBox().bigEnd(2);
        Real z_lev_top = z_levels_h[khi_lev];
        Real z_lev_bot = z_levels_h[domlo_z];
        AMREX_ALWAYS_ASSERT(z_lev_top > z_lev_bot);
 
        Gpu::DeviceVector<Real> z_levels_d(z_levels_h.size());
        Gpu::copy(Gpu::hostToDevice, z_levels_h.begin(), z_levels_h.end(), z_levels_d.begin());
        auto const& z_lev = z_levels_d.data();
 
        //
        // Correct the ghost nodes along with the valid ones.  On a fine level whose
        // grids do not cover the domain, the lateral ghost nodes have no valid owner
        // anywhere, so the FillBoundary below cannot bring them up to date: left
        // uncorrected they would keep the interpolated mesh above the surface while
        // their k=0 node already holds the fine terrain, and the first layer of such a
        // column inverts wherever the fine terrain stands more than a layer thickness
        // above the interpolated one.  Both fields the correction reads are defined on
        // those nodes -- but only inside the domain, so stop there and leave the nodes
        // outside it to the clamp below, which is what gives them a mesh at all.
        //
        const Box domain_nd = convert(domain, IntVect(1,1,1));
 
        for (MFIter mfi(z_phys_nd, TilingIfNotGPU()); mfi.isValid(); ++mfi)
        {
            const Box bx = mfi.growntilebox(ngrowVect) & domain_nd;
 
            Array4<Real      > const& z_arr    = z_phys_nd.array(mfi);
            Array4<Real const> const& z_interp = z_phys_interp->const_array(mfi);
 
            ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
            {
                if (k > domlo_z && k <= khi_lev) {
                    // How much taller the fine terrain is than the interpolated surface
                    Real delta_terrain = z_arr(i,j,domlo_z) - z_interp(i,j,domlo_z);
 
                    Real decay = (z_lev_top - z_lev[k]) / (z_lev_top - z_lev_bot);
                    decay = amrex::max(Real(0.0), amrex::min(Real(1.0), decay));
 
                    z_arr(i,j,k) = z_interp(i,j,k) + decay * delta_terrain;
                }
            });
        } // mfi
 
        // Give every ghost node that does have a valid owner that owner's value, rather
        // than the one this box computed for itself, before the fill below reads them
        z_phys_nd.FillBoundary(geom.periodicity());
    }
 
    //
    // Extend the mesh to the nodes outside the lateral domain boundary by copying the
    // nearest column inside it, which is the value the BTF branch gives those nodes.
    // Like that branch we clamp in a periodic direction too, even though FillBoundary
    // below will replace what we write there with the periodic image: that way every
    // node we read here is one inside the domain, whose value we know is already set.
    //
    int imin = domlo_x; int imax = domhi_x;
    int jmin = domlo_y; int jmax = domhi_y;
 
    for (MFIter mfi(z_phys_nd, TilingIfNotGPU()); mfi.isValid(); ++mfi)
    {
        Box gbx = mfi.growntilebox(ngrowVect);
        Array4<Real> const& z_arr = z_phys_nd.array(mfi);
 
        // Note that we only write nodes outside the domain and only read nodes inside
        // it, so no node is both read and written here.
        ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
        {
            int ii = amrex::max(amrex::min(i,imax),imin);
            int jj = amrex::max(amrex::min(j,jmax),jmin);
            if (ii != i || jj != j) {
                z_arr(i,j,k) = z_arr(ii,jj,k);
            }
        });
    } // mfi
 
    z_phys_nd.FillBoundary(geom.periodicity());
 
    for (MFIter mfi(z_phys_nd, TilingIfNotGPU()); mfi.isValid(); ++mfi)
    {
        const Box& bx = mfi.validbox();
        Box gbx = mfi.growntilebox(ngrowVect);
 
        if (bx.smallEnd(2) == domlo_z)
        {
            Array4<Real> const& z_arr = z_phys_nd.array(mfi);
 
            // Fill lateral boundaries below the bottom surface
            ParallelFor(makeSlab(gbx,2,domlo_z), [=] AMREX_GPU_DEVICE (int i, int j, int)
            {
                z_arr(i,j,domlo_z-1) = two*z_arr(i,j,domlo_z) - z_arr(i,j,domlo_z+1);
            });
        }
    } // mfi
} // init_fine_terrain_grid

Referenced by init_which_terrain_grid().

Here is the call graph for this function:

Here is the caller graph for this function:

init_which_terrain_grid()

void init_which_terrain_grid	(	int	lev,
		Geometry const &	geom,
		MultiFab &	z_phys_nd,
		Vector< Real > const &	z_levels_h,
		FineTerrain	fine_terrain,
		MultiFab const *	z_phys_interp
	)

Initialize the terrain grid using the selected model (BTF, STF, or Sullivan TF).

Parameters

[in]	lev	Level in the AMR hierarchy.
[in]	geom	Geometry defining the domain.
[out]	z_phys_nd	MultiFab to be filled with nodal physical heights.
[in]	z_levels_h	Reference z-levels for grid stretching.

{
    // User-selected method from inputs file (BTF default)
    ParmParse pp("erf");
    int terrain_smoothing = 0;
    pp.queryAdd("terrain_smoothing", terrain_smoothing);
 
    //
    // The STF and Sullivan TF transformations need the full column from the surface to
    // the domain top, so on a fine level we build the mesh from the one interpolated
    // from the coarse level instead.  Only a caller that has done that interpolation
    // (and told us so) may take this path; everyone else still aborts.
    //
    if (lev > 0 && terrain_smoothing != 0) {
        if (fine_terrain == FineTerrain::None) {
            Abort("terrain_smoothing = " + std::to_string(terrain_smoothing) +
                  " on level > 0 requires setting erf.amr_terrain_refinement (transform or interpolate)");
        }
        init_fine_terrain_grid(lev, geom, z_phys_nd, z_levels_h, fine_terrain, z_phys_interp);
        return;
    }
 
    // Number of ghost cells
    int ngrow = z_phys_nd.nGrow();
    IntVect ngrowVect = z_phys_nd.nGrowVect();
 
    const Box& domain = geom.Domain();
    int domlo_x = domain.smallEnd(0); int domhi_x = domain.bigEnd(0) + 1;
    int domlo_y = domain.smallEnd(1); int domhi_y = domain.bigEnd(1) + 1;
    int domlo_z = domain.smallEnd(2); int domhi_z = domain.bigEnd(2) + 1;
 
    int imin = domlo_x; // if (geom.isPeriodic(0)) imin -= z_phys_nd.nGrowVect()[0];
    int jmin = domlo_y; // if (geom.isPeriodic(1)) jmin -= z_phys_nd.nGrowVect()[1];
 
    int imax = domhi_x; // if (geom.isPeriodic(0)) imax += z_phys_nd.nGrowVect()[0];
    int jmax = domhi_y; // if (geom.isPeriodic(1)) jmax += z_phys_nd.nGrowVect()[1];
 
    int nz = static_cast<int>(z_levels_h.size());
    Real z_top = z_levels_h[nz-1];
 
    // The STF and Sullivan transforms measure heights and terrain from the bottom of the
    // domain, which is z_levels_h[0] = geometry.prob_lo[2]; the BTF one does so already
    Real z_bot = z_levels_h[0];
    Real depth = z_top - z_bot;
 
    Gpu::DeviceVector<Real> z_levels_d;
    z_levels_d.resize(nz);
    Gpu::copy(Gpu::hostToDevice, z_levels_h.begin(), z_levels_h.end(), z_levels_d.begin());
 
    switch(terrain_smoothing) {
    case 0: // BTF Method
    {
        for ( MFIter mfi(z_phys_nd, TilingIfNotGPU()); mfi.isValid(); ++mfi )
        {
            // Note that this box is nodal because it is based on z_phys_nd
            const Box&  bx = mfi.validbox();
 
            int k0 = bx.smallEnd()[2];
 
            // Grown box with corrected ghost cells at top
            Box gbx = mfi.growntilebox(ngrowVect);
 
            if (bx.smallEnd(2) == domlo_z) {
                gbx.setSmall(2,domlo_z);
            } else {
                gbx.growLo(2,-1);
            }
 
            // Note that we don't overwrite the values at the high end of the box
            // regardless of whether the box reaches the top of the domain or not.
            // In the case of lev > 0, this ensures that the fine nodes at the top
            // of a fine box are those that are interpolated from the coarse grid
            if (bx.bigEnd(2) == domhi_z) {
                gbx.setBig(2,domhi_z);
            } else {
                gbx.growHi(2,-1);
                if (gbx.bigEnd(2) > domhi_z) {
                    gbx.setBig(2,domhi_z);
                }
            }
 
            Array4<Real> const& z_arr = z_phys_nd.array(mfi);
            auto const&         z_lev = z_levels_d.data();
 
            //
            // Vertical grid stretching using BTF model from p2163 of Klemp2011
            // z_levels are only defined from k = dom_lo to dom_hi (nodal)
            //
            ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
            {
                int ii = amrex::max(amrex::min(i,imax),imin);
                int jj = amrex::max(amrex::min(j,jmax),jmin);
 
                //
                // Start with flat z_lev set either with uniform cell size or specified z_levels
                // If k0 = 0 then z_arr at k0 has already been filled from the terrain data
                // If k0 > 0 then z_arr at k0 has already been filled from interpolation
                //
                // A box that does not start at the bottom (e.g. a refined box starting at
                // k0 = 1) has ghost nodes below the domain, where z_lev is not defined;
                // extend the levels there by reflection about the bottom
                Real z         = (k >= domlo_z) ? z_lev[k] : two*z_lev[domlo_z] - z_lev[2*domlo_z-k];
                Real z_sfc     = z_arr(ii,jj,k0);
                Real z_lev_sfc = z_lev[k0];
 
                z_arr(i,j,k) = ( (z_sfc - z_lev_sfc) * z_top +
                                 (z_top - z_sfc    ) * z     ) / (z_top - z_lev_sfc);
            });
        } // mfi
 
        z_phys_nd.FillBoundary(geom.periodicity());
 
        for ( MFIter mfi(z_phys_nd, TilingIfNotGPU()); mfi.isValid(); ++mfi )
        {
            // Note that this box is nodal because it is based on z_phys_nd
            const Box&  bx = mfi.validbox();
            Box gbx = mfi.growntilebox(ngrowVect);
 
            int k0 = bx.smallEnd()[2];
 
            if (k0 == 0)
            {
                Array4<Real> const& z_arr = z_phys_nd.array(mfi);
 
                // Fill lateral boundaries below the bottom surface
                ParallelFor(makeSlab(gbx,2,0), [=] AMREX_GPU_DEVICE (int i, int j, int)
                {
                    z_arr(i,j,-1) = two*z_arr(i,j,0) - z_arr(i,j,1);
                });
            }
        } // mfi
        break;
    } // case 0
 
    case 1: // STF Method
    {
        // Get MultiFab spanning domain with 1 level of ghost cells
        MultiFab h_mf(    z_phys_nd.boxArray(), z_phys_nd.DistributionMap(), 1, ngrow+1);
        MultiFab h_mf_old(z_phys_nd.boxArray(), z_phys_nd.DistributionMap(), 1, ngrow+1);
 
        // Save max height for smoothing
        Real h_m;
 
        // Create 2D MF without allocation
        MultiFab mf2d;
        {
            BoxList bl2d = h_mf.boxArray().boxList();
            for (auto& b : bl2d) { b.setRange(2,b.smallEnd(2)); }
            BoxArray ba2d(std::move(bl2d));
            mf2d = MultiFab(ba2d, h_mf.DistributionMap(), 1, ngrow, MFInfo().SetAlloc(false));
        }
 
        // Get MultiArray4s from the multifabs
        MultiArray4<Real> const& ma_h_s     = h_mf.arrays();
        MultiArray4<Real> const& ma_h_s_old = h_mf_old.arrays();
        MultiArray4<Real> const& ma_z_phys  = z_phys_nd.arrays();
 
        // Bottom boundary
        int k0 = domlo_z;
 
        // Get max value
        h_m = ParReduce(TypeList<ReduceOpMax>{}, TypeList<Real>{}, mf2d, IntVect(ngrow,ngrow,0),
                    [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int) noexcept
                        -> GpuTuple<Real>
                {
                  // Get Array4s
                  const auto & h     = ma_h_s[box_no];
                  const auto & z_arr = ma_z_phys[box_no];
 
                  int ii = amrex::max(amrex::min(i,imax),imin);
                  int jj = amrex::max(amrex::min(j,jmax),jmin);
 
                  // Fill the lateral boundaries
                  z_arr(i,j,k0) = z_arr(ii,jj,k0);
 
                  // Populate h with the terrain height above the bottom of the domain
                  h(i,j,k0) = z_arr(i,j,k0) - z_bot;
 
                  // Return height for max
                  return { h(i,j,k0) };
                });
        amrex::ParallelDescriptor::ReduceRealMax(h_m);
 
        if (h_m < std::numeric_limits<Real>::epsilon()) h_m = Real(1e-16);
 
        // Fill ghost cells (neglects domain boundary if not periodic)
        h_mf.FillBoundary(geom.periodicity());
 
        // Make h_mf copy for old values
        MultiFab::Copy(h_mf_old, h_mf,0,0,1,h_mf_old.nGrow());
 
        // Minimum allowed fractional grid spacing
        Real gamma_m = myhalf;
        pp.queryAdd("terrain_gamma_m", gamma_m);
        Real z_H     = Real(2.44)*h_m/(1-gamma_m); // Klemp2011 Eqn. 11
 
        // Levels at which the minimum-thickness criterion forced us to stop smoothing
        // early, reported once below: a terrain that needs this at many levels is at the
        // edge of what STF can represent on this vertical grid.
        int n_capped = 0;
 
        // Populate h_mf at k>0 with h_s, solving in ordered 2D slices
        for (int k = domlo_z+1; k <= domhi_z; k++) // skip terrain level
        {
            auto const& z_lev_h = z_levels_h.data();
 
            Real zz       = z_lev_h[k]   - z_bot;
            Real zz_minus = z_lev_h[k-1] - z_bot;
 
            // Hybrid attenuation profile, Klemp2011 Eqn. 9
            Real A;
            Real foo = std::cos((PI/2)*(zz/z_H));
            if(zz < z_H) { A = foo*foo*foo*foo*foo*foo; } // A controls rate of return to atm
            else         { A = 0; }
            Real foo_minus = std::cos((PI/2)*(zz_minus/z_H));
            Real A_minus;
            if(zz_minus < z_H) { A_minus = foo_minus*foo_minus*foo_minus*foo_minus*foo_minus*foo_minus; } // A controls rate of return to atm
            else               { A_minus = 0; }
 
            unsigned maxIter = 50; // M_k in paper
            unsigned iter    = 0;
            Real threshold   = gamma_m;
            Real diff        = bogus_large_value;
            while (iter < maxIter)
            {
 
                diff = ParReduce(TypeList<ReduceOpMin>{}, TypeList<Real>{}, mf2d, IntVect(ngrow,ngrow,0),
                    [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int) noexcept
                        -> GpuTuple<Real>
                {
                    const auto & h_s     = ma_h_s[box_no];
                    const auto & h_s_old = ma_h_s_old[box_no];
 
                    Real beta_k = Real(0.2)*std::min(zz/(2*h_m),Real(1.0)); //smoothing coefficient (Eqn. 8)
 
                    // Clip indices for ghost-cells
                    int ii = amrex::min(amrex::max(i,domlo_x),domhi_x);
                    int jj = amrex::min(amrex::max(j,domlo_y),domhi_y);
 
                    if (iter == 0) {
                        h_s(i,j,k) = h_s_old(i,j,k-1) + beta_k*(h_s_old(ii+1,jj  ,k-1)
                                                              + h_s_old(ii-1,jj  ,k-1)
                                                              + h_s_old(ii  ,jj+1,k-1)
                                                              + h_s_old(ii  ,jj-1,k-1) - 4*h_s_old(ii,jj,k-1));
                    }
                    else {
                        h_s(i,j,k) = h_s_old(i,j,k  ) + beta_k*(h_s_old(ii+1,jj  ,k  )
                                                              + h_s_old(ii-1,jj  ,k  )
                                                              + h_s_old(ii  ,jj+1,k  )
                                                              + h_s_old(ii  ,jj-1,k  ) - 4*h_s_old(ii,jj,k  ));
                    }
 
                    // Minimum vertical grid spacing condition (Klemp2011 Eqn. 7)
                    return { (zz + A * h_s(i,j,k) - (zz_minus + A_minus * h_s(i,j,k-1))) / (zz - zz_minus) };
 
                }); //ParReduce
 
                ParallelDescriptor::ReduceRealMin(diff);
 
                //
                // Each pass smooths h_s at THIS level while the level below is already
                // fixed, so over a peak it always thins the layer: the passes have to
                // stop before Eqn. 7 is violated, not after.  Measuring only once a pass
                // has been committed is not enough -- nothing then stops the pass that
                // breaks the criterion from being the one we keep, which is how a column
                // over steep terrain ended up at 0.21 of its nominal layer thickness
                // against gamma_m = 0.5.  Reject that pass and restore what preceded it.
                //
                // At iter 0 there is no earlier pass at this level to fall back to, so
                // carry the level below up unsmoothed instead.  That leaves only the
                // A-decay term, whose worst case is exactly what z_H = 2.44 h_m/(1-gamma_m)
                // was chosen to bound: max|dA/dz| = 2.439/z_H, so h_m max|dA/dz| dz equals
                // (1-gamma_m) dz and the thickness ratio lands on gamma_m.
                //
                if (diff < threshold) {
                    ++n_capped;
                    if (iter == 0) {
                        for ( MFIter mfi(h_mf, TilingIfNotGPU()); mfi.isValid(); ++mfi )
                        {
                            Box xybx = mfi.growntilebox(ngrow);
                            xybx.setRange(2,0);
                            Array4<Real> const& h_s = h_mf.array(mfi);
                            ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int) noexcept
                            {
                                h_s(i,j,k) = h_s(i,j,k-1);
                            });
                        }
                        MultiFab::Copy(h_mf_old, h_mf,0,0,1,h_mf_old.nGrow());
                        h_mf_old.FillBoundary(geom.periodicity());
                    } else {
                        // h_mf_old still holds the last pass that satisfied the criterion
                        MultiFab::Copy(h_mf, h_mf_old,0,0,1,h_mf.nGrow());
                    }
                    break;
                }
 
                MultiFab::Copy(h_mf_old, h_mf,0,0,1,h_mf_old.nGrow());
 
                iter++;
 
                // Fill ghost points
                h_mf_old.FillBoundary(geom.periodicity());
 
            } // while
 
            auto const& z_lev_d = z_levels_d.data();
 
            // Populate z_phys_nd by solving z_arr(i,j,k) = z + A*h_s(i,j,k)
            for ( MFIter mfi(z_phys_nd, TilingIfNotGPU()); mfi.isValid(); ++mfi )
            {
                // Grown box with no z range
                Box xybx = mfi.growntilebox(ngrow);
                xybx.setRange(2,0);
 
                Array4<Real> const& h_s   = h_mf_old.array(mfi);
                Array4<Real> const& z_arr = z_phys_nd.array(mfi);
 
                ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int) {
 
                    // Location of nodes
                    Real z = z_lev_d[k];
 
                    // STF model from p2164 of Klemp2011 (Eqn. 4)
                    z_arr(i,j,k) = z + A*h_s(i,j,k);
 
                    // Fill below the bottom surface
                    if (k == 1) {
                        z_arr(i,j,k0-1) = two*z_arr(i,j,k0) - z_arr(i,j,k);
                    }
                });
            } // mfi
        } // k
 
            if (n_capped > 0) {
                amrex::Print() << "NOTE: the smoothed terrain-following grid at level " << lev
                               << " hit its minimum layer-thickness floor (erf.terrain_gamma_m = "
                               << gamma_m << ") at " << n_capped << " of "
                               << (domhi_z - domlo_z) << " levels, so the terrain smoothing was "
                               << "capped there.  Raising erf.terrain_gamma_m, coarsening the "
                               << "vertical grid or using erf.terrain_smoothing = 0 would give "
                               << "this terrain more room.\n";
            }
 
            Gpu::streamSynchronize();
 
            break;
        } // case 1
 
        case 2: // Sullivan TF Method
        {
            int k0    = 0;
 
            for ( MFIter mfi(z_phys_nd, TilingIfNotGPU()); mfi.isValid(); ++mfi )
            {
                // Grown box with corrected ghost cells at top
                Box gbx = mfi.growntilebox(ngrow);
                gbx.setRange(2,domlo_z,domhi_z+1);
 
                Array4<Real> const& z_arr = z_phys_nd.array(mfi);
                auto const&         z_lev = z_levels_d.data();
 
                ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
                {
                    // Vertical grid stretching
                    Real z =  z_lev[k];
 
                    int ii = amrex::max(amrex::min(i,imax),imin);
                    int jj = amrex::max(amrex::min(j,jmax),jmin);
 
                    // Fill levels using model from Sullivan et. al. 2014
                    int omega = 3; //Used to adjust how rapidly grid lines level out. omega=1 is BTF!
                    z_arr(i,j,k) = z + static_cast<Real>(std::pow((one - (z - z_bot)/depth),omega) * (z_arr(ii,jj,k0) - z_bot));
 
                    // Fill lateral boundaries and below the bottom surface
                    if (k == k0) {
                        z_arr(i,j,k0  ) = z_arr(ii,jj,k0);
                    }
                    if (k == 1) {
                        z_arr(i,j,k0-1) = two*z_arr(ii,jj,k0) - z_arr(i,j,k);
                    }
                });
            } // mfi
            break;
        } // case 2
 
        case 3: // Debugging Test Method -- applies Sullivan TF starting at k = 1 so that domain does not change size
        {
            int k0    = 0;
 
            for ( MFIter mfi(z_phys_nd, TilingIfNotGPU()); mfi.isValid(); ++mfi )
            {
                // Grown box with corrected ghost cells at top
                Box gbx = mfi.growntilebox(ngrow);
                gbx.setRange(2,domlo_z,domhi_z+1);
 
                Array4<Real> const& z_arr = z_phys_nd.array(mfi);
                auto const&         z_lev = z_levels_d.data();
 
                ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
                {
                    // Vertical grid stretching
                    Real z =  z_lev[k];
 
                    int ii = amrex::max(amrex::min(i,imax),imin);
                    int jj = amrex::max(amrex::min(j,jmax),jmin);
 
                    // Fill values outside the lateral boundaries and below the bottom surface (necessary if init_type = "WRFInput")
                    if (k == k0+1) {
                        z_arr(i,j,k) = z + (z_arr(ii,jj,k0) - z_bot);
                    } else {
                        // Fill levels using model from Sullivan et. al. 2014
                        int omega = 3; //Used to adjust how rapidly grid lines level out. omega=1 is BTF!
                        z_arr(i,j,k) = z + static_cast<Real>(std::pow((one - (z - z_bot)/depth),omega) * (z_arr(ii,jj,k0) - z_bot));
                    }
                });
                gbx.setBig(2,0);
                ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
                {
                    z_arr(i,j,k  ) = z_bot;
                    z_arr(i,j,k-1) = two*z_bot - z_arr(i,j,k+1);
                });
            } // mfi
            break;
        } // case 3
    } //switch
}

Referenced by make_terrain_fitted_coords().

Here is the call graph for this function:

Here is the caller graph for this function:

join_boxes_stacked_in_z()

BoxArray join_boxes_stacked_in_z

(

const BoxArray &

ba

)

The cells of a (cell-centered) BoxArray as boxes that are never stacked on each other.

Wherever a box of ba sits on the top face of another box of ba, the run of cells in z that the two belong to becomes a single box. Boxes stacked with different lateral extents are first cut at each other's lateral edges; the joined pieces are then merged in x and y and chopped to no more than the largest lateral size of the boxes they came from. Returns ba itself when no two boxes of ba are stacked.

{
    AMREX_ALWAYS_ASSERT(ba.ixType().cellCentered());
 
    const int nboxes = static_cast<int>(ba.size());
 
    // Gather the boxes into stacks: link every box to the boxes that sit on its top face
    std::vector<int> stack_of(nboxes);
    std::iota(stack_of.begin(), stack_of.end(), 0);
    auto find_stack = [&stack_of] (int i) {
        while (stack_of[i] != i) {
            stack_of[i] = stack_of[stack_of[i]];
            i = stack_of[i];
        }
        return i;
    };
 
    bool any_stacked = false;
    for (int i = 0; i < nboxes; ++i) {
        Box above(ba[i]);
        above.setRange(2, ba[i].bigEnd(2)+1);
        for (const auto& isect : ba.intersections(above)) {
            stack_of[find_stack(isect.first)] = find_stack(i);
            any_stacked = true;
        }
    }
    if (!any_stacked) { return ba; }
 
    // std::map keeps the order, and so the joined boxes, the same on every rank
    std::map<int, std::vector<Box>> stacks;
    for (int i = 0; i < nboxes; ++i) {
        stacks[find_stack(i)].push_back(ba[i]);
    }
 
    BoxList bl_joined;
    for (const auto& stack : stacks)
    {
        const std::vector<Box>& boxes = stack.second;
        if (boxes.size() == 1) {
            bl_joined.push_back(boxes[0]);
            continue;
        }
 
        // Cut the stack at every lateral edge of its boxes
        std::array<std::vector<int>,2> cuts;
        IntVect max_size(1);
        for (const auto& b : boxes) {
            for (int idim = 0; idim < 2; ++idim) {
                cuts[idim].push_back(b.smallEnd(idim));
                cuts[idim].push_back(b.bigEnd(idim)+1);
                max_size[idim] = std::max(max_size[idim], b.length(idim));
            }
        }
        for (auto& c : cuts) {
            std::sort(c.begin(), c.end());
            c.erase(std::unique(c.begin(), c.end()), c.end());
        }
        auto cut_index = [&cuts] (int idim, int coord) {
            return static_cast<int>(std::lower_bound(cuts[idim].begin(), cuts[idim].end(), coord)
                                    - cuts[idim].begin());
        };
 
        // The z extents of the boxes over each piece, keyed by the piece's cut indices
        std::map<std::pair<int,int>, std::vector<std::pair<int,int>>> zruns_of_piece;
        for (const auto& b : boxes) {
            for (int ix = cut_index(0, b.smallEnd(0)); ix < cut_index(0, b.bigEnd(0)+1); ++ix) {
                for (int iy = cut_index(1, b.smallEnd(1)); iy < cut_index(1, b.bigEnd(1)+1); ++iy) {
                    zruns_of_piece[std::make_pair(ix,iy)].emplace_back(b.smallEnd(2), b.bigEnd(2));
                }
            }
        }
 
        // Over each piece, join the extents that touch into one box per run of cells
        BoxList bl_stack;
        for (auto& piece : zruns_of_piece) {
            const int ix = piece.first.first;
            const int iy = piece.first.second;
            auto& zruns = piece.second;
            std::sort(zruns.begin(), zruns.end());
            std::size_t n_first = 0;
            for (std::size_t n = 1; n <= zruns.size(); ++n) {
                if (n == zruns.size() || zruns[n].first != zruns[n-1].second+1) {
                    bl_stack.push_back(Box(IntVect(cuts[0][ix],     cuts[1][iy],     zruns[n_first].first),
                                           IntVect(cuts[0][ix+1]-1, cuts[1][iy+1]-1, zruns[n-1].second)));
                    n_first = n;
                }
            }
        }
 
        // Merging in x and y joins pieces with the same z extent, so it cannot stack them again
        bl_stack.simplify();
        max_size[2] = bl_stack.minimalBox().length(2);
        bl_stack.maxSize(max_size);
        bl_joined.join(bl_stack);
    }
 
    return BoxArray(std::move(bl_joined));
}

Referenced by SurfaceLayer::define_pblh_columns().

Here is the call graph for this function:

Here is the caller graph for this function:

make_areas()

void make_areas	(	const Geometry &	geom,
		MultiFab &	z_phys_nd,
		MultiFab &	ax,
		MultiFab &	ay,
		MultiFab &	az
	)

Computation of area fractions on faces

{
    const auto* dx = geom.CellSize();
    Real dzInv = one/dx[2];
 
    // Domain valid box (z_nd is nodal)
    const Box& domain = geom.Domain();
    int domlo_z = domain.smallEnd(2);
 
    // The z-faces are always full when using terrain-fitted coordinates
    az.setVal(1.0);
 
    //
    // x-areas
    //
#ifdef _OPENMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
    for ( MFIter mfi(ax, TilingIfNotGPU()); mfi.isValid(); ++mfi )
    {
        Box gbx = mfi.growntilebox(ax.nGrow());
        if (gbx.smallEnd(2) < domlo_z) {
            gbx.setSmall(2,domlo_z);
        }
 
        Array4<Real const> z_nd = z_phys_nd.const_array(mfi);
        Array4<Real      > ax_arr = ax.array(mfi);
        ParallelFor(gbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
               ax_arr(i, j, k) = Real(.5) * dzInv * (
                       z_nd(i,j,k+1) + z_nd(i,j+1,k+1) - z_nd(i,j,k) - z_nd(i,j+1,k));
        });
    }
 
    //
    // y-areas
    //
#ifdef _OPENMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
    for ( MFIter mfi(ay, TilingIfNotGPU()); mfi.isValid(); ++mfi )
    {
        Box gbx = mfi.growntilebox(ay.nGrow());
        if (gbx.smallEnd(2) < domlo_z) {
            gbx.setSmall(2,domlo_z);
        }
 
        Array4<Real const> z_nd = z_phys_nd.const_array(mfi);
        Array4<Real      > ay_arr = ay.array(mfi);
        ParallelFor(gbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
               ay_arr(i, j, k) = Real(.5) * dzInv * (
                       z_nd(i,j,k+1) + z_nd(i+1,j,k+1) - z_nd(i,j,k) - z_nd(i+1,j,k));
        });
    }
 
    ax.FillBoundary(geom.periodicity());
    ay.FillBoundary(geom.periodicity());
    az.FillBoundary(geom.periodicity());
}

Referenced by ERF::update_terrain_arrays().

Here is the call graph for this function:

Here is the caller graph for this function:

make_J()

void make_J	(	const Geometry &	geom,
		MultiFab &	z_phys_nd,
		MultiFab &	detJ_cc
	)

Computation of detJ at cell-center

{
    const auto *dx = geom.CellSize();
    Real dzInv = one/dx[2];
 
    // Domain valid box (z_nd is nodal)
    const Box& domain = geom.Domain();
    int domlo_z = domain.smallEnd(2);
 
    // Number of ghost cells
    int ngrow= detJ_cc.nGrow();
 
#ifdef _OPENMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
    for ( MFIter mfi(detJ_cc, TilingIfNotGPU()); mfi.isValid(); ++mfi )
    {
        Box gbx = mfi.growntilebox(ngrow);
        if (gbx.smallEnd(2) < domlo_z) {
            gbx.setSmall(2,domlo_z);
        }
        Array4<Real const> z_nd = z_phys_nd.const_array(mfi);
        Array4<Real      > detJ = detJ_cc.array(mfi);
        ParallelFor(gbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
           detJ(i, j, k) = Real(.25) * dzInv * (
                   z_nd(i,j,k+1) + z_nd(i+1,j,k+1) + z_nd(i,j+1,k+1) + z_nd(i+1,j+1,k+1)
                  -z_nd(i,j,k  ) - z_nd(i+1,j,k  ) - z_nd(i,j+1,k  ) - z_nd(i+1,j+1,k  ) );
        });
    }
    detJ_cc.FillBoundary(geom.periodicity());
}

Referenced by ERF::update_terrain_arrays().

Here is the call graph for this function:

Here is the caller graph for this function:

make_terrain_fitted_coords()

void make_terrain_fitted_coords	(	int	lev,
		const Geometry &	geom,
		MultiFab &	z_phys_nd,
		Vector< Real > const &	z_levels_h,
		GpuArray< ERF_BC, AMREX_SPACEDIM *2 > &	phys_bc_type,
		FineTerrain	fine_terrain,
		MultiFab const *	z_phys_interp
	)

Computation of the terrain grid from BTF, STF, or Sullivan TF model

{
    const Box& domain = geom.Domain();
 
    int domlo_z = domain.smallEnd(2);
    int domhi_z = domain.bigEnd(2) + 1;
 
    // Just in case ...
    z_phys_nd.setDomainBndry(bogus_large_value,0,1,geom);
 
    // ****************************************************************************
 
    if (lev == 0) {
        BoxArray ba(z_phys_nd.boxArray());
        bool all_boxes_touch_bottom = true;
        for (int i = 0; i < ba.size(); i++) {
            if (ba[i].smallEnd(2) != domlo_z) {
                all_boxes_touch_bottom = false;
            }
        }
 
        if (all_boxes_touch_bottom) {
            init_which_terrain_grid(lev, geom, z_phys_nd, z_levels_h, fine_terrain, z_phys_interp);
        } else {
 
            BoxArray ba_new(domain);
            ChopGrids2D(ba_new, domain, ParallelDescriptor::NProcs());
 
            DistributionMapping dm_new(ba_new);
            ba_new.surroundingNodes();
 
            MultiFab z_phys_nd_new(ba_new, dm_new, 1, z_phys_nd.nGrowVect());
 
            z_phys_nd_new.ParallelCopy(z_phys_nd,0,0,1,z_phys_nd.nGrowVect(),z_phys_nd.nGrowVect());
 
            init_which_terrain_grid(lev, geom, z_phys_nd_new, z_levels_h, fine_terrain, z_phys_interp);
 
            z_phys_nd.ParallelCopy(z_phys_nd_new,0,0,1,z_phys_nd.nGrowVect(),z_phys_nd.nGrowVect());
        }
    } else { // lev > 0
        //
        // BTF maps each column of a box between the box's lowest node and the domain top.
        // A box that sits on another box of this level must continue the columns of the box
        // below it: its own lowest node holds only the mesh interpolated from the coarse
        // level, so starting from there would make the mesh above it follow the coarse
        // terrain and depend on where the grids are split in z.  As on level 0 we therefore
        // build the mesh on boxes that hold whole columns and copy it back.  The bottom box
        // of a column, such as a refined patch aloft, still starts from its own lowest node.
        //
        ParmParse pp("erf");
        int terrain_smoothing = 0;
        pp.query("terrain_smoothing", terrain_smoothing);
 
        const ColumnGrids* cg = (terrain_smoothing == 0)
            ? &column_grids(lev, geom, amrex::convert(z_phys_nd.boxArray(), IntVect::TheCellVector()))
            : nullptr;
 
        if (cg == nullptr || !cg->stacked) {
            init_which_terrain_grid(lev, geom, z_phys_nd, z_levels_h, fine_terrain, z_phys_interp);
        } else {
            const IntVect& ng = z_phys_nd.nGrowVect();
            MultiFab z_phys_nd_col(cg->ba_col_nd, cg->dm_col, 1, ng);
 
            z_phys_nd_col.ParallelCopy(z_phys_nd,0,0,1,ng,ng);
 
            init_which_terrain_grid(lev, geom, z_phys_nd_col, z_levels_h, fine_terrain, z_phys_interp);
 
            // Copy the ghost nodes too, then give every node the value of a box that holds it
            // as valid, rather than one a neighbouring box computed for its ghost node
            z_phys_nd.ParallelCopy(z_phys_nd_col,0,0,1,ng,ng);
            z_phys_nd.ParallelCopy(z_phys_nd_col,0,0,1,IntVect(0),ng);
        }
    }
 
    //
    // A node on the high side of a periodic direction is the same node as its image on the
    // low side, and a box that holds one as valid is not touched by FillBoundary -- so give
    // both copies one value, or the faces on either side of the node see different metrics
    // and the mass flux through them differs.  make_terrain_surface_periodic already does
    // this for a surface built by init_terrain_surface (and that is what the immersed
    // boundary reads), but not for one that comes from wrfinput, met_em, a restart or the
    // interpolation from a coarser level; this covers all of them, and the whole column.
    // The copy kept is that of the lower-numbered box, which for the usual grids (and
    // always within one box) is the low side, as in make_terrain_surface_periodic.
    //
    // Only the lateral periodicity: in a domain that is also periodic in z, the top node of
    // a column is the image of the bottom one, but its height is not -- syncing them would
    // give the top node the bottom's height and the top cell a negative thickness.
    //
    if (geom.isPeriodic(0) || geom.isPeriodic(1)) {
        const Box& dom = geom.Domain();
        const Periodicity lateral_period(IntVect(AMREX_D_DECL(geom.isPeriodic(0) ? dom.length(0) : 0,
                                                              geom.isPeriodic(1) ? dom.length(1) : 0,
                                                              0)));
        z_phys_nd.OverrideSync(lateral_period);
    }
 
    //
    // Fill ghost layers and corners (including periodic) -- no matter what level
    //
    z_phys_nd.FillBoundary(geom.periodicity());
 
    if (phys_bc_type[Orientation(0,Orientation::low )] == ERF_BC::symmetry ||
        phys_bc_type[Orientation(0,Orientation::high)] == ERF_BC::symmetry ||
        phys_bc_type[Orientation(1,Orientation::low )] == ERF_BC::symmetry ||
        phys_bc_type[Orientation(1,Orientation::high)] == ERF_BC::symmetry) {
 
        const auto& dom_lo = lbound(convert(domain,IntVect(1,1,1)));
        const auto& dom_hi = ubound(convert(domain,IntVect(1,1,1)));
 
        //
        // Mirror the whole ghost region, not just its first layer: the base state carries
        // ComputeGhostCells+1 layers and the state stencils reach several cells past the
        // boundary, so mirroring only the first would leave the rest of the ghost region
        // holding whatever the terrain construction extrapolated into it.
        //
        const IntVect& ngv = z_phys_nd.nGrowVect();
 
        for (MFIter mfi(z_phys_nd,true); mfi.isValid(); ++mfi) {
            //
            // Whether this tile reaches a domain boundary is a question about its valid
            // nodes, so ask it of the tilebox.  The grown tilebox begins a ghost width short
            // of the boundary and so never equals it -- asking it of the grown box is what
            // kept this mirroring from ever running.  The grown box is still what we sweep,
            // so that the ghost rows in the other two directions are covered as well.
            //
            const Box& vbx = mfi.tilebox();
            const Box& bx  = mfi.growntilebox();
            const Array4< Real> z_nd_arr = z_phys_nd.array(mfi);
            if (phys_bc_type[Orientation(0,Orientation::low)] == ERF_BC::symmetry && vbx.smallEnd(0) == dom_lo.x) {
                const int ng = ngv[0];
                ParallelFor(makeSlab(bx,0,1), [=] AMREX_GPU_DEVICE (int , int j, int k)
                {
                    for (int n = 1; n <= ng; ++n) {
                        z_nd_arr(dom_lo.x-n,j,k) = z_nd_arr(dom_lo.x+n,j,k);
                    }
                });
            }
            if (phys_bc_type[Orientation(0,Orientation::high)] == ERF_BC::symmetry && vbx.bigEnd(0) == dom_hi.x) {
                const int ng = ngv[0];
                ParallelFor(makeSlab(bx,0,1), [=] AMREX_GPU_DEVICE (int , int j, int k)
                {
                    for (int n = 1; n <= ng; ++n) {
                        z_nd_arr(dom_hi.x+n,j,k) = z_nd_arr(dom_hi.x-n,j,k);
                    }
                });
            }
            if (phys_bc_type[Orientation(1,Orientation::low)] == ERF_BC::symmetry && vbx.smallEnd(1) == dom_lo.y) {
                const int ng = ngv[1];
                ParallelFor(makeSlab(bx,1,1), [=] AMREX_GPU_DEVICE (int i, int  , int k)
                {
                    for (int n = 1; n <= ng; ++n) {
                        z_nd_arr(i,dom_lo.y-n,k) = z_nd_arr(i,dom_lo.y+n,k);
                    }
                });
            }
            if (phys_bc_type[Orientation(1,Orientation::high)] == ERF_BC::symmetry && vbx.bigEnd(1) == dom_hi.y) {
                const int ng = ngv[1];
                ParallelFor(makeSlab(bx,1,1), [=] AMREX_GPU_DEVICE (int i, int  , int k)
                {
                    for (int n = 1; n <= ng; ++n) {
                        z_nd_arr(i,dom_hi.y+n,k) = z_nd_arr(i,dom_hi.y-n,k);
                    }
                });
            }
        }
    }
 
    //********************************************************************************
    // Populate domain boundary cells in z-direction
    //********************************************************************************
    int ngrow = z_phys_nd.nGrow();
 
    for ( MFIter mfi(z_phys_nd, TilingIfNotGPU()); mfi.isValid(); ++mfi )
    {
        // Only set values above top of domain if this box reaches that far
        Box nd_bx = mfi.tilebox();
 
        // Note that domhi_z is already nodal in the z-direction
        if (nd_bx.bigEnd(2) >= domhi_z) {
            // Grown box with no z range
            Box bx_zhi = mfi.growntilebox(ngrow);
            bx_zhi.setSmall(2,domhi_z+1);
            Array4<Real> const& z_arr = z_phys_nd.array(mfi);
 
            // Extrapolate top layer
            ParallelFor(bx_zhi, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
                z_arr(i,j,k) = z_arr(i,j,domhi_z)
                             + (k-domhi_z) * (z_arr(i,j,domhi_z) - z_arr(i,j,domhi_z-1));
            });
        }
    }
} // make_terrain_fitted_coords

Referenced by ERF::init_zphys(), and ERF::remake_zphys().

Here is the call graph for this function:

Here is the caller graph for this function:

make_zcc()

void make_zcc	(	const Geometry &	geom,
		MultiFab &	z_phys_nd,
		MultiFab &	z_phys_cc
	)

Computation of z_phys at cell-center

{
#ifdef _OPENMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
    for ( MFIter mfi(z_phys_cc, TilingIfNotGPU()); mfi.isValid(); ++mfi )
    {
        Box gbx = mfi.growntilebox();
        Array4<Real const> z_nd = z_phys_nd.const_array(mfi);
        Array4<Real      > z_cc = z_phys_cc.array(mfi);
        ParallelFor(gbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
           z_cc(i, j, k) = Real(.125) * ( z_nd(i,j,k  ) + z_nd(i+1,j,k  ) + z_nd(i,j+1,k  ) + z_nd(i+1,j+1,k  )
                                         +z_nd(i,j,k+1) + z_nd(i+1,j,k+1) + z_nd(i,j+1,k+1) + z_nd(i+1,j+1,k+1) );
       });
    }
    z_phys_cc.FillBoundary(geom.periodicity());
}

Referenced by ERF::post_timestep(), and ERF::update_terrain_arrays().

Here is the call graph for this function:

Here is the caller graph for this function:

validate_flat_terrain()

void validate_flat_terrain	(	int	lev,
		const MultiFab &	z_phys_nd,
		const Vector< Real > &	z_levels_h
	)

{
    // The comparison is against heights that have been round-tripped through z_phys_nd, so the
    // threshold has to track the working precision: a hard-wired 1e-10 is far below single-
    // precision round-off at domain-top magnitudes and would abort on a genuinely flat mesh.
    // Keep the original (double-precision) threshold, floored at a few hundred ulps.
    const Real scale = std::max({Real(1.0), std::abs(z_levels_h.front()), std::abs(z_levels_h.back())});
    const Real rel_tol = std::max(Real(1.0e-10), Real(100.0) * std::numeric_limits<Real>::epsilon());
    const Real tolerance = rel_tol * scale;
 
    MultiFab expected(z_phys_nd.boxArray(), z_phys_nd.DistributionMap(), 1, 0);
    Gpu::DeviceVector<Real> z_levels_d(z_levels_h.size());
    Gpu::copy(Gpu::hostToDevice, z_levels_h.begin(), z_levels_h.end(), z_levels_d.begin());
    const Real* z_levels = z_levels_d.data();
 
    for (MFIter mfi(expected); mfi.isValid(); ++mfi) {
        const Box& bx = mfi.validbox();
        auto const& expected_arr = expected.array(mfi);
        ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
            expected_arr(i,j,k) = z_levels[k];
        });
    }
 
    MultiFab error(expected.boxArray(), expected.DistributionMap(), 1, 0);
    MultiFab::Copy(error, expected, 0, 0, 1, 0);
    MultiFab::Subtract(error, z_phys_nd, 0, 0, 1, 0);
    const Real max_error = error.norminf(0, 0, false);
 
    if (max_error > tolerance) {
        Abort("erf.flat_terrain = true, but z_phys_nd is not horizontally "
              "flat or does not match the configured vertical spacing at "
              "level " + std::to_string(lev) + ".");
    }
}

Referenced by ERF::init_zphys().

Here is the call graph for this function:

Here is the caller graph for this function:

validate_flat_top()

void validate_flat_top	(	int	lev,
		const Geometry &	geom,
		const MultiFab &	z_phys_nd
	)

Abort unless the top of the domain is flat, i.e. z_phys_nd is the same at every node of the top slab. The top boundary is a horizontal plane: erf_slow_rhs_pre sets Omega = rho w on the top face, the anelastic projection converts it with OmegaFromW, and the two only agree when the terrain slopes vanish there. Every mesh builder pins the top node to z_top, so this guards against a mesh (custom, read from a file or a checkpoint) that does not.

Parameters

[in]	lev	Level index.
[in]	geom	Geometry information.
[in]	z_phys_nd	Nodal physical height field.

{
    const int k_top = geom.Domain().bigEnd(2) + 1;
 
    ReduceOps<ReduceOpMin, ReduceOpMax> reduce_op;
    ReduceData<Real, Real> reduce_data(reduce_op);
    using ReduceTuple = typename decltype(reduce_data)::Type;
 
    for (MFIter mfi(z_phys_nd); mfi.isValid(); ++mfi) {
        const Box& bx = mfi.validbox();
        if (bx.bigEnd(2) < k_top) { continue; }
        auto const& z_arr = z_phys_nd.const_array(mfi);
        reduce_op.eval(makeSlab(bx,2,k_top), reduce_data,
        [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple {
            return {z_arr(i,j,k), z_arr(i,j,k)};
        });
    }
 
    ReduceTuple hv = reduce_data.value(reduce_op);
    Real z_top_min = amrex::get<0>(hv);
    Real z_top_max = amrex::get<1>(hv);
    ParallelDescriptor::ReduceRealMin(z_top_min);
    ParallelDescriptor::ReduceRealMax(z_top_max);
 
    // No box at this level reaches the top of the domain
    if (z_top_min > z_top_max) { return; }
 
    // Same precision-aware tolerance as validate_flat_terrain
    const Real scale = std::max({Real(1.0), std::abs(z_top_min), std::abs(z_top_max)});
    const Real rel_tol = std::max(Real(1.0e-10), Real(100.0) * std::numeric_limits<Real>::epsilon());
 
    if (z_top_max - z_top_min > rel_tol * scale) {
        Print() << "Top of the domain at level " << lev << " ranges from "
                << z_top_min << " to " << z_top_max << std::endl;
        Abort("The top of the domain must be flat: z_phys_nd varies along the top boundary at level "
              + std::to_string(lev) + ".");
    }
}

Referenced by ERF::update_terrain_arrays().

Here is the caller graph for this function:

which_fine_terrain()

FineTerrain which_fine_terrain

(

)

Translate erf.amr_terrain_refinement into a FineTerrain.

{
    ParmParse pp("erf");
    std::string amr_terrain_refinement = "interpolate";
    pp.query("amr_terrain_refinement", amr_terrain_refinement);
 
    if (amr_terrain_refinement == "interpolate") {
        return FineTerrain::Interpolate;
    } else if (amr_terrain_refinement == "transform") {
        return FineTerrain::Transform;
    }
 
    Abort("erf.amr_terrain_refinement = " + amr_terrain_refinement +
          " is not recognized; it must be \"interpolate\" or \"transform\"");
    return FineTerrain::Interpolate;
}

Referenced by ERF::init_zphys(), and ERF::remake_zphys().

Here is the call graph for this function:

Here is the caller graph for this function: