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"
#include "volumes_and_areas.hpp"
#include "flux_surface_average.hpp"
#include "magnetic_equil_files.hpp"

// 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;
            // Except plane_phi
            m.rplane.plane_phi = rplane.plane_phi;
            m.lplane.plane_phi = lplane.plane_phi;
        }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
        m.fsa_matrices = fsa_matrices.template mirror<Device2>();
        if(flux_surfaces_h.size()>0){
            m.flux_surfaces_h = HostArray<VertexList>(flux_surfaces_h.size());
            for(int isurf=0; isurf<flux_surfaces_h.size(); isurf++) m.flux_surfaces_h(isurf) = flux_surfaces_h(isurf);
        }

        return m;
    }

    void psi_surf_region_check(const MagneticField<DeviceType>& magnetic_field, const PlaneFiles& right_plane_files);

    void write_to_file(const DomainDecomposition<DeviceType>& pol_decomp, MagneticEquilFiles::Ptr equil_files_ptr, bool is_stellarator) 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() 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_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;

    template<class Device2>
    KOKKOS_INLINE_FUNCTION bool node_is_in_region_1_or_2(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 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;

    KOKKOS_INLINE_FUNCTION int get_ip(const int inode) 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
    FluxSurfaceAverage<Device> fsa_matrices;
    HostArray<VertexList> flux_surfaces_h;
};

#include "grid.tpp"
#endif