anuga.Quantity

class anuga.Quantity(domain: Generic_Domain, vertex_values: ArrayLike | None = None, name: str | None = None, register: bool = False, qty_type: str | None = None)

Bases: object

Class Quantity - Implements values at each triangular element

__init__(domain: Generic_Domain, vertex_values: ArrayLike | None = None, name: str | None = None, register: bool = False, qty_type: str | None = None) → None

Create Quantity object

Parameters:

• domain – Associated domain structure. Required.

• vertex_values – N x 3 array of values at each vertex for each element. Default None

• name (str) – Provides a way to refer to a created quantity

• register – Register a quantity

• qty_type (str) –

  Memory allocation strategy. One of:

  • 'evolved' (default) — centroid, edge, explicit/semi-implicit update, centroid backup. x/y-gradient and phi are lazy (allocated on first access). Use for time-stepped conserved quantities (stage, xmomentum, ymomentum). 56 bytes per triangle eager; gradients/phi only allocated if accessed (e.g. by erosion operators).

  • 'edge_diagnostic' — centroid and edge values only. Use for static or derived diagnostic fields (elevation, height, velocity). 32 bytes per triangle. Gradients lazy.

  • 'centroid_only' — centroid values only. Use for scalar parameters that are never edge-interpolated (e.g. friction). 8 bytes per triangle.

  • 'coordinate' — centroid, edge, and eager vertex values. Use for mesh coordinate quantities (x, y). 56 bytes per triangle.

  If None (default), the type is looked up in domain._quantity_type_map[name], falling back to 'evolved' for unknown names.

Usage:

>>> Quantity(domain, name="newQ", register=True)
>>> Quantity(domain, name="tracer", register=True, qty_type='centroid_only')

If vertex_values are None Create array of zeros compatible with domain. Otherwise check that it is compatible with dimensions of domain. Otherwise raise an exception

For Quantities that need to be saved during checkpointing, set register=True. Registered Quantities can be found in the dictionary domain.quantities (note, other Quantities can exist).

Methods

__init__(domain[, vertex_values, name, ...])	Create Quantity object
backup_centroid_values()
bound_vertices_below_by_constant(bound)
bound_vertices_below_by_quantity(quantity)
compute_gradients()
compute_local_gradients()
extrapolate_first_order()	Extrapolate conserved quantities from centroid to vertices and edges for each volume using first order scheme.
extrapolate_second_order()
extrapolate_second_order_and_limit_by_edge()
extrapolate_second_order_and_limit_by_vertex()
get_beta()	Get default beta value for limiting
get_extremum_index([mode, indices])	Return index for maximum or minimum value of quantity (on centroids)
get_gradients()	Provide gradients.
get_integral([full_only, region, indices])	Compute the integral of quantity across entire domain, or over a region.
get_interpolated_values(interpolation_points)	Get values at interpolation points
get_maximum_index([indices])	See get extreme index for details
get_maximum_location([indices])	Return location of maximum value of quantity (on centroids)
get_maximum_value([indices])	Return maximum value of quantity (on centroids)
get_minimum_index([indices])	See get extreme index for details
get_minimum_location([indices])	Return location of minimum value of quantity (on centroids)
get_minimum_value([indices])	Return minimum value of quantity (on centroids)
get_name()
get_values([interpolation_points, location, ...])	Get values for quantity
get_vertex_values([xy, smooth, ...])	Return vertex values like an OBJ format i.e. one value per node.
interpolate()	Compute interpolated values at edges and centroid Pre-condition: vertex_values have been set
interpolate_from_edges_to_vertices()
interpolate_from_vertices_to_edges()
interpolate_old()	Compute interpolated values at edges and centroid Pre-condition: vertex_values have been set
limit()
limit_edges_by_all_neighbours()
limit_edges_by_neighbour()
limit_vertices_by_all_neighbours()
maximum(other)	Max of self with anything that could populate a quantity
minimum(other)	Max of self with anything that could populate a quantity
plot_quantity([filename, draw])
save_centroid_data_to_csv([filename])
save_data_to_dem([filename])
save_to_array([cellsize, NODATA_value, ...])	Interpolate quantity to an array
saxpy_centroid_values(a, b)
set_beta(beta)	Set default beta value for limiting
set_boundary_values([numeric])	Set boundary values
set_boundary_values_from_edges()	Set boundary values by simply extrapolating from the cells
set_name([name])
set_values([numeric, quantity, function, ...])	Set values for quantity based on different sources.
set_values_from_array(values[, location, ...])	Set values for quantity
set_values_from_constant(X, location, ...)	Set quantity values from specified constant X
set_values_from_file(filename, ...[, ...])	Set quantity based on arbitrary points in a points file using attribute_name selects name of attribute present in file.
set_values_from_function(f[, location, ...])	Set values for quantity using specified function.
set_values_from_geospatial_data(...[, ...])	Set values based on geo referenced geospatial data object.
set_values_from_lat_long_grid_file(filename)	Read Digital model from the following ASCII format (.asc or .grd)
set_values_from_points(points, values, ...)	Set quantity values from arbitray data points using fit_interpolate.fit
set_values_from_quantity(q, location, ...)	Set quantity values from specified quantity instance q
set_values_from_tif_file(filename[, ...])
set_values_from_utm_grid_file(filename[, ...])	Read Digital Elevation model from the following ASCII format (.asc, .grd or .dem)
set_values_from_utm_raster(raster[, ...])
set_vertex_values(A[, indices, use_cache, ...])	Set vertex values for all unique vertices based on input array A which has one entry per unique vertex, i.e. one value for each row in array self.domain.nodes.
smooth_vertex_values([use_cache, verbose])	Smooths vertex values.
update(timestep)

Attributes

counter
vertex_values

extrapolate_first_order()

Extrapolate conserved quantities from centroid to vertices and edges for each volume using first order scheme.

get_beta() → float

Get default beta value for limiting

get_extremum_index(mode: str | None = None, indices: list[int] | ndarray | None = None) → int

Return index for maximum or minimum value of quantity (on centroids)

Optional arguments:

mode is either ‘max’(default) or ‘min’. indices is the set of element ids that the operation applies to.

Usage:

i = Q.get_extreme_index()

Notes

We do not seek the extremum at vertices as each vertex can have multiple values - one for each triangle sharing it.

If there are multiple cells with same maximum value, the first cell encountered in the triangle array is returned.

get_gradients() → tuple[ndarray, ndarray]

Provide gradients. Use compute_gradients first.

get_integral(full_only: bool = True, region=None, indices: list[int] | ndarray | None = None) → float

Compute the integral of quantity across entire domain, or over a region. Eg my_region = anuga.Region(polygon = user_polygon)

Q.get_integral(region = my_region)

get_interpolated_values(interpolation_points: ArrayLike | GeospatialData, use_cache: bool = False, verbose: bool = False) → num.ndarray

Get values at interpolation points

The argument interpolation points must be given as either a list of absolute UTM coordinates or a geospatial data object.

get_maximum_index(indices: list[int] | ndarray | None = None) → int

See get extreme index for details

get_maximum_location(indices: list[int] | ndarray | None = None) → tuple[float, float]

Return location of maximum value of quantity (on centroids)

Optional argument:

indices is the set of element ids that the operation applies to.

Usage:

x, y = get_maximum_location()

Notes

We do not seek the maximum at vertices as each vertex can have multiple values - one for each triangle sharing it.

If there are multiple cells with same maximum value, the first cell encountered in the triangle array is returned.

get_maximum_value(indices: list[int] | ndarray | None = None) → float

Return maximum value of quantity (on centroids)

Optional argument:

indices is the set of element ids that the operation applies to.

Usage:

v = get_maximum_value()

Note, we do not seek the maximum at vertices as each vertex can have multiple values - one for each triangle sharing it

get_minimum_index(indices: list[int] | ndarray | None = None) → int

See get extreme index for details

get_minimum_location(indices: list[int] | ndarray | None = None) → tuple[float, float]

Return location of minimum value of quantity (on centroids)

Optional argument:

indices is the set of element ids that the operation applies to.

Usage:

x, y = get_minimum_location()

Notes

We do not seek the maximum at vertices as each vertex can have multiple values - one for each triangle sharing it.

If there are multiple cells with same maximum value, the first cell encountered in the triangle array is returned.

get_minimum_value(indices: list[int] | ndarray | None = None) → float

Return minimum value of quantity (on centroids)

Optional argument:

