Geomodeling benchmark: the “Moureze”-Model

This model is part of a geomodeling benchmaring effort. More information (and, hopefully, publication) coming.

import os

These two lines are necessary only if gempy is not installed

# Importing gempy
import gempy as gp
import gempy_viewer as gpv

# Aux imports
import numpy as np
import pandas as pd

from gempy_engine.config import AvailableBackends

Loading surface points from repository:

With pandas we can do it directly from the web and with the right args we can directly tidy the data in gempy style:

data_path = os.path.abspath('../../data/input_data/Moureze')
Moureze_points = pd.read_csv(
    filepath_or_buffer=data_path + '/Moureze_Points.csv',
    sep=';',
    names=['X', 'Y', 'Z', 'G_x', 'G_y', 'G_z', '_'],
    header=0,
)

Sections_EW = pd.read_csv(
    filepath_or_buffer=data_path + '/Sections_EW.csv',
    sep=';',
    names=['X', 'Y', 'Z', 'ID', '_'], header=1
).dropna()

Sections_NS = pd.read_csv(
    filepath_or_buffer=data_path + '/Sections_NS.csv',
    sep=';',
    names=['X', 'Y', 'Z', 'ID', '_'], header=1
).dropna()

Extracting the orientatins:

mask_surfpoints = Moureze_points['G_x'] < -9999
surface_points = Moureze_points[mask_surfpoints][::10]
orientations = Moureze_points[~mask_surfpoints][::10]

Giving an arbitrary value name to the surface

surface_points['surface'] = '0'
orientations['surface'] = '0'

surface_points.tail()

	X	Y	Z	G_x	G_y	G_z	_	surface
3370	224.000000	65.335228	-120.0000	-99999.0	-99999.0	-99999.0	0.670953	0
3386	232.000000	20.000000	-85.1014	-99999.0	-99999.0	-99999.0	0.700194	0
3396	135.939468	154.308319	-154.0210	-99999.0	-99999.0	-99999.0	0.351625	0
3411	296.000000	378.000000	-156.4050	-99999.0	-99999.0	-99999.0	0.852947	0
3424	91.960609	151.813431	-83.2593	-99999.0	-99999.0	-99999.0	0.950689	0

orientations.tail()

	X	Y	Z	G_x	G_y	G_z	_	surface
3186	139.936676	207.712677	-154.766	-0.034360	0.068784	0.997040	0.050783	0
3254	109.926086	47.649914	-135.337	0.351890	-0.766284	-0.537571	0.942560	0
3315	158.341888	211.320877	-147.312	0.053863	0.143364	0.988203	0.053091	0
3371	0.734351	90.000000	-100.000	0.756724	0.533905	0.377245	0.170960	0
3409	47.970890	129.885971	-132.010	-0.453523	-0.849487	0.269609	0.016933	0

Data initialization:

Suggested size of the axis-aligned modeling box:

Origin: -5 -5 -200

Maximum: 305 405 -50

Suggested resolution: 2m (grid size 156 x 206 x 76)

Only using one orientation because otherwhise it gets a mess

Number voxels

np.array([156, 206, 76]).prod()

np.int64(2442336)

resolution_requ = [156, 206, 76]
resolution = [77, 103, 38]
resolution_low = [45, 51, 38]


surface_points_table: gp.data.SurfacePointsTable = gp.data.SurfacePointsTable.from_arrays(
    x=surface_points['X'].values,
    y=surface_points['Y'].values,
    z=surface_points['Z'].values,
    names=surface_points['surface'].values.astype(str)
)

orientations_table: gp.data.OrientationsTable = gp.data.OrientationsTable.from_arrays(
    x=orientations['X'].values,
    y=orientations['Y'].values,
    z=orientations['Z'].values,
    G_x=orientations['G_x'].values,
    G_y=orientations['G_y'].values,
    G_z=orientations['G_z'].values,
    names=orientations['surface'].values.astype(str),
    name_id_map=surface_points_table.name_id_map  # ! Make sure that ids and names are shared
)

structural_frame: gp.data.StructuralFrame = gp.data.StructuralFrame.from_data_tables(
    surface_points=surface_points_table,
    orientations=orientations_table
)

geo_model: gp.data.GeoModel = gp.create_geomodel(
    project_name='Moureze',
    extent=[-5, 305, -5, 405, -200, -50],
    # resolution=resolution_low,
    refinement=5,
    structural_frame=structural_frame
)

Now we can see how the data looks so far:

gpv.plot_2d(geo_model, direction='y')

<gempy_viewer.modules.plot_2d.visualization_2d.Plot2D object at 0x7fbc9ec7a3e0>

The default range is always the diagonal of the extent. Since in this model data is very close we will need to reduce the range to 5-10% of that value:

geo_model.interpolation_options.kernel_options.range *= 0.2
geo_model.interpolation_options.evaluation_options.verbose = True
geo_model.interpolation_options.evaluation_options.octree_error_threshold = 1.5
geo_model.interpolation_options.evaluation_options.number_octree_levels_surface = 5

gp.compute_model(
    gempy_model=geo_model,
    engine_config=gp.data.GemPyEngineConfig(
        use_gpu=False,
        dtype='float32',
        backend=AvailableBackends.PYTORCH
    )
)

• 
• 

Setting Backend To: AvailableBackends.PYTORCH
Number of voxels marked by stats: 437 of torch.Size([512]).
 Number of voxels marked by corners : 408
Total voxels: 437
Dense Grid would be 512 voxels
Chunking done: 30 chunks
Number of voxels marked by stats: 3053 of torch.Size([3736]).
 Number of voxels marked by corners : 1989
Total voxels: 3053
Dense Grid would be 4096 voxels
Chunking done: 25 chunks
Chunking done: 197 chunks
Chunking done: 35 chunks

Solutions: 5 Octree Levels, 1 DualContouringMeshes

Time

300k voxels 3.5k points

• Nvidia 2080: 500 ms ± 1.3 ms per loop (mean ± std. dev. of 7 runs, 1 loop each), Memory 1 Gb

• CPU 14.2 s ± 82.4 ms per loop (mean ± std. dev. of 7 runs, 1 loop each), Memory: 1.3 Gb

2.4 M voxels, 3.5k points

• CPU 2min 33s ± 216 ms per loop (mean ± std. dev. of 7 runs, 1 loop each) Memory: 1.3 GB

• Nvidia 2080: 1.92 s ± 6.74 ms per loop (mean ± std. dev. of 7 runs, 1 loop each) 1 Gb

2.4 M voxels, 3.5k points 3.5 k orientations

• Nvidia 2080: 2.53 s ± 1.31 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

gpv.plot_2d(geo_model, cell_number='mid', series_n=0, show_scalar=True)

<gempy_viewer.modules.plot_2d.visualization_2d.Plot2D object at 0x7fbc9ec7aec0>

gpv.plot_2d(geo_model, cell_number='mid', show_data=True, direction='y')

<gempy_viewer.modules.plot_2d.visualization_2d.Plot2D object at 0x7fbcb5a16890>

sphinx_gallery_thumbnail_number = 4

gpv.plot_3d(geo_model)

<gempy_viewer.modules.plot_3d.vista.GemPyToVista object at 0x7fbc9ec7add0>