ERF_NCWpsFile.H File Reference

#include <sstream>
#include <string>
#include <atomic>
#include <cfenv>
#include "AMReX_FArrayBox.H"
#include "AMReX_IArrayBox.H"
#include "ERF_EpochTime.H"
#include "ERF_NCInterface.H"

Include dependency graph for ERF_NCWpsFile.H:

This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Classes
struct	NDArray< DataType >

Typedefs
using	PlaneVector = amrex::Vector< amrex::FArrayBox >

Enumerations
enum class	NC_Data_Dims_Type {   Time_SL_SN_WE , Time_BT_SN_WE , Time_SN_WE , Time_BT ,   Time_SL , Time , Time_BdyWidth_BT_SN , Time_BdyWidth_BT_WE ,   Time_BdyWidth_SN , Time_BdyWidth_WE }

Functions
int	BuildFABsFromWRFBdyFile (const std::string &fname, amrex::Vector< amrex::Vector< amrex::FArrayBox >> &bdy_data_xlo, amrex::Vector< amrex::Vector< amrex::FArrayBox >> &bdy_data_xhi, amrex::Vector< amrex::Vector< amrex::FArrayBox >> &bdy_data_ylo, amrex::Vector< amrex::Vector< amrex::FArrayBox >> &bdy_data_yhi)
template<typename DType >
void	ReadTimeSliceFromNetCDFFile (const std::string &fname, const int tidx, amrex::Vector< std::string > names, amrex::Vector< NDArray< DType > > &arrays, amrex::Vector< int > &success)
template<typename DType >
void	ReadNetCDFFile (const std::string &fname, amrex::Vector< std::string > names, amrex::Vector< NDArray< DType > > &arrays, amrex::Vector< int > &success)
template<class FAB , typename DType >
void	fill_fab_from_arrays (int iv, const amrex::Box &domain, amrex::Vector< NDArray< float >> &nc_arrays, const std::string &var_name, NC_Data_Dims_Type &NC_dim_type, FAB &temp)
template<class FAB , typename DType >
void	BuildFABsFromNetCDFFile (const amrex::Box &domain, const std::string &fname, amrex::Vector< std::string > nc_var_names, amrex::Vector< enum NC_Data_Dims_Type > NC_dim_types, amrex::Vector< FAB * > fab_vars, amrex::Vector< int > &success, bool keep_on_host=false)

Typedef Documentation

PlaneVector

using PlaneVector = amrex::Vector<amrex::FArrayBox>

Enumeration Type Documentation

NC_Data_Dims_Type

enum NC_Data_Dims_Type

strong

Enumerator
Time_SL_SN_WE
Time_BT_SN_WE
Time_SN_WE
Time_BT
Time_SL
Time
Time_BdyWidth_BT_SN
Time_BdyWidth_BT_WE
Time_BdyWidth_SN
Time_BdyWidth_WE

                             {
    Time_SL_SN_WE,       // Time, Soil Layers, South-North, West-East
    Time_BT_SN_WE,
    Time_SN_WE,
    Time_BT,
    Time_SL,             // Time, Soil layers
    Time,
    Time_BdyWidth_BT_SN,
    Time_BdyWidth_BT_WE,
    Time_BdyWidth_SN,
    Time_BdyWidth_WE
};

Function Documentation

BuildFABsFromNetCDFFile()

template<class FAB , typename DType >

void BuildFABsFromNetCDFFile	(	const amrex::Box &	domain,
		const std::string &	fname,
		amrex::Vector< std::string >	nc_var_names,
		amrex::Vector< enum NC_Data_Dims_Type >	NC_dim_types,
		amrex::Vector< FAB * >	fab_vars,
		amrex::Vector< int > &	success,
		bool	keep_on_host = false
	)

Function to read NetCDF variables and fill the corresponding Array4's

Parameters

domain	Box specifying the domain covered by the NetCDF fields
fname	Name of the NetCDF file to be read
nc_var_names	Variable names in the NetCDF file
NC_dim_types	NetCDF data dimension types
fab_vars	Fab data we are to fill
success	Success flags for each requested variable
keep_on_host	If true, leave each fab in the (pinned, on GPU) host memory it was broadcast into rather than copying it to the device. Each fab spans the whole domain on every rank, so a caller that only uses it as the source of copies into its own MultiFabs should set this to avoid a full-domain device allocation per variable.

