Underworld¶
Underworld is a Python friendly version of the Underworld code which provides a programmable and flexible front end to all the functionality of the code running in a parallel HPC environment. This gives signficant advantages to the user, with access to the power of Python libraries for setup of complex problems, analysis at runtime, problem steering, and multi physics coupling. While Underworld2 embraces Jupyter Notebooks as the preferred modelling environment, only standard Python is required.
The Underworld2 development team is based in Melbourne, Australia at the University of Melbourne and at Monash University led by Louis Moresi. We would like to acknowledge AuScope Simulation, Analysis and Modelling for providing long term funding which has made the project possible. Additional funding for specific improvements and additional functionality has come from the Australian Research Council. The Python toolkit was funded by the NeCTAR eresearch_tools program. Underworld was originally developed in collaboration with the Victorian Partnership for Advanced Computing.
Privacy¶
Note that basic usage metrics are dispatched when you use Underworld. We do this to help assess the usage of our code which is important in justifying our funding. To opt out, set the UW_NO_USAGE_METRICS environment variable. See PRIVACY.md for full details.
Bedtime reading¶
These papers explain the theory and implementation for Underworld. The code itself can also be cited via the zenodo DOI. There is a master DOI for all releases and releases after 2.6.0 are automatically given a DOI under the master. If you are using a development branch and wish to obtain a DOI for your specific version we ask that you contact us to make an interim release under the master DOI.
Moresi, L., Dufour, F., and Muhlhaus, H.B., 2002, Mantle convection modeling with viscoelastic/brittle lithosphere: Numerical methodology and plate tectonic modeling: Pure And Applied Geophysics, v. 159, no. 10, p. 2335–2356, doi: 10.1007/s00024-002-8738-3.
Moresi, L., Dufour, F., and Muhlhaus, H.B., 2003, A Lagrangian integration point finite element method for large deformation modeling of viscoelastic geomaterials: Journal of Computational Physics, v. 184, no. 2, p. 476–497.
Moresi, L., Quenette, S., Lemiale, V., Mériaux, C., Appelbe, W., Mühlhaus, 2007, Computational approaches to studying non-linear dynamics of the crust and mantle: Phys. Earth Planet. Inter, v. 163, p. 69–82, doi: 10.1016/j.pepi.2007.06.009.
Geodynamics - mathematical background¶
The simplest template set of equations for solid-Earth dynamics cover mass, momentum and heat conservation in a highly viscous fluid allowing for additional effects due to elasticity and plasticity. The Stokes momentum equation neglects inertia but includes an additional term on the right hand side that represents stress history associated with an explicit treatment of viscoelasticity.[1] [2] [3]
\(\tau\) is the deviatoric stress, \(p\) represents the pressure, \(\rho\) is density, \(T\) is the temperature, \(C\) is a concentration intended to represent changes in composition.
At pressures in planetary interiors, silicate minerals are weakly compressible and this is generally considered as a perturbation to an incompressible flow ignoring bulk viscosity and only considering the long-term elastic resistance to volume change. For the purposes of explaining the formulation, the incompressible constraint equation on the velocity, \(u\) is sufficient.
The thermal evolution of the system expresses the balance between heat transport by fluid motion, thermal diffusion and internal heat generation. Additional terms can be included to account for heating due to viscous dissipation, for example, but do not change the overall character of the conservation equation.
The most significant feature of this system is the spontaneous appearance of boundary layers where horizontal advection and vertical diffusion are approximately balanced. By contrast, compositional variations are characterised by a much smaller, usually negligible, rate of diffusion:
The thermal and compositional variations couple to the momentum equation through their effect on density. The Boussinesq approximation [6], accounts for the buoyancy forces while neglecting the associated volume change allowing us to assume incompressibility. If the non-diffusive, compositional variation represents a smoothly varying concentration, then the density can be written as [4]
In the case where C represents a state with discrete steps (e.g.a phase change or immiscible fluids), it is common to let \(\rho_0\) take discrete values and assume \(\alpha_C=0\).
The final requirement is a constitutive relationship for the momentum equation that links the stress to the velocity unknown. Rheology is one of the defining aspects of the dynamics of the mantle, particularly in the cooler parts of the upper boundary layer where elasticity, non-linearity, and brittle behaviour plays a significant role. A general constitutive law can be expressed as:
where \(\mu\) is the elastic shear modulus and \(\eta\) is the shear viscosity (both of which may vary with temperature and composition). \(\Lambda\) is a structural tensor that represents the orientation of the plastic deformation relative to the applied stress and \(\lambda\) is a scalar multiplier that is computed to satisfy the stress conditions at yield [5]. Typically, \(\eta\) varies by several tens of orders of magnitude over the typical temperature ranges expected between the Earth’s surface and interior. The orientation tensor and the yield stress are usually modelled to include a simple damage evolution that relates to the work expended in deforming the material at yield.
- [1]:
- Moresi, L. N., F. Dufour, and H. B. Muhlhaus (2002), Mantle convection modeling with viscoelastic/brittle lithosphere: Numerical methodology and plate tectonic modeling, Pure And Applied Geophysics, 159(10), 2335–2356, doi:10.1007/s00024-002-8738-3.
- [2]:
- Moresi, L. N., F. Dufour, and H. B. Muhlhaus (2003), A Lagrangian integration point finite element method for large deformation modeling of viscoelastic geomaterials, Journal of Computational Physics, 184(2), 476–497, doi:10.1016/S0021-9991(02)00031-1.
- [3]:
- Farrington, R. J., L. N. Moresi, and F. A. Capitanio (2014), The role of viscoelasticity in subducting plates, Geochemistry, Geophysics, Geosystems, 15(11), 4291–4304, doi:10.1002/2014GC005507.
- [4]:
- van Keken, P. E., S. D. King, H. Schmeling, U. R. Christensen, D. Neumeister, and M. P. Doin (1997), A comparison of methods for the modeling of thermochemical convection, J. Geophys. Res., 102(B10), 22,477–22,495.
- [5]:
- Moresi, H. B. Muhlhaus, V. Lemiale, and D. A. May (2007), Incompressible viscous formulations for deformation and yielding of the lithosphere, Geological Society London Special Publications, 282(1), 457–472, doi:10.1144/SP282.19.
- [6]:
- Boussinesq, J. (1903), The ́orie analytique de la chaleur, Vol. 2, Gauthier-Villars, Paris (Reproduction Bibliothe ́que Nationale de France, 1995).
Numerical methods - background¶
Description¶
Underworld is a Lagrangian integration point finite element code. This is a modernization of the original particle-in-cell concept from the 1960s in which a structured mesh and an unstructured particle swarm co-exist. The mesh is used to solve diffusion-dominated parts of the problem and the particle swarm is used to track advected quantities. In the finite element context, the mapping from mesh to particles is through the usual basis functions of the elements and the mapping from particles to mesh is through the integration scheme used to build up the stiffness matrices etc.
The applications of the method are mainly in modelling of complex fluids where very large strains occur but the material also has a memory of the entire strain / strain-rate history. In geosciences this occurs due to the visco-elasticity of rocks at lithospheric temperature and their tendency to develop fabric (lattice preferred orientation and stress/strain-dependent grain size). Problems with material interfaces which undergo severe distortion during the deformation are also naturally handled by this method provided there is no slip on the interface.
Background¶
The method has been published in detail in Moresi et al (2002, 2003)[1]. These papers dealt exclusively with 2D applications but in recent years, we have introduced a number of improvements in the method to enable us to scale the problem to 3D. For example we developed a fast discrete Voronoi method to compute the integration weights of the particle-to-mesh mapping efficiently [2]. We have also concentrated on extremely robust solvers / preconditioners which are necessary because the material variations and geometrical complexity are both large and unpredictable at the start of the simulation.
The benefit of this approach is associated with the separation of the computational mesh from the swarm of points which track the history. This allows us to retain a much more structured computational mesh than the deformation / material history would otherwise allow. We can take full advantage of the most efficient geometrical multigrid solvers and there is no need to preserve structure during any remeshing operations we undertake (for example if we do need to track a free surface or an internal interface). Although there are several complexities introduced by enforcing this separation, we find that the benefits, for our particular class of problems, are significant.
Implementation and parallelism¶
Underworld is implemented using the StGermain framework . This provides the essential infrastructure to manage i/o, meshes, particle swarms, finite element operations, in a parallel (domain decomposition, message passing) environment. The numerical solvers are based around the PETSc software suite which focuses on delivering good parallel scalability. Good scalability results have been achieved for over 10000 core simulations.
- [1]:
- Moresi, L. N., F. Dufour, and H. B. Muhlhaus (2003), A Lagrangian integration point finite element method for large deformation modeling of viscoelastic geomaterials, Journal of Computational Physics, 184(2), 476–497, doi:10.1016/S0021-9991(02)00031-1.
- [2]:
- Velić, M., D. A. May, and L. N. Moresi (2009), A fast robust algorithm for computing discrete voronoi diagrams, Journal of Mathematical Modelling and …, doi:10.1007/s10852-008-9097-6.
- [3]
- Moresi, F. Dufour, and H. B. Muhlhaus. A Lagrangian integration point finite element method for large deformation modeling of viscoelastic geomaterials. Journal Of Computational Physics, 184:476–497, 2003.
Underworld Installation¶
Docker¶
Docker is a type of lightweight virtualisation, and is the preferred method for Underworld usage on personal computers. You will first need to install Docker on your system, and then you may install Underworld via Docker. Docker may be driven from the command line, but new users may wish to use the Docker Kitematic GUI instead for ease. Simply search for ‘underworldcode/underworld2’ within Kitematic, and then click ‘CREATE’ to launch a container. You will eventually wish to modify your container settings (again through Kitematic) to enable local folder volume mapping, which will allow you to access your local drives within your container.
For Linux users, and those who prefer the command line, the following minimal command should be sufficient to access the Underworld2 Jupyter Notebook examples:
docker run -p 8888:8888 underworldcode/underworld2
Navigate to localhost:8888 to see the notebooks. A number of useful commands are provided within the Underworld cheat-sheet.
Windows users should note that for Windows 10 Home you should install Docker Toolbox, while for Windows 10 Professional you should install Docker for Windows. The Docker Toolbox edition utilised VirtualBox for virtualisation, and therefore to access any running Jupyter servers you must browse to the virtual machine address (instead of localhost). To find the VM address, you will generally execute
docker-machine ip default
but note that this will only work correctly from the Docker Quickstart Terminal.
Native¶
For installation on HPC facilities and if you would like a native local build, you will need to download, install and compile the Underworld code and relevant dependencies. A native build can be difficult in non-standard environments and we are currently investigating HPC deployment of Docker containers.
For information on how to build, please refer to the top level COMPILE.md file. Instructions for HPC builds may be found at docs/install_guides. You may also find useful build related information at the Underworld blog pages.
User Guide¶
Getting started¶
Welcome to using Underworld!¶
Underworld is a Python library for the development of long time scale earth process models. The Underworld Python interface is designed to facilitate interactive and intuitive model development. To this end, we embrace Jupyter Notebooks as the preferred development environment, although the standard Python interpreter is the only requirement. Underworld utilises MPI parallelisation to allow large simulation across HPC facilities. The Python interface was partly funded by the NeCTAR eResearch_tools program.
Underworld includes the glucifer package for interactive
visualisation within the Jupter Notebook environment, and for seamless
visualisation capabilities across parallel simulations. glucifer
leverages the LavaVu engine for
rendering capabilities.
Resources¶
There are numerous resources that you might find useful for learning and using Underworld:
- underworld2.readthedocs.io: The page where most of our primary documentation is published, including this user guide, and the API documentation. It also links out to most other resources.
- github.com/underworldcode/underworld2:
Our code repository. You can keep track of changes to the codes base
here, and you can also browse most of our documention directly within
the Github interface. Don’t forget to
StarandWatchour project if you find it useful! - github.com/underworldcode/underworld2/issues: Please use our issue tracker to report any bugs or difficulties you encounter. This includes any general usage questions as well as technial issues. If you wish to submit a suspected bug, please include a minimal example which allows us to reproduce the issue. It is also the place to post any feature requests you may have!
- The User Guide: A more focused look at the various aspects of
Underworld modelling. Note that each section is a self container
Jupyter Notebook document, available in the Underworld repository
(
docs/user_guide). It is also published at the Underworld ReadTheDocs page. - Examples: These notebooks go through the entire Underworld
workflow for constructing and solving geophysics models. These models
demonstrate Underworld current best usage practises, and are
guaranteed to operate correctly for each Underworld release. They are
available in the repository at
docs/examples. - API Documentation: The Underworld API documentation is published
at our ReadTheDocs page, along with the glucifer API documentation.
Note that the API documentation is embedded in the Python
implementation as docstrings, and is therefore also available
directy via the Python
help()built-in. More usefully, this information is accessible within Jupyter Notebooks via the tooltips shortcutShift-Tab(when in edit mode). - Underworld Models Library: A repository of Underworld models. The library includes models which reproduce publication results, tutorials, examples and usage tidbits. Note that these models are not explicitly maintained, and so may not operate against the latest version of Underworld. If you scroll to the end of each model, it should state which version of Underworld the repository model was successfully ran against.
Required Python Skills¶
To use Underworld successfully, you will need to have an understanding of the following Python constructs:
- Basic Python, such as importing modules, and syntax structure (indenting, functions, etc).
- Containers such as dictionaries, lists and tuples.
- Flow control (loops, if-else conditionals, etc).
- Python objects, object methods, object attributes, object lifecycles.
Most beginner (or intermediate) Python tutorials should cover these concepts. Also useful, though not strictly necessary, is some familiarity with the following:
- Exception handling (for dealing with errors that might occur).
- Context managers (for mesh and swarm deformations).
- Operator overloading.
Note that Underworld heavily leverages the numpy Python numerical
library for all heavy data access and manipulation. All Underworld
objects that record heavy data will expose their data via the data
attribute, which is actually a handle to a numpy array. As such,
familiary with numpy usage paradigms is a must, and more advanced usage
patterns (array slicing, advanced indexing, etc) will become important
as your models increase in complexity.
Similarly, Underworld uses h5py for all heavy data disk IO. H5py is
a Python library which provides a Python interface to writing/reading
HDF5 format files. While not strictly required, more advanced users
will certainly find having some familiarity with the h5py libary
useful, both for directly querying files Underworld has generated, and
also for writing their own files (in preference to CSV for example).
Jupyter Notebooks¶
Jupyter Notebooks is the recommended environment for most model
development. In Underworld we utilise notebooks to provide inline
visualisation of your model configurations, allowing you to quickly see
your results, modify as required, and then regenerate and repeat.
Equally important is the tooltips and autocomplete functionality
provided within the notebooks. To access tooltips, use shift-tab
while you have the text editing cursor located within an Underworld
object (or you can write a question mark after the object, and execute
the cell). For autocomplete, after you type a few letters you can press
tab to be provided with all possible completion options. Using these
tools is essential to rapid and frustration free model development,
especially for new users.
If you are new to Jupyter Notebooks, you should familiarise yourself with the notebook environment first. Also, remember the Help menu bar option provides useful references and keyboard shortcuts.
How to get help¶
If you encounter issues or suspect a bug, please create a ticket using the issue tracker on github.
A quick demo¶
Let’s do a quick run through of setting up some basic Underworld objects.
# Import underworld
import underworld as uw
# Create a mesh:
mesh = uw.mesh.FeMesh_Cartesian( elementRes = (8, 8),
minCoord = (0., 0.),
maxCoord = (2., 1.))
# Next we create a mesh variable:
temperatureField = mesh.add_variable( nodeDofCount=1 )
# Let's initialise the variable with some data:
for index, coord in enumerate(mesh.data):
temperatureField.data[index] = coord[1] # set the temperature to be the vertical (y) coordinate
# Finally we will plot the temperature field using ``gLucifer`` after importing the gLucifer module.
import glucifer
fig = glucifer.Figure(figsize=(800,400))
fig.append( glucifer.objects.Surface(mesh, temperatureField, colours="blue white red") )
fig.append( glucifer.objects.Mesh(mesh) )
fig.show()
Typically we might then setup boundary conditions, particle swarms, rheology and systems to be solved.
All of these topics are discussed in the following sections of the user guide.
The Mesh¶
The finite element mesh is a fundamental construct for Underworld
modelling. It will generally determine your domain geometry, and the
resolution of the finite element system. For parallel simulations, the
mesh topology will also determine the domain decomposition for the
problem. Currently underworld only provides curvilinear mesh
capabilities.
Overview:¶
- Creating mesh objects.
- Element types.
- Deforming the mesh.
- Loading and saving the mesh.
- Special sets.
- Mesh variables
- Setting values on a mesh variables.
- Gradients of mesh variable fields.
- Loading and saving mesh variable data.
Keywords: mesh variables, finite elements, load, save, initial conditions
Creating the mesh¶
First create an 2x2 element mesh. By default the mesh will be of
rectangular geometry, with domain extents specified via the minCoord
and maxCoord constructor parameters.
# This example creates a mesh, and visualises it.
import underworld as uw
import glucifer
mesh = uw.mesh.FeMesh_Cartesian( elementRes = (2, 2),
minCoord = (0.0, 0.0),
maxCoord = (2.0, 1.0) )
# visualising the result
figMesh = glucifer.Figure(figsize=(800,400))
figMesh.append( glucifer.objects.Mesh(mesh, nodeNumbers=True, pointsize=10) )
figMesh.show()
As is usually the pattern for data carrying objects in Underworld,
underlying data is available as numpy arrays via the data mesh
instances attribute. In the case of mesh objects, this data corresponds
to mesh node coordinate information.
# This example simply prints mesh node information.
import underworld as uw
mesh = uw.mesh.FeMesh_Cartesian(elementRes=(2, 2))
for item in mesh.data:
print(item)
[ 0. 0.]
[ 0.5 0. ]
[ 1. 0.]
[ 0. 0.5]
[ 0.5 0.5]
[ 1. 0.5]
[ 0. 1.]
[ 0.5 1. ]
[ 1. 1.]
Note that data in underworld is always published as arrays of shape
(num_items, count), where num_items is the number of items in
the dataset (in this case, 3x3 nodes, so 9), and count is the
size of each data item (so 1 for scalar data, 2 or 3 for
vectors in 2- and 3-dimensions respectively, etc).
In parallel, each process will be assigned responsibility for a
subsection of the domain, and as such the length (num_items) of the
arrays will generally be different across processes. As such, while it
might be tempting to use numpy.reshape to reshape your arrays to
reflect the cartesian mesh topology, this approach will not work in
parallel and is therefore not usually recommended. You should instead
traverse data arrays as 1-d lists (each item of size count).
Following this pattern will usually result in models which can safely be
run in serial or parallel.
Element Types¶
Underworld supports two primary element types: a linear element type
(Q1) and a quadratic element type (Q2). Interpolation using
these element types will provide continuous results across element
boundaries, though note that interpolant derivatives will be
discontinous across boundaries. Underworld also supports mixed element
types, where a secondary element type may also be specified. This
formulation is critical for stable solutions to Stokes type problems,
with velocity fields being represented by the primary element type, and
pressure fields by the secondary type. Secondary element types are not
continuous across element boundaries. Supported secondary types are
constant (dQ0), linear (dPc1) and bilinear (dQ1). Refer to
http://femtable.org for further details on element types.
Mesh element types are specified as a textual argument to the
elementType constructor parameter. A primary type must always be
provided (Q1 or Q2), and where a secondary type is required, the
primary/secondary types are provided as a pair (Q1/dQ0, Q2/dPc1
or Q2/dQ1). For mixed type elements, two mesh objects are generated
by the constructor, one for the primary and one for the secondary
elements. The primary mesh will be the object returned directly by the
constructor, with the secondary mesh available via the primary mesh’s
subMesh attribute.
Note that where constructor parameters are omitted, default values are
used. Defaults values are baked into the underworld API and may be
determined in Jupyter by pressing shift-tab with the cursor
inside the constructor parenthesis, or by using the Python help()
function (with the class/function of interest as the argument). Default
values are useful for quickly mocking up models, but may change without
notice so should not be relied upon.
Element boundaries (and therefore domain geometry and resolution) is entirely determined by the primary mesh, with the location of the secondary mesh nodes slaved to the primary mesh nodes. Secondary mesh nodes are visualised below.
# Create a mesh with submesh import underworld as uw import glucifer mesh = uw.mesh.FeMesh_Cartesian( elementType="Q2/dpc1", maxCoord=(2.0, 1.0) ) submesh = mesh.subMesh figMesh = glucifer.Figure(figsize=(800,400)) figMesh.append( glucifer.objects.Mesh(mesh) ) figMesh.append( glucifer.objects.Mesh(submesh, nodeNumbers=True, pointsize=10) ) figMesh.show()
Deforming the mesh¶
It is often desirable to deform the mesh to either conform the domain to some geometry, or to implement mesh refinement by bunching nodes in certain domain regions. The user is free to modify the mesh as desired, though care must be taken to ensure elements are not improperly transformed (Jacobians should not become singular) and to avoid tangling (face normals should be consistent).
To deform the mesh, the user will directly modify the node coordinates
via the usual data array. Note however that modifications are only
possible from within the deform_mesh() context manager, as otherwise
the array is read-only.
# A trivial deformation
import underworld as uw
import glucifer
mesh = uw.mesh.FeMesh_Cartesian(elementRes=(2, 2), maxCoord=(2.0, 1.0) )
with mesh.deform_mesh():
mesh.data[4][0] += 0.05
figMesh = glucifer.Figure(figsize=(800,400))
figMesh.append( glucifer.objects.Mesh(mesh, nodeNumbers=True, pointsize=10) )
figMesh.show()
Note that in the above example, for simplicity, we’ve simply moved an arbitrary mesh node. Usually you will not want to construct your logic in ways which rely on the node numbers themselves, as this will usually not be consistent across parallel simulations, and/or as you vary model resolution.
We now look at a more realistic example where we deform the mesh to increase the resolution at the top of the domain. There are countless ways this may be achieved, but here we use the simple transformation \(z \rightarrow \sqrt{z}\):
# A more useful deformation
import underworld as uw
import glucifer
mesh = uw.mesh.FeMesh_Cartesian(elementRes=(32,32), maxCoord=(2.0, 1.0))
with mesh.deform_mesh():
for index, coord in enumerate(mesh.data):
mesh.data[index][1] = mesh.data[index][1]**0.5
figMesh = glucifer.Figure(figsize=(800,400))
figMesh.append(glucifer.objects.Mesh(mesh))
figMesh.show()
Saving and loading the mesh¶
Your mesh data may be loaded or saved using the mesh instance load()
and save() method respectively. All data in Underworld is written as
HDF5 format binary files. There are numerous tools available for
interrogating HDF5 files, but in Underworld we primarily leverage the
h5py Python package. We recommend that users also familiarise
themselves with this package to allow querying of Underworld generated
data, and also as an excellent option for directly writing their own
simulation data.
You will wish to save data for analysis purposes, but also for
simulation restart. Note that Underworld does not provide any explicit
restart functionality, and you will instead need to use object
load() methods to return each object to their previous states. In
Underworld, the pattern is always to first recreate vanilla equivalent
objects, and then reload the data. So in the case of the mesh above,
you must first recreate a standard mesh object, and it must be of
identical resolution and element type. Finally you will reload the data.
# In this example we'll walk through some of the mesh loading mechanisms
import underworld as uw
mesh = uw.mesh.FeMesh_Cartesian()
# trivially deform
with mesh.deform_mesh():
mesh.data[3]+=0.1
# save
mesh.save('deformedMesh.h5')
# create second mesh
mesh2 = uw.mesh.FeMesh_Cartesian()
# confirm that it is different to first
import numpy as np
result = np.allclose(mesh.data,mesh2.data)
if result==True: # raise error if result not as expected
raise RuntimeError("Mesh objects should not match. Something wrong has happenend")
# now load data onto second mesh, and re-test
mesh2.load('deformedMesh.h5')
result = np.allclose(mesh.data,mesh2.data)
if result==False: # raise error if result not as expected
raise RuntimeError("Mesh objects should match. Something wrong has happenend")
# let's remove the saved file as it is no longer required.
if uw.rank()==0:
import os;
os.remove('deformedMesh.h5')
Special sets¶
Special sets are sets of nodes which are in some way special to a given
mesh object. For instance, for the cartesian mesh you will often wish to
tie special values to the walls to form your problem boundary
conditions. The list of special sets is provided via the specialSets
list object. The sets are named with respect to the (I,J,K)
cartesian indexing of the mesh, with the I index usually used for
the horizontal (left/right), J for the vertical, and K for
horizontal (front/back). Note that this is just a matter of convention,
and you are free to chose how the cartesian coordinates map within your
model. So, for example, MinI_VertexSet specifies the set of nodes
belonings to what we would normally consider to be the left wall, and
similarly MaxI_VertexSet the right wall. We also provide the alias
(Left_VertexSet, Right_VertexSet, etc) for the most common use
case.
import underworld as uw
mesh = uw.mesh.FeMesh_Cartesian(elementRes=(2,2))
print("Available special sets for this mesh are:")
print(mesh.specialSets.keys())
leftset = mesh.specialSets['MinI_VertexSet']
print("\nIndices in the left set wall set:")
print(leftset)
Available special sets for this mesh are:
['MaxI_VertexSet', 'Top_VertexSet', 'Left_VertexSet', 'MinI_VertexSet', 'AllWalls_VertexSet', 'Bottom_VertexSet', 'Right_VertexSet', 'MinJ_VertexSet', 'MaxJ_VertexSet', 'Empty']
Indices in the left set wall set:
FeMesh_IndexSet([0, 3, 6])
You can confirm that these are indeed the left nodes with reference to
the mesh visualisation at the top of the page. Note that the special set
indices always specify parallel local indices. In almost all
instances, you will interact with only local data and identifiers. This
pattern is critical to successful parallel operation, and is usually a
natural way of constructing your models. You should also note that the
above visualisation shows the global node numbering, though for these
serial demonstrations the values are identical. Your checkpointed mesh
(and mesh variable) data will however be ordered according to its global
numbering. If for some reason you need to determine the global
identifiers, they are stored in the data_elgId array.
Mesh Variables¶
Mesh variables are used to encode spatial data across the mesh. Underworld mesh variables assign data for each node of the mesh, so for example you might assign a temperature value at each mesh node. Mesh variables also leverage the mesh element shape functions to form interpolations within elements, allowing our discrete datasets to have continuous analogues formed.
Mesh variables may be used for any number of purposes when constructing models, but most importantly they form the unknowns for the finite element numerical systems you will be solving.
We create a mesh variable by using the mesh’s add_variable() method,
and specifying the nodeDofCount (node degree of freedom count). So
for example, setting nodeDofCount=1 yields a scalar mesh variable,
having a single data value recorded at each node of the mesh, as you may
require for perhaps a temperature field.
Setting values on the MeshVariable¶
As you will find with most Underworld data carrying objects, you will
access the data directly via the data attribute. Our data is stored
there as a numpy array, and you may therefore accessed or modified it
using standard Numpy operations.
Let’s initialise a variable with a function based on its spatial coordinates
This variable may be used to represent a temperature field. We will walk
over the mesh vertex data to access coordinate information. As the
temperature field array is (by design) a 1-1 map to the vertex array, we
know the index of the mesh vertex will be the index of the associated
temperature datum. Note that the Python built-in enumerate() acts
to return both the index of a piece of data in a list (or array),
and the data itself.
# create a mesh & meshvariable, and then initialise with data.
import underworld as uw
import glucifer
import math
mesh = uw.mesh.FeMesh_Cartesian( maxCoord=(2., 1.) )
# now create our mesh variable
temperatureField = mesh.add_variable( nodeDofCount=1 )
pi = math.pi
for index, coord in enumerate(mesh.data):
temperatureField.data[index] = 100.*math.exp(1.-coord[1])
# vis results
fig = glucifer.Figure(figsize=(800,400))
fig.append( glucifer.objects.Surface(mesh, temperatureField) )
fig.append( glucifer.objects.Mesh(mesh) )
fig.show()
To describe for example a velocity, a vector mesh variable is required.
We set nodeDofCount=2 accordingly, and initialise to:
# vector fields
import underworld as uw
import glucifer
import math
mesh = uw.mesh.FeMesh_Cartesian( maxCoord=(2., 1.) )
velocityField = mesh.add_variable( nodeDofCount=2 )
coordmid = (1., 0.5)
for index, coord in enumerate(mesh.data):
vx = coordmid[1] - coord[1]
vz = coord[0] - coordmid[0]
velocityField.data[index] = (vx, vz)
# visualise both the velocity and temperature
fig = glucifer.Figure(figsize=(800,400))
fig.append( glucifer.objects.VectorArrows(mesh, velocityField, scaling=0.2, arrowHead=0.2) )
fig.show()
Mesh variables gradients¶
Mesh variable gradients are calculated according to:
Note that gradients derived in this way are not generally continuous across element boundaries. Note that we access the different gradient components via the square bracket operator, with ordering:
or for a vector field:
The gradient of the field is accessible via the fn_gradient
attribute on the mesh variable.
Loading and saving variables¶
The procedure for loading and saving variables is very similar to that
used for the mesh, and uses the usual save() and load() methods.
For model restarts, the normal pattern applies, with the user required
to generate unmodified objects first, and then reloading the required
data.
# mesh variable loading & saving
import underworld as uw
import numpy as np
mesh = uw.mesh.FeMesh_Cartesian()
meshvariable = mesh.add_variable(1)
# this might be a pattern you use in your models.
# set `restart` to required value, but note that
# you *must* first run with `True` to ensure that
# a file has been created!
restart=False
restartFile = "meshvariable.h5"
if not restart:
# not restarting, so init as required.
meshvariable.data[:] = 1.234
else:
meshvariable.load(restartFile)
# confirm that data is correct at this point!
result = np.allclose(meshvariable.data,1.234)
if result==False: # raise error if result not as expected
raise RuntimeError("Mesh variable does not contain expected values.")
# now save. generally this will occur after
# you've done something interesting, but we
# don't have time for that in this example.
ignore = meshvariable.save(restartFile)
# let's remove the saved file as it is no longer required.
if uw.rank()==0:
import os;
os.remove(restartFile)
Swarms¶
In Underworld, the swarm is an object that defines a collection of particles. Swarms are usually used to track the Lagrangian quantities required for a given model, such as a material type identifier, or the plastic strain of an advecting parcel of fluid. Swarm may be used passively to record information as they advect, or actively where the values recorded on particles actually feed into rheologies or forces.
Swarms of particles may:
- Advect through the mesh according to a user specified velocity.
- Store arbitrary data on a per-particle basis.
- Freely cross process boundaries in parallel simulations.
The user is free to create as many swarms as required and each swarm may
contain an arbitrary number of particles. The data layout for any
given swarm is identical across all its particles, though different
swarms may have different layouts. For example, SwarmA may contain 100
particles, and each particle may encode an int and a float,
while SwarmB may contain 15 particles, each particle encoding three
float values. Particles may also be added to a swarm at any stage
either directly or through population control mechanisms. It is also
possible to delete particles, though currently this is only possible
through indirect means.
Overview¶
- Creating a swarm object and adding particles.
- Moving particles.
- Swarm variables.
- Shapes with particle swarms.
- Saving and loading swarms.
Keywords: swarms, particles, shapes.
Creating a swarm object and adding particles¶
Creating a swarm is very simple, but you will first require a mesh. The mesh will inform the swarm of your model geometry. In particular, your model domain is generally determined by the mesh, and particles may be deleted if they leave the mesh. For parallel simulations, the mesh also determines the problem partitioning across multiple processes, and therefore the local process domain.
Note that your newly created swarm will be empty! It is just a container for particles, and you will need to explicitly add particles. To add particles, you will either use a layout object, or generate a numpy array with the required particle coordinates.
import underworld as uw import glucifer import numpy as np mesh = uw.mesh.FeMesh_Cartesian(maxCoord=(2.,1.)) # create a mesh! swarm = uw.swarm.Swarm( mesh=mesh ) # now let's add particles. first create a layout object. swarmLayout = uw.swarm.layouts.PerCellSpaceFillerLayout( swarm=swarm, particlesPerCell=20 ) # now populate using layout. swarm.populate_using_layout( layout=swarmLayout ) # vis fig = glucifer.Figure(figsize=(800,400)) fig.append( glucifer.objects.Points(swarm=swarm, pointSize=5, colourBar=False) ) fig.append( glucifer.objects.Mesh(mesh) ) fig.show()
Layouts are usually the best option where you want to populate the
entire domain with particles. This is the usual use case where you
intend on using particles dynamically. Note also that layouts may only
be used to populate swarm objects which are currently empty. There are
numerous layouts available in the underworld.swarm.layouts
submodule.
Swarms can also be populated using numpy arrays that specify particle coordinates. This is how you will add particles when you wish to explicitly specify their coordinates and the total number of particles. The usual use case for this is adding passive swarms for analytic purposes, such as tracking interfaces, or extracting data at fixed coordinates throughout your simulations.
Let’s create a new swarm and add particles explicitly.
import underworld as uw
import glucifer
import numpy as np
mesh = uw.mesh.FeMesh_Cartesian(maxCoord=(2.,1.))
swarmCustom = uw.swarm.Swarm( mesh=mesh, particleEscape=True )
# create the array
swarmCoords = np.array([ [0.4,0.2], [0.8,0.4],[1.2,0.6],[1.6,0.8],[3.6,1.8]])
# use the array to add particles at the specified coordinates.
result = swarmCustom.add_particles_with_coordinates(swarmCoords)
print("Returned array: {}.".format(result))
fig = glucifer.Figure(figsize=(800,400))
fig.append( glucifer.objects.Points(swarm=swarmCustom, pointSize=10, colourBar=False) )
fig.append( glucifer.objects.Mesh(mesh))
fig.show()
Returned array: [ 0 1 2 3 -1].
Note the array returned above by the add_particles_with_coordinates
method. It specifies the local identifiers for the added particles.
Coordinates which are not within the (local) domain will be ignored and
are signified with a -1. For parallel simulation, the domain is
partitioned across all processes, and therefore coordinates which are
ignored on one process (ie, no particle created) may be consumed on
another process (ie, particle created).
Moving particles¶
To move particles, you have two options. You can explicitly move
particles using the data numpy array on the swarm (or equivalently
via particleCoordinates.data), or you can use the swarm advector
object (underworld.Systems.SwarmAdvector). The latter is a natural
choice when advecting a swarm with a velocity field that you have either
defined, or is the result of a Stokes solve. It is the usual way in
which particles are advected.
Let us start with a very basic constant velocity field with which to advect particles
# note that this cell requires the previous cell to be executed first. vel = mesh.add_variable(2) vel.data[:] = (1.,0.) advector = uw.systems.SwarmAdvector(vel,swarmCustom) dt = advector.get_max_dt() advector.integrate(dt) fig.show()
Note that the swarm advector requires a MeshVariable class velocity
field object, and does not yet support the more general Function
class.
To modify particle coordinates directly, we must use the
deform_swarm context manager, much as the deform_mesh manager is
required to modify the mesh. We will first move a single particle:
# note that this cell requires the previous cell to be executed first.
with swarmCustom.deform_swarm():
swarmCustom.data[1] = (0.4,0.5)
fig.show()
Underworld also supports particle deletion, though currently this is
only possible by moving particles outside the domain. Note that the
swarm object must also be created with particleEscape=True.
# note that this cell requires the previous cell to be executed first.
with swarmCustom.deform_swarm():
swarmCustom.data[1] = (9999,9999.)
fig.show()
Swarm variables¶
You may add data storing variables to the swarm via the
add_variable() swarm method. Each particle will record a value of
the given type and count:
import underworld as uw
import glucifer
import numpy as np
mesh = uw.mesh.FeMesh_Cartesian(maxCoord=(2.,1.))
swarmCustom = uw.swarm.Swarm( mesh=mesh, particleEscape=True )
swarmCoords = np.array([ [0.4,0.2], [0.8,0.4],[1.2,0.6],[1.6,0.8],[3.6,1.8]])
result = swarmCustom.add_particles_with_coordinates(swarmCoords)
swarmVariable = swarmCustom.add_variable(dataType='double', count=1)
swarmVariable.data[0] = 1.
swarmVariable.data[1] = 10.
swarmVariable.data[2] = 100.
swarmVariable.data[3] = 1000.
fig = glucifer.Figure(figsize=(800,400))
fig.append( glucifer.objects.Points( swarm=swarmCustom, pointSize=20,
fn_colour=swarmVariable, colourBar = True,
colours="red green blue", logScale=True) )
fig.append( glucifer.objects.Mesh(mesh, opacity=0.25))
fig.show()
Shapes with particle swarms¶
This example will demonstrate the creation of geometry using swarm
variables. Particle geometries often form means to define initial
fluid/material distribution in models. You will see later how once you
have defined a material identifier on each particle, you can then
delegate different behaviours based on this identifier using
Function objects.
First we’ll create a swarm, and then add a materialIndex swarm
variable to record the material type. We want to configure our
materialIndex such that it defines a circle. We will use a Python
loop to traverse all the particles, and for particles who’s coordinates
fall inside our circle shape, we will set their materialIndex value
to 1. Note that we have already initialised the materialIndex
values to 0 above, so we only need to modify the values of particles
within the circle.
import underworld as uw
import glucifer
import numpy as np
mesh = uw.mesh.FeMesh_Cartesian( elementRes = (64, 64),
maxCoord = (2., 1.) )
swarm = uw.swarm.Swarm( mesh=mesh )
# add a data variable which will store an index to determine material
materialIndex = swarm.add_variable( dataType="int", count=1 )
# populate our swarm across the mesh domain
swarmLayout = uw.swarm.layouts.PerCellSpaceFillerLayout( swarm=swarm, particlesPerCell=20 )
swarm.populate_using_layout( layout=swarmLayout )
materialIndex.data[:] = 0
# our circles parameters
circleRadius = 0.1
circleCentre = (1., 0.5)
# the particle loop
for index, coord in enumerate(swarm.particleCoordinates.data):
x = coord[0]
z = coord[1]
xx = x - circleCentre[0]
zz = z - circleCentre[1]
condition = (xx*xx + zz*zz < circleRadius**2)
if(condition == True): # inside the circle
materialIndex.data[index] = 1
# vis
fig = glucifer.Figure(figsize=(800,400))
fig.append( glucifer.objects.Points( swarm=swarm, fn_colour=materialIndex, colours='blue red',
colourBar=True, pointSize=2.0 ) )
fig.show()
More efficient ways to achieve this are presented in the Functions section of the user guide, though for operations that are usually only performed once in a simulation (such as initialisation), we recommend using whichever method you find most natural, as the costs are usually negligible.
Saving and load a swarms and swarm variables¶
As with other data types, the save() method is used to save
Swarm and SwarmVariable objects. Similarly, the load()
method is used to load Swarm and SwarmVariable objects. Note
that it is necessary to load the Swarm object before the
corresponding SwarmVariable object can be loaded. Ie, you cannot
load directly onto an existing SwarmVariable without first loading a
Swarm. This is to ensure that both objects are in a compatible state
(specifically, the correct number of particles exists and the data
ordering is identical).
Analogous to the process for mesh objects, where you wish to restart your simulation using saved data, you will first create an empty swarm, and then load the required swarm data. Similarly, once you have a swarm, you will create the required swarm variable, and then load its data.
# note that this cell requires the previous cell to be executed first.
ignore = swarm.save("SwarmWithCircle.h5")
ignore = materialIndex.save("SwarmWithCircle.materialIndex.h5")
swarmCopy = uw.swarm.Swarm( mesh=mesh )
swarmCopy.load("SwarmWithCircle.h5")
# The swarm geometry is now loaded but the swarmVariables are not restored.
# We have to explicitly create and re-populate any variables we have previously saved
materialIndexCopy = swarmCopy.add_variable("int",1)
materialIndexCopy.load("SwarmWithCircle.materialIndex.h5")
fig = glucifer.Figure(figsize=(800,400))
fig.append( glucifer.objects.Points( swarm=swarmCopy, fn_colour=materialIndexCopy, colours='blue red',
colourBar=True, pointSize=2.0 ) )
# Cleanup
if uw.rank()==0:
import os
os.remove( "SwarmWithCircle.h5" )
os.remove( "SwarmWithCircle.materialIndex.h5" )
fig.show()
Note that the above methods return SavedFileData objects which are
used required if you wish to create XDMF files and can be ignored
otherwise. Although it is not necessary to do so, we record them to the
ignore instance as they are not required here. See the Utilities
section of the user guide for further information on generating XDMF
files.
Functions¶
Function class objects provide the building blocks for mathematical
expression within Underworld2. The primary aim of this class is to
enable a natural description of mathematics through the Python syntax so
that users may quickly and accurately prototype model behaviour.
Functions are used extensively across the Underworld2 API and provide a
unified interface to Underworld2 discrete objects (SwarmVariable and
MeshVariable objects).
The user is encouraged to drill down interactively into the function submodules to discover available functionality, or alternatively browse the API reference reference at the underworld documentation site.
Overview:
- A simple example.
- Usage basics.
- Module overview.
- The
evaluate()method. - The
inputfunction. - Branching functions.
Keywords: functions, swarms, meshvariables, materials
A Simple Example¶
Let us define a function which we might use as a variable heat conductivity for a thermal problem. It will take the following temperature dependent form:
import underworld as uw
from underworld import function as fn
import glucifer
# first create a mesh and variable
mesh = uw.mesh.FeMesh_Cartesian(elementRes=(64,64), maxCoord=(2.,1.))
tempVar = mesh.add_variable(1)
# init the temp variable
for index,coord in enumerate(mesh.data):
tempVar.data[index] = coord[1]
# and now define the function.
fn_k = 5. + 8.*fn.math.exp(5.*tempVar) # a
# vis
fig = glucifer.Figure(figsize=(800,400))
fig.append(glucifer.objects.Surface(mesh,fn_k)) # b
fig.show() # c
Let’s deconstruct the Function defined at step (a) above:
fn_k = 5. # 1
+ 8.* # 2
fn.math.exp( # 3
5.*temp ) # 4
Things to note at the positions above:
- You can directly use Python native numerical objects. Under the hood,
the native Python float object created here ‘
5.’ will be automatically converted to an Underworld2Constanttype function. Note that arithmetic operations only currently support float type objects, and an exception will be thrown for other types. For this reason, you often need to be careful where you use Python natives (for instance, using ‘5’ here instead of ‘5.’ would result in an error). - Again, the native ‘
8.’ will be automatically converted. The addition operator here will be automatically converted to an UnderworldAdditionoperation through operator overloading. Likewise for the multiplication operation. - Note that for an exponential function, we need to use the Underworld
provided
fn.math.expfunction, not the Pythonmathmoduleexpfunction. - Here the argument (
5.*temp) is itself aFunction, andFunctioncompounding applies. Importantly, note that theMeshVariableis used directly in the arithmetic, and this is possible because it is also aFunctionclass object (more on this soon).
At step (b), we provide the function to the visualisation object
(Surface). Function objects are expected in many places across
the Underworld2 API. Finally, at step (c), the actually function
evaluation occurs, though under the hood it is a two step process. In
the first step, tests are performed to ensure that the provided function
is compatible with the required operation. Compatibility is dependent on
a number of factors which we discuss below (‘Conformal input & output
checking’), and if there are any problems an exception will be raised
here to notify the user. The second step is the required function
evaluation, and will generally be the most compuationally expensive
phase.
Usage basics¶
Underworld data objects¶
As seen in the example above, Underworld data objects (MeshVariable
and SwarmVariable types) may be used directly within functions, as
they are indeed themselves Function objects (via Python multiple
inheritance). The Function class provides a uniform interface to these
objects and is (largely) agnostic to the underlying data discretisation,
instead providing mechanisms for evaluation at arbitrary coordinates.
Note that true arbitrary coordinate evaluation is not possible for
SwarmVariable objects, as they are purely discrete and do not
(currently) have supporting interpolation functions. The special case
inputs provided by objects of the FunctionInput class may be used
for SwarmVariable objects (more on this soon).
The following simple example demonstrates querying a MeshVariable
object.
import underworld as uw
from underworld import function as fn
mesh = uw.mesh.FeMesh_Cartesian(elementRes=(64,64), maxCoord=(2.,1.))
sinvar = mesh.add_variable(1)
# initialise with a sine perturbation
import numpy as np
sinvar.data[:,0] = np.sin( np.pi*mesh.data[:,0] )
# use the `evaluate()` method to perform query.. more on this below
result = sinvar.evaluate( (0.25,0.) )
print("Evaluation result: {}".format(result))
import math
if not np.allclose( result, math.sqrt(2.)/2. ):
raise RuntimeError("Error! Expected result was not obtained.")
Evaluation result: [[ 0.70710678]]
Elementary algebraic operations¶
Use the Python equivalents (+,-,*,/) directly with
your Function objects! Function objects are operator overloaded
to facilitate this. Note however that only functions which return
floating point type values are compatible with elementary operations
currently, and an exception will be raised otherwise.
Let’s have a play with some mesh variables initialised to constants.
import underworld as uw
from underworld import function as fn
import numpy as np
mesh = uw.mesh.FeMesh_Cartesian(elementRes=(64,64), maxCoord=(2.,1.))
# create some more mesh variables and initialise
two_var = mesh.add_variable(1)
six_var = mesh.add_variable(1)
two_var.data[:] = 2.
six_var.data[:] = 6.
# create some functions via the Python operators
fn_plus = two_var + six_var
fn_minus = two_var - six_var
fn_div = two_var / six_var
fn_times = two_var * six_var
# check results.. evaluate anywhere as our mesh variables are constant
coord = (0.1234, 0.5678)
resultplus = fn_plus.evaluate( coord )
resultminus = fn_minus.evaluate( coord )
resultdiv = fn_div.evaluate( coord )
resulttimes = fn_times.evaluate( coord )
print("Addition result : {}".format(resultplus))
print("Subtraction result : {}".format(resultminus))
print("Division result : {}".format(resultdiv))
print("Multiplication result: {}".format(resulttimes))
# run tests
if not ( np.allclose(resultplus, 8.) and
np.allclose(resultminus,-4.) and
np.allclose(resultdiv, 1./3.) and
np.allclose(resulttimes, 12.) ):
raise RuntimeError("Error! Expected result was not obtained.")
Addition result : [[ 8.]]
Subtraction result : [[-4.]]
Division result : [[ 0.33333333]]
Multiplication result: [[ 12.]]
Convenience conversions¶
Python elementary types (int/floats/etc) may be used directly with
Underworld Function objects in algebraic operations. Likewise, Python
tuples (or lists) of Underworld Functions are automatically converted
into Function objects which return vector results composed of the
tuple/list entries. It is often important to remember that the
elementary algebraic operations are only compatible with float type
objects, and therefore you will often need to write ‘5.’ (which will
be converted to a double precision float object) instead of ‘5’
(which will be converted to a integer object).
Note that conversions can only occur automatically where the Python
object (tuple or primary math object) comes in contact with an
Underworld function, whereby conversion occurs by virtue of the object
overloading. However, occasionally we may wish to explicitly perform
conversions, generally where we wish to use the evaluate() method or
when we wish to utilise overloading. The convert() static method on
the Function class provides this functionality.
import underworld as uw
from underworld import function as fn
import numpy as np
import math
mesh = uw.mesh.FeMesh_Cartesian(elementRes=(64,64), maxCoord=(2.,1.))
# create a new mesh variable and initialise
cosvar = mesh.add_variable(1)
cosvar.data[:,0] = np.cos( np.pi*mesh.data[:,0] )
sinvar = mesh.add_variable(1)
sinvar.data[:,0] = np.sin( np.pi*mesh.data[:,0] )
# create function. note that `2.` and the `1.` are automatically converted
sin2 = 2.*sinvar*cosvar + 1.
# evaluate somewhere and then test
import random
coord = ( random.random(), 0. )
result = sin2.evaluate(coord)
expected = math.sin(2.*math.pi*coord[0]) + 1. # via double angle formula
if not np.allclose(result, expected, rtol=5e-2, atol=5e-2):
print(coord, expected, result)
raise RuntimeError("Error! Expected result was not obtained.")
# also let's create a vector function on the fly.
# first create a vector mesh variable
vecvar = mesh.add_variable(2)
vecvar.data[:] = (1.,1.)
# now the function definition, evaluation, and test
fn_vec = vecvar + (sinvar, cosvar)
result = fn_vec.evaluate(coord)
expected = (math.sin(math.pi*coord[0]) + 1., math.cos(math.pi*coord[0]) + 1.)
if not np.allclose(result, expected, rtol=5e-2, atol=5e-2):
print(result, expected)
raise RuntimeError("Error! Expected result was not obtained.")
# now let's create a Python tuple from UW functions.
vec_as_py_tuple = (sinvar, cosvar)
print("`vec_as_py_tuple` type is: {}".format(type(vec_as_py_tuple)))
# The following will not work! This is because it is still a Python
# native `tuple` object, and as such does not have an `evaluate()` method.
# vec_as_py_tuple.evaluate()
# However, we can convert it to a UW function, and then use `evaulate()`.
vec_as_uw_fn = fn.Function.convert(vec_as_py_tuple)
print("`vec_as_uw_fn` type is: {}".format(type(vec_as_uw_fn)))
# this is better
print("evaluate: {}".format(vec_as_uw_fn.evaluate(coord)))
vec_as_py_tuple type is: <type 'tuple'> vec_as_uw_fn type is: <class 'underworld.function._function.add'> evaluate: [[ 0.91933336 0.39229551]]
Basic mathematical functions¶
Basic functions (such as sin() and exp()) are provided by the
underworld.function.math module. Note that the Python math
module is not compatible with Underworld2 functions (you must use
our math module). Operator overloads are also provided from the
indexing operator ([]) and the power operator (**).
We will construct some more double angle formula here. While in the
previous example, we used Numpy to initialise a mesh variable object to
construct sin and cos like functions, here all mathematical
operations will be performed by Underworld Function objects.
import underworld as uw
from underworld import function as fn
import random
import numpy as np
import math
# trig funcs
sin = fn.math.sin()
cos = fn.math.cos()
tan = fn.math.tan()
# double angle formula
sin_2theta = 2.*sin*cos
cos_2theta = 1. - 2.*sin**2
tan_2theta = (2.*tan)/(1.-tan**2)
# get somewhere to evaluate
theta = random.random()
# do things check out?
if not ( np.allclose( sin_2theta.evaluate(theta), math.sin(2*theta) ) and
np.allclose( cos_2theta.evaluate(theta), math.cos(2*theta) ) and
np.allclose( tan_2theta.evaluate(theta), math.tan(2*theta) ) ):
raise RuntimeError("Error! Expected result was not obtained.")
Relational and logical functions¶
Relational functions are constructed via the Python relational operators
(<,<=,>,>=). Underworld functions for AND,
OR and XOR logical operations are also available, and these
overload the Python bitwise operators (&,|,^) (though
they do not perform bitwise operations). These functions will all
return boolean results.
import underworld as uw
from underworld import function as fn
import numpy as np
import math
import glucifer
mesh = uw.mesh.FeMesh_Cartesian(elementRes=(64,64), maxCoord=(2.,1.))
# define a logical function for inside a unit circle.
# we will use the `input()` function as a proxy for the coordinate (see below).
radius = 1.
coord = fn.input()
inside_circle = ( coord[0]**2 + coord[1]**2 < radius**2 )
# test at some locations
if not ( np.allclose( inside_circle.evaluate( (0. ,0. ) ), True ) and
np.allclose( inside_circle.evaluate( (1. ,1. ) ), False ) and
np.allclose( inside_circle.evaluate( (0.5,0.5) ), True ) and
np.allclose( inside_circle.evaluate( (0.9,0.5) ), False ) ):
raise RuntimeError("Error! Expected result was not obtained.")
# now something a bit more complex.. first, some circles
radius = 0.25
offset = (1.0,0.5)
deltax = radius*math.cos(math.pi/4.)
deltay = radius*math.sin(math.pi/4.)
coord1 = fn.input() - offset
circle1 = ( fn.math.dot(coord1,coord1) < radius**2 )
coord2 = coord1 - ( deltax, deltay)
circle2 = ( fn.math.dot(coord2,coord2) < radius**2 )
coord3 = coord1 - (-deltax, deltay)
circle3 = ( fn.math.dot(coord3,coord3) < radius**2 )
coord4 = coord1 - ( deltax,-deltay)
circle4 = ( fn.math.dot(coord4,coord4) < radius**2 )
coord5 = coord1 - (-deltax,-deltay)
circle5 = ( fn.math.dot(coord5,coord5) < radius**2 )
# now create a cross.. note the use of the OR operator
cross = ( (fn.math.abs(coord1[0])<radius/2.) | (fn.math.abs(coord1[1])<radius/2.) )
# visualise the XOR of these shapes
fig = glucifer.Figure(figsize=(800,400))
fig.append(glucifer.objects.Surface(mesh, circle1^circle2^circle3^circle4^circle5^cross, resolution=600, onMesh=False, colours=['white','blue']))
fig.show()
Conformal input & output checking¶
When you define your functions, they are input agnostic in the sense that you describe your function without declaring the type of argument that will be used when the function is eventually evaluated. At evaluation time, checks are performed to ensure the argument is compatible with the provided function. Likewise, the output returned by the evaluated function will be checked to ensure it is of the required form. If a check fails, an exception is thrown.
Incompatible function input¶
The example below demonstrates an incompatible input type. We will use
the tempVar object, and attempt to evaluate it with a scalar input.
This is not a valid operation, as tempVar is a MeshVariable type
object and therefore can only be successfully evaluated at a valid
domain coordinate, and therefore requires a vector input (the
coordinate!). We will also use the min() function to set a lower
bound (of zero) on the function’s returned results.
import underworld as uw
from underworld import function as fn
mesh = uw.mesh.FeMesh_Cartesian()
tempVar = mesh.add_variable(1)
positiveTemp = fn.misc.min(0.,tempVar)
try:
positiveTemp.evaluate((0.1,))
except RuntimeError as e:
print("RuntimeError: "+str(e))
RuntimeError: Issue utilising function of class 'MeshVariable' constructed at:
Line 4 of notebook cell 7:
tempVar = mesh.add_variable(1)
Error message:
Function input dimensionality (1) does not appear to match mesh variable dimensionality (2).
Let’s deconstruct this error message. The second part specifies the
actual problem that was encountered, in this case the input being a
scalar where a vector was required. The first part tells you where the
function that cannot be evaluated was defined. This is useful
information because often you will define a function in a very different
place to where you will eventually evaluate it. More importantly,
complex functions are usually defined through numerous Python calls (as
in this case where we introduce min()), so we try to flag to you the
actual subfunction that was problematic (in this case tempVar, not
min).
Functions can of course only be evaluated at a point within their
domain. So for MeshVariable objects, you may only evaluate at a
coordinate within the mesh domain. Furthermore, in parallel, you may
only evaluate within the mesh domain local to the process:
# note that this cell requires you to run the previous cell first
try:
positiveTemp.evaluate((-1.0,1.0)) # but the failure will occur at evaluation time.
except ValueError as e:
print("RuntimeError: "+str(e))
RuntimeError: Issue utilising function of class 'MeshVariable' constructed at:
Line 4 of notebook cell 7:
tempVar = mesh.add_variable(1)
Error message:
FeVariable interpolation at location (-1, 1) does not appear to be valid.
Location is probably outside local domain.
SwarmVariable objects are purely discrete (they do not carry
interpolation functions as do their MeshVariable counterpart), so
their domain is the set of all particles coordinates. As such, an
arbitrary coordinate cannot be used as an input. However due to
algorithmic considerations, it is also not possible to feed in the
coordinate of a particle for evaluation, and instead you must use the
special FunctionInput class for evaluations, which we will explore
shortly.
import underworld as uw
from underworld import function as fn
mesh = uw.mesh.FeMesh_Cartesian()
swarm = uw.swarm.Swarm(mesh)
svar = swarm.add_variable('double',1)
swarm.populate_using_layout(uw.swarm.layouts.PerCellSpaceFillerLayout(swarm,1))
try:
svar.evaluate((0.2,0.3))
except RuntimeError as e:
print("RuntimeError: "+str(e))
RuntimeError: Issue utilising function of class 'SwarmVariable' constructed at:
Line 5 of notebook cell 9:
svar = swarm.add_variable('double',1)
Error message:
Unable to evaluate function using provided input.
Note that your function is constructed using a SwarmVariable. We do not
currently support interpolation of SwarmVariables, so where you wish to
evaluate() a SwarmVariable based function, you may only do so using
the swarm itself as the input to the function evaluation, resulting in
evaluation at each particle location. Note also that for evaluation,
you must use the Swarm object from which the SwarmVariable was derived.
Alternatively, you may wish to project your function onto a MeshVariable
(using a MeshVariable_Projection object) and then perform interpolation
on the MeshVariable object.
Incompatible function output¶
All functions require inputs of a certain form, and also return results
of a (generally different) form. The object that is performing the
evaluation then utilises these outputs for its given task. In the
following example, we are trying to visualise a function using the
Surface visualisation object. This object generates a raster image
representing the provided function, with a single colour mapped to a
given scalar function result. As such, while the function is
successfully evaluated in the following, it returns a vector result
where a scalar is required.
# the `Surface` object can only visualise scalar objects
import underworld as uw
from underworld import function as fn
import glucifer
mesh = uw.mesh.FeMesh_Cartesian()
vecvar = mesh.add_variable(2)
fig = glucifer.Figure()
try:
fig.append(glucifer.objects.Surface(mesh,vecvar)) # failure will occur here
except RuntimeError as e:
print("RuntimeError: "+str(e))
RuntimeError: Provided function must return a scalar result.
Function are dynamic¶
Functions are dynamic in the sense that they do not capture a static representation of their constituent functions at definition time, but instead simply refer to these functions at evaluation time. This means that if any of the constituent functions change (such as when you update your velocity field), these changes are reflected immediately:
import underworld as uw from underworld import function as fn mesh = uw.mesh.FeMesh_Cartesian() # create mesh var and init to 1. meshvar = mesh.add_variable(1) meshvar.data[:] = 1. # now create a function which utilises the mesh var and evaluate evaluation_coord = (0.1,0.2) fn_test = meshvar * 5. result_before = fn_test.evaluate(evaluation_coord) # now modify the meshvar and evaluate the fn again meshvar.data[:] = -1. result_after = fn_test.evaluate(evaluation_coord) print(result_before,result_after)
(array([[ 5.]]), array([[-5.]]))
The evaluate() Method¶
Once you have created functions, you will pass these into various
objects within the Underworld2 API, however you will often also wish to
directly evaluate the functions you have created for analytic or testing
purposes. The evaluate() method provides this capability. Here are
some basic examples using the function created previously:
import underworld as uw
from underworld import function as fn
coord = fn.input() # more on the input function soon
fn_k = 5. + 8.*fn.math.exp(5.*coord[1])
# evaluate at single location, provide as a coordinate tuple or list.
result = fn_k.evaluate( (0.3,0.2) )
print( "Single evaluation result = \n{}\n".format(result))
# evaluate at a set of locations, provide these as a numpy array.
import numpy as np
count = 10
# create an empty array
locations = np.zeros( (count,2) )
# specify evaluation coodinates
locations[:,0] = 0.5
locations[:,1] = np.linspace(0.,1.,count)
# evaluate
result = fn_k.evaluate(locations)
print( "Multi evaluation result = \n{}".format(result))
Single evaluation result =
[[ 26.74625463]]
Multi evaluation result =
[[ 13. ]
[ 18.94327199]
[ 29.30185422]
[ 47.3559204 ]
[ 78.82251482]
[ 133.66592538]
[ 229.25299916]
[ 395.8525702 ]
[ 686.22046174]
[ 1192.30527282]]
The FunctionInput class¶
A further type of input to function evaluation are FunctionInput
class objects. These are shortcuts to their underlying data, but
leverage the fundamental object data for higher efficiency evaluations,
and sometimes as a necessity to facilitate the evaluation. Currently
provided FunctionInput classes include:
underworld.mesh.FeMesh: Evaluation at all mesh vertex coordinates.underworld.mesh.FeMesh_IndexSet: Evaluation at the mesh vertex coordinates within the set.underworld.swarm.Swarm: Evaluation at all swarm particle coordinates.underworld.swarm.VoronoiIntegrationSwarm: Evaluation at all voronoi swarm particle coordinates.
The above behave as FunctionInput classes by way of the multiple
inheritence mechanisms of Python. For example, mesh objects (ie, objects
of class FeMesh) are also FunctionInput objects.
Note that you will get identical evaluation results using the entire mesh vertex Numpy array as evaluation input, or using the mesh directly as an input:
# note that this cell requires the previous cell to have been executed
mesh = uw.mesh.FeMesh_Cartesian(elementRes=(64,64), maxCoord=(2.,1.))
# first lets confirm that mesh is indeed a `FunctionInput` object
if not isinstance(mesh, fn.FunctionInput):
raise RuntimeError("Error! The mesh does not appear to be an instance of the `FunctionInput` class.")
# now evaluate at all mesh vertices:
results_using_functioninput = fn_k.evaluate(mesh)
# likewise, let's do the numpy equivalent
results_using_numpy = fn_k.evaluate(mesh.data)
# confirm identical results
if not np.allclose( results_using_functioninput, results_using_numpy ):
raise RuntimeError("Error! Results differ where they should be the same.")
# Note that the `FeMesh_IndexSet` which contains all mesh vertices should *also* return identical results:
# let's create the 'empty' set
allindices = mesh.specialSets['Empty']
# now invert to obtain the 'full' set
allindices.invert()
# evaluate
results_using_indexset = fn_k.evaluate(allindices)
# again confirm identical results
if not np.allclose( results_using_functioninput, results_using_indexset ):
raise RuntimeError("Error! Results differ where they should be the same.")
While all the above methods yield the same numerical results, note that
there are efficiency differences. When you provide a Numpy array as
input to a MeshVariable object evaluation, it is first necessary to
determine which element the evaluation coordinate belongs to, which can
be an expensive operation (in particular for deformed mesh). Once the
owning element is determined, the interpolation itself needs to be
calculated. However, when using the mesh object directly as a function
input, the evaluation leverages higher level object information to
directly extract the nodal values (ie, owning element and interpolation
are not required). Consider the following timing results:
import underworld as uw mesh = uw.mesh.FeMesh_Cartesian(elementRes=(64,64), maxCoord=(2.,1.)) var = mesh.add_variable(1) # first evaluate directly using mesh %timeit var.evaluate(mesh)
1000 loops, best of 3: 650 µs per loop
import underworld as uw mesh = uw.mesh.FeMesh_Cartesian(elementRes=(64,64), maxCoord=(2.,1.)) var = mesh.add_variable(1) # now evaluate via numpy array %timeit var.evaluate(mesh.data)
100 loops, best of 3: 3.43 ms per loop
import underworld as uw
mesh = uw.mesh.FeMesh_Cartesian(elementRes=(64,64), maxCoord=(2.,1.))
var = mesh.add_variable(1)
# now flag the mesh as deformed. note that we are not
# moving any mesh vertices, but simply flagging the
# mesh as deformed.
with mesh.deform_mesh():
pass
%timeit var.evaluate(mesh.data)
10 loops, best of 3: 80.6 ms per loop
So we see that there isn’t a dramatic difference until the mesh is deformed. Note also that as we haven’t truly deformed the mesh, the observed result will generally be the best case scenario.
The other important FunctionInput objects are the swarm based items.
In the case of functions which utilise SwarmVariable objects, these
are the only option for evaluation from within Python. This is
because particles are inherently discrete, and therefore evaluation at
arbitrary locations is not directly possible. In some case where you
supply SwarmVariable functions to be used within Underworld
operations (such as a viscosity for the Stokes system), a nearest
neighbour calculation is implicitly utilised, however these mechanisms
are not (yet) exposed for calls to the evaluate() method. Future
releases will provide greater functionality for operating over
SwarmVariable functions. Again, simply pass the swarm to the
evaluate function to utilse the swarm FunctionInput behaviour:
import underworld as uw
from underworld import function as fn
import numpy as np
fn_k = 5. + 8.*fn.math.exp(5.*fn.input()[1])
mesh = uw.mesh.FeMesh_Cartesian(elementRes=(64,64), maxCoord=(2.,1.))
swarm = uw.swarm.Swarm(mesh)
swarm.populate_using_layout( uw.swarm.layouts.PerCellSpaceFillerLayout(swarm,20) )
# add variable to use for conductivity
k_var = swarm.add_variable('double',1)
# now initialise these particles to have values which correspond to those
# provided by fn_k. we can use either `swarm` or `swarm.data` as inputs
# here and the results in this case will be identical.
k_var.data[:] = fn_k.evaluate(swarm)
# note that k_var is itself of the `Function` class
if not isinstance(k_var, fn.Function):
raise RuntimeError("Error! Swarm variable does not appear to be of `Function` class")
# now try to evaluate at arbitrary location
# k_var.evaluate( (0.3,0.2) ) # this will raise a RuntimeError!
# now evaluate using swarm
results_using_swarm_k_var = k_var.evaluate( swarm )
# also evaluate fn_k using swarm
results_using_swarm_fn_k = fn_k.evaluate( swarm )
if not np.allclose( results_using_swarm_fn_k, results_using_swarm_fn_k ):
raise RuntimeError("Error! These arrays should have identical values.")
The input function¶
The function.input class simply provides the identity function. It
returns whatever values are passed to it!
import underworld as uw
from underworld import function as fn
import numpy as np
mesh = uw.mesh.FeMesh_Cartesian(elementRes=(64,64), maxCoord=(2.,1.))
# create the input class/function
fn_input = fn.input()
# confirm behaviour
print(fn_input.evaluate( 2. ) )
print(fn_input.evaluate( ((0.2), (0.3)) ))
print(fn_input.evaluate( ((0.2), (0.3), (0.4)) ))
print(fn_input.evaluate( mesh.data[0:5] ))
if not np.allclose(mesh.data, fn_input.evaluate(mesh.data)):
raise RuntimeError("Error! These arrays should have identical values.")
[[ 2.]]
[[ 0.2 0.3]]
[[ 0.2 0.3 0.4]]
[[ 0. 0. ]
[ 0.03125 0. ]
[ 0.0625 0. ]
[ 0.09375 0. ]
[ 0.125 0. ]]
The input function is most often used when you wish to construct
functions which operate over a particular coordinate axis:
import underworld as uw from underworld import function as fn import glucifer import numpy as np fn_input = fn.input() mesh = uw.mesh.FeMesh_Cartesian(elementRes=(64,64), maxCoord=(2.,1.)) # Create sinusoidal function in the horizontal direction. # Note that the square brackets operator is overloaded to extract # a certain axis' value from vector functions. fn_sin = fn.math.sin(8.*np.pi*fn_input[0]) # take a look fig = glucifer.Figure(figsize=(800,400)) fig.append(glucifer.objects.Surface(mesh,fn_sin)) fig.show()
We also provide an alias to this class in function.coord. This is
strictly a Python alias, and is identical in all ways (except name!) to
function.input. It is provided because we often use the input
function to operate on coordinate values (as above), but it is worth
noting that function inputs are certainly not restricted to coordinates
or vectors, and in these cases it would be misleading to use the
coord function.
Branching Functions¶
Branching functions are functions which nest other functions, and select which function to execute based on some condition (also expressed as a function!). The most common use case for branching functions are for defining materials within your model which exhibit different behaviours (based on their associated function). Note that Underworld2 does not include an explict ‘material’ construct (unlike Underworld1), however the equivalent functionality may be obtained through branching functions.
Two branching functions are currently provided:
function.branching.map: Key/value type mapping of behavior.
function.branching.conditional: if/elif type behaviour.
Refer to respective API documentation for further information on these
function classes. The map function is usually used to construct
material type behaviours. Note that the map function provides a
subset of the conditional function behaviour, though it has
performance advantages.
Let us construct a basic model with material type behaviours:
import underworld as uw
from underworld import function as fn
import glucifer
import numpy as np
mesh = uw.mesh.FeMesh_Cartesian(elementRes=(64,32),minCoord=(-2.0, -1.0), maxCoord=(2.0, 1.0))
swarm = uw.swarm.Swarm(mesh)
swarm.populate_using_layout(uw.swarm.layouts.PerCellSpaceFillerLayout(swarm,80))
# add a variable which will act as the key function for the map function.
material_index = swarm.add_variable("int",1)
# add the following for convenience and clarity.
outside_circle = 0
inside_circle = 1
# initialise material_index to be outside_circle everywhere
material_index.data[:] = outside_circle
# now set to inside_circle where inside circle!
r2 = 0.8
for index, position in enumerate(swarm.particleCoordinates.data):
if position[0]**2 + position[1]**2 < r2:
material_index.data[index] = inside_circle
# create mapped behaviour dictionary
fn_z = fn.input()[1]
fn_sin = fn.math.sin(8.*np.pi*fn.input()[0])
map_dict = { outside_circle:fn_z,
inside_circle:fn_sin }
# create function
fn_map = fn.branching.map( fn_key = material_index,
mapping = map_dict )
# viz
fig = glucifer.Figure(figsize=(800,400))
fig.append( glucifer.objects.Points(swarm,fn_map, pointSize=1.5) )
fig.show()
Alternatively, we can achieve an identical result using the
conditional function, and its evaluate() method:
# note that the previous cell needs to be executed before this.
# the 'position' variable (created above) retains a reference to the
# particleCoordinates numpy array. this inteferes with adding a new
# swarm variable, so we must delete it.
del position
# add a new variable first (we will test later for equality)
material_index2 = swarm.add_variable("int",1)
# first define a function which returns true if inside the circle
coord = fn.input()
fn_in_circle = fn.math.dot(coord,coord) < r2 # note the overload of the '<' operator
# now create the conditional
fn_conditional = fn.branching.conditional( ( (fn_in_circle, inside_circle),
( True, outside_circle) ) )
# use the evaluate, writing the results out to the material_index2
material_index2.data[:] = fn_conditional.evaluate(swarm)
# check to ensure that we have identical results
if not (material_index2.data == material_index.data).all():
raise RuntimeError("Error! These arrays should have identical values.")
Systems¶
Numerous dynamic systems are implemented in Underworld. They may be
found within the systems submodule. For specific information on the
different system classes, the user is directed to the API reference
available at the Underworld documentation site:
http://underworld2.readthedocs.io/
We will consider here the basic workflow for creating and configuring an Underworld numerical system. The general process is as follows:
- Create your mesh.
- Create any required field(s) on the mesh (such as a temperature field).
- Create any required boundary condition objects.
- Create function objects to define any required physical quantities.
- Create system.
- Create solver and solve system.
Note that for systems that are solved many times (for perhaps time-stepping), you will generally only create the system (and solver) once, and then solve it numerous times.
Boundary Conditions¶
Boundary conditions form part of your known information for the system you are creating. In Underworld, these knowns are either of the Dirichlet (fixed value) or Neumann (fixed gradient) type, and are applied on a per node and per degree of freedom (DOF) basis. Note also that the entire domain boundary must be piecemeal either Dirichlet or Neumann type, but no section of boundary can be both Dirichlet and Neumann. Sections of the boundary that do not have any BCs explicitly set by the user will implicitly be of Neumann type (with zero gradient).
To give this discussion a mathematical grounding, let’s consider the heat equation. The strong form of the boundary value problem is
where, \(\alpha\) is the diffusivity, \(u\) is the temperature, \(f\) is a source term, \(g\) is the Dirichlet condition, and \(h\) is a Neumann condition. The problem boundary, \(\Gamma\), admits the decomposition \(\Gamma=\Gamma_g\cup\Gamma_h\) where \(\emptyset=\Gamma_g\cap\Gamma_h\).
The Dirichlet Condition¶
The Dirichlet condition part of the problem specification requires that
our solution temperature \(u\) is equal to the user specified
function \(g\) across the user specified boundary \(\Gamma_g\).
In Underworld, the Dirichlet boundary \(\Gamma_g\) is specified
through the use of IndexSet objects. These objects simply contain
lists of mesh nodes (or rather the node indices) which the user wishes
to flag as boundary nodes. A number of commonly used IndexSet
objects are provided directly through the mesh object’s specialSets
dictionary:
import underworld as uw mesh = uw.mesh.FeMesh_Cartesian() mesh.specialSets.keys()
['MaxI_VertexSet',
'Top_VertexSet',
'Left_VertexSet',
'MinI_VertexSet',
'AllWalls_VertexSet',
'Bottom_VertexSet',
'Right_VertexSet',
'MinJ_VertexSet',
'MaxJ_VertexSet',
'Empty']
IndexSet objects are discussed further in the mesh section of this
user guide. Importantly, note that the specialSets objects are
simply provided for convenience, and the user is indeed free to create
any set they wish. As an example, we can easily reproduce a set provided
by the specialSets dictionary by explicitly adding the required
indices:
import underworld as uw
import numpy as np
mesh = uw.mesh.FeMesh_Cartesian()
# first option, use provided convenience set.
leftwall = mesh.specialSets['Left_VertexSet']
# second option, create empty set and add required indices.
leftwall_handcrafted = mesh.specialSets['Empty']
for index,coord in enumerate(mesh.data):
# grab only left wall indices
if coord[0]<0.01:
leftwall_handcrafted.add(index)
# confirm sets are identical
if len(leftwall) != len(leftwall_handcrafted):
raise RuntimeError("Sets should have identical size!")
if not np.allclose(leftwall,leftwall_handcrafted):
raise RuntimeError("Sets should contain identical values!")
For the Dirichlet condition, the user specified function \(g\) which
provides the required values along the Dirichlet boundary is implicitly
set by directly modifying the solution MeshVariable at the required
nodes. Note that you can set or modify your required \(g\) function
at any stage before the solve operation.
Let’s put all these pieces together to generate a basic heat equation solution. For this system, we’ll set the bottom boundary to have a temperature of \(u=373\), and the top boundary to have a temperature \(u=273\). We won’t set anything for the side walls, but note that this will implicitly yield a zero temperature gradient (ie heat flux) Neumann condition.
import underworld as uw
import glucifer
# 1. Create a mesh.
mesh = uw.mesh.FeMesh_Cartesian(maxCoord=(2.,1.))
# 2. Create the required field(s).
tempField = mesh.add_variable(1)
# 3. Create required boundary condition object.
# 3.a First create the required IndexSet objects to specify the required boundary.
# Note shorthand addition operation.
botSet = mesh.specialSets['Bottom_VertexSet']
topSet = mesh.specialSets[ 'Top_VertexSet']
totSet = botSet + topSet
# 3.b Now we can create the Dirichlet condition object.
condition = uw.conditions.DirichletCondition(tempField,totSet)
# 3.c Let's set the g function now by directly modifying the temperatureField object.
# Note that this operation can be done at any stage!
# Note also the shorthand we use to set the required values.
tempField.data[botSet] = 373.
tempField.data[topSet] = 273.
# 4. Let's go ahead and set a simple diffusivity.
alpha = 1.
# 5. Create the required system.
heatSys = uw.systems.SteadyStateHeat(temperatureField=tempField,
fn_diffusivity=alpha,
conditions=condition)
# 6. Create solver and solve.
solver = uw.systems.Solver(heatSys)
solver.solve()
# vis.
fig = glucifer.Figure(figsize=(800,400))
fig.append(glucifer.objects.Surface(mesh, tempField, colours="cubelaw2"))
fig.show()
Note in the above visualisation that we have satisfied our model requirements:
- User specified Dirichlet condition is satisfied \((u=373 \text{ on } \Gamma_{\text{bottom}},\,\, u=273 \text{ on } \Gamma_{\text{top}})\).
- Inspection suggests side boundaries have zero normal temperature gradient, but the example in the Neumann section below explicitly sets the zero flux Neumann condition and yields identical results.
As mentioned, the user is free to flag any mesh nodes as knowns of the
system, including nodes internal to the domain. They are also free to
change their \(g\) specification (ie, the values set on the
MeshVariable boundaries) at any point, however if they wish to
change the boundary specification (\(\Gamma_g\)), it is current
necessary to create a new system object. Note also that while the user
may freely set values on any node of the MeshVariable object, for
instantaneous systems only values set on BC-flagged nodes will be
retained after the solve() operation, with the solver then delegated
to determine all other nodal values. We will demonstrate these
behaviours in the following example:
import underworld as uw
import glucifer
res=32
mesh = uw.mesh.FeMesh_Cartesian(elementRes=(res,res),minCoord=(-0.5,-0.5), maxCoord=(0.5,0.5))
tempField = mesh.add_variable(1)
# create bc sets. first, allwalls.
allWallSet = mesh.specialSets['AllWalls_VertexSet']
# now grab central patch
centralSet = mesh.specialSets['Empty']
for index,coord in enumerate(mesh.data):
if max(coord*coord)<0.03:
centralSet.add(index)
condition = uw.conditions.DirichletCondition(tempField,centralSet+allWallSet)
tempField.data[allWallSet] = 0.
tempField.data[centralSet] = 1.
# for fun, let's also set a random non-bc node to a random value.
# the value will be replaced at solve time by the solver.
tempField.data[res+2] = -999999.
heatSys = uw.systems.SteadyStateHeat(temperatureField=tempField,
fn_diffusivity=1.,
conditions=condition)
solver = uw.systems.Solver(heatSys)
solver.solve()
# grab a copy of solved data for later use, and then let's
# flip the config and re-solve immediately
tempField_orig = tempField.copy(deepcopy=True)
tempField.data[allWallSet] = 1.
tempField.data[centralSet] = 0.
solver.solve()
# vis
fig1 = glucifer.Figure(figsize=(400,400), title="original")
fig1.append(glucifer.objects.Surface(mesh, tempField_orig, colours="ocean", colourBar=False, onMesh=False))
fig2 = glucifer.Figure(figsize=(400,400), title="flipped")
fig2.append(glucifer.objects.Surface(mesh, tempField , colours="ocean", colourBar=False, onMesh=False))
glucifer.Figure.show_grid([fig1,fig2])
 |  |
