Skip to content

Conversation

@piyueh
Copy link

@piyueh piyueh commented Apr 11, 2019

Thank you all for developing this awesome software! We expanded GeoClaw to simulate the hydrocarbon overland flow released from ruptured pipelines. All extensions and modifications are located inside the folder src/2d/shallow/landspill. None of the original files in src/2d/shallow is modified, except the changes made in the commit 042d04d. In this way, we can minimize the possibility of introducing new bugs into GeoClaw. Also, examples are in examples/landspill.

All source code of our work had been reviewed by @mandli at our fork of GeoClaw. See the following PRs in the forked repository: PR8, PR9, PR10, PR14, PR16, and PR17. We use object-oriented design to make things more modularized.

The following is an overview of our work:

1. Point sources

The releasing points of fluid from pipelines are modeled as point sources. The code is inside src/2d/shallow/landspill/point_source. A point source is represented by PointSource class. And a collection of point sources is represented by PointSourceCollection. We simply add point sources to the right-hand-side of the mass equation in src2.f90.

See PR8 for more details.

2. Darcy-Weisbach friction models

We prefer Darcy-Weisbach friction models over Manning's. The source code is in src/2d/shallow/landspill/darcy_weisbach. There are many different models to calculate the friction coefficient in Darcy-Weiscath friction. So we use an abstract class, DarcyWeisbachBase, to rule how all different models should be implemented and follow a unified interface. All models for calculating coefficients should be derived from this abstract class. Currently, we have implemented 7 different models: DarcyWeisbachNull, DarcyWeisbachConstant, DarcyWeisbachBlockConstants, DarcyWeisbachCells, DarcyWeisbachChurchill, DarcyWeisbachTwoRegimes, and DarcyWeisbachThreeRegimes.

The abstract-class design requires using pointers in the main function. But pointers seem to be not welcomed in Fortran, so we use another class, DarcyWeisbach, to encapsulate the pointers. And users of this class will not have to deal with pointers in their main Fortran code.

The Darcy-Weisbach friction model replaces Manning's model in src2.f90. See PR9 for more information.

3. Intersection with On-land hydrological features

When hydrocarbon fluid contact with on-land hydrological features (e.g., rivers, lakes, ponds, etc.), we need to record the location of contact lines, contact times, and the volumes flowing into these water bodies. And after an overland flow simulation, we can use these data as the input to other solvers that simulate hydrological transport. The source code is in src/2d/shallow/landspill/hydrologic.

Users provide hydrological data through ESRI ASCII files, and one ASCII file is represented by an instance of HydroFeature class. The class, HydroFeatureCollection, represents a collection of several files. Our implementation uses a concept similar to CSR sparse-matrix format to describe the locations of water bodies. Cells with zero values are land while those with values greater than zero are water bodies. Different values may refer to different types of water bodies. Now, all types of water bodies are handled in the same way, but in the future, different types of water bodies may be processed in different ways.

See PR10.

4. Temperature-dependent viscosity with the Lewis-Squires viscosity correction

Many models for the Darcy-Weisbach coefficients require the viscosity of the working fluid. We use the Lewis-Squire viscosity correction to get temperature-dependent viscosity. Currently, however, due to that ambient temperature and fluid temperature do not change during a simulation, so the viscosity is actually a constant from the beginning to the end of a simulation. We determine the viscosity based on the ambient temperature provided by users.

See PR14.

5. Evaporation models

We again utilize abstract-class design to implement evaporation models. Currently, there are only two models: Fingas' log model and Fingas' square root model, which are implemented in the classes EvapFingas1996Log and EvapFingas1996SQRT, respectively. The class EvapNull is used when evaporation is disabled, and EvapBase is the abstract class. Again, we use a top-level class EvapModel to encapsulate pointers. We hope this kind of design can make it easier to add more models in the future.

We treat the evaporation as a source (a negative source) in the right-hand side of the mass equation. And this is done in src2.f90.

See PR16 for more information.

6. Miscellaneous and minor optimization

Two slight optimizations are done in PR17. One is to ignore totally dry grid patches during time-marching. The other one is just to remove unnecessary (to oil spill applications) refining criteria in flag2refinre2.f90 so that the code can be cleaner for reading. The optimized version of the two source files is located in src/2d/shallow/landspill/optimized.

The only two modifications in the original GeoClaw source code are also in this PR (see 042d04d). The first one adds a numerical parameter to update.f90. This parameter controls how we update a coarse cell when the averaged fluid level of its child cells is below the coarse cell's topographical elevation. Originally in GeoClaw, when the averaged fluid level of child cells is lower than the parent cell's elevation, the solver chooses to maintain the steady-state condition rather than chooses to conserve the mass, so the parent cell will have zero fluid volume and will lose some fluid mass. But this causes some trouble for overland flows when there are canyon-like features in the topography. See the tests and the short report in this repository. After discussing this issue with @mandli, we decided to add a parameter to control the fluid volume of a coarse cell in this situation. The default value of this parameter is zero, and this will behave like the original GeoClaw. In our overland-flow applications, we set this parameter to dry_tol so that a coarse cell will always have some fluid in it as long as one of its child cells has fluid. This may break the steady-state condition, but it saves us from losing significant fluid mass in our applications.

