Stellarator-Tools
vmec_quantities.hpp
Go to the documentation of this file.
1 //******************************************************************************
4 //******************************************************************************
5 
6 #ifndef vmec_quantities_hpp
7 #define vmec_quantities_hpp
8 
9 #include <string>
10 #include <netcdf.h>
11 #include <array>
12 
13 #include "grid_quantity.hpp"
14 
15 //------------------------------------------------------------------------------
19 //------------------------------------------------------------------------------
21 public:
37  const double signj;
55  const double curtor;
58 
59 //------------------------------------------------------------------------------
65 //------------------------------------------------------------------------------
66  vmec_quantities(const std::string &wout_file) :
69  drds(to_prime(r)),
70  z(load_fourier<full_grid, sine> (wout_file, "zmns")),
71  dzds(to_prime(z)),
72  phipf(load<full_grid> (wout_file, "phipf")),
73  chipf(load<full_grid> (wout_file, "chipf")),
74  signj(load_scalar<int> (wout_file, "signgs")),
83  curtor(load_scalar<double> (wout_file, "ctor")),
84  p(load<full_grid> (wout_file, "presf")) {}
85 
86 //------------------------------------------------------------------------------
90 //------------------------------------------------------------------------------
92  std::vector<double> temp(100);
93 
94  for (size_t i = 0; i < 100; i++) {
95  const double s = static_cast<double> (i)/99.0;
96 
97  temp[i] = cos(s);
98  }
99 
100  return vmec_grid<full_grid> (temp);
101  }
102 
103 //------------------------------------------------------------------------------
109 //------------------------------------------------------------------------------
110  template<typename TYPE>
111  static double load_scalar(const std::string &wout_file,
112  const std::string &name) {
113  int ncid;
114  nc_open(wout_file.c_str(), NC_NOWRITE, &ncid);
115 
116  int varid;
117  TYPE temp;
118  nc_inq_varid(ncid, name.c_str(), &varid);
119  nc_get_var(ncid, varid, &temp);
120 
121  nc_close(ncid);
122 
123  return static_cast<double> (temp);
124  }
125 
126 //------------------------------------------------------------------------------
132 //------------------------------------------------------------------------------
133  template<class GIRD_CLASS>
134  static vmec_grid<GIRD_CLASS> load(const std::string &wout_file,
135  const std::string &name) {
136  int ncid;
137  nc_open(wout_file.c_str(), NC_NOWRITE, &ncid);
138 
139  int varid;
140  int ns;
141  nc_inq_varid(ncid, "ns", &varid);
142  nc_get_var(ncid, varid, &ns);
143 
144  std::vector<double> temp(ns);
145  nc_inq_varid(ncid, name.c_str(), &varid);
146  nc_get_var(ncid, varid, temp.data());
147 
148  nc_close(ncid);
149 
150  return vmec_grid<GIRD_CLASS> (temp);
151  }
152 
153 //------------------------------------------------------------------------------
159 //------------------------------------------------------------------------------
160  template<class GIRD_CLASS, class PARITY>
162  const std::string &name) {
163  int ncid;
164  nc_open(wout_file.c_str(), NC_NOWRITE, &ncid);
165 
166  int varid;
167  int ns;
168  nc_inq_varid(ncid, "ns", &varid);
169  nc_get_var(ncid, varid, &ns);
170 
171  int mnmax;
172  nc_inq_varid(ncid, "mnmax", &varid);
173  nc_get_var(ncid, varid, &mnmax);
174 
175  std::vector<double> xm(mnmax);
176  nc_inq_varid(ncid, "xm", &varid);
177  nc_get_var(ncid, varid, xm.data());
178 
179  std::vector<double> xn(mnmax);
180  nc_inq_varid(ncid, "xn", &varid);
181  nc_get_var(ncid, varid, xn.data());
182 
183  vmec_quantity<GIRD_CLASS> quantity;
184  for (size_t i = 0; i < mnmax; i++) {
185  std::vector<double> temp(ns);
186  nc_inq_varid(ncid, name.c_str(), &varid);
187 
188  const std::array<size_t, 2> start = {0, i};
189  const std::array<size_t, 2> end = {static_cast<size_t> (ns), 1};
190  nc_get_vara(ncid, varid,
191  start.data(), end.data(),
192  temp.data());
193 
194  quantity.push_back(vmec_grid<GIRD_CLASS> (temp));
195  }
196 
197  nc_close(ncid);
198 
199  return vmec_fourier<GIRD_CLASS, PARITY> (quantity, xm, xn);
200  }
201 
202 //------------------------------------------------------------------------------
208 //------------------------------------------------------------------------------
209  template<class GIRD_CLASS, class PARITY>
211  const std::string &name) {
212  int ncid;
213  nc_open(wout_file.c_str(), NC_NOWRITE, &ncid);
214 
215  int varid;
216  int ns;
217  nc_inq_varid(ncid, "ns", &varid);
218  nc_get_var(ncid, varid, &ns);
219 
220  int mnmax;
221  nc_inq_varid(ncid, "mnmax_nyq", &varid);
222  nc_get_var(ncid, varid, &mnmax);
223 
224  std::vector<double> xm(mnmax);
225  nc_inq_varid(ncid, "xm_nyq", &varid);
226  nc_get_var(ncid, varid, xm.data());
227 
228  std::vector<double> xn(mnmax);
229  nc_inq_varid(ncid, "xn_nyq", &varid);
230  nc_get_var(ncid, varid, xn.data());
231 
232  vmec_quantity<GIRD_CLASS> quantity;
233  for (size_t i = 0; i < mnmax; i++) {
234  std::vector<double> temp(ns);
235  nc_inq_varid(ncid, name.c_str(), &varid);
236 
237  const std::array<size_t, 2> start = {0, i};
238  const std::array<size_t, 2> end = {static_cast<size_t> (ns), 1};
239  nc_get_vara(ncid, varid,
240  start.data(), end.data(),
241  temp.data());
242 
243  quantity.push_back(vmec_grid<GIRD_CLASS> (temp));
244  }
245 
246  nc_close(ncid);
247 
248  return vmec_fourier<GIRD_CLASS, PARITY> (quantity, xm, xn);
249  }
250 
251 //------------------------------------------------------------------------------
256 //------------------------------------------------------------------------------
257  template<class PARITY>
259  vmec_quantity<half_grid> half_quanity;
260 
261  for (const vmec_grid<full_grid> &q : vmec.quantity) {
262  half_quanity.push_back(vmec_grid<half_grid> (q.grid.get_prime()));
263  }
264 
265  return vmec_fourier<half_grid, PARITY> (half_quanity, vmec.m, vmec.n);
266  }
267 
268 //------------------------------------------------------------------------------
273 //------------------------------------------------------------------------------
274  template<class PARITY>
276  vmec_quantity<full_grid> full_quanity;
277 
278  for (size_t i = 0, e = vmec.m.size(); i < e; i++) {
279  full_quanity.push_back(vmec_grid<full_grid> (vmec.quantity[i].grid.get_prime(vmec.m[i])));
280  }
281 
282  return vmec_fourier<full_grid, PARITY> (full_quanity, vmec.m, vmec.n);
283  }
284 
285 //------------------------------------------------------------------------------
292 //------------------------------------------------------------------------------
293  double get_r(const double s, const double u, const double v) const {
294  return r.get(s, u, v);
295  }
296 
297 //------------------------------------------------------------------------------
304 //------------------------------------------------------------------------------
305  double get_r_prime(const double s, const double u, const double v) const {
306  return 2.0*sqrt(s)*drds.get(s, u, v);
307  }
308 
309 //------------------------------------------------------------------------------
316 //------------------------------------------------------------------------------
317  double get_dr_du(const double s, const double u, const double v) const {
318  return r.get_du(s, u, v);
319  }
320 
321 //------------------------------------------------------------------------------
328 //------------------------------------------------------------------------------
329  double get_dr_dv(const double s, const double u, const double v) const {
330  return r.get_dv(s, u, v);
331  }
332 
333 //------------------------------------------------------------------------------
340 //------------------------------------------------------------------------------
341  double get_z(const double s, const double u, const double v) const {
342  return z.get(s, u, v);
343  }
344 
345 //------------------------------------------------------------------------------
352 //------------------------------------------------------------------------------
353  double get_z_prime(const double s, const double u, const double v) const {
354  return 2.0*sqrt(s)*dzds.get(s, u, v);
355  }
356 
357 //------------------------------------------------------------------------------
364 //------------------------------------------------------------------------------
365  double get_dz_du(const double s, const double u, const double v) const {
366  return z.get_du(s, u, v);
367  }
368 
369 //------------------------------------------------------------------------------
376 //------------------------------------------------------------------------------
377  double get_dz_dv(const double s, const double u, const double v) const {
378  return z.get_dv(s, u, v);
379  }
380 
381 //------------------------------------------------------------------------------
386 //------------------------------------------------------------------------------
387  double get_phipf(const double s) const {
388  return signj*sqrt(s)*phipf.get(s)/M_PI;
389  }
390 
391 //------------------------------------------------------------------------------
396 //------------------------------------------------------------------------------
397  double get_chipf(const double s) const {
398  return signj*sqrt(s)*chipf.get(s)/M_PI;
399  }
400 
401 //------------------------------------------------------------------------------
408 //------------------------------------------------------------------------------
409  double get_bsups(const double s, const double u, const double v) const {
410  return 0.0;
411  }
412 
413 //------------------------------------------------------------------------------
420 //------------------------------------------------------------------------------
421  double get_bsupu(const double s, const double u, const double v) const {
422  return bsupu.get(s, u, v);
423  }
424 
425 //------------------------------------------------------------------------------
432 //------------------------------------------------------------------------------
433  double get_bsupv(const double s, const double u, const double v) const {
434  return bsupv.get(s, u, v);
435  }
436 
437 //------------------------------------------------------------------------------
444 //------------------------------------------------------------------------------
445  double get_bsubs(const double s, const double u, const double v) const {
446  return bsubs.get(s, u, v);
447  }
448 
449 //------------------------------------------------------------------------------
456 //------------------------------------------------------------------------------
457  double get_bsubu(const double s, const double u, const double v) const {
458  return bsubu.get(s, u, v);
459  }
460 
461 //------------------------------------------------------------------------------
468 //------------------------------------------------------------------------------
469  double get_bsubv(const double s, const double u, const double v) const {
470  return bsubv.get(s, u, v);
471  }
472 
473 //------------------------------------------------------------------------------
480 //------------------------------------------------------------------------------
481  double get_jksups(const double s, const double u, const double v) const {
482  return 0.0;
483  }
484 
485 //------------------------------------------------------------------------------
492 //------------------------------------------------------------------------------
493  double get_jksupu(const double s, const double u, const double v) const {
494  return 2.0*sqrt(s)*jksupu.get(s, u, v);
495  }
496 
497 //------------------------------------------------------------------------------
504 //------------------------------------------------------------------------------
505  double get_jksupv(const double s, const double u, const double v) const {
506  return 2.0*sqrt(s)*jksupv.get(s, u, v);
507  }
508 
509 //------------------------------------------------------------------------------
516 //------------------------------------------------------------------------------
517  double get_jbsups(const double s, const double u, const double v) const {
518  return 0.0;
519  }
520 
521 //------------------------------------------------------------------------------
528 //------------------------------------------------------------------------------
529  double get_jbsupu(const double s, const double u, const double v) const {
530  return 2.0*sqrt(s)*j.get(s, u, v)*get_bsupu(s, u, v);
531  }
532 
533 //------------------------------------------------------------------------------
540 //------------------------------------------------------------------------------
541  double get_jbsupv(const double s, const double u, const double v) const {
542  return 2.0*sqrt(s)*j.get(s, u, v)*get_bsupv(s, u, v);
543  }
544 
545 //------------------------------------------------------------------------------
552 //------------------------------------------------------------------------------
553  double get_jacobian(const double s, const double u, const double v) const {
554  return 2.0*sqrt(s)*j.get(s, u, v);
555  }
556 
557 //------------------------------------------------------------------------------
561 //------------------------------------------------------------------------------
562  double get_curtor() const {
563  return curtor;
564  }
565 
566 //------------------------------------------------------------------------------
573 //------------------------------------------------------------------------------
574  double get_dbsupu_du(const double s, const double u, const double v) const {
575  return 2.0*sqrt(s)*j.get(s, u, v)*bsupu.get_du(s, u, v) + 2.0*sqrt(s)*j.get_du(s, u, v)*bsupu.get(s, u, v);
576  }
577 
578 //------------------------------------------------------------------------------
585 //------------------------------------------------------------------------------
586  double get_dbsupv_dv(const double s, const double u, const double v) const {
587  return 2.0*sqrt(s)*j.get(s, u, v)*bsupv.get_dv(s, u, v) + 2.0*sqrt(s)*j.get_dv(s, u, v)*bsupv.get(s, u, v);
588  }
589 
590 //------------------------------------------------------------------------------
597 //------------------------------------------------------------------------------
598  double get_divb(const double s, const double u, const double v) const {
599  return get_dbsupu_du(s,u,v) + get_dbsupv_dv(s,u,v);
600  }
601 
602 //------------------------------------------------------------------------------
607 //------------------------------------------------------------------------------
608  double get_pressure(const double s) const {
609  return p.get(s);
610  }
611 
612 //------------------------------------------------------------------------------
617 //------------------------------------------------------------------------------
618  double get_test(const double s) const {
619  return test_full.get(s);
620  }
621 };
622 
623 #endif
Cosine Parity function interface.
Definition: parity.hpp:62
A full grid quantity.
Definition: grid_quantity.hpp:49
half_grid get_prime() const
Get a radial derivative at a s position.
Definition: grid_quantity.cpp:59
A half grid quantity.
Definition: grid_quantity.hpp:80
Sine Parity function interface.
Definition: parity.hpp:36
double get_du(const double s, const double u, const double v) const
Get a poloidal derivative at a radial s position.
Definition: grid_quantity.hpp:241
double get(const double s, const double u, const double v) const
Get a value at a radial s position.
Definition: grid_quantity.hpp:219
double get_dv(const double s, const double u, const double v) const
Get a toroidal derivative at a radial s position.
Definition: grid_quantity.hpp:263
const GRID_CLASS grid
Radial grid buffer.
Definition: grid_quantity.hpp:115
double get(const double s) const
Get a value at a radial s position.
Definition: grid_quantity.hpp:137
Vmec quantities.
Definition: vmec_quantities.hpp:20
double get_bsubs(const double s, const double u, const double v) const
Get bsubs quantity.
Definition: vmec_quantities.hpp:445
double get_r(const double s, const double u, const double v) const
Get R quantity.
Definition: vmec_quantities.hpp:293
const vmec_fourier< full_grid, cosine > jksupv
JK^v.
Definition: vmec_quantities.hpp:53
double get_test(const double s) const
Get test quantity.
Definition: vmec_quantities.hpp:618
double get_divb(const double s, const double u, const double v) const
Get divergence of B.
Definition: vmec_quantities.hpp:598
const vmec_grid< full_grid > phipf
Radial derivative of toroidal flux.
Definition: vmec_quantities.hpp:33
vmec_quantities(const std::string &wout_file)
Vmec quantities.
Definition: vmec_quantities.hpp:66
double get_bsubu(const double s, const double u, const double v) const
Get bsubu quantity.
Definition: vmec_quantities.hpp:457
const vmec_fourier< full_grid, cosine > jksupu
JK^u.
Definition: vmec_quantities.hpp:51
static vmec_fourier< GIRD_CLASS, PARITY > load_fourier_nyq(const std::string &wout_file, const std::string &name)
Factory method to load a nyquest fourier vmec quantity.
Definition: vmec_quantities.hpp:210
const vmec_fourier< full_grid, cosine > r
r
Definition: vmec_quantities.hpp:25
double get_chipf(const double s) const
Get chipf quantity.
Definition: vmec_quantities.hpp:397
const vmec_fourier< half_grid, sine > bsubs
B_s.
Definition: vmec_quantities.hpp:45
double get_jksupv(const double s, const double u, const double v) const
Get jksupv quantity.
Definition: vmec_quantities.hpp:505
const vmec_fourier< half_grid, cosine > drds
drds
Definition: vmec_quantities.hpp:27
const vmec_grid< full_grid > p
Pressure.
Definition: vmec_quantities.hpp:57
double get_jksups(const double s, const double u, const double v) const
Get jksups quantity.
Definition: vmec_quantities.hpp:481
static double load_scalar(const std::string &wout_file, const std::string &name)
Factory method to load a scalar vmec quantity.
Definition: vmec_quantities.hpp:111
double get_bsupv(const double s, const double u, const double v) const
Get bsupv quantity.
Definition: vmec_quantities.hpp:433
double get_dbsupv_dv(const double s, const double u, const double v) const
Get dbsupvdv quantity.
Definition: vmec_quantities.hpp:586
const vmec_fourier< half_grid, cosine > bsupu
B^u.
Definition: vmec_quantities.hpp:41
static vmec_fourier< half_grid, PARITY > to_prime(const vmec_fourier< full_grid, PARITY > vmec)
Convert full grid quantity to half grid primed.
Definition: vmec_quantities.hpp:258
double get_z(const double s, const double u, const double v) const
Get Z quantity.
Definition: vmec_quantities.hpp:341
double get_dz_dv(const double s, const double u, const double v) const
Get dZdv quantity.
Definition: vmec_quantities.hpp:377
double get_jksupu(const double s, const double u, const double v) const
Get jksupu quantity.
Definition: vmec_quantities.hpp:493
const vmec_fourier< full_grid, sine > z
z
Definition: vmec_quantities.hpp:29
double get_bsupu(const double s, const double u, const double v) const
Get bsupu quantity.
Definition: vmec_quantities.hpp:421
double get_curtor() const
Get total toroidal current.
Definition: vmec_quantities.hpp:562
double get_dz_du(const double s, const double u, const double v) const
Get dZdu quantity.
Definition: vmec_quantities.hpp:365
static vmec_grid< full_grid > make_full()
Factory method to make a test quantity.
Definition: vmec_quantities.hpp:91
const vmec_grid< full_grid > test_full
test function
Definition: vmec_quantities.hpp:23
const double signj
Sign of vmec jacobian.
Definition: vmec_quantities.hpp:37
static vmec_fourier< full_grid, PARITY > to_prime(const vmec_fourier< half_grid, PARITY > vmec)
Convert full grid quantity to half grid primed.
Definition: vmec_quantities.hpp:275
static vmec_grid< GIRD_CLASS > load(const std::string &wout_file, const std::string &name)
Factory method to load a vmec quantity.
Definition: vmec_quantities.hpp:134
double get_bsups(const double s, const double u, const double v) const
Get bsups quantity.
Definition: vmec_quantities.hpp:409
double get_phipf(const double s) const
Get phipf quantity.
Definition: vmec_quantities.hpp:387
const vmec_grid< full_grid > chipf
Radial derivative of poloidal flux.
Definition: vmec_quantities.hpp:35
const vmec_fourier< half_grid, cosine > j
Jacobian.
Definition: vmec_quantities.hpp:39
double get_jbsupv(const double s, const double u, const double v) const
Get jbsupv quantity.
Definition: vmec_quantities.hpp:541
double get_r_prime(const double s, const double u, const double v) const
Get R' quantity.
Definition: vmec_quantities.hpp:305
double get_dbsupu_du(const double s, const double u, const double v) const
Get dbsupudu quantity.
Definition: vmec_quantities.hpp:574
double get_jbsups(const double s, const double u, const double v) const
Get jbsups quantity.
Definition: vmec_quantities.hpp:517
double get_dr_du(const double s, const double u, const double v) const
Get dRdu quantity.
Definition: vmec_quantities.hpp:317
const vmec_fourier< half_grid, sine > dzds
dzds
Definition: vmec_quantities.hpp:31
double get_jbsupu(const double s, const double u, const double v) const
Get jbsupu quantity.
Definition: vmec_quantities.hpp:529
double get_dr_dv(const double s, const double u, const double v) const
Get dRdv quantity.
Definition: vmec_quantities.hpp:329
double get_jacobian(const double s, const double u, const double v) const
Get jacobian.
Definition: vmec_quantities.hpp:553
double get_pressure(const double s) const
Get pressure quantity.
Definition: vmec_quantities.hpp:608
const double curtor
Total toroidal current.
Definition: vmec_quantities.hpp:55
const vmec_fourier< half_grid, cosine > bsubv
B_v.
Definition: vmec_quantities.hpp:49
double get_bsubv(const double s, const double u, const double v) const
Get bsubv quantity.
Definition: vmec_quantities.hpp:469
const vmec_fourier< half_grid, cosine > bsupv
B^v.
Definition: vmec_quantities.hpp:43
static vmec_fourier< GIRD_CLASS, PARITY > load_fourier(const std::string &wout_file, const std::string &name)
Factory method to load a fourier vmec quantity.
Definition: vmec_quantities.hpp:161
const vmec_fourier< half_grid, cosine > bsubu
B_u.
Definition: vmec_quantities.hpp:47
double get_z_prime(const double s, const double u, const double v) const
Get Z' quantity.
Definition: vmec_quantities.hpp:353
Contains classes to interpolate full and half grid quanities.
std::vector< vmec_grid< GRID_CLASS > > vmec_quantity
Type for vmec fourier quantities.
Definition: grid_quantity.hpp:183
character(len=siesta_namelist_name_length) wout_file
Filename of the VMEC woutfile.
Definition: siesta_namelist.f90:282