Forward Simulation

Simulation Class

The problem is a partial differential equation of the form:

\[c(m, u) = 0\]

Here, \(m\) is the model and u is the field (or fields). Given the model, \(m\), we can calculate the fields \(u(m)\), however, the data we collect is a subset of the fields, and can be defined by a linear projection, \(P\).

\[d_\text{pred} = P u(m)\]

For the inverse problem, we are interested in how changing the model transforms the data, as such we can take write the Taylor expansion:

\[Pu(m + hv) = Pu(m) + hP\frac{\partial u(m)}{\partial m} v + \mathcal{O}(h^2 \left\| v \right\| )\]

We can linearize and define the sensitivity matrix as:

\[J = P\frac{\partial u}{\partial m}\]

The sensitivity matrix, and it’s transpose will be used in the inverse problem to (locally) find how model parameters change the data, and optimize!

Working with the general PDE, \(c(m, u) = 0\), where m is the model and u is the field, the sensitivity is defined as:

\[J = P\frac{\partial u}{\partial m}\]

We can take the derivative of the PDE:

\[\nabla_m c(m, u) \partial m + \nabla_u c(m, u) \partial u = 0\]

If the forward problem is invertible, then we can rearrange for \(\frac{\partial u}{\partial m}\):

\[J = - P \left( \nabla_u c(m, u) \right)^{-1} \nabla_m c(m, u)\]

This can often be computed given a vector (i.e. \(J(v)\)) rather than stored, as \(J\) is a large dense matrix.

The API

Simulation

class SimPEG.simulation.BaseSimulation(*args, **kwargs)[source]

BaseSimulation is the base class for all geophysical forward simulations in SimPEG.

Required Properties:

  • counter (Counter): A SimPEG.utils.Counter object, an instance of Counter

  • mesh (BaseMesh): a discretize mesh instance, an instance of BaseMesh

  • sensitivity_path (String): path to store the sensitivty, a unicode string, Default: ./sensitivity/

  • None

  • solver_opts (Dictionary): solver options as a kwarg dict, a dictionary

  • survey (BaseSurvey): a survey object, an instance of BaseSurvey

Optional Properties:

  • model (Model): Inversion model., a numpy array of <class ‘float’>, <class ‘int’> with shape (*, *) or (*)

property mesh

mesh (BaseMesh): a discretize mesh instance, an instance of BaseMesh

property survey

survey (BaseSurvey): a survey object, an instance of BaseSurvey

property counter

counter (Counter): A SimPEG.utils.Counter object, an instance of Counter

property sensitivity_path

sensitivity_path (String): path to store the sensitivty, a unicode string, Default: ./sensitivity/

property solver_opts

solver_opts (Dictionary): solver options as a kwarg dict, a dictionary

deleteTheseOnModelUpdate = []
clean_on_model_update = []

List of matrix names to have their factors cleared on a model update

property Solver

Solver has been deprecated. See simulation.solver for documentation

property solverOpts

solverOpts has been deprecated. See solver_opts for documentation

property solver

None

fields(m=None)[source]

u = fields(m) The field given the model. :param numpy.ndarray m: model :rtype: numpy.ndarray :return: u, the fields

dpred(m, f=None)[source]

Create the projected data from a model. The fields, f, (if provided) will be used for the predicted data instead of recalculating the fields (which may be expensive!).

\[d_\text{pred} = P(f(m))\]

Where P is a projection of the fields onto the data space.

Jvec(m, v, f=None)[source]

Jv = Jvec(m, v, f=None) Effect of J(m) on a vector v. :param numpy.ndarray m: model :param numpy.ndarray v: vector to multiply :param Fields f: fields :rtype: numpy.ndarray :return: Jv

Jtvec(m, v, f=None)[source]

Jtv = Jtvec(m, v, f=None) Effect of transpose of J(m) on a vector v. :param numpy.ndarray m: model :param numpy.ndarray v: vector to multiply :param Fields f: fields :rtype: numpy.ndarray :return: JTv

Jvec_approx(m, v, f=None)[source]

Approximate effect of J(m) on a vector v :param numpy.ndarray m: model :param numpy.ndarray v: vector to multiply :param Fields f: fields :rtype: numpy.ndarray :return: approxJv

Jtvec_approx(m, v, f=None)[source]

Approximate effect of transpose of J(m) on a vector v. :param numpy.ndarray m: model :param numpy.ndarray v: vector to multiply :param Fields f: fields :rtype: numpy.ndarray :return: JTv

residual(m, dobs, f=None)[source]

The data residual:

\[\mu_\text{data} = \mathbf{d}_\text{pred} - \mathbf{d}_\text{obs}\]
Parameters
Return type

numpy.ndarray

Returns

data residual

make_synthetic_data(m, relative_error=0.05, noise_floor=0.0, f=None, add_noise=False, **kwargs)[source]

Make synthetic data given a model, and a standard deviation. :param numpy.ndarray m: geophysical model :param numpy.ndarray relative_error: standard deviation :param numpy.ndarray noise_floor: noise floor :param numpy.ndarray f: fields for the given model (if pre-calculated)

pair(survey)[source]

Deprecated pairing method. Please use simulation.survey=survey instead

class SimPEG.simulation.BaseTimeSimulation(*args, **kwargs)[source]

Base class for a time domain simulation

Required Properties:

  • counter (Counter): A SimPEG.utils.Counter object, an instance of Counter

  • mesh (BaseMesh): a discretize mesh instance, an instance of BaseMesh

  • sensitivity_path (String): path to store the sensitivty, a unicode string, Default: ./sensitivity/

  • None

  • solver_opts (Dictionary): solver options as a kwarg dict, a dictionary

  • survey (BaseSurvey): a survey object, an instance of BaseSurvey

  • t0 (Float): Origin of the time discretization, a float, Default: 0.0

  • time_steps (TimeStepArray):

    Sets/gets the time steps for the time domain simulation. You can set as an array of dt’s or as a list of tuples/floats. Tuples must be length two with […, (dt, repeat), …] For example, the following setters are the same:

    sim.time_steps = [(1e-6, 3), 1e-5, (1e-4, 2)]
sim.time_steps = np.r_[1e-6,1e-6,1e-6,1e-5,1e-4,1e-4]

    , an array or list of tuples specifying the mesh tensor of <class ‘float’> with shape (*)

Optional Properties:

  • model (Model): Inversion model., a numpy array of <class ‘float’>, <class ‘int’> with shape (*, *) or (*)

property time_steps

time_steps (TimeStepArray): Sets/gets the time steps for the time domain simulation. You can set as an array of dt’s or as a list of tuples/floats. Tuples must be length two with […, (dt, repeat), …] For example, the following setters are the same:

sim.time_steps = [(1e-6, 3), 1e-5, (1e-4, 2)]
sim.time_steps = np.r_[1e-6,1e-6,1e-6,1e-5,1e-4,1e-4]

, an array or list of tuples specifying the mesh tensor of <class ‘float’> with shape (*)

property t0

t0 (Float): Origin of the time discretization, a float, Default: 0.0

property time_mesh
property nT
property times

Modeling times

property timeSteps

timeSteps has been deprecated. See time_steps for documentation

property timeMesh

time_mesh.timeMesh has been deprecated. See time_mesh for documentation

dpred(m, f=None)[source]

Create the projected data from a model. The fields, f, (if provided) will be used for the predicted data instead of recalculating the fields (which may be expensive!).

\[d_\text{pred} = P(f(m))\]

Where P is a projection of the fields onto the data space.

class SimPEG.simulation.LinearSimulation(*args, **kwargs)[source]

Class for a linear simulation of the form

\[d = Gm\]

where \(d\) is a vector of the data, G is the simulation matrix and \(m\) is the model. Inherit this class to build a linear simulation.

Required Properties:

  • counter (Counter): A SimPEG.utils.Counter object, an instance of Counter

  • mesh (BaseMesh): a discretize mesh instance, an instance of BaseMesh

  • sensitivity_path (String): path to store the sensitivty, a unicode string, Default: ./sensitivity/

  • None

  • solver_opts (Dictionary): solver options as a kwarg dict, a dictionary

  • survey (BaseSurvey): a survey object, an instance of BaseSurvey

Optional Properties:

  • linear_model (PhysicalProperty): The model for a linear problem, a physical property

  • model (Model): Inversion model., a numpy array of <class ‘float’>, <class ‘int’> with shape (*, *) or (*)

  • model_map (Mapping): Mapping of The model for a linear problem to the inversion model., a SimPEG Map

Other Properties:

  • model_deriv (Derivative): Derivative of The model for a linear problem wrt the model.

property linear_model

The model for a linear problem

property model_map

Mapping of The model for a linear problem to the inversion model.

property model_deriv

Derivative of The model for a linear problem wrt the model.

property mesh

mesh (BaseMesh): a discretize mesh instance, an instance of BaseMesh

property survey

survey (BaseSurvey): a survey object, an instance of BaseSurvey

property G
fields(m)[source]

u = fields(m) The field given the model. :param numpy.ndarray m: model :rtype: numpy.ndarray :return: u, the fields

dpred(m, f=None)[source]

Create the projected data from a model. The fields, f, (if provided) will be used for the predicted data instead of recalculating the fields (which may be expensive!).

\[d_\text{pred} = P(f(m))\]

Where P is a projection of the fields onto the data space.

getJ(m, f=None)[source]
Jvec(m, v, f=None)[source]

Jv = Jvec(m, v, f=None) Effect of J(m) on a vector v. :param numpy.ndarray m: model :param numpy.ndarray v: vector to multiply :param Fields f: fields :rtype: numpy.ndarray :return: Jv

Jtvec(m, v, f=None)[source]

Jtv = Jtvec(m, v, f=None) Effect of transpose of J(m) on a vector v. :param numpy.ndarray m: model :param numpy.ndarray v: vector to multiply :param Fields f: fields :rtype: numpy.ndarray :return: JTv

class SimPEG.simulation.TimeStepArray(*args, **kwargs)[source]
class_info = 'an array or list of tuples specifying the mesh tensor'
validate(instance, value)[source]

Determine if array is valid based on shape and dtype

Fields

class SimPEG.fields.Fields(*args, **kwargs)[source]
Fancy Field Storage

Required Properties:

  • aliasFields (Dictionary):

    a dictionary of the aliased fields with [alias, location, function], e.g. {“b”:[“e”,”F”,lambda(F,e,ind)]} , a dictionary

  • knownFields (Dictionary):

    a dictionary with the names of the know fields and their location on a mesh e.g. {“e”: “E”, “phi”: “CC”} , a dictionary

  • simulation (BaseSimulation): a SimPEG simulation, an instance of BaseSimulation

property knownFields

knownFields (Dictionary): a dictionary with the names of the know fields and their location on a mesh e.g. {“e”: “E”, “phi”: “CC”} , a dictionary

property aliasFields

aliasFields (Dictionary): a dictionary of the aliased fields with [alias, location, function], e.g. {“b”:[“e”,”F”,lambda(F,e,ind)]} , a dictionary

dtype

alias of builtins.float

property simulation

simulation (BaseSimulation): a SimPEG simulation, an instance of BaseSimulation

property mesh
property survey
startup()[source]
property approxSize

The approximate cost to storing all of the known fields.

class SimPEG.fields.TimeFields(*args, **kwargs)[source]
Fancy Field Storage for time domain problems
fields = TimeFields(simulation=simulation, knownFields={'phi':'CC'})
fields[:,'phi', timeInd] = phi
print(fields[src0,'phi'])

Required Properties:

  • aliasFields (Dictionary):

    a dictionary of the aliased fields with [alias, location, function], e.g. {“b”:[“e”,”F”,lambda(F,e,ind)]} , a dictionary

  • knownFields (Dictionary):

    a dictionary with the names of the know fields and their location on a mesh e.g. {“e”: “E”, “phi”: “CC”} , a dictionary

  • simulation (BaseTimeSimulation): a SimPEG time simulation, an instance of BaseTimeSimulation

property simulation

simulation (BaseTimeSimulation): a SimPEG time simulation, an instance of BaseTimeSimulation

Survey

class SimPEG.survey.BaseSurvey(*args, **kwargs)[source]

Survey holds the sources and receivers for a survey

Required Properties:

  • counter (Counter): A SimPEG counter object, an instance of Counter

  • source_list (a list of BaseSrc): A list of sources for the survey, a list (each item is an instance of BaseSrc)

property counter

counter (Counter): A SimPEG counter object, an instance of Counter

property source_list

source_list (a list of BaseSrc): A list of sources for the survey, a list (each item is an instance of BaseSrc)

getSourceIndex(sources)[source]
property nD

Number of data

property vnD

Vector number of data

property nSrc

Number of Sources

property srcList

srcList has been deprecated. See source_list for documentation

dpred(m=None, f=None)[source]
makeSyntheticData(m, std=None, f=None, force=False, **kwargs)[source]
pair(simulation)[source]
class SimPEG.survey.BaseTimeSurvey(*args, **kwargs)[source]

Required Properties:

  • counter (Counter): A SimPEG counter object, an instance of Counter

  • source_list (a list of BaseSrc): A list of sources for the survey, a list (each item is an instance of BaseSrc)

property unique_times
property times

unique_times.times has been deprecated. See unique_times for documentation

Data

class SimPEG.data.Data(*args, **kwargs)[source]

Data storage. This class keeps track of observed data, relative error of those data and the noise floor.

data = Data(survey, dobs=dobs, relative_error=relative, noise_floor=floor)

or

data = Data(survey, dobs=dobs, standard_deviation=standard_deviation)

Required Properties:

  • dobs (Array):

    Vector of the observed data. The data can be set using the survey parameters:

    data = Data(survey)
for src in survey.source_list:
    for rx in src.receiver_list:
        data[src, rx] = datum

    , a list or numpy array of <class ‘float’>, <class ‘int’> with shape (*)

  • noise_floor (UncertaintyArray):

    Noise floor of the data. This can be set using an array of the same size as the data (e.g. if you want to assign a different noise floor to each datum) or as a scalar if you would like to assign a the same noise floor to all data.

    The standard_deviation is constructed as follows:

    relative_error * np.abs(dobs) + noise_floor

    For example, if you set

    data = Data(survey, dobs=dobs)
data.noise_floor = 1e-10

    then the contribution to the standard_deviation is equal to

    data.noise_floor

    , An array that can be set by a scalar value or numpy array of <class ‘float’>, <class ‘int’> with shape (*)

  • relative_error (UncertaintyArray):

    Relative error of the data. This can be set using an array of the same size as the data (e.g. if you want to assign a different relative error to each datum) or as a scalar if you would like to assign a the same relative error to all data.

    The standard_deviation is constructed as follows:

    relative_error * np.abs(dobs) + noise_floor

    For example, if you set

    data = Data(survey, dobs=dobs)
data.relative_error = 0.05

    then the contribution to the standard_deviation is equal to

    data.relative_error * np.abs(data.dobs)

    , An array that can be set by a scalar value or numpy array of <class ‘float’>, <class ‘int’> with shape (*)

  • survey (BaseSurvey): a SimPEG survey object, an instance of BaseSurvey

property survey

survey (BaseSurvey): a SimPEG survey object, an instance of BaseSurvey

property dobs

dobs (Array): Vector of the observed data. The data can be set using the survey parameters:

data = Data(survey)
for src in survey.source_list:
    for rx in src.receiver_list:
        data[src, rx] = datum

, a list or numpy array of <class ‘float’>, <class ‘int’> with shape (*)

property relative_error

relative_error (UncertaintyArray): Relative error of the data. This can be set using an array of the same size as the data (e.g. if you want to assign a different relative error to each datum) or as a scalar if you would like to assign a the same relative error to all data.

The standard_deviation is constructed as follows:

relative_error * np.abs(dobs) + noise_floor

For example, if you set

data = Data(survey, dobs=dobs)
data.relative_error = 0.05

then the contribution to the standard_deviation is equal to

data.relative_error * np.abs(data.dobs)

, An array that can be set by a scalar value or numpy array of <class ‘float’>, <class ‘int’> with shape (*)

property noise_floor

noise_floor (UncertaintyArray): Noise floor of the data. This can be set using an array of the same size as the data (e.g. if you want to assign a different noise floor to each datum) or as a scalar if you would like to assign a the same noise floor to all data.

The standard_deviation is constructed as follows:

relative_error * np.abs(dobs) + noise_floor

For example, if you set

data = Data(survey, dobs=dobs)
data.noise_floor = 1e-10

then the contribution to the standard_deviation is equal to

data.noise_floor

, An array that can be set by a scalar value or numpy array of <class ‘float’>, <class ‘int’> with shape (*)

property standard_deviation

Data standard deviations. If a relative error and noise floor are provided, the standard_deviation is

data.standard_deviation = (
    data.relative_error*np.abs(data.dobs) +
    data.noise_floor
)

otherwise, the standard_deviation can be set directly

data.standard_deviation = 0.05 * np.absolute(self.dobs) + 1e-12

Note that setting the standard_deviation directly will clear the relative_error and set the value to the noise_floor property.

property nD
property shape
property index_dictionary

Dictionary of data indices by sources and receivers. To set data using survey parameters:

data = Data(survey)
for src in survey.source_list:
    for rx in src.receiver_list:
        index = data.index_dictionary[src][rx]
        data.dobs[index] = datum
tovec()[source]
fromvec(v)[source]
property std

std has been deprecated. See relative_error for documentation

property eps

eps has been deprecated. See noise_floor for documentation

class SimPEG.data.SyntheticData(*args, **kwargs)[source]

Data class for synthetic data. It keeps track of observed and clean data

Required Properties:

  • dclean (Array):

    Vector of the clean synthetic data. The data can be set using the survey parameters:

    data = Data(survey)
for src in survey.source_list:
    for rx in src.receiver_list:
        index = data.index_dictionary(src, rx)
        data.dclean[indices] = datum

    , a list or numpy array of <class ‘float’>, <class ‘int’> with shape (*)

  • dobs (Array):

    Vector of the observed data. The data can be set using the survey parameters:

    data = Data(survey)
for src in survey.source_list:
    for rx in src.receiver_list:
        data[src, rx] = datum

    , a list or numpy array of <class ‘float’>, <class ‘int’> with shape (*)

  • noise_floor (UncertaintyArray):

    Noise floor of the data. This can be set using an array of the same size as the data (e.g. if you want to assign a different noise floor to each datum) or as a scalar if you would like to assign a the same noise floor to all data.

    The standard_deviation is constructed as follows:

    relative_error * np.abs(dobs) + noise_floor

    For example, if you set

    data = Data(survey, dobs=dobs)
data.noise_floor = 1e-10

    then the contribution to the standard_deviation is equal to

    data.noise_floor

    , An array that can be set by a scalar value or numpy array of <class ‘float’>, <class ‘int’> with shape (*)

  • relative_error (UncertaintyArray):

    Relative error of the data. This can be set using an array of the same size as the data (e.g. if you want to assign a different relative error to each datum) or as a scalar if you would like to assign a the same relative error to all data.

    The standard_deviation is constructed as follows:

    relative_error * np.abs(dobs) + noise_floor

    For example, if you set

    data = Data(survey, dobs=dobs)
data.relative_error = 0.05

    then the contribution to the standard_deviation is equal to

    data.relative_error * np.abs(data.dobs)

    , An array that can be set by a scalar value or numpy array of <class ‘float’>, <class ‘int’> with shape (*)

  • survey (BaseSurvey): a SimPEG survey object, an instance of BaseSurvey

property dclean

dclean (Array): Vector of the clean synthetic data. The data can be set using the survey parameters:

data = Data(survey)
for src in survey.source_list:
    for rx in src.receiver_list:
        index = data.index_dictionary(src, rx)
        data.dclean[indices] = datum

, a list or numpy array of <class ‘float’>, <class ‘int’> with shape (*)

class SimPEG.data.UncertaintyArray(*args, **kwargs)[source]
class_info = 'An array that can be set by a scalar value or numpy array'
validate(instance, value)[source]

Determine if array is valid based on shape and dtype

Sources

class SimPEG.survey.BaseSrc(*args, **kwargs)[source]

SimPEG Source Object

Required Properties:

  • receiver_list (a list of BaseRx): receiver list, a list (each item is an instance of BaseRx)

Optional Properties:

  • location (SourceLocationArray): Location of the source [x, y, z] in 3D, a 1D array denoting the source location of <class ‘float’>, <class ‘int’> with shape (*)

property loc

loc has been deprecated. See location for documentation

property rxList

rxList has been deprecated. See receiver_list for documentation

getReceiverIndex(receiver)[source]
property nD

Number of data

property vnD

Vector number of data

property receiver_list

receiver_list (a list of BaseRx): receiver list, a list (each item is an instance of BaseRx)

property location

location (SourceLocationArray): Location of the source [x, y, z] in 3D, a 1D array denoting the source location of <class ‘float’>, <class ‘int’> with shape (*)

class SimPEG.survey.SourceLocationArray(*args, **kwargs)[source]
class_info = 'a 1D array denoting the source location'
validate(instance, value)[source]

Determine if array is valid based on shape and dtype

Receivers

class SimPEG.survey.BaseRx(*args, **kwargs)[source]

SimPEG Receiver Object

Required Properties:

  • locations (RxLocationArray): Locations of the receivers (nRx x nDim), an array of receiver locations of <class ‘float’>, <class ‘int’> with shape (*, *)

  • projGLoc (StringChoice): Projection grid location, default is CC, any of “CC”, “Fx”, “Fy”, “Fz”, “Ex”, “Ey”, “Ez”, “N”, Default: CC

  • storeProjections (Boolean): Store calls to getP (organized by mesh), a boolean, Default: True

property projGLoc

projGLoc (StringChoice): Projection grid location, default is CC, any of “CC”, “Fx”, “Fy”, “Fz”, “Ex”, “Ey”, “Ez”, “N”, Default: CC

property storeProjections

storeProjections (Boolean): Store calls to getP (organized by mesh), a boolean, Default: True

property locations

locations (RxLocationArray): Locations of the receivers (nRx x nDim), an array of receiver locations of <class ‘float’>, <class ‘int’> with shape (*, *)

property locs

locs has been deprecated. See locations for documentation

property nD

Number of data in the receiver.

getP(mesh, projGLoc=None)[source]

Returns the projection matrices as a list for all components collected by the receivers.

Note

Projection matrices are stored as a dictionary listed by meshes.

eval(**kwargs)[source]
evalDeriv(**kwargs)[source]
class SimPEG.survey.BaseTimeRx(*args, **kwargs)[source]

SimPEG Receiver Object for time-domain simulations

Required Properties:

  • locations (RxLocationArray): Locations of the receivers (nRx x nDim), an array of receiver locations of <class ‘float’>, <class ‘int’> with shape (*, *)

  • projGLoc (StringChoice): Projection grid location, default is CC, any of “CC”, “Fx”, “Fy”, “Fz”, “Ex”, “Ey”, “Ez”, “N”, Default: CC

  • projTLoc (StringChoice): location on the time mesh where the data are projected from, either “N” or “CC”, Default: N

  • storeProjections (Boolean): Store calls to getP (organized by mesh), a boolean, Default: True

  • times (Array): times where the recievers measure data, a list or numpy array of <class ‘float’>, <class ‘int’> with shape (*)

property projTLoc

projTLoc (StringChoice): location on the time mesh where the data are projected from, either “N” or “CC”, Default: N

property times

times (Array): times where the recievers measure data, a list or numpy array of <class ‘float’>, <class ‘int’> with shape (*)

property nD

Number of data in the receiver.

getSpatialP(mesh)[source]

Returns the spatial projection matrix.

Note

This is not stored in memory, but is created on demand.

getTimeP(timeMesh)[source]

Returns the time projection matrix.

Note

This is not stored in memory, but is created on demand.

getP(mesh, timeMesh)[source]

Returns the projection matrices as a list for all components collected by the receivers.

Note

Projection matrices are stored as a dictionary (mesh, timeMesh) if storeProjections is True

class SimPEG.survey.RxLocationArray(*args, **kwargs)[source]
class_info = 'an array of receiver locations'
validate(instance, value)[source]

Determine if array is valid based on shape and dtype