From a240440c7aa0e4170a78370fd2227f7e40ff8eca Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dav=C3=ADd=20Brakenhoff?= Date: Tue, 24 Dec 2024 14:33:46 +0100 Subject: [PATCH 01/10] Add PRT support - add prt submodule with PRT packages (missing FMI at the moment) --- nlmod/prt/__init__.py | 3 + nlmod/prt/prt.py | 215 ++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 218 insertions(+) create mode 100644 nlmod/prt/__init__.py create mode 100644 nlmod/prt/prt.py diff --git a/nlmod/prt/__init__.py b/nlmod/prt/__init__.py new file mode 100644 index 00000000..e7d3a18b --- /dev/null +++ b/nlmod/prt/__init__.py @@ -0,0 +1,3 @@ +# ruff: noqa: F401 F403 +from . import prt +from .prt import * diff --git a/nlmod/prt/prt.py b/nlmod/prt/prt.py new file mode 100644 index 00000000..d77987b1 --- /dev/null +++ b/nlmod/prt/prt.py @@ -0,0 +1,215 @@ +import logging +from pathlib import Path + +import flopy as fp + +from ..gwf.gwf import _dis, _disv, _set_record +from ..util import _get_value_from_ds_attr, _get_value_from_ds_datavar + +logger = logging.getLogger(__name__) + + +def prt(ds, sim, modelname=None, save_flows=True, **kwargs): + """Create particle tracking model from the model dataset. + + Parameters + ---------- + ds : xarray.Dataset + dataset with model data. Should have the dimension 'time' and the + attributes: model_name, mfversion, model_ws, time_units, start, + perlen, nstp, tsmult + sim : flopy MFSimulation + simulation object. + modelname : str + name of the particle tracking model + + Returns + ------- + prt : flopy ModflowPrt + particle tracking model object + """ + # start creating model + logger.info("creating mf6 PRT") + + if modelname is None: + modelname = f"{ds.model_name}_prt" + model_nam_file = f"{modelname}.nam" + + prt = fp.mf6.ModflowPrt( + sim, + modelname=modelname, + model_nam_file=model_nam_file, + save_flows=save_flows, + **kwargs, + ) + return prt + + +def dis(ds, prt, length_units="METERS", pname="dis", **kwargs): + """Create discretisation package from the model dataset. + + Parameters + ---------- + ds : xarray.Dataset + dataset with model data. + prt : flopy ModflowPrt + particle tracking object + length_units : str, optional + length unit. The default is 'METERS'. + pname : str, optional + package name + + Returns + ------- + dis : flopy ModflowPrtdis + discretisation package. + """ + return _dis(ds, prt, length_units, pname, **kwargs) + + +def disv(ds, prt, length_units="METERS", pname="disv", **kwargs): + """Create discretisation vertices package from the model dataset. + + Parameters + ---------- + ds : xarray.Dataset + dataset with model data. + prt : flopy ModflowPrt + particle tracking object + length_units : str, optional + length unit. The default is 'METERS'. + pname : str, optional + package name + + Returns + ------- + disv : flopy ModflowPrtdisv + disv package + """ + return _disv(ds, prt, length_units, pname, **kwargs) + + +def mip(ds, prt, porosity=None, **kwargs): + """Create matrix input package for groundwater transport model. + + Parameters + ---------- + ds : xarray.Dataset + dataset with model data. + prt : ModflowPrt + modflow particle tracking object + porosity : str, float, array + porosity array/attribute name, value or array. The default is None. + + Returns + ------- + mip : flopy ModflowGwtmip + matrix input package + """ + logger.info("creating mf6 MIP") + + # NOTE: attempting to look for porosity in attributes first, then data variables. + # If both are defined, the attribute value will be used. The log message in this + # case is not entirely correct. This is something we may need to sort out, and + # also think about the order we do this search. + if isinstance(porosity, str): + value = None + else: + value = porosity + porosity = _get_value_from_ds_attr( + ds, "porosity", porosity, value=value, warn=False + ) + porosity = _get_value_from_ds_datavar(ds, "porosity", porosity) + return fp.mf6.ModflowPrtmip(prt, pname="mip", porosity=porosity, **kwargs) + + +def prp(ds, prt, packagedata, perioddata, pname="prp", **kwargs): + """Create particle release package for groundwater transport model. + + Parameters + ---------- + ds : xarray.Dataset + dataset with model data + prt : ModflowPrt + modflow particle tracking object + releasepts : list + list of tuples containing release point data + + Returns + ------- + prp : flopy ModflowGwtprp + particle release package + """ + logger.info("creating mf6 PRP") + prp_track_file = kwargs.pop("prp_track_file", f"{ds.model_name}.prp.trk") + prp = fp.mf6.ModflowPrtprp( + prt, + pname=pname, + filename=f"{ds.model_name}_prt.prp", + nreleasepts=len(packagedata), + packagedata=packagedata, + perioddata=perioddata, + track_filerecord=[prp_track_file], + stop_at_weak_sink=kwargs.pop("stop_at_weak_sink", False), + boundnames=kwargs.pop("boundnames", False), + exit_solve_tolerance=kwargs.pop("exit_solve_tolerance", 1e-5), + **kwargs, + ) + return prp + + +def fmi(ds, prt, packagedata=None, **kwargs): + """Create flow model interface package for groundwater transport model. + + Parameters + ---------- + ds : xarray.Dataset + dataset with model data + prt : ModflowPrt + modflow particle tracking object + packagedata : list, optional + list of tuples containing package name and file path to heads and budget files. + The default is None, which will derive heads and buget file names from the model + dataset. + **kwargs : dict + additional keyword arguments passed to flopy.mf6.ModflowPrtfmi() + + Returns + ------- + fmi : flopy ModflowGwtfmi + flow model interface package + + """ + logger.info("creating mf6 FMI") + + if packagedata is None: + gwf_head_file = Path(ds.model_ws) / f"{ds.model_name}.hds" + gwf_budget_file = Path(ds.model_ws) / f"{ds.model_name}.cbc" + packagedata = [ + ("GWFHEAD", gwf_head_file.absolute().resolve()), + ("GWFBUDGET", gwf_budget_file.absolute().resolve()), + ] + + fmi = fp.mf6.ModflowPrtfmi(prt, packagedata=packagedata, pname="fmi", **kwargs) + return fmi + + +def oc(ds, prt, save_budget=True, print_budget=False, **kwargs): + logger.info("creating mf6 OC") + + budget_filerecord = kwargs.pop("budget_filerecord", [f"{ds.model_name}_prt.cbc"]) + track_filerecord = kwargs.pop("track_filerecord", [f"{ds.model_name}_prt.trk"]) + + saverecord = _set_record(False, save_budget, output="budget") + printrecord = _set_record(False, print_budget, output="concentration") + oc = fp.mf6.ModflowPrtoc( + prt, + pname="oc", + filename=f"{ds.model_name}_prt.oc", + saverecord=saverecord, + printrecord=printrecord, + budget_filerecord=budget_filerecord, + track_filerecord=track_filerecord, + **kwargs, + ) + return oc From fb37bf179e343d02befef43da76193d03ea7991c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dav=C3=ADd=20Brakenhoff?= Date: Tue, 24 Dec 2024 14:34:17 +0100 Subject: [PATCH 02/10] modify _dis(v) pkgs to support Prtdis(v) packages --- nlmod/gwf/gwf.py | 133 +++++++++++++++++++++-------------------------- 1 file changed, 58 insertions(+), 75 deletions(-) diff --git a/nlmod/gwf/gwf.py b/nlmod/gwf/gwf.py index 215383a1..241232f1 100644 --- a/nlmod/gwf/gwf.py +++ b/nlmod/gwf/gwf.py @@ -114,46 +114,34 @@ def _dis(ds, model, length_units="METERS", pname="dis", **kwargs): idomain = get_idomain(ds).data if model.model_type == "gwf6": - dis = flopy.mf6.ModflowGwfdis( - model, - pname=pname, - length_units=length_units, - xorigin=xorigin, - yorigin=yorigin, - angrot=angrot, - nlay=ds.sizes["layer"], - nrow=ds.sizes["y"], - ncol=ds.sizes["x"], - delr=get_delr(ds), - delc=get_delc(ds), - top=ds["top"].data, - botm=ds["botm"].data, - idomain=idomain, - filename=f"{ds.model_name}.dis", - **kwargs, - ) + filename = f"{ds.model_name}.dis" + klass = flopy.mf6.ModflowGwfdis elif model.model_type == "gwt6": - dis = flopy.mf6.ModflowGwtdis( - model, - pname=pname, - length_units=length_units, - xorigin=xorigin, - yorigin=yorigin, - angrot=angrot, - nlay=ds.sizes["layer"], - nrow=ds.sizes["y"], - ncol=ds.sizes["x"], - delr=get_delr(ds), - delc=get_delc(ds), - top=ds["top"].data, - botm=ds["botm"].data, - idomain=idomain, - filename=f"{ds.model_name}_gwt.dis", - **kwargs, - ) + filename = (f"{ds.model_name}_gwt.dis",) + klass = flopy.mf6.ModflowGwtdis + elif model.model_type == "prt6": + filename = f"{ds.model_name}_prt.dis" + klass = flopy.mf6.ModflowPrtdis else: raise ValueError("Unknown model type.") - + dis = klass( + model, + pname=pname, + length_units=length_units, + xorigin=xorigin, + yorigin=yorigin, + angrot=angrot, + nlay=ds.sizes["layer"], + nrow=ds.sizes["y"], + ncol=ds.sizes["x"], + delr=get_delr(ds), + delc=get_delc(ds), + top=ds["top"].data, + botm=ds["botm"].data, + idomain=idomain, + filename=filename, + **kwargs, + ) return dis @@ -213,45 +201,34 @@ def _disv(ds, model, length_units="METERS", pname="disv", **kwargs): cell2d = grid.get_cell2d_from_ds(ds) idomain = get_idomain(ds).data if model.model_type == "gwf6": - disv = flopy.mf6.ModflowGwfdisv( - model, - idomain=idomain, - xorigin=xorigin, - yorigin=yorigin, - length_units=length_units, - angrot=angrot, - nlay=len(ds.layer), - ncpl=len(ds.icell2d), - nvert=len(ds.iv), - top=ds["top"].data, - botm=ds["botm"].data, - vertices=vertices, - cell2d=cell2d, - pname=pname, - **kwargs, - ) + klass = flopy.mf6.ModflowGwfdisv + filename = f"{ds.model_name}.disv" elif model.model_type == "gwt6": - disv = flopy.mf6.ModflowGwtdisv( - model, - idomain=idomain, - xorigin=xorigin, - yorigin=yorigin, - length_units=length_units, - angrot=angrot, - nlay=len(ds.layer), - ncpl=len(ds.icell2d), - nvert=len(ds.iv), - top=ds["top"].data, - botm=ds["botm"].data, - vertices=vertices, - cell2d=cell2d, - pname=pname, - filename=f"{ds.model_name}_gwt.disv", - **kwargs, - ) + klass = flopy.mf6.ModflowGwtdisv + filename = f"{ds.model_name}.disv" + elif model.model_type == "prt6": + klass = flopy.mf6.ModflowPrtdisv + filename = (f"{ds.model_name}_prt.disv",) else: raise ValueError("Unknown model type.") - + disv = klass( + model, + idomain=idomain, + xorigin=xorigin, + yorigin=yorigin, + length_units=length_units, + angrot=angrot, + nlay=len(ds.layer), + ncpl=len(ds.icell2d), + nvert=len(ds.iv), + top=ds["top"].data, + botm=ds["botm"].data, + vertices=vertices, + cell2d=cell2d, + pname=pname, + filename=filename, + **kwargs, + ) if "angrot" in ds.attrs and ds.attrs["angrot"] != 0.0: model.modelgrid.set_coord_info(xoff=xorigin, yoff=yorigin, angrot=angrot) @@ -875,7 +852,7 @@ def _set_record(out, budget, output="head"): return record -def buy(ds, gwf, pname="buy", **kwargs): +def buy(ds, gwf, pname="buy", gwt_modelname=None, **kwargs): """Create buoyancy package from model dataset. Parameters @@ -886,6 +863,9 @@ def buy(ds, gwf, pname="buy", **kwargs): groundwaterflow object. pname : str, optional package name, by default "buy" + gwt_modelname : str, optional + name of the groundwater transport model, by default None, + which uses gwf modelname with '_gwt' appended. Returns ------- @@ -903,6 +883,7 @@ def buy(ds, gwf, pname="buy", **kwargs): "BUY package requires a groundwater transport model. " "Set 'transport' to True in model dataset." ) + logger.info("creating mf6 BUY") drhodc = _get_value_from_ds_attr( ds, "drhodc", attr="drhodc", value=kwargs.pop("drhodc", None) @@ -914,7 +895,9 @@ def buy(ds, gwf, pname="buy", **kwargs): ds, "denseref", attr="denseref", value=kwargs.pop("denseref", None) ) - pdata = [(0, drhodc, crhoref, f"{ds.model_name}_gwt", "CONCENTRATION")] + if gwt_modelname is None: + gwt_modelname = f"{ds.model_name}_gwt" + pdata = [(0, drhodc, crhoref, gwt_modelname, "CONCENTRATION")] buy = flopy.mf6.ModflowGwfbuy( gwf, From 2d00838a87011a0b5d93d29aff44e1a3c67d956a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dav=C3=ADd=20Brakenhoff?= Date: Tue, 24 Dec 2024 14:34:28 +0100 Subject: [PATCH 03/10] improve docstring --- nlmod/mfoutput/mfoutput.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/nlmod/mfoutput/mfoutput.py b/nlmod/mfoutput/mfoutput.py index 40b83041..9b1c8f45 100644 --- a/nlmod/mfoutput/mfoutput.py +++ b/nlmod/mfoutput/mfoutput.py @@ -252,8 +252,7 @@ def _get_budget_da( def _get_flopy_data_object( var, ds=None, gwml=None, fname=None, grb_file=None, **kwargs ): - """Get modflow HeadFile or CellBudgetFile object, containg heads, budgets or - concentrations. + """Get flopy data object, containing heads, budgets or concentrations. Provide one of ds, gwf or fname. From 56c2789ccabda62d157c9dce52556dd3ea739c4b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dav=C3=ADd=20Brakenhoff?= Date: Tue, 24 Dec 2024 14:35:38 +0100 Subject: [PATCH 04/10] allow package_to_nodes to be used without Modpath model - deprecate mpf, add ibound kwarg - add structured kwarg to pg_from_pd (needed for structured models) --- nlmod/modpath/modpath.py | 45 +++++++++++++++++++++++++++++----------- 1 file changed, 33 insertions(+), 12 deletions(-) diff --git a/nlmod/modpath/modpath.py b/nlmod/modpath/modpath.py index fe1a135d..21d32bb5 100644 --- a/nlmod/modpath/modpath.py +++ b/nlmod/modpath/modpath.py @@ -2,6 +2,7 @@ import logging import numbers import os +import warnings from shutil import copyfile import flopy @@ -10,7 +11,7 @@ from packaging.version import parse as parse_version from .. import util -from ..dims.grid import get_row_col_from_xy, get_icell2d_from_xy +from ..dims.grid import get_icell2d_from_xy, get_row_col_from_xy logger = logging.getLogger(__name__) @@ -105,7 +106,7 @@ def get_node_vertex(lay, icell2d, shape): return lay * shape[1] + icell2d -def package_to_nodes(gwf, package_name, mpf): +def package_to_nodes(gwf, package_name, mpf=None, ibound=None): """Return a list of nodes from the cells with certain boundary conditions. Parameters @@ -114,8 +115,8 @@ def package_to_nodes(gwf, package_name, mpf): Groundwater flow model. package_name : str name of the package. - mpf : flopy.modpath.mp7.Modpath7 - modpath object. + ibound : array + array indicating active cells Raises ------ @@ -127,6 +128,13 @@ def package_to_nodes(gwf, package_name, mpf): nodes : list of ints node numbers corresponding to the cells with a certain boundary condition. """ + if mpf is not None: + warnings.warn( + "The 'mpf' parameter is deprecated and will be removed in a future version. " + "Please pass 'ibound' directly.", + DeprecationWarning, + ) + ibound = mpf.ib gwf_package = gwf.get_package(package_name) if hasattr(gwf_package, "stress_period_data"): pkg_cid = gwf_package.stress_period_data.array[0]["cellid"] @@ -138,17 +146,28 @@ def package_to_nodes(gwf, package_name, mpf): ) nodes = [] for cid in pkg_cid: - if len(mpf.ib.shape) == 3: - if mpf.ib[cid[0], cid[1], cid[2]] > 0: - nodes.append(get_node_structured(cid[0], cid[1], cid[2], mpf.ib.shape)) + if ibound is None: + if gwf.modelgrid.grid_type == "structured": + nodes.append( + get_node_structured(cid[0], cid[1], cid[2], gwf.modelgrid.shape) + ) + elif gwf.modelgrid.grid_type == "vertex": + nodes.append(get_node_vertex(cid[0], cid[1], gwf.modelgrid.shape)) + else: + raise NotImplementedError( + "only structured and vertex grids are supported" + ) + elif len(ibound.shape) == 3: + if ibound[cid[0], cid[1], cid[2]] > 0: + nodes.append(get_node_structured(cid[0], cid[1], cid[2], ibound.shape)) else: - if mpf.ib[cid[0], cid[1]] > 0: - nodes.append(get_node_vertex(cid[0], cid[1], mpf.ib.shape)) + if ibound[cid[0], cid[1]] > 0: + nodes.append(get_node_vertex(cid[0], cid[1], ibound.shape)) return nodes def layer_to_nodes(mpf, modellayer): - """Get the nodes of all cells in one ore more model layer(s). + """Get the nodes of all cells in one or more model layer(s). Parameters ---------- @@ -403,7 +422,7 @@ def pg_from_fdt(nodes, divisions=3): return pg -def pg_from_pd(nodes, localx=0.5, localy=0.5, localz=0.5): +def pg_from_pd(nodes, localx=0.5, localy=0.5, localz=0.5, structured=False): """Create a particle group using the ParticleData. Parameters @@ -429,6 +448,8 @@ def pg_from_pd(nodes, localx=0.5, localy=0.5, localz=0.5): be provided for each partloc. If localz is None, a value of 0.5 (center of the cell) will be used (default is None). A localz value of 1.0 indicates the top of a cell. + structured : bool, optional + if True, assumes structured model grid. Returns ------- @@ -437,7 +458,7 @@ def pg_from_pd(nodes, localx=0.5, localy=0.5, localz=0.5): """ p = flopy.modpath.ParticleData( partlocs=nodes, - structured=False, + structured=structured, localx=localx, localy=localy, localz=localz, From 24abb67b2a34d047be7429abfcc3229a6e590240 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dav=C3=ADd=20Brakenhoff?= Date: Tue, 24 Dec 2024 14:36:03 +0100 Subject: [PATCH 05/10] import sim submodule --- nlmod/sim/__init__.py | 1 + 1 file changed, 1 insertion(+) diff --git a/nlmod/sim/__init__.py b/nlmod/sim/__init__.py index 57b71d59..6cc70a87 100644 --- a/nlmod/sim/__init__.py +++ b/nlmod/sim/__init__.py @@ -1,2 +1,3 @@ # ruff: noqa: F403 +from . import sim from .sim import * From b8d898dffcb5f5a5bb53582d9c5b1cb0fc5719c7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dav=C3=ADd=20Brakenhoff?= Date: Tue, 24 Dec 2024 14:36:56 +0100 Subject: [PATCH 06/10] improve solvers - allow directly registering solution pkgs if model is passed - add EMS package --- nlmod/sim/sim.py | 28 +++++++++++++++++++++++++--- 1 file changed, 25 insertions(+), 3 deletions(-) diff --git a/nlmod/sim/sim.py b/nlmod/sim/sim.py index 6480dbca..2f9c21aa 100644 --- a/nlmod/sim/sim.py +++ b/nlmod/sim/sim.py @@ -216,8 +216,8 @@ def tdis(ds, sim, pname="tdis", nstp="nstp", tsmult="tsmult", **kwargs): return tdis -def ims(sim, complexity="MODERATE", pname="ims", **kwargs): - """Create IMS package. +def ims(sim, complexity="MODERATE", pname="ims", model=None, **kwargs): + """Create implicit model solution (IMS) package. Parameters ---------- @@ -245,9 +245,31 @@ def ims(sim, complexity="MODERATE", pname="ims", **kwargs): complexity=complexity, **kwargs, ) - + if model is not None: + register_solution_package(sim, model, ims) return ims +def ems(sim, pname="ems", model=None, **kwargs): + """Create explicit model solution (EMS) package. + + Parameters + ---------- + sim : flopy MFSimulation + simulation object. + pname : str, optional + package name + + """ + ems = flopy.mf6.ModflowEms(sim, pname=pname, **kwargs) + if model is not None: + register_solution_package(sim, model, ems) + return ems + + def register_ims_package(sim, model, ims): sim.register_ims_package(ims, [model.name]) + + +def register_solution_package(sim, model, solver): + sim.register_solution_package(solver, [model.name]) From 8ca772885f31aa382b09166d2b96f466580cc50f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dav=C3=ADd=20Brakenhoff?= Date: Tue, 24 Dec 2024 14:37:03 +0100 Subject: [PATCH 07/10] import prt --- nlmod/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nlmod/__init__.py b/nlmod/__init__.py index 5582151b..05b33ca0 100644 --- a/nlmod/__init__.py +++ b/nlmod/__init__.py @@ -3,7 +3,7 @@ NLMOD_DATADIR = os.path.join(os.path.dirname(__file__), "data") -from . import dims, gis, gwf, gwt, modpath, plot, read, sim, util +from . import dims, gis, gwf, gwt, modpath, plot, prt, read, sim, util from .dims import base, get_ds, grid, layers, resample, time, to_model_ds from .util import download_mfbinaries from .version import __version__, show_versions From 20cd7ffa90ed1f5f4be4c9ce80515c64bd817f9c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dav=C3=ADd=20Brakenhoff?= Date: Tue, 24 Dec 2024 14:37:12 +0100 Subject: [PATCH 08/10] add test notebook --- docs/examples/18_particle_tracking_prt.ipynb | 168 +++++++++++++++++++ 1 file changed, 168 insertions(+) create mode 100644 docs/examples/18_particle_tracking_prt.ipynb diff --git a/docs/examples/18_particle_tracking_prt.ipynb b/docs/examples/18_particle_tracking_prt.ipynb new file mode 100644 index 00000000..39d6bae3 --- /dev/null +++ b/docs/examples/18_particle_tracking_prt.ipynb @@ -0,0 +1,168 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Particle Tracking with PRT" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from pathlib import Path\n", + "\n", + "import flopy as fp\n", + "import numpy as np\n", + "import xarray as xr\n", + "\n", + "import nlmod" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "nlmod.util.get_color_logger(\"INFO\")\n", + "nlmod.show_versions()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "model_ws = Path(\"./scratch_model\")\n", + "model_name = \"from_scratch\"\n", + "\n", + "# create new subdir for PRT model\n", + "prt_ws = Path(model_ws) / \"prt\"" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "ds = xr.open_dataset(model_ws / f\"{model_name}.nc\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Create a new Simulation object\n", + "simprt = nlmod.sim.sim(ds, sim_ws=prt_ws)\n", + "tdis = nlmod.sim.tdis(ds, simprt)\n", + "\n", + "# Add PRT model\n", + "prt = nlmod.prt.prt(ds, simprt)\n", + "\n", + "# DIS: discretization package\n", + "dis = nlmod.prt.dis(ds, prt)\n", + "\n", + "# MIP: matrix input package(?)\n", + "mip = nlmod.prt.mip(ds, prt, porosity=0.3)\n", + "\n", + "# PRP: particle release package\n", + "pd = fp.modpath.ParticleData(\n", + " [\n", + " (k, i, j)\n", + " for i in np.arange(0, ds.sizes[\"y\"], 11)\n", + " for j in np.arange(0, ds.sizes[\"x\"], 11)\n", + " for k in [0, 1, 2]\n", + " ],\n", + " structured=True,\n", + ")\n", + "releasepts = list(pd.to_prp(prt.modelgrid, global_xy=False))\n", + "prp = nlmod.prt.prp(ds, prt, packagedata=releasepts, perioddata={0: [\"FIRST\"]})\n", + "\n", + "# FMI: flow model interface\n", + "gwf_budget_file = Path(ds.model_ws) / f\"{model_name}.cbc\"\n", + "gwf_head_file = Path(ds.model_ws) / f\"{model_name}.hds\"\n", + "fmi = fp.mf6.ModflowPrtfmi(\n", + " prt,\n", + " packagedata=[\n", + " (\"GWFHEAD\", gwf_head_file.absolute().resolve()),\n", + " (\"GWFBUDGET\", gwf_budget_file.absolute().resolve()),\n", + " ],\n", + ")\n", + "\n", + "# OC: output control\n", + "oc = nlmod.prt.oc(ds, prt)\n", + "\n", + "# EMS: explicit model solution\n", + "ems = nlmod.sim.ems(simprt, model=prt)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "simprt.write_simulation()\n", + "simprt.run_simulation()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Load pathline data" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from flopy.utils.prtfile import PathlineFile\n", + "\n", + "trkfile = PathlineFile(\n", + " Path(prt.model_ws) / f\"{prt.name}.trk\",\n", + " header_filename=Path(prt.model_ws) / f\"{prt.name}.trk.hdr\",\n", + ")\n", + "trk = trkfile.get_alldata()\n", + "trk_df = trkfile.get_dataframe()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Plot results" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "pmv = fp.plot.PlotMapView(prt)\n", + "# loop to avoid different tracks get connected in plot\n", + "for i in range(len(trk)):\n", + " pmv.plot_pathline(trk[i:i+1], layer=\"all\", colors=\"k\", lw=1.0)\n", + "pmv.plot_endpoint(trk_df, color=\"r\", marker=\".\", direction=\"starting\", zorder=10);" + ] + } + ], + "metadata": { + "language_info": { + "name": "python" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} From f784fa97e1d93a30faf9716ac25b63dcf47e34ff Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dav=C3=ADd=20Brakenhoff?= Date: Tue, 24 Dec 2024 14:44:33 +0100 Subject: [PATCH 09/10] save all flows for prt example --- docs/examples/00_model_from_scratch.ipynb | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/docs/examples/00_model_from_scratch.ipynb b/docs/examples/00_model_from_scratch.ipynb index 81842350..0107e9d5 100644 --- a/docs/examples/00_model_from_scratch.ipynb +++ b/docs/examples/00_model_from_scratch.ipynb @@ -125,7 +125,9 @@ "ims = nlmod.sim.ims(sim, complexity=\"SIMPLE\")\n", "gwf = nlmod.gwf.gwf(ds, sim)\n", "dis = nlmod.gwf.dis(ds, gwf)\n", - "npf = nlmod.gwf.npf(ds, gwf)\n", + "npf = nlmod.gwf.npf(\n", + " ds, gwf, save_flows=True, save_specific_discharge=True, save_saturation=True\n", + ")\n", "ic = nlmod.gwf.ic(ds, gwf, starting_head=1.0)\n", "oc = nlmod.gwf.oc(ds, gwf, save_head=True)" ] From 73c7c4e2ba0b653046b3367f9695707bc8efb926 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dav=C3=ADd=20Brakenhoff?= Date: Thu, 26 Dec 2024 17:15:18 +0100 Subject: [PATCH 10/10] fix gwtdis fname, should be str not tuple --- nlmod/gwf/gwf.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nlmod/gwf/gwf.py b/nlmod/gwf/gwf.py index 241232f1..f09b2b30 100644 --- a/nlmod/gwf/gwf.py +++ b/nlmod/gwf/gwf.py @@ -117,7 +117,7 @@ def _dis(ds, model, length_units="METERS", pname="dis", **kwargs): filename = f"{ds.model_name}.dis" klass = flopy.mf6.ModflowGwfdis elif model.model_type == "gwt6": - filename = (f"{ds.model_name}_gwt.dis",) + filename = f"{ds.model_name}_gwt.dis" klass = flopy.mf6.ModflowGwtdis elif model.model_type == "prt6": filename = f"{ds.model_name}_prt.dis"