esys.downunder Package

Classes

class esys.downunder.AbstractMinimizer(J=None, m_tol=0.0001, J_tol=None, imax=300)

Bases: object

Base class for function minimization methods.

getCostFunction()

return the cost function to be minimized

Return type:CostFunction
getOptions()

Returns a dictionary of minimizer-specific options.

getResult()

Returns the result of the minimization.

logSummary()

Outputs a summary of the completed minimization process to the logger.

run(x0)

Executes the minimization algorithm for f starting with the initial guess x0.

Returns:the result of the minimization
setCallback(callback)

Sets a callback function to be called after every iteration. The arguments to the function are: (k, x, Jx, g_Jxx), where k is the current iteration, x is the current estimate, Jx=f(x) and g_Jxx=grad f(x).

setCostFunction(J)

set the cost function to be minimized

Parameters:J (CostFunction) – the cost function to be minimized
setMaxIterations(imax)

Sets the maximum number of iterations before the minimizer terminates.

setOptions(**opts)

Sets minimizer-specific options. For a list of possible options see getOptions().

setTolerance(m_tol=0.0001, J_tol=None)

Sets the tolerance for the stopping criterion. The minimizer stops when an appropriate norm is less than m_tol.

class esys.downunder.AcousticVelocityMapping(V_prior, Q_prior)

Bases: esys.downunder.mappings.Mapping

Maps a p-velocity and Q-index to slowness square sigma=(V*(1-i*1/(2*Q))^{-2} in the form sigma=e^{Mr+m[0])}*( cos(Mi+m[1])) + i * sin(Mi+m[1]) where

getDerivative(m)

returns the value for the derivative of the mapping for m

getInverse(s)

returns the value of the inverse of the mapping for s

getTypicalDerivative()

returns a typical value for the derivative

getValue(m)

returns the value of the mapping for m

class esys.downunder.AcousticWaveForm(domain, omega, w, data, F, coordinates=None, fixAtBottom=False, tol=1e-08, saveMemory=True, scaleF=True)

Bases: esys.downunder.forwardmodels.ForwardModel

Forward Model for acoustic waveform inversion in the frequence domain It defines a cost function:

defect = 1/2 integrate( ( w * ( a * u - data ) ) ** 2 )

where w are weighting factors, data are the measured data (as a 2-comp vector of real and imaginary part) for real frequency omega, and u is the coresponding result produced by the forward model. u (as a 2-comp vector) is the solution of the complex Helmholtz equation for frequency omega, source F and complex, inverse, squared p-velocity sigma:

  • -u_{ii} - omega**2 * sigma * u = F

It is assumed that the exact scale of source F is unknown and the scaling factor a of F is calculated by minimizing the defect

getArguments(sigma)

Returns precomputed values shared by getDefect() and getGradient().

Parameters:sigma (Data of shape (2,)) – a suggestion for complex 1/V**2
Returns:solution, uTar, uTai, uTu
Return type:Data of shape (2,), 3 x float
getCoordinateTransformation()

returns the coordinate transformation being used

Return type:CoordinateTransformation
getDefect(sigma, u, uTar, uTai, uTu)

Returns the defect value.

Parameters:
  • sigma (Data of shape (2,)) – a suggestion for complex 1/V**2
  • u (Data of shape (2,)) – a u vector

:param uTar : = integrate( w * (data[0]*u[0]+data[1]*u[1])) :type uTar: float :param uTai : = integrate( w * (data[1]*u[0]-data[0]*u[1])) :type uTa: float :param uTu : = integrate( w * (u,u)) :type uTu: float

Return type:float
getDomain()

Returns the domain of the forward model.

Return type:Domain
getGradient(sigma, u, uTar, uTai, uTu)

Returns the gradient of the defect with respect to density.

Parameters:
  • sigma (Data of shape (2,)) – a suggestion for complex 1/V**2
  • u (Data of shape (2,)) – a u vector

:param uTar : = integrate( w * (data[0]*u[0]+data[1]*u[1])) :type uTar: float :param uTai : = integrate( w * (data[1]*u[0]-data[0]*u[1])) :type uTa: float :param uTu : = integrate( w * (u,u)) :type uTu: float

getSourceScaling(u)

returns the scaling factor s required to rescale source F to minimize defect |s * u- data|^2 :param u: value of pressure solution (real and imaginary part) :type u: Data of shape (2,) :rtype: complex

getSurvey(index=None)

Returns the pair (data, weight)

If argument index is ignored.

rescaleWeights(scale=1.0, sigma_scale=1.0)

rescales the weights such that

integrate( ( w omega**2 * sigma_scale * data * ((1/L_j)**2)*-1) +1 )/(data*omega**2 * ((1/L_j)**2)**-1) * sigma_scale )=scale*

Parameters:
  • scale (positive float) – scale of data weighting factors
  • sigma_scale (Scalar) – scale of 1/vp**2 velocity.
setUpPDE()

Return the underlying PDE.

Return type:LinearPDE
class esys.downunder.BoundedRangeMapping(s_min=0, s_max=1)

Bases: esys.downunder.mappings.Mapping

Maps an unbounded parameter to a bounded range. The mapping is smooth and continuous.

getDerivative(m)

returns the value for the derivative of the mapping for m

getInverse(s)

returns the value of the inverse of the mapping for s

getTypicalDerivative()

returns a typical value for the derivative

getValue(m)

returns the value of the mapping for m

class esys.downunder.CartesianReferenceSystem(name='CARTESIAN')

Bases: esys.downunder.coordinates.ReferenceSystem

Identifies the Cartesian coordinate system

createTransformation(domain)

creates an appropriate coordinate transformation on a given domain

Parameters:domain (esys.escript.AbstractDomain) – domain of transformation
Return type:SpatialCoordinateTransformation
getName()

returns the name of the reference system

isCartesian()

returns if the reference system is Cartesian

Return type:bool
isTheSame(other)

test if argument other defines the same reference system

Parameters:other (ReferenceSystem) – a second reference system
Returns:True if other is a CartesianReferenceSystem instance.
Return type:bool
Note:every two CartesianReferenceSystem instances are considered as being the same.
class esys.downunder.CostFunction

Bases: object

A function f(x) that can be minimized (base class).

Example of usage:

cf=DerivedCostFunction()
# ... calculate x ...
args=cf.getArguments(x) # this could be potentially expensive!
f=cf.getValue(x, *args)
# ... it could be required to update x without using the gradient...
# ... but then ...
gf=cf.getGradient(x, *args)

The class distinguishes between the representation of the solution x (x-type) and the gradients (r-type).

Variables:provides_inverse_Hessian_approximation – This member should be set to True in subclasses that provide a valid implementation of getInverseHessianApproximation()
getArguments(x)

returns precalculated values that are shared in the calculation of f(x) and grad f(x) and the Hessian operator. The default implementation returns an empty tuple.

Note

The tuple returned by this call will be passed back to this CostFunction in other calls(eg: getGradient). Its contents are not specified at this level because no code, other than the CostFunction which created it, will be interacting with it. That is, the implementor can put whatever information they find useful in it.

Parameters:x (x-type) – location of derivative
Return type:tuple
getDualProduct(x, r)

returns the dual product of x and r

Return type:float
getGradient(x, *args)

returns the gradient of f at x using the precalculated values for x.

Parameters:
  • x (x-type) – location of derivative
  • args – pre-calculated values for x from getArguments()
Return type:

r-type

getInverseHessianApproximation(x, r, *args)

returns an approximative evaluation p of the inverse of the Hessian operator of the cost function for a given gradient r at a given location x: H(x) p = r

Parameters:
  • x (x-type) – location of Hessian operator to be evaluated
  • r (r-type) – a given gradient
  • args – pre-calculated values for x from getArguments()
Return type:

x-type

Note:

In general it is assumed that the Hessian H(x) needs to be calculated in each call for a new location x. However, the solver may suggest that this is not required, typically when the iteration is close to completeness.

Note:

Subclasses that implement this method should set the class variable provides_inverse_Hessian_approximation to True to enable the solver to call this method.

getNorm(x)

returns the norm of x

Return type:float
getValue(x, *args)

returns the value f(x) using the precalculated values for x.

Parameters:x (x-type) – a solution approximation
Return type:float
provides_inverse_Hessian_approximation = False
updateHessian()

notifies the class that the Hessian operator needs to be updated. This method is called by the solver class.

class esys.downunder.DataSource(reference_system=None, tags=[])

Bases: object

A class that provides survey data for the inversion process. This is an abstract base class that implements common functionality. Methods to be overwritten by subclasses are marked as such. This class assumes 2D data which is mapped to a slice of a 3D domain. For other setups override the methods as required.

ACOUSTIC = 2
GRAVITY = 0
MAGNETIC = 1
getDataExtents()

returns a tuple of tuples ( (x0, y0), (nx, ny), (dx, dy) ), where

  • x0, y0 = coordinates of data origin
  • nx, ny = number of data points in x and y
  • dx, dy = spacing of data points in x and y

This method must be implemented in subclasses.

getDataType()

Returns the type of survey data managed by this source. Subclasses must return GRAVITY or MAGNETIC or ACOUSTIC as appropriate.

getHeightScale()

returns the height scale factor to convert from meters to the appropriate units of the reference system used.

Return type:float
getReferenceSystem()

returns the reference coordinate system

Return type:ReferenceSystem
getSubsamplingFactor()

Returns the subsampling factor that was set via setSubsamplingFactor (see there).

getSurveyData(domain, origin, NE, spacing)

This method is called by the DomainBuilder to retrieve the survey data as Data objects on the given domain.

Subclasses should return one or more Data objects with survey data interpolated on the given escript domain. The exact return type depends on the type of data.

Parameters:
  • domain (esys.escript.Domain) – the escript domain to use
  • origin (tuple or list) – the origin coordinates of the domain
  • NE (tuple or list) – the number of domain elements in each dimension
  • spacing (tuple or list) – the cell sizes (node spacing) in the domain
getTags()

returns the list of tags

Return type:list
getUtmZone()

All data source coordinates are converted to UTM (Universal Transverse Mercator) in order to have useful domain extents. Subclasses should implement this method and return the UTM zone number of the projected coordinates.

hasTag(tag)

returns true if the data set has tag tag

Return type:bool
setSubsamplingFactor(f)

Sets the data subsampling factor (default=1).

The factor is applied in all dimensions. For example a 2D dataset with 300 x 150 data points will be reduced to 150 x 75 when a subsampling factor of 2 is used. This becomes important when adding data of varying resolution to a DomainBuilder.

class esys.downunder.DensityMapping(domain, z0=None, rho0=None, drho=None, beta=None)

Bases: esys.downunder.mappings.LinearMapping

Density mapping with depth weighting

rho = rho0 + drho * ( (x_2 - z0)/l_z)^(beta/2) ) * m

getDerivative(m)

returns the value for the derivative of the mapping for m

getInverse(p)

returns the value of the inverse of the mapping for s

getTypicalDerivative()

returns a typical value for the derivative

getValue(m)

returns the value of the mapping for m

class esys.downunder.DomainBuilder(dim=3, reference_system=None)

Bases: object

This class is responsible for constructing an escript Domain object with suitable extents and resolution for survey data (DataSource objects) that are added to it.

The domain covers a region above and below the Earth surface. The East-West direction is used as the x- or longitudinal or x[0] direction, the North-South direction is used as the y- or latitudinal or x[1] direction, the vertical direction is denoted by z or radial or x[2] direction. The corresponding terms are used synonymously.

addSource(source)

Adds a survey data provider to the domain builder. An exception is raised if the domain has already been built or if the UTM zone of source does not match the UTM zone of sources already added to the domain builder (see Inversion Cookbook for more information). An exception is also raised if the dimensionality of the data source is incompatible with this domain builder. That is, the dimensionality of the data must be one less than the dimensionality of the domain (specified in the constructor).

Parameters:source (DataSource) – The data source to be added. Its reference system needs to match the reference system of the DomainBuilder.
fixDensityBelow(depth=None)

Defines the depth below which the density anomaly is set to a given value. If no value is given zero is assumed.

Parameters:depth (float) – depth below which the density is fixed. If not set, no constraint at depth is applied.
fixSusceptibilityBelow(depth=None)

Defines the depth below which the susceptibility anomaly is set to a given value. If no value is given zero is assumed.

Parameters:depth (float) – depth below which the susceptibility is fixed. If not set, no constraint at depth is applied.
fixVelocityBelow(depth=None)

Defines the depth below which the velocity and Q index is set to a given value. If no value is given zero is assumed.

Parameters:depth (float) – depth below which the velocity is fixed. If not set, no constraint at depth is applied.
getBackgroundMagneticFluxDensity()

Returns the background magnetic flux density.

getDomain()

Returns a domain that spans the data area plus padding.

The domain is created the first time this method is called, subsequent calls return the same domain so anything that affects the domain (such as padding) needs to be set beforehand.

Returns:The escript domain for this data source
Return type:esys.escript.Domain
getGravitySurveys()

Returns a list of gravity surveys, see getSurveys for details.

getMagneticSurveys()

Returns a list of magnetic surveys, see getSurveys for details.

getReferenceSystem()

returns the reference coordinate system

Return type:ReferenceSystem
getSetDensityMask()

Returns the density mask data object which is non-zero for cells whose density value is fixed, zero otherwise.

getSetSusceptibilityMask()

Returns the susceptibility mask data object which is non-zero for cells whose susceptibility value is fixed, zero otherwise.

getSurveys(datatype, tags=None)

Returns a list of Data objects for all surveys of type datatype available to this domain builder. If a list of tags is given only data sources whose tag matching the tag list are returned

Returns:List of surveys which are tuples (anomaly,error).
Return type:list
getTags()

returns a list of all tags beeing used by the attached data sources. the list may be empty.

setBackgroundMagneticFluxDensity(B)

Sets the background magnetic flux density B=(B_East, B_North, B_Vertical)

setElementPadding(pad_x=None, pad_y=None, pad_lat=None, pad_lon=None)

Sets the amount of padding around the dataset in number of elements (cells).

When the domain is constructed pad_x (pad_y) elements are added on each side of the x- (y-) dimension. The arguments must be non-negative.

Parameters:
  • pad_x (int) – Padding per side in x direction (default: no padding)
  • pad_y (int) – Padding per side in y direction (default: no padding)
Note:

pad_y is ignored for 2-dimensional datasets.

setFractionalPadding(pad_x=None, pad_y=None, pad_lat=None, pad_lon=None)

Sets the amount of padding around the dataset as a fraction of the dataset side lengths.

For example, calling setFractionalPadding(0.2, 0.1) with a data source of size 10x20 will result in the padded data set size 14x24 (10*(1+2*0.2), 20*(1+2*0.1))

Parameters:
  • pad_x (float) – Padding per side in x direction (default: no padding)
  • pad_y (float) – Padding per side in y direction (default: no padding)
  • pad_lat (float) – Padding per side in latitudinal direction (default: no padding)
  • pad_lon (float) – Padding per side in longitudinal direction (default: no padding)
Note:

pad_y is ignored for 2-dimensional domains.

setGeoPadding(pad_lat=None, pad_lon=None)

Sets the amount of padding around the dataset in longitude and latitude.

The final domain size will be the extent in the latitudinal (in longitudinal) direction of the dataset plus twice the value of pad_lat (pad_lon). The arguments must be non-negative.

Parameters:
  • pad_lat (float in units of degree) – Padding per side in latitudinal direction (default: 0)
  • pad_lon (float in units of degree) – Padding per side in longitudinal direction (default: 0)
Note:

pad_lon is ignored for 2-dimensional domains.

Note:

this function can only be used if the reference system is not Cartesian

setPadding(pad_x=None, pad_y=None, pad_lat=None, pad_lon=None)

Sets the amount of padding around the dataset in absolute length units.

The final domain size will be the length in x (in y) of the dataset plus twice the value of pad_x (pad_y). The arguments must be non-negative.

Parameters:
  • pad_x (float in units of length (meter)) – Padding per side in x direction (default: no padding)
  • pad_y (float in units of length (meter)) – Padding per side in y direction (default: no padding)
Note:

pad_y is ignored for 2-dimensional domains.

Note:

this function can only be used if the reference system is Cartesian

setVerticalExtents(depth=40000.0, air_layer=10000.0, num_cells=25)

This method sets the target domain parameters for the vertical dimension.

Parameters:
  • depth (float) – Depth of the domain (in meters)
  • air_layer (float) – Depth of the layer above sea level (in meters)
  • num_cells (int) – Number of domain elements for the entire vertical dimension
class esys.downunder.ErMapperData(data_type, headerfile, datafile=None, altitude=0.0, error=None, scale_factor=None, null_value=None, reference_system=None)

Bases: esys.downunder.datasources.DataSource

Data Source for ER Mapper raster data. Note that this class only accepts a very specific type of ER Mapper data input and will raise an exception if other data is found.

ACOUSTIC = 2
GRAVITY = 0
MAGNETIC = 1
getDataExtents()

returns ( (x0, y0), (nx, ny), (dx, dy) )

getDataType()
getHeightScale()

returns the height scale factor to convert from meters to the appropriate units of the reference system used.

Return type:float
getReferenceSystem()

returns the reference coordinate system

Return type:ReferenceSystem
getSubsamplingFactor()

Returns the subsampling factor that was set via setSubsamplingFactor (see there).

getSurveyData(domain, origin, NE, spacing)
getTags()

returns the list of tags

Return type:list
getUtmZone()
hasTag(tag)

returns true if the data set has tag tag

Return type:bool
setSubsamplingFactor(f)

Sets the data subsampling factor (default=1).

The factor is applied in all dimensions. For example a 2D dataset with 300 x 150 data points will be reduced to 150 x 75 when a subsampling factor of 2 is used. This becomes important when adding data of varying resolution to a DomainBuilder.

class esys.downunder.ForwardModel

Bases: object

An abstract forward model that can be plugged into a cost function. Subclasses need to implement getDefect(), getGradient(), and possibly getArguments() and ‘getCoordinateTransformation’.

getArguments(x)
getCoordinateTransformation()
getDefect(x, *args)
getGradient(x, *args)
class esys.downunder.ForwardModelWithPotential(domain, w, data, coordinates=None, fixPotentialAtBottom=False, tol=1e-08)

Bases: esys.downunder.forwardmodels.ForwardModel

Base class for a forward model using a potential such as magnetic or gravity. It defines a cost function:

defect = 1/2 sum_s integrate( ( weight_i[s] * ( r_i - data_i[s] ) ) ** 2 )

where s runs over the survey, weight_i are weighting factors, data_i are the data, and r_i are the results produced by the forward model. It is assumed that the forward model is produced through postprocessing of the solution of a potential PDE.

getArguments(x)
getCoordinateTransformation()

returns the coordinate transformation being used

Return type:CoordinateTransformation
getDefect(x, *args)
getDefectGradient(result)
getDomain()

Returns the domain of the forward model.

Return type:Domain
getGradient(x, *args)
getPDE()

Return the underlying PDE.

Return type:LinearPDE
getSurvey(index=None)

Returns the pair (data_index, weight_index), where data_i is the data of survey i, weight_i is the weighting factor for survey i. If index is None, all surveys will be returned in a pair of lists.

class esys.downunder.GeodeticCoordinateTransformation(domain, reference=<esys.downunder.coordinates.GeodeticReferenceSystem object at 0x7fd148998590>)

Bases: esys.downunder.coordinates.SpatialCoordinateTransformation

A geodetic coordinate transformation

getDomain()

returns the domain of the coordinate transformation.

Return type:esys.escript.AbstractDomain
getReferenceSystem()

returns the reference system used to to define the coordinate transformation

Return type:ReferenceSystem
getScalingFactors()

returns the scaling factors

Return type:esys.escript.Vector
getVolumeFactor()

returns the volume factor for the coordinate transformation

Return type:esys.escript.Scalar
isCartesian()

returns True if the scaling factors (and the volume factor) are equal to 1

Return type:bool
isTheSame(other)

test if argument other defines the same coordinate transformation

Parameters:other (SpatialCoordinateTransformation) – a second coordinate transformation
Returns:True if other defines then same coordinate transformation
Return type:bool
class esys.downunder.GeodeticReferenceSystem(a=6378137.0, f=0.0033528106647474805, angular_unit=0.017453292519943295, height_unit=1000.0, name='WGS84')

Bases: esys.downunder.coordinates.ReferenceSystem

Identifies a Geodetic coordinate system

createTransformation(domain)

creates an appropriate coordinate transformation on a given domain

Parameters:domain (esys.escript.AbstractDomain) – domain of transformation
Return type:SpatialCoordinateTransformation
getAngularUnit()

returns the angular unit

getFlattening()

returns the flattening

getHeightUnit()

returns the height unit

getName()

returns the name of the reference system

getSemiMajorAxis()

returns the length of semi major axis

getSemiMinorAxis()

returns the length of semi minor axis

isCartesian()

returns if the reference system is Cartesian

Return type:bool
isTheSame(other)

test if other defines the same reference system

Parameters:other (ReferenceSystem) – a second reference system
Returns:True if other defines then same reference system
Return type:bool

Note

two GeodeticReferenceSystem are considered to be the same if the use the same semi major axis, the same flattening and the same angular unit.

class esys.downunder.GravityInversion(solverclass=None, debug=False, self_demagnetization=False)

Bases: esys.downunder.inversions.InversionDriver

Driver class to perform an inversion of Gravity (Bouguer) anomaly data. The class uses the standard Regularization class for a single level set function, DensityMapping mapping, and the gravity forward model GravityModel.

fixGravityPotentialAtBottom(status=False)

indicates to fix the gravity potential at the bottom to zero (in addition to the top)

Parameters:status (bool) – if True gravity potential at the bottom is set to zero
fixMagneticPotentialAtBottom(status=True)

indicates to fix the magnetic potential at the bottom to zero (in addition to the top)

Parameters:status (bool) – if True magnetic potential at the bottom is set to zero
getCostFunction()

returns the domain of the inversion

Return type:‘InversionCostFunction’
getDomain()

returns the domain of the inversion

Return type:Domain
getLevelSetFunction()

returns the level set function as solution of the optimization problem

Return type:Data
getSolver()

The solver to be used in the inversion process. See the minimizers module for available solvers. By default, the L-BFGS minimizer is used.

Return type:‘AbstractMinimizer’.
isSetUp()

returns True if the inversion is set up and is ready to run.

Return type:bool
run()

This function runs the inversion.

Returns:physical parameters as result of the inversion
Return type:list of physical parameters or a physical parameter
setCostFunction(costfunction)

sets the cost function of the inversion. This function needs to be called before the inversion iteration can be started.

Parameters:costfunction (‘InversionCostFunction’) – domain of the inversion
setInitialGuess(rho=None)

set the initial guess rho for density the inversion iteration. If no rho present then an appropriate initial guess is chosen.

Parameters:rho (Scalar) – initial value for the density anomaly.
setSolverCallback(callback=None)

Sets the callback function which is called after every solver iteration.

setSolverMaxIterations(maxiter=None)

Sets the maximum number of solver iterations to run. If maxiter is reached the iteration is terminated and MinimizerMaxIterReached is thrown.

Parameters:maxiter (positive int) – maximum number of iteration steps.
setSolverTolerance(m_tol=None, J_tol=None)

Sets the error tolerance for the solver. An acceptable solution is considered to be found once the tolerance is reached.

Parameters:
  • m_tol (float or None) – tolerance for changes to level set function. If None changes to the level set function are not checked for convergence during iteration.
  • J_tol (float or None) – tolerance for changes to cost function. If None changes to the cost function are not checked for convergence during iteration.
Note:

if both arguments are equal to None the default setting m_tol=1e-4, J_tol=None is used.

setup(domainbuilder, rho0=None, drho=None, z0=None, beta=None, w0=None, w1=None, rho_at_depth=None)

Sets up the inversion parameters from a DomainBuilder.

Parameters:
  • domainbuilder (DomainBuilder) – Domain builder object with gravity source(s)
  • rho0 (float or Scalar) – reference density, see DensityMapping. If not specified, zero is used.
  • drho (float or Scalar) – density scale, see DensityMapping. If not specified, 2750kg/m^3 is used.
  • z0 (float or Scalar) – reference depth for depth weighting, see DensityMapping. If not specified, zero is used.
  • beta (float or Scalar) – exponent for depth weighting, see DensityMapping. If not specified, zero is used.
  • w0 (Scalar or float) – weighting factor for level set term regularization. If not set zero is assumed.
  • w1 (Vector or list of float) – weighting factor for the gradient term in the regularization. If not set zero is assumed
  • rho_at_depth (float or None) – value for density at depth, see DomainBuilder.
siloWriterCallback(k, m, Jm, g_Jm)

callback function that can be used to track the solution

Parameters:
  • k – iteration count
  • m – current m approximation
  • Jm – value of cost function
  • g_Jm – gradient of f at x
class esys.downunder.GravityModel(domain, w, g, gravity_constant=6.6742e-11, coordinates=None, fixPotentialAtBottom=False, tol=1e-08)

Bases: esys.downunder.forwardmodels.ForwardModelWithPotential

Forward Model for gravity inversion as described in the inversion cookbook.

getArguments(rho)

Returns precomputed values shared by getDefect() and getGradient().

Parameters:rho (Scalar) – a suggestion for the density distribution
Returns:gravity potential and corresponding gravity field.
Return type:Scalar, Vector
getCoordinateTransformation()

returns the coordinate transformation being used

Return type:CoordinateTransformation
getDefect(rho, phi, gravity_force)

Returns the value of the defect

Parameters:
  • rho (Scalar) – density distribution
  • phi (Scalar) – corresponding potential
  • gravity_force (Vector) – gravity force
Return type:

float

getDefectGradient(result)
getDomain()

Returns the domain of the forward model.

Return type:Domain
getGradient(rho, phi, gravity_force)

Returns the gradient of the defect with respect to density.

Parameters:
  • rho (Scalar) – density distribution
  • phi (Scalar) – corresponding potential
  • gravity_force (Vector) – gravity force
Return type:

Scalar

getPDE()

Return the underlying PDE.

Return type:LinearPDE
getPotential(rho)

Calculates the gravity potential for a given density distribution.

Parameters:rho (Scalar) – a suggestion for the density distribution
Returns:gravity potential
Return type:Scalar
getSurvey(index=None)

Returns the pair (data_index, weight_index), where data_i is the data of survey i, weight_i is the weighting factor for survey i. If index is None, all surveys will be returned in a pair of lists.

rescaleWeights(scale=1.0, rho_scale=1.0)

rescales the weights such that

sum_s integrate( ( w_i[s] *g_i[s]) (w_j[s]*1/L_j) * L**2 * 4*pi*G*rho_scale )=scale

Parameters:
  • scale (positive float) – scale of data weighting factors
  • rho_scale (Scalar) – scale of density.
class esys.downunder.HTIWave(domain, v_p, v_s, wavelet, source_tag, source_vector=[1.0, 0.0, 0.0], eps=0.0, gamma=0.0, delta=0.0, rho=1.0, dt=None, u0=None, v0=None, absorption_zone=300.0, absorption_cut=0.01, lumping=True, useFast=False)

Bases: esys.downunder.seismic.WaveBase

Solving the HTI wave equation (along the x_0 axis)

Note:In case of a two dimensional domain a horizontal domain is considered, i.e. the depth component is dropped.
getTimeStepSize()
update(t)

returns the solution for the next time marker t which needs to greater than the time marker from the previous call.

class esys.downunder.InversionCostFunction(regularization, mappings, forward_models)

Bases: esys.downunder.costfunctions.MeteredCostFunction

Class to define cost function J(m) for inversion with one or more forward models based on a multi-valued level set function m:

J(m) = J_reg(m) + sum_f mu_f * J_f(p)

where J_reg(m) is the regularization and cross gradient component of the cost function applied to a level set function m, J_f(p) are the data defect cost functions involving a physical forward model using the physical parameter(s) p and mu_f is the trade-off factor for model f.

A forward model depends on a set of physical parameters p which are constructed from components of the level set function m via mappings.

Example 1 (single forward model):
m=Mapping() f=ForwardModel() J=InversionCostFunction(Regularization(), m, f)
Example 2 (two forward models on a single valued level set)

m0=Mapping() m1=Mapping() f0=ForwardModel() f1=ForwardModel()

J=InversionCostFunction(Regularization(), mappings=[m0, m1], forward_models=[(f0, 0), (f1,1)])

Example 3 (two forward models on 2-valued level set)

m0=Mapping() m1=Mapping() f0=ForwardModel() f1=ForwardModel()

J=InversionCostFunction(Regularization(self.numLevelSets=2), mappings=[(m0,0), (m1,0)], forward_models=[(f0, 0), (f1,1)])

Variables:provides_inverse_Hessian_approximation – if true the class provides an approximative inverse of the Hessian operator.
createLevelSetFunction(*props)

returns an instance of an object used to represent a level set function initialized with zeros. Components can be overwritten by physical properties props. If present entries must correspond to the mappings arguments in the constructor. Use None for properties for which no value is given.

getArguments(x)

returns precalculated values that are shared in the calculation of f(x) and grad f(x) and the Hessian operator

Note

The tuple returned by this call will be passed back to this CostFunction in other calls(eg: getGradient). Its contents are not specified at this level because no code, other than the CostFunction which created it, will be interacting with it. That is, the implementor can put whatever information they find useful in it.

Parameters:x (x-type) – location of derivative
Return type:tuple
getComponentValues(m, *args)
getDomain()

returns the domain of the cost function

Return type:Domain
getDualProduct(x, r)

returns the dual product of x and r

Return type:float
getForwardModel(idx=None)

returns the idx-th forward model.

Parameters:idx (int) – model index. If cost function contains one model only idx can be omitted.
getGradient(x, *args)

returns the gradient of f at x using the precalculated values for x.

Parameters:
  • x (x-type) – location of derivative
  • args – pre-calculated values for x from getArguments()
Return type:

