Skip to content

Commit

Permalink
feat(PRT): introduce Particle Tracking (PRT) model (#1389)
Browse files Browse the repository at this point in the history
At this time PRT supports particle tracking on all DIS and some but not all DISV grids. DISV cells containing adjacent parallel faces are not properly handled and may cause crashes, this includes quad-refined grids as well as DISV grids generated by FloPy's voronoi module.

A GWF-PRT exchange is included, which can be used to couple GWF and PRT models in the same simulation. PRT can also consume flows from a previously run GWF model via FMI.

Support is planned but not yet implemented for:

* configurable particle mass
* local mass conservation — conserved between cells but not yet guaranteed at the subcell level
* exchanging particles between models
* DISU grids (longer term)
* solving models in parallel (longer term)

Docs/dfns and tests make up about half the changeset, the rest is fortran source. Input/output format is likely to change without notice until PRT reaches a release-ready state.

PRT's primary author is Alden Provost. Some code is adapted from MODPATH and originally due to David Pollock.

Co-authored-by: Alden Provost <[email protected]>
  • Loading branch information
wpbonelli and aprovost-usgs authored Feb 23, 2024
1 parent ab6d3f8 commit 61b0883
Show file tree
Hide file tree
Showing 150 changed files with 19,759 additions and 151 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
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 61b0883

Please sign in to comment.