custEM.core package

Submodules

custEM.core.model_base module

@author: Rochlitz.R

class custEM.core.model_base.MOD(mod_name, mesh_name, approach, debug_level=20, fenics_debug_level=50, mute=False, **main_kwargs)[source]

Bases: object

Main class for 3D CSEM modeling using FEniCS. After an instance of this mainclass is created, all modeling steps are conducted within this environement. Instances within the mainclass ‘MOD’ are:

Methods

  • solve_main_problem()

    function to solve the weak formulation of the main system, defined in the FE instance; note that iterative solvers are theoretically considered but not working properly yet

  • init_third_party_parameters()

    adjust a few dolfin and numpy parameters; so far, no need for customization

  • init_default_MOD_parameters()

    initialize default parameters which can be overwritten by specifying keyword arguments when calling the MOD instance

  • init_instances()

    initialize all required ‘sub’-classes of the MOD class

  • handle_exception()

    overwrite excepthook from sys module for customized logger

handle_exception(exc_type, exc_value, exc_traceback)[source]

Overwrite excepthook from sys module to enable raising uncaught exceptions by the customized self.logger from the logging module and write them to the “mod_name + _debug.log” file for debugging.

Required arguments

  • see description of Python exceptions and corresponding traceback

init_default_model_parameters()[source]

Initalizes default parameters of the MOD class which will be updated if keyword arguments are set when calling the MOD instance. The parameters are described in the MOD class documentation.

init_instances()[source]

Initalizes instances within MOD (FS, MP, PP, FE, Solver, IB). For more information, we refer to the main documentation of the sub-modules.

init_third_party_parameters(fenics_debug_level)[source]

Set some dolfin and print parameters. Might contain more specifications with custom choices in future.

Required arguments

  • fenics_debug_level, type int

    specified log level in FEniCS, default is 50, most levels refer to the ones of the standard Python library logging:

    10 or “DEBUG”, subdry 13 or “TRACE”, what’s happening (in detail) 16 or “PROGRESS”, what’s happening (broadly) 20 or “INFO”, information of general interest 30 or “WARNING”, things that may go boom later 40 or “ERROR”, things that go boom 50 or “CRITICAL”, ciritcal errors between go boom and beyond

solve_main_problem(bc=None, sym=True, out_of_core=False, convert=True, method='fgmres', delete_factorization=True, auto_interpolate=False, calculation_frequencies=[], secondary_potential=False, inv_dir=None, interpolation_quantities=['E_t', 'H_t'], remote_station=None, **post_proc_kwargs)[source]

Assemble and solve linear systems of equations. Some arguments are piped to the PostProcessingTD or PostProcessingFD classes with the post_proc_kwargs dictionary during the post-processing phase, if automated postprocessing is enabled (convert=True). For the documentation of valid keyword arguments for the PP.convert_results() method, it is referred to the post-processing classes.

Keyword arguments

  • solver = ‘MUMPS’, type int

    use direct solver MUMPS, no alternatives so far

  • bc = None, type str

    default boundary conditions (BC) are implicit Neumann bc; alternatively, ‘ZeroDirichlet’, short ‘ZD’, can be choosen or inhomogeneous Dirichlet conditions ‘ID’ in case of the frequency-domain E-field approach

  • sym = 1, 2 or True, type int or bool

    so far, all implemented systems of equations are symmetric, this parameter may be False in other FE formulations 1 stands for positive definite symmetric matrices. 2 stands for gerneral structurally symmetric matrices

  • out_of_core = False, type bool

    set True to enable ‘out_of_core’ option of direct solver MUMPS if main memory is exceeded. WARNING! may be REALY slow

  • convert = True, type bool

    automatically converts solution function U to real and imarginary parts of electric and magnetic fields and saves ‘xml’ and ‘pvd’ files. If False, functions convert_results() and export_all_results in the postprocessing instance PP can be called manually with custom conversion and export parameters

  • delete_factorization = True, type bool

    if False, keep solver instance with factorization of main system for later usage; usually not required

  • calculation_frequencies = [], type list of int

    only compute results for frequencies specified in this list; useful to recompute reminaing results after a random MUMPS crash

  • auto_interpolate = False, type bool

    set True to automatically interpolate solutions in parallel after the solution and conversion phase has been finished for a frequency (useful in multi-frequency simulations)

  • interpolate_quantities = [‘E_t’, ‘H_t’], type list of str

    list of quantities for parallel interpolation on the fly, choose from ‘E_t’, ‘H_t’, ‘E_s’, or ‘H_s’

  • method = tfqmr, type str

    if iterative solvers are considered in future, choose solver here

  • pc = ilu, type str

    if iterative solvers are considered in future, choose preconditioner here

  • secondary_potential = False, type bool

    set True to use secondary potential formulation for calculating DC fields, only relevant for DC approach, can also be specified before by build_var_form() call

  • post_proc_kwargs

    pipe post-processing keyword arugments to post-processing instance if automated conversion and export is enabled