r-type

getInverseHessianApproximation(x, r, *args)

returns an approximative evaluation p of the inverse of the Hessian operator of the cost function for a given gradient r at a given location x: H(x) p = r

Note

In general it is assumed that the Hessian H(x) needs to be calculate in each call for a new location x. However, the solver may suggest that this is not required, typically when the iteration is close to completeness.

Parameters:
  • x (x-type) – location of Hessian operator to be evaluated.
  • r (r-type) – a given gradient
  • args – pre-calculated values for x from getArguments()
Return type:

x-type

getNorm(x)

returns the norm of x

Return type:float
getNumTradeOffFactors()

returns the number of trade-off factors being used including the trade-off factors used in the regularization component.

Return type:int
getProperties(m, return_list=False)

returns a list of the physical properties from a given level set function m using the mappings of the cost function.

Parameters:
  • m (Data) – level set function
  • return_list (bool) – if True a list is returned.
Return type:

list of Data

getRegularization()

returns the regularization

Return type:Regularization
getTradeOffFactors(mu=None)

returns a list of the trade-off factors.

Return type:list of float
getTradeOffFactorsModels()

returns the trade-off factors for the forward models

Return type:float or list of float
getValue(x, *args)

returns the value f(x) using the precalculated values for x.

Parameters:x (x-type) – a solution approximation
Return type:float
provides_inverse_Hessian_approximation = True
resetCounters()

resets all statistical counters

setTradeOffFactors(mu=None)

sets the trade-off factors for the forward model and regularization terms.

Parameters:mu (list of float) – list of trade-off factors.
setTradeOffFactorsModels(mu=None)

sets the trade-off factors for the forward model components.

Parameters:mu (float in case of a single model or a list of float with the length of the number of models.) – list of the trade-off factors. If not present ones are used.
setTradeOffFactorsRegularization(mu=None, mu_c=None)

sets the trade-off factors for the regularization component of the cost function, see Regularization for details.

Parameters:
  • mu – trade-off factors for the level-set variation part
  • mu_c – trade-off factors for the cross gradient variation part
updateHessian()

notifies the class that the Hessian operator needs to be updated.

class esys.downunder.InversionDriver(solverclass=None, debug=False, self_demagnetization=False)

Bases: object

Base class for running an inversion

fixGravityPotentialAtBottom(status=False)

indicates to fix the gravity potential at the bottom to zero (in addition to the top)

Parameters:status (bool) – if True gravity potential at the bottom is set to zero
fixMagneticPotentialAtBottom(status=True)

indicates to fix the magnetic potential at the bottom to zero (in addition to the top)

Parameters:status (bool) – if True magnetic potential at the bottom is set to zero
getCostFunction()

returns the domain of the inversion

Return type:‘InversionCostFunction’
getDomain()

returns the domain of the inversion

Return type:Domain
getLevelSetFunction()

returns the level set function as solution of the optimization problem

Return type:Data
getSolver()

The solver to be used in the inversion process. See the minimizers module for available solvers. By default, the L-BFGS minimizer is used.

Return type:‘AbstractMinimizer’.
isSetUp()

returns True if the inversion is set up and is ready to run.

Return type:bool
run()

This function runs the inversion.

Returns:physical parameters as result of the inversion
Return type:list of physical parameters or a physical parameter
setCostFunction(costfunction)

sets the cost function of the inversion. This function needs to be called before the inversion iteration can be started.

Parameters:costfunction (‘InversionCostFunction’) – domain of the inversion
setInitialGuess(*props)

Sets the initial guess for the inversion iteration. By default zero is used.

setSolverCallback(callback=None)

Sets the callback function which is called after every solver iteration.

setSolverMaxIterations(maxiter=None)

Sets the maximum number of solver iterations to run. If maxiter is reached the iteration is terminated and MinimizerMaxIterReached is thrown.

Parameters:maxiter (positive int) – maximum number of iteration steps.
setSolverTolerance(m_tol=None, J_tol=None)

Sets the error tolerance for the solver. An acceptable solution is considered to be found once the tolerance is reached.

Parameters:
  • m_tol (float or None) – tolerance for changes to level set function. If None changes to the level set function are not checked for convergence during iteration.
  • J_tol (float or None) – tolerance for changes to cost function. If None changes to the cost function are not checked for convergence during iteration.
Note:

if both arguments are equal to None the default setting m_tol=1e-4, J_tol=None is used.

setup(*args, **k_args)

returns True if the inversion is set up and ready to run.

class esys.downunder.JointGravityMagneticInversion(solverclass=None, debug=False, self_demagnetization=False)

Bases: esys.downunder.inversions.InversionDriver

Driver class to perform a joint inversion of Gravity (Bouguer) and magnetic anomaly data. The class uses the standard Regularization class for two level set functions with cross-gradient correlation, DensityMapping and SusceptibilityMapping mappings, the gravity forward model GravityModel and the linear magnetic forward model MagneticModel.

DENSITY = 0
SUSCEPTIBILITY = 1
fixGravityPotentialAtBottom(status=False)

indicates to fix the gravity potential at the bottom to zero (in addition to the top)

Parameters:status (bool) – if True gravity potential at the bottom is set to zero
fixMagneticPotentialAtBottom(status=True)

indicates to fix the magnetic potential at the bottom to zero (in addition to the top)

Parameters:status (bool) – if True magnetic potential at the bottom is set to zero
getCostFunction()

returns the domain of the inversion

Return type:‘InversionCostFunction’
getDomain()

returns the domain of the inversion

Return type:Domain
getLevelSetFunction()

returns the level set function as solution of the optimization problem

Return type:Data
getSolver()

The solver to be used in the inversion process. See the minimizers module for available solvers. By default, the L-BFGS minimizer is used.

Return type:‘AbstractMinimizer’.
isSetUp()

returns True if the inversion is set up and is ready to run.

Return type:bool
run()

This function runs the inversion.

Returns:physical parameters as result of the inversion
Return type:list of physical parameters or a physical parameter
setCostFunction(costfunction)

sets the cost function of the inversion. This function needs to be called before the inversion iteration can be started.

Parameters:costfunction (‘InversionCostFunction’) – domain of the inversion
setInitialGuess(rho=None, k=None)

set the initial guess rho for density and k for susceptibility for the inversion iteration.

Parameters:
  • rho (Scalar) – initial value for the density anomaly.
  • k (Scalar) – initial value for the susceptibility anomaly.
setSolverCallback(callback=None)

Sets the callback function which is called after every solver iteration.

setSolverMaxIterations(maxiter=None)

Sets the maximum number of solver iterations to run. If maxiter is reached the iteration is terminated and MinimizerMaxIterReached is thrown.

Parameters:maxiter (positive int) – maximum number of iteration steps.
setSolverTolerance(m_tol=None, J_tol=None)

Sets the error tolerance for the solver. An acceptable solution is considered to be found once the tolerance is reached.

Parameters:
  • m_tol (float or None) – tolerance for changes to level set function. If None changes to the level set function are not checked for convergence during iteration.
  • J_tol (float or None) – tolerance for changes to cost function. If None changes to the cost function are not checked for convergence during iteration.
Note:

if both arguments are equal to None the default setting m_tol=1e-4, J_tol=None is used.

setup(domainbuilder, rho0=None, drho=None, rho_z0=None, rho_beta=None, k0=None, dk=None, k_z0=None, k_beta=None, w0=None, w1=None, w_gc=None, rho_at_depth=None, k_at_depth=None)

Sets up the inversion from an instance domainbuilder of a DomainBuilder. Gravity and magnetic data attached to the domainbuilder are considered in the inversion. If magnetic data are given as scalar it is assumed that values are collected in direction of the background magnetic field.

Parameters:
  • domainbuilder (DomainBuilder) – Domain builder object with gravity source(s)
  • rho0 (float or Scalar) – reference density, see DensityMapping. If not specified, zero is used.
  • drho (float or Scalar) – density scale, see DensityMapping. If not specified, 2750kg/m^3 is used.
  • rho_z0 (float or Scalar) – reference depth for depth weighting for density, see DensityMapping. If not specified, zero is used.
  • rho_beta (float or Scalar) – exponent for depth weighting for density, see DensityMapping. If not specified, zero is used.
  • k0 (float or Scalar) – reference susceptibility, see SusceptibilityMapping. If not specified, zero is used.
  • dk (float or Scalar) – susceptibility scale, see SusceptibilityMapping. If not specified, 2750kg/m^3 is used.
  • k_z0 (float or Scalar) – reference depth for depth weighting for susceptibility, see SusceptibilityMapping. If not specified, zero is used.
  • k_beta (float or Scalar) – exponent for depth weighting for susceptibility, see SusceptibilityMapping. If not specified, zero is used.
  • w0 (Data or ndarray of shape (2,)) – weighting factors for level set term regularization, see Regularization. If not set zero is assumed.
  • w1 (Data or ndarray of shape (2,DIM)) – weighting factor for the gradient term in the regularization see Regularization. If not set zero is assumed
  • w_gc (Scalar or float) – weighting factor for the cross gradient term in the regularization, see Regularization. If not set one is assumed
  • k_at_depth (float or None) – value for susceptibility at depth, see DomainBuilder.
  • rho_at_depth (float or None) – value for density at depth, see DomainBuilder.
siloWriterCallback(k, m, Jm, g_Jm)

callback function that can be used to track the solution

Parameters:
  • k – iteration count
  • m – current m approximation
  • Jm – value of cost function
  • g_Jm – gradient of f at x
class esys.downunder.LinearMapping(a=1, p0=0)

Bases: esys.downunder.mappings.Mapping

Maps a parameter by a linear transformation p = a * m + p0

getDerivative(m)

returns the value for the derivative of the mapping for m

getInverse(p)

returns the value of the inverse of the mapping for s

getTypicalDerivative()

returns a typical value for the derivative

getValue(m)

returns the value of the mapping for m

class esys.downunder.MagneticInversion(solverclass=None, debug=False, self_demagnetization=False)

Bases: esys.downunder.inversions.InversionDriver

Driver class to perform an inversion of magnetic anomaly data. The class uses the standard Regularization class for a single level set function, SusceptibilityMapping mapping and the linear magnetic forward model MagneticModel.

fixGravityPotentialAtBottom(status=False)

indicates to fix the gravity potential at the bottom to zero (in addition to the top)

Parameters:status (bool) – if True gravity potential at the bottom is set to zero
fixMagneticPotentialAtBottom(status=True)

indicates to fix the magnetic potential at the bottom to zero (in addition to the top)

Parameters:status (bool) – if True magnetic potential at the bottom is set to zero
getCostFunction()

returns the domain of the inversion

Return type:‘InversionCostFunction’
getDomain()

returns the domain of the inversion

Return type:Domain
getLevelSetFunction()

returns the level set function as solution of the optimization problem

Return type:Data
getSolver()

The solver to be used in the inversion process. See the minimizers module for available solvers. By default, the L-BFGS minimizer is used.

Return type:‘AbstractMinimizer’.
isSetUp()

returns True if the inversion is set up and is ready to run.

Return type:bool
run()

This function runs the inversion.

Returns:physical parameters as result of the inversion
Return type:list of physical parameters or a physical parameter
setCostFunction(costfunction)

sets the cost function of the inversion. This function needs to be called before the inversion iteration can be started.

Parameters:costfunction (‘InversionCostFunction’) – domain of the inversion
setInitialGuess(k=None)

set the initial guess k for susceptibility for the inversion iteration. If no k present then an appropriate initial guess is chosen.

Parameters:k (Scalar) – initial value for the susceptibility anomaly.
setSolverCallback(callback=None)

Sets the callback function which is called after every solver iteration.

setSolverMaxIterations(maxiter=None)

Sets the maximum number of solver iterations to run. If maxiter is reached the iteration is terminated and MinimizerMaxIterReached is thrown.

Parameters:maxiter (positive int) – maximum number of iteration steps.
setSolverTolerance(m_tol=None, J_tol=None)

Sets the error tolerance for the solver. An acceptable solution is considered to be found once the tolerance is reached.

Parameters:
  • m_tol (float or None) – tolerance for changes to level set function. If None changes to the level set function are not checked for convergence during iteration.
  • J_tol (float or None) – tolerance for changes to cost function. If None changes to the cost function are not checked for convergence during iteration.
Note:

if both arguments are equal to None the default setting m_tol=1e-4, J_tol=None is used.

setup(domainbuilder, k0=None, dk=None, z0=None, beta=None, w0=None, w1=None, k_at_depth=None)

Sets up the inversion from a DomainBuilder. If magnetic data are given as scalar it is assumed that values are collected in direction of the background magnetic field.

Parameters:
  • domainbuilder (DomainBuilder) – Domain builder object with gravity source(s)
  • k0 (float or Scalar) – reference susceptibility, see SusceptibilityMapping. If not specified, zero is used.
  • dk (float or Scalar) – susceptibility scale, see SusceptibilityMapping. If not specified, 2750kg/m^3 is used.
  • z0 (float or Scalar) – reference depth for depth weighting, see SusceptibilityMapping. If not specified, zero is used.
  • beta (float or Scalar) – exponent for depth weighting, see SusceptibilityMapping. If not specified, zero is used.
  • w0 (Scalar or float) – weighting factor for level set term regularization. If not set zero is assumed.
  • w1 (Vector or list of float) – weighting factor for the gradient term in the regularization. If not set zero is assumed
  • k_at_depth (float or None) – value for susceptibility at depth, see DomainBuilder.
siloWriterCallback(k, m, Jm, g_Jm)

callback function that can be used to track the solution

Parameters:
  • k – iteration count
  • m – current m approximation
  • Jm – value of cost function
  • g_Jm – gradient of f at x
class esys.downunder.MagneticModel(domain, w, B, background_magnetic_flux_density, coordinates=None, fixPotentialAtBottom=False, tol=1e-08)

Bases: esys.downunder.forwardmodels.ForwardModelWithPotential

Forward Model for magnetic inversion as described in the inversion cookbook.

getArguments(k)

Returns precomputed values shared by getDefect() and getGradient().

Parameters:k (Scalar) – susceptibility
Returns:scalar magnetic potential and corresponding magnetic field
Return type:Scalar, Vector
getCoordinateTransformation()

returns the coordinate transformation being used

Return type:CoordinateTransformation
getDefect(k, phi, magnetic_flux_density)

Returns the value of the defect.

Parameters:
  • k (Scalar) – susceptibility
  • phi (Scalar) – corresponding potential
  • magnetic_flux_density (Vector) – magnetic field
Return type:

float

getDefectGradient(result)
getDomain()

Returns the domain of the forward model.

Return type:Domain
getGradient(k, phi, magnetic_flux_density)

Returns the gradient of the defect with respect to susceptibility.

Parameters:
  • k (Scalar) – susceptibility
  • phi (Scalar) – corresponding potential
  • magnetic_flux_density (Vector) – magnetic field
Return type:

Scalar

getPDE()

Return the underlying PDE.

Return type:LinearPDE
getPotential(k)

Calculates the magnetic potential for a given susceptibility.

Parameters:k (Scalar) – susceptibility
Returns:magnetic potential
Return type:Scalar
getSurvey(index=None)

Returns the pair (data_index, weight_index), where data_i is the data of survey i, weight_i is the weighting factor for survey i. If index is None, all surveys will be returned in a pair of lists.

rescaleWeights(scale=1.0, k_scale=1.0)

rescales the weights such that

sum_s integrate( ( w_i[s] *B_i[s]) (w_j[s]*1/L_j) * L**2 * (background_magnetic_flux_density_j[s]*1/L_j) * k_scale )=scale

Parameters:
  • scale (positive float) – scale of data weighting factors
  • k_scale (Scalar) – scale of susceptibility.
class esys.downunder.Mapping(*args)

Bases: object

An abstract mapping class to map level set functions m to physical parameters p.

getDerivative(m)

returns the value for the derivative of the mapping for m

getInverse(s)

returns the value of the inverse of the mapping for physical parameter p

getTypicalDerivative()

returns a typical value for the derivative

getValue(m)

returns the value of the mapping for m

class esys.downunder.MeteredCostFunction

Bases: esys.downunder.costfunctions.CostFunction

This an intrumented version of the CostFunction class. The function calls update statistical information. The actual work is done by the methods with corresponding name and a leading underscore. These functions need to be overwritten for a particular cost function implementation.

getArguments(x)

returns precalculated values that are shared in the calculation of f(x) and grad f(x) and the Hessian operator

Note

The tuple returned by this call will be passed back to this CostFunction in other calls(eg: getGradient). Its contents are not specified at this level because no code, other than the CostFunction which created it, will be interacting with it. That is, the implementor can put whatever information they find useful in it.

Parameters:x (x-type) – location of derivative
Return type:tuple
getDualProduct(x, r)

returns the dual product of x and r

Return type:float
getGradient(x, *args)

returns the gradient of f at x using the precalculated values for x.

Parameters:
  • x (x-type) – location of derivative
  • args – pre-calculated values for x from getArguments()
Return type:

r-type

getInverseHessianApproximation(x, r, *args)

returns an approximative evaluation p of the inverse of the Hessian operator of the cost function for a given gradient r at a given location x: H(x) p = r

Note

In general it is assumed that the Hessian H(x) needs to be calculate in each call for a new location x. However, the solver may suggest that this is not required, typically when the iteration is close to completeness.

Parameters:
  • x (x-type) – location of Hessian operator to be evaluated.
  • r (r-type) – a given gradient
  • args – pre-calculated values for x from getArguments()
Return type:

x-type

getNorm(x)

returns the norm of x

Return type:float
getValue(x, *args)

returns the value f(x) using the precalculated values for x.

Parameters:x (x-type) – a solution approximation
Return type:float
provides_inverse_Hessian_approximation = False
resetCounters()

resets all statistical counters

updateHessian()

notifies the class that the Hessian operator needs to be updated. This method is called by the solver class.

class esys.downunder.MinimizerBFGS(J=None, m_tol=0.0001, J_tol=None, imax=300)

Bases: esys.downunder.minimizers.AbstractMinimizer

Minimizer that uses the Broyden-Fletcher-Goldfarb-Shanno method.

getCostFunction()

return the cost function to be minimized

Return type:CostFunction
getOptions()
getResult()

Returns the result of the minimization.

logSummary()

Outputs a summary of the completed minimization process to the logger.

run(x)
setCallback(callback)

Sets a callback function to be called after every iteration. The arguments to the function are: (k, x, Jx, g_Jxx), where k is the current iteration, x is the current estimate, Jx=f(x) and g_Jxx=grad f(x).

setCostFunction(J)

set the cost function to be minimized

Parameters:J (CostFunction) – the cost function to be minimized
setMaxIterations(imax)

Sets the maximum number of iterations before the minimizer terminates.

setOptions(**opts)
setTolerance(m_tol=0.0001, J_tol=None)

Sets the tolerance for the stopping criterion. The minimizer stops when an appropriate norm is less than m_tol.

class esys.downunder.MinimizerException

Bases: exceptions.Exception

This is a generic exception thrown by a minimizer.

args
message
class esys.downunder.MinimizerIterationIncurableBreakDown

Bases: esys.downunder.minimizers.MinimizerException

Exception thrown if the iteration scheme encountered an incurable breakdown.

args
message
class esys.downunder.MinimizerLBFGS(J=None, m_tol=0.0001, J_tol=None, imax=300)

Bases: esys.downunder.minimizers.AbstractMinimizer

Minimizer that uses the limited-memory Broyden-Fletcher-Goldfarb-Shanno method.

getCostFunction()

return the cost function to be minimized

Return type:CostFunction
getOptions()
getResult()

Returns the result of the minimization.

logSummary()

Outputs a summary of the completed minimization process to the logger.

run(x)
Parameters:x (Data) – Level set function representing our initial guess
Returns:Level set function representing the solution
Return type:Data
setCallback(callback)

Sets a callback function to be called after every iteration. The arguments to the function are: (k, x, Jx, g_Jxx), where k is the current iteration, x is the current estimate, Jx=f(x) and g_Jxx=grad f(x).

setCostFunction(J)

set the cost function to be minimized

Parameters:J (CostFunction) – the cost function to be minimized
setMaxIterations(imax)

Sets the maximum number of iterations before the minimizer terminates.

setOptions(**opts)
setTolerance(m_tol=0.0001, J_tol=None)

Sets the tolerance for the stopping criterion. The minimizer stops when an appropriate norm is less than m_tol.

class esys.downunder.MinimizerMaxIterReached

Bases: esys.downunder.minimizers.MinimizerException

Exception thrown if the maximum number of iteration steps is reached.

args
message
class esys.downunder.MinimizerNLCG(J=None, m_tol=0.0001, J_tol=None, imax=300)

Bases: esys.downunder.minimizers.AbstractMinimizer

Minimizer that uses the nonlinear conjugate gradient method (Fletcher-Reeves variant).

getCostFunction()

return the cost function to be minimized

Return type:CostFunction
getOptions()

Returns a dictionary of minimizer-specific options.

getResult()

Returns the result of the minimization.

logSummary()

Outputs a summary of the completed minimization process to the logger.

run(x)
setCallback(callback)

Sets a callback function to be called after every iteration. The arguments to the function are: (k, x, Jx, g_Jxx), where k is the current iteration, x is the current estimate, Jx=f(x) and g_Jxx=grad f(x).

setCostFunction(J)

set the cost function to be minimized

Parameters:J (CostFunction) – the cost function to be minimized
setMaxIterations(imax)

Sets the maximum number of iterations before the minimizer terminates.

setOptions(**opts)

Sets minimizer-specific options. For a list of possible options see getOptions().

setTolerance(m_tol=0.0001, J_tol=None)

Sets the tolerance for the stopping criterion. The minimizer stops when an appropriate norm is less than m_tol.

class esys.downunder.NetCdfData(data_type, filename, altitude=0.0, data_variable=None, error=None, scale_factor=None, null_value=None, reference_system=None)

Bases: esys.downunder.datasources.DataSource

Data Source for gridded netCDF data that use CF/COARDS conventions.

ACOUSTIC = 2
GRAVITY = 0
MAGNETIC = 1
getDataExtents()

returns ( (x0, y0), (nx, ny), (dx, dy) )

getDataType()
getHeightScale()

returns the height scale factor to convert from meters to the appropriate units of the reference system used.

Return type:float
getReferenceSystem()

returns the reference coordinate system

Return type:ReferenceSystem
getSubsamplingFactor()

Returns the subsampling factor that was set via setSubsamplingFactor (see there).

getSurveyData(domain, origin, NE, spacing)
getTags()

returns the list of tags

Return type:list
getUtmZone()
hasTag(tag)

returns true if the data set has tag tag

Return type:bool
setSubsamplingFactor(f)

Sets the data subsampling factor (default=1).

The factor is applied in all dimensions. For example a 2D dataset with 300 x 150 data points will be reduced to 150 x 75 when a subsampling factor of 2 is used. This becomes important when adding data of varying resolution to a DomainBuilder.

class esys.downunder.NumpyData(data_type, data, error=1.0, length=1000.0, null_value=-1.0, tags=, []origin=None)

Bases: esys.downunder.datasources.DataSource

ACOUSTIC = 2
GRAVITY = 0
MAGNETIC = 1
getDataExtents()
getDataType()
getHeightScale()

returns the height scale factor to convert from meters to the appropriate units of the reference system used.

Return type:float
getReferenceSystem()

returns the reference coordinate system

Return type:ReferenceSystem
getSubsamplingFactor()

Returns the subsampling factor that was set via setSubsamplingFactor (see there).

getSurveyData(domain, origin, NE, spacing)
getTags()

returns the list of tags

Return type:list
getUtmZone()

returns a dummy UTM zone since this class does not use real coordinate values.

hasTag(tag)

returns true if the data set has tag tag

Return type:bool
setSubsamplingFactor(f)

Sets the data subsampling factor (default=1).

The factor is applied in all dimensions. For example a 2D dataset with 300 x 150 data points will be reduced to 150 x 75 when a subsampling factor of 2 is used. This becomes important when adding data of varying resolution to a DomainBuilder.

class esys.downunder.ReferenceSystem(name='none')

Bases: object

Generic identifier for coordinate systems.

createTransformation(domain)

creates an appropriate coordinate transformation on a given domain

Note

needs to be overwritten by a particular reference system

Parameters:domain (esys.escript.AbstractDomain) – domain of transformation
Return type:SpatialCoordinateTransformation
getName()

returns the name of the reference system

isCartesian()

returns if the reference system is Cartesian

Note

needs to be overwritten by a particular reference system

Return type:bool
isTheSame(other)

test if argument other defines the same reference system

Parameters:other (ReferenceSystem) – a second reference system
Returns:True if other defines the same reference system
Return type:bool

Note

needs to be overwritten by a particular reference system

class esys.downunder.Regularization(domain, numLevelSets=1, w0=None, w1=None, wc=None, location_of_set_m=<esys.escriptcore.escriptcpp.Data object at 0x7fd13d1723c0>, useDiagonalHessianApproximation=False, tol=1e-08, coordinates=None, scale=None, scale_c=None)

Bases: esys.downunder.costfunctions.CostFunction

The regularization term for the level set function m within the cost function J for an inversion:

J(m)=1/2 * sum_k integrate( mu[k] * ( w0[k] * m_k**2 * w1[k,i] * m_{k,i}**2) + sum_l<k mu_c[l,k] wc[l,k] * | curl(m_k) x curl(m_l) |^2

where w0[k], w1[k,i] and wc[k,l] are non-negative weighting factors and mu[k] and mu_c[l,k] are trade-off factors which may be altered during the inversion. The weighting factors are normalized such that their integrals over the domain are constant:

integrate(w0[k] + inner(w1[k,:],1/L[:]**2))=scale[k] volume(domain)* integrate(wc[l,k]*1/L**4)=scale_c[k] volume(domain) *

getArguments(m)
getCoordinateTransformation()

returns the coordinate transformation being used

Return type:CoordinateTransformation
getDomain()

returns the domain of the regularization term

Return type:Domain
getDualProduct(m, r)

returns the dual product of a gradient represented by X=r[1] and Y=r[0] with a level set function m:

Y_i*m_i + X_ij*m_{i,j}
Return type:float
getGradient(m, grad_m)

returns the gradient of the cost function J with respect to m.

Note:This implementation returns Y_k=dPsi/dm_k and X_kj=dPsi/dm_kj
getInverseHessianApproximation(m, r, grad_m)
getNorm(m)

returns the norm of m.

Parameters:m (Data) – level set function
Return type:float
getNumLevelSets()

returns the number of level set functions

Return type:int
getNumTradeOffFactors()

returns the number of trade-off factors being used.

Return type:int
getPDE()

returns the linear PDE to be solved for the Hessian Operator inverse

Return type:LinearPDE
getValue(m, grad_m)

returns the value of the cost function J with respect to m. This equation is specified in the inversion cookbook.

Return type:float
provides_inverse_Hessian_approximation = False
setTradeOffFactors(mu=None)

sets the trade-off factors for the level-set variation and the cross-gradient.

Parameters:mu (list of float or `numpy.array`) – new values for the trade-off factors where values mu[:numLevelSets] are the trade-off factors for the level-set variation and the remaining values for the cross-gradient part with mu_c[l,k]=mu[numLevelSets+l+((k-1)*k)/2] (l<k). If no values for mu are given ones are used. Values must be positive.
setTradeOffFactorsForCrossGradient(mu_c=None)

sets the trade-off factors for the cross-gradient terms.

Parameters:mu_c (float, list of float or numpy.array) – new values for the trade-off factors for the cross-gradient terms. Values must be positive. If no value is given ones are used. Only value mu_c[l,k] for l<k are used.
setTradeOffFactorsForVariation(mu=None)

sets the trade-off factors for the level-set variation part.

Parameters:mu (float, list of float or `numpy.array`) – new values for the trade-off factors. Values must be positive.
updateHessian()

notifies the class to recalculate the Hessian operator.

class esys.downunder.Ricker(f_dom=40, t_dom=None)

Bases: esys.downunder.seismic.Wavelet

The Ricker Wavelet w=f(t)

getAcceleration(t)

get the acceleration f’‘(t) at time t

getCenter()

return value of wavelet center

getTimeScale()

returns the time scale which is the inverse of the largest freqence with a significant spectral component.

getValue(t)

get value of wavelet at time t

class esys.downunder.SeismicSource(x, y=0.0, omega=0.0, elevation=0.0, power=None, orientation=None)

Bases: object

describes a seimic source by location (x,y), frequency omega, power (if known) and orientation (if known). this class is used to tag seismic data sources.

getElevation()

return elevation of source :rtype: float

getFrequency()

return frequency of source :rtype: float

getLocation()

return location of source :rtype: tuple of float

getOrientation()

return power of source orientation at frequency :rtype: vector type object or None

getPower()

return power of source at frequency :rtype: complex or None

class esys.downunder.SelfDemagnetizationModel(domain, w, B, background_magnetic_flux_density, coordinates=None, fixPotentialAtBottom=False, tol=1e-08)

Bases: esys.downunder.forwardmodels.ForwardModelWithPotential

Forward Model for magnetic inversion with self-demagnetization as described in the inversion cookbook.

getArguments(k)

Returns precomputed values shared by getDefect() and getGradient().

Parameters:k (Scalar) – susceptibility
Returns:scalar magnetic potential and corresponding magnetic field
Return type:Scalar, Vector
getCoordinateTransformation()

returns the coordinate transformation being used

Return type:CoordinateTransformation
getDefect(k, phi, grad_phi, magnetic_flux_density)

Returns the value of the defect.

Parameters:
  • k (Scalar) – susceptibility
  • phi (Scalar) – corresponding potential
  • magnetic_flux_density (Vector) – magnetic field
Return type:

float

getDefectGradient(result)
getDomain()

Returns the domain of the forward model.

Return type:Domain
getGradient(k, phi, grad_phi, magnetic_flux_density)

Returns the gradient of the defect with respect to susceptibility.

Parameters:
  • k (Scalar) – susceptibility
  • phi (Scalar) – corresponding potential
  • magnetic_flux_density (Vector) – magnetic field
Return type:

Scalar

getPDE()

Return the underlying PDE.

Return type:LinearPDE
getPotential(k)

Calculates the magnetic potential for a given susceptibility.

Parameters:k (Scalar) – susceptibility
Returns:magnetic potential
Return type:Scalar
getSurvey(index=None)

Returns the pair (data_index, weight_index), where data_i is the data of survey i, weight_i is the weighting factor for survey i. If index is None, all surveys will be returned in a pair of lists.

rescaleWeights(scale=1.0, k_scale=1.0)

rescales the weights such that

sum_s integrate( ( w_i[s] *B_i[s]) (w_j[s]*1/L_j) * L**2 * (background_magnetic_flux_density_j[s]*1/L_j) * k_scale )=scale

Parameters:
  • scale (positive float) – scale of data weighting factors
  • k_scale (Scalar) – scale of susceptibility.
class esys.downunder.SimpleSEGYWriter(receiver_group=None, source=0.0, sampling_interval=0.004, text='some seimic data')

Bases: object

as simple writer for 2D and 3D seimic lines in particular for synthetic data

Typical usage:

from esys.escript import unitsSI as U sw=SimpleSEGYWriter([0.,100*U.m,200*U,m,300.], source=200*U.m, sampling_interval=4*U.msec) while n < 10:

sw.addRecord([i*2., i*0.67, i**2, -i*7])

sw.write(‘example.segy’)

Note:the writer uses obspy
COORDINATE_SCALE = 1000.0
addRecord(record)

adds a record to the traces. a time difference of sample_interval between two records is assumed. The record mast be a list of as many values as given receivers or a float if a single receiver is used.

Parameters:record – list of tracks to be added to the record.
getSamplingInterval()

returns the sampling interval in seconds.

write(filename)

writes to segy file

Parameters:filename – file name
Note:the function uses the obspy module.
class esys.downunder.SmoothAnomaly(lx, ly, lz, x, y, depth, v_inner=None, v_outer=None)

Bases: esys.downunder.datasources.SourceFeature

A source feature in the form of a blob (roughly gaussian).

getMask(x)
getValue(x)
class esys.downunder.SonicHTIWave(domain, v_p, wavelet, source_tag, source_vector=[1.0, 0.0], eps=0.0, delta=0.0, azimuth=0.0, dt=None, p0=None, v0=None, absorption_zone=300.0, absorption_cut=0.01, lumping=True)

Bases: esys.downunder.seismic.WaveBase

Solving the HTI wave equation (along the x_0 axis) with azimuth (rotation around verticle axis) under the assumption of zero shear wave velocities The unknowns are the transversal (along x_0) and vertial stress (Q, P)

Note:In case of a two dimensional domain the second spatial dimenion is depth.
getTimeStepSize()
update(t)

returns the solution for the next time marker t which needs to greater than the time marker from the previous call.

class esys.downunder.SonicWave(domain, v_p, wavelet, source_tag, dt=None, p0=None, p0_t=None, absorption_zone=300.0, absorption_cut=0.01, lumping=True)

Bases: esys.downunder.seismic.WaveBase

Solving the sonic wave equation

p_tt = (v_p**2 * p_i)_i + f(t) * delta_s where (p-) velocity v_p.

f(t) is wavelet acting at a point source term at positon s

getTimeStepSize()
update(t)

returns the solution for the next time marker t which needs to greater than the time marker from the previous call.

class esys.downunder.SpatialCoordinateTransformation(domain, reference=<esys.downunder.coordinates.CartesianReferenceSystem object at 0x7fd149603410>)

Bases: object

Defines an orthogonal coordinate transformation from a domain into the Cartesian domain using a coordinate transformation.

The default implementation is the identity transformation (i.e. no changes are applied to the domain). Overwrite the appropriate methods to define other coordinate systems.

getDomain()

returns the domain of the coordinate transformation.

Return type:esys.escript.AbstractDomain
getReferenceSystem()

returns the reference system used to to define the coordinate transformation

Return type:ReferenceSystem
getScalingFactors()

returns the scaling factors

Return type:esys.escript.Vector
getVolumeFactor()

returns the volume factor for the coordinate transformation

Return type:esys.escript.Scalar
isCartesian()

returns True if the scaling factors (and the volume factor) are equal to 1

Return type:bool
isTheSame(other)

test if argument other defines the same coordinate transformation

Parameters:other (SpatialCoordinateTransformation) – a second coordinate transformation
Returns:True if other defines then same coordinate transformation
Return type:bool
class esys.downunder.StrongJointGravityMagneticInversion(solverclass=None, debug=False, self_demagnetization=False)

Bases: esys.downunder.inversions.InversionDriver

Driver class to perform a joint inversion of Gravity (Bouguer) and magnetic anomaly date with the assumption that there is a functional relationship between density and susceptibility.

The class uses the standard Regularization class for a single level set function,`DensityMapping` and SusceptibilityMapping mappings, the gravity forward model GravityModel and the linear magnetic forward model MagneticModel.

DENSITY = 0
SUSCEPTIBILITY = 1
fixGravityPotentialAtBottom(status=False)

indicates to fix the gravity potential at the bottom to zero (in addition to the top)

Parameters:status (bool) – if True gravity potential at the bottom is set to zero
fixMagneticPotentialAtBottom(status=True)

indicates to fix the magnetic potential at the bottom to zero (in addition to the top)

Parameters:status (bool) – if True magnetic potential at the bottom is set to zero
getCostFunction()

returns the domain of the inversion

Return type:‘InversionCostFunction’
getDomain()

returns the domain of the inversion

Return type:Domain
getLevelSetFunction()

returns the level set function as solution of the optimization problem

Return type:Data
getSolver()

The solver to be used in the inversion process. See the minimizers module for available solvers. By default, the L-BFGS minimizer is used.

Return type:‘AbstractMinimizer’.
isSetUp()

returns True if the inversion is set up and is ready to run.

Return type:bool
run()

This function runs the inversion.

Returns:physical parameters as result of the inversion
Return type:list of physical parameters or a physical parameter
setCostFunction(costfunction)

sets the cost function of the inversion. This function needs to be called before the inversion iteration can be started.

Parameters:costfunction (‘InversionCostFunction’) – domain of the inversion
setInitialGuess(rho=None, k=None)

set the initial guess rho for density and k for susceptibility for the inversion iteration.

Parameters:
  • rho (Scalar) – initial value for the density anomaly.
  • k (Scalar) – initial value for the susceptibility anomaly.
setSolverCallback(callback=None)

Sets the callback function which is called after every solver iteration.

setSolverMaxIterations(maxiter=None)

Sets the maximum number of solver iterations to run. If maxiter is reached the iteration is terminated and MinimizerMaxIterReached is thrown.

Parameters:maxiter (positive int) – maximum number of iteration steps.
setSolverTolerance(m_tol=None, J_tol=None)

Sets the error tolerance for the solver. An acceptable solution is considered to be found once the tolerance is reached.

Parameters:
  • m_tol (float or None) – tolerance for changes to level set function. If None changes to the level set function are not checked for convergence during iteration.
  • J_tol (float or None) – tolerance for changes to cost function. If None changes to the cost function are not checked for convergence during iteration.
Note:

if both arguments are equal to None the default setting m_tol=1e-4, J_tol=None is used.

setup(domainbuilder, rho0=None, drho=None, rho_z0=None, rho_beta=None, k0=None, dk=None, k_z0=None, k_beta=None, w0=None, w1=None, w_gc=None, rho_at_depth=None, k_at_depth=None)

Sets up the inversion from an instance domainbuilder of a DomainBuilder. Gravity and magnetic data attached to the domainbuilder are considered in the inversion. If magnetic data are given as scalar it is assumed that values are collected in direction of the background magnetic field.

Parameters:
  • domainbuilder (DomainBuilder) – Domain builder object with gravity source(s)
  • rho0 (float or Scalar) – reference density, see DensityMapping. If not specified, zero is used.
  • drho (float or Scalar) – density scale, see DensityMapping. If not specified, 2750kg/m^3 is used.
  • rho_z0 (float or Scalar) – reference depth for depth weighting for density, see DensityMapping. If not specified, zero is used.
  • rho_beta (float or Scalar) – exponent for depth weighting for density, see DensityMapping. If not specified, zero is used.
  • k0 (float or Scalar) – reference susceptibility, see SusceptibilityMapping. If not specified, zero is used.
  • dk (float or Scalar) – susceptibility scale, see SusceptibilityMapping. If not specified, 2750kg/m^3 is used.
  • k_z0 (float or Scalar) – reference depth for depth weighting for susceptibility, see SusceptibilityMapping. If not specified, zero is used.
  • k_beta (float or Scalar) – exponent for depth weighting for susceptibility, see SusceptibilityMapping. If not specified, zero is used.
  • w0 (Scalar or float) – weighting factor for level set term regularization. If not set zero is assumed.
  • w1 (Data or ndarray of shape (DIM,)) – weighting factor for the gradient term in the regularization see Regularization. If not set zero is assumed
  • w_gc (Scalar or float) – weighting factor for the cross gradient term in the regularization, see Regularization. If not set one is assumed
  • k_at_depth (float or None) – value for susceptibility at depth, see DomainBuilder.
  • rho_at_depth (float or None) – value for density at depth, see DomainBuilder.
siloWriterCallback(k, m, Jm, g_Jm)

callback function that can be used to track the solution

Parameters:
  • k – iteration count
  • m – current m approximation
  • Jm – value of cost function
  • g_Jm – gradient of f at x
class esys.downunder.SusceptibilityMapping(domain, z0=None, k0=None, dk=None, beta=None)

Bases: esys.downunder.mappings.LinearMapping

Susceptibility mapping with depth weighting

k = k0 + dk * ( (x_2 - z0)/l_z)^(beta/2) ) * m

getDerivative(m)

returns the value for the derivative of the mapping for m

getInverse(p)

returns the value of the inverse of the mapping for s

getTypicalDerivative()

