Skip to content

Commit

Permalink
feat(PRT): add particle tracking model
Browse files Browse the repository at this point in the history
Co-authored-by: Alden Provost <[email protected]>
  • Loading branch information
wpbonelli and aprovost-usgs committed Feb 23, 2024
1 parent fa29361 commit e7315bd
Show file tree
Hide file tree
Showing 153 changed files with 19,756 additions and 153 deletions.
159 changes: 159 additions & 0 deletions .github/common/update_fortran_style.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,159 @@
import argparse
import re
from contextlib import nullcontext
from itertools import repeat
from pathlib import Path
from typing import Iterator, Optional
from warnings import warn

from fprettify.fparse_utils import InputStream

INTENT_PATTERN = re.compile(r".*(intent\(.+\)).*")


def get_intent(s) -> Optional[str]:
result = INTENT_PATTERN.match(s)
return result.group(1) if result else None


def get_param(s) -> bool:
return "parameter" in s


def get_comments(comments) -> Iterator[str]:
for comment in comments:
if not any(comment):
continue
yield comment.rstrip()


class Transforms:
@staticmethod
def separate_lines(path, overwrite=False):
"""Variables defined on separate lines"""

flines = []
with open(path, "r") as f:
stream = InputStream(f)
while 1:
line, comments, lines = stream.next_fortran_line()
if not lines:
break
line = line.rstrip()
parts = line.rpartition("::")
comments = " " + "".join(get_comments(comments))
if not parts[1] or "procedure" in parts[0]:
for l in lines:
flines.append(l.rstrip())
continue

nspaces = len(lines[0]) - len(lines[0].lstrip())
prefix = "".join(repeat(" ", nspaces))
vtype = parts[0].split(",")[0].strip()
split = parts[2].split(",")
intent = get_intent(parts[0])
param = get_param(parts[0])

if not line:
continue
if (len(parts[0]) == 0 and len(parts[1]) == 0) or (
"(" in parts[2] or ")" in parts[2]
):
flines.append(prefix + line + comments)
elif len(split) == 1:
flines.append(prefix + line + comments)
elif param:
flines.append(prefix + line + comments)
else:
for s in split:
if s.strip() == "&":
continue
l = prefix + vtype
if intent:
l += f", {intent}"
l += f" :: {s.strip()}"
flines.append(l + comments)

with open(path, "w") if overwrite else nullcontext() as f:

def write(line):
if overwrite:
f.write(line + "\n")
else:
print(line)

for line in flines:
write(line)

@staticmethod
def no_return_statements(path, overwrite=False):
"""Remove return statements at the end of routines"""
# todo
pass

@staticmethod
def no_empty_comments(path, overwrite=False):
"""Remove comments on lines with only whitespace"""
# todo
pass


def reformat(path, overwrite, separate_lines, no_return_statements, no_empty_comments):
if separate_lines:
Transforms.separate_lines(path, overwrite=overwrite)
if no_return_statements:
Transforms.no_return_statements(path, overwrite=overwrite)
warn("--no-return not implemented yet")
if no_empty_comments:
Transforms.no_empty_comments(path, overwrite=overwrite)
warn("--no-empty-comments not implemented yet")


if __name__ == "__main__":
parser = argparse.ArgumentParser(
"""
Modify MODFLOW 6 Fortran source code, either writing to stdout or
overwriting the input file. Options are provided for several code
styles.
"""
)
parser.add_argument(
"-i", "--input", help="path to input file" # todo: or directory
)
parser.add_argument(
"-f",
"--force",
action="store_true",
default=False,
required=False,
help="overwrite/reformat files",
)
parser.add_argument(
"--separate-lines",
action="store_true",
default=True,
required=False,
help="define variables on separate lines",
)
parser.add_argument(
"--no-return_statements",
action="store_true",
default=False,
required=False,
help="no return statements at the end of routines",
)
parser.add_argument(
"--no-empty-comments",
action="store_true",
default=False,
required=False,
help="no empty comments",
)
args = parser.parse_args()
reformat(
path=Path(args.input).expanduser().absolute(),
overwrite=args.force,
separate_lines=args.separate_lines,
no_return_statements=args.no_return_statements,
no_empty_comments=args.no_empty_comments,
)
6 changes: 6 additions & 0 deletions .github/workflows/ci.yml
Original file line number Diff line number Diff line change
Expand Up @@ -336,6 +336,12 @@ jobs:
with:
repository: MODFLOW-USGS/modflow6-testmodels
path: modflow6-testmodels

- name: Checkout modflow6-examples
uses: actions/checkout@v4
with:
repository: MODFLOW-USGS/modflow6-examples
path: modflow6-examples

- name: Setup Micromamba
uses: mamba-org/setup-micromamba@v1
Expand Down
2 changes: 1 addition & 1 deletion CITATION.cff
Original file line number Diff line number Diff line change
Expand Up @@ -55,7 +55,7 @@ authors:
alias: mjreno
- family-names: Bonelli
given-names: Wesley P.
alias: w-bonelli
alias: wpbonelli
affiliation: U.S. Geological Survey
orcid: https://orcid.org/0000-0002-2665-5078
- family-names: Boyce
Expand Down
4 changes: 3 additions & 1 deletion README.md
Original file line number Diff line number Diff line change
Expand Up @@ -67,12 +67,14 @@ Instructions for building definition files for new packages are summarized in [d

MODFLOW is a popular open-source groundwater flow model distributed by the U.S. Geological Survey. For 30 years, the MODFLOW program has been widely used by academic researchers, private consultants, and government scientists to accurately, reliably, and efficiently simulate groundwater flow. With time, growing interest in surface and groundwater interactions, local refinement with nested and unstructured grids, karst groundwater flow, solute transport, and saltwater intrusion, has led to the development of numerous MODFLOW versions. Although these MODFLOW versions are often based on the core version (presently MODFLOW-2005), there are often incompatibilities that restrict their use with one another. In many cases, development of these alternative versions has been challenging due to the underlying MODFLOW structure, which was designed for the simulation with a single groundwater flow model using a rectilinear grid.

MODFLOW 6 is the latest core version of MODFLOW. It synthesizes many of the capabilities available in MODFLOW-2005, MODFLOW-NWT, and MODFLOW-LGR. MODFLOW 6 was built on a new object-oriented framework that allows new packages and models to be added, and allows any number of models to be run simultaneously in a single simulation. Model may be coupled sequentially, such as for flow and transport, or the models may be tightly coupled at the matrix level, such as for multiple flow models. MODFLOW 6 presently contains two types of hydrologic models, the Groundwater Flow (GWF) Model and the Groundwater Transport (GWT) Model.
MODFLOW 6 is the latest core version of MODFLOW. It synthesizes many of the capabilities available in MODFLOW-2005, MODFLOW-NWT, and MODFLOW-LGR. MODFLOW 6 was built on a new object-oriented framework that allows new packages and models to be added, and allows any number of models to be run simultaneously in a single simulation. Model may be coupled sequentially, such as for flow and transport, or the models may be tightly coupled at the matrix level, such as for multiple flow models. MODFLOW 6 presently contains three types of hydrologic models, the Groundwater Flow (GWF) Model, the Groundwater Transport (GWT) Model, and the Particle Tracking (PRT) Model.

The Groundwater Flow (GWF) Model was the first model to be released in MODFLOW 6. It supports regular MODFLOW grids consisting of layers, rows, and columns, but it also supports more flexible grids that may conform to irregular boundaries or have increased resolution in areas of interest. The GWF Model consists of the original MODFLOW stress packages (CHD, WEL, DRN, RIV, GHB, RCH, and EVT) and four advanced stress packages (MAW, SFR, LAK, and UZF), which have been distilled from their predecessors to contain the most commonly used capabilities. MODFLOW 6 contains a new Water Mover (MVR) Package that can transfer water from provider packages to receiver packages. Providers can be many of the stress and advanced stress packages; receivers can be any of the advanced stress packages. This new capability makes it possible to route water between lakes and streams, route rejected infiltration into a nearby stream, or augment lakes using groundwater pumped from wells, for example. To modernize user interaction with the program, the MODFLOW 6 input structure was redesigned. Within package input files, information is divided into blocks, and informative keywords are used to label numeric data and activate options. This new input structure was designed to make it easier for users to adjust simulation options in an intuitive manner, reduce user input errors, and allow new capabilities to be added without causing problems with backward compatibility.

The GWT model for MODFLOW 6 simulates three-dimensional transport of a single solute species in flowing groundwater. The GWT Model solves the solute transport equation using numerical methods and a generalized CVFD approach, which can be used with regular MODFLOW grids or with unstructured grids. The GWT Model is designed to work with most of the new capabilities released with the GWF Model, including the Newton flow formulation, unstructured grids, advanced packages, and the movement of water between packages. The GWF and GWT Models operate simultaneously during a MODFLOW 6 simulation to represent coupled groundwater flow and solute transport. The GWT Model can also run separately from a GWF Model by reading the heads and flows saved by a previously run GWF Model. The GWT model is also capable of working with the flows from another groundwater flow model, as long as the flows from that model can be written in the correct form to flow and head files.

The Particle Tracking (PRT) Model simulates three-dimensional particle trajectories in flowing groundwater. The PRT Model can operate simultaneously with a GWF model via an exchange, or can consume GWF outputs via Flow Model Interface (FMI). The PRT Model solves structured DIS grids analytically and unstructured DISV grids semi-analytically. Tracking delegates from the model domain to individual cells, and to sub-components of cells for DISV grids. On structured grids the approach reduces to Pollock's method used in MODPATH 7. On DISV grids, polygonal cells are decomposed into triangles, within which particle exit faces/times are solved numerically. Track data are recorded at the boundaries between spatial volumes (cells or subcells) and time segments (timesteps or stress periods), and at other relevant times, e.g. release and termination. Events to record may be configured by the user. Though each particle's motion may be computed independently, parallel solving is not yet supported. Particle exchange between models is also planned but not yet supported. Particle mass is conserved at the cell level, but may not be conserved at the subcell level for DISV grids.


## How to Cite MODFLOW 6

Expand Down
1 change: 1 addition & 0 deletions autotest/TestGeomUtil.f90
Original file line number Diff line number Diff line change
Expand Up @@ -292,6 +292,7 @@ subroutine test_point_in_polygon_irr(error)
deallocate (face_pts)
end subroutine test_point_in_polygon_irr

!> @brief Test 2D skew
subroutine test_skew(error)
type(error_type), allocatable, intent(out) :: error
real(DP) :: v(2)
Expand Down
1 change: 0 additions & 1 deletion autotest/TestList.f90
Original file line number Diff line number Diff line change
Expand Up @@ -112,7 +112,6 @@ subroutine test_get_next_previous_item_reset(error)
type(ListType), pointer :: list
type(IntNodeType), pointer :: n1, n2, n3
class(*), pointer :: p
integer(I4B) :: i

allocate (list)
allocate (n1)
Expand Down
44 changes: 22 additions & 22 deletions autotest/TestMathUtil.f90
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@ module TestMathUtil
use testdrive, only: check, error_type, new_unittest, test_failed, &
to_string, unittest_type
use MathUtilModule, only: f1d, is_close, mod_offset, &
zeroch, zerotest, zeroin
zero_ch, zero_test, zero_br
implicit none
private
public :: collect_mathutil
Expand All @@ -19,12 +19,12 @@ subroutine collect_mathutil(testsuite)
test_is_close_symmetric_near_0), &
new_unittest("mod_offset", &
test_mod_offset), &
new_unittest("zeroch", &
test_zeroch), &
new_unittest("zeroin", &
test_zeroin), &
new_unittest("zerotest", &
test_zerotest) &
new_unittest("zero_ch", &
test_zero_ch), &
new_unittest("zero_br", &
test_zero_br), &
new_unittest("zero_test", &
test_zero_test) &
]
end subroutine collect_mathutil

