##~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~##
## ##
## This file forms part of the Underworld geophysics modelling application. ##
## ##
## For full license and copyright information, please refer to the LICENSE.md file ##
## located at the project root, or contact the authors. ##
## ##
##~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~##
import underworld as uw
import underworld._stgermain as _stgermain
from . import sle
import underworld.libUnderworld as libUnderworld
import numpy
[docs]class SteadyStateDarcyFlow(_stgermain.StgCompoundComponent):
"""
This class provides functionality for a discrete representation
of the steady state darcy flow equation.
The class uses a standard Galerkin finite element method to
construct a system of linear equations which may then be solved
using an object of the underworld.system.Solver class.
The underlying element types are determined by the supporting
mesh used for the 'pressureField'.
The strong form of the given boundary value problem, for :math:`f`,
:math:`q` and :math:`h` given, is
.. math::
\\begin{align}
q_i =& \\kappa \\, ( -u_{,i} + S_i ) & \\\\
q_{i,i} =& \\: f & \\text{ in } \\Omega \\\\
u =& \\: q & \\text{ on } \\Gamma_q \\\\
-q_i n_i =& \\: h & \\text{ on } \\Gamma_h \\\\
\\end{align}
where,
* :math:`\\kappa` is the diffusivity, :math:`u` is the pressure,
* :math:`S` is a flow body-source, for example due to gravity,
* :math:`f` is a source term, :math:`q` is the Dirichlet condition, and
* :math:`h` is a Neumann condition.
The problem boundary, :math:`\\Gamma`, admits the decomposition
:math:`\\Gamma=\\Gamma_q\\cup\\Gamma_h` where
:math:`\\emptyset=\\Gamma_q\\cap\\Gamma_h`. The equivalent weak form is:
.. math::
-\\int_{\\Omega} w_{,i} \\, q_i \\, d \\Omega = \\int_{\\Omega} w \\, f \\, d\\Omega + \\int_{\\Gamma_h} w \\, h \\, d \\Gamma
where we must find :math:`u` which satisfies the above for all :math:`w`
in some variational space.
Parameters
----------
pressureField : underworld.mesh.MeshVariable
The solution field for pressure.
fn_diffusivity : underworld.function.Function
The function that defines the diffusivity across the domain.
fn_bodyforce : underworld.function.Function
A function that defines the flow body-force across the domain,
for example gravity. Must be a vector. Optional.
voronoi_swarm : underworld.swarm.Swarm
A swarm with just one particle within each cell should be provided.
This avoids the evaluation of the velocity on nodes and inaccuracies
arising from diffusivity changes within cells.
If a swarm is provided, voronoi type numerical integration is
utilised. The provided swarm is used as the basis for the voronoi
integration. If no voronoi_swarm is provided, Gauss integration
is used.
conditions : underworld.conditions.SystemCondition
Numerical conditions to impose on the system. This should be supplied as
the condition itself, or a list object containing the conditions.
gauss_swarm : underworld.swarm.GaussIntegrationSwarm
If provided this gauss_swarm will be used for (gaussian) numerical
integration rather than a default gauss integration swarm that is
automatically generated and dependent on the element order of the mesh.
NB: if a voronoi_swarm is defined it OVERRIDES this gauss_swarm as the
preferred integration swarm (quadrature method).
velocityField : underworld.mesh.MeshVariable
Solution field for darcy flow velocity. Optional.
swarmVarVelocity : undeworld.swarm.SwarmVariable
If a swarm variable is provided, the velocity calculated on the swarm will be stored.
This is the most representative velocity data object, as the velocity calculation occurs
on the swarm, away from mesh nodes. Optional.
Notes
-----
Constructor must be called collectively by all processes.
"""
_objectsDict = { "_system" : "SystemLinearEquations" }
_selfObjectName = "_system"
def __init__(self, pressureField, fn_diffusivity, fn_bodyforce=None, voronoi_swarm=None, conditions=[], velocityField=None, swarmVarVelocity=None, _removeBCs=True, gauss_swarm=None, **kwargs):
if not isinstance( pressureField, uw.mesh.MeshVariable):
raise TypeError( "Provided 'pressureField' must be of 'MeshVariable' class." )
self._pressureField = pressureField
try:
_fn_diffusivity = uw.function.Function.convert(fn_diffusivity)
except Exception as e:
raise uw._prepend_message_to_exception(e, "Exception encountered. Note that provided 'fn_diffusivity' must be of or convertible to 'Function' class.\nEncountered exception message:\n")
if not fn_bodyforce:
if pressureField.mesh.dim == 2:
fn_bodyforce = (0.,0.)
else:
fn_bodyforce = (0.,0.,0.)
try:
_fn_bodyforce = uw.function.Function.convert(fn_bodyforce)
except Exception as e:
raise uw._prepend_message_to_exception(e, "Exception encountered. Note that provided 'fn_bodyforce' must be of or convertible to 'Function' class.\nEncountered exception message:\n")
if voronoi_swarm and not isinstance(voronoi_swarm, uw.swarm.Swarm):
raise TypeError( "Provided 'swarm' must be of 'Swarm' class." )
self._swarm = voronoi_swarm
if len(numpy.unique(voronoi_swarm.owningCell.data[:])) != len(voronoi_swarm.owningCell.data[:]):
import warnings
warnings.warn("It is not advised to fill any cell with more than one particle, as the Q1 shape function cannot capture material interfaces. Use at your own risk.")
if velocityField and not isinstance( velocityField, uw.mesh.MeshVariable):
raise TypeError( "Provided 'velocityField' must be of 'MeshVariable' class." )
self._velocityField = velocityField
if swarmVarVelocity and not isinstance(swarmVarVelocity, uw.swarm.SwarmVariable):
raise TypeError("Provided 'swarmVarVelocity' must be of 'SwarmVariable' class.")
self._swarmVarVelocity = swarmVarVelocity
if voronoi_swarm and pressureField.mesh.elementType=='Q2':
import warnings
warnings.warn("Voronoi integration may yield unsatisfactory results for Q2 element types. Q2 element types can also result in spurious velocity calculations.")
if not isinstance( _removeBCs, bool):
raise TypeError( "Provided '_removeBCs' must be of type bool." )
self._removeBCs = _removeBCs
# error check dcs 'dirichlet conditions' and ncs 'neumann cond.' FeMesh_IndexSets
nbc = None
mesh = pressureField.mesh
if not isinstance(conditions,(list,tuple)):
conditionslist = []
conditionslist.append(conditions)
conditions = conditionslist
for cond in conditions:
if not isinstance( cond, uw.conditions.SystemCondition ):
raise TypeError( "Provided 'conditions' must be a list 'SystemCondition' objects." )
# set the bcs on here
if type( cond ) == uw.conditions.DirichletCondition:
if cond.variable == pressureField:
libUnderworld.StgFEM.FeVariable_SetBC( pressureField._cself, cond._cself )
elif type( cond ) == uw.conditions.NeumannCondition:
nbc=cond
else:
raise RuntimeError("Can't decide on input condition")
# setup the gauss integration swarm
if gauss_swarm != None:
if type(gauss_swarm) != uw.swarm.GaussIntegrationSwarm:
raise RuntimeError( "Provided 'gauss_swarm' must be a GaussIntegrationSwarm object" )
intswarm = gauss_swarm
else:
intswarm = uw.swarm.GaussIntegrationSwarm(mesh)
# we will use voronoi if that has been requested by the user, else use
# gauss integration.
if self._swarm:
intswarm = self._swarm._voronoi_swarm
# need to ensure voronoi is populated now, as assembly terms will call
# initial test functions which may require a valid voronoi swarm
self._swarm._voronoi_swarm.repopulate()
# build the equation numbering for the temperature field discretisation
tEqNums = self._tEqNums = sle.EqNumber( pressureField, removeBCs=self._removeBCs )
# create solutions vector
self._solutionVector = sle.SolutionVector(pressureField, tEqNums )
libUnderworld.StgFEM.SolutionVector_LoadCurrentFeVariableValuesOntoVector( self._solutionVector._cself )
# create force vectors
self._fvector = sle.AssembledVector(pressureField, tEqNums )
# and matrices
self._kmatrix = sle.AssembledMatrix( self._solutionVector, self._solutionVector, rhs=self._fvector )
self._kMatTerm = sle.MatrixAssemblyTerm_NA_i__NB_i__Fn( integrationSwarm=intswarm,
assembledObject=self._kmatrix,
fn=_fn_diffusivity)
self._forceVecTerm = sle.VectorAssemblyTerm_NA_i__Fn_i( integrationSwarm=intswarm,
assembledObject=self._fvector,
fn=fn_bodyforce*fn_diffusivity)
# search for neumann conditions
for cond in conditions:
if isinstance( cond, uw.conditions.NeumannCondition ):
#NOTE many NeumannConditions can be used but the _sufaceFluxTerm only records the last
### -VE flux because of Energy_SLE_Solver ###
negativeCond = uw.conditions.NeumannCondition( fn_flux=-1.0*cond.fn_flux,
variable=cond.variable,
indexSetsPerDof=cond.indexSet )
self._surfaceFluxTerm = sle.VectorSurfaceAssemblyTerm_NA__Fn__ni(
assembledObject = self._fvector,
surfaceGaussPoints = 2,
nbc = negativeCond )
super(SteadyStateDarcyFlow, self).__init__(**kwargs)
def _setup(self):
uw.libUnderworld.StGermain.Stg_ObjectList_Append( self._cself.stiffnessMatrices, self._kmatrix._cself )
uw.libUnderworld.StGermain.Stg_ObjectList_Append( self._cself.forceVectors, self._fvector._cself )
uw.libUnderworld.StGermain.Stg_ObjectList_Append( self._cself.solutionVectors, self._solutionVector._cself )
def _add_to_stg_dict(self,componentDictionary):
# call parents method
super(SteadyStateDarcyFlow,self)._add_to_stg_dict(componentDictionary)
def solve_velocityField(self):
fnVel = (-self._pressureField.fn_gradient + self.fn_bodyforce) * self._kMatTerm.fn
if self._swarmVarVelocity:
self._swarmVarVelocity.data[:] = fnVel.evaluate(self._swarm)
if self._velocityField:
velproj = uw.utils.MeshVariable_Projection(self._velocityField,fnVel,self._swarm)
velproj.solve()
@property
def fn_diffusivity(self):
"""
The diffusivity function. You may change this function directly via this
property.
"""
return self._kMatTerm.fn
@fn_diffusivity.setter
def fn_diffusivity(self, value):
self._kMatTerm.fn = value
@property
def fn_bodyforce(self):
"""
The heating function. You may change this function directly via this
property.
"""
return self._forceVecTerm.fn / self._kMatTerm.fn
@fn_bodyforce.setter
def fn_bodyforce(self, value):
self._forceVecTerm.fn = value * self._kMatTerm.fn