returns a typical value for the derivative

getValue(m)

returns the value of the mapping for m

class esys.downunder.SyntheticData(data_type, n_length=1, n_depth=1, depth_offset=0.0, depth=None, amplitude=None, DIM=2, number_of_elements=10, length=1000.0, B_b=None, data_offset=0, full_knowledge=False, s=0.0)

Bases: esys.downunder.datasources.SyntheticDataBase

Defines synthetic gravity/magnetic data based on harmonic property anomaly

rho = amplitude * sin(n_depth * pi /depth * (z+depth_offset)) * sin(n_length * pi /length * (x - shift) )

for all x and z<=0. For z>0 rho = 0.

ACOUSTIC = 2
GRAVITY = 0
MAGNETIC = 1
getDataExtents()

returns the lateral data extend of the data set

getDataType()

returns the data type

getHeightScale()

returns the height scale factor to convert from meters to the appropriate units of the reference system used.

Return type:float
getReferenceProperty(domain=None)

Returns the reference Data object that was used to generate the gravity/susceptibility anomaly data.

getReferenceSystem()

returns the reference coordinate system

Return type:ReferenceSystem
getSubsamplingFactor()

Returns the subsampling factor that was set via setSubsamplingFactor (see there).

getSurveyData(domain, origin, number_of_elements, spacing)

returns the survey data placed on a given domain.

Parameters:
  • domain (Domain) – domain on which the data are to be placed
  • origin (list of float) – origin of the domain
  • number_of_elements (list of int) – number of elements (or cells) in each spatial direction used to span the domain
  • spacing (list of float) – cell size in each spatial direction
Returns:

observed gravity field or magnetic flux density for each cell in the domain and for each cell an indicator 1/0 if the data are valid or not.

Return type:

pair of Scalar

getTags()

returns the list of tags

Return type:list
getUtmZone()

returns a dummy UTM zone since this class does not use real coordinate values.

hasTag(tag)

returns true if the data set has tag tag

Return type:bool
setSubsamplingFactor(f)

Sets the data subsampling factor (default=1).

The factor is applied in all dimensions. For example a 2D dataset with 300 x 150 data points will be reduced to 150 x 75 when a subsampling factor of 2 is used. This becomes important when adding data of varying resolution to a DomainBuilder.

class esys.downunder.SyntheticDataBase(data_type, DIM=2, number_of_elements=10, length=1000.0, B_b=None, data_offset=0, full_knowledge=False)

Bases: esys.downunder.datasources.DataSource

Base class to define reference data based on a given property distribution (density or susceptibility). Data are collected from a square region of vertical extent length on a grid with number_of_elements cells in each direction.

The synthetic data are constructed by solving the appropriate forward problem. Data can be sampled with an offset from the surface at z=0 or using the entire subsurface region.

ACOUSTIC = 2
GRAVITY = 0
MAGNETIC = 1
getDataExtents()

returns the lateral data extend of the data set

getDataType()

returns the data type

getHeightScale()

returns the height scale factor to convert from meters to the appropriate units of the reference system used.

Return type:float
getReferenceProperty(domain=None)

Returns the reference Data object that was used to generate the gravity/susceptibility anomaly data.

Returns:the density or susceptibility anomaly used to create the survey data
Note:it can be assumed that in the first call the domain argument is present so the actual anomaly data can be created. In subsequent calls this may not be true.
Note:method needs to be overwritten
getReferenceSystem()

returns the reference coordinate system

Return type:ReferenceSystem
getSubsamplingFactor()

Returns the subsampling factor that was set via setSubsamplingFactor (see there).

getSurveyData(domain, origin, number_of_elements, spacing)

returns the survey data placed on a given domain.

Parameters:
  • domain (Domain) – domain on which the data are to be placed
  • origin (list of float) – origin of the domain
  • number_of_elements (list of int) – number of elements (or cells) in each spatial direction used to span the domain
  • spacing (list of float) – cell size in each spatial direction
Returns:

observed gravity field or magnetic flux density for each cell in the domain and for each cell an indicator 1/0 if the data are valid or not.

Return type:

pair of Scalar

getTags()

returns the list of tags

Return type:list
getUtmZone()

returns a dummy UTM zone since this class does not use real coordinate values.

hasTag(tag)

returns true if the data set has tag tag

Return type:bool
setSubsamplingFactor(f)

Sets the data subsampling factor (default=1).

The factor is applied in all dimensions. For example a 2D dataset with 300 x 150 data points will be reduced to 150 x 75 when a subsampling factor of 2 is used. This becomes important when adding data of varying resolution to a DomainBuilder.

class esys.downunder.SyntheticFeatureData(data_type, features, DIM=2, number_of_elements=10, length=1000.0, B_b=None, data_offset=0, full_knowledge=False)

Bases: esys.downunder.datasources.SyntheticDataBase

Uses a list of SourceFeature objects to define synthetic anomaly data.

ACOUSTIC = 2
GRAVITY = 0
MAGNETIC = 1
getDataExtents()

returns the lateral data extend of the data set

getDataType()

returns the data type

getHeightScale()

returns the height scale factor to convert from meters to the appropriate units of the reference system used.

Return type:float
getReferenceProperty(domain=None)

Returns the reference Data object that was used to generate the gravity/susceptibility anomaly data.

getReferenceSystem()

returns the reference coordinate system

Return type:ReferenceSystem
getSubsamplingFactor()

Returns the subsampling factor that was set via setSubsamplingFactor (see there).

getSurveyData(domain, origin, number_of_elements, spacing)

returns the survey data placed on a given domain.

Parameters:
  • domain (Domain) – domain on which the data are to be placed
  • origin (list of float) – origin of the domain
  • number_of_elements (list of int) – number of elements (or cells) in each spatial direction used to span the domain
  • spacing (list of float) – cell size in each spatial direction
Returns:

observed gravity field or magnetic flux density for each cell in the domain and for each cell an indicator 1/0 if the data are valid or not.

Return type:

pair of Scalar

getTags()

returns the list of tags

Return type:list
getUtmZone()

returns a dummy UTM zone since this class does not use real coordinate values.

hasTag(tag)

returns true if the data set has tag tag

Return type:bool
setSubsamplingFactor(f)

Sets the data subsampling factor (default=1).

The factor is applied in all dimensions. For example a 2D dataset with 300 x 150 data points will be reduced to 150 x 75 when a subsampling factor of 2 is used. This becomes important when adding data of varying resolution to a DomainBuilder.

class esys.downunder.TTIWave(domain, v_p, v_s, wavelet, source_tag, source_vector=[0.0, 1.0], eps=0.0, delta=0.0, theta=0.0, rho=1.0, dt=None, u0=None, v0=None, absorption_zone=300.0, absorption_cut=0.01, lumping=True)

Bases: esys.downunder.seismic.WaveBase

Solving the 2D TTI wave equation with

sigma_xx= c11*e_xx + c13*e_zz + c15*e_xz sigma_zz= c13*e_xx + c33*e_zz + c35*e_xz sigma_xz= c15*e_xx + c35*e_zz + c55*e_xz

the coefficicients c11, c13, etc are calculated from the tompsen parameters eps and delta and the tilt theta

Note:currently only the 2D case is supported.
getTimeStepSize()
update(t)

returns the solution for the next time marker t which needs to greater than the time marker from the previous call.

class esys.downunder.VTIWave(domain, v_p, v_s, wavelet, source_tag, source_vector=[0.0, 0.0, 1.0], eps=0.0, gamma=0.0, delta=0.0, rho=1.0, dt=None, u0=None, v0=None, absorption_zone=300.0, absorption_cut=0.01, lumping=True, useFast=False)

Bases: esys.downunder.seismic.WaveBase

Solving the VTI wave equation

Note:In case of a two dimensional domain the second spatial dimenion is depth.
getTimeStepSize()
update(t)

returns the solution for the next time marker t which needs to greater than the time marker from the previous call.

class esys.downunder.WaveBase(dt, u0, v0, t0=0.0)

Bases: object

Base for wave propagation using the Verlet scheme.

u_tt = A(t,u), u(t=t0)=u0, u_t(t=t0)=v0

with a given acceleration force as function of time.

a_n=A(t_{n-1}) v_n=v_{n-1} + dt * a_n u_n=u_{n-1} + dt * v_n

getTimeStepSize()
update(t)

returns the solution for the next time marker t which needs to greater than the time marker from the previous call.

Functions

esys.downunder.CartesianCoordinateTransformation(domain, reference=<esys.downunder.coordinates.CartesianReferenceSystem object at 0x7fd148998450>)
esys.downunder.GRS80ReferenceSystem()

returns the GeodeticReferenceSystem for the GRS80 Ellipsoid eg. used by Geocentric Datum of Australia GDA94

esys.downunder.SphericalReferenceSystem(R=6378137.0)

returns the GeodeticReferenceSystem of a sphere :param R: sphere radius :type R: positive double

esys.downunder.WGS84ReferenceSystem()

returns the GeodeticReferenceSystem for the WGS84 Ellipsoid

esys.downunder.createAbsorbtionLayerFunction(x, absorption_zone=300.0, absorption_cut=0.01, top_absorbation=False)

creating a distribution which is one in the interior of the domain of x and is falling down to the value ‘absorption_cut’ over a margain of thickness ‘absorption_zone’ toward each boundary except the top of the domain.

Parameters:
  • x (Data) – location of points in the domain
  • absorption_zone – thickness of the aborption zone
  • absorption_cut – value of decay function on domain boundary
Returns:

function on ‘x’ which is one in the iterior and decays to almost zero over a margin toward the boundary.

esys.downunder.makeTranformation(domain, coordinates=None)

returns a SpatialCoordinateTransformation for the given domain

Parameters:
Returns:

the spatial coordinate system for the given domain of the specified reference system coordinates. If coordinates is already spatial coordinate system based on the riven domain coordinates is returned. Otherwise an appropriate spatial coordinate system is created.

Return type:

SpatialCoordinateTransformation

Others

  • __builtins__
  • __copyright__
  • __doc__
  • __file__
  • __license__
  • __name__
  • __nodocorecursion
  • __package__
  • __path__
  • __url__