Extension: Geometries
#include "dg/geometries/geometries.h"
Loading...
Searching...
No Matches
make_field.h
Go to the documentation of this file.
1#include "magnetic_field.h"
2#include "solovev.h"
3#include "guenter.h"
4#include "polynomial.h"
5#include "toroidal.h"
6#include "taylor.h"
7#include "sheath.h"
8#include "modified.h"
9#include <dg/file/json_utilities.h>
10
15namespace dg{
16namespace geo{
17
75inline TokamakMagneticField createMagneticField( dg::file::WrappedJsonValue gs)
76{
77 std::string e = gs.get( "equilibrium", "solovev" ).asString();
79 try{
80 equi = str2equilibrium.at( e);
81 }catch ( std::out_of_range& err)
82 {
83 std::string message = "ERROR: Key \"" + e
84 + "\" not valid in field:\n\t"
85 + gs.access_string() + "\"equilibrium\" \n";
86 throw std::out_of_range(message);
87 }
88 switch( equi){
90 {
92 return createPolynomialField( gp);
93 }
95 {
96 double R0 = gs.get( "R_0", 10.0).asDouble();
97 return createToroidalField( R0);
98 }
100 {
101 double I0 = gs.get( "I_0", 20.0).asDouble();
102 double R0 = gs.get( "R_0", 10.0).asDouble();
103 return createGuenterField( R0, I0);
104 }
106 {
107 double I0 = gs.get( "I_0", 20.0).asDouble();
108 double R0 = gs.get( "R_0", 10.0).asDouble();
109 double a = gs.get( "a", 1.0).asDouble();
110 double b = gs.get( "b", 1.0).asDouble();
111 return createCircularField( R0, I0, a, b);
112 }
114 {
115 solovev::Parameters gp( gs);
116 return createTaylorField( gp);
117 }
118 default:
119 {
120 solovev::Parameters gp( gs);
121 return createSolovevField( gp);
122 }
123 }
124}
125
126
127//TODO The place where these functions are somewhat tested is geometry_diag.cpp Formalise?
128
130namespace detail{
131
132inline void transform_psi( TokamakMagneticField mag, double& psi0, double& alpha0, double& sign0)
133{
134 double RO=mag.R0(), ZO=0.;
135 dg::geo::findOpoint( mag.get_psip(), RO, ZO);
136 double psipO = mag.psip()( RO, ZO);
137 double wall_psi0p = (1.-psi0*psi0)*psipO;
138 double wall_alpha0p = -(2.*psi0+alpha0)*alpha0*psipO;
139 psi0 = wall_psi0p + sign0*wall_alpha0p/2.;
140 alpha0 = fabs( wall_alpha0p/2.);
141 sign0 = sign0*((psipO>0)-(psipO<0));
142}
143
144inline void createModifiedField(
146 dg::file::WrappedJsonValue jsmod,
147 modifier& mod,
148 CylindricalFunctorsLvl2& mod_psip,
149 CylindricalFunctor& wall, CylindricalFunctor& transition)
150{
151 std::string m = jsmod.get( "type", "heaviside" ).asString();
152 mod = modifier::heaviside;
153 description desc = mag.params().getDescription();
154 try{
155 mod = str2modifier.at( m);
156 }catch ( std::out_of_range& err)
157 {
158 std::string message = "ERROR: Key \"" + m
159 + "\" not valid in field:\n\t"
160 + jsmod.access_string() + "\"type\" \n";
161 throw std::out_of_range(message);
162 }
163 double RX1 = 0., ZX1 = 0., RX2 = 0., ZX2 = 0.;
164 double RO=mag.R0(), ZO=0.;
165 if( desc != description::none and desc != description::centeredX)
166 dg::geo::findOpoint( mag.get_psip(), RO, ZO);
167 if ( desc == description::standardX or desc == description::doubleX)
168 {
169 // Find first X-point
170 RX1 = mag.R0()-1.1*mag.params().triangularity()*mag.params().a();
171 ZX1 = -1.1*mag.params().elongation()*mag.params().a();
172 dg::geo::findXpoint( mag.get_psip(), RX1, ZX1);
173 }
174 if ( desc == description::doubleX)
175 {
176 // Find second X-point
177 RX2 = mag.R0()-1.1*mag.params().triangularity()*mag.params().a();
178 ZX2 = +1.1*mag.params().elongation()*mag.params().a();
179 dg::geo::findXpoint( mag.get_psip(), RX2, ZX2);
180 }
181 switch (mod) {
182 default: //none
183 {
184 wall = mod::DampingRegion( mod::nowhere, mag.psip(), 0, 1, -1);
185 transition = mod::MagneticTransition( mod::nowhere, mag.psip(),
186 0, 1, -1);
187 mod_psip = mag.get_psip();
188 break;
189 }
190 case modifier::heaviside:
191 {
192 double psi0 = jsmod.get( "boundary", 1.1 ).asDouble();
193 double alpha = jsmod.get( "alpha", 0.2 ).asDouble();
194 double sign = +1;
195 if( desc == description::standardX || desc == description::standardO ||
196 desc == description::doubleX)
197 detail::transform_psi( mag, psi0, alpha, sign);
198 else
199 sign = jsmod.get( "sign", -1. ).asDouble();
200
201 mod_psip = mod::createPsip( mod::everywhere, mag.get_psip(), psi0, alpha, sign);
202 wall = mod::DampingRegion( mod::everywhere, mag.psip(), psi0, alpha, -sign);
203 transition = mod::MagneticTransition( mod::everywhere, mag.psip(), psi0, alpha, sign);
204 break;
205 }
206 case modifier::sol_pfr:
207 {
208 double psi0 = jsmod["boundary"].get( 0, 1.1 ).asDouble(); //given in rho_p
209 double alpha0 = jsmod["alpha"].get(0, 0.2 ).asDouble();
210 double psi1 = jsmod["boundary"].get(1, 0.97 ).asDouble();
211 double alpha1 = jsmod[ "alpha"].get(1, 0.2 ).asDouble();
212 switch( desc){
213 // TODO Maybe this can be unified with modifier::sol_pfr_2X
214 case description::standardX:
215 {
216 double sign0 = +1., sign1 = -1.;
217 detail::transform_psi( mag, psi0, alpha0, sign0);
218 detail::transform_psi( mag, psi1, alpha1, sign1);
219 CylindricalFunctorsLvl2 mod0_psip;
220 mod0_psip = mod::createPsip(
221 mod::everywhere, mag.get_psip(), psi0, alpha0, sign0);
222 mod_psip = mod::createPsip(
223 mod::HeavisideZ( ZX1, -1), mod0_psip, psi1, alpha1, sign1);
224
225 CylindricalFunctor wall0 = mod::DampingRegion(
226 mod::everywhere, mag.psip(), psi0, alpha0, -sign0);
227 CylindricalFunctor transition0 = mod::MagneticTransition(
228 mod::everywhere, mag.psip(), psi0, alpha0, sign0);
229 CylindricalFunctor wall1 = mod::DampingRegion(
230 mod::HeavisideZ(ZX1, -1), mag.psip(), psi1, alpha1, -sign1);
231 CylindricalFunctor transition1 = mod::MagneticTransition(
232 mod::HeavisideZ(ZX1, -1), mag.psip(), psi1, alpha1, sign1);
233 wall = mod::SetUnion( wall0, wall1);
234 transition = mod::SetUnion( transition0, transition1);
235 break;
236 }
237 case description::doubleX:
238 {
239 double sign0 = +1., sign1 = -1.;
240 detail::transform_psi( mag, psi0, alpha0, sign0);
241 detail::transform_psi( mag, psi1, alpha1, sign1);
242 CylindricalFunctorsLvl2 mod0_psip, mod1_psip;
243 mod0_psip = mod::createPsip(
244 mod::everywhere, mag.get_psip(), psi0, alpha0, sign0);
245 mod1_psip = mod::createPsip(
246 mod::HeavisideZ( ZX1, -1), mod0_psip, psi1, alpha1, sign1);
247 mod_psip = mod::createPsip(
248 mod::HeavisideZ( ZX2, +1), mod1_psip, psi1, alpha1, sign1);
249 CylindricalFunctor wall0 = mod::DampingRegion(
250 mod::everywhere, mag.psip(), psi0, alpha0, -sign0);
251 CylindricalFunctor wall1 = mod::DampingRegion(
252 mod::HeavisideZ(ZX1, -1), mag.psip(), psi1, alpha1, -sign1);
253 CylindricalFunctor wall2 = mod::DampingRegion(
254 mod::HeavisideZ(ZX2, +1), mag.psip(), psi1, alpha1, -sign1);
255 CylindricalFunctor transition0 = mod::MagneticTransition(
256 mod::everywhere, mag.psip(), psi0, alpha0, sign0);
257 CylindricalFunctor transition1 = mod::MagneticTransition(
258 mod::HeavisideZ(ZX1, -1), mag.psip(), psi1, alpha1, sign1);
259 CylindricalFunctor transition2 = mod::MagneticTransition(
260 mod::HeavisideZ(ZX2, +1), mag.psip(), psi1, alpha1, sign1);
261 transition = mod::SetUnion( mod::SetUnion( transition0, transition1), transition2);
262 wall = mod::SetUnion( mod::SetUnion( wall0, wall1), wall2);
263 break;
264 }
265 default:
266 {
267 double sign0 = jsmod[ "sign"].get( 0, -1. ).asDouble();
268 double sign1 = jsmod[ "sign"].get( 1, +1. ).asDouble();
269 CylindricalFunctorsLvl2 mod0_psip;
270 mod0_psip = mod::createPsip(
271 mod::everywhere, mag.get_psip(), psi0, alpha0, sign0);
272 mod_psip = mod::createPsip(
273 mod::everywhere, mod0_psip, psi1, alpha1, sign1);
274 CylindricalFunctor wall0 = mod::DampingRegion(
275 mod::everywhere, mag.psip(), psi0, alpha0, sign0);
276 CylindricalFunctor transition0 = mod::MagneticTransition(
277 mod::everywhere, mag.psip(), psi0, alpha0, sign0);
278 CylindricalFunctor wall1 = mod::DampingRegion(
279 mod::everywhere, mag.psip(), psi1, alpha1, sign1);
280 CylindricalFunctor transition1 = mod::MagneticTransition(
281 mod::everywhere, mag.psip(), psi1, alpha1, sign1);
282 wall = mod::SetUnion( wall0, wall1);
283 transition = mod::SetUnion( transition0, transition1);
284 break;
285 }
286 }
287 break;
288 }
289 case modifier::sol_pfr_2X:
290 {
291 if( desc != description::doubleX)
292 throw Error( Message(_ping_) << "Description must be doubleX");
293 constexpr unsigned num = 4;
294 std::vector<double> psi(num), alpha(num);
295 for( unsigned u=0; u<num; u++)
296 {
297 psi[u] = jsmod["boundary"].get( u, 1.0) .asDouble();
298 alpha[u] = jsmod["alpha"].get( u, 0.1) .asDouble();
299 }
300 std::vector<double> sign = {+1,-1.,+1.,-1};
301 for( unsigned u=0; u<num; u++)
302 detail::transform_psi( mag, psi[u], alpha[u], sign[u]);
303
304 std::vector<std::function< bool( double,double)>> mods = {
305 mod::RightSideOf( {RX1,ZX1},{RO,ZO},{RX2,ZX2}),
306 mod::Above( {RX1, ZX1}, {RO,ZO}, false), // = Below
307 mod::RightSideOf( {RX2,ZX2},{RO,ZO},{RX1,ZX1}),
308 mod::Above( {RX2, ZX2}, {RO, ZO}, false) };
309
310 CylindricalFunctorsLvl2 mod_psips[num];
311 mod_psips[0] = mod::createPsip( mods[0], mag.get_psip(), psi[0],
312 alpha[0], sign[0]);
313 for( unsigned u=1; u<num; u++)
314 mod_psips[u] = mod::createPsip( mods[u], mod_psips[u-1], psi[u],
315 alpha[u], sign[u]);
316 mod_psip = mod_psips[num-1];
317
318 CylindricalFunctor walls[num], transitions[num];
319 for( unsigned u=0; u<4; u++)
320 {
321 walls[u] = mod::DampingRegion( mods[u], mag.psip(), psi[u],
322 alpha[u], -sign[u]);
323 transitions[u] = mod::MagneticTransition( mods[u], mag.psip(),
324 psi[u], alpha[u], sign[u]);
325 }
326 transition = mod::SetUnion( transitions[0], transitions[1]);
327 wall = mod::SetUnion( walls[0], walls[1]);
328 for( unsigned u=2; u<4; u++)
329 {
330 transition = mod::SetUnion( transition, transitions[u]);
331 wall = mod::SetUnion( wall, walls[u]);
332 }
333 break;
334 }
335 }
336}
337}//namespace detail
339
342
422 dg::file::WrappedJsonValue gs, dg::file::WrappedJsonValue jsmod,
423 CylindricalFunctor& wall, CylindricalFunctor& transition)
424{
427 const MagneticFieldParameters& inp = mag.params();
428
429 modifier mod;
430 detail::createModifiedField( mag, jsmod, mod, mod_psip, wall, transition);
431 description desc = inp.getDescription();
432 equilibrium equi = inp.getEquilibrium();
433 MagneticFieldParameters mod_params{ inp.a(), inp.elongation(),
434 inp.triangularity(), equi, mod, desc};
435 switch( equi)
436 {
438 {
439 solovev::Parameters gp( gs);
440 return TokamakMagneticField( gp.R_0,mod_psip,
441 solovev::createIpol( gp, mod_psip), mod_params);
442 }
443 default:
444 {
445 return TokamakMagneticField( mag.R0(), mod_psip,
446 mag.get_ipol(), mod_params);
447 }
448 }
449}
450
451
455 dg::file::WrappedJsonValue jsmod)
456{
457 CylindricalFunctor wall, transition;
459 modifier mod;
460 detail::createModifiedField( mag, jsmod, mod, mod_psip, wall, transition);
461 return wall;
462}
464inline CylindricalFunctor createWallRegion( dg::file::WrappedJsonValue gs,
465 dg::file::WrappedJsonValue jsmod)
466{
467 return createWallRegion( createMagneticField(gs), jsmod);
468}
469
470
537 dg::file::WrappedJsonValue jsmod, TokamakMagneticField mag,
538 CylindricalFunctor wall, dg::Grid2d& sheath_walls,
539 CylindricalFunctor& sheath)
540{
541 double R0 = sheath_walls.x0(), R1 = sheath_walls.x1();
542 double Z0 = sheath_walls.y0(), Z1 = sheath_walls.y1();
543 Grid1d gR1d ( R0, R1, 1, 100);
544 Grid1d gZ1d ( Z0, Z1, 1, 100);
545 std::array<bool,2> sheathR = {false,false}, sheathZ = {false,false};
546 for ( unsigned i=0; i<100; i++)
547 {
548 if( wall( R0+i*gR1d.h(), Z0) == 0)
549 sheathR[0] = true;
550 if( wall( R0+i*gR1d.h(), Z1) == 0)
551 sheathR[1] = true;
552 if( wall( R0, Z0 + i*gZ1d.h()) == 0)
553 sheathZ[0] = true;
554 if( wall( R1, Z0 + i*gZ1d.h()) == 0)
555 sheathZ[1] = true;
556 }
557 if( false == sheathR[0]) Z0 = -1e10;
558 if( false == sheathR[1]) Z1 = 1e10;
559 if( false == sheathZ[0]) R0 = -1e10;
560 if( false == sheathZ[1]) R1 = 1e10;
561 sheath_walls = dg::Grid2d( R0, R1, Z0, Z1, 1,1,1);
562 double boundary = jsmod.get( "boundary", 0.0625 ).asDouble(); // 1/16
563 double alpha = jsmod.get( "alpha", 0.015625 ).asDouble(); // 1/64
565 mag), sheath_walls, (-boundary-1e-3)*2.0*M_PI,
566 1e-6, "phi", mod::SOLRegion(mag,wall));
568 mag), sheath_walls, (+boundary+1e-3)*2.0*M_PI,
569 1e-6, "phi", mod::SOLRegion(mag,wall));
570 dg::PolynomialHeaviside polyM( -boundary*2.*M_PI + alpha*M_PI, alpha*M_PI, +1);
571 dg::PolynomialHeaviside polyP( boundary*2.*M_PI - alpha*M_PI, alpha*M_PI, -1);
572 auto sheathM = dg::compose( polyM, distM); //positive (because distance)
573 auto sheathP = dg::compose( polyP, distP);
574 sheath = mod::SetUnion( sheathM, sheathP);
575 sheath = mod::SetIntersection( mod::SetNot( wall), sheath);
576}
578
579} //namespace geo
580}//namespace dg
#define M_PI
dg::geo::TokamakMagneticField createCircularField(double R0, double I0, double a=1, double b=1)
Definition toroidal.h:136
auto compose(UnaryOp f, Functor g)
RealCylindricalFunctor< double > CylindricalFunctor
Most of the times we use double.
Definition fluxfunctions.h:50
TokamakMagneticField createMagneticField(dg::file::WrappedJsonValue gs)
Create a Magnetic field based on the given parameters.
Definition make_field.h:75
dg::RealGrid< double, 2 > Grid2d
dg::geo::TokamakMagneticField createGuenterField(double R_0, double I_0)
Create a Guenter Magnetic field.
Definition guenter.h:158
modifier
How flux-function is modified.
Definition magnetic_field.h:36
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
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
@ 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
void findXpoint(const CylindricalFunctorsLvl2 &psi, double &RC, double &ZC)
This function finds X-points of psi.
Definition fluxfunctions.h:353
int findOpoint(const CylindricalFunctorsLvl2 &psi, double &RC, double &ZC)
This function finds O-points of psi.
Definition fluxfunctions.h:335
dg::geo::TokamakMagneticField createPolynomialField(dg::geo::polynomial::Parameters gp)
Create a Polynomial Magnetic field.
Definition polynomial.h:175
dg::geo::CylindricalFunctorsLvl1 createIpol(const Parameters &gp, const CylindricalFunctorsLvl1 &psip)
Definition solovev.h:354
dg::geo::TokamakMagneticField createSolovevField(dg::geo::solovev::Parameters gp)
Create a Solovev Magnetic field.
Definition solovev.h:375
dg::geo::TokamakMagneticField createTaylorField(dg::geo::solovev::Parameters gp)
Create a Taylor Magnetic field.
Definition taylor.h:308
dg::geo::TokamakMagneticField createToroidalField(double R0)
Create a Toroidal Magnetic field.
Definition toroidal.h:118
void createSheathRegion(dg::file::WrappedJsonValue jsmod, TokamakMagneticField mag, CylindricalFunctor wall, dg::Grid2d &sheath_walls, CylindricalFunctor &sheath)
Create the sheath region where fieldlines intersect the boundary.
Definition make_field.h:536
CylindricalFunctor createWallRegion(dg::geo::TokamakMagneticField mag, dg::file::WrappedJsonValue jsmod)
A convenience function call for dg::geo::createModifiedField that ignores the transition parameter an...
Definition make_field.h:454
TokamakMagneticField createModifiedField(dg::file::WrappedJsonValue gs, dg::file::WrappedJsonValue jsmod, CylindricalFunctor &wall, CylindricalFunctor &transition)
Modify Magnetic Field and create wall above or below certain Psi values according to given parameters...
Definition make_field.h:421
real_type h(unsigned u=0) const
real_type x0() const
real_type y0() const
real_type y1() const
real_type x1() const
This struct bundles a function and its first and second derivatives.
Definition fluxfunctions.h:222
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
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
A tokamak field as given by R0, Psi and Ipol plus Meta-data like shape and equilibrium.
Definition magnetic_field.h:165
const CylindricalFunctorsLvl2 & get_psip() const
Definition magnetic_field.h:200
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
double R0() const
Definition magnetic_field.h:180
const CylindricalFunctor & psip() const
, where R, Z are cylindrical coordinates
Definition magnetic_field.h:182
Distance to wall along fieldline in phi or s coordinate
Definition sheath.h:32
The default predicate for sheath integration.
Definition modified.h:425
Definition modified.h:497
Definition modified.h:517
Definition modified.h:477
Constructs and display geometric parameters for the polynomial fields.
Definition polynomial_parameters.h:43
Constructs and display geometric parameters for the solovev and taylor fields.
Definition solovev_parameters.h:43
double R_0
major tokamak radius
Definition solovev_parameters.h:45