XGC1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
grid_setup.hpp
Go to the documentation of this file.
1 #ifndef GRID_SETUP_HPP
2 #define GRID_SETUP_HPP
3 
4 #include "globals.hpp"
5 #include "grid.hpp"
6 
7 extern "C" void my_qsort(int left,int right,int* idx,double* psi_surf2_tmp,int npsi_surf2);
8 
9 extern "C" void init_guess_xtable_fort( int* nguess_2, double* grid_guess_max, double* grid_guess_min, double* grid_inv_guess_d, int* ntriangle, int* nnodes, const RZPair* grid_x,
10  Vertex* grid_nd, int* guess_xtable, int* guess_count, int* guess_list_len_in);
11 
12 
13 extern "C" void init_guess_list_fort( int* nguess_2, double* grid_guess_max, double* grid_guess_min, double* grid_inv_guess_d, int* ntriangle, int* nnodes, const RZPair* grid_x,
14  Vertex* grid_nd, int* guess_list, int* guess_xtable, int* guess_count, int* guess_list_len_in);
15 
16 // Provides setup and checks for an analytic example of a grid
17 
18 namespace{
19 
20 // Mapping matrix init
21 void setup_mapping(const View<Vertex*, CLayout, HostType> nd, const View<RZPair*, CLayout, HostType> x, View<VertexMap*, CLayout, HostType>& mapping){
22  Kokkos::parallel_for("setup_mapping", Kokkos::RangePolicy<HostExSpace>(0,mapping.extent(0)), KOKKOS_LAMBDA( const int i ){
23  double dx1r=x[nd[i].vertex[0]-1].r - x[nd[i].vertex[2]-1].r;
24  double dx1z=x[nd[i].vertex[0]-1].z - x[nd[i].vertex[2]-1].z;
25  double dx2r=x[nd[i].vertex[1]-1].r - x[nd[i].vertex[2]-1].r;
26  double dx2z=x[nd[i].vertex[1]-1].z - x[nd[i].vertex[2]-1].z;
27 
28  double det = 1.0/( dx1r*dx2z - dx2r*dx1z);
29 
30  mapping[i].vertex[0].r = dx2z*det;
31  mapping[i].vertex[1].r =-dx2r*det;
32  mapping[i].vertex[0].z =-dx1z*det;
33  mapping[i].vertex[1].z = dx1r*det;
34 
35 
36  mapping[i].vertex[2].r = x[nd[i].vertex[2]-1].r;
37  mapping[i].vertex[2].z = x[nd[i].vertex[2]-1].z;
38  });
39 }
40 
41 // These are the regions designated in the .node file
43  Inside=0,
44  OnWall
45 };
46 
47 // These are the regions used in XGC
48 enum Region{
49  RegionOne=1,
50  RegionTwo,
51  PrivateRegion,
52  Wall=100 // temporary, since 100 is the value used elsewhere
53 };
54 
55 void set_psi(){
56  //for(int i = 0; i<psi.size(); i++){
57  //grid%psi(i)=psi_interpol(grid%x(1,i),grid%x(2,i),0,0);
58 }
59 
60 template<class Device, class GridDevice>
61 void set_rgn(const MagneticField<Device>& bfield, const View<RZPair*,CLayout,GridDevice>& gx_h,
62  const View<double*,CLayout,GridDevice>& psi_h, View<int*, CLayout, GridDevice>& rgn_h){
63 
64  int nnode = gx_h.extent(0);
65  View<RZPair*,CLayout,Device> gx("gx",nnode);
66  View<double*,CLayout,Device> psi("psi",nnode);
67  View<int*,CLayout,Device> rgn("rgn",nnode);
68  Kokkos::deep_copy(gx, gx_h);
69  Kokkos::deep_copy(psi, psi_h);
70  Kokkos::deep_copy(rgn, rgn_h);
71 
72  Kokkos::parallel_for("set_rgn", Kokkos::RangePolicy<typename Device::execution_space>(0,rgn.extent(0)), KOKKOS_LAMBDA( const int i ){
73  if(rgn(i)==Inside){
74  if(bfield.equil.is_in_region_1_or_2(gx(i).r,gx(i).z,psi(i))){
75  if(psi(i) > bfield.equil.xpt_psi - bfield.equil.epsil_psi){
76  rgn(i)=RegionTwo;
77  }else{
78  rgn(i)=RegionOne;
79  }
80  }else{
81  rgn(i)=PrivateRegion;
82  }
83  } else {
84  rgn(i)=Wall;
85  }
86  });
87  Kokkos::fence();
88 
89  Kokkos::deep_copy(gx_h, gx);
90  Kokkos::deep_copy(psi_h, psi);
91  Kokkos::deep_copy(rgn_h, rgn);
92 }
93 
94 // Counts number of wall nodes
95 int get_n_wall_nodes(const View<int*, HostType>& rgn){
96  int count=0;
97  for(int i = 0; i<rgn.size(); i++){
98  if(rgn[i]==Wall) count++;
99  }
100 
101  if (count == 0){
102  // No wall nodes found --> print warning and return
103  if(is_rank_zero()) printf("\ninit_wall: No wall nodes found!\n");
104  }
105 
106  return count;
107 }
108 
109 void setup_wall(const View<int*, CLayout, HostType>& rgn, View<int*, CLayout, HostType>& wall_nodes, View<int*, CLayout, HostType>& node_to_wall){
110  int count=0;
111  for(int i = 0; i<rgn.size(); i++){
112  if(rgn[i]==Wall){
113  wall_nodes[count]=i + 1; // 1-indexed
114  node_to_wall[i]=count + 1; // 1-indexed
115  count++;
116  }else{
117  node_to_wall[i]=0;
118  }
119  }
120 }
121 
122 void order_wall(const View<RZPair*,CLayout,HostType>& gx, const View<int*, CLayout, HostType>& rgn, const View<Vertex*, CLayout, HostType>& nodes, const View<int*,CLayout,DeviceType>& num_t_node_d, const View<int**,CLayout,DeviceType>& tr_node_d, const View<int**,CLayout,HostType>& adj, double z_axis, const View<int*,CLayout,HostType>& wall_nodes, const View<int*,CLayout,HostType>& node_to_wall){
123  // Copy num_t_node and tr_node to host
124  View<int*,CLayout,HostType> num_t_node("num_t_node", num_t_node_d.layout());
125  View<int**,CLayout,HostType> tr_node("tr_node", tr_node_d.layout());
126  Kokkos::deep_copy(num_t_node, num_t_node_d);
127  Kokkos::deep_copy(tr_node, tr_node_d);
128 
129  int nwall = wall_nodes.size();
130 
131  // Allocate temporary array for the re-ordered wall nodes
132  View<int*, CLayout, HostType> ordered_wall_nodes("ordered_wall_nodes", wall_nodes.layout());
133  int lastnode = -1;
134  int thisnode = -1;
135  int nextnode = wall_nodes(0);
136  for(int w=0; w<nwall; w++){
137  if (nextnode==lastnode && is_rank_zero()) printf("WARNING: nextnode and lastnode the same!\n");
138  ordered_wall_nodes(w) = nextnode;
139  lastnode = thisnode;
140  thisnode = nextnode;
141  nextnode = -1;
142 
143  // loop through triangles touching this node
144  for(int j=0; j<num_t_node(thisnode-1); j++){
145  int itr = tr_node(j,thisnode-1);
146 
147  // Find index of wall_node in triangle
148  int iv = nodes(itr).vertex[0]==thisnode ? 0 :
149  nodes(itr).vertex[1]==thisnode ? 1 : 2;
150 
151  // Loop through nodes of this triangle
152  for(int k=0; k<3; k++){
153  int candidate = nodes(itr).vertex[k];
154  if(rgn(candidate-1)==Wall &&
155  candidate!=thisnode && // Skip triangle node of the wall_node in question
156  // adj matrix setup so index 0 -> edge 12, index 1 -> edge 02, index 2 -> edge 01.
157  adj(itr, 3 - (iv + k))==-1 && // Not a boundary if there is an adjacent triangle
158  candidate!=lastnode){ // Skip the lastnode entered (move forward) if after first node
159 
160  // Candidate is the correct next node
161  nextnode = candidate;
162  break;
163  }
164  }
165  if(nextnode!=-1) break; // next node has been found
166  }
167  }
168 
169  if (nextnode != ordered_wall_nodes(0) && is_rank_zero()) printf("WARNING: wall_nodes ordered doesnt form closed surface\n");
170 
171  // First node should be right above the inboard midplane, and last node right below
172  // Find inner-most intersection of z-axis
173  int crossing_segment = -1;
174  double min_r_crossing;
175  bool clockwise;
176  for(int w=0; w<nwall; w++){
177  // Use next point to define segment
178  int w2 = (w+1)%nwall;
179  double z1 = gx(ordered_wall_nodes(w)-1).z;
180  double z2 = gx(ordered_wall_nodes(w2)-1).z;
181  // Check if z-axis is crossed
182  bool z_axis_crossed = (z1>=z_axis && z2<z_axis) || (z1<z_axis && z2>=z_axis);
183  if(z_axis_crossed){
184  double r1 = gx(ordered_wall_nodes(w)-1).r;
185  double r2 = gx(ordered_wall_nodes(w2)-1).r;
186  double dr = r2-r1;
187  double dz = z2-z1; // z_axis_crossed condition guarantees this is not exactly zero
188  double r_crossing = r1 + dr*(z_axis-z1)/dz;
189  // Bound r_crossing to reside on the segment: this could be an issue if dz!= but is so small
190  // there is a numerical issue. (It may also be impossible.)
191  r_crossing = std::max(std::min(r1, r2), r_crossing);
192  r_crossing = std::min(std::max(r1, r2), r_crossing);
193  if(crossing_segment==-1 || r_crossing<min_r_crossing){
194  crossing_segment = w;
195  min_r_crossing = r_crossing;
196  clockwise = (z2>z1);
197  }
198  }
199  }
200 
201  // Edge case: wall does not intersect z-axis. In this case, keep the existing first node.
202  // Because tr_node order is not deterministic, some MPI ranks may choose clockwise/counterclockwise differently
203  // So choose an arbitrary but consistent order now to make sure wall_nodes is the same for every rank
204  if(crossing_segment==-1){
205  if(is_rank_zero()) printf("WARNING: Wall does not appear to intersect z-axis.\n");
206 
207  crossing_segment = 0;
208  int p0 = crossing_segment;
209  int p1 = (crossing_segment+1)%nwall; // Second point in current ordering
210  int pm1 = (nwall+crossing_segment-1)%nwall; // Final point in current ordering
211  double r0 = gx(ordered_wall_nodes(p0)-1).r;
212  double r1 = gx(ordered_wall_nodes(p1)-1).r;
213  double rm1 = gx(ordered_wall_nodes(pm1)-1).r;
214  double z0 = gx(ordered_wall_nodes(p0)-1).z;
215  double z1 = gx(ordered_wall_nodes(p1)-1).z;
216  double zm1 = gx(ordered_wall_nodes(pm1)-1).z;
217 
218  // Determine direction based on looking at the 2 neighboring points. This may put us in the wrong global
219  // order, but the important thing is to be consistent.
220  double det = (r1-r0)*(zm1-z0)-(rm1-r0)*(z1-z0);
221  if(det!=0.0){
222  clockwise = (det<0.0);
223  }else{
224  if(r1!=rm1){
225  // If the three points are colinear, choose direction based on which point has lower r
226  clockwise = (r1<rm1);
227  }else{
228  if(z1!=zm1){
229  clockwise = (z1<zm1);
230  }else{
231  if(nwall<=2){
232  clockwise = true; // This doesn't matter if there are only 1 or 2 wall nodes
233  }else{
234  exit_XGC("\nError: Two wall nodes are the same.\n");
235  }
236  }
237  }
238  }
239  }
240 
241  if(clockwise){
242  // If already clockwise, we just need to cyclically shift the wall
243  int first_point = (crossing_segment+1)%nwall;
244  for(int w=0; w<nwall; w++){
245  wall_nodes(w) = ordered_wall_nodes((first_point+w)%nwall);
246  }
247  }else{
248  // If counterclockwise, we need to reverse the order
249  int first_point = crossing_segment;
250  for(int w=0; w<nwall; w++){
251  wall_nodes(w) = ordered_wall_nodes((nwall+first_point-w)%nwall);
252  }
253  }
254 
255  // Reset node_to_wall
256  Kokkos::deep_copy(node_to_wall, 0);
257  for(int i = 0; i<nwall; i++){
258  node_to_wall(wall_nodes(i)-1) = i+1; // 1-indexed
259  }
260 }
261 
262 
263 template<class Device, class GridDevice>
264 View<double*,CLayout,GridDevice> setup_grid_psi(const MagneticField<Device>& magnetic_field, const View<RZPair*,CLayout,GridDevice>& gx_h){
265  int nnode = gx_h.extent(0);
266  View<RZPair*,CLayout,Device> gx("gx",nnode);
267  Kokkos::deep_copy(gx, gx_h);
268 
269  // Calculate psi at grid points
270  View<double*,CLayout,Device> psi("psi", nnode);
271  Kokkos::parallel_for("psi_calc", Kokkos::RangePolicy<typename Device::execution_space>(0, nnode), KOKKOS_LAMBDA( const int idx ){
272  double psi_local;
273  magnetic_field.psi_bicub.der_zero(gx(idx).r,gx(idx).z, psi_local);
274  psi(idx) = max(1e-70, psi_local);
275  });
276 
277  // Check whether the first vertex is the magnetic axis
278  // If so, set psi exactly to 0. This is needed for 1D
279  // psi search on irregular grid
280  double psi_zero_eps = 5.0e-3;
281  Kokkos::parallel_for("psi_calc_check_first", Kokkos::RangePolicy<typename Device::execution_space>(0, 1), KOKKOS_LAMBDA( const int idx ){
282  double dr = gx(idx).r - magnetic_field.equil.axis_r;
283  double dz = gx(idx).z - magnetic_field.equil.axis_z;
284  if(sqrt(dr*dr + dz*dz) < psi_zero_eps){
285  psi(idx) = 0.0;
286  }
287  });
288 
289  Kokkos::fence();
290 
291  View<double*,CLayout,GridDevice> psi_h("psi", nnode);
292  Kokkos::deep_copy(psi_h, psi);
293 
294  return psi_h;
295 }
296 
297 // Returns volume-weighted average of node psi values along each flux surface
298 template<class Device>
299 View<double*,CLayout,HostType> get_psi_surf(const MagneticField<Device>& magnetic_field,
300  const View<double*,CLayout,HostType>& psi_h,
301  const View<int*, CLayout, HostType>& nnodes_on_surface_h,
302  const View<int**,CLayout,HostType>& surf_idx_h,
303  const View<double*,CLayout,HostType>& node_vol_h){
304 
305  // Copy host views to device
306  View<double*,CLayout,Device> psi(NoInit("psi"), psi_h.layout());
307  View<int*, CLayout,Device> nnodes_on_surface(NoInit("nnodes_on_surface"),nnodes_on_surface_h.layout());
308  View<int**,CLayout,Device> surf_idx(NoInit("surf_idx"),surf_idx_h.layout());
309  View<double*,CLayout,Device> node_vol(NoInit("node_vol"),node_vol_h.layout());
310  Kokkos::deep_copy(psi, psi_h);
311  Kokkos::deep_copy(nnodes_on_surface, nnodes_on_surface_h);
312  Kokkos::deep_copy(node_vol, node_vol_h);
313  Kokkos::deep_copy(surf_idx, surf_idx_h);
314 
315  // Allocate output view
316  View<double*,CLayout,Device> psi_surf(NoInit("psi_surf"), nnodes_on_surface.size());
317 
318  // Loop over surfaces
319  Kokkos::parallel_for("get_psi_surf", Kokkos::RangePolicy<typename Device::execution_space>(0,psi_surf.size()), KOKKOS_LAMBDA( const int i ){
320  double surf_vol = 0.0;
321  double psi_sum = 0.0;
322  for(int j=0; j<nnodes_on_surface(i); j++){
323  int k=surf_idx(i, j)-1; // 1-indexed
324 
325  // Calculate psi_surf by taking a volume-weighted average of psi values along
326  // the nodes of the surface
327  psi_sum += psi(k)*node_vol(k);
328  surf_vol += node_vol(k);
329  }
330  psi_surf(i) = psi_sum/surf_vol;
331  });
332  Kokkos::fence();
333 
334  View<double*,CLayout,HostType> psi_surf_h(NoInit("psi_surf"), psi_surf.layout());
335  Kokkos::deep_copy(psi_surf_h, psi_surf);
336  return psi_surf_h;
337 }
338 
339 /*
340  int nsurfaces;
341  std::vector<int> nnodes_on_surface;
342  int sum_nnodes_on_surfaces;
343  std::vector<int> nodes_of_surfaces;
344  */
345 template<class Device>
346 View<double*,CLayout,HostType> setup_psi_surf2(const MagneticField<Device>& magnetic_field,
347  const View<int*, CLayout, HostType>& nnodes_on_surface_h,
348  const View<int**,CLayout, HostType>& surf_idx_h,
349  const View<double*,CLayout,HostType>& psi_surf_h,
350  const View<RZPair*,CLayout,HostType>& gx_h){
351  const double use_surf_thresh=3.0e-2;
352  int nnode = gx_h.extent(0);
353  View<RZPair*,CLayout,Device> gx("gx",nnode);
354  View<int*, CLayout, Device> nnodes_on_surface(NoInit("nnodes_on_surface"),nnodes_on_surface_h.layout());
355  View<int**,CLayout, Device> surf_idx(NoInit("surf_idx"),surf_idx_h.layout());
356  Kokkos::deep_copy(nnodes_on_surface, nnodes_on_surface_h);
357  Kokkos::deep_copy(surf_idx, surf_idx_h);
358  Kokkos::deep_copy(gx, gx_h);
359 
360  int nsurfaces = nnodes_on_surface.size();
361 
362  int npsi_surf2=0;
363  View<bool*,CLayout,Device> use_surf("use_surf", nsurfaces);
364  Kokkos::parallel_reduce("get_psi_surf2", Kokkos::RangePolicy<typename Device::execution_space>(0,psi_surf_h.size()), KOKKOS_LAMBDA( const int i, int& npsi_surf2_l){
365  double zmin=1e10;
366  for(int j=0; j<nnodes_on_surface(i); j++){
367  int k=surf_idx(i, j)-1; // 1-indexed
368  // Determine which flux surfaces to use for 1D diagnosis
369  // Only use flux surfaces that cross outer midplane
370  if ( (std::abs(gx(k).z-magnetic_field.equil.axis_z) < zmin) && (gx(k).r-(magnetic_field.equil.axis_r - 5.0e-5) >= 0.0) ){
371  zmin=std::abs(gx(k).z -magnetic_field.equil.axis_z);
372  }
373  }
374 
375  if (zmin <= use_surf_thresh){
376  use_surf(i) = true;
377  npsi_surf2_l++;
378  }else{
379  use_surf(i) = false;
380  }
381  }, npsi_surf2);
382  Kokkos::fence();
383 
384  // Bring results to host
385  View<bool*,CLayout,HostType> use_surf_h("use_surf", nsurfaces);
386  Kokkos::deep_copy(use_surf_h, use_surf);
387 
388  // Populate psi_surf2 with desired surfaces from psi_surf
389  View<double*, CLayout, HostType> psi_surf2_h(NoInit("psi_surf2"),npsi_surf2);
390  int j=0;
391  for(int i=0; i<nsurfaces; i++){
392  if (use_surf_h(i)){
393  psi_surf2_h(j) = psi_surf_h(i);
394  j++;
395  }
396  }
397 
398  // Make sure surfaces are in ascending order
399  View<int*, HostType> idx(NoInit("idx"),npsi_surf2);
400  for(int i=0; i<npsi_surf2; i++) idx(i) = i + 1; // idx is one-indexed because it is reordered in Fortran
401  int left=1;
402  int right=npsi_surf2;
403  my_qsort(left,right,idx.data(),psi_surf2_h.data(),npsi_surf2);
404  View<double*, HostType> tmp(NoInit("tmp"),npsi_surf2);
405  for(int i=0; i<npsi_surf2; i++) tmp(i)=psi_surf2_h(idx(i) - 1); // -1 since idx is one-indexed
406  for(int i=0; i<npsi_surf2; i++) psi_surf2_h(i)=tmp(i);
407 
408  return psi_surf2_h;
409 }
410 
411 // Calculate and store magnetic field at nodes
412 template<class Device, class GridDevice>
413 View<double*,CLayout,GridDevice> setup_grid_bfield(const MagneticField<Device>& magnetic_field, const View<RZPair*,CLayout,GridDevice>& gx_h){
414  int nnode = gx_h.extent(0);
415  View<RZPair*,CLayout,Device> gx("gx",nnode);
416  Kokkos::deep_copy(gx, gx_h);
417 
418  // Calculate psi at grid points
419  View<double*,CLayout,Device> bfield("bfield", nnode);
420  Kokkos::parallel_for("node_bfield", Kokkos::RangePolicy<typename Device::execution_space>(0,nnode), KOKKOS_LAMBDA( const int idx ){
421  double br, bz, bphi;
422  magnetic_field.bvec_interpol(gx(idx).r,gx(idx).z,0.0,br,bz,bphi);
423  bfield(idx) = sqrt(br*br + bz*bz + bphi*bphi);
424  });
425  Kokkos::fence();
426 
427  View<double*,CLayout,GridDevice> bfield_h("bfield", nnode);
428  Kokkos::deep_copy(bfield_h, bfield);
429 
430  return bfield_h;
431 }
432 
433 
434 /*void plot_grid(Grid<HostType>& grid){
435  int nx = 80;
436  int ny = 40;
437  double dr = (equil.max_r - equil.min_r)/nx;
438  double dz = (equil.max_z - equil.min_z)/ny;
439  printf("\n");
440  std::vector<int> prevrow(nx+1,-1);
441 
442  for (int j=ny;j>=0;j--){
443  printf("\n");
444  int prevnum = -1;
445  for (int i=0;i<=nx;i++){
446  double r = equil.min_r + i*dr;
447  double z = equil.min_z + j*dz;
448  int num = get_analytic_psi(equil,r,z)*5.0;
449  char mychar = 48+num;
450  if (prevnum!=num || prevrow[i]!=num) printf("%c",mychar);
451  else printf(" ");
452  prevnum = num;
453  prevrow[i] = num;
454  }
455  }
456  printf("\n");
457 }
458 */
459 
460 }
461 #endif
KOKKOS_INLINE_FUNCTION void bvec_interpol(double r, double z, double phi, double &br, double &bz, double &bphi) const
Definition: magnetic_field.tpp:204
bool is_rank_zero()
Definition: globals.hpp:27
void init_guess_xtable_fort(int *nguess_2, double *grid_guess_max, double *grid_guess_min, double *grid_inv_guess_d, int *ntriangle, int *nnodes, const RZPair *grid_x, Vertex *grid_nd, int *guess_xtable, int *guess_count, int *guess_list_len_in)
void my_qsort(int left, int right, int *idx, double *psi_surf2_tmp, int npsi_surf2)
Definition: magnetic_field.hpp:12
Equilibrium equil
The object containing information about the magnetic equilibrium.
Definition: magnetic_field.hpp:32
double axis_z
z coordinate of axis
Definition: equil.hpp:90
Definition: grid_structs.hpp:28
idx
Definition: diag_f0_df_port1.hpp:32
FileRegion
Definition: grid_setup.hpp:42
double axis_r
r coordinate of axis
Definition: equil.hpp:89
void init_guess_list_fort(int *nguess_2, double *grid_guess_max, double *grid_guess_min, double *grid_inv_guess_d, int *ntriangle, int *nnodes, const RZPair *grid_x, Vertex *grid_nd, int *guess_list, int *guess_xtable, int *guess_count, int *guess_list_len_in)
Definition: grid_structs.hpp:24
Region
Definition: grid_setup.hpp:48
void exit_XGC(std::string msg)
Definition: globals.hpp:37
Definition: magnetic_field.F90:1
void parallel_for(const std::string name, int n_ptl, Function func, Option option, HostAoSoA aosoa_h, DeviceAoSoA aosoa_d)
Definition: streamed_parallel_for.hpp:252
Bicub< Device > psi_bicub
The object for interpolating psi (magnetic flux surfaces)
Definition: magnetic_field.hpp:30
double epsil_psi
Not sure?
Definition: equil.hpp:80
Kokkos::ViewAllocateWithoutInitializing NoInit
Definition: space_settings.hpp:68
KOKKOS_INLINE_FUNCTION bool is_in_region_1_or_2(double r, double z, double psi) const
Definition: equil.tpp:82
double xpt_psi
Psi coordinate of 1st X-point.
Definition: equil.hpp:79