diff --git a/landlab/components/__init__.py b/landlab/components/__init__.py
index e9c3396245..7a607f4edd 100644
--- a/landlab/components/__init__.py
+++ b/landlab/components/__init__.py
@@ -60,6 +60,7 @@
from .priority_flood_flow_router import PriorityFloodFlowRouter
from .profiler import ChannelProfiler, Profiler, TrickleDownProfiler
from .radiation import Radiation
+from .river_flow_dynamics import river_flow_dynamics
from .sink_fill import SinkFiller, SinkFillerBarnes
from .soil_moisture import SoilInfiltrationGreenAmpt, SoilMoisture
from .space import Space, SpaceLargeScaleEroder
@@ -138,6 +139,7 @@
PrecipitationDistribution,
Profiler,
Radiation,
+ river_flow_dynamics,
SedDepEroder,
SedimentPulserAtLinks,
SedimentPulserEachParcel,
diff --git a/landlab/components/river_flow_dynamics/__init__.py b/landlab/components/river_flow_dynamics/__init__.py
new file mode 100644
index 0000000000..ba4d1f1711
--- /dev/null
+++ b/landlab/components/river_flow_dynamics/__init__.py
@@ -0,0 +1,3 @@
+from .river_flow_dynamics import river_flow_dynamics
+
+__all__ = ["river_flow_dynamics"]
diff --git a/landlab/components/river_flow_dynamics/river_flow_dynamics.py b/landlab/components/river_flow_dynamics/river_flow_dynamics.py
new file mode 100644
index 0000000000..3d83bc3c63
--- /dev/null
+++ b/landlab/components/river_flow_dynamics/river_flow_dynamics.py
@@ -0,0 +1,1906 @@
+"""river_flow_dynamics.py
+
+This component runs a semi-implicit, semi-Lagrangian
+finite-volume approximation to the depth-averaged shallow
+water equations of Casulli and Cheng (1992) and related work.
+
+Written by Sebastian Bernal and Angel Monsalve.
+
+Last updated: August 11, 2023
+
+Examples
+--------
+
+Simulation of River Flow Dynamics
+
+This example demonstrates how to simulate river flow dynamics using the Landlab library.
+
+First, let's import necessary libraries such as NumPy, Matplotlib, and Landlab components.
+
+>>> import numpy as np
+>>> import matplotlib.pyplot as plt
+>>> from landlab import RasterModelGrid
+>>> from landlab.components import river_flow_dynamics
+>>> from landlab.io import read_esri_ascii
+>>> from landlab.plot.imshow import imshow_grid
+
+Now, create a rectangular grid for flow dynamics calculations with specified dimensions and cell size.
+
+>>> nRows = 20
+>>> nCols = 60
+>>> cellSize = 0.1
+>>> grid = RasterModelGrid((nRows,nCols), xy_spacing=(cellSize,cellSize))
+
+Defining the Topography: Set up the initial topographic elevation for the grid, creating a basic rectangular channel with a slope of 0.01.
+
+>>> te = grid.add_zeros("topographic__elevation", at="node")
+>>> te += 0.059 - 0.01*grid.x_of_node
+>>> te[grid.y_of_node > 1.5] = 1.0
+>>> te[grid.y_of_node < 0.5] = 1.0
+
+We can use `imshow_grid` to visualize the elevation profile.
+Let's take a look at the grid's middle longitudinal section.
+
+>>> middleBedProfile = np.reshape(te,(nRows,nCols))[10,:]
+>>> np.round(middleBedProfile,3)
+array([ 0.059, 0.058, 0.057, 0.056, 0.055, 0.054, 0.053, 0.052,
+ 0.051, 0.05 , 0.049, 0.048, 0.047, 0.046, 0.045, 0.044,
+ 0.043, 0.042, 0.041, 0.04 , 0.039, 0.038, 0.037, 0.036,
+ 0.035, 0.034, 0.033, 0.032, 0.031, 0.03 , 0.029, 0.028,
+ 0.027, 0.026, 0.025, 0.024, 0.023, 0.022, 0.021, 0.02 ,
+ 0.019, 0.018, 0.017, 0.016, 0.015, 0.014, 0.013, 0.012,
+ 0.011, 0.01 , 0.009, 0.008, 0.007, 0.006, 0.005, 0.004,
+ 0.003, 0.002, 0.001, -0. ])
+
+Instantiating the component: Initialize the river_flow_dynamics component with specified parameters, including time step and Manning's roughness coefficient.
+
+The grid will need some data to run the river_flow_dynamics component.
+To check the names of the required inputs use the *input_var_names*
+class property.
+
+>>> river_flow_dynamics.input_var_names
+('surface_water__depth',
+ 'surface_water__elevation',
+ 'surface_water__velocity',
+ 'topographic__elevation')
+
+Also, if we want to know where these fields are mapped use:
+
+>>> river_flow_dynamics.var_mapping
+(('surface_water__depth', 'node'),
+ ('surface_water__elevation', 'node'),
+ ('surface_water__velocity', 'link'),
+ ('topographic__elevation', 'node'))
+
+Create fields of data for each of these input variables. The channel is
+empty at the beginning of the simulation, so we create the water depth field.
+
+>>> h = grid.add_zeros("surface_water__depth", at="node")
+
+Water velocity is zero everywhere since there is no water yet.
+
+>>> vel = grid.add_zeros("surface_water__velocity", at="link")
+
+Calculating the initial water surface elevation from water depth and topographic elevation.
+
+>>> wse = grid.add_zeros("surface_water__elevation", at="node")
+>>> wse += h + te
+
+Then, we specify the nodes at which water is entering into the domain, and also the associated links.
+These will be the inlet boundary conditions for water depth and velocity.
+In this case, water flows from left to right at 0.5 m depth, with a velocity of 0.45 m/s.
+
+>>> fixed_entry_nodes = np.arange(300,910,60)
+>>> fixed_entry_links = grid.links_at_node[fixed_entry_nodes][:,0]
+
+And we set the fixed values in these entry nodes/links.
+
+>>> entry_nodes_h_values = np.full(11, 0.5)
+>>> entry_links_vel_values = np.full(11, 0.45)
+
+We instantiate river_flow_dynamics with the arguments we defined previously:
+
+>>> rfd = river_flow_dynamics(grid, dt=0.1, mannings_n=0.012, fixed_entry_nodes=fixed_entry_nodes,
+... fixed_entry_links=fixed_entry_links,entry_nodes_h_values=entry_nodes_h_values,
+... entry_links_vel_values=entry_links_vel_values)
+
+And we run the simulation for 100 timesteps (10 seconds):
+
+>>> n_timesteps = 100
+>>> for timestep in range(n_timesteps):
+... rfd.run_one_step()
+
+Let's see what the flow depth at the center of the channel is after 10 seconds.
+
+>>> np.reshape(grid['node']["surface_water__depth"],(nRows,nCols))[10,:]
+array([ 0.5 , 0.49119496, 0.48070117, 0.47309464, 0.46764464,
+ 0.46370528, 0.46076177, 0.45840315, 0.45619896, 0.45419369,
+ 0.45250159, 0.45106348, 0.44972981, 0.4484113 , 0.44702824,
+ 0.445566 , 0.44414898, 0.44279566, 0.44144646, 0.44007219,
+ 0.43868299, 0.43726816, 0.43582451, 0.4344022 , 0.43300387,
+ 0.43159008, 0.43012984, 0.42867334, 0.42723278, 0.42577057,
+ 0.42426757, 0.42277662, 0.42127323, 0.41977445, 0.41825745,
+ 0.41674194, 0.41522539, 0.41368181, 0.41214349, 0.41060916,
+ 0.40905859, 0.40748906, 0.40591859, 0.40435914, 0.40275767,
+ 0.40113684, 0.39954249, 0.39795401, 0.39635525, 0.39474439,
+ 0.39313631, 0.39155279, 0.38993334, 0.38828883, 0.38664499,
+ 0.3849828 , 0.38333126, 0.381655 , 0.37977436, 0.37881807])
+
+And the velocity at links laong the center of the channel
+
+>>> linksAtCenter = grid.links_at_node[np.array(np.arange(600,660))][:-1,0]
+>>> grid['link']["surface_water__velocity"][linksAtCenter]
+array([ 0.45 , 0.59396861, 0.69574417, 0.75731793, 0.79261482,
+ 0.8123394 , 0.82700801, 0.84314289, 0.85567049, 0.85979031,
+ 0.85838515, 0.85658338, 0.85716582, 0.85944389, 0.86308335,
+ 0.86497462, 0.86463743, 0.86367451, 0.86314567, 0.86565779,
+ 0.86981418, 0.8726041 , 0.87338014, 0.87403681, 0.87494614,
+ 0.87549414, 0.87571505, 0.87786193, 0.87953196, 0.8794075 ,
+ 0.88030366, 0.88296806, 0.88376619, 0.88519116, 0.88610576,
+ 0.88744898, 0.88909536, 0.8898242 , 0.89135825, 0.89318975,
+ 0.89358086, 0.89405077, 0.89630327, 0.89764043, 0.8976985 ,
+ 0.89895516, 0.90065785, 0.90219923, 0.90227029, 0.90168375,
+ 0.90460332, 0.90538973, 0.9036968 , 0.90546049, 0.90666316,
+ 0.90682276, 0.90754984, 0.90761995, 0.91305765])
+
+"""
+
+import numpy as np
+import scipy as sp
+
+from landlab import Component, FieldError
+
+
+class river_flow_dynamics(Component):
+ """Simulate surface fluid flow based on Casulli and Cheng (1992).
+
+ Landlab component that simulates surface fluid flow using the Casulli and Cheng (1992)
+ approximations of the 2D shallow water equations.
+
+ This components calculates water depth and velocity across the raster grid, given
+ a certain input discharge.
+
+ References
+ ----------
+ **Required Software Citation(s) Specific to this Component**
+
+ None Listed
+
+ **Additional References**
+
+ Casulli, V., Cheng, R.T. (1992). “Semi-implicit finite difference methods for
+ three-dimensional shallow water flow”. International Journal for Numerical Methods
+ in Fluids. 15: 629-648.
+ https://doi.org/10.1002/fld.1650150602
+
+ """
+
+ _name = "river_flow_dynamics"
+
+ _unit_agnostic = False
+
+ _info = {
+ "surface_water__depth": {
+ "dtype": float,
+ "intent": "inout",
+ "optional": False,
+ "units": "m",
+ "mapping": "node",
+ "doc": "Depth of water on the surface",
+ },
+ "surface_water__velocity": {
+ "dtype": float,
+ "intent": "inout",
+ "optional": False,
+ "units": "m/s",
+ "mapping": "link",
+ "doc": "Speed of water flow above the surface",
+ },
+ "surface_water__elevation": {
+ "dtype": float,
+ "intent": "inout",
+ "optional": False,
+ "units": "m",
+ "mapping": "node",
+ "doc": "Water surface elevation at time N",
+ },
+ "topographic__elevation": {
+ "dtype": float,
+ "intent": "in",
+ "optional": False,
+ "units": "m",
+ "mapping": "node",
+ "doc": "Land surface topographic elevation",
+ },
+ }
+
+ def __init__(
+ self,
+ grid,
+ dt=0.01, # Sets the time step (s)
+ eddy_viscosity=1e-4, # ddy viscosity coefficient
+ mannings_n=0.012, # Manning's n
+ threshold_depth=0.01, # Sets the wet/dry threshold
+ theta=0.5, # Degree of 'implicitness' of the solution
+ fixed_entry_nodes=[], # Node IDs where flow enters the domain
+ fixed_entry_links=[], # Link IDs where flow enters the domain
+ entry_nodes_h_values=[], # Water depth at nodes where flow enters the domain
+ entry_links_vel_values=[], # Water velocity at links where flow enters the domain
+ pcg_tolerance=1e-05, # Preconditioned Conjugate Gradient convergence tolerance
+ pcg_max_iterations=None, # Preconditioned Conjugate Gradient max iterations
+ surface_water__elevation_at_N_1=None, # Sets the surface water elevation at previous time
+ surface_water__elevation_at_N_2=None, # Sets the surface water elevation at previous previous time
+ surface_water__velocity_at_N_1=None, # Sets the speed of water flow above the surface at time N-1
+ ):
+ """Simulate the vertical-averaged surface fluid flow
+
+ Landlab component that simulates surface fluid flow using the Casulli and Cheng (1992)
+ approximations of the 2D shallow water equations.
+
+ This components calculates water depth and velocity across the raster grid, given
+ a certain input discharge.
+
+ Parameters
+ ----------
+ grid : RasterModelGrid
+ A grid.
+ dt : float, optional
+ Time step in seconds. If not given, it is calculated from CFL condition.
+ eddy_viscosity : float, optional
+ Eddy viscosity coefficient. Default = 1e-4 :math:`m^2 / s`
+ mannings_n : float or array_like, optional
+ Manning's roughness coefficient. Default = 0.012 :math:`s / m^1/3`
+ threshold_depth : float, optional
+ Threshold at which a cell is considered wet. Default = 0.01 m
+ theta : float, optional
+ Degree of 'implicitness' of the solution, ranging between 0.5 and 1.0. Default = 0.5.
+ When it is equal to 0.5, the approximation is centered in time.
+ When it is equal to 1.0, the approximation is fully implicit.
+ fixed_entry_nodes : array_like or None, optional
+ Node IDs where flow enters the domain (Dirichlet boundary condition).
+ If not provided, the water already present in the domain is not renewed.
+ fixed_entry_links : array_like or None, optional
+ Link IDs where flow enters the domain (Dirichlet boundary condition).
+ If not provided, the water already present in the domain is not renewed.
+ entry_nodes_h_values : array_like, optional
+ Water depth values at nodes where flow enters the domain (Dirichlet boundary condition).
+ If not provided, the water already present in the domain is not renewed.
+ entry_links_vel_values : array_like, optional
+ Water velocity values at links where flow enters the domain (Dirichlet boundary condition).
+ If not provided, the water already present in the domain is not renewed.
+ pcg_tolerance : float, optional
+ Preconditioned Conjugate Gradient method argument. Tolerance for convergence.
+ Default = 1e-05
+ pcg_max_iterations : integer, optional
+ Preconditioned Conjugate Gradient method argument. Maximum number of iterations.
+ Iteration will stop after maxiter steps even if the specified tolerance has not been achieved.
+ Default = None
+ surface_water__elevation_at_N_1: float, optional
+ Water surface elevation at time N-1.
+ units: m, mapping: node
+ surface_water__elevation_at_N_2: float, optional
+ Water surface elevation at time N-2.
+ units: m, mapping: node
+ surface_water__velocity_at_N_1: float, optional
+ Speed of water flow above the surface at time N-1",
+ units: m/2, mapping: link
+ """
+ super().__init__(grid)
+
+ # User inputs
+ self._dt = dt
+ self._eddy_viscosity = eddy_viscosity
+ self._g = sp.constants.g
+ self._mannings_n = mannings_n
+ self._threshold_depth = threshold_depth
+ self._theta = theta
+ self._fixed_entry_nodes = fixed_entry_nodes
+ self._fixed_entry_links = fixed_entry_links
+ self._entry_nodes_h_values = entry_nodes_h_values
+ self._entry_links_vel_values = entry_links_vel_values
+ self._pcg_tolerance = pcg_tolerance
+ self._pcg_max_iterations = pcg_max_iterations
+
+ # Getting topography for further calculations
+ self._additional_z = 10 # To set the virtual reference elevation (z=0)
+ self._max_elevation = self._grid["node"]["topographic__elevation"].max()
+ self._z = self._grid["node"]["topographic__elevation"]
+ self._z = (self._z.max() + self._additional_z) - self._z
+
+ # Creating fields if they don't exist
+ try:
+ self._grid["node"]["surface_water__depth"] = grid.add_zeros(
+ "surface_water__depth",
+ at="node",
+ units=self._info["surface_water__depth"]["units"],
+ )
+ except FieldError:
+ self._grid["node"]["surface_water__depth"] = grid.at_node[
+ "surface_water__depth"
+ ]
+
+ try:
+ self._grid["link"]["surface_water__velocity"] = grid.add_zeros(
+ "surface_water__velocity",
+ at="link",
+ units=self._info["surface_water__velocity"]["units"],
+ )
+ except FieldError:
+ self._grid["link"]["surface_water__velocity"] = grid.at_link[
+ "surface_water__velocity"
+ ]
+
+ try:
+ self._grid["node"]["surface_water__elevation"] = grid.add_zeros(
+ "surface_water__elevation",
+ at="node",
+ units=self._info["surface_water__elevation"]["units"],
+ )
+ self._grid["node"]["surface_water__elevation"] = (
+ self._grid["node"]["surface_water__depth"] - self._z
+ )
+ except FieldError:
+ self._grid["node"]["surface_water__elevation"] = grid.at_node[
+ "surface_water__elevation"
+ ]
+
+ if surface_water__elevation_at_N_1 is None:
+ self._surface_water__elevation_at_N_1 = np.zeros(self._grid.number_of_nodes)
+ else:
+ if surface_water__elevation_at_N_1.size > 0:
+ if (
+ surface_water__elevation_at_N_1.shape[0]
+ == self._grid.number_of_nodes
+ ):
+ self._surface_water__elevation_at_N_1 = (
+ surface_water__elevation_at_N_1
+ )
+ else:
+ raise ValueError(
+ "surface_water__elevation_at_N_1 \
+ does not have the same dimensions of the grid's nodes"
+ )
+
+ if surface_water__elevation_at_N_2 is None:
+ self._surface_water__elevation_at_N_2 = np.zeros(self._grid.number_of_nodes)
+ else:
+ if surface_water__elevation_at_N_2.size > 0:
+ if (
+ surface_water__elevation_at_N_2.shape[0]
+ == self._grid.number_of_nodes
+ ):
+ self._surface_water__elevation_at_N_2 = (
+ surface_water__elevation_at_N_2
+ )
+ else:
+ raise ValueError(
+ "surface_water__elevation_at_N_2 \
+ does not have the same dimensions of the grid's nodes"
+ )
+
+ if surface_water__velocity_at_N_1 is None:
+ self._surface_water__velocity_at_N_1 = np.zeros(self._grid.number_of_links)
+ else:
+ if surface_water__velocity_at_N_1.size > 0:
+ if (
+ surface_water__velocity_at_N_1.shape[0]
+ == self._grid.number_of_links
+ ):
+ self._surface_water__velocity_at_N_1 = (
+ surface_water__velocity_at_N_1
+ )
+ else:
+ raise ValueError(
+ "surface_water__velocity_at_N_1 \
+ does not have the same dimensions of the grid's links"
+ )
+
+ # Assigning a class variable to the fields
+ self._h = self._grid["node"]["surface_water__depth"]
+ self._vel = self._grid["link"]["surface_water__velocity"]
+ self._vel_at_N_1 = self._surface_water__velocity_at_N_1
+ self._eta = self._grid["node"]["surface_water__elevation"] - (
+ self._max_elevation + self._additional_z
+ )
+ self._eta_at_N_1 = self._surface_water__elevation_at_N_1 - (
+ self._max_elevation + self._additional_z
+ )
+ self._eta_at_N_2 = self._surface_water__elevation_at_N_2 - (
+ self._max_elevation + self._additional_z
+ )
+
+ # Getting nodes and links at each edge
+ self._nodes_at_right_edge = self._grid.nodes_at_right_edge
+ self._nodes_at_top_edge = self._grid.nodes_at_top_edge
+ self._nodes_at_left_edge = self._grid.nodes_at_left_edge
+ self._nodes_at_bottom_edge = self._grid.nodes_at_bottom_edge
+
+ # links_at_right_edge = self._grid.links_at_node[self._nodes_at_right_edge][:, 2]
+ # links_at_top_edge = self._grid.links_at_node[self._nodes_at_top_edge][:, 3]
+ # links_at_left_edge = self._grid.links_at_node[self._nodes_at_left_edge][:, 0]
+ # links_at_bottom_edge = self._grid.links_at_node[self._nodes_at_bottom_edge][:, 1]
+
+ self._nx = len(self._nodes_at_top_edge)
+ self._ny = len(self._nodes_at_left_edge)
+
+ self._dx = self._grid.dx
+ self._dy = self._grid.dy
+
+ # Open boundary conditions (water can leave the domain at everywhere, only limited by topography)
+ self._grid.status_at_node[
+ self._nodes_at_left_edge
+ ] = self._grid.BC_NODE_IS_FIXED_VALUE
+ self._grid.status_at_node[
+ self._nodes_at_right_edge
+ ] = self._grid.BC_NODE_IS_FIXED_VALUE
+ self._grid.status_at_node[
+ self._nodes_at_bottom_edge
+ ] = self._grid.BC_NODE_IS_FIXED_VALUE
+ self._grid.status_at_node[
+ self._nodes_at_top_edge
+ ] = self._grid.BC_NODE_IS_FIXED_VALUE
+
+ # Identifying node and link ids for later use.
+ self._core_nodes = self._grid.core_nodes
+ self._corner_nodes = self._grid.corner_nodes
+ self._active_links = self._grid.active_links
+ self._horizontal_links = self._grid.horizontal_links
+ self._vertical_links = self._grid.vertical_links
+ self._number_of_nodes = self._grid.number_of_nodes
+
+ self._adjacent_nodes_at_corner_nodes = np.array(
+ [
+ [
+ self._nodes_at_top_edge[-2],
+ self._nodes_at_right_edge[-2],
+ ], # Top right
+ [self._nodes_at_top_edge[1], self._nodes_at_left_edge[-2]], # Top left
+ [
+ self._nodes_at_left_edge[1],
+ self._nodes_at_bottom_edge[1],
+ ], # Bottom left
+ [
+ self._nodes_at_right_edge[1],
+ self._nodes_at_bottom_edge[-2],
+ ], # Bottom right
+ ]
+ )
+
+ # Updating open boundary nodes/links
+ self._open_boundary_nodes = self._grid.boundary_nodes
+ self._open_boundary_links = np.unique(
+ self._grid.links_at_node[self._open_boundary_nodes]
+ )
+
+ self._open_boundary_nodes = np.setdiff1d(
+ self._open_boundary_nodes, self._fixed_entry_nodes
+ )
+ self._open_boundary_links = np.setdiff1d(
+ self._open_boundary_links, self._fixed_entry_links
+ )
+
+ self._fixed_corner_nodes = np.setdiff1d(
+ self._corner_nodes, self._open_boundary_nodes
+ )
+ self._open_corner_nodes = np.setdiff1d(
+ self._corner_nodes, self._fixed_corner_nodes
+ )
+
+ self._open_boundary_nodes = np.setdiff1d(
+ self._open_boundary_nodes, self._open_corner_nodes
+ )
+
+ # Using fixed entry nodes/links only when they exist
+ if len(self._fixed_entry_nodes) > 0:
+ self._fixed_nodes_exist = True
+ else:
+ self._fixed_nodes_exist = False
+
+ if len(self._fixed_entry_links) > 0:
+ self._fixed_links_exist = True
+ else:
+ self._fixed_links_exist = False
+
+ # Updating grid fixed values according to the user input
+ if self._fixed_nodes_exist is True:
+ self._h[self._fixed_entry_nodes] = entry_nodes_h_values
+ self._eta[self._fixed_entry_nodes] = (
+ entry_nodes_h_values - self._z[self._fixed_entry_nodes]
+ )
+ if self._fixed_links_exist is True:
+ self._vel[self._fixed_entry_links] = entry_links_vel_values
+
+ # Mapping node values at links
+ self._z_at_links = self._grid.map_mean_of_link_nodes_to_link(self._z)
+ self._h_at_links = self._grid.map_mean_of_link_nodes_to_link(self._h)
+ self._eta_at_links = self._h_at_links - self._z_at_links
+
+ # Passing values to the time step N
+ self._h_at_N = self._h.copy()
+ self._h_at_N_at_links = self._grid.map_mean_of_link_nodes_to_link(self._h_at_N)
+
+ self._vel_at_N = self._vel.copy()
+ self._eta_at_N = self._eta.copy()
+
+ # Boolean for wet nodes/links
+ self._wet_nodes = np.where(self._h_at_N >= self._threshold_depth, True, False)
+ self._wet_links = np.where(
+ self._h_at_N_at_links >= self._threshold_depth, True, False
+ )
+
+ # Defining some functions
+ def find_nearest_link(self, x_coordinates, y_coordinates, objective_links="all"):
+ """Link nearest a point.
+
+ Find the index to the link nearest the given x, y coordinates.
+ Returns the indices of the links nearest the given coordinates.
+
+ """
+ # Defining the set of links that are going to be used
+ if objective_links == "all":
+ objective_links = np.arange(self._grid.number_of_links)
+ elif objective_links == "horizontal":
+ objective_links = self._horizontal_links
+ elif objective_links == "vertical":
+ objective_links = self._vertical_links
+ # if (objective_links == "all") END
+
+ # Coordinates of all the RasterModelGrid links
+ x_of_objective_links = np.unique(self._grid.xy_of_link[objective_links][:, 0])
+ y_of_objective_links = np.unique(self._grid.xy_of_link[objective_links][:, 1])
+
+ # Getting the closest link-coordinate to the exit point
+ tempCalc1 = np.repeat(x_coordinates, len(x_of_objective_links)).reshape(
+ len(x_coordinates), len(x_of_objective_links)
+ )
+ tempCalc2 = np.tile(x_of_objective_links, len(x_coordinates)).reshape(
+ len(x_coordinates), len(x_of_objective_links)
+ )
+ indices = abs(tempCalc2 - tempCalc1).argmin(axis=1)
+ nearest_x = x_of_objective_links[indices]
+
+ tempCalc1 = np.repeat(y_coordinates, len(y_of_objective_links)).reshape(
+ len(y_coordinates), len(y_of_objective_links)
+ )
+ tempCalc2 = np.tile(y_of_objective_links, len(y_coordinates)).reshape(
+ len(y_coordinates), len(y_of_objective_links)
+ )
+ indices = abs(tempCalc2 - tempCalc1).argmin(axis=1)
+ nearest_y = y_of_objective_links[indices]
+
+ # Getting the closest link to link
+ tempCalc1 = np.repeat(
+ nearest_x, len(self._grid.xy_of_link[objective_links][:, 0])
+ ).reshape(len(nearest_x), len(self._grid.xy_of_link[objective_links][:, 0]))
+ tempCalc2 = np.tile(
+ self._grid.xy_of_link[objective_links][:, 0], len(x_coordinates)
+ ).reshape(len(x_coordinates), len(self._grid.xy_of_link[objective_links][:, 0]))
+ tempB1 = tempCalc1 == tempCalc2
+
+ tempCalc1 = np.repeat(
+ nearest_y, len(self._grid.xy_of_link[objective_links][:, 1])
+ ).reshape(len(nearest_y), len(self._grid.xy_of_link[objective_links][:, 1]))
+ tempCalc2 = np.tile(
+ self._grid.xy_of_link[objective_links][:, 1], len(y_coordinates)
+ ).reshape(len(y_coordinates), len(self._grid.xy_of_link[objective_links][:, 1]))
+ tempB2 = tempCalc1 == tempCalc2
+
+ tempCalc3 = (
+ np.repeat(objective_links, len(x_coordinates))
+ .reshape((len(objective_links), len(y_coordinates)))
+ .T
+ )
+ nearest_link = tempCalc3[tempB1 * tempB2]
+
+ return nearest_link.astype(int)
+
+ def find_adjacent_links_at_link(self, current_link, objective_links="horizontal"):
+ """Get adjacent links to the link.
+
+ This function finds the links at right, above, left and below the given link.
+ Similar purpose to the "adjacent_nodes_at_node" function.
+ Return the adjacent links in as many rows as given links.
+ Link IDs are returned as columns in clock-wise order starting from East (E, N, W, S).
+
+ """
+
+ # Defining the set of links that are going to be used
+ if objective_links == "horizontal":
+ objective_links = self._horizontal_links
+ reshape_pair = (self._ny, self._nx - 1)
+ elif objective_links == "vertical":
+ objective_links = self._vertical_links
+ reshape_pair = (self._ny - 1, self._nx)
+ # if (objective_links == "horizontal") END
+
+ # Coordinates of the current link
+ x_of_current_links = self._grid.xy_of_link[current_link][:, 0]
+ y_of_current_links = self._grid.xy_of_link[current_link][:, 1]
+
+ # Coordinates of all the RasterModelGrid links
+ x_of_objective_links = np.unique(self._grid.xy_of_link[objective_links][:, 0])
+ y_of_objective_links = np.unique(self._grid.xy_of_link[objective_links][:, 1])
+
+ # Getting links that share the same y-coordinates
+ # The following matrices are built to be compared to each other.
+ # tempCalc1 repeats "y_of_current_links" for every x-coordinate in "objective_links": cols = "y_of_current_links"
+ # tempCalc2 repeats "y_of_objective_links" for every x-coordinate in the "current_link": rows = "y_of_objective_links"
+ # tempCalc3 give us the index to extract all the objective links that are located in the same row than the current link: rows = [0, 1, 2, ...]
+ tempCalc1 = np.repeat(y_of_current_links, len(y_of_objective_links)).reshape(
+ len(y_of_current_links), len(y_of_objective_links)
+ )
+ tempCalc2 = np.tile(y_of_objective_links, len(y_of_current_links)).reshape(
+ len(y_of_current_links), len(y_of_objective_links)
+ )
+ tempCalc3 = (
+ np.repeat(np.arange(len(y_of_objective_links)), len(y_of_current_links))
+ .reshape(len(y_of_objective_links), len(y_of_current_links))
+ .T
+ )
+
+ indices = tempCalc3[tempCalc1 == tempCalc2]
+ links_at_same_rows = objective_links.reshape(reshape_pair)[indices, :]
+ links_at_same_rows = np.append(
+ np.array([-np.ones_like(current_link)]).T, links_at_same_rows, axis=1
+ )
+ links_at_same_rows = np.append(
+ links_at_same_rows, np.array([-np.ones_like(current_link)]).T, axis=1
+ )
+
+ # Getting links that share the same x-coordinates
+ # The following matrices are built to be compared to each other.
+ # tempCalc1 repeats "x_of_current_links" for every x-coordinate in "objective_links": cols = "x_of_current_links"
+ # tempCalc2 repeats "x_of_objective_links" for every x-coordinate in the "current_link": rows = "x_of_objective_links"
+ # tempCalc3 give us the index to extract all the objective links that are located in the same row than the current link: rows = [0, 1, 2, ...]
+ tempCalc1 = np.repeat(x_of_current_links, len(x_of_objective_links)).reshape(
+ len(x_of_current_links), len(x_of_objective_links)
+ )
+ tempCalc2 = np.tile(x_of_objective_links, len(x_of_current_links)).reshape(
+ len(x_of_current_links), len(x_of_objective_links)
+ )
+ tempCalc3 = (
+ np.repeat(np.arange(len(x_of_objective_links)), len(x_of_current_links))
+ .reshape(len(x_of_objective_links), len(x_of_current_links))
+ .T
+ )
+
+ indices = tempCalc3[tempCalc1 == tempCalc2]
+ links_at_same_cols = objective_links.reshape(reshape_pair)[:, indices].T
+ links_at_same_cols = np.append(
+ np.array([-np.ones_like(current_link)]).T, links_at_same_cols, axis=1
+ )
+ links_at_same_cols = np.append(
+ links_at_same_cols, np.array([-np.ones_like(current_link)]).T, axis=1
+ )
+
+ # Extracing the adjacent links to current link (E,N,W,S)
+ adjacent_links_at_link = np.zeros((current_link.shape[0], 4))
+
+ # Rows (E,W)
+ tempCalc1 = np.repeat(current_link, links_at_same_rows.shape[1]).reshape(
+ current_link.shape[0], links_at_same_rows.shape[1]
+ )
+ tempCalc2 = (
+ np.repeat(np.arange(links_at_same_rows.shape[1]), current_link.shape[0])
+ .reshape(links_at_same_rows.shape[1], current_link.shape[0])
+ .T
+ )
+ tempCalc3 = tempCalc2[tempCalc1 == links_at_same_rows]
+
+ adjacent_links_at_link[:, 0] = links_at_same_rows[
+ (range(links_at_same_rows.shape[0])), (tempCalc3 + 1)
+ ]
+ adjacent_links_at_link[:, 2] = links_at_same_rows[
+ (range(links_at_same_rows.shape[0])), (tempCalc3 - 1)
+ ]
+
+ # Cols (N,S)
+ tempCalc1 = np.repeat(current_link, links_at_same_cols.shape[1]).reshape(
+ current_link.shape[0], links_at_same_cols.shape[1]
+ )
+ tempCalc2 = (
+ np.repeat(np.arange(links_at_same_cols.shape[1]), current_link.shape[0])
+ .reshape(links_at_same_cols.shape[1], current_link.shape[0])
+ .T
+ )
+ tempCalc3 = tempCalc2[tempCalc1 == links_at_same_cols]
+
+ adjacent_links_at_link[:, 1] = links_at_same_cols[
+ (range(links_at_same_cols.shape[0])), (tempCalc3 + 1)
+ ]
+ adjacent_links_at_link[:, 3] = links_at_same_cols[
+ (range(links_at_same_cols.shape[0])), (tempCalc3 - 1)
+ ]
+
+ return adjacent_links_at_link.astype(int)
+
+ def path_line_tracing(self):
+ """ " Path line tracing algorithm.
+
+ This function implements the semi-analytical path line tracing method of Pollock (1988).
+
+ The semi-analytical path line tracing method was developed for particle tracking in ground
+ water flow models. The assumption that each directional velocity component varies linearly
+ in its coordinate directions within each computational volume or cell underlies the method.
+ Linear variation allows the derivation of an analytical expression for the path line of a
+ particle across a volume.
+
+ Given an initial point located at each volume faces of the domain, particle trayectories are
+ traced backwards on time. Then, this function returns the departure point of the particle at
+ the beginning of the time step.
+
+ """
+
+ # Calculating the partial time-step TAUx, TAUy, dt - sum_partial_times
+ sum_partial_times = np.zeros_like(self._u_vel_of_particle)
+ remaining_time = self._dt - sum_partial_times
+ keep_tracing = np.where(remaining_time > 0, True, False)
+
+ while np.any(remaining_time > 0):
+ # Using the previous exit point as the new entry point
+ self._x_of_particle = np.where(keep_tracing, self._x_at_exit_point, self._x_of_particle)
+ self._y_of_particle = np.where(
+ keep_tracing , self._y_at_exit_point, self._y_of_particle
+ )
+
+ # Checking if the particles departs (backwards) from a link position (True)
+ tempBx = np.isin(
+ self._x_of_particle, self._grid.xy_of_link[self._active_links][:, 0]
+ ) # Particles located on horizontal-links/vertical-faces
+ tempBy = np.isin(
+ self._y_of_particle, self._grid.xy_of_link[self._active_links][:, 1]
+ ) # Particles located on vertical-links/horizontal-faces
+
+ # True, particles depart from link positions.
+ # False, particles depart from random locations inside a cell
+ tempBxy = tempBx + tempBy
+
+ # Getting surrounding links for particles located inside a cell
+ tempCalc1 = np.append(
+ np.array([self._x_of_particle]), np.array([self._y_of_particle]), axis=0
+ )
+ tempCalc2 = self._grid.find_nearest_node(tempCalc1, mode="raise")
+ temp_links_from_node = self._grid.links_at_node[tempCalc2]
+ nodes_from_particle = tempCalc2
+
+ # Getting surrounding links for particles located at link positions
+ tempCalc1 = np.where(
+ self._u_vel_of_particle >= 0,
+ np.array([self._x_of_particle]) - self._dx / 10,
+ np.array([self._x_of_particle]) + self._dx / 10,
+ )
+ tempCalc2 = np.where(
+ self._v_vel_of_particle >= 0,
+ np.array([self._y_of_particle]) - self._dy / 10,
+ np.array([self._y_of_particle]) + self._dy / 10,
+ )
+ tempCalc3 = np.append(tempCalc1, tempCalc2, axis=0)
+ tempCalc4 = self._grid.find_nearest_node(tempCalc3, mode="raise")
+ temp_links_from_link = self._grid.links_at_node[tempCalc4]
+ nodes_from_particle = np.where(
+ tempBxy , tempCalc4, nodes_from_particle
+ )
+
+ # Getting links around particle
+ tempBxy = np.tile(tempBxy, 4).reshape(4, len(tempBxy)).T
+ links_at_particle = np.where(
+ tempBxy , temp_links_from_link, temp_links_from_node
+ )
+
+ # Defining links based on velocity direction
+ link_at_x2 = np.where(
+ self._u_vel_of_particle >= 0,
+ links_at_particle[:, 0],
+ links_at_particle[:, 2],
+ )
+ link_at_x1 = np.where(
+ self._u_vel_of_particle >= 0,
+ links_at_particle[:, 2],
+ links_at_particle[:, 0],
+ )
+ link_at_y2 = np.where(
+ self._v_vel_of_particle >= 0,
+ links_at_particle[:, 1],
+ links_at_particle[:, 3],
+ )
+ link_at_y1 = np.where(
+ self._v_vel_of_particle >= 0,
+ links_at_particle[:, 3],
+ links_at_particle[:, 1],
+ )
+
+ x_at_x2 = np.where(
+ self._u_vel_of_particle >= 0,
+ self._grid.x_of_node[nodes_from_particle] + self._dx / 2,
+ self._grid.x_of_node[nodes_from_particle] - self._dx / 2,
+ )
+ x_at_x1 = np.where(
+ self._u_vel_of_particle >= 0,
+ self._grid.x_of_node[nodes_from_particle] - self._dx / 2,
+ self._grid.x_of_node[nodes_from_particle] + self._dx / 2,
+ )
+ y_at_y2 = np.where(
+ self._v_vel_of_particle >= 0,
+ self._grid.y_of_node[nodes_from_particle] + self._dy / 2,
+ self._grid.y_of_node[nodes_from_particle] - self._dy / 2,
+ )
+ y_at_y1 = np.where(
+ self._v_vel_of_particle >= 0,
+ self._grid.y_of_node[nodes_from_particle] - self._dy / 2,
+ self._grid.y_of_node[nodes_from_particle] + self._dy / 2,
+ )
+
+ # Getting velocity around the particle
+ u_vel_at_x2 = np.where(
+ link_at_x2 >= 0, self._vel_at_N[link_at_x2], self._vel_at_N[link_at_x1]
+ )
+ u_vel_at_x1 = np.where(
+ link_at_x1 >= 0, self._vel_at_N[link_at_x1], self._vel_at_N[link_at_x2]
+ )
+ v_vel_at_y2 = np.where(
+ link_at_y2 >= 0, self._vel_at_N[link_at_y2], self._vel_at_N[link_at_y1]
+ )
+ v_vel_at_y1 = np.where(
+ link_at_y1 >= 0, self._vel_at_N[link_at_y1], self._vel_at_N[link_at_y2]
+ )
+
+ # Calculating gradients for path line tracing
+ gradient_x_direction = (u_vel_at_x2 - u_vel_at_x1) / self._dx
+ gradient_y_direction = (v_vel_at_y2 - v_vel_at_y1) / self._dy
+
+ # Calculating entry velocity for each particle
+ self._u_vel_of_particle = u_vel_at_x2 - gradient_x_direction * (
+ x_at_x2 - self._x_of_particle
+ )
+ self._v_vel_of_particle = v_vel_at_y2 - gradient_y_direction * (
+ y_at_y2 - self._y_of_particle
+ )
+ self._u_vel_of_particle = np.where(
+ self._u_vel_of_particle < 1e-10, 0, self._u_vel_of_particle
+ )
+ self._v_vel_of_particle = np.where(
+ self._v_vel_of_particle < 1e-10, 0, self._v_vel_of_particle
+ )
+
+ ### Calculation accoss x-direction
+ # Avoiding divisions by zero
+ tempCalc1 = np.where(
+ self._u_vel_of_particle == 0, 9999, self._u_vel_of_particle
+ )
+ tempCalc2 = np.where(u_vel_at_x1 == 0, 9999, u_vel_at_x1)
+ tempCalc3 = np.where(gradient_x_direction == 0, 9999, gradient_x_direction)
+ TAUx = (1 / tempCalc3) * np.log(abs(tempCalc1 / tempCalc2))
+
+ # Calculation when gradient is equal to zero
+ tempCalc4 = abs((self._x_of_particle - x_at_x1) / tempCalc2)
+ TAUx = np.where(gradient_x_direction == 0, tempCalc4, TAUx)
+
+ # Calculation when:
+ # a) Uxp/Ux1 = 1,
+ # b) Uxp,Vyp = 0,
+ # c) Ux1,Vy1 = 0, and
+ # d) Uxp/Ux1, Vxp/Vy1 = -1
+ tempCalc5 = self._u_vel_of_particle / tempCalc2
+ TAUx = np.where(tempCalc5 == 1, tempCalc4, TAUx)
+ TAUx = np.where(self._u_vel_of_particle == 0, remaining_time, TAUx)
+ TAUx = np.where(u_vel_at_x1 == 0, remaining_time, TAUx)
+ TAUx = np.where(tempCalc5 < 0, remaining_time, TAUx)
+ TAUx = np.where(TAUx > self._dt, self._dt, TAUx)
+ TAUx = np.where(TAUx < 0, 0, TAUx)
+
+ ### Calculation across y-direction
+ # Avoiding divisions by zero
+ tempCalc1 = np.where(
+ self._v_vel_of_particle == 0, 9999, self._v_vel_of_particle
+ )
+ tempCalc2 = np.where(v_vel_at_y1 == 0, 9999, v_vel_at_y1)
+ tempCalc3 = np.where(gradient_y_direction == 0, 9999, gradient_y_direction)
+ TAUy = (1 / tempCalc3) * np.log(abs(tempCalc1 / tempCalc2))
+
+ # Calculation when gradient is equal to zero
+ tempCalc4 = abs((self._y_of_particle - y_at_y1) / tempCalc2)
+ TAUy = np.where(gradient_y_direction == 0, tempCalc4, TAUy)
+
+ # Calculation when
+ # a) Vyp/Vy1 = 1,
+ # b) Uxp,Vyp = 0,
+ # c) Ux1,Vy1 = 0, and
+ # d) Uxp/Ux1, Vxp/Vy1 = -1
+ tempCalc5 = self._v_vel_of_particle / tempCalc2
+ TAUy = np.where(tempCalc5 == 1, tempCalc4, TAUy)
+ TAUy = np.where(self._v_vel_of_particle == 0, remaining_time, TAUy)
+ TAUy = np.where(v_vel_at_y1 == 0, remaining_time, TAUy)
+ TAUy = np.where(tempCalc5 < 0, remaining_time, TAUy)
+ TAUy = np.where(TAUy > self._dt, self._dt, TAUy)
+ TAUy = np.where(TAUy < 0, 0, TAUy)
+
+ # Obtaining TAU = min(TAUx, TAUy, (dt - sum_partial_times))
+ TAUx = abs(TAUx)
+ TAUy = abs(TAUy)
+ TAU = np.array((TAUx, TAUy, remaining_time)).min(axis=0)
+ # TAU = np.where(TAU < 1e-10, 0, TAU)
+
+ # Calculating exit point Xe, Ye
+ tempCalc1 = np.where(gradient_x_direction == 0, 9999, gradient_x_direction)
+ tempCalc2 = np.where(gradient_y_direction == 0, 9999, gradient_y_direction)
+
+ # Exit point Xe (tempCalc3) and Ye (tempCalc4)
+ tempCalc3 = x_at_x2 - (1 / tempCalc1) * (
+ u_vel_at_x2
+ - self._u_vel_of_particle / np.exp(gradient_x_direction * TAU)
+ )
+ tempCalc4 = y_at_y2 - (1 / tempCalc2) * (
+ v_vel_at_y2
+ - self._v_vel_of_particle / np.exp(gradient_y_direction * TAU)
+ )
+
+ tempCalc3 = np.where(
+ gradient_x_direction == 0,
+ self._x_of_particle - u_vel_at_x2 * TAU,
+ tempCalc3,
+ )
+ tempCalc4 = np.where(
+ gradient_y_direction == 0,
+ self._y_of_particle - v_vel_at_y2 * TAU,
+ tempCalc4,
+ )
+
+ tempCalc3 = np.where(
+ self._u_vel_of_particle == 0, self._x_of_particle, tempCalc3
+ )
+ tempCalc4 = np.where(
+ self._v_vel_of_particle == 0, self._y_of_particle, tempCalc4
+ )
+
+ self._x_at_exit_point = np.where(
+ keep_tracing , tempCalc3, self._x_at_exit_point
+ )
+ self._y_at_exit_point = np.where(
+ keep_tracing , tempCalc4, self._y_at_exit_point
+ )
+
+ # Updating sum of partial time-steps, TAU
+ sum_partial_times = np.where(
+ keep_tracing , sum_partial_times + TAU, self._dt
+ )
+
+ # Checking remaining_time == 0 (dt = sum_partial_times)
+ remaining_time = np.where(
+ remaining_time == 0, 0, self._dt - sum_partial_times
+ )
+
+ # Correcting entry velocity
+ tempCalc1 = np.where(self._u_vel_of_particle == 0, 1, 0)
+ tempCalc2 = np.where(self._v_vel_of_particle == 0, 1, 0)
+ remaining_time = np.where((tempCalc1 * tempCalc2) == 1, 0, remaining_time)
+
+ # Correction for static particles
+ remaining_time = np.where(
+ abs(self._x_of_particle - self._x_at_exit_point) < 1e-7,
+ 0,
+ remaining_time,
+ )
+ remaining_time = np.where(
+ abs(self._y_of_particle - self._y_at_exit_point) < 1e-7,
+ 0,
+ remaining_time,
+ )
+
+ # Stop tracing if a particle hits the boundary
+ # Keep tracing for all particles
+ tempCalc1 = np.repeat(True, len(keep_tracing))
+ # Particle hits the left edge?
+ tempCalc2 = np.isin(
+ self._x_at_exit_point, self._grid.x_of_node[self._nodes_at_left_edge]
+ )
+ # If above True, stop tracing for that particle
+ tempCalc1 = np.where(tempCalc2 , False, tempCalc1)
+ # Particle hits the right edge?
+ tempCalc2 = np.isin(
+ self._x_at_exit_point, self._grid.x_of_node[self._nodes_at_right_edge]
+ )
+ # If above True, stop tracing for that particle
+ tempCalc1 = np.where(tempCalc2 , False, tempCalc1)
+ # Particle hits the top edge?
+ tempCalc2 = np.isin(
+ self._y_at_exit_point, self._grid.y_of_node[self._nodes_at_top_edge]
+ )
+ # If above True, stop tracing for that particle
+ tempCalc1 = np.where(tempCalc2 , False, tempCalc1)
+ # Particle hits the bottom edge?
+ tempCalc2 = np.isin(
+ self._y_at_exit_point, self._grid.y_of_node[self._nodes_at_bottom_edge]
+ )
+ # If above True, stop tracing for that particle
+ tempCalc1 = np.where(tempCalc2 , False, tempCalc1)
+ # Where particles reach the boundary, remaining time is equal to zero
+ #remaining_time = np.where(not tempCalc1, 0, remaining_time)
+ remaining_time = np.where(~tempCalc1, 0, remaining_time)
+
+
+ # Updating on particles that need to traced backwards
+ keep_tracing = np.where(remaining_time > 0, True, False)
+
+ # WHILE "np.any(remaining_time > 0)" END
+ # DEF "path_line_tracing(self)" END
+
+ def run_one_step(self):
+ """Calculate water depth and water velocity for a time period dt."""
+
+ # Getting velocity as U,V components
+ self._u_vel = self._vel_at_N[self._horizontal_links]
+ self._v_vel = self._vel_at_N[self._vertical_links]
+
+ # Calculating Chezy coefficient
+ self._chezy_at_nodes = self._h_at_N ** (1 / 6) / self._mannings_n
+ self._chezy_at_links = self._h_at_N_at_links ** (1 / 6) / self._mannings_n
+
+ # Computing V-velocity (vertical links) at U-velocity positions (horizontal links)
+ tempCalc1 = self._grid.map_mean_of_horizontal_links_to_node(self._vel_at_N)
+ self._u_vel_at_v_links = np.mean(
+ tempCalc1[self._grid.nodes_at_link[self._vertical_links]], axis=1
+ )
+
+ # Computing U-velocity (horizontal links) at V-velocity positions (vertical links)
+ tempCalc1 = self._grid.map_mean_of_vertical_links_to_node(self._vel_at_N)
+ self._v_vel_at_u_links = np.mean(
+ tempCalc1[self._grid.nodes_at_link[self._horizontal_links]], axis=1
+ )
+
+ """ Setting A-faces """
+ # Setting A-faces
+ self._a_links = np.zeros_like(self._vel_at_N)
+
+ # Setting dry links equal to 1 to avoid divisions by zero
+ tempCalc1 = np.where(self._wet_links , self._chezy_at_links, 1)
+
+ # Computing A-faces
+ self._a_links[self._horizontal_links] = (
+ self._h_at_N_at_links[self._horizontal_links]
+ + self._g
+ * self._dt
+ * (
+ self._vel_at_N[self._horizontal_links] ** 2
+ + self._v_vel_at_u_links**2
+ )
+ ** (1 / 2)
+ / tempCalc1[self._horizontal_links]
+ )
+ self._a_links[self._vertical_links] = (
+ self._h_at_N_at_links[self._vertical_links]
+ + self._g
+ * self._dt
+ * (self._vel_at_N[self._vertical_links] ** 2 + self._u_vel_at_v_links**2)
+ ** (1 / 2)
+ / tempCalc1[self._vertical_links]
+ )
+
+ # Using only wet-link values, and setting dry links equal to 1 to avoid divisions by zero
+ self._a_links = np.where(self._wet_links , self._a_links, 1)
+
+ """ Path line tracing
+ U-velocity, x-direction, horizontal links
+ """
+ # Getting the initial particle location at each volume faces
+ tempB1 = [i in self._horizontal_links for i in self._active_links]
+ self._x_of_particle = self._grid.xy_of_link[:, 0][self._active_links][tempB1]
+ self._y_of_particle = self._grid.xy_of_link[:, 1][self._active_links][tempB1]
+
+ # Getting the initial particle velocity
+ tempB2 = [i in self._active_links for i in self._horizontal_links]
+ self._u_vel_of_particle = self._u_vel[tempB2]
+ self._v_vel_of_particle = self._v_vel_at_u_links[tempB2]
+
+ # Getting a first 'exit' point to begin the loop
+ self._x_at_exit_point = self._x_of_particle
+ self._y_at_exit_point = self._y_of_particle
+
+ # Calculating path line backwards on time
+ self.path_line_tracing()
+
+ """ Bicuatradic interpolation
+ U-velocity, x-direction, around (p,q) location
+ """
+ self._UsL = np.zeros_like(self._u_vel)
+
+ # Getting V-velocity at U-links
+ temp_Vvel = np.zeros_like(self._vel_at_N)
+ temp_Vvel[self._horizontal_links] = self._v_vel_at_u_links
+
+ # Getting links around the particle and defining downstream direction based on velocity
+ nearest_link_to_particle = self.find_nearest_link(
+ self._x_at_exit_point, self._y_at_exit_point, objective_links="horizontal"
+ )
+ adjacent_links_to_particle = self.find_adjacent_links_at_link(
+ nearest_link_to_particle, objective_links="horizontal"
+ )
+
+ link_at_B2 = nearest_link_to_particle
+ link_at_A2 = adjacent_links_to_particle[:, 1] # 1: N, top
+ link_at_C2 = adjacent_links_to_particle[:, 3] # 3: S, bottom
+ link_at_A2 = np.where(temp_Vvel[link_at_B2] >= 0, link_at_A2, link_at_C2)
+ link_at_C2 = np.where(temp_Vvel[link_at_B2] >= 0, link_at_C2, link_at_A2)
+ # We avoid "-1" links close to the boundary
+ link_at_A2 = np.where(link_at_A2 >= 0, link_at_A2, link_at_B2)
+ link_at_C2 = np.where(link_at_C2 >= 0, link_at_C2, link_at_B2)
+
+ # Getting the surrounding links to every particle from closest link
+ # 0: E, Right
+ # 2: W, Left
+ link_at_A1 = self.find_adjacent_links_at_link(
+ link_at_A2, objective_links="horizontal"
+ )[:, 0]
+ link_at_A3 = self.find_adjacent_links_at_link(
+ link_at_A2, objective_links="horizontal"
+ )[:, 2]
+
+ # Selecting downstream link based on velocity direction
+ link_at_A1 = np.where(self._vel_at_N[link_at_B2] >= 0, link_at_A1, link_at_A3)
+ link_at_A3 = np.where(self._vel_at_N[link_at_B2] >= 0, link_at_A3, link_at_A1)
+
+ link_at_B1 = self.find_adjacent_links_at_link(
+ link_at_B2, objective_links="horizontal"
+ )[
+ :, 0
+ ] # 0: E, Right
+ link_at_B3 = self.find_adjacent_links_at_link(
+ link_at_B2, objective_links="horizontal"
+ )[
+ :, 2
+ ] # 2: W, Left
+
+ # Selecting downstream link based on velocity direction
+ link_at_B1 = np.where(self._vel_at_N[link_at_B2] >= 0, link_at_B1, link_at_B3)
+ link_at_B3 = np.where(self._vel_at_N[link_at_B2] >= 0, link_at_B3, link_at_B1)
+
+ link_at_C1 = self.find_adjacent_links_at_link(
+ link_at_C2, objective_links="horizontal"
+ )[
+ :, 0
+ ] # 0: E, Right
+ link_at_C3 = self.find_adjacent_links_at_link(
+ link_at_C2, objective_links="horizontal"
+ )[
+ :, 2
+ ] # 2: W, Left
+
+ # Selecting downstream link based on velocity direction
+ link_at_C1 = np.where(self._vel_at_N[link_at_B2] >= 0, link_at_C1, link_at_C3)
+ link_at_C3 = np.where(self._vel_at_N[link_at_B2] >= 0, link_at_C3, link_at_C1)
+
+ # Getting velocity around the particle
+ vel_at_A1 = np.where(
+ link_at_A1 >= 0, self._vel_at_N[link_at_A1], self._vel_at_N[link_at_A2]
+ )
+ vel_at_A2 = self._vel_at_N[link_at_A2]
+ vel_at_A3 = np.where(
+ link_at_A3 >= 0, self._vel_at_N[link_at_A3], self._vel_at_N[link_at_A2]
+ )
+
+ vel_at_B1 = np.where(
+ link_at_B1 >= 0, self._vel_at_N[link_at_B1], self._vel_at_N[link_at_B2]
+ )
+ vel_at_B2 = self._vel_at_N[link_at_B2]
+ vel_at_B3 = np.where(
+ link_at_B3 >= 0, self._vel_at_N[link_at_B3], self._vel_at_N[link_at_B2]
+ )
+
+ vel_at_C1 = np.where(
+ link_at_C1 >= 0, self._vel_at_N[link_at_C1], self._vel_at_N[link_at_C2]
+ )
+ vel_at_C2 = self._vel_at_N[link_at_C2]
+ vel_at_C3 = np.where(
+ link_at_C3 >= 0, self._vel_at_N[link_at_C3], self._vel_at_N[link_at_C2]
+ )
+
+ # Getting coordinates around the particle
+ x_at_2 = self._grid.xy_of_link[link_at_B2][:, 0]
+ x_at_1 = np.where(
+ self._vel_at_N[link_at_B2] >= 0, x_at_2 + self._dx, x_at_2 - self._dx
+ )
+ x_at_3 = np.where(
+ self._vel_at_N[link_at_B2] >= 0, x_at_2 - self._dx, x_at_2 + self._dx
+ )
+
+ y_at_B = self._grid.xy_of_link[link_at_B2][:, 1]
+ y_at_A = np.where(
+ temp_Vvel[link_at_B2] >= 0, y_at_B + self._dy, y_at_B - self._dy
+ )
+ y_at_C = np.where(
+ temp_Vvel[link_at_B2] >= 0, y_at_B - self._dy, y_at_B + self._dy
+ )
+
+ # Calculating the weights W(i,j) for k around x-direction
+ W1 = (
+ (self._x_at_exit_point - x_at_2)
+ * (self._x_at_exit_point - x_at_3)
+ / ((x_at_1 - x_at_2) * (x_at_1 - x_at_3))
+ )
+ W2 = (
+ (self._x_at_exit_point - x_at_1)
+ * (self._x_at_exit_point - x_at_3)
+ / ((x_at_2 - x_at_1) * (x_at_2 - x_at_3))
+ )
+ W3 = (
+ (self._x_at_exit_point - x_at_1)
+ * (self._x_at_exit_point - x_at_2)
+ / ((x_at_3 - x_at_1) * (x_at_3 - x_at_2))
+ )
+
+ # Interpolation by row around 'x_at_exit_point'
+ A = W1 * vel_at_A1 + W2 * vel_at_A2 + W3 * vel_at_A3
+ B = W1 * vel_at_B1 + W2 * vel_at_B2 + W3 * vel_at_B3
+ C = W1 * vel_at_C1 + W2 * vel_at_C2 + W3 * vel_at_C3
+
+ # Calculating the weghts W(i,j) for l around y-direction
+ W1 = (
+ (self._y_at_exit_point - y_at_B)
+ * (self._y_at_exit_point - y_at_C)
+ / ((y_at_A - y_at_B) * (y_at_A - y_at_C))
+ )
+ W2 = (
+ (self._y_at_exit_point - y_at_A)
+ * (self._y_at_exit_point - y_at_C)
+ / ((y_at_B - y_at_A) * (y_at_B - y_at_C))
+ )
+ W3 = (
+ (self._y_at_exit_point - y_at_A)
+ * (self._y_at_exit_point - y_at_B)
+ / ((y_at_C - y_at_A) * (y_at_C - y_at_B))
+ )
+
+ # Calculating UsL by bicuadratic interpolation
+ self._UsL[tempB2] = W1 * A + W2 * B + W3 * C
+
+ """ Computing viscous terms
+ U-located particles
+ """
+ # Central difference scheme around x- and y- direction for U-located particles
+ self._Uvis = np.zeros_like(self._u_vel)
+
+ tempCalc1 = (
+ self._eddy_viscosity
+ * self._dt
+ * (vel_at_B3 - 2 * vel_at_B2 + vel_at_B1)
+ / (self._dx**2)
+ )
+ tempCalc2 = (
+ self._eddy_viscosity
+ * self._dt
+ * (vel_at_C2 - 2 * vel_at_B2 + vel_at_A2)
+ / (self._dy**2)
+ )
+
+ self._Uvis[tempB2] = tempCalc1 + tempCalc2
+
+ """ Path line tracing
+ V-velocity, y-direction, vertical links
+ """
+ # Getting the initial particle location at each volume faces
+ tempB1 = [j in self._vertical_links for j in self._active_links]
+ self._x_of_particle = self._grid.xy_of_link[:, 0][self._active_links][tempB1]
+ self._y_of_particle = self._grid.xy_of_link[:, 1][self._active_links][tempB1]
+
+ # Getting the initial particle velocity
+ tempB2 = [j in self._active_links for j in self._vertical_links]
+ self._v_vel_of_particle = self._v_vel[tempB2]
+ self._u_vel_of_particle = self._u_vel_at_v_links[tempB2]
+
+ # Getting a first 'exit' point to begin the loop
+ self._x_at_exit_point = self._x_of_particle
+ self._y_at_exit_point = self._y_of_particle
+
+ # Calculating path line backwards on time
+ self.path_line_tracing()
+
+ """ Bicuatradic interpolation
+ V-velocity, y-direction, around (p,q) location
+ """
+ self._VsL = np.zeros_like(self._v_vel)
+
+ # Getting V-velocity at U-links
+ temp_Uvel = np.zeros_like(self._vel_at_N)
+ temp_Uvel[self._vertical_links] = self._u_vel_at_v_links
+
+ # Getting links around the particle and defining downstream direction based on velocity
+ nearest_link_to_particle = self.find_nearest_link(
+ self._x_at_exit_point, self._y_at_exit_point, objective_links="vertical"
+ )
+ adjacent_links_to_particle = self.find_adjacent_links_at_link(
+ nearest_link_to_particle, objective_links="vertical"
+ )
+
+ link_at_B2 = nearest_link_to_particle
+ link_at_A2 = adjacent_links_to_particle[:, 1] # 1: N, top
+ link_at_C2 = adjacent_links_to_particle[:, 3] # 3: S, bottom
+ link_at_A2 = np.where(self._vel_at_N[link_at_B2] >= 0, link_at_A2, link_at_C2)
+ link_at_C2 = np.where(self._vel_at_N[link_at_B2] >= 0, link_at_C2, link_at_A2)
+ link_at_A2 = np.where(link_at_A2 >= 0, link_at_A2, link_at_B2)
+ link_at_C2 = np.where(link_at_C2 >= 0, link_at_C2, link_at_B2)
+ # We avoid "-1" links close to the boundary
+
+ # Getting the surrounding links to every particle from closest link
+ link_at_A1 = self.find_adjacent_links_at_link(
+ link_at_A2, objective_links="vertical"
+ )[
+ :, 0
+ ] # 0: E, Left
+ link_at_A3 = self.find_adjacent_links_at_link(
+ link_at_A2, objective_links="vertical"
+ )[
+ :, 2
+ ] # 2: W, Right
+
+ # Selecting downstream link based on velocity direction
+ link_at_A1 = np.where(temp_Uvel[link_at_B2] >= 0, link_at_A1, link_at_A3)
+ link_at_A3 = np.where(temp_Uvel[link_at_B2] >= 0, link_at_A3, link_at_A1)
+
+ link_at_B1 = self.find_adjacent_links_at_link(
+ link_at_B2, objective_links="vertical"
+ )[
+ :, 0
+ ] # 0: E, Left
+ link_at_B3 = self.find_adjacent_links_at_link(
+ link_at_B2, objective_links="vertical"
+ )[
+ :, 2
+ ] # 2: W, Right
+
+ # Selecting downstream link based on velocity direction
+ link_at_B1 = np.where(temp_Uvel[link_at_B2] >= 0, link_at_B1, link_at_B3)
+ link_at_B3 = np.where(temp_Uvel[link_at_B2] >= 0, link_at_B3, link_at_B1)
+
+ link_at_C1 = self.find_adjacent_links_at_link(
+ link_at_C2, objective_links="vertical"
+ )[
+ :, 0
+ ] # 0: E, Left
+ link_at_C3 = self.find_adjacent_links_at_link(
+ link_at_C2, objective_links="vertical"
+ )[
+ :, 2
+ ] # 2: W, Right
+
+ # Selecting downstream link based on velocity direction
+ link_at_C1 = np.where(temp_Uvel[link_at_B2] >= 0, link_at_C1, link_at_C3)
+ link_at_C3 = np.where(temp_Uvel[link_at_B2] >= 0, link_at_C3, link_at_C1)
+
+ # Getting velocity around the particle
+ vel_at_A1 = np.where(
+ link_at_A1 >= 0, self._vel_at_N[link_at_A1], self._vel_at_N[link_at_A2]
+ )
+ vel_at_A2 = self._vel_at_N[link_at_A2]
+ vel_at_A3 = np.where(
+ link_at_A3 >= 0, self._vel_at_N[link_at_A3], self._vel_at_N[link_at_A2]
+ )
+
+ vel_at_B1 = np.where(
+ link_at_B1 >= 0, self._vel_at_N[link_at_B1], self._vel_at_N[link_at_B2]
+ )
+ vel_at_B2 = self._vel_at_N[link_at_B2]
+ vel_at_B3 = np.where(
+ link_at_B3 >= 0, self._vel_at_N[link_at_B3], self._vel_at_N[link_at_B2]
+ )
+
+ vel_at_C1 = np.where(
+ link_at_C1 >= 0, self._vel_at_N[link_at_C1], self._vel_at_N[link_at_C2]
+ )
+ vel_at_C2 = self._vel_at_N[link_at_C2]
+ vel_at_C3 = np.where(
+ link_at_C3 >= 0, self._vel_at_N[link_at_C3], self._vel_at_N[link_at_C2]
+ )
+
+ # Getting coordinates around the particle
+ x_at_2 = self._grid.xy_of_link[link_at_B2][:, 0]
+ x_at_1 = np.where(
+ temp_Uvel[link_at_B2] >= 0, x_at_2 + self._dx, x_at_2 - self._dx
+ )
+ x_at_3 = np.where(
+ temp_Uvel[link_at_B2] >= 0, x_at_2 - self._dx, x_at_2 + self._dx
+ )
+
+ y_at_B = self._grid.xy_of_link[link_at_B2][:, 1]
+ y_at_A = np.where(
+ self._vel_at_N[link_at_B2] >= 0, y_at_B + self._dy, y_at_B - self._dy
+ )
+ y_at_C = np.where(
+ self._vel_at_N[link_at_B2] >= 0, y_at_B - self._dy, y_at_B + self._dy
+ )
+
+ # Calculating the weights W(i,j) for k around x-direction
+ W1 = (
+ (self._x_at_exit_point - x_at_2)
+ * (self._x_at_exit_point - x_at_3)
+ / ((x_at_1 - x_at_2) * (x_at_1 - x_at_3))
+ )
+ W2 = (
+ (self._x_at_exit_point - x_at_1)
+ * (self._x_at_exit_point - x_at_3)
+ / ((x_at_2 - x_at_1) * (x_at_2 - x_at_3))
+ )
+ W3 = (
+ (self._x_at_exit_point - x_at_1)
+ * (self._x_at_exit_point - x_at_2)
+ / ((x_at_3 - x_at_1) * (x_at_3 - x_at_2))
+ )
+
+ # Interpolation by row around 'y_at_exit_point'
+ A = W1 * vel_at_A1 + W2 * vel_at_A2 + W3 * vel_at_A3
+ B = W1 * vel_at_B1 + W2 * vel_at_B2 + W3 * vel_at_B3
+ C = W1 * vel_at_C1 + W2 * vel_at_C2 + W3 * vel_at_C3
+
+ # Calculating the weghts W(i,j) for l around y-direction
+ W1 = (
+ (self._y_at_exit_point - y_at_B)
+ * (self._y_at_exit_point - y_at_C)
+ / ((y_at_A - y_at_B) * (y_at_A - y_at_C))
+ )
+ W2 = (
+ (self._y_at_exit_point - y_at_A)
+ * (self._y_at_exit_point - y_at_C)
+ / ((y_at_B - y_at_A) * (y_at_B - y_at_C))
+ )
+ W3 = (
+ (self._y_at_exit_point - y_at_A)
+ * (self._y_at_exit_point - y_at_B)
+ / ((y_at_C - y_at_A) * (y_at_C - y_at_B))
+ )
+
+ # Calculating VsL by bicuadratic interpolation
+ self._VsL[tempB2] = W1 * A + W2 * B + W3 * C
+
+ """ Computing viscous terms
+ V-located particles
+ """
+ # Central difference scheme around x- and y- direction for V-located particles
+ self._Vvis = np.zeros_like(self._v_vel)
+
+ tempCalc1 = (
+ self._eddy_viscosity
+ * self._dt
+ * (vel_at_B3 - 2 * vel_at_B2 + vel_at_B1)
+ / (self._dx**2)
+ )
+ tempCalc2 = (
+ self._eddy_viscosity
+ * self._dt
+ * (vel_at_C2 - 2 * vel_at_B2 + vel_at_A2)
+ / (self._dy**2)
+ )
+
+ self._Vvis[tempB2] = tempCalc1 + tempCalc2
+
+ """ Computing advective terms F(U,V)
+ """
+ # Computing advective terms (FU, FV)
+ self._f_vel = np.zeros_like(self._vel_at_N)
+
+ # Adding semi-lagrangian and viscous terms
+ tempCalc1 = self._UsL + self._Uvis
+ tempCalc2 = self._VsL + self._Vvis
+
+ # Including the results according with links directions
+ self._f_vel[self._horizontal_links] = tempCalc1
+ self._f_vel[self._vertical_links] = tempCalc2
+
+ """ Setting G-faces
+ """
+ # Setting G-faces
+ self._g_links = np.zeros_like(self._vel_at_N)
+
+ # Computing G-faces
+ self._g_links = self._h_at_N_at_links * self._f_vel - self._h_at_N_at_links * (
+ 1 - self._theta
+ ) * self._g * self._dt / self._dx * self._grid.calc_diff_at_link(self._eta_at_N)
+
+ # Using only wet-link values, and setting dry links equal to 0 to avoid using wrong values
+ self._g_links = np.where(self._wet_links , self._g_links, 0)
+ self._g_links = np.where(
+ self._grid.status_at_link == 4, 0, self._g_links
+ ) # link is active (0), fixed (2), inactive (4)
+
+ """ Solving semi-implicit scheme with PCG method
+ """
+ # Building the system of equations 'A*x=b'
+ A = np.zeros(
+ (self._number_of_nodes, self._number_of_nodes)
+ ) # Full 'A' matrix with all nodes on it
+ b = np.zeros(self._number_of_nodes) # Full 'b' vector with all nodes on it
+
+ # Getting surrounding locations for core nodes
+ adjacent_nodes = self._grid.adjacent_nodes_at_node[
+ self._core_nodes
+ ] # East, North, West, South
+ adjacent_links = self._grid.links_at_node[
+ self._core_nodes
+ ] # East, North, West, South
+ nodes_location = np.append(
+ adjacent_nodes, np.array([self._core_nodes]).T, axis=1
+ ) # East, North, West, South, Center
+
+ # Boolean to differentiate between core and boundary nodes
+ tempB1 = np.isin(nodes_location, self._core_nodes)
+ # Core node if tempB1 == True, boundary node if tempB1 == False
+ tempB2 = ~tempB1
+ # Boundary node if tempB2 == True, core node if tempB2 == False
+
+ ## Building 'b' for the right-hand side of the system
+ # Calculating 'delta' to build 'b'
+ tempCalc1 = (
+ self._h_at_N_at_links[adjacent_links] * self._vel_at_N[adjacent_links]
+ )
+ tempCalc2 = (
+ self._eta_at_N[self._core_nodes]
+ - (1 - self._theta)
+ * self._dt
+ / self._dx
+ * (tempCalc1[:, 0] - tempCalc1[:, 2])
+ - (1 - self._theta)
+ * self._dt
+ / self._dy
+ * (tempCalc1[:, 1] - tempCalc1[:, 3])
+ )
+
+ tempCalc1 = (
+ self._h_at_N_at_links[adjacent_links]
+ * self._g_links[adjacent_links]
+ / self._a_links[adjacent_links]
+ )
+ b[self._core_nodes] = (
+ tempCalc2
+ - self._theta * self._dt / self._dx * (tempCalc1[:, 0] - tempCalc1[:, 2])
+ - self._theta * self._dt / self._dy * (tempCalc1[:, 1] - tempCalc1[:, 3])
+ )
+
+ ## Building 'A' for the left-hand side of the system
+ # Calculating coefficients for the system of equations
+ tempCalc1 = (
+ self._h_at_N_at_links[adjacent_links] ** 2 / self._a_links[adjacent_links]
+ )
+ coefficients = [
+ -tempCalc1[:, 0] * (self._g * self._theta * self._dt / self._dx) ** 2,
+ -tempCalc1[:, 1] * (self._g * self._theta * self._dt / self._dy) ** 2,
+ -tempCalc1[:, 2] * (self._g * self._theta * self._dt / self._dx) ** 2,
+ -tempCalc1[:, 3] * (self._g * self._theta * self._dt / self._dy) ** 2,
+ 1
+ + (tempCalc1[:, 0] + tempCalc1[:, 2])
+ * (self._g * self._theta * self._dt / self._dx) ** 2
+ + (tempCalc1[:, 1] + tempCalc1[:, 3])
+ * (self._g * self._theta * self._dt / self._dy) ** 2,
+ ]
+ coefficients = np.array(coefficients).T
+
+ ## For loop iterates through every row of 'nodes_location' to:
+ # a) find the node number (row in A) and choose an equation, and
+ # b) find the columns associated with adjacent nodes
+ for row in range(nodes_location.shape[0]):
+ # Getting the current node
+ current_node = nodes_location[row, 4]
+ # Getting the row associated with the current node
+ current_rows_in_A = np.array([current_node])
+ # Getting the columns associated with surrounding nodes
+ current_cols_in_A = nodes_location[row, :]
+
+ # Filling A matrix with coefficients associated with surrounding nodes
+ A[np.ix_(current_rows_in_A, current_cols_in_A)] = (
+ coefficients[row, :] * tempB1[row, :]
+ )
+
+ # Adding known terms (boundary nodes) to the right-hand side of the equation
+ b[current_rows_in_A] = b[current_rows_in_A] - sum(
+ self._eta_at_N[nodes_location[row, :]]
+ * coefficients[row, :]
+ * tempB2[row, :]
+ )
+ # for END
+
+ # Extracting only core nodes to be solved
+ left_hand_side = A[np.ix_(self._core_nodes, self._core_nodes)]
+ right_hand_side = b[self._core_nodes]
+
+ # Applying PCG method to 'LHS*eta=RHS' using np.diag() as a preconditioner for 'LHS'
+ # Preconditioned conjugated gradient output flag:
+ # 0 : successful exit
+ # >0 : convergence to tolerance not achieved, number of iterations
+ # <0 : illegal input or breakdown
+ # Alternative preconditiner: Li = np.linalg.cholesky(left_hand_side)
+ Mi = np.diag(np.diag(left_hand_side))
+ pcg_results = sp.sparse.linalg.cg(
+ left_hand_side,
+ right_hand_side,
+ M=Mi,
+ tol=self._pcg_tolerance,
+ maxiter=self._pcg_max_iterations,
+ atol=0
+ )
+
+ # Getting the new water surface elevation
+ self._eta = np.zeros_like(self._eta_at_N)
+ self._eta[self._core_nodes] = pcg_results[0]
+
+ """ Boundary conditions
+ Radiation Boundary Conditions of Roed & Smedstad (1984) applied on open boundaries
+ Water surface elevation
+ """
+ ## Updating the new WSE ('eta') with the fixed nodes values
+ if self._fixed_nodes_exist is True:
+ self._eta[self._fixed_entry_nodes] = (
+ self._entry_nodes_h_values - self._z[self._fixed_entry_nodes]
+ )
+
+ ## Getting the 1-line-upstream nodes from boundary nodes
+ tempCalc1 = self._grid.active_adjacent_nodes_at_node[self._open_boundary_nodes]
+ open_boundary_nodes_1_backwards = np.extract(tempCalc1 >= 0, tempCalc1)
+
+ ## Getting the 2-line-upstream nodes from boundary nodes
+ # Looking for these nodes
+ tempCalc1 = np.tile(self._open_boundary_nodes, (4, 1)).T
+ tempCalc2 = self._grid.active_adjacent_nodes_at_node[
+ open_boundary_nodes_1_backwards
+ ] # Surrounding nodes to 1-line-upstream
+
+ # Getting the node positions to extract them from all the surrounding nodes
+ tempB1 = np.tile([0, 1, 2, 3], (len(self._open_boundary_nodes), 1))
+ tempB2 = tempB1[tempCalc1 == tempCalc2]
+
+ # It gives me the node indices to extract
+ # folowing the face direction
+ tempB1 = np.where(tempB2 == 0, 2, tempB2)
+ tempB1 = np.where(tempB2 == 1, 3, tempB1)
+ tempB1 = np.where(tempB2 == 2, 0, tempB1)
+ tempB1 = np.where(tempB2 == 3, 1, tempB1)
+
+ open_boundary_nodes_2_backwards = tempCalc2[
+ [range(tempCalc2.shape[0])], tempB1
+ ][0]
+
+ ## Getting WSE at different locations
+ eta_at_N_at_B = self._eta_at_N[self._open_boundary_nodes]
+ eta_at_N_at_B_1 = self._eta_at_N[open_boundary_nodes_1_backwards]
+ eta_at_N_1_at_B_1 = self._eta_at_N_1[open_boundary_nodes_1_backwards]
+ eta_at_N_1_at_B_2 = self._eta_at_N_1[open_boundary_nodes_2_backwards]
+
+ ## Computing boundary condition
+ tempCalc1 = eta_at_N_at_B_1 - eta_at_N_1_at_B_1
+ tempCalc2 = eta_at_N_1_at_B_1 - eta_at_N_1_at_B_2
+ tempCalc3 = np.where(tempCalc2 == 0, 1, tempCalc2)
+
+ Ce = np.where(tempCalc2 == 0, 0, tempCalc1 / tempCalc3 * (-self._dx / self._dt))
+ Ce = np.where(tempCalc2 == 0, 0, Ce)
+
+ # eta[open_boundary_nodes] = tempCalc1/tempCalc2
+ self._eta[self._open_boundary_nodes] = np.where(
+ Ce >= 0, eta_at_N_at_B_1, eta_at_N_at_B
+ )
+ self._eta = np.where(
+ abs(self._eta) > abs(self._z), -self._z, self._eta
+ ) # Correcting WSE below topographic elevation
+
+ self._eta_at_links = self._grid.map_mean_of_link_nodes_to_link(self._eta)
+
+ # Corner nodes treatment
+ self._eta[self._corner_nodes] = np.mean(
+ self._eta[self._adjacent_nodes_at_corner_nodes], axis=1
+ )
+
+ """ Updating water velocity
+ """
+ # tempB1 : Considering only wet cells
+ # tempB2 : Cells with elevation below the water surface elevation
+ tempB1 = np.where(
+ abs(self._eta[self._grid.node_at_link_head])
+ <= abs(self._z[self._grid.node_at_link_head] - self._threshold_depth),
+ 1,
+ 0,
+ )
+ tempB2 = np.where(
+ abs(self._z[self._grid.node_at_link_head])
+ > abs(self._eta[self._grid.node_at_link_tail]),
+ 1,
+ 0,
+ )
+ tempB1 = tempB1 + tempB2
+ tempB1 = np.where(tempB1 > 1, 1, 0)
+
+ # Updating water velocity
+ tempCalc1 = (
+ self._theta
+ * self._g
+ * self._dt
+ / self._dx
+ * (self._grid.calc_diff_at_link(self._eta))
+ * self._h_at_N_at_links
+ / self._a_links
+ )
+ self._vel = self._g_links / self._a_links - tempCalc1 * tempB1
+
+ # Only updating velocity on wet cells
+ self._vel = np.where(self._wet_links , self._vel, 0)
+
+ """ Boundary conditions
+ Radiation Boundary Conditions of Roed & Smedstad (1984) applied on open boundaries
+ Water velocity
+ """
+ ## Updating the new Velocity with the fixed links values
+ if self._fixed_links_exist is True:
+ self._vel[self._fixed_entry_links] = self._entry_links_vel_values
+
+ ## Getting the boundary links
+ tempB1 = [i in self._open_boundary_links for i in self._active_links]
+ open_boundary_active_links = self._grid.active_links[tempB1]
+
+ ## Getting the 1-line-upstream links from boundary links
+ tempCalc1 = np.tile(open_boundary_active_links, (4, 1)).T
+ tempCalc2 = self._grid.links_at_node[open_boundary_nodes_1_backwards]
+
+ # Getting the link positions to extract them from all the surrounding nodes
+ tempB1 = np.tile([0, 1, 2, 3], (len(self._open_boundary_nodes), 1))
+ tempB2 = tempB1[tempCalc1 == tempCalc2]
+
+ # It gives me the link indices to extract
+ # folowing the face direction
+ tempB1 = np.where(tempB2 == 0, 2, tempB2)
+ tempB1 = np.where(tempB2 == 1, 3, tempB1)
+ # tempB1 is where the target link is located
+ tempB1 = np.where(tempB2 == 2, 0, tempB1)
+ # tempB2 is where the upstream link is located
+ tempB1 = np.where(tempB2 == 3, 1, tempB1)
+
+ open_boundary_active_links_1_backwards = tempCalc2[
+ [range(tempCalc2.shape[0])], tempB1
+ ][0]
+
+ ### Getting the 2-line-upstream links from boundary nodes
+ tempCalc1 = np.tile(open_boundary_active_links_1_backwards, (4, 1)).T
+ tempCalc2 = self._grid.links_at_node[open_boundary_nodes_2_backwards]
+
+ # Getting the link positions to extract them from all the surrounding nodes
+ tempB1 = np.tile([0, 1, 2, 3], (len(self._open_boundary_nodes), 1))
+ tempB2 = tempB1[(tempCalc1 == tempCalc2)]
+
+ # It gives me the link indices to extract
+ # folowing the face direction
+ tempB1 = np.where(tempB2 == 0, 2, tempB2)
+ tempB1 = np.where(tempB2 == 1, 3, tempB1)
+ # tempB1 is where the target link is located
+ tempB1 = np.where(tempB2 == 2, 0, tempB1)
+ # tempB2 is where the upstream link is located
+ tempB1 = np.where(tempB2 == 3, 1, tempB1)
+
+ open_boundary_active_links_2_backwards = tempCalc2[
+ [range(tempCalc2.shape[0])], tempB1
+ ][0]
+
+ ## Getting water velocity at different locations
+ vel_at_N_at_B = self._vel_at_N[open_boundary_active_links]
+ vel_at_N_at_B_1 = self._vel_at_N[open_boundary_active_links_1_backwards]
+ vel_at_N_1_at_B_1 = self._vel_at_N_1[open_boundary_active_links_1_backwards]
+ vel_at_N_1_at_B_2 = self._vel_at_N_1[open_boundary_active_links_2_backwards]
+
+ ## Computing boundary condition
+ tempCalc1 = vel_at_N_at_B_1 - vel_at_N_1_at_B_1
+ tempCalc2 = vel_at_N_1_at_B_1 - vel_at_N_1_at_B_2
+ tempCalc3 = np.where(tempCalc2 == 0, 1, tempCalc2)
+
+ Ce = np.where(tempCalc2 == 0, 0, tempCalc1 / tempCalc3 * (-self._dx / self._dt))
+ Ce = np.where(tempCalc2 == 0, 0, Ce)
+
+ self._vel[open_boundary_active_links] = np.where(
+ Ce >= 0, vel_at_N_at_B_1, vel_at_N_at_B
+ )
+
+ """ Updating water depth at links
+ """
+ # Using only values where the WSE is above the topographic elevation
+ tempB1 = np.where(abs(self._eta) <= abs(self._z - self._threshold_depth), 1, 0)
+
+ # Updating water depth at links
+ tempCalc1 = (
+ self._z_at_links + self._eta[self._grid.node_at_link_head]
+ ) * tempB1[self._grid.node_at_link_head]
+ tempCalc2 = (
+ self._z_at_links + self._eta[self._grid.node_at_link_tail]
+ ) * tempB1[self._grid.node_at_link_tail]
+ tempCalc3 = np.zeros_like(self._h_at_N_at_links)
+
+ self._h_at_links = np.array((tempCalc1, tempCalc2, tempCalc3)).max(axis=0)
+
+ # Applying boundary condition at links
+ self._h_at_links[self._open_boundary_links] = (
+ self._z_at_links[self._open_boundary_links]
+ + self._eta_at_links[self._open_boundary_links]
+ )
+
+ # Wet cells threshold
+ self._h_at_links = np.where(
+ self._h_at_links < self._threshold_depth, 0, self._h_at_links
+ )
+
+ # Updating wet links
+ self._wet_links = np.where(
+ self._h_at_links >= self._threshold_depth, True, False
+ )
+ self._vel = self._vel * self._wet_links
+
+ """ Calculating average water depth at nodes
+ """
+ # If a node is dry, using only surrounding links such that 'WSE' is above 'z'
+ # If a node is wet, using all surrounding links even if 'WSE' is below 'z' (jumps)
+
+ # Checking surrounding wet links
+ surrounding_links = self._grid.links_at_node[self._core_nodes]
+
+ # Checking whether the core node is wet (T) or dry (F)
+ tempB1 = abs(self._eta[self._core_nodes]) < abs(
+ self._z[self._core_nodes] - self._threshold_depth
+ )
+
+ # Checking whether surrounding links are wet (T) or dry (F)
+ tempB2 = self._wet_links[surrounding_links]
+
+ # Checking whether surrounding 'WSE' links are above (T) or below (F) 'z' at nodes
+ tempB3 = (
+ abs(self._eta_at_links[surrounding_links])
+ < abs(self._z[self._core_nodes] - self._threshold_depth)[:, None]
+ )
+
+ # Getting the number of wet links around each core node, satisfying tempB2, and avoiding divisions by zero
+ tempCalc2 = np.sum(tempB2 * 1, axis=1)
+ tempCalc2 = np.where(tempCalc2 == 0, -9999, tempCalc2)
+
+ # Getting the number of wet links around each core node, satisfying tempB3, and avoiding divisions by zero
+ tempCalc3 = np.sum(tempB2 * tempB3 * 1, axis=1)
+ tempCalc3 = np.where(tempCalc3 == 0, -9999, tempCalc3)
+
+ # Updating water depth
+ # h = h_at_N - rmg.calc_net_flux_at_node(h_at_links*vel) # (influx if negative)
+ self._h[self._core_nodes] = np.where(
+ tempCalc3 > 0,
+ np.sum(self._h_at_links[surrounding_links] * tempB2 * tempB3, axis=1)
+ / tempCalc3,
+ 0,
+ ) # Dry nodes, tempB1 == False
+
+ ### Updating boundary nodes
+ if self._fixed_nodes_exist is True:
+ self._h[self._fixed_entry_nodes] = self._entry_nodes_h_values
+ self._h[self._open_boundary_nodes] = (
+ self._eta[self._open_boundary_nodes] + self._z[self._open_boundary_nodes]
+ )
+ self._h = np.where(self._h < self._threshold_depth, 0, self._h)
+
+ # Corner nodes treatment
+ self._h[self._corner_nodes] = np.mean(
+ self._h[self._adjacent_nodes_at_corner_nodes], axis=1
+ )
+
+ # Updating wet nodes
+ self._wet_nodes = np.where(self._h >= self._threshold_depth, True, False)
+
+ """ Storing values in the grid
+ """
+ self._grid.at_node["surface_water__depth"] = self._h
+ self._grid.at_link["surface_water__velocity"] = self._vel
+ self._grid.at_node["surface_water__elevation"] = self._eta + (
+ self._max_elevation + self._additional_z
+ )
+
+ self._grid.at_link["surface_water__velocity_at_N-1"] = self._vel_at_N
+ self._grid.at_node["surface_water__elevation_at_N-1"] = self._eta_at_N + (
+ self._max_elevation + self._additional_z
+ )
+ self._grid.at_node["surface_water__elevation_at_N-2"] = self._eta_at_N_1 + (
+ self._max_elevation + self._additional_z
+ )
+
+ """ Storing values at previous time steps
+ """
+ self._eta_at_N = self._eta.copy()
+ self._eta_at_N_1 = self._eta_at_N.copy()
+ self._eta_at_N_2 = self._eta_at_N_1.copy()
+
+ self._h_at_N = self._h.copy()
+ self._h_at_N_at_links = self._h_at_links.copy()
+
+ self._vel_at_N = self._vel.copy()
+ self._vel_at_N_1 = self._vel_at_N.copy()
diff --git a/notebooks/tutorials/river_flow_dynamics/river_flow_dynamics_tutorial.ipynb b/notebooks/tutorials/river_flow_dynamics/river_flow_dynamics_tutorial.ipynb
new file mode 100644
index 0000000000..e0e4dccd8f
--- /dev/null
+++ b/notebooks/tutorials/river_flow_dynamics/river_flow_dynamics_tutorial.ipynb
@@ -0,0 +1,997 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ ""
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# 2D Surface Water Flow component\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "
\n",
+ "For more Landlab tutorials, click here: https://landlab.readthedocs.io/en/latest/user_guide/tutorials.html\n",
+ ""
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Overview\n",
+ "\n",
+ "This notebook demonstrate the usage of the `river flow dynamics` Landlab component. The component runs a semi-implicit, semi-Lagrangian finite-volume approximation to the depth-averaged 2D shallow-water equations of Casulli and Cheng (1992) and related work.\n",
+ "\n",
+ "### Theory\n",
+ "\n",
+ "The depth-averaged 2D shallow-water equations are the simplification of the Navier-Stokes equations, which correspond to the balance of momentum and mass in the fluid. It is possible to simplify these equations by assuming a well-mixed water column and a small water depth to width ratio, where a vertical integration results in depth-averaged equations. These require boundary conditions at the top and bottom of the water column, which are provided by the wind stress and the Manning-Chezy formula, respectively:\n",
+ "\n",
+ "$$\n",
+ "\\frac{\\partial U}{\\partial t}\n",
+ "+ U\\frac{\\partial U}{\\partial x} + V\\frac{\\partial U}{\\partial y}\n",
+ "= \n",
+ "- g\\frac{\\partial \\eta}{\\partial x}\n",
+ "+ \\epsilon\\left(\\frac{\\partial^2 U}{\\partial x^2} + \\frac{\\partial^2 U}{\\partial y^2}\\right)\n",
+ "+ \\frac{\\gamma_T(U_a - U)}{H} - g\\frac{\\sqrt{U^2 + V^2}}{Cz^2}U + \\mathbf{f}V\n",
+ "$$\n",
+ "\n",
+ "$$\n",
+ "\\frac{\\partial V}{\\partial t}\n",
+ "+ U\\frac{\\partial V}{\\partial x} + V\\frac{\\partial V}{\\partial y}\n",
+ "= \n",
+ "- g\\frac{\\partial \\eta}{\\partial y}\n",
+ "+ \\epsilon\\left(\\frac{\\partial^2 V}{\\partial x^2} + \\frac{\\partial^2 V}{\\partial y^2}\\right)\n",
+ "+ \\frac{\\gamma_T(V_a - V)}{H} - g\\frac{\\sqrt{U^2 + V^2}}{Cz^2}V + \\mathbf{f}U\n",
+ "$$\n",
+ "\n",
+ "$$\n",
+ "\\frac{\\partial \\eta}{\\partial t}\n",
+ "+ \\frac{\\partial (HU)}{\\partial x} + \\frac{\\partial (HV)}{\\partial y}\n",
+ "= 0\n",
+ "$$\n",
+ "\n",
+ "where $U$ is the water velocity in the $x$-direction, $V$ is the water velocity in the $y$-direction, $H$ is the water depth, $\\eta$ is the water surface elevation, $Cz$ is the Chezy friction coefficient, and $t$ is time. For the constants $g$ is the gravity acceleration, $\\epsilon$ is the horizontal eddy viscosity, $\\mathbf{f}$ is the Coriolis parameter, $\\gamma_T$ is the wind stress coefficient, and $U_a$ and $V_a$ are the prescribed wind velocities.\n",
+ "\n",
+ "### Numerical representation\n",
+ "\n",
+ "A semi-implicit, semi-Lagrangian, finite volume numerical approximation represents the depth averaged, 2D shallow-water equations described before. The water surface elevation, $\\eta$, is defined at the center of each computational volume (nodes). Water depth, $H$, and velocity components, $U$ and $V$, are defined at the midpoint of volume faces (links). The finite volume structure provides a control volume representation that is inherently mass conservative.\n",
+ "\n",
+ "The combination of a semi-implciit water surface elevation solution and a semi-Lagrangian representation of advection provides the advantages of a stable solution and of time steps that exceed the CFL criterion. In the semi-implicit process, $\\eta$ in the momentum equations, and the velocity divergence in the continuity equation, are trated implicitly. The advective terms in the momentum equations, are discretized explicitly. See the cited literature for more details.\n",
+ "\n",
+ "### The component\n",
+ "\n",
+ "Import the needed libraries:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import numpy as np\n",
+ "import matplotlib.pyplot as plt\n",
+ "from landlab import RasterModelGrid\n",
+ "from landlab.components import river_flow_dynamics\n",
+ "from landlab.io import read_esri_ascii\n",
+ "from landlab.plot.imshow import imshow_grid"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Information about the component\n",
+ "\n",
+ "Using the class name as argument for the `help` function returns descriptions of the various methods and parameters."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "metadata": {
+ "scrolled": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Help on class river_flow_dynamics in module landlab.components.river_flow_dynamics.river_flow_dynamics:\n",
+ "\n",
+ "class river_flow_dynamics(landlab.core.model_component.Component)\n",
+ " | river_flow_dynamics(grid, dt=0.01, eddy_viscosity=0.0001, mannings_n=0.012, threshold_depth=0.01, theta=0.5, fixed_entry_nodes=[], fixed_entry_links=[], entry_nodes_h_values=[], entry_links_vel_values=[], pcg_tolerance=1e-05, pcg_max_iterations=None, surface_water__elevation_at_N_1=None, surface_water__elevation_at_N_2=None, surface_water__velocity_at_N_1=None)\n",
+ " | \n",
+ " | Simulate surface fluid flow based on Casulli and Cheng (1992).\n",
+ " | \n",
+ " | Landlab component that simulates surface fluid flow using the Casulli and Cheng (1992)\n",
+ " | approximations of the 2D shallow water equations.\n",
+ " | \n",
+ " | This components calculates water depth and velocity across the raster grid, given\n",
+ " | a certain input discharge.\n",
+ " | \n",
+ " | References\n",
+ " | ----------\n",
+ " | **Required Software Citation(s) Specific to this Component**\n",
+ " | \n",
+ " | None Listed\n",
+ " | \n",
+ " | **Additional References**\n",
+ " | \n",
+ " | Casulli, V., Cheng, R.T. (1992). “Semi-implicit finite difference methods for\n",
+ " | three-dimensional shallow water flow”. International Journal for Numerical Methods\n",
+ " | in Fluids. 15: 629-648.\n",
+ " | https://doi.org/10.1002/fld.1650150602\n",
+ " | \n",
+ " | Method resolution order:\n",
+ " | river_flow_dynamics\n",
+ " | landlab.core.model_component.Component\n",
+ " | builtins.object\n",
+ " | \n",
+ " | Methods defined here:\n",
+ " | \n",
+ " | __init__(self, grid, dt=0.01, eddy_viscosity=0.0001, mannings_n=0.012, threshold_depth=0.01, theta=0.5, fixed_entry_nodes=[], fixed_entry_links=[], entry_nodes_h_values=[], entry_links_vel_values=[], pcg_tolerance=1e-05, pcg_max_iterations=None, surface_water__elevation_at_N_1=None, surface_water__elevation_at_N_2=None, surface_water__velocity_at_N_1=None)\n",
+ " | Simulate the vertical-averaged surface fluid flow\n",
+ " | \n",
+ " | Landlab component that simulates surface fluid flow using the Casulli and Cheng (1992)\n",
+ " | approximations of the 2D shallow water equations.\n",
+ " | \n",
+ " | This components calculates water depth and velocity across the raster grid, given\n",
+ " | a certain input discharge.\n",
+ " | \n",
+ " | Parameters\n",
+ " | ----------\n",
+ " | grid : RasterModelGrid\n",
+ " | A grid.\n",
+ " | dt : float, optional\n",
+ " | Time step in seconds. If not given, it is calculated from CFL condition.\n",
+ " | eddy_viscosity : float, optional\n",
+ " | Eddy viscosity coefficient. Default = 1e-4 :math:`m^2 / s`\n",
+ " | mannings_n : float or array_like, optional\n",
+ " | Manning's roughness coefficient. Default = 0.012 :math:`s / m^1/3`\n",
+ " | threshold_depth : float, optional\n",
+ " | Threshold at which a cell is considered wet. Default = 0.01 m\n",
+ " | theta : float, optional\n",
+ " | Degree of 'implicitness' of the solution, ranging between 0.5 and 1.0. Default = 0.5.\n",
+ " | When it is equal to 0.5, the approximation is centered in time.\n",
+ " | When it is equal to 1.0, the approximation is fully implicit.\n",
+ " | fixed_entry_nodes : array_like or None, optional\n",
+ " | Node IDs where flow enters the domain (Dirichlet boundary condition).\n",
+ " | If not provided, the water already present in the domain is not renewed.\n",
+ " | fixed_entry_links : array_like or None, optional\n",
+ " | Link IDs where flow enters the domain (Dirichlet boundary condition).\n",
+ " | If not provided, the water already present in the domain is not renewed.\n",
+ " | entry_nodes_h_values : array_like, optional\n",
+ " | Water depth values at nodes where flow enters the domain (Dirichlet boundary condition).\n",
+ " | If not provided, the water already present in the domain is not renewed.\n",
+ " | entry_links_vel_values : array_like, optional\n",
+ " | Water velocity values at links where flow enters the domain (Dirichlet boundary condition).\n",
+ " | If not provided, the water already present in the domain is not renewed.\n",
+ " | pcg_tolerance : float, optional\n",
+ " | Preconditioned Conjugate Gradient method argument. Tolerance for convergence.\n",
+ " | Default = 1e-05\n",
+ " | pcg_max_iterations : integer, optional\n",
+ " | Preconditioned Conjugate Gradient method argument. Maximum number of iterations.\n",
+ " | Iteration will stop after maxiter steps even if the specified tolerance has not been achieved.\n",
+ " | Default = None\n",
+ " | surface_water__elevation_at_N_1: float, optional\n",
+ " | Water surface elevation at time N-1.\n",
+ " | units: m, mapping: node\n",
+ " | surface_water__elevation_at_N_2: float, optional\n",
+ " | Water surface elevation at time N-2.\n",
+ " | units: m, mapping: node\n",
+ " | surface_water__velocity_at_N_1: float, optional\n",
+ " | Speed of water flow above the surface at time N-1\",\n",
+ " | units: m/2, mapping: link\n",
+ " | \n",
+ " | find_adjacent_links_at_link(self, current_link, objective_links='horizontal')\n",
+ " | Get adjacent links to the link.\n",
+ " | \n",
+ " | This function finds the links at right, above, left and below the given link.\n",
+ " | Similar purpose to the \"adjacent_nodes_at_node\" function.\n",
+ " | Return the adjacent links in as many rows as given links.\n",
+ " | Link IDs are returned as columns in clock-wise order starting from East (E, N, W, S).\n",
+ " | \n",
+ " | find_nearest_link(self, x_coordinates, y_coordinates, objective_links='all')\n",
+ " | Link nearest a point.\n",
+ " | \n",
+ " | Find the index to the link nearest the given x, y coordinates.\n",
+ " | Returns the indices of the links nearest the given coordinates.\n",
+ " | \n",
+ " | path_line_tracing(self)\n",
+ " | \" Path line tracing algorithm.\n",
+ " | \n",
+ " | This function implements the semi-analytical path line tracing method of Pollock (1988).\n",
+ " | \n",
+ " | The semi-analytical path line tracing method was developed for particle tracking in ground\n",
+ " | water flow models. The assumption that each directional velocity component varies linearly\n",
+ " | in its coordinate directions within each computational volume or cell underlies the method.\n",
+ " | Linear variation allows the derivation of an analytical expression for the path line of a\n",
+ " | particle across a volume.\n",
+ " | \n",
+ " | Given an initial point located at each volume faces of the domain, particle trayectories are\n",
+ " | traced backwards on time. Then, this function returns the departure point of the particle at\n",
+ " | the beginning of the time step.\n",
+ " | \n",
+ " | run_one_step(self)\n",
+ " | Calculate water depth and water velocity for a time period dt.\n",
+ " | \n",
+ " | ----------------------------------------------------------------------\n",
+ " | Methods inherited from landlab.core.model_component.Component:\n",
+ " | \n",
+ " | initialize_optional_output_fields(self)\n",
+ " | Create fields for a component based on its optional field outputs,\n",
+ " | if declared in _optional_var_names.\n",
+ " | \n",
+ " | This method will create new fields (without overwrite) for any\n",
+ " | fields output by the component as optional. New fields are\n",
+ " | initialized to zero. New fields are created as arrays of floats,\n",
+ " | unless the component also contains the specifying property\n",
+ " | _var_type.\n",
+ " | \n",
+ " | initialize_output_fields(self, values_per_element=None)\n",
+ " | Create fields for a component based on its input and output var\n",
+ " | names.\n",
+ " | \n",
+ " | This method will create new fields (without overwrite) for any fields\n",
+ " | output by, but not supplied to, the component. New fields are\n",
+ " | initialized to zero. Ignores optional fields. New fields are created as\n",
+ " | arrays of floats, unless the component specifies the variable type.\n",
+ " | \n",
+ " | Parameters\n",
+ " | ----------\n",
+ " | values_per_element: int (optional)\n",
+ " | On occasion, it is necessary to create a field that is of size\n",
+ " | (n_grid_elements, values_per_element) instead of the default size\n",
+ " | (n_grid_elements,). Use this keyword argument to acomplish this\n",
+ " | task.\n",
+ " | \n",
+ " | ----------------------------------------------------------------------\n",
+ " | Class methods inherited from landlab.core.model_component.Component:\n",
+ " | \n",
+ " | from_path(grid, path) from builtins.type\n",
+ " | Create a component from an input file.\n",
+ " | \n",
+ " | Parameters\n",
+ " | ----------\n",
+ " | grid : ModelGrid\n",
+ " | A landlab grid.\n",
+ " | path : str or file_like\n",
+ " | Path to a parameter file, contents of a parameter file, or\n",
+ " | a file-like object.\n",
+ " | \n",
+ " | Returns\n",
+ " | -------\n",
+ " | Component\n",
+ " | A newly-created component.\n",
+ " | \n",
+ " | var_definition(name) from builtins.type\n",
+ " | Get a description of a particular field.\n",
+ " | \n",
+ " | Parameters\n",
+ " | ----------\n",
+ " | name : str\n",
+ " | A field name.\n",
+ " | \n",
+ " | Returns\n",
+ " | -------\n",
+ " | tuple of (*name*, *description*)\n",
+ " | A description of each field.\n",
+ " | \n",
+ " | var_help(name) from builtins.type\n",
+ " | Print a help message for a particular field.\n",
+ " | \n",
+ " | Parameters\n",
+ " | ----------\n",
+ " | name : str\n",
+ " | A field name.\n",
+ " | \n",
+ " | var_loc(name) from builtins.type\n",
+ " | Location where a particular variable is defined.\n",
+ " | \n",
+ " | Parameters\n",
+ " | ----------\n",
+ " | name : str\n",
+ " | A field name.\n",
+ " | \n",
+ " | Returns\n",
+ " | -------\n",
+ " | str\n",
+ " | The location ('node', 'link', etc.) where a variable is defined.\n",
+ " | \n",
+ " | var_type(name) from builtins.type\n",
+ " | Returns the dtype of a field (float, int, bool, str...).\n",
+ " | \n",
+ " | Parameters\n",
+ " | ----------\n",
+ " | name : str\n",
+ " | A field name.\n",
+ " | \n",
+ " | Returns\n",
+ " | -------\n",
+ " | dtype\n",
+ " | The dtype of the field.\n",
+ " | \n",
+ " | var_units(name) from builtins.type\n",
+ " | Get the units of a particular field.\n",
+ " | \n",
+ " | Parameters\n",
+ " | ----------\n",
+ " | name : str\n",
+ " | A field name.\n",
+ " | \n",
+ " | Returns\n",
+ " | -------\n",
+ " | str\n",
+ " | Units for the given field.\n",
+ " | \n",
+ " | ----------------------------------------------------------------------\n",
+ " | Static methods inherited from landlab.core.model_component.Component:\n",
+ " | \n",
+ " | __new__(cls, *args, **kwds)\n",
+ " | Create and return a new object. See help(type) for accurate signature.\n",
+ " | \n",
+ " | ----------------------------------------------------------------------\n",
+ " | Readonly properties inherited from landlab.core.model_component.Component:\n",
+ " | \n",
+ " | cite_as\n",
+ " | Citation information for component.\n",
+ " | \n",
+ " | Return required software citation, if any. An empty string indicates\n",
+ " | that no citations other than the standard Landlab package citations are\n",
+ " | needed for the component.\n",
+ " | \n",
+ " | Citations are provided in BibTeX format.\n",
+ " | \n",
+ " | Returns\n",
+ " | -------\n",
+ " | cite_as\n",
+ " | \n",
+ " | coords\n",
+ " | Return the coordinates of nodes on grid attached to the\n",
+ " | component.\n",
+ " | \n",
+ " | definitions\n",
+ " | Get a description of each field.\n",
+ " | \n",
+ " | Returns\n",
+ " | -------\n",
+ " | tuple of (*name*, *description*)\n",
+ " | A description of each field.\n",
+ " | \n",
+ " | grid\n",
+ " | Return the grid attached to the component.\n",
+ " | \n",
+ " | input_var_names\n",
+ " | Names of fields that are used by the component.\n",
+ " | \n",
+ " | Returns\n",
+ " | -------\n",
+ " | tuple of str\n",
+ " | Tuple of field names.\n",
+ " | \n",
+ " | name\n",
+ " | Name of the component.\n",
+ " | \n",
+ " | Returns\n",
+ " | -------\n",
+ " | str\n",
+ " | Component name.\n",
+ " | \n",
+ " | optional_var_names\n",
+ " | Names of fields that are optionally provided by the component, if\n",
+ " | any.\n",
+ " | \n",
+ " | Returns\n",
+ " | -------\n",
+ " | tuple of str\n",
+ " | Tuple of field names.\n",
+ " | \n",
+ " | output_var_names\n",
+ " | Names of fields that are provided by the component.\n",
+ " | \n",
+ " | Returns\n",
+ " | -------\n",
+ " | tuple of str\n",
+ " | Tuple of field names.\n",
+ " | \n",
+ " | shape\n",
+ " | Return the grid shape attached to the component, if defined.\n",
+ " | \n",
+ " | unit_agnostic\n",
+ " | Whether the component is unit agnostic.\n",
+ " | \n",
+ " | If True, then the component is unit agnostic. Under this condition a\n",
+ " | user must still provide consistent units across all input arguments,\n",
+ " | keyword arguments, and fields. However, when ``unit_agnostic`` is True\n",
+ " | the units specified can be interpreted as dimensions.\n",
+ " | \n",
+ " | When False, then the component requires inputs in the specified units.\n",
+ " | \n",
+ " | Returns\n",
+ " | -------\n",
+ " | bool\n",
+ " | \n",
+ " | units\n",
+ " | Get the units for all field values.\n",
+ " | \n",
+ " | Returns\n",
+ " | -------\n",
+ " | tuple or str\n",
+ " | Units for each field.\n",
+ " | \n",
+ " | var_mapping\n",
+ " | Location where variables are defined.\n",
+ " | \n",
+ " | Returns\n",
+ " | -------\n",
+ " | tuple of (name, location)\n",
+ " | Tuple of variable name and location ('node', 'link', etc.) pairs.\n",
+ " | \n",
+ " | ----------------------------------------------------------------------\n",
+ " | Data descriptors inherited from landlab.core.model_component.Component:\n",
+ " | \n",
+ " | __dict__\n",
+ " | dictionary for instance variables (if defined)\n",
+ " | \n",
+ " | __weakref__\n",
+ " | list of weak references to the object (if defined)\n",
+ " | \n",
+ " | current_time\n",
+ " | Current time.\n",
+ " | \n",
+ " | Some components may keep track of the current time. In this case, the\n",
+ " | ``current_time`` attribute is incremented. Otherwise it is set to None.\n",
+ " | \n",
+ " | Returns\n",
+ " | -------\n",
+ " | current_time\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "help(river_flow_dynamics)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Examples\n",
+ "\n",
+ "-- --\n",
+ "\n",
+ "### Example 1: Flow in a rectangular channel 6.0 m long\n",
+ "\n",
+ "This first basic example illustrates water flowing through a rectangular channel 1.0 $m$ wide and 6.0 $m$ long. Our channel is made in concrete, so we choose a Manning's roughness coefficient equal to 0.012 $s/m^\\frac{1}{3}$, and it has a slope of 0.01 $m/m$.\n",
+ "\n",
+ "We specify some basic parameters such as the grid resolution, time step duration, number of time steps, and the domain dimensions by specifying the number of columns and rows. "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Basic parameters\n",
+ "n = 0.012 # Manning's roughness coefficient, [s/m^(1/3)]\n",
+ "S0 = 0.01 # Channel slope [m/m]\n",
+ "\n",
+ "# Simulation parameters\n",
+ "n_timesteps = 100 # Number of timesteps\n",
+ "dt = 0.1 # Timestep duration, [s]\n",
+ "nrows = 20 # Number of node rows\n",
+ "ncols = 60 # Number of node cols\n",
+ "dx = 0.1 # Node spacing in the x-direction, [m]\n",
+ "dy = 0.1 # Node spacing in the y-direction, [m]"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Create the grid:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Create and set up the grid\n",
+ "grid = RasterModelGrid((nrows, ncols), xy_spacing=(dx, dy))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Create the elevation field and define the topography to represent our rectangular channel:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 6,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# The grid represents a basic rectangular channel with slope equal to 0.01 m/m\n",
+ "te = grid.add_zeros(\"topographic__elevation\", at=\"node\")\n",
+ "te += 1.0 - S0*grid.x_of_node\n",
+ "te[grid.y_of_node > 1.5] = 2.5\n",
+ "te[grid.y_of_node < 0.5] = 2.5"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "We show a top view of the domain:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh8AAAGFCAYAAABdSJFpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAqlklEQVR4nO3df3BV9Z3/8dclMQlq7nWD5JckGlulLCjyTdwSlF8CcYOTkQ470toxVMGvGUNQU5w1uF9EyzbbLrrRDQRo0ZSKyGj45TSLZNYmEYUpwUSotSzU1ARMRJjtvSTVhCTn+wfNHSO/7s29OefcnOej8/njHs7nnHcu1bx9vz/nc1yGYRgCAAAwyQirAwAAAM5C8gEAAExF8gEAAExF8gEAAExF8gEAAExF8gEAAExF8gEAAEwVbXUAAAA43VdffaXu7u6QrxMTE6O4uLgwRDS0SD4AALDQV199pYyMDLW3t4d8reTkZDU3N9s+ASH5AADAQt3d3Wpvb1dra6vcbvegr+Pz+ZSWlqbu7m6SDwAAcHlud7zc7vgQrhA5b0sh+QAAwAYMw1Aor1uLpFe18bQLAAAwFZUPAABswVBorZPIqXyQfAAAYAvOST5ouwAAAFNR+QAAwAactOCU5AMAAFug7QIAADAkqHwAAGALzql8kHwAAGADrPkAAAAmc07lgzUfAADAVFQ+AACwBedUPkg+AACwASet+aDtAgAATEXlAwAAW6DtAgAATOWc5IO2CwAADlRaWqrbb79d8fHxSkxM1Lx583TkyJGA57/33nuKjo7WbbfdFvS9ST4AALCB/gWnoYxg1NXVqbCwUPv371dNTY16enqUk5Ojzs7Oy871er3Kz8/XrFmzBvWzuoxIWh4LAMAw4/P55PF4dOpUi9xud0jXufbadHm93kFd54svvlBiYqLq6uo0bdq0S577/e9/XzfddJOioqK0Y8cONTU1BXUvKh8AAAwjPp9vwOjq6gpontfrlSQlJCRc8rxXXnlFf/rTn/TMM88MOkaSDwAAbMEIw5DS0tLk8Xj8o7S09PJ3NgwVFxfrzjvv1IQJEy563tGjR/XUU09p8+bNio4e/DMrPO0CAIANhGuTsdbW1gFtl9jY2MvOXbJkiQ4dOqS9e/de9Jze3l7df//9evbZZ3XzzTcPOk6JNR8AAFiqf83HF180h7zmY/TojKDXfBQVFWnHjh2qr69XRkbGRc/7y1/+or/7u79TVFSU/1hfX58Mw1BUVJT27Nmju+66K6B7UvkAAMCBDMNQUVGRtm/frtra2ksmHpLkdrt1+PDhAcfWrl2rd955R2+++eZl538dyQcAALZgSOoLcX7gCgsL9dprr2nnzp2Kj49Xe3u7JMnj8WjkyJGSpJKSEp04cUKbNm3SiBEjzlsPkpiYqLi4uEuuE7kQFpwCAGADZu/zUVFRIa/XqxkzZiglJcU/tm7d6j+nra1NLS0t4f5RWfMBAICV+td8nDz5J7nd8SFc54wSE7816H0+zETbBQAAW3DOu11IPgAAsAXnJB+s+QAAAKai8gEAgA2Ea5OxSEDyAQCALTin7ULyAQCALTgn+WDNBwAAMBWVDwAAbIA1HwAAwGS0XQAAAIYElQ8AAGyAtgsAADBZn0J7q20oc81F2wUAAJiKygcAALbgnAWnJB8AANiAk9Z80HYBAACmovIBAIAt0HYBAACmIvkAAAAmYs0HAADAEKHyAQCAbURO9SIUJB8AANgCO5wCAAAMCSofAADYgJMWnJJ8AABgC8551Ja2CwAAMBWVDwAAbME5lQ+SDwAAbMBJaz5ouwAAAFNR+QAAwBZouwAAAFM5Z5Mxkg8AAGyANR8AAABDhMoHAAC2wJoPAABgKuckH7RdAACAqah8AABgA4bRJ8MY/BMrocw1G5UPAABswQjDCFxpaaluv/12xcfHKzExUfPmzdORI0cuOWfbtm2aM2eORo8eLbfbrezsbL399ttB3Vci+QAAwJHq6upUWFio/fv3q6amRj09PcrJyVFnZ+dF59TX12vOnDmqrq7WwYMHNXPmTOXl5amxsTGoe7uMSHowGACAYcbn88nj8ejPf35HbvfVIVynQzfccJe8Xq/cbnfQ87/44gslJiaqrq5O06ZNC3je+PHjtWDBAq1YsSLgOaz5AADAFsLztIvP5xtwNDY2VrGxsZed7fV6JUkJCQkB37Gvr09nzpwJao5E2wUAgGElLS1NHo/HP0pLSy87xzAMFRcX684779SECRMCvtfzzz+vzs5O3XfffUHFSOUDAABbMBTa+1nOVT5aW1sHtF0CqXosWbJEhw4d0t69ewO+25YtW7Ry5Urt3LlTiYmJQUVK8gEAgA2E690ubrc7qDUfRUVF2rVrl+rr6zVmzJiA5mzdulWLFi3SG2+8odmzZwcdK8kHAAC2YO4Op4ZhqKioSNu3b1dtba0yMjICmrdlyxY99NBD2rJli+65557BBEryAQCAExUWFuq1117Tzp07FR8fr/b2dkmSx+PRyJEjJUklJSU6ceKENm3aJOlc4pGfn68XX3xRkydP9s8ZOXKkPB5PwPdmwSkAALZg7iZjFRUV8nq9mjFjhlJSUvxj69at/nPa2trU0tLi/7x+/Xr19PSosLBwwJzHHnssqHtT+QAAwBZC21492MWqgawvqaysHPC5trY2qHtcDJUPAABgKiofAADYgrkLTq1E8gEAgC04J/mg7QIAAExF5QMAABswjNAWnIa2WNVcJB8AANgCbRcAAIAhQeUDAABbcE7lg+QDAAAbYM0HAAAwmXMqH6z5AAAApqLyAQCADRiGEdD7Vi41P1KQfAAAYAuGgn053PnzIwNtFwAAYCoqHwAA2IJzFpySfAAAYANOWvNB2wUAAJiKygcAAHZgGOdGKPMjBJUPAABgKpIPAABgKtouAADYgJMWnJJ8AABgB8550pbkAwAAO3BS5YM1HwAAwFRUPgAAsAPaLgAAwFTs8wEAADA0qHwAAGADDip8kHwAAGALDso+aLsAAABTkXwAAABT0XYBAMAG2GQMAABgiFD5AADADthkDAAAmMpBT7uQfAAAYAMOyj1Y8wEAgBOVlpbq9ttvV3x8vBITEzVv3jwdOXLksvPq6uqUmZmpuLg43XjjjVq3bl3Q9yb5AADADvpLH6GMINTV1amwsFD79+9XTU2Nenp6lJOTo87OzovOaW5u1ty5czV16lQ1NjZq+fLlWrp0qaqqqoK6t8uIpGdzAAAYZnw+nzwejz76YJPir75y0Nc50/FXjf8/+fJ6vXK73UHP/+KLL5SYmKi6ujpNmzbtguf88z//s3bt2qWPP/7Yf6ygoEAffvih9u3bF/C9HLfmo6+vT5999pni4+PlcrmsDgcAYGOGYejMmTNKTU3ViBGR0Szw+XwDPsfGxio2Nvay87xeryQpISHhoufs27dPOTk5A47dfffd2rhxo86ePasrrrgioBgdl3x89tlnSktLszoMAEAEaW1t1ZgxY4b2JmFacfrN33HPPPOMVq5ceZmphoqLi3XnnXdqwoQJFz2vvb1dSUlJA44lJSWpp6dHp06dUkpKSkChOi75iI+Pl3Tu/0iDKUsBAJzD5/MpLS3N/7tjSIVpn49v/n4LpOqxZMkSHTp0SHv37r3sud/sGvSv3gimm+C45KP/y3G73SQfAICARFKbPtjfb0VFRdq1a5fq6+svW91JTk5We3v7gGMnT55UdHS0Ro0aFfA9HZd8AABgR8bf/hfK/KDONwwVFRVp+/btqq2tVUZGxmXnZGdn66233hpwbM+ePcrKygp4vYfEo7YAANiDEYYRhMLCQr366qt67bXXFB8fr/b2drW3t+vLL7/0n1NSUqL8/Hz/54KCAn366acqLi7Wxx9/rJdfflkbN27UsmXLgro3yQcAAA5UUVEhr9erGTNmKCUlxT+2bt3qP6etrU0tLS3+zxkZGaqurlZtba1uu+02/eQnP9FLL72k+fPnB3Vv2i4AANiByfurB7LNV2Vl5XnHpk+frg8++CCoe30TyQcAADbgpHe7kHwAAGAHDso+WPMBAABMReUDAAAbcFDhg+QDAABbcFD2QdsFAACYisoHAAB2EKZ3u0QCkg8AAGzAMIyA9t641PxIQdsFAACYisoHAAB2QNsFAACYibYLAADAEKHyAQCAHdB2AQAApnLQJmMkHwAA2ICDcg/WfAAAAHNR+QAAwA4MhVj6CFskQ47kAwAAO3DQglPaLgAAwFSWJh/19fXKy8tTamqqXC6XduzYccnza2tr5XK5zht//OMfzQkYAIAhYsjwbzQ2qBFBpQ9L2y6dnZ2aOHGiHnzwQc2fPz/geUeOHJHb7fZ/Hj16dND3zsvMVHRUVNDzwqWvr8+yeweip7fX6hAuqdfm8fXY+O/X7t/d2a4uq0O4pLPd3VaHcFF8d+Fn6j/JDmq7WJp85ObmKjc3N+h5iYmJuuaaa8IfEAAAGHIRueZj0qRJSklJ0axZs/Tb3/72kud2dXXJ5/MNGAAA2E4oLZdQNwkxWUQlHykpKdqwYYOqqqq0bds2jR07VrNmzVJ9ff1F55SWlsrj8fhHWlqaiREDABAgIwwjQkTUo7Zjx47V2LFj/Z+zs7PV2tqq1atXa9q0aRecU1JSouLiYv9nn89HAgIAgIUiKvm4kMmTJ+vVV1+96J/HxsYqNjbWxIgAABgEB+2vHvHJR2Njo1JSUqwOAwCAkDgo97A2+ejo6NCxY8f8n5ubm9XU1KSEhASlp6erpKREJ06c0KZNmyRJZWVluuGGGzR+/Hh1d3fr1VdfVVVVlaqqqqz6EQAACA8HZR+WJh8NDQ2aOXOm/3P/2oyFCxeqsrJSbW1tamlp8f95d3e3li1bphMnTmjkyJEaP368fvOb32ju3Lmmxw4AAAbHZRgRlCqFgc/nk8fj0bRvf5tNxi6BTcZCwyZjg8dGWYPHdxd+fZKOS/J6vQM2twyn/t9LDdVrdfVVIwd9nY7OL5U199EhjTVcIn7NBwAAw4KD2i4Rtc8HAACIfFQ+AACwA97tAgAAzOTfJj2E+ZGCtgsAADAVlQ8AAOyABacAAMBM/blHKCNY9fX1ysvLU2pqqlwul3bs2HHZOZs3b9bEiRN15ZVXKiUlRQ8++KBOnz4d1H1JPgAAcKjOzk5NnDhR5eXlAZ2/d+9e5efna9GiRfroo4/0xhtv6MCBA1q8eHFQ96XtAgCAHVjQdsnNzVVubm7A5+/fv1833HCDli5dKknKyMjQI488op///OdB3ZfKBwAAdmCEYejcjqlfH11h3Pl2ypQpOn78uKqrq2UYhj7//HO9+eabuueee4K6DskHAAA20P+obShDktLS0uTxePyjtLQ0bDFOmTJFmzdv1oIFCxQTE6Pk5GRdc801+s///M+grkPyAQDAMNLa2iqv1+sfJSUlYbv2H/7wBy1dulQrVqzQwYMHtXv3bjU3N6ugoCCo67DmAwAAOwjTDqdut3vIXixXWlqqO+64Q08++aQk6dZbb9VVV12lqVOnatWqVUpJSQnoOiQfAADYQCTscPrXv/5V0dEDU4eov70hPpj703YBAMChOjo61NTUpKamJklSc3Ozmpqa1NLSIkkqKSlRfn6+//y8vDxt27ZNFRUV+uSTT/Tee+9p6dKl+od/+AelpqYGfF8qHwAA2IXJm5Q2NDRo5syZ/s/FxcWSpIULF6qyslJtbW3+RESSfvSjH+nMmTMqLy/Xj3/8Y11zzTW666679LOf/Syo+5J8AABgBxbs8zFjxoxLtksqKyvPO1ZUVKSioqKg7/V1tF0AAICpqHwAAGADkbDgNFxIPgAAsIMwPWobCWi7AAAAU1H5AADABpzUdrG08lFfX6+8vDylpqbK5XJpx44dl51TV1enzMxMxcXF6cYbb9S6deuGPlAAAIZaXxhGhLA0+ejs7NTEiRNVXl4e0PnNzc2aO3eupk6dqsbGRi1fvlxLly5VVVXVEEcKAMAQ63/UNpQRISxtu+Tm5io3Nzfg89etW6f09HSVlZVJksaNG6eGhgatXr1a8+fPH6IoAQBAOAVc+Th+/PhQxhGQffv2KScnZ8Cxu+++Ww0NDTp79uwF53R1dcnn8w0YAADYjYMKH4EnHxMmTNCvf/3roYzlstrb25WUlDTgWFJSknp6enTq1KkLziktLZXH4/GPtLQ0M0IFACA4Dso+Ak4+fvrTn6qwsFDz58/X6dOnhzKmS3K5XAM+96/u/ebxfiUlJfJ6vf7R2to65DECAICLCzj5ePTRR/Xhhx/qf//3fzV+/Hjt2rVrKOO6oOTkZLW3tw84dvLkSUVHR2vUqFEXnBMbGyu32z1gAABgNw4qfAS34DQjI0PvvPOOysvLNX/+fI0bN07R0QMv8cEHH4Q1wK/Lzs7WW2+9NeDYnj17lJWVpSuuuGLI7gsAwJCz4MVyVgn6aZdPP/1UVVVVSkhI0L333nte8hGMjo4OHTt2zP+5ublZTU1NSkhIUHp6ukpKSnTixAlt2rRJklRQUKDy8nIVFxfr4Ycf1r59+7Rx40Zt2bJl0DEAAABzBZU5/OIXv9CPf/xjzZ49W7///e81evTokG7e0NCgmTNn+j8XFxdLkhYuXKjKykq1tbWppaXF/+cZGRmqrq7WE088oTVr1ig1NVUvvfQSj9kCACIflY/z/eM//qN+97vfqby8XPn5+WG5+YwZMy65HWxlZeV5x6ZPnz6krR0AAKxg9J0bocyPFAEnH729vTp06JDGjBkzlPEAAIBhLuDko6amZijjAADA2QyF2HYJWyRDjrfaAgBgAw5a8kHyAQCALTgo+7D0rbYAAMB5qHwAAGAHDqp8kHwAAGADhhHio7aRk3vQdgEAAOai8gEAgB3QdgEAAGZyUO5B2wUAAJiLygcAAHbgoNIHyQcAADZgGMYlX7YayPxIQdsFAACYisoHAAB20Pe3Ecr8CEHyAQCAHbDmAwAAmMlBuQdrPgAAgLmofAAAYAd9xrkRyvwIQfIBAIAN8KgtAADAECH5AADADowwjCDV19crLy9Pqampcrlc2rFjx2XndHV16emnn9b111+v2NhYfetb39LLL78c1H0tTz7Wrl2rjIwMxcXFKTMzU+++++5Fz62trZXL5Tpv/PGPfzQxYgAAhkD/4y6hjCB1dnZq4sSJKi8vD3jOfffdp//+7//Wxo0bdeTIEW3ZskXf+c53grqvpWs+tm7dqscff1xr167VHXfcofXr1ys3N1d/+MMflJ6eftF5R44ckdvt9n8ePXq0GeECADCs5ObmKjc3N+Dzd+/erbq6On3yySdKSEiQJN1www1B39fSyscLL7ygRYsWafHixRo3bpzKysqUlpamioqKS85LTExUcnKyf0RFRZkUMQAAQ6TPkBHC6H/axefzDRhdXV1hC3HXrl3KysrSz3/+c1133XW6+eabtWzZMn355ZdBXcey5KO7u1sHDx5UTk7OgOM5OTl6//33Lzl30qRJSklJ0axZs/Tb3/72kud2dXWd9xcBAIDthGnNR1pamjwej3+UlpaGLcRPPvlEe/fu1e9//3tt375dZWVlevPNN1VYWBjUdSxru5w6dUq9vb1KSkoacDwpKUnt7e0XnJOSkqINGzYoMzNTXV1d+vWvf61Zs2aptrZW06ZNu+Cc0tJSPfvss2GPHwAAO2ptbR2wNCE2NjZs1+7r65PL5dLmzZvl8Xgkneti/NM//ZPWrFmjkSNHBnQdy/f5cLlcAz4bhnHesX5jx47V2LFj/Z+zs7PV2tqq1atXXzT5KCkpUXFxsf+zz+dTWlpaGCIHACB8DIW4z8ffSh9ut3tA8hFOKSkpuu666/yJhySNGzdOhmHo+PHjuummmwK6jmVtl2uvvVZRUVHnVTlOnjx5XjXkUiZPnqyjR49e9M9jY2P9fxFD+RcCAEBI+sIwhtgdd9yhzz77TB0dHf5j//M//6MRI0ZozJgxAV/HsspHTEyMMjMzVVNTo+9973v+4zU1Nbr33nsDvk5jY6NSUlKCvv9bBw+SiAAALsnn8w34r/yhZMUOpx0dHTp27Jj/c3Nzs5qampSQkKD09HSVlJToxIkT2rRpkyTp/vvv109+8hM9+OCDevbZZ3Xq1Ck9+eSTeuihhwJuuUgWt12Ki4v1wAMPKCsrS9nZ2dqwYYNaWlpUUFAgSef90GVlZbrhhhs0fvx4dXd369VXX1VVVZWqqqqs/DEAAIhIDQ0Nmjlzpv9z/zKFhQsXqrKyUm1tbWppafH/+dVXX62amhoVFRUpKytLo0aN0n333adVq1YFdV9Lk48FCxbo9OnTeu6559TW1qYJEyaourpa119/vSSd90N3d3dr2bJlOnHihEaOHKnx48frN7/5jebOnWvVjwAAQHhY8GK5GTNmXLJiUllZed6x73znO6qpqQn6Xl/nMiLpTTRh0F9C83q9tF0AAJdkxu+M/nvU/ORfdFVc3KCv0/nVV5rz/1ZFxO83y7dXBwAAzmL5o7YAAMCaBadWIfkAAMAOLFjzYRXaLgAAwFRUPgAAsAHDODdCmR8pSD4AALAD2i4AAABDg8oHAAA2wNMuAADAXKG+HM6EF8uFC8kHAAA2cG7BaSiVjzAGM8RY8wEAAExF5QMAADswQnzaJYJKHyQfAADYgfG3Ecr8CEHbBQAAmIrKBwAANsCjtgAAwFROSj5ouwAAAFNR+QAAwAYctMcYyQcAAHZA2wUAAGCIUPkAAMAGnFT5IPkAAMAGHLTHmPVtl7Vr1yojI0NxcXHKzMzUu+++e8nz6+rqlJmZqbi4ON14441at26dSZECADB0+isfoYxIYWnysXXrVj3++ON6+umn1djYqKlTpyo3N1ctLS0XPL+5uVlz587V1KlT1djYqOXLl2vp0qWqqqoyOXIAADBYliYfL7zwghYtWqTFixdr3LhxKisrU1pamioqKi54/rp165Senq6ysjKNGzdOixcv1kMPPaTVq1ebHDkAAOFF5cME3d3dOnjwoHJycgYcz8nJ0fvvv3/BOfv27Tvv/LvvvlsNDQ06e/bsBed0dXXJ5/MNGAAA2A3JhwlOnTql3t5eJSUlDTielJSk9vb2C85pb2+/4Pk9PT06derUBeeUlpbK4/H4R1paWnh+AAAAMCiWLzh1uVwDPhuGcd6xy51/oeP9SkpK5PV6/aO1tTXEiAEACD8jDCNSWPao7bXXXquoqKjzqhwnT548r7rRLzk5+YLnR0dHa9SoURecExsbq9jYWP/n/mSF9gsA4HL6f1eY0tIItXUSQW0Xy5KPmJgYZWZmqqamRt/73vf8x2tqanTvvfdecE52drbeeuutAcf27NmjrKwsXXHFFQHd98yZM5JE+wUAELAzZ87I4/FYHcawYekmY8XFxXrggQeUlZWl7OxsbdiwQS0tLSooKJB0rmVy4sQJbdq0SZJUUFCg8vJyFRcX6+GHH9a+ffu0ceNGbdmyJeB7pqamqrW1VfHx8XK5XPL5fEpLS1Nra6vcbveQ/JyRgO/hHL6Hc/gezuF7OMfJ34NhGDpz5oxSU1NNuRc7nJpgwYIFOn36tJ577jm1tbVpwoQJqq6u1vXXXy9JamtrG7DnR0ZGhqqrq/XEE09ozZo1Sk1N1UsvvaT58+cHfM8RI0ZozJgx5x13u92O+4fqQvgezuF7OIfv4Ry+h3Oc+j2YVfFw0g6nlm+v/uijj+rRRx+94J9VVlaed2z69On64IMPhjgqAAAwVCxPPgAAAG0XR4mNjdUzzzwz4IkYJ+J7OIfv4Ry+h3P4Hs7hezCHk5IPlxFJ0QIAMMz4fD55PB5tKSzUlSEkeH/t6tIP1qyR1+sNeG1OfX29/v3f/10HDx5UW1ubtm/frnnz5gU097333tP06dM1YcIENTU1BRWr5ZuMAQAAa3R2dmrixIkqLy8Pap7X61V+fr5mzZo1qPs6vu0CAIAd9P1thDI/WLm5ucrNzQ163iOPPKL7779fUVFR2rFjR9DzqXwAAGAD4Xqx3DdfptrV1RXWOF955RX96U9/0jPPPDPoa5B8AAAwjKSlpQ14oWppaWnYrn306FE99dRT2rx5s6KjB988oe0CAIANhOtpl2/uRBuup5R6e3t1//3369lnn9XNN98c0rUcXflYu3atMjIyFBcXp8zMTL377rtWh2S6+vp65eXlKTU1VS6Xa1C9u+GgtLRUt99+u+Lj45WYmKh58+bpyJEjVodluoqKCt16663+nSyzs7P1X//1X1aHZanS0lK5XC49/vjjVodiupUrV8rlcg0YycnJVoc1bIWr7dL/z2//CFfycebMGTU0NGjJkiWKjo5WdHS0nnvuOX344YeKjo7WO++8E/C1HJt8bN26VY8//riefvppNTY2aurUqcrNzR2wnbsTDHal83BTV1enwsJC7d+/XzU1Nerp6VFOTo46OzutDs1UY8aM0b/927+poaFBDQ0Nuuuuu3Tvvffqo48+sjo0Sxw4cEAbNmzQrbfeanUolhk/frza2tr84/Dhw1aHBIu43W4dPnxYTU1N/lFQUKCxY8eqqalJ3/3udwO+lmPbLi+88IIWLVqkxYsXS5LKysr09ttvq6KiIqz9Mbsb7Ern4Wb37t0DPr/yyitKTEzUwYMHNW3aNIuiMl9eXt6Az//6r/+qiooK7d+/X+PHj7coKmt0dHTohz/8oX7xi19o1apVVodjmejoaKodJrHi3S4dHR06duyY/3Nzc7OampqUkJCg9PT0AS94HTFihCZMmDBgfmJiouLi4s47fjmOrHx0d3fr4MGDysnJGXA8JydH77//vkVRwU68Xq8kKSEhweJIrNPb26vXX39dnZ2dys7Otjoc0xUWFuqee+7R7NmzrQ7FUkePHlVqaqoyMjL0/e9/X5988onVIQ1b4Wq7BKOhoUGTJk3SpEmTJJ172/ykSZO0YsUKSee/4DVcHFn5OHXqlHp7e5WUlDTgeFJSktrb2y2KCnZhGIaKi4t15513Bp3NDweHDx9Wdna2vvrqK1199dXavn27/v7v/97qsEz1+uuv64MPPtCBAwesDsVS3/3ud7Vp0ybdfPPN+vzzz7Vq1SpNmTJFH330kUaNGmV1eMOOodC2SB/MzBkzZlzynhd6wevXrVy5UitXrgz6vo5MPvq5XK4Bnw3DOO8YnGfJkiU6dOiQ9u7da3Uolujv3/7lL39RVVWVFi5cqLq6OsckIK2trXrssce0Z88excXFWR2Opb7ekr3llluUnZ2tb33rW/rVr36l4uJiCyNDpHNk8nHttdcqKirqvCrHyZMnz6uGwFmKioq0a9cu1dfXa8yYMVaHY4mYmBh9+9vfliRlZWXpwIEDevHFF7V+/XqLIzPHwYMHdfLkSWVmZvqP9fb2qr6+XuXl5erq6lJUVJSFEVrnqquu0i233KKjR49aHcqw1GcY6guh8hHKXLM5cs1HTEyMMjMzVVNTM+B4TU2NpkyZYlFUsJJhGFqyZIm2bdumd955RxkZGVaHZBuGYYR9h0Q7mzVr1nkr+rOysvTDH/5QTU1Njk08JKmrq0sff/yxUlJSrA5lWDLCMCKFIysf0rlFNQ888ICysrKUnZ2tDRs2qKWlRQUFBVaHZqrLrXR2isLCQr322mvauXOn4uPj/VUxj8ejkSNHWhydeZYvX67c3FylpaXpzJkzev3111VbW3ve00DDWXx8/Hlrfa666iqNGjXKcWuAli1bpry8PKWnp+vkyZNatWqVfD6fFi5caHVoiHCOTT4WLFig06dP67nnnlNbW5smTJig6upqXX/99VaHZqqGhgbNnDnT/7m/j7tw4cLLLjQaTioqKiSdW3z1da+88op+9KMfmR+QRT7//HM98MADamtrk8fj0a233qrdu3drzpw5VocGCxw/flw/+MEPdOrUKY0ePVqTJ0/W/v37HffvSbOEa4fTSOAyIilaAACGGZ/PJ4/How2LF+vKmJhBX+ev3d36v7/8pbxe74Dt1e3IkWs+AACAdRzbdgEAwFZCbLsoghoZJB8AANgAj9oCAAAMESofAADYgJOediH5AADABqx4q61VSD4AALABJ1U+WPMBAABMRfIBOFRvb6+mTJmi+fPnDzju9XqVlpamf/mXf7EoMsCZ+p92CWVECpIPwKGioqL0q1/9Srt379bmzZv9x4uKipSQkKAVK1ZYGB3gPP1tl1BGpGDNB+BgN910k0pLS1VUVKSZM2fqwIEDev311/W73/1OMSFs8wwAl0LyAThcUVGRtm/frvz8fB0+fFgrVqzQbbfdZnVYgOM4acEpyQfgcC6XSxUVFRo3bpxuueUWPfXUU1aHBDiSk5IP1nwA0Msvv6wrr7xSzc3NOn78uNXhABjmSD4Ah9u3b5/+4z/+Qzt37lR2drYWLVoUUf8FBQwXfWEYkYLkA3CwL7/8UgsXLtQjjzyi2bNn65e//KUOHDig9evXWx0a4DhOetqF5ANwsKeeekp9fX362c9+JklKT0/X888/ryeffFJ//vOfrQ0OwLBF8gE4VF1dndasWaPKykpdddVV/uMPP/ywpkyZQvsFMJmTKh887QI41PTp09XT03PBP3v77bdNjgZAqLuURtIOpyQfAADYAI/aAgAADBEqHwAA2EGo6zYiqPJB8gEAgA04ac0HbRcAAGAqKh8AANiAodAWjUZO3YPkAwAAW+BpFwAAgCFC5QMAABtw0oJTkg8AAGyAtgsAAMAQofIBAIANOKntQuUDAAAbsOKttvX19crLy1NqaqpcLpd27NhxyfO3bdumOXPmaPTo0XK73crOzh7UiyhJPgAAsIH+ykcoI1idnZ2aOHGiysvLAzq/vr5ec+bMUXV1tQ4ePKiZM2cqLy9PjY2NQd2XtgsAAA6Vm5ur3NzcgM8vKysb8PmnP/2pdu7cqbfeekuTJk0K+DokHwAA2EC4nnbx+XwDjsfGxio2Njak2C6mr69PZ86cUUJCQlDzaLsAAGAD4Wq7pKWlyePx+EdpaemQxfz888+rs7NT9913X1DzqHwAADCMtLa2yu12+z8PVdVjy5YtWrlypXbu3KnExMSg5pJ8AABgA+Fqu7jd7gHJx1DYunWrFi1apDfeeEOzZ88Oej7JBwAANmCEuM+HWTucbtmyRQ899JC2bNmie+65Z1DXIPkAAMChOjo6dOzYMf/n5uZmNTU1KSEhQenp6SopKdGJEye0adMmSecSj/z8fL344ouaPHmy2tvbJUkjR46Ux+MJ+L4sOAUAwAas2GSsoaFBkyZN8j8mW1xcrEmTJmnFihWSpLa2NrW0tPjPX79+vXp6elRYWKiUlBT/eOyxx4K6L5UPAABswIrt1WfMmHHJpKWysnLA59ra2qDvcSFUPgAAgKmofAAAYAPhetolEpB8AABgA056qy3JBwAANuCk5IM1HwAAwFRUPgAAsAHWfAAAAFPRdgEAABgiVD4AALABJ1U+SD4AALABwzBk9PWFND9S0HYBAACmovIBAIAN0HYBAACmMkJMPmi7AAAAXASVDwAAbIC2CwAAMFVfX5/6QnjaJZS5ZiP5AADABpy0vTprPgAAgKmofAAAYAOs+QAAAKZy0poP2i4AAMBUVD4AALAB2i4AAMBUTko+aLsAAABTUfkAAMAGnLTglOQDAAAb6DMM9dJ2AQAACD8qHwAA2ABtFwAAYKo+hdY6iZzUg+QDAABb6OvrU5/LFdL8SMGaDwAAYCoqHwAA2EBvX59GhFD56I2gygfJBwAANsAOpwAAAEOEygcAADZA2wUAAJiKp10AAMCwV19fr7y8PKWmpsrlcmnHjh2XnVNXV6fMzEzFxcXpxhtv1Lp164K+L8kHAAA20NfXp94QxmAqH52dnZo4caLKy8sDOr+5uVlz587V1KlT1djYqOXLl2vp0qWqqqoK6r60XQAAsIHevj65TF7zkZubq9zc3IDPX7dundLT01VWViZJGjdunBoaGrR69WrNnz8/4OtQ+QAAYBjx+XwDRldXV9iuvW/fPuXk5Aw4dvfdd6uhoUFnz54N+DokHwAA2EAoLZf+IUlpaWnyeDz+UVpaGrYY29vblZSUNOBYUlKSenp6dOrUqYCvQ9sFAAAbCNfTLq2trXK73f7jsbGxIcf2dd9sDRl/29wsmJYRyQcAADbQE+Kjsv3z3W73gOQjnJKTk9Xe3j7g2MmTJxUdHa1Ro0YFfB3aLgAAICDZ2dmqqakZcGzPnj3KysrSFVdcEfB1SD4AALCB3t7ekEewOjo61NTUpKamJknnHqVtampSS0uLJKmkpET5+fn+8wsKCvTpp5+quLhYH3/8sV5++WVt3LhRy5YtC+q+tF0AALCBULdHH8z8hoYGzZw50/+5uLhYkrRw4UJVVlaqra3Nn4hIUkZGhqqrq/XEE09ozZo1Sk1N1UsvvRTUY7aS5DKMCHoNHgAAw4zP55PH49H/ue46RY0YfEOit69PH5w4Ia/XO2RrPsKFygcAADbQ29srhVAP4MVyAAAgKGe7utQXYuUjUrDgFAAAmIrKBwAANnD27NmQNhnrjaAlnCQfAADYQE9XV2g7nEZQ8kHbBQAAmIrKBwAANtB19mxIFYHIWW5K8gEAgKViYmKUnJysz77xzpTBSE5OVkxMTBiiGlpsMgYAgMW++uordXd3h3ydmJgYxcXFhSGioUXyAQAATMWCUwAAYCqSDwAAYCqSDwAAYCqSDwAAYCqSDwAAYCqSDwAAYCqSDwAAYKr/D+1vf7OUHoSJAAAAAElFTkSuQmCC",
+ "text/plain": [
+ "