convected_helmholtz_3d_hf_fmm.hpp

Go to the documentation of this file.

 
#ifndef NIHU_CONVECTED_HELMHOLTZ_3D_HF_FMM_HPP_INCLUDED
#define NIHU_CONVECTED_HELMHOLTZ_3D_HF_FMM_HPP_INCLUDED
 
#include "cluster_tree.hpp"
#include "fmm_operator.hpp"
#include "helmholtz_3d_hf_cluster.h"
#include "helmholtz_3d_hf_shift.h"
#include "m2l_indices.hpp"
#include "operator_with_wave_number.hpp"
 
#include "library/convected_helmholtz_3d_kernel.hpp"
#include "library/convected_helmholtz_3d_singular_integrals.hpp"
#include "library/convected_helmholtz_3d_nearly_singular_integrals.hpp"
 
#include <boost/math/constants/constants.hpp>
#include <boost/math/special_functions/hankel.hpp>
#include <boost/math/special_functions/legendre.hpp>
 
#include <vector>
 
namespace NiHu
{
namespace fmm
{
 
template <class WaveNumber>
class convected_operator
    : public operator_with_wave_number<WaveNumber>
{
public:
    typedef Eigen::Matrix<double, 3, 3> dmatrix3_t;
    typedef WaveNumber wave_number_t;
 
    convected_operator(wave_number_t const &k, wave_number_t const &kappa, dmatrix3_t const &T)
        : operator_with_wave_number<wave_number_t>(k)
        , m_kappa(kappa)
        , m_T(T)   
    {
    }
 
    dmatrix3_t get_T() const
    {
        return m_T;
    }
 
    wave_number_t const &get_kappa() const
    {
        return m_kappa;
    }
 
private:
    wave_number_t const m_kappa;
    dmatrix3_t const m_T;
};
 
template <class WaveNumber>
class convected_helmholtz_3d_hf_fmm
{
public:
    using wave_number_t = WaveNumber;
    using cluster_t = helmholtz_3d_hf_cluster;
    using cluster_tree_t = cluster_tree<cluster_t>;
    using bounding_box_t = typename cluster_t::bounding_box_t;
    using location_t = typename bounding_box_t::location_t;
    using dmatrix3_t = Eigen::Matrix<double, 3, 3>;
    using cvector_t = Eigen::Matrix<std::complex<double>, Eigen::Dynamic, 1>;
 
    convected_helmholtz_3d_hf_fmm(wave_number_t const &wn, location_t const &mach_vector)
        : m_wave_number(wn)
        , m_mach_vector(mach_vector)
        , m_mach_number(mach_vector.norm())
        , m_kappa(m_wave_number / (1. - m_mach_number * m_mach_number))
        , m_accuracy(3.0)
    {
        // compute transform matrix
        location_t e = m_mach_vector.normalized();
        dmatrix3_t E = e * e.transpose();
        m_T = (dmatrix3_t::Identity() - E) * std::sqrt(1 - m_mach_number*m_mach_number) + E;
    }
 
    static size_t compute_expansion_length(double drel, double C)
    {
        using namespace boost::math::double_constants;
        auto kd = two_pi * drel;
        return size_t(std::ceil(kd + C * std::log(kd + pi)));
    }
 
    void set_accuracy(double C)
    {
        m_accuracy = C;
    }
 
    void init_level_data(cluster_tree_t const &tree)
    {
        using namespace boost::math::double_constants;
        double lambda_scaled = two_pi / std::real(m_kappa);
        m_level_data_vector.clear();
        m_level_data_vector.resize(tree.get_n_levels());
 
        // set the expansion length for each level
        for (size_t i = 0; i < tree.get_n_levels(); ++i)
        {
            auto &ld = m_level_data_vector[i];
            // get the diameter from the first cluster on the level
            auto idx = tree.level_begin(i);
            auto d = tree[idx].get_bounding_box().get_diameter();
            auto L = compute_expansion_length(d / lambda_scaled, m_accuracy);
            ld.set_expansion_length(L);
        }
 
        // compute interpolation matrices
        for (size_t i = 0; i < tree.get_n_levels(); ++i)
        {
            auto &ld = m_level_data_vector[i];
            auto const &Sto = ld.get_unit_sphere();
            if (i >= 3)
                ld.set_interp_dn(interpolator(m_level_data_vector[i - 1].get_unit_sphere(), Sto));
            if (i < tree.get_n_levels() - 1 && i >= 2)
                ld.set_interp_up(interpolator(m_level_data_vector[i + 1].get_unit_sphere(), Sto));
        }
    } // end of function init_level_data
 
