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
#include "perturbed_B_field.hpp"
#include "get_dAdt_guess.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::dAdt_guess>>;
};
 
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
extern "C" void set_rho_ff_pointers(int n, Field<VarType::Scalar,PhiInterpType::None>* dpot, Field<VarType::Scalar,PhiInterpType::Planes>* dpot_ff);
extern "C" void set_rho_pointers(int n, Field<VarType::Scalar,PhiInterpType::None>* dpot);
 
// 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> dAdt_guess;
 
    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; 
 
    // Host View for the adjustable loop voltage current drive
    GridField<HostType,VarType::Scalar,PhiInterpType::None,TorType::OnePlane,KinType::DriftKin> loop_voltage_h;
 
    // Solver boundaries; These are (shallow) copies so that they aren't needed as arguments to get_potential_grid etc. Should be removed in refactor
    Boundary AbdH;
    Boundary pbd0;
 
    // 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) : //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
      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)
    {
        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);
        std::string em_pullback_str = nlr.get<std::string>("sml_em_pullback_method","Electrostatic");
        if(em_pullback_mode==4) em_pullback_str = "IdealMHD";
        PullbackMethod em_pullback_method = (em_pullback_str=="Electrostatic" ? PullbackMethod::Electrostatic : PullbackMethod::IdealMHD);
 
        if (!em_mixed_variable)
            exit_XGC("\nError: sml_em_mixed_variable is explicitly set to .false. which is incompatible with EXPLICIT_EM\n");
 
 
        dpot_es = View<double*,CLayout, HostType>("dpot_es", nnode);
        dAh = GridField<HostType,vec2d_if_axisym<PIT_GLOBAL>(),PIT_GLOBAL,TorType::OnePlane,KinType::GyroKin>("dAh",nphi,nrho,nnode);
        Ah = GridField<HostType,VarType::Scalar,PIT_GLOBAL,TorType::OnePlane,KinType::GyroKin>("Ah",nphi,nrho,nnode);
        Ah_h = GridField<HostType,VarType::Scalar,PhiInterpType::None,TorType::MultiplePlanes,KinType::DriftKin>("Ah_h",4,nnode); // -1:2 in fortran
        Ah_h_backup = GridField<HostType,VarType::Scalar,PhiInterpType::None,TorType::OnePlane,KinType::DriftKin>("Ah_h_backup",nnode); // -1:2 in fortran
        As_h = GridField<HostType,VarType::Scalar,PhiInterpType::None,TorType::MultiplePlanes,KinType::DriftKin>("As_h",4,nnode); // -1:2 in fortran
        As_h_backup = GridField<HostType,VarType::Scalar,PhiInterpType::None,TorType::MultiplePlanes,KinType::DriftKin>("As_h_backup",4,nnode); // -1:2 in fortran
#endif
 
        nlr.use_namelist("performance_param", NLReader::NamelistReader::Optional);
        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.
 
#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_method==PullbackMethod::IdealMHD){
                dAdt_guess = GridField<HostType,VarType::Scalar,PhiInterpType::Planes,TorType::OnePlane,KinType::GyroKin>("dAdt_guess",nphi,nrho,nnode);
            }
        }
#endif
 
        // Set pointers in Fortran to rho arrays
#ifndef NO_FORTRAN_MODULES
#ifdef XGC1
        set_rho_ff_pointers(nnode, dpot_h.data(), dpot_ff_h.data());
#else
        set_rho_pointers(nnode, dpot_h.data());
#endif
#endif
 
        // Initialize loop voltage
        nlr.use_namelist("src_param");
        const double initial_loop_vol = nlr.get<double>("src_loop_voltage",0.0); // Initial value for loop voltage
                                                                                 // Positive loop voltage corresponds to an electron torque opposite to the toroidal
                                                                                 // direction (clockwise looking down on the torus, i.e. in negative Z-direction)
 
        for (int inode=0; inode<nnode; inode++){
             loop_voltage_h(inode) = initial_loop_vol;
        }
    }
 
 
    /**************************************************/
    /************  Copy to GridFieldPack  *************/
    /**************************************************/
    template<MarkerType MT, MagneticFieldMode MFM>
    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, 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>::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);
 
        if constexpr(MT==MarkerType::ReducedDeltaF && MFM==MagneticFieldMode::Electrostatic){
            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
        }
 
        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, 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);
            if constexpr(MFM==MagneticFieldMode::Electromagnetic){
                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, 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, gpg_tmp, my_subview(As_h.f, 2), pot0_em_h, gfpack.template get<Label::dAs>(), gfpack.template get<Label::As>());
            }
        }else{
            pfg.gather_phi_ff_on_device(E.f, gfpack.template get<Label::E>().f);
            if constexpr(MFM==MagneticFieldMode::Electromagnetic){
                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);
            }
        }
        if constexpr(MFM==MagneticFieldMode::Electromagnetic){
            if(dAdt_guess.f.size()>0) pfg.gather_phi_ff_on_device(dAdt_guess.f, gfpack.template get<Label::dAdt_guess>().f);
        }
        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, MagneticFieldMode MFM>
    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, bool near_field=false) const{
        using gfp_push_type = GridFieldPack<Device, typename gfpack_push_type<KinType::GyroKin, PIT_GLOBAL, MT, MFM>::t>;
        gfp_push_type gfpack(turb_efield);
 
        GPTLstart("Copy_rho_ff_fields_to_device");
        grid_field_copy(gfpack.template get<Label::E>(), E);
        if constexpr(MFM==MagneticFieldMode::Electromagnetic){
            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::dAdt_guess>(), dAdt_guess);
        }
        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, bool near_field=false) const{
        if(!sml.explicit_electromagnetic){
            if(markertype_in==MarkerType::ReducedDeltaF){
                if(kintype_in==DriftKin){
                    return copy_push_fields_to_device_drift_kin_XGC1<MarkerType::ReducedDeltaF, MagneticFieldMode::Electrostatic>(sml, pol_decomp, magnetic_field, grid, near_field);
                }else{
                    return copy_push_fields_to_device_gyro_kin_XGC1<MarkerType::ReducedDeltaF, MagneticFieldMode::Electrostatic>(sml, pol_decomp, magnetic_field, grid, near_field);
                }
            }else{
                if(kintype_in==DriftKin){
                    return copy_push_fields_to_device_drift_kin_XGC1<MarkerType::TotalF, MagneticFieldMode::Electrostatic>(sml, pol_decomp, magnetic_field, grid, near_field);
                }else{
                    return copy_push_fields_to_device_gyro_kin_XGC1<MarkerType::TotalF, MagneticFieldMode::Electrostatic>(sml, pol_decomp, magnetic_field, grid, near_field);
                }
            }
        }else{
            if(markertype_in==MarkerType::ReducedDeltaF){
                if(kintype_in==DriftKin){
                    return copy_push_fields_to_device_drift_kin_XGC1<MarkerType::ReducedDeltaF, MagneticFieldMode::Electromagnetic>(sml, pol_decomp, magnetic_field, grid, near_field);
                }else{
                    return copy_push_fields_to_device_gyro_kin_XGC1<MarkerType::ReducedDeltaF, MagneticFieldMode::Electromagnetic>(sml, pol_decomp, magnetic_field, grid, near_field);
                }
            }else{
                if(kintype_in==DriftKin){
                    return copy_push_fields_to_device_drift_kin_XGC1<MarkerType::TotalF, MagneticFieldMode::Electromagnetic>(sml, pol_decomp, magnetic_field, grid, near_field);
                }else{
                    return copy_push_fields_to_device_gyro_kin_XGC1<MarkerType::TotalF, MagneticFieldMode::Electromagnetic>(sml, pol_decomp, magnetic_field, grid, 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, 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, MagneticFieldMode::Electrostatic>::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, 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(grid, 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 PerturbedBField<DeviceType>& perturbed_B_field, 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;
        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);
 
        // 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);
 
        bool need_para_grad = sml.explicit_electromagnetic && sml.em_mixed_variable && sml.em_pullback_method==PullbackMethod::IdealMHD;
        View<double**,CLayout, DeviceType> para_grad;
        if(need_para_grad){
            // dAdt_guess takes parallel gradient as an input, but E is already gyroaveraged and field following
            // So save the simple parallel gradient here
            para_grad = View<double**,CLayout, DeviceType>(NoInit("para_grad"), 2, grid.nnode);
            gpg_field_args.request_para_grad(para_grad);
        }
 
        // Get potentials
        get_field_grad(grid, gpg_field_args, tmp);
 
        GPTLstop("GET_POT_GRAD_PHI");
 
        if(sml.explicit_electromagnetic){
            GPTLstart("GET_EM_PARA");
            if(sml.em_mixed_variable && sml.em_pullback_method==PullbackMethod::IdealMHD){
                get_dAdt_guess(sml, grid, pol_decomp, magnetic_field, smoothing, perturbed_B_field, tmp, para_grad, AbdH, dAdt_guess);
            }
            GPTLstop("GET_EM_PARA");
 
 
            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);
            gpg_field_args.request_gradient(dAh);
            get_field_grad(grid, 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);
                gpg_field_args.request_gradient(dAs);
                get_field_grad(grid, 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(grid, gpg_field_args, tmp);
 
        gpg_field_args = FieldArgs(Ez_B2_h);
        gpg_field_args.request_gradient(dEz_B2);
        get_field_grad(grid, gpg_field_args, tmp);
 
        gpg_field_args = FieldArgs(u2_E_h);
        gpg_field_args.request_gradient(du2_E);
        get_field_grad(grid, 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