anuga.Domain

class anuga.Domain(coordinates: _Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | bool | int | float | complex | str | bytes | _NestedSequence[bool | int | float | complex | str | bytes] | None = None, vertices: _Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | bool | int | float | complex | str | bytes | _NestedSequence[bool | int | float | complex | str | bytes] | None = None, boundary: dict | None = None, tagged_elements: dict | None = None, geo_reference=None, use_inscribed_circle: bool = False, mesh_filename: str | None = None, use_cache: bool = False, verbose: bool = False, conserved_quantities: list[str] | None = None, evolved_quantities: list[str] | None = None, other_quantities: list[str] | None = None, full_send_dict: dict | None = None, ghost_recv_dict: dict | None = None, starttime: float = 0, processor: int = 0, numproc: int = 1, number_of_full_nodes: int | None = None, number_of_full_triangles: int | None = None, ghost_layer_width: int = 2, **kwargs)[source]

Bases: Generic_Domain

Object which encapulates the shallow water model

This class is a specialization of class Generic_Domain from module generic_domain.py consisting of methods specific to the Shallow Water Wave Equation

Shallow Water Wave Equation

\[U_t + E_x + G_y = S\]

where

\[U = [w, uh, vh]^T\]
\[E = [uh, u^2h + gh^2/2, uvh]\]
\[G = [vh, uvh, v^2h + gh^2/2]\]

S represents source terms forcing the system (e.g. gravity, friction, wind stress, …)

and _t, _x, _y denote the derivative with respect to t, x and y respectively.

The quantities are

symbol

variable name

explanation

x

x

horizontal distance from origin [m]

y

y

vertical distance from origin [m]

z

elevation

elevation of bed on which flow is modelled [m]

h

height

water height above z [m]

w

stage

absolute water level, w = z+h [m]

u

speed in the x direction [m/s]

v

speed in the y direction [m/s]

uh

xmomentum

momentum in the x direction [m^2/s]

vh

ymomentum

momentum in the y direction [m^2/s]

eta

mannings friction coefficient [to appear]

nu

wind stress coefficient [to appear]

The conserved quantities are w, uh, vh

__init__(coordinates: _Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | bool | int | float | complex | str | bytes | _NestedSequence[bool | int | float | complex | str | bytes] | None = None, vertices: _Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | bool | int | float | complex | str | bytes | _NestedSequence[bool | int | float | complex | str | bytes] | None = None, boundary: dict | None = None, tagged_elements: dict | None = None, geo_reference=None, use_inscribed_circle: bool = False, mesh_filename: str | None = None, use_cache: bool = False, verbose: bool = False, conserved_quantities: list[str] | None = None, evolved_quantities: list[str] | None = None, other_quantities: list[str] | None = None, full_send_dict: dict | None = None, ghost_recv_dict: dict | None = None, starttime: float = 0, processor: int = 0, numproc: int = 1, number_of_full_nodes: int | None = None, number_of_full_triangles: int | None = None, ghost_layer_width: int = 2, **kwargs) → None[source]

Instantiate a shallow water domain.

Parameters:

  • coordinates – vertex locations for the mesh

  • vertices – vertex indices defining the triangles of the mesh

  • boundary – boundaries of the mesh

Methods

__init__([coordinates, vertices, boundary, ...])

Instantiate a shallow water domain.

add_quantity(name, *args, **kwargs)

Add values to a named quantity

apply_fractional_steps()

Override to sync GPU data before fractional step operators run.

apply_protection_against_isolated_degenerate_timesteps()

backup_conserved_quantities()

boundary_statistics([quantities, tags])

Output statistics about boundary forcing at each timestep

build_tagged_elements_dictionary(*args, **kwargs)

centroid_norm(quantity, normfunc)

Calculate the norm of the centroid values of a specific quantity, using normfunc.

check_integrity()

Run integrity checks on shallow water domain.

compute_boundary_flows()

Compute boundary flows at current timestep.

compute_fluxes()

Compute fluxes and timestep suitable for all volumes in domain.

compute_forcing_terms()

If there are any forcing functions driving the system they should be defined in Domain subclass and appended to the list self.forcing_terms

compute_total_volume()

Compute total volume (m^3) of water in entire domain

conserved_values_to_evolved_values(q_cons, ...)

Needs to be overridden by Domain subclass

create_quantity_from_expression(expression)

Create new quantity from other quantities using arbitrary expression.

distribute_edges_to_vertices()

Distribute edge values to vertices.

distribute_to_edges()

extrapolate centroid values edges

distribute_to_vertices_and_edges([...])

extrapolate centroid values to vertices and edges

dump_triangulation([filename])

Get vertex coordinates, partition full and ghost triangles based on self.tri_full_flag

evolve([yieldstep, outputstep, finaltime, ...])

Evolve method from Domain class.

evolve_one_euler_step(yieldstep, finaltime)

One Euler Time Step Q^{n+1} = E(h) Q^n

evolve_one_rk2_step(yieldstep, finaltime)

One 2nd order RK timestep Q^{n+1} = 0.5 Q^n + 0.5 E(h)^2 Q^n

evolve_one_rk3_step(yieldstep, finaltime)

One 3rd order RK timestep Q^(1) = 3/4 Q^n + 1/4 E(h)^2 Q^n (at time t^n + h/2) Q^{n+1} = 1/3 Q^n + 2/3 E(h) Q^(1) (at time t^{n+1})

evolve_to_end([finaltime])

Iterate evolve all the way to the end.

get_CFL()

get CFL

get_algorithm_parameters()

Get the standard parameter that are currently set (as a dictionary)

get_area(*args, **kwargs)

get_areas(*args, **kwargs)

get_beta()

Get default beta for limiting.

get_boundary_flux_integral()

Compute the boundary flux integral.

get_boundary_polygon(*args, **kwargs)

get_boundary_tags(*args, **kwargs)

get_centroid_coordinates(*args, **kwargs)

get_centroid_transmissive_bc()

Get value of centroid_transmissive_bc flag.

get_cfl()

get CFL

get_compute_fluxes_method()

Get method for computing fluxes.

get_conserved_quantities(vol_id[, vertex, edge])

Get conserved quantities at volume vol_id.

get_datadir()

get_datetime([timestamp])

Retrieve datetime corresponding to current timestamp wrt to domain timezone

get_disconnected_triangles(*args, **kwargs)

get_edge_midpoint_coordinate(*args, **kwargs)

get_edge_midpoint_coordinates(*args, **kwargs)

get_energy_through_cross_section(polyline[, ...])

Obtain average energy head [m] across specified cross section.

get_evolve_max_timestep()

Set default max_timestep for evolving.

get_evolve_min_timestep()

Set default max_timestep for evolving.

get_evolve_starttime()

get_evolved_quantities(vol_id[, vertex, edge])

Get evolved quantities at volume vol_id.

get_extent(*args, **kwargs)

get_flow_algorithm()

Get method used for timestepping and spatial discretisation

get_flow_through_cross_section(polyline[, ...])

Get the total flow through an arbitrary poly line.

get_fractional_step_volume_integral()

Compute the integrated flows from fractional steps.

get_full_centroid_coordinates(*args, **kwargs)

get_full_nodes(*args, **kwargs)

get_full_triangles(*args, **kwargs)

get_full_vertex_coordinates(*args, **kwargs)

get_georeference(*args, **kwargs)

get_global_name()

get_hemisphere()

get_interpolation_object(*args, **kwargs)

get_intersecting_segments(*args, **kwargs)

get_inv_tri_map()

get_lone_vertices(*args, **kwargs)

get_maximum_inundation_elevation([indices, ...])

Return highest elevation where h > 0

get_maximum_inundation_location([indices])

Return location of highest elevation where h > 0

get_minimum_allowed_height()

get_minimum_storable_height()

get_multiprocessor_mode()

Get multiprocessor mode

get_name()

get_nodes(*args, **kwargs)

get_normal(*args, **kwargs)

get_number_of_full_triangles(*args, **kwargs)

get_number_of_nodes(*args, **kwargs)

get_number_of_triangles(*args, **kwargs)

get_number_of_triangles_per_node(*args, **kwargs)

get_quantity(name[, location, indices])

Get pointer to quantity object.

get_quantity_names()

Get a list of all the quantity names that this domain is aware of.

get_radii(*args, **kwargs)

get_relative_time()

Set internal relative time

get_starttime([datetime])

return starttime, either as timestamp, or as a datetime

get_store()

Get whether data saved to sww file.

get_store_centroids()

Get whether data saved to sww file.

get_tagged_elements(*args, **kwargs)

get_time()

Get the absolute model time (seconds).

get_timestep()

get current timestep (seconds).

get_timestepping_method()

get_timezone()

Retrieve current domain timezone

get_tri_map()

get_triangle_containing_point(*args, **kwargs)

get_triangles(*args, **kwargs)

get_triangles_and_vertices_per_node(*args, ...)

get_triangles_inside_polygon(*args, **kwargs)

get_unique_vertices(*args, **kwargs)

get_using_centroid_averaging()

get_using_discontinuous_elevation()

Return boolean indicating whether algorithm is using dicontinuous elevation

get_vertex_coordinate(*args, **kwargs)

get_vertex_coordinates(*args, **kwargs)

get_water_volume()

get_wet_elements([indices, minimum_height])

Return indices for elements where h > minimum_allowed_height

get_zone()

get zone for domain Geo_reference

initialise_storage()

Create and initialise self.writer object for storing data.

log_operator_timestepping_statistics()

maximum_quantity(name, *args, **kwargs)

max of values to a named quantity

minimum_quantity(name, *args, **kwargs)

min of values to a named quantity

print_algorithm_parameters()

Print the standard parameters that are curently set (as a dictionary)

print_boundary_statistics([quantities, tags])

print_operator_statistics()

print_operator_timestepping_statistics()

print_statistics(*args, **kwargs)

print_timestepping_statistics(*args, **kwargs)

Print time stepping statistics.

print_volumetric_balance_statistics()

protect_against_infinitesimal_and_negative_heights()

Clean up the stage and momentum values to ensure non-negative heights

quantity_statistics([precision])

Return string with statistics about quantities for printing or logging

report_cells_with_small_local_timestep([...])

Convenience function to print the locations of cells with a small local timestep.

report_water_volume_statistics([verbose, ...])

Compute the volume, boundary flux integral, fractional step volume integral, and their difference

saxpy_conserved_quantities(a, b[, c])

set_CFL([cfl])

Set CFL parameter, warn if greater than 2.0

set_beta(beta)

Shorthand to assign one constant value [0,2] to all limiters.

set_betas(beta_w, beta_w_dry, beta_uh, ...)

Assign beta values in the range [0,2] to all limiters.

set_boundary(boundary_map)

Associate boundary objects with tagged segments (see base class).

set_centroid_transmissive_bc(flag)

Set behaviour of the transmissive boundary condition, namely calculate the BC using the centroid value of neighbouring cell or the calculated edge value.

set_collect_max_quantities([...])

Create (or return existing) Collect_max_quantities_operator on this domain.

set_cfl([cfl])

Set CFL parameter, warn if greater than 2.0

set_checkpointing([checkpoint, ...])

Set up checkpointing.

set_compute_fluxes_method([flag])

Set method for computing fluxes.

set_datadir(name)

set_default_order(n)

Set default (spatial) order to either 1 or 2.

set_evolve_max_timestep(max_timestep)

Set default max_timestep for evolving.

set_evolve_min_timestep(min_timestep)

Set default min_timestep for evolving.

set_evolve_starttime(time)

set_extrapolate_velocity([flag])

Extrapolation routine uses momentum by default, can change to velocity extrapolation which seems to work better.

set_fixed_flux_timestep([flux_timestep])

Disable variable timestepping and manually set a fixed flux_timestep

set_flow_algorithm([algorithm])

Set combination of slope limiting and time stepping

set_fractional_step_operator(operator)

set_georeference(*args, **kwargs)

set_gpu_interface()

set_hemisphere(hemisphere)

set_institution(institution)

set_low_froude([low_froude])

For low Froude problems the standard flux calculations can lead to excessive damping.

set_maximum_allowed_speed(maximum_allowed_speed)

Set the maximum particle speed that is allowed in water shallower than minimum_allowed_height.

set_minimum_allowed_height(...)

Set minimum depth that will be recognised in the numerical scheme.

set_minimum_storable_height(...)

Set the minimum depth that will be written to an SWW file.

set_multiprocessor_mode([multiprocessor_mode])

Set multiprocessor mode (legacy integer API).

set_name([name, timestamp])

Assign a name to this simulation.

set_omp_num_threads([omp_num_threads, verbose])

Set the OpenMP thread count (process-wide).

set_plotter(*args, **kwargs)

Set the plotter for this domain

set_points_file_block_line_size(...)

set_quantities_to_be_monitored(q[, polygon, ...])

Specify which quantities will be monitored for extrema.

set_quantities_to_be_stored(q)

Specify which quantities will be stored in the SWW file.

set_quantity(name, *args, **kwargs)

Set values for named quantity

set_quantity_vertices_dict(quantity_dict)

Set values for named quantities.

set_relative_time([time])

Set internal relative time

set_sloped_mannings_function([flag])

Set mannings friction function to use the sloped wetted area.

set_starttime([timestamp])

Set the starttime for the evolution

set_store([flag])

Set whether data saved to sww file.

set_store_centroids([flag])

Set whether centroid data is saved to sww file.

set_store_vertices_smoothly([flag, reduction])

Decide whether vertex values should be stored smoothly (one value per vertex) or uniquely as computed in the model (False).

set_store_vertices_uniquely([flag, reduction])

Decide whether vertex values should be stored uniquely as computed in the model (True) or whether they should be reduced to one value per vertex using self.reduction (False).

set_tag_region(*args, **kwargs)

Set quantities based on a regional tag.

set_time([time])

Set the model time (seconds).

set_timestepping_method(timestepping_method)

set_timezone([tz])

Set timezone for domain

set_use_kinematic_viscosity([flag])

set_use_optimise_dry_cells([flag])

Try to optimize calculations where region is dry

set_using_centroid_averaging([flag])

Set flag to use centroid averaging in output of smoothed vertex values.

set_using_discontinuous_elevation([flag])

Set flag to show whether compute flux algorithm is allowing discontinuous elevation.

set_zone(zone)

Set zone for domain.

statistics(*args, **kwargs)

store_timestep()

Store time dependent quantities and time.

sww_merge(*args, **kwargs)

Dummy function for sequential algorithms where the sww produced is the final products.

timestepping_statistics([track_speeds, ...])

Return string with time stepping statistics for printing or logging

tripcolor(*args, **kwargs)

triplot(*args, **kwargs)

update_boundary()

Go through list of boundary objects and update boundary values for all conserved quantities on boundary.

update_boundary_old()

Go through list of boundary objects and update boundary values for all conserved quantities on boundary.

update_boundary_old_2()

Go through list of boundary objects and update boundary values for all conserved quantities on boundary.

update_centroids_of_momentum_from_velocity()

Calculate the centroid value of x and y momentum from height and velocities.

update_centroids_of_velocities_and_height()

Calculate the centroid values of velocities and height based on the values of the quantities stage and x and y momentum

update_conserved_quantities()

Update vectors of conserved quantities using previously computed fluxes and specified forcing functions.

update_extrema()

Update extrema if requested by set_quantities_to_be_monitored.

update_ghosts([quantities])

Override to use GPU ghost exchange when in GPU mode.

update_other_quantities()

There may be a need to calculate some of the other quantities based on the new values of conserved quantities

update_special_conditions()

update_timestep(yieldstep, finaltime)

Calculate the next timestep to take

volumetric_balance_statistics()

Create volumetric balance report suitable for printing or logging.

write_boundary_statistics([quantities, tags])

write_time([track_speeds])

Attributes

flux_timestep

Current flux timestep [s]

g

Gravitational acceleration [m/s^2]

timestep

Current timestep [s]
apply_fractional_steps()[source]

Override to sync GPU data before fractional step operators run.

Boyd culvert operators are handled via GPUCulvertManager (batched, only 2 GPU sync points) instead of the per-operator Python loop.

check_integrity() → None[source]

Run integrity checks on shallow water domain.

compute_boundary_flows() → tuple[dict[str, float], float, float][source]

Compute boundary flows at current timestep.

Computes the total inflow and outflow across the domain boundary, as well as the flow across each tagged boundary segment.

Returns:

  • boundary_flows (dict) – Flow rates [m^3/s] for each boundary tag

  • total_boundary_inflow (float) – Total inflow across boundary [m^3/s]

  • total_boundary_outflow (float) – Total outflow across boundary [m^3/s]

Notes

These calculations are only approximate since they don’t use the flux calculation used in evolve. For exact computation, see get_boundary_flux_integral.

compute_capabilities() → dict[source]

Report which compute backends this build/run supports.

Returns a dict with:

'gpu_offload'     : bool - process can offload mode-2 to a GPU device
                           (build supports it, device present, offload
                           not disabled); see set_gpu_offload()
'num_gpu_devices' : int  - number of offload devices visible
'mpi'             : bool - gpu_ext built with C MPI ('unified' parallel ok)
'modes'           : list - per-domain modes available ('unified' only
                           when the gpu_ext extension is importable)
compute_fluxes()[source]

Compute fluxes and timestep suitable for all volumes in domain.

Compute total flux for each conserved quantity using “flux_function”

Fluxes across each edge are scaled by edgelengths and summed up Resulting flux is then scaled by area and stored in explicit_update for each of the three conserved quantities stage, xmomentum and ymomentum

The maximal allowable speed computed by the flux_function for each volume is converted to a timestep that must not be exceeded. The minimum of those is computed as the next overall timestep.

Post conditions:

domain.explicit_update is reset to computed flux values domain.flux_timestep is set to the largest step satisfying all volumes.

This wrapper calls the underlying C version of compute fluxes

compute_forcing_terms()[source]

If there are any forcing functions driving the system they should be defined in Domain subclass and appended to the list self.forcing_terms

compute_total_volume() → float[source]

Compute total volume (m^3) of water in entire domain

diagnose_timestep(threshold_dt: float | None = None, top_n: int = 5, threshold_depth: float | None = None, verbose: bool = True) → dict[source]

Diagnose CFL-limited timestep and NaN/inf in conserved quantities.

Computes the CFL-limiting dt per centroid (radius / wave_speed) and reports the cells driving the global minimum. Also scans stage, xmomentum and ymomentum for NaN / inf, which is the usual signature of an ADER-2 (or any) timestepping blow-up. Intended to be called from inside an evolve loop (e.g. each yieldstep) to localise where and when the solver is about to fail.

Parameters:

  • threshold_dt (float, optional) – If given, a WARNING is emitted whenever the global min local-dt falls below this value. Useful to catch the moment dt collapses.

  • top_n (int) – Number of worst (smallest-dt) cells per rank to include in the returned dict and printed report.

  • threshold_depth (float, optional) – Depths below this are clamped when computing speed, to avoid division-by-tiny. Defaults to self.minimum_allowed_height.

  • verbose (bool) – If True, print a per-rank report.

Returns:

Keys:

time              relative simulation time
current_timestep  self.timestep (last accepted global dt)
cfl               self.CFL setting
local_min_dt      smallest cell-local CFL dt on this rank
global_min_dt     smallest cell-local CFL dt across all ranks
local_max_speed   max wave_speed on this rank
global_max_speed  max wave_speed across all ranks
nan_counts        {'stage': n, 'xmomentum': n, 'ymomentum': n}
first_nan_cell    centroid xy of the first NaN/inf cell on
                  this rank, or None
worst_cells       list of dicts (top_n) with keys
                  {triangle, x, y, depth, speed, local_dt}

Return type:

dict

distribute_edges_to_vertices()[source]

Distribute edge values to vertices.

This is a wrapper for the C implementation of the distribution from edges to vertices.

distribute_to_edges()[source]

extrapolate centroid values edges

distribute_to_vertices_and_edges(distribute_to_vertices=True)[source]

extrapolate centroid values to vertices and edges

evolve(yieldstep: float | None = None, outputstep: float | None = None, finaltime: float | DateTime | None = None, duration: float | None = None, skip_initial_step: bool = False) → Iterator[float][source]

Evolve method from Domain class.

Parameters:

  • yieldstep (float, optional) – Yield every yieldstep time period

  • outputstep (float, optional) – Output to sww file every outputstep time period. outputstep should be an integer multiple of yieldstep.

  • finaltime (float or datetime, optional) – Evolve until finaltime (can be a float in seconds or a datetime object)

  • duration (float, optional) – Evolve for a time of length duration (seconds)

  • skip_initial_step (bool, optional) – Can be used to restart a simulation (not often used).

Notes

If outputstep is None, the output to sww file happens every yieldstep. If yieldstep is None then simply evolve to finaltime or for a duration.

evolve_one_ader2_step(yieldstep, finaltime)[source]

One ADER-2 timestep using the local Cauchy-Kovalewski predictor.

Q^{n+1} = Q^n + dt * R(Q^{n+1/2})

Uses the fused edge predictor: edge values are shifted to Q^{n+1/2} in-place while centroid values remain at Q^n, eliminating the second extrapolation pass and the backup/saxpy restore pattern.

Single-flux-call variant: the previous step’s CFL timestep is reused for the predictor half-advance. The first step bootstraps with dt=0 (Euler) to establish the initial CFL timestep.

Cost: 1 flux call + 1 extrapolation + 1 edge C-K predictor (after step 1). Accuracy: 2nd-order in space and time.

evolve_one_euler_step(yieldstep, finaltime)[source]

One Euler Time Step Q^{n+1} = E(h) Q^n

Does not assume that centroid values have been extrapolated to vertices and edges

evolve_one_rk2_step(yieldstep, finaltime)[source]

One 2nd order RK timestep Q^{n+1} = 0.5 Q^n + 0.5 E(h)^2 Q^n

Does not assume that centroid values have been extrapolated to vertices and edges

evolve_one_rk3_step(yieldstep, finaltime)[source]

One 3rd order RK timestep Q^(1) = 3/4 Q^n + 1/4 E(h)^2 Q^n (at time t^n + h/2) Q^{n+1} = 1/3 Q^n + 2/3 E(h) Q^(1) (at time t^{n+1})

Does not assume that centroid values have been extrapolated to vertices and edges

property flux_timestep: float

Current flux timestep [s]

property g: float

Gravitational acceleration [m/s^2]

get_algorithm_parameters() → dict[source]

Get the standard parameter that are currently set (as a dictionary)

get_boundary_flux_integral() → float[source]

Compute the boundary flux integral.

Should work in parallel

get_compute_fluxes_method() → str[source]

Get method for computing fluxes.

See set_compute_fluxes_method for possible choices.

get_compute_mode() → str[source]

Return the active per-domain compute mode: ‘legacy’ or ‘unified’.

get_datetime(timestamp: float | None = None) → DateTime[source]

Retrieve datetime corresponding to current timestamp wrt to domain timezone

param: timestamp: return datetime corresponding to given timestamp

get_energy_through_cross_section(polyline: _Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | bool | int | float | complex | str | bytes | _NestedSequence[bool | int | float | complex | str | bytes], kind: str = 'total', verbose: bool = False) → float[source]

Obtain average energy head [m] across specified cross section.

Inputs:
polyline: Representation of desired cross section - it may contain

multiple sections allowing for complex shapes. Assume absolute UTM coordinates. Format [[x0, y0], [x1, y1], …]

kind: Select which energy to compute.

Options are ‘specific’ and ‘total’ (default)

Output:

E: Average energy [m] across given segments for all stored times.

The average velocity is computed for each triangle intersected by the polyline and averaged weighted by segment lengths.

The typical usage of this function would be to get average energy of flow in a channel, and the polyline would then be a cross section perpendicular to the flow.

#FIXME (Ole) - need name for this energy reflecting that its dimension is [m].

get_flow_algorithm() → str[source]

Get method used for timestepping and spatial discretisation

get_flow_through_cross_section(polyline: _Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | bool | int | float | complex | str | bytes | _NestedSequence[bool | int | float | complex | str | bytes], verbose: bool = False) → float[source]

Get the total flow through an arbitrary poly line.

This is a run-time equivalent of the function with same name in sww_interrogate.py

Input:
polyline: Representation of desired cross section - it may contain

multiple sections allowing for complex shapes. Assume absolute UTM coordinates. Format [[x0, y0], [x1, y1], …]

Output:

Q: Total flow [m^3/s] across given segments.

get_fractional_step_volume_integral() → float[source]

Compute the integrated flows from fractional steps.

This requires that the fractional step operators update the fractional_step_volume_integral.

Should work in parallel

get_global_max_speed() → float[source]

Return maximum speed across all MPI ranks.

Usage:

max_speed = get_global_max_speed()

Note: This performs MPI reduction, so all ranks get the same result.

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

Return maximum stage value across all MPI ranks.

Optional argument:

indices: set of element ids that the operation applies to

Usage:

max_stage = get_global_max_stage()

Note: This performs MPI reduction, so all ranks get the same result.

get_global_wet_element_count(indices: list[int] | ndarray | None = None, minimum_height: float | None = None) → int[source]

Return total number of wet elements across all MPI ranks.

Optional arguments:

indices: set of element ids that the operation applies to minimum_height: threshold for considering an element wet

Usage:

count = get_global_wet_element_count()

Note: This performs MPI reduction, so all ranks get the same result.

get_maximum_inundation_elevation(indices: list[int] | ndarray | None = None, minimum_height: float | None = None) → float[source]

Return highest elevation where h > 0

Optional argument:

indices is the set of element ids that the operation applies to. minimum_height for testing h > minimum_height

Usage:

q = get_maximum_inundation_elevation()

Note, centroid values are used for this operation

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

Return location of highest elevation where h > 0

Optional argument:

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

Usage:

q = get_maximum_inundation_location()

Note, centroid values are used for this operation

get_multiprocessor_mode() → int[source]

Get multiprocessor mode

  1. openmp (in development)

  2. gpu/mpi (in development)

get_starttime(datetime: bool = False) → float | DateTime[source]

return starttime, either as timestamp, or as a datetime

get_store() → bool[source]

Get whether data saved to sww file.

get_store_centroids() → bool[source]

Get whether data saved to sww file.

get_timezone() → ZoneInfoType[source]

Retrieve current domain timezone

get_wet_elements(indices: list[int] | ndarray | None = None, minimum_height: float | None = None) → ndarray[source]

Return indices for elements where h > minimum_allowed_height

Optional argument:

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

Usage:

indices = get_wet_elements()

Note, centroid values are used for this operation

property gpu

Shortcut to access the GPU interface.

Returns the gpu_interface object which provides FLOP counters and kernel wrappers. Raises AttributeError if GPU mode is not enabled.

initialise_storage() → None[source]

Create and initialise self.writer object for storing data. Also, save x,y and bed elevation

load_balance_statistics(minimum_height: float | None = None) → dict[source]

Return load balance statistics for this domain (single-rank version).

For a parallel domain use Parallel_domain.load_balance_statistics() which gathers across all MPI ranks via Allgather. This serial version returns a dict with length-1 arrays so the interface is identical.

Parameters:

minimum_height (float, optional) – Depth threshold for “wet” classification. Defaults to anuga.config.minimum_allowed_height.

Returns:

Keys and shapes are the same as the parallel version:

n_full           int[1]   total triangle count
n_ghost          int[1]   0 (no ghost triangles in serial)
n_wet_full       int[1]   wet triangle count
wet_fraction     float[1] n_wet_full / n_full
ghost_fraction   float[1] 0.0
wall_time        float[1] total wall time since evolve() started
comm_time        float[1] 0.0
reduce_wait_time float[1] 0.0
compute_time     float[1] same as wall_time
imbalance_ratio  float    1.0
wet_compute_corr float    nan

Return type:

dict

print_algorithm_parameters() → None[source]

Print the standard parameters that are curently set (as a dictionary)

print_load_balance_statistics(minimum_height: float | None = None) → None[source]

Print a load balance summary to stdout.

For a single-process domain this just reports wet fraction and triangle count. The parallel override prints a per-rank table.

Parameters:

minimum_height (float, optional) – Passed through to load_balance_statistics().

print_timestepping_statistics(*args, **kwargs) → None[source]

Print time stepping statistics.

Parameters:

  • time_units (str, optional) – Time units for reporting. Options are ‘sec’, ‘min’, ‘hr’, ‘day’.

  • datetime (bool, optional) – Flag to use timestamp or datetime.

  • track_speed (bool, optional) – Optional boolean keyword that decides whether to report location of smallest timestep as well as a histogram and percentile report.

  • relative_time (bool, optional) – Flag to report relative time instead of absolute time.

  • triangle_id (int, optional) – Can be used to specify a particular triangle rather than the one with the largest speed.

protect_against_infinitesimal_and_negative_heights()[source]

Clean up the stage and momentum values to ensure non-negative heights

report_cells_with_small_local_timestep(threshold_depth: float | None = None) → None[source]

Convenience function to print the locations of cells with a small local timestep.

Computations are at cell centroids

Useful in models with complex meshes, to find ways to speed up the model

report_water_volume_statistics(verbose: bool = True, returnStats: bool = False) → list[float] | None[source]

Compute the volume, boundary flux integral, fractional step volume integral, and their difference

If verbose, print a summary If returnStats, return a list with the volume statistics

set_beta(beta: float) → None[source]

Shorthand to assign one constant value [0,2] to all limiters. 0 Corresponds to first order, where as larger values make use of the second order scheme.

set_betas(beta_w: float, beta_w_dry: float, beta_uh: float, beta_uh_dry: float, beta_vh: float, beta_vh_dry: float) → None[source]

Assign beta values in the range [0,2] to all limiters. 0 Corresponds to first order, where as larger values make use of the second order scheme.

set_boundary(boundary_map)[source]

Associate boundary objects with tagged segments (see base class).

Mode-2 (‘unified’) captures a device-side boundary classification and the per-edge Dirichlet/Time/File/… mappings when the GPU interface is first built. A later set_boundary() — e.g. switching a tag from Reflective to Dirichlet partway through a simulation — must invalidate those caches, otherwise the stale boundary keeps being applied on the device and the new condition is silently ignored.

set_checkpointing(checkpoint: bool = True, checkpoint_dir: str = 'CHECKPOINTS', checkpoint_step: int = 10, checkpoint_time: float | None = None) → None[source]

Set up checkpointing.

param checkpoint: Default = True. Set to False will turn off checkpointing param checkpoint_dir: Where to store checkpointing files param checkpoint_step: Save checkpoint files after this many yieldsteps param checkpoint_time: If set, over-rides checkpoint_step. save checkpoint files after this amount of walltime

set_collect_max_quantities(update_frequency=1, collection_start_time=0.0, velocity_zero_height=None, store_to_sww=True) → Collect_max_quantities_operator[source]

Create (or return existing) Collect_max_quantities_operator on this domain.

Tracks running maxima of stage, depth, speed, and momentum magnitude (||(uh, vh)||) over the simulation. Call once before domain.evolve().

Parameters:

  • update_frequency (int) – Update maxima every this many timesteps (default 1).

  • collection_start_time (float) – Only collect after this simulation time (default 0).

  • velocity_zero_height (float or None) – Zero velocity below this depth; defaults to minimum_allowed_height.

  • store_to_sww (bool) – If True (default), write running maxima to the SWW file every yield step as centroid quantities max_stage_c, max_depth_c, max_speed_c, max_uh_c.

Return type:

Collect_max_quantities_operator

set_compute_fluxes_method(flag: str = 'original') → None[source]

Set method for computing fluxes.

Currently

original wb_1 wb_2 wb_3 tsunami DE

set_compute_mode(mode: str = 'unified', verbose: bool = False) → None[source]

Select this domain’s compute mode (per-domain).

This is a per-domain setting — different domains in one script may use different modes. Whether ‘unified’ uses a GPU is a separate, process-wide decision (see set_gpu_offload()), because OpenMP target offload is a process-level runtime setting, not a per-domain one.

Under MPI, ‘unified’ requires a gpu_ext built with MPI; otherwise this falls back to ‘legacy’ (whose Python MPI exchange is correct in parallel) with a rank-0 warning. The active mode is recorded in self.compute_mode; the original request in self.requested_compute_mode.

Parameters:

mode ({'legacy', 'unified'}) –

  • 'legacy' — mode 1: the sw_domain_openmp_ext solver with serial-Python fractional-step operators.

  • 'unified' — mode 2: the unified sw_domain_gpu_ext C kernels (solver and operators). Runs CPU-multicore by default; offloads to a GPU only when GPU offload is enabled process-wide via anuga.set_gpu_offload() on a GPU-capable build.

set_extrapolate_velocity(flag: bool = True) → None[source]

Extrapolation routine uses momentum by default, can change to velocity extrapolation which seems to work better.

set_flow_algorithm(algorithm: str = 'DE0') → None[source]

Set combination of slope limiting and time stepping

Currently

DE0 DE1 DE2 DE0_7 DE1_7

set_low_froude(low_froude: int = 0) → None[source]

For low Froude problems the standard flux calculations can lead to excessive damping. Set low_froude to 1 or 2 for flux calculations which minimize the damping in this case.

set_maximum_allowed_speed(maximum_allowed_speed: float) → None[source]

Set the maximum particle speed that is allowed in water shallower than minimum_allowed_height.

maximum_allowed_speed

This is useful for controlling speeds in very thin layers of water and at the same time allow some movement avoiding pooling of water.

set_minimum_allowed_height(minimum_allowed_height: float) → None[source]

Set minimum depth that will be recognised in the numerical scheme.

minimum_allowed_height minimum allowed depth in meters

The parameter H0 (Minimal height for flux computation) is also set by this function.

set_minimum_storable_height(minimum_storable_height: float) → None[source]

Set the minimum depth that will be written to an SWW file.

minimum_storable_height minimum allowed SWW depth is in meters

This is useful for removing thin water layers that seems to be caused by friction creep.

set_multiprocessor_mode(multiprocessor_mode: int = 1) → None[source]

Set multiprocessor mode (legacy integer API).

  1. openmp - Python RK loop (use_c_rk_loop=False)

  2. gpu/mpi - C RK loop (use_c_rk_loop=True)

Thin wrapper over set_compute_mode(): 1 maps to 'legacy', 2 to 'unified'. Whether ‘unified’ offloads to a GPU is a separate, process-wide choice — see anuga.set_gpu_offload(). New code should prefer set_compute_mode().

set_omp_num_threads(omp_num_threads: int | None = None, verbose: bool = True) → None[source]

Set the OpenMP thread count (process-wide).

OpenMP thread count is a process-level setting, not per-domain. This is a thin wrapper that delegates to the module-level anuga.set_omp_num_threads(); prefer that in new code. Kept for backward compatibility, and records self.omp_num_threads.

set_plotter(*args, **kwargs)[source]

Set the plotter for this domain

set_points_file_block_line_size(points_file_block_line_size: int) → None[source]
set_quantities_to_be_stored(q: dict[str, int] | list[str] | None) → None[source]

Specify which quantities will be stored in the SWW file.

q must be either:

  • a dictionary with quantity names

  • a list of quantity names (for backwards compatibility)

  • None

The format of the dictionary is as follows

quantity_name: flag where flag must be either 1 or 2. If flag is 1, the quantity is considered static and will be stored once at the beginning of the simulation in a 1D array.

If flag is 2, the quantity is considered time dependent and it will be stored at each yieldstep by appending it to the appropriate 2D array in the sww file.

If q is None, storage will be switched off altogether.

Once the simulation has started and thw sww file opened, this function will have no effect.

The format, where q is a list of names is for backwards compatibility only. It will take the specified quantities to be time dependent and assume ‘elevation’ to be static regardless.

set_quantity(name: str, *args, **kwargs) → None[source]

Set values for named quantity

We have to do something special for ‘elevation’ otherwise pass through to generic set_quantity

Mode-2 (‘unified’): the device holds the authoritative centroid state once the GPU interface exists, so a host-only write has to be mirrored to it. That is handled one level down, in Quantity.set_values(), via the _notify_*_host_quantity_write() hooks below — which also covers callers that reach a Quantity directly and bypass this method.

set_sloped_mannings_function(flag: bool = True) → None[source]

Set mannings friction function to use the sloped wetted area.

The flag is tested in the python wrapper mannings_friction_implicit

set_starttime(timestamp: float | DateTime = 0.0) → None[source]

Set the starttime for the evolution

Parameters:

timestamp – Either a float or a datetime object

Essentially we use unix time as our absolute time. So time = 0 corresponds to Jan 1st 1970 UTC

Use naive datetime which will be localized to the domain timezone or or use zoneinfo.ZoneInfo to set the timezone of datetime. Don’t use the tzinfo argument of datetime to set timezone as this does not work!

Example

Without setting timezone for the domain and the starttime then time calculations are all based on UTC. Note the timestamp, which is time in seconds from 1st Jan 1970 UTC.

>>> import anuga
>>> from zoneinfo import ZoneInfo
>>> from datetime import datetime
>>>
>>> domain = anuga.rectangular_cross_domain(10,10)
>>> dt = datetime(2021,3,21,18,30)
>>> domain.set_starttime(dt)
>>> print(domain.get_datetime(), 'TZ', domain.get_timezone(), 'Timestamp: ', domain.get_time())
2021-03-21 18:30:00+00:00 TZ UTC Timestamp:  1616351400.0

Example

Setting timezone for the domain, then naive datetime will be localizes to the domain timezone. Note the timestamp, which is time in seconds from 1st Jan 1970 UTC.

>>> import anuga
>>> from zoneinfo import ZoneInfo
>>> from datetime import datetime
>>>
>>> domain = anuga.rectangular_cross_domain(10,10)
>>> AEST = ZoneInfo('Australia/Sydney')
>>> domain.set_timezone(AEST)
>>>
>>> dt = datetime(2021,3,21,18,30)
>>> domain.set_starttime(dt)
>>> print(domain.get_datetime(), 'TZ', domain.get_timezone(), 'Timestamp: ', domain.get_time())
2021-03-21 18:30:00+11:00 TZ Australia/Sydney Timestamp:  1616311800.0

Example

Setting timezone for the domain, and setting the timezone for the datetime. Note the timestamp, which is time in seconds from 1st Jan 1970 UTC is the same as the previous example.

>>> import anuga
>>> from zoneinfo import ZoneInfo
>>> from datetime import datetime
>>>
>>> domain = anuga.rectangular_cross_domain(10,10)
>>>
>>> ACST = ZoneInfo('Australia/Adelaide')
>>> domain.set_timezone(ACST)
>>>
>>> AEST = ZoneInfo('Australia/Sydney')
>>> dt = datetime(2021,3,21,18,30, tzinfo=AEST)
>>>
>>> domain.set_starttime(dt)
>>> print(domain.get_datetime(), 'TZ', domain.get_timezone(), 'Timestamp: ', domain.get_time())
2021-03-21 18:00:00+10:30 TZ Australia/Adelaide Timestamp:  1616311800.0
set_store(flag: bool = True) → None[source]

Set whether data saved to sww file.

set_store_centroids(flag: bool = True) → None[source]

Set whether centroid data is saved to sww file.

set_store_vertices_smoothly(flag: bool = True, reduction: Callable | None = None) → None[source]

Decide whether vertex values should be stored smoothly (one value per vertex) or uniquely as computed in the model (False).

set_store_vertices_uniquely(flag: bool = True, reduction: Callable | None = None) → None[source]

Decide whether vertex values should be stored uniquely as computed in the model (True) or whether they should be reduced to one value per vertex using self.reduction (False).

set_timezone(tz: str | ZoneInfoType | None = None) → None[source]

Set timezone for domain

Parameters:

tz – either a timezone object or string

We recommend using the ZoneInfo provided by zoneinfo. Default is ZoneInfo(‘UTC’)

Example: Set default timezone UTC

>>> domain.set_timezone()

Example: Set timezone using tsdata string

>>> domain.set_timezone('Australia/Syndey')

Example: Set timezone using ZoneInfo timezone

>>> from zoneinfo import ZoneInfo
>>> new_tz = ZoneInfo('Australia/Sydney')
>>> domain.set_timezone(new_tz)
set_use_optimise_dry_cells(flag: bool = True) → None[source]

Try to optimize calculations where region is dry

store_timestep() → None[source]

Store time dependent quantities and time.

Precondition:

self.writer has been initialised

sww_merge(*args, **kwargs) → None[source]

Dummy function for sequential algorithms where the sww produced is the final products.

For parallel runs, a similarly named routine in parallel_shallow_water will merge all the sub domain sww files into a global sww file

Parameters:

  • verbose (bool) – Flag to produce more output

  • delete_old (bool) – Flag to delete sub domain sww files after creating global sww file

property timestep: float

Current timestep [s]

timestepping_statistics(track_speeds: bool = False, triangle_id: int | None = None, relative_time: bool = False, time_unit: str = 'sec', datetime: bool = False) → str[source]

Return string with time stepping statistics for printing or logging

Parameters:

  • time_units (str, optional) – Time units for reporting. Options are ‘sec’, ‘min’, ‘hr’, ‘day’.

  • datetime (bool, optional) – Flag to use timestamp or datetime.

  • track_speeds (bool, optional) – Optional boolean keyword that decides whether to report location of smallest timestep as well as a histogram and percentile report.

  • relative_time (bool, optional) – Flag to report relative time instead of absolute time.

  • triangle_id (int, optional) – Can be used to specify a particular triangle rather than the one with the largest speed.

Returns:

Formatted string with time stepping statistics.

Return type:

str

update_boundary()[source]

Go through list of boundary objects and update boundary values for all conserved quantities on boundary. It is assumed that the ordering of conserved quantities is consistent between the domain and the boundary object, i.e. the jth element of vector q must correspond to the jth conserved quantity in domain.

update_centroids_of_momentum_from_velocity()[source]

Calculate the centroid value of x and y momentum from height and velocities.

This method computes the centroid values of x and y momentum (xmomentum and ymomentum) by multiplying the centroid velocities by the centroid height values. The method assumes that the centroids of height and velocities are already up to date.

This is particularly useful for kinematic viscosity calculations where momentum values at cell centroids are required.

The method updates: - xmomentum.centroid_values: product of xvelocity and height at centroids - ymomentum.centroid_values: product of yvelocity and height at centroids

After updating centroid values, the method distributes these values to vertices and edges via distribute_to_vertices_and_edges().

Notes

This method modifies the centroid_values arrays in-place for both xmomentum and ymomentum quantities.

See also

distribute_to_vertices_and_edges

Distribute centroid values to vertices and edges

update_centroids_of_velocities_and_height()[source]

Calculate the centroid values of velocities and height based on the values of the quantities stage and x and y momentum

Assumes that stage and momentum are up to date

Useful for kinematic viscosity calculations

update_conserved_quantities()[source]

Update vectors of conserved quantities using previously computed fluxes and specified forcing functions.

update_domain_c_struct()[source]

Update the C domain structure from the Python Domain object.

update_ghosts(quantities=None)[source]

Override to use GPU ghost exchange when in GPU mode.

update_other_quantities()[source]

There may be a need to calculate some of the other quantities based on the new values of conserved quantities

update_timestep(yieldstep, finaltime)[source]

Calculate the next timestep to take

property use_c_rk2_loop

use use_c_rk_loop instead.

Type:

Deprecated

volumetric_balance_statistics() → str[source]

Create volumetric balance report suitable for printing or logging.