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
224namespace equilibrium {
226 static std::mutex sync;
227
228//******************************************************************************
229// Equilibrium interface
230//******************************************************************************
231//------------------------------------------------------------------------------
236//------------------------------------------------------------------------------
237 template<jit::float_scalar T, bool SAFE_MATH=false>
238 class generic {
239 protected:
241 const std::vector<T> ion_masses;
243 const std::vector<uint8_t> ion_charges;
244
245 public:
246//------------------------------------------------------------------------------
251//------------------------------------------------------------------------------
252 generic(const std::vector<T> &masses,
253 const std::vector<uint8_t> &charges) :
254 ion_masses(masses), ion_charges(charges) {
255 assert(ion_masses.size() == ion_charges.size() &&
256 "Masses and charges need the same number of elements.");
257 }
258
259//------------------------------------------------------------------------------
261//------------------------------------------------------------------------------
262 virtual ~generic() {}
263
264//------------------------------------------------------------------------------
268//------------------------------------------------------------------------------
269 size_t get_num_ion_species() const {
270 return ion_masses.size();
271 }
272
273//------------------------------------------------------------------------------
278//------------------------------------------------------------------------------
279 T get_ion_mass(const size_t index) const {
280 return ion_masses.at(index);
281 }
282
283//------------------------------------------------------------------------------
288//------------------------------------------------------------------------------
289 uint8_t get_ion_charge(const size_t index) const {
290 return ion_charges.at(index);
291 }
292
293//------------------------------------------------------------------------------
300//------------------------------------------------------------------------------
301 virtual graph::shared_leaf<T, SAFE_MATH>
302 get_electron_density(graph::shared_leaf<T, SAFE_MATH> x,
303 graph::shared_leaf<T, SAFE_MATH> y,
304 graph::shared_leaf<T, SAFE_MATH> z) = 0;
305
306//------------------------------------------------------------------------------
314//------------------------------------------------------------------------------
315 virtual graph::shared_leaf<T, SAFE_MATH>
316 get_ion_density(const size_t index,
317 graph::shared_leaf<T, SAFE_MATH> x,
318 graph::shared_leaf<T, SAFE_MATH> y,
319 graph::shared_leaf<T, SAFE_MATH> z) = 0;
320
321//------------------------------------------------------------------------------
328//------------------------------------------------------------------------------
329 virtual graph::shared_leaf<T, SAFE_MATH>
330 get_electron_temperature(graph::shared_leaf<T, SAFE_MATH> x,
331 graph::shared_leaf<T, SAFE_MATH> y,
332 graph::shared_leaf<T, SAFE_MATH> z) = 0;
333
334//------------------------------------------------------------------------------
342//------------------------------------------------------------------------------
343 virtual graph::shared_leaf<T, SAFE_MATH>
344 get_ion_temperature(const size_t index,
345 graph::shared_leaf<T, SAFE_MATH> x,
346 graph::shared_leaf<T, SAFE_MATH> y,
347 graph::shared_leaf<T, SAFE_MATH> z) = 0;
348
349//------------------------------------------------------------------------------
356//------------------------------------------------------------------------------
357 virtual graph::shared_vector<T, SAFE_MATH>
358 get_magnetic_field(graph::shared_leaf<T, SAFE_MATH> x,
359 graph::shared_leaf<T, SAFE_MATH> y,
360 graph::shared_leaf<T, SAFE_MATH> z) = 0;
361
362//------------------------------------------------------------------------------
369//------------------------------------------------------------------------------
370 virtual graph::shared_leaf<T, SAFE_MATH>
371 get_characteristic_field(const size_t device_number=0) = 0;
372
373//------------------------------------------------------------------------------
380//------------------------------------------------------------------------------
381 virtual graph::shared_vector<T, SAFE_MATH>
389
390//------------------------------------------------------------------------------
397//------------------------------------------------------------------------------
398 virtual graph::shared_vector<T, SAFE_MATH>
406
407//------------------------------------------------------------------------------
414//------------------------------------------------------------------------------
415 virtual graph::shared_vector<T, SAFE_MATH>
423
424//------------------------------------------------------------------------------
431//------------------------------------------------------------------------------
432 virtual graph::shared_leaf<T, SAFE_MATH>
438
439//------------------------------------------------------------------------------
446//------------------------------------------------------------------------------
447 virtual graph::shared_leaf<T, SAFE_MATH>
453
454//------------------------------------------------------------------------------
461//------------------------------------------------------------------------------
462 virtual graph::shared_leaf<T, SAFE_MATH>
468 };
469
471 template<jit::float_scalar T, bool SAFE_MATH=false>
472 using shared = std::shared_ptr<generic<T, SAFE_MATH>>;
473
474//******************************************************************************
475// No Magnetic equilibrium.
476//******************************************************************************
477//------------------------------------------------------------------------------
482//------------------------------------------------------------------------------
483 template<jit::float_scalar T, bool SAFE_MATH=false>
484 class no_magnetic_field : public generic<T, SAFE_MATH> {
485 public:
486//------------------------------------------------------------------------------
488//------------------------------------------------------------------------------
489 no_magnetic_field() :
490 generic<T, SAFE_MATH> ({3.34449469E-27}, {1}) {}
491
492//------------------------------------------------------------------------------
499//------------------------------------------------------------------------------
500 virtual graph::shared_leaf<T, SAFE_MATH>
508
509//------------------------------------------------------------------------------
517//------------------------------------------------------------------------------
518 virtual graph::shared_leaf<T, SAFE_MATH>
519 get_ion_density(const size_t index,
520 graph::shared_leaf<T, SAFE_MATH> x,
521 graph::shared_leaf<T, SAFE_MATH> y,
522 graph::shared_leaf<T, SAFE_MATH> z) final {
523 return graph::constant<T, SAFE_MATH> (static_cast<T> (1.0E19)) *
524 (graph::constant<T, SAFE_MATH> (static_cast<T> (0.1))*x +
525 graph::one<T, SAFE_MATH> ());
526 }
527
528//------------------------------------------------------------------------------
535//------------------------------------------------------------------------------
536 virtual graph::shared_leaf<T, SAFE_MATH>
542
543//------------------------------------------------------------------------------
551//------------------------------------------------------------------------------
552 virtual graph::shared_leaf<T, SAFE_MATH>
553 get_ion_temperature(const size_t index,
554 graph::shared_leaf<T, SAFE_MATH> x,
555 graph::shared_leaf<T, SAFE_MATH> y,
556 graph::shared_leaf<T, SAFE_MATH> z) final {
557 return graph::constant<T, SAFE_MATH> (static_cast<T> (1000.0));
558 }
559
560//------------------------------------------------------------------------------
567//------------------------------------------------------------------------------
568 virtual graph::shared_vector<T, SAFE_MATH>
575
576//------------------------------------------------------------------------------
583//------------------------------------------------------------------------------
584 virtual graph::shared_leaf<T, SAFE_MATH>
585 get_characteristic_field(const size_t device_number=0) final {
586 return graph::one<T, SAFE_MATH> ();
587 }
588 };
589
590//------------------------------------------------------------------------------
597//------------------------------------------------------------------------------
598 template<jit::float_scalar T, bool SAFE_MATH=false>
599 shared<T, SAFE_MATH> make_no_magnetic_field() {
600 return std::make_shared<no_magnetic_field<T, SAFE_MATH>> ();
601 }
602
603//******************************************************************************
604// Slab equilibrium.
605//******************************************************************************
606//------------------------------------------------------------------------------
611//------------------------------------------------------------------------------
612 template<jit::float_scalar T, bool SAFE_MATH=false>
613 class slab : public generic<T, SAFE_MATH> {
614 public:
615//------------------------------------------------------------------------------
617//------------------------------------------------------------------------------
618 slab() :
619 generic<T, SAFE_MATH> ({3.34449469E-27}, {1}) {}
620
621//------------------------------------------------------------------------------
628//------------------------------------------------------------------------------
629 virtual graph::shared_leaf<T, SAFE_MATH>
635
636//------------------------------------------------------------------------------
644//------------------------------------------------------------------------------
645 virtual graph::shared_leaf<T, SAFE_MATH>
646 get_ion_density(const size_t index,
647 graph::shared_leaf<T, SAFE_MATH> x,
648 graph::shared_leaf<T, SAFE_MATH> y,
649 graph::shared_leaf<T, SAFE_MATH> z) final {
650 return graph::constant<T, SAFE_MATH> (static_cast<T> (1.0E19));
651 }
652
653//------------------------------------------------------------------------------
660//------------------------------------------------------------------------------
661 virtual graph::shared_leaf<T, SAFE_MATH>
667
668//------------------------------------------------------------------------------
676//------------------------------------------------------------------------------
677 virtual graph::shared_leaf<T, SAFE_MATH>
678 get_ion_temperature(const size_t index,
679 graph::shared_leaf<T, SAFE_MATH> x,
680 graph::shared_leaf<T, SAFE_MATH> y,
681 graph::shared_leaf<T, SAFE_MATH> z) final {
682 return graph::constant<T, SAFE_MATH> (static_cast<T> (1000.0));
683 }
684
685//------------------------------------------------------------------------------
692//------------------------------------------------------------------------------
693 virtual graph::shared_vector<T, SAFE_MATH>
699
700//------------------------------------------------------------------------------
707//------------------------------------------------------------------------------
708 virtual graph::shared_leaf<T, SAFE_MATH>
709 get_characteristic_field(const size_t device_number=0) final {
710 return graph::one<T, SAFE_MATH> ();
711 }
712 };
713
714//------------------------------------------------------------------------------
721//------------------------------------------------------------------------------
722 template<jit::float_scalar T, bool SAFE_MATH=false>
723 shared<T, SAFE_MATH> make_slab() {
724 return std::make_shared<slab<T, SAFE_MATH>> ();
725 }
726
727//******************************************************************************
728// Slab density equilibrium.
729//******************************************************************************
730//------------------------------------------------------------------------------
735//------------------------------------------------------------------------------
736 template<jit::float_scalar T, bool SAFE_MATH=false>
737 class slab_density : public generic<T, SAFE_MATH> {
738 public:
739//------------------------------------------------------------------------------
741//------------------------------------------------------------------------------
742 slab_density() :
743 generic<T, SAFE_MATH> ({3.34449469E-27}, {1}) {}
744
745//------------------------------------------------------------------------------
752//------------------------------------------------------------------------------
753 virtual graph::shared_leaf<T, SAFE_MATH>
761
762//------------------------------------------------------------------------------
770//------------------------------------------------------------------------------
771 virtual graph::shared_leaf<T, SAFE_MATH>
772 get_ion_density(const size_t index,
773 graph::shared_leaf<T, SAFE_MATH> x,
774 graph::shared_leaf<T, SAFE_MATH> y,
775 graph::shared_leaf<T, SAFE_MATH> z) final {
776 return graph::constant<T, SAFE_MATH> (static_cast<T> (1.0E19)) *
777 (graph::constant<T, SAFE_MATH> (static_cast<T> (0.1))*x +
778 graph::one<T, SAFE_MATH> ());
779 }
780
781//------------------------------------------------------------------------------
788//------------------------------------------------------------------------------
789 virtual graph::shared_leaf<T, SAFE_MATH>
795
796//------------------------------------------------------------------------------
804//------------------------------------------------------------------------------
805 virtual graph::shared_leaf<T, SAFE_MATH>
806 get_ion_temperature(const size_t index,
807 graph::shared_leaf<T, SAFE_MATH> x,
808 graph::shared_leaf<T, SAFE_MATH> y,
809 graph::shared_leaf<T, SAFE_MATH> z) final {
810 return graph::constant<T, SAFE_MATH> (static_cast<T> (1000.0));
811 }
812
813//------------------------------------------------------------------------------
820//------------------------------------------------------------------------------
821 virtual graph::shared_vector<T, SAFE_MATH>
828
829//------------------------------------------------------------------------------
836//------------------------------------------------------------------------------
837 virtual graph::shared_leaf<T, SAFE_MATH>
838 get_characteristic_field(const size_t device_number=0) final {
839 return graph::one<T, SAFE_MATH> ();
840 }
841 };
842
843//------------------------------------------------------------------------------
850//------------------------------------------------------------------------------
851 template<jit::float_scalar T, bool SAFE_MATH=false>
852 shared<T, SAFE_MATH> make_slab_density() {
853 return std::make_shared<slab_density<T, SAFE_MATH>> ();
854 }
855
856//******************************************************************************
857// Slab field gradient equilibrium.
858//******************************************************************************
859//------------------------------------------------------------------------------
864//------------------------------------------------------------------------------
865 template<jit::float_scalar T, bool SAFE_MATH=false>
866 class slab_field : public generic<T, SAFE_MATH> {
867 public:
868//------------------------------------------------------------------------------
870//------------------------------------------------------------------------------
871 slab_field() :
872 generic<T, SAFE_MATH> ({3.34449469E-27}, {1}) {}
873
874//------------------------------------------------------------------------------
881//------------------------------------------------------------------------------
882 virtual graph::shared_leaf<T, SAFE_MATH>
890
891//------------------------------------------------------------------------------
899//------------------------------------------------------------------------------
900 virtual graph::shared_leaf<T, SAFE_MATH>
907
908//------------------------------------------------------------------------------
915//------------------------------------------------------------------------------
916 virtual graph::shared_leaf<T, SAFE_MATH>
924
925//------------------------------------------------------------------------------
933//------------------------------------------------------------------------------
934 virtual graph::shared_leaf<T, SAFE_MATH>
941
942//------------------------------------------------------------------------------
949//------------------------------------------------------------------------------
950 virtual graph::shared_vector<T, SAFE_MATH>
956
957//------------------------------------------------------------------------------
964//------------------------------------------------------------------------------
965 virtual graph::shared_leaf<T, SAFE_MATH>
966 get_characteristic_field(const size_t device_number=0) final {
967 return graph::one<T, SAFE_MATH> ();
968 }
969 };
970
971//------------------------------------------------------------------------------
978//------------------------------------------------------------------------------
979 template<jit::float_scalar T, bool SAFE_MATH=false>
980 shared<T, SAFE_MATH> make_slab_field() {
981 return std::make_shared<slab_field<T, SAFE_MATH>> ();
982 }
983//******************************************************************************
984// Guassian density with a uniform magnetic field.
985//******************************************************************************
986//------------------------------------------------------------------------------
991//------------------------------------------------------------------------------
992 template<jit::float_scalar T, bool SAFE_MATH=false>
993 class guassian_density : public generic<T, SAFE_MATH> {
994 public:
995//------------------------------------------------------------------------------
997//------------------------------------------------------------------------------
998 guassian_density() :
999 generic<T, SAFE_MATH> ({3.34449469E-27}, {1}) {}
1000
1001//------------------------------------------------------------------------------
1008//------------------------------------------------------------------------------
1009 virtual graph::shared_leaf<T, SAFE_MATH>
1016
1017//------------------------------------------------------------------------------
1025//------------------------------------------------------------------------------
1026 virtual graph::shared_leaf<T, SAFE_MATH>
1027 get_ion_density(const size_t index,
1028 graph::shared_leaf<T, SAFE_MATH> x,
1029 graph::shared_leaf<T, SAFE_MATH> y,
1030 graph::shared_leaf<T, SAFE_MATH> z) final {
1031 return graph::constant<T, SAFE_MATH> (static_cast<T> (1.0E19)) *
1032 graph::exp((x*x + y*y)/graph::constant<T, SAFE_MATH> (static_cast<T> (-0.2)));
1033 }
1034
1035//------------------------------------------------------------------------------
1042//------------------------------------------------------------------------------
1043 virtual graph::shared_leaf<T, SAFE_MATH>
1049
1050//------------------------------------------------------------------------------
1058//------------------------------------------------------------------------------
1059 virtual graph::shared_leaf<T, SAFE_MATH>
1060 get_ion_temperature(const size_t index,
1061 graph::shared_leaf<T, SAFE_MATH> x,
1062 graph::shared_leaf<T, SAFE_MATH> y,
1063 graph::shared_leaf<T, SAFE_MATH> z) final {
1064 return graph::constant<T, SAFE_MATH> (static_cast<T> (1000.0));
1065 }
1066
1067//------------------------------------------------------------------------------
1074//------------------------------------------------------------------------------
1075 virtual graph::shared_vector<T, SAFE_MATH>
1082
1083//------------------------------------------------------------------------------
1090//------------------------------------------------------------------------------
1091 virtual graph::shared_leaf<T, SAFE_MATH>
1092 get_characteristic_field(const size_t device_number=0) final {
1093 return graph::one<T, SAFE_MATH> ();
1094 }
1095 };
1096
1097//------------------------------------------------------------------------------
1104//------------------------------------------------------------------------------
1105 template<jit::float_scalar T, bool SAFE_MATH=false>
1106 shared<T, SAFE_MATH> make_guassian_density() {
1107 return std::make_shared<guassian_density<T, SAFE_MATH>> ();
1108 }
1109
1110//------------------------------------------------------------------------------
1121//------------------------------------------------------------------------------
1122 template<jit::float_scalar T, bool SAFE_MATH=false>
1123 graph::shared_leaf<T, SAFE_MATH> build_1D_spline(graph::output_nodes<T, SAFE_MATH> c,
1124 graph::shared_leaf<T, SAFE_MATH> x,
1125 const T scale,
1126 const T offset) {
1127 auto c3 = c[3]/(scale*scale*scale);
1128 auto c2 = c[2]/(scale*scale) - static_cast<T> (3.0)*offset*c[3]/(scale*scale*scale);
1129 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);
1130 auto c0 = c[0] - offset*c[1]/scale + offset*offset*c[2]/(scale*scale) - offset*offset*offset*c[3]/(scale*scale*scale);
1131
1132 return graph::fma(graph::fma(graph::fma(c3, x, c2), x, c1), x, c0);
1133 }
1134
1135//******************************************************************************
1136// 2D EFIT equilibrium.
1137//******************************************************************************
1138//------------------------------------------------------------------------------
1146//------------------------------------------------------------------------------
1147 template<jit::float_scalar T, bool SAFE_MATH=false>
1148 class efit final : public generic<T, SAFE_MATH> {
1149 private:
1151 const T psimin;
1153 const T dpsi;
1154
1155// Temperature spline coefficients.
1157 const backend::buffer<T> te_c0;
1159 const backend::buffer<T> te_c1;
1161 const backend::buffer<T> te_c2;
1163 const backend::buffer<T> te_c3;
1165 graph::shared_leaf<T, SAFE_MATH> te_scale;
1166
1167// Density spline coefficients.
1169 const backend::buffer<T> ne_c0;
1171 const backend::buffer<T> ne_c1;
1173 const backend::buffer<T> ne_c2;
1175 const backend::buffer<T> ne_c3;
1177 graph::shared_leaf<T, SAFE_MATH> ne_scale;
1178
1179// Pressure spline coefficients.
1181 const backend::buffer<T> pres_c0;
1183 const backend::buffer<T> pres_c1;
1185 const backend::buffer<T> pres_c2;
1187 const backend::buffer<T> pres_c3;
1189 graph::shared_leaf<T, SAFE_MATH> pres_scale;
1190
1192 const T rmin;
1194 const T dr;
1196 const T zmin;
1198 const T dz;
1199
1200// Fpol spline coefficients.
1202 const backend::buffer<T> fpol_c0;
1204 const backend::buffer<T> fpol_c1;
1206 const backend::buffer<T> fpol_c2;
1208 const backend::buffer<T> fpol_c3;
1209
1210// Psi spline coefficients.
1212 const size_t num_cols;
1214 const backend::buffer<T> c00;
1216 const backend::buffer<T> c01;
1218 const backend::buffer<T> c02;
1220 const backend::buffer<T> c03;
1222 const backend::buffer<T> c10;
1224 const backend::buffer<T> c11;
1226 const backend::buffer<T> c12;
1228 const backend::buffer<T> c13;
1230 const backend::buffer<T> c20;
1232 const backend::buffer<T> c21;
1234 const backend::buffer<T> c22;
1236 const backend::buffer<T> c23;
1238 const backend::buffer<T> c30;
1240 const backend::buffer<T> c31;
1242 const backend::buffer<T> c32;
1244 const backend::buffer<T> c33;
1245
1246// Cached values.
1248 graph::shared_leaf<T, SAFE_MATH> x_cache;
1250 graph::shared_leaf<T, SAFE_MATH> y_cache;
1252 graph::shared_leaf<T, SAFE_MATH> z_cache;
1253
1255 graph::shared_leaf<T, SAFE_MATH> ne_cache;
1257 graph::shared_leaf<T, SAFE_MATH> ni_cache;
1259 graph::shared_leaf<T, SAFE_MATH> te_cache;
1261 graph::shared_leaf<T, SAFE_MATH> ti_cache;
1262
1264 graph::shared_vector<T, SAFE_MATH> b_cache;
1265
1267 graph::shared_leaf<T, SAFE_MATH> psi_cache;
1268
1269//------------------------------------------------------------------------------
1279//------------------------------------------------------------------------------
1280 graph::shared_leaf<T, SAFE_MATH>
1281 build_psi(graph::shared_leaf<T, SAFE_MATH> r,
1282 const T r_scale,
1283 const T r_offset,
1284 graph::shared_leaf<T, SAFE_MATH> z,
1285 const T z_scale,
1286 const T z_offset) {
1287 auto c00_temp = graph::piecewise_2D(c00, num_cols, r, r_scale, r_offset, z, z_scale, z_offset);
1288 auto c01_temp = graph::piecewise_2D(c01, num_cols, r, r_scale, r_offset, z, z_scale, z_offset);
1289 auto c02_temp = graph::piecewise_2D(c02, num_cols, r, r_scale, r_offset, z, z_scale, z_offset);
1290 auto c03_temp = graph::piecewise_2D(c03, num_cols, r, r_scale, r_offset, z, z_scale, z_offset);
1291
1292 auto c10_temp = graph::piecewise_2D(c10, num_cols, r, r_scale, r_offset, z, z_scale, z_offset);
1293 auto c11_temp = graph::piecewise_2D(c11, num_cols, r, r_scale, r_offset, z, z_scale, z_offset);
1294 auto c12_temp = graph::piecewise_2D(c12, num_cols, r, r_scale, r_offset, z, z_scale, z_offset);
1295 auto c13_temp = graph::piecewise_2D(c13, num_cols, r, r_scale, r_offset, z, z_scale, z_offset);
1296
1297 auto c20_temp = graph::piecewise_2D(c20, num_cols, r, r_scale, r_offset, z, z_scale, z_offset);
1298 auto c21_temp = graph::piecewise_2D(c21, num_cols, r, r_scale, r_offset, z, z_scale, z_offset);
1299 auto c22_temp = graph::piecewise_2D(c22, num_cols, r, r_scale, r_offset, z, z_scale, z_offset);
1300 auto c23_temp = graph::piecewise_2D(c23, num_cols, r, r_scale, r_offset, z, z_scale, z_offset);
1301
1302 auto c30_temp = graph::piecewise_2D(c30, num_cols, r, r_scale, r_offset, z, z_scale, z_offset);
1303 auto c31_temp = graph::piecewise_2D(c31, num_cols, r, r_scale, r_offset, z, z_scale, z_offset);
1304 auto c32_temp = graph::piecewise_2D(c32, num_cols, r, r_scale, r_offset, z, z_scale, z_offset);
1305 auto c33_temp = graph::piecewise_2D(c33, num_cols, r, r_scale, r_offset, z, z_scale, z_offset);
1306
1307 auto r_norm = (r - r_offset)/r_scale;
1308
1309 auto c0 = build_1D_spline({c00_temp, c01_temp, c02_temp, c03_temp}, z, z_scale, z_offset);
1310 auto c1 = build_1D_spline({c10_temp, c11_temp, c12_temp, c13_temp}, z, z_scale, z_offset);
1311 auto c2 = build_1D_spline({c20_temp, c21_temp, c22_temp, c23_temp}, z, z_scale, z_offset);
1312 auto c3 = build_1D_spline({c30_temp, c31_temp, c32_temp, c33_temp}, z, z_scale, z_offset);
1313
1314 return ((c3*r_norm + c2)*r_norm + c1)*r_norm + c0;
1315 }
1316
1317//------------------------------------------------------------------------------
1325//------------------------------------------------------------------------------
1326 void set_cache(graph::shared_leaf<T, SAFE_MATH> x,
1327 graph::shared_leaf<T, SAFE_MATH> y,
1328 graph::shared_leaf<T, SAFE_MATH> z) {
1329 if (!x->is_match(x_cache) ||
1330 !y->is_match(y_cache) ||
1331 !z->is_match(z_cache)) {
1332 x_cache = x;
1333 y_cache = y;
1334 z_cache = z;
1335
1336 auto r = graph::sqrt(x*x + y*y);
1337
1338 psi_cache = build_psi(r, dr, rmin, z, dz, zmin);
1339
1340 auto n0_temp = graph::piecewise_1D(ne_c0, psi_cache, dpsi, psimin);
1341 auto n1_temp = graph::piecewise_1D(ne_c1, psi_cache, dpsi, psimin);
1342 auto n2_temp = graph::piecewise_1D(ne_c2, psi_cache, dpsi, psimin);
1343 auto n3_temp = graph::piecewise_1D(ne_c3, psi_cache, dpsi, psimin);
1344
1345 ne_cache = ne_scale*build_1D_spline({n0_temp, n1_temp, n2_temp, n3_temp}, psi_cache, dpsi, psimin);
1346
1347 auto t0_temp = graph::piecewise_1D(te_c0, psi_cache, dpsi, psimin);
1348 auto t1_temp = graph::piecewise_1D(te_c1, psi_cache, dpsi, psimin);
1349 auto t2_temp = graph::piecewise_1D(te_c2, psi_cache, dpsi, psimin);
1350 auto t3_temp = graph::piecewise_1D(te_c3, psi_cache, dpsi, psimin);
1351
1352 te_cache = te_scale*build_1D_spline({t0_temp, t1_temp, t2_temp, t3_temp}, psi_cache, dpsi, psimin);
1353
1354 auto p0_temp = graph::piecewise_1D(pres_c0, psi_cache, dpsi, psimin);
1355 auto p1_temp = graph::piecewise_1D(pres_c1, psi_cache, dpsi, psimin);
1356 auto p2_temp = graph::piecewise_1D(pres_c2, psi_cache, dpsi, psimin);
1357 auto p3_temp = graph::piecewise_1D(pres_c3, psi_cache, dpsi, psimin);
1358
1359 auto pressure = pres_scale*build_1D_spline({p0_temp, p1_temp, p2_temp, p3_temp}, psi_cache, dpsi, psimin);
1360
1361 auto q = graph::constant<T, SAFE_MATH> (static_cast<T> (1.60218E-19));
1362
1363 ni_cache = te_cache;
1364 ti_cache = (pressure - ne_cache*te_cache*q)/(ni_cache*q);
1365
1366 auto phi = graph::atan(x, y);
1367
1368 auto br = psi_cache->df(z)/r;
1369
1370 auto b0_temp = graph::piecewise_1D(fpol_c0, psi_cache, dpsi, psimin);
1371 auto b1_temp = graph::piecewise_1D(fpol_c1, psi_cache, dpsi, psimin);
1372 auto b2_temp = graph::piecewise_1D(fpol_c2, psi_cache, dpsi, psimin);
1373 auto b3_temp = graph::piecewise_1D(fpol_c3, psi_cache, dpsi, psimin);
1374
1375 auto bp = build_1D_spline({b0_temp, b1_temp, b2_temp, b3_temp}, psi_cache, dpsi, psimin)/r;
1376
1377 auto bz = -psi_cache->df(r)/r;
1378
1379 auto cos = graph::cos(phi);
1380 auto sin = graph::sin(phi);
1381
1382 b_cache = graph::vector(br*cos - bp*sin,
1383 br*sin + bp*cos,
1384 bz);
1385 }
1386 }
1387
1388 public:
1389//------------------------------------------------------------------------------
1434//------------------------------------------------------------------------------
1435 efit(const T psimin,
1436 const T dpsi,
1437 const backend::buffer<T> te_c0,
1438 const backend::buffer<T> te_c1,
1439 const backend::buffer<T> te_c2,
1440 const backend::buffer<T> te_c3,
1441 graph::shared_leaf<T, SAFE_MATH> te_scale,
1442 const backend::buffer<T> ne_c0,
1443 const backend::buffer<T> ne_c1,
1444 const backend::buffer<T> ne_c2,
1445 const backend::buffer<T> ne_c3,
1446 graph::shared_leaf<T, SAFE_MATH> ne_scale,
1447 const backend::buffer<T> pres_c0,
1448 const backend::buffer<T> pres_c1,
1449 const backend::buffer<T> pres_c2,
1450 const backend::buffer<T> pres_c3,
1451 graph::shared_leaf<T, SAFE_MATH> pres_scale,
1452 const T rmin,
1453 const T dr,
1454 const T zmin,
1455 const T dz,
1456 const backend::buffer<T> fpol_c0,
1457 const backend::buffer<T> fpol_c1,
1458 const backend::buffer<T> fpol_c2,
1459 const backend::buffer<T> fpol_c3,
1460 const size_t num_cols,
1461 const backend::buffer<T> c00,
1462 const backend::buffer<T> c01,
1463 const backend::buffer<T> c02,
1464 const backend::buffer<T> c03,
1465 const backend::buffer<T> c10,
1466 const backend::buffer<T> c11,
1467 const backend::buffer<T> c12,
1468 const backend::buffer<T> c13,
1469 const backend::buffer<T> c20,
1470 const backend::buffer<T> c21,
1471 const backend::buffer<T> c22,
1472 const backend::buffer<T> c23,
1473 const backend::buffer<T> c30,
1474 const backend::buffer<T> c31,
1475 const backend::buffer<T> c32,
1476 const backend::buffer<T> c33) :
1477 generic<T, SAFE_MATH> ({3.34449469E-27} ,{1}),
1478 psimin(psimin), dpsi(dpsi), num_cols(num_cols),
1479 te_c0(te_c0), te_c1(te_c1), te_c2(te_c2), te_c3(te_c3), te_scale(te_scale),
1480 ne_c0(te_c0), ne_c1(te_c1), ne_c2(ne_c2), ne_c3(ne_c3), ne_scale(ne_scale),
1481 pres_c0(pres_c0), pres_c1(pres_c1), pres_c2(pres_c2), pres_c3(pres_c3),
1482 pres_scale(pres_scale), rmin(rmin), dr(dr), zmin(zmin), dz(dz),
1483 fpol_c0(fpol_c0), fpol_c1(fpol_c1), fpol_c2(fpol_c2), fpol_c3(fpol_c3),
1484 c00(c00), c01(c01), c02(c02), c03(c03),
1485 c10(c10), c11(c11), c12(c12), c13(c13),
1486 c20(c20), c21(c21), c22(c22), c23(c23),
1487 c30(c30), c31(c31), c32(c32), c33(c33) {
1488 auto zero = graph::zero<T, SAFE_MATH> ();
1489 x_cache = zero;
1490 y_cache = zero;
1491 z_cache = zero;
1492 }
1493
1494//------------------------------------------------------------------------------
1501//------------------------------------------------------------------------------
1502 virtual graph::shared_leaf<T, SAFE_MATH>
1509
1510//------------------------------------------------------------------------------
1518//------------------------------------------------------------------------------
1519 virtual graph::shared_leaf<T, SAFE_MATH>
1520 get_ion_density(const size_t index,
1521 graph::shared_leaf<T, SAFE_MATH> x,
1522 graph::shared_leaf<T, SAFE_MATH> y,
1523 graph::shared_leaf<T, SAFE_MATH> z) {
1524 set_cache(x, y, z);
1525 return ni_cache;
1526 }
1527
1528//------------------------------------------------------------------------------
1535//------------------------------------------------------------------------------
1536 virtual graph::shared_leaf<T, SAFE_MATH>
1543
1544//------------------------------------------------------------------------------
1552//------------------------------------------------------------------------------
1553 virtual graph::shared_leaf<T, SAFE_MATH>
1554 get_ion_temperature(const size_t index,
1555 graph::shared_leaf<T, SAFE_MATH> x,
1556 graph::shared_leaf<T, SAFE_MATH> y,
1557 graph::shared_leaf<T, SAFE_MATH> z) {
1558 set_cache(x, y, z);
1559 return ti_cache;
1560 }
1561
1562//------------------------------------------------------------------------------
1569//------------------------------------------------------------------------------
1570 virtual graph::shared_vector<T, SAFE_MATH>
1577
1578//------------------------------------------------------------------------------
1585//------------------------------------------------------------------------------
1586 virtual graph::shared_leaf<T, SAFE_MATH>
1587 get_characteristic_field(const size_t device_number=0) final {
1588 auto x_axis = graph::variable<T, SAFE_MATH> (1, "x");
1589 auto y_axis = graph::variable<T, SAFE_MATH> (1, "y");
1590 auto z_axis = graph::variable<T, SAFE_MATH> (1, "z");
1591 x_axis->set(static_cast<T> (1.7));
1592 y_axis->set(static_cast<T> (0.0));
1593 z_axis->set(static_cast<T> (0.0));
1594 auto b_vec = get_magnetic_field(x_axis, y_axis, z_axis);
1595 auto b_mod = b_vec->length();
1596
1597 graph::input_nodes<T, SAFE_MATH> inputs {
1598 graph::variable_cast(x_axis),
1599 graph::variable_cast(y_axis),
1600 graph::variable_cast(z_axis)
1601 };
1602
1603 workflow::manager<T, SAFE_MATH> work(device_number);
1604 solver::newton(work, {
1605 x_axis, z_axis
1606 }, inputs, (psi_cache - psimin)/dpsi, graph::shared_random_state<T, SAFE_MATH> (), static_cast<T> (1.0E-30), 1000, static_cast<T> (0.1));
1607 work.add_item(inputs, {b_mod}, {},
1608 graph::shared_random_state<T, SAFE_MATH> (),
1609 "bmod_at_axis", inputs.back()->size());
1610 work.compile();
1611 work.run();
1612
1613 T result;
1614 work.copy_to_host(b_mod, &result);
1615
1616 return graph::constant<T, SAFE_MATH> (result);
1617 }
1618 };
1619
1620//------------------------------------------------------------------------------
1628//------------------------------------------------------------------------------
1629 template<jit::float_scalar T, bool SAFE_MATH=false>
1630 shared<T, SAFE_MATH> make_efit(const std::string &spline_file) {
1631 int ncid;
1632 sync.lock();
1633 nc_open(spline_file.c_str(), NC_NOWRITE, &ncid);
1634
1635// Load scalar quantities.
1636 int varid;
1637
1638 double rmin_value;
1639 nc_inq_varid(ncid, "rmin", &varid);
1640 nc_get_var(ncid, varid, &rmin_value);
1641
1642 double dr_value;
1643 nc_inq_varid(ncid, "dr", &varid);
1644 nc_get_var(ncid, varid, &dr_value);
1645
1646 double zmin_value;
1647 nc_inq_varid(ncid, "zmin", &varid);
1648 nc_get_var(ncid, varid, &zmin_value);
1649
1650 double dz_value;
1651 nc_inq_varid(ncid, "dz", &varid);
1652 nc_get_var(ncid, varid, &dz_value);
1653
1654 double psimin_value;
1655 nc_inq_varid(ncid, "psimin", &varid);
1656 nc_get_var(ncid, varid, &psimin_value);
1657
1658 double dpsi_value;
1659 nc_inq_varid(ncid, "dpsi", &varid);
1660 nc_get_var(ncid, varid, &dpsi_value);
1661
1662 double pres_scale_value;
1663 nc_inq_varid(ncid, "pres_scale", &varid);
1664 nc_get_var(ncid, varid, &pres_scale_value);
1665
1666 double ne_scale_value;
1667 nc_inq_varid(ncid, "ne_scale", &varid);
1668 nc_get_var(ncid, varid, &ne_scale_value);
1669
1670 double te_scale_value;
1671 nc_inq_varid(ncid, "te_scale", &varid);
1672 nc_get_var(ncid, varid, &te_scale_value);
1673
1674// Load 1D quantities.
1675 int dimid;
1676
1677 size_t numr;
1678 nc_inq_dimid(ncid, "numr", &dimid);
1679 nc_inq_dimlen(ncid, dimid, &numr);
1680
1681 size_t numpsi;
1682 nc_inq_dimid(ncid, "numpsi", &dimid);
1683 nc_inq_dimlen(ncid, dimid, &numpsi);
1684
1685 std::vector<double> fpol_c0_buffer(numpsi);
1686 std::vector<double> fpol_c1_buffer(numpsi);
1687 std::vector<double> fpol_c2_buffer(numpsi);
1688 std::vector<double> fpol_c3_buffer(numpsi);
1689
1690 nc_inq_varid(ncid, "fpol_c0", &varid);
1691 nc_get_var(ncid, varid, fpol_c0_buffer.data());
1692 nc_inq_varid(ncid, "fpol_c1", &varid);
1693 nc_get_var(ncid, varid, fpol_c1_buffer.data());
1694 nc_inq_varid(ncid, "fpol_c2", &varid);
1695 nc_get_var(ncid, varid, fpol_c2_buffer.data());
1696 nc_inq_varid(ncid, "fpol_c3", &varid);
1697 nc_get_var(ncid, varid, fpol_c3_buffer.data());
1698
1699// Load psi grids.
1700 size_t numz;
1701 nc_inq_dimid(ncid, "numz", &dimid);
1702 nc_inq_dimlen(ncid, dimid, &numz);
1703
1704 std::vector<double> psi_c00_buffer(numz*numr);
1705 std::vector<double> psi_c01_buffer(numz*numr);
1706 std::vector<double> psi_c02_buffer(numz*numr);
1707 std::vector<double> psi_c03_buffer(numz*numr);
1708 std::vector<double> psi_c10_buffer(numz*numr);
1709 std::vector<double> psi_c11_buffer(numz*numr);
1710 std::vector<double> psi_c12_buffer(numz*numr);
1711 std::vector<double> psi_c13_buffer(numz*numr);
1712 std::vector<double> psi_c20_buffer(numz*numr);
1713 std::vector<double> psi_c21_buffer(numz*numr);
1714 std::vector<double> psi_c22_buffer(numz*numr);
1715 std::vector<double> psi_c23_buffer(numz*numr);
1716 std::vector<double> psi_c30_buffer(numz*numr);
1717 std::vector<double> psi_c31_buffer(numz*numr);
1718 std::vector<double> psi_c32_buffer(numz*numr);
1719 std::vector<double> psi_c33_buffer(numz*numr);
1720
1721 nc_inq_varid(ncid, "psi_c00", &varid);
1722 nc_get_var(ncid, varid, psi_c00_buffer.data());
1723 nc_inq_varid(ncid, "psi_c01", &varid);
1724 nc_get_var(ncid, varid, psi_c01_buffer.data());
1725 nc_inq_varid(ncid, "psi_c02", &varid);
1726 nc_get_var(ncid, varid, psi_c02_buffer.data());
1727 nc_inq_varid(ncid, "psi_c03", &varid);
1728 nc_get_var(ncid, varid, psi_c03_buffer.data());
1729 nc_inq_varid(ncid, "psi_c10", &varid);
1730 nc_get_var(ncid, varid, psi_c10_buffer.data());
1731 nc_inq_varid(ncid, "psi_c11", &varid);
1732 nc_get_var(ncid, varid, psi_c11_buffer.data());
1733 nc_inq_varid(ncid, "psi_c12", &varid);
1734 nc_get_var(ncid, varid, psi_c12_buffer.data());
1735 nc_inq_varid(ncid, "psi_c13", &varid);
1736 nc_get_var(ncid, varid, psi_c13_buffer.data());
1737 nc_inq_varid(ncid, "psi_c20", &varid);
1738 nc_get_var(ncid, varid, psi_c20_buffer.data());
1739 nc_inq_varid(ncid, "psi_c21", &varid);
1740 nc_get_var(ncid, varid, psi_c21_buffer.data());
1741 nc_inq_varid(ncid, "psi_c22", &varid);
1742 nc_get_var(ncid, varid, psi_c22_buffer.data());
1743 nc_inq_varid(ncid, "psi_c23", &varid);
1744 nc_get_var(ncid, varid, psi_c23_buffer.data());
1745 nc_inq_varid(ncid, "psi_c30", &varid);
1746 nc_get_var(ncid, varid, psi_c30_buffer.data());
1747 nc_inq_varid(ncid, "psi_c31", &varid);
1748 nc_get_var(ncid, varid, psi_c31_buffer.data());
1749 nc_inq_varid(ncid, "psi_c32", &varid);
1750 nc_get_var(ncid, varid, psi_c32_buffer.data());
1751 nc_inq_varid(ncid, "psi_c33", &varid);
1752 nc_get_var(ncid, varid, psi_c33_buffer.data());
1753
1754 std::vector<double> pressure_c0_buffer(numpsi);
1755 std::vector<double> pressure_c1_buffer(numpsi);
1756 std::vector<double> pressure_c2_buffer(numpsi);
1757 std::vector<double> pressure_c3_buffer(numpsi);
1758
1759 nc_inq_varid(ncid, "pressure_c0", &varid);
1760 nc_get_var(ncid, varid, pressure_c0_buffer.data());
1761 nc_inq_varid(ncid, "pressure_c1", &varid);
1762 nc_get_var(ncid, varid, pressure_c1_buffer.data());
1763 nc_inq_varid(ncid, "pressure_c2", &varid);
1764 nc_get_var(ncid, varid, pressure_c2_buffer.data());
1765 nc_inq_varid(ncid, "pressure_c3", &varid);
1766 nc_get_var(ncid, varid, pressure_c3_buffer.data());
1767
1768 std::vector<double> te_c0_buffer(numpsi);
1769 std::vector<double> te_c1_buffer(numpsi);
1770 std::vector<double> te_c2_buffer(numpsi);
1771 std::vector<double> te_c3_buffer(numpsi);
1772
1773 nc_inq_varid(ncid, "te_c0", &varid);
1774 nc_get_var(ncid, varid, te_c0_buffer.data());
1775 nc_inq_varid(ncid, "te_c1", &varid);
1776 nc_get_var(ncid, varid, te_c1_buffer.data());
1777 nc_inq_varid(ncid, "te_c2", &varid);
1778 nc_get_var(ncid, varid, te_c2_buffer.data());
1779 nc_inq_varid(ncid, "te_c3", &varid);
1780 nc_get_var(ncid, varid, te_c3_buffer.data());
1781
1782 std::vector<double> ne_c0_buffer(numpsi);
1783 std::vector<double> ne_c1_buffer(numpsi);
1784 std::vector<double> ne_c2_buffer(numpsi);
1785 std::vector<double> ne_c3_buffer(numpsi);
1786
1787 nc_inq_varid(ncid, "ne_c0", &varid);
1788 nc_get_var(ncid, varid, ne_c0_buffer.data());
1789 nc_inq_varid(ncid, "ne_c1", &varid);
1790 nc_get_var(ncid, varid, ne_c1_buffer.data());
1791 nc_inq_varid(ncid, "ne_c2", &varid);
1792 nc_get_var(ncid, varid, ne_c2_buffer.data());
1793 nc_inq_varid(ncid, "ne_c3", &varid);
1794 nc_get_var(ncid, varid, ne_c3_buffer.data());
1795
1796 nc_close(ncid);
1797 sync.unlock();
1798
1799 auto rmin = static_cast<T> (rmin_value);
1800 auto dr = static_cast<T> (dr_value);
1801 auto zmin = static_cast<T> (zmin_value);
1802 auto dz = static_cast<T> (dz_value);
1803 auto psimin = static_cast<T> (psimin_value);
1804 auto dpsi = static_cast<T> (dpsi_value);
1805 auto pres_scale = graph::constant<T, SAFE_MATH> (static_cast<T> (pres_scale_value));
1806 auto ne_scale = graph::constant<T, SAFE_MATH> (static_cast<T> (ne_scale_value));
1807 auto te_scale = graph::constant<T, SAFE_MATH> (static_cast<T> (te_scale_value));
1808
1809 const auto fpol_c0 = backend::buffer(std::vector<T> (fpol_c0_buffer.begin(), fpol_c0_buffer.end()));
1810 const auto fpol_c1 = backend::buffer(std::vector<T> (fpol_c1_buffer.begin(), fpol_c1_buffer.end()));
1811 const auto fpol_c2 = backend::buffer(std::vector<T> (fpol_c2_buffer.begin(), fpol_c2_buffer.end()));
1812 const auto fpol_c3 = backend::buffer(std::vector<T> (fpol_c3_buffer.begin(), fpol_c3_buffer.end()));
1813
1814 const auto c00 = backend::buffer(std::vector<T> (psi_c00_buffer.begin(), psi_c00_buffer.end()));
1815 const auto c01 = backend::buffer(std::vector<T> (psi_c01_buffer.begin(), psi_c01_buffer.end()));
1816 const auto c02 = backend::buffer(std::vector<T> (psi_c02_buffer.begin(), psi_c02_buffer.end()));
1817 const auto c03 = backend::buffer(std::vector<T> (psi_c03_buffer.begin(), psi_c03_buffer.end()));
1818 const auto c10 = backend::buffer(std::vector<T> (psi_c10_buffer.begin(), psi_c10_buffer.end()));
1819 const auto c11 = backend::buffer(std::vector<T> (psi_c11_buffer.begin(), psi_c11_buffer.end()));
1820 const auto c12 = backend::buffer(std::vector<T> (psi_c12_buffer.begin(), psi_c12_buffer.end()));
1821 const auto c13 = backend::buffer(std::vector<T> (psi_c13_buffer.begin(), psi_c13_buffer.end()));
1822 const auto c20 = backend::buffer(std::vector<T> (psi_c20_buffer.begin(), psi_c20_buffer.end()));
1823 const auto c21 = backend::buffer(std::vector<T> (psi_c21_buffer.begin(), psi_c21_buffer.end()));
1824 const auto c22 = backend::buffer(std::vector<T> (psi_c22_buffer.begin(), psi_c22_buffer.end()));
1825 const auto c23 = backend::buffer(std::vector<T> (psi_c23_buffer.begin(), psi_c23_buffer.end()));
1826 const auto c30 = backend::buffer(std::vector<T> (psi_c30_buffer.begin(), psi_c30_buffer.end()));
1827 const auto c31 = backend::buffer(std::vector<T> (psi_c31_buffer.begin(), psi_c31_buffer.end()));
1828 const auto c32 = backend::buffer(std::vector<T> (psi_c32_buffer.begin(), psi_c32_buffer.end()));
1829 const auto c33 = backend::buffer(std::vector<T> (psi_c33_buffer.begin(), psi_c33_buffer.end()));
1830
1831 const auto pres_c0 = backend::buffer(std::vector<T> (pressure_c0_buffer.begin(), pressure_c0_buffer.end()));
1832 const auto pres_c1 = backend::buffer(std::vector<T> (pressure_c1_buffer.begin(), pressure_c1_buffer.end()));
1833 const auto pres_c2 = backend::buffer(std::vector<T> (pressure_c2_buffer.begin(), pressure_c2_buffer.end()));
1834 const auto pres_c3 = backend::buffer(std::vector<T> (pressure_c3_buffer.begin(), pressure_c3_buffer.end()));
1835
1836 const auto te_c0 = backend::buffer(std::vector<T> (te_c0_buffer.begin(), te_c0_buffer.end()));
1837 const auto te_c1 = backend::buffer(std::vector<T> (te_c1_buffer.begin(), te_c1_buffer.end()));
1838 const auto te_c2 = backend::buffer(std::vector<T> (te_c2_buffer.begin(), te_c2_buffer.end()));
1839 const auto te_c3 = backend::buffer(std::vector<T> (te_c3_buffer.begin(), te_c3_buffer.end()));
1840
1841 const auto ne_c0 = backend::buffer(std::vector<T> (ne_c0_buffer.begin(), ne_c0_buffer.end()));
1842 const auto ne_c1 = backend::buffer(std::vector<T> (ne_c1_buffer.begin(), ne_c1_buffer.end()));
1843 const auto ne_c2 = backend::buffer(std::vector<T> (ne_c2_buffer.begin(), ne_c2_buffer.end()));
1844 const auto ne_c3 = backend::buffer(std::vector<T> (ne_c3_buffer.begin(), ne_c3_buffer.end()));
1845
1846 return std::make_shared<efit<T, SAFE_MATH>> (psimin, dpsi,
1847 te_c0, te_c1, te_c2, te_c3, te_scale,
1848 ne_c0, ne_c1, ne_c2, ne_c3, ne_scale,
1849 pres_c0, pres_c1, pres_c2, pres_c3, pres_scale,
1850 rmin, dr, zmin, dz,
1851 fpol_c0, fpol_c1, fpol_c2, fpol_c3, numz,
1852 c00, c01, c02, c03,
1853 c10, c11, c12, c13,
1854 c20, c21, c22, c23,
1855 c30, c31, c32, c33);
1856 }
1857
1858//******************************************************************************
1859// 3D VMEC equilibrium.
1860//******************************************************************************
1861//------------------------------------------------------------------------------
1868//------------------------------------------------------------------------------
1869 template<jit::float_scalar T, bool SAFE_MATH=false>
1870 class vmec final : public generic<T, SAFE_MATH> {
1871 private:
1873 const T sminh;
1875 const T sminf;
1877 const T ds;
1879 graph::shared_leaf<T, SAFE_MATH> signj;
1880
1881// Poloidal flux coefficients.
1883 const backend::buffer<T> chi_c0;
1885 const backend::buffer<T> chi_c1;
1887 const backend::buffer<T> chi_c2;
1889 const backend::buffer<T> chi_c3;
1890
1891// Toroidal flux coefficients.
1892 graph::shared_leaf<T, SAFE_MATH> dphi;
1893
1894// Radial coefficients.
1896 const std::vector<backend::buffer<T>> rmnc_c0;
1898 const std::vector<backend::buffer<T>> rmnc_c1;
1900 const std::vector<backend::buffer<T>> rmnc_c2;
1902 const std::vector<backend::buffer<T>> rmnc_c3;
1903
1904// Vertical coefficients.
1906 const std::vector<backend::buffer<T>> zmns_c0;
1908 const std::vector<backend::buffer<T>> zmns_c1;
1910 const std::vector<backend::buffer<T>> zmns_c2;
1912 const std::vector<backend::buffer<T>> zmns_c3;
1913
1914// Lambda coefficients.
1916 const std::vector<backend::buffer<T>> lmns_c0;
1918 const std::vector<backend::buffer<T>> lmns_c1;
1920 const std::vector<backend::buffer<T>> lmns_c2;
1922 const std::vector<backend::buffer<T>> lmns_c3;
1923
1925 const backend::buffer<T> xm;
1927 const backend::buffer<T> xn;
1928
1929// Cached values.
1931 graph::shared_leaf<T, SAFE_MATH> s_cache;
1933 graph::shared_leaf<T, SAFE_MATH> u_cache;
1935 graph::shared_leaf<T, SAFE_MATH> v_cache;
1937 graph::shared_leaf<T, SAFE_MATH> x_cache;
1939 graph::shared_leaf<T, SAFE_MATH> y_cache;
1941 graph::shared_leaf<T, SAFE_MATH> z_cache;
1942
1944 graph::shared_vector<T, SAFE_MATH> esups_cache;
1946 graph::shared_vector<T, SAFE_MATH> esupu_cache;
1948 graph::shared_vector<T, SAFE_MATH> esupv_cache;
1949
1951 graph::shared_vector<T, SAFE_MATH> bvec_cache;
1952
1953//------------------------------------------------------------------------------
1959//------------------------------------------------------------------------------
1960 graph::shared_vector<T, SAFE_MATH>
1961 get_esubs(graph::shared_leaf<T, SAFE_MATH> r,
1962 graph::shared_leaf<T, SAFE_MATH> z) {
1963 auto cosv = graph::cos(v_cache);
1964 auto sinv = graph::sin(v_cache);
1965 auto one = graph::one<T, SAFE_MATH> ();
1966 auto zero = graph::zero<T, SAFE_MATH> ();
1967
1968 auto m = graph::matrix(graph::vector(cosv, -sinv, zero),
1969 graph::vector(sinv, cosv, zero),
1970 graph::vector(zero, zero, one ));
1971 return m->dot(graph::vector(r->df(s_cache),
1972 zero,
1973 z->df(s_cache)));
1974 }
1975
1976//------------------------------------------------------------------------------
1982//------------------------------------------------------------------------------
1983 graph::shared_vector<T, SAFE_MATH>
1984 get_esubu(graph::shared_leaf<T, SAFE_MATH> r,
1985 graph::shared_leaf<T, SAFE_MATH> z) {
1986 auto cosv = graph::cos(v_cache);
1987 auto sinv = graph::sin(v_cache);
1988 auto one = graph::one<T, SAFE_MATH> ();
1989 auto zero = graph::zero<T, SAFE_MATH> ();
1990
1991 auto m = graph::matrix(graph::vector(cosv, -sinv, zero),
1992 graph::vector(sinv, cosv, zero),
1993 graph::vector(zero, zero, one ));
1994 return m->dot(graph::vector(r->df(u_cache),
1995 zero,
1996 z->df(u_cache)));
1997 }
1998
1999//------------------------------------------------------------------------------
2005//------------------------------------------------------------------------------
2006 graph::shared_vector<T, SAFE_MATH>
2007 get_esubv(graph::shared_leaf<T, SAFE_MATH> r,
2008 graph::shared_leaf<T, SAFE_MATH> z) {
2009 auto cosv = graph::cos(v_cache);
2010 auto sinv = graph::sin(v_cache);
2011 auto one = graph::one<T, SAFE_MATH> ();
2012 auto zero = graph::zero<T, SAFE_MATH> ();
2013
2014 auto m = graph::matrix(graph::vector(cosv, -sinv, zero),
2015 graph::vector(sinv, cosv, zero),
2016 graph::vector(zero, zero, one ));
2017 return m->dot(graph::vector(r->df(v_cache),
2018 r,
2019 z->df(v_cache)));
2020 }
2021
2022//------------------------------------------------------------------------------
2031//------------------------------------------------------------------------------
2032 graph::shared_leaf<T, SAFE_MATH>
2033 get_jacobian(graph::shared_vector<T, SAFE_MATH> esub_s,
2034 graph::shared_vector<T, SAFE_MATH> esub_u,
2035 graph::shared_vector<T, SAFE_MATH> esub_v) {
2036 return esub_s->dot(esub_u->cross(esub_v));
2037 }
2038
2039//------------------------------------------------------------------------------
2044//------------------------------------------------------------------------------
2045 graph::shared_leaf<T, SAFE_MATH>
2046 get_chi(graph::shared_leaf<T, SAFE_MATH> s) {
2047 auto c0_temp = graph::piecewise_1D(chi_c0, s, ds, sminf);
2048 auto c1_temp = graph::piecewise_1D(chi_c1, s, ds, sminf);
2049 auto c2_temp = graph::piecewise_1D(chi_c2, s, ds, sminf);
2050 auto c3_temp = graph::piecewise_1D(chi_c3, s, ds, sminf);
2051
2052 return build_1D_spline({c0_temp, c1_temp, c2_temp, c3_temp}, s, ds, sminf);
2053 }
2054
2055//------------------------------------------------------------------------------
2060//------------------------------------------------------------------------------
2061 graph::shared_leaf<T, SAFE_MATH>
2062 get_phi(graph::shared_leaf<T, SAFE_MATH> s) {
2063 return signj*dphi*s;
2064 }
2065
2066//------------------------------------------------------------------------------
2074//------------------------------------------------------------------------------
2075 void set_cache(graph::shared_leaf<T, SAFE_MATH> s,
2076 graph::shared_leaf<T, SAFE_MATH> u,
2077 graph::shared_leaf<T, SAFE_MATH> v) {
2078 if (!s->is_match(s_cache) ||
2079 !u->is_match(u_cache) ||
2080 !v->is_match(v_cache)) {
2081 s_cache = s;
2082 u_cache = u;
2083 v_cache = v;
2084
2085 auto s_norm_f = (s - sminf)/ds;
2086
2087 auto zero = graph::zero<T, SAFE_MATH> ();
2088 auto r = zero;
2089 auto z = zero;
2090 auto l = zero;
2091
2092 for (size_t i = 0, ie = xm.size(); i < ie; i++) {
2093 auto rmnc_c0_temp = graph::piecewise_1D(rmnc_c0[i], s, ds, sminf);
2094 auto rmnc_c1_temp = graph::piecewise_1D(rmnc_c1[i], s, ds, sminf);
2095 auto rmnc_c2_temp = graph::piecewise_1D(rmnc_c2[i], s, ds, sminf);
2096 auto rmnc_c3_temp = graph::piecewise_1D(rmnc_c3[i], s, ds, sminf);
2097
2098 auto zmns_c0_temp = graph::piecewise_1D(zmns_c0[i], s, ds, sminf);
2099 auto zmns_c1_temp = graph::piecewise_1D(zmns_c1[i], s, ds, sminf);
2100 auto zmns_c2_temp = graph::piecewise_1D(zmns_c2[i], s, ds, sminf);
2101 auto zmns_c3_temp = graph::piecewise_1D(zmns_c3[i], s, ds, sminf);
2102
2103 auto lmns_c0_temp = graph::piecewise_1D(lmns_c0[i], s, ds, sminh);
2104 auto lmns_c1_temp = graph::piecewise_1D(lmns_c1[i], s, ds, sminh);
2105 auto lmns_c2_temp = graph::piecewise_1D(lmns_c2[i], s, ds, sminh);
2106 auto lmns_c3_temp = graph::piecewise_1D(lmns_c3[i], s, ds, sminh);
2107
2108 auto rmnc = build_1D_spline({rmnc_c0_temp, rmnc_c1_temp, rmnc_c2_temp, rmnc_c3_temp},
2109 s, ds, sminf);
2110 auto zmns = build_1D_spline({zmns_c0_temp, zmns_c1_temp, zmns_c2_temp, zmns_c3_temp},
2111 s, ds, sminf);
2112 auto lmns = build_1D_spline({lmns_c0_temp, lmns_c1_temp, lmns_c2_temp, lmns_c3_temp},
2113 s, ds, sminh);
2114
2115 auto m = graph::constant<T, SAFE_MATH> (xm[i]);
2116 auto n = graph::constant<T, SAFE_MATH> (xn[i]);
2117
2118 auto sinmn = graph::sin(m*u - n*v);
2119
2120 r = r + rmnc*graph::cos(m*u - n*v);
2121 z = z + zmns*sinmn;
2122 l = l + lmns*sinmn;
2123 }
2124
2125 x_cache = r*graph::cos(v);
2126 y_cache = r*graph::sin(v);
2127 z_cache = z;
2128
2129 auto esubs = get_esubs(r, z);
2130 auto esubu = get_esubu(r, z);
2131 auto esubv = get_esubv(r, z);
2132
2133 auto jacobian = get_jacobian(esubs, esubu, esubv);
2134
2135 esups_cache = esubu->cross(esubv)/jacobian;
2136 esupu_cache = esubv->cross(esubs)/jacobian;
2137 esupv_cache = esubs->cross(esubu)/jacobian;
2138
2139 auto phip = get_phi(s)->df(s);
2140 auto jbsupu = get_chi(s_norm_f)->df(s) - phip*l->df(v);
2141 auto jbsupv = phip*(1.0 + l->df(u));
2142 bvec_cache = (jbsupu*esubu + jbsupv*esubv)/jacobian;
2143 }
2144 }
2145
2146//------------------------------------------------------------------------------
2151//------------------------------------------------------------------------------
2152 graph::shared_leaf<T, SAFE_MATH>
2153 get_profile(graph::shared_leaf<T, SAFE_MATH> s) {
2154 return graph::pow((1.0 - graph::pow(graph::sqrt(s*s), 1.5)), 2.0);
2155 }
2156
2157 public:
2158//------------------------------------------------------------------------------
2184//------------------------------------------------------------------------------
2185 vmec(const T sminh,
2186 const T sminf,
2187 const T ds,
2188 graph::shared_leaf<T, SAFE_MATH> dphi,
2189 graph::shared_leaf<T, SAFE_MATH> signj,
2190 const backend::buffer<T> chi_c0,
2191 const backend::buffer<T> chi_c1,
2192 const backend::buffer<T> chi_c2,
2193 const backend::buffer<T> chi_c3,
2194 const std::vector<backend::buffer<T>> rmnc_c0,
2195 const std::vector<backend::buffer<T>> rmnc_c1,
2196 const std::vector<backend::buffer<T>> rmnc_c2,
2197 const std::vector<backend::buffer<T>> rmnc_c3,
2198 const std::vector<backend::buffer<T>> zmns_c0,
2199 const std::vector<backend::buffer<T>> zmns_c1,
2200 const std::vector<backend::buffer<T>> zmns_c2,
2201 const std::vector<backend::buffer<T>> zmns_c3,
2202 const std::vector<backend::buffer<T>> lmns_c0,
2203 const std::vector<backend::buffer<T>> lmns_c1,
2204 const std::vector<backend::buffer<T>> lmns_c2,
2205 const std::vector<backend::buffer<T>> lmns_c3,
2206 const backend::buffer<T> xm,
2207 const backend::buffer<T> xn) :
2208 generic<T, SAFE_MATH> ({3.34449469E-27} ,{1}),
2209 sminh(sminh), sminf(sminf), ds(ds), dphi(dphi), signj(signj),
2210 chi_c0(chi_c0), chi_c1(chi_c1), chi_c2(chi_c2), chi_c3(chi_c3),
2211 rmnc_c0(rmnc_c0), rmnc_c1(rmnc_c1), rmnc_c2(rmnc_c2), rmnc_c3(rmnc_c3),
2212 zmns_c0(zmns_c0), zmns_c1(zmns_c1), zmns_c2(zmns_c2), zmns_c3(zmns_c3),
2213 lmns_c0(lmns_c0), lmns_c1(lmns_c1), lmns_c2(lmns_c2), lmns_c3(lmns_c3),
2214 xm(xm), xn(xn) {
2215 auto zero = graph::zero<T, SAFE_MATH> ();
2216 s_cache = zero;
2217 u_cache = zero;
2218 v_cache = zero;
2219 }
2220
2221//------------------------------------------------------------------------------
2228//------------------------------------------------------------------------------
2229 virtual graph::shared_vector<T, SAFE_MATH>
2230 get_esup1(graph::shared_leaf<T, SAFE_MATH> s,
2231 graph::shared_leaf<T, SAFE_MATH> u,
2232 graph::shared_leaf<T, SAFE_MATH> v) {
2233 set_cache(s, u, v);
2234 return esups_cache;
2235 }
2236
2237//------------------------------------------------------------------------------
2244//------------------------------------------------------------------------------
2245 virtual graph::shared_vector<T, SAFE_MATH>
2246 get_esup2(graph::shared_leaf<T, SAFE_MATH> s,
2247 graph::shared_leaf<T, SAFE_MATH> u,
2248 graph::shared_leaf<T, SAFE_MATH> v) {
2249 set_cache(s, u, v);
2250 return esupu_cache;
2251 }
2252
2253//------------------------------------------------------------------------------
2260//------------------------------------------------------------------------------
2261 virtual graph::shared_vector<T, SAFE_MATH>
2262 get_esup3(graph::shared_leaf<T, SAFE_MATH> s,
2263 graph::shared_leaf<T, SAFE_MATH> u,
2264 graph::shared_leaf<T, SAFE_MATH> v) {
2265 set_cache(s, u, v);
2266 return esupv_cache;
2267 }
2268
2269//------------------------------------------------------------------------------
2276//------------------------------------------------------------------------------
2277 virtual graph::shared_leaf<T, SAFE_MATH>
2278 get_electron_density(graph::shared_leaf<T, SAFE_MATH> s,
2279 graph::shared_leaf<T, SAFE_MATH> u,
2280 graph::shared_leaf<T, SAFE_MATH> v) {
2281 auto ne_scale = graph::constant<T, SAFE_MATH> (static_cast<T> (1.0E19));
2282 return ne_scale*get_profile(s);
2283 }
2284
2285//------------------------------------------------------------------------------
2293//------------------------------------------------------------------------------
2294 virtual graph::shared_leaf<T, SAFE_MATH>
2301
2302//------------------------------------------------------------------------------
2309//------------------------------------------------------------------------------
2310 virtual graph::shared_leaf<T, SAFE_MATH>
2311 get_electron_temperature(graph::shared_leaf<T, SAFE_MATH> s,
2312 graph::shared_leaf<T, SAFE_MATH> u,
2313 graph::shared_leaf<T, SAFE_MATH> v) {
2314 auto te_scale = graph::constant<T, SAFE_MATH> (static_cast<T> (1000.0));
2315 return te_scale*get_profile(s);
2316 }
2317
2318//------------------------------------------------------------------------------
2326//------------------------------------------------------------------------------
2327 virtual graph::shared_leaf<T, SAFE_MATH>
2334
2335//------------------------------------------------------------------------------
2342//------------------------------------------------------------------------------
2343 virtual graph::shared_vector<T, SAFE_MATH>
2344 get_magnetic_field(graph::shared_leaf<T, SAFE_MATH> s,
2345 graph::shared_leaf<T, SAFE_MATH> u,
2346 graph::shared_leaf<T, SAFE_MATH> v) {
2347 set_cache(s, u, v);
2348 return bvec_cache;
2349 }
2350
2351//------------------------------------------------------------------------------
2358//------------------------------------------------------------------------------
2359 virtual graph::shared_leaf<T, SAFE_MATH>
2360 get_characteristic_field(const size_t device_number=0) final {
2361 auto s_axis = graph::zero<T, SAFE_MATH> ();
2362 auto u_axis = graph::zero<T, SAFE_MATH> ();
2363 auto v_axis = graph::zero<T, SAFE_MATH> ();
2364 auto b_vec = get_magnetic_field(s_axis, u_axis, v_axis);
2365 return b_vec->length();
2366 }
2367
2368//------------------------------------------------------------------------------
2375//------------------------------------------------------------------------------
2376 virtual graph::shared_leaf<T, SAFE_MATH>
2377 get_x(graph::shared_leaf<T, SAFE_MATH> s,
2378 graph::shared_leaf<T, SAFE_MATH> u,
2379 graph::shared_leaf<T, SAFE_MATH> v) {
2380 set_cache(s, u, v);
2381 return x_cache;
2382 }
2383
2384//------------------------------------------------------------------------------
2391//------------------------------------------------------------------------------
2392 virtual graph::shared_leaf<T, SAFE_MATH>
2393 get_y(graph::shared_leaf<T, SAFE_MATH> s,
2394 graph::shared_leaf<T, SAFE_MATH> u,
2395 graph::shared_leaf<T, SAFE_MATH> v) {
2396 set_cache(s, u, v);
2397 return y_cache;
2398 }
2399
2400//------------------------------------------------------------------------------
2407//------------------------------------------------------------------------------
2408 virtual graph::shared_leaf<T, SAFE_MATH>
2409 get_z(graph::shared_leaf<T, SAFE_MATH> s,
2410 graph::shared_leaf<T, SAFE_MATH> u,
2411 graph::shared_leaf<T, SAFE_MATH> v) {
2412 set_cache(s, u, v);
2413 return z_cache;
2414 }
2415 };
2416
2417//------------------------------------------------------------------------------
2425//------------------------------------------------------------------------------
2426 template<jit::float_scalar T, bool SAFE_MATH=false>
2427 shared<T, SAFE_MATH> make_vmec(const std::string &spline_file) {
2428 int ncid;
2429 sync.lock();
2430 nc_open(spline_file.c_str(), NC_NOWRITE, &ncid);
2431
2432// Load scalar quantities.
2433 int varid;
2434
2435 double sminf_value;
2436 nc_inq_varid(ncid, "sminf", &varid);
2437 nc_get_var(ncid, varid, &sminf_value);
2438
2439 double sminh_value;
2440 nc_inq_varid(ncid, "sminh", &varid);
2441 nc_get_var(ncid, varid, &sminh_value);
2442
2443 double ds_value;
2444 nc_inq_varid(ncid, "ds", &varid);
2445 nc_get_var(ncid, varid, &ds_value);
2446
2447 double dphi_value;
2448 nc_inq_varid(ncid, "dphi", &varid);
2449 nc_get_var(ncid, varid, &dphi_value);
2450
2451 double signj_value;
2452 nc_inq_varid(ncid, "signj", &varid);
2453 nc_get_var(ncid, varid, &signj_value);
2454
2455// Load 1D quantities.
2456 int dimid;
2457
2458 size_t numsf;
2459 nc_inq_dimid(ncid, "numsf", &dimid);
2460 nc_inq_dimlen(ncid, dimid, &numsf);
2461
2462 std::vector<double> chi_c0_buffer(numsf);
2463 std::vector<double> chi_c1_buffer(numsf);
2464 std::vector<double> chi_c2_buffer(numsf);
2465 std::vector<double> chi_c3_buffer(numsf);
2466
2467 nc_inq_varid(ncid, "chi_c0", &varid);
2468 nc_get_var(ncid, varid, chi_c0_buffer.data());
2469 nc_inq_varid(ncid, "chi_c1", &varid);
2470 nc_get_var(ncid, varid, chi_c1_buffer.data());
2471 nc_inq_varid(ncid, "chi_c2", &varid);
2472 nc_get_var(ncid, varid, chi_c2_buffer.data());
2473 nc_inq_varid(ncid, "chi_c3", &varid);
2474 nc_get_var(ncid, varid, chi_c3_buffer.data());
2475
2476// Load 2D quantities.
2477 size_t numsh;
2478 nc_inq_dimid(ncid, "numsh", &dimid);
2479 nc_inq_dimlen(ncid, dimid, &numsh);
2480
2481 size_t nummn;
2482 nc_inq_dimid(ncid, "nummn", &dimid);
2483 nc_inq_dimlen(ncid, dimid, &nummn);
2484
2485 std::vector<std::vector<double>> rmnc_c0_buffer(nummn, std::vector<double> (numsf));
2486 std::vector<std::vector<double>> rmnc_c1_buffer(nummn, std::vector<double> (numsf));
2487 std::vector<std::vector<double>> rmnc_c2_buffer(nummn, std::vector<double> (numsf));
2488 std::vector<std::vector<double>> rmnc_c3_buffer(nummn, std::vector<double> (numsf));
2489
2490 nc_inq_varid(ncid, "rmnc_c0", &varid);
2491 for (size_t i = 0; i < nummn; i++) {
2492 const array<size_t, 2> start = {i, 0};
2493 const array<size_t, 2> count = {1, numsf};
2494 nc_get_vara(ncid, varid, start.data(), count.data(),
2495 rmnc_c0_buffer[i].data());
2496 }
2497 nc_inq_varid(ncid, "rmnc_c1", &varid);
2498 for (size_t i = 0; i < nummn; i++) {
2499 const array<size_t, 2> start = {i, 0};
2500 const array<size_t, 2> count = {1, numsf};
2501 nc_get_vara(ncid, varid, start.data(), count.data(),
2502 rmnc_c1_buffer[i].data());
2503 }
2504 nc_inq_varid(ncid, "rmnc_c2", &varid);
2505 for (size_t i = 0; i < nummn; i++) {
2506 const array<size_t, 2> start = {i, 0};
2507 const array<size_t, 2> count = {1, numsf};
2508 nc_get_vara(ncid, varid, start.data(), count.data(),
2509 rmnc_c2_buffer[i].data());
2510 }
2511 nc_inq_varid(ncid, "rmnc_c3", &varid);
2512 for (size_t i = 0; i < nummn; i++) {
2513 const array<size_t, 2> start = {i, 0};
2514 const array<size_t, 2> count = {1, numsf};
2515 nc_get_vara(ncid, varid, start.data(), count.data(),
2516 rmnc_c3_buffer[i].data());
2517 }
2518
2519 std::vector<std::vector<double>> zmns_c0_buffer(nummn, std::vector<double> (numsf));
2520 std::vector<std::vector<double>> zmns_c1_buffer(nummn, std::vector<double> (numsf));
2521 std::vector<std::vector<double>> zmns_c2_buffer(nummn, std::vector<double> (numsf));
2522 std::vector<std::vector<double>> zmns_c3_buffer(nummn, std::vector<double> (numsf));
2523
2524 nc_inq_varid(ncid, "zmns_c0", &varid);
2525 for (size_t i = 0; i < nummn; i++) {
2526 const array<size_t, 2> start = {i, 0};
2527 const array<size_t, 2> count = {1, numsf};
2528 nc_get_vara(ncid, varid, start.data(), count.data(),
2529 zmns_c0_buffer[i].data());
2530 }
2531 nc_inq_varid(ncid, "zmns_c1", &varid);
2532 for (size_t i = 0; i < nummn; i++) {
2533 const array<size_t, 2> start = {i, 0};
2534 const array<size_t, 2> count = {1, numsf};
2535 nc_get_vara(ncid, varid, start.data(), count.data(),
2536 zmns_c1_buffer[i].data());
2537 }
2538 nc_inq_varid(ncid, "zmns_c2", &varid);
2539 for (size_t i = 0; i < nummn; i++) {
2540 const array<size_t, 2> start = {i, 0};
2541 const array<size_t, 2> count = {1, numsf};
2542 nc_get_vara(ncid, varid, start.data(), count.data(),
2543 zmns_c2_buffer[i].data());
2544 }
2545 nc_inq_varid(ncid, "zmns_c3", &varid);
2546 for (size_t i = 0; i < nummn; i++) {
2547 const array<size_t, 2> start = {i, 0};
2548 const array<size_t, 2> count = {1, numsf};
2549 nc_get_vara(ncid, varid, start.data(), count.data(),
2550 zmns_c3_buffer[i].data());
2551 }
2552
2553 std::vector<std::vector<double>> lmns_c0_buffer(nummn, std::vector<double> (numsh));
2554 std::vector<std::vector<double>> lmns_c1_buffer(nummn, std::vector<double> (numsh));
2555 std::vector<std::vector<double>> lmns_c2_buffer(nummn, std::vector<double> (numsh));
2556 std::vector<std::vector<double>> lmns_c3_buffer(nummn, std::vector<double> (numsh));
2557
2558 nc_inq_varid(ncid, "lmns_c0", &varid);
2559 for (size_t i = 0; i < nummn; i++) {
2560 const array<size_t, 2> start = {i, 0};
2561 const array<size_t, 2> count = {1, numsh};
2562 nc_get_vara(ncid, varid, start.data(), count.data(),
2563 lmns_c0_buffer[i].data());
2564 }
2565 nc_inq_varid(ncid, "lmns_c1", &varid);
2566 for (size_t i = 0; i < nummn; i++) {
2567 const array<size_t, 2> start = {i, 0};
2568 const array<size_t, 2> count = {1, numsh};
2569 nc_get_vara(ncid, varid, start.data(), count.data(),
2570 lmns_c1_buffer[i].data());
2571 }
2572 nc_inq_varid(ncid, "lmns_c2", &varid);
2573 for (size_t i = 0; i < nummn; i++) {
2574 const array<size_t, 2> start = {i, 0};
2575 const array<size_t, 2> count = {1, numsh};
2576 nc_get_vara(ncid, varid, start.data(), count.data(),
2577 lmns_c2_buffer[i].data());
2578 }
2579 nc_inq_varid(ncid, "lmns_c3", &varid);
2580 for (size_t i = 0; i < nummn; i++) {
2581 const array<size_t, 2> start = {i, 0};
2582 const array<size_t, 2> count = {1, numsh};
2583 nc_get_vara(ncid, varid, start.data(), count.data(),
2584 lmns_c3_buffer[i].data());
2585 }
2586
2587 std::vector<double> xm_buffer(nummn);
2588 nc_inq_varid(ncid, "xm", &varid);
2589 nc_get_var(ncid, varid, xm_buffer.data());
2590
2591 std::vector<double> xn_buffer(nummn);
2592 nc_inq_varid(ncid, "xn", &varid);
2593 nc_get_var(ncid, varid, xn_buffer.data());
2594
2595 nc_close(ncid);
2596 sync.unlock();
2597
2598 auto sminf = static_cast<T> (sminf_value);
2599 auto sminh = static_cast<T> (sminh_value);
2600 auto ds = static_cast<T> (ds_value);
2601 auto dphi = graph::constant<T, SAFE_MATH> (static_cast<T> (dphi_value));
2602 auto signj = graph::constant<T, SAFE_MATH> (static_cast<T> (signj_value));
2603
2604 const backend::buffer<T> chi_c0(std::vector<T> (chi_c0_buffer.begin(), chi_c0_buffer.end()));
2605 const backend::buffer<T> chi_c1(std::vector<T> (chi_c1_buffer.begin(), chi_c1_buffer.end()));
2606 const backend::buffer<T> chi_c2(std::vector<T> (chi_c2_buffer.begin(), chi_c2_buffer.end()));
2607 const backend::buffer<T> chi_c3(std::vector<T> (chi_c3_buffer.begin(), chi_c3_buffer.end()));
2608
2609 std::vector<backend::buffer<T>> rmnc_c0(nummn);
2610 std::vector<backend::buffer<T>> rmnc_c1(nummn);
2611 std::vector<backend::buffer<T>> rmnc_c2(nummn);
2612 std::vector<backend::buffer<T>> rmnc_c3(nummn);
2613
2614 std::vector<backend::buffer<T>> zmns_c0(nummn);
2615 std::vector<backend::buffer<T>> zmns_c1(nummn);
2616 std::vector<backend::buffer<T>> zmns_c2(nummn);
2617 std::vector<backend::buffer<T>> zmns_c3(nummn);
2618
2619 std::vector<backend::buffer<T>> lmns_c0(nummn);
2620 std::vector<backend::buffer<T>> lmns_c1(nummn);
2621 std::vector<backend::buffer<T>> lmns_c2(nummn);
2622 std::vector<backend::buffer<T>> lmns_c3(nummn);
2623
2624 for (size_t i = 0; i < nummn; i++) {
2625 rmnc_c0[i] = backend::buffer(std::vector<T> (rmnc_c0_buffer[i].begin(), rmnc_c0_buffer[i].end()));
2626 rmnc_c1[i] = backend::buffer(std::vector<T> (rmnc_c1_buffer[i].begin(), rmnc_c1_buffer[i].end()));
2627 rmnc_c2[i] = backend::buffer(std::vector<T> (rmnc_c2_buffer[i].begin(), rmnc_c2_buffer[i].end()));
2628 rmnc_c3[i] = backend::buffer(std::vector<T> (rmnc_c3_buffer[i].begin(), rmnc_c3_buffer[i].end()));
2629
2630 zmns_c0[i] = backend::buffer(std::vector<T> (zmns_c0_buffer[i].begin(), zmns_c0_buffer[i].end()));
2631 zmns_c1[i] = backend::buffer(std::vector<T> (zmns_c1_buffer[i].begin(), zmns_c1_buffer[i].end()));
2632 zmns_c2[i] = backend::buffer(std::vector<T> (zmns_c2_buffer[i].begin(), zmns_c2_buffer[i].end()));
2633 zmns_c3[i] = backend::buffer(std::vector<T> (zmns_c3_buffer[i].begin(), zmns_c3_buffer[i].end()));
2634
2635 lmns_c0[i] = backend::buffer(std::vector<T> (lmns_c0_buffer[i].begin(), lmns_c0_buffer[i].end()));
2636 lmns_c1[i] = backend::buffer(std::vector<T> (lmns_c1_buffer[i].begin(), lmns_c1_buffer[i].end()));
2637 lmns_c2[i] = backend::buffer(std::vector<T> (lmns_c2_buffer[i].begin(), lmns_c2_buffer[i].end()));
2638 lmns_c3[i] = backend::buffer(std::vector<T> (lmns_c3_buffer[i].begin(), lmns_c3_buffer[i].end()));
2639 }
2640
2641 const backend::buffer<T> xm(std::vector<T> (xm_buffer.begin(), xm_buffer.end()));
2642 const backend::buffer<T> xn(std::vector<T> (xn_buffer.begin(), xn_buffer.end()));
2643
2644 return std::make_shared<vmec<T, SAFE_MATH>> (sminh, sminf, ds, dphi, signj,
2645 chi_c0, chi_c1, chi_c2, chi_c3,
2646 rmnc_c0, rmnc_c1, rmnc_c2, rmnc_c3,
2647 zmns_c0, zmns_c1, zmns_c2, zmns_c3,
2648 lmns_c0, lmns_c1, lmns_c2, lmns_c3,
2649 xm, xn);
2650 }
2651}
2652
2653#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:1148
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:1520
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:1571
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:1537
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:1435
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:1503
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:1554
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:1587
equilibrium::generic
Class representing a generic equilibrium.
Definition equilibrium.hpp:238
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:399
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:433
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:448
equilibrium::generic::ion_charges
const std::vector< uint8_t > ion_charges
Ion charge for each species.
Definition equilibrium.hpp:243
equilibrium::generic::~generic
virtual ~generic()
Destructor.
Definition equilibrium.hpp:262
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:416
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:463
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:289
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:382
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:269
equilibrium::generic::get_ion_mass
T get_ion_mass(const size_t index) const
Get the mass for an ion species.
Definition equilibrium.hpp:279
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:241
equilibrium::guassian_density
Guassian density with uniform magnetic field equilibrium.
Definition equilibrium.hpp:993
equilibrium::guassian_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:1010
equilibrium::guassian_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:1092
equilibrium::guassian_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:1076
equilibrium::guassian_density::guassian_density
guassian_density()
Construct a guassian density with uniform magnetic field.
Definition equilibrium.hpp:998
equilibrium::guassian_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:1044
equilibrium::guassian_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:1027
equilibrium::guassian_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:1060
equilibrium::no_magnetic_field
Uniform density with no magnetic field equilibrium.
Definition equilibrium.hpp:484
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:569
equilibrium::no_magnetic_field::no_magnetic_field
no_magnetic_field()
Construct a linear density with no magnetic field.
Definition equilibrium.hpp:489
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:553
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:519
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:537
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:501
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:585
equilibrium::slab_density
Vary density with uniform magnetic field equilibrium.
Definition equilibrium.hpp:737
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:806
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:822
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:838
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:790
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:772
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:754
equilibrium::slab_density::slab_density
slab_density()
Construct a guassian density with uniform magnetic field.
Definition equilibrium.hpp:742
equilibrium::slab_field
Vary density with uniform magnetic field equilibrium.
Definition equilibrium.hpp:866
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:883
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:901
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:917
equilibrium::slab_field::slab_field
slab_field()
Construct a guassian density with uniform magnetic field.
Definition equilibrium.hpp:871
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:966
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:951
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:935
equilibrium::slab
Uniform density with varying magnetic field equilibrium.
Definition equilibrium.hpp:613
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:709
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:662
equilibrium::slab::slab
slab()
Construct a guassian density with uniform magnetic field.
Definition equilibrium.hpp:618
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:694
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:646
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:630
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:678
equilibrium::vmec
3D VMEC equilibrium.
Definition equilibrium.hpp:1870
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:2278
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:2262
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:2377
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:2393
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:2360
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:2246
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:2311
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:2230
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:2328
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:2295
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:2344
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:2409
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:2185
workflow::manager
Class representing a workflow manager.
Definition workflow.hpp:171
workflow::manager::run
void run()
Run work items.
Definition workflow.hpp:285
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:315
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:221
workflow::manager::compile
void compile()
Compile the workflow items.
Definition workflow.hpp:262
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:224
equilibrium::make_slab_field
shared< T, SAFE_MATH > make_slab_field()
Convenience function to build a slab density equilibrium.
Definition equilibrium.hpp:980
equilibrium::make_slab_density
shared< T, SAFE_MATH > make_slab_density()
Convenience function to build a slab density equilibrium.
Definition equilibrium.hpp:852
equilibrium::make_vmec
shared< T, SAFE_MATH > make_vmec(const std::string &spline_file)
Convenience function to build an VMEC equilibrium.
Definition equilibrium.hpp:2427
equilibrium::shared
std::shared_ptr< generic< T, SAFE_MATH > > shared
Convenience type alias for shared equilibria.
Definition equilibrium.hpp:472
equilibrium::make_guassian_density
shared< T, SAFE_MATH > make_guassian_density()
Convenience function to build a guassian density equilibrium.
Definition equilibrium.hpp:1106
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:599
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:1123
equilibrium::make_slab
shared< T, SAFE_MATH > make_slab()
Convenience function to build a slab equilibrium.
Definition equilibrium.hpp:723
equilibrium::make_efit
shared< T, SAFE_MATH > make_efit(const std::string &spline_file)
Convenience function to build an EFIT equilibrium.
Definition equilibrium.hpp:1630
graph::zero
constexpr shared_leaf< T, SAFE_MATH > zero()
Forward declare for zero.
Definition node.hpp:994
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 convience function.
Definition piecewise.hpp:563
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:802
graph::shared_vector
std::shared_ptr< vector_quantity< T, SAFE_MATH > > shared_vector
Convenience type for shared vector quantities.
Definition vector.hpp:119
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 convience function.
Definition piecewise.hpp:1281
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:133
graph::exp
shared_leaf< T, SAFE_MATH > exp(shared_leaf< T, SAFE_MATH > x)
Define exp convience 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:272
graph::input_nodes
std::vector< shared_variable< T, SAFE_MATH > > input_nodes
Convenience type alias for a vector of inputs.
Definition node.hpp:1730
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:1746
graph::sin
shared_leaf< T, SAFE_MATH > sin(shared_leaf< T, SAFE_MATH > x)
Define sine convience function.
Definition trigonometry.hpp:229
graph::sqrt
shared_leaf< T, SAFE_MATH > sqrt(shared_leaf< T, SAFE_MATH > x)
Define sqrt convience 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:413
graph::shared_leaf
std::shared_ptr< leaf_node< T, SAFE_MATH > > shared_leaf
Convenience type alias for shared leaf nodes.
Definition node.hpp:673
graph::cos
shared_leaf< T, SAFE_MATH > cos(shared_leaf< T, SAFE_MATH > x)
Define cosine convience function.
Definition trigonometry.hpp:481
graph::output_nodes
std::vector< shared_leaf< T, SAFE_MATH > > output_nodes
Convenience type alias for a vector of output nodes.
Definition node.hpp:688
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 tolarance=1.0E-30, const size_t max_iterations=1000, const T step=1.0)
Determine the value of vars to minimze the loss function.
Definition newton.hpp:34
newton.hpp
Sets up the kernel for a newtons method.
piecewise.hpp
Piecewise nodes.
trigonometry.hpp
Trigonometry functions.
vector.hpp
Defines vectors of graphs.