Graph Framework
Loading...
Searching...
No Matches
equilibrium.hpp
Go to the documentation of this file.
1//------------------------------------------------------------------------------
6//------------------------------------------------------------------------------
207//------------------------------------------------------------------------------
208
209#ifndef equilibrium_h
210#define equilibrium_h
211
212#include <mutex>
213
214#include <netcdf.h>
215
216#include "vector.hpp"
217#include "trigonometry.hpp"
218#include "piecewise.hpp"
219#include "math.hpp"
220#include "arithmetic.hpp"
221#include "newton.hpp"
222#include "output.hpp"
223
225namespace equilibrium {
226//******************************************************************************
227// Equilibrium interface
228//******************************************************************************
229//------------------------------------------------------------------------------
234//------------------------------------------------------------------------------
235 template<jit::float_scalar T, bool SAFE_MATH=false>
236 class generic {
237 protected:
239 const std::vector<T> ion_masses;
241 const std::vector<uint8_t> ion_charges;
242
243 public:
244//------------------------------------------------------------------------------
249//------------------------------------------------------------------------------
250 generic(const std::vector<T> &masses,
251 const std::vector<uint8_t> &charges) :
252 ion_masses(masses), ion_charges(charges) {
253 assert(ion_masses.size() == ion_charges.size() &&
254 "Masses and charges need the same number of elements.");
255 }
256
257//------------------------------------------------------------------------------
259//------------------------------------------------------------------------------
260 virtual ~generic() {}
261
262//------------------------------------------------------------------------------
266//------------------------------------------------------------------------------
267 size_t get_num_ion_species() const {
268 return ion_masses.size();
269 }
270
271//------------------------------------------------------------------------------
276//------------------------------------------------------------------------------
277 T get_ion_mass(const size_t index) const {
278 return ion_masses.at(index);
279 }
280
281//------------------------------------------------------------------------------
286//------------------------------------------------------------------------------
287 uint8_t get_ion_charge(const size_t index) const {
288 return ion_charges.at(index);
289 }
290
291//------------------------------------------------------------------------------
298//------------------------------------------------------------------------------
299 virtual graph::shared_leaf<T, SAFE_MATH>
300 get_electron_density(graph::shared_leaf<T, SAFE_MATH> x,
301 graph::shared_leaf<T, SAFE_MATH> y,
302 graph::shared_leaf<T, SAFE_MATH> z) = 0;
303
304//------------------------------------------------------------------------------
312//------------------------------------------------------------------------------
313 virtual graph::shared_leaf<T, SAFE_MATH>
314 get_ion_density(const size_t index,
315 graph::shared_leaf<T, SAFE_MATH> x,
316 graph::shared_leaf<T, SAFE_MATH> y,
317 graph::shared_leaf<T, SAFE_MATH> z) = 0;
318
319//------------------------------------------------------------------------------
326//------------------------------------------------------------------------------
327 virtual graph::shared_leaf<T, SAFE_MATH>
328 get_electron_temperature(graph::shared_leaf<T, SAFE_MATH> x,
329 graph::shared_leaf<T, SAFE_MATH> y,
330 graph::shared_leaf<T, SAFE_MATH> z) = 0;
331
332//------------------------------------------------------------------------------
340//------------------------------------------------------------------------------
341 virtual graph::shared_leaf<T, SAFE_MATH>
342 get_ion_temperature(const size_t index,
343 graph::shared_leaf<T, SAFE_MATH> x,
344 graph::shared_leaf<T, SAFE_MATH> y,
345 graph::shared_leaf<T, SAFE_MATH> z) = 0;
346
347//------------------------------------------------------------------------------
354//------------------------------------------------------------------------------
355 virtual graph::shared_vector<T, SAFE_MATH>
356 get_magnetic_field(graph::shared_leaf<T, SAFE_MATH> x,
357 graph::shared_leaf<T, SAFE_MATH> y,
358 graph::shared_leaf<T, SAFE_MATH> z) = 0;
359
360//------------------------------------------------------------------------------
367//------------------------------------------------------------------------------
368 virtual graph::shared_leaf<T, SAFE_MATH>
369 get_characteristic_field(const size_t device_number=0) = 0;
370
371//------------------------------------------------------------------------------
378//------------------------------------------------------------------------------
379 virtual graph::shared_vector<T, SAFE_MATH>
387
388//------------------------------------------------------------------------------
395//------------------------------------------------------------------------------
396 virtual graph::shared_vector<T, SAFE_MATH>
404
405//------------------------------------------------------------------------------
412//------------------------------------------------------------------------------
413 virtual graph::shared_vector<T, SAFE_MATH>
421
422//------------------------------------------------------------------------------
429//------------------------------------------------------------------------------
430 virtual graph::shared_leaf<T, SAFE_MATH>
436
437//------------------------------------------------------------------------------
444//------------------------------------------------------------------------------
445 virtual graph::shared_leaf<T, SAFE_MATH>
451
452//------------------------------------------------------------------------------
459//------------------------------------------------------------------------------
460 virtual graph::shared_leaf<T, SAFE_MATH>
466 };
467
469 template<jit::float_scalar T, bool SAFE_MATH=false>
470 using shared = std::shared_ptr<generic<T, SAFE_MATH>>;
471
472//******************************************************************************
473// No Magnetic equilibrium.
474//******************************************************************************
475//------------------------------------------------------------------------------
480//------------------------------------------------------------------------------
481 template<jit::float_scalar T, bool SAFE_MATH=false>
482 class no_magnetic_field : public generic<T, SAFE_MATH> {
483 public:
484//------------------------------------------------------------------------------
486//------------------------------------------------------------------------------
487 no_magnetic_field() :
488 generic<T, SAFE_MATH> ({3.34449469E-27}, {1}) {}
489
490//------------------------------------------------------------------------------
497//------------------------------------------------------------------------------
498 virtual graph::shared_leaf<T, SAFE_MATH>
506
507//------------------------------------------------------------------------------
515//------------------------------------------------------------------------------
516 virtual graph::shared_leaf<T, SAFE_MATH>
517 get_ion_density(const size_t index,
518 graph::shared_leaf<T, SAFE_MATH> x,
519 graph::shared_leaf<T, SAFE_MATH> y,
520 graph::shared_leaf<T, SAFE_MATH> z) final {
521 return graph::constant<T, SAFE_MATH> (static_cast<T> (1.0E19)) *
522 (graph::constant<T, SAFE_MATH> (static_cast<T> (0.1))*x +
523 graph::one<T, SAFE_MATH> ());
524 }
525
526//------------------------------------------------------------------------------
533//------------------------------------------------------------------------------
534 virtual graph::shared_leaf<T, SAFE_MATH>
540
541//------------------------------------------------------------------------------
549//------------------------------------------------------------------------------
550 virtual graph::shared_leaf<T, SAFE_MATH>
551 get_ion_temperature(const size_t index,
552 graph::shared_leaf<T, SAFE_MATH> x,
553 graph::shared_leaf<T, SAFE_MATH> y,
554 graph::shared_leaf<T, SAFE_MATH> z) final {
555 return graph::constant<T, SAFE_MATH> (static_cast<T> (1000.0));
556 }
557
558//------------------------------------------------------------------------------
565//------------------------------------------------------------------------------
566 virtual graph::shared_vector<T, SAFE_MATH>
573
574//------------------------------------------------------------------------------
581//------------------------------------------------------------------------------
582 virtual graph::shared_leaf<T, SAFE_MATH>
583 get_characteristic_field(const size_t device_number=0) final {
584 return graph::one<T, SAFE_MATH> ();
585 }
586 };
587
588//------------------------------------------------------------------------------
595//------------------------------------------------------------------------------
596 template<jit::float_scalar T, bool SAFE_MATH=false>
597 shared<T, SAFE_MATH> make_no_magnetic_field() {
598 return std::make_shared<no_magnetic_field<T, SAFE_MATH>> ();
599 }
600
601//******************************************************************************
602// Slab equilibrium.
603//******************************************************************************
604//------------------------------------------------------------------------------
609//------------------------------------------------------------------------------
610 template<jit::float_scalar T, bool SAFE_MATH=false>
611 class slab : public generic<T, SAFE_MATH> {
612 public:
613//------------------------------------------------------------------------------
615//------------------------------------------------------------------------------
616 slab() :
617 generic<T, SAFE_MATH> ({3.34449469E-27}, {1}) {}
618
619//------------------------------------------------------------------------------
626//------------------------------------------------------------------------------
627 virtual graph::shared_leaf<T, SAFE_MATH>
633
634//------------------------------------------------------------------------------
642//------------------------------------------------------------------------------
643 virtual graph::shared_leaf<T, SAFE_MATH>
644 get_ion_density(const size_t index,
645 graph::shared_leaf<T, SAFE_MATH> x,
646 graph::shared_leaf<T, SAFE_MATH> y,
647 graph::shared_leaf<T, SAFE_MATH> z) final {
648 return graph::constant<T, SAFE_MATH> (static_cast<T> (1.0E19));
649 }
650
651//------------------------------------------------------------------------------
658//------------------------------------------------------------------------------
659 virtual graph::shared_leaf<T, SAFE_MATH>
665
666//------------------------------------------------------------------------------
674//------------------------------------------------------------------------------
675 virtual graph::shared_leaf<T, SAFE_MATH>
676 get_ion_temperature(const size_t index,
677 graph::shared_leaf<T, SAFE_MATH> x,
678 graph::shared_leaf<T, SAFE_MATH> y,
679 graph::shared_leaf<T, SAFE_MATH> z) final {
680 return graph::constant<T, SAFE_MATH> (static_cast<T> (1000.0));
681 }
682
683//------------------------------------------------------------------------------
690//------------------------------------------------------------------------------
691 virtual graph::shared_vector<T, SAFE_MATH>
697
698//------------------------------------------------------------------------------
705//------------------------------------------------------------------------------
706 virtual graph::shared_leaf<T, SAFE_MATH>
707 get_characteristic_field(const size_t device_number=0) final {
708 return graph::one<T, SAFE_MATH> ();
709 }
710 };
711
712//------------------------------------------------------------------------------
719//------------------------------------------------------------------------------
720 template<jit::float_scalar T, bool SAFE_MATH=false>
721 shared<T, SAFE_MATH> make_slab() {
722 return std::make_shared<slab<T, SAFE_MATH>> ();
723 }
724
725//******************************************************************************
726// Slab density equilibrium.
727//******************************************************************************
728//------------------------------------------------------------------------------
733//------------------------------------------------------------------------------
734 template<jit::float_scalar T, bool SAFE_MATH=false>
735 class slab_density : public generic<T, SAFE_MATH> {
736 public:
737//------------------------------------------------------------------------------
739//------------------------------------------------------------------------------
740 slab_density() :
741 generic<T, SAFE_MATH> ({3.34449469E-27}, {1}) {}
742
743//------------------------------------------------------------------------------
750//------------------------------------------------------------------------------
751 virtual graph::shared_leaf<T, SAFE_MATH>
759
760//------------------------------------------------------------------------------
768//------------------------------------------------------------------------------
769 virtual graph::shared_leaf<T, SAFE_MATH>
770 get_ion_density(const size_t index,
771 graph::shared_leaf<T, SAFE_MATH> x,
772 graph::shared_leaf<T, SAFE_MATH> y,
773 graph::shared_leaf<T, SAFE_MATH> z) final {
774 return graph::constant<T, SAFE_MATH> (static_cast<T> (1.0E19)) *
775 (graph::constant<T, SAFE_MATH> (static_cast<T> (0.1))*x +
776 graph::one<T, SAFE_MATH> ());
777 }
778
779//------------------------------------------------------------------------------
786//------------------------------------------------------------------------------
787 virtual graph::shared_leaf<T, SAFE_MATH>
793
794//------------------------------------------------------------------------------
802//------------------------------------------------------------------------------
803 virtual graph::shared_leaf<T, SAFE_MATH>
804 get_ion_temperature(const size_t index,
805 graph::shared_leaf<T, SAFE_MATH> x,
806 graph::shared_leaf<T, SAFE_MATH> y,
807 graph::shared_leaf<T, SAFE_MATH> z) final {
808 return graph::constant<T, SAFE_MATH> (static_cast<T> (1000.0));
809 }
810
811//------------------------------------------------------------------------------
818//------------------------------------------------------------------------------
819 virtual graph::shared_vector<T, SAFE_MATH>
826
827//------------------------------------------------------------------------------
834//------------------------------------------------------------------------------
835 virtual graph::shared_leaf<T, SAFE_MATH>
836 get_characteristic_field(const size_t device_number=0) final {
837 return graph::one<T, SAFE_MATH> ();
838 }
839 };
840
841//------------------------------------------------------------------------------
848//------------------------------------------------------------------------------
849 template<jit::float_scalar T, bool SAFE_MATH=false>
850 shared<T, SAFE_MATH> make_slab_density() {
851 return std::make_shared<slab_density<T, SAFE_MATH>> ();
852 }
853
854//******************************************************************************
855// Slab field gradient equilibrium.
856//******************************************************************************
857//------------------------------------------------------------------------------
862//------------------------------------------------------------------------------
863 template<jit::float_scalar T, bool SAFE_MATH=false>
864 class slab_field : public generic<T, SAFE_MATH> {
865 public:
866//------------------------------------------------------------------------------
868//------------------------------------------------------------------------------
869 slab_field() :
870 generic<T, SAFE_MATH> ({3.34449469E-27}, {1}) {}
871
872//------------------------------------------------------------------------------
879//------------------------------------------------------------------------------
880 virtual graph::shared_leaf<T, SAFE_MATH>
888
889//------------------------------------------------------------------------------
897//------------------------------------------------------------------------------
898 virtual graph::shared_leaf<T, SAFE_MATH>
905
906//------------------------------------------------------------------------------
913//------------------------------------------------------------------------------
914 virtual graph::shared_leaf<T, SAFE_MATH>
922
923//------------------------------------------------------------------------------
931//------------------------------------------------------------------------------
932 virtual graph::shared_leaf<T, SAFE_MATH>
939
940//------------------------------------------------------------------------------
947//------------------------------------------------------------------------------
948 virtual graph::shared_vector<T, SAFE_MATH>
954
955//------------------------------------------------------------------------------
962//------------------------------------------------------------------------------
963 virtual graph::shared_leaf<T, SAFE_MATH>
964 get_characteristic_field(const size_t device_number=0) final {
965 return graph::one<T, SAFE_MATH> ();
966 }
967 };
968
969//------------------------------------------------------------------------------
976//------------------------------------------------------------------------------
977 template<jit::float_scalar T, bool SAFE_MATH=false>
978 shared<T, SAFE_MATH> make_slab_field() {
979 return std::make_shared<slab_field<T, SAFE_MATH>> ();
980 }
981//******************************************************************************
982// Gaussian density with a uniform magnetic field.
983//******************************************************************************
984//------------------------------------------------------------------------------
989//------------------------------------------------------------------------------
990 template<jit::float_scalar T, bool SAFE_MATH=false>
991 class gaussian_density : public generic<T, SAFE_MATH> {
992 public:
993//------------------------------------------------------------------------------
995//------------------------------------------------------------------------------
996 gaussian_density() :
997 generic<T, SAFE_MATH> ({3.34449469E-27}, {1}) {}
998
999//------------------------------------------------------------------------------
1006//------------------------------------------------------------------------------
1007 virtual graph::shared_leaf<T, SAFE_MATH>
1014
1015//------------------------------------------------------------------------------
1023//------------------------------------------------------------------------------
1024 virtual graph::shared_leaf<T, SAFE_MATH>
1025 get_ion_density(const size_t index,
1026 graph::shared_leaf<T, SAFE_MATH> x,
1027 graph::shared_leaf<T, SAFE_MATH> y,
1028 graph::shared_leaf<T, SAFE_MATH> z) final {
1029 return graph::constant<T, SAFE_MATH> (static_cast<T> (1.0E19)) *
1030 graph::exp((x*x + y*y)/graph::constant<T, SAFE_MATH> (static_cast<T> (-0.2)));
1031 }
1032
1033//------------------------------------------------------------------------------
1040//------------------------------------------------------------------------------
1041 virtual graph::shared_leaf<T, SAFE_MATH>
1047
1048//------------------------------------------------------------------------------
1056//------------------------------------------------------------------------------
1057 virtual graph::shared_leaf<T, SAFE_MATH>
1058 get_ion_temperature(const size_t index,
1059 graph::shared_leaf<T, SAFE_MATH> x,
1060 graph::shared_leaf<T, SAFE_MATH> y,
1061 graph::shared_leaf<T, SAFE_MATH> z) final {
1062 return graph::constant<T, SAFE_MATH> (static_cast<T> (1000.0));
1063 }
1064
1065//------------------------------------------------------------------------------
1072//------------------------------------------------------------------------------
1073 virtual graph::shared_vector<T, SAFE_MATH>
1080
1081//------------------------------------------------------------------------------
1088//------------------------------------------------------------------------------
1089 virtual graph::shared_leaf<T, SAFE_MATH>
1090 get_characteristic_field(const size_t device_number=0) final {
1091 return graph::one<T, SAFE_MATH> ();
1092 }
1093 };
1094
1095//------------------------------------------------------------------------------
1102//------------------------------------------------------------------------------
1103 template<jit::float_scalar T, bool SAFE_MATH=false>
1104 shared<T, SAFE_MATH> make_gaussian_density() {
1105 return std::make_shared<gaussian_density<T, SAFE_MATH>> ();
1106 }
1107
1108//------------------------------------------------------------------------------
1119//------------------------------------------------------------------------------
1120 template<jit::float_scalar T, bool SAFE_MATH=false>
1121 graph::shared_leaf<T, SAFE_MATH> build_1D_spline(graph::output_nodes<T, SAFE_MATH> c,
1122 graph::shared_leaf<T, SAFE_MATH> x,
1123 const T scale,
1124 const T offset) {
1125 auto c3 = c[3]/(scale*scale*scale);
1126 auto c2 = c[2]/(scale*scale) - static_cast<T> (3.0)*offset*c[3]/(scale*scale*scale);
1127 auto c1 = c[1]/scale - static_cast<T> (2.0)*offset*c[2]/(scale*scale) + static_cast<T> (3.0)*offset*offset*c[3]/(scale*scale*scale);
1128 auto c0 = c[0] - offset*c[1]/scale + offset*offset*c[2]/(scale*scale) - offset*offset*offset*c[3]/(scale*scale*scale);
1129
1130 return graph::fma(graph::fma(graph::fma(c3, x, c2), x, c1), x, c0);
1131 }
1132
1133//******************************************************************************
1134// 2D EFIT equilibrium.
1135//******************************************************************************
1136//------------------------------------------------------------------------------
1144//------------------------------------------------------------------------------
1145 template<jit::float_scalar T, bool SAFE_MATH=false>
1146 class efit final : public generic<T, SAFE_MATH> {
1147 private:
1149 const T psimin;
1151 const T dpsi;
1152
1153// Temperature spline coefficients.
1155 const backend::buffer<T> te_c0;
1157 const backend::buffer<T> te_c1;
1159 const backend::buffer<T> te_c2;
1161 const backend::buffer<T> te_c3;
1163 graph::shared_leaf<T, SAFE_MATH> te_scale;
1164
1165// Density spline coefficients.
1167 const backend::buffer<T> ne_c0;
1169 const backend::buffer<T> ne_c1;
1171 const backend::buffer<T> ne_c2;
1173 const backend::buffer<T> ne_c3;
1175 graph::shared_leaf<T, SAFE_MATH> ne_scale;
1176
1177// Pressure spline coefficients.
1179 const backend::buffer<T> pres_c0;
1181 const backend::buffer<T> pres_c1;
1183 const backend::buffer<T> pres_c2;
1185 const backend::buffer<T> pres_c3;
1187 graph::shared_leaf<T, SAFE_MATH> pres_scale;
1188
1190 const T rmin;
1192 const T dr;
1194 const T zmin;
1196 const T dz;
1197
1198// Fpol spline coefficients.
1200 const backend::buffer<T> fpol_c0;
1202 const backend::buffer<T> fpol_c1;
1204 const backend::buffer<T> fpol_c2;
1206 const backend::buffer<T> fpol_c3;
1207
1208// Psi spline coefficients.
1210 const size_t num_cols;
1212 const backend::buffer<T> c00;
1214 const backend::buffer<T> c01;
1216 const backend::buffer<T> c02;
1218 const backend::buffer<T> c03;
1220 const backend::buffer<T> c10;
1222 const backend::buffer<T> c11;
1224 const backend::buffer<T> c12;
1226 const backend::buffer<T> c13;
1228 const backend::buffer<T> c20;
1230 const backend::buffer<T> c21;
1232 const backend::buffer<T> c22;
1234 const backend::buffer<T> c23;
1236 const backend::buffer<T> c30;
1238 const backend::buffer<T> c31;
1240 const backend::buffer<T> c32;
1242 const backend::buffer<T> c33;
1243
1244// Cached values.
1246 graph::shared_leaf<T, SAFE_MATH> x_cache;
1248 graph::shared_leaf<T, SAFE_MATH> y_cache;
1250 graph::shared_leaf<T, SAFE_MATH> z_cache;
1251
1253 graph::shared_leaf<T, SAFE_MATH> ne_cache;
1255 graph::shared_leaf<T, SAFE_MATH> ni_cache;
1257 graph::shared_leaf<T, SAFE_MATH> te_cache;
1259 graph::shared_leaf<T, SAFE_MATH> ti_cache;
1260
1262 graph::shared_vector<T, SAFE_MATH> b_cache;
1263
1265 graph::shared_leaf<T, SAFE_MATH> psi_cache;
1266
1267//------------------------------------------------------------------------------
1277//------------------------------------------------------------------------------
1278 graph::shared_leaf<T, SAFE_MATH>
1279 build_psi(graph::shared_leaf<T, SAFE_MATH> r,
1280 const T r_scale,
1281 const T r_offset,
1282 graph::shared_leaf<T, SAFE_MATH> z,
1283 const T z_scale,
1284 const T z_offset) {
1285 auto c00_temp = graph::piecewise_2D(c00, num_cols, r, r_scale, r_offset, z, z_scale, z_offset);
1286 auto c01_temp = graph::piecewise_2D(c01, num_cols, r, r_scale, r_offset, z, z_scale, z_offset);
1287 auto c02_temp = graph::piecewise_2D(c02, num_cols, r, r_scale, r_offset, z, z_scale, z_offset);
1288 auto c03_temp = graph::piecewise_2D(c03, num_cols, r, r_scale, r_offset, z, z_scale, z_offset);
1289
1290 auto c10_temp = graph::piecewise_2D(c10, num_cols, r, r_scale, r_offset, z, z_scale, z_offset);
1291 auto c11_temp = graph::piecewise_2D(c11, num_cols, r, r_scale, r_offset, z, z_scale, z_offset);
1292 auto c12_temp = graph::piecewise_2D(c12, num_cols, r, r_scale, r_offset, z, z_scale, z_offset);
1293 auto c13_temp = graph::piecewise_2D(c13, num_cols, r, r_scale, r_offset, z, z_scale, z_offset);
1294
1295 auto c20_temp = graph::piecewise_2D(c20, num_cols, r, r_scale, r_offset, z, z_scale, z_offset);
1296 auto c21_temp = graph::piecewise_2D(c21, num_cols, r, r_scale, r_offset, z, z_scale, z_offset);
1297 auto c22_temp = graph::piecewise_2D(c22, num_cols, r, r_scale, r_offset, z, z_scale, z_offset);
1298 auto c23_temp = graph::piecewise_2D(c23, num_cols, r, r_scale, r_offset, z, z_scale, z_offset);
1299
1300 auto c30_temp = graph::piecewise_2D(c30, num_cols, r, r_scale, r_offset, z, z_scale, z_offset);
1301 auto c31_temp = graph::piecewise_2D(c31, num_cols, r, r_scale, r_offset, z, z_scale, z_offset);
1302 auto c32_temp = graph::piecewise_2D(c32, num_cols, r, r_scale, r_offset, z, z_scale, z_offset);
1303 auto c33_temp = graph::piecewise_2D(c33, num_cols, r, r_scale, r_offset, z, z_scale, z_offset);
1304
1305 auto r_norm = (r - r_offset)/r_scale;
1306
1307 auto c0 = build_1D_spline({c00_temp, c01_temp, c02_temp, c03_temp}, z, z_scale, z_offset);
1308 auto c1 = build_1D_spline({c10_temp, c11_temp, c12_temp, c13_temp}, z, z_scale, z_offset);
1309 auto c2 = build_1D_spline({c20_temp, c21_temp, c22_temp, c23_temp}, z, z_scale, z_offset);
1310 auto c3 = build_1D_spline({c30_temp, c31_temp, c32_temp, c33_temp}, z, z_scale, z_offset);
1311
1312 return ((c3*r_norm + c2)*r_norm + c1)*r_norm + c0;
1313 }
1314
1315//------------------------------------------------------------------------------
1323//------------------------------------------------------------------------------
1324 void set_cache(graph::shared_leaf<T, SAFE_MATH> x,
1325 graph::shared_leaf<T, SAFE_MATH> y,
1326 graph::shared_leaf<T, SAFE_MATH> z) {
1327 if (!x->is_match(x_cache) ||
1328 !y->is_match(y_cache) ||
1329 !z->is_match(z_cache)) {
1330 x_cache = x;
1331 y_cache = y;
1332 z_cache = z;
1333
1334 auto r = graph::sqrt(x*x + y*y);
1335
1336 psi_cache = build_psi(r, dr, rmin, z, dz, zmin);
1337
1338 auto n0_temp = graph::piecewise_1D(ne_c0, psi_cache, dpsi, psimin);
1339 auto n1_temp = graph::piecewise_1D(ne_c1, psi_cache, dpsi, psimin);
1340 auto n2_temp = graph::piecewise_1D(ne_c2, psi_cache, dpsi, psimin);
1341 auto n3_temp = graph::piecewise_1D(ne_c3, psi_cache, dpsi, psimin);
1342
1343 ne_cache = ne_scale*build_1D_spline({n0_temp, n1_temp, n2_temp, n3_temp}, psi_cache, dpsi, psimin);
1344
1345 auto t0_temp = graph::piecewise_1D(te_c0, psi_cache, dpsi, psimin);
1346 auto t1_temp = graph::piecewise_1D(te_c1, psi_cache, dpsi, psimin);
1347 auto t2_temp = graph::piecewise_1D(te_c2, psi_cache, dpsi, psimin);
1348 auto t3_temp = graph::piecewise_1D(te_c3, psi_cache, dpsi, psimin);
1349
1350 te_cache = te_scale*build_1D_spline({t0_temp, t1_temp, t2_temp, t3_temp}, psi_cache, dpsi, psimin);
1351
1352 auto p0_temp = graph::piecewise_1D(pres_c0, psi_cache, dpsi, psimin);
1353 auto p1_temp = graph::piecewise_1D(pres_c1, psi_cache, dpsi, psimin);
1354 auto p2_temp = graph::piecewise_1D(pres_c2, psi_cache, dpsi, psimin);
1355 auto p3_temp = graph::piecewise_1D(pres_c3, psi_cache, dpsi, psimin);
1356
1357 auto pressure = pres_scale*build_1D_spline({p0_temp, p1_temp, p2_temp, p3_temp}, psi_cache, dpsi, psimin);
1358
1359 auto q = graph::constant<T, SAFE_MATH> (static_cast<T> (1.60218E-19));
1360
1361 ni_cache = te_cache;
1362 ti_cache = (pressure - ne_cache*te_cache*q)/(ni_cache*q);
1363
1364 auto phi = graph::atan(x, y);
1365
1366 auto br = psi_cache->df(z)/r;
1367
1368 auto b0_temp = graph::piecewise_1D(fpol_c0, psi_cache, dpsi, psimin);
1369 auto b1_temp = graph::piecewise_1D(fpol_c1, psi_cache, dpsi, psimin);
1370 auto b2_temp = graph::piecewise_1D(fpol_c2, psi_cache, dpsi, psimin);
1371 auto b3_temp = graph::piecewise_1D(fpol_c3, psi_cache, dpsi, psimin);
1372
1373 auto bp = build_1D_spline({b0_temp, b1_temp, b2_temp, b3_temp}, psi_cache, dpsi, psimin)/r;
1374
1375 auto bz = -psi_cache->df(r)/r;
1376
1377 auto cos = graph::cos(phi);
1378 auto sin = graph::sin(phi);
1379
1380 b_cache = graph::vector(br*cos - bp*sin,
1381 br*sin + bp*cos,
1382 bz);
1383 }
1384 }
1385
1386 public:
1387//------------------------------------------------------------------------------
1432//------------------------------------------------------------------------------
1433 efit(const T psimin,
1434 const T dpsi,
1435 const backend::buffer<T> te_c0,
1436 const backend::buffer<T> te_c1,
1437 const backend::buffer<T> te_c2,
1438 const backend::buffer<T> te_c3,
1439 graph::shared_leaf<T, SAFE_MATH> te_scale,
1440 const backend::buffer<T> ne_c0,
1441 const backend::buffer<T> ne_c1,
1442 const backend::buffer<T> ne_c2,
1443 const backend::buffer<T> ne_c3,
1444 graph::shared_leaf<T, SAFE_MATH> ne_scale,
1445 const backend::buffer<T> pres_c0,
1446 const backend::buffer<T> pres_c1,
1447 const backend::buffer<T> pres_c2,
1448 const backend::buffer<T> pres_c3,
1449 graph::shared_leaf<T, SAFE_MATH> pres_scale,
1450 const T rmin,
1451 const T dr,
1452 const T zmin,
1453 const T dz,
1454 const backend::buffer<T> fpol_c0,
1455 const backend::buffer<T> fpol_c1,
1456 const backend::buffer<T> fpol_c2,
1457 const backend::buffer<T> fpol_c3,
1458 const size_t num_cols,
1459 const backend::buffer<T> c00,
1460 const backend::buffer<T> c01,
1461 const backend::buffer<T> c02,
1462 const backend::buffer<T> c03,
1463 const backend::buffer<T> c10,
1464 const backend::buffer<T> c11,
1465 const backend::buffer<T> c12,
1466 const backend::buffer<T> c13,
1467 const backend::buffer<T> c20,
1468 const backend::buffer<T> c21,
1469 const backend::buffer<T> c22,
1470 const backend::buffer<T> c23,
1471 const backend::buffer<T> c30,
1472 const backend::buffer<T> c31,
1473 const backend::buffer<T> c32,
1474 const backend::buffer<T> c33) :
1475 generic<T, SAFE_MATH> ({3.34449469E-27} ,{1}),
1476 psimin(psimin), dpsi(dpsi), num_cols(num_cols),
1477 te_c0(te_c0), te_c1(te_c1), te_c2(te_c2), te_c3(te_c3), te_scale(te_scale),
1478 ne_c0(te_c0), ne_c1(te_c1), ne_c2(ne_c2), ne_c3(ne_c3), ne_scale(ne_scale),
1479 pres_c0(pres_c0), pres_c1(pres_c1), pres_c2(pres_c2), pres_c3(pres_c3),
1480 pres_scale(pres_scale), rmin(rmin), dr(dr), zmin(zmin), dz(dz),
1481 fpol_c0(fpol_c0), fpol_c1(fpol_c1), fpol_c2(fpol_c2), fpol_c3(fpol_c3),
1482 c00(c00), c01(c01), c02(c02), c03(c03),
1483 c10(c10), c11(c11), c12(c12), c13(c13),
1484 c20(c20), c21(c21), c22(c22), c23(c23),
1485 c30(c30), c31(c31), c32(c32), c33(c33) {
1486 auto zero = graph::zero<T, SAFE_MATH> ();
1487 x_cache = zero;
1488 y_cache = zero;
1489 z_cache = zero;
1490 }
1491
1492//------------------------------------------------------------------------------
1499//------------------------------------------------------------------------------
1500 virtual graph::shared_leaf<T, SAFE_MATH>
1507
1508//------------------------------------------------------------------------------
1516//------------------------------------------------------------------------------
1517 virtual graph::shared_leaf<T, SAFE_MATH>
1518 get_ion_density(const size_t index,
1519 graph::shared_leaf<T, SAFE_MATH> x,
1520 graph::shared_leaf<T, SAFE_MATH> y,
1521 graph::shared_leaf<T, SAFE_MATH> z) {
1522 set_cache(x, y, z);
1523 return ni_cache;
1524 }
1525
1526//------------------------------------------------------------------------------
1533//------------------------------------------------------------------------------
1534 virtual graph::shared_leaf<T, SAFE_MATH>
1541
1542//------------------------------------------------------------------------------
1550//------------------------------------------------------------------------------
1551 virtual graph::shared_leaf<T, SAFE_MATH>
1552 get_ion_temperature(const size_t index,
1553 graph::shared_leaf<T, SAFE_MATH> x,
1554 graph::shared_leaf<T, SAFE_MATH> y,
1555 graph::shared_leaf<T, SAFE_MATH> z) {
1556 set_cache(x, y, z);
1557 return ti_cache;
1558 }
1559
1560//------------------------------------------------------------------------------
1567//------------------------------------------------------------------------------
1568 virtual graph::shared_vector<T, SAFE_MATH>
1575
1576//------------------------------------------------------------------------------
1583//------------------------------------------------------------------------------
1584 virtual graph::shared_leaf<T, SAFE_MATH>
1585 get_characteristic_field(const size_t device_number=0) final {
1586 auto x_axis = graph::variable<T, SAFE_MATH> (1, "x");
1587 auto y_axis = graph::variable<T, SAFE_MATH> (1, "y");
1588 auto z_axis = graph::variable<T, SAFE_MATH> (1, "z");
1589 x_axis->set(static_cast<T> (1.7));
1590 y_axis->set(static_cast<T> (0.0));
1591 z_axis->set(static_cast<T> (0.0));
1592 auto b_vec = get_magnetic_field(x_axis, y_axis, z_axis);
1593 auto b_mod = b_vec->length();
1594
1595 graph::input_nodes<T, SAFE_MATH> inputs {
1596 graph::variable_cast(x_axis),
1597 graph::variable_cast(y_axis),
1598 graph::variable_cast(z_axis)
1599 };
1600
1601 workflow::manager<T, SAFE_MATH> work(device_number);
1602 solver::newton(work, {
1603 x_axis, z_axis
1604 }, inputs, (psi_cache - psimin)/dpsi, graph::shared_random_state<T, SAFE_MATH> (), static_cast<T> (1.0E-30), 1000, static_cast<T> (0.1));
1605 work.add_item(inputs, {b_mod}, {},
1606 graph::shared_random_state<T, SAFE_MATH> (),
1607 "bmod_at_axis", inputs.back()->size());
1608 work.compile();
1609 work.run();
1610
1611 T result;
1612 work.copy_to_host(b_mod, &result);
1613
1614 return graph::constant<T, SAFE_MATH> (result);
1615 }
1616 };
1617
1618//------------------------------------------------------------------------------
1626//------------------------------------------------------------------------------
1627 template<jit::float_scalar T, bool SAFE_MATH=false>
1628 shared<T, SAFE_MATH> make_efit(const std::string &spline_file) {
1629 int ncid;
1630 output::sync.lock();
1631 nc_open(spline_file.c_str(), NC_NOWRITE, &ncid);
1632
1633// Load scalar quantities.
1634 int varid;
1635
1636 double rmin_value;
1637 nc_inq_varid(ncid, "rmin", &varid);
1638 nc_get_var(ncid, varid, &rmin_value);
1639
1640 double dr_value;
1641 nc_inq_varid(ncid, "dr", &varid);
1642 nc_get_var(ncid, varid, &dr_value);
1643
1644 double zmin_value;
1645 nc_inq_varid(ncid, "zmin", &varid);
1646 nc_get_var(ncid, varid, &zmin_value);
1647
1648 double dz_value;
1649 nc_inq_varid(ncid, "dz", &varid);
1650 nc_get_var(ncid, varid, &dz_value);
1651
1652 double psimin_value;
1653 nc_inq_varid(ncid, "psimin", &varid);
1654 nc_get_var(ncid, varid, &psimin_value);
1655
1656 double dpsi_value;
1657 nc_inq_varid(ncid, "dpsi", &varid);
1658 nc_get_var(ncid, varid, &dpsi_value);
1659
1660 double pres_scale_value;
1661 nc_inq_varid(ncid, "pres_scale", &varid);
1662 nc_get_var(ncid, varid, &pres_scale_value);
1663
1664 double ne_scale_value;
1665 nc_inq_varid(ncid, "ne_scale", &varid);
1666 nc_get_var(ncid, varid, &ne_scale_value);
1667
1668 double te_scale_value;
1669 nc_inq_varid(ncid, "te_scale", &varid);
1670 nc_get_var(ncid, varid, &te_scale_value);
1671
1672// Load 1D quantities.
1673 int dimid;
1674
1675 size_t numr;
1676 nc_inq_dimid(ncid, "numr", &dimid);
1677 nc_inq_dimlen(ncid, dimid, &numr);
1678
1679 size_t numpsi;
1680 nc_inq_dimid(ncid, "numpsi", &dimid);
1681 nc_inq_dimlen(ncid, dimid, &numpsi);
1682
1683 std::vector<double> fpol_c0_buffer(numpsi);
1684 std::vector<double> fpol_c1_buffer(numpsi);
1685 std::vector<double> fpol_c2_buffer(numpsi);
1686 std::vector<double> fpol_c3_buffer(numpsi);
1687
1688 nc_inq_varid(ncid, "fpol_c0", &varid);
1689 nc_get_var(ncid, varid, fpol_c0_buffer.data());
1690 nc_inq_varid(ncid, "fpol_c1", &varid);
1691 nc_get_var(ncid, varid, fpol_c1_buffer.data());
1692 nc_inq_varid(ncid, "fpol_c2", &varid);
1693 nc_get_var(ncid, varid, fpol_c2_buffer.data());
1694 nc_inq_varid(ncid, "fpol_c3", &varid);
1695 nc_get_var(ncid, varid, fpol_c3_buffer.data());
1696
1697// Load psi grids.
1698 size_t numz;
1699 nc_inq_dimid(ncid, "numz", &dimid);
1700 nc_inq_dimlen(ncid, dimid, &numz);
1701
1702 std::vector<double> psi_c00_buffer(numz*numr);
1703 std::vector<double> psi_c01_buffer(numz*numr);
1704 std::vector<double> psi_c02_buffer(numz*numr);
1705 std::vector<double> psi_c03_buffer(numz*numr);
1706 std::vector<double> psi_c10_buffer(numz*numr);
1707 std::vector<double> psi_c11_buffer(numz*numr);
1708 std::vector<double> psi_c12_buffer(numz*numr);
1709 std::vector<double> psi_c13_buffer(numz*numr);
1710 std::vector<double> psi_c20_buffer(numz*numr);
1711 std::vector<double> psi_c21_buffer(numz*numr);
1712 std::vector<double> psi_c22_buffer(numz*numr);
1713 std::vector<double> psi_c23_buffer(numz*numr);
1714 std::vector<double> psi_c30_buffer(numz*numr);
1715 std::vector<double> psi_c31_buffer(numz*numr);
1716 std::vector<double> psi_c32_buffer(numz*numr);
1717 std::vector<double> psi_c33_buffer(numz*numr);
1718
1719 nc_inq_varid(ncid, "psi_c00", &varid);
1720 nc_get_var(ncid, varid, psi_c00_buffer.data());
1721 nc_inq_varid(ncid, "psi_c01", &varid);
1722 nc_get_var(ncid, varid, psi_c01_buffer.data());
1723 nc_inq_varid(ncid, "psi_c02", &varid);
1724 nc_get_var(ncid, varid, psi_c02_buffer.data());
1725 nc_inq_varid(ncid, "psi_c03", &varid);
1726 nc_get_var(ncid, varid, psi_c03_buffer.data());
1727 nc_inq_varid(ncid, "psi_c10", &varid);
1728 nc_get_var(ncid, varid, psi_c10_buffer.data());
1729 nc_inq_varid(ncid, "psi_c11", &varid);
1730 nc_get_var(ncid, varid, psi_c11_buffer.data());
1731 nc_inq_varid(ncid, "psi_c12", &varid);
1732 nc_get_var(ncid, varid, psi_c12_buffer.data());
1733 nc_inq_varid(ncid, "psi_c13", &varid);
1734 nc_get_var(ncid, varid, psi_c13_buffer.data());
1735 nc_inq_varid(ncid, "psi_c20", &varid);
1736 nc_get_var(ncid, varid, psi_c20_buffer.data());
1737 nc_inq_varid(ncid, "psi_c21", &varid);
1738 nc_get_var(ncid, varid, psi_c21_buffer.data());
1739 nc_inq_varid(ncid, "psi_c22", &varid);
1740 nc_get_var(ncid, varid, psi_c22_buffer.data());
1741 nc_inq_varid(ncid, "psi_c23", &varid);
1742 nc_get_var(ncid, varid, psi_c23_buffer.data());
1743 nc_inq_varid(ncid, "psi_c30", &varid);
1744 nc_get_var(ncid, varid, psi_c30_buffer.data());
1745 nc_inq_varid(ncid, "psi_c31", &varid);
1746 nc_get_var(ncid, varid, psi_c31_buffer.data());
1747 nc_inq_varid(ncid, "psi_c32", &varid);
1748 nc_get_var(ncid, varid, psi_c32_buffer.data());
1749 nc_inq_varid(ncid, "psi_c33", &varid);
1750 nc_get_var(ncid, varid, psi_c33_buffer.data());
1751
1752 std::vector<double> pressure_c0_buffer(numpsi);
1753 std::vector<double> pressure_c1_buffer(numpsi);
1754 std::vector<double> pressure_c2_buffer(numpsi);
1755 std::vector<double> pressure_c3_buffer(numpsi);
1756
1757 nc_inq_varid(ncid, "pressure_c0", &varid);
1758 nc_get_var(ncid, varid, pressure_c0_buffer.data());
1759 nc_inq_varid(ncid, "pressure_c1", &varid);
1760 nc_get_var(ncid, varid, pressure_c1_buffer.data());
1761 nc_inq_varid(ncid, "pressure_c2", &varid);
1762 nc_get_var(ncid, varid, pressure_c2_buffer.data());
1763 nc_inq_varid(ncid, "pressure_c3", &varid);
1764 nc_get_var(ncid, varid, pressure_c3_buffer.data());
1765
1766 std::vector<double> te_c0_buffer(numpsi);
1767 std::vector<double> te_c1_buffer(numpsi);
1768 std::vector<double> te_c2_buffer(numpsi);
1769 std::vector<double> te_c3_buffer(numpsi);
1770
1771 nc_inq_varid(ncid, "te_c0", &varid);
1772 nc_get_var(ncid, varid, te_c0_buffer.data());
1773 nc_inq_varid(ncid, "te_c1", &varid);
1774 nc_get_var(ncid, varid, te_c1_buffer.data());
1775 nc_inq_varid(ncid, "te_c2", &varid);
1776 nc_get_var(ncid, varid, te_c2_buffer.data());
1777 nc_inq_varid(ncid, "te_c3", &varid);
1778 nc_get_var(ncid, varid, te_c3_buffer.data());
1779
1780 std::vector<double> ne_c0_buffer(numpsi);
1781 std::vector<double> ne_c1_buffer(numpsi);
1782 std::vector<double> ne_c2_buffer(numpsi);
1783 std::vector<double> ne_c3_buffer(numpsi);
1784
1785 nc_inq_varid(ncid, "ne_c0", &varid);
1786 nc_get_var(ncid, varid, ne_c0_buffer.data());
1787 nc_inq_varid(ncid, "ne_c1", &varid);
1788 nc_get_var(ncid, varid, ne_c1_buffer.data());
1789 nc_inq_varid(ncid, "ne_c2", &varid);
1790 nc_get_var(ncid, varid, ne_c2_buffer.data());
1791 nc_inq_varid(ncid, "ne_c3", &varid);
1792 nc_get_var(ncid, varid, ne_c3_buffer.data());
1793
1794 nc_close(ncid);
1795 output::sync.unlock();
1796
1797 auto rmin = static_cast<T> (rmin_value);
1798 auto dr = static_cast<T> (dr_value);
1799 auto zmin = static_cast<T> (zmin_value);
1800 auto dz = static_cast<T> (dz_value);
1801 auto psimin = static_cast<T> (psimin_value);
1802 auto dpsi = static_cast<T> (dpsi_value);
1803 auto pres_scale = graph::constant<T, SAFE_MATH> (static_cast<T> (pres_scale_value));
1804 auto ne_scale = graph::constant<T, SAFE_MATH> (static_cast<T> (ne_scale_value));
1805 auto te_scale = graph::constant<T, SAFE_MATH> (static_cast<T> (te_scale_value));
1806
1807 const auto fpol_c0 = backend::buffer(std::vector<T> (fpol_c0_buffer.begin(), fpol_c0_buffer.end()));
1808 const auto fpol_c1 = backend::buffer(std::vector<T> (fpol_c1_buffer.begin(), fpol_c1_buffer.end()));
1809 const auto fpol_c2 = backend::buffer(std::vector<T> (fpol_c2_buffer.begin(), fpol_c2_buffer.end()));
1810 const auto fpol_c3 = backend::buffer(std::vector<T> (fpol_c3_buffer.begin(), fpol_c3_buffer.end()));
1811
1812 const auto c00 = backend::buffer(std::vector<T> (psi_c00_buffer.begin(), psi_c00_buffer.end()));
1813 const auto c01 = backend::buffer(std::vector<T> (psi_c01_buffer.begin(), psi_c01_buffer.end()));
1814 const auto c02 = backend::buffer(std::vector<T> (psi_c02_buffer.begin(), psi_c02_buffer.end()));
1815 const auto c03 = backend::buffer(std::vector<T> (psi_c03_buffer.begin(), psi_c03_buffer.end()));
1816 const auto c10 = backend::buffer(std::vector<T> (psi_c10_buffer.begin(), psi_c10_buffer.end()));
1817 const auto c11 = backend::buffer(std::vector<T> (psi_c11_buffer.begin(), psi_c11_buffer.end()));
1818 const auto c12 = backend::buffer(std::vector<T> (psi_c12_buffer.begin(), psi_c12_buffer.end()));
1819 const auto c13 = backend::buffer(std::vector<T> (psi_c13_buffer.begin(), psi_c13_buffer.end()));
1820 const auto c20 = backend::buffer(std::vector<T> (psi_c20_buffer.begin(), psi_c20_buffer.end()));
1821 const auto c21 = backend::buffer(std::vector<T> (psi_c21_buffer.begin(), psi_c21_buffer.end()));
1822 const auto c22 = backend::buffer(std::vector<T> (psi_c22_buffer.begin(), psi_c22_buffer.end()));
1823 const auto c23 = backend::buffer(std::vector<T> (psi_c23_buffer.begin(), psi_c23_buffer.end()));
1824 const auto c30 = backend::buffer(std::vector<T> (psi_c30_buffer.begin(), psi_c30_buffer.end()));
1825 const auto c31 = backend::buffer(std::vector<T> (psi_c31_buffer.begin(), psi_c31_buffer.end()));
1826 const auto c32 = backend::buffer(std::vector<T> (psi_c32_buffer.begin(), psi_c32_buffer.end()));
1827 const auto c33 = backend::buffer(std::vector<T> (psi_c33_buffer.begin(), psi_c33_buffer.end()));
1828
1829 const auto pres_c0 = backend::buffer(std::vector<T> (pressure_c0_buffer.begin(), pressure_c0_buffer.end()));
1830 const auto pres_c1 = backend::buffer(std::vector<T> (pressure_c1_buffer.begin(), pressure_c1_buffer.end()));
1831 const auto pres_c2 = backend::buffer(std::vector<T> (pressure_c2_buffer.begin(), pressure_c2_buffer.end()));
1832 const auto pres_c3 = backend::buffer(std::vector<T> (pressure_c3_buffer.begin(), pressure_c3_buffer.end()));
1833
1834 const auto te_c0 = backend::buffer(std::vector<T> (te_c0_buffer.begin(), te_c0_buffer.end()));
1835 const auto te_c1 = backend::buffer(std::vector<T> (te_c1_buffer.begin(), te_c1_buffer.end()));
1836 const auto te_c2 = backend::buffer(std::vector<T> (te_c2_buffer.begin(), te_c2_buffer.end()));
1837 const auto te_c3 = backend::buffer(std::vector<T> (te_c3_buffer.begin(), te_c3_buffer.end()));
1838
1839 const auto ne_c0 = backend::buffer(std::vector<T> (ne_c0_buffer.begin(), ne_c0_buffer.end()));
1840 const auto ne_c1 = backend::buffer(std::vector<T> (ne_c1_buffer.begin(), ne_c1_buffer.end()));
1841 const auto ne_c2 = backend::buffer(std::vector<T> (ne_c2_buffer.begin(), ne_c2_buffer.end()));
1842 const auto ne_c3 = backend::buffer(std::vector<T> (ne_c3_buffer.begin(), ne_c3_buffer.end()));
1843
1844 return std::make_shared<efit<T, SAFE_MATH>> (psimin, dpsi,
1845 te_c0, te_c1, te_c2, te_c3, te_scale,
1846 ne_c0, ne_c1, ne_c2, ne_c3, ne_scale,
1847 pres_c0, pres_c1, pres_c2, pres_c3, pres_scale,
1848 rmin, dr, zmin, dz,
1849 fpol_c0, fpol_c1, fpol_c2, fpol_c3, numz,
1850 c00, c01, c02, c03,
1851 c10, c11, c12, c13,
1852 c20, c21, c22, c23,
1853 c30, c31, c32, c33);
1854 }
1855
1856//******************************************************************************
1857// 3D VMEC equilibrium.
1858//******************************************************************************
1859//------------------------------------------------------------------------------
1866//------------------------------------------------------------------------------
1867 template<jit::float_scalar T, bool SAFE_MATH=false>
1868 class vmec final : public generic<T, SAFE_MATH> {
1869 private:
1871 const T sminh;
1873 const T sminf;
1875 const T ds;
1877 graph::shared_leaf<T, SAFE_MATH> signj;
1878
1879// Poloidal flux coefficients.
1881 const backend::buffer<T> chi_c0;
1883 const backend::buffer<T> chi_c1;
1885 const backend::buffer<T> chi_c2;
1887 const backend::buffer<T> chi_c3;
1888
1889// Toroidal flux coefficients.
1890 graph::shared_leaf<T, SAFE_MATH> dphi;
1891
1892// Radial coefficients.
1894 const std::vector<backend::buffer<T>> rmnc_c0;
1896 const std::vector<backend::buffer<T>> rmnc_c1;
1898 const std::vector<backend::buffer<T>> rmnc_c2;
1900 const std::vector<backend::buffer<T>> rmnc_c3;
1901
1902// Vertical coefficients.
1904 const std::vector<backend::buffer<T>> zmns_c0;
1906 const std::vector<backend::buffer<T>> zmns_c1;
1908 const std::vector<backend::buffer<T>> zmns_c2;
1910 const std::vector<backend::buffer<T>> zmns_c3;
1911
1912// Lambda coefficients.
1914 const std::vector<backend::buffer<T>> lmns_c0;
1916 const std::vector<backend::buffer<T>> lmns_c1;
1918 const std::vector<backend::buffer<T>> lmns_c2;
1920 const std::vector<backend::buffer<T>> lmns_c3;
1921
1923 const backend::buffer<T> xm;
1925 const backend::buffer<T> xn;
1926
1927// Cached values.
1929 graph::shared_leaf<T, SAFE_MATH> s_cache;
1931 graph::shared_leaf<T, SAFE_MATH> u_cache;
1933 graph::shared_leaf<T, SAFE_MATH> v_cache;
1935 graph::shared_leaf<T, SAFE_MATH> x_cache;
1937 graph::shared_leaf<T, SAFE_MATH> y_cache;
1939 graph::shared_leaf<T, SAFE_MATH> z_cache;
1940
1942 graph::shared_vector<T, SAFE_MATH> esups_cache;
1944 graph::shared_vector<T, SAFE_MATH> esupu_cache;
1946 graph::shared_vector<T, SAFE_MATH> esupv_cache;
1947
1949 graph::shared_vector<T, SAFE_MATH> bvec_cache;
1950
1951//------------------------------------------------------------------------------
1957//------------------------------------------------------------------------------
1958 graph::shared_vector<T, SAFE_MATH>
1959 get_esubs(graph::shared_leaf<T, SAFE_MATH> r,
1960 graph::shared_leaf<T, SAFE_MATH> z) {
1961 auto cosv = graph::cos(v_cache);
1962 auto sinv = graph::sin(v_cache);
1963 auto one = graph::one<T, SAFE_MATH> ();
1964 auto zero = graph::zero<T, SAFE_MATH> ();
1965
1966 auto m = graph::matrix(graph::vector(cosv, -sinv, zero),
1967 graph::vector(sinv, cosv, zero),
1968 graph::vector(zero, zero, one ));
1969 return m->dot(graph::vector(r->df(s_cache),
1970 zero,
1971 z->df(s_cache)));
1972 }
1973
1974//------------------------------------------------------------------------------
1980//------------------------------------------------------------------------------
1981 graph::shared_vector<T, SAFE_MATH>
1982 get_esubu(graph::shared_leaf<T, SAFE_MATH> r,
1983 graph::shared_leaf<T, SAFE_MATH> z) {
1984 auto cosv = graph::cos(v_cache);
1985 auto sinv = graph::sin(v_cache);
1986 auto one = graph::one<T, SAFE_MATH> ();
1987 auto zero = graph::zero<T, SAFE_MATH> ();
1988
1989 auto m = graph::matrix(graph::vector(cosv, -sinv, zero),
1990 graph::vector(sinv, cosv, zero),
1991 graph::vector(zero, zero, one ));
1992 return m->dot(graph::vector(r->df(u_cache),
1993 zero,
1994 z->df(u_cache)));
1995 }
1996
1997//------------------------------------------------------------------------------
2003//------------------------------------------------------------------------------
2004 graph::shared_vector<T, SAFE_MATH>
2005 get_esubv(graph::shared_leaf<T, SAFE_MATH> r,
2006 graph::shared_leaf<T, SAFE_MATH> z) {
2007 auto cosv = graph::cos(v_cache);
2008 auto sinv = graph::sin(v_cache);
2009 auto one = graph::one<T, SAFE_MATH> ();
2010 auto zero = graph::zero<T, SAFE_MATH> ();
2011
2012 auto m = graph::matrix(graph::vector(cosv, -sinv, zero),
2013 graph::vector(sinv, cosv, zero),
2014 graph::vector(zero, zero, one ));
2015 return m->dot(graph::vector(r->df(v_cache),
2016 r,
2017 z->df(v_cache)));
2018 }
2019
2020//------------------------------------------------------------------------------
2029//------------------------------------------------------------------------------
2030 graph::shared_leaf<T, SAFE_MATH>
2031 get_jacobian(graph::shared_vector<T, SAFE_MATH> esub_s,
2032 graph::shared_vector<T, SAFE_MATH> esub_u,
2033 graph::shared_vector<T, SAFE_MATH> esub_v) {
2034 return esub_s->dot(esub_u->cross(esub_v));
2035 }
2036
2037//------------------------------------------------------------------------------
2042//------------------------------------------------------------------------------
2043 graph::shared_leaf<T, SAFE_MATH>
2044 get_chi(graph::shared_leaf<T, SAFE_MATH> s) {
2045 auto c0_temp = graph::piecewise_1D(chi_c0, s, ds, sminf);
2046 auto c1_temp = graph::piecewise_1D(chi_c1, s, ds, sminf);
2047 auto c2_temp = graph::piecewise_1D(chi_c2, s, ds, sminf);
2048 auto c3_temp = graph::piecewise_1D(chi_c3, s, ds, sminf);
2049
2050 return build_1D_spline({c0_temp, c1_temp, c2_temp, c3_temp}, s, ds, sminf);
2051 }
2052
2053//------------------------------------------------------------------------------
2058//------------------------------------------------------------------------------
2059 graph::shared_leaf<T, SAFE_MATH>
2060 get_phi(graph::shared_leaf<T, SAFE_MATH> s) {
2061 return signj*dphi*s;
2062 }
2063
2064//------------------------------------------------------------------------------
2072//------------------------------------------------------------------------------
2073 void set_cache(graph::shared_leaf<T, SAFE_MATH> s,
2074 graph::shared_leaf<T, SAFE_MATH> u,
2075 graph::shared_leaf<T, SAFE_MATH> v) {
2076 if (!s->is_match(s_cache) ||
2077 !u->is_match(u_cache) ||
2078 !v->is_match(v_cache)) {
2079 s_cache = s;
2080 u_cache = u;
2081 v_cache = v;
2082
2083 auto s_norm_f = (s - sminf)/ds;
2084
2085 auto zero = graph::zero<T, SAFE_MATH> ();
2086 auto r = zero;
2087 auto z = zero;
2088 auto l = zero;
2089
2090 for (size_t i = 0, ie = xm.size(); i < ie; i++) {
2091 auto rmnc_c0_temp = graph::piecewise_1D(rmnc_c0[i], s, ds, sminf);
2092 auto rmnc_c1_temp = graph::piecewise_1D(rmnc_c1[i], s, ds, sminf);
2093 auto rmnc_c2_temp = graph::piecewise_1D(rmnc_c2[i], s, ds, sminf);
2094 auto rmnc_c3_temp = graph::piecewise_1D(rmnc_c3[i], s, ds, sminf);
2095
2096 auto zmns_c0_temp = graph::piecewise_1D(zmns_c0[i], s, ds, sminf);
2097 auto zmns_c1_temp = graph::piecewise_1D(zmns_c1[i], s, ds, sminf);
2098 auto zmns_c2_temp = graph::piecewise_1D(zmns_c2[i], s, ds, sminf);
2099 auto zmns_c3_temp = graph::piecewise_1D(zmns_c3[i], s, ds, sminf);
2100
2101 auto lmns_c0_temp = graph::piecewise_1D(lmns_c0[i], s, ds, sminh);
2102 auto lmns_c1_temp = graph::piecewise_1D(lmns_c1[i], s, ds, sminh);
2103 auto lmns_c2_temp = graph::piecewise_1D(lmns_c2[i], s, ds, sminh);
2104 auto lmns_c3_temp = graph::piecewise_1D(lmns_c3[i], s, ds, sminh);
2105
2106 auto rmnc = build_1D_spline({rmnc_c0_temp, rmnc_c1_temp, rmnc_c2_temp, rmnc_c3_temp},
2107 s, ds, sminf);
2108 auto zmns = build_1D_spline({zmns_c0_temp, zmns_c1_temp, zmns_c2_temp, zmns_c3_temp},
2109 s, ds, sminf);
2110 auto lmns = build_1D_spline({lmns_c0_temp, lmns_c1_temp, lmns_c2_temp, lmns_c3_temp},
2111 s, ds, sminh);
2112
2113 auto m = graph::constant<T, SAFE_MATH> (xm[i]);
2114 auto n = graph::constant<T, SAFE_MATH> (xn[i]);
2115
2116 auto sinmn = graph::sin(m*u - n*v);
2117
2118 r = r + rmnc*graph::cos(m*u - n*v);
2119 z = z + zmns*sinmn;
2120 l = l + lmns*sinmn;
2121 }
2122
2123 x_cache = r*graph::cos(v);
2124 y_cache = r*graph::sin(v);
2125 z_cache = z;
2126
2127 auto esubs = get_esubs(r, z);
2128 auto esubu = get_esubu(r, z);
2129 auto esubv = get_esubv(r, z);
2130
2131 auto jacobian = get_jacobian(esubs, esubu, esubv);
2132
2133 esups_cache = esubu->cross(esubv)/jacobian;
2134 esupu_cache = esubv->cross(esubs)/jacobian;
2135 esupv_cache = esubs->cross(esubu)/jacobian;
2136
2137 auto phip = get_phi(s)->df(s);
2138 auto jbsupu = get_chi(s_norm_f)->df(s) - phip*l->df(v);
2139 auto jbsupv = phip*(1.0 + l->df(u));
2140 bvec_cache = (jbsupu*esubu + jbsupv*esubv)/jacobian;
2141 }
2142 }
2143
2144//------------------------------------------------------------------------------
2149//------------------------------------------------------------------------------
2150 graph::shared_leaf<T, SAFE_MATH>
2151 get_profile(graph::shared_leaf<T, SAFE_MATH> s) {
2152 return graph::pow((1.0 - graph::pow(graph::sqrt(s*s), 1.5)), 2.0);
2153 }
2154
2155 public:
2156//------------------------------------------------------------------------------
2182//------------------------------------------------------------------------------
2183 vmec(const T sminh,
2184 const T sminf,
2185 const T ds,
2186 graph::shared_leaf<T, SAFE_MATH> dphi,
2187 graph::shared_leaf<T, SAFE_MATH> signj,
2188 const backend::buffer<T> chi_c0,
2189 const backend::buffer<T> chi_c1,
2190 const backend::buffer<T> chi_c2,
2191 const backend::buffer<T> chi_c3,
2192 const std::vector<backend::buffer<T>> rmnc_c0,
2193 const std::vector<backend::buffer<T>> rmnc_c1,
2194 const std::vector<backend::buffer<T>> rmnc_c2,
2195 const std::vector<backend::buffer<T>> rmnc_c3,
2196 const std::vector<backend::buffer<T>> zmns_c0,
2197 const std::vector<backend::buffer<T>> zmns_c1,
2198 const std::vector<backend::buffer<T>> zmns_c2,
2199 const std::vector<backend::buffer<T>> zmns_c3,
2200 const std::vector<backend::buffer<T>> lmns_c0,
2201 const std::vector<backend::buffer<T>> lmns_c1,
2202 const std::vector<backend::buffer<T>> lmns_c2,
2203 const std::vector<backend::buffer<T>> lmns_c3,
2204 const backend::buffer<T> xm,
2205 const backend::buffer<T> xn) :
2206 generic<T, SAFE_MATH> ({3.34449469E-27} ,{1}),
2207 sminh(sminh), sminf(sminf), ds(ds), dphi(dphi), signj(signj),
2208 chi_c0(chi_c0), chi_c1(chi_c1), chi_c2(chi_c2), chi_c3(chi_c3),
2209 rmnc_c0(rmnc_c0), rmnc_c1(rmnc_c1), rmnc_c2(rmnc_c2), rmnc_c3(rmnc_c3),
2210 zmns_c0(zmns_c0), zmns_c1(zmns_c1), zmns_c2(zmns_c2), zmns_c3(zmns_c3),
2211 lmns_c0(lmns_c0), lmns_c1(lmns_c1), lmns_c2(lmns_c2), lmns_c3(lmns_c3),
2212 xm(xm), xn(xn) {
2213 auto zero = graph::zero<T, SAFE_MATH> ();
2214 s_cache = zero;
2215 u_cache = zero;
2216 v_cache = zero;
2217 }
2218
2219//------------------------------------------------------------------------------
2226//------------------------------------------------------------------------------
2227 virtual graph::shared_vector<T, SAFE_MATH>
2228 get_esup1(graph::shared_leaf<T, SAFE_MATH> s,
2229 graph::shared_leaf<T, SAFE_MATH> u,
2230 graph::shared_leaf<T, SAFE_MATH> v) {
2231 set_cache(s, u, v);
2232 return esups_cache;
2233 }
2234
2235//------------------------------------------------------------------------------
2242//------------------------------------------------------------------------------
2243 virtual graph::shared_vector<T, SAFE_MATH>
2244 get_esup2(graph::shared_leaf<T, SAFE_MATH> s,
2245 graph::shared_leaf<T, SAFE_MATH> u,
2246 graph::shared_leaf<T, SAFE_MATH> v) {
2247 set_cache(s, u, v);
2248 return esupu_cache;
2249 }
2250
2251//------------------------------------------------------------------------------
2258//------------------------------------------------------------------------------
2259 virtual graph::shared_vector<T, SAFE_MATH>
2260 get_esup3(graph::shared_leaf<T, SAFE_MATH> s,
2261 graph::shared_leaf<T, SAFE_MATH> u,
2262 graph::shared_leaf<T, SAFE_MATH> v) {
2263 set_cache(s, u, v);
2264 return esupv_cache;
2265 }
2266
2267//------------------------------------------------------------------------------
2274//------------------------------------------------------------------------------
2275 virtual graph::shared_leaf<T, SAFE_MATH>
2276 get_electron_density(graph::shared_leaf<T, SAFE_MATH> s,
2277 graph::shared_leaf<T, SAFE_MATH> u,
2278 graph::shared_leaf<T, SAFE_MATH> v) {
2279 auto ne_scale = graph::constant<T, SAFE_MATH> (static_cast<T> (1.0E19));
2280 return ne_scale*get_profile(s);
2281 }
2282
2283//------------------------------------------------------------------------------
2291//------------------------------------------------------------------------------
2292 virtual graph::shared_leaf<T, SAFE_MATH>
2299
2300//------------------------------------------------------------------------------
2307//------------------------------------------------------------------------------
2308 virtual graph::shared_leaf<T, SAFE_MATH>
2309 get_electron_temperature(graph::shared_leaf<T, SAFE_MATH> s,
2310 graph::shared_leaf<T, SAFE_MATH> u,
2311 graph::shared_leaf<T, SAFE_MATH> v) {
2312 auto te_scale = graph::constant<T, SAFE_MATH> (static_cast<T> (1000.0));
2313 return te_scale*get_profile(s);
2314 }
2315
2316//------------------------------------------------------------------------------
2324//------------------------------------------------------------------------------
2325 virtual graph::shared_leaf<T, SAFE_MATH>
2332
2333//------------------------------------------------------------------------------
2340//------------------------------------------------------------------------------
2341 virtual graph::shared_vector<T, SAFE_MATH>
2342 get_magnetic_field(graph::shared_leaf<T, SAFE_MATH> s,
2343 graph::shared_leaf<T, SAFE_MATH> u,
2344 graph::shared_leaf<T, SAFE_MATH> v) {
2345 set_cache(s, u, v);
2346 return bvec_cache;
2347 }
2348
2349//------------------------------------------------------------------------------
2356//------------------------------------------------------------------------------
2357 virtual graph::shared_leaf<T, SAFE_MATH>
2358 get_characteristic_field(const size_t device_number=0) final {
2359 auto s_axis = graph::zero<T, SAFE_MATH> ();
2360 auto u_axis = graph::zero<T, SAFE_MATH> ();
2361 auto v_axis = graph::zero<T, SAFE_MATH> ();
2362 auto b_vec = get_magnetic_field(s_axis, u_axis, v_axis);
2363 return b_vec->length();
2364 }
2365
2366//------------------------------------------------------------------------------
2373//------------------------------------------------------------------------------
2374 virtual graph::shared_leaf<T, SAFE_MATH>
2375 get_x(graph::shared_leaf<T, SAFE_MATH> s,
2376 graph::shared_leaf<T, SAFE_MATH> u,
2377 graph::shared_leaf<T, SAFE_MATH> v) {
2378 set_cache(s, u, v);
2379 return x_cache;
2380 }
2381
2382//------------------------------------------------------------------------------
2389//------------------------------------------------------------------------------
2390 virtual graph::shared_leaf<T, SAFE_MATH>
2391 get_y(graph::shared_leaf<T, SAFE_MATH> s,
2392 graph::shared_leaf<T, SAFE_MATH> u,
2393 graph::shared_leaf<T, SAFE_MATH> v) {
2394 set_cache(s, u, v);
2395 return y_cache;
2396 }
2397
2398//------------------------------------------------------------------------------
2405//------------------------------------------------------------------------------
2406 virtual graph::shared_leaf<T, SAFE_MATH>
2407 get_z(graph::shared_leaf<T, SAFE_MATH> s,
2408 graph::shared_leaf<T, SAFE_MATH> u,
2409 graph::shared_leaf<T, SAFE_MATH> v) {
2410 set_cache(s, u, v);
2411 return z_cache;
2412 }
2413 };
2414
2415//------------------------------------------------------------------------------
2423//------------------------------------------------------------------------------
2424 template<jit::float_scalar T, bool SAFE_MATH=false>
2425 shared<T, SAFE_MATH> make_vmec(const std::string &spline_file) {
2426 int ncid;
2427 output::sync.lock();
2428 nc_open(spline_file.c_str(), NC_NOWRITE, &ncid);
2429
2430// Load scalar quantities.
2431 int varid;
2432
2433 double sminf_value;
2434 nc_inq_varid(ncid, "sminf", &varid);
2435 nc_get_var(ncid, varid, &sminf_value);
2436
2437 double sminh_value;
2438 nc_inq_varid(ncid, "sminh", &varid);
2439 nc_get_var(ncid, varid, &sminh_value);
2440
2441 double ds_value;
2442 nc_inq_varid(ncid, "ds", &varid);
2443 nc_get_var(ncid, varid, &ds_value);
2444
2445 double dphi_value;
2446 nc_inq_varid(ncid, "dphi", &varid);
2447 nc_get_var(ncid, varid, &dphi_value);
2448
2449 double signj_value;
2450 nc_inq_varid(ncid, "signj", &varid);
2451 nc_get_var(ncid, varid, &signj_value);
2452
2453// Load 1D quantities.
2454 int dimid;
2455
2456 size_t numsf;
2457 nc_inq_dimid(ncid, "numsf", &dimid);
2458 nc_inq_dimlen(ncid, dimid, &numsf);
2459
2460 std::vector<double> chi_c0_buffer(numsf);
2461 std::vector<double> chi_c1_buffer(numsf);
2462 std::vector<double> chi_c2_buffer(numsf);
2463 std::vector<double> chi_c3_buffer(numsf);
2464
2465 nc_inq_varid(ncid, "chi_c0", &varid);
2466 nc_get_var(ncid, varid, chi_c0_buffer.data());
2467 nc_inq_varid(ncid, "chi_c1", &varid);
2468 nc_get_var(ncid, varid, chi_c1_buffer.data());
2469 nc_inq_varid(ncid, "chi_c2", &varid);
2470 nc_get_var(ncid, varid, chi_c2_buffer.data());
2471 nc_inq_varid(ncid, "chi_c3", &varid);
2472 nc_get_var(ncid, varid, chi_c3_buffer.data());
2473
2474// Load 2D quantities.
2475 size_t numsh;
2476 nc_inq_dimid(ncid, "numsh", &dimid);
2477 nc_inq_dimlen(ncid, dimid, &numsh);
2478
2479 size_t nummn;
2480 nc_inq_dimid(ncid, "nummn", &dimid);
2481 nc_inq_dimlen(ncid, dimid, &nummn);
2482
2483 std::vector<std::vector<double>> rmnc_c0_buffer(nummn, std::vector<double> (numsf));
2484 std::vector<std::vector<double>> rmnc_c1_buffer(nummn, std::vector<double> (numsf));
2485 std::vector<std::vector<double>> rmnc_c2_buffer(nummn, std::vector<double> (numsf));
2486 std::vector<std::vector<double>> rmnc_c3_buffer(nummn, std::vector<double> (numsf));
2487
2488 nc_inq_varid(ncid, "rmnc_c0", &varid);
2489 for (size_t i = 0; i < nummn; i++) {
2490 const array<size_t, 2> start = {i, 0};
2491 const array<size_t, 2> count = {1, numsf};
2492 nc_get_vara(ncid, varid, start.data(), count.data(),
2493 rmnc_c0_buffer[i].data());
2494 }
2495 nc_inq_varid(ncid, "rmnc_c1", &varid);
2496 for (size_t i = 0; i < nummn; i++) {
2497 const array<size_t, 2> start = {i, 0};
2498 const array<size_t, 2> count = {1, numsf};
2499 nc_get_vara(ncid, varid, start.data(), count.data(),
2500 rmnc_c1_buffer[i].data());
2501 }
2502 nc_inq_varid(ncid, "rmnc_c2", &varid);
2503 for (size_t i = 0; i < nummn; i++) {
2504 const array<size_t, 2> start = {i, 0};
2505 const array<size_t, 2> count = {1, numsf};
2506 nc_get_vara(ncid, varid, start.data(), count.data(),
2507 rmnc_c2_buffer[i].data());
2508 }
2509 nc_inq_varid(ncid, "rmnc_c3", &varid);
2510 for (size_t i = 0; i < nummn; i++) {
2511 const array<size_t, 2> start = {i, 0};
2512 const array<size_t, 2> count = {1, numsf};
2513 nc_get_vara(ncid, varid, start.data(), count.data(),
2514 rmnc_c3_buffer[i].data());
2515 }
2516
2517 std::vector<std::vector<double>> zmns_c0_buffer(nummn, std::vector<double> (numsf));
2518 std::vector<std::vector<double>> zmns_c1_buffer(nummn, std::vector<double> (numsf));
2519 std::vector<std::vector<double>> zmns_c2_buffer(nummn, std::vector<double> (numsf));
2520 std::vector<std::vector<double>> zmns_c3_buffer(nummn, std::vector<double> (numsf));
2521
2522 nc_inq_varid(ncid, "zmns_c0", &varid);
2523 for (size_t i = 0; i < nummn; i++) {
2524 const array<size_t, 2> start = {i, 0};
2525 const array<size_t, 2> count = {1, numsf};
2526 nc_get_vara(ncid, varid, start.data(), count.data(),
2527 zmns_c0_buffer[i].data());
2528 }
2529 nc_inq_varid(ncid, "zmns_c1", &varid);
2530 for (size_t i = 0; i < nummn; i++) {
2531 const array<size_t, 2> start = {i, 0};
2532 const array<size_t, 2> count = {1, numsf};
2533 nc_get_vara(ncid, varid, start.data(), count.data(),
2534 zmns_c1_buffer[i].data());
2535 }
2536 nc_inq_varid(ncid, "zmns_c2", &varid);
2537 for (size_t i = 0; i < nummn; i++) {
2538 const array<size_t, 2> start = {i, 0};
2539 const array<size_t, 2> count = {1, numsf};
2540 nc_get_vara(ncid, varid, start.data(), count.data(),
2541 zmns_c2_buffer[i].data());
2542 }
2543 nc_inq_varid(ncid, "zmns_c3", &varid);
2544 for (size_t i = 0; i < nummn; i++) {
2545 const array<size_t, 2> start = {i, 0};
2546 const array<size_t, 2> count = {1, numsf};
2547 nc_get_vara(ncid, varid, start.data(), count.data(),
2548 zmns_c3_buffer[i].data());
2549 }
2550
2551 std::vector<std::vector<double>> lmns_c0_buffer(nummn, std::vector<double> (numsh));
2552 std::vector<std::vector<double>> lmns_c1_buffer(nummn, std::vector<double> (numsh));
2553 std::vector<std::vector<double>> lmns_c2_buffer(nummn, std::vector<double> (numsh));
2554 std::vector<std::vector<double>> lmns_c3_buffer(nummn, std::vector<double> (numsh));
2555
2556 nc_inq_varid(ncid, "lmns_c0", &varid);
2557 for (size_t i = 0; i < nummn; i++) {
2558 const array<size_t, 2> start = {i, 0};
2559 const array<size_t, 2> count = {1, numsh};
2560 nc_get_vara(ncid, varid, start.data(), count.data(),
2561 lmns_c0_buffer[i].data());
2562 }
2563 nc_inq_varid(ncid, "lmns_c1", &varid);
2564 for (size_t i = 0; i < nummn; i++) {
2565 const array<size_t, 2> start = {i, 0};
2566 const array<size_t, 2> count = {1, numsh};
2567 nc_get_vara(ncid, varid, start.data(), count.data(),
2568 lmns_c1_buffer[i].data());
2569 }
2570 nc_inq_varid(ncid, "lmns_c2", &varid);
2571 for (size_t i = 0; i < nummn; i++) {
2572 const array<size_t, 2> start = {i, 0};
2573 const array<size_t, 2> count = {1, numsh};
2574 nc_get_vara(ncid, varid, start.data(), count.data(),
2575 lmns_c2_buffer[i].data());
2576 }
2577 nc_inq_varid(ncid, "lmns_c3", &varid);
2578 for (size_t i = 0; i < nummn; i++) {
2579 const array<size_t, 2> start = {i, 0};
2580 const array<size_t, 2> count = {1, numsh};
2581 nc_get_vara(ncid, varid, start.data(), count.data(),
2582 lmns_c3_buffer[i].data());
2583 }
2584
2585 std::vector<double> xm_buffer(nummn);
2586 nc_inq_varid(ncid, "xm", &varid);
2587 nc_get_var(ncid, varid, xm_buffer.data());
2588
2589 std::vector<double> xn_buffer(nummn);
2590 nc_inq_varid(ncid, "xn", &varid);
2591 nc_get_var(ncid, varid, xn_buffer.data());
2592
2593 nc_close(ncid);
2594 output::sync.unlock();
2595
2596 auto sminf = static_cast<T> (sminf_value);
2597 auto sminh = static_cast<T> (sminh_value);
2598 auto ds = static_cast<T> (ds_value);
2599 auto dphi = graph::constant<T, SAFE_MATH> (static_cast<T> (dphi_value));
2600 auto signj = graph::constant<T, SAFE_MATH> (static_cast<T> (signj_value));
2601
2602 const backend::buffer<T> chi_c0(std::vector<T> (chi_c0_buffer.begin(), chi_c0_buffer.end()));
2603 const backend::buffer<T> chi_c1(std::vector<T> (chi_c1_buffer.begin(), chi_c1_buffer.end()));
2604 const backend::buffer<T> chi_c2(std::vector<T> (chi_c2_buffer.begin(), chi_c2_buffer.end()));
2605 const backend::buffer<T> chi_c3(std::vector<T> (chi_c3_buffer.begin(), chi_c3_buffer.end()));
2606
2607 std::vector<backend::buffer<T>> rmnc_c0(nummn);
2608 std::vector<backend::buffer<T>> rmnc_c1(nummn);
2609 std::vector<backend::buffer<T>> rmnc_c2(nummn);
2610 std::vector<backend::buffer<T>> rmnc_c3(nummn);
2611
2612 std::vector<backend::buffer<T>> zmns_c0(nummn);
2613 std::vector<backend::buffer<T>> zmns_c1(nummn);
2614 std::vector<backend::buffer<T>> zmns_c2(nummn);
2615 std::vector<backend::buffer<T>> zmns_c3(nummn);
2616
2617 std::vector<backend::buffer<T>> lmns_c0(nummn);
2618 std::vector<backend::buffer<T>> lmns_c1(nummn);
2619 std::vector<backend::buffer<T>> lmns_c2(nummn);
2620 std::vector<backend::buffer<T>> lmns_c3(nummn);
2621
2622 for (size_t i = 0; i < nummn; i++) {
2623 rmnc_c0[i] = backend::buffer(std::vector<T> (rmnc_c0_buffer[i].begin(), rmnc_c0_buffer[i].end()));
2624 rmnc_c1[i] = backend::buffer(std::vector<T> (rmnc_c1_buffer[i].begin(), rmnc_c1_buffer[i].end()));
2625 rmnc_c2[i] = backend::buffer(std::vector<T> (rmnc_c2_buffer[i].begin(), rmnc_c2_buffer[i].end()));
2626 rmnc_c3[i] = backend::buffer(std::vector<T> (rmnc_c3_buffer[i].begin(), rmnc_c3_buffer[i].end()));
2627
2628 zmns_c0[i] = backend::buffer(std::vector<T> (zmns_c0_buffer[i].begin(), zmns_c0_buffer[i].end()));
2629 zmns_c1[i] = backend::buffer(std::vector<T> (zmns_c1_buffer[i].begin(), zmns_c1_buffer[i].end()));
2630 zmns_c2[i] = backend::buffer(std::vector<T> (zmns_c2_buffer[i].begin(), zmns_c2_buffer[i].end()));
2631 zmns_c3[i] = backend::buffer(std::vector<T> (zmns_c3_buffer[i].begin(), zmns_c3_buffer[i].end()));
2632
2633 lmns_c0[i] = backend::buffer(std::vector<T> (lmns_c0_buffer[i].begin(), lmns_c0_buffer[i].end()));
2634 lmns_c1[i] = backend::buffer(std::vector<T> (lmns_c1_buffer[i].begin(), lmns_c1_buffer[i].end()));
2635 lmns_c2[i] = backend::buffer(std::vector<T> (lmns_c2_buffer[i].begin(), lmns_c2_buffer[i].end()));
2636 lmns_c3[i] = backend::buffer(std::vector<T> (lmns_c3_buffer[i].begin(), lmns_c3_buffer[i].end()));
2637 }
2638
2639 const backend::buffer<T> xm(std::vector<T> (xm_buffer.begin(), xm_buffer.end()));
2640 const backend::buffer<T> xn(std::vector<T> (xn_buffer.begin(), xn_buffer.end()));
2641
2642 return std::make_shared<vmec<T, SAFE_MATH>> (sminh, sminf, ds, dphi, signj,
2643 chi_c0, chi_c1, chi_c2, chi_c3,
2644 rmnc_c0, rmnc_c1, rmnc_c2, rmnc_c3,
2645 zmns_c0, zmns_c1, zmns_c2, zmns_c3,
2646 lmns_c0, lmns_c1, lmns_c2, lmns_c3,
2647 xm, xn);
2648 }
2649}
2650
2651#endif /* equilibrium_h */
arithmetic.hpp
Basic arithmetic operations.
backend::buffer
Class representing a generic buffer.
Definition backend.hpp:29
backend::buffer::size
size_t size() const
Get size of the buffer.
Definition backend.hpp:116
equilibrium::efit
2D EFIT equilibrium.
Definition equilibrium.hpp:1146
equilibrium::efit::get_ion_density
virtual graph::shared_leaf< T, SAFE_MATH > get_ion_density(const size_t index, graph::shared_leaf< T, SAFE_MATH > x, graph::shared_leaf< T, SAFE_MATH > y, graph::shared_leaf< T, SAFE_MATH > z)
Get the ion density.
Definition equilibrium.hpp:1518
equilibrium::efit::get_magnetic_field
virtual graph::shared_vector< T, SAFE_MATH > get_magnetic_field(graph::shared_leaf< T, SAFE_MATH > x, graph::shared_leaf< T, SAFE_MATH > y, graph::shared_leaf< T, SAFE_MATH > z)
Get the magnetic field.
Definition equilibrium.hpp:1569
equilibrium::efit::get_electron_temperature
virtual graph::shared_leaf< T, SAFE_MATH > get_electron_temperature(graph::shared_leaf< T, SAFE_MATH > x, graph::shared_leaf< T, SAFE_MATH > y, graph::shared_leaf< T, SAFE_MATH > z)
Get the electron temperature.
Definition equilibrium.hpp:1535
equilibrium::efit::efit
efit(const T psimin, const T dpsi, const backend::buffer< T > te_c0, const backend::buffer< T > te_c1, const backend::buffer< T > te_c2, const backend::buffer< T > te_c3, graph::shared_leaf< T, SAFE_MATH > te_scale, const backend::buffer< T > ne_c0, const backend::buffer< T > ne_c1, const backend::buffer< T > ne_c2, const backend::buffer< T > ne_c3, graph::shared_leaf< T, SAFE_MATH > ne_scale, const backend::buffer< T > pres_c0, const backend::buffer< T > pres_c1, const backend::buffer< T > pres_c2, const backend::buffer< T > pres_c3, graph::shared_leaf< T, SAFE_MATH > pres_scale, const T rmin, const T dr, const T zmin, const T dz, const backend::buffer< T > fpol_c0, const backend::buffer< T > fpol_c1, const backend::buffer< T > fpol_c2, const backend::buffer< T > fpol_c3, const size_t num_cols, const backend::buffer< T > c00, const backend::buffer< T > c01, const backend::buffer< T > c02, const backend::buffer< T > c03, const backend::buffer< T > c10, const backend::buffer< T > c11, const backend::buffer< T > c12, const backend::buffer< T > c13, const backend::buffer< T > c20, const backend::buffer< T > c21, const backend::buffer< T > c22, const backend::buffer< T > c23, const backend::buffer< T > c30, const backend::buffer< T > c31, const backend::buffer< T > c32, const backend::buffer< T > c33)
Construct a EFIT equilibrium.
Definition equilibrium.hpp:1433
equilibrium::efit::get_electron_density
virtual graph::shared_leaf< T, SAFE_MATH > get_electron_density(graph::shared_leaf< T, SAFE_MATH > x, graph::shared_leaf< T, SAFE_MATH > y, graph::shared_leaf< T, SAFE_MATH > z)
Get the electron density.
Definition equilibrium.hpp:1501
equilibrium::efit::get_ion_temperature
virtual graph::shared_leaf< T, SAFE_MATH > get_ion_temperature(const size_t index, graph::shared_leaf< T, SAFE_MATH > x, graph::shared_leaf< T, SAFE_MATH > y, graph::shared_leaf< T, SAFE_MATH > z)
Get the ion temperature.
Definition equilibrium.hpp:1552
equilibrium::efit::get_characteristic_field
virtual graph::shared_leaf< T, SAFE_MATH > get_characteristic_field(const size_t device_number=0) final
Get the characteristic field.
Definition equilibrium.hpp:1585
equilibrium::gaussian_density
Gaussian density with uniform magnetic field equilibrium.
Definition equilibrium.hpp:991
equilibrium::gaussian_density::get_ion_temperature
virtual graph::shared_leaf< T, SAFE_MATH > get_ion_temperature(const size_t index, graph::shared_leaf< T, SAFE_MATH > x, graph::shared_leaf< T, SAFE_MATH > y, graph::shared_leaf< T, SAFE_MATH > z) final
Get the ion temperature.
Definition equilibrium.hpp:1058
equilibrium::gaussian_density::gaussian_density
gaussian_density()
Construct a gaussian density with uniform magnetic field.
Definition equilibrium.hpp:996
equilibrium::gaussian_density::get_ion_density
virtual graph::shared_leaf< T, SAFE_MATH > get_ion_density(const size_t index, graph::shared_leaf< T, SAFE_MATH > x, graph::shared_leaf< T, SAFE_MATH > y, graph::shared_leaf< T, SAFE_MATH > z) final
Get the ion density.
Definition equilibrium.hpp:1025
equilibrium::gaussian_density::get_electron_temperature
virtual graph::shared_leaf< T, SAFE_MATH > get_electron_temperature(graph::shared_leaf< T, SAFE_MATH > x, graph::shared_leaf< T, SAFE_MATH > y, graph::shared_leaf< T, SAFE_MATH > z) final
Get the electron temperature.
Definition equilibrium.hpp:1042
equilibrium::gaussian_density::get_magnetic_field
virtual graph::shared_vector< T, SAFE_MATH > get_magnetic_field(graph::shared_leaf< T, SAFE_MATH > x, graph::shared_leaf< T, SAFE_MATH > y, graph::shared_leaf< T, SAFE_MATH > z) final
Get the magnetic field.
Definition equilibrium.hpp:1074
equilibrium::gaussian_density::get_characteristic_field
virtual graph::shared_leaf< T, SAFE_MATH > get_characteristic_field(const size_t device_number=0) final
Get the characteristic field.
Definition equilibrium.hpp:1090
equilibrium::gaussian_density::get_electron_density
virtual graph::shared_leaf< T, SAFE_MATH > get_electron_density(graph::shared_leaf< T, SAFE_MATH > x, graph::shared_leaf< T, SAFE_MATH > y, graph::shared_leaf< T, SAFE_MATH > z) final
Get the electron density.
Definition equilibrium.hpp:1008
equilibrium::generic
Class representing a generic equilibrium.
Definition equilibrium.hpp:236
equilibrium::generic::get_esup2
virtual graph::shared_vector< T, SAFE_MATH > get_esup2(graph::shared_leaf< T, SAFE_MATH > x1, graph::shared_leaf< T, SAFE_MATH > x2, graph::shared_leaf< T, SAFE_MATH > x3)
Get the contravariant basis vector in the x2 direction.
Definition equilibrium.hpp:397
equilibrium::generic::get_x
virtual graph::shared_leaf< T, SAFE_MATH > get_x(graph::shared_leaf< T, SAFE_MATH > x1, graph::shared_leaf< T, SAFE_MATH > x2, graph::shared_leaf< T, SAFE_MATH > x3)
Get the x position.
Definition equilibrium.hpp:431
equilibrium::generic::get_y
virtual graph::shared_leaf< T, SAFE_MATH > get_y(graph::shared_leaf< T, SAFE_MATH > x1, graph::shared_leaf< T, SAFE_MATH > x2, graph::shared_leaf< T, SAFE_MATH > x3)
Get the y position.
Definition equilibrium.hpp:446
equilibrium::generic::ion_charges
const std::vector< uint8_t > ion_charges
Ion charge for each species.
Definition equilibrium.hpp:241
equilibrium::generic::~generic
virtual ~generic()
Destructor.
Definition equilibrium.hpp:260
equilibrium::generic::get_ion_temperature
virtual graph::shared_leaf< T, SAFE_MATH > get_ion_temperature(const size_t index, graph::shared_leaf< T, SAFE_MATH > x, graph::shared_leaf< T, SAFE_MATH > y, graph::shared_leaf< T, SAFE_MATH > z)=0
Get the ion temperature.
equilibrium::generic::get_electron_density
virtual graph::shared_leaf< T, SAFE_MATH > get_electron_density(graph::shared_leaf< T, SAFE_MATH > x, graph::shared_leaf< T, SAFE_MATH > y, graph::shared_leaf< T, SAFE_MATH > z)=0
Get the electron density.
equilibrium::generic::get_esup3
virtual graph::shared_vector< T, SAFE_MATH > get_esup3(graph::shared_leaf< T, SAFE_MATH > x1, graph::shared_leaf< T, SAFE_MATH > x2, graph::shared_leaf< T, SAFE_MATH > x3)
Get the contravariant basis vector in the x3 direction.
Definition equilibrium.hpp:414
equilibrium::generic::get_z
virtual graph::shared_leaf< T, SAFE_MATH > get_z(graph::shared_leaf< T, SAFE_MATH > x1, graph::shared_leaf< T, SAFE_MATH > x2, graph::shared_leaf< T, SAFE_MATH > x3)
Get the z position.
Definition equilibrium.hpp:461
equilibrium::generic::get_electron_temperature
virtual graph::shared_leaf< T, SAFE_MATH > get_electron_temperature(graph::shared_leaf< T, SAFE_MATH > x, graph::shared_leaf< T, SAFE_MATH > y, graph::shared_leaf< T, SAFE_MATH > z)=0
Get the electron temperature.
equilibrium::generic::get_ion_charge
uint8_t get_ion_charge(const size_t index) const
Get the charge for an ion species.
Definition equilibrium.hpp:287
equilibrium::generic::get_characteristic_field
virtual graph::shared_leaf< T, SAFE_MATH > get_characteristic_field(const size_t device_number=0)=0
Get the characteristic field.
equilibrium::generic::get_esup1
virtual graph::shared_vector< T, SAFE_MATH > get_esup1(graph::shared_leaf< T, SAFE_MATH > x1, graph::shared_leaf< T, SAFE_MATH > x2, graph::shared_leaf< T, SAFE_MATH > x3)
Get the contravariant basis vector in the x1 direction.
Definition equilibrium.hpp:380
equilibrium::generic::get_magnetic_field
virtual graph::shared_vector< T, SAFE_MATH > get_magnetic_field(graph::shared_leaf< T, SAFE_MATH > x, graph::shared_leaf< T, SAFE_MATH > y, graph::shared_leaf< T, SAFE_MATH > z)=0
Get the magnetic field.
equilibrium::generic::get_num_ion_species
size_t get_num_ion_species() const
Get the number of ion species.
Definition equilibrium.hpp:267
equilibrium::generic::get_ion_mass
T get_ion_mass(const size_t index) const
Get the mass for an ion species.
Definition equilibrium.hpp:277
equilibrium::generic::get_ion_density
virtual graph::shared_leaf< T, SAFE_MATH > get_ion_density(const size_t index, graph::shared_leaf< T, SAFE_MATH > x, graph::shared_leaf< T, SAFE_MATH > y, graph::shared_leaf< T, SAFE_MATH > z)=0
Get the ion density.
equilibrium::generic::ion_masses
const std::vector< T > ion_masses
Ion masses for each species.
Definition equilibrium.hpp:239
equilibrium::no_magnetic_field
Uniform density with no magnetic field equilibrium.
Definition equilibrium.hpp:482
equilibrium::no_magnetic_field::get_magnetic_field
virtual graph::shared_vector< T, SAFE_MATH > get_magnetic_field(graph::shared_leaf< T, SAFE_MATH > x, graph::shared_leaf< T, SAFE_MATH > y, graph::shared_leaf< T, SAFE_MATH > z) final
Get the magnetic field.
Definition equilibrium.hpp:567
equilibrium::no_magnetic_field::no_magnetic_field
no_magnetic_field()
Construct a linear density with no magnetic field.
Definition equilibrium.hpp:487
equilibrium::no_magnetic_field::get_ion_temperature
virtual graph::shared_leaf< T, SAFE_MATH > get_ion_temperature(const size_t index, graph::shared_leaf< T, SAFE_MATH > x, graph::shared_leaf< T, SAFE_MATH > y, graph::shared_leaf< T, SAFE_MATH > z) final
Get the ion temperature.
Definition equilibrium.hpp:551
equilibrium::no_magnetic_field::get_ion_density
virtual graph::shared_leaf< T, SAFE_MATH > get_ion_density(const size_t index, graph::shared_leaf< T, SAFE_MATH > x, graph::shared_leaf< T, SAFE_MATH > y, graph::shared_leaf< T, SAFE_MATH > z) final
Get the ion density.
Definition equilibrium.hpp:517
equilibrium::no_magnetic_field::get_electron_temperature
virtual graph::shared_leaf< T, SAFE_MATH > get_electron_temperature(graph::shared_leaf< T, SAFE_MATH > x, graph::shared_leaf< T, SAFE_MATH > y, graph::shared_leaf< T, SAFE_MATH > z) final
Get the electron temperature.
Definition equilibrium.hpp:535
equilibrium::no_magnetic_field::get_electron_density
virtual graph::shared_leaf< T, SAFE_MATH > get_electron_density(graph::shared_leaf< T, SAFE_MATH > x, graph::shared_leaf< T, SAFE_MATH > y, graph::shared_leaf< T, SAFE_MATH > z) final
Get the electron density.
Definition equilibrium.hpp:499
equilibrium::no_magnetic_field::get_characteristic_field
virtual graph::shared_leaf< T, SAFE_MATH > get_characteristic_field(const size_t device_number=0) final
Get the characteristic field.
Definition equilibrium.hpp:583
equilibrium::slab_density
Vary density with uniform magnetic field equilibrium.
Definition equilibrium.hpp:735
equilibrium::slab_density::get_ion_temperature
virtual graph::shared_leaf< T, SAFE_MATH > get_ion_temperature(const size_t index, graph::shared_leaf< T, SAFE_MATH > x, graph::shared_leaf< T, SAFE_MATH > y, graph::shared_leaf< T, SAFE_MATH > z) final
Get the ion temperature.
Definition equilibrium.hpp:804
equilibrium::slab_density::get_magnetic_field
virtual graph::shared_vector< T, SAFE_MATH > get_magnetic_field(graph::shared_leaf< T, SAFE_MATH > x, graph::shared_leaf< T, SAFE_MATH > y, graph::shared_leaf< T, SAFE_MATH > z) final
Get the magnetic field.
Definition equilibrium.hpp:820
equilibrium::slab_density::get_characteristic_field
virtual graph::shared_leaf< T, SAFE_MATH > get_characteristic_field(const size_t device_number=0) final
Get the characteristic field.
Definition equilibrium.hpp:836
equilibrium::slab_density::get_electron_temperature
virtual graph::shared_leaf< T, SAFE_MATH > get_electron_temperature(graph::shared_leaf< T, SAFE_MATH > x, graph::shared_leaf< T, SAFE_MATH > y, graph::shared_leaf< T, SAFE_MATH > z) final
Get the electron temperature.
Definition equilibrium.hpp:788
equilibrium::slab_density::get_ion_density
virtual graph::shared_leaf< T, SAFE_MATH > get_ion_density(const size_t index, graph::shared_leaf< T, SAFE_MATH > x, graph::shared_leaf< T, SAFE_MATH > y, graph::shared_leaf< T, SAFE_MATH > z) final
Get the ion density.
Definition equilibrium.hpp:770
equilibrium::slab_density::get_electron_density
virtual graph::shared_leaf< T, SAFE_MATH > get_electron_density(graph::shared_leaf< T, SAFE_MATH > x, graph::shared_leaf< T, SAFE_MATH > y, graph::shared_leaf< T, SAFE_MATH > z) final
Get the electron density.
Definition equilibrium.hpp:752
equilibrium::slab_density::slab_density
slab_density()
Construct a gaussian density with uniform magnetic field.
Definition equilibrium.hpp:740
equilibrium::slab_field
Vary density with uniform magnetic field equilibrium.
Definition equilibrium.hpp:864
equilibrium::slab_field::get_electron_density
virtual graph::shared_leaf< T, SAFE_MATH > get_electron_density(graph::shared_leaf< T, SAFE_MATH > x, graph::shared_leaf< T, SAFE_MATH > y, graph::shared_leaf< T, SAFE_MATH > z) final
Get the electron density.
Definition equilibrium.hpp:881
equilibrium::slab_field::get_ion_density
virtual graph::shared_leaf< T, SAFE_MATH > get_ion_density(const size_t index, graph::shared_leaf< T, SAFE_MATH > x, graph::shared_leaf< T, SAFE_MATH > y, graph::shared_leaf< T, SAFE_MATH > z) final
Get the ion density.
Definition equilibrium.hpp:899
equilibrium::slab_field::get_electron_temperature
virtual graph::shared_leaf< T, SAFE_MATH > get_electron_temperature(graph::shared_leaf< T, SAFE_MATH > x, graph::shared_leaf< T, SAFE_MATH > y, graph::shared_leaf< T, SAFE_MATH > z) final
Get the electron temperature.
Definition equilibrium.hpp:915
equilibrium::slab_field::slab_field
slab_field()
Construct a gaussian density with uniform magnetic field.
Definition equilibrium.hpp:869
equilibrium::slab_field::get_characteristic_field
virtual graph::shared_leaf< T, SAFE_MATH > get_characteristic_field(const size_t device_number=0) final
Get the characteristic field.
Definition equilibrium.hpp:964
equilibrium::slab_field::get_magnetic_field
virtual graph::shared_vector< T, SAFE_MATH > get_magnetic_field(graph::shared_leaf< T, SAFE_MATH > x, graph::shared_leaf< T, SAFE_MATH > y, graph::shared_leaf< T, SAFE_MATH > z) final
Get the magnetic field.
Definition equilibrium.hpp:949
equilibrium::slab_field::get_ion_temperature
virtual graph::shared_leaf< T, SAFE_MATH > get_ion_temperature(const size_t index, graph::shared_leaf< T, SAFE_MATH > x, graph::shared_leaf< T, SAFE_MATH > y, graph::shared_leaf< T, SAFE_MATH > z) final
Get the ion temperature.
Definition equilibrium.hpp:933
equilibrium::slab
Uniform density with varying magnetic field equilibrium.
Definition equilibrium.hpp:611
equilibrium::slab::get_characteristic_field
virtual graph::shared_leaf< T, SAFE_MATH > get_characteristic_field(const size_t device_number=0) final
Get the characteristic field.
Definition equilibrium.hpp:707
equilibrium::slab::get_electron_temperature
virtual graph::shared_leaf< T, SAFE_MATH > get_electron_temperature(graph::shared_leaf< T, SAFE_MATH > x, graph::shared_leaf< T, SAFE_MATH > y, graph::shared_leaf< T, SAFE_MATH > z) final
Get the electron temperature.
Definition equilibrium.hpp:660
equilibrium::slab::slab
slab()
Construct a gaussian density with uniform magnetic field.
Definition equilibrium.hpp:616
equilibrium::slab::get_magnetic_field
virtual graph::shared_vector< T, SAFE_MATH > get_magnetic_field(graph::shared_leaf< T, SAFE_MATH > x, graph::shared_leaf< T, SAFE_MATH > y, graph::shared_leaf< T, SAFE_MATH > z) final
Get the magnetic field.
Definition equilibrium.hpp:692
equilibrium::slab::get_ion_density
virtual graph::shared_leaf< T, SAFE_MATH > get_ion_density(const size_t index, graph::shared_leaf< T, SAFE_MATH > x, graph::shared_leaf< T, SAFE_MATH > y, graph::shared_leaf< T, SAFE_MATH > z) final
Get the ion density.
Definition equilibrium.hpp:644
equilibrium::slab::get_electron_density
virtual graph::shared_leaf< T, SAFE_MATH > get_electron_density(graph::shared_leaf< T, SAFE_MATH > x, graph::shared_leaf< T, SAFE_MATH > y, graph::shared_leaf< T, SAFE_MATH > z) final
Get the electron density.
Definition equilibrium.hpp:628
equilibrium::slab::get_ion_temperature
virtual graph::shared_leaf< T, SAFE_MATH > get_ion_temperature(const size_t index, graph::shared_leaf< T, SAFE_MATH > x, graph::shared_leaf< T, SAFE_MATH > y, graph::shared_leaf< T, SAFE_MATH > z) final
Get the ion temperature.
Definition equilibrium.hpp:676
equilibrium::vmec
3D VMEC equilibrium.
Definition equilibrium.hpp:1868
equilibrium::vmec::get_electron_density
virtual graph::shared_leaf< T, SAFE_MATH > get_electron_density(graph::shared_leaf< T, SAFE_MATH > s, graph::shared_leaf< T, SAFE_MATH > u, graph::shared_leaf< T, SAFE_MATH > v)
Get the electron density.
Definition equilibrium.hpp:2276
equilibrium::vmec::get_esup3
virtual graph::shared_vector< T, SAFE_MATH > get_esup3(graph::shared_leaf< T, SAFE_MATH > s, graph::shared_leaf< T, SAFE_MATH > u, graph::shared_leaf< T, SAFE_MATH > v)
Get the contravariant basis vector in the V direction.
Definition equilibrium.hpp:2260
equilibrium::vmec::get_x
virtual graph::shared_leaf< T, SAFE_MATH > get_x(graph::shared_leaf< T, SAFE_MATH > s, graph::shared_leaf< T, SAFE_MATH > u, graph::shared_leaf< T, SAFE_MATH > v)
Get the x position.
Definition equilibrium.hpp:2375
equilibrium::vmec::get_y
virtual graph::shared_leaf< T, SAFE_MATH > get_y(graph::shared_leaf< T, SAFE_MATH > s, graph::shared_leaf< T, SAFE_MATH > u, graph::shared_leaf< T, SAFE_MATH > v)
Get the y position.
Definition equilibrium.hpp:2391
equilibrium::vmec::get_characteristic_field
virtual graph::shared_leaf< T, SAFE_MATH > get_characteristic_field(const size_t device_number=0) final
Get the characteristic field.
Definition equilibrium.hpp:2358
equilibrium::vmec::get_esup2
virtual graph::shared_vector< T, SAFE_MATH > get_esup2(graph::shared_leaf< T, SAFE_MATH > s, graph::shared_leaf< T, SAFE_MATH > u, graph::shared_leaf< T, SAFE_MATH > v)
Get the contravariant basis vector in the U direction.
Definition equilibrium.hpp:2244
equilibrium::vmec::get_electron_temperature
virtual graph::shared_leaf< T, SAFE_MATH > get_electron_temperature(graph::shared_leaf< T, SAFE_MATH > s, graph::shared_leaf< T, SAFE_MATH > u, graph::shared_leaf< T, SAFE_MATH > v)
Get the electron temperature.
Definition equilibrium.hpp:2309
equilibrium::vmec::get_esup1
virtual graph::shared_vector< T, SAFE_MATH > get_esup1(graph::shared_leaf< T, SAFE_MATH > s, graph::shared_leaf< T, SAFE_MATH > u, graph::shared_leaf< T, SAFE_MATH > v)
Get the contravariant basis vector in the S direction.
Definition equilibrium.hpp:2228
equilibrium::vmec::get_ion_temperature
virtual graph::shared_leaf< T, SAFE_MATH > get_ion_temperature(const size_t index, graph::shared_leaf< T, SAFE_MATH > s, graph::shared_leaf< T, SAFE_MATH > u, graph::shared_leaf< T, SAFE_MATH > v)
Get the ion temperature.
Definition equilibrium.hpp:2326
equilibrium::vmec::get_ion_density
virtual graph::shared_leaf< T, SAFE_MATH > get_ion_density(const size_t index, graph::shared_leaf< T, SAFE_MATH > s, graph::shared_leaf< T, SAFE_MATH > u, graph::shared_leaf< T, SAFE_MATH > v)
Get the ion density.
Definition equilibrium.hpp:2293
equilibrium::vmec::get_magnetic_field
virtual graph::shared_vector< T, SAFE_MATH > get_magnetic_field(graph::shared_leaf< T, SAFE_MATH > s, graph::shared_leaf< T, SAFE_MATH > u, graph::shared_leaf< T, SAFE_MATH > v)
Get the magnetic field.
Definition equilibrium.hpp:2342
equilibrium::vmec::get_z
virtual graph::shared_leaf< T, SAFE_MATH > get_z(graph::shared_leaf< T, SAFE_MATH > s, graph::shared_leaf< T, SAFE_MATH > u, graph::shared_leaf< T, SAFE_MATH > v)
Get the z position.
Definition equilibrium.hpp:2407
equilibrium::vmec::vmec
vmec(const T sminh, const T sminf, const T ds, graph::shared_leaf< T, SAFE_MATH > dphi, graph::shared_leaf< T, SAFE_MATH > signj, const backend::buffer< T > chi_c0, const backend::buffer< T > chi_c1, const backend::buffer< T > chi_c2, const backend::buffer< T > chi_c3, const std::vector< backend::buffer< T > > rmnc_c0, const std::vector< backend::buffer< T > > rmnc_c1, const std::vector< backend::buffer< T > > rmnc_c2, const std::vector< backend::buffer< T > > rmnc_c3, const std::vector< backend::buffer< T > > zmns_c0, const std::vector< backend::buffer< T > > zmns_c1, const std::vector< backend::buffer< T > > zmns_c2, const std::vector< backend::buffer< T > > zmns_c3, const std::vector< backend::buffer< T > > lmns_c0, const std::vector< backend::buffer< T > > lmns_c1, const std::vector< backend::buffer< T > > lmns_c2, const std::vector< backend::buffer< T > > lmns_c3, const backend::buffer< T > xm, const backend::buffer< T > xn)
Construct a EFIT equilibrium.
Definition equilibrium.hpp:2183
workflow::manager
Class representing a workflow manager.
Definition workflow.hpp:215
workflow::manager::run
void run()
Run work items.
Definition workflow.hpp:359
workflow::manager::copy_to_host
void copy_to_host(graph::shared_leaf< T, SAFE_MATH > &node, T *destination)
Copy contexts of buffer to host.
Definition workflow.hpp:389
workflow::manager::add_item
void add_item(graph::input_nodes< T, SAFE_MATH > in, graph::output_nodes< T, SAFE_MATH > out, graph::map_nodes< T, SAFE_MATH > maps, graph::shared_random_state< T, SAFE_MATH > state, const std::string name, const size_t size)
Add a workflow item.
Definition workflow.hpp:271
workflow::manager::compile
void compile()
Compile the workflow items.
Definition workflow.hpp:336
assert
subroutine assert(test, message)
Assert check.
Definition f_binding_test.f90:38
math.hpp
Defined basic math functions.
equilibrium
Name space for equilibrium models.
Definition equilibrium.hpp:225
equilibrium::make_slab_field
shared< T, SAFE_MATH > make_slab_field()
Convenience function to build a slab density equilibrium.
Definition equilibrium.hpp:978
equilibrium::make_gaussian_density
shared< T, SAFE_MATH > make_gaussian_density()
Convenience function to build a gaussian density equilibrium.
Definition equilibrium.hpp:1104
equilibrium::make_slab_density
shared< T, SAFE_MATH > make_slab_density()
Convenience function to build a slab density equilibrium.
Definition equilibrium.hpp:850
equilibrium::make_vmec
shared< T, SAFE_MATH > make_vmec(const std::string &spline_file)
Convenience function to build an VMEC equilibrium.
Definition equilibrium.hpp:2425
equilibrium::shared
std::shared_ptr< generic< T, SAFE_MATH > > shared
Convenience type alias for shared equilibria.
Definition equilibrium.hpp:470
equilibrium::make_no_magnetic_field
shared< T, SAFE_MATH > make_no_magnetic_field()
Convenience function to build a no magnetic field equilibrium.
Definition equilibrium.hpp:597
equilibrium::build_1D_spline
graph::shared_leaf< T, SAFE_MATH > build_1D_spline(graph::output_nodes< T, SAFE_MATH > c, graph::shared_leaf< T, SAFE_MATH > x, const T scale, const T offset)
Build a 1D spline.
Definition equilibrium.hpp:1121
equilibrium::make_slab
shared< T, SAFE_MATH > make_slab()
Convenience function to build a slab equilibrium.
Definition equilibrium.hpp:721
equilibrium::make_efit
shared< T, SAFE_MATH > make_efit(const std::string &spline_file)
Convenience function to build an EFIT equilibrium.
Definition equilibrium.hpp:1628
graph::zero
constexpr shared_leaf< T, SAFE_MATH > zero()
Forward declare for zero.
Definition node.hpp:986
graph::piecewise_1D
shared_leaf< T, SAFE_MATH > piecewise_1D(const backend::buffer< T > &d, shared_leaf< T, SAFE_MATH > x, const T scale, const T offset)
Define piecewise_1D convenience function.
Definition piecewise.hpp:570
graph::atan
shared_leaf< T, SAFE_MATH > atan(shared_leaf< T, SAFE_MATH > l, shared_leaf< T, SAFE_MATH > r)
Build arctan node.
Definition trigonometry.hpp:808
graph::shared_vector
std::shared_ptr< vector_quantity< T, SAFE_MATH > > shared_vector
Convenience type for shared vector quantities.
Definition vector.hpp:142
graph::pow
shared_leaf< T, SAFE_MATH > pow(shared_leaf< T, SAFE_MATH > l, shared_leaf< T, SAFE_MATH > r)
Build power node.
Definition math.hpp:1352
graph::piecewise_2D
shared_leaf< T, SAFE_MATH > piecewise_2D(const backend::buffer< T > &d, const size_t n, shared_leaf< T, SAFE_MATH > x, const T x_scale, const T x_offset, shared_leaf< T, SAFE_MATH > y, const T y_scale, const T y_offset)
Define piecewise_2D convenience function.
Definition piecewise.hpp:1312
graph::vector
shared_vector< T, SAFE_MATH > vector(shared_leaf< T, SAFE_MATH > x, shared_leaf< T, SAFE_MATH > y, shared_leaf< T, SAFE_MATH > z)
Build a shared vector quantity.
Definition vector.hpp:156
graph::exp
shared_leaf< T, SAFE_MATH > exp(shared_leaf< T, SAFE_MATH > x)
Define exp convenience function.
Definition math.hpp:544
graph::shared_random_state
std::shared_ptr< random_state_node< T, SAFE_MATH > > shared_random_state
Convenience type alias for shared sqrt nodes.
Definition random.hpp:263
graph::input_nodes
std::vector< shared_variable< T, SAFE_MATH > > input_nodes
Convenience type alias for a vector of inputs.
Definition node.hpp:1711
graph::fma
shared_leaf< T, SAFE_MATH > fma(shared_leaf< T, SAFE_MATH > l, shared_leaf< T, SAFE_MATH > m, shared_leaf< T, SAFE_MATH > r)
Build fused multiply add node.
Definition arithmetic.hpp:5202
graph::variable_cast
shared_variable< T, SAFE_MATH > variable_cast(shared_leaf< T, SAFE_MATH > x)
Cast to a variable node.
Definition node.hpp:1727
graph::sin
shared_leaf< T, SAFE_MATH > sin(shared_leaf< T, SAFE_MATH > x)
Define sine convenience function.
Definition trigonometry.hpp:229
graph::sqrt
shared_leaf< T, SAFE_MATH > sqrt(shared_leaf< T, SAFE_MATH > x)
Define sqrt convenience function.
Definition math.hpp:279
graph::matrix
shared_matrix< T, SAFE_MATH > matrix(shared_vector< T, SAFE_MATH > r1, shared_vector< T, SAFE_MATH > r2, shared_vector< T, SAFE_MATH > r3)
Build a shared vector quantity.
Definition vector.hpp:436
graph::shared_leaf
std::shared_ptr< leaf_node< T, SAFE_MATH > > shared_leaf
Convenience type alias for shared leaf nodes.
Definition node.hpp:676
graph::cos
shared_leaf< T, SAFE_MATH > cos(shared_leaf< T, SAFE_MATH > x)
Define cosine convenience function.
Definition trigonometry.hpp:487
graph::output_nodes
std::vector< shared_leaf< T, SAFE_MATH > > output_nodes
Convenience type alias for a vector of output nodes.
Definition node.hpp:691
solver::newton
void newton(workflow::manager< T, SAFE_MATH > &work, graph::output_nodes< T, SAFE_MATH > vars, graph::input_nodes< T, SAFE_MATH > inputs, graph::shared_leaf< T, SAFE_MATH > func, graph::shared_random_state< T, SAFE_MATH > state, const T tolerance=1.0E-30, const size_t max_iterations=1000, const T step=1.0)
Determine the value of vars to minimize the loss function.
Definition newton.hpp:34
newton.hpp
Sets up the kernel for a newtons method.
output.hpp
Implements output files in a netcdf format.
piecewise.hpp
Piecewise nodes.
trigonometry.hpp
Trigonometry functions.
vector.hpp
Defines vectors of graphs.