diff --git a/.gitignore b/.gitignore index fbc9b2b7..11b18283 100644 --- a/.gitignore +++ b/.gitignore @@ -1,5 +1,7 @@ build/ docs/examples/*.ipynb +docs/examples/quickstart/* +docs/examples/twri/* venv/ .DS_Store .ruff_cache diff --git a/docs/examples/advanced_packages.py b/docs/examples/advanced_packages.py deleted file mode 100644 index e69de29b..00000000 diff --git a/docs/examples/image/qs.png b/docs/examples/image/qs.png deleted file mode 100644 index 0d64e08c..00000000 Binary files a/docs/examples/image/qs.png and /dev/null differ diff --git a/docs/examples/image/qsprj.png b/docs/examples/image/qsprj.png deleted file mode 100644 index f71b239a..00000000 Binary files a/docs/examples/image/qsprj.png and /dev/null differ diff --git a/docs/examples/image/quickstart.png b/docs/examples/image/quickstart.png index 9e4daba5..d775d78e 100644 Binary files a/docs/examples/image/quickstart.png and b/docs/examples/image/quickstart.png differ diff --git a/docs/examples/qsplot.py b/docs/examples/qsplot.py deleted file mode 100644 index 02d82045..00000000 --- a/docs/examples/qsplot.py +++ /dev/null @@ -1,111 +0,0 @@ -import os - -import cartopy.crs as ccrs -import flopy -import matplotlib.pyplot as plt -import numpy as np -import xarray as xr -from matplotlib.collections import PolyCollection - -QS_ROOT = os.path.abspath(os.path.join(__file__, "..")) -modelname = "mymodel" -ws = os.path.join(QS_ROOT, "quickstart_data") - - -def _pcolormesh(grid, da, ax): - vertices = [] - for face in grid.iverts: - f = [] - for i in face: - f.append(grid.verts[i]) - vertices.append(f) - vertices = np.array(vertices) - - collection = PolyCollection(vertices) - collection.set_array(da.to_numpy().ravel()) - p = ax.add_collection(collection, autolim=False) - - xmin, xmax, ymin, ymax = grid.extent - ax.set_xlim(xmin, xmax) - ax.set_ylim(ymin, ymax) - cverts = (xmin, ymin), (xmax, ymax) - ax.update_datalim(cverts) - - return p - - -# create budget reader -bpth = os.path.join(ws, f"{modelname}.bud") -bobj = flopy.utils.CellBudgetFile(bpth, precision="double") - -# set specific discharge -spdis = bobj.get_data(text="DATA-SPDIS")[0] - -# create head reader -hpth = os.path.join(ws, f"{modelname}.hds") -hobj = flopy.utils.HeadFile(hpth, precision="double") - -# set heads -heads = hobj.get_alldata() - -# create grb reader -grbpth = os.path.join(ws, f"{modelname}.dis.grb") -grbobj = flopy.mf6.utils.MfGrdFile(grbpth) -grid = None - -# flow vector component arrays -uflow = None -vflow = None - -# create grid object -if "NROW" in grbobj._datadict and "NCOL" in grbobj._datadict: - grid = flopy.discretization.StructuredGrid.from_binary_grid_file(grbpth) - u = [] - v = [] - for r in spdis: - u.append(r[3]) - v.append(r[4]) - uflow = np.array(u).reshape(grid.nrow, grid.ncol) - vflow = np.array(v).reshape(grid.nrow, grid.ncol) -else: - raise Exception("unsupported discretization") - -# set data coordinate arrays -xcrs = grid.xycenters[0] -ycrs = grid.xycenters[1] - -# create qx, qy dataarrys -qx = xr.DataArray(uflow, dims=("y", "x"), coords={"x": xcrs, "y": ycrs}) -qy = xr.DataArray(vflow, dims=("y", "x"), coords={"x": xcrs, "y": ycrs}) - -# create head data array -da = xr.DataArray( - heads[0][0], - dims=("y", "x"), - coords={"x": xcrs, "y": ycrs}, - name="head", -) - -# quickstart mesh with contours and flow arrows -fig, ax = plt.subplots() -_pcolormesh(grid, da, ax) -qz = np.sin(np.sqrt(uflow**2 + vflow**2)) -plt.quiver(xcrs, ycrs, qx, qy, color="w") -cbar = plt.colorbar(label=da.name) -plt.contour(xcrs, ycrs, qz, 4, cmap="viridis") -# plt.contourf(xcrs, ycrs, qz, 5, cmap='viridis', alpha=0.4) -qs_pth = os.path.join(QS_ROOT, "image", "qs.png") -fig.savefig(qs_pth) - -# projection example with cartopy -fig, ax = plt.subplots(subplot_kw={"projection": ccrs.PlateCarree()}) -da.plot(ax=ax, transform=ccrs.PlateCarree()) -ax.coastlines() -ax.stock_img() -ax.set_extent([-5, 15, -4, 14], crs=ccrs.PlateCarree()) -gln = ax.gridlines(draw_labels=True) -gln.top_labels = False -gln.right_labels = False -ax.set_title("Head") -qsprj_pth = os.path.join(QS_ROOT, "image", "qsprj.png") -fig.savefig(qsprj_pth) diff --git a/docs/examples/quickstart.py b/docs/examples/quickstart.py index ad2e5228..3e64c1b4 100644 --- a/docs/examples/quickstart.py +++ b/docs/examples/quickstart.py @@ -11,6 +11,7 @@ name = "quickstart" workspace = Path(__file__).parent / name +workspace.mkdir(exist_ok=True) time = Time(perlen=[1.0], nstp=[1]) grid = StructuredGrid( nlay=1, @@ -25,7 +26,7 @@ gwf_name = "mymodel" ims = Ims(parent=sim, models=[gwf_name]) # temporary hack gwf = Gwf(parent=sim, name=gwf_name, save_flows=True, dis=grid) -npf = Npf(parent=gwf, save_specific_discharge=True) +npf = Npf(parent=gwf, print_flows=True, save_flows=True, save_specific_discharge=True) chd = Chd( parent=gwf, head={0: {(0, 0, 0): 1.0, (0, 9, 9): 0.0}}, diff --git a/docs/examples/quickstart/mfsim.lst b/docs/examples/quickstart/mfsim.lst deleted file mode 100644 index 5a80e3b9..00000000 --- a/docs/examples/quickstart/mfsim.lst +++ /dev/null @@ -1,248 +0,0 @@ - MODFLOW 6 - U.S. GEOLOGICAL SURVEY MODULAR HYDROLOGIC MODEL - VERSION 6.6.1 02/10/2025 - - MODFLOW 6 compiled Feb 14 2025 13:38:43 with Intel(R) Fortran Intel(R) 64 - Compiler Classic for applications running on Intel(R) 64, Version 2021.7.0 - Build 20220726_000000 - -This software has been approved for release by the U.S. Geological -Survey (USGS). Although the software has been subjected to rigorous -review, the USGS reserves the right to update the software as needed -pursuant to further analysis and review. No warranty, expressed or -implied, is made by the USGS or the U.S. Government as to the -functionality of the software and related material nor shall the -fact of release constitute any such warranty. Furthermore, the -software is released on condition that neither the USGS nor the U.S. -Government shall be held liable for any damages resulting from its -authorized or unauthorized use. Also refer to the USGS Water -Resources Software User Rights Notice for complete use, copyright, -and distribution information. - - -As a work of the United States Government, this USGS product is -in the public domain within the United States. You can copy, -modify, distribute, and perform the work, even for commercial -purposes, all without asking permission. Additionally, USGS -waives copyright and related rights in the work worldwide -through CC0 1.0 Universal Public Domain Dedication -(https://creativecommons.org/publicdomain/zero/1.0/). - -The following GNU Lesser General Public License (LGPL) libraries -are used in this USGS product: - - SPARSKIT version 2.0 - ilut, luson, and qsplit - (https://www-users.cse.umn.edu/~saad/software/SPARSKIT/) - - RCM - Reverse Cuthill McKee Ordering - (https://people.math.sc.edu/Burkardt/f_src/rcm/rcm.html) - - BLAS - Basic Linear Algebra Subprograms Level 1 - (https://people.math.sc.edu/Burkardt/f_src/blas1_d/blas1_d.html) - - SPARSEKIT - Sparse Matrix Utility Package - amux, dperm, dvperm, rperm, and cperm - (https://people.sc.fsu.edu/~jburkardt/f77_src/sparsekit/sparsekit.html) - -The following BSD-3 License libraries are used in this USGS product: - - Modern Fortran DAG Library - Copyright (c) 2018, Jacob Williams - All rights reserved. - (https://github.com/jacobwilliams/daglib) - -MODFLOW 6 compiler options: -Imac/obj_mf6 -O2 -no-heap-arrays -fpe0 --traceback -Qdiag-disable:7416 -Qdiag-disable:7025 -Qdiag-disable:5268 -fpp --module mac/mod_mf6/ -c -o mac/obj_mf6/compilerversion.o - -System command used to initiate simulation: -/Users/wbonelli/dev/flopy/.venv/bin/mf6 - -MODFLOW was compiled using uniform precision. - -Real Variables - KIND: 8 - TINY (smallest non-zero value): 2.225074-308 - HUGE (largest value): 1.797693+308 - PRECISION: 15 - SIZE IN BITS: 64 - -Integer Variables - KIND: 4 - HUGE (largest value): 2147483647 - SIZE IN BITS: 32 - -Long Integer Variables - KIND: 8 - HUGE (largest value): 9223372036854775807 - SIZE IN BITS: 64 - -Logical Variables - KIND: 4 - SIZE IN BITS: 32 - - - OPENED mfsim.nam - FILE TYPE:NAM6 UNIT 1001 STATUS:OLD - FORMAT:FORMATTED ACCESS:SEQUENTIAL - ACTION:READ - - # File generated by Flopy version 3.10.0.dev2 on 04/04/2025 at 09:49:59. - - OPENED mymodel.tdis - FILE TYPE:TDIS6 UNIT 1002 STATUS:OLD - FORMAT:FORMATTED ACCESS:SEQUENTIAL - ACTION:READ - - # File generated by Flopy version 3.10.0.dev2 on 04/04/2025 at 09:49:59. - - OPENED mymodel.nam - FILE TYPE:GWF6 UNIT 1003 STATUS:OLD - FORMAT:FORMATTED ACCESS:SEQUENTIAL - ACTION:READ - - # File generated by Flopy version 3.10.0.dev2 on 04/04/2025 at 09:49:59. - - OPENED mymodel.dis - FILE TYPE:DIS6 UNIT 1004 STATUS:OLD - FORMAT:FORMATTED ACCESS:SEQUENTIAL - ACTION:READ - - # File generated by Flopy version 3.10.0.dev2 on 04/04/2025 at 09:49:59. - - OPENED mymodel.npf - FILE TYPE:NPF6 UNIT 1005 STATUS:OLD - FORMAT:FORMATTED ACCESS:SEQUENTIAL - ACTION:READ - - # File generated by Flopy version 3.10.0.dev2 on 04/04/2025 at 09:49:59. - - OPENED mymodel.ic - FILE TYPE:IC6 UNIT 1006 STATUS:OLD - FORMAT:FORMATTED ACCESS:SEQUENTIAL - ACTION:READ - - # File generated by Flopy version 3.10.0.dev2 on 04/04/2025 at 09:49:59. - - OPENED mymodel.oc - FILE TYPE:OC6 UNIT 1007 STATUS:OLD - FORMAT:FORMATTED ACCESS:SEQUENTIAL - ACTION:READ - - - OPENED mymodel.chd - FILE TYPE:CHD6 UNIT 1008 STATUS:OLD - FORMAT:FORMATTED ACCESS:SEQUENTIAL - ACTION:READ - - # File generated by Flopy version 3.10.0.dev2 on 04/04/2025 at 09:49:59. - - READING SIMULATION OPTIONS - MAXIMUM NUMBER OF ERRORS THAT WILL BE STORED IS 1000 - END OF SIMULATION OPTIONS - - READING SIMULATION TIMING - - TDIS -- TEMPORAL DISCRETIZATION PACKAGE, - VERSION 1 : 11/13/2014 - INPUT READ FROM MEMPATH: __INPUT__/SIM/TDIS - PROCESSING TDIS OPTIONS - SIMULATION TIME UNIT IS UNDEFINED - END OF TDIS OPTIONS - PROCESSING TDIS DIMENSIONS - 1 STRESS PERIOD(S) IN SIMULATION - END OF TDIS DIMENSIONS - PROCESSING TDIS PERIODDATA - - - STRESS PERIOD LENGTH TIME STEPS MULTIPLIER FOR DELT - ---------------------------------------------------------------------------- - 1 1.000000 1 1.000 - END OF TDIS PERIODDATA - END OF SIMULATION TIMING - - READING SIMULATION MODELS - GWF6 model 1 will be created - END OF SIMULATION MODELS - - READING SIMULATION EXCHANGES - END OF SIMULATION EXCHANGES - - READING SOLUTIONGROUP - - Creating solution: SLN_1 - - OPENED mymodel.ims - FILE TYPE:IMS UNIT 1010 STATUS:OLD - FORMAT:FORMATTED ACCESS:SEQUENTIAL - ACTION:READ - - END OF SOLUTIONGROUP - -PROCESSING MODEL CONNECTIONS - - IMS -- ITERATIVE MODEL SOLUTION PACKAGE, VERSION 6, 4/28/2017 - INPUT READ FROM UNIT 1010 - # File generated by Flopy version 3.10.0.dev2 on 04/04/2025 at 09:49:59. - - PROCESSING IMS OPTIONS - END OF IMS OPTIONS - ***UNDER-RELAXATION WILL NOT BE USED*** - - ***IMS LINEAR SOLVER WILL BE USED*** - - IMSLINEAR -- UNSTRUCTURED LINEAR SOLUTION PACKAGE, VERSION 8, 04/28/2017 - A symmetric matrix will be solved - - OUTER ITERATION CONVERGENCE CRITERION (DVCLOSE) = 0.100000E-02 - MAXIMUM NUMBER OF OUTER ITERATIONS (MXITER) = 25 - SOLVER PRINTOUT INDEX (IPRIMS) = 0 - NONLINEAR ITERATION METHOD (NONLINMETH) = 0 - LINEAR SOLUTION METHOD (LINMETH) = 1 - - SOLUTION BY THE CONJUGATE-GRADIENT METHOD - ------------------------------------------------------------------ - MAXIMUM OF 25 CALLS OF SOLUTION ROUTINE - MAXIMUM OF 50 INTERNAL ITERATIONS PER CALL TO SOLUTION ROUTINE - LINEAR ACCELERATION METHOD = CG - MATRIX PRECONDITIONING TYPE = INCOMPLETE LU - MATRIX SCALING APPROACH = NO SCALING - MATRIX REORDERING APPROACH = ORIGINAL ORDERING - NUMBER OF ORTHOGONALIZATIONS = 0 - HEAD CHANGE CRITERION FOR CLOSURE = 0.10000E-02 - RESIDUAL CHANGE CRITERION FOR CLOSURE = 0.10000E+00 - RESIDUAL CONVERGENCE OPTION = 0 - RESIDUAL CONVERGENCE NORM = INFINITY NORM - RELAXATION FACTOR = 0.00000E+00 - - - - - Solving: Stress period: 1 Time step: 1 - -1 - STRESS PERIOD NO. 1, LENGTH = 1.000000 - ------------------------------------------ - NUMBER OF TIME STEPS = 1 - MULTIPLIER FOR DELT = 1.000 - INITIAL TIME STEP SIZE = 1.000000 - - MEMORY MANAGER TOTAL STORAGE BY DATA TYPE, IN KILOBYTES - ------------------------------- - ALLOCATED - DATA TYPE MEMORY - ------------------------------- - Character 6.8550000 - Logical 4.80000000E-02 - Integer 18.988000 - Real 43.160000 - ------------------------------- - Total 69.051000 - Virtual 0.0000000 - ------------------------------- - - - - Run end date and time (yyyy/mm/dd hh:mm:ss): 2025/04/04 9:50:03 - Elapsed run time: 0.057 Seconds - Normal termination of simulation. diff --git a/docs/examples/quickstart/mfsim.nam b/docs/examples/quickstart/mfsim.nam deleted file mode 100644 index ffbd1173..00000000 --- a/docs/examples/quickstart/mfsim.nam +++ /dev/null @@ -1,19 +0,0 @@ -# File generated by Flopy version 3.10.0.dev2 on 04/04/2025 at 09:49:59. -BEGIN options -END options - -BEGIN timing - TDIS6 mymodel.tdis -END timing - -BEGIN models - gwf6 mymodel.nam mymodel -END models - -BEGIN exchanges -END exchanges - -BEGIN solutiongroup 1 - ims6 mymodel.ims mymodel -END solutiongroup 1 - diff --git a/docs/examples/quickstart/mymodel.bud b/docs/examples/quickstart/mymodel.bud deleted file mode 100644 index 0938f1d2..00000000 Binary files a/docs/examples/quickstart/mymodel.bud and /dev/null differ diff --git a/docs/examples/quickstart/mymodel.chd b/docs/examples/quickstart/mymodel.chd deleted file mode 100644 index 6710879b..00000000 --- a/docs/examples/quickstart/mymodel.chd +++ /dev/null @@ -1,13 +0,0 @@ -# File generated by Flopy version 3.10.0.dev2 on 04/04/2025 at 09:49:59. -BEGIN options -END options - -BEGIN dimensions - MAXBOUND 2 -END dimensions - -BEGIN period 1 - 1 1 1 1.00000000E+00 - 1 10 10 0.00000000E+00 -END period 1 - diff --git a/docs/examples/quickstart/mymodel.dis b/docs/examples/quickstart/mymodel.dis deleted file mode 100644 index 9fea5c3e..00000000 --- a/docs/examples/quickstart/mymodel.dis +++ /dev/null @@ -1,21 +0,0 @@ -# File generated by Flopy version 3.10.0.dev2 on 04/04/2025 at 09:49:59. -BEGIN options -END options - -BEGIN dimensions - NLAY 1 - NROW 10 - NCOL 10 -END dimensions - -BEGIN griddata - delr - CONSTANT 1.00000000 - delc - CONSTANT 1.00000000 - top - CONSTANT 1.00000000 - botm - CONSTANT 0.00000000 -END griddata - diff --git a/docs/examples/quickstart/mymodel.dis.grb b/docs/examples/quickstart/mymodel.dis.grb deleted file mode 100644 index e8373212..00000000 Binary files a/docs/examples/quickstart/mymodel.dis.grb and /dev/null differ diff --git a/docs/examples/quickstart/mymodel.hds b/docs/examples/quickstart/mymodel.hds deleted file mode 100644 index 29824fc5..00000000 Binary files a/docs/examples/quickstart/mymodel.hds and /dev/null differ diff --git a/docs/examples/quickstart/mymodel.ic b/docs/examples/quickstart/mymodel.ic deleted file mode 100644 index 34fff374..00000000 --- a/docs/examples/quickstart/mymodel.ic +++ /dev/null @@ -1,9 +0,0 @@ -# File generated by Flopy version 3.10.0.dev2 on 04/04/2025 at 09:49:59. -BEGIN options -END options - -BEGIN griddata - strt - CONSTANT 1.00000000 -END griddata - diff --git a/docs/examples/quickstart/mymodel.ims b/docs/examples/quickstart/mymodel.ims deleted file mode 100644 index 64506566..00000000 --- a/docs/examples/quickstart/mymodel.ims +++ /dev/null @@ -1,4 +0,0 @@ -# File generated by Flopy version 3.10.0.dev2 on 04/04/2025 at 09:49:59. -BEGIN options -END options - diff --git a/docs/examples/quickstart/mymodel.lst b/docs/examples/quickstart/mymodel.lst deleted file mode 100644 index 652e8be7..00000000 --- a/docs/examples/quickstart/mymodel.lst +++ /dev/null @@ -1,217 +0,0 @@ - MODFLOW 6 - U.S. GEOLOGICAL SURVEY MODULAR HYDROLOGIC MODEL - GROUNDWATER FLOW MODEL (GWF) - VERSION 6.6.1 02/10/2025 - - MODFLOW 6 compiled Feb 14 2025 13:38:43 with Intel(R) Fortran Intel(R) 64 - Compiler Classic for applications running on Intel(R) 64, Version 2021.7.0 - Build 20220726_000000 - -This software has been approved for release by the U.S. Geological -Survey (USGS). Although the software has been subjected to rigorous -review, the USGS reserves the right to update the software as needed -pursuant to further analysis and review. No warranty, expressed or -implied, is made by the USGS or the U.S. Government as to the -functionality of the software and related material nor shall the -fact of release constitute any such warranty. Furthermore, the -software is released on condition that neither the USGS nor the U.S. -Government shall be held liable for any damages resulting from its -authorized or unauthorized use. Also refer to the USGS Water -Resources Software User Rights Notice for complete use, copyright, -and distribution information. - - -As a work of the United States Government, this USGS product is -in the public domain within the United States. You can copy, -modify, distribute, and perform the work, even for commercial -purposes, all without asking permission. Additionally, USGS -waives copyright and related rights in the work worldwide -through CC0 1.0 Universal Public Domain Dedication -(https://creativecommons.org/publicdomain/zero/1.0/). - -The following GNU Lesser General Public License (LGPL) libraries -are used in this USGS product: - - SPARSKIT version 2.0 - ilut, luson, and qsplit - (https://www-users.cse.umn.edu/~saad/software/SPARSKIT/) - - RCM - Reverse Cuthill McKee Ordering - (https://people.math.sc.edu/Burkardt/f_src/rcm/rcm.html) - - BLAS - Basic Linear Algebra Subprograms Level 1 - (https://people.math.sc.edu/Burkardt/f_src/blas1_d/blas1_d.html) - - SPARSEKIT - Sparse Matrix Utility Package - amux, dperm, dvperm, rperm, and cperm - (https://people.sc.fsu.edu/~jburkardt/f77_src/sparsekit/sparsekit.html) - -The following BSD-3 License libraries are used in this USGS product: - - Modern Fortran DAG Library - Copyright (c) 2018, Jacob Williams - All rights reserved. - (https://github.com/jacobwilliams/daglib) - -MODFLOW 6 compiler options: -Imac/obj_mf6 -O2 -no-heap-arrays -fpe0 --traceback -Qdiag-disable:7416 -Qdiag-disable:7025 -Qdiag-disable:5268 -fpp --module mac/mod_mf6/ -c -o mac/obj_mf6/compilerversion.o - -System command used to initiate simulation: -/Users/wbonelli/dev/flopy/.venv/bin/mf6 - -MODFLOW was compiled using uniform precision. - -Real Variables - KIND: 8 - TINY (smallest non-zero value): 2.225074-308 - HUGE (largest value): 1.797693+308 - PRECISION: 15 - SIZE IN BITS: 64 - -Integer Variables - KIND: 4 - HUGE (largest value): 2147483647 - SIZE IN BITS: 32 - -Long Integer Variables - KIND: 8 - HUGE (largest value): 9223372036854775807 - SIZE IN BITS: 64 - -Logical Variables - KIND: 4 - SIZE IN BITS: 32 - - NAMEFILE OPTIONS: - FLOWS WILL BE SAVED TO BUDGET FILE SPECIFIED IN OUTPUT CONTROL - END NAMEFILE OPTIONS: - - DIS -- STRUCTURED GRID DISCRETIZATION PACKAGE, VERSION 2 : 3/27/2014 - INPUT READ FROM MEMPATH: __INPUT__/MYMODEL/DIS - - - NPF -- NODE PROPERTY FLOW PACKAGE, VERSION 1, 3/30/2015 INPUT READ FROM MEMPATH: __INPUT__/MYMODEL/NPF - - - IC -- Initial Conditions Package, Version 8, 3/28/2015 input read from mempath: __INPUT__/MYMODEL/IC - - - Setting Discretization Options - End Setting Discretization Options - - Setting Discretization Dimensions - NLAY = 1 - NROW = 10 - NCOL = 10 - End Setting Discretization Dimensions - - Setting Discretization Griddata - DELR set from input file - DELC set from input file - TOP set from input file - BOTM set from input file - End Setting Discretization Griddata - - Setting NPF Options - Specific discharge will be calculated at cell centers and written to DATA-SPDIS in budget file when requested. - End Setting NPF Options - - Setting NPF Griddata - ICELLTYPE set from input file - K set from input file - K33 not provided. Setting K33 = K. - K22 not provided. Setting K22 = K. - End Setting NPF Griddata - - - CHD -- CHD PACKAGE, VERSION 8, 2/22/2014 INPUT READ FROM MEMPATH: __INPUT__/MYMODEL/CHD_0 - - PROCESSING CHD BASE OPTIONS - END OF CHD BASE OPTIONS - - PROCESSING CHD BASE DIMENSIONS - MAXBOUND = 2 - END OF CHD BASE DIMENSIONS - STRT set from input file - BINARY GRID INFORMATION WILL BE WRITTEN TO: - UNIT NUMBER: 1011 - FILE NAME: mymodel.dis.grb - - OPENED mymodel.dis.grb - FILE TYPE:DATA(BINARY) UNIT 1011 STATUS:REPLACE - FORMAT:UNFORMATTED ACCESS:STREAM - ACTION:READWRITE - - THE LAST TIME STEP WILL BE PRINTED - THE LAST TIME STEP WILL BE PRINTED - # File generated by Flopy version 3.10.0.dev2 on 04/04/2025 at 09:49:59. - - PROCESSING OC OPTIONS - - BUDGET INFORMATION WILL BE WRITTEN TO: - UNIT NUMBER: 1012 - FILE NAME: mymodel.bud - - OPENED mymodel.bud - FILE TYPE:DATA(BINARY) UNIT 1012 STATUS:REPLACE - FORMAT:UNFORMATTED ACCESS:STREAM - ACTION:READWRITE - - HEAD INFORMATION WILL BE WRITTEN TO: - UNIT NUMBER: 1013 - FILE NAME: mymodel.hds - - OPENED mymodel.hds - FILE TYPE:DATA(BINARY) UNIT 1013 STATUS:REPLACE - FORMAT:UNFORMATTED ACCESS:STREAM - ACTION:READWRITE - - END OF OC OPTIONS - -start timestep kper="1" kstp="1" mode="normal" - - - BEGIN READING OUTPUT CONTROL FOR STRESS PERIOD 1 - ALL TIME STEPS WILL BE SAVED - ALL TIME STEPS WILL BE SAVED - - END READING OUTPUT CONTROL FOR STRESS PERIOD 1 - UBDSV1 SAVING FLOW-JA-FACE ON UNIT 1012 AT TIME STEP 1, STRESS PERIOD 1 - UBDSV06 SAVING DATA-SPDIS IN MODEL MYMODEL PACKAGE NPF CONNECTED TO MODEL MYMODEL PACKAGE NPF ON UNIT 1012 AT TIME STEP 1, STRESS PERIOD 1 - UBDSV06 SAVING CHD IN MODEL MYMODEL PACKAGE MYMODEL CONNECTED TO MODEL MYMODEL PACKAGE CHD_0 ON UNIT 1012 AT TIME STEP 1, STRESS PERIOD 1 - - HEAD WILL BE SAVED ON UNIT 1013 AT END OF TIME STEP 1, STRESS PERIOD 1 - - - VOLUME BUDGET FOR ENTIRE MODEL AT END OF TIME STEP 1, STRESS PERIOD 1 - --------------------------------------------------------------------------------------------------- - - CUMULATIVE VOLUME L**3 RATES FOR THIS TIME STEP L**3/T PACKAGE NAME - ------------------ ------------------------ ---------------- - - IN: IN: - --- --- - CHD = 0.3321 CHD = 0.3321 CHD_0 - - TOTAL IN = 0.3321 TOTAL IN = 0.3321 - - OUT: OUT: - ---- ---- - CHD = 0.3320 CHD = 0.3320 CHD_0 - - TOTAL OUT = 0.3320 TOTAL OUT = 0.3320 - - IN - OUT = 1.2232E-05 IN - OUT = 1.2232E-05 - - PERCENT DISCREPANCY = 0.00 PERCENT DISCREPANCY = 0.00 - - - - - TIME SUMMARY AT END OF TIME STEP 1 IN STRESS PERIOD 1 - TIME STEP LENGTH = 1.00000 - STRESS PERIOD TIME = 1.00000 - TOTAL SIMULATION TIME = 1.00000 - -end timestep - diff --git a/docs/examples/quickstart/mymodel.nam b/docs/examples/quickstart/mymodel.nam deleted file mode 100644 index f9c98ed6..00000000 --- a/docs/examples/quickstart/mymodel.nam +++ /dev/null @@ -1,13 +0,0 @@ -# File generated by Flopy version 3.10.0.dev2 on 04/04/2025 at 09:49:59. -BEGIN options - SAVE_FLOWS -END options - -BEGIN packages - DIS6 mymodel.dis dis - IC6 mymodel.ic ic - NPF6 mymodel.npf npf - CHD6 mymodel.chd chd_0 - OC6 mymodel.oc oc -END packages - diff --git a/docs/examples/quickstart/mymodel.npf b/docs/examples/quickstart/mymodel.npf deleted file mode 100644 index 208c46b1..00000000 --- a/docs/examples/quickstart/mymodel.npf +++ /dev/null @@ -1,12 +0,0 @@ -# File generated by Flopy version 3.10.0.dev2 on 04/04/2025 at 09:49:59. -BEGIN options - SAVE_SPECIFIC_DISCHARGE -END options - -BEGIN griddata - icelltype - CONSTANT 0 - k - CONSTANT 1.00000000 -END griddata - diff --git a/docs/examples/quickstart/mymodel.oc b/docs/examples/quickstart/mymodel.oc deleted file mode 100644 index fda2173d..00000000 --- a/docs/examples/quickstart/mymodel.oc +++ /dev/null @@ -1,11 +0,0 @@ -# File generated by Flopy version 3.10.0.dev2 on 04/04/2025 at 09:49:59. -BEGIN options - BUDGET FILEOUT mymodel.bud - HEAD FILEOUT mymodel.hds -END options - -BEGIN period 1 - SAVE HEAD ALL - SAVE BUDGET ALL -END period 1 - diff --git a/docs/examples/quickstart/mymodel.tdis b/docs/examples/quickstart/mymodel.tdis deleted file mode 100644 index 810a7b01..00000000 --- a/docs/examples/quickstart/mymodel.tdis +++ /dev/null @@ -1,12 +0,0 @@ -# File generated by Flopy version 3.10.0.dev2 on 04/04/2025 at 09:49:59. -BEGIN options -END options - -BEGIN dimensions - NPER 1 -END dimensions - -BEGIN perioddata - 1.00000000 1 1.00000000 -END perioddata - diff --git a/docs/examples/quickstart/quickstart.png b/docs/examples/quickstart/quickstart.png deleted file mode 100644 index 9e4daba5..00000000 Binary files a/docs/examples/quickstart/quickstart.png and /dev/null differ diff --git a/docs/examples/quickstart_expanded.py b/docs/examples/quickstart_expanded.py index 1d53e24f..2149446e 100644 --- a/docs/examples/quickstart_expanded.py +++ b/docs/examples/quickstart_expanded.py @@ -31,16 +31,16 @@ ws = "./mymodel" name = "mymodel" -sim = Simulation(name=name, workspace=ws, exe="mf6") -tdis = Tdis(sim) -gwf = Gwf(sim, name=name, save_flows=True) -dis = Dis(gwf, nrow=10, ncol=10) -ic = Ic(gwf) -npf = Npf(gwf, save_specific_discharge=True) +sim = Simulation(name=name, workspace=ws) +tdis = Tdis(parent=sim) +gwf = Gwf(parent=sim, name=name, save_flows=True) +dis = Dis(parent=gwf, nrow=10, ncol=10) +ic = Ic(parent=gwf) +npf = Npf(parent=gwf, save_specific_discharge=True) # CHD. first, nothing but builtins. chd = Chd( - gwf, + parent=gwf, # 1) tuples (like flopy3) stress_period_data=[[(0, 0, 0), 1.0], [(0, 9, 9), 0.0]], # @@ -90,7 +90,7 @@ budget_file = name + ".bud" head_file = name + ".hds" oc = Oc( - gwf, + parent=gwf, budget_filerecord=budget_file, head_filerecord=head_file, # 1) tuples (like flopy3) diff --git a/docs/examples/stress_packages.py b/docs/examples/stress_packages.py deleted file mode 100644 index e69de29b..00000000 diff --git a/docs/examples/twri.py b/docs/examples/twri.py new file mode 100644 index 00000000..3f743e86 --- /dev/null +++ b/docs/examples/twri.py @@ -0,0 +1,184 @@ +# # TWRI +# +# Import dependencies. + +from pathlib import Path + +import numpy as np + +import flopy4 + +# Timing +time = flopy4.mf6.utils.time.Time.from_timestamps( + ["2000-01-01", "2000-01-02", "2000-01-03", "2000-01-04"] +) +nper = time.nper + +# Grid +nlay = 3 +nrow = 15 +ncol = 15 +shape = (nlay, nrow, ncol) +nodes = np.prod(shape) +delr = np.ones((ncol,)) * 5000.0 +delc = np.ones((nrow,)) * 5000.0 +idomain = np.ones(shape, dtype=int) +top = np.ones((nrow, ncol), dtype=float) * 200.0 +bottom = np.stack([np.full((nrow, ncol), val) for val in [-200.0, -300.0, -450.0]]) +grid = flopy4.mf6.utils.grid.StructuredGrid( + nlay=nlay, nrow=nrow, ncol=ncol, top=top, botm=bottom, delr=delr, delc=delc, idomain=idomain +) +dims = {"nper": nper, **dict(grid.dataset.sizes)} # TODO: temporary + +# Grid discretization +# TODO: xorigin, yorigin +dis = flopy4.mf6.gwf.Dis.from_grid(grid=grid) + +# Constant head boundary on the left +chd = flopy4.mf6.gwf.Chd( + head={"*": {(k, i, 0): 0.0 for k in range(nlay - 1) for i in range(nrow)}}, + print_input=True, + print_flows=True, + save_flows=True, + dims=dims, +) + +# Drain in the center left of the model +elevation = [0.0, 0.0, 10.0, 20.0, 30.0, 50.0, 70.0, 90.0, 100.0] +conductance = 1.0 +drn = flopy4.mf6.gwf.Drn( + elev={"*": {(0, 7, j + 1): elevation[j] for j in range(9)}}, + cond={"*": {(0, 7, j + 1): conductance for j in range(9)}}, + print_input=True, + print_flows=True, + save_flows=True, + dims=dims, +) + +# Initial conditions +ic = flopy4.mf6.gwf.Ic(strt=0.0, dims=dims) + +# Node properties +icelltype = np.stack([np.full((nrow, ncol), val) for val in [1, 0, 0]]) +k = np.stack([np.full((nrow, ncol), val) for val in [1.0e-3, 1.0e-4, 2.0e-4]]) +k33 = np.stack([np.full((nrow, ncol), val) for val in [2.0e-8, 2.0e-8, 2.0e-8]]) +npf = flopy4.mf6.gwf.Npf( + # TODO: no need for reshaping once array structuring converter is done + icelltype=icelltype.reshape((nodes,)), + k=k.reshape((nodes,)), + k33=k33.reshape((nodes,)), + cvoptions=flopy4.mf6.gwf.Npf.CvOptions(dewatered=True), + perched=True, + save_flows=True, + dims=dims, +) + +# Storage +sto = flopy4.mf6.gwf.Sto( + storagecoefficient=False, + ss=1.0e-5, + sy=0.15, + steady_state=[True, False, False], + iconvert=0, + dims=dims, +) + +# Uniform recharge on the top layer +rch_rate = np.stack(np.full((nlay, nrow, ncol), flopy4.mf6.constants.FILL_DNODATA)) +rate = np.repeat(np.expand_dims(rch_rate, axis=0), repeats=nper, axis=0) +rate[0, 0, :] = 3.0e-8 +rch = flopy4.mf6.gwf.Rch(recharge=rate.reshape(nper, -1), dims=dims) + +# Output control +# TODO: show both ways to set up the Oc package, strings +# and proper record types? and/or with perioddata param? +oc = flopy4.mf6.gwf.Oc( + budget_file="gwf.bud", + head_file="gwf.hds", + save_head={0: "all"}, + save_budget={0: "all"}, + dims=dims, +) + +# Wells scattered throughout the model +wel_q = -5.0 +wel_nodes = [ + [2, 4, 10, -5.0], + [1, 3, 5, -5.0], + [1, 5, 11, -5.0], + [0, 8, 7, -5.0], + [0, 8, 9, -5.0], + [0, 8, 11, -5.0], + [0, 8, 13, -5.0], + [0, 10, 7, -5.0], + [0, 10, 9, -5.0], + [0, 10, 11, -5.0], + [0, 10, 13, -5.0], + [0, 12, 7, -5.0], + [0, 12, 9, -5.0], + [0, 12, 11, -5.0], + [0, 12, 13, -5.0], +] +wel = flopy4.mf6.gwf.Wel( + q={"*": {(layer, row, col): wel_q for layer, row, col, wel_q in wel_nodes}}, + dims=dims, +) + +# Flow model +gwf = flopy4.mf6.gwf.Gwf( + dis=grid, + ic=ic, + npf=npf, + sto=sto, + oc=oc, + chd=[chd], + rch=[rch], + wel=[wel], + drn=[drn], + dims=dims, +) + +# Solver +ims = flopy4.mf6.Ims( + print_option="summary", + outer_dvclose=1.0e-4, + outer_maximum=500, + under_relaxation=None, + inner_dvclose=1.0e-4, + inner_rclose=0.001, + inner_maximum=100, + linear_acceleration="cg", + scaling_method=None, + reordering_method=None, + relaxation_factor=0.97, + models=["gwf"], +) + +# TDIS +tdis = flopy4.mf6.simulation.Tdis.from_time(time) + +# Create workspace +workspace = Path(__file__).parent / "twri" +workspace.mkdir(parents=True, exist_ok=True) + +# Create simulation +sim = flopy4.mf6.simulation.Simulation( + name="twri", + tdis=tdis, + models={"gwf": gwf}, + solutions={"ims": ims}, + workspace=workspace, +) + +# Write input files and run the simulation +sim.write() +sim.run() # assumes the ``mf6`` executable is available on your PATH. + +# Load head results +head = flopy4.mf6.utils.open_hds( + workspace / f"{gwf.name}.hds", + workspace / f"{gwf.name}.dis.grb", +) + +# Plot head results +head.isel(layer=0, time=0).plot.contourf() diff --git a/flopy4/__init__.py b/flopy4/__init__.py index e69de29b..b7fbce76 100644 --- a/flopy4/__init__.py +++ b/flopy4/__init__.py @@ -0,0 +1,3 @@ +from flopy4 import mf6 + +__all__ = ["mf6"] diff --git a/flopy4/mf6/__init__.py b/flopy4/mf6/__init__.py index e92f50cc..31ce43a9 100644 --- a/flopy4/mf6/__init__.py +++ b/flopy4/mf6/__init__.py @@ -5,12 +5,25 @@ from tomli import load as load_toml from tomli_w import dump as dump_toml +# Import submodules to make them accessible via flopy4.mf6.* +from flopy4.mf6 import gwf, simulation, solution, utils from flopy4.mf6.codec import dump as dump_mf6 from flopy4.mf6.codec import load as load_mf6 from flopy4.mf6.component import Component from flopy4.mf6.converter import structure, unstructure +from flopy4.mf6.ims import Ims +from flopy4.mf6.simulation import Simulation +from flopy4.mf6.tdis import Tdis from flopy4.uio import DEFAULT_REGISTRY +__all__ = ["gwf", "simulation", "solution", "utils", "Ims", "Tdis", "Simulation"] + + +class WriteError(Exception): + """An error occurred while writing a component.""" + + pass + def _load_mf6(path: Path) -> Component: with open(path, "r") as fp: @@ -29,17 +42,26 @@ def _load_toml(path: Path) -> Component: def _write_mf6(component: Component) -> None: with open(component.path, "w") as fp: - dump_mf6(unstructure(component), fp) + data = unstructure(component) + try: + dump_mf6(data, fp) + except Exception as e: + raise WriteError( + f"Failed to write MF6 format file for component '{component.name}' " # type: ignore + f"of type {component.__class__.__name__}" + ) from e def _write_json(component: Component) -> None: with open(component.path, "w") as fp: - dump_json(unstructure(component), fp, indent=4) + data = unstructure(component) + dump_json(data, fp, indent=4) def _write_toml(component: Component) -> None: with open(component.path, "wb") as fp: - dump_toml(unstructure(component), fp) + data = unstructure(component) + dump_toml(data, fp) DEFAULT_REGISTRY.register_loader(Component, "mf6", _load_mf6) diff --git a/flopy4/mf6/codec/writer/filters.py b/flopy4/mf6/codec/writer/filters.py index 8260b4f2..69010135 100644 --- a/flopy4/mf6/codec/writer/filters.py +++ b/flopy4/mf6/codec/writer/filters.py @@ -10,7 +10,7 @@ from flopy4.mf6.constants import FILL_DNODATA -ArrayHow = Literal["constant", "internal", "external"] +ArrayHow = Literal["constant", "internal", "external", "layered internal"] def array_how(value: xr.DataArray) -> ArrayHow: @@ -26,7 +26,11 @@ def array_how(value: xr.DataArray) -> ArrayHow: return "external" if value.max() == value.min(): return "constant" - return "internal" + if value.ndim <= 2: + return "internal" + if value.ndim == 3: + return "layered internal" + raise ValueError(f"Arrays with ndim > 3 are not supported, got ndim={value.ndim}") def array2const(value: xr.DataArray) -> Scalar: @@ -95,7 +99,7 @@ def array2string(value: NDArray) -> str: buffer = StringIO() value = np.asarray(value) if value.ndim > 2: - raise ValueError("Only 1D and 2D arrays are supported.") + raise ValueError(f"Only 1D and 2D arrays are supported, got ndim={value.ndim}") if value.ndim == 1: # add an axis to 1d arrays so np.savetxt writes elements on 1 line value = value[None] @@ -103,7 +107,7 @@ def array2string(value: NDArray) -> str: format = ( "%d" if np.issubdtype(value.dtype, np.integer) - else "%f" + else "%.9f" if np.issubdtype(value.dtype, np.floating) else "%s" ) diff --git a/flopy4/mf6/codec/writer/templates/macros.jinja b/flopy4/mf6/codec/writer/templates/macros.jinja index 45069913..6f6e1a7e 100644 --- a/flopy4/mf6/codec/writer/templates/macros.jinja +++ b/flopy4/mf6/codec/writer/templates/macros.jinja @@ -37,8 +37,16 @@ {% for layer in value -%} {{ inset }}CONSTANT {{ layer|array2const }} {%- endfor %} +{% elif how == "layered internal" %} +{% for layer in value %} + +{{ inset }}INTERNAL +{% for chunk in layer|array2chunks -%} +{{ (2 * inset) ~ chunk|array2string }} +{%- endfor %} +{%- endfor %} {% elif how == "internal" %} -INTERNAL +{{ inset }}INTERNAL {% for chunk in value|array2chunks -%} {{ (2 * inset) ~ chunk|array2string }} {%- endfor %} diff --git a/flopy4/mf6/config.py b/flopy4/mf6/config.py index 58fd9b6b..14da6705 100644 --- a/flopy4/mf6/config.py +++ b/flopy4/mf6/config.py @@ -1,3 +1,3 @@ # TODO use https://environ-config.readthedocs.io/en/stable/? -SPARSE_THRESHOLD = 1000 +SPARSE_THRESHOLD = 100000 diff --git a/flopy4/mf6/context.py b/flopy4/mf6/context.py index f3d38228..2f31e4a0 100644 --- a/flopy4/mf6/context.py +++ b/flopy4/mf6/context.py @@ -7,11 +7,12 @@ from flopy4.mf6.component import Component from flopy4.mf6.constants import MF6 from flopy4.mf6.spec import field +from flopy4.utils import to_path @xattree class Context(Component, ABC): - workspace: Path = field(default=None) + workspace: Path = field(default=None, converter=to_path) def __attrs_post_init__(self): super().__attrs_post_init__() diff --git a/flopy4/mf6/converter/unstructure.py b/flopy4/mf6/converter/unstructure.py index 0a8f076c..76a3f95f 100644 --- a/flopy4/mf6/converter/unstructure.py +++ b/flopy4/mf6/converter/unstructure.py @@ -112,12 +112,18 @@ def oc_setting_data(rec): match value.dtype: case np.bool: # supports boolean dataarrays, e.g. STO steady_state and transient - dat = {kper: "" for kper in range(value.sizes["nper"]) if value.values[kper]} # type: ignore - data[name] = dat + # e.g. steady_state to steady-state, why is't this the dfn name? + fname = name.replace("_", "-") # type: ignore + dat = {kper: "" for kper in range(value.sizes["nper"]) if value.values[kper]} + data[fname] = dat case np.dtypes.StringDType(): # supports string dataarrays, e.g. OC save_budget, save_head fname = name.replace("_", " ") - dat = {kper: value.values[kper] for kper in range(value.sizes["nper"])} + dat = { + kper: value.values[kper] + for kper in range(value.sizes["nper"]) + if value.values[kper].lower() in ["first", "last", "steps", "all"] + } data[fname] = dat case object(): # supports object dataararys, e.g. OC PrintSaveSetting @@ -139,6 +145,8 @@ def oc_setting_data(rec): def unstructure_component(value: Component) -> dict[str, Any]: + from flopy4.mf6.constants import FILL_DNODATA + blockspec = dict(sorted(value.dfn.blocks.items(), key=block_sort_key)) # type: ignore blocks: dict[str, dict[str, Any]] = {} xatspec = xattree.get_xatspec(type(value)) @@ -236,18 +244,23 @@ def unstructure_component(value: Component) -> dict[str, Any]: # setup indexed period blocks, combine arrays into datasets for kper, block in period_blocks.items(): - blocks[f"period {kper + 1}"] = {} + key = f"period {kper + 1}" for arr_name, val in block.items(): - match block[arr_name]: - case str(): - # non data period parameters have their period - # write key set in the _hack_period_non_numeric - # routine - blocks[f"period {kper + 1}"][arr_name] = val - case xr.DataArray(): - blocks[f"period {kper + 1}"]["period"] = xr.Dataset( - block, coords=block[arr_name].coords - ) + if np.any(val != FILL_DNODATA): + # don't create the block (so it isn't written) + # unless there is data to write + if key not in blocks: + blocks[key] = {} + match block[arr_name]: + case str(): + # non data period parameters have their period + # write key set in the _hack_period_non_numeric + # routine + blocks[f"period {kper + 1}"][arr_name] = val + case xr.DataArray(): + blocks[f"period {kper + 1}"]["period"] = xr.Dataset( + block, coords=block[arr_name].coords + ) # combine "perioddata" block arrays (tdis, ats) into datasets # so they render as lists. temp hack TODO do this generically diff --git a/flopy4/mf6/gwf/__init__.py b/flopy4/mf6/gwf/__init__.py index 14ad6512..8c69b9a1 100644 --- a/flopy4/mf6/gwf/__init__.py +++ b/flopy4/mf6/gwf/__init__.py @@ -13,6 +13,7 @@ from flopy4.mf6.gwf.ic import Ic from flopy4.mf6.gwf.npf import Npf from flopy4.mf6.gwf.oc import Oc +from flopy4.mf6.gwf.rch import Rch from flopy4.mf6.gwf.sto import Sto from flopy4.mf6.gwf.wel import Wel from flopy4.mf6.model import Model @@ -20,7 +21,7 @@ from flopy4.mf6.utils import open_cbc, open_hds from flopy4.utils import to_path -__all__ = ["Gwf", "Chd", "Dis", "Drn", "Ic", "Npf", "Oc", "Sto", "Wel"] +__all__ = ["Gwf", "Chd", "Dis", "Drn", "Ic", "Npf", "Oc", "Sto", "Wel", "Rch"] def convert_grid(value): @@ -80,6 +81,7 @@ def budget(self): chd: list[Chd] = field(block="packages") wel: list[Wel] = field(block="packages") drn: list[Drn] = field(block="packages") + rch: list[Rch] = field(block="packages") output: Output = attrs.field( default=attrs.Factory(lambda self: Gwf.Output(self), takes_self=True) ) diff --git a/flopy4/mf6/simulation.py b/flopy4/mf6/simulation.py index 86a58498..2ddebde5 100644 --- a/flopy4/mf6/simulation.py +++ b/flopy4/mf6/simulation.py @@ -34,11 +34,14 @@ def default_filename(self) -> str: def __attrs_post_init__(self): super().__attrs_post_init__() if self.filename != "mfsim.nam": - warn( - "Simulation filename must be 'mfsim.nam'.", - UserWarning, - ) + if self.filename is not None: + warn( + "Simulation filename must be 'mfsim.nam'.", + UserWarning, + ) self.filename = "mfsim.nam" + for model in self.models.values(): + model.workspace = self.workspace @property def time(self) -> Time: diff --git a/flopy4/mf6/utils/__init__.py b/flopy4/mf6/utils/__init__.py index be3986c4..03b73cae 100644 --- a/flopy4/mf6/utils/__init__.py +++ b/flopy4/mf6/utils/__init__.py @@ -1,4 +1,7 @@ +# Import submodules to make them accessible via flopy4.mf6.utils.* +from flopy4.mf6.utils import grid, time + from .cbc_reader import open_cbc from .heads_reader import open_hds -__all__ = ["open_hds", "open_cbc"] +__all__ = ["open_hds", "open_cbc", "grid", "time"] diff --git a/flopy4/mf6/utils/grid.py b/flopy4/mf6/utils/grid.py index 4bd3fb5b..fa28f21b 100644 --- a/flopy4/mf6/utils/grid.py +++ b/flopy4/mf6/utils/grid.py @@ -2,16 +2,139 @@ from typing import Any import numpy as np +import xarray as xr from attrs import fields -from flopy.discretization import StructuredGrid +from flopy.discretization import StructuredGrid as LegacyStructuredGrid +from xarray.core.indexes import PandasIndex from flopy4.mf6.constants import FILL_DNODATA +class StructuredGrid(LegacyStructuredGrid): + """Extend flopy3's StructuredGrid""" + + # TODO de-duplicate array/dataset setup. no need to do it in both __init__ and properties + + def __init__(self, *args, **kwargs): + super().__init__(*args, **kwargs) + self._dims_coords = { + "nlay": "k", + "nrow": "i", + "ncol": "j", + "nodes": "node", + } + self._coords = { + "k": xr.DataArray(np.arange(self.nlay, dtype=int), dims=("nlay",)), + "i": xr.DataArray(np.arange(self.nrow, dtype=int), dims=("nrow",)), + "j": xr.DataArray(np.arange(self.ncol, dtype=int), dims=("ncol",)), + "node": xr.DataArray(np.arange(self.nnodes, dtype=int), dims=("nodes",)), + } + data_vars = { + "delr": self.delr, + "delc": self.delc, + "delz": self.delz, + "top": self.top, + "botm": self.botm, + "idomain": self.idomain, + } + self._dataset = ( + xr.Dataset({k: v for k, v in data_vars.items() if v is not None}, coords=self._coords) + .set_xindex("k", PandasIndex) + .set_xindex("i", PandasIndex) + .set_xindex("j", PandasIndex) + .set_xindex("node", PandasIndex) + ) + + @property + def dataset(self) -> xr.Dataset: + return self._dataset + + @property + def delc(self): + if self.__delc is None: + return None + dims = ("ncol",) + coord_name = self._dims_coords[dims[0]] + coords = coords = {coord_name: self._coords[coord_name]} + return xr.DataArray(super().delc, coords=coords, dims=dims).set_xindex( + coord_name, PandasIndex + ) + + @property + def delr(self): + if self.__delr is None: + return None + dims = ("nrow",) + coord_name = self._dims_coords[dims[0]] + coords = {coord_name: self._coords[coord_name]} + return xr.DataArray(super().delr, coords=coords, dims=dims).set_xindex( + coord_name, PandasIndex + ) + + @property + def delz(self): + dims = ("nlay", "nrow", "ncol") + coord_names = ( + self._dims_coords[dims[0]], + self._dims_coords[dims[1]], + self._dims_coords[dims[2]], + ) + coords = {coord_name: self._coords[coord_name] for coord_name in coord_names} + return ( + xr.DataArray(super().delz, coords=coords, dims=dims) + .set_xindex(coord_names[0], PandasIndex) + .set_xindex(coord_names[1], PandasIndex) + .set_xindex(coord_names[2], PandasIndex) + ) + + @property + def top(self): + dims = ("nrow", "ncol") + coord_names = (self._dims_coords[dims[0]], self._dims_coords[dims[1]]) + coords = {coord_name: self._coords[coord_name] for coord_name in coord_names} + return ( + xr.DataArray(super().top, coords=coords, dims=dims) + .set_xindex(coord_names[0], PandasIndex) + .set_xindex(coord_names[1], PandasIndex) + ) + + @property + def botm(self): + dims = ("nlay", "nrow", "ncol") + coord_names = ( + self._dims_coords[dims[0]], + self._dims_coords[dims[1]], + self._dims_coords[dims[2]], + ) + coords = {coord_name: self._coords[coord_name] for coord_name in coord_names} + return ( + xr.DataArray(super().botm, coords=coords, dims=dims) + .set_xindex(coord_names[0], PandasIndex) + .set_xindex(coord_names[1], PandasIndex) + .set_xindex(coord_names[2], PandasIndex) + ) + + @property + def idomain(self): + dims = ("nlay", "nrow", "ncol") + coord_names = ( + self._dims_coords[dims[0]], + self._dims_coords[dims[1]], + self._dims_coords[dims[2]], + ) + coords = {coord_name: self._coords[coord_name] for coord_name in coord_names} + return ( + xr.DataArray(super().idomain, coords=coords, dims=dims) + .set_xindex(coord_names[0], PandasIndex) + .set_xindex(coord_names[1], PandasIndex) + .set_xindex(coord_names[2], PandasIndex) + ) + + def get_coords(grid: StructuredGrid) -> dict[str, Any]: # unpack tuples xmin, xmax, ymin, ymax = grid.extent - dx, dy = (grid.delr, -grid.delc) + dx, dy = (grid.delr, -grid.delc) # type: ignore coords: collections.OrderedDict[str, Any] = collections.OrderedDict() # from cell size to x and y coordinates if isinstance(dx, (int, float, np.int_)): # equidistant diff --git a/test/conftest.py b/test/conftest.py index cdc59fa9..79316b7c 100644 --- a/test/conftest.py +++ b/test/conftest.py @@ -11,7 +11,6 @@ PROJ_ROOT_PATH = Path(__file__).parents[1] DOCS_PATH = PROJ_ROOT_PATH / "docs" EXAMPLES_PATH = DOCS_PATH / "examples" -EXCLUDED_EXAMPLES = [] def _get_dfn_path(tmp_path_factory): @@ -59,13 +58,3 @@ def patch_macos_ci_matplotlib(): import matplotlib matplotlib.use("agg") - - -def pytest_generate_tests(metafunc): - if "example_script" in metafunc.fixturenames: - scripts = { - file.name: file - for file in sorted(EXAMPLES_PATH.glob("*example.py")) - if file.stem not in EXCLUDED_EXAMPLES - } - metafunc.parametrize("example_script", scripts.values(), ids=scripts.keys()) diff --git a/test/test_examples.py b/test/test_examples.py index 5a008066..46142fe5 100644 --- a/test/test_examples.py +++ b/test/test_examples.py @@ -3,8 +3,20 @@ import sys import pytest -from modflow_devtools.markers import requires_exe -from modflow_devtools.misc import cd, run_cmd +from conftest import EXAMPLES_PATH +from modflow_devtools.misc import run_cmd + +EXCLUDE = ["quickstart_expanded"] + + +def pytest_generate_tests(metafunc): + if "example_script" in metafunc.fixturenames: + scripts = { + file.name: file + for file in sorted(EXAMPLES_PATH.glob("*.py")) + if file.stem not in EXCLUDE + } + metafunc.parametrize("example_script", scripts.values(), ids=scripts.keys()) @pytest.mark.slow @@ -12,20 +24,3 @@ def test_scripts(example_script): args = [sys.executable, example_script] stdout, stderr, retcode = run_cmd(*args, verbose=True) assert not retcode, stdout + stderr - - -@pytest.mark.slow -@requires_exe("jupytext") -def test_notebooks(example_script, tmp_path): - args = [ - "jupytext", - "--from", - "py", - "--to", - "ipynb", - "--execute", - example_script, - ] - with cd(tmp_path): - out, err, ret = run_cmd(*args, verbose=True) - assert not ret, out + err diff --git a/test/test_mf6_codec.py b/test/test_mf6_codec.py index 796c3333..84db7e65 100644 --- a/test/test_mf6_codec.py +++ b/test/test_mf6_codec.py @@ -2,6 +2,8 @@ from pprint import pprint +import pytest + from flopy4.mf6.codec import dumps, loads from flopy4.mf6.converter import COMPONENT_CONVERTER @@ -71,7 +73,7 @@ def test_dumps_sto(): print("STO dump:") print(dumped) assert "BEGIN PERIOD 1\n TRANSIENT" in dumped - assert "BEGIN PERIOD 2\n STEADY_STATE" in dumped + assert "BEGIN PERIOD 2\n STEADY-STATE" in dumped assert "BEGIN PERIOD 3\n TRANSIENT" in dumped assert dumped @@ -148,11 +150,12 @@ def test_dumps_oc2(): pprint(loaded) -def test_dumps_dis(): +@pytest.fixture +def dis_with_constant_arrays(): from flopy4.mf6.gwf import Dis - dis = Dis( - nlay=1, + return Dis( + nlay=2, nrow=10, ncol=10, delr=100.0, @@ -161,6 +164,9 @@ def test_dumps_dis(): length_units="feet", ) + +def test_dumps_dis_with_constant_arrays(dis_with_constant_arrays): + dis = dis_with_constant_arrays dumped = dumps(COMPONENT_CONVERTER.unstructure(dis)) print("DIS dump:") print(dumped) @@ -171,11 +177,26 @@ def test_dumps_dis(): pprint(loaded) assert ["LENGTH_UNITS", "feet"] in loaded["OPTIONS"] - assert loaded["DIMENSIONS"] == [["NLAY", 1], ["NCOL", 10], ["NROW", 10]] + assert loaded["DIMENSIONS"] == [["NLAY", 2], ["NCOL", 10], ["NROW", 10]] assert ["DELR"] in loaded["GRIDDATA"] assert ["DELC"] in loaded["GRIDDATA"] +def test_dumps_dis_with_layered_arrays(dis_with_constant_arrays): + dis = dis_with_constant_arrays + dis.delr[0] = 101.0 + dis.botm[0, 0, 0] = -1.0 # 3d array will force layered output + dumped = dumps(COMPONENT_CONVERTER.unstructure(dis)) + print("DIS dump:") + print(dumped) + assert dumped + assert "BOTM LAYERED" in dumped + + loaded = loads(dumped) + print("DIS load:") + pprint(loaded) + + def test_dumps_tdis(): from flopy4.mf6.tdis import Tdis from flopy4.mf6.utils.time import Time diff --git a/test/test_mf6_component.py b/test/test_mf6_component.py index e1a753b6..03cd9ff1 100644 --- a/test/test_mf6_component.py +++ b/test/test_mf6_component.py @@ -216,13 +216,13 @@ def test_init_sim_explicit_dims(): def test_init_big_sim(): # if size over threshold, arrays should be sparse time = Time(perlen=[1.0], nstp=[1], tsmult=[1.0]) - grid = StructuredGrid(nlay=1, nrow=100, ncol=100) + grid = StructuredGrid(nlay=1, nrow=1000, ncol=1000) sim = Simulation(tdis=time) gwf = Gwf(parent=sim, dis=grid) ic = Ic(parent=gwf) oc = Oc(parent=gwf) npf = Npf(parent=gwf) - chd = Chd(parent=gwf, head={"*": {(0, 0, 0): 1.0, (0, 99, 99): 0.0}}) + chd = Chd(parent=gwf, head={"*": {(0, 0, 0): 1.0, (0, 999, 999): 0.0}}) assert sim.models["gwf"] is gwf assert isinstance(sim.data, DataTree) @@ -231,11 +231,11 @@ def test_init_big_sim(): assert gwf.oc is oc assert gwf.npf is npf assert gwf.chd[0] is chd - assert np.array_equal(sim.models["gwf"].npf.k, np.ones(10000)) - assert np.array_equal(sim.models["gwf"].npf.data.k, np.ones(10000)) + assert np.array_equal(sim.models["gwf"].npf.k, np.ones(1000000)) + assert np.array_equal(sim.models["gwf"].npf.data.k, np.ones(1000000)) assert chd.head[0, 0].item() == 1.0 - assert chd.head[0, 9999].item() == 0.0 - assert np.array_equal(chd.head[0, 1:9999].data.todense(), np.full((9998,), FILL_DNODATA)) + assert chd.head[0, 999999].item() == 0.0 + assert np.array_equal(chd.head[0, 1:999999].data.todense(), np.full((999998,), FILL_DNODATA)) assert np.array_equal(chd.head.data.todense(), chd.data.head.data.todense()) assert np.array_equal( chd.head.data.todense(),