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 righthandsides (two source polarizations implemented as boundary conditions on the top
of the model domain) for naturalsource modeling. assemble()
 assembly function of left and righthand sides for secondaryfield formulations

class
custEM.fem.assembly.
SecondaryFieldAssembler
(FE, bc)[source]¶ Bases:
object
Class that assembles the system matrix and the righthandsides for secondaryfield formulations based on primary fields and delta_sigma.
 assemble()
 assembly function of left and righthand sides for secondaryfield formulations

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 righthand side vector for total field formulations.
 assemble()
 general assembly function of left and righthand sides for totalfield 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 AV formulation on nodal elements
 add_topo(),
 modify zvalues of coordinates according to specified topography

add_topo
(tx, tx_topo=None)[source]¶ Apply topography information to the zvalues 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 dofcoordinate 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 righthand side assembly of source currents manually and let dolfin assemble the system matrix.
 fi = 0, type int
 iteration index over frequencies for specifying the lefthandside

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, counterclockwise in a righthanded coordinate system

find_and_fill_start_stop_dofs
(bb, current, ti=0, fi=0)[source]¶ Required for the AV 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 righthand 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.
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 pyhedinternal multiprocessing from a serial process to ensure correct dof ordering.
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 pyhedinternal multiprocessing from a serial process.
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 ‘FEniCSway’, e.g. predefined anomaly instances from the file misc/anomaly_expressions.py or userdefined dolfin expressions can be added to the model, each anomaly gets a new marker

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 elementwise scalar parameter function (isotropic parameters)
 build_dg_tensor_func()
 build elementwise tensor parameter function (anisotropic parameters)
 build_complex_sigma_func()
 build complexvalued sigma functions (only isotropic) for modeling induced polarization effects
 import mesh (see init_mesh documentation)
 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 = elementwise scalar discontinuous FunctionSpace
 order 0, 1 value per cell
 S_dgt = elementwise tensorvalued 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)
 initalize empty solution Function U on mixed solution space
 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 AV approach, which requires Efield 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 Hfield conversion, provides no benefits so far but might be reconsidered 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 complexvalued conductivities for simulating induced polarization effects with a ColeCole 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))

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, singlevalued functions for isotropic values and tensorvalued functions for anisotropic values
 print_model_parameters()
 print updated list of model parameters
 load_mesh_parameteres()
 update meshrelated model parameters from ‘meshparameters’ 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 timedomain and for DC modeling, the first dimension of the current data array corresponds to n_tx; in frequencydomain 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 (righthand sides of the system of linear equations), calculated automatically
 sigma_air = 1e8, type float or list
 homogeneous airspace conductivity, for modeling a fullspace, set sigma_air to the same value as sigma_ground
 sigma_ground = 1e2, 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 secondaryfield/potential formulations; arbitrary general anisotropic values are supported for sigma_ground and sigma_0; delta_sigma is evaluated internally
 sigma = [1e2, 1e8], 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 anomalydomains 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 threelayered 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 * 1e7, 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 * 1e12, 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 ColeCole model for modeling induced polarization effects
 ip_c = 1., type float or list of floats
 c values of ColeCole model for modeling induced polarization effects
 ip_m = 0., type float or list of floats
 m values of ColeCole model for modeling induced polarization effects
 n_layers = 1, type int
 number of subsurface layers, automatically imported from meshparameter file
 layer_depths = [], type list of float
 depths of subsurfacelayer interfaces, automatically imported from meshparameter 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 meshparameter 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.
custEM.fem.fem_utils module¶
@author: Rochlitz.R

class
custEM.fem.fem_utils.
ApproachBaseFD
[source]¶ Bases:
object
Initialize FE variables for all frequencydomain 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 meshparameter 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 meshparameter file; otherwise, the geometry parameters will be used to define the tx
 n_tx = 1, type int
number of transmitters (righthand 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 meshparameter 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 righthand 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 ColeCole 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 = 1e2, 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

class
custEM.fem.fem_utils.
ApproachBaseTD
[source]¶ Bases:
object
Initialize FE variables for all timedomain 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 meshparameter file in the MP instance
 calc_static_magnetic_field()
 calculate static magnetic field with the BiotSavart law
 biot_savart()
 vectorized implementation of BiotSavart law
 interpolate_static()
 interpolate static magentic field on specified interpolation mesh for independent postprocessing 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 postprocessing 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 timedomain 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 timedomain 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 shutoff of shuton response

biot_savart
(rx, tx)[source]¶ Vectorized implementation of the BiotSavart 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 shutoff 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 BiotSavart 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
 nondistributed 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

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 lefthand side of variational formulation
 rhs_form_total()
 build righthand side of variational formulation for secondary potential approach
 rhs_form_secondary()
 build righthand 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 frequencydomain 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 righthand 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

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.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.


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).

class
custEM.fem.frequency_domain_approaches.
E_vector
(FS, MP)[source]¶ Bases:
custEM.fem.fem_utils.ApproachBaseFD
FE implementation of Efield 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_full
()[source]¶ Lefthand 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.


class
custEM.fem.frequency_domain_approaches.
F_U_mixed
(FS, MP)[source]¶ Bases:
custEM.fem.fem_utils.ApproachBaseFD
FE implementation of TauOmega approach on mixed elements (Mitsuhata, 2004), called FU 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.


class
custEM.fem.frequency_domain_approaches.
F_U_nodal
(FS, MP)[source]¶ Bases:
custEM.fem.fem_utils.ApproachBaseFD
FE implementation of FU (TauOmega) 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 geoEM applications.

class
custEM.fem.frequency_domain_approaches.
H_vector
(FS, MP)[source]¶ Bases:
custEM.fem.fem_utils.ApproachBaseFD
FE implementation of Hfild 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]¶ Lefthand side contributions for full (not quasistatic) Hfield approach, in development, not fully tested yet.

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 primaryfield solution for a HED in a fullspace
 calc_vmd_fullspace_field()
 calculate analytic primaryfield solution for a VMD in a fullspace
 calc_layered_earth_field()
 calculate semianalytic primaryfield 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 finitelength 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 finitelength 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 edgemidpoints 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 recaculation of exisiting primary field with an identical parametrization
 pyhed_interp = False, type bool
 set True, if the alternative interpolationbased primary field calculation technique of the pyhedmodule 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 semianalytical results at positions with singular fields values
 dipole_z = 0., type float
 shift zvalue 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 semianalytic primary field solution with the pyhed software for layeredearth 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

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
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 Fouriertransform based approach. It calculates several tens of frequencydomain solutions with the total Efield approach and converts them to time domain afterwards.
 build_var_form()
 set up attributes required in the frequencydomain variational formulations
 initialize frequencies()
 initialize frequencies for the frequencydomain computations automatically if not given explicitly as keyword argument
 calc_frequency_solutions()
 calculate frequencydomain solutions with total Efield approach
 transform_to_time_domain()
 transform frequencydomain 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 Fouriertransform 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 frequencydomain 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
 frequencydomain data, only one component
 modus, type int
 modus switch to calculate either shuton, shutoff or impulse response results
 dc_level = None, type float
 used to apply explicitly a calculated DC level for shutoff Efield transients calculations; otherwise, the DC level will be approximated by the fields values of the lowest frequency

initalize_frequencies
()[source]¶ Initialize frequencies for the frequencydomain 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 frequencydomain solutions to time domain using the fast Hankel transformation routine.
 interp_mesh, type str
 interpolation mesh for the chosen frequencydomain 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 shuton field transients or not if False
 shut_off = True, type str
 transform to shutoff 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 timestepping method, e.g., (Um, 2010). Users can use either a first or secondorder discretization in time and, with cautioan, a formulation considering displacement currents which needs some more tests.
 time_stepping()
 runs the complete timestepping solution process as well as Hfield 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 curlcurl matrices as well righthand side (source) vector
 init_initial_condition()
 initialize initial condition for timestepping approach
 first_order_IE()
 perform firstorder implicit Euler timestepping approach
 second_order_IE()
 perform secondorder implicit Euler timestepping approach
 first_order_IE_full()
 perform firstorder implicit Euler timestepping 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 timestepping process
 source_type = ‘djdt’, type str
 alternatively, use ‘j’ to calculate Bfields 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 rampfree source functions are used
 stepping_times = None, type numpy array
 specify custom calculation times for the timestepping 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 backsubstitution 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 = 1e2, type float
 relative width of dirac function, not recommended to be changed
 t_tol = 1e2, 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 Efield timestepping 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]¶ Firstorder backward Euler timestepping 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]¶ Firstorder backward Euler timestepping 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 timestepping.
 ini_con, type dolfin function
 use customized initial fields, not supported yet
 order = 1, type int
 specify first or second order timestepping scheme

init_times
()[source]¶ Initialize times for timestepping. 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 timestepping process.
 dt, type array
 vector of timesteps

class
custEM.fem.time_domain_approaches.
RationalArnoldi
(FS, MP)[source]¶ Bases:
custEM.fem.fem_utils.ApproachBaseTD
FE implementation of rational Arnoldi (Krylovsubspace) method after Börner (2015).
 rational_arnoldi()
 runs the complete rationalarnoldi solution process as well as Hfield 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 curlcurl matrices as well as righthand 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_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.
Module contents¶
fem¶
Submodules:
 ** frequency_domain_approaches**: Efield, Hfield, AVmixed, AVnodal, FUmixed
 ** time_domain_approaches**: Implicit Euler, Inverse FourierTransform 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/righthandsides using total field approaches