probes.h

Go to the documentation of this file.

#pragma once
 
 
#include "dg/topology/interpolation.h"
#ifdef MPI_VERSION
#include "dg/topology/mpi_projection.h"
#include "nc_mpi_file.h"
#endif //MPI_VERSION
 
#include "probes_params.h"
#include "nc_file.h"
#include "records.h"
 
namespace dg
{
namespace file
{
template<class NcFile, class Topology>
struct Probes
{
    Probes() = default;
    Probes(
        NcFile& file,
        const Topology& grid,
        const ProbesParams& params // do nothing if probes is false
        ) : m_file(&file)
    {
        m_probes = params.probes;
        if( !params.probes) return;
 
        if ( params.coords.size() != grid.ndim())
            throw std::runtime_error( "Need "+std::to_string(grid.ndim())+" values in coords!");
        unsigned num_pins = params.get_coords_sizes();
        // params.coords is empty on ranks other than master
 
// TODO We could think about distributing the coords among ranks using grid.contains ...
        m_probe_interpolate = dg::create::interpolation( params.coords, grid,
            grid.get_bc(), "dg");
        // Create helper storage probe variables
        m_simple_probes = create_probes_vec( params.coords[0], grid);
        m_resultH = dg::evaluate( dg::zero, grid);
 
        file.def_grp( "probes");
        m_grp = file.get_current_path() / "probes";
        file.set_grp( m_grp);
        file.put_att( {"format", params.format});
        file.def_dim( "pdim", num_pins);
        file.template def_dimvar_as<double>( "ptime", NC_UNLIMITED, {{"axis" , "T"}});
        for( unsigned i=0; i<params.coords.size(); i++)
        {
            auto probes_vec = create_probes_vec( params.coords[0], grid);
            file.defput_var( params.coords_names[i], {"pdim"},
                {{"long_name" , "Coordinate variable for probe position"}},
                {probes_vec}, probes_vec);
        }
        file.set_grp("..");
    }
 
    template<class ListClass, class ...Params>
    void static_write( const ListClass& records, Params&& ... ps)
    {
        if(!m_probes) return;
        using Result = get_first_argument_type_t<typename ListClass::value_type::Signature>;
        auto result = dg::construct<Result>( m_resultH);
        static_assert( dg::is_vector_v<Result>, "Result must be a vector type");
        m_file->set_grp( m_grp);
 
        for ( auto& record : records)
        {
            record.function( result, std::forward<Params>(ps)...);
            dg::assign( result, m_resultH);
            dg::blas2::symv( m_probe_interpolate, m_resultH, m_simple_probes);
 
            m_file->defput_var( record.name, {"pdim"}, record.atts,
                {m_simple_probes}, m_simple_probes);
        }
        m_file->set_grp("..");
    }
 
    template<class ListClass, class ...Params>
    void buffer( double time, const ListClass& probe_list, Params&& ... ps)
    {
        if(!m_probes) return;
        using Result = get_first_argument_type_t<typename ListClass::value_type::Signature>;
        auto result = dg::construct<Result>( m_resultH);
        m_time_intern.push_back(time);
        if( m_simple_probes_intern.empty())
            init<dg::get_value_type<Result>>( probe_list);
 
        for( auto& record : probe_list)
        {
            record.function( result, std::forward<Params>(ps)...);
            dg::assign( result, m_resultH);
            dg::blas2::symv( m_probe_interpolate, m_resultH, m_simple_probes);
            m_simple_probes_intern.at(record.name).push_back(m_simple_probes);
        }
    }
 
    void flush()
    {
        if(!m_probes) return;
        // else we write the internal buffer to file
        m_file->set_grp( m_grp);
        for( unsigned j=0; j<m_time_intern.size(); j++)
        {
            m_file->put_var( "ptime", {m_probe_start}, m_time_intern[j]);
            for( auto& field : m_simple_probes_intern)
                m_file->put_var( field.first, {m_probe_start, field.second[j]},
                    field.second[j]);
            m_probe_start ++;
        }
        // clear internal buffers
        m_time_intern.clear();
        for( auto& field : m_simple_probes_intern)
            field.second.clear();
        m_file->set_grp( "..");
    }
 
    template<class ListClass, class ...Params>
    void write( double time, const ListClass& probe_list, Params&& ... ps)
    {
        buffer( time, probe_list, std::forward<Params>(ps)...);
        flush();
    }
 
    private:
    template<class T=Topology, std::enable_if_t<dg::is_vector_v<typename
    T::host_vector, dg::SharedVectorTag>, bool> = true>
    auto create_probes_vec( const dg::HVec& coord, const T& grid)
    {
        return typename Topology::host_vector(coord);
    }
    template<class T=Topology, std::enable_if_t<dg::is_vector_v<typename
    T::host_vector, dg::MPIVectorTag>, bool> = true>
    auto create_probes_vec( const dg::HVec& coord, const T& grid)
    {
        return typename Topology::host_vector(coord, grid.communicator());
    }
    template<class value_type, class ListClass>
    void init( const ListClass& probe_list)
    {
        m_file->set_grp( m_grp);
        for( auto& record : probe_list)
        {
            m_simple_probes_intern[record.name] = {}; // empty vectors
            m_file->template def_var_as<value_type>( record.name, {"ptime", "pdim"});
            m_file->put_atts( record.atts, record.name);
        }
        m_file->set_grp( "..");
    }
    bool m_probes = false;
    NcFile* m_file;
    std::filesystem::path m_grp;
    int m_probe_grp_id;
    typename Topology::host_vector m_simple_probes;
#ifdef MPI_VERSION
    std::conditional_t<
        dg::is_vector_v<typename Topology::host_vector, dg::SharedVectorTag>,
        dg::IHMatrix_t< typename Topology::value_type>,
        dg::MIHMatrix_t<typename Topology::value_type> > m_probe_interpolate;
#else
    dg::IHMatrix_t< typename Topology::value_type> m_probe_interpolate;
#endif
    std::map<std::string, std::vector<typename Topology::host_vector>> m_simple_probes_intern;
    std::vector<double> m_time_intern;
    //Container m_resultD;
    typename Topology::host_vector m_resultH;
    size_t m_probe_start = 0; // always point to where we currently can write
};
 
} //namespace file
}//namespace dg