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:
FS.DOM:
custEM.fem.fem_base.Domains
Solver:
custEM.core.solvers.Solver
- FE - depending on the time- or frequency-domain modeling approach:
custEM.fem.time_domain_approaches.FourierTransformBased.
- IB - depending on time- or frequency-domain modeling:
- PP - depending on post-processing or pre-processing:
custEM.core.post_proc.PostProcessing
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