XGCa
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
magnetic_equil_files.hpp
Go to the documentation of this file.
1 #ifndef MAGNETIC_EQUIL_FILES_HPP
2 #define MAGNETIC_EQUIL_FILES_HPP
3 
4 #include "space_settings.hpp"
5 #include "NamelistReader.hpp"
6 #include "globals.hpp"
7 #include "range_view.hpp"
8 #include "rz_bounds.hpp"
9 #include "file_reader.hpp"
10 
11 extern "C" void xgc_read_eq_m3dc1_dims(int* mpsi, int* mr, int* mz, double* min_r, double* max_r, double* min_z, double* max_z);
12 extern "C" void xgc_read_eq_m3dc1(int mpsi, int mr, int mz, double min_r, double max_r, double min_z, double max_z,
13  double* axis_r, double* axis_z, double* axis_b, double* x_psi, double* x_r, double* x_z,
14  double* psi_grid, double* rgrid, double* zgrid, double* I, double* psi_rz);
15 
17  int mpsi;
18  int mr;
19  int mz;
21  double axis_r;
22  double axis_z;
23  double axis_b;
24  double x_psi;
25  double x_r;
26  double x_z;
27  View<double*, CLayout, HostType> psi_grid;
28  View<double*, HostType> rgrid;
29  View<double*, HostType> zgrid;
30  View<double*, CLayout, HostType> I;
31  View<double**, CLayout, HostType> psi_rz;
32 
33  // Reads and communicates eq dimensions
34  void read_eq_dims(FileReader& file_reader){
36  file_reader.skip_line(); // Skip first line (Header)
37  file_reader.read(&mr, 1);
38  file_reader.read(&mz, 1);
39  file_reader.read(&mpsi, 1);
40  file_reader.read(&bounds.min_r, 1);
41  file_reader.read(&bounds.max_r, 1);
42  file_reader.read(&bounds.min_z, 1);
43  file_reader.read(&bounds.max_z, 1);
44  }
45 
46  void read_eq(FileReader& file_reader){
47  file_reader.read(&axis_r, 1);
48  file_reader.read(&axis_z, 1);
49  file_reader.read(&axis_b, 1);
50  file_reader.read(&x_psi, 1);
51  file_reader.read(&x_r, 1);
52  file_reader.read(&x_z, 1);
53 
54  file_reader.read(psi_grid.data(), psi_grid.size());
55  file_reader.read(I.data(), I.size());
56  file_reader.read(psi_rz.data(), psi_rz.size());
57 
58  int end_flag;
59  file_reader.read(&end_flag, 1);
60  if(end_flag!=-1){
61  exit_XGC("\nError reading eqd file: Did not encounter end flag (-1) where expected. The file format may be incorrect.\n");
62  }
63  }
64 
65  // Read the files
66  MagneticEquilFiles(bool read_m3dc1, const std::string& eq_filename){
67  // file_reader if using simple eqd file
68  FileReader file_reader;
69 
70  // First read in the dimensions
71  if(is_rank_zero()){
72  if(read_m3dc1){
74  }else{
75  file_reader.open(eq_filename);
76  read_eq_dims(file_reader);
77  }
78  }
79 
80  // Broadcast dimensions to other processes
81 #ifdef USE_MPI
82  int root_rank = 0;
83  MPI_Bcast(&mr, 1, MPI_INT, root_rank, SML_COMM_WORLD);
84  MPI_Bcast(&mz, 1, MPI_INT, root_rank, SML_COMM_WORLD);
85  MPI_Bcast(&mpsi, 1, MPI_INT, root_rank, SML_COMM_WORLD);
86 
87  MPI_Bcast(&bounds, 4, MPI_DOUBLE, root_rank, SML_COMM_WORLD);
88 #endif
89 
90  // Allocate views
91  psi_grid = View<double*, CLayout, HostType>("psi_grid", mpsi);
92  I = View<double*, CLayout, HostType>("I", mpsi);
93  psi_rz = View<double**, CLayout, HostType>("psi_rz", mz, mr);
94 
95  // Set rgrid and zgrid to be linear
98 
99  // Read in data
100  if(is_rank_zero()){
101  if(read_m3dc1){
103  psi_grid.data(), rgrid.data(), zgrid.data(), I.data(), psi_rz.data());
104  }else{
105  read_eq(file_reader);
106  }
107  }
108 
109  // Broadcast values to other processes
110 #ifdef USE_MPI
111  MPI_Bcast(&axis_r,1, MPI_DOUBLE, root_rank, SML_COMM_WORLD);
112  MPI_Bcast(&axis_z,1, MPI_DOUBLE, root_rank, SML_COMM_WORLD);
113  MPI_Bcast(&axis_b,1, MPI_DOUBLE, root_rank, SML_COMM_WORLD);
114 
115  MPI_Bcast(&x_psi, 1, MPI_DOUBLE, root_rank, SML_COMM_WORLD);
116  MPI_Bcast(&x_r, 1, MPI_DOUBLE, root_rank, SML_COMM_WORLD);
117  MPI_Bcast(&x_z, 1, MPI_DOUBLE, root_rank, SML_COMM_WORLD);
118 
119  MPI_Bcast(psi_grid.data(), psi_grid.size(), MPI_DOUBLE, root_rank, SML_COMM_WORLD);
120  MPI_Bcast(I.data(), I.size(), MPI_DOUBLE, root_rank, SML_COMM_WORLD);
121  MPI_Bcast(psi_rz.data(), psi_rz.size(), MPI_DOUBLE, root_rank, SML_COMM_WORLD);
122 #endif
123  }
124 
125  // Analytic constructor
127 
128  if(is_rank_zero()) printf("\nCreating a circular analytic magnetic equilibrium.\n");
129 
130  axis_r = 1.5;
131  if(nlr.namelist_present("analytic_grid_param")){
132  nlr.use_namelist("analytic_grid_param");
133  axis_r = nlr.get<double>("axis_r", 1.5);
134  }
135 
136  axis_z = 0.0;
137  axis_b = 1.0;
138 
139  //nlr.use_namelist("analytic_eq_files_param");
140  mpsi=100;
141  mr=100;
142  mz=100;
143  bounds.min_r = axis_r - 1.0;
144  bounds.max_r = axis_r + 1.0;
145  bounds.min_z = axis_z - 1.0;
146  bounds.max_z = axis_z + 1.0;
147  x_psi = 2.0; // outside bounds
148  x_r = axis_r;
149  x_z = bounds.min_z - 1.0; // outside bounds
150 
151  // Set rgrid and zgrid to be linear
154  psi_grid = create_range_view("psi_grid",0.0, 1.0, mpsi);
155 
156  // Allocate views
157  I = View<double*, CLayout, HostType>("I", mpsi);
158  Kokkos::deep_copy(I, 10000.0); // High safety factor for now - make input later
159 
160  // Set up psi_grid - concentric linear psi
161  psi_rz = View<double**, CLayout, HostType>("psi_rz", mz, mr);
162  for (int j=0; j<zgrid.size(); j++){
163  for (int i=0; i<rgrid.size(); i++){
164  double rp = (rgrid(i) - axis_r);
165  double zp = (zgrid(j) - axis_z);
166  psi_rz(j,i) = sqrt(rp*rp + zp*zp);
167  }
168  }
169  }
170 
171  /* Returns MagneticEquilFiles object using the file-reader constructor. This may be updated so the input file
172  * can specify an analytic grid.
173  * */
175  nlr.use_namelist("sml_param");
176 
177  if(nlr.get<bool>("sml_create_analytic_equil", false)==true) return MagneticEquilFiles(nlr);
178 
179  // Invert psi values if appropriate
180  bool bt_sign = nlr.get<double>("sml_bt_sign", -1.0);
181  bool bp_sign = nlr.get<double>("sml_bp_sign", 1.0);
182 
183  // Read file names and path from input file
184  bool read_eqd = nlr.get<bool>("sml_read_eq", true);
185  bool read_m3dc1 = nlr.get<bool>("sml_read_eq_m3dc1", false);
186  // Ignore read_eqd if read_m3dc1 is true
187 
188  if(!(read_eqd || read_m3dc1)) exit_XGC("Error: must use either sml_read_eq or sml_read_eq_m3dc1");
189 
190  std::string input_file_dir = nlr.get<std::string>("sml_input_file_dir","./");
191 
192  nlr.use_namelist("eq_param");
193  std::string eq_filename = nlr.get<std::string>("eq_filename","example_file.eqd");
194 
195  // Add path to file name
196  eq_filename = input_file_dir + "/" + eq_filename;
197 
198  // Read file
199  auto eq_file = MagneticEquilFiles(read_m3dc1, eq_filename);
200 
201  eq_file.x_psi *= bp_sign;
202  eq_file.axis_b *= bt_sign;
203 
204  for(int i=0; i<eq_file.psi_grid.size(); i++) eq_file.psi_grid(i) *= bp_sign;
205  for(int i=0; i<eq_file.I.size(); i++) eq_file.I(i) *= bt_sign;
206  for(int i=0; i<eq_file.psi_rz.extent(0); i++)
207  for(int j=0; j<eq_file.psi_rz.extent(1); j++)
208  eq_file.psi_rz(i,j) *= bp_sign;
209 
210  // Return result
211  return eq_file;
212  }
213 
214 };
215 
216 #endif
View< double *, HostType > rgrid
Definition: magnetic_equil_files.hpp:28
double axis_b
Definition: magnetic_equil_files.hpp:23
Definition: magnetic_equil_files.hpp:16
bool is_rank_zero()
Definition: globals.hpp:27
MPI_Comm SML_COMM_WORLD
Definition: my_mpi.cpp:4
double x_psi
Definition: magnetic_equil_files.hpp:24
T get(const string &param, const T default_val, int val_ind=0)
Definition: NamelistReader.hpp:373
double x_z
Definition: magnetic_equil_files.hpp:26
Definition: rz_bounds.hpp:4
View< double *, HostType > create_range_view(std::string name, double x_min, double x_max, int n)
Definition: range_view.hpp:6
Definition: NamelistReader.hpp:193
void read_eq_dims(FileReader &file_reader)
Definition: magnetic_equil_files.hpp:34
void read(T *var, int n)
Definition: file_reader.hpp:28
double axis_z
Definition: magnetic_equil_files.hpp:22
View< double *, HostType > zgrid
Definition: magnetic_equil_files.hpp:29
MagneticEquilFiles(NLReader::NamelistReader &nlr)
Definition: magnetic_equil_files.hpp:126
double x_r
Definition: magnetic_equil_files.hpp:25
int mr
Definition: magnetic_equil_files.hpp:18
void use_namelist(const string &namelist)
Definition: NamelistReader.hpp:355
double max_z
Definition: rz_bounds.hpp:8
View< double **, CLayout, HostType > psi_rz
Definition: magnetic_equil_files.hpp:31
subroutine xgc_read_eq_m3dc1_dims(mpsi, mr, mz, min_r, max_r, min_z, max_z)
Definition: read.F90:22
double axis_r
Definition: magnetic_equil_files.hpp:21
MagneticEquilFiles(bool read_m3dc1, const std::string &eq_filename)
Definition: magnetic_equil_files.hpp:66
double min_r
Definition: rz_bounds.hpp:5
int mpsi
Definition: magnetic_equil_files.hpp:17
void read_eq(FileReader &file_reader)
Definition: magnetic_equil_files.hpp:46
Definition: file_reader.hpp:6
static MagneticEquilFiles eq_files_from_namelist(NLReader::NamelistReader &nlr)
Definition: magnetic_equil_files.hpp:174
subroutine xgc_read_eq_m3dc1(mpsi, mr, mz, min_r, max_r, min_z, max_z, axis_r, axis_z, axis_b, x_psi, x_r, x_z, psi_grid, rgrid, zgrid, eq_I, psi_rz)
Definition: read.F90:54
bool namelist_present(const string &namelist)
Definition: NamelistReader.hpp:351
View< double *, CLayout, HostType > I
Definition: magnetic_equil_files.hpp:30
void exit_XGC(std::string msg)
Definition: globals.hpp:37
View< double *, CLayout, HostType > psi_grid
Definition: magnetic_equil_files.hpp:27
RZBounds bounds
Definition: magnetic_equil_files.hpp:20
double min_z
Definition: rz_bounds.hpp:7
void skip_line()
Definition: file_reader.hpp:62
int mz
Definition: magnetic_equil_files.hpp:19
void open(const std::string &input)
Definition: file_reader.hpp:13
double max_r
Definition: rz_bounds.hpp:6