perturbed_B_field.hpp

Go to the documentation of this file.

#ifndef PTB_3DB_HPP
#define PTB_3DB_HPP
#include "NamelistReader.hpp"
#include "grid.hpp"
#include "grid_field.hpp"
#include "smoothing.hpp"
 
extern "C" void get_ptb_3db_arrays(Field<VarType::Vector,PhiInterpType::None>** ptb_3db_bfield_re_vac, Field<VarType::Vector,PhiInterpType::None>** ptb_3db_bfield_im_vac);
 
// Perturbed magnetic field class
template<class Device>
class PerturbedBField {
    public:
 
    // Field variables
    bool rampup_vac;      
    int num_ntor;         
    int rampup_time;      
    int start_time;       
    int mode;             
    bool full_spec_on;    
 
    double rampup_fac;    
    double rampup_fac0;   
 
    // RMP penetration calculation
    int dasdt_opt;          
    int mstep_es;           
    int mstep_em;           
    int es_to_em_dt_ratio;  
 
    bool active; 
 
    // Field arrays
    View<int*,CLayout,Device> ntor;             
    View<Field<VarType::Vector,PhiInterpType::None>**,CLayout,Device> bfield_im_vac;
    View<Field<VarType::Vector,PhiInterpType::None>**,CLayout,Device> bfield_re_vac;
 
    private:
    View<double**,CLayout,HostType> jpar0; // Backup of total parallel current for sliding time average
    View<double***,CLayout,HostType> cden0; // Backup of total charge for sliding time average
    View<double**,CLayout,HostType> As_vac; // Vacuum RMP field
    View<double***,CLayout,HostType> Ah_rmp; // Backup of Ah for RMP penetration; used to approximate dAs/dt
 
    public:
 
    // Constructor
    PerturbedBField(NLReader::NamelistReader& nlr, const DomainDecomposition<DeviceType>& pol_decomp, const MagneticField<DeviceType>& magnetic_field, const Grid<DeviceType>& grid, const Smoothing& smoothing, bool explicit_em);
 
    // Null constructor if ptb_3db is off
    PerturbedBField()
      : mode(0), active(false) {}
 
    void initialize(NLReader::NamelistReader& nlr, const DomainDecomposition<DeviceType>& pol_decomp, const MagneticField<DeviceType>& magnetic_field, const Grid<DeviceType>& grid, const Smoothing& smoothing, bool explicit_em);
 
    void read_3db_xgc_mesh(const std::string& input_dir, const std::string& filename, const Grid<DeviceType>& grid);
 
    void copy_bfield_from_fortran() const{
#ifndef NO_FORTRAN_MODULES
        Field<VarType::Vector,PhiInterpType::None>* ptb_3db_bfield_re_vac;
        Field<VarType::Vector,PhiInterpType::None>* ptb_3db_bfield_im_vac;
        get_ptb_3db_arrays(&ptb_3db_bfield_re_vac, &ptb_3db_bfield_im_vac);
 
        array_deep_copy(bfield_re_vac, ptb_3db_bfield_re_vac);
        array_deep_copy(bfield_im_vac, ptb_3db_bfield_im_vac);
#endif
    }
 
    void copy_to_device_if_on(int gstep){
        if(is_triggered(gstep)){
            copy_bfield_from_fortran();
            active = true;
        }else{
            active = false;
        }
    }
 
    void save_Ah(const GridField<HostType,VarType::Scalar,PhiInterpType::None,TorType::MultiplePlanes,KinType::DriftKin>& Ah) const;
    void replace_B_grad_phi(const View<double**,CLayout,HostType>& view, int gstep, double dt) const;
    void ramp_up_As(int ipc, int gstep, double dt, double dt_full, const View<double*,CLayout,HostType>& As_local);
    void update_rampup_and_write(bool explicit_electromagnetic, int gstep, double time);
    void set_to_As_vac(const GridField<HostType,VarType::Scalar,PhiInterpType::None,TorType::MultiplePlanes,KinType::DriftKin>& As) const;
    void subtract_As_vac_fac(const View<double*,CLayout,HostType>& rhs2) const;
 
    void save_rhs_poisson(int iflag, int ipc, bool em_es_step, const View<double*,CLayout,HostType>& rhs) const;
    void save_rhs_ampere(const Simulation<DeviceType>& sml, const MagneticField<DeviceType>& magnetic_field, const Grid<DeviceType>& grid, const DomainDecomposition<DeviceType>& pol_decomp, Plasma& plasma, const GridField<HostType,VarType::Scalar,PhiInterpType::None,TorType::MultiplePlanes,KinType::DriftKin>& dpot_h, const View<double*,CLayout,HostType>& rhs, int ipc) const;
 
    void write_checkpoint_files(const Grid<DeviceType>& grid, const DomainDecomposition<DeviceType>& pol_decomp, const XGC_IO_Stream& stream) const;
    void read_checkpoint_files(const Grid<DeviceType>& grid, const DomainDecomposition<DeviceType>& pol_decomp, const XGC_IO_Stream& stream, const XGC_IO_Stream& em_stream);
 
    KOKKOS_INLINE_FUNCTION bool is_triggered(int gstep) const;
    KOKKOS_INLINE_FUNCTION bool rampup_is_triggered(int gstep) const;
 
    // Get perturbed b-field
    template<Order OT>
    KOKKOS_INLINE_FUNCTION void get_delta_b(const Grid<Device> &grid, const Simd<double>& phi, const SimdGridWeights<OT, PhiInterpType::None>& grid_wts, SimdVector &tdb) const;
 
};
 
#include "perturbed_B_field.tpp"
#endif