functors.h

Go to the documentation of this file.

#pragma once
 
#include <cmath>
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
#include <vector>
#include <functional>
#include <thrust/random/linear_congruential_engine.h>
#include <thrust/random/normal_distribution.h>
#include <thrust/random/uniform_real_distribution.h>
#include "blas1.h"
#include "topology/grid.h"
#include "topology/evaluation.h"
#include "topology/functions.h"
namespace dg
{
 
//Everything that is quite basic and simple
 
 
template <class T = double>
struct PLUS
{
    PLUS( T value): x_(value){}
    DG_DEVICE
    T operator()( T x)const{ return x + x_;}
    private:
    T x_;
};
 
template< class T = double >
struct EXP
{
    DG_DEVICE T operator() ( T x) const
    {
        return exp(x);
    }
};
 
template < class T = double>
struct LN
{
    DG_DEVICE T operator() (const T& x) const
    {
        return log(x);
    }
 
};
 
template < class T = double>
struct SQRT
{
    DG_DEVICE T operator() (T x) const
    {
        return sqrt(x);
    }
};
 
struct Square
{
    template<class T>
    DG_DEVICE T operator()( T x) const{ return x*x;}
};
 
template < class T = double>
struct InvSqrt
{
    DG_DEVICE T operator() (T x) const
    {
        return 1./sqrt(x);
    }
};
 
template <class T = double>
struct INVERT
{
    DG_DEVICE T operator()( T x)const{ return 1./x;}
};
 
template <class T = double>
struct ABS
{
    DG_DEVICE T operator()(T x)const{ return fabs(x);}
};
 
template <class T = double>
struct Sign
{
    DG_DEVICE T operator()(T x)const{ return (T(0) < x) - (x < T(0));}
};
 
template <class T = double>
struct AbsMax
{
    DG_DEVICE T operator() ( T x, T y) const
    {
        T absx = x>0 ? x : -x;
        T absy = y>0 ? y : -y;
        return absx > absy ? absx : absy;
    }
};
 
template <class T = double>
struct AbsMin
{
    DG_DEVICE T operator() (T x, T y) const
    {
        T absx = x<0 ? -x : x;
        T absy = y<0 ? -y : y;
        return absx < absy ? absx : absy;
    }
};
 
template <class T = double>
struct POSVALUE
{
    DG_DEVICE T operator()( T x)const{
        if (x >= 0.0) return x;
        return 0.0;
    }
};
 
template <class T= double>
struct MOD
{
    MOD( T m): m_m(m){}
 
    DG_DEVICE
    T operator()( T x)const{
        return (fmod(x,m_m) < 0 ) ? (m_m + fmod(x,m_m)) : fmod(x,m_m);
    }
    private:
    T m_m;
 
};
 
template <class T>
struct ISNFINITE
{
#ifdef __CUDACC__
    DG_DEVICE bool operator()(T x){ return !isfinite(x);}
#else
    bool operator()( T x){ return !std::isfinite(x);}
#endif
};
 
template <class T>
struct ISNSANE
{
#ifdef __CUDACC__
    DG_DEVICE bool operator()(T x){
        if( !isfinite(x))
            return true;
        if( x > 1e100 || x < -1e100)
            return true;
        return false;
    }
#else
    bool operator()( T x){
        if( !std::isfinite(x))
            return true;
        if( x > 1e100 || x < -1e100)
            return true;
        return false;
    }
#endif
};
 
struct MinMod
{
#ifdef __CUDACC__
    template < class T>
    DG_DEVICE T operator()( T x1, T x2) const
    {
        if( x1 > 0 && x2 > 0)
            return min(x1,x2);
        else if( x1 < 0 && x2 < 0)
            return max(x1,x2);
        return 0.;
    }
#else
    template < class T>
    T operator()( T x1, T x2) const
    {
        if( x1 > 0 && x2 > 0)
            return std::min(x1,x2);
        else if( x1 < 0 && x2 < 0)
            return std::max(x1,x2);
        return 0.;
    }
#endif
    template<class T>
    DG_DEVICE T operator() ( T x1, T x2, T x3)const
    {
        return this-> operator()( this-> operator()( x1, x2), x3);
    }
};
 
template<class T>
struct UniformRealDistribution
{
    UniformRealDistribution( T a, T b): m_rng(), m_dist(a,b){}
 