indices is the set of element ids that the operation applies to.

Usage:

v = get_minimum_value()

See get_maximum_value for more details.

get_values(interpolation_points: ArrayLike | GeospatialData | None = None, location: str = 'vertices', indices: list[int] | num.ndarray | None = None, use_cache: bool = False, verbose: bool = False) → num.ndarray

Get values for quantity

Extract values for quantity as a numeric array.

Inputs:

interpolation_points: List of x, y coordinates where value is

sought (using interpolation). If points are given, values of location and indices are ignored. Assume either absolute UTM coordinates or geospatial data object.

location: Where values are to be stored.

Permissible options are: vertices, edges, centroids and unique vertices. Default is ‘vertices’

The returned values will have the leading dimension equal to length of the indices list or N (all values) if indices is None.

In case of location == ‘centroids’ the dimension of returned values will be a list or a numerical array of length N, N being the number of elements.

In case of location == ‘vertices’ or ‘edges’ the dimension of returned values will be of dimension Nx3

In case of location == ‘unique vertices’ the average value at each vertex will be returned and the dimension of returned values will be a 1d array of length “number of vertices”

Indices is the set of element ids that the operation applies to. If indices is None (or omitted) all elements are returned as a copy of the relevant array. If performance is critical, use arrays domain.centroid_values, domain.vertex_values and domain.edge_values directly.

The values will be stored in elements following their internal ordering.

get_vertex_values(xy: bool = True, smooth: bool | None = None, centroid_averaging: bool | None = None, precision: type | None = None) → tuple[ndarray, ndarray] | tuple[ndarray, ndarray, ndarray, ndarray]

Return vertex values like an OBJ format i.e. one value per node.

The vertex values are returned as one sequence in the 1D float array A. If requested the coordinates will be returned in 1D arrays X and Y.

The connectivity is represented as an integer array, V, of dimension Mx3, where M is the number of triangles. Each row has three indices defining the triangle and they correspond to elements in the arrays X, Y and A.

If smooth is True, vertex values corresponding to one common coordinate set will be smoothed by taking the average of vertex values for each node. In this case vertex coordinates will be de-duplicated corresponding to the original nodes as obtained from the method general_mesh.get_nodes()

If no smoothings is required, vertex coordinates and values will be aggregated as a concatenation of values at vertices 0, vertices 1 and vertices 2. This corresponds to the node coordinates obtained from the method general_mesh.get_vertex_coordinates()

Calling convention:

if xy is True:
    X, Y, A, V = get_vertex_values()
else:
    A, V = get_vertex_values()

interpolate()

Compute interpolated values at edges and centroid Pre-condition: vertex_values have been set

interpolate_old()

Compute interpolated values at edges and centroid Pre-condition: vertex_values have been set

maximum(other)

Max of self with anything that could populate a quantity

E.g other can be a constant, an array, a function, another quantity (except for a filename or points, attributes (for now)) - see set_values for details

minimum(other)

Max of self with anything that could populate a quantity

E.g other can be a constant, an array, a function, another quantity (except for a filename or points, attributes (for now)) - see set_values for details

save_to_array(cellsize=None, NODATA_value=-9999.0, smooth=None, easting_min=None, easting_max=None, northing_min=None, northing_max=None, origin=None, verbose=False)

Interpolate quantity to an array

set_beta(beta: float) → None

Set default beta value for limiting

set_boundary_values(numeric: _Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | bool | int | float | complex | str | bytes | _NestedSequence[bool | int | float | complex | str | bytes] | Callable = 0.0) → None

Set boundary values

set_boundary_values_from_edges() → None

Set boundary values by simply extrapolating from the cells

set_values(numeric: ArrayLike | Callable | Quantity | float | str | None = None, quantity: Quantity | None = None, function: Callable | None = None, geospatial_data: GeospatialData | None = None, filename: str | None = None, raster: tuple | None = None, attribute_name: str | None = None, alpha: float | str | None = None, location: str = 'vertices', polygon: ArrayLike | None = None, indices: list[int] | num.ndarray | None = None, smooth: bool = False, verbose: bool = False, use_cache: bool = False) → None

Set values for quantity based on different sources.

numeric:

Compatible list, numeric array (see below) or constant. If callable it will treated as a function (see below) If instance of another Quantity it will be treated as such. If geo_spatial object it will be treated as such

quantity:

Another quantity (compatible quantity, e.g. obtained as a linear combination of quantities)

function:

Any callable object that takes two 1d arrays x and y each of length N and returns an array also of length N. The function will be evaluated at points determined by location and indices in the underlying mesh.

geospatial_data:

Arbitrary geo spatial dataset in the form of the class Geospatial_data. Mesh points are populated using fit_interpolate.fit fitting

filename:

Name of a points file (extension .pts, .csv, .txt or .xya) or dem file (ext .dem, .asc, .grd or .tif) containing data points and attributes for use with fit_interpolate.fit.

raster:

A class or a tuple (x,y,Z)

attribute_name:

If specified, any array matching that name will be used. from file or geospatial_data. Otherwise a default will be used.

alpha:

Smoothing parameter to be used with fit_interpolate.fit. See module fit_interpolate.fit for further details about alpha. Alpha will only be used with points, values or filename. Otherwise it will be ignored.

location: Where values are to be stored.

Permissible options are: vertices, edges, centroids Default is ‘vertices’

In case of location == ‘centroids’ the dimension values must be a list of a numerical array of length N, N being the number of elements. Otherwise it must be of dimension Nx3

The values will be stored in elements following their internal ordering.

If location is ‘unique vertices’ indices refers the set of node ids that the operation applies to. If location is not ‘unique vertices’ indices refers the set of triangle ids that the operation applies to.

If selected location is vertices, values for centroid and edges will be assigned interpolated values. In any other case, only values for the specified locations will be assigned and the others will be left undefined.

polygon: Restrict update of quantity to locations that fall

inside polygon. Polygon works by selecting indices and calling set_values recursively. Polygon mode has only been implemented for constant values so far.

indices: Restrict update of quantity to locations that are

identified by indices (e.g. node ids if location is ‘unique vertices’ or triangle ids otherwise).

verbose: True means that output to stdout is generated

use_cache: True means that caching of intermediate results is

attempted for fit_interpolate.fit.

Exactly one of the arguments numeric, quantity, function or filename must be present.

set_values_from_array(values, location='vertices', indices=None, use_cache=False, verbose=False)

Set values for quantity

values: numeric array location: Where values are to be stored. Permissible options are: vertices, centroid, unique vertices Default is ‘vertices’

indices - if this action is carried out on a subset of elements or unique vertices The element/unique vertex indices are specified here.

In case of location == ‘centroid’ the dimension values must be a list of a numerical array of length N, N being the number of elements.

Otherwise it must be of dimension Nx3

The values will be stored in elements following their internal ordering.

If selected location is vertices, values for centroid and edges will be assigned interpolated values. In any other case, only values for the specified locations will be assigned and the others will be left undefined.

set_values_from_constant(X, location, indices, verbose)

Set quantity values from specified constant X

set_values_from_file(filename, attribute_name, alpha, location, indices, verbose=False, use_cache=False, max_read_lines=None)

Set quantity based on arbitrary points in a points file using attribute_name selects name of attribute present in file. If attribute_name is not specified, use first available attribute as defined in geospatial_data.

set_values_from_function(f, location='vertices', indices=None, use_cache=False, verbose=False)

Set values for quantity using specified function.

Input:

f: x, y -> z Function where x, y and z are arrays
location: Where values are to be stored.
          Permissible options are: vertices, centroid,
          unique vertices
          Default is "vertices"
indices:

set_values_from_geospatial_data(geospatial_data, alpha, location, indices, verbose=False, use_cache=False)

Set values based on geo referenced geospatial data object.

set_values_from_lat_long_grid_file(filename, location='centroids', indices=None, verbose=False)

Read Digital model from the following ASCII format (.asc or .grd)

Example: ncols 41 nrows 41 xllcorner 140 yllcorner -30 cellsize 0.025 NODATA_value -9999 28.6 28.6 28.6 28.6 28.7 28.7 28.7 28.7 28.7 28.7 ….

This file would represent raster data from lower left corner at 140 long and -30 lat over to upper right corner 141 long -29 lat.

Here cellsize = 0.025 represents 1/40 of a degree

Parameters:

• filename (str) – name of input file in asc format

• location (str) – vertices or centroids, interpolation onto these locations

• indices – None or a list of indices where interploation occurs

• northern (bool) – Flag to specify northern or southern hemisphere

• verbose (bool) – level of printed feedback

set_values_from_points(points, values, alpha, location, indices, data_georef=None, verbose=False, use_cache=False)

Set quantity values from arbitray data points using fit_interpolate.fit

set_values_from_quantity(q, location, indices, verbose)

Set quantity values from specified quantity instance q

Location is ignored - vertices will always be used here.

set_values_from_utm_grid_file(filename, location='vertices', indices=None, verbose=False)

Read Digital Elevation model from the following ASCII format (.asc, .grd or .dem)

Example: ncols 3121 nrows 1800 xllcorner 722000 yllcorner 5893000 cellsize 25 NODATA_value -9999 138.3698 137.4194 136.5062 135.5558 ……….

An accompanying file with same basename but extension .prj must exist and is used to fix the UTM zone, datum, false northings and eastings.

The prj format is assumed to be as

Projection UTM Zone 56 Datum WGS84 Zunits NO Units METERS Spheroid WGS84 Xshift 0.0000000000 Yshift 10000000.0000000000 Parameters

set_vertex_values(A: _Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | bool | int | float | complex | str | bytes | _NestedSequence[bool | int | float | complex | str | bytes], indices: list[int] | ndarray | None = None, use_cache: bool = False, verbose: bool = False) → None

Set vertex values for all unique vertices based on input array A which has one entry per unique vertex, i.e. one value for each row in array self.domain.nodes.

indices is the list of vertex_id’s that will be set.

This function is used by set_values_from_array

smooth_vertex_values(use_cache: bool = False, verbose: bool = False) → None

Smooths vertex values.

Memory layout (qty_type)

Each Quantity allocates a subset of its internal arrays based on its type. The type is resolved in the following priority order:

1. The explicit qty_type keyword argument passed to Quantity.__init__.

2. domain._quantity_type_map[name] — the domain’s per-name override table.

3. 'evolved' — the default, which allocates every array for full backward compatibility.

The four built-in types and their memory footprint (assuming float64, 8 bytes per element, N triangles) are:

qty_type	Arrays allocated	Bytes / triangle
'evolved'	centroid, edge, explicit update, semi-implicit update, centroid backup. x/y-gradient and phi are lazy (allocated on first access).	56 N eager
'edge_diagnostic'	centroid, edge. x/y-gradient and phi are lazy.	32 N
'centroid_only'	centroid only	8 N
'coordinate'	centroid, edge, vertex (eager)	56 N

vertex_values is lazily allocated for all types except 'coordinate': the backing array is created on first access and is transparent to existing code. Gradient arrays (x_gradient, y_gradient) and phi are also lazy for all types — they are allocated on first access, which occurs only when explicitly called (e.g. by the erosion operators), not during a normal DE timestep.

Default types used by the shallow-water Domain:

Quantity	Type	Reason
stage, xmomentum, ymomentum	'evolved'	Time-stepped conserved quantities — need update arrays
elevation	'edge_diagnostic'	Static field; gradients lazy (allocated only if erosion operators call compute_local_gradients)
height, xvelocity, yvelocity	'edge_diagnostic'	Derived fields; edge values needed for output, no update machinery
friction	'centroid_only'	Scalar parameter; never edge-interpolated
x, y	'coordinate'	Mesh coordinates; vertex values required immediately at construction

For a 10-quantity shallow-water domain the optimised layout uses ~368 bytes/triangle versus 800 bytes/triangle if every quantity were 'evolved' — a saving of ~54 %.

Example — adding a custom tracer with minimal memory:

import anuga
from anuga.abstract_2d_finite_volumes.quantity import Quantity

domain = anuga.rectangular_cross_domain(100, 100)

# Option 1: pass qty_type directly
tracer = Quantity(domain, name='tracer', register=True,
                  qty_type='centroid_only')
tracer.set_values(0.0)

# Option 2: register via the domain's type map, then create with Quantity()
# Note: set_quantity() and add_quantity() operate on *existing* quantities
# only — they do not create a new one.  Use Quantity(..., register=True).
domain._quantity_type_map['salinity'] = 'edge_diagnostic'
salinity = Quantity(domain, name='salinity', register=True)
salinity.set_values(0.0)