profile.hpp

Go to the documentation of this file.

#ifndef PROFILE_HPP
#define PROFILE_HPP

#include "magnetic_field.hpp"

// Putting this in a namespace since Shape names sound so generic
namespace Eq {

enum Shape{
    Constant=0,
    TanH,
    Linear,
    Shape3,
    Shape4,
    Shape5,
    Shape6,
    CustomSpline,
    CustomLinear
};

//parameters for linear equal distance interpolation of profiles
//maybe need some compiler flags for generality. Currently hard coded to 1000 and true.
const int ftn_lin_n=1000; // 1000 data point for linear interpolation

// data structure for linear interpolation -- used in eq_ftn_type
template<class Device>
class CustomLinShape{
    int n;
    double p_min,p_del;
    Kokkos::View<double*,Kokkos::LayoutRight,Device> v;

    public:

    CustomLinShape(double p_min_in, double p_max_in, double* v_in)
        : n(ftn_lin_n),
          p_min(p_min_in),
          p_del((p_max_in - p_min_in)/(ftn_lin_n-1)),
          v("v",ftn_lin_n)
    {
        array_deep_copy(v,v_in);
    }

    CustomLinShape(){}

    KOKKOS_INLINE_FUNCTION double interp(double x) const{
        double xn=(x-p_min)/p_del;
        xn=max(min(xn,double(n)),0.0);
        int ip=int(xn);
        ip=min(ip,n-2);
        double wp=1.0-(xn-double(ip));
        return v(ip)*wp+v(ip+1)*(1.0-wp);
    }

    KOKKOS_INLINE_FUNCTION double dinterp(double x) const{
        double xn=(x-p_min)/p_del;
        xn=max(min(xn,double(n)),0.0);
        int ip=int(xn);
        ip=min(ip,n-2);
        return (v(ip+1)-v(ip))/p_del;
    }
};

// data structure for equilbirum profile
template<class Device>
class Profile{
    Shape shape;
    double inx[5];
    double iny[4];
    double sv[6];
    double p_min;
    double p_max;

    // for linear interpolation
    CustomLinShape<Device> lin;

    // get value of profile at the given coordinates
    KOKKOS_INLINE_FUNCTION double value_at_coords(const MagneticField<DeviceType> &b_field, double psi,double r,double z) const {
        if(shape==Constant){
            return iny[0];
        } else if(shape==TanH){
            return sv[1]*tanh((inx[0]-psi)*sv[2]) + sv[0];
        } else if(shape==Linear){
            return sv[0]*psi+sv[1];
        } else if(shape==Shape5){
            if(abs(psi/b_field.equil.xpt_psi)<1.0e-3) psi = 1.0e-3*b_field.equil.xpt_psi;
            double tmp2=sv[2]*(inx[0]-sqrt(inx[0]*psi)); // z
            double tmp3=exp(tmp2); // expz
            double tmp4=exp(-tmp2); // expmz
            // A * ( (1+z*slope)*expz - expmz )/(expz + expmz) + B
            return sv[1]*((1+tmp2*sv[3])*tmp3-tmp4)/(tmp3+tmp4) + sv[0];
        } else if(shape==CustomLinear){
            psi=min(max(psi,p_min),p_max);
            return lin.interp(psi);
        } else{ // Shouldnt get here
            return 0.0;
        }
    }

    public:

    // Profile from Fortran
    Profile(double* inx_in, double* iny_in, double* sv_in, double p_min_in, double p_max_in, int shape_int, double* v_in )
        : p_min(p_min_in),
          p_max(p_max_in)
    {
        for (int i=0; i<5; i++) inx[i]=inx_in[i];
        for (int i=0; i<4; i++) iny[i]=iny_in[i];
        for (int i=0; i<6; i++) sv[i]=sv_in[i];

        // Currently, only a few shapes available
        if(shape_int==0){
            shape=Constant;
        } else if(shape_int==1){
            shape=TanH;
        } else if(shape_int==2){
            shape=Linear;
        } else if(shape_int==5){
            shape=Shape5;
        } else if(shape_int==-1){ // shape==-1 is also CustomSpline; this is caught in the init right now
            shape=CustomLinear;
            lin = CustomLinShape<Device>(p_min, p_max, v_in);
        } else {
            printf("invalid shape number in eq_ftn_setup: %d",shape_int);
            exit(0);
        }
    }

    // Constructor for simple linear profile
    Profile(double val_at_zero_psi, double slope)
     : shape(Linear)
    {
        for (int i=0; i<5; i++) inx[i]=0.0;
        for (int i=0; i<4; i++) iny[i]=0.0;
        for (int i=0; i<6; i++) sv[i]=0.0;
        sv[0] = slope;
        sv[1] = val_at_zero_psi;
    }

    // Profile for analytic testing
    Profile(Shape shape_in)
     : shape(shape_in),
       p_min(0.0),
       p_max(0.0)
    {
        for (int i=0; i<5; i++) inx[i]=0.0;
        for (int i=0; i<4; i++) iny[i]=0.0;
        for (int i=0; i<6; i++) sv[i]=0.0;
        if(shape==Constant){
            iny[0] = 2.0; // value
        } else if(shape==TanH){
            inx[0] = 1.0;
            sv[0] = 3.0; // half height
            sv[1] = 3.0; // value at psi=inx[0]
            sv[2] = 5.0; // width
        } else if(shape==Linear){
            sv[0] = 2.0; // slope
            sv[1] = 3.0; // value at psi=0
        } else if(shape==Shape5){
            inx[0] = 2.5e-1;
            sv[0] = 5.5e2;
            sv[1] = 4.5e2;
            sv[2] = 2.1;
            sv[3] = 8.0e-3;
        } else if(shape==CustomLinear){
            double v_in[ftn_lin_n];
            p_min = 0.0;
            p_max = 1.0;
            for (int i = 0; i<ftn_lin_n; i++){
                v_in[i] = 1.0+i/(ftn_lin_n-1.0);
            }
            lin = CustomLinShape<Device>(p_min, p_max, v_in);
        }
    }

    Profile(){}

    // get value of profile, including decay if out of range
    KOKKOS_INLINE_FUNCTION double value(const MagneticField<DeviceType> &b_field, double psi_in,double r,double z) const {
        bool is_in_rgn12 = b_field.equil.is_in_region_1_or_2(r,z,psi_in);
        if(is_in_rgn12 && (psi_in<=b_field.outpsi)){
            return value_at_coords(b_field,psi_in,r,z);
        } else if(is_in_rgn12 && (psi_in>b_field.outpsi)){
            // Get function value at psi=b_field.outpsi
            // let ftn decay exponentially to eq_out_decay_factor
            // times the value at b_field.outpsi with a decay
            // length of eq_out_decay_width
            double psi=b_field.outpsi;
            double tmp2=value_at_coords(b_field,psi,b_field.equil.axis_r,b_field.equil.axis_z);
            double tmp3=tmp2*(1.0-b_field.equil.out_decay_factor);
            return tmp3*exp(-abs(psi_in-psi)/b_field.equil.out_decay_width)+b_field.equil.out_decay_factor*tmp2;
        } else if(!is_in_rgn12 && z<b_field.equil.xpt_z){
            // For primary private flux region
            double psi=b_field.equil.xpt_psi;
            double tmp2=value_at_coords(b_field,psi,b_field.equil.axis_r,b_field.equil.axis_z);
            double tmp3=tmp2*(1.0-b_field.equil.priv_flux_decay_factor);
            return tmp3*exp(-abs(psi_in-psi)/b_field.equil.priv_flux_decay_width)+b_field.equil.priv_flux_decay_factor*tmp2;
        } else if(!is_in_rgn12 && z>b_field.equil.xpt2_z && b_field.equil.set_xpt2){
            // For primary private flux region
            double psi=b_field.equil.xpt2_psi;
            double tmp2=value_at_coords(b_field,psi,b_field.equil.axis_r,b_field.equil.axis_z);
            double tmp3=tmp2*(1.0-b_field.equil.priv_flux_decay_factor);
            return tmp3*exp(-abs(psi_in-psi)/b_field.equil.priv_flux_decay_width)+b_field.equil.priv_flux_decay_factor*tmp2;
        } else {
            // No way to get here? - ALS
            double psi=b_field.equil.xpt_psi;
            return value_at_coords(b_field,psi,r,z);
        }
    }

    // get derivative of profile
    KOKKOS_INLINE_FUNCTION double slope(const MagneticField<DeviceType> &b_field, double psi_in,double r,double z) const {
        // region searching and enforcing region 3 value to x-point value
        //rh return derivative=0 for psi values larger than sml_outpsi
        //rh This is not a serious constraint since eq_dftn is called only in
        //rh the particle push to get the turbulence drive for delta-f simulations
        if(!b_field.equil.is_in_region_1_or_2(r,z,psi_in) || (psi_in>b_field.outpsi)){
           return 0.0;
        }

        if(shape==Constant){
           return 0.0;
        } else if(shape==TanH){
           double tmp2 = 1.0/cosh((inx[0]-psi_in)*sv[2]);
           return -sv[2]*sv[1]*tmp2*tmp2;
        } else if(shape==Linear){
           return sv[0];
        } else if(shape==Shape5){
           double tmp2=sv[2]*(inx[0]-sqrt(inx[0]*psi_in)); // z
           double tmp3=exp(tmp2); // expz
           // A * ( (1+z*slope)*expz - expmz )/(expz + expmz) + B
           if(abs(psi_in/b_field.equil.xpt_psi)<1.0e-3){
             return 0.0;
           } else {
             return -inx[0]*sv[1]*sv[2]*tmp3*tmp3*(4.0+sv[3]*(1.0+tmp3*tmp3+2.0*tmp2))
                        /(2.0*(1.0+tmp3*tmp3)*(1.0+tmp3*tmp3)*(sqrt(inx[0]*psi_in)));
           }
        } else if(shape==CustomLinear){
            if(p_min<psi_in && psi_in<p_max){
                return lin.dinterp(psi_in);
            } else {
                return 0.0;
            }
        } else{ // Shouldnt get here
            return 0.0;
        }
    }

};
}
#endif