Extension: Geometries
#include "dg/geometries/geometries.h"
Loading...
Searching...
No Matches
magnetic_field.h
Go to the documentation of this file.
1#pragma once
2
3#include <map>
4#include "fluxfunctions.h"
5
10namespace dg
11{
12namespace geo
13{
23
25enum class equilibrium
26{
27 solovev,
28 taylor,
30 guenter,
31 toroidal,
33};
35enum class modifier
36{
37 none,
38 heaviside,
39 sol_pfr,
41 // TODO There should be the "circular" parameter from feltor and should there be a "heavisideX"?
42};
49enum class description
50{
51 standardO,
52 standardX,
53 doubleX,
54 none,
55 square,
57};
59inline const std::map<std::string, equilibrium> str2equilibrium{
60 {"solovev", equilibrium::solovev},
61 {"taylor", equilibrium::taylor},
62 {"polynomial", equilibrium::polynomial},
63 {"guenter", equilibrium::guenter},
64 {"toroidal", equilibrium::toroidal},
65 {"circular", equilibrium::circular}
66};
67inline const std::map<std::string, modifier> str2modifier{
68 {"none", modifier::none},
69 {"heaviside", modifier::heaviside},
70 {"sol_pfr", modifier::sol_pfr},
71 {"sol_pfr_2X", modifier::sol_pfr_2X}
72};
73inline const std::map<std::string, description> str2description{
74 {"standardO", description::standardO},
75 {"standardX", description::standardX},
76 {"doubleX", description::doubleX},
77 {"square", description::square},
78 {"none", description::none},
79 {"centeredX", description::centeredX}
80};
82
83//Meta-data about magnetic fields
84//
94{
99 m_a = 1, m_elongation = 1, m_triangularity = 0;
100 m_equilibrium = equilibrium::toroidal;
101 m_modifier = modifier::none;
102 m_description = description::none;
103 }
115 equilibrium equ, modifier mod, description des): m_a(a),
116 m_elongation(elongation),
117 m_triangularity( triangularity),
118 m_equilibrium( equ),
119 m_modifier(mod), m_description( des){}
125 double a() const{return m_a;}
131 double elongation() const{return m_elongation;}
137 double triangularity() const{return m_triangularity;}
139 equilibrium getEquilibrium() const{return m_equilibrium;}
141 modifier getModifier() const{return m_modifier;}
143 description getDescription() const{return m_description;}
144 private:
145 double m_a,
146 m_elongation,
147 m_triangularity;
148 equilibrium m_equilibrium;
149 modifier m_modifier;
150 description m_description;
151};
152
165{
170 ): m_R0(R0), m_psip(psip), m_ipol(ipol), m_params(gp){}
171 void set( double R0, const CylindricalFunctorsLvl2& psip, const
173 {
174 m_R0=R0;
175 m_psip=psip;
176 m_ipol=ipol;
177 m_params = gp;
178 }
180 double R0()const {return m_R0;}
182 const CylindricalFunctor& psip()const{return m_psip.f();}
184 const CylindricalFunctor& psipR()const{return m_psip.dfx();}
186 const CylindricalFunctor& psipZ()const{return m_psip.dfy();}
188 const CylindricalFunctor& psipRR()const{return m_psip.dfxx();}
190 const CylindricalFunctor& psipRZ()const{return m_psip.dfxy();}
192 const CylindricalFunctor& psipZZ()const{return m_psip.dfyy();}
194 const CylindricalFunctor& ipol()const{return m_ipol.f();}
196 const CylindricalFunctor& ipolR()const{return m_ipol.dfx();}
198 const CylindricalFunctor& ipolZ()const{return m_ipol.dfy();}
199
200 const CylindricalFunctorsLvl2& get_psip() const{return m_psip;}
201 const CylindricalFunctorsLvl1& get_ipol() const{return m_ipol;}
207 const MagneticFieldParameters& params() const{return m_params;}
208
209 private:
210 double m_R0;
214};
215
217inline CylindricalFunctorsLvl1 periodify( const CylindricalFunctorsLvl1& in, double R0, double R1, double Z0, double Z1, bc bcx, bc bcy)
218{
219 return CylindricalFunctorsLvl1(
220 Periodify( in.f(), R0, R1, Z0, Z1, bcx, bcy),
221 Periodify( in.dfx(), R0, R1, Z0, Z1, inverse(bcx), bcy),
222 Periodify( in.dfy(), R0, R1, Z0, Z1, bcx, inverse(bcy)));
223}
224inline CylindricalFunctorsLvl2 periodify( const CylindricalFunctorsLvl2& in, double R0, double R1, double Z0, double Z1, bc bcx, bc bcy)
225{
226 return CylindricalFunctorsLvl2(
227 Periodify( in.f(), R0, R1, Z0, Z1, bcx, bcy),
228 Periodify( in.dfx(), R0, R1, Z0, Z1, inverse(bcx), bcy),
229 Periodify( in.dfy(), R0, R1, Z0, Z1, bcx, inverse(bcy)),
230 Periodify( in.dfxx(), R0, R1, Z0, Z1, bcx, bcy),
231 Periodify( in.dfxy(), R0, R1, Z0, Z1, inverse(bcx), inverse(bcy)),
232 Periodify( in.dfyy(), R0, R1, Z0, Z1, bcx, bcy));
233}
235
250inline TokamakMagneticField periodify( const TokamakMagneticField& mag, double R0, double R1, double Z0, double Z1, dg::bc bcx, dg::bc bcy)
251{
252 return TokamakMagneticField( mag.R0(),
253 periodify( mag.get_psip(), R0, R1, Z0, Z1, bcx, bcy),
254 //what if Dirichlet BC in the current? Won't that generate a NaN?
255 periodify( mag.get_ipol(), R0, R1, Z0, Z1, bcx, bcy), mag.params());
256}
257
259struct LaplacePsip : public aCylindricalFunctor<LaplacePsip>
260{
261 LaplacePsip( const TokamakMagneticField& mag): m_mag(mag) { }
262 double do_compute(double R, double Z) const
263 {
264 return m_mag.psipR()(R,Z)/R+ m_mag.psipRR()(R,Z) + m_mag.psipZZ()(R,Z);
265 }
266 private:
268};
269
271struct Bmodule : public aCylindricalFunctor<Bmodule>
272{
273 Bmodule( const TokamakMagneticField& mag): m_mag(mag) { }
274 double do_compute(double R, double Z) const
275 {
276 double psipR = m_mag.psipR()(R,Z), psipZ = m_mag.psipZ()(R,Z), ipol = m_mag.ipol()(R,Z);
277 return m_mag.R0()/R*sqrt(ipol*ipol+psipR*psipR +psipZ*psipZ);
278 }
279 private:
281};
282
290struct InvB : public aCylindricalFunctor<InvB>
291{
292 InvB( const TokamakMagneticField& mag): m_mag(mag){ }
293 double do_compute(double R, double Z) const
294 {
295 double psipR = m_mag.psipR()(R,Z), psipZ = m_mag.psipZ()(R,Z), ipol = m_mag.ipol()(R,Z);
296 return R/(m_mag.R0()*sqrt(ipol*ipol + psipR*psipR +psipZ*psipZ)) ;
297 }
298 private:
300};
301
309struct LnB : public aCylindricalFunctor<LnB>
310{
311 LnB(const TokamakMagneticField& mag): m_mag(mag) { }
312 double do_compute(double R, double Z) const
313 {
314 double psipR = m_mag.psipR()(R,Z), psipZ = m_mag.psipZ()(R,Z), ipol = m_mag.ipol()(R,Z);
315 return log(m_mag.R0()/R*sqrt(ipol*ipol + psipR*psipR +psipZ*psipZ)) ;
316 }
317 private:
319};
320
331struct BR: public aCylindricalFunctor<BR>
332{
333 BR(const TokamakMagneticField& mag): m_invB(mag), m_mag(mag) { }
334 double do_compute(double R, double Z) const
335 {
336 double Rn = R/m_mag.R0();
337 double invB = m_invB(R,Z);
338 return -1./R/invB + invB/Rn/Rn*(m_mag.ipol()(R,Z)*m_mag.ipolR()(R,Z) + m_mag.psipR()(R,Z)*m_mag.psipRR()(R,Z) + m_mag.psipZ()(R,Z)*m_mag.psipRZ()(R,Z));
339 }
340 private:
341 InvB m_invB;
343};
344
353struct BZ: public aCylindricalFunctor<BZ>
354{
355 BZ(const TokamakMagneticField& mag ): m_mag(mag), m_invB(mag) { }
356 double do_compute(double R, double Z) const
357 {
358 double Rn = R/m_mag.R0();
359 return (m_invB(R,Z)/Rn/Rn)*(m_mag.ipol()(R,Z)*m_mag.ipolZ()(R,Z) + m_mag.psipR()(R,Z)*m_mag.psipRZ()(R,Z) + m_mag.psipZ()(R,Z)*m_mag.psipZZ()(R,Z));
360 }
361 private:
363 InvB m_invB;
364};
365
370struct CurvatureNablaBR: public aCylindricalFunctor<CurvatureNablaBR>
371{
372 CurvatureNablaBR(const TokamakMagneticField& mag, int sign): m_invB(mag), m_bZ(mag) {
373 if( sign >0)
374 m_sign = +1.;
375 else
376 m_sign = -1;
377 }
378 double do_compute( double R, double Z) const
379 {
380 return -m_sign*m_invB(R,Z)*m_invB(R,Z)*m_bZ(R,Z);
381 }
382 private:
383 double m_sign;
384 InvB m_invB;
385 BZ m_bZ;
386};
387
392struct CurvatureNablaBZ: public aCylindricalFunctor<CurvatureNablaBZ>
393{
394 CurvatureNablaBZ( const TokamakMagneticField& mag, int sign): m_invB(mag), m_bR(mag) {
395 if( sign >0)
396 m_sign = +1.;
397 else
398 m_sign = -1;
399 }
400 double do_compute( double R, double Z) const
401 {
402 return m_sign*m_invB(R,Z)*m_invB(R,Z)*m_bR(R,Z);
403 }
404 private:
405 double m_sign;
406 InvB m_invB;
407 BR m_bR;
408};
409
414struct CurvatureKappaR: public aCylindricalFunctor<CurvatureKappaR>
415{
417 CurvatureKappaR( const TokamakMagneticField& mag, int sign = +1){ }
418 double do_compute( double R, double Z) const
419 {
420 return 0.;
421 }
422 private:
423};
424
429struct CurvatureKappaZ: public aCylindricalFunctor<CurvatureKappaZ>
430{
431 CurvatureKappaZ( const TokamakMagneticField& mag, int sign): m_invB(mag) {
432 if( sign >0)
433 m_sign = +1.;
434 else
435 m_sign = -1;
436 }
437 double do_compute( double R, double Z) const
438 {
439 return -m_sign*m_invB(R,Z)/R;
440 }
441 private:
442 double m_sign;
443 InvB m_invB;
444};
445
450struct DivCurvatureKappa: public aCylindricalFunctor<DivCurvatureKappa>
451{
452 DivCurvatureKappa( const TokamakMagneticField& mag, int sign): m_invB(mag), m_bZ(mag){
453 if( sign >0)
454 m_sign = +1.;
455 else
456 m_sign = -1;
457 }
458 double do_compute( double R, double Z) const
459 {
460 return m_sign*m_bZ(R,Z)*m_invB(R,Z)*m_invB(R,Z)/R;
461 }
462 private:
463 double m_sign;
464 InvB m_invB;
465 BZ m_bZ;
466};
467
472struct DivCurvatureNablaB: public aCylindricalFunctor<DivCurvatureNablaB>
473{
474 DivCurvatureNablaB( const TokamakMagneticField& mag, int sign): m_div(mag, sign){ }
475 double do_compute( double R, double Z) const
476 {
477 return -m_div(R,Z);
478 }
479 private:
480 DivCurvatureKappa m_div;
481};
485struct TrueCurvatureNablaBR: public aCylindricalFunctor<TrueCurvatureNablaBR>
486{
487 TrueCurvatureNablaBR(const TokamakMagneticField& mag): m_R0(mag.R0()), m_mag(mag), m_invB(mag), m_bZ(mag) { }
488 double do_compute( double R, double Z) const
489 {
490 double invB = m_invB(R,Z), ipol = m_mag.ipol()(R,Z);
491 return -invB*invB*invB*ipol*m_R0/R*m_bZ(R,Z);
492 }
493 private:
494 double m_R0;
496 InvB m_invB;
497 BZ m_bZ;
498};
499
503struct TrueCurvatureNablaBZ: public aCylindricalFunctor<TrueCurvatureNablaBZ>
504{
505 TrueCurvatureNablaBZ(const TokamakMagneticField& mag): m_R0(mag.R0()), m_mag(mag), m_invB(mag), m_bR(mag) { }
506 double do_compute( double R, double Z) const
507 {
508 double invB = m_invB(R,Z), ipol = m_mag.ipol()(R,Z);
509 return invB*invB*invB*ipol*m_R0/R*m_bR(R,Z);
510 }
511 private:
512 double m_R0;
514 InvB m_invB;
515 BR m_bR;
516};
517
521struct TrueCurvatureNablaBP: public aCylindricalFunctor<TrueCurvatureNablaBP>
522{
523 TrueCurvatureNablaBP(const TokamakMagneticField& mag): m_mag(mag), m_invB(mag),m_bR(mag), m_bZ(mag) { }
524 double do_compute( double R, double Z) const
525 {
526 double invB = m_invB(R,Z);
527 return m_mag.R0()*invB*invB*invB/R/R*(m_mag.psipZ()(R,Z)*m_bZ(R,Z) + m_mag.psipR()(R,Z)*m_bR(R,Z));
528 }
529 private:
531 InvB m_invB;
532 BR m_bR;
533 BZ m_bZ;
534};
535
537struct TrueCurvatureKappaR: public aCylindricalFunctor<TrueCurvatureKappaR>
538{
539 TrueCurvatureKappaR( const TokamakMagneticField& mag):m_mag(mag), m_invB(mag), m_bZ(mag){ }
540 double do_compute( double R, double Z) const
541 {
542 double invB = m_invB(R,Z);
543 return m_mag.R0()*invB*invB/R*(m_mag.ipolZ()(R,Z) - m_mag.ipol()(R,Z)*invB*m_bZ(R,Z));
544 }
545 private:
547 InvB m_invB;
548 BZ m_bZ;
549};
550
552struct TrueCurvatureKappaZ: public aCylindricalFunctor<TrueCurvatureKappaZ>
553{
554 TrueCurvatureKappaZ( const TokamakMagneticField& mag):m_mag(mag), m_invB(mag), m_bR(mag){ }
555 double do_compute( double R, double Z) const
556 {
557 double invB = m_invB(R,Z);
558 return m_mag.R0()*invB*invB/R*( - m_mag.ipolR()(R,Z) + m_mag.ipol()(R,Z)*invB*m_bR(R,Z));
559 }
560 private:
562 InvB m_invB;
563 BR m_bR;
564};
566struct TrueCurvatureKappaP: public aCylindricalFunctor<TrueCurvatureKappaP>
567{
568 TrueCurvatureKappaP( const TokamakMagneticField& mag):m_mag(mag), m_invB(mag), m_bR(mag), m_bZ(mag){ }
569 double do_compute( double R, double Z) const
570 {
571 double invB = m_invB(R,Z);
572 return m_mag.R0()*invB*invB/R/R*(
573 + invB*m_mag.psipZ()(R,Z)*m_bZ(R,Z) + invB *m_mag.psipR()(R,Z)*m_bR(R,Z)
574 + m_mag.psipR()(R,Z)/R - m_mag.psipRR()(R,Z) - m_mag.psipZZ()(R,Z));
575 }
576 private:
578 InvB m_invB;
579 BR m_bR;
580 BZ m_bZ;
581};
582
584struct TrueDivCurvatureKappa: public aCylindricalFunctor<TrueDivCurvatureKappa>
585{
586 TrueDivCurvatureKappa( const TokamakMagneticField& mag): m_mag(mag), m_invB(mag), m_bR(mag), m_bZ(mag){}
587 double do_compute( double R, double Z) const
588 {
589 double invB = m_invB(R,Z);
590 return m_mag.R0()*invB*invB*invB/R*( m_mag.ipolR()(R,Z)*m_bZ(R,Z) - m_mag.ipolZ()(R,Z)*m_bR(R,Z) );
591 }
592 private:
594 InvB m_invB;
595 BR m_bR;
596 BZ m_bZ;
597};
598
600struct TrueDivCurvatureNablaB: public aCylindricalFunctor<TrueDivCurvatureNablaB>
601{
603 double do_compute( double R, double Z) const {
604 return - m_div(R,Z);
605 }
606 private:
608};
609
615struct GradLnB: public aCylindricalFunctor<GradLnB>
616{
617 GradLnB( const TokamakMagneticField& mag): m_mag(mag), m_invB(mag), m_bR(mag), m_bZ(mag) { }
618 double do_compute( double R, double Z) const
619 {
620 double invB = m_invB(R,Z);
621 return m_mag.R0()*invB*invB*(m_bR(R,Z)*m_mag.psipZ()(R,Z)-m_bZ(R,Z)*m_mag.psipR()(R,Z))/R ;
622 }
623 private:
625 InvB m_invB;
626 BR m_bR;
627 BZ m_bZ;
628};
635struct Divb: public aCylindricalFunctor<Divb>
636{
637 Divb( const TokamakMagneticField& mag): m_gradLnB(mag) { }
638 double do_compute( double R, double Z) const
639 {
640 return -m_gradLnB(R,Z);
641 }
642 private:
643 GradLnB m_gradLnB;
644};
645
647struct BFieldP: public aCylindricalFunctor<BFieldP>
648{
649 BFieldP( const TokamakMagneticField& mag): m_mag(mag){}
650 double do_compute( double R, double Z) const
651 {
652 return m_mag.R0()*m_mag.ipol()(R,Z)/R/R;
653 }
654 private:
655
657};
658
660struct BFieldR: public aCylindricalFunctor<BFieldR>
661{
662 BFieldR( const TokamakMagneticField& mag): m_mag(mag){}
663 double do_compute( double R, double Z) const
664 {
665 return m_mag.R0()/R*m_mag.psipZ()(R,Z);
666 }
667 private:
669
670};
671
673struct BFieldZ: public aCylindricalFunctor<BFieldZ>
674{
675 BFieldZ( const TokamakMagneticField& mag): m_mag(mag){}
676 double do_compute( double R, double Z) const
677 {
678 return -m_mag.R0()/R*m_mag.psipR()(R,Z);
679 }
680 private:
682};
683
685struct BFieldT: public aCylindricalFunctor<BFieldT>
686{
687 BFieldT( const TokamakMagneticField& mag): m_R0(mag.R0()), m_fieldR(mag), m_fieldZ(mag){}
688 double do_compute(double R, double Z) const
689 {
690 double r2 = (R-m_R0)*(R-m_R0) + Z*Z;
691 return m_fieldR(R,Z)*(-Z/r2) + m_fieldZ(R,Z)*(R-m_R0)/r2;
692 }
693 private:
694 double m_R0;
695 BFieldR m_fieldR;
696 BFieldZ m_fieldZ;
697};
698
700struct BHatR: public aCylindricalFunctor<BHatR>
701{
702 BHatR( const TokamakMagneticField& mag): m_mag(mag), m_invB(mag){ }
703 double do_compute( double R, double Z) const
704 {
705 return m_invB(R,Z)*m_mag.R0()/R*m_mag.psipZ()(R,Z);
706 }
707 private:
709 InvB m_invB;
710
711};
712
714struct BHatZ: public aCylindricalFunctor<BHatZ>
715{
716 BHatZ( const TokamakMagneticField& mag): m_mag(mag), m_invB(mag){ }
717 double do_compute( double R, double Z) const
718 {
719 return -m_invB(R,Z)*m_mag.R0()/R*m_mag.psipR()(R,Z);
720 }
721 private:
723 InvB m_invB;
724};
725
727struct BHatP: public aCylindricalFunctor<BHatP>
728{
729 BHatP( const TokamakMagneticField& mag): m_mag(mag), m_invB(mag){ }
730 double do_compute( double R, double Z) const
731 {
732 return m_invB(R,Z)*m_mag.R0()*m_mag.ipol()(R,Z)/R/R;
733 }
734 private:
736 InvB m_invB;
737};
738
739
749 if( sign > 0)
750 return CylindricalVectorLvl1( Constant(0), Constant(0), [](double x, double y){ return 1./x;}, Constant(0), Constant(0));
751 return CylindricalVectorLvl1( Constant(0), Constant(0), [](double x, double y){ return -1./x;}, Constant(0), Constant(0));
752}
773 return CylindricalVectorLvl0( CurvatureKappaR(mag, sign), CurvatureKappaZ(mag, sign), Constant(0));
774}
805
806
808struct BHatRR: public aCylindricalFunctor<BHatRR>
809{
810 BHatRR( const TokamakMagneticField& mag): m_invB(mag), m_br(mag), m_mag(mag){}
811 double do_compute( double R, double Z) const
812 {
813 double psipZ = m_mag.psipZ()(R,Z);
814 double psipRZ = m_mag.psipRZ()(R,Z);
815 double binv = m_invB(R,Z);
816 return -psipZ*m_mag.R0()*binv/R/R + psipRZ*binv*m_mag.R0()/R
817 -psipZ*m_mag.R0()/R*binv*binv*m_br(R,Z);
818 }
819 private:
820 InvB m_invB;
821 BR m_br;
823};
825struct BHatRZ: public aCylindricalFunctor<BHatRZ>
826{
827 BHatRZ( const TokamakMagneticField& mag): m_invB(mag), m_bz(mag), m_mag(mag){}
828 double do_compute( double R, double Z) const
829 {
830 double psipZ = m_mag.psipZ()(R,Z);
831 double psipZZ = m_mag.psipZZ()(R,Z);
832 double binv = m_invB(R,Z);
833 return m_mag.R0()/R*( psipZZ*binv -binv*binv*m_bz(R,Z)*psipZ );
834 }
835 private:
836 InvB m_invB;
837 BZ m_bz;
839};
841struct BHatZR: public aCylindricalFunctor<BHatZR>
842{
843 BHatZR( const TokamakMagneticField& mag): m_invB(mag), m_br(mag), m_mag(mag){}
844 double do_compute( double R, double Z) const
845 {
846 double psipR = m_mag.psipR()(R,Z);
847 double psipRR = m_mag.psipRR()(R,Z);
848 double binv = m_invB(R,Z);
849 return +psipR*m_mag.R0()*binv/R/R - psipRR*binv*m_mag.R0()/R
850 +psipR*m_mag.R0()/R*binv*binv*m_br(R,Z);
851 }
852 private:
853 InvB m_invB;
854 BR m_br;
856};
858struct BHatZZ: public aCylindricalFunctor<BHatZZ>
859{
860 BHatZZ( const TokamakMagneticField& mag): m_invB(mag), m_bz(mag), m_mag(mag){}
861 double do_compute( double R, double Z) const
862 {
863 double psipR = m_mag.psipR()(R,Z);
864 double psipRZ = m_mag.psipRZ()(R,Z);
865 double binv = m_invB(R,Z);
866 return -m_mag.R0()/R*( psipRZ*binv -binv*binv*m_bz(R,Z)*psipR );
867 }
868 private:
869 InvB m_invB;
870 BZ m_bz;
872};
874struct BHatPR: public aCylindricalFunctor<BHatPR>
875{
876 BHatPR( const TokamakMagneticField& mag): m_mag(mag), m_invB(mag), m_br(mag){ }
877 double do_compute( double R, double Z) const
878 {
879 double binv = m_invB(R,Z);
880 double ipol = m_mag.ipol()(R,Z);
881 double ipolR = m_mag.ipolR()(R,Z);
882 return -binv*binv*m_br(R,Z)*m_mag.R0()*ipol/R/R
883 - 2./R/R/R*binv*m_mag.R0()*ipol
884 + binv *m_mag.R0()/R/R*ipolR;
885 }
886 private:
888 InvB m_invB;
889 BR m_br;
890};
892struct BHatPZ: public aCylindricalFunctor<BHatPZ>
893{
894 BHatPZ( const TokamakMagneticField& mag): m_mag(mag), m_invB(mag), m_bz(mag){ }
895 double do_compute( double R, double Z) const
896 {
897 double binv = m_invB(R,Z);
898 double ipol = m_mag.ipol()(R,Z);
899 double ipolZ = m_mag.ipolZ()(R,Z);
900 return -binv*binv*m_bz(R,Z)*m_mag.R0()*ipol/R/R
901 + binv *m_mag.R0()/R/R*ipolZ;
902 }
903 private:
905 InvB m_invB;
906 BZ m_bz;
907};
908
910struct DivVVP: public aCylindricalFunctor<DivVVP>
911{
912 DivVVP( const TokamakMagneticField& mag): m_mag(mag),
913 m_bhatP(mag){ }
914 double do_compute( double R, double Z) const
915 {
916 double ipol = m_mag.ipol()(R,Z), ipolR = m_mag.ipolR()(R,Z),
917 ipolZ = m_mag.ipolZ()(R,Z);
918 double psipR = m_mag.psipR()(R,Z), psipZ = m_mag.psipZ()(R,Z);
919 double bphi = m_bhatP(R,Z);
920 return -(psipZ*(ipolR/R - 2.*ipol/R/R) - ipolZ/R*psipR)/
921 (ipol*ipol + psipR*psipR + psipZ*psipZ)/bphi/bphi;
922 }
923 private:
925 BHatP m_bhatP;
926};
927
935 return CylindricalVectorLvl1( BHatR(mag), BHatZ(mag), BHatP(mag),
936 Divb(mag), DivVVP(mag)
937 );
938}
939
946struct RhoP: public aCylindricalFunctor<RhoP>
947{
948 RhoP( const TokamakMagneticField& mag): m_mag(mag){
949 double RO = m_mag.R0(), ZO = 0;
950 try{
951 findOpoint( mag.get_psip(), RO, ZO);
952 m_psipmin = m_mag.psip()(RO, ZO);
953 } catch ( dg::Error& err)
954 {
955 m_psipmin = 1.;
957 m_psipmin = -10;
958 }
959 }
960 double do_compute( double R, double Z) const
961 {
962 return sqrt( 1.-m_mag.psip()(R,Z)/m_psipmin ) ;
963 }
964 private:
965 double m_psipmin;
967
968};
969
972{
974 double do_compute( double R, double Z) const
975 {
976 double psipR = m_mag.psipR()(R,Z), psipZ = m_mag.psipZ()(R,Z), ipol = m_mag.ipol()(R,Z);
977 double psip2 = psipR*psipR+psipZ*psipZ;
978 if( psip2 == 0)
979 psip2 = 1e-16;
980 return (ipol*ipol + psip2)/R/R/psip2;
981 }
982 private:
984};
985
989struct WallDirection : public dg::geo::aCylindricalFunctor<WallDirection>
990{
999 vertical, std::vector<double> horizontal) : m_vertical(vertical),
1000 m_horizontal(horizontal), m_BR( mag), m_BZ(mag){}
1008 dg::Grid2d walls) : m_vertical({walls.x0(), walls.x1()}),
1009 m_horizontal({walls.y0(), walls.y1()}), m_BR( mag), m_BZ(mag){}
1010 double do_compute ( double R, double Z) const
1011 {
1012 std::vector<double> v_dist(1,1e100), h_dist(1,1e100);
1013 for( auto v : m_vertical)
1014 v_dist.push_back( R-v );
1015 for( auto h : m_horizontal)
1016 h_dist.push_back( Z-h );
1017 double v_min = *std::min_element( v_dist.begin(), v_dist.end(),
1018 [](double a, double b){ return fabs(a) < fabs(b);} );
1019 double h_min = *std::min_element( h_dist.begin(), h_dist.end(),
1020 [](double a, double b){ return fabs(a) < fabs(b);} );
1021 if( fabs(v_min) < fabs(h_min) ) // if vertical R wall is closer
1022 {
1023 double br = m_BR( R,Z);
1024 return v_min*br < 0 ? +1 : -1;
1025 }
1026 else //horizontal Z wall is closer
1027 {
1028 double bz = m_BZ( R,Z);
1029 return h_min*bz < 0 ? +1 : -1;
1030 }
1031 }
1032 private:
1033 std::vector<double> m_vertical, m_horizontal;
1034 dg::geo::BFieldR m_BR;
1035 dg::geo::BFieldZ m_BZ;
1036};
1038
1039} //namespace geo
1040} //namespace dg
1041
bc inverse(bc bound)
modifier
How flux-function is modified.
Definition magnetic_field.h:36
CylindricalVectorLvl0 createCurvatureNablaB(const TokamakMagneticField &mag, int sign)
Approximate curvature vector field (CurvatureNablaBR, CurvatureNablaBZ, Constant(0))
Definition magnetic_field.h:761
CylindricalVectorLvl0 createGradPsip(const TokamakMagneticField &mag)
Gradient Psip vector field (PsipR, PsipZ, 0)
Definition magnetic_field.h:802
CylindricalVectorLvl0 createTrueCurvatureNablaB(const TokamakMagneticField &mag)
True curvature vector field (TrueCurvatureNablaBR, TrueCurvatureNablaBZ, TrueCurvatureNablaBP)
Definition magnetic_field.h:792
CylindricalVectorLvl0 createTrueCurvatureKappa(const TokamakMagneticField &mag)
True curvature vector field (TrueCurvatureKappaR, TrueCurvatureKappaZ, TrueCurvatureKappaP)
Definition magnetic_field.h:782
CylindricalVectorLvl0 createCurvatureKappa(const TokamakMagneticField &mag, int sign)
Approximate curvature vector field (CurvatureKappaR, CurvatureKappaZ, Constant(0))
Definition magnetic_field.h:772
equilibrium
How flux-function is computed. Decides how to construct magnetic field.
Definition magnetic_field.h:26
description
How flux function looks like. Decider on whether and what flux aligned grid to construct.
Definition magnetic_field.h:50
TokamakMagneticField periodify(const TokamakMagneticField &mag, double R0, double R1, double Z0, double Z1, dg::bc bcx, dg::bc bcy)
Use dg::geo::Periodify to periodify every function in the magnetic field.
Definition magnetic_field.h:250
CylindricalVectorLvl1 createBHat(const TokamakMagneticField &mag)
Contravariant components of the magnetic unit vector field and its Divergence and derivative in cylin...
Definition magnetic_field.h:934
CylindricalVectorLvl1 createEPhi(int sign)
Contravariant components of the unit vector field (0, 0, ) and its Divergence and derivative (0,...
Definition magnetic_field.h:748
@ none
no modification
@ sol_pfr
Psip is dampened in the SOL and PFR regions but not in the closed field line region.
@ sol_pfr_2X
Psip is dampened in the SOL and PFR regions of each of 2 X-points but not in the closed field line re...
@ heaviside
Psip is dampened to a constant outside a critical value.
@ solovev
dg::geo::solovev::Psip
@ taylor
dg::geo::taylor::Psip
@ polynomial
dg::geo::polynomial::Psip
@ guenter
dg::geo::guenter::Psip
@ toroidal
dg::geo::createToroidalField
@ circular
dg::geo::circular::Psip
@ square
closed flux surfaces centered around an O-point and bordered by a square with four X-points in the co...
@ none
no shaping: Purely toroidal magnetic field
@ standardO
closed flux surfaces centered around an O-point located near (R_0, 0); flux-aligned grids can be cons...
@ doubleX
closed flux surfaces centered around an O-point located near (R_0, 0) and bordered by a separatrix wi...
@ centeredX
one X-point in the middle, no O-point, only open flux surfaces, X-grids cannot be constructed
@ standardX
closed flux surfaces centered around an O-point located near (R_0, 0) and bordered by a separatrix wi...
int findOpoint(const CylindricalFunctorsLvl2 &psi, double &RC, double &ZC)
This function finds O-points of psi.
Definition fluxfunctions.h:336
real_type x0() const
real_type x1() const
Definition magnetic_field.h:648
double do_compute(double R, double Z) const
Definition magnetic_field.h:650
BFieldP(const TokamakMagneticField &mag)
Definition magnetic_field.h:649
Definition magnetic_field.h:661
BFieldR(const TokamakMagneticField &mag)
Definition magnetic_field.h:662
double do_compute(double R, double Z) const
Definition magnetic_field.h:663
Definition magnetic_field.h:686
BFieldT(const TokamakMagneticField &mag)
Definition magnetic_field.h:687
double do_compute(double R, double Z) const
Definition magnetic_field.h:688
Definition magnetic_field.h:674
double do_compute(double R, double Z) const
Definition magnetic_field.h:676
BFieldZ(const TokamakMagneticField &mag)
Definition magnetic_field.h:675
Definition magnetic_field.h:728
double do_compute(double R, double Z) const
Definition magnetic_field.h:730
BHatP(const TokamakMagneticField &mag)
Definition magnetic_field.h:729
Definition magnetic_field.h:875
BHatPR(const TokamakMagneticField &mag)
Definition magnetic_field.h:876
double do_compute(double R, double Z) const
Definition magnetic_field.h:877
Definition magnetic_field.h:893
double do_compute(double R, double Z) const
Definition magnetic_field.h:895
BHatPZ(const TokamakMagneticField &mag)
Definition magnetic_field.h:894
Definition magnetic_field.h:701
BHatR(const TokamakMagneticField &mag)
Definition magnetic_field.h:702
double do_compute(double R, double Z) const
Definition magnetic_field.h:703
Definition magnetic_field.h:809
BHatRR(const TokamakMagneticField &mag)
Definition magnetic_field.h:810
double do_compute(double R, double Z) const
Definition magnetic_field.h:811
Definition magnetic_field.h:826
double do_compute(double R, double Z) const
Definition magnetic_field.h:828
BHatRZ(const TokamakMagneticField &mag)
Definition magnetic_field.h:827
Definition magnetic_field.h:715
BHatZ(const TokamakMagneticField &mag)
Definition magnetic_field.h:716
double do_compute(double R, double Z) const
Definition magnetic_field.h:717
Definition magnetic_field.h:842
BHatZR(const TokamakMagneticField &mag)
Definition magnetic_field.h:843
double do_compute(double R, double Z) const
Definition magnetic_field.h:844
Definition magnetic_field.h:859
double do_compute(double R, double Z) const
Definition magnetic_field.h:861
BHatZZ(const TokamakMagneticField &mag)
Definition magnetic_field.h:860
Definition magnetic_field.h:332
double do_compute(double R, double Z) const
Definition magnetic_field.h:334
BR(const TokamakMagneticField &mag)
Definition magnetic_field.h:333
Definition magnetic_field.h:354
BZ(const TokamakMagneticField &mag)
Definition magnetic_field.h:355
double do_compute(double R, double Z) const
Definition magnetic_field.h:356
Definition magnetic_field.h:272
Bmodule(const TokamakMagneticField &mag)
Definition magnetic_field.h:273
double do_compute(double R, double Z) const
Definition magnetic_field.h:274
Definition fluxfunctions.h:114
Approximate .
Definition magnetic_field.h:415
CurvatureKappaR(const TokamakMagneticField &mag, int sign=+1)
Definition magnetic_field.h:417
CurvatureKappaR()
Definition magnetic_field.h:416
double do_compute(double R, double Z) const
Definition magnetic_field.h:418
Approximate .
Definition magnetic_field.h:430
double do_compute(double R, double Z) const
Definition magnetic_field.h:437
CurvatureKappaZ(const TokamakMagneticField &mag, int sign)
Definition magnetic_field.h:431
Approximate .
Definition magnetic_field.h:371
CurvatureNablaBR(const TokamakMagneticField &mag, int sign)
Definition magnetic_field.h:372
double do_compute(double R, double Z) const
Definition magnetic_field.h:378
Approximate .
Definition magnetic_field.h:393
CurvatureNablaBZ(const TokamakMagneticField &mag, int sign)
Definition magnetic_field.h:394
double do_compute(double R, double Z) const
Definition magnetic_field.h:400
This struct bundles a function and its first derivatives.
Definition fluxfunctions.h:186
const CylindricalFunctor & dfx() const
Definition fluxfunctions.h:209
const CylindricalFunctor & f() const
Definition fluxfunctions.h:207
const CylindricalFunctor & dfy() const
Definition fluxfunctions.h:211
This struct bundles a function and its first and second derivatives.
Definition fluxfunctions.h:223
const CylindricalFunctor & dfxy() const
Definition fluxfunctions.h:255
const CylindricalFunctor & dfy() const
Definition fluxfunctions.h:251
const CylindricalFunctor & dfx() const
Definition fluxfunctions.h:249
const CylindricalFunctor & dfxx() const
Definition fluxfunctions.h:253
const CylindricalFunctor & f() const
Definition fluxfunctions.h:247
const CylindricalFunctor & dfyy() const
Definition fluxfunctions.h:257
Definition fluxfunctions.h:416
This struct bundles a vector field and its divergence.
Definition fluxfunctions.h:444
Approximate .
Definition magnetic_field.h:451
DivCurvatureKappa(const TokamakMagneticField &mag, int sign)
Definition magnetic_field.h:452
double do_compute(double R, double Z) const
Definition magnetic_field.h:458
Approximate .
Definition magnetic_field.h:473
double do_compute(double R, double Z) const
Definition magnetic_field.h:475
DivCurvatureNablaB(const TokamakMagneticField &mag, int sign)
Definition magnetic_field.h:474
Definition magnetic_field.h:911
double do_compute(double R, double Z) const
Definition magnetic_field.h:914
DivVVP(const TokamakMagneticField &mag)
Definition magnetic_field.h:912
Definition magnetic_field.h:636
double do_compute(double R, double Z) const
Definition magnetic_field.h:638
Divb(const TokamakMagneticField &mag)
Definition magnetic_field.h:637
Definition magnetic_field.h:616
GradLnB(const TokamakMagneticField &mag)
Definition magnetic_field.h:617
double do_compute(double R, double Z) const
Definition magnetic_field.h:618
Inertia factor .
Definition magnetic_field.h:972
Hoo(dg::geo::TokamakMagneticField mag)
Definition magnetic_field.h:973
double do_compute(double R, double Z) const
Definition magnetic_field.h:974
Definition magnetic_field.h:291
InvB(const TokamakMagneticField &mag)
Definition magnetic_field.h:292
double do_compute(double R, double Z) const
Definition magnetic_field.h:293
Definition magnetic_field.h:260
double do_compute(double R, double Z) const
Definition magnetic_field.h:262
LaplacePsip(const TokamakMagneticField &mag)
Definition magnetic_field.h:261
Definition magnetic_field.h:310
double do_compute(double R, double Z) const
Definition magnetic_field.h:312
LnB(const TokamakMagneticField &mag)
Definition magnetic_field.h:311
Meta-data about the magnetic field in particular the flux function.
Definition magnetic_field.h:94
double a() const
The minor radius.
Definition magnetic_field.h:125
modifier getModifier() const
the way the flux function is modified
Definition magnetic_field.h:141
MagneticFieldParameters(double a, double elongation, double triangularity, equilibrium equ, modifier mod, description des)
Constructor.
Definition magnetic_field.h:114
description getDescription() const
how the flux function looks
Definition magnetic_field.h:143
double elongation() const
Definition magnetic_field.h:131
double triangularity() const
Definition magnetic_field.h:137
equilibrium getEquilibrium() const
the way the flux function is computed
Definition magnetic_field.h:139
MagneticFieldParameters()
Default values are for a Toroidal field.
Definition magnetic_field.h:98
Definition magnetic_field.h:947
double do_compute(double R, double Z) const
Definition magnetic_field.h:960
RhoP(const TokamakMagneticField &mag)
Definition magnetic_field.h:948
A tokamak field as given by R0, Psi and Ipol plus Meta-data like shape and equilibrium.
Definition magnetic_field.h:165
const CylindricalFunctor & psipR() const
, where R, Z are cylindrical coordinates
Definition magnetic_field.h:184
const CylindricalFunctorsLvl2 & get_psip() const
Definition magnetic_field.h:200
const CylindricalFunctor & psipRR() const
, where R, Z are cylindrical coordinates
Definition magnetic_field.h:188
TokamakMagneticField()
as long as the field stays empty the access functions are undefined
Definition magnetic_field.h:167
const CylindricalFunctor & ipol() const
the current
Definition magnetic_field.h:194
const CylindricalFunctor & ipolR() const
Definition magnetic_field.h:196
const CylindricalFunctor & psipZZ() const
, where R, Z are cylindrical coordinates
Definition magnetic_field.h:192
const CylindricalFunctorsLvl1 & get_ipol() const
Definition magnetic_field.h:201
const MagneticFieldParameters & params() const
Access Meta-data of the field.
Definition magnetic_field.h:207
const CylindricalFunctor & psipZ() const
, where R, Z are cylindrical coordinates
Definition magnetic_field.h:186
double R0() const
Definition magnetic_field.h:180
const CylindricalFunctor & ipolZ() const
Definition magnetic_field.h:198
void set(double R0, const CylindricalFunctorsLvl2 &psip, const CylindricalFunctorsLvl1 &ipol, MagneticFieldParameters gp)
Definition magnetic_field.h:171
const CylindricalFunctor & psipRZ() const
, where R, Z are cylindrical coordinates
Definition magnetic_field.h:190
const CylindricalFunctor & psip() const
, where R, Z are cylindrical coordinates
Definition magnetic_field.h:182
TokamakMagneticField(double R0, const CylindricalFunctorsLvl2 &psip, const CylindricalFunctorsLvl1 &ipol, MagneticFieldParameters gp)
Definition magnetic_field.h:168
True .
Definition magnetic_field.h:567
TrueCurvatureKappaP(const TokamakMagneticField &mag)
Definition magnetic_field.h:568
double do_compute(double R, double Z) const
Definition magnetic_field.h:569
True .
Definition magnetic_field.h:538
double do_compute(double R, double Z) const
Definition magnetic_field.h:540
TrueCurvatureKappaR(const TokamakMagneticField &mag)
Definition magnetic_field.h:539
True .
Definition magnetic_field.h:553
TrueCurvatureKappaZ(const TokamakMagneticField &mag)
Definition magnetic_field.h:554
double do_compute(double R, double Z) const
Definition magnetic_field.h:555
True .
Definition magnetic_field.h:522
double do_compute(double R, double Z) const
Definition magnetic_field.h:524
TrueCurvatureNablaBP(const TokamakMagneticField &mag)
Definition magnetic_field.h:523
True .
Definition magnetic_field.h:486
double do_compute(double R, double Z) const
Definition magnetic_field.h:488
TrueCurvatureNablaBR(const TokamakMagneticField &mag)
Definition magnetic_field.h:487
True .
Definition magnetic_field.h:504
double do_compute(double R, double Z) const
Definition magnetic_field.h:506
TrueCurvatureNablaBZ(const TokamakMagneticField &mag)
Definition magnetic_field.h:505
True .
Definition magnetic_field.h:585
double do_compute(double R, double Z) const
Definition magnetic_field.h:587
TrueDivCurvatureKappa(const TokamakMagneticField &mag)
Definition magnetic_field.h:586
True .
Definition magnetic_field.h:601
double do_compute(double R, double Z) const
Definition magnetic_field.h:603
TrueDivCurvatureNablaB(const TokamakMagneticField &mag)
Definition magnetic_field.h:602
Determine if poloidal field points towards or away from the nearest wall.
Definition magnetic_field.h:990
WallDirection(dg::geo::TokamakMagneticField mag, std::vector< double > vertical, std::vector< double > horizontal)
Allocate lines.
Definition magnetic_field.h:998
WallDirection(dg::geo::TokamakMagneticField mag, dg::Grid2d walls)
Allocate lines.
Definition magnetic_field.h:1007
double do_compute(double R, double Z) const
Definition magnetic_field.h:1010
Represent functions written in cylindrical coordinates that are independent of the angle phi serving ...
Definition fluxfunctions.h:66