custEM.fem package

Submodules

custEM.fem.assembly module

@author: Rochlitz.R

class custEM.fem.assembly.MTAssembler(FE, bc)

Bases: object

Class that assembles the system matrix and the right-hand-sides (two source polarizations implemented as boundary conditions on the top

of the model domain) for natural-source modeling.

Methods

• assemble()

  assembly function of left- and right-hand sides for secondary-field formulations

assemble(fi=0)

class custEM.fem.assembly.SecondaryFieldAssembler(FE, bc)

Bases: object

Class that assembles the system matrix and the right-hand-sides for secondary-field formulations based on primary fields and delta_sigma.

Methods

• assemble()

  assembly function of left- and right-hand sides for secondary-field formulations

assemble(fi=0)

Let dolfin assemble the system matrix and right-hand side vectors.

Keyword arguments

• fi = 0, type int

  iteration index over frequencies for specifying the left-hand-side

class custEM.fem.assembly.TotalFieldAssembler(FE, bc, td_dc=False, force_p1=False)

Bases: object

Class for setting the assembled contribution of the current density on transmitter dof in the right-hand side vector for total field formulations.

Methods

• assemble()

  general assembly function of left- and right-hand sides for total-field formulations; calls all the other methods

• find_dof_coords()

  utility function searching for all Nedelec dof and whose coordinates which belong to the source defined on edges

• find_cells()

  for each dof which belongs to the source, search for an element which contains the former

• find_directions()

  evaluate the local direction of the Nedelec dof in the element which contains it and switch the signto match the global direction

• find_and_fill_start_stop_dof()

  required for the A_Phi approach on mixed elements to set divergence term on the start and stop point of grounded wire sources

• find_for_hed(), find_for_vmd()

  utility functions called by find_dof_coords()

• find_for_line(), find_for_loop_r(), find_for_path()

  utility functions called by find_dof_coords()

• append_lists()

  utility function to collect the transmitter dof

• eval_target_dir()

  utility function to find the global direction of an edge

• test_An_total(),

  method for testing total A-V formulation on nodal elements

• add_topo(),

  modify z-values of coordinates according to specified topography

add_topo(tx, tx_topo=None)

Apply topography information to the z-values of the transmitter coordinates.

Required arguements

• tx, type list of coordinates

  list of transmitter coordinates

Keyword arguments

• tx_topo = None, type str or function

  user might specify a custom topography reference here, not recommended, only implemented for tests

append_lists(idx, tdir=1.0, comp=0)

Utility function appending dof, coordinates, directions and component switches to lists during the dof-coordinate search.

Required arguements

• idx, type list of int

  integers specifying the transmitter dof

Keyword arguments

• tdir, type float

  either positive or negative 1

• comp = 0, type int

  either 0 or 1 for x or y component

assemble(fi=0)

Conduct right-hand side assembly of source currents manually and let dolfin assemble the system matrix.

Keyword arguments

• fi = 0, type int

  iteration index over frequencies for specifying the left-hand-side

eval_target_dir(a, b)

Evaluate global direction of each transmitter segment.

Required arguements

• a, b, type array/list of len 3

  start and stop coordinates to find the global positive or negative target direction, counter-clockwise in a right-handed coordinate system

find_and_fill_start_stop_dof(bb, current, ti=0, fi=0)

Required for the A-V approach on mixed elements and the DC approach to set divergence term on the start and stop point of grounded wire sources.

Required arguements

• bb, dolfin function

  assembled right-hand side function

• current, type float,

  source current for the corresponsing tx

Keyword arguments

• ti = 0, type int

  index of the corresponding tx

find_cells()

For each dof which belongs to the source, search for an element which contains this dof.

find_directions()

Evaluate the local direction of the Nedelec dof in the element, which contains it and switch the sign to match the global direction.

find_dof_coords()

Utility function searching for all Nedelec dof on edges and the corresponding coordinates which belong to the transmitter.

find_for_hed()

Find source dof and coordinates for HED source.

find_for_line()

Find source dof and coordinates for “line” source.

find_for_loop_r()

Find source dof and coordinates for “loop_r” source.

find_for_path()

Find source dof and coordinates for arbitrary grounded or ungrounded transmitters described by a path of coordinates.

find_for_vmd()

Find source dof and coordinates for VMD source.

test_An_total(fi)

Experimental total A-Phi nodal approach, use with caution and contact the authors for advise.

custEM.fem.calc_primary_boundary_fields module

@author: Rochlitz.R

class custEM.fem.calc_primary_boundary_fields.Pyhed_calc(mesh_name, results_dir, pf_type, file_format)

Bases: object

Primary field calculation class which starts calculating primary fields with pyhed-internal multiprocessing from a serial process to ensure correct dof ordering.

calc_fields()

convertCoordinates(gimli, dolfin)

Input: arr1, arr2: two coordinate lists of same shape (n, 3) which contains the same coordinates but in a diffrent order. output: arr1_arr2, arr2_arr1: index arrays which converts coordinates from input1 to input2 and from input2 to input1.

read_h5(var, fname)

read_serial_calculation_parameters()

write_h5(var, fname)

custEM.fem.calc_primary_fields module

@author: Rochlitz.R

class custEM.fem.calc_primary_fields.Pyhed_calc(mesh_name, results_dir, pf_type, file_format)

Bases: object

Primary field calculation class which starts calculating primary fields with pyhed-internal multiprocessing from a serial process.

calc_fields()

convertCoordinates(gimli, dolfin)

input: arr1, arr2: two coordinate lists of same shape (n, 3) which contains the same coordinates but in a diffrent order. output: arr1_arr2, arr2_arr1: index arrays which converts coordinates from input1 to input2 and from input2 to input1.

read_h5(var, fname)

read_serial_calculation_parameters()

write_h5(var, fname)

custEM.fem.fem_base module

@author: Rochlitz.R

class custEM.fem.fem_base.Domains(mesh, n_dom, domain_func=None, domain_vals=None)

Bases: object

Domain class called internally from FunctionSpace instance.

Methods

• add_anomaly()

  DEPRECATED - should not be used anymore! Add an anomaly domain in the ‘FEniCS-way’, e.g. predefined anomaly instances from the file misc/anomaly_expressions.py or user-defined dolfin expressions can be added to the model, each anomaly gets a new marker

add_anomaly(instance)

DEPRECATED! Do not use anymore without contacting the authors. (from custEM.misc import anomaly_expressions as ae) (MOD.FS.DOM.add_anomaly(ae.Plate(origin=[0.0, -100.0, -400.0])))

class custEM.fem.fem_base.FunctionSpaces(MP, p, test_mode, self_mode, ned_kind, dg_interp, load_existing, reuse_fs=None, overwrite_mesh=False, mesh_only=False)

Bases: object

FunctionSpaces class called from MOD instance. Used for initializing the required function spaces for the FE assembly, boundary condtions, and primary fields for all secondary field formulations.

Methods

• init_mesh()

  import mesh and domain information from .h5 file (default) in parallel or from .xml file in serial

• add_primary_fields()

  initialize PrimaryField instance, calculate or import primary fields

• add_dirichlet_bc()

  initalize dirichlet boundary conditions on the model boundary

• build_dg_func()

  build element-wise scalar parameter function (isotropic parameters)

• build_dg_tensor_func()

  build element-wise tensor parameter function (anisotropic parameters)

• build_complex_sigma_func()

  build complex-valued sigma functions (only isotropic) for modeling induced polarization effects

Description of initialization order and FunctionSpaces

1. import mesh (see init_mesh documentation)

2. build FEniCS FunctionSpaces depending on the mesh, polynomial order p and the chosen approach

   • S = Nodal (LaGrange) FunctionSpace

     order p and dimension 1

   • S_dg = element-wise scalar discontinuous FunctionSpace

     order 0, 1 value per cell

   • S_dgt = element-wise tensor-valued discontinuous FunctionSpace

     order 0, 9 values per cell

   • V = Nedelec Space

     order p and dimension 3

   • V_cg = Nodal (LaGrange) VectorFunctionSpace

     order p and dimension 3

   • W = mixed Nedelec Space

     order p and dimension 6 for coupled real and imaginary parts

   • W_cg = mixed (LaGrange) VectorFunctionSpace

     order p and dimension 6 for coupled real and imaginary parts

   • M = mixed solution FunctionSpace

     order p and dimension 6 or 8 (dependent on appraoch)

3. initalize empty solution Function U on mixed solution space

4. if a secondary field formulation is used, a PrimaryField instance PF is added to the FunctionSpace class later on by calling the add_Primary_Field() method automatically during the FE.build_var_form() call

add_dirichlet_bc(MP=None, dof_re=None, dof_im=None, bc_H=False)

Initalize Dichichlet BoundaryConditions on the solution space, only for internal usage. The keyword arguments are set automatically to handle different cases.

Keyword arguments

• MP = None, type class

  only used to apply inhomogeneous Dirichlet conditions in A-V approach, which requires E-field conversion to A depending on MP.omega

• dof_re = None, type list

  list of boundary dof of “real” part of FE system to apply inhomogeneous Dirichlet conditions

• dof_im = None, type list

  list of boundary dof of “imaginary” part of FE system to apply inhomogeneous Dirichlet conditions

• bc_H = False, type bool

  flag for testing inhomogeneous Dirichlet conditions for H-field conversion, provides no benefits so far but might be re-considered for other purposes in future

add_inner_dirichlet_bc(logger)

Initalize Dichichlet BoundaryConditions for inner edges of the mesh, used to specify “perfect conductor” infrastructure.

Required arguments ———-_——-

• inner_dof, type list

  list to specify the correct dof corresponding to the infrastructure implemented as edges in the mesh.

add_primary_fields(FE, MP, calc=True, **kwargs)

Add a PrimaryField instance to the FunctionSpace instance, called by FE if the ‘Secondary’ field formulation is used.

Required arguments

• FE, type class

  FE approach instance (e.g., E_vector)

• MP, type class

  ModelParameters instance

Keyword arguments

• calc = True, type bool

  enable primary fields calculation or PF class initialization only

• further kwargs

  see full list in the description of the PrimaryField class

build_complex_sigma_func(omega, sigma, tau, m, c)

Calculate and assign complex-valued conductivities for simulating induced polarization effects with a Cole-Cole model discontinous function spaces. This only works for isotropic conductivities for now.

Required arguments

• omega, type float

  angular frequency

• sigma, type list

  list of conductivity values

• tau, type list

  list of ip_tau values, see MP class documentation

• m, type list

  list of ip_m values, see MP class documentation

• c, type list

  list of ip_c values, see MP class documentation

build_dg_func(vals, inverse=False)

Assign isotropic conductivity values to discontinous function spaces.

Required arguments

• vals, type list (of lists)

  conductivity values specified in MP instance

Keyword arguments

• inverse = False, type bool

  return either sigma or 1/sigma (=inv(sigma))

build_dg_tensor_func(vals, inverse=False, logger=None)

Assign anisotropic conductivity values to discontinous tensor function spaces.

Required arguments

• vals, type list (of lists)

  conductivity values specified in MP instance

Keyword arguments

• inverse = False, type bool

  return either sigma or inv(sigma)

init_mesh(test_mode, MP, self_mode, overwrite_mesh)

Import the requested mesh and identify domain markers, if available.

Required arguments

• see documentation of MOD class

class custEM.fem.fem_base.ModelParameters(logger, str_args, test_mode, solvingapproach)

Bases: object

ModelParameter class called internally from MOD instance. This class contains, holds, and updates all physical and modeling parameters during the FE simulation. All model parameters can be updated with the update_model_parameters() method. A description of all available parameters of the ModelParameters class is listed below.

Methods

• update_model_parameters()

  update any supported physical or modeling parameter

• calc_delta_sigma()

  calculate delta_sigma as sigma_ground - sigma_0 and take reshape lists to account for any anisotropic value descriptions

• build_dg_parameter_functions()

  build discontinuous Parameter functions for all physical parameters, single-valued functions for isotropic values and tensor-valued functions for anisotropic values

• print_model_parameters()

  print updated list of model parameters

• load_mesh_parameteres()

  update mesh-related model parameters from ‘mesh-parameters’ json file

• check_parameter_dx_conformity()

  check conformity of number of domains and number of physical parameter values

Description of physical parameters

• f = 10., type float

  frequency (Hz)

• omega = 1e1, type float

  angular frequency, omega = 2 * pi * f

• frequencies = [], type list

  specify multiple frequencies (Hz)

• omegas = [], type list

  specify multiple angular frequencies, omegas = 2 * pi * frequencies

• n_freqs = 1, type int

  number of frequencies, calculated automatically

• current = 1., type float

  source current I (A) in the current carrying wire, variable J is deprecated but still used internally, current is the new alias for it

• currents = [], type 1D or 2D array or list of lists

  unique source currents I (A) for multiple Tx or multiple frequencies; In time-domain and for DC modeling, the first dimension of the current data array corresponds to n_tx; in frequency-domain modeling, the first dimension of the current data array corresponds to n_freqs and the second dimension of the current data array correponds to n_tx; if not specifed, current = 1. is used for all frequencies and Tx

• m_moment = 1., type float

  magnetic dipole moment for calculating background fields of a vmd source in a fullspace analytically

• n_tx = 1, type int

  number of transmitters (right-hand sides of the system of linear equations), calculated automatically

• sigma_air = 1e-8, type float or list

  homogeneous air-space conductivity, for modeling a fullspace, set sigma_air to the same value as sigma_ground

• sigma_ground = 1e-2, type list of float or list of list of floats

  subsurface conductivities corresponding to all domain markers (except 0 which corresponds to sigma_air); anisotropy can be specified by a list of three values (main diagonal) or a list of six values (upper triangle of 3x3 conductivity tensor) for each domain; a list of one value is handled identical to providing this single value as float instead of embedding it into a list

• sigma_0 = [], type list of float or list of list of floats

  background condcutivities for all sigma_ground, only required for secondary-field/potential formulations; arbitrary general anisotropic values are supported for sigma_ground and sigma_0; delta_sigma is evaluated internally

• sigma = [1e-2, 1e-8], type list of float

  list of all conductivities, the first two values in sigma are always sigma_ground (first layer) and sigma_air; if these values are updated later on, sigma will be updated internally; further values are related to subsurface layers or anomaly-domains in the model (see sigma_anom)

• delta_sigma = [0., 0.], type list of float

  conductivity differences, needed for secondary field formulations, updated internally

• anomaly_layer_markers = [], type list of int

  define the reference subsurface layer (from top to depth) for each anomaly to calculate delta_sigma. For instance, if two anomalies are located in the first and third layer of a three-layered earth, the correct choice would be [1, 3]; used for secondary field approaches only, this functionality will be updated by implementing an automated mapping in future

• mu = mu_0 = 4 * pi * 1e-7, type float

  magentic permeability (mu_r * mu_0), default is vacuum permeability

• mu_r = 1., type float or list of floats

  relative magnetic permeability

• eps = eps_0 = 8.85 * 1e-12, type float

  electric permittivity (eps_r * eps_0), default is vacuum permittivity

• eps_r = 1., type float or list of floats

  values of relative electric permittivity

• ip_tau = 1., type float or list of floats

  tau values of Cole-Cole model for modeling induced polarization effects

• ip_c = 1., type float or list of floats

  c values of Cole-Cole model for modeling induced polarization effects

• ip_m = 0., type float or list of floats

  m values of Cole-Cole model for modeling induced polarization effects

• system_it = False, type bool

  used to change system for iterative solver. Default is false. Set true if iterative solver is used

Description of mesh parameters

• n_layers = 1, type int

  number of subsurface layers, automatically imported from mesh-parameter file

• layer_depths = [], type list of float

  depths of subsurface-layer interfaces, automatically imported from mesh-parameter file

• layer_thicknesses = [], type list of float

  subsurface layer thicknesses, calculated automatically from layer_depths attribute

• topo = None, type str

  definition of DEM- or synthetic topography, imported automatically from mesh-parameter file.

build_dg_parameter_functions(FS, pf_EH=None, eps=False)

Build discontinuous isotropic or anisotropic parameter functions.

calc_delta_sigma()

Calculate delta sigma using sigma_ground and the corresponding background conductivities sigma_0.

check_parameter_dx_conformity()

Evaulate, if number of domains matches number of given physical parameter values.

load_mesh_parameters()

Import mesh parameters such as the number of domains or topography information.

print_model_parameters()

Print model parameters with adjusted style, e.g., no lengthly lists.

update_model_parameters(**mod_kwargs)

Updates all non-default physical parameters if called. A list of available parameters is given in the documentation of the ModelingParameters class.

custEM.fem.fem_utils module

@author: Rochlitz.R

class custEM.fem.fem_utils.ApproachBaseFD

Bases: object

Initialize FE variables for all frequency-domain approaches. All variables can be updated by specifying the corresponding keyword arguments in the build_var_form() method.

Methods

• update_fem_parameters()

  update valid parameters of the FE and PF classes and print updated class attributes

• update_tx_parameters()

  if s_type=’auto’, update Tx information from imported mesh-parameter file in the MP instance

Description of FE parameters

