From baf4300049be8b677d934da48ac15dfbf094c6fb Mon Sep 17 00:00:00 2001 From: jdhughes-usgs Date: Tue, 9 Jul 2024 21:59:19 -0500 Subject: [PATCH] fix(csub): fix a few memory issues in csub package (#1935) * fix z-displacement issue with inactive cells when saving z-displacement data * initialize flag_string * add warning if saving delay bed convergence data and delay beds are not included in a simulation * add a simple csub test that covers the changes (inactive cells, z-displacement output, and csub convergence output w/no delay beds) --- autotest/test_gwf_csub_zdisp02.py | 280 +++++++++++++++++++++++++ doc/ReleaseNotes/develop.tex | 4 + doc/mf6io/mf6ivar/dfn/gwf-csub.dfn | 2 +- src/Model/GroundWaterFlow/gwf-csub.f90 | 159 ++++++++------ 4 files changed, 377 insertions(+), 68 deletions(-) create mode 100644 autotest/test_gwf_csub_zdisp02.py diff --git a/autotest/test_gwf_csub_zdisp02.py b/autotest/test_gwf_csub_zdisp02.py new file mode 100644 index 00000000000..a5a506ae1b3 --- /dev/null +++ b/autotest/test_gwf_csub_zdisp02.py @@ -0,0 +1,280 @@ +import pathlib as pl + +import flopy +import numpy as np +import pytest + +from framework import TestFramework + +cases = ["csub_zdisp02"] +htol = [None for _ in range(len(cases))] +dtol = 1e-3 +budtol = 1e-2 + +# static model data +# temporal discretization +nper = 31 +perlen = [1.0] + [365.2500000 for _ in range(nper - 1)] +nstp = [1] + [6 for _ in range(nper - 1)] +tsmult = [1.0] + [1.3 for _ in range(nper - 1)] +tdis_rc = [] +for i in range(nper): + tdis_rc.append((perlen[i], nstp[i], tsmult[i])) + +# spatial discretization data +nlay, nrow, ncol = 3, 6, 6 +shape3d = (nlay, nrow, ncol) +size3d = nlay * nrow * ncol +delr, delc = 1000.0, 1000.0 +top = 0.0 +botm = [-100, -150.0, -350.0] +zthick = [top - botm[0], botm[0] - botm[1], botm[1] - botm[2]] +strt = 0.0 + +# create idomain and ibound +idomain = np.ones(shape3d, dtype=np.int32) +idomain[0, :, :] = np.array( + [ + [0, 0, 0, 1, 1, 1], + [0, 0, 1, 1, 1, 1], + [0, 1, 1, 1, 1, 1], + [1, 1, 1, 1, 1, 0], + [1, 1, 1, 1, 0, 0], + [1, 1, 1, 0, 0, 0], + ] +) +idomain[1, :, :] = np.flip(idomain[0], axis=1) +idomain[2] = idomain[0] * idomain[1] + +# calculate hk +hk1fact = 1.0 / 50.0 +hk1 = 0.5 * hk1fact +hk = [20.0, hk1, 5.0] + +# calculate vka +vka = [1e6, 7.5e-5, 1e6] + +# layer 1 is convertible +laytyp = [1, 0, 0] + +# solver options +nouter, ninner = 500, 300 +hclose, rclose, relax = 1e-9, 1e-6, 1.0 +newtonoptions = "NEWTON" +imsla = "BICGSTAB" + + +# pumping well data +spd_wel = [(2, 3, 2, -1400.0)] +maxwel = len(spd_wel) + +# storage and compaction data +ss = [0.0, 0.0, 0.0] +void = 0.82 +theta = void / (1.0 + void) + +# static ibc and sub data +sgm = 0.0 +sgs = 0.0 +omega = 1.0 + +# no delay bed data +lnd = [0, 1, 2] +hc = -7.0 +thicknd0 = [15.0, 50.0, 30.0] +ccnd0 = [6e-4, 3e-4, 6e-4] +crnd0 = [6e-6, 3e-6, 6e-6] +sfv = [] +sfe = [] +for k in range(nlay): + sfv.append(ccnd0[k] * thicknd0[k]) + sfe.append(crnd0[k] * thicknd0[k]) + +# no delay bed packagedata entries +sub6 = [] +cdelays = "nodelay" +ibcno = 0 +for kdx, k in enumerate(lnd): + for i in range(nrow): + for j in range(ncol): + # skip inactive cells + if idomain[k, i, j] == 0: + continue + tag = f"{k + 1:02d}_{i + 1:02d}_{j + 1:02d}" + # create nodelay entry + b = thicknd0[kdx] + d = ( + ibcno, + (k, i, j), + cdelays, + hc, + b, + 1.0, + ccnd0[kdx], + crnd0[kdx], + theta, + 999.0, + -999.0, + tag, + ) + sub6.append(d) + ibcno += 1 + +# create coarse-grained material storage +ske_scaled = [] +for k in range(nlay): + sst = (zthick[k] - thicknd0[k]) * ss[k] / zthick[k] + ske_scaled.append(sst) + + +def build_models(idx, test): + name = cases[idx] + + # build MODFLOW 6 files + ws = test.workspace + sim = flopy.mf6.MFSimulation( + sim_name=name, version="mf6", exe_name="mf6", sim_ws=ws + ) + + flopy.mf6.ModflowTdis( + sim, time_units="DAYS", nper=nper, perioddata=tdis_rc + ) + + flopy.mf6.ModflowIms( + sim, + print_option="SUMMARY", + outer_dvclose=hclose, + outer_maximum=nouter, + under_relaxation="NONE", + inner_maximum=ninner, + inner_dvclose=hclose, + rcloserecord=rclose, + linear_acceleration=imsla, + scaling_method="NONE", + reordering_method="NONE", + relaxation_factor=relax, + ) + + # create gwf model + gwf = flopy.mf6.ModflowGwf( + sim, modelname=name, newtonoptions=newtonoptions, save_flows=True + ) + + flopy.mf6.ModflowGwfdis( + gwf, + nlay=nlay, + nrow=nrow, + ncol=ncol, + delr=delr, + delc=delc, + top=top, + botm=botm, + idomain=idomain, + ) + + flopy.mf6.ModflowGwfic(gwf, strt=strt, filename=f"{name}.ic") + + flopy.mf6.ModflowGwfnpf( + gwf, save_flows=False, icelltype=laytyp, k=hk, k33=vka + ) + + flopy.mf6.ModflowGwfsto( + gwf, + save_flows=False, + iconvert=laytyp, + ss=0, + sy=0, + storagecoefficient=None, + steady_state={0: True}, + transient={1: True}, + ) + + copth = f"{name}.compaction.gridbin" + zopth = f"{name}.zdisplacement.gridbin" + cvgpth = f"{name}.csub.convergence.csv" + flopy.mf6.ModflowGwfcsub( + gwf, + boundnames=True, + head_based=True, + update_material_properties=True, + specified_initial_interbed_state=True, + effective_stress_lag=False, + save_flows=True, + compaction_filerecord=copth, + zdisplacement_filerecord=zopth, + package_convergence_filerecord=cvgpth, + ninterbeds=len(sub6), + beta=0.0, + cg_ske_cr=ss, + packagedata=sub6, + ) + + flopy.mf6.ModflowGwfwel( + gwf, + print_input=True, + print_flows=True, + maxbound=maxwel, + stress_period_data={1: spd_wel}, + ) + + flopy.mf6.ModflowGwfoc( + gwf, + budget_filerecord=f"{name}.cbc", + head_filerecord=f"{name}.hds", + saverecord=[("HEAD", "ALL"), ("BUDGET", "ALL")], + printrecord=[("HEAD", "LAST"), ("BUDGET", "ALL")], + ) + + return sim, None + + +def check_output(idx, test): + with open(pl.Path(test.workspace) / "mfsim.lst", "r") as f: + lines = f.readlines() + tag = ( + "1. Package convergence data is requested but delay interbeds are not" + ) + found = False + for line in lines[::-1]: + if tag in line: + found = True + break + assert found, f"WARNING: '{tag}' not found in mfsim.lst" + + sim = flopy.mf6.MFSimulation.load(sim_ws=test.workspace) + gwf = sim.get_model() + idomain = gwf.dis.idomain.array + + times = gwf.csub.output.compaction().get_times() + totim = times[-1] + + compaction = gwf.csub.output.compaction().get_data(totim=totim) + compaction[compaction == 1e30] = 0.0 + zdis_calc = np.zeros(compaction.shape, dtype=float) + zdis_calc[2] = compaction[2] + for k in (1, 0): + zdis_calc[k] = zdis_calc[k + 1] + compaction[k] + zdis_calc[idomain < 1] = 0.0 + + zpth = pl.Path(test.workspace) / "csub_zdisp02.zdisplacement.gridbin" + zobj = flopy.utils.HeadFile(zpth, text="CSUB-ZDISPLACE") + zdis = zobj.get_data(totim=totim) + zdis[zdis == 1e30] = 0.0 + + assert np.allclose( + zdis, zdis_calc + ), "Calculated z-displacement is not equal to simulated z-displacement" + + +@pytest.mark.slow +@pytest.mark.parametrize("idx, name", enumerate(cases)) +def test_mf6model(idx, name, function_tmpdir, targets): + test = TestFramework( + name=name, + workspace=function_tmpdir, + targets=targets, + build=lambda t: build_models(idx, t), + check=lambda t: check_output(idx, t), + htol=htol[idx], + ) + test.run() diff --git a/doc/ReleaseNotes/develop.tex b/doc/ReleaseNotes/develop.tex index d179f99a402..9e159f3d060 100644 --- a/doc/ReleaseNotes/develop.tex +++ b/doc/ReleaseNotes/develop.tex @@ -29,6 +29,10 @@ \item The PRT model's cell face flows were improperly combined with boundary flows; for cell faces with active neighbors, the face flow replaced any boundary flows (likely a rare situation because IFLOWFACE is typically applied to faces on the boundary of the active domain). The face flow calculation has been corrected. \item A bad pointer reference has been fixed in the PRT model. Depending on the combination of platform and compiler used to build the program, this could result in crashes (e.g. divide by zero errors) or in silent failures to properly pass particles downward between vertically adjacent cells, leading to early termination. The latter could occur as a consequence of undefined behavior which prevented detection of situations when a particle should exit a cell through its bottom face. \item For a UZF setup with multiple UZF objects in a cell, the auxmultname option would cause the program to issue an unnecessary error. The program was corrected to remove the unnecessary error. The revised program was tested to ensure that the assignment of multiple UZF objects per cell works when it should and fails when the cumulative area of the UZF objects within a cell exceed the total area for the cell. + \item An array out of bounds error for Z-displacement output in inactive cells has been fixed in the CSUB package. Depending on the combination of platform and compiler used to build the program, this could result in crashes. + \item Initialize a few uninitialized variables in the CSUB package. Depending on the combination of platform and compiler used to build the program, this could result in different program behaviour. + \item Add a warning if saving convergence data for the CSUB package when delay beds are not included in a simulation. Also modified the convergence output data so that it is clear that DVMAX and DSTORAGEMAX data and locations are not calculated in this case (by writing `-\,-' to the output file). + % \item xxx \end{itemize} diff --git a/doc/mf6io/mf6ivar/dfn/gwf-csub.dfn b/doc/mf6io/mf6ivar/dfn/gwf-csub.dfn index 954240daaf3..421456b2f18 100644 --- a/doc/mf6io/mf6ivar/dfn/gwf-csub.dfn +++ b/doc/mf6io/mf6ivar/dfn/gwf-csub.dfn @@ -413,7 +413,7 @@ reader urword tagged true optional false longname package_convergence keyword -description keyword to specify that record corresponds to the package convergence comma spaced values file. +description keyword to specify that record corresponds to the package convergence comma spaced values file. Package convergence data is for delay interbeds. A warning message will be issued if package convergence data is requested but delay interbeds are not included in the package. block options name package_convergence_filename diff --git a/src/Model/GroundWaterFlow/gwf-csub.f90 b/src/Model/GroundWaterFlow/gwf-csub.f90 index f6a56bc07c8..1f9831f4c47 100644 --- a/src/Model/GroundWaterFlow/gwf-csub.f90 +++ b/src/Model/GroundWaterFlow/gwf-csub.f90 @@ -7,7 +7,7 @@ !! !< module GwfCsubModule - use KindModule, only: I4B, DP + use KindModule, only: I4B, DP, LGP use ConstantsModule, only: DPREC, DZERO, DEM20, DEM15, DEM10, DEM8, DEM7, & DEM6, DEM4, DP9, DHALF, DEM1, DONE, DTWO, DTHREE, & DGRAVITY, DTEN, DHUNDRED, DNODATA, DHNOFLO, & @@ -80,7 +80,7 @@ module GwfCsubModule character(len=LENAUXNAME), dimension(:), & pointer, contiguous :: auxname => null() !< vector of auxname ! -- logical scalars - logical, pointer :: lhead_based => null() !< logical variable indicating if head-based solution + logical(LGP), pointer :: lhead_based => null() !< logical variable indicating if head-based solution ! -- integer scalars integer(I4B), pointer :: istounit => null() !< unit number of storage package integer(I4B), pointer :: istrainib => null() !< unit number of interbed strain output @@ -303,6 +303,9 @@ module GwfCsubModule procedure, private :: csub_delay_assemble_fc procedure, private :: csub_delay_assemble_fn procedure, private :: csub_delay_head_check + + ! methods for tables + procedure, private :: csub_initialize_tables ! ! -- methods for observations procedure, public :: csub_obs_supported @@ -367,7 +370,7 @@ subroutine csub_ar(this, dis, ibound) class(DisBaseType), pointer, intent(in) :: dis !< model discretization integer(I4B), dimension(:), pointer, contiguous :: ibound !< model ibound array ! -- local variables - logical :: isfound, endOfBlock + logical(LGP) :: isfound, endOfBlock character(len=:), allocatable :: line character(len=LINELENGTH) :: keyword character(len=20) :: cellid @@ -419,6 +422,9 @@ subroutine csub_ar(this, dis, ibound) ! ! - observation data call this%obs%obs_ar() + + ! setup tables + call this%csub_initialize_tables() ! ! -- terminate if errors dimensions block data if (count_errors() > 0) then @@ -610,8 +616,8 @@ subroutine read_options(this) character(len=:), allocatable :: line character(len=MAXCHARLEN) :: fname character(len=LENAUXNAME), dimension(:), allocatable :: caux - logical :: isfound - logical :: endOfBlock + logical(LGP) :: isfound + logical(LGP) :: endOfBlock integer(I4B) :: n integer(I4B) :: lloc integer(I4B) :: istart @@ -1048,7 +1054,7 @@ subroutine csub_read_dimensions(this) ! -- local variables character(len=LENBOUNDNAME) :: keyword integer(I4B) :: ierr - logical :: isfound, endOfBlock + logical(LGP) :: isfound, endOfBlock ! -- format ! ! -- initialize dimensions to -1 @@ -1421,8 +1427,8 @@ subroutine csub_read_packagedata(this) character(len=10) :: text character(len=LENBOUNDNAME) :: bndName character(len=7) :: cdelay - logical :: isfound - logical :: endOfBlock + logical(LGP) :: isfound + logical(LGP) :: endOfBlock integer(I4B) :: ival integer(I4B) :: n integer(I4B) :: nn @@ -2514,8 +2520,8 @@ subroutine csub_rp(this) character(len=LINELENGTH) :: title character(len=LINELENGTH) :: text character(len=20) :: cellid - logical :: isfound - logical :: endOfBlock + logical(LGP) :: isfound + logical(LGP) :: endOfBlock integer(I4B) :: jj integer(I4B) :: ierr integer(I4B) :: node @@ -3001,6 +3007,60 @@ subroutine csub_fn(this, kiter, hold, hnew, matrix_sln, idxglo, rhs) return end subroutine csub_fn + !> @ brief Initialize optional tables + !! + !! Subroutine to initialize optional tables. Tables include: + !! o delay interbeds convergence tables + !! + !< + subroutine csub_initialize_tables(this) + class(GwfCsubType) :: this + + character(len=LINELENGTH) :: tag + integer(I4B) :: ntabrows + integer(I4B) :: ntabcols + + if (this%ipakcsv > 0) then + if (this%ndelaybeds < 1) then + write (warnmsg, '(a,1x,3a)') & + 'Package convergence data is requested but delay interbeds', & + 'are not included in package (', & + trim(adjustl(this%packName)), ').' + call store_warning(warnmsg) + end if + + ntabrows = 1 + ntabcols = 9 + + ! setup table + call table_cr(this%pakcsvtab, this%packName, '') + call this%pakcsvtab%table_df(ntabrows, ntabcols, this%ipakcsv, & + lineseparator=.FALSE., separator=',', & + finalize=.FALSE.) + + ! add columns to package csv + tag = 'total_inner_iterations' + call this%pakcsvtab%initialize_column(tag, 10, alignment=TABLEFT) + tag = 'totim' + call this%pakcsvtab%initialize_column(tag, 10, alignment=TABLEFT) + tag = 'kper' + call this%pakcsvtab%initialize_column(tag, 10, alignment=TABLEFT) + tag = 'kstp' + call this%pakcsvtab%initialize_column(tag, 10, alignment=TABLEFT) + tag = 'nouter' + call this%pakcsvtab%initialize_column(tag, 10, alignment=TABLEFT) + tag = 'dvmax' + call this%pakcsvtab%initialize_column(tag, 15, alignment=TABLEFT) + tag = 'dvmax_loc' + call this%pakcsvtab%initialize_column(tag, 15, alignment=TABLEFT) + tag = 'dstoragemax' + call this%pakcsvtab%initialize_column(tag, 15, alignment=TABLEFT) + tag = 'dstoragemax_loc' + call this%pakcsvtab%initialize_column(tag, 15, alignment=TABLEFT) + end if + + end subroutine csub_initialize_tables + !> @ brief Final convergence check !! !! Final convergence check for the CSUB package. The final convergence @@ -3030,13 +3090,10 @@ subroutine csub_cc(this, innertot, kiter, iend, icnvgmod, nodes, & character(len=LENPAKLOC), intent(inout) :: cpak !< string location of the maximum change in csub package integer(I4B), intent(inout) :: ipak !< node with the maximum change in csub package real(DP), intent(inout) :: dpak !< maximum change in csub package - ! -- local variables - character(len=LINELENGTH) :: tag + ! local variables character(len=LENPAKLOC) :: cloc integer(I4B) :: icheck integer(I4B) :: ipakfail - integer(I4B) :: ntabrows - integer(I4B) :: ntabcols integer(I4B) :: ib integer(I4B) :: node integer(I4B) :: idelay @@ -3066,57 +3123,17 @@ subroutine csub_cc(this, innertot, kiter, iend, icnvgmod, nodes, & ipakfail = 0 locdhmax = 0 locrmax = 0 + ifirst = 1 dhmax = DZERO rmax = DZERO - ifirst = 1 ! ! -- additional checks to see if convergence needs to be checked ! -- no convergence check for steady-state stress periods if (this%gwfiss /= 0) then icheck = 0 else - ! - ! -- if not saving package convergence data on check convergence if - ! the model is considered converged - if (this%ipakcsv == 0) then - if (icnvgmod == 0) then - icheck = 0 - end if - else - ! - ! -- header for package csv - if (.not. associated(this%pakcsvtab)) then - ! - ! -- determine the number of columns and rows - ntabrows = 1 - ntabcols = 9 - ! - ! -- setup table - call table_cr(this%pakcsvtab, this%packName, '') - call this%pakcsvtab%table_df(ntabrows, ntabcols, this%ipakcsv, & - lineseparator=.FALSE., separator=',', & - finalize=.FALSE.) - ! - ! -- add columns to package csv - tag = 'total_inner_iterations' - call this%pakcsvtab%initialize_column(tag, 10, alignment=TABLEFT) - tag = 'totim' - call this%pakcsvtab%initialize_column(tag, 10, alignment=TABLEFT) - tag = 'kper' - call this%pakcsvtab%initialize_column(tag, 10, alignment=TABLEFT) - tag = 'kstp' - call this%pakcsvtab%initialize_column(tag, 10, alignment=TABLEFT) - tag = 'nouter' - call this%pakcsvtab%initialize_column(tag, 10, alignment=TABLEFT) - tag = 'dvmax' - call this%pakcsvtab%initialize_column(tag, 15, alignment=TABLEFT) - tag = 'dvmax_loc' - call this%pakcsvtab%initialize_column(tag, 15, alignment=TABLEFT) - tag = 'dstoragemax' - call this%pakcsvtab%initialize_column(tag, 15, alignment=TABLEFT) - tag = 'dstoragemax_loc' - call this%pakcsvtab%initialize_column(tag, 15, alignment=TABLEFT) - end if + if (icnvgmod == 0) then + icheck = 0 end if end if ! @@ -3214,10 +3231,17 @@ subroutine csub_cc(this, innertot, kiter, iend, icnvgmod, nodes, & call this%pakcsvtab%add_term(kper) call this%pakcsvtab%add_term(kstp) call this%pakcsvtab%add_term(kiter) - call this%pakcsvtab%add_term(dhmax) - call this%pakcsvtab%add_term(locdhmax) - call this%pakcsvtab%add_term(rmax) - call this%pakcsvtab%add_term(locrmax) + if (this%ndelaybeds > 0) then + call this%pakcsvtab%add_term(dhmax) + call this%pakcsvtab%add_term(locdhmax) + call this%pakcsvtab%add_term(rmax) + call this%pakcsvtab%add_term(locrmax) + else + call this%pakcsvtab%add_term('--') + call this%pakcsvtab%add_term('--') + call this%pakcsvtab%add_term('--') + call this%pakcsvtab%add_term('--') + end if ! ! -- finalize the package csv if (iend == 1) then @@ -3793,7 +3817,7 @@ subroutine csub_ot_dv(this, idvfl, idvprint) ! -- fill buff with data from buffusr do nodeu = 1, this%dis%nodesuser node = this%dis%get_nodenumber_idx1(nodeu, 1) - if (node /= 0) then + if (node > 0) then this%buff(node) = this%buffusr(nodeu) end if end do @@ -5765,11 +5789,11 @@ subroutine csub_delay_sln(this, ib, hcell, update) class(GwfCsubType), intent(inout) :: this integer(I4B), intent(in) :: ib !< interbed number real(DP), intent(in) :: hcell !< current head in a cell - logical, intent(in), optional :: update !< optional logical variable indicating - !! if the maximum head change variable - !! in a delay bed should be updated + logical(LGP), intent(in), optional :: update !< optional logical variable indicating + !! if the maximum head change variable + !! in a delay bed should be updated ! -- local variables - logical :: lupdate + logical(LGP) :: lupdate integer(I4B) :: n integer(I4B) :: icnvg integer(I4B) :: iter @@ -7512,10 +7536,11 @@ subroutine csub_process_obsID(obsrv, dis, inunitobs, iout) integer(I4B) :: icol, istart, istop character(len=LINELENGTH) :: string character(len=LENBOUNDNAME) :: bndname - logical :: flag_string + logical(LGP) :: flag_string ! ! -- initialize variables string = obsrv%IDstring + flag_string = .TRUE. ! ! -- Extract reach number from string and store it. ! If 1st item is not an integer(I4B), it should be a