3#include "dg/algorithm.h"
24template<
class real_type>
34 template<
class BinaryFunctor>
42 real_type
operator()( real_type R, real_type Z, real_type phi)
const{
46 std::function<real_type(real_type,real_type)> m_f;
64template<
class Derived>
77 const Derived& underlying =
static_cast<const Derived&
>(*this);
78 return underlying.do_compute(R,Z);
91 const Derived& underlying =
static_cast<const Derived&
>(*this);
92 return underlying.do_compute(R,Z);
102 aCylindricalFunctor(
const aCylindricalFunctor&){}
106 aCylindricalFunctor& operator=(
const aCylindricalFunctor&){
return *
this;}
131 ZCutter(
double ZX,
int sign = +1): m_heavi( ZX, sign){}
165 m_g( R0, R1, Z0, Z1, 3, 10, 10, bcx, bcy), m_f(functor)
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);
213 std::array<CylindricalFunctor,3> p_;
235 f0(f,fx,fy), f1(fxx,fxy,fyy)
242 f0.reset( f,fx,fy), f1.reset(fxx,fxy,fyy);
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;
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);
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);
298 while( (eps < eps_old || eps > 1e-7) && eps > 1e-10 && counter < 100)
301 XN[0] = X[0] - Dinv*(psipZZ*psipR - psipRZ*psipZ);
302 XN[1] = X[1] - Dinv*(-psipRZ*psipR + psipRR*psipZ);
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);
314 if ( counter >= 100 || D0 == 0|| std::isnan( Dinv) )
316 RC = X[0], ZC = X[1];
317 if( Dinv > 0 && psipRR > 0)
319 if( Dinv > 0 && psipRR < 0)
339 if( point == 3 || point == 0 )
385 p_{{ chi_xx,chi_xy,chi_yy,divChiX,divChiY}}
410 std::array<CylindricalFunctor,5> p_;
437 std::array<CylindricalFunctor,3> p_;
453 ): f0{v_x, v_y, v_z},
454 m_div(div), m_divvvz(divvvz) {}
463 f0.reset( v_x,v_y,v_z);
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);
510 return sqrt(m_s(R,Z));
528template<
class Geometry3d>
532 using host_vector =
typename Geometry3d::host_vector;
534 std::array<host_vector,3> bt;
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;
561template<
class Geometry3d>
565 using host_vector =
typename Geometry3d::host_vector;
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
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