anuga.Domain

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

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]

BEDLOAD_PARAMETER_SETS = {'wong_parker_eq23': (4.93, 1.6, 0.047), 'wong_parker_eq24': (3.97, 1.5, 0.0495)}: Bedload parameter sets for [K-1], keyed by name. Which of the two Wong & Parker relations FG21 used is still open, so both are provided and neither is privileged beyond the default.

add_sediment_fraction(name, diameter, rho_s=2650.0, tau_c_star=0.04, d_star=1.0, beta=None, initial_concentration=0.0, reference_height=None, **settling_kwargs)[source]

Register one suspended sediment fraction and return its index.

Everything here is a property of THIS fraction. The domain-wide parameters live on initialize_sediment_operator(); in particular there is no rho_w here, because there is one fluid.

domain.add_sediment_fraction('sand', diameter=2.0e-4)
domain.add_sediment_fraction('silt', diameter=2.0e-5, tau_c_star=0.11)

A sediment fraction is a tracer with settling parameters attached, so it inherits the transport, boundary and conservation machinery described under Passive Tracers, and takes the tracer slot of the same index. Do not interleave add_tracer() and add_sediment_fraction on the same domain if you rely on that correspondence.

If the domain has no sediment operator yet, one is created with the default domain-wide parameters.

Parameters:

name (str) – Identifier for this fraction, e.g. ‘sand’. Also its tracer name.
diameter (float) – Grain diameter in metres.
rho_s (float, optional) – Sediment particle density, kg/m3. Default 2650 (quartz). Enters as the submerged specific gravity R = rho_s/rho_w - 1.
tau_c_star (float, optional) – Critical Shields stress for entrainment [E-1]. Default 0.04. Setting it to 0 disables entrainment for this fraction, leaving deposition only.
d_star (float, optional) – Ratio of near-bed to depth-averaged concentration in [D-1]. Default 1.0, the well-mixed limit. Ignored when the domain’s near_bed mode is ‘rouse’, which computes it per cell.
beta (float, optional) – Edge reconstruction limiter for this fraction’s tracer.
initial_concentration (float or array-like, optional) – Initial c_s; seeds m = h*c consistently.
reference_height (float, optional) – a in [S-4], in metres. Defaults to 2*diameter.
**settling_kwargs – Passed to settling_velocity() – e.g. shape=’natural’.

Returns:

The index of this fraction, which is also its tracer index.

Return type:

int

See also

initialize_sediment_operator: the domain-wide parameters.

add_tracer(name, beta=None, initial_value=0.0)[source]

Register a passive tracer named name and return its index.

A tracer is a depth-averaged concentration c advected with the water flux. The conserved variable is m = h*c; c is derived from it each substep, exactly as height is derived from stage.

Parameters:

name (str) – Identifier for the tracer. Must be unique on this domain.
beta (float, optional) – Edge-reconstruction limiter coefficient. 0 selects first order, > 0 a limited second-order reconstruction. This is a single value shared by every tracer (the C struct carries one beta_tracer scalar), so passing a value that disagrees with an already-registered tracer is an error rather than a silent last-writer-wins. Defaults to leaving the current value alone.
initial_value (float or array-like, optional) – Initial concentration c. A scalar fills the domain; an array must have one value per cell. m = h*c is seeded consistently.

Returns:

The tracer’s index, i.e. its row in the (ns, …) arrays.

Return type:

int

Notes

Registering a tracer reallocates every tracer array, so any reference held to one of them beforehand becomes stale. Add every tracer before seeding values, or re-fetch via get_tracer.

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.

Operator ORDER matters: fractional-step operators mutate stage in sequence, so running them in a different order gives different answers. Mode 1 runs them in registration order. Mode 2 must too, or the same script silently disagrees with itself across compute modes — and mode-1-vs-mode-2 agreement is the oracle we use to validate GPU work, so it must not depend on the order a user happened to construct their operators in.

The batch is therefore fired at the position of the first culvert in the list, not before the loop. When the culverts are registered together — which is what every real script does (towradgi registers 22 in a row) — that is exactly registration order, and the batching is untouched, so this costs nothing.

