ERF_WeatherDataInterpolation.cpp File Reference

#include <filesystem>
#include <stdexcept>
#include "ERF.H"
#include "ERF_ReadCustomBinaryIC.H"
#include "ERF_Interpolation_Bilinear.H"

Include dependency graph for ERF_WeatherDataInterpolation.cpp:

Macros
#define	ERF_WEATHERDATAINTERPOLATION_H_

Enumerations
enum class	BoundType { Lo , Hi }
enum class	MultiFabType { CC , NC }

Functions
void	fill_weather_data_multifab (MultiFab &mf, const Geometry &geom_weather, const int nx, const int ny, const int nz, const Vector< Real > &latvec_h, const Vector< Real > &lonvec_h, const Vector< Real > &zvec_h, const Vector< Real > &rho_h, const Vector< Real > &uvel_h, const Vector< Real > &vvel_h, const Vector< Real > &wvel_h, const Vector< Real > &theta_h, const Vector< Real > &qv_h, const Vector< Real > &qc_h, const Vector< Real > &qr_h)
void	PlotMultiFab (const MultiFab &mf, const Geometry &geom_mf, const std::string plotfilename, MultiFabType mftype)
IntVect	find_bound_idx (const Real &x, const Real &y, const Real &z, const BoxList &bl_weather, const Geometry &geom_weather, BoundType bound_type)

Macro Definition Documentation

ERF_WEATHERDATAINTERPOLATION_H_

#define ERF_WEATHERDATAINTERPOLATION_H_

Enumeration Type Documentation

BoundType

enum BoundType

strong

Enumerator
Lo
Hi

{ Lo, Hi };

MultiFabType

enum MultiFabType

strong

Enumerator
CC
NC

{ CC, NC };

Function Documentation

fill_weather_data_multifab()

void fill_weather_data_multifab	(	MultiFab &	mf,
		const Geometry &	geom_weather,
		const int	nx,
		const int	ny,
		const int	nz,
		const Vector< Real > &	latvec_h,
		const Vector< Real > &	lonvec_h,
		const Vector< Real > &	zvec_h,
		const Vector< Real > &	rho_h,
		const Vector< Real > &	uvel_h,
		const Vector< Real > &	vvel_h,
		const Vector< Real > &	wvel_h,
		const Vector< Real > &	theta_h,
		const Vector< Real > &	qv_h,
		const Vector< Real > &	qc_h,
		const Vector< Real > &	qr_h
	)

{
    const int nx_d = nx;
    const int ny_d = ny;
    const int nz_d = nz;
    const int tot_size = nx*ny*nz;
 
    amrex::Gpu::DeviceVector<Real> latvec_d(nx*ny), lonvec_d(nx*ny), zvec_d(nz);
    amrex::Gpu::DeviceVector<Real> rho_d(tot_size), uvel_d(tot_size),
                                   vvel_d(tot_size), wvel_d(tot_size), theta_d(tot_size),
                                   qv_d(tot_size), qc_d(tot_size), qr_d(tot_size);
 
    amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, latvec_h.begin(), latvec_h.end(), latvec_d.begin());
    amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, lonvec_h.begin(), lonvec_h.end(), lonvec_d.begin());
    amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, zvec_h.begin(), zvec_h.end(), zvec_d.begin());
    amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, rho_h.begin(), rho_h.end(), rho_d.begin());
    amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, theta_h.begin(), theta_h.end(), theta_d.begin());
    amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, uvel_h.begin(), uvel_h.end(), uvel_d.begin());
    amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, vvel_h.begin(), vvel_h.end(), vvel_d.begin());
    amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, wvel_h.begin(), wvel_h.end(), wvel_d.begin());
    amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, qv_h.begin(), qv_h.end(), qv_d.begin());
    amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, qc_h.begin(), qc_h.end(), qc_d.begin());
    amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, qr_h.begin(), qr_h.end(), qr_d.begin());
 
    Real* latvec_d_ptr  = latvec_d.data();
    Real* lonvec_d_ptr  = lonvec_d.data();
    Real* zvec_d_ptr  = zvec_d.data();
    Real* rho_d_ptr   = rho_d.data();
    Real* uvel_d_ptr  = uvel_d.data();
    Real* vvel_d_ptr  = vvel_d.data();
    Real* wvel_d_ptr  = wvel_d.data();
    Real* theta_d_ptr = theta_d.data();
    Real* qv_d_ptr = qv_d.data();
    Real* qc_d_ptr = qc_d.data();
    Real* qr_d_ptr = qr_d.data();
 
    const auto prob_lo  = geom_weather.ProbLoArray();
    const auto dx       = geom_weather.CellSizeArray();
 
    for (MFIter mfi(mf, TilingIfNotGPU()); mfi.isValid(); ++mfi)
    {
        const Box& bx = mfi.nodaltilebox();
        Array4<Real> const& arr = mf.array(mfi);
 
        ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
        {
            const Real z = prob_lo[2] + k * dx[2];
            int kloc = -1;
            for (int kk = 0; kk < nz_d; kk++) {
                 if (zvec_d_ptr[kk] > z) {
                     kloc = kk;
                     break;
                 }
            }
            // Replace this with logic using your coarse input data
            int idx1 = get_single_index(i,j,kloc-1,nx_d,ny_d);
            int idx2 = get_single_index(i,j,kloc,nx_d,ny_d);
            Real dz = zvec_d_ptr[kloc] - zvec_d_ptr[kloc-1];
            Real fac = (z-zvec_d_ptr[kloc-1])/dz;
            arr(i,j,k,0) = rho_d_ptr[idx1]   + fac*(rho_d_ptr[idx2]-rho_d_ptr[idx1]);
            arr(i,j,k,1) = uvel_d_ptr[idx1]  + fac*(uvel_d_ptr[idx2]-uvel_d_ptr[idx1]);
            arr(i,j,k,2) = vvel_d_ptr[idx1]  + fac*(vvel_d_ptr[idx2]-vvel_d_ptr[idx1]);
            arr(i,j,k,3) = wvel_d_ptr[idx1]  + fac*(wvel_d_ptr[idx2]-wvel_d_ptr[idx1]);
            arr(i,j,k,4) = theta_d_ptr[idx1] + fac*(theta_d_ptr[idx2]-theta_d_ptr[idx1]);
            arr(i,j,k,5) = qv_d_ptr[idx1]    + fac*(qv_d_ptr[idx2]-qv_d_ptr[idx1]);
            arr(i,j,k,6) = qc_d_ptr[idx1]    + fac*(qc_d_ptr[idx2]-qc_d_ptr[idx1]);
            arr(i,j,k,7) = qr_d_ptr[idx1]    + fac*(qr_d_ptr[idx2]-qr_d_ptr[idx1]);
            idx1 = get_single_index(i,j,0,nx_d,ny_d);
            arr(i,j,k,8) = latvec_d_ptr[idx1];
            arr(i,j,k,9) = lonvec_d_ptr[idx1];
        });
    }
}