• s_type = ‘auto’, type str

  alternatives: hed (HED source, treated identical to line tx in halfspace

  and as infinitesimal source in fullspace with analytical solution)

  vmd (VMD source, treated as infinitesimal source in fullspace

  with analytical solution)

  heds (multiple HED sources, treated as infinitesimal sources

  in fullspace with analytical solution; specify origins with a list of 3D coordinated in the tx variable)

  vmds (multiple VMD sources, treated as infinitesimal sources

  in fullspace with analytical solution; specify origins with a list of 3D coordinated in the tx variable)

  line (real finite length dipole source) loop_c (circular loop source) loop_r (rectangular loop source) path (arbitrary sources)

• tx = None, type list

  internal usage only; if s_type=auto, this variable will be updated by the information in the mesh-parameter file; otherwise, the geometry parameters will be used to define the tx

• n_tx = 1, type int

  number of transmitters (right-hand sides of system of linear equations)

• grounding = False,

  internal usage only and only relevant if s_type=auto (this variable will be updated by the information in the mesh-parameter file)

• origin = [0.0, 0.0, 0.0], type list of len (3)

  origin of hed or vmd Tx; only valid if s_type is not auto and n_tx=1

• length = 1., type float

  length of hed Tx; only valid if s_type is not auto and n_tx=1

• self.azimuth = 0., type float

  azimuth of hed or vmd Tx in degree; only valid if s_type is not auto and n_tx=1

• start = [-100, 0.0, 0.0], type list of len(3)

  start coordinate of line Tx or lower left corner of loop_r Tx; only valid if s_type is not auto and n_tx=1; z_value is ignored and always set to zero

• stop = [100.0, 0.0, 0.0], type list of len(3)

  stop coordinate of line Tx or upper right corner of loop_r Tx; only valid if s_type is not auto and n_tx=1; z_value is ignored and always set to zero

• radius = 100., type float

  radius of loop_c Tx; only valid if s_type is not auto and n_tx=1

• pf_EH_flag = ‘E’, type str

  set to ‘H’ for using primary magentic fields to calculate the right-hand side source contributions in secondary field formulations

• bc = None, type str

  specify alternative boundary conditions (see FS class)

• quasi_static = True, type bool

  use diffusive (quasi static) approximation of the Helmholtz equation without displacement currents, False only implemented in the E_vector approach

• ip = False, type bool

  set True to enable modeling of induced polarization effects by specifying Cole-Cole model parameters (see MP class), only implemented in the E_vector approach

• tx_topo = None, type str or function

  specify topography (from file or as function z=f(x, y)) for searching the vertical Tx coordinates on a manually chosen topography for the total field approaches, not recommendend to be changed

• a_tol = 1e-2, type float

  find correct source dof locations for total field formulations; might be changed but the default choice is almost always sufficient

• assembly_time = 0., type float

  used internally, adds all assembly times to this variable during the to be stored by the profiler later on

update_fem_parameters(fem_kwargs)

Update and print FE and PF class attributes for frequency-domain approaches by overwriting the corresponding keyword arguments.

Required arguments

fem_kwargs, type dict

dictionary of keyword arguments for build_var_form() method

update_tx_parameters()

Update Tx information from imported mesh-parameter file in the MP instance if s_type=’auto’.

class custEM.fem.fem_utils.ApproachBaseTD

Bases: object

Initialize FE variables for all time-domain approaches. All variables can be updated by specifying the corresponding keyword arguments in the build_var_form() method.

Methods

• update_fem_parameters()

  update valid parameters of the FE and PF classes and print updated class attributes

• update_tx_parameters()

  if s_type=’auto’, update Tx information from imported mesh-parameter file in the MP instance

• calc_static_magnetic_field()

  calculate static magnetic field with the Biot-Savart law

• biot_savart()

  vectorized implementation of Biot-Savart law

• interpolate_static()

  interpolate static magentic field on specified interpolation mesh for independent post-processing purposes

• calc_dc_field()

  calculate direct current electric field with either total or secondary potential formulation of the common potential approach

• interpolate_dc()

  interpolate direct current field on specified interpolation mesh for independent post-processing purposes

• export_interpolation_data()

  export interpolated static magentic field or direct current electric field results as numpy files and, if specified, in .pvd file format

Description of FE parameters

Here, only parameters relevant for the time-domain approaches are listed. For all other parameters, it is referred to the documentation of the ApproachBaseFD class.

Note that further relevant parameters that can be updated with the build_var_form_method() are initialized by each time-domain approach class individually. It is referred to the corresponding documentation in the time_domain_approaches.py file.

• log_t_min = -6, type int/float

  specify start observation time of logarithmically divided interval for computing transients; defines the start time as 10^log_t_min

• log_t_max = 0, type int/float

  specify stop observation time of logarithmically divided interval for computing transients; defines the stop time as 10^log_t_max

• n_log_steps = 61, type int

  number of steps to logarithmically divide the time interval of the transient, equal to the number of results exported at the corresponding times

• times = None, type list

  manually define export times of the transients

• shut_off = True, type bool

  calculate either transient transmitter shut-off of shut-on response

biot_savart(rx, tx)

Vectorized implementation of the Biot-Savart law.

Required arguments

• rx, type list of coordinates

  receiver coordinates

• tx, type list of coordinates

  list of transmitter coordinates describing the transmitter path

calc_dc_field(dc_mesh=None, secondary_potential=False, store=False, td_dc=True)

Calculate direct current electric fields with either secondary or total field formulation of potential approach for obtaining shut-off transients.

Keyword arguments

• secondary_potential = False, type bool

  use either total or secondary potential formulation

• store = False, type bool

  save static magnetic fields as HDF5 file

calc_static_magnetic_field(store=False)

Calculate static magentic fields with the Biot-Savart law.

Keyword arguments

• store = False, type bool

  save static magnetic fields as HDF5 file

export_interpolation_data(field, V_serial, ti, export_pvd, interp_mesh)

Export interpolated results as complex numpy arrays and, if specified, as .pvd files to be viewed in Paraview.

Required arguments

• field, type dolfin Function

  dolfin function containing the field values

• V_serial, type dolin FunctionSpace

  non-distributed solution function space

• ti, type int

  number of the Tx to specify output file name

• export_pvd, type bool

  aside from numpy array, save interpolated fields also as .pvd file for Paraview

• interp_mesh, type str

  interpolation mesh name

interpolate_dc(interp_mesh, interp_p=None, export_pvd=True, interp_dir=None)

Customized interpolation function to interpolate direct current fields on specified interpolation mesh.

Required arguments

• interp_mesh, type str

  name of target interpolation mesh

Keyword arguments

• interp_p = None, type int

  default is interp_p=1, not reasonable to be changed

• export_pvd = True, type bool

  save fields as .pvd file for Paraview or not

• interp_dir = None, type str

  specify custom interpolation directory, not recommended

interpolate_static(interp_mesh, interp_p=None, export_pvd=True, interp_dir=None)

Customized interpolation function to interpolate static magnetic fields on specified interpolation mesh.

Required arguments

• interp_mesh, type str

  name of target interpolation mesh

Keyword arguments

• interp_p = None, type int

  default is interp_p=1, not reasonable to be changed

• export_pvd = True, type bool

  save fields as .pvd file for Paraview or not

• interp_dir = None, type str

  specify custom interpolation directory, not recommended

update_fem_parameters(fem_kwargs)

Update and print FE class attributes for time-domain approaches by overwriting the corresponding keyword arguments.

Required arguments

• fem_kwargs, type dict

  dictionary of keyword arguments for build_var_form() method

update_tx_parameters()

Update Tx information from imported mesh-parameter file in the MP instance if s_type=’auto’.

class custEM.fem.fem_utils.DC(FS, MP)

Bases: ApproachBaseFD

Implementation of common potential approach for caluclating direct current (DC) electric fields.

Methods

• build_var_form()

  build variational formulation with FEniCS

• lhs_form()

  build left-hand side of variational formulation

• rhs_form_total()

  build right-hand side of variational formulation for secondary potential approach

• rhs_form_secondary()

  build right-hand side of variational formulation for total potential approach

• calc_primary_potential()

  calculate primary potential for 1D backround resistivity distribution analytically

• calc_dc_field()

  customoized method for solving the DC FE problem

Description of DC parameters

Most parameters are equivalent to the parameters for the frequency-domain approaches in the ApproachBaseFD. Please refer to the corresponding documentaion. The only DC class relevant keyword argument to be updated with the build_var_form() method is listed below.

• secondary_potential = False, type bool

  use either total or secondary potential formulation

build_var_form(**fem_kwargs)

Update FE parameters and build variational formulations for the matrix assembly. For the description of valid keyword arguments, it is referred to the ApproachBaseFD and DC class documentations.

calc_dc_field(PP, profiler, secondary_potential=None, delete_factorization=True)

Assemble the system matrix and right-hand sides and solve the resulting linear system of equations for either the total or secondary potential approach. The electric DC field is also derived after the solution of the potential FE problem.

Keyword arguments

• secondary_potential = False, type bool

  use either total or secondary potential formulation, can be updated during calling this method is required, but should be already specified before

calc_primary_potential()

Calculate primary potential of static current analytically.

lhs_form()

Build left-hand side of variational formulation for DC approach.

rhs_form_secondary()

Build right-hand side of variational formulation for secondary DC potential approach.

rhs_form_total()

Build right-hand side of variational formulation for total potential approach.

custEM.fem.fem_utils.add_sigma_vals_to_list(vals, inverse=False)

Rearrange conductivity values to appropriate list format, e.g., convert isotropic values given as float to list of length 1.

Required arguments

• vals, type list (of lists)

  conductivity values specified in MP instance

Keyword arguments

• inverse = False, type bool

  return either sigma or 1/sigma (inv(sigma))

custEM.fem.fem_utils.check_sigma_vals(MP)

Evaulate, if anomalies are included or halfspace or fullspace model is chosen.

Required arguments

• MP, type class

  ModelParameters instance

custEM.fem.fem_utils.check_source(s_type, tx, pf_flag='E', approach=None)

Evaulate, which source type is chosen.

Required arguments

• s_type, type int

  source type (see ApproachBaseFD documentation)

• tx, type list of arrays

  list of transmitter coordinates for each Tx

Keywprd arguments

• pf_flag = ‘E’, type str

  either ‘E’ or ‘H’

custEM.fem.frequency_domain_approaches module

@author: Rochlitz.R

class custEM.fem.frequency_domain_approaches.A_V_mixed(FS, MP)

Bases: ApproachBaseFD

FE implementation of mixed potential approach, see (Ansari, 2014) or a modified formulation by (Schwarzbach, 2009).

build_var_form(**fem_kwargs)

Print FE parameters and initialize variational formulations for assembly. Valid keyword arguments are described in the ApproachBaseFD and PrimaryFields classes.

lhs_forms()

Left-hand side formulation as the basis for the system matrix assembly.

rhs_form_secondary(fi=0)

Right-hand side formulation for secondary potential field approach.

rhs_form_total()

Right-hand side formulation for total field approach. Right-hand side vector will be filled later by the TotalFieldAssembler instance.

class custEM.fem.frequency_domain_approaches.A_V_nodal(FS, MP)

Bases: ApproachBaseFD

FE implementation of nodal potential approach, see (Badea, 2001) or (Puzyref, 2013).

build_var_form(**fem_kwargs)

Print FE parameters and initialize variational formulations for assembly. Valid keyword arguments are described in the ApproachBaseFD and PrimaryFields classes.

lhs_form()

Left-hand side formulation as the basis for the system matrix assembly.

rhs_form_secondary(fi=0)

Right-hand side formulation for secondary field approach.

rhs_form_total()

Right-hand side formulation for total field approach. Right-hand side vector will be filled later by the TotalFieldAssembler instance.

class custEM.fem.frequency_domain_approaches.E_vector(FS, MP)

Bases: ApproachBaseFD

FE implementation of E-field approach (Grayver, 2014, Schwarzbach 2009).

build_var_form(**fem_kwargs)

Print FE parameters and initialize variational formulations for assembly. Valid keyword arguments are described in the ApproachBaseFD and PrimaryFields classes.

lhs_forms(system_it)

Left-hand side formulation as the basis for the system matrix assembly.

lhs_forms_full()

Left-hand side formulation as the basis for the system matrix assembly. Here, electric permittivities and induced polarization parameters are considered in addition to the conductivity values.

rhs_form_secondary(fi=0)

Right-hand side formulation for secondary electric field approach.

rhs_form_total(system_it)

Right-hand side formulation for total field approach. Right-hand side vector will be filled later by the TotalFieldAssembler instance.

class custEM.fem.frequency_domain_approaches.F_U_mixed(FS, MP)

Bases: ApproachBaseFD

FE implementation of Tau-Omega approach on mixed elements (Mitsuhata, 2004), called F-U potential approach in custEM.

build_var_form(**fem_kwargs)

Print FE parameters and initialize variational formulations for assembly. Valid keyword arguments are described in the ApproachBaseFD and PrimaryFields classes.

lhs_forms()

Left-hand side formulation as the basis for the system matrix assembly.

rhs_form_secondary(fi=0)

Right-hand side formulation for secondary potential field approach.

rhs_form_total()

Right-hand side formulation for total field approach. Right-hand side vector will be filled later by the TotalFieldAssembler instance.

class custEM.fem.frequency_domain_approaches.F_U_nodal(FS, MP)

Bases: ApproachBaseFD

FE implementation of F-U (Tau-Omega) approach on nodal elements. This approach was not implemented yet as it requires a homogeneous conductivity distribution (at least in our derivation of the equations) in the complete computational domain. Hence, it is irrelevant for geo-EM applications.

class custEM.fem.frequency_domain_approaches.H_vector(FS, MP)

Bases: ApproachBaseFD

FE implementation of H-fild approach (Rochlitz, 2019).

build_var_form(add_primary=True, **fem_kwargs)

Print FE parameters and initialize variational formulations for assembly. Valid keyword arguments are described in the ApproachBaseFD and PrimaryFields classes.

lhs_form_full()

Left-hand side contributions for full (not quasi-static) H-field approach, in development, not fully tested yet.

lhs_forms(system_it)

Left-hand side formulation as the basis for the system matrix assembly.

rhs_form_secondary(fi=0)

Right-hand side formulation for secondary magentic field approach.

rhs_form_total()

Right-hand side formulation for total field approach. Right-hand side vector will be filled later by the TotalFieldAssembler instance.

custEM.fem.primary_fields module

@author: Rochlitz.R

class custEM.fem.primary_fields.PrimaryField(FS, FE, MP, **pf_kwargs)

Bases: ApproachBaseFD

PrimaryField class for computing or loading existing primary fields for secondary field formulations. Called by FE.build_var_form() via the FunctionSpace class method add_Primary_Field().

For the documentation of the the *comet toolbox (pyhed module), it is referred to the documentation of this software.

Methods

• export_primary_fields()

  save calcualted primary fields

• load_primary_fields()

  import existing primary fields; if not existing, start calculation

• calc_hed_fullspace_field()

  calculate analytic primary-field solution for a HED in a fullspace

• calc_vmd_fullspace_field()

  calculate analytic primary-field solution for a VMD in a fullspace

• calc_layered_earth_field()

  calculate semi-analytic primary-field solutions for layered earth geometries with pyhed

• init_primary_field_name()

  generate an encrypted primary field data file name based on the Tx configuration to avoid recomputation of exisiting indentical setups (same mesh, same Tx geometry, same physical parameters)

• write_serial_calculation_parameters()

  store PF class parameters in a file for the reimport of serial calculation scripts

Description of PF parameters

Note that further relevant parameters are initialized by the ApproachBaseFD class. They are also valid keyword arguments of the build_var_form() method and can be changed accordingly.

• max_length = None, type float

  maximum length of the single HED segments used to approximate a finite-length bipole or loop transmitter; overwrites n_segs if both parameters are specified

• n_segs = None, type int

  total number of segments (single HED transmitters) for approximating finite-length bipole or loop transmitters; overwritten by max_length if both parameters are specified

• pf_type = ‘default’, type str

  primary field type, alternatives: fullspace or custom

• pf_name = None, type str

  an alternative name for the primary fields of the model can be defined; it is not recommended to change the default name

• export_pvd = False, type bool

  flag controlling if ‘pvd’ - files are exported (‘True’) or not

• eval_tx = False, type bool

  evaulate Tx coordinates on corresponding edge-midpoints between nodes and use them for the primary field calculation instead of the node coordinates, also provides edge lengths and orientations

• force_primary = False, type bool

  force re-caculation of exisiting primary field with an identical parametrization

• pyhed_interp = False, type bool

  set True, if the alternative interpolation-based primary field calculation technique of the pyhed-module should be used

• pf_EH_flag = ‘E’, type str

  either E or H to specify, which type of primary fields should be used for the FE caluclations later on

• procs_per_proc = 1, type int

  change, if more processes should be used to calculate primary fields with ‘pyhed’ than specified via the mpirun command (e.g., -n 12) for sloving the main FEm problem. For instance, a value of 3 would lead to 36 processes used for the primary field calculation.

