ERF_LatLonMap.cpp File Reference

#include "ERF_LatLonMap.H"
#include "ERF_EarthRelativeWind.H"
#include "ERF_NumericalConstants.H"
#include <algorithm>
#include <AMReX_Arena.H>
#include <AMReX_GpuDevice.H>
#include <AMReX_MFIter.H>
#include <AMReX_ParallelDescriptor.H>

Include dependency graph for ERF_LatLonMap.cpp:

Functions
void	pack_latlon (MultiFab &latlon, const MultiFab &lat, const MultiFab &lon)
void	grid_rotation_from_latlon (const Array4< const Real > &ll, const Box &dom, int i, int j, Real &cos_alpha, Real &sin_alpha)
void	fill_grid_rotation (MultiFab &cos_alpha, MultiFab &sin_alpha, const MultiFab &lat, const MultiFab &lon, const Geometry &geom)
LatLonStatus	locate_latlon_on_grid (const Array4< const Real > &ll, const Box &dom, const GpuArray< Real, AMREX_SPACEDIM > &problo, const GpuArray< Real, AMREX_SPACEDIM > &dx, Real req_lat, Real req_lon, LatLonLocation &loc)

Function Documentation

fill_grid_rotation()

void fill_grid_rotation	(	MultiFab &	cos_alpha,
		MultiFab &	sin_alpha,
		const MultiFab &	lat,
		const MultiFab &	lon,
		const Geometry &	geom
	)

{
    const Box& dom = geom.Domain();
    const int ilo = dom.smallEnd(0);
    const int ihi = dom.bigEnd(0);
 
    // Ghost cells outside the domain are left at zero: the stencil never reads them
    const IntVect ng(1,1,0);
    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(lat.nGrowVect().allGE(ng) && lon.nGrowVect().allGE(ng),
        "fill_grid_rotation: lat/lon need one ghost cell in x and y");
    MultiFab ll(cos_alpha.boxArray(), cos_alpha.DistributionMap(), 2, ng);
    ll.setVal(Real(0.0));
    const Periodicity period(IntVect(geom.isPeriodic(0) ? dom.length(0) : 0,
                                     geom.isPeriodic(1) ? dom.length(1) : 0, 0));
    ll.ParallelCopy(lat, 0, 0, 1, ng, ng, period);
    ll.ParallelCopy(lon, 0, 1, 1, ng, ng, period);
 
    for (MFIter mfi(cos_alpha, TilingIfNotGPU()); mfi.isValid(); ++mfi) {
        const Box& bx = mfi.tilebox();
        const Array4<const Real>& ll_arr = ll.const_array(mfi);
        const Array4<Real>& c_arr = cos_alpha.array(mfi);
        const Array4<Real>& s_arr = sin_alpha.array(mfi);
        ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
        {
            erf_earth_wind::grid_rotation_at(ll_arr, ilo, ihi, i, j, c_arr(i,j,k), s_arr(i,j,k));
        });
    }
}

Referenced by ERF::FillPlot2DVars(), and TopographicShading::initialize_geometry().

Here is the call graph for this function:

Here is the caller graph for this function:

grid_rotation_from_latlon()

void grid_rotation_from_latlon	(	const Array4< const Real > &	ll,
		const Box &	dom,
		int	i,
		int	j,
		Real &	cos_alpha,
		Real &	sin_alpha
	)

{
    erf_earth_wind::grid_rotation_at(ll, dom.smallEnd(0), dom.bigEnd(0), i, j,
                                     cos_alpha, sin_alpha);
}

Referenced by locate_latlon_on_grid(), and LatLonMap::rotation_at().

Here is the call graph for this function:

Here is the caller graph for this function:

locate_latlon_on_grid()

LatLonStatus locate_latlon_on_grid	(	const Array4< const Real > &	ll,
		const Box &	dom,
		const GpuArray< Real, AMREX_SPACEDIM > &	problo,
		const GpuArray< Real, AMREX_SPACEDIM > &	dx,
		Real	req_lat,
		Real	req_lon,
		LatLonLocation &	loc
	)

{
    const int ilo = dom.smallEnd(0), ihi = dom.bigEnd(0);
    const int jlo = dom.smallEnd(1), jhi = dom.bigEnd(1);
 
    const Real cosfac = std::cos(req_lat * PI / Real(180.0));
    auto dist2 = [&](int i, int j) {
        const Real dlat = ll(i,j,0,0) - req_lat;
        const Real dlon = wrap_longitude_difference(ll(i,j,0,1) - req_lon) * cosfac;
        return dlat*dlat + dlon*dlon;
    };
 
    int  bi = ilo, bj = jlo;
    Real best = dist2(ilo,jlo);
    for (int j = jlo; j <= jhi; ++j) {
        for (int i = ilo; i <= ihi; ++i) {
            const Real d = dist2(i,j);
            if (d < best) { best = d; bi = i; bj = j; }
        }
    }
 
    loc.near_lat = ll(bi,bj,0,0);
    loc.near_lon = ll(bi,bj,0,1);
 
    // One linear solve in index space inverts the map: over a cell the
    // projection is linear to well below a cell width.
    const GradientStencil s = gradient_stencil(dom, bi, bj);
 
    const Real dlat_di = (ll(s.ip,bj,0,0) - ll(s.im,bj,0,0)) * s.inv_di;
    const Real dlon_di = wrap_longitude_difference(ll(s.ip,bj,0,1) - ll(s.im,bj,0,1)) * s.inv_di;
    const Real dlat_dj = (ll(bi,s.jp,0,0) - ll(bi,s.jm,0,0)) * s.inv_dj;
    const Real dlon_dj = wrap_longitude_difference(ll(bi,s.jp,0,1) - ll(bi,s.jm,0,1)) * s.inv_dj;
 
    const Real rlat = req_lat - ll(bi,bj,0,0);
    const Real rlon = wrap_longitude_difference(req_lon - ll(bi,bj,0,1));
    Real di = Real(0.0), dj = Real(0.0);
    if (!invert_latlon_offset(rlat, rlon, dlat_di, dlon_di, dlat_dj, dlon_dj, di, dj)) {
        return LatLonStatus::Degenerate;
    }
 
    // The solve places the point relative to the nearest grid point, so more
    // than a cell away means the request cannot be trusted: either it is
    // outside the domain, or the map is too poorly conditioned near it for one
    // linear solve.
    if (std::abs(di) > Real(1.0) || std::abs(dj) > Real(1.0)) {
        return LatLonStatus::TooFar;
    }
 
    loc.x   = problo[0] + (Real(bi) + Real(0.5) + di) * dx[0];
    loc.y   = problo[1] + (Real(bj) + Real(0.5) + dj) * dx[1];
    loc.lat = ll(bi,bj,0,0) + dlat_di*di + dlat_dj*dj;
    // Resolving a point from a nearest grid point on the other side of the
    // antimeridian steps the longitude just past +-180; report the longitude
    // that names the place, not the one that runs off the end of the range
    loc.lon = wrap_longitude_difference(ll(bi,bj,0,1) + dlon_di*di + dlon_dj*dj);
 
    grid_rotation_from_latlon(ll, dom, bi, bj, loc.cos_alpha, loc.sin_alpha);
 
    return LatLonStatus::Ok;
}

Referenced by LatLonMap::locate().

Here is the call graph for this function:

Here is the caller graph for this function:

pack_latlon()

void pack_latlon	(	MultiFab &	latlon,
		const MultiFab &	lat,
		const MultiFab &	lon
	)

{
    for (MFIter mfi(latlon, TilingIfNotGPU()); mfi.isValid(); ++mfi) {
        const Box& bx = mfi.tilebox();
        const Array4<Real>&       ll  = latlon.array(mfi);
        const Array4<const Real>& lat_arr = lat.const_array(mfi);
        const Array4<const Real>& lon_arr = lon.const_array(mfi);
        ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
        {
            ll(i,j,k,0) = lat_arr(i,j,0);
            ll(i,j,k,1) = lon_arr(i,j,0);
        });
    }
}

Referenced by LatLonMap::LatLonMap().

Here is the call graph for this function:

Here is the caller graph for this function: