From d572d0c154a8abe43d35eb8c4645df7cf925cd52 Mon Sep 17 00:00:00 2001 From: emorway-usgs Date: Mon, 25 Nov 2024 16:18:56 -0800 Subject: [PATCH] fix(gwe-uze): stop simulation when multiple uze objects exist in one cell --- autotest/test_gwe_uze_auxmult.py | 470 ++++++++++++++++++++++++ doc/ReleaseNotes/develop.tex | 1 + doc/mf6io/gwe/uze.tex | 2 +- src/Model/GroundWaterEnergy/gwe-uze.f90 | 68 +++- src/Model/GroundWaterFlow/gwf-uzf.f90 | 2 +- 5 files changed, 538 insertions(+), 5 deletions(-) create mode 100644 autotest/test_gwe_uze_auxmult.py diff --git a/autotest/test_gwe_uze_auxmult.py b/autotest/test_gwe_uze_auxmult.py new file mode 100644 index 00000000000..defd5e7d890 --- /dev/null +++ b/autotest/test_gwe_uze_auxmult.py @@ -0,0 +1,470 @@ +""" +A simple test of confirm that when the auxmultname option is active in UZF, +UZE will stop with an error msg. Multiple UZF objects in a cell is not +supported in UZE. Two variations of the same test problem are scripted; the +first test combines GWF and GWE into a single simulation and should error out. +The second test splits GWF and GWE into two different simulations where the +simulation containing the GWE model should error out. +""" + +import flopy +import numpy as np + +cases = ["uzauxmlt-err", "uzauxmlt-fmi"] + +nouter, ninner = 100, 300 +hclose, rclose, relax = 1e-9, 1e-3, 0.97 + +uzarea_data = { + "uzf_pkdat": [ + [0, (0, 0, 3), 1, -1, 1.0, 1e-5, 0.2, 0.3, 0.25, 3.5, "uzf01"], + [1, (0, 0, 3), 1, -1, 1.0, 1e-5, 0.2, 0.3, 0.25, 3.5, "uzf01"], + [2, (0, 0, 4), 1, -1, 1.0, 1e-5, 0.2, 0.3, 0.25, 3.5, "uzf02"], + [3, (0, 0, 4), 1, -1, 1.0, 1e-5, 0.2, 0.3, 0.25, 3.5, "uzf02"], + [4, (0, 0, 5), 1, -1, 1.0, 1e-5, 0.2, 0.3, 0.25, 3.5, "uzf03"], + [5, (0, 0, 5), 1, -1, 1.0, 1e-5, 0.2, 0.3, 0.25, 3.5, "uzf03"], + [6, (0, 0, 6), 1, -1, 1.0, 1e-5, 0.2, 0.3, 0.25, 3.5, "uzf04"], + [7, (0, 0, 6), 1, -1, 1.0, 1e-5, 0.2, 0.3, 0.25, 3.5, "uzf04"], + ], + "auxmultval": 0.5, +} + +nlay, nrow, ncol = 3, 1, 10 +nper = 1 +perlen = [1.0] +nstp = [1] +tsmult = [1.0] + +delr = 100.0 +delc = 100.0 +strt = -9.0 + +strt_temp = 1.0 +dispersivity = 1.0 +cpw = 4180.0 +rhow = 1000.0 +cps = 800.0 +rhos = 2500.0 +prsity = 0.2 + + +def build_gwf(gwf, gwfname): + flopy.mf6.ModflowGwfdis( + gwf, + nlay=nlay, + nrow=nrow, + ncol=ncol, + delr=delr, + delc=delc, + top=0.0, + botm=[-10.0, -20.0, -30.0], + ) + + # initial conditions + flopy.mf6.ModflowGwfic(gwf, strt=strt) + + # node property flow + flopy.mf6.ModflowGwfnpf( + gwf, + save_flows=True, + save_specific_discharge=True, + save_saturation=True, + icelltype=1, + k=1.0e-4, + k22=1.0e-4, + k33=1.0e-5, + ) + + # aquifer storage + flopy.mf6.ModflowGwfsto( + gwf, + save_flows=True, + iconvert=1, + ss=1e-5, + sy=prsity, + transient=True, + ) + + # chd files + chdval = -9.0 + iface = 0 + temperature = 10.0 + chdspd = { + 0: [ + [(0, 0, 0), chdval, iface, temperature, "object0"], + [(0, 0, ncol - 1), chdval, iface, temperature, "object0"], + ] + } + flopy.mf6.ModflowGwfchd( + gwf, + auxiliary=["iface", "temperature"], + boundnames=True, + print_input=True, + save_flows=True, + stress_period_data=chdspd, + pname="CHD", + ) + + # transient uzf info + # ifno cellid landflg ivertcn surfdp vks thtr thts thti eps [bndnm] + uz_pkdat = uzarea_data["uzf_pkdat"] + + finf = 1.0 + extdp = 0.0 + extwc = 0.0 + pet = 0.0 + zero = 0.0 + auxmultval = uzarea_data["auxmultval"] + uzf_spd = { + 0: [ + [i, finf, pet, extdp, extwc, zero, zero, zero, auxmultval] + for i in np.arange(len(uz_pkdat)) + ] + } + + flopy.mf6.ModflowGwfuzf( + gwf, + print_input=True, + print_flows=True, + save_flows=True, + boundnames=True, + ntrailwaves=7, + nwavesets=40, + auxiliary="multiplier", + auxmultname="multiplier", + package_convergence_filerecord=f"{gwfname}.UzfConvergence.csv", + wc_filerecord=f"{gwfname}.wc", + nuzfcells=len(uz_pkdat), + packagedata=uz_pkdat, + perioddata=uzf_spd, + budget_filerecord=f"{gwfname}.uzf.bud", + pname="UZF", + filename=f"{gwfname}.uzf", + ) + + # output control + flopy.mf6.ModflowGwfoc( + gwf, + budget_filerecord=f"{gwfname}.cbc", + head_filerecord=f"{gwfname}.hds", + headprintrecord=[("COLUMNS", 10, "WIDTH", 15, "DIGITS", 6, "GENERAL")], + saverecord=[("HEAD", "ALL"), ("BUDGET", "ALL")], + printrecord=[("HEAD", "ALL"), ("BUDGET", "ALL")], + filename=f"{gwfname}.oc", + ) + + return gwf + + +def build_gwe(gwe, gwename, fmi=False): + # Instantiating MODFLOW 6 transport discretization package + flopy.mf6.ModflowGwedis( + gwe, + nogrb=True, + nlay=nlay, + nrow=nrow, + ncol=ncol, + delr=delr, + delc=delc, + top=0.0, + botm=[-10.0, -20.0, -30.0], + pname="DIS", + filename=f"{gwename}.dis", + ) + + # Instantiating MODFLOW 6 transport initial concentrations + flopy.mf6.ModflowGweic( + gwe, + strt=strt_temp, + pname="IC", + filename=f"{gwename}.ic", + ) + + # Instantiating MODFLOW 6 transport advection package + flopy.mf6.ModflowGweadv(gwe, scheme="TVD", pname="ADV", filename=f"{gwename}.adv") + + # Instantiating MODFLOW 6 transport dispersion package + flopy.mf6.ModflowGwecnd( + gwe, + xt3d_off=False, + alh=dispersivity, + ath1=dispersivity, + ktw=0.5918 * 86400, + kts=0.2700 * 86400, + pname="CND", + filename=f"{gwename}.cnd", + ) + + # Instantiating MODFLOW 6 transport mass storage package + flopy.mf6.ModflowGweest( + gwe, + save_flows=True, + porosity=prsity, + heat_capacity_water=cpw, + density_water=rhow, + heat_capacity_solid=cps, + density_solid=rhos, + pname="EST", + filename=f"{gwename}.est", + ) + + # Instantiating MODFLOW 6 transport source-sink mixing package + srctype = "AUX" + auxname = "TEMPERATURE" + pname = ["CHD"] + # Inpput to SSM is: + sources = [[itm, srctype, auxname] for itm in pname] + + flopy.mf6.ModflowGwessm( + gwe, + sources=sources, + pname="SSM", + filename=f"{gwename}.ssm", + ) + + # Instantiating MODFLOW 6 energy transport source-sink mixing package + uz_pkdat = uzarea_data["uzf_pkdat"] + uzepackagedata = [(ct, 1.0) for ct, iuz in enumerate(uz_pkdat)] + uzeperioddata = { + 0: [[ct, "INFILTRATION", 1.0] for ct, itm in enumerate(uz_pkdat)], + } + + flopy.mf6.ModflowGweuze( + gwe, + flow_package_name="UZF", + boundnames=False, + save_flows=True, + print_input=True, + print_flows=True, + print_temperature=True, + temperature_filerecord=gwename + ".uze.bin", + budget_filerecord=gwename + ".uze.bud", + packagedata=uzepackagedata, + uzeperioddata=uzeperioddata, + pname="UZE", + filename=f"{gwename}.uze", + ) + + # Instantiate MODFLOW 6 heat transport output control package + flopy.mf6.ModflowGweoc( + gwe, + pname="OC", + budget_filerecord=f"{gwename}.cbc", + temperature_filerecord=f"{gwename}.ucn", + temperatureprintrecord=[("COLUMNS", 10, "WIDTH", 15, "DIGITS", 6, "GENERAL")], + saverecord=[("TEMPERATURE", "ALL"), ("BUDGET", "ALL")], + printrecord=[("TEMPERATURE", "ALL"), ("BUDGET", "ALL")], + filename=f"{gwename}.oc", + ) + + if fmi: + gwfname = gwename.replace("gwe", "gwf") + pd = [ + ("GWFHEAD", gwfname + ".hds", None), + ("GWFBUDGET", gwfname + ".cbc", None), + ("UZF", gwfname + ".uzf.bud", None), + ] + fmi = flopy.mf6.ModflowGwefmi(gwe, packagedata=pd) + + return gwe + + +def build_single_sim(idx, ws, exe): + name = cases[idx] + gwfname = "gwf-" + name + gwename = "gwe-" + name + + tdis_rc = [] + for i in range(nper): + tdis_rc.append((perlen[i], nstp[i], tsmult[i])) + + # build MODFLOW 6 files + sim = flopy.mf6.MFSimulation(sim_name=name, version="mf6", exe_name=exe, sim_ws=ws) + + # create tdis package + flopy.mf6.ModflowTdis(sim, time_units="DAYS", nper=nper, perioddata=tdis_rc) + + # create gwf model + gwf = flopy.mf6.ModflowGwf( + sim, modelname=gwfname, newtonoptions="NEWTON", save_flows=True + ) + + # build out gwf model + gwf = build_gwf(gwf, gwfname) + + # create iterative model solution and register the gwf model with it + ims = flopy.mf6.ModflowIms( + sim, + print_option="SUMMARY", + complexity="MODERATE", + outer_dvclose=hclose, + outer_maximum=nouter, + under_relaxation="DBD", + inner_maximum=ninner, + inner_dvclose=hclose, + rcloserecord=rclose, + linear_acceleration="BICGSTAB", + scaling_method="NONE", + reordering_method="NONE", + relaxation_factor=relax, + ) + sim.register_ims_package(ims, [gwf.name]) + + # ---------------------------------- + # Instantiating MODFLOW 6 GWE model + # ---------------------------------- + gwe = flopy.mf6.ModflowGwe(sim, modelname=gwename, model_nam_file=f"{gwename}.nam") + gwe.name_file.save_flows = True + + # build out gwe model + gwe = build_gwe(gwe, gwename) + + imsgwe = 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="BICGSTAB", + scaling_method="NONE", + reordering_method="NONE", + relaxation_factor=relax, + filename=f"{gwename}.ims", + ) + sim.register_ims_package(imsgwe, [gwe.name]) + + # Instantiate Gwf-Gwe Exchange package + flopy.mf6.ModflowGwfgwe( + sim, + exgtype="GWF6-GWE6", + exgmnamea=gwfname, + exgmnameb=gwename, + filename=f"{gwename}.gwfgwe", + ) + + return sim + + +def run_single_sim(dir, exe): + idx = 0 + sim = build_single_sim(idx, dir, exe) + sim.write_simulation() + success, buff = sim.run_simulation(silent=False) + errmsg = f"simulation should terminate with error message, but " f"did not.\n{buff}" + assert not success, errmsg + + +def build_gwf_sim_only(idx, dir, exe): + name = cases[idx] + gwfname = "gwf-" + name + + tdis_rc = [] + for i in range(nper): + tdis_rc.append((perlen[i], nstp[i], tsmult[i])) + + # build MODFLOW 6 files + sim = flopy.mf6.MFSimulation(sim_name=name, version="mf6", exe_name=exe, sim_ws=dir) + + # create tdis package + flopy.mf6.ModflowTdis(sim, time_units="DAYS", nper=nper, perioddata=tdis_rc) + + # create gwf model + gwf = flopy.mf6.ModflowGwf( + sim, modelname=gwfname, newtonoptions="NEWTON", save_flows=True + ) + + # build out gwf model + gwf = build_gwf(gwf, gwfname) + + # create iterative model solution and register the gwf model with it + ims = flopy.mf6.ModflowIms( + sim, + print_option="SUMMARY", + complexity="MODERATE", + outer_dvclose=hclose, + outer_maximum=nouter, + under_relaxation="DBD", + inner_maximum=ninner, + inner_dvclose=hclose, + rcloserecord=rclose, + linear_acceleration="BICGSTAB", + scaling_method="NONE", + reordering_method="NONE", + relaxation_factor=relax, + ) + sim.register_ims_package(ims, [gwf.name]) + + return sim + + +def build_gwe_sim_only(idx, dir, exe): + name = cases[idx] + gwename = "gwe-" + name + + tdis_rc = [] + for i in range(nper): + tdis_rc.append((perlen[i], nstp[i], tsmult[i])) + + # build MODFLOW 6 files + sim = flopy.mf6.MFSimulation(sim_name=name, version="mf6", exe_name=exe, sim_ws=dir) + + # create tdis package + flopy.mf6.ModflowTdis(sim, time_units="DAYS", nper=nper, perioddata=tdis_rc) + + # create gwf model + gwe = flopy.mf6.ModflowGwe( + sim, modelname=gwename, model_nam_file=f"{gwename}.nam", save_flows=True + ) + gwe = build_gwe(gwe, gwename, fmi=True) + + imsgwe = 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="BICGSTAB", + scaling_method="NONE", + reordering_method="NONE", + relaxation_factor=relax, + filename=f"{gwename}.ims", + ) + sim.register_ims_package(imsgwe, [gwe.name]) + + return sim + + +def run_separate_sims(dir, exe): + idx = 1 + sim = build_gwf_sim_only(idx, dir, exe) + sim.write_simulation() + success, buff = sim.run_simulation(silent=False) + + # So long as the flow simulation ran successfully, then run GWE' + errmsg = "GWF only simulation did not run and should have" + assert success, errmsg + sim = build_gwe_sim_only(idx, dir, exe) + sim.write_simulation() + success, buff = sim.run_simulation(silent=False) + errmsg = ( + f"GWE simulation should terminate with error message, but " f"did not.\n{buff}" + ) + assert not success, errmsg + + +def test_auxmultname(function_tmpdir, targets): + # Whether running the GWF & GWE model in a single sim or + # in separate sims via FMI, both should error out, which is + # how the assertion statements are constructed. + # Start with check where both GWF and GWE are in same simulation + run_single_sim(str(function_tmpdir), targets["mf6"]) + + # Next, separate GWF and GWE models into 2 different simulations + # and ensure MF quits with error message + run_separate_sims(str(function_tmpdir), targets["mf6"]) diff --git a/doc/ReleaseNotes/develop.tex b/doc/ReleaseNotes/develop.tex index b82e6075d56..397dd3614e9 100644 --- a/doc/ReleaseNotes/develop.tex +++ b/doc/ReleaseNotes/develop.tex @@ -47,6 +47,7 @@ \item The GWF-GWF Exchange has been fixed to support the configuration of the Mover package (MVR) also in cases where the number of exchanges equals zero (NEXG = 0). In earlier versions this has caused MODFLOW to terminate with an error. \item A PRT bug was fixed in which array-based input read from the RCH (recharge) or EVT (evapotranspiration) packages could fail to be processed correctly by the PRT FMI (flow model interface) package, causing a crash. \item When the SQUARE\_GWET option was invoked in the UZF options block, evapotranspiration from the water table (GWET) was calculated incorrectly. Instead of acting as a sink, the calculated evapotranspiration flux was added as a source of water. The applied fix ensures that groundwater evapotranspiration is removed from the water table and as a result the GWET values are accumulated as outflows in the budget table. + \item The UZE package will instantiate one UZE object for each UZF object specified in the flow model. UZF will support specification of multiple UZF objects within a given cell; however, GWE will not work properly where multiple UZE objects exist in a cell. Code was added to check for multiple UZE objects in a grid cell and will stop the simulation with an appropriate message when this condition is discovered. % \item xxx \end{itemize} diff --git a/doc/mf6io/gwe/uze.tex b/doc/mf6io/gwe/uze.tex index 1e65d5bbd11..255348d9a28 100644 --- a/doc/mf6io/gwe/uze.tex +++ b/doc/mf6io/gwe/uze.tex @@ -1,6 +1,6 @@ Unsaturated Zone Energy Transport (UZE) Package information is read from the file that is specified by ``UZE6'' as the file type. There can be as many UZE Packages as necessary for a GWE model. Each UZE Package is designed to work with flows from a corresponding GWF UZF Package. By default \mf uses the UZE package name to determine which UZF Package corresponds to the UZE Package. Therefore, the package name of the UZE Package (as specified in the GWE name file) must match with the name of the corresponding UZF Package (as specified in the GWF name file). Alternatively, the name of the flow package can be specified using the FLOW\_PACKAGE\_NAME keyword in the options block. The GWE UZE Package cannot be used without a corresponding GWF UZF Package. -The UZE Package does not have a dimensions block; instead, dimensions for the UZE Package are set using the dimensions from the corresponding UZF Package. For example, the UZF Package requires specification of the number of cells (NUZFCELLS). UZE sets the number of UZE cells equal to NUZFCELLS. Therefore, the PACKAGEDATA block below must have NUZFCELLS entries in it. +The UZE Package does not have a dimensions block; instead, dimensions for the UZE Package are set using the dimensions from the corresponding UZF Package. For example, the UZF Package requires specification of the number of cells (NUZFCELLS). UZE sets the number of UZE cells equal to NUZFCELLS. Therefore, the PACKAGEDATA block below must have NUZFCELLS entries in it. Because of how UZE is formulated, only one UZE object per grid cell is permitted. If more than one UZE object is discovered in a grid cell, GWE will exit with an error message. To disable multiple UZE objects in a grid cell, remove the AUXMULTNAME keyword from the OPTIONS block of the corresponding UZF input file. \vspace{5mm} \subsubsection{Structure of Blocks} diff --git a/src/Model/GroundWaterEnergy/gwe-uze.f90 b/src/Model/GroundWaterEnergy/gwe-uze.f90 index a6bf10dce25..f7e2058dd65 100644 --- a/src/Model/GroundWaterEnergy/gwe-uze.f90 +++ b/src/Model/GroundWaterEnergy/gwe-uze.f90 @@ -82,6 +82,8 @@ module GweUzeModule procedure :: bnd_ac => uze_ac procedure :: bnd_mc => uze_mc procedure :: get_mvr_depvar + procedure, private :: stop_if_multi_uzobj_in_cell + procedure, private :: too_many_uzobjs end type GweUzeType @@ -156,22 +158,28 @@ subroutine find_uze_package(this) ! -- dummy class(GweUzeType) :: this ! -- local - character(len=LINELENGTH) :: errmsg class(BndType), pointer :: packobj - integer(I4B) :: ip, icount + integer(I4B) :: ip, icount, nuz integer(I4B) :: nbudterm + character(len=LINELENGTH) :: errmsg logical :: found ! ! -- Initialize found to false, and error later if flow package cannot ! be found found = .false. ! + nuz = 1 + ! ! -- If user is specifying flows in a binary budget file, then set up ! the budget file reader, otherwise set a pointer to the flow package ! budobj if (this%fmi%flows_from_file) then call this%fmi%set_aptbudobj_pointer(this%flowpackagename, this%flowbudptr) - if (associated(this%flowbudptr)) found = .true. + if (associated(this%flowbudptr)) then + found = .true. + nuz = this%flowbudptr%budterm(nuz)%maxlist + call this%stop_if_multi_uzobj_in_cell(nuz) + end if ! else if (associated(this%fmi%gwfbndlist)) then @@ -190,6 +198,11 @@ subroutine find_uze_package(this) this%flowbudptr => packobj%budobj end select end if + ! + ! -- Error if option for multiple UZF objects per cell is in use + if (packobj%iauxmultcol > 0) then + call this%too_many_uzobjs() + end if if (found) exit end do end if @@ -1293,6 +1306,55 @@ subroutine uze_bd_obs(this, obstypeid, jj, v, found) end select end subroutine uze_bd_obs + !> @brief Stop simulation when multiple UZF objects discovered in cell + !! + !! GWE equilibrates the temperature of a UZE object with the host cell. + !! This is accomplished through the use of a residual term in the energy + !! balance of each cell. As a result, GWE, specifically UZE, does not + !! support multiple UZE objects within a given cell, requiring the code + !! to exit with an appropriate message. + !< + subroutine stop_if_multi_uzobj_in_cell(this, nuz) + ! -- modules + use ConstantsModule, only: IZERO + ! -- dummy + class(GweUzeType) :: this + integer(I4B) :: nuz + ! -- local + integer(I4B) :: n, iloc + integer(I4B), dimension(:), allocatable :: dup_chk + character(len=LINELENGTH) :: errmsg + ! + allocate (dup_chk(nuz)) + do n = 1, nuz + dup_chk(n) = IZERO + end do + ! + ! -- add explanation + do n = 1, nuz + if (this%flowbudptr%budterm(1)%id2(n) /= IZERO) then + call this%too_many_uzobjs() + end if + iloc = this%flowbudptr%budterm(1)%id2(n) + dup_chk(iloc) = 1 + end do + end subroutine stop_if_multi_uzobj_in_cell + + !> @brief Print and store error msg for too many UZF/UZE objs in cell + !< + subroutine too_many_uzobjs(this) + ! -- dummy + class(GweUzeType) :: this + ! -- local + character(len=LINELENGTH) :: errmsg + ! + write (errmsg, '(a)') 'UZE does not support the use of multiple UZF & + &objects in a cell. Check use of AUXMULTNAME & + &option in UZF package.' + call store_error(errmsg) + call this%parser%StoreErrorUnit() + end subroutine too_many_uzobjs + !> @brief Sets the stress period attributes for keyword use. !< subroutine uze_set_stressperiod(this, itemno, keyword, found) diff --git a/src/Model/GroundWaterFlow/gwf-uzf.f90 b/src/Model/GroundWaterFlow/gwf-uzf.f90 index 33f7ab1cd1e..d8aedabd490 100644 --- a/src/Model/GroundWaterFlow/gwf-uzf.f90 +++ b/src/Model/GroundWaterFlow/gwf-uzf.f90 @@ -2032,7 +2032,7 @@ subroutine read_cell_properties(this) end do ! ! -- do an initial evaluation of the sum of uzfarea relative to the - ! GWF cell area in the case that there is more than one UZF cell + ! GWF cell area in the case that there is more than one UZF object ! in a GWF cell and a auxmult value is not being applied to the ! calculate the UZF cell area from the GWF cell area. if (this%imaxcellcnt > 1 .and. this%iauxmultcol < 1) then