fourier_filter.hpp

Go to the documentation of this file.

#ifndef FOURIER_FILTER_HPP
#define FOURIER_FILTER_HPP
 
#include "grid.hpp"
#include "domain_decomposition.hpp"
#include "plasma.hpp"
#include "complex.hpp"
#include "FFT_wrapper.hpp"
#include "spline_1d_interpolation.hpp"
 
/* Evaluate a damping factor for the boundary condition of the Fourier filter functions
 * @param[in]      psi          Input poloidal flux
 * @param[in]      inpsi        Inner filter boundary
 * @param[in]      outpsi       Outer filter boundary 
 * @param[in]      bd_width     Boundary width
 */
KOKKOS_INLINE_FUNCTION double get_damp_fac(double psi, double inpsi, double outpsi, double bd_width){
    double damp_fac = 0.5*(  tanh(2.0*(psi-inpsi )/bd_width)
                           - tanh(2.0*(psi-outpsi)/bd_width) );
    return (damp_fac>1.0e-6 ? damp_fac : 0.0);
}
 
class FourierFilter{
    public:
    enum class ScrapeOffLayerMethod{
        BlendWindow=0,
        RetainSineModes,
        RetainAllModes
    };
 
    enum class FilterMode{
        None=0,
        SingleN,
        M_Range,
        N_M_Range,
        N_Range_M_Resonant
    };
 
    struct Options{
        double inpsi;
        double outpsi;
        double bd_width;
        int minm;
        int maxm;
        bool is_spectral;
        bool is_resonant;
        int ntor_real;
        ScrapeOffLayerMethod sol_method;
 
        // For n-m filter
        bool div_mix;
        bool bands_on;
        int nmin;
        int nmax;
        int nmin_real;
        int nmax_real;
 
        FilterMode mode;
 
        // For single n
        int single_n;
    };
 
    private:
    int nsurf; //
    View<bool*,CLayout,HostType> run_filter;
    int max_n_points;
    View<bool*,CLayout,HostType> is_periodic;
    View<int*,CLayout,DeviceType> nnodes_on_surf;
    View<int*,CLayout,DeviceType> nfreq_on_surf;
    View<int*,CLayout,HostType> nuni_on_surf;
    View<int*,CLayout,DeviceType> node_offsets;
    View<int*,CLayout,HostType> node_offsets_h;
    View<int*,CLayout,DeviceType> freq_offsets;
    View<int*,CLayout,HostType> uni_offsets;
    View<int*,CLayout,DeviceType> node_surf_map;
    View<int*,CLayout,DeviceType> freq_surf_map;
    //View<int*,CLayout,HostType> uni_surf_map;
 
    View<double*,CLayout,DeviceType> inv_npt;
 
    public: // Because helper function needs them
    View<int*,CLayout,DeviceType> m_limit;
    View<double*,CLayout,DeviceType> qsafety;
    View<int*, CLayout, DeviceType> m_max_surf;
    int mres_q;
 
    private:
 
    View<int**,CLayout,DeviceType> m_minmax;
 
    // FFT
    std::vector<FFTWrapper> fft_plan;
 
    // Mapping
    View<int*,CLayout,DeviceType> mapping_all;
    View<double*,CLayout,DeviceType> theta_all;
    View<double*,CLayout,HostType> theta_all_h;
 
    View<UniformRange*,CLayout,HostType> theta_uniform_range_all;
 
    // Window
    View<double*,CLayout,DeviceType> window_all;
 
    // Spline precomputation
    std::vector<SplinePrecomputations<HostType>> spline_precomputation;
 
    // Resonance filter
    int b_sign;
    int wedge_n;
    int nplanes;
    int get_dir(int ntor_real) const;
 
    View<double*,CLayout,DeviceType> psi_surf;
 
    // Scratch
    View<double**, CLayout, DeviceType> x_ordered_all;
    View<double**, CLayout, HostType> x_ordered_h_view;
    std::vector<std::vector<View<double*, CLayout, HostType>>> x_ordered_h;
    View<double**, CLayout, HostType> x_uniform_h_view;
    std::vector<std::vector<View<double*, CLayout, HostType>>> x_uniform_h;
    View<Complex**,CLayout,DeviceType> spec_all;
    View<Complex**,CLayout,HostType> spec_h_view;
    std::vector<std::vector<View<Complex*, CLayout, HostType>>> spec_h;
 
    // Non-aligned vertices
#ifdef USE_GPU
    View<int*, CLayout, DeviceType> non_aligned_vert;
    View<int**, CLayout, DeviceType> non_aligned_ff_idx;
    View<double**, CLayout, DeviceType> non_aligned_ff_w;
#endif
    View<int*, CLayout, HostType> grid_non_aligned_vert;
    View<int**, CLayout, HostType> grid_non_aligned_ff_idx;
    View<double**, CLayout, HostType> grid_non_aligned_ff_w;
 
    void fft_execute(int isurf, const View<double*,CLayout,HostType>& x, const View<Complex*,CLayout,HostType>& spec) const;
 
    void set_nonaligned(const View<double**,CLayout,HostType>& x) const;
#ifdef USE_GPU
    void set_nonaligned(const View<double**,CLayout,DeviceType>& x) const;
#endif
 
    void m_range_one_surf(int isurf, const View<double**,CLayout,HostType>& x,const Options& options) const;
    void m_range_all_surfs(const View<double**,CLayout,DeviceType>& x,const Options& options) const;
 
    int get_n_points(int n_nodes) const;
    int get_n_frequencies(int n_points) const;
 
    public:
 
    FourierFilter(){}
 
    FourierFilter(NLReader::NamelistReader& nlr, const MagneticField<DeviceType>& magnetic_field, const Grid<DeviceType>& grid, Plasma& plasma);
 
 
#ifdef USE_GPU
    void m_range(const DomainDecomposition<DeviceType>& pol_decomp, const Options& options, const View<double*,CLayout,HostType>& x) const;
    void m_range(const DomainDecomposition<DeviceType>& pol_decomp, const Options& options, const View<double**,CLayout,HostType>& x) const;
#endif
    void m_range(const DomainDecomposition<DeviceType>& pol_decomp, const Options& options, const View<double*,CLayout,DeviceType>& x) const;
    void m_range(const DomainDecomposition<DeviceType>& pol_decomp, const Options& options, const View<double**,CLayout,DeviceType>& x) const;
 
    static void single_n(const Grid<DeviceType>& grid, const DomainDecomposition<DeviceType>& pol_decomp, const View<double*,CLayout,HostType>& view, const Options& options);
 
    static void single_n_flex(const Grid<DeviceType>& grid, const DomainDecomposition<DeviceType>& pol_decomp, const View<double*,CLayout,HostType>& view_in, const View<double**,CLayout,HostType>& result, const Options& options, int ntor_num);
 
    void n_m_range(int nplanes, const View<double*,CLayout,HostType>& view, const Options& options) const;
 
    void apply(const Grid<DeviceType>& grid, const DomainDecomposition<DeviceType>& pol_decomp, const Options& options, const View<double*,CLayout,HostType>& view) const;
 
    void destroy();
};
 
 
#endif