Manual

This manual provides reference documentation to PETGEM from a user’s and developer’s perspective.

How generate documentation

PETGEM is documented using Sphinx and LaTeX. The documentation source are in the doc directory.

The following steps summarize how generate PETGEM documentation.

  1. Install Sphinx and LaTeX.
  2. If is necessary, edit the rst files in doc/ directory using your favorite text editor. Nothing fancy is needed since ReST format is really simple.
  3. (Re) generate the PETGEM documentation as follows:
  • By using setup.py script:

    $ python3 setup.py output_format
    

    since PETGEM actually support html and pdf formats, valid options for output_format = [pdfdocs, htmldocs].

  • By using Sphinx commands directly:

    $ cd doc
    $ make output_format
    

    where output_format = [html, latexpdf].

Coding style

PETGEM’s coding style is based on PEP-0008 guidelines. Main guidelines are the following:

  • 79 characteres per line.

  • 4 spaces per indentation level.

  • Variables are lower case meanwhile constants are upper case.

  • Comments convention for functions is as follows:

    def function(arg1, arg2):
        ''' This is a function.
    
            :param int arg1: array of dimensions ...
            :param str arg2: string that ...
        '''
    
  • The use of inline comments is sparingly.

  • Use of lowercase to name functions. Furthermore, functions names have following form: <action>_<subject>(), e.g. compute_matrix().

  • Use of whole words instead abbreviations, examples:

    • Yes: solve_system(), compute_edges(), compute_matrix().
    • No: solve(), compedges(), compmatrix().

PETGEM design

PETGEM use a code structure for the Nédelec FE algorithm that emphasizes good parallel scalability, which is crucial in the multi-core era. Furthermore, it’s modularity should simplify the process of reaching the best possible performance in terms of percentage of the peak amount of floating point operations provided by the architecture. The code structure is modular, simple and flexible which allows exploiting not just our own modules but also third party libraries. Therefore, the software stack includes interfaces to external suites of data structures and libraries that contain most of the necessary building blocks needed for programming large scale numerical applications, i.e. sparse matrices, vectors, iterative and direct solvers or domain decomposition. As a result, the code is compact and flexible whose main UML diagrams are described as follows.

Component diagram

Main components of PETGEM are described in Figure 7.1. Pre-processing, processing and post-processing phases are included in Figure 7.1.

PETGEM: component diagram.

Figure 7.1. PETGEM: class diagram.

Class diagram

Main classes for PETGEM are described in Figure 7.2. Among all, the kernel class is the most important since it manage and control the others classes, as consequence, kernel class is the start point for any modelling.

PETGEM: class diagram.

Figure 7.2. PETGEM: class diagram.

Sequence diagram

Figure 7.3 describe the sequence for a standard CSEM forward modelling.

PETGEM: sequence diagram.

Figure 7.3. PETGEM: sequence diagram.

PETGEM directory structure

This subsection is dedicated to list and describe the PETGEM directory structure.

Top directory structure
Name Description
doc/ Source files for PETGEM documentation
examples/ CSEM FM examples directory
petgem/ Source code
tests/ PETGEM tests directory
utils/ Set of Python/Matlab scripts for pre-processing and post-processing
AUTHORS Authors file
builder.py Build helper for PETGEM setup script. The original version of this script was adapted from NiPy project
config_site.py PETGEM configuration file
INSTALL Installation and configuration instructions
kernel.py Heart of the code. This file contains the entire work-flow for a CSEM FM
LICENSE License file
Makefile Makefile with main PETGEM setup commands
MANIFEST.in File with exact list of files to include in PETGEM distribution
README Readme file
setup.py Main PETGEM setup script, it is based on Python setup-module
VERSION File with PETGEM version

The PETGEM source code is petgem/, which has the following contents:

petgem/ directory structure
Name Description
base/ Common classes and functions for init process.
decomposition/ Modules for domain decomposition (required for parallel execution of PETGEM)
efem/ General modules and classes for describing a CSEM FM by using EFEM of lowest order in tetrahedral meshes
mesh/ Interface to import mesh files
monitoring/ Modules for PETGEM performance monitoring
parallel/ Modules for supporting parallel assembling and solution of CSEM FM
postprocessing/ Modules for postprocessing stage
solver/ Interface to various internal/external solvers

Running a simulation

The following command should be run in the top-level directory of the PETGEM source tree. PETGEM kernel is invoked as follows:

$ mpirun -n <# of MPI processes> python3 kernel.py -options_file path/petsc.opts path/input_file.py

where petsc.opts is the PETSc options file and input_file.py is the PETGEM parameters file, whose main structure is explained in the following section.

Parameters file description

By definition, any model should include: physical parameters, a mesh file, source and receivers files, computational issues (solver type, solver tolerance, preconditioner) and an output file format.

In sake of simplicity and in order to avoid a specific parser, the PETGEM parameters file is defined as a Pyhton dictionary. As consequence, the dictionary name and his key names MUST NOT BE changed. The content of this file is the following:

 modelling = {
  # -------------------------------
  # ----------- General -----------
  # -------------------------------

  # -------------------------------
  # ----- Pyshical parameters -----
  # -------------------------------
  # Source
  # Source frequency. Type: float
  # Optional: NO
  'freq': 1.0,

  # Source position(x, y, z). Type: float
  # Optional: NO
  'src_pos': [1750.0, 1750.0, -975.0],

  # Source orientation. Type: int
  # 1 = X-directed source
  # 2 = Y-directed source
  # 3 = Z-directed source
  # Optional: NO
  'src_direc': 1,

  # Source current. Type: float
  # Optional: NO
  'src_current': 1.0,

  # Source length. Type: float
  # Optional: NO
  'src_length': 1.0,

  # Conductivity model. Type: str
  # Optional: NO
  'sigma_file': 'examples/WHAM/elemsSigma.dat',

  # Background conductivity. Type: float
  # Optional: NO
  'sigma_background': 3.33,

  # -------------------------------
  # ------- Mesh information ------
  # -------------------------------
  # Mesh files

  # Nodes spatial coordinates. Type: str
  # Optional: NO
  'nodes_file': 'examples/WHAM/nodes.dat',

  # Elements-nodes connectivity. Type: str
  # Optional: NO
  'elemsN_file': 'examples/WHAM/elemsN.dat',

  # Elements-edges connectivity. Type: str
  # Optional: NO
  'elemsE_file': 'examples/WHAM/elemsE.dat',

  # Edges-nodes connectivity. Type: str
  # Optional: NO
  'edgesN_file': 'examples/WHAM/edgesN.dat',

  # nnz for matrix allocation (PETSc)
  'nnz_file': 'examples/WHAM/nnz.dat',

  # Boundary-edges. Type: str
  # Optional: NO
  'bEdges_file': 'examples/WHAM/bEdges.dat',

  # -------------------------------
  # ------------ Solver -----------
  # -------------------------------
  # Solver options must be set in
  # examples/WHAM/petsc_solver.opts

  # -------------------------------
  # ------------ Output -----------
  # -------------------------------
  # Name of the file that contains the receivers position. Type: str
  # Optional: NO
  'receivers_file': 'examples/WHAM/receivers.dat',
}

The PETGEM parameters file is divided into 4 sections, namely, physical parameters (frequency, source position, source current, source length, conductivity model, background conductivity), mesh information (file path of nodal spatial coordinates, nodal element connectivity, edge element connectivity, edges nodes connectivity, and sparsity structure for PETSc matrix allocation), solver parameters (solver type, tolerance, maximum number of iterations), results information (receivers position file path).

Mesh formats

PETGEM support following mesh formats:

  • STL format from NETGEN.
  • STL format from Gambit.
  • Msh format from Gmsh.
  • Matlab format: matrix of nodal coordinates and element connectivity.
  • Files in Alya format (dom.dat and geo.dat).

Aforementioned formats must be pre-processed using the petgem_preprocessing.m Matlab script. This script is included in utils/Preprocessing/petgem_preprocessing.m of the PETGEM source tree. From a Matlab terminal, petgem_preprocessing.m is invoked has follows:

$ petgem_preprocessing('mesh_file', 'mesh_format', sigma_domain_per_subdomain, 'receivers_file')

See Utilities section for more details about petgem_preprocessing.m.

Available solvers

This section describes solvers available in PETGEM from user’s perspective. Direct as well as iterative solvers and preconditioners are supported through an interface to third party libraries (PETSc, MUMPs).

PETGEM uses petsc4py package in order to support the Krylov Subspace Package (KSP) from PETSc. The object KSP provides an easy-to-use interface to the combination of a parallel Krylov iterative method and a preconditioner (PC) or a sequential direct solver. As result, PETGEM users can set various solver options and preconditioner options at runtime via the PETSc options database.

Simulations in parallel

In FEM or EFEM simulations, the need for efficient algorithms for assembling the matrices may be crucial, especially when the DOF is considerably large. This is the case for real scenarios of 3D CSEM forward modelling because the required accuracy. In such situation, assembly process remains a critical portion of the code optimization since solution of linear systems which asymptotically dominates in large-scale computing, could be done with linear solvers such as PETSc, MUMPs, PARDISO). In fact, in PETGEM V1.0, the system assembly takes around $40%$ of the total time.

The classical assembly in FEM or EFEM programming is based on a loop over the elements. Different techniques and algorithms have been proposed and nowadays is possible performing these computations at the same time, i.e., to compute them in parallel. However, parallel programming is not a trivial task in most programming languages, and demands a strong theoretical knowledge about the hardware architecture. Fortunately, Python presents user friendly solutions for parallel computations, namely, mpi4py. The mpi4py package is an open source software that provides bindings of the MPI standard for the Python programming language, allowing any Python code to exploit multiple processors architectures.

On top of that, Figure 7.4 depicts shows an upper view of the matrix assembly and solution using the mpi4py and petsc4py package in PETGEM. The first step is to partition the work-load into subdomains. This task can be done by Metis library, which makes load over processes balanced. After domain partition, subdomains are read and assigned to MPI tasks and the elemental matrices are calculated concurrently. These local contributions are then accumulated into the global matrix system. The process for global vector assembly is similar.

Subsequently, the system is ready to be solved. PETGEM uses the Krylov Subspace Package (KSP) from PETSc through the petsc4py package. The object KSP provides an easy-to-use interface to the combination of a parallel Krylov iterative method and a preconditioner (PC) or a sequential direct solver. As result, PETGEM users can set various solver options and preconditioner options at runtime via the PETSc options database. Since PETSc knows which portions of the matrix and vectors are locally owned by each processor, the post-processing task is also completed in parallel following the numerical scheme described in CSEM problem section.

All petsc4py classes and methods are called from the PETGEM kernel in a manner that allows a parallel matrix and parallel vectors to be created automatically when the code is run on many processors. Similarly, if only one processor is specified the code will run in a sequential mode. Although petsc4py allows control the way in which the matrices and vectors to be split across the processors on the architecture, PETGEM simply let petsc4py decide the local sizes in sake of computational flexibility. However, this can be modified in an easy way without any extra coding required.

Parallel scheme for assembly in PETGEM V1.0.

Figure 7.4. Parallel scheme for assembly and solution in PETGEM using 4 MPI tasks. Here the elemental matrices computation is done in parallel. After calculations the global system is built and solved in parallel using the petsc4py and mpi4py packages.

Visualization of results

Once a solution of CSEM forward modelling has been obtained, it should be post-processed by using a visualization program. PETGEM does not do the visualization by itself, but it generates output files (vtk format) with the final results. Figure 7.5 shows an example of PETGEM output for the modelling described in Examples section. Figure 7.5 was obtained using Paraview.

PETGEM vtk output.

