Evaluating the Helmholtz kernel using the Fast Multipole Method

Table of Contents

• Introduction
• Code
  • Type Definitions
  • Source and receiver nodes
  • Cluster tree
  • FMM initializaion
  • FMM operator manipulations
  • FMM matrix and MVP

Introduction

In this example we compute the acoustic field radiated by a set of 2D acoustic monopoles. This problem involves the evaluation of the 2d Helmholtz kernel between source and receiver locations and adding the sources' contributions. We solve this matrix multiplication by means of the wideband fast multipole method for the 2D Helmholtz kernel.

Code

The below sections present the full code needed to compute the MVP (matrix vector product).

Type Definitions

A few types are defined first fo convenience, including the cluster type, cluster tree type and physical location type.

typedef double wave_number_t;
typedef NiHu::fmm::helmholtz_2d_wb_fmm<wave_number_t> fmm_t;
typedef fmm_t::cluster_t cluster_t;
typedef NiHu::fmm::cluster_tree<cluster_t> cluster_tree_t;
typedef cluster_t::location_t location_t;

Source and receiver nodes

We define random source and receiver locations:

    size_t N_src = 5000;    // number of source nodes
    location_t *src_begin = new location_t[N_src];
    for (size_t i = 0; i < N_src; ++i)
        src_begin[i].setRandom();
    size_t N_rec = 3000;    // number of receiver nodes
    location_t *rec_begin = new location_t[N_rec];
    for (size_t i = 0; i < N_rec; ++i)
    {
        rec_begin[i].setRandom();
        rec_begin[i](0) += 2.0;
    }

Cluster tree

A cluster tree is built that contains all the source and receiver nodes. The clusters are divided until they contain less than 10 nodes.

    cluster_tree_t tree(src_begin, src_begin + N_src,
        rec_begin, rec_begin + N_rec,
        NiHu::fmm::divide_num_nodes(10));
    
    NiHu::fmm::interaction_lists lists(tree);

FMM initializaion

We initialize further data structures needed to compute the fmm operators. This is done by initializing the fmm object, and initializing its level data vector using the cluster tree. The level data vector is attached to the tree by linking the appropriate level data structures to each cluster.

    wave_number_t k = 2.0;
    fmm_t fmm(k);
 
    fmm.init_level_data(tree);
 
    for (size_t c = 0; c < tree.get_n_clusters(); ++c)
        tree[c].set_p_level_data(&fmm.get_level_data(tree[c].get_level()));

FMM operator manipulations

We get the fmm operators from the fmm object, and combine them into a collection for later convenience. Two functors are defined, first to assign the nodes and the cluster tree to the operators, and a second to assign the interaction lists to them.

Using the two functors, all the fmm operators are easily indexed and precomputed.

    auto ops = NiHu::fmm::create_fmm_operator_collection(
        fmm.create_p2p<0, 0>(),
        fmm.create_p2m<0>(),
        fmm.create_p2l<0>(),
        fmm.create_m2p<0>(),
        fmm.create_l2p<0>(),
        fmm.create_m2m(),
        fmm.create_l2l(),
        fmm.create_m2l()
    );
 
    auto p2p = fmm.create_p2p<0, 0>();
    std::cout << p2p(*rec_begin, *src_begin) << std::endl;
 
    auto idx_fctr = NiHu::fmm::create_indexed_functor(
        src_begin, src_begin + N_src,
        rec_begin, rec_begin + N_rec,
        tree);
 
    auto pre_fctr = NiHu::fmm::create_precompute_functor(tree, lists);
 
    auto pre_collection = ops.transform(idx_fctr).transform(pre_fctr);

FMM matrix and MVP

Finally, an fmm matrix is instantiated from the precomputed operators, a source vector is allocated, and the response is computed by simple matrix vector product.

    auto mat = NiHu::fmm::create_fmm_matrix(pre_collection, tree, lists);
 
    Eigen::Matrix<std::complex<double>, Eigen::Dynamic, 1> xct(N_src, 1);
    xct.setConstant(1.0);
 
    auto resp = mat * xct;