    helmholtz_3d_hf_level_data const &get_level_data(size_t idx) const
    {
        return m_level_data_vector[idx];
    }
 
    helmholtz_3d_hf_level_data &get_level_data(size_t idx)
    {
        return m_level_data_vector[idx];
    }
 
    class m2m
        : public convected_operator<wave_number_t>
        , public fmm_operator<m2m_tag>
    {
    public:
        using wn_base_t = convected_operator<wave_number_t>;
        using result_t = helmholtz_3d_hf_fmm_upshift;
        using cluster_t = convected_helmholtz_3d_hf_fmm::cluster_t;
 
        m2m(wave_number_t const &wave_number, wave_number_t const &kappa, dmatrix3_t const &T)
            : wn_base_t(wave_number, kappa, T)
        {
        }
 
        static size_t unique_idx(cluster_t const &to, cluster_t const &from)
        {
            return bounding_box_t::dist2idx(
                from.get_bounding_box().get_center(),
                to.get_bounding_box().get_center());
        }
 
        result_t operator()(cluster_t const &to, cluster_t const &from) const
        {
            using namespace std::complex_literals;
            auto const &Sto = to.get_level_data().get_unit_sphere();
            location_t TD = this->get_T() * (to.get_bounding_box().get_center()
                - from.get_bounding_box().get_center());
            auto const &kappa = this->get_kappa();
            cvector_t shift = exp((-1.i * kappa * Sto.get_s().transpose() * TD).array());
            return result_t(shift, to.get_level_data());
        }
    };
 
 
    class l2l
        : public convected_operator<wave_number_t>
        , public fmm_operator<l2l_tag>
    {
    public:
        using wn_base_t = convected_operator<WaveNumber>;
        using result_t = helmholtz_3d_hf_fmm_downshift;
        using cluster_t = convected_helmholtz_3d_hf_fmm::cluster_t;
 
        l2l(wave_number_t const &wave_number, wave_number_t const &kappa, dmatrix3_t const &T)
            : wn_base_t(wave_number, kappa, T)
        {
        }
 
        static size_t unique_idx(cluster_t const &to, cluster_t const &from)
        {
            return bounding_box_t::dist2idx(
                to.get_bounding_box().get_center(),
                from.get_bounding_box().get_center());
        }
 
        result_t operator()(cluster_t const &to, cluster_t const &from) const
        {
            using namespace std::complex_literals;
#if L2L_SHIFT_FIRST
            auto const &Sfrom = from.get_level_data().get_unit_sphere();
#else
            auto const &Sto = to.get_level_data().get_unit_sphere();
#endif
            location_t TD = this->get_T() * (to.get_bounding_box().get_center()
                - from.get_bounding_box().get_center());
            auto const &kappa = this->get_kappa();
#if L2L_SHIFT_FIRST
            cvector_t shift = exp((-1.i * kappa * Sfrom.get_s().transpose() * TD).array());
#else
            cvector_t shift = exp((-1.i * kappa * Sto.get_s().transpose() * TD).array());
#endif
            return result_t(shift, to.get_level_data());
        }
    };
 
 
    class m2l
        : public convected_operator<wave_number_t>
        , public fmm_operator<m2l_tag>
    {
    public:
        using convected_base_t = convected_operator<wave_number_t>;
        using result_t = Eigen::DiagonalMatrix<std::complex<double>, Eigen::Dynamic>;
        using cluster_t = convected_helmholtz_3d_hf_fmm::cluster_t;
 
        m2l(wave_number_t const &wave_number, wave_number_t const &kappa, dmatrix3_t const &T)
            : convected_base_t(wave_number, kappa, T)
        {
        }
 
        static size_t unique_idx(cluster_t const &to, cluster_t const &from)
        {
            return m2l_indices<3>::eval(to.get_bounding_box(), from.get_bounding_box());
        }
 
        result_t operator()(cluster_t const &to, cluster_t const &from) const
        {
            auto const &kappa = this->get_kappa();
 
            auto const &X = to.get_bounding_box().get_center();
            auto const &Y = from.get_bounding_box().get_center();
            location_t TDvec = this->get_T() * (X - Y);
 
            auto L = to.get_level_data().get_expansion_length();
            auto const &s = to.get_level_data().get_unit_sphere().get_s();
 
            return m2l_matrix_impl(TDvec, s, kappa, L).asDiagonal();
        }
 
    private:
        static cvector_t m2l_matrix_impl(location_t const &TDvec,
            Eigen::Matrix<double, 3, Eigen::Dynamic> const &s,
            wave_number_t const &kappa,
            size_t L)
        {
            using namespace boost::math::double_constants;
            using namespace std::complex_literals;
            using boost::math::legendre_p;
            using boost::math::sph_hankel_1;
 
            double TD = TDvec.norm();
            location_t TUvec = TDvec / TD;
            auto N = s.cols();
 
            Eigen::Matrix<double, 1, Eigen::Dynamic> x = TUvec.transpose() * s;
            auto z = -kappa * TD;
 
            cvector_t M2L = cvector_t::Zero(N, 1);
 
            for (size_t l = 0; l <= L; ++l)
            {
                auto h = sph_hankel_1(l, z);
                auto c = double(2 * l + 1) * std::pow(1.i, l) * h;
                for (int i = 0; i < s.cols(); ++i)
                    M2L(i) += c * legendre_p(int(l), x(0, i));
            }
            M2L *= -1.i * kappa / (4.0 * pi) / (4.0 * pi);
 
            return M2L;
        }
    };
 
    class p2m_SLP
        : public convected_operator<wave_number_t>
        , public fmm_operator<p2m_tag>
    {
    public:
        using test_input_t = cluster_t;
        using trial_input_t = typename NiHu::convected_helmholtz_3d_SLP_kernel<wave_number_t>::trial_input_t;
        using result_t = Eigen::Matrix<std::complex<double>, Eigen::Dynamic, 1>;
 
        p2m_SLP(wave_number_t const &k, wave_number_t const &kappa, dmatrix3_t const &T, location_t const &mach_vector)
            : convected_operator<wave_number_t>(k, kappa, T)
            , m_mach_vector(mach_vector)
        {
        }
 
        size_t rows(test_input_t const &to) const
        {
            return to.get_level_data().get_unit_sphere().get_num_directions();
        }
 
        result_t operator()(test_input_t const &to, trial_input_t const &tri) const
        {
            return eval(to, tri);
        }
 
    private:
        result_t eval(test_input_t const &to, trial_input_t const &tri) const
        {
            using namespace std::complex_literals;
            auto const &Y = to.get_bounding_box().get_center();
            auto const &y = tri.get_x();
            auto const &s = to.get_level_data().get_unit_sphere().get_s();
            auto const &kappa = this->get_kappa();
            auto const &Mvec = this->m_mach_vector;
            location_t Td = this->get_T() * (Y - y);
            double M_dot_ny = Mvec.dot(tri.get_unit_normal());
            return Eigen::exp(-1.i * kappa * (s.transpose() * Td).array())
                * std::exp(-1.i * kappa * Mvec.dot(y))
                * (1. - M_dot_ny*M_dot_ny);
        }
 
        location_t m_mach_vector;
    };
 
    class p2m_DLP
        : public convected_operator<wave_number_t>
        , public fmm_operator<p2m_tag>
    {
    public:
        using test_input_t = cluster_t;
        using trial_input_t = typename NiHu::convected_helmholtz_3d_DLP_kernel<wave_number_t>::trial_input_t;
        using result_t = Eigen::Matrix<std::complex<double>, Eigen::Dynamic, 1>;
 
        p2m_DLP(wave_number_t const &k, wave_number_t const &kappa, dmatrix3_t const &T, location_t const &mach_vector)
            : convected_operator<wave_number_t>(k, kappa, T)
            , m_mach_vector(mach_vector)
        {
        }
 
        size_t rows(test_input_t const &to) const
        {
            return to.get_level_data().get_unit_sphere().get_num_directions();
        }
 
        result_t operator()(test_input_t const &to, trial_input_t const &tri) const
        {
            return eval(to, tri);
        }
 
    private:
        result_t eval(test_input_t const &to, trial_input_t const &tri) const
        {
            using namespace std::complex_literals;
            auto const &Y = to.get_bounding_box().get_center();
            auto const &y = tri.get_x();
            auto const &ny = tri.get_unit_normal();
            auto const &s = to.get_level_data().get_unit_sphere().get_s();
            auto const &kappa = this->get_kappa();
            auto const &Mvec = this->m_mach_vector;
            double M2 = Mvec.dot(Mvec);
            auto const &T = this->get_T();
            location_t Td = T * (Y - y);
            return Eigen::exp(-1.i * kappa * (s.transpose() * Td).array()) * std::exp(-1.i * kappa * Mvec.dot(y))
                * (1.i * kappa) *
                ((s.transpose() * (T * ny)).array() + Mvec.dot(ny) * (1. - M2 - (s.transpose()*Mvec).array()));
        }
 
        location_t m_mach_vector;
    };
 
    class p2l
        : public convected_operator<wave_number_t>
        , public fmm_operator<p2l_tag>
    {
    public:
        using test_input_t = cluster_t;
        using trial_input_t = typename NiHu::convected_helmholtz_3d_SLP_kernel<wave_number_t>::trial_input_t;
        using result_t = Eigen::Matrix<std::complex<double>, Eigen::Dynamic, 1>;
 
        p2l(wave_number_t const &k, wave_number_t const &kappa, dmatrix3_t const &T, location_t const &mach_vector)
            : convected_operator<wave_number_t>(k, kappa, T)
            , m_mach_vector(mach_vector)
        {
        }
 
        size_t rows(test_input_t const &to) const
        {
            return to.get_level_data().get_unit_sphere().get_num_directions();
        }
 
        result_t operator()(test_input_t const &to, trial_input_t const &tri) const
        {
            return eval(to, tri);
        }
 
    private:
        result_t eval(test_input_t const &to, trial_input_t const &tri) const
        {
            throw std::logic_error("Unimplemented p2l operator");
        }
 
        location_t m_mach_vector;
    };
 
    class m2p
        : public convected_operator<wave_number_t>
        , public fmm_operator<m2p_tag>
    {
    public:
        using test_input_t = typename NiHu::convected_helmholtz_3d_SLP_kernel<wave_number_t>::test_input_t;
        using trial_input_t = cluster_t;
        using result_t = Eigen::Matrix<std::complex<double>, Eigen::Dynamic, 1>;
 
        m2p(wave_number_t const &k, wave_number_t const &kappa, dmatrix3_t const &T, location_t const &mach_vector)
            : convected_operator<wave_number_t>(k, kappa, T)
            , m_mach_vector(mach_vector)
        {
        }
 
        size_t cols(trial_input_t const &from) const
        {
            return from.get_level_data().get_unit_sphere().get_num_directions();
        }
 
        result_t operator()(test_input_t const &tsi, trial_input_t const &from) const
        {
            return eval(tsi, from);
        }
 
    private:
        result_t eval(test_input_t const &tsi, trial_input_t const &from) const
        {
            throw std::logic_error("Unimplemented m2p operator");
        }
 
        location_t m_mach_vector;
    };
 
    class l2p
        : public convected_operator<wave_number_t>
        , public fmm_operator<l2p_tag>
    {
    public:
        using test_input_t = typename NiHu::convected_helmholtz_3d_SLP_kernel<wave_number_t>::test_input_t;
        using trial_input_t = cluster_t;
        using result_t = Eigen::Matrix<std::complex<double>, 1, Eigen::Dynamic>;
 
        l2p(wave_number_t const &k, wave_number_t const &kappa, dmatrix3_t const &T, location_t const &mach_vector)
            : convected_operator<wave_number_t>(k, kappa, T)
            , m_mach_vector(mach_vector)
        {
        }
 
        size_t cols(trial_input_t const &from) const
        {
            return from.get_level_data().get_unit_sphere().get_num_directions();
        }
 
        result_t operator()(test_input_t const &tsi, trial_input_t const &from) const
        {
            return eval(tsi, from);
        }
 
    private:
        result_t eval(test_input_t const &tsi, trial_input_t const &from) const
        {
            using namespace std::complex_literals;
            auto const &X = from.get_bounding_box().get_center();
            auto const &x = tsi.get_x();
            auto const &S = from.get_level_data().get_unit_sphere();
            auto const &s = S.get_s();
            auto const &w = S.get_w();
            auto const &kappa = this->get_kappa();
            auto const &Mvec = this->m_mach_vector;
            location_t Td = this->get_T() * (x - X);
            return Eigen::exp(-1.i * kappa * (s.transpose() * Td).array()) * w.array() * std::exp(1.i * kappa * Mvec.dot(x));
        }
 
        location_t m_mach_vector;
    };
 
