Skip to content

Optimize gravity forward of prism layer - #674

Merged
santisoler merged 15 commits into
mainfrom
prism-layer-memory-fix
Jul 13, 2026
Merged

Optimize gravity forward of prism layer#674
santisoler merged 15 commits into
mainfrom
prism-layer-memory-fix

Conversation

@santisoler

@santisoler santisoler commented Apr 21, 2026

Copy link
Copy Markdown
Member

Re-implement the gravity method of the prism layer: rewrite the forward modeling computation of the gravity fields of a prism layer through a custom Numba function that iterates directly over the 1D horizontal coordinates of the layer. Replace the old implementation that used to build the whole set of prism boundaries and call the prism_gravity function. This new implementation doesn't require to store all those prisms in memory, significantly lowering the memory requirements, especially for large problems. Remove obsolete private function and method in layer.py. Also, reorganize prism forward modelling into a submodule: group together prism-related pieces of code inside a new harmonica._forward/prisms module, organized in a similar way as the one for ellipsoids.

Relevant issues/PRs:

This issue has been found by @aguspesce, who noticed a significant amount of memory being used while forward modeling a large prism layer.

Group together prism-related pieces of code inside a new
`harmonica._forward/prisms` module, organized in a similar way as the
one for `ellipsoids`.
Reimplement the forward modeling computation of the gravity fields of
a prism layer through a custom Numba function that iterates directly
over the 1D horizontal coordinates of the layer. Replace the old
implementation that used to build the whole set of prism boundaries and
call the `prism_gravity` function. This new implementation doesn't
require to store all those prisms in memory, significantly lowering the
memory requirements, especially for large problems.
@santisoler santisoler added this to the v0.8.0 milestone Apr 21, 2026
@santisoler
santisoler marked this pull request as ready for review July 13, 2026 18:01
@santisoler

santisoler commented Jul 13, 2026

Copy link
Copy Markdown
Member Author

I benchmarked this branch against main (e60e717) with this script, running the forward of a big prism layer on single computation point (just to speed up things, since the memory issue is coming not from the number of observation points, but from the allocation of prisms boundaries):

import numpy as np
import harmonica as hm
import bordado as bd

# Build prism layer
region = (-10e3, 10e3, -10e3, 10e3)
shape = (5000, 5000)
non_dimensional_coords = 100.0
easting, northing, surface = bd.grid_coordinates(
    region=region, shape=shape, non_dimensional_coords=non_dimensional_coords
)

properties = {"density": np.ones_like(easting)}
prism_layer = hm.prism_layer(
    (easting, northing), surface=surface, reference=0, properties=properties
)
print(prism_layer)

# Build observation point
coordinates = (0.0, 0.0, 500.0)

# Benchmark performance
prism_layer.prism_layer.gravity(coordinates, field="g_z")

I ran this script using:

command time -v python script.py

On main I obtained a peak usage of memory of:

Maximum resident set size (kbytes): 4505612

While on this branch (eeb94dd) I get:

Maximum resident set size (kbytes): 1877456

This shows a significant reduction in memory usage.

@santisoler
santisoler merged commit 92bd344 into main Jul 13, 2026
14 checks passed
@santisoler
santisoler deleted the prism-layer-memory-fix branch July 13, 2026 18:23
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant