openPMD_io.hh

Go to the documentation of this file.

#ifndef RAYTRACE_API_HH
#define RAYTRACE_API_HH
#include "ray.hh"
#include <openPMD/openPMD.hpp> // openPMD C++ API
#include <string>

#include <exception>

namespace raytracing {
#define ITER 1

class openPMD_io {
    // Auxiliary classes that are private, so not part of the public API
private:
    class Rays {
        // public to make it available to openPMD_io methods
    public:
        template <typename T> class Record {
            std::vector<T> _vals;
            T _min, _max;

        public:
            Record(): _vals(), _min(), _max() { clear(); }
            const std::vector<T>& vals(void) const { return _vals; };
            std::vector<T>& vals(void) { return _vals; };
            T min(void) const { return _min; };
            T max(void) const { return _max; };

            // used when filling before writing
            void push_back(T val) {
                _vals.push_back(val);
                if (_min > val) (_min) = val;
                if (_max < val) (_max) = val;
            }

            // this is used when reading from the openPMD file
            void store(const T* vec, size_t n, float min, float max) {
                //              clear();
                _vals.insert(_vals.begin(), vec, vec + n);
                _min = min;
                _max = max;
            }

            void clear_chunk(void) { _vals.clear(); }

            void clear(void) {
                std::numeric_limits<T> lim;
                _max = lim.min();
                _min = lim.max();
                _vals.clear();
            }
            const T operator[](size_t i) const { return _vals[i]; }; // cannot modify
        }; // end of Record class

        //------------------------------ public memebers
    public:
        // I make the Records public members to avoid writing methods to access them
        //
        // the 3D components are in separate vectors because this is the way the openPMD API
        // wants them to be
        Record<float> _x, _y, _z,                   // position
                _dx, _dy, _dz,                      // direction (vx^2+vy^2+vz^2) = 1
                _sx, _sy, _sz,                      // non-photon polarization
                _sPolAx, _sPolAy, _sPolAz, _sPolPh, // photon s-polarization amplitude
                _pPolAx, _pPolAy, _pPolAz, _pPolPh, // photon p-polarization amplitude
                _wavelength,                        // wavelength
                _time, _weight;                     // ray time, weight
        Record<unsigned long long int> _id;         // id
        Record<int> _status;                        // alive status

        //------------------------------ private memebers
    private:
        size_t _size = 0; // number of stored rays = min(chunk_size, remaining rays to read)
        size_t _read = 0; // current index when reading

        //------------------------------ public methods
    public:
        Rays();

        void push(const Ray& this_ray);

        Ray pop(bool next = true);

        void clear_chunk(void) {
            _size = 0;
            _read = 0;

            _x.clear_chunk();
            _y.clear_chunk();
            _z.clear_chunk();

            _dx.clear_chunk();
            _dy.clear_chunk();
            _dz.clear_chunk();

            _sx.clear_chunk();
            _sy.clear_chunk();
            _sz.clear_chunk();

            _sPolAx.clear_chunk();
            _sPolAy.clear_chunk();
            _sPolAz.clear_chunk();
            _pPolAx.clear_chunk();
            _pPolAy.clear_chunk();
            _pPolAz.clear_chunk();

            _sPolPh.clear_chunk();
            _pPolPh.clear_chunk();

            _wavelength.clear_chunk();

            _time.clear_chunk();
            _weight.clear_chunk();

            _id.clear_chunk();
            _status.clear_chunk();
        }
        void clear(void) {
            _size = 0;
            _read = 0;

            _x.clear();
            _y.clear();
            _z.clear();

            _dx.clear();
            _dy.clear();
            _dz.clear();

            _sx.clear();
            _sy.clear();
            _sz.clear();

            _sPolAx.clear();
            _sPolAy.clear();
            _sPolAz.clear();
            _pPolAx.clear();
            _pPolAy.clear();
            _pPolAz.clear();

            _sPolPh.clear();
            _pPolPh.clear();

            _wavelength.clear();

            _time.clear();
            _weight.clear();

            _id.clear();
            _status.clear();
        };

        size_t size() const { return _size; };

        void size(size_t s) {
            _size = s;
            if (_x.vals().size() != s)
                throw std::runtime_error(
                        "size of stored Rays and one of its records is different");
        };

        bool is_chunk_finished(void) { return _read == _size; }
    }; // end of Rays class

public:
    explicit openPMD_io(
            const std::string& filename, 
            const std::string mc_code_name =
                    "", 
            const std::string mc_code_version = "", 
            const std::string instrument_name = "", 
            const std::string name_current_component =
                    "" 
    );

    /***************************************************************/

    void init_write(std::string particle_species,  
                    unsigned long long int n_rays, 

                    unsigned int iter = 1 
    );

    void
    init_rays(std::string particle_species, unsigned long long int n_rays, unsigned int iter);

    void trace_write(Ray this_ray);

    void save_write(void);

    /***************************************************************/

    unsigned long long int
    init_read(std::string pdgId,                 
              unsigned int iter             = 1, 
              unsigned long long int n_rays = 0, 
              unsigned int repeat =
                      1 
    );

    Ray trace_read(void);

    void set_gravity_direction(float x, float y, float z);

    void get_gravity_direction(float* x, float* y, float* z);
    void get_horizontal_direction(float* x, float* y, float* z);
    void set_horizontal_direction(float* x, float* y, float* z);
private:
    void load_chunk(void);

private:
    // parameters defined at construction
    std::string _name;
    std::string _mc_code_name;
    std::string _mc_code_version;
    std::string _instrument_name;
    std::string _name_current_component;
    unsigned int _i_repeat, _n_repeat;
    unsigned long long int _nrays, _max_allowed_rays;

    // internal usage
    //  openPMD::Access _access_mode;
    openPMD::Offset _offset;
    bool _isWriteMode;
    std::unique_ptr<openPMD::Series> _series;
    Rays _rays;
    Ray _last_ray;
    unsigned int _iter;
    std::string _particle_species;

    //------------------------------ set of helper methods
    inline openPMD::Iteration& iter_pmd(unsigned int iter) { return _series->iterations[iter]; }

    // returns the current particle species from the current iteration
    inline openPMD::ParticleSpecies& rays_pmd(void) {
        auto i = iter_pmd(_iter);
        return i.particles[_particle_species];
    }

    // returns the given particle species from the current iteration
    inline openPMD::ParticleSpecies& rays_pmd(std::string particle_species) {
        _particle_species = particle_species;
        auto i            = iter_pmd(_iter);
        return i.particles[particle_species];
    }

    void init_ray_prop(std::string name,          
                       openPMD::Dataset& dataset, 
                       bool isScalar,             
                       std::map<openPMD::UnitDimension, double> const& dims =
                               {{openPMD::UnitDimension::L, 0.}}, 
                       double unitSI = 0.);                       
};

} // namespace raytracing
#endif