custEM.core.post_proc_fd module

@author: Rochlitz.R

class custEM.core.post_proc_fd.PostProcessingFD(FE, export_domains, solvingapproach)[source]

Bases: object

Frequency-domain PostProcessing class for computing and exporting electric and magnetic fields, called from MOD instance.

Methods

  • convert_results()

    convert main solution to other field quantities, e.g., E to H

  • export_E_fields()

    export total and/or secondary electric field data

  • export_H_fields()

    export total and/or secondary magnetic field data

  • export_A_fields()

    export total and/or secondary potential ‘field’ data

  • export_all_results()

    export all available EM-field and/or potential data

  • save_field()

    utility function for writing data on hard drive

  • convert_H_to_E()

    convert magnetic fields obtained with magnetic field or potential approaches to electric fields

  • dc_post_proc()

    conduct post-porcessing for direct current DC modeling approach

convert_H_to_E(from_potential=False)[source]

Compute electric on basis of magnetic fields if enabled.

convert_results(fi=0, convert_to_H=True, convert_to_E=False, export_cg=False, export_pvd=True, export_nedelec=True, **dummy_kwargs)[source]

Automatically convert results from E-fields to H-fields or vice versa or from potentials to E or H, depending on the utilized approach.

Keyword arguments

  • fi = 0, type int

    integer specifiying to which frequency the current solution belongs

  • convert_to_H = True, type bool

    calculate H-fields from E-field or potential approaches

  • convert_to_E = True, type bool

    calculate E-fields, if H-field or F-U approach was used

  • export_cg = False, type bool

    set True for exporting calculated data on Lagrange spaces

  • export_pvd = True, type bool

    set True to export fields in .pvd format for Paraview

  • export_nedelec = True, type bool

    set True for exporting calculated data on a Nedelec spaces

dc_post_proc(export_cg=True, export_pvd=True, export_nedelec=True)[source]

Post-processing for direct-current (DC) approach.

Required arguments

  • export_cg, type bool

    specify if fields in data format should be exported directly after conversion

  • export_pvd, type bool

    specify if fields are also exported in .pvd format for Paraview

  • export_nedelec = True, type bool

    set True for exporting caluclated data on a Nedelec space

export_A_fields(export_pvd, export_cg, export_nedelec)[source]

Export potentials using specified data (xml/h5) and/or pvd format from either Nedelec space or Lagrange VectorFunctionSpace.

Required arguments

  • export_cg, type bool

    specify if fields in data format should be exported directly after conversion

  • export_pvd, type bool

    specify if fields are also exported in .pvd format for Paraview

  • export_nedelec = True, type bool

    set True for exporting caluclated data on a Nedelec space

export_E_fields(export_pvd, export_cg, export_nedelec)[source]

Export electric fields using specified data (xml/h5) and/or pvd format from either Nedelec space or Lagrange VectorFunctionSpace.

Required arguments

  • export_pvd, type bool

    specify if fields are also exported in .pvd format for Paraview

  • export_cg, type bool

    specify if fields in data format should be exported directly after conversion

  • export_nedelec = True, type bool

    set True for exporting caluclated data on a NedelecSpace

export_H_fields(export_pvd, export_cg, export_nedelec)[source]

Export magnetic fields using specified data (xml/h5) and/or pvd format from either Nedelec space or Lagrange VectorFunctionSpace.

Required arguments

  • export_pvd, type bool

    specify if fields are also exported in .pvd format for Paraview

  • export_cg, type bool

    specify if fields in data format should be exported directly after conversion

  • export_nedelec = True, type bool

    set True for exporting caluclated data on a Nedelec space

export_all_results(quantities='EAH', export_cg=False, export_pvd=True, export_nedelec=True)[source]

Export a selection of calculated electric, magnetic or potential fields using specified data (xml/h5) and/or pvd format.

Keyword arguments

  • quantities = None, type str

    if export_all is False, specify which data should be exported, use a combination of E, H and/or A combined in one string

  • export_cg = False, type bool

    specify if fields in data format should be exported directly after conversion

  • export_pvd = True, type bool

    specify if fields are also exported in .pvd format for Paraview

  • export_nedelec = True, type bool

    set True for exporting caluclated data on a Nedelec space

save_field(quant, q1, q2, export_pvd, export_cg, export_nedelec)[source]

Store fields on hard disc drive.

Required arguments

  • quant, type str

    string specifying the field quantitiy to be exported

  • q1, q2, type dolfin Function

    functions containing the results

  • export_cg, type bool

    specify if fields in data format should be exported directly after conversion

  • export_pvd, type bool

    specify if fields are also exported in .pvd format for Paraview

  • export_nedelec = True, type bool

    set True for exporting caluclated data on a Nedelec space

custEM.core.post_proc_td module

@author: Rochlitz.R

class custEM.core.post_proc_td.PostProcessingTD(FE, export_domains)[source]

Bases: object

Time-domain PostProcessing class for computing and exporting electric and magnetic fields, called from MOD instance. Note that interpolated results are currently only stored as numpy arrays. A method to export time-dependent solutions as .pvd file for Paraview might be implemented in future if required.

Methods

  • convert_results()

    convert frequency-domain solutions of FT approach to time domain, derive H from E in IE and RA approaches

  • export_td_data()

    save total electric and/or magnetic fields on hard disc drive

  • interpolate_frequency_solutions()

    interpolate frequency-domain results in FT approach

  • reorder_results()

    reorder interpolated frequency-domain results in FT approach

  • merge_td_results()

    combine all interpolated time-domain results in a common data structure and export all time-domain results of an interpolation mesh as numpy arrays (pvd files not supported yet)

convert_results(convert_to_H=True, export_cg=False, export_pvd=False, export_nedelec=True, interp_mesh=None, interpolate_fd=True, reorder_fd=True, transform_td=True)[source]

Automatically convert results from E-fields to H-fields.

Keyword arguments

  • convert_to_H = True, type bool

    calculate H-fields from E-fields

  • export_cg = False, type bool

    set True for exporting calculated data on Lagrange spaces

  • export_pvd = True, type bool

    set True to export fields in .pvd format for Paraview

  • export_nedelec = True, type bool

    set True for exporting calculated data on a Nedelec spaces

  • interp_mesh = None, type str

    specify interpolation mesh for automated FT approach post-processing

  • interpolate_fd = True, type bool

    enable automated interpolation of frequency-domain results in FT approach

  • reorder_fd = True, type bool

    enable automated reordering of interpolated frequency-domain results in FT approach

  • transform_td = True, type bool

    enable automated transformation of reordered frequency-domain results in FT approach

export_td_data(export_cg, export_pvd, export_nedelec)[source]

Export electric and or magnetic fields using specified data (xml/h5) and/or pvd format from Nedelec spaces.

Required arguments

  • export_pvd, type bool

    specify if fields are also exported in .pvd format for Paraview

  • export_cg, type bool

    specify if fields in data format should be exported directly after conversion

  • export_nedelec = True, type bool

    set True for exporting caluclated data on a NedelecSpace

interpolate_frequency_solutions(interp_mesh, EH)[source]

Interpolate frequency-domain data on specified interpolation mesh.

Required arguments

  • interp_mesh, type str

    target interpolation mesh

  • EH, type str

    by default, EH=’EH’, that means that both, electric and magnetic field results, are interpolated

merge_td_results(interp_mesh, EH='EH', interp_dir=None)[source]

Resort time-dependent results on interpolation meshes in a common data format for all time-domain modeling approaches and export data as numpy arrays.

Required arguments

  • interp_mesh, type str

    target interpolation mesh

Keyword arguments

  • EH, type str

    by default, EH=’EH’, that means that both, electric and magnetic field results, are sorted and exported

reorder_results(interp_mesh, EH)[source]

Reorder frequency-domain data station-wise by importing all results and storing them in one file per station, containing the values for all chosen frequencies.

Required arguments

  • interp_mesh, type str

    target interpolation mesh

  • EH, type str

    by default, EH=’EH’, that means that both, electric and magnetic field results, are reordered

custEM.core.pre_proc module

@author: Rochlitz.R

class custEM.core.pre_proc.PreProcessing(FS, MP, import_freq=None, field_selection='all', self_mode=False, fs_type='None')[source]

Bases: object

PreProcessing class called from MOD instance.

Methods

  • import_all_results()

    import results of an existing model - all quantities

  • import_selected_results()

    import results of an existing model - selected quantities

  • load()

    basic import function for data files

  • read_h5()

    utility import function for h5 data files

import_all_results(file_format)[source]

Load all existing quantities from already calculated solutions.

Required arguments

  • file_format, type str

    either ‘h5’ or ‘xml’, set within model paramers in MOD instance

import_selected_results(selection, file_format)[source]

Load selected quantities from already calculated solutions.

Required arguments:

  • selection, type str

    string that must contain a combination of E_t, E_s, H_t, H_s, A_t or A_s added to a continuous string

  • file_format, type str

    either ‘h5’ or ‘xml’, set within model paramers in MOD instance

load(quantity, file_format, switch, ri='data')[source]

Import data in either xml or h5 format.

Required arguments

  • quantitiy, type str

    quantitiy, e.g., E_t_real_cg, which should be imported

  • file_format, type str

    either ‘h5’ or ‘xml’, set within model paramers in MOD instance

  • switch, type int

    internally used switch, either 1 to use Lagrange VectorFunctionSpace or 2 for Nedelec space

Keyword arguments

-ri = ‘data’, type str

specify a function to be imported from the HDF5 data container, if it contains multiple results; used internally to distinguish between real and imag fields

read_h5(f_name, switch, ri)[source]

Read data files in ‘h5’ format.

Required arguments

  • fname, type str

    file name to import from, containing also the export path

  • switch, type int

    internally used switch, either 1 to use Lagrange VectorFunctionSpace or 2 for Nedelec space

-ri = ‘data’, type str

specify a function to be imported from the HDF5 data container, if it contains multiple results; used internally to distinguish between real and imag fields

custEM.core.solvers module

@author: Rochlitz.R

class custEM.core.solvers.Solver(FS, FE, mumps_debug=False, serial_ordering=False)[source]

Bases: object

Solver class called internally from MOD instance.

Methods

  • solve_var_form_default()

    solve variaitonal formulation using default FEniCS solver, no support of symmetry, memory limits

  • solve_system_default()

    solve assembled system of equations with default FEniCS solver, no support of symmetry, memory limits

  • solve_system_mumps()

    use MUMPS solver to solve system of equations, no memory limits and symmetry support

  • call_mumps()

    call MUMPS solver within solve_system_mumps() method

  • solve_var_form_iter()

    use iteative solvers to solve variational formulation, in development

  • solve_system_iter()

    use iteative solvers to solve system of equations, in development

BoundaryConditions

  • default choice of boundary conditions are implicit Neumann BC

  • for zero Dirichlet BC, change bc argument in the ‘solve’ functions or during assembly for all total field approaches

build_PRESB_Jacobian(Asys, Hh, B, A, inner_tol=0.1)[source]

Sets up the iterative solver and system to be solved. (GCR preconditioned with PRESB adjusted for the computation of entries of the Jacobian)

inner_tol = 1e-1 (default): set tolerance for inner solver

(in inv_base.py)

Function returns required variables (matrices, vectors and inner solver) for GCR method

Implemented and commented by M. Weiss (January 2024)

call_mumps(A, b, u, sym, out_of_core, first_call=True)[source]

Call MUMPS solver.

Required arguments

  • A = left-hand-side matrix, type PETScMatrix

    assembled LHS matrix

  • b = right-hand-side vector, type PETScVector

    assembled RHS vector

  • u = solution vector, type PETScVector

    empty solution vector to be filled

  • see the definitions of the solve_system_mumps() method for the description of the other input parameters

solve_PRESB_Jacobian(ksp, f1, f2, g, h, Ag, H, A11, Asystem, Aw, w, r, DD, AD, b, fs=None, flag_mt=False, outer_tol=0.0001, sol=None)[source]

