Skip to content
Snippets Groups Projects
particles_system_builder.hpp 13.17 KiB
#ifndef PARTICLES_SYSTEM_BUILDER_HPP
#define PARTICLES_SYSTEM_BUILDER_HPP

#include <string>

#include "abstract_particles_system.hpp"
#include "particles_system.hpp"
#include "particles_input_hdf5.hpp"
#include "particles_interp_spline.hpp"

#include "field.hpp"
#include "kspace.hpp"



//////////////////////////////////////////////////////////////////////////////
///
/// Double template "for"
///
//////////////////////////////////////////////////////////////////////////////

namespace Template_double_for_if{

template <class RetType,
          class IterType1, IterType1 CurrentIter1,
          class IterType2, const IterType2 CurrentIter2, const IterType2 iterTo2, const IterType2 IterStep2,
          class Func, bool IsNotOver, typename... Args>
struct For2{
    static RetType evaluate(IterType2 value2, Args... args){
        if(CurrentIter2 == value2){
            return std::move(Func::template instanciate<CurrentIter1, CurrentIter2>(args...));
        }
        else{
            return std::move(For2<RetType,
                                        IterType1, CurrentIter1,
                                        IterType2, CurrentIter2+IterStep2, iterTo2, IterStep2,
                                        Func, (CurrentIter2+IterStep2 < iterTo2), Args...>::evaluate(value2, args...));
        }
    }
};

template <class RetType,
          class IterType1, IterType1 CurrentIter1,
          class IterType2, const IterType2 CurrentIter2, const IterType2 iterTo2, const IterType2 IterStep2,
          class Func, typename... Args>
struct For2<RetType,
                  IterType1, CurrentIter1,
                  IterType2, CurrentIter2, iterTo2, IterStep2,
                  Func, false, Args...>{
    static RetType evaluate(IterType2 value2, Args... args){
        std::cout << __FUNCTION__ << "[ERROR] template values for loop 2 " << value2 << " does not exist\n";
        return RetType();
    }
};

template <class RetType,
          class IterType1, const IterType1 CurrentIter1, const IterType1 iterTo1, const IterType1 IterStep1,
          class IterType2, const IterType2 IterFrom2, const IterType2 iterTo2, const IterType2 IterStep2,
          class Func, bool IsNotOver, typename... Args>
struct For1{
    static RetType evaluate(IterType1 value1, IterType2 value2, Args... args){
        if(CurrentIter1 == value1){
            return std::move(For2<RetType,
                                        IterType1, CurrentIter1,
                                        IterType2, IterFrom2, iterTo2, IterStep2,
                                        Func, (IterFrom2<iterTo2), Args...>::evaluate(value2, args...));
        }
        else{
            return std::move(For1<RetType,
                              IterType1, CurrentIter1+IterStep1, iterTo1, IterStep1,
                              IterType2, IterFrom2, iterTo2, IterStep2,
                              Func, (CurrentIter1+IterStep1 < iterTo1), Args...>::evaluate(value1, value2, args...));
        }
    }
};

template <class RetType,
          class IterType1, const IterType1 IterFrom1, const IterType1 iterTo1, const IterType1 IterStep1,
          class IterType2, const IterType2 IterFrom2, const IterType2 iterTo2, const IterType2 IterStep2,
          class Func, typename... Args>
struct For1<RetType,
                IterType1, IterFrom1, iterTo1, IterStep1,
                IterType2, IterFrom2, iterTo2, IterStep2,
                Func, false, Args...>{
    static RetType evaluate(IterType1 value1, IterType2 value2, Args... args){
        std::cout << __FUNCTION__ << "[ERROR] template values for loop 1 " << value1 << " does not exist\n";
        return RetType();
    }
};

template <class RetType,
          class IterType1, const IterType1 IterFrom1, const IterType1 iterTo1, const IterType1 IterStep1,
          class IterType2, const IterType2 IterFrom2, const IterType2 iterTo2, const IterType2 IterStep2,
          class Func, typename... Args>
inline RetType evaluate(IterType1 value1, IterType2 value2, Args... args){
    return std::move(For1<RetType,
            IterType1, IterFrom1, iterTo1, IterStep1,
            IterType2, IterFrom2, iterTo2, IterStep2,
            Func, (IterFrom1<iterTo1), Args...>::evaluate(value1, value2, args...));
}

}


