Export a geological model from GemPy to use in MOOSE

import gempy as gp

Creating a geological model

The procedure of generating a geological model is presented in detail in Chapter 1-1 of the GemPy tutorials, so it will only be briefly presented here

Initiate a model

geo_model = gp.create_model('tutorial_moose_exp')
data_path = 'https://raw.githubusercontent.com/cgre-aachen/gempy_data/master/'

Import data from CSV-files with setting the resolution and model extent

gp.init_data(geo_model, [0, 2000., 0, 2000., 0, 2000.], [50, 50, 80],
             path_o=data_path + "/data/input_data/tut_chapter1/simple_fault_model_orientations.csv",
             path_i=data_path + "/data/input_data/tut_chapter1/simple_fault_model_points.csv",
             default_values=True)

Out:

Active grids: ['regular']

tutorial_moose_exp  2021-04-18 11:39

present the units and series

geo_model.surfaces

	surface	series	order_surfaces	color	id
0	Shale	Default series	1	#015482	1
1	Sandstone_1	Default series	2	#9f0052	2
2	Siltstone	Default series	3	#ffbe00	3
3	Sandstone_2	Default series	4	#728f02	4
4	Main_Fault	Default series	5	#443988	5
5	basement	Basement	1	#ff3f20	6

combine units in series and make two series, as the fault needs its own

gp.map_stack_to_surfaces(geo_model,
                         {"Fault_Series": 'Main_Fault',
                          "Strat_Series": ('Sandstone_2', 'Siltstone', 'Shale', 'Sandstone_1', 'basement')},
                         remove_unused_series=True)

# set the fault series to be fault object
geo_model.set_is_fault(['Fault_Series'], change_color=False)

	order_series	BottomRelation	isActive	isFault	isFinite
Fault_Series	1	Fault	True	True	False
Strat_Series	2	Erosion	True	False	False

check whether series were assigned correctly

geo_model.surfaces

	surface	series	order_surfaces	color	id
4	Main_Fault	Fault_Series	1	#443988	1
0	Shale	Strat_Series	1	#015482	2
1	Sandstone_1	Strat_Series	2	#9f0052	3
2	Siltstone	Strat_Series	3	#ffbe00	4
3	Sandstone_2	Strat_Series	4	#728f02	5
5	basement	Strat_Series	5	#ff3f20	6

Model generation

After loading in the data, we set it up for interpolation and compute the model.

set up interpolator

gp.set_interpolator(geo_model,
                    compile_theano=True,
                    theano_optimizer='fast_compile',
                    verbose=[])

Out:

Setting kriging parameters to their default values.
Compiling theano function...
Level of Optimization:  fast_compile
Device:  cpu
Precision:  float64
Number of faults:  1
Compilation Done!
Kriging values:
                     values
range               3464.1
$C_o$            285714.29
drift equations     [3, 3]

<gempy.core.interpolator.InterpolatorModel object at 0x7fcb8b722e50>

compute the model

gp.compute_model(geo_model, compute_mesh=False);

Out:

Lithology ids
  [6. 6. 6. ... 2. 2. 2.]

have a look at the data and computed model

gp.plot_3d(geo_model)

Out:

<gempy.plot.vista.GemPyToVista object at 0x7fcb89256100>

Exporting the Model to MOOSE

The voxel-model above already is the same as a model discretized in a hexahedral grid, so my immediately be used as input in a simulation tool, e.g. MOOSE. For this, we need to access to the unit IDs assigned to each voxel in GemPy. The array containing these IDs is called lith_block.

ids = geo_model.solutions.lith_block
print(ids)

Out:

[6. 6. 6. ... 2. 2. 2.]

This array has the shape of (x,) and would be immediately useful, if GemPy and the chosen simulation code would populate a grid in the same way. Of course, however, that is not the case. This is why we have to restructure the lith_block array, so it can be read correctly by MOOSE.

model resolution

nx, ny, nz = geo_model.grid.regular_grid.resolution

# model extent
xmin, xmax, ymin, ymax, zmin, zmax = geo_model.grid.regular_grid.extent

These two parameters are important to, a) restructure lith_block, and b) write the input file for MOOSE correctly. For a), we need to reshape lith_block again to its three dimensions and re-flatten it in a MOOSE-conform way.

reshape to 3D array

units = ids.reshape((nx, ny, nz))
# flatten MOOSE conform
units = units.flatten('F')

The importance of nx, ny, nz is apparent from the cell above. But what about xmin, …, zmax?

A MOOSE input-file for mesh generation has the following syntax:

[MeshGenerators]
  [./gmg]
    type = GeneratedMeshGenerator
    dim = 3
    nx = 50
    ny = 50
    nz = 80
    xmin = 0.0
    xmax = 2000.0
    yim = 0.0
    ymax = 2000.0
    zmin = 0.0
    zmax = 2000.0
    block_id = '1 2 3 4 5 6'
    block_name = 'Main_Fault Sandstone_2 Siltstone Shale Sandstone_1 basement'
  [../]

  [./subdomains]
    type = ElementSubdomainIDGenerator
    input = gmg
    subdomain_ids = ' ' # here you paste the transformed lith_block vector
  [../]
[]

[Mesh]
  type = MeshGeneratorMesh
[]

So these parameters are required inputs in the [MeshGenerators] object in the MOOSE input file. GemPy has a method to directly create such an input file, stored in gempy.utils.export.py.

The following cell shows how to call the method:

# sphinx_gallery_thumbnail_path = '_static/GemPy_model_combined.png'
import gempy.utils.export as export
export.export_moose_input(geo_model, path='')

Out:

Successfully exported geological model as moose input to ./

This method automatically stores a file geo_model_units_moose_input.i at the specified path. Either this input file could be extended with parameters to directly run a simulation, or it is used just for creating a mesh. In the latter case, the next step would be, to run the compiled MOOSE executable witch the optional flag --mesh-only.

E.g. with using the PorousFlow module:

$path_to_moose/moose/modules/porous_flow/porous_flow-opt -i pct_voxel_mesh.i --mesh-only

How to compile MOOSE is described in their documentation.

The now generated mesh with the name geo_model_units_moose_input_in.e can be used as input for another MOOSE input file, which contains the main simulation parameters. To call the file with the grid, the following part has to be added in the MOOSE simulation input file:

[Mesh]
  file = geo_model_units_moose_input_in.e
[]

---

The final output of the simulation may also be such an .e, which can, for instance, be opened with paraview. A simulated temperature field (purely conductive) of the created model would look like this: