Extension: Geometries
#include "dg/geometries/geometries.h"
Loading...
Searching...
No Matches
fluxfunctions.h
Go to the documentation of this file.
1#pragma once
2#include <functional>
3#include "dg/algorithm.h"
4
5namespace dg
6{
7namespace geo
8{
9
12
24template<class real_type>
26{
34 template<class BinaryFunctor>
35 RealCylindricalFunctor( BinaryFunctor f):
36 m_f(f) {}
38 real_type operator()( real_type R, real_type Z) const{
39 return m_f(R,Z);
40 }
42 real_type operator()( real_type R, real_type Z, real_type phi) const{
43 return m_f(R,Z);
44 }
45 private:
46 std::function<real_type(real_type,real_type)> m_f;
47};
48
51
64template<class Derived>
66{
75 double operator()(double R, double Z) const
76 {
77 const Derived& underlying = static_cast<const Derived&>(*this);
78 return underlying.do_compute(R,Z);
79 }
89 double operator()(double R, double Z, double phi)const
90 {
91 const Derived& underlying = static_cast<const Derived&>(*this);
92 return underlying.do_compute(R,Z);
93 }
94#ifndef __CUDACC__ //nvcc below 10 has problems with the following construct
95 //This trick avoids that classes inherit from the wrong Base:
96 private:
97 friend Derived;
102 aCylindricalFunctor(const aCylindricalFunctor&){}
106 aCylindricalFunctor& operator=(const aCylindricalFunctor&){return *this;}
107#endif //__CUDACC__
108};
109
113struct Constant: public aCylindricalFunctor<Constant>
114{
115 Constant(double c):c_(c){}
116 double do_compute(double R,double Z)const{return c_;}
117 private:
118 double c_;
119};
129struct ZCutter : public aCylindricalFunctor<ZCutter>
130{
131 ZCutter(double ZX, int sign = +1): m_heavi( ZX, sign){}
132 double do_compute(double R, double Z) const {
133 return m_heavi(Z);
134 }
135 private:
136 dg::Heaviside m_heavi;
137};
138
143struct Periodify : public aCylindricalFunctor<Periodify>
144{
151 Periodify( CylindricalFunctor functor, dg::Grid2d g): m_g( g), m_f(functor) {}
163 Periodify( CylindricalFunctor functor, double R0, double R1, double Z0,
164 double Z1, dg::bc bcx, dg::bc bcy):
165 m_g( R0, R1, Z0, Z1, 3, 10, 10, bcx, bcy), m_f(functor)
166 {}
167 double do_compute( double R, double Z) const
168 {
169 bool negative = false;
170 dg::create::detail::shift( negative, R, m_g.bcx(), m_g.x0(), m_g.x1());
171 dg::create::detail::shift( negative, Z, m_g.bcy(), m_g.y0(), m_g.y1());
172 if( negative) return -m_f(R,Z);
173 return m_f( R, Z);
174 }
175 private:
176 dg::Grid2d m_g;
178};
179
186{
201 {
202 p_[0] = f;
203 p_[1] = fx;
204 p_[2] = fy;
205 }
207 const CylindricalFunctor& f()const{return p_[0];}
209 const CylindricalFunctor& dfx()const{return p_[1];}
211 const CylindricalFunctor& dfy()const{return p_[2];}
212 private:
213 std::array<CylindricalFunctor,3> p_;
214};
215
216
223{
241 {
242 f0.reset( f,fx,fy), f1.reset(fxx,fxy,fyy);
243 }
245 operator CylindricalFunctorsLvl1 ()const {return f0;}
247 const CylindricalFunctor& f()const{return f0.f();}
249 const CylindricalFunctor& dfx()const{return f0.dfx();}
251 const CylindricalFunctor& dfy()const{return f0.dfy();}
253 const CylindricalFunctor& dfxx()const{return f1.f();}
255 const CylindricalFunctor& dfxy()const{return f1.dfx();}
257 const CylindricalFunctor& dfyy()const{return f1.dfy();}
258 private:
260};
261
262
279inline int findCriticalPoint( const CylindricalFunctorsLvl2& psi, double& RC, double& ZC)
280{
281 std::array<double, 2> X{ {0,0} }, XN(X), X_OLD(X);
282 X[0] = RC, X[1] = ZC;
283 double eps = 1e10, eps_old= 2e10;
284 unsigned counter = 0; //safety measure to avoid deadlock
285 double psipRZ = psi.dfxy()(X[0], X[1]);
286 double psipRR = psi.dfxx()(X[0], X[1]), psipZZ = psi.dfyy()(X[0],X[1]);
287 double psipR = psi.dfx()(X[0], X[1]), psipZ = psi.dfy()(X[0], X[1]);
288 double D0 = (psipZZ*psipRR - psipRZ*psipRZ);
289 if(D0 == 0) // try to change initial guess slightly if we are very lucky
290 {
291 X[0] *= 1.0001, X[1]*=1.0001;
292 psipRZ = psi.dfxy()(X[0], X[1]);
293 psipRR = psi.dfxx()(X[0], X[1]), psipZZ = psi.dfyy()(X[0],X[1]);
294 psipR = psi.dfx()(X[0], X[1]), psipZ = psi.dfy()(X[0], X[1]);
295 D0 = (psipZZ*psipRR - psipRZ*psipRZ);
296 }
297 double Dinv = 1./D0;
298 while( (eps < eps_old || eps > 1e-7) && eps > 1e-10 && counter < 100)
299 {
300 //newton iteration
301 XN[0] = X[0] - Dinv*(psipZZ*psipR - psipRZ*psipZ);
302 XN[1] = X[1] - Dinv*(-psipRZ*psipR + psipRR*psipZ);
303 XN.swap(X);
304 eps = sqrt( (X[0]-X_OLD[0])*(X[0]-X_OLD[0]) + (X[1]-X_OLD[1])*(X[1]-X_OLD[1]));
305 X_OLD = X; eps_old= eps;
306 psipRZ = psi.dfxy()(X[0], X[1]);
307 psipRR = psi.dfxx()(X[0], X[1]), psipZZ = psi.dfyy()(X[0],X[1]);
308 psipR = psi.dfx()(X[0], X[1]), psipZ = psi.dfy()(X[0], X[1]);
309 D0 = (psipZZ*psipRR - psipRZ*psipRZ);
310 Dinv = 1./D0;
311 if( D0 == 0) break;
312 counter++;
313 }
314 if ( counter >= 100 || D0 == 0|| std::isnan( Dinv) )
315 return 0;
316 RC = X[0], ZC = X[1];
317 if( Dinv > 0 && psipRR > 0)
318 return 1; //local minimum
319 if( Dinv > 0 && psipRR < 0)
320 return 2; //local maximum
321 //if( Dinv < 0)
322 return 3; //saddle point
323}
324
336inline int findOpoint( const CylindricalFunctorsLvl2& psi, double& RC, double& ZC)
337{
338 int point = findCriticalPoint( psi, RC, ZC);
339 if( point == 3 || point == 0 )
340 throw dg::Error(dg::Message(_ping_)<<"There is no O-point near "<<RC<<" "<<ZC);
341 return point;
342}
343
354inline void findXpoint( const CylindricalFunctorsLvl2& psi, double& RC, double& ZC)
355{
356 int point = findCriticalPoint( psi, RC, ZC);
357 if( point != 3)
358 throw dg::Error(dg::Message(_ping_)<<"There is no X-point near "<<RC<<" "<<ZC);
359}
360
361
365{
370 reset( Constant(1), Constant(0), Constant(1), Constant(0), Constant(0));
371 }
384 CylindricalFunctor divChiX, CylindricalFunctor divChiY) :
385 p_{{ chi_xx,chi_xy,chi_yy,divChiX,divChiY}}
386 {
387 }
391 CylindricalFunctor divChiY)
392 {
393 p_[0] = chi_xx;
394 p_[1] = chi_xy;
395 p_[2] = chi_yy;
396 p_[3] = divChiX;
397 p_[4] = divChiY;
398 }
400 const CylindricalFunctor& xx()const{return p_[0];}
402 const CylindricalFunctor& xy()const{return p_[1];}
404 const CylindricalFunctor& yy()const{return p_[2];}
406 const CylindricalFunctor& divX()const{return p_[3];}
408 const CylindricalFunctor& divY()const{return p_[4];}
409 private:
410 std::array<CylindricalFunctor,5> p_;
411};
412
416{
425 {
426 p_[0] = v_x;
427 p_[1] = v_y;
428 p_[2] = v_z;
429 }
431 const CylindricalFunctor& x()const{return p_[0];}
433 const CylindricalFunctor& y()const{return p_[1];}
435 const CylindricalFunctor& z()const{return p_[2];}
436 private:
437 std::array<CylindricalFunctor,3> p_;
438};
439
444{
452 CylindricalFunctor divvvz
453 ): f0{v_x, v_y, v_z},
454 m_div(div), m_divvvz(divvvz) {}
460 CylindricalFunctor divvvz
461 )
462 {
463 f0.reset( v_x,v_y,v_z);
464 m_div = div;
465 m_divvvz = divvvz;
466 }
468 operator CylindricalVectorLvl0 ()const {return f0;}
470 const CylindricalFunctor& x()const{return f0.x();}
472 const CylindricalFunctor& y()const{return f0.y();}
474 const CylindricalFunctor& z()const{return f0.z();}
476 const CylindricalFunctor& div()const{return m_div;}
478 const CylindricalFunctor& divvvz()const{return m_divvvz;}
479 private:
481 CylindricalFunctor m_div, m_divvvz;
482};
483
487struct ScalarProduct : public aCylindricalFunctor<ScalarProduct>
488{
490 double do_compute( double R, double Z) const
491 {
492 return m_v.x()(R,Z)*m_w.x()(R,Z)
493 + m_v.y()(R,Z)*m_w.y()(R,Z)
494 + m_v.z()(R,Z)*m_w.z()(R,Z);
495 }
496 private:
497 CylindricalVectorLvl0 m_v, m_w;
498};
499
505struct SquareNorm : public aCylindricalFunctor<SquareNorm>
506{
508 double do_compute( double R, double Z) const
509 {
510 return sqrt(m_s(R,Z));
511 }
512 private:
513 ScalarProduct m_s;
514};
515
516
528template<class Geometry3d>
530 const dg::geo::CylindricalVectorLvl0& bhat, const Geometry3d& g)
531{
532 using host_vector = typename Geometry3d::host_vector;
534 std::array<host_vector,3> bt;
535 dg::pushForward( bhat.x(), bhat.y(), bhat.z(), bt[0], bt[1], bt[2], g);
536 std::vector<host_vector> chi(6, dg::evaluate( dg::zero,g));
537 dg::blas1::pointwiseDot( bt[0], bt[0], chi[0]);
538 dg::blas1::pointwiseDot( bt[0], bt[1], chi[1]);
539 dg::blas1::pointwiseDot( bt[0], bt[2], chi[2]);
540 dg::blas1::pointwiseDot( bt[1], bt[1], chi[3]);
541 dg::blas1::pointwiseDot( bt[1], bt[2], chi[4]);
542 dg::blas1::pointwiseDot( bt[2], bt[2], chi[5]);
543 t.idx(0,0) = 0, t.idx(0,1) = t.idx(1,0) = 1,
544 t.idx(0,2) = t.idx(2,0) = 2;
545 t.idx(1,1) = 3, t.idx(1,2) = t.idx(2,1) = 4;
546 t.idx(2,2) = 5;
547 t.values() = chi;
548 return t;
549}
561template<class Geometry3d>
563 const dg::geo::CylindricalVectorLvl0& bhat, const Geometry3d& g)
564{
565 using host_vector = typename Geometry3d::host_vector;
567 dg::SparseTensor<host_vector> m = g.metric();
568 dg::blas1::axpby( 1., m.value(0,0), -1., t.values()[0]);
569 dg::blas1::axpby( 1., m.value(0,1), -1., t.values()[1]);
570 dg::blas1::axpby( 1., m.value(0,2), -1., t.values()[2]);
571 dg::blas1::axpby( 1., m.value(1,1), -1., t.values()[3]);
572 dg::blas1::axpby( 1., m.value(1,2), -1., t.values()[4]);
573 dg::blas1::axpby( 1., m.value(2,2), -1., t.values()[5]);
574 return t;
575}
576
578}//namespace geo
579}//namespace dg
#define _ping_
DG_DEVICE T zero(T x, Ts ...xs)
void axpby(value_type alpha, const ContainerType1 &x, value_type1 beta, ContainerType &y)
void pointwiseDot(value_type alpha, const ContainerType1 &x1, const ContainerType2 &x2, value_type1 beta, ContainerType &y)
auto evaluate(Functor &&f, const Topology &g)
dg::SparseTensor< typename Geometry3d::host_vector > createAlignmentTensor(const dg::geo::CylindricalVectorLvl0 &bhat, const Geometry3d &g)
Definition fluxfunctions.h:529
dg::SparseTensor< typename Geometry3d::host_vector > createProjectionTensor(const dg::geo::CylindricalVectorLvl0 &bhat, const Geometry3d &g)
Definition fluxfunctions.h:562
void findXpoint(const CylindricalFunctorsLvl2 &psi, double &RC, double &ZC)
This function finds X-points of psi.
Definition fluxfunctions.h:354
int findOpoint(const CylindricalFunctorsLvl2 &psi, double &RC, double &ZC)
This function finds O-points of psi.
Definition fluxfunctions.h:336
int findCriticalPoint(const CylindricalFunctorsLvl2 &psi, double &RC, double &ZC)
This function finds critical points of psi (any point with vanishing gradient, including the X-point ...
Definition fluxfunctions.h:279
void pushForward(const Functor1 &vR, const Functor2 &vZ, const Functor3 &vPhi, container &vx, container &vy, container &vz, const Geometry &g)
int idx(unsigned i, unsigned j) const
std::vector< container > & values()
const container & value(size_t i, size_t j) const
real_type x0() const
dg::bc bcy() const
real_type y0() const
real_type y1() const
real_type x1() const
dg::bc bcx() const
Definition fluxfunctions.h:114
double do_compute(double R, double Z) const
Definition fluxfunctions.h:116
Constant(double c)
Definition fluxfunctions.h:115
This struct bundles a function and its first derivatives.
Definition fluxfunctions.h:186
void reset(CylindricalFunctor f, CylindricalFunctor fx, CylindricalFunctor fy)
copy given functors
Definition fluxfunctions.h:200
CylindricalFunctorsLvl1(CylindricalFunctor f, CylindricalFunctor fx, CylindricalFunctor fy)
Construct with given functors.
Definition fluxfunctions.h:196
const CylindricalFunctor & dfx() const
Definition fluxfunctions.h:209
const CylindricalFunctor & f() const
Definition fluxfunctions.h:207
const CylindricalFunctor & dfy() const
Definition fluxfunctions.h:211
CylindricalFunctorsLvl1()
the access functions are undefined as long as the class remains empty
Definition fluxfunctions.h:188
This struct bundles a function and its first and second derivatives.
Definition fluxfunctions.h:223
const CylindricalFunctor & dfxy() const
Definition fluxfunctions.h:255
const CylindricalFunctor & dfy() const
Definition fluxfunctions.h:251
CylindricalFunctorsLvl2(CylindricalFunctor f, CylindricalFunctor fx, CylindricalFunctor fy, CylindricalFunctor fxx, CylindricalFunctor fxy, CylindricalFunctor fyy)
Construct with given functors.
Definition fluxfunctions.h:232
CylindricalFunctorsLvl2()
the access functions are undefined as long as the class remains empty
Definition fluxfunctions.h:225
const CylindricalFunctor & dfx() const
Definition fluxfunctions.h:249
const CylindricalFunctor & dfxx() const
Definition fluxfunctions.h:253
void reset(CylindricalFunctor f, CylindricalFunctor fx, CylindricalFunctor fy, CylindricalFunctor fxx, CylindricalFunctor fxy, CylindricalFunctor fyy)
Replace with given Functors.
Definition fluxfunctions.h:238
const CylindricalFunctor & f() const
Definition fluxfunctions.h:247
const CylindricalFunctor & dfyy() const
Definition fluxfunctions.h:257
Definition fluxfunctions.h:365
void reset(CylindricalFunctor chi_xx, CylindricalFunctor chi_xy, CylindricalFunctor chi_yy, CylindricalFunctor divChiX, CylindricalFunctor divChiY)
replace with given functors
Definition fluxfunctions.h:389
const CylindricalFunctor & yy() const
yy component
Definition fluxfunctions.h:404
const CylindricalFunctor & divY() const
is the y-component of the divergence of the tensor
Definition fluxfunctions.h:408
CylindricalSymmTensorLvl1(CylindricalFunctor chi_xx, CylindricalFunctor chi_xy, CylindricalFunctor chi_yy, CylindricalFunctor divChiX, CylindricalFunctor divChiY)
Copy given functors.
Definition fluxfunctions.h:382
CylindricalSymmTensorLvl1()
Initialize with the identity tensor.
Definition fluxfunctions.h:369
const CylindricalFunctor & xy() const
xy component
Definition fluxfunctions.h:402
const CylindricalFunctor & divX() const
is the x-component of the divergence of the tensor
Definition fluxfunctions.h:406
const CylindricalFunctor & xx() const
xy component
Definition fluxfunctions.h:400
Definition fluxfunctions.h:416
const CylindricalFunctor & z() const
z-component of the vector
Definition fluxfunctions.h:435
CylindricalVectorLvl0(CylindricalFunctor v_x, CylindricalFunctor v_y, CylindricalFunctor v_z)
Copy given Functors.
Definition fluxfunctions.h:419
const CylindricalFunctor & y() const
y-component of the vector
Definition fluxfunctions.h:433
void reset(CylindricalFunctor v_x, CylindricalFunctor v_y, CylindricalFunctor v_z)
replace with given functors
Definition fluxfunctions.h:423
const CylindricalFunctor & x() const
x-component of the vector
Definition fluxfunctions.h:431
CylindricalVectorLvl0()
Definition fluxfunctions.h:417
This struct bundles a vector field and its divergence.
Definition fluxfunctions.h:444
void reset(CylindricalFunctor v_x, CylindricalFunctor v_y, CylindricalFunctor v_z, CylindricalFunctor div, CylindricalFunctor divvvz)
replace with given functors
Definition fluxfunctions.h:456
const CylindricalFunctor & y() const
y-component of the vector
Definition fluxfunctions.h:472
const CylindricalFunctor & x() const
x-component of the vector
Definition fluxfunctions.h:470
const CylindricalFunctor & div() const
Definition fluxfunctions.h:476
CylindricalVectorLvl1(CylindricalFunctor v_x, CylindricalFunctor v_y, CylindricalFunctor v_z, CylindricalFunctor div, CylindricalFunctor divvvz)
Copy given Functors.
Definition fluxfunctions.h:448
const CylindricalFunctor & divvvz() const
Definition fluxfunctions.h:478
const CylindricalFunctor & z() const
z-component of the vector
Definition fluxfunctions.h:474
CylindricalVectorLvl1()
the access functions are undefined as long as the class remains empty
Definition fluxfunctions.h:446
This function extends another function beyond the grid boundaries.
Definition fluxfunctions.h:144
Periodify(CylindricalFunctor functor, dg::Grid2d g)
Construct from grid.
Definition fluxfunctions.h:151
Periodify(CylindricalFunctor functor, double R0, double R1, double Z0, double Z1, dg::bc bcx, dg::bc bcy)
provide 2d grid boundaries by hand
Definition fluxfunctions.h:163
double do_compute(double R, double Z) const
Definition fluxfunctions.h:167
Inject both 2d and 3d operator() to a 2d functor.
Definition fluxfunctions.h:26
RealCylindricalFunctor()
Definition fluxfunctions.h:27
real_type operator()(real_type R, real_type Z, real_type phi) const
Definition fluxfunctions.h:42
real_type operator()(real_type R, real_type Z) const
Definition fluxfunctions.h:38
RealCylindricalFunctor(BinaryFunctor f)
Construct from any binary functor.
Definition fluxfunctions.h:35
Return scalar product of two vector fields .
Definition fluxfunctions.h:488
double do_compute(double R, double Z) const
Definition fluxfunctions.h:490
ScalarProduct(CylindricalVectorLvl0 v, CylindricalVectorLvl0 w)
Definition fluxfunctions.h:489
Return norm of scalar product of two vector fields .
Definition fluxfunctions.h:506
SquareNorm(CylindricalVectorLvl0 v, CylindricalVectorLvl0 w)
Definition fluxfunctions.h:507
double do_compute(double R, double Z) const
Definition fluxfunctions.h:508
Definition fluxfunctions.h:130
double do_compute(double R, double Z) const
Definition fluxfunctions.h:132
ZCutter(double ZX, int sign=+1)
Definition fluxfunctions.h:131
Represent functions written in cylindrical coordinates that are independent of the angle phi serving ...
Definition fluxfunctions.h:66
double operator()(double R, double Z) const
do_compute(R,Z)
Definition fluxfunctions.h:75
double operator()(double R, double Z, double phi) const
do_compute(R,Z)
Definition fluxfunctions.h:89