    template< class ...Ts>
DG_DEVICE T operator()( Ts... )
    {
        return m_dist(m_rng); //what happens on GPU? Does every thread get the same?
    }
    private:
    thrust::minstd_rand m_rng;
    thrust::uniform_real_distribution<T> m_dist;
};
 
template<class T>
using RandomNumbers = UniformRealDistribution<T>;
 
struct VanLeer
{
    template<class T>
    DG_DEVICE T operator()( T x1, T x2) const
    {
        if( x1*x2 <= 0)
            return 0.;
        return 2.*x1*x2/(x1+x2);
    }
};
 
struct Upwind
{
    template<class T>
    DG_DEVICE T operator()( T velocity, T backward, T forward) const{
        if( velocity >= 0)
            return backward;
        else
            return forward;
    }
};
 
struct UpwindProduct
{
    template<class T>
    DG_DEVICE T operator()( T velocity, T backward, T forward)const{
        return velocity*m_up(velocity, backward, forward);
    }
    private:
    Upwind m_up;
};
 
template<class Limiter>
struct SlopeLimiter
{
    SlopeLimiter() {}
    SlopeLimiter( Limiter l ) : m_l( l){}
    template<class T>
    DG_DEVICE T operator()( T v, T gm, T g0, T gp, T hm, T hp ) const{
        if( v >= 0)
            return +hm*m_l( g0, gm);
        else
            return -hp*m_l( gp, g0);
    }
    private:
    Limiter m_l;
};
 
template<class Limiter>
struct SlopeLimiterProduct
{
    SlopeLimiterProduct() {}
    SlopeLimiterProduct( Limiter l ) : m_s( l){}
    template<class T>
    DG_DEVICE T operator()( T v, T gm, T g0, T gp, T hm, T hp ) const{
        return v*m_s(v,gm,g0,gp,hm,hp);
    }
    private:
    SlopeLimiter<Limiter> m_s;
};
 
 
 
//
 
struct Iris
{
    Iris( double psi_min, double psi_max ):
        m_psimin(psi_min), m_psimax(psi_max) { }
    DG_DEVICE
    double operator()(double psi)const
    {
        if( psi > m_psimax) return 0.;
        if( psi < m_psimin) return 0.;
        return 1.;
    }
    private:
    double m_psimin, m_psimax;
};
struct Pupil
{
    Pupil( double psimax):
        psimax_(psimax) { }
    DG_DEVICE
    double operator()(double psi)const
    {
        if( psi > psimax_) return 0.;
        return 1.;
    }
    private:
    double psimax_;
};
struct PsiPupil
{
    PsiPupil(double psimax):
        psimax_(psimax){ }
    DG_DEVICE
    double operator()(double psi)const
    {
        if( psi > psimax_) return psimax_;
        return  psi;
    }
    private:
    double psimax_;
};
struct Heaviside
{
 
    Heaviside( double xb, int sign = +1):
        m_xb(xb), m_s(sign){ }
 
    DG_DEVICE
    double operator()(double x)const
    {
        if( (x < m_xb && m_s == 1) || (x > m_xb && m_s == -1)) return 0.;
        return 1.;
    }
    private:
    double m_xb;
    int m_s;
};
 
 
struct Distance
{
    Distance( double x0, double y0): m_x0(x0), m_y0(y0){}
    DG_DEVICE
    double operator()(double x, double y){
        return sqrt( (x-m_x0)*(x-m_x0) + (y-m_y0)*(y-m_y0));
    }
    private:
    double m_x0, m_y0;
};
struct Line{
    Line(double x0, double y0, double x1, double y1) :
        m_x0(x0), m_y0(y0), m_x1(x1), m_y1(y1){}
    double operator()(double x){
        return m_y1*(x-m_x0)/(m_x1-m_x0) + m_y0*(x-m_x1)/(m_x0-m_x1);
    }
    private:
    double m_x0, m_y0, m_x1, m_y1;
};
 
struct LinearX
{
    LinearX( double a, double b):a_(a), b_(b){}
    DG_DEVICE
    double operator()(double x)const{ return a_*x+b_;}
    DG_DEVICE
    double operator()( double x, double)const{ return this->operator()(x);}
    DG_DEVICE
    double operator()( double x, double, double)const{ return this->operator()(x);}
   private:
    double a_,b_;
};
struct LinearY
{
    LinearY( double a, double b):a_(a), b_(b){}
    DG_DEVICE
    double operator()( double, double y, double)const { return a_*y+b_;}
    DG_DEVICE
    double operator()( double, double y)const{ return a_*y+b_;}
  private:
    double a_,b_;
};
struct LinearZ
{
    LinearZ( double a, double b):a_(a), b_(b){}
    DG_DEVICE
    double operator()( double, double, double z)const{ return a_*z+b_;}
  private:
    double a_,b_;
};
struct Gaussian
{
    Gaussian( double x0, double y0, double sigma_x, double sigma_y, double amp)
        : m_x0(x0), m_y0(y0), m_sigma_x(sigma_x), m_sigma_y(sigma_y), m_amp(amp){
            assert( m_sigma_x != 0  &&  "sigma_x must not be 0 in Gaussian");
            assert( m_sigma_y != 0  &&  "sigma_y must not be 0 in Gaussian");
    }
    DG_DEVICE
    double operator()(double x, double y) const
    {
        return  m_amp*
                   exp( -((x-m_x0)*(x-m_x0)/2./m_sigma_x/m_sigma_x +
                          (y-m_y0)*(y-m_y0)/2./m_sigma_y/m_sigma_y) );
    }
    DG_DEVICE
    double operator()(double x, double y, double) const
    {
        return  this->operator()(x,y);
    }
  private:
    double  m_x0, m_y0, m_sigma_x, m_sigma_y, m_amp;
 
};
 
struct Cauchy
{
    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){
        assert( sigma_x != 0  &&  "sigma_x must be !=0 in Cauchy");
        assert( sigma_y != 0  &&  "sigma_y must be !=0 in Cauchy");
    }
    DG_DEVICE
    double operator()(double x, double y )const{
        double xbar = (x-x0_)/sigmaX_;
        double ybar = (y-y0_)/sigmaY_;
        if( xbar*xbar + ybar*ybar < 1.)
            return amp_*exp( 1. +  1./( xbar*xbar + ybar*ybar -1.) );
        return 0.;
    }
    bool inside( double x, double y)const
    {
        double xbar = (x-x0_)/sigmaX_;
        double ybar = (y-y0_)/sigmaY_;
        if( xbar*xbar + ybar*ybar < 1.)
            return true;
        return false;
    }
 