So far we have only considered scalar systems having only a single
degree of freedom per mesh node. For the Stokes system however we solve
for a velocity unknown having 2 or 3 DOFs per mesh node (for 2D and 3D
simulations respectively). While the general process is unchanged, you
will now be required to specify boundary conditions on a per-DOF basis.
The way this is handled in Underworld is to have the user provide an
IndexSet object for each DOF. So for the previous examples, only a
single IndexSet was required for the single DOF of the temperature
field, while for a velocity field we require 2 (or 3) IndexSet
objects, one for the \(V_x\) DOF, and another for the \(V_y\)
DOF (and a further for the \(V_z\) DOF in 3D).
Let’s consider the basic no-slip configuration first:
import underworld as uw
import glucifer
# 1. Create a mesh.
mesh = uw.mesh.FeMesh_Cartesian(elementRes=(32,32))
# 2. Create the required field(s).
velField = mesh.add_variable(2)
preField = mesh.subMesh.add_variable(1)
# 3. Create required boundary condition object.
allWalls = mesh.specialSets['AllWalls_VertexSet']
# Note that we now pass in a tuple of IndexSet objects, representing
# the nodes to flag for (V_x, V_y). For no-slip, we flag both DOFs for
# each node, so therefore the index sets must be identical.
# Note also that we generally do not set any conditions on the pressure
# unknowns.
condition = uw.conditions.DirichletCondition(velField, (allWalls,allWalls) )
# Init the velocity field to zero everywhere. Therefore g=(0,0).
velField.data[:] = (0.,0.)
# 4. Let's create a simple function for the viscosity & bodyforce.
viscosity = 1.
coord = uw.function.coord()
bodyforce = (0.,-coord[0])
# 5. Create the required system.
stokesSys = uw.systems.Stokes(velocityField=velField,
pressureField=preField,
fn_viscosity=viscosity,
fn_bodyforce=bodyforce,
conditions=condition )
# 6. Create solver and solve
solver = uw.systems.Solver(stokesSys)
solver.solve()
# Grab copy of solution.
noslipvel = velField.copy(deepcopy=True)
# As before, we can modify g at any point.
# Create a mobile lid config:
topWall = mesh.specialSets['Top_VertexSet']
velField.data[topWall] = (1.,0.)
# also drop body force
stokesSys.fn_bodyforce = (0.,0.)
solver.solve()
# vis
fignoslip = glucifer.Figure(figsize=(400,400), title="no slip")
fignoslip.append(glucifer.objects.VectorArrows(mesh, noslipvel))
fignoslip.append(glucifer.objects.Surface(mesh, uw.function.math.dot(noslipvel,noslipvel), colours="isorainbow"))
figtraction = glucifer.Figure(figsize=(400,400), title="mobile lid")
figtraction.append(glucifer.objects.VectorArrows(mesh, velField))
figtraction.append(glucifer.objects.Surface(mesh, uw.function.math.dot(velField,velField), colours="isorainbow"))
glucifer.Figure.show_grid([fignoslip,figtraction])
 |  |
