electric_field.hpp

Go to the documentation of this file.

#ifndef ELECTRIC_FIELD_HPP
#define ELECTRIC_FIELD_HPP
#include "space_settings.hpp"
#include "NamelistReader.hpp"
#include "grid.hpp"
#include "grid_field_pack.hpp"
#include "timer_macro.hpp"
#include "smoothing.hpp"
#include "get_potential_grad.hpp"
#include "plane_field_gatherer.hpp"
#include "field_following_coordinates.hpp"
#ifdef SONIC_GK
#include "get_sonic_fields.hpp"
#endif
#ifdef EXPLICIT_EM
#include "globals.hpp"
#endif
#include "bias_potential.hpp"
#include "simple00_solver.hpp"

template<KinType KT, PhiInterpType PIT, MarkerType MT, MagneticFieldMode MFM>
struct gfpack_push_type;

template<KinType KT> struct getGyroVecFldType;
template<> struct getGyroVecFldType<KinType::DriftKin>{ using t = EfieldType; };
template<> struct getGyroVecFldType<KinType::GyroKin>{ using t = EfieldGyroType; };

template<KinType KT> struct getGyroScaFldType;
template<> struct getGyroScaFldType<KinType::DriftKin>{ using t = PotType; };
template<> struct getGyroScaFldType<KinType::GyroKin>{ using t = PotGyroType; };

// XGC1

template<KinType KT, MarkerType MT>
struct gfpack_push_type<KT, PhiInterpType::Planes, MT, MagneticFieldMode::Electromagnetic> {
    using t = Pack<Labeled<typename getGyroVecFldType<KT>::t, Label::E>,
                   Labeled<typename getGyroVecFldType<KT>::t, Label::dAh>,
                   Labeled<typename getGyroVecFldType<KT>::t, Label::dAs>,
                   Labeled<typename getGyroScaFldType<KT>::t, Label::Ah>,
                   Labeled<typename getGyroScaFldType<KT>::t, Label::As>,
                   Labeled<typename getGyroScaFldType<KT>::t, Label::Epar_em>>;
};

template<KinType KT, MarkerType MT>
struct gfpack_push_type<KT, PhiInterpType::Planes, MT, MagneticFieldMode::Electrostatic> {
    using t = Pack<Labeled<typename getGyroVecFldType<KT>::t, Label::E>>;
};

template<KinType KT>
struct gfpack_push_type<KT, PhiInterpType::Planes, MarkerType::ReducedDeltaF, MagneticFieldMode::Electrostatic> {
    using t = Pack<Labeled<typename getGyroVecFldType<KT>::t, Label::E>,
                   Labeled<E00Type, Label::E00>,
                   Labeled<PotType, Label::ddpotdt>>;
};

// XGCA
template<KinType KT, MarkerType MT>
struct gfpack_push_type<KT, PhiInterpType::None, MT, MagneticFieldMode::Electrostatic> {
# ifdef SONIC_GK
    using t = Pack<Labeled<Efield2DGyroType, Label::E>,
                   Labeled<E00GyroType, Label::dEr_B2>,
                   Labeled<E00GyroType, Label::dEz_B2>,
                   Labeled<E00GyroType, Label::du2_E>>;
# else
    using t = Pack<Labeled<Efield2DGyroType, Label::E>>;
# endif
};


// Fortran calls
#ifdef XGC1
#  ifdef EXPLICIT_EM
extern "C" void set_rho_ff_pointers(int n, Field<VarType::Scalar,PhiInterpType::None>* dpot, Field<VarType::Scalar,PhiInterpType::Planes>* dpot_ff,
      Field<VarType::Scalar,PhiInterpType::None>* Ah, Field<VarType::Scalar,PhiInterpType::None>* As, Field<VarType::Scalar,PhiInterpType::None>* As_backup);
#  else
extern "C" void set_rho_ff_pointers(int n, Field<VarType::Scalar,PhiInterpType::None>* dpot, Field<VarType::Scalar,PhiInterpType::Planes>* dpot_ff);
#  endif
#else
extern "C" void set_rho_pointers(int n, Field<VarType::Scalar,PhiInterpType::None>* dpot);
#endif

// Electric field class
template<class Device>
struct ElectricField {

    // Dimensions
    int nnode;
    int nphi;

    // Gyromatrix info (to be moved soon)
    int nrho;

    // Optimization options
    bool calculate_phi_ff_on_device;

