diff --git a/requirements.txt b/requirements.txt index a912259e..4ead18e6 100644 --- a/requirements.txt +++ b/requirements.txt @@ -3,3 +3,4 @@ numpy xarray dask landlab>=1.7 +bmipy diff --git a/terrainbento/__init__.py b/terrainbento/__init__.py index dc24c899..b56b31fa 100644 --- a/terrainbento/__init__.py +++ b/terrainbento/__init__.py @@ -9,53 +9,54 @@ have been designed to modify model boundary conditions during a model run. """ +from ._version import get_versions +from .base_class import ( + ErosionModel, + StochasticErosionModel, + TwoLithologyErosionModel, +) +from .boundary_handlers import ( + CaptureNodeBaselevelHandler, + GenericFuncBaselevelHandler, + NotCoreNodeBaselevelHandler, + PrecipChanger, + SingleNodeBaselevelHandler, +) from .clock import Clock - -from .base_class import ErosionModel -from .base_class import StochasticErosionModel -from .base_class import TwoLithologyErosionModel - +from .derived_models import ( + Basic, + BasicCh, + BasicChRt, + BasicChRtTh, + BasicChSa, + BasicCv, + BasicDd, + BasicDdHy, + BasicDdRt, + BasicDdSt, + BasicDdVs, + BasicHy, + BasicHyRt, + BasicHySa, + BasicHySt, + BasicHyVs, + BasicRt, + BasicRtSa, + BasicRtTh, + BasicRtVs, + BasicSa, + BasicSaVs, + BasicSt, + BasicStTh, + BasicStVs, + BasicTh, + BasicThVs, + BasicVs, +) from .model_template import ModelTemplate - -from .boundary_handlers import PrecipChanger -from .boundary_handlers import SingleNodeBaselevelHandler -from .boundary_handlers import CaptureNodeBaselevelHandler -from .boundary_handlers import NotCoreNodeBaselevelHandler -from .boundary_handlers import GenericFuncBaselevelHandler - -from .precipitators import UniformPrecipitator, RandomPrecipitator +from .precipitators import RandomPrecipitator, UniformPrecipitator from .runoff_generators import SimpleRunoff -from .derived_models import Basic -from .derived_models import BasicTh -from .derived_models import BasicDd -from .derived_models import BasicHy -from .derived_models import BasicCh -from .derived_models import BasicSt -from .derived_models import BasicVs -from .derived_models import BasicSa -from .derived_models import BasicRt -from .derived_models import BasicCv -from .derived_models import BasicDdHy -from .derived_models import BasicStTh -from .derived_models import BasicDdSt -from .derived_models import BasicHySt -from .derived_models import BasicThVs -from .derived_models import BasicDdVs -from .derived_models import BasicHyVs -from .derived_models import BasicStVs -from .derived_models import BasicHySa -from .derived_models import BasicChSa -from .derived_models import BasicSaVs -from .derived_models import BasicRtTh -from .derived_models import BasicDdRt -from .derived_models import BasicHyRt -from .derived_models import BasicChRt -from .derived_models import BasicRtVs -from .derived_models import BasicRtSa -from .derived_models import BasicChRtTh - - __all__ = [ "Clock", "ModelTemplate", @@ -100,7 +101,6 @@ "TwoLithologyErosionModel", ] -from ._version import get_versions __version__ = get_versions()["version"] del get_versions diff --git a/terrainbento/_version.py b/terrainbento/_version.py index cef31606..fddb9a1a 100644 --- a/terrainbento/_version.py +++ b/terrainbento/_version.py @@ -1,4 +1,3 @@ - # This file helps to compute a version number in source trees obtained from # git-archive tarball (such as those provided by githubs download-from-tag # feature). Distribution tarballs (built by setup.py sdist) and build diff --git a/terrainbento/base_class/erosion_model.py b/terrainbento/base_class/erosion_model.py index 45ad4572..e6a22920 100644 --- a/terrainbento/base_class/erosion_model.py +++ b/terrainbento/base_class/erosion_model.py @@ -6,6 +6,7 @@ import sys import time as tm +import cfunits import dask import numpy as np import six @@ -16,7 +17,7 @@ from landlab.components import FlowAccumulator, NormalFault from landlab.graph import Graph from landlab.io.netcdf import write_raster_netcdf -from terrainbento import Clock +from terrainbento.bmi import BmiModel from terrainbento.boundary_handlers import ( CaptureNodeBaselevelHandler, GenericFuncBaselevelHandler, @@ -24,6 +25,7 @@ PrecipChanger, SingleNodeBaselevelHandler, ) +from terrainbento.clock import Clock from terrainbento.precipitators import RandomPrecipitator, UniformPrecipitator from terrainbento.runoff_generators import SimpleRunoff @@ -126,7 +128,7 @@ def _setup_boundary_handlers(grid, name, params): return boundary_handler -class ErosionModel(object): +class ErosionModel(BmiModel): """Base class providing common functionality for terrainbento models. @@ -139,13 +141,28 @@ class ErosionModel(object): It is expected that a derived model will define an **__init__** and a **run_one_step** method. If desired, the derived model can overwrite the - existing **run_for**, **run**, and **finalize** methods. + existing **run_for** and **run** methods. The following at-node fields must be specified in the grid: - ``topographic__elevation`` """ - _required_fields = ["topographic__elevation"] + _name = "ErosionModel" + + _input_var_names = ("topographic__elevation",) + + _output_var_names = ("topographic__elevation",) + + _var_info = { + "topographic__elevation": { + "type": float, + "units": "m", + "at": "node", + "description": "Land surface topographic elevation", + } + } + + _param_info = {} @classmethod def from_file(cls, file_like): @@ -160,6 +177,7 @@ def from_file(cls, file_like): Examples -------- >>> from six import StringIO + >>> from terrainbento import ErosionModel >>> filelike = StringIO(''' ... grid: ... grid: @@ -179,18 +197,18 @@ def from_file(cls, file_like): 1.0 >>> model.clock.stop 200.0 - >>> model.grid.shape + >>> model._grid.shape (4, 5) """ # first get contents. try: contents = file_like.read() - except AttributeError: # was a str + except AttributeError: if os.path.isfile(file_like): with open(file_like, "r") as fp: contents = fp.read() else: - contents = file_like # not tested + contents = file_like # then parse contents. params = yaml.safe_load(contents) @@ -204,7 +222,7 @@ def from_dict(cls, params, output_writers=None): The input parameter dictionary portion associated with the "grid" keword will be passed directly to the Landlab - `create_grid `_. + `create_grid `_. function. Parameters @@ -237,7 +255,7 @@ def from_dict(cls, params, output_writers=None): 1.0 >>> model.clock.stop 200.0 - >>> model.grid.shape + >>> model._grid.shape (4, 5) """ cls._validate(params) @@ -370,29 +388,30 @@ def __init__( recommend that you look at the terrainbento tutorials for examples of usage. """ + self.clock = clock + self._grid = grid + # call ErosionModel's base class init + super(ErosionModel, self).__init__() + flow_accumulator_kwargs = flow_accumulator_kwargs or {} boundary_handlers = boundary_handlers or {} output_writers = output_writers or {} fields = fields or ["topographic__elevation"] + # type checking if isinstance(clock, Clock) is False: raise ValueError("Provided Clock is not valid.") if isinstance(grid, ModelGrid) is False: raise ValueError("Provided Grid is not valid.") - # save the grid, clock, and parameters. - self.grid = grid - self.clock = clock - - # first pass of verifying fields - self._verify_fields(self._required_fields) - # save reference to elevation + self._verify_fields(ErosionModel._input_var_names) + self.z = grid.at_node["topographic__elevation"] - self.grid.add_zeros("node", "cumulative_elevation_change") + self._grid.add_zeros("node", "cumulative_elevation_change") - self.grid.add_field( + self._grid.add_field( "node", "initial_topographic__elevation", self.z.copy() ) @@ -405,9 +424,6 @@ def __init__( output_interval = clock.stop self.output_interval = output_interval - # instantiate model time. - self._model_time = 0. - # instantiate container for computational timestep: self._compute_time = [tm.time()] @@ -417,7 +433,7 @@ def __init__( # verify that precipitator is valid if precipitator is None: - precipitator = UniformPrecipitator(self.grid) + precipitator = UniformPrecipitator(self._grid) else: if isinstance(precipitator, _VALID_PRECIPITATORS) is False: raise ValueError("Provided value for precipitator not valid.") @@ -425,7 +441,7 @@ def __init__( # verify that runoff_generator is valid if runoff_generator is None: - runoff_generator = SimpleRunoff(self.grid) + runoff_generator = SimpleRunoff(self._grid) else: if isinstance(runoff_generator, _VALID_RUNOFF_GENERATORS) is False: raise ValueError( @@ -443,18 +459,21 @@ def __init__( flow_director == "FlowDirectorSteepest" ): self.flow_accumulator = FlowAccumulator( - self.grid, + self._grid, routing="D4", depression_finder=depression_finder, **flow_accumulator_kwargs ) else: self.flow_accumulator = FlowAccumulator( - self.grid, + self._grid, depression_finder=depression_finder, **flow_accumulator_kwargs ) + # first pass of verifying fields + self._verify_fields(self._input_var_names) + ################################################################### # Boundary Conditions and Output Writers ################################################################### @@ -469,24 +488,21 @@ def __init__( self.output_writers = output_writers + # todo: check field units + def _verify_fields(self, required_fields): """Verify all required fields are present.""" for field in required_fields: - if field not in self.grid.at_node: + if field not in self._grid.at_node: raise ValueError( "Required field {field} not present.".format(field=field) ) - @property - def model_time(self): - """Return current time of model integration in model time units.""" - return self._model_time - def calculate_cumulative_change(self): """Calculate cumulative node-by-node changes in elevation.""" - self.grid.at_node["cumulative_elevation_change"][:] = ( - self.grid.at_node["topographic__elevation"] - - self.grid.at_node["initial_topographic__elevation"] + self._grid.at_node["cumulative_elevation_change"][:] = ( + self._grid.at_node["topographic__elevation"] + - self._grid.at_node["initial_topographic__elevation"] ) def create_and_move_water(self, step): @@ -512,21 +528,24 @@ def write_output(self): self._output_files.append(filename) try: write_raster_netcdf( - filename, self.grid, names=self.output_fields, format="NETCDF4" + filename, + self._grid, + names=self.output_fields, + format="NETCDF4", ) except NotImplementedError: graph = Graph.from_dict( { - "y_of_node": self.grid.y_of_node, - "x_of_node": self.grid.x_of_node, - "nodes_at_link": self.grid.nodes_at_link, + "y_of_node": self._grid.y_of_node, + "x_of_node": self._grid.x_of_node, + "nodes_at_link": self._grid.nodes_at_link, } ) for field_name in self.output_fields: graph._ds.__setitem__( - field_name, ("node", self.grid.at_node[field_name]) + field_name, ("node", self._grid.at_node[field_name]) ) graph.to_netcdf(path=filename, mode="w", format="NETCDF4") @@ -539,18 +558,14 @@ def finalize__run_one_step(self, step): This base-class method increments model time and updates boundary conditions. """ - # calculate model time - self._model_time += step + # update model time + self.clock.advance(step) # Update boundary conditions self.update_boundary_conditions(step) - def finalize(self): - """Finalize model. - - This base-class method does nothing. Derived classes can - override it to run any required finalization steps. - """ + def finalize_terrainbento_run(self): + """Finalize erosion models.""" pass def run_for(self, step, runtime): @@ -566,7 +581,7 @@ def run_for(self, step, runtime): runtime : float Total duration for which to run model. """ - elapsed_time = 0. + elapsed_time = 0.0 keep_running = True while keep_running: if elapsed_time + step >= runtime: @@ -590,19 +605,18 @@ def run(self): self._itters.append(0) self.write_output() self.iteration = 1 - time_now = self._model_time - while time_now < self.clock.stop: + while self.clock.time < self.clock.stop: next_run_pause = min( - time_now + self.output_interval, self.clock.stop + self.clock.time + self.output_interval, self.clock.stop ) - self.run_for(self.clock.step, next_run_pause - time_now) - time_now = self._model_time + self.run_for(self.clock.step, next_run_pause - self.clock.time) + self.iteration += 1 self._itters.append(self.iteration) self.write_output() # now that the model is finished running, execute finalize. - self.finalize() + self.finalize_terrainbento_run() def _ensure_precip_runoff_are_vanilla(self, vsa_precip=False): """Ensure only default versions of precipitator/runoff are used. diff --git a/terrainbento/base_class/stochastic_erosion_model.py b/terrainbento/base_class/stochastic_erosion_model.py index 17c66726..a5732c7b 100644 --- a/terrainbento/base_class/stochastic_erosion_model.py +++ b/terrainbento/base_class/stochastic_erosion_model.py @@ -16,7 +16,7 @@ class StochasticErosionModel(ErosionModel): - """Base class for stochastic-precipitation terrainbento models. + r"""Base class for stochastic-precipitation terrainbento models. A **StochasticErosionModel** inherits from :py:class:`ErosionModel` and provides functionality needed by all @@ -60,11 +60,11 @@ class StochasticErosionModel(ErosionModel): routing method. It then performs one of two options, depending on the user's choice of ``opt_stochastic_duration`` (True or False). - If the user requests stochastic duration, the model iterates through a sequence - of storm and interstorm periods. Storm depth is drawn at random from a gamma - distribution, and storm duration from an exponential distribution; storm - intensity is then depth divided by duration. This sequencing is implemented by - overriding the run_for method. + If the user requests stochastic duration, the model iterates through a + sequence of storm and interstorm periods. Storm depth is drawn at random + from a gamma distribution, and storm duration from an exponential + distribution; storm intensity is then depth divided by duration. This + sequencing is implemented by overriding the run_for method. If the user does not request stochastic duration (indicated by setting ``opt_stochastic_duration`` to ``False``), then the default @@ -84,18 +84,22 @@ class StochasticErosionModel(ErosionModel): R = P - I (1 - \exp ( -P / I )) - where :math:`I` is the soil infiltration capacity. At the sub-grid scale, soil - infiltration capacity is assumed to have an exponential distribution of which - :math:`I` is the mean. Hence, there are always some spots within any given grid cell - that will generate runoff. This approach yields a smooth transition from - near-zero runoff (when :math:`I>>P`) to :math:`R \\approx P` - (when :math:`P>>I`), without a "hard threshold." + where :math:`I` is the soil infiltration capacity. At the sub-grid scale, + soil infiltration capacity is assumed to have an exponential distribution + of which :math:`I` is the mean. Hence, there are always some spots within + any given grid cell that will generate runoff. This approach yields a + smooth transition from near-zero runoff (when :math:`I>>P`) to + :math:`R \\approx P` (when :math:`P>>I`), without a "hard threshold." The following at-node fields must be specified in the grid: - ``topographic__elevation`` """ - _required_fields = ["topographic__elevation"] + _name = "StochasticErosionModel" + + _input_var_names = ("topographic__elevation",) + + _output_var_names = ("topographic__elevation",) def __init__( self, @@ -201,8 +205,8 @@ def __init__( # initialize record for storms. Depending on how this model is run # (stochastic time, number_time_steps>1, more manually) the step may # change. Thus, rather than writing routines to reconstruct the time - # series of precipitation from the step could change based on users use, - # we"ll record this with the model run instead of re-running. + # series of precipitation from the step could change based on users + # use, we"ll record this with the model run instead of re-running. # make this the non-default option. @@ -229,8 +233,8 @@ def calc_runoff_and_discharge(self): runoff = 0 # pragma: no cover else: runoff = self.rain_rate - self.grid.at_node["surface_water__discharge"][:] = ( - runoff * self.grid.at_node["drainage_area"] + self._grid.at_node["surface_water__discharge"][:] = ( + runoff * self._grid.at_node["drainage_area"] ) return runoff @@ -358,13 +362,13 @@ def handle_water_erosion(self, step, flooded): if self.record_rain: # save record into the rain record self.record_rain_event( - self.model_time, step, self.rain_rate, runoff + self.clock.time, step, self.rain_rate, runoff ) elif self.opt_stochastic_duration and self.rain_rate <= 0.0: # calculate and record the time with no rain: if self.record_rain: - self.record_rain_event(self.model_time, step, 0, 0) + self.record_rain_event(self.clock.time, step, 0, 0) elif not self.opt_stochastic_duration: @@ -382,7 +386,7 @@ def handle_water_erosion(self, step, flooded): self.eroder.run_one_step(dt_water, flooded_nodes=flooded) # save record into the rain record if self.record_rain: - event_start_time = self.model_time + (i * dt_water) + event_start_time = self.clock.time + (i * dt_water) self.record_rain_event( event_start_time, dt_water, self.rain_rate, runoff ) @@ -396,12 +400,12 @@ def handle_water_erosion(self, step, flooded): # if dry time is greater than zero, record. if dt_dry > 0: - event_start_time = self.model_time + ( + event_start_time = self.clock.time + ( self.n_sub_steps * dt_water ) self.record_rain_event(event_start_time, dt_dry, 0.0, 0.0) - def finalize(self): + def finalize_terrainbento_run(self): """Finalize stochastic erosion models. The finalization step of stochastic erosion models in @@ -581,11 +585,11 @@ def write_exceedance_frequency_file( # calculate the predictions for 10, 25, and 100 year event based on # the analytical form of the exceedance function. - event_intervals = np.array([10., 25, 100.]) + event_intervals = np.array([10.0, 25, 100.0]) # calculate the probability of each event based on the number of years # and the number of wet days per year. - daily_distribution_exceedance_probabilities = 1. / ( + daily_distribution_exceedance_probabilities = 1.0 / ( nwet * event_intervals ) @@ -601,8 +605,8 @@ def write_exceedance_frequency_file( # po = P * (- ln (P(p>po))) ^ (1 / c) expected_rainfall = self.scale_factor * ( - -1. * np.log(daily_distribution_exceedance_probabilities) - ) ** (1. / self.shape_factor) + -1.0 * np.log(daily_distribution_exceedance_probabilities) + ) ** (1.0 / self.shape_factor) exceedance_file.write("\n\nSection 2: Theoretical Predictions\n") @@ -706,7 +710,7 @@ def write_exceedance_frequency_file( ) * ( np.exp( - -1. + -1.0 * (expected_rainfall / self.scale_factor) ** self.shape_factor ) @@ -779,7 +783,7 @@ def write_exceedance_frequency_file( ] = selected_wet_day_totals.max() # calculate the distribution percentiles associated with each interval - event_percentiles = (1. - (1. / event_intervals)) * 100. + event_percentiles = (1.0 - (1.0 / event_intervals)) * 100.0 # calculated the event magnitudes associated with the percentiles. event_magnitudes = np.percentile( diff --git a/terrainbento/base_class/two_lithology_erosion_model.py b/terrainbento/base_class/two_lithology_erosion_model.py index 86b583c9..b42fd025 100644 --- a/terrainbento/base_class/two_lithology_erosion_model.py +++ b/terrainbento/base_class/two_lithology_erosion_model.py @@ -26,10 +26,14 @@ class TwoLithologyErosionModel(ErosionModel): - ``lithology_contact__elevation`` """ - _required_fields = [ + _name = "TwoLithologyErosionModel" + + _input_var_names = ( "topographic__elevation", "lithology_contact__elevation", - ] + ) + + _output_var_names = ("topographic__elevation",) def __init__( self, @@ -40,7 +44,7 @@ def __init__( water_erodibility_lower=0.0001, water_erodibility_upper=0.001, regolith_transport_parameter=0.1, - contact_zone__width=1., + contact_zone__width=1.0, **kwargs ): """ @@ -68,7 +72,7 @@ def __init__( super(TwoLithologyErosionModel, self).__init__(clock, grid, **kwargs) # verify correct fields are present. - self._verify_fields(self._required_fields) + self._verify_fields(self._input_var_names) self.m = m_sp self.n = n_sp @@ -86,7 +90,7 @@ def __init__( self.rock_erody = self.K_rock self.till_erody = self.K_till - self.rock_till_contact = self.grid.at_node[ + self.rock_till_contact = self._grid.at_node[ "lithology_contact__elevation" ] @@ -95,10 +99,10 @@ def _setup_rock_and_till(self): erodibility.""" # Create field for erodibility - self.erody = self.grid.add_zeros("node", "substrate__erodibility") + self.erody = self._grid.add_zeros("node", "substrate__erodibility") # Create array for erodibility weighting function - self.erody_wt = np.zeros(self.grid.number_of_nodes) + self.erody_wt = np.zeros(self._grid.number_of_nodes) # Set values correctly self._update_erodywt() @@ -109,15 +113,15 @@ def _setup_rock_and_till_with_threshold(self): erodibility.""" # Create field for erodibility - self.erody = self.grid.add_zeros("node", "substrate__erodibility") + self.erody = self._grid.add_zeros("node", "substrate__erodibility") # Create field for threshold values - self.threshold = self.grid.add_zeros( + self.threshold = self._grid.add_zeros( "node", "water_erosion_rule__threshold" ) # Create array for erodibility weighting function - self.erody_wt = np.zeros(self.grid.number_of_nodes) + self.erody_wt = np.zeros(self._grid.number_of_nodes) # set values correctly self._update_erodywt() @@ -125,7 +129,7 @@ def _setup_rock_and_till_with_threshold(self): def _update_erodywt(self): # Update the erodibility weighting function (this is "F") - core = self.grid.core_nodes + core = self._grid.core_nodes if self.contact_width > 0.0: self.erody_wt[core] = 1.0 / ( 1.0 diff --git a/terrainbento/bmi/__init__.py b/terrainbento/bmi/__init__.py new file mode 100644 index 00000000..eb63d9ad --- /dev/null +++ b/terrainbento/bmi/__init__.py @@ -0,0 +1,4 @@ +"""The Basic Model Interface.""" +from .bmi import BmiModel + +__all__ = ["BmiModel"] diff --git a/terrainbento/bmi/bmi.py b/terrainbento/bmi/bmi.py new file mode 100644 index 00000000..5cd70671 --- /dev/null +++ b/terrainbento/bmi/bmi.py @@ -0,0 +1,631 @@ +# -*- coding: utf-8 -*- + +import numpy + +from bmipy import Bmi + + +class BmiModel(Bmi): + def finalize(self) -> None: + """Perform tear-down tasks for the model. + + Perform all tasks that take place after exiting the model's time + loop. This typically includes deallocating memory, closing files and + printing reports. + """ + raise NotImplementedError("finalize") + + def get_component_name(self) -> str: + """Name of the component. + + Returns + ------- + str + The name of the component. + """ + raise NotImplementedError("get_component_name") + + def get_current_time(self) -> float: + """Current time of the model. + + Returns + ------- + float + The current model time. + """ + raise NotImplementedError("get_current_time") + + def get_end_time(self) -> float: + """End time of the model. + + Returns + ------- + float + The maximum model time. + """ + raise NotImplementedError("get_end_time") + + def get_grid_edge_count(self, grid: int) -> int: + """Get the number of edges in the grid. + + Parameters + ---------- + grid : int + A grid identifier. + + Returns + ------- + int + The total number of grid edges. + """ + raise NotImplementedError("get_grid_edge_count") + + def get_grid_edge_nodes( + self, grid: int, edge_nodes: numpy.ndarray + ) -> numpy.ndarray: + """Get the edge-node connectivity. + + Parameters + ---------- + grid : int + A grid identifier. + edge_nodes : ndarray of int, shape *(2 x nnodes,)* + A numpy array to place the edge-node connectivity. For each edge, + connectivity is given as node at edge tail, followed by node at + edge head. + + Returns + ------- + ndarray of int + The input numpy array that holds the edge-node connectivity. + """ + raise NotImplementedError("get_grid_edge_nodes") + + def get_grid_face_count(self, grid: int) -> int: + """Get the number of faces in the grid. + + Parameters + ---------- + grid : int + A grid identifier. + + Returns + ------- + int + The total number of grid faces. + """ + raise NotImplementedError("get_grid_face_count") + + def get_grid_face_nodes( + self, grid: int, face_nodes: numpy.ndarray + ) -> numpy.ndarray: + """Get the face-node connectivity. + + Parameters + ---------- + grid : int + A grid identifier. + face_nodes : ndarray of int + A numpy array to place the face-node connectivity. For each face, + the nodes (listed in a counter-clockwise direction) that form the + boundary of the face. + + Returns + ------- + ndarray of int + The input numpy array that holds the face-node connectivity. + """ + raise NotImplementedError("get_grid_face_nodes") + + def get_grid_node_count(self, grid: int) -> int: + """Get the number of nodes in the grid. + + Parameters + ---------- + grid : int + A grid identifier. + + Returns + ------- + int + The total number of grid nodes. + """ + raise NotImplementedError("get_grid_node_count") + + def get_grid_nodes_per_face( + self, grid: int, nodes_per_face: numpy.ndarray + ) -> numpy.ndarray: + """Get the number of nodes for each face. + + Parameters + ---------- + grid : int + A grid identifier. + nodes_per_face : ndarray of int, shape *(nfaces,)* + A numpy array to place the number of edges per face. + + Returns + ------- + ndarray of int + The input numpy array that holds the number of nodes per edge. + """ + raise NotImplementedError("get_grid_nodes_per_face") + + def get_grid_origin( + self, grid: int, origin: numpy.ndarray + ) -> numpy.ndarray: + """Get coordinates for the lower-left corner of the computational grid. + + Parameters + ---------- + grid : int + A grid identifier. + origin : ndarray of float, shape *(ndim,)* + A numpy array to hold the coordinates of the lower-left corner of + the grid. + + Returns + ------- + ndarray of float + The input numpy array that holds the coordinates of the grid's + lower-left corner. + """ + raise NotImplementedError("get_grid_origin") + + def get_grid_rank(self, grid: int) -> int: + """Get number of dimensions of the computational grid. + + Parameters + ---------- + grid : int + A grid identifier. + + Returns + ------- + int + Rank of the grid. + """ + raise NotImplementedError("get_grid_rank") + + def get_grid_shape(self, grid: int, shape: numpy.ndarray) -> numpy.ndarray: + """Get dimensions of the computational grid. + + Parameters + ---------- + grid : int + A grid identifier. + shape : ndarray of int, shape *(ndim,)* + A numpy array into which to place the shape of the grid. + + Returns + ------- + ndarray of int + The input numpy array that holds the grid's shape. + """ + raise NotImplementedError("get_grid_shape") + + def get_grid_size(self, grid: int) -> int: + """Get the total number of elements in the computational grid. + + Parameters + ---------- + grid : int + A grid identifier. + + Returns + ------- + int + Size of the grid. + """ + raise NotImplementedError("get_grid_size") + + def get_grid_spacing( + self, grid: int, spacing: numpy.ndarray + ) -> numpy.ndarray: + """Get distance between nodes of the computational grid. + + Parameters + ---------- + grid : int + A grid identifier. + spacing : ndarray of float, shape *(ndim,)* + A numpy array to hold the spacing between grid rows and columns. + + Returns + ------- + ndarray of float + The input numpy array that holds the grid's spacing. + """ + raise NotImplementedError("get_grid_spacing") + + def get_grid_type(self, grid: int) -> str: + """Get the grid type as a string. + + Parameters + ---------- + grid : int + A grid identifier. + + Returns + ------- + str + Type of grid as a string. + """ + raise NotImplementedError("get_grid_type") + + def get_grid_x(self, grid: int, x: numpy.ndarray) -> numpy.ndarray: + """Get coordinates of grid nodes in the x direction. + + Parameters + ---------- + grid : int + A grid identifier. + x : ndarray of float, shape *(nrows,)* + A numpy array to hold the x-coordinates of the grid node columns. + + Returns + ------- + ndarray of float + The input numpy array that holds the grid's column x-coordinates. + """ + raise NotImplementedError("get_grid_x") + + def get_grid_y(self, grid: int, y: numpy.ndarray) -> numpy.ndarray: + """Get coordinates of grid nodes in the y direction. + + Parameters + ---------- + grid : int + A grid identifier. + y : ndarray of float, shape *(ncols,)* + A numpy array to hold the y-coordinates of the grid node rows. + + Returns + ------- + ndarray of float + The input numpy array that holds the grid's row y-coordinates. + """ + raise NotImplementedError("get_grid_y") + + def get_grid_z(self, grid: int, z: numpy.ndarray) -> numpy.ndarray: + """Get coordinates of grid nodes in the z direction. + + Parameters + ---------- + grid : int + A grid identifier. + z : ndarray of float, shape *(nlayers,)* + A numpy array to hold the z-coordinates of the grid nodes layers. + + Returns + ------- + ndarray of float + The input numpy array that holds the grid's layer z-coordinates. + """ + raise NotImplementedError("get_grid_z") + + def get_input_var_names(self): + """List of a model's input variables. + + Input variable names must be CSDMS Standard Names, also known + as *long variable names*. + + Returns + ------- + list of str + The input variables for the model. + + Notes + ----- + Standard Names enable the CSDMS framework to determine whether + an input variable in one model is equivalent to, or compatible + with, an output variable in another model. This allows the + framework to automatically connect components. + + Standard Names do not have to be used within the model. + """ + raise NotImplementedError("get_input_var_names") + + def get_output_var_names(self): + """List of a model's output variables. + + Output variable names must be CSDMS Standard Names, also known + as *long variable names*. + + Returns + ------- + list of str + The output variables for the model. + """ + raise NotImplementedError("get_output_var_names") + + def get_start_time(self) -> float: + """Start time of the model. + + Model times should be of type float. + + Returns + ------- + float + The model start time. + """ + raise NotImplementedError("get_start_time") + + def get_time_step(self) -> float: + """Current time step of the model. + + The model time step should be of type float. + + Returns + ------- + float + The time step used in model. + """ + raise NotImplementedError("get_time_step") + + def get_time_units(self) -> str: + """Time units of the model. + + Returns + ------- + float + The model time unit; e.g., `days` or `s`. + + Notes + ----- + CSDMS uses the UDUNITS standard from Unidata. + """ + raise NotImplementedError("get_time_units") + + def get_value(self, name: str, dest: numpy.ndarray) -> numpy.ndarray: + """Get a copy of values of the given variable. + + This is a getter for the model, used to access the model's + current state. It returns a *copy* of a model variable, with + the return type, size and rank dependent on the variable. + + Parameters + ---------- + name : str + An input or output variable name, a CSDMS Standard Name. + dest : ndarray + A numpy array into which to place the values. + + Returns + ------- + ndarray + The same numpy array that was passed as an input buffer. + """ + raise NotImplementedError("get_value") + + def get_value_at_indices( + self, name: str, dest: numpy.ndarray, inds: numpy.ndarray + ) -> numpy.ndarray: + """Get values at particular indices. + + Parameters + ---------- + name : str + An input or output variable name, a CSDMS Standard Name. + dest : ndarray + A numpy array into which to place the values. + indices : array_like + The indices into the variable array. + + Returns + ------- + array_like + Value of the model variable at the given location. + """ + raise NotImplementedError("get_value_at_indices") + + def get_value_ptr(self, name: str) -> numpy.ndarray: + """Get a reference to values of the given variable. + + This is a getter for the model, used to access the model's + current state. It returns a reference to a model variable, + with the return type, size and rank dependent on the variable. + + Parameters + ---------- + name : str + An input or output variable name, a CSDMS Standard Name. + + Returns + ------- + array_like + A reference to a model variable. + """ + raise NotImplementedError("get_value_ptr") + + def get_var_grid(self, name: str) -> int: + """Get grid identifier for the given variable. + + Parameters + ---------- + name : str + An input or output variable name, a CSDMS Standard Name. + + Returns + ------- + int + The grid identifier. + """ + raise NotImplementedError("get_var_grid") + + def get_var_itemsize(self, name: str) -> int: + """Get memory use for each array element in bytes. + + Parameters + ---------- + name : str + An input or output variable name, a CSDMS Standard Name. + + Returns + ------- + int + Item size in bytes. + """ + raise NotImplementedError("get_var_itemsize") + + def get_var_location(self, name: str) -> str: + """Get the grid element type that the a given variable is defined on. + + The grid topology can be composed of *nodes*, *edges*, and *faces*. + + *node* + A point that has a coordinate pair or triplet: the most + basic element of the topology. + + *edge* + A line or curve bounded by two *nodes*. + + *face* + A plane or surface enclosed by a set of edges. In a 2D + horizontal application one may consider the word “polygon”, + but in the hierarchy of elements the word “face” is most common. + + Parameters + ---------- + name : str + An input or output variable name, a CSDMS Standard Name. + + Returns + ------- + str + The grid location on which the variable is defined. Must be one of + `"node"`, `"edge"`, or `"face"`. + + Notes + ----- + CSDMS uses the `ugrid conventions`_ to define unstructured grids. + + .. _ugrid conventions: http://ugrid-conventions.github.io/ugrid-conventions + """ + raise NotImplementedError("get_var_location") + + def get_var_nbytes(self, name: str) -> int: + """Get size, in bytes, of the given variable. + + Parameters + ---------- + name : str + An input or output variable name, a CSDMS Standard Name. + + Returns + ------- + int + The size of the variable, counted in bytes. + """ + raise NotImplementedError("get_var_nbytes") + + def get_var_type(self, name: str) -> str: + """Get data type of the given variable. + + Parameters + ---------- + name : str + An input or output variable name, a CSDMS Standard Name. + + Returns + ------- + str + The Python variable type; e.g., ``str``, ``int``, ``float``. + """ + raise NotImplementedError("get_var_type") + + def get_var_units(self, name: str) -> str: + """Get units of the given variable. + + Standard unit names, in lower case, should be used, such as + ``meters`` or ``seconds``. Standard abbreviations, like ``m`` for + meters, are also supported. For variables with compound units, + each unit name is separated by a single space, with exponents + other than 1 placed immediately after the name, as in ``m s-1`` + for velocity, ``W m-2`` for an energy flux, or ``km2`` for an + area. + + Parameters + ---------- + name : str + An input or output variable name, a CSDMS Standard Name. + + Returns + ------- + str + The variable units. + + Notes + ----- + CSDMS uses the `UDUNITS`_ standard from Unidata. + + .. _UDUNITS: http://www.unidata.ucar.edu/software/udunits + """ + raise NotImplementedError("get_var_units") + + def initialize(self, config_file: str) -> None: + """Perform startup tasks for the model. + + Perform all tasks that take place before entering the model's time + loop, including opening files and initializing the model state. Model + inputs are read from a text-based configuration file, specified by + `filename`. + + Parameters + ---------- + config_file : str, optional + The path to the model configuration file. + + Notes + ----- + Models should be refactored, if necessary, to use a + configuration file. CSDMS does not impose any constraint on + how configuration files are formatted, although YAML is + recommended. A template of a model's configuration file + with placeholder values is used by the BMI. + """ + raise NotImplementedError("initialize") + + def set_value(self, name: str, values: numpy.ndarray) -> None: + """Specify a new value for a model variable. + + This is the setter for the model, used to change the model's + current state. It accepts, through *src*, a new value for a + model variable, with the type, size and rank of *src* + dependent on the variable. + + Parameters + ---------- + var_name : str + An input or output variable name, a CSDMS Standard Name. + src : array_like + The new value for the specified variable. + """ + raise NotImplementedError("set_value") + + def set_value_at_indices( + self, name: str, inds: numpy.ndarray, src: numpy.ndarray + ) -> None: + """Specify a new value for a model variable at particular indices. + + Parameters + ---------- + var_name : str + An input or output variable name, a CSDMS Standard Name. + indices : array_like + The indices into the variable array. + src : array_like + The new value for the specified variable. + """ + raise NotImplementedError("set_value_at_indices") + + def update(self) -> None: + """Advance model state by one time step. + + Perform all tasks that take place within one pass through the model's + time loop. This typically includes incrementing all of the model's + state variables. If the model's state variables don't change in time, + then they can be computed by the :func:`initialize` method and this + method can return with no action. + """ + raise NotImplementedError("update") diff --git a/terrainbento/boundary_handlers/__init__.py b/terrainbento/boundary_handlers/__init__.py index f0c1ee39..965d6758 100644 --- a/terrainbento/boundary_handlers/__init__.py +++ b/terrainbento/boundary_handlers/__init__.py @@ -1,10 +1,10 @@ """Classes to assist with boundary conditions in the terrainbento package.""" from .capture_node_baselevel_handler import CaptureNodeBaselevelHandler -from .not_core_node_baselevel_handler import NotCoreNodeBaselevelHandler -from .single_node_baselevel_handler import SingleNodeBaselevelHandler from .generic_function_baselevel_handler import GenericFuncBaselevelHandler +from .not_core_node_baselevel_handler import NotCoreNodeBaselevelHandler from .precip_changer import PrecipChanger +from .single_node_baselevel_handler import SingleNodeBaselevelHandler __all__ = [ "CaptureNodeBaselevelHandler", diff --git a/terrainbento/clock/clock.py b/terrainbento/clock/clock.py index 9f5e0de7..71fd9fe6 100644 --- a/terrainbento/clock/clock.py +++ b/terrainbento/clock/clock.py @@ -1,5 +1,6 @@ """Clock sets the run duration and timestep in terrainbento model runs.""" +import cfunits import yaml @@ -60,7 +61,7 @@ def from_dict(cls, params): """ return cls(**params) - def __init__(self, start=0., step=10., stop=100.): + def __init__(self, start=0.0, step=10.0, stop=100.0, units="day"): """ Parameters ---------- @@ -70,7 +71,8 @@ def __init__(self, start=0., step=10., stop=100.): Model stop time. Default is 100. step : float, optional Model time step. Default is 10. - + units : str, optional + Default is "day". Examples -------- >>> from terrainbento import Clock @@ -94,9 +96,42 @@ def __init__(self, start=0., step=10., stop=100.): 2400.0 >>> clock.step 200.0 + + There are two ways to advance the model time stored in `clock.time`. + + >>> clock.time + 0.0 + + First, to advance by the step size provided, use `advance`: + + >>> clock.advance() + >>> clock.time + 200.0 + + To advance by an arbitrary time pass the `dt` value to `advance`: + + >>> clock.advance(2.) + >>> clock.time + 202.0 + + It is also possible to change the timestep. + + >>> clock.step = 18. + >>> clock.advance() + >>> clock.time + 220.0 + + And to change the stop time. + + >>> clock.stop = 1000. + >>> clock.stop + 1000.0 """ + # verify that unit is a valid CFUNITS + # raise ValueError() + try: - self.start = float(start) + self._start = float(start) except ValueError: msg = ( "Clock: Required parameter *start* is " @@ -105,7 +140,7 @@ def __init__(self, start=0., step=10., stop=100.): raise ValueError(msg) try: - self.step = float(step) + self._step = float(step) except ValueError: msg = ( "Clock: Required parameter *step* is " @@ -114,7 +149,7 @@ def __init__(self, start=0., step=10., stop=100.): raise ValueError(msg) try: - self.stop = float(stop) + self._stop = float(stop) except ValueError: msg = ( "Clock: Required parameter *stop* is " @@ -125,3 +160,51 @@ def __init__(self, start=0., step=10., stop=100.): if self.start > self.stop: msg = "Clock: *start* is larger than *stop*." raise ValueError(msg) + + self._time = 0.0 + + def advance(self, dt=None): + """Advance the time stepper by one time step or a provided value. + + Parameters + ---------- + dt : float, optional + Model time step. Default is to use the step provided at + instantiation. + """ + step = dt or self._step + self._time += step + if self._stop is not None and self._time > self._stop: + raise StopIteration() + + @property + def time(self): + """Current time.""" + return self._time + + @property + def start(self): + """Start time.""" + return self._start + + @property + def stop(self): + """Stop time.""" + return self._stop + + @stop.setter + def stop(self, new_val): + """Change the stop time.""" + if self._time > new_val: + raise ValueError("") + self._stop = new_val + + @property + def step(self): + """Time Step.""" + return self._step + + @step.setter + def step(self, new_val): + """Change the time step.""" + self._step = new_val diff --git a/terrainbento/derived_models/__init__.py b/terrainbento/derived_models/__init__.py index 1be57abf..d75d29c3 100644 --- a/terrainbento/derived_models/__init__.py +++ b/terrainbento/derived_models/__init__.py @@ -1,63 +1,63 @@ """Derived models in the terrainbento package.""" -from .model_basic import Basic -from .model_basicTh import BasicTh -from .model_basicDd import BasicDd -from .model_basicHy import BasicHy +from .model_basic import Basic from .model_basicCh import BasicCh -from .model_basicSt import BasicSt -from .model_basicVs import BasicVs -from .model_basicSa import BasicSa -from .model_basicRt import BasicRt +from .model_basicChRt import BasicChRt +from .model_basicChRtTh import BasicChRtTh +from .model_basicChSa import BasicChSa from .model_basicCv import BasicCv - +from .model_basicDd import BasicDd from .model_basicDdHy import BasicDdHy -from .model_basicStTh import BasicStTh +from .model_basicDdRt import BasicDdRt from .model_basicDdSt import BasicDdSt -from .model_basicHySt import BasicHySt -from .model_basicThVs import BasicThVs from .model_basicDdVs import BasicDdVs -from .model_basicHyVs import BasicHyVs -from .model_basicStVs import BasicStVs +from .model_basicHy import BasicHy +from .model_basicHyRt import BasicHyRt from .model_basicHySa import BasicHySa -from .model_basicChSa import BasicChSa -from .model_basicSaVs import BasicSaVs +from .model_basicHySt import BasicHySt +from .model_basicHyVs import BasicHyVs +from .model_basicRt import BasicRt +from .model_basicRtSa import BasicRtSa from .model_basicRtTh import BasicRtTh -from .model_basicDdRt import BasicDdRt -from .model_basicHyRt import BasicHyRt -from .model_basicChRt import BasicChRt from .model_basicRtVs import BasicRtVs -from .model_basicRtSa import BasicRtSa - -from .model_basicChRtTh import BasicChRtTh +from .model_basicSa import BasicSa +from .model_basicSaVs import BasicSaVs +from .model_basicSt import BasicSt +from .model_basicStTh import BasicStTh +from .model_basicStVs import BasicStVs +from .model_basicTh import BasicTh +from .model_basicThVs import BasicThVs +from .model_basicVs import BasicVs -__all__ = [ - "Basic", - "BasicTh", - "BasicDd", - "BasicHy", - "BasicCh", - "BasicSt", - "BasicVs", - "BasicSa", - "BasicRt", - "BasicCv", - "BasicDdHy", - "BasicStTh", - "BasicDdSt", - "BasicHySt", - "BasicThVs", - "BasicDdVs", - "BasicStVs", - "BasicHySa", - "BasicHyVs", - "BasicChSa", - "BasicSaVs", - "BasicRtTh", - "BasicDdRt", - "BasicHyRt", - "BasicChRt", - "BasicRtVs", - "BasicRtSa", - "BasicChRtTh", +MODELS = [ + Basic, + BasicTh, + BasicDd, + BasicHy, + BasicCh, + BasicSt, + BasicVs, + BasicSa, + BasicRt, + BasicCv, + BasicDdHy, + BasicStTh, + BasicDdSt, + BasicHySt, + BasicThVs, + BasicDdVs, + BasicStVs, + BasicHySa, + BasicHyVs, + BasicChSa, + BasicSaVs, + BasicRtTh, + BasicDdRt, + BasicHyRt, + BasicChRt, + BasicRtVs, + BasicRtSa, + BasicChRtTh, ] + +__all__ = [cls.__name__ for cls in MODELS] diff --git a/terrainbento/derived_models/model_basic.py b/terrainbento/derived_models/model_basic.py index 0f518853..277d1f7d 100644 --- a/terrainbento/derived_models/model_basic.py +++ b/terrainbento/derived_models/model_basic.py @@ -40,7 +40,51 @@ class Basic(ErosionModel): - ``topographic__elevation`` """ - _required_fields = ["topographic__elevation"] + _name = "Basic" + + _input_var_names = ("topographic__elevation", "water__unit_flux_in") + + _output_var_names = ("topographic__elevation",) + + _var_info = { + "topographic__elevation": { + "type": float, + "units": "m", + "at": "node", + "long_name": "Land surface topographic elevation", + "_FillValue": -999.0, + }, + "water__unit_flux_in": { + "type": float, + "units": "m/s", + "at": "node", + "long_name": ( + "Volume water per area per time input to each node." + ), + }, + } + + _param_info = { + "water_erodability": { + "type": float, + "units": None, # should I put in dependence on M? + "default": 0.001, # these values assume yr instead of second + "range": [1e-7, 1e-2], # these values assume yr instead of second + "long_name": "", + "definition": ( + "Water erodibility coefficient. Units depend on " + "values of m_sp and n_sp." + ), + }, + "regolith_transport_parameter": { + "type": float, + "units": "m/(s^2)", + "default": 0.1, # these values assume yr instead of second + "range": [1e-3, 1e-1], # these values assume yr instead of second + "long_name": "", + "definition": ("Hillslope diffusivity."), + }, + } def __init__( self, @@ -98,7 +142,7 @@ def __init__( will just run it one step. >>> model.run_one_step(1.) - >>> model.model_time + >>> model.clock.time 1.0 """ @@ -106,7 +150,7 @@ def __init__( super(Basic, self).__init__(clock, grid, **kwargs) # verify correct fields are present. - self._verify_fields(self._required_fields) + self._verify_fields(self._input_var_names) # Get Parameters: self.m = m_sp @@ -117,7 +161,7 @@ def __init__( # Instantiate a FastscapeEroder component self.eroder = FastscapeEroder( - self.grid, + self._grid, K_sp=self.K, m_sp=self.m, n_sp=self.n, @@ -126,7 +170,7 @@ def __init__( # Instantiate a LinearDiffuser component self.diffuser = LinearDiffuser( - self.grid, linear_diffusivity=self.regolith_transport_parameter + self._grid, linear_diffusivity=self.regolith_transport_parameter ) def run_one_step(self, step): diff --git a/terrainbento/derived_models/model_basicCh.py b/terrainbento/derived_models/model_basicCh.py index 77ed5e7d..1fb088d8 100644 --- a/terrainbento/derived_models/model_basicCh.py +++ b/terrainbento/derived_models/model_basicCh.py @@ -48,7 +48,11 @@ class BasicCh(ErosionModel): - ``topographic__elevation`` """ - _required_fields = ["topographic__elevation"] + _name = "BasicCh" + + _input_var_names = ("topographic__elevation", "water__unit_flux_in") + + _output_var_names = ("topographic__elevation",) def __init__( self, @@ -113,7 +117,7 @@ def __init__( will just run it one step. >>> model.run_one_step(1.) - >>> model.model_time + >>> model.clock.time 1.0 """ @@ -122,7 +126,7 @@ def __init__( super(BasicCh, self).__init__(clock, grid, **kwargs) # verify correct fields are present. - self._verify_fields(self._required_fields) + self._verify_fields(self._input_var_names) # Get Parameters and convert units if necessary: self.m = m_sp @@ -133,7 +137,7 @@ def __init__( # Instantiate a FastscapeEroder component self.eroder = FastscapeEroder( - self.grid, + self._grid, K_sp=self.K, m_sp=self.m, n_sp=self.n, @@ -142,7 +146,7 @@ def __init__( # Instantiate a NonLinearDiffuser component self.diffuser = TaylorNonLinearDiffuser( - self.grid, + self._grid, linear_diffusivity=regolith_transport_parameter, slope_crit=critical_slope, nterms=number_of_taylor_terms, diff --git a/terrainbento/derived_models/model_basicChRt.py b/terrainbento/derived_models/model_basicChRt.py index c2d02c42..0e1c626e 100644 --- a/terrainbento/derived_models/model_basicChRt.py +++ b/terrainbento/derived_models/model_basicChRt.py @@ -68,10 +68,15 @@ class BasicChRt(TwoLithologyErosionModel): - ``lithology_contact__elevation`` """ - _required_fields = [ + _name = "BasicChRt" + + _input_var_names = ( "topographic__elevation", "lithology_contact__elevation", - ] + "water__unit_flux_in", + ) + + _output_var_names = ("topographic__elevation",) def __init__( self, @@ -140,7 +145,7 @@ def __init__( will just run it one step. >>> model.run_one_step(1.) - >>> model.model_time + >>> model.clock.time 1.0 """ # Call ErosionModel"s init @@ -151,7 +156,7 @@ def __init__( # Instantiate a FastscapeEroder component self.eroder = FastscapeEroder( - self.grid, + self._grid, m_sp=self.m, n_sp=self.n, K_sp=self.erody, @@ -160,7 +165,7 @@ def __init__( # Instantiate a LinearDiffuser component self.diffuser = TaylorNonLinearDiffuser( - self.grid, + self._grid, linear_diffusivity=self.regolith_transport_parameter, slope_crit=critical_slope, nterms=number_of_taylor_terms, diff --git a/terrainbento/derived_models/model_basicChRtTh.py b/terrainbento/derived_models/model_basicChRtTh.py index 95f5ccf8..12def074 100644 --- a/terrainbento/derived_models/model_basicChRtTh.py +++ b/terrainbento/derived_models/model_basicChRtTh.py @@ -78,10 +78,15 @@ class BasicChRtTh(TwoLithologyErosionModel): - ``lithology_contact__elevation`` """ - _required_fields = [ + _name = "BasicChRtTh" + + _input_var_names = ( "topographic__elevation", "lithology_contact__elevation", - ] + "water__unit_flux_in", + ) + + _output_var_names = ("topographic__elevation",) def __init__( self, @@ -158,7 +163,7 @@ def __init__( will just run it one step. >>> model.run_one_step(1.) - >>> model.model_time + >>> model.clock.time 1.0 """ # Call ErosionModel"s init @@ -168,7 +173,7 @@ def __init__( raise ValueError("Model only supports n equals 1.") # verify correct fields are present. - self._verify_fields(self._required_fields) + self._verify_fields(self._input_var_names) # Save the threshold values for rock and till self.rock_thresh = water_erosion_rule_lower__threshold @@ -179,7 +184,7 @@ def __init__( # Instantiate a StreamPowerSmoothThresholdEroder component self.eroder = StreamPowerSmoothThresholdEroder( - self.grid, + self._grid, K_sp=self.erody, threshold_sp=self.threshold, m_sp=self.m, @@ -189,7 +194,7 @@ def __init__( # Instantiate a LinearDiffuser component self.diffuser = TaylorNonLinearDiffuser( - self.grid, + self._grid, linear_diffusivity=self.regolith_transport_parameter, slope_crit=critical_slope, nterms=number_of_taylor_terms, diff --git a/terrainbento/derived_models/model_basicChSa.py b/terrainbento/derived_models/model_basicChSa.py index a9e73791..c1c10073 100644 --- a/terrainbento/derived_models/model_basicChSa.py +++ b/terrainbento/derived_models/model_basicChSa.py @@ -74,7 +74,11 @@ class BasicChSa(ErosionModel): - ``soil__depth`` """ - _required_fields = ["topographic__elevation", "soil__depth"] + _name = "BasicChSa" + + _input_var_names = ("topographic__elevation", "soil__depth") + + _output_var_names = ("topographic__elevation", "soil__depth") def __init__( self, @@ -149,7 +153,7 @@ def __init__( will just run it one step. >>> model.run_one_step(10) - >>> model.model_time + >>> model.clock.time 10.0 """ @@ -158,20 +162,20 @@ def __init__( super(BasicChSa, self).__init__(clock, grid, **kwargs) # verify correct fields are present. - self._verify_fields(self._required_fields) + self._verify_fields(self._input_var_names) self.m = m_sp self.n = n_sp self.K = water_erodibility # Create bedrock elevation field - soil_thickness = self.grid.at_node["soil__depth"] - bedrock_elev = self.grid.add_zeros("node", "bedrock__elevation") + soil_thickness = self._grid.at_node["soil__depth"] + bedrock_elev = self._grid.add_zeros("node", "bedrock__elevation") bedrock_elev[:] = self.z - soil_thickness # Instantiate a FastscapeEroder component self.eroder = FastscapeEroder( - self.grid, + self._grid, K_sp=self.K, m_sp=self.m, n_sp=self.n, @@ -180,14 +184,14 @@ def __init__( # Instantiate a weathering component self.weatherer = ExponentialWeatherer( - self.grid, + self._grid, soil_production__maximum_rate=soil_production__maximum_rate, soil_production__decay_depth=soil_production__decay_depth, ) # Instantiate a soil-transport component self.diffuser = DepthDependentTaylorDiffuser( - self.grid, + self._grid, linear_diffusivity=regolith_transport_parameter, slope_crit=critical_slope, soil_transport_decay_depth=soil_transport_decay_depth, @@ -251,8 +255,8 @@ def run_one_step(self, step): # into bedrock has occurred, the bedrock elevation will be higher than # the actual elevation, so we simply re-set bedrock elevation to the # lower of itself or the current elevation. - b = self.grid.at_node["bedrock__elevation"] - b[:] = np.minimum(b, self.grid.at_node["topographic__elevation"]) + b = self._grid.at_node["bedrock__elevation"] + b[:] = np.minimum(b, self._grid.at_node["topographic__elevation"]) # Calculate regolith-production rate self.weatherer.calc_soil_prod_rate() diff --git a/terrainbento/derived_models/model_basicCv.py b/terrainbento/derived_models/model_basicCv.py index eb0c1fbc..49fd141f 100644 --- a/terrainbento/derived_models/model_basicCv.py +++ b/terrainbento/derived_models/model_basicCv.py @@ -48,7 +48,11 @@ class BasicCv(ErosionModel): - ``topographic__elevation`` """ - _required_fields = ["topographic__elevation"] + _name = "BasicCv" + + _input_var_names = ("topographic__elevation", "water__unit_flux_in") + + _output_var_names = ("topographic__elevation",) def __init__( self, @@ -114,7 +118,7 @@ def __init__( will just run it one step. >>> model.run_one_step(1.) - >>> model.model_time + >>> model.clock.time 1.0 """ @@ -122,7 +126,7 @@ def __init__( super(BasicCv, self).__init__(clock, grid, **kwargs) # verify correct fields are present. - self._verify_fields(self._required_fields) + self._verify_fields(self._input_var_names) self.m = m_sp self.n = n_sp @@ -143,7 +147,7 @@ def __init__( # Instantiate a FastscapeEroder component self.eroder = FastscapeEroder( - self.grid, + self._grid, K_sp=K[0], m_sp=self.m, n_sp=self.n, @@ -152,7 +156,7 @@ def __init__( # Instantiate a LinearDiffuser component self.diffuser = LinearDiffuser( - self.grid, linear_diffusivity=regolith_transport_parameter + self._grid, linear_diffusivity=regolith_transport_parameter ) def run_one_step(self, step): @@ -193,7 +197,7 @@ def run_one_step(self, step): )[0] # Update erosion based on climate - self.eroder.K = float(self.K_through_time(self.model_time)) + self.eroder.K = float(self.K_through_time(self.clock.time)) # Do some erosion (but not on the flooded nodes) self.eroder.run_one_step(step, flooded_nodes=flooded) diff --git a/terrainbento/derived_models/model_basicDd.py b/terrainbento/derived_models/model_basicDd.py index 37cd59f3..a779a92b 100644 --- a/terrainbento/derived_models/model_basicDd.py +++ b/terrainbento/derived_models/model_basicDd.py @@ -56,7 +56,11 @@ class BasicDd(ErosionModel): - ``topographic__elevation`` """ - _required_fields = ["topographic__elevation"] + _name = "BasicDd" + + _input_var_names = ("topographic__elevation", "water__unit_flux_in") + + _output_var_names = ("topographic__elevation",) def __init__( self, @@ -122,14 +126,14 @@ def __init__( will just run it one step. >>> model.run_one_step(1.) - >>> model.model_time + >>> model.clock.time 1.0 """ # Call ErosionModel"s init super(BasicDd, self).__init__(clock, grid, **kwargs) # verify correct fields are present. - self._verify_fields(self._required_fields) + self._verify_fields(self._input_var_names) # Get Parameters and convert units if necessary: self.m = m_sp @@ -144,14 +148,14 @@ def __init__( self.threshold_value = water_erosion_rule__threshold # Create a field for the (initial) erosion threshold - self.threshold = self.grid.add_zeros( + self.threshold = self._grid.add_zeros( "node", "water_erosion_rule__threshold" ) self.threshold[:] = self.threshold_value # Instantiate a FastscapeEroder component self.eroder = StreamPowerSmoothThresholdEroder( - self.grid, + self._grid, m_sp=self.m, n_sp=self.n, K_sp=self.K, @@ -166,7 +170,7 @@ def __init__( # Instantiate a LinearDiffuser component self.diffuser = LinearDiffuser( - self.grid, linear_diffusivity=regolith_transport_parameter + self._grid, linear_diffusivity=regolith_transport_parameter ) def update_erosion_threshold_values(self): @@ -191,9 +195,9 @@ def update_erosion_threshold_values(self): # The second line handles the case where there is growth, in which case # we want the threshold to stay at its initial value rather than # getting smaller. - cum_ero = self.grid.at_node["cumulative_elevation_change"] + cum_ero = self._grid.at_node["cumulative_elevation_change"] cum_ero[:] = ( - self.z - self.grid.at_node["initial_topographic__elevation"] + self.z - self._grid.at_node["initial_topographic__elevation"] ) self.threshold[:] = self.threshold_value - ( self.thresh_change_per_depth * cum_ero diff --git a/terrainbento/derived_models/model_basicDdHy.py b/terrainbento/derived_models/model_basicDdHy.py index f0580fcf..e5ebb442 100644 --- a/terrainbento/derived_models/model_basicDdHy.py +++ b/terrainbento/derived_models/model_basicDdHy.py @@ -68,7 +68,11 @@ class BasicDdHy(ErosionModel): - ``topographic__elevation`` """ - _required_fields = ["topographic__elevation"] + _name = "BasicDdHy" + + _input_var_names = ("topographic__elevation", "water__unit_flux_in") + + _output_var_names = ("topographic__elevation",) def __init__( self, @@ -79,7 +83,7 @@ def __init__( water_erodibility=0.0001, regolith_transport_parameter=0.1, water_erosion_rule__threshold=0.01, - water_erosion_rule__thresh_depth_derivative=0., + water_erosion_rule__thresh_depth_derivative=0.0, settling_velocity=0.001, sediment_porosity=0.3, fraction_fines=0.5, @@ -150,14 +154,14 @@ def __init__( will just run it one step. >>> model.run_one_step(1.) - >>> model.model_time + >>> model.clock.time 1.0 """ # Call ErosionModel"s init super(BasicDdHy, self).__init__(clock, grid, **kwargs) # verify correct fields are present. - self._verify_fields(self._required_fields) + self._verify_fields(self._input_var_names) # Get Parameters and convert units if necessary: self.m = m_sp @@ -166,14 +170,14 @@ def __init__( self.sp_crit = water_erosion_rule__threshold # Create a field for the (initial) erosion threshold - self.threshold = self.grid.add_zeros( + self.threshold = self._grid.add_zeros( "node", "water_erosion_rule__threshold" ) self.threshold[:] = self.sp_crit # starting value # Instantiate an ErosionDeposition component self.eroder = ErosionDeposition( - self.grid, + self._grid, K=self.K, F_f=fraction_fines, phi=sediment_porosity, @@ -192,7 +196,7 @@ def __init__( # Instantiate a LinearDiffuser component self.diffuser = LinearDiffuser( - self.grid, linear_diffusivity=regolith_transport_parameter + self._grid, linear_diffusivity=regolith_transport_parameter ) def run_one_step(self, step): @@ -234,9 +238,9 @@ def run_one_step(self, step): )[0] # Calculate cumulative erosion and update threshold - cum_ero = self.grid.at_node["cumulative_elevation_change"] + cum_ero = self._grid.at_node["cumulative_elevation_change"] cum_ero[:] = ( - self.z - self.grid.at_node["initial_topographic__elevation"] + self.z - self._grid.at_node["initial_topographic__elevation"] ) self.threshold[:] = self.sp_crit - ( self.thresh_change_per_depth * cum_ero diff --git a/terrainbento/derived_models/model_basicDdRt.py b/terrainbento/derived_models/model_basicDdRt.py index e5d343c4..11e0444b 100644 --- a/terrainbento/derived_models/model_basicDdRt.py +++ b/terrainbento/derived_models/model_basicDdRt.py @@ -74,17 +74,22 @@ class BasicDdRt(TwoLithologyErosionModel): - ``lithology_contact__elevation`` """ - _required_fields = [ + _name = "BasicDdRt" + + _input_var_names = ( "topographic__elevation", "lithology_contact__elevation", - ] + "water__unit_flux_in", + ) + + _output_var_names = ("topographic__elevation",) def __init__( self, clock, grid, water_erosion_rule__threshold=0.01, - water_erosion_rule__thresh_depth_derivative=0., + water_erosion_rule__thresh_depth_derivative=0.0, **kwargs ): """ @@ -147,14 +152,14 @@ def __init__( will just run it one step. >>> model.run_one_step(1.) - >>> model.model_time + >>> model.clock.time 1.0 """ # Call ErosionModel"s init super(BasicDdRt, self).__init__(clock, grid, **kwargs) # verify correct fields are present. - self._verify_fields(self._required_fields) + self._verify_fields(self._input_var_names) if float(self.n) != 1.0: raise ValueError("Model only supports n equals 1.") @@ -165,14 +170,14 @@ def __init__( self._setup_rock_and_till() # Create a field for the (initial) erosion threshold - self.threshold = self.grid.add_zeros( + self.threshold = self._grid.add_zeros( "node", "water_erosion_rule__threshold" ) self.threshold[:] = self.threshold_value # Instantiate a StreamPowerSmoothThresholdEroder component self.eroder = StreamPowerSmoothThresholdEroder( - self.grid, + self._grid, K_sp=self.erody, m_sp=self.m, n_sp=self.n, @@ -187,7 +192,7 @@ def __init__( # Instantiate a LinearDiffuser component self.diffuser = LinearDiffuser( - self.grid, linear_diffusivity=self.regolith_transport_parameter + self._grid, linear_diffusivity=self.regolith_transport_parameter ) def _update_erosion_threshold_values(self): @@ -199,9 +204,9 @@ def _update_erosion_threshold_values(self): # The second line handles the case where there is growth, in which case # we want the threshold to stay at its initial value rather than # getting smaller. - cum_ero = self.grid.at_node["cumulative_elevation_change"] + cum_ero = self._grid.at_node["cumulative_elevation_change"] cum_ero[:] = ( - self.z - self.grid.at_node["initial_topographic__elevation"] + self.z - self._grid.at_node["initial_topographic__elevation"] ) self.threshold[:] = self.threshold_value - ( self.thresh_change_per_depth * cum_ero diff --git a/terrainbento/derived_models/model_basicDdSt.py b/terrainbento/derived_models/model_basicDdSt.py index 8f4ecdd0..0760e2cd 100644 --- a/terrainbento/derived_models/model_basicDdSt.py +++ b/terrainbento/derived_models/model_basicDdSt.py @@ -67,7 +67,11 @@ class BasicDdSt(StochasticErosionModel): - ``topographic__elevation`` """ - _required_fields = ["topographic__elevation"] + _name = "BasicDdSt" + + _input_var_names = ("topographic__elevation",) + + _output_var_names = ("topographic__elevation",) def __init__( self, @@ -78,7 +82,7 @@ def __init__( water_erodibility=0.0001, regolith_transport_parameter=0.1, water_erosion_rule__threshold=0.01, - water_erosion_rule__thresh_depth_derivative=0., + water_erosion_rule__thresh_depth_derivative=0.0, infiltration_capacity=1.0, **kwargs ): @@ -135,7 +139,7 @@ def __init__( will just run it one step. >>> model.run_one_step(1.) - >>> model.model_time + >>> model.clock.time 1.0 """ @@ -143,7 +147,7 @@ def __init__( super(BasicDdSt, self).__init__(clock, grid, **kwargs) # verify correct fields are present. - self._verify_fields(self._required_fields) + self._verify_fields(self._input_var_names) # Get Parameters: self.m = m_sp @@ -165,14 +169,14 @@ def __init__( self.flow_accumulator.run_one_step() # Create a field for the (initial) erosion threshold - self.threshold = self.grid.add_zeros( + self.threshold = self._grid.add_zeros( "node", "water_erosion_rule__threshold" ) self.threshold[:] = self.threshold_value # Instantiate a FastscapeEroder component self.eroder = StreamPowerSmoothThresholdEroder( - self.grid, + self._grid, m_sp=self.m, n_sp=self.n, K_sp=self.K, @@ -182,14 +186,14 @@ def __init__( # Instantiate a LinearDiffuser component self.diffuser = LinearDiffuser( - self.grid, linear_diffusivity=regolith_transport_parameter + self._grid, linear_diffusivity=regolith_transport_parameter ) def update_threshold_field(self): """Update the threshold based on cumulative erosion depth.""" - cum_ero = self.grid.at_node["cumulative_elevation_change"] + cum_ero = self._grid.at_node["cumulative_elevation_change"] cum_ero[:] = ( - self.z - self.grid.at_node["initial_topographic__elevation"] + self.z - self._grid.at_node["initial_topographic__elevation"] ) self.threshold[:] = self.threshold_value - ( self.thresh_change_per_depth * cum_ero diff --git a/terrainbento/derived_models/model_basicDdVs.py b/terrainbento/derived_models/model_basicDdVs.py index 8a88145e..c3a80f08 100644 --- a/terrainbento/derived_models/model_basicDdVs.py +++ b/terrainbento/derived_models/model_basicDdVs.py @@ -67,7 +67,11 @@ class BasicDdVs(ErosionModel): - ``soil__depth`` """ - _required_fields = ["topographic__elevation", "soil__depth"] + _name = "BasicDdVs" + + _input_var_names = ("topographic__elevation", "soil__depth") + + _output_var_names = ("topographic__elevation",) def __init__( self, @@ -78,7 +82,7 @@ def __init__( water_erodibility=0.0001, regolith_transport_parameter=0.1, water_erosion_rule__threshold=0.01, - water_erosion_rule__thresh_depth_derivative=0., + water_erosion_rule__thresh_depth_derivative=0.0, hydraulic_conductivity=0.1, **kwargs ): @@ -137,7 +141,7 @@ def __init__( will just run it one step. >>> model.run_one_step(1.) - >>> model.model_time + >>> model.clock.time 1.0 """ @@ -148,7 +152,7 @@ def __init__( self._ensure_precip_runoff_are_vanilla(vsa_precip=True) # verify correct fields are present. - self._verify_fields(self._required_fields) + self._verify_fields(self._input_var_names) if float(n_sp) != 1.0: raise ValueError("Model BasicDdVs only supports n = 1.") @@ -159,17 +163,17 @@ def __init__( self.threshold_value = water_erosion_rule__threshold # Get the effective-area parameter - self._Kdx = hydraulic_conductivity * self.grid.dx + self._Kdx = hydraulic_conductivity * self._grid.dx # Create a field for the (initial) erosion threshold - self.threshold = self.grid.add_zeros( + self.threshold = self._grid.add_zeros( "node", "water_erosion_rule__threshold" ) self.threshold[:] = self.threshold_value # Instantiate a FastscapeEroder component self.eroder = StreamPowerSmoothThresholdEroder( - self.grid, + self._grid, use_Q="surface_water__discharge", K_sp=self.K, m_sp=self.m, @@ -184,27 +188,27 @@ def __init__( # Instantiate a LinearDiffuser component self.diffuser = LinearDiffuser( - self.grid, linear_diffusivity=regolith_transport_parameter + self._grid, linear_diffusivity=regolith_transport_parameter ) def _calc_effective_drainage_area(self): """Calculate and store effective drainage area.""" - area = self.grid.at_node["drainage_area"] - slope = self.grid.at_node["topographic__steepest_slope"] - cores = self.grid.core_nodes + area = self._grid.at_node["drainage_area"] + slope = self._grid.at_node["topographic__steepest_slope"] + cores = self._grid.core_nodes sat_param = ( self._Kdx - * self.grid.at_node["soil__depth"] - / self.grid.at_node["rainfall__flux"] + * self._grid.at_node["soil__depth"] + / self._grid.at_node["rainfall__flux"] ) eff_area = area[cores] * ( np.exp(-sat_param[cores] * slope[cores] / area[cores]) ) - self.grid.at_node["surface_water__discharge"][cores] = eff_area + self._grid.at_node["surface_water__discharge"][cores] = eff_area def run_one_step(self, step): """Advance model **BasicVs** for one time-step of duration step. @@ -249,7 +253,7 @@ def run_one_step(self, step): )[0] # Zero out effective area in flooded nodes - self.grid.at_node["surface_water__discharge"][flooded] = 0.0 + self._grid.at_node["surface_water__discharge"][flooded] = 0.0 # Set the erosion threshold. # @@ -258,9 +262,9 @@ def run_one_step(self, step): # The second line handles the case where there is growth, in which case # we want the threshold to stay at its initial value rather than # getting smaller. - cum_ero = self.grid.at_node["cumulative_elevation_change"] + cum_ero = self._grid.at_node["cumulative_elevation_change"] cum_ero[:] = ( - self.z - self.grid.at_node["initial_topographic__elevation"] + self.z - self._grid.at_node["initial_topographic__elevation"] ) self.threshold[:] = self.threshold_value - ( self.thresh_change_per_depth * cum_ero diff --git a/terrainbento/derived_models/model_basicHy.py b/terrainbento/derived_models/model_basicHy.py index bfce72cc..aef94546 100644 --- a/terrainbento/derived_models/model_basicHy.py +++ b/terrainbento/derived_models/model_basicHy.py @@ -50,7 +50,11 @@ class BasicHy(ErosionModel): - ``topographic__elevation`` """ - _required_fields = ["topographic__elevation"] + _name = "BasicHy" + + _input_var_names = ("topographic__elevation", "water__unit_flux_in") + + _output_var_names = ("topographic__elevation",) def __init__( self, @@ -124,7 +128,7 @@ def __init__( will just run it one step. >>> model.run_one_step(1.) - >>> model.model_time + >>> model.clock.time 1.0 """ @@ -133,7 +137,7 @@ def __init__( super(BasicHy, self).__init__(clock, grid, **kwargs) # verify correct fields are present. - self._verify_fields(self._required_fields) + self._verify_fields(self._input_var_names) # Get Parameters self.m = m_sp @@ -142,7 +146,7 @@ def __init__( # Instantiate a Space component self.eroder = ErosionDeposition( - self.grid, + self._grid, K=self.K, phi=sediment_porosity, F_f=fraction_fines, @@ -155,7 +159,7 @@ def __init__( # Instantiate a LinearDiffuser component self.diffuser = LinearDiffuser( - self.grid, linear_diffusivity=regolith_transport_parameter + self._grid, linear_diffusivity=regolith_transport_parameter ) def run_one_step(self, step): diff --git a/terrainbento/derived_models/model_basicHyRt.py b/terrainbento/derived_models/model_basicHyRt.py index 818db1bc..2c5e63ce 100644 --- a/terrainbento/derived_models/model_basicHyRt.py +++ b/terrainbento/derived_models/model_basicHyRt.py @@ -69,10 +69,15 @@ class BasicHyRt(TwoLithologyErosionModel): - ``lithology_contact__elevation`` """ - _required_fields = [ + _name = "BasicHyRt" + + _input_var_names = ( "topographic__elevation", "lithology_contact__elevation", - ] + "water__unit_flux_in", + ) + + _output_var_names = ("topographic__elevation",) def __init__( self, @@ -150,25 +155,25 @@ def __init__( will just run it one step. >>> model.run_one_step(1.) - >>> model.model_time + >>> model.clock.time 1.0 """ # Call ErosionModel"s init super(BasicHyRt, self).__init__(clock, grid, **kwargs) # verify correct fields are present. - self._verify_fields(self._required_fields) + self._verify_fields(self._input_var_names) # Save the threshold values for rock and till - self.rock_thresh = 0. - self.till_thresh = 0. + self.rock_thresh = 0.0 + self.till_thresh = 0.0 # Set up rock-till boundary and associated grid fields. self._setup_rock_and_till_with_threshold() # Instantiate an ErosionDeposition ("hybrid") component self.eroder = ErosionDeposition( - self.grid, + self._grid, K="substrate__erodibility", F_f=fraction_fines, phi=sediment_porosity, @@ -181,7 +186,7 @@ def __init__( # Instantiate a LinearDiffuser component self.diffuser = LinearDiffuser( - self.grid, linear_diffusivity=self.regolith_transport_parameter + self._grid, linear_diffusivity=self.regolith_transport_parameter ) def run_one_step(self, step): diff --git a/terrainbento/derived_models/model_basicHySa.py b/terrainbento/derived_models/model_basicHySa.py index b395486f..3120d6b1 100644 --- a/terrainbento/derived_models/model_basicHySa.py +++ b/terrainbento/derived_models/model_basicHySa.py @@ -76,7 +76,11 @@ class BasicHySa(ErosionModel): - ``soil__depth`` """ - _required_fields = ["topographic__elevation", "soil__depth"] + _name = "BasicHySa" + + _input_var_names = ("topographic__elevation", "soil__depth") + + _output_var_names = ("topographic__elevation", "soil__depth") def __init__( self, @@ -166,7 +170,7 @@ def __init__( will just run it one step. >>> model.run_one_step(1.) - >>> model.model_time + >>> model.clock.time 1.0 """ @@ -174,10 +178,10 @@ def __init__( super(BasicHySa, self).__init__(clock, grid, **kwargs) # verify correct fields are present. - self._verify_fields(self._required_fields) + self._verify_fields(self._input_var_names) - soil_thickness = self.grid.at_node["soil__depth"] - bedrock_elev = self.grid.add_zeros("node", "bedrock__elevation") + soil_thickness = self._grid.at_node["soil__depth"] + bedrock_elev = self._grid.add_zeros("node", "bedrock__elevation") bedrock_elev[:] = self.z - soil_thickness self.m = m_sp @@ -187,7 +191,7 @@ def __init__( # Instantiate a SPACE component self.eroder = Space( - self.grid, + self._grid, K_sed=self.K_sed, K_br=self.K_br, sp_crit_br=sp_crit_br, @@ -204,20 +208,20 @@ def __init__( # Instantiate diffusion and weathering components self.diffuser = DepthDependentDiffuser( - self.grid, + self._grid, linear_diffusivity=regolith_transport_parameter, soil_transport_decay_depth=soil_transport_decay_depth, ) self.weatherer = ExponentialWeatherer( - self.grid, + self._grid, soil_production__maximum_rate=soil_production__maximum_rate, soil_production__decay_depth=soil_production__decay_depth, ) - self.grid.at_node["soil__depth"][:] = ( - self.grid.at_node["topographic__elevation"] - - self.grid.at_node["bedrock__elevation"] + self._grid.at_node["soil__depth"][:] = ( + self._grid.at_node["topographic__elevation"] + - self._grid.at_node["bedrock__elevation"] ) def run_one_step(self, step): @@ -273,8 +277,8 @@ def run_one_step(self, step): # into bedrock has occurred, the bedrock elevation will be higher than # the actual elevation, so we simply re-set bedrock elevation to the # lower of itself or the current elevation. - b = self.grid.at_node["bedrock__elevation"] - b[:] = np.minimum(b, self.grid.at_node["topographic__elevation"]) + b = self._grid.at_node["bedrock__elevation"] + b[:] = np.minimum(b, self._grid.at_node["topographic__elevation"]) # Calculate regolith-production rate self.weatherer.calc_soil_prod_rate() @@ -290,10 +294,10 @@ def run_one_step(self, step): def check_stability(self): """Check model stability and exit if unstable.""" - fields = self.grid.at_node.keys() + fields = self._grid.at_node.keys() for f in fields: - if np.any(np.isnan(self.grid.at_node[f])) or np.any( - np.isinf(self.grid.at_node[f]) + if np.any(np.isnan(self._grid.at_node[f])) or np.any( + np.isinf(self._grid.at_node[f]) ): raise SystemExit( "terrainbento ModelHySa: Model became unstable" diff --git a/terrainbento/derived_models/model_basicHySt.py b/terrainbento/derived_models/model_basicHySt.py index 6b0d371d..72fdf7e3 100644 --- a/terrainbento/derived_models/model_basicHySt.py +++ b/terrainbento/derived_models/model_basicHySt.py @@ -56,7 +56,11 @@ class BasicHySt(StochasticErosionModel): - ``topographic__elevation`` """ - _required_fields = ["topographic__elevation"] + _name = "BasicHySt" + + _input_var_names = ("topographic__elevation",) + + _output_var_names = ("topographic__elevation",) def __init__( self, @@ -132,7 +136,7 @@ def __init__( will just run it one step. >>> model.run_one_step(1.) - >>> model.model_time + >>> model.clock.time 1.0 """ @@ -140,7 +144,7 @@ def __init__( super(BasicHySt, self).__init__(clock, grid, **kwargs) # verify correct fields are present. - self._verify_fields(self._required_fields) + self._verify_fields(self._input_var_names) # Get Parameters: self.m = m_sp @@ -156,7 +160,7 @@ def __init__( # Instantiate an ErosionDeposition component self.eroder = ErosionDeposition( - self.grid, + self._grid, K=self.K, F_f=fraction_fines, phi=sediment_porosity, @@ -169,7 +173,7 @@ def __init__( # Instantiate a LinearDiffuser component self.diffuser = LinearDiffuser( - self.grid, linear_diffusivity=regolith_transport_parameter + self._grid, linear_diffusivity=regolith_transport_parameter ) def run_one_step(self, step): diff --git a/terrainbento/derived_models/model_basicHyVs.py b/terrainbento/derived_models/model_basicHyVs.py index c8339872..efe2b77a 100644 --- a/terrainbento/derived_models/model_basicHyVs.py +++ b/terrainbento/derived_models/model_basicHyVs.py @@ -62,7 +62,11 @@ class BasicHyVs(ErosionModel): - ``soil__depth`` """ - _required_fields = ["topographic__elevation", "soil__depth"] + _name = "BasicHyVs" + + _input_var_names = ("topographic__elevation", "soil__depth") + + _output_var_names = ("topographic__elevation",) def __init__( self, @@ -140,7 +144,7 @@ def __init__( will just run it one step. >>> model.run_one_step(1.) - >>> model.model_time + >>> model.clock.time 1.0 """ @@ -152,18 +156,18 @@ def __init__( self._ensure_precip_runoff_are_vanilla(vsa_precip=True) # verify correct fields are present. - self._verify_fields(self._required_fields) + self._verify_fields(self._input_var_names) self.m = m_sp self.n = n_sp self.K = water_erodibility # Get the effective-area parameter - self._Kdx = hydraulic_conductivity * self.grid.dx + self._Kdx = hydraulic_conductivity * self._grid.dx # Instantiate a SPACE component self.eroder = ErosionDeposition( - self.grid, + self._grid, K=self.K, F_f=fraction_fines, phi=sediment_porosity, @@ -176,27 +180,27 @@ def __init__( # Instantiate a LinearDiffuser component self.diffuser = LinearDiffuser( - self.grid, linear_diffusivity=regolith_transport_parameter + self._grid, linear_diffusivity=regolith_transport_parameter ) def _calc_effective_drainage_area(self): """Calculate and store effective drainage area.""" - area = self.grid.at_node["drainage_area"] - slope = self.grid.at_node["topographic__steepest_slope"] - cores = self.grid.core_nodes + area = self._grid.at_node["drainage_area"] + slope = self._grid.at_node["topographic__steepest_slope"] + cores = self._grid.core_nodes sat_param = ( self._Kdx - * self.grid.at_node["soil__depth"] - / self.grid.at_node["rainfall__flux"] + * self._grid.at_node["soil__depth"] + / self._grid.at_node["rainfall__flux"] ) eff_area = area[cores] * ( np.exp(-sat_param[cores] * slope[cores] / area[cores]) ) - self.grid.at_node["surface_water__discharge"][cores] = eff_area + self._grid.at_node["surface_water__discharge"][cores] = eff_area def run_one_step(self, step): """Advance model **BasicVs** for one time-step of duration step. @@ -241,7 +245,7 @@ def run_one_step(self, step): )[0] # Zero out effective area in flooded nodes - self.grid.at_node["surface_water__discharge"][flooded] = 0.0 + self._grid.at_node["surface_water__discharge"][flooded] = 0.0 # Do some erosion # (if we're varying K through time, update that first) diff --git a/terrainbento/derived_models/model_basicRt.py b/terrainbento/derived_models/model_basicRt.py index 6c4e77fd..e403877c 100644 --- a/terrainbento/derived_models/model_basicRt.py +++ b/terrainbento/derived_models/model_basicRt.py @@ -61,10 +61,15 @@ class BasicRt(TwoLithologyErosionModel): - ``lithology_contact__elevation`` """ - _required_fields = [ + _name = "BasicRt" + + _input_var_names = ( "topographic__elevation", "lithology_contact__elevation", - ] + "water__unit_flux_in", + ) + + _output_var_names = ("topographic__elevation",) def __init__(self, clock, grid, **kwargs): """ @@ -121,7 +126,7 @@ def __init__(self, clock, grid, **kwargs): will just run it one step. >>> model.run_one_step(1.) - >>> model.model_time + >>> model.clock.time 1.0 """ @@ -129,14 +134,14 @@ def __init__(self, clock, grid, **kwargs): super(BasicRt, self).__init__(clock, grid, **kwargs) # verify correct fields are present. - self._verify_fields(self._required_fields) + self._verify_fields(self._input_var_names) # Set up rock-till boundary and associated grid fields. self._setup_rock_and_till() # Instantiate a FastscapeEroder component self.eroder = FastscapeEroder( - self.grid, + self._grid, K_sp=self.erody, m_sp=self.m, n_sp=self.n, @@ -145,7 +150,7 @@ def __init__(self, clock, grid, **kwargs): # Instantiate a LinearDiffuser component self.diffuser = LinearDiffuser( - self.grid, linear_diffusivity=self.regolith_transport_parameter + self._grid, linear_diffusivity=self.regolith_transport_parameter ) def run_one_step(self, step): diff --git a/terrainbento/derived_models/model_basicRtSa.py b/terrainbento/derived_models/model_basicRtSa.py index daebea97..21f90a47 100644 --- a/terrainbento/derived_models/model_basicRtSa.py +++ b/terrainbento/derived_models/model_basicRtSa.py @@ -89,11 +89,16 @@ class BasicRtSa(TwoLithologyErosionModel): - ``soil__depth`` """ - _required_fields = [ + _name = "BasicRtSa" + + _input_var_names = ( "topographic__elevation", "lithology_contact__elevation", + "water__unit_flux_in", "soil__depth", - ] + ) + + _output_var_names = ("topographic__elevation", "soil__depth") def __init__( self, @@ -165,7 +170,7 @@ def __init__( will just run it one step. >>> model.run_one_step(1.) - >>> model.model_time + >>> model.clock.time 1.0 """ @@ -173,33 +178,33 @@ def __init__( super(BasicRtSa, self).__init__(clock, grid, **kwargs) # verify correct fields are present. - self._verify_fields(self._required_fields) + self._verify_fields(self._input_var_names) # Set up rock-till boundary and associated grid fields. self._setup_rock_and_till() # Instantiate a FastscapeEroder component self.eroder = FastscapeEroder( - self.grid, + self._grid, K_sp=self.erody, m_sp=self.m, n_sp=self.n, discharge_name="surface_water__discharge", ) - soil_thickness = self.grid.at_node["soil__depth"] - bedrock_elev = self.grid.add_zeros("node", "bedrock__elevation") + soil_thickness = self._grid.at_node["soil__depth"] + bedrock_elev = self._grid.add_zeros("node", "bedrock__elevation") bedrock_elev[:] = self.z - soil_thickness # Instantiate diffusion and weathering components self.diffuser = DepthDependentDiffuser( - self.grid, + self._grid, linear_diffusivity=self.regolith_transport_parameter, soil_transport_decay_depth=soil_transport_decay_depth, ) self.weatherer = ExponentialWeatherer( - self.grid, + self._grid, soil_production__maximum_rate=soil_production__maximum_rate, soil_production__decay_depth=soil_production__decay_depth, ) @@ -258,8 +263,8 @@ def run_one_step(self, step): # into bedrock has occurred, the bedrock elevation will be higher than # the actual elevation, so we simply re-set bedrock elevation to the # lower of itself or the current elevation. - b = self.grid.at_node["bedrock__elevation"] - b[:] = np.minimum(b, self.grid.at_node["topographic__elevation"]) + b = self._grid.at_node["bedrock__elevation"] + b[:] = np.minimum(b, self._grid.at_node["topographic__elevation"]) # Calculate regolith-production rate self.weatherer.calc_soil_prod_rate() diff --git a/terrainbento/derived_models/model_basicRtTh.py b/terrainbento/derived_models/model_basicRtTh.py index ede48517..6f5f8d5e 100644 --- a/terrainbento/derived_models/model_basicRtTh.py +++ b/terrainbento/derived_models/model_basicRtTh.py @@ -72,17 +72,22 @@ class BasicRtTh(TwoLithologyErosionModel): - ``lithology_contact__elevation`` """ - _required_fields = [ + _name = "BasicRtTh" + + _input_var_names = ( "topographic__elevation", "lithology_contact__elevation", - ] + "water__unit_flux_in", + ) + + _output_var_names = ("topographic__elevation",) def __init__( self, clock, grid, - water_erosion_rule_upper__threshold=1., - water_erosion_rule_lower__threshold=1., + water_erosion_rule_upper__threshold=1.0, + water_erosion_rule_lower__threshold=1.0, **kwargs ): """ @@ -145,7 +150,7 @@ def __init__( will just run it one step. >>> model.run_one_step(1.) - >>> model.model_time + >>> model.clock.time 1.0 """ @@ -156,7 +161,7 @@ def __init__( raise ValueError("Model only supports n equals 1.") # verify correct fields are present. - self._verify_fields(self._required_fields) + self._verify_fields(self._input_var_names) # Save the threshold values for rock and till self.rock_thresh = water_erosion_rule_lower__threshold @@ -167,7 +172,7 @@ def __init__( # Instantiate a StreamPowerSmoothThresholdEroder component self.eroder = StreamPowerSmoothThresholdEroder( - self.grid, + self._grid, K_sp=self.erody, threshold_sp=self.threshold, m_sp=self.m, @@ -177,7 +182,7 @@ def __init__( # Instantiate a LinearDiffuser component self.diffuser = LinearDiffuser( - self.grid, linear_diffusivity=self.regolith_transport_parameter + self._grid, linear_diffusivity=self.regolith_transport_parameter ) def run_one_step(self, step): diff --git a/terrainbento/derived_models/model_basicRtVs.py b/terrainbento/derived_models/model_basicRtVs.py index 13ef0b1d..ba1f427b 100644 --- a/terrainbento/derived_models/model_basicRtVs.py +++ b/terrainbento/derived_models/model_basicRtVs.py @@ -72,11 +72,16 @@ class BasicRtVs(TwoLithologyErosionModel): - ``soil__depth`` """ - _required_fields = [ + _name = "BasicRtVs" + + _input_var_names = ( "topographic__elevation", "lithology_contact__elevation", + "water__unit_flux_in", "soil__depth", - ] + ) + + _output_var_names = ("topographic__elevation",) def __init__(self, clock, grid, hydraulic_conductivity=0.1, **kwargs): """ @@ -136,7 +141,7 @@ def __init__(self, clock, grid, hydraulic_conductivity=0.1, **kwargs): will just run it one step. >>> model.run_one_step(1.) - >>> model.model_time + >>> model.clock.time 1.0 """ @@ -147,17 +152,17 @@ def __init__(self, clock, grid, hydraulic_conductivity=0.1, **kwargs): self._ensure_precip_runoff_are_vanilla() # verify correct fields are present. - self._verify_fields(self._required_fields) + self._verify_fields(self._input_var_names) # Set up rock-till boundary and associated grid fields. self._setup_rock_and_till() # Get the effective-area parameter - self._Kdx = hydraulic_conductivity * self.grid.dx + self._Kdx = hydraulic_conductivity * self._grid.dx # Instantiate a FastscapeEroder component self.eroder = FastscapeEroder( - self.grid, + self._grid, K_sp=self.erody, m_sp=self.m, n_sp=self.n, @@ -166,7 +171,7 @@ def __init__(self, clock, grid, hydraulic_conductivity=0.1, **kwargs): # Instantiate a LinearDiffuser component self.diffuser = LinearDiffuser( - self.grid, linear_diffusivity=self.regolith_transport_parameter + self._grid, linear_diffusivity=self.regolith_transport_parameter ) def _calc_effective_drainage_area(self): @@ -182,21 +187,21 @@ def _calc_effective_drainage_area(self): drainage area, :math:`R_r` is the runoff ratio, and :math:`\alpha` is the saturation parameter. """ - area = self.grid.at_node["drainage_area"] - slope = self.grid.at_node["topographic__steepest_slope"] - cores = self.grid.core_nodes + area = self._grid.at_node["drainage_area"] + slope = self._grid.at_node["topographic__steepest_slope"] + cores = self._grid.core_nodes sat_param = ( self._Kdx - * self.grid.at_node["soil__depth"] - / self.grid.at_node["rainfall__flux"] + * self._grid.at_node["soil__depth"] + / self._grid.at_node["rainfall__flux"] ) eff_area = area[cores] * ( np.exp(-sat_param[cores] * slope[cores] / area[cores]) ) - self.grid.at_node["surface_water__discharge"][cores] = eff_area + self._grid.at_node["surface_water__discharge"][cores] = eff_area def run_one_step(self, step): """Advance model **BasicRtVs** for one time-step of duration step. @@ -245,7 +250,7 @@ def run_one_step(self, step): )[0] # Zero out effective area in flooded nodes - self.grid.at_node["surface_water__discharge"][flooded] = 0.0 + self._grid.at_node["surface_water__discharge"][flooded] = 0.0 # Update the erodibility field self._update_erodibility_field() diff --git a/terrainbento/derived_models/model_basicSa.py b/terrainbento/derived_models/model_basicSa.py index 60b60ade..0d652ebc 100644 --- a/terrainbento/derived_models/model_basicSa.py +++ b/terrainbento/derived_models/model_basicSa.py @@ -70,7 +70,11 @@ class BasicSa(ErosionModel): - ``soil__depth`` """ - _required_fields = ["topographic__elevation", "soil__depth"] + _name = "BasicSa" + + _input_var_names = ("topographic__elevation", "soil__depth") + + _output_var_names = ("topographic__elevation", "soil__depth") def __init__( self, @@ -138,7 +142,7 @@ def __init__( will just run it one step. >>> model.run_one_step(1.) - >>> model.model_time + >>> model.clock.time 1.0 """ @@ -147,7 +151,7 @@ def __init__( super(BasicSa, self).__init__(clock, grid, **kwargs) # verify correct fields are present. - self._verify_fields(self._required_fields) + self._verify_fields(self._input_var_names) # Get Parameters and convert units if necessary: self.m = m_sp @@ -156,26 +160,26 @@ def __init__( # Instantiate a FastscapeEroder component self.eroder = FastscapeEroder( - self.grid, + self._grid, K_sp=self.K, m_sp=self.m, n_sp=self.n, discharge_name="surface_water__discharge", ) - soil_thickness = self.grid.at_node["soil__depth"] - bedrock_elev = self.grid.add_zeros("node", "bedrock__elevation") + soil_thickness = self._grid.at_node["soil__depth"] + bedrock_elev = self._grid.add_zeros("node", "bedrock__elevation") bedrock_elev[:] = self.z - soil_thickness # Instantiate diffusion and weathering components self.diffuser = DepthDependentDiffuser( - self.grid, + self._grid, linear_diffusivity=regolith_transport_parameter, soil_transport_decay_depth=soil_transport_decay_depth, ) self.weatherer = ExponentialWeatherer( - self.grid, + self._grid, soil_production__maximum_rate=soil_production__maximum_rate, soil_production__decay_depth=soil_production__decay_depth, ) @@ -235,8 +239,8 @@ def run_one_step(self, step): # into bedrock has occurred, the bedrock elevation will be higher than # the actual elevation, so we simply re-set bedrock elevation to the # lower of itself or the current elevation. - b = self.grid.at_node["bedrock__elevation"] - b[:] = np.minimum(b, self.grid.at_node["topographic__elevation"]) + b = self._grid.at_node["bedrock__elevation"] + b[:] = np.minimum(b, self._grid.at_node["topographic__elevation"]) # Calculate regolith-production rate self.weatherer.calc_soil_prod_rate() diff --git a/terrainbento/derived_models/model_basicSaVs.py b/terrainbento/derived_models/model_basicSaVs.py index 81a845de..6867bbc1 100644 --- a/terrainbento/derived_models/model_basicSaVs.py +++ b/terrainbento/derived_models/model_basicSaVs.py @@ -71,7 +71,11 @@ class BasicSaVs(ErosionModel): - ``soil__depth`` """ - _required_fields = ["topographic__elevation", "soil__depth"] + _name = "BasicSaVs" + + _input_var_names = ("topographic__elevation", "soil__depth") + + _output_var_names = ("topographic__elevation", "soil__depth") def __init__( self, @@ -142,7 +146,7 @@ def __init__( will just run it one step. >>> model.run_one_step(1.) - >>> model.model_time + >>> model.clock.time 1.0 """ @@ -153,23 +157,23 @@ def __init__( self._ensure_precip_runoff_are_vanilla(vsa_precip=True) # verify correct fields are present. - self._verify_fields(self._required_fields) + self._verify_fields(self._input_var_names) # Get Parameters and convert units if necessary: self.m = m_sp self.n = n_sp self.K = water_erodibility - soil_thickness = self.grid.at_node["soil__depth"] - bedrock_elev = self.grid.add_zeros("node", "bedrock__elevation") + soil_thickness = self._grid.at_node["soil__depth"] + bedrock_elev = self._grid.add_zeros("node", "bedrock__elevation") bedrock_elev[:] = self.z - soil_thickness # Get the effective-area parameter - self._Kdx = hydraulic_conductivity * self.grid.dx + self._Kdx = hydraulic_conductivity * self._grid.dx # Instantiate a FastscapeEroder component self.eroder = FastscapeEroder( - self.grid, + self._grid, discharge_name="surface_water__discharge", K_sp=self.K, m_sp=self.m, @@ -178,34 +182,34 @@ def __init__( # Instantiate a DepthDependentDiffuser component self.diffuser = DepthDependentDiffuser( - self.grid, + self._grid, linear_diffusivity=regolith_transport_parameter, soil_transport_decay_depth=soil_transport_decay_depth, ) self.weatherer = ExponentialWeatherer( - self.grid, + self._grid, soil_production__maximum_rate=soil_production__maximum_rate, soil_production__decay_depth=soil_production__decay_depth, ) def _calc_effective_drainage_area(self): """Calculate and store effective drainage area.""" - area = self.grid.at_node["drainage_area"] - slope = self.grid.at_node["topographic__steepest_slope"] - cores = self.grid.core_nodes + area = self._grid.at_node["drainage_area"] + slope = self._grid.at_node["topographic__steepest_slope"] + cores = self._grid.core_nodes sat_param = ( self._Kdx - * self.grid.at_node["soil__depth"] - / self.grid.at_node["rainfall__flux"] + * self._grid.at_node["soil__depth"] + / self._grid.at_node["rainfall__flux"] ) eff_area = area[cores] * ( np.exp(-sat_param[cores] * slope[cores] / area[cores]) ) - self.grid.at_node["surface_water__discharge"][cores] = eff_area + self._grid.at_node["surface_water__discharge"][cores] = eff_area def run_one_step(self, step): """Advance model **BasicVs** for one time-step of duration step. @@ -250,7 +254,7 @@ def run_one_step(self, step): )[0] # Zero out effective area in flooded nodes - self.grid.at_node["surface_water__discharge"][flooded] = 0.0 + self._grid.at_node["surface_water__discharge"][flooded] = 0.0 # Do some erosion (but not on the flooded nodes) # (if we're varying K through time, update that first) @@ -267,8 +271,8 @@ def run_one_step(self, step): # into bedrock has occurred, the bedrock elevation will be higher than # the actual elevation, so we simply re-set bedrock elevation to the # lower of itself or the current elevation. - b = self.grid.at_node["bedrock__elevation"] - b[:] = np.minimum(b, self.grid.at_node["topographic__elevation"]) + b = self._grid.at_node["bedrock__elevation"] + b[:] = np.minimum(b, self._grid.at_node["topographic__elevation"]) # Calculate regolith-production rate self.weatherer.calc_soil_prod_rate() diff --git a/terrainbento/derived_models/model_basicSt.py b/terrainbento/derived_models/model_basicSt.py index d895ce29..d852807d 100644 --- a/terrainbento/derived_models/model_basicSt.py +++ b/terrainbento/derived_models/model_basicSt.py @@ -45,7 +45,11 @@ class BasicSt(StochasticErosionModel): - ``topographic__elevation`` """ - _required_fields = ["topographic__elevation"] + _name = "BasicSt" + + _input_var_names = ("topographic__elevation",) + + _output_var_names = ("topographic__elevation",) def __init__( self, @@ -105,7 +109,7 @@ def __init__( will just run it one step. >>> model.run_one_step(1.) - >>> model.model_time + >>> model.clock.time 1.0 """ @@ -113,7 +117,7 @@ def __init__( super(BasicSt, self).__init__(clock, grid, **kwargs) # verify correct fields are present. - self._verify_fields(self._required_fields) + self._verify_fields(self._input_var_names) # Get Parameters: self.m = m_sp @@ -129,7 +133,7 @@ def __init__( # Instantiate a FastscapeEroder component self.eroder = FastscapeEroder( - self.grid, + self._grid, K_sp=self.K, m_sp=self.m, n_sp=self.n, @@ -138,7 +142,7 @@ def __init__( # Instantiate a LinearDiffuser component self.diffuser = LinearDiffuser( - self.grid, linear_diffusivity=regolith_transport_parameter + self._grid, linear_diffusivity=regolith_transport_parameter ) def run_one_step(self, step): diff --git a/terrainbento/derived_models/model_basicStTh.py b/terrainbento/derived_models/model_basicStTh.py index 47940db0..36159010 100644 --- a/terrainbento/derived_models/model_basicStTh.py +++ b/terrainbento/derived_models/model_basicStTh.py @@ -49,7 +49,11 @@ class BasicStTh(StochasticErosionModel): - ``topographic__elevation`` """ - _required_fields = ["topographic__elevation"] + _name = "BasicStTh" + + _input_var_names = ("topographic__elevation",) + + _output_var_names = ("topographic__elevation",) def __init__( self, @@ -113,7 +117,7 @@ def __init__( will just run it one step. >>> model.run_one_step(1.) - >>> model.model_time + >>> model.clock.time 1.0 """ @@ -121,7 +125,7 @@ def __init__( super(BasicStTh, self).__init__(clock, grid, **kwargs) # verify correct fields are present. - self._verify_fields(self._required_fields) + self._verify_fields(self._input_var_names) # Get Parameters: self.m = m_sp @@ -140,7 +144,7 @@ def __init__( # Instantiate a FastscapeEroder component self.eroder = StreamPowerSmoothThresholdEroder( - self.grid, + self._grid, K_sp=self.K, m_sp=self.m, n_sp=self.n, @@ -150,7 +154,7 @@ def __init__( # Instantiate a LinearDiffuser component self.diffuser = LinearDiffuser( - self.grid, linear_diffusivity=regolith_transport_parameter + self._grid, linear_diffusivity=regolith_transport_parameter ) def run_one_step(self, step): diff --git a/terrainbento/derived_models/model_basicStVs.py b/terrainbento/derived_models/model_basicStVs.py index 079f4493..8fc8bb98 100644 --- a/terrainbento/derived_models/model_basicStVs.py +++ b/terrainbento/derived_models/model_basicStVs.py @@ -60,7 +60,11 @@ class BasicStVs(StochasticErosionModel): - ``topographic__elevation`` """ - _required_fields = ["topographic__elevation"] + _name = "BasicStVs" + + _input_var_names = ("topographic__elevation",) + + _output_var_names = ("topographic__elevation",) def __init__( self, @@ -123,7 +127,7 @@ def __init__( will just run it one step. >>> model.run_one_step(1.) - >>> model.model_time + >>> model.clock.time 1.0 """ @@ -131,20 +135,20 @@ def __init__( super(BasicStVs, self).__init__(clock, grid, **kwargs) # verify correct fields are present. - self._verify_fields(self._required_fields) + self._verify_fields(self._input_var_names) # Get Parameters: self.m = m_sp self.n = n_sp self.K = water_erodibility - soil_thickness = self.grid.at_node["soil__depth"] + soil_thickness = self._grid.at_node["soil__depth"] # instantiate rain generator self.instantiate_rain_generator() # Add a field for subsurface discharge - self.qss = self.grid.add_zeros("node", "subsurface_water__discharge") + self.qss = self._grid.add_zeros("node", "subsurface_water__discharge") # Get the transmissivity parameter # transmissivity is hydraulic condiuctivity times soil thickness @@ -153,14 +157,14 @@ def __init__( if np.any(self.trans) <= 0.0: raise ValueError("BasicStVs: Transmissivity must be > 0") - self.tlam = self.trans * self.grid._dx # assumes raster + self.tlam = self.trans * self._grid._dx # assumes raster # Run flow routing and lake filler self.flow_accumulator.run_one_step() # Instantiate a FastscapeEroder component self.eroder = FastscapeEroder( - self.grid, + self._grid, K_sp=self.K, m_sp=self.m, n_sp=self.m, @@ -169,24 +173,24 @@ def __init__( # Instantiate a LinearDiffuser component self.diffuser = LinearDiffuser( - self.grid, linear_diffusivity=regolith_transport_parameter + self._grid, linear_diffusivity=regolith_transport_parameter ) def calc_runoff_and_discharge(self): """Calculate runoff rate and discharge; return runoff.""" # Here"s the total (surface + subsurface) discharge - pa = self.rain_rate * self.grid.at_node["drainage_area"] + pa = self.rain_rate * self._grid.at_node["drainage_area"] # slope > 0 active_nodes = np.where( - self.grid.at_node["topographic__steepest_slope"] > 0.0 + self._grid.at_node["topographic__steepest_slope"] > 0.0 )[0] # Transmissivity x lambda x slope = subsurface discharge capacity tls = ( self.tlam[active_nodes] - * self.grid.at_node["topographic__steepest_slope"][active_nodes] + * self._grid.at_node["topographic__steepest_slope"][active_nodes] ) # Subsurface discharge: zero where slope is flat @@ -197,9 +201,9 @@ def calc_runoff_and_discharge(self): # # Note that roundoff errors can sometimes produce a tiny negative # value when qss and pa are close; make sure these are set to 0 - self.grid.at_node["surface_water__discharge"][:] = pa - self.qss - self.grid.at_node["surface_water__discharge"][ - self.grid.at_node["surface_water__discharge"] < 0.0 + self._grid.at_node["surface_water__discharge"][:] = pa - self.qss + self._grid.at_node["surface_water__discharge"][ + self._grid.at_node["surface_water__discharge"] < 0.0 ] = 0.0 return np.nan diff --git a/terrainbento/derived_models/model_basicTh.py b/terrainbento/derived_models/model_basicTh.py index e0bd52e7..143ce481 100644 --- a/terrainbento/derived_models/model_basicTh.py +++ b/terrainbento/derived_models/model_basicTh.py @@ -44,7 +44,11 @@ class BasicTh(ErosionModel): - ``topographic__elevation`` """ - _required_fields = ["topographic__elevation"] + _name = "BasicTh" + + _input_var_names = ("topographic__elevation", "water__unit_flux_in") + + _output_var_names = ("topographic__elevation",) def __init__( self, @@ -106,7 +110,7 @@ def __init__( will just run it one step. >>> model.run_one_step(1.) - >>> model.model_time + >>> model.clock.time 1.0 """ @@ -114,7 +118,7 @@ def __init__( super(BasicTh, self).__init__(clock, grid, **kwargs) # verify correct fields are present. - self._verify_fields(self._required_fields) + self._verify_fields(self._input_var_names) # Get Parameters and convert units if necessary: self.m = m_sp @@ -126,7 +130,7 @@ def __init__( # Instantiate a FastscapeEroder component self.eroder = StreamPowerSmoothThresholdEroder( - self.grid, + self._grid, K_sp=self.K, m_sp=self.m, n_sp=self.n, @@ -136,7 +140,7 @@ def __init__( # Instantiate a LinearDiffuser component self.diffuser = LinearDiffuser( - self.grid, linear_diffusivity=regolith_transport_parameter + self._grid, linear_diffusivity=regolith_transport_parameter ) def run_one_step(self, step): diff --git a/terrainbento/derived_models/model_basicThVs.py b/terrainbento/derived_models/model_basicThVs.py index af7fe674..bd655145 100644 --- a/terrainbento/derived_models/model_basicThVs.py +++ b/terrainbento/derived_models/model_basicThVs.py @@ -55,7 +55,11 @@ class BasicThVs(ErosionModel): - ``soil__depth`` """ - _required_fields = ["topographic__elevation", "soil__depth"] + _name = "BasicThVs" + + _input_var_names = ("topographic__elevation", "soil__depth") + + _output_var_names = ("topographic__elevation",) def __init__( self, @@ -121,7 +125,7 @@ def __init__( will just run it one step. >>> model.run_one_step(1.) - >>> model.model_time + >>> model.clock.time 1.0 """ @@ -132,7 +136,7 @@ def __init__( self._ensure_precip_runoff_are_vanilla(vsa_precip=True) # verify correct fields are present. - self._verify_fields(self._required_fields) + self._verify_fields(self._input_var_names) self.m = m_sp self.n = n_sp @@ -142,11 +146,11 @@ def __init__( raise ValueError("Model only supports n = 1.") # Get the effective-area parameter - self._Kdx = hydraulic_conductivity * self.grid.dx + self._Kdx = hydraulic_conductivity * self._grid.dx # Instantiate a FastscapeEroder component self.eroder = StreamPowerSmoothThresholdEroder( - self.grid, + self._grid, K_sp=self.K, m_sp=self.m, n_sp=self.n, @@ -156,27 +160,27 @@ def __init__( # Instantiate a LinearDiffuser component self.diffuser = LinearDiffuser( - self.grid, linear_diffusivity=regolith_transport_parameter + self._grid, linear_diffusivity=regolith_transport_parameter ) def _calc_effective_drainage_area(self): """Calculate and store effective drainage area.""" - area = self.grid.at_node["drainage_area"] - slope = self.grid.at_node["topographic__steepest_slope"] - cores = self.grid.core_nodes + area = self._grid.at_node["drainage_area"] + slope = self._grid.at_node["topographic__steepest_slope"] + cores = self._grid.core_nodes sat_param = ( self._Kdx - * self.grid.at_node["soil__depth"] - / self.grid.at_node["rainfall__flux"] + * self._grid.at_node["soil__depth"] + / self._grid.at_node["rainfall__flux"] ) eff_area = area[cores] * ( np.exp(-sat_param[cores] * slope[cores] / area[cores]) ) - self.grid.at_node["surface_water__discharge"][cores] = eff_area + self._grid.at_node["surface_water__discharge"][cores] = eff_area def run_one_step(self, step): """Advance model **BasicThVs** for one time-step of duration step. @@ -221,7 +225,7 @@ def run_one_step(self, step): )[0] # Zero out effective area in flooded nodes - self.grid.at_node["surface_water__discharge"][flooded] = 0.0 + self._grid.at_node["surface_water__discharge"][flooded] = 0.0 # Do some erosion (but not on the flooded nodes) # (if we're varying K through time, update that first) diff --git a/terrainbento/derived_models/model_basicVs.py b/terrainbento/derived_models/model_basicVs.py index 4f83c680..72d249db 100644 --- a/terrainbento/derived_models/model_basicVs.py +++ b/terrainbento/derived_models/model_basicVs.py @@ -52,7 +52,11 @@ class BasicVs(ErosionModel): - ``soil__depth`` """ - _required_fields = ["topographic__elevation", "soil__depth"] + _name = "BasicVs" + + _input_var_names = ("topographic__elevation", "soil__depth") + + _output_var_names = ("topographic__elevation",) def __init__( self, @@ -114,7 +118,7 @@ def __init__( will just run it one step. >>> model.run_one_step(1.) - >>> model.model_time + >>> model.clock.time 1.0 """ @@ -126,7 +130,7 @@ def __init__( self._ensure_precip_runoff_are_vanilla(vsa_precip=True) # verify correct fields are present. - self._verify_fields(self._required_fields) + self._verify_fields(self._input_var_names) # Get Parameters: self.m = m_sp @@ -134,16 +138,16 @@ def __init__( self.K = water_erodibility # Add a field for effective drainage area - self.grid.at_node["surface_water__discharge"] = self.grid.add_zeros( + self._grid.at_node["surface_water__discharge"] = self._grid.add_zeros( "node", "effective_drainage_area" ) # Get the effective-area parameter - self._Kdx = hydraulic_conductivity * self.grid.dx + self._Kdx = hydraulic_conductivity * self._grid.dx # Instantiate a FastscapeEroder component self.eroder = FastscapeEroder( - self.grid, + self._grid, discharge_name="surface_water__discharge", K_sp=self.K, m_sp=self.m, @@ -152,26 +156,26 @@ def __init__( # Instantiate a LinearDiffuser component self.diffuser = LinearDiffuser( - self.grid, linear_diffusivity=regolith_transport_parameter + self._grid, linear_diffusivity=regolith_transport_parameter ) def _calc_effective_drainage_area(self): """Calculate and store effective drainage area.""" - area = self.grid.at_node["drainage_area"] - slope = self.grid.at_node["topographic__steepest_slope"] - cores = self.grid.core_nodes + area = self._grid.at_node["drainage_area"] + slope = self._grid.at_node["topographic__steepest_slope"] + cores = self._grid.core_nodes sat_param = ( self._Kdx - * self.grid.at_node["soil__depth"] - / self.grid.at_node["rainfall__flux"] + * self._grid.at_node["soil__depth"] + / self._grid.at_node["rainfall__flux"] ) eff_area = area[cores] * ( np.exp(-sat_param[cores] * slope[cores] / area[cores]) ) - self.grid.at_node["surface_water__discharge"][cores] = eff_area + self._grid.at_node["surface_water__discharge"][cores] = eff_area def run_one_step(self, step): """Advance model **BasicVs** for one time-step of duration step. @@ -216,7 +220,7 @@ def run_one_step(self, step): )[0] # Zero out effective area in flooded nodes - self.grid.at_node["surface_water__discharge"][flooded] = 0.0 + self._grid.at_node["surface_water__discharge"][flooded] = 0.0 # Do some erosion (but not on the flooded nodes) # (if we're varying K through time, update that first) diff --git a/terrainbento/model_template/model_template.py b/terrainbento/model_template/model_template.py index f8787f1c..93b00c30 100644 --- a/terrainbento/model_template/model_template.py +++ b/terrainbento/model_template/model_template.py @@ -33,7 +33,9 @@ class ModelTemplate(ErosionModel): # The model must inherit from either 3. Description of the governing equation of the model. """ - _required_fields = ["topographic__elevation"] + _input_var_names = ("topographic__elevation", "water__unit_flux_in") + + _output_var_names = ("topographic__elevation",) def __init__( self, @@ -88,7 +90,7 @@ def __init__( # Do not change any additional parts of this line. # verify correct fields are present. - self._verify_fields(self._required_fields) + self._verify_fields(self._input_var_names) # put all actions needed to initialize the model below this line. diff --git a/terrainbento/precipitators/__init__.py b/terrainbento/precipitators/__init__.py index 9f293e9d..38e9160c 100644 --- a/terrainbento/precipitators/__init__.py +++ b/terrainbento/precipitators/__init__.py @@ -1,6 +1,6 @@ """Precipitators in the terrainbento package.""" -from .uniform_precipitation import UniformPrecipitator from .random_precipitation import RandomPrecipitator +from .uniform_precipitation import UniformPrecipitator __all__ = ["UniformPrecipitator", "RandomPrecipitator"] diff --git a/terrainbento/precipitators/random_precipitation.py b/terrainbento/precipitators/random_precipitation.py index 0abde95e..f3564c17 100644 --- a/terrainbento/precipitators/random_precipitation.py +++ b/terrainbento/precipitators/random_precipitation.py @@ -56,7 +56,7 @@ def __init__(self, grid, distribution="uniform", **kwargs): grid.add_ones("node", "rainfall__flux") self.function = np.random.__dict__[distribution] self._kwargs = kwargs - self.run_one_step(0.) + self.run_one_step(0.0) def run_one_step(self, step): """Run **RandomPrecipitator** forward by duration ``step``""" diff --git a/terrainbento/utilities/__init__.py b/terrainbento/utilities/__init__.py index 2220a581..115e5072 100644 --- a/terrainbento/utilities/__init__.py +++ b/terrainbento/utilities/__init__.py @@ -1,4 +1,3 @@ - from terrainbento.utilities.file_compare import filecmp __all__ = ["filecmp"] diff --git a/tests/conftest.py b/tests/conftest.py index a33f025c..82bd2ace 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -30,7 +30,7 @@ def Kt(): @pytest.fixture() def grid_1(): - grid = RasterModelGrid((3, 21), xy_spacing=100.) + grid = RasterModelGrid((3, 21), xy_spacing=100.0) grid.set_closed_boundaries_at_grid_edges(False, True, False, True) grid.add_zeros("node", "topographic__elevation") grid.add_ones("node", "soil__depth") @@ -40,37 +40,37 @@ def grid_1(): @pytest.fixture() def grid_2(): - grid = RasterModelGrid((8, 20), xy_spacing=100.) + grid = RasterModelGrid((8, 20), xy_spacing=100.0) grid.set_closed_boundaries_at_grid_edges(False, True, False, True) grid.add_zeros("node", "topographic__elevation") grid.add_ones("node", "soil__depth") lith = grid.add_zeros("node", "lithology_contact__elevation") lith[:80] = 10 - lith[80:] = -10000. + lith[80:] = -10000.0 return grid @pytest.fixture() def grid_3(): - grid = RasterModelGrid((21, 3), xy_spacing=100.) + grid = RasterModelGrid((21, 3), xy_spacing=100.0) grid.set_closed_boundaries_at_grid_edges(False, True, False, True) grid.add_zeros("node", "topographic__elevation") grid.add_ones("node", "soil__depth") lith = grid.add_zeros("node", "lithology_contact__elevation") - lith[grid.core_nodes[:9]] = -100000. - lith[grid.core_nodes[9:]] = 100000. + lith[grid.core_nodes[:9]] = -100000.0 + lith[grid.core_nodes[9:]] = 100000.0 return grid @pytest.fixture() def grid_4(): - grid = RasterModelGrid((3, 21), xy_spacing=10.) + grid = RasterModelGrid((3, 21), xy_spacing=10.0) grid.set_closed_boundaries_at_grid_edges(False, True, False, True) grid.add_zeros("node", "topographic__elevation") grid.add_ones("node", "soil__depth") lith = grid.add_zeros("node", "lithology_contact__elevation") - lith[grid.core_nodes[:9]] = -100000. - lith[grid.core_nodes[9:]] = 100000. + lith[grid.core_nodes[:9]] = -100000.0 + lith[grid.core_nodes[9:]] = 100000.0 return grid @@ -81,14 +81,14 @@ def grid_5(): grid.add_ones("node", "soil__depth") lith = grid.add_zeros("node", "lithology_contact__elevation") lith[:27] = -30 - lith[27:] = 10. - lith[grid.boundary_nodes] = -9999. + lith[27:] = 10.0 + lith[grid.boundary_nodes] = -9999.0 return grid @pytest.fixture def almost_default_grid(): - grid = RasterModelGrid((4, 5), xy_spacing=100.) + grid = RasterModelGrid((4, 5), xy_spacing=100.0) grid.add_zeros("node", "topographic__elevation") return grid @@ -103,49 +103,49 @@ def simple_square_grid(): @pytest.fixture() def clock_simple(): - clock_simple = Clock(step=1000., stop=5.1e6) + clock_simple = Clock(step=1000.0, stop=5.1e6) return clock_simple @pytest.fixture() def clock_02(): - clock_02 = Clock.from_dict({"step": 10., "stop": 1000.}) + clock_02 = Clock.from_dict({"step": 10.0, "stop": 1000.0}) return clock_02 @pytest.fixture() def clock_04(): - clock_04 = Clock.from_dict({"step": 10., "stop": 100000.}) + clock_04 = Clock.from_dict({"step": 10.0, "stop": 100000.0}) return clock_04 @pytest.fixture() def clock_05(): - clock_05 = Clock.from_dict({"step": 10., "stop": 200.}) + clock_05 = Clock.from_dict({"step": 10.0, "stop": 200.0}) return clock_05 @pytest.fixture() def clock_06(): - clock_06 = Clock.from_dict({"step": 1., "stop": 3.}) + clock_06 = Clock.from_dict({"step": 1.0, "stop": 3.0}) return clock_06 @pytest.fixture() def clock_07(): - clock_07 = Clock.from_dict({"step": 10., "stop": 10000.}) + clock_07 = Clock.from_dict({"step": 10.0, "stop": 10000.0}) return clock_07 @pytest.fixture() def clock_08(): - clock_08 = Clock(step=1., stop=20.) + clock_08 = Clock(step=1.0, stop=50000) return clock_08 @pytest.fixture() def clock_09(): - clock_09 = Clock(step=2., stop=200.) + clock_09 = Clock(step=2.0, stop=60000.0) return clock_09 diff --git a/tests/test_bad_n_sp.py b/tests/test_bad_n_sp.py index 98aa1a77..61a02c01 100644 --- a/tests/test_bad_n_sp.py +++ b/tests/test_bad_n_sp.py @@ -1,4 +1,3 @@ - import pytest from terrainbento import ( diff --git a/tests/test_base_class_erosion_model_boundary_handers.py b/tests/test_base_class_erosion_model_boundary_handers.py index 4d27ec04..a1dda850 100644 --- a/tests/test_base_class_erosion_model_boundary_handers.py +++ b/tests/test_base_class_erosion_model_boundary_handers.py @@ -53,7 +53,7 @@ def test_single_node_blh_with_closed_boundaries( "boundary_handlers": {"SingleNodeBaselevelHandler": snblh}, } model = Basic(**params) - assert model.grid.status_at_node[3] == FIXED_VALUE_BOUNDARY + assert model._grid.status_at_node[3] == FIXED_VALUE_BOUNDARY def test_pass_two_boundary_handlers(clock_simple, simple_square_grid, U): @@ -76,13 +76,13 @@ def test_pass_two_boundary_handlers(clock_simple, simple_square_grid, U): truth = np.zeros(model.z.size) truth[0] -= U - truth[model.grid.core_nodes] += U + truth[model._grid.core_nodes] += U assert_array_equal(model.z, truth) status_at_node = np.zeros(model.z.size) - status_at_node[model.grid.boundary_nodes] = CLOSED_BOUNDARY + status_at_node[model._grid.boundary_nodes] = CLOSED_BOUNDARY status_at_node[0] = FIXED_VALUE_BOUNDARY - assert_array_equal(model.grid.status_at_node, status_at_node) + assert_array_equal(model._grid.status_at_node, status_at_node) def test_generic_bch(clock_simple, simple_square_grid): @@ -102,15 +102,15 @@ def test_generic_bch(clock_simple, simple_square_grid): # assertion tests assert "GenericFuncBaselevelHandler" in model.boundary_handlers - assert_array_equal(np.where(bh.nodes_to_lower)[0], model.grid.core_nodes) + assert_array_equal(np.where(bh.nodes_to_lower)[0], model._grid.core_nodes) - step = 10. + step = 10.0 model.run_one_step(step) - dzdt = -(model.grid.x_of_node + model.grid.y_of_node) - truth_z = -1. * dzdt * step + dzdt = -(model._grid.x_of_node + model._grid.y_of_node) + truth_z = -1.0 * dzdt * step assert_array_equal( - model.z[model.grid.core_nodes], truth_z[model.grid.core_nodes] + model.z[model._grid.core_nodes], truth_z[model._grid.core_nodes] ) @@ -134,9 +134,9 @@ def test_capture_node(clock_simple, simple_square_grid): # assertion tests assert "CaptureNodeBaselevelHandler" in model.boundary_handlers assert model.z[1] == 0 - model.run_one_step(10.) + model.run_one_step(10.0) assert model.z[1] == 0 model.run_one_step(10) - assert model.z[1] == -30. + assert model.z[1] == -30.0 model.run_one_step(10) - assert model.z[1] == -31. + assert model.z[1] == -31.0 diff --git a/tests/test_base_class_erosion_model_inputs.py b/tests/test_base_class_erosion_model_inputs.py index 6f96f3b5..b9dc389a 100644 --- a/tests/test_base_class_erosion_model_inputs.py +++ b/tests/test_base_class_erosion_model_inputs.py @@ -64,17 +64,17 @@ def test_input_file(tmpdir, inputs_yaml): em = ErosionModel.from_file("./params.yaml") - assert isinstance(em.grid, HexModelGrid) - assert em.grid.number_of_nodes == 56 + assert isinstance(em._grid, HexModelGrid) + assert em._grid.number_of_nodes == 56 for field in at_node_fields: - assert field in em.grid.at_node + assert field in em._grid.at_node assert isinstance(em.flow_accumulator, FlowAccumulator) is True assert em.flow_accumulator.flow_director._name == "FlowDirectorSteepest" assert em.boundary_handlers == {} assert em.output_writers == {} assert em.save_first_timestep is True assert em._out_file_name == "terrainbento_output" - assert em._model_time == 0. + assert em.clock.time == 0.0 def test_parameters(clock_simple): @@ -96,17 +96,17 @@ def test_parameters(clock_simple): } em = ErosionModel.from_dict(params) - assert isinstance(em.grid, HexModelGrid) - assert em.grid.number_of_nodes == 56 + assert isinstance(em._grid, HexModelGrid) + assert em._grid.number_of_nodes == 56 for field in at_node_fields: - assert field in em.grid.at_node + assert field in em._grid.at_node assert isinstance(em.flow_accumulator, FlowAccumulator) is True assert em.flow_accumulator.flow_director._name == "FlowDirectorSteepest" assert em.boundary_handlers == {} assert em.output_writers == {} assert em.save_first_timestep is True assert em._out_file_name == "terrainbento_output" - assert em._model_time == 0. + assert em.clock.time == 0.0 def test_string(tmpdir, inputs_yaml): @@ -118,14 +118,14 @@ def test_string(tmpdir, inputs_yaml): contents = f.read() em = ErosionModel.from_file(contents) - assert isinstance(em.grid, HexModelGrid) - assert em.grid.number_of_nodes == 56 + assert isinstance(em._grid, HexModelGrid) + assert em._grid.number_of_nodes == 56 for field in at_node_fields: - assert field in em.grid.at_node + assert field in em._grid.at_node assert isinstance(em.flow_accumulator, FlowAccumulator) is True assert em.flow_accumulator.flow_director._name == "FlowDirectorSteepest" assert em.boundary_handlers == {} assert em.output_writers == {} assert em.save_first_timestep is True assert em._out_file_name == "terrainbento_output" - assert em._model_time == 0. + assert em.clock.time == 0.0 diff --git a/tests/test_base_class_erosion_model_misc.py b/tests/test_base_class_erosion_model_misc.py index 18251733..75fb9d06 100644 --- a/tests/test_base_class_erosion_model_misc.py +++ b/tests/test_base_class_erosion_model_misc.py @@ -8,10 +8,10 @@ def test_calc_cumulative_erosion(clock_simple, grid_1): em = ErosionModel(grid=grid_1, clock=clock_simple) - assert np.array_equiv(em.z, 0.) is True - em.z += 1. + assert np.array_equiv(em.z, 0.0) is True + em.z += 1.0 em.calculate_cumulative_change() assert ( - np.array_equiv(em.grid.at_node["cumulative_elevation_change"], 1.) + np.array_equiv(em._grid.at_node["cumulative_elevation_change"], 1.0) is True ) diff --git a/tests/test_base_class_erosion_model_output_writers.py b/tests/test_base_class_erosion_model_output_writers.py index 91644428..891a57bb 100644 --- a/tests/test_base_class_erosion_model_output_writers.py +++ b/tests/test_base_class_erosion_model_output_writers.py @@ -13,26 +13,26 @@ def output_writer_function_a(model): - average_elevation = np.mean(model.z[model.grid.core_nodes]) - with open("ow_func_a." + str(model.model_time) + ".txt", "w") as f: + average_elevation = np.mean(model.z[model._grid.core_nodes]) + with open("ow_func_a." + str(model.clock.time) + ".txt", "w") as f: f.write(str(average_elevation)) def output_writer_function_b(model): - minimum_elevation = np.min(model.z[model.grid.core_nodes]) - with open("ow_func_b." + str(model.model_time) + ".txt", "w") as f: + minimum_elevation = np.min(model.z[model._grid.core_nodes]) + with open("ow_func_b." + str(model.clock.time) + ".txt", "w") as f: f.write(str(minimum_elevation)) class output_writer_class_a(object): def __init__(self, model): self.model = model - self.change = model.grid.at_node["cumulative_elevation_change"] + self.change = model._grid.at_node["cumulative_elevation_change"] def run_one_step(self): - average_change = np.mean(self.change[self.model.grid.core_nodes]) + average_change = np.mean(self.change[self.model._grid.core_nodes]) with open( - "ow_class_a." + str(self.model.model_time) + ".txt", "w" + "ow_class_a." + str(self.model.clock.time) + ".txt", "w" ) as f: f.write(str(average_change)) @@ -40,12 +40,12 @@ def run_one_step(self): class output_writer_class_b(object): def __init__(self, model): self.model = model - self.change = model.grid.at_node["cumulative_elevation_change"] + self.change = model._grid.at_node["cumulative_elevation_change"] def run_one_step(self): - min_change = np.min(self.change[self.model.grid.core_nodes]) + min_change = np.min(self.change[self.model._grid.core_nodes]) with open( - "ow_class_b." + str(self.model.model_time) + ".txt", "w" + "ow_class_b." + str(self.model.clock.time) + ".txt", "w" ) as f: f.write(str(min_change)) diff --git a/tests/test_base_class_erosion_model_run.py b/tests/test_base_class_erosion_model_run.py index b0992aca..9f2f8c31 100644 --- a/tests/test_base_class_erosion_model_run.py +++ b/tests/test_base_class_erosion_model_run.py @@ -10,8 +10,8 @@ def test_run_for(tmpdir, basic_inputs_yaml): fp.write(basic_inputs_yaml) model = Basic.from_file("./params.yaml") model._out_file_name = "run_for_output" - model.run_for(10., 100.) - assert model.model_time == 100. + model.run_for(10.0, 100.0) + assert model.clock.time == 100.0 def test_finalize(tmpdir, basic_inputs_yaml): @@ -19,7 +19,7 @@ def test_finalize(tmpdir, basic_inputs_yaml): with open("params.yaml", "w") as fp: fp.write(basic_inputs_yaml) model = Basic.from_file("./params.yaml") - model.finalize() + model.finalize_terrainbento_run() def test_run(tmpdir, basic_inputs_yaml): @@ -29,5 +29,5 @@ def test_run(tmpdir, basic_inputs_yaml): model = Basic.from_file("./params.yaml") model._out_file_name = "run_output" model.run() - assert model.model_time == 200. + assert model.clock.time == 200.0 model.remove_output_netcdfs() diff --git a/tests/test_base_class_stochastic_erosion_model.py b/tests/test_base_class_stochastic_erosion_model.py index 75624df5..ebd032a4 100644 --- a/tests/test_base_class_stochastic_erosion_model.py +++ b/tests/test_base_class_stochastic_erosion_model.py @@ -52,9 +52,9 @@ def test_run_stochastic_opt_true(clock_04, grid_1): "water_erodibility": 0.01, "regolith_transport_parameter": 0.1, "infiltration_capacity": 0.0, - "mean_storm_duration": 2., - "mean_interstorm_duration": 3., - "mean_storm_depth": 1., + "mean_storm_duration": 2.0, + "mean_interstorm_duration": 3.0, + "mean_storm_depth": 1.0, "random_seed": 1234, } @@ -101,7 +101,7 @@ def test_run_stochastic_opt_false(clock_05, grid_1): "water_erodibility": 0.01, "regolith_transport_parameter": 0.1, "infiltration_capacity": 0.0, - "rainfall__mean_rate": 1., + "rainfall__mean_rate": 1.0, "rainfall_intermittency_factor": 0.1, "rainfall__shape_factor": 0.6, "number_of_sub_time_steps": 1, @@ -109,8 +109,9 @@ def test_run_stochastic_opt_false(clock_05, grid_1): } model = BasicSt(**params) + model.clock.stop = 20000.0 assert model.opt_stochastic_duration is False - model.run_for(model.clock.step, 10000.) + model.run_for(model.clock.step, 10000.0) rainfall_rate = np.asarray(model.rain_record["rainfall_rate"]) event_duration = np.asarray(model.rain_record["event_duration"]) @@ -121,7 +122,7 @@ def test_run_stochastic_opt_false(clock_05, grid_1): assert ( np.array_equiv( dry_times, - model.clock.step * (1. - params["rainfall_intermittency_factor"]), + model.clock.step * (1.0 - params["rainfall_intermittency_factor"]), ) is True ) @@ -151,9 +152,9 @@ def test_reset_random_seed_stochastic_duration_true(clock_simple, grid_1): "water_erodibility": 0.01, "regolith_transport_parameter": 0.1, "infiltration_capacity": 0.0, - "mean_storm_duration": 2., - "mean_interstorm_duration": 3., - "mean_storm_depth": 1., + "mean_storm_duration": 2.0, + "mean_interstorm_duration": 3.0, + "mean_storm_depth": 1.0, "random_seed": 0, } @@ -203,7 +204,7 @@ def test_reset_random_seed_stochastic_duration_false(clock_05, grid_1): "water_erodibility": 0.01, "regolith_transport_parameter": 0.1, "infiltration_capacity": 0.0, - "rainfall__mean_rate": 1., + "rainfall__mean_rate": 1.0, "rainfall_intermittency_factor": 0.1, "rainfall__shape_factor": 0.6, "number_of_sub_time_steps": 1, @@ -242,7 +243,7 @@ def test_float_number_of_sub_time_steps(clock_05, grid_1): "water_erodibility": 0.01, "regolith_transport_parameter": 0.1, "infiltration_capacity": 0.0, - "rainfall__mean_rate": 1., + "rainfall__mean_rate": 1.0, "rainfall_intermittency_factor": 0.1, "rainfall__shape_factor": 0.6, "number_of_sub_time_steps": 1.5, @@ -264,7 +265,7 @@ def test_run_opt_false_with_changer(clock_06, grid_1, precip_defaults): "water_erodibility": 0.01, "regolith_transport_parameter": 0.1, "infiltration_capacity": 0.0, - "rainfall__mean_rate": 1., + "rainfall__mean_rate": 1.0, "rainfall_intermittency_factor": 0.5, "rainfall__shape_factor": 0.65, "number_of_sub_time_steps": 1, @@ -317,7 +318,7 @@ def test_not_specifying_record_rain(clock_05, grid_1): "water_erodibility": 0.01, "regolith_transport_parameter": 0.1, "infiltration_capacity": 0.0, - "rainfall__mean_rate": 1., + "rainfall__mean_rate": 1.0, "rainfall_intermittency_factor": 0.1, "rainfall__shape_factor": 0.6, "number_of_sub_time_steps": 1, @@ -345,7 +346,7 @@ def test_finalize_opt_duration_stochastic_false_too_short(clock_05, grid_1): "water_erodibility": 0.01, "regolith_transport_parameter": 0.1, "infiltration_capacity": 0.0, - "rainfall__mean_rate": 1., + "rainfall__mean_rate": 1.0, "rainfall_intermittency_factor": 0.1, "rainfall__shape_factor": 0.6, "number_of_sub_time_steps": 1, @@ -356,7 +357,7 @@ def test_finalize_opt_duration_stochastic_false_too_short(clock_05, grid_1): model.reset_random_seed() model.run_for(model.clock.step, model.clock.stop) with pytest.raises(RuntimeError): - model.finalize() + model.finalize_terrainbento_run() os.remove("storm_sequence.txt") @@ -372,7 +373,7 @@ def test_finalize_opt_duration_stochastic_false_no_rain(clock_07, grid_1): "water_erodibility": 0.01, "regolith_transport_parameter": 0.1, "infiltration_capacity": 0.0, - "rainfall__mean_rate": 1., + "rainfall__mean_rate": 1.0, "rainfall_intermittency_factor": 0.0, "rainfall__shape_factor": 0.6, "number_of_sub_time_steps": 1, @@ -382,7 +383,7 @@ def test_finalize_opt_duration_stochastic_false_no_rain(clock_07, grid_1): model.reset_random_seed() model.run_for(model.clock.step, model.clock.stop) with pytest.raises(ValueError): - model.finalize() + model.finalize_terrainbento_run() def test_finalize_opt_duration_stochastic_false(clock_07, grid_1): @@ -396,7 +397,7 @@ def test_finalize_opt_duration_stochastic_false(clock_07, grid_1): "water_erodibility": 0.01, "regolith_transport_parameter": 0.1, "infiltration_capacity": 0.0, - "rainfall__mean_rate": 1., + "rainfall__mean_rate": 1.0, "rainfall_intermittency_factor": 0.1, "rainfall__shape_factor": 0.6, "number_of_sub_time_steps": 1, @@ -405,7 +406,7 @@ def test_finalize_opt_duration_stochastic_false(clock_07, grid_1): model = BasicSt(**params) model.reset_random_seed() model.run_for(model.clock.step, model.clock.stop) - model.finalize() + model.finalize_terrainbento_run() # assert that these are correct truth_file = os.path.join( @@ -433,16 +434,16 @@ def test_finalize_opt_duration_stochastic_true(clock_07, grid_1): "water_erodibility": 0.01, "regolith_transport_parameter": 0.1, "infiltration_capacity": 0.0, - "mean_storm_duration": 2., - "mean_interstorm_duration": 3., - "mean_storm_depth": 1., + "mean_storm_duration": 2.0, + "mean_interstorm_duration": 3.0, + "mean_storm_depth": 1.0, "random_seed": 1234, } model = BasicSt(**params) model.reset_random_seed() model.run_for(model.clock.step, model.clock.stop) - model.finalize() + model.finalize_terrainbento_run() # assert that these are correct truth_file = os.path.join( @@ -463,14 +464,14 @@ def test_runoff_equals_zero(clock_07, grid_1): "n_sp": 1.0, "water_erodibility": 0.01, "regolith_transport_parameter": 0.1, - "infiltration_capacity": 100000., + "infiltration_capacity": 100000.0, "rainfall__mean_rate": 0.0, "rainfall_intermittency_factor": 0.1, - "rainfall__shape_factor": 1., + "rainfall__shape_factor": 1.0, "number_of_sub_time_steps": 1, "random_seed": 1234, } model = BasicSt(**params) - model.run_one_step(1.) + model.run_one_step(1.0) runoff = model.calc_runoff_and_discharge() assert runoff == 0 diff --git a/tests/test_base_class_two_lithology_base.py b/tests/test_base_class_two_lithology_base.py index a04f1123..f36d77a0 100644 --- a/tests/test_base_class_two_lithology_base.py +++ b/tests/test_base_class_two_lithology_base.py @@ -14,8 +14,8 @@ def test_no_contact_zone_width(clock_simple, grid_5): model = TwoLithologyErosionModel(**params) model._setup_rock_and_till() - truth = np.ones(model.grid.size("node")) - truth[model.grid.core_nodes[14:]] = 0.0 + truth = np.ones(model._grid.size("node")) + truth[model._grid.core_nodes[14:]] = 0.0 assert_array_equal(model.erody_wt, truth) @@ -24,7 +24,7 @@ def test_contact_zone_width(clock_simple, grid_5): params = { "grid": grid_5, "clock": clock_simple, - "contact_zone__width": 10., + "contact_zone__width": 10.0, } model = TwoLithologyErosionModel(**params) @@ -32,16 +32,16 @@ def test_contact_zone_width(clock_simple, grid_5): truth = np.array( [ - 0., - 0., - 0., - 0., - 0., - 0., - 0., - 0., - 0., - 0., + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, 0.95257413, 0.95257413, 0.95257413, @@ -49,8 +49,8 @@ def test_contact_zone_width(clock_simple, grid_5): 0.95257413, 0.95257413, 0.95257413, - 0., - 0., + 0.0, + 0.0, 0.95257413, 0.95257413, 0.95257413, @@ -58,8 +58,8 @@ def test_contact_zone_width(clock_simple, grid_5): 0.95257413, 0.95257413, 0.95257413, - 0., - 0., + 0.0, + 0.0, 0.26894142, 0.26894142, 0.26894142, @@ -67,8 +67,8 @@ def test_contact_zone_width(clock_simple, grid_5): 0.26894142, 0.26894142, 0.26894142, - 0., - 0., + 0.0, + 0.0, 0.26894142, 0.26894142, 0.26894142, @@ -76,16 +76,16 @@ def test_contact_zone_width(clock_simple, grid_5): 0.26894142, 0.26894142, 0.26894142, - 0., - 0., - 0., - 0., - 0., - 0., - 0., - 0., - 0., - 0., + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, ] ) diff --git a/tests/test_base_level_not_core_node_baselevel_handler.py b/tests/test_base_level_not_core_node_baselevel_handler.py index 40544bdb..d46bea8f 100644 --- a/tests/test_base_level_not_core_node_baselevel_handler.py +++ b/tests/test_base_level_not_core_node_baselevel_handler.py @@ -28,7 +28,7 @@ def test_hex(): not_closed = mg.status_at_node == 0 # closed should have been downdropped 10*0.1 - assert_array_equal(z[closed], -1. * np.ones(np.sum(closed))) + assert_array_equal(z[closed], -1.0 * np.ones(np.sum(closed))) # not closed should have stayed the same assert_array_equal(z[not_closed], np.zeros(np.sum(not_closed))) diff --git a/tests/test_clock.py b/tests/test_clock.py index 9164fec3..9b1eda80 100644 --- a/tests/test_clock.py +++ b/tests/test_clock.py @@ -11,9 +11,9 @@ def test_from_file(tmpdir, clock_yaml): with open("params.yaml", "w") as fp: fp.write(clock_yaml) clock = Clock.from_file("./params.yaml") - assert clock.start == 1. - assert clock.stop == 11. - assert clock.step == 2. + assert clock.start == 1.0 + assert clock.stop == 11.0 + assert clock.step == 2.0 @pytest.mark.parametrize("bad", ["start", "stop", "step"]) diff --git a/tests/test_create_precipitator.py b/tests/test_create_precipitator.py index 527490d5..915c1421 100644 --- a/tests/test_create_precipitator.py +++ b/tests/test_create_precipitator.py @@ -1,4 +1,3 @@ - import pytest from terrainbento import ErosionModel diff --git a/tests/test_linear_diffusion.py b/tests/test_linear_diffusion.py index 74a00712..aff8be8d 100644 --- a/tests/test_linear_diffusion.py +++ b/tests/test_linear_diffusion.py @@ -52,20 +52,20 @@ def test_stochastic_linear_diffusion(clock_simple, grid_1, U, Model): for _ in range(nts): model.run_one_step(1000) reference_node = 9 - predicted_z = model.z[model.grid.core_nodes[reference_node]] - ( - U / (2. * params["regolith_transport_parameter"]) + predicted_z = model.z[model._grid.core_nodes[reference_node]] - ( + U / (2.0 * params["regolith_transport_parameter"]) ) * ( ( - model.grid.x_of_node - - model.grid.x_of_node[model.grid.core_nodes[reference_node]] + model._grid.x_of_node + - model._grid.x_of_node[model._grid.core_nodes[reference_node]] ) ** 2 ) # assert actual and predicted elevations are the same. assert_array_almost_equal( - predicted_z[model.grid.core_nodes], - model.z[model.grid.core_nodes], + predicted_z[model._grid.core_nodes], + model.z[model._grid.core_nodes], decimal=2, ) @@ -105,20 +105,20 @@ def test_diffusion_only(clock_simple, grid_1, U, Model): for _ in range(nts): model.run_one_step(1000) reference_node = 9 - predicted_z = model.z[model.grid.core_nodes[reference_node]] - ( - U / (2. * params["regolith_transport_parameter"]) + predicted_z = model.z[model._grid.core_nodes[reference_node]] - ( + U / (2.0 * params["regolith_transport_parameter"]) ) * ( ( - model.grid.x_of_node - - model.grid.x_of_node[model.grid.core_nodes[reference_node]] + model._grid.x_of_node + - model._grid.x_of_node[model._grid.core_nodes[reference_node]] ) ** 2 ) # assert actual and predicted elevations are the same. assert_array_almost_equal( - predicted_z[model.grid.core_nodes], - model.z[model.grid.core_nodes], + predicted_z[model._grid.core_nodes], + model.z[model._grid.core_nodes], decimal=2, ) @@ -147,19 +147,19 @@ def test_rock_till_linear_diffusion(clock_simple, grid_1, U, Model): for _ in range(nts): model.run_one_step(1000) reference_node = 9 - predicted_z = model.z[model.grid.core_nodes[reference_node]] - ( - U / (2. * params["regolith_transport_parameter"]) + predicted_z = model.z[model._grid.core_nodes[reference_node]] - ( + U / (2.0 * params["regolith_transport_parameter"]) ) * ( ( - model.grid.x_of_node - - model.grid.x_of_node[model.grid.core_nodes[reference_node]] + model._grid.x_of_node + - model._grid.x_of_node[model._grid.core_nodes[reference_node]] ) ** 2 ) # assert actual and predicted elevations are the same. assert_array_almost_equal( - predicted_z[model.grid.core_nodes], - model.z[model.grid.core_nodes], + predicted_z[model._grid.core_nodes], + model.z[model._grid.core_nodes], decimal=2, ) diff --git a/tests/test_model_basicCh.py b/tests/test_model_basicCh.py index 61ac1e96..9174ec4b 100644 --- a/tests/test_model_basicCh.py +++ b/tests/test_model_basicCh.py @@ -43,6 +43,6 @@ def test_diffusion_only(clock_09, grid_4): p_roots = np.roots(p) predicted_slope = np.abs(np.real(p_roots[-1])) actual_slope = np.abs( - model.grid.at_node["topographic__steepest_slope"][39] + model._grid.at_node["topographic__steepest_slope"][39] ) assert_array_almost_equal(actual_slope, predicted_slope, decimal=3) diff --git a/tests/test_model_basicChRt_and_basicChRtTh.py b/tests/test_model_basicChRt_and_basicChRtTh.py index 2ebe0ace..e05dd2ba 100644 --- a/tests/test_model_basicChRt_and_basicChRtTh.py +++ b/tests/test_model_basicChRt_and_basicChRtTh.py @@ -45,6 +45,6 @@ def test_diffusion_only(clock_09, grid_4, Model): # print(predicted_slope) actual_slope = np.abs( - model.grid.at_node["topographic__steepest_slope"][39] + model._grid.at_node["topographic__steepest_slope"][39] ) assert_array_almost_equal(actual_slope, predicted_slope, decimal=3) diff --git a/tests/test_model_basicChSa.py b/tests/test_model_basicChSa.py index 515d9f22..99a91b34 100644 --- a/tests/test_model_basicChSa.py +++ b/tests/test_model_basicChSa.py @@ -39,7 +39,7 @@ def test_diffusion_only(clock_08, grid_4): model.run_one_step(clock_08.step) # test steady state soil depth - actual_depth = model.grid.at_node["soil__depth"][30] + actual_depth = model._grid.at_node["soil__depth"][30] predicted_depth = -soil_production_decay_depth * np.log( U / max_soil_production_rate ) @@ -61,6 +61,6 @@ def test_diffusion_only(clock_08, grid_4): p_roots = np.roots(p) predicted_slope = np.abs(np.real(p_roots[-1])) actual_slope = np.abs( - model.grid.at_node["topographic__steepest_slope"][39] + model._grid.at_node["topographic__steepest_slope"][39] ) assert_array_almost_equal(actual_slope, predicted_slope, decimal=3) diff --git a/tests/test_model_basicDd.py b/tests/test_model_basicDd.py index 36cd76df..fe6bb320 100644 --- a/tests/test_model_basicDd.py +++ b/tests/test_model_basicDd.py @@ -7,13 +7,13 @@ from terrainbento import BasicDd, NotCoreNodeBaselevelHandler -@pytest.mark.parametrize("m_sp", [1. / 3, 0.5]) -@pytest.mark.parametrize("n_sp", [1.]) +@pytest.mark.parametrize("m_sp", [1.0 / 3, 0.5]) +@pytest.mark.parametrize("n_sp", [1.0]) @pytest.mark.parametrize( "depression_finder", [None, "DepressionFinderAndRouter"] ) @pytest.mark.parametrize("threshold", [0.1]) -@pytest.mark.parametrize("thresh_change_per_depth", [0.]) +@pytest.mark.parametrize("thresh_change_per_depth", [0.0]) def test_steady_Ksp_no_precip_changer_no_thresh_change( clock_simple, grid_2, @@ -32,7 +32,7 @@ def test_steady_Ksp_no_precip_changer_no_thresh_change( params = { "grid": grid_2, "clock": clock_simple, - "regolith_transport_parameter": 0., + "regolith_transport_parameter": 0.0, "water_erodibility": K, "m_sp": m_sp, "n_sp": n_sp, @@ -51,23 +51,23 @@ def test_steady_Ksp_no_precip_changer_no_thresh_change( # note that since we have a smooth threshold, we do not have a true # analytical solution, but a bracket within wich we expect the actual # slopes to fall. - actual_slopes = model.grid.at_node["topographic__steepest_slope"] - actual_areas = model.grid.at_node["surface_water__discharge"] + actual_slopes = model._grid.at_node["topographic__steepest_slope"] + actual_areas = model._grid.at_node["surface_water__discharge"] predicted_slopes_upper = ( (U + threshold) / (K * (actual_areas ** m_sp)) - ) ** (1. / n_sp) + ) ** (1.0 / n_sp) predicted_slopes_lower = ((U + 0.0) / (K * (actual_areas ** m_sp))) ** ( - 1. / n_sp + 1.0 / n_sp ) # assert actual and predicted slopes are in the correct range for the # slopes. assert np.all( - actual_slopes[model.grid.core_nodes[1:-1]] - > predicted_slopes_lower[model.grid.core_nodes[1:-1]] + actual_slopes[model._grid.core_nodes[1:-1]] + > predicted_slopes_lower[model._grid.core_nodes[1:-1]] ) assert np.all( - actual_slopes[model.grid.core_nodes[1:-1]] - < predicted_slopes_upper[model.grid.core_nodes[1:-1]] + actual_slopes[model._grid.core_nodes[1:-1]] + < predicted_slopes_upper[model._grid.core_nodes[1:-1]] ) diff --git a/tests/test_model_basicDdHy.py b/tests/test_model_basicDdHy.py index e36bef96..c5ab0767 100644 --- a/tests/test_model_basicDdHy.py +++ b/tests/test_model_basicDdHy.py @@ -8,7 +8,7 @@ @pytest.mark.parametrize("threshold", [0.0, 0.00001]) -@pytest.mark.parametrize("m_sp,n_sp", [(1. / 3, 2. / 3.), (0.5, 1.0)]) +@pytest.mark.parametrize("m_sp,n_sp", [(1.0 / 3, 2.0 / 3.0), (0.5, 1.0)]) @pytest.mark.parametrize( "depression_finder", [None, "DepressionFinderAndRouter"] ) @@ -35,7 +35,7 @@ def test_stream_DdHy( params = { "grid": grid_1, "clock": clock_simple, - "regolith_transport_parameter": 0., + "regolith_transport_parameter": 0.0, "water_erodibility": K, "settling_velocity": v_sc, "sediment_porosity": phi, @@ -55,24 +55,24 @@ def test_stream_DdHy( model.run_one_step(10) # construct actual and predicted slopes - actual_slopes = model.grid.at_node["topographic__steepest_slope"] - actual_areas = model.grid.at_node["surface_water__discharge"] + actual_slopes = model._grid.at_node["topographic__steepest_slope"] + actual_areas = model._grid.at_node["surface_water__discharge"] predicted_slopes = np.power( ((U * v_sc) / (K * np.power(actual_areas, m_sp))) + (U / (K * np.power(actual_areas, m_sp))), - 1. / n_sp, + 1.0 / n_sp, ) # assert actual and predicted slopes are the same. assert_array_almost_equal( - actual_slopes[model.grid.core_nodes[1:-1]], - predicted_slopes[model.grid.core_nodes[1:-1]], + actual_slopes[model._grid.core_nodes[1:-1]], + predicted_slopes[model._grid.core_nodes[1:-1]], decimal=4, ) # assert actual and predicted slopes are the same. assert_array_almost_equal( - actual_slopes[model.grid.core_nodes[1:-1]], - predicted_slopes[model.grid.core_nodes[1:-1]], + actual_slopes[model._grid.core_nodes[1:-1]], + predicted_slopes[model._grid.core_nodes[1:-1]], decimal=4, ) diff --git a/tests/test_model_basicDdRt.py b/tests/test_model_basicDdRt.py index a499a2a9..dcd8d152 100644 --- a/tests/test_model_basicDdRt.py +++ b/tests/test_model_basicDdRt.py @@ -7,13 +7,13 @@ from terrainbento import BasicDdRt, NotCoreNodeBaselevelHandler -@pytest.mark.parametrize("m_sp", [1. / 3, 0.5]) -@pytest.mark.parametrize("n_sp", [1.]) +@pytest.mark.parametrize("m_sp", [1.0 / 3, 0.5]) +@pytest.mark.parametrize("n_sp", [1.0]) @pytest.mark.parametrize( "depression_finder", [None, "DepressionFinderAndRouter"] ) @pytest.mark.parametrize("threshold", [0.1]) -@pytest.mark.parametrize("thresh_change_per_depth", [0.]) +@pytest.mark.parametrize("thresh_change_per_depth", [0.0]) def test_steady_Ksp_no_precip_changer_no_thresh_change( clock_simple, grid_2, @@ -33,7 +33,7 @@ def test_steady_Ksp_no_precip_changer_no_thresh_change( params = { "grid": grid_2, "clock": clock_simple, - "regolith_transport_parameter": 0., + "regolith_transport_parameter": 0.0, "water_erodibility_lower": Kr, "water_erodibility_upper": Kt, "water_erosion_rule__threshold": threshold, @@ -49,10 +49,10 @@ def test_steady_Ksp_no_precip_changer_no_thresh_change( for _ in range(100): model.run_one_step(1000) - actual_slopes = model.grid.at_node["topographic__steepest_slope"] - actual_areas = model.grid.at_node["drainage_area"] - rock_predicted_slopes = (U / (Kr * (actual_areas ** m_sp))) ** (1. / n_sp) - till_predicted_slopes = (U / (Kt * (actual_areas ** m_sp))) ** (1. / n_sp) + actual_slopes = model._grid.at_node["topographic__steepest_slope"] + actual_areas = model._grid.at_node["drainage_area"] + rock_predicted_slopes = (U / (Kr * (actual_areas ** m_sp))) ** (1.0 / n_sp) + till_predicted_slopes = (U / (Kt * (actual_areas ** m_sp))) ** (1.0 / n_sp) # assert actual slopes are steeper than simple stream power prediction assert np.all(actual_slopes[22:37] > rock_predicted_slopes[22:37]) diff --git a/tests/test_model_basicDdSt.py b/tests/test_model_basicDdSt.py index 1ae8fc7d..058e286e 100644 --- a/tests/test_model_basicDdSt.py +++ b/tests/test_model_basicDdSt.py @@ -7,13 +7,13 @@ from terrainbento import BasicDdSt, NotCoreNodeBaselevelHandler -@pytest.mark.parametrize("m_sp", [1. / 3, 0.5]) -@pytest.mark.parametrize("n_sp", [1.]) +@pytest.mark.parametrize("m_sp", [1.0 / 3, 0.5]) +@pytest.mark.parametrize("n_sp", [1.0]) @pytest.mark.parametrize( "depression_finder", [None, "DepressionFinderAndRouter"] ) @pytest.mark.parametrize("threshold", [0.1]) -@pytest.mark.parametrize("thresh_change_per_depth", [0.]) +@pytest.mark.parametrize("thresh_change_per_depth", [0.0]) def test_steady_Ksp_no_precip_changer_no_thresh_change( clock_simple, grid_2, @@ -32,7 +32,7 @@ def test_steady_Ksp_no_precip_changer_no_thresh_change( params = { "grid": grid_2, "clock": clock_simple, - "regolith_transport_parameter": 0., + "regolith_transport_parameter": 0.0, "water_erodibility": K, "m_sp": m_sp, "n_sp": n_sp, @@ -51,23 +51,23 @@ def test_steady_Ksp_no_precip_changer_no_thresh_change( # note that since we have a smooth threshold, we do not have a true # analytical solution, but a bracket within wich we expect the actual # slopes to fall. - actual_slopes = model.grid.at_node["topographic__steepest_slope"] - actual_areas = model.grid.at_node["surface_water__discharge"] + actual_slopes = model._grid.at_node["topographic__steepest_slope"] + actual_areas = model._grid.at_node["surface_water__discharge"] predicted_slopes_upper = ( (U + threshold) / (K * (actual_areas ** m_sp)) - ) ** (1. / n_sp) + ) ** (1.0 / n_sp) predicted_slopes_lower = ((U + 0.0) / (K * (actual_areas ** m_sp))) ** ( - 1. / n_sp + 1.0 / n_sp ) # assert actual and predicted slopes are in the correct range for the # slopes. assert np.all( - actual_slopes[model.grid.core_nodes[1:-1]] - > predicted_slopes_lower[model.grid.core_nodes[1:-1]] + actual_slopes[model._grid.core_nodes[1:-1]] + > predicted_slopes_lower[model._grid.core_nodes[1:-1]] ) assert np.all( - actual_slopes[model.grid.core_nodes[1:-1]] - < predicted_slopes_upper[model.grid.core_nodes[1:-1]] + actual_slopes[model._grid.core_nodes[1:-1]] + < predicted_slopes_upper[model._grid.core_nodes[1:-1]] ) diff --git a/tests/test_model_basicDdVs.py b/tests/test_model_basicDdVs.py index 4cb81262..a43bc9ef 100644 --- a/tests/test_model_basicDdVs.py +++ b/tests/test_model_basicDdVs.py @@ -8,13 +8,13 @@ from terrainbento import BasicDdVs, NotCoreNodeBaselevelHandler -@pytest.mark.parametrize("m_sp", [1. / 3, 0.5]) -@pytest.mark.parametrize("n_sp", [1.]) +@pytest.mark.parametrize("m_sp", [1.0 / 3, 0.5]) +@pytest.mark.parametrize("n_sp", [1.0]) @pytest.mark.parametrize( "depression_finder", [None, "DepressionFinderAndRouter"] ) @pytest.mark.parametrize("threshold", [0.1]) -@pytest.mark.parametrize("thresh_change_per_depth", [0.]) +@pytest.mark.parametrize("thresh_change_per_depth", [0.0]) def test_steady_Ksp_no_precip_changer_no_thresh_change( clock_simple, grid_2, @@ -33,7 +33,7 @@ def test_steady_Ksp_no_precip_changer_no_thresh_change( params = { "grid": grid_2, "clock": clock_simple, - "regolith_transport_parameter": 0., + "regolith_transport_parameter": 0.0, "water_erodibility": K, "water_erosion_rule__threshold": threshold, "water_erosion_rule__thresh_depth_derivative": thresh_change_per_depth, @@ -49,25 +49,25 @@ def test_steady_Ksp_no_precip_changer_no_thresh_change( for _ in range(200): model.run_one_step(1000) - actual_slopes = model.grid.at_node["topographic__steepest_slope"] - actual_areas = model.grid.at_node["surface_water__discharge"] + actual_slopes = model._grid.at_node["topographic__steepest_slope"] + actual_areas = model._grid.at_node["surface_water__discharge"] predicted_slopes_upper = ( (U + threshold) / (K * (actual_areas ** m_sp)) - ) ** (1. / n_sp) + ) ** (1.0 / n_sp) predicted_slopes_lower = ((U + 0.0) / (K * (actual_areas ** m_sp))) ** ( - 1. / n_sp + 1.0 / n_sp ) # assert actual and predicted slopes are in the correct range for the # slopes. assert np.all( - actual_slopes[model.grid.core_nodes[1:-1]] - > predicted_slopes_lower[model.grid.core_nodes[1:-1]] + actual_slopes[model._grid.core_nodes[1:-1]] + > predicted_slopes_lower[model._grid.core_nodes[1:-1]] ) assert np.all( - actual_slopes[model.grid.core_nodes[1:-1]] - < predicted_slopes_upper[model.grid.core_nodes[1:-1]] + actual_slopes[model._grid.core_nodes[1:-1]] + < predicted_slopes_upper[model._grid.core_nodes[1:-1]] ) @@ -81,7 +81,7 @@ def test_Aeff(clock_simple, grid_2, K, U): params = { "grid": grid_2, "clock": clock_simple, - "regolith_transport_parameter": 0., + "regolith_transport_parameter": 0.0, "water_erodibility": K, "water_erosion_rule__threshold": threshold, "hydraulic_conductivity": 0.01, @@ -96,8 +96,8 @@ def test_Aeff(clock_simple, grid_2, K, U): for _ in range(200): model.run_one_step(1000) - actual_slopes = model.grid.at_node["topographic__steepest_slope"] - actual_areas = model.grid.at_node["drainage_area"] + actual_slopes = model._grid.at_node["topographic__steepest_slope"] + actual_areas = model._grid.at_node["drainage_area"] alpha = ( params["hydraulic_conductivity"] @@ -111,7 +111,9 @@ def test_Aeff(clock_simple, grid_2, K, U): # assert aeff internally calculated correclty assert_array_almost_equal( - model.grid.at_node["surface_water__discharge"][model.grid.core_nodes], - A_eff_predicted[model.grid.core_nodes], + model._grid.at_node["surface_water__discharge"][ + model._grid.core_nodes + ], + A_eff_predicted[model._grid.core_nodes], decimal=1, ) diff --git a/tests/test_model_basicHyRt.py b/tests/test_model_basicHyRt.py index fa33aaf3..6e6b6264 100644 --- a/tests/test_model_basicHyRt.py +++ b/tests/test_model_basicHyRt.py @@ -7,7 +7,7 @@ from terrainbento import BasicHyRt, NotCoreNodeBaselevelHandler -@pytest.mark.parametrize("m_sp,n_sp", [(1. / 3, 2. / 3.), (0.5, 1.0)]) +@pytest.mark.parametrize("m_sp,n_sp", [(1.0 / 3, 2.0 / 3.0), (0.5, 1.0)]) @pytest.mark.parametrize( "depression_finder", [None, "DepressionFinderAndRouter"] ) @@ -27,7 +27,7 @@ def test_channel_erosion( params = { "grid": grid_2, "clock": clock_simple, - "regolith_transport_parameter": 0., + "regolith_transport_parameter": 0.0, "water_erodibility_upper": Kt, "water_erodibility_lower": Kr, "settling_velocity": v_sc, @@ -46,17 +46,17 @@ def test_channel_erosion( model.run_one_step(10) # construct actual and predicted slopes - actual_slopes = model.grid.at_node["topographic__steepest_slope"] - actual_areas = model.grid.at_node["drainage_area"] + actual_slopes = model._grid.at_node["topographic__steepest_slope"] + actual_areas = model._grid.at_node["drainage_area"] rock_predicted_slopes = np.power( ((U * v_sc) / (Kr * np.power(actual_areas, m_sp))) + (U / (Kr * np.power(actual_areas, m_sp))), - 1. / n_sp, + 1.0 / n_sp, ) till_predicted_slopes = np.power( ((U * v_sc) / (Kt * np.power(actual_areas, m_sp))) + (U / (Kt * np.power(actual_areas, m_sp))), - 1. / n_sp, + 1.0 / n_sp, ) # assert actual and predicted slopes are the same for rock and till diff --git a/tests/test_model_basicHySa.py b/tests/test_model_basicHySa.py index ef918cc9..cb88ebc0 100644 --- a/tests/test_model_basicHySa.py +++ b/tests/test_model_basicHySa.py @@ -10,7 +10,7 @@ from terrainbento import BasicHySa, NotCoreNodeBaselevelHandler, PrecipChanger -@pytest.mark.parametrize("m_sp,n_sp", [(1. / 3, 2. / 3.), (0.5, 1.0)]) +@pytest.mark.parametrize("m_sp,n_sp", [(1.0 / 3, 2.0 / 3.0), (0.5, 1.0)]) @pytest.mark.parametrize( "depression_finder", [None, "DepressionFinderAndRouter"] ) @@ -36,7 +36,7 @@ def test_channel_erosion( params = { "grid": grid_1, "clock": clock_simple, - "regolith_transport_parameter": 0., + "regolith_transport_parameter": 0.0, "water_erodibility_rock": K_rock_sp, "water_erodibility_sediment": K_sed_sp, "sp_crit_br": sp_crit_br, @@ -61,18 +61,18 @@ def test_channel_erosion( model.run_one_step(10) # construct actual and predicted slopes - actual_slopes = model.grid.at_node["topographic__steepest_slope"] - actual_areas = model.grid.at_node["surface_water__discharge"] + actual_slopes = model._grid.at_node["topographic__steepest_slope"] + actual_areas = model._grid.at_node["surface_water__discharge"] predicted_slopes = np.power( ((U * v_sc) / (K_sed_sp * np.power(actual_areas, m_sp))) + (U / (K_rock_sp * np.power(actual_areas, m_sp))), - 1. / n_sp, + 1.0 / n_sp, ) # assert actual and predicted slopes are the same. assert_array_almost_equal( - actual_slopes[model.grid.core_nodes[1:-1]], - predicted_slopes[model.grid.core_nodes[1:-1]], + actual_slopes[model._grid.core_nodes[1:-1]], + predicted_slopes[model._grid.core_nodes[1:-1]], decimal=4, ) @@ -88,7 +88,7 @@ def test_with_precip_changer( params = { "grid": grid_1, "clock": clock_simple, - "regolith_transport_parameter": 0., + "regolith_transport_parameter": 0.0, "water_erodibility_rock": 0.001, "water_erodibility_sediment": 0.01, "boundary_handlers": {"PrecipChanger": precip_changer}, diff --git a/tests/test_model_basicHyVs.py b/tests/test_model_basicHyVs.py index 22a0bcca..75c168cd 100644 --- a/tests/test_model_basicHyVs.py +++ b/tests/test_model_basicHyVs.py @@ -8,7 +8,7 @@ from terrainbento import BasicHyVs, NotCoreNodeBaselevelHandler -@pytest.mark.parametrize("m_sp,n_sp", [(1. / 3, 2. / 3.), (0.5, 1.0)]) +@pytest.mark.parametrize("m_sp,n_sp", [(1.0 / 3, 2.0 / 3.0), (0.5, 1.0)]) @pytest.mark.parametrize( "depression_finder", [None, "DepressionFinderAndRouter"] ) @@ -26,12 +26,12 @@ def test_channel_erosion( params = { "grid": grid_1, "clock": clock_simple, - "regolith_transport_parameter": 0., + "regolith_transport_parameter": 0.0, "water_erodibility": K, "settling_velocity": v_sc, "sediment_porosity": phi, "fraction_fines": F_f, - "hydraulic_conductivity": 0., + "hydraulic_conductivity": 0.0, "solver": solver, "m_sp": m_sp, "n_sp": n_sp, @@ -45,17 +45,17 @@ def test_channel_erosion( model.run_one_step(5) # construct actual and predicted slopes - actual_slopes = model.grid.at_node["topographic__steepest_slope"] - actual_areas = model.grid.at_node["surface_water__discharge"] + actual_slopes = model._grid.at_node["topographic__steepest_slope"] + actual_areas = model._grid.at_node["surface_water__discharge"] predicted_slopes = np.power( ((U * v_sc) / (K * np.power(actual_areas, m_sp))) + (U / (K * np.power(actual_areas, m_sp))), - 1. / n_sp, + 1.0 / n_sp, ) # assert actual and predicted slopes are the same. assert_array_almost_equal( - actual_slopes[model.grid.core_nodes[1:-1]], - predicted_slopes[model.grid.core_nodes[1:-1]], + actual_slopes[model._grid.core_nodes[1:-1]], + predicted_slopes[model._grid.core_nodes[1:-1]], decimal=4, ) diff --git a/tests/test_model_basicHy_and_basicHySt.py b/tests/test_model_basicHy_and_basicHySt.py index 019017bc..6c40d3b7 100644 --- a/tests/test_model_basicHy_and_basicHySt.py +++ b/tests/test_model_basicHy_and_basicHySt.py @@ -11,7 +11,7 @@ "Model,param_name", [(BasicHy, "water_erodibility"), (BasicHySt, "water_erodibility")], ) -@pytest.mark.parametrize("m_sp,n_sp", [(1. / 3, 2. / 3.), (0.5, 1.0)]) +@pytest.mark.parametrize("m_sp,n_sp", [(1.0 / 3, 2.0 / 3.0), (0.5, 1.0)]) @pytest.mark.parametrize( "depression_finder", [None, "DepressionFinderAndRouter"] ) @@ -38,7 +38,7 @@ def test_channel_erosion( params = { "grid": grid_1, "clock": clock_simple, - "regolith_transport_parameter": 0., + "regolith_transport_parameter": 0.0, param_name: K, "settling_velocity": v_sc, "sediment_porosity": phi, @@ -56,17 +56,17 @@ def test_channel_erosion( model.run_one_step(10) # construct actual and predicted slopes - actual_slopes = model.grid.at_node["topographic__steepest_slope"] - actual_areas = model.grid.at_node["surface_water__discharge"] + actual_slopes = model._grid.at_node["topographic__steepest_slope"] + actual_areas = model._grid.at_node["surface_water__discharge"] predicted_slopes = np.power( ((U * v_sc) / (K * np.power(actual_areas, m_sp))) + (U / (K * np.power(actual_areas, m_sp))), - 1. / n_sp, + 1.0 / n_sp, ) # assert actual and predicted slopes are the same. assert_array_almost_equal( - actual_slopes[model.grid.core_nodes[1:-1]], - predicted_slopes[model.grid.core_nodes[1:-1]], + actual_slopes[model._grid.core_nodes[1:-1]], + predicted_slopes[model._grid.core_nodes[1:-1]], decimal=4, ) diff --git a/tests/test_model_basicRtSa_basicSa_basicSaVs.py b/tests/test_model_basicRtSa_basicSa_basicSaVs.py index 4a8b1ac7..edffe338 100644 --- a/tests/test_model_basicRtSa_basicSa_basicSaVs.py +++ b/tests/test_model_basicRtSa_basicSa_basicSaVs.py @@ -29,7 +29,7 @@ def test_diffusion_only(clock_simple, grid_4, Model, water_params): regolith_transport_parameter = 1.0 soil_transport_decay_depth = 0.5 - grid_4.at_node["soil__depth"][:] = 0. + grid_4.at_node["soil__depth"][:] = 0.0 ncnblh = NotCoreNodeBaselevelHandler( grid_4, modify_core_nodes=True, lowering_rate=-U @@ -55,18 +55,18 @@ def test_diffusion_only(clock_simple, grid_4, Model, water_params): dx = grid_4.dx # test steady state soil depthf - actual_depth = model.grid.at_node["soil__depth"][28] + actual_depth = model._grid.at_node["soil__depth"][28] predicted_depth = -soil_production_decay_depth * np.log( U / max_soil_production_rate ) assert_array_almost_equal(actual_depth, predicted_depth, decimal=3) # test steady state slope - actual_profile = model.grid.at_node["topographic__elevation"][21:42] + actual_profile = model._grid.at_node["topographic__elevation"][21:42] - domain = np.arange(0, max(model.grid.node_x + dx), dx) + domain = np.arange(0, max(model._grid.node_x + dx), dx) - half_domain = np.arange(0, max(domain) / 2. + dx, dx) + half_domain = np.arange(0, max(domain) / 2.0 + dx, dx) one_minus_h_hstar = 1 - np.exp( -predicted_depth / soil_transport_decay_depth @@ -75,7 +75,7 @@ def test_diffusion_only(clock_simple, grid_4, Model, water_params): half_domain_z = ( -half_domain ** 2 * U - / (regolith_transport_parameter * 2. * one_minus_h_hstar) + / (regolith_transport_parameter * 2.0 * one_minus_h_hstar) ) steady_z_profile = np.concatenate( diff --git a/tests/test_model_basicRtTh.py b/tests/test_model_basicRtTh.py index 4548e689..2f705fa9 100644 --- a/tests/test_model_basicRtTh.py +++ b/tests/test_model_basicRtTh.py @@ -27,7 +27,7 @@ def test_steady_Ksp_no_precip_changer( params = { "grid": grid_2, "clock": clock_simple, - "regolith_transport_parameter": 0., + "regolith_transport_parameter": 0.0, "water_erodibility_lower": Kr, "water_erodibility_upper": Kt, "water_erosion_rule_upper__threshold": Tt, @@ -43,25 +43,25 @@ def test_steady_Ksp_no_precip_changer( for _ in range(200): model.run_one_step(1000) - actual_slopes = model.grid.at_node["topographic__steepest_slope"] - actual_areas = model.grid.at_node["surface_water__discharge"] + actual_slopes = model._grid.at_node["topographic__steepest_slope"] + actual_areas = model._grid.at_node["surface_water__discharge"] # note that since we have a smooth threshold, we do not have a true # analytical solution, but a bracket within wich we expect the actual # slopes to fall. rock_predicted_slopes_upper = ( (U + Tr) / (Kr * (actual_areas ** m_sp)) - ) ** (1. / n_sp) + ) ** (1.0 / n_sp) till_predicted_slopes_upper = ( (U + Tt) / (Kt * (actual_areas ** m_sp)) - ) ** (1. / n_sp) + ) ** (1.0 / n_sp) rock_predicted_slopes_lower = ( - (U + 0.) / (Kr * (actual_areas ** m_sp)) - ) ** (1. / n_sp) + (U + 0.0) / (Kr * (actual_areas ** m_sp)) + ) ** (1.0 / n_sp) till_predicted_slopes_lower = ( - (U + 0.) / (Kt * (actual_areas ** m_sp)) - ) ** (1. / n_sp) + (U + 0.0) / (Kt * (actual_areas ** m_sp)) + ) ** (1.0 / n_sp) # assert actual and predicted slopes are the same for rock and till # portions. diff --git a/tests/test_model_basicRtVs.py b/tests/test_model_basicRtVs.py index bef9ddec..8def4ca4 100644 --- a/tests/test_model_basicRtVs.py +++ b/tests/test_model_basicRtVs.py @@ -4,7 +4,7 @@ from terrainbento import BasicRtVs, NotCoreNodeBaselevelHandler -@pytest.mark.parametrize("m_sp,n_sp", [(1. / 3, 2. / 3.), (0.5, 1.0)]) +@pytest.mark.parametrize("m_sp,n_sp", [(1.0 / 3, 2.0 / 3.0), (0.5, 1.0)]) @pytest.mark.parametrize( "depression_finder", [None, "DepressionFinderAndRouter"] ) @@ -21,7 +21,7 @@ def test_steady_Kss_no_precip_changer( params = { "grid": grid_2, "clock": clock_simple, - "regolith_transport_parameter": 0., + "regolith_transport_parameter": 0.0, "water_erodibility_lower": Kr, "water_erodibility_upper": Kt, "hydraulic_conductivity": hydraulic_conductivity, @@ -36,10 +36,10 @@ def test_steady_Kss_no_precip_changer( for _ in range(100): model.run_one_step(1000) - actual_slopes = model.grid.at_node["topographic__steepest_slope"] - actual_areas = model.grid.at_node["surface_water__discharge"] - rock_predicted_slopes = (U / (Kr * (actual_areas ** m_sp))) ** (1. / n_sp) - till_predicted_slopes = (U / (Kt * (actual_areas ** m_sp))) ** (1. / n_sp) + actual_slopes = model._grid.at_node["topographic__steepest_slope"] + actual_areas = model._grid.at_node["surface_water__discharge"] + rock_predicted_slopes = (U / (Kr * (actual_areas ** m_sp))) ** (1.0 / n_sp) + till_predicted_slopes = (U / (Kt * (actual_areas ** m_sp))) ** (1.0 / n_sp) # assert actual and predicted slopes are the same for rock and till. assert_array_almost_equal( diff --git a/tests/test_model_basicSaVs.py b/tests/test_model_basicSaVs.py index 3c4d8a90..e7f1c746 100644 --- a/tests/test_model_basicSaVs.py +++ b/tests/test_model_basicSaVs.py @@ -7,7 +7,7 @@ from terrainbento import BasicSaVs, NotCoreNodeBaselevelHandler -@pytest.mark.parametrize("m_sp,n_sp", [(1. / 3, 2. / 3.), (0.5, 1.0)]) +@pytest.mark.parametrize("m_sp,n_sp", [(1.0 / 3, 2.0 / 3.0), (0.5, 1.0)]) @pytest.mark.parametrize( "depression_finder", [None, "DepressionFinderAndRouter"] ) @@ -22,7 +22,7 @@ def test_detachment_steady_no_precip_changer( params = { "grid": grid_1, "clock": clock_simple, - "regolith_transport_parameter": 0., + "regolith_transport_parameter": 0.0, "water_erodibility": 0.001, "depression_finder": depression_finder, "m_sp": m_sp, @@ -36,14 +36,14 @@ def test_detachment_steady_no_precip_changer( model.run_one_step(1000) # construct actual and predicted slopes - actual_slopes = model.grid.at_node["topographic__steepest_slope"] - actual_areas = model.grid.at_node["surface_water__discharge"] + actual_slopes = model._grid.at_node["topographic__steepest_slope"] + actual_areas = model._grid.at_node["surface_water__discharge"] predicted_slopes = ( U / (params["water_erodibility"] * (actual_areas ** params["m_sp"])) - ) ** (1. / params["n_sp"]) + ) ** (1.0 / params["n_sp"]) # assert actual and predicted slopes are the same. assert_array_almost_equal( - actual_slopes[model.grid.core_nodes[1:-1]], - predicted_slopes[model.grid.core_nodes[1:-1]], + actual_slopes[model._grid.core_nodes[1:-1]], + predicted_slopes[model._grid.core_nodes[1:-1]], ) diff --git a/tests/test_model_basicStVs.py b/tests/test_model_basicStVs.py index bc0a2b0e..321242b9 100644 --- a/tests/test_model_basicStVs.py +++ b/tests/test_model_basicStVs.py @@ -8,7 +8,7 @@ def test_bad_transmiss(grid_2, clock_simple): params = { "grid": grid_2, "clock": clock_simple, - "hydraulic_conductivity": 0., + "hydraulic_conductivity": 0.0, } with pytest.raises(ValueError): @@ -30,7 +30,7 @@ def test_steady_without_stochastic_duration( params = { "grid": grid_2, "clock": clock_simple, - "regolith_transport_parameter": 0., + "regolith_transport_parameter": 0.0, "water_erodibility": K, "m_sp": m_sp, "n_sp": n_sp, @@ -50,10 +50,10 @@ def test_steady_without_stochastic_duration( model.run_one_step(1.0) # construct actual and predicted slopes - ic = model.grid.core_nodes[1:-1] # "inner" core nodes - actual_slopes = model.grid.at_node["topographic__steepest_slope"] - actual_areas = model.grid.at_node["surface_water__discharge"] - predicted_slopes = (U / (K * (actual_areas ** m_sp))) ** (1. / n_sp) + ic = model._grid.core_nodes[1:-1] # "inner" core nodes + actual_slopes = model._grid.at_node["topographic__steepest_slope"] + actual_areas = model._grid.at_node["surface_water__discharge"] + predicted_slopes = (U / (K * (actual_areas ** m_sp))) ** (1.0 / n_sp) assert_array_almost_equal( actual_slopes[ic], predicted_slopes[ic], decimal=4 ) diff --git a/tests/test_model_basicSt_basicStTh.py b/tests/test_model_basicSt_basicStTh.py index 363ba8ac..c76ae9da 100644 --- a/tests/test_model_basicSt_basicStTh.py +++ b/tests/test_model_basicSt_basicStTh.py @@ -83,7 +83,7 @@ def test_steady_without_stochastic_duration( m = 1.0 n = 1.0 - grid = RasterModelGrid((3, 6), xy_spacing=100.) + grid = RasterModelGrid((3, 6), xy_spacing=100.0) grid.set_closed_boundaries_at_grid_edges(False, True, False, True) grid.add_zeros("node", "topographic__elevation") s = grid.add_zeros("node", "soil__depth") @@ -97,7 +97,7 @@ def test_steady_without_stochastic_duration( params = { "grid": grid, "clock": clock_simple, - "regolith_transport_parameter": 0., + "regolith_transport_parameter": 0.0, "water_erodibility": K, "m_sp": m, "n_sp": n, @@ -119,9 +119,9 @@ def test_steady_without_stochastic_duration( model.run_one_step(1.0) # construct actual and predicted slopes - ic = model.grid.core_nodes[1:-1] # "inner" core nodes - actual_slopes = model.grid.at_node["topographic__steepest_slope"][ic] - actual_areas = model.grid.at_node["drainage_area"][ic] + ic = model._grid.core_nodes[1:-1] # "inner" core nodes + actual_slopes = model._grid.at_node["topographic__steepest_slope"][ic] + actual_areas = model._grid.at_node["drainage_area"][ic] predicted_slopes = 2 * U / (K * (actual_areas)) # assert actual and predicted slopes are the same. @@ -138,7 +138,7 @@ def test_stochastic_duration_rainfall_means(): m = 1.0 n = 1.0 - grid = RasterModelGrid((3, 6), xy_spacing=100.) + grid = RasterModelGrid((3, 6), xy_spacing=100.0) grid.set_closed_boundaries_at_grid_edges(True, False, True, False) grid.add_zeros("node", "topographic__elevation") @@ -151,7 +151,7 @@ def test_stochastic_duration_rainfall_means(): params = { "grid": grid, "clock": clock, - "regolith_transport_parameter": 0., + "regolith_transport_parameter": 0.0, "water_erodibility": K, "m_sp": m, "n_sp": n, diff --git a/tests/test_model_basicTh.py b/tests/test_model_basicTh.py index 08fa4467..f1f88db1 100644 --- a/tests/test_model_basicTh.py +++ b/tests/test_model_basicTh.py @@ -21,7 +21,7 @@ def test_steady_Kss_no_precip_changer( params = { "grid": grid_2, "clock": clock_simple, - "regolith_transport_parameter": 0., + "regolith_transport_parameter": 0.0, "m_sp": m_sp, "n_sp": n_sp, "water_erodibility": K, @@ -35,23 +35,23 @@ def test_steady_Kss_no_precip_changer( for _ in range(200): model.run_one_step(1000) - actual_slopes = model.grid.at_node["topographic__steepest_slope"] - actual_areas = model.grid.at_node["surface_water__discharge"] + actual_slopes = model._grid.at_node["topographic__steepest_slope"] + actual_areas = model._grid.at_node["surface_water__discharge"] predicted_slopes_upper = ( (U + threshold) / (K * (actual_areas ** m_sp)) - ) ** (1. / n_sp) + ) ** (1.0 / n_sp) predicted_slopes_lower = ((U + 0.0) / (K * (actual_areas ** m_sp))) ** ( - 1. / n_sp + 1.0 / n_sp ) # assert actual and predicted slopes are in the correct range for the # slopes. assert np.all( - actual_slopes[model.grid.core_nodes[1:-1]] - > predicted_slopes_lower[model.grid.core_nodes[1:-1]] + actual_slopes[model._grid.core_nodes[1:-1]] + > predicted_slopes_lower[model._grid.core_nodes[1:-1]] ) assert np.all( - actual_slopes[model.grid.core_nodes[1:-1]] - < predicted_slopes_upper[model.grid.core_nodes[1:-1]] + actual_slopes[model._grid.core_nodes[1:-1]] + < predicted_slopes_upper[model._grid.core_nodes[1:-1]] ) diff --git a/tests/test_model_basicThVs.py b/tests/test_model_basicThVs.py index c7d639c6..06a46ff4 100644 --- a/tests/test_model_basicThVs.py +++ b/tests/test_model_basicThVs.py @@ -25,7 +25,7 @@ def test_steady_Kss_no_precip_changer( params = { "grid": grid_2, "clock": clock_simple, - "regolith_transport_parameter": 0., + "regolith_transport_parameter": 0.0, "water_erodibility": K, "hydraulic_conductivity": hydraulic_conductivity, "m_sp": m_sp, @@ -40,23 +40,23 @@ def test_steady_Kss_no_precip_changer( for _ in range(200): model.run_one_step(1000) - actual_slopes = model.grid.at_node["topographic__steepest_slope"] - actual_areas = model.grid.at_node["surface_water__discharge"] + actual_slopes = model._grid.at_node["topographic__steepest_slope"] + actual_areas = model._grid.at_node["surface_water__discharge"] predicted_slopes_upper = ( (U + threshold) / (K * (actual_areas ** m_sp)) - ) ** (1. / n_sp) + ) ** (1.0 / n_sp) predicted_slopes_lower = ((U + 0.0) / (K * (actual_areas ** m_sp))) ** ( - 1. / n_sp + 1.0 / n_sp ) # assert actual and predicted slopes are in the correct range for the # slopes. assert np.all( - actual_slopes[model.grid.core_nodes[1:-1]] - > predicted_slopes_lower[model.grid.core_nodes[1:-1]] + actual_slopes[model._grid.core_nodes[1:-1]] + > predicted_slopes_lower[model._grid.core_nodes[1:-1]] ) assert np.all( - actual_slopes[model.grid.core_nodes[1:-1]] - < predicted_slopes_upper[model.grid.core_nodes[1:-1]] + actual_slopes[model._grid.core_nodes[1:-1]] + < predicted_slopes_upper[model._grid.core_nodes[1:-1]] ) diff --git a/tests/test_model_basicVs.py b/tests/test_model_basicVs.py index 20317044..c2d0b909 100644 --- a/tests/test_model_basicVs.py +++ b/tests/test_model_basicVs.py @@ -6,7 +6,7 @@ from terrainbento import BasicVs, NotCoreNodeBaselevelHandler -@pytest.mark.parametrize("m_sp,n_sp", [(1. / 3, 2. / 3.), (0.5, 1.0)]) +@pytest.mark.parametrize("m_sp,n_sp", [(1.0 / 3, 2.0 / 3.0), (0.5, 1.0)]) @pytest.mark.parametrize( "depression_finder", [None, "DepressionFinderAndRouter"] ) @@ -19,7 +19,7 @@ def test_detachment_steady_no_precip_changer( params = { "grid": grid_1, "clock": clock_simple, - "regolith_transport_parameter": 0., + "regolith_transport_parameter": 0.0, "water_erodibility": 0.001, "hydraulic_conductivity": 1.0, "depression_finder": depression_finder, @@ -33,14 +33,14 @@ def test_detachment_steady_no_precip_changer( model.run_one_step(1000) # construct actual and predicted slopes - actual_slopes = model.grid.at_node["topographic__steepest_slope"] - actual_areas = model.grid.at_node["surface_water__discharge"] + actual_slopes = model._grid.at_node["topographic__steepest_slope"] + actual_areas = model._grid.at_node["surface_water__discharge"] predicted_slopes = ( U / (params["water_erodibility"] * (actual_areas ** params["m_sp"])) - ) ** (1. / params["n_sp"]) + ) ** (1.0 / params["n_sp"]) # assert actual and predicted slopes are the same. assert_array_almost_equal( - actual_slopes[model.grid.core_nodes[1:-1]], - predicted_slopes[model.grid.core_nodes[1:-1]], + actual_slopes[model._grid.core_nodes[1:-1]], + predicted_slopes[model._grid.core_nodes[1:-1]], ) diff --git a/tests/test_precip_changer.py b/tests/test_precip_changer.py index d31de732..e61cd2a8 100644 --- a/tests/test_precip_changer.py +++ b/tests/test_precip_changer.py @@ -53,7 +53,7 @@ def test_simple_precip_changer( params = { "grid": grid_1, "clock": clock_simple, - "regolith_transport_parameter": 0., + "regolith_transport_parameter": 0.0, "water_erodibility": K, "boundary_handlers": {"PrecipChanger": precip_changer}, } @@ -98,10 +98,10 @@ def test_rock_till_precip_changer( model._update_erodibility_field() assert ( - np.array_equiv(model.eroder.K[model.grid.core_nodes[:8]], Kt) is True + np.array_equiv(model.eroder.K[model._grid.core_nodes[:8]], Kt) is True ) assert ( - np.array_equiv(model.eroder.K[model.grid.core_nodes[10:]], Kr) is True + np.array_equiv(model.eroder.K[model._grid.core_nodes[10:]], Kr) is True ) assert "PrecipChanger" in model.boundary_handlers @@ -109,10 +109,10 @@ def test_rock_till_precip_changer( model.run_one_step(1.0) assert_array_almost_equal( - model.eroder.K[model.grid.core_nodes[:8]], + model.eroder.K[model._grid.core_nodes[:8]], Kt * precip_testing_factor * np.ones((8)), ) assert_array_almost_equal( - model.eroder.K[model.grid.core_nodes[10:]], + model.eroder.K[model._grid.core_nodes[10:]], Kr * precip_testing_factor * np.ones((9)), ) diff --git a/tests/test_simple_detachment_limited.py b/tests/test_simple_detachment_limited.py index 639bcf57..6347f219 100644 --- a/tests/test_simple_detachment_limited.py +++ b/tests/test_simple_detachment_limited.py @@ -17,8 +17,8 @@ @pytest.mark.parametrize("Model", [BasicRt, BasicChRt, BasicRtSa]) -@pytest.mark.parametrize("m_sp", [1. / 3, 0.5, 0.75, 0.25]) -@pytest.mark.parametrize("n_sp", [2. / 3., 1.]) +@pytest.mark.parametrize("m_sp", [1.0 / 3, 0.5, 0.75, 0.25]) +@pytest.mark.parametrize("n_sp", [2.0 / 3.0, 1.0]) @pytest.mark.parametrize( "depression_finder", [None, "DepressionFinderAndRouter"] ) @@ -31,7 +31,7 @@ def test_rock_till_steady_no_precip_changer( params = { "grid": grid_2, "clock": clock_simple, - "regolith_transport_parameter": 0., + "regolith_transport_parameter": 0.0, "water_erodibility_lower": Kr, "water_erodibility_upper": Kt, "depression_finder": depression_finder, @@ -46,10 +46,10 @@ def test_rock_till_steady_no_precip_changer( model.run_one_step(1000) # construct actual and predicted slopes - actual_slopes = model.grid.at_node["topographic__steepest_slope"] - actual_areas = model.grid.at_node["surface_water__discharge"] - rock_predicted_slopes = (U / (Kr * (actual_areas ** m_sp))) ** (1. / n_sp) - till_predicted_slopes = (U / (Kt * (actual_areas ** m_sp))) ** (1. / n_sp) + actual_slopes = model._grid.at_node["topographic__steepest_slope"] + actual_areas = model._grid.at_node["surface_water__discharge"] + rock_predicted_slopes = (U / (Kr * (actual_areas ** m_sp))) ** (1.0 / n_sp) + till_predicted_slopes = (U / (Kt * (actual_areas ** m_sp))) ** (1.0 / n_sp) # assert actual and predicted slopes are the same for rock and till # portions. @@ -63,8 +63,8 @@ def test_rock_till_steady_no_precip_changer( @pytest.mark.parametrize("Model", [BasicRtTh, BasicChRtTh]) -@pytest.mark.parametrize("m_sp", [1. / 3, 0.5, 0.75, 0.25]) -@pytest.mark.parametrize("n_sp", [1.]) +@pytest.mark.parametrize("m_sp", [1.0 / 3, 0.5, 0.75, 0.25]) +@pytest.mark.parametrize("n_sp", [1.0]) @pytest.mark.parametrize( "depression_finder", [None, "DepressionFinderAndRouter"] ) @@ -77,7 +77,7 @@ def test_rock_till_steady_no_precip_changer_ChRtTh( params = { "grid": grid_2, "clock": clock_simple, - "regolith_transport_parameter": 0., + "regolith_transport_parameter": 0.0, "water_erodibility_lower": Kr, "water_erodibility_upper": Kt, "depression_finder": depression_finder, @@ -94,10 +94,10 @@ def test_rock_till_steady_no_precip_changer_ChRtTh( model.run_one_step(1000) # construct actual and predicted slopes - actual_slopes = model.grid.at_node["topographic__steepest_slope"] - actual_areas = model.grid.at_node["surface_water__discharge"] - rock_predicted_slopes = (U / (Kr * (actual_areas ** m_sp))) ** (1. / n_sp) - till_predicted_slopes = (U / (Kt * (actual_areas ** m_sp))) ** (1. / n_sp) + actual_slopes = model._grid.at_node["topographic__steepest_slope"] + actual_areas = model._grid.at_node["surface_water__discharge"] + rock_predicted_slopes = (U / (Kr * (actual_areas ** m_sp))) ** (1.0 / n_sp) + till_predicted_slopes = (U / (Kt * (actual_areas ** m_sp))) ** (1.0 / n_sp) # assert actual and predicted slopes are the same for rock and till # portions. @@ -113,8 +113,8 @@ def test_rock_till_steady_no_precip_changer_ChRtTh( @pytest.mark.parametrize( "Model", [Basic, BasicCv, BasicCh, BasicChSa, BasicSa] ) -@pytest.mark.parametrize("m_sp", [1. / 3, 0.5, 0.75, 0.25]) -@pytest.mark.parametrize("n_sp", [2. / 3., 1.]) +@pytest.mark.parametrize("m_sp", [1.0 / 3, 0.5, 0.75, 0.25]) +@pytest.mark.parametrize("n_sp", [2.0 / 3.0, 1.0]) @pytest.mark.parametrize( "depression_finder", [None, "DepressionFinderAndRouter"] ) @@ -127,7 +127,7 @@ def test_detachment_steady_no_precip_changer( params = { "grid": grid_1, "clock": clock_simple, - "regolith_transport_parameter": 0., + "regolith_transport_parameter": 0.0, "water_erodibility": 0.001, "depression_finder": depression_finder, "m_sp": m_sp, @@ -140,14 +140,14 @@ def test_detachment_steady_no_precip_changer( model.run_one_step(1000) # construct actual and predicted slopes - actual_slopes = model.grid.at_node["topographic__steepest_slope"] - actual_areas = model.grid.at_node["surface_water__discharge"] + actual_slopes = model._grid.at_node["topographic__steepest_slope"] + actual_areas = model._grid.at_node["surface_water__discharge"] predicted_slopes = ( U / (params["water_erodibility"] * (actual_areas ** params["m_sp"])) - ) ** (1. / params["n_sp"]) + ) ** (1.0 / params["n_sp"]) # assert actual and predicted slopes are the same. assert_array_almost_equal( - actual_slopes[model.grid.core_nodes[1:-1]], - predicted_slopes[model.grid.core_nodes[1:-1]], + actual_slopes[model._grid.core_nodes[1:-1]], + predicted_slopes[model._grid.core_nodes[1:-1]], ) diff --git a/tests/test_vsa_aeff.py b/tests/test_vsa_aeff.py index f9fce8f1..48b325bf 100644 --- a/tests/test_vsa_aeff.py +++ b/tests/test_vsa_aeff.py @@ -23,7 +23,7 @@ def test_Aeff(clock_simple, grid_2, K, U, Model): params = { "grid": grid_2, "clock": clock_simple, - "regolith_transport_parameter": 0., + "regolith_transport_parameter": 0.0, "water_erodibility": K, "hydraulic_conductivity": 0.02, "m_sp": m_sp, @@ -38,8 +38,8 @@ def test_Aeff(clock_simple, grid_2, K, U, Model): model.run_one_step(10) # construct actual and predicted slopes - actual_slopes = model.grid.at_node["topographic__steepest_slope"] - actual_areas = model.grid.at_node["drainage_area"] + actual_slopes = model._grid.at_node["topographic__steepest_slope"] + actual_areas = model._grid.at_node["drainage_area"] alpha = ( params["hydraulic_conductivity"] @@ -54,7 +54,9 @@ def test_Aeff(clock_simple, grid_2, K, U, Model): # assert aeff internally calculated correclty assert_array_almost_equal( - model.grid.at_node["surface_water__discharge"][model.grid.core_nodes], - A_eff_predicted[model.grid.core_nodes], + model._grid.at_node["surface_water__discharge"][ + model._grid.core_nodes + ], + A_eff_predicted[model._grid.core_nodes], decimal=1, ) diff --git a/tests/test_wrap_as_bmi.py b/tests/test_wrap_as_bmi.py new file mode 100644 index 00000000..20d52469 --- /dev/null +++ b/tests/test_wrap_as_bmi.py @@ -0,0 +1,22 @@ +import sys +import warnings + +import cfunits +import pytest + +from terrainbento.derived_models import MODELS + + +@pytest.mark.parametrize("Model", MODELS) +def test_wrap_as_bmi(Model): + + # verify that all models BMI attributes are the correct values in the + # docstring information. + + # verify that all the grid information is correct. + + # verify that the model can run as expected. + + # verify that all units are either None or UDUNITS compatible. + + pass