flux_surface_average.hpp

Go to the documentation of this file.

#ifndef FLUX_SURFACE_AVERAGE_HPP
#define FLUX_SURFACE_AVERAGE_HPP

#include "plane.hpp"
#include "volumes_and_areas.hpp"

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

    Matrix<Device> to_flux_surface;
    Matrix<Device> from_flux_surface;
    View<double*, CLayout, Device> surf_vol;
    View<double*, CLayout, Device> surf_vol_3D;

    void wall_nodes_strikep_method(NLReader::NamelistReader& nlr, const MagneticField<DeviceType>& magnetic_field, const Plane<HostType>& plane, const PlaneFiles& plane_files, const View<double*,CLayout,HostType>& psi_surf, const View<int*,CLayout,DeviceType>& num_t_node, const View<int**,CLayout,DeviceType>& tr_node);

    public:

    FluxSurfaceAverage(){}

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

        m.to_flux_surface = to_flux_surface.template mirror<Device2>();
        m.from_flux_surface = from_flux_surface.template mirror<Device2>();

        m.surf_vol = my_mirror_view(surf_vol, Device2());
        mirror_copy(m.surf_vol, surf_vol);

        m.surf_vol_3D = my_mirror_view(surf_vol_3D, Device2());
        mirror_copy(m.surf_vol_3D, surf_vol_3D);

        return m;
    }

    FluxSurfaceAverage(NLReader::NamelistReader& nlr, const DomainDecomposition<DeviceType>& pol_decomp, const MagneticField<DeviceType>& magnetic_field, const Plane<HostType>& plane, const PlaneFiles& plane_files, const VolumesAndAreas& volumes_and_areas, const View<int*,CLayout,DeviceType>& num_t_node, const View<int**,CLayout,DeviceType>& tr_node);

    View<double*, CLayout, Device> to_surf(const View<double*, CLayout, Device>& input) const;

    void from_surf(const View<double*, CLayout, Device>& input, const View<double*, CLayout, Device>& output) const;

    View<double*, CLayout, Device> average(const View<double*, CLayout, Device>& array) const;

    View<double*, CLayout, Device> average(const View<double*, CLayout, Device>& array, const View<double*, CLayout, Device>& output_surf) const;

    void average(const View<double*, CLayout, Device>& array, const View<double*, CLayout, Device>& output_surf, const View<double*, CLayout, Device>& tmp00) const;

    void average(const View<double**, CLayout, Device>& array) const;

    View<double*, CLayout, Device> average_3D(const DomainDecomposition<DeviceType>& pol_decomp, const View<double*, CLayout, Device>& array) const;
    View<double*, CLayout, Device> average_3D(const DomainDecomposition<DeviceType>& pol_decomp, const View<double*, CLayout, Device>& array, const View<double*, CLayout, Device>& output_array) const;

    int nsurfs() const;

    // Return index of psi surface on which the input vertex falls. Returns -1 for non-aligned vertices
    KOKKOS_INLINE_FUNCTION int get_node_surf_index(int inode) const{
        if(to_flux_surface.nelement(inode)==1){ // Non-aligned vertices do not have a mapping to surface so nelement is 0 for them
            return to_flux_surface.eindex(inode,0) - 1; // Map to psi_surf
        }else{
            return -1;
        }
    }
};

#endif