The other change to GeoClaw is another numerical parameter in flag2refine2.f90. This parameter controls the criteria to refine a cell according to the fluid depth. The default value of this parameter is dry_tol, and it behaves like the original GeoClaw. In our applications, we'll set this parameter to zero, and the finest mesh will be applied to wherever there's fluid. The area of overland flow is usually small compared to the computational domain, so using the finest mesh to cover the whole flow area is computationally acceptable.

7. Examples

There are 5 examples along with their README in examples/landspill:

The first one, inclined_plane is a validation test. We compare the result with the experimental data from Lister (1992). This test takes a very long time to run due to the very fine mesh. Be cautious.

utah_dem and utah_dem_hydro_feat simulate the overland flow of Maya crude somewhere close to Salt Lake City. utah_dem represents a case where no on-land hydrological features present, while utah_dem_hydro_feat has a pond near the rupture point.

utah_dem_opt and utah_dem_hydro_feat_opt are the same as utah_dem and utah_dem_hydro_feat, except they use optimized version of update.f90 and flag2refine2.f90.

piyueh and others added 30 commits September 12, 2018 19:12
- initialize landspill module and set_landspill
- python interface for specifying point source data in setrun.py
- use instrinsic function new_line in generic write
- add set_landspill to amr2.f90
- modify Makefile.geoclaw
- add a function to apply point sources
- fix a bug in cell_id subroutine of point source instances
- add function call for applying point sources to src2.f90
This type calculate the friction coefficients on the fly based on local
Reynolds number and three models corresponding to laminar, smooth
turbulent, and fully rough turbulent regimes.

- modify the signature of get_coefficient member function
- add kinematic viscosity to land_spill module
- add the 4th derived Darcy-Weisbach implementation
…blems

- remove misused idint function
- intrinsic math function: replace single precision version with double
  precision version
- rename three_models with three_regimes
- fix the wrong signature of get_coefficient in DarcyWeisbach class
- a new class PointSourceCollection to encapsulate point sources
- rename class point_source -> PointSource
- rewrite I/O interface
- put point-source-related source files to a sub-folder
- calculation of Reynolds number
- definition of starting Re for full turbulence
- DarcyWeisbachThreeModels => DarcyWeisbachThreeRegimes
- semi-implicit firction term
- add flexibility for FFLAGS in the Makefile of the example
- remove `make .data` and `make .exe` as `make .output` will do that
- add python 2 support to download_topo.py and fix the bug of
  `/usr/bin/env`
- set up module_setup to true after module initialization.
- replace `write(*, *)` with `print *,`
- check if the list of end times provided by users is sorted or not in
  `data.py`
- explicitly triger refinement around point sources when t = t0
- reduce the initial depth at the location of the point sources ot
  0.9*dry_tolerance
piyueh and others added 21 commits January 15, 2019 18:29
Capability of setting extra physical properties
Landspill module and its objects do not support lat-long coordinate do
far. But when these objects are null objects (e.g., zero point source,
zero hydrological feature, etc), lat-long coordinate should still be
allowed.
If a point source is located at the boundary grid line of a grid patch,
rounding errors may cause the cell_id function not able to find correct
cell indices of the point source. In this fix, we add a tolerance to the
calculation.
* add settings of hidden parameters to setrun.py
* modify Makefile files
* revise README.md files
* add a new example utah_dem_hydro_feat_opt
Optimization and landspill-specific modification
fix #18: row indices not initialized correctly
remove dockerfiles for PR reviewers
@mandli mandli added this to the 5.7.0 milestone Aug 16, 2019
@mandli mandli requested a review from rjleveque August 29, 2019 23:48
@mandli
Copy link
Member

mandli commented Feb 20, 2020

@piyueh in case you are wondering we are looking at this right now.

@piyueh
Copy link
Author

piyueh commented Feb 20, 2020

Thank you very much!

@piyueh
Copy link
Author

piyueh commented Dec 10, 2020

@mandli Hi, I'm thinking about re-working on this PR:

  1. Rebase the branch onto the latest commit of the master branch (or onto v5.7.1), so the conflict can be resolved when merging
  2. Squash minor commits to reduce number of commits and make the commit history more clear
  3. move the big examples to the apps repository

Is it OK if I close this PR and make a new one later? Thanks!

@mandli
Copy link
Member

mandli commented Dec 10, 2020

@piyueh that sounds good. Sorry we never got back to this!

@piyueh piyueh closed this Dec 10, 2020
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants