# esys.escript.models Package

## Classes

- `DarcyFlow`

- `FaultSystem`

- `IncompressibleIsotropicFlowCartesian`

- `LevelSet`

- `MaxIterReached`

- `Mountains`

- `PowerLaw`

- `Rheology`

- `StokesProblemCartesian`

- `SubSteppingException`

- `TemperatureCartesian`

- `Tracer`

class esys.escript.models.DarcyFlow(domain, useReduced=False, solver='POST', verbose=False, w=1.0)

  solves the problem

  u_i+k_{ij}*p_{,j} = g_i u_{i,i} = f

  where p represents the pressure and u the Darcy flux. k represents the permeability,

  Variables

    - EVAL – direct pressure gradient evaluation for flux

    - POST – global postprocessing of flux by solving the PDE K_{ij} u_j + (w * K * l u_{k,k})_{,i}= - p_{,j} + K_{ij} g_j where l is the length scale, K is the inverse of the permeability tensor, and w is a positive weighting factor.

    - SMOOTH – global smoothing by solving the PDE K_{ij} u_j= - p_{,j} + K_{ij} g_j

  __init__(domain, useReduced=False, solver='POST', verbose=False, w=1.0)

    initializes the Darcy flux problem.

    Parameters

      - domain (`Domain`) – domain of the problem

      - useReduced (`bool`) – uses reduced oreder on flux and pressure

      - solver (in [`DarcyFlow.EVAL`, `DarcyFlow.POST`, `DarcyFlow.SMOOTH` ]) – solver method

      - verbose (`bool`) – if `True` some information on the iteration progress are printed.

      - w (`float`) – weighting factor for `DarcyFlow.POST` solver

  EVAL = 'EVAL'
  POST = 'POST'
  SIMPLE = 'EVAL'
  SMOOTH = 'SMOOTH'
  getFlux(p, u0=None)

    returns the flux for a given pressure `p` where the flux is equal to `u0` on locations where `location_of_fixed_flux` is positive (see `setValue`). Notice that `g` is used, see `setValue`.

    Parameters

      - p (scalar value on the domain (e.g. `escript.Data`).) – pressure.

      - u0 (vector values on the domain (e.g. `escript.Data`) or `None`) – flux on the locations of the domain marked be `location_of_fixed_flux`.

    Returns

      flux

    Return type

      `escript.Data`

  getSolverOptionsFlux()

    Returns the solver options used to solve the flux problems :return: `SolverOptions`

  getSolverOptionsPressure()

    Returns the solver options used to solve the pressure problems :return: `SolverOptions`

  setSolverOptionsFlux(options=None)

    Sets the solver options used to solve the flux problems If `options` is not present, the options are reset to default :param options: `SolverOptions`

  setSolverOptionsPressure(options=None)

    Sets the solver options used to solve the pressure problems If `options` is not present, the options are reset to default

    Parameters

      options – `SolverOptions`

    Note

      if the adaption of subtolerance is choosen, the tolerance set by `options` will be overwritten before the solver is called.

  setValue(f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None)

    assigns values to model parameters

    Parameters

      - f (scalar value on the domain (e.g. `escript.Data`)) – volumetic sources/sinks

      - g (vector values on the domain (e.g. `escript.Data`)) – flux sources/sinks

      - location_of_fixed_pressure (scalar value on the domain (e.g. `escript.Data`)) – mask for locations where pressure is fixed

      - location_of_fixed_flux (vector values on the domain (e.g. `escript.Data`)) – mask for locations where flux is fixed.

      - permeability (scalar or symmetric tensor values on the domain (e.g. `escript.Data`)) – permeability tensor. If scalar `s` is given the tensor with `s` on the main diagonal is used.

    Note

      the values of parameters which are not set by calling `setValue` are not altered.

    Note

      at any point on the boundary of the domain the pressure (`location_of_fixed_pressure` >0) or the normal component of the flux (`location_of_fixed_flux[i]>0`) if direction of the normal is along the x_i axis.

  solve(u0, p0)

    solves the problem.

    Parameters

      - u0 (vector value on the domain (e.g. `escript.Data`).) – initial guess for the flux. At locations in the domain marked by `location_of_fixed_flux` the value of `u0` is kept unchanged.

      - p0 (scalar value on the domain (e.g. `escript.Data`).) – initial guess for the pressure. At locations in the domain marked by `location_of_fixed_pressure` the value of `p0` is kept unchanged.

    Returns

      flux and pressure

    Return type

      `tuple` of `escript.Data`.