//////////////////////////////////////////////////////////////////////////////
///
/// Builder Functions
///
//////////////////////////////////////////////////////////////////////////////

template <class field_rnumber, field_backend be, class particles_rnumber>
struct particles_system_build_container {
    template <const int interpolation_size, const int spline_mode>
    static std::unique_ptr<abstract_particles_system<particles_rnumber>> instanciate(
             const field<field_rnumber, be, THREE>* fs_field, // (field object)
             const kspace<be, SMOOTH>* fs_kk, // (kspace object, contains dkx, dky, dkz)
             const int nsteps, // to check coherency between parameters and hdf input file (nb rhs)
             const int nparticles, // to check coherency between parameters and hdf input file
             const std::string& fname_input, // particles input filename
            const std::string& inDatanameState, const std::string& inDatanameRhs, // input dataset names
             MPI_Comm mpi_comm){

        // The size of the field grid (global size) all_size seems
        std::array<size_t,3> field_grid_dim;
        field_grid_dim[IDX_X] = fs_field->rlayout->sizes[FIELD_IDX_X];// nx
        field_grid_dim[IDX_Y] = fs_field->rlayout->sizes[FIELD_IDX_Y];// nx
        field_grid_dim[IDX_Z] = fs_field->rlayout->sizes[FIELD_IDX_Z];// nz

        // The size of the local field grid (the field nodes that belong to current process)
        std::array<size_t,3> local_field_dims;
        local_field_dims[IDX_X] = fs_field->rlayout->subsizes[FIELD_IDX_X];
        local_field_dims[IDX_Y] = fs_field->rlayout->subsizes[FIELD_IDX_Y];
        local_field_dims[IDX_Z] = fs_field->rlayout->subsizes[FIELD_IDX_Z];

        // The offset of the local field grid
        std::array<size_t,3> local_field_offset;
        local_field_offset[IDX_X] = fs_field->rlayout->starts[FIELD_IDX_X];
        local_field_offset[IDX_Y] = fs_field->rlayout->starts[FIELD_IDX_Y];
        local_field_offset[IDX_Z] = fs_field->rlayout->starts[FIELD_IDX_Z];

        // Retreive split from fftw to know processes that have no work
        int my_rank, nb_processes;
        AssertMpi(MPI_Comm_rank(mpi_comm, &my_rank));
        AssertMpi(MPI_Comm_size(mpi_comm, &nb_processes));

        const int split_step = (int(field_grid_dim[IDX_Z])+nb_processes-1)/nb_processes;
        const int nb_processes_involved = (int(field_grid_dim[IDX_Z])+split_step-1)/split_step;

        assert((my_rank < nb_processes_involved && local_field_dims[IDX_Z] != 0)
               || (nb_processes_involved <= my_rank && local_field_dims[IDX_Z] == 0));
        assert(nb_processes_involved <= int(field_grid_dim[IDX_Z]));

        // Make the idle processes starting from the limit (and not 0 as set by fftw)
        if(nb_processes_involved <= my_rank){
            local_field_offset[IDX_Z] = field_grid_dim[IDX_Z];
        }

        // Ensure that 1D partitioning is used
        {
            assert(local_field_offset[IDX_X] == 0);
            assert(local_field_offset[IDX_Y] == 0);
            assert(local_field_dims[IDX_X] == field_grid_dim[IDX_X]);
            assert(local_field_dims[IDX_Y] == field_grid_dim[IDX_Y]);

            assert(my_rank >= nb_processes_involved || ((my_rank == 0 && local_field_offset[IDX_Z] == 0)
                   || (my_rank != 0 && local_field_offset[IDX_Z] != 0)));
            assert(my_rank >= nb_processes_involved || ((my_rank == nb_processes_involved-1 && local_field_offset[IDX_Z]+local_field_dims[IDX_Z] == field_grid_dim[IDX_Z])
                   || (my_rank != nb_processes_involved-1 && local_field_offset[IDX_Z]+local_field_dims[IDX_Z] != field_grid_dim[IDX_Z])));
        }

        // The offset of the local field grid
        std::array<size_t,3> local_field_mem_size;
        local_field_mem_size[IDX_X] = fs_field->rmemlayout->subsizes[FIELD_IDX_X];
        local_field_mem_size[IDX_Y] = fs_field->rmemlayout->subsizes[FIELD_IDX_Y];
        local_field_mem_size[IDX_Z] = fs_field->rmemlayout->subsizes[FIELD_IDX_Z];

        // The spatial box size (all particles should be included inside)
        std::array<particles_rnumber,3> spatial_box_width;
        spatial_box_width[IDX_X] = 4 * acos(0) / (fs_kk->dkx);
        spatial_box_width[IDX_Y] = 4 * acos(0) / (fs_kk->dky);
        spatial_box_width[IDX_Z] = 4 * acos(0) / (fs_kk->dkz);

        // The distance between two field nodes in z
        std::array<particles_rnumber,3> spatial_partition_width;
        spatial_partition_width[IDX_X] = spatial_box_width[IDX_X]/particles_rnumber(field_grid_dim[IDX_X]);
        spatial_partition_width[IDX_Y] = spatial_box_width[IDX_Y]/particles_rnumber(field_grid_dim[IDX_Y]);
        spatial_partition_width[IDX_Z] = spatial_box_width[IDX_Z]/particles_rnumber(field_grid_dim[IDX_Z]);
        // The spatial interval of the current process
        const particles_rnumber my_spatial_low_limit_z = particles_rnumber(local_field_offset[IDX_Z])*spatial_partition_width[IDX_Z];
        const particles_rnumber my_spatial_up_limit_z = particles_rnumber(local_field_offset[IDX_Z]+local_field_dims[IDX_Z])*spatial_partition_width[IDX_Z];

        // Create the particles system
        particles_system<particles_rnumber, field_rnumber, particles_interp_spline<particles_rnumber, interpolation_size,spline_mode>, interpolation_size>* part_sys
         = new particles_system<particles_rnumber, field_rnumber, particles_interp_spline<particles_rnumber, interpolation_size,spline_mode>, interpolation_size>(field_grid_dim,
                                                                                                   spatial_box_width,
                                                                                                   spatial_partition_width,
                                                                                                   my_spatial_low_limit_z,
                                                                                                   my_spatial_up_limit_z,
                                                                                                   fs_field->get_rdata(),
                                                                                                   local_field_dims,
                                                                                                   local_field_offset,
                                                                                                   local_field_mem_size,
                                                                                                   mpi_comm);

        // Load particles from hdf5
        particles_input_hdf5<particles_rnumber, 3,3> generator(mpi_comm, fname_input,
                                            inDatanameState, inDatanameRhs, my_spatial_low_limit_z, my_spatial_up_limit_z);

        // Ensure parameters match the input file
        if(generator.getNbRhs() != nsteps){
            std::runtime_error(std::string("Nb steps is ") + std::to_string(nsteps)
                               + " in the parameters but " + std::to_string(generator.getNbRhs()) + " in the particles file.");
        }
        // Ensure parameters match the input file
        if(generator.getTotalNbParticles() != nparticles){
            std::runtime_error(std::string("Nb particles is ") + std::to_string(nparticles)
                               + " in the parameters but " + std::to_string(generator.getTotalNbParticles()) + " in the particles file.");
        }

        // Load the particles and move them to the particles system
        part_sys->init(generator);

        assert(part_sys->getNbRhs() == nsteps);

        // Return the created particles system
        return std::unique_ptr<abstract_particles_system<particles_rnumber>>(part_sys);
    }
};


template <class field_rnumber, field_backend be, class particles_rnumber = double>
inline std::unique_ptr<abstract_particles_system<particles_rnumber>> particles_system_builder(
        const field<field_rnumber, be, THREE>* fs_field, // (field object)
        const kspace<be, SMOOTH>* fs_kk, // (kspace object, contains dkx, dky, dkz)
        const int nsteps, // to check coherency between parameters and hdf input file (nb rhs)
        const int nparticles, // to check coherency between parameters and hdf input file
        const std::string& fname_input, // particles input filename
        const std::string& inDatanameState, const std::string& inDatanameRhs, // input dataset names
        const int interpolation_size,
        const int spline_mode,
        MPI_Comm mpi_comm){
    return Template_double_for_if::evaluate<std::unique_ptr<abstract_particles_system<particles_rnumber>>,
                       int, 1, 7, 1, // interpolation_size
                       int, 0, 3, 1, // spline_mode
                       particles_system_build_container<field_rnumber,be,particles_rnumber>>(
                           interpolation_size, // template iterator 1
                           spline_mode, // template iterator 2
                           fs_field,fs_kk, nsteps, nparticles, fname_input, inDatanameState, inDatanameRhs, mpi_comm);
}


#endif