In [1]:

import numpy as np
import exafmm.laplace as laplace

In [2]:

laplace.__doc__

Out[2]:

"exafmm's submodule for Laplace kernel"

1. create sources and targets¶

init_sources() takes in two arguments: coordinates of source bodies src_coords and charges (weights) of source bodies src_charges. Both should be numpy.ndarray. It returns a list of sources.

init_targets() only requires an array of coordinates trg_coords from input and returns a list of targets.

In [3]:

nsrcs = 100000
ntrgs = 200000

# generate random positions for particles
src_coords = np.random.random((nsrcs, 3))
trg_coords = np.random.random((ntrgs, 3))

# generate random charges for sources
src_charges = np.random.random(nsrcs)

In [4]:

# create a list of source instances
sources = laplace.init_sources(src_coords, src_charges)

# create a list of target instances
targets = laplace.init_targets(trg_coords)

sources and targets are two lists, the type of each element is exafmm.laplace.Body.

In [5]:

print(type(sources))
print(type(sources[0]), type(targets[0]))

<class 'list'>
<class 'exafmm.laplace.Body'> <class 'exafmm.laplace.Body'>

2. create a LaplaceFmm instance¶

To run FMM, we need to create a LaplaceFmm instance. The constructor takes in 3 integers, the expansion order p, max number of bodies per leaf ncrit, and the tree depth depth, and a string filename, the file name of the pre-computation matrix. Its default value is laplace_d_p[$p].dat. This file will be created during setup() call.

For now, the topology of the adaptive tree is only determined by ncrit. Changing depth does not affect the tree structure. We keep this interface depth for future use (a parameter that controls the domain decomposition when MPI is enabled).

In [6]:

fmm = laplace.LaplaceFmm(p=10, ncrit=200, filename="test_file.dat")

3. setup fmm¶

Given sources, targets and fmm, the function setup() handles three tasks:

build the tree
build interaction lists
pre-compute (or load) invariant matrices.

It returns an octree, whose type is exafmm.laplace.Tree.

In [7]:

tree = laplace.setup(sources, targets, fmm)

In [8]:

print(type(tree))

<class 'exafmm.laplace.Tree'>

The pre-computation file should be generated already.

In [9]:

ls *.dat

test_file.dat

4. evaluate¶

evaluate() triggers the evaluation and returns the potential and gradient as an numpy.ndarray of shape (ntrgs, 4). The $i$ -th row of trg_values starts with the potential value of the $i$ -th target, followed by its three gradient values.

In [10]:

trg_values = laplace.evaluate(tree, fmm)

In [11]:

trg_values[6]

Out[11]:

array([ 7298.79741989, -9117.70340717, -2316.54604176, -1612.97018776])

Set verbose to True to show timings.

In [12]:

laplace.clear_values(tree)
trg_values = laplace.evaluate(tree, fmm, True)

5. check accuracy (optional)¶

fmm.verify(tree.leafs) returns L2-norm the relative error of potential and gradient in a list, compared with the values calculated from direct method.

In [13]:

fmm.verify(tree.leafs)

Out[13]:

[1.599971542163853e-09, 1.5686554179176571e-07]

6. update charges of sources and run FMM iteratively¶

In [14]:

niters = 5

for i in range(niters):
    print('-'*10 + ' iteration {} '.format(i) + '-'*10)  # print divider between iterations
    
    src_charges = np.random.random(nsrcs)          # generate new random charges
    laplace.update_charges(tree, src_charges)      # update charges
    laplace.clear_values(tree)                     # clear values
    trg_values = laplace.evaluate(tree, fmm)       # evaluate potentials

    print("Error: ", fmm.verify(tree.leafs))       # check accuracy

---------- iteration 0 ----------
Error:  [1.618751446467221e-09, 1.586287849001751e-07]
---------- iteration 1 ----------
Error:  [1.6754867707208437e-09, 1.609071569273411e-07]
---------- iteration 2 ----------
Error:  [1.6481379491975236e-09, 1.6048968403703656e-07]
---------- iteration 3 ----------
Error:  [1.627084147780005e-09, 1.5683217437032964e-07]
---------- iteration 4 ----------
Error:  [1.5953579903620995e-09, 1.5678691639592906e-07]

In [ ]: