loop_voltage.hpp

Go to the documentation of this file.

#ifndef LOOP_VOLTAGE_HPP
#define LOOP_VOLTAGE_HPP

#include "grid.hpp"
#include "globals.hpp"


template<class Device,PhiInterpType PIT>
KOKKOS_INLINE_FUNCTION void add_loop_voltage(const PushControls& push_controls, const Grid<Device>& grid, const MagneticField<Device>& magnetic_field, const GridFieldPack<Device,PIT>& gfpack, const Species<Device>& species, const SimdParticles& part, const Simd<double>& B_mag, const Simd<double>& psi, SimdVector& efield, SimdGridWeights<Order::One, PIT>& grid_wts){

    Simd<double> loop_voltage;

//#ifdef XGC1
//    // Must use non-field-following coordinates
//    SimdGridWeights<Order::One, PhiInterpType::None> grid_wts_no_ff;
//    grid.get_grid_weights(magnetic_field, part.ph.v(), grid_wts_no_ff);
//    gfpack.loop_voltage.gather(grid,grid_wts_no_ff,loop_voltage)
//#else
    gfpack.loop_voltage.gather(grid,grid_wts,part.ph.phi,loop_voltage);
//#endif

    for (int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
        double cmrhoB = species.c_m*part.ph.rho[i_simd]*B_mag[i_simd];
        double f0_T_ev = species.eq_temp.value(magnetic_field, psi[i_simd],part.ph.r[i_simd],part.ph.z[i_simd]);

        // Add loop voltage if particle is within specified psi range and isn't beyond the energy threshold
        constexpr double LOOP_VOLTAGE_HIGH_E_LIMIT_FACTOR = 16.0;
        bool loop_voltage_on = (psi[i_simd]<push_controls.loop_voltage_psimax) &&
                               (magnetic_field.equil.is_in_region_1_or_2(part.ph.r[i_simd],part.ph.z[i_simd],psi[i_simd])) &&
                               ((part.ct.mu[i_simd]*2.0*B_mag[i_simd] + cmrhoB*cmrhoB*species.mass)*J_2_EV<(LOOP_VOLTAGE_HIGH_E_LIMIT_FACTOR*f0_T_ev));
        //if(loop_voltage_on) efield.phi[i_simd] += push_controls.loop_voltage/(TWOPI*part.ph.r[i_simd]);
        if(loop_voltage_on) efield.phi[i_simd] += loop_voltage[i_simd]/(TWOPI*part.ph.r[i_simd]);
    }
}

#endif