Figure 7.5. PETGEM vtk output.

Utilities

A set of Matlab functions to preprocessing and postprocessing included. These scripts are located in utils/Preprocessing/ and utils/Postprocessing/ respectively. A more detailed explanation for Preprocessing is the following:

Preprocessing
Name Description Parameters
export_matrix_hdf5.m Export a matrix to a hdf5 file Matrix to be exported and out file name
petgem_preprocessing.m Build input files for PETGEM (PETSc format) File mesh to be loaded, mesh format, conductivity model (layer model), receivers file name

On the other hand, the Postprocessing utils are the following:

Postprocessing
Name Description Parameters
petgem_import.m Import PETGEM output to Matlab PETGEM output file name. This file contains 3 arrays: Primary field (\(Ep\)), Secondary field (\(Es\)) and Total field (\(Et\)) for each receiver.

Examples

This section includes a step-by-step walk-through of the process to solve a simple CSEM forward modelling. The typical process to solve a problem using PETGEM is followed: a model is meshed, PETGEM input files are created, a parameters file is drafted, PETGEM is run to solve the modelling and finally the results of the analysis are visualised.

Problem statement: Isotropic model

In order to explain the CSEM forward modelling using PETGEM, here is considered the canonical model by Weiss2006 which consists in four-layers: 1000 m thick seawater (3.3 \(S/m\)), 1000 m thick sediments (1 \(S/m\)), 100 m thick oil (0.01 \(S/m\)) and 1400 m thick sediments (1 \(S/m\)). The computational domain is a \([0,3500]^3\) m cube. For this model, a 1 Hz x-directed dipole source is located at \(z=975\), \(x=1750\), \(y=1750\). The receivers are placed in-line to the source position and along its orientation, directly above the seafloor (\(z = 990\)).

Meshing

PETGEM V1.0 is based on tetrahedral meshes of lowest order. Therefore, Figure 7.6 shows a 3D view of the model with its unstructured tetrahedral mesh for the halfspace \(y>1750\), with the color scale representing the electrical conductivity \(\sigma\) for each layer. Mesh in Figure 7.6 have been obtained using Gambit.

In-line canonical off-shore hydrocarbon model with its unstructured tetrahedral mesh for :math:`y>150`

Figure 7.6. In-line canonical off-shore hydrocarbon model with its unstructured tetrahedral mesh for \(y>1750\). The color scale represents the electrical conductivity \(\sigma\) for each layer.

PETGEM preprocessing

The next step in the process is the PETGEM input files construction. These files can be created using the utils/Preprocessing/petgem_preprocessing.m Matlab script. For this modelling, petgem_preprocessing.m should be invoked has follows:

$ petgem_preprocessing('examples/WHAM/Mesh/Mesh_WHAM.neu', 'Gambit', [1/.3 1 1/100 1], 'examples/WHAM/Mesh/Receivers_WHAM')

petgem_preprocessing.m will create 8 files in binary PETSc format: nodes.dat (nodal spatial coordinates), elemsN.dat (nodal elements connectivity), elemsE.dat (edges elements connectivity), edgesN.dat (edges nodes connectivity), bEdges.dat (edges boundary array), nnz.dat (sparsity pattern for matrix allocation), receivers.dat (receivers position).

Parameters file for PETGEM