    // Field variables
    bool turb_efield; 
    bool E00_efield;   
    bool n0_m_efield;

    // Host Views
    GridField<HostType,VarType::Vector2D,PIT_GLOBAL,TorType::OnePlane,KinType::DriftKin> E00_ff_h; 
    GridField<HostType,VarType::Scalar,PhiInterpType::Planes,TorType::OnePlane,KinType::DriftKin> ddpotdt; 

    //For EM
    GridField<HostType,VarType::Scalar,PhiInterpType::None,TorType::MultiplePlanes,KinType::DriftKin> Ah_h;
    GridField<HostType,VarType::Scalar,PhiInterpType::None,TorType::OnePlane,KinType::DriftKin> Ah_h_backup;
    GridField<HostType,vec2d_if_axisym<PIT_GLOBAL>(),PIT_GLOBAL,TorType::OnePlane,KinType::GyroKin> dAh;
    GridField<HostType,VarType::Scalar,PIT_GLOBAL,TorType::OnePlane,KinType::GyroKin> Ah;
    GridField<HostType,VarType::Scalar,PhiInterpType::None,TorType::MultiplePlanes,KinType::DriftKin> As_h;
    GridField<HostType,VarType::Scalar,PhiInterpType::None,TorType::MultiplePlanes,KinType::DriftKin> As_h_backup;
    GridField<HostType,vec2d_if_axisym<PIT_GLOBAL>(),PIT_GLOBAL,TorType::OnePlane,KinType::GyroKin> dAs;
    GridField<HostType,VarType::Scalar,PIT_GLOBAL,TorType::OnePlane,KinType::GyroKin> As;
    GridField<HostType,VarType::Scalar,PhiInterpType::None,TorType::MultiplePlanes,KinType::DriftKin> Ah_cv_h;
    GridField<HostType,VarType::Scalar,PhiInterpType::Planes,TorType::OnePlane,KinType::DriftKin> Ah_cv;
    GridField<HostType,VarType::Scalar,PhiInterpType::Planes,TorType::OnePlane,KinType::GyroKin> Epar_em;

    GridField<HostType,VarType::Scalar,PIT_GLOBAL,TorType::OnePlane,KinType::GyroKin> pot;
    GridField<HostType,vec2d_if_axisym<PIT_GLOBAL>(),PIT_GLOBAL,TorType::OnePlane,KinType::GyroKin> E;
//#ifdef XGC1
    GridField<HostType,VarType::Scalar,PhiInterpType::Planes,TorType::OnePlane,KinType::DriftKin> dpot_ff_h;
//#else
//#  ifdef SONIC_GK
    GridField<HostType,VarType::Vector2D,PIT_GLOBAL,TorType::OnePlane,KinType::GyroKin> dEr_B2; 
    GridField<HostType,VarType::Vector2D,PIT_GLOBAL,TorType::OnePlane,KinType::GyroKin> dEz_B2; 
    GridField<HostType,VarType::Vector2D,PIT_GLOBAL,TorType::OnePlane,KinType::GyroKin> du2_E ; 
//#  endif
//#endif
    GridField<HostType,VarType::Scalar,PhiInterpType::None,TorType::MultiplePlanes,KinType::DriftKin> dpot_h;
    GridField<HostType,VarType::Scalar,PhiInterpType::None,TorType::OnePlane,KinType::DriftKin> pot0_h;
    GridField<HostType,VarType::Scalar,PhiInterpType::None,TorType::OnePlane,KinType::DriftKin> dpot_n0_h;
    View<double***,CLayout, HostType> save_dpot0; // Used for adiabatic response
    View<double***,CLayout, HostType> save_dpot; // Used for adiabatic response
    View<double***,CLayout, HostType> dpotsave; // Used for split-weight scheme
    View<double*,CLayout, HostType> dpot_es; // Used for a diagnostic

    View<double*,CLayout, HostType> pot0m; 
    View<double*,CLayout, HostType> pot00_1d; 

    Simple00Solver simple00;

    // Host View for the adjustable loop voltage current drive
    GridField<HostType,VarType::Scalar,PhiInterpType::None,TorType::OnePlane,KinType::DriftKin> loop_voltage_h;
    GridField<HostType,VarType::Scalar,PhiInterpType::None,TorType::OnePlane,KinType::DriftKin> jpara_diff_previous_h;
    GridField<HostType,VarType::Scalar,PhiInterpType::None,TorType::OnePlane,KinType::DriftKin> jpara_diff_integral_h;

    // Bias potential
    BiasPotential bias_potential;

    // Solver boundaries; These are lists of nodes that are EXCLUDED
    Boundary jbdH;
    Boundary AbdH;
    Boundary cbd0;
    Boundary pbd0;
    Boundary cbdH;
    Boundary pbdH;

    // Default constructor
    ElectricField(){}

    // Makes sense to use the Grid in the constructor here, but the Grid still copies from fortran so until that changes,
    // it's simpler to not use it here so that ElectricField can be defined before the restart, which uses psn%As
    ElectricField(NLReader::NamelistReader& nlr, const Grid<DeviceType>& grid, int nrho_in, bool es_reduced_deltaf_setup=false, double dt=1.0, double psi_norm=1.0) : //Grid<HostTypeType>& grid) :
      nnode(grid.nnode),
      nphi(grid.nplanes),
      nrho(nrho_in),
      // Host Views
#ifdef XGC1
      pot("pot",nphi,nrho,nnode),
      E("E",nphi,nrho,nnode),
      save_dpot("save_dpot", 2, 2, nnode), // View. dims: RK2, ff, nnode
      save_dpot0("save_dpot0", 2, 2, nnode), // View. dims: RK2, ff, nnode
#  ifdef EXPLICIT_EM
      dpot_es("dpot_es", nnode),
      dAh("dAh",nphi,nrho,nnode),
      Ah("Ah",nphi,nrho,nnode),
      Ah_h("Ah_h",4,nnode), // -1:2 in fortran
      Ah_h_backup("Ah_h_backup",nnode), // -1:2 in fortran
      As_h("As_h",4,nnode), // -1:2 in fortran
      As_h_backup("As_h_backup",4,nnode), // -1:2 in fortran
#  endif
      dpot_h("dpot_h",4,nnode), // -1:2 in fortran
      dpot_ff_h("dpot_ff_h",nnode),
      dpot_n0_h("dpot_n0_h",nnode),
#else
      E("E",nrho, nnode),
      dpot_h("dpot_h",1,nnode),
#  ifdef SONIC_GK
      dEr_B2("dEr_B2",nrho, nnode),
      dEz_B2("dEz_B2",nrho, nnode),
      du2_E("du2_E",nrho, nnode),
#  endif
      dpot_n0_h("dpot_n0_h",nnode),
#endif
      pot0m("pot0m", nnode),
      pot00_1d("pot00_1d", grid.psi00.n),
      pot0_h("pot0_h",nnode),
      loop_voltage_h("loop_voltage_h",nnode),
      jpara_diff_previous_h("jpara_diff_previous_h",nnode),
      jpara_diff_integral_h("jpara_diff_integral_h",nnode),
      bias_potential(nlr, dt, psi_norm)
    {
        if(es_reduced_deltaf_setup){
            E00_ff_h = GridField<HostType,VarType::Vector2D,PIT_GLOBAL,TorType::OnePlane,KinType::DriftKin>("E00_ff_h",nnode);
            ddpotdt = GridField<HostType,VarType::Scalar,PhiInterpType::Planes,TorType::OnePlane,KinType::DriftKin>("ddpotdt",nphi,nrho,nnode);
            dpotsave = View<double***,CLayout, HostType>("dpotsave", 2, 2, nnode); // View. dims: RK2, ff, nnode
        }

        nlr.use_namelist("sml_param");
        turb_efield = nlr.get<bool>("sml_turb_efield",true); // E-field calculated only with \f$<\phi>\f$ if .false.
                                                             // psn%dpot will still contain all (n=0,|m|>0) and (|n|>0,m) modes
        E00_efield = nlr.get<bool>("sml_00_efield",true); // Flux-surface averaged potential not used for calculating the electric field if .false.
        n0_m_efield = nlr.get<bool>("sml_0m_efield",true); // The axisymmetric potential variation \f$\langle \phi-\langle\phi\rangle\rangle_T\f$ is
                                              // set to zero if .false.

#ifdef EXPLICIT_EM
        // Pass these (or sml) as arguments
        bool em_mixed_variable    = nlr.get<bool>("sml_em_mixed_variable",true);
        bool em_control_variate   = nlr.get<bool>("sml_em_control_variate",false);
        int  em_pullback_mode     = nlr.get<int>("sml_em_pullback_mode",3);

        if (!em_mixed_variable)
            exit_XGC("\nError: sml_em_mixed_variable is explicitly set to .false. which is incompatible with EXPLICIT_EM\n");
#endif

        if(nlr.namelist_present("performance_param")){
            nlr.use_namelist("performance_param"); 
            calculate_phi_ff_on_device = nlr.get<bool>("calculate_phi_ff_on_device", false); // With this option, potential is communicated and the electric field is calculated on every MPI process. This means additional computation, but less communication. This tradeoff is probably beneficial on all current flagship platforms, but this has not been verified.
        }else{
            calculate_phi_ff_on_device = false;
        }

#ifdef EXPLICIT_EM
        if (em_control_variate){
            Ah_cv_h = GridField<HostType,VarType::Scalar,PhiInterpType::None,TorType::MultiplePlanes,KinType::DriftKin>("Ah_cv_h",4,nnode); // -1:2 in fortran
            Ah_cv = GridField<HostType,VarType::Scalar,PhiInterpType::Planes,TorType::OnePlane,KinType::DriftKin>("Ah_cv",nphi,nrho,nnode);
        }

        if (em_mixed_variable){
            As = GridField<HostType,VarType::Scalar,PIT_GLOBAL,TorType::OnePlane,KinType::GyroKin>("As",nphi,nrho,nnode);
            dAs = GridField<HostType,vec2d_if_axisym<PIT_GLOBAL>(),PIT_GLOBAL,TorType::OnePlane,KinType::GyroKin>("dAs",nphi,nrho,nnode);
            if (em_pullback_mode==4){
                Epar_em = GridField<HostType,VarType::Scalar,PhiInterpType::Planes,TorType::OnePlane,KinType::GyroKin>("Epar_em",nphi,nrho,nnode);
            }
        }
#endif

        // Set pointers in Fortran to rho arrays
#ifdef XGC1
#  ifdef EXPLICIT_EM
        set_rho_ff_pointers(nnode, dpot_h.data(), dpot_ff_h.data(), Ah_h.data(), As_h.data(), As_h_backup.data());
#  else
        set_rho_ff_pointers(nnode, dpot_h.data(), dpot_ff_h.data());
#  endif
#else
        set_rho_pointers(nnode, dpot_h.data());
#endif
        // Initialize loop voltage
        nlr.use_namelist("src_param");
        const double initial_loop_vol = nlr.get<double>("src_loop_voltage",0.0);
        //
        for (int inode=0; inode<nnode; inode++){
             loop_voltage_h(inode) = initial_loop_vol;
        }
    }