    double dx( double x, double y )const{
        double xbar = (x-x0_)/sigmaX_;
        double ybar = (y-y0_)/sigmaY_;
        double temp = sigmaX_*(xbar*xbar + ybar*ybar  - 1.);
        return -2.*(x-x0_)*this->operator()(x,y)/temp/temp;
    }
    double dxx( double x, double y)const{
        double temp = sigmaY_*sigmaY_*(x-x0_)*(x-x0_) + sigmaX_*sigmaX_*((y-y0_)*(y-y0_) - sigmaY_*sigmaY_);
        double bracket = sigmaX_*sigmaX_*((y-y0_)*(y-y0_)-sigmaY_*sigmaY_)*sigmaX_*sigmaX_*((y-y0_)*(y-y0_)-sigmaY_*sigmaY_)
            -3.*sigmaY_*sigmaY_*sigmaY_*sigmaY_*(x-x0_)*(x-x0_)*(x-x0_)*(x-x0_)
            -2.*sigmaY_*sigmaY_*sigmaX_*sigmaX_*(x-x0_)*(x-x0_)*(y-y0_)*(y-y0_);
        return -2.*sigmaX_*sigmaX_*sigmaY_*sigmaY_*sigmaY_*sigmaY_*this->operator()(x,y)*bracket/temp/temp/temp/temp;
    }
    double dy( double x, double y)const{
        double xbar = (x-x0_)/sigmaX_;
        double ybar = (y-y0_)/sigmaY_;
        double temp = sigmaY_*(xbar*xbar + ybar*ybar  - 1.);
        return -2.*(y-y0_)*this->operator()(x,y)/temp/temp;
    }
    double dyy( double x, double y)const{
        double temp = sigmaX_*sigmaX_*(y-y0_)*(y-y0_) + sigmaY_*sigmaY_*((x-x0_)*(x-x0_) - sigmaX_*sigmaX_);
        double bracket = sigmaY_*sigmaY_*((x-x0_)*(x-x0_)-sigmaX_*sigmaX_)*sigmaY_*sigmaY_*((x-x0_)*(x-x0_)-sigmaX_*sigmaX_)
            -3.*sigmaX_*sigmaX_*sigmaX_*sigmaX_*(y-y0_)*(y-y0_)*(y-y0_)*(y-y0_)
            -2.*sigmaX_*sigmaX_*sigmaY_*sigmaY_*(y-y0_)*(y-y0_)*(x-x0_)*(x-x0_);
        return -2.*sigmaY_*sigmaY_*sigmaX_*sigmaX_*sigmaX_*sigmaX_*this->operator()(x,y)*bracket/temp/temp/temp/temp;
    }
    double dxy( double x, double y )const{
        double xbar = (x-x0_)/sigmaX_;
        double ybar = (y-y0_)/sigmaY_;
        double temp = (xbar*xbar + ybar*ybar  - 1.);
        return 8.*xbar*ybar*this->operator()(x,y)/temp/temp/temp/sigmaX_/sigmaY_
            + 4.*xbar*ybar*this->operator()(x,y)/temp/temp/temp/temp/sigmaX_/sigmaY_
;
    }
    private:
    double x0_, y0_, sigmaX_, sigmaY_, amp_;
};
 
struct CauchyX
{
    CauchyX( double x0,  double sigma_x, double amp): x0_(x0), sigmaX_(sigma_x),  amp_(amp){
        assert( sigma_x != 0  &&  "sigma_x must be !=0 in Cauchy");
    }
    DG_DEVICE
    double operator()(double x, double )const{
        double xbar = (x-x0_)/sigmaX_;
        if( xbar*xbar  < 1.)
            return amp_*exp( 1. +  1./( xbar*xbar -1.) );
        return 0.;
    }
    bool inside( double x, double)const
    {
        double xbar = (x-x0_)/sigmaX_;
        if( xbar*xbar < 1.)
            return true;
        return false;
    }
    private:
    double x0_, sigmaX_,  amp_;
};
 
struct Gaussian3d
{
    Gaussian3d( double x0, double y0, double z0, double sigma_x, double sigma_y, double sigma_z, double amp)
        : 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){
            assert( m_sigma_x != 0  &&  "sigma_x must be !=0 in Gaussian3d");
            assert( m_sigma_y != 0  &&  "sigma_y must be !=0 in Gaussian3d");
            assert( m_sigma_z != 0  &&  "sigma_z must be !=0 in Gaussian3d");
    }
    DG_DEVICE
    double operator()(double x, double y) const
    {
        return  m_amp*
                   exp( -((x-m_x0)*(x-m_x0)/2./m_sigma_x/m_sigma_x +
                          (y-m_y0)*(y-m_y0)/2./m_sigma_y/m_sigma_y) );
    }
    DG_DEVICE
    double operator()(double x, double y, double z) const
    {
        return  m_amp*
                exp( -((x-m_x0)*(x-m_x0)/2./m_sigma_x/m_sigma_x +
                       (z-m_z0)*(z-m_z0)/2./m_sigma_z/m_sigma_z +
                       (y-m_y0)*(y-m_y0)/2./m_sigma_y/m_sigma_y) );
    }
  private:
    double  m_x0, m_y0, m_z0, m_sigma_x, m_sigma_y, m_sigma_z, m_amp;
 
};
struct GaussianX
{
    GaussianX( double x0, double sigma_x, double amp)
        :m_x0(x0), m_sigma_x(sigma_x), m_amp(amp){
            assert( m_sigma_x != 0  &&  "sigma_x must be !=0 in GaussianX");
    }
    DG_DEVICE
    double operator()(double x) const
    {
        return  m_amp* exp( -((x-m_x0)*(x-m_x0)/2./m_sigma_x/m_sigma_x ));
    }
    DG_DEVICE
    double operator()(double x, double) const
    {
        return this->operator()(x);
    }
    DG_DEVICE
    double operator()(double x, double, double) const
    {
        return this->operator()(x);
    }
  private:
    double  m_x0, m_sigma_x, m_amp;
 
};
struct GaussianY
{
    GaussianY( double y0, double sigma_y, double amp)
        : m_y0(y0), m_sigma_y(sigma_y), m_amp(amp){
            assert( m_sigma_y != 0  &&  "sigma_x must be !=0 in GaussianY");
    }
    DG_DEVICE
    double operator()(double, double y) const
    {
        return  m_amp*exp( -((y-m_y0)*(y-m_y0)/2./m_sigma_y/m_sigma_y) );
    }
  private:
    double  m_y0, m_sigma_y, m_amp;
 
};
struct GaussianZ
{
    GaussianZ( double z0, double sigma_z, double amp)
        : m_z0(z0), m_sigma_z(sigma_z), m_amp(amp){
            assert( m_sigma_z != 0  &&  "sigma_z must be !=0 in GaussianZ");
    }
    DG_DEVICE
    double operator()( double z) const
    {
        return  m_amp*exp( -((z-m_z0)*(z-m_z0)/2./m_sigma_z/m_sigma_z) );
    }
    DG_DEVICE
    double operator()(double, double, double z) const
    {
        return  m_amp*exp( -((z-m_z0)*(z-m_z0)/2./m_sigma_z/m_sigma_z) );
    }
  private:
    double  m_z0, m_sigma_z, m_amp;
 
};
struct IslandXY
{
     IslandXY( double lambda, double eps):lambda_(lambda), eps_(eps){
         assert( lambda != 0 && "Lambda parameter in IslandXY must not be zero!");
     }
    DG_DEVICE
    double operator()( double x, double y)const{ return lambda_*log(cosh(x/lambda_)+eps_*cos(y/lambda_));}
  private:
    double lambda_,eps_;
};
struct SinXSinY
{
    SinXSinY( double amp, double bamp, double kx, double ky):amp_(amp), bamp_(bamp),kx_(kx),ky_(ky){}
    DG_DEVICE
    double operator()( double x, double y)const{ return bamp_+amp_*sin(x*kx_)*sin(y*ky_);}
  private:
    double amp_,bamp_,kx_,ky_;
};
struct CosXCosY
{
    CosXCosY( double amp, double bamp, double kx, double ky):amp_(amp), bamp_(bamp),kx_(kx),ky_(ky){}
    DG_DEVICE
    double operator()( double x, double y)const{ return bamp_+amp_*cos(x*kx_)*cos(y*ky_);}
  private:
    double amp_,bamp_,kx_,ky_;
};
struct SinXCosY
{
    SinXCosY( double amp, double bamp, double kx, double ky):amp_(amp), bamp_(bamp),kx_(kx),ky_(ky){}
    DG_DEVICE
    double operator()( double x, double y)const{ return bamp_+amp_*sin(x*kx_)*cos(y*ky_);}
  private:
    double amp_,bamp_,kx_,ky_;
};
struct SinX
{
    SinX( double amp, double bamp, double kx):amp_(amp), bamp_(bamp),kx_(kx){}
    DG_DEVICE
    double operator()( double x)const{ return bamp_+amp_*sin(x*kx_);}
    DG_DEVICE
    double operator()( double x, double)const{ return this->operator()(x);}
    DG_DEVICE
    double operator()( double x, double, double)const{ return this->operator()(x);}
  private:
    double amp_,bamp_,kx_;
};
struct SinY
{
    SinY( double amp, double bamp, double ky):amp_(amp), bamp_(bamp),ky_(ky){}
    DG_DEVICE
    double operator()( double, double y)const{ return bamp_+amp_*sin(y*ky_);}
  private:
    double amp_,bamp_,ky_;
};
struct CosY
{
    CosY( double amp, double bamp, double ky):amp_(amp), bamp_(bamp),ky_(ky){}
    DG_DEVICE
    double operator()( double, double y)const{ return bamp_+amp_*cos(y*ky_);}
  private:
    double amp_,bamp_,ky_;
};
struct InvCoshXsq
{
    InvCoshXsq( double amp, double kx):m_amp(amp), m_kx(kx){}
    DG_DEVICE
    double operator()( double x)const{ return m_amp/cosh(x*m_kx)/cosh(x*m_kx);}
    DG_DEVICE
    double operator()( double x, double)const{ return this->operator()(x);}
    DG_DEVICE
    double operator()( double x, double, double)const{ return this->operator()(x);}
  private:
    double m_amp, m_kx;
};
struct SinProfX
{
    SinProfX( double amp, double bamp, double kx):m_amp(amp), m_bamp(bamp),m_kx(kx){}
    DG_DEVICE
    double operator()( double x)const{ return m_bamp+m_amp*(1.-sin(x*m_kx));}
    DG_DEVICE
    double operator()( double x, double)const{ return this->operator()(x);}
    DG_DEVICE
    double operator()( double x, double, double)const{ return this->operator()(x);}
  private:
    double m_amp, m_bamp, m_kx;
};
struct ExpProfX
{
    ExpProfX( double amp, double bamp, double ln):m_amp(amp),m_bamp(bamp),m_ln(ln){
        assert( ln!=0 && "ln parameter must be != 0 in ExpProfX!");
    }
    DG_DEVICE
    double operator()( double x)const{ return m_bamp+m_amp*exp(-x/m_ln);}
    DG_DEVICE
    double operator()( double x, double)const{ return this->operator()(x);}
    DG_DEVICE
    double operator()( double x, double, double)const{ return this->operator()(x);}
  private:
    double m_amp, m_bamp, m_ln;
};
 
 
struct GaussianDamping
{
    GaussianDamping( double psimax, double alpha):
        m_psimax(psimax), m_alpha(alpha) {
            assert( alpha!= 0 && "Damping width in GaussianDamping must not be zero");
        }
    DG_DEVICE
    double operator()(double psi)const
    {
        if( psi > m_psimax + 4.*m_alpha) return 0.;
        if( psi < m_psimax) return 1.;
        return exp( -( psi-m_psimax)*( psi-m_psimax)/2./m_alpha/m_alpha);
    }
    private:
    double m_psimax, m_alpha;
};
struct TanhProfX {
    TanhProfX(double xb, double width, int sign =1,double bgamp = 0.,
        double profamp = 1.) :
        xb_(xb),w_(width), s_(sign),bga_(bgamp),profa_(profamp)  {
            assert( width != 0&& "Width in TanhProfX must not be zero!");
    }
    DG_DEVICE
    double operator() (double x)const
    {
        return profa_*0.5*(1.+s_*tanh((x-xb_)/w_))+bga_;
    }
    DG_DEVICE
    double operator()( double x, double)const{ return this->operator()(x);}
    DG_DEVICE
    double operator()( double x, double, double)const{ return this->operator()(x);}
    private:
    double xb_;
    double w_;
    int s_;
    double bga_;
    double profa_;
};
 
