# At which grid size does FMM become 'worth it'?

**URL:** <https://bempp.discourse.group/t/at-which-grid-size-does-fmm-become-worth-it/50>\
**Category:** Bempp\
**Created:** [2 December 2020 09:31 UTC](https://bempp.discourse.group/t/at-which-grid-size-does-fmm-become-worth-it/50 "2020-12-02T09:31:02Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![Rasmus](https://avatars.discourse-cdn.com/v4/letter/r/7993a0/32.png) [@Rasmus](https://bempp.discourse.group/u/Rasmus)\
**Post date:** [2 December 2020 09:31 UTC](https://bempp.discourse.group/t/at-which-grid-size-does-fmm-become-worth-it/50/1 "2020-12-02T09:31:02Z")

</div>

Hi everyone. Thanks a lot for all the help so far.

I’m (still) doing Helmholtz scattering on a rigid object. I’m using somewhat large grid’s, and my simulations takes a long time to compute.

```auto
import time
import bempp.api
import numpy as np

grid = bempp.api.import_grid('some_grid.stl')

k = 100 # wave number 
space = bempp.api.function_space(grid, "P", 1) # piecewise linear needed due to the hypersingular operator

I = bempp.api.operators.boundary.sparse.identity(space, space, space)
M = bempp.api.operators.boundary.helmholtz.double_layer(space, space, space,k,assembler= "fmm") # double layer potential
E = bempp.api.operators.boundary.helmholtz.hypersingular(space, space, space, k,assembler= "fmm")
lhs = 0.5 * I - M - 1j/k * E

# point source far away from the boundary
r = 2
theta = 45 * np.pi/180
phi = 20 * np.pi/180
r_s = r * np.array([np.cos(theta) * np.cos(phi) , np.sin(theta) * np.cos(phi) , np.sin(phi)])

@bempp.api.complex_callable
def f(x, n, domain_index, result):
     dist = np.sqrt((r_s[0] - x[0])**2 + (r_s[1] - x[1])**2 + (r_s[2] - x[2])**2)
     result[0] = np.exp(1j * k * dist)/dist
rhs = bempp.api.GridFunction(space, fun=f)

t1 = time.time()
print('Solving equation system')
sol, info = bempp.api.linalg.gmres(lhs, rhs,tol=1E-4)

print('solve time ')
print(time.time() - t1)
print('Solver returned with ')
print(info)
print('\n\n'

```

I’m currently running the above on a cloud computer with 4 VCPU’s and 32 GB Ram.  
For a grid with ~20.000 elements the system logged a solve time of

```auto
solve time 
6649.893636465073
Solver returned with 
0

```

So a little under two hours.

Now if i instead don’t specify the assembler like

```auto
M = bempp.api.operators.boundary.helmholtz.double_layer(space, space, space,k) 
E = bempp.api.operators.boundary.helmholtz.hypersingular(space, space, space, k)

```

then the program finished in about half the time

```auto
solve time 
3013.474852323532
Solver returned with 
0

```

Which leads me to belive that either I’m utilizing the FMM wrong, or that it is only worth it for huge grids if the overhead is still greater than dense assembly for grids with over 20.000 elements (the first option being much more likely).  
In theory FMM should scale with O[\log (N) N] instead of O[N^2] right?

I don’t know the significance of this but if i interrupt the computations the last funtion call is almost always

```auto
  File "/home/rasmus/anaconda3/lib/python3.8/site-packages/Bempp_cl-0.2.2-py3.8.egg/bempp/api/fmm/exafmm.py", line 160, in evaluate
    result = self._module.evaluate(self._tree, self._fmm)

```

Am i doing something stupid in my code? Even if there is no solution i would still like to know when the FMM is useful or not, so i can write better simulations in the future.

Best regards  
Rasmus

---

<div class="post-metadata">

**Author:** ![Xavier\_V](https://avatars.discourse-cdn.com/v4/letter/x/ad7895/32.png) [@Xavier\_V](https://bempp.discourse.group/u/Xavier_V)\
**Post date:** [2 December 2020 17:40 UTC](https://bempp.discourse.group/t/at-which-grid-size-does-fmm-become-worth-it/50/2 "2020-12-02T17:40:07Z")

</div>

Hi Rasmus,  
I will try to help you through some lines of investigation. You are right twice the time for the FMM solver for 20k elements seems somewhat “weird”.

First of all I am wondering if the number of elements is sufficient. In other words, what is your mesh density criteria per wavelength. In addition is it possible for you to get a reference solution thanks to a dense assembler.

In any case, I would check first the convergence of the solution by displaying the residuals for both assemblers. Normally they are quite similar. If the solution is struggling to converge with the FMM so you will most likely have to increase the expansion order; 5 by default (a very important parameter to ensure a reliable expansion of FMM functions). The drawback is that leads to increase your memory requirements.

Best Regards.

Xavier V.

---

<div class="post-metadata">

**Author:** ![Rasmus](https://avatars.discourse-cdn.com/v4/letter/r/7993a0/32.png) [@Rasmus](https://bempp.discourse.group/u/Rasmus)\
**Post date:** [3 December 2020 09:17 UTC](https://bempp.discourse.group/t/at-which-grid-size-does-fmm-become-worth-it/50/3 "2020-12-03T09:17:44Z")

</div>

Hi Xavier,  
Thanks a lot for taking the time.

> [@Xavier\_V](#):
>
> First of all I am wondering if the number of elements is sufficient. In other words, what is your mesh density criteria per wavelength.

I’m trying to aim for ~8-10 elements per wavelength, and i would need a finer mesh than the one i tested the above with eventually (50-60k elements). But i would prefer to get the program to run a bit faster before feeding it with larger meshes. However the 10 elements per wavelength criteria should be satisfied for the tested wavenuber k = 100.

> [@Xavier\_V](#):
>
> In addition is it possible for you to get a reference solution thanks to a dense assembler.

I will have a look at it. What would i do if the solutions are dissimilar?

> [@Xavier\_V](#):
>
> If the solution is struggling to converge with the FMM so you will most likely have to increase the expansion order; 5 by default (a very important parameter to ensure a reliable expansion of FMM functions). The drawback is that leads to increase your memory

The compute cluster at my university are very ‘generous’ with memory so my current virtual machine has about 64 GB of RAM, of which only about 30% gets utilized, so Increasing memory usage is not a problem if it can increase convergence.

I tried looking at the BEMPP documentation and I can’t find how to increase the expansion order of the FMM. How do i change this parameter?

---

<div class="post-metadata">

**Author:** ![Rasmus](https://avatars.discourse-cdn.com/v4/letter/r/7993a0/32.png) [@Rasmus](https://bempp.discourse.group/u/Rasmus)\
**Post date:** [4 December 2020 10:26 UTC](https://bempp.discourse.group/t/at-which-grid-size-does-fmm-become-worth-it/50/4 "2020-12-04T10:26:44Z")

</div>

So i have done some investigations, and have turned to logging everything in a log file.

What i have found for meshes around 20k elements, is that the convergence of GMRES seems to be very slow, for both dense assembly and FMM. generally the number of iterations and values of the residuals seem to be exactly the same. For high values of k the iteration count gets well into the 100’s. I came back this morning from having a test script run overnight, and the iteration counter was over 500.

I have played a little around with changing the expansion order by using

```auto
bempp.api.GLOBAL_PARAMETERS.fmm.expansion_order = X

```

at the default value of 5 each FMM call takes about 5-6 seconds

```auto
TIMING:bempp:Start operation: Evaluating Fmm.
TIMING:bempp:Start operation: Calling ExaFMM.
TIMING:bempp:Finished Operation: Calling ExaFMM.: 5.457153558731079s
TIMING:bempp:Start operation: Singular Corrections Evaluator
TIMING:bempp:Finished Operation: Singular Corrections Evaluator: 0.18956685066223145s
TIMING:bempp:Finished Operation: Evaluating Fmm.: 5.652909994125366s

```

For one GMRES iteration, 12 of these callse to FMM were logged.

by increasing it to say 6, my memory consumption increases a little but so do the EXAFMM calls

```auto
TIMING:bempp:Start operation: Evaluating Fmm.
TIMING:bempp:Start operation: Calling ExaFMM.
TIMING:bempp:Finished Operation: Calling ExaFMM.: 6.3357932567596436s
TIMING:bempp:Start operation: Singular Corrections Evaluator
TIMING:bempp:Finished Operation: Singular Corrections Evaluator: 0.22159862518310547s
TIMING:bempp:Finished Operation: Evaluating Fmm.: 6.563088893890381s

```

Which is slower per call, but now only 9 calls were made, so the total time should be similar.  
The residual also seem to converge at the same rate, here for k \approx 18, which is the low end of my frequency range of interest. (0.96 , 0.87 , 0.77 , 0.64, 0.56 , 0.52 , 0.48 , etc)

Extrapolating the curve, assuming an exponential decay, it seems it will take about 40 iterations to reach a tolerance of 1e-2 and around 70 for 1e-3.

Without using FMM, the GMRES goes a lot quicker like

```auto
INFO:bempp:GMRES Iteration 1 with residual 0.9608518350984073
INFO:bempp:GMRES Iteration 2 with residual 0.8774929392680748
INFO:bempp:GMRES Iteration 3 with residual 0.7720222395249918
INFO:bempp:GMRES Iteration 4 with residual 0.6491716894707288
...
INFO:bempp:GMRES Iteration 55 with residual 0.01023279292849304
INFO:bempp:GMRES Iteration 56 with residual 0.009307566374188364
INFO:bempp:GMRES finished in 56 iterations and took 5.87E+01 sec.

```

So in this example the non-FMM GMRES finishes in one minute, where the FMM GMRES takes around 1 minute for a _single_ GMRES iteration.

Meanwhile the place were FMM is actually faster than dense is when assembling the matrices with

```auto
TIMING:bempp:assemble_sparse : 3.238e+00s
INFO:bempp:OpenCL Device set to: pthread-AMD EPYC Processor (with IBPB)
TIMING:bempp:Start operation: Initialising Exafmm.
TIMING:bempp:Finished Operation: Initialising Exafmm.: 2.3232381343841553s
TIMING:bempp:Start operation: Singular assembler:helmholtz_double_layer_boundary:opencl
TIMING:bempp:Finished Operation: Singular assembler:helmholtz_double_layer_boundary:opencl: 5.6644837856292725s
DEBUG:bempp:Using cached Fmm Interface.
TIMING:bempp:Start operation: Singular assembler:helmholtz_hypersingular_boundary:opencl
TIMING:bempp:Finished Operation: Singular assembler:helmholtz_hypersingular_boundary:opencl: 4.836869716644287s
INFO:bempp:Starting GMRES iteration

```

Compared to

```auto
TIMING:bempp:Start operation: Regular assembler:helmholtz_double_layer_boundary:opencl
INFO:bempp:OpenCL Device set to: pthread-AMD EPYC Processor (with IBPB)
TIMING:bempp:Finished Operation: Regular assembler:helmholtz_double_layer_boundary:opencl: 279.1030373573303s
TIMING:bempp:Start operation: Singular assembler:helmholtz_double_layer_boundary:opencl
TIMING:bempp:Finished Operation: Singular assembler:helmholtz_double_layer_boundary:opencl: 5.6845927238464355s
TIMING:bempp:Start operation: Regular assembler:helmholtz_hypersingular_boundary:opencl
TIMING:bempp:Finished Operation: Regular assembler:helmholtz_hypersingular_boundary:opencl: 271.2179608345032s
TIMING:bempp:Start operation: Singular assembler:helmholtz_hypersingular_boundary:opencl
TIMING:bempp:Finished Operation: Singular assembler:helmholtz_hypersingular_boundary:opencl: 4.796561241149902s
INFO:bempp:Starting GMRES iteration

```

So it seems like it is obviously the GMRES solver that slows down using FMM, but i don’t understand why it should be over 50 times slower for similar convergence.

So my conclusion is that i’m probably using FMM/EXAFMM wrong.

I know that the variables

```auto
bempp.api.GLOBAL_PARAMETERS.fmm.dense_evaluation
bempp.api.GLOBAL_PARAMETERS.fmm.ncrit
bempp.api.GLOBAL_PARAMETERS.fmm.depth

```

are used in the EXAFMM calls, but i don’t know enough about FMM to know what values these should be.

Any of you know what is going on here?

---

<div class="post-metadata">

**Author:** ![timo](https://yyz2.discourse-cdn.com/free1/user_avatar/bempp.discourse.group/timo/32/7_2.png) [@timo](https://bempp.discourse.group/u/timo)\
**Post date:** [7 January 2021 13:41 UTC](https://bempp.discourse.group/t/at-which-grid-size-does-fmm-become-worth-it/50/5 "2021-01-07T13:41:15Z")

</div>

Hi Rasmus,

FMM performance is quite a complex topic. First of all, consider a simple single-layer Helmholtz or Laplace operator. Then one matvec requires one FMM pass. While the FMM scales close to linearly, the algorithm itself is quite complex. So even for a few ten thousand elements, as soon as you have the matrix available as dense assembly the matvec may be faster in dense mode due to the much simpler data transport. Obviously, there will be a crossover point and you will soon run out of memory anyway if you continue dense.

The next issue is the number of passes for a single operator. The conjugate double layer also only requires one pass, while the double layer operator requires three passess and the Helmholtz hypersingular requires six passes (Laplace hypersingular needs three passes). This further slows things down.

We use FMM for some very large problems, which cannot be represented in dense mode (think hundred thousands to millions of unknowns). But in personal experience, whenever a dense representation is possible and initial assembly is not far too slow, then that is preferable to FMM.