grid_files.hpp

Go to the documentation of this file.

#ifndef GRID_FILES_HPP
#define GRID_FILES_HPP
 
#ifdef USE_MPI
#include "my_mpi.hpp"
#endif
#include "file_reader.hpp"
#include "magnetic_field.hpp"
#include "grid_structs.hpp"
#include "space_settings.hpp"
#include "NamelistReader.hpp"
 
class PlaneFilesInterface{
    public:
      virtual ~PlaneFilesInterface() = default;
    
      // Functions to query data from PlaneFiles
      [[nodiscard]] virtual int get_nnodes() const = 0;
      [[nodiscard]] virtual View<RZPair*, HostType> get_x() const = 0;
      [[nodiscard]] virtual View<int*, HostType> get_rgn() const = 0;
 
      [[nodiscard]] virtual int get_ntriangles() const = 0;
      [[nodiscard]] virtual View<Vertex*, HostType> get_nodes() const = 0;
 
      [[nodiscard]] virtual View<int*, HostType> get_i_x() const = 0;
 
      [[nodiscard]] virtual int get_nsurfaces() const = 0;
      [[nodiscard]] virtual int get_nsurfaces1() const = 0;
      [[nodiscard]] virtual int get_nsurfaces2() const = 0;
      [[nodiscard]] virtual int get_nsurfaces3a() const = 0;
      [[nodiscard]] virtual int get_nsurfaces3b() const = 0;
      [[nodiscard]] virtual View<int*, HostType> get_i_surf_separatrices() const = 0;
      [[nodiscard]] virtual View<int*, HostType> get_nnodes_on_surface() const = 0;
      [[nodiscard]] virtual int get_sum_nnodes_on_surfaces() const = 0;
      [[nodiscard]] virtual int get_num_non_aligned ()const = 0;
      [[nodiscard]] virtual View<int*, HostType> get_non_aligned_vert() const = 0;
      [[nodiscard]] virtual View<int*, HostType> get_non_aligned_nsurf() const = 0;
      [[nodiscard]] virtual View<int**, CLayout, HostType> get_non_aligned_surf_idx() const = 0;
      [[nodiscard]] virtual View<int**, CLayout, HostType> get_surf_idx() const = 0;
 
      [[nodiscard]] virtual View<int*, HostType> pack_scalars(int max_nodes_on_surface, int n_xpts, int n_separatrices) const = 0;
      virtual void unpack_scalars(const View<int*, HostType>& scalars, int& max_nodes_on_surface, int& n_xpts, int& n_separatrices) = 0;
 
      [[nodiscard]] virtual int n_specified_xpts() const = 0;
      [[nodiscard]] virtual View<Equil::XPoint*, HostType> get_xpts() const = 0;
      [[nodiscard]] virtual RZPair get_axis() const = 0;
      [[nodiscard]] virtual RZPair get_reference_point() const = 0;
      [[nodiscard]] virtual View<double*, HostType> get_psiN() const = 0; 
};
 
class PlaneFilesBase: public PlaneFilesInterface{
    public: 
      View<int*, HostType> pack_scalars(int max_nodes_on_surface, int n_xpts, int n_separatrices) const override;
      void unpack_scalars(const View<int*, HostType>& scalars, int& max_nodes_on_surface, int& n_xpts, int& n_separatrices) override;
   
    protected:
      // These are the regions designated in the .node file
        enum FileRegion{
          Inside=0,
          OnWall
      };
 
      // Contents of .node
      int nnodes;
      View<RZPair*, HostType> x;
      View<int*, HostType> rgn;
 
      // Contents of .ele
      int ntriangles;
      View<Vertex*, HostType> nodes;
 
      // Contents of .flx.aif
      View<int*, HostType> i_x;
 
      int nsurfaces;
      int nsurfaces1, nsurfaces2, nsurfaces3a, nsurfaces3b;
      View<int*, HostType> i_surf_separatrices;
      View<int*, HostType> nnodes_on_surface;
      int sum_nnodes_on_surfaces;
      int num_non_aligned;
      View<int*, HostType> non_aligned_vert;
      View<int*, HostType> non_aligned_nsurf;
      View<int**, CLayout, HostType> non_aligned_surf_idx;
      View<int**, CLayout, HostType> surf_idx;
 
