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 problemuseReduced (
bool
) – uses reduced oreder on flux and pressuresolver (in [
DarcyFlow.EVAL
,DarcyFlow.POST
,DarcyFlow.SMOOTH
]) – solver methodverbose (
bool
) – ifTrue
some information on the iteration progress are printed.w (
float
) – weighting factor forDarcyFlow.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 tou0
on locations wherelocation_of_fixed_flux
is positive (seesetValue
). Notice thatg
is used, seesetValue
.- Parameters
p (scalar value on the domain (e.g.
escript.Data
).) – pressure.u0 (vector values on the domain (e.g.
escript.Data
) orNone
) – flux on the locations of the domain marked belocation_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/sinksg (vector values on the domain (e.g.
escript.Data
)) – flux sources/sinkslocation_of_fixed_pressure (scalar value on the domain (e.g.
escript.Data
)) – mask for locations where pressure is fixedlocation_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 scalars
is given the tensor withs
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 bylocation_of_fixed_flux
the value ofu0
is kept unchanged.p0 (scalar value on the domain (e.g.
escript.Data
).) – initial guess for the pressure. At locations in the domain marked bylocation_of_fixed_pressure
the value ofp0
is kept unchanged.
- Returns
flux and pressure
- Return type
tuple
ofescript.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 lengthl
. The strikes and the length is used to define a polyline with pointsV[i]
such thatV[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 lengthls
. The strikes and the length is used to define a polyline with pointsV[i]
such thatV[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 3Dls[i]
==0 is allowed.In case of a 3D model a fault plane is defined through a dip
dips
and depthdepths
. From the dip and the depth the polylinebottom
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]
andv3=bottom[i+1]
The segment is parametrized byw0
andw1
withw0_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
offloat
) – list of strikes. This is the angle of the fault segment direction with x0 axis. Right hand rule applies.ls (
list
offloat
) – list of fault lengths. In the case of a 3D fault a segment may have length 0.V0 (
list
ornumpy.array
with 2 or 3 components.V0[2]
must be zero.) – start point of the faulttag (
float
orstr
) – the tag of the fault. If faulttag
already exists it is overwritten.dips (
list
offloat
) – list of dip angles. Right hand rule around strike direction applies.depths (
list
offloat
) – list of segment depth. Value mut be positive in the 3D case.w0_offsets (
list
offloat
orNone
) –w0_offsets[i]
defines the offset of the segmenti
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
andtop
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
orstr
: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
orstr
: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
orstr
: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
orstr
: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
offloat
-
getMaxValue
(f, tol=1.4901161193847656e-08)¶ returns the tag of the fault of where
f
takes the maximum value and aLocator
object which can be used to collect values fromData
class objects at the location where the minimum is taken.- Parameters
f (
escript.Data
) – a distribution of valuestol (
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 aLocator
object which can be used to collect values fromData
class objects at the location where the minimum is taken.- Parameters
f (
escript.Data
) – a distribution of valuestol (
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 coordinatesx
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 tolerancetol
.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 (asmask=m
) together with their locationp
in the fault coordinate system.- Parameters
x (
escript.Data
object ornumpy.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 ornumpy.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
orstr
: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 faulttag
.- Parameters
x (
escript.Data
object ornumpy.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
ofnumpy.array
.
-
getStrikes
(tag=None)¶ - Returns
the strike of the segements in fault
tag
- Return type
list
offloat
-
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
orstr
) – 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
offloat
-
getW0Range
(tag=None)¶ returns the range of the parameterization in
w0
:rtype: twofloat
-
getW1Range
(tag=None)¶ returns the range of the parameterization in
w1
:rtype: twofloat
-
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 RADshift (
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 domainstress (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 timenumMaterials (
int
) – number of materialsverbose (
bool
) – ifTrue
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
andp
-
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 everyreinit_after
stepsmooth – 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
wherephi<0
andparam_pos
wherephi>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
wherephi<0
andparam_pos
wherephi>0
which is smoothed over a lengthsmoothing_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 givenflux
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
ifpositiveSide
andphi(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)¶ Initialize self. See help(type(self)) for accurate signature.
-
-
class
esys.escript.models.
Mountains
(domain, eps=0.01)¶ The Mountains class is defined by the following equations:
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.
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 materialsverbose (
bool
) – ifTrue
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 materialid
is returned.- Returns
the list of the viscosities for all matrials is returned. If
id
is present only the viscosity for materialid
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 materialid
is returned.- Returns
the list of the powers for all matrials is returned. If
id
is present only the power for materialid
is returned.
-
getTauT
(id=None)¶ returns the transition stress
- Parameters
id (
int
) – if present, the transition stress for materialid
is returned.- Returns
the list of the transition stresses for all matrials is returned. If
id
is present only the transition stress for materialid
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 ideta_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 domainstress (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 velocityv
:- 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 velocityv
- 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 ofD
.- 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
orData
-
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
ofData
s
-
setDeviatoricStrain
(D=None)¶ set deviatoric strain
- Parameters
D (
Data
of rank 2) – new deviatoric strain. IfD
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 markv (vector
Data
) – new current velocityp (scalar
Data
) – new deviatoric stressstress (
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 inFunctionSpace
Function
or similar) – external forcefixed_u_mask (
Vector
object onFunctionSpace
Solution
or similar) – mask of locations with fixed velocity.eta (
Scalar
object onFunctionSpace
Function
or similar) – viscositysurface_stress (
Vector
object onFunctionSpace
FunctionOnBoundary
or similar) – normal surface stressstress (
Tensor
object onFunctionSpace
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 inFunctionSpace
Function
or similar) – external forcefixed_u_mask (
Vector
object onFunctionSpace
Solution
or similar) – mask of locations with fixed velocity.eta (
Scalar
object onFunctionSpace
Function
or similar) – viscositysurface_stress (
Vector
object onFunctionSpace
FunctionOnBoundary
or similar) – normal surface stressstress (
Tensor
object onFunctionSpace
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
andp
:note: This method can be overwritten by a subclass. UsesetStokesEquation
to set new values to model parameters.
-
class
esys.escript.models.
SubSteppingException
¶ Thrown if the L{Mountains} class uses substepping.
-
__init__
(*args, **kwargs)¶ Initialize self. See help(type(self)) for accurate signature.
-
-
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 onFunction
) – value for coefficientM
M_reduced (any type that can be cast to a
Data
object onFunction
) – value for coefficientM_reduced
A (any type that can be cast to a
Data
object onFunction
) – value for coefficientA
A_reduced (any type that can be cast to a
Data
object onReducedFunction
) – value for coefficientA_reduced
B (any type that can be cast to a
Data
object onFunction
) – value for coefficientB
B_reduced (any type that can be cast to a
Data
object onReducedFunction
) – value for coefficientB_reduced
C (any type that can be cast to a
Data
object onFunction
) – value for coefficientC
C_reduced (any type that can be cast to a
Data
object onReducedFunction
) – value for coefficientC_reduced
D (any type that can be cast to a
Data
object onFunction
) – value for coefficientD
D_reduced (any type that can be cast to a
Data
object onReducedFunction
) – value for coefficientD_reduced
X (any type that can be cast to a
Data
object onFunction
) – value for coefficientX
X_reduced (any type that can be cast to a
Data
object onReducedFunction
) – value for coefficientX_reduced
Y (any type that can be cast to a
Data
object onFunction
) – value for coefficientY
Y_reduced (any type that can be cast to a
Data
object onReducedFunction
) – value for coefficientY_reduced
m (any type that can be cast to a
Data
object onFunctionOnBoundary
) – value for coefficientm
m_reduced (any type that can be cast to a
Data
object onFunctionOnBoundary
) – value for coefficientm_reduced
d (any type that can be cast to a
Data
object onFunctionOnBoundary
) – value for coefficientd
d_reduced (any type that can be cast to a
Data
object onReducedFunctionOnBoundary
) – value for coefficientd_reduced
y (any type that can be cast to a
Data
object onFunctionOnBoundary
) – value for coefficienty
d_contact (any type that can be cast to a
Data
object onFunctionOnContactOne
orFunctionOnContactZero
) – value for coefficientd_contact
d_contact_reduced (any type that can be cast to a
Data
object onReducedFunctionOnContactOne
orReducedFunctionOnContactZero
) – value for coefficientd_contact_reduced
y_contact (any type that can be cast to a
Data
object onFunctionOnContactOne
orFunctionOnContactZero
) – value for coefficienty_contact
y_contact_reduced (any type that can be cast to a
Data
object onReducedFunctionOnContactOne
orReducedFunctionOnContactZero
) – value for coefficienty_contact_reduced
d_dirac (any type that can be cast to a
Data
object onDiracDeltaFunctions
) – value for coefficientd_dirac
y_dirac (any type that can be cast to a
Data
object onDiracDeltaFunctions
) – value for coefficienty_dirac
r (any type that can be cast to a
Data
object onSolution
orReducedSolution
depending on whether reduced order is used for the solution) – values prescribed to the solution at the locations of constraintsq (any type that can be cast to a
Data
object onSolution
orReducedSolution
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=0Typical 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 onFunction
) – value for coefficientM
M_reduced (any type that can be cast to a
Data
object onFunction
) – value for coefficientM_reduced
A (any type that can be cast to a
Data
object onFunction
) – value for coefficientA
A_reduced (any type that can be cast to a
Data
object onReducedFunction
) – value for coefficientA_reduced
B (any type that can be cast to a
Data
object onFunction
) – value for coefficientB
B_reduced (any type that can be cast to a
Data
object onReducedFunction
) – value for coefficientB_reduced
C (any type that can be cast to a
Data
object onFunction
) – value for coefficientC
C_reduced (any type that can be cast to a
Data
object onReducedFunction
) – value for coefficientC_reduced
D (any type that can be cast to a
Data
object onFunction
) – value for coefficientD
D_reduced (any type that can be cast to a
Data
object onReducedFunction
) – value for coefficientD_reduced
X (any type that can be cast to a
Data
object onFunction
) – value for coefficientX
X_reduced (any type that can be cast to a
Data
object onReducedFunction
) – value for coefficientX_reduced
Y (any type that can be cast to a
Data
object onFunction
) – value for coefficientY
Y_reduced (any type that can be cast to a
Data
object onReducedFunction
) – value for coefficientY_reduced
m (any type that can be cast to a
Data
object onFunctionOnBoundary
) – value for coefficientm
m_reduced (any type that can be cast to a
Data
object onFunctionOnBoundary
) – value for coefficientm_reduced
d (any type that can be cast to a
Data
object onFunctionOnBoundary
) – value for coefficientd
d_reduced (any type that can be cast to a
Data
object onReducedFunctionOnBoundary
) – value for coefficientd_reduced
y (any type that can be cast to a
Data
object onFunctionOnBoundary
) – value for coefficienty
d_contact (any type that can be cast to a
Data
object onFunctionOnContactOne
orFunctionOnContactZero
) – value for coefficientd_contact
d_contact_reduced (any type that can be cast to a
Data
object onReducedFunctionOnContactOne
orReducedFunctionOnContactZero
) – value for coefficientd_contact_reduced
y_contact (any type that can be cast to a
Data
object onFunctionOnContactOne
orFunctionOnContactZero
) – value for coefficienty_contact
y_contact_reduced (any type that can be cast to a
Data
object onReducedFunctionOnContactOne
orReducedFunctionOnContactZero
) – value for coefficienty_contact_reduced
d_dirac (any type that can be cast to a
Data
object onDiracDeltaFunctions
) – value for coefficientd_dirac
y_dirac (any type that can be cast to a
Data
object onDiracDeltaFunctions
) – value for coefficienty_dirac
r (any type that can be cast to a
Data
object onSolution
orReducedSolution
depending on whether reduced order is used for the solution) – values prescribed to the solution at the locations of constraintsq (any type that can be cast to a
Data
object onSolution
orReducedSolution
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
-