sources.hpp

Go to the documentation of this file.

#ifndef SOURCES_HPP
#define SOURCES_HPP

#include "step_trigger.hpp"
#include "monte_carlo_collisions.hpp"
#include "col_grid.hpp"
#include "diffusion.hpp"
#include "radiation.hpp"
#include "pellet.hpp"
#include "thermal_bath.hpp"

enum class CollisionMode{
    None = 0,
    MonteCarloNonConserving, // monte-carlo non-conserving
    NonLinearFokkerPlanckLandau, // fully non-linear multi-species Fokker-Planck-Landau [Hager JCP 315, 644-660 (2016)]
    PETScLandauAMR // PETSc multi-species Landau with adaptive mesh refinement (ongoing implementation)
};

struct Sources{
    private:

    StepTrigger step_trigger; // Triggers when f_source is called

    public:

    int f_source_period; 

    // List of sources (other collisions-related members will be folded into a child class of Source)
    Source neutrals;
    Diffusion diffusion;
    Radiation radiation;
    Source heat_and_torque;
    Source collisions;
    Source current_drive;
    Pellet pellet;
    ThermalBath thermal_bath;

    CollisionMode col_mode; 
    MonteCarloCollider<DeviceType> col_mc; 
    CollisionGrid<DeviceType> col_grid; 

    Sources(NLReader::NamelistReader& nlr, const DomainDecomposition<DeviceType>& pol_decomp, const MagneticField<DeviceType>& magnetic_field, const Grid<DeviceType>& grid, const VelocityGrid& vgrid, Plasma& plasma, bool overwrite_existing_files) {
        nlr.use_namelist("sml_param");
        f_source_period = nlr.get<int>("sml_f_source_period", 1);

        step_trigger = StepTrigger(f_source_period);

        // Set up for each source:
        const int CALL_AT_F_SOURCE = 1; // Set these triggers so that these sources are called every time f_source is called
        bool sml_radiation = nlr.get<bool>("sml_radiation",false); 
                                                                   // data base (must be configured with the rad_param namelist)
        bool sml_neutral = nlr.get<bool>("sml_neutral",false);
        bool sml_source = nlr.get<bool>("sml_source",false);
        bool sml_current_drive = nlr.get<bool>("sml_current_drive_on",false);
        bool sml_thermal_bath_on = nlr.get<bool>("sml_thermal_bath_on",false);
        if (sml_radiation) sml_neutral = true; // radiation overrides neutrals (should throw error instead)

        if(sml_radiation){
            radiation = Radiation(nlr, grid, magnetic_field);
        }

        // Temporary: 'f0_grid' is true if every species is either adiabatic or dynamic_f0
        bool f0_grid = plasma.true_for_all_species([&](const Species<DeviceType>& species){return (species.is_adiabatic || species.dynamic_f0);});

        if(sml_neutral){
            if(!f0_grid) exit_XGC("\nERROR: sml_neutrals can't (yet) be used unless all non-adiabatic species have ptl_dynamic_f0 = .true.\n");

            nlr.use_namelist("neu_param");
            int start_time = nlr.get<int>("neu_start_time",2); 
            neutrals.step_trigger = StepTrigger(CALL_AT_F_SOURCE, start_time);
        }

        if(sml_source){
            if(!f0_grid) exit_XGC("\nERROR: sml_source can't (yet) be used unless all non-adiabatic species have ptl_dynamic_f0 = .true.\n");

            nlr.use_namelist("src_param");
            int start_time = nlr.get<int>("src_start_time",1); 
            heat_and_torque.step_trigger = StepTrigger(CALL_AT_F_SOURCE, start_time);
        }

        if(sml_current_drive){
            nlr.use_namelist("src_param");
            int start_time = nlr.get<int>("src_current_drive_start_time",1);
            current_drive.step_trigger = StepTrigger(CALL_AT_F_SOURCE, start_time);
        }

        if(sml_thermal_bath_on){
            thermal_bath = ThermalBath(nlr);
        }

        nlr.use_namelist("diff_param");
        bool diff_on = nlr.get<bool>("diff_on", false); 
        if(diff_on){
            nlr.use_namelist("sml_param");
            bool electron_on = nlr.get<bool>("sml_electron_on",false);
            // Prevent use of anomalous transport model when electrons are adiabatic
            if (!electron_on) {
                if(is_rank_zero()) printf("\ndiff_on=TRUE is not compatible with sml_electron_on=FALSE");
                exit_XGC("\nThe anomalous diffusion model requires kinetic electrons.\n");
            }
            nlr.use_namelist("diff_param");
            // If diffusion is on, call diffusion constructor
            diffusion = Diffusion(nlr, magnetic_field);

            int start_time = nlr.get<int>("diff_start_time", 1); 
            diffusion.step_trigger = StepTrigger(CALL_AT_F_SOURCE, start_time);
        }

        nlr.use_namelist("src_param");
        bool pellet_on = nlr.get<bool>("src_pellet_on", false); 
        if(pellet_on){
            // If pellet is on, call pellet constructor
            pellet = Pellet(nlr, grid, plasma, magnetic_field);
        }

        // Collisions
        nlr.use_namelist("col_param");
        int col_mode_in = nlr.get<int>("col_mode", 0);
        if(col_mode_in==0) col_mode = CollisionMode::None;
        else if(col_mode_in==1) col_mode = CollisionMode::MonteCarloNonConserving;
        else if(col_mode_in==4) col_mode = CollisionMode::NonLinearFokkerPlanckLandau;
        else if(col_mode_in==5) col_mode = CollisionMode::PETScLandauAMR;
        else{
            if(col_mode_in==2){
                if(is_rank_zero()) printf("\ncol_mode==2 is currently not supported. If you are reviving it, dont forget to address this potential issue: https://github.com/PrincetonUniversity/XGC-Devel/issues/266");
            }
            exit_XGC("\nInvalid col_mode specified. Acceptable options are 0, 1, 4, or 5.\n");
        }

        if(col_mode==CollisionMode::NonLinearFokkerPlanckLandau){
            col_grid = CollisionGrid<DeviceType>(nlr, magnetic_field, grid, overwrite_existing_files);
        }else if(col_mode==CollisionMode::MonteCarloNonConserving){
            col_mc = MonteCarloCollider<DeviceType>(nlr, pol_decomp, magnetic_field, grid, plasma.n_nonadiabatic_species);
        }

        if(col_mode!=CollisionMode::None){
            collisions.step_trigger = StepTrigger(CALL_AT_F_SOURCE, 0); // start time is zero
        }
    }

    // Close streams for diagnostics associated with a source. This could be part of a destructor.
    void close_streams(){
        if(col_mode==CollisionMode::NonLinearFokkerPlanckLandau && col_grid.diag_on){
            col_grid.io_stream->Close();
        }
    }

    /* Returns whether the sources should be executed on this step
     * */
    bool is_triggered(int step){
        return step_trigger.is_triggered(step);
    }
};

#endif