Solves the system for multiple RHS using GCR preconditioned with PRESB to obtain the entries of the Jacobian (sensitivities)

outer_tol = 1e-4 (default): set tolerance for outer solver

(in inv_base.py)

Function returns sensitivites for one receiver (equivalent to a column of the Jacobian)

Function writes out convergence history (i.e. the residual for each iteration) to outerrelresPID.txt, where PID is the process ID

Implemented and commented by M. Weiss (January 2024)

solve_system_default(A, b)[source]

Use default petsc_lu solver to solve an assembled system of equations.

Required arguments

  • A = left-hand-side matrix, type PETScMatrix

    assembled LHS matrix

  • b = right-hand-side vector, type PETScVector

    assembled RHS vector

solve_system_iter_BD(Asys, Hh, b, B, A, inner, method, flag_mt, fs=None, serial_ordering=False, outer_tol=1e-08, inner_tol=0.001, maxiter=2000)[source]

Use generalised conjugate residual (GCR) method preconditioned by block-diagonal preconditioner to solve the system of equation (based on Grayver & Buerg (2014) Robust and scalable 3-D geo-electromagnetic modelling approach using the finite element method and Grayver & Kolev (2015) Large-scale 3D geoelectromagnetic modeling using parallel adaptive high-order finite element method)

The procedure is of “inner-outer” type, where the inner solver is a FGMRES method preconditioned with the multi-grid based auxiliary-space preconditioner AMS.

The implementation is largely based on available functions of the PETSc and hypre libraries.

The convergence history of the outer solver (i.e. the residual for each iteration) is written to outerrelres.txt.

Implemented and commented by M. Weiss (January 2024)

Required arguments (assembly of matrices in .py)

  • Asys = 2 by 2 system block matrix, type PETScMatrix

  • H = matrix to be solved by inner solver, type PETScMatrix

  • b = right-hand-side vector, type PETScVector

  • A = diagonal block of system matrix, type PETScMatrix

  • B = off-diagonal block of system matrix, type PETScMatrix

  • inner = ‘iterative’ or ‘directly’ depending on whether the inner systems

    are to be solved by an iterative or direct method (first argument in bundled input variable solvingapproach required to

    initialise a custEM model)

  • method = ‘fgmres’ by default, alternatively ‘cg’ or ‘gcr’

    (set when calling solve_main_problem)

  • flag_mt = False (default) or True (automatically set if MT problem is solved)

  • fs = None, type str or dolfin FunctionSpace

    define FunctionSpace for solution function(s)

  • serial_ordering = False (default) or True

    to use serial ordering of METIS instead of parallel ordering by SCOTCH; leads to slight increase of solution time for MUMPS but significantly reduces chance of potential MUMPS crashes

Keyword arguments

  • outer_tol = 1e-8, type float

    tolerance of outer GCR solver

  • inner_tol = 1e-3, type float

    tolerance of inner solver

  • maxiter = 2000, type int

    maximum number of iterations

solve_system_iter_Hfield(H, RHS, fs=None, method='fgmres')[source]

Computes the magnetic fields from the electric fields using the FGMRES method preconditioned with the auxiliary-space solver AMS.

Presents an alternative to deriving the magnetic fields by the direct solver MUMPS.

Can be employed by setting the 4th argument of the bundled input variable solvingapproach (when initialising a custem model) to ‘H_iterative’.

The convergence history of the solver (i.e. the residual for each iteration) is written to Hfielditer.txt.

Implemented and commented by M. Weiss (January 2024)

solve_system_iter_PRESB(Asys, Hh, b, B, A, inner, method, flag_mt, fs=None, serial_ordering=False, outer_tol=1e-08, inner_tol=0.001, maxiter=2000)[source]

Use generalised conjugate residual (GCR) method preconditioned by PRESB to solve the system of equation (based on Weiss et al. (2023) Iterative solution methods for 3D controlled-source electromagnetic modelling of geophysical exploration scenarios)

The procedure is of “inner-outer” type, where the inner solver is a FGMRES method preconditioned with the multi-grid based auxiliary-space preconditioner AMS.

The implementation is largely based on available functions of the PETSc and hypre libraries.

The convergence history of the outer solver (i.e. the residual for each iteration) is written to outerrelres.txt.

Implemented and commented by M. Weiss (January 2024)

