gyro_avg_mat.hpp

Go to the documentation of this file.

#ifndef GYRO_AVG_MAT_HPP
#define GYRO_AVG_MAT_HPP

#include "space_settings.hpp"
#include "matrix.hpp"
#include "task_group.hpp"
#include "host_array.hpp"
#include "grid.hpp"
#include "domain_decomposition.hpp"

enum class GyroCommMethod{
    PlaneComm=0,
    PackedPlaneComm,
    SelfComm
};

template<class Device>
class GyroAverageMatrices{
    public: // for testing
    HostArray<Matrix<Device>> matrices;
    HostArray<Matrix<HostType>> matrices_h;

    int nrho;
    double rhomax;
    double inv_drho;
    int nphi;
    GyroCommMethod gyro_comm_method;
    bool resident_on_device;

    TaskGroup task_group;

    public:

    // Default constructor used if species is drift-kinetic
    GyroAverageMatrices()
      : nrho(0),
        gyro_comm_method(GyroCommMethod::PlaneComm)
    {}

    // Construct the gyromatrices used on this MPI rank
    GyroAverageMatrices(NLReader::NamelistReader& nlr, const Grid<DeviceType>& grid, const DomainDecomposition<DeviceType>& pol_decomp, const MagneticField<DeviceType>& magnetic_field, bool is_field_following, int n_larmor_pts, int nrho_in, double rhomax_in);

    // Copy matrices to device
    void copy_matrices_to_device(){
        for(int i=0; i<matrices.size(); i++){
            matrices(i) = matrices_h(i).template mirror<Device>();
        }
    }

    void copy_to_device_if_not_resident(){
        if(!resident_on_device){
            copy_matrices_to_device();
        }
    }

    void deallocate_device_matrices_if_not_resident(){
        if(!resident_on_device){
            for(int i=0; i<matrices.size(); i++){
                matrices(i) = Matrix<Device>();
            }
        }
    }

    GyroCommMethod read_gyro_comm_method(NLReader::NamelistReader& nlr) const;

#ifdef USE_MPI
    MPI_Comm choose_comm(const DomainDecomposition<DeviceType>& pol_decomp, int ntasks) const;
#endif

    // Could be a constructor of a new GyroMatrix class (which could be a derived Matrix class)
    void init_one_gyro_avg_mat(const Grid<DeviceType>& grid, const MagneticField<DeviceType>& magnetic_field, Matrix<HostType>& matrix, int iphi, int irho, bool is_field_following, int n_larmor_pts);

    void gyroaverage_field(const Grid<DeviceType>& grid, const View<double***,CLayout,DeviceType>& field_rho,
#ifdef USE_MPI
                           const MPI_Comm& plane_comm,
#endif
                           const View<double**,CLayout,HostType>& field) const;

    void get_pot_mat_mult(const View<double***,CLayout, DeviceType>& field_rho, const View<double****,CLayout, DeviceType>& gradient_rho, const View<double***,CLayout, DeviceType>& E_para_em, bool potential_requested, bool gradient_requested, bool E_para_em_requested) const;

    void broadcast_potentials(const View<double***,CLayout, DeviceType>& potential_rho, const View<double****,CLayout, DeviceType>& gradient_rho, const View<double***,CLayout, DeviceType>& E_para_em, bool potential_requested, bool gradient_requested, bool E_para_em_requested) const;

    void get_gyroaveraged(const View<double***,CLayout, DeviceType>& field_rho, const View<double****,CLayout, DeviceType>& gradient_rho, const View<double***,CLayout, DeviceType>& E_para_em, bool potential_requested, bool gradient_requested, bool E_para_em_requested) const;

    void gyroaverage_moments(const View<double**, CLayout, DeviceType>& density_rho,
                                                      const View<double**, CLayout, DeviceType>& flow_rho,
                                                      const View<double**, CLayout, DeviceType>& energy_rho,
                                                      const View<double**, CLayout, DeviceType>& density_f0g_rho,
                                                      const View<double*, CLayout, DeviceType>& density,
                                                      const View<double*, CLayout, DeviceType>& flow,
                                                      const View<double*, CLayout, DeviceType>& energy, 
                                                      const View<double*, CLayout, DeviceType>& density_f0g
#ifdef USE_MPI
                                                      , const MPI_Comm& plane_comm
#endif
                                                      );


    int find_a_plane_rank_with_last_gyromatrix() const;

    void iterative_gyroaverage(Matrix<DeviceType>& gyromatrix, bool print_residual, View<double*, CLayout, DeviceType>& view) const;

    void smooth_with_last_gyroaverage(bool print_residual, View<double*, CLayout, DeviceType>& delta_n) const;
};

#endif