The setup necessary to achieve a free-slip configuration is more complex. In this scenario, we wish to treat velocity components normal to the boundaries as knowns of the system, tethering this component to zero if we wish to prevent in/outflows. The tangential component should be left free to allow slip. Therefore, for a basic free-slip configuration on all walls of the standard cartesian geometry, we require:
- Left/right wall nodes: \(V_x = 0\), \(V_y\) free.
- Top/bottom wall nodes: \(V_x\) free, \(V_y = 0\).
This translates into flagging all the \(V_x\) nodes on the left/right wall, and all the \(V_y\) nodes on the top/bottom walls:
import underworld as uw
import glucifer
# 1. Create a mesh.
mesh = uw.mesh.FeMesh_Cartesian(elementRes=(32,32))
# 2. Create the required field(s).
velField = mesh.add_variable(2)
preField = mesh.subMesh.add_variable(1)
# 3. Create required boundary condition object.
leftrightWalls = mesh.specialSets['Left_VertexSet'] + mesh.specialSets['Right_VertexSet']
topbottWalls = mesh.specialSets[ 'Top_VertexSet'] + mesh.specialSets['Bottom_VertexSet']
condition = uw.conditions.DirichletCondition(velField, (leftrightWalls,topbottWalls) )
# Init the velocity field to zero everywhere.
velField.data[:] = (0.,0.)
# 4. Let's create a simple function for the viscosity & bodyforce.
viscosity = 1.
coord = uw.function.coord()
bodyforce = (0.,-coord[0])
# 5. Create the required system.
stokesSys = uw.systems.Stokes(velocityField=velField,
pressureField=preField,
fn_viscosity=viscosity,
fn_bodyforce=bodyforce,
conditions=condition )
# 6. Create solver and solve
solver = uw.systems.Solver(stokesSys)
solver.solve()
# Grab copy of solution.
noslipvel = velField.copy(deepcopy=True)
# As before, we can modify g at any point.
# Create inflow/outflow condition
for index in leftrightWalls:
if mesh.data[index][1] < 0.3:
velField.data[index] = (1.,0.)
# ramp up body force
stokesSys.fn_bodyforce = 100.*stokesSys.fn_bodyforce
solver.solve()
# vis
figfreeslip = glucifer.Figure(figsize=(400,400), title="free slip")
figfreeslip.append(glucifer.objects.VectorArrows(mesh, noslipvel))
figfreeslip.append(glucifer.objects.Surface(mesh, uw.function.math.dot(noslipvel,noslipvel), colours="gebco"))
figinflow = glucifer.Figure(figsize=(400,400), title="free slip with in/out flow")
figinflow.append(glucifer.objects.VectorArrows(mesh, velField))
figinflow.append(glucifer.objects.Surface(mesh, uw.function.math.dot(velField,velField), colours="gebco"))
glucifer.Figure.show_grid([figfreeslip,figinflow])
 |  |
Note in the second example above that while you are free to create any
inflow/outflow condition, unless you have otherwise specified, the
system will be configured for incompessibility, and therefore any inflow
must exactly balance outflow. For complex geometries, ensuring this
balance is not always possible given the discrete nature of the models,
and you might instead consider allowing some part of your domain to be
compressible via the fn_one_on_lambda Stokes constructor
parameter (refer to associated API reference for concise details).
If you encounter solver difficulties while using inflow/outflow, it is often the case that you are violating volume conservation, and are therefore trying to find a solution that cannot exist.
The Neumann Condition¶
The general process for specifying a Neumann condition is almost
identical to that for the Dirichlet condition. As before, you will
specify the required boundary using IndexSet objects. However, as
you are no longer fixing a value for the unknown itself (temperature,
velocity, etc), but instead a derived quantity (heat flux, stress) you
are now required to provide a Function object to specify the
necessary \(h\) function. Remember that MeshVariable objects are
also Function objects, so you are free to create another
MeshVariable to define \(h\) if that makes most sense. Refer to
the Function section of this guide for further details.
import underworld as uw
import glucifer
# 1. Create a mesh.
mesh = uw.mesh.FeMesh_Cartesian(elementRes=(32,32))
# 2. Create the required field(s).
tempField = mesh.add_variable(1)
# 3.a Create the Dirichlet Condition
botSet = mesh.specialSets['Bottom_VertexSet']
topSet = mesh.specialSets[ 'Top_VertexSet']
totSet = botSet + topSet
cond_dirichlet = uw.conditions.DirichletCondition(tempField,totSet)
tempField.data[botSet] = 373.
tempField.data[topSet] = 273.
# 3.b Create the Neumann condition
sidewalls = mesh.specialSets['Left_VertexSet'] + mesh.specialSets['Right_VertexSet']
# Create this constant value function explicitly so we can
# modify it later.
flux = uw.function.Function.convert(0.)
cond_neumann = uw.conditions.NeumannCondition(tempField,sidewalls,flux)
# 4. Diffusivity
alpha = 1.
# 5. Create the required system.
heatSys = uw.systems.SteadyStateHeat(temperatureField=tempField,
fn_diffusivity=alpha,
conditions=[cond_dirichlet,cond_neumann])
# 6. Create solver and solve.
solver = uw.systems.Solver(heatSys)
solver.solve()
# grab copy
nofluxTemp = tempField.copy(deepcopy=True)
# now set flux and re-solve
flux.value = -100.
solver.solve()
# vis.
fig_noflux = glucifer.Figure(figsize=(400,400), title="No Lateral Flux")
fig_noflux.append(glucifer.objects.Surface(mesh, nofluxTemp, colours="cubelaw2"))
fig_flux = glucifer.Figure(figsize=(400,400), title="Lateral Flux")
fig_flux.append(glucifer.objects.Surface(mesh, tempField, colours="cubelaw2"))
glucifer.Figure.show_grid([fig_noflux,fig_flux])
 |  |
Initial Conditions¶
For instantaneous systems (SteadyStateHeat,
SteadyStateDarcyFlow, Stokes) you will not generally be required
to set initial conditions, although care must be taken for non-linear
systems.
For the AdvectionDiffusion system you will be required to configure
your initial variable field as necessary. See the Mesh section of the
user guide for examples on initialising MeshVariable objects.
Naturally if your model utilises particles swarms, they will also need
to be configured (see the Swarm section of the user guide).
Defining Physical Quantities¶
You will create an Underworld Function object to define your
physical quantities. Refer to the Function section of the user guide for
details.
Non-Linearity¶
Underworld currently provides a basic Picard based iterative scheme for solving non-linear problems. Note, depending on the form of your non-linearity, a reasonable initial field configuration may be required to avoid solver convergence difficulty or division by zero type issues.
If your model has a non-linearity, the solver will flag you to this at
solve time, and you must explicitly set the solve() parameter
nonLinearIterate to True (Underworld will iterate for a
solution) or False (Underworld will not iterate, but will perform a
single solve). The False option might be used if you would like to
explicitly handle the non-linearity, or perform some analysis or
modifications after each iteration. An alternative option to enable
per-iteration interaction is offered by the creating a callback function
(passed in via the callback_post_solve paramater to the solve()
method).
import underworld as uw
import glucifer
# 1. Create a mesh.
mesh = uw.mesh.FeMesh_Cartesian(elementRes=(32,32))
# 2. Create the required field(s).
tempField = mesh.add_variable(1)
# 3.a Create the Dirichlet Condition
botSet = mesh.specialSets['Bottom_VertexSet']
topSet = mesh.specialSets[ 'Top_VertexSet']
totSet = botSet + topSet
cond_dirichlet = uw.conditions.DirichletCondition(tempField,totSet)
tempField.data[botSet] = 1.
tempField.data[topSet] = 0.
# 4. Diffusivity
alpha = uw.function.math.exp(tempField)
# 5. Create the required system.
heatSys = uw.systems.SteadyStateHeat(temperatureField=tempField,
fn_diffusivity=alpha,
conditions=[cond_dirichlet])
# 6. Create solver and solve.
solver = uw.systems.Solver(heatSys)
# THIS WILL NOT WORK!
# solver.solve()
# THIS WILL
solver.solve(nonLinearIterate=True)
# Use a callback
rms_temp = uw.utils.Integral(tempField*tempField,mesh)
import numpy as np
def mycallback():
print("Average temp is {}".format(np.sqrt(rms_temp.evaluate())))
# reset temp field first for demo
tempField.data[:] = 0.
tempField.data[botSet] = 1.
solver.solve(nonLinearIterate=True, callback_post_solve=mycallback)
Average temp is [ 0.58463673]
Average temp is [ 0.64936446]
Average temp is [ 0.6476043]
Average temp is [ 0.6476043]
Utilities¶
Overview:
- Integrals.
- Checkpointing.
- Generating XDMF files.
Keywords: checkpointing, utilities, volume integrals, surface integrals, xdmf
Integral Class¶
The Integral class constructs the volume integral
for some function \(f_i\) (specified by a Function object), over
some domain \(V\) (specified by an FeMesh object), or the
surface integral
for some surface \(\Gamma\) (specified via an IndexSet object on
the mesh).
# setup some objects for demonstration import underworld as uw from underworld import function as fn import glucifer import math # setup required objects mesh = uw.mesh.FeMesh_Cartesian(minCoord=(-1,-1), elementRes=(32,32)) temperatureField = mesh.add_variable( 1 ) velocityField = mesh.add_variable( 2 ) # init temperatureField.data[:,0] = -mesh.data[:,1]/2. + 0.5 velocityField.data[:,0] = -mesh.data[:,1] velocityField.data[:,1] = mesh.data[:,0] # viz fig1 = glucifer.Figure() velmagfield = uw.function.math.sqrt( uw.function.math.dot( velocityField, velocityField ) ) fig1.append( glucifer.objects.VectorArrows(mesh, velocityField, arrowHead=0.2, scaling=0.1) ) fig1.append( glucifer.objects.Surface( mesh, temperatureField ) ) fig1.show()
In the following example, we will calculate a root mean square velocity defined as:
The same result can be achieved through a number of paths.
- Call the
integrate()method on aFunctionobject. - Call the
integrate()method on anFeMeshobject. - Create an
Integralclass object, and call itsintegrate()method.
Note that all three methods result in identical calculations.
# note that the previous cell needs to have been executed before this one.
# define required function
vdotv = fn.math.dot( velocityField, velocityField )
# evaluate area (domain) integrals on the function
v2sum = mesh.integrate( vdotv ) # ... option 2
volume = mesh.integrate( 1.0 )
# finally, calculate RMS
v_rms = math.sqrt( v2sum[0] )/volume[0]
print('RMS velocity = {0:.3f}'.format(v_rms))
# option 3
# create integral objects, passing in functions
v2sum_integral = uw.utils.Integral( mesh=mesh, fn=vdotv )
volume_integral = uw.utils.Integral( mesh=mesh, fn=1. )
# evaluate integrals
v2sum = v2sum_integral.evaluate()
volume = volume_integral.evaluate()
# finally, calculate RMS
v_rms = math.sqrt( v2sum[0] )/volume[0]
print('RMS velocity = {0:.3f}'.format(v_rms))
RMS velocity = 0.408
RMS velocity = 0.408
To evaluate an integral over a subdomain the
fn.branching.conditional class may be useful:
import underworld as uw
from underworld import function as fn
import glucifer
mesh = uw.mesh.FeMesh_Cartesian(minCoord=(-1,-1), elementRes=(32,32))
# create circle function
radius = 1.
coord = fn.coord()
fn_sphere = fn.math.dot( coord, coord ) < radius**2
# setup a function that is 1 if the coordinates, are inside the circle, and zero otherwise.
conditions = [ ( fn_sphere , 1.0),
( True , 0.0) ]
kernelFunction = fn.branching.conditional( conditions )
# create and evaluate integral
volume = mesh.integrate( kernelFunction )
print('Area from integral = {0:6.8e}'.format(volume[0]))
fig = glucifer.Figure()
fig.append( glucifer.objects.Mesh(mesh) )
fig.append( glucifer.objects.Surface(mesh, kernelFunction))
fig.show()
Area from integral = 3.14559221e+00
In the following example, we will calculate the Nusselt number defined by
where \(h\) is the height of the domain, and \(\Gamma_t\) and \(\Gamma_b\) are the top and bottom surfaces respectively.
NB. Surface integrals must still be implemented using the old
implementation mesh of creating a uw.utils.Integral object
import underworld as uw
from underworld import function as fn
import glucifer
import math
mesh = uw.mesh.FeMesh_Cartesian(minCoord=(-1,-1), elementRes=(32,32))
temperatureField = mesh.add_variable( 1 )
temperatureField.data[:,0] = -mesh.data[:,1]/2. + 0.5
nuTop = uw.utils.Integral( fn=temperatureField.fn_gradient[1],
mesh=mesh, integrationType='Surface',
surfaceIndexSet=mesh.specialSets["MaxJ_VertexSet"])
nuBottom = uw.utils.Integral( fn=temperatureField,
mesh=mesh, integrationType='Surface',
surfaceIndexSet=mesh.specialSets["MinJ_VertexSet"])
nu = - nuTop.evaluate()[0]/nuBottom.evaluate()[0]
print('Nusselt number = {0:.6f}'.format(nu))
Nusselt number = 0.500000
Checkpointing¶
Checkpointing is the process of saving sufficient data to facilitate restarting your simulations at a later stage. Note that we do not provide explicit checkpointing functionality, but instead provide the tools required for the loading and saving of heavy data. Which data items are required for restart will depend on the systems you have used and how you have constructed your models.
The following Underworld data structures have load/save functionality:
* SwarmVariables * Swarm * MeshVariables * Mesh
All files are saved in HDF5 format.
Note: When saving a SwarmVariable, if you wish to reload the
SwarmVariable data at a later stage, you must also save the
Swarm for the corresponding state (generally, the same timestep).
This requirement is due to the population control mechanisms swarms
generally used for swarms, and also due to particles crossing process
boundaries. When you come to reload the Swarm and SwarmVariable,
you must load the Swarm first. Note again that the Swarm and
SwarmVariable must be of corresponding state for successful reload.
Swarms and Swarm Variables¶
Below a Swarm and a SwarmVariable, are created and saved to
disk, then a new swarm loads the data from disk.
- Save the swarm data to disk using the
save()method on theSwarmandSwarmVariableobjects. Note the handle object that is returned from thesave()method. This is currently used forxdmf()operation, see below. - Load the swarm data from disk using the
load()method on theSwarmandSwarmVariableobjects
import underworld as uw
from underworld import function as fn
import glucifer
import numpy as np
mesh = uw.mesh.FeMesh_Cartesian(minCoord=(-1,-1), elementRes=(32,32))
swarm1 = uw.swarm.Swarm(mesh)
swarm1var = swarm1.add_variable(dataType='int', count=1)
layout = uw.swarm.layouts.PerCellSpaceFillerLayout(swarm1,particlesPerCell=5)
swarm1.populate_using_layout(layout)
# set random data to save
swarm1var.data[:,0] = np.random.rand(len(swarm1var.data))
ignore = swarm1.save('swarm.h5')
ignore = swarm1var.save('swarmvar.h5')
# new swarm
swarm2 = uw.swarm.Swarm(mesh)
swarm2var = swarm2.add_variable(dataType='int', count=1)
swarm2.load('swarm.h5')
swarm2var.load('swarmvar.h5')
if not np.allclose( swarm2var.data[:], swarm1var.data[:] ):
raise RuntimeError("Something went wrong, the swarms are not identical.")
if not np.allclose( swarm2.particleCoordinates.data[:], swarm1.particleCoordinates.data[:]):
raise RuntimeError("Something went wrong, the swarms variables are not identical.")
if uw.rank() == 0:
import os
os.remove('swarm.h5')
os.remove('swarmvar.h5')
Mesh and Mesh Variables¶
The MeshVariable object behaves similarly to SwarmVariable
object with the save() and load() functionality.
import underworld as uw
from underworld import function as fn
import glucifer
import numpy as np
mesh = uw.mesh.FeMesh_Cartesian(minCoord=(-1,-1), elementRes=(32,32))
# create var and set random data
meshvar = mesh.add_variable(1)
meshvar.data[:,0] = np.random.rand(len(meshvar.data))
ignore = mesh.save('mesh.h5')
ignore = meshvar.save('meshvar.h5')
newvar = mesh.add_variable(1)
newvar.load('meshvar.h5')
# ensure expected result
if not np.allclose( newvar.data, meshvar.data ):
raise RuntimeError("Something went wrong, the mesh variables are not identical.")
if uw.rank() == 0:
import os
os.remove('mesh.h5')
os.remove('meshvar.h5')
Writing to XDMF files¶
The XDMF file specifies how your model data (MeshVariable
and/or SwarmVariable objects) and geometry (Mesh and/or
Swarm objects) is layed out within the HDF5 files you have saved to
disk. It is a light (or sidecar) file written in XML which accompanies
your heavy HDF5 files to allow applications such as
ParaView to visualise your models.
In Underworld, whenever you call a save() operation, we return
handler objects which specify the information required to create the
XDMF files. If you do not wish to create XDMF files, you can ignore
these handler objects. Otherwise, you need to retain them, and then pass
them into the xdmf() methods so that it has the information it
requires. Note that any object that has a save() method, will also
have an xdmf() method.
We generate XDMF files in the following example. We will also print an XDMF file, but note that you generally should not need to manually view or modify these files. For details on writing XDMF files in a dynamical context, please see the Rayleigh-Taylor example.
import underworld as uw
from underworld import function as fn
import glucifer
import numpy as np
# mesh
mesh = uw.mesh.FeMesh_Cartesian(elementType='Q1', elementRes=(4,4)) # 16 elements total
meshvar = mesh.add_variable(1) # 25 meshvar nodes total
xdmf_info_mesh = mesh.save('mesh.h5')
xdmf_info_meshvar = meshvar.save('meshvar.h5')
# swarm
swarm1 = uw.swarm.Swarm(mesh)
swarm1var = swarm1.add_variable(dataType='int', count=1)
layout = uw.swarm.layouts.PerCellSpaceFillerLayout(swarm1,particlesPerCell=5) # 80 particles total
swarm1.populate_using_layout(layout)
xdmf_info_swarm = swarm1.save('swarm.h5')
xdmf_info_swarmvar = swarm1var.save('swarmvar.h5')
# xdmf
meshvar.xdmf('meshvar.xdmf', xdmf_info_meshvar, "MyField", xdmf_info_mesh, "TheMesh", modeltime=0.0)
swarm1var.xdmf('swarmvar.xdmf', xdmf_info_swarmvar, "SwarmVariable", xdmf_info_swarm, "TheSwarm", modeltime=0.1)
# print an xdmf file
# printfile = 'meshvar.xdmf'
printfile = 'swarmvar.xdmf'
with open(printfile, 'r') as f:
print(f.read())
if uw.rank() == 0:
import os
os.remove('mesh.h5')
os.remove('meshvar.h5')
os.remove('swarm.h5')
os.remove('swarmvar.h5')
os.remove('meshvar.xdmf')
os.remove('swarmvar.xdmf')
<?xml version="1.0" ?>
<Xdmf xmlns:xi="http://www.w3.org/2001/XInclude" Version="2.0">
<Domain>
<Grid Name="TheSwarm" GridType="Uniform">
<Time Value="0.1" />
<Topology Type="POLYVERTEX" NodesPerElement="80"> </Topology>
<Geometry Type="XY">
<DataItem Format="HDF" NumberType="Float" Precision="8" Dimensions="80 2">swarm.h5:/data</DataItem>
</Geometry>
<Attribute Type="Scalar" Center="Node" Name="SwarmVariable">
<DataItem Format="HDF" NumberType="Float" Precision="8" Dimensions="80 1">swarmvar.h5:/data</DataItem>
</Attribute>
</Grid>
</Domain>
</Xdmf>
Visualisation¶
The glucifer module provides visualisation capabilities for Underworld modelling. It provides a higher level interface to the rendering capabilities provided by LavaVu, but also performs all the required collation of parallel data back to the root process, where it is then rendered.
- Creating figures.
- Drawing objects within figures.
- Saving figures to a file.
- Advanced figure control.
- The Interactive viewer
The Figure¶
The Figure class is the base container object for your glucifer
visualisations. It provides the canvas to which you will add the
renderings from your drawing objects. Use the show() method to
render your blank canvas:
import underworld as uw import glucifer newfigure = glucifer.Figure(facecolour="lightgrey", title="Empty",axis=True, figsize=(800,400)) newfigure.show()
Drawing Objects¶
Drawing objects are the items that are rendered within a Figure.
Available drawing objects (and associated documentation) may be found
within the glucifer.objects submodule or at the glucifer API
documentation page. We consider a few of the drawing objects below.
objects.Mesh¶
Render mesh geometry and node indices. Note, the append() method is
used to attach this drawing object to a figure object lists of drawing
objects.
import underworld as uw import glucifer fig = glucifer.Figure() # create mesh and display it mesh = uw.mesh.FeMesh_Cartesian( elementRes=(2,2), minCoord=(0.,0.), maxCoord=(2.,1.) ) fig = glucifer.Figure(figsize=(800,400)) fig.append( glucifer.objects.Mesh( mesh, nodeNumbers=True ) ) fig.show()
objects.Surface¶
This object will render a provided field across the faces of a mesh object.
import underworld as uw import glucifer fig = glucifer.Figure() mesh = uw.mesh.FeMesh_Cartesian( elementRes=(2,2), minCoord=(0.,0.), maxCoord=(2.,1.) ) # create an object with a single value at each mesh point fevar = mesh.add_variable( 1 ) # give the variable some values fevar.data[:] = 0. fevar.data[0] = 10. fevar.data[4] = 30. fevar.data[8] = 10. fig = glucifer.Figure ( edgecolour="black", figsize=(800,400) ) fig.append( glucifer.objects.Surface( mesh, fevar, colours="red yellow green", onMesh=False ) ) fig.show()
objects.VectorArrows¶
This object will draw an array of vector arrows across the image using
the provided vector field to determine their direction. Check
help(VectorArrows) for the full options of the VectorArrows drawing
object.
import underworld as uw import glucifer fig = glucifer.Figure() mesh = uw.mesh.FeMesh_Cartesian( elementRes=(2,2), minCoord=(0.,0.), maxCoord=(2.,1.) ) # create an object with a single value at each mesh point fevar = mesh.add_variable(2) # initialise the vector field to something fevar.data[:] = [ 0., 0.] fevar.data[0] = [ 1., 1.] fevar.data[1] = [ 0., 1.] fevar.data[2] = [-1., 1.] fevar.data[3] = [ 1., 0.] fevar.data[5] = [-1., 0.] fevar.data[6] = [ 1.,-1.] fevar.data[7] = [ 0.,-1.] fevar.data[8] = [-1.,-1.] fig = glucifer.Figure( figsize=(800,400), edgecolour="black" ) fig.append( glucifer.objects.VectorArrows( mesh, fevar, scaling=0.1, arrowHead=0.2 ) ) fig.show()
glucifer.Points¶
This object will draw a swarm of points using the provided Underworld swarm for the point locations and an Underworld swarm variable for the point colours (or size or opacity).
import underworld as uw import glucifer fig = glucifer.Figure() mesh = uw.mesh.FeMesh_Cartesian( elementRes=(2,2), minCoord=(0.,0.), maxCoord=(2.,1.) ) swarm = uw.swarm.Swarm( mesh=mesh ) layout = uw.swarm.layouts.PerCellSpaceFillerLayout( swarm=swarm, particlesPerCell=200 ) swarm.populate_using_layout( layout ) # render points, and also colour them by depth fig = glucifer.Figure( figsize=(800,400), edgecolour="blue" ) fig.append(glucifer.objects.Points( swarm=swarm, fn_colour=uw.function.input()[1], pointSize=5, )) fig.show()
Multiple drawing objects can, of course, be layered upon each other to
build up complex images. This is achieved by making multiple calls
append():
import underworld as uw import glucifer fig = glucifer.Figure() mesh = uw.mesh.FeMesh_Cartesian( elementRes=(2,2), minCoord=(0.,0.), maxCoord=(2.,1.) ) # create vec field, setting it to coordinate vector vec = mesh.add_variable(2) vec.data[:] = mesh.data[:] fig = glucifer.Figure( figsize=(800,400), edgecolour="black" ) fig.append(glucifer.objects.VectorArrows( mesh, vec, scaling=0.1, arrowHead=0.2, opacity=0.6) ) fig.append(glucifer.objects.Surface( mesh, uw.function.math.dot(vec,vec), colours="red yellow green") ) fig.append(glucifer.objects.Mesh(mesh)) fig.show()
Saving Results¶
To output results to raster files (such as PNG), simply use the
savefig() method.
import underworld as uw
import glucifer
fig = glucifer.Figure()
figfile = fig.save_image("savedfigure")
import glob
import os
# do we have a file?
if figfile:
print(glob.glob( figfile ))
# cleanup
os.remove( figfile )
['savedfigure.png']
Advanced Figure Control¶
Various imaging parameters are provided directly through the
glucifer API, and the user is encouraged to explore these options
interactively and also through the API documentation. Note that for full
control of all image parameters, the user should consult the LavaVu
Documentation. In
particular, the LavaVu Property
Reference
is useful, as is the colour maps example
notebook.
We have a quick look at some available options below.
import underworld as uw
import glucifer
import numpy as np
# Create figure, set the margin around the edges of the plot and apply rulers with tick labels.
fig = glucifer.Figure(figsize=(1000,500),
title="Test Figure",
quality=2,
margin=0.06,
rulers=True,
rulerticks=5,
rulerwidth=0.5,
fontscale=0.7)
# create a var and init with some data
mesh = uw.mesh.FeMesh_Cartesian(elementRes=(64,32), minCoord=(-10.,-5), maxCoord=(10.,5.))
meshvar = mesh.add_variable(1)
meshvar.data[:,0] = mesh.data[:,0]**2 + mesh.data[:,1]**2
# create surf object.
# note that colours specification. for further details and options,
# consult the LavaVu colour maps example notebook.
surf = glucifer.objects.Surface( mesh, meshvar, colours=' red (20)white \
(21)black (40)black \
(41)white (60)white \
(61)black (80)black \
(81)white (100)white \
(101)black (120)black \
(121)white (140)white ',
onMesh=False, resolution=500, range=[0.,125])
# go ahead and add object to figure
fig.append(surf)
# we can further modify the colour bar associated with a drawing object by
# acessing it directly.
# set the number of extra tick marks (besides start/finish)
surf.colourBar["ticks"] = 6
# # set values for extra ticks
surf.colourBar["tickvalues"] = [20, 40, 60, 80, 100, 120]
# set size
surf.colourBar["size"] = [0.75,0.02]
# set alignment. Can be left/right (vertical colour bar), bottom/top (horizontal colour bar)
surf.colourBar["align"] = "right"
# set position
surf.colourBar["position"] = 0.
# also let's add an annotation
fig.draw.label("Label In Centre", (0,0,0), scaling=10.)
fig.show()
Interactivity & Scripting¶
An interactive viewer can be opened from inside the notebook by importing the “lavavu” module and using the figure’s window() method. Once this interactive viewer is open we can either shift the view using the mouse, or we can send commands to shift the view from inside the notebook.
Note that interactive usage is only possible from within a notebook, and will not work from the ReadTheDocs html widgets. You can access the notebook document from which this page is derived from the Underworld github page.
import underworld as uw
import glucifer
mesh = uw.mesh.FeMesh_Cartesian()
var = mesh.add_variable(2)
var.data[:] = mesh.data[:]
fig = glucifer.Figure( edgecolour="black", quality=3 )
fig.Surface( mesh, uw.function.math.dot(var,var), colours="red yellow green", colourBar = False )
fig.VectorArrows( mesh, var, scaling=0.1, arrowHead=0.2)
if not glucifer.lavavu:
raise RuntimeError("LavaVu not available.") #Stop notebook here if no vis enabled
lv = fig.window()
lv.rotate('y', 45)
lv.redisplay()
We can retrieve a list of objects from a viewer by inspecting the “objects” property. References to these objects can be used to modify the object appearance, either directly or using controls.
lv.reset() #Restore camera to default view print lv.objects surf = lv.objects["ScalarField_0"] surf["opacity"] = 0.75 lv.redisplay()
Interactive controls can be created to adjust visualisation parameters, here we control the opacity of the previously retrieved “surf” object and the global vector scaling parameter, changes will be reflected in the viewer window above:
surf.control.Range(property="opacity", range=(0,1), command="reload")
lv.control.Range("scalevectors", value=1., range=(0.1, 2.5), step=0.1, label="Scale Arrows")
lv.control.show()
Stokes Solver¶
We want to solve the following Stokes block system.
If we apply Gaussian elimination to the above as a 2x2 block matrix system we can write this as:
where \(S=G^{T}K^{-1}G-C\) is the Schur complement and \(\hat{h}=G^{T}K^{-1}f -h\).
This system is now solved first for the pressure using the Schur complement matrix, \(S\). Then a backsolve for the velocity gives the complete solution.
Note that wherever \(K^{-1}\) appears, the inverse is never explicitly calculated but is achieved via a PETSc solve method. While solving for the pressure, there are necessarily solves using \(K\) inside of the matrix \(S\). We often refer to these as ‘inner’ solves.
Basic usage of the Stokes solver class involves being able to easily set up the inner solves in a few different ways (Setting up the pressure solve is more advanced).
To illustrate some basic usage let’s set up a simple problem to solve.
import underworld as uw
from underworld import function as fn
res=128
mesh = uw.mesh.FeMesh_Cartesian(elementRes=(res,res))
velocityField = mesh.add_variable(2)
pressureField = mesh.subMesh.add_variable(1)
velocityField.data[:] = (0.,0.)
pressureField.data[:] = 0.
# We are going to make use of one of the existing analytic solutions so that we may easily
# obtain functions for a viscosity profile and forcing terms.
# Exact solution solCx with defaults
sol = fn.analytic.SolCx()
stokesSystem = uw.systems.Stokes(velocityField,pressureField,sol.fn_viscosity,sol.fn_bodyforce,conditions=sol.get_bcs(velocityField))
# Now we create a solver.
solver=uw.systems.Solver(stokesSystem)
# The Stokes solver will use multigrid as a preconditioner along with
# PETSc's 'fgmres' Krylov method by default for the 'inner' solve.
solver.solve()
# Now let's set up the inner solve to do a direct solve.
# Note that the `lu` direct solver will not work in parallel.
solver.set_inner_method("lu")
solver.solve()
Let’s run underworld’s help function on the solver.configure function.
help(solver.set_inner_method)
Help on method set_inner_method in module underworld.systems._bsscr:
set_inner_method(self, solve_type='mg') method of underworld.systems._bsscr.StokesSolver instance
Configure velocity/inner solver (A11 PETSc prefix).
Available options below. Note that associated solver software
(for example mumps) must be installed.
- mg : Geometric multigrid (default).
- nomg : Disables multigrid.
- lu : LU direct solver (serial only).
- mumps : MUMPS parallel direct solver.
- superludist : SuperLU parallel direct solver.
- superlu : SuperLU direct solver (serial only).
We can see all the of the options that are configured in the solver
using the list() functions for each component of the solver. These
are the most useful ones.
print("System Level Options:")
solver.options.main.list()
print("")
print("Schur Complement Solve Options:")
solver.options.scr.list() # Specifics for the
print("")
print("Inner (velocity) Solve Options:")
solver.options.A11.list() # Specifics for the inner (velocity) solve
print("")
print("Multigrid (where enabled) Options:")
solver.options.mg.list() # Options relevant to multigrid (if chosen)
System Level Options:
('remove_constant_pressure_null_space', False)
('ksp_k2_type', 'NULL')
('change_backsolve', False)
('penalty', 0.0)
('pc_type', 'none')
('force_correction', True)
('k_scale_only', True)
('Q22_pc_type', 'uw')
('change_A11rhspresolve', False)
('ksp_type', 'bsscr')
('rescale_equations', False)
('restore_K', False)
Schur Complement Solve Options:
('ksp_type', 'fgmres')
('ksp_rtol', 1e-05)
Inner (velocity) Solve Options:
('pc_type', 'lu')
('_mg_active', False)
('ksp_type', 'preonly')
Multigrid (where enabled) Options:
('active', False)
('levels', 8)
Further information about options is available via the help() Python
function:
help(solver.options.A11)
The options generally follow PETSc naming conventions.
A useful trick is to be able to imitate the classic “penalty method” which can be very efficient with modest-sized (2D) problems.
In the penalty method, we add a term to the weak form of the Stokes equation which penalises \(\lambda | \nabla \cdot \mathbf{u}| > 0\) and where \(\lambda\) is a sufficiently large constant to enforce the constraint. Typically \(10^7\) is considered sufficient.
The problem with this method is that the condition number of the system is severely compromised by adding the penalty term and standard iterative methods do not work well.
Our solvers have been configured with the penalty term and, for sufficiently robust choices of the inner solver, this can help solve problems faster (by improving pressure convergence).
An indestructible solver like lu or mumps (Mumps is a direct
solve that will work in parallel) can use very large penalties. Hence we
can recreate the penalty method as follows (though it still follows the
pattern of the Schur complement solver, while the classical method takes
some shortcuts).
solver.set_inner_method("mumps")
solver.options.scr.ksp_type="cg"
solver.set_penalty(1.0e7)
solver.solve()
solver.print_stats()
[1;35m
Pressure iterations: 4
Velocity iterations: 1 (presolve)
Velocity iterations: 4 (pressure solve)
Velocity iterations: 1 (backsolve)
Velocity iterations: 6 (total solve)
SCR RHS setup time: 4.8789e-01
SCR RHS solve time: 5.6970e-03
Pressure setup time: 1.0583e-03
Pressure solve time: 2.6979e-02
Velocity setup time: 8.5831e-06 (backsolve)
Velocity solve time: 5.4877e-03 (backsolve)
Total solve time : 5.7443e-01
Velocity solution min/max: 0.0000e+00/0.0000e+00
Pressure solution min/max: 0.0000e+00/0.0000e+00
[00m
Now let’s go back to using multigrid. We can use a penalty here too, but the gigantic numbers won’t work.
solver.set_inner_method("mg")
solver.set_penalty(1.0)
solver.solve()
solver.print_stats()
[1;35m
Pressure iterations: 17
Velocity iterations: 6 (presolve)
Velocity iterations: 88 (pressure solve)
Velocity iterations: 5 (backsolve)
Velocity iterations: 99 (total solve)
SCR RHS setup time: 2.8383e-02
SCR RHS solve time: 8.8551e-02
Pressure setup time: 1.5616e-03
Pressure solve time: 8.6988e-01
Velocity setup time: 9.5367e-06 (backsolve)
Velocity solve time: 3.8740e-02 (backsolve)
Total solve time : 1.1070e+00
Velocity solution min/max: 0.0000e+00/0.0000e+00
Pressure solution min/max: 0.0000e+00/0.0000e+00
[00m
underworld module¶
Underworld2 is a python-friendly version of the Underworld geodynamics code which provides a programmable and flexible front end to all the functionality of the code running in a parallel HPC environment. This gives signficant advantages to the user, with access to the power of python libraries for setup of complex problems, analysis at runtime, problem steering, and coupling of multiple problems.
Underworld2 is integrated with the literate programming environment of the jupyter notebook system for tutorials and as a teaching tool for solid Earth geoscience.
Underworld is an open-source, particle-in-cell finite element code tuned for large-scale geodynamics simulations. The numerical algorithms allow the tracking of history information through the high-strain deformation associated with fluid flow (for example, transport of the stress tensor in a viscoelastic, convecting medium, or the advection of fine-scale damage parameters by the large-scale flow). The finite element mesh can be static or dynamic, but it is not contrained to move in lock-step with the evolving geometry of the fluid. This hybrid approach is very well suited to complex fluids which is how the solid Earth behaves on a geological timescale.
Module Summary¶
submodules:¶
underworld.function module¶
The function module contains the Function class, and related classes.
Function objects are constructed in python, but evaluated in C for efficiency. They provide a high level interface for users to compose model behaviour (such as viscosity), as well as a natural interface by which discrete data (such as meshvariables) may be utilised.
The branching module provides functions which provide branching behaviour. Typically, these functions will select other user provided functions when certain conditions are met (with the condition also described by a function!).
underworld.function.branching.map |
This function performs a map to other functions. |
underworld.function.branching.conditional |
This function provides ‘if/elif’ type conditional behaviour. |
-
class
underworld.function.branching.map(fn_key=None, mapping=None, fn_default=None, *args, **kwargs)[source]¶ Bases:
underworld.function._function.FunctionThis function performs a map to other functions. The user provides a python dictionary which maps unsigned integers keys to underworld functions. The user must also provide a key function. At evaluation time, the key function is evaluated first, with the outcome determining which function should finally be evaluated to return a value.
For a set of value functions \(\{f_{v_0},f_{v_1},\ldots,f_{v_n}\}\), corresponding keys \(\{k_0,k_1,\ldots,k_n\}\), and key function \(f_{k}\), we have:
\[\begin{split}f(\mathbf{r})= \begin{cases} f_{v_0}(\mathbf{r}), & \text{if } f_{k}(\mathbf{r}) = k_0\\ f_{v_1}(\mathbf{r}), & \text{if } f_{k}(\mathbf{r}) = k_1\\ ... \\ f_{v_n}(\mathbf{r}), & \text{if } f_{k}(\mathbf{r}) = k_n\\ f_{d} (\mathbf{r}), & \text{otherwise} \end{cases}\end{split}\]As stated, the keys must be unsigned integers. The key function need not return an unsigned integer, but whatever value it returns will be cast to an unsigned integer so caution is advised.
The default function is optional, but if none is provided, and the key function evaluates to a value which is not within the user provide set of keys, an exception will be thrown.
Parameters: - fn_key (underworld.function.Function (or convertible)) – Function which returns integer key values. This function will be evaluated first to determine which function from the mapping is to be used.
- mapping (dict(Function)) – Python dictionary providing a mapping from unsigned integer ‘key’ values to underworld ‘value’ functions. Note that the provided ‘value’ functions must return values of type ‘double’.
- fn_default (underworld.function.Function (or convertible) (optional)) – Default function to be utilised when the key (returned by fn_key function) does not correspond to any key value in the mapping dictionary.
The following example sets different function behaviour inside and outside of a unit sphere. The unit sphere is represented by particles which record a swarm variable to determine if they are or not inside the sphere.
Example
Setup mesh, swarm, swarmvariable & populate
>>> import underworld as uw >>> import underworld.function as fn >>> import numpy as np >>> mesh = uw.mesh.FeMesh_Cartesian(elementRes=(8,8),minCoord=(-1.0, -1.0), maxCoord=(1.0, 1.0)) >>> swarm = uw.swarm.Swarm(mesh) >>> svar = swarm.add_variable("int",1) >>> swarm.populate_using_layout(uw.swarm.layouts.PerCellSpaceFillerLayout(swarm,20))
For all particles in unit circle, set svar to 1
>>> svar.data[:] = 0 >>> for index, position in enumerate(swarm.particleCoordinates.data): ... if position[0]**2 + position[1]**2 < 1.: ... svar.data[index] = 1
Create a function which reports the value ‘1.’ inside the sphere, and ‘0.’ otherwise. Note that while we have only used constant value functions here, you can use any object of the class Function.
>>> fn_map = fn.branching.map(fn_key=svar, mapping={0: 0., 1:1.}) >>> np.allclose(np.pi, uw.utils.Integral(fn_map,mesh).evaluate(),rtol=2e-2) True
Alternatively, we could utilise the default function to achieve the same result.
>>> fn_map = fn.branching.map(fn_key=svar, mapping={1: 1.}, fn_default=0.) >>> np.allclose(np.pi, uw.utils.Integral(fn_map,mesh).evaluate(),rtol=2e-2) True
-
class
underworld.function.branching.conditional(clauses, *args, **kwargs)[source]¶ Bases:
underworld.function._function.FunctionThis function provides ‘if/elif’ type conditional behaviour.
The user provides a list of tuples, with each tuple being of the form (fn_condition, fn_resultant). Effectively, each tuple provides a clause within the if/elif statement.
When evaluated, the function traverses the clauses, stopping at the first fn_condition which returns ‘true’. It then executes the corresponding fn_resultant and returns the results.
If none of the provided clauses return a ‘True’ result, an exception is raised.
For a set of condition functions { fc_0, fc_1, … ,fc_n }, and corresponding resultant functions { fr_0, fr_1, … ,fr_n }, we have for a provided input f_in:
if fc_0(f_in) : return fr_0(f_in) elif fc_1(f_in) : return fr_1(f_in) ... elif fc_n(f_in) : return fr_n(f_in) else : raise RuntimeError("Reached end of conditional statement. At least one of the clause conditions must evaluate to 'True'." );
Parameters: clauses (list) – list of tuples, with each tuple being of the form (fn_condition, fn_resultant). Example
The following example uses functions to represent a unit circle. Here a conditional function report back the value ‘1.’ inside the sphere (as per the first condition), and ‘0.’ otherwise.
>>> import underworld as uw >>> import underworld.function as fn >>> import numpy as np >>> mesh = uw.mesh.FeMesh_Cartesian(elementRes=(16,16),minCoord=(-1.0, -1.0), maxCoord=(1.0, 1.0)) >>> circleFn = fn.coord()[0]**2 + fn.coord()[1]**2 >>> fn_conditional = fn.branching.conditional( [ (circleFn < 1., 1. ), ( True, 0. ) ] ) >>> np.allclose(np.pi, uw.utils.Integral(fn_conditional,mesh).evaluate(),rtol=1e-2) True
This module provides functions which raise an exception when given conditions are encountered during function evaluations. Exception functions never modify query data.
underworld.function.exception.CustomException |
This function allows you to set custom exceptions within your model. |
underworld.function.exception.SafeMaths |
This function checks if any of the following have been encountered during the evaluation of its argument function: |
-
class
underworld.function.exception.CustomException(fn_input, fn_condition, fn_print=None, *args, **kwargs)[source]¶ Bases:
underworld.function._function.FunctionThis function allows you to set custom exceptions within your model. You must pass it two functions: the first function is the pass through function, the second function is the required condition. You may also pass in a optional third function whose output will be printed if the condition evaluates to False.
A CustomException function will perform the following logic:
- Evaluate the condition function.
- If it evaluates to False, an exception is thrown and the simulation is halted. If a print function is provided, it will be evaluated and its results will be included in the exception message.
- If it evaluates to True, the pass through function is evaluated with the result then being return.
Parameters: - fn_passthrough (underworld.function.Function) – The pass through function
- fn_condition (underworld.function.Function) – The condition function
- fn_print (underworld.function.Function) – The print function (optional).
Example
>>> import underworld as uw >>> import underworld.function as fn >>> one = fn.misc.constant(1.) >>> passing_one = fn.exception.CustomException( one, (one < 2.) ) >>> passing_one.evaluate() array([[ 1.]]) >>> failing_one = fn.exception.CustomException( one, (one > 2.) ) >>> failing_one.evaluate() Traceback (most recent call last): ... RuntimeError: Issue utilising function of class 'CustomException' constructed at: --- CONSTRUCTION TIME STACK --- Error message: CustomException condition function has evaluated to False for current input!
Now with printing
>>> failing_one_by_five = fn.exception.CustomException( one, (one*5. > 20.), one*5. ) >>> failing_one_by_five.evaluate() Traceback (most recent call last): ... RuntimeError: Issue utilising function of class 'CustomException' constructed at: --- CONSTRUCTION TIME STACK --- Error message: CustomException condition function has evaluated to False for current input! Print function returns the following values (cast to double precision): ( 5 )
-
class
underworld.function.exception.SafeMaths(fn, *args, **kwargs)[source]¶ Bases:
underworld.function._function.FunctionThis function checks if any of the following have been encountered during the evaluation of its argument function:
- Divide by zero
- Invalid domain was used for evaluation
- Value overflow errors
- Value underflow errors
If any of the above are encountered, an exception is thrown at the conclusion of the argument function evaluation.
Parameters: fn (underworld.function.Function) – The function that is subject to the testing. Example
>>> import underworld as uw >>> import underworld.function as fn >>> one = fn.misc.constant(1.) >>> zero = fn.misc.constant(0.) >>> fn_dividebyzero = one/zero >>> safedividebyzero = fn.exception.SafeMaths(fn_dividebyzero) >>> safedividebyzero.evaluate() Traceback (most recent call last): ... RuntimeError: Issue utilising function of class 'SafeMaths' constructed at: --- CONSTRUCTION TIME STACK --- Error message: Floating point exception(s) encountered while evaluating SafeMaths argument function: Divide by zero
This module provides functions relating to tensor operations.
All Underworld2 functions return 1d array type objects. For tensor objects, the following convention is used:
Full tensors:
- 2D:
- \[\begin{split}\left[ a_{00}, a_{01}, \\ a_{10}, a_{11} \right]\end{split}\]
- 3D:
- \[\begin{split}\left[ a_{00}, a_{01}, a_{02}, \\ a_{10}, a_{11}, a_{12}, \\ a_{20}, a_{21}, a_{22} \right]\end{split}\]
Symmetric tensors:
- 2D:
- \[\left[ a_{00}, a_{11}, a_{01} \right]\]
- 3D:
- \[\left[ a_{00}, a_{11}, a_{22}, a_{01}, a_{02}, a_{12} \right]\]
underworld.function.tensor.symmetric |
This function calculates the symmetric part of a tensor and returns it as a symmetric tensor. |
underworld.function.tensor.deviatoric |
This function calculates the deviatoric stress tensor from the provided symmetric tensor. |
underworld.function.tensor.antisymmetric |
This function calculates the anti-symmetric part of a tensor, returning it as a tensor. |
underworld.function.tensor.second_invariant |
This function calculates the second invariant of (symmetric)tensor provided by the subject function. |
-
class
underworld.function.tensor.symmetric(fn, *args, **kwargs)[source]¶ Bases:
underworld.function._function.FunctionThis function calculates the symmetric part of a tensor and returns it as a symmetric tensor. The function generated by this class returns objects of type SymmetricTensorType.
\[v_{ij} = \tfrac{1}{2} ( u_{ij} + u_{ji} )\]Parameters: fn (underworld.function.Function) – The function which provides the required tensor. This function must return objects of type TensorType.
-
class
underworld.function.tensor.deviatoric(fn, *args, **kwargs)[source]¶ Bases:
underworld.function._function.FunctionThis function calculates the deviatoric stress tensor from the provided symmetric tensor. The function generated by this class returns objects of type SymmetricTensorType.
\[\tau_{ij} = \sigma_{ij} - \frac{\sigma_{kk}}{\delta_{ll}}\delta_{ij}\]Parameters: fn (underworld.function.Function) – The function which provides the required stress symmetric tensor. This function must return objects of type SymmetricTensorType.
-
class
underworld.function.tensor.antisymmetric(fn, *args, **kwargs)[source]¶ Bases:
underworld.function._function.FunctionThis function calculates the anti-symmetric part of a tensor, returning it as a tensor. The function generated by this class returns objects of type TensorType.
\[v_{ij} = \tfrac{1}{2} ( u_{ij} - u_{ji} )\]Parameters: fn (underworld.function.Function) – The function which provides the required tensor. This function must return objects of type TensorType.
-
class
underworld.function.tensor.second_invariant(fn, *args, **kwargs)[source]¶ Bases:
underworld.function._function.FunctionThis function calculates the second invariant of (symmetric)tensor provided by the subject function. The function generated by this class returns objects of type ScalarType.
\[u = \sqrt{ \tfrac{1}{2} u_{ij} u_{ij} }\]Parameters: fn (underworld.function.Function) – The function which provides the required tensor. This function must return objects of type TensorType or SymmetricTensorType.
Miscellaneous functions.
underworld.function.misc.max |
Returns the maximum of the results returned from its two argument function. |
underworld.function.misc.constant |
This function returns a constant value. |
underworld.function.misc.min |
Returns the minimum of the results returned from its two argument function. |
-
class
underworld.function.misc.max(fn1, fn2, **kwargs)[source]¶ Bases:
underworld.function._function.FunctionReturns the maximum of the results returned from its two argument function.
Parameters: - fn1 (underworld.function.Function) – First argument function. Function must return a float type.
- fn2 (underworld.function.Function) – Second argument function. Function must return a float type.
Example
>>> import underworld as uw >>> import underworld.function as fn >>> import numpy as np >>> testpoints = np.array(([[ 0.0], [0.2], [0.4], [0.6], [0.8], [1.01], [1.2], [1.4], [1.6], [1.8], [2.0],]))
Create which return identical results via different paths:
>>> fn_x = fn.input()[0] >>> fn_x_minus_one = fn_x - 1. >>> fn_one_minus_x = 1. - fn_x
Here we use ‘max’ and ‘min’ functions:
>>> fn_max = fn.misc.max(fn_one_minus_x,fn_x_minus_one) >>> fn_min = fn.misc.min(fn_one_minus_x,fn_x_minus_one)
Here we use the conditional functions:
>>> fn_conditional_max = fn.branching.conditional( ( ( fn_x <= 1., fn_one_minus_x ), ( fn_x > 1., fn_x_minus_one ) )) >>> fn_conditional_min = fn.branching.conditional( ( ( fn_x >= 1., fn_one_minus_x ), ( fn_x < 1., fn_x_minus_one ) ))
They should return identical results:
>>> np.allclose(fn_max.evaluate(testpoints),fn_conditional_max.evaluate(testpoints)) True >>> np.allclose(fn_min.evaluate(testpoints),fn_conditional_min.evaluate(testpoints)) True
-
class
underworld.function.misc.constant(value, *args, **kwargs)[source]¶ Bases:
underworld.function._function.FunctionThis function returns a constant value.
Parameters: value (int,float,bool, iterable) – The value the function should return. Note that iterable objects which contain valid types are permitted, but must be homogeneous in their type. Example
>>> import underworld as uw >>> import underworld.function as fn >>> fn_const = fn.misc.constant( 3 ) >>> fn_const.evaluate(0.) # eval anywhere for constant array([[3]], dtype=int32) >>> fn_const = fn.misc.constant( (3,2,1) ) >>> fn_const.evaluate(0.) # eval anywhere for constant array([[3, 2, 1]], dtype=int32) >>> fn_const = fn.misc.constant( 3. ) >>> fn_const.evaluate(0.) # eval anywhere for constant array([[ 3.]]) >>> fn_const = fn.misc.constant( (3.,2.,1.) ) >>> fn_const.evaluate(0.) # eval anywhere for constant array([[ 3., 2., 1.]]) >>> fn_const = fn.misc.constant( True ) >>> fn_const.evaluate(0.) # eval anywhere for constant array([[ True]], dtype=bool) >>> fn_const = fn.misc.constant( (True,False,True) ) >>> fn_const.evaluate(0.) # eval anywhere for constant array([[ True, False, True]], dtype=bool)
-
value¶ constant value this function returns
Type: value
-
-
class
underworld.function.misc.min(fn1, fn2, **kwargs)[source]¶ Bases:
underworld.function._function.FunctionReturns the minimum of the results returned from its two argument function.
Parameters: - fn1 (underworld.function.Function) – First argument function. Function must return a float type.
- fn2 (underworld.function.Function) – Second argument function. Function must return a float type.
Example
See the example provided for ‘max’ function.
This module provides a suite of models which satisfy the Stokes system of equations.
All models are considered across a unit square (or cube) domain, and utilise (unless otherwise stated) free-slip conditions on all boundaries.
Each model object provides a set of Underworld Functions for description of physical quantities such as velocity, pressure and viscosity.
For numerical validation in Underworld, we construct a Stokes system with appropriate domain and boundary conditions. Viscosity and body forces are set directly using corresponding Functions provided by the solution object. Generated numerical solution for velocity and pressure (or derivated quantities) may then be compared with exact solutions provided by solution objects.
underworld.function.analytic.SolDB3d |
SolDB2d and solDB3d from: |
underworld.function.analytic.SolM |
|
underworld.function.analytic.SolH |
Density step profile in (x,y). |
underworld.function.analytic.SolCx |
Viscosity step profile in x, trigonometric density profile. |
underworld.function.analytic.SolB |
Trigonometric/hyperbolic body forcing. |
underworld.function.analytic.SolC |
Discontinuous body forcing. |
underworld.function.analytic.SolDA |
Columnar density profile in x, and viscosity step in z. |
underworld.function.analytic.SolA |
Trigonometric body forcing. |
underworld.function.analytic.SolKz |
The boundary conditions are free-slip everywhere on a unit domain. |
underworld.function.analytic.SolDB2d |
SolDB2d and solDB3d from: |
underworld.function.analytic.SolKx |
The boundary conditions are free-slip everywhere on a unit domain. |
underworld.function.analytic.SolNL |
SolNL requires tighter solver tolerances and/or a direct solve for best results. |
-
class
underworld.function.analytic.SolDB3d(Beta=4.0, *args, **kwargs)[source]¶ Bases:
underworld.function.analytic._SolBaseFixedBcSolDB2d and solDB3d from:
Dohrmann, C.R., Bochev, P.B., A stabilized finite element method for the Stokes problem based on polynomial pressure projections, Int. J. Numer. Meth. Fluids 46, 183-201 (2004).
-
class
underworld.function.analytic.SolM(eta_0=1.0, n_x=1, n_z=1, r=1.5, *args, **kwargs)[source]¶ Bases:
underworld.function.analytic._SolBaseFreeSlipBc
-
class
underworld.function.analytic.SolH(sigma_0=1.0, x_c=0.5, y_c=0.5, eta_0=1.0, nmodes=30, *args, **kwargs)[source]¶ Bases:
underworld.function.analytic._SolBaseFreeSlipBcDensity step profile in (x,y). Constant viscosity.
Parameters: - sigma_0 (float) – Perturbation strength factor.
- x_c (float) – Step position (in x).
- y_c (float) – Step position (in y).
- eta_0 (float) – Viscosity.
- nmodes (int) – Number of Fourier modes used when evaluating analytic solution.
-
class
underworld.function.analytic.SolCx(n_x=1, eta_A=1.0, eta_B=100000.0, x_c=0.75, *args, **kwargs)[source]¶ Bases:
underworld.function.analytic._SolBaseFreeSlipBcViscosity step profile in x, trigonometric density profile.
Parameters: - n_x (unsigned) – Wavenumber parameter (in x).
- eta_A (float) – Viscosity of region A.
- eta_B (float) – Viscosity of region B.
- eta_c (float) – Viscosity step location.
-
class
underworld.function.analytic.SolB(sigma_0=1.0, n_x=1, n_z=1.5, eta_0=1.0, *args, **kwargs)[source]¶ Bases:
underworld.function.analytic._SolBaseFreeSlipBcTrigonometric/hyperbolic body forcing. Isoviscous.
Parameters: - sigma_0 (float) – Perturbation strength factor.
- n_x (int) – Wavenumber parameter (in x).
- n_z (float) – Wavenumber parameter (in z).
- eta_0 (float) – Viscosity.
-
class
underworld.function.analytic.SolC(sigma_0=1.0, x_c=0.5, eta_0=1.0, nmodes=200, *args, **kwargs)[source]¶ Bases:
underworld.function.analytic._SolBaseFreeSlipBcDiscontinuous body forcing. Isoviscous.
Parameters: - sigma_0 (float) – Perturbation strength factor.
- x_c (float) – Perturbation step location.
- eta_0 (float) – Viscosity.
- nmodes (int) – Number of Fourier modes used when evaluating analytic solution.
Notes
This solution is quiet slow to evaluate due to large number of Fourier terms required. Number of terms is hard code in solC.c.
-
class
underworld.function.analytic.SolDA(sigma_0=1.0, x_c=0.375, x_w=0.25, eta_A=1.0, eta_B=10.0, z_c=0.75, nmodes=200, *args, **kwargs)[source]¶ Bases:
underworld.function.analytic._SolBaseFreeSlipBcColumnar density profile in x, and viscosity step in z.
Parameters: - sigma_0 (float) – Perturbation strength factor.
- x_c (float) – Centre of column.
- x_w (float) – Width of column.
- eta_A (float) – Viscosity of region A.
- eta_B (float) – Viscosity of region B.
- z_c (float) – Viscosity step location.
- nmodes (int) – Number of Fourier modes used when evaluating analytic solution.
-
class
underworld.function.analytic.SolA(sigma_0=1.0, n_x=1, n_z=1.0, eta_0=1.0, *args, **kwargs)[source]¶ Bases:
underworld.function.analytic._SolBaseFreeSlipBcTrigonometric body forcing. Isoviscous.
Parameters: - sigma_0 (float) – Perturbation strength factor.
- n_x (int) – Wavenumber parameter (in x).
- n_z (float) – Wavenumber parameter (in z).
- eta_0 (float) – Viscosity.
-
class
underworld.function.analytic.SolKz(sigma_0=1.0, n_x=1, n_z=1.0, B=1.1512925465, *args, **kwargs)[source]¶ Bases:
underworld.function.analytic._SolBaseFreeSlipBcThe boundary conditions are free-slip everywhere on a unit domain. The viscosity varies exponentially in the z direction and is given by \(\eta = \exp (2 B z)\). The flow is driven by the following density perturbation:
Parameters: - sigma_0 (float) – Perturbation strength factor.
- n_x (int) – Wavenumber parameter (in x).
- n_z (float) – Wavenumber parameter (in z).
- B (float) – Viscosity parameter.
-
class
underworld.function.analytic.SolDB2d(*args, **kwargs)[source]¶ Bases:
underworld.function.analytic._SolBaseFixedBcSolDB2d and solDB3d from:
Dohrmann, C.R., Bochev, P.B., A stabilized finite element method for the Stokes problem based on polynomial pressure projections, Int. J. Numer. Meth. Fluids 46, 183-201 (2004).
Check get_bcs() for BC setup.
-
class
underworld.function.analytic.SolKx(sigma_0=1.0, n_x=1, n_z=1.0, B=1.1512925465, *args, **kwargs)[source]¶ Bases:
underworld.function.analytic._SolBaseFreeSlipBcThe boundary conditions are free-slip everywhere on a unit domain. The viscosity varies exponentially in the x direction and is given by \(\eta = \exp (2 B x)\). The flow is driven by the following density perturbation:
Parameters: - sigma_0 (float) – Perturbation strength factor.
- n_x (int) – Wavenumber parameter (in x).
- n_z (float) – Wavenumber parameter (in z).
- B (float) – Viscosity parameter.
-
class
underworld.function.analytic.SolNL(eta_0=1.0, n_z=1, r=1.5, *args, **kwargs)[source]¶ Bases:
underworld.function.analytic._SolBaseSolNL requires tighter solver tolerances and/or a direct solve for best results. Need to check in with Caesar Dandenonensis as to the origins of this solution.
Check get_bcs() for BC setup.
-
get_bcs(velVar)[source]¶ ( Fixed,Fixed ) conditions left/right. ( Free, Fixed ) conditions top/bottom.
All fixed DOFs set to analytic soln values.
Parameters: velVar (underworld.mesh.MeshVariable) – The velocity variable is required to construct the BC object. Returns: The BC object. It should be passed in to the system being constructed. Return type: underworld.conditions.SystemCondition
-
This module includes shape type functions. Shape functions generally define some geometric object, and return boolean values to indicate whether the queried locations are inside or outside the shape.
underworld.function.shape.Polygon |
This function creates a polygon shape. |
-
class
underworld.function.shape.Polygon(vertices, fn=None, *args, **kwargs)[source]¶ Bases:
underworld.function._function.FunctionThis function creates a polygon shape. Note that this is strictly a 2d shape, and the third dimension of any query will be ignored. You may create a box type function if you wish to limit the shape extent in the third dimension.
You will need to use rotations to orient the polygon in other directions. Rotations functions will be available shortly (hopefully!).
Parameters: - vertices (np.ndarray) – This array provides the vertices for the polygon. Note that the order of the vertices is important. The polygon is defined by a piecewise linear edge joining the vertices in the order provided by the array. The final vertex and the initial vertex are joined to complete the polygon.
- fn (underworld.function.Function, default=None) – This is the input function. Generally it will not be required, but you may need to use (for example) to transform the incoming coordinates.
Example
In this example we will create a triangle shape and test some points.
>>> import underworld as uw >>> import numpy as np
Create the array to define the triangle, and the function
>>> vertex_array = np.array( [(0.0,0.0),(0.5,1.0),(1.0,0.0)] ) >>> polyfn = uw.function.shape.Polygon(vertex_array)
Create some test points, and do a test evaluation
>>> test_array = np.array( [(0.0,0.9),(0.5,0.5),(0.9,0.2)] ) >>> polyfn.evaluate(test_array) array([[False], [ True], [False]], dtype=bool)
This module contains functions relating to rheological operations.
underworld.function.rheology.stress_limiting_viscosity |
Returns a viscosity value which effectively limits the maximum fluid stress. |
-
class
underworld.function.rheology.stress_limiting_viscosity(fn_stress, fn_stresslimit, fn_inputviscosity, *args, **kwargs)[source]¶ Bases:
underworld.function._function.FunctionReturns a viscosity value which effectively limits the maximum fluid stress. Where the stress invariant (as calculated using the provided fn_stress) is greater than the stress limit (as provided by the fn_stresslimit), the returned viscosity will affect a fluid stress at the stress limit. Otherwise, fn_inputviscosity is passed through.
Parameters: - fn_stress (underworld.function.Function) – Function which returns the current stress in the fluid. Function should return a symmetric tensor of floating point values.
- fn_stresslimit (underworld.function.Function) – Function which defines the stress limit. Function should return a scalar floating point value.
- fn_inputviscosity (underworld.function.Function) – Function which defines the non-yielded viscosity value. Function should return a scalar floating point value.
Example
Lets setup a simple shear type configuration but with a viscosity that increase vertically:
>>> import underworld as uw >>> import underworld.function as fn >>> mesh = uw.mesh.FeMesh_Cartesian(elementRes=(16,16), periodic=(True,False)) >>> velVar = uw.mesh.MeshVariable(mesh,2) >>> pressVar = uw.mesh.MeshVariable(mesh.subMesh,1)
Simple shear boundary conditions:
>>> bot_nodes = mesh.specialSets["MinJ_VertexSet"] >>> top_nodes = mesh.specialSets["MaxJ_VertexSet"] >>> bc = uw.conditions.DirichletCondition(velVar, (top_nodes+bot_nodes,top_nodes+bot_nodes)) >>> velVar.data[bot_nodes.data] = (-0.5,0.) >>> velVar.data[top_nodes.data] = ( 0.5,0.)
Vertically increasing exponential viscosity:
>>> fn_visc = 1. >>> stokesSys = uw.systems.Stokes(velVar,pressVar,fn_visc,conditions=[bc,])
Solve:
>>> solver = uw.systems.Solver(stokesSys) >>> solver.solve()
Use the min_max function to determine a maximum stress:
>>> fn_stress = 2.*fn_visc*uw.function.tensor.symmetric( velVar.fn_gradient ) >>> fn_minmax_inv = fn.view.min_max(fn.tensor.second_invariant(fn_stress)) >>> ignore = fn_minmax_inv.evaluate(mesh) >>> import numpy as np >>> np.allclose(fn_minmax_inv.max_global(), 1.0, rtol=1e-05) True
Now lets set the limited viscosity. Note that the system is now non-linear.
>>> fn_visc_limited = fn.rheology.stress_limiting_viscosity(fn_stress,0.5,fn_visc) >>> stokesSys.fn_viscosity = fn_visc_limited >>> solver.solve(nonLinearIterate=True)
Now check the stress:
>>> fn_stress = 2.*fn_visc_limited*uw.function.tensor.symmetric( velVar.fn_gradient ) >>> fn_minmax_inv = fn.view.min_max(fn.tensor.second_invariant(fn_stress)) >>> ignore = fn_minmax_inv.evaluate(mesh) >>> np.allclose(fn_minmax_inv.max_global(), 0.5, rtol=1e-05) True
This module provides math functions. All functions take functions (or convertibles) as arguments. These functions effectively wrap to the c++ standard template library equivalent. For example, the ‘exp’ class generates a function with uses std::exp(double).
All functions operate on and return ‘double’ type data (or ‘float’ from python).
underworld.function.math.pow |
Power function. |
underworld.function.math.cosh |
Computes the hyperbolic cosine of its argument function. |
underworld.function.math.acosh |
Computes the inverse hyperbolic cosine of its argument function. |
underworld.function.math.tan |
Computes the tangent of its argument function (measured in radians). |
underworld.function.math.asin |
Computes the principal value of the arc sine of x, expressed in radians. |
underworld.function.math.log |
Computes the natural logarithm of its argument function. |
underworld.function.math.atanh |
Computes the inverse hyperbolic tangent of its argument function. |
underworld.function.math.sqrt |
Computes the square root of its argument function. |
underworld.function.math.abs |
Computes the absolute value of its argument function. |
underworld.function.math.log10 |
Computes the base 10 logarithm of its argument function. |
underworld.function.math.sin |
Computes the sine of its argument function (measured in radians). |
underworld.function.math.asinh |
Computes the inverse hyperbolic sine of its argument function. |
underworld.function.math.log2 |
Computes the base 2 logarithm of its argument function. |
underworld.function.math.atan |
Computes the principal value of the arc tangent of x, expressed in radians. |
underworld.function.math.sinh |
Computes the hyperbolic sine of its argument function. |
underworld.function.math.cos |
Computes the cosine of its argument function (measured in radians). |
underworld.function.math.tanh |
Computes the hyperbolic tangent of its argument function. |
underworld.function.math.erf |
Computes the error function of its argument function. |
underworld.function.math.erfc |
Computes the complementary error function of its argument function. |
underworld.function.math.exp |
Computes the exponent of its argument function. |
underworld.function.math.acos |
Computes the principal value of the arc cosine of x, expressed in radians. |
underworld.function.math.dot |
Dot product function. |
-
class
underworld.function.math.pow(fn1, fn2, **kwargs)[source]¶ Bases:
underworld.function._function.FunctionPower function. Raises fn1 to the power of fn2.
Parameters: - fn1 (underworld.function.Function (or convertible)) – The base function.
- fn2 (underworld.function.Function (or convertible)) – The power function.
Example
>>> from . import _systemmath as math >>> import numpy as np >>> func = pow(_uw.function.input(),3.) >>> np.allclose( func.evaluate(2.), math.pow(2.,3.) ) True
-
class
underworld.function.math.cosh(fn=None, *args, **kwargs)[source]¶ Bases:
underworld.function._function.FunctionComputes the hyperbolic cosine of its argument function.
Parameters: fn (underworld.function.Function (or convertible)) – Optionally provided for function composition. Example
>>> from . import _systemmath as math >>> import numpy as np >>> func = cosh() >>> np.allclose( func.evaluate(0.1234), math.cosh(0.1234) ) True
-
class
underworld.function.math.acosh(fn=None, *args, **kwargs)[source]¶ Bases:
underworld.function._function.FunctionComputes the inverse hyperbolic cosine of its argument function.
Parameters: fn (underworld.function.Function (or convertible)) – Optionally provided for function composition. Example
>>> from . import _systemmath as math >>> import numpy as np >>> func = acosh() >>> np.allclose( func.evaluate(5.1234), math.acosh(5.1234) ) True
-
class
underworld.function.math.tan(fn=None, *args, **kwargs)[source]¶ Bases:
underworld.function._function.FunctionComputes the tangent of its argument function (measured in radians).
Parameters: fn (underworld.function.Function (or convertible)) – Optionally provided for function composition. Example
>>> from . import _systemmath as math >>> import numpy as np >>> func = tan() >>> np.allclose( func.evaluate(0.1234), math.tan(0.1234) ) True
-
class
underworld.function.math.asin(fn=None, *args, **kwargs)[source]¶ Bases:
underworld.function._function.FunctionComputes the principal value of the arc sine of x, expressed in radians.
Parameters: fn (underworld.function.Function (or convertible)) – Optionally provided for function composition. Example
>>> from . import _systemmath as math >>> import numpy as np >>> func = asin() >>> np.allclose( func.evaluate(0.1234), math.asin(0.1234) ) True
-
class
underworld.function.math.log(fn=None, *args, **kwargs)[source]¶ Bases:
underworld.function._function.FunctionComputes the natural logarithm of its argument function.
Parameters: fn (underworld.function.Function (or convertible)) – Optionally provided for function composition. Example
>>> from . import _systemmath as math >>> import numpy as np >>> func = log() >>> np.allclose( func.evaluate(0.1234), math.log(0.1234) ) True
-
class
underworld.function.math.atanh(fn=None, *args, **kwargs)[source]¶ Bases:
underworld.function._function.FunctionComputes the inverse hyperbolic tangent of its argument function.
Parameters: fn (underworld.function.Function (or convertible)) – Optionally provided for function composition. Example
>>> from . import _systemmath as math >>> import numpy as np >>> func = atanh() >>> np.allclose( func.evaluate(0.1234), math.atanh(0.1234) ) True
-
class
underworld.function.math.sqrt(fn=None, *args, **kwargs)[source]¶ Bases:
underworld.function._function.FunctionComputes the square root of its argument function.
Parameters: fn (underworld.function.Function (or convertible)) – Optionally provided for function composition. Example
>>> from . import _systemmath as math >>> import numpy as np >>> func = sqrt() >>> np.allclose( func.evaluate(0.1234), math.sqrt(0.1234) ) True
-
class
underworld.function.math.abs(fn=None, *args, **kwargs)[source]¶ Bases:
underworld.function._function.FunctionComputes the absolute value of its argument function.
Parameters: fn (underworld.function.Function (or convertible)) – Optionally provided for function composition. Example
>>> from . import _systemmath as math >>> import numpy as np >>> func = abs() >>> np.allclose( func.evaluate(-0.1234), math.fabs(0.1234) ) True
-
class
underworld.function.math.log10(fn=None, *args, **kwargs)[source]¶ Bases:
underworld.function._function.FunctionComputes the base 10 logarithm of its argument function.
Parameters: fn (underworld.function.Function (or convertible)) – Optionally provided for function composition. Example
>>> from . import _systemmath as math >>> import numpy as np >>> func = log10() >>> np.allclose( func.evaluate(0.1234), math.log10(0.1234) ) True
-
class
underworld.function.math.sin(fn=None, *args, **kwargs)[source]¶ Bases:
underworld.function._function.FunctionComputes the sine of its argument function (measured in radians).
Parameters: fn (underworld.function.Function (or convertible)) – Optionally provided for function composition. Example
>>> from . import _systemmath as math >>> import numpy as np >>> func = sin() >>> np.allclose( func.evaluate(0.1234), math.sin(0.1234) ) True
-
class
underworld.function.math.asinh(fn=None, *args, **kwargs)[source]¶ Bases:
underworld.function._function.FunctionComputes the inverse hyperbolic sine of its argument function.
Parameters: fn (underworld.function.Function (or convertible)) – Optionally provided for function composition. Example
>>> from . import _systemmath as math >>> import numpy as np >>> func = asinh() >>> np.allclose( func.evaluate(5.1234), math.asinh(5.1234) ) True
-
class
underworld.function.math.log2(fn=None, *args, **kwargs)[source]¶ Bases:
underworld.function._function.FunctionComputes the base 2 logarithm of its argument function.
Parameters: fn (underworld.function.Function (or convertible)) – Optionally provided for function composition. Example
>>> from . import _systemmath as math >>> import numpy as np >>> func = log2() >>> np.allclose( func.evaluate(0.1234), math.log(0.1234,2) ) True
-
class
underworld.function.math.atan(fn=None, *args, **kwargs)[source]¶ Bases:
underworld.function._function.FunctionComputes the principal value of the arc tangent of x, expressed in radians.
Parameters: fn (underworld.function.Function (or convertible)) – Optionally provided for function composition. Example
>>> from . import _systemmath as math >>> import numpy as np >>> func = atan() >>> np.allclose( func.evaluate(0.1234), math.atan(0.1234) ) True
-
class
underworld.function.math.sinh(fn=None, *args, **kwargs)[source]¶ Bases:
underworld.function._function.FunctionComputes the hyperbolic sine of its argument function.
Parameters: fn (underworld.function.Function (or convertible)) – Optionally provided for function composition. Example
>>> from . import _systemmath as math >>> import numpy as np >>> func = sinh() >>> np.allclose( func.evaluate(0.1234), math.sinh(0.1234) ) True
-
class
underworld.function.math.cos(fn=None, *args, **kwargs)[source]¶ Bases:
underworld.function._function.FunctionComputes the cosine of its argument function (measured in radians).
Parameters: fn (underworld.function.Function (or convertible)) – Optionally provided for function composition. Example
>>> from . import _systemmath as math >>> import numpy as np >>> func = cos() >>> np.allclose( func.evaluate(0.1234), math.cos(0.1234) ) True
-
class
underworld.function.math.tanh(fn=None, *args, **kwargs)[source]¶ Bases:
underworld.function._function.FunctionComputes the hyperbolic tangent of its argument function.
Parameters: fn (underworld.function.Function (or convertible)) – Optionally provided for function composition. Example
>>> from . import _systemmath as math >>> import numpy as np >>> func = tanh() >>> np.allclose( func.evaluate(0.1234), math.tanh(0.1234) ) True
-
class
underworld.function.math.erf(fn=None, *args, **kwargs)[source]¶ Bases:
underworld.function._function.FunctionComputes the error function of its argument function.
Parameters: fn (underworld.function.Function (or convertible)) – Optionally provided for function composition. Example
>>> from . import _systemmath as math >>> import numpy as np >>> func = erf() >>> np.allclose( func.evaluate(0.1234), math.erf(0.1234) ) True
-
class
underworld.function.math.erfc(fn=None, *args, **kwargs)[source]¶ Bases:
underworld.function._function.FunctionComputes the complementary error function of its argument function.
Parameters: fn (underworld.function.Function (or convertible)) – Optionally provided for function composition. Example
>>> from . import _systemmath as math >>> import numpy as np >>> func = erfc() >>> np.allclose( func.evaluate(0.1234), math.erfc(0.1234) ) True
-
class
underworld.function.math.exp(fn=None, *args, **kwargs)[source]¶ Bases:
underworld.function._function.FunctionComputes the exponent of its argument function.
Parameters: fn (underworld.function.Function (or convertible)) – Optionally provided for function composition. Example
>>> from . import _systemmath as math >>> import numpy as np >>> func = exp() >>> np.allclose( func.evaluate(0.1234), math.exp(0.1234) ) True
-
class
underworld.function.math.acos(fn=None, *args, **kwargs)[source]¶ Bases:
underworld.function._function.FunctionComputes the principal value of the arc cosine of x, expressed in radians.
Parameters: fn (underworld.function.Function (or convertible)) – Optionally provided for function composition. Example
>>> from . import _systemmath as math >>> import numpy as np >>> func = acos() >>> np.allclose( func.evaluate(0.1234), math.acos(0.1234) ) True
-
class
underworld.function.math.dot(fn1, fn2, **kwargs)[source]¶ Bases:
underworld.function._function.FunctionDot product function. Returns fn1.fn2. Argument functions must return values of identical size.
Parameters: - fn1 (underworld.function.Function (or convertible)) – Argument function 1.
- fn2 (underworld.function.Function (or convertible)) – Argument function 2.
Example
>>> from . import _systemmath as math >>> import numpy as np >>> input1 = (2.,3.,4.) >>> input2 = (5.,6.,7.) >>> func = dot( input1, input2 )
The function is constant, so evaluate anywhere:
>>> np.allclose( func.evaluate(0.), np.dot(input1,input2) ) True
This module includes functions which provide views into the results of function queries. These functions never modify query data.
underworld.function.view.min_max |
This function records the min & max result from a queried function. |
-
class
underworld.function.view.min_max(fn, fn_norm=None, fn_auxiliary=None, *args, **kwargs)[source]¶ Bases:
underworld.function._function.FunctionThis function records the min & max result from a queried function.
Note that this function simply records the min/max values encountered when it is evaluated. Therefore, if it has not been evaluated at all, the values returned via one of its methods (‘min_local’, ‘min_global’, etc ) will simply be initialisation values.
For vector input types, this function will report on the magnitude of the vector.
Parameters: - fn (underworld.function.Function) – The primary function. If fn_norm is not provided, this is used to calculate the min_max. Results from this function are always passed back.
- fn_norm (underworld.function.Function) – This function returns a norm like quantity by which the min and max are determined. For example, where the primary function is a vector quantity, this function might calculate the magnitude of that vector. This function must return a scalar result, and must be provided where the primary function is non-scalar. See the example below for usage.
- fn_auxiliary (underworld.function.Function) – An auxiliary function which is evaluated at the location of the min/max. For example, often the coordinate where the min/max values occur are required, and so the user may pass in fn.input() as the auxiliary function to achieve this
Example
Create a simple function which returns two times its input:
>>> import underworld as uw >>> import underworld.function as fn >>> import numpy as np >>> fn_simple = fn.input()[0]*2.
Let’s wrap it with a min_max function:
>>> fn_minmax_simple = fn.view.min_max(fn_simple)
Now do an evaluation:
>>> fn_minmax_simple.evaluate(5.) array([[ 10.]])
Since there’s only been one evaluation, min and max values should be identical:
>>> fn_minmax_simple.min_global() 10.0 >>> fn_minmax_simple.max_global() 10.0
Do another evaluation:
>>> fn_minmax_simple.evaluate(-3.) array([[-6.]])
Now check min and max again:
>>> fn_minmax_simple.min_global() -6.0 >>> fn_minmax_simple.max_global() 10.0
Note that if we only evaluate the subject function, no min/max values are recorded:
>>> fn_simple.evaluate(3000.) array([[ 6000.]]) >>> fn_minmax_simple.max_global() 10.0
Also note that for vector valued subject function, fn_norm must be provided:
>>> fn_vec = fn.input()*(1.,1.) >>> fn_vec_mm = fn.view.min_max(fn_vec) >>> fn_vec_mm.evaluate( 2. ) Traceback (most recent call last): ... RuntimeError: Issue utilising function of class 'min_max' constructed at: --- CONSTRUCTION TIME STACK --- Error message: Argument function does not return scalar results. You must also provide a function which calculates the required norm like quantity via the `fn_norm` parameter.
>>> fn_vec_mm = fn.view.min_max(fn_vec, fn_norm=fn.math.dot(fn_vec,fn_vec)) >>> fn_vec_mm.evaluate( 2. ) array([[ 2., 2.]]) >>> fn_vec_mm.max_global() 8.0 >>> fn_vec_mm.evaluate( -1. ) array([[-1., -1.]]) >>> fn_vec_mm.min_global() 2.0 >>> fn_vec_mm.max_global() 8.0
To obtain the min/max values across a MeshVariable object, you will need to evaluate the function across all nodes of the MeshVariable:
>>> mesh = uw.mesh.FeMesh_Cartesian() >>> meshvariable = uw.mesh.MeshVariable( mesh, 1 ) >>> meshvariable.data[:] = np.random.randint(100,size=meshvariable.data.shape) # init with random data >>> fn_mv = fn.view.min_max(meshvariable) # create min_max view wrapper >>> ignore = fn_mv.evaluate(mesh) # this call will evaluate at all nodes >>> np.allclose(fn_mv.min_local(),meshvariable.data.min()) True >>> np.allclose(fn_mv.max_local(),meshvariable.data.max()) True
Note that when operating in parallel, the min_global() and max_global() methods are a good option for extracting discrete object global min/max values, as the numpy views will only report the local min/max values.
Also note that since min_max views only record results as they are evaluated, if the underlying subject function min/max values change, this will not be recorded by the min_max view until its evaluate encounters the new min/max values:
>>> meshvariable.data[3] = 1000 # change some random value >>> np.allclose(fn_mv.max_local(),meshvariable.data.max()) # check again, it should be false False >>> ignore = fn_mv.evaluate(mesh) # evaluate across all nodes again >>> np.allclose(fn_mv.max_local(),meshvariable.data.max()) # check again True
Similarly, the view’s min/max values are only updated when smaller/larger min/max values are encountered. So, if the underlying subject function’s maximum (for example) is reduced, the view will not record this if its currently stored value exceeds the new maximum. A call to reset() is required:
>>> fn_mv.max_local() 1000.0 >>> meshvariable.data[3] = 500 # reduce max >>> ignore = fn_mv.evaluate(mesh) # evaluate across all nodes again >>> fn_mv.max_local() # note that it still records old value 1000.0 >>> fn_mv.reset() # now re-init view's min/max >>> ignore = fn_mv.evaluate(mesh) # evaluate across all nodes again >>> fn_mv.max_local() # it should now record new value 500.0
The auxiliary function allows you to obtain secondary information at the function minimum. One common use case would be to obtain a location where the min/max was obtained:
>>> fn_mv = fn.view.min_max(meshvariable, fn_auxiliary=fn.input()) >>> meshvariable.data[1] = 1000.0 # set second node to have the highest value >>> ignore = fn_mv.evaluate(mesh) >>> fn_mv.max_global() 1000.0 >>> np.allclose( mesh.data[1], fn_mv.max_global_auxiliary() ) # ensure max is obtain at required mesh node. True
-
max_global()[source]¶ Returns the maximum value encountered across all processes.
Notes
This method must be called by collectively all processes.
Returns: double Return type: maximum value
-
max_global_auxiliary()[source]¶ Returns the results of the auxiliary function evaluated at the location corresponding to the primary function maximum. This method considers results across all processes (ie, globally).
Notes
This method must be called by collectively all processes.
Returns: FunctionIO Return type: value at global maximum.
-
max_local()[source]¶ Returns the max value encountered locally on the current process.
Returns: double Return type: maximum value
-
max_local_auxiliary()[source]¶ Returns the results of the auxiliary function evaluated at the location corresponding to the primary function maximum. This method only considers results on the current process.
Returns: FunctionIO Return type: value at local maximum.
-
max_rank()[source]¶ Returns the rank where the maximum occurs. Note that this method will return -1 until max_global has been called.
Returns: int Return type: rank
-
min_global()[source]¶ Returns the minimum value encountered across all processes.
Notes
This method must be called by collectively all processes.
Returns: double Return type: minimum value
-
min_global_auxiliary()[source]¶ Returns the results of the auxiliary function evaluated at the location corresponding to the primary function minimum. This method considers results across all processes (ie, globally).
Notes
This method must be called by collectively all processes.
Returns: FunctionIO Return type: value at global minimum.
-
min_local()[source]¶ Returns the minimum value encountered locally on the current process.
Returns: double Return type: minimum value
-
min_local_auxiliary()[source]¶ Returns the results of the auxiliary function evaluated at the location corresponding to the primary function minimum. This method only considers results on the current process.
Returns: FunctionIO Return type: value at local minimum.
underworld.function.Function |
Objects which inherit from this class provide user definable functions within Underworld. |
underworld.function.FunctionInput |
Objects that inherit from this class are able to act as inputs to function evaluation from python. |
underworld.function.coord |
alias of underworld.function._function.input |
underworld.function.input |
This class generates a function which simply passes through its input. |
-
class
underworld.function.Function(argument_fns, **kwargs)[source]¶ Bases:
underworld._stgermain.LeftOverParamsCheckerObjects which inherit from this class provide user definable functions within Underworld.
Functions aim to achieve a number of goals: * Provide a natural interface for mathematical behaviour description within python. * Provide a high level interface to Underworld discrete objects. * Allow discrete objects to be used in combination with continuous objects. * Handle the evaluation of discrete objects in the most efficient manner. * Perform all heavy calculations at the C-level for efficiency. * Provide an interface for users to evaluate functions directly within python, utilising numpy arrays for input/output.
-
__add__(other)[source]¶ Operator overloading for ‘+’ operation:
Fn3 = Fn1 + Fn2
Creates a new function Fn3 which performs additions of Fn1 and Fn2.
Returns: fn.add Return type: Add function Examples
>>> import misc >>> import numpy as np >>> three = misc.constant(3.) >>> four = misc.constant(4.) >>> np.allclose( (three + four).evaluate(0.), [[ 7.]] ) # note we can evaluate anywhere because it's a constant True
-
__and__(other)[source]¶ Operator overloading for ‘&’ operation:
Fn3 = Fn1 & Fn2
Creates a new function Fn3 which returns a bool result for the operation.
Returns: fn.logical_and Return type: AND function Examples
>>> import misc >>> trueFn = misc.constant(True) >>> falseFn = misc.constant(False) >>> (trueFn & falseFn).evaluate() array([[False]], dtype=bool)
Notes
The ‘&’ operator in python is usually used for bitwise ‘and’ operations, with the ‘and’ operator used for boolean type operators. It is not possible to overload the ‘and’ operator in python, so instead the bitwise equivalent has been utilised.
-
__div__(other)[source]¶ Operator overloading for ‘/’ operation:
Fn3 = Fn1 / Fn2
Creates a new function Fn3 which returns the quotient of Fn1 and Fn2.
Returns: fn.divide Return type: Divide function Examples
>>> import misc >>> import numpy as np >>> two = misc.constant(2.) >>> four = misc.constant(4.) >>> np.allclose( (four/two).evaluate(0.), [[ 2.]] ) # note we can evaluate anywhere because it's a constant True
-
__ge__(other)[source]¶ Operator overloading for ‘>=’ operation:
Fn3 = Fn1 >= Fn2
Creates a new function Fn3 which returns a bool result for the relation.
Returns: fn.greater_equal Return type: Greater than or equal to function Examples
>>> import misc >>> import numpy as np >>> two = misc.constant(2.) >>> (two >= two).evaluate() array([[ True]], dtype=bool)
-
__getitem__(index)[source]¶ Operator overloading for ‘[]’ operation:
FnComponent = Fn[0]
Creates a new function FnComponent which returns the required component of Fn.
Returns: fn.at Return type: component function Examples
>>> import misc >>> fn = misc.constant((2.,3.,4.)) >>> np.allclose( fn[1].evaluate(0.), [[ 3.]] ) # note we can evaluate anywhere because it's a constant True
-
__gt__(other)[source]¶ Operator overloading for ‘>’ operation:
Fn3 = Fn1 > Fn2
Creates a new function Fn3 which returns a bool result for the relation.
Returns: fn.greater Return type: Greater than function Examples
>>> import misc >>> import numpy as np >>> two = misc.constant(2.) >>> four = misc.constant(4.) >>> (two > four).evaluate() array([[False]], dtype=bool)
-
__le__(other)[source]¶ Operator overloading for ‘<=’ operation:
Fn3 = Fn1 <= Fn2
Creates a new function Fn3 which returns a bool result for the relation.
Returns: fn.less_equal Return type: Less than or equal to function Examples
>>> import misc >>> import numpy as np >>> two = misc.constant(2.) >>> (two <= two).evaluate() array([[ True]], dtype=bool)
-
__lt__(other)[source]¶ Operator overloading for ‘<’ operation:
Fn3 = Fn1 < Fn2
Creates a new function Fn3 which returns a bool result for the relation.
Returns: fn.less Return type: Less than function Examples
>>> import misc >>> import numpy as np >>> two = misc.constant(2.) >>> four = misc.constant(4.) >>> (two < four).evaluate() array([[ True]], dtype=bool)
-
__mul__(other)[source]¶ Operator overloading for ‘*’ operation:
Fn3 = Fn1 * Fn2
Creates a new function Fn3 which returns the product of Fn1 and Fn2.
Returns: fn.multiply Return type: Multiply function Examples
>>> import misc >>> import numpy as np >>> three = misc.constant(3.) >>> four = misc.constant(4.) >>> np.allclose( (three*four).evaluate(0.), [[ 12.]] ) # note we can evaluate anywhere because it's a constant True
-
__neg__()[source]¶ Operator overloading for unary ‘-‘.
FnNeg = -Fn
Creates a new function FnNeg which is the negative of Fn.
Returns: fn.multiply Return type: Negative function Examples
>>> import misc >>> import numpy as np >>> four = misc.constant(4.) >>> np.allclose( (-four).evaluate(0.), [[ -4.]] ) # note we can evaluate anywhere because it's a constant True
-
__or__(other)[source]¶ Operator overloading for ‘|’ operation:
Fn3 = Fn1 | Fn2
Creates a new function Fn3 which returns a bool result for the operation.
Returns: fn.logical_or Return type: OR function Examples
>>> import misc >>> trueFn = misc.constant(True) >>> falseFn = misc.constant(False) >>> (trueFn | falseFn).evaluate() array([[ True]], dtype=bool)
Notes
The ‘|’ operator in python is usually used for bitwise ‘or’ operations, with the ‘or’ operator used for boolean type operators. It is not possible to overload the ‘or’ operator in python, so instead the bitwise equivalent has been utilised.
-
__pow__(other)[source]¶ Operator overloading for ‘**’ operation:
Fn3 = Fn1 ** Fn2
Creates a new function Fn3 which returns Fn1 to the power of Fn2.
Returns: fn.math.pow Return type: Power function Examples
>>> import misc >>> import numpy as np >>> two = misc.constant(2.) >>> four = misc.constant(4.) >>> np.allclose( (two**four).evaluate(0.), [[ 16.]] ) # note we can evaluate anywhere because it's a constant True
-
__radd__(other)¶ Operator overloading for ‘+’ operation:
Fn3 = Fn1 + Fn2
Creates a new function Fn3 which performs additions of Fn1 and Fn2.
Returns: fn.add Return type: Add function Examples
>>> import misc >>> import numpy as np >>> three = misc.constant(3.) >>> four = misc.constant(4.) >>> np.allclose( (three + four).evaluate(0.), [[ 7.]] ) # note we can evaluate anywhere because it's a constant True
-
__rmul__(other)¶ Operator overloading for ‘*’ operation:
Fn3 = Fn1 * Fn2
Creates a new function Fn3 which returns the product of Fn1 and Fn2.
Returns: fn.multiply Return type: Multiply function Examples
>>> import misc >>> import numpy as np >>> three = misc.constant(3.) >>> four = misc.constant(4.) >>> np.allclose( (three*four).evaluate(0.), [[ 12.]] ) # note we can evaluate anywhere because it's a constant True
-
__rsub__(other)[source]¶ Operator overloading for ‘-‘ operation. Right hand version.
Fn3 = Fn1 - Fn2
Creates a new function Fn3 which performs subtraction of Fn2 from Fn1.
Returns: fn.subtract Return type: RHS subtract function Examples
>>> import misc >>> import numpy as np >>> four = misc.constant(4.) >>> np.allclose( (5. - four).evaluate(0.), [[ 1.]] ) # note we can evaluate anywhere because it's a constant True
-
__sub__(other)[source]¶ Operator overloading for ‘-‘ operation:
Fn3 = Fn1 - Fn2
Creates a new function Fn3 which performs subtraction of Fn2 from Fn1.
Returns: fn.subtract Return type: Subtract function Examples
>>> import misc >>> import numpy as np >>> three = misc.constant(3.) >>> four = misc.constant(4.) >>> np.allclose( (three - four).evaluate(0.), [[ -1.]] ) # note we can evaluate anywhere because it's a constant True
-
__xor__(other)[source]¶ Operator overloading for ‘^’ operation:
Fn3 = Fn1 ^ Fn2
Creates a new function Fn3 which returns a bool result for the operation.
Returns: fn.logical_xor Return type: XOR function Examples
>>> import misc >>> trueFn = misc.constant(True) >>> falseFn = misc.constant(False) >>> (trueFn ^ falseFn).evaluate() array([[ True]], dtype=bool) >>> (trueFn ^ trueFn).evaluate() array([[False]], dtype=bool) >>> (falseFn ^ falseFn).evaluate() array([[False]], dtype=bool)
Notes
The ‘^’ operator in python is usually used for bitwise ‘xor’ operations, however here we always use the logical version, with the operation inputs cast to their bool equivalents before the operation.
-
static
convert(obj)[source]¶ This method will attempt to convert the provided input into an equivalent underworld function. If the provided input is already of Function type, it is immediately returned. Likewise, if the input is of None type, it is also returned.
Parameters: obj (fn_like) – The object to be converted. Note that if obj is of type None or Function, it is simply returned immediately. Where obj is of type int/float/double, a Constant type function is returned which evaluates to the provided object’s value. Where obj is of type list/tuple, a function will be returned which evaluates to a vector of the provided list/tuple’s values (where possible). Returns: Return type: Fn.Function or None. Examples
>>> import underworld as uw >>> import underworld.function as fn
>>> fn_const = fn.Function.convert( 3 ) >>> fn_const.evaluate(0.) # eval anywhere for constant array([[3]], dtype=int32)
>>> fn_const == fn.Function.convert( fn_const ) True
>>> fn.Function.convert( None )
>>> fn1 = fn.input() >>> fn2 = 10.*fn.input() >>> fn3 = 100.*fn.input() >>> vec = (fn1,fn2,fn3) >>> fn_vec = fn.Function.convert(vec) >>> fn_vec.evaluate([3.]) array([[ 3., 30., 300.]])
-
evaluate(inputData=None, inputType=None)[source]¶ This method performs evaluate of a function at the given input(s).
It accepts floats, lists, tuples, numpy arrays, or any object which is of class FunctionInput. lists/tuples must contain floats only.
FunctionInput class objects are shortcuts to their underlying data, often with performance advantages, and sometimes they are the only valid input type (such as using Swarm objects as an inputs to SwarmVariable evaluation). Objects of class FeMesh, Swarm, FeMesh_IndexSet and VoronoiIntegrationSwarm are also of class FunctionInput. See the Function section of the user guide for more information.
Results are returned as numpy array.
Parameters: - inputData (float, list, tuple, ndarray, underworld.function.FunctionInput) – The input to the function. The form of this input must be appropriate for the function being evaluated, or an exception will be thrown. Note that if no input is provided, function will be evaluated at 0.
- inputType (str) – Specifies the type the provided data represents. Acceptable values are ‘scalar’, ‘vector’, ‘symmetrictensor’, ‘tensor’, ‘array’.
Returns: ndarray
Return type: array of results
Examples
>>> from . import _systemmath as math >>> import underworld.function.math as fnmath >>> sinfn = fnmath.sin()
Single evaluation:
>>> np.allclose( sinfn.evaluate(math.pi/4.), [[ 0.5*math.sqrt(2.)]] ) True
Multiple evaluations
>>> input = (0.,math.pi/4.,2.*math.pi) >>> np.allclose( sinfn.evaluate(input), [[ 0., 0.5*math.sqrt(2.), 0.]] ) True
Single MeshVariable evaluations
>>> mesh = uw.mesh.FeMesh_Cartesian() >>> var = uw.mesh.MeshVariable(mesh,1) >>> import numpy as np >>> var.data[:,0] = np.linspace(0,1,len(var.data)) >>> result = var.evaluate( (0.2,0.5 ) ) >>> np.allclose( result, np.array([[ 0.45]]) ) True
Numpy input MeshVariable evaluation
>>> # evaluate at a set of locations.. provide these as a numpy array. >>> count = 10 >>> # create an empty array >>> locations = np.zeros( (count,2)) >>> # specify evaluation coodinates >>> locations[:,0] = 0.5 >>> locations[:,1] = np.linspace(0.,1.,count) >>> # evaluate >>> result = var.evaluate(locations) >>> np.allclose( result, np.array([[ 0.08333333], [ 0.17592593], [ 0.26851852], [ 0.36111111], [ 0.4537037 ], [ 0.5462963 ], [ 0.63888889], [ 0.73148148], [ 0.82407407], [ 0.91666667]]) ) True
Using the mesh object as a FunctionInput
>>> np.allclose( var.evaluate(mesh), var.evaluate(mesh.data)) True
-
evaluate_global(inputData, inputType=None)[source]¶ This method attempts to evalute inputData across all processes, and then consolide the results on the root processor. This is most useful where you wish to evalute your functions using global coordinates which may span processes in a parallel simulation.
Note that this method does not currently support ‘FunctionInput’ class input data.
Due to the communications required for this method, a significant performance overhead may be encountered. The standard evaluate method should be used instead wherever possible.
Please see evaluate method for parameter details.
Notes
This method must be called collectively by all processes.
Returns: - Only the root process gets the final results array. All other processes
- are returned None.
-
integrate(mesh)[source]¶ Perform an integral of this underworld function over the given mesh
Parameters: mesh (uw.mesh.FeMesh_Cartesian) – Domain to perform integral over. Examples
>>> mesh = uw.mesh.FeMesh_Cartesian(minCoord=(0.0,0.0), maxCoord=(1.0,2.0)) >>> fn_1 = uw.function.misc.constant(2.0) >>> np.allclose( fn_1.integrate( mesh )[0], 4 ) True
>>> fn_2 = uw.function.misc.constant(2.0) * (0.5, 1.0) >>> np.allclose( fn_2.integrate( mesh ), [2,4] ) True
-
-
class
underworld.function.FunctionInput(*args, **kwargs)[source]¶ Bases:
underworld._stgermain.LeftOverParamsCheckerObjects that inherit from this class are able to act as inputs to function evaluation from python.
-
underworld.function.coord¶ alias of
underworld.function._function.input
-
class
underworld.function.input(*args, **kwargs)[source]¶ Bases:
underworld.function._function.FunctionThis class generates a function which simply passes through its input. It is the identity function. It is often useful when construct functions where the input itself needs to be accessed, such as to extract a particular component.
For example, you may wish to use this function when you wish to extract a particular coordinate component for manipulation. For this reason, we also provide an alias to this class called ‘coord’.
Returns: fn.input Return type: the input function Examples
Here we see the input function simply passing through its input.
>>> infunc = input() >>> np.allclose( infunc.evaluate( (1.,2.,3.) ), [ 1., 2., 3.] ) True
Often this behaviour is useful when we want to construct a function which operates on only a particular coordinate, such as a depth dependent density. We may wish to extract the z coordinate (in 2d):
>>> zcoord = input()[1] >>> baseDensity = 1. >>> density = baseDensity - 0.01*zcoord >>> testCoord1 = (0.1,0.4) >>> testCoord2 = (0.9,0.4) >>> np.allclose( density.evaluate( testCoord1 ), density.evaluate( testCoord2 ) ) True
underworld.container module¶
Implementation relating to container objects.
underworld.container.IndexSet |
The IndexSet class provides a set type container for integer values. |
underworld.container.ObjectifiedIndexSet |
This class simply adds an object to IndexSet data. |
-
class
underworld.container.IndexSet(size, fromObject=None)[source]¶ Bases:
objectThe IndexSet class provides a set type container for integer values. The underlying implementation is designed for memory efficiency. Index insertion and removal is a constant time operation.
Parameters: - size (int) – The size of the IndexSet. Note that the size corresponds to the maximum index value (plus 1) the set is required to hold, NOT the number of elements in the set. See IndexSet.size docstring for more information.
- fromObject (iterable, array_like, IndexSet. Optional.) – If provided, an IndexSet will be constructed using provided object’s data. See ‘add’ method for more details on acceptable objects. If not provided, empty set is generated.
Examples
You can add items via the constructor:
>>> someSet = uw.container.IndexSet( 15, [3,14,2] ) >>> someSet IndexSet([ 2, 3, 14])
Alternatively, create an empty set and add items as necessary:
>>> someSet = uw.container.IndexSet( 15 ) >>> someSet IndexSet([]) >>> someSet.add(3) >>> someSet IndexSet([3]) >>> someSet.add( [2,11] ) >>> someSet IndexSet([ 2, 3, 11])
Python operators are overloaded for convenience. Check class method details for full details.
-
AND(indices)[source]¶ Logical AND operation performed with provided IndexSet.
Parameters: indices (IndexSet) – IndexSet for which AND operation is performed. Note that provided set must be of type IndexSet. Example
>>> someSet1 = uw.container.IndexSet( 15, [3,9,10] ) >>> someSet2 = uw.container.IndexSet( 15, [1,9,12] ) >>> someSet1.AND(someSet2) >>> someSet1 IndexSet([9])
-
__add__(other)[source]¶ Operator overloading for
C = A + BCreates a new set C, then adds indices from A and B.
Returns: indexSet – The new set (C). Return type: IndexSet Example
>>> someSet1 = uw.container.IndexSet( 15, [3,9,10] ) >>> someSet2 = uw.container.IndexSet( 15, [1,9,12] ) >>> someSet1 + someSet2 IndexSet([ 1, 3, 9, 10, 12])
-
__and__(other)[source]¶ Operator overloading for
C = A & BCreates a new set C, then adds indices from A, and performs AND logic with B.
Returns: indexSet – The new set (C). Return type: IndexSet Example
>>> someSet1 = uw.container.IndexSet( 15, [3,9,10] ) >>> someSet2 = uw.container.IndexSet( 15, [1,9,12] ) >>> someSet1 & someSet2 IndexSet([9])
-
__contains__(index)[source]¶ Check if item is in IndexSet.
Parameters: index (unsigned int) – Check if index is in IndexSet. Returns: inSet – True if item is in set, False otherwise. Return type: bool Example
>>> someSet = uw.container.IndexSet( 15, [3,9,10] ) >>> 3 in someSet True
-
__deepcopy__(memo)[source]¶ Custom deepcopy routine required because python won’t know how to copy memory owned by stgermain.
-
__iadd__(other)[source]¶ Operator overloading for
A += BAdds indices from A and B.
Example
>>> someSet1 = uw.container.IndexSet( 15, [3,9,10] ) >>> someSet2 = uw.container.IndexSet( 15, [1,9,12] ) >>> someSet1 += someSet2 >>> someSet1 IndexSet([ 1, 3, 9, 10, 12])
-
__iand__(other)[source]¶ Operator overloading for
A &= BPerforms logical AND operation with A and B. Results are stored in A.
Example
>>> someSet1 = uw.container.IndexSet( 15, [3,9,10] ) >>> someSet2 = uw.container.IndexSet( 15, [1,9,12] ) >>> someSet1 &= someSet2 >>> someSet1 IndexSet([9])
-
__ior__(other)[source]¶ Operator overloading for
A |= BPerforms logical OR operation with A and B. Results are stored in A.
Example
>>> someSet1 = uw.container.IndexSet( 15, [3,9,10] ) >>> someSet2 = uw.container.IndexSet( 15, [1,9,12] ) >>> someSet1 |= someSet2 >>> someSet1 IndexSet([ 1, 3, 9, 10, 12])
-
__isub__(other)[source]¶ Operator overloading for
A -= BRemoves from A indices in B.
Example
>>> someSet1 = uw.container.IndexSet( 15, [3,9,10] ) >>> someSet2 = uw.container.IndexSet( 15, [1,9,12] ) >>> someSet1 -= someSet2 >>> someSet1 IndexSet([ 3, 10])
-
__len__()[source]¶ Overload for Python len usage.
Returns: int – Returns the total number of members this set contains. Return type: member count Example
>>> someSet = uw.container.IndexSet( 15, [3,9,10] ) >>> len(someSet) 3
-
__or__(other)[source]¶ Operator overloading for
C = A | BCreates a new set C, then adds indices from A, and performs OR logic with B.
Returns: indexSet – The new set (C). Return type: IndexSet Example
>>> someSet1 = uw.container.IndexSet( 15, [3,9,10] ) >>> someSet2 = uw.container.IndexSet( 15, [1,9,12] ) >>> someSet1 | someSet2 IndexSet([ 1, 3, 9, 10, 12])
-
__sub__(other)[source]¶ Operator overloading for
C = A - BCreates a new set C, then adds indices from A, and removes those from B.
Returns: indexSet – The new set (C). Return type: IndexSet Example
>>> someSet1 = uw.container.IndexSet( 15, [3,9,10] ) >>> someSet2 = uw.container.IndexSet( 15, [1,9,12] ) >>> someSet1 - someSet2 IndexSet([ 3, 10])
-
add(indices)[source]¶ Add item(s) to IndexSet.
Parameters: indices (unsigned int, ndarray, IndexSet, iterable object.) – Index or indices to be added to the IndexSet. Ensure value(s) are integer and non-negative. An iterable object may also be provided, with numpy arrays and IndexSets being significantly more efficient. Example
Create an empty set and add items as necessary:
>>> someSet = uw.container.IndexSet( 15 ) >>> someSet.add(3) >>> someSet IndexSet([3]) >>> 3 in someSet True >>> someSet.add([5,3,7,8]) >>> someSet IndexSet([3, 5, 7, 8]) >>> someSet.add(np.array([10,11,3])) >>> someSet IndexSet([ 3, 5, 7, 8, 10, 11])
-
addAll()[source]¶ Set all indices of set to added.
Example
>>> someSet = uw.container.IndexSet( 5 ) >>> someSet IndexSet([]) >>> someSet.addAll() >>> someSet IndexSet([0, 1, 2, 3, 4])
-
clear()[source]¶ Clear set. ie, set all indices to not included.
Example
>>> someSet = uw.container.IndexSet( 5, [1,2,3] ) >>> someSet IndexSet([1, 2, 3]) >>> someSet.clear() >>> someSet IndexSet([])
-
data¶ Returns the set members as a numpy array.
Note that only a numpy copy of the set is returned, and modifying this array is disabled (and would have no effect).
Returns: Array containing IndexSet members. Return type: numpy.ndarray (uint32) Example
>>> someSet = uw.container.IndexSet( 15, [3,9,10] ) >>> someSet.data array([ 3, 9, 10], dtype=uint32)
-
invert()[source]¶ Inverts the index set in place.
Example
>>> someSet = uw.container.IndexSet( 15, [1,3,5,7,9,11,13] ) >>> someSet.invert() >>> someSet IndexSet([ 0, 2, 4, 6, 8, 10, 12, 14])
-
remove(indices)[source]¶ Remove item(s) from IndexSet.
Parameters: indices (unsigned int, ndarray, iterable object) – Index or indices to be removed from the IndexSet. Ensure value(s) are integer and non-negative. An iterable object may also be provided, with numpy arrays being significantly more efficient. Note that the ‘remove’ method can not be provided with an IndexSet object, as the ‘add’ object can. Example
>>> someSet = uw.container.IndexSet( 15, [3,9,10] ) >>> someSet IndexSet([ 3, 9, 10]) >>> someSet.remove(3) >>> 3 in someSet False >>> someSet IndexSet([ 9, 10]) >>> someSet.remove([9,10]) >>> someSet IndexSet([])
-
size¶ The size of the IndexSet. Note that the size corresponds to the maximum index value (plus 1) the set is required to hold, NOT the number of elements in the set. So for example, a size of 16, would result in an IndexSet which can retain values between 0 and 15 (inclusive). Note also that the the IndexSet will require ( size/8 + 1 ) bytes of memory storage.
-
class
underworld.container.ObjectifiedIndexSet(object=None, *args, **kwargs)[source]¶ Bases:
underworld.container._indexset.IndexSetThis class simply adds an object to IndexSet data. Usually this object will be the object for which the IndexSet data relates to.. For example, we can attach a Mesh object to an IndexSet containing mesh vertices.
-
__init__(object=None, *args, **kwargs)[source]¶ Class initialiser
Parameters: - object (any, default=None) – Object to tether to data
- parent classes for further parameters. (See) –
Returns: objectifiedIndexSet
Return type:
-
object¶ Object for which IndexSet data relates.
-
underworld.utils module¶
Various utility classes & functions.
underworld.utils.is_kernel |
Function to determine if the script is being run in an ipython or jupyter notebook or in a regular python interpreter. |
underworld.utils.SavedFileData |
A class used to define saved data. |
underworld.utils.Integral |
The Integral class constructs the volume integral |
underworld.utils.MeshVariable_Projection |
This class provides functionality for projecting data from any underworld function onto a provided mesh variable. |
-
class
underworld.utils.SavedFileData(pyobj, filename)[source]¶ Bases:
objectA class used to define saved data.
Parameters: - pyobj (object) – python object saved data relates to.
- filename (str) – filename for saved data, full path
-
class
underworld.utils.Integral(fn, mesh=None, integrationType='volume', surfaceIndexSet=None, integrationSwarm=None, **kwargs)[source]¶ Bases:
underworld._stgermain.StgCompoundComponentThe Integral class constructs the volume integral
\[F_{i} = \int_V \, f_i(\mathbf{x}) \, \mathrm{d} V\]for some function \(f_i\) (specified by a Function object), over some domain \(V\) (specified by an FeMesh object), or the surface integral
\[F_{i} = \oint_{\Gamma} \, f_i(\mathbf{x}) \, \mathrm{d}\Gamma\]for some surface \(\Gamma\) (specified via an IndexSet object on the mesh).
Parameters: - fn (uw.function.Function) – Function to be integrated.
- mesh (uw.mesh.FeMesh) – The mesh over which integration is performed.
- integrationType (str) – Type of integration to perform. Options are “volume” or “surface”.
- surfaceIndexSet (uw.mesh.FeMesh_IndexSet) – Must be provided where integrationType is “surface”. This IndexSet determines which surface is to be integrated over. Note that surface integration over interior nodes is not currently supported.
- integrationSwarm (uw.swarm.IntegrationSwarm (optional)) – User provided integration swarm.
Notes
Constructor must be called by collectively all processes.
Example
Calculate volume of mesh:
>>> import underworld as uw >>> mesh = uw.mesh.FeMesh_Cartesian(minCoord=(0.,0.), maxCoord=(1.,1.)) >>> volumeIntegral = uw.utils.Integral(fn=1.,mesh=mesh) >>> np.allclose( 1., volumeIntegral.evaluate(), rtol=1e-8) True
Calculate surface area of mesh:
>>> surfaceIntegral = uw.utils.Integral(fn=1., mesh=mesh, integrationType='surface', surfaceIndexSet=mesh.specialSets["AllWalls_VertexSet"]) >>> np.allclose( 4., surfaceIntegral.evaluate(), rtol=1e-8) True
-
evaluate()[source]¶ Perform integration.
Notes
Method must be called collectively by all processes.
Returns: result – Integration result. For vector integrals, a vector is returned. Return type: list of floats
-
maskFn¶ The integration mask used where surface integration is performed.
-
class
underworld.utils.MeshVariable_Projection(meshVariable=None, fn=None, voronoi_swarm=None, type=0, **kwargs)[source]¶ Bases:
underworld._stgermain.StgCompoundComponentThis class provides functionality for projecting data from any underworld function onto a provided mesh variable.
For the variable \(\bf{U} = \bf{u}_a N_a\) and function \(F\), we wish to determine appropriate values for \(\bf{u}_a\) such that \(\bf{U} \simeq F\).
Two projection methods are supported; weighted averages and weighted residuals. Generally speaking, weighted averages provide robust low order results, while weighted residuals give higher accuracy but spurious results for difficult functions \(F\).
The weighted average method is defined as:
\[\bf{u}_a = \frac{\int_{\Omega} \bf{F} N_a \partial\Omega }{\int_{\Omega} N_a \partial\Omega }\]The weighted residual method constructs an SLE which is then solved to determine the unknowns:
\[\bf{u}_a\int_{\Omega} N_a N_b \partial\Omega = \int_{\Omega} \bf{F} N_b \partial\Omega\]Parameters: - meshVariable (underworld.mesh.MeshVariable) – The variable you wish to project the function onto.
- fn (underworld.function.Function) – The function you wish to project.
- voronoi_swarm (underworld.swarm.Swarm) – Optional. If a voronoi_swarm is provided, voronoi type integration is utilised to integrate across elements. The provided swarm is used as the basis for the voronoi integration. If no swarm is provided, Gauss integration is used.
- type (int, default=0) – Projection type. 0:weighted average, 1:weighted residual
Notes
Constructor must be called collectively by all processes.
Examples
>>> import underworld as uw >>> import numpy as np >>> mesh = uw.mesh.FeMesh_Cartesian() >>> U = uw.mesh.MeshVariable( mesh, 1 )
Lets cast a constant value onto this mesh variable
>>> const = 1.23456 >>> projector = uw.utils.MeshVariable_Projection( U, const, type=0 ) >>> np.allclose(U.data, const) False >>> projector.solve() >>> np.allclose(U.data, const) True
Now cast mesh coordinates onto a vector variable
>>> U_coord = uw.mesh.MeshVariable( mesh, 2 ) >>> projector = uw.utils.MeshVariable_Projection( U_coord, uw.function.coord(), type=1 ) >>> projector.solve() >>> np.allclose(U_coord.data, mesh.data) True
Project one mesh variable onto another
>>> U_copy = uw.mesh.MeshVariable( mesh, 2 ) >>> projector = uw.utils.MeshVariable_Projection( U_copy, U_coord, type=1 ) >>> projector.solve() >>> np.allclose(U_copy.data, U_coord.data) True
Project the coords to the submesh (usually the constant mesh)
>>> U_submesh = uw.mesh.MeshVariable( mesh.subMesh, 2 ) >>> projector = uw.utils.MeshVariable_Projection( U_submesh, U_coord, type=1 ) >>> projector.solve() >>> np.allclose(U_submesh.data, mesh.subMesh.data) True
Create swarm, then project particle owning elements onto mesh
>>> U_submesh = uw.mesh.MeshVariable( mesh.subMesh, 1 ) >>> swarm = uw.swarm.Swarm(mesh) >>> swarm.populate_using_layout(uw.swarm.layouts.PerCellSpaceFillerLayout(swarm,4)) >>> projector = uw.utils.MeshVariable_Projection( U_submesh, swarm.owningCell, type=1 ) >>> projector.solve() >>> np.allclose(U_submesh.data, mesh.data_elgId) True
underworld.scaling module¶
The scaling module provides units and scaling capabilities.
underworld.scaling.non_dimensionalise |
Non-dimensionalize (scale) provided quantity. |
underworld.scaling.dimensionalise |
Dimensionalise a value. |
underworld.scaling.get_coefficients |
Returns the global scaling dictionary. |
-
underworld.scaling.non_dimensionalise(dimValue)[source]¶ Non-dimensionalize (scale) provided quantity.
This function uses pint to perform a dimension analysis and return a value scaled according to a set of scaling coefficients.
Parameters: dimValue (pint quantity) – Returns: - float (The scaled value.)
- Example
- ——–
- >>> import underworld as uw
- >>> u = uw.scaling.units
- >>> # Characteristic values of the system
- >>> half_rate = 0.5 * u.centimeter / u.year
- >>> model_height = 600e3 * u.meter
- >>> refViscosity = 1e24 * u.pascal * u.second
- >>> surfaceTemp = 0. * u.kelvin
- >>> baseModelTemp = 1330. * u.kelvin
- >>> baseCrustTemp = 550. * u.kelvin
- >>> KL_meters = model_height
- >>> KT_seconds = KL_meters / half_rate
- >>> KM_kilograms = refViscosity * KL_meters * KT_seconds
- >>> Kt_degrees = (baseModelTemp - surfaceTemp)
- >>> K_substance = 1. * u.mole
- >>> scaling_coefficients = uw.scaling.get_coefficients()
- >>> scaling_coefficients[“[time]”] = KT_seconds
- >>> scaling_coefficients[“[length]”] = KL_meters
- >>> scaling_coefficients[“[mass]”] = KM_kilograms
- >>> scaling_coefficients[“[temperature]”] = Kt_degrees
- >>> scaling_coefficients[“[substance]”] -= K_substance
- >>> # Get a scaled value
- >>> gravity = uw.scaling.non_dimensionalise(9.81 * u.meter / u.second**2)
-
underworld.scaling.dimensionalise(value, units)[source]¶ Dimensionalise a value.
Parameters: - value (float, int) – The value to be assigned units.
- units (pint units) – The units to be assigned.
Returns: pint quantity
Return type: dimensionalised value.
Example
>>> import underworld as uw >>> A = uw.scaling.dimensionalise(1.0, u.metre)
underworld.swarm module¶
This module contains routines relating to swarm type objects.
This module contains classes for populating swarms with particles across the domain.
underworld.swarm.layouts.PerCellSpaceFillerLayout |
This layout fills the domain with particles in a quasi-random pattern. |
underworld.swarm.layouts.GlobalSpaceFillerLayout |
This layout fills the domain with particles in a quasi-random pattern. |
underworld.swarm.layouts.PerCellGaussLayout |
This layout populates the domain with particles located at gauss locations within each element of the swarm’s associated finite element mesh. |
underworld.swarm.layouts.PerCellRandomLayout |
This layout fills the domain with particles in a random (per element) pattern. |
-
class
underworld.swarm.layouts.PerCellSpaceFillerLayout(swarm, particlesPerCell, **kwargs)[source]¶ Bases:
underworld.swarm.layouts._PerCellMeshParticleLayoutThis layout fills the domain with particles in a quasi-random pattern. It utilises sobol sequences to generate per element particle locations which are more uniform than that achieved by a purely random generator.
Parameters: - swarm (underworld.swarm.Swarm) – The swarm this layout will act upon
- particlesPerCell (int) – The number of particles per element that this layout will generate.
Example
>>> import underworld as uw >>> mesh = uw.mesh.FeMesh_Cartesian('Q1/dQ0', (1,1), (0.,0.), (1.,1.)) >>> swarm = uw.swarm.Swarm(mesh) >>> layout = uw.swarm.layouts.PerCellSpaceFillerLayout(swarm,particlesPerCell=4) >>> swarm.populate_using_layout(layout) >>> swarm.particleLocalCount 4 >>> swarm.particleCoordinates.data array([[ 0.5 , 0.5 ], [ 0.25 , 0.75 ], [ 0.75 , 0.25 ], [ 0.375, 0.625]])
-
class
underworld.swarm.layouts.GlobalSpaceFillerLayout(swarm, particlesPerCell, **kwargs)[source]¶ Bases:
underworld.swarm.layouts._ParticleLayoutAbstractThis layout fills the domain with particles in a quasi-random pattern. It utilises sobol sequences to generate global particle locations which are more uniform than that achieved by a purely random generator. This layout is mostly useful where populating particles across a rectangular domain.
Parameters: - swarm (underworld.swarm.Swarm) – The swarm this layout will act upon
- particlesPerCell (float) – The average number of particles per element that this layout will generate.
Example
>>> import underworld as uw >>> mesh = uw.mesh.FeMesh_Cartesian('Q1/dQ0', (1,1), (0.,0.), (1.,1.)) >>> swarm = uw.swarm.Swarm(mesh) >>> layout = uw.swarm.layouts.GlobalSpaceFillerLayout(swarm,particlesPerCell=4) >>> swarm.populate_using_layout(layout) >>> swarm.particleLocalCount 4 >>> swarm.particleCoordinates.data array([[ 0.5 , 0.5 ], [ 0.25 , 0.75 ], [ 0.75 , 0.25 ], [ 0.375, 0.625]])
-
class
underworld.swarm.layouts.PerCellGaussLayout(swarm, gaussPointCount, **kwargs)[source]¶ Bases:
underworld.swarm.layouts._ParticleLayoutAbstractThis layout populates the domain with particles located at gauss locations within each element of the swarm’s associated finite element mesh.
Parameters: - swarm (underworld.swarm.Swarm) – The swarm this layout will act upon
- gaussPointCount (int) – Per cell, the number of gauss points in each dimensional direction. Must take an int value between 1 and 5 inclusive.
Example
>>> import underworld as uw >>> # choose mesh to coincide with global element >>> mesh = uw.mesh.FeMesh_Cartesian('Q1/dQ0', (1,1), (-1.,-1.), (1.,1.)) >>> swarm = uw.swarm.Swarm(mesh) >>> layout = uw.swarm.layouts.PerCellGaussLayout(swarm,gaussPointCount=2) >>> swarm.populate_using_layout(layout) >>> swarm.particleLocalCount 4 >>> swarm.particleCoordinates.data array([[-0.57735027, -0.57735027], [ 0.57735027, -0.57735027], [-0.57735027, 0.57735027], [ 0.57735027, 0.57735027]]) >>> import math >>> # lets check one of these gauss points >>> ( swarm.particleCoordinates.data[3][0] - math.sqrt(1./3.) ) < 1.e-10 True
-
class
underworld.swarm.layouts.PerCellRandomLayout(swarm, particlesPerCell, seed=13, **kwargs)[source]¶ Bases:
underworld.swarm.layouts._PerCellMeshParticleLayoutThis layout fills the domain with particles in a random (per element) pattern.
Parameters: - swarm (underworld.swarm.Swarm) – The swarm this layout will act upon
- particlesPerCell (int) – The number of particles per element that this layout will generate.
- seed (int) – Seed for random generator. Default is 13.
Example
>>> import underworld as uw >>> mesh = uw.mesh.FeMesh_Cartesian('Q1/dQ0', (1,1), (0.,0.), (1.,1.)) >>> swarm = uw.swarm.Swarm(mesh) >>> layout = uw.swarm.layouts.PerCellRandomLayout(swarm,particlesPerCell=4) >>> swarm.populate_using_layout(layout) >>> swarm.particleLocalCount 4 >>> swarm.particleCoordinates.data array([[ 0.24261743, 0.67115852], [ 0.16116546, 0.70790335], [ 0.73160516, 0.08792286], [ 0.71953113, 0.15966135]])
underworld.swarm.PopulationControl |
This class implements swarm population control mechanism. |
underworld.swarm.SwarmVariable |
The SwarmVariable class allows users to add data to swarm particles. The data |
underworld.swarm.GaussIntegrationSwarm |
Integration swarm which creates particles within an element at the Gauss points. |
underworld.swarm.VoronoiIntegrationSwarm |
Class for an IntegrationSwarm that maps to another Swarm |
underworld.swarm.Swarm |
The Swarm class supports particle like data structures. |
underworld.swarm.IntegrationSwarm |
Abstract class definition for IntegrationSwarms. |
underworld.swarm.GaussBorderIntegrationSwarm |
Integration swarm which creates particles within the boundary faces of an element, at the Gauss points. |
underworld.swarm.SwarmAbstract |
The SwarmAbstract class supports particle like data structures. |
-
class
underworld.swarm.PopulationControl(swarm, deleteThreshold=0.006, splitThreshold=0.25, maxDeletions=0, maxSplits=3, aggressive=False, aggressiveThreshold=0.8, particlesPerCell=None, **kwargs)[source]¶ Bases:
underworld._stgermain.LeftOverParamsCheckerThis class implements swarm population control mechanism. Population control acts on a per element basic, with a discrete voronoi algorithm is used to determine where particles should be added or removed.
Parameters: - swarm (underworld.swarm.Swarm) – The swarm for which population control should occur.
- deleteThreshold (float) – Particle volume fraction threshold below which particle is deleted. i.e if (particleVolume/elementVolume)<deleteThreshold, then the particle is deleted.
- splitThreshold (float) – Particle volume fraction threshold above which particle is split. i.e if (particleVolume/elementVolume)>splitThreshold, then the particle is split.
- maxDeletions (int) – maximum number of particle deletions per cell
- maxSplits (int) – maximum number of particles splits per cell
- aggressive (bool) – When enabled, this option will invoke aggressive population control in elements where particle counts drop below some threshold.
- aggressiveThreshold (float) – lower cell particle population threshold beyond which aggressive population control occurs. i.e if (cellParticleCount/particlesPerCell)<aggressiveThreshold, then aggressive pop control will occur. Note that this option is only valid if ‘aggressive’ is enabled.
- particlesPerCell (int) – This is the desired number of particles each element should contain. Note that this option is only valid if ‘aggressive’ is enabled.
Example
This simple example generates a swarm, then applies population control to split particles.
>>> import underworld as uw >>> import numpy as np >>> mesh = uw.mesh.FeMesh_Cartesian() >>> swarm = uw.swarm.Swarm(mesh) >>> swarm.populate_using_layout(uw.swarm.layouts.PerCellGaussLayout(swarm,4)) >>> population_control = uw.swarm.PopulationControl(swarm,deleteThreshold=0.,splitThreshold=0.,maxDeletions=0,maxSplits=9999) >>> population_control.repopulate() >>> swarm.particleGlobalCount 512
-
class
underworld.swarm.SwarmVariable(swarm, dataType, count, writeable=True, **kwargs)[source]¶ Bases:
underworld._stgermain.StgClass,underworld.function._function.Function- The SwarmVariable class allows users to add data to swarm particles. The data
- can be of type “char”, “short”, “int”, “long, “float” or “double”.
Note that the swarm allocates one block of contiguous memory for all the particles. The per particle variable datums is then interlaced across this memory block.
The recommended practise is to add all swarm variables before populating the swarm to avoid costly reallocations.
Swarm variables should be added via the add_variable swarm method.
Parameters: - swarm (underworld.swarm.Swarm) – The swarm of particles for which we wish to add the variable
- dataType (str) – The data type for the variable. Available types are “char”, “short”, “int”, “long”, “float” or “double”.
- count (unsigned) – The number of values to be stored for each particle.
- writeable (bool) – Signifies if the variable should be writeable.
-
count¶ Returns: Number of data items for this variable stored on each particle. Return type: int
-
data¶ Returns: Numpy proxy array to underlying variable data. Note that the returned array is a proxy for all the local particle data. As numpy arrays are simply proxys to the underlying memory structures, no data copying is required. Return type: numpy.ndarray Example
>>> # create mesh >>> mesh = uw.mesh.FeMesh_Cartesian( elementType='Q1/dQ0', elementRes=(16,16), minCoord=(0.,0.), maxCoord=(1.,1.) ) >>> # create empty swarm >>> swarm = uw.swarm.Swarm(mesh) >>> # add a variable >>> svar = swarm.add_variable("int",1) >>> # add particles >>> swarm.populate_using_layout(uw.swarm.layouts.PerCellGaussLayout(swarm,2)) >>> swarm.particleLocalCount 1024 >>> len(svar.data) # should be the same as particle local count 1024 >>> swarm.owningCell.data # check particle owning cells/elements. array([[ 0], [ 0], [ 0], ..., [255], [255], [255]], dtype=int32)
>>> # particle coords >>> swarm.particleCoordinates.data[0] array([ 0.0132078, 0.0132078]) >>> # move the particle >>> with swarm.deform_swarm(): ... swarm.particleCoordinates.data[0] = [0.2,0.2] >>> swarm.particleCoordinates.data[0] array([ 0.2, 0.2])
-
dataType¶ Returns: Data type for variable. Supported types are ‘char’, ‘short’, ‘int’, ‘long’, ‘float’ and ‘double’. Return type: str
-
data_shadow¶ Returns: Numpy proxy array to underlying variable shadow data. Return type: numpy.ndarray Example
Refer to example provided for ‘data’ property(/method).
-
load(filename)[source]¶ Load the swarm variable from disk. This must be called after the swarm.load().
Parameters: filename (str) – The filename for the saved file. Relative or absolute paths may be used, but all directories must exist. Notes
This method must be called collectively by all processes.
Example
Refer to example provided for ‘save’ method.
-
save(filename)[source]¶ Save the swarm variable to disk.
Parameters: - filename (str) – The filename for the saved file. Relative or absolute paths may be used, but all directories must exist.
- swarmHandle (uw.utils.SavedFileData , optional) – The saved swarm file handle. If provided, a reference to the swarm file is made. Currently this doesn’t provide any extra functionality.
Returns: Data object relating to saved file. This only needs to be retained if you wish to create XDMF files and can be ignored otherwise.
Return type: Notes
This method must be called collectively by all processes.
Example
First create the swarm, populate, then add a variable:
>>> mesh = uw.mesh.FeMesh_Cartesian( elementType='Q1/dQ0', elementRes=(16,16), minCoord=(0.,0.), maxCoord=(1.,1.) ) >>> swarm = uw.swarm.Swarm(mesh) >>> swarm.populate_using_layout(uw.swarm.layouts.PerCellGaussLayout(swarm,2)) >>> svar = swarm.add_variable("int",1)
Write something to variable
>>> import numpy as np >>> svar.data[:,0] = np.arange(swarm.particleLocalCount)
Save to a file:
>>> ignoreMe = swarm.save("saved_swarm.h5") >>> ignoreMe = svar.save("saved_swarm_variable.h5")
Now let’s try and reload. First create a new swarm and swarm variable, and then load both:
>>> clone_swarm = uw.swarm.Swarm(mesh) >>> clone_svar = clone_swarm.add_variable("int",1) >>> clone_swarm.load("saved_swarm.h5") >>> clone_svar.load("saved_swarm_variable.h5")
Now check for equality:
>>> import numpy as np >>> np.allclose(svar.data,clone_svar.data) True
>>> # clean up: >>> if uw.rank() == 0: ... import os; ... os.remove( "saved_swarm.h5" ) ... os.remove( "saved_swarm_variable.h5" )
-
swarm¶ Returns: The swarm this variable belongs to. Return type: underworld.swarm.Swarm
-
xdmf(filename, varSavedData, varname, swarmSavedData, swarmname, modeltime=0.0)[source]¶ Creates an xdmf file, filename, associating the varSavedData file on the swarmSavedData file
Notes
xdmf contain 2 files: an .xml and a .h5 file. See http://www.xdmf.org/index.php/Main_Page This method only needs to be called by the master process, all other processes return quietly.
Parameters: - filename (str) – The output path to write the xdmf file. Relative or absolute paths may be used, but all directories must exist.
- varname (str) – The xdmf name to give the swarmVariable
- swarmname (str) – The xdmf name to give the swarm
- swarmSavedData (underworld.utils.SaveFileData) – Handler returned for saving a swarm. underworld.swarm.Swarm.save(xxx)
- varSavedData (underworld.utils.SavedFileData) – Handler returned from saving a SwarmVariable. underworld.swarm.SwarmVariable.save(xxx)
- modeltime (float (default 0.0)) – The time recorded in the xdmf output file
Example
First create the swarm and add a variable:
>>> mesh = uw.mesh.FeMesh_Cartesian( elementType='Q1/dQ0', elementRes=(16,16), minCoord=(0.,0.), maxCoord=(1.,1.) ) >>> swarm = uw.swarm.Swarm( mesh=mesh ) >>> swarmLayout = uw.swarm.layouts.PerCellGaussLayout(swarm,2) >>> swarm.populate_using_layout( layout=swarmLayout ) >>> swarmVar = swarm.add_variable( dataType="int", count=1 )
Write something to variable
>>> import numpy as np >>> swarmVar.data[:,0] = np.arange(swarmVar.data.shape[0])
Save mesh and var to a file:
>>> swarmDat = swarm.save("saved_swarm.h5") >>> swarmVarDat = swarmVar.save("saved_swarmvariable.h5")
Now let’s create the xdmf file
>>> swarmVar.xdmf("TESTxdmf", swarmVarDat, "var1", swarmDat, "MrSwarm" )
Does file exist?
>>> import os >>> if uw.rank() == 0: os.path.isfile("TESTxdmf.xdmf") True
>>> # clean up: >>> if uw.rank() == 0: ... import os; ... os.remove( "saved_swarm.h5" ) ... os.remove( "saved_swarmvariable.h5" ) ... os.remove( "TESTxdmf.xdmf" )
-
class
underworld.swarm.GaussIntegrationSwarm(mesh, particleCount=None, **kwargs)[source]¶ Bases:
underworld.swarm._integration_swarm.IntegrationSwarmIntegration swarm which creates particles within an element at the Gauss points.
Parameters: - mesh (underworld.mesh.FeMesh) – The FeMesh the swarm is supported by. See Swarm.mesh property docstring for further information.
- particleCount (unsigned. Default is 3, unless Q2 mesh which takes default 5.) – Number of gauss particles in each direction. Must take value in [1,5].
-
class
underworld.swarm.VoronoiIntegrationSwarm(swarm, **kwargs)[source]¶ Bases:
underworld.swarm._integration_swarm.IntegrationSwarm,underworld.function._function.FunctionInputClass for an IntegrationSwarm that maps to another Swarm
Note that this swarm can act as a function input. In this capacity, the fundamental function input type is the FEMCoordinate (ie, the particle local coordinate, the owning mesh, and the owning element). This input can be reduced to the global coordinate when returned within python. The FEMCoordinate particle representation is useful when deforming a mesh, as it is possible to deform the mesh, and then use the FEMCoordinate to reset the particles within the moved mesh.
Parameters: swarm (underworld.swarm.Swarm) – The PIC integration swarm maps to this user provided swarm. Example
This simple example checks that the true global coordiante, and that derived from the local coordinate, are close to equal. Note that the VoronoiIntegrationSwarm uses a voronoi centroid algorithm so we do not expect particle to exactly coincide.
>>> import underworld as uw >>> import numpy as np >>> mesh = uw.mesh.FeMesh_Cartesian() >>> swarm = uw.swarm.Swarm(mesh) >>> swarm.populate_using_layout(uw.swarm.layouts.PerCellGaussLayout(swarm,4)) >>> vswarm = uw.swarm.VoronoiIntegrationSwarm(swarm) >>> vswarm.repopulate() >>> np.allclose(swarm.particleCoordinates.data, uw.function.input().evaluate(vswarm),atol=1e-1) True
-
repopulate(weights_calculator=None)[source]¶ This method repopulates the voronoi swarm using the provided global swarm. The weights are also recalculated.
Parameters: weights_calculator (underworld.swarm.Weights) – The weights calculator for the Voronoi swarm. If none is provided, a default DVCWeights calculator is used.
-
-
class
underworld.swarm.Swarm(mesh, particleEscape=False, **kwargs)[source]¶ Bases:
underworld.swarm._swarmabstract.SwarmAbstract,underworld.function._function.FunctionInput,underworld._stgermain.SaveThe Swarm class supports particle like data structures. Each instance of this class will store a set of unique particles. In this context, particles are data structures which store a location variable, along with any other variables the user requests.
Parameters: - mesh (underworld.mesh.FeMesh) – The FeMesh the swarm is supported by. See Swarm.mesh property docstring for further information.
- particleEscape (bool) – If set to true, particles are deleted when they leave the domain. This may occur during particle advection, or when the mesh is deformed.
Example
Create a swarm with some variables:
>>> # First we need a mesh: >>> mesh = uw.mesh.FeMesh_Cartesian( elementType='Q1/dQ0', elementRes=(16,16), minCoord=(0.,0.), maxCoord=(1.,1.) ) >>> # Create empty swarm: >>> swarm = uw.swarm.Swarm(mesh) >>> # Add a variable: >>> svar = swarm.add_variable("char",1) >>> # Add another: >>> svar = swarm.add_variable("double",3) >>> # Now use a layout to fill with particles >>> swarm.particleLocalCount 0 >>> layout = uw.swarm.layouts.PerCellGaussLayout(swarm,2) >>> swarm.populate_using_layout(layout) >>> swarm.particleLocalCount 1024 >>> swarm.data[0] array([ 0.0132078, 0.0132078]) >>> swarm.owningCell.data[0] array([0], dtype=int32)
With particleEscape enabled, particles which are no longer within the mesh domain are deleted.
>>> mesh = uw.mesh.FeMesh_Cartesian( elementType='Q1/dQ0', elementRes=(16,16), minCoord=(0.,0.), maxCoord=(1.,1.) ) >>> swarm = uw.swarm.Swarm(mesh, particleEscape=True) >>> swarm.particleLocalCount 0 >>> layout = uw.swarm.layouts.PerCellGaussLayout(swarm,2) >>> swarm.populate_using_layout(layout) >>> swarm.particleGlobalCount 1024 >>> with mesh.deform_mesh(): ... mesh.data[:] += (0.5,0.) >>> swarm.particleGlobalCount 512
Alternatively, moving the particles:
>>> mesh = uw.mesh.FeMesh_Cartesian( elementType='Q1/dQ0', elementRes=(16,16), minCoord=(0.,0.), maxCoord=(1.,1.) ) >>> swarm = uw.swarm.Swarm(mesh, particleEscape=True) >>> swarm.particleLocalCount 0 >>> layout = uw.swarm.layouts.PerCellGaussLayout(swarm,2) >>> swarm.populate_using_layout(layout) >>> swarm.particleGlobalCount 1024 >>> with swarm.deform_swarm(): ... swarm.data[:] -= (0.5,0.) >>> swarm.particleGlobalCount 512
-
add_particles_with_coordinates(coordinatesArray)[source]¶ This method adds particles to the swarm using particle coordinates provided using a numpy array.
Note that particles with coordinates NOT local to the current processor will be reject/ignored.
Parameters: coordinatesArray (numpy.ndarray) – The numpy array containing the coordinate of the new particles. Array is expected to take shape n*dim, where n is the number of new particles, and dim is the dimensionality of the swarm’s supporting mesh. Returns: Array containing the local index of the added particles. Rejected particles are denoted with an index of -1. Return type: numpy.ndarray Example
>>> mesh = uw.mesh.FeMesh_Cartesian( elementType='Q1/dQ0', elementRes=(4,4), minCoord=(0.,0.), maxCoord=(1.,1.) ) >>> swarm = uw.swarm.Swarm(mesh) >>> import numpy as np >>> arr = np.zeros((5,2)) >>> arr[0] = [0.1,0.1] >>> arr[1] = [0.2,0.1] >>> arr[2] = [0.1,0.2] >>> arr[3] = [-0.1,-0.1] >>> arr[4] = [0.8,0.8] >>> swarm.add_particles_with_coordinates(arr) array([ 0, 1, 2, -1, 3], dtype=int32) >>> swarm.particleLocalCount 4 >>> swarm.data array([[ 0.1, 0.1], [ 0.2, 0.1], [ 0.1, 0.2], [ 0.8, 0.8]])
-
deform_swarm(**kwds)[source]¶ Any particle location modifications must occur within this python context manager. This is necessary as it is critical that certain internal objects are updated when particle locations are modified.
Parameters: update_owners (bool, default=True) – If this is set to False, particle ownership (which element owns a particular particle) is not updated at the conclusion of the context manager. This is often necessary when both the mesh and particles are advecting simutaneously. Example
>>> mesh = uw.mesh.FeMesh_Cartesian( elementType='Q1/dQ0', elementRes=(16,16), minCoord=(0.,0.), maxCoord=(1.,1.) ) >>> swarm = uw.swarm.Swarm(mesh) >>> layout = uw.swarm.layouts.PerCellGaussLayout(swarm,2) >>> swarm.populate_using_layout(layout) >>> swarm.data[0] array([ 0.0132078, 0.0132078])
Attempted modification without using deform_swarm() should fail:
>>> swarm.data[0] = [0.2,0.2] Traceback (most recent call last): ... ValueError: assignment destination is read-only
Within the deform_swarm() context manager, modification is allowed:
>>> with swarm.deform_swarm(): ... swarm.data[0] = [0.2,0.2] >>> swarm.data[0] array([ 0.2, 0.2])
-
fn_particle_found()[source]¶ This function returns True where a particle is able to be found using the provided input to function evaluation.
Returns: The function object. Return type: underworld.function.Function Example
Setup some things:
>>> import numpy as np >>> mesh = uw.mesh.FeMesh_Cartesian(elementRes=(32,32)) >>> swarm = uw.swarm.Swarm(mesh, particleEscape=True) >>> layout = uw.swarm.layouts.PerCellGaussLayout(swarm,2) >>> swarm.populate_using_layout(layout)
Now, evaluate the fn_particle_found function on the swarm.. all should be true
>>> fn_pf = swarm.fn_particle_found() >>> fn_pf.evaluate(swarm).all() True
Evalute at arbitrary coord… should return False
>>> fn_pf.evaluate( (0.3,0.9) ) array([[False]], dtype=bool)
Now, lets get rid of all particles outside of a circle, and look to obtain pi/4. First eject particles:
>>> with swarm.deform_swarm(): ... for ind,coord in enumerate(swarm.data): ... if np.dot(coord,coord)>1.: ... swarm.data[ind] = (99999.,99999.)
Now integrate and test
>>> incirc = uw.function.branching.conditional( ( (fn_pf,1.),(True,0.) ) ) >>> np.isclose(uw.utils.Integral(incirc, mesh).evaluate(),np.pi/4.,rtol=2e-2) array([ True], dtype=bool)
-
load(filename, try_optimise=True, verbose=False)[source]¶ Load a swarm from disk. Note that this must be called before any SwarmVariable members are loaded.
Parameters: - filename (str) – The filename for the saved file. Relative or absolute paths may be used.
- try_optimise (bool, Default=True) – Will speed up the swarm load time but warning - this algorithm assumes the previously saved swarm data was made on an identical mesh and mesh partitioning (number of processors) with respect to the current mesh. If this isn’t the case then the reloaded particle ordering will be broken, leading to an invalid swarms. One can disable this optimisation and revert to a brute force algorithm, much slower, by setting this option to False.
- verbose (bool) – Prints a swarm load progress bar.
Notes
This method must be called collectively by all processes.
Example
Refer to example provided for ‘save’ method.
-
particleGlobalCount¶ Returns: The global number (across all processes) of particles in the swarm. Return type: int
-
save(filename)[source]¶ Save the swarm to disk.
Parameters: filename (str) – The filename for the saved file. Relative or absolute paths may be used, but all directories must exist. Returns: Data object relating to saved file. This only needs to be retained if you wish to create XDMF files and can be ignored otherwise. Return type: underworld.utils.SavedFileData Notes
This method must be called collectively by all processes.
Example
First create the swarm, and populate with layout:
>>> mesh = uw.mesh.FeMesh_Cartesian( elementType='Q1/dQ0', elementRes=(16,16), minCoord=(0.,0.), maxCoord=(1.,1.) ) >>> swarm = uw.swarm.Swarm(mesh) >>> swarm.populate_using_layout(uw.swarm.layouts.PerCellGaussLayout(swarm,2))
Save to a file:
>>> ignoreMe = swarm.save("saved_swarm.h5")
Now let’s try and reload. First create an empty swarm, and then load:
>>> clone_swarm = uw.swarm.Swarm(mesh) >>> clone_swarm.load( "saved_swarm.h5" )
Now check for equality:
>>> import numpy as np >>> np.allclose(swarm.data,clone_swarm.data) True
>>> # clean up: >>> if uw.rank() == 0: ... import os; ... os.remove( "saved_swarm.h5" )
-
shadow_particles_fetch()[source]¶ When called, neighbouring processor particles which have coordinates within the current processor’s shadow zone will be communicated to the current processor. Ie, the processor shadow zone is populated using particles that are owned by neighbouring processors. After this method has been called, particle shadow data is available via the data_shadow handles of SwarmVariable objects. This data is read only.
Note that you will need to call this whenever neighbouring information has potentially changed, for example after swarm advection, or after you have modified a SwarmVariable object.
Any existing shadow information will be discarded when this is called.
Notes
This method must be called collectively by all processes.
-
update_particle_owners()[source]¶ This routine will update particles owners after particles have been moved. This is both in terms of the cell/element the the particle resides within, and also in terms of the parallel processor decomposition (particles belonging on other processors will be sent across).
Users should not generally need to call this as it will be called automatically at the conclusion of a deform_swarm() block.
Notes
This method must be called collectively by all processes.
Example
>>> mesh = uw.mesh.FeMesh_Cartesian( elementType='Q1/dQ0', elementRes=(16,16), minCoord=(0.,0.), maxCoord=(1.,1.) ) >>> swarm = uw.swarm.Swarm(mesh) >>> swarm.populate_using_layout(uw.swarm.layouts.PerCellGaussLayout(swarm,2)) >>> swarm.data[0] array([ 0.0132078, 0.0132078]) >>> swarm.owningCell.data[0] array([0], dtype=int32) >>> with swarm.deform_swarm(): ... swarm.data[0] = [0.1,0.1] >>> swarm.owningCell.data[0] array([17], dtype=int32)
-
class
underworld.swarm.IntegrationSwarm(mesh, **kwargs)[source]¶ Bases:
underworld.swarm._swarmabstract.SwarmAbstractAbstract class definition for IntegrationSwarms.
All IntegrationSwarms have the following SwarmVariables from this class:
- localCoordVariable : double (number of particle, dimensions)
- For local element coordinates of the particle
- weightVariable : double (number of particles)
- For the integration weight of each particle
-
particleWeights¶ Returns: Swarm variable recording the weight of the swarm particles. Return type: underworld.swarm.SwarmVariable
-
class
underworld.swarm.GaussBorderIntegrationSwarm(mesh, particleCount=None, **kwargs)[source]¶ Bases:
underworld.swarm._integration_swarm.GaussIntegrationSwarmIntegration swarm which creates particles within the boundary faces of an element, at the Gauss points.
See parent class for parameters.
-
class
underworld.swarm.SwarmAbstract(mesh, **kwargs)[source]¶ Bases:
underworld._stgermain.StgCompoundComponentThe SwarmAbstract class supports particle like data structures. Each instance of this class will store a set of unique particles. In this context, particles are data structures which store a location variable, along with any other variables the user requests.
Parameters: mesh (underworld.mesh.FeMesh) – The FeMesh the swarm is supported by. See Swarm.mesh property docstring for further information. -
add_variable(dataType, count)[source]¶ Add a variable to each particle in this swarm. Variables can be added at any point. Removal of variables is however not currently supported. See help(SwarmVariable) for further information.
Parameters: - dataType (str) – The data type for the variable. Available types are “char”, “short”, “int”, “float” or “double”.
- count (unsigned) – The number of values to be stored for each particle.
Returns: The newly created swarm variable.
Return type: Example
>>> # first we need a mesh >>> mesh = uw.mesh.FeMesh_Cartesian( elementType='Q1/dQ0', elementRes=(16,16), minCoord=(0.,0.), maxCoord=(1.,1.) ) >>> # create swarm >>> swarm = uw.swarm.Swarm(mesh) >>> # add a variable >>> svar = swarm.add_variable("char",1) >>> # add another >>> svar = swarm.add_variable("double",3) >>> # add some particles >>> swarm.populate_using_layout(uw.swarm.layouts.PerCellGaussLayout(swarm,2)) >>> # add another variable >>> svar = swarm.add_variable("double",5)
-
data¶ Returns: Numpy proxy array to underlying particle coordinate data. Note that this is an alias to swarm.particleCoordinates.data. Check SwarmVariable.data for further info. Return type: numpy.ndarray
-
globalId¶ Returns: Swarm variable recording a particle global identifier. Not yet implemented. Return type: underworld.swarm.SwarmVariable
-
mesh¶ Returns: Supporting FeMesh for this Swarm. All swarms are required to be supported by mesh (or similar) objects, which provide the data structures necessary for efficient particle locating/tracking, as well as the necessary mechanism for the swarm parallel decomposition. Return type: underworld.mesh.FeMesh
-
owningCell¶ Returns: Swarm variable recording the owning cell of the swarm particles. This will usually correspond to the owning element local id. Return type: underworld.swarm.SwarmVariable
-
particleCoordinates¶ Returns: Swarm variable recording the coordinates of the swarm particles. Return type: underworld.swarm.SwarmVariable
-
particleLocalCount¶ Returns: Number of swarm particles within the current processes local space. Return type: int
-
populate_using_layout(layout)[source]¶ This method uses the provided layout to populate the swarm with particles. Usually layouts add particles across the entire domain. Available layouts may be found in the swarm.layouts module. Note that Layouts can only currently be used on empty swarms. Also note that all numpy arrays associated with swarm variables must be deleted before a layout can be applied.
Parameters: layout (underworld.swarm.layouts._ParticleLayoutAbstract) – The layout which determines where particles are created and added. Example
>>> # first we need a mesh >>> mesh = uw.mesh.FeMesh_Cartesian( elementType='Q1/dQ0', elementRes=(16,16), minCoord=(0.,0.), maxCoord=(1.,1.) ) >>> # create swarm >>> swarm = uw.swarm.Swarm(mesh) >>> # add populate >>> swarm.populate_using_layout(uw.swarm.layouts.PerCellGaussLayout(swarm,2))
-
stateId¶ Returns: Swarm state identifier. This is incremented whenever the swarm is modified. Return type: int
-
variables¶ Returns: List of swarm variables associated with this swarm. Return type: list
-
underworld.mesh module¶
Implementation relating to meshing.
underworld.mesh.FeMesh_IndexSet |
This class ties the FeMesh instance to an index set, and stores other metadata relevant to the set. |
underworld.mesh.FeMesh |
The FeMesh class provides the geometry and topology of a finite element discretised domain. |
underworld.mesh.FeMesh_Cartesian |
This class generates a finite element mesh which is topologically cartesian and geometrically regular. |
underworld.mesh.MeshVariable |
The MeshVariable class generates a variable supported by a finite element mesh. |
-
class
underworld.mesh.FeMesh_IndexSet(object, topologicalIndex=None, *args, **kwargs)[source]¶ Bases:
underworld.container._indexset.ObjectifiedIndexSet,underworld.function._function.FunctionInputThis class ties the FeMesh instance to an index set, and stores other metadata relevant to the set.
Parameters: - object (underworld.mesh.FeMesh) – The FeMesh instance from which the IndexSet was extracted.
- topologicalIndex (int) – Mesh topological index for which the IndexSet relates. See docstring for further info.
Example
>>> amesh = uw.mesh.FeMesh_Cartesian( elementType='Q1/dQ0', elementRes=(4,4), minCoord=(0.,0.), maxCoord=(1.,1.) ) >>> iset = uw.libUnderworld.StgDomain.RegularMeshUtils_CreateGlobalMaxISet( amesh._mesh ) >>> uw.mesh.FeMesh_IndexSet( amesh, topologicalIndex=0, size=amesh.nodesGlobal, fromObject=iset ) FeMesh_IndexSet([ 4, 9, 14, 19, 24])
-
topologicalIndex¶ Returns: - The topological index for the indices. The mapping is:
- 0 - vertex 1 - edge 2 - face 3 - volume
Return type: int
-
class
underworld.mesh.FeMesh(elementType, generator=None, **kwargs)[source]¶ Bases:
underworld._stgermain.StgCompoundComponent,underworld.function._function.FunctionInputThe FeMesh class provides the geometry and topology of a finite element discretised domain. The FeMesh is implicitly parallel. Some aspects may be local or global, but this is generally handled automatically.
A number of element types are supported.
Parameters: - elementType (str) – Element type for FeMesh. See FeMesh.elementType docstring for further info.
- generator (underworld.mesh.MeshGenerator) – Generator object which builds the FeMesh. See FeMesh.generator docstring for further info.
-
add_post_deform_function(function)[source]¶ Adds a function function to be executed after mesh deformation is applied.
Parameters: function (FunctionType) – Python (not underworld) function to be executed. Closures should be used where parameters are required.
-
add_pre_deform_function(function)[source]¶ Adds a function function to be executed before mesh deformation is applied.
Parameters: function (FunctionType) – Python (not underworld) function to be executed. Closures should be used where parameters are required.
-
add_variable(nodeDofCount, dataType='double', **kwargs)[source]¶ Creates and returns a mesh variable using the discretisation of the given mesh.
To set / read nodal values, use the numpy interface via the ‘data’ property.
Parameters: - dataType (string) – The data type for the variable. Note that only ‘double’ type variables are currently supported.
- nodeDofCount (int) – Number of degrees of freedom per node the variable will have
Returns: The newly created mesh variable.
Return type: Example
>>> linearMesh = uw.mesh.FeMesh_Cartesian( elementType='Q1/dQ0', elementRes=(16,16), minCoord=(0.,0.), maxCoord=(1.,1.) ) >>> scalarFeVar = linearMesh.add_variable( nodeDofCount=1, dataType="double" ) >>> q0field = linearMesh.subMesh.add_variable( 1 ) # adds variable to secondary elementType discretisation
-
data¶ Numpy proxy array proxy to underlying object vertex data. Note that the returned array is a proxy for all the local vertices, and it is provided as 1d list.
As these arrays are simply proxys to the underlying memory structures, no data copying is required.
Note that this property returns a read-only numpy array as default. If you wish to modify mesh vertex locations, you are required to use the deform_mesh context manager.
If you are modifying the mesh, remember to modify any submesh associated with it accordingly.
Returns: The data proxy array. Return type: numpy.ndarray Example
>>> import underworld as uw >>> someMesh = uw.mesh.FeMesh_Cartesian( elementType='Q1', elementRes=(2,2), minCoord=(-1.,-1.), maxCoord=(1.,1.) ) >>> someMesh.data.shape (9, 2)
You can retrieve individual vertex locations
>>> someMesh.data[1] array([ 0., -1.])
You can modify these locations directly, but take care not to tangle the mesh! Mesh modifications must occur within the deform_mesh context manager.
>>> with someMesh.deform_mesh(): ... someMesh.data[1] = [0.1,-1.1] >>> someMesh.data[1] array([ 0.1, -1.1])
-
data_elementNodes¶ Returns: Array specifying the nodes (global node id) for a given element (local element id). NOTE: Length is local size. Return type: numpy.ndarray
-
data_elgId¶ Returns: Array specifying global element ids. Length is domain size, (local+shadow). Return type: numpy.ndarray
-
data_nodegId¶ Returns: Array specifying global node ids. Length is domain size, (local+shadow). Return type: numpy.ndarray
-
deform_mesh(**kwds)[source]¶ Any mesh deformation must occur within this python context manager. Note that certain algorithms may be switched to their irregular mesh equivalents (if not already set this way). This may have performance implications.
Any submesh will also be appropriately updated on return from the context manager, as will various mesh metrics.
Parameters: isRegular (bool) – The general assumption is that the deformed mesh will no longer be regular (orthonormal), and more general but less efficient algorithms will be selected via this context manager. To over-ride this behaviour, set this parameter to True. Example
>>> import underworld as uw >>> someMesh = uw.mesh.FeMesh_Cartesian() >>> with someMesh.deform_mesh(): ... someMesh.data[0] = [0.1,0.1] >>> someMesh.data[0] array([ 0.1, 0.1])
-
elementType¶ Returns: Element type for FeMesh. Supported types are “Q2”, “Q1”, “dQ1”, “dPc1” and “dQ0”. Return type: str
-
elementsDomain¶ Returns: Returns the number of domain (local+shadow) elements on the mesh Return type: int Example
>>> someMesh = uw.mesh.FeMesh_Cartesian( elementType='Q1', elementRes=(2,2), minCoord=(-1.,-1.), maxCoord=(1.,1.) ) >>> someMesh.elementsDomain 4
-
elementsGlobal¶ Returns: Returns the number of global elements on the mesh Return type: int Example
>>> someMesh = uw.mesh.FeMesh_Cartesian( elementType='Q1', elementRes=(2,2), minCoord=(-1.,-1.), maxCoord=(1.,1.) ) >>> someMesh.elementsGlobal 4
-
elementsLocal¶ Returns: Returns the number of local elements on the mesh Return type: int Example
>>> someMesh = uw.mesh.FeMesh_Cartesian( elementType='Q1', elementRes=(2,2), minCoord=(-1.,-1.), maxCoord=(1.,1.) ) >>> someMesh.elementsLocal 4
-
generator¶ Getter/Setter for the mesh MeshGenerator object.
Returns: Object which builds the mesh. Note that the mesh itself may be a generator, in which case this property will return the mesh object iself. Return type: underworld.mesh.MeshGenerator
-
load(filename)[source]¶ Load the mesh from disk.
Parameters: filename (str) – The filename for the saved file. Relative or absolute paths may be used, but all directories must exist. Notes
This method must be called collectively by all processes.
If the file data array is the same length as the current mesh global size, it is assumed the file contains compatible data. Note that this may not be the case where for example where you have saved using a 2*8 resolution mesh, then loaded using an 8*2 resolution mesh.
Provided files must be in hdf5 format, and use the correct schema.
Example
Refer to example provided for ‘save’ method.
-
nodesDomain¶ Returns: Returns the number of domain (local+shadow) nodes on the mesh. Return type: int
-
nodesGlobal¶ Returns: Returns the number of global nodes on the mesh Return type: int Example
>>> someMesh = uw.mesh.FeMesh_Cartesian( elementType='Q1', elementRes=(2,2), minCoord=(-1.,-1.), maxCoord=(1.,1.) ) >>> someMesh.nodesGlobal 9
-
nodesLocal¶ Returns: Returns the number of local nodes on the mesh. Return type: int
-
reset()[source]¶ Reset the mesh.
Templated mesh (such as the DQ0 mesh) will be reset according to the current state of their geometryMesh.
Other mesh (such as the Q1 & Q2) will be reset to their post-construction state.
Notes
This method must be called collectively by all processes.
-
save(filename)[source]¶ Save the mesh to disk
Parameters: filename (string) – The name of the output file. Returns: Data object relating to saved file. This only needs to be retained if you wish to create XDMF files and can be ignored otherwise. Return type: underworld.utils.SavedFileData Notes
This method must be called collectively by all processes.
Example
First create the mesh:
>>> mesh = uw.mesh.FeMesh_Cartesian( elementType='Q1/dQ0', elementRes=(16,16), minCoord=(0.,0.), maxCoord=(1.,1.) )
Save to a file (note that the ‘ignoreMe’ object isn’t really required):
>>> ignoreMe = mesh.save("saved_mesh.h5")
Now let’s try and reload. First create new mesh (note the different spatial size):
>>> clone_mesh = uw.mesh.FeMesh_Cartesian( elementType='Q1/dQ0', elementRes=(16,16), minCoord=(0.,0.), maxCoord=(1.5,1.5) )
Confirm clone mesh is different from original mesh:
>>> import numpy as np >>> np.allclose(mesh.data,clone_mesh.data) False
Now reload using saved file:
>>> clone_mesh.load("saved_mesh.h5")
Now check for equality:
>>> np.allclose(mesh.data,clone_mesh.data) True
>>> # clean up: >>> if uw.rank() == 0: ... import os; ... os.remove( "saved_mesh.h5" )
-
specialSets¶ Returns: This dictionary stores a set of special data sets relating to mesh objects. Return type: dict Example
>>> import underworld as uw >>> someMesh = uw.mesh.FeMesh_Cartesian( elementType='Q1', elementRes=(2,2), minCoord=(0.,0.), maxCoord=(1.,1.) ) >>> someMesh.specialSets.keys() ['MaxI_VertexSet', 'Top_VertexSet', 'Left_VertexSet', 'MinI_VertexSet', 'AllWalls_VertexSet', 'Bottom_VertexSet', 'Right_VertexSet', 'MinJ_VertexSet', 'MaxJ_VertexSet', 'Empty'] >>> someMesh.specialSets["MinJ_VertexSet"] FeMesh_IndexSet([0, 1, 2])
-
class
underworld.mesh.FeMesh_Cartesian(elementType='Q1/dQ0', elementRes=(4, 4), minCoord=(0.0, 0.0), maxCoord=(1.0, 1.0), periodic=None, partitioned=True, **kwargs)[source]¶ Bases:
underworld.mesh._mesh.FeMesh,underworld.mesh._mesh.CartesianMeshGeneratorThis class generates a finite element mesh which is topologically cartesian and geometrically regular. It is possible to directly build a dual mesh by passing a pair of element types to the constructor.
Refer to parent classes for parameters beyond those below.
Parameters: - elementType (str) – Mesh element type. Note that this class allows the user to (optionally) provide a pair of elementTypes for which a dual mesh will be created. The submesh is accessible through the ‘subMesh’ property. The primary mesh itself is the object returned by this constructor.
- elementRes (list,tuple) – List or tuple of ints specifying mesh resolution. See CartesianMeshGenerator.elementRes docstring for further information.
- minCoord (list, tuple) – List or tuple of floats specifying minimum mesh location. See CartesianMeshGenerator.minCoord docstring for further information.
- maxCoord (list, tuple) – List or tuple of floats specifying maximum mesh location. See CartesianMeshGenerator.maxCoord docstring for further information.
- periodic (list, tuple) – List or tuple of bools, specifying mesh periodicity in each direction.
- partitioned (bool) – If false, the mesh is not partitioned across entire processor pool. Instead mesh is entirely owned by processor which generated it.
Examples
To create a linear mesh:
>>> import underworld as uw >>> someMesh = uw.mesh.FeMesh_Cartesian( elementType='Q1', elementRes=(16,16), minCoord=(0.,0.), maxCoord=(1.,1.) ) >>> someMesh.dim 2 >>> someMesh.elementRes (16, 16)
Alternatively, you can create a linear/constant dual mesh
>>> someDualMesh = uw.mesh.FeMesh_Cartesian( elementType='Q1/dQ0', elementRes=(16,16), minCoord=(0.,0.), maxCoord=(1.,1.) ) >>> someDualMesh.elementType 'Q1' >>> subMesh = someDualMesh.subMesh >>> subMesh.elementType 'DQ0'
To set / read vertex coords, use the numpy interface via the ‘data’ property.
-
integrate(fn)[source]¶ Perform a domain integral of the given underworld function over this mesh
Parameters: mesh (uw.mesh.FeMesh_Cartesian) – Domain to perform integral over. Examples
>>> mesh = uw.mesh.FeMesh_Cartesian(minCoord=(0.0,0.0), maxCoord=(1.0,2.0)) >>> fn_1 = uw.function.misc.constant(2.0) >>> np.allclose( mesh.integrate( fn_1 )[0], 4 ) True
>>> fn_2 = uw.function.misc.constant(2.0) * (0.5, 1.0) >>> np.allclose( mesh.integrate( fn_2 ), [2,4] ) True
-
subMesh¶ Returns: Returns the submesh where the object is a dual mesh, or None otherwise. Return type: underworld.mesh.FeMesh
-
class
underworld.mesh.MeshVariable(mesh, nodeDofCount, dataType='double', **kwargs)[source]¶ Bases:
underworld._stgermain.StgCompoundComponent,underworld.function._function.Function,underworld._stgermain.Save,underworld._stgermain.LoadThe MeshVariable class generates a variable supported by a finite element mesh.
To set / read nodal values, use the numpy interface via the ‘data’ property.
Parameters: - mesh (underworld.mesh.FeMesh) – The supporting mesh for the variable.
- dataType (string) – The data type for the variable. Note that only ‘double’ type variables are currently supported.
- nodeDofCount (int) – Number of degrees of freedom per node the variable will have.
See property docstrings for further information.
Example
For example, to create a scalar meshVariable:
>>> linearMesh = uw.mesh.FeMesh_Cartesian( elementType='Q1/dQ0', elementRes=(16,16), minCoord=(0.,0.), maxCoord=(1.,1.) ) >>> scalarFeVar = uw.mesh.MeshVariable( mesh=linearMesh, nodeDofCount=1, dataType="double" )
or a vector meshvariable can be created:
>>> vectorFeVar = uw.mesh.MeshVariable( mesh=linearMesh, nodeDofCount=3, dataType="double" )
-
copy(deepcopy=False)[source]¶ This method returns a copy of the meshvariable.
Parameters: deepcopy (bool) – If True, the variable’s data is also copied into new variable. Returns: The mesh variable copy. Return type: underworld.mesh.MeshVariable Example
>>> mesh = uw.mesh.FeMesh_Cartesian() >>> var = uw.mesh.MeshVariable(mesh,2) >>> import math >>> var.data[:] = (math.pi,math.exp(1.)) >>> varCopy = var.copy() >>> varCopy.mesh == var.mesh True >>> varCopy.nodeDofCount == var.nodeDofCount True >>> import numpy as np >>> np.allclose(var.data,varCopy.data) False >>> varCopy2 = var.copy(deepcopy=True) >>> np.allclose(var.data,varCopy2.data) True
-
data¶ Numpy proxy array to underlying variable data. Note that the returned array is a proxy for all the local nodal data, and is provided as 1d list. It is possible to change the shape of this numpy array to reflect the cartesian topology (where appropriate), though again only the local portion of the decomposed domain will be available, and the shape will not necessarily be identical on all processors.
As these arrays are simply proxys to the underlying memory structures, no data copying is required.
Returns: The proxy array. Return type: numpy.ndarray Example
>>> linearMesh = uw.mesh.FeMesh_Cartesian( elementType='Q1/dQ0', elementRes=(16,16), minCoord=(0.,0.), maxCoord=(1.,1.) ) >>> scalarFeVar = uw.mesh.MeshVariable( mesh=linearMesh, nodeDofCount=1, dataType="double" ) >>> scalarFeVar.data.shape (289, 1)
You can retrieve individual nodal values
>>> scalarFeVar.data[100] array([ 0.])
Likewise you can modify nodal values
>>> scalarFeVar.data[100] = 15.333 >>> scalarFeVar.data[100] array([ 15.333])
-
dataType¶ Returns: Data type for variable. Supported types are ‘double’. Return type: str
-
fn_gradient¶ Returns a Function for the gradient field of this meshvariable.
Note that for a scalar variable T, the gradient function returns an array of the form:
\[[ \frac{\partial T}{\partial x}, \frac{\partial T}{\partial y}, \frac{\partial T}{\partial z} ]\]and for a vector variable v:
\[[ \frac{\partial v_x}{\partial x}, \frac{\partial v_x}{\partial y}, \frac{\partial v_x}{\partial z}, \frac{\partial v_y}{\partial x}, \frac{\partial v_y}{\partial y}, \frac{\partial v_y}{\partial z}, \frac{\partial v_z}{\partial x}, \frac{\partial v_z}{\partial y}, \frac{\partial v_z}{\partial z} ]\]Returns: The gradient function. Return type: underworld.function.Function
-
load(filename, interpolate=False)[source]¶ Load the MeshVariable from disk.
Parameters: - filename (str) – The filename for the saved file. Relative or absolute paths may be used, but all directories must exist.
- interpolate (bool) – Set to True to interpolate a file containing different resolution data. Note that a temporary MeshVariable with the file data will be build on each processor. Also note that the temporary MeshVariable can only be built if its corresponding mesh file is available. Also note that the supporting mesh mush be regular.
Notes
This method must be called collectively by all processes.
If the file data array is the same length as the current variable global size, it is assumed the file contains compatible data. Note that this may not be the case where for example where you have saved using a 2*8 resolution mesh, then loaded using an 8*2 resolution mesh.
Provided files must be in hdf5 format, and use the correct schema.
Example
Refer to example provided for ‘save’ method.
-
mesh¶ Returns: Supporting FeMesh for this MeshVariable. Return type: underworld.mesh.FeMesh
-
nodeDofCount¶ Returns: Degrees of freedom on each mesh node that this variable provides. Return type: int
-
save(filename, meshHandle=None)[source]¶ Save the MeshVariable to disk.
Parameters: - filename (string) – The name of the output file. Relative or absolute paths may be used, but all directories must exist.
- meshHandle (uw.utils.SavedFileData , optional) – The saved mesh file handle. If provided, a link is created within the mesh variable file to this saved mesh file. Important for checkpoint when the mesh deforms.
Notes
This method must be called collectively by all processes.
Returns: Data object relating to saved file. This only needs to be retained if you wish to create XDMF files and can be ignored otherwise. Return type: underworld.utils.SavedFileData Example
First create the mesh add a variable:
>>> mesh = uw.mesh.FeMesh_Cartesian( elementType='Q1/dQ0', elementRes=(16,16), minCoord=(0.,0.), maxCoord=(1.,1.) ) >>> var = uw.mesh.MeshVariable( mesh=mesh, nodeDofCount=1, dataType="double" )
Write something to variable
>>> import numpy as np >>> var.data[:,0] = np.arange(var.data.shape[0])
Save to a file (note that the ‘ignoreMe’ object isn’t really required):
>>> ignoreMe = var.save("saved_mesh_variable.h5")
Now let’s try and reload.
>>> clone_var = uw.mesh.MeshVariable( mesh=mesh, nodeDofCount=1, dataType="double" ) >>> clone_var.load("saved_mesh_variable.h5")
Now check for equality:
>>> np.allclose(var.data,clone_var.data) True
>>> # clean up: >>> if uw.rank() == 0: ... import os; ... os.remove( "saved_mesh_variable.h5" )
-
syncronise()[source]¶ This method is often necessary when Underworld is operating in parallel.
It will syncronise the mesh variable so that it is consistent with it’s parallel neighbours. Specifically, the shadow space of each process obtains the required data from neighbouring processes.
-
xdmf(filename, fieldSavedData, varname, meshSavedData, meshname, modeltime=0.0)[source]¶ Creates an xdmf file, filename, associating the fieldSavedData file on the meshSavedData file
Notes
xdmf contain 2 files: an .xml and a .h5 file. See http://www.xdmf.org/index.php/Main_Page This method only needs to be called by the master process, all other processes return quiely.
Parameters: - filename (str) – The output path to write the xdmf file. Relative or absolute paths may be used, but all directories must exist.
- varname (str) – The xdmf name to give the field
- meshSavedData (underworld.utils.SaveFileData) – Handler returned for saving a mesh. underworld.mesh.save(xxx)
- meshname (str) – The xdmf name to give the mesh
- fieldSavedData (underworld.SavedFileData) – Handler returned from saving a field. underworld.mesh.save(xxx)
- modeltime (float) – The time recorded in the xdmf output file
Example
First create the mesh add a variable:
>>> mesh = uw.mesh.FeMesh_Cartesian( elementType='Q1/dQ0', elementRes=(16,16), minCoord=(0.,0.), maxCoord=(1.,1.) ) >>> var = uw.mesh.MeshVariable( mesh=mesh, nodeDofCount=1, dataType="double" )
Write something to variable
>>> import numpy as np >>> var.data[:,0] = np.arange(var.data.shape[0])
Save mesh and var to a file:
>>> meshDat = mesh.save("saved_mesh.h5") >>> varDat = var.save("saved_mesh_variable.h5")
Now let’s create the xdmf file
>>> var.xdmf("TESTxdmf", varDat, "var1", meshDat, "meshie" )
Does file exist?
>>> import os >>> if uw.rank() == 0: os.path.isfile("TESTxdmf.xdmf") True
Clean up:
>>> if uw.rank() == 0: ... import os; ... os.remove( "saved_mesh_variable.h5" ) ... os.remove( "saved_mesh.h5" ) ... os.remove( "TESTxdmf.xdmf" )
underworld.systems module¶
This module contains routines relating to differential system.
-
class
underworld.systems.sle.ConstitutiveMatrixTerm(fn_visc1=None, fn_visc2=None, fn_director=None, **kwargs)[source]¶ Bases:
underworld.systems.sle._assemblyterm.MatrixAssemblyTerm
-
class
underworld.systems.sle.PreconditionerMatrixTerm(assembledObject=None, **kwargs)[source]¶ Bases:
underworld.systems.sle._assemblyterm.MatrixAssemblyTerm
-
class
underworld.systems.sle.AssembledVector(meshVariable, eqNum, **kwargs)[source]¶ Bases:
underworld.systems.sle._svector.SolutionVectorVector object, generally assembled as a result of the FEM framework.
See parent class for parameters.
-
petscVector¶ Underlying PETSc vector object.
Type: petscVector (swig petsc vector)
-
-
class
underworld.systems.sle.AssemblyTerm(integrationSwarm, extraInfo=None, **kwargs)[source]¶ Bases:
underworld._stgermain.StgCompoundComponent
-
class
underworld.systems.sle.AssembledMatrix(rowVector, colVector, rhs=None, rhs_T=None, assembleOnNodes=False, **kwargs)[source]¶ Bases:
underworld._stgermain.StgCompoundComponentMatrix object, generally assembled as a result of the FEM framework.
Parameters: - meshVariableRow (underworld.mesh.MeshVariable) – MeshVariable object for matrix row.
- meshVariableCol (underworld.mesh.MeshVariable) – MeshVariable object for matrix column.
-
class
underworld.systems.sle.VectorAssemblyTerm(assembledObject, **kwargs)[source]¶ Bases:
underworld.systems.sle._assemblyterm.AssemblyTerm
-
class
underworld.systems.sle.MatrixAssemblyTerm_NA_i__NB_i__Fn(fn, **kwargs)[source]¶ Bases:
underworld.systems.sle._assemblyterm.MatrixAssemblyTerm
-
class
underworld.systems.sle.LumpedMassMatrixVectorTerm(assembledObject, **kwargs)[source]¶ Bases:
underworld.systems.sle._assemblyterm.VectorAssemblyTerm
-
class
underworld.systems.sle.VectorAssemblyTerm_NA_i__Fn_i(fn, mesh=None, **kwargs)[source]¶ Bases:
underworld.systems.sle._assemblyterm.VectorAssemblyTerm
-
class
underworld.systems.sle.GradientStiffnessMatrixTerm(assembledObject=None, **kwargs)[source]¶ Bases:
underworld.systems.sle._assemblyterm.MatrixAssemblyTerm
-
class
underworld.systems.sle.MatrixAssemblyTerm_NA__NB__Fn(fn, mesh, **kwargs)[source]¶ Bases:
underworld.systems.sle._assemblyterm.MatrixAssemblyTerm
-
class
underworld.systems.sle.EqNumber(meshVariable, removeBCs=True, **kwargs)[source]¶ Bases:
underworld._stgermain.StgClassThe SolutionVector manages the numerical solution vectors used by Underworld’s equation systems. Interface between meshVariables and systems.
Parameters: meshVariable (uw.mesh.MeshVariable) – MeshVariable object for which this equation numbering corresponds. Example
>>> linearMesh = uw.mesh.FeMesh_Cartesian( elementType='Q1/dQ0', elementRes=(4,4), minCoord=(0.,0.), maxCoord=(1.,1.) ) >>> tField = uw.mesh.MeshVariable( linearMesh, 1 ) >>> teqNum = uw.systems.sle.EqNumber( tField )
-
class
underworld.systems.sle.AdvDiffResidualVectorTerm(velocityField, diffusivity, sourceTerm, **kwargs)[source]¶ Bases:
underworld.systems.sle._assemblyterm.VectorAssemblyTerm
-
class
underworld.systems.sle.SolutionVector(meshVariable, eqNumber, **kwargs)[source]¶ Bases:
underworld._stgermain.StgCompoundComponentThe SolutionVector manages the numerical solution vectors used by Underworld’s equation systems. Interface between meshVariables and systems.
Parameters: - meshVariable (underworld.mesh.MeshVariable) – MeshVariable object for which this SLE vector corresponds.
- eqNumber (underworld.systems.sle.EqNumber) – Equation numbering object corresponding to this vector.
Example
>>> linearMesh = uw.mesh.FeMesh_Cartesian( elementType='Q1/dQ0', elementRes=(4,4), minCoord=(0.,0.), maxCoord=(1.,1.) ) >>> tField = uw.mesh.MeshVariable( linearMesh, 1 ) >>> eqNum = uw.systems.sle.EqNumber( tField ) >>> sVector = uw.systems.sle.SolutionVector(tField, eqNum )
-
class
underworld.systems.sle.MatrixAssemblyTerm(assembledObject=None, **kwargs)[source]¶ Bases:
underworld.systems.sle._assemblyterm.AssemblyTerm
-
class
underworld.systems.sle.VectorSurfaceAssemblyTerm_NA__Fn__ni(nbc, integrationSwarm=None, surfaceGaussPoints=2, **kwargs)[source]¶ Bases:
underworld.systems.sle._assemblyterm.VectorAssemblyTermBuild an assembly term for a surface integral.
Parameters: - nbc (underworld.conditions.NeumannCondition) – See uw.conditions.NeumannCondition for details
- integrationSwarm (underworld.swarm.GaussBorderIntegrationSwarm) – Optional integration swarm to be used for numerical integration.
- surfaceGaussPoints (int) – The number of quadrature points per element face to use in surface integration. Will be used to create a GaussBorderIntegrationSwarm in the case the ‘swarm’ input is ‘None’.
-
class
underworld.systems.sle.VectorAssemblyTerm_NA__Fn(fn, mesh=None, **kwargs)[source]¶ Bases:
underworld.systems.sle._assemblyterm.VectorAssemblyTerm
-
class
underworld.systems.sle.VectorAssemblyTerm_NA_j__Fn_ij(fn, mesh=None, **kwargs)[source]¶ Bases:
underworld.systems.sle._assemblyterm.VectorAssemblyTermBuild an assembly term for a spatial gradient, used for the viscoelastic force term.
Parameters: - fn (underworld.function.Function) – Function is a vector of size 3 (dim=2) or 6 (dim=3) representing a symetric tensor
- mesh (uw.mesh.FeMesh_Cartesian) –
underworld.systems.Solver |
This method simply returns a necessary solver for the provided system. |
underworld.systems.Stokes |
This class provides functionality for a discrete representation of the Stokes flow equations. |
underworld.systems.HeatSolver |
Steady State Heat Equation Solver. |
underworld.systems.SteadyStateHeat |
This class provides functionality for a discrete representation of the steady state heat equation. |
underworld.systems.StokesSolver |
The Block Stokes Schur Complement Solver: This solves the saddle-point system |
underworld.systems.SteadyStateDarcyFlow |
This class provides functionality for a discrete representation of the steady state darcy flow equation. |
underworld.systems.TimeIntegration |
Abstract class for integrating numerical objects (fields, swarms, etc.) in time. |
underworld.systems.AdvectionDiffusion |
This class provides functionality for a discrete representation of an advection-diffusion equation. |
underworld.systems.SwarmAdvector |
Objects of this class advect a swarm through time using the provided velocity field. |
-
class
underworld.systems.Stokes(velocityField, pressureField, fn_viscosity, fn_bodyforce=None, fn_one_on_lambda=None, fn_source=None, voronoi_swarm=None, conditions=[], _removeBCs=True, _fn_viscosity2=None, _fn_director=None, fn_stresshistory=None, _fn_stresshistory=None, _fn_v0=None, _fn_p0=None, _callback_post_solve=None, **kwargs)[source]¶ Bases:
underworld._stgermain.StgCompoundComponentThis class provides functionality for a discrete representation of the Stokes flow equations.
Specifically, the class uses a mixed 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 ‘velocityField’ and ‘pressureField’ parameters.
The strong form of the given boundary value problem, for \(f\), \(g\) and \(h\) given, is
\[\begin{split}\begin{align} \sigma_{ij,j} + f_i =& \: 0 & \text{ in } \Omega \\ u_{k,k} + \frac{p}{\lambda} =& \: H & \text{ in } \Omega \\ u_i =& \: g_i & \text{ on } \Gamma_{g_i} \\ \sigma_{ij}n_j =& \: h_i & \text{ on } \Gamma_{h_i} \\ \end{align}\end{split}\]where,
- \(\sigma_{i,j}\) is the stress tensor
- \(u_i\) is the velocity,
- \(p\) is the pressure,
- \(f_i\) is a body force,
- \(\lambda\) is pseudo compressibility factor,
- \(H\) is the compressible equation source term,
- \(g_i\) are the velocity boundary conditions (DirichletCondition)
- \(h_i\) are the traction boundary conditions (NeumannCondition).
The problem boundary, \(\Gamma\), admits the decompositions \(\Gamma=\Gamma_{g_i}\cup\Gamma_{h_i}\) where \(\emptyset=\Gamma_{g_i}\cap\Gamma_{h_i}\). The equivalent weak form is:
\[\int_{\Omega} w_{(i,j)} \sigma_{ij} \, d \Omega = \int_{\Omega} w_i \, f_i \, d\Omega + \sum_{j=1}^{n_{sd}} \int_{\Gamma_{h_j}} w_i \, h_i \, d \Gamma\]where we must find \(u\) which satisfies the above for all \(w\) in some variational space.
Parameters: - velocityField (underworld.mesh.MeshVariable) – Variable used to record system velocity.
- pressureField (underworld.mesh.MeshVariable) – Variable used to record system pressure.
- fn_viscosity (underworld.function.Function) – Function which reports a viscosity value. Function must return scalar float values.
- fn_bodyforce (underworld.function.Function, Default = None) – Function which reports a body force for the system. Function must return float values of identical dimensionality to the provided velocity variable.
- fn_one_on_lambda (underworld.function.Function, Default = None) – Pseudo-compressibility factor. Note that non-zero values are incompatible with the ‘penalty’ stokes solver. Ensure a ‘penalty’ equal to 0 is used if this function is non-zero. By default this is the case.
- fn_source (underworld.function.Function, Default = None) – Mass source term. Check fn_one_on_lambda for usage caveats.
- fn_stresshistory (underworld.function.Function, Default = None) – Function which defines the stress history term used for viscoelasticity. Function is a vector of size 3 (dim=2) or 6 (dim=3) representing a symetric tensor.
- voronoi_swarm (underworld.swarm.Swarm) – If a voronoi_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.
Notes
Constructor must be called by collectively all processes.
-
eqResiduals¶ Returns the stokes flow equations’ residuals from the latest solve. Residual calculations use the matrices and vectors of the discretised problem. The residuals correspond to the momentum equation and the continuity equation.
Returns: r1 is the momentum equation residual r2 is the continuity equation residual Return type: (r1, r2) - 2 tuple of doubles Notes
This method must be called collectively by all processes.
-
fn_bodyforce¶ The body force function. You may change this function directly via this property.
-
fn_one_on_lambda¶ A bulk viscosity parameter
-
fn_source¶ The volumetric source term function. You may change this function directly via this property.
-
fn_viscosity¶ The viscosity function. You may change this function directly via this property.
-
stokes_callback¶ Return the callback function used by this system
-
class
underworld.systems.HeatSolver(heatSLE, **kwargs)[source]¶ Bases:
underworld._stgermain.StgCompoundComponentSteady State Heat Equation Solver.
-
configure(solve_type='')[source]¶ Configure velocity/inner solver (A11 PETSc prefix).
solve_type can be one of:
- mumps : MUMPS parallel direct solver.
- superludist : SuperLU parallel direct solver.
- superlu : SuperLU direct solver (serial only).
- lu : LU direct solver (serial only).
-
solve(nonLinearIterate=None, nonLinearTolerance=0.01, nonLinearMaxIterations=500, callback_post_solve=None, **kwargs)[source]¶ Solve the HeatEq system
Parameters: - nonLinearIterate (bool) – True will perform non linear iterations iterations, False (or 0) will not
- nonLinearTolerance (float, Default=1.0e-2) – Relative tolerance criterion for the change in the velocity field
- nonLinearMaxIterations (int, Default=500) – Maximum number of non linear iteration to perform
- callback_post_sovle (func, Default=None) – Optional callback function to be performed at the end of a linear solve iteration. Commonly this will be used to perform operations between non linear iterations, for example, calibrating the solution or removing the system null space.
-
-
class
underworld.systems.SteadyStateHeat(temperatureField, fn_diffusivity, fn_heating=0.0, voronoi_swarm=None, conditions=[], _removeBCs=True, **kwargs)[source]¶ Bases:
underworld._stgermain.StgCompoundComponentThis class provides functionality for a discrete representation of the steady state heat 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 ‘temperatureField’.
The strong form of the given boundary value problem, for \(f\), \(h\) and \(h\) given, is
\[\begin{split}\begin{align} q_i =& - \alpha \, u_{,i} & \\ q_{i,i} =& \: f & \text{ in } \Omega \\ u =& \: g & \text{ on } \Gamma_g \\ -q_i n_i =& \: h & \text{ on } \Gamma_h \\ \end{align}\end{split}\]where, \(\alpha\) is the diffusivity, \(u\) is the temperature, \(f\) is a source term, \(g\) is the Dirichlet condition, and \(h\) is a Neumann condition. The problem boundary, \(\Gamma\), admits the decomposition \(\Gamma=\Gamma_g\cup\Gamma_h\) where \(\emptyset=\Gamma_g\cap\Gamma_h\). The equivalent weak form is:
\[-\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 \(u\) which satisfies the above for all \(w\) in some variational space.
Parameters: - temperatureField (underworld.mesh.MeshVariable) – The solution field for temperature.
- fn_diffusivity (underworld.function.Function) – The function that defines the diffusivity across the domain.
- fn_heating (underworld.function.Function) – A function that defines the heating across the domain. Optional.
- voronoi_swarm (underworld.swarm.Swarm) – If a voronoi_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.
Notes
Constructor must be called collectively by all processes.
Example
Setup a basic thermal system:
>>> linearMesh = uw.mesh.FeMesh_Cartesian( elementType='Q1/dQ0', elementRes=(4,4), minCoord=(0.,0.), maxCoord=(1.,1.) ) >>> tField = uw.mesh.MeshVariable( linearMesh, 1 ) >>> topNodes = linearMesh.specialSets["MaxJ_VertexSet"] >>> bottomNodes = linearMesh.specialSets["MinJ_VertexSet"] >>> tbcs = uw.conditions.DirichletCondition(tField, topNodes + bottomNodes) >>> tField.data[topNodes.data] = 0.0 >>> tField.data[bottomNodes.data] = 1.0 >>> tSystem = uw.systems.SteadyStateHeat(temperatureField=tField, fn_diffusivity=1.0, conditions=[tbcs])
Example with non diffusivity:
>>> k = tField + 1.0 >>> tSystem = uw.systems.SteadyStateHeat(temperatureField=tField, fn_diffusivity=k, conditions=[tbcs]) >>> solver = uw.systems.Solver(tSystem) >>> solver.solve() Traceback (most recent call last): ... RuntimeError: Nonlinearity detected. Diffusivity function depends on the temperature field provided to the system. Please set the 'nonLinearIterate' solve parameter to 'True' or 'False' to continue. >>> solver.solve(nonLinearIterate=True)
-
fn_diffusivity¶ The diffusivity function. You may change this function directly via this property.
-
fn_heating¶ The heating function. You may change this function directly via this property.
-
class
underworld.systems.StokesSolver(stokesSLE, **kwargs)[source]¶ Bases:
underworld._stgermain.StgCompoundComponentThe Block Stokes Schur Complement Solver: This solves the saddle-point system
\[\begin{split}\begin{bmatrix} K & G \\ G^T & C \end{bmatrix} \begin{bmatrix} u \\ p \end{bmatrix} = \begin{bmatrix}f \\ h \end{bmatrix}\end{split}\]via a Schur complement method.
We first solve:
(1)¶\[S p= G^T K^{-1} f - h,\]where \(S = G^T K^{-1} G-C\)
Then we backsolve for the velocity:
(2)¶\[K u = f - G p.\]The effect of \(K^{-1}\) in (1) is obtained via a KSPSolve in PETSc. This has the prefix ‘A11’ (often called the ‘inner’ solve)
The solve in (1) for the pressure has prefix ‘scr’.
Assuming the returned solver is called ‘solver’, it is possible to configure these solves individually via the solver.options.A11 and solver.options.scr dictionaries.
Try help(solver.options.A11) for some details.
Common configurations are provided via the set_inner_method() method.
help(solver.set_inner_method) for more.
For more advanced configurations use the solver.options.A11/scr dictionaries directly.
help(solver.options) to see more.
-
set_inner_method(solve_type='mg')[source]¶ Configure velocity/inner solver (A11 PETSc prefix).
Available options below. Note that associated solver software (for example mumps) must be installed.
- mg : Geometric multigrid (default).
- nomg : Disables multigrid.
- lu : LU direct solver (serial only).
- mumps : MUMPS parallel direct solver.
- superludist : SuperLU parallel direct solver.
- superlu : SuperLU direct solver (serial only).
-
set_mg_levels(levels)[source]¶ Set the number of multigrid levels manually. It is set automatically by default.
-
set_penalty(penalty)[source]¶ By setting the penalty, the Augmented Lagrangian Method is used as the solve. This method is not recommended for normal use as there is additional memory and cpu overhead. This method can often help improve convergence issues for problems with large viscosity contrasts that are having trouble converging.
A penalty of roughly 0.1 of the maximum viscosity contrast is not a bad place to start as a rule of thumb. (check notes/paper)
-
solve(nonLinearIterate=None, nonLinearTolerance=0.01, nonLinearKillNonConvergent=False, nonLinearMinIterations=1, nonLinearMaxIterations=500, callback_post_solve=None, print_stats=False, reinitialise=True, **kwargs)[source]¶ Solve the stokes system
Parameters: - nonLinearIterate (bool) – True will perform non linear iterations iterations, False (or 0) will not
- nonLinearTolerance (float, Default=1.0e-2) – Relative tolerance criterion for the change in the velocity field
- nonLinearMaxIterations (int, Default=500) – Maximum number of non linear iteration to perform
- callback_post_sovle (func, Default=None) – Optional callback function to be performed at the end of a linear solve iteration. Commonly this will be used to perform operations between non linear iterations, for example, calibrating the pressure solution or removing the system null space.
- print_stats (bool, Default=False) – Print out solver iteration and timing counts per solver
- reinitialise (bool, Default=True,) – Rebuild the system discretisation storage (location matrix/petsc mats & vecs) and repopulate, if available, the stokes voronio swarm before the system is solved.
-
-
class
underworld.systems.SteadyStateDarcyFlow(pressureField, fn_diffusivity, fn_bodyforce=None, voronoi_swarm=None, conditions=[], velocityField=None, swarmVarVelocity=None, _removeBCs=True, **kwargs)[source]¶ Bases:
underworld._stgermain.StgCompoundComponentThis 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 \(f\), \(q\) and \(h\) given, is
\[\begin{split}\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}\end{split}\]where,
- \(\kappa\) is the diffusivity, \(u\) is the pressure,
- \(S\) is a flow body-source, for example due to gravity,
- \(f\) is a source term, \(q\) is the Dirichlet condition, and
- \(h\) is a Neumann condition.
The problem boundary, \(\Gamma\), admits the decomposition \(\Gamma=\Gamma_q\cup\Gamma_h\) where \(\emptyset=\Gamma_q\cap\Gamma_h\). The equivalent weak form is:
\[-\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 \(u\) which satisfies the above for all \(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.
- 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.
-
fn_bodyforce¶ The heating function. You may change this function directly via this property.
-
fn_diffusivity¶ The diffusivity function. You may change this function directly via this property.
-
class
underworld.systems.TimeIntegration(order, **kwargs)[source]¶ Bases:
underworld._stgermain.StgCompoundComponentAbstract class for integrating numerical objects (fields, swarms, etc.) in time.
The integration algorithm is a modified Runge Kutta method that only evaluates midpoint information varying in space - using only the present timestep solution. The order of the integration used can be 1,2,4
Parameters: order (int {1,2,4}) – Defines the numerical order ‘in space’ of the Runge Kutta like integration scheme. -
dt¶ Time integrator timestep size.
-
time¶ Time integrator time value.
-
-
class
underworld.systems.AdvectionDiffusion(phiField, phiDotField, velocityField, fn_diffusivity, fn_sourceTerm=None, conditions=[], _allow_non_q1=False, **kwargs)[source]¶ Bases:
underworld._stgermain.StgCompoundComponentThis class provides functionality for a discrete representation of an advection-diffusion equation.
The class uses the Streamline Upwind Petrov Galerkin SUPG method to integrate through time.
\[\frac{\partial\phi}{\partial t} + {\bf u } \cdot \nabla \phi= \nabla { ( k \nabla \phi ) } + H\]Parameters: - phiField (underworld.mesh.MeshVariable) – The concentration field, typically the temperature field
- phiDotField (underworld.mesh.MeshVariable) – A MeshVariable that defines the initial time derivative of the phiField. Typically 0 at the beginning of a model, e.g. phiDotField.data[:]=0 When using a phiField loaded from disk one should also load the phiDotField to ensure the solving method has the time derivative information for a smooth restart. No dirichlet conditions are required for this field as the phiField degrees of freedom map exactly to this field’s dirichlet conditions, the value of which ought to be 0 for constant values of phi.
- velocityField (underworld.mesh.MeshVariable) – The velocity field.
- fn_diffusivity (underworld.function.Function) – A function that defines the diffusivity within the domain.
- fn_sourceTerm (underworld.function.Function) – A function that defines the heating within the domain. Optional.
- 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.
Notes
Constructor must be called by collectively all processes.
-
class
underworld.systems.SwarmAdvector(velocityField, swarm, order=2, **kwargs)[source]¶ Bases:
underworld.systems._timeintegration.TimeIntegrationObjects of this class advect a swarm through time using the provided velocity field.
Parameters: - velocityField (underworld.mesh.MeshVariable) – The MeshVariable field used for evaluating the velocity field that advects the swarm particles
- swarm (underworld.swarm.Swarm) – Particle swarm that will be advected by the given velocity field
-
integrate(dt, update_owners=True)[source]¶ Integrate the associated swarm in time, by dt, using the velocityfield that is associated with this class
Parameters: - dt (double) – The timestep to use in the intergration
- update_owners (bool) – If this is set to False, particle ownership (which element owns a particular particle) is not updated after advection. This is often necessary when both the mesh and particles are advecting simutaneously.
Example
>>> import underworld as uw >>> import numpy as np >>> from underworld import function as fn >>> dim=2; >>> elementMesh = uw.mesh.FeMesh_Cartesian(elementType="Q1/dQ0", elementRes=(9,9), minCoord=(-1.,-1.), maxCoord=(1.,1.)) >>> velocityField = uw.mesh.MeshVariable( mesh=elementMesh, nodeDofCount=dim ) >>> swarm = uw.swarm.Swarm(mesh=elementMesh) >>> particle = np.zeros((1,2)) >>> particle[0] = [0.2,-0.2] >>> swarm.add_particles_with_coordinates(particle) array([0], dtype=int32) >>> velocityField.data[:]=[1.0,1.0] >>> swarmAdvector = uw.systems.SwarmAdvector(velocityField=velocityField, swarm=swarm, order=2) >>> dt=swarmAdvector.get_max_dt() >>> swarmAdvector.integrate(dt) >>> np.allclose(swarm.particleCoordinates.data[0], [ 0.27856742, -0.12143258], rtol=1e-4) True
underworld.timing module¶
This module implements some high level timing operations for Underworld, allowing users to determine how walltime is divided between different Underworld API calls. Note that this module only records timing for Underworld API calls, and has no way of knowing how much time has been spent elsewhere (such as numpy, scipy etc). The total runtime is also recorded which gives users an indication of how much time is spent outside Underworld.
Timing routines enabled by this module should introduce negligible computational overhead.
Only the root process records timing information.
Note that to utilise timing routines, you must first set the ‘UW_ENABLE_TIMING’ environment variable, and this must be done before you call import underworld.
>>> import os
>>> os.environ["UW_ENABLE_TIMING"] = "1"
>>> import underworld as uw
>>> uw.timing.start()
>>> someMesh = uw.mesh.FeMesh_Cartesian()
>>> with someMesh.deform_mesh():
... someMesh.data[0] = [0.1,0.1]
>>> uw.timing.stop()
>>> # uw.print_table() # This will print the data.
>>> # Commented out as not doctest friendly.
underworld.timing.log_result |
Allows the user to manually add entries to data. |
underworld.timing.reset |
Reset timing data. |
underworld.timing.get_data |
Returns dict with timing data. |
underworld.timing.stop |
Call this function to stop recording timing data. |
underworld.timing.print_table |
Print timing results to stdout or to a provided file. |
underworld.timing.start |
Call this function to start recording timing data. |
-
underworld.timing.log_result(time, name)[source]¶ Allows the user to manually add entries to data.
Parameters: - time (float) – Time spent.
- name (str) – Name to record to dataset. Note that the current stack information is generated internally and recorded.
-
underworld.timing.reset()[source]¶ Reset timing data. Note that this function calls stop(), and the user must call start() to resume recording timing data.
-
underworld.timing.get_data(group_by='line_routine')[source]¶ Returns dict with timing data.
Parameters: group_by (str) – Reported timing data is grouped according to the following options: “line” : Calling line of code. “routine” : Class routine. “line_routine”: Line&routine form an individual timing group.
-
underworld.timing.stop()[source]¶ Call this function to stop recording timing data. Note that this is automatically called when print_table() is called.
-
underworld.timing.print_table(group_by='line_routine', sort_by='total', display_fraction=0.95, float_precision='.3f', output_file=None, **kwargs)[source]¶ Print timing results to stdout or to a provided file. Call this function stops timing.
Parameters: - group_by (str) – See get_data() function
- sort_by (str) – Data is sorted according to: “total” : Total time allocated to any group. “average” : Average time attributed to any group
- display_fraction (float) – Set this option to cull insignificant (short time) results.
- output_file (str) – File to record table to. If none provided, outputs to stdout.
- **kwargs – Any extra kwargs are passed to tabulate module (if installed). This allows you to tweak the output format. Consule the tabulate module instructions for details.
underworld.conditions module¶
Implementation relating to system conditions.
underworld.conditions.NeumannCondition |
This class defines Neumann conditions for a differential equation. |
underworld.conditions.DirichletCondition |
The DirichletCondition class provides the required functionality to imposed Dirichlet conditions on your differential equation system. |
underworld.conditions.SystemCondition |
-
class
underworld.conditions.NeumannCondition(variable, indexSetsPerDof=None, fn_flux=None)[source]¶ Bases:
underworld.conditions._conditions.SystemConditionThis class defines Neumann conditions for a differential equation. Neumann conditions specifiy a field’s flux along a boundary.
As such the user specifices the field’s flux as a uw.Function and the nodes where this flux is to be applied - similar to uw.conditions.DirichletCondtion
Parameters: - fn_flux (underworld.function.Function) – Function which determines flux values.
- variable (underworld.mesh.MeshVariable) – The variable that describes the discretisation (mesh & DOFs) for ‘indexSetsPerDof’
- indexSetsPerDof (list, tuple, IndexSet) – The index set(s) which flag nodes/DOFs as Neumann conditions. Note that the user must provide an index set for each degree of freedom of the variable above. So for a vector variable of rank 2 (say Vx & Vy), two index sets must be provided (say VxDofSet, VyDofSet).
Example
Basic setup and usage of Neumann conditions:
>>> linearMesh = uw.mesh.FeMesh_Cartesian( elementType='Q1/dQ0', elementRes=(4,4), minCoord=(0.,0.), maxCoord=(1.,1.) ) >>> velocityField = uw.mesh.MeshVariable( linearMesh, 2 ) >>> velocityField.data[:] = [0.,0.] # set velocity zero everywhere, which will of course include the boundaries. >>> myFunc = (uw.function.coord()[1],0.0) >>> bottomWall = linearMesh.specialSets["MinJ_VertexSet"] >>> tractionBC = uw.conditions.NeumannCondition(variable=velocityField, fn_flux=myFunc, indexSetsPerDof=(None,bottomWall) )
-
fn_flux¶ Get the underworld.Function that defines the flux
-
class
underworld.conditions.DirichletCondition(variable, indexSetsPerDof)[source]¶ Bases:
underworld.conditions._conditions.SystemConditionThe DirichletCondition class provides the required functionality to imposed Dirichlet conditions on your differential equation system.
The user is simply required to flag which nodes/DOFs should be considered by the system to be a Dirichlet condition. The values at the Dirichlet nodes/DOFs is then left untouched by the system.
Parameters: - variable (underworld.mesh.MeshVariable) – This is the variable for which the Dirichlet condition applies.
- indexSetsPerDof (list, tuple, IndexSet) – The index set(s) which flag nodes/DOFs as Dirichlet conditions. Note that the user must provide an index set for each degree of freedom of the variable. So for a vector variable of rank 2 (say Vx & Vy), two index sets must be provided (say VxDofSet, VyDofSet).
Notes
Note that it is necessary for the user to set the required value on the variable, possibly via the numpy interface.
Constructor must be called collectively all processes.
Example
Basic setup and usage of Dirichlet conditions:
>>> linearMesh = uw.mesh.FeMesh_Cartesian( elementType='Q1/dQ0', elementRes=(4,4), minCoord=(0.,0.), maxCoord=(1.,1.) ) >>> velocityField = uw.mesh.MeshVariable( linearMesh, 2 ) >>> velocityField.data[:] = [0.,0.] # set velocity zero everywhere, which will of course include the boundaries. >>> IWalls = linearMesh.specialSets["MinI_VertexSet"] + linearMesh.specialSets["MaxI_VertexSet"] # get some wall index sets >>> JWalls = linearMesh.specialSets["MinJ_VertexSet"] + linearMesh.specialSets["MaxJ_VertexSet"] >>> freeSlipBC = uw.conditions.DirichletCondition(velocityField, (IWalls,JWalls) ) # this will give free slip sides >>> noSlipBC = uw.conditions.DirichletCondition(velocityField, (IWalls+JWalls,IWalls+JWalls) ) # this will give no slip sides
functions:¶
underworld.matplotlib_inline |
This function simply enables Jupyter Notebook inlined matplotlib results. |
underworld.nProcs |
Returns the number of processes being utilised by the simulation. |
underworld.rank |
Returns the rank of the current process. |
underworld.barrier |
Creates an MPI barrier. |
Module Details¶
functions:¶
-
underworld.matplotlib_inline()[source]¶ This function simply enables Jupyter Notebook inlined matplotlib results. This function should be called at the start of your notebooks as a replacement for the Jupyter Notebook %matplotlib inline magic. It provides the same functionality, however it allows notebooks to be converted to python without having to explicitly remove these calls.
-
underworld.nProcs()[source]¶ Returns the number of processes being utilised by the simulation.
Returns: Number of processors. Return type: unsigned
glucifer module¶
The glucifer module provides visualisation algorithms for Underworld.
Visualisation data is generated in parallel, with each processes generating the necessary data for its part of the domain. This data is written into a data file.
Actual rendering is performed in serial using the LavaVu rendering engine.
glucifer provides many flexible rendering options, including client-server based operation for remote usage. Users may choose to renderer outputs to raster images, or save a database file for later rendering. For those working in the Jupyter environment, glucifer will inline rendered images or even interactive webgl frames (still experimental).
Module Summary¶
submodules:¶
glucifer.objects module¶
glucifer.objects.ColourBar |
The ColourBar drawing object draws a colour bar for the provided colour map. |
glucifer.objects.VectorArrows |
This drawing object class draws vector arrows corresponding to the provided vector field. |
glucifer.objects.Points |
This drawing object class draws a swarm of points. |
glucifer.objects.Sampler |
The Sampler class provides functionality for sampling a field at a number of provided vertices. |
glucifer.objects.Surface |
This drawing object class draws a surface using the provided scalar field. |
glucifer.objects.Volume |
This drawing object class draws a volume using the provided scalar field. |
glucifer.objects.Mesh |
This drawing object class draws a mesh. |
glucifer.objects.Contours |
This drawing object class draws contour lines in a cross section using the provided scalar field. |
glucifer.objects.ColourMap |
The ColourMap class provides functionality for mapping colours to numerical values. |
glucifer.objects.CrossSection |
This drawing object class defines a cross-section plane, derived classes plot data over this cross section |
glucifer.objects.Drawing |
This is the base class for all drawing objects but can also be instantiated as is for direct/custom drawing. |
glucifer.objects.IsoSurface |
This drawing object class draws an isosurface using the provided scalar field. |
-
class
glucifer.objects.ColourBar(colourMap, *args, **kwargs)[source]¶ Bases:
glucifer.objects.DrawingThe ColourBar drawing object draws a colour bar for the provided colour map.
Parameters: colourMap (glucifer.objects.ColourMap) – Colour map for which the colour bar will be drawn.
-
class
glucifer.objects.VectorArrows(mesh, fn, resolution=[16, 16, 16], autoscale=True, *args, **kwargs)[source]¶ Bases:
glucifer.objects._GridSampler3DThis drawing object class draws vector arrows corresponding to the provided vector field.
See parent class for further parameter details. Also see property docstrings.
Parameters: - mesh (underworld.mesh.FeMesh) – Mesh over which vector arrows are rendered.
- fn (underworld.function.Function) – Function used to determine vectors to render. Function should return a vector of floats/doubles of appropriate dimensionality.
- arrowHead (float) – The size of the head of the arrow. If > 1.0 is ratio to arrow radius If in range [0.,1.] is ratio to arrow length
- scaling (float) – Scaling for entire arrow.
- autoscale (bool) – Scaling based on field min/max
- glyphs (int) – Type of glyph to render for vector arrow. 0: Line, 1 or more: 3d arrow, higher number => better quality.
- resolution (list(unsigned)) – Number of samples in the I,J,K directions.
-
class
glucifer.objects.Points(swarm, fn_colour=None, fn_mask=None, fn_size=None, colourVariable=None, colourBar=True, *args, **kwargs)[source]¶ Bases:
glucifer.objects.DrawingThis drawing object class draws a swarm of points.
See parent class for further parameter details. Also see property docstrings.
Parameters: - swarm (underworld.swarm.Swarm) – Swarm which provides locations for point rendering.
- fn_colour (underworld.function.Function) – Function used to determine colour to render particle. This function should return float/double values.
- fn_mask (underworld.function.Function) – Function used to determine if a particle should be rendered. This function should return bool values.
- fn_size (underworld.function.Function) – Function used to determine size to render particle. This function should return float/double values.
-
class
glucifer.objects.Sampler(mesh, fn, *args, **kwargs)[source]¶ Bases:
glucifer.objects.DrawingThe Sampler class provides functionality for sampling a field at a number of provided vertices.
Parameters: - vertices (list,array) – List of vertices to sample the field at, either a list or numpy array
- mesh (underworld.mesh.FeMesh) – Mesh over which the values are sampled
- fn (underworld.function.Function) – Function used to get the sampled values.
-
class
glucifer.objects.Surface(mesh, fn, drawSides='xyzXYZ', colourBar=True, onMesh=True, *args, **kwargs)[source]¶ Bases:
glucifer.objects.CrossSectionThis drawing object class draws a surface using the provided scalar field.
See parent class for further parameter details. Also see property docstrings.
Parameters: - mesh (underworld.mesh.FeMesh) – Mesh over which cross section is rendered.
- fn (underworld.function.Function) – Function used to determine values to render.
- drawSides (str) – Sides (x,y,z,X,Y,Z) for which the surface should be drawn. For example, “xyzXYZ” would render the provided function across all surfaces of the domain in 3D. In 2D, this object always renders across the entire domain.
-
class
glucifer.objects.Volume(mesh, fn, resolution=[64, 64, 64], colourBar=True, *args, **kwargs)[source]¶ Bases:
glucifer.objects._GridSampler3DThis drawing object class draws a volume using the provided scalar field.
See parent class for further parameter details. Also see property docstrings.
Parameters: - mesh (underworld.mesh.FeMesh) – Mesh over which object is rendered.
- fn (underworld.function.Function) – Function used to determine colour values. Function should return a vector of floats/doubles of appropriate dimensionality.
- resolution (list(unsigned)) – Number of samples in the I,J,K directions.
-
class
glucifer.objects.Mesh(mesh, nodeNumbers=False, segmentsPerEdge=1, *args, **kwargs)[source]¶ Bases:
glucifer.objects.DrawingThis drawing object class draws a mesh.
See parent class for further parameter details. Also see property docstrings.
Parameters: - mesh (underworld.mesh.FeMesh) – Mesh to render.
- nodeNumbers (bool) – Bool to determine whether global node numbers should be rendered.
- segmentsPerEdge (unsigned) – Number of segments to render per cell/element edge. For higher order mesh, more segments are useful to render mesh curvature correctly.
-
class
glucifer.objects.Contours(mesh, fn, labelFormat='', unitScaling=1.0, interval=0.33, limits=(0.0, 0.0), *args, **kwargs)[source]¶ Bases:
glucifer.objects.CrossSectionThis drawing object class draws contour lines in a cross section using the provided scalar field.
See parent class for further parameter details. Also see property docstrings.
Parameters: - mesh (underworld.mesh.FeMesh) – Mesh over which cross section is rendered.
- fn (underworld.function.Function) – Function used to determine values to render.
- labelFormat (str) – Format string (printf style) used to print a contour label, eg: ” %g K”
- unitScaling – Scaling factor to apply to value when printing labels
- interval (float) – Interval between contour lines
- limits (tuple, list) – User defined minimum and maximum limits for the contours. Provided as a tuple/list of floats (minValue, maxValue). If none is provided, the limits will be determined automatically.
-
class
glucifer.objects.ColourMap(colours='diverge', valueRange=None, logScale=False, discrete=False, **kwargs)[source]¶ Bases:
underworld._stgermain.StgCompoundComponentThe ColourMap class provides functionality for mapping colours to numerical values.
Parameters: - colours (str, list) – List of colours to use for drawing object colour map. Provided as a string or as a list of strings. Example, “red blue”, or [“red”, “blue”]
- valueRange (tuple, list) – User defined value range to apply to colour map. Provided as a tuple of floats (minValue, maxValue). If none is provided, the value range will be determined automatically.
- logScale (bool) – Bool to determine if the colourMap should use a logarithmic scale.
- discrete (bool) – Bool to determine if a discrete colour map should be used. Discrete colour maps do not interpolate between colours and instead use nearest neighbour for colouring.
-
class
glucifer.objects.CrossSection(mesh, fn, crossSection='', resolution=[100, 100, 1], colourBar=True, offsetEdges=None, onMesh=False, *args, **kwargs)[source]¶ Bases:
glucifer.objects.DrawingThis drawing object class defines a cross-section plane, derived classes plot data over this cross section
See parent class for further parameter details. Also see property docstrings.
Parameters: - mesh (underworld.mesh.FeMesh) – Mesh over which cross section is rendered.
- fn (underworld.function.Function) – Function used to determine values to render.
- crossSection (str) – Cross Section definition, eg. z=0.
- resolution (list(unsigned)) – Surface sampling resolution.
- onMesh (boolean) – Sample the mesh nodes directly, as opposed to sampling across a regular grid. This flag should be used in particular where a mesh has been deformed.
-
crossSection¶ Cross Section definition, eg;: z=0.
Type: crossSection (str)
-
class
glucifer.objects.Drawing(name=None, colours=None, colourMap='', colourBar=False, valueRange=None, logScale=False, discrete=False, *args, **kwargs)[source]¶ Bases:
underworld._stgermain.StgCompoundComponentThis is the base class for all drawing objects but can also be instantiated as is for direct/custom drawing.
Note that the defaults here are often overridden by the child objects.
Parameters: - colours (str, list.) – See ColourMap class docstring for further information
- colourMap (glucifer.objects.ColourMap) – A ColourMap object for the object to use. This should not be specified if ‘colours’ is specified.
- opacity (float) – Opacity of object. If provided, must take values from 0. to 1.
- colourBar (bool) – Bool to determine if a colour bar should be rendered.
- valueRange (tuple, list) – See ColourMap class docstring for further information
- logScale (bool) – See ColourMap class docstring for further information
- discrete (bool) – See ColourMap class docstring for further information
-
colourBar¶ return colour bar of drawing object, create if doesn’t yet exist.
Type: colourBar (object)
-
colourMap¶ return colour map of drawing object
Type: colourMap (object)
-
label(text, pos=(0.0, 0.0, 0.0), font='sans', scaling=1)[source]¶ Writes a label string
Parameters: - text (str) – label text.
- pos (tuple) – X,Y,Z position to place the label.
- font (str) – label font (small/fixed/sans/serif/vector).
- scaling (float) – label font scaling (for “vector” font only).
-
line(start=(0.0, 0.0, 0.0), end=(0.0, 0.0, 0.0))[source]¶ Draws a line
Parameters: - start (tuple) – X,Y,Z position to start line
- end (tuple) – X,Y,Z position to end line
-
class
glucifer.objects.IsoSurface(mesh, fn, fn_colour=None, resolution=[64, 64, 64], colourBar=True, isovalue=None, *args, **kwargs)[source]¶ Bases:
glucifer.objects.VolumeThis drawing object class draws an isosurface using the provided scalar field.
See parent class for further parameter details. Also see property docstrings.
Parameters: - mesh (underworld.mesh.FeMesh) – Mesh over which object is rendered.
- fn (underworld.function.Function) – Function used to determine surface position. Function should return a vector of floats/doubles.
- fn_colour (underworld.function.Function) – Function used to determine colour of surface.
- resolution (list(unsigned)) – Number of samples in the I,J,K directions.
- isovalue (float) – Isovalue to plot.
- isovalues (list of float) – List of multiple isovalues to plot.
classes:¶
glucifer.Store |
The Store class provides a database which stores gLucifer drawing objects as they are rendered in figures. |
glucifer.Figure |
The Figure class provides a window within which gLucifer drawing objects may be rendered. |
Module Details¶
classes:¶
-
class
glucifer.Store(filename=None, split=False, compress=True, **kwargs)[source]¶ Bases:
underworld._stgermain.StgCompoundComponentThe Store class provides a database which stores gLucifer drawing objects as they are rendered in figures. It also provides associated routines for saving and reloading this database to external files
In addition to parameter specification below, see property docstrings for further information.
Parameters: - filename (str) – Filename to use for a disk database, default is to create a temporary database filename.
- split (bool) – Set to true to write a separate database file for each timestep visualised
- view (bool) – Set to true and pass filename if loading a saved database for revisualisation
- compress (bool) – Set to true to enable database compression.
Example
Create a database:
>>> import glucifer >>> store = glucifer.Store()
Optionally provide a filename so you don’t need to call save later (no extension required)
>>> store = glucifer.Store('myvis')
Pass to figures when creating them (Providing a name allows you to revisualise the figure from the name)
>>> fig = glucifer.Figure(store, name="myfigure")
When figures are rendered with show() or save(imgname), they are saved to storage If you don’t need to render an image but still want to store the figure to view later, just call save() without a filename
>>> fig.save()
Save the database (only necessary if no filename provided when created)
>>> dbfile = store.save("myvis")
-
class
glucifer.Figure(store=None, name=None, figsize=None, boundingBox=None, facecolour='white', edgecolour='black', title='', axis=False, quality=1, *args, **kwargs)[source]¶ Bases:
dictThe Figure class provides a window within which gLucifer drawing objects may be rendered. It also provides associated routines for image generation and visualisation.
In addition to parameter specification below, see property docstrings for further information.
Parameters: - store (glucifer.Store) – Database to collect visualisation data, this may be shared among figures to collect their data into a single file.
- name (str) – Name of this figure, optional, used for revisualisation of stored figures.
- resolution (tuple) – Image resolution provided as a tuple.
- boundingBox (tuple) – Tuple of coordinate tuples defining figure bounding box. For example ( (0.1,0.1), (0.9,0.9) )
- facecolour (str) – Background colour for figure.
- edgecolour (str) – Edge colour for figure.
- title (str) – Figure title.
- axis (bool) – Bool to determine if figure axis should be drawn.
- quality (unsigned) – Antialiasing oversampling quality. For a value of 2, the image will be rendered at twice the resolution, and then downsampled. Setting this to 1 disables antialiasing, values higher than 3 are not recommended..
- properties (str) – Further properties to set on the figure.
Example
Create a figure:
>>> import glucifer >>> fig = glucifer.Figure()
We need a mesh
>>> import underworld as uw >>> mesh = uw.mesh.FeMesh_Cartesian()
Add drawing objects:
>>> fig.append(glucifer.objects.Surface( mesh, 1.))
Draw image. Note that if called from within a Jupyter notebook, image will be rendered inline. Otherwise, image will be saved to disk.
>>> fig.show()
Save the image
>>> imgfile = fig.save("test_image")
Clean up:
>>> if imgfile: ... import os; ... os.remove( imgfile )
-
append(drawingObject)[source]¶ Add a drawing object to the figure.
Parameters: drawingObject (glucifer.objects.Drawing) – The drawing object to add to the figure
-
axis¶ Axis enabled if true.
Type: axis
-
edgecolour¶ colour of figure border.
Type: edgecolour
-
facecolour¶ colour of face background.
Type: facecolour
-
objects¶ list of objects to be drawn within the figure.
Type: objects
-
properties¶ visual properties of viewport, passed to LavaVu to control rendering output of figure.
When using the property setter, new properties are set, overwriting any duplicate keys but keeping existing values otherwise.
Type: properties (dict)
-
resolution¶ size of window in pixels.
Type: resolution (tuple(int,int))
-
save(filename=None, size=(0, 0), type='Image')[source]¶ Saves the generated image to the provided filename or the figure to the database.
Parameters: - filename (str) – Filename to save file to. May include an absolute or relative path.
- (tuple(int,int)) (size) – If omitted, simply saves the figure data without generating an image
- type (str) – Type of visualisation to save (‘Image’ or ‘WebGL’).
Returns: filename – The final filename (including extension) used to save the image will be returned. Note that only the root process will return this filename. All other processes will not return anything.
Return type: str
-
script(cmd=None)[source]¶ Append to or get contents of the saved script.
Parameters: cmd (str) – Command to add to script.
-
send_command(cmd, retry=True)[source]¶ Run command on an open viewer instance.
Parameters: cmd (str) – Command to send to open viewer.
-
show(type='Image')[source]¶ Shows the generated image inline within an ipython notebook.
Parameters: - type (str) – Type of visualisation to display (‘Image’ or ‘WebGL’).
- IPython is installed, displays the result image or WebGL content inline (If) –
- IPython is not installed, this method will call the default image/web (If) –
- routines to save the result with a default filename in the current directory (output) –
-
static
show_grid(*rows)[source]¶ Shows a set of Figure objects in a grid. Note that this method currently only supports rendering images within a Jupyter Notebook, and saving gridded images to a file is not currently supported.
Parameters: rows (set of tuples) – Each provided tuple represents a row of Figures, and should only contain Figure class objects. Example
Create a bunch of figures: >>> import glucifer >>> fig1 = glucifer.Figure() >>> fig2 = glucifer.Figure() >>> fig3 = glucifer.Figure() >>> fig4 = glucifer.Figure() >>> fig5 = glucifer.Figure() >>> fig6 = glucifer.Figure()
We need a mesh >>> import underworld as uw >>> mesh = uw.mesh.FeMesh_Cartesian()
Add drawing objects as usual: >>> r = uw.function.input() >>> fig1.append(glucifer.objects.Surface( mesh, 1.)) >>> fig2.append(glucifer.objects.Mesh( mesh )) >>> fig3.append(glucifer.objects.Mesh( mesh, nodeNumbers=True )) >>> fig4.append(glucifer.objects.Surface( mesh, r[0])) >>> fig5.append(glucifer.objects.Surface( mesh, r[1])) >>> fig6.append(glucifer.objects.VectorArrows( mesh, r ))
Draw images in a grid. Note that in a Jupyter notebook, this will render the image within the notebook, though it will not be rendered in this example. Also note that show_grid() is a static method, and so we call it directly as below (instead of as a method on a Figure instance).
>>> glucifer.Figure.show_grid( (fig1,fig2,fig3), ... (fig4,fig5,fig6) ) <IPython.core.display.HTML object>
The above should generate a 2x3 (row x col) grid. For a 3x2 grid we would instead call:
>>> glucifer.Figure.show_grid( (fig1,fig2), ... (fig3,fig4), ... (fig5,fig6) ) <IPython.core.display.HTML object>
-
step¶ current timestep
Type: step (int)
-
title¶ a title for the image.
Type: title