.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "Tutorials\04_POC-HFD_calculation.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note Click :ref:`here ` to download the full example code .. rst-class:: sphx-glr-example-title .. _sphx_glr_Tutorials_04_POC-HFD_calculation.py: 04 - Analysis of SHEMAT-Suite models ========================== The created geological models with gempy were exported as SHEMAT-Suite input files. `SHEMAT-Suite `_ [1] is a code for solving coupled heat transport in porous media. It is written in fortran and uses a finite differences scheme in a hexahedral grid. In this example, we will load a heat transport simulation from the base POC model we created in "Geological model creation and gravity simulation". We will demonstrate methods contained in OpenWF for loading the result file, displaying the parameters it contains and how to visualize these parameters. Finally, we will calculate the conductive heat flow and plot it. .. GENERATED FROM PYTHON SOURCE LINES 11-19 .. code-block:: default # import some libraries import numpy as np import sys sys.path.append('../../') import OpenWF.postprocessing as pp import matplotlib.pyplot as plt .. GENERATED FROM PYTHON SOURCE LINES 20-26 Simulation results from SHEMAT-Suite can be written in different file formats, such as VTK (Visualization Toolkit) HDF5 (Hierarchical Data Format), or PLT (TecPlot). In addition, optional outputs, such as ASCII-Files with comparison of simulated values to measured ones can be provided. Further, a status log of the simulation and other meta-files. A full coverage of possible output files is provided in the `SHEMAT-Suite wiki `_. In this tutorial, we will work with HDF5 files and VTK files. The majority of methods in OpenWF are tailored towards HDF5 files, which are smaller than their VTK relatives. However, there exists a powerful visualization code for python which builds upon vtk, called pyvista. We will briefly showcase its capabilities at the end of this tutorial. .. GENERATED FROM PYTHON SOURCE LINES 28-41 Load HDF5 file -------------- From the base POC model, we created a SHEMAT-Suite input file. This was then executed with the compiled SHEMAT-Suite code. As basic information: we look at conductive heat transport, i.e. no fluid flow, and heat transport is described by Fourier's law of heat conduction :math:`q = - \lambda \nabla T`. At the base of the model, the heat flow boundary condition is set to 72 mW/m:math:`^2`. OpenWF has a built in method for loading HDF5 files, though reading a file is a one-liner using the library ``h5py``. .. code-block:: python fid = h5py.File('../../models/SHEMAT-Suite_output/SHEMAT_PCT_base_model_final.h5') The file can be loaded in different states, among others for 'r' for read, 'a' for append, 'w' for write, etc. The ``read_hdf`` method in OpenWF lets the user also choose the state to load the HDF5 file. .. GENERATED FROM PYTHON SOURCE LINES 41-46 .. code-block:: default model_path = '../../models/SHEMAT-Suite_output/SHEMAT_PCT_base_model_temp_final.h5' fid = pp.read_hdf_file(model_path, write=False) .. GENERATED FROM PYTHON SOURCE LINES 47-48 To check the parameters stored in the HDF5 file, you can query the loaded h5py file for its keys, i.e. the "labels" of the data boxes stored in the HDF5 file. .. GENERATED FROM PYTHON SOURCE LINES 48-51 .. code-block:: default fid.keys() .. rst-class:: sphx-glr-script-out Out: .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 52-54 As some of these acronyms can have no meaning to new users, we implemented a method, specifically for SHEMAT-Suite generated HDF5 files to present information about the stored parameters: .. GENERATED FROM PYTHON SOURCE LINES 54-57 .. code-block:: default pp.available_parameters(fid) .. rst-class:: sphx-glr-script-out Out: .. code-block:: none {'comp': 'compressibility', 'delx': 'discretization in x direction in meter', 'dely': 'discretization in y direction in meter', 'delz': 'discretization in z direction in meter', 'df': '?', 'ec': '?', 'head': 'hydraulic potential in meter', 'itemp_bcd': '?', 'itemp_bcn': '?', 'kx': 'log-permeability (square meter) in x direction', 'ky': 'log-permeability (square meter) in y direction', 'kz': 'log-permeability (square meter) in z direction', 'lc': '?', 'lx': 'thermal conductivity in x direction in watt per meter and kelvin', 'ly': 'thermal conductivity in y direction in watt per meter and kelvin', 'lz': 'thermal conductivity in z direction in watt per meter and kelvin', 'por': 'porosity', 'pres': 'pressure in megapascal', 'q': '?', 'rc': '?', 'rhof': 'density water in kilogram per cubic meter', 'temp': 'temperature in degrees celsius', 'temp_bcd': 'temperature dirichlet boundary condition in degrees celsius', 'temp_bcn': 'temperature neumann boundary condition in degrees celsius', 'uindex': 'rock unit index - geological unit present in the cell', 'visf': 'fluid viscosity', 'vx': 'velocity in x direction in meters per second', 'vy': 'velocity in y direction in meters per second', 'vz': 'velocity in z direction in meters per second', 'x': 'x coordinate in meters', 'y': 'y coordinate in meters', 'z': 'z coordinate in meters'} .. GENERATED FROM PYTHON SOURCE LINES 58-61 The postprocessing in OpenWF has methods for quickly displaying the parameters in each of the model dimensions in a 2D slice. For instance, we will look at a profile through the model parallel to the y direction, thus showing a crosssection of the model. In lines, it shows the interfaces between different geological units, and the specified parameter as a colored contour field. .. GENERATED FROM PYTHON SOURCE LINES 61-66 .. code-block:: default # plot slice temperature fig = plt.figure(figsize=[15,7]) pp.plot_slice(model_path, parameter='temp', direction='y', cell_number=25, model_depth=6500.) .. image:: /Tutorials/images/sphx_glr_04_POC-HFD_calculation_001.png :alt: temp,y-direction, cell 25 :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 67-73 .. code-block:: default # plot slice fluid density fig = plt.figure(figsize=[15,7]) pp.plot_slice(model_path, parameter='rhof', direction='y', cell_number=25, model_depth=6500) .. image:: /Tutorials/images/sphx_glr_04_POC-HFD_calculation_002.png :alt: rhof,y-direction, cell 25 :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 74-79 Heat flow estimation -------------------- SHEMAT-Suite does not provide the heat flow for HDF5 files. It does, however, store it in the vtk output. To also have the heat flow for HDF5 files, we provide a method for calculating it. In the future, it may be directly written in HDF5 files by SHEMAT-Suite. The method calculates the conductive heat flow in all model dimensions per default. The user can specify, if only one direction should be yielded by the method. .. GENERATED FROM PYTHON SOURCE LINES 79-82 .. code-block:: default qx, qy, qz = pp.calc_cond_hf(fid) .. GENERATED FROM PYTHON SOURCE LINES 83-89 Maybe here is a good point to talk about the dimensions and according directions in the HDF5 file. We see above, that qz has three dimensions, one with 60 entries, one with 50 entries, and one with 100 entries. These is also the cell discretization in z, y, and x direction. That is, in an HDF5 file from SHEMAT-Suite, we have the dimensions [z,y,x], so here qz[z,y,x]. The three variables now contain the heat flow in x, y, z direction for each cell in the model. Let's have a look at a horizontal slice through the model center and the heat flow in z-direction. Remembering the notation for directions, and seeing that in z-direction, we have 60 cells, the horizontal slice would reflect ``qz[29,:,:]``, i.e. all entries in x- and y-direction at index 29 of the z-direction. In the HDF5-file, we further count from the bottom, so ``qz[0,:,:]`` is the deepest slice, ``qz[-1,:,:]`` the shallowest. .. GENERATED FROM PYTHON SOURCE LINES 89-105 .. code-block:: default # Get the model dimensions in x, y, z x = fid['x'][0,0,:] y = fid['y'][0,:,0] z = fid['z'][:,0,0] cell_number = 29 ui_cs = fid['uindex'][cell_number,:,:] fig = plt.figure(figsize=[15,7]) cs = plt.contourf(x,y,qz[cell_number]*1000, 20, cmap='viridis') plt.contour(x,y,ui_cs, colors='k') plt.colorbar(cs, label='HF mW/m$^2$') plt.xlabel('X [m]') plt.ylabel('Y [m]'); .. image:: /Tutorials/images/sphx_glr_04_POC-HFD_calculation_003.png :alt: 04 POC HFD calculation :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out Out: .. code-block:: none Text(100.83333333333333, 0.5, 'Y [m]') .. GENERATED FROM PYTHON SOURCE LINES 106-108 Next to calculating the heatflow in each cell, we implemented a method to calculate it over a specified interval, e.g. over the depth interval of -4000 m to -3000 m, so providing the average heat flow over this depth interval. .. GENERATED FROM PYTHON SOURCE LINES 108-124 .. code-block:: default deeper = -4000 shallower = -3000 mid_depth = deeper - (deeper - shallower) / 2 qz_int = pp.calc_cond_hf_over_interval(fid, depth_interval=[deeper,shallower], model_depth=6500) fig = plt.figure(figsize=[15,7]) cs = plt.contourf(x,y,qz_int*1000, 20, cmap='viridis') plt.contour(x,y,ui_cs, colors='k') plt.colorbar(cs, label='HF mW/m$^2$') plt.title(f"Heat flow, at {mid_depth} m a.s.l. depth, over a {np.abs(deeper-shallower)} m depth interval") plt.xlabel('X [m]') plt.ylabel('Y [m]') plt.show(); .. image:: /Tutorials/images/sphx_glr_04_POC-HFD_calculation_004.png :alt: Heat flow, at -3500.0 m a.s.l. depth, over a 1000 m depth interval :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 125-128 VTK and pyvista --------------- `Pyvista `_ [2] is a python library for working with 3D meshes and providing an interface for VTK files. .. GENERATED FROM PYTHON SOURCE LINES 128-134 .. code-block:: default import pyvista as pv pv.set_plot_theme("document") simulation = pv.read('../../models/SHEMAT-Suite_output/SHEMAT_PCT_base_model_temp_final.vtk') .. GENERATED FROM PYTHON SOURCE LINES 135-136 This line loads the VTK file. For information about its content, we can simply call the variable: .. GENERATED FROM PYTHON SOURCE LINES 136-141 .. code-block:: default simulation simulation.plot(); .. image:: /Tutorials/images/sphx_glr_04_POC-HFD_calculation_005.png :alt: 04 POC HFD calculation :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out Out: .. code-block:: none [(49464.563148517394, 42464.563148517394, 39214.563148517394), (14000.0, 7000.0, 3750.0), (0.0, 0.0, 1.0)] .. GENERATED FROM PYTHON SOURCE LINES 142-143 The vtk file has a couple of scalar values stored (seen in the table with data arrays). We can switch the active scalars to temperature for example using: .. GENERATED FROM PYTHON SOURCE LINES 143-148 .. code-block:: default simulation.set_active_scalars('temp') simulation.plot(); .. image:: /Tutorials/images/sphx_glr_04_POC-HFD_calculation_006.png :alt: 04 POC HFD calculation :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out Out: .. code-block:: none [(49464.563148517394, 42464.563148517394, 39214.563148517394), (14000.0, 7000.0, 3750.0), (0.0, 0.0, 1.0)] .. GENERATED FROM PYTHON SOURCE LINES 149-153 References ---------- | [1] Keller, J., Rath, V., Bruckmann, J., Mottaghy, D., Clauser, C., Wolf, A., Seidler, R., Bücker, H.M., Klitzsch, N. (2020). SHEMAT-Suite: An open-source code for simulating flow, heat and species transport in porous media. SoftwareX, 12, 100533. | [2] Sullivan, C., & Kaszynski, A. (2019). PyVista: 3D plotting and mesh analysis through a streamlined interface for the Visualization Toolkit (VTK). Journal of Open Source Software, 4(37), 1450. .. rst-class:: sphx-glr-timing **Total running time of the script:** ( 0 minutes 11.709 seconds) .. _sphx_glr_download_Tutorials_04_POC-HFD_calculation.py: .. only :: html .. container:: sphx-glr-footer :class: sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: 04_POC-HFD_calculation.py <04_POC-HFD_calculation.py>` .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: 04_POC-HFD_calculation.ipynb <04_POC-HFD_calculation.ipynb>` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_