class esys.escript.models.FaultSystem(dim=3)

  The FaultSystem class defines a system of faults in the Earth’s crust.

  A fault system is defined by set of faults index by a tag. Each fault is defined by a starting point V0 and a list of strikes `strikes` and length `l`. The strikes and the length is used to define a polyline with points `V[i]` such that

  - `V[0]=V0`

  - `V[i]=V[i]+ls[i]*array(cos(strikes[i]),sin(strikes[i]),0)`

  So `strikes` defines the angle between the direction of the fault segment and the x0 axis. ls[i]==0 is allowed.

  In case of a 3D model a fault plane is defined through a dip and depth.

  The class provides a mechanism to parametrise each fault with the domain [0,w0_max] x [0, w1_max] (to [0,w0_max] in the 2D case).

  __init__(dim=3)

    Sets up the fault system

    Parameters

      dim (`int` of value 2 or 3) – spatial dimension

  MIN_DEPTH_ANGLE = 0.1
  NOTAG = '__NOTAG__'
  addFault(strikes, ls, V0=[0.0, 0.0, 0.0], tag=None, dips=None, depths=None, w0_offsets=None, w1_max=None)

    adds a new fault to the fault system. The fault is named by `tag`.

    The fault is defined by a starting point V0 and a list of `strikes` and length `ls`. The strikes and the length is used to define a polyline with points `V[i]` such that

    - `V[0]=V0`

    - `V[i]=V[i]+ls[i]*array(cos(strikes[i]),sin(strikes[i]),0)`

    So `strikes` defines the angle between the direction of the fault segment and the x0 axis. In 3D `ls[i]` ==0 is allowed.

    In case of a 3D model a fault plane is defined through a dip `dips` and depth `depths`. From the dip and the depth the polyline `bottom` of the bottom of the fault is computed.

    Each segment in the fault is decribed by the for vertices `v0=top[i]`, `v1==top[i+1]`, `v2=bottom[i]` and `v3=bottom[i+1]` The segment is parametrized by `w0` and `w1` with `w0_offsets[i]<=w0<=w0_offsets[i+1]` and `-w1_max<=w1<=0`. Moreover

    - `(w0,w1)=(w0_offsets[i] , 0)->v0`

    - `(w0,w1)=(w0_offsets[i+1], 0)->v1`

    - `(w0,w1)=(w0_offsets[i] , -w1_max)->v2`

    - `(w0,w1)=(w0_offsets[i+1], -w1_max)->v3`

    If no `w0_offsets` is given,

    - `w0_offsets[0]=0`

    - `w0_offsets[i]=w0_offsets[i-1]+L[i]`

    where `L[i]` is the length of the segments on the top in 2D and in the middle of the segment in 3D.

    If no `w1_max` is given, the average fault depth is used.

    Parameters

      - strikes (`list` of `float`) – list of strikes. This is the angle of the fault segment direction with x0 axis. Right hand rule applies.

      - ls (`list` of `float`) – list of fault lengths. In the case of a 3D fault a segment may have length 0.

      - V0 (`list` or `numpy.array` with 2 or 3 components. `V0[2]` must be zero.) – start point of the fault

      - tag (`float` or `str`) – the tag of the fault. If fault `tag` already exists it is overwritten.

      - dips (`list` of `float`) – list of dip angles. Right hand rule around strike direction applies.

      - depths (`list` of `float`) – list of segment depth. Value mut be positive in the 3D case.

      - w0_offsets (`list` of `float` or `None`) – `w0_offsets[i]` defines the offset of the segment `i` in the fault to be used in the parametrization of the fault. If not present the cumulative length of the fault segments is used.

      - w1_max (`float`) – the maximum value used for parametrization of the fault in the depth direction. If not present the mean depth of the fault segments is used.

    Note

      In the three dimensional case the lists `dip` and `top` must have the same length.

  getBottomPolyline(tag=None)

    returns the list of the vertices defining the bottom of the fault `tag` :param tag: the tag of the fault :type tag: `float` or `str` :return: the list of vertices. In the 2D case None is returned.

  getCenterOnSurface()

    returns the center point of the fault system at the surface :rtype: `numpy.array`

  getDepthVectors(tag=None)

    returns the list of the depth vector at top vertices in fault `tag`. :param tag: the tag of the fault :type tag: `float` or `str` :return: the list of segment depths. In the 2D case None is returned.

  getDepths(tag=None)

    returns the list of the depths of the segements in fault `tag`. :param tag: the tag of the fault :type tag: `float` or `str` :return: the list of segment depths. In the 2D case None is returned.

  getDim()

    returns the spatial dimension :rtype: `int`

  getDips(tag=None)

    returns the list of the dips of the segements in fault `tag` :param tag: the tag of the fault :type tag: `float` or `str` :return: the list of segment dips. In the 2D case None is returned.

  getLengths(tag=None)
    Returns

      the lengths of segments in fault `tag`

    Return type

      `list` of `float`

  getMaxValue(f, tol=1.4901161193847656e-08)

    returns the tag of the fault of where `f` takes the maximum value and a `Locator` object which can be used to collect values from `Data` class objects at the location where the minimum is taken.

    Parameters

      - f (`escript.Data`) – a distribution of values

      - tol (`tol`) – relative tolerance used to decide if point is on fault

    Returns

      the fault tag the maximum is taken, and a `Locator` object to collect the value at location of maximum value.

  getMediumDepth(tag=None)

    returns the medium depth of fault `tag` :rtype: `float`

  getMinValue(f, tol=1.4901161193847656e-08)

    returns the tag of the fault of where `f` takes the minimum value and a `Locator` object which can be used to collect values from `Data` class objects at the location where the minimum is taken.

    Parameters

      - f (`escript.Data`) – a distribution of values

      - tol (`tol`) – relative tolerance used to decide if point is on fault

    Returns

      the fault tag the minimum is taken, and a `Locator` object to collect the value at location of minimum value.

  getOrientationOnSurface()

    returns the orientation of the fault system in RAD on the surface around the fault system center :rtype: `float`

  getParametrization(x, tag=None, tol=1.4901161193847656e-08, outsider=None)

    returns the parametrization of the fault `tag` in the fault system. In fact the values of the parametrization for at given coordinates `x` is returned. In addition to the value of the parametrization a mask is returned indicating if the given location is on the fault with given tolerance `tol`.

    Typical usage of the this method is

    dom=Domain(..) x=dom.getX() fs=FaultSystem() fs.addFault(tag=3,…) p, m=fs.getParametrization(x, outsider=0,tag=3) saveDataCSV(‘x.csv’,p=p, x=x, mask=m)

    to create a file with the coordinates of the points in `x` which are on the fault (as `mask=m`) together with their location `p` in the fault coordinate system.

    Parameters

      - x (`escript.Data` object or `numpy.array`) – location(s)

      - tag – the tag of the fault

      - tol (`float`) – relative tolerance to check if location is on fault.

      - outsider (`float`) – value used for parametrization values outside the fault. If not present an appropriate value is choosen.

    Returns

      the coordinates `x` in the coordinate system of the fault and a mask indicating coordinates in the fault by 1 (0 elsewhere)

    Return type

      `escript.Data` object or `numpy.array`

  getSegmentNormals(tag=None)

    returns the list of the normals of the segments in fault `tag` :param tag: the tag of the fault :type tag: `float` or `str` :return: the list of vectors normal to the segments. In the 2D case None is returned.

  getSideAndDistance(x, tag=None)

    returns the side and the distance at `x` from the fault `tag`.

    Parameters

      - x (`escript.Data` object or `numpy.array`) – location(s)

      - tag – the tag of the fault

    Returns

      the side of `x` (positive means to the right of the fault, negative to the left) and the distance to the fault. Note that a value zero for the side means that that the side is undefined.

  getStart(tag=None)

    returns the starting point of fault `tag` :rtype: `numpy.array`.

  getStrikeVectors(tag=None)
    Returns

      the strike vectors of fault `tag`

    Return type

      `list` of `numpy.array`.

  getStrikes(tag=None)
    Returns

      the strike of the segements in fault `tag`

    Return type

      `list` of `float`

  getTags()

    returns a list of the tags used by the fault system :rtype: `list`

  getTopPolyline(tag=None)

    returns the polyline used to describe fault tagged by `tag`

    Parameters

      tag (`float` or `str`) – the tag of the fault

    Returns

      the list of vertices defining the top of the fault. The coordinates are `numpy.array`.

  getTotalLength(tag=None)
    Returns

      the total unrolled length of fault `tag`

    Return type

      `float`

  getW0Offsets(tag=None)

    returns the offsets for the parametrization of fault `tag`.

    Returns

      the offsets in the parametrization

    Return type

      `list` of `float`

  getW0Range(tag=None)

    returns the range of the parameterization in `w0` :rtype: two `float`

  getW1Range(tag=None)

    returns the range of the parameterization in `w1` :rtype: two `float`

  transform(rot=0, shift=array([0., 0., 0.]))

    applies a shift and a consecutive rotation in the x0x1 plane.

    Parameters

      - rot (`float`) – rotation angle in RAD

      - shift (`numpy.array` of size 2 or 3) – shift vector to be applied before rotation

class esys.escript.models.IncompressibleIsotropicFlowCartesian(domain, stress=0, v=0, p=0, t=0, numMaterials=1, verbose=True)

  This class implements the rheology of an isotropic Kelvin material.

  Typical usage:

      sp = IncompressibleIsotropicFlowCartesian(domain, stress=0, v=0)
      sp.initialize(...)
      v,p = sp.solve()

  Note

    This model has been used in the self-consistent plate-mantle model proposed in Hans-Bernd Muhlhaus and Klaus Regenauer-Lieb: “Towards a self-consistent plate mantle model that includes elasticity: simple benchmarks and application to basic modes of convection”, see doi: 10.1111/j.1365-246X.2005.02742.x

  __init__(domain, stress=0, v=0, p=0, t=0, numMaterials=1, verbose=True)

    Initializes the model.

    Parameters

      - domain (`Domain`) – problem domain

      - stress (a tensor value/field of order 2) – initial (deviatoric) stress

      - v (a vector value/field) – initial velocity field

      - p (a scalar value/field) – initial pressure

      - t (`float`) – initial time

      - numMaterials (`int`) – number of materials

      - verbose (`bool`) – if `True` some information is printed.

  getCurrentEtaEff()

    returns the effective viscosity used in the last iteration step of the last time step.

  initialize(F=None, f=None, fixed_v_mask=None, v_boundary=None, restoration_factor=None)

    sets external forces and velocity constraints

    Parameters

      - F (vector value/field) – external force

      - f (vector value/field on boundary) – surface force

      - fixed_v_mask (vector value/field) – location of constraints maked by positive values

      - v_boundary (vector value/field) – value of velocity at location of constraints

      - restoration_factor (scalar values/field) – factor for normal restoration force

    Note

      Only changing parameters need to be specified.

  update(dt, iter_max=10, eta_iter_max=20, verbose=False, usePCG=True, max_correction_steps=100)

    Updates stress, velocity and pressure for time increment dt.

    Parameters

      - dt – time increment

      - iter_max – maximum number of iteration steps in the incompressible solver

      - eta_iter_max – maximum number of iteration steps in the incompressible solver

      - verbose – prints some infos in the incompressible solver

  updateStokesEquation(v, p)

    updates the underlying Stokes equation to consider dependencies from `v` and `p`

class esys.escript.models.LevelSet(phi, reinit_max=10, reinitialize_after=1, smooth=2.0, useReducedOrder=False)

  The level set method tracking an interface defined by the zero contour of the level set function phi which defines the signed distance of a point x from the interface. The contour phi(x)=0 defines the interface.

  It is assumed that phi(x)<0 defines the volume of interest, phi(x)>0 the outside world.

  __init__(phi, reinit_max=10, reinitialize_after=1, smooth=2.0, useReducedOrder=False)

    Sets up the level set method.

    Parameters

      - phi – the initial level set function

      - reinit_max – maximum number of reinitialization steps

      - reinitialize_after – `phi` is reinitialized every `reinit_after` step

      - smooth – smoothing width

  getAdvectionSolverOptions()

    Returns the solver options for the interface advective.

  getDomain()

    Returns the domain.

  getH()

    Returns the mesh size.

  getInterface(phi=None, smoothing_width=None)

    creates a characteristic function which is 1 over the over the length 2*h*smoothing_width around the interface and zero elsewhere

  getJumpingParameter(param_neg=- 1, param_pos=1, phi=None)

    Creates a function with `param_neg` where `phi<0` and `param_pos` where `phi>0` (no smoothing).

    Parameters

      - param_neg – value of parameter on the negative side (phi<0)

      - param_pos – value of parameter on the positive side (phi>0)

      - phi – level set function to be used. If not present the current level set is used.

  getLevelSetFunction()

    Returns the level set function.

  getReinitializationSolverOptions()

    Returns the options of the solver for the reinitialization

  getSmoothedJump(phi=None, smoothing_width=None)

    Creates a smooth interface from -1 to 1 over the length 2*h*smoothing_width where -1 is used where the level set is negative and 1 where the level set is 1.

  getSmoothedParameter(param_neg=- 1, param_pos=1, phi=None, smoothing_width=None)

    Creates a smoothed function with `param_neg` where `phi<0` and `param_pos` where `phi>0` which is smoothed over a length `smoothing_width` across the interface.

    Parameters

      smoothing_width – width of the smoothing zone relative to mesh size. If not present the initial value of `smooth` is used.

  getTimeStepSize(flux)

    Returns a new `dt` for a given `flux` using the Courant condition.

    Parameters

      flux – flux field

  getVolume()

    Returns the volume of the phi(x)<0 region.

  makeCharacteristicFunction(contour=0, phi=None, positiveSide=True, smoothing_width=None)

    Makes a smooth characteristic function of the region `phi(x)>contour` if `positiveSide` and `phi(x)<contour` otherwise.

    Parameters

      - phi – level set function to be used. If not present the current level set is used.

      - smoothing_width – width of the smoothing zone relative to mesh size. If not present the initial value of `smooth` is used.

  update(dt)

    Updates the level set function.

    Parameters

      dt – time step forward

class esys.escript.models.MaxIterReached

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

  __init__(*args, **kwargs)

class esys.escript.models.Mountains(domain, eps=0.01)

  The Mountains class is defined by the following equations:

  1. eps*w_i,aa+w_i,33=0 where 0<=eps<<1 and a=1,2 and w_i is the extension of the surface velocity where w_i(x_3=1)=v_i.

  2. Integration of topography PDE using Taylor-Galerkin upwinding to stabilize the advection terms H^(t+dt)=H^t+dt*w_3+w_hat*dt*[(div(w_hat*H^t)+w_3)+(dt/2)+H^t], where w_hat=w*[1,1,0], dt<0.5*d/max(w_i), d is a characteristic element size; H(x_3=1)=lambda (?)

  __init__(domain, eps=0.01)

    Sets up the level set method.

    Parameters

      - domain – the domain where the mountains is used

      - eps – the smoothing parameter for (1)

  getDomain()

    Returns the domain.

  getSafeTimeStepSize()

    Returns the time step value.

    Return type

      `float`

  getSolverOptionsForSmooting()

    returns the solver options for the smoothing/extrapolation

  getSolverOptionsForUpdate()

    returns the solver options for the topograthy update

  getTopography()

    returns the current topography. :rtype: scalar `Data`

  getVelocity()

    returns the smoothed/extrapolated velocity :rtype: vector `Data`

  setTopography(H=None)

    set the topography to H where H defines the vertical displacement. H is defined for the entire domain.

    Parameters

      H (scalar) – the topography. If None zero is used.

  setVelocity(v=None)

    set a new velocity. v is define on the entire domain but only the surface values are used.

    Parameters

      v (vector) – velocity field. If None zero is used.

  update(dt=None, allow_substeps=True)

    Sets a new W and updates the H function.

    Parameters

      dt (positve `float` which is less or equal than the safe time step size.) – time step forward. If None the save time step size is used.

class esys.escript.models.PowerLaw(numMaterials=1, verbose=False)

  this implements the power law for a composition of a set of materials where the viscosity eta of each material is given by a power law relationship of the form

  eta=eta_N*(tau/tau_t)**(1./power-1.)

  where tau is equivalent stress and eta_N, tau_t and power are given constant. Moreover an elastic component can be considered. Moreover tau meets the Drucker-Prager type yield condition

  tau <= tau_Y + friction * pressure

  where gamma_dot is the equivalent.

  __init__(numMaterials=1, verbose=False)

    initializes a power law

    Parameters

      - numMaterials (`int`) – number of materials

      - verbose (`bool`) – if `True` some information is printed.

  getElasticShearModulus()

    returns the elastic shear modulus.

    Returns

      elastic shear modulus

  getEtaEff(gamma_dot, eta0=None, pressure=None, dt=None, iter_max=30)

    returns the effective viscosity eta_eff such that

    tau=eta_eff * gamma_dot

    by solving a non-linear problem for tau.

    Parameters

      - gamma_dot – equivalent strain gamma_dot

      - eta0 – initial guess for the effective viscosity (e.g from a previous time step). If not present, an initial guess is calculated.

      - pressure – pressure used to calculate yield condition

      - dt (positive `float` if present) – time step size. only needed if elastic component is considered.

      - iter_max (`int`) – maximum number of iteration steps.

    Returns

      effective viscosity.

  getEtaN(id=None)

    returns the viscosity

    Parameters

      id (`int`) – if present, the viscosity for material `id` is returned.

    Returns

      the list of the viscosities for all matrials is returned. If `id` is present only the viscosity for material `id` is returned.

  getEtaTolerance()

    returns the relative tolerance for the effectice viscosity.

    Returns

      relative tolerance

    Return type

      positive `float`

  getFriction()

    returns the friction coefficient

    Returns

      friction coefficient

  getNumMaterials()

    returns the numebr of materials

    Returns

      number of materials

    Return type

      `int`

  getPower(id=None)

    returns the power in the power law

    Parameters

      id (`int`) – if present, the power for material `id` is returned.

    Returns

      the list of the powers for all matrials is returned. If `id` is present only the power for material `id` is returned.

  getTauT(id=None)

    returns the transition stress

    Parameters

      id (`int`) – if present, the transition stress for material `id` is returned.

    Returns

      the list of the transition stresses for all matrials is returned. If `id` is present only the transition stress for material `id` is returned.

  getTauY()

    returns the yield stress

    Returns

      the yield stress

  setDruckerPragerLaw(tau_Y=None, friction=None)

    Sets the parameters for the Drucker-Prager model.

    Parameters

      - tau_Y – yield stress

      - friction – friction coefficient

  setElasticShearModulus(mu=None)

    Sets the elastic shear modulus.

    Parameters

      mu – elastic shear modulus

  setEtaTolerance(rtol=0.0001)

    sets the relative tolerance for the effectice viscosity.

    Parameters

      rtol (positive `float`) – relative tolerance

  setPowerLaw(eta_N, id=0, tau_t=1, power=1)

    Sets the power-law parameters for material id

    Parameters

      - id (`int`) – material id

      - eta_N – viscosity for tau=tau_t

      - tau_t – transition stress

      - power – power law coefficient

  setPowerLaws(eta_N, tau_t, power)

    Sets the parameters of the power-law for all materials.

    Parameters

      - eta_N – list of viscosities for tau=tau_t

      - tau_t – list of transition stresses

      - power – list of power law coefficient

  validMaterialId(id=0)

    checks if a given material id is valid

    Parameters

      id (`int`) – a material id

    Returns

      `True` is the id is valid

    Return type

      `bool`

class esys.escript.models.Rheology(domain, stress=None, v=None, p=None, t=0, verbose=True)

  General framework to implement a rheology

  __init__(domain, stress=None, v=None, p=None, t=0, verbose=True)

    Initializes the rheology

    Parameters

      - domain (`Domain`) – problem domain

      - stress (a tensor value/field of order 2) – initial (deviatoric) stress

      - v (a vector value/field) – initial velocity field

      - p (a scalar value/field) – initial pressure

      - t (`float`) – initial time

  checkVerbose()

    Returns True if verbose is switched on

    Returns

      value of verbosity flag

    Return type

      `bool`

  getDeviatoricStrain(v=None)

    Returns deviatoric strain of current velocity or if `v` is present the deviatoric strain of velocity `v`:

    Parameters

      v (`Data` of rank 1) – a velocity field

    Returns

      deviatoric strain of the current velocity field or if `v` is present the deviatoric strain of velocity `v`

    Return type

      `Data` of rank 2

  getDeviatoricStress()

    Returns current deviatoric stress.

    Returns

      current deviatoric stress

    Return type

      `Data` of rank 2

  getDomain()

    returns the domain.

    Returns

      the domain

    Return type

      `Domain`

  getForce()

    Returns the external force

    Returns

      external force

    Return type

      `Data`

  getGammaDot(D=None)

    Returns current second invariant of deviatoric strain rate or if `D` is present the second invariant of `D`.

    Parameters

      D (`Data` of rank 0) – deviatoric strain rate tensor

    Returns

      second invariant of deviatoric strain

    Return type

      scalar `Data`

  getPressure()

    Returns current pressure.

    Returns

      current stress

    Return type

      scalar `Data`

  getRestorationFactor()

    Returns the restoring force factor

    Returns

      restoring force factor

    Return type

      `float` or `Data`

  getStress()

    Returns current stress.

    Returns

      current stress

    Return type

      `Data` of rank 2

  getSurfaceForce()

    Returns the surface force

    Returns

      surface force

    Return type

      `Data`

  getTau()

    Returns current second invariant of deviatoric stress

    Returns

      second invariant of deviatoric stress

    Return type

      scalar `Data`

  getTime()

    Returns current time.

    Returns

      current time

    Return type

      `float`

  getVelocity()

    Returns current velocity.

    Returns

      current velocity

    Return type

      vector `Data`

  getVelocityConstraint()

    Returns the constraint for the velocity as a pair of the mask of the location of the constraint and the values.

    Returns

      the locations of fixed velocity and value of velocities at these locations

    Return type

      `tuple` of `Data` s

  setDeviatoricStrain(D=None)

    set deviatoric strain

    Parameters

      D (`Data` of rank 2) – new deviatoric strain. If `D` is not present the current velocity is used.

  setDeviatoricStress(stress)

    Sets the current deviatoric stress

    Parameters

      stress (`Data` of rank 2) – new deviatoric stress

  setExternals(F=None, f=None, fixed_v_mask=None, v_boundary=None, restoration_factor=None)

    sets external forces and velocity constraints

    Parameters

      - F (vector value/field) – external force

      - f (vector value/field on boundary) – surface force

      - fixed_v_mask (vector value/field) – location of constraints maked by positive values

      - v_boundary (vector value/field) – value of velocity at location of constraints

      - restoration_factor (scalar values/field) – factor for normal restoration force

    Note

      Only changing parameters need to be specified.

  setGammaDot(gammadot=None)

    set the second invariant of deviatoric strain rate. If `gammadot` is not present zero is used.

    Parameters

      gammadot (`Data` of rank 1) – second invariant of deviatoric strain rate.

  setPressure(p)

    Sets current pressure. :param p: new deviatoric stress :type p: scalar `Data`

  setStatus(t, v, p, stress)

    Resets the current status given by pressure p and velocity v.

    Parameters

      - t (`float`) – new time mark

      - v (vector `Data`) – new current velocity

      - p (scalar `Data`) – new deviatoric stress

      - stress (`Data` of rank 2) – new deviatoric stress

  setTime(t=0.0)

    Updates current time.

    Parameters

      t (`float`) – new time mark

  setVelocity(v)

    Sets current velocity.

    Parameters

      v (vector `Data`) – new current velocity