Required arguments (assembly of matrices in .py)

  • Asys = 2 by 2 system block matrix, type PETScMatrix

  • H = matrix to be solved by inner solver, type PETScMatrix

  • b = right-hand-side vector, type PETScVector

  • A = diagonal block of system matrix, type PETScMatrix

  • B = off-diagonal block of system matrix, type PETScMatrix

  • inner = ‘iterative’ or ‘directly’ depending on whether the inner systems

    are to be solved by an iterative or direct method (first argument in bundled input variable solvingapproach required to

    initialise a custEM model)

  • method = ‘fgmres’ by default, alternatively ‘cg’ or ‘gcr’

    (set when calling solve_main_problem)

  • flag_mt = False (default) or True (automatically set if MT problem is solved)

  • fs = None, type str or dolfin FunctionSpace

    define FunctionSpace for solution function(s)

  • serial_ordering = False (default) or True

    to use serial ordering of METIS instead of parallel ordering by SCOTCH; leads to slight increase of solution time for MUMPS but significantly reduces chance of potential MUMPS crashes

Keyword arguments

  • outer_tol = 1e-8, type float

    tolerance of outer GCR solver

  • inner_tol = 1e-3, type float

    tolerance of inner solver

  • maxiter = 2000, type int

    maximum number of iterations

solve_system_mumps(A, b, fs=None, sym=False, out_of_core=False, first_call=True, tx_selection=None)[source]

Use MUMPS solver to solve an assembled linear system of equations.

Required arguments

  • A = left-hand-side matrix, type PETScMatrix

    assembled LHS matrix

  • b = right-hand-side vector, type PETScVector

    assembled RHS vector

Keyword arguments

  • fs = None, type str or dolfin FunctionSpace

    define FunctionSpace for solution function(s)

  • sym = False, type int or bool

    set to 1 for positive definite symmetric matrix set to 2 or True for general structurally symmetric matrix

  • out_of_core = False, type bool

    set to True if MUMPS should be allowed to use disc space for the solution process if the memory is completely exploited

  • tx_selection = None, type list of int

    specify if system should be solved for not all but only selected right-hand sides; so far used to solve DC system only for grounded Tx if there are also ungrounded ones in the same model

  • first_call = True, type bool

    specify if this is the first call of a solution based on the same system matrix, re-use factorization if first_call=**False**

solve_var_form_default(L, R, U, bc=None)[source]

Use default petsc_lu solver to solve a LinearVariationProblem.

Required arguments

  • L = left-hand-side, type UFL form

    not assembled! UFL form of the left-hand-side of the system of equations

  • R = right-hand-side, type UFL form

    not assembled! UFL form of the right-hand-side of the system of equations

  • U = Solution space, type dolfin FunctionSpace

    mixed solution FunctionSpace

Keyword arguments

  • bc = None, type str

    so far ZeroDirichlet or short ZD bc and implicit Neumann bc by using None are supported

solve_var_form_iter(L, R, U, bc, sym=False, method='gmres', pc='petsc_amg', abs_tol=1e-06, rel_tol=1e-06, maxiter=2000)[source]

Use Krylov solver to solve a LinearVariationProblem.

Required arguments

  • L = left-hand-side, type UFL form

    not assembled! UFL form of the left-hand-side of the system of equations

  • R = right-hand-side, type UFL form

    not assembled! UFL form of the right-hand-side of the system of equations

  • U = Solution space, type dolfin FunctionSpace

    mixed solution FunctionSpace

  • bc = boundary conditions, type str

    so far ‘ZeroDirichlet’ or short ‘ZD’ bc and implicit Neumann bc by using None are supported

Keyword arguments

  • sym = False, type bool

    set to True if a symmetric system of equations is solved

  • method = ‘gmres’, type str

    iterative solution method

  • pc = ‘petsc_amg’, type str

    prconditioning type

  • abs_tol = 1e-6, type float

    absolute tolereance

  • rel_tol = 1e-6, type float

    relative tolerance

  • maxiter = 2000, type int

    maximum number of iterations

Module contents

core

Submodules:

  • model_base for defining main model class

  • pre_proc for defining pre-processing class

  • post_proc_fd for defining frequency-domain post-processing class

  • post_proc_td for defining time-domain post-processing class

  • solver for solving linear systems of equations