diff --git a/.github/workflows/notebooks.yml b/.github/workflows/notebooks.yml index 1e5fdcf..e30e8e1 100644 --- a/.github/workflows/notebooks.yml +++ b/.github/workflows/notebooks.yml @@ -33,5 +33,5 @@ jobs: - name: Install dependencies run: pip install nox - - name: Build documentation + - name: Test notebooks run: nox -s test-notebooks diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index c968a56..7657d14 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -20,7 +20,7 @@ jobs: strategy: matrix: os: [ubuntu-latest, macos-latest, windows-latest] - python-version: ["3.9", "3.10", "3.11"] + python-version: ["3.10", "3.11"] steps: - uses: actions/checkout@v2 diff --git a/docs/reference/cli/index.md b/docs/reference/cli/index.md index 2e868de..d8e492c 100644 --- a/docs/reference/cli/index.md +++ b/docs/reference/cli/index.md @@ -1,4 +1,4 @@ -# ``sequence`` +# Command-line interface The *Sequence* model can be run from the command line using the *sequence* program. diff --git a/docs/reference/cli/input_files.md b/docs/reference/cli/input_files.md index 71d77b9..e0b73a1 100644 --- a/docs/reference/cli/input_files.md +++ b/docs/reference/cli/input_files.md @@ -58,11 +58,34 @@ spacing = 1000.0 ``` In this case we have a grid that represents a 1D profile that consists of -500 columns (i.e. vertical stacks) of sediment (the *n_cols* parameter). +100 columns (i.e. vertical stacks) of sediment (the *n_cols* parameter). The *spacing* parameter is the width of each of your sediment stacks in meters. Thus, the length of you domain is the product of the number of columns with -the spacing (that is, for this example, 500 * 100 m or 50 km). +the spacing (that is, for this example, 100 * 1000 m or 100 km). + +You can also run *Sequence* in a quasi-2d mode in which *Sequence* models parallel +cross-shore profiles. Each profile will have its own sediment supply that is +transported into the ocean and then can be diffused between profiles. To set up +such a case, the grid section can, for example, be modified in the following way, + +```toml +[sequence.grid] +shape = [3, 100] +spacing = [10000.0, 1000.0] +``` +In this case *Sequence* will create three parallel profiles (with an along-shore +width of 10000.0 m). + +:::{note} +The following is equivalent to the 1D example above, + +```toml +[sequence.grid] +shape = [1, 100] +spacing = [1.0, 1000.0] +``` +::: (the-output-section)= diff --git a/docs/tutorials/cli/sealevel.md b/docs/tutorials/cli/sealevel.md index 8922851..e17ee09 100644 --- a/docs/tutorials/cli/sealevel.md +++ b/docs/tutorials/cli/sealevel.md @@ -27,7 +27,7 @@ we replace this section with the following, ```toml [sequence.sea_level] -filepath = sealevel.csv +filepath = "sealevel.csv" ``` With this configuration, *Sequence* will read sea level values from this diff --git a/news/78.feature b/news/78.feature new file mode 100644 index 0000000..dd4171d --- /dev/null +++ b/news/78.feature @@ -0,0 +1,3 @@ + +Added a quasi-2d mode for *Sequence* that keeps track of multiple cross-shore +rows. diff --git a/sequence/_grid.py b/sequence/_grid.py index a422259..deb4008 100644 --- a/sequence/_grid.py +++ b/sequence/_grid.py @@ -15,26 +15,62 @@ class SequenceModelGrid(RasterModelGrid): """Create a Landlab ModelGrid for use with Sequence.""" - def __init__(self, n_cols: int, spacing: float = 100.0): + def __init__( + self, shape: int | tuple[int, int], spacing: float | tuple[float, float] = 100.0 + ): """Create a *Landlab* :class:`~landlab.grid.base.ModelGrid` for use with *Sequence*. Parameters ---------- - n_cols : int - The number of columns. - spacing : float, optional - The spacing between columns. + shape : int or tuple of int + The number of columns in the cross-shore direction or, if a ``tuple``, + ``(n_rows, n_cols)`` where rows are in the along-shore direction and + columns in the cross-shore direction. + spacing : float or tuple of float, optional + The spacing between columns or, if ``len(shape) == 2``, + ``(row_spacing, col_spacing)``. Examples -------- >>> from sequence import SequenceModelGrid - >>> grid = SequenceModelGrid(500, spacing=10.0) - >>> grid.shape - (3, 500) - >>> grid.spacing - (1.0, 10.0) + >>> grid = SequenceModelGrid(5, spacing=10.0) + >>> grid.y_of_row + array([ 0., 1., 2.]) + >>> grid.x_of_column + array([ 0., 10., 20., 30., 40.]) + + >>> grid = SequenceModelGrid((3, 5), spacing=(10000.0, 10.0)) + >>> grid.y_of_row + array([ 0., 10000., 20000., 30000., 40000.]) + >>> grid.x_of_column + array([ 0., 10., 20., 30., 40.]) """ - super().__init__((3, n_cols), xy_spacing=(spacing, 1.0)) + row_col_shape: list[int] = np.atleast_1d(np.asarray(shape, dtype=int)) + row_col_spacing: list[float] = np.atleast_1d(np.asarray(spacing, dtype=float)) + + if (len(row_col_shape) == 1 and len(row_col_spacing) != 1) or ( + len(row_col_shape) != 1 and len(row_col_spacing) != len(row_col_shape) + ): + raise ValueError( + f"spacing dimension ({len(row_col_spacing)}) does not match shape" + f" dimension ({len(row_col_shape)})" + ) + + if len(row_col_shape) == 1: + n_rows, n_cols = 1, row_col_shape[0] + elif len(row_col_shape) == 2: + n_rows, n_cols = row_col_shape + else: + raise ValueError( + f"invalid number of dimensions for grid ({len(row_col_shape)})" + ) + + if len(row_col_shape) == 1: + row_spacing, col_spacing = 1.0, row_col_spacing[0] + elif len(row_col_shape) == 2: + row_spacing, col_spacing = row_col_spacing + + super().__init__((n_rows + 2, n_cols), xy_spacing=(col_spacing, row_spacing)) self.status_at_node[self.nodes_at_top_edge] = self.BC_NODE_IS_CLOSED self.status_at_node[self.nodes_at_bottom_edge] = self.BC_NODE_IS_CLOSED @@ -42,6 +78,30 @@ def __init__(self, n_cols: int, spacing: float = 100.0): self.at_node["sediment_deposit__thickness"] = self.zeros(at="node") self.at_grid["sea_level__elevation"] = 0.0 + # self.new_field_location("row", size=int(n_rows + 2)) + self.new_field_location("row", size=int(n_rows)) + + @property + def x_of_column(self) -> NDArray: + """X-coordinate for each column of the grid.""" + return self.x_of_node[self.nodes_at_top_edge] + + @property + def number_of_rows(self) -> int: + """Number of rows of nodes.""" + return self.shape[0] + # return self.shape[0] - 2 + + @property + def number_of_columns(self) -> int: + """Number of columns of nodes.""" + return self.shape[1] + + @property + def y_of_row(self) -> NDArray: + """Y-coordinate for each row of the grid.""" + return self.y_of_node[self.nodes_at_left_edge] + def get_profile(self, name: str) -> NDArray: """Return the values of a field along the grid's profile. @@ -55,7 +115,8 @@ def get_profile(self, name: str) -> NDArray: values : ndarray The values of the field located at the middle row of nodes. """ - return self.at_node[name].reshape(self.shape)[1] + return self.at_node[name].reshape(self.shape)[1:-1] + # return self.at_node[name].reshape(self.shape)[row] @classmethod def from_toml(cls, filepath: os.PathLike[str]) -> "SequenceModelGrid": @@ -85,19 +146,36 @@ def from_dict(cls, params: dict) -> "SequenceModelGrid": params : dict A dictionary that contains the parameters needed to create the grid. + + Examples + -------- + >>> from sequence import SequenceModelGrid + >>> params = {"shape": 5, "spacing": 10.0} + >>> grid = SequenceModelGrid.from_dict(params) + >>> grid.y_of_row + array([ 0., 1., 2.]) + >>> grid.x_of_column + array([ 0., 10., 20., 30., 40.]) + + >>> params = {"shape": (3, 5), "spacing": (10000.0, 10.0)} + >>> grid = SequenceModelGrid.from_dict(params) + >>> grid.y_of_row + array([ 0., 10000., 20000., 30000., 40000.]) + >>> grid.x_of_column + array([ 0., 10., 20., 30., 40.]) """ - if "n_cols" in params: - n_cols = params["n_cols"] - elif "shape" in params: - n_cols = params["shape"][1] + if "shape" in params: + shape = params["shape"] + elif "n_cols" in params: + shape = params["n_cols"] else: - raise KeyError("n_cols") + raise KeyError("shape") if "spacing" in params: spacing = params["spacing"] elif "xy_spacing" in params: - spacing = np.broadcast_to(params["xy_spacing"], 2)[0] + spacing = np.atleast_1d(params["xy_spacing"])[::-1] else: raise KeyError("spacing") - return cls(n_cols, spacing=spacing) + return cls(shape, spacing=spacing) diff --git a/sequence/cli.py b/sequence/cli.py index bc91e20..87b7d78 100644 --- a/sequence/cli.py +++ b/sequence/cli.py @@ -17,7 +17,7 @@ from numpy.typing import ArrayLike from tqdm import tqdm -from .errors import MissingRequiredVariable +from .errors import OutputValidationError from .input_reader import TimeVaryingConfig from .logging import LoggingHandler from .plot import plot_file, plot_layers @@ -381,12 +381,20 @@ def setup(set: str) -> None: @sequence.command() -@click.option("--set", multiple=True, help="Set model parameters") +@click.option("--set", "-S", multiple=True, help="Set model parameters") +@click.argument( + "netcdf_file", + type=click.Path( + exists=True, file_okay=True, dir_okay=False, path_type=pathlib.Path + ), + nargs=1, +) @click.pass_context -def plot(ctx: Any, set: str) -> None: +def plot(ctx: Any, set: str, netcdf_file) -> None: """Plot a Sequence output file.""" verbose = ctx.parent.params["verbose"] folder = pathlib.Path.cwd() + path_to_file = folder / netcdf_file config = PLOT_KEYWORDS.copy() @@ -398,7 +406,6 @@ def plot(ctx: Any, set: str) -> None: ) config.update(**_load_params_from_strings(set)) - if verbose and len(config) > 0: logger.info( os.linesep.join( @@ -409,14 +416,11 @@ def plot(ctx: Any, set: str) -> None: ) ) - logger.info(f"Plotting {folder / 'sequence.nc'}") + logger.info(f"Plotting {path_to_file!s}") try: - plot_file(folder / "sequence.nc", **config) - except MissingRequiredVariable as error: - logger.error( - f"{folder / 'sequence.nc'}: output file is missing a required variable " - f"({error})" - ) + plot_file(path_to_file, **config) + except OutputValidationError as error: + logger.error(f"{path_to_file!s}: output file is invalid ({error!s})") raise click.Abort() from error diff --git a/sequence/errors.py b/sequence/errors.py index 92ee018..c6b6a45 100644 --- a/sequence/errors.py +++ b/sequence/errors.py @@ -29,15 +29,33 @@ def __str__(self) -> str: return self._msg -class MissingRequiredVariable(SequenceError): +class OutputValidationError(SequenceError): + """Raise this exception if there is something wrong with an output file.""" + + pass + + +class InvalidRowError(OutputValidationError): + """Raise this exception if a required variable is missing from an output file.""" + + def __init__(self, row: int, n_rows: int): + self._row = row + self._n_rows = n_rows + + def __str__(self) -> str: + """Return error message that includes the requested row.""" + return f"row {self._row!r} is out of bounds for output file with {self._n_rows} rows" + + +class MissingRequiredVariable(OutputValidationError): """Raise this exception if a required variable is missing from an output file.""" def __init__(self, name: str): self._name = name - def __str_(self) -> str: + def __str__(self) -> str: """Return error message that includes the name of the missing variable.""" - return f"{self._name!r}" + return f"requireed variable ({self._name!r}) is missing from output file" class ParameterMismatchError(SequenceError): diff --git a/sequence/fluvial.py b/sequence/fluvial.py index 6469a00..d7b4bb6 100644 --- a/sequence/fluvial.py +++ b/sequence/fluvial.py @@ -122,7 +122,6 @@ def run_one_step(self, dt: float) -> None: dt : float The time step to update the component by. """ - # Upstream boundary conditions */ mud_vol = self.sediment_load * (1.0 - self.sand_frac) / self.sand_frac sand_vol = self.sediment_load qs = ( @@ -142,18 +141,22 @@ def run_one_step(self, dt: float) -> None: # channel_width = sand_vol * self.basin_width / qs / 31536000. channel_width = sand_vol / qs - x = self.grid.x_of_node.reshape(self.grid.shape)[1] + # x = self.grid.x_of_node.reshape(self.grid.shape)[1] + x = self.grid.x_of_column elevation = self.grid.get_profile("topographic__elevation") + middle_row = elevation.shape[0] // 2 shore = find_shoreline( - x, elevation, sea_level=self.grid.at_grid["sea_level__elevation"] + x, + elevation[middle_row], + sea_level=self.grid.at_grid["sea_level__elevation"], ) land = x < shore water = x >= shore # land = self.grid.x_of_node[self.grid.node_at_cell] < shore # slope = np.gradient(z[1, land]) / self.grid.dx - slope = np.gradient(elevation) / self.grid.dx + slope = np.gradient(elevation[middle_row]) / self.grid.dx # slp = np.zeros(1) # slp[0] = self.plain_slope # slope = np.concatenate((slp,slop)) @@ -163,7 +166,7 @@ def run_one_step(self, dt: float) -> None: # ) # channel_depth = np.zeros(self.grid.number_of_cells) - channel_depth = np.zeros_like(elevation) # use upsteam slope + channel_depth = np.zeros_like(elevation[middle_row]) # use upsteam slope # channel_depth[0] = ( # (self.sand_density - 1000.) / 1000. * self.sand_grain / self.plain_slope channel_depth[land] = ( @@ -178,13 +181,19 @@ def run_one_step(self, dt: float) -> None: # width_cb = channel_width/epsilon # Original: r_cb = (model.new_height[i]-model.height[i]+model.d_sl); - r_cb = self.grid.get_profile("sediment_deposit__thickness") * epsilon - dz = self.grid.get_profile("bedrock_surface__increment_of_elevation").copy() + r_cb = ( + self.grid.get_profile("sediment_deposit__thickness")[middle_row] * epsilon + ) + dz = self.grid.get_profile("bedrock_surface__increment_of_elevation")[ + middle_row + ].copy() # original: r_b = model.thickness[i]; - r_b = self.grid.get_profile("sediment_deposit__thickness").copy() + r_b = self.grid.get_profile("sediment_deposit__thickness")[middle_row].copy() # r_fp = np.zeros(self.grid.shape[1]) - r_fp = np.zeros_like(elevation) - percent_sand = self.grid.get_profile("delta_sediment_sand__volume_fraction") + r_fp = np.zeros_like(elevation[middle_row]) + percent_sand = self.grid.get_profile("delta_sediment_sand__volume_fraction")[ + middle_row + ] percent_sand.fill(0.0) for i in np.where(land)[0]: @@ -265,7 +274,8 @@ def run_one_step(self, dt: float) -> None: water_depth = ( self.grid.at_grid["sea_level__elevation"] - - self.grid.at_node["topographic__elevation"] + # - self.grid.at_node["topographic__elevation"] + - self.grid.get_profile("topographic__elevation")[middle_row] ) add_mud = np.zeros(self.grid.shape[1]) taper = 1 @@ -287,7 +297,9 @@ def run_one_step(self, dt: float) -> None: if add_mud[i] < 0.0: add_mud[i] = 0.0 - sediment_thickness = self.grid.get_profile("sediment_deposit__thickness") + sediment_thickness = self.grid.get_profile("sediment_deposit__thickness")[ + middle_row + ] sediment_thickness[water] += add_mud[water] np.divide( diff --git a/sequence/netcdf.py b/sequence/netcdf.py index 8c8a9a9..e230275 100644 --- a/sequence/netcdf.py +++ b/sequence/netcdf.py @@ -1,4 +1,5 @@ """Write a `SequenceModelGrid` to a NetCDF file.""" +import contextlib import os import warnings from collections import defaultdict @@ -6,6 +7,7 @@ from typing import Any, Iterable, Optional, Sequence, Union import netCDF4 as nc +import numpy as np from numpy.typing import NDArray from ._grid import SequenceModelGrid @@ -43,15 +45,28 @@ def _create_grid_dimension( ) -> Any: """Create grid dimensions for a netcdf file.""" if ids is None: - # ids = Ellipsis ids = slice(None) - if at not in ("node", "link", "patch", "corner", "face", "cell", "grid"): + if at not in ( + "node", + "link", + "patch", + "corner", + "face", + "cell", + "grid", + "row", + "column", + ): raise ValueError(f"unknown grid location {at}") if at not in root.dimensions: if at == "grid": number_of_elements = 1 + elif at == "row": + number_of_elements = len(getattr(grid, f"y_of_{at}")[ids]) + elif at == "column": + number_of_elements = len(getattr(grid, f"x_of_{at}")[ids]) else: number_of_elements = len(getattr(grid, f"xy_of_{at}")[ids]) root.createDimension(at, number_of_elements) @@ -95,9 +110,16 @@ def _set_grid_coordinates( if at == "grid": return - coords = getattr(grid, f"xy_of_{at}") - root.variables[f"x_of_{at}"][:] = coords[ids, 0] - root.variables[f"y_of_{at}"][:] = coords[ids, 1] + if at == "row": + coords = getattr(grid, f"y_of_{at}") + root.variables[f"y_of_{at}"][:] = coords[ids] + elif at == "column": + coords = getattr(grid, f"x_of_{at}") + root.variables[f"x_of_{at}"][:] = coords[ids] + else: + coords = getattr(grid, f"xy_of_{at}") + root.variables[f"x_of_{at}"][:] = coords[ids, 0] + root.variables[f"y_of_{at}"][:] = coords[ids, 1] def _create_field( @@ -139,6 +161,12 @@ def _set_field( ) names = names & set(grid[at]) + if at == "grid": + with contextlib.suppress(KeyError): + names.remove("x_of_shelf_edge") + with contextlib.suppress(KeyError): + names.remove("x_of_shore") + _create_field(root, grid, at=at, names=names) if "time" in root.dimensions: @@ -149,7 +177,6 @@ def _set_field( else: values = grid[at][name] root.variables[_netcdf_var_name(name, at)][n_times - 1, :] = values - else: for name in names: root.variables[_netcdf_var_name(name, at)][:] = grid[at][name][ids] @@ -165,16 +192,21 @@ def _create_layers( netcdf_name = _netcdf_var_name(name, "layer") if netcdf_name not in root.variables: root.createVariable( - netcdf_name, _netcdf_type(grid.event_layers[name]), ("layer", "cell") + netcdf_name, + _netcdf_type(grid.event_layers[name]), + ("layer", "row", "column"), ) netcdf_name = _netcdf_var_name("thickness", "layer") if netcdf_name not in root.variables: - root.createVariable(netcdf_name, "f8", ("layer", "cell")) + root.createVariable(netcdf_name, "f8", ("layer", "row", "column")) def _set_layers( - root: Any, grid: SequenceModelGrid, names: Optional[Iterable[str]] = None + root: Any, + grid: SequenceModelGrid, + names: Optional[Iterable[str]] = None, + ids: Optional[Union[slice, Iterable[int]]] = None, ) -> None: """Set values for variables at a grid layers.""" if isinstance(names, str): @@ -188,12 +220,15 @@ def _set_layers( n_layers = grid.event_layers.number_of_layers for name in names: - root.variables[_netcdf_var_name(name, "layer")][ - :n_layers, : - ] = grid.event_layers[name][:] - root.variables[_netcdf_var_name("thickness", "layer")][ - :n_layers, : - ] = grid.event_layers.dz[:] + layers = grid.event_layers[name][:, ids].reshape( + (-1, grid.shape[0] - 2, grid.shape[1] - 2) + ) + root.variables[_netcdf_var_name(name, "layer")][:n_layers, :, :] = layers + + dz = grid.event_layers.dz[:, ids].reshape( + (-1, grid.shape[0] - 2, grid.shape[1] - 2) + ) + root.variables[_netcdf_var_name("thickness", "layer")][:n_layers, :, :] = dz def to_netcdf( @@ -203,7 +238,7 @@ def to_netcdf( format: str = "NETCDF4", time: float = 0.0, at: Optional[Union[str, Sequence[str]]] = None, - ids: Optional[Union[dict[str, Iterable[int]], int, Iterable[int], slice]] = None, + ids: Optional[dict[str, NDArray[np.integer]] | int | Iterable[int] | slice] = None, names: Optional[ Union[dict[str, Optional[Iterable[str]]], str, Iterable[str]] ] = None, @@ -267,6 +302,12 @@ def to_netcdf( else: names_dict[loc] = list(names_) + with contextlib.suppress(ValueError): + names_dict["grid"].remove("x_of_shore") + with contextlib.suppress(ValueError): + names_dict["grid"].remove("x_of_shelf_edge") + names_dict["row"] = ["x_of_shore", "x_of_shelf_edge"] + ids_dict: dict[str, Union[slice, Iterable[int]]] = defaultdict(lambda: slice(None)) if not isinstance(ids, dict): for loc in at: @@ -286,6 +327,8 @@ def to_netcdf( root.createDimension("time", None) root.createVariable("time", "f8", ("time",)) + _set_grid_coordinates(root, grid, at="row", ids=ids_dict["row"]) + _set_grid_coordinates(root, grid, at="column", ids=ids_dict["column"]) for loc in at: _set_grid_coordinates(root, grid, at=loc, ids=ids_dict[loc]) @@ -294,7 +337,15 @@ def to_netcdf( for loc in at: _set_field(root, grid, at=loc, ids=ids_dict[loc], names=names_dict[loc]) + _set_field( + root, + grid, + at="row", + ids=ids_dict["row"], + names=["x_of_shore", "x_of_shelf_edge"], + ) + if with_layers: - _set_layers(root, grid, names=None) + _set_layers(root, grid, names=None, ids=ids_dict.get("cell", None)) root.close() diff --git a/sequence/output_writer.py b/sequence/output_writer.py index 108dd98..9d816f9 100644 --- a/sequence/output_writer.py +++ b/sequence/output_writer.py @@ -4,6 +4,7 @@ from os import PathLike from typing import Iterable, Optional, Union +import numpy as np from landlab import Component from ._grid import SequenceModelGrid @@ -51,9 +52,9 @@ def __init__( self.filepath = filepath if rows is not None: - self._nodes = grid.nodes[(rows,)].flatten() + self._rows = np.asarray(rows) - 1 else: - self._nodes = None + self._rows = np.arange(grid.shape[0] - 2) self._time = 0.0 self._step_count = 0 @@ -74,7 +75,10 @@ def run_one_step(self, dt: Optional[float] = None) -> None: mode="a", time=self._time, names={"node": self.fields}, - ids={"node": self._nodes}, + ids={ + "row": self._rows, + "column": np.arange(self.grid.shape[1] - 2), + }, ) self._time += dt self._step_count += 1 diff --git a/sequence/plot.py b/sequence/plot.py index dcd2f59..4f78fd0 100644 --- a/sequence/plot.py +++ b/sequence/plot.py @@ -11,7 +11,7 @@ from scipy.interpolate import interp1d from ._grid import SequenceModelGrid -from .errors import MissingRequiredVariable +from .errors import InvalidRowError, MissingRequiredVariable def plot_layers( @@ -171,32 +171,42 @@ def plot_layers( plt.show() -def plot_grid(grid: SequenceModelGrid, **kwds: Any) -> None: +def plot_grid(grid: SequenceModelGrid, row: Optional[int] = None, **kwds: Any) -> None: """Plot a :class:`~SequenceModelGrid`. Parameters ---------- grid : SequenceModelGrid The grid to plot. + row : int, optional + The row of the grid to plot. If not provided, plot the middle row. See Also -------- :func:`plot_layers` : Plot layers from a 2D array of elevations. :func:`plot_file` : Plot a `SequenceModelGrid`'s layers from a file. """ + if row is None: + row = (grid.number_of_rows - 2) // 2 + elevation_at_layer = ( - grid.at_node["bedrock_surface__elevation"][grid.node_at_cell] + grid.at_layer.z + grid.get_profile("bedrock_surface__elevation")[row, 1:-1] + grid.at_layer.z ) - x_of_stack = grid.x_of_node[grid.node_at_cell] - x_of_shore = grid.at_layer_grid["x_of_shore"].flatten() - x_of_shelf_edge = grid.at_layer_grid["x_of_shelf_edge"].flatten() + x_of_stack = grid.x_of_column[1:-1] + + x_of_shore = grid.at_layer_row["x_of_shore"] + x_of_shelf_edge = grid.at_layer_row["x_of_shelf_edge"] time_at_layer_grid = grid.at_layer_grid["age"].flatten() time_at_layer = grid.at_layer["age"] - x_of_shore = interp1d(time_at_layer_grid, x_of_shore)(time_at_layer[:, 0]) - x_of_shelf_edge = interp1d(time_at_layer_grid, x_of_shelf_edge)(time_at_layer[:, 0]) + x_of_shore = interp1d(time_at_layer_grid, x_of_shore[:, row].squeeze())( + time_at_layer[:, 0] + ) + x_of_shelf_edge = interp1d(time_at_layer_grid, x_of_shelf_edge[:, row].squeeze())( + time_at_layer[:, 0] + ) kwds.setdefault("title", f"time = {time_at_layer[-1, 0]} years") @@ -209,13 +219,17 @@ def plot_grid(grid: SequenceModelGrid, **kwds: Any) -> None: ) -def plot_file(filename: Union[str, PathLike], **kwds: Any) -> None: +def plot_file( + filename: Union[str, PathLike], row: Optional[int] = None, **kwds: Any +) -> None: """Plot a `SequenceModelGrid` from a *Sequence* output file. Parameters ---------- filename : path-like Path to the file to plot. + row : int, optional + Row to plot. If not provided, plot the middle row. See Also -------- @@ -224,25 +238,45 @@ def plot_file(filename: Union[str, PathLike], **kwds: Any) -> None: """ kwds.setdefault("title", f"{filename}") with xr.open_dataset(filename) as ds: + if "row" not in ds.dims: + raise MissingRequiredVariable("row") + + if row is None: + row = ds.dims["row"] // 2 + elif (row >= ds.dims["row"]) or (row < -ds.dims["row"]): + raise InvalidRowError(row, ds.dims["row"]) + try: - thickness_at_layer = ds["at_layer:thickness"] - x_of_shore = ds["at_grid:x_of_shore"].data.squeeze() - x_of_shelf_edge = ds["at_grid:x_of_shelf_edge"].data.squeeze() - bedrock = ds["at_node:bedrock_surface__elevation"].data.squeeze() + thickness_at_layer = ds["at_layer:thickness"][:, row, :] + x_of_shore = ds["at_row:x_of_shore"].data + x_of_shelf_edge = ds["at_row:x_of_shelf_edge"].data + bedrock = ( + ds["at_node:bedrock_surface__elevation"] + .data.reshape((-1, ds.dims["row"] + 2, ds.dims["column"] + 2)) + .squeeze() + ) time = ds["time"] - time_at_layer = ds["at_layer:age"] + time_at_layer = ds["at_layer:age"].data.reshape( + (-1, ds.dims["row"], ds.dims["column"]) + ) except KeyError as err: raise MissingRequiredVariable(str(err)) from err try: - x_of_stack = ds["x_of_cell"].data.squeeze() + x_of_stack = ( + ds["x_of_cell"] + .data.reshape((ds.dims["row"], ds.dims["column"]))[row, :] + .squeeze() + ) except KeyError: x_of_stack = np.arange(ds.dims["cell"]) - elevation_at_layer = bedrock[-1, 1:-1] + np.cumsum(thickness_at_layer, axis=0) + elevation_at_layer = bedrock[-1, row, 1:-1] + np.cumsum( + thickness_at_layer, axis=0 + ) - x_of_shore = interp1d(time, x_of_shore)(time_at_layer[:, 0]) - x_of_shelf_edge = interp1d(time, x_of_shelf_edge)(time_at_layer[:, 0]) + x_of_shore = interp1d(time, x_of_shore[:, row])(time_at_layer[:, row, 0]) + x_of_shelf_edge = interp1d(time, x_of_shelf_edge[:, row])(time_at_layer[:, row, 0]) plot_layers( elevation_at_layer, diff --git a/sequence/sediment_flexure.py b/sequence/sediment_flexure.py index 3dacc29..853fc58 100644 --- a/sequence/sediment_flexure.py +++ b/sequence/sediment_flexure.py @@ -30,9 +30,11 @@ def __init__(self, grid, isostasytime: Optional[float] = 7000.0, **kwds: dict): else: self._isostasy_time = self.validate_isostasy_time(isostasytime) - super().__init__(grid, rows=1, method="flexure", **kwds) + active_rows = np.arange(1, grid.shape[0] - 1) - self._subsidence_pool = np.zeros(grid.shape[1], dtype=float) + super().__init__(grid, rows=active_rows, method="flexure", **kwds) + + self._subsidence_pool = np.zeros((len(active_rows), grid.shape[1]), dtype=float) @staticmethod def validate_isostasy_time(time: float) -> float: @@ -103,6 +105,8 @@ def calc_dynamic_deflection( """ isostasy_fraction = self.calc_isostasy_fraction(self.isostasy_time, dt) + isostasy_fraction /= self.grid.shape[0] - 2 + self._subsidence_pool[:] += isostatic_deflection deflection = self._subsidence_pool[:] * isostasy_fraction self._subsidence_pool[:] -= deflection diff --git a/sequence/sequence.py b/sequence/sequence.py index 38825b6..5a9bd4b 100644 --- a/sequence/sequence.py +++ b/sequence/sequence.py @@ -78,6 +78,7 @@ def __init__( self._n_archived_layers = 0 self.grid.at_layer_grid = EventLayers(1) + self.grid.at_layer_row = EventLayers(grid.number_of_rows) if "bedrock_surface__elevation" not in self.grid.at_node: self.grid.add_field( @@ -225,6 +226,14 @@ def add_layer(self, dz_at_cell: NDArray[np.floating]) -> None: x_of_shelf_edge = self.grid.at_grid["x_of_shelf_edge"] except KeyError: x_of_shelf_edge = np.nan + try: + x_of_shores = self.grid.at_row["x_of_shore"] + except KeyError: + x_of_shores = np.full(self.grid.number_of_rows, np.nan) + try: + x_of_shelf_edges = self.grid.at_row["x_of_shelf_edge"] + except KeyError: + x_of_shelf_edges = np.full(self.grid.number_of_rows, np.nan) self.grid.at_node["topographic__elevation"][ self.grid.node_at_cell @@ -238,6 +247,13 @@ def add_layer(self, dz_at_cell: NDArray[np.floating]) -> None: x_of_shore=x_of_shore, x_of_shelf_edge=x_of_shelf_edge, ) + self.grid.at_layer_row.add( + 1.0, + age=self.time, + sea_level=self.grid.at_grid["sea_level__elevation"], + x_of_shore=x_of_shores, + x_of_shelf_edge=x_of_shelf_edges, + ) if ( self.grid.event_layers.number_of_layers - self._n_archived_layers diff --git a/sequence/sequence_model.py b/sequence/sequence_model.py index 68d389e..fc23b25 100644 --- a/sequence/sequence_model.py +++ b/sequence/sequence_model.py @@ -39,13 +39,12 @@ class SequenceModel: "flexure", ) DEFAULT_PARAMS = { - "grid": {"n_cols": 100, "spacing": 1000.0}, + "grid": {"shape": (1, 100), "spacing": (1.0, 1000.0)}, "clock": {"start": 0.0, "stop": 600000.0, "step": 100.0}, "output": { "interval": 10, "filepath": "sequence.nc", "clobber": True, - "rows": [1], "fields": ["sediment_deposit__thickness", "bedrock_surface__elevation"], }, "submarine_diffusion": { @@ -116,8 +115,8 @@ def __init__( if output is not None: self._components["output"] = OutputWriter(self._grid, **output) - self.grid.at_grid["x_of_shore"] = np.nan - self.grid.at_grid["x_of_shelf_edge"] = np.nan + self.grid.at_row["x_of_shore"][:] = np.nan + self.grid.at_row["x_of_shelf_edge"][:] = np.nan self.grid.at_grid["sea_level__elevation"] = 0.0 self._n_archived_layers = 0 diff --git a/sequence/shoreline.py b/sequence/shoreline.py index fa59340..6bf4725 100644 --- a/sequence/shoreline.py +++ b/sequence/shoreline.py @@ -44,7 +44,7 @@ class ShorelineFinder(Component): "intent": "out", "optional": False, "units": "m", - "mapping": "grid", + "mapping": "row", "doc": "Position of the shore line", }, "x_of_shelf_edge": { @@ -52,7 +52,7 @@ class ShorelineFinder(Component): "intent": "out", "optional": False, "units": "m", - "mapping": "grid", + "mapping": "row", "doc": "Position of the shelf edge", }, } @@ -71,24 +71,31 @@ def __init__(self, grid: SequenceModelGrid, alpha: float = 0.0005): super().__init__(grid) self._alpha = alpha - self.grid.at_grid["x_of_shore"] = 0.0 - self.grid.at_grid["x_of_shelf_edge"] = 0.0 + self.grid.at_row["x_of_shore"] = np.zeros(self.grid.shape[0] - 2) + self.grid.at_row["x_of_shelf_edge"] = np.zeros(self.grid.shape[0] - 2) - x = self.grid.x_of_node[self.grid.node_at_cell] - z = self.grid.at_node["topographic__elevation"][self.grid.node_at_cell] - dz = self.grid.at_node["sediment_deposit__thickness"][self.grid.node_at_cell] + x = self.grid.x_of_column + z = self.grid.get_profile("topographic__elevation") + dz = self.grid.get_profile("sediment_deposit__thickness") sea_level = self.grid.at_grid["sea_level__elevation"] - x_of_shore = find_shoreline(x, z, sea_level=sea_level) - try: - x_of_shelf_edge = find_shelf_edge( - x, dz, x_of_shore=x_of_shore, alpha=self._alpha - ) - except ShelfEdgeError: - x_of_shelf_edge = np.nan + x_of_shores = np.empty(z.shape[0], dtype=float) + x_of_shelf_edges = np.empty(z.shape[0], dtype=float) + + for row in range(z.shape[0]): + x_of_shore = find_shoreline(x, z[row], sea_level=sea_level) + try: + x_of_shelf_edge = find_shelf_edge( + x, dz[row], x_of_shore=x_of_shore, alpha=self._alpha + ) + except ShelfEdgeError: + x_of_shelf_edge = np.nan + + x_of_shores[row] = x_of_shore + x_of_shelf_edges[row] = x_of_shelf_edge - self.grid.at_grid["x_of_shore"] = x_of_shore - self.grid.at_grid["x_of_shelf_edge"] = x_of_shelf_edge + self.grid.at_row["x_of_shore"][:] = x_of_shores + self.grid.at_row["x_of_shelf_edge"][:] = x_of_shelf_edges @property def alpha(self) -> float: @@ -97,23 +104,30 @@ def alpha(self) -> float: def update(self) -> None: """Update the component one time step to find the new shoreline.""" - x = self.grid.x_of_node[self.grid.node_at_cell] - z = self.grid.at_node["topographic__elevation"][self.grid.node_at_cell] - dz = self.grid.at_node["sediment_deposit__thickness"][self.grid.node_at_cell] + x = self.grid.x_of_column + z = self.grid.get_profile("topographic__elevation") + dz = self.grid.get_profile("sediment_deposit__thickness") sea_level = self.grid.at_grid["sea_level__elevation"] - x_of_shore = find_shoreline(x, z, sea_level=sea_level) - try: - x_of_shelf_edge = find_shelf_edge( - x, dz, x_of_shore=x_of_shore, alpha=self._alpha - ) - except ShelfEdgeError: - x_of_shelf_edge = np.nan - else: - assert x_of_shelf_edge > x_of_shore + x_of_shores = np.empty(z.shape[0], dtype=float) + x_of_shelf_edges = np.empty(z.shape[0], dtype=float) - self.grid.at_grid["x_of_shore"] = x_of_shore - self.grid.at_grid["x_of_shelf_edge"] = x_of_shelf_edge + for row in range(z.shape[0]): + x_of_shore = find_shoreline(x, z[row], sea_level=sea_level) + try: + x_of_shelf_edge = find_shelf_edge( + x, dz[row], x_of_shore=x_of_shore, alpha=self._alpha + ) + except ShelfEdgeError: + x_of_shelf_edge = np.nan + else: + assert x_of_shelf_edge > x_of_shore + + x_of_shores[row] = x_of_shore + x_of_shelf_edges[row] = x_of_shelf_edge + + self.grid.at_row["x_of_shore"][:] = x_of_shores + self.grid.at_row["x_of_shelf_edge"][:] = x_of_shelf_edges def run_one_step(self, dt: Optional[float] = None) -> None: """Update the component on time step. @@ -163,11 +177,10 @@ def find_shelf_edge_by_curvature( return x[np.argmin(curvature)] -# def find_shelf_edge(grid, x, wd, x_of_shore, sea_level=0.0, alpha = 0.0005): def find_shelf_edge( - x: list[float], - dz: list[float], - x_of_shore: float = 0.0, + x: NDArray[np.floating], + dz: NDArray[np.floating], + x_of_shore: NDArray[np.floating] | float = 0.0, alpha: float = 0.0005, ) -> float: """Find the shelf edge based on deposit thickness. @@ -218,7 +231,7 @@ def find_shoreline( z: NDArray[np.floating], sea_level: float = 0.0, kind: str = "cubic", -) -> float: +) -> NDArray[np.floating] | float: """Find the shoreline of a profile. Parameters @@ -269,18 +282,27 @@ def find_shoreline( x, z = np.asarray(x), np.asarray(z) - try: - index_at_shore = find_shoreline_index(x, z, sea_level=sea_level) - except ShorelineError: - if z[0] < sea_level: - x_of_shoreline = x[0] + z = np.atleast_2d(z) + n_rows = z.shape[0] + + x_of_shorelines = np.empty(n_rows, dtype=float) + for row in range(n_rows): + try: + index_at_shore = find_shoreline_index(x, z[row], sea_level=sea_level) + except ShorelineError: + if z[row, 0] < sea_level: + x_of_shoreline = x[0] + else: + x_of_shoreline = x[-1] else: - x_of_shoreline = x[-1] - else: - func = interpolate.interp1d(x, z - sea_level, kind=kind) - x_of_shoreline = bisect(func, x[index_at_shore - 1], x[index_at_shore]) + func = interpolate.interp1d(x, z[row] - sea_level, kind=kind) + x_of_shoreline = bisect(func, x[index_at_shore - 1], x[index_at_shore]) + x_of_shorelines[row] = x_of_shoreline - return x_of_shoreline + if n_rows == 1: + return x_of_shorelines.item() + else: + return x_of_shorelines def _find_shoreline_polyfit( @@ -330,7 +352,7 @@ def _find_shoreline_polyfit( def find_shoreline_index( x: NDArray[np.floating], z: NDArray[np.floating], sea_level: float = 0.0 -) -> int: +) -> NDArray[np.integer] | int: """Find the landward-index of the shoreline. Parameters @@ -381,14 +403,25 @@ def find_shoreline_index( Traceback (most recent call last): sequence.errors.ShorelineError: No shoreline found. The profile is all above sea level """ - (below_water,) = np.where(z < sea_level) - - if len(below_water) == 0 or len(below_water) == len(x): - raise ShorelineError( - "No shoreline found. The profile is all {} sea level".format( - "above" if len(below_water) == 0 else "below" + z = np.atleast_2d(z) + n_rows = z.shape[0] + + index_of_shoreline = np.full(n_rows, -1) + for row in range(n_rows): + (below_water,) = np.where(z[row] < sea_level) + + if len(below_water) == 0 or len(below_water) == len(x): + raise ShorelineError( + "No shoreline found. The profile is all {} sea level".format( + "above" if len(below_water) == 0 else "below" + ) ) - ) - return None + # return None + else: + index_of_shoreline[row] = below_water[0] + # return below_water[0] + + if n_rows == 1: + return index_of_shoreline.item() else: - return below_water[0] + return index_of_shoreline diff --git a/sequence/submarine.py b/sequence/submarine.py index 4284b38..fff093f 100644 --- a/sequence/submarine.py +++ b/sequence/submarine.py @@ -168,10 +168,6 @@ def sediment_load(self, value: float) -> None: self._ksh = self._load / self._plain_slope self.grid.at_grid["sediment_load"] = self._load - # @property - # def k0(self): - # return self._k0 - @property def load0(self) -> float: """Return the sediment load entering the profile.""" @@ -196,7 +192,9 @@ def sea_level(self) -> float: def sea_level(self, sea_level: float) -> None: self.grid.at_grid["sea_level__elevation"] = sea_level - def calc_diffusion_coef(self, x_of_shore: float) -> NDArray[np.floating]: + def calc_diffusion_coef( + self, x_of_shore: NDArray[np.floating] | float + ) -> NDArray[np.floating]: """Calculate and store diffusion coefficient values. Parameters @@ -210,10 +208,10 @@ def calc_diffusion_coef(self, x_of_shore: float) -> NDArray[np.floating]: sea level and three below, two of which are in deep water (below the default 60 m wave base). - >>> from landlab import RasterModelGrid + >>> from sequence import SequenceModelGrid >>> import numpy as np - >>> grid = RasterModelGrid((3, 6), xy_spacing=200.0) + >>> grid = SequenceModelGrid((1, 6), spacing=(1.0, 200.0)) >>> z = grid.add_zeros("topographic__elevation", at="node") >>> z[6:12] = np.array([3., 3., 1., -1., -85., -85.]) >>> z.reshape((3, 6)) @@ -232,36 +230,74 @@ def calc_diffusion_coef(self, x_of_shore: float) -> NDArray[np.floating]: >>> diffusion_coef is grid.at_node["kd"] True """ - sea_level = self.grid.at_grid["sea_level__elevation"] - water_depth = sea_level - self._grid.at_node["topographic__elevation"] - - under_water = water_depth > 0.0 - deep_water = water_depth > self._wave_base - land = ~under_water - - k = self.grid.at_node["kd"] - - x = self.grid.x_of_node - b = (self._shoreface_height * self._alpha + self._shelf_slope) * self.grid.dx - - k[under_water] = ( - self._load - * ((x[under_water] - x_of_shore) + self.grid.dx) - / (water_depth[under_water] + b) - ) - - k[deep_water] *= np.exp( - -(water_depth[deep_water] - self._wave_base) / self._wave_base - ) + x_of_shore = np.atleast_1d(x_of_shore) - self._load = self._load0 * (1 + sea_level * self._load_sl) - self._ksh = self._load / self._plain_slope - if self._basin_width > 0.0: - k[land] = self._ksh * (self._basin_width + x[land]) / self._basin_width - else: - k[land] = self._ksh + sea_level = self.grid.at_grid["sea_level__elevation"] - return k + water_depth = sea_level - self._grid.get_profile("topographic__elevation") + k = self.grid.get_profile("kd") + x = self.grid.x_of_column + + n_rows = k.shape[0] + + assert len(x_of_shore) == n_rows + + for row in range(n_rows): + under_water = water_depth[row] > 0.0 + deep_water = water_depth[row] > self._wave_base + land = ~under_water + + b = ( + self._shoreface_height * self._alpha + self._shelf_slope + ) * self.grid.dx + + k[row, under_water] = ( + self._load + * ((x[under_water] - x_of_shore[row]) + self.grid.dx) + / (water_depth[row, under_water] + b) + ) + + k[row, deep_water] *= np.exp( + -(water_depth[row, deep_water] - self._wave_base) / self._wave_base + ) + + self._load = self._load0 * (1 + sea_level * self._load_sl) + self._ksh = self._load / self._plain_slope + + if self._basin_width > 0.0: + k[row, land] = ( + self._ksh * (self._basin_width + x[land]) / self._basin_width + ) + else: + k[row, land] = self._ksh + + # TODO: modify diffusion outside of the channel row. + if row != n_rows // 2: + k[row, land] *= 0.5 + + # if row == n_rows // 2: + # if self._basin_width > 0.0: + # k[row, land] = ( + # self._ksh * (self._basin_width + x[land]) / self._basin_width + # ) + # else: + # k[row, land] = self._ksh + # else: # outside of the channel with low diffusivity + # if self._basin_width > 0.0: + # k[row, land] = self._ksh * (self.grid.dx + x[land]) / self._basin_width + # else: + # k[row, land] = self._ksh * (self.grid.dx) / self._basin_width + + # if self._basin_width > 0.0: + # k[land] = self._ksh * (self._basin_width + x[land]) / self._basin_width + # else: + # k[land] = self._ksh + + k = self.grid.at_node["kd"].reshape(self.grid.shape) + k[0, :] = k[1, :] + k[-1, :] = k[-1, :] + + return self.grid.at_node["kd"] def run_one_step(self, dt: float) -> None: """Advance component one time step. @@ -272,20 +308,22 @@ def run_one_step(self, dt: float) -> None: Time step to advance component by. """ shore = find_shoreline( - self.grid.x_of_node[self.grid.node_at_cell], - self.grid.at_node["topographic__elevation"][self.grid.node_at_cell], + self.grid.x_of_column, + self.grid.get_profile("topographic__elevation"), sea_level=self.grid.at_grid["sea_level__elevation"], ) self.calc_diffusion_coef(shore) # set elevation at upstream boundary to ensure proper sediment influx - x = self.grid.x_of_node.reshape(self.grid.shape) + x = self.grid.x_of_column z = self._grid.at_node["topographic__elevation"].reshape(self.grid.shape) # k = self._grid.at_node["kd"].reshape(self.grid.shape) # z[1, 0] = z[1,1] + self._load / k[1, 0] * (x[1,1]-x[1,0]) - z[1, 0] = z[1, 1] + self._plain_slope * (x[1, 1] - x[1, 0]) - # self._load/self._load0) + # z[1, 0] = z[1, 1] + self._plain_slope * (x[1, 1] - x[1, 0]) + z[1, 0] = z[1, 1] + self._plain_slope * (x[1] - x[0]) + + z[1:-1, 0] = z[1:-1, 1] + self._plain_slope * (x[1] - x[0]) z_before = self.grid.at_node["topographic__elevation"].copy() diff --git a/tests/test_grid.py b/tests/test_grid.py index 44edb85..a84b309 100644 --- a/tests/test_grid.py +++ b/tests/test_grid.py @@ -22,13 +22,11 @@ def test_from_dict(): def test_from_dict_old_style(): expected = SequenceModelGrid(500, spacing=5.0) - actual = SequenceModelGrid.from_dict( - {"shape": (13, 500), "xy_spacing": (5.0, 25.0)} - ) + actual = SequenceModelGrid.from_dict({"shape": (1, 500), "xy_spacing": (5.0, 1.0)}) assert actual.shape == expected.shape assert actual.spacing == expected.spacing - actual = SequenceModelGrid.from_dict({"shape": (13, 500), "xy_spacing": 5.0}) + actual = SequenceModelGrid.from_dict({"shape": 500, "xy_spacing": 5.0}) assert actual.shape == expected.shape assert actual.spacing == expected.spacing diff --git a/tests/test_netcdf.py b/tests/test_netcdf.py index 7c75f24..12c4e27 100644 --- a/tests/test_netcdf.py +++ b/tests/test_netcdf.py @@ -94,7 +94,7 @@ def test_with_layers(tmpdir): grid = SequenceModelGrid(4) grid.event_layers.add(10.0, age=0.0, water_depth=np.arange(2)) with tmpdir.as_cwd(): - to_netcdf(grid, "test.nc") + to_netcdf(grid, "test.nc", ids={"row": [1], "column": [1, 2]}) ds = xr.open_dataset("test.nc") assert ds.dims["layer"] == 1 assert np.all(ds["at_layer:thickness"] == np.array([[10.0, 10.0]])) diff --git a/tests/test_output_writer.py b/tests/test_output_writer.py index dac021b..6869eed 100644 --- a/tests/test_output_writer.py +++ b/tests/test_output_writer.py @@ -1,13 +1,13 @@ import numpy as np import xarray as xr -from landlab import RasterModelGrid from pytest import approx, raises +from sequence import SequenceModelGrid from sequence.output_writer import OutputWriter def test_no_fields(tmpdir): - grid = RasterModelGrid((3, 4)) + grid = SequenceModelGrid((1, 4), spacing=(1.0, 1.0)) with tmpdir.as_cwd(): writer = OutputWriter(grid, "test.nc") writer.run_one_step() @@ -19,7 +19,7 @@ def test_no_fields(tmpdir): def test_default_output_interval(tmpdir): - grid = RasterModelGrid((3, 4)) + grid = SequenceModelGrid((1, 4), spacing=(1.0, 1.0)) with tmpdir.as_cwd(): writer = OutputWriter(grid, "test.nc") for _ in range(4): @@ -32,7 +32,7 @@ def test_default_output_interval(tmpdir): def test_output_interval(tmpdir): - grid = RasterModelGrid((3, 4)) + grid = SequenceModelGrid((1, 4), spacing=(1.0, 1.0)) with tmpdir.as_cwd(): writer = OutputWriter(grid, "test.nc", interval=2) for _ in range(3): @@ -45,7 +45,7 @@ def test_output_interval(tmpdir): def test_default_clobber_is_false(tmpdir): - grid = RasterModelGrid((3, 4)) + grid = SequenceModelGrid((1, 4), spacing=(1.0, 1.0)) with tmpdir.as_cwd(): OutputWriter(grid, "test.nc").run_one_step() with raises(RuntimeError): @@ -53,7 +53,7 @@ def test_default_clobber_is_false(tmpdir): def test_clobber_keyword(tmpdir): - grid = RasterModelGrid((3, 4)) + grid = SequenceModelGrid((1, 4), spacing=(1.0, 1.0)) with tmpdir.as_cwd(): writer = OutputWriter(grid, "test.nc") for _ in range(10): @@ -70,7 +70,7 @@ def test_clobber_keyword(tmpdir): def test_default_is_all_rows(tmpdir): - grid = RasterModelGrid((3, 4)) + grid = SequenceModelGrid((1, 4), spacing=(1.0, 1.0)) with tmpdir.as_cwd(): OutputWriter(grid, "test.nc").run_one_step() with xr.open_dataset("test.nc") as ds: @@ -78,14 +78,15 @@ def test_default_is_all_rows(tmpdir): def test_rows_keyword(tmpdir): - grid = RasterModelGrid((3, 4)) + grid = SequenceModelGrid((1, 4), spacing=(1.0, 1.0)) with tmpdir.as_cwd(): OutputWriter(grid, "test.nc", rows=(1,)).run_one_step() with xr.open_dataset("test.nc") as ds: - assert np.all(ds.x_of_node == approx([0, 1, 2, 3])) - assert np.all(ds.y_of_node == approx(1.0)) + assert np.all(ds.x_of_column == approx([0, 1])) + assert np.all(ds.y_of_row == approx(0.0)) - OutputWriter(grid, "test.nc", rows=((0, 2),), clobber=True).run_one_step() + # OutputWriter(grid, "test.nc", rows=((0, 2),), clobber=True).run_one_step() + OutputWriter(grid, "test.nc", rows=(0, 2), clobber=True).run_one_step() with xr.open_dataset("test.nc") as ds: - assert np.all(ds.x_of_node == approx([0, 1, 2, 3, 0, 1, 2, 3])) - assert np.all(ds.y_of_node == approx([0, 0, 0, 0, 2, 2, 2, 2])) + assert np.all(ds.x_of_column == approx([0, 1])) + assert np.all(ds.y_of_row == approx([2, 1])) diff --git a/tests/test_sequence.py b/tests/test_sequence.py index b04b5fc..9ecdbe1 100644 --- a/tests/test_sequence.py +++ b/tests/test_sequence.py @@ -70,15 +70,15 @@ def test_sequence_x_of_shore_missing(grid): def test_sequence_x_of_shore(grid): seq = Sequence(grid, components=[ShorelineFinder(grid)]) - assert "x_of_shore" in grid.at_grid - assert "x_of_shelf_edge" in grid.at_grid + assert "x_of_shore" in grid.at_row + assert "x_of_shelf_edge" in grid.at_row assert grid.at_layer_grid.number_of_layers == 1 - assert grid.at_layer_grid["x_of_shore"][-1] == pytest.approx(20000.0) - assert grid.at_layer_grid["x_of_shelf_edge"][-1] > 20000.0 + assert np.all(grid.at_layer_row["x_of_shore"][-1, :] == pytest.approx(20000.0)) + assert np.all(grid.at_layer_row["x_of_shelf_edge"][-1, :] > 20000.0) seq.update() assert grid.at_layer_grid.number_of_layers == 2 - assert grid.at_layer_grid["x_of_shore"][-1] == pytest.approx(20000.0) - assert grid.at_layer_grid["x_of_shelf_edge"][-1] > 20000.0 + assert np.all(grid.at_layer_row["x_of_shore"][-1, :] == pytest.approx(20000.0)) + assert np.all(grid.at_layer_row["x_of_shelf_edge"][-1, :] > 20000.0) diff --git a/tests/test_sequence_model/marmara.yaml b/tests/test_sequence_model/marmara.yaml index ae28f82..c3593e1 100644 --- a/tests/test_sequence_model/marmara.yaml +++ b/tests/test_sequence_model/marmara.yaml @@ -15,8 +15,7 @@ output: interval: 10 filepath: marmara.nc clobber: True - rows: - - 1 + rows: [0, 1, 2] fields: - sediment_deposit__thickness - kd diff --git a/tests/test_subsidence.py b/tests/test_subsidence.py index 7c7a9af..ce89d4d 100644 --- a/tests/test_subsidence.py +++ b/tests/test_subsidence.py @@ -73,7 +73,7 @@ def test_linear_subsidence(tmpdir): subsidence.run_one_step(dt=1.0) assert_array_equal( grid.get_profile("bedrock_surface__increment_of_elevation"), - [0.0, 10.0, 20.0, 30.0, 40.0], + [[0.0, 10.0, 20.0, 30.0, 40.0]], ) @@ -91,7 +91,7 @@ def test_add_to_existing_subsidence(tmpdir): subsidence.run_one_step(dt=1.0) assert_array_equal( grid.get_profile("bedrock_surface__increment_of_elevation"), - [1.0, 11.0, 21.0, 31.0, 41.0], + [[1.0, 11.0, 21.0, 31.0, 41.0]], )