Extension: Geometries
#include "dg/geometries/geometries.h"
utilities.h
Go to the documentation of this file.
1#pragma once
2#include "fluxfunctions.h"
3
5namespace dg
6{
7namespace geo
8{
9
11namespace flux{
12
19struct FieldRZYT
20{
21 //x0 and y0 sould be O-point to define angle with respect to O-point
22 FieldRZYT( const CylindricalFunctorsLvl1& psip, const CylindricalFunctorsLvl1& ipol, double R0, double Z0): R_0_(R0), Z_0_(Z0), psip_(psip), ipol_(ipol){}
23 void operator()(double t, const std::array<double,3>& y, std::array<double,3>& yp) const
24 {
25 double psipR = psip_.dfx()(y[0], y[1]), psipZ = psip_.dfy()(y[0],y[1]);
26 double ipol=ipol_.f()(y[0], y[1]);
27 yp[0] = psipZ;//fieldR
28 yp[1] = -psipR;//fieldZ
29 yp[2] =ipol/y[0];
30 double r2 = (y[0]-R_0_)*(y[0]-R_0_) + (y[1]-Z_0_)*(y[1]-Z_0_);
31 double fieldT = psipZ*(y[1]-Z_0_)/r2 + psipR*(y[0]-R_0_)/r2;
32 yp[0] /= fieldT;
33 yp[1] /= fieldT;
34 yp[2] /= fieldT;
35 }
36 private:
37 double R_0_, Z_0_;
38 CylindricalFunctorsLvl1 psip_, ipol_;
39};
40
41struct FieldRZYZ
42{
43 FieldRZYZ( const CylindricalFunctorsLvl1& psip, const CylindricalFunctorsLvl1& ipol):psip_(psip), ipol_(ipol) {}
44 void operator()(double t, const std::array<double,3>& y, std::array<double,3>& yp) const
45 {
46 double psipR = psip_.dfx()(y[0], y[1]), psipZ = psip_.dfy()(y[0],y[1]);
47 double ipol=ipol_.f()(y[0], y[1]);
48 yp[0] = psipZ;//fieldR
49 yp[1] = -psipR;//fieldZ
50 yp[2] = ipol/y[0]; //fieldYbar
51 yp[0] /= yp[1];
52 yp[2] /= yp[1];
53 yp[1] = 1.;
54 }
55 private:
56 CylindricalFunctorsLvl1 psip_, ipol_;
57};
58
64struct FieldRZY
65{
66 FieldRZY( const CylindricalFunctorsLvl2& psip, const CylindricalFunctorsLvl1& ipol): f_(1.), psip_(psip), ipol_(ipol){}
67 void set_f(double f){ f_ = f;}
68 void operator()(double t, const std::array<double,2>& y, std::array<double,2>& yp) const
69 {
70 double psipR = psip_.dfx()(y[0], y[1]), psipZ = psip_.dfy()(y[0],y[1]);
71 double ipol=ipol_.f()(y[0], y[1]);
72 double fnorm = y[0]/ipol/f_;
73 yp[0] = (psipZ)*fnorm;
74 yp[1] = -(psipR)*fnorm;
75 }
76 private:
77 double f_;
78 CylindricalFunctorsLvl1 psip_,ipol_;
79};
80
92struct FieldRZYRYZY
93{
94 FieldRZYRYZY(const CylindricalFunctorsLvl2& psip, const CylindricalFunctorsLvl1& ipol): f_(1.), f_prime_(1), psip_(psip), ipol_(ipol){}
95 void set_f( double new_f){ f_ = new_f;}
96 void set_fp( double new_fp){ f_prime_ = new_fp;}
97 void initialize( double R0, double Z0, double& yR, double& yZ)
98 {
99 double psipR = psip_.dfx()(R0, Z0), psipZ = psip_.dfy()(R0,Z0);
100 double psip2 = (psipR*psipR+ psipZ*psipZ);
101 double fnorm =R0/ipol_.f()(R0,Z0)/f_; //=Rq/I
102 yR = -psipZ/psip2/fnorm;
103 yZ = +psipR/psip2/fnorm;
104 }
105 void derive( double R0, double Z0, double& xR, double& xZ)
106 {
107 xR = +f_*psip_.dfx()(R0, Z0);
108 xZ = +f_*psip_.dfy()(R0, Z0);
109 }
110
111 void operator()(double t, const std::array<double,4>& y, std::array<double,4>& yp) const
112 {
113 double psipR = psip_.dfx()(y[0], y[1]), psipZ = psip_.dfy()(y[0],y[1]);
114 double psipRR = psip_.dfxx()(y[0], y[1]), psipRZ = psip_.dfxy()(y[0],y[1]), psipZZ = psip_.dfyy()(y[0],y[1]);
115 double ipol=ipol_.f()(y[0], y[1]);
116 double ipolR=ipol_.dfx()(y[0], y[1]);
117 double ipolZ=ipol_.dfy()(y[0], y[1]);
118 double fnorm =y[0]/ipol/f_; //=R/(I/q)
119
120 yp[0] = -(psipZ)*fnorm;
121 yp[1] = +(psipR)*fnorm;
122 yp[2] = (+psipRZ*y[2]- psipRR*y[3])*fnorm + f_prime_/f_*psipR + ipolR/ipol - 1./y[0];
123 yp[3] = (-psipRZ*y[3]+ psipZZ*y[2])*fnorm + f_prime_/f_*psipZ + ipolZ/ipol;
124
125 }
126 private:
127 double f_, f_prime_;
128 CylindricalFunctorsLvl2 psip_;
129 CylindricalFunctorsLvl1 ipol_;
130};
131
132}//namespace flux
133namespace ribeiro{
134
135struct FieldRZYT
136{
137 //x0 and y0 sould be O-point to define angle with respect to O-point
138 FieldRZYT( const CylindricalFunctorsLvl1& psip, double R0, double Z0, const CylindricalSymmTensorLvl1& chi = CylindricalSymmTensorLvl1() ): R_0_(R0), Z_0_(Z0), psip_(psip), chi_(chi){}
139 void operator()(double t, const std::array<double,3>& y, std::array<double,3>& yp) const
140 {
141 double psipR = psip_.dfx()(y[0], y[1]), psipZ = psip_.dfy()(y[0],y[1]);
142 double psip2 = chi_.xx()(y[0], y[1])*psipR*psipR
143 + 2.*chi_.xy()(y[0], y[1])*psipR*psipZ
144 + chi_.yy()(y[0], y[1])*psipZ*psipZ;
145 yp[0] = -psipZ;//fieldR
146 yp[1] = +psipR;//fieldZ
147 //yp[2] = 1; //volume
148 //yp[2] = sqrt(psip2); //equalarc
149 yp[2] = psip2; //ribeiro
150 //yp[2] = psip2*sqrt(psip2); //separatrix
151 double r2 = (y[0]-R_0_)*(y[0]-R_0_) + (y[1]-Z_0_)*(y[1]-Z_0_);
152 double fieldT = psipZ*(y[1]-Z_0_)/r2 + psipR*(y[0]-R_0_)/r2;
153 yp[0] /= fieldT;
154 yp[1] /= fieldT;
155 yp[2] /= fieldT;
156 }
157 private:
158 double R_0_, Z_0_;
159 CylindricalFunctorsLvl1 psip_;
160 CylindricalSymmTensorLvl1 chi_;
161};
162
163struct FieldRZYZ
164{
165 FieldRZYZ( const CylindricalFunctorsLvl1& psip, const CylindricalSymmTensorLvl1& chi = CylindricalSymmTensorLvl1() ): psip_(psip), chi_(chi){}
166 void operator()(double t, const std::array<double,3>& y, std::array<double,3>& yp) const
167 {
168 double psipR = psip_.dfx()(y[0], y[1]), psipZ = psip_.dfy()(y[0],y[1]);
169 double psip2 = chi_.xx()(y[0], y[1])*psipR*psipR
170 + 2.*chi_.xy()(y[0], y[1])*psipR*psipZ
171 + chi_.yy()(y[0], y[1])*psipZ*psipZ;
172 yp[0] = -psipZ;//fieldR
173 yp[1] = psipR;//fieldZ
174 //yp[2] = 1.0; //volume
175 //yp[2] = sqrt(psip2); //equalarc
176 yp[2] = psip2; //ribeiro
177 //yp[2] = psip2*sqrt(psip2); //separatrix
178 yp[0] /= yp[1];
179 yp[2] /= yp[1];
180 yp[1] = 1.;
181 }
182 private:
183 CylindricalFunctorsLvl1 psip_;
184 CylindricalSymmTensorLvl1 chi_;
185};
186
187struct FieldRZY
188{
189 FieldRZY( const CylindricalFunctorsLvl1& psip, const CylindricalSymmTensorLvl1& chi = CylindricalSymmTensorLvl1() ): f_(1.), psip_(psip), chi_(chi){}
190 void set_f(double f){ f_ = f;}
191 void operator()(double t, const std::array<double,2>& y, std::array<double,2>& yp) const
192 {
193 double psipR = psip_.dfx()(y[0], y[1]), psipZ = psip_.dfy()(y[0],y[1]);
194 double psip2 = chi_.xx()(y[0], y[1])*psipR*psipR
195 + 2.*chi_.xy()(y[0], y[1])*psipR*psipZ
196 + chi_.yy()(y[0], y[1])*psipZ*psipZ;
197 //yp[0] = +psipZ/f_;//volume
198 //yp[1] = -psipR/f_;//volume
199 //yp[0] = +psipZ/sqrt(psip2)/f_;//equalarc
200 //yp[1] = -psipR/sqrt(psip2)/f_;//equalarc
201 yp[0] = -psipZ/psip2/f_;//ribeiro
202 yp[1] = +psipR/psip2/f_;//ribeiro
203 //yp[0] = +psipZ/psip2/sqrt(psip2)/f_;//separatrix
204 //yp[1] = -psipR/psip2/sqrt(psip2)/f_;//separatrix
205 }
206 private:
207 double f_;
208 CylindricalFunctorsLvl1 psip_;
209 CylindricalSymmTensorLvl1 chi_;
210};
211
212
213struct FieldRZYRYZY
214{
215 FieldRZYRYZY( const CylindricalFunctorsLvl2& psip):psip_(psip){ f_ = f_prime_ = 1.;}
216 void set_f( double new_f){ f_ = new_f;}
217 void set_fp( double new_fp){ f_prime_ = new_fp;}
218 void initialize( double R0, double Z0, double& yR, double& yZ)
219 {
220 yR = -f_*psip_.dfy()(R0, Z0);
221 yZ = +f_*psip_.dfx()(R0, Z0);
222 }
223 void derive( double R0, double Z0, double& xR, double& xZ)
224 {
225 xR = +f_*psip_.dfx()(R0, Z0);
226 xZ = +f_*psip_.dfy()(R0, Z0);
227 }
228
229 void operator()(double t, const std::array<double,4>& y, std::array<double,4>& yp) const
230 {
231 double psipR = psip_.dfx()(y[0], y[1]), psipZ = psip_.dfy()(y[0],y[1]);
232 double psipRR = psip_.dfxx()(y[0], y[1]), psipRZ = psip_.dfxy()(y[0],y[1]), psipZZ = psip_.dfyy()(y[0],y[1]);
233 double psip2 = (psipR*psipR+ psipZ*psipZ);
234
235 yp[0] = -psipZ/f_/psip2;
236 yp[1] = +psipR/f_/psip2;
237 yp[2] = ( + psipRZ*y[2] - psipRR*y[3] )/f_/psip2
238 + f_prime_/f_* psipR + 2.*(psipR*psipRR + psipZ*psipRZ)/psip2 ;
239 yp[3] = (-psipRZ*y[3] + psipZZ*y[2])/f_/psip2
240 + f_prime_/f_* psipZ + 2.*(psipR*psipRZ + psipZ*psipZZ)/psip2;
241 }
242 private:
243 double f_, f_prime_;
244 CylindricalFunctorsLvl2 psip_;
245};
246}//namespace ribeiro
247namespace equalarc{
248
249
250struct FieldRZYT
251{
252 //x0 and y0 sould be O-point to define angle with respect to O-point
253 FieldRZYT( const CylindricalFunctorsLvl1& psip, double R0, double Z0, const CylindricalSymmTensorLvl1& chi = CylindricalSymmTensorLvl1() ): R_0_(R0), Z_0_(Z0), psip_(psip), chi_(chi){}
254 void operator()(double t, const std::array<double,3>& y, std::array<double,3>& yp) const
255 {
256 double psipR = psip_.dfx()(y[0], y[1]), psipZ = psip_.dfy()(y[0],y[1]);
257 double psip2 = chi_.xx()(y[0], y[1])*psipR*psipR
258 + 2.*chi_.xy()(y[0], y[1])*psipR*psipZ
259 + chi_.yy()(y[0], y[1])*psipZ*psipZ;
260 yp[0] = -psipZ;//fieldR
261 yp[1] = +psipR;//fieldZ
262 //yp[2] = 1; //volume
263 yp[2] = sqrt(psip2); //equalarc
264 //yp[2] = psip2; //ribeiro
265 //yp[2] = psip2*sqrt(psip2); //separatrix
266 double r2 = (y[0]-R_0_)*(y[0]-R_0_) + (y[1]-Z_0_)*(y[1]-Z_0_);
267 double fieldT = psipZ*(y[1]-Z_0_)/r2 + psipR*(y[0]-R_0_)/r2; //fieldT
268 yp[0] /= fieldT;
269 yp[1] /= fieldT;
270 yp[2] /= fieldT;
271 }
272 private:
273 double R_0_, Z_0_;
274 CylindricalFunctorsLvl1 psip_;
275 CylindricalSymmTensorLvl1 chi_;
276};
277
278struct FieldRZYZ
279{
280 FieldRZYZ( const CylindricalFunctorsLvl1& psip, const CylindricalSymmTensorLvl1& chi = CylindricalSymmTensorLvl1() ): psip_(psip), chi_(chi){}
281 void operator()(double t, const std::array<double,3>& y, std::array<double,3>& yp) const
282 {
283 double psipR = psip_.dfx()(y[0], y[1]), psipZ = psip_.dfy()(y[0],y[1]);
284 double psip2 = chi_.xx()(y[0], y[1])*psipR*psipR
285 + 2.*chi_.xy()(y[0], y[1])*psipR*psipZ
286 + chi_.yy()(y[0], y[1])*psipZ*psipZ;
287 yp[0] = -psipZ;//fieldR
288 yp[1] = +psipR;//fieldZ
289 //yp[2] = 1.0; //volume
290 yp[2] = sqrt(psip2); //equalarc
291 //yp[2] = psip2; //ribeiro
292 //yp[2] = psip2*sqrt(psip2); //separatrix
293 yp[0] /= yp[1];
294 yp[2] /= yp[1];
295 yp[1] = 1.;
296 }
297 private:
298 CylindricalFunctorsLvl1 psip_;
299 CylindricalSymmTensorLvl1 chi_;
300};
301
302struct FieldRZY
303{
304 FieldRZY( const CylindricalFunctorsLvl1& psip, const CylindricalSymmTensorLvl1& chi = CylindricalSymmTensorLvl1() ): f_(1.), psip_(psip), chi_(chi){}
305 void set_f(double f){ f_ = f;}
306 void operator()(double t, const std::array<double,2>& y, std::array<double,2>& yp) const
307 {
308 double psipR = psip_.dfx()(y[0], y[1]), psipZ = psip_.dfy()(y[0],y[1]);
309 double psip2 = chi_.xx()(y[0], y[1])*psipR*psipR
310 + 2.*chi_.xy()(y[0], y[1])*psipR*psipZ
311 + chi_.yy()(y[0], y[1])*psipZ*psipZ;
312 //yp[0] = +psipZ/f_;//volume
313 //yp[1] = -psipR/f_;//volume
314 yp[0] = -psipZ/sqrt(psip2)/f_;//equalarc
315 yp[1] = +psipR/sqrt(psip2)/f_;//equalarc
316 //yp[0] = -psipZ/psip2/f_;//ribeiro
317 //yp[1] = +psipR/psip2/f_;//ribeiro
318 //yp[0] = +psipZ/psip2/sqrt(psip2)/f_;//separatrix
319 //yp[1] = -psipR/psip2/sqrt(psip2)/f_;//separatrix
320 }
321 private:
322 double f_;
323 CylindricalFunctorsLvl1 psip_;
324 CylindricalSymmTensorLvl1 chi_;
325};
326
327
328struct FieldRZYRYZY
329{
330 FieldRZYRYZY( const CylindricalFunctorsLvl2& psip): psip_(psip){ f_ = f_prime_ = 1.;}
331 void set_f( double new_f){ f_ = new_f;}
332 void set_fp( double new_fp){ f_prime_ = new_fp;}
333 void initialize( double R0, double Z0, double& yR, double& yZ)
334 {
335 double psipR = psip_.dfx()(R0, Z0), psipZ = psip_.dfy()(R0,Z0);
336 double psip2 = (psipR*psipR+ psipZ*psipZ);
337 yR = -f_*psipZ/sqrt(psip2);
338 yZ = +f_*psipR/sqrt(psip2);
339 }
340 void derive( double R0, double Z0, double& xR, double& xZ)
341 {
342 xR = +f_*psip_.dfx()(R0, Z0);
343 xZ = +f_*psip_.dfy()(R0, Z0);
344 }
345
346 void operator()(double t, const std::array<double,4>& y, std::array<double,4>& yp) const
347 {
348 double psipR = psip_.dfx()(y[0], y[1]), psipZ = psip_.dfy()(y[0],y[1]);
349 double psipRR = psip_.dfxx()(y[0], y[1]), psipRZ = psip_.dfxy()(y[0],y[1]), psipZZ = psip_.dfyy()(y[0],y[1]);
350 double psip2 = (psipR*psipR+ psipZ*psipZ);
351
352 yp[0] = -psipZ/f_/sqrt(psip2);
353 yp[1] = +psipR/f_/sqrt(psip2);
354 yp[2] = ( +psipRZ*y[2] - psipRR *y[3])/f_/sqrt(psip2)
355 + f_prime_/f_* psipR + (psipR*psipRR + psipZ*psipRZ)/psip2 ;
356 yp[3] = ( -psipRZ*y[3] + psipZZ*y[2])/f_/sqrt(psip2)
357 + f_prime_/f_* psipZ + (psipR*psipRZ + psipZ*psipZZ)/psip2;
358 }
359 private:
360 double f_, f_prime_;
361 CylindricalFunctorsLvl2 psip_;
362};
363
364}//namespace equalarc
365
366struct FieldRZtau
367{
368 FieldRZtau(const CylindricalFunctorsLvl1& psip, const CylindricalSymmTensorLvl1& chi = CylindricalSymmTensorLvl1()): psip_(psip), chi_(chi){}
369 void operator()(double t, const std::array<double,2>& y, std::array<double,2>& yp) const
370 {
371 double psipR = psip_.dfx()(y[0], y[1]), psipZ = psip_.dfy()(y[0],y[1]);
372 double chiRR = chi_.xx()(y[0], y[1]),
373 chiRZ = chi_.xy()(y[0], y[1]),
374 chiZZ = chi_.yy()(y[0], y[1]);
375 double psip2 = chiRR*psipR*psipR + 2.*chiRZ*psipR*psipZ + chiZZ*psipZ*psipZ;
376 yp[0] = (chiRR*psipR + chiRZ*psipZ)/psip2;
377 yp[1] = (chiRZ*psipR + chiZZ*psipZ)/psip2;
378 }
379 private:
380 CylindricalFunctorsLvl1 psip_;
381 CylindricalSymmTensorLvl1 chi_;
382};
383
384struct HessianRZtau
385{
386 HessianRZtau( const CylindricalFunctorsLvl2& psip): norm_(false), quad_(1), psip_(psip){}
387 // if true goes into positive Z - direction and X else
388 void set_quadrant( int quadrant) {quad_ = quadrant;}
389 void set_norm( bool normed) {norm_ = normed;}
390 void operator()(double t, const std::array<double,2>& y, std::array<double,2>& yp) const
391 {
392 double psipRZ = psip_.dfxy()(y[0], y[1]);
393 if( psipRZ == 0)
394 {
395 if( quad_ == 0) { yp[0] = 1; yp[1] = 0; }
396 else if( quad_ == 1) { yp[0] = 0; yp[1] = 1; }
397 else if( quad_ == 2) { yp[0] = -1; yp[1] = 0; }
398 else if( quad_ == 3) { yp[0] = 0; yp[1] = -1; }
399 }
400 else
401 {
402 double psipRR = psip_.dfxx()(y[0], y[1]), psipZZ = psip_.dfyy()(y[0],y[1]);
403 double T = psipRR + psipZZ;
404 double D = psipZZ*psipRR - psipRZ*psipRZ;
405 double L1 = 0.5*T+sqrt( 0.25*T*T-D); // > 0
406 double L2 = 0.5*T-sqrt( 0.25*T*T-D); // < 0; D = L1*L2
407 if ( quad_ == 0){ yp[0] = L1 - psipZZ; yp[1] = psipRZ;}
408 else if ( quad_ == 1){ yp[0] = -psipRZ; yp[1] = psipRR - L2;}
409 else if ( quad_ == 2){ yp[0] = psipZZ - L1; yp[1] = -psipRZ;}
410 else if ( quad_ == 3){ yp[0] = +psipRZ; yp[1] = L2 - psipRR;}
411 }
412 if( norm_)
413 {
414 double vgradpsi = yp[0]*psip_.dfx()(y[0],y[1]) + yp[1]*psip_.dfy()(y[0],y[1]);
415 yp[0] /= vgradpsi, yp[1] /= vgradpsi;
416 }
417 else
418 {
419 double norm = sqrt(yp[0]*yp[0]+yp[1]*yp[1]);
420 yp[0]/= norm, yp[1]/= norm;
421 }
422
423 }
424 private:
425 bool norm_;
426 int quad_;
427 CylindricalFunctorsLvl2 psip_;
428};
429
430struct MinimalCurve
431{
432 MinimalCurve(const CylindricalFunctorsLvl1& psip): norm_(false),
433 psip_(psip){}
434 void set_norm( bool normed) {norm_ = normed;}
435 void operator()(double t, const std::array<double,4>& y, std::array<double,4>& yp) const
436 {
437 double psipR = psip_.dfx()(y[0], y[1]), psipZ = psip_.dfy()(y[0], y[1]);
438 yp[0] = y[2];
439 yp[1] = y[3];
440 //double psipRZ = psipRZ_(y[0], y[1]), psipR = psipR_(y[0], y[1]), psipZ = psipZ_(y[0], y[1]), psipRR=psipRR_(y[0], y[1]), psipZZ=psipZZ_(y[0], y[1]);
441 //double D2 = psipRR*y[2]*y[2] + 2.*psipRZ*y[2]*y[3] + psipZZ*y[3]*y[3];
442 //double grad2 = psipR*psipR+psipZ*psipZ;
443 //yp[2] = D2/(1.+grad2) * psipR ;
444 //yp[3] = D2/(1.+grad2) * psipZ ;
445 if( psip_.f()(y[0], y[1]) < 0)
446 {
447 yp[2] = -10.*psipR;
448 yp[3] = -10.*psipZ;
449 }
450 else
451 {
452 yp[2] = 10.*psipR;
453 yp[3] = 10.*psipZ;
454 }
455
456 if( norm_)
457 {
458 double vgradpsi = y[2]*psipR + y[3]*psipZ;
459 yp[0] /= vgradpsi, yp[1] /= vgradpsi, yp[2] /= vgradpsi, yp[3] /= vgradpsi;
460 }
461 }
462 private:
463 bool norm_;
464 CylindricalFunctorsLvl1 psip_;
465};
466
468
469
470namespace detail
471{
472//compute psi(x) and f(x) for given discretization of x and a fpsiMinv functor
473//doesn't integrate over the x-point
474//returns psi_1
475template <class FieldFinv>
476void construct_psi_values( FieldFinv fpsiMinv,
477 const double psi_0, const double psi_1, const double x_0, const thrust::host_vector<double>& x_vec, const double x_1,
478 thrust::host_vector<double>& psi_x,
479 thrust::host_vector<double>& f_x_, bool verbose = false)
480{
481 f_x_.resize( x_vec.size()), psi_x.resize( x_vec.size());
482 double begin(psi_0), end(begin), temp(begin);
483 unsigned N = 1;
484 double eps = 1e10, eps_old=2e10;
485 if(verbose)std::cout << "In psi function:\n";
486 double x0=x_0, x1 = psi_1>psi_0? x_vec[0]:-x_vec[0];
487 dg::SinglestepTimeloop<double> odeint( dg::RungeKutta<double>(
488 "Feagin-17-8-10",0.), fpsiMinv);
489 while( (eps < eps_old || eps > 1e-8) && eps > 1e-14) //1e-8 < eps < 1e-14
490 {
491 eps_old = eps;
492 x0 = x_0, x1 = x_vec[0];
493 if( psi_1<psi_0) x1*=-1;
494 odeint.integrate_steps( x0, begin, x1, end, N);
495 psi_x[0] = end; fpsiMinv(0.,end,temp); f_x_[0] = temp;
496 for( unsigned i=1; i<x_vec.size(); i++)
497 {
498 temp = end;
499 x0 = x_vec[i-1], x1 = x_vec[i];
500 if( psi_1<psi_0) x0*=-1, x1*=-1;
501 odeint.integrate_steps( x0, temp, x1, end, N);
502 psi_x[i] = end; fpsiMinv(0.,end,temp); f_x_[i] = temp;
503 }
504 temp = end;
505 odeint.integrate_steps(x1, temp, psi_1>psi_0?x_1:-x_1, end,N);
506 double psi_1_numerical = end;
507 eps = fabs( psi_1_numerical-psi_1);
508 if(verbose)std::cout << "Effective Psi error is "<<eps<<" with "<<N<<" steps\n";
509 N*=2;
510 }
511
512}
513
514//compute the vector of r and z - values that form one psi surface
515//assumes that the initial line is perpendicular
516template <class Fpsi, class FieldRZYRYZY>
517void compute_rzy(Fpsi fpsi, FieldRZYRYZY fieldRZYRYZY,
518 double psi, const thrust::host_vector<double>& y_vec,
519 thrust::host_vector<double>& r,
520 thrust::host_vector<double>& z,
521 thrust::host_vector<double>& yr,
522 thrust::host_vector<double>& yz,
523 thrust::host_vector<double>& xr,
524 thrust::host_vector<double>& xz,
525 double& R_0, double& Z_0, double& f, double& fp, bool verbose = false )
526{
527 thrust::host_vector<double> r_old(y_vec.size(), 0), r_diff( r_old), yr_old(r_old), xr_old(r_old);
528 thrust::host_vector<double> z_old(y_vec.size(), 0), z_diff( z_old), yz_old(r_old), xz_old(z_old);
529 r.resize( y_vec.size()), z.resize(y_vec.size()), yr.resize(y_vec.size()), yz.resize(y_vec.size()), xr.resize(y_vec.size()), xz.resize(y_vec.size());
530
531 //now compute f and starting values
532 std::array<double, 4> begin{ {0,0,0,0} }, end(begin), temp(begin);
533 const double f_psi = fpsi.construct_f( psi, begin[0], begin[1]);
534 fieldRZYRYZY.set_f(f_psi);
535 double fprime = fpsi.f_prime( psi);
536 fieldRZYRYZY.set_fp(fprime);
537 fieldRZYRYZY.initialize( begin[0], begin[1], begin[2], begin[3]);
538 R_0 = begin[0], Z_0 = begin[1];
539 if(verbose)std::cout <<f_psi<<" "<<" "<< begin[0] << " "<<begin[1]<<"\t";
540 unsigned steps = 1;
541 double eps = 1e10, eps_old=2e10;
542 using Vec = std::array<double,4>;
543 dg::SinglestepTimeloop<Vec> odeint( dg::RungeKutta<Vec>( "Feagin-17-8-10",
544 begin), fieldRZYRYZY);
545 while( eps < eps_old)
546 {
547 //begin is left const
548 eps_old = eps, r_old = r, z_old = z, yr_old = yr, yz_old = yz, xr_old = xr, xz_old = xz;
549 odeint.integrate_steps( 0, begin, y_vec[0], end, steps);
550 r[0] = end[0], z[0] = end[1], yr[0] = end[2], yz[0] = end[3];
551 fieldRZYRYZY.derive( r[0], z[0], xr[0], xz[0]);
552 //std::cout <<end[0]<<" "<< end[1] <<"\n";
553 for( unsigned i=1; i<y_vec.size(); i++)
554 {
555 temp = end;
556 odeint.integrate_steps( y_vec[i-1], temp, y_vec[i], end, steps);
557 r[i] = end[0], z[i] = end[1], yr[i] = end[2], yz[i] = end[3];
558 fieldRZYRYZY.derive( r[i], z[i], xr[i], xz[i]);
559 }
560 //compute error in R,Z only
561 dg::blas1::axpby( 1., r, -1., r_old, r_diff);
562 dg::blas1::axpby( 1., z, -1., z_old, z_diff);
563 double er = dg::blas1::dot( r_diff, r_diff);
564 double ez = dg::blas1::dot( z_diff, z_diff);
565 double ar = dg::blas1::dot( r, r);
566 double az = dg::blas1::dot( z, z);
567 eps = sqrt( er + ez)/sqrt(ar+az);
568 //std::cout << "rel. error is "<<eps<<" with "<<steps<<" steps\n";
569 steps*=2;
570 }
571 r = r_old, z = z_old, yr = yr_old, yz = yz_old, xr = xr_old, xz = xz_old;
572 f = f_psi;
573
574}
575
576} //namespace detail
577} //namespace geo
578} //namespace dg
580
dg::blas1::axpby
void axpby(get_value_type< ContainerType > alpha, const ContainerType1 &x, get_value_type< ContainerType > beta, ContainerType &y)
dg::blas1::dot
get_value_type< ContainerType1 > dot(const ContainerType1 &x, const ContainerType2 &y)
coo3d::yz
@ yz
coo3d::xz
@ xz
coo3d::z
@ z
coo2d::y
@ y