verify partitioning results for anisotropy in k

[Manangka]

In GitLab by @luitjansl on Nov 16, 2023, 13:51

The following is a starting point for a test with anisotropy in k. After partitioning the simulation and running the simulation, we should get the same results as when running the simulation without partitioning. But this is not the case. We should investigate if the same happens in flopy and if not, why.

@pytest.mark.usefixtures("transient_twri_model")
@pytest.mark.parametrize(
    "partition_name",
    ["diagonal_1", "diagonal_2", "four_squares", "intrusion", "island"],
)
def test_partitioning_strong_k_anisotropy(
    tmp_path: Path, transient_twri_model: Modflow6Simulation, partition_name: str
):
    simulation = transient_twri_model
    idomain = simulation["GWF_1"].domain
    k = transient_twri_model["GWF_1"]["npf"].dataset["k"]
    transient_twri_model["GWF_1"]["npf"].dataset["k22"]= k *100

    partitioning_arrays = setup_partitioning_arrays(idomain.isel(layer=0))
    # run the original example, so without partitioning, and save the simulation results
    orig_dir = tmp_path / "original"
    simulation.write(orig_dir, binary=False, use_absolute_paths=True)
    #-----------------------------
    flopy_sim =flopy.mf6.MFSimulation.load( sim_ws=orig_dir, verbosity_level=1,)    

    mf_splitter = Mf6Splitter(flopy_sim)
    flopy_split_sim =  mf_splitter.split_model(partitioning_arrays[partition_name])
    flopy_split_dir = tmp_path / "flopy_split"   
    flopy_split_sim.set_sim_path(flopy_split_dir)
    flopy_split_sim.write_simulation(silent=False)
    flopy_split_sim.run_simulation(silent=False)    
    #-----------------------------
    simulation.run()

    orig_head = imod.mf6.open_hds(
        orig_dir / "GWF_1/GWF_1.hds",
        orig_dir / "GWF_1/dis.dis.grb",
    )

    # partition the simulation, run it, and save the (merged) results

    split_simulation = simulation.split(partitioning_arrays[partition_name])

    split_simulation.write(tmp_path, binary=False)
    split_simulation.run()

    head = merge_heads(tmp_path, split_simulation)
    _ = merge_balances(tmp_path, split_simulation)

    # compare the head result of the original simulation with the result of the partitioned simulation
    np.testing.assert_allclose(head.values, orig_head.values, rtol=1e-4, atol=1e-4)

[Manangka]

In GitLab by @luitjansl on Nov 16, 2023, 15:12

tests were done with the twri model. It had anisotropy added to it (K22 was set to 20*K in the npf package). It was partitioned in 4 different ways, and it was partitioned using flopy and imod-python. FIrst it was established that without partitioning, flopy and imod-python give the same results.

Then the difference between the unpartitioned head and the partitioned head was computed for the last timestep in the simulation.

The following table gives the max and mean difference for the head results. Note that these results are exactly equal for flopy and imod-python, except for the "four squares"partitioning. There, the flopy error is larger than the imod-python error. So the story was considered validated.	Partition boundary	Flopy split max	Flopy split mean	Imod python split max	Imod split mean
Diagonal_1	0.5191836798861	-2.220987858	0.5191836798861	-2.220987858
Island	4.494228506469	-0.02097133	4.494228506469	-0.02097133
intrusion	0.86444	0.0078825572697	0.86444	0.0078825572697
Four squares	12.096101364	-20.828593	12.096101364	3.88502683

Note: when imod-python creates a partition file, it describes the exchange from the partition with the highest partition number to the partition with the lowest partition number. So from partition 2 to partition 1. Flopy does this the other way around and gives the partition from 1 to 2. we changed imod python to do the same thing as flopy, to facilitate comparisons.