diagnostics.hpp

Go to the documentation of this file.

#ifndef DIAGNOSTICS_HPP
#define DIAGNOSTICS_HPP

#include "push_diagnostic.hpp"
#include "heat_diagnostics.hpp"
#include "pseudo_inverse_diag.hpp"
#include "tracer_diag.hpp"
#include "loop_voltage_diag.hpp"
#include "diag_f0_df.hpp"
#include "diag_3d_f0_f.hpp"
#include "diag_f0.hpp"
#include "diag_1d.hpp"
#include "diag_3d.hpp"
#include "particle_stream.hpp"
#include "diag_diffusion_profiles.hpp"
#include <mpi.h>

/* Contains all the diagnostics */
struct Diagnostics{
    PushDiagnostic<DeviceType> push_diag;
    HeatDiagnostics<DeviceType> heat_diag;
    PseudoInvDiagnostics pseudo_inv_diag;
    TracerDiagnostics tracer_diag;
    LoopVolDiagnostics loop_vol_diag;
    ParticleStream particle_stream;
    DiagF0DF diag_f0_df;
    DiagF0 diag_f0;
    Diag3DF0F diag_3d_f0_f;
    Diag1D diag_1d;
    Diag3D diag_3d;
    DiffusionProfilesDiag diag_diff_prof;


    Diagnostics(){}

    Diagnostics(NLReader::NamelistReader& nlr, const Simulation<DeviceType>& sml,
            const MagneticField<DeviceType>& magnetic_field, const Plasma& plasma,
            const DomainDecomposition<DeviceType>& pol_decomp,
            const Grid<DeviceType>& grid, const VelocityGrid& vgrid, int f_source_period,
            bool overwrite_existing){

        nlr.use_namelist("diag_param");
        int diag_1d_period = nlr.get<int>("diag_1d_period",10);
        int diag_3d_on = nlr.get<bool>("diag_3d_on",true);

        bool heat_diag_is_on = nlr.get<bool>("diag_heat_on", false); // Whether to write heat diagnostic.

        // Construct each diagnostic
        if(heat_diag_is_on){
            heat_diag.init(nlr, magnetic_field, grid, plasma.nspecies);
        }

        // Check if some species are deltaf
        bool use_deltaf_diag = plasma.true_for_some_species([&](const Species<DeviceType>& species){return (species.f_analytic_shape!=FAnalyticShape::None && !species.is_adiabatic);});
        push_diag = PushDiagnostic<DeviceType>(nlr, use_deltaf_diag, grid.npsi_surf2, plasma.n_nonadiabatic_species, sml.nthreads);

        diag_f0_df.init(nlr, sml, grid, plasma, f_source_period);

        diag_f0.init(nlr, f_source_period, sml.is_XGCa);

        diag_3d_f0_f.init(nlr, magnetic_field, grid, f_source_period);

#ifdef GB2024
        particle_stream = ParticleStream(nlr, magnetic_field);
#endif
#ifndef NO_PETSC
        if (vgrid.pseudo_inv_on){
            pseudo_inv_diag.init(nlr, f_source_period, "pseudo_inv_diag");
        }
        if (sml.current_drive_on){
            loop_vol_diag.init(nlr, f_source_period,"loop_vol_diag");
        }
#endif
#ifdef ESC_PTL
        tracer_diag.init("tracer_diag");
#endif


        diag_1d.init(nlr, pol_decomp, grid, magnetic_field, plasma, use_deltaf_diag);

        if(diag_3d_on){
            diag_3d.init(nlr, grid, diag_1d_period);
        }

        nlr.use_namelist("diff_param");
        bool diag_diff_prof_on = nlr.get<bool>("diff_update_on",false);
        if (diag_diff_prof_on){
            diag_diff_prof.init(nlr, sml, grid, plasma, f_source_period);
        }

        // Open file streams for diagnostics
        open_streams(overwrite_existing);
    }

    void open_streams(bool overwrite_existing){
        XGC_IO_Mode mode = (overwrite_existing ? XGC_IO_Mode::Write : XGC_IO_Mode::Append);
        if(heat_diag.is_on){
            if(heat_diag.mode==1) heat_diag.open_stream("xgc.heatdiag.bp", mode);
            else heat_diag.open_stream("xgc.heatdiag2.bp", mode);
        }
#ifdef XGCA
        // XGC1 uses a different file for each step. This is kept for back compatibility, but should be consolidated.
        if(diag_f0_df.is_on){
            diag_f0_df.open_stream("xgc.fsourcediag.bp", mode);
        }
#endif
        if(pseudo_inv_diag.is_on){
            pseudo_inv_diag.open_stream("xgc.pseudo_inv.bp", mode);
        }
        if(loop_vol_diag.is_on){
#ifdef USE_MPI
            loop_vol_diag.open_stream("xgc.loop_vol.bp", mode, MPI_COMM_SELF);
#else
            loop_vol_diag.open_stream("xgc.loop_vol.bp", mode);
#endif
        }
        if(tracer_diag.is_on){
            tracer_diag.open_stream("xgc.particle_stream.bp", XGC_IO_Mode::Write);
        }
        if(diag_1d.is_on){
            diag_1d.open_stream("xgc.oneddiag.bp", mode);
        }
#ifdef XGC1
        // Steps within the same file are only used if XGC1 and adios_stage_3d
        // This is kept for back compatibility, but should be consolidated
        if(diag_3d.is_on && diag_3d.adios_stage){
            diag_3d.open_stream("xgc.3d.bp", mode);
        }

        if(diag_f0.is_on && diag_f0.adios_stage){
            diag_f0.open_stream("xgc.f0.bp", mode);
        }

        if(diag_3d_f0_f.is_on && diag_3d_f0_f.adios_stage){
            diag_3d_f0_f.open_stream("xgc.f3d.bp", mode);
        }

        if(diag_diff_prof.is_on){
#ifdef USE_MPI
            diag_diff_prof.open_stream("xgc.diffusion_profiles.bp", mode, MPI_COMM_SELF);
#else
            diag_diff_prof.open_stream("xgc.diffusion_profiles.bp", mode);
#endif
        }
#endif
    }

    void close_streams(){
#ifdef XGCA
        if(diag_f0_df.is_on) diag_f0_df.close_stream();
#endif
        if(heat_diag.is_on) heat_diag.close_stream();
        if(pseudo_inv_diag.is_on) pseudo_inv_diag.close_stream();
        if(loop_vol_diag.is_on) loop_vol_diag.close_stream();
        if(tracer_diag.is_on) tracer_diag.close_stream();
        if(diag_1d.is_on) diag_1d.close_stream();
#ifdef XGC1
        if(diag_3d.is_on && diag_3d.adios_stage) diag_3d.close_stream();
        if(diag_f0.is_on && diag_f0.adios_stage) diag_f0.close_stream();
        if(diag_3d_f0_f.is_on && diag_3d_f0_f.adios_stage) diag_3d_f0_f.close_stream();
#endif
        if(diag_diff_prof.is_on) diag_diff_prof.close_stream();
    }
};

#endif