esys.escript.models Package

Classes

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

optionsSolverOptions

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