custEM.fem package

Submodules

custEM.fem.assembly module

@author: Rochlitz.R

class custEM.fem.assembly.MTAssembler(FE, bc)[source]

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.
  • assemble()
    assembly function of left- and right-hand sides for secondary-field formulations
assemble(fi=0)[source]
class custEM.fem.assembly.SecondaryFieldAssembler(FE, bc)[source]

Bases: object

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

  • assemble()
    assembly function of left- and right-hand sides for secondary-field formulations
assemble(fi=0)[source]

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

  • 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)[source]

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.

  • 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 dofs 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_dofs()
    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)[source]

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

  • tx, type list of coordinates
    list of transmitter coordinates
  • 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)[source]

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

  • idx, type list of int
    integers specifying the transmitter dof
  • tdir, type float
    either positive or negative 1
  • comp = 0, type int
    either 0 or 1 for x or y component
assemble(fi=0)[source]

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

  • fi = 0, type int
    iteration index over frequencies for specifying the left-hand-side
eval_target_dir(a, b)[source]

Evaluate global direction of each transmitter segment.

  • 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_dofs(bb, current, ti=0, fi=0)[source]

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.

  • bb, dolfin function
    assembled right-hand side function
  • current, type float,
    source current for the corresponsing tx
  • ti = 0, type int
    index of the corresponding tx
find_cells()[source]

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

find_directions()[source]

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()[source]

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

find_for_hed()[source]

Find source dof and coordinates for HED source.

find_for_line()[source]

Find source dof and coordinates for “line” source.

find_for_loop_r()[source]

Find source dof and coordinates for “loop_r” source.

find_for_path()[source]

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

find_for_vmd()[source]

Find source dof and coordinates for VMD source.

test_An_total(fi)[source]

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)[source]

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()[source]
convertCoordinates(gimli, dolfin)[source]

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)[source]
read_serial_calculation_parameters()[source]
write_h5(var, fname)[source]

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)[source]

Bases: object

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

calc_fields()[source]
convertCoordinates(gimli, dolfin)[source]

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)[source]
read_serial_calculation_parameters()[source]
write_h5(var, fname)[source]

custEM.fem.fem_base module

@author: Rochlitz.R

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

Bases: object

Domain class called internally from FunctionSpace instance.

  • 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)[source]

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)[source]

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.

  • 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
  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)[source]

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

  • 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_primary_fields(FE, MP, calc=True, **kwargs)[source]

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

  • FE, type class
    FE approach instance (e.g., E_vector)
  • MP, type class
    ModelParameters instance
  • 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)[source]

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.

  • 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)[source]

Assign isotropic conductivity values to discontinous function spaces.

  • vals, type list (of lists)
    conductivity values specified in MP instance
  • inverse = False, type bool
    return either sigma or 1/sigma (=inv(sigma))
build_dg_tensor_func(vals, inverse=False, logger=None)[source]

Assign anisotropic conductivity values to discontinous tensor function spaces.

  • vals, type list (of lists)
    conductivity values specified in MP instance
  • inverse = False, type bool
    return either sigma or inv(sigma)
init_mesh(test_mode, MP, self_mode, overwrite_mesh)[source]

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

  • see documentation of MOD class
class custEM.fem.fem_base.ModelParameters(logger, str_args, test_mode)[source]

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.

  • 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
  • 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
  • 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)[source]

Build discontinuous isotropic or anisotropic parameter functions.

calc_delta_sigma()[source]

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

check_parameter_dx_conformity()[source]

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

load_mesh_parameters()[source]

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

print_model_parameters()[source]

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

update_model_parameters(**mod_kwargs)[source]

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[source]

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.

  • 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
  • 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)[source]

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

fem_kwargs, type dict
dictionary of keyword arguments for build_var_form() method
update_tx_parameters()[source]

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

class custEM.fem.fem_utils.ApproachBaseTD[source]

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.

  • 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

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)[source]

Vectorized implementation of the Biot-Savart law.

  • 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)[source]

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

  • 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)[source]

Calculate static magentic fields with the Biot-Savart law.

  • store = False, type bool
    save static magnetic fields as HDF5 file
export_interpolation_data(field, V_serial, ti, export_pvd, interp_mesh)[source]

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

  • 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)[source]

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

  • interp_mesh, type str
    name of target interpolation mesh
  • 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)[source]

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

  • interp_mesh, type str
    name of target interpolation mesh
  • 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)[source]

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

  • fem_kwargs, type dict
    dictionary of keyword arguments for build_var_form() method
update_tx_parameters()[source]

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

class custEM.fem.fem_utils.DC(FS, MP)[source]

Bases: custEM.fem.fem_utils.ApproachBaseFD

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

  • 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

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)[source]

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)[source]

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.

  • 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()[source]

Calculate primary potential of static current analytically.

lhs_form()[source]

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

rhs_form_secondary()[source]

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

rhs_form_total()[source]

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

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

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

  • vals, type list (of lists)
    conductivity values specified in MP instance
  • inverse = False, type bool
    return either sigma or 1/sigma (inv(sigma))
custEM.fem.fem_utils.check_sigma_vals(MP)[source]

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

  • MP, type class
    ModelParameters instance
custEM.fem.fem_utils.check_source(s_type, tx, pf_flag='E', approach=None)[source]

Evaulate, which source type is chosen.

  • s_type, type int
    source type (see ApproachBaseFD documentation)
  • tx, type list of arrays
    list of transmitter coordinates for each Tx
  • 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)[source]

Bases: custEM.fem.fem_utils.ApproachBaseFD

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

build_var_form(**fem_kwargs)[source]

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

lhs_forms()[source]

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

rhs_form_secondary(fi=0)[source]

Right-hand side formulation for secondary potential field approach.

rhs_form_total()[source]

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)[source]

Bases: custEM.fem.fem_utils.ApproachBaseFD

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

build_var_form(**fem_kwargs)[source]

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

lhs_form()[source]

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

rhs_form_secondary(fi=0)[source]

Right-hand side formulation for secondary field approach.

rhs_form_total()[source]

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)[source]

Bases: custEM.fem.fem_utils.ApproachBaseFD

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

build_var_form(**fem_kwargs)[source]

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

lhs_forms()[source]

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

lhs_forms_full()[source]

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)[source]

Right-hand side formulation for secondary electric field approach.

rhs_form_total()[source]

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)[source]

Bases: custEM.fem.fem_utils.ApproachBaseFD

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

build_var_form(**fem_kwargs)[source]

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

lhs_forms()[source]

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

rhs_form_secondary(fi=0)[source]

Right-hand side formulation for secondary potential field approach.

rhs_form_total()[source]

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)[source]

Bases: custEM.fem.fem_utils.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)[source]

Bases: custEM.fem.fem_utils.ApproachBaseFD

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

build_var_form(add_primary=True, **fem_kwargs)[source]

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

lhs_form_full()[source]

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

lhs_forms()[source]

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

rhs_form_secondary(fi=0)[source]

Right-hand side formulation for secondary magentic field approach.

rhs_form_total()[source]

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)[source]

Bases: custEM.fem.fem_utils.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.

  • 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

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)[source]

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

  • V_cg, type dolfin VectorFunctionSpace
    dolfin Lagrange vector function space
calc_layered_earth_field(V_cg, boundary_fields, ti)[source]

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

  • 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)[source]

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

  • V_cg, type dolfin VectorFunctionSpace
    dolfin Lagrange vector function space
calc_vmd_fullspace_field(V_cg)[source]

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

  • V_cg, type dolfin VectorFunctionSpace
    dolfin Lagrange vector function space
export_primary_fields()[source]

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

init_primary_field_name(boundary_fields)[source]

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

  • 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)[source]

Import existing or start calculation of primary fields.

  • 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()[source]

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

  • 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)[source]

Bases: custEM.fem.fem_utils.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.

  • 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
  • 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)[source]

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)[source]

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)[source]

Default fast Hankel transformation routine with 80 filter coefficients.

  • 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
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()[source]

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)[source]

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

  • 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’.
  • 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)[source]

Bases: custEM.fem.fem_utils.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.

  • 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
  • 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)[source]

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)[source]

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

  • u_0, type dolfin function
    initial solution at first time step
first_order_IE_full(u_0, u_1)[source]

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

  • 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)[source]

Calculate gaussian pulse.

  • 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)[source]

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

  • ini_con, type dolfin function
    use customized initial fields, not supported yet
  • order = 1, type int
    specify first or second order time-stepping scheme
init_times()[source]

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)[source]

Print some information before starting the time-stepping process.

  • dt, type array
    vector of time-steps
second_order_IE(u_0, u_1)[source]

Second-order backward Euler time-stepping scheme.

  • u_0, type dolfin function
    initial solution at first time step
  • u_1, type dolfin function
    intital solution at second time step
time_stepping()[source]

Initialize time-discfretization and variational formulations for assembly.

class custEM.fem.time_domain_approaches.RationalArnoldi(FS, MP)[source]

Bases: custEM.fem.fem_utils.ApproachBaseTD

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

  • 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
  • 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()[source]

Build the Krylov subspace basis.

build_var_form(**fem_kwargs)[source]

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)[source]

Calculate matrix exponentials with eigenvalues and eigenvectors.

  • A, E, type numpy ndarray
    input matrices for the calculation
get_poles_from_table()[source]

Define poles for ration Arnoldi method.

project_solution()[source]

Project subspace solution back to original solution space.

rational_arnoldi()[source]

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