If culverts are interleaved with other operators, the later ones are still pulled forward to the first culvert’s slot (the C batches all registered culverts in one gather/compute/scatter cycle, and there is no unbatched mode-2 culvert path). We warn in that case rather than diverge silently. See issue #192.

check_integrity() → None[source]: Run integrity checks on shallow water domain.

check_tracer_conservation(name)[source]

Return (change_in_mass, boundary_flux_integral, discrepancy).

For a domain with no tracer source terms the first two should agree, so the discrepancy is the conservation error. On a CLOSED domain the flux integral is zero and this reduces to “did the mass stay constant”.

Returns absolute quantities, not a relative error: for a tracer that is zero almost everywhere a relative measure is meaningless, and the caller knows the scale that matters.

COLLECTIVE in parallel: every rank must call it.

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

erodible_thickness()[source]

Remaining erodible thickness per centroid [m], or None if no base.

elevation - sediment_z_base, clipped at zero. Zero means the cell has reached bedrock and will not erode further.

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_active_set_stats()[source]

Return (mean_active_fraction, n_rebuilds) once evolution has run.

mean_active_fraction is 1.0 whenever the active set never engaged (flag off, non-GPU compute mode, MPI fallback, or before the first rebuild).

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

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

openmp (in development)
gpu/mpi (in development)

get_sediment_names()[source]: Return the registered sediment class names, in index order.

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_tracer(name)[source]: Return the concentration c of tracer name, one value per cell.

get_tracer_boundary(name, tag)[source]: Return the boundary concentrations of name on tag.

get_tracer_boundary_flux_integral(name)[source]

Net tracer mass that has crossed the domain boundary, since t=0.

Positive is INTO the domain, matching the sign of the water balance. The tracer analogue of get_boundary_flux_integral, and the other half of a mass budget:

mass(t) - mass(0) == boundary_flux_integral(t)

to within the timestepping error, for a domain with no source terms. check_tracer_conservation does exactly that comparison.

The kernel accumulates the per-edge boundary flux into a per-cell array each substep; tracer_flux_integral_operator applies the timestepping method’s own weights and zeroes it, exactly as boundary_flux_integral_operator does for water.

Reduced across ranks, matching get_boundary_flux_integral for water: each rank counts only the boundary edges of the cells it owns, so the per-rank integrals sum to the whole-domain figure. COLLECTIVE.

get_tracer_index(name)[source]: Return the row index of tracer name.

get_tracer_mass(name)[source]

Total mass of tracer name in the domain: the integral of m = h*c.

The tracer analogue of get_water_volume. Summed over owned cells only and reduced across ranks, so it is the whole-domain figure in parallel as well as in serial.

COLLECTIVE in parallel, like get_water_volume: every rank must call it, or the ones that do will block.

get_tracer_names()[source]: Return the registered tracer names, in index order.

get_use_active_set() → bool[source]: Whether active-set stepping is requested (see set_use_active_set).

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

initialize_sediment_operator(porosity=None, c_max=None, c_pack=None, bed_evolution=None, rho_w=None, description=None, label=None, logging=False, verbose=False)[source]

Switch sediment transport on, and return the operator.

This is the entry point for sediment transport. It takes the DOMAIN-WIDE parameters – the ones that describe the run as a whole rather than any one fraction – and creates the single Sediment_transport_operator that carries the bed exchange.

Grain sizes are added separately, with add_sediment_fraction():

domain.initialize_sediment_operator(porosity=0.28, rho_w=1000.0)
domain.add_sediment_fraction('sand', diameter=2.0e-4)
domain.add_sediment_fraction('silt', diameter=2.0e-5)

The split is the point. A parameter belongs to exactly one of the two calls, so there is never a question of which call wins: porosity and rho_w are properties of the run, diameter and rho_s are properties of a fraction.

ONE OPERATOR PER DOMAIN. Calling this twice returns the same operator, applying any parameters given the second time – the kernel makes a single pass over every registered fraction, so a second operator in the fractional-step list would apply the bed exchange twice per step.