    class p2p_SLP
        : public fmm_operator<p2p_tag>
    {
    public:
        typedef NiHu::convected_helmholtz_3d_SLP_kernel<wave_number_t> kernel_t;
        typedef typename kernel_t::result_t result_t;
        typedef typename kernel_t::test_input_t test_input_t;
        typedef typename kernel_t::trial_input_t trial_input_t;
 
        p2p_SLP(wave_number_t const &k, location_t const &mvec)
            : m_kernel(k, mvec)
        {
        }
 
        result_t operator()(test_input_t const &tsi, trial_input_t const &tri) const
        {
            return m_kernel(tsi, tri);
        }
 
        kernel_t const &get_kernel() const
        {
            return m_kernel;
        }
 
    private:
        kernel_t m_kernel;
    };
 
    class p2p_DLP
        : public fmm_operator<p2p_tag>
    {
    public:
        typedef NiHu::convected_helmholtz_3d_DLP_kernel<wave_number_t> kernel_t;
        typedef typename kernel_t::result_t result_t;
        typedef typename kernel_t::test_input_t test_input_t;
        typedef typename kernel_t::trial_input_t trial_input_t;
 
        p2p_DLP(wave_number_t const &k, location_t const &mvec)
            : m_kernel(k, mvec)
        {
        }
 
        result_t operator()(test_input_t const &tsi, trial_input_t const &tri) const
        {
            return m_kernel(tsi, tri);
        }
 
        kernel_t const &get_kernel() const
        {
            return m_kernel;
        }
 
    private:
        kernel_t m_kernel;
    };
 
    m2m create_m2m() const
    {
        return m2m(m_wave_number, m_kappa, m_T);
    }
 
    l2l create_l2l() const
    {
        return l2l(m_wave_number, m_kappa, m_T);
    }
 
    m2l create_m2l() const
    {
        return m2l(m_wave_number, m_kappa, m_T);
    }
 
    p2l create_p2l() const
    {
        return p2l(m_wave_number, m_kappa, m_T, m_mach_vector);
    }
 
    m2p create_m2p() const
    {
        return m2p(m_wave_number, m_kappa, m_T, m_mach_vector);
    }
 
    p2m_SLP create_p2m_SLP() const
    {
        return p2m_SLP(m_wave_number, m_kappa, m_T, m_mach_vector);
    }
 
    p2m_DLP create_p2m_DLP() const
    {
        return p2m_DLP(m_wave_number, m_kappa, m_T, m_mach_vector);
    }
 
    l2p create_l2p() const
    {
        return l2p(m_wave_number, m_kappa, m_T, m_mach_vector);
    }
 
    p2p_SLP create_p2p_SLP() const
    {
        return p2p_SLP(m_wave_number, m_mach_vector);
    }
 
    p2p_DLP create_p2p_DLP() const
    {
        return p2p_DLP(m_wave_number, m_mach_vector);
    }
 
private:
    wave_number_t m_wave_number;    // the wave number
    location_t m_mach_vector;       // the mach number vector (v / c)
    double m_mach_number;           // the scalar mach number
    wave_number_t m_kappa;          // the upscaled wave_number
    dmatrix3_t m_T;                 // the coordinate transform matrix
    double m_accuracy;              // the expansion accuracy
    std::vector<helmholtz_3d_hf_level_data> m_level_data_vector;
};
 
}
}
 
#endif // NIHU_CONVECTED_HELMHOLTZ_3D_HF_FMM_HPP_INCLUDED