{
    int ioproc = amrex::ParallelDescriptor::IOProcessorNumber();  // I/O rank
 
    amrex::Vector<NDArray<float>> nc_arrays(nc_var_names.size());
 
    if (amrex::ParallelDescriptor::IOProcessor())
    {
        ReadNetCDFFile(fname, nc_var_names, nc_arrays, success);
    }
 
    amrex::ParallelDescriptor::Bcast(success.dataPtr(), success.size(), ioproc);
 
    for (int iv = 0; iv < nc_var_names.size(); iv++)
    {
        if (success[iv] == 1) {
            FAB tmp;
            if (amrex::ParallelDescriptor::IOProcessor()) {
                fill_fab_from_arrays<FAB,DType>(iv, domain, nc_arrays, nc_var_names[iv], NC_dim_types[iv], tmp);
            }
 
            int ncomp = tmp.nComp();
            amrex::Box box = tmp.box();
 
            amrex::ParallelDescriptor::Bcast(&box,   1, ioproc);
            amrex::ParallelDescriptor::Bcast(&ncomp, 1, ioproc);
 
            if (!amrex::ParallelDescriptor::IOProcessor()) {
#ifdef AMREX_USE_GPU
                tmp.resize(box,ncomp,amrex::The_Pinned_Arena());
#else
                tmp.resize(box,ncomp);
#endif
            }
 
            amrex::ParallelDescriptor::Bcast(tmp.dataPtr(), tmp.size(), ioproc);
 
            // Shift box by the domain lower corner
            amrex::Box  fab_bx = tmp.box();
            amrex::Dim3 dom_lb = lbound(domain);
            fab_bx += amrex::IntVect(dom_lb.x,dom_lb.y,dom_lb.z);
            if (keep_on_host) {
                // Hand over the host buffer rather than copying it
                *fab_vars[iv] = std::move(tmp);
                fab_vars[iv]->shift(amrex::IntVect(dom_lb.x,dom_lb.y,dom_lb.z));
                continue;
            }
            // fab_vars points to data on device
            fab_vars[iv]->resize(fab_bx,1);
#ifdef AMREX_USE_GPU
            amrex::Gpu::copy(amrex::Gpu::hostToDevice,
                             tmp.dataPtr(), tmp.dataPtr() + tmp.size(),
                             fab_vars[iv]->dataPtr());
#else
            // Provided by BaseFab inheritance through FArrayBox
            fab_vars[iv]->copy(tmp,tmp.box(),0,fab_bx,0,1);
#endif
        } // success
    } // iv
}

Here is the call graph for this function:

BuildFABsFromWRFBdyFile()

int BuildFABsFromWRFBdyFile	(	const std::string &	fname,
		amrex::Vector< amrex::Vector< amrex::FArrayBox >> &	bdy_data_xlo,
		amrex::Vector< amrex::Vector< amrex::FArrayBox >> &	bdy_data_xhi,
		amrex::Vector< amrex::Vector< amrex::FArrayBox >> &	bdy_data_ylo,
		amrex::Vector< amrex::Vector< amrex::FArrayBox >> &	bdy_data_yhi
	)

fill_fab_from_arrays()

template<class FAB , typename DType >

void fill_fab_from_arrays	(	int	iv,
		const amrex::Box &	domain,
		amrex::Vector< NDArray< float >> &	nc_arrays,
		const std::string &	var_name,
		NC_Data_Dims_Type &	NC_dim_type,
		FAB &	temp
	)

Helper function for reading data from NetCDF file into a provided FAB.

Parameters

iv	Index for which variable we are going to fill
domain	Box specifying the horizontal and vertical extent to read
nc_arrays	Arrays of data from NetCDF file
var_name	Variable name
NC_dim_type	Dimension type for the variable as stored in the NetCDF file
temp	FAB where we store the variable data from the NetCDF Arrays

