tricub.hpp

Go to the documentation of this file.

#ifndef TRICUB_HPP
#define TRICUB_HPP
#include "space_settings.hpp"
#include "simd.hpp"
#include "my_mirror_view.hpp"

struct TricubCoeff {
    double coeff[4][4][4];
};

// Tricubic spline class
template<class Device>
class Tricub {
    // This line enables private access between the host and device implementations of this class
    // i.e. Makes the copy from host to device more straightforward
    template<class Device2> friend class Tricub;

    // enums to read indices more easily
    enum{ J0 = 0, J1 = 1, J2 = 2, J3 = 3 };
    enum{ K0 = 0, K1 = 1, K2 = 2, K3 = 3 };

    // bicub variables
    int nr, nz, np;            
    double rmin, zmin, pmin;     
    double dr_inv, dz_inv, dp_inv; 
    Kokkos::View<double*,Kokkos::LayoutRight,Device> rc; 
    Kokkos::View<double*,Kokkos::LayoutRight,Device> zc; 
    Kokkos::View<double*,Kokkos::LayoutRight,Device> pc; 
    Kokkos::View<TricubCoeff***,Kokkos::LayoutRight,Device> acoeff_all;  

    KOKKOS_INLINE_FUNCTION void eval_tricub_0(double x, double y, double z, double xc, double yc, double zc, const double (&acoeff)[4][4][4],
     double &f000) const;
    static KOKKOS_INLINE_FUNCTION double interp_1d(const double (&coeff1d)[4], double dx);
    static KOKKOS_INLINE_FUNCTION double slope_1d(const double (&coeff1d)[4], double dx);

    static KOKKOS_INLINE_FUNCTION double interp_2d(const double (&coeff2d)[4][4], double dx, double dy);

    KOKKOS_INLINE_FUNCTION void eval_tricub_1(double x, double y, double z, double xc, double yc, double zc, const double (&acoeff)[4][4][4],
     double &f000, double &f100, double &f010, double &f001) const;

    public:

    // Constructor from grid on r,z
    Tricub(const View<double*, HostType>& r, const View<double*, HostType>& z, const View<double*, HostType>& phi, const View<double***, CLayout, HostType>& val, bool keep_positive=false);

    // Create a mirror with a different device type
    template<class Device2>
    inline Tricub<Device2> mirror() const{
        Tricub<Device2> m;
        m.nr = nr;
        m.nz = nz;
        m.np = np;
        m.rmin = rmin;
        m.dr_inv = dr_inv;
        m.zmin = zmin;
        m.dz_inv = dz_inv;
        m.pmin = pmin;
        m.dp_inv = dp_inv;
        m.rc = my_mirror_view(rc, Device2());
        mirror_copy(m.rc, rc);
        m.zc = my_mirror_view(zc, Device2());
        mirror_copy(m.zc, zc);
        m.pc = my_mirror_view(pc, Device2());
        mirror_copy(m.pc, pc);
        m.acoeff_all = my_mirror_view(acoeff_all, Device2());
        mirror_copy(m.acoeff_all, acoeff_all);

        return m;
    }

    // Default constructor
    Tricub(){}

    // Get values
    KOKKOS_INLINE_FUNCTION void der_zero(double x, double y, double z, double &f000) const;

    // Get first derivatives
    KOKKOS_INLINE_FUNCTION void der_one(double x, double y, double z, double &f000, double &f100, double &f010, double &f001) const;
};

#include "tricub.tpp"

#endif