zgimbutas / fmmlib3d

Fast Multipole Method (FMM) library in R^3
BSD 3-Clause "New" or "Revised" License
39 stars 18 forks source link

Unexposed evaluation grid #2

Closed GavinBarill closed 4 years ago

GavinBarill commented 6 years ago

Hello, Thank you for this library. I am primarily interested in using this library for evaluating the influence of some point-dipole sources, with a goal of root finding an isovalue of the charge field. Because I'll be using root finding, I want to query the value of the field at multiple places in each iteration of my root-convergence. As such, for performance, I would be nice if I could save on the octree pre-computation. Ideally, I could save the local expansions of the finest grid level, so that I can evaluate the field each iteration without repeating pre-computation. Although I don't know what all of my evaluation points will be at the start of my root finding, I am able to put a bounding box within which I can guarantee my evaluation points will stay inside.

For example, here's some high level pseudocode of the ideal kind of interface for my problem:

evaluation_grid = precomputation(point_sources, bounding_box) for i = 0 ... until convergence function_values = fmmevaluation(evaluation_grid,evaluation_point) process function_values to create the next evaluation_point end

Personally I'm a complete novice at Fortran, but I'm good at c++ and matlab. Could you please provide any guidance on modifying your library to expose the pre-computation?

Thank you, Gavin

zgimbutas commented 6 years ago

Assuming that you are evaluating Laplace/Helmholtz interactions, the routines that need to be modified are in lfmm3dpart.f or hfmm3dpart.f, respectively. The oct-tree is stored in the wlists array (variable: real 8, allocatable :: wlists(:), the actual tree creation is done in fmm3dpartree routine, which is a driver for internal routine d3tstrcr); the multipole and local expansion due to sources are stored in wrmlexp array (variable: real *8, allocatable :: wrmlexp(:), with indexing defined in iaddr array, which is stored in w array).

Now, what you want to do is to rewrite fmm3dparttargmain routine to disable steps 6, 7, and 8, that perform evaluation on targets, and return w, wlists, and wrmlexp arrays from fmm3dparttarg routine. Currently, we are using Fortran 77 calling conventions for user convenience, returning pointers to dynamically allocated arrays will require using interface statements from Fortran 90 standard. Pointers in Fortran 90 are somewhat different beasts than in C++, but one can return something like them into C++ by using Fortran 90 to C compatibility packages. After that, an additional routine needs to be created to perform missing processing from steps 6, 7, and 8, by assigning targets to the oct-tree.

I don’t exactly know where your targets points are going to be located, there could be some extra problem here that will need to be addressed: the oct-tree is currently heavily pruned for efficiency reasons. Only oct-tree boxes that contain the initial sources and targets are preserved, therefore, there could be a situation where the new targets are outside the oct-tree while being inside the main bounding box. If that is the case, an extra tree processing will be necessary. If you are going to do the target evaluation multiple times, a relatively easy fix for this problem would be to switch on ifempty flag inside the oct-tree generation routine d3tstrcr to ifempty=1, forcing generation of oct-tree boxes and the corresponding local and multipole expansions, tiling the whole computational box. One note - you will still need full access to the oct-tree to finish your calculations: the local expansions describe the far field only due to well-separated boxes (step 7), you will need to add the direct contribution of neighboring boxes inside the oct-tree (step 8), as well as the contribution of small neighboring boxes via multipole expansions (step 6).

On Jan 4, 2018, at 2:40 PM, GavinBarill notifications@github.com wrote:

Hello, Thank you for this library. I am primarily interested in using this library for evaluating the influence of some point-dipole sources, with a goal of root finding an isovalue of the charge field. Because I'll be using root finding, I want to query the value of the field at multiple places in each iteration of my root-convergence. As such, for performance, I would be nice if I could save on the octree pre-computation. Ideally, I could save the local expansions of the finest grid level, so that I can evaluate the field each iteration without repeating pre-computation. Although I don't know what all of my evaluation points will be at the start of my root finding, I am able to put a bounding box within which I can guarantee my evaluation points will stay inside.

For example, here's some high level pseudocode of the ideal kind of interface for my problem:

evaluation_grid = precomputation(point_sources, bounding_box) for i = 0 ... until convergence function_values = fmmevaluation(evaluation_grid,evaluation_point) process function_values to create the next evaluation_point end

Personally I'm a complete novice at Fortran, but I'm good at c++ and matlab. Could you please provide any guidance on modifying your library to expose the pre-computation?

Thank you, Gavin

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/zgimbutas/fmmlib3d/issues/2, or mute the thread https://github.com/notifications/unsubscribe-auth/AGydmnKevLRo8BrE8ci_8qc40JenvxG3ks5tHUVqgaJpZM4RTvl7.