Adaptive Forward Modelling

[JonasLiebsch]

During my bachelor thesis I worked on gravity forward modelling using prisms. Since I focused on a terrain correction, the datasets used are very massive and caused long computation times. To accelerate the modelling, I decided to use an adaptive approach for the forward modelling of a regular spaced grid.
The idea is to use a high resolution only where it is needed.
To achieve that goal the algorithm performs the following steps:

• Create n resolution levels, level n has the highest resolution (lowest grid spacing) and is the input data which is used to compute the other levels. With each lower resolution level the grid spacing increases by a factor of 2.
• The actual forward modelling algorithm will start on level 0 (lowest resolution) and compute the gravity of each prism.
• This is repeated for level 1 with the higher resolution.
• Afterwards it compares the computed gravity values of the one large prism (level 0) and the four smaller prisms (level 1). If the difference is below a predefined threshold the value for that area is final and added to the overall result.
• In case the difference is higher than the threshold, the computation is repeated for level 2.
• Now level 1 and 2 get compared to decide where a higher resolution is needed.
• And so on ...
• The Process is finished if there is no more prism splitting or level n is reached

In my tests the algorithm was able to cut the runtime of a Terrain-/ Bouguer Correction on the Svalbard archipelago by a factor of 100 to 200, while creating errors smaller than 0.01 mgal (of course lower/higher thresholds could lead to more/less accurate results).
The algorithm was also useful to perform a simple isostatic correction, assuming airy isostasy.

Would such a feature be of interest for Harmonica?
How could this feature fit in the current structure of Harmonica? Perhaps a separate function besides the existing prism function and maybe a separate function for Bouger-/Terrain Correction which consider lakes and ice?
Of course the algorithm could be implemented similarly for Teseroids.
I am completely new to the Harmonica project and Git, so I am happy for any advice.

Are you willing to help implement and maintain this feature? Yes

[welcome]

👋 Thanks for opening your first issue here! Please make sure you filled out the template with as much detail as possible.

You might also want to take a look at our Contributing Guide and Code of Conduct.

[santisoler]

Hi @JonasLiebsch !
Sorry for the huge delay of this response 🙈

This is a very interesting algorithm, and we would love to have some way to reduce the computation time of terrain corrections and prisms layers forward modelling in Harmonica. Since yesterday we have a new function to create prisms layers as xarray.Datasets and an accessor class that enables to compute their gravity fields: https://www.fatiando.org/harmonica/dev/api/generated/harmonica.prisms_layer.html#harmonica.prisms_layer

I think a method like the one you developed should be accessible though that .gravity() method.

Perhaps a separate function besides the existing prism function and maybe a separate function for Bouger-/Terrain Correction which consider lakes and ice?

We decided not to include a function or class to specifically compute terrain corrections. Trying to design a piece of code that could be generally applied would be almost impossible, but we will encourage users to use prisms layers to define their own terrain models and use them to compute their gravitational fields, something like this: https://www.fatiando.org/harmonica/dev/gallery/forward/prisms_topo_gravity.html#sphx-glr-gallery-forward-prisms-topo-gravity-py

Regarding your algorithm, I have a question:
Isn't a way to avoid computing the gravitational effect of prisms of each level until you reach the level that fits your needs?
From my experience, the most computer demanding task is the actual forward modelling, so I would try to reduce the number of forward computations, even though the ones for lower resolutions will be faster than for the subsequent levels.

Here's my own though: what if instead of considering reducing the resolution of the whole layer of prisms, we reduce the resolution only on areas far from the computation point, and keep high resolution prisms near it?

What do you think?