Hydrologic features

[piyueh]

Removing working fluid when encountering hydrological features

If the working fluid is not water, then when it flows into water bodies/hydrological features, another solver will be used for this hydrographic transport. And the working fluid flowing into water bodies requires to be removed from the overland flow solver, which in this case is GeoClaw. This PR implements the following items to achieve this:

1. Reading hydrological features: Users will have to provide raster files of the hydrological features. The number of files can be arbitrary, but currently, only Esri ASCII format is accepted. The setup is through the HydroFeatureData object in setrun.py. For example (Note, in the future, after merging back to GeoClaw and Clawpack repos, the registration of HydroFeatureData will happen in the initialization of rundata, so users won't have to do import and rundata.add_data manually in setrun.py.):

   from clawpack.geoclaw.data import HydroFeatureData                          
   rundata.add_data(HydroFeatureData(), 'hydro_feature_data')                  
   hydro_feature_data = rundata.hydro_feature_data                             
   hydro_feature_data.files.append("./hydro_feature1.asc")                     
   hydro_feature_data.files.append("./hydro_feature2.asc")                     
   hydro_feature_data.files.append("./hydro_feature3.asc")                     

2. Identifying which cells belong to hydro features: I store this information in aux array. The setaux.f90 is in charge of initializing this information whenever there's a new grid created. The tricky part is that users will have to set up the num_aux correctly in their setrun.py. For the sake of defensive programming, the solver will check if the number of aux datasets is correct or not at the initialization stage of the hydro feature class.

3. Removing the fluid flowing into hydro features: this is done in src2.f90.

4. A new example for hydrological features: See examples/landspill/utah_dem_hydro_feat.

Miscellaneous

Some items/bugs unrelated to hydrological features are also implemented/fixed in this PR.

1. Wrong order of reading each row in ESRI ASCII format.
2. Using pyclaw.Solution object in the post-processing Python scripts for the solution I/O.
3. Changing the output file format to binary.
4. A new post-processing script to create a single NetCDF file for GIS software: writenc.py.

[piyueh]

@mandli This PR is ready for review. The Docker image: https://hub.docker.com/r/barbagroup/landspill/ with tag PR10. Thank you so much!!

[piyueh]

@mandli More commits have been added to this PR. The new commits are about tracing the locations where the working fluid is removed and the volumes removed from these locations.

I didn't use the fixed grid monitors or gauges. What I want to record is the total volumes removed, so this requires some modification in the original GeoClaw source code, which is I try to avoid.

The tracer implemented in this PR is a sparse matrix in which the non-zero elements represent the locations of the boundaries of hydro features. The module creates a virtual Cartesian grid with domain size equal to the computational domain and with cell size equal to the finest cell size in AMR setting. A 2D Cartesian grid can be represented as a matrix. This matrix records the volume removed from each cell. The removal of fluid can only happen at the boundaries of hydrological features, so very few elements in this matrix (i.e., cells in this virtual grid) are non-zero. And hence I use a sparse matrix to represent this grid and store the removed volumes at hydrological boundaries.

[piyueh]

Thanks for your review, @mandli! Since this PR is not blocking my current development, we can keep it open and keep finding a new name for the functionality.

[piyueh]

Thank you @mandli! Yes, I believe renaming should be easy.