field_aligned_local_fields.hpp

Go to the documentation of this file.

#ifndef FIELD_ALIGNED_LOCAL_FIELDS_HPP
#define FIELD_ALIGNED_LOCAL_FIELDS_HPP

#include "globals.hpp"
#include "push_controls.hpp"
#include "grid.hpp"
#include "gyro_radius.hpp"
#include "electric_field.hpp"
#include "grid_field_pack.hpp"


/* This class manages how to obtain the local grid field data during the push. It either:
 *  - a) performs a gather operation when the method fields_at_point is used, or 
 *  - b) performs a gather operation when initialized, stores those values and
 *       returns them when the fields_at_point method is used.
 * Option b) is used in XGC1 for ions, because the ion location in an intermediate RK step
 * might not have a valid associated field since it could be outside of the toroidal sector
 * for which field information is present.
 * The generic declaration uses Option a) */
template<KinType KT, PhiInterpType PIT, typename GFPT>
class FieldAlignedLocalFields{
    // No members

    public:

    // Nothing to do in constructor
    template<class Device>
    KOKKOS_INLINE_FUNCTION FieldAlignedLocalFields(const Grid<Device> &grid, const PushControls &push_controls, const Species<Device> &species,
            const MagneticField<Device> &magnetic_field, const GridFieldPack<Device, GFPT> &gfpack, SimdParticles &part, SimdGridWeights<Order::One, PIT>& grid_wts){}

    /* Gather fields if calculating for electrons */
    template<class Device>
    KOKKOS_INLINE_FUNCTION void fields_at_point(const GridFieldPack<Device, GFPT> fields, const PushControls &push_controls, const Grid<Device> &grid,
            const SimdVector &B, const SimdVector2D &gradpsi, const SimdGridWeights<Order::One, PIT>& grid_wts,
            SimdGyroWeights<KT>& rho_wts, LocalFields<GFPT>& fld) const{
        fields.fields_at_point(push_controls, grid, B,gradpsi,grid_wts,rho_wts,  fld);
    }
};

/* If using XGC1, then the class specialization on KinType GyroKin stores the field aligned
 * local field calculated in the constructor, and converts it to cylindrical coordinates in
 * fields_at_point. */
template<typename GFPT>
class FieldAlignedLocalFields<GyroKin, PhiInterpType::Planes, GFPT>{
    // Simd vectors are hard-coded to (r, z, phi), but here they represent field-aligned coordinates instead
    // "phi" --> parallel to B-field
    LocalFields<GFPT> fa_fld;

    /* Convert stored electric/magnetic fields to field-aligned coordinates from cylindrical coordinates */
    KOKKOS_INLINE_FUNCTION void convert_cyl_to_mag(const SimdVector &B, const SimdVector2D &gradpsi, const LocalFields<GFPT>& fld){
        const auto& fld_E = fld.template get<Label::E>(); // Alias (not a copy)
#  ifdef EXPLICIT_EM
        const auto& fld_dAh = fld.template get<Label::dAh>(); // Alias (not a copy)
        const auto& fld_dAs = fld.template get<Label::dAs>(); // Alias (not a copy)
#  endif
        auto& E = fa_fld.template get<Label::E>(); // Alias (not a copy)
#  ifdef EXPLICIT_EM
        auto& dAh = fa_fld.template get<Label::dAh>(); // Alias (not a copy)
        auto& dAs = fa_fld.template get<Label::dAs>(); // Alias (not a copy)
#  endif
        for (int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
            double Bmag = B.magnitude(i_simd);
            double Bp = B.poloidal_magnitude(i_simd);
            double grad2 = gradpsi.r[i_simd]*gradpsi.r[i_simd] + gradpsi.z[i_simd]*gradpsi.z[i_simd];

            // get E-field in magnetic field-aligned basis
            E.r[i_simd]   = (fld_E.r[i_simd]  *gradpsi.r[i_simd] + fld_E.z[i_simd]  *gradpsi.z[i_simd])/grad2;
            E.z[i_simd]   = (fld_E.r[i_simd]  *B.r[i_simd]       + fld_E.z[i_simd]  *B.z[i_simd]      )/Bp;
            E.phi[i_simd] = (E.z[i_simd]*Bp       + fld_E.phi[i_simd]*B.phi[i_simd]                   )/Bmag;    // parallel field

#  ifdef EXPLICIT_EM
            // grad(A_h)
            dAh.r[i_simd] = (fld_dAh.r[i_simd]*gradpsi.r[i_simd] + fld_dAh.z[i_simd]*gradpsi.z[i_simd])/grad2;
            dAh.z[i_simd] = (fld_dAh.r[i_simd]*B.r[i_simd]     + fld_dAh.z[i_simd]*B.z[i_simd]    )/Bp;
            dAh.phi[i_simd] = (dAh.z[i_simd] *Bp       + fld_dAh.phi[i_simd]*B.phi[i_simd]  )/Bmag;

            // These variables are for the mixed-variables formulation.
            // Make sure that they are initialized to 0 if Hamiltonian
            // formulation is used.
            dAs.r[i_simd] = (fld_dAs.r[i_simd]*gradpsi.r[i_simd] + fld_dAs.z[i_simd]*gradpsi.z[i_simd])/grad2;
            dAs.z[i_simd] = (fld_dAs.r[i_simd]*B.r[i_simd]     + fld_dAs.z[i_simd]*B.z[i_simd]    )/Bp;
            dAs.phi[i_simd] = (dAs.z[i_simd] *Bp       + fld_dAs.phi[i_simd]*B.phi[i_simd]  )/Bmag;
#  endif
        }
    }

    /* Convert stored electric/magnetic fields from field-aligned coordinates to cylindrical coordinates */
    KOKKOS_INLINE_FUNCTION void convert_mag_to_cyl(const SimdVector &B, const SimdVector2D &gradpsi, LocalFields<GFPT>& fld) const{
        auto& fld_E = fld.template get<Label::E>(); // Alias (not a copy)
#  ifdef EXPLICIT_EM
        auto& fld_dAh = fld.template get<Label::dAh>(); // Alias (not a copy)
        auto& fld_dAs = fld.template get<Label::dAs>(); // Alias (not a copy)
#  endif
        const auto& E = fa_fld.template get<Label::E>(); // Alias (not a copy)
#  ifdef EXPLICIT_EM
        const auto& dAh = fa_fld.template get<Label::dAh>(); // Alias (not a copy)
        const auto& dAs = fa_fld.template get<Label::dAs>(); // Alias (not a copy)
#  endif
        for (int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
            double Bmag = B.magnitude(i_simd);
            double Bp = B.poloidal_magnitude(i_simd);
            double Br_norm = B.r[i_simd] / Bp;
            double Bz_norm = B.z[i_simd] / Bp;

            fld_E.r[i_simd]   = E.r[i_simd]*gradpsi.r[i_simd]+E.z[i_simd]*Br_norm;
            fld_E.z[i_simd]   = E.r[i_simd]*gradpsi.z[i_simd]+E.z[i_simd]*Bz_norm;
            fld_E.phi[i_simd] = (E.phi[i_simd]*Bmag - fld_E.r[i_simd]*B.r[i_simd] - fld_E.z[i_simd]*B.z[i_simd])/B.phi[i_simd];

#  ifdef EXPLICIT_EM
            fld_dAh.r[i_simd] = dAh.r[i_simd]*gradpsi.r[i_simd]+dAh.z[i_simd]*Br_norm;
            fld_dAh.z[i_simd] = dAh.r[i_simd]*gradpsi.z[i_simd]+dAh.z[i_simd]*Bz_norm;
            fld_dAh.phi[i_simd] = (dAh.phi[i_simd]*Bmag - fld_dAh.r[i_simd]*B.r[i_simd] - fld_dAh.z[i_simd]*B.z[i_simd])/B.phi[i_simd];

            // The following are for the mixed-variable formulation
            // Make sure that no all variables are initialized
            fld_dAs.r[i_simd] = dAs.r[i_simd]*gradpsi.r[i_simd]+dAs.z[i_simd]*Br_norm;
            fld_dAs.z[i_simd] = dAs.r[i_simd]*gradpsi.z[i_simd]+dAs.z[i_simd]*Bz_norm;
            fld_dAs.phi[i_simd] = (dAs.phi[i_simd]*Bmag - fld_dAs.r[i_simd]*B.r[i_simd] - fld_dAs.z[i_simd]*B.z[i_simd])/B.phi[i_simd];
#  endif
        }
    }

    public:

    /* Gather electric/magnetic fields and convert to field-aligned coordinates for storage */
    template<class Device>
    KOKKOS_INLINE_FUNCTION FieldAlignedLocalFields(const Grid<Device> &grid, const PushControls &push_controls, const Species<Device> &species,
            const MagneticField<Device> &magnetic_field, const GridFieldPack<Device, GFPT> &gfpack, SimdParticles &part, SimdGridWeights<Order::One, PhiInterpType::Planes>& grid_wts){
        // Declare arrays returned from B-field
        SimdVector B;
        SimdVector jacb[3];
        Simd<double> psi;
        SimdVector2D gradpsi;
        Simd<bool> rz_outside;

        // obtain B-field information
        magnetic_field.field(part.ph.v(), B, jacb, psi, gradpsi, rz_outside);
        part.deactivate(rz_outside);

        // Get the current field following grid triangle and p weights
        grid.get_grid_weights(magnetic_field, part.ph.v(), psi, grid_wts);

        // calculate rho weights
        SimdGyroWeights<GyroKin> rho_wts(grid, magnetic_field, species, part.ph.v(), part.ct.mu);

        // get electric field 
        LocalFields<GFPT> fld;
        gfpack.fields_at_point(push_controls, grid, B,gradpsi,grid_wts,rho_wts,  fld);

        // Convert to field-aligned coordinates and store
        convert_cyl_to_mag(B, gradpsi, fld);

        // These fields are not used for ions:
//#if defined(DELTAF_CONV) && !defined(EXPLICIT_EM)
//        fa_fld.template get<Label::E00>() = fld.template get<Label::E00>(); // Should E00 also be field-aligned?
//        fa_fld.template get<Label::ddpotdt>() = fld.template get<Label::ddpotdt>();
//#endif
        // These get stored without any conversion; Consolidate field_aligned_local_fields with local_fields objects to avoid this copy
#ifdef EXPLICIT_EM
        fa_fld.template get<Label::Ah>() = fld.template get<Label::Ah>();
        fa_fld.template get<Label::As>() = fld.template get<Label::As>();
        fa_fld.template get<Label::Epar_em>() = fld.template get<Label::Epar_em>();
#endif
    }

    /* Return the stored field values for XGC1 ions */
    template<class Device>
    KOKKOS_INLINE_FUNCTION void fields_at_point(const GridFieldPack<Device, GFPT> fields, const PushControls &push_controls, const Grid<Device> &grid,
            const SimdVector &B, const SimdVector2D &gradpsi, const SimdGridWeights<Order::One, PhiInterpType::Planes>& grid_wts,
            SimdGyroWeights<GyroKin>& rho_wts, LocalFields<GFPT>& fld) const{
        convert_mag_to_cyl(B, gradpsi, fld);

        // These fields are not used for ions:
//#if defined(DELTAF_CONV) && !defined(EXPLICIT_EM)
//        fld.template get<Label::E00>()=fa_fld.template get<Label::E00>(); // Should E00 also be field-aligned?
//        fld.template get<Label::ddpotdt>() = fa_fld.template get<Label::ddpotdt>();
//#endif
        // These get stored without any conversion; Consolidate field_aligned_local_fields with local_fields objects to avoid this copy
#ifdef EXPLICIT_EM
        fld.template get<Label::Ah>() = fa_fld.template get<Label::Ah>();
        fld.template get<Label::As>() = fa_fld.template get<Label::As>();
        fld.template get<Label::Epar_em>() = fa_fld.template get<Label::Epar_em>();
#endif
    }
};


#endif