esys.downunder.dcresistivityforwardmodeling Package

Classes

class esys.downunder.dcresistivityforwardmodeling.DcResistivityForward

This class allows for the solution of dc resistivity forward problems via the calculation of a primary and secondary potential. Conductivity values are to be provided for the primary problem which is a homogeneous half space of a chosen conductivity and for the secondary problem which typically varies it conductivity spatially across the domain. The primary potential acts as a reference point typically based of some know reference conductivity however any value will suffice. The primary potential is implemented to avoid the use of dirac delta functions.

__init__()

This is a skeleton class for all the other forward modeling classes.

checkBounds()
getApparentResistivity()
getElectrodes()

retuns the list of electrodes with locations

getPotential()
class esys.downunder.dcresistivityforwardmodeling.DipoleDipoleSurvey(domain, primaryConductivity, secondaryConductivity, current, a, n, midPoint, directionVector, numElectrodes)

DipoleDipoleSurvey forward modeling

__init__(domain, primaryConductivity, secondaryConductivity, current, a, n, midPoint, directionVector, numElectrodes)

This is a skeleton class for all the other forward modeling classes.

getApparentResistivityPrimary()
getApparentResistivitySecondary()
getApparentResistivityTotal()
getPotential()

Returns 3 list each made up of a number of list containing primary, secondary and total potentials diferences. Each of the lists contain a list for each value of n.

class esys.downunder.dcresistivityforwardmodeling.LinearPDE(domain, numEquations=None, numSolutions=None, isComplex=False, debug=False)

This class is used to define a general linear, steady, second order PDE for an unknown function u on a given domain defined through a Domain object.

For a single PDE having a solution with a single component the linear PDE is defined in the following form:

-(grad(A[j,l]+A_reduced[j,l])*grad(u)[l]+(B[j]+B_reduced[j])u)[j]+(C[l]+C_reduced[l])*grad(u)[l]+(D+D_reduced)=-grad(X+X_reduced)[j,j]+(Y+Y_reduced)

where grad(F) denotes the spatial derivative of F. Einstein’s summation convention, ie. summation over indexes appearing twice in a term of a sum performed, is used. The coefficients A, B, C, D, X and Y have to be specified through Data objects in Function and the coefficients A_reduced, B_reduced, C_reduced, D_reduced, X_reduced and Y_reduced have to be specified through Data objects in ReducedFunction. It is also allowed to use objects that can be converted into such Data objects. A and A_reduced are rank two, B, C, X, B_reduced, C_reduced and X_reduced are rank one and D, D_reduced, Y and Y_reduced are scalar.

The following natural boundary conditions are considered:

n[j]*((A[i,j]+A_reduced[i,j])*grad(u)[l]+(B+B_reduced)[j]*u)+(d+d_reduced)*u=n[j]*(X[j]+X_reduced[j])+y

where n is the outer normal field. Notice that the coefficients A, A_reduced, B, B_reduced, X and X_reduced are defined in the PDE. The coefficients d and y are each a scalar in FunctionOnBoundary and the coefficients d_reduced and y_reduced are each a scalar in ReducedFunctionOnBoundary.

Constraints for the solution prescribe the value of the solution at certain locations in the domain. They have the form

u=r where q>0

r and q are each scalar where q is the characteristic function defining where the constraint is applied. The constraints override any other condition set by the PDE or the boundary condition.

The PDE is symmetrical if

A[i,j]=A[j,i] and B[j]=C[j] and A_reduced[i,j]=A_reduced[j,i] and B_reduced[j]=C_reduced[j]

For a system of PDEs and a solution with several components the PDE has the form

