.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples/co2_leakage_along_a_fault/2_generate_mesh_and_incon_files.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code .. rst-class:: sphx-glr-example-title .. _sphx_glr_examples_co2_leakage_along_a_fault_2_generate_mesh_and_incon_files.py: Generate MESH and INCON files ============================= In this example, we assume that the mesh has already been generated by a third-party software (here Gmsh via :mod:`pygmsh`). .. GENERATED FROM PYTHON SOURCE LINES 10-11 First, we import :mod:`numpy` and :mod:`toughio`. .. GENERATED FROM PYTHON SOURCE LINES 11-15 .. code-block:: Python import numpy as np import toughio .. GENERATED FROM PYTHON SOURCE LINES 18-19 A supported mesh can be read using the function :func:`toughio.read_mesh` that returns a :class:`toughio.Mesh` object. .. GENERATED FROM PYTHON SOURCE LINES 19-22 .. code-block:: Python mesh = toughio.read_mesh("mesh.vtu") .. GENERATED FROM PYTHON SOURCE LINES 25-27 (Optional) Although physical group names have been defined previously, most meshes do not export this information, including VTU. We can use the method :meth:`toughio.Mesh.add_material` to manually add material names. .. GENERATED FROM PYTHON SOURCE LINES 27-35 .. code-block:: Python mesh.add_material("UPPAQ", 1) mesh.add_material("CAPRO", 2) mesh.add_material("CENAQ", 3) mesh.add_material("BASAQ", 4) mesh.add_material("FAULT", 5) mesh.add_material("BOUND", 6) .. GENERATED FROM PYTHON SOURCE LINES 38-39 The mesh used in this sample problem is 2D and has been defined in the XY plane, but the points have 3D coordinates (with zeros as 3rd dimension for every cells). To make it 3D in the XZ plane, we swap the 2nd and 3rd dimensions, and then extrude the mesh by 1 meter along the Y axis (2nd dimension). .. GENERATED FROM PYTHON SOURCE LINES 39-43 .. code-block:: Python mesh.points[:, [1, 2]] = mesh.points[:, [2, 1]] mesh.extrude_to_3d(height=1.0, axis=1) .. GENERATED FROM PYTHON SOURCE LINES 46-49 (Optional) Before going any further, it is good practice to first check the quality of the mesh generated. TOUGH does not use any geometrical coordinate system and assumes that the line connecting a cell with its neighbor is orthogonal to their common interface. :mod:`toughio` provides a mesh property that measures the quality of a cell as the average absolute cosine angle between the line connecting a cell with its neighbor and the normal vector of the common interface. The mesh used in this example is of rather good quality. Bad quality cells are located at the boundaries of the model and mostly belong to the material `"BOUND"`. As this material is only used to impose Dirichlet boundary conditions in TOUGH, these bad quality cells will not impact the simulation outputs. .. GENERATED FROM PYTHON SOURCE LINES 49-66 .. code-block:: Python import pyvista pyvista.set_plot_theme("document") p = pyvista.Plotter(window_size=(800, 800)) p.add_mesh( mesh=mesh.to_pyvista(), scalars=mesh.qualities, scalar_bar_args={"title": "Average cell quality"}, clim=(0.0, 1.0), cmap="RdBu", show_scalar_bar=True, show_edges=True, ) p.view_xz() p.show() .. image-sg:: /examples/co2_leakage_along_a_fault/images/sphx_glr_2_generate_mesh_and_incon_files_001.png :alt: 2 generate mesh and incon files :srcset: /examples/co2_leakage_along_a_fault/images/sphx_glr_2_generate_mesh_and_incon_files_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 69-70 (Optional) Usually, a simple distribution plot is enough to rapidly assess the quality of a mesh. .. GENERATED FROM PYTHON SOURCE LINES 70-75 .. code-block:: Python import seaborn ax = seaborn.displot(mesh.qualities[mesh.materials != "BOUND"], kind="hist") .. image-sg:: /examples/co2_leakage_along_a_fault/images/sphx_glr_2_generate_mesh_and_incon_files_002.png :alt: 2 generate mesh and incon files :srcset: /examples/co2_leakage_along_a_fault/images/sphx_glr_2_generate_mesh_and_incon_files_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 78-80 We start by defining the boundary conditions. :mod:`toughio` recognizes the cell data key `"boundary_condition"` and automatically imposes Dirichlet boundary conditions to cells that have any value other than 0 in this cell data array. In this example, we simply set 1 to cells that belong to the group `"BOUND"` and 0 to others. .. GENERATED FROM PYTHON SOURCE LINES 80-85 .. code-block:: Python materials = mesh.materials bcond = (materials == "BOUND").astype(int) mesh.add_cell_data("boundary_condition", bcond) .. GENERATED FROM PYTHON SOURCE LINES 88-89 Initial conditions can be defined as a cell data array associated to key `"initial_condition"` where each column of the array corresponds to a primary variable. Note that :mod:`toughio` will not write any initial condition value that is lower than the threshold flag -1.0e9. .. GENERATED FROM PYTHON SOURCE LINES 89-98 .. code-block:: Python centers = mesh.centers incon = np.full((mesh.n_cells, 4), -1.0e9) incon[:, 0] = 1.0e5 - 9810.0 * centers[:, 2] incon[:, 1] = 0.05 incon[:, 2] = 0.0 incon[:, 3] = 10.0 - 0.025 * centers[:, 2] mesh.add_cell_data("initial_condition", incon) .. GENERATED FROM PYTHON SOURCE LINES 101-102 :mod:`toughio` also recognizes the cell data keys `"porosity"` and `"permeability"` in case we want to initialize porosity and/or permeability fields (e.g., if well logs data are available). Like boundary and initial conditions, we only have to associate new cell data arrays to keys `"porosity"` and/or `"permeability"`. The way these arrays are generated does not matter, they can be the results of simple interpolations (e.g., with :mod:`scipy`) or more advanced geostatistical interpolations (e.g., with :mod:`pykrige`). .. GENERATED FROM PYTHON SOURCE LINES 106-108 We can now write the `MESH` and `INCON` files by calling the method :meth:`toughio.Mesh.write_tough`. Additionally, we can also pickle the final mesh for later use (reading a pickle file is much faster than reading any mesh format). .. GENERATED FROM PYTHON SOURCE LINES 108-112 .. code-block:: Python mesh.write_tough("MESH", incon=True) mesh.write("mesh.pickle") .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 5.553 seconds) .. _sphx_glr_download_examples_co2_leakage_along_a_fault_2_generate_mesh_and_incon_files.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: 2_generate_mesh_and_incon_files.ipynb <2_generate_mesh_and_incon_files.ipynb>` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: 2_generate_mesh_and_incon_files.py <2_generate_mesh_and_incon_files.py>` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_