Discontinuous Galerkin Library
#include "dg/algorithm.h"
Loading...
Searching...
No Matches
functors.h
Go to the documentation of this file.
1#pragma once
2
3#include <cmath>
5#ifndef M_PI
6#define M_PI 3.14159265358979323846
7#endif
8#include <vector>
9#include <random>
10#include <functional>
11#include "blas1.h"
12#include "topology/grid.h"
13#include "topology/evaluation.h"
14#include "topology/functions.h"
18namespace dg
19{
20
23//Everything that is quite basic and simple
24
25
27template <class T = double>
28struct PLUS
29{
35 PLUS( T value): x_(value){}
37 T operator()( T x)const{ return x + x_;}
38 private:
39 T x_;
40};
41
43template< class T = double >
44struct EXP
45{
46 DG_DEVICE T operator() ( T x) const
47 {
48 return exp(x);
49 }
50};
51
53template < class T = double>
54struct LN
55{
56 DG_DEVICE T operator() (const T& x) const
57 {
58 return log(x);
59 }
60
61};
62
64template < class T = double>
65struct SQRT
66{
68 {
69 return sqrt(x);
70 }
71};
72
74struct Square
75{
76 template<class T>
77 DG_DEVICE T operator()( T x) const{ return x*x;}
78};
79
81template < class T = double>
82struct InvSqrt
83{
85 {
86 return 1./sqrt(x);
87 }
88};
89
91template <class T = double>
92struct INVERT
93{
94 DG_DEVICE T operator()( T x)const{ return 1./x;}
95};
96
98template <class T = double>
99struct ABS
100{
101 DG_DEVICE T operator()(T x)const{ return fabs(x);}
102};
103
112template <class T = double>
113struct Sign
114{
115 DG_DEVICE T operator()(T x)const{ return (T(0) < x) - (x < T(0));}
116};
117
119template <class T = double>
120struct AbsMax
121{
122 DG_DEVICE T operator() ( T x, T y) const
123 {
124 T absx = x>0 ? x : -x;
125 T absy = y>0 ? y : -y;
126 return absx > absy ? absx : absy;
127 }
128};
129
131template <class T = double>
132struct AbsMin
133{
134 DG_DEVICE T operator() (T x, T y) const
135 {
136 T absx = x<0 ? -x : x;
137 T absy = y<0 ? -y : y;
138 return absx < absy ? absx : absy;
139 }
140};
141
150template <class T = double>
152{
153 DG_DEVICE T operator()( T x)const{
154 if (x >= 0.0) return x;
155 return 0.0;
156 }
157};
158
165template <class T= double>
166struct MOD
167{
173 MOD( T m): m_m(m){}
174
176 T operator()( T x)const{
177 return (fmod(x,m_m) < 0 ) ? (m_m + fmod(x,m_m)) : fmod(x,m_m);
178 }
179 private:
180 T m_m;
181
182};
183
196template <class T>
198{
199#ifdef __CUDACC__
200 DG_DEVICE bool operator()(T x){ return !isfinite(x);}
201#else
202 bool operator()( T x){ return !std::isfinite(x);}
203#endif
204};
205
222template <class T>
224{
225#ifdef __CUDACC__
226 DG_DEVICE bool operator()(T x){
227 if( !isfinite(x))
228 return true;
229 if( x > 1e100 || x < -1e100)
230 return true;
231 return false;
232 }
233#else
234 bool operator()( T x){
235 if( !std::isfinite(x))
236 return true;
237 if( x > 1e100 || x < -1e100)
238 return true;
239 return false;
240 }
241#endif
242};
243
255struct MinMod
256{
258#ifdef __CUDACC__
259 template < class T>
260 DG_DEVICE T operator()( T x1, T x2) const
261 {
262 if( x1 > 0 && x2 > 0)
263 return min(x1,x2);
264 else if( x1 < 0 && x2 < 0)
265 return max(x1,x2);
266 return 0.;
267 }
268#else
269 template < class T>
270 T operator()( T x1, T x2) const
271 {
272 if( x1 > 0 && x2 > 0)
273 return std::min(x1,x2);
274 else if( x1 < 0 && x2 < 0)
275 return std::max(x1,x2);
276 return 0.;
277 }
278#endif
280 template<class T>
281 DG_DEVICE T operator() ( T x1, T x2, T x3)const
282 {
283 return this-> operator()( this-> operator()( x1, x2), x3);
284 }
285};
286
296{
297 template<class T>
298 DG_DEVICE T operator()( T x1, T x2) const
299 {
300 if( x1*x2 <= 0)
301 return 0.;
302 return 2.*x1*x2/(x1+x2);
303 }
304};
305
312struct Upwind
313{
314 template<class T>
315 DG_DEVICE T operator()( T velocity, T backward, T forward) const{
316 if( velocity >= 0)
317 return backward;
318 else
319 return forward;
320 }
321};
322
330{
331 template<class T>
332 DG_DEVICE T operator()( T velocity, T backward, T forward)const{
333 return velocity*m_up(velocity, backward, forward);
334 }
335 private:
336 Upwind m_up;
337};
338
348template<class Limiter>
350{
352 SlopeLimiter( Limiter l ) : m_l( l){}
353 template<class T>
354 DG_DEVICE T operator()( T v, T gm, T g0, T gp, T hm, T hp ) const{
355 if( v >= 0)
356 return +hm*m_l( g0, gm);
357 else
358 return -hp*m_l( gp, g0);
359 }
360 private:
361 Limiter m_l;
362};
363
373template<class Limiter>
375{
377 SlopeLimiterProduct( Limiter l ) : m_s( l){}
378 template<class T>
379 DG_DEVICE T operator()( T v, T gm, T g0, T gp, T hm, T hp ) const{
380 return v*m_s(v,gm,g0,gp,hm,hp);
381 }
382 private:
384};
386
387
393
396//
397
405struct Iris
406{
407 Iris( double psi_min, double psi_max ):
408 m_psimin(psi_min), m_psimax(psi_max) { }
410 double operator()(double psi)const
411 {
412 if( psi > m_psimax) return 0.;
413 if( psi < m_psimin) return 0.;
414 return 1.;
415 }
416 private:
417 double m_psimin, m_psimax;
418};
426struct Pupil
427{
428 Pupil( double psimax):
429 psimax_(psimax) { }
431 double operator()(double psi)const
432 {
433 if( psi > psimax_) return 0.;
434 return 1.;
435 }
436 private:
437 double psimax_;
438};
447{
448 PsiPupil(double psimax):
449 psimax_(psimax){ }
451 double operator()(double psi)const
452 {
453 if( psi > psimax_) return psimax_;
454 return psi;
455 }
456 private:
457 double psimax_;
458};
468{
469
478 Heaviside( double xb, int sign = +1):
479 m_xb(xb), m_s(sign){ }
480
482 double operator()(double x)const
483 {
484 if( (x < m_xb && m_s == 1) || (x > m_xb && m_s == -1)) return 0.;
485 return 1.;
486 }
487 private:
488 double m_xb;
489 int m_s;
490};
491
492
497{
498 Distance( double x0, double y0): m_x0(x0), m_y0(y0){}
500 double operator()(double x, double y){
501 return sqrt( (x-m_x0)*(x-m_x0) + (y-m_y0)*(y-m_y0));
502 }
503 private:
504 double m_x0, m_y0;
505};
512struct Line{
513 Line(double x0, double y0, double x1, double y1) :
514 m_x0(x0), m_y0(y0), m_x1(x1), m_y1(y1){}
515 double operator()(double x){
516 return m_y1*(x-m_x0)/(m_x1-m_x0) + m_y0*(x-m_x1)/(m_x0-m_x1);
517 }
518 private:
519 double m_x0, m_y0, m_x1, m_y1;
520};
521
527{
534 LinearX( double a, double b):a_(a), b_(b){}
536 double operator()(double x)const{ return a_*x+b_;}
538 double operator()( double x, double)const{ return this->operator()(x);}
540 double operator()( double x, double, double)const{ return this->operator()(x);}
541 private:
542 double a_,b_;
543};
549{
556 LinearY( double a, double b):a_(a), b_(b){}
558 double operator()( double, double y, double)const { return a_*y+b_;}
560 double operator()( double, double y)const{ return a_*y+b_;}
561 private:
562 double a_,b_;
563};
569{
576 LinearZ( double a, double b):a_(a), b_(b){}
578 double operator()( double, double, double z)const{ return a_*z+b_;}
579 private:
580 double a_,b_;
581};
589{
599 Gaussian( double x0, double y0, double sigma_x, double sigma_y, double amp)
600 : m_x0(x0), m_y0(y0), m_sigma_x(sigma_x), m_sigma_y(sigma_y), m_amp(amp){
601 assert( m_sigma_x != 0 && "sigma_x must not be 0 in Gaussian");
602 assert( m_sigma_y != 0 && "sigma_y must not be 0 in Gaussian");
603 }
616 double operator()(double x, double y) const
617 {
618 return m_amp*
619 exp( -((x-m_x0)*(x-m_x0)/2./m_sigma_x/m_sigma_x +
620 (y-m_y0)*(y-m_y0)/2./m_sigma_y/m_sigma_y) );
621 }
633 double operator()(double x, double y, double) const
634 {
635 return this->operator()(x,y);
636 }
637 private:
638 double m_x0, m_y0, m_sigma_x, m_sigma_y, m_amp;
639
640};
641
653struct Cauchy
654{
664 Cauchy( double x0, double y0, double sigma_x, double sigma_y, double amp): x0_(x0), y0_(y0), sigmaX_(sigma_x), sigmaY_(sigma_y), amp_(amp){
665 assert( sigma_x != 0 && "sigma_x must be !=0 in Cauchy");
666 assert( sigma_y != 0 && "sigma_y must be !=0 in Cauchy");
667 }
669 double operator()(double x, double y )const{
670 double xbar = (x-x0_)/sigmaX_;
671 double ybar = (y-y0_)/sigmaY_;
672 if( xbar*xbar + ybar*ybar < 1.)
673 return amp_*exp( 1. + 1./( xbar*xbar + ybar*ybar -1.) );
674 return 0.;
675 }
676 bool inside( double x, double y)const
677 {
678 double xbar = (x-x0_)/sigmaX_;
679 double ybar = (y-y0_)/sigmaY_;
680 if( xbar*xbar + ybar*ybar < 1.)
681 return true;
682 return false;
683 }
684
685 double dx( double x, double y )const{
686 double xbar = (x-x0_)/sigmaX_;
687 double ybar = (y-y0_)/sigmaY_;
688 double temp = sigmaX_*(xbar*xbar + ybar*ybar - 1.);
689 return -2.*(x-x0_)*this->operator()(x,y)/temp/temp;
690 }
691 double dxx( double x, double y)const{
692 double temp = sigmaY_*sigmaY_*(x-x0_)*(x-x0_) + sigmaX_*sigmaX_*((y-y0_)*(y-y0_) - sigmaY_*sigmaY_);
693 double bracket = sigmaX_*sigmaX_*((y-y0_)*(y-y0_)-sigmaY_*sigmaY_)*sigmaX_*sigmaX_*((y-y0_)*(y-y0_)-sigmaY_*sigmaY_)
694 -3.*sigmaY_*sigmaY_*sigmaY_*sigmaY_*(x-x0_)*(x-x0_)*(x-x0_)*(x-x0_)
695 -2.*sigmaY_*sigmaY_*sigmaX_*sigmaX_*(x-x0_)*(x-x0_)*(y-y0_)*(y-y0_);
696 return -2.*sigmaX_*sigmaX_*sigmaY_*sigmaY_*sigmaY_*sigmaY_*this->operator()(x,y)*bracket/temp/temp/temp/temp;
697 }
698 double dy( double x, double y)const{
699 double xbar = (x-x0_)/sigmaX_;
700 double ybar = (y-y0_)/sigmaY_;
701 double temp = sigmaY_*(xbar*xbar + ybar*ybar - 1.);
702 return -2.*(y-y0_)*this->operator()(x,y)/temp/temp;
703 }
704 double dyy( double x, double y)const{
705 double temp = sigmaX_*sigmaX_*(y-y0_)*(y-y0_) + sigmaY_*sigmaY_*((x-x0_)*(x-x0_) - sigmaX_*sigmaX_);
706 double bracket = sigmaY_*sigmaY_*((x-x0_)*(x-x0_)-sigmaX_*sigmaX_)*sigmaY_*sigmaY_*((x-x0_)*(x-x0_)-sigmaX_*sigmaX_)
707 -3.*sigmaX_*sigmaX_*sigmaX_*sigmaX_*(y-y0_)*(y-y0_)*(y-y0_)*(y-y0_)
708 -2.*sigmaX_*sigmaX_*sigmaY_*sigmaY_*(y-y0_)*(y-y0_)*(x-x0_)*(x-x0_);
709 return -2.*sigmaY_*sigmaY_*sigmaX_*sigmaX_*sigmaX_*sigmaX_*this->operator()(x,y)*bracket/temp/temp/temp/temp;
710 }
711 double dxy( double x, double y )const{
712 double xbar = (x-x0_)/sigmaX_;
713 double ybar = (y-y0_)/sigmaY_;
714 double temp = (xbar*xbar + ybar*ybar - 1.);
715 return 8.*xbar*ybar*this->operator()(x,y)/temp/temp/temp/sigmaX_/sigmaY_
716 + 4.*xbar*ybar*this->operator()(x,y)/temp/temp/temp/temp/sigmaX_/sigmaY_
717;
718 }
719 private:
720 double x0_, y0_, sigmaX_, sigmaY_, amp_;
721};
722
735{
743 CauchyX( double x0, double sigma_x, double amp): x0_(x0), sigmaX_(sigma_x), amp_(amp){
744 assert( sigma_x != 0 && "sigma_x must be !=0 in Cauchy");
745 }
747 double operator()(double x, double )const{
748 double xbar = (x-x0_)/sigmaX_;
749 if( xbar*xbar < 1.)
750 return amp_*exp( 1. + 1./( xbar*xbar -1.) );
751 return 0.;
752 }
753 bool inside( double x, double)const
754 {
755 double xbar = (x-x0_)/sigmaX_;
756 if( xbar*xbar < 1.)
757 return true;
758 return false;
759 }
760 private:
761 double x0_, sigmaX_, amp_;
762};
763
771{
783 Gaussian3d( double x0, double y0, double z0, double sigma_x, double sigma_y, double sigma_z, double amp)
784 : m_x0(x0), m_y0(y0), m_z0(z0), m_sigma_x(sigma_x), m_sigma_y(sigma_y), m_sigma_z(sigma_z), m_amp(amp){
785 assert( m_sigma_x != 0 && "sigma_x must be !=0 in Gaussian3d");
786 assert( m_sigma_y != 0 && "sigma_y must be !=0 in Gaussian3d");
787 assert( m_sigma_z != 0 && "sigma_z must be !=0 in Gaussian3d");
788 }
801 double operator()(double x, double y) const
802 {
803 return m_amp*
804 exp( -((x-m_x0)*(x-m_x0)/2./m_sigma_x/m_sigma_x +
805 (y-m_y0)*(y-m_y0)/2./m_sigma_y/m_sigma_y) );
806 }
820 double operator()(double x, double y, double z) const
821 {
822 return m_amp*
823 exp( -((x-m_x0)*(x-m_x0)/2./m_sigma_x/m_sigma_x +
824 (z-m_z0)*(z-m_z0)/2./m_sigma_z/m_sigma_z +
825 (y-m_y0)*(y-m_y0)/2./m_sigma_y/m_sigma_y) );
826 }
827 private:
828 double m_x0, m_y0, m_z0, m_sigma_x, m_sigma_y, m_sigma_z, m_amp;
829
830};
838{
846 GaussianX( double x0, double sigma_x, double amp)
847 :m_x0(x0), m_sigma_x(sigma_x), m_amp(amp){
848 assert( m_sigma_x != 0 && "sigma_x must be !=0 in GaussianX");
849 }
851 double operator()(double x) const
852 {
853 return m_amp* exp( -((x-m_x0)*(x-m_x0)/2./m_sigma_x/m_sigma_x ));
854 }
856 double operator()(double x, double) const
857 {
858 return this->operator()(x);
859 }
861 double operator()(double x, double, double) const
862 {
863 return this->operator()(x);
864 }
865 private:
866 double m_x0, m_sigma_x, m_amp;
867
868};
876{
884 GaussianY( double y0, double sigma_y, double amp)
885 : m_y0(y0), m_sigma_y(sigma_y), m_amp(amp){
886 assert( m_sigma_y != 0 && "sigma_x must be !=0 in GaussianY");
887 }
899 double operator()(double, double y) const
900 {
901 return m_amp*exp( -((y-m_y0)*(y-m_y0)/2./m_sigma_y/m_sigma_y) );
902 }
903 private:
904 double m_y0, m_sigma_y, m_amp;
905
906};
914{
922 GaussianZ( double z0, double sigma_z, double amp)
923 : m_z0(z0), m_sigma_z(sigma_z), m_amp(amp){
924 assert( m_sigma_z != 0 && "sigma_z must be !=0 in GaussianZ");
925 }
937 double operator()( double z) const
938 {
939 return m_amp*exp( -((z-m_z0)*(z-m_z0)/2./m_sigma_z/m_sigma_z) );
940 }
952 double operator()(double, double, double z) const
953 {
954 return m_amp*exp( -((z-m_z0)*(z-m_z0)/2./m_sigma_z/m_sigma_z) );
955 }
956 private:
957 double m_z0, m_sigma_z, m_amp;
958
959};
965{
972 IslandXY( double lambda, double eps):lambda_(lambda), eps_(eps){
973 assert( lambda != 0 && "Lambda parameter in IslandXY must not be zero!");
974 }
983 double operator()( double x, double y)const{ return lambda_*log(cosh(x/lambda_)+eps_*cos(y/lambda_));}
984 private:
985 double lambda_,eps_;
986};
991{
1000 SinXSinY( double amp, double bamp, double kx, double ky):amp_(amp), bamp_(bamp),kx_(kx),ky_(ky){}
1009 DG_DEVICE
1010 double operator()( double x, double y)const{ return bamp_+amp_*sin(x*kx_)*sin(y*ky_);}
1011 private:
1012 double amp_,bamp_,kx_,ky_;
1013};
1019{
1028 CosXCosY( double amp, double bamp, double kx, double ky):amp_(amp), bamp_(bamp),kx_(kx),ky_(ky){}
1037 DG_DEVICE
1038 double operator()( double x, double y)const{ return bamp_+amp_*cos(x*kx_)*cos(y*ky_);}
1039 private:
1040 double amp_,bamp_,kx_,ky_;
1041};
1047{
1056 SinXCosY( double amp, double bamp, double kx, double ky):amp_(amp), bamp_(bamp),kx_(kx),ky_(ky){}
1065 DG_DEVICE
1066 double operator()( double x, double y)const{ return bamp_+amp_*sin(x*kx_)*cos(y*ky_);}
1067 private:
1068 double amp_,bamp_,kx_,ky_;
1069};
1074struct SinX
1075{
1083 SinX( double amp, double bamp, double kx):amp_(amp), bamp_(bamp),kx_(kx){}
1084 DG_DEVICE
1085 double operator()( double x)const{ return bamp_+amp_*sin(x*kx_);}
1086 DG_DEVICE
1087 double operator()( double x, double)const{ return this->operator()(x);}
1088 DG_DEVICE
1089 double operator()( double x, double, double)const{ return this->operator()(x);}
1090 private:
1091 double amp_,bamp_,kx_;
1092};
1097struct SinY
1098{
1106 SinY( double amp, double bamp, double ky):amp_(amp), bamp_(bamp),ky_(ky){}
1107 DG_DEVICE
1108 double operator()( double, double y)const{ return bamp_+amp_*sin(y*ky_);}
1109 private:
1110 double amp_,bamp_,ky_;
1111};
1116struct CosY
1117{
1125 CosY( double amp, double bamp, double ky):amp_(amp), bamp_(bamp),ky_(ky){}
1126 DG_DEVICE
1127 double operator()( double, double y)const{ return bamp_+amp_*cos(y*ky_);}
1128 private:
1129 double amp_,bamp_,ky_;
1130};
1136{
1143 InvCoshXsq( double amp, double kx):m_amp(amp), m_kx(kx){}
1144 DG_DEVICE
1145 double operator()( double x)const{ return m_amp/cosh(x*m_kx)/cosh(x*m_kx);}
1146 DG_DEVICE
1147 double operator()( double x, double)const{ return this->operator()(x);}
1148 DG_DEVICE
1149 double operator()( double x, double, double)const{ return this->operator()(x);}
1150 private:
1151 double m_amp, m_kx;
1152};
1158{
1166 SinProfX( double amp, double bamp, double kx):m_amp(amp), m_bamp(bamp),m_kx(kx){}
1167 DG_DEVICE
1168 double operator()( double x)const{ return m_bamp+m_amp*(1.-sin(x*m_kx));}
1169 DG_DEVICE
1170 double operator()( double x, double)const{ return this->operator()(x);}
1171 DG_DEVICE
1172 double operator()( double x, double, double)const{ return this->operator()(x);}
1173 private:
1174 double m_amp, m_bamp, m_kx;
1175};
1181{
1189 ExpProfX( double amp, double bamp, double ln):m_amp(amp),m_bamp(bamp),m_ln(ln){
1190 assert( ln!=0 && "ln parameter must be != 0 in ExpProfX!");
1191 }
1192 DG_DEVICE
1193 double operator()( double x)const{ return m_bamp+m_amp*exp(-x/m_ln);}
1194 DG_DEVICE
1195 double operator()( double x, double)const{ return this->operator()(x);}
1196 DG_DEVICE
1197 double operator()( double x, double, double)const{ return this->operator()(x);}
1198 private:
1199 double m_amp, m_bamp, m_ln;
1200};
1201
1202
1215{
1216 GaussianDamping( double psimax, double alpha):
1217 m_psimax(psimax), m_alpha(alpha) {
1218 assert( alpha!= 0 && "Damping width in GaussianDamping must not be zero");
1219 }
1220 DG_DEVICE
1221 double operator()(double psi)const
1222 {
1223 if( psi > m_psimax + 4.*m_alpha) return 0.;
1224 if( psi < m_psimax) return 1.;
1225 return exp( -( psi-m_psimax)*( psi-m_psimax)/2./m_alpha/m_alpha);
1226 }
1227 private:
1228 double m_psimax, m_alpha;
1229};
1247 TanhProfX(double xb, double width, int sign =1,double bgamp = 0.,
1248 double profamp = 1.) :
1249 xb_(xb),w_(width), s_(sign),bga_(bgamp),profa_(profamp) {
1250 assert( width != 0&& "Width in TanhProfX must not be zero!");
1251 }
1252 DG_DEVICE
1253 double operator() (double x)const
1254 {
1255 return profa_*0.5*(1.+s_*tanh((x-xb_)/w_))+bga_;
1256 }
1257 DG_DEVICE
1258 double operator()( double x, double)const{ return this->operator()(x);}
1259 DG_DEVICE
1260 double operator()( double x, double, double)const{ return this->operator()(x);}
1261 private:
1262 double xb_;
1263 double w_;
1264 int s_;
1265 double bga_;
1266 double profa_;
1267};
1268
1290 PolynomialHeaviside(double xb, double a, int sign = +1) :
1291 x0(xb), a(a), m_s(sign){
1292 assert( a!=0 && "PolynomialHeaviside width must not be zero");
1293 }
1294 DG_DEVICE
1295 double operator() (double x)const
1296 {
1297 if( m_s == -1) x = 2*x0-x; //mirror
1298 if ( x < x0-a) return 0;
1299 if ( x > x0+a) return 1;
1300 return ((16.*a*a*a - 29.*a*a*(x - x0)
1301 + 20.*a*(x - x0)*(x - x0)
1302 - 5.*(x - x0)*(x-x0)*(x-x0))
1303 *(a + x - x0)*(a + x - x0)
1304 *(a + x - x0)*(a + x - x0))/(32.*a*a*a * a*a*a*a);
1305 }
1306 DG_DEVICE
1307 double operator()( double x, double)const{ return this->operator()(x);}
1308 DG_DEVICE
1309 double operator()( double x, double, double)const{ return this->operator()(x);}
1310 private:
1311 double x0, a;
1312 int m_s;
1313};
1314
1342 PolynomialRectangle(double xl, double al, double xr, double ar) :
1343 m_hl( xl, al, +1), m_hr( xr, ar, -1) {
1344 assert( xl < xr && "left boundary must be left of right boundary");
1345 }
1346 DG_DEVICE
1347 double operator() (double x)const
1348 {
1349 return m_hl(x)*m_hr(x);
1350 }
1351 DG_DEVICE
1352 double operator()( double x, double)const{ return this->operator()(x);}
1353 DG_DEVICE
1354 double operator()( double x, double, double)const{ return this->operator()(x);}
1355 private:
1356 PolynomialHeaviside m_hl, m_hr;
1357};
1358
1380 IPolynomialHeaviside(double xb, double a, int sign = +1) :
1381 x0(xb), a(a), m_s(sign){
1382 assert( a!=0 && "IPolynomialHeaviside width must not be zero");
1383 }
1384 DG_DEVICE
1385 double operator() (double x)const
1386 {
1387 if( m_s == -1) x = 2*x0-x; //mirror
1388 double result;
1389 if ( x < x0-a) result = x0;
1390 else if ( x > x0+a) result = x;
1391 else
1392 result = x0 + ((35.* a*a*a - 47.* a*a*(x - x0) + 25.*a*(x - x0)*(x-x0)
1393 - 5.*(x - x0)*(x-x0)*(x-x0))
1394 *(a+x-x0)*(a+x-x0)*(a+x-x0)*(a+x-x0)*(a+x-x0))
1395 /(256.*a*a*a * a*a*a*a);
1396 if ( m_s == +1) return result;
1397 return 2*x0 - result;
1398
1399 }
1400 DG_DEVICE
1401 double operator()( double x, double)const{ return this->operator()(x);}
1402 DG_DEVICE
1403 double operator()( double x, double, double)const{ return this->operator()(x);}
1404 private:
1405 double x0, a;
1406 int m_s;
1407};
1408
1432 DPolynomialHeaviside(double xb, double a, int = +1) :
1433 x0(xb), a(a){
1434 assert( a!=0 && "DPolynomialHeaviside width must not be zero");
1435 }
1436 DG_DEVICE
1437 double operator() (double x)const
1438 {
1439 if ( (x < x0-a) || (x > x0+a)) return 0;
1440 return (35.*(a+x-x0)*(a+x-x0)*(a+x-x0)*(a-x+x0)*(a-x+x0)*(a-x+x0))
1441 /(32.*a*a*a * a*a*a*a);
1442 }
1443 DG_DEVICE
1444 double operator()( double x, double)const{ return this->operator()(x);}
1445 DG_DEVICE
1446 double operator()( double x, double, double)const{ return this->operator()(x);}
1447 private:
1448 double x0, a;
1449};
1450
1465{
1479 ExponentialFilter( double alpha, double eta_c, unsigned order, unsigned n):
1480 m_alpha(alpha), m_etac(eta_c), m_s(order), m_n(n) {}
1481 double operator()( unsigned i) const
1482 {
1483 double eta = (double)i/(double)(m_n-1);
1484 if( m_n == 1) eta = 0.;
1485 if( eta < m_etac)
1486 return 1.;
1487 if( eta <= 1.)
1488 return exp( -m_alpha*pow( (eta-m_etac)/(1.-m_etac), 2*m_s));
1489 return 0;
1490 }
1491 private:
1492 double m_alpha, m_etac;
1493 unsigned m_s, m_n;
1494};
1495
1510struct Lamb
1511{
1520 Lamb( double x0, double y0, double R, double U):R_(R), U_(U), x0_(x0), y0_(y0)
1521 {
1522 gamma_ = 3.83170597020751231561;
1523 lambda_ = gamma_/R;
1524#ifdef _MSC_VER
1525 j_ = _j0(gamma_);
1526#else
1527 j_ = j0( gamma_);
1528#endif
1529 //std::cout << r_ <<u_<<x0_<<y0_<<lambda_<<gamma_<<j_<<std::endl;
1530 }
1539 DG_DEVICE
1540 double operator() (double x, double y)const
1541 {
1542 double radius = sqrt( (x-x0_)*(x-x0_) + (y-y0_)*(y-y0_));
1543 double theta = atan2( (y-y0_),(x-x0_));
1544
1545 if( radius <= R_)
1546#ifdef _MSC_VER
1547 return 2.*lambda_*U_*_j1(lambda_*radius)/j_*cos( theta);
1548#else
1549 return 2.*lambda_*U_*j1( lambda_*radius)/j_*cos( theta);
1550#endif
1551 return 0;
1552 }
1560 double enstrophy( ) const { return M_PI*U_*U_*gamma_*gamma_;}
1561
1568 double energy() const { return 2.*M_PI*R_*R_*U_*U_;}
1569 private:
1570 double R_, U_, x0_, y0_, lambda_, gamma_, j_;
1571};
1572
1595{
1606 Vortex( double x0, double y0, unsigned state,
1607 double R, double u_dipole, double kz = 0):
1608 x0_(x0), y0_(y0), s_(state), R_(R), u_d( u_dipole), kz_(kz){
1609 g_[0] = 3.831896621;
1610 g_[1] = -3.832353624;
1611 g_[2] = 7.016;
1612 b_[0] = 0.03827327723;
1613 b_[1] = 0.07071067810 ;
1614 b_[2] = 0.07071067810 ;
1615 }
1632 DG_DEVICE
1633 double operator()( double x, double y)const
1634 {
1635 double r = sqrt( (x-x0_)*(x-x0_)+(y-y0_)*(y-y0_));
1636 double theta = atan2( y-y0_, x-x0_);
1637 double beta = b_[s_];
1638 double norm = 1.2965125;
1639
1640 if( r/R_<=1.)
1641 return u_d*(
1642 r *( 1 +beta*beta/g_[s_]/g_[s_] )
1643#ifdef _MSC_VER
1644 - R_* beta*beta/g_[s_]/g_[s_] *_j1(g_[s_]*r/R_)/_j1(g_[s_])
1645#else
1646 - R_ * beta*beta/g_[s_]/g_[s_] * j1(g_[s_]*r/R_)/ j1(g_[s_])
1647#endif
1648 )*cos(theta)/norm;
1649 return u_d * R_* bessk1(beta*r/R_)/bessk1(beta)*cos(theta)/norm;
1650 // TODO Can these be replaced by std::cyl_bessel_k? Not sure what to do on device though
1651 }
1669 DG_DEVICE
1670 double operator()( double x, double y, double z)const
1671 {
1672 return this->operator()(x,y)*cos(kz_*z);
1673 }
1674 private:
1675 // Returns the modified Bessel function K1(x) for positive real x.
1676 DG_DEVICE
1677 double bessk1(double x)const
1678 {
1679 double y,ans;
1680 if (x <= 2.0)
1681 {
1682 y=x*x/4.0;
1683 ans = (log(x/2.0)*bessi1(x))+(1.0/x)*(1.0+y*(0.15443144 +
1684 y*(-0.67278579+y*(-0.18156897+y*(-0.1919402e-1 +
1685 y*(-0.110404e-2+y*(-0.4686e-4)))))));
1686 }
1687 else
1688 {
1689 y=2.0/x;
1690 ans = (exp(-x)/sqrt(x))*(1.25331414+y*(0.23498619 +
1691 y*(-0.3655620e-1+y*(0.1504268e-1+y*(-0.780353e-2 +
1692 y*(0.325614e-2+y*(-0.68245e-3)))))));
1693 }
1694 return ans;
1695 }
1696 //Returns the modified Bessel function I1(x) for any real x.
1697 DG_DEVICE
1698 double bessi1(double x) const
1699 {
1700 double ax,ans;
1701 double y;
1702 if ((ax=fabs(x)) < 3.75)
1703 {
1704 y=x/3.75;
1705 y*=y;
1706 ans = ax*(0.5+y*(0.87890594+y*(0.51498869+y*(0.15084934 +
1707 y*(0.2658733e-1+y*(0.301532e-2+y*0.32411e-3))))));
1708 }
1709 else
1710 {
1711 y=3.75/ax;
1712 ans = 0.2282967e-1+y*(-0.2895312e-1+y*(0.1787654e-1 -
1713 y*0.420059e-2)); ans=0.39894228+y*(-0.3988024e-1+
1714 y*(-0.362018e-2 +y*(0.163801e-2+y*(-0.1031555e-1+y*ans))));
1715 ans *= (exp(ax)/sqrt(ax));
1716 }
1717 return x < 0.0 ? -ans : ans;
1718 }
1719 double x0_, y0_;
1720 unsigned s_;
1721 double R_, b_[3], u_d;
1722 double g_[3];
1723 double kz_;
1724};
1725
1744struct BathRZ{
1756 BathRZ( unsigned N_kR, unsigned N_kZ, double R_min, double Z_min, double gamma, double L_E, double amp) :
1757 N_kR_(N_kR), N_kZ_(N_kZ),
1758 R_min_(R_min), Z_min_(Z_min),
1759 gamma_(gamma), L_E_(L_E) , amp_(amp),
1760 kvec( N_kR_*N_kZ_, 0), sqEkvec(kvec), unif1(kvec), unif2(kvec),
1761 normal1(kvec), normal2(kvec), alpha(kvec), theta(kvec)
1762 {
1763 double N_kR2=(double)(N_kR_*N_kR_);
1764 double N_kZ2=(double)(N_kZ_*N_kZ_);
1765 double N_k= sqrt(N_kR2+N_kZ2);
1766
1767 norm_=sqrt(2./(double)N_kR_/(double)N_kZ_);
1768 double tpi=2.*M_PI, tpi2=tpi*tpi;
1769 double k0= tpi*L_E_/N_k;
1770 double N_kRh = N_kR_/2.;
1771 double N_kZh = N_kZ_/2.;
1772
1773 std::minstd_rand generator;
1774 std::normal_distribution<double> ndistribution( 0.0, 1.0); // ( mean, stddev)
1775 std::uniform_real_distribution<double> udistribution(0.0,tpi); //between [0 and 2pi)
1776 for (unsigned j=1;j<=N_kZ_;j++)
1777 {
1778 double kZ2=tpi2*(j-N_kZh)*(j-N_kZh)/(N_kZ2);
1779 for (unsigned i=1;i<=N_kR_;i++)
1780 {
1781 double kR2=tpi2*(i-N_kRh)*(i-N_kRh)/(N_kR2);
1782 int z=(j-1)*(N_kR_)+(i-1);
1783 kvec[z]= sqrt(kR2 + kZ2); //radial k number
1784 sqEkvec[z]=pow(kvec[z]*4.*k0/(kvec[z]+k0)/(kvec[z]+k0),gamma_/2.); //Energie in k space with max at 1.
1785 unif1[z]=cos(udistribution(generator));
1786 unif2[z]=sin(udistribution(generator));
1787 normal1[z]=ndistribution(generator);
1788 normal2[z]=ndistribution(generator);
1789 alpha[z]=sqrt(normal1[z]*normal1[z]+normal2[z]*normal2[z]);
1790 theta[z]=atan2(normal2[z],normal1[z]);
1791 }
1792 }
1793
1794 }
1818 double operator()(double R, double Z)const
1819 {
1820 double f, kappa, RR, ZZ;
1821 RR=R-R_min_;
1822 ZZ=Z-Z_min_;
1823 f=0.;
1824 for (unsigned j=0;j<N_kZ_;j++)
1825 {
1826 for (unsigned i=0;i<N_kR_;i++)
1827 {
1828 int z=j*N_kR_+i;
1829 kappa= RR*unif1[z]+ZZ*unif2[z];
1830 f+= sqEkvec[z]*alpha[z]*cos(kvec[z]*kappa+theta[z]);
1831 }
1832 }
1833 return amp_*norm_*f;
1834 }
1859 double operator()(double R, double Z, double)const {
1860 double f, kappa;
1861 double RR, ZZ;
1862 RR=R-R_min_;
1863 ZZ=Z-Z_min_;
1864 f=0;
1865 for (unsigned j=0;j<N_kZ_;j++)
1866 {
1867 for (unsigned i=0;i<N_kR_;i++)
1868 {
1869 int z=(j)*(N_kR_)+(i);
1870 kappa= RR*unif1[z]+ZZ*unif2[z];
1871 f+= sqEkvec[z]*alpha[z]*cos(kvec[z]*kappa+theta[z]);
1872 }
1873 }
1874 return amp_*norm_*f;
1875 }
1876 private:
1877 unsigned N_kR_,N_kZ_;
1878 double R_min_, Z_min_;
1879 double gamma_, L_E_;
1880 double amp_;
1881 double norm_;
1882 std::vector<double> kvec;
1883 std::vector<double> sqEkvec;
1884 std::vector<double> unif1, unif2, normal1,normal2,alpha,theta;
1885};
1886
1893{
1895 Horner2d(): m_c( 1, 1), m_M(1), m_N(1){}
1896
1904 Horner2d( std::vector<double> c, unsigned M, unsigned N): m_c(c), m_M(M), m_N(N){}
1905 double operator()( double x, double y) const
1906 {
1907 std::vector<double> cx( m_M);
1908 for( unsigned i=0; i<m_M; i++)
1909 cx[i] = horner( &m_c[i*m_N], m_N, y);
1910 return horner( &cx[0], m_M, x);
1911 }
1912 private:
1913 double horner( const double * c, unsigned M, double x) const
1914 {
1915 double b = c[M-1];
1916 for( unsigned i=0; i<M-1; i++)
1917 b = c[M-2-i] + b*x;
1918 return b;
1919 }
1920 std::vector<double> m_c;
1921 unsigned m_M, m_N;
1922};
1923
1924
1929template <class container = thrust::host_vector<double> >
1931{
1937 Histogram(const dg::Grid1d& g1d, const std::vector<double>& in) :
1938 g1d_(g1d),
1939 in_(in),
1940 binwidth_(g1d_.h()),
1941 count_(dg::evaluate(dg::zero,g1d_))
1942 {
1943 for (unsigned j=0;j<in_.size();j++)
1944 {
1945 unsigned bin =floor( (in_[j]-g1d_.x0())/binwidth_ );
1946 bin = std::max(bin,(unsigned) 0);
1947 bin = std::min(bin,(unsigned)(g1d_.size()-1));
1948 count_[bin ]+=1.;
1949 }
1950 //Normalize
1951 unsigned Ampmax = (unsigned)thrust::reduce( count_.begin(), count_.end(),0., thrust::maximum<double>() );
1952 dg::blas1::scal(count_,1./Ampmax);
1953
1954 }
1955
1961 double binwidth() {return binwidth_;}
1969 double operator()(double x)const
1970 {
1971 unsigned bin = floor((x-g1d_.x0())/binwidth_+0.5);
1972 bin = std::max(bin,(unsigned) 0);
1973 bin = std::min(bin,(unsigned)(g1d_.size()-1));
1974 return count_[bin];
1975 }
1976
1977 private:
1978 dg::Grid1d g1d_;
1979 const std::vector<double> in_;
1980 double binwidth_;
1981 container count_;
1982};
1983
1988template <class container = thrust::host_vector<double> >
1990{
1997 Histogram2D(const dg::Grid2d& g2d, const std::vector<double>& inx,const std::vector<double>& iny) :
1998 g2d_(g2d),
1999 inx_(inx),
2000 iny_(iny),
2001 binwidthx_(g2d_.hx()),
2002 binwidthy_(g2d_.hy()),
2003 count_(dg::evaluate(dg::zero,g2d_))
2004 {
2005
2006 for (unsigned j=0;j<iny_.size();j++)
2007 {
2008 unsigned biny =floor((iny_[j]-g2d_.y0())/binwidthy_) ;
2009 biny = std::max(biny,(unsigned) 0);
2010 biny = std::min(biny,(unsigned)(g2d_.Ny()-1));
2011
2012 unsigned binx =floor((inx_[j]-g2d_.x0())/binwidthx_) ;
2013 binx = std::max(binx,(unsigned) 0);
2014 binx = std::min(binx,(unsigned)(g2d_.Nx()-1));
2015 count_[biny*g2d_.Nx()+binx ]+=1.;
2016
2017 }
2018 //Normalize
2019 unsigned Ampmax = (unsigned)thrust::reduce( count_.begin(), count_.end(),0.,thrust::maximum<double>() );
2020 dg::blas1::scal(count_, 1./Ampmax);
2021
2022 }
2023
2032 double operator()(double x, double y)const
2033 {
2034 unsigned binx = floor((x-g2d_.x0())/binwidthx_+0.5) ;
2035 binx = std::max(binx,(unsigned) 0);
2036 binx = std::min(binx,(unsigned)(g2d_.Nx()-1));
2037 unsigned biny = floor((y-g2d_.y0())/binwidthy_+0.5) ;
2038 biny = std::max(biny,(unsigned) 0);
2039 biny = std::min(biny,(unsigned)(g2d_.Ny()-1));
2040 return count_[biny*g2d_.Nx()+binx ];
2041
2042 }
2043 private:
2044 dg::Grid2d g2d_;
2045 const std::vector<double> inx_,iny_;
2046 double binwidthx_,binwidthy_;
2047 container count_;
2048};
2049
2050
2051
2059{
2066 WallDistance( std::vector<double> vertical, std::vector<double> horizontal) :
2067 m_vertical(vertical), m_horizontal( horizontal) {}
2073 WallDistance( dg::Grid2d walls) : m_vertical({walls.x0(), walls.x1()}),
2074 m_horizontal({walls.y0(), walls.y1()}){}
2078 double operator() (double R, double Z) const
2079 {
2080 std::vector<double> dist( 1, 1e100); //fill in at least one (large) number in case vectors are empty)
2081 for( auto v : m_vertical)
2082 dist.push_back(fabs( R-v));
2083 for( auto h : m_horizontal)
2084 dist.push_back(fabs( Z-h));
2085 return *std::min_element( dist.begin(), dist.end());
2086 }
2087 private:
2088 std::vector<double> m_vertical;
2089 std::vector<double> m_horizontal;
2090};
2091
2092
2094} //namespace dg
2095
Function discretization routines.
Some utility functions for the dg::evaluate routines.
#define M_PI
M_PI is non-standard ... so MSVC complains.
Definition functors.h:6
base topology classes
DG_DEVICE T zero(T, Ts ...)
This enum can be used in dg::evaluate.
Definition functions.h:19
void scal(ContainerType &x, value_type alpha)
Definition blas1.h:263
@ z
z direction
@ backward
backward derivative (cell to the left and current cell)
Definition enums.h:99
@ forward
forward derivative (cell to the right and current cell)
Definition enums.h:98
@ y
y direction
@ x
x direction
auto evaluate(Functor &&f, const Topology &g)
Evaluate a function on grid coordinates
Definition evaluation.h:74
#define DG_DEVICE
Expands to __host__ __device__ if compiled with nvcc else is empty.
Definition dg_doc.h:378
const double m
const double alpha
const double n
This is the namespace for all functions and classes defined and used by the discontinuous Galerkin li...
double amp
const double beta
Definition functors.h:100
DG_DEVICE T operator()(T x) const
Definition functors.h:101
Definition functors.h:121
DG_DEVICE T operator()(T x, T y) const
Definition functors.h:122
Definition functors.h:133
DG_DEVICE T operator()(T x, T y) const
Definition functors.h:134
Definition functors.h:1744
double operator()(double R, double Z) const
Return the value of the Bath.
Definition functors.h:1818
BathRZ(unsigned N_kR, unsigned N_kZ, double R_min, double Z_min, double gamma, double L_E, double amp)
Functor returning a random field in the RZ-plane or in the first RZ-plane.
Definition functors.h:1756
double operator()(double R, double Z, double) const
Return the value of the Bath.
Definition functors.h:1859
Definition functors.h:654
double dxx(double x, double y) const
Definition functors.h:691
Cauchy(double x0, double y0, double sigma_x, double sigma_y, double amp)
A blob that drops to zero.
Definition functors.h:664
double dxy(double x, double y) const
Definition functors.h:711
double dx(double x, double y) const
Definition functors.h:685
DG_DEVICE double operator()(double x, double y) const
Definition functors.h:669
double dy(double x, double y) const
Definition functors.h:698
double dyy(double x, double y) const
Definition functors.h:704
bool inside(double x, double y) const
Definition functors.h:676
Definition functors.h:735
bool inside(double x, double) const
Definition functors.h:753
CauchyX(double x0, double sigma_x, double amp)
A 1D-blob that drops to zero.
Definition functors.h:743
DG_DEVICE double operator()(double x, double) const
Definition functors.h:747
Definition functors.h:1019
CosXCosY(double amp, double bamp, double kx, double ky)
Construct.
Definition functors.h:1028
DG_DEVICE double operator()(double x, double y) const
Return profile.
Definition functors.h:1038
Definition functors.h:1117
DG_DEVICE double operator()(double, double y) const
Definition functors.h:1127
CosY(double amp, double bamp, double ky)
Construct.
Definition functors.h:1125
The derivative of PolynomialHeaviside approximates delta(x)
Definition functors.h:1421
DG_DEVICE double operator()(double x) const
Definition functors.h:1437
DG_DEVICE double operator()(double x, double, double) const
Definition functors.h:1446
DPolynomialHeaviside(double xb, double a, int=+1)
Construct with xb, width and sign.
Definition functors.h:1432
DG_DEVICE double operator()(double x, double) const
Definition functors.h:1444
Definition functors.h:497
Distance(double x0, double y0)
Definition functors.h:498
DG_DEVICE double operator()(double x, double y)
Definition functors.h:500
Definition functors.h:45
DG_DEVICE T operator()(T x) const
Definition functors.h:46
Definition functors.h:1181
DG_DEVICE double operator()(double x) const
Definition functors.h:1193
DG_DEVICE double operator()(double x, double, double) const
Definition functors.h:1197
ExpProfX(double amp, double bamp, double ln)
Construct with three coefficients.
Definition functors.h:1189
DG_DEVICE double operator()(double x, double) const
Definition functors.h:1195
Definition functors.h:1465
ExponentialFilter(double alpha, double eta_c, unsigned order, unsigned n)
Create exponential filter .
Definition functors.h:1479
double operator()(unsigned i) const
Definition functors.h:1481
Definition functors.h:771
DG_DEVICE double operator()(double x, double y) const
Return a 2d Gaussian.
Definition functors.h:801
DG_DEVICE double operator()(double x, double y, double z) const
Return the value of the Gaussian.
Definition functors.h:820
Gaussian3d(double x0, double y0, double z0, double sigma_x, double sigma_y, double sigma_z, double amp)
Functor returning a Gaussian.
Definition functors.h:783
Definition functors.h:1215
DG_DEVICE double operator()(double psi) const
Definition functors.h:1221
GaussianDamping(double psimax, double alpha)
Definition functors.h:1216
Definition functors.h:589
DG_DEVICE double operator()(double x, double y) const
Return the value of the Gaussian.
Definition functors.h:616
Gaussian(double x0, double y0, double sigma_x, double sigma_y, double amp)
Functor returning a Gaussian.
Definition functors.h:599
DG_DEVICE double operator()(double x, double y, double) const
Return the value of the Gaussian.
Definition functors.h:633
Definition functors.h:838
DG_DEVICE double operator()(double x, double, double) const
Definition functors.h:861
DG_DEVICE double operator()(double x) const
Definition functors.h:851
DG_DEVICE double operator()(double x, double) const
Definition functors.h:856
GaussianX(double x0, double sigma_x, double amp)
A Gaussian in x.
Definition functors.h:846
Definition functors.h:876
GaussianY(double y0, double sigma_y, double amp)
Functor returning a gaussian.
Definition functors.h:884
DG_DEVICE double operator()(double, double y) const
Return the value of the gaussian.
Definition functors.h:899
Definition functors.h:914
DG_DEVICE double operator()(double z) const
Return the value of the gaussian.
Definition functors.h:937
GaussianZ(double z0, double sigma_z, double amp)
Functor returning a gaussian.
Definition functors.h:922
DG_DEVICE double operator()(double, double, double z) const
Return the value of the gaussian.
Definition functors.h:952
Definition functors.h:468
DG_DEVICE double operator()(double x) const
Definition functors.h:482
Heaviside(double xb, int sign=+1)
Construct with xb and sign.
Definition functors.h:478
Compute a histogram on a 2D grid.
Definition functors.h:1990
Histogram2D(const dg::Grid2d &g2d, const std::vector< double > &inx, const std::vector< double > &iny)
Construct a histogram from number of bins and an input vector.
Definition functors.h:1997
double operator()(double x, double y) const
Access computed histogram.
Definition functors.h:2032
Compute a histogram on a 1D grid.
Definition functors.h:1931
Histogram(const dg::Grid1d &g1d, const std::vector< double > &in)
Construct a histogram from number of bins and an input vector.
Definition functors.h:1937
double binwidth()
get binwidth
Definition functors.h:1961
double operator()(double x) const
Access computed histogram.
Definition functors.h:1969
Definition functors.h:1893
Horner2d()
Initialize 1 coefficient to 1.
Definition functors.h:1895
double operator()(double x, double y) const
Definition functors.h:1905
Horner2d(std::vector< double > c, unsigned M, unsigned N)
Initialize coefficients and dimensions.
Definition functors.h:1904
Definition functors.h:93
DG_DEVICE T operator()(T x) const
Definition functors.h:94
The integral of PolynomialHeaviside approximates xH(x)
Definition functors.h:1372
DG_DEVICE double operator()(double x, double) const
Definition functors.h:1401
IPolynomialHeaviside(double xb, double a, int sign=+1)
Construct with xb, width and sign.
Definition functors.h:1380
DG_DEVICE double operator()(double x) const
Definition functors.h:1385
DG_DEVICE double operator()(double x, double, double) const
Definition functors.h:1403
Definition functors.h:198
bool operator()(T x)
Definition functors.h:202
Definition functors.h:224
bool operator()(T x)
Definition functors.h:234
Definition functors.h:1136
DG_DEVICE double operator()(double x, double, double) const
Definition functors.h:1149
DG_DEVICE double operator()(double x, double) const
Definition functors.h:1147
DG_DEVICE double operator()(double x) const
Definition functors.h:1145
InvCoshXsq(double amp, double kx)
Construct with two coefficients.
Definition functors.h:1143
Definition functors.h:83
DG_DEVICE T operator()(T x) const
Definition functors.h:84
Definition functors.h:406
DG_DEVICE double operator()(double psi) const
Definition functors.h:410
Iris(double psi_min, double psi_max)
Definition functors.h:407
Definition functors.h:965
DG_DEVICE double operator()(double x, double y) const
Return profile.
Definition functors.h:983
IslandXY(double lambda, double eps)
Construct Island.
Definition functors.h:972
Definition functors.h:55
DG_DEVICE T operator()(const T &x) const
Definition functors.h:56
Definition functors.h:1511
double enstrophy() const
The total enstrophy of the dipole.
Definition functors.h:1560
Lamb(double x0, double y0, double R, double U)
Functor returning a Lamb-dipole.
Definition functors.h:1520
DG_DEVICE double operator()(double x, double y) const
Return the value of the dipole.
Definition functors.h:1540
double energy() const
The total energy of the dipole.
Definition functors.h:1568
Definition functors.h:512
double operator()(double x)
Definition functors.h:515
Line(double x0, double y0, double x1, double y1)
Definition functors.h:513
Definition functors.h:527
DG_DEVICE double operator()(double x) const
Definition functors.h:536
LinearX(double a, double b)
Construct with two coefficients.
Definition functors.h:534
DG_DEVICE double operator()(double x, double, double) const
Definition functors.h:540
DG_DEVICE double operator()(double x, double) const
Definition functors.h:538
Definition functors.h:549
LinearY(double a, double b)
Construct with two coefficients.
Definition functors.h:556
DG_DEVICE double operator()(double, double y) const
Definition functors.h:560
DG_DEVICE double operator()(double, double y, double) const
Definition functors.h:558
Definition functors.h:569
LinearZ(double a, double b)
Construct with two coefficients.
Definition functors.h:576
DG_DEVICE double operator()(double, double, double z) const
Definition functors.h:578
x mod m > 0 ? x mod m : x mod m + m
Definition functors.h:167
MOD(T m)
Construct from modulo.
Definition functors.h:173
DG_DEVICE T operator()(T x) const
Definition functors.h:176
Definition functors.h:256
T operator()(T x1, T x2) const
Definition functors.h:270
Definition functors.h:29
PLUS(T value)
Construct.
Definition functors.h:35
DG_DEVICE T operator()(T x) const
Definition functors.h:37
Definition functors.h:152
DG_DEVICE T operator()(T x) const
Definition functors.h:153
Definition functors.h:1281
DG_DEVICE double operator()(double x, double, double) const
Definition functors.h:1309
DG_DEVICE double operator()(double x, double) const
Definition functors.h:1307
DG_DEVICE double operator()(double x) const
Definition functors.h:1295
PolynomialHeaviside(double xb, double a, int sign=+1)
Construct with xb, width and sign.
Definition functors.h:1290
Definition functors.h:1333
DG_DEVICE double operator()(double x, double, double) const
Definition functors.h:1354
DG_DEVICE double operator()(double x, double) const
Definition functors.h:1352
DG_DEVICE double operator()(double x) const
Definition functors.h:1347
PolynomialRectangle(double xl, double al, double xr, double ar)
Construct with xb, width and sign.
Definition functors.h:1342
Definition functors.h:447
PsiPupil(double psimax)
Definition functors.h:448
DG_DEVICE double operator()(double psi) const
Definition functors.h:451
Definition functors.h:427
Pupil(double psimax)
Definition functors.h:428
DG_DEVICE double operator()(double psi) const
Definition functors.h:431
Definition functors.h:66
DG_DEVICE T operator()(T x) const
Definition functors.h:67
Definition functors.h:114
DG_DEVICE T operator()(T x) const
Definition functors.h:115
Definition functors.h:1158
SinProfX(double amp, double bamp, double kx)
Construct.
Definition functors.h:1166
DG_DEVICE double operator()(double x) const
Definition functors.h:1168
DG_DEVICE double operator()(double x, double, double) const
Definition functors.h:1172
DG_DEVICE double operator()(double x, double) const
Definition functors.h:1170
Definition functors.h:1047
SinXCosY(double amp, double bamp, double kx, double ky)
Construct.
Definition functors.h:1056
DG_DEVICE double operator()(double x, double y) const
Return profile.
Definition functors.h:1066
Definition functors.h:1075
DG_DEVICE double operator()(double x) const
Definition functors.h:1085
DG_DEVICE double operator()(double x, double, double) const
Definition functors.h:1089
DG_DEVICE double operator()(double x, double) const
Definition functors.h:1087
SinX(double amp, double bamp, double kx)
Construct.
Definition functors.h:1083
Definition functors.h:991
SinXSinY(double amp, double bamp, double kx, double ky)
Construct.
Definition functors.h:1000
DG_DEVICE double operator()(double x, double y) const
Return profile.
Definition functors.h:1010
Definition functors.h:1098
SinY(double amp, double bamp, double ky)
Construct.
Definition functors.h:1106
DG_DEVICE double operator()(double, double y) const
Definition functors.h:1108
Definition functors.h:350
DG_DEVICE T operator()(T v, T gm, T g0, T gp, T hm, T hp) const
Definition functors.h:354
SlopeLimiter()
Definition functors.h:351
SlopeLimiter(Limiter l)
Definition functors.h:352
Definition functors.h:375
DG_DEVICE T operator()(T v, T gm, T g0, T gp, T hm, T hp) const
Definition functors.h:379
SlopeLimiterProduct(Limiter l)
Definition functors.h:377
SlopeLimiterProduct()
Definition functors.h:376
Definition functors.h:75
DG_DEVICE T operator()(T x) const
Definition functors.h:77
Definition functors.h:1236
DG_DEVICE double operator()(double x, double) const
Definition functors.h:1258
DG_DEVICE double operator()(double x, double, double) const
Definition functors.h:1260
DG_DEVICE double operator()(double x) const
Definition functors.h:1253
TanhProfX(double xb, double width, int sign=1, double bgamp=0., double profamp=1.)
Construct with xb, width and sign.
Definition functors.h:1247
Definition functors.h:313
DG_DEVICE T operator()(T velocity, T backward, T forward) const
Definition functors.h:315
Definition functors.h:330
DG_DEVICE T operator()(T velocity, T backward, T forward) const
Definition functors.h:332
Definition functors.h:296
DG_DEVICE T operator()(T x1, T x2) const
Definition functors.h:298
Definition functors.h:1595
DG_DEVICE double operator()(double x, double y) const
Definition functors.h:1633
Vortex(double x0, double y0, unsigned state, double R, double u_dipole, double kz=0)
Definition functors.h:1606
DG_DEVICE double operator()(double x, double y, double z) const
Definition functors.h:1670
Shortest Distance to a collection of vertical and horizontal lines.
Definition functors.h:2059
WallDistance(dg::Grid2d walls)
Allocate lines.
Definition functors.h:2073
WallDistance(std::vector< double > vertical, std::vector< double > horizontal)
Allocate lines.
Definition functors.h:2066
double operator()(double R, double Z) const
Distance to closest wall in a box.
Definition functors.h:2078
real_type x0() const
Equivalent to p(0)
Definition grid.h:285
unsigned size() const
The total number of points.
Definition grid.h:532
unsigned Nx() const
Equivalent to N(0)
Definition grid.h:334
real_type y0() const
Equivalent to p(2)
Definition grid.h:291
real_type x1() const
Equivalent to p(1)
Definition grid.h:288
unsigned Ny() const
Equivalent to N(1)
Definition grid.h:337