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 "field_weights.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>
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, PIT> &gfpack, SimdParticles &part, Simd<int>& itr){}

    /* Gather fields if calculating for electrons */
    template<class Device>
    KOKKOS_INLINE_FUNCTION void fields_at_point(const GridFieldPack<Device, PIT> fields, const PushControls &push_controls, const Grid<Device> &grid,
            const SimdVector &B, const SimdVector2D &gradpsi, const Simd<double>& fld_phi, const Simd<int>& itr, const SimdGridVec &p,
            SimdGyroRadius<KT>& rho, LocalFields& fld) const{
        fields.fields_at_point(push_controls, grid, B,gradpsi,fld_phi,itr,p,rho,  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<>
class FieldAlignedLocalFields<GyroKin, PhiInterpType::Planes>{
    // Simd vectors are hard-coded to (r, z, phi), but here they represent field-aligned coordinates instead
    // "phi" --> parallel to B-field
    SimdVector E;
#  ifdef EXPLICIT_EM
    SimdVector dAh;
    SimdVector dAs;
#  endif

    // Members that don't get adjusted but need to be stored
#ifdef DELTAF_CONV
    SimdVector E00;
    Simd<double> ddpotdt;
#endif
#ifdef EXPLICIT_EM
    Simd<double> Ah;
    Simd<double> As;
    Simd<double> Epar_em;
#endif

    /* 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& fld){
        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& fld) const{
        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, PhiInterpType::Planes> &gfpack, SimdParticles &part, Simd<int>& itr){
        // Declare arrays returned from B-field
        SimdVector B;
        SimdVector jacb[3];
        Simd<double> psi;
        SimdVector2D gradpsi;
        SimdVector tdb;
        Simd<bool> rz_outside;

        // obtain B-field information
        magnetic_field.field(part.ph.x(), B, jacb, psi, gradpsi, tdb, rz_outside);
        remove_particles(part.gid,rz_outside);

        // Get toroidal angle
        Simd<double> fld_phi = part.ph.phi;
        grid.wedge_modulo_phi(fld_phi);

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

        // calculate rho
        SimdGyroRadius<GyroKin> rho(grid, magnetic_field, species, part.ph.x(), fld_phi, part.ct.mu, itr);

        // get electric field 
        LocalFields fld;
        gfpack.fields_at_point(push_controls, grid, B,gradpsi,fld_phi,itr,p,rho,  fld);

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

        // These get stored without any conversion; Consolidate field_aligned_local_fields with local_fields objects to avoid this copy
#ifdef DELTAF_CONV
        E00.r=fld.E00.r; // Should E00 also be field-aligned?
        E00.z=fld.E00.z;
        E00.phi=fld.E00.phi;
        ddpotdt = fld.ddpotdt;
#endif
#ifdef EXPLICIT_EM
        Ah = fld.Ah;
        As = fld.As;
        Epar_em = fld.Epar_em;
#endif
    }

    /* Return the stored field values for XGC1 ions */
    template<class Device>
    KOKKOS_INLINE_FUNCTION void fields_at_point(const GridFieldPack<Device, PhiInterpType::Planes> fields, const PushControls &push_controls, const Grid<Device> &grid,
            const SimdVector &B, const SimdVector2D &gradpsi, const Simd<double>& fld_phi, const Simd<int>& itr, const SimdGridVec &p,
            SimdGyroRadius<GyroKin>& rho, LocalFields& fld) const{
        convert_mag_to_cyl(B, gradpsi, fld);

        // These get stored without any conversion; Consolidate field_aligned_local_fields with local_fields objects to avoid this copy
#ifdef DELTAF_CONV
        fld.E00.r=E00.r; // Should E00 also be field-aligned?
        fld.E00.z=E00.z;
        fld.E00.phi=E00.phi;
        fld.ddpotdt = ddpotdt;
#endif
#ifdef EXPLICIT_EM
        fld.Ah = Ah;
        fld.As = As;
        fld.Epar_em = Epar_em;
#endif
    }
};


#endif