{
    // bottom_top, south_north, west_east (these can be staggered or unstaggered).
    // Initialised because amrex::Abort is not [[noreturn]], so the compiler
    // cannot tell that the final else never falls through.
    int ns1 = 0, ns2 = 0, ns3 = 0;
    if (NC_dim_type == NC_Data_Dims_Type::Time_BT) {
        ns1 = nc_arrays[iv].get_vshape()[1];
        ns2 = 1;
        ns3 = 1;
        // amrex::Print() << "TYPE BT " << ns1 << std::endl;
    } else if (NC_dim_type == NC_Data_Dims_Type::Time_SN_WE) {
        ns1 = 1;
        ns2 = nc_arrays[iv].get_vshape()[1];
        ns3 = nc_arrays[iv].get_vshape()[2];
        // amrex::Print() << "TYPE SN WE " << ns2 << " " << ns3 << std::endl;
    } else if (NC_dim_type == NC_Data_Dims_Type::Time_BT_SN_WE) {
        ns1 = nc_arrays[iv].get_vshape()[1];
        ns2 = nc_arrays[iv].get_vshape()[2];
        ns3 = nc_arrays[iv].get_vshape()[3];
        // amrex::Print() << "TYPE BT SN WE " << ns1 << " " << ns2 << " " << ns3 << std::endl;
    } else if (NC_dim_type == NC_Data_Dims_Type::Time_SL_SN_WE) {
        ns1 = nc_arrays[iv].get_vshape()[1];
        ns2 = nc_arrays[iv].get_vshape()[2];
        ns3 = nc_arrays[iv].get_vshape()[3];
        // amrex::Print() << "TYPE SL SN WE " << ns1 << " " << ns2 << " " << ns3 << std::endl;
    } else if (NC_dim_type == NC_Data_Dims_Type::Time_SL) {
        ns1 = nc_arrays[iv].get_vshape()[1];
        ns2 = 1;
        ns3 = 1;
        // amrex::Print() << "TYPE SL " << ns1 << std::endl;
    }
    else {
        amrex::Abort("fill_fab_from_arrays: variable " + var_name +
                     " has NC_Data_Dims_Type " + std::to_string(static_cast<int>(NC_dim_type)) +
                     ", which this reader does not handle (only Time_BT, Time_SN_WE,"
                     " Time_BT_SN_WE, Time_SL_SN_WE and Time_SL)");
    }
 
    // allow fewer z levels to be read in
    // - ns1 may be the bottom_top or bottom_top_stag dim
    // - domain.bigEnd(2) is the input number of cells in z
    // - khi may be nz==bottom_top, nz < bottom_top, or 0 for 2D fields
    int khi = std::min(domain.bigEnd(2), ns1-1);
 
    // TODO:  The box will only start at (0,0,0) at level 0 -- we need to generalize this
    amrex::Box in_box(amrex::IntVect(0,0,0), amrex::IntVect(ns3-1,ns2-1,ns1-1)); // ns1 may or may not be staggered
    amrex::Box my_box(amrex::IntVect(0,0,0), amrex::IntVect(ns3-1,ns2-1,khi));
 
    if (var_name == "PH" || var_name == "PHB") {
        my_box.setType(amrex::IndexType(amrex::IntVect(0,0,1)));
        my_box.setBig(2, khi+1);
    }
    else if (var_name == "U" || var_name == "UU" || var_name == "MAPFAC_U") {
        my_box.setType(amrex::IndexType(amrex::IntVect(1,0,0)));
    }
    else if (var_name == "V" || var_name == "VV" || var_name == "MAPFAC_V") {
        my_box.setType(amrex::IndexType(amrex::IntVect(0,1,0)));
    }
    else if (var_name == "W" || var_name == "WW") {
        my_box.setType(amrex::IndexType(amrex::IntVect(0,0,1)));
        my_box.setBig(2, khi+1);
    }
    else if (var_name == "TSLB" || var_name == "SMOIS" || var_name == "SH2O" || var_name == "ZS" || var_name == "DZS") {
        // soil variables are staggered in z
        my_box.setType(amrex::IndexType(amrex::IntVect(0,0,1)));
    }
 
    amrex::Arena* Arena_Used = amrex::The_Arena();
#ifdef AMREX_USE_GPU
    // Make sure temp lives on CPU since nc_arrays lives on CPU only
    Arena_Used = amrex::The_Pinned_Arena();
#endif
    temp.resize(my_box,1, Arena_Used);
    amrex::Array4<DType> fab_arr = temp.array();
 
    int ioff = temp.box().smallEnd()[0];
    int joff = temp.box().smallEnd()[1];
    int kmax = temp.box().bigEnd()[2]; // highest z level to be read in
    //amrex::Print() << var_name << " domain khi, kmax, ns1 = " << khi << " " << kmax << " " << ns1 << std::endl;
 
    auto num_pts = in_box.numPts();
 
    for (int n(0); n < num_pts; ++n) {
        int k  = n / (ns2*ns3);
        if (k > kmax) continue;
        int j  = (n - k*(ns2*ns3)) / ns3 + joff;
        int i  = n - k*(ns2*ns3) - (j-joff) * ns3 + ioff;
        fab_arr(i,j,k,0) = static_cast<DType>(*(nc_arrays[iv].get_data()+n));
    }
}

ReadNetCDFFile()

template<typename DType >

void ReadNetCDFFile	(	const std::string &	fname,
		amrex::Vector< std::string >	names,
		amrex::Vector< NDArray< DType > > &	arrays,
		amrex::Vector< int > &	success
	)

