easy_dims.h

Go to the documentation of this file.

#pragma once
 
#include <netcdf.h>
#include "thrust/host_vector.h"
 
#include "dg/runge_kutta.h" // for dg::is_same
#include "dg/topology/grid.h"
#include "dg/topology/gridX.h"
#include "dg/topology/evaluation.h"
#ifdef MPI_VERSION
#include "dg/topology/mpi_grid.h"
#include "dg/topology/mpi_evaluation.h"
#endif //MPI_VERSION
 
#include "nc_error.h"
#include "easy_atts.h"
 
namespace dg
{
namespace file
{
 
 
namespace detail
{
 
template<class real_type>
std::vector<dg::RealGrid1d<real_type> > grids( const RealGrid0d<real_type>& ) { return {};}
template<class real_type>
std::vector<dg::RealGrid1d<real_type> > grids( const RealGrid1d<real_type>& g ) { return {g};}
template<class real_type>
std::vector<dg::RealGrid1d<real_type> > grids( const aRealTopology2d<real_type>& g ) { return {g.gy(), g.gx()};}
template<class real_type>
std::vector<dg::RealGrid1d<real_type> > grids( const aRealTopology3d<real_type>& g ) { return {g.gz(), g.gy(), g.gx()};}
template<class real_type>
std::vector<dg::RealGrid1d<real_type> > grids( const aRealTopologyX2d<real_type>& g ) { return grids(g.grid());}
template<class real_type>
std::vector<dg::RealGrid1d<real_type> > grids( const aRealTopologyX3d<real_type>& g ) { return grids(g.grid());;}
 
template<class real_type>
std::vector<std::string> dim_names( const RealGrid0d<real_type>& ) { return {};}
template<class real_type>
std::vector<std::string> dim_names( const RealGrid1d<real_type>& g ) { return {"x"};}
template<class real_type>
std::vector<std::string> dim_names( const aRealTopology2d<real_type>& g ) { return {"y", "x"};}
template<class real_type>
std::vector<std::string> dim_names( const aRealTopology3d<real_type>& g ) { return {"z", "y", "x"};}
template<class real_type>
std::vector<std::string> dim_names( const aRealTopologyX2d<real_type>& g ) { return {"y", "x"};}
template<class real_type>
std::vector<std::string> dim_names( const aRealTopologyX3d<real_type>& g ) { return {"z", "y", "x"};}
template<class real_type>
std::vector<std::string> axis_names( const RealGrid0d<real_type>& ) { return {};}
template<class real_type>
std::vector<std::string> axis_names( const RealGrid1d<real_type>& g ) { return {"X"};}
template<class real_type>
std::vector<std::string> axis_names( const aRealTopology2d<real_type>& g ) { return {"Y", "X"};}
template<class real_type>
std::vector<std::string> axis_names( const aRealTopology3d<real_type>& g ) { return {"Z", "Y", "X"};}
template<class real_type>
std::vector<std::string> axis_names( const aRealTopologyX2d<real_type>& g ) { return {"Y", "X"};}
template<class real_type>
std::vector<std::string> axis_names( const aRealTopologyX3d<real_type>& g ) { return {"Z", "Y", "X"};}
 
 
inline void assign_defaults( std::vector<std::string>& name_dims, const std::vector<std::string>& default_names)
{
    if( name_dims.empty())
        name_dims = default_names;
    if( name_dims.size() != default_names.size())
        throw std::runtime_error( "Number of given dimension names "+std::to_string(name_dims.size())+"does not match required number "+std::to_string(default_names.size())+"\n");
}
} // namespace detail
 
 
template<class T>
bool check_real_time( int ncid, const char* name, int* dimID, int* tvarID)
{
    // TODO Axis attribute check is still missing
    file::NC_Error_Handle err;
    // Check that it exists
    int retval = nc_inq_dimid( ncid, name, dimID);
    if( retval)
        return false;
    // Check that it is unlimited
    int num_unlim;
    err = nc_inq_unlimdims( ncid, &num_unlim, NULL);
    std::vector<int> unlim_dims(num_unlim);
    err = nc_inq_unlimdims( ncid, &num_unlim, &unlim_dims[0]);
    if( std::find( unlim_dims.begin(), unlim_dims.end(), *dimID) == std::end( unlim_dims))
        throw std::runtime_error( "Dimension "+std::string(name)+" already defined but not unlimited!\n");
    //size_t length;
    //err = nc_inq_dimlen( ncid, *dimID, &length);
    //if( length != 0)
    //    throw std::runtime_error( "Unlimited dimension "+std::string(name)+" already has values!\n");
    // Now check if the dimension variable exists already
    err = nc_inq_varid( ncid, name, tvarID);
    int xtype;
    int ndims;
    int dimids[10];
    err = nc_inq_var( ncid, *tvarID, NULL, &xtype, &ndims, dimids, NULL);
    if( xtype != detail::getNCDataType<T>() || ndims != 1 || dimids[0] != *dimID)
        throw std::runtime_error( "Unlimited variable "+std::string(name)+" has wrong type or wrong dimensions!\n");
    return true;
}
 
template<class T>
inline int define_real_time( int ncid, const char* name, int* dimID, int* tvarID, bool full_check = false)
{
    if( full_check)
        if( check_real_time<T>( ncid, name, dimID, tvarID))
            return NC_NOERR;
    int retval;
    if( (retval = nc_def_dim( ncid, name, NC_UNLIMITED, dimID)) ){ return retval;}
    if( (retval = nc_def_var( ncid, name, detail::getNCDataType<T>(), 1, dimID, tvarID))){return retval;}
    std::string t = "time since start"; //needed for paraview to recognize timeaxis
    // Update: Actually paraview also recognizes time from the "T" "axis" without "unit"
    std::string axis = "T";
    if( (retval = nc_put_att_text(ncid, *tvarID, "axis", axis.size(), axis.data())) ){return retval;}
    if( (retval = nc_put_att_text(ncid, *tvarID, "units", t.size(), t.data())) ){ return retval;}
    return retval;
}
 
inline int define_time( int ncid, const char* name, int* dimID, int* tvarID)
{
    return define_real_time<double>( ncid, name, dimID, tvarID);
}
 
 
inline int define_limited_time( int ncid, const char* name, int size, int* dimID, int* tvarID)
{
    int retval;
    if( (retval = nc_def_dim( ncid, name, size, dimID)) ){ return retval;}
    if( (retval = nc_def_var( ncid, name, NC_DOUBLE, 1, dimID, tvarID))){return retval;}
    std::string t = "time since start"; //needed for paraview to recognize timeaxis
    std::string axis = "T";
    if( (retval = nc_put_att_text(ncid, *tvarID, "axis", axis.size(), axis.data())) ){return retval;}
    if( (retval = nc_put_att_text(ncid, *tvarID, "units", t.size(), t.data())) ){ return retval;}
    return retval;
}
 
template<class T>
bool check_dimension( int ncid, int* dimID, const dg::RealGrid1d<T>& g, std::string name_dim = "x", std::string axis = "X")
{
    // TODO Axis attribute check is still missing
    file::NC_Error_Handle err;
    // check if the dimension exists already:
    int retval = nc_inq_dimid( ncid, name_dim.data(), dimID);
    if( retval)
        return false;
    size_t length;
    retval = nc_inq_dimlen( ncid, *dimID, &length);
    thrust::host_vector<T> points = g.abscissas();
    if( length != points.size())
        throw std::runtime_error( "Length of dimension "+name_dim+" does not match grid!\n");
    // Now check if the dimension variable exists already
    int varID;
    err = nc_inq_varid( ncid, name_dim.data(), &varID);
    int xtype;
    int ndims;
    int dimids[10];
    err = nc_inq_var( ncid, varID, NULL, &xtype, &ndims, dimids, NULL);
    if( xtype != detail::getNCDataType<T>() || ndims != 1 || dimids[0] != *dimID)
        throw std::runtime_error( "Dimension variable "+name_dim+" has wrong type or wrong dimensions!\n");
    thrust::host_vector<T> data( points);
    err = nc_get_var( ncid, varID, data.data());
    for( unsigned i=0; i<data.size(); i++)
        if( !dg::is_same( data[i], points[i]))
            throw std::runtime_error( "Dimension variable "+name_dim+" has values different from grid!\n");
    return true;
}
 
template<class T>
inline int define_dimension( int ncid, int* dimID, const dg::RealGrid1d<T>& g, std::string name_dim = "x", std::string axis = "X", bool full_check = false)
{
    if( full_check)
        if( check_dimension( ncid, dimID, g, name_dim, axis))
            return NC_NOERR;
    int retval;
    std::string long_name = name_dim+"-coordinate in Computational coordinate system";
    thrust::host_vector<T> points = g.abscissas();
    if( (retval = nc_def_dim( ncid, name_dim.data(), points.size(), dimID))){ return retval;}
    int varID;
    if( (retval = nc_def_var( ncid, name_dim.data(), detail::getNCDataType<T>(), 1, dimID, &varID))){return retval;}
    if( (retval = nc_put_var( ncid, varID, points.data())) ){ return retval;}
    if( (retval = nc_put_att_text( ncid, varID, "axis", axis.size(), axis.data())) ){return retval;}
    retval = nc_put_att_text( ncid, varID, "long_name", long_name.size(), long_name.data());
    return retval;
}
 
template<class Topology, std::enable_if_t<dg::is_vector_v<typename Topology::host_vector, dg::SharedVectorTag>,bool > = true>
int define_dimensions( int ncid, int *dimsIDs, const Topology& g, std::vector<std::string> name_dims = {}, bool full_check = false)
{
    int retval = NC_NOERR;
    auto grids = detail::grids( g);
    auto default_names = detail::dim_names(g);
    auto axis_names = detail::axis_names(g);
    detail::assign_defaults( name_dims, default_names);
    for ( unsigned i=0; i<grids.size(); i++)
    {
        retval = define_dimension( ncid, &dimsIDs[i], grids[i], name_dims[i], axis_names[i], full_check);
        if( retval)
            return retval;
    }
    return retval;
}
 
template<class Topology, std::enable_if_t<dg::is_vector_v<typename Topology::host_vector, dg::SharedVectorTag>,bool > = true>
int define_dimensions( int ncid, int *dimsIDs, int* tvarID, const Topology& g, std::vector<std::string> name_dims = {}, bool full_check = false)
{
    auto default_names = detail::dim_names(g);
    default_names.insert( default_names.begin(), "time");
    detail::assign_defaults( name_dims, default_names);
 
    int retval = define_real_time<typename Topology::value_type>( ncid, name_dims[0].data(), &dimsIDs[0], tvarID, full_check);
    if(retval)
        return retval;
    return define_dimensions( ncid, &dimsIDs[1], g, {name_dims.begin()+1, name_dims.end()}, full_check);
}
 
template<class Topology, std::enable_if_t<dg::is_vector_v<typename Topology::host_vector, dg::SharedVectorTag>,bool > = true>
bool check_dimensions( int ncid, int *dimsIDs, const Topology& g, std::vector<std::string> name_dims = {})
{
    auto grids = detail::grids( g);
    auto default_names = detail::dim_names(g);
    auto axis_names = detail::axis_names(g);
    detail::assign_defaults( name_dims, default_names);
    for ( unsigned i=0; i<grids.size(); i++)
    {
        if (!check_dimension( ncid, &dimsIDs[i], grids[i], name_dims[i], axis_names[i]))
            return false;
    }
    return true;
}
 
template<class Topology, std::enable_if_t<dg::is_vector_v<typename Topology::host_vector, dg::SharedVectorTag>,bool > = true>
bool check_dimensions( int ncid, int *dimsIDs, int* tvarID, const Topology& g, std::vector<std::string> name_dims = {})
{
    auto default_names = detail::dim_names(g);
    default_names.insert( default_names.begin(), "time");
    detail::assign_defaults( name_dims, default_names);
    if ( !check_real_time<typename Topology::value_type>( ncid, name_dims[0].data(), &dimsIDs[0], tvarID))
        return false;
    return check_dimensions( ncid, &dimsIDs[1], g, {name_dims.begin()+1, name_dims.end()});
}
 
 
 
 
template<class T>
inline int define_limtime_xy( int ncid, int* dimsIDs, int size, int* tvarID, const dg::aRealTopology2d<T>& g, std::vector<std::string> name_dims = {"time", "y", "x"})
{
    int retval;
    retval = define_limited_time( ncid, name_dims[0].data(), size, &dimsIDs[0], tvarID);
    if(retval)
        return retval;
    return define_dimensions( ncid, &dimsIDs[1], g, {name_dims[1], name_dims[2]});
}
 
 
#ifdef MPI_VERSION
template<class MPITopology, std::enable_if_t<dg::is_vector_v<typename MPITopology::host_vector, dg::MPIVectorTag>,bool > = true>
inline int define_dimensions( int ncid, int* dimsIDs, const MPITopology& g, std::vector<std::string> name_dims = {}, bool full_check = false)
{
    int rank;
    MPI_Comm_rank( MPI_COMM_WORLD, &rank);
    if( rank==0) return define_dimensions( ncid, dimsIDs, g.global(), name_dims, full_check);
    else
        return NC_NOERR;
}
template<class MPITopology, std::enable_if_t<dg::is_vector_v<typename MPITopology::host_vector, dg::MPIVectorTag>,bool > = true>
inline int define_dimensions( int ncid, int* dimsIDs, int* tvarID, const MPITopology& g, std::vector<std::string> name_dims = {}, bool full_check = false)
{
    int rank;
    MPI_Comm_rank( MPI_COMM_WORLD, &rank);
    if( rank==0) return define_dimensions( ncid, dimsIDs, tvarID, g.global(), name_dims, full_check);
    else
        return NC_NOERR;
}
template<class MPITopology, std::enable_if_t<dg::is_vector_v<typename MPITopology::host_vector, dg::MPIVectorTag>,bool > = true>
inline bool check_dimensions( int ncid, int* dimsIDs, const MPITopology& g, std::vector<std::string> name_dims = {})
{
    // all processes can read NetCDF in parallel by default
    return check_dimensions( ncid, dimsIDs, g.global(), name_dims);
}
template<class MPITopology, std::enable_if_t<dg::is_vector_v<typename MPITopology::host_vector, dg::MPIVectorTag>,bool > = true>
inline bool check_dimensions( int ncid, int* dimsIDs, int* tvarID, const MPITopology& g, std::vector<std::string> name_dims = {})
{
    // all processes can read NetCDF in parallel by default
    return check_dimensions( ncid, dimsIDs, tvarID, g.global(), name_dims);
}
#endif //MPI_VERSION
 
} //namespace file
} //namespace dg