struct PolynomialHeaviside {
    PolynomialHeaviside(double xb, double a, int sign = +1) :
        x0(xb), a(a), m_s(sign){
            assert( a!=0 && "PolynomialHeaviside width must not be zero");
    }
    DG_DEVICE
    double operator() (double x)const
    {
        if( m_s == -1) x = 2*x0-x; //mirror
        if ( x < x0-a) return 0;
        if ( x > x0+a) return 1;
        return ((16.*a*a*a - 29.*a*a*(x - x0)
               + 20.*a*(x - x0)*(x - x0)
               - 5.*(x - x0)*(x-x0)*(x-x0))
               *(a + x - x0)*(a + x - x0)
               *(a + x - x0)*(a + x - x0))/(32.*a*a*a * a*a*a*a);
    }
    DG_DEVICE
    double operator()( double x, double)const{ return this->operator()(x);}
    DG_DEVICE
    double operator()( double x, double, double)const{ return this->operator()(x);}
    private:
    double x0, a;
    int m_s;
};
 
struct PolynomialRectangle {
    PolynomialRectangle(double xl, double al, double xr, double ar) :
        m_hl( xl, al, +1), m_hr( xr, ar, -1) {
        assert( xl < xr && "left boundary must be left of right boundary");
    }
    DG_DEVICE
    double operator() (double x)const
    {
        return m_hl(x)*m_hr(x);
    }
    DG_DEVICE
    double operator()( double x, double)const{ return this->operator()(x);}
    DG_DEVICE
    double operator()( double x, double, double)const{ return this->operator()(x);}
    private:
    PolynomialHeaviside m_hl, m_hr;
};
 
struct IPolynomialHeaviside {
    IPolynomialHeaviside(double xb, double a, int sign = +1) :
        x0(xb), a(a), m_s(sign){
            assert( a!=0 && "IPolynomialHeaviside width must not be zero");
        }
    DG_DEVICE
    double operator() (double x)const
    {
        if( m_s == -1) x = 2*x0-x; //mirror
        double result;
        if ( x < x0-a) result =  x0;
        else if ( x > x0+a) result =  x;
        else
            result =  x0 + ((35.* a*a*a - 47.* a*a*(x - x0) + 25.*a*(x - x0)*(x-x0)
                - 5.*(x - x0)*(x-x0)*(x-x0))
                *(a+x-x0)*(a+x-x0)*(a+x-x0)*(a+x-x0)*(a+x-x0))
            /(256.*a*a*a * a*a*a*a);
        if ( m_s == +1) return result;
        return 2*x0 - result;
 
    }
    DG_DEVICE
    double operator()( double x, double)const{ return this->operator()(x);}
    DG_DEVICE
    double operator()( double x, double, double)const{ return this->operator()(x);}
    private:
    double x0, a;
    int m_s;
};
 
struct DPolynomialHeaviside {
    DPolynomialHeaviside(double xb, double a, int = +1) :
        x0(xb), a(a){
            assert( a!=0 && "DPolynomialHeaviside width must not be zero");
    }
    DG_DEVICE
    double operator() (double x)const
    {
        if ( (x < x0-a) || (x > x0+a)) return 0;
        return (35.*(a+x-x0)*(a+x-x0)*(a+x-x0)*(a-x+x0)*(a-x+x0)*(a-x+x0))
            /(32.*a*a*a * a*a*a*a);
    }
    DG_DEVICE
    double operator()( double x, double)const{ return this->operator()(x);}
    DG_DEVICE
    double operator()( double x, double, double)const{ return this->operator()(x);}
    private:
    double x0, a;
};
 
