Extension: Geometries
#include "dg/geometries/geometries.h"
testfunctors.h
Go to the documentation of this file.
1#pragma once
2
3#include "dg/functors.h"
4#include "magnetic_field.h"
5
10namespace dg
11{
12namespace geo
13{
15//For testing purposes only
19
20//exp(R-R_0)exp(Z)cos^2 (phi)
21struct TestFunctionPsi2
22{
23 TestFunctionPsi2( const TokamakMagneticField& c, double a):R_0(c.R0()), a(a), c_(c){}
24 double operator()(double R, double Z, double phi) const {
25 return exp((R-R_0)/a)*exp(Z/a)*cos(phi)*cos(phi);
26 }
27 double dR( double R, double Z, double phi) const{
28 return exp((R-R_0)/a)*exp(Z/a)*cos(phi)*cos(phi)/a;
29 }
30 double dRR( double R, double Z, double phi) const{
31 return exp((R-R_0)/a)*exp(Z/a)*cos(phi)*cos(phi)/a/a;
32 }
33 double dRZ( double R, double Z, double phi) const{
34 return exp((R-R_0)/a)*exp(Z/a)*cos(phi)*cos(phi)/a/a;
35 }
36 double dZ( double R, double Z, double phi) const{
37 return exp((R-R_0)/a)*exp(Z/a)*cos(phi)*cos(phi)/a;
38 }
39 double dZZ( double R, double Z, double phi) const{
40 return exp((R-R_0)/a)*exp(Z/a)*cos(phi)*cos(phi)/a/a;
41 }
42 double dP( double R, double Z, double phi) const{
43 return exp((R-R_0)/a)*exp(Z/a)*(-2.*sin(phi)*cos(phi));
44 }
45 double dRP( double R, double Z, double phi) const{
46 return exp((R-R_0)/a)*exp(Z/a)*(-2.*sin(phi)*cos(phi))/a;
47 }
48 double dZP( double R, double Z, double phi) const{
49 return exp((R-R_0)/a)*exp(Z/a)*(-2.*sin(phi)*cos(phi))/a;
50 }
51 double dPP( double R, double Z, double phi) const{
52 return exp((R-R_0)/a)*exp(Z/a)*(2.*sin(phi)*sin(phi)-2.*cos(phi)*cos(phi));
53 }
54 private:
55 double R_0, a;
56 TokamakMagneticField c_;
57};
58
59// (cos(M_PI*(R-R_0))+1)*(cos(M_PI*Z/2.)+1)*sin(phi);
60struct TestFunctionDirNeu{
61 TestFunctionDirNeu( const TokamakMagneticField& c){
62 R_0 = c.R0();
63 }
64 double operator()(double R, double Z, double phi)const{
65 return (cos(M_PI*(R-R_0))+1.)*(cos(M_PI*Z)+1.)*sin(phi);
66 }
67 double dR(double R, double Z, double phi)const{
68 return -M_PI*sin(M_PI*(R-R_0))*(cos(M_PI*Z)+1.)*sin(phi);
69 }
70 double dZ(double R, double Z, double phi)const{
71 return -M_PI*(cos(M_PI*(R-R_0))+1.)*sin(M_PI*Z)*sin(phi);
72 }
73 double dP(double R, double Z, double phi)const{
74 return (cos(M_PI*(R-R_0))+1.)*(cos(M_PI*Z)+1.)*cos(phi);
75 }
76 double dRR(double R, double Z, double phi)const{
77 return -M_PI*M_PI*cos(M_PI*(R-R_0))*(cos(M_PI*Z)+1.)*sin(phi);
78 }
79 double dZZ(double R, double Z, double phi)const{
80 return -M_PI*M_PI*(cos(M_PI*(R-R_0))+1.)*cos(M_PI*Z)*sin(phi);
81 }
82 double dPP(double R, double Z, double phi)const{
83 return -(cos(M_PI*(R-R_0))+1.)*(cos(M_PI*Z)+1.)*sin(phi);
84 }
85 double dRZ(double R, double Z, double phi)const{
86 return M_PI*M_PI*sin(M_PI*(R-R_0))*sin(M_PI*Z)*sin(phi);
87 }
88 double dRP(double R, double Z, double phi)const{
89 return -M_PI*sin(M_PI*(R-R_0))*(cos(M_PI*Z)+1.)*cos(phi);
90 }
91 double dZP(double R, double Z, double phi)const{
92 return -M_PI*(cos(M_PI*(R-R_0))+1.)*sin(M_PI*Z)*cos(phi);
93 }
94 private:
95 double R_0;
96};
97
99// b \nabla f
100template<class Function>
101struct DsFunction
102{
103 DsFunction( const TokamakMagneticField& c, Function f): f_(f), c_(c),
104 bhatR_(c), bhatZ_(c), bhatP_(c){}
105 double operator()(double R, double Z, double phi) const {
106 return bhatR_(R,Z)*f_.dR(R,Z,phi) +
107 bhatZ_(R,Z)*f_.dZ(R,Z,phi) +
108 bhatP_(R,Z)*f_.dP(R,Z,phi);
109 }
110 private:
111 Function f_;
112 TokamakMagneticField c_;
113 dg::geo::BHatR bhatR_;
114 dg::geo::BHatZ bhatZ_;
115 dg::geo::BHatP bhatP_;
116};
117//\nabla( b f)
118template<class Function>
119struct DsDivFunction
120{
121 DsDivFunction( const TokamakMagneticField& c, Function f):
122 f_(f), dsf_(c,f), divb_(c){}
123 double operator()(double R, double Z, double phi) const {
124 return f_(R,Z,phi)*divb_(R,Z) + dsf_(R,Z,phi);
125 }
126 private:
127 Function f_;
128 DsFunction<Function> dsf_;
129 dg::geo::Divb divb_;
130};
131
132//2nd derivative \nabla_\parallel^2
133template<class Function>
134struct DssFunction
135{
136 DssFunction( TokamakMagneticField c, Function f):f_(f), c_(c),
137 bhatR_(c), bhatZ_(c), bhatP_(c),
138 bhatRR_(c), bhatZR_(c), bhatPR_(c),
139 bhatRZ_(c), bhatZZ_(c), bhatPZ_(c){}
140 double operator()(double R, double Z, double phi) const {
141 double bhatR = bhatR_(R,Z), bhatZ = bhatZ_(R,Z), bhatP = bhatP_(R,Z);
142 double bhatRR = bhatRR_(R,Z), bhatZR = bhatZR_(R,Z), bhatPR = bhatPR_(R,Z);
143 double bhatRZ = bhatRZ_(R,Z), bhatZZ = bhatZZ_(R,Z), bhatPZ = bhatPZ_(R,Z);
144 double fR = f_.dR(R,Z,phi), fZ = f_.dZ(R,Z,phi), fP = f_.dP(R,Z,phi);
145 double fRR = f_.dRR(R,Z,phi), fRZ = f_.dRZ(R,Z,phi), fZZ = f_.dZZ(R,Z,phi);
146 double fRP = f_.dRP(R,Z,phi), fZP = f_.dZP(R,Z,phi), fPP = f_.dPP(R,Z,phi);
147 double gradbhatR = bhatR*bhatRR+bhatZ*bhatRZ,
148 gradbhatZ = bhatR*bhatZR+bhatZ*bhatZZ,
149 gradbhatP = bhatR*bhatPR+bhatZ*bhatPZ;
150 return bhatR*bhatR*fRR + bhatZ*bhatZ*fZZ + bhatP*bhatP*fPP
151 +2.*(bhatR*bhatZ*fRZ + bhatR*bhatP*fRP + bhatZ*bhatP*fZP)
152 + gradbhatR*fR + gradbhatZ*fZ + gradbhatP*fP;
153 }
154 private:
155 Function f_;
156 TokamakMagneticField c_;
157 dg::geo::BHatR bhatR_;
158 dg::geo::BHatZ bhatZ_;
159 dg::geo::BHatP bhatP_;
160 dg::geo::BHatRR bhatRR_;
161 dg::geo::BHatZR bhatZR_;
162 dg::geo::BHatPR bhatPR_;
163 dg::geo::BHatRZ bhatRZ_;
164 dg::geo::BHatZZ bhatZZ_;
165 dg::geo::BHatPZ bhatPZ_;
166};
167
168//positive Laplacian \Delta_\parallel
169template<class Function>
170struct DsDivDsFunction
171{
172 DsDivDsFunction( const TokamakMagneticField& c,Function f): dsf_(c,f), dssf_(c,f), divb_(c){}
173 double operator()(double R, double Z, double phi) const {
174 return divb_(R,Z)*dsf_(R,Z,phi) + dssf_(R,Z,phi);
175 }
176 private:
177 DsFunction<Function> dsf_;
178 DssFunction<Function> dssf_;
179 dg::geo::Divb divb_;
180};
181
182//positive perp Laplacian \Delta_\perp
183template<class Function>
184struct DPerpFunction
185{
186 DPerpFunction( const TokamakMagneticField& c, Function f): f_(f), dsf_(c,f){}
187 double operator()(double R, double Z, double phi) const {
188 return f_.dR(R,Z,phi)/R + f_.dRR(R,Z,phi) + f_.dZZ(R,Z,phi) + f_.dPP(R,Z,phi)/R/R - dsf_(R,Z,phi);
189 }
190 private:
191 Function f_;
192 DsDivDsFunction<Function> dsf_;
193};
194
195template<class Function>
196struct OMDsDivDsFunction
197{
198 OMDsDivDsFunction( const TokamakMagneticField& c,Function f): f_(f), df_(c,f){}
199 double operator()(double R, double Z, double phi) const {
200 return f_(R,Z,phi)-df_(R,Z,phi);
201 }
202 private:
203 Function f_;
204 DsDivDsFunction<Function> df_;
205};
206
207template<class Function>
208struct Variation
209{
210 Variation( Function f): f_(f){}
211 double operator()(double R, double Z, double phi) const {
212 return sqrt(f_.dR(R,Z,phi)*f_.dR(R,Z,phi) + f_.dZ(R,Z,phi)*f_.dZ(R,Z,phi));
213 }
214 private:
215 Function f_;
216};
217
219template<class DS, class container>
220void callDS( DS& ds, std::string name, const container& in, container& out,
221unsigned max_iter = 1e4, double eps = 1e-6)
222{
223 if( name == "forward") ds.ds( dg::forward, in, out);
224 else if( name == "backward") ds.ds( dg::backward, in, out);
225 else if( name == "forward2") ds.forward2( 1., in, 0., out);
226 else if( name == "backward2") ds.backward2( 1., in, 0., out);
227 else if( name == "centered") ds.ds( dg::centered, in, out);
228 else if( name == "dss") ds.dss( in, out);
229 else if( name == "centered_bc_along")
230 ds.centered_bc_along_field( 1., in, 0., out, ds.fieldaligned().bcx(), {0,0});
231 else if( name == "dss_bc_along")
232 ds.dss_bc_along_field( 1., in, 0., out, ds.fieldaligned().bcx(), {0,0});
233 else if( name == "centered") ds.ds( dg::centered, in, out);
234 else if( name == "dss") ds.dss( in, out);
235 else if( name == "divForward") ds.div( dg::forward, in, out);
236 else if( name == "divBackward") ds.div( dg::backward, in, out);
237 else if( name == "divCentered") ds.div( dg::centered, in, out);
238 else if( name == "directLap") {
239 ds.dssd( 1., in, 0., out);
240 }
241 else if( name == "directLap_bc_along") {
242 ds.dssd_bc_along_field( 1., in, 0., out, ds.fieldaligned().bcx(), {0,0});
243 }
244 else if( name == "invCenteredLap"){
245 //dg::LGMRES<container> invert( in, 30,3,10000);
246 dg::BICGSTABl<container> invert( in, 30000,3);
247 dg::Timer t;
248 t.tic();
249 double precond = 1.;
250 unsigned number = invert.solve( [&](const auto& x, auto& y){
251 // y = ( 1 - D) x
252 dg::blas2::symv( ds, x, y);
253 dg::blas1::axpby( 1., x, -1., y, y);
254 }, out, in, precond, ds.weights(), eps);
255 t.toc();
256#ifdef MPI_VERSION
257 int rank;
258 MPI_Comm_rank(MPI_COMM_WORLD, &rank);
259 if(rank==0)
260 {
261#endif //MPI
262 std::cout << "#Number of BICGSTABl iterations: "<<number<<"\n";
263 std::cout << "#Took : "<<t.diff()<<"\n";
264#ifdef MPI_VERSION
265 }
266#endif //MPI
267 return;
268 }
269
270}
272
273//psi * cos(theta)
274struct FuncDirPer
275{
276 FuncDirPer( const TokamakMagneticField& c, double psi_0, double psi_1, double k):
277 R_0_(c.R0()), psi0_(psi_0), psi1_(psi_1), k_(k), c_(c) {}
278 double operator()(double R, double Z) const {
279 double psip = c_.psip()(R,Z);
280 double result = (psip-psi0_)*(psip-psi1_)*cos(k_*theta(R,Z));
281 return 0.1*result;
282 }
283 double operator()(double R, double Z, double phi) const {
284 return operator()(R,Z);
285 }
286 double dR( double R, double Z)const
287 {
288 double psip = c_.psip()(R,Z), psipR = c_.psipR()(R,Z), theta_ = k_*theta(R,Z);
289 double result = (2.*psip*psipR - (psi0_+psi1_)*psipR)*cos(theta_)
290 - (psip-psi0_)*(psip-psi1_)*sin(theta_)*k_*thetaR(R,Z);
291 return 0.1*result;
292 }
293 double dRR( double R, double Z)const
294 {
295 double psip = c_.psip()(R,Z), psipR = c_.psipR()(R,Z), theta_=k_*theta(R,Z), thetaR_=k_*thetaR(R,Z);
296 double psipRR = c_.psipRR()(R,Z);
297 double result = (2.*(psipR*psipR + psip*psipRR) - (psi0_+psi1_)*psipRR)*cos(theta_)
298 - 2.*(2.*psip*psipR-(psi0_+psi1_)*psipR)*sin(theta_)*thetaR_
299 - (psip-psi0_)*(psip-psi1_)*(k_*thetaRR(R,Z)*sin(theta_)+cos(theta_)*thetaR_*thetaR_);
300 return 0.1*result;
301
302 }
303 double dZ( double R, double Z)const
304 {
305 double psip = c_.psip()(R,Z), psipZ = c_.psipZ()(R,Z), theta_=k_*theta(R,Z);
306 double result = (2*psip*psipZ - (psi0_+psi1_)*psipZ)*cos(theta_)
307 - (psip-psi0_)*(psip-psi1_)*sin(theta_)*k_*thetaZ(R,Z);
308 return 0.1*result;
309 }
310 double dZZ( double R, double Z)const
311 {
312 double psip = c_.psip()(R,Z), psipZ = c_.psipZ()(R,Z), theta_=k_*theta(R,Z), thetaZ_=k_*thetaZ(R,Z);
313 double psipZZ = c_.psipZZ()(R,Z);
314 double result = (2.*(psipZ*psipZ + psip*psipZZ) - (psi0_+psi1_)*psipZZ)*cos(theta_)
315 - 2.*(2.*psip*psipZ-(psi0_+psi1_)*psipZ)*sin(theta_)*thetaZ_
316 - (psip-psi0_)*(psip-psi1_)*(k_*thetaZZ(R,Z)*sin(theta_) + cos(theta_)*thetaZ_*thetaZ_ );
317 return 0.1*result;
318 }
319 private:
320 double theta( double R, double Z) const {
321 double dR = R-R_0_;
322 if( Z >= 0)
323 return acos( dR/sqrt( dR*dR + Z*Z));
324 else
325 return 2.*M_PI-acos( dR/sqrt( dR*dR + Z*Z));
326 }
327 double thetaR( double R, double Z) const {
328 double dR = R-R_0_;
329 return -Z/(dR*dR+Z*Z);
330 }
331 double thetaZ( double R, double Z) const {
332 double dR = R-R_0_;
333 return dR/(dR*dR+Z*Z);
334 }
335 double thetaRR( double R, double Z) const {
336 double dR = R-R_0_;
337 return 2*Z*dR/(dR*dR+Z*Z)/(dR*dR+Z*Z);
338 }
339 double thetaZZ( double R, double Z) const { return -thetaRR(R,Z);}
340 double R_0_;
341 double psi0_, psi1_, k_;
342 const TokamakMagneticField c_;
343};
344// Variation of FuncDirPer
345struct VariationDirPer
346{
347 VariationDirPer( dg::geo::TokamakMagneticField mag, double psi_0, double psi_1): m_f(mag, psi_0, psi_1,4. ){}
348 double operator()(double R, double Z, double phi) const {
349 return this->operator()(R,Z);}
350
351 double operator()(double R, double Z) const {
352 return m_f.dR( R,Z)*m_f.dR(R,Z) + m_f.dZ(R,Z)*m_f.dZ(R,Z);
353 }
354 private:
355 dg::geo::FuncDirPer m_f;
356};
357
358
359//takes the magnetic field as chi
360struct EllipticDirPerM
361{
362 EllipticDirPerM( const TokamakMagneticField& c, double psi_0, double psi_1, double k): func_(c, psi_0, psi_1, k), bmod_(c), br_(c), bz_(c) {}
363 double operator()(double R, double Z, double phi) const {
364 return operator()(R,Z);}
365 double operator()(double R, double Z) const {
366 double bmod = bmod_(R,Z), br = br_(R,Z), bz = bz_(R,Z);
367 return -(br*func_.dR(R,Z) + bz*func_.dZ(R,Z) + bmod*(func_.dRR(R,Z) + func_.dZZ(R,Z) ));
368
369 }
370 private:
371 FuncDirPer func_;
372 dg::geo::Bmodule bmod_;
373 dg::geo::BR br_;
374 dg::geo::BZ bz_;
375};
376
377//Blob function
378struct FuncDirNeu
379{
380 FuncDirNeu( const TokamakMagneticField& c, double psi_0, double psi_1, double R_blob, double Z_blob, double sigma_blob, double amp_blob):
381 psi0_(psi_0), psi1_(psi_1),
382 cauchy_(R_blob, Z_blob, sigma_blob, sigma_blob, amp_blob){}
383
384 double operator()(double R, double Z, double phi) const {
385 return operator()(R,Z);}
386 double operator()(double R, double Z) const {
387 return cauchy_(R,Z);
388 //double psip = psip_(R,Z);
389 //return (psip-psi0_)*(psip-psi1_)+cauchy_(R,Z);
390 //return (psip-psi0_)*(psip-psi1_);
391 }
392 double dR( double R, double Z)const
393 {
394 return cauchy_.dx(R,Z);
395 //double psip = psip_(R,Z), psipR = psipR_(R,Z);
396 //return (2.*psip-psi0_-psi1_)*psipR + cauchy_.dx(R,Z);
397 //return (2.*psip-psi0_-psi1_)*psipR;
398 }
399 double dRR( double R, double Z)const
400 {
401 return cauchy_.dxx(R,Z);
402 //double psip = psip_(R,Z), psipR = psipR_(R,Z);
403 //double psipRR = psipRR_(R,Z);
404 //return (2.*(psipR*psipR + psip*psipRR) - (psi0_+psi1_)*psipRR)+cauchy_.dxx(R,Z);
405 //return (2.*(psipR*psipR + psip*psipRR) - (psi0_+psi1_)*psipRR);
406
407 }
408 double dZ( double R, double Z)const
409 {
410 return cauchy_.dy(R,Z);
411 //double psip = psip_(R,Z), psipZ = psipZ_(R,Z);
412 //return (2*psip-psi0_-psi1_)*psipZ+cauchy_.dy(R,Z);
413 //return (2*psip-psi0_-psi1_)*psipZ;
414 }
415 double dZZ( double R, double Z)const
416 {
417 return cauchy_.dyy(R,Z);
418 //double psip = psip_(R,Z), psipZ = psipZ_(R,Z);
419 //double psipZZ = psipZZ_(R,Z);
420 //return (2.*(psipZ*psipZ + psip*psipZZ) - (psi0_+psi1_)*psipZZ)+cauchy_.dyy(R,Z);
421 //return (2.*(psipZ*psipZ + psip*psipZZ) - (psi0_+psi1_)*psipZZ);
422 }
423 private:
424 double psi0_, psi1_;
425 dg::Cauchy cauchy_;
426};
427
428
429//takes the magnetic field multiplied by (1+0.5sin(theta)) as chi
430struct BmodTheta
431{
432 BmodTheta( const TokamakMagneticField& c): R_0_(c.R0()), bmod_(c){}
433 double operator()(double R,double Z, double phi) const{
434 return operator()(R,Z);}
435 double operator()(double R,double Z) const{
436 return bmod_(R,Z)*(1.+0.5*sin(theta(R,Z)));
437 }
438 private:
439 double theta( double R, double Z) const {
440 double dR = R-R_0_;
441 if( Z >= 0)
442 return acos( dR/sqrt( dR*dR + Z*Z));
443 else
444 return 2.*M_PI-acos( dR/sqrt( dR*dR + Z*Z));
445 }
446 double R_0_;
447 dg::geo::Bmodule bmod_;
448
449};
450
451//take BmodTheta as chi
452struct EllipticDirNeuM
453{
454 EllipticDirNeuM( const TokamakMagneticField& c, double psi_0, double psi_1, double R_blob, double Z_blob, double sigma_blob, double amp_blob): R_0_(c.R0()),
455 func_(c, psi_0, psi_1, R_blob, Z_blob, sigma_blob,amp_blob), bmod_(c), br_(c), bz_(c) {}
456 double operator()(double R, double Z) const {
457 double bmod = bmod_(R,Z), br = br_(R,Z), bz = bz_(R,Z), theta_ = theta(R,Z);
458 double chi = bmod*(1.+0.5*sin(theta_));
459 double chiR = br*(1.+0.5*sin(theta_)) + bmod*0.5*cos(theta_)*thetaR(R,Z);
460 double chiZ = bz*(1.+0.5*sin(theta_)) + bmod*0.5*cos(theta_)*thetaZ(R,Z);
461 return -(chiR*func_.dR(R,Z) + chiZ*func_.dZ(R,Z) + chi*( func_.dRR(R,Z) + func_.dZZ(R,Z) ));
462
463 }
464 double operator()(double R, double Z, double phi) const {
465 return operator()(R,Z);
466 }
467 private:
468 double theta( double R, double Z) const {
469 double dR = R-R_0_;
470 if( Z >= 0)
471 return acos( dR/sqrt( dR*dR + Z*Z));
472 else
473 return 2.*M_PI-acos( dR/sqrt( dR*dR + Z*Z));
474 }
475 double thetaR( double R, double Z) const {
476 double dR = R-R_0_;
477 return -Z/(dR*dR+Z*Z);
478 }
479 double thetaZ( double R, double Z) const {
480 double dR = R-R_0_;
481 return dR/(dR*dR+Z*Z);
482 }
483 double R_0_;
484 FuncDirNeu func_;
485 dg::geo::Bmodule bmod_;
486 dg::geo::BR br_;
487 dg::geo::BZ bz_;
488};
489
490//the psi surfaces
491struct FuncXDirNeu
492{
493 FuncXDirNeu( const TokamakMagneticField& c, double psi_0, double psi_1):
494 c_(c), psi0_(psi_0), psi1_(psi_1){}
495
496 double operator()(double R, double Z, double phi) const {
497 return operator()(R,Z);}
498 double operator()(double R, double Z) const {
499 double psip = c_.psip()(R,Z);
500 return (psip-psi0_)*(psip-psi1_);
501 }
502 double dR( double R, double Z)const
503 {
504 double psip = c_.psip()(R,Z), psipR = c_.psipR()(R,Z);
505 return (2.*psip-psi0_-psi1_)*psipR;
506 }
507 double dRR( double R, double Z)const
508 {
509 double psip = c_.psip()(R,Z), psipR = c_.psipR()(R,Z);
510 double psipRR = c_.psipRR()(R,Z);
511 return (2.*(psipR*psipR + psip*psipRR) - (psi0_+psi1_)*psipRR);
512
513 }
514 double dZ( double R, double Z)const
515 {
516 double psip = c_.psip()(R,Z), psipZ = c_.psipZ()(R,Z);
517 return (2*psip-psi0_-psi1_)*psipZ;
518 }
519 double dZZ( double R, double Z)const
520 {
521 double psip = c_.psip()(R,Z), psipZ = c_.psipZ()(R,Z);
522 double psipZZ = c_.psipZZ()(R,Z);
523 return (2.*(psipZ*psipZ + psip*psipZZ) - (psi0_+psi1_)*psipZZ);
524 }
525 private:
526 TokamakMagneticField c_;
527 double psi0_, psi1_;
528};
529
530//take Bmod as chi
531struct EllipticXDirNeuM
532{
533 EllipticXDirNeuM( const TokamakMagneticField& c, double psi_0, double psi_1): R_0_(c.R0()),
534 func_(c, psi_0, psi_1), bmod_(c), br_(c), bz_(c) {}
535 double operator()(double R, double Z) const {
536 double bmod = bmod_(R,Z), br = br_(R,Z), bz = bz_(R,Z);
537 //double chi = 1e4+bmod; //bmod can be zero for a Taylor state(!)
538 double chi = bmod; //bmod for solovev state
539 double chiR = br;
540 double chiZ = bz;
541 return -(chiR*func_.dR(R,Z) + chiZ*func_.dZ(R,Z) + chi*( func_.dRR(R,Z) + func_.dZZ(R,Z) ));
542 //return -( func_.dRR(R,Z) + func_.dZZ(R,Z) );
543
544 }
545 double operator()(double R, double Z, double phi) const {
546 return operator()(R,Z);
547 }
548 private:
549 double R_0_;
550 FuncXDirNeu func_;
551 dg::geo::Bmodule bmod_;
552 dg::geo::BR br_;
553 dg::geo::BZ bz_;
554};
555
556//take Blob and chi=1
557struct EllipticBlobDirNeuM
558{
559 EllipticBlobDirNeuM( const TokamakMagneticField& c, double psi_0, double psi_1, double R_blob, double Z_blob, double sigma_blob, double amp_blob):
560 func_(c, psi_0, psi_1, R_blob, Z_blob, sigma_blob, amp_blob){}
561 double operator()(double R, double Z) const {
562 return -( func_.dRR(R,Z) + func_.dZZ(R,Z) );
563 }
564 double operator()(double R, double Z, double phi) const {
565 return operator()(R,Z);
566 }
567 private:
568 double R_0_;
569 FuncDirNeu func_;
570};
571
572struct EllipticDirSimpleM
573{
574 EllipticDirSimpleM( const TokamakMagneticField& c, double psi_0, double psi_1, double R_blob, double Z_blob, double sigma_blob, double amp_blob): func_(c, psi_0, psi_1, R_blob, Z_blob, sigma_blob, amp_blob) {}
575 double operator()(double R, double Z, double phi) const {
576 return -(( 1./R*func_.dR(R,Z) + func_.dRR(R,Z) + func_.dZZ(R,Z) ));
577
578 }
579 private:
580 FuncDirNeu func_;
581};
582
585//
586} //namespace functors
587} //namespace dg
M_PI
#define M_PI
dg::blas1::axpby
void axpby(get_value_type< ContainerType > alpha, const ContainerType1 &x, get_value_type< ContainerType > beta, ContainerType &y)
dg::blas2::symv
void symv(MatrixType &&M, const ContainerType1 &x, ContainerType2 &y)
dg::backward
backward
dg::forward
forward
dg::centered
centered
dg::invert
dg::Operator< T > invert(const dg::Operator< T > &in)
magnetic_field.h
dg::Timer::diff
double diff() const
dg::Timer::toc
void toc()
dg::Timer::tic
void tic()
dg::geo::BHatP
Definition: magnetic_field.h:725
dg::geo::BHatPR
Definition: magnetic_field.h:872
dg::geo::BHatPZ
Definition: magnetic_field.h:890
dg::geo::BHatR
Definition: magnetic_field.h:698
dg::geo::BHatRR
Definition: magnetic_field.h:806
dg::geo::BHatRZ
Definition: magnetic_field.h:823
dg::geo::BHatZ
Definition: magnetic_field.h:712
dg::geo::BHatZR
Definition: magnetic_field.h:839
dg::geo::BHatZZ
Definition: magnetic_field.h:856
dg::geo::BR
Definition: magnetic_field.h:329
dg::geo::BZ
Definition: magnetic_field.h:351
dg::geo::Bmodule
Definition: magnetic_field.h:269
dg::geo::Divb
Definition: magnetic_field.h:633
dg::geo::TokamakMagneticField
A tokamak field as given by R0, Psi and Ipol plus Meta-data like shape and equilibrium.
Definition: magnetic_field.h:162