The parameters file structure for PETGEM is well documented in Parameters file description section. The parameters file used for this example follows:

 modelling = {
  # -------------------------------
  # ----------- General -----------
  # -------------------------------

  # -------------------------------
  # ----- Pyshical parameters -----
  # -------------------------------
  # Source
  # Source frequency. Type: float
  # Optional: NO
  'freq': 1.0,

  # Source position(x, y, z). Type: float
  # Optional: NO
  'src_pos': [1750.0, 1750.0, -975.0],

  # Source orientation. Type: int
  # 1 = X-directed source
  # 2 = Y-directed source
  # 3 = Z-directed source
  # Optional: NO
  'src_direc': 1,

  # Source current. Type: float
  # Optional: NO
  'src_current': 1.0,

  # Source length. Type: float
  # Optional: NO
  'src_length': 1.0,

  # Conductivity model. Type: str
  # Optional: NO
  'sigma_file': 'examples/WHAM/elemsSigma.dat',

  # Background conductivity. Type: float
  # Optional: NO
  'sigma_background': 3.33,

  # -------------------------------
  # ------- Mesh information ------
  # -------------------------------
  # Mesh files

  # Nodes spatial coordinates. Type: str
  # Optional: NO
  'nodes_file': 'examples/WHAM/nodes.dat',

  # Elements-nodes connectivity. Type: str
  # Optional: NO
  'elemsN_file': 'examples/WHAM/elemsN.dat',

  # Elements-edges connectivity. Type: str
  # Optional: NO
  'elemsE_file': 'examples/WHAM/elemsE.dat',

  # Edges-nodes connectivity. Type: str
  # Optional: NO
  'edgesN_file': 'examples/WHAM/edgesN.dat',

  # nnz for matrix allocation (PETSc)
  'nnz_file': 'examples/WHAM/nnz.dat',

  # Boundary-edges. Type: str
  # Optional: NO
  'bEdges_file': 'examples/WHAM/bEdges.dat',

  # -------------------------------
  # ------------ Solver -----------
  # -------------------------------
  # Solver options must be set in
  # examples/WHAM/petsc_solver.opts

  # -------------------------------
  # ------------ Output -----------
  # -------------------------------
  # Name of the file that contains the receivers position. Type: str
  # Optional: NO
  'receivers_file': 'examples/WHAM/receivers.dat',
}

Note that you may wish to change the location of the input files to somewhere on your drive. By default PETGEM will create the output directory at same level where the parameters file is located. For this example and following the PETSc documentation, the solver options have been configured in the file as petsc.opts follows:

# Solver options for PETSc
-ksp_type gmres
-pc_type sor
-ksp_rtol 1e-8
-ksp_converged_reason
-log_summary

That’s it, we are now ready to solve the modelling.

Running PETGEM

To run the simulation, the following command should be run in the top-level directory of the PETGEM source tree:

$ mpirun -n 16 python3 kernel.py -options_file examples/WHAM/petsc.opts examples/WHAM/params_file.py

PETGEM solves the problem and outputs the solution to the output path (examples/WHAM/Output/). The output files will be PETSc binary format. By default PETGEM save the electric field components in different files:

  • EpX.dat, EpY.dat, EpZ.dat: primary electric field components.
  • EsX.dat, EsY.dat, EsZ.dat: secondary electric field components.
  • EtX.dat, EtY.dat, EtZ.dat: total electric field components.

PETGEM postprocessing

Once the simulation has ended, PETGEM output can be imported to Matlab by using the script utils/Postprocessing/petgem_import.m. Hence, in order to import the x-component of the total electric field for this modelling, petgem_preprocessing.m should be invoked has follows:

$ [EtX] = petgem_import('examples/WHAM/Output/EtX.dat')

The dimension of arrays of \(EtX\) is determined by the number of receivers. Once the arrays are imported, feel free to handling and ploting. Figure 7.7 shows a comparison of the x-component of total electric field between PETGEM results and the quasi-analytical solution obtained with the WHAM tool. In Figure 7.7 it is easy to see the effect of the imperfect absorbing boundaries which can be mitigated by enlargening the domain with element sizes increasing logarithmically outwards from the zone of interest. The total electric field in Figure 7.7 was calculated using a mesh with \(\approx12\) millions of edges (degrees of freedom).

In-line canonical off-shore hydrocarbon model with its unstructured tetrahedral mesh for :math:`y>150`

Figure 7.7. Total electric field comparative for x-component between PETGEM V1.0 and WHAM. The effect of our imperfect absorbing boundaries can be mitigated by enlargening the computational domain with element sizes increasing logarithmically outwards from the zone of interest.