XGC1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
velocity_grid.hpp
Go to the documentation of this file.
1 #ifndef VELOCITY_GRID_HPP
2 #define VELOCITY_GRID_HPP
3 
4 #include "vgrid_weights.hpp"
5 #include "lagrange_weights.hpp"
6 
7 // Specifies the dimensions of the velocity grid
8 struct VelocityGrid{
9  int nvp;
10  double vp_max;
11  double dvp;
12 
13  int nmu;
14  double smu_max;
15  double dsmu;
16 
17  double inv_mu0_factor = 1.0/3.0;
18 
19  int nvr;
20  int nvz;
21 
23 
25 
27  // Read in basic velocity grid inputs
28  nlr.use_namelist("f0_param");
29  nvp = nlr.get<int>("f0_nvp", 15);
30  nmu = nlr.get<int>("f0_nmu", 31);
31  vp_max = nlr.get<double>("f0_vp_max", 3.0);
32  smu_max = nlr.get<double>("f0_smu_max", 3.0);
33 
34  // Derived values
35  dvp = vp_max/nvp;
36  dsmu = smu_max/nmu;
37  nvr = nmu+1;
38  nvz = nvp*2+1;
39 
40  // Pseudo inverse
41  pseudo_inv_on = nlr.get<bool>("f0_velocity_interp_use_pseudo_inv", false);
42  if(pseudo_inv_on){
43  element_order = nlr.get<int>("f0_velocity_interp_pseudo_inv_order", 2);
44 #ifdef NO_PETSC
45  exit_XGC("\nError: Cannot run with f0_velocity_interp_use_pseudo_inv = true without using PETSc.\n");
46 #endif
47  if (element_order < 0){
48  exit_XGC("\nError: f0_velocity_interp_pseudo_inv_order must be 0 or larger.\n");
49  }
51  std::string msg = "\nError: Cannot use f0_velocity_interp_pseudo_inv_order larger than ";
52  msg += std::to_string(LAGRANGE_MAX_ORDER);
53  msg += ", if you really want to use that large order then modify LAGRANGE_MAX_ORDER in the code.\n";
54  exit_XGC(msg);
55  }
56  if (element_order > 0){
57  if(nmu % element_order != 0){
58  exit_XGC("\nError: f0_nmu must be divisible by f0_velocity_interp_pseudo_inv_order when f0_velocity_interp_use_pseudo_inv = true.\n");
59  }
60  if(2*nvp % element_order != 0){
61  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");
62  }
63  }
64  }else{
65  element_order = 1;
66  }
67  }
68 
69  // Returns parallel velocity given a grid index
70  KOKKOS_INLINE_FUNCTION double vpar_n(int ivz) const{
71  return (ivz - nvp)*dvp;
72  }
73 
74  // Returns perpendicular velocity given a grid index
75  KOKKOS_INLINE_FUNCTION double vperp_n(int ivr) const{
76  return ivr*dsmu;
77  }
78 
79  // Determine vgrid weights from input smu, vp and vgrid_vol
80  KOKKOS_INLINE_FUNCTION VGridWeights get_weights(double smu, double vp) const{ // double inv_vgrid_vol
81  // Calculate linear weights for velocity grid
82  LinearWeights vr_wt(smu, 1.0/dsmu);
83  LinearWeights vz_wt(vp + vp_max, 1.0/dvp); // Shift so minimum is 0.0
84 
85  //rh For f0_f: conversion from mu-vpar to sqrt(mu)-vpar or vperp-vpar
86  //vr_wt.w[0] *= inv_vgrid_vol * (vr_wt.i==0 ? 2.0 : 1.0);
87  //vr_wt.w[1] *= inv_vgrid_vol * (vr_wt.i==n_vr()-2 ? 2.0 : 1.0);
88 
89  if(vr_wt.i>=nmu || vz_wt.i>=2*nvp || vz_wt.i<0){
90  // Return invalid grid if out of bounds
91  return VGridWeights();
92  }else{
93  return VGridWeights(vr_wt, vz_wt);
94  }
95  }
96 
97  // Divide by dsmu, accounting for difference at edges
98  KOKKOS_INLINE_FUNCTION void weight_by_inv_dsmu(VGridWeights& vgrid_wts) const{
99  double inv_smu0 = 1.0/(dsmu*(vgrid_wts.i_vr==0 ? inv_mu0_factor : vgrid_wts.i_vr));
100  double inv_smu1 = 1.0/(dsmu*(vgrid_wts.i_vr==nmu-1 ? nmu-inv_mu0_factor : vgrid_wts.i_vr+1));
101 
102  vgrid_wts.w_00 *= inv_smu0;
103  vgrid_wts.w_01 *= inv_smu0;
104  vgrid_wts.w_10 *= inv_smu1;
105  vgrid_wts.w_11 *= inv_smu1;
106  }
107 
108  // Multiply by volume, accounting for difference at edges
109  KOKKOS_INLINE_FUNCTION void weight_by_vol(double vol, VGridWeights& vgrid_wts) const{
110  double vol0 = vol*(vgrid_wts.i_vr==0 ? 0.5 : 1.0);
111  double vol1 = vol*(vgrid_wts.i_vr==nmu-1 ? 0.5 : 1.0);
112 
113  vgrid_wts.w_00 *= vol0;
114  vgrid_wts.w_01 *= vol0;
115  vgrid_wts.w_10 *= vol1;
116  vgrid_wts.w_11 *= vol1;
117  }
118 
119  // Divide by volume, accounting for difference at edges
120  KOKKOS_INLINE_FUNCTION void weight_by_inv_vol(double inv_vol, VGridWeights& vgrid_wts) const{
121  double inv_vol0 = inv_vol*(vgrid_wts.i_vr==0 ? 2.0 : 1.0);
122  double inv_vol1 = inv_vol*(vgrid_wts.i_vr==nmu-1 ? 2.0 : 1.0);
123 
124  vgrid_wts.w_00 *= inv_vol0;
125  vgrid_wts.w_01 *= inv_vol0;
126  vgrid_wts.w_10 *= inv_vol1;
127  vgrid_wts.w_11 *= inv_vol1;
128  }
129 };
130 #endif
int element_order
velocity grid finite element order (0 = nearest neighbor), (1 = linear), (2 = quadratic), (3 = cubic), ....
Definition: velocity_grid.hpp:22
double smu_max
max mu
Definition: velocity_grid.hpp:14
int nmu
n points in mu (not including zero)
Definition: velocity_grid.hpp:13
int nvp
n points in parallel velocity (not including zero)
Definition: velocity_grid.hpp:9
T get(const string &param, const T default_val, int val_ind=0)
Definition: NamelistReader.hpp:373
bool pseudo_inv_on
whether pseudo-inverse interpolation is used in velocity space
Definition: velocity_grid.hpp:24
double w_00
Definition: vgrid_weights.hpp:10
Definition: linear_weights.hpp:8
Definition: velocity_grid.hpp:8
KOKKOS_INLINE_FUNCTION void weight_by_inv_vol(double inv_vol, VGridWeights &vgrid_wts) const
Definition: velocity_grid.hpp:120
Definition: NamelistReader.hpp:193
int nvr
full grid size (including zero)
Definition: velocity_grid.hpp:19
double w_01
Definition: vgrid_weights.hpp:11
KOKKOS_INLINE_FUNCTION double vpar_n(int ivz) const
Definition: velocity_grid.hpp:70
double w_11
Definition: vgrid_weights.hpp:13
void use_namelist(const string &namelist)
Definition: NamelistReader.hpp:355
#define LAGRANGE_MAX_ORDER
Definition: lagrange_weights.hpp:18
int i
Definition: linear_weights.hpp:9
double vp_max
max parallel velocity
Definition: velocity_grid.hpp:10
int i_vr
Definition: vgrid_weights.hpp:8
KOKKOS_INLINE_FUNCTION void weight_by_vol(double vol, VGridWeights &vgrid_wts) const
Definition: velocity_grid.hpp:109
void exit_XGC(std::string msg)
Definition: globals.hpp:37
VelocityGrid(NLReader::NamelistReader &nlr)
Definition: velocity_grid.hpp:26
Definition: vgrid_weights.hpp:7
double dsmu
grid spacing in mu
Definition: velocity_grid.hpp:15
double inv_mu0_factor
Set value of lowest mu in grid –&gt; 1/mu0_factor.
Definition: velocity_grid.hpp:17
int nvz
full grid size (including negative and zero)
Definition: velocity_grid.hpp:20
KOKKOS_INLINE_FUNCTION VGridWeights get_weights(double smu, double vp) const
Definition: velocity_grid.hpp:80
double dvp
grid spacing in parallel velocity
Definition: velocity_grid.hpp:11
double w_10
Definition: vgrid_weights.hpp:12
KOKKOS_INLINE_FUNCTION void weight_by_inv_dsmu(VGridWeights &vgrid_wts) const
Definition: velocity_grid.hpp:98
KOKKOS_INLINE_FUNCTION double vperp_n(int ivr) const
Definition: velocity_grid.hpp:75