25 View<RZPair*, HostType>
x;
26 View<int*, HostType>
rgn;
33 View<int*, HostType>
i_x;
47 Kokkos::parallel_for(
"transpose_x_and_rgn", Kokkos::RangePolicy<HostExSpace>(0,tmp_nodes.extent(0)), KOKKOS_LAMBDA(
const int i ){
49 x(i).r = tmp_nodes(i,1);
50 x(i).z = tmp_nodes(i,2);
51 rgn(i) = tmp_nodes(i,3);
56 Kokkos::parallel_for(
"transpose_x_and_rgn", Kokkos::RangePolicy<HostExSpace>(0,tmp_ele.extent(0)), KOKKOS_LAMBDA(
const int i ){
58 nodes(i).vertex[0] = tmp_ele(i,1);
59 nodes(i).vertex[1] = tmp_ele(i,2);
60 nodes(i).vertex[2] = tmp_ele(i,3);
70 PlaneFiles(
const std::string& node_file,
const std::string& element_file,
const std::string& flx_aif_file,
const MPI_Comm& comm){
74 MPI_Comm_rank(comm, &my_rank);
80 int max_nodes_on_surface;
84 if(my_rank==ROOT_RANK){
89 file_reader.
open(node_file);
94 View<double**, CLayout, HostType> tmp_nodes(
NoInit(
"tmp_nodes"),
nnodes,4);
95 file_reader.
read(tmp_nodes.data(), tmp_nodes.size());
103 file_reader.
open(element_file);
109 file_reader.
read(tmp_ele.data(), tmp_ele.size());
120 file_reader.
open(flx_aif_file);
122 if(old_flx_aif_format){
126 i_x = View<int*, HostType>(
"i_x",0);
135 file_reader.
read(&n_xpts, 1);
136 i_x = View<int*, HostType>(
"i_x", n_xpts);
140 file_reader.
read(&i_xpt1, 1);
141 file_reader.
read(&i_xpt2, 1);
144 if(n_xpts==1 && i_xpt2 != -1)
exit_XGC(
"Error in flx.aif file: n_xpts=1, but i_xpt2 is specified.\n");
145 if(n_xpts==2 && i_xpt2 <= 0)
exit_XGC(
"Error in flx.aif file: n_xpts=2, but i_xpt2 is not specified.\n");
147 if(n_xpts>=1)
i_x(0) = i_xpt1 - 1;
148 if(n_xpts>=2)
i_x(1) = i_xpt2 - 1;
172 max_nodes_on_surface = 0;
176 surf_idx = View<int**, CLayout, HostType>(
"surf_idx",
nsurfaces, max_nodes_on_surface);
184 int should_be_negative_one;
185 file_reader.
read(&should_be_negative_one, 1);
186 if(should_be_negative_one != -1){
187 printf(
"\nExpected -1 at end of section in surface file, encountered %d instead\n", should_be_negative_one);
188 exit_XGC(
"Error reading surface file.\n");
201 if(old_flx_aif_format){
211 file_reader.
read(&should_be_negative_one, 1);
212 if(should_be_negative_one != -1){
213 printf(
"\nExpected -1 at end of the surface file, encountered %d instead\n", should_be_negative_one);
214 exit_XGC(
"Error reading surface file.\n");
222 View<int*, HostType> scalars(
"scalars", n_scalars);
226 scalars(3) = max_nodes_on_surface;
229 scalars(6) = n_separatrices;
237 MPI_Bcast(scalars.data(), scalars.size(), MPI_INT, ROOT_RANK, comm);
240 if(my_rank!=ROOT_RANK){
244 nsurfaces = scalars(2);
245 max_nodes_on_surface = scalars(3);
246 num_non_aligned = scalars(4);
248 n_separatrices = scalars(6);
252 nsurfaces3b = scalars(10);
261 i_x = View<int*, HostType>(
"i_x", n_xpts);
265 non_aligned_surf_idx = View<int**, CLayout, HostType>(
NoInit(
"non_aligned_surf_idx"),
num_non_aligned, 4);
272 MPI_Bcast(
x.data(),
x.size()*2, MPI_DOUBLE, ROOT_RANK, comm);
273 MPI_Bcast( rgn.data(), rgn.size(), MPI_INT, ROOT_RANK, comm);
274 MPI_Bcast(
nodes.data(),
nodes.size()*3, MPI_INT, ROOT_RANK, comm);
276 if(
i_x.size()>0) MPI_Bcast(
i_x.data(),
i_x.size(), MPI_INT, ROOT_RANK, comm);
279 MPI_Bcast( non_aligned_nsurf.data(), non_aligned_nsurf.size(), MPI_INT, ROOT_RANK, comm);
280 MPI_Bcast( non_aligned_surf_idx.data(), non_aligned_surf_idx.size(), MPI_INT, ROOT_RANK, comm);
309 enum { NE = 0,
N, NW, SW, S, SE };
312 int max_nodes_on_surface = 6*nshells;
323 for(
int ish = 1; ish<=nshells; ish++){
324 int nodes_in_shell = 6*ish;
325 int first_inner_shell_node = ish==1 ? 0 : inode - 6*(ish-1);
326 int first_node = inode;
327 int inner_shell_node = first_inner_shell_node;
329 double r_track = raxis + rscale*ish;
330 double z_track = zaxis;
331 for (
int inodesh = 0; inodesh < nodes_in_shell; inodesh++){
332 x(inode).r = r_track;
333 x(inode).z = z_track;
339 nodes(itri).vertex[0] = inner_shell_node;
340 nodes(itri).vertex[1] = inode-1;
341 nodes(itri).vertex[2] = inode;
344 nodes(itri).vertex[0] = inner_shell_node;
345 nodes(itri).vertex[1] = (inodesh == nodes_in_shell - 1) ? first_inner_shell_node : (inner_shell_node+1);
346 nodes(itri).vertex[2] = inode;
353 int hex_region = inodesh/ish;
354 double r_move, z_move;
356 {r_move = -0.5; z_move = 1.0;}
357 else if (hex_region ==
N)
358 {r_move = -1.0; z_move = 0.0;}
359 else if (hex_region == NW)
360 {r_move = -0.5; z_move = -1.0;}
361 else if (hex_region == SW)
362 {r_move = 0.5; z_move = -1.0;}
363 else if (hex_region == S)
364 {r_move = 1.0; z_move = 0.0;}
365 else if (hex_region == SE)
366 {r_move = 0.5; z_move = 1.0;}
367 r_track += rscale*r_move;
368 z_track += zscale*z_move;
373 nodes(itri).vertex[0] = first_inner_shell_node;
374 nodes(itri).vertex[1] = inode-1;
375 nodes(itri).vertex[2] = first_node;
383 nodes(i).vertex[0]++;
384 nodes(i).vertex[1]++;
385 nodes(i).vertex[2]++;
390 int n_separatrices = 0;
391 i_x = View<int*, HostType>(
"i_x", n_xpts);
397 non_aligned_surf_idx = View<int**, CLayout, HostType>(
NoInit(
"non_aligned_surf_idx"),
num_non_aligned, 4);
406 return sqrt((a.
r-b.
r)*(a.
r-b.
r) + (a.
z-b.
z)*(a.
z-b.
z));
430 int max_nodes_on_surface = 6*nshells;
437 for(
int ish = 1; ish<=nshells; ish++){
438 int nodes_in_shell = 6*ish;
439 int nodes_in_inner_shell = ish==1 ? 1 : 6*(ish-1);
440 int first_inner_shell_node = ish==1 ? 0 : inode - 6*(ish-1);
441 int first_node = inode;
442 int inner_shell_node = first_inner_shell_node;
444 for (
int inodesh = 0; inodesh <= nodes_in_shell; inodesh++){
446 if(inodesh<nodes_in_shell){
447 double angle =
TWOPI*double(inodesh)/double(nodes_in_shell);
448 x(inode).r = raxis + ish*rscale*cos(angle);
449 x(inode).z = zaxis + ish*zscale*sin(angle);
457 int outer_node = (inodesh < nodes_in_shell) ? inode : first_node;
458 int prev_outer_node = inode-1;
461 int next_inner_shell_node = inner_shell_node+1;
462 if(next_inner_shell_node-first_inner_shell_node==nodes_in_inner_shell) next_inner_shell_node = first_inner_shell_node;
467 if(distance1>distance2 && next_inner_shell_node!=inner_shell_node){
468 if(itri>=
nodes.size()) {printf(
"\nMore elements than expected in circular grid generation. Dynamically resize?"); exit(1); }
469 nodes(itri).vertex[0] = inner_shell_node;
470 nodes(itri).vertex[1] = next_inner_shell_node;
471 nodes(itri).vertex[2] = prev_outer_node;
472 inner_shell_node = next_inner_shell_node;
476 if(itri>=
nodes.size()) {printf(
"\nMore elements than expected in circular grid generation. Dynamically resize?"); exit(1); }
477 nodes(itri).vertex[0] = inner_shell_node;
478 nodes(itri).vertex[1] = prev_outer_node;
479 nodes(itri).vertex[2] = outer_node;
483 if(inodesh<nodes_in_shell){
491 if(itri!=
nodes.size()) {printf(
"\nFewer elements than expected in circular grid generation. Dynamically resize?"); exit(1); }
495 nodes(i).vertex[0]++;
496 nodes(i).vertex[1]++;
497 nodes(i).vertex[2]++;
502 int n_separatrices = 0;
503 i_x = View<int*, HostType>(
"i_x", n_xpts);
514 non_aligned_surf_idx = View<int**, CLayout, HostType>(
NoInit(
"non_aligned_surf_idx"),
num_non_aligned, 4);
524 if(nnodes<1)
return 0;
525 else return std::ceil(sqrt(12*nnodes - 3)/6 - 0.5);
539 :
nnodes(1 + 3*(nshells+1)*nshells),
544 nsurfaces(nshells+1),
563 :
nnodes(1 + 3*(nshells+1)*nshells),
568 nsurfaces(nshells+1),
589 double triangle_width = width/(2*nshells);
590 double triangle_height = height/(2*nshells);
594 double scale = std::min(width,height)/(2*nshells);
604 View<Equil::XPoint*, HostType> xpts(
NoInit(
"xpts"),
i_x.size());
605 for(
int i = 0; i<
i_x.size(); i++){
606 xpts(i).r =
x(
i_x(i)).r;
607 xpts(i).z =
x(
i_x(i)).z;
640 if(nlr.
get<
bool>(
"sml_create_analytic_grid",
false)==
true){
647 nshells = nlr.
get<
int>(
"nsurfaces", 30);
648 axis_r = nlr.
get<
double>(
"axis_r", 1.5);
651 if(
is_rank_zero()) printf(
"\nCreating a circular analytic grid with %d surfaces and %d nodes.\n", nshells, 1 + 3*(nshells+1)*nshells);
656 double scale = 0.999/nshells;
657 double rscale = scale;
658 double zscale = scale;
666 std::string node_file = nlr.
get<std::string>(
"sml_node_file",
"example_file.node");
667 std::string ele_file = nlr.
get<std::string>(
"sml_ele_file",
"example_file.ele");
668 std::string surf_file = nlr.
get<std::string>(
"sml_surf_file",
"example_file.flx.aif");
669 std::string input_file_dir = nlr.
get<std::string>(
"sml_input_file_dir",
"./");
672 node_file = input_file_dir +
"/" + node_file;
673 ele_file = input_file_dir +
"/" + ele_file;
674 surf_file = input_file_dir +
"/" + surf_file;
678 int ifile = iplane*2;
680 std::string midplane_node_file = node_file +
"." + std::to_string(ifile);
681 std::string midplane_ele_file = ele_file +
"." + std::to_string(ifile);
682 std::string midplane_surf_file = surf_file +
"." + std::to_string(ifile);
688 node_file = node_file +
"." + std::to_string(ifile);
689 ele_file = ele_file +
"." + std::to_string(ifile);
690 surf_file = surf_file +
"." + std::to_string(ifile);
Definition: grid_files.hpp:290
PlaneFiles(int nshells, double raxis, double zaxis, double rscale, double zscale, GridCreateOpt grid_opt=Hexagonal)
Definition: grid_files.hpp:538
FileRegion
Definition: grid_files.hpp:18
static constexpr GridCreateOpt Hexagonal
Definition: grid_files.hpp:615
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:44
MPI_Comm plane_comm
Definition: my_mpi.hpp:38
T get(const string ¶m, const T default_val, int val_ind=0)
Definition: NamelistReader.hpp:373
static int nshells_from_nnodes(int nnodes)
Definition: grid_files.hpp:519
static constexpr GridCreateOpt Circular
Definition: grid_files.hpp:616
int n_planes() const
Definition: grid_files.hpp:702
int my_intpl_rank
Definition: my_mpi.hpp:49
void set_fortran_flx_aif_format(int is_flx_aif_format_int)
PlaneFiles()
Definition: grid_files.hpp:65
View< int *, HostType > non_aligned_vert
Definition: grid_files.hpp:41
Definition: NamelistReader.hpp:193
Definition: magnetic_field.hpp:12
View< int *, HostType > non_aligned_nsurf
Definition: grid_files.hpp:42
int ntriangles
Definition: grid_files.hpp:29
RZBounds bounds
Simulation boundary.
Definition: magnetic_field.hpp:23
Definition: grid_files.hpp:613
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:40
void construct_circular_grid(int nshells, double raxis, double zaxis, double rscale, double zscale)
Definition: grid_files.hpp:418
View< Equil::XPoint *, HostType > get_xpts() const
Definition: grid_files.hpp:603
int nsurfaces3a
Definition: grid_files.hpp:36
View< int *, HostType > nnodes_on_surface
Definition: grid_files.hpp:38
int n_specified_xpts() const
Definition: grid_files.hpp:599
Definition: grid_files.hpp:289
GridFiles(int nshells, double raxis, double zaxis, double rscale, double zscale, GridCreateOpt grid_opt=Hexagonal)
Definition: grid_files.hpp:624
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:36
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:46
std::vector< PlaneFiles > plane_files_vec
Definition: grid_files.hpp:618
Definition: my_mpi.hpp:19
GridFiles(MagneticField< HostType > &magnetic_field, int nshells, GridCreateOpt grid_opt=Hexagonal)
Definition: grid_files.hpp:628
void use_namelist(const string &namelist)
Definition: NamelistReader.hpp:355
View< int *, HostType > i_x
Definition: grid_files.hpp:33
double max_z
Definition: rz_bounds.hpp:8
Definition: grid_files.hpp:20
double axis_r
r coordinate of axis
Definition: equil.hpp:89
double min_r
Definition: rz_bounds.hpp:5
Definition: col_grid.hpp:35
int nsurfaces1
Definition: grid_files.hpp:36
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:55
const PlaneFiles & operator[](int i) const
Definition: grid_files.hpp:621
GridCreateOpt
Definition: grid_files.hpp:288
View< int **, CLayout, HostType > non_aligned_surf_idx
Definition: grid_files.hpp:43
bool namelist_present(const string &namelist)
Definition: NamelistReader.hpp:351
int nnodes
Definition: grid_files.hpp:24
void construct_hex_grid(int nshells, double raxis, double zaxis, double rscale, double zscale)
Definition: grid_files.hpp:302
int nsurfaces3b
Definition: grid_files.hpp:36
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:562
static int nshells_from_nnodes(int nnodes)
Definition: grid_files.hpp:711
Definition: grid_files.hpp:16
double node_to_node_dist(RZPair a, RZPair b)
Definition: grid_files.hpp:405
void clear()
Definition: grid_files.hpp:706
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:26
int nsurfaces
Definition: grid_files.hpp:35
void skip(int n)
Definition: file_reader.hpp:45
View< RZPair *, HostType > x
Definition: grid_files.hpp:25
Kokkos::ViewAllocateWithoutInitializing NoInit
Definition: space_settings.hpp:68
constexpr double TWOPI
Definition: constants.hpp:9
Definition: grid_files.hpp:19
View< int *, HostType > i_surf_separatrices
Definition: grid_files.hpp:37
View< Vertex *, HostType > nodes
Definition: grid_files.hpp:30
int sum_nnodes_on_surfaces
Definition: grid_files.hpp:39
void open(const std::string &input)
Definition: file_reader.hpp:13
double max_r
Definition: rz_bounds.hpp:6