class esys.escript.models.StokesProblemCartesian(domain, **kwargs)

  solves

  > -(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j}

  >   u_{i,i}=0

  > u=0 where fixed_u_mask>0 eta*(u_{i,j}+u_{j,i})*n_j-p*n_i=surface_stress +stress_{ij}n_j

  if surface_stress is not given 0 is assumed.

  typical usage:

  > sp=StokesProblemCartesian(domain) sp.setTolerance() sp.initialize(…) v,p=sp.solve(v0,p0) sp.setStokesEquation(…) # new values for some parameters v1,p1=sp.solve(v,p)

  __init__(domain, **kwargs)

    initialize the Stokes Problem

    The approximation spaces used for velocity (=Solution(domain)) and pressure (=ReducedSolution(domain)) must be LBB complient, for instance using quadratic and linear approximation on the same element or using linear approximation with macro elements for the pressure.

    Parameters

      domain (`Domain`) – domain of the problem.

  Bv(v, tol)

    returns inner product of element p and div(v)

    Parameters

      v – a residual

    Returns

      inner product of element p and div(v)

    Return type

      `float`

  getDV(p, v, tol)

    return the value for v for a given p

    Parameters

      - p – a pressure

      - v – a initial guess for the value v to return.

    Returns

      dv given as Adv=(f-Av-B^*p)

  getSolverOptionsDiv()

    returns the solver options for solving the equation to project the divergence of the velocity onto the function space of presure.

    Return type

      `SolverOptions`

  getSolverOptionsPressure()

    returns the solver options used solve the equation for pressure. :rtype: `SolverOptions`

  getSolverOptionsVelocity()

    returns the solver options used solve the equation for velocity.

    Return type

      `SolverOptions`

  initialize(f=<esys.escriptcore.escriptcpp.Data object>, fixed_u_mask=<esys.escriptcore.escriptcpp.Data object>, eta=1, surface_stress=<esys.escriptcore.escriptcpp.Data object>, stress=<esys.escriptcore.escriptcpp.Data object>, restoration_factor=0)

    assigns values to the model parameters

    Parameters

      - f (`Vector` object in `FunctionSpace` `Function` or similar) – external force

      - fixed_u_mask (`Vector` object on `FunctionSpace` `Solution` or similar) – mask of locations with fixed velocity.

      - eta (`Scalar` object on `FunctionSpace` `Function` or similar) – viscosity

      - surface_stress (`Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar) – normal surface stress

      - stress (`Tensor` object on `FunctionSpace` `Function` or similar) – initial stress

  inner_p(p0, p1)

    Returns inner product of p0 and p1

    Parameters

      - p0 – a pressure

      - p1 – a pressure

    Returns

      inner product of p0 and p1

    Return type

      `float`

  inner_pBv(p, Bv)

    returns inner product of element p and Bv=-div(v)

    Parameters

      - p – a pressure increment

      - Bv – a residual

    Returns

      inner product of element p and Bv=-div(v)

    Return type

      `float`

  norm_Bv(Bv)

    Returns Bv (overwrite).

    Return type

      equal to the type of p

    Note

      boundary conditions on p should be zero!

  norm_v(v)

    returns the norm of v

    Parameters

      v – a velovity

    Returns

      norm of v

    Return type

      non-negative `float`

  setSolverOptionsDiv(options=None)

    set the solver options for solving the equation to project the divergence of the velocity onto the function space of presure.

    Parameters

      options (`SolverOptions`) – new solver options

  setSolverOptionsPressure(options=None)

    set the solver options for solving the equation for pressure. :param options: new solver options :type options: `SolverOptions`

  setSolverOptionsVelocity(options=None)

    set the solver options for solving the equation for velocity.

    Parameters

      options (`SolverOptions`) – new solver options

  setStokesEquation(f=None, fixed_u_mask=None, eta=None, surface_stress=None, stress=None, restoration_factor=None)

    assigns new values to the model parameters.

    Parameters

      - f (`Vector` object in `FunctionSpace` `Function` or similar) – external force

      - fixed_u_mask (`Vector` object on `FunctionSpace` `Solution` or similar) – mask of locations with fixed velocity.

      - eta (`Scalar` object on `FunctionSpace` `Function` or similar) – viscosity

      - surface_stress (`Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar) – normal surface stress

      - stress (`Tensor` object on `FunctionSpace` `Function` or similar) – initial stress

  solve_AinvBt(p, tol)

    Solves Av=B^*p with accuracy `tol`

    Parameters

      p – a pressure increment

    Returns

      the solution of Av=B^*p

    Note

      boundary conditions on v should be zero!

  solve_prec(Bv, tol)

    applies preconditioner for for BA^{-1}B^* to Bv with accuracy `self.getSubProblemTolerance()`

    Parameters

      Bv – velocity increment

    Returns

      p=P(Bv) where P^{-1} is an approximation of BA^{-1}B^ * )

    Note

      boundary conditions on p are zero.

  updateStokesEquation(v, p)

    updates the Stokes equation to consider dependencies from `v` and `p` :note: This method can be overwritten by a subclass. Use `setStokesEquation` to set new values to model parameters.