struct ExponentialFilter
{
    ExponentialFilter( double alpha, double eta_c, unsigned order, unsigned n):
        m_alpha(alpha), m_etac(eta_c), m_s(order), m_n(n) {}
    double operator()( unsigned i) const
    {
        double eta = (double)i/(double)(m_n-1);
        if( m_n == 1) eta = 0.;
        if( eta < m_etac)
            return 1.;
        if( eta <= 1.)
            return exp( -m_alpha*pow( (eta-m_etac)/(1.-m_etac), 2*m_s));
        return 0;
    }
    private:
    double m_alpha, m_etac;
    unsigned m_s, m_n;
};
 
struct Lamb
{
    Lamb(  double x0, double y0, double R, double U):R_(R), U_(U), x0_(x0), y0_(y0)
    {
        gamma_ = 3.83170597020751231561;
        lambda_ = gamma_/R;
#ifdef _MSC_VER
        j_ = _j0(gamma_);
#else
        j_ = j0( gamma_);
#endif
        //std::cout << r_ <<u_<<x0_<<y0_<<lambda_<<gamma_<<j_<<std::endl;
    }
    DG_DEVICE
    double operator() (double x, double y)const
    {
        double radius = sqrt( (x-x0_)*(x-x0_) + (y-y0_)*(y-y0_));
        double theta = atan2( (y-y0_),(x-x0_));
 
        if( radius <= R_)
#ifdef _MSC_VER
            return 2.*lambda_*U_*_j1(lambda_*radius)/j_*cos( theta);
#else
            return 2.*lambda_*U_*j1( lambda_*radius)/j_*cos( theta);
#endif
        return 0;
    }
    double enstrophy( ) const { return M_PI*U_*U_*gamma_*gamma_;}
 
    double energy() const { return 2.*M_PI*R_*R_*U_*U_;}
  private:
    double R_, U_, x0_, y0_, lambda_, gamma_, j_;
};
 
struct Vortex
{
    Vortex( double x0, double y0, unsigned state,
          double R,  double u_dipole, double kz = 0):
        x0_(x0), y0_(y0), s_(state),  R_(R), u_d( u_dipole), kz_(kz){
        g_[0] = 3.831896621;
        g_[1] = -3.832353624;
        g_[2] = 7.016;
        b_[0] = 0.03827327723;
        b_[1] = 0.07071067810 ;
        b_[2] = 0.07071067810 ;
    }
    DG_DEVICE
    double operator()( double x, double y)const
    {
        double r = sqrt( (x-x0_)*(x-x0_)+(y-y0_)*(y-y0_));
        double theta = atan2( y-y0_, x-x0_);
        double beta = b_[s_];
        double norm = 1.2965125;
 
        if( r/R_<=1.)
            return u_d*(
                      r *( 1 +beta*beta/g_[s_]/g_[s_] )
#ifdef _MSC_VER
                    - R_*  beta*beta/g_[s_]/g_[s_] *_j1(g_[s_]*r/R_)/_j1(g_[s_])
#else
                    - R_ * beta*beta/g_[s_]/g_[s_] * j1(g_[s_]*r/R_)/ j1(g_[s_])
#endif
                    )*cos(theta)/norm;
        return u_d * R_* bessk1(beta*r/R_)/bessk1(beta)*cos(theta)/norm;
        // TODO Can these be replaced by std::cyl_bessel_k? Not sure what to do on device though
    }
    DG_DEVICE
    double operator()( double x, double y, double z)const
    {
        return this->operator()(x,y)*cos(kz_*z);
    }
    private:
    // Returns the modified Bessel function K1(x) for positive real x.
    DG_DEVICE
    double bessk1(double x)const
    {
        double y,ans;
        if (x <= 2.0)
        {
            y=x*x/4.0;
            ans = (log(x/2.0)*bessi1(x))+(1.0/x)*(1.0+y*(0.15443144 +
                       y*(-0.67278579+y*(-0.18156897+y*(-0.1919402e-1 +
                       y*(-0.110404e-2+y*(-0.4686e-4)))))));
        }
        else
        {
            y=2.0/x;
            ans = (exp(-x)/sqrt(x))*(1.25331414+y*(0.23498619 +
                      y*(-0.3655620e-1+y*(0.1504268e-1+y*(-0.780353e-2 +
                      y*(0.325614e-2+y*(-0.68245e-3)))))));
        }
        return ans;
    }
    //Returns the modified Bessel function I1(x) for any real x.
    DG_DEVICE
    double bessi1(double x) const
    {
        double ax,ans;
        double y;
        if ((ax=fabs(x)) < 3.75)
        {
            y=x/3.75;
            y*=y;
            ans = ax*(0.5+y*(0.87890594+y*(0.51498869+y*(0.15084934 +
                       y*(0.2658733e-1+y*(0.301532e-2+y*0.32411e-3))))));
        }
        else
        {
            y=3.75/ax;
            ans = 0.2282967e-1+y*(-0.2895312e-1+y*(0.1787654e-1 -
                      y*0.420059e-2)); ans=0.39894228+y*(-0.3988024e-1+
                      y*(-0.362018e-2 +y*(0.163801e-2+y*(-0.1031555e-1+y*ans))));
            ans *= (exp(ax)/sqrt(ax));
        }
        return x < 0.0 ? -ans : ans;
    }
    double x0_, y0_;
    unsigned s_;
    double R_, b_[3], u_d;
    double g_[3];
    double kz_;
};
 
struct BathRZ{
    BathRZ( unsigned N_kR, unsigned N_kZ, double R_min, double Z_min, double gamma, double L_E, double amp) :
        N_kR_(N_kR), N_kZ_(N_kZ),
        R_min_(R_min), Z_min_(Z_min),
        gamma_(gamma), L_E_(L_E) , amp_(amp),
        kvec( N_kR_*N_kZ_, 0), sqEkvec(kvec), unif1(kvec), unif2(kvec),
        normal1(kvec), normal2(kvec), alpha(kvec), theta(kvec)
    {
        double N_kR2=(double)(N_kR_*N_kR_);
        double N_kZ2=(double)(N_kZ_*N_kZ_);
        double N_k= sqrt(N_kR2+N_kZ2);
 
        norm_=sqrt(2./(double)N_kR_/(double)N_kZ_);
        double tpi=2.*M_PI, tpi2=tpi*tpi;
        double k0= tpi*L_E_/N_k;
        double N_kRh = N_kR_/2.;
        double N_kZh = N_kZ_/2.;
 
        thrust::minstd_rand generator;
        thrust::normal_distribution<double> ndistribution( 0.0, 1.0); // ( mean, stddev)
        thrust::uniform_real_distribution<double> udistribution(0.0,tpi); //between [0 and 2pi)
        for (unsigned j=1;j<=N_kZ_;j++)
        {
            double kZ2=tpi2*(j-N_kZh)*(j-N_kZh)/(N_kZ2);
            for (unsigned i=1;i<=N_kR_;i++)
            {
                double kR2=tpi2*(i-N_kRh)*(i-N_kRh)/(N_kR2);
                int z=(j-1)*(N_kR_)+(i-1);
                kvec[z]= sqrt(kR2 + kZ2);  //radial k number
                sqEkvec[z]=pow(kvec[z]*4.*k0/(kvec[z]+k0)/(kvec[z]+k0),gamma_/2.); //Energie in k space with max at 1.
                unif1[z]=cos(udistribution(generator));
                unif2[z]=sin(udistribution(generator));
                normal1[z]=ndistribution(generator);
                normal2[z]=ndistribution(generator);
                alpha[z]=sqrt(normal1[z]*normal1[z]+normal2[z]*normal2[z]);
                theta[z]=atan2(normal2[z],normal1[z]);
            }
        }
 
    }
    double operator()(double R, double Z)const
    {
        double f, kappa, RR, ZZ;
        RR=R-R_min_;
        ZZ=Z-Z_min_;
        f=0.;
        for (unsigned j=0;j<N_kZ_;j++)
        {
            for (unsigned i=0;i<N_kR_;i++)
            {
                int z=j*N_kR_+i;
                kappa= RR*unif1[z]+ZZ*unif2[z];
                f+= sqEkvec[z]*alpha[z]*cos(kvec[z]*kappa+theta[z]);
            }
        }
        return amp_*norm_*f;
    }
    double operator()(double R, double Z, double)const {
        double f, kappa;
        double  RR, ZZ;
        RR=R-R_min_;
        ZZ=Z-Z_min_;
        f=0;
        for (unsigned j=0;j<N_kZ_;j++)
        {
            for (unsigned i=0;i<N_kR_;i++)
            {
                int z=(j)*(N_kR_)+(i);
                kappa= RR*unif1[z]+ZZ*unif2[z];
                f+= sqEkvec[z]*alpha[z]*cos(kvec[z]*kappa+theta[z]);
            }
        }
        return amp_*norm_*f;
    }
  private:
    unsigned N_kR_,N_kZ_;
    double R_min_, Z_min_;
    double gamma_, L_E_;
    double amp_;
    double norm_;
    std::vector<double> kvec;
    std::vector<double> sqEkvec;
    std::vector<double> unif1, unif2, normal1,normal2,alpha,theta;
};
 
namespace detail
{
inline double horner( const double * c, unsigned M, double x)
{
    double b = c[M-1];
    for( unsigned i=0; i<M-1; i++)
        b = c[M-2-i] + b*x;
    return b;
}
 
} // namespace detail
 
struct Horner1d
{
    Horner1d(): m_c( 1, 1), m_M(1), m_prev( {0,1}){}
 
    Horner1d( const std::vector<double>& c): m_c(c), m_M(c.size()), m_prev( {1e300,1e300}){
        if( 0 == m_M )
        {
            // Make safe for zero coefficients
            m_c.resize(1, 0.);
            m_M = 1;
        }
    }
    double operator()( double x) const
    {
        if( m_prev[0] == x)
            return m_prev[1];
        m_prev[1] = detail::horner( &m_c[0], m_M, x);
        m_prev[0] = x;
        return m_prev[1];
    }
    private:
    std::vector<double> m_c;
    unsigned m_M;
    mutable std::array<double,2> m_prev;
};
 
struct Horner2d
{
    Horner2d(): m_c( 1, 1), m_cx( 1,1), m_M(1), m_N(1), m_prev( {0,0,1}){}
 
    Horner2d( const std::vector<double>& c, unsigned M, unsigned N): m_c(c), m_cx( M), m_M(M), m_N(N), m_prev( {1e300,1e300,1e300}){
        if( 0 == M || 0 == N)
        {
            // Make safe for zero coefficients
            m_c.resize(1, 0.);
            m_M = m_N = 1;
            m_cx.resize(1, 0.);
        }
    }
    double operator()( double x, double y) const
    {
        // Optimization rationale: Horner is typically repeatedly evaluated at
        // the same point (as it typically happens through the
        // dg::geo::polynomial::Psip classes in connection with magnetic field
        // quantities) or with constant y and varying x (such as when
        // dg::evaluate is called on Horner2d). In those cases significant
        // computation can be saved by storing the result of the previous call.
        if( m_prev[1] == y && m_prev[0] == x)
            return m_prev[2];
        if( m_prev[1] != y )
            for( unsigned i=0; i<m_M; i++)
                m_cx[i] = detail::horner( &m_c[i*m_N], m_N, y);
        m_prev[2] = detail::horner( &m_cx[0], m_M, x);
        m_prev[0] = x, m_prev[1] = y;
        return m_prev[2];
    }
    private:
    std::vector<double> m_c;
    mutable std::vector<double> m_cx;
    unsigned m_M, m_N;
    mutable std::array<double,3> m_prev;
};
 
template <class container = thrust::host_vector<double> >
struct Histogram
{
    Histogram(const dg::Grid1d& g1d, const std::vector<double>& in) :
    g1d_(g1d),
    in_(in),
    binwidth_(g1d_.h()),
    count_(dg::evaluate(dg::zero,g1d_))
    {
        for (unsigned j=0;j<in_.size();j++)
        {
            unsigned bin =floor( (in_[j]-g1d_.x0())/binwidth_ );
            bin = std::max(bin,(unsigned) 0);
            bin = std::min(bin,(unsigned)(g1d_.size()-1));
            count_[bin ]+=1.;
        }
        //Normalize
        unsigned Ampmax = (unsigned)thrust::reduce( count_.begin(), count_.end(),0.,   thrust::maximum<double>()  );
        dg::blas1::scal(count_,1./Ampmax);
 
    }
 
    double binwidth() {return binwidth_;}
    double operator()(double x)const
    {
        unsigned bin = floor((x-g1d_.x0())/binwidth_+0.5);
        bin = std::max(bin,(unsigned) 0);
        bin = std::min(bin,(unsigned)(g1d_.size()-1));
        return count_[bin];
    }
 
    private:
    dg::Grid1d g1d_;
    const std::vector<double> in_;
    double binwidth_;
    container  count_;
};
 
template <class container = thrust::host_vector<double> >
struct Histogram2D
{
    Histogram2D(const dg::Grid2d& g2d, const std::vector<double>& inx,const std::vector<double>& iny) :
    g2d_(g2d),
    inx_(inx),
    iny_(iny),
    binwidthx_(g2d_.hx()),
    binwidthy_(g2d_.hy()),
    count_(dg::evaluate(dg::zero,g2d_))
    {
 
        for (unsigned j=0;j<iny_.size();j++)
        {
            unsigned biny =floor((iny_[j]-g2d_.y0())/binwidthy_) ;
            biny = std::max(biny,(unsigned) 0);
            biny = std::min(biny,(unsigned)(g2d_.Ny()-1));
 
            unsigned binx =floor((inx_[j]-g2d_.x0())/binwidthx_) ;
            binx = std::max(binx,(unsigned) 0);
            binx = std::min(binx,(unsigned)(g2d_.Nx()-1));
            count_[biny*g2d_.Nx()+binx ]+=1.;
 
        }
        //Normalize
        unsigned Ampmax =  (unsigned)thrust::reduce( count_.begin(),   count_.end(),0.,thrust::maximum<double>()  );
        dg::blas1::scal(count_,  1./Ampmax);
 
    }
 
    double operator()(double x, double y)const
    {
        unsigned binx = floor((x-g2d_.x0())/binwidthx_+0.5) ;
        binx = std::max(binx,(unsigned) 0);
        binx = std::min(binx,(unsigned)(g2d_.Nx()-1));
        unsigned biny = floor((y-g2d_.y0())/binwidthy_+0.5) ;
        biny = std::max(biny,(unsigned) 0);
        biny = std::min(biny,(unsigned)(g2d_.Ny()-1));
        return count_[biny*g2d_.Nx()+binx ];
 
    }
    private:
    dg::Grid2d g2d_;
    const std::vector<double> inx_,iny_;
    double binwidthx_,binwidthy_;
    container count_;
};
 
 
 
struct WallDistance
{
    WallDistance( std::vector<double> vertical, std::vector<double> horizontal) :
        m_vertical(vertical), m_horizontal( horizontal) {}
    WallDistance( dg::Grid2d walls) : m_vertical({walls.x0(), walls.x1()}),
        m_horizontal({walls.y0(), walls.y1()}){}
    double operator() (double R, double Z) const
    {
        std::vector<double> dist( 1, 1e100); //fill in at least one (large) number in case vectors are empty)
        for( auto v : m_vertical)
            dist.push_back(fabs( R-v));
        for( auto h : m_horizontal)
            dist.push_back(fabs( Z-h));
        return *std::min_element( dist.begin(), dist.end());
    }
    private:
    std::vector<double> m_vertical;
    std::vector<double> m_horizontal;
};
 
 
} //namespace dg