Calling it is optional in the simplest case: add_sediment_fraction() will create the operator with default domain-wide parameters if none exists. The two may be called in either order.

The closure choices – which shear, erosion, deposition and bedload laws to use – have their own setters, because each carries its own parameters and validation: set_shear_closure(), set_sediment_friction(), set_bed_material(), set_deposition() and set_bedload().

Parameters:

porosity (float, optional) – Bed porosity lambda in [G-4]. Default 0.30.
c_max (float, optional) – [L-2], the ceiling on depth-averaged concentration. Default 0.30.
c_pack (float, optional) – [L-4], maximum packing bounding near-bed concentration.
bed_evolution (bool, optional) – True (default) evolves the bed; False is the fixed-bed stage.
rho_w (float, optional) – Water density, kg/m3. Default 1000.
description – Passed to the operator; see Operator.
label – Passed to the operator; see Operator.
logging – Passed to the operator; see Operator.
verbose – Passed to the operator; see Operator.

Returns:

The domain’s sediment operator.

Return type:

Sediment_transport_operator

See also

add_sediment_fraction: register one sediment fraction.
set_sediment_parameters: change the domain-wide parameters later.
sediment_summary: print the complete active configuration.

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

property omp_num_threads: int

The process-wide OpenMP thread count (read-only view).

OpenMP thread count is a process-level setting, not per-domain, so this reflects the live value set by anuga.set_omp_num_threads() for every domain in the session — including domains constructed before the call. Assigning to it (domain.omp_num_threads = n) is kept for backward compatibility and sets the count process-wide.

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

sediment_summary()[source]

Return the complete active sediment configuration as text.

Every choice that affects the answer, with its units and the spec label it implements. Worth printing at the top of any run: the module has enough switches that “which erosion law was that?” is a real question six months later, and several of the choices are physics statements rather than tuning (spec 4.1.1).

set_angle_of_repose(angle=None, relax=1.0, max_sweeps=50)[source]

Relax bed slopes steeper than angle by moving material downslope.

Spec 7, from FG21 §2.2.4. Where the centroid-to-centroid bed slope exceeds the critical angle, material is diffused downslope until it does not.

domain.set_angle_of_repose(35.0) # degrees; FG21 use 35 domain.set_angle_of_repose(None) # off again (the default)

Parameters:

angle (float) – Critical angle in DEGREES, in (0, 90). None or 0 disables it.
relax (float) – Relaxation in (0, 1], default 1.0. The kernel already divides by the edge count for stability, so 1.0 is the fastest STABLE setting rather than an aggressive one – measured 793 sweeps to converge an over-steep cone against 2400+ at 0.3. Lower it only if you see something pathological.
max_sweeps (int) – Hard cap on sweeps per timestep, default 50. Reaching it is reported, because it means the bed may still be over-steep.

Notes

FG21 are explicit that this is a numerical heuristic, not physics: real bed slope failures are advective. It exists to stop the rest of the model breaking on over-steep slopes. It has a side effect worth knowing before you switch it on – it limits the steepness of canyon walls and knickpoints, and so suppresses knickpoint retreat that may be real. That is why it is off by default.

Mass is conserved: material removed from an over-steep cell is deposited on its neighbour, never discarded. This is the one place this module differs sharply from sanddune_erosion_operator, which lowers an over-steep cell and lets the material vanish.

Respects [L-5]: a cell cannot slump away material it is not allowed to lose, so a locked cell or one at its base stays put and its neighbours relax around it.

On the sweep count. This is an explicit diffusion solve, so convergence from a badly over-steep bed is slow: an over-steep cone needed 793 sweeps to reach the critical angle from cold. That is not what the per-timestep cap is sized for. In a running model the bed is already near-relaxed and each step needs a handful of sweeps; the cap is there for the pathological case, and hitting it is not fatal – progress carries over, so the bed keeps relaxing on subsequent steps. It is reported so that you know relaxation is lagging rather than finished.

If you START from a bed steeper than the critical angle, expect the cap to be hit on the first steps. Either accept that it settles over the first few, or raise max_sweeps for that run.

set_bed_material(material='noncohesive', tau_crit=0.088, K_e=None)[source]

