diff --git a/autotest/prt/prt_test_utils.py b/autotest/prt/prt_test_utils.py index e61805325f0..9d45895faa0 100644 --- a/autotest/prt/prt_test_utils.py +++ b/autotest/prt/prt_test_utils.py @@ -1,20 +1,24 @@ import os from types import SimpleNamespace -from typing import Optional, Tuple, Union +from typing import Optional, Tuple import flopy import numpy as np -import pandas as pd +import matplotlib as mpl -def all_equal(col, val): - a = col.to_numpy() +def all_equal(series, val): + a = series.to_numpy() return a[0] == val and (a[0] == a).all() def get_gwf_sim( name, ws, mf6 ) -> Tuple[flopy.mf6.MFSimulation, SimpleNamespace]: + """ + Simple GWF simulation for use/modification by PRT tests + """ + # test case context ctx = SimpleNamespace( nlay=1, @@ -30,15 +34,12 @@ def get_gwf_sim( # mp7 release points (cell-local coordinates) releasepts_mp7=[ # node number, localx, localy, localz - # (0-based indexing converted to 1-based for mp7 by flopy) (0, float(f"0.{i + 1}"), float(f"0.{i + 1}"), 0.5) for i in range(9) ], - # expected release points in PRT format; use flopy - # to convert mp7 to prt format then check equality + # PRT release points (cell indices and global coordinates both required) releasepts_prt=[ # particle index, k, i, j, x, y, z - # (0-based indexing converted to 1-based for mf6 by flopy) [i, 0, 0, 0, float(f"0.{i + 1}"), float(f"9.{i + 1}"), 0.5] for i in range(9) ], @@ -119,7 +120,7 @@ def check_track_data( track_hdr: os.PathLike, track_csv: os.PathLike, ): - """Check that track data written to binary, CSV, and budget files are equal.""" + """Check that binary and CSV track files are equal.""" # get dtype from ascii header file dt = get_track_dtype(track_hdr) @@ -135,8 +136,7 @@ def check_track_data( data_bin.shape == data_csv.shape ), f"Binary and CSV track data shapes do not match: {data_bin.shape} != {data_csv.shape}" - # check particle tracks written to all output files are equal - # check each column separately to avoid: + # check each column separately to avoid # TypeError: The DType could not be promoted by for k in data_bin.dtype.names: if k == "name": @@ -184,7 +184,8 @@ def get_model_name(name, mdl): def get_track_dtype(path: os.PathLike): - """Get the dtype of the track data recarray from the ascii header file.""" + """Read a numpy dtype describing particle track + data format from the ascii track header file.""" hdr_lns = open(path).readlines() hdr_lns_spl = [[ll.strip() for ll in l.split(",")] for l in hdr_lns] @@ -192,6 +193,11 @@ def get_track_dtype(path: os.PathLike): def get_ireason_code(output_event): + """ + Map output event to PRT ireason code specifing + the reason a particle track datum was recorded. + """ + return ( 0 if output_event == "RELEASE" @@ -208,6 +214,10 @@ def get_ireason_code(output_event): def get_partdata(grid, rpts): + """ + Make a flopy.modpath.ParticleData from the given grid and release points. + """ + if grid.grid_type == "structured": return flopy.modpath.ParticleData( partlocs=[grid.get_lrc(p[0])[0] for p in rpts], @@ -230,26 +240,55 @@ def get_partdata(grid, rpts): ) -def get_perioddata(name, periods=1, fraction=None) -> Optional[dict]: - if "reft" in name: - return None - opt = [ - "FIRST" - if "frst" in name - else "ALL" - if "all" in name - else ("STEPS", 1) - if "stps" in name - else None - ] - if opt[0] is None: - raise ValueError(f"Invalid period option: {name}") - if fraction is not None: - opt.append(("FRACTION", fraction)) - return {i: opt for i in range(periods)} - - def has_default_boundnames(data): name = [int(n.partition("0")[2]) for n in data["name"].to_numpy()] irpt = data["irpt"].to_numpy() return np.array_equal(name, irpt) + + +def plot_nodes_and_vertices( + gwf, mg, ibd, ncpl, ax, xmin=None, xmax=None, ymin=None, ymax=None +): + """ + Plot cell nodes and vertices (and IDs) on a zoomed inset + """ + + ax.set_aspect("equal") + xlim = False + ylim = False + if xmin is not None and xmax is not None: + ax.set_xlim([xmin, xmax]) + xlim = True + if ymin is not None and ymax is not None: + ax.set_ylim([ymin, ymax]) + ylim = True + + # create map view plot + pmv = flopy.plot.PlotMapView(gwf, ax=ax) + v = pmv.plot_grid(lw=0.5, edgecolor="black") + t = ax.set_title("Node and vertex indices (one-based)\n", fontsize=14) + ax.set_xlim([xmin, xmax]) + ax.set_ylim([ymin, ymax]) + + # plot vertices + verts = mg.verts + ax.plot(verts[:, 0], verts[:, 1], "bo") + for i in range(ncpl): + x, y = verts[i, 0], verts[i, 1] + ax.annotate(str(i + 1), verts[i, :], color="b") + + # plot nodes + xc, yc = mg.get_xcellcenters_for_layer(0), mg.get_ycellcenters_for_layer(0) + for i in range(ncpl): + x, y = xc[i], yc[i] + ax.plot(x, y, "o", color="grey") + ax.annotate(str(i + 1), (x, y), color="grey") + + # create legend + ax.legend( + handles=[ + mpl.patches.Patch(color="blue", label="vertex"), + mpl.patches.Patch(color="grey", label="node"), + ], + loc="upper left", + ) diff --git a/autotest/prt/test_prt_exg01.py b/autotest/prt/test_prt_exg01.py index ed80291bd93..cd1555c1385 100644 --- a/autotest/prt/test_prt_exg01.py +++ b/autotest/prt/test_prt_exg01.py @@ -25,32 +25,14 @@ from flopy.plot.plotutil import to_mp7_pathlines from flopy.utils import PathlineFile from flopy.utils.binaryfile import HeadFile -from prt_test_utils import (all_equal, check_budget_data, check_track_data, - get_gwf_sim, has_default_boundnames) +from prt_test_utils import check_budget_data, check_track_data, get_gwf_sim from framework import TestFramework from simulation import TestSimulation -# simulation/model names simname = "prtexg01" - -# test cases ex = [simname, f"{simname}bnms"] -# release points -releasepts = [ - # index, k, i, j, x, y, z - # (0-based indexing converted to 1-based for mf6 by flopy) - [i, 0, 0, 0, float(f"0.{i + 1}"), float(f"9.{i + 1}"), 0.5] - for i in range(9) -] -releasepts_mp7 = [ - # node number, localx, localy, localz - # (0-based indexing converted to 1-based for mp7 by flopy) - (0, float(f"0.{i + 1}"), float(f"0.{i + 1}"), 0.5) - for i in range(9) -] - # model names def get_model_name(idx, mdl): @@ -80,9 +62,9 @@ def build_sim(idx, ws, mf6): # create prp package rpts = ( - [r + [str(r[0] + 1)] for r in releasepts] + [r + [str(r[0] + 1)] for r in ctx.releasepts_prt] if "bnms" in name - else releasepts + else ctx.releasepts_prt ) flopy.mf6.ModflowPrtprp( prt, @@ -137,10 +119,10 @@ def build_sim(idx, ws, mf6): def build_mp7_sim(ctx, idx, ws, mp7, gwf): partdata = flopy.modpath.ParticleData( - partlocs=[p[0] for p in releasepts_mp7], - localx=[p[1] for p in releasepts_mp7], - localy=[p[2] for p in releasepts_mp7], - localz=[p[3] for p in releasepts_mp7], + partlocs=[p[0] for p in ctx.releasepts_mp7], + localx=[p[1] for p in ctx.releasepts_mp7], + localy=[p[2] for p in ctx.releasepts_mp7], + localz=[p[3] for p in ctx.releasepts_mp7], timeoffset=0, drape=0, ) diff --git a/autotest/prt/test_prt_fmi01.py b/autotest/prt/test_prt_fmi01.py index 461c89921d9..e5f1197fb25 100644 --- a/autotest/prt/test_prt_fmi01.py +++ b/autotest/prt/test_prt_fmi01.py @@ -36,14 +36,17 @@ from flopy.plot.plotutil import to_mp7_pathlines from flopy.utils import PathlineFile from flopy.utils.binaryfile import HeadFile -from prt_test_utils import (all_equal, check_budget_data, check_track_data, - get_gwf_sim, get_model_name, get_partdata, - has_default_boundnames) +from prt_test_utils import ( + all_equal, + check_budget_data, + check_track_data, + get_gwf_sim, + get_model_name, + get_partdata, + has_default_boundnames, +) -# simulation name simname = "prtfmi01" - -# test cases ex = [simname, f"{simname}saws"] diff --git a/autotest/prt/test_prt_fmi02.py b/autotest/prt/test_prt_fmi02.py index 846a00ba250..17fea05db5d 100644 --- a/autotest/prt/test_prt_fmi02.py +++ b/autotest/prt/test_prt_fmi02.py @@ -36,13 +36,14 @@ from flopy.plot.plotutil import to_mp7_pathlines from flopy.utils import PathlineFile from flopy.utils.binaryfile import HeadFile -from prt_test_utils import (check_budget_data, check_track_data, get_gwf_sim, - get_model_name, get_partdata) +from prt_test_utils import ( + check_budget_data, + check_track_data, + get_gwf_sim, + get_model_name, +) -# simulation/model names simname = "prtfmi02" - -# test cases ex = [ f"{simname}all", f"{simname}rel", @@ -50,18 +51,14 @@ f"{simname}tstp", f"{simname}wksk", ] - -# release points releasepts_prt = { "a": [ # index, k, i, j, x, y, z - # (0-based indexing converted to 1-based for mf6 by flopy) [i, 0, 0, 0, float(f"0.{i + 1}"), float(f"9.{i + 1}"), 0.5] for i in range(4) ], "b": [ # index, k, i, j, x, y, z - # (0-based indexing converted to 1-based for mf6 by flopy) [i, 0, 0, 0, float(f"0.{i + 5}"), float(f"9.{i + 5}"), 0.5] for i in range(5) ], @@ -69,13 +66,11 @@ releasepts_mp7 = { "a": [ # node number, localx, localy, localz - # (0-based indexing converted to 1-based for mf6 by flopy) (0, float(f"0.{i + 1}"), float(f"0.{i + 1}"), 0.5) for i in range(4) ], "b": [ # node number, localx, localy, localz - # (0-based indexing converted to 1-based for mf6 by flopy) (0, float(f"0.{i + 5}"), float(f"0.{i + 5}"), 0.5) for i in range(5) ], diff --git a/autotest/prt/test_prt_fmi03.py b/autotest/prt/test_prt_fmi03.py index 90b197485af..c9884fb308e 100644 --- a/autotest/prt/test_prt_fmi03.py +++ b/autotest/prt/test_prt_fmi03.py @@ -37,16 +37,15 @@ from flopy.utils import PathlineFile from flopy.utils.binaryfile import HeadFile from matplotlib.collections import LineCollection -from prt_test_utils import (check_budget_data, check_track_data, get_gwf_sim, - get_model_name) +from prt_test_utils import ( + check_budget_data, + check_track_data, + get_gwf_sim, + get_model_name, +) -# simulation name simname = "prtfmi03" - -# test cases ex = [f"{simname}_l1", f"{simname}_l2"] - -# izone stopzone_cells = [(0, 1, 8), (0, 8, 1)] @@ -57,22 +56,6 @@ def create_izone(nlay, nrow, ncol): return izone -# release points -# todo: define for mp7 first, then use flopy utils to convert to global coords for mf6 prt -releasepts_mp7 = [ - # node number, localx, localy, localz - # (0-based indexing converted to 1-based for mp7 by flopy) - (0, float(f"0.{i + 1}"), float(f"0.{i + 1}"), 0.5) - for i in range(9) -] -releasepts_prt = [ - # particle index, k, i, j, x, y, z - # (0-based indexing converted to 1-based for mf6 by flopy) - (i, 0, 0, 0, float(f"0.{i + 1}"), float(f"9.{i + 1}"), 0.5) - for i in range(9) -] - - def build_prt_sim(ctx, name, ws, mf6): # create simulation sim = flopy.mf6.MFSimulation( @@ -122,8 +105,8 @@ def build_prt_sim(ctx, name, ws, mf6): prt, pname="prp1", filename=f"{prtname}_1.prp", - nreleasepts=len(releasepts_prt), - packagedata=releasepts_prt, + nreleasepts=len(ctx.releasepts_prt), + packagedata=ctx.releasepts_prt, perioddata={0: ["FIRST"]}, istopzone=1, ) @@ -163,10 +146,10 @@ def build_prt_sim(ctx, name, ws, mf6): def build_mp7_sim(ctx, name, ws, mp7, gwf): partdata = flopy.modpath.ParticleData( - partlocs=[p[0] for p in releasepts_mp7], - localx=[p[1] for p in releasepts_mp7], - localy=[p[2] for p in releasepts_mp7], - localz=[p[3] for p in releasepts_mp7], + partlocs=[p[0] for p in ctx.releasepts_mp7], + localx=[p[1] for p in ctx.releasepts_mp7], + localy=[p[2] for p in ctx.releasepts_mp7], + localz=[p[3] for p in ctx.releasepts_mp7], timeoffset=0, drape=0, ) diff --git a/autotest/prt/test_prt_fmi04.py b/autotest/prt/test_prt_fmi04.py index 961e36ff55d..9e1ae1d2173 100644 --- a/autotest/prt/test_prt_fmi04.py +++ b/autotest/prt/test_prt_fmi04.py @@ -35,29 +35,17 @@ from flopy.plot.plotutil import to_mp7_pathlines from flopy.utils import PathlineFile from flopy.utils.binaryfile import HeadFile -from prt_test_utils import (check_budget_data, check_track_data, get_gwf_sim, - get_ireason_code, get_model_name) +from prt_test_utils import ( + check_budget_data, + check_track_data, + get_gwf_sim, + get_ireason_code, + get_model_name, +) -# simulation name simname = "prtfmi04" - -# test cases ex = [simname, f"{simname}saws"] -# release points -releasepts_mp7 = [ - # node number, localx, localy, localz - # (0-based indexing converted to 1-based for mp7 by flopy) - (0, float(f"0.{i + 1}"), float(f"0.{i + 1}"), 0.5) - for i in range(9) -] -releasepts_prt = [ - # particle index, k, i, j, x, y, z - # (0-based indexing converted to 1-based for mf6 by flopy) - (i, 0, 0, 0, float(f"0.{i + 1}"), float(f"9.{i + 1}"), 0.5) - for i in range(9) -] - def build_prt_sim(ctx, name, ws, mf6): # output file names @@ -106,8 +94,8 @@ def build_prt_sim(ctx, name, ws, mf6): prt, pname="prp1", filename=f"{prtname}_1.prp", - nreleasepts=len(releasepts_prt), - packagedata=releasepts_prt, + nreleasepts=len(ctx.releasepts_prt), + packagedata=ctx.releasepts_prt, perioddata={0: ["FIRST"]}, stop_at_weak_sink="saws" in name, ) @@ -145,10 +133,10 @@ def build_mp7_sim(ctx, name, ws, mp7, gwf): mp7_pathline_file = f"{mp7name}.mppth" partdata = flopy.modpath.ParticleData( - partlocs=[p[0] for p in releasepts_mp7], - localx=[p[1] for p in releasepts_mp7], - localy=[p[2] for p in releasepts_mp7], - localz=[p[3] for p in releasepts_mp7], + partlocs=[p[0] for p in ctx.releasepts_mp7], + localx=[p[1] for p in ctx.releasepts_mp7], + localy=[p[2] for p in ctx.releasepts_mp7], + localz=[p[3] for p in ctx.releasepts_mp7], timeoffset=0, drape=0, ) diff --git a/autotest/prt/test_prt_fmi05.py b/autotest/prt/test_prt_fmi05.py index 5f0be45e705..ccc83c3a725 100644 --- a/autotest/prt/test_prt_fmi05.py +++ b/autotest/prt/test_prt_fmi05.py @@ -37,14 +37,16 @@ from flopy.plot.plotutil import to_mp7_pathlines from flopy.utils import PathlineFile from flopy.utils.binaryfile import HeadFile -from prt_test_utils import (all_equal, check_budget_data, check_track_data, - get_gwf_sim, get_model_name, get_perioddata, - has_default_boundnames) +from prt_test_utils import ( + all_equal, + check_budget_data, + check_track_data, + get_gwf_sim, + get_model_name, + get_partdata, +) -# simulation name simname = "prtfmi05" - -# test cases ex = [ # options block options f"{simname}reft", # REFERENCETIME 0.5 @@ -54,34 +56,24 @@ f"{simname}stps", # STEPS 1 FRACTION 0.5 ] -# release points in mp7 format (using local coordinates) -releasepts_mp7 = [ - # node number, localx, localy, localz - # (0-based indexing converted to 1-based for mp7 by flopy) - (0, float(f"0.{i + 1}"), float(f"0.{i + 1}"), 0.5) - for i in range(9) -] - -# expected release points in PRT format; below we will use flopy -# to convert from mp7 to prt format and make sure they are equal -releasepts_prt = [ - # particle index, k, i, j, x, y, z - # (0-based indexing converted to 1-based for mf6 by flopy) - (i, 0, 0, 0, float(f"0.{i + 1}"), float(f"9.{i + 1}"), 0.5) - for i in range(9) -] - -def get_partdata(grid): - return flopy.modpath.ParticleData( - partlocs=[grid.get_lrc(p[0])[0] for p in releasepts_mp7], - structured=True, - localx=[p[1] for p in releasepts_mp7], - localy=[p[2] for p in releasepts_mp7], - localz=[p[3] for p in releasepts_mp7], - timeoffset=0, - drape=0, - ) +def get_perioddata(name, periods=1, fraction=None) -> Optional[dict]: + if "reft" in name: + return None + opt = [ + "FIRST" + if "frst" in name + else "ALL" + if "all" in name + else ("STEPS", 1) + if "stps" in name + else None + ] + if opt[0] is None: + raise ValueError(f"Invalid period option: {name}") + if fraction is not None: + opt.append(("FRACTION", fraction)) + return {i: opt for i in range(periods)} def build_prt_sim(ctx, name, ws, mf6, fraction=None): @@ -119,11 +111,11 @@ def build_prt_sim(ctx, name, ws, mf6, fraction=None): flopy.mf6.ModflowPrtmip(prt, pname="mip", porosity=ctx.porosity) # convert mp7 particledata to prt release points - partdata = get_partdata(prt.modelgrid) + partdata = get_partdata(prt.modelgrid, ctx.releasepts_mp7) releasepts = list(partdata.to_prp(prt.modelgrid)) # check release points match expectation - assert np.allclose(releasepts_prt, releasepts) + assert np.allclose(ctx.releasepts_prt, releasepts) # create prp package prp_track_file = f"{prtname}.prp.trk" @@ -178,7 +170,7 @@ def build_prt_sim(ctx, name, ws, mf6, fraction=None): def build_mp7_sim(ctx, name, ws, mp7, gwf): # convert mp7 particledata to prt release points - partdata = get_partdata(gwf.modelgrid) + partdata = get_partdata(gwf.modelgrid, ctx.releasepts_mp7) # create modpath 7 simulation mp7name = get_model_name(name, "mp7") diff --git a/autotest/prt/test_prt_fmi06.py b/autotest/prt/test_prt_fmi06.py index 30cab70ec2a..1f26ee6d868 100644 --- a/autotest/prt/test_prt_fmi06.py +++ b/autotest/prt/test_prt_fmi06.py @@ -4,7 +4,6 @@ from pathlib import Path from pprint import pformat -from typing import Optional import flopy import matplotlib.cm as cm @@ -16,16 +15,17 @@ from flopy.utils import PathlineFile from flopy.utils.binaryfile import HeadFile from flopy.utils.gridutil import get_disv_kwargs -from prt_test_utils import (all_equal, check_budget_data, check_track_data, - get_partdata, has_default_boundnames) +from prt_test_utils import ( + all_equal, + check_budget_data, + check_track_data, + get_partdata, + has_default_boundnames, + plot_nodes_and_vertices, +) -# simulation name simname = "prtfmi06" - -# test cases -ex = [ - f"{simname}a", -] +ex = [f"{simname}", f"{simname}bprp"] # model info nlay = 1 @@ -57,10 +57,9 @@ botm, ) -# release points in mp7 format (using local coordinates) +# release points in mp7 format releasepts_mp7 = [ # node number, localx, localy, localz - # (0-based indexing converted to 1-based for mp7 by flopy) (i * 10, 0.5, 0.5, 0.5) for i in range(10) ] @@ -156,8 +155,9 @@ def build_gwf_sim(idx, dir, mf6): def build_prt_sim(idx, ws, mf6): # create simulation + name = ex[idx] sim = flopy.mf6.MFSimulation( - sim_name=ex[idx], + sim_name=name, exe_name=mf6, version="mf6", sim_ws=ws, @@ -181,6 +181,9 @@ def build_prt_sim(idx, ws, mf6): # convert mp7 particledata to prt release points partdata = get_partdata(prt.modelgrid, releasepts_mp7) releasepts = list(partdata.to_prp(prt.modelgrid)) + if "bprp" in name: + # wrong cell index, point is in cell (0, 0) + releasepts[0] = (0, (0, 1), 0.5, 9.5, 22.5) # create prp package prp_track_file = f"{prtname}.prp.trk" @@ -194,7 +197,7 @@ def build_prt_sim(idx, ws, mf6): perioddata={0: ["FIRST"]}, track_filerecord=[prp_track_file], trackcsv_filerecord=[prp_track_csv_file], - stop_at_weak_sink="saws" in prtname, + stop_at_weak_sink=False, boundnames=True, ) @@ -281,10 +284,19 @@ def test_mf6model(idx, name, function_tmpdir, targets): gwfsim = build_gwf_sim(idx, ws, targets.mf6) prtsim = build_prt_sim(idx, ws, targets.mf6) - # run mf6 models - for sim in [gwfsim, prtsim]: - sim.write_simulation() - success, buff = sim.run_simulation(report=True) + # run gwf model + gwfsim.write_simulation() + success, buff = gwfsim.run_simulation(report=True) + assert success, pformat(buff) + + # run prt model + prtsim.write_simulation() + success, buff = prtsim.run_simulation(report=True) + if "bprp" in name: + assert not success, pformat(buff) + assert any("Error: release point" in l for l in buff) + return + else: assert success, pformat(buff) # extract mf6 models @@ -369,6 +381,10 @@ def test_mf6model(idx, name, function_tmpdir, targets): pmv.plot_grid() pmv.plot_array(hds[0], alpha=0.2) pmv.plot_vector(qx, qy, normalize=True, color="white") + # set zoom area + # xmin, xmax = 2050, 4800 + # ymin, ymax = 5200, 7550 + plot_nodes_and_vertices(gwf, mg, None, mg.ncpl, ax[0]) mf6_plines = mf6_pls.groupby(["iprp", "irpt", "trelease"]) for ipl, ((iprp, irpt, trelease), pl) in enumerate(mf6_plines): pl.plot( diff --git a/autotest/prt/test_prt_fmi07.py b/autotest/prt/test_prt_fmi07.py new file mode 100644 index 00000000000..31c02e61010 --- /dev/null +++ b/autotest/prt/test_prt_fmi07.py @@ -0,0 +1,135 @@ +""" +Tests ability to run a GWF model then a PRT model +in separate simulations via flow model interface, +with release points improperly mapped to cell IDs +(expect failures). + +The grid is a 10x10 square with a single layer, +the same flow system shown on the FloPy readme. + +Particles are released from the top left cell. +""" + + +from pprint import pformat + +import flopy +import pytest +from prt_test_utils import ( + get_gwf_sim, + get_model_name, + get_partdata, +) + +simname = "prtfmi07" +ex = [simname] + + +def build_prt_sim(ctx, name, ws, mf6): + # create simulation + sim = flopy.mf6.MFSimulation( + sim_name=name, + exe_name=mf6, + version="mf6", + sim_ws=ws, + ) + + # create tdis package + flopy.mf6.modflow.mftdis.ModflowTdis( + sim, + pname="tdis", + time_units="DAYS", + nper=ctx.nper, + perioddata=[(ctx.perlen, ctx.nstp, ctx.tsmult)], + ) + + # create prt model + prtname = get_model_name(name, "prt") + prt = flopy.mf6.ModflowPrt(sim, modelname=prtname) + + # create prt discretization + flopy.mf6.modflow.mfgwfdis.ModflowGwfdis( + prt, + pname="dis", + nlay=ctx.nlay, + nrow=ctx.nrow, + ncol=ctx.ncol, + ) + + # create mip package + flopy.mf6.ModflowPrtmip(prt, pname="mip", porosity=ctx.porosity) + + # convert mp7 to prt release points and check against expectation + partdata = get_partdata(prt.modelgrid, ctx.releasepts_mp7) + coords = partdata.to_coords(prt.modelgrid) + # bad cell indices! + releasepts = [(i, 0, 1, 1, c[0], c[1], c[2]) for i, c in enumerate(coords)] + + # create prp package + prp_track_file = f"{prtname}.prp.trk" + prp_track_csv_file = f"{prtname}.prp.trk.csv" + flopy.mf6.ModflowPrtprp( + prt, + pname="prp1", + filename=f"{prtname}_1.prp", + nreleasepts=len(releasepts), + packagedata=releasepts, + perioddata={0: ["FIRST"]}, + track_filerecord=[prp_track_file], + trackcsv_filerecord=[prp_track_csv_file], + stop_at_weak_sink="saws" in prtname, + boundnames=True, + ) + + # create output control package + prt_track_file = f"{prtname}.trk" + prt_track_csv_file = f"{prtname}.trk.csv" + flopy.mf6.ModflowPrtoc( + prt, + pname="oc", + track_filerecord=[prt_track_file], + trackcsv_filerecord=[prt_track_csv_file], + ) + + # create the flow model interface + gwfname = get_model_name(name, "gwf") + gwf_budget_file = f"{gwfname}.bud" + gwf_head_file = f"{gwfname}.hds" + flopy.mf6.ModflowPrtfmi( + prt, + packagedata=[ + ("GWFHEAD", gwf_head_file), + ("GWFBUDGET", gwf_budget_file), + ], + ) + + # add explicit model solution + ems = flopy.mf6.ModflowEms( + sim, + pname="ems", + filename=f"{prtname}.ems", + ) + sim.register_solution_package(ems, [prt.name]) + + return sim + + +@pytest.mark.parametrize("name", ex) +def test_mf6model(name, function_tmpdir, targets): + # workspace + ws = function_tmpdir + + # build mf6 models + gwfsim, ctx = get_gwf_sim(name, ws, targets.mf6) + prtsim = build_prt_sim(ctx, name, ws, targets.mf6) + + # run gwf models + gwfsim.write_simulation() + success, buff = gwfsim.run_simulation(report=True) + assert success, pformat(buff) + + # run prt model (expect failure) + prtsim.write_simulation() + success, buff = prtsim.run_simulation(report=True) + assert not success, pformat(buff) + assert any("Error: release point" in l for l in buff) diff --git a/distribution/makedefaults b/distribution/makedefaults new file mode 100644 index 00000000000..ce27db970e6 --- /dev/null +++ b/distribution/makedefaults @@ -0,0 +1,135 @@ +# makedefaults created by pymake (version 1.2.8) for the 'mf6' executable. + +# determine OS +ifeq ($(OS), Windows_NT) + detected_OS = Windows + OS_macro = -D_WIN32 +else + detected_OS = $(shell sh -c 'uname 2>/dev/null || echo Unknown') + ifeq ($(detected_OS), Darwin) + OS_macro = -D__APPLE__ + else + OS_macro = -D__LINUX__ + endif +endif + +# Define the directories for the object and module files +# and the executable and its path. +BINDIR = ../bin +OBJDIR = obj_temp +MODDIR = mod_temp +INCSWITCH = -I $(OBJDIR) +MODSWITCH = -J $(MODDIR) + +# define os dependent program name +ifeq ($(detected_OS), Windows) + PROGRAM = $(BINDIR)/mf6.exe +else ifeq ($(detected_OS), Darwin) + PROGRAM = $(BINDIR)/mf6 +else + PROGRAM = $(BINDIR)/mf6 +endif + +# use GNU compilers if defined compilers do not exist +ifeq ($(detected_OS), Windows) + WHICH = where +else + WHICH = which +endif +ifeq (, $(shell $(WHICH) $(FC))) + FC = gfortran +endif + +# set fortran compiler to gfortran if it is f77 +ifeq ($(FC), f77) + FC = gfortran + # set c compiler to gcc if not passed on the command line + ifneq ($(origin CC), "command line") + ifneq ($(CC), gcc) + CC = gcc + endif + endif +endif + +# set the optimization level (OPTLEVEL) if not defined +OPTLEVEL ?= -O2 + +# set the fortran flags +ifeq ($(detected_OS), Windows) + ifeq ($(FC), gfortran) + FFLAGS ?= -static -fbacktrace -ffpe-summary=overflow -ffpe-trap=overflow,zero,invalid $(OS_macro) -fall-intrinsics -Wtabs -Wline-truncation -Wunused-label -Wunused-variable -pedantic -std=f2008 -Wcharacter-truncation -cpp + endif +else + ifeq ($(FC), gfortran) + FFLAGS ?= -fbacktrace -ffpe-summary=overflow -ffpe-trap=overflow,zero,invalid $(OS_macro) -fall-intrinsics -Wtabs -Wline-truncation -Wunused-label -Wunused-variable -pedantic -std=f2008 -Wcharacter-truncation -cpp + endif + ifeq ($(FC), $(filter $(FC), ifort mpiifort)) + FFLAGS ?= -no-heap-arrays -fpe0 -traceback -fpp + MODSWITCH = -module $(MODDIR) + endif + ifeq ($(FC), $(filter $(FC), ftn)) + FFLAGS ?= -h noheap_allocate + endif +endif + +# set the ldflgs +ifeq ($(detected_OS), Windows) + ifeq ($(FC), $(filter $(FC), gfortran)) + LDFLAGS ?= -static -static-libgfortran -static-libgcc -static-libstdc++ -lm + endif +else + ifeq ($(FC), gfortran) + LDFLAGS ?= -lc + endif + ifeq ($(FC), $(filter $(FC), ifort mpiifort)) + LDFLAGS ?= -lc + endif + ifeq ($(FC), $(filter $(FC), ftn)) + LDFLAGS ?= -lc + endif +endif + +# check for Windows error condition +ifeq ($(detected_OS), Windows) + ifeq ($(FC), $(filter $(FC), ifort mpiifort)) + WINDOWSERROR = $(FC) + endif +endif + +# Define task functions +# Create the bin directory and compile and link the program +all: windowscheck makedirs | $(PROGRAM) + +# test for windows error +windowscheck: +ifdef WINDOWSERROR + $(error cannot use makefile on windows with $(WINDOWSERROR)) +endif + +# Make the bin directory for the executable +makedirs: + mkdir -p $(BINDIR) + mkdir -p $(MODDIR) + +# Write selected compiler settings +.PHONY: settings +settings: + @echo "Optimization level: $(OPTLEVEL)" + @echo "Fortran compiler: $(FC)" + @echo "Fortran flags: $(FFLAGS)" + @echo "Linker: $(FC)" + @echo "SYSLIBS: $(LDFLAGS)" + +# Clean the object and module files and the executable +.PHONY: clean +clean: + -rm -rf $(OBJDIR) + -rm -rf $(MODDIR) + -rm -rf $(PROGRAM) + +# Clean the object and module files +.PHONY: cleanobj +cleanobj: + -rm -rf $(OBJDIR) + -rm -rf $(MODDIR) + diff --git a/distribution/makefile b/distribution/makefile new file mode 100644 index 00000000000..a1ac2689cc0 --- /dev/null +++ b/distribution/makefile @@ -0,0 +1,385 @@ +# makefile created by pymake (version 1.2.8) for the 'mf6' executable. + + +include ./makedefaults + +# Define the source file directories +SOURCEDIR1=../src +SOURCEDIR2=../src/Exchange +SOURCEDIR3=../src/Distributed +SOURCEDIR4=../src/Solution +SOURCEDIR5=../src/Solution/LinearMethods +SOURCEDIR6=../src/Solution/ParticleTracker +SOURCEDIR7=../src/Solution/PETSc +SOURCEDIR8=../src/Timing +SOURCEDIR9=../src/Utilities +SOURCEDIR10=../src/Utilities/Idm +SOURCEDIR11=../src/Utilities/Idm/selector +SOURCEDIR12=../src/Utilities/Idm/mf6blockfile +SOURCEDIR13=../src/Utilities/TimeSeries +SOURCEDIR14=../src/Utilities/Memory +SOURCEDIR15=../src/Utilities/OutputControl +SOURCEDIR16=../src/Utilities/ArrayRead +SOURCEDIR17=../src/Utilities/Libraries +SOURCEDIR18=../src/Utilities/Libraries/rcm +SOURCEDIR19=../src/Utilities/Libraries/blas +SOURCEDIR20=../src/Utilities/Libraries/sparskit2 +SOURCEDIR21=../src/Utilities/Libraries/daglib +SOURCEDIR22=../src/Utilities/Libraries/sparsekit +SOURCEDIR23=../src/Utilities/Vector +SOURCEDIR24=../src/Utilities/Matrix +SOURCEDIR25=../src/Utilities/Observation +SOURCEDIR26=../src/Model +SOURCEDIR27=../src/Model/Connection +SOURCEDIR28=../src/Model/ParticleTracking +SOURCEDIR29=../src/Model/GroundWaterTransport +SOURCEDIR30=../src/Model/ModelUtilities +SOURCEDIR31=../src/Model/GroundWaterFlow +SOURCEDIR32=../src/Model/TransportModel +SOURCEDIR33=../src/Model/Geometry + +VPATH = \ +${SOURCEDIR1} \ +${SOURCEDIR2} \ +${SOURCEDIR3} \ +${SOURCEDIR4} \ +${SOURCEDIR5} \ +${SOURCEDIR6} \ +${SOURCEDIR7} \ +${SOURCEDIR8} \ +${SOURCEDIR9} \ +${SOURCEDIR10} \ +${SOURCEDIR11} \ +${SOURCEDIR12} \ +${SOURCEDIR13} \ +${SOURCEDIR14} \ +${SOURCEDIR15} \ +${SOURCEDIR16} \ +${SOURCEDIR17} \ +${SOURCEDIR18} \ +${SOURCEDIR19} \ +${SOURCEDIR20} \ +${SOURCEDIR21} \ +${SOURCEDIR22} \ +${SOURCEDIR23} \ +${SOURCEDIR24} \ +${SOURCEDIR25} \ +${SOURCEDIR26} \ +${SOURCEDIR27} \ +${SOURCEDIR28} \ +${SOURCEDIR29} \ +${SOURCEDIR30} \ +${SOURCEDIR31} \ +${SOURCEDIR32} \ +${SOURCEDIR33} + +.SUFFIXES: .f90 .F90 .o + +OBJECTS = \ +$(OBJDIR)/kind.o \ +$(OBJDIR)/Constants.o \ +$(OBJDIR)/SimVariables.o \ +$(OBJDIR)/genericutils.o \ +$(OBJDIR)/defmacro.o \ +$(OBJDIR)/compilerversion.o \ +$(OBJDIR)/ArrayHandlers.o \ +$(OBJDIR)/version.o \ +$(OBJDIR)/Message.o \ +$(OBJDIR)/Sim.o \ +$(OBJDIR)/OpenSpec.o \ +$(OBJDIR)/InputOutput.o \ +$(OBJDIR)/TableTerm.o \ +$(OBJDIR)/Table.o \ +$(OBJDIR)/MemoryHelper.o \ +$(OBJDIR)/CharString.o \ +$(OBJDIR)/Memory.o \ +$(OBJDIR)/List.o \ +$(OBJDIR)/LongLineReader.o \ +$(OBJDIR)/DevFeature.o \ +$(OBJDIR)/MemoryList.o \ +$(OBJDIR)/TimeSeriesRecord.o \ +$(OBJDIR)/BlockParser.o \ +$(OBJDIR)/MemoryManager.o \ +$(OBJDIR)/TimeSeries.o \ +$(OBJDIR)/ats.o \ +$(OBJDIR)/TimeSeriesLink.o \ +$(OBJDIR)/TimeSeriesFileList.o \ +$(OBJDIR)/tdis.o \ +$(OBJDIR)/HashTable.o \ +$(OBJDIR)/VectorBase.o \ +$(OBJDIR)/Sparse.o \ +$(OBJDIR)/DisvGeom.o \ +$(OBJDIR)/ArrayReaders.o \ +$(OBJDIR)/TimeSeriesManager.o \ +$(OBJDIR)/SmoothingFunctions.o \ +$(OBJDIR)/MemoryManagerExt.o \ +$(OBJDIR)/MatrixBase.o \ +$(OBJDIR)/ListReader.o \ +$(OBJDIR)/Connections.o \ +$(OBJDIR)/InputDefinition.o \ +$(OBJDIR)/TimeArray.o \ +$(OBJDIR)/ObsOutput.o \ +$(OBJDIR)/DiscretizationBase.o \ +$(OBJDIR)/simnamidm.o \ +$(OBJDIR)/prt1idm.o \ +$(OBJDIR)/prt1mip1idm.o \ +$(OBJDIR)/prt1disv1idm.o \ +$(OBJDIR)/prt1dis1idm.o \ +$(OBJDIR)/gwt1idm.o \ +$(OBJDIR)/gwt1ic1idm.o \ +$(OBJDIR)/gwt1dsp1idm.o \ +$(OBJDIR)/gwt1disv1idm.o \ +$(OBJDIR)/gwt1disu1idm.o \ +$(OBJDIR)/gwt1dis1idm.o \ +$(OBJDIR)/gwt1cnc1idm.o \ +$(OBJDIR)/gwf3wel8idm.o \ +$(OBJDIR)/gwf3riv8idm.o \ +$(OBJDIR)/gwf3rch8idm.o \ +$(OBJDIR)/gwf3rcha8idm.o \ +$(OBJDIR)/gwf3npf8idm.o \ +$(OBJDIR)/gwf3idm.o \ +$(OBJDIR)/gwf3ic8idm.o \ +$(OBJDIR)/gwf3ghb8idm.o \ +$(OBJDIR)/gwf3evt8idm.o \ +$(OBJDIR)/gwf3evta8idm.o \ +$(OBJDIR)/gwf3drn8idm.o \ +$(OBJDIR)/gwf3disv8idm.o \ +$(OBJDIR)/gwf3disu8idm.o \ +$(OBJDIR)/gwf3dis8idm.o \ +$(OBJDIR)/gwf3chd8idm.o \ +$(OBJDIR)/TimeArraySeries.o \ +$(OBJDIR)/ObsOutputList.o \ +$(OBJDIR)/Observe.o \ +$(OBJDIR)/CellDefn.o \ +$(OBJDIR)/IdmSimDfnSelector.o \ +$(OBJDIR)/IdmPrtDfnSelector.o \ +$(OBJDIR)/IdmGwtDfnSelector.o \ +$(OBJDIR)/IdmGwfDfnSelector.o \ +$(OBJDIR)/TimeArraySeriesLink.o \ +$(OBJDIR)/ObsUtility.o \ +$(OBJDIR)/ObsContainer.o \ +$(OBJDIR)/UtilMisc.o \ +$(OBJDIR)/GlobalData.o \ +$(OBJDIR)/BudgetFileReader.o \ +$(OBJDIR)/IdmDfnSelector.o \ +$(OBJDIR)/TimeArraySeriesManager.o \ +$(OBJDIR)/PackageMover.o \ +$(OBJDIR)/Obs3.o \ +$(OBJDIR)/NumericalPackage.o \ +$(OBJDIR)/Budget.o \ +$(OBJDIR)/Particle.o \ +$(OBJDIR)/BudgetTerm.o \ +$(OBJDIR)/sort.o \ +$(OBJDIR)/SfrCrossSectionUtils.o \ +$(OBJDIR)/SourceCommon.o \ +$(OBJDIR)/BoundaryPackage.o \ +$(OBJDIR)/TernaryUtil.o \ +$(OBJDIR)/Ternary.o \ +$(OBJDIR)/TrackData.o \ +$(OBJDIR)/VirtualBase.o \ +$(OBJDIR)/STLVecInt.o \ +$(OBJDIR)/BaseModel.o \ +$(OBJDIR)/PackageBudget.o \ +$(OBJDIR)/HeadFileReader.o \ +$(OBJDIR)/BudgetObject.o \ +$(OBJDIR)/PrintSaveManager.o \ +$(OBJDIR)/SfrCrossSectionManager.o \ +$(OBJDIR)/dag_module.o \ +$(OBJDIR)/BoundaryPackageExt.o \ +$(OBJDIR)/TernarySolveTrack.o \ +$(OBJDIR)/SubcellTri.o \ +$(OBJDIR)/Method.o \ +$(OBJDIR)/SubcellRect.o \ +$(OBJDIR)/VirtualDataLists.o \ +$(OBJDIR)/VirtualDataContainer.o \ +$(OBJDIR)/SimStages.o \ +$(OBJDIR)/NumericalModel.o \ +$(OBJDIR)/FlowModelInterface.o \ +$(OBJDIR)/OutputControlData.o \ +$(OBJDIR)/gwf3ic8.o \ +$(OBJDIR)/Xt3dAlgorithm.o \ +$(OBJDIR)/gwf3tvbase8.o \ +$(OBJDIR)/gwf3sfr8.o \ +$(OBJDIR)/gwf3riv8.o \ +$(OBJDIR)/gwf3maw8.o \ +$(OBJDIR)/mf6lists.o \ +$(OBJDIR)/gwf3lak8.o \ +$(OBJDIR)/GwfVscInputData.o \ +$(OBJDIR)/gwf3ghb8.o \ +$(OBJDIR)/gwf3drn8.o \ +$(OBJDIR)/IndexMap.o \ +$(OBJDIR)/MethodSubcellTernary.o \ +$(OBJDIR)/MethodSubcellPollock.o \ +$(OBJDIR)/VirtualModel.o \ +$(OBJDIR)/BaseExchange.o \ +$(OBJDIR)/tsp1fmi1.o \ +$(OBJDIR)/GwtSpc.o \ +$(OBJDIR)/OutputControl.o \ +$(OBJDIR)/tsp1ic1.o \ +$(OBJDIR)/TspAdvOptions.o \ +$(OBJDIR)/UzfCellGroup.o \ +$(OBJDIR)/Xt3dInterface.o \ +$(OBJDIR)/gwf3tvk8.o \ +$(OBJDIR)/gwf3vsc8.o \ +$(OBJDIR)/GwfNpfOptions.o \ +$(OBJDIR)/InterfaceMap.o \ +$(OBJDIR)/SeqVector.o \ +$(OBJDIR)/ImsLinearSettings.o \ +$(OBJDIR)/ConvergenceSummary.o \ +$(OBJDIR)/ArrayReaderBase.o \ +$(OBJDIR)/MethodSubcellPool.o \ +$(OBJDIR)/CellPoly.o \ +$(OBJDIR)/CellRectQuad.o \ +$(OBJDIR)/CellRect.o \ +$(OBJDIR)/CellWithNbrs.o \ +$(OBJDIR)/NumericalExchange.o \ +$(OBJDIR)/tsp1ssm1.o \ +$(OBJDIR)/tsp1oc1.o \ +$(OBJDIR)/tsp1obs1.o \ +$(OBJDIR)/tsp1mvt1.o \ +$(OBJDIR)/tsp1adv1.o \ +$(OBJDIR)/gwf3disv8.o \ +$(OBJDIR)/gwf3disu8.o \ +$(OBJDIR)/gwf3dis8.o \ +$(OBJDIR)/gwf3uzf8.o \ +$(OBJDIR)/tsp1apt1.o \ +$(OBJDIR)/gwt1mst1.o \ +$(OBJDIR)/GwtDspOptions.o \ +$(OBJDIR)/gwf3npf8.o \ +$(OBJDIR)/gwf3tvs8.o \ +$(OBJDIR)/GwfStorageUtils.o \ +$(OBJDIR)/Mover.o \ +$(OBJDIR)/GwfMvrPeriodData.o \ +$(OBJDIR)/ims8misc.o \ +$(OBJDIR)/GwfBuyInputData.o \ +$(OBJDIR)/VirtualSolution.o \ +$(OBJDIR)/SparseMatrix.o \ +$(OBJDIR)/LinearSolverBase.o \ +$(OBJDIR)/ims8reordering.o \ +$(OBJDIR)/ModflowInput.o \ +$(OBJDIR)/Integer2dReader.o \ +$(OBJDIR)/MethodCellTernary.o \ +$(OBJDIR)/MethodCellPollockQuad.o \ +$(OBJDIR)/MethodCellPollock.o \ +$(OBJDIR)/MethodCellPassToBot.o \ +$(OBJDIR)/VirtualExchange.o \ +$(OBJDIR)/GridSorting.o \ +$(OBJDIR)/DisConnExchange.o \ +$(OBJDIR)/CsrUtils.o \ +$(OBJDIR)/tsp1.o \ +$(OBJDIR)/gwt1uzt1.o \ +$(OBJDIR)/gwt1src1.o \ +$(OBJDIR)/gwt1sft1.o \ +$(OBJDIR)/gwt1mwt1.o \ +$(OBJDIR)/gwt1lkt1.o \ +$(OBJDIR)/gwt1ist1.o \ +$(OBJDIR)/gwt1dsp1.o \ +$(OBJDIR)/gwt1cnc1.o \ +$(OBJDIR)/gwf3api8.o \ +$(OBJDIR)/gwf3wel8.o \ +$(OBJDIR)/gwf3rch8.o \ +$(OBJDIR)/gwf3sto8.o \ +$(OBJDIR)/gwf3oc8.o \ +$(OBJDIR)/gwf3obs8.o \ +$(OBJDIR)/gwf3mvr8.o \ +$(OBJDIR)/gwf3hfb8.o \ +$(OBJDIR)/gwf3csub8.o \ +$(OBJDIR)/gwf3buy8.o \ +$(OBJDIR)/GhostNode.o \ +$(OBJDIR)/gwf3evt8.o \ +$(OBJDIR)/gwf3chd8.o \ +$(OBJDIR)/RouterBase.o \ +$(OBJDIR)/ImsLinearSolver.o \ +$(OBJDIR)/ims8base.o \ +$(OBJDIR)/StructVector.o \ +$(OBJDIR)/IdmLogger.o \ +$(OBJDIR)/DefinitionSelect.o \ +$(OBJDIR)/InputLoadType.o \ +$(OBJDIR)/Integer1dReader.o \ +$(OBJDIR)/Double2dReader.o \ +$(OBJDIR)/Double1dReader.o \ +$(OBJDIR)/ExplicitModel.o \ +$(OBJDIR)/prt1fmi1.o \ +$(OBJDIR)/GeomUtil.o \ +$(OBJDIR)/MethodCellPool.o \ +$(OBJDIR)/CellUtil.o \ +$(OBJDIR)/GridConnection.o \ +$(OBJDIR)/DistributedVariable.o \ +$(OBJDIR)/gwt1.o \ +$(OBJDIR)/gwf3.o \ +$(OBJDIR)/SerialRouter.o \ +$(OBJDIR)/Timer.o \ +$(OBJDIR)/LinearSolverFactory.o \ +$(OBJDIR)/ims8linear.o \ +$(OBJDIR)/BaseSolution.o \ +$(OBJDIR)/StructArray.o \ +$(OBJDIR)/BoundInputContext.o \ +$(OBJDIR)/AsciiInputLoadType.o \ +$(OBJDIR)/LayeredArrayReader.o \ +$(OBJDIR)/TrackingModel.o \ +$(OBJDIR)/prt1prp1.o \ +$(OBJDIR)/prt1oc1.o \ +$(OBJDIR)/prt1obs1.o \ +$(OBJDIR)/prt1mip.o \ +$(OBJDIR)/MethodDisv.o \ +$(OBJDIR)/MethodDis.o \ +$(OBJDIR)/SpatialModelConnection.o \ +$(OBJDIR)/GwtInterfaceModel.o \ +$(OBJDIR)/GwtGwtExchange.o \ +$(OBJDIR)/GwfInterfaceModel.o \ +$(OBJDIR)/GwfGwfExchange.o \ +$(OBJDIR)/RouterFactory.o \ +$(OBJDIR)/NumericalSolution.o \ +$(OBJDIR)/MappedMemory.o \ +$(OBJDIR)/StressListInput.o \ +$(OBJDIR)/StressGridInput.o \ +$(OBJDIR)/LoadMf6File.o \ +$(OBJDIR)/prt1.o \ +$(OBJDIR)/ExplicitSolution.o \ +$(OBJDIR)/GwtGwtConnection.o \ +$(OBJDIR)/GwfGwfConnection.o \ +$(OBJDIR)/VirtualDataManager.o \ +$(OBJDIR)/Mapper.o \ +$(OBJDIR)/IdmMf6File.o \ +$(OBJDIR)/ModelPackageInput.o \ +$(OBJDIR)/VirtualGwtModel.o \ +$(OBJDIR)/VirtualGwtExchange.o \ +$(OBJDIR)/VirtualGwfModel.o \ +$(OBJDIR)/VirtualGwfExchange.o \ +$(OBJDIR)/SolutionGroup.o \ +$(OBJDIR)/SolutionFactory.o \ +$(OBJDIR)/GwfPrtExchange.o \ +$(OBJDIR)/GwfGwtExchange.o \ +$(OBJDIR)/RunControl.o \ +$(OBJDIR)/SourceLoad.o \ +$(OBJDIR)/ModelPackageInputs.o \ +$(OBJDIR)/SimulationCreate.o \ +$(OBJDIR)/RunControlFactory.o \ +$(OBJDIR)/IdmLoad.o \ +$(OBJDIR)/ConnectionBuilder.o \ +$(OBJDIR)/comarg.o \ +$(OBJDIR)/mf6core.o \ +$(OBJDIR)/BaseGeometry.o \ +$(OBJDIR)/mf6.o \ +$(OBJDIR)/StringList.o \ +$(OBJDIR)/MemorySetHandler.o \ +$(OBJDIR)/ilut.o \ +$(OBJDIR)/sparsekit.o \ +$(OBJDIR)/rcm.o \ +$(OBJDIR)/blas1_d.o \ +$(OBJDIR)/Iunit.o \ +$(OBJDIR)/RectangularGeometry.o \ +$(OBJDIR)/CircularGeometry.o + +# Define the objects that make up the program +$(PROGRAM) : $(OBJECTS) + -$(FC) $(OPTLEVEL) -o $@ $(OBJECTS) $(LDFLAGS) + +$(OBJDIR)/%.o : %.f90 + @mkdir -p $(@D) + $(FC) $(OPTLEVEL) $(FFLAGS) -c $< -o $@ $(INCSWITCH) $(MODSWITCH) + +$(OBJDIR)/%.o : %.F90 + @mkdir -p $(@D) + $(FC) $(OPTLEVEL) $(FFLAGS) -c $< -o $@ $(INCSWITCH) $(MODSWITCH) + diff --git a/environment.yml b/environment.yml index 9b735965ae6..dd3ebc9fa04 100644 --- a/environment.yml +++ b/environment.yml @@ -19,12 +19,11 @@ dependencies: - shapely - pip - pip: - - https://github.com/w-bonelli/flopy/archive/prt-utils.zip + - https://github.com/wpbonelli/flopy/archive/prt-utils.zip - git+https://github.com/modflowpy/pymake.git - git+https://github.com/Deltares/xmipy.git - git+https://github.com/MODFLOW-USGS/modflowapi.git - modflow-devtools - - flopy[test,optional] - fortls - pytest - pytest-cases diff --git a/make/makefile b/make/makefile index 55971897078..7364de98f43 100644 --- a/make/makefile +++ b/make/makefile @@ -150,6 +150,7 @@ $(OBJDIR)/gwf3chd8idm.o \ $(OBJDIR)/TimeArraySeries.o \ $(OBJDIR)/ObsOutputList.o \ $(OBJDIR)/Observe.o \ +$(OBJDIR)/CellDefn.o \ $(OBJDIR)/IdmSimDfnSelector.o \ $(OBJDIR)/IdmPrtDfnSelector.o \ $(OBJDIR)/IdmGwtDfnSelector.o \ @@ -157,6 +158,8 @@ $(OBJDIR)/IdmGwfDfnSelector.o \ $(OBJDIR)/TimeArraySeriesLink.o \ $(OBJDIR)/ObsUtility.o \ $(OBJDIR)/ObsContainer.o \ +$(OBJDIR)/UtilMisc.o \ +$(OBJDIR)/GlobalData.o \ $(OBJDIR)/BudgetFileReader.o \ $(OBJDIR)/IdmDfnSelector.o \ $(OBJDIR)/TimeArraySeriesManager.o \ @@ -164,14 +167,18 @@ $(OBJDIR)/PackageMover.o \ $(OBJDIR)/Obs3.o \ $(OBJDIR)/NumericalPackage.o \ $(OBJDIR)/Budget.o \ -$(OBJDIR)/CellDefn.o \ +$(OBJDIR)/Particle.o \ $(OBJDIR)/BudgetTerm.o \ $(OBJDIR)/sort.o \ $(OBJDIR)/SfrCrossSectionUtils.o \ $(OBJDIR)/SourceCommon.o \ $(OBJDIR)/BoundaryPackage.o \ -$(OBJDIR)/UtilMisc.o \ -$(OBJDIR)/GlobalData.o \ +$(OBJDIR)/TernaryUtil.o \ +$(OBJDIR)/Ternary.o \ +$(OBJDIR)/TrackData.o \ +$(OBJDIR)/VirtualBase.o \ +$(OBJDIR)/STLVecInt.o \ +$(OBJDIR)/BaseModel.o \ $(OBJDIR)/PackageBudget.o \ $(OBJDIR)/HeadFileReader.o \ $(OBJDIR)/BudgetObject.o \ @@ -179,7 +186,14 @@ $(OBJDIR)/PrintSaveManager.o \ $(OBJDIR)/SfrCrossSectionManager.o \ $(OBJDIR)/dag_module.o \ $(OBJDIR)/BoundaryPackageExt.o \ -$(OBJDIR)/Particle.o \ +$(OBJDIR)/TernarySolveTrack.o \ +$(OBJDIR)/SubcellTri.o \ +$(OBJDIR)/Method.o \ +$(OBJDIR)/SubcellRect.o \ +$(OBJDIR)/VirtualDataLists.o \ +$(OBJDIR)/VirtualDataContainer.o \ +$(OBJDIR)/SimStages.o \ +$(OBJDIR)/NumericalModel.o \ $(OBJDIR)/FlowModelInterface.o \ $(OBJDIR)/OutputControlData.o \ $(OBJDIR)/gwf3ic8.o \ @@ -193,12 +207,11 @@ $(OBJDIR)/gwf3lak8.o \ $(OBJDIR)/GwfVscInputData.o \ $(OBJDIR)/gwf3ghb8.o \ $(OBJDIR)/gwf3drn8.o \ -$(OBJDIR)/TernaryUtil.o \ -$(OBJDIR)/Ternary.o \ -$(OBJDIR)/TrackData.o \ -$(OBJDIR)/VirtualBase.o \ -$(OBJDIR)/STLVecInt.o \ -$(OBJDIR)/BaseModel.o \ +$(OBJDIR)/IndexMap.o \ +$(OBJDIR)/MethodSubcellTernary.o \ +$(OBJDIR)/MethodSubcellPollock.o \ +$(OBJDIR)/VirtualModel.o \ +$(OBJDIR)/BaseExchange.o \ $(OBJDIR)/tsp1fmi1.o \ $(OBJDIR)/GwtSpc.o \ $(OBJDIR)/OutputControl.o \ @@ -209,15 +222,17 @@ $(OBJDIR)/Xt3dInterface.o \ $(OBJDIR)/gwf3tvk8.o \ $(OBJDIR)/gwf3vsc8.o \ $(OBJDIR)/GwfNpfOptions.o \ -$(OBJDIR)/TernarySolveTrack.o \ -$(OBJDIR)/SubcellTri.o \ -$(OBJDIR)/Method.o \ -$(OBJDIR)/SubcellRect.o \ -$(OBJDIR)/VirtualDataLists.o \ -$(OBJDIR)/VirtualDataContainer.o \ -$(OBJDIR)/SimStages.o \ -$(OBJDIR)/NumericalModel.o \ -$(OBJDIR)/IndexMap.o \ +$(OBJDIR)/InterfaceMap.o \ +$(OBJDIR)/SeqVector.o \ +$(OBJDIR)/ImsLinearSettings.o \ +$(OBJDIR)/ConvergenceSummary.o \ +$(OBJDIR)/ArrayReaderBase.o \ +$(OBJDIR)/MethodSubcellPool.o \ +$(OBJDIR)/CellPoly.o \ +$(OBJDIR)/CellRectQuad.o \ +$(OBJDIR)/CellRect.o \ +$(OBJDIR)/CellWithNbrs.o \ +$(OBJDIR)/NumericalExchange.o \ $(OBJDIR)/tsp1ssm1.o \ $(OBJDIR)/tsp1oc1.o \ $(OBJDIR)/tsp1obs1.o \ @@ -237,15 +252,20 @@ $(OBJDIR)/Mover.o \ $(OBJDIR)/GwfMvrPeriodData.o \ $(OBJDIR)/ims8misc.o \ $(OBJDIR)/GwfBuyInputData.o \ -$(OBJDIR)/MethodSubcellTernary.o \ -$(OBJDIR)/MethodSubcellPollock.o \ -$(OBJDIR)/VirtualModel.o \ -$(OBJDIR)/BaseExchange.o \ -$(OBJDIR)/InterfaceMap.o \ -$(OBJDIR)/SeqVector.o \ -$(OBJDIR)/ImsLinearSettings.o \ -$(OBJDIR)/ConvergenceSummary.o \ -$(OBJDIR)/ArrayReaderBase.o \ +$(OBJDIR)/VirtualSolution.o \ +$(OBJDIR)/SparseMatrix.o \ +$(OBJDIR)/LinearSolverBase.o \ +$(OBJDIR)/ims8reordering.o \ +$(OBJDIR)/ModflowInput.o \ +$(OBJDIR)/Integer2dReader.o \ +$(OBJDIR)/MethodCellTernary.o \ +$(OBJDIR)/MethodCellPollockQuad.o \ +$(OBJDIR)/MethodCellPollock.o \ +$(OBJDIR)/MethodCellPassToBot.o \ +$(OBJDIR)/VirtualExchange.o \ +$(OBJDIR)/GridSorting.o \ +$(OBJDIR)/DisConnExchange.o \ +$(OBJDIR)/CsrUtils.o \ $(OBJDIR)/tsp1.o \ $(OBJDIR)/gwt1uzt1.o \ $(OBJDIR)/gwt1src1.o \ @@ -268,28 +288,6 @@ $(OBJDIR)/gwf3buy8.o \ $(OBJDIR)/GhostNode.o \ $(OBJDIR)/gwf3evt8.o \ $(OBJDIR)/gwf3chd8.o \ -$(OBJDIR)/MethodSubcellPool.o \ -$(OBJDIR)/CellPoly.o \ -$(OBJDIR)/CellRectQuad.o \ -$(OBJDIR)/CellRect.o \ -$(OBJDIR)/CellWithNbrs.o \ -$(OBJDIR)/NumericalExchange.o \ -$(OBJDIR)/VirtualSolution.o \ -$(OBJDIR)/SparseMatrix.o \ -$(OBJDIR)/LinearSolverBase.o \ -$(OBJDIR)/ims8reordering.o \ -$(OBJDIR)/ModflowInput.o \ -$(OBJDIR)/Integer2dReader.o \ -$(OBJDIR)/gwt1.o \ -$(OBJDIR)/gwf3.o \ -$(OBJDIR)/MethodCellTernary.o \ -$(OBJDIR)/MethodCellPollockQuad.o \ -$(OBJDIR)/MethodCellPollock.o \ -$(OBJDIR)/MethodCellPassToBot.o \ -$(OBJDIR)/VirtualExchange.o \ -$(OBJDIR)/GridSorting.o \ -$(OBJDIR)/DisConnExchange.o \ -$(OBJDIR)/CsrUtils.o \ $(OBJDIR)/RouterBase.o \ $(OBJDIR)/ImsLinearSolver.o \ $(OBJDIR)/ims8base.o \ @@ -302,11 +300,13 @@ $(OBJDIR)/Double2dReader.o \ $(OBJDIR)/Double1dReader.o \ $(OBJDIR)/ExplicitModel.o \ $(OBJDIR)/prt1fmi1.o \ -$(OBJDIR)/ModelPackageInput.o \ +$(OBJDIR)/GeomUtil.o \ $(OBJDIR)/MethodCellPool.o \ $(OBJDIR)/CellUtil.o \ $(OBJDIR)/GridConnection.o \ $(OBJDIR)/DistributedVariable.o \ +$(OBJDIR)/gwt1.o \ +$(OBJDIR)/gwf3.o \ $(OBJDIR)/SerialRouter.o \ $(OBJDIR)/Timer.o \ $(OBJDIR)/LinearSolverFactory.o \ @@ -321,7 +321,6 @@ $(OBJDIR)/prt1prp1.o \ $(OBJDIR)/prt1oc1.o \ $(OBJDIR)/prt1obs1.o \ $(OBJDIR)/prt1mip.o \ -$(OBJDIR)/ModelPackageInputs.o \ $(OBJDIR)/MethodDisv.o \ $(OBJDIR)/MethodDis.o \ $(OBJDIR)/SpatialModelConnection.o \ @@ -335,13 +334,14 @@ $(OBJDIR)/MappedMemory.o \ $(OBJDIR)/StressListInput.o \ $(OBJDIR)/StressGridInput.o \ $(OBJDIR)/LoadMf6File.o \ -$(OBJDIR)/ExplicitSolution.o \ $(OBJDIR)/prt1.o \ +$(OBJDIR)/ExplicitSolution.o \ $(OBJDIR)/GwtGwtConnection.o \ $(OBJDIR)/GwfGwfConnection.o \ $(OBJDIR)/VirtualDataManager.o \ $(OBJDIR)/Mapper.o \ $(OBJDIR)/IdmMf6File.o \ +$(OBJDIR)/ModelPackageInput.o \ $(OBJDIR)/VirtualGwtModel.o \ $(OBJDIR)/VirtualGwtExchange.o \ $(OBJDIR)/VirtualGwfModel.o \ @@ -352,6 +352,7 @@ $(OBJDIR)/GwfPrtExchange.o \ $(OBJDIR)/GwfGwtExchange.o \ $(OBJDIR)/RunControl.o \ $(OBJDIR)/SourceLoad.o \ +$(OBJDIR)/ModelPackageInputs.o \ $(OBJDIR)/SimulationCreate.o \ $(OBJDIR)/RunControlFactory.o \ $(OBJDIR)/IdmLoad.o \ @@ -367,7 +368,6 @@ $(OBJDIR)/sparsekit.o \ $(OBJDIR)/rcm.o \ $(OBJDIR)/blas1_d.o \ $(OBJDIR)/Iunit.o \ -$(OBJDIR)/GeomUtil.o \ $(OBJDIR)/RectangularGeometry.o \ $(OBJDIR)/CircularGeometry.o diff --git a/src/Model/ParticleTracking/prt1.f90 b/src/Model/ParticleTracking/prt1.f90 index bf69e823b68..69be12c6649 100644 --- a/src/Model/ParticleTracking/prt1.f90 +++ b/src/Model/ParticleTracking/prt1.f90 @@ -123,6 +123,9 @@ module PrtModule data PRT_MULTIPKG/'PRP6 ', ' ', ' ', ' ', ' ', & ! 5 &45*' '/ ! 50 + ! -- size of supported model package arrays + integer(I4B), parameter :: NIUNIT_PRT = PRT_NBASEPKG + PRT_NMULTIPKG + contains !> @brief Create a new particle tracking model object @@ -207,7 +210,6 @@ end subroutine prt_cr subroutine prt_df(this) ! -- modules use PrtPrpModule, only: PrtPrpType - use ModelPackageInputsModule, only: NIUNIT_PRT ! -- dummy class(PrtModelType) :: this ! -- local @@ -458,7 +460,7 @@ subroutine prt_cq(this, icnvg, isuppress_output) ! conc solution. do ip = 1, this%bndlist%Count() packobj => GetBndFromList(this%bndlist, ip) - call packobj%bnd_cf(reset_mover=.false.) ! kluge note: cf call probably not needed + call packobj%bnd_cf() ! kluge note: cf call probably not needed call packobj%bnd_cq(this%x, this%flowja) end do ! @@ -962,8 +964,8 @@ subroutine allocate_arrays(this) end subroutine allocate_arrays !> @brief Create boundary condition packages for this model - subroutine package_create(this, filtyp, ipakid, ipaknum, pakname, inunit, & - iout) + subroutine package_create(this, filtyp, ipakid, ipaknum, pakname, mempath, & + inunit, iout) ! -- modules use ConstantsModule, only: LINELENGTH use PrtPrpModule, only: prp_create @@ -975,6 +977,7 @@ subroutine package_create(this, filtyp, ipakid, ipaknum, pakname, inunit, & integer(I4B), intent(in) :: ipakid integer(I4B), intent(in) :: ipaknum character(len=*), intent(in) :: pakname + character(len=*), intent(in) :: mempath integer(I4B), intent(in) :: inunit integer(I4B), intent(in) :: iout ! -- local @@ -985,9 +988,9 @@ subroutine package_create(this, filtyp, ipakid, ipaknum, pakname, inunit, & ! -- This part creates the package object select case (filtyp) case ('PRP6') - this%nprp = this%nprp + 1 ! kluge? + this%nprp = this%nprp + 1 call prp_create(packobj, ipakid, ipaknum, inunit, iout, this%name, & - pakname, this%fmi) + pakname, mempath, this%fmi) case ('API6') call api_create(packobj, ipakid, ipaknum, inunit, iout, this%name, pakname) case default @@ -1225,8 +1228,8 @@ subroutine create_bndpkgs(this, bndpkgs, pkgtypes, pkgnames, & bndptype = pkgtype end if ! - call this%package_create(pkgtype, ipakid, ipaknum, pkgname, inunit, & - this%iout) + call this%package_create(pkgtype, ipakid, ipaknum, pkgname, mempath, & + inunit, this%iout) ipakid = ipakid + 1 ipaknum = ipaknum + 1 end do @@ -1333,8 +1336,8 @@ subroutine create_packages(this) call create_methodDis(this%methoddis) call create_methodDisv(this%methoddisv) ! call create_methodDisu(this%methodDisu) - call create_methodCellPool() ! kluge - call create_methodSubcellPool() ! kluge + call create_methodCellPool() + call create_methodSubcellPool() ! ! -- Create packages that are tied directly to model call mip_cr(this%mip, this%name, mempathmip, this%inmip, this%iout, this%dis) diff --git a/src/Model/ParticleTracking/prt1prp1.f90 b/src/Model/ParticleTracking/prt1prp1.f90 index baf49cd06fb..949ee5755e4 100644 --- a/src/Model/ParticleTracking/prt1prp1.f90 +++ b/src/Model/ParticleTracking/prt1prp1.f90 @@ -19,6 +19,7 @@ module PrtPrpModule use ArrayHandlersModule, only: expandarray use GlobalDataModule use TrackModule, only: TrackControlType + use GeomUtilModule, only: point_in_polygon implicit none @@ -87,7 +88,7 @@ module PrtPrpModule !> @brief Create a new particle release package subroutine prp_create(packobj, id, ibcnum, inunit, iout, namemodel, & - pakname, fmi) + pakname, mempath, fmi) ! -- dummy class(BndType), pointer :: packobj integer(I4B), intent(in) :: id @@ -96,17 +97,22 @@ subroutine prp_create(packobj, id, ibcnum, inunit, iout, namemodel, & integer(I4B), intent(in) :: iout character(len=*), intent(in) :: namemodel character(len=*), intent(in) :: pakname + character(len=*), intent(in) :: mempath type(PrtFmiType), pointer :: fmi ! -- local type(PrtPrpType), pointer :: prpobj + ! -- formats + character(len=*), parameter :: fmtheader = & + "(1x, /1x, 'PRP -- Particle Release Point package,', & + &' INPUT READ FROM MEMPATH: ', A, /)" ! ! -- allocate the object and assign values to object variables allocate (prpobj) packobj => prpobj ! ! -- create name and memory path - call packobj%set_names(ibcnum, namemodel, pakname, ftype) - packobj%text = text + call packobj%set_names(ibcnum, namemodel, pakname, ftype, mempath) + prpobj%text = text ! ! -- allocate scalars call prpobj%prp_allocate_scalars() @@ -121,9 +127,16 @@ subroutine prp_create(packobj, id, ibcnum, inunit, iout, namemodel, & packobj%ncolbnd = 4 packobj%iscloc = 1 ! - ! -- Store pointer to flow model interface + ! -- store pointer to flow model interface prpobj%fmi => fmi ! + ! -- check if prp is enabled + if (inunit > 0) then + ! + ! -- Print a message identifying the node property flow package. + write (iout, fmtheader) mempath + end if + ! ! -- return return end subroutine prp_create @@ -359,8 +372,11 @@ subroutine prp_ad(this) integer(I4B) :: i, n, ic, icu, icpl, irow, icol, ilay integer(I4B) :: nps, np real(DP) :: rls_frac, trelease, tstop + real(DP) :: x, y logical(LGP) :: isRelease character(len=LENBOUNDNAME) :: bndName = '' + character(len=LINELENGTH) :: errmsg + real(DP), allocatable :: polyverts(:, :) isRelease = .false. ! @@ -454,7 +470,7 @@ subroutine prp_ad(this) if (this%stoptime < tstop) tstop = this%stoptime end if - ! -- Compute starting location's (user) node number and layer number + ! -- get user node number and layer number for starting location icu = this%dis%get_nodeuser(ic) select type (dis => this%dis) ! kluge??? type is (GwfDisType) @@ -463,22 +479,30 @@ subroutine prp_ad(this) call get_jk(icu, dis%ncpl, dis%nlay, icpl, ilay) end select - ! -- todo: check that location is within the specified cell + ! -- make sure the release point is in the specified cell + x = this%x(nps) + y = this%y(nps) + call this%dis%get_polyverts(ic, polyverts) + if (.not. point_in_polygon(x, y, polyverts)) then + write (errmsg, '(a,g0,a,g0,a,i0)') & + 'Error: release point (', x, ', ', y, ') is not in cell ', icu + call store_error(errmsg, terminate=.false.) + call store_error_unit(this%inunit, terminate=.true.) + end if - ! -- Set particle boundname + ! -- set particle boundname if (size(this%boundname) /= 0) & bndName = this%boundname(nps) - ! -- Add particle to particle list (todo: factor out a routine?) - ! The routine should branch based on whether this is an exchange - ! PRP or a normal PRP. If exchange PRP, particle identity info - ! should have been injected into the PRP object by the exchange - ! (e.g. this%imdl, this%iprp, this%irpt, this%trelease arrays). - ! If normal PRP, imdl and iprp should be set from the pointers - ! set on this PRP from the PRT model, and irpt and trelease set - ! as below. - this%partlist%x(np) = this%x(nps) - this%partlist%y(np) = this%y(nps) + ! -- add particle to particle list + ! -- todo: factor out a routine when implementing exchange... + ! should branch depending on exchange PRP or a normal PRP. + ! if exchange PRP, particle identities should be passed in + ! (e.g. this%imdl, this%iprp, this%irpt, this%trelease). + ! if normal PRP, imdl and iprp should be set from pointers + ! provided to PRP by PRT model; irpt and trelease as below. + this%partlist%x(np) = x + this%partlist%y(np) = y this%partlist%z(np) = this%z(nps) this%partlist%trelease(np) = trelease this%partlist%tstop(np) = tstop diff --git a/src/Solution/ParticleTracker/Method.f90 b/src/Solution/ParticleTracker/Method.f90 index c566556cace..a2e2a850375 100644 --- a/src/Solution/ParticleTracker/Method.f90 +++ b/src/Solution/ParticleTracker/Method.f90 @@ -9,6 +9,7 @@ module MethodModule private public :: MethodType + public :: get_iatop type, abstract :: MethodType ! private @@ -179,4 +180,19 @@ subroutine update(this, particle, celldefn) end subroutine update + !> @brief Get the index corresponding to top elevation of a cell in the grid. + !! This is -1 if the cell is in the top layer and 1 otherwise. + function get_iatop(ncpl, icu) result(iatop) + integer, intent(in) :: ncpl, icu + integer :: iatop + ! + if (icu .le. ncpl) then + iatop = -1 + else + iatop = 1 + end if + ! + return + end function get_iatop + end module MethodModule diff --git a/src/Solution/ParticleTracker/MethodDis.f90 b/src/Solution/ParticleTracker/MethodDis.f90 index f7c8e4b534f..a5347bdf047 100644 --- a/src/Solution/ParticleTracker/MethodDis.f90 +++ b/src/Solution/ParticleTracker/MethodDis.f90 @@ -2,9 +2,8 @@ module MethodDisModule use KindModule, only: DP, I4B use ConstantsModule, only: DONE - use MethodModule + use MethodModule, only: MethodType, get_iatop use MethodCellPoolModule - ! use CellModule use CellDefnModule use CellRectModule use ParticleModule @@ -32,16 +31,16 @@ module MethodDisModule procedure, public :: apply => apply_mGD ! applies the DIS-grid method procedure, public :: pass => pass_mGD ! passes the particle to the next cell procedure, public :: loadsub => loadsub_mGD ! loads the cell method - procedure, private :: get_npolyverts ! returns the number of polygon vertices for a cell in the grid - procedure, private :: get_iatop ! returns index used to locate top elevation of cell in the grid procedure, private :: get_top ! returns top elevation based on index iatop + + ! todo: maybe separate LoadCellDefn module? procedure, public :: load_cellDefn ! loads cell definition from the grid - procedure, public :: load_cellDefn_facenbr ! loads face neighbors to a cell object - procedure, public :: load_cellDefn_polyverts ! loads polygon vertices to a cell object - procedure, public :: load_cellDefn_flows ! loads flows to a cell object - procedure, public :: load_cellDefn_ispv180 ! loads 180-degree vertex indicator to a cell object - procedure, public :: load_cellDefn_basic ! loads basic components to a cell object from its grid - procedure, public :: addBoundaryFlows_cellRect ! adds BoundaryFlows from the grid to the faceflow array of a rectangular cell + procedure, private :: load_cellDefn_basic ! loads basic components to a cell object from its grid + procedure, private :: load_cellDefn_polyverts ! loads polygon vertices to a cell object + procedure, private :: load_cellDefn_facenbr ! loads face neighbors to a cell object + procedure, private :: load_cellDefn_ispv180 ! loads 180-degree vertex indicator to a cell object + procedure, private :: load_cellDefn_flows ! loads flows to a cell object + procedure, private :: addBoundaryFlows_cellRect ! adds BoundaryFlows from the grid to the faceflow array of a rectangular cell end type MethodDisType contains @@ -117,14 +116,15 @@ subroutine loadsub_mGD(this, particle, levelNext, submethod) ic = particle%iTrackingDomain(levelNext) ! kluge note: is cell number always known coming in? ! -- load cellDefn call this%load_cellDefn(ic, this%cellRect%cellDefn) + + ! -- If cell is active but dry, select and initialize + ! -- pass-to-bottom method and set cell method pointer if (this%fmi%ibdgwfsat0(ic) == 0) then ! kluge note: use cellDefn%sat == DZERO here instead? - ! -- Cell is active but dry, so select and initialize pass-to-bottom - ! -- cell method and set cell method pointer call methodCellPassToBot%init(particle, this%cellRect%cellDefn, & this%trackctl) submethod => methodCellPassToBot else - ! -- load rectangular cell + ! -- load rectangular cell (todo: refactor into separate routine) select type (dis => this%fmi%dis) ! kluge??? type is (GwfDisType) icu = dis%get_nodeuser(ic) @@ -156,8 +156,8 @@ subroutine loadsub_mGD(this, particle, levelNext, submethod) term = factor / areaz this%cellRect%vz1 = this%cellRect%cellDefn%faceflow(6) * term this%cellRect%vz2 = -this%cellRect%cellDefn%faceflow(7) * term - ! -- Select and initialize Pollock's cell method and set cell method - ! -- pointer + + ! -- Select and initialize Pollock's method and set method pointer call methodCellPollock%init(particle, this%cellRect, this%trackctl) submethod => methodCellPollock end if @@ -226,8 +226,6 @@ subroutine pass_mGD(this, particle) particle%iTrackingDomainBoundary(2) = inface if (inface < 5) then ! -- Map z between cells - ! iatop = this%get_iatop(ic) - ! top = this%get_top(iatop) topfrom = this%cellRect%cellDefn%top botfrom = this%cellRect%cellDefn%bot zrel = (z - botfrom) / (topfrom - botfrom) @@ -264,44 +262,6 @@ subroutine apply_mGD(this, particle, tmax) ! end subroutine apply_mGD - !> @brief Return the number of polygon vertices for a cell in the grid - function get_npolyverts(this, ic) result(npolyverts) - use GwfDisModule - implicit none - ! -- dummy - class(MethodDisType), intent(inout) :: this - integer, intent(in) :: ic - ! -- result - integer :: npolyverts - ! - npolyverts = 4 - ! - return - end function get_npolyverts - - !> @brief Get the index used to locate top elevation of a cell in the grid - function get_iatop(this, ic) result(iatop) - use GwfDisModule - implicit none - ! -- dummy - class(MethodDisType), intent(inout) :: this - integer, intent(in) :: ic - ! -- result - integer :: iatop - ! -- local - integer :: ncpl, icu - ! - ncpl = this%fmi%dis%get_ncpl() - icu = this%fmi%dis%get_nodeuser(ic) - if (icu .le. ncpl) then - iatop = -1 ! kluge note: just returns -1 if in top layer and 1 otherwise - else - iatop = 1 - end if - ! - return - end function get_iatop - !> @brief Returns a top elevation based on index iatop function get_top(this, iatop) result(top) implicit none @@ -340,15 +300,16 @@ subroutine load_cellDefn(this, ic, cellDefn) ! -- Load 180-degree indicator call this%load_cellDefn_ispv180(cellDefn) ! - ! -- Load flows - ! -- (assumes face neighbors already loaded) + ! -- Load flows (assumes face neighbors already loaded) call this%load_cellDefn_flows(cellDefn) ! return ! end subroutine load_cellDefn - !> @brief Loads basic cell definition components from the grid + !> @brief Loads basic cell definition components from the grid. + !! These include cell index, vertex count, (ia)top and botm, + !! porosity, retfactor, and izone. subroutine load_cellDefn_basic(this, ic, cellDefn) implicit none ! -- dummy @@ -359,15 +320,16 @@ subroutine load_cellDefn_basic(this, ic, cellDefn) integer :: iatop real(DP) :: top, bot, sat ! + ! -- cell index cellDefn%icell = ic ! - ! -- Load npolyverts - cellDefn%npolyverts = this%get_npolyverts(ic) + ! -- number of polygon vertices + cellDefn%npolyverts = 4 ! rectangular cell always has 4 vertices ! - ! -- Load iatop, top, and bot - iatop = this%get_iatop(ic) + ! -- iatop, top, and bot + iatop = get_iatop(this%fmi%dis%get_ncpl(), & + this%fmi%dis%get_nodeuser(ic)) cellDefn%iatop = iatop - ! cellDefn%top = this%get_top(iatop) top = this%fmi%dis%top(ic) bot = this%fmi%dis%bot(ic) sat = this%fmi%gwfsat(ic) @@ -376,7 +338,7 @@ subroutine load_cellDefn_basic(this, ic, cellDefn) cellDefn%bot = bot cellDefn%sat = sat ! - ! -- Load porosity, retfactor, and izone + ! -- porosity, retfactor, and izone cellDefn%porosity = this%porosity(ic) cellDefn%retfactor = this%retfactor(ic) cellDefn%izone = this%izone(ic) @@ -386,10 +348,7 @@ subroutine load_cellDefn_basic(this, ic, cellDefn) end subroutine load_cellDefn_basic !> @brief Load polygon vertices into cell definition from the grid - !! - !! kluge note: are polyverts even needed for MethodDis??? - !! - !< + !! Assumes cell index and number of vertices are already loaded. subroutine load_cellDefn_polyverts(this, cellDefn) use InputOutputModule ! kluge use GwfDisModule ! kluge??? @@ -404,8 +363,6 @@ subroutine load_cellDefn_polyverts(this, cellDefn) double precision :: cellx, celly, dxhalf, dyhalf ! ic = cellDefn%icell - ! - if (cellDefn%npolyverts .eq. 0) cellDefn%npolyverts = this%get_npolyverts(ic) ! kluge note: just set it to 4 ??? npolyverts = cellDefn%npolyverts ! ! -- Load polygon vertices. Note that the polyverts array @@ -440,13 +397,14 @@ subroutine load_cellDefn_polyverts(this, cellDefn) cellDefn%polyvert(5)%ivert = cellDefn%polyvert(1)%ivert cellDefn%polyvert(5)%x = cellDefn%polyvert(1)%x cellDefn%polyvert(5)%y = cellDefn%polyvert(1)%y - end select ! ... kluge??? + end select ! return ! end subroutine load_cellDefn_polyverts - !> @brief Loads face neighbors to cell definition from the grid + !> @brief Loads face neighbors to cell definition from the grid. + !! Assumes cell index and number of vertices are already loaded. subroutine load_cellDefn_facenbr(this, cellDefn) use InputOutputModule use GwfDisModule @@ -462,8 +420,6 @@ subroutine load_cellDefn_facenbr(this, cellDefn) integer :: ncpl ! ic = cellDefn%icell - ! - if (cellDefn%npolyverts .eq. 0) cellDefn%npolyverts = this%get_npolyverts(ic) npolyverts = cellDefn%npolyverts ! ! -- Load face neighbors. Note that the facenbr array @@ -516,12 +472,9 @@ subroutine load_cellDefn_facenbr(this, cellDefn) ! end subroutine load_cellDefn_facenbr - !> @brief Load flows into the cell definition - !! - !! Loads polygon face, top and bottom, and net distributed - !! flows from the grid into the cell definition. - !! - !> + !> @brief Load flows into the cell definition. + !! These include face flows and net distributed flows. + !! Assumes cell index and number of vertices are already loaded. subroutine load_cellDefn_flows(this, cellDefn) implicit none ! -- dummy @@ -531,8 +484,6 @@ subroutine load_cellDefn_flows(this, cellDefn) integer :: ic, npolyverts, m, n ! ic = cellDefn%icell - ! - if (cellDefn%npolyverts .eq. 0) cellDefn%npolyverts = this%get_npolyverts(ic) npolyverts = cellDefn%npolyverts ! ! -- Load face flows. Note that the faceflow array @@ -574,7 +525,8 @@ subroutine load_cellDefn_flows(this, cellDefn) ! end subroutine load_cellDefn_flows - !> @brief Add BoundaryFlows to the cell definition faceflow array + !> @brief Add boundary flows to the cell definition faceflow array. + !! Assumes cell index and number of vertices are already loaded. subroutine addBoundaryFlows_cellRect(this, cellDefn) implicit none ! -- dummy @@ -606,13 +558,9 @@ subroutine addBoundaryFlows_cellRect(this, cellDefn) ! end subroutine addBoundaryFlows_cellRect - !> @brief Load vertex 180-degree indicator array into cell definition - !! - !! Loads 180-degree vertex indicator array to cell - !! definition and sets flags that indicate how cell - !! can be represented. Todo: rename? - !! - !< + !> @brief Load 180-degree vertex indicator array and set flags + !! indicating how cell can be represented (kluge: latter needed?). + !! Assumes cell index and number of vertices are already loaded. subroutine load_cellDefn_ispv180(this, cellDefn) implicit none ! -- dummy @@ -621,16 +569,15 @@ subroutine load_cellDefn_ispv180(this, cellDefn) ! -- local integer :: npolyverts ! - ! -- Set up polyverts if not done previously - if (cellDefn%npolyverts .eq. 0) call this%load_cellDefn_polyverts(cellDefn) npolyverts = cellDefn%npolyverts ! ! -- Load 180-degree indicator. Note that the ispv180 array ! -- does not get reallocated if it is already allocated ! -- to a size greater than or equal to npolyverts+1. - ! -- Also, set flags that indicate how cell can be represented. call allocate_as_needed(cellDefn%ispv180, npolyverts + 1) cellDefn%ispv180(1:npolyverts + 1) = .false. + ! + ! -- Set flags that indicate how cell can be represented. cellDefn%canBeCellRect = .true. cellDefn%canBeCellRectQuad = .false. ! diff --git a/src/Solution/ParticleTracker/MethodDisv.f90 b/src/Solution/ParticleTracker/MethodDisv.f90 index a78ed5d08c5..5ca6890c36a 100644 --- a/src/Solution/ParticleTracker/MethodDisv.f90 +++ b/src/Solution/ParticleTracker/MethodDisv.f90 @@ -2,7 +2,7 @@ module MethodDisvModule use KindModule, only: DP, I4B use ConstantsModule, only: DONE - use MethodModule + use MethodModule, only: MethodType, get_iatop use MethodCellPoolModule use CellDefnModule use CellPolyModule @@ -33,20 +33,18 @@ module MethodDisvModule procedure, public :: loadsub => loadsub_mGDv ! loads the cell method procedure, public :: mapToNbrCell ! maps a location on the cell face to the shared face of a neighbor procedure, private :: get_npolyverts ! returns the number of polygon vertices for a cell in the grid - procedure, private :: get_iatop ! returns index used to locate top elevation of cell in the grid procedure, private :: get_top ! returns top elevation based on index iatop + + ! todo: maybe separate LoadCellDefn module? procedure, public :: load_cellDefn ! loads cell definition from the grid - procedure, public :: load_cellDefn_facenbr ! loads face neighbors to a cell object - procedure, public :: load_cellDefn_polyverts ! loads polygon vertices to a cell object - procedure, public :: load_cellDefn_flows ! loads flows to a cell object - procedure, public :: load_cellDefn_ispv180 ! loads 180-degree vertex indicator to a cell object - procedure, public :: load_cellDefn_basic ! loads basic components to a cell object from its grid - ! adds BoundaryFlows from grid to faceflow array of a rectangular cell - procedure, public :: addBoundaryFlows_cellRect - ! adds BoundaryFlows from the grid to the faceflow array of a rectangular-quad cell - procedure, public :: addBoundaryFlows_cellRectQuad - ! adds BoundaryFlows from the grid to the faceflow array of a polygonal cell - procedure, public :: addBoundaryFlows_cellPoly + procedure, private :: load_cellDefn_basic ! loads basic components to a cell object from its grid + procedure, private :: load_cellDefn_polyverts ! loads polygon vertices to a cell object + procedure, private :: load_cellDefn_facenbr ! loads face neighbors to a cell object + procedure, private :: load_cellDefn_ispv180 ! loads 180-degree vertex indicator to a cell object + procedure, private :: load_cellDefn_flows ! loads flows to a cell object + procedure, private :: addBoundaryFlows_cellRect ! adds BoundaryFlows from the grid to the faceflow array of a rectangular cell + procedure, private :: addBoundaryFlows_cellRectQuad ! adds BoundaryFlows from the grid to the faceflow array of a rectangular-quad cell + procedure, private :: addBoundaryFlows_cellPoly ! adds BoundaryFlows from the grid to the faceflow array of a polygonal cell end type MethodDisvType contains @@ -256,8 +254,6 @@ subroutine mapToNbrCell(this, cellDefnin, inface, z) end if end do ! -- Map z between cells - ! iatop = this%get_iatop(ic) - ! top = this%get_top(iatop) topfrom = cellDefnin%top botfrom = cellDefnin%bot zrel = (z - botfrom) / (topfrom - botfrom) @@ -314,29 +310,6 @@ function get_npolyverts(this, ic) result(npolyverts) return end function get_npolyverts - !> @brief Get index used to locate top elevation of a cell in the grid - function get_iatop(this, ic) result(iatop) - use GwfDisvModule ! kluge??? - implicit none - ! -- dummy - class(MethodDisvType), intent(inout) :: this - integer, intent(in) :: ic - ! -- result - integer :: iatop - ! -- local - integer :: ncpl, icu - ! - ncpl = this%fmi%dis%get_ncpl() - icu = this%fmi%dis%get_nodeuser(ic) - if (icu .le. ncpl) then - iatop = -1 ! kluge note: just returns -1 if in top layer and 1 otherwise - else - iatop = 1 - end if - ! - return - end function get_iatop - !> @brief Get top elevation based on index iatop !! kluge note: not needed??? function get_top(this, iatop) result(top) @@ -376,8 +349,7 @@ subroutine load_cellDefn(this, ic, cellDefn) ! -- Load 180-degree indicator call this%load_cellDefn_ispv180(cellDefn) ! - ! -- Load flows - ! -- (assumes face neighbors already loaded) + ! -- Load flows (assumes face neighbors already loaded) call this%load_cellDefn_flows(cellDefn) ! return @@ -385,6 +357,8 @@ subroutine load_cellDefn(this, ic, cellDefn) end subroutine load_cellDefn !> @brief Loads basic cell definition components from the grid + !! These include cell index, vertex count, (ia)top and botm, + !! porosity, retfactor, and izone. subroutine load_cellDefn_basic(this, ic, cellDefn) implicit none ! -- dummy @@ -396,13 +370,15 @@ subroutine load_cellDefn_basic(this, ic, cellDefn) integer :: iatop real(DP) :: top, bot, sat ! + ! -- cell index cellDefn%icell = ic ! - ! -- Load npolyverts + ! -- number of polygon vertices cellDefn%npolyverts = this%get_npolyverts(ic) ! - ! -- Load iatop, top, and bot - iatop = this%get_iatop(ic) + ! -- iatop, top, and bot + iatop = get_iatop(this%fmi%dis%get_ncpl(), & + this%fmi%dis%get_nodeuser(ic)) cellDefn%iatop = iatop top = this%fmi%dis%top(ic) bot = this%fmi%dis%bot(ic) @@ -422,6 +398,7 @@ subroutine load_cellDefn_basic(this, ic, cellDefn) end subroutine load_cellDefn_basic !> @brief Loads polygon vertices to cell definition from the grid + !! Assumes cell index and number of vertices are already loaded. subroutine load_cellDefn_polyverts(this, cellDefn) use GwfDisvModule implicit none @@ -433,8 +410,6 @@ subroutine load_cellDefn_polyverts(this, cellDefn) integer :: ncpl ! kluge??? ! ic = cellDefn%icell - ! - if (cellDefn%npolyverts .eq. 0) cellDefn%npolyverts = this%get_npolyverts(ic) npolyverts = cellDefn%npolyverts ! ! -- Load polygon vertices. Note that the polyverts array @@ -466,6 +441,7 @@ subroutine load_cellDefn_polyverts(this, cellDefn) end subroutine load_cellDefn_polyverts !> @brief Loads face neighbors to cell definition from the grid + !! Assumes cell index and number of vertices are already loaded. subroutine load_cellDefn_facenbr(this, cellDefn) use InputOutputModule use GwfDisvModule @@ -480,8 +456,6 @@ subroutine load_cellDefn_facenbr(this, cellDefn) integer :: ncpl ! ic = cellDefn%icell - ! - if (cellDefn%npolyverts .eq. 0) cellDefn%npolyverts = this%get_npolyverts(ic) npolyverts = cellDefn%npolyverts ! ! -- Load face neighbors. Note that the facenbr array @@ -568,12 +542,9 @@ subroutine shared_edgeface(ivlist1, ivlist2, iedgeface) end do outerloop end subroutine shared_edgeface - !> @brief Load flows into the cell from the grid - !! - !! Load polygon face, top and bottom, and net distributed - !! flows to cell definition from the grid. - !! - !< + !> @brief Load flows into the cell definition. + !! These include face flows and net distributed flows. + !! Assumes cell index and number of vertices are already loaded. subroutine load_cellDefn_flows(this, cellDefn) implicit none ! -- dummy @@ -583,8 +554,6 @@ subroutine load_cellDefn_flows(this, cellDefn) integer :: ic, npolyverts, m, n ! ic = cellDefn%icell - ! - if (cellDefn%npolyverts .eq. 0) cellDefn%npolyverts = this%get_npolyverts(ic) npolyverts = cellDefn%npolyverts ! ! -- Load face flows. Note that the faceflow array @@ -624,7 +593,8 @@ subroutine load_cellDefn_flows(this, cellDefn) ! end subroutine load_cellDefn_flows - !> @brief Load boundary flows from the grid into a rectangular cell + !> @brief Load boundary flows from the grid into a rectangular cell. + !! Assumes cell index and number of vertices are already loaded. subroutine addBoundaryFlows_cellRect(this, cellDefn) implicit none ! -- dummy @@ -660,7 +630,8 @@ subroutine addBoundaryFlows_cellRect(this, cellDefn) ! end subroutine addBoundaryFlows_cellRect - !> @brief Load boundary flows from the grid into rectangular quadcell + !> @brief Load boundary flows from the grid into rectangular quadcell. + !! Assumes cell index and number of vertices are already loaded. subroutine addBoundaryFlows_cellRectQuad(this, cellDefn) implicit none ! -- dummy @@ -727,7 +698,8 @@ subroutine addBoundaryFlows_cellRectQuad(this, cellDefn) ! end subroutine addBoundaryFlows_cellRectQuad - !> @brief Load boundary flows from the grid into a polygonal cell + !> @brief Load boundary flows from the grid into a polygonal cell. + !! Assumes cell index and number of vertices are already loaded. subroutine addBoundaryFlows_cellPoly(this, cellDefn) implicit none ! -- dummy @@ -760,20 +732,15 @@ subroutine addBoundaryFlows_cellPoly(this, cellDefn) ! end subroutine addBoundaryFlows_cellPoly - !> @brief Load 180-degree vertex indicator array into cell. - !! - !! Loads 180-degree vertex indicator array to cell - !! definition and sets flags that indicate how cell - !! can be represented - !! - !< + !> @brief Load 180-degree vertex indicator array and set flags + !! indicating how cell can be represented (kluge: latter needed?). + !! Assumes cell index and number of vertices are already loaded. subroutine load_cellDefn_ispv180(this, cellDefn) ! kluge note: rename??? implicit none ! -- dummy class(MethodDisvType), intent(inout) :: this type(CellDefnType), pointer, intent(inout) :: cellDefn ! -- local - ! integer :: ig,ic,npolyverts,iapnbr,iapv180 integer :: npolyverts, m, m0, m1, m2 integer :: ic integer :: num90, num180, numacute @@ -783,9 +750,6 @@ subroutine load_cellDefn_ispv180(this, cellDefn) ! kluge note: rename??? logical last180 ! ic = cellDefn%icell - ! - ! -- Set up polyverts if not done previously - if (cellDefn%npolyverts .eq. 0) call this%load_cellDefn_polyverts(cellDefn) npolyverts = cellDefn%npolyverts ! ! -- Load 180-degree indicator. Note that the ispv180 array diff --git a/src/Utilities/GeomUtil.f90 b/src/Utilities/GeomUtil.f90 index f6014f56c23..93d7a30a6e2 100644 --- a/src/Utilities/GeomUtil.f90 +++ b/src/Utilities/GeomUtil.f90 @@ -1,5 +1,5 @@ module GeomUtilModule - use KindModule, only: I4B, DP + use KindModule, only: I4B, DP, LGP implicit none private public :: between, point_in_polygon @@ -12,8 +12,12 @@ logical function between(x, a, b) end function between !> @brief Check if a point is within a polygon. - !! Vertices and edge points are considered in. - !! Reference: https://stackoverflow.com/a/63436180/6514033 + !! + !! Vertices and edge points are considered in the polygon. + !! + !! Adapted from https://stackoverflow.com/a/63436180/6514033, + !! modified to check for approximate floating point equality. + !< logical function point_in_polygon(x, y, poly) ! dummy real(DP), intent(in) :: x !< x point coordinate @@ -33,23 +37,27 @@ logical function point_in_polygon(x, y, poly) xb = poly(1, ii) yb = poly(2, ii) - if ((x == xa .and. y == ya) .or. (x == xb .and. y == yb)) then + if ((x == xa .and. y == ya) .or. & + (x == xb .and. y == yb)) then ! on vertex point_in_polygon = .true. exit - else if (ya == yb .and. y == ya .and. between(x, xa, xb)) then + else if (ya == yb .and. & + y == ya .and. & + between(x, xa, xb)) then ! on horizontal edge point_in_polygon = .true. exit else if (between(y, ya, yb)) then - if ((y == ya .and. yb >= ya) .or. (y == yb .and. ya >= yb)) then + if ((y == ya .and. yb >= ya) .or. & + (y == yb .and. ya >= yb)) then xa = xb ya = yb cycle end if ! cross product c = (xa - x) * (yb - y) - (xb - x) * (ya - y) - if (c == 0) then + if (c == 0.0_DP) then ! on edge point_in_polygon = .true. exit