diff --git a/autotest/test_par_gwf01.py b/autotest/test_par_gwf01.py index 57c38874892..a3c460baf4b 100644 --- a/autotest/test_par_gwf01.py +++ b/autotest/test_par_gwf01.py @@ -190,6 +190,7 @@ def get_model(idx, dir): exgmnameb=name_right, exchangedata=gwfgwf_data, auxiliary=["ANGLDEGX", "CDIST"], + print_input=True, ) return sim diff --git a/autotest/test_par_gwf_disu.py b/autotest/test_par_gwf_disu.py new file mode 100644 index 00000000000..da77695e6c2 --- /dev/null +++ b/autotest/test_par_gwf_disu.py @@ -0,0 +1,238 @@ +""" +Test for parallel MODFLOW running on two cpus. +Each case contains two coupled DISU models which are +generated from their DIS counterparts: + +1d: (nlay,nrow,ncol) = (1,1,5), +2d: (nlay,nrow,ncol) = (1,5,5), +3d: (nlay,nrow,ncol) = (5,5,5), + +constant head boundaries left=1.0, right=10.0. +The result should be a uniform flow field. +""" + +import os + +import flopy +import numpy as np +import pytest +from flopy.utils.gridutil import get_disu_kwargs + +from framework import TestFramework + +cases = [ + "par_gwf_disu1d", + "par_gwf_disu2d", + "par_gwf_disu3d", +] +dis_shape = [(1, 1, 5), (1, 1, 5), (1, 1, 5), (1, 5, 5), (5, 5, 5)] + +# global convenience... +name_left = "leftmodel" +name_right = "rightmodel" + +# solver data +nouter, ninner = 100, 300 +hclose, rclose, relax = 10e-9, 1e-3, 0.97 + + +def get_model(idx, dir): + name = cases[idx] + + # parameters and spd + # tdis + nper = 1 + tdis_rc = [] + for i in range(nper): + tdis_rc.append((1.0, 1, 1)) + + # model spatial discretization + nlay = dis_shape[idx][0] + nrow = dis_shape[idx][1] + ncol = dis_shape[idx][2] + + # cell spacing + delr = 100.0 + delr_arr = delr * np.ones(ncol) + delc = 100.0 + delc_arr = delc * np.ones(nrow) + + # shift + shift_x = 5 * delr + shift_y = 0.0 + + # top/bot of the aquifer + tops = [0.0, -100.0, -200.0, -300.0, -400.0, -500.0] + + # conversion to DISU + disukwargs = get_disu_kwargs( + nlay, + nrow, + ncol, + delr_arr, + delc_arr, + tops[0], + tops[1 : nlay + 1], + return_vertices=True, + ) + + # hydraulic conductivity + k11 = 1.0 + + # boundary stress period data + h_left = 1.0 + h_right = 10.0 + + # initial head + h_start = 0.0 + + sim = flopy.mf6.MFSimulation( + sim_name=name, + version="mf6", + exe_name="mf6", + sim_ws=dir, + ) + + tdis = flopy.mf6.ModflowTdis( + sim, time_units="DAYS", nper=nper, perioddata=tdis_rc + ) + + ims = flopy.mf6.ModflowIms( + sim, + print_option="ALL", + outer_dvclose=hclose, + outer_maximum=nouter, + inner_maximum=ninner, + inner_dvclose=hclose, + rcloserecord=rclose, + linear_acceleration="BICGSTAB", + relaxation_factor=relax, + pname="ims", + ) + + # submodel on the left: + left_chd = [ + [(ilay * ncol * nrow + irow * ncol), h_left] + for irow in range(nrow) + for ilay in range(nlay) + ] + chd_spd_left = {0: left_chd} + + gwf = flopy.mf6.ModflowGwf(sim, modelname=name_left, save_flows=True) + disu = flopy.mf6.ModflowGwfdisu(gwf, **disukwargs) + ic = flopy.mf6.ModflowGwfic(gwf, strt=h_start) + npf = flopy.mf6.ModflowGwfnpf( + gwf, + save_specific_discharge=False, # let's skip angledegx + save_flows=True, + icelltype=0, + k=k11, + ) + chd = flopy.mf6.ModflowGwfchd(gwf, stress_period_data=chd_spd_left) + oc = flopy.mf6.ModflowGwfoc( + gwf, + head_filerecord=f"{name_left}.hds", + budget_filerecord=f"{name_left}.cbc", + headprintrecord=[("COLUMNS", 10, "WIDTH", 15, "DIGITS", 6, "GENERAL")], + saverecord=[("HEAD", "LAST"), ("BUDGET", "LAST")], + ) + + # submodel on the right: + right_chd = [ + [(ilay * ncol * nrow + irow * ncol + ncol - 1), h_right] + for irow in range(nrow) + for ilay in range(nlay) + ] + chd_spd_right = {0: right_chd} + + gwf = flopy.mf6.ModflowGwf(sim, modelname=name_right, save_flows=True) + disukwargs["xorigin"] = shift_x + disukwargs["yorigin"] = shift_y + disu = flopy.mf6.ModflowGwfdisu(gwf, **disukwargs) + ic = flopy.mf6.ModflowGwfic(gwf, strt=h_start) + npf = flopy.mf6.ModflowGwfnpf( + gwf, + save_specific_discharge=False, # let's skip angledegx + save_flows=True, + icelltype=0, + k=k11, + ) + chd = flopy.mf6.ModflowGwfchd(gwf, stress_period_data=chd_spd_right) + oc = flopy.mf6.ModflowGwfoc( + gwf, + head_filerecord=f"{name_right}.hds", + budget_filerecord=f"{name_right}.cbc", + headprintrecord=[("COLUMNS", 10, "WIDTH", 15, "DIGITS", 6, "GENERAL")], + saverecord=[("HEAD", "LAST"), ("BUDGET", "LAST")], + ) + + # exchangedata + angldegx = 0.0 + cdist = delr + gwfgwf_data = [ + [ + (ilay * ncol * nrow + irow * ncol + ncol - 1), + (ilay * ncol * nrow + irow * ncol), + 1, + delr / 2.0, + delr / 2.0, + delc, + angldegx, + cdist, + ] + for irow in range(nrow) + for ilay in range(nlay) + ] + + gwfgwf = flopy.mf6.ModflowGwfgwf( + sim, + exgtype="GWF6-GWF6", + nexg=len(gwfgwf_data), + exgmnamea=name_left, + exgmnameb=name_right, + exchangedata=gwfgwf_data, + auxiliary=["ANGLDEGX", "CDIST"], + print_input=True, + ) + + return sim + + +def build_models(idx, test): + sim = get_model(idx, test.workspace) + return sim, None + + +def check_output(idx, test): + pass + # two coupled models with a uniform flow field, + # here we assert the known head values at the + # cell centers + fpth = os.path.join(test.workspace, f"{name_left}.hds") + hds = flopy.utils.HeadFile(fpth) + heads_left = hds.get_data().flatten() + fpth = os.path.join(test.workspace, f"{name_right}.hds") + hds = flopy.utils.HeadFile(fpth) + heads_right = hds.get_data().flatten() + np.testing.assert_array_almost_equal( + heads_left[0:5], [1.0, 2.0, 3.0, 4.0, 5.0] + ) + np.testing.assert_array_almost_equal( + heads_right[0:5], [6.0, 7.0, 8.0, 9.0, 10.0] + ) + + +@pytest.mark.parallel +@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), + compare=None, + parallel=True, + ncpus=2, + ) + test.run() diff --git a/src/Distributed/VirtualModel.f90 b/src/Distributed/VirtualModel.f90 index be66b02b794..bdb2a30fdf4 100644 --- a/src/Distributed/VirtualModel.f90 +++ b/src/Distributed/VirtualModel.f90 @@ -45,6 +45,8 @@ module VirtualModelModule type(VirtualDbl1dType), pointer :: x => null() type(VirtualDbl1dType), pointer :: x_old => null() type(VirtualInt1dType), pointer :: ibound => null() + ! Base Model fields + type(VirtualIntType), pointer :: idsoln => null() contains ! public procedure :: create => vm_create @@ -116,6 +118,8 @@ subroutine init_virtual_data(this) call this%set(this%x%base(), 'X', '', MAP_NODE_TYPE) call this%set(this%x_old%base(), 'XOLD', '', MAP_NODE_TYPE) call this%set(this%ibound%base(), 'IBOUND', '', MAP_NODE_TYPE) + ! Base model + call this%set(this%idsoln%base(), 'IDSOLN', '', MAP_ALL_TYPE) end subroutine init_virtual_data @@ -128,6 +132,7 @@ subroutine vm_prepare_stage(this, stage) if (stage == STG_AFT_MDL_DF) then + call this%map(this%idsoln%base(), (/STG_AFT_MDL_DF/)) call this%map(this%con_ianglex%base(), (/STG_AFT_MDL_DF/)) call this%map(this%dis_ndim%base(), (/STG_AFT_MDL_DF/)) call this%map(this%dis_nodes%base(), (/STG_AFT_MDL_DF/)) @@ -250,6 +255,7 @@ subroutine allocate_data(this) allocate (this%x) allocate (this%x_old) allocate (this%ibound) + allocate (this%idsoln) end subroutine allocate_data @@ -285,6 +291,8 @@ subroutine deallocate_data(this) deallocate (this%x) deallocate (this%x_old) deallocate (this%ibound) + ! Base model + deallocate (this%idsoln) end subroutine deallocate_data diff --git a/src/Exchange/DisConnExchange.f90 b/src/Exchange/DisConnExchange.f90 index 760e0f32f41..3429c5bf963 100644 --- a/src/Exchange/DisConnExchange.f90 +++ b/src/Exchange/DisConnExchange.f90 @@ -194,8 +194,8 @@ subroutine source_dimensions(this, iout) return end subroutine source_dimensions - !> @brief - !< + !> @brief Returns reduced node number from user + !< specified cell id. function noder(this, model, cellid, iout) ! -- modules use GeomUtilModule, only: get_node @@ -226,11 +226,11 @@ end function noder !> @brief !< - function cellstr(this, model, cellid, iout) + function cellstr(this, ndim, cellid, iout) ! -- modules ! -- dummy class(DisConnExchangeType) :: this !< instance of exchange object - class(NumericalModelType), pointer, intent(in) :: model + integer(I4B) :: ndim !< model DIS dimension integer(I4B), dimension(:), pointer, intent(in) :: cellid integer(I4B), intent(in) :: iout !< the output file unit character(len=20) :: cellstr @@ -243,7 +243,7 @@ function cellstr(this, model, cellid, iout) ! cellstr = '' ! - select case (model%dis%ndim) + select case (ndim) case (1) write (cellstr, fmtndim1) cellid(1) case (2) @@ -274,6 +274,7 @@ subroutine source_data(this, iout) real(DP), dimension(:), contiguous, pointer :: hwva real(DP), dimension(:, :), contiguous, pointer :: auxvar type(CharacterStringType), dimension(:), contiguous, pointer :: boundname + integer(I4B) :: ndim1, ndim2 character(len=20) :: cellstr1, cellstr2 character(len=2) :: cnfloat integer(I4B) :: nerr, iaux @@ -292,6 +293,8 @@ subroutine source_data(this, iout) call mem_setptr(hwva, 'HWVA', this%input_mempath) call mem_setptr(auxvar, 'AUX', this%input_mempath) call mem_setptr(boundname, 'BOUNDNAME', this%input_mempath) + ndim1 = size(cellidm1, dim=1) + ndim2 = size(cellidm2, dim=1) ! write (iout, '(1x,a)') 'PROCESSING EXCHANGEDATA' ! @@ -316,7 +319,7 @@ subroutine source_data(this, iout) ! if (associated(this%model1)) then ! - ! -- Determine user node number + ! -- Determine reduced node number nodem1 = this%noder(this%model1, cellidm1(:, iexg), iout) this%nodem1(iexg) = nodem1 ! @@ -326,7 +329,7 @@ subroutine source_data(this, iout) ! if (associated(this%model2)) then ! - ! -- Determine user node number + ! -- Determine reduced node number nodem2 = this%noder(this%model2, cellidm2(:, iexg), iout) this%nodem2(iexg) = nodem2 ! @@ -348,8 +351,8 @@ subroutine source_data(this, iout) ! ! -- Write the data to listing file if requested if (this%iprpak /= 0) then - cellstr1 = this%cellstr(this%model1, cellidm1(:, iexg), iout) - cellstr2 = this%cellstr(this%model2, cellidm2(:, iexg), iout) + cellstr1 = this%cellstr(ndim1, cellidm1(:, iexg), iout) + cellstr2 = this%cellstr(ndim2, cellidm2(:, iexg), iout) if (this%inamedbound == 0) then write (iout, fmtexgdata) trim(cellstr1), trim(cellstr2), & this%ihc(iexg), this%cl1(iexg), this%cl2(iexg), & @@ -367,7 +370,7 @@ subroutine source_data(this, iout) ! -- Check to see if nodem1 is outside of active domain if (associated(this%model1)) then if (nodem1 <= 0) then - cellstr1 = this%cellstr(this%model1, cellidm1(:, iexg), iout) + cellstr1 = this%cellstr(ndim1, cellidm1(:, iexg), iout) write (errmsg, *) & trim(adjustl(this%model1%name))// & ' Cell is outside active grid domain ('// & @@ -379,7 +382,7 @@ subroutine source_data(this, iout) ! -- Check to see if nodem2 is outside of active domain if (associated(this%model2)) then if (nodem2 <= 0) then - cellstr2 = this%cellstr(this%model2, cellidm2(:, iexg), iout) + cellstr2 = this%cellstr(ndim2, cellidm2(:, iexg), iout) write (errmsg, *) & trim(adjustl(this%model2%name))// & ' Cell is outside active grid domain ('// & diff --git a/src/Exchange/GwfGwfExchange.f90 b/src/Exchange/GwfGwfExchange.f90 index e48b7be02dc..4b933374c53 100644 --- a/src/Exchange/GwfGwfExchange.f90 +++ b/src/Exchange/GwfGwfExchange.f90 @@ -218,15 +218,13 @@ subroutine gwf_gwf_df(this) write (iout, '(/a,a)') ' Creating exchange: ', this%name ! ! -- Ensure models are in same solution - if (associated(this%gwfmodel1) .and. associated(this%gwfmodel2)) then - if (this%gwfmodel1%idsoln /= this%gwfmodel2%idsoln) then - call store_error('Two models are connected in a GWF '// & - 'exchange but they are in different solutions. '// & - 'GWF models must be in same solution: '// & - trim(this%gwfmodel1%name)//' '// & - trim(this%gwfmodel2%name)) - call store_error_filename(this%filename) - end if + if (this%v_model1%idsoln%get() /= this%v_model2%idsoln%get()) then + call store_error('Two models are connected in a GWF '// & + 'exchange but they are in different solutions. '// & + 'GWF models must be in same solution: '// & + trim(this%v_model1%name)//' '// & + trim(this%v_model2%name)) + call store_error_filename(this%filename) end if ! ! -- source options @@ -282,7 +280,7 @@ subroutine validate_exchange(this) logical(LGP) :: has_k22, has_spdis, has_vsc ! ! Periodic boundary condition in exchange don't allow XT3D (=interface model) - if (associated(this%model1, this%model2)) then + if (this%v_model1 == this%v_model2) then if (this%ixt3d > 0) then write (errmsg, '(3a)') 'GWF-GWF exchange ', trim(this%name), & ' is a periodic boundary condition which cannot'// & diff --git a/src/Exchange/GwtGwtExchange.f90 b/src/Exchange/GwtGwtExchange.f90 index 220294a96a4..c4833aeef7d 100644 --- a/src/Exchange/GwtGwtExchange.f90 +++ b/src/Exchange/GwtGwtExchange.f90 @@ -207,15 +207,13 @@ subroutine gwt_gwt_df(this) write (iout, '(/a,a)') ' Creating exchange: ', this%name ! ! -- Ensure models are in same solution - if (associated(this%gwtmodel1) .and. associated(this%gwtmodel2)) then - if (this%gwtmodel1%idsoln /= this%gwtmodel2%idsoln) then - call store_error('Two models are connected in a GWT '// & - 'exchange but they are in different solutions. '// & - 'GWT models must be in same solution: '// & - trim(this%gwtmodel1%name)//' '// & - trim(this%gwtmodel2%name)) - call store_error_filename(this%filename) - end if + if (this%v_model1%idsoln%get() /= this%v_model2%idsoln%get()) then + call store_error('Two models are connected in a GWT '// & + 'exchange but they are in different solutions. '// & + 'GWT models must be in same solution: '// & + trim(this%v_model1%name)//' '// & + trim(this%v_model2%name)) + call store_error_filename(this%filename) end if ! ! -- source options @@ -270,7 +268,7 @@ subroutine validate_exchange(this) end if ! ! Periodic boundary condition in exchange don't allow XT3D (=interface model) - if (associated(this%model1, this%model2)) then + if (this%v_model1 == this%v_model2) then if (this%ixt3d > 0) then write (errmsg, '(3a)') 'GWT-GWT exchange ', trim(this%name), & ' is a periodic boundary condition which cannot'// & diff --git a/src/Model/Connection/GridConnection.f90 b/src/Model/Connection/GridConnection.f90 index ae4b45e075f..03a61f5ae3d 100644 --- a/src/Model/Connection/GridConnection.f90 +++ b/src/Model/Connection/GridConnection.f90 @@ -1056,9 +1056,7 @@ subroutine buildInterfaceMap(this) ! first get the participating models call model_ids%init() do i = 1, this%nrOfCells - if (.not. model_ids%contains(this%idxToGlobal(i)%v_model%id)) then - call model_ids%push_back(this%idxToGlobal(i)%v_model%id) - end if + call model_ids%push_back_unique(this%idxToGlobal(i)%v_model%id) end do ! initialize the map