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)¶
- 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)¶
- 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