{
    AMREX_ASSERT(arrays.size() == names.size());
 
    if (amrex::ParallelDescriptor::IOProcessor())
    {
        auto prev_fpe_except = amrex::disableFPExcept(amrex::FPExcept::all);
 
        auto ncf = ncutils::NCFile::open(fname, NC_CLOBBER | NC_NETCDF4);
 
        /*
        // get the dimension information
        int Time               = static_cast<int>(ncf.dim("Time").len());
        int DateStrLen         = static_cast<int>(ncf.dim("DateStrLen").len());
        int west_east          = static_cast<int>(ncf.dim("west_east").len());
        int south_north        = static_cast<int>(ncf.dim("south_north").len());
        int bottom_top         = static_cast<int>(ncf.dim("bottom_top").len());
        int bottom_top_stag    = static_cast<int>(ncf.dim("bottom_top_stag").len());
        int west_east_stag     = static_cast<int>(ncf.dim("west_east_stag").len());
        int south_north_stag   = static_cast<int>(ncf.dim("south_north_stag").len());
        int bdy_width          = static_cast<int>(ncf.dim("bdy_width").len());
        */
 
        // amrex::Print() << "Reading the dimensions from the netcdf file " << "\n";
 
        for (auto n=0; n<arrays.size(); ++n)
        {
            std::string vname_to_write = names[n];
            std::string vname_to_read  = names[n];
 
            if (vname_to_read.substr(0,2) == "R_") {
                vname_to_read  = names[n+4]; // This allows us to read "T" instead -- we will over-write this later
            }
 
            success[n] = ncf.has_var(vname_to_read);
            if (success[n] == 0) {
                amrex::Print() << " Skipping " << vname_to_read << std::endl;
            } else {
                amrex::Print() << " Reading  " << vname_to_read << std::endl;
            }
 
            if (success[n] == 1) {
 
                std::vector<std::string> dimnames = ncf.var(vname_to_read).dimnames();
                AMREX_ALWAYS_ASSERT(dimnames[0] == "Time" || dimnames[0] == "time");
 
                std::vector<size_t> shape = ncf.var(vname_to_read).shape();
                arrays[n]                 = NDArray<DType>(vname_to_read,shape);
                DType* dataPtr            = arrays[n].get_data();
 
                std::vector<size_t> start(shape.size(), 0);
 
#if 0
                auto numPts               = arrays[n].ndim();
                amrex::Print() << "NetCDF Variable name = " << vname_to_read << std::endl;
                amrex::Print() << "numPts read from NetCDF file/var = " << numPts << std::endl;
                amrex::Print() << "Dims in var = " << ncf.var(vname_to_read).ndim() << std::endl;
                amrex::Print()  << "Dim names = (";
                for (auto &dim:dimnames)
                    amrex::Print() << dim << ", " ;
                amrex::Print()  << ")" << std::endl;
                amrex::Print()  << "Dims of the variable = (";
                for (auto &dim:shape)
                    amrex::Print() << dim << ", " ;
                amrex::Print()  << ")" << std::endl;
#endif
 
                ncf.var(vname_to_read).get(dataPtr, start, shape);
 
            } // has_var
        }
        ncf.close();
 
        std::feclearexcept(FE_ALL_EXCEPT);
        amrex::setFPExcept(prev_fpe_except);
    }
}

Referenced by BuildFABsFromNetCDFFile().

Here is the call graph for this function:

Here is the caller graph for this function:

ReadTimeSliceFromNetCDFFile()

template<typename DType >

void ReadTimeSliceFromNetCDFFile	(	const std::string &	fname,
		const int	tidx,
		amrex::Vector< std::string >	names,
		amrex::Vector< NDArray< DType > > &	arrays,
		amrex::Vector< int > &	success
	)

{
    amrex::Print() << "Reading time slice " << tidx << " from " << fname << std::endl;
    AMREX_ASSERT(arrays.size() == names.size());
 
    if (amrex::ParallelDescriptor::IOProcessor())
    {
        auto ncf = ncutils::NCFile::open(fname, NC_CLOBBER | NC_NETCDF4);
 
        int ntimes = ncf.dim("Time").len();
        AMREX_ALWAYS_ASSERT((tidx >= 0) && (tidx < ntimes));
 
        for (auto n=0; n<arrays.size(); ++n)
        {
            std::string vname_to_write = names[n];
            std::string vname_to_read  = names[n];
            if (vname_to_read.substr(0,2) == "R_") {
                vname_to_read  = names[n+4]; // This allows us to read "T" instead -- we will over-write this later
            }
 
            success[n] = ncf.has_var(vname_to_read);
 
            if (success[n] == 1)
            {
                std::vector<std::string> dimnames = ncf.var(vname_to_read).dimnames();
                AMREX_ALWAYS_ASSERT(dimnames[0] == "Time");
 
                std::vector<size_t> count = ncf.var(vname_to_read).shape();
                std::vector<size_t> start(count.size(), 0);
                start[0] = tidx;
                count[0] = 1;
 
                arrays[n] = NDArray<DType>(vname_to_read, count);
                DType* dataPtr = arrays[n].get_data();
 
                ncf.var(vname_to_read).get(dataPtr, start, count);
            } // has_var
        }
        ncf.close();
    }
}

Here is the call graph for this function: