NiHu  2.0
Evaluating the Helmholtz kernel using the Fast Multipole Method

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 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,

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;
NiHu::fmm::chebyshev_cluster
Cluster class of the black box fmm.
Definition: chebyshev_cluster.hpp:28
NiHu::fmm::helmholtz_2d_wb_fmm
the 2d wide band helmholetz fmm
Definition: helmholtz_2d_wb_fmm.hpp:43
NiHu::fmm::black_box_fmm
Black box FMM for a smooth kernel.
Definition: black_box_fmm.hpp:61
NiHu::fmm::cluster_tree< cluster_t >
NiHu::fmm::divide_num_nodes
class representing a cluster division based on number of nodes
Definition: divide.hpp:74
NiHu::fmm::create_fmm_matrix
auto create_fmm_matrix(fmm_operator_collection< CollOps... > const &collection, cluster_tree< Cluster > const &tree, interaction_lists const &lists)
Factory function to create fmm_matrix from an operator collection.
Definition: fmm_matrix.hpp:628
NiHu::fmm::interaction_lists
class storing the different interaction lists of a tree
Definition: lists.hpp:22