XGCa
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator 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 
14 
15 struct GridFiles{
16  // Contents of .node
17  int nnodes;
18  View<RZPair*, HostType> x;
19  View<int*, HostType> rgn;
20 
21  // Contents of .ele
23  View<Vertex*, HostType> nodes;
24 
25  // Contents of .flx.aif
26  int nsurfaces;
27  View<int*, HostType> nnodes_on_surface;
29  View<int*, HostType> nodes_of_surfaces;
30 
31  void convert_to_x_rgn_format(const View<double**, CLayout, HostType>& tmp_nodes, View<RZPair*, HostType>& x, View<int*, HostType>& rgn){
32  Kokkos::parallel_for("transpose_x_and_rgn", Kokkos::RangePolicy<HostExSpace>(0,tmp_nodes.extent(0)), KOKKOS_LAMBDA( const int i ){
33  // tmp_nodes[i,0]; <-- first column is ignored
34  x(i).r = tmp_nodes(i,1);
35  x(i).z = tmp_nodes(i,2);
36  rgn(i) = tmp_nodes(i,3);
37  });
38  }
39 
40  void convert_to_node_format(const View<int**, CLayout, HostType>& tmp_ele, View<Vertex*, HostType>& nodes){
41  Kokkos::parallel_for("transpose_x_and_rgn", Kokkos::RangePolicy<HostExSpace>(0,tmp_ele.extent(0)), KOKKOS_LAMBDA( const int i ){
42  // tmp_ele[i,0]; <-- first column is ignored
43  nodes(i).vertex[0] = tmp_ele(i,1);
44  nodes(i).vertex[1] = tmp_ele(i,2);
45  nodes(i).vertex[2] = tmp_ele(i,3);
46  });
47  }
48 
49  // Read the files
50  GridFiles(const std::string& node_file, const std::string& element_file, const std::string& flx_aif_file){
51  // Only one MPI rank reads the files, then broadcasts to the others
52  if(is_rank_zero()){
53 
54  FileReader file_reader;
55 
57  file_reader.open(node_file);
58  file_reader.read(&nnodes, 1); // Get nnodes from file
59  file_reader.skip(3); // Skip 3 unused values in file
60 
61  // Use dummy array since reformatting is needed
62  View<double**, CLayout, HostType> tmp_nodes(NoInit("tmp_nodes"),nnodes,4); // 4 entries per node: index, r, z, region
63  file_reader.read(tmp_nodes.data(), tmp_nodes.size());
64 
65  // Reformat to x and rgn
66  x = View<RZPair*, HostType>(NoInit("x"),nnodes);
67  rgn = View<int*, HostType>(NoInit("rgn"),nnodes);
68  convert_to_x_rgn_format(tmp_nodes, x, rgn);
69 
71  file_reader.open(element_file);
72  file_reader.read(&ntriangles, 1); // Get nnodes from file
73  file_reader.skip(2); // Skip 2 unused values in file
74 
75  // Use dummy array since reformatting is needed
76  View<int**, CLayout, HostType> tmp_ele(NoInit("tmp_ele"),ntriangles,4); // 4 entries per node: index, r, z, region
77  file_reader.read(tmp_ele.data(), tmp_ele.size());
78 
79  // Reformat since first column is ignored
80  nodes = View<Vertex*, HostType>("nodes",ntriangles);
81  convert_to_node_format(tmp_ele, nodes);
82 
84 /*
85 #ifndef NEW_FLX_AIF
86  exit_XGC("\nError: C++ grid file reader can't read old flux files (i.e. NEW_FLX_AIF must be On)");
87 #endif
88  file_reader.open(flx_aif_file);
89  file_reader.skip(3); // Skip 3 unused values in file
90  int nsurfaces1, nsurfaces2, nsurfaces3a, nsurfaces3b;
91  file_reader.read(&nsurfaces1, 1); // Get number of surfaces
92  file_reader.read(&nsurfaces2, 1); // Get number of surfaces
93  file_reader.read(&nsurfaces3a, 1); // Get number of surfaces
94  file_reader.read(&nsurfaces3b, 1); // Get number of surfaces
95  // nsurfaces is actually the sum of these 4 inputs
96  nsurfaces = nsurfaces1 + nsurfaces2 + nsurfaces3a + nsurfaces3b;
97 
98  file_reader.skip(2); // Skip 2 unused values in file
99 
100  // Read in nnodes on each surface
101  nnodes_on_surface = std::vector<int>(nsurfaces);
102  file_reader.read(nnodes_on_surface.data(), nnodes_on_surface.size());
103 
104  // Loop through surfaces to calculate sum total of nodes on surfaces
105  sum_nnodes_on_surfaces = 0;
106  for(int i=0; i<nsurfaces; i++) sum_nnodes_on_surfaces += nnodes_on_surface[i];
107 
108  // Allocate nodes_of_surfaces
109  nodes_of_surfaces = std::vector<int>(sum_nnodes_on_surfaces);
110 
111  // Read in values
112  file_reader.read(nodes_of_surfaces.data(), nodes_of_surfaces.size());
113 */
114  }
115 
116  /*** Broadcast grid contents to other MPI ranks ***/
117 
118  // First broadcast scalars so we know the sizes of the arrays
119  View<int*, HostType> scalars("scalars", 4);
120  scalars(0) = nnodes;
121  scalars(1) = ntriangles;
122 // scalars(2) = nsurfaces;
123 // scalars(3) = sum_nnodes_on_surfaces;
124 
125  // Broadcast file size to other processes
126 #ifdef USE_MPI
127  int root_rank = 0;
128  MPI_Bcast(scalars.data(), scalars.size(), MPI_INT, root_rank, SML_COMM_WORLD);
129 #endif
130 
131  if(!is_rank_zero()){
132  // Set scalars on non-root ranks
133  nnodes = scalars(0);
134  ntriangles = scalars(1);
135 // nsurfaces = scalars(2);
136 // sum_nnodes_on_surfaces = scalars(3);
137 
138  // Allocate vectors on non-root ranks
139  Kokkos::resize(x, nnodes);
140  Kokkos::resize(rgn, nnodes);
141  Kokkos::resize(nodes, ntriangles);
142 // nnodes_on_surface.resize(nsurfaces);
143 // nodes_of_surfaces.resize(sum_nnodes_on_surfaces);
144  }
145 
146  // Broadcast grid arrays
147 #ifdef USE_MPI
148  MPI_Bcast( x.data(), x.size()*2, MPI_DOUBLE, root_rank, SML_COMM_WORLD);
149  MPI_Bcast( rgn.data(), rgn.size(), MPI_INT, root_rank, SML_COMM_WORLD);
150  MPI_Bcast( nodes.data(), nodes.size()*3, MPI_INT, root_rank, SML_COMM_WORLD);
151 // MPI_Bcast( nnodes_on_surface.data(), nnodes_on_surface.size(), MPI_INT, root_rank, SML_COMM_WORLD);
152 // MPI_Bcast( nodes_of_surfaces.data(), nodes_of_surfaces.size(), MPI_INT, root_rank, SML_COMM_WORLD);
153 #endif
154  }
155 
156  // Options: hex vs circular
160  };
161 
171  void construct_hex_grid(int nshells, double raxis, double zaxis, double rscale, double zscale) {
172  // First, build construction of what .node and .ele files would look like
173 
174  // Build hexagonal sample grid of concentric shells
175  // Minimum 1 shell, which produces 7 nodes and 6 triangles
176 
177  // Grid construction uses this enum for directions
178  enum { NE = 0, N, NW, SW, S, SE };
179 
180  // Work from center outward
181  x(0).r = raxis;
182  x(0).z = zaxis;
183  rgn(0) = 0; // not wall
184  int inode = 1;
185  int itri = 0;
186  for(int ish = 1; ish<=nshells; ish++){
187  int nodes_in_shell = 6*ish;
188  int first_inner_shell_node = ish==1 ? 0 : inode - 6*(ish-1);
189  int first_node = inode;
190  int inner_shell_node = first_inner_shell_node;
191  // Starting position:
192  double r_track = raxis + rscale*ish;
193  double z_track = zaxis;
194  for (int inodesh = 0; inodesh < nodes_in_shell; inodesh++){
195  x(inode).r = r_track;
196  x(inode).z = z_track;
197  rgn(inode) = (ish==nshells? OnWall : Inside); // last shell is wall
198 
199  if (inodesh>0){
200  nodes(itri).vertex[0] = inner_shell_node;
201  nodes(itri).vertex[1] = inode-1;
202  nodes(itri).vertex[2] = inode;
203  itri++;
204  if (inodesh%ish>0){ // corners only add one triangle
205  nodes(itri).vertex[0] = inner_shell_node;
206  nodes(itri).vertex[1] = (inodesh == nodes_in_shell - 1) ? first_inner_shell_node : (inner_shell_node+1);
207  nodes(itri).vertex[2] = inode;
208  inner_shell_node++;
209  itri++;
210  }
211  }
212 
213  // Determine position of next node
214  int hex_region = inodesh/ish;
215  double r_move, z_move;
216  if(hex_region == NE)
217  {r_move = -0.5; z_move = 1.0;}
218  else if (hex_region == N)
219  {r_move = -1.0; z_move = 0.0;}
220  else if (hex_region == NW)
221  {r_move = -0.5; z_move = -1.0;}
222  else if (hex_region == SW)
223  {r_move = 0.5; z_move = -1.0;}
224  else if (hex_region == S)
225  {r_move = 1.0; z_move = 0.0;}
226  else if (hex_region == SE)
227  {r_move = 0.5; z_move = 1.0;}
228  r_track += rscale*r_move;
229  z_track += zscale*z_move;
230  inode++;
231  }
232 
233  // Add last element to close ring
234  nodes(itri).vertex[0] = first_inner_shell_node;
235  nodes(itri).vertex[1] = inode-1;
236  nodes(itri).vertex[2] = first_node;
237  itri++;
238  }
239 
240  // Shift node list to 1-indexing (input files are 1-indexed)
241  for(int i = 0; i<ntriangles; i++){
242  nodes(i).vertex[0]++;
243  nodes(i).vertex[1]++;
244  nodes(i).vertex[2]++;
245  }
246 
247  // nsurfaces is 0 because the hex grid is not surface-aligned
248  nsurfaces=0;
249  }
250 
257  return sqrt((a.r-b.r)*(a.r-b.r) + (a.z-b.z)*(a.z-b.z));
258  }
259 
269  void construct_circular_grid(int nshells, double raxis, double zaxis, double rscale, double zscale) {
270  // First, build construction of what .node and .ele files would look like
271 
272  // Build circular sample grid of concentric shells
273  // Minimum 1 shell, which produces 7 nodes and 6 triangles
274 
275  // Work from center outward
276  x(0).r = raxis;
277  x(0).z = zaxis;
278  rgn(0) = 0; // not wall
279 
280  nnodes_on_surface(0) = 1;
281  nodes_of_surfaces(0) = 0;
282  int inode = 1;
283  int itri = 0;
284  for(int ish = 1; ish<=nshells; ish++){
285  int nodes_in_shell = 6*ish;
286  int nodes_in_inner_shell = ish==1 ? 1 : 6*(ish-1);
287  int first_inner_shell_node = ish==1 ? 0 : inode - 6*(ish-1);
288  int first_node = inode;
289  int inner_shell_node = first_inner_shell_node;
290  // Starting position:
291  for (int inodesh = 0; inodesh <= nodes_in_shell; inodesh++){
292  // Create node
293  if(inodesh<nodes_in_shell){
294  double angle = TWOPI*double(inodesh)/double(nodes_in_shell);
295  x(inode).r = raxis + ish*rscale*cos(angle);
296  x(inode).z = zaxis + ish*zscale*sin(angle);
297  rgn(inode) = (ish==nshells? OnWall : Inside); // last shell is wall
298  }
299 
300  // Create element(s)
301  if (inodesh>0){
302  int outer_node = (inodesh < nodes_in_shell) ? inode : first_node;
303  int prev_outer_node = inode-1;
304 
305  // First, check cross-distances
306  int next_inner_shell_node = inner_shell_node+1;
307  if(next_inner_shell_node-first_inner_shell_node==nodes_in_inner_shell) next_inner_shell_node = first_inner_shell_node;
308  double distance1 = node_to_node_dist(x(outer_node), x(inner_shell_node));
309  double distance2 = node_to_node_dist(x(prev_outer_node), x(next_inner_shell_node));
310 
311  // The next node is far away enough that we need to add two elements
312  if(distance1>distance2 && next_inner_shell_node!=inner_shell_node){
313  if(itri>=nodes.size()) {printf("\nMore elements than expected in circular grid generation. Dynamically resize?"); exit(1); }
314  nodes(itri).vertex[0] = inner_shell_node;
315  nodes(itri).vertex[1] = next_inner_shell_node;
316  nodes(itri).vertex[2] = prev_outer_node;
317  inner_shell_node = next_inner_shell_node;
318  itri++;
319  }
320 
321  if(itri>=nodes.size()) {printf("\nMore elements than expected in circular grid generation. Dynamically resize?"); exit(1); }
322  nodes(itri).vertex[0] = inner_shell_node;
323  nodes(itri).vertex[1] = prev_outer_node;
324  nodes(itri).vertex[2] = outer_node;
325  itri++;
326  }
327 
328  if(inodesh<nodes_in_shell){
329  nodes_of_surfaces(inode) = inode;
330  inode++;
331  }
332  }
333 
334  nnodes_on_surface(ish) = nodes_in_shell;
335  }
336 
337  if(itri!=nodes.size()) {printf("\nFewer elements than expected in circular grid generation. Dynamically resize?"); exit(1); }
338 
339  // Shift node list to 1-indexing (input files are 1-indexed)
340  for(int i = 0; i<ntriangles; i++){
341  nodes(i).vertex[0]++;
342  nodes(i).vertex[1]++;
343  nodes(i).vertex[2]++;
344  }
345  }
346 
347 
358  GridFiles(int nshells, double raxis, double zaxis, double rscale, double zscale, GridCreateOpt grid_opt=Hexagonal)
359  : nnodes(1 + 3*(nshells+1)*nshells),
360  x(NoInit("x"),nnodes),
361  rgn(NoInit("rgn"),nnodes),
362  ntriangles(6*nshells*nshells),
363  nodes(NoInit("nodes"),ntriangles),
364  nsurfaces(nshells+1),
365  nnodes_on_surface(NoInit("nnodes_on_surface"),nsurfaces),
367  nodes_of_surfaces(NoInit("nodes_of_surfaces"),nnodes)
368  {
369  if(grid_opt==Hexagonal){
370  construct_hex_grid(nshells, raxis, zaxis, rscale, zscale);
371  }else{
372  construct_circular_grid(nshells, raxis, zaxis, rscale, zscale);
373  }
374  }
375 
384  : nnodes(1 + 3*(nshells+1)*nshells),
385  x(NoInit("x"),nnodes),
386  rgn(NoInit("rgn"),nnodes),
387  ntriangles(6*nshells*nshells),
388  nodes(NoInit("nodes"),ntriangles),
389  nsurfaces(nshells+1),
390  nnodes_on_surface(NoInit("nnodes_on_surface"),nsurfaces),
392  nodes_of_surfaces(NoInit("nodes_of_surfaces"),nnodes)
393  {
394  // Magnetic axis will be the center of the grid
395  double raxis = magnetic_field.equil.axis_r;
396  double zaxis = magnetic_field.equil.axis_z;
397 
398  // Since the axis is not necessarily centered in the field map, take the shorter
399  // distance to the edge to use as the dimensions of the grid
400  double width = 2*std::min(magnetic_field.bd_max_r - magnetic_field.equil.axis_r,
401  magnetic_field.equil.axis_r - magnetic_field.bd_min_r);
402 
403  double height = 2*std::min(magnetic_field.bd_max_z - magnetic_field.equil.axis_z,
404  magnetic_field.equil.axis_z - magnetic_field.bd_min_z);
405 
406  // Make sure the grid is within the psi bounds
407  width *= 0.99;
408  height *= 0.99;
409 
410  if(grid_opt==Hexagonal){
411  double triangle_width = width/(2*nshells);
412  double triangle_height = height/(2*nshells);
413  construct_hex_grid(nshells, raxis, zaxis, triangle_width, triangle_height);
414  }else{
415  // Circular grid; take the smaller of the width or height to make sure it fits
416  double scale = std::min(width,height)/(2*nshells);
417  construct_circular_grid(nshells, raxis, zaxis, scale, scale);
418  }
419  }
420 
421  /* Returns GridFiles object using the file-reader constructor. This may be updated so the input file
422  * can specify an analytic grid.
423  * */
425  // Read file names and path from input file
426  nlr.use_namelist("sml_param");
427  std::string node_file = nlr.get<std::string>("sml_node_file","example_file.node");
428  std::string ele_file = nlr.get<std::string>("sml_ele_file","example_file.ele");
429  std::string surf_file = nlr.get<std::string>("sml_surf_file","example_file.flx.aif");
430  std::string input_file_dir = nlr.get<std::string>("sml_input_file_dir","./");
431 
432  // Add path to file names
433  node_file = input_file_dir + "/" + node_file;
434  ele_file = input_file_dir + "/" + ele_file;
435  surf_file = input_file_dir + "/" + surf_file;
436 
437  // Read in grid files and return result
438  return GridFiles(node_file, ele_file, surf_file);
439  }
440 
441 };
442 
443 #endif
double bd_max_z
Maximum z (of what exactly?)
Definition: magnetic_field.hpp:23
View< int *, HostType > nodes_of_surfaces
Definition: grid_files.hpp:29
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:31
bool is_rank_zero()
Definition: globals.hpp:26
MPI_Comm SML_COMM_WORLD
Definition: my_mpi.cpp:4
T get(const string &param, const T default_val, int val_ind=0)
Definition: NamelistReader.hpp:353
View< Vertex *, HostType > nodes
Definition: grid_files.hpp:23
int nsurfaces
Definition: grid_files.hpp:26
Definition: NamelistReader.hpp:163
Definition: magnetic_field.hpp:9
Definition: grid_files.hpp:15
static GridFiles grid_files_from_namelist(NLReader::NamelistReader &nlr)
Definition: grid_files.hpp:424
void read(T *var, int n)
Definition: file_reader.hpp:28
View< int *, HostType > nnodes_on_surface
Definition: grid_files.hpp:27
void construct_circular_grid(int nshells, double raxis, double zaxis, double rscale, double zscale)
Definition: grid_files.hpp:269
View< int *, HostType > rgn
Definition: grid_files.hpp:19
GridFiles(int nshells, double raxis, double zaxis, double rscale, double zscale, GridCreateOpt grid_opt=Hexagonal)
Definition: grid_files.hpp:358
double z
Definition: grid_structs.hpp:30
Definition: grid_structs.hpp:28
double bd_min_r
Minimum r (of what exactly?)
Definition: magnetic_field.hpp:20
View< RZPair *, HostType > x
Definition: grid_files.hpp:18
Definition: grid_files.hpp:159
GridFiles(MagneticField< HostType > &magnetic_field, int nshells, GridCreateOpt grid_opt=Hexagonal)
Definition: grid_files.hpp:383
int ntriangles
Definition: grid_files.hpp:22
void use_namelist(const string &namelist)
Definition: NamelistReader.hpp:322
Definition: grid_files.hpp:158
double bd_min_z
Minimum z (of what exactly?)
Definition: magnetic_field.hpp:22
void convert_to_node_format(const View< int **, CLayout, HostType > &tmp_ele, View< Vertex *, HostType > &nodes)
Definition: grid_files.hpp:40
Definition: col_grid.hpp:35
Definition: file_reader.hpp:6
constexpr double TWOPI
Definition: globals.hpp:123
void construct_hex_grid(int nshells, double raxis, double zaxis, double rscale, double zscale)
Definition: grid_files.hpp:171
double r
Definition: grid_structs.hpp:29
Equilibrium< Device > equil
The object containing information about the magnetic equilibrium.
Definition: magnetic_field.hpp:32
Definition: magnetic_field.F90:1
int sum_nnodes_on_surfaces
Definition: grid_files.hpp:28
GridCreateOpt
Definition: grid_files.hpp:157
GridFiles(const std::string &node_file, const std::string &element_file, const std::string &flx_aif_file)
Definition: grid_files.hpp:50
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
int nnodes
Definition: grid_files.hpp:17
void skip(int n)
Definition: file_reader.hpp:45
Kokkos::ViewAllocateWithoutInitializing NoInit
Definition: space_settings.hpp:61
double node_to_node_dist(RZPair a, RZPair b)
Definition: grid_files.hpp:256
void open(const std::string &input)
Definition: file_reader.hpp:13
double bd_max_r
Maximum r (of what exactly?)
Definition: magnetic_field.hpp:21