Referenced by ERF::FillWeatherDataMultiFab().

Here is the call graph for this function:

Here is the caller graph for this function:

find_bound_idx()

IntVect find_bound_idx	(	const Real &	x,
		const Real &	y,
		const Real &	z,
		const BoxList &	bl_weather,
		const Geometry &	geom_weather,
		BoundType	bound_type
	)

{
    const auto prob_lo_weather  = geom_weather.ProbLoArray();
    const auto dx_weather       = geom_weather.CellSizeArray();
 
    int i, j, k;
 
    if (bound_type == BoundType::Lo) {
        i = static_cast<int>(std::floor((x - prob_lo_weather[0]) / dx_weather[0]));
        j = static_cast<int>(std::floor((y - prob_lo_weather[1]) / dx_weather[1]));
        k = static_cast<int>(std::floor((z - prob_lo_weather[2]) / dx_weather[2]));
    } else { // BoundType::Hi
        i = static_cast<int>(std::ceil((x - prob_lo_weather[0]) / dx_weather[0]));
        j = static_cast<int>(std::ceil((y - prob_lo_weather[1]) / dx_weather[1]));
        k = static_cast<int>(std::ceil((z - prob_lo_weather[2]) / dx_weather[2]));
    }
 
    IntVect idx(i, j, k);
 
    for (const auto& b : bl_weather) {
        if (b.contains(idx)) {
            return idx;
        }
    }
 
    amrex::Abort("Bound index not found in any box in BoxList!");
    return IntVect::TheZeroVector(); // unreachable if Abort
}

Referenced by ERF::InterpWeatherDataOntoMesh().

Here is the caller graph for this function:

PlotMultiFab()

void PlotMultiFab	(	const MultiFab &	mf,
		const Geometry &	geom_mf,
		const std::string	plotfilename,
		MultiFabType	mftype
	)

{
 
    Vector<std::string> varnames = {
    "rho", "uvel", "vvel", "wvel", "theta", "qv", "qc", "qr", "latitude", "longitude"
    }; // Customize variable names
 
    const Real time = 0.0;
 
 
    // Assume weather_mf is nodal in all directions
    if(mftype == MultiFabType::NC) {
        BoxArray cba = mf.boxArray();
        cba = amrex::convert(mf.boxArray(), IntVect::TheCellVector());
 
        MultiFab cc_mf(cba, mf.DistributionMap(),
               mf.nComp(), 0);
 
        amrex::average_node_to_cellcenter(cc_mf, 0, mf, 0, mf.nComp());
 
        WriteSingleLevelPlotfile(
            plotfilename,
            cc_mf,
            varnames,
            geom_mf,
            time,
            0 // level
        );
    } else {
        WriteSingleLevelPlotfile(
            plotfilename,
            mf,
            varnames,
            geom_mf,
            time,
            0 // level
        );
    }
}