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"

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

    int nrho;
    double rhomax;
    double inv_drho;
    int nphi;
    TaskGroup task_group;

    // Access individual matrices (may not be necessary eventually)
    Matrix<Device> operator [](int i) const {return matrices(i);}
    Matrix<Device>& operator [](int i) {return matrices(i);}

    public:

    GyroAverageMatrices(){} // For testing

    // 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);

    // Create a mirror with a different device type
    template<class Device2>
    inline GyroAverageMatrices<Device2> mirror() const{
        GPTLstart("GET_POT_GRAD_GYROMAT_SETUP");
        GyroAverageMatrices<Device2> m;
        const std::string label ("v_outer");
        m.matrices = HostArray<Matrix<Device2>>(matrices.size());
        for(int i=0; i<matrices.size(); i++){
            m.matrices(i) = matrices(i).template mirror<Device2>();
        }

        m.nrho = nrho;
        m.rhomax = rhomax;
        m.inv_drho = inv_drho;
        m.nphi = nphi;
        m.task_group = task_group; // Can't be used from a device kernel, but no harm in sending it

        GPTLstop("GET_POT_GRAD_GYROMAT_SETUP");

        return m;
    }

    // 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<Device>& matrix, int iphi, int irho, bool is_field_following, int n_larmor_pts);

    void gyroaverage_field(const Grid<DeviceType>& grid, View<double***,CLayout,HostType>& field_rho
#ifdef USE_MPI
                           , const MPI_Comm& plane_comm
#endif
                           ) 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
                                                      ) const;


    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