plane.hpp

Go to the documentation of this file.

#ifndef PLANE_HPP
#define PLANE_HPP

#include "xgc_io.hpp"
#include "magnetic_field.hpp"
#include "grid_structs.hpp"
#include "grid_weights.hpp"
#include "grid_files.hpp"
#include "guess_table.hpp"
#include "uniform_range.hpp"
#include "guess_list_1d.hpp"
#include "gradient_matrices.hpp"


// Plane class 
template<class Device>
class Plane {
    // 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 Plane;

    public:

    // Constructor from namelist
    Plane(NLReader::NamelistReader& nlr, const PlaneFiles& plane_files, const MagneticField<DeviceType>& magnetic_field, double plane_index_dbl, bool use_as_fortran_grid=false);

    // Default constructor
    Plane(){}

    // Create a mirror with a different device type
    template<class Device2>
    inline Plane<Device2> mirror() const{
        Plane<Device2> m;

        m.ntriangle = ntriangle;
        m.nnode = nnode;
        m.nplanes = nplanes;
        m.plane_phi = plane_phi;
        m.wedge_angle = wedge_angle;
        m.delta_phi = delta_phi;
        m.inv_delta_phi = inv_delta_phi;

        m.psi = my_mirror_view(psi, Device2());
        mirror_copy(m.psi, psi);
        m.bfield = my_mirror_view(bfield, Device2());
        mirror_copy(m.bfield, bfield);
        m.basis = my_mirror_view(basis, Device2());
        mirror_copy(m.basis, basis);
        m.nodes = my_mirror_view(nodes, Device2());
        mirror_copy(m.nodes, nodes);
        m.mapping = my_mirror_view(mapping, Device2());
        mirror_copy(m.mapping, mapping);

        m.guess = guess.template mirror<Device2>();

        m.psi_surf = my_mirror_view(psi_surf, Device2());
        mirror_copy(m.psi_surf, psi_surf);
        m.psi_surf_h = psi_surf_h; // Shallow copy since this stays on host
        m.npsi_surf2 = npsi_surf2;
        m.minval_psi_surf2 = minval_psi_surf2;
        m.maxval_psi_surf2 = maxval_psi_surf2;
        m.psi_surf2 = my_mirror_view(psi_surf2, Device2());
        mirror_copy(m.psi_surf2, psi_surf2);
        m.psi_map_surf_to_surf2 = my_mirror_view(psi_map_surf_to_surf2, Device2());
        mirror_copy(m.psi_map_surf_to_surf2, psi_map_surf_to_surf2);
        m.psi_map_surf2_to_surf = my_mirror_view(psi_map_surf2_to_surf, Device2());
        mirror_copy(m.psi_map_surf2_to_surf, psi_map_surf2_to_surf);

        m.psi00 = psi00;
        m.psi_guess = psi_guess.template mirror<Device2>();

        m.rgn = my_mirror_view(rgn, Device2());
        mirror_copy(m.rgn, rgn);
        m.gx = my_mirror_view(gx, Device2());
        mirror_copy(m.gx, gx);
        // This isn't used on GPU right now
        //m.adj = my_mirror_view(adj, Device2());
        //mirror_copy(m.adj, adj);

        m.nwall = nwall;
        m.wall_nodes = my_mirror_view(wall_nodes, Device2());
        mirror_copy(m.wall_nodes, wall_nodes);
        m.node_to_wall = my_mirror_view(node_to_wall, Device2());
        mirror_copy(m.node_to_wall, node_to_wall);

        m.eps_flux_surface = eps_flux_surface;

#ifdef NO_FORTRAN_MODULES
        m.gradient_matrices_h = gradient_matrices_h;
#endif

        return m;
    }

    void setup_1d_flux_surface_grid(const PlaneFiles& plane_files, const MagneticField<DeviceType>& magnetic_field, const View<double*,CLayout,HostType>& node_vol_h);

    KOKKOS_INLINE_FUNCTION void t_coeff(const SimdVector2D &x, const Simd<int>& itr, SimdGridVec &p) const;

    KOKKOS_INLINE_FUNCTION void t_coeff_mod(const MagneticField<Device> &magnetic_field, const double r, const double z, const double psiin, const int itr, double (&p)[3] ) const;

    KOKKOS_INLINE_FUNCTION void t_coeff_mod(const MagneticField<Device> &magnetic_field, const SimdVector2D &xy, const Simd<double>& psiin, const Simd<int>& itr, SimdGridVec &p ) const;

    KOKKOS_INLINE_FUNCTION void search_tr2_no_precheck(const double r, const double z, int& itr, double (&p)[3]) const;

    KOKKOS_INLINE_FUNCTION void search_tr2( const SimdVector2D &xy, Simd<int>& itr, SimdGridVec &pout ) const;

    KOKKOS_INLINE_FUNCTION void search_tr_check_guess(const SimdVector2D &x, const Simd<int>& old_itr,Simd<int>& itr,SimdGridVec &p ) const;

    KOKKOS_INLINE_FUNCTION void psi_search(double psi,double &wp,int &ip) const;

    KOKKOS_INLINE_FUNCTION void follow_field_to_nearest_midplane(const MagneticField<Device> &magnetic_field, const SimdVector2D &x, const Simd<double>& phi, SimdVector2D &xff) const;

    KOKKOS_INLINE_FUNCTION void get_grid_weights_ff(const MagneticField<Device> &magnetic_field, const SimdVector &v, const Simd<double>& psi_in,SimdVector2D &xff,SimdGridWeights<Order::One, PhiInterpType::Planes>& grid_wts) const;

    KOKKOS_INLINE_FUNCTION void get_grid_weights_plane_ff(const MagneticField<Device> &magnetic_field, const SimdVector &v, double phi_plane, SimdGridWeights<Order::One, PhiInterpType::Planes>& grid_wts) const;

    KOKKOS_INLINE_FUNCTION void get_grid_weights_no_ff(const MagneticField<Device> &magnetic_field, const SimdVector2D &x, const Simd<double>& psi_in,SimdGridWeights<Order::One, PhiInterpType::None>& grid_wts) const;

    template<PhiInterpType PIT>
    KOKKOS_INLINE_FUNCTION void nearest_node(const SimdGridWeights<Order::One, PIT>& grid_wts, SimdGridWeights<Order::Zero, PIT>& grid_wts0) const;

    template<PhiInterpType PIT>
    KOKKOS_INLINE_FUNCTION void get_grid_weights(const MagneticField<Device> &magnetic_field, const SimdVector &v, const Simd<double>& psi, SimdVector2D &xff,SimdGridWeights<Order::One, PIT>& grid_wts) const;

    template<PhiInterpType PIT>
    KOKKOS_INLINE_FUNCTION void get_grid_weights(const MagneticField<Device> &magnetic_field, const SimdVector &v, SimdVector2D &xff, SimdGridWeights<Order::One, PIT>& grid_wts) const;

    template<PhiInterpType PIT>
    KOKKOS_INLINE_FUNCTION void get_grid_weights(const MagneticField<Device> &magnetic_field, const SimdVector &v, const Simd<double>& psi_in, SimdGridWeights<Order::Zero, PIT>& grid_wts) const;

    template<PhiInterpType PIT>
    KOKKOS_INLINE_FUNCTION void get_grid_weights(const MagneticField<Device> &magnetic_field, const SimdVector &v, const Simd<double>& psi_in, SimdGridWeights<Order::One, PIT>& grid_wts) const;

    template<PhiInterpType PIT>
    KOKKOS_INLINE_FUNCTION void get_grid_weights(const MagneticField<Device> &magnetic_field, const SimdVector &v, SimdGridWeights<Order::Zero, PIT>& grid_wts) const;

    template<PhiInterpType PIT>
    KOKKOS_INLINE_FUNCTION void get_grid_weights(const MagneticField<Device> &magnetic_field, const SimdVector &v, SimdGridWeights<Order::One, PIT>& grid_wts) const;

    KOKKOS_INLINE_FUNCTION void get_wall_index(const Simd<bool>& just_left_the_grid,const SimdVector2D &x, const SimdGridWeights<Order::One, PIT_GLOBAL>& grid_wts, Simd<int>& widx) const;

    KOKKOS_INLINE_FUNCTION void wedge_modulo_phi(Simd<double>& phi_mod) const;

    KOKKOS_INLINE_FUNCTION double wedge_modulo_phi(double phi) const;

    KOKKOS_INLINE_FUNCTION bool node_is_inside_psi_range(const MagneticField<Device> &magnetic_field, const int node) const;

    KOKKOS_INLINE_FUNCTION int get_plane_index(double phi) const;

    KOKKOS_INLINE_FUNCTION int get_node_index(int triangle_index, int tri_vertex_index) const;

    KOKKOS_INLINE_FUNCTION double node_area_to_volume(const MagneticField<Device> &magnetic_field, double area, int node_index) const;

    KOKKOS_INLINE_FUNCTION void get_triangle_area_and_volume(const MagneticField<Device> &magnetic_field, int i, double& area, double& volume) const;

    KOKKOS_INLINE_FUNCTION double get_r_center_of_mass(int itr) const;

    KOKKOS_INLINE_FUNCTION bool node_is_in_region_1_or_2_no_wall(const int inode) const;

    KOKKOS_INLINE_FUNCTION bool node_is_in_region_2_or_3_no_wall(const int inode) const;

    KOKKOS_INLINE_FUNCTION bool node_is_in_private_region_no_wall(const int inode) const;

    KOKKOS_INLINE_FUNCTION bool node_is_in_included_region(const int inode, const bool exclude_private_region) const;

    KOKKOS_INLINE_FUNCTION bool node_is_in_region_1_or_2(const MagneticField<Device> &magnetic_field, const int inode) const;

    template<class Device2>
    KOKKOS_INLINE_FUNCTION bool node_is_in_region_3b(const MagneticField<Device2> &magnetic_field, const int inode) const;

    KOKKOS_INLINE_FUNCTION bool node_is_on_wall(const int inode) const;

    template<class Device2>
    KOKKOS_INLINE_FUNCTION bool node_is_in_region_3a_no_wall(const MagneticField<Device2> &magnetic_field, const int inode) const;

    template<class Device2>
    KOKKOS_INLINE_FUNCTION bool node_is_in_region_3b_no_wall(const MagneticField<Device2> &magnetic_field, const int inode) const;

    KOKKOS_INLINE_FUNCTION int uses_rz_basis(const int inode) const;

    KOKKOS_INLINE_FUNCTION double get_r(const int inode) const;

    KOKKOS_INLINE_FUNCTION double get_dist2_from_node(const int inode, const double r, const double z) const;

    KOKKOS_INLINE_FUNCTION double get_dist_from_node(const int inode, const double r, const double z) const;

    KOKKOS_INLINE_FUNCTION void get_rz_coordinates(const int inode, double& r, double& z) const;

    KOKKOS_INLINE_FUNCTION void get_rz_coordinates(const Simd<int>& grid_inds, SimdVector2D& x) const;

    KOKKOS_INLINE_FUNCTION void get_rz_coordinates(const int itr, const double (&p)[3], double& r, double& z) const;

    KOKKOS_INLINE_FUNCTION void get_rz_coordinates(const Simd<int>& itr, const SimdGridVec& p, SimdVector2D& x) const;

    KOKKOS_INLINE_FUNCTION int get_nearest_node(const Simd<int>& itr, const SimdGridVec& p, int i_simd) const;

    KOKKOS_INLINE_FUNCTION int get_nearest_node(const int itr, const double (&p)[3]) const;

    KOKKOS_INLINE_FUNCTION double get_nearest_midplane(double phi) const;

    KOKKOS_INLINE_FUNCTION void set_gradient_mat_triangle(const View<int*,CLayout,DeviceType>& num_t_node,
        const View<int**,CLayout,DeviceType>& tr_node, const View<double*,CLayout,DeviceType>& tr_area,
        const View<double***, CLayout, DeviceType>& unit_vecs,
        const int i, const GradientMatrices<Device>& gradient_matrices) const;

    KOKKOS_INLINE_FUNCTION int get_characteristic_length(const View<int*,CLayout,DeviceType>& num_t_node,
        const View<int**,CLayout,DeviceType>& tr_node, const int i,
        double& dist_triangle, double& dist_psi, double& dist_theta) const;

    KOKKOS_INLINE_FUNCTION void set_grad_matrix_from_psi_theta(bool is_psi_dir, int itr_pos, int itr_neg, double (&p_pos)[3], double (&p_neg)[3], double dl_pos, double dl_neg, int inode, const Matrix<Device>& matrix) const;

    KOKKOS_INLINE_FUNCTION RZPair get_wall_rz(int i_wall) const;

    void write_to_file(const XGC_IO_Stream& stream, const DomainDecomposition<DeviceType>& pol_decomp, bool is_stellarator, const std::string& prefix) const;
    
    void draw_ascii_grid(MagneticField<Device> magnetic_field, double rmin, double rmax, double dr, double zmin, double zmax, double dz){
        printf("\n\n******** PLOT OF GRID ******* \n");
        Kokkos::View<int*,HostType> save_vert("save_vert", 1+((rmax-rmin)/dr));
        for (double zplot = zmax; zplot>=zmin; zplot-=dz){
            printf("\n");
            int save_num = 0;
            int rsave_ind = -1;
            for (double rplot = rmin; rplot<=rmax; rplot+=dr){
                rsave_ind++;
                SimdVector v;
                v.r = rplot;
                v.z = zplot;
                v.phi = 0.001;

                // Get itr (triangle index)
                SimdGridVec p;
                Simd<int> itr = -1;
                get_grid_weights(magnetic_field, v, itr, p);
                if(itr[0]!=save_num || itr[0]!=save_vert[rsave_ind]){
                    if(itr[0]<10 && itr[0]>=0) printf("%d  ",itr[0]);
                    else printf("%d ",itr[0]);
                    save_num = itr[0];
                    save_vert[rsave_ind] = itr[0];
                } else {printf("   ");}
            }
        }
    }

    int ntriangle;         
    int nnode;             
    int nplanes;           
    double plane_phi;      
    double wedge_angle;    
    double delta_phi;      
    double inv_delta_phi;  

    Kokkos::View<double*,Kokkos::LayoutRight,Device> psi;   
    Kokkos::View<double*,Kokkos::LayoutRight,Device> bfield;   

    private:
    Kokkos::View<int*,   Kokkos::LayoutRight,Device> basis; 
    Kokkos::View<Vertex*,Kokkos::LayoutRight,Device> nodes;
    Kokkos::View<VertexMap*,Kokkos::LayoutRight,Device> mapping;
    Kokkos::View<int*,Kokkos::LayoutRight,Device> rgn;
    Kokkos::View<RZPair*,Kokkos::LayoutRight,Device> gx;
    GuessTable<Device> guess;

    public:
    int npsi_surf2;
    double minval_psi_surf2;
    double maxval_psi_surf2;
    Kokkos::View<double*,Kokkos::LayoutRight,Device> psi_surf;
    Kokkos::View<double*,Kokkos::LayoutRight,HostType> psi_surf_h;
    Kokkos::View<double*,Kokkos::LayoutRight,Device> psi_surf2;
    View<int*, CLayout, Device> psi_map_surf_to_surf2;
    View<int*, CLayout, Device> psi_map_surf2_to_surf;
    UniformRange psi00;

    private:
    GuessList1D<Device> psi_guess;

    View<int**,CLayout,Device> adj; 

    public:
    int nwall;

    //private:
    Kokkos::View<int*,Kokkos::LayoutRight,Device> wall_nodes;
    Kokkos::View<int*,Kokkos::LayoutRight,Device> node_to_wall;

    double eps_flux_surface; 

#ifdef NO_FORTRAN_MODULES
    GradientMatrices<HostType> gradient_matrices_h;
#endif
};

#include "plane.tpp"
#endif