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