Expand Down Expand Up @@ -177,67 +177,67 @@ pure function sine(bet) result(s)
s = sin(bet)
end function sine

subroutine test_zeroch(error)
subroutine test_zero_ch(error)
type(error_type), allocatable, intent(out) :: error
real(DP), parameter :: pi = 4 * atan(1.0_DP)
real(DP) :: z
procedure(f1d), pointer :: f

f => sine

z = zeroch(-1.0_DP, 1.0_DP, f, 0.001_DP)
z = zero_ch(-1.0_DP, 1.0_DP, f, 0.001_DP)
call check(error, is_close(z, 0.0_DP, atol=1d-6), &
'expected 0, got: '//to_string(z))

z = zeroch(-4.0_DP, -1.0_DP, f, 0.001_DP)
z = zero_ch(-4.0_DP, -1.0_DP, f, 0.001_DP)
call check(error, is_close(z, -pi, atol=1d-6), &
'expected -pi, got: '//to_string(z))

z = zeroch(1.0_DP, 4.0_DP, f, 0.001_DP)
z = zero_ch(1.0_DP, 4.0_DP, f, 0.001_DP)
call check(error, is_close(z, pi, atol=1d-6), &
'expected pi, got: '//to_string(z))
end subroutine test_zeroch
end subroutine test_zero_ch

subroutine test_zeroin(error)
subroutine test_zero_br(error)
type(error_type), allocatable, intent(out) :: error
real(DP), parameter :: pi = 4 * atan(1.0_DP)
real(DP) :: z
procedure(f1d), pointer :: f

f => sine

z = zeroin(-1.0_DP, 1.0_DP, f, 0.001_DP)
z = zero_br(-1.0_DP, 1.0_DP, f, 0.001_DP)
call check(error, is_close(z, 0.0_DP, atol=1d-6), &
'expected 0, got: '//to_string(z))

z = zeroin(-4.0_DP, -1.0_DP, f, 0.001_DP)
z = zero_br(-4.0_DP, -1.0_DP, f, 0.001_DP)
call check(error, is_close(z, -pi, atol=1d-6), &
'expected -pi, got: '//to_string(z))

z = zeroin(1.0_DP, 4.0_DP, f, 0.001_DP)
z = zero_br(1.0_DP, 4.0_DP, f, 0.001_DP)
call check(error, is_close(z, pi, atol=1d-6), &
'expected pi, got: '//to_string(z))
end subroutine test_zeroin
end subroutine test_zero_br

subroutine test_zerotest(error)
subroutine test_zero_test(error)
type(error_type), allocatable, intent(out) :: error
real(DP), parameter :: pi = 4 * atan(1.0_DP)
real(DP) :: z
procedure(f1d), pointer :: f

f => sine

z = zerotest(-1.0_DP, 1.0_DP, f, 0.001_DP)
z = zero_test(-1.0_DP, 1.0_DP, f, 0.001_DP)
call check(error, is_close(z, 0.0_DP, atol=1d-6), &
'expected 0, got: '//to_string(z))

z = zerotest(-4.0_DP, -1.0_DP, f, 0.001_DP)
z = zero_test(-4.0_DP, -1.0_DP, f, 0.001_DP)
call check(error, is_close(z, -pi, atol=1d-6), &
'expected -pi, got: '//to_string(z))

z = zerotest(1.0_DP, 4.0_DP, f, 0.001_DP)
z = zero_test(1.0_DP, 4.0_DP, f, 0.001_DP)
call check(error, is_close(z, pi, atol=1d-6), &
'expected pi, got: '//to_string(z))
end subroutine test_zerotest
end subroutine test_zero_test

end module TestMathUtil
Loading

0 comments on commit e7315bd

Please sign in to comment.