import math
from collections.abc import Callable
from collections.abc import Iterator
from numbers import Integral
import numpy as np
from numpy.typing import ArrayLike
from numpy.typing import NDArray
from requireit import require_between
from requireit import require_instance
from requireit import require_positive
from landlab.field.errors import FieldError
from landlab.field.graph_field import GraphFields
[docs]
def iter_time_steps(duration: float, *, dt: float | None = None) -> Iterator[float]:
"""Yield fixed-size time steps that evenly span a requested duration.
Split *duration* into equally-sized substeps, so that no substeps
are longer that *dt*.
Parameters
----------
duration : float
Total amount of time to advance.
dt : float, optional
Maximum time-step size. If not given, use *duration* as a
single time step.
Yields
------
float
The next time-step size.
Raises
------
ValueError
If *duration* is negative, or if *dt* is not finite and
positive.
Examples
--------
>>> from landlab.core.component_utils import iter_time_steps
>>> list(iter_time_steps(10.0, dt=2.5))
[2.5, 2.5, 2.5, 2.5]
A duration that doesn't divide evenly is split into equal substeps,
each no longer than *dt*, rather than leaving a short final step.
>>> list(iter_time_steps(10.0, dt=3.0))
[2.5, 2.5, 2.5, 2.5]
If *dt* isn't given, *duration* is used as a single time step.
>>> list(iter_time_steps(5.0))
[5.0]
"""
duration = require_between(
duration, 0.0, math.inf, inclusive_max=False, name="duration"
)
if duration == 0.0:
return
dt = duration if dt is None else dt
dt = require_between(
dt, 0.0, math.inf, inclusive_min=False, inclusive_max=False, name="dt"
)
n_steps = math.ceil(duration / dt)
step = duration / n_steps
for _ in range(n_steps):
yield step
[docs]
def iter_adaptive_time_steps(
duration: float,
*,
calc_dt: Callable[[], float | None],
max_steps: int | None = None,
rtol: float = 1e-12,
) -> Iterator[float]:
"""Yield adaptive time steps that advance up to a requested duration.
Repeatedly call *calc_dt* to obtain the next stable time-step size,
capping each step so that the total does not exceed *duration*.
Note that iteration may stop within the tolerance specified by *rtol*.
Parameters
----------
duration : float
Total amount of time to advance.
calc_dt : callable
Called with no arguments before each substep to obtain the
current stable time-step size. A return value of ``None`` signals
that iteration should be stopped before *duration* is reached.
A return value of `inf` advances to *duration*.
max_steps : int, optional
Maximum number of substeps to yield before raising a
``RuntimeError``.
rtol : float, optional
Stop once the remaining time is no greater than
``rtol * duration``. Consequently, the yielded time steps may
sum to slightly less than *duration*.
Yields
------
float
The next time-step size.
Raises
------
ValueError
If *duration* or *rtol* are out of range.
RuntimeError
If *max_steps* is exceeded, or if a returned step is either
invalid or is too small, relative to the elapsed time, to make
further progress.
Examples
--------
>>> from landlab.core.component_utils import iter_adaptive_time_steps
>>> steps = iter([2.0, 2.0, 2.0, 1.0])
>>> list(iter_adaptive_time_steps(7.0, calc_dt=lambda: next(steps)))
[2.0, 2.0, 2.0, 1.0]
Return ``None`` from *calc_dt* to stop before *duration* is reached.
>>> steps = iter([2.0, 2.0, None])
>>> list(iter_adaptive_time_steps(10.0, calc_dt=lambda: next(steps)))
[2.0, 2.0]
>>> list(iter_adaptive_time_steps(10.0, calc_dt=lambda: 3.0))
[3.0, 3.0, 3.0, 1.0]
"""
duration = require_between(
duration, 0.0, math.inf, inclusive_max=False, name="duration"
)
rtol = require_between(rtol, 0.0, 1.0, inclusive_max=False, name="rtol")
if max_steps is not None:
require_instance(max_steps, Integral, name="max_steps")
require_positive(max_steps, name="max_steps")
if duration == 0.0:
return
tol = rtol * duration
elapsed = 0.0
n_steps = 0
while duration - elapsed > tol:
if max_steps is not None and n_steps >= max_steps:
raise RuntimeError(f"unable to reach {duration!r} in {max_steps} substeps")
step = calc_dt()
if step is None:
break
if step <= 0.0 or math.isnan(step):
raise RuntimeError("step must be positive or None")
this_dt = min(step, duration - elapsed)
new_elapsed = elapsed + this_dt
if new_elapsed == elapsed:
raise RuntimeError(
"time step is too small relative to the elapsed time"
" to make further progress due to floating-point precision"
)
yield this_dt
elapsed = new_elapsed
n_steps += 1
type FieldLike = str | ArrayLike
[docs]
def validate_field(
value: FieldLike,
*,
grid: GraphFields | None = None,
at: str | None = None,
) -> str | NDArray:
"""Validate a field name or array-like value.
Check that *value* is either the name of an existing field, or an
array-like object that could represent values located on a grid.
Unlike :func:`resolve_field`, a field name is *not* resolved to its
values, so the returned value can be stored and later passed to
:func:`resolve_field` to pick up the field's current
values at the time it's needed. *value* is converted to a *numpy*
array, but, unlike :func:`resolve_field`, it is not broadcast.
Parameters
----------
value : str or array_like
A field name, or an array of values.
grid : GraphFields, optional
Grid used to validate *value* against. If not given, *value* is
not checked against a grid (and, if a field name, is returned
unchecked).
at : str, optional
Name of the group (e.g. "node", "link") that *value* is defined
on. Must be given if, and only if, *grid* is given.
Returns
-------
str or ndarray
*value* unchanged, if it is a field name, otherwise *value* as an
array (not copied, if *value* is already an ndarray).
Raises
------
ValueError
If only one of *grid* and *at* is given, or if *value* is
array-like (and not a scalar) but its first dimension does not
match the number of *at* elements of *grid*.
landlab.field.errors.FieldError
If *value* is a field name that does not exist in *grid* at *at*.
Examples
--------
>>> from landlab.core.component_utils import validate_field
>>> from landlab import RasterModelGrid
>>> grid = RasterModelGrid((3, 4))
>>> _ = grid.add_ones("foo", at="node")
>>> validate_field("foo", grid=grid, at="node")
'foo'
>>> validate_field(range(12), grid=grid, at="node")
array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11])
>>> validate_field(42, grid=grid, at="node")
array(42)
"""
if (grid is None) != (at is None):
raise ValueError("grid and at must both be given, or neither")
if isinstance(value, str):
if grid is not None and value not in grid[at]:
raise FieldError(f"{value!r} is not a field at {at!r}")
return value
value = np.asarray(value)
if at is not None:
size = grid.size(at)
if value.ndim > 0 and value.shape[0] != size:
raise ValueError(
f"first dimension of value, {value.shape}, does not match the"
f" number of {at!r} elements of grid, {size}"
)
return value
[docs]
def resolve_field(
field: FieldLike,
*,
grid: GraphFields,
at: str,
) -> NDArray:
"""Resolve a field name or array to an array of values.
If *field* is a field name, look up and return its current values on
*grid*. If *field* is 0-dimensional, broadcast it to the number of
*at* elements of *grid*. Otherwise, return *field* unchanged. Call
this each time you need the values of a *field* previously validated with
:func:`validate_field`, so that a field name always resolves to that
field's current values, even if they have changed since *field* was
validated.
Parameters
----------
field : str or array_like
A field name, or an array of values, as returned by
:func:`validate_field`.
grid : GraphFields
Grid to look up *field* on, or broadcast it against, as needed.
at : str
Name of the group (e.g. "node", "link") that *field* is defined
on.
Returns
-------
ndarray
The current values of *field*. If *field* is 0-dimensional, this
is a read-only, broadcast view rather than a newly allocated array.
If *field* is already an *numpy* array, it is returned unchanged
(not copied). Otherwise, a new array is created.
Raises
------
landlab.field.errors.GroupError
If *at* is not a group on *grid*.
landlab.field.errors.FieldError
If *field* is a field name that does not exist in *grid* at *at*.
Examples
--------
>>> from landlab.core.component_utils import resolve_field
>>> from landlab import RasterModelGrid
>>> grid = RasterModelGrid((3, 4))
>>> foo = grid.add_ones("foo", at="node")
>>> resolve_field("foo", grid=grid, at="node") is foo
True
>>> resolve_field(range(12), grid=grid, at="node")
array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11])
>>> resolve_field(42, grid=grid, at="node")
array([42, 42, 42, 42, 42, 42, 42, 42, 42, 42, 42, 42])
"""
if isinstance(field, str):
return grid.field_values(field, at=at)
field = np.asarray(field)
return np.broadcast_to(field, (grid.size(at),)) if field.ndim == 0 else field