Extension: Geometries
#include "dg/geometries/geometries.h"
simple_orthogonal.h
Go to the documentation of this file.
1#pragma once
2
3#include "dg/algorithm.h"
4#include "generator.h"
5#include "utilities.h"
6#include "adaption.h"
7
8
9namespace dg
10{
11namespace geo
12{
14namespace orthogonal
15{
16
17namespace detail
18{
19
20//This leightweights struct and its methods finds the initial R and Z values and the coresponding f(\psi) as
21//good as it can, i.e. until machine precision is reached
22struct Fpsi
23{
24
25 //firstline = 0 -> conformal, firstline = 1 -> equalarc
26 Fpsi( const CylindricalFunctorsLvl1& psi, const CylindricalSymmTensorLvl1& chi, double x0, double y0, int firstline):
27 psip_(psi), fieldRZYTconf_(psi, x0, y0, chi),fieldRZYTequl_(psi, x0, y0, chi), fieldRZtau_(psi, chi)
28 {
29 X_init = x0, Y_init = y0;
30 while( fabs( psi.dfx()(X_init, Y_init)) <= 1e-10 && fabs( psi.dfy()( X_init, Y_init)) <= 1e-10)
31 X_init += 1.;
32 m_firstline = firstline;
33 }
34 //finds the starting points for the integration in y direction
35 void find_initial( double psi, double& R_0, double& Z_0)
36 {
37 unsigned N = 50;
38 std::array<double, 2> begin2d{ {0,0} }, end2d(begin2d), end2d_old(begin2d);
39 begin2d[0] = end2d[0] = end2d_old[0] = X_init;
40 begin2d[1] = end2d[1] = end2d_old[1] = Y_init;
41 double eps = 1e10, eps_old = 2e10;
42 while( (eps < eps_old || eps > 1e-7) && eps > 1e-14)
43 {
44 eps_old = eps; end2d_old = end2d;
45 N*=2;
46 double psi0 = psip_.f()(X_init, Y_init);
47 using Vec = std::array<double,2>;
48 dg::SinglestepTimeloop<Vec> odeint( dg::RungeKutta<Vec>(
49 "Feagin-17-8-10", {0,0}), fieldRZtau_);
50 odeint.integrate_steps( psi0, begin2d, psi, end2d, N);
51 eps = sqrt( (end2d[0]-end2d_old[0])*(end2d[0]-end2d_old[0]) +
52 (end2d[1]-end2d_old[1])*(end2d[1]-end2d_old[1]));
53 }
54 X_init = R_0 = end2d_old[0], Y_init = Z_0 = end2d_old[1];
55 }
56
57 //compute f for a given psi between psi0 and psi1
58 double construct_f( double psi, double& R_0, double& Z_0)
59 {
60 find_initial( psi, R_0, Z_0);
61 std::array<double, 3> begin{ {0,0,0} }, end(begin);
62 begin[0] = R_0, begin[1] = Z_0;
63 double eps = 1e10, eps_old = 2e10;
64 unsigned N = 50;
65 while( (eps < eps_old || eps > 1e-7)&& eps > 1e-14)
66 {
67 eps_old = eps; N*=2;
68 using Vec = std::array<double,3>;
69 dg::RungeKutta<Vec> rk( "Feagin-17-8-10", begin);
70 dg::SinglestepTimeloop<Vec> odeint;
71 if( m_firstline == 0)
72 odeint = dg::SinglestepTimeloop<Vec>( rk,
73 fieldRZYTconf_);
74 if( m_firstline == 1)
75 odeint = dg::SinglestepTimeloop<Vec>( rk,
76 fieldRZYTequl_);
77 odeint.integrate_steps( 0., begin, 2.*M_PI, end, N);
78 eps = sqrt( (end[0]-begin[0])*(end[0]-begin[0]) + (end[1]-begin[1])*(end[1]-begin[1]));
79 }
80 //std::cout << "\t error "<<eps<<" with "<<N<<" steps\t";
81 double f_psi = 2.*M_PI/end[2];
82 //std::cout <<"f_psi: "<<f_psi << " end "<<end[2] <<"\n";
83 return f_psi;
84 }
85 double operator()( double psi)
86 {
87 double R_0, Z_0;
88 return construct_f( psi, R_0, Z_0);
89 }
90
91 private:
92 int m_firstline;
93 double X_init, Y_init;
94 CylindricalFunctorsLvl1 psip_;
95 CylindricalSymmTensorLvl1 chi_;
96 dg::geo::ribeiro::FieldRZYT fieldRZYTconf_;
97 dg::geo::equalarc::FieldRZYT fieldRZYTequl_;
98 dg::geo::FieldRZtau fieldRZtau_;
99
100};
101
102//compute the vector of r and z - values that form one psi surface
103//assumes y_0 = 0
104template<class real_type>
105void compute_rzy( const CylindricalFunctorsLvl1& psi, const CylindricalSymmTensorLvl1& chi,
106 const thrust::host_vector<real_type>& y_vec,
107 thrust::host_vector<real_type>& r,
108 thrust::host_vector<real_type>& z,
109 real_type R_0, real_type Z_0, real_type f_psi, int mode )
110{
111
112 thrust::host_vector<real_type> r_old(y_vec.size(), 0), r_diff( r_old);
113 thrust::host_vector<real_type> z_old(y_vec.size(), 0), z_diff( z_old);
114 r.resize( y_vec.size()), z.resize(y_vec.size());
115 std::array<real_type,2> begin{ {0,0} }, end(begin), temp(begin);
116 begin[0] = R_0, begin[1] = Z_0;
117 //std::cout <<"f_psi "<<f_psi<<" "<<" "<< begin[0] << " "<<begin[1]<<"\t";
118 dg::geo::ribeiro::FieldRZY fieldRZYconf(psi, chi);
119 dg::geo::equalarc::FieldRZY fieldRZYequi(psi, chi);
120 fieldRZYconf.set_f(f_psi);
121 fieldRZYequi.set_f(f_psi);
122 unsigned steps = 1;
123 real_type eps = 1e10, eps_old=2e10;
124 while( (eps < eps_old||eps > 1e-7) && eps > 1e-14)
125 {
126 //begin is left const
127 eps_old = eps, r_old = r, z_old = z;
128 using Vec = std::array<real_type,2>;
129 dg::RungeKutta<Vec> rk ( "Feagin-17-8-10", {0,0});
130 dg::SinglestepTimeloop<Vec> odeint;
131 if( mode == 0)
132 odeint = dg::SinglestepTimeloop<Vec>( rk, fieldRZYconf);
133 if( mode == 1)
134 odeint = dg::SinglestepTimeloop<Vec>( rk, fieldRZYequi);
135 odeint.integrate_steps( 0., begin, y_vec[0], end, steps);
136 r[0] = end[0], z[0] = end[1];
137 for( unsigned i=1; i<y_vec.size(); i++)
138 {
139 temp = end;
140 odeint.integrate_steps( y_vec[i-1], temp, y_vec[i], end, steps);
141 r[i] = end[0], z[i] = end[1];
142 }
143 //compute error in R,Z only
144 dg::blas1::axpby( 1., r, -1., r_old, r_diff);
145 dg::blas1::axpby( 1., z, -1., z_old, z_diff);
146 real_type er = dg::blas1::dot( r_diff, r_diff);
147 real_type ez = dg::blas1::dot( z_diff, z_diff);
148 real_type ar = dg::blas1::dot( r, r);
149 real_type az = dg::blas1::dot( z, z);
150 eps = sqrt( er + ez)/sqrt(ar+az);
151 steps*=2;
152 }
153 r = r_old, z = z_old;
154
155}
156
157//This struct computes -2pi/f with a fixed number of steps for all psi
158//and provides the Nemov algorithm for orthogonal grid
159struct Nemov
160{
161 Nemov( const CylindricalFunctorsLvl2 psi, const CylindricalSymmTensorLvl1& chi, double f0, int mode):
162 f0_(f0), mode_(mode),
163 psip_(psi), chi_(chi), lapPsi_(psi,chi) { }
164 void initialize(
165 const thrust::host_vector<double>& r_init, //1d intial values
166 const thrust::host_vector<double>& z_init, //1d intial values
167 thrust::host_vector<double>& h_init) //,
168 // thrust::host_vector<double>& hr_init,
169 // thrust::host_vector<double>& hz_init)
170 {
171 unsigned size = r_init.size();
172 h_init.resize( size);//, hr_init.resize( size), hz_init.resize( size);
173 for( unsigned i=0; i<size; i++)
174 {
175 if(mode_ == 0)
176 h_init[i] = f0_;
177 if(mode_ == 1)
178 {
179 double x = r_init[i], y = z_init[i];
180 double psipR = psip_.dfx()(x, y), psipZ = psip_.dfy()(x,y);
181 double chiRR = chi_.xx()(x, y),
182 chiRZ = chi_.xy()(x, y),
183 chiZZ = chi_.yy()(x, y);
184 double psip2 = chiRR*psipR*psipR + 2.*chiRZ*psipR*psipZ + chiZZ*psipZ*psipZ;
185 h_init[i] = f0_/sqrt(psip2); //equalarc
186 }
187 //double laplace = psipRR_(r_init[i], z_init[i]) +
188 //psipZZ_(r_init[i], z_init[i]);
189 //hr_init[i] = -f0_*laplace/psip2*psipR;
190 //hz_init[i] = -f0_*laplace/psip2*psipZ;
191 }
192 }
193
194 void operator()(double t, const std::array<thrust::host_vector<double>,3 >& y, std::array<thrust::host_vector<double>,3>& yp)
195 {
196 //y[0] = R, y[1] = Z, y[2] = h, y[3] = hr, y[4] = hz
197 unsigned size = y[0].size();
198 for( unsigned i=0; i<size; i++)
199 {
200 double xx = y[0][i], yy = y[1][i];
201 double psipR = psip_.dfx()(xx, yy), psipZ = psip_.dfy()(xx,yy);
202 double chiRR = chi_.xx()(xx, yy),
203 chiRZ = chi_.xy()(xx, yy),
204 chiZZ = chi_.yy()(xx, yy);
205 double psip2 = chiRR*psipR*psipR + 2.*chiRZ*psipR*psipZ + chiZZ*psipZ*psipZ;
206 yp[0][i] = (chiRR*psipR + chiRZ*psipZ)/psip2/f0_;
207 yp[1][i] = (chiRZ*psipR + chiZZ*psipZ)/psip2/f0_;
208 yp[2][i] = y[2][i]*( - lapPsi_(xx,yy) )/psip2/f0_;
209 //yp[3][i] = ( -(2.*psipRR+psipZZ)*y[3][i] - psipRZ*y[4][i] - laplacePsipR_(y[0][i], y[1][i])*y[2][i])/psip2; //wrong with monitor metric!!
210 //yp[4][i] = ( -psipRZ*y[3][i] - (2.*psipZZ+psipRR)*y[4][i] - laplacePsipZ_(y[0][i], y[1][i])*y[2][i])/psip2; //wrong with monitor metric!!
211 }
212 }
213 private:
214 double f0_;
215 int mode_;
216 CylindricalFunctorsLvl2 psip_;
217 CylindricalSymmTensorLvl1 chi_;
218 dg::geo::detail::LaplaceChiPsi lapPsi_;
219};
220
221template<class Nemov>
222void construct_rz( Nemov nemov,
223 double x_0, //the x value that corresponds to the first psi surface
224 const thrust::host_vector<double>& x_vec, //1d x values
225 const thrust::host_vector<double>& r_init, //1d intial values of the first psi surface
226 const thrust::host_vector<double>& z_init, //1d intial values of the first psi surface
227 thrust::host_vector<double>& r,
228 thrust::host_vector<double>& z,
229 thrust::host_vector<double>& h
230 )
231{
232 unsigned N = 1;
233 double eps = 1e10, eps_old=2e10;
234 std::array<thrust::host_vector<double>,3> begin;
235 thrust::host_vector<double> h_init( r_init.size(), 0.);
236 nemov.initialize( r_init, z_init, h_init);
237 begin[0] = r_init, begin[1] = z_init, begin[2] = h_init;
238 //now we have the starting values
239 std::array<thrust::host_vector<double>,3> end(begin), temp(begin);
240 unsigned sizeX = x_vec.size(), sizeY = r_init.size();
241 unsigned size2d = x_vec.size()*r_init.size();
242 r.resize(size2d), z.resize(size2d), h.resize(size2d);
243 double x0=x_0, x1 = x_vec[0];
244 thrust::host_vector<double> r_new(r_init), r_old(r_init), r_diff(r_init);
245 thrust::host_vector<double> z_new(z_init), z_old(z_init), z_diff(z_init);
246 thrust::host_vector<double> h_new(h_init); //, h_old(h_init), h_diff(h_init);
247 for( unsigned i=0; i<sizeX; i++)
248 {
249 N = 1;
250 eps = 1e10, eps_old=2e10;
251 begin = temp;
252 while( (eps < eps_old || eps > 1e-6) && eps > 1e-13)
253 {
254 r_old = r_new, z_old = z_new; eps_old = eps;
255 //h_old = h_new;
256 temp = begin;
258 x0 = i==0?x_0:x_vec[i-1], x1 = x_vec[i];
260 using Vec = std::array<thrust::host_vector<double>,3>;
261 dg::RungeKutta<Vec> rk( "Feagin-17-8-10", temp);
262 dg::SinglestepTimeloop<Vec>( rk, nemov).integrate_steps( x0, temp,
263 x1, end, N);
264 for( unsigned j=0; j<sizeY; j++)
265 {
266 r_new[j] = end[0][j], z_new[j] = end[1][j];
267 h_new[j] = end[2][j];
268 }
270 temp = end;
271 dg::blas1::axpby( 1., r_new, -1., r_old, r_diff);
272 dg::blas1::axpby( 1., z_new, -1., z_old, z_diff);
273 //dg::blas1::axpby( 1., h_new, -1., h_old, h_diff);
274 //dg::blas1::pointwiseDot( h_diff, h_diff, h_diff);
275 //dg::blas1::pointwiseDivide( h_diff, h_old, h_diff); // h is always > 0
276 dg::blas1::pointwiseDot( r_diff, r_diff, r_diff);
277 dg::blas1::pointwiseDot( 1., z_diff, z_diff, 1., r_diff);
278 //dg::blas1::axpby( 1., h_diff, 1., r_diff);
279 try{
280 eps = sqrt( dg::blas1::dot( r_diff, 1.)/sizeY); //should be relative to the interpoint distances
281 //double eps_h = sqrt( dg::blas1::dot( h_diff, 1.)/sizeY);
282 //std::cout << "Effective Relative diff-h error is "<<eps_h<<" with "<<N<<" steps\n";
283 } catch ( dg::Error& )
284 {
285 eps = eps_old;
286 r_new = r_old , z_new = z_old;
287 //h_new = h_old;
288 }
289 //std::cout << "Effective Absolute diff-r error is "<<eps<<" with "<<N<<" steps\n";
290 N*=2;
291 if( (eps > eps_old && N > 1024 && eps > 1e-6) || N > 64000)
292 throw dg::Error(dg::Message(_ping_) <<
293 "Grid generator encountered loss of convergence integrating from x = "
294 <<x0<<" to x = "<<x1);
295 }
296 for( unsigned j=0; j<sizeY; j++)
297 {
298 unsigned idx = sizeX*j+i;
299 r[idx] = r_new[j], z[idx] = z_new[j], h[idx] = h_new[j];
300 }
301 }
302
303}
304
305} //namespace detail
306
307}//namespace orthogonal
309
310
327struct SimpleOrthogonal : public aGenerator2d
328{
342 SimpleOrthogonal(const CylindricalFunctorsLvl2& psi, double psi_0, double
343 psi_1, double x0, double y0, double psi_firstline, int mode =0
344 ):
345 SimpleOrthogonal( psi, CylindricalSymmTensorLvl1(), psi_0, psi_1, x0, y0,
346 psi_firstline, mode)
347 {
348 m_orthogonal = true;
349 }
364 SimpleOrthogonal(const CylindricalFunctorsLvl2& psi, const
365 CylindricalSymmTensorLvl1& chi, double psi_0, double psi_1,
366 double x0, double y0, double psi_firstline, int mode = 0
367 ):
368 psi_(psi), chi_(chi)
369 {
370 assert( psi_1 != psi_0);
371 m_firstline = mode;
372 orthogonal::detail::Fpsi fpsi(psi, chi, x0, y0, mode);
373 f0_ = fabs( fpsi.construct_f( psi_firstline, R0_, Z0_));
374 if( psi_1 < psi_0) f0_*=-1;
375 lz_ = f0_*(psi_1-psi_0);
376 m_orthogonal = false;
377 m_zeta_first = f0_*(psi_firstline-psi_0);
378 }
379
385 double f0() const{return f0_;}
386 virtual SimpleOrthogonal* clone() const override final{return new SimpleOrthogonal(*this);}
387
388 private:
389 // length of zeta-domain (f0*(psi_1-psi_0))
390 virtual double do_width() const override final{return lz_;}
391 virtual double do_height() const override final{return 2.*M_PI;}
392 virtual bool do_isOrthogonal() const override final{return m_orthogonal;}
393 virtual void do_generate(
394 const thrust::host_vector<double>& zeta1d,
395 const thrust::host_vector<double>& eta1d,
396 thrust::host_vector<double>& x,
397 thrust::host_vector<double>& y,
398 thrust::host_vector<double>& zetaX,
399 thrust::host_vector<double>& zetaY,
400 thrust::host_vector<double>& etaX,
401 thrust::host_vector<double>& etaY) const override final
402 {
403 thrust::host_vector<double> r_init, z_init;
404 orthogonal::detail::compute_rzy( psi_, chi_, eta1d, r_init, z_init,
405 R0_, Z0_, f0_, m_firstline);
406 orthogonal::detail::Nemov nemov(psi_, chi_, f0_, m_firstline);
407 thrust::host_vector<double> h;
408 orthogonal::detail::construct_rz(nemov, m_zeta_first, zeta1d, r_init,
409 z_init, x, y, h);
410 unsigned size = x.size();
411 for( unsigned idx=0; idx<size; idx++)
412 {
413 double psipR = psi_.dfx()(x[idx], y[idx]);
414 double psipZ = psi_.dfy()(x[idx], y[idx]);
415 double chiRR = chi_.xx()( x[idx], y[idx]),
416 chiRZ = chi_.xy()( x[idx], y[idx]),
417 chiZZ = chi_.yy()( x[idx], y[idx]);
418 zetaX[idx] = f0_*psipR;
419 zetaY[idx] = f0_*psipZ;
420 etaX[idx] = -h[idx]*(chiRZ*psipR + chiZZ*psipZ);
421 etaY[idx] = +h[idx]*(chiRR*psipR + chiRZ*psipZ);
422 }
423 }
424 CylindricalFunctorsLvl2 psi_;
425 CylindricalSymmTensorLvl1 chi_;
426 double f0_, lz_, R0_, Z0_;
427 int m_firstline;
428 double m_zeta_first;
429 bool m_orthogonal;
430};
431
432}//namespace geo
433}//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::blas1::pointwiseDot
void pointwiseDot(get_value_type< ContainerType > alpha, const ContainerType1 &x1, const ContainerType2 &x2, get_value_type< ContainerType > beta, ContainerType &y)
dg::blas1::dot
get_value_type< ContainerType1 > dot(const ContainerType1 &x, const ContainerType2 &y)
coo3d::z
@ z
coo2d::y
@ y
coo2d::x
@ x
dg::SinglestepTimeloop::integrate_steps
void integrate_steps(value_type t0, const container_type &u0, value_type t1, container_type &u1, unsigned steps)
dg::geo::CylindricalFunctorsLvl2
This struct bundles a function and its first and second derivatives.
Definition: fluxfunctions.h:219
dg::geo::CylindricalFunctorsLvl2::dfy
const CylindricalFunctor & dfy() const
Definition: fluxfunctions.h:247
dg::geo::CylindricalFunctorsLvl2::dfx
const CylindricalFunctor & dfx() const
Definition: fluxfunctions.h:245
dg::geo::CylindricalSymmTensorLvl1
Definition: fluxfunctions.h:361
dg::geo::CylindricalSymmTensorLvl1::yy
const CylindricalFunctor & yy() const
yy component
Definition: fluxfunctions.h:400
dg::geo::CylindricalSymmTensorLvl1::xy
const CylindricalFunctor & xy() const
xy component
Definition: fluxfunctions.h:398
dg::geo::CylindricalSymmTensorLvl1::xx
const CylindricalFunctor & xx() const
xy component
Definition: fluxfunctions.h:396
dg::geo::SimpleOrthogonal
Generate a simple orthogonal grid.
Definition: simple_orthogonal.h:328
dg::geo::SimpleOrthogonal::SimpleOrthogonal
SimpleOrthogonal(const CylindricalFunctorsLvl2 &psi, const CylindricalSymmTensorLvl1 &chi, double psi_0, double psi_1, double x0, double y0, double psi_firstline, int mode=0)
Construct a simple orthogonal grid.
Definition: simple_orthogonal.h:364
dg::geo::SimpleOrthogonal::f0
double f0() const
The grid constant.
Definition: simple_orthogonal.h:385
dg::geo::SimpleOrthogonal::SimpleOrthogonal
SimpleOrthogonal(const CylindricalFunctorsLvl2 &psi, double psi_0, double psi_1, double x0, double y0, double psi_firstline, int mode=0)
Construct a simple orthogonal grid.
Definition: simple_orthogonal.h:342
dg::geo::SimpleOrthogonal::clone
virtual SimpleOrthogonal * clone() const override final
Abstract clone method that returns a copy on the heap.
Definition: simple_orthogonal.h:386
dg::geo::aRealGenerator2d
The abstract generator base class.
Definition: generator.h:20