# -*- coding: utf-8 -*-
"""
@author: Rochlitz.R
"""
import pygimli as pg
pg.core.setThreadCount(1)
import dolfin as df
import os
import custEM as ce
from custEM.misc import logger_print as lp
from custEM.misc import export_config_file
from custEM.misc import export_resource_file
from custEM.misc import check_approach_and_file_format
from custEM.misc import get_logger
from custEM.misc import max_mem
import numpy as np
import sys
import time
import logging
[docs]
class MOD:
"""
Main class for 3D CSEM modeling using FEniCS. After an instance of this
mainclass is created, all modeling steps are conducted within this
environement. Instances within the mainclass 'MOD' are:
- FS: :class:`custEM.fem.fem_base.FunctionSpaces`
- FS.PF: :class:`custEM.fem.primary_fields.PrimaryField`
- FS.DOM: :class:`custEM.fem.fem_base.Domains`
- MP: :class:`custEM.fem.fem_base.ModelParameters`
- Solver: :class:`custEM.core.solvers.Solver`
- FE - depending on the time- or frequency-domain modeling approach:
* :class:`custEM.fem.fem_utils.ApproachBaseTD`
* :class:`custEM.fem.time_domain_approaches.ImplicitEuler`
* :class:`custEM.fem.time_domain_approaches.FourierTransformBased.`
* :class:`custEM.fem.time_domain_approaches.RationalArnoldi`
* :class:`custEM.fem.fem_utils.ApproachBaseFD`
* :class:`custEM.fem.frequency_domain_approaches.E_vector`
* :class:`custEM.fem.frequency_domain_approaches.H_vector`
* :class:`custEM.fem.frequency_domain_approaches.A_V_nodal`
* :class:`custEM.fem.frequency_domain_approaches.A_V_mixed`
* :class:`custEM.fem.frequency_domain_approaches.F_U_mixed`
* :class:`custEM.fem.fem_utils.DC`
- IB - depending on time- or frequency-domain modeling:
* :class:`custEM.post.interp_tools_td.TimeDomainInterpolator`
* :class:`custEM.post.interp_tools_fd.FrequencyDomainInterpolator`
- PP - depending on post-processing or pre-processing:
* :class:`custEM.core.post_proc.PostProcessing`
* :class:`custEM.core.pre_proc.PreProcessing`
Methods
-------
- solve_main_problem()
function to solve the weak formulation of the main system,
defined in the **FE** instance; note that iterative solvers are
theoretically considered but not working properly yet
- init_third_party_parameters()
adjust a few dolfin and numpy parameters; so far, no need for
customization
- init_default_MOD_parameters()
initialize default parameters which can be overwritten by
specifying keyword arguments when calling the **MOD** instance
- init_instances()
initialize all required 'sub'-classes of the **MOD** class
- handle_exception()
overwrite excepthook from sys module for customized logger
"""
def __init__(self, mod_name, mesh_name, approach, debug_level=20,
fenics_debug_level=50, mute=False, **main_kwargs):
"""
Initializes Mesh, Functionspaces, ModelParameters, FE kernel, Solvers,
Interpolation base, Post- or Pre-Processing kernel. Note that the
*fenics_debug_level* influences only the debugging level of FEniCS.
In custEM, most prints can be suppressed with the *mute* argument,
but different
log levels were not implemented yet (planned in future versions).
Required arguments
------------------
- mod_name, type str
name of the model
- mesh_name, type str
name of input mesh without suffix
- approach, type str
* either **E_t**, **H_t**, **Am_t**, **An_t** or
**E_s**, **H_s**, **Am_s**, **Fm_s**, **An_s** for frequency-domain
modeling; Note: **H_t** and **An_t** are theoretically considered
but not working yet
* either **E_IE**, **E_FT**, or **E_RA** for time-domain modeling
* **DC** for direct current field simulations
* **MT** for magentotelluric modeling (in development)
Keyword arguments
-----------------
- debug_level = 20 = "INFO", type int or str
specify level of custEM debug logger, levels refer to the ones of
the standard Python library *logging*:
10 or "DEBUG", diagnose problems, show detailed information
20 or "INFO", confirm that things are working as expected
30 or "WARNING", indicate if anyhing unexpected happens
40 or "ERROR", show serious problems
50 or "CRITICAL", show serious errors which force to quit
- fenics_debug_level = 50, type int
specified log level in FEniCS, default is *50* to keep FEniCS as
quiet as possible and return only critical errors.
- mute = False, type bool
set **True**, if unimportant prints should be muted
- p = 1, type int
polynomial order, alternatively change to order **2** (or **3**)
- r_dir = 'results', type str
results directory, can be changed here or as general default
variable **your_out_dir** in the file *misc/paths.dat*
- m_dir = 'meshes', type str
mesh directory, same possibilites as for *r_dir*
- para_dir = m_dir + '/para', type str
directory of mesh parameter files
- out_dir = r_dir + '/' + approach, type str
redefine to **your_output_directory** for custom output directory
- profiler = True, type bool
automatically write logger information, model configuration and
resource demands to file
- load_exisiting = False, type bool
set **True** to load results form an existing model, readable
data are stored in the "PP" (preprocessing) instance
- overwrite_results = False, type bool
set **True** to overwrite results of a model with the same name
on the same mesh
- overwrite_mesh = False, type bool
set **True** to allow automated conversion of a never version of a
created mesh with the same name
- serial_ordering = False, type bool
use serial ordering of *METIS* instead parallel ordering of
*SCOTCH*; leads to slight increase of solution time for *MUMPS*
but significantly reduces chance of potential *MUMPS* crashes
- field_selection = 'all', type str
specify a number of selected quantities to speed up and avoid
unnecessary imports. The string can contain "E_t", "E_s", "H_t",
"H_s", (potentially "A_t" and "A_s") to specify the quantities
- test_mode = False, type bool
set **True** for test mode with a 6 * (20 x 20 x 20) cells
UnitCubeMesh (dimensions of 20 x 20 x 20 km). If **True**, a
specified **mesh** is always ignored
- import_freq = None, type int
integer specifiying which frequency id should be imported;
if *None*, a preprocessing instance can be initialized in parallel
without importing any fileds to speed up the intializtion process;
if *n_freqs* for this model is 1, the value is
irrelevant as there is only one frequency
- file_format = 'h5', type str
h5 for parallel I/O, alternatively change value
to **'xml'** for dolfin-xml format
- fs_type = 'None', type str
set **Ned/ned** or **CG/cg** to restrict import of functions
on either Nedelec or CG Funcionspaces (if they exist)
- dg_interpolation = False, type bool
set **True** if fields should be interpolation from Nedelec-Spaces
on discontinuous Lagrange-VectorFunctionSpaces during
post-processing/interpolation; use with caution!, this is usually
not necessary
- export domains = True, type bool
set **False** to not export the mesn with domain information
- mumps_debug = False, type bool
set **True** to enable debugging of MUMPS solver by adjusting
the MUMPS log level
- reuse_fs = None, type bool
reuse the function spaces of an existing model (with the same
mesh), currently used for internal usage of the DC class to copy
the dof information for calculating DC fields in time-domain
problems; note that this flag was required since FEniCS does NOT
distribute the same mesh identically in subsequent **MOD** calls
- mesh_only = False, type bool
set True to use the implemented routines for mesh and survey
parameter import only without initializing further instances,
mainly used as mesh-convert utility in the inv_utils
- self_mode = False, type bool
the complete simulation will be conducted in serial on the process
intitializing **MOD**, used mainly for the serial interpolation
for internal usage
- ned_kind = '1st', type str
used only for experiments, alternatives would be "2nd" or "div",
do not change!
- mpi_cw/mpi_cs, type MPI communicators
MPI comm. attributes for internal usage, FEniCS version dependent
- mpi_rank, type MPI communicators
MPI rank of procs for internal usage, FEniCS version dependent
"""
self.mute = mute
self.mod_name = mod_name
self.mesh_name = mesh_name
self.approach = approach
self.init_default_model_parameters()
# not the nicest implementation for initialization, but the logger
# already requires the data directory structure to be able to write
# the log file
if 'r_dir' in main_kwargs:
self.r_dir = main_kwargs['r_dir']
if 'm_dir' in main_kwargs:
self.m_dir = main_kwargs['m_dir']
if 'profiler' in main_kwargs:
self.profiler = main_kwargs['profiler']
self.out_dir = (self.r_dir + '/' + self.approach +
'/' + self.mesh_name)
ce.misc.make_directories(self.r_dir, self.m_dir, self.approach)
if df.MPI.size(self.mpi_cw) > 1 and df.MPI.rank(self.mpi_cw) == 0:
if not os.path.isdir(self.out_dir):
os.makedirs(self.out_dir)
else:
if df.MPI.size(self.mpi_cw) == 1:
if not os.path.isdir(self.out_dir):
os.makedirs(self.out_dir)
else:
pass
self.logger = get_logger(debug_level,
self.profiler,
self.out_dir + '/' + mod_name)
sys.excepthook = self.handle_exception
if self.load_existing:
lp(self.logger, 20, 'Model class initialization:',
pre_dash=False, barrier=False)
else:
lp(self.logger, 20, 'Model class initialization:', barrier=False)
# support old syntax for a while
if 'overwrite' in main_kwargs:
lp(self.logger, 30,
'Warning! The use of overwrite is deprecated and will '
'be disabled\nin future custEM versions.\nUse the'
'*overwrite_results* and *overwrite_mesh* flags '
'instead.\nContinuing with *overwrite* as alias for '
'*overwrite_results* ...', post_dash=True)
self.__dict__['overwrite_results'] = \
main_kwargs.pop('overwrite')
for key in main_kwargs:
if key not in self.__dict__:
lp(self.logger, 50,
'Error! Unknown keyword argument for MOD set', key)
lp(self.logger, 50, "... let's stop "
"before something unexpected happens ...")
raise SystemExit
if 'solvingapproach' in main_kwargs:
self.solvingapproach = main_kwargs['solvingapproach']
if 'iterative' in main_kwargs:
self.iterative = main_kwargs['iterative']
if self.iterative == False:
self.solvingapproach=['direct','MUMPS','directl','H_direct']
elif self.iterative == True:
self.solvingapproach=['iterative','PRESB','iterative','H_iterative']
else:
lp(self.logger, 50,
'Error! Unknown keyword argument for MOD set', key)
lp(self.logger, 50, "... let's stop "
"before something unexpected happens ...")
raise SystemExit
self.__dict__.update(main_kwargs)
if 'para_dir' not in main_kwargs:
self.para_dir = self.m_dir + '/para'
self.init_third_party_parameters(fenics_debug_level)
if self.test_mode is not False:
self.mod_name = 'TEST_MODE_' + self.mod_name
lp(self.logger, 20, 'Notice: Prefix "TEST_MODE_" added to the '
'given *mod_name*', self.mod_name)
if not self.self_mode:
logger = self.logger
lp(self.logger, 20, '... initializing sub-module '
'instances ...', pre_dash=False)
else: # for nested serial interpolation
logger = int(self.logger.level)
if logger < 21:
print('... initializing sub-module instances ...')
ce.misc.check_if_model_exists(self.out_dir + '/' + self.mod_name,
self.overwrite_results,
self.load_existing, logger)
self.init_instances()
[docs]
def solve_main_problem(self, bc=None, sym=True,
out_of_core=False, convert=True,
method='fgmres', delete_factorization=True,
auto_interpolate=False, calculation_frequencies=[],
secondary_potential=False, inv_dir=None,
interpolation_quantities=['E_t', 'H_t'],
remote_station=None,
**post_proc_kwargs):
"""
Assemble and solve linear systems of equations. Some arguments are
piped to the **PostProcessingTD** or **PostProcessingFD** classes
with the *post_proc_kwargs* dictionary during the post-processing
phase, if automated postprocessing is enabled (*convert=True*).
For the documentation of valid keyword arguments for the
*PP.convert_results()* method, it is referred to the post-processing
classes.
Keyword arguments
-----------------
- solver = 'MUMPS', type int
use direct solver MUMPS, no alternatives so far
- bc = None, type str
default boundary conditions (BC) are implicit Neumann bc;
alternatively, **'ZeroDirichlet'**, short **'ZD'**, can be choosen
or inhomogeneous Dirichlet conditions **'ID'** in case of the
frequency-domain *E-field* approach
- sym = 1, 2 or True, type int or bool
so far, all implemented systems of equations are symmetric,
this parameter may be **False** in other FE formulations
1 stands for positive definite symmetric matrices.
2 stands for gerneral structurally symmetric matrices
- out_of_core = False, type bool
set **True** to enable *'out_of_core'* option of direct solver
MUMPS if main memory is exceeded. WARNING! may be REALY slow
- convert = True, type bool
automatically converts solution function *U* to real and
imarginary parts of electric and magnetic fields and saves 'xml'
and 'pvd' files. If **False**, functions *convert_results()*
and *export_all_results* in the postprocessing instance *PP* can
be called manually with custom conversion and export parameters
- delete_factorization = True, type bool
if **False**, keep solver instance with factorization of main
system for later usage; usually not required
- calculation_frequencies = [], type list of int
only compute results for frequencies specified in this list;
useful to recompute reminaing results after a random MUMPS crash
- auto_interpolate = False, type bool
set **True** to automatically interpolate solutions in parallel
after the solution and conversion phase has been finished for a
frequency (useful in multi-frequency simulations)
- interpolate_quantities = ['E_t', 'H_t'], type list of str
list of quantities for parallel interpolation on the fly, choose
from *'E_t'*, *'H_t'*, *'E_s'*, or *'H_s'*
- method = tfqmr, type str
if iterative solvers are considered in future, choose solver here
- pc = ilu, type str
if iterative solvers are considered in future, choose
preconditioner here
- secondary_potential = False, type bool
set True to use secondary potential formulation for calculating
DC fields, only relevant for **DC** approach, can also be specified
before by build_var_form() call
- post_proc_kwargs
pipe post-processing keyword arugments to post-processing instance
if automated conversion and export is enabled
"""
#if solver.lower() != 'mumps':
# lp(self.logger, 50,
# 'Error! For performance and implementation reasons, '
# 'only the direct solver MUMPS is supported so far. '
# 'Aborting ...')
# raise SystemExit
if self.solvingapproach[0]=='direct':
solver='MUMPS'
elif self.solvingapproach[0]=='iterative':
if self.solvingapproach[1]=='PRESB':
solver='PRESB'
if self.solvingapproach[2]=='directly':
inner='directly'
elif self.solvingapproach[2]=='iterative':
inner='iterative'
else:
lp(self.logger, 50, "Unknown flag. Abort simulation!")
raise SystemExit
elif self.solvingapproach[1]=='BD':
solver='BD'
if self.solvingapproach[2]=='directly':
inner='directly'
elif self.solvingapproach[2]=='iterative':
inner='iterative'
else:
lp(self.logger, 50, "Unknown flag. Abort simulation!")
raise SystemExit
else:
lp(self.logger, 50, "Unknown flag. Abort simulation!")
raise SystemExit
else:
lp(self.logger, 50, "Unknown flag. Abort simulation!")
raise SystemExit
if bc is None:
bc = self.FE.bc
if self.MP.dc:
self.FE.calc_dc_field(self.PP, self.profiler,
secondary_potential=secondary_potential,
delete_factorization=delete_factorization)
if convert:
lp(self.logger, 20,
' - automated post-processing and conversion '
'enabled - ')
self.PP.convert_results(**post_proc_kwargs)
elif self.MP.td:
if self.profiler:
export_config_file(self.PP)
if self.approach == 'E_IE':
self.FE.time_stepping()
elif self.approach == 'E_RA':
self.FE.rational_arnoldi()
elif self.approach == 'E_FT':
self.FE.calc_frequency_solutions(auto_interpolate,
calculation_frequencies,
**post_proc_kwargs)
if convert:
if self.approach in ['E_IE', 'E_RA']:
t0 = time.time()
self.PP.convert_results(**post_proc_kwargs)
if self.approach in ['E_IE', 'E_RA']:
self.FS.solution_time += (time.time() - t0)
else:
flag_mt=False
if not delete_factorization:
self.factorizations = []
for fi in range(self.MP.n_freqs):
if fi == 0:
if 'mt' in self.approach.lower():
Assembler = ce.fem.MTAssembler(self.FE, bc)
Assembler.assemble(fi=fi)
flag_mt=True
elif '_t' in self.approach:
Assembler = ce.fem.TotalFieldAssembler(self.FE, bc)
Assembler.assemble(fi=fi)
elif '_s' in self.approach:
Assembler = ce.fem.SecondaryFieldAssembler(self.FE, bc)
Assembler.assemble(fi=fi)
# used to recompute if MUMPS crashed during H-field conversion
if len(calculation_frequencies) > 0 and \
fi not in calculation_frequencies:
continue
if fi != 0:
Assembler.assemble(fi=fi)
if self.profiler:
if not os.path.isfile(self.PP.export_dir + "_config.json"):
export_config_file(self.PP)
if fi == 0:
lp(self.logger, 20,
'Solution of main FE system(s) of equations:')
if len(calculation_frequencies) > 0:
lp(self.logger, 30,
' - solving system(s) for ' +
str(len(calculation_frequencies)) +
' frequencies - ', pre_dash=False)
else:
lp(self.logger, 30,
' - solving system(s) for ' +
str(self.MP.n_freqs) +
' frequencies - ', pre_dash=False)
else:
lp(self.logger, 30,
'... solving and post-processing for frequency number'
' --> ' + str(fi) + ' ...')
if solver == 'MUMPS':
if '_s' in self.approach and not self.FS.anom_flag:
self.FS.U = [df.Function(self.FS.M) for ti in
range(self.FE.n_tx)]
# need to check this below, seems to make sense
# print('hacked no anomaly solve')
# self.Solver.solve_system_mumps(
# Assembler.A, Assembler.b,
# sym=sym, out_of_core=out_of_core)
else:
self.Solver.solve_system_mumps(
Assembler.A, Assembler.b,
sym=sym, out_of_core=out_of_core)
if convert:
if fi == 0:
lp(self.logger, 20,
' - automated post-processing and conversion '
'enabled - ', pre_dash=False)
max_mem(logger=self.logger)
self.PP.convert_results(fi=fi, **post_proc_kwargs)
if auto_interpolate:
self.IB.auto_interpolate_parallel(
interpolation_quantities, fi=fi)
if hasattr(self, 'INV'):
if fi == 0:
lp(self.PP.MP.logger, 20,
'Sensitivity Calculation (Jacobi Matrix):')
self.INV.explicit_jacobian(self.Solver, fi,
self.jacobian, inv_dir,
remote_station)
if delete_factorization:
if hasattr(self.Solver, 'solver'):
del self.Solver.solver
else:
self.factorizations.append(self.Solver.solver)
elif solver == 'PRESB':
if '_s' in self.approach and not self.FS.anom_flag:
self.FS.U = [df.Function(self.FS.M) for ti in
range(self.FE.n_tx)]
# need to check this below, seems to make sense
# print('hacked no anomaly solve')
# self.Solver.solve_system_mumps(
# Assembler.A, Assembler.b,
# sym=sym, out_of_core=out_of_core)
else:
self.Solver.solve_system_iter_PRESB(Assembler.Asys, Assembler.Hh,
Assembler.b,
Assembler.Aoffdiag,
Assembler.Adiag,inner,
method, flag_mt)
if convert:
if fi == 0:
lp(self.logger, 20,
' - automated post-processing and conversion '
'enabled - ', pre_dash=False)
max_mem(logger=self.logger)
self.PP.convert_results(fi=fi, **post_proc_kwargs)
if auto_interpolate:
self.IB.auto_interpolate_parallel(
interpolation_quantities, fi=fi)
if hasattr(self, 'INV'):
if fi == 0:
lp(self.PP.MP.logger, 20,
'Sensitivity Calculation (Jacobi Matrix):')
self.INV.explicit_jacobian(self.Solver, fi,
self.jacobian, inv_dir,
self.solvingapproach[0])
if delete_factorization:
if hasattr(self.Solver, 'solver'):
del self.Solver.solver
else:
self.factorizations.append(self.Solver.solver)
elif solver == 'BD':
if '_s' in self.approach and not self.FS.anom_flag:
self.FS.U = [df.Function(self.FS.M) for ti in
range(self.FE.n_tx)]
# need to check this below, seems to make sense
# print('hacked no anomaly solve')
# self.Solver.solve_system_mumps(
# Assembler.A, Assembler.b,
# sym=sym, out_of_core=out_of_core)
else:
self.Solver.solve_system_iter_BD(Assembler.Asys, Assembler.Hh,
Assembler.b,
Assembler.Aoffdiag,
Assembler.Adiag,inner,
method, flag_mt)
if convert:
if fi == 0:
lp(self.logger, 20,
' - automated post-processing and conversion '
'enabled - ', pre_dash=False)
max_mem(logger=self.logger)
self.PP.convert_results(fi=fi, **post_proc_kwargs)
if auto_interpolate:
self.IB.auto_interpolate_parallel(
interpolation_quantities, fi=fi)
if hasattr(self, 'INV'):
if fi == 0:
lp(self.PP.MP.logger, 20,
'Sensitivity Calculation (Jacobi Matrix):')
self.INV.explicit_jacobian(self.Solver, fi,
self.jacobian, inv_dir,
self.solvingapproach[0])
if delete_factorization:
if hasattr(self.Solver, 'solver'):
del self.Solver.solver
else:
self.factorizations.append(self.Solver.solver)
# elif solver == 'default':
# self.Solver.solve_system_default(Assembler.A,
# Assembler.b)
# elif solver == 'iter':
# self.Solver.solve_system_iter(Assembler.A, Assembler.b,
# method=method, pc=pc)
# use warning log level to suppress unnecessary prints if
# debugging level is not specified explicitly
# make code silent for repeated frequency calculations
if fi == 0 and self.MP.n_freqs != 1:
if self.logger.level > 11 and self.logger.level < 30:
self.logger.handlers[0].setLevel(30)
# switch back to specified debug level of main model
self.logger.handlers[0].setLevel(self.logger.level)
if self.profiler:
max_mem(logger=self.logger)
export_resource_file(self.PP)
[docs]
def init_third_party_parameters(self, fenics_debug_level):
"""
Set some dolfin and print parameters. Might contain more specifications
with custom choices in future.
Required arguments
------------------
- fenics_debug_level, type int
specified log level in FEniCS, default is *50*, most levels refer
to the ones of the standard Python library *logging*:
10 or "DEBUG", subdry
13 or "TRACE", what's happening (in detail)
16 or "PROGRESS", what's happening (broadly)
20 or "INFO", information of general interest
30 or "WARNING", things that may go boom later
40 or "ERROR", things that go boom
50 or "CRITICAL", ciritcal errors between go boom and beyond
"""
df.parameters['form_compiler']['optimize'] = True
df.parameters['form_compiler']['cpp_optimize'] = True
logging.getLogger('FFC').setLevel(fenics_debug_level)
logging.getLogger('UFL').setLevel(fenics_debug_level)
df.set_log_level(fenics_debug_level)
df.parameters["allow_extrapolation"] = True
# np.set_printoptions(edgeitems=5,
# threshold=11)
sys.stdout.flush()
[docs]
def init_default_model_parameters(self):
"""
Initalizes default parameters of the **MOD** class which will be
updated if keyword arguments are set when calling the MOD instance.
The parameters are described in the **MOD** class documentation.
"""
self.r_dir, self.m_dir = ce.misc.read_paths(
os.path.dirname(ce.__file__) + '/misc/paths.dat')
self.mpi_cw = df.MPI.comm_world
self.mpi_cs = df.MPI.comm_self
self.mpi_rank = df.MPI.rank(self.mpi_cw)
self.profiler = True
self.test_mode = False
self.file_format = None
self.ned_kind = '1st'
self.p = 1
self.overwrite_results = False
self.overwrite_mesh = False
self.load_existing = False
self.field_selection = 'None'
self.import_freq = None
self.fs_type = 'None'
self.mute = False
self.self_mode = False
self.para_dir = self.m_dir + '/para'
self.export_domains = True
self.dg_interpolation = False
self.mumps_debug = False
self.serial_ordering = False
self.reuse_fs = None
self.jacobian = None
self.mesh_only = False
self.iterative = False
self.solvingapproach=['direct','MUMPS','directly','H_direct']
[docs]
def init_instances(self):
"""
Initalizes instances within MOD (**FS, MP, PP, FE, Solver, IB**).
For more information, we refer to the main documentation of the
sub-modules.
"""
# don't know if this if statement is still required, but it works
if 'DC' in self.approach.upper() and self.approach != 'DC':
self.approch = 'DC'
self.file_format, dc_flag, mt_flag, td_flag, fd_flag = \
check_approach_and_file_format(
self.approach, self.file_format,
os.path.dirname(ce.__file__) + '/misc/', self.logger)
# initialize model parameters
self.MP = ce.fem.ModelParameters(self.logger,
[self.mod_name,
self.mesh_name,
self.approach,
self.file_format,
self.m_dir,
self.r_dir,
self.para_dir,
self.out_dir,
self.mute,
self.mpi_cw,
self.mpi_cs,
self.mpi_rank,
dc_flag,
mt_flag,
td_flag,
fd_flag,
self.mumps_debug,
self.serial_ordering,
self.self_mode],
self.test_mode,
self.solvingapproach[0])
# initialize function spaces
self.FS = ce.fem.FunctionSpaces(self.MP, self.p, self.test_mode,
self.self_mode, self.ned_kind,
self.dg_interpolation,
self.load_existing, self.reuse_fs,
self.overwrite_mesh, self.mesh_only)
if self.mesh_only:
return
if self.load_existing:
self.PP = ce.core.PreProcessing(
self.FS, self.MP, self.import_freq,
field_selection=self.field_selection,
self_mode=self.self_mode,
fs_type=self.fs_type)
else:
# DC field FE modeling
if dc_flag:
self.FE = ce.fem.fem_utils.DC(self.FS, self.MP)
self.Solver = ce.core.Solver(self.FS, self.FE,
mumps_debug=self.mumps_debug)
self.PP = ce.core.PostProcessingFD(self.FE,
self.export_domains)
# MT FE modeling with E-field approach
elif mt_flag:
self.FE = ce.fem.frequency_domain_approaches.E_vector(
self.FS, self.MP)
self.Solver = ce.core.Solver(self.FS, self.FE,
mumps_debug=self.mumps_debug)
self.PP = ce.core.PostProcessingFD(self.FE,
self.export_domains,
self.solvingapproach[3])
# Time domain FE approaches
elif td_flag:
td = ce.fem.time_domain_approaches
if self.approach == 'E_IE':
self.FE = td.ImplicitEuler(self.FS, self.MP)
self.time_stepping = self.FE.time_stepping
elif self.approach == 'E_RA':
self.FE = td.RationalArnoldi(self.FS, self.MP)
elif self.approach == 'E_FT':
self.FE = td.FourierTransformBased(self.FS, self.MP)
self.PP = ce.core.PostProcessingTD(self.FE,
self.export_domains)
# Frequency domain FE approaches
elif fd_flag:
fd = ce.fem.frequency_domain_approaches
if 'E' in self.approach:
self.FE = fd.E_vector(self.FS, self.MP)
elif 'H' in self.approach:
self.FE = fd.H_vector(self.FS, self.MP)
elif 'Am' in self.approach:
self.FE = fd.A_V_mixed(self.FS, self.MP)
elif 'An' in self.approach:
self.FE = fd.A_V_nodal(self.FS, self.MP)
elif'Fm' in self.approach:
self.FE = fd.F_U_mixed(self.FS, self.MP)
elif'Fn' in self.approach: # dummy, not implemented yet
self.FE = fd.F_U_nodal(self.FS, self.MP)
self.Solver = ce.core.Solver(
self.FS, self.FE, mumps_debug=self.mumps_debug,
serial_ordering=self.serial_ordering)
self.PP = ce.core.PostProcessingFD(self.FE,
self.export_domains,
self.solvingapproach[3])
# Initalize interpolation instances
if td_flag:
self.IB = ce.post.TimeDomainInterpolator(
self.PP, self.dg_interpolation,
self.self_mode)
else:
self.IB = ce.post.FrequencyDomainInterpolator(
self.PP, self.dg_interpolation,
self.self_mode, self.fs_type)
if self.jacobian is not None:
self.INV = ce.inv.InversionBase(self.PP, self.IB)
[docs]
def handle_exception(self, exc_type, exc_value, exc_traceback):
"""
Overwrite *excepthook* from *sys* module to enable raising uncaught
exceptions by the customized *self.logger* from the *logging* module
and write them to the "*mod_name* + _debug.log" file for debugging.
Required arguments
------------------
- see description of Python exceptions and corresponding traceback
"""
if issubclass(exc_type, KeyboardInterrupt):
sys.__excepthook__(exc_type, exc_value, exc_traceback)
return
self.logger.critical("Uncaught exception", exc_info=(
exc_type, exc_value, exc_traceback))