XGC1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
grid_files.hpp
Go to the documentation of this file.
1 #ifndef GRID_FILES_HPP
2 #define GRID_FILES_HPP
3 
4 #ifdef USE_MPI
5 #include "my_mpi.hpp"
6 #endif
7 #include "file_reader.hpp"
8 #include "magnetic_field.hpp"
9 #include "grid_structs.hpp"
10 #include "grid_setup.hpp"
11 #include "space_settings.hpp"
12 #include "NamelistReader.hpp"
13 #include "check_aif_flx_format.hpp"
14 
15 extern "C" void set_fortran_flx_aif_format(int is_flx_aif_format_int);
16 
17 struct PlaneFiles{
18  // Contents of .node
19  int nnodes;
20  View<RZPair*, HostType> x;
21  View<int*, HostType> rgn;
22 
23  // Contents of .ele
25  View<Vertex*, HostType> nodes;
26 
27  // Contents of .flx.aif
28  View<int*, HostType> i_x;
29 
30  int nsurfaces;
32  View<int*, HostType> i_surf_separatrices;
33  View<int*, HostType> nnodes_on_surface;
36  View<int*, HostType> non_aligned_vert;
37  View<int*, HostType> non_aligned_nsurf;
38  View<int**, CLayout, HostType> non_aligned_surf_idx;
39  View<int**, CLayout, HostType> surf_idx;
40 
41  void convert_to_x_rgn_format(const View<double**, CLayout, HostType>& tmp_nodes, View<RZPair*, HostType>& x, View<int*, HostType>& rgn){
42  Kokkos::parallel_for("transpose_x_and_rgn", Kokkos::RangePolicy<HostExSpace>(0,tmp_nodes.extent(0)), KOKKOS_LAMBDA( const int i ){
43  // tmp_nodes[i,0]; <-- first column is ignored
44  x(i).r = tmp_nodes(i,1);
45  x(i).z = tmp_nodes(i,2);
46  rgn(i) = tmp_nodes(i,3);
47  });
48  }
49 
50  void convert_to_node_format(const View<int**, CLayout, HostType>& tmp_ele, View<Vertex*, HostType>& nodes){
51  Kokkos::parallel_for("transpose_x_and_rgn", Kokkos::RangePolicy<HostExSpace>(0,tmp_ele.extent(0)), KOKKOS_LAMBDA( const int i ){
52  // tmp_ele[i,0]; <-- first column is ignored
53  nodes(i).vertex[0] = tmp_ele(i,1);
54  nodes(i).vertex[1] = tmp_ele(i,2);
55  nodes(i).vertex[2] = tmp_ele(i,3);
56  });
57  }
58 
59  // Default constructor, to use to deallocate the file arrays after use
61 
62 // Should remove this USE_MPI at some point
63 #ifdef USE_MPI
64  // Read the files
65  PlaneFiles(const std::string& node_file, const std::string& element_file, const std::string& flx_aif_file, const MPI_Comm& comm){
66  // MPI setup
67  int my_rank;
68 #ifdef USE_MPI
69  MPI_Comm_rank(comm, &my_rank);
70 #else
71  my_rank = 0;
72 #endif
73  int ROOT_RANK = 0;
74 
75  int max_nodes_on_surface;
76  int n_xpts;
77  int n_separatrices;
78  // Only one MPI rank reads the files, then broadcasts to the others
79  if(my_rank==ROOT_RANK){
80 
81  FileReader file_reader;
82 
84  file_reader.open(node_file);
85  file_reader.read(&nnodes, 1); // Get nnodes from file
86  file_reader.skip(3); // Skip 3 unused values in file
87 
88  // Use dummy array since reformatting is needed
89  View<double**, CLayout, HostType> tmp_nodes(NoInit("tmp_nodes"),nnodes,4); // 4 entries per node: index, r, z, region
90  file_reader.read(tmp_nodes.data(), tmp_nodes.size());
91 
92  // Reformat to x and rgn
93  x = View<RZPair*, HostType>(NoInit("x"),nnodes);
94  rgn = View<int*, HostType>(NoInit("rgn"),nnodes);
95  convert_to_x_rgn_format(tmp_nodes, x, rgn);
96 
98  file_reader.open(element_file);
99  file_reader.read(&ntriangles, 1); // Get nnodes from file
100  file_reader.skip(2); // Skip 2 unused values in file
101 
102  // Use dummy array since reformatting is needed
103  View<int**, CLayout, HostType> tmp_ele(NoInit("tmp_ele"),ntriangles,4); // 4 entries per node: index, r, z, region
104  file_reader.read(tmp_ele.data(), tmp_ele.size());
105 
106  // Reformat since first column is ignored
107  nodes = View<Vertex*, HostType>("nodes",ntriangles);
108  convert_to_node_format(tmp_ele, nodes);
109 
111  if(nnodes==1) exit_XGC("\nGrid must have more than 1 nnode!\n"); // This is assumed in the flx format check
112  bool old_flx_aif_format = is_old_flx_format(flx_aif_file);
113  set_fortran_flx_aif_format(old_flx_aif_format ? 1 : 0);
114 
115  file_reader.open(flx_aif_file);
116 
117  if(old_flx_aif_format){
118  file_reader.read(&nsurfaces, 1); // Get number of surfaces
119  n_xpts = 0;
120  n_separatrices = 0;
121  i_x = View<int*, HostType>("i_x",0); // No xpoints specified
122  i_surf_separatrices = View<int*, HostType>(NoInit("i_surf_separatrices"), 0); // No separatrices specified
123 
124  // These don't get read or used by the old format
125  nsurfaces1 = 0;
126  nsurfaces2 = 0;
127  nsurfaces3a = 0;
128  nsurfaces3b = 0;
129  }else{
130  file_reader.read(&n_xpts, 1); // Get number of x-points
131  i_x = View<int*, HostType>("i_x", n_xpts);
132 
133  // Current flx.aif files list 2 x-point indices, regardless of n_xpts
134  int i_xpt1, i_xpt2;
135  file_reader.read(&i_xpt1, 1); // Get index of 1st x-point
136  file_reader.read(&i_xpt2, 1); // Get index of 2nd x-point
137 
138  // Check for consistency
139  if(n_xpts==1 && i_xpt2 != -1) exit_XGC("Error in flx.aif file: n_xpts=1, but i_xpt2 is specified.\n");
140  if(n_xpts==2 && i_xpt2 <= 0) exit_XGC("Error in flx.aif file: n_xpts=2, but i_xpt2 is not specified.\n");
141 
142  if(n_xpts>=1) i_x(0) = i_xpt1 - 1; // switch to zero-index
143  if(n_xpts>=2) i_x(1) = i_xpt2 - 1;
144 
145  // Next are the nsurface variables
146  file_reader.read(&nsurfaces1, 1); // Get number of surfaces
147  file_reader.read(&nsurfaces2, 1); // Get number of surfaces
148  file_reader.read(&nsurfaces3a, 1); // Get number of surfaces
149  file_reader.read(&nsurfaces3b, 1); // Get number of surfaces
150 
151  // nsurfaces is actually the sum of these 4 inputs
153 
154  // Get surface index of separatrices
155  // Hard-coded to 2, but generalize the read just in case
156  n_separatrices = 2;
157  i_surf_separatrices = View<int*, HostType>(NoInit("i_surf_separatrices"), n_separatrices);
158  file_reader.read(&i_surf_separatrices(0), 1);
159  file_reader.read(&i_surf_separatrices(1), 1);
160  }
161 
162  // Read in nnodes on each surface
163  nnodes_on_surface = View<int*, HostType>(NoInit("nnodes_on_surface"), nsurfaces);
164  file_reader.read(nnodes_on_surface.data(), nnodes_on_surface.size());
165 
166  // Loop through surfaces to find max number of nodes in a surface
167  max_nodes_on_surface = 0;
168  for(int i=0; i<nnodes_on_surface.size(); i++) max_nodes_on_surface = std::max(nnodes_on_surface(i), max_nodes_on_surface);
169 
170  // Allocate surf_idx and initialize to 0
171  surf_idx = View<int**, CLayout, HostType>("surf_idx",nsurfaces, max_nodes_on_surface);
172 
173  // Read in values
174  for(int i=0; i<nnodes_on_surface.size(); i++){
175  file_reader.read(&surf_idx(i,0), nnodes_on_surface(i));
176  }
177 
178  // There is a "-1" in the file here to mark the end of the section
179  int should_be_negative_one;
180  file_reader.read(&should_be_negative_one, 1);
181  if(should_be_negative_one != -1){
182  printf("\nExpected -1 at end of section in surface file, encountered %d instead\n", should_be_negative_one);
183  exit_XGC("Error reading surface file.\n");
184  }
185 
186  file_reader.read(&num_non_aligned, 1);
187  printf(" - reading num_non_aligned: %d",num_non_aligned);
188 
189  // Read the non-aligned vertex data
190  non_aligned_vert = View<int*, HostType>(NoInit("non_aligned_vert"), num_non_aligned);
191  non_aligned_nsurf = View<int*, HostType>(NoInit("non_aligned_nsurf"), num_non_aligned);
192  // Initialize to 0
193  non_aligned_surf_idx = View<int**, CLayout, HostType>("non_aligned_surf_idx", num_non_aligned, 4);
194  for(int i=0; i<num_non_aligned; i++){
195  file_reader.read(&non_aligned_vert(i), 1);
196  if(old_flx_aif_format){
197  file_reader.read(&non_aligned_nsurf(i), 1);
198  }else{
199  // New format does not contain non_aligned_nsurf
200  non_aligned_nsurf(i)=2;
201  }
202  file_reader.read(&non_aligned_surf_idx(i,0), non_aligned_nsurf(i));
203  }
204 
205  // Another "-1" at the end of the file
206  file_reader.read(&should_be_negative_one, 1);
207  if(should_be_negative_one != -1){
208  printf("\nExpected -1 at end of the surface file, encountered %d instead\n", should_be_negative_one);
209  exit_XGC("Error reading surface file.\n");
210  }
211  }
212 
213  /*** Broadcast grid contents to other MPI ranks ***/
214 
215  // First broadcast scalars so we know the sizes of the arrays
216  int n_scalars = 11;
217  View<int*, HostType> scalars("scalars", n_scalars);
218  scalars(0) = nnodes;
219  scalars(1) = ntriangles;
220  scalars(2) = nsurfaces;
221  scalars(3) = max_nodes_on_surface;
222  scalars(4) = num_non_aligned;
223  scalars(5) = n_xpts;
224  scalars(6) = n_separatrices;
225  scalars(7) = nsurfaces1;
226  scalars(8) = nsurfaces2;
227  scalars(9) = nsurfaces3a;
228  scalars(10) = nsurfaces3b;
229 
230  // Broadcast file size to other processes
231 #ifdef USE_MPI
232  MPI_Bcast(scalars.data(), scalars.size(), MPI_INT, ROOT_RANK, comm);
233 #endif
234 
235  if(my_rank!=ROOT_RANK){
236  // Set scalars on non-root ranks
237  nnodes = scalars(0);
238  ntriangles = scalars(1);
239  nsurfaces = scalars(2);
240  max_nodes_on_surface = scalars(3);
241  num_non_aligned = scalars(4);
242  n_xpts = scalars(5);
243  n_separatrices = scalars(6);
244  nsurfaces1 = scalars(7);
245  nsurfaces2 = scalars(8);
246  nsurfaces3a = scalars(9);
247  nsurfaces3b = scalars(10);
248 
250  // node and ele views
251  x = View<RZPair*, HostType>(NoInit("x"),nnodes);
252  rgn = View<int*, HostType>(NoInit("rgn"),nnodes);
253  nodes = View<Vertex*, HostType>("nodes",ntriangles);
254 
255  // flx.aif views
256  i_x = View<int*, HostType>("i_x", n_xpts);
257  i_surf_separatrices = View<int*, HostType>(NoInit("i_surf_separatrices"), n_separatrices);
258  non_aligned_vert = View<int*, HostType>(NoInit("non_aligned_vert"), num_non_aligned);
259  non_aligned_nsurf = View<int*, HostType>(NoInit("non_aligned_nsurf"), num_non_aligned);
260  non_aligned_surf_idx = View<int**, CLayout, HostType>(NoInit("non_aligned_surf_idx"), num_non_aligned, 4);
261  nnodes_on_surface = View<int*, HostType>(NoInit("nnodes_on_surface"), nsurfaces);
262  surf_idx = View<int**, CLayout, HostType>(NoInit("surf_idx"),nsurfaces, max_nodes_on_surface);
263  }
264 
265  // Broadcast grid arrays
266 #ifdef USE_MPI
267  MPI_Bcast( x.data(), x.size()*2, MPI_DOUBLE, ROOT_RANK, comm);
268  MPI_Bcast( rgn.data(), rgn.size(), MPI_INT, ROOT_RANK, comm);
269  MPI_Bcast( nodes.data(), nodes.size()*3, MPI_INT, ROOT_RANK, comm);
270 
271  if(i_x.size()>0) MPI_Bcast( i_x.data(), i_x.size(), MPI_INT, ROOT_RANK, comm);
272  if(i_surf_separatrices.size()>0) MPI_Bcast( i_surf_separatrices.data(), i_surf_separatrices.size(), MPI_INT, ROOT_RANK, comm);
273  MPI_Bcast( non_aligned_vert.data(), non_aligned_vert.size(), MPI_INT, ROOT_RANK, comm);
274  MPI_Bcast( non_aligned_nsurf.data(), non_aligned_nsurf.size(), MPI_INT, ROOT_RANK, comm);
275  MPI_Bcast( non_aligned_surf_idx.data(), non_aligned_surf_idx.size(), MPI_INT, ROOT_RANK, comm);
276  MPI_Bcast( nnodes_on_surface.data(), nnodes_on_surface.size(), MPI_INT, ROOT_RANK, comm);
277  MPI_Bcast( surf_idx.data(), surf_idx.size(), MPI_INT, ROOT_RANK, comm);
278 #endif
279  }
280 #endif
281 
282  // Options: hex vs circular
286  };
287 
297  void construct_hex_grid(int nshells, double raxis, double zaxis, double rscale, double zscale) {
298  // First, build construction of what .node and .ele files would look like
299 
300  // Build hexagonal sample grid of concentric shells
301  // Minimum 1 shell, which produces 7 nodes and 6 triangles
302 
303  // Grid construction uses this enum for directions
304  enum { NE = 0, N, NW, SW, S, SE };
305 
306  // Node indices of each surface
307  int max_nodes_on_surface = 6*nshells;
308  surf_idx = View<int**, CLayout, HostType>(NoInit("surf_idx"),nsurfaces, max_nodes_on_surface);
309 
310  // Work from center outward
311  x(0).r = raxis;
312  x(0).z = zaxis;
313  rgn(0) = 0; // not wall
314  nnodes_on_surface(0) = 1;
315  surf_idx(0, 0) = 1; // Since surf_idx is 1-indexed
316  int inode = 1;
317  int itri = 0;
318  for(int ish = 1; ish<=nshells; ish++){
319  int nodes_in_shell = 6*ish;
320  int first_inner_shell_node = ish==1 ? 0 : inode - 6*(ish-1);
321  int first_node = inode;
322  int inner_shell_node = first_inner_shell_node;
323  // Starting position:
324  double r_track = raxis + rscale*ish;
325  double z_track = zaxis;
326  for (int inodesh = 0; inodesh < nodes_in_shell; inodesh++){
327  x(inode).r = r_track;
328  x(inode).z = z_track;
329  rgn(inode) = (ish==nshells? OnWall : Inside); // last shell is wall
330 
331  surf_idx(ish, inodesh) = inode+1; // Since surf_idx is 1-indexed
332 
333  if (inodesh>0){
334  nodes(itri).vertex[0] = inner_shell_node;
335  nodes(itri).vertex[1] = inode-1;
336  nodes(itri).vertex[2] = inode;
337  itri++;
338  if (inodesh%ish>0){ // corners only add one triangle
339  nodes(itri).vertex[0] = inner_shell_node;
340  nodes(itri).vertex[1] = (inodesh == nodes_in_shell - 1) ? first_inner_shell_node : (inner_shell_node+1);
341  nodes(itri).vertex[2] = inode;
342  inner_shell_node++;
343  itri++;
344  }
345  }
346 
347  // Determine position of next node
348  int hex_region = inodesh/ish;
349  double r_move, z_move;
350  if(hex_region == NE)
351  {r_move = -0.5; z_move = 1.0;}
352  else if (hex_region == N)
353  {r_move = -1.0; z_move = 0.0;}
354  else if (hex_region == NW)
355  {r_move = -0.5; z_move = -1.0;}
356  else if (hex_region == SW)
357  {r_move = 0.5; z_move = -1.0;}
358  else if (hex_region == S)
359  {r_move = 1.0; z_move = 0.0;}
360  else if (hex_region == SE)
361  {r_move = 0.5; z_move = 1.0;}
362  r_track += rscale*r_move;
363  z_track += zscale*z_move;
364  inode++;
365  }
366 
367  // Add last element to close ring
368  nodes(itri).vertex[0] = first_inner_shell_node;
369  nodes(itri).vertex[1] = inode-1;
370  nodes(itri).vertex[2] = first_node;
371  itri++;
372 
373  nnodes_on_surface(ish) = nodes_in_shell;
374  }
375 
376  // Shift node list to 1-indexing (input files are 1-indexed)
377  for(int i = 0; i<ntriangles; i++){
378  nodes(i).vertex[0]++;
379  nodes(i).vertex[1]++;
380  nodes(i).vertex[2]++;
381  }
382 
383  // Field geometry
384  int n_xpts = 0; // Assume no X-points for now
385  int n_separatrices = 0; // Assume no separatrices for now
386  i_x = View<int*, HostType>("i_x", n_xpts);
387  i_surf_separatrices = View<int*, HostType>(NoInit("i_surf_separatrices"), n_separatrices);
388 
389  // Assume no non-aligned surfaces for now
390  num_non_aligned = 0;
391  non_aligned_nsurf = View<int*, HostType>(NoInit("non_aligned_nsurf"), num_non_aligned);
392  non_aligned_surf_idx = View<int**, CLayout, HostType>(NoInit("non_aligned_surf_idx"), num_non_aligned, 4);
393  }
394 
401  return sqrt((a.r-b.r)*(a.r-b.r) + (a.z-b.z)*(a.z-b.z));
402  }
403 
413  void construct_circular_grid(int nshells, double raxis, double zaxis, double rscale, double zscale) {
414  // First, build construction of what .node and .ele files would look like
415 
416  // Build circular sample grid of concentric shells
417  // Minimum 1 shell, which produces 7 nodes and 6 triangles
418 
419  // Work from center outward
420  x(0).r = raxis;
421  x(0).z = zaxis;
422  rgn(0) = 0; // not wall
423 
424  // Node indices of each surface
425  int max_nodes_on_surface = 6*nshells;
426  surf_idx = View<int**, CLayout, HostType>(NoInit("surf_idx"),nsurfaces, max_nodes_on_surface);
427 
428  nnodes_on_surface(0) = 1;
429  surf_idx(0, 0) = 1; // Since surf_idx is 1-indexed
430  int inode = 1;
431  int itri = 0;
432  for(int ish = 1; ish<=nshells; ish++){
433  int nodes_in_shell = 6*ish;
434  int nodes_in_inner_shell = ish==1 ? 1 : 6*(ish-1);
435  int first_inner_shell_node = ish==1 ? 0 : inode - 6*(ish-1);
436  int first_node = inode;
437  int inner_shell_node = first_inner_shell_node;
438  // Starting position:
439  for (int inodesh = 0; inodesh <= nodes_in_shell; inodesh++){
440  // Create node
441  if(inodesh<nodes_in_shell){
442  double angle = TWOPI*double(inodesh)/double(nodes_in_shell);
443  x(inode).r = raxis + ish*rscale*cos(angle);
444  x(inode).z = zaxis + ish*zscale*sin(angle);
445  rgn(inode) = (ish==nshells? OnWall : Inside); // last shell is wall
446 
447  surf_idx(ish, inodesh) = inode+1; // Since surf_idx is 1-indexed
448  }
449 
450  // Create element(s)
451  if (inodesh>0){
452  int outer_node = (inodesh < nodes_in_shell) ? inode : first_node;
453  int prev_outer_node = inode-1;
454 
455  // First, check cross-distances
456  int next_inner_shell_node = inner_shell_node+1;
457  if(next_inner_shell_node-first_inner_shell_node==nodes_in_inner_shell) next_inner_shell_node = first_inner_shell_node;
458  double distance1 = node_to_node_dist(x(outer_node), x(inner_shell_node));
459  double distance2 = node_to_node_dist(x(prev_outer_node), x(next_inner_shell_node));
460 
461  // The next node is far away enough that we need to add two elements
462  if(distance1>distance2 && next_inner_shell_node!=inner_shell_node){
463  if(itri>=nodes.size()) {printf("\nMore elements than expected in circular grid generation. Dynamically resize?"); exit(1); }
464  nodes(itri).vertex[0] = inner_shell_node;
465  nodes(itri).vertex[1] = next_inner_shell_node;
466  nodes(itri).vertex[2] = prev_outer_node;
467  inner_shell_node = next_inner_shell_node;
468  itri++;
469  }
470 
471  if(itri>=nodes.size()) {printf("\nMore elements than expected in circular grid generation. Dynamically resize?"); exit(1); }
472  nodes(itri).vertex[0] = inner_shell_node;
473  nodes(itri).vertex[1] = prev_outer_node;
474  nodes(itri).vertex[2] = outer_node;
475  itri++;
476  }
477 
478  if(inodesh<nodes_in_shell){
479  inode++;
480  }
481  }
482 
483  nnodes_on_surface(ish) = nodes_in_shell;
484  }
485 
486  if(itri!=nodes.size()) {printf("\nFewer elements than expected in circular grid generation. Dynamically resize?"); exit(1); }
487 
488  // Shift node list to 1-indexing (input files are 1-indexed)
489  for(int i = 0; i<ntriangles; i++){
490  nodes(i).vertex[0]++;
491  nodes(i).vertex[1]++;
492  nodes(i).vertex[2]++;
493  }
494 
495  // Field geometry
496  int n_xpts = 0; // Assume no X-points for now
497  int n_separatrices = 0; // Assume no separatrices for now
498  i_x = View<int*, HostType>("i_x", n_xpts);
499  i_surf_separatrices = View<int*, HostType>(NoInit("i_surf_separatrices"), n_separatrices);
500 
502  nsurfaces2 = 0;
503  nsurfaces3a = 0;
504  nsurfaces3b = 0;
505 
506  // Assume no non-aligned surfaces for now
507  num_non_aligned = 0;
508  non_aligned_nsurf = View<int*, HostType>(NoInit("non_aligned_nsurf"), num_non_aligned);
509  non_aligned_surf_idx = View<int**, CLayout, HostType>(NoInit("non_aligned_surf_idx"), num_non_aligned, 4);
510  non_aligned_vert = View<int*, HostType>(NoInit("non_aligned_vert"), num_non_aligned);
511  }
512 
513  // Returns the number of shells needed to create a grid with at least nnodes nodes.
514  static int nshells_from_nnodes(int nnodes){
515  // nnodes = 1 + 3*(nshells+1)*nshells
516  // 0 = 3s^2 + 3s + 1 - n
517  // s = (-3 +/- sqrt(9 + 4*3*(n-1)))/6
518  // s = sqrt(9 + 12*(n-1))/6 - 0.5
519  if(nnodes<1) return 0;
520  else return std::ceil(sqrt(12*nnodes - 3)/6 - 0.5);
521  }
522 
533  PlaneFiles(int nshells, double raxis, double zaxis, double rscale, double zscale, GridCreateOpt grid_opt=Hexagonal)
534  : nnodes(1 + 3*(nshells+1)*nshells),
535  x(NoInit("x"),nnodes),
536  rgn(NoInit("rgn"),nnodes),
537  ntriangles(6*nshells*nshells),
538  nodes(NoInit("nodes"),ntriangles),
539  nsurfaces(nshells+1),
540  nnodes_on_surface(NoInit("nnodes_on_surface"),nsurfaces),
542  {
543  if(grid_opt==Hexagonal){
544  construct_hex_grid(nshells, raxis, zaxis, rscale, zscale);
545  }else{
546  construct_circular_grid(nshells, raxis, zaxis, rscale, zscale);
547  }
548  }
549 
558  : nnodes(1 + 3*(nshells+1)*nshells),
559  x(NoInit("x"),nnodes),
560  rgn(NoInit("rgn"),nnodes),
561  ntriangles(6*nshells*nshells),
562  nodes(NoInit("nodes"),ntriangles),
563  nsurfaces(nshells+1),
564  nnodes_on_surface(NoInit("nnodes_on_surface"),nsurfaces),
566  {
567  // Magnetic axis will be the center of the grid
568  double raxis = magnetic_field.equil.axis_r;
569  double zaxis = magnetic_field.equil.axis_z;
570 
571  // Since the axis is not necessarily centered in the field map, take the shorter
572  // distance to the edge to use as the dimensions of the grid
573  double width = 2*std::min(magnetic_field.bounds.max_r - magnetic_field.equil.axis_r,
574  magnetic_field.equil.axis_r - magnetic_field.bounds.min_r);
575 
576  double height = 2*std::min(magnetic_field.bounds.max_z - magnetic_field.equil.axis_z,
577  magnetic_field.equil.axis_z - magnetic_field.bounds.min_z);
578 
579  // Make sure the grid is within the psi bounds
580  width *= 0.99;
581  height *= 0.99;
582 
583  if(grid_opt==Hexagonal){
584  double triangle_width = width/(2*nshells);
585  double triangle_height = height/(2*nshells);
586  construct_hex_grid(nshells, raxis, zaxis, triangle_width, triangle_height);
587  }else{
588  // Circular grid; take the smaller of the width or height to make sure it fits
589  double scale = std::min(width,height)/(2*nshells);
590  construct_circular_grid(nshells, raxis, zaxis, scale, scale);
591  }
592  }
593 
594  int n_specified_xpts() const{
595  return i_x.size();
596  }
597 
598  View<Equil::XPoint*, HostType> get_xpts() const{
599  View<Equil::XPoint*, HostType> xpts(NoInit("xpts"), i_x.size());
600  for(int i = 0; i<i_x.size(); i++){
601  xpts(i).r = x(i_x(i)).r;
602  xpts(i).z = x(i_x(i)).z;
603  }
604  return xpts;
605  }
606 };
607 
608 struct GridFiles{
612 
613  std::vector<PlaneFiles> plane_files_vec;
614 
615  // Allows array-like access pattern
616  const PlaneFiles& operator [](int i) const {return plane_files_vec[i];}
617 
618  // For test grids
619  GridFiles(int nshells, double raxis, double zaxis, double rscale, double zscale, GridCreateOpt grid_opt=Hexagonal){
620  plane_files_vec.push_back(PlaneFiles(nshells, raxis, zaxis, rscale, zscale, grid_opt));
621  }
622 
624  plane_files_vec.push_back(PlaneFiles(magnetic_field, nshells, grid_opt));
625  }
626 
627 #ifdef USE_MPI
628  /* Returns PlaneFiles object using the file-reader constructor. This may be updated so the input file
629  * can specify an analytic grid.
630  * */
631  GridFiles(NLReader::NamelistReader& nlr, bool is_stellarator, const MyMPI& my_mpi){
632  // Read file names and path from input file
633  nlr.use_namelist("sml_param");
634 
635  if(nlr.get<bool>("sml_create_analytic_grid", false)==true){
636  // Keep analytic grid options simple: number of surfaces, and axis_r
637  int nshells = 30;
638  double axis_r = 1.5;
639  if(nlr.namelist_present("analytic_grid_param")){
640  nlr.use_namelist("analytic_grid_param");
641 
642  nshells = nlr.get<int>("nsurfaces", 30);
643  axis_r = nlr.get<double>("axis_r", 1.5);
644  }
645 
646  if(is_rank_zero()) printf("\nCreating a circular analytic grid with %d surfaces and %d nodes.\n", nshells, 1 + 3*(nshells+1)*nshells);
647 
648  double axis_z = 0.0;
649 
650  // Leave a margin so that all vertices are inside equilibrium
651  double scale = 0.999/nshells; // Height of each triangle
652  double rscale = scale; // rscale = zscale for circular
653  double zscale = scale;
654  plane_files_vec.push_back(PlaneFiles(nshells, axis_r, axis_z, rscale, zscale, GridCreateOpt::Circular));
655 
656  if(is_stellarator){
657  // Push back a copy of the grid since stellarator needs two
658  plane_files_vec.push_back(PlaneFiles(nshells, axis_r, axis_z, rscale, zscale, GridCreateOpt::Circular));
659  }
660  }else{
661  std::string node_file = nlr.get<std::string>("sml_node_file","example_file.node");
662  std::string ele_file = nlr.get<std::string>("sml_ele_file","example_file.ele");
663  std::string surf_file = nlr.get<std::string>("sml_surf_file","example_file.flx.aif");
664  std::string input_file_dir = nlr.get<std::string>("sml_input_file_dir","./");
665 
666  // Add path to file names
667  node_file = input_file_dir + "/" + node_file;
668  ele_file = input_file_dir + "/" + ele_file;
669  surf_file = input_file_dir + "/" + surf_file;
670 
671  if(is_stellarator){
672  int iplane = my_mpi.my_intpl_rank;
673  int ifile = iplane*2; // Midplane is .0, Edge is .1, next midplane is .2, etc.
674  // Make some assumption about what these files are temporarily called
675  std::string midplane_node_file = node_file + "." + std::to_string(ifile);
676  std::string midplane_ele_file = ele_file + "." + std::to_string(ifile);
677  std::string midplane_surf_file = surf_file + "." + std::to_string(ifile);
678 
679  // Read in midplane grid files (one rank per plane reads)
680  plane_files_vec.push_back(PlaneFiles(midplane_node_file, midplane_ele_file, midplane_surf_file, my_mpi.plane_comm));
681 
682  ifile += 1; // shift to edge plane
683  node_file = node_file + "." + std::to_string(ifile);
684  ele_file = ele_file + "." + std::to_string(ifile);
685  surf_file = surf_file + "." + std::to_string(ifile);
686 
687  // Read in edgeplane grid files (one rank per plane reads)
688  plane_files_vec.push_back(PlaneFiles(node_file, ele_file, surf_file, my_mpi.plane_comm));
689  }else{
690  // Read in grid files of single plane
691  plane_files_vec.push_back(PlaneFiles(node_file, ele_file, surf_file, SML_COMM_WORLD));
692  }
693  }
694  }
695 #endif
696 
697  void clear(){
698  plane_files_vec.clear();
699  }
700 
701  // Returns the number of shells needed to create a grid with at least nnodes nodes.
702  static int nshells_from_nnodes(int nnodes){
703  return PlaneFiles::nshells_from_nnodes(nnodes);
704  }
705 
706 };
707 
708 #endif
Definition: grid_files.hpp:285
PlaneFiles(int nshells, double raxis, double zaxis, double rscale, double zscale, GridCreateOpt grid_opt=Hexagonal)
Definition: grid_files.hpp:533
static constexpr GridCreateOpt Hexagonal
Definition: grid_files.hpp:610
bool is_rank_zero()
Definition: globals.hpp:27
MPI_Comm SML_COMM_WORLD
Definition: my_mpi.cpp:4
View< int **, CLayout, HostType > surf_idx
Definition: grid_files.hpp:39
MPI_Comm plane_comm
Definition: my_mpi.hpp:38
T get(const string &param, const T default_val, int val_ind=0)
Definition: NamelistReader.hpp:373
static int nshells_from_nnodes(int nnodes)
Definition: grid_files.hpp:514
static constexpr GridCreateOpt Circular
Definition: grid_files.hpp:611
int my_intpl_rank
Definition: my_mpi.hpp:49
void set_fortran_flx_aif_format(int is_flx_aif_format_int)
Definition: my_mpi.F90:1
PlaneFiles()
Definition: grid_files.hpp:60
View< int *, HostType > non_aligned_vert
Definition: grid_files.hpp:36
Definition: NamelistReader.hpp:193
Definition: magnetic_field.hpp:12
View< int *, HostType > non_aligned_nsurf
Definition: grid_files.hpp:37
int ntriangles
Definition: grid_files.hpp:24
RZBounds bounds
Simulation boundary.
Definition: magnetic_field.hpp:23
Definition: grid_files.hpp:608
Equilibrium equil
The object containing information about the magnetic equilibrium.
Definition: magnetic_field.hpp:32
void read(T *var, int n)
Definition: file_reader.hpp:28
int num_non_aligned
Definition: grid_files.hpp:35
void construct_circular_grid(int nshells, double raxis, double zaxis, double rscale, double zscale)
Definition: grid_files.hpp:413
View< Equil::XPoint *, HostType > get_xpts() const
Definition: grid_files.hpp:598
int nsurfaces3a
Definition: grid_files.hpp:31
View< int *, HostType > nnodes_on_surface
Definition: grid_files.hpp:33
int n_specified_xpts() const
Definition: grid_files.hpp:594
Definition: grid_files.hpp:284
GridFiles(int nshells, double raxis, double zaxis, double rscale, double zscale, GridCreateOpt grid_opt=Hexagonal)
Definition: grid_files.hpp:619
double z
Definition: grid_structs.hpp:30
double axis_z
z coordinate of axis
Definition: equil.hpp:90
Definition: grid_structs.hpp:28
int nsurfaces2
Definition: grid_files.hpp:31
void convert_to_x_rgn_format(const View< double **, CLayout, HostType > &tmp_nodes, View< RZPair *, HostType > &x, View< int *, HostType > &rgn)
Definition: grid_files.hpp:41
std::vector< PlaneFiles > plane_files_vec
Definition: grid_files.hpp:613
Definition: my_mpi.hpp:19
GridFiles(MagneticField< HostType > &magnetic_field, int nshells, GridCreateOpt grid_opt=Hexagonal)
Definition: grid_files.hpp:623
void use_namelist(const string &namelist)
Definition: NamelistReader.hpp:355
View< int *, HostType > i_x
Definition: grid_files.hpp:28
double max_z
Definition: rz_bounds.hpp:8
double axis_r
r coordinate of axis
Definition: equil.hpp:89
double min_r
Definition: rz_bounds.hpp:5
Definition: col_grid.hpp:34
int nsurfaces1
Definition: grid_files.hpp:31
Definition: file_reader.hpp:6
void convert_to_node_format(const View< int **, CLayout, HostType > &tmp_ele, View< Vertex *, HostType > &nodes)
Definition: grid_files.hpp:50
const PlaneFiles & operator[](int i) const
Definition: grid_files.hpp:616
GridCreateOpt
Definition: grid_files.hpp:283
View< int **, CLayout, HostType > non_aligned_surf_idx
Definition: grid_files.hpp:38
bool namelist_present(const string &namelist)
Definition: NamelistReader.hpp:351
int nnodes
Definition: grid_files.hpp:19
void construct_hex_grid(int nshells, double raxis, double zaxis, double rscale, double zscale)
Definition: grid_files.hpp:297
int nsurfaces3b
Definition: grid_files.hpp:31
double r
Definition: grid_structs.hpp:29
void exit_XGC(std::string msg)
Definition: globals.hpp:37
Definition: magnetic_field.F90:1
PlaneFiles(MagneticField< HostType > &magnetic_field, int nshells, GridCreateOpt grid_opt=Hexagonal)
Definition: grid_files.hpp:557
static int nshells_from_nnodes(int nnodes)
Definition: grid_files.hpp:702
Definition: grid_files.hpp:17
double node_to_node_dist(RZPair a, RZPair b)
Definition: grid_files.hpp:400
void clear()
Definition: grid_files.hpp:697
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
double min_z
Definition: rz_bounds.hpp:7
View< int *, HostType > rgn
Definition: grid_files.hpp:21
int nsurfaces
Definition: grid_files.hpp:30
void skip(int n)
Definition: file_reader.hpp:45
View< RZPair *, HostType > x
Definition: grid_files.hpp:20
Kokkos::ViewAllocateWithoutInitializing NoInit
Definition: space_settings.hpp:68
constexpr double TWOPI
Definition: constants.hpp:9
bool is_old_flx_format(const std::string &filename)
Definition: check_aif_flx_format.hpp:13
View< int *, HostType > i_surf_separatrices
Definition: grid_files.hpp:32
View< Vertex *, HostType > nodes
Definition: grid_files.hpp:25
int sum_nnodes_on_surfaces
Definition: grid_files.hpp:34
void open(const std::string &input)
Definition: file_reader.hpp:13
double max_r
Definition: rz_bounds.hpp:6