class esys.escript.models.SubSteppingException

  Thrown if the L{Mountains} class uses substepping.

  __init__(*args, **kwargs)

class esys.escript.models.TemperatureCartesian(domain, **kwargs)

  Represents and solves the temperature advection-diffusion problem

  rhocp(T_{,t} + v_i T_{,i} - ( k T_{,i})_i = Q

  k T_{,i}*n_i=surface_flux and T_{,t} = 0 where `given_T_mask`>0.

  If surface_flux is not given 0 is assumed.

  Typical usage:

      sp = TemperatureCartesian(domain)
      sp.setTolerance(1.e-4)
      t = 0
      T = ...
      sp.setValues(rhocp=...,  v=..., k=..., given_T_mask=...)
      sp.setInitialTemperature(T)
      while t < t_end:
          sp.setValue(Q=...)
          T = sp.getTemperature(dt)
          t += dt

  __init__(domain, **kwargs)

    Initializes the temperature advection-diffusion problem.

    Parameters

      domain – domain of the problem

    Note

      the approximation order is switched to reduced if the approximation order is nnot linear (equal to 1).

  getTemperature(dt, **kwargs)

    Same as `getSolution`.

  setInitialTemperature(T)

    Same as `setInitialSolution`.

  setValue(rhocp=None, v=None, k=None, Q=None, surface_flux=None, given_T_mask=None)

    Sets new values to coefficients.

    Parameters

      - coefficients – new values assigned to coefficients

      - M (any type that can be cast to a `Data` object on `Function`) – value for coefficient `M`

      - M_reduced (any type that can be cast to a `Data` object on `Function`) – value for coefficient `M_reduced`

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

      - m (any type that can be cast to a `Data` object on `FunctionOnBoundary`) – value for coefficient `m`

      - m_reduced (any type that can be cast to a `Data` object on `FunctionOnBoundary`) – value for coefficient `m_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 the location of constraints

    Raises

      IllegalCoefficient – if an unknown coefficient keyword is used

class esys.escript.models.Tracer(domain, useBackwardEuler=False, **kwargs)

  Represents and solves the tracer problem

  C_{,t} + v_i C_{,i} - ( k T_{,i})_i) = 0

  C_{,t} = 0 where `given_C_mask`>0. C_{,i}*n_i=0

  Typical usage:

      sp = Tracer(domain)
      sp.setTolerance(1.e-4)
      t = 0
      T = ...
      sp.setValues(given_C_mask=...)
      sp.setInitialTracer(C)
      while t < t_end:
          sp.setValue(v=...)
          dt.getSaveTimeStepSize()
          C = sp.getTracer(dt)
          t += dt

  __init__(domain, useBackwardEuler=False, **kwargs)

    Initializes the Tracer advection problem

    Parameters

      - domain – domain of the problem

      - useBackwardEuler (`bool`) – if set the backward Euler scheme is used. Otherwise the Crank-Nicholson scheme is applied. Not that backward Euler scheme will return a safe time step size which is practically infinity as the scheme is unconditional unstable. So other measures need to be applied to control the time step size. The Crank-Nicholson scheme provides a higher accuracy but requires to limit the time step size to be stable.

    Note

      the approximation order is switched to reduced if the approximation order is nnot linear (equal to 1).

  getTracer(dt, **kwargs)

    Same as `getSolution`.

  setInitialTracer(C)

    Same as `setInitialSolution`.

  setValue(v=None, given_C_mask=None, k=None)

    Sets new values to coefficients.

    Parameters

      - coefficients – new values assigned to coefficients

      - M (any type that can be cast to a `Data` object on `Function`) – value for coefficient `M`

      - M_reduced (any type that can be cast to a `Data` object on `Function`) – value for coefficient `M_reduced`

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

      - m (any type that can be cast to a `Data` object on `FunctionOnBoundary`) – value for coefficient `m`

      - m_reduced (any type that can be cast to a `Data` object on `FunctionOnBoundary`) – value for coefficient `m_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 the location of constraints

    Raises

      IllegalCoefficient – if an unknown coefficient keyword is used

## Functions

## Others

## Packages