-grad((A[i,j,k,l]+A_reduced[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k])[j]+(C[i,k,l]+C_reduced[i,k,l])*grad(u[k])[l]+(D[i,k]+D_reduced[i,k]*u[k] =-grad(X[i,j]+X_reduced[i,j])[j]+Y[i]+Y_reduced[i]

A and A_reduced are of rank four, B, B_reduced, C and C_reduced are each of rank three, D, D_reduced, X_reduced and X are each of rank two and Y and Y_reduced are of rank one. The natural boundary conditions take the form:

n[j]*((A[i,j,k,l]+A_reduced[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k])+(d[i,k]+d_reduced[i,k])*u[k]=n[j]*(X[i,j]+X_reduced[i,j])+y[i]+y_reduced[i]

The coefficient d is of rank two and y is of rank one both in FunctionOnBoundary. The coefficients d_reduced is of rank two and y_reduced is of rank one both in ReducedFunctionOnBoundary.

Constraints take the form

u[i]=r[i] where q[i]>0

r and q are each rank one. Notice that at some locations not necessarily all components must have a constraint.

The system of PDEs is symmetrical if

  • A[i,j,k,l]=A[k,l,i,j]

  • A_reduced[i,j,k,l]=A_reduced[k,l,i,j]

  • B[i,j,k]=C[k,i,j]

  • B_reduced[i,j,k]=C_reduced[k,i,j]

  • D[i,k]=D[i,k]

  • D_reduced[i,k]=D_reduced[i,k]

  • d[i,k]=d[k,i]

  • d_reduced[i,k]=d_reduced[k,i]

LinearPDE also supports solution discontinuities over a contact region in the domain. To specify the conditions across the discontinuity we are using the generalised flux J which, in the case of a system of PDEs and several components of the solution, is defined as

J[i,j]=(A[i,j,k,l]+A_reduced[[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k]-X[i,j]-X_reduced[i,j]

For the case of single solution component and single PDE J is defined as

J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[j]+(B[i]+B_reduced[i])*u-X[i]-X_reduced[i]

In the context of discontinuities n denotes the normal on the discontinuity pointing from side 0 towards side 1 calculated from FunctionSpace.getNormal of FunctionOnContactZero. For a system of PDEs the contact condition takes the form

n[j]*J0[i,j]=n[j]*J1[i,j]=(y_contact[i]+y_contact_reduced[i])- (d_contact[i,k]+d_contact_reduced[i,k])*jump(u)[k]

where J0 and J1 are the fluxes on side 0 and side 1 of the discontinuity, respectively. jump(u), which is the difference of the solution at side 1 and at side 0, denotes the jump of u across discontinuity along the normal calculated by jump. The coefficient d_contact is of rank two and y_contact is of rank one both in FunctionOnContactZero or FunctionOnContactOne. The coefficient d_contact_reduced is of rank two and y_contact_reduced is of rank one both in ReducedFunctionOnContactZero or ReducedFunctionOnContactOne. In case of a single PDE and a single component solution the contact condition takes the form

n[j]*J0_{j}=n[j]*J1_{j}=(y_contact+y_contact_reduced)-(d_contact+y_contact_reduced)*jump(u)

In this case the coefficient d_contact and y_contact are each scalar both in FunctionOnContactZero or FunctionOnContactOne and the coefficient d_contact_reduced and y_contact_reduced are each scalar both in ReducedFunctionOnContactZero or ReducedFunctionOnContactOne.

Typical usage:

p = LinearPDE(dom)
p.setValue(A=kronecker(dom), D=1, Y=0.5)
u = p.getSolution()
__init__(domain, numEquations=None, numSolutions=None, isComplex=False, debug=False)

Initializes a new linear PDE.

Parameters:
  • domain (Domain) – domain of the PDE

  • numEquations – number of equations. If None the number of equations is extracted from the PDE coefficients.

  • numSolutions – number of solution components. If None the number of solution components is extracted from the PDE coefficients.

  • debug – if True debug information is printed

checkSymmetry(verbose=True)

Tests the PDE for symmetry.

Parameters:

verbose (bool) – if set to True or not present a report on coefficients which break the symmetry is printed.

Returns:

True if the PDE is symmetric

Return type:

bool

Note:

This is a very expensive operation. It should be used for degugging only! The symmetry flag is not altered.

createOperator()

Returns an instance of a new operator.

getFlux(u=None)

Returns the flux J for a given u.

J[i,j]=(A[i,j,k,l]+A_reduced[A[i,j,k,l]]*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])u[k]-X[i,j]-X_reduced[i,j]

or

J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[l]+(B[j]+B_reduced[j])u-X[j]-X_reduced[j]

Parameters:

u (Data or None) – argument in the flux. If u is not present or equals None the current solution is used.

Returns:

flux

Return type:

Data

getRequiredOperatorType()

Returns the system type which needs to be used by the current set up.

getResidual(u=None)

Returns the residual of u or the current solution if u is not present.

Parameters:

u (Data or None) – argument in the residual calculation. It must be representable in self.getFunctionSpaceForSolution(). If u is not present or equals None the current solution is used.

Returns:

residual of u

Return type:

Data

getSolution()

Returns the solution of the PDE.

Returns:

the solution

Return type:

Data

getSystem()

Returns the operator and right hand side of the PDE.

Returns:

the discrete version of the PDE

Return type:

tuple of Operator and Data

insertConstraint(rhs_only=False)

Applies the constraints defined by q and r to the PDE.

Parameters:

rhs_only (bool) – if True only the right hand side is altered by the constraint

setValue(**coefficients)

Sets new values to coefficients.

Parameters:
  • coefficients – new values assigned to coefficients

  • A (any type that can be cast to a Data object on Function) – value for coefficient A

  • A_reduced (any type that can be cast to a Data object on ReducedFunction) – value for coefficient A_reduced

  • B (any type that can be cast to a Data object on Function) – value for coefficient B

  • B_reduced (any type that can be cast to a Data object on ReducedFunction) – value for coefficient B_reduced

  • C (any type that can be cast to a Data object on Function) – value for coefficient C

  • C_reduced (any type that can be cast to a Data object on ReducedFunction) – value for coefficient C_reduced

  • D (any type that can be cast to a Data object on Function) – value for coefficient D

  • D_reduced (any type that can be cast to a Data object on ReducedFunction) – value for coefficient D_reduced

  • X (any type that can be cast to a Data object on Function) – value for coefficient X

  • X_reduced (any type that can be cast to a Data object on ReducedFunction) – value for coefficient X_reduced

  • Y (any type that can be cast to a Data object on Function) – value for coefficient Y

  • Y_reduced (any type that can be cast to a Data object on ReducedFunction) – value for coefficient Y_reduced

  • d (any type that can be cast to a Data object on FunctionOnBoundary) – value for coefficient d

  • d_reduced (any type that can be cast to a Data object on ReducedFunctionOnBoundary) – value for coefficient d_reduced

  • y (any type that can be cast to a Data object on FunctionOnBoundary) – value for coefficient y

  • d_contact (any type that can be cast to a Data object on FunctionOnContactOne or FunctionOnContactZero) – value for coefficient d_contact

  • d_contact_reduced (any type that can be cast to a Data object on ReducedFunctionOnContactOne or ReducedFunctionOnContactZero) – value for coefficient d_contact_reduced

  • y_contact (any type that can be cast to a Data object on FunctionOnContactOne or FunctionOnContactZero) – value for coefficient y_contact

  • y_contact_reduced (any type that can be cast to a Data object on ReducedFunctionOnContactOne or ReducedFunctionOnContactZero) – value for coefficient y_contact_reduced

  • d_dirac (any type that can be cast to a Data object on DiracDeltaFunctions) – value for coefficient d_dirac

  • y_dirac (any type that can be cast to a Data object on DiracDeltaFunctions) – value for coefficient y_dirac

  • r (any type that can be cast to a Data object on Solution or ReducedSolution depending on whether reduced order is used for the solution) – values prescribed to the solution at the locations of constraints

  • q (any type that can be cast to a Data object on Solution or ReducedSolution depending on whether reduced order is used for the representation of the equation) – mask for location of constraints

Raises:

IllegalCoefficient – if an unknown coefficient keyword is used

class esys.downunder.dcresistivityforwardmodeling.Locator(where, x=array([0., 0., 0.]))

Locator provides access to the values of data objects at a given spatial coordinate x.

In fact, a Locator object finds the sample in the set of samples of a given function space or domain which is closest to the given point x.

__init__(where, x=array([0., 0., 0.]))

Initializes a Locator to access values in Data objects on the Doamin or FunctionSpace for the sample point which is closest to the given point x.

Parameters:
  • where (escript.FunctionSpace) – function space

  • x (numpy.ndarray or list of numpy.ndarray) – location(s) of the Locator

getFunctionSpace()

Returns the function space of the Locator.

getId(item=None)

Returns the identifier of the location.

getValue(data)

Returns the value of data at the Locator if data is a Data object otherwise the object is returned.

getX()

Returns the exact coordinates of the Locator.

setValue(data, v)

Sets the value of the data at the Locator.

class esys.downunder.dcresistivityforwardmodeling.PoleDipoleSurvey(domain, primaryConductivity, secondaryConductivity, current, a, n, midPoint, directionVector, numElectrodes)

Forward model class for a poledipole setup

__init__(domain, primaryConductivity, secondaryConductivity, current, a, n, midPoint, directionVector, numElectrodes)
Parameters:
  • domain (Domain) – domain of the model

  • primaryConductivity (data) – preset primary conductivity data object

  • secondaryConductivity (data) – preset secondary conductivity data object

  • current (float or int) – amount of current to be injected at the current electrode

  • a (list) – the spacing between current and potential electrodes

  • n (float or int) – multiple of spacing between electrodes. typicaly interger

  • midPoint – midPoint of the survey, as a list containing x,y coords

  • directionVector – two element list specifying the direction the survey should extend from the midpoint

  • numElectrodes (int) – the number of electrodes to be used in the survey must be a multiple of 2 for polepole survey:

getApparentResistivityPrimary()
getApparentResistivitySecondary()
getApparentResistivityTotal()
getPotential()

Returns 3 list each made up of a number of list containing primary, secondary and total potentials diferences. Each of the lists contain a list for each value of n.

class esys.downunder.dcresistivityforwardmodeling.PolePoleSurvey(domain, primaryConductivity, secondaryConductivity, current, a, midPoint, directionVector, numElectrodes)

Forward model class for a polepole setup

__init__(domain, primaryConductivity, secondaryConductivity, current, a, midPoint, directionVector, numElectrodes)
Parameters:
  • domain (Domain) – domain of the model

  • primaryConductivity (data) – preset primary conductivity data object

  • secondaryConductivity (data) – preset secondary conductivity data object

  • current (float or int) – amount of current to be injected at the current electrode

  • a (list) – the spacing between current and potential electrodes

  • midPoint – midPoint of the survey, as a list containing x,y coords

  • directionVector – two element list specifying the direction the survey should extend from the midpoint

  • numElectrodes (int) – the number of electrodes to be used in the survey must be a multiple of 2 for polepole survey:

getApparentResistivityPrimary()
getApparentResistivitySecondary()
getApparentResistivityTotal()
getPotential()

returns a list containing 3 lists one for each the primary, secondary and total potential.

class esys.downunder.dcresistivityforwardmodeling.SchlumbergerSurvey(domain, primaryConductivity, secondaryConductivity, current, a, n, midPoint, directionVector, numElectrodes)

Schlumberger survey forward calculation

__init__(domain, primaryConductivity, secondaryConductivity, current, a, n, midPoint, directionVector, numElectrodes)

This is a skeleton class for all the other forward modeling classes.

getApparentResistivity(delPhiList)
getElectrodeDict()

retuns the electrode dictionary

getPotentialAnalytic()

Returns 3 list each made up of a number of list containing primary, secondary and total potentials diferences. Each of the lists contain a list for each value of n.

getPotentialNumeric()

Returns 3 list each made up of a number of list containing primary, secondary and total potentials diferences. Each of the lists contain a list for each value of n.

getSourcesSamples()

return a list of tuples of sample locations followed by a list of tuples of source locations.

class esys.downunder.dcresistivityforwardmodeling.WennerSurvey(domain, primaryConductivity, secondaryConductivity, current, a, midPoint, directionVector, numElectrodes)

WennerSurvey forward calculation

__init__(domain, primaryConductivity, secondaryConductivity, current, a, midPoint, directionVector, numElectrodes)
Parameters:
  • domain (Domain) – domain of the model

  • primaryConductivity (data) – preset primary conductivity data object

  • secondaryConductivity (data) – preset secondary conductivity data object

  • current (float or int) – amount of current to be injected at the current electrode

  • a (list) – the spacing between current and potential electrodes

  • midPoint – midPoint of the survey, as a list containing x,y coords

  • directionVector – two element list specifying the direction the survey should extend from the midpoint

  • numElectrodes (int) – the number of electrodes to be used in the survey must be a multiple of 2 for polepole survey

getApparentResistivityPrimary()
getApparentResistivitySecondary()
getApparentResistivityTotal()
getPotential()

returns a list containing 3 lists one for each the primary, secondary and total potential.

esys.downunder.dcresistivityforwardmodeling.xrange

alias of range

Functions

esys.downunder.dcresistivityforwardmodeling.saveSilo(filename, domain=None, write_meshdata=False, time=0.0, cycle=0, **data)

Writes Data objects and their mesh to a file using the SILO file format.

Example:

temp=Scalar(..)
v=Vector(..)
saveSilo("solution.silo", temperature=temp, velocity=v)

temp and v are written to “solution.silo” where temp is named “temperature” and v is named “velocity”.

Parameters:
  • filename (str) – name of the output file (‘.silo’ is added if required)

  • domain (escript.Domain) – domain of the Data objects. If not specified, the domain of the given Data objects is used.

  • write_meshdata (bool) – whether to save mesh-related data such as element identifiers, ownership etc. This is mainly useful for debugging.

  • time (float) – the timestamp to save within the file

  • cycle (int) – the cycle (or timestep) of the data

  • <name> – writes the assigned value to the Silo file using <name> as identifier

Note:

All data objects have to be defined on the same domain but they may be defined on separate FunctionSpace s.

Others

  • pi

Packages