globals.hpp

Go to the documentation of this file.

#ifndef GLOBALS_HPP
#define GLOBALS_HPP
#include <limits.h>
#include <string>
#include <cassert>
#include "space_settings.hpp"
#include "array_deep_copy.hpp"
#include "access_add.hpp"
#include "simd.hpp"
#ifdef USE_MPI
#include "my_mpi.hpp"
#endif
#include "constants.hpp"

/* Returns max number of omp threads, or 1 if no omp
 */
inline int get_num_cpu_threads(){
#ifdef USE_OMP
    return omp_get_max_threads();
#else
    return 1;
#endif
}

/* Return true if MPI rank is zero
 * */
inline bool is_rank_zero(){
#ifdef USE_MPI
    return SML_COMM_RANK==0;
#else
    return true;
#endif
}

/* Safely abort and exit the code
 * */
inline void exit_XGC(std::string msg){
    printf("%s",msg.c_str());
    fflush(stdout);
#ifdef USE_MPI
    MPI_Abort(SML_COMM_WORLD, 1);
#else
    exit(1);
#endif
}

/* Asserts condition
 * (Callable from Device)
 * (C++ has no standardized assert with a message)
 * */
KOKKOS_INLINE_FUNCTION void assert_XGC(bool cond, const char* msg){
    if(!cond){
        DEVICE_PRINTF("%s",msg);
        assert(false);
    }
}

/* Check that two integers multiplied together doesn't cause an overflow
 * */
inline bool causes_multiplication_overflow(int a,int b){
    int c = a*b;
    if(b==0) return false; // No overflow
    return !(c/b == a);
}

/* Check that two integers added together doesn't cause an overflow
 * */
inline bool causes_addition_overflow(int a,int b){
    if ( ((b > 0) && (a > INT_MAX - b))    // a + b would overflow
       ||((b < 0) && (a < INT_MIN - b)) ){ // a + b would underflow
        return true;
    } else {
        return false;
    }
}

enum class Order{
    Zero,
    One,
    Two
};

enum SpeciesType {
    ELECTRON = 0,
    MAIN_ION = 1
};

enum KinType{
    DriftKin=0,
    GyroKin=1
};

/* PhiInterpType specifies whether the field uses two values in order to interpolate in the phi direction
 * */
enum class PhiInterpType{
    Planes,
    None
};

// Eventually, nearly the entire code will be templated on PhiInterpType, which is XGC1 vs XGCa
// As an intermediate step, it is convenient to specify the template in one global location
#ifdef XGC1
constexpr PhiInterpType PIT_GLOBAL = PhiInterpType::Planes;
#else
constexpr PhiInterpType PIT_GLOBAL = PhiInterpType::None;
#endif

/* What type of marker weight algorithm to use, reduced delta-f or total-f
 */
enum class MarkerType{
    ReducedDeltaF,
    FullF,
    TotalF,
    None
};
enum class FAnalyticShape{
    Maxwellian,
    SlowingDown,
    None
};
enum class WeightEvoEq{
    Direct,
    PDE,
    None
};

enum class MagneticFieldMode{
    Electromagnetic,
    Electrostatic
};
#ifdef EXPLICIT_EM
constexpr MagneticFieldMode MFM_GLOBAL = MagneticFieldMode::Electromagnetic;
#else
constexpr MagneticFieldMode MFM_GLOBAL = MagneticFieldMode::Electrostatic;
#endif

enum class BFieldSymmetry{
    Tokamak,
    Stellarator
};
#ifdef STELLARATOR
constexpr BFieldSymmetry BFS_GLOBAL = BFieldSymmetry::Stellarator;
#else
constexpr BFieldSymmetry BFS_GLOBAL = BFieldSymmetry::Tokamak;
#endif

/* Whether to use cylindrical limit geometry with periodic boundary function
 */
enum class GeometryType{
    Toroidal,
    CylindricalLimit
};
#ifdef CYLINDRICAL
constexpr GeometryType GEOMETRY = GeometryType::CylindricalLimit;
#else
constexpr GeometryType GEOMETRY = GeometryType::Toroidal;
#endif

template<GeometryType GT>
KOKKOS_INLINE_FUNCTION double geometry_switch(double a, double b);

template<>
KOKKOS_INLINE_FUNCTION double geometry_switch<GeometryType::Toroidal>(double a, double b){
    return a;
}

template<>
KOKKOS_INLINE_FUNCTION double geometry_switch<GeometryType::CylindricalLimit>(double a, double b){
    return b;
}


template<GeometryType GT>
KOKKOS_INLINE_FUNCTION constexpr bool use_toroidal_terms();

template<>
KOKKOS_INLINE_FUNCTION constexpr bool use_toroidal_terms<GeometryType::CylindricalLimit>(){
    return 0;
}

template<>
KOKKOS_INLINE_FUNCTION constexpr bool use_toroidal_terms<GeometryType::Toroidal>(){
    return 1;
}

/***/

enum ParticlePhase {
    PIR = 0, 
    PIZ,     
    PIP,     
    PIRHO,   
    PIW1,    
    PIW2,    
    PTL_NPHASE
};

enum ParticleConsts {
    PIM = 0, 
    PIW0,    
    PIF0,    
#ifdef PTL_G2
    PIG2,    
#endif
    PTL_NCONST
};

// Divide two integers, then round up
KOKKOS_INLINE_FUNCTION int divide_and_round_up(int a, int b){
    return (a+b-1)/b;
}

// Positive modulo. i.e. positive_modulo(-1,3) will return 2, whereas (-1)%3 returns -1.
KOKKOS_INLINE_FUNCTION unsigned positive_modulo( int value, unsigned m) {
    int mod = value % (int)m;
    if (mod < 0) {
        mod += m;
    }
    return mod;
}

/* Returns the offset of an evenly distributed set of objects into even subsets. If the set can be distributed evenly,
 * then this operation is simply the local subset times the result of the division. If not, then the first k
 * subsets will have one extra object, where k is the remainder.
 *
 * n_obj      is the number of objects that need to be divided
 * n_subsets  is the number of ranks the objects are being divided into
 * i_subset   is the index of this particular subset
 *
 */
inline long long int offsets_of_even_distribution(long long int n_obj, long long int n_subsets, long long int i_subset){
    long long int obj_per_subset = n_obj/n_subsets;
    long long int remainder = n_obj%n_subsets;
    if (i_subset<remainder){
        obj_per_subset += 1;
        return i_subset*obj_per_subset;
    }else{
        return i_subset*obj_per_subset + remainder;
    }
}

/* Returns the count of an evenly distributed set of objects into even subsets. If the set can be distributed evenly,
 * then this operation is simply the result of the division. If not, then the first k
 * subsets will have one extra object, where k is the remainder.
 *
 * n_obj      is the number of objects that need to be divided
 * n_subsets  is the number of ranks the objects are being divided into
 * i_subset   is the index of this particular subset
 *
 */
inline long long int counts_of_even_distribution(long long int n_obj, long long int n_subsets, long long int i_subset){
    long long int obj_per_subset = n_obj/n_subsets;
    long long int remainder = n_obj%n_subsets;
    if (i_subset<remainder){
        obj_per_subset += 1;
    }
    return obj_per_subset;
}

// Converts integer to string and adds leading zeros
inline std::string formatted_int2str(int input, int n_digits){
    std::string string_no_leading = std::to_string(input); // no leading zeros

    int initial_length = string_no_leading.length();
    int zeros_to_add = n_digits - std::min(n_digits, initial_length);

    // Add leading zeros
    return std::string(zeros_to_add, '0') + string_no_leading;
}

// is_same_type is equivalent to std::is_same, but can be used in cuda kernels
template<class T, class U> struct is_same_type{static constexpr bool val=false;};
template<class T         > struct is_same_type<T, T>{static constexpr bool val=true;};

// Defined and initialized in my_mpi.cpp, reset in sml.tpp
extern bool global_debug_flag;
extern bool global_perf_barriers_flag;

#endif