ERF_InitTurbPert.cpp File Reference

#include <ERF.H>
#include <ERF_TerrainSurfaceSlab.H>
#include <ERF_PerturbationWind.H>
#include <AMReX_MultiFabUtil.H>
#include <ERF_TileNoZ.H>
#include <ERF_ProbCommon.H>

Include dependency graph for ERF_InitTurbPert.cpp:

Functions
PerturbationWindSample	sample_perturbation_wind (const Box &region, Real height_agl, const Geometry &geometry, const MultiFab &state, const MultiFab &xvel, const MultiFab &yvel, const MultiFab *z_cc, const MultiFab *z_nd)

Function Documentation

sample_perturbation_wind()

PerturbationWindSample sample_perturbation_wind	(	const Box &	region,
		Real	height_agl,
		const Geometry &	geometry,
		const MultiFab &	state,
		const MultiFab &	xvel,
		const MultiFab &	yvel,
		const MultiFab *	z_cc,
		const MultiFab *	z_nd
	)

{
    // Keep only the edge strip, with the source grids' ownership. A vertical
    // halo exchanges cell-centered winds/heights across stacked boxes without
    // allocating a full-volume copy or depending on source ghost values.
    BoxList boxes;
    Vector<int> source_indices, owners;
    const auto& source_ba = state.boxArray();
    const auto& source_owners = state.DistributionMap().ProcessorMap();
    for (int b = 0; b < source_ba.size(); ++b) {
        const Box bx = source_ba[b] & region;
        if (!bx.ok()) { continue; }
        boxes.push_back(bx);
        source_indices.push_back(b);
        owners.push_back(source_owners[b]);
    }
    if (boxes.isEmpty()) { return {}; }
    const BoxArray ba(boxes);
    const DistributionMapping dm(owners);
    // Components are u, v, cell-center height, bottom height, top height (AGL).
    MultiFab samples(ba, dm, 5, IntVect(0,0,1));
    samples.setVal(bogus_large_value);
    const int klo = geometry.Domain().smallEnd(2);
    const auto dx = geometry.CellSizeArray();
    auto ground = perturbation_ground(z_nd, samples, geometry);
    for (MFIter mfi(samples); mfi.isValid(); ++mfi) {
        const int src = source_indices[mfi.index()];
        const auto ua = xvel.const_array(src);
        const auto va = yvel.const_array(src);
        const auto cc = z_cc ? z_cc->const_array(src) : Array4<const Real>{};
        const auto nd = z_nd ? z_nd->const_array(src) : Array4<const Real>{};
        const auto ga = ground ? ground->const_array(mfi) : Array4<const Real>{};
        const auto a = samples.array(mfi);
        ParallelFor(mfi.validbox(), [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
            Real zlo = (k-klo)*dx[2];
            Real zhi = zlo + dx[2];
            Real z = myhalf*(zlo+zhi);
            if (nd) {
                const Real zs = fourth*(ga(i,j,klo) + ga(i+1,j,klo) +
                                        ga(i,j+1,klo) + ga(i+1,j+1,klo));
                zlo = fourth*(nd(i,j,k) + nd(i+1,j,k) +
                              nd(i,j+1,k) + nd(i+1,j+1,k)) - zs;
                zhi = fourth*(nd(i,j,k+1) + nd(i+1,j,k+1) +
                              nd(i,j+1,k+1) + nd(i+1,j+1,k+1)) - zs;
                z = cc ? cc(i,j,k) - zs : myhalf*(zlo+zhi);
            }
            a(i,j,k,0) = myhalf*(ua(i,j,k)+ua(i+1,j,k));
            a(i,j,k,1) = myhalf*(va(i,j,k)+va(i,j+1,k));
            a(i,j,k,2) = z;
            a(i,j,k,3) = zlo;
            a(i,j,k,4) = zhi;
        });
    }
    samples.FillBoundary(geometry.periodicity());
 
    Real sums[3] = {zero, zero, zero};
    for (MFIter mfi(samples); mfi.isValid(); ++mfi) {
        const auto a = samples.const_array(mfi);
        ReduceOps<ReduceOpSum, ReduceOpSum, ReduceOpSum> op;
        ReduceData<Real, Real, Real> data(op);
        using Tuple = typename decltype(data)::Type;
        op.eval(mfi.validbox(), data, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept -> Tuple {
            const Real z = a(i,j,k,2);
            const Real znext = a(i,j,k+1,2);
            Real u = a(i,j,k,0);
            Real v = a(i,j,k,1);
            if (znext != bogus_large_value && znext > z && height_agl >= z && height_agl < znext) {
                const Real weight = (height_agl-z)/(znext-z);
                u += weight*(a(i,j,k+1,0)-u);
                v += weight*(a(i,j,k+1,1)-v);
            } else if ((k == klo && height_agl >= a(i,j,k,3) && height_agl < z) ||
                       (znext == bogus_large_value && height_agl >= z && height_agl <= a(i,j,k,4))) {
                // Constant extension within the first/last half-cell only.
            } else {
                return {Real(zero), Real(zero), Real(zero)};
            }
            return {u, v, Real(one)};
        });
        const auto values = data.value(op);
        sums[0] += amrex::get<0>(values);
        sums[1] += amrex::get<1>(values);
        sums[2] += amrex::get<2>(values);
    }
    ParallelAllReduce::Sum(sums, 3, ParallelContext::CommunicatorSub());
    Box footprint = region;
    footprint.setRange(2, klo);
    if (sums[2] != Real(footprint.numPts()) || sums[2] == zero) {
        // Report the actual sampling failure: an incomplete count does not
        // necessarily mean the ancestor's domain is too shallow. These
        // collective reductions only run when a source level is rejected.
        const Real min_height = samples.min(2);
        const Real max_height = samples.max(2);
        Print() << "Automatic perturbation wind sampling: region=" << region
                << " target_height_agl=" << height_agl << " m"
                << " samples=" << sums[2]
                << " expected_columns=" << footprint.numPts()
                << " cell_height_agl_min=" << min_height << " m"
                << " cell_height_agl_max=" << max_height << " m"
                << " source_boxes=" << ba.size() << "\n";
        return {};
    }
    const Real u = sums[0]/sums[2];
    const Real v = sums[1]/sums[2];
    if (!std::isfinite(u) || !std::isfinite(v)) {
        Abort("Automatic perturbations: non-finite wind at the CPM sampling height");
    }
    return {u, v, true};
}

Referenced by ERF::turbPert_update_automatic(), and ERF::turbPert_update_uniform_wind().

Here is the call graph for this function:

Here is the caller graph for this function: