profile.hpp

Go to the documentation of this file.

#ifndef PROFILE_HPP
#define PROFILE_HPP

#include "file_reader.hpp"
#include "magnetic_field.hpp"
#ifdef USE_MPI
#include "my_mpi.hpp"
#endif

extern "C" void init_ftn_lin(int num, double* psi, double* var, double psi_min, double dpsi, int n, double* v);

// 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
};

// data structure for linear interpolation
template<class Device>
class CustomLinShape{
    static constexpr int n = 1000; // 1000 data point for linear interpolation
    double p_min;
    double p_max;
    double p_del;
    Kokkos::View<double*,Kokkos::LayoutRight,Device> v;

    public:

    CustomLinShape<Device>(double psi_normalization, const std::string& filename)
       : v(NoInit("v"),n)
    {
        // Only one MPI rank reads the files, then broadcasts to the others
        int num;
        View<double*, CLayout, HostType> psi;
        View<double*, CLayout, HostType> var;
        if(is_rank_zero()){
            // Open file
            FileReader file_reader;
            file_reader.open(filename);
            file_reader.read(&num, 1);
            if(num<=0) exit_XGC("\nError in profile ftn init : invalid number of data\n");

            psi = View<double*, CLayout, HostType>(NoInit("psi"), num);
            var = View<double*, CLayout, HostType>(NoInit("var"), num);

            // Read data
            for(int i=0; i<num; i++){
                file_reader.read(&psi(i), 1);
                file_reader.read(&var(i), 1);
            }

            // Check end of file
            int flag;
            file_reader.read(&flag, 1);
            if(flag!=-1) exit_XGC("\nError in profile ftn init : invalid number of data or ending flag -1 is not set correctly\n");
        }

        // Broadcast file size to other processes
#ifdef USE_MPI
        int root_rank = 0;
        MPI_Bcast(&num, 1, MPI_INT, root_rank, SML_COMM_WORLD);
#endif
        if(!is_rank_zero()){
            psi = View<double*, CLayout, HostType>(NoInit("psi"), num);
            var = View<double*, CLayout, HostType>(NoInit("var"), num);
        }
        // Broadcast arrays
#ifdef USE_MPI
        MPI_Bcast( psi.data(), psi.size(), MPI_DOUBLE, root_rank, SML_COMM_WORLD);
        MPI_Bcast( var.data(), var.size(), MPI_DOUBLE, root_rank, SML_COMM_WORLD);
#endif

        // Normalization
        for(int i=0; i<num; i++){
            psi(i) *= psi_normalization;
        }

        // save psi range
        p_min = psi(0);
        p_max = psi(num-1);
        p_del = (p_max - p_min)/(n-1);

        // Set value of v on host
        View<double*,CLayout,HostType> v_h(NoInit("v_h"), v.layout());
#ifdef PSPLINE
        init_ftn_lin(num, psi.data(), var.data(), p_min, p_del, n, v_h.data());
#else
        exit_XGC("\nError: CustomLinear profile shape requires PSPLINE.\n");
#endif

        // Copy to device
        Kokkos::deep_copy(v, v_h);
    }

    // Constructor provided an array of values and psis
    CustomLinShape<Device>(const View<double*, CLayout, HostType>& psi, const View<double*, CLayout, HostType>& var)
      : v(NoInit("v"),n)
    {
        int num = var.size();

        // save psi range
        p_min = psi(0);
        p_max = psi(num-1);
        p_del = (p_max - p_min)/(n-1);

        // Set value of v on host
        View<double*,CLayout,HostType> v_h(NoInit("v_h"), v.layout());
#ifdef PSPLINE
        init_ftn_lin(num, psi.data(), var.data(), p_min, p_del, n, v_h.data());
#else
        exit_XGC("\nError: CustomLinear profile shape requires PSPLINE.\n");
#endif

        // Copy to device
        Kokkos::deep_copy(v, v_h);
    }

    // Constructor with linear value for testing
    CustomLinShape<Device>(double p_min_in, double p_max_in)
      : v(NoInit("v"),n),
        p_min(p_min_in),
        p_max(p_max_in),
        p_del((p_max_in-p_min_in)/(n-1))
    {
        View<double*,CLayout,HostType> v_h(NoInit("v_h"), v.layout());
        for (int i = 0; i<n; i++){
            v_h(i) = 1.0 + i/(n-1.0);
        }

        // Copy to device
        Kokkos::deep_copy(v, v_h);
    }

    CustomLinShape(){}

    KOKKOS_INLINE_FUNCTION double interp(double x) const{
        // Bound by psi limits - this seems redundant with the bounds immediately following
        x = min(max(x,p_min),p_max);
        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{
        if(p_min<x && x<p_max){
            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;
        } else {
            return 0.0;
        }
    }
};

// Data structure for equilibrium profile
template<class Device>
class Profile{
    public:
    static constexpr int N_X = 5; // SK: N_X and N_Y can be reduced to 3
    static constexpr int N_Y = 4;
    static constexpr int N_V = 6;

    private:

    Shape shape;
    double inx[N_X];
    double iny[N_Y];
    double sv[N_V];

    // for linear interpolation
    CustomLinShape<Device> lin;

    // get value of profile at the given coordinates
    KOKKOS_INLINE_FUNCTION double value_at_coords(double psi_norm, double psi) 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)<1.0e-3*psi_norm) psi = 1.0e-3*psi_norm;
            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){
            return lin.interp(psi);
        } else{ // Shouldnt get here
            return 0.0;
        }
    }

    public:

    // This function converts from the integer currently provided in the namelist
    // into the enum
    static Shape shape_from_int(int shape_int){
        if(shape_int==0){
            return Constant;
        } else if(shape_int==1){
            return TanH;
        } else if(shape_int==2){
            return Linear;
        } else if(shape_int==5){
            return Shape5;
        } else if(shape_int==-1){ // shape==-1 is also CustomSpline; this is caught in the init right now
            return CustomLinear;
        } else {
            printf("Invalid shape number in profile setup: %d",shape_int);
            exit_XGC("\nError: Invalid shape number\n");
            return Constant; // Not reached
        }
    }

    Profile(double psi_normalization, Shape shape_in, const std::vector<double>& inx_in, const std::vector<double>& iny_in, const std::string& filename)
      : shape(shape_in)
    {
        // Copy vector inputs into member C-style arrays
        for(int i=0; i<inx_in.size(); i++) inx[i] = inx_in[i];
        for(int i=0; i<iny_in.size(); i++) iny[i] = iny_in[i];
        

        // Precalculations depending on profile shape
        if(shape==Constant){
            // No additional setup needed
        }else if(shape==TanH){
            sv[0]=0.5*(iny[0]+iny[1]);
            sv[1]=0.5*(iny[0]-iny[1]);
            sv[2]=2/inx[1];
        }else if(shape==Linear){
            sv[0]=(iny[1]-iny[0])/ inx[1]; // slope
            sv[1]=iny[0]-sv[0]*inx[0];     // offset
        }else if(shape==Shape5){
            // A ( (1+z slope) e^z - e^-z )/(e^z + e^-z) + B
            // z= 4 ( center - sqrt(center psi) )/ width
            // iny(1) - edge_val : ped_top
            // iny(2) - out_val : ped bottom
            // inx(1) - center ,  inx(2) - width
            // iny(3) - core(axis)_val , inx(3) - core(axis)_psi
            sv[0]=0.5*(iny[0]+iny[1]);
            sv[1]=0.5*(iny[0]-iny[1]);    // height
            sv[2]=4/inx[1];                   // width inverse
            sv[3]=(iny[2]-iny[0]) / (sqrt(inx[2]) - sqrt(inx[0]))
                  /( - (sv[1]+1.0e-99) * sv[2] * sqrt(inx[0]) );  // slope
        }else if(shape==CustomLinear){
            lin = CustomLinShape<Device>(psi_normalization, filename);
        }else{
            exit_XGC("\nError: Invalid profile shape provided.\n");
        }
    }

    // Constructor that takes a pre-existing custom profile x and var
    Profile(const View<double*, CLayout, HostType>& psi, const View<double*, CLayout, HostType>& var)
      : shape(CustomLinear),
        lin(CustomLinShape<Device>(psi, var))
    {}

    // 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)
    {
        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){
            // Linear value for testing
            double p_min = 0.0;
            double p_max = 1.0;
            lin = CustomLinShape<Device>(p_min, p_max);
        }
    }

    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.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_norm(),psi_in);
        } 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_norm(),psi);
            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_norm(),psi);
            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_norm(),psi);
            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_norm(),psi);
        }
    }

    // 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.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){
            return lin.dinterp(psi_in);
        } else{ // Shouldnt get here
            return 0.0;
        }
    }

    void get_shape_inx_iny(int& shape_int, std::vector<double>& inx_out, std::vector<double>& iny_out) {

        shape_int = get_shape_int();

        for(int i=0; i<N_X; i++) inx_out[i] = inx[i];
        for(int i=0; i<N_Y; i++) iny_out[i] = iny[i];
    }



    // This function returns shape int from shape as the integer provied in the namelist
    int get_shape_int() {
        if(shape==Constant){
            return 0;
       }else if(shape==TanH){
            return 1;
        }else if(shape==Linear){
            return 2;
        }else if(shape==Shape5){
            return 5;
        }else if(shape==CustomLinear){
            return -1;
        }else{
            exit_XGC("\nError: Invalid profile shape provided.\n");
            return -1; // To turn off void return warning message
        }
    }

};
}
#endif