fix(par): specific discharge for parallel simulations (#2479)

* - add spdis check to virtual gwe and gwt model: only synchronize spdis when it is available

* - python format

* - add check and error for index map out of bounds

* - add release note

* - fixed typo in note

* - forgot about map being present but empty

* - remove commented code

* - check: map can be empty
- fix: map is zero-based

Author: mjr-deltares

autotest/test_gwt_adv03_gwtgwt.py

@@ -0,0 +1,317 @@
+"""
+Test gwt-gwt coupling for advection without dispersion.
+"""
+
+import flopy
+import pytest
+from framework import TestFramework
+
+cases = ["adv03_gwtgwt"]
+scheme = "tvd"
+
+# solver settings
+nouter, ninner = 100, 300
+hclose, rclose, relax = 1e-6, 1e-6, 1.0
+
+g_delr = 1.0
+
+
+def get_gwf_model(sim, gwfname, gwfpath, modelshape, chdspd=None, welspd=None):
+    nlay, nrow, ncol, xshift, yshift = modelshape
+    delr = g_delr
+    delc = 1.0
+    top = 1.0
+    botm = [0.0]
+    strt = 1.0
+    hk = 1.0
+    laytyp = 0
+
+    gwf = flopy.mf6.ModflowGwf(
+        sim,
+        modelname=gwfname,
+        save_flows=True,
+    )
+
+    dis = flopy.mf6.ModflowGwfdis(
+        gwf,
+        nlay=nlay,
+        nrow=nrow,
+        ncol=ncol,
+        delr=delr,
+        delc=delc,
+        top=top,
+        botm=botm,
+        xorigin=xshift,
+        yorigin=yshift,
+    )
+
+    # initial conditions
+    ic = flopy.mf6.ModflowGwfic(gwf, strt=strt)
+
+    # node property flow
+    npf = flopy.mf6.ModflowGwfnpf(
+        gwf,
+        icelltype=laytyp,
+        k=hk,
+        save_specific_discharge=False,
+    )
+
+    # chd files
+    if chdspd is not None:
+        chd = flopy.mf6.modflow.mfgwfchd.ModflowGwfchd(
+            gwf,
+            stress_period_data=chdspd,
+            save_flows=False,
+            pname="CHD-1",
+        )
+
+    # wel files
+    if welspd is not None:
+        wel = flopy.mf6.ModflowGwfwel(
+            gwf,
+            print_input=True,
+            print_flows=True,
+            stress_period_data=welspd,
+            save_flows=False,
+            auxiliary="CONCENTRATION",
+            pname="WEL-1",
+        )
+
+    # output control
+    oc = flopy.mf6.ModflowGwfoc(
+        gwf,
+        budget_filerecord=f"{gwfname}.cbc",
+        head_filerecord=f"{gwfname}.hds",
+        headprintrecord=[("COLUMNS", 10, "WIDTH", 15, "DIGITS", 6, "GENERAL")],
+        saverecord=[("HEAD", "LAST"), ("BUDGET", "LAST")],
+        printrecord=[("HEAD", "LAST"), ("BUDGET", "LAST")],
+    )
+    gwf.set_model_relative_path(gwfpath)
+    return gwf
+
+
+def get_gwt_model(sim, gwtname, gwtpath, modelshape, scheme, sourcerecarray=None):
+    nlay, nrow, ncol, xshift, yshift = modelshape
+    delr = g_delr
+    delc = 1.0
+    top = 1.0
+    botm = [0.0]
+    strt = 1.0
+    hk = 1.0
+    laytyp = 0
+
+    gwt = flopy.mf6.MFModel(
+        sim,
+        model_type="gwt6",
+        modelname=gwtname,
+    )
+    gwt.name_file.save_flows = True
+
+    dis = flopy.mf6.ModflowGwtdis(
+        gwt,
+        nlay=nlay,
+        nrow=nrow,
+        ncol=ncol,
+        delr=delr,
+        delc=delc,
+        top=top,
+        botm=botm,
+        xorigin=xshift,
+        yorigin=yshift,
+    )
+
+    # initial conditions
+    ic = flopy.mf6.ModflowGwtic(gwt, strt=0.0)
+
+    # advection
+    adv = flopy.mf6.ModflowGwtadv(gwt, scheme=scheme)
+
+    # mass storage and transfer
+    mst = flopy.mf6.ModflowGwtmst(gwt, porosity=0.1)
+
+    # sources
+    ssm = flopy.mf6.ModflowGwtssm(gwt, sources=sourcerecarray)
+
+    # output control
+    oc = flopy.mf6.ModflowGwtoc(
+        gwt,
+        budget_filerecord=f"{gwtname}.cbc",
+        concentration_filerecord=f"{gwtname}.ucn",
+        concentrationprintrecord=[("COLUMNS", 10, "WIDTH", 15, "DIGITS", 6, "GENERAL")],
+        saverecord=[("CONCENTRATION", "LAST"), ("BUDGET", "LAST")],
+        printrecord=[("CONCENTRATION", "LAST"), ("BUDGET", "LAST")],
+    )
+
+    gwt.set_model_relative_path(gwtpath)
+    return gwt
+
+
+def build_models(idx, test):
+    # temporal discretization
+    nper = 1
+    perlen = [5.0]
+    nstp = [20]
+    tsmult = [1.0]
+    tdis_rc = []
+    for i in range(nper):
+        tdis_rc.append((perlen[i], nstp[i], tsmult[i]))
+
+    # build MODFLOW 6 files
+    ws = test.workspace
+    sim = flopy.mf6.MFSimulation(sim_name=ws, version="mf6", exe_name="mf6", sim_ws=ws)
+    # create tdis package
+    tdis = flopy.mf6.ModflowTdis(
+        sim, time_units="DAYS", nper=nper, perioddata=tdis_rc, pname="sim.tdis"
+    )
+
+    # grid information
+    nlay, nrow, ncol = 1, 2, 50
+
+    # Create gwf1 model
+    welspd = {0: [[(0, 0, 0), 1.0, 1.0], [(0, 1, 0), 1.0, 1.0]]}
+    chdspd = None
+    gwf1 = get_gwf_model(
+        sim,
+        "flow1",
+        "flow1",
+        (nlay, nrow, ncol, 0.0, 0.0),
+        chdspd=chdspd,
+        welspd=welspd,
+    )
+
+    # Create gwf2 model
+    welspd = None
+    chdspd = {0: [[(0, 0, ncol - 1), 0.0000000], [(0, 1, ncol - 1), 0.0000000]]}
+    gwf2 = get_gwf_model(
+        sim,
+        "flow2",
+        "flow2",
+        (nlay, nrow, ncol, ncol * g_delr, 0.0),
+        chdspd=chdspd,
+        welspd=welspd,
+    )
+
+    # gwf-gwf with interface model enabled
+    gwfgwf_data = [
+        [(0, 0, ncol - 1), (0, 0, 0), 1, 0.5, 0.5, 1.0, 0.0, 1.0],
+        [(0, 1, ncol - 1), (0, 1, 0), 1, 0.5, 0.5, 1.0, 0.0, 1.0],
+    ]
+    gwfgwf = flopy.mf6.ModflowGwfgwf(
+        sim,
+        exgtype="GWF6-GWF6",
+        nexg=len(gwfgwf_data),
+        exgmnamea=gwf1.name,
+        exgmnameb=gwf2.name,
+        exchangedata=gwfgwf_data,
+        auxiliary=["ANGLDEGX", "CDIST"],
+        filename="flow1_flow2.gwfgwf",
+    )
+
+    # create iterative model solution and register the gwf model with it
+    imsgwf = 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="CG",
+        scaling_method="NONE",
+        reordering_method="NONE",
+        relaxation_factor=relax,
+        filename="flow.ims",
+    )
+    sim.register_ims_package(imsgwf, [gwf1.name, gwf2.name])
+
+    # Create gwt model
+    sourcerecarray = [("WEL-1", "AUX", "CONCENTRATION")]
+    gwt1 = get_gwt_model(
+        sim,
+        "transport1",
+        "transport1",
+        (nlay, nrow, ncol, 0.0, 0.0),
+        scheme,
+        sourcerecarray=sourcerecarray,
+    )
+
+    # Create gwt model
+    sourcerecarray = None
+    gwt2 = get_gwt_model(
+        sim,
+        "transport2",
+        "transport2",
+        (nlay, nrow, ncol, ncol * g_delr, 0.0),
+        scheme,
+        sourcerecarray=sourcerecarray,
+    )
+
+    # Create GWT GWT exchange
+    gwt1gwt2 = flopy.mf6.ModflowGwtgwt(
+        sim,
+        exgtype="GWT6-GWT6",
+        gwfmodelname1=gwf1.name,
+        gwfmodelname2=gwf2.name,
+        adv_scheme=scheme,
+        nexg=len(gwfgwf_data),
+        exgmnamea=gwt1.name,
+        exgmnameb=gwt2.name,
+        exchangedata=gwfgwf_data,
+        auxiliary=["ANGLDEGX", "CDIST"],
+        filename="transport1_transport2.gwtgwt",
+    )
+
+    # GWF GWT exchange
+    gwfgwt1 = flopy.mf6.ModflowGwfgwt(
+        sim,
+        exgtype="GWF6-GWT6",
+        exgmnamea="flow1",
+        exgmnameb="transport1",
+        filename="flow1_transport1.gwfgwt",
+    )
+    gwfgwt2 = flopy.mf6.ModflowGwfgwt(
+        sim,
+        exgtype="GWF6-GWT6",
+        exgmnamea="flow2",
+        exgmnameb="transport2",
+        filename="flow2_transport2.gwfgwt",
+    )
+
+    # create iterative model solution and register the gwt model with it
+    imsgwt = 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="transport.ims",
+    )
+    sim.register_ims_package(imsgwt, [gwt1.name, gwt2.name])
+
+    return sim, None
+
+
+def check_output(idx, test):
+    # for now, just check that simulation finishes correctly
+    return
+
+
+@pytest.mark.parametrize("idx, name", enumerate(cases))
+@pytest.mark.developmode
+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),
+    )
+    test.run()

autotest/test_par_gwt_adv03.py

@@ -0,0 +1,28 @@
+"""
+This test reuses the simulation data and config in
+test_gwt_adv03_gwtgwt.py and runs it in parallel mode.
+"""
+
+import pytest
+from framework import TestFramework
+
+cases = ["par_adv03_gwtgwt"]
+
+
+@pytest.mark.parallel
+@pytest.mark.developmode
+@pytest.mark.parametrize("idx, name", enumerate(cases))
+def test_mf6model(idx, name, function_tmpdir, targets):
+    from test_gwt_adv03_gwtgwt import build_models, check_output
+
+    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()

doc/ReleaseNotes/develop.toml

@@ -88,3 +88,8 @@ section = "fixes"
 subsection = "solution"
 description = "Fixes a bug first included in version 6.4.2 that has changed the convergence check of the numerical solution. In cases where the linear solver failed to converge, or when it did converge but the STRICT option was active and convergence did not occur on the first iteration, the solution was erroneously considered for overall convergence."
 
+[[items]]
+section = "fixes"
+subsection = "parallel"
+description = "Fixes a memory exception that has occurred when running parallel solute transport on Windows without the dispersion package. This has no effect on the simulation results other than causing a premature termination in some cases."
+