XGCa
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
col_grid.hpp
Go to the documentation of this file.
1 #ifndef COL_GRID_HPP
2 #define COL_GRID_HPP
3 
4 #include <vector>
5 #include <stdio.h>
6 #include <Kokkos_Core.hpp>
7 #include <Kokkos_DualView.hpp>
8 
10 #include "velocity_grid.hpp"
11 #include "vgrid_distribution.hpp"
12 
13 #include "col_grid_matrix.hpp"
14 #include "col_species.hpp"
15 #include "col_vgrids.hpp"
16 
17 extern "C" int get_col_range_setup(int inode);
18 
19 //> XGC's multi-species Fokker-Planck-Landau collision operator
20 // see R. Hager et al., Journal of Computational Physics 315 (2016), pp. 644-660
21 // E. S. Yoon et al., Physics of Plasmas 21 (2014), 032503
22 //
23 // The headers of the subroutines in this module refer to the corresponding parts of Hager et al.
24 // where applicable in order to make the source code easier to understand.
25 
26 
27 // Enum for derivatives of ED array
28 namespace ED{
29 enum{
30  RR=0,
31  RZ,
32  ZZ,
33  ER,
34  EZ,
35  N
36 };
37 }
38 
39 // <SOLVER-related>
40 const int vpic_inner_iter_max=20;
41 
42 // set use_superlu = .true. to use sparse direct SuperLU solver
43 // set use_superlu = .false. to use lapack band solver
44 const bool use_superlu=false;
45 
46 namespace Convergence{
52  Normal=1,
53 
54 };
55 
56 inline bool is_okay(int convergence_status){
57  return (convergence_status!=Convergence::MaxedIter);
58 }
59 
60 class Moments{
61  // Values
62  double w_sum; //> sum of energy over species
63  double p_sum; //> sum of momentum over species
64  double min_p_thres; //> momentum error criterion
65 
66  // Deltas
67  double dw_sum; //> sum of energy change over species (i.e. error in energy)
68  double dp_sum; //> sum of momentum change over species (i.e. error in momentum)
69  double dn_n_max; //> maximum fraction of density changes among species (i.e. maximum relative error in density in all species )
71  bool en_exit_ok;
73 
74  //> Checks the convergence criterion (i.e. conservation of mass, energy and momentum)
75  template<class Device>
76  inline void eval(int isp, int mesh_ind, const CollisionVelocityGrids<Device>& cvg, const CollisionSpecies<Device>& cs,double& dw,double& dp,double& dn_n) const{
77  double vpic_dn = 0.0;
78  double vpic_dw = 0.0;
79  double vpic_dp = 0.0;
80 
81  int gri = cs.s.h_view(isp, mesh_ind).grid_index;
82 
83  for (int index_i = 0; index_i<cvg.mesh_z.extent(2); index_i++){
84  double mesh_z2 = cvg.mesh_z(mesh_ind,gri,index_i)*cvg.mesh_z(mesh_ind,gri,index_i);
85  for (int index_j = 0; index_j<cvg.mesh_r.extent(2); index_j++){
86  double mesh_r2 = cvg.mesh_r(mesh_ind,gri,index_j)*cvg.mesh_r(mesh_ind,gri,index_j);
87  double vpic_dfc = (cs.pdf_np1.h_view(mesh_ind,isp,index_i,index_j) - cs.pdf_n.h_view(mesh_ind,isp,index_i,index_j)) * cvg.vol_h(mesh_ind,gri,index_j);
88 
89  vpic_dn += vpic_dfc;
90  vpic_dw += vpic_dfc*(mesh_z2 + mesh_r2);
91  vpic_dp += vpic_dfc*cvg.mesh_z(mesh_ind,gri,index_i);
92  }
93  }
94 
95  dw = vpic_dw*cs.s.h_view(isp,mesh_ind).mass;
96  dp = vpic_dp*cs.s.h_view(isp,mesh_ind).mass;
97  dn_n = std::abs(vpic_dn/cs.s.h_view(isp,mesh_ind).dens);
98  }
99 
100  public:
101 
102  template<class Device>
103  inline void set(int mesh_ind, const CollisionSpecies<Device>& col_spall){
104  w_sum = 0;
105  p_sum = 0;
106  min_p_thres = 1e99; //just large value to be reinitialized with a small value
107  for (int spi = 0; spi<col_spall.n(); spi++){
108  w_sum += col_spall.s.h_view(spi,mesh_ind).ens;
109  p_sum += col_spall.s.h_view(spi,mesh_ind).mom;
110  min_p_thres = std::min( col_spall.s.h_view(spi,mesh_ind).mass*col_spall.s.h_view(spi,mesh_ind).dens*sqrt(col_spall.s.h_view(spi,mesh_ind).numeric_vth2),
111  min_p_thres );
112  }
113  }
114 
115  template<class Device>
116  inline bool evaluate_conservation(int mesh_ind, const CollisionVelocityGrids<Device>& col_vgrids, const CollisionSpecies<Device>& col_spall){
117  dw_sum = 0.0;
118  dp_sum = 0.0;
119  dn_n_max = 0.0;
120  dens_exit_ok = true;
121 
122  for (int spi = 0; spi<col_spall.n(); spi++){
123  // Re-compute conservation errors after correction
124  double dw,dp,dn_n;
125  eval(spi, mesh_ind, col_vgrids, col_spall,dw,dp,dn_n);// moments evalulation from difference betw. pdf_n and pdf_np1
126  dw_sum += dw;
127  dp_sum += dp;
128  if(dn_n_max<dn_n) dn_n_max = dn_n;
129  dens_exit_ok = dens_exit_ok && (dn_n<1e-10);
130  }
131 
132  en_exit_ok = std::abs(dw_sum/w_sum) <= 1e-7;
133  // bool mom_exit_ok = std::abs(col_dp_sum)/(std::max(std::abs(col_p_sum),min_p_thres)) <= (col_f_vp_max/col_f_nvz);// v_par = max_v_par of simulation domain / vth
134  // ALS : momentum tolerance changed to the tolerance in XGC's former collisions module:
135  // For momentum tolerance use main ion mass and electron density as normalization (when there are more than 2 species).
136  double mass_use = (col_spall.s.h_view(0,mesh_ind).is_electron ? col_spall.s.h_view(1,mesh_ind).mass : col_spall.s.h_view(0,mesh_ind).mass);
137  mom_exit_ok = std::abs(dp_sum)/(std::max(std::abs(p_sum),1e-3*mass_use*col_spall.s.h_view(0,mesh_ind).dens)) <= 1.0e-7;
138 
139  return (dens_exit_ok && en_exit_ok && mom_exit_ok);
140  }
141 };
142 
143 }
144 
145 
146 template<class Device>
147 struct TmpColData{
148 
149  // These views are both on CPU and GPU - dual views
150  Kokkos::View<double***,Device> gammac_spall;
151  Kokkos::DualView<int[4][4],Device> mat_pos_rel_indx; // Shouldn't need to be a dual view
152 
153  Kokkos::View<double***,Device> fhalf;
154  Kokkos::View<double***,Device> dfdr;
155  Kokkos::View<double***,Device> dfdz;
156  Kokkos::View<double****,Device> EDs;
157  Kokkos::View<double******,Device> M_ab;
158  Kokkos::View<double*****,Device> M_s;
159 
160  TmpColData<Device>(int nvr, int nvz, int sp_num, int n_vgrids, int mb_n_nodes)
161  : // Note: layout left
162  gammac_spall("gammac_spall",mb_n_nodes,sp_num,sp_num),
163  fhalf("fhalf",mb_n_nodes,nvz-1,nvr-1),
164  dfdr("dfdr",mb_n_nodes,nvz-1,nvr-1),
165  dfdz("dfdz",mb_n_nodes,nvz-1,nvr-1),
166  EDs("EDs",mb_n_nodes,nvz-1,nvr-1,ED::N),
167  M_s("M_s",mb_n_nodes,n_vgrids, (nvr-1)*(nvz-1), ED::N, nvr-1),
168  M_ab("M_ab",mb_n_nodes,n_vgrids, n_vgrids-1, (nvr-1)*(nvz-1), 3, (nvr-1)*(nvz-1)),
169  mat_pos_rel_indx("mat_pos_rel_indx")
170  {
171  // Set up mat_pos_rel_indx and send to GPU
172  int mat_loc[16] = {4,5,7,8,
173  1,2,4,5,
174  3,4,6,7,
175  0,1,3,4};
176 
177  for (int i=0;i<16;i++)
178  mat_pos_rel_indx.h_view(i/4,i%4)=mat_loc[i];
179 
180 #ifdef USE_GPU
181  Kokkos::deep_copy(mat_pos_rel_indx.d_view, mat_pos_rel_indx.h_view);
182 #endif
183  }
184 };
185 
186 
187 template<class Device>
189 
190  public:
191 
193 
194  double* node_cost; // pointer to array of timings for collisions
195 
196  Kokkos::View<bool*, Kokkos::HostSpace> in_range;
197 
199 
201 
203 
204  //> Set up global variables of the collision operator
208  // These defaults have not been systematically optimized
209 #ifdef USE_GPU
210  int default_batch_size = 64;
211 #else
212  int default_batch_size = 1;//omp_get_max_threads();
213 #endif
214  std::string linalg_backend_str;
215  if(nlr.namelist_present("performance_param")){
216  nlr.use_namelist("performance_param");
217  batch_size = nlr.get<int>("mesh_batch_size", default_batch_size);
218  linalg_backend_str = nlr.get<std::string>("collisions_solver", "lapack");
219  }else{
220  batch_size = default_batch_size;
221  linalg_backend_str = "lapack";
222  }
223 
225  }
226 
227  //> Updates some parameters that may change in the course of the simulation
228  void update_params(bool sml_symmetric_f, int nnode, int first_node, double* node_cost_ptr){
229  // This may not need to be updated
230  node_cost = node_cost_ptr;
231  symmetric_f = sml_symmetric_f;
232 
233  in_range = Kokkos::View<bool*, Kokkos::HostSpace>("in_range",nnode);
234 
235  // nodes outside [sml_inpsi,sml_outpsi] and in private flux region are skipped
236  for (int i=0;i<nnode;i++){
237  int ftr_node = i + first_node;
238  in_range(i) = (get_col_range_setup(ftr_node)==1);
239  }
240  }
241 
242  // Core init
243  void core_init(int isp, int mesh_ind, const CollisionVelocityGrids<Device>& col_vgrids, const CollisionSpecies<Device>& col_spall ) const;
244 
245  // Get maxw_fac
246  static KOKKOS_INLINE_FUNCTION double get_maxw_fac(double mesh_dr, double mesh_r, double numeric_vth2) ;
247 
248  // Get numeric v_thermal equilibrium for a given species against a given grid
249  View<double*,Device> get_numeric_v_thermal_equil(int mb_n_nodes, int spi, int grj, const CollisionVelocityGrids<Device>& col_vgrids, const CollisionSpecies<Device>& col_spall) const;
250 
251  // Delta init
252  void core_delta_init(int mb_n_nodes, int gri, int grj, int spi, CollisionVelocityGrids<Device>& col_vgrids, const CollisionSpecies<Device>& col_spall) const;
253 
254  // LU_matrix_ftn
255  static KOKKOS_INLINE_FUNCTION void LU_matrix_ftn(int mesh_ind, int gri, int grj, int spi, int cell_i,int cell_j,const CollisionVelocityGrids<Device>& col_vgrids,int mprl_col, int mat_pos,double coeff1,double coeff2, const TmpColData<Device>& tcd, const Kokkos::View<int**,Device>& index_map_LU_d,typename Collisions::GridMatrix<Device>::values_array_t LU_values);
256 
257  // LU_matrix
258  void LU_matrix(int mb_n_nodes, int gri, int grj, int spi, const CollisionVelocityGrids<Device>& col_vgrids,const TmpColData<Device>& tcd, Collisions::GridMatrix<Device>* const mtx) const;
259 
260  // the core loop including the picard loop is now inside this function
261  void picard_loop(int vpic_inner_iter_max, const CollisionVelocityGrids<Device>& col_vgrids, CollisionSpecies<Device>& col_spall, TmpColData<Device>& tcd, int mb_n_nodes, Kokkos::View<int*,HostType>& convergence_status, Kokkos::View<int*,HostType>& n_iterations_performed, Kokkos::View<Convergence::Moments*,HostType>& moments) const;
262 
263  // dispatch for E_and_D calc
264  void E_and_D(int mb_n_nodes, int gri, int grj, const CollisionVelocityGrids<Device>& col_vgrids, TmpColData<Device>& tcd) const;
265 
266  // E_and_D_s
267  static KOKKOS_INLINE_FUNCTION void E_and_D_s(int idx, int mb_n_nodes, int nvrm1, int nvzm1, const TmpColData<Device>& tcd, int gri) ;
268 
269  // dispatch for angle_avg calc
270  void angle_avg(int mb_n_nodes, int gri, int grj, CollisionVelocityGrids<Device>& col_vgrids, TmpColData<Device>& tcd) const;
271 
272  // angle_avg_s
273  static KOKKOS_INLINE_FUNCTION void angle_avg_s(int idx, int mb_n_nodes, int nvrm1, int nvzm1, const CollisionVelocityGrids<Device>& col_vgrids, const TmpColData<Device>& tcd, int gri);
274 
275  // E_and_D_ab
276  static KOKKOS_INLINE_FUNCTION void E_and_D_ab(int idx, int mb_n_nodes, int nvrm1, int nvzm1, const CollisionVelocityGrids<Device>& col_vgrids, const TmpColData<Device>& tcd, int gri, int grj);
277 
278  // angle_avg_ab
279  static KOKKOS_INLINE_FUNCTION void angle_avg_ab(int idx, int mb_n_nodes, int nvrm1, int nvzm1, const CollisionVelocityGrids<Device>& col_vgrids, const TmpColData<Device>& tcd, int gri, int grj);
280 
281  void f_df(int mb_n_nodes, const CollisionVelocityGrids<Device>& col_vgrids, const CollisionSpecies<Device>& col_spall,int spi,int grj,TmpColData<Device>& tcd) const;
282 
283  // Main collision algorithm
284  Kokkos::View<int*,HostType> core(CollisionVelocityGrids<Device>& col_vgrids, CollisionSpecies<Device>& col_spall, TmpColData<Device>& tcd, int mb_n_nodes) const;
285 
286 
287  // Carries out collisions
288  void collision(std::vector<Species<DeviceType>> &all_species, const DomainDecomposition<DeviceType>& pol_decomp, const VGridDistribution<HostType>& f0_f, Kokkos::View<double****,Kokkos::LayoutRight, Kokkos::HostSpace>& df0g_tmp, double dt) const;
289 };
290 
291 #include "col_grid.tpp"
292 
293 #endif
double min_p_thres
Definition: col_grid.hpp:64
Definition: col_grid_matrix.hpp:23
double w_sum
Definition: col_grid.hpp:62
void collision(std::vector< Species< DeviceType >> &all_species, const DomainDecomposition< DeviceType > &pol_decomp, const VGridDistribution< HostType > &f0_f, Kokkos::View< double ****, Kokkos::LayoutRight, Kokkos::HostSpace > &df0g_tmp, double dt) const
Definition: col_grid.tpp:25
Definition: col_grid.hpp:50
const bool use_superlu
Definition: col_grid.hpp:44
Kokkos::View< double ******, Device > M_ab
Definition: col_grid.hpp:157
static KOKKOS_INLINE_FUNCTION void LU_matrix_ftn(int mesh_ind, int gri, int grj, int spi, int cell_i, int cell_j, const CollisionVelocityGrids< Device > &col_vgrids, int mprl_col, int mat_pos, double coeff1, double coeff2, const TmpColData< Device > &tcd, const Kokkos::View< int **, Device > &index_map_LU_d, typename Collisions::GridMatrix< Device >::values_array_t LU_values)
Definition: col_grid.tpp:859
int batch_size
Definition: col_grid.hpp:192
Kokkos::View< double ***, Kokkos::LayoutRight, HostType > mesh_z
Definition: col_vgrids.hpp:21
T get(const string &param, const T default_val, int val_ind=0)
Definition: NamelistReader.hpp:353
Kokkos::View< double *****, Device > M_s
Definition: col_grid.hpp:158
Definition: col_grid.hpp:188
Kokkos::DualView< double ****, Kokkos::LayoutRight, Device > pdf_np1
Definition: col_species.hpp:90
void set(int mesh_ind, const CollisionSpecies< Device > &col_spall)
Definition: col_grid.hpp:103
bool symmetric_f
Definition: col_grid.hpp:198
Kokkos::View< double ***, Device > fhalf
Definition: col_grid.hpp:153
Collisions::LinAlgBackend labackend
Definition: col_grid.hpp:200
double p_sum
Definition: col_grid.hpp:63
Kokkos::View< value_type ***, Kokkos::LayoutRight, Device, Kokkos::MemoryTraits< Kokkos::Unmanaged >> values_array_t
Definition: col_grid_matrix.hpp:31
Kokkos::View< double ***, Device > gammac_spall
Definition: col_grid.hpp:150
Definition: NamelistReader.hpp:163
Kokkos::View< double ***, Kokkos::LayoutRight, HostType > vol_h
Definition: col_vgrids.hpp:23
bool en_exit_ok
Definition: col_grid.hpp:71
Definition: col_grid.hpp:30
Kokkos::DualView< CollisionSpeciesScalars **, Device > s
Definition: col_species.hpp:87
double dn_n_max
Definition: col_grid.hpp:69
CollisionGrid()
Definition: col_grid.hpp:202
Kokkos::View< bool *, Kokkos::HostSpace > in_range
Definition: col_grid.hpp:196
Kokkos::View< int *, HostType > core(CollisionVelocityGrids< Device > &col_vgrids, CollisionSpecies< Device > &col_spall, TmpColData< Device > &tcd, int mb_n_nodes) const
Definition: col_grid.tpp:129
ConvergenceStatus
Definition: col_grid.hpp:47
Definition: col_vgrids.hpp:13
void core_delta_init(int mb_n_nodes, int gri, int grj, int spi, CollisionVelocityGrids< Device > &col_vgrids, const CollisionSpecies< Device > &col_spall) const
Definition: col_grid.tpp:798
View< double *, Device > get_numeric_v_thermal_equil(int mb_n_nodes, int spi, int grj, const CollisionVelocityGrids< Device > &col_vgrids, const CollisionSpecies< Device > &col_spall) const
Definition: col_grid.tpp:766
static KOKKOS_INLINE_FUNCTION void E_and_D_ab(int idx, int mb_n_nodes, int nvrm1, int nvzm1, const CollisionVelocityGrids< Device > &col_vgrids, const TmpColData< Device > &tcd, int gri, int grj)
Definition: col_grid.tpp:536
Definition: col_grid.hpp:147
void LU_matrix(int mb_n_nodes, int gri, int grj, int spi, const CollisionVelocityGrids< Device > &col_vgrids, const TmpColData< Device > &tcd, Collisions::GridMatrix< Device > *const mtx) const
Definition: col_grid.tpp:904
const int vpic_inner_iter_max
Definition: col_grid.hpp:40
static KOKKOS_INLINE_FUNCTION void E_and_D_s(int idx, int mb_n_nodes, int nvrm1, int nvzm1, const TmpColData< Device > &tcd, int gri)
Definition: col_grid.tpp:705
void core_init(int isp, int mesh_ind, const CollisionVelocityGrids< Device > &col_vgrids, const CollisionSpecies< Device > &col_spall) const
Definition: col_grid.tpp:377
double dw_sum
Definition: col_grid.hpp:67
Definition: col_grid.hpp:52
idx
Definition: diag_f0_df_port1.hpp:32
Kokkos::DualView< double ****, Kokkos::LayoutRight, Device > pdf_n
Definition: col_species.hpp:89
Definition: col_grid.hpp:51
void use_namelist(const string &namelist)
Definition: NamelistReader.hpp:322
Definition: col_grid.hpp:60
Definition: col_grid.hpp:32
void eval(int isp, int mesh_ind, const CollisionVelocityGrids< Device > &cvg, const CollisionSpecies< Device > &cs, double &dw, double &dp, double &dn_n) const
Definition: col_grid.hpp:76
void E_and_D(int mb_n_nodes, int gri, int grj, const CollisionVelocityGrids< Device > &col_vgrids, TmpColData< Device > &tcd) const
Definition: col_grid.tpp:512
Kokkos::View< double ***, Device > dfdz
Definition: col_grid.hpp:155
Definition: col_grid.hpp:35
static KOKKOS_INLINE_FUNCTION void angle_avg_s(int idx, int mb_n_nodes, int nvrm1, int nvzm1, const CollisionVelocityGrids< Device > &col_vgrids, const TmpColData< Device > &tcd, int gri)
Definition: col_grid.tpp:589
bool dens_exit_ok
Definition: col_grid.hpp:70
bool namelist_present(const string &namelist)
Definition: NamelistReader.hpp:312
Definition: col_grid.hpp:31
bool evaluate_conservation(int mesh_ind, const CollisionVelocityGrids< Device > &col_vgrids, const CollisionSpecies< Device > &col_spall)
Definition: col_grid.hpp:116
Definition: col_grid.hpp:49
int get_col_range_setup(int inode)
Kokkos::View< double ****, Device > EDs
Definition: col_grid.hpp:156
Definition: col_grid.hpp:33
bool is_okay(int convergence_status)
Definition: col_grid.hpp:56
double * node_cost
Definition: col_grid.hpp:194
LinAlgBackend
Definition: col_grid_matrix.hpp:216
int n() const
Definition: col_species.hpp:175
Kokkos::View< double ***, Device > dfdr
Definition: col_grid.hpp:154
static KOKKOS_INLINE_FUNCTION void angle_avg_ab(int idx, int mb_n_nodes, int nvrm1, int nvzm1, const CollisionVelocityGrids< Device > &col_vgrids, const TmpColData< Device > &tcd, int gri, int grj)
Definition: col_grid.tpp:442
Definition: col_grid.hpp:48
bool mom_exit_ok
Definition: col_grid.hpp:72
void picard_loop(int vpic_inner_iter_max, const CollisionVelocityGrids< Device > &col_vgrids, CollisionSpecies< Device > &col_spall, TmpColData< Device > &tcd, int mb_n_nodes, Kokkos::View< int *, HostType > &convergence_status, Kokkos::View< int *, HostType > &n_iterations_performed, Kokkos::View< Convergence::Moments *, HostType > &moments) const
Definition: col_grid.tpp:258
Definition: species.hpp:66
Definition: col_species.hpp:83
Kokkos::DualView< int[4][4], Device > mat_pos_rel_indx
Definition: col_grid.hpp:151
void angle_avg(int mb_n_nodes, int gri, int grj, CollisionVelocityGrids< Device > &col_vgrids, TmpColData< Device > &tcd) const
Definition: col_grid.tpp:424
double dp_sum
Definition: col_grid.hpp:68
void f_df(int mb_n_nodes, const CollisionVelocityGrids< Device > &col_vgrids, const CollisionSpecies< Device > &col_spall, int spi, int grj, TmpColData< Device > &tcd) const
Definition: col_grid.tpp:654
Kokkos::View< double ***, Kokkos::LayoutRight, HostType > mesh_r
Definition: col_vgrids.hpp:20
static KOKKOS_INLINE_FUNCTION double get_maxw_fac(double mesh_dr, double mesh_r, double numeric_vth2)
Definition: col_grid.tpp:750
void update_params(bool sml_symmetric_f, int nnode, int first_node, double *node_cost_ptr)
Definition: col_grid.hpp:228
CollisionGrid(NLReader::NamelistReader &nlr)
Definition: col_grid.hpp:205
Definition: col_grid.hpp:34