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(){}
 
    ElectricField(NLReader::NamelistReader& nlr, const Grid<DeviceType>& grid, int nrho_in, bool es_reduced_deltaf_setup=false);
 
    void write_checkpoint_files(const Grid<DeviceType>& grid, const DomainDecomposition<DeviceType>& pol_decomp, bool loop_voltage_on, bool use_split_weight_scheme, bool em, const XGC_IO_Stream& stream);
 
    void read_checkpoint_files(const Grid<DeviceType>& grid, const DomainDecomposition<DeviceType>& pol_decomp, bool loop_voltage_on, bool use_split_weight_scheme, bool em, const XGC_IO_Stream& stream, const XGC_IO_Stream& em_stream);
 
 
    /**************************************************/
    /************  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