    /**************************************************/
    /************  Copy to GridFieldPack  *************/
    /**************************************************/
    template<MarkerType MT>
    GridFieldPackPtr copy_push_fields_to_device_drift_kin_XGC1(const Simulation<DeviceType>& sml, const DomainDecomposition<DeviceType>& pol_decomp, const MagneticField<DeviceType>& magnetic_field, const Grid<DeviceType>& grid, Smoothing& smoothing, bool near_field=false) const{

        // Set push type here for now
        using gfp_push_type = GridFieldPack<Device, typename gfpack_push_type<KinType::DriftKin, PIT_GLOBAL, MT, MFM_GLOBAL>::t>;
        gfp_push_type gfpack(turb_efield);

        GPTLstart("GAT_FIELD_INFO");
        if(pol_decomp.decompose_fields){
            if(near_field){
                // Gathering the "near field", defined as the field decomposition region that contains this rank's
                // global domain decomposition
                int local_pid = pol_decomp.field_decomp.find_domain_owner(pol_decomp.plane_index, grid.nplanes, pol_decomp.node_offset, grid.nnode);
                gfpack.phi_offset = pol_decomp.field_decomp.all_first_plane(local_pid);
                gfpack.node_offset = pol_decomp.field_decomp.all_first_node(local_pid);
            }else{
                gfpack.phi_offset = pol_decomp.field_decomp.first_plane;
                gfpack.node_offset = pol_decomp.field_decomp.first_node;
            }
        }
        PlaneFieldGatherer pfg(pol_decomp, grid, near_field);

#ifndef EXPLICIT_EM
        if constexpr(MT==MarkerType::ReducedDeltaF){
            if(ddpotdt.f.size()>0) pfg.gather_phi_ff_on_device(ddpotdt.f, gfpack.template get<Label::ddpotdt>().f);
            grid_field_copy(gfpack.template get<Label::E00>(), E00_ff_h); // No gather_phi needed due to axisymmetry
        }
#endif

        if(calculate_phi_ff_on_device){
            // Allocate temporary views
            GPTLstart("GET_POT_GRAD_ALLOC_VIEWS");
            int n_input_field_planes = pot0_h.nphi();
            bool gradient_requested = true;
            bool gyroaverage_requested = false;
            GetPotentialGradTemp<DeviceType, DeviceType> gpg_tmp(sml, grid, pol_decomp, magnetic_field, pbd0, n_input_field_planes, gradient_requested, gyroaverage_requested);
            GPTLstop("GET_POT_GRAD_ALLOC_VIEWS");

            // Whether to ignore dpot
            bool ignore_poloidal_dpot = !(turb_efield);
            bool use_field00 = E00_efield;
            bool potential_is_requested = false; // Not used in push
            GridField<Device, VarType::Scalar,PhiInterpType::Planes, TorType::MultiplePlanes, KinType::DriftKin> unused_pot("unused_pot",0,0);
            pfg.calculate_phi_ff_on_device(sml,grid, magnetic_field, smoothing, gpg_tmp, my_subview(dpot_h.f, 2), pot0_h, gfpack.template get<Label::E>(), unused_pot, potential_is_requested, use_field00, ignore_poloidal_dpot);
#  ifdef EXPLICIT_EM
            GridField<HostType,VarType::Scalar,PhiInterpType::None,TorType::OnePlane,KinType::DriftKin> pot0_em_h("pot0_em_h",0); // unused
            pfg.calculate_phi_ff_on_device(sml,grid, magnetic_field, smoothing, gpg_tmp, my_subview(Ah_h.f, 2), pot0_em_h, gfpack.template get<Label::dAh>(), gfpack.template get<Label::Ah>());
            if(As.f.size()>0) pfg.calculate_phi_ff_on_device(sml,grid, magnetic_field, smoothing, gpg_tmp, my_subview(As_h.f, 2), pot0_em_h, gfpack.template get<Label::dAs>(), gfpack.template get<Label::As>());
#  endif
        }else{
            pfg.gather_phi_ff_on_device(E.f, gfpack.template get<Label::E>().f);
#  ifdef EXPLICIT_EM
            pfg.gather_phi_ff_on_device(dAh.f, gfpack.template get<Label::dAh>().f);
            pfg.gather_phi_ff_on_device(Ah.f, gfpack.template get<Label::Ah>().f);
            if(dAs.f.size()>0) pfg.gather_phi_ff_on_device(dAs.f, gfpack.template get<Label::dAs>().f);
            if(As.f.size()>0) pfg.gather_phi_ff_on_device(As.f, gfpack.template get<Label::As>().f);
#  endif
        }
#  ifdef EXPLICIT_EM
        if(Epar_em.f.size()>0) pfg.gather_phi_ff_on_device(Epar_em.f, gfpack.template get<Label::Epar_em>().f);
#  endif
        GPTLstop("GAT_FIELD_INFO");

        FieldFollowingCoordinates ff(grid.half_plane_ff);
        // Copy axisymmetric loop voltage to temporary device-view
        ScalarGridField loop_voltage("loop_voltage",loop_voltage_h.nnode());
        Kokkos::deep_copy(loop_voltage.f,loop_voltage_h.f);
        gfpack.loop_voltage = GridField<Device, VarType::Scalar,PhiInterpType::Planes, TorType::OnePlane, KinType::DriftKin>("gfpack_loop_voltage",loop_voltage.nnode());
        ff.cnvt_grid_real2ff(grid,loop_voltage,gfpack.loop_voltage);

        return std::make_unique<gfp_push_type>(gfpack);
    }

    template<MarkerType MT>
    GridFieldPackPtr copy_push_fields_to_device_gyro_kin_XGC1(const Simulation<DeviceType>& sml, const DomainDecomposition<DeviceType>& pol_decomp, const MagneticField<DeviceType>& magnetic_field, const Grid<DeviceType>& grid, Smoothing& smoothing, bool near_field=false) const{
        using gfp_push_type = GridFieldPack<Device, typename gfpack_push_type<KinType::GyroKin, PIT_GLOBAL, MT, MFM_GLOBAL>::t>;
        gfp_push_type gfpack(turb_efield);

        GPTLstart("Copy_rho_ff_fields_to_device");
        grid_field_copy(gfpack.template get<Label::E>(), E);
#  ifdef EXPLICIT_EM
        grid_field_copy(gfpack.template get<Label::dAs>(), dAs);
        grid_field_copy(gfpack.template get<Label::dAh>(), dAh);
        grid_field_copy(gfpack.template get<Label::As>(), As);
        grid_field_copy(gfpack.template get<Label::Ah>(), Ah);
        grid_field_copy(gfpack.template get<Label::Epar_em>(), Epar_em);
#  endif
        GPTLstop("Copy_rho_ff_fields_to_device");

        // Loop voltage is identical for gyrokinetic and drift-kinetic species,
        // but different between XGC1 and XGCa
        FieldFollowingCoordinates ff(grid.half_plane_ff);
        // Copy axisymmetric loop voltage to temporary device-view
        ScalarGridField loop_voltage("loop_voltage",loop_voltage_h.nnode());
        Kokkos::deep_copy(loop_voltage.f,loop_voltage_h.f);
        gfpack.loop_voltage = GridField<Device, VarType::Scalar,PhiInterpType::Planes, TorType::OnePlane, KinType::DriftKin>("gfpack_loop_voltage",loop_voltage.nnode());
        ff.cnvt_grid_real2ff(grid,loop_voltage,gfpack.loop_voltage);

        return std::make_unique<gfp_push_type>(gfpack);
    }

#ifdef XGC1
    // Resizes allocations for large arrays that have a phi dimension (i.e used for electron push)
    GridFieldPackPtr copy_push_fields_to_device(KinType kintype_in, MarkerType markertype_in, const Simulation<DeviceType>& sml, const DomainDecomposition<DeviceType>& pol_decomp, const MagneticField<DeviceType>& magnetic_field, const Grid<DeviceType>& grid, Smoothing& smoothing, bool near_field=false) const{
        if(markertype_in==MarkerType::ReducedDeltaF){
            if(kintype_in==DriftKin){
                return copy_push_fields_to_device_drift_kin_XGC1<MarkerType::ReducedDeltaF>(sml, pol_decomp, magnetic_field, grid, smoothing, near_field);
            }else{
                return copy_push_fields_to_device_gyro_kin_XGC1<MarkerType::ReducedDeltaF>(sml, pol_decomp, magnetic_field, grid, smoothing, near_field);
            }
        }else{
            if(kintype_in==DriftKin){
                return copy_push_fields_to_device_drift_kin_XGC1<MarkerType::TotalF>(sml, pol_decomp, magnetic_field, grid, smoothing, near_field);
            }else{
                return copy_push_fields_to_device_gyro_kin_XGC1<MarkerType::TotalF>(sml, pol_decomp, magnetic_field, grid, smoothing, near_field);
            }
        }
    }
#else
    // Resizes allocations for large arrays that have a phi dimension (i.e used for electron push)
    GridFieldPackPtr copy_push_fields_to_device(KinType kintype_in, MarkerType markertype_in, const Simulation<DeviceType>& sml, const DomainDecomposition<DeviceType>& pol_decomp, const MagneticField<DeviceType>& magnetic_field, const Grid<DeviceType>& grid, Smoothing& smoothing, bool near_field=false) const{
        // Re-initialize new gridpack - deallocates whatever Views were already allocated
        using gfp_push_type = GridFieldPack<Device, typename gfpack_push_type<KinType::GyroKin, PIT_GLOBAL, MarkerType::TotalF, MFM_GLOBAL>::t>;
        gfp_push_type gfpack(turb_efield);

        grid_field_copy(gfpack.template get<Label::E>(), E);
#  ifdef SONIC_GK
        grid_field_copy(gfpack.template get<Label::dEr_B2>, dEr_B2);
        grid_field_copy(gfpack.template get<Label::dEz_B2>, dEz_B2);
        grid_field_copy(gfpack.template get<Label::du2_E>, du2_E);
#  endif
        // Copy axisymmetric loop voltage to temporary device-view
        ScalarGridField loop_voltage("loop_voltage",loop_voltage_h.nnode());
        Kokkos::deep_copy(loop_voltage.f,loop_voltage_h.f);
        grid_field_copy(gfpack.loop_voltage,loop_voltage);

        //grid_field_copy(gfpack.loop_voltage,loop_voltage_h);    
        return std::make_unique<gfp_push_type>(gfpack);
    }
#endif

