ecohydrologyteam / clearwater-riverine Goto Github PK

The transport equation includes $n + \theta$, a temporal weighting factor on the concentrations in the advection and diffusion terms. The current Clearwater-riverine model assumes that $\theta = 1$ (i.e., a first-order fully implicit Backward Euler scheme), but users should be able to specify a $\theta$ of their choice (e.g. $\theta = 0.5$ for a second-order Crank-Nicholson scheme). In general, $X^{n+\theta} = \theta X^{n+1} + (1-\theta)X^n$ (Generalized Euler scheme).

This issue can be closed when:

1. Users can specify $\theta$
2. A comparison has been run to determine how changing $theta$ impacts accuracy for models of various timesteps (e.g., how does a 5 second model's mass balance change when $\theta =1$ versus $\theta = 0.5$?).

The model that we worked with for all of our testing had both face flow and cell volume printed to output (by default, this information is NOT printed, but we turned it on). If those values are not included, the face flow and cell volume should be calculated by the Clearwater Riverine model.

However, in some recent testing with a different model, @jrutyna found that there were errors thrown when the the face flow and cell volume outputs were not turned on. Therefore, this error handling and/or calculation is not properly implemented and should be addressed.

The boundary_conditions() method is the slowest method. Explore possibilities for speeding this calculation up. I believe the slowest piece of this method is the following line:

boundary_df['Time Index'] = [np.where(self.mesh.time == i)[0][0] for i in boundary_df.Datetime]

Some xarray methods could likely be leveraged to improve this code.

@jrutyna encountered the following mass balance issue:

I conducted the 100 everywhere initial condition test and boundary condition test on the Ohio River model and it also is not balancing mass...it loses about 0.5% by the time it reaches the downstream boundary. I think it is related to the boundary conditions.

I believe this issue is associated with ghost cell volumes. In our initial testing models, we did not have adequate information to calculate ghost cell volumes, so we approximated those volumes using the values flowing in/out of the cell. This would over-approximate the ghost cell volumes.

An exploration of RAS output of a different model now shows that the ghost cell volumes can be printed. We should remove the ghost cell approximation calculations (or only make them triggered if required).

As described in the The Hitchhiker’s Guide to Python section on Structuring Your Project:

Repository structure is a crucial part of your project’s architecture.

Following their guidance, and, looking at other similar repositories, I propose that we have the following directory structure:

• clearwater_riverine
• docs
• examples
  • data
  • dev_sandbox
  • temp
• tests

Here's the intended purpose of each of these directories:

• clearwater_riverine is package of source code that users will import.
• docs is for all documentation, including for sphinx.
• examples is for Jupyter notebooks to teach users how to use the package. If designed as a tutorial, notebooks might be numbered in the sequence they should be used (i.e. 1_Intro.ipynb, etc.).
• examples/data is for data files required to run the example notebooks.
• examples/dev_sandbox is for exploratory work by the development team as they develop new capabilities, including for informal testing.
• examples/temp is for temporary output files that users might create during the demo. This subdirectory will be in the .gitignore, so that output files are not saved to the repo.
• tests is for unit testing scripts (i.e. test_basic.py) designed for automated testing. Example models do not belong in this folder, nor do Jupyter Notebooks.Testing scripts in this folder may access use case datasets in examples/data.

@TSteissberg, @sjordan29, @ptomasula, what do you all think?

Continue building out unit tests with pytest in this file for unit testing of Clearwater Riverine.

• Use codecov to ensure complete coverage.
• Set up workflow to run unit tests on commit
• Set up workflow to ensure all unit tests pass before merging to main

We are going to compare E. Coli concentrations from an EFDC model to the output from our RAS2D water quality model. This will involve the following steps:

1. Use EFDC boundary conditions to define water quality boundary conditions [internal and external?] in the RAS2D water quality model.
2. Figure out how to compare downstream boundary conditions with results from RAS2D. The two models have different spatial resolutions.

In place of certain print statements, switch to logging information so we can better control when and how information is displayed. Will need to research if logging is compatible with numba.

@sjordan29, @ptomasula, and I decided today that it made sense to separate Input/Output functions from the model class and into it's own class.

To get started, @ptomasula suggested that we might use the Factory Method creational design pattern to create class methods (using the @classmethod decorator (see https://pynative.com/python-class-method/) to create grid objects and initial/boundary condition objects from the RAS2D model HDF5 output and other files.

Another approach might be to implement implicit interfaces using Protocol Classes, which enables static type checking (via MyPy) of functions. I'm not sure if this is mutually exclusive with Factory Method approach or not.

Creating this issue to track our exploration of different approaches to solve this.

@jrutyna & @aufdenkampe, just merged #54 to address #47, so I think we are ready to assess whether it'd be possible to have a longer time step than 1 second to decrease model run time.

@jrutyna - can you set up a few models with various timestep lengths (maybe ranging from 2-30 second timesteps depending on what you think is reasonable?). How long would the total model run need to be for you to feel comfortable assessing the impact? For example, could we do a 1-hour or 3-hour model run rather than a day or several day model run for this classification?

Let me know if you want to do the analysis between model runs or if you want me to take care of that. I can use your post-processing scripts if you want me to compare and I can present the results here, or you can take over. Whatever works best with your schedule this week!

Use clearwater-modules-python to demo TSM. This issue can be closed when there is an initial implementation of TSM for the October deadline. Future issues may focus on improving the integration between Clearwater-riverine and Clearwater-modules-python.

I used memory_profiler for instantiation of Clearwater-riverine and running a few timesteps. The following were interesting findings:

• the calculate_required_parameters method is very memory intensive, because transport calculations require a few key variables with dimensions of time by the number of edges, which get very large as the time dimension grows.
• turning off printing intermediate variables in clearwater-modules helped with memory management

Next steps:

• Delete unnecessary intermediate variables from Clearwater-riverine
• Investigate if we can remove FLOW_ACROSS_FACE variable via #60 to reduce memory-intensive redundancy
• Leverage lazy loading. See this pangeo discussion and the reading/writing xarray docs.

Develop a test case in the tests folder in which initial conditions and boundary conditions are all set to 100 (arbitrary units).

In the Box Model example, we find that the .prj file says there are SI units.

However, when we parse the attributes from the HDF output file, we find Imperial units (feet, feet per second, etc.).

@CTaylorLTI - do you have any idea of which units are correct? Are there different units when you are setting up the model vs. in model output?

@jrutyna noticed that there was a mass balance issue in plan 241. Upon investigating it was determined that the mass flux variables for the downstream boundary was only using the boundary condition of 25 mg/L for this calculation and not using the concentration values in the cells next to the boundary. @sjordan29 should investigate further to determine how to fix this issue.

@drucinski suggested via email to @tslawecki, @TSteissberg, and me that we might be able to this repo and associated open-access paper. Perhaps there is something in here that we might use, or maybe this is just a motivation to publish this work in 2023.

Shabani, A., Woznicki, S.A., Mehaffey, M., Butcher, J., Wool, T.A. and Whung, P.Y., 2021. A coupled hydrodynamic (HEC‐RAS 2D) and water quality model (WASP) for simulating flood‐induced soil, sediment, and contaminant transport. Journal of Flood Risk Management, 14(4), p.e12747.

• https://onlinelibrary.wiley.com/doi/10.1111/jfr3.12747

• Git repo with python script and control file format:
  https://github.com/AfshinShabani/External-Coupler-

My quick takes are:

• This is a single 600-line script that runs sequentially, not software.
• There are only 5 commits, from a single person, made around the time that the publication went to press
• There are some blocks of code that could be useful, if/when we want to:
  • Read/write WASP files
  • Convert a finer scale RAS grid into a courser scale water quality grid

cc: @sjordan29

While trying to clear up some memory for Clearwater-riverine, I found an interesting discrepancy between EDGE_VELOCITY and FLOW_ACROSS_FACE. There are instances when the EDGE_VELOCITY is equal to zero but the FLOW_ACROSS_FACE is nonzero. I suspect this is due to the coarser resolution of RAS output than the actual model timestep; i.e., maybe EDGE_VELOCITY pulls from the instantaneous velocity whereas the FLOW_ACROSS_FACE is the aggregate value of any sub-second flows across the face in RAS. I observe this in Plan 37.

@jrutyna- do you know if this is the case?

Implications:

• Advection Coefficient: this is calculated as follows, which was a remnant of some old logic that actually exposed this issue. I suspect we'll just want the advection coefficient to be equal the flow across the face, and not set this to zero if the edge velocity is 0.

mesh[FLOW_ACROSS_FACE] * np.sign(abs(mesh[EDGE_VELOCITY]))

• Flow Across Face: if we do change the advection coefficient as proposed above, then we could delete this variable (since it should be identical to the advection coefficient) in the matrix and free up some memory.

For eventual use with pytest for automated testing.

Every boundary (ghost) cell should include a boundary condition. This issue can be resolved when there is a mechanism in place to throw an error if a boundary condition is not specified. The error should be informative, telling users which boundaries are missing.

cProfile results (profile.prof), explored in SnakeViz identify bottlenecks in the code. For clearwater-riverine, the main bottleneck is in

clearwater_riverine\transport.py:159(update), and specifically within the code that updates the left hand side matrix to solve the transport equation. Within that function, time is primarily sent in xarray.

This issue can be closed when some low-hanging performance gains are made.

Hypothesis
Variables are read from RAS HDF output to an xarray, but then those values must be put into a numpy format for compatibility with the sparse matrix solver. I think these data could be stored outside of an xarray and only saved to the xarray if requested by the user, which would reduce conversion in and out of xarray format.

• Explore how data is stored for variables in clearwater-modules and implement a similar approach.
• Engage with #57 for clearwater-modules on skipping writing intermediate arrays to xarray since this will be a way to speed up both of these models.

@aufdenkampe, @jrutyna - some questions about boundary conditions to keep in mind as you review the code.

Ghost Cells Receiving Flow
In RAS2D output, all ghost cells (boundary cells) have a volume of 0. However, this results in an error in the sparse matrix solver because nothing can leave the mesh, even when it should (e.g. in our box model.) Therefore, we must calculate the 'volume' in the ghost cells as the total flow across the edge between the real cell and ghost cell in a given timestep: implemented here. These ghost volumes are then inserted into the sparse LHS matrix on the diagonals (implemented here), so that material can move out of the grid when there is flow across an edge to a ghost cell.

Ghost Cells Flowing Into the Mesh
In order to set boundary conditions of a concentration of a certain contaminant/ghost cell, we also need to consider the volume of that cell when solving our sparse matrix. I calculate the volume of these ghost cells in a similar way as above, where I calculate teh flow across the edge between the real cell and the ghost cell in a given timestep, but now I assign that volume to the prior timestep for flow coming into the model mesh (implemented here).

Putting this information in the Sparse Matrix
I know I need to consider ghost cells flowing into the mesh on the right hand side of the sparse matrix problem in order to have any boundary condition concentrations flow into the mesh. But how should I handle the ghost cells receiving flow in this right hand side? Should the concentration reset to 0 at every timestep? Should users be defining ambient boundary conditions for every timestep? Currently, I include both of these ghost cell volumes in the [RHS class] (https://github.com/EnvironmentalSystems/ClearWater-riverine/blob/main/examples/ras2dwq.py#L648-L649)

TLDR

• How should boundary conditions be handled for ghost cells?
• Is there any difference in how boundary conditions should be handled for ghost cells that are flowing into the mesh vs receiving flow from the mesh?

Linking to TSM (and likely other clearwater modules) requires wetted surface area, which is not included in the RAS output. We need to:

• Add method to calculate wetted surface area for a mesh
• Make that an optional output when initializing (i.e., don't need it if just doing transport).

The model produces poor results when the Courant Condition is violated. HEC-RAS can run at a fine timestep where the timestep where the condition is not violated (e.g., 1 second) but print to a coarser timestep where the condition would be violated (e.g., 5 minutes). We should investigate:

1. If it can be read from the RAS HDF file (accurately to the output timestep, not the model timestep)
2. If not, if it can be calculated from output that is printed to the RAS HDF file
3. Or estimated in some other way

This issue can be resolved when there is a mechanism within Clearwater-Riverine to check if the Courant condition is violated, and throw a warning to the user if it is.

It will be fastest to run a number of constituents at once rather than re-running clearwater riverine multiple times. This issue can be closed when there is a functional Constituents class that manages each individual constituent.

We have increased the domain size of sumwere creek. Re-run the 100 everywhere test at 5, 10, 15 second timesteps and look at percent error in the output. This will help speed up demonstrations.

The HEC-RAS HDF file does not contain the coordinate system of the model grid. Users must therefore define the CRS in the plot() method. However, if a user has the model grid shapefile, it may be possible to parse the coordinate system from that file. We should explore the possibility of allowing users to pass the path to the RAS model shapefile to parse the CRS (in addition to allowing for them to manually define it as a string).

The _calc_ghost_cell_volumes() function in the utilities module has too much responsibility and we should seek to find ways to reduce the complexity. We could move some of the for loops into separate functions or find ways to combine them.

Wetting and drying cells leads to spikes in concentration. This issue can be resolved when we have developed a way to manage these spikes.

Thanks to Todd for digging in on the proper way to calculate cell volumes and face areas, originally outlined in this notebook and transferred to a python module here

• There appear to be a lot of for loops, so these functions should be reviewed to see if there is any opportunity for vectorization to improve performance.
• Explore best practice on function inputs (better to have fewer general inputs (e.g., one dataframe) or more explicit inputs (e.g., columns of that dataframe as lists).
• Figure out if face_areas_1 or face_areas_2 should be used
  • Detailed in: #4.

Users define initial and boundary conditions, which then must be properly converted for integration into the water quality simulation, which uses whatever units are used in the HEC-RAS model. Currently, the model assumes initial and boundary conditions are in the same units, and then users define their mass, volume, and unit conversions in the simulate_wq() method, after defining their initial and boundary conditions. Future versions of Clearwater Riverine should refine how these units are handled.

In resolving this issue we should do the following:

1. Define Conceptual Framework: how much freedom do we want to give users? Do we want to handle all unit conversions, or make users do the leg work on their own to input defined units? If we want to handle certain unit conversions, what are common units and conversions that we could contain within the module? Should we make users input initial and boundary conditions in the same units, or should we allow for different units and handle conversions? Answering these kinds of questions will involve consulting with water quality modelers to better understand how these situations are handled in other water quality models.
2. Move Unit Conversions to Appropriate Methods: Currently, any input units are defined and converted in the simulate_wq() method. We should consider moving this to the initial_conditions() and/or boundary_conditions(). The ultimate configuration will depend on the results of the conceptual framework defined in step 1.

We should add the following features to the Clearwater plotting methods:

• Play button
• Option to save to HTML or GIF embedded within the methods
• Allow users to pass additional plotting arguments (clim, cmap, height, width, etc.)
• Should we allow users to plot the variable of their choice (e.g., volume, load) rather than forcing them to plot concentration? It would have to be data at the cell level.
• Are there any other types of plots -- either maps (contour, rasterized) or pulling data (creating a timeseries from a point) -- that would be useful to integrate into the Clearwater Riverine package?

Need to integrate unit handling in the clearwater-riverine model with the post processing utilities created by @jrutyna.

• Current implementation takes user input units and converts them to model units for calculations, and then converts them back to whatever the user entered

Decision points

• Handling imperial units: convert imperial units to metric units right at the start?
• When and how to convert units from inputs (initial conditions and boundary conditions): require users to enter unit conversion factor when initializing model to an internal modeling unit (e.g., mg/L)
• Amount of flexibility to provide users (vs. deciding output units).

This issue can be closed when the decisions above have been made and the units in clearwater-riverine's transport module and postproc_utils modules are aligned.

The current Clearwater Riverine model uses Upwind Differencing. Users should be able to specify which advection scheme they want (linear, linear upwind, or upwind). See this online textbook for more details on these schemes and this video for a helpful overview.

This issue can be closed when users can specify an advection scheme when they run the model.

Currently, the hdf.py module (see line 30) is using the "Time Date Stamp" table from the hdf5 file. This is okay to use unless the HEC-RAS output is set to sub-second intervals then the "Time Date Stamp (ms)" table should be used. Need to investigate if this has any adverse impacts elsewhere in the repository before making this change.

As discussed in #49 here, converting from binary times in the RAS HDF file is a slow point in model instantiation.

A brief analysis found that leveraging pandas yielded faster results than np.vectorize(). Additional solutions should be assessed to help speed up this step.

In the utilities module, there are two functions (_compute_face_areas() and _compute_cell_volumes()) that are O(n^2), which could impact performance. These functions should be analyzed to see if we can achieve any gains in performance. There are also opportunities to refine these functions and/or split into smaller functions to reduce complexity and improve readability, since the nested logic within these functions makes it difficult to track, test, and debug.

Neither of these functions are used in current testing, because face areas and cell volumes are printed to the RAS output. We should develop a test HEC-RAS case that has these outputs turned off so that we can compare performance with and without these outputs.

Continue the work in #14, but now for improving the handling of reading other types of inputs:

• Initial Conditions
• Boundary Conditions

All examples we are working with (Box model, Muncie, Ohio River) have a single 2D flow area. The name of that flow area is parsed using the parse_project_name function, which simply grabs the first 2D Flow Area name:

def parse_project_name(infile: h5py._hl.files.File) -> str:
    '''Parse the name of a project'''
    project_name = infile['Geometry/2D Flow Areas/Attributes'][()][0][0].decode('UTF-8')
    return project_name

@jrutyna pointed out that RAS models can have multiple 2D flow areas within a single project. The code is not currently configured to handle this kind of situation. We do not need to address this in the short term, but we will eventually need to figure out how we want to handle projects with multiple flow areas (e.g., looping though all of them, allowing users to select a single area, some combination of those options, etc.).

In the current Clearwater Riverine version, CSV files for boundary conditions must contain very specific information: the boundary condition name (as named in the HEC-RAS model) and the concentration of that boundary cell for each timestep in the HEC-RAS model output. If a boundary cell / timestep combo is not in the CSV, the concentration is assumed to be 0. However, users would likely find some additional capabilities useful:

• Boundary condition interpolation: users will likely have water quality data that is sporadic, and that may not line up exactly with the HEC-RAS timesteps. It would therefore be useful for this package to contain methods that could handle these scenarios and interpolate water quality to align with the appropriate modeling timestep.
• Handling spatial data: Allow users to input monitoring locations (lat/long or shapefile) and develop logic to align that information with the HEC-RAS boundary cells.

Issue #6 identified an issue with the way HEC-RAS mapped boundary condition lines to external faces in its hdf5 results file; a bug report has been sent to the HEC-RAS development team regarding this issue. A work around solution was implemented for issue #6 for the example simulation comparison between an EFDC model of the Ohio River and the same system in HEC-RAS/Clear Water.

A more robust solution will be needed to handle these problematic boundary condition lines if the HEC-RAS development team does not resolve the issue or before they release the update. The two proposed solutions that @jrutyna identified in issue #6 that could be implemented within the Clear Water framework are:

1. Develop a similar geo-processing tool that HEC-RAS performs to identify external faces from the boundary condition line, or
2. Develop a method from the face flow results to identify which external faces have flows coming in/out of the model domain and use these faces as a lookup in the Geometry/Boundary Condition Lines/External Faces table

@TSteissberg - FYI, making this a formal issue so we can track our process and decision making.

Edge vertical areas can be calculated in multiple ways - either using information from the first cell or the second cell in the "Faces Cell Index" dataset, as outlined in this notebook. One could also potentially use alternative methods (e.g., a weighted average of the values from these two cells).

We will need to review all of these potential choices, and possibly conduct a sensitivity analysis to see how different choices impact the results.

Add pre-commit and update repository to pass all pre-commits.

This is calculated using the _calc_sum_coeff_to_diffusion_term function and _sum_vals, which currently loops over time. As the number of timesteps is increasing, this is becoming extremely slow (>3 min to instantiate a model). We should investigate how to vectorize this calculation over the time dimension.

The plot() method of the ClearwaterRiverine class throws a ShapelyDeprecationWarning. We currently ignore that warning with warnings.filterwarnings("ignore", category=ShapelyDeprecationWarning), but we should ultimately work to understand and resolve the underlying root cause.

It makes sense to work toward similar if not identical development environments for both ClearWater-Riverine and ClearWater-modules-python. This will help with this next phase of integration, even though it will be less important once we properly package each library.

Thanks for the suggestion @jrutyna!

cc: @xaviernogueira, @sjordan29

Matrices can be reduced in size when there are dry cells to potentially reduce memory and solve time.

Currently, if there is a dry cell, a placeholder value is placed in the diagonal of the matrix so that the matrix can be inverted. However, a more elegant solution would be to remove dry (unnecessary) rows/columns from the matrix entirely. We would need to assign the matrix row and column indices to their corresponding cell indices in the mesh and reassign the solved concentration values to the proper cells after solving.

Furthermore, if a cell is dry over the entire simulation horizon, we could potentially remove it from the simulation entirely and free up memory (both in the matrix and in the model mesh itself) to help with #61.

@aufdenkampe, @jrutyna, @xaviernogueira - as I'm thinking about #43, I was wondering if it might be useful to leverage a config file for Clearwater-Riverine. I feel like this approach could have a few advantages:

• Easy to set up for modelers who might not have as much Python experience
• Simple way to turn on / off modules (are they included in the config or not?)
• Easier to add flexibility in the number of general constituents to include
• More transparency on what information is needed for each of the modules from a user standpoint.

I'm envisioning clearwater looking for a config.yml file that looks something like this, which would be the only input for instantiating a Clearwater-Riverine model:

- name: Sumwere Creek Plan 10
- description: Demo model
- path: ./sumwere_creek_p10.hdf
- diffusion_coefficient: 0.1
- modules:
  - gsm:
    - constituent-A:
      - initial_conditions: ./constituent_a_initial_conditions.csv
      - boundary_conditions: ./constituent_a_boundary_conditions.csv
      - order: 0
    - constituent-B:
      - initial_conditions: ./constituent_b_initial_conditions.csv
      - boundary_conditions: ./constituent_b_boundary_conditions.csv
      - order: 0
  - tsm:
    - initial_conditions: ./tsm_initial_conditions.csv
    - boundary_conditions: ./tsm_boundary_conditions.csv
    - meteo_parameters: ./tsm_meteo_parameters.csv
    - temp_parameters: ./tsm_temp_parameters.csv

To me, this feels like a more organized approach for setting up a Clearwater Riverine model than having a ton of optional inputs when instantiating the model, but I'm curious to hear what others thing from both a programming and modelling standpoint.

When a new user tries to clone the entire repository onto their machine, we run into a git lfs bandwidth error. We need to remove any HDF file that is not used for unit testing.

@aufdenkampe - can you link to the steps you took to delete certain large files from the entire git history when splitting this repository from the original Clearwater repo?

It appears that the PyData/sparse library more natively integrates with Xarray and Dask than the scipy.sparce subpackage, because it was designed from the ground up to natively use arrays instead of matrices. See:

• https://sparse.pydata.org/en/stable/roadmap.html
• https://tutorial.xarray.dev/overview/xarray-in-45-min.html#alternative-array-types
• https://docs.dask.org/en/latest/array-sparse.html

It's possible that this is being solved by the recent refactoring of scipy.sparce, but Xarray and Dask are already integrated with PyData/sparse.

It might be worth exploring if migrating to PyData/sparse would provide any benefits to performance or code simplicity.

On the other hand, development doesn't appear very active, so it's possible that momentum has shifted away from this effort. See: https://sparse.pydata.org/en/stable/changelog.html

Check unit conversions for mg/L inputs.

This issue can be closed when inputting 100 mg/L as units for mass balance test and ~100 mg/L pops out as the result.