6#define M_PI 3.14159265358979323846
27template <
class T =
double>
35 PLUS( T value): x_(value){}
43template<
class T =
double >
53template <
class T =
double>
64template <
class T =
double>
81template <
class T =
double>
91template <
class T =
double>
98template <
class T =
double>
112template <
class T =
double>
119template <
class T =
double>
124 T absx =
x>0 ?
x : -
x;
125 T absy =
y>0 ?
y : -
y;
126 return absx > absy ? absx : absy;
131template <
class T =
double>
136 T absx =
x<0 ? -
x :
x;
137 T absy =
y<0 ? -
y :
y;
138 return absx < absy ? absx : absy;
150template <
class T =
double>
154 if (
x >= 0.0)
return x;
165template <
class T=
double>
177 return (fmod(
x,m_m) < 0 ) ? (m_m + fmod(
x,m_m)) : fmod(
x,m_m);
229 if(
x > 1e100 ||
x < -1e100)
235 if( !std::isfinite(
x))
237 if(
x > 1e100 ||
x < -1e100)
262 if( x1 > 0 && x2 > 0)
264 else if( x1 < 0 && x2 < 0)
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);
283 return this->
operator()( this->
operator()( x1, x2), x3);
302 return 2.*x1*x2/(x1+x2);
348template<
class Limiter>
356 return +hm*m_l( g0, gm);
358 return -hp*m_l( gp, g0);
373template<
class Limiter>
380 return v*m_s(v,gm,g0,gp,hm,hp);
407 Iris(
double psi_min,
double psi_max ):
408 m_psimin(psi_min), m_psimax(psi_max) { }
412 if( psi > m_psimax)
return 0.;
413 if( psi < m_psimin)
return 0.;
417 double m_psimin, m_psimax;
433 if( psi > psimax_)
return 0.;
453 if( psi > psimax_)
return psimax_;
479 m_xb(xb), m_s(sign){ }
484 if( (
x < m_xb && m_s == 1) || (
x > m_xb && m_s == -1))
return 0.;
498 Distance(
double x0,
double y0): m_x0(x0), m_y0(y0){}
501 return sqrt( (
x-m_x0)*(
x-m_x0) + (
y-m_y0)*(
y-m_y0));
513 Line(
double x0,
double y0,
double x1,
double y1) :
514 m_x0(x0), m_y0(y0), m_x1(x1), m_y1(y1){}
516 return m_y1*(
x-m_x0)/(m_x1-m_x0) + m_y0*(
x-m_x1)/(m_x0-m_x1);
519 double m_x0, m_y0, m_x1, m_y1;
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");
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) );
638 double m_x0, m_y0, m_sigma_x, m_sigma_y, m_amp;
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");
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.) );
678 double xbar = (
x-x0_)/sigmaX_;
679 double ybar = (
y-y0_)/sigmaY_;
680 if( xbar*xbar + ybar*ybar < 1.)
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;
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;
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;
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;
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_
720 double x0_, y0_, sigmaX_, sigmaY_, amp_;
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");
748 double xbar = (
x-x0_)/sigmaX_;
750 return amp_*exp( 1. + 1./( xbar*xbar -1.) );
755 double xbar = (
x-x0_)/sigmaX_;
761 double x0_, sigmaX_, amp_;
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");
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) );
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) );
828 double m_x0, m_y0, m_z0, m_sigma_x, m_sigma_y, m_sigma_z, m_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");
853 return m_amp* exp( -((
x-m_x0)*(
x-m_x0)/2./m_sigma_x/m_sigma_x ));
866 double m_x0, m_sigma_x, m_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");
901 return m_amp*exp( -((
y-m_y0)*(
y-m_y0)/2./m_sigma_y/m_sigma_y) );
904 double m_y0, m_sigma_y, m_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");
939 return m_amp*exp( -((
z-m_z0)*(
z-m_z0)/2./m_sigma_z/m_sigma_z) );
954 return m_amp*exp( -((
z-m_z0)*(
z-m_z0)/2./m_sigma_z/m_sigma_z) );
957 double m_z0, m_sigma_z, m_amp;
972 IslandXY(
double lambda,
double eps):lambda_(lambda), eps_(eps){
973 assert( lambda != 0 &&
"Lambda parameter in IslandXY must not be zero!");
983 double operator()(
double x,
double y)
const{
return lambda_*log(cosh(
x/lambda_)+eps_*cos(
y/lambda_));}
1000 SinXSinY(
double amp,
double bamp,
double kx,
double ky):amp_(
amp), bamp_(bamp),kx_(kx),ky_(ky){}
1010 double operator()(
double x,
double y)
const{
return bamp_+amp_*sin(
x*kx_)*sin(
y*ky_);}
1012 double amp_,bamp_,kx_,ky_;
1028 CosXCosY(
double amp,
double bamp,
double kx,
double ky):amp_(
amp), bamp_(bamp),kx_(kx),ky_(ky){}
1038 double operator()(
double x,
double y)
const{
return bamp_+amp_*cos(
x*kx_)*cos(
y*ky_);}
1040 double amp_,bamp_,kx_,ky_;
1056 SinXCosY(
double amp,
double bamp,
double kx,
double ky):amp_(
amp), bamp_(bamp),kx_(kx),ky_(ky){}
1066 double operator()(
double x,
double y)
const{
return bamp_+amp_*sin(
x*kx_)*cos(
y*ky_);}
1068 double amp_,bamp_,kx_,ky_;
1083 SinX(
double amp,
double bamp,
double kx):amp_(
amp), bamp_(bamp),kx_(kx){}
1091 double amp_,bamp_,kx_;
1106 SinY(
double amp,
double bamp,
double ky):amp_(
amp), bamp_(bamp),ky_(ky){}
1110 double amp_,bamp_,ky_;
1125 CosY(
double amp,
double bamp,
double ky):amp_(
amp), bamp_(bamp),ky_(ky){}
1129 double amp_,bamp_,ky_;
1145 double operator()(
double x)
const{
return m_amp/cosh(
x*m_kx)/cosh(
x*m_kx);}
1166 SinProfX(
double amp,
double bamp,
double kx):m_amp(
amp), m_bamp(bamp),m_kx(kx){}
1168 double operator()(
double x)
const{
return m_bamp+m_amp*(1.-sin(
x*m_kx));}
1174 double m_amp, m_bamp, m_kx;
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!");
1199 double m_amp, m_bamp, m_ln;
1217 m_psimax(psimax), m_alpha(
alpha) {
1218 assert(
alpha!= 0 &&
"Damping width in GaussianDamping must not be zero");
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);
1228 double m_psimax, m_alpha;
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!");
1255 return profa_*0.5*(1.+s_*tanh((
x-xb_)/w_))+bga_;
1291 x0(xb), a(a), m_s(sign){
1292 assert( a!=0 &&
"PolynomialHeaviside width must not be zero");
1297 if( m_s == -1)
x = 2*x0-
x;
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);
1343 m_hl( xl, al, +1), m_hr( xr, ar, -1) {
1344 assert( xl < xr &&
"left boundary must be left of right boundary");
1349 return m_hl(
x)*m_hr(
x);
1381 x0(xb), a(a), m_s(sign){
1382 assert( a!=0 &&
"IPolynomialHeaviside width must not be zero");
1387 if( m_s == -1)
x = 2*x0-
x;
1389 if (
x < x0-a) result = x0;
1390 else if (
x > x0+a) result =
x;
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;
1434 assert( a!=0 &&
"DPolynomialHeaviside width must not be zero");
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);
1480 m_alpha(
alpha), m_etac(eta_c), m_s(order), m_n(
n) {}
1483 double eta = (double)i/(
double)(m_n-1);
1484 if( m_n == 1) eta = 0.;
1488 return exp( -m_alpha*pow( (eta-m_etac)/(1.-m_etac), 2*m_s));
1492 double m_alpha, m_etac;
1520 Lamb(
double x0,
double y0,
double R,
double U):R_(R), U_(U), x0_(x0), y0_(y0)
1522 gamma_ = 3.83170597020751231561;
1542 double radius = sqrt( (
x-x0_)*(
x-x0_) + (
y-y0_)*(
y-y0_));
1543 double theta = atan2( (
y-y0_),(
x-x0_));
1547 return 2.*lambda_*U_*_j1(lambda_*radius)/j_*cos( theta);
1549 return 2.*lambda_*U_*j1( lambda_*radius)/j_*cos( theta);
1570 double R_, U_, x0_, y0_, lambda_, gamma_, j_;
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;
1612 b_[0] = 0.03827327723;
1613 b_[1] = 0.07071067810 ;
1614 b_[2] = 0.07071067810 ;
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;
1644 - R_*
beta*
beta/g_[s_]/g_[s_] *_j1(g_[s_]*r/R_)/_j1(g_[s_])
1646 - R_ *
beta*
beta/g_[s_]/g_[s_] * j1(g_[s_]*r/R_)/ j1(g_[s_])
1649 return u_d * R_* bessk1(
beta*r/R_)/bessk1(
beta)*cos(theta)/norm;
1677 double bessk1(
double x)
const
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)))))));
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)))))));
1698 double bessi1(
double x)
const
1702 if ((ax=fabs(x)) < 3.75)
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))))));
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));
1717 return x < 0.0 ? -ans : ans;
1721 double R_, b_[3], u_d;
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)
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);
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.;
1773 std::minstd_rand generator;
1774 std::normal_distribution<double> ndistribution( 0.0, 1.0);
1775 std::uniform_real_distribution<double> udistribution(0.0,tpi);
1776 for (
unsigned j=1;j<=N_kZ_;j++)
1778 double kZ2=tpi2*(j-N_kZh)*(j-N_kZh)/(N_kZ2);
1779 for (
unsigned i=1;i<=N_kR_;i++)
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);
1784 sqEkvec[
z]=pow(kvec[
z]*4.*k0/(kvec[
z]+k0)/(kvec[
z]+k0),gamma_/2.);
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]);
1820 double f, kappa, RR, ZZ;
1824 for (
unsigned j=0;j<N_kZ_;j++)
1826 for (
unsigned i=0;i<N_kR_;i++)
1829 kappa= RR*unif1[
z]+ZZ*unif2[
z];
1830 f+= sqEkvec[
z]*alpha[
z]*cos(kvec[
z]*kappa+theta[
z]);
1833 return amp_*norm_*f;
1865 for (
unsigned j=0;j<N_kZ_;j++)
1867 for (
unsigned i=0;i<N_kR_;i++)
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]);
1874 return amp_*norm_*f;
1877 unsigned N_kR_,N_kZ_;
1878 double R_min_, Z_min_;
1879 double gamma_, L_E_;
1882 std::vector<double> kvec;
1883 std::vector<double> sqEkvec;
1884 std::vector<double> unif1, unif2, normal1,normal2,alpha,theta;
1904 Horner2d( std::vector<double> c,
unsigned M,
unsigned N): m_c(c), m_M(M), m_N(N){}
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);
1913 double horner(
const double * c,
unsigned M,
double x)
const
1916 for(
unsigned i=0; i<M-1; i++)
1920 std::vector<double> m_c;
1929template <
class container = thrust::host_vector<
double> >
1940 binwidth_(g1d_.h()),
1943 for (
unsigned j=0;j<in_.size();j++)
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));
1951 unsigned Ampmax = (unsigned)thrust::reduce( count_.begin(), count_.end(),0., thrust::maximum<double>() );
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));
1979 const std::vector<double> in_;
1988template <
class container = thrust::host_vector<
double> >
2001 binwidthx_(g2d_.hx()),
2002 binwidthy_(g2d_.hy()),
2006 for (
unsigned j=0;j<iny_.size();j++)
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));
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.;
2019 unsigned Ampmax = (unsigned)thrust::reduce( count_.begin(), count_.end(),0.,thrust::maximum<double>() );
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 ];
2045 const std::vector<double> inx_,iny_;
2046 double binwidthx_,binwidthy_;
2066 WallDistance( std::vector<double> vertical, std::vector<double> horizontal) :
2067 m_vertical(vertical), m_horizontal( horizontal) {}
2074 m_horizontal({walls.y0(), walls.y1()}){}
2080 std::vector<double> dist( 1, 1e100);
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());
2088 std::vector<double> m_vertical;
2089 std::vector<double> m_horizontal;
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
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
@ 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
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
This is the namespace for all functions and classes defined and used by the discontinuous Galerkin li...
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
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
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
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
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
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
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
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