Select the erosion law by naming the BED MATERIAL (spec 4.1.1).

‘noncohesive’ (default) – sand, gravel, boulders. Shields entrainment via Smith & McLean / Parker, [E-1]/[E-2], with a critical Shields stress per fraction (tau_c_star on add_sediment_fraction).

‘partheniades’ – [E-4], E = K_p (tau_b - tau_c)/tau_c, the form RDycore-sediment uses. K_e here is the Partheniades coefficient as a mass flux in kg m-2 s-1 (RDy26 use 1e-4), a different quantity from [E-5]’s erodibility, and is divided internally by the class density.

‘cohesive’ – silt, clay, cohesive bank material. Excess DIMENSIONAL shear via Hanson & Simon [E-3], E = K_e (tau_b - tau_c), with the jet-test erodibility [E-5] K_e = 0.2e-6 / sqrt(tau_c) unless K_e is given explicitly. tau_crit is in pascals, not Shields units; aS16 use 0.088.

Notes

[E-1] and [E-3] are not competing formulations of the same physics – they describe different sediment, calibrated from different experiments. Spec 4.1.1 puts it plainly: choosing between them is a statement about the bed, and getting it wrong is a physics error, not a tuning error. That is why this method is named for the material rather than for the equation.

set_bedload(formula='wong_parker_eq24', K=None, m=None, tau_c_star=None)[source]

Enable bedload transport [K-1]-[K-4] and its bed evolution [G-5].

Parameters:

formula (str) – ‘wong_parker_eq24’ (default; K`=3.97, `m`=1.5, `tau_c*`=0.0495), `’wong_parker_eq23’ (4.93, 1.60, 0.0470), ‘engelund_hansen’, or ‘off’.
K (float, optional) – Override the chosen set’s values individually.
m (float, optional) – Override the chosen set’s values individually.
tau_c_star (float, optional) – Override the chosen set’s values individually.

Notes

Engelund & Hansen is a TOTAL LOAD relation [K-5]: suspension is already inside it, so running the suspended source alongside it double counts. Spec 6 calls this out as a critical usage rule, and this method enforces it rather than warning – selecting ‘engelund_hansen’ turns the suspended exchange off, and it has no threshold by construction.

Which Wong & Parker relation FG21 used, their Eq 23 or Eq 24, is an open item. The two differ enough to matter, so the parameters are exposed: resolving it is a change of default, not an edit.

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_deposition(law='d_star', tau_d=0.0, near_bed='constant', reference_height_floor=0.01)[source]

Select the deposition law and its near-bed treatment (spec 4.4).

Parameters:

law ({'d_star', 'threshold'}) – ‘d_star’ (default) – [D-1], D = d*(Z) c v_s. ‘threshold’ – [D-2], D = v_s c (1 - tau_b/tau_d) for tau_b < tau_d, else zero; RDycore-sediment’s form.
tau_d (float) – Critical deposition stress in pascals, for ‘threshold’. tau_d = 0 disables deposition entirely – the hook RDycore’s passive-transport benchmarks rely on.
near_bed ({'constant', 'rouse'}) – How d* in [D-1] is obtained. ‘constant’ uses the per-fraction value given to add_sediment_fraction (default 1.0, the well-mixed limit of P14/P13). ‘rouse’ evaluates the fitted [S-4] profile per cell from the local Rouse number.
reference_height_floor (float) – The van Rijn-style floor a >= floor * h, used only by ‘rouse’. Default 0.01. Set to 0 to reach anugaSed’s regime, which applies no floor – see spec 12, D4b, where this is the largest single divergence from them.

set_erodible_base(elevation=None, depth=None)[source]

Set the non-erodible base – bedrock – below which no erosion acts.

[L-5]. Without this the bed is bottomless: erosion lowers it for as long as the flow has the strength to, which is the right default for an alluvial channel and wrong wherever the erodible layer is finite – a scoured reach over an outcrop, a lined culvert, a dam apron, a soil layer of known depth over rock.

The base is a per-CENTROID elevation, not a scalar, because bedrock is a surface. Give it either way round:

domain.set_erodible_base(elevation=z_rock)  # absolute [m]
domain.set_erodible_base(depth=0.5)         # 0.5 m below the
                                            # elevation set so far
domain.set_erodible_base()                  # remove the base

Parameters:

elevation (float or array (n,)) – Base elevation, in the same datum as the domain’s elevation quantity. Scalar broadcasts.
depth (float or array (n,)) – Erodible thickness below the CURRENT bed. The base is recorded as an elevation at the moment of the call, so later changes to the elevation quantity do not move it. Scalar broadcasts.
neither (Give exactly one. With)
is (the base is removed and the bed)
again. (bottomless)

Notes

The limiter acts on the SOURCE, never by clamping elevation. Erosion is scaled back to what the remaining thickness can supply, and the sediment that is not eroded never enters the water column, so the budget still closes exactly. Clamping z afterwards would leave suspended sediment that came from nowhere.

Where several classes compete for the last of the material, they are scaled by a shared proportional factor rather than served in order: the bed carries no per-class stratigraphy, so no class has a better claim, and registration order must not change the answer. Deposition is never scaled – it is what replenishes the bed.

Bedload is limited too, and stays exactly conservative while it is: the limit applies to the transport vector and to whole edges, both of which the two cells sharing an edge see identically.

set_erodible_region(region=None, polygon=None, center=None, radius=None, indices=None, erodible=True)[source]

Restrict erosion to part of the domain (or lock part of it).

The bed is erodible everywhere by default. Give a region to say otherwise:

domain.set_erodible_region(polygon=breach)     # ONLY here erodes
domain.set_erodible_region(polygon=apron,      # everywhere BUT here
                           erodible=False)
domain.set_erodible_region(my_region)          # a Region object
domain.set_erodible_region()                   # remove the restriction

Parameters:

region (Region) – An already-built anuga.abstract_2d_finite_volumes.region.Region. This is the general form: Region also understands line=, poly= and expand_polygon=, which have no keyword here, so build one and pass it when you need them. It must belong to THIS domain – its indices mean nothing on another mesh, and one built elsewhere is rejected rather than quietly mis-selecting.
polygon (list of [x, y]) – Region boundary, as for the region-based operators.
center ([x, y], float) – A circular region instead.
radius ([x, y], float) – A circular region instead.
indices (array of int) – Triangle ids directly. Overrides the geometric arguments.
erodible (bool) – True (default): the region named is the ONLY erodible part. False: the region named is the only LOCKED part.
the (The keyword arguments are the same ones Erosion_operator and)
take (other region-based operators)
same (and are resolved by the)
class (Region)
selects (so a polygon that selects a set of cells there)
here. (the same set)

Notes

A locked cell is held at the elevation it has WHEN THIS IS CALLED, by giving it zero erodible thickness – the region restriction is [L-5] with the layer set to nothing, not a separate mechanism. So call it after the elevation is set.

Locked means it cannot be SCOURED. Sediment may still settle onto it, which is what a concrete apron or a rock bar does in the field, and that new material is erodible again – it is above the base. If you want a cell that neither erodes nor accretes, that is not this.

Composes with set_erodible_base(): the base sets how DEEP erosion may go, the region sets WHERE it may happen, and setting one leaves the other in place. Where they disagree the stricter wins.

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

openmp - Python RK loop (use_c_rk_loop=False)
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_negative_volume_warning_fraction(fraction: float) → None[source]

Set the fraction of the total domain water volume that must be added by clamping negative-depth cells to zero depth in a single timestep before update_conserved_quantities() emits a “possible loss of conservation” warning.

Clamping a few cells by a near-zero depth is normal in wetting/drying and involves negligible volume, so the default (see anuga.config.negative_volume_warning_fraction) avoids warning on almost every step. Set to 0.0 to warn whenever any volume is added.

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.

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_sediment_friction(mode='constant', k_s=None, r_d=2.0, r_br=2.0, sigma_br=None, bed='sand', grain_size=None)[source]

Select the friction closure used by the sediment kernel (spec 3.3).

This affects only tau_b in the sediment source term; the hydrodynamic friction operator is untouched.

mode=’constant’ (default): f_c = g n^2 h^(-1/3) [T-6] with n from the domain’s own friction quantity.
mode=’larsen_lamb’: Manning-Strickler [T-13]-[T-15]: k_s = r_d r_br sigma_br and n = k_s^(1/6) / (8.1 sqrt(g)), uniform in space and time. Pass either k_s directly or sigma_br (one standard deviation of bedrock elevation). LL16 measured sigma_br ~ 5 m at Moses Coulee, giving k_s = 20 m and n = 0.065 – a site-measured quantity, not a universal default.
mode=’wilson’: f_c from bed type and relative submergence [T-8]-[T-10]. bed is ‘sand’ (uses D50), ‘gravel’ or ‘boulder’ (both use D84); pass the percentile as grain_size.

Notes

`larsen_lamb` and `wilson` are megaflood/planetary parameterisations, not general-purpose flood closures. W04 is a Mars outflow channel study under Martian gravity; LL16 is the Channeled Scablands. Neither is calibrated for ordinary river or urban flood modelling, which is ANUGA’s main use. For standard flood work keep constant with an n from conventional tables.

set_sediment_parameters(porosity=None, c_max=None, c_pack=None, bed_evolution=None, rho_w=None)[source]

Set the scalar sediment parameters, with validation.

Everything here is a physical property of the run, not a numerical knob. All are optional; only what you pass is changed.

Parameters:

porosity (float) – Bed porosity lambda in [G-4]. The sediment VOLUME leaving suspension is (1-lambda) dz, the rest being pore space filled from the water column. Default 0.30; LM15 use 0.28.
c_max (float) – [L-2], the ceiling on depth-averaged volumetric concentration. Default 0.30 (FG21); aS16 use 0.20.
c_pack (float) – [L-4], maximum packing bounding the NEAR-BED concentration c_b = d* c. Default 0.65, the same constant that bounds E* in [E-1]. Only bites when d* != 1.
bed_evolution (bool) – Spec 2.4’s coupling stage. True (default) evolves the bed via [G-4]/[G-5]; False is the FIXED BED of Phase 3, which is RDycore v1.0’s configuration and what the analytic constant-depth deposition solutions assume.
rho_w (float) – Water density, used to form the dimensional bed shear stress.

set_shear_closure(closure='quadratic_drag')[source]

Select how bed shear stress is obtained (spec 3.1 / 3.4).

‘quadratic_drag’ (default) – [T-1], tau_b = rho f_c |v|^2. Makes no equilibrium assumption.

‘depth_slope’ – [T-7], tau_b = rho g h S with S the bed slope magnitude, as aSM16 Eqs 6-7 and hence anugaSed. This is the steady uniform (normal) flow approximation: it assumes the energy slope equals the bed slope and the flow is locally in equilibrium.

‘energy_slope’ – [T-7e], the same tau_b = rho g h S with S the free-surface slope magnitude instead. Under the shallow-water assumption the free surface is the energy grade line, so this drops [T-7]’s equilibrium assumption and uses the slope actually driving the flow. Prefer it wherever the bed slope is not a good proxy for the energy slope: backwater, a pool-riffle sequence, a bed that is flat but drawing down, a dam break. It is also what the older Bed_shear_erosion_operator used (EN_slope), so it is the closure to pick when reproducing a model built on that operator.

Notes

Spec 3.4 recommends [T-1] and keeps [T-7] only for reproducing published anugaSed results, for three reasons: normal-flow equilibrium is exactly what fails in the dam-breach and outburst floods this work targets; S should be the energy slope, not the bed slope (substituting the energy slope into [T-7] recovers [T-1] identically); and the domain-global slope clamp anugaSed applies has no counterpart in their own manual.

This does not reproduce anugaSed exactly. Their code additionally divides the elevation gradient by a domain-mean cell size and applies S <- min(S, mean(S)/2). Neither is in aSM16; the first is dimensionally inconsistent and the second is the undocumented clamp of divergence D1a. [T-7] here follows the manual, not the code.

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_tracer(name, values)[source]

Set the concentration c of tracer name.

Seeds both c and the conserved m = h*c, which must agree: setting c alone leaves the conserved variable stale and the first substep overwrites c from it.

Because m = h*c, this reads the CURRENT depth: set stage and elevation before the tracer, not after.

Called before the run starts, this also rebases the conservation baseline (see check_tracer_conservation); called mid-run it does not.

set_tracer_boundary(name, tag, value)[source]

Concentration that tracer name brings in across boundary tag.

Only used where water FLOWS IN. The flux kernel picks the upwind value edge by edge from the sign of the water flux:

inflow (n.U into the domain) -> this boundary value outflow (n.U out of the domain) -> the interior edge value

so there is nothing to set for an outflow, and nothing that can be set: prescribing a concentration on an outflow would over-determine the advection, and the kernel would ignore it anyway. A boundary that only ever lets water out therefore needs no call at all, and one that alternates gets Dirichlet-on-inflow / transmissive-on-outflow automatically – the characteristic condition, without a switch.

Parameters:

name (str) – A registered tracer.
tag (str) – A boundary tag, as used by set_boundary. On a distributed sub-domain a tag this rank owns no part of is silently ignored, so the same call works on every rank; in serial an unknown tag is an error.
value (float, array-like, or callable) – A scalar applies to every edge carrying tag. An array must have one value per edge of that tag, ordered as domain.tag_boundary_cells[tag]. A callable is evaluated as value(t) at each timestep and must return a scalar or such an array, for a time-varying inflow.

Notes

UNSET BOUNDARIES BRING c = 0. The array is zero-filled at add_tracer, so with no call here an inflow carries clean water. That is a modelling assumption, not a neutral default: for salinity or suspended sediment it is usually wrong, and it is invisible in the output. It is preserved because changing it would silently alter existing results, but set it explicitly wherever water enters.

set_tracer_source(name, values)[source]

Prescribe an external source S_ms for tracer name, in m/s.

[G-3]’s optional external supply: hillslope yield, tributary load, rainfall washoff. Same units as E and D.

Applied after the [L-1]/[L-2] limiters. Those bound the bed exchange by what the bed and water column can supply; an external source is neither, so clipping it there would be wrong – and would also make a manufactured solution impossible to impose exactly.

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

Enable active-set stepping: skip cells that are provably unchangeable this step (dry, with an all-dry 2-ring neighbourhood).

Only the GPU/unified compute path honours the flag (rk2/DE1, ader2/DE_ader2 and euler/DE0 stepping; DE2/rk3 runs full steps with a notice). It is also ignored, with a notice, on domains with passive tracers or sediment classes, riverwalls, or sloped Manning friction: the active path uses the single-solve scatter flux kernel, which carries none of those. It also switches the flux kernel to single-solve scatter mode, whose results differ from the default cell-based kernel only at floating-point roundoff. Serial domains only for now – under MPI the flag is ignored with a warning.

On mostly-dry flood domains (dam break, levee breach, storm surge) this is the largest single speedup measured for the GPU path. Water added by rate operators (rainfall) is never lost – newly wetted cells activate on the next step – but widespread rain activates the whole mesh and the speedup degrades toward zero, hence the warning when both are in play.

set_use_optimise_dry_cells(flag: bool = True) → None[source]: Try to optimize calculations where region is dry

settling_velocity(diameter, rho_s=2650.0, rho_w=1000.0, nu=1e-06, C1=18.0, C2=0.4)[source]

Settling velocity v_s from Ferguson & Church (2004), spec [S-1].

v_s = R g d^2 / ( C1 nu + sqrt(0.75 C2 R g d^3) )

Smooth across the Stokes/turbulent transition and branch-free, which is why the spec prefers it over the Dietrich (1982) polynomial fit.

C1 = 18, C2 = 0.4 are the smooth-sphere constants; use 1.0, 1.1 for natural irregular grains. R = rho_s/rho_w - 1 is submerged specific gravity.

Verified against the spec: d = 4.5e-5 m quartz gives 1.75e-3 m/s, the value P13 report for Ferguson & Church.

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.