Source code¶
Here you can find the documentation with direct links to the actual AeoLiS code. You can click on the green [source] button next to the classes and modules below to access the specific source code. You can use ctr-f to look for a specific functionality or variable. It still may be a bit difficult to browse through, in addition you can find an overview of all module code here
Model classes¶
The AeoLiS model is based on two main model classes:
AeoLiS
and
AeoLiSRunner
. The former is the actual,
low-level, BMI-compatible class that implements the basic model
functions and numerical schemes. The latter is a convenience class
that implements a time loop, netCDF4 output, a progress indicator and
a callback function that allows the used to interact with the model
during runtime.
The additional WindGenerator
class to generate
realistic wind time series is available from the same module.
AeoLiS¶
- class model.AeoLiS(configfile: str)[source]¶
AeoLiS model class
AeoLiS is a process-based model for simulating supply-limited aeolian sediment transport. This model class is compatible with the Basic Model Interface (BMI) and provides basic model operations, like initialization, time stepping, finalization and data exchange. For higher level operations, like a progress indicator and netCDF4 output is refered to the AeoLiS model runner class, see
AeoLiSRunner
.Examples
>>> with AeoLiS(configfile='aeolis.txt') as model: >>> while model.get_current_time() <= model.get_end_time(): >>> model.update()
>>> model = AeoLiS(configfile='aeolis.txt') >>> model.initialize() >>> zb = model.get_var('zb') >>> model.set_var('zb', zb + 1) >>> for i in range(10): >>> model.update(60.) # step 60 seconds forward >>> model.finalize()
- __init__(configfile: str) None [source]¶
Initialize class
- Parameters:
configfile – Path to model configuration file. See
read_configfile()
.
- crank_nicolson() Any [source]¶
Convenience function for semi-implicit solver based on Crank-Nicolson scheme
See also
- static dimensions(var: str | None = None) Tuple | dict [source]¶
Static method that returns named dimensions of all spatial grids
- Parameters:
var (optional) – Name of spatial grid
- Returns:
Tuple with named dimensions of requested spatial grid or dictionary with all named dimensions of all spatial grids. Returns nothing if requested spatial grid is not defined.
- Return type:
tuple or dict
- euler_backward() Any [source]¶
Convenience function for implicit solver based on Euler backward scheme
See also
- euler_forward() Any [source]¶
Convenience function for explicit solver based on Euler forward scheme
See also
- get_var(var: str) ndarray | int | float | str | list [source]¶
Returns spatial grid or model configuration parameter
If the given variable name matches with a spatial grid, the spatial grid is returned. If not, the given variable name is matched with a model configuration parameter. If a match is found, the parameter value is returned. Otherwise, nothing is returned.
- Parameters:
var – Name of spatial grid or model configuration parameter
- Returns:
Spatial grid or model configuration parameter
- Return type:
np.ndarray or int, float, str or list
Examples
>>> # returns bathymetry grid ... model.get_var('zb')
>>> # returns simulation duration ... model.get_var('tstop')
See also
- get_var_name(i: int) str | int [source]¶
Returns name of spatial grid by index (in alphabetical order)
- Parameters:
i – Index of spatial grid
- Returns:
Name of spatial grid or -1 in case index exceeds the number of grids
- Return type:
str or -1
- get_var_rank(var: str) int [source]¶
Returns rank of spatial grid
- Parameters:
var – Name of spatial grid
- Return type:
Rank of spatial grid or -1 if not found
- get_var_shape(var: str) Tuple | str [source]¶
Returns shape of spatial grid
- Parameters:
var – Name of spatial grid
- Return type:
Dimensions of spatial grid or -1 if not found
- get_var_type(var: str) str | int [source]¶
Returns data type of variable in spatial grid
- Parameters:
var – Name of spatial grid
- Return type:
Variable type of spatial grid or -1 if not found
- initialize() None [source]¶
Read model configuration file, initialize parameters and spatial grids dictionary, and load bathymetry and bed composition.
- set_timestep(dt: float = -1.0) bool [source]¶
Determine optimal time step
If no time step is given the optimal time step is determined. For explicit numerical schemes the time step is based in the Courant-Frierichs-Lewy (CFL) condition. For implicit numerical schemes the time step specified in the model configuration file is used. Alternatively, a preferred time step is given that is maximized by the CFL condition in case of an explicit numerical scheme.
Returns True except when:
1. No time step could be determined, for example when there is no wind and the numerical scheme is explicit. In this case the time step is set arbitrarily to one second.
2. Or when the time step is smaller than -1. In this case the time is updated with the absolute value of the time step, but no model execution is performed. This funcionality can be used to skip fast-forward in time.
- Parameters:
df – Preferred time step
- Return type:
False if determination of time step was unsuccessful, True otherwise
- set_var(var: str, val: ndarray | int | float | str | list) None [source]¶
Sets spatial grid or model configuration parameter
If the given variable name matches with a spatial grid, the spatial grid is set. If not, the given variable name is matched with a model configuration parameter. If a match is found, the parameter value is set. Otherwise, nothing is set.
- Parameters:
var – Name of spatial grid or model configuration parameter
val – Spatial grid or model configuration parameter
Examples
>>> # set bathymetry grid ... model.set_var('zb', np.array([[0.,0., ... ,0.]]))
>>> # set simulation duration ... model.set_var('tstop', 3600.)
See also
- set_var_index(i: int, val: ndarray) None [source]¶
Set spatial grid by index (in alphabetical order)
- Parameters:
i – Index of spatial grid
val – Spatial grid
- set_var_slice() None [source]¶
Overwrite the values in variable name with data from var, in the range (start:start+count). Start, count can be integers for rank 1, and can be tuples of integers for higher ranks. For some implementations it can be equivalent and more efficient to do: get_var(name)[start[0]:start[0]+count[0], …, start[n]:start[n]+count[n]] = var
- solve(alpha: float = 0.5, beta: float = 1.0) dict [source]¶
Implements the explicit Euler forward, implicit Euler backward and semi-implicit Crank-Nicolson numerical schemes
Determines weights of sediment fractions, sediment pickup and instantaneous sediment concentration. Returns a partial spatial grid dictionary that can be used to update the global spatial grid dictionary.
- Parameters:
alpha – Implicitness coefficient (0.0 for Euler forward, 1.0 for Euler backward or 0.5 for Crank-Nicolson, default=0.5)
beta – Centralization coefficient (1.0 for upwind or 0.5 for centralized, default=1.0)
- Return type:
Partial spatial grid dictionary
Examples
>>> model.s.update(model.solve(alpha=1., beta=1.) # euler backward
>>> model.s.update(model.solve(alpha=.5, beta=1.) # crank-nicolson
- solve_pieter(alpha: float = 0.5, beta: float = 1.0) dict [source]¶
Implements the explicit Euler forward, implicit Euler backward and semi-implicit Crank-Nicolson numerical schemes
Determines weights of sediment fractions, sediment pickup and instantaneous sediment concentration. Returns a partial spatial grid dictionary that can be used to update the global spatial grid dictionary.
- Parameters:
alpha – Implicitness coefficient (0.0 for Euler forward, 1.0 for Euler backward or 0.5 for Crank-Nicolson, default=0.5)
beta (float, optional) – Centralization coefficient (1.0 for upwind or 0.5 for centralized, default=1.0)
- Return type:
Partial spatial grid dictionary
Examples
>>> model.s.update(model.solve(alpha=1., beta=1.) # euler backward
>>> model.s.update(model.solve(alpha=.5, beta=1.) # crank-nicolson
- update(dt: float = -1) None [source]¶
Time stepping function
Takes a single step in time. Interpolates wind and hydrodynamic time series to the current time, updates the soil moisture, mixes the bed due to wave action, computes wind velocity threshold and the equilibrium sediment transport concentration. Subsequently runs one of the available numerical schemes to compute the instantaneous sediment concentration and pickup for the next time step and updates the bed accordingly.
For explicit schemes the time step is maximized by the Courant-Friedrichs-Lewy (CFL) condition. See
set_timestep()
.- Parameters:
dt (optional) – Time step in seconds. The time step specified in the model configuration file is used in case dt is smaller than zero. For explicit numerical schemes the time step is maximized by the CFL confition.
AeoLiSRunner¶
- class model.AeoLiSRunner(configfile: str = 'aeolis.txt')[source]¶
AeoLiS model runner class
This runner class is a convenience class for the BMI-compatible AeoLiS model class (
AeoLiS()
). It implements a time loop, a progress indicator and netCDF4 output. It also provides the definition of a callback function that can be used to interact with the AeoLiS model during runtime.The command-line function
aeolis
is available that uses this class to start an AeoLiS model run.Examples
>>> # run with default settings ... AeoLiSRunner().run()
>>> AeoLiSRunner(configfile='aeolis.txt').run()
>>> model = AeoLiSRunner(configfile='aeolis.txt') >>> model.run(callback=lambda model: model.set_var('zb', zb))
>>> model.run(callback='bar.py:add_bar')
See also
console.aeolis
- __init__(configfile: str = 'aeolis.txt') None [source]¶
Initialize class
Reads model configuration file without parsing all referenced files for the progress indicator and netCDF output. If no configuration file is given, the default settings are used.
- Parameters:
configfile (str, optional) – Model configuration file. See
read_configfile()
.
- get_statistic(var: str, stat: str = 'avg') ndarray | None [source]¶
Return statistic of spatial grid
- Parameters:
var (str) – Name of spatial grid
stat (str) – Name of statistic (avg, sum, var, min or max)
- Return type:
Statistic of spatial grid
- get_var(var: str, clear: bool = True) ndarray | int | float | str | list [source]¶
Returns spatial grid, statistic or model configuration parameter
Overloads the
get_var()
function and extends it with the functionality to return statistics on spatial grids by adding a postfix to the variable name (e.g. Ct_avg). Supported statistics are avg, sum, var, min and max.- Parameters:
var (str) – Name of spatial grid or model configuration parameter. Spatial grid name can be extended with a postfix to request a statistic (_avg, _sum, _var, _min or _max).
clear (bool) – Clear output statistics afterwards.
- Returns:
Spatial grid, statistic or model configuration parameter
- Return type:
np.ndarray or int, float, str or list
Examples
>>> # returns average sediment concentration ... model.get_var('Ct_avg')
>>> # returns variance in wave height ... model.get_var('Hs_var')
See also
- initialize() None [source]¶
Initialize model
Overloads the
initialize()
function, but also initializes output statistics.
- load_hotstartfiles() None [source]¶
Load model state from hotstart files
Hotstart files are plain text representations of model state variables that can be used to hotstart the (partial) model state. Hotstart files should have the name of the model state variable it contains and have the extension .hotstart. Hotstart files differ from restart files in that restart files contain entire model states and are pickled Python objects.
See also
- load_restartfile(restartfile: str) bool [source]¶
Load model state from restart file
- Parameters:
restartfile (str) – Path to previously written restartfile.
- Return type:
True if model state from restartfile is loaded successfully
- output_clear() None [source]¶
Clears output statistics dictionary
Creates a matrix for minimum, maximum, variance and summed values for each output variable and sets the time step counter to zero.
- output_update() None [source]¶
Updates output statistics dictionary
Updates matrices with minimum, maximum, variance and summed values for each output variable with current spatial grid values and increases time step counter with one.
- output_write() None [source]¶
Appends output to netCDF4 output file
If the time since the last output is equal or larger than the set output interval, append current output to the netCDF4 output file. Computes the average and variance values based on available output statistics and clear output statistics dictionary.
- parse_callback(callback)[source]¶
Parses callback definition and returns function
The callback function can be specified in two formats:
As a native Python function
As a string refering to a Python script and function, separated by a colon (e.g.
example/callback.py:function
)
- Parameters:
callback (str or function) – Callback definition
- Returns:
Python callback function
- Return type:
function
- print_progress(fraction: float = 0.01, min_interval: float = 1.0, max_interval: float = 60.0) None [source]¶
Print progress to screen
- Parameters:
fraction (optional) – Fraction of simulation at which to print progress (default: 1%)
min_interval (optional) – Minimum time in seconds between subsequent progress prints (default: 1s)
max_interval (optional) – Maximum time in seconds between subsequent progress prints (default: 60s)
- run(callback=None, restartfile: str | None = None) None [source]¶
Start model time loop
Changes current working directory to the model directory, prints model configuration parameters and progress indicator to the screen, writes netCDF4 output and calls a callback function upon request.
- Parameters:
callback (str or function) – The callback function is called at the start of every single time step and takes the AeoLiS model object as input. The callback function can be used to interact with the model during simulation (e.g. update the bed with new measurements). See for syntax
parse_callback()
.restartfile (str) – Path to previously written restartfile. The model state is loaded from this file after initialization of the model.
See also
- update(dt: float = -1) None [source]¶
Time stepping function
Overloads the
update()
function, but also updates output statistics and clears output statistics upon request.- Parameters:
dt (float, optional) – Time step in seconds.
WindGenerator¶
- class model.WindGenerator(mean_speed: float = 9.0, max_speed: float = 30.0, dt: float = 60.0, n_states: int = 30, shape: float = 2.0, scale: float = 2.0)[source]¶
Wind velocity time series generator
Generates a random wind velocity time series with given mean and maximum wind speed, duration and time resolution. The wind velocity time series is generated using a Markov Chain Monte Carlo (MCMC) approach based on a Weibull distribution. The wind time series can be written to an AeoLiS-compatible wind input file assuming a constant wind direction of zero degrees.
The command-line function
aeolis-wind
is available that uses this class to generate AeoLiS wind input files.Examples
>>> wind = WindGenerator(mean_speed=10.).generate(duration=24*3600.) >>> wind.write_time_series('wind.txt') >>> wind.plot() >>> wind.hist()
See also
console.wind
- __getitem__(s: int) ndarray [source]¶
Retrieves item from list of wind speeds as an array
Parameters:¶
- s:
index in a list
- __init__(mean_speed: float = 9.0, max_speed: float = 30.0, dt: float = 60.0, n_states: int = 30, shape: float = 2.0, scale: float = 2.0) None [source]¶
- __weakref__¶
list of weak references to the object (if defined)
Physics modules¶
Bathymetry and bed composition¶
- bed.arrange_layers(m, dm, d, nl, ix_ero, ix_dep)[source]¶
Arranges mass redistrubution between layers. This function is called in the bed.update fucntion to speed up code using numba
- Parameters:
m (array) – mass in layers
dm (array) – total mass exchanged between layers derrived from pickup
d (array) – normalized mass in layers
nl (int) – number of layers
ix_dep (array) – cells for deposition
ix_ero (array) – cells for erosion
- Return type:
m
- bed.initialize(s, p)[source]¶
Initialize bathymetry and bed composition
Initialized bathymetry, computes cell sizes and orientation, bed layer thickness and bed composition.
- Parameters:
s (dict) – Spatial grids
p (dict) – Model configuration parameters
- Returns:
Spatial grids
- Return type:
dict
- bed.mixtoplayer(s, p)[source]¶
Mix grain size distribution in top layer of the bed.
Simulates mixing of the top layers of the bed by wave action. The wave action is represented by a local wave height maximized by a maximum wave hieght over depth ratio
gamma
. The mixing depth is a fraction of the local wave height indicated byfacDOD
. The mixing depth is used to compute the number of bed layers that should be included in the mixing. The grain size distribution in these layers is then replaced by the average grain size distribution over these layers.- Parameters:
s (dict) – Spatial grids
p (dict) – Model configuration parameters
- Returns:
Spatial grids
- Return type:
dict
- bed.prevent_negative_mass(m, dm, pickup)[source]¶
Handle situations in which negative mass may occur due to numerics
Negative mass may occur by moving sediment to lower layers down to accomodate deposition of sediments. In particular two cases are important:
A net deposition cell has some erosional fractions.
In this case the top layer mass is reduced according to the existing sediment distribution in the layer to accomodate deposition of fresh sediment. If the erosional fraction is subtracted afterwards, negative values may occur. Therefore the erosional fractions are subtracted from the top layer beforehand in this function. An equal mass of deposition fractions is added to the top layer in order to keep the total layer mass constant. Subsequently, the distribution of the sediment to be moved to lower layers is determined and the remaining deposits are accomodated.
Deposition is larger than the total mass in a layer.
In this case a non-uniform distribution in the bed may also lead to negative values as the abundant fractions are reduced disproportionally as sediment is moved to lower layers to accomodate the deposits. This function fills the top layers entirely with fresh deposits and moves the existing sediment down such that the remaining deposits have a total mass less than the total bed layer mass. Only the remaining deposits are fed to the routine that moves sediment through the layers.
- Parameters:
m (np.ndarray) – Sediment mass in bed (nx*ny, nl, nf)
dm (np.ndarray) – Total sediment mass exchanged between layers (nx*ny, nf)
pickup (np.ndarray) – Sediment pickup (nx*ny, nf)
- Returns:
np.ndarray – Sediment mass in bed (nx*ny, nl, nf)
np.ndarray – Total sediment mass exchanged between layers (nx*ny, nf)
np.ndarray – Sediment pickup (nx*ny, nf)
Note
The situations handled in this function can also be prevented by reducing the time step, increasing the layer mass or increasing the adaptation time scale.
- bed.update(s, p)[source]¶
Update bathymetry and bed composition
Update bed composition by moving sediment fractions between bed layers. The total mass in a single bed layer does not change as sediment removed from a layer is repleted with sediment from underlying layers. Similarly, excess sediment added in a layer is moved to underlying layers in order to keep the layer mass constant. The lowest bed layer exchanges sediment with an infinite sediment source that follows the original grain size distribution as defined in the model configuration file by
grain_size
andgrain_dist
. The bathymetry is updated following the cummulative erosion/deposition over the fractions ifbedupdate
isTrue
.- Parameters:
s (dict) – Spatial grids
p (dict) – Model configuration parameters
- Returns:
Spatial grids
- Return type:
dict
Wind velocity and direction¶
- wind.calculate_z0(p, s)[source]¶
Calculate z0 according to chosen roughness method
The z0 is required for the calculation of the shear velocity. Here, z0 is calculated based on a user-defined method. The constant method defines the value of z0 as equal to k (z0 = ks). This was implemented to ensure backward compatibility and does not follow the definition of Nikuradse (z0 = k / 30). For following the definition of Nikuradse use the method constant_nikuradse. The mean_grainsize_initial method uses the intial mean grain size ascribed to the bed (grain_dist and grain_size in the input file) to calculate the z0. The median_grainsize_adaptive bases the z0 on the median grain size (D50) in the surface layer in every time step. The resulting z0 is variable accross the domain (x,y). The strypsteen_vanrijn method is based on the roughness calculation in their paper.
- Parameters:
s (dict) – Spatial grids
p (dict) – Model configuration parameters
- Returns:
z0
- Return type:
array
- wind.compute_shear1d(s, p)[source]¶
Compute wind shear perturbation for given free-flow wind speed on computational grid. based on same implementation in Duna
- wind.interpolate(s, p, t)[source]¶
Interpolate wind velocity and direction to current time step
Interpolates the wind time series for velocity and direction to the current time step. The cosine and sine of the direction angle are interpolated separately to prevent zero-crossing errors. The wind velocity is decomposed in two grid components based on the orientation of each individual grid cell. In case of a one-dimensional model only a single positive component is used.
- Parameters:
s (dict) – Spatial grids
p (dict) – Model configuration parameters
t (float) – Current time
- Returns:
Spatial grids
- Return type:
dict
- class shear.WindShear(x, y, z, dx, dy, L, l, z0, buffer_width, buffer_relaxation=None)[source]¶
Class for computation of 2DH wind shear perturbations over a topography.
The class implements a 2D FFT solution to the wind shear perturbation on curvilinear grids. As the FFT solution is only defined on an equidistant rectilinear grid with circular boundary conditions that is aligned with the wind direction, a rotating computational grid is automatically defined for the computation. The computational grid is extended in all directions using a logistic sigmoid function as to ensure full coverage of the input grid for all wind directions, circular boundaries and preservation of the alongshore uniformity. An extra buffer distance can be used as to minimize the disturbence from the borders in the input grid. The results are interpolated back to the input grid when necessary.
Frequencies related to wave lengths smaller than a computational grid cell are filtered from the 2D spectrum of the topography using a logistic sigmoid tapering. The filtering aims to minimize the disturbance as a result of discontinuities in the topography that may physically exists, but cannot be solved for in the computational grid used.
Example
>>> w = WindShear(x, y, z) >>> w(u0=10., udir=30.).add_shear(taux, tauy)
Notes
To do:
- Actual resulting values are still to be compared with the results
from Kroy et al. (2002)
Grid interpolation can still be optimized
Separation bubble is still to be improved
- add_shear()[source]¶
Add wind shear perturbations to a given wind shear field
- Parameters:
taux (numpy.ndarray) – Wind shear in x-direction
tauy (numpy.ndarray) – Wind shear in y-direction
- Returns:
taux (numpy.ndarray) – Wind shear including perturbations in x-direction
tauy (numpy.ndarray) – Wind shear including perturbations in y-direction
- compute_shear(u0, nfilter=(1.0, 2.0))[source]¶
Compute wind shear perturbation for given free-flow wind speed on computational grid
- Parameters:
u0 (float) – Free-flow wind speed
nfilter (2-tuple) – Wavenumber range used for logistic sigmoid filter. See
filter_highfrequencies()
- filter_highfrequenies(kx, ky, hs, nfilter=(1, 2))[source]¶
Filter high frequencies from a 2D spectrum
A logistic sigmoid filter is used to taper higher frequencies from the 2D spectrum. The range over which the sigmoid runs from 0 to 1 with a precision
p
is given by the 2-tuplenfilter
. The range is defined as wavenumbers in terms of gridcells, i.e. a value 1 corresponds to a wave with lengthdx
.- Parameters:
kx (numpy.ndarray) – Wavenumbers in x-direction
ky (numpy.ndarray) – Wavenumbers in y-direction
hs (numpy.ndarray) – 2D spectrum
nfilter (2-tuple) – Wavenumber range used for logistic sigmoid filter
p (float) – Precision of sigmoid range definition
- Returns:
hs – Filtered 2D spectrum
- Return type:
numpy.ndarray
- static get_exact_grid(xmin, xmax, ymin, ymax, dx, dy)[source]¶
Returns a grid with given gridsizes approximately within given bounding box
- get_separation()[source]¶
Returns difference in height between z-coordinate of the separation polynomial and of the bed level
- Returns:
hsep – Height of seperation bubble
- Return type:
numpy.ndarray
- get_shear()[source]¶
Returns wind shear perturbation field
- Returns:
taux (numpy.ndarray) – Wind shear perturbation in x-direction
tauy (numpy.ndarray) – Wind shear perturbation in y-direction
- plot(ax=None, cmap='Reds', stride=10, computational_grid=False, **kwargs)[source]¶
Plot wind shear perturbation
- Parameters:
ax (matplotlib.pyplot.Axes, optional) – Axes to plot onto
cmap (matplotlib.cm.Colormap or string, optional) – Colormap for topography (default: Reds)
stride (int, optional) – Stride to apply to wind shear vectors (default: 10)
computational_grid (bool, optional) – Plot on computational grid rather than input grid (default: False)
kwargs (dict) – Additional arguments to
matplotlib.pyplot.quiver()
- Returns:
ax – Axes used for plotting
- Return type:
matplotlib.pyplot.Axes
- static rotate(x, y, alpha, origin=(0, 0))[source]¶
Rotate a matrix over given angle around given origin
Wind velocity threshold¶
- threshold.compute(s, p)[source]¶
Compute wind velocity threshold based on bed surface properties
Computes wind velocity threshold based on grain size fractions, bed slope, soil moisture content, air humidity, the presence of roughness elements and a non-erodible layer. All bed surface properties increase the current wind velocity threshold, except for the grain size fractions. Therefore, the computation is initialized by the grain size fractions and subsequently altered by the other bed surface properties.
- Parameters:
s (dict) – Spatial grids
p (dict) – Model configuration parameters
- Returns:
Spatial grids
- Return type:
dict
See also
compute_grainsize
,compute_bedslope
,compute_moisture
,compute_humidity
,compute_sheltering
,non_erodible
- threshold.compute_bedslope(s, p)[source]¶
Modify wind velocity threshold based on bed slopes following Dyer (1986)
- Parameters:
s (dict) – Spatial grids
p (dict) – Model configuration parameters
- Returns:
Spatial grids
- Return type:
dict
- threshold.compute_grainsize(s, p)[source]¶
Compute wind velocity threshold based on grain size fractions following Bagnold (1937)
- Parameters:
s (dict) – Spatial grids
p (dict) – Model configuration parameters
- Returns:
Spatial grids
- Return type:
dict
- threshold.compute_moisture(s, p)[source]¶
Modify wind velocity threshold based on soil moisture content following Belly (1964) or Hotta (1984)
- Parameters:
s (dict) – Spatial grids
p (dict) – Model configuration parameters
- Returns:
Spatial grids
- Return type:
dict
- threshold.compute_salt(s, p)[source]¶
Modify wind velocity threshold based on salt content following Nickling and Ecclestone (1981)
- Parameters:
s (dict) – Spatial grids
p (dict) – Model configuration parameters
- Returns:
Spatial grids
- Return type:
dict
- threshold.compute_sheltering(s, p)[source]¶
Modify wind velocity threshold based on the presence of roughness elements following Raupach (1993)
Raupach (1993) presents the following amplification factor for the shear velocity threshold due to the presence of roughness elements.
\[R_t = \frac{u_{*,th,s}}{u_{*,th,r}} = \sqrt{\frac{\tau_s''}{\tau}} = \frac{1}{\sqrt{\left( 1 - m \sigma \lambda \right) \left( 1 + m \beta \lambda \right)}}\]\(m\) is a constant smaller or equal to unity that accounts for the difference between the average stress on the bed surface \(\tau_s\) and the maximum stress on the bed surface \(\tau_s''\).
\(\beta\) is the stress partition coefficient defined as the ratio between the drag coefficient of the roughness element itself \(C_r\) and the drag coefficient of the bare surface without roughness elements \(C_s\).
\(\sigma\) is the shape coefficient defined as the basal area divided by the frontal area: \(\frac{A_b}{A_f}\). For hemispheres \(\sigma = 2\), for spheres \(\sigma = 1\).
\(\lambda\) is the roughness density defined as the number of elements per surface area \(\frac{n}{S}\) multiplied by the frontal area of a roughness element \(A_f\), also known as the frontal area index:
\[\lambda = \frac{n b h}{S} = \frac{n A_f}{S}\]If multiplied by \(\sigma\) the equation simplifies to the mass fraction of non-erodible elements:
\[\sigma \lambda = \frac{n A_b}{S} = \sum_{k=n_0}^{n_k} \hat{w}^{\mathrm{bed}}_k\]where \(k\) is the fraction index, \(n_0\) is the smallest non-erodible fraction, \(n_k\) is the largest non-erodible fraction and \(\hat{w}^{\mathrm{bed}}_k\) is the mass fraction of sediment fraction \(k\). It is assumed that the fractions are ordered by increasing size.
Substituting the derivation in the Raupach (1993) equation gives the formulation implemented in this function:
\[u_{*,th,r} = u_{*,th,s} * \sqrt{\left( 1 - m \sum_{k=n_0}^{n_k} \hat{w}^{\mathrm{bed}}_k \right) \left( 1 + m \frac{\beta}{\sigma} \sum_{k=n_0}^{n_k} \hat{w}^{\mathrm{bed}}_k \right)}\]- Parameters:
s (dict) – Spatial grids
p (dict) – Model configuration parameters
- Returns:
Spatial grids
- Return type:
dict
Tides, meteorology and soil moisture content¶
- hydro.Boussinesq(GW, DSWL, ds, GW_stat, K_gw, ne_gw, D_gw, shl_ix, bound, zb, process_seepage_face)[source]¶
Add description
- hydro.evaporation(s, p, met)[source]¶
Compute evaporation according to the Penman equation (Shuttleworth, 1993)
- Parameters:
s (dict) – Spatial grids
p (dict) – Model configuration parameters
met (dict) – meteorologial parameters T: Temperature, degrees Celsius Q : Global radiation, MJ/m2/d P : Atmospheric pressure, kPa U: Relative humidity, %
- Returns:
Evaporation (mm/day)
- Return type:
float
- hydro.interpolate(s, p, t)[source]¶
Interpolate hydrodynamic and meteorological conditions to current time step
Interpolates the hydrodynamic and meteorological time series to the current time step, if available. Meteorological parameters are stored as dictionary rather than a single value.
- Parameters:
s (dict) – Spatial grids
p (dict) – Model configuration parameters
t (float) – Current time
- Returns:
Spatial grids
- Return type:
dict
- hydro.saturation_pressure(T)[source]¶
Compute saturation pressure based on air temperature, Tetens equation
- Parameters:
T (float) – Air temperature in degrees Celcius
- Returns:
Saturation pressure
- Return type:
float
- hydro.update(s, p, dt, t)[source]¶
Update soil moisture content
Updates soil moisture content in all cells. The soil moisure content is computed either with the infiltration-method or surface_moist method. The infiltration method accounts for surface moisture as a function of runup and the subsequent infiltration and evaporation. The surface_moist method takes into account the effect of wave runup, precipitation, evaporation, infiltration, and capillary rise from the groundwater table.
- Parameters:
s (dict) – Spatial grids
p (dict) – Model configuration parameters
dt (float) – Current time step
- Returns:
Spatial grids
- Return type:
dict
Sediment transport¶
- transport.compute_weights(s, p)[source]¶
Compute weights for sediment fractions
Multi-fraction sediment transport needs to weigh the transport of each sediment fraction to prevent the sediment transport to increase with an increasing number of sediment fractions. The weighing is not uniform over all sediment fractions, but depends on the sediment availibility in the air and the bed and the bed interaction parameter
bi
.- Parameters:
s (dict) – Spatial grids
p (dict) – Model configuration parameters
- Returns:
Array with weights for each sediment fraction
- Return type:
numpy.ndarray
- transport.constant_grainspeed(s, p)[source]¶
Define saltation velocity u [m/s]
- Parameters:
s (dict) – Spatial grids
p (dict) – Model configuration parameters
- Returns:
Spatial grids
- Return type:
dict
- transport.duran_grainspeed(s, p)[source]¶
Compute grain speed according to Duran 2007 (p. 42)
- Parameters:
s (dict) – Spatial grids
p (dict) – Model configuration parameters
- Returns:
Spatial grids
- Return type:
dict
- transport.equilibrium(s, p)[source]¶
Compute equilibrium sediment concentration following Bagnold (1937)
- Parameters:
s (dict) – Spatial grids
p (dict) – Model configuration parameters
- Returns:
Spatial grids
- Return type:
dict
- transport.renormalize_weights(w, ix)[source]¶
Renormalizes weights for sediment fractions
Renormalizes weights for sediment fractions such that the sum of all weights is unity. To ensure that the erosion of specific fractions does not exceed the sediment availibility in the bed, the normalization only modifies the weights with index equal or larger than
ix
.- Parameters:
w (numpy.ndarray) – Array with weights for each sediment fraction
ix (int) – Minimum index to be modified
- Returns:
Array with weights for each sediment fraction
- Return type:
numpy.ndarray
Avalanching¶
- avalanching.angele_of_repose(s, p)[source]¶
Determine the dynamic and static angle of repose.
Both the critical dynamic and static angle of repose are spatial varying and depend on surface moisture content and roots of present vegetation and ….
- Parameters:
s (dict) – Spatial grids
p (dict) – Model configuration parameters
- Returns:
Spatial grids
- Return type:
dict
- avalanching.avalanche(s, p)[source]¶
Avalanching occurs if bed slopes exceed critical slopes.
Simulates the process of avalanching that is triggered by the exceedence of a critical static slope
theta_stat
by the bed slope. The iteration stops if the bed slope does not exceed the dynamic critical slopetheta_dyn
.- Parameters:
s (dict) – Spatial grids
p (dict) – Model configuration parameters
- Returns:
Spatial grids
- Return type:
dict
Vegetation¶
Fences¶
- fences.calc_fence_shear(igrid, cgrid, udir, x0, y0, p)[source]¶
Compute wind shear for given wind speed and direction
- Parameters:
u0 (float) – Free-flow wind speed
udir (float) – Wind direction in degrees
process_separattion –
- fences.compute_fenceshear(igrid, cgrid, udir, p)[source]¶
Compute wind shear perturbation for given free-flow wind speed on computational grid
- Parameters:
u0 (float) – Free-flow wind speed
nfilter (2-tuple) – Wavenumber range used for logistic sigmoid filter. See
filter_highfrequencies()
- fences.get_exact_grid(xmin, xmax, ymin, ymax, dx, dy)[source]¶
Returns a grid with given gridsizes approximately within given bounding box
- fences.populate_computational_grid(igrid, cgrid, alpha, x0, y0)[source]¶
Interpolate input topography to computational grid
Adds and fills buffer zone around the initial grid and rotates the computational grid to current wind direction. The computational grid is filled by interpolating the input topography and initial wind induced shear stresses to it.
- Parameters:
alpha (float) – Rotation angle in degrees
Marine Erosion¶
- erosion.run_ph12(s, p, t)[source]¶
Calculates bed level change due to dune erosion
Calculates bed level change due to dune erosion accoording to Palmsten and Holman (2012).
- Parameters:
s (dict) – Spatial grids
p (dict) – Model configuration parameters
t (float) – Model time
- Returns:
Spatial grids
- Return type:
dict
Helper modules¶
Input/Output¶
- inout.check_configuration(p)[source]¶
Check model configuration validity
Checks if required parameters are set and if the references files for bathymetry, wind, tide and meteorological input are valid. Throws an error if one or more requirements are not met.
- Parameters:
p (dict) – Model configuration dictionary with parsed files
See also
- inout.parse_value(val, parse_files=True, force_list=False)[source]¶
Casts a string to the most appropriate variable type
- Parameters:
val (str) – String representation of value
parse_files (bool) – If True, files referred to by string parameters are parsed by
numpy.loadtxt
force_list – If True, interpret the value as a list, even if it consists of a single value
- Returns:
Casted value
- Return type:
misc
Examples
>>> type(parse_value('T')) bool >>> type(parse_value('F')) bool >>> type(parse_value('123')) int >>> type(parse_value('123.2')) float >>> type(parse_value('euler_forward')) str >>> type(parse_value('')) NoneType >>> type(parse_value('zb zs Ct')) numpy.ndarray >>> type(parse_value('zb', force_list=True)) numpy.ndarray >>> type(parse_value('0.1 0.2 0.3')[0]) float >>> type(parse_value('wind.txt'), parse_files=True) numpy.ndarray >>> type(parse_value('wind.txt'), parse_files=False) str
- inout.read_configfile(configfile, parse_files=True, load_defaults=True)[source]¶
Read model configuration file
Updates default model configuration based on a model configuration file. The model configuration file should be a text file with one parameter on each line. The parameter name and value are seperated by an equal sign (=). Any lines that start with a percent sign (%) or do not contain an equal sign are omitted.
Parameters are casted into the best matching variable type. If the variable type is
str
it is optionally interpreted as a filename. If the corresponding file is found it is parsed using thenumpy.loadtxt
function.- Parameters:
configfile (str) – Model configuration file
parse_files (bool) – If True, files referred to by string parameters are parsed
load_defaults (bool) – If True, default settings are loaded and overwritten by the settings from the configuration file
- Returns:
Dictionary with casted and optionally parsed model configuration parameters
- Return type:
dict
See also
- inout.visualize_grid(s, p)[source]¶
Create figures and tables for the user to check whether the grid-input is correctly interpreted
- inout.visualize_spatial(s, p)[source]¶
Create figures and tables for the user to check whether the input is correctly interpreted
- inout.visualize_timeseries(p, t)[source]¶
Create figures and tables for the user to check whether the timeseries-input is correctly interpreted
- inout.write_configfile(configfile, p=None)[source]¶
Write model configuration file
Writes model configuration to file. If no model configuration is given, the default configuration is written to file. Any parameters with a name ending with _file and holding a matrix are treated as separate files. The matrix is then written to an ASCII file using the
numpy.savetxt
function and the parameter value is replaced by the name of the ASCII file.- Parameters:
configfile (str) – Model configuration file
p (dict, optional) – Dictionary with model configuration parameters
- Returns:
Dictionary with casted and optionally parsed model configuration parameters
- Return type:
dict
See also
netCDF4 output¶
- netcdf.append(outputfile, variables)[source]¶
Append variables to existing netCDF4 output file
Increments the time axis length with one and appends the provided spatial grids along the time axis. The
variables
dictionary should at least have thetime
field indicating the current simulation time. The CF time bounds are updated accordingly.- Parameters:
outputfile (str) – Name of netCDF4 output file
variables (dict) – Dictionary with spatial grids and time
Examples
>>> netcdf.append('aeolis.nc', {'time', 3600., ... 'Ct', np.array([[0.,0., ... ,0.]]), ... 'Cu', np.array([[1.,1., ... ,1.]]))
See also
- netcdf.dump(outputfile, dumpfile, var='mass', ix=-1)[source]¶
Dumps time slice from netCDF4 output file to ASCII file
This function can be used to use a specific time slice from a netCDF4 output file as input file for another AeoLiS model run. For example, the bed composition from a spinup run can be used as initial composition for other runs reducing the spinup time.
- Parameters:
outputfile (str) – Name of netCDF4 output file
dumpfile (str) – Name of ASCII dump file
var (str, optional) – Name of spatial grid to be dumped (default: mass)
ix (int) – Time slice index to be dumped (default: -1)
Examples
>>> # use bedcomp_file = bedcomp.txt in model configuration file ... netcdf.dump('aeolis.nc', 'bedcomp.txt', var='mass')
- netcdf.initialize(outputfile, outputvars, s, p, dimensions)[source]¶
Create empty CF-compatible netCDF4 output file
- Parameters:
outputfile (str) – Name of netCDF4 output file
outputvars (dictionary) – Spatial grids to be written to netCDF4 output file
s (dict) – Spatial grids
p (dict) – Model configuration parameters
dimensions (dict) – Dictionary that specifies a tuple with the named dimensions for each spatial grid (e.g. (‘ny’, ‘nx’, ‘nfractions’))
Examples
>>> netcdf.initialize('aeolis.nc', ... ['Ct', 'Cu', 'zb'], ... ['avg', 'max'], ... s, p, {'Ct':('ny','nx','nfractions'), ... 'Cu':('ny','nx','nfractions'), ... 'zb':('ny','nx')})
- netcdf.parse_metadata(outputvars)[source]¶
Parse metadata from constants.py
Parses the Python comments in constants.py to extract meta data, like units, for the model state variables that can be used as netCDF4 meta data.
- Parameters:
outputvars (dictionary) – Spatial grids to be written to netCDF4 output file
- Returns:
meta – Dictionary with meta data for the output variables
- Return type:
dict
Plotting¶
Command-line tools¶
Debugging¶
Miscellaneous¶
- utils.apply_mask(arr, mask)[source]¶
Apply complex mask
The real part of the complex mask is multiplied with the input array. Subsequently the imaginary part is added and the result returned.
The shape of the mask is assumed to match the first few dimensions of the input array. If the input array is larger than the mask, the mask is repeated for any additional dimensions.
- Parameters:
arr (numpy.ndarray) – Array or matrix to which the mask needs to be applied
mask (numpy.ndarray) – Array or matrix with complex mask values
- Returns:
arr – Array or matrix to which the mask is applied
- Return type:
numpy.ndarray
- utils.calc_grain_size(p, s, percent)[source]¶
Calculate grain size characteristics based on mass in each fraction
Calculate grain size distribution for each cell based on weight distribution over the fractions. Interpolates to the requested percentage in the grain size distribution. For example, percent=50 will result in calculation of the D50. Calculation is only executed for the top layer
- Parameters:
s (dict) – Spatial grids
p (dict) – Model configuration parameters
percent (float) – Requested percentage in grain size dsitribution
- Returns:
grain size per grid cell
- Return type:
array
- utils.calc_mean_grain_size(p, s)[source]¶
Calculate mean grain size based on mass in each fraction
Calculate mean grain size for each cell based on weight distribution over the fractions. Calculation is only executed for the top layer.
- Parameters:
s (dict) – Spatial grids
p (dict) – Model configuration parameters
percent (float) – Requested percentage in grain size dsitribution
- Returns:
mean grain size per grid cell
- Return type:
array
- utils.format_log(msg, ncolumns=2, **props)[source]¶
Format log message into columns
Prints log message and additional data into a column format that fits into a 70 character terminal.
- Parameters:
msg (str) – Main log message
ncolumns (int) – Number of columns
props (key/value pairs) – Properties to print in column format
- Returns:
Formatted log message
- Return type:
str
Note
Properties names starting with
min
,max
ornr
are respectively replaced bymin.
,max.
or#
.
- utils.interp_array(x: ndarray, xp: ndarray, fp: ndarray, circular: bool = False, **kwargs: dict) ndarray [source]¶
Interpolate multiple time series at once
- Parameters:
x (array_like) – The x-coordinates of the interpolated values.
xp (1-D sequence of floats) – The x-coordinates of the data points, must be increasing.
fp (2-D sequence of floats) – The y-coordinates of the data points, same length as
xp
.circular (bool) – Use the
interp_circular()
function rather than thenumpy.interp()
function.kwargs (dict) – Keyword options to the
numpy.interp()
function
- Returns:
The interpolated values, same length as second dimension of
fp
.- Return type:
ndarray
- utils.interp_circular(x: ndarray, xp: ndarray, fp: ndarray, **kwargs) ndarray [source]¶
One-dimensional linear interpolation.
Returns the one-dimensional piecewise linear interpolant to a function with given values at discrete data-points. Values beyond the limits of
x
are interpolated in circular manner. For example, a value ofx > x.max()
evaluates asf(x-x.max())
assuming thatx.max() - x < x.max()
.- Parameters:
x (array_like) – The x-coordinates of the interpolated values.
xp (1-D sequence of floats) – The x-coordinates of the data points, must be increasing.
fp (1-D sequence of floats) – The y-coordinates of the data points, same length as
xp
.kwargs (dict) – Keyword options to the
numpy.interp()
function
- Returns:
y – The interpolated values, same shape as
x
.- Return type:
{float, ndarray}
- Raises:
ValueError – If
xp
andfp
have different length
- utils.normalize(x: ndarray, ref: ndarray | None = None, axis: int = 0, fill: float = 0.0)[source]¶
Normalize array
Normalizes an array to make it sum to unity over a specific axis. The procedure is safe for dimensions that sum to zero. These dimensions return the
fill
value instead.- Parameters:
x (array_like) – The array to be normalized
ref (array_like, optional) – Alternative normalization reference, if not specified, the sum of x is used
axis (optional) – The normalization axis (default: 0)
fill (optional) – The return value for all-zero dimensions (default: 0.)
- utils.prevent_tiny_negatives(x: ndarray, max_error: float = 1e-10, replacement: float = 0.0) ndarray [source]¶
Replace tiny negative values in array
- Parameters:
x (np.ndarray) – Array with potential tiny negative values
max_error (float) – Maximum absolute value to be replaced
replacement (float) – Replacement value
- Returns:
Array with tiny negative values removed
- Return type:
np.ndarray
- utils.print_value(val, fill='<novalue>')[source]¶
Construct a string representation from an arbitrary value
- Parameters:
val (misc) – Value to be represented as string
fill (str, optional) – String representation used in case no value is given
- Returns:
String representation of value
- Return type:
str
- utils.rotate(x, y, alpha, origin=(0, 0))[source]¶
Rotate a matrix over given angle around given origin
- gridparams.initialize(s: Dict, p: Dict) Tuple[Dict, Dict] [source]¶
Initializes a grid and a checks if the grid is properly defined based on the x and y input files. Grid distance and grid surface are also initialized, and then broadcasted to the s-dictionary for use by other modules.
- Parameters:
s (dict) – Spatial grids
p (dict) – Model configuration parameters
- Returns:
Initialized spatial grids as a tuple containing the s-dictionary and the p-dictionary.
- Return type:
tuple [dict, dict]
Note
Sierd’s favorite function is:
aeolis.bed.prevent_tiny_negatives