    void get_Ah_rho_ff(const Simulation<DeviceType>& sml, const Grid<DeviceType>& grid, const DomainDecomposition<DeviceType>& pol_decomp, const MagneticField<DeviceType>& magnetic_field, Smoothing& smoothing, Species<DeviceType>& species){
        GPTLstart("GET_AH_RHO_FF_EXCL_DESTR"); // Timer excludes the destructors which are non-negligible

        // Allocate temporary views
        GPTLstart("GET_AH_RHO_FF_ALLOC_VIEWS");
        int n_input_field_planes = Ah_h.nphi();
        bool gradient_requested = false;
        bool gyroaverage_requested = true;
        species.gyro_avg_matrices.copy_to_device_if_not_resident();
        GetPotentialGradTemp<DeviceType, HostType> tmp(sml, grid, pol_decomp, magnetic_field, pbd0, n_input_field_planes, gradient_requested, gyroaverage_requested, species.gyro_avg_matrices);
        GPTLstop("GET_AH_RHO_FF_ALLOC_VIEWS");

        /****************************************/
        // Calculation of the Hamiltonian vector potential A_h
        /****************************************/
        using FieldArgs = GetPotGradFieldArgs<HostType, HostType, vec2d_if_axisym<PIT_GLOBAL>(), PIT_GLOBAL, TorType::OnePlane, KinType::GyroKin>;
        FieldArgs gpg_field_args(Ah_h.unmanaged());
        gpg_field_args.request_potential(Ah);
        get_field_grad(sml, grid, pol_decomp, magnetic_field, smoothing, gpg_field_args, tmp);

        GPTLstop("GET_AH_RHO_FF_EXCL_DESTR");

        species.gyro_avg_matrices.deallocate_device_matrices_if_not_resident();
    }

    void get_potential_grad(const Simulation<DeviceType>& sml, const Grid<DeviceType>& grid, const DomainDecomposition<DeviceType>& pol_decomp, const MagneticField<DeviceType>& magnetic_field, Smoothing& smoothing, const View<double*, CLayout, HostType>& spitzer_resistivity, Species<DeviceType>& species){
        GPTLstart("GET_POT_GRAD_EXCL_DESTR"); // Timer excludes the destructors which are non-negligible

        // Allocate temporary views
        GPTLstart("GET_POT_GRAD_ALLOC_VIEWS");
        int n_input_field_planes = dpot_h.nphi();
        bool gradient_requested = true;
        bool gyroaverage_requested = true;
#ifdef MULTI_RATE
        gyroaverage_requested = false; // This is a hack. This variable controls the size of the tmp arrays for get_pot_grad,
                                       // but the rho_ff arrays will still exist and the rho=0 values will be correct
#endif
        if(gyroaverage_requested) species.gyro_avg_matrices.copy_to_device_if_not_resident();
        GetPotentialGradTemp<DeviceType, HostType> tmp(sml, grid, pol_decomp, magnetic_field, pbd0, n_input_field_planes, gradient_requested, gyroaverage_requested, species.gyro_avg_matrices);
        GPTLstop("GET_POT_GRAD_ALLOC_VIEWS");

        /****************************************/
        // Calculation of the electric potential E
        /****************************************/
        GPTLstart("GET_POT_GRAD_PHI");

        using FieldArgs = GetPotGradFieldArgs<HostType, HostType, vec2d_if_axisym<PIT_GLOBAL>(), PIT_GLOBAL, TorType::OnePlane, KinType::GyroKin>;
        FieldArgs gpg_field_args(dpot_h.unmanaged(), !(turb_efield));

        // Set up E00
        if(E00_efield){
            gpg_field_args.field00 = Field00<DeviceType>(pot0_h, sml.grad_psitheta);
            gpg_field_args.field00.calculate_gradient(grid, tmp.grad_matrices);
        }
        // Write results of E00
        if(sml.reduced_deltaf && sml.electron_on) gpg_field_args.field00.set_field(grid, tmp.ff, E00_ff_h.f);

        // Set up EMPar
        if(sml.explicit_electromagnetic && (sml.em_mixed_variable && sml.em_pullback_mode==4)) gpg_field_args.request_para_em(Epar_em, sml.em_pullback_dampfac, spitzer_resistivity, AbdH, gyroaverage_requested);

        // Request potential and/or gradient outputs
        // XGC1 algorithm uses potential but XGCa doesnt
        if(!sml.is_XGCa) gpg_field_args.request_potential(pot);
        gpg_field_args.request_gradient(E);

        // Get potentials
        get_field_grad(sml, grid, pol_decomp, magnetic_field, smoothing, gpg_field_args, tmp);

        GPTLstop("GET_POT_GRAD_PHI");

        if(sml.explicit_electromagnetic){
            GPTLstart("GET_POT_GRAD_APAR");

            /****************************************/
            // Calculation of the Hamiltonian vector potential A_h
            /****************************************/
            gpg_field_args = FieldArgs(Ah_h.unmanaged());
            gpg_field_args.request_potential(Ah);
#ifndef MULTI_RATE
            gpg_field_args.request_gradient(dAh);
#endif
            get_field_grad(sml, grid, pol_decomp, magnetic_field, smoothing, gpg_field_args, tmp);

            /****************************************/
            // Calculation of the symplectic vector potential A_s
            /****************************************/
            if (sml.em_mixed_variable){
                gpg_field_args = FieldArgs(As_h.unmanaged());
                gpg_field_args.request_potential(As);
#ifndef MULTI_RATE
                gpg_field_args.request_gradient(dAs);
#endif
                get_field_grad(sml, grid, pol_decomp, magnetic_field, smoothing, gpg_field_args, tmp);
            }

            /****************************************/
            // Calculation of the Hamiltonian vector potential A_h for the control-variate method
            /****************************************/
            if (sml.em_control_variate){
                TIMER("GET_POT_APAR_CV",
                    get_field_Ah_cv_ff(grid, pol_decomp, Ah_cv_h, Ah_cv) );
            }
            GPTLstop("GET_POT_GRAD_APAR");
        }

#ifdef SONIC_GK
        GPTLstart("GET_POT_GRAD_SONIC");
        View<double*,CLayout,HostType> Er_B2_h("Er_B2_h",grid.nnode); 
        View<double*,CLayout,HostType> Ez_B2_h("Ez_B2_h",grid.nnode); 
        View<double*,CLayout,HostType> u2_E_h("u2_E_h",grid.nnode); 

        TIMER("GET_POT_SONIC",
        get_sonic_fields(grid, E, Er_B2_h, Ez_B2_h, u2_E_h) );

        gpg_field_args = FieldArgs(Er_B2_h);
        gpg_field_args.request_gradient(dEr_B2);
        get_field_grad(sml, grid, pol_decomp, magnetic_field, smoothing, gpg_field_args, tmp);

        gpg_field_args = FieldArgs(Ez_B2_h);
        gpg_field_args.request_gradient(dEz_B2);
        get_field_grad(sml, grid, pol_decomp, magnetic_field, smoothing, gpg_field_args, tmp);

        gpg_field_args = FieldArgs(u2_E_h);
        gpg_field_args.request_gradient(du2_E);
        get_field_grad(sml, grid, pol_decomp, magnetic_field, smoothing, gpg_field_args, tmp);

        GPTLstop("GET_POT_GRAD_SONIC");
#endif

        species.gyro_avg_matrices.deallocate_device_matrices_if_not_resident();

        GPTLstop("GET_POT_GRAD_EXCL_DESTR");
    }

    static MemoryPrediction estimate_memory_usage(NLReader::NamelistReader& nlr, const Grid<DeviceType> &grid){
        nlr.use_namelist("sml_param");
        int field_n_rho = nlr.get<int>("sml_grid_nrho", 6);

#ifdef EXPLICIT_EM
        int n_fields = 3 + 1; // E, dAs, dAh, + scalars As,Ah
#else
        int n_fields = 1;
#endif
        double host_field_GB_per_vertex = (53 + 8*grid.nplanes + 11*(field_n_rho + 1))*8*BYTES_TO_GB;
        double cpu_memory_usage = grid.nnode*host_field_GB_per_vertex;

        double device_field_GB_per_vertex = std::max(n_fields*8*grid.nplanes, n_fields*8*(field_n_rho+1))*8*BYTES_TO_GB;
        double gpu_memory_usage = grid.nnode*device_field_GB_per_vertex;

        MemoryPrediction memory_prediction{"Collisions", cpu_memory_usage, gpu_memory_usage};

        return memory_prediction;
    }
};

#endif