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"

// 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, 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),
#  if defined(DELTAF_CONV) && !defined(EXPLICIT_EM)
      E00_ff_h("E00_ff_h",nnode),
      ddpotdt("ddpotdt",nphi,nrho,nnode),
      dpotsave("dpotsave", 2, 2, nnode), // View. dims: RK2, ff, nnode
#  endif
      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)
    {
        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  *************/
    /**************************************************/

    // Resizes allocations for large arrays that have a phi dimension (i.e used for electron push)
    template<PhiInterpType PIT>
    GridFieldPack<Device, PIT> copy_push_fields_to_device(KinType kintype_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
        GridFieldPack<Device, PIT> gfpack(kintype_in, turb_efield);
        //
        // 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);
        //
#ifdef XGC1
        // Copy field data to device (or shallow copy if no device)
        if(kintype_in==DriftKin){
            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 defined(DELTAF_CONV) && !defined(EXPLICIT_EM)
            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,PIT, 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");
        }else{
            GPTLstart("Copy_rho_ff_fields_to_device");
            grid_field_copy(gfpack.template get<Label::E_gyro>(), E);
#  ifdef EXPLICIT_EM
            grid_field_copy(gfpack.template get<Label::dAs_gyro>(), dAs);
            grid_field_copy(gfpack.template get<Label::dAh_gyro>(), dAh);
            grid_field_copy(gfpack.template get<Label::As_gyro>(), As);
            grid_field_copy(gfpack.template get<Label::Ah_gyro>(), Ah);
            grid_field_copy(gfpack.template get<Label::Epar_em_gyro>(), 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);
        gfpack.loop_voltage = GridField<Device, VarType::Scalar,PIT, TorType::OnePlane, KinType::DriftKin>("gfpack_loop_voltage",loop_voltage.nnode());
        ff.cnvt_grid_real2ff(grid,loop_voltage,gfpack.loop_voltage);
#else
        grid_field_copy(gfpack.E_rho, E);
#  ifdef SONIC_GK
        grid_field_copy(gfpack.dEr_B2_rho, dEr_B2);
        grid_field_copy(gfpack.dEz_B2_rho, dEz_B2);
        grid_field_copy(gfpack.du2_E_rho, du2_E);
#  endif
        grid_field_copy(gfpack.loop_voltage,loop_voltage);
#endif
        //grid_field_copy(gfpack.loop_voltage,loop_voltage_h);    
        return gfpack;
    }

    /* Set of functions used to allocate/deallocate gpu memory for Ah during particle scatters */
    // Resizes allocation of Ah_cv
    template<PhiInterpType PIT>
    GridFieldPack<Device, PIT> copy_Ah_cv_to_device_as_Ah(KinType kintype_in, const DomainDecomposition<DeviceType>& pol_decomp, const Grid<DeviceType>& grid) const{
        GridFieldPack<Device, PIT> gfpack(kintype_in, turb_efield);
        // No need to gather Ah_phi_ff, because copy_Ah_cv_to_device_as_Ah is only used in scatter, which already assumes particles are
        // on the correct plane!
#ifdef EXPLICIT_EM
        constexpr int ZERO_PHI = 0;
        gfpack.template get<Label::Ah>().f = View<Field<VarType::Scalar,PIT_GLOBAL>**, CLayout,DeviceType>("gfpack_view",grid.nplanes, grid.nnode);
        Kokkos::deep_copy(my_subview(gfpack.template get<Label::Ah>().f,ZERO_PHI), Ah_cv.f);
#endif
        return gfpack;
    }

    template<PhiInterpType PIT>
    GridFieldPack<Device, PIT> copy_Ah_ff_to_device(KinType kintype_in, const DomainDecomposition<DeviceType>& pol_decomp, const Grid<DeviceType>& grid) const{
        GridFieldPack<Device, PIT> gfpack(kintype_in, turb_efield);
#ifdef EXPLICIT_EM
        if(kintype_in==DriftKin){
            // No need to gather Ah_phi_ff, because copy_Ah_ff_to_device is only used in scatter and pullback, which already assume particles are
            // on the correct plane!
            constexpr int ZERO_GYRO = 0;
            constexpr int ZERO_PHI = 0;
            GridField<HostType,VarType::Scalar,PIT_GLOBAL,TorType::OnePlane,KinType::DriftKin> Ah_tmp("Ah_tmp",1,1,grid.nnode);;
            Kokkos::parallel_for("Ah_ff_rhozero", Kokkos::RangePolicy<HostExSpace>( 0, grid.nnode), [=] (const int inode){
                Ah_tmp(inode) = Ah(inode,ZERO_GYRO);
            });
            gfpack.template get<Label::Ah>().f = View<Field<VarType::Scalar,PIT_GLOBAL>**, CLayout,DeviceType>("gfpack_view",grid.nplanes, grid.nnode);
            Kokkos::deep_copy(my_subview(gfpack.template get<Label::Ah>().f,ZERO_PHI), Ah_tmp.f);
        }else{
            grid_field_copy(gfpack.template get<Label::Ah_gyro>(), Ah);
        }
#endif
        return gfpack;
    }


    // Resizes allocation of dpot_ff and copies to device
    template<PhiInterpType PIT>
    GridFieldPack<Device, PIT> copy_dpot_to_device(KinType kintype_in) const{
        GridFieldPack<Device, PIT> gfpack(kintype_in, turb_efield);
#ifdef XGC1
        grid_field_copy(gfpack.dpot_ff, dpot_ff_h);
#else
        const int PLANE_ZERO=0;
        grid_field_copy(gfpack.dpot, dpot_h.subfield(PLANE_ZERO));
#endif
        return gfpack;
    }

    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