• pf_p = 1, type int

  polynomials order of Lagrange VectorFunctionSpace which is used to calculate the primary fields; could be changes to 2, but this only leads to significantly increased primary field computation times and barely more accurate results

• ph_log_level = 0, type int

  log level of the pyhed software

• pyhed_drop_tol = 1., type bool

  internal drop tolerance for approximating the semi-analytical results at positions with singular fields values

• dipole_z = 0., type float

  shift z-value of source to the subsurface to enable calculations of primary fields for marine setups; use with caution as this option was not tested sufficiently

• Assembler = None, type class

  Assembler class that might be passed to the PF class for custom purposes

calc_hed_fullspace_field(V_cg)

Calculate analytic primary field solution for a HED in a fullspace.

Required arguments

• V_cg, type dolfin VectorFunctionSpace

  dolfin Lagrange vector function space

calc_layered_earth_field(V_cg, boundary_fields, ti)

Calculate semi-analytic primary field solution with the pyhed software for layered-earth geometries.

Required arguments

• V_cg, type dolfin VectorFunctionSpace

  dolfin Lagrange vector function space

• boundary_fields, type bool

  only for internal usage, used to calculate primary fields only on the domain boundaries to apply inhomogeneous Dirichlet BC

• ti, type int

  number of the Tx if multiple Tx are specified

calc_mt_fullspace_field(V_cg)

Calculate analytic primary field solution for a HED in a fullspace.

Required arguments

• V_cg, type dolfin VectorFunctionSpace

  dolfin Lagrange vector function space

calc_vmd_fullspace_field(V_cg)

Calculate analytic primary field solution for a VMD in a fullspace.

Required arguments

• V_cg, type dolfin VectorFunctionSpace

  dolfin Lagrange vector function space

export_primary_fields()

Saves calculated primary fields in dedicated HDF5 or .pvd files.

init_primary_field_name(boundary_fields)

Initialize an encoded, unique name for the files to store primary fields based on the mesh, transmitter characteristics and others.

Required arguments

• boundary_fields = False, type bool

  only for internal usage, used to calculate primary fields only on the domain boundaries to apply inhomogeneous Dirichlet BC

load_primary_fields(boundary_fields=False, fi=0)

Import existing or start calculation of primary fields.

Keyword arguments

• boundary_fields = False, type bool

  only for internal usage, used to calculate primary fields only on the domain boundaries to apply inhomogeneous Dirichlet BC

• fi = 0, type int

  iteration index over frequencies

write_serial_calculation_parameters()

Export geometry and transmitter setup for seperate primary field calculations in serial.

Required arguments

• path, type str

  current working directory, switch to this directory after the export of the parameter file

custEM.fem.time_domain_approaches module

@author: Rochlitz.R

class custEM.fem.time_domain_approaches.FourierTransformBased(FS, MP)

Bases: ApproachBaseTD

FE implementation of inverse Fourier-transform based approach. It calculates several tens of frequency-domain solutions with the total E-field approach and converts them to time domain afterwards.

Methods

• build_var_form()

  set up attributes required in the frequency-domain variational formulations

• initialize frequencies()

  initialize frequencies for the frequency-domain computations automatically if not given explicitly as keyword argument

• calc_frequency_solutions()

  calculate frequency-domain solutions with total E-field approach

• transform_to_time_domain()

  transform frequency-domain solutions to time domain

• fht_hahnstein_80()

  use default set of 80 filter coefficients to perform the fast Hankel transformation

RA approach specific parameters

• n_mult = 6, type int

  multiplyer of poles, do not change

• n_poles = 2, type int

  number of poles to be used

• pole_shift = None, type float

  shift poles in time; the default covered range is about 4 to 5 decades from the earliest times on; this values is adjusted automatically and it is not recommended to change it

• dc_field = None, type dolfin Function

  the DC field is automatically calculated on the fly; for special purposes, it might be computed externally and passed here

build_var_form(**fem_kwargs)

Initialize parameters for the Fourier-transform based approach and print FE parameters. A description of the keyword arguments is given in the ApproachBaseTD and **FourierTransformBased**classes.

calc_frequency_solutions(auto_interpolate, repeat_frequencies, **solve_kwargs)

Calculate frequency-domain solutions. Keyword arguments are piped to this method by the build_var_form() method.

fht_hanstein_80(fd_data, modus, dc_level=None)

Default fast Hankel transformation routine with 80 filter coefficients.

Required arguments

• fd_data, type numpy array

  frequency-domain data, only one component

• modus, type int

  modus switch to calculate either shut-on, shut-off or impulse response results

Keyword arguments

dc_level = None, type float

used to apply explicitly a calculated DC level for shut-off E-field transients calculations; otherwise, the DC level will be approximated by the fields values of the lowest frequency

initalize_frequencies()

Initialize frequencies for the frequency-domain solutions automatically based on the default 80 FHT filter coefficients.

transform_to_time_domain(interp_mesh, EH, shut_on=True, shut_off=True, impulse=True)

Transform frequency-domain solutions to time domain using the fast Hankel transformation routine.

Required arguments

• interp_mesh, type str

  interpolation mesh for the chosen frequency-domain results which should be converted to time domain

• EH, type str

  default is ‘EH’, that means both, electric and magentic fields are converted. Users can choose only ‘E’ or ‘H’.

Keyword arguments

• shut_on = True, type str

  transform to shut-on field transients or not if False

• shut_off = True, type str

  transform to shut-off field transients or not if False

• impulse = True, type str

  transform to impulse (dB/dt) repsonse transients or not if False

class custEM.fem.time_domain_approaches.ImplicitEuler(FS, MP)

Bases: ApproachBaseTD

FE implementation of Implicit Euler time-stepping method, e.g., (Um, 2010). Users can use either a first- or second-order discretization in time and, with cautioan, a formulation considering displacement currents which needs some more tests.

Methods

• time_stepping()

  runs the complete time-stepping solution process as well as H-field conversion and export methods, FE modeling and time stepping related keyword arguments can be passed to this method.

• build_var_form()

  set up variational formulation and assemble mass and curl-curl matrices as well right-hand side (source) vector

• init_initial_condition()

  initialize initial condition for time-stepping approach

• first_order_IE()

  perform first-order implicit Euler time-stepping approach

• second_order_IE()

  perform second-order implicit Euler time-stepping approach

• first_order_IE_full()

  perform first-order implicit Euler time-stepping approach including displacement currents

• init_times()

  initalize internally used time vector for the time stepping

• gaussian()

  define gaussian source pulse for impulse response modeling

• print_stepping_info()

  print some information before starting the time-stepping process

IE approach specific parameters

• source_type = ‘djdt’, type str

  alternatively, use ‘j’ to calculate B-fields instead of dB/dt

• source_func = None, type numpy array

  specify custom source (transmitter) function over time; needs to correspond to suitable source_times

• source_times = None, type numpy array

  specify custom source_times if no default ramp-free source functions are used

• stepping_times = None, type numpy array

  specify custom calculation times for the time-stepping after the source pulse

• be_order = 1, type int

  order of discretization in time for backward Euler method, alterntively set to 2

• n_lin_steps = 30, type int

  number of linear steps between each logarithmic time step, each linear step requires only an inexpensive back-substitution for solving the system of linear equations

• n_on = 101, type int

  number of linear time steps for creating source signal, not recommended to be changed

• t0_j = -10., type float

  specify start time (negative) for specifying the duration of the source signal ramp if the source type is j

• pulse_width = 1e-2, type float

  relative width of dirac function, not recommended to be changed

• t_tol = 1e-2, type float

  tolerance to check if a calculated time step is linear or not (if not linear, a new matrix factorization will be computed)

build_var_form(**fem_kwargs)

Build variational formulation for E-field time-stepping approach and print FE parameters. A description of the keyword arguments is given in the ApproachBaseTD and **ImplicitEuler**classes.

first_order_IE(u_0)

First-order backward Euler time-stepping scheme using diffusive approximation (no displacement currents).

Required arguments

• u_0, type dolfin function

  initial solution at first time step

first_order_IE_full(u_0, u_1)

First-order backward Euler time-stepping scheme for full formulation considering displacement currents.

Required arguments

• u_0, type dolfin function

  initial solution at first time step

• u_1, type dolfin function

  intital solution at second time step

gaussian(x, x0, alpha)

Calculate gaussian pulse.

Required arguments

• x, type numpy array

  time values

• x0, type float

  central value of times values

• alpha, type float

  stretching parameter of Dirac pulse function

init_initial_condition(ini_con, order=1)

Initialize the initial fields u_0 and u1 (if second order scheme or full formulation) for the time-stepping.

Required arguments

• ini_con, type dolfin function

  use customized initial fields, not supported yet

Keyword arguments

• order = 1, type int

  specify first or second order time-stepping scheme

init_times()

Initialize times for time-stepping. Depending on the source_type, initial time steps before are added automatically to incorporate the source signal before the observation times are calculated.

print_stepping_info(dt)

Print some information before starting the time-stepping process.

Required arguments

• dt, type array

  vector of time-steps

second_order_IE(u_0, u_1)

Second-order backward Euler time-stepping scheme.

Required arguments

• u_0, type dolfin function

  initial solution at first time step

• u_1, type dolfin function

  intital solution at second time step

time_stepping()

Initialize time-discfretization and variational formulations for assembly.

class custEM.fem.time_domain_approaches.RationalArnoldi(FS, MP)

Bases: ApproachBaseTD

FE implementation of rational Arnoldi (Krylov-subspace) method after Börner (2015).

Methods

• rational_arnoldi()

  runs the complete rational-arnoldi solution process as well as H-field conversion and export methods, FE modeling and time stepping related keyword arguments can be passed to this method

• build_var_form()

  set up variational formulation and assemble mass and curl-curl matrices as well as right-hand side (source) vector

• get_poles_from_table()

  define poles used in rational Arnoldi method

• build_rational_arnoldi_basis()

  build Krylov subspace basis

• expm_ra()

  utility function to calculate matrix exponentials

• project_solution()

  calculate solution on Krylov subspace and project it back

RA approach specific parameters

• n_mult = 6, type int

  multiplyer of poles, do not change

• n_poles = 2, type int

  number of poles to be used

• pole_shift = None, type float

  shift poles in time; the default covered range is about 4 to 5 decades from the earliest times on; this values is adjusted automatically and it is not recommended to change it

• dc_field = None, type dolfin Function

  the DC field is automatically calculated on the fly; for special purposes, it might be computed externally and passed here

build_krylov_subspace_basis()

Build the Krylov subspace basis.

build_var_form(**fem_kwargs)

Build variational formulation for rational Arnoldi approach and print FE parameters. A description of the keyword arguments is given in the ApproachBaseTD and RationalArnoldi classes.

expm_ra(A, E)

Calculate matrix exponentials with eigenvalues and eigenvectors.

Required arguments

• A, E, type numpy ndarray

  input matrices for the calculation

get_poles_from_table()

Define poles for ration Arnoldi method.

project_solution()

Project subspace solution back to original solution space.

rational_arnoldi()

Conduct rational Arnoldi modeling.

Module contents

fem

Submodules:

• ** frequency_domain_approaches**: E-field, H-field, A-V-mixed, A-V-nodal, F-U-mixed

• ** time_domain_approaches**: Implicit Euler, Inverse Fourier-Transform based, Rational Arnoldi

• fem_base for setting up the finite element kernel

• fem_utils for defining approach base class and related utility functions

• primary_fields for setting up primary fields

• assembly of left-/right-hand-sides using total field approaches

---