velocity_grid.hpp

Go to the documentation of this file.

#ifndef VELOCITY_GRID_HPP
#define VELOCITY_GRID_HPP

#include "lagrange_weights.hpp"

// Specifies the dimensions of the velocity grid
struct VelocityGrid{
    int nvp;        
    double vp_max;  
    double dvp;     

    int nmu;        
    double smu_max; 
    double dsmu;    

    double inv_mu0_factor = 1.0/3.0; 

    int nvr;        
    int nvz;        

    int element_order; 

    bool pseudo_inv_on; 

    VelocityGrid(NLReader::NamelistReader& nlr){
        // Read in basic velocity grid inputs
        nlr.use_namelist("f0_param");
        nvp = nlr.get<int>("f0_nvp", 15);
        nmu = nlr.get<int>("f0_nmu", 31);
        vp_max = nlr.get<double>("f0_vp_max", 3.0);
        smu_max = nlr.get<double>("f0_smu_max", 3.0);

        // Derived values
        dvp = vp_max/nvp;
        dsmu = smu_max/nmu;
        nvr = nmu+1;
        nvz = nvp*2+1;

        // Pseudo inverse
        pseudo_inv_on = nlr.get<bool>("f0_velocity_interp_use_pseudo_inv", false);
        if(pseudo_inv_on){
            element_order = nlr.get<int>("f0_velocity_interp_pseudo_inv_order", 2);
#ifdef NO_PETSC
            exit_XGC("\nError: Cannot run with f0_velocity_interp_use_pseudo_inv = true without using PETSc.\n");
#endif
            if (element_order < 0){
                exit_XGC("\nError: f0_velocity_interp_pseudo_inv_order must be 0 or larger.\n");
            }
            else if (element_order > LAGRANGE_MAX_ORDER){
                std::string msg = "\nError: Cannot use f0_velocity_interp_pseudo_inv_order larger than ";
                msg += std::to_string(LAGRANGE_MAX_ORDER);
                msg += ", if you really want to use that large order then modify LAGRANGE_MAX_ORDER in the code.\n";
                exit_XGC(msg);
            }
            if (element_order > 0){
                if(nmu % element_order != 0){
                    exit_XGC("\nError: f0_nmu must be divisible by f0_velocity_interp_pseudo_inv_order when f0_velocity_interp_use_pseudo_inv = true.\n");
                }
                if(2*nvp % element_order != 0){
                    exit_XGC("\nError: 2*f0_nvp must be divisible by f0_velocity_interp_pseudo_inv_order when f0_velocity_interp_use_pseudo_inv = true.\n");
                }
            }
        }else{
            element_order = 1;
        }
    }
};
#endif