Model of Ales, France: Plotting Sections and Maps

Explanation

This model is generally quite unstable and requires float64 precision to determine a solution. The lack of data in one corner for the TRIAS and LIAS series causes the model to bend unrealistically, eroding the CARBO layer, which disappears in that section.

import gempy as gp
import gempy_viewer as gpv
import os

Setting up paths

cwd = os.getcwd()
data_path = os.path.join(cwd, 'examples') if 'examples' not in cwd else os.path.join(cwd, '../..')

path_interf = os.path.join(data_path, "data/input_data/AlesModel/2018_interf.csv")
path_orient = os.path.join(data_path, "data/input_data/AlesModel/2018_orient_clust_n_init5_0.csv")
path_dem = os.path.join(data_path, "data/input_data/AlesModel/_cropped_DEM_coarse.tif")

Creating the geological model

geo_model: gp.data.GeoModel = gp.create_geomodel(
    project_name='Claudius',
    extent=[729550.0, 751500.0, 1913500.0, 1923650.0, -1800.0, 800.0],
    resolution=None,
    refinement=6,
    importer_helper=gp.data.ImporterHelper(
        path_to_orientations=path_orient,
        path_to_surface_points=path_interf,
    )
)

Setting up the section grid

gp.set_section_grid(
    grid=geo_model.grid,
    section_dict={
            'section1': ([732000, 1916000], [745000, 1916000], [200, 150])
    }
)

Active grids: GridTypes.NONE|SECTIONS|OCTREE

	start	stop	resolution	dist
section1	[732000, 1916000]	[745000, 1916000]	[200, 150]	13000.0

Sorting lithologies

gp.map_stack_to_surfaces(
    gempy_model=geo_model,
    mapping_object={
            'fault_left': 'fault_left',
            'fault_right': 'fault_right',
            'fault_lr': 'fault_lr',
            'Trias_Series': ('TRIAS', 'LIAS'),
            'Carbon_Series': 'CARBO',
            'Basement_Series': 'basement'
    },
    remove_unused_series=True
)

Could not find element 'basement' in any group.

Structural Groups:	StructuralGroup: Name: fault_left Structural Relation: StackRelationType.ERODE Elements: StructuralElement: Name: fault_left StructuralGroup: Name: fault_right Structural Relation: StackRelationType.ERODE Elements: StructuralElement: Name: fault_right StructuralGroup: Name: fault_lr Structural Relation: StackRelationType.ERODE Elements: StructuralElement: Name: fault_lr StructuralGroup: Name: Trias_Series Structural Relation: StackRelationType.ERODE Elements: StructuralElement: Name: TRIAS StructuralElement: Name: LIAS StructuralGroup: Name: Carbon_Series Structural Relation: StackRelationType.ERODE Elements: StructuralElement: Name: CARBO
Fault Relations:	fault_left fault_righ... fault_lr Trias_Seri... Carbon_Ser... fault_left fault_right fault_lr Trias_Series Carbon_Series
	True False

Changing colors

geo_model.structural_frame.get_element_by_name("LIAS").color = "#015482"
geo_model.structural_frame.get_element_by_name("TRIAS").color = "#9f0052"
geo_model.structural_frame.get_element_by_name("CARBO").color = "#ffbe00"

Plotting the 2D model

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

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

Plotting section traces

gpv.plot_section_traces(geo_model)

<function plot_section_traces at 0x7fbd043bb010>

Setting faults

gp.set_is_fault(
    frame=geo_model.structural_frame,
    fault_groups=[
            geo_model.structural_frame.get_group_by_name('fault_left'),
            geo_model.structural_frame.get_group_by_name('fault_right'),
            geo_model.structural_frame.get_group_by_name('fault_lr')
    ],
    change_color=True
)

Structural Groups:	StructuralGroup: Name: fault_left Structural Relation: StackRelationType.FAULT Elements: StructuralElement: Name: fault_left StructuralGroup: Name: fault_right Structural Relation: StackRelationType.FAULT Elements: StructuralElement: Name: fault_right StructuralGroup: Name: fault_lr Structural Relation: StackRelationType.FAULT Elements: StructuralElement: Name: fault_lr StructuralGroup: Name: Trias_Series Structural Relation: StackRelationType.ERODE Elements: StructuralElement: Name: TRIAS StructuralElement: Name: LIAS StructuralGroup: Name: Carbon_Series Structural Relation: StackRelationType.ERODE Elements: StructuralElement: Name: CARBO
Fault Relations:	fault_left fault_righ... fault_lr Trias_Seri... Carbon_Ser... fault_left fault_right fault_lr Trias_Series Carbon_Series
	True False

Setting topography from file

gp.set_topography_from_file(
    grid=geo_model.grid,
    filepath=path_dem,
    crop_to_extent=[729550.0, 751500.0, 1913500.0, 1923650.0]
)

# Plotting 3D model with topography
gpv.plot_3d(geo_model, show_topography=True, ve=1, image=True)

