grid.hpp

Go to the documentation of this file.

#ifndef GRID_HPP
#define GRID_HPP

#include "magnetic_field.hpp"
#include "bounds.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 "plane.hpp"
#include "ff_projection.hpp"
#include "domain_decomposition.hpp"

extern "C" void init_fortran_vols_and_areas(double* tr_area_h, double* tr_vol_h,
#ifdef XGC1
                    double* node_vol_ff_h,
#endif
                    double* node_vol_nearest_h, double* node_vol_h, double* node_area_h, double* psi00_h, int npsi00);

// Volume- and Area-related Host Views. Should find a different location for these so that these unused views aren't sent to GPU
struct VolumesAndAreas{
    View<double*,CLayout,HostType> tr_area_h;
    View<double*,CLayout,HostType> tr_vol_h;
    View<double**,CLayout,HostType> node_vol_ff_h; // Default to midplane, for now
    View<double**,CLayout,HostType> node_vol_ff_l_h; // Field following volume for left plane nodes
    View<double**,CLayout,HostType> node_vol_ff_r_h; // Field following volume for right plane nnodes
    View<double*,CLayout,HostType> node_vol_nearest_h;
    View<double*,CLayout,HostType> node_vol_h;
    View<double*,CLayout,HostType> node_area_h;
    View<double*,CLayout,HostType> vol00_h;

    VolumesAndAreas(){}

    VolumesAndAreas(int ntriangle_in, int nnode_in, int npsi00_in)
      : tr_area_h(NoInit("tr_area"), ntriangle_in),
        tr_vol_h(NoInit("tr_vol"), ntriangle_in),
#ifdef XGC1
        node_vol_ff_h(NoInit("node_vol_ff"),2, nnode_in),
#ifdef STELLARATOR
        node_vol_ff_l_h(NoInit("node_vol_ff_l"),2, nnode_in),
        node_vol_ff_r_h(NoInit("node_vol_ff_r"),2, nnode_in),
#else
        node_vol_ff_l_h(node_vol_ff_h), // Shallow copy of midplane node_vol_ff
        node_vol_ff_r_h(node_vol_ff_h),
#endif
#endif
        node_vol_nearest_h(NoInit("node_vol_nearest"),nnode_in),
        node_vol_h(NoInit("node_vol"),nnode_in),
        node_area_h(NoInit("node_area"),nnode_in),
        vol00_h(NoInit("vol00_h"), npsi00_in)
    {}

    void set_fortran_ptrs(){
#ifndef NO_FORTRAN_MODULES
        int npsi00 = vol00_h.size();
        init_fortran_vols_and_areas(tr_area_h.data(), tr_vol_h.data(),
#ifdef XGC1
                node_vol_ff_h.data(),
#endif
                node_vol_nearest_h.data(), node_vol_h.data(), node_area_h.data(),
                vol00_h.data(), npsi00);
#endif
    }
};

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

    public:

    // Constructor from namelist
    Grid(NLReader::NamelistReader& nlr, const GridFiles& grid_files, const DomainDecomposition<DeviceType>& pol_decomp, const MagneticField<DeviceType>& magnetic_field);

    // Default constructor
    Grid(){}

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

        m.axisymmetric = axisymmetric;

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

        m.npsi_surf2 = npsi_surf2;
        m.minval_psi_surf2 = minval_psi_surf2;
        m.maxval_psi_surf2 = maxval_psi_surf2;
        m.psi00 = psi00;

        if(axisymmetric){
            m.midplane = midplane.template mirror<Device2>();
            m.rplane = m.midplane; // Shallow copy of midplane if axisymmetric
            m.lplane = m.midplane;
        }else{
            m.midplane = midplane.template mirror<Device2>();
            m.rplane = rplane.template mirror<Device2>();
            m.lplane = lplane.template mirror<Device2>();
        }

        // Set these views to be their midplane equivalents (shallow copy)
        m.psi = m.midplane.psi;
        m.bfield = m.midplane.bfield;
        m.psi_surf2 = m.midplane.psi_surf2;
        m.wall_nodes = m.midplane.wall_nodes;

        m.volumes_and_areas = volumes_and_areas; // Doesn't require mirror view since everything is on host

        m.half_plane_ff = half_plane_ff;
        m.ff_rplane_to_neighbors = ff_rplane_to_neighbors;
        m.ff_lplane_to_neighbors = ff_lplane_to_neighbors;
        m.ff_to_midplane = ff_to_midplane;

#ifdef NO_FORTRAN_MODULES
        m.gradient_matrices_h = gradient_matrices_h;
#endif

        return m;
    }

    void write_to_file() const;

    KOKKOS_INLINE_FUNCTION void psi_search(double psi,double &wp,int &ip) 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;

    Bounds get_local_wedge_bounds(int plane_index) 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 bool node_is_in_region_1_or_2_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;

    KOKKOS_INLINE_FUNCTION int uses_rz_basis(const int inode) const;

    KOKKOS_INLINE_FUNCTION double get_r(const int inode) 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 RZPair get_wall_rz(int i_wall) const;

    bool axisymmetric;

    // Planes
    Plane<Device> midplane;
    Plane<Device> lplane;
    Plane<Device> rplane;

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

    int npsi_surf2;
    double minval_psi_surf2;
    double maxval_psi_surf2;
    UniformRange psi00;
    Kokkos::View<double*,Kokkos::LayoutRight,Device> psi;   
    Kokkos::View<double*,Kokkos::LayoutRight,Device> bfield;   
    Kokkos::View<double*,Kokkos::LayoutRight,Device> psi_surf2;
    Kokkos::View<int*,Kokkos::LayoutRight,Device> wall_nodes;

    // Volume- and Area-related Host Views. Should find a different location for these so that these unused views aren't sent to GPU
    VolumesAndAreas volumes_and_areas;

    // Projections - reside in host memory for now
    Projection<HostType> half_plane_ff;
    Projection<HostType> ff_rplane_to_neighbors;
    Projection<HostType> ff_lplane_to_neighbors;
    Projection<HostType> ff_to_midplane;


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

#include "grid.tpp"
#endif