      View<double*, HostType> psiN; 
};
 
struct PlaneFiles: public PlaneFilesBase{
    
    void convert_to_x_rgn_format(const View<double**, CLayout, HostType>& tmp_nodes, View<RZPair*, HostType>& x, View<int*, HostType>& rgn, bool is_stellarator);
 
    void convert_to_node_format(const View<int**, CLayout, HostType>& tmp_ele, View<Vertex*, HostType>& nodes);
 
    // Default constructor, to use to deallocate the file arrays after use
    PlaneFiles(){}
 
// Should remove this USE_MPI at some point
#ifdef USE_MPI
 
    // Copy data from a PlaneFiles object on a neighboring rank
    PlaneFiles(const MyMPI& my_mpi, const PlaneFilesInterface& plane_files);
 
    // Read the files
    PlaneFiles(const std::string& node_file, const std::string& element_file, const std::string& flx_aif_file, bool is_stellarator, const MPI_Comm& comm);
#endif
 
    // Options: hex vs circular
    enum GridCreateOpt{
        Hexagonal=0,
        Circular=1
    };
 
    void construct_hex_grid(int nshells, double raxis, double zaxis, double rscale, double zscale);
 
    double node_to_node_dist(RZPair a, RZPair b);
 
    void construct_circular_grid(int nshells, double raxis, double zaxis, double rscale, double zscale);
 
    static int nshells_from_nnodes(int nnodes);
 
    PlaneFiles(int nshells, double raxis, double zaxis, double rscale, double zscale, GridCreateOpt grid_opt=Hexagonal);
 
    PlaneFiles(MagneticField<HostType>& magnetic_field, int nshells, GridCreateOpt grid_opt=Hexagonal);
 
    // Functions
    int get_nnodes() const override;
    View<RZPair*, HostType> get_x() const override;
    View<int*, HostType> get_rgn() const override;
 
    int get_ntriangles() const override;
    View<Vertex*, HostType> get_nodes() const override;
 
    View<int*, HostType> get_i_x() const override; 
 
    int get_nsurfaces() const override;
    int get_nsurfaces1() const override;
    int get_nsurfaces2() const override;
    int get_nsurfaces3a() const override;
    int get_nsurfaces3b() const override;
    View<int*, HostType> get_i_surf_separatrices() const override;
    View<int*, HostType> get_nnodes_on_surface() const override; 
    int get_sum_nnodes_on_surfaces() const override;
    int get_num_non_aligned() const override;
    View<int*, HostType> get_non_aligned_vert() const override;
    View<int*, HostType> get_non_aligned_nsurf() const override;
    View<int**, CLayout, HostType> get_non_aligned_surf_idx() const override;
    View<int**, CLayout, HostType> get_surf_idx() const override;
 
    int n_specified_xpts() const override;
    View<Equil::XPoint*, HostType> get_xpts() const override;
    RZPair get_axis() const override;
    RZPair get_reference_point() const override;
    View<double*, HostType> get_psiN() const override;
};
 
struct GridFiles{
    using GridCreateOpt = PlaneFiles::GridCreateOpt;
    static constexpr GridCreateOpt Hexagonal = PlaneFiles::Hexagonal;
    static constexpr GridCreateOpt Circular = PlaneFiles::Circular;
 
    std::vector<std::unique_ptr<PlaneFilesInterface>> plane_files_vec;
 
    // For test grids
    GridFiles(int nshells, double raxis, double zaxis, double rscale, double zscale, GridCreateOpt grid_opt=Hexagonal);
 
    GridFiles(MagneticField<HostType>& magnetic_field, int nshells, GridCreateOpt grid_opt=Hexagonal);
 
    std::string ifile_int2str(int ifile) const;
 
#ifdef USE_MPI
    GridFiles(NLReader::NamelistReader& nlr, bool is_stellarator, const MyMPI& my_mpi);
#endif
 
    const PlaneFilesInterface& lplane() const;
 
    const PlaneFilesInterface& midplane() const;
 
    const PlaneFilesInterface& rplane() const;
 
    int n_planes() const;
 
    void clear();
 
    static int nshells_from_nnodes(int nnodes);
};
 
#endif