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"

// Temporary views, preallocated for reuse by each OMP thread
struct FilterTmpViews{
    View<double*,CLayout,HostType> tmp;
    View<int*,CLayout,HostType> node_mapping;
    View<double*,CLayout,HostType> tmp2;
    View<double*,CLayout,HostType> theta;
    View<double*,CLayout,HostType> win;
    View<Complex**,CLayout,HostType> spec;
                
    FilterTmpViews(){}
            
    FilterTmpViews(int n_dim2, int n_max)
      : tmp(NoInit("tmp"), n_max),
        node_mapping(NoInit("node_mapping"), n_max),
        tmp2(NoInit("tmp2"), n_max),
        theta(NoInit("theta"), n_max),
        win(NoInit("win"), n_max),
        spec(NoInit("spec"), n_dim2, n_max){}

    // Create a view that is the size needed for a specific surface
    static const View<double*,CLayout,HostType> view(const View<double*,CLayout,HostType>& in, int n){
        return View<double*,CLayout,HostType, Kokkos::MemoryTraits<Kokkos::Unmanaged>>(in.data(), n);
    }
    
    static const View<int*,CLayout,HostType> view(const View<int*,CLayout,HostType>& in, int n){
        return View<int*,CLayout,HostType, Kokkos::MemoryTraits<Kokkos::Unmanaged>>(in.data(), n);
    }

    static const View<Complex**,CLayout,HostType> view(const View<Complex**,CLayout,HostType>& in, int n_dim2, int n){
        return View<Complex**,CLayout,HostType, Kokkos::MemoryTraits<Kokkos::Unmanaged>>(in.data(), n_dim2, n);
    }
};

class FourierFilter{
    public:
    enum class ScrapeOffLayerMethod{
        BlendWindow=0,
        RetainSineModes,
        RetainAllModes
    };

    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;
    };

    private:
    View<bool*,CLayout,HostType> run_filter;
    int mres_q;
    int max_n_points;
    View<int*,CLayout,HostType> m_limit;
    View<bool*,CLayout,HostType> is_periodic;
    View<int*,CLayout,HostType> nnodes_on_surf;

    // FFT
    std::vector<FFTWrapper> fft_plan;

    // Mapping
    std::vector<View<int*,CLayout,HostType>> mapping;
    std::vector<View<double*,CLayout,HostType>> theta;

    // Window
    std::vector<View<double*,CLayout,HostType>> window;

    // Spline precomputation
    std::vector<SplinePrecomputations<HostType>> spline_precomputation;

    // Resonance filter
    int b_sign;

    // Scratch
    std::vector<View<double**, CLayout, HostType>> x_ordered_all;
    std::vector<View<Complex**,CLayout,HostType>> spec_all;
    std::vector<View<double**, CLayout, HostType>> x_uniform_all;

    void get_m_min_max(const Options& options, int isurf, double qsafety, int& minm, int& maxm) const;

    void fft_execute(int isurf, const View<double*,CLayout,HostType>& x, const View<Complex*,CLayout,HostType>& spec) const;

    void m_range_one_surf(const Grid<DeviceType>& grid, int isurf, const View<double**,CLayout,HostType>& 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);

    void m_range(const Grid<DeviceType>& grid, const DomainDecomposition<DeviceType>& pol_decomp, const Options& options, const View<double*,CLayout,HostType>& x) const;

    void m_range(const Grid<DeviceType>& grid, const DomainDecomposition<DeviceType>& pol_decomp, const Options& options, const View<double**,CLayout,HostType>& x) const;

    void destroy();
};


void fourier_filter(const Grid<DeviceType>& grid, const DomainDecomposition<DeviceType>& pol_decomp, const View<double*,CLayout,HostType>& view);

#endif