ERF_LatLonMap.cpp File Reference

#include "ERF_LatLonMap.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)
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

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
	)

{
    const GradientStencil s = gradient_stencil(dom, i, j);
 
    // Displacement on the sphere per cell in i: east and north components, in
    // degrees of latitude (the common factor, the Earth's radius, cancels)
    const Real lat   = ll(i,j,0,0);
    const Real north = (ll(s.ip,j,0,0) - ll(s.im,j,0,0)) * s.inv_di;
    const Real east  = wrap_longitude_difference(ll(s.ip,j,0,1) - ll(s.im,j,0,1)) * s.inv_di
                     * std::cos(lat * PI / Real(180.0));
    const Real norm  = std::sqrt(east*east + north*north);
 
    if (norm > Real(0.0)) {
        cos_alpha = east  / norm;
        sin_alpha = north / norm;
    } else {
        // A domain one cell wide in x has no i direction to measure
        cos_alpha = Real(1.0);
        sin_alpha = Real(0.0);
    }
}

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: