Extension: Geometries
#include "dg/geometries/geometries.h"
flux.h
Go to the documentation of this file.
1#pragma once
2
3#include "dg/algorithm.h"
4#include "fluxfunctions.h"
5#include "ribeiro.h"
6
7
8namespace dg
9{
10namespace geo
11{
12
14namespace flux
15{
16namespace detail
17{
18
19//This leightweights struct and its methods finds the initial R and Z values and the coresponding f(\psi) as
20//good as it can, i.e. until machine precision is reached
21//Note that f(psi) = 1/q(psi) (The safety factor)
22struct Fpsi
23{
24
25 //firstline = 0 -> conformal, firstline = 1 -> equalarc
26 Fpsi( const CylindricalFunctorsLvl2& psip, const CylindricalFunctorsLvl1& ipol, double x0, double y0, bool verbose = false):
27 psip_(psip), fieldRZYT_(psip, ipol, x0, y0), fieldRZtau_(psip),m_verbose(verbose)
28 {
29 //Find O-point
30 double R_O = x0, Z_O = y0;
31 dg::geo::findOpoint( psip, R_O, Z_O);
32 //define angle with respect to O-point
33 fieldRZYT_ = dg::geo::flux::FieldRZYT(psip, ipol, R_O, Z_O);
34 X_init = x0, Y_init = y0;
35 while( fabs( psip.dfx()(X_init, Y_init)) <= 1e-10 && fabs( psip.dfy()( X_init, Y_init)) <= 1e-10)
36 X_init += 1.;
37 }
38 //finds the starting points for the integration in y direction
39 void find_initial( double psi, double& R_0, double& Z_0)
40 {
41 unsigned N = 50;
42 std::array<double, 2> begin2d{ {0,0} }, end2d(begin2d), end2d_old(begin2d);
43 if(m_verbose)std::cout << "In init function\n";
44 begin2d[0] = end2d[0] = end2d_old[0] = X_init;
45 begin2d[1] = end2d[1] = end2d_old[1] = Y_init;
46 double eps = 1e10, eps_old = 2e10;
47 using Vec = std::array<double,2>;
48 dg::SinglestepTimeloop<Vec> odeint( dg::RungeKutta<Vec>( "Feagin-17-8-10",
49 begin2d), fieldRZtau_);
50 while( (eps < eps_old || eps > 1e-7) && eps > 1e-14)
51 {
52 eps_old = eps; end2d_old = end2d;
53 N*=2; odeint.integrate_steps( psip_.f()(X_init, Y_init), begin2d,
54 psi, end2d, N);
55 eps = sqrt( (end2d[0]-end2d_old[0])*(end2d[0]-end2d_old[0]) + (end2d[1]-end2d_old[1])*(end2d[1]-end2d_old[1]));
56 }
57 X_init = R_0 = end2d_old[0], Y_init = Z_0 = end2d_old[1];
58 if(m_verbose)std::cout << "In init function error: psi(R,Z)-psi0: "<<psip_.f()(X_init, Y_init)-psi<<"\n";
59 }
60
61 //compute f for a given psi between psi0 and psi1
62 double construct_f( double psi, double& R_0, double& Z_0)
63 {
64 find_initial( psi, R_0, Z_0);
65 std::array<double,3> begin{ {0,0,0} }, end(begin), end_old(begin);
66 begin[0] = R_0, begin[1] = Z_0;
67 double eps = 1e10, eps_old = 2e10;
68 unsigned N = 50;
69 unsigned nan_counter = 0;
70 using Vec = std::array<double,3>;
71 dg::SinglestepTimeloop<Vec> odeint( dg::RungeKutta<Vec>( "Feagin-17-8-10",
72 begin), fieldRZYT_);
73 while( (eps < eps_old || eps > 1e-7)&& eps > 1e-14)
74 {
75 eps_old = eps, end_old = end; N*=2;
76 odeint.integrate_steps( 0., begin, 2*M_PI, end, N);
77 eps = sqrt( (end[0]-begin[0])*(end[0]-begin[0]) +
78 (end[1]-begin[1])*(end[1]-begin[1]));
79 if(m_verbose)std::cout << "\t error "<<eps<<" with "<<N<<" steps\t";
80 if( std::isnan( eps) && nan_counter < 4) eps = 1e10, end = end_old, nan_counter++;
81 }
82 if(m_verbose)std::cout << "\t error "<<eps<<" with "<<N<<" steps\t";
83 if(m_verbose)std::cout <<end_old[2] << " "<<end[2] <<"\n";
84 double f_psi = 2.*M_PI/end[2]; //this actually is 1/q the safety factor
85 return f_psi;
86 }
87
88 double operator()( double psi)
89 {
90 double R_0, Z_0;
91 return construct_f( psi, R_0, Z_0);
92 }
93 double f_prime( double psi)
94 {
95 //compute fprime
96 double deltaPsi = fabs(psi)/100.;
97 double fofpsi[4];
98 fofpsi[1] = operator()(psi-deltaPsi);
99 fofpsi[2] = operator()(psi+deltaPsi);
100 double fprime = (-0.5*fofpsi[1]+0.5*fofpsi[2])/deltaPsi, fprime_old;
101 double eps = 1e10, eps_old=2e10;
102 while( eps < eps_old)
103 {
104 deltaPsi /=2.;
105 fprime_old = fprime;
106 eps_old = eps;
107 fofpsi[0] = fofpsi[1], fofpsi[3] = fofpsi[2];
108 fofpsi[1] = operator()(psi-deltaPsi);
109 fofpsi[2] = operator()(psi+deltaPsi);
110 //reuse previously computed fpsi for current fprime
111 fprime = (+ 1./12.*fofpsi[0]
112 - 2./3. *fofpsi[1]
113 + 2./3. *fofpsi[2]
114 - 1./12.*fofpsi[3]
115 )/deltaPsi;
116 eps = fabs((fprime - fprime_old)/fprime);
117 //std::cout << "fprime "<<fprime<<" rel error fprime is "<<eps<<" delta psi "<<deltaPsi<<"\n";
118 }
119 return fprime_old;
120 }
121
122 private:
123 double X_init, Y_init;
124 CylindricalFunctorsLvl1 psip_;
125 dg::geo::flux::FieldRZYT fieldRZYT_;
126 dg::geo::FieldRZtau fieldRZtau_;
127 bool m_verbose;
128
129};
130
131} //namespace detail
132}//namespace flux
134
144struct FluxGenerator : public aGenerator2d
145{
159 FluxGenerator( const CylindricalFunctorsLvl2& psi, const CylindricalFunctorsLvl1& ipol, double psi_0, double psi_1, double x0, double y0, int mode=0, bool verbose = false):
160 psi_(psi), ipol_(ipol), mode_(mode), m_verbose( verbose)
161 {
162 psi0_ = psi_0, psi1_ = psi_1;
163 assert( psi_1 != psi_0);
164 if( mode==0)
165 {
166 flux::detail::Fpsi fpsi(psi, ipol, x0, y0, m_verbose);
167 f0_ = fabs( fpsi.construct_f( psi_0, x0_, y0_));
168 }
169 else
170 {
171 ribeiro::detail::Fpsi fpsi(psi, x0, y0, mode, m_verbose);
172 f0_ = fabs( fpsi.construct_f( psi_0, x0_, y0_));
173 }
174 if( psi_1 < psi_0) f0_*=-1;
175 lx_ = f0_*(psi_1-psi_0);
176 x0_=x0, y0_=y0, psi0_=psi_0, psi1_=psi_1;
177 if(m_verbose)std::cout << "lx = "<<lx_<<"\n";
178 }
179
180 virtual FluxGenerator* clone() const override final{return new FluxGenerator(*this);}
181
182 private:
183 // length of zeta-domain (f0*(psi_1-psi_0))
184 virtual double do_width() const override final{return lx_;}
185 virtual double do_height() const override final{return 2.*M_PI;}
186 virtual void do_generate(
187 const thrust::host_vector<double>& zeta1d,
188 const thrust::host_vector<double>& eta1d,
189 thrust::host_vector<double>& x,
190 thrust::host_vector<double>& y,
191 thrust::host_vector<double>& zetaX,
192 thrust::host_vector<double>& zetaY,
193 thrust::host_vector<double>& etaX,
194 thrust::host_vector<double>& etaY) const override final
195 {
196 //compute psi(x) for a grid on x and call construct_rzy for all psi
197 thrust::host_vector<double> psi_x(zeta1d);
198 for( unsigned i=0; i<psi_x.size(); i++)
199 psi_x[i] = zeta1d[i]/f0_ +psi0_;
200
201 if(m_verbose)std::cout << "In grid function:"<<std::endl;
202 flux::detail::Fpsi fpsi(psi_, ipol_, x0_, y0_, m_verbose);
203 dg::geo::flux::FieldRZYRYZY fieldRZYRYZY(psi_, ipol_);
204 ribeiro::detail::Fpsi fpsiRibeiro(psi_, x0_, y0_, mode_, m_verbose);
205 dg::geo::equalarc::FieldRZYRYZY fieldRZYRYZYequalarc(psi_);
206 thrust::host_vector<double> fx_;
207 fx_.resize( zeta1d.size());
208 thrust::host_vector<double> f_p(fx_);
209 unsigned Nx = zeta1d.size(), Ny = eta1d.size();
210 for( unsigned i=0; i<zeta1d.size(); i++)
211 {
212 thrust::host_vector<double> ry, zy;
213 thrust::host_vector<double> yr, yz, xr, xz;
214 double R0, Z0;
215 if(mode_==0)dg::geo::detail::compute_rzy( fpsi, fieldRZYRYZY, psi_x[i], eta1d, ry, zy, yr, yz, xr, xz, R0, Z0, fx_[i], f_p[i], m_verbose);
216 if(mode_==1)dg::geo::detail::compute_rzy( fpsiRibeiro, fieldRZYRYZYequalarc, psi_x[i], eta1d, ry, zy, yr, yz, xr, xz, R0, Z0, fx_[i], f_p[i], m_verbose);
217 for( unsigned j=0; j<Ny; j++)
218 {
219 x[j*Nx+i] = ry[j], y[j*Nx+i] = zy[j];
220 etaX[j*Nx+i] = yr[j], etaY[j*Nx+i] = yz[j];
221 zetaX[j*Nx+i] = xr[j]/fx_[i]*f0_, zetaY[j*Nx+i] = xz[j]/fx_[i]*f0_;
222 }
223 }
224 }
225 CylindricalFunctorsLvl2 psi_;
226 CylindricalFunctorsLvl1 ipol_;
227 double f0_, lx_, x0_, y0_, psi0_, psi1_;
228 int mode_;
229 bool m_verbose;
230};
231
237struct RibeiroFluxGenerator : public aGenerator2d
238{
250 RibeiroFluxGenerator( const CylindricalFunctorsLvl2& psi, double psi_0, double psi_1, double x0, double y0, int mode=0, bool verbose = false):
251 psip_(psi), mode_(mode), m_verbose(verbose)
252 {
253 psi0_ = psi_0, psi1_ = psi_1;
254 assert( psi_1 != psi_0);
255 ribeiro::detail::Fpsi fpsi(psi, x0, y0, mode, m_verbose);
256 f0_ = fabs( fpsi.construct_f( psi_0, x0_, y0_));
257 if( psi_1 < psi_0) f0_*=-1;
258 lx_ = f0_*(psi_1-psi_0);
259 x0_=x0, y0_=y0, psi0_=psi_0, psi1_=psi_1;
260 if(m_verbose)std::cout << "lx = "<<lx_<<"\n";
261 }
262 virtual RibeiroFluxGenerator* clone() const{return new RibeiroFluxGenerator(*this);}
263
264 private:
265 //length of zeta-domain (f0*(psi_1-psi_0))
266 virtual double do_width() const{return lx_;}
267 virtual double do_height() const{return 2.*M_PI;}
268 virtual void do_generate(
269 const thrust::host_vector<double>& zeta1d,
270 const thrust::host_vector<double>& eta1d,
271 thrust::host_vector<double>& x,
272 thrust::host_vector<double>& y,
273 thrust::host_vector<double>& zetaX,
274 thrust::host_vector<double>& zetaY,
275 thrust::host_vector<double>& etaX,
276 thrust::host_vector<double>& etaY) const
277 {
278 //compute psi(x) for a grid on x and call construct_rzy for all psi
279 thrust::host_vector<double> psi_x(zeta1d);
280 for( unsigned i=0; i<psi_x.size(); i++)
281 psi_x[i] = zeta1d[i]/f0_ +psi0_;
282
283 ribeiro::detail::Fpsi fpsi(psip_, x0_, y0_, mode_, m_verbose);
284 dg::geo::ribeiro::FieldRZYRYZY fieldRZYRYZYribeiro(psip_);
285 dg::geo::equalarc::FieldRZYRYZY fieldRZYRYZYequalarc(psip_);
286 thrust::host_vector<double> fx_;
287 fx_.resize( zeta1d.size());
288 thrust::host_vector<double> f_p(fx_);
289 unsigned Nx = zeta1d.size(), Ny = eta1d.size();
290 for( unsigned i=0; i<zeta1d.size(); i++)
291 {
292 thrust::host_vector<double> ry, zy;
293 thrust::host_vector<double> yr, yz, xr, xz;
294 double R0, Z0;
295 if(mode_==0)dg::geo::detail::compute_rzy( fpsi, fieldRZYRYZYribeiro, psi_x[i], eta1d, ry, zy, yr, yz, xr, xz, R0, Z0, fx_[i], f_p[i], m_verbose);
296 if(mode_==1)dg::geo::detail::compute_rzy( fpsi, fieldRZYRYZYequalarc, psi_x[i], eta1d, ry, zy, yr, yz, xr, xz, R0, Z0, fx_[i], f_p[i], m_verbose);
297 for( unsigned j=0; j<Ny; j++)
298 {
299 x[j*Nx+i] = ry[j], y[j*Nx+i] = zy[j];
300 etaX[j*Nx+i] = yr[j], etaY[j*Nx+i] = yz[j];
301 zetaX[j*Nx+i] = xr[j]/fx_[i]*f0_, zetaY[j*Nx+i] = xz[j]/fx_[i]*f0_;
302 }
303 }
304 }
305 CylindricalFunctorsLvl2 psip_;
306 double f0_, lx_, x0_, y0_, psi0_, psi1_;
307 int mode_;
308 bool m_verbose;
309};
310}//namespace geo
311}//namespace dg
M_PI
#define M_PI
coo3d::yz
@ yz
coo3d::xz
@ xz
coo2d::y
@ y
coo2d::x
@ x
dg::geo::findOpoint
static int findOpoint(const CylindricalFunctorsLvl2 &psi, double &RC, double &ZC)
This function finds O-points of psi.
Definition: fluxfunctions.h:332
dg::geo::CylindricalFunctorsLvl1
This struct bundles a function and its first derivatives.
Definition: fluxfunctions.h:182
dg::geo::CylindricalFunctorsLvl2
This struct bundles a function and its first and second derivatives.
Definition: fluxfunctions.h:219
dg::geo::FluxGenerator
A symmetry flux generator.
Definition: flux.h:145
dg::geo::FluxGenerator::clone
virtual FluxGenerator * clone() const override final
Abstract clone method that returns a copy on the heap.
Definition: flux.h:180
dg::geo::FluxGenerator::FluxGenerator
FluxGenerator(const CylindricalFunctorsLvl2 &psi, const CylindricalFunctorsLvl1 &ipol, double psi_0, double psi_1, double x0, double y0, int mode=0, bool verbose=false)
Construct a symmetry flux grid generator.
Definition: flux.h:159
dg::geo::RibeiroFluxGenerator
Same as the Ribeiro class just but uses psi as a flux label directly.
Definition: flux.h:238
dg::geo::RibeiroFluxGenerator::clone
virtual RibeiroFluxGenerator * clone() const
Abstract clone method that returns a copy on the heap.
Definition: flux.h:262
dg::geo::RibeiroFluxGenerator::RibeiroFluxGenerator
RibeiroFluxGenerator(const CylindricalFunctorsLvl2 &psi, double psi_0, double psi_1, double x0, double y0, int mode=0, bool verbose=false)
Construct a flux aligned grid generator.
Definition: flux.h:250
dg::geo::aRealGenerator2d
The abstract generator base class.
Definition: generator.h:20