charge.hpp

Go to the documentation of this file.

#ifndef CHARGE_HPP
#define CHARGE_HPP
#include "space_settings.hpp"
#include "electric_field.hpp"

extern "C" void set_electron_charge_pointers(int nnode, int npsi00, int phi_ind0, double* density_h, double* density0_h, double* den00_1d_h);
extern "C" void set_electron_charge_pointers_em(int nnode, double* jpar_h, double* jpar0_h);
extern "C" void set_ion_charge_pointers(int nnode, int npsi00, int phi_ind0, double* density_h, double* density0_h, double* den00_1d_h);
extern "C" void set_ion_charge_pointers_em(int nnode, double* jpar_h, double* jpar0_h);
extern "C" void set_cden00_1d_fort(int npsi00, double* cden00_1d_h);

// Charge class
template<class Device, KinType KT>
class Charge {
    public:

    bool is_electrons;
    bool use_current; 

    // Host views
    View<double**,CLayout,HostType> density_h;
    View<double*,CLayout,HostType> density0_h;
    View<double*,CLayout,HostType> den00_1d_h;

    View<double**,CLayout,HostType> jpar_h;
    View<double*,CLayout,HostType> jpar0_h;
    
    View<double**,CLayout,HostType> den_f0_h;

    // ions only
    View<double***,CLayout,HostType> den_rho_ff_h;
    View<double***,CLayout,HostType> jpar_rho_ff_h;

    // Device views
    GridField<Device, VarType::Scalar,PIT_GLOBAL, TorType::OnePlane, KT, SCATTER_TYPE_GLOBAL> density; 
    GridField<Device, VarType::Scalar,PIT_GLOBAL, TorType::OnePlane, KT, SCATTER_TYPE_GLOBAL> jpar;    

    // Default constructor
    Charge(){}

    // Constructor
    Charge(bool use_current, int nnode, int nphi, int nrho, int npsi00, bool is_electrons)
      : is_electrons(is_electrons),
        use_current(use_current),
          // Host views
          density_h("density_h", nphi, nnode),
          density0_h("density0_h", nnode),
          den00_1d_h("den00_1d_h", npsi00),

          jpar_h("jpar_h", nphi, use_current ? nnode : 0),
          jpar0_h("jpar0_h", use_current ? nnode : 0),

          den_f0_h("den_f0_h", nrho+1, nnode),

          // ions only
          den_rho_ff_h("den_rho_ff_h", nrho+1, nphi, nnode),
          jpar_rho_ff_h("jpar_rho_ff_h", nrho+1, nphi, use_current ? nnode : 0),

          // Device views
          jpar("jpar", 1, nrho, use_current ? nnode : 0), // nphi=1
          density("density", 1, nrho, nnode) // nphi=1
    {
#ifndef NO_FORTRAN_MODULES
        int phi_ind0 = 2 - nphi; // Starting index for fortran density = 1 for XGCA, 0 for XGC1
        if(is_electrons){
            set_electron_charge_pointers(nnode, npsi00, phi_ind0, density_h.data(), density0_h.data(), den00_1d_h.data());
            if(use_current) set_electron_charge_pointers_em(nnode, jpar_h.data(), jpar0_h.data());
        }else{
            set_ion_charge_pointers(nnode, npsi00, phi_ind0, density_h.data(), density0_h.data(), den00_1d_h.data());
            if(use_current) set_ion_charge_pointers_em(nnode, jpar_h.data(), jpar0_h.data());
        }
#endif
    }

    void reset_to_zero(){
        density.reset_to_zero();
        if(use_current) jpar.reset_to_zero();
    }
};

struct Charges{
    Charge<DeviceType, KinType::DriftKin> electron;
    Charge<DeviceType, KinType::GyroKin> ion;

    View<double*,CLayout,HostType> cden00_1d_h;

    Charges(const Grid<DeviceType>& grid, const Simulation<DeviceType>& sml, int nrho){
        int nphi = (sml.is_XGCa ? 1 : 2);
        electron = Charge<DeviceType,KinType::DriftKin>(sml.explicit_electromagnetic, grid.nnode, nphi, 0, grid.psi00.n, true);
        ion = Charge<DeviceType,KinType::GyroKin>(sml.explicit_electromagnetic, grid.nnode, nphi, nrho, grid.psi00.n, false);
        cden00_1d_h = View<double*,CLayout,HostType>("cden00_1d_h", grid.psi00.n);
#ifndef NO_FORTRAN_MODULES
        set_cden00_1d_fort(grid.psi00.n, cden00_1d_h.data());
#endif
    }
};

#endif