Skip to content

Commit

Permalink
fix(par): bug with print_input in exchange (#1561)
Browse files Browse the repository at this point in the history
- fix bug with exchanges and `print_input` in parallel mode
- add test for `print_input`
- add test for parallel disu coupling
  • Loading branch information
mjr-deltares authored Jan 19, 2024
1 parent 69f2713 commit 2100373
Show file tree
Hide file tree
Showing 7 changed files with 278 additions and 34 deletions.
1 change: 1 addition & 0 deletions autotest/test_par_gwf01.py
Original file line number Diff line number Diff line change
Expand Up @@ -190,6 +190,7 @@ def get_model(idx, dir):
exgmnameb=name_right,
exchangedata=gwfgwf_data,
auxiliary=["ANGLDEGX", "CDIST"],
print_input=True,
)

return sim
Expand Down
238 changes: 238 additions & 0 deletions autotest/test_par_gwf_disu.py
Original file line number Diff line number Diff line change
@@ -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()
8 changes: 8 additions & 0 deletions src/Distributed/VirtualModel.f90
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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

Expand All @@ -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/))
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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

Expand Down
25 changes: 14 additions & 11 deletions src/Exchange/DisConnExchange.f90
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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)
Expand Down Expand Up @@ -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
Expand All @@ -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'
!
Expand All @@ -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
!
Expand All @@ -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
!
Expand All @@ -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), &
Expand All @@ -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 ('// &
Expand All @@ -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 ('// &
Expand Down
Loading

0 comments on commit 2100373

Please sign in to comment.