velocity_grid.hpp

Go to the documentation of this file.

#ifndef VELOCITY_GRID_HPP
#define VELOCITY_GRID_HPP

#include "vgrid_weights.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() = default; // Default constructor for dummy object.
    
    VelocityGrid(NLReader::NamelistReader& nlr){
        // Read in basic velocity grid inputs
        nlr.use_namelist("f0_param");
        nvp = nlr.get<int>("f0_nvp", 15); // Resolution of velocity space in parallel direcion. Total parallel grid points is 2*f0_nvp+1
        nmu = nlr.get<int>("f0_nmu", 31); // Resolution of velocity space in perpendicular direction. Total perpendicular grid points is f0_nmu+1
        vp_max = nlr.get<double>("f0_vp_max", 3.0); // Maximum parallel velocity of velocity space grid, normalized to thermal velocity
        smu_max = nlr.get<double>("f0_smu_max", 3.0); // Maximum perpendicular velocity of velocity space grid, normalized to thermal velocity

        // 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); // Whether to use pseudo-inverse interpolation between particles and velocity grid
        if(pseudo_inv_on){
            element_order = nlr.get<int>("f0_velocity_interp_pseudo_inv_order", 2); // 1 (linear) or 2 (quadratic) pseudo-inverse interpolation, only used if f0_velocity_interp_use_pseudo_inv = .true.
#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;
        }
    }

    // Returns parallel velocity given a grid index
    KOKKOS_INLINE_FUNCTION double vpar_n(int ivz) const{
        return (ivz - nvp)*dvp;
    }

    // Returns perpendicular velocity given a grid index
    KOKKOS_INLINE_FUNCTION double vperp_n(int ivr) const{
        return ivr*dsmu;
    }

    // Returns factor to handle velocity space volume on edge of velocity grid
    KOKKOS_INLINE_FUNCTION double mu_vol_fac(int ivr) const{
        return (ivr==0 || ivr==nvr-1) ? 0.5 : 1.0;
    }

    // Determine vgrid weights from input smu, vp and vgrid_vol
    KOKKOS_INLINE_FUNCTION VGridWeights get_weights(double smu, double vp) const{ // double inv_vgrid_vol
        // Calculate linear weights for velocity grid
        LinearWeights vr_wt(smu, 1.0/dsmu);
        LinearWeights vz_wt(vp + vp_max, 1.0/dvp); // Shift so minimum is 0.0

        //rh For f0_f: conversion from mu-vpar to sqrt(mu)-vpar or vperp-vpar
        //vr_wt.w[0] *= inv_vgrid_vol * (vr_wt.i==0 ? 2.0 : 1.0);
        //vr_wt.w[1] *= inv_vgrid_vol * (vr_wt.i==n_vr()-2 ? 2.0 : 1.0);

        if(vr_wt.i>=nmu || vz_wt.i>=2*nvp || vz_wt.i<0){
            // Return invalid grid if out of bounds
            return VGridWeights();
        }else{
            return VGridWeights(vr_wt, vz_wt);
        }
    }

    // Divide by dsmu, accounting for difference at edges
    KOKKOS_INLINE_FUNCTION void weight_by_inv_dsmu(VGridWeights& vgrid_wts) const{
        double inv_smu0 = 1.0/(dsmu*(vgrid_wts.i_vr==0     ? inv_mu0_factor     : vgrid_wts.i_vr));
        double inv_smu1 = 1.0/(dsmu*(vgrid_wts.i_vr==nmu-1 ? nmu-inv_mu0_factor : vgrid_wts.i_vr+1));

        vgrid_wts.w_00 *= inv_smu0;
        vgrid_wts.w_01 *= inv_smu0;
        vgrid_wts.w_10 *= inv_smu1;
        vgrid_wts.w_11 *= inv_smu1;
    }

    // Multiply by volume, accounting for difference at edges
    KOKKOS_INLINE_FUNCTION void weight_by_vol(double vol, VGridWeights& vgrid_wts) const{
        double vol0 = vol*(vgrid_wts.i_vr==0     ? 0.5 : 1.0);
        double vol1 = vol*(vgrid_wts.i_vr==nmu-1 ? 0.5 : 1.0);

        vgrid_wts.w_00 *= vol0;
        vgrid_wts.w_01 *= vol0;
        vgrid_wts.w_10 *= vol1;
        vgrid_wts.w_11 *= vol1;
    }

    // Divide by volume, accounting for difference at edges
    KOKKOS_INLINE_FUNCTION void weight_by_inv_vol(double inv_vol, VGridWeights& vgrid_wts) const{
        double inv_vol0 = inv_vol*(vgrid_wts.i_vr==0     ? 2.0 : 1.0);
        double inv_vol1 = inv_vol*(vgrid_wts.i_vr==nmu-1 ? 2.0 : 1.0);

        vgrid_wts.w_00 *= inv_vol0;
        vgrid_wts.w_01 *= inv_vol0;
        vgrid_wts.w_10 *= inv_vol1;
        vgrid_wts.w_11 *= inv_vol1;
    }
};
#endif