Active grids: GridTypes.NONE|SECTIONS|TOPOGRAPHY|OCTREE

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

Getting the Carbon Series

carbo = geo_model.structural_frame.get_group_by_name("Carbon_Series")

Modifying interpolation options for better model fitting

geo_model.interpolation_options.number_octree_levels_surface = 4
geo_model.interpolation_options.kernel_options.range = 0.8
gp.modify_surface_points(
    geo_model=geo_model,
    elements_names=["CARBO", "LIAS", "TRIAS"],
    nugget=0.005
)

Structural Groups:	StructuralGroup: Name: fault_left Structural Relation: StackRelationType.FAULT Elements: StructuralElement: Name: fault_left StructuralGroup: Name: fault_right Structural Relation: StackRelationType.FAULT Elements: StructuralElement: Name: fault_right StructuralGroup: Name: fault_lr Structural Relation: StackRelationType.FAULT Elements: StructuralElement: Name: fault_lr StructuralGroup: Name: Trias_Series Structural Relation: StackRelationType.ERODE Elements: StructuralElement: Name: TRIAS StructuralElement: Name: LIAS StructuralGroup: Name: Carbon_Series Structural Relation: StackRelationType.ERODE Elements: StructuralElement: Name: CARBO
Fault Relations:	fault_left fault_righ... fault_lr Trias_Seri... Carbon_Ser... fault_left fault_right fault_lr Trias_Series Carbon_Series
	True False

Displaying the structural frame

print(geo_model.structural_frame)
geo_model.structural_frame

StructuralFrame(
        structural_groups=[
StructuralGroup(
        name=fault_left,
        structural_relation=StackRelationType.FAULT,
        elements=[
Element(
        name=fault_left,
        color=#527682,
        is_active=True
)
]
),
StructuralGroup(
        name=fault_right,
        structural_relation=StackRelationType.FAULT,
        elements=[
Element(
        name=fault_right,
        color=#527682,
        is_active=True
)
]
),
StructuralGroup(
        name=fault_lr,
        structural_relation=StackRelationType.FAULT,
        elements=[
Element(
        name=fault_lr,
        color=#527682,
        is_active=True
)
]
),
StructuralGroup(
        name=Trias_Series,
        structural_relation=StackRelationType.ERODE,
        elements=[
Element(
        name=TRIAS,
        color=#9f0052,
        is_active=True
),
Element(
        name=LIAS,
        color=#015482,
        is_active=True
)
]
),
StructuralGroup(
        name=Carbon_Series,
        structural_relation=StackRelationType.ERODE,
        elements=[
Element(
        name=CARBO,
        color=#ffbe00,
        is_active=True
)
]
)
],
        fault_relations=
[[False, False, False,  True,  True],
 [False, False, False,  True,  True],
 [False, False, False,  True,  True],
 [False, False, False, False, False],
 [False, False, False, False, False]],

Structural Groups:	StructuralGroup: Name: fault_left Structural Relation: StackRelationType.FAULT Elements: StructuralElement: Name: fault_left StructuralGroup: Name: fault_right Structural Relation: StackRelationType.FAULT Elements: StructuralElement: Name: fault_right StructuralGroup: Name: fault_lr Structural Relation: StackRelationType.FAULT Elements: StructuralElement: Name: fault_lr StructuralGroup: Name: Trias_Series Structural Relation: StackRelationType.ERODE Elements: StructuralElement: Name: TRIAS StructuralElement: Name: LIAS StructuralGroup: Name: Carbon_Series Structural Relation: StackRelationType.ERODE Elements: StructuralElement: Name: CARBO
Fault Relations:	fault_left fault_righ... fault_lr Trias_Seri... Carbon_Ser... fault_left fault_right fault_lr Trias_Series Carbon_Series
	True False

Explanation of model characteristics and adjustments This model has characteristics that make it difficult to get the right default values: - It is large, and we want high resolution - Some series have a large conditional number (i.e., the model input is not very stable) To address these issues: - Reduce the chunk size during evaluation to trade speed for memory - Reduce the std of the error parameter in octree refinement, which evaluates fewer voxels but may leave some without refinement Enable debugging options to help tune these parameters.

Setting verbose and condition number options for debugging

geo_model.interpolation_options.evaluation_options.verbose = True
geo_model.interpolation_options.kernel_options.compute_condition_number = True

Observations and parameter adjustments The octree refinement is making the octree grid almost dense, and smaller chunks are needed to avoid running out of memory. Adjusting parameters accordingly:

geo_model.interpolation_options.evaluation_options.octree_error_threshold = 0.5
geo_model.interpolation_options.evaluation_options.evaluation_chunk_size = 50_000

Computing the model with the adjusted settings

geo_model.interpolation_options.mesh_extraction = False
_ = gp.compute_model(
    geo_model,
    engine_config=gp.data.GemPyEngineConfig(
        backend=gp.data.AvailableBackends.PYTORCH,
        use_gpu=True,
        dtype="float64"
    )
)

• 
• 
• 
• 
• 
• 
• 
• 
• 
• 
• 
• 
• 
• 
• 

Setting Backend To: AvailableBackends.PYTORCH
Condition number: 2917573.7407638184.
Chunking done: 19 chunks
Condition number: 1174481.1989654184.
Chunking done: 9 chunks
Condition number: 1642063.3117937546.
Chunking done: 11 chunks
Condition number: 1586673.9519253937.
Chunking done: 248 chunks
Condition number: 1122611.6721707566.
Chunking done: 198 chunks
Chunking done: 32 chunks
Chunking done: 26 chunks
Number of voxels marked by stats: 136 of torch.Size([512]).
 Number of voxels marked by corners : 328
Total voxels: 136
Dense Grid would be 512 voxels
Number of voxels marked by stats: 129 of torch.Size([512]).
 Number of voxels marked by corners : 328
Total voxels: 186
Dense Grid would be 512 voxels
Number of voxels marked by stats: 152 of torch.Size([512]).
 Number of voxels marked by corners : 328
Total voxels: 242
Dense Grid would be 512 voxels
Number of voxels marked by stats: 136 of torch.Size([512]).
 Number of voxels marked by corners : 328
Total voxels: 314
Dense Grid would be 512 voxels
Number of voxels marked by stats: 170 of torch.Size([512]).
 Number of voxels marked by corners : 328
Total voxels: 367
Dense Grid would be 512 voxels
Chunking done: 27 chunks
Chunking done: 22 chunks
Chunking done: 16 chunks
Chunking done: 8 chunks
Chunking done: 9 chunks
Chunking done: 205 chunks
Chunking done: 164 chunks
Number of voxels marked by stats: 1014 of torch.Size([3360]).
 Number of voxels marked by corners : 1656
Total voxels: 1014
Dense Grid would be 4096 voxels
Number of voxels marked by stats: 976 of torch.Size([3360]).
 Number of voxels marked by corners : 1656
Total voxels: 1429
Dense Grid would be 4096 voxels
Number of voxels marked by stats: 1135 of torch.Size([3360]).
 Number of voxels marked by corners : 1656
Total voxels: 1852
Dense Grid would be 4096 voxels
Number of voxels marked by stats: 1051 of torch.Size([3360]).
 Number of voxels marked by corners : 1656
Total voxels: 2390
Dense Grid would be 4096 voxels
Number of voxels marked by stats: 1247 of torch.Size([3360]).
 Number of voxels marked by corners : 1656
Total voxels: 2741
Dense Grid would be 4096 voxels
Chunking done: 14 chunks
Chunking done: 7 chunks
Chunking done: 8 chunks
Chunking done: 184 chunks
Chunking done: 147 chunks
Chunking done: 108 chunks
Chunking done: 51 chunks
Chunking done: 62 chunks
Chunking done: 1459 chunks
Chunking done: 1167 chunks
Number of voxels marked by stats: 7902 of torch.Size([24048]).
 Number of voxels marked by corners : 8203
Total voxels: 7902
Dense Grid would be 32768 voxels
Number of voxels marked by stats: 7629 of torch.Size([24048]).
 Number of voxels marked by corners : 8203
Total voxels: 11219
Dense Grid would be 32768 voxels
Number of voxels marked by stats: 8784 of torch.Size([24048]).
 Number of voxels marked by corners : 8203
Total voxels: 14630
Dense Grid would be 32768 voxels
Number of voxels marked by stats: 8070 of torch.Size([24048]).
 Number of voxels marked by corners : 8203
Total voxels: 18800
Dense Grid would be 32768 voxels
Number of voxels marked by stats: 9624 of torch.Size([24048]).
 Number of voxels marked by corners : 8203
Total voxels: 21430
Dense Grid would be 32768 voxels
Chunking done: 102 chunks
Chunking done: 48 chunks
Chunking done: 59 chunks
Chunking done: 1378 chunks
Chunking done: 1103 chunks

Plotting the 2D model with and without topography

gpv.plot_2d(geo_model, show_topography=False, section_names=['topography'], show_lith=True)
gpv.plot_2d(geo_model, cell_number=[4], direction=['y'], show_topography=True, show_data=True)
gpv.plot_2d(geo_model, cell_number=[-4], direction=['y'], show_topography=True, show_data=True)

• 
• 
• 

/home/leguark/TeamCity/work/3a8738c25f60c3c9/venv/lib/python3.10/site-packages/gempy_viewer/API/_plot_2d_sections_api.py:106: UserWarning: Section contacts not implemented yet. We need to pass scalar field for the sections grid
  warnings.warn(

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

Setting thumbnail number for Sphinx-Gallery sphinx_gallery_thumbnail_number = -1

gpv.plot_3d(geo_model, show_lith=True, show_topography=True, kwargs_plot_structured_grid={'opacity': 0.8})

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