From c6a8700f0997376497f88a59989a353b9b46a547 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Thu, 8 Oct 2026 11:49:48 +0200 Subject: [PATCH 1/7] Implement UTVD advection scheme --- imod/mf6/adv.py | 35 +++++++++++++++++++++++++++-- imod/tests/test_mf6/test_mf6_adv.py | 9 ++++++++ 2 files changed, 42 insertions(+), 2 deletions(-) diff --git a/imod/mf6/adv.py b/imod/mf6/adv.py index e2dd0438e..b39c6da02 100644 --- a/imod/mf6/adv.py +++ b/imod/mf6/adv.py @@ -135,8 +135,9 @@ class AdvectionCentral(AdvectionBase): class AdvectionTVD(AdvectionBase): """ - An implicit second order TVD scheme. More expensive than upstream - weighting but more robust. + An implicit second order TVD scheme. More expensive than upstream weighting + but more robust. This scheme is not advised for unstructured grids, for that + see: :class:`AdvectionUTVD` Parameters ---------- @@ -159,3 +160,33 @@ class AdvectionTVD(AdvectionBase): """ _scheme = "TVD" + + +class AdvectionUTVD(AdvectionBase): + """ + An implicit second order TVD scheme for unstructured grids. `This scheme + seems to perform best in preserving sharp fronts for both structured and + unstructured grids. + `_ + Requires MODFLOW 6.7.0 or higher. + + Parameters + ---------- + ats_percel: float, optional + Fractional cell distance submitted by the ADV Package to the + :class:`imod.mf6.AdaptiveTimeStepping` (ATS) package. If ``ats_percel`` + is specified and the ATS Package is active, a time step calculation will + be made for each cell based on flow through the cell and cell + properties. The largest time step will be calculated such that the + advective fractional cell distance (``ats_percel``) is not exceeded for + any active cell in the grid. This time-step constraint will be submitted + to the ATS Package, perhaps with constraints submitted by other + packages, in the calculation of the time step. ``ats_percel`` must be + greater than zero. If a value of zero is specified for ``ats_percel`` + the program will automatically reset it to an internal no data value to + indicate that time steps should not be subject to this constraint. + validate: bool, optional + Validate the package upon initialization. Defaults to True. + """ + + _scheme = "UTVD" diff --git a/imod/tests/test_mf6/test_mf6_adv.py b/imod/tests/test_mf6/test_mf6_adv.py index 053d0ff64..effbf2549 100644 --- a/imod/tests/test_mf6/test_mf6_adv.py +++ b/imod/tests/test_mf6/test_mf6_adv.py @@ -34,6 +34,15 @@ def test_advection_TVD(): assert actual == expected +def test_advection_UTVD(): + directory = pathlib.Path("mymodel") + globaltimes = [np.datetime64("2000-01-01")] + adv = imod.mf6.AdvectionUTVD() + actual = adv._render(directory, "adv", globaltimes, True) + expected = "begin options\n scheme UTVD\nend options" + assert actual == expected + + def test_advection_ats_percel(): directory = pathlib.Path("mymodel") globaltimes = [np.datetime64("2000-01-01")] From a402e260a32c9d2ce7a6ed9d35e067cfdff834f6 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Thu, 8 Oct 2026 11:52:04 +0200 Subject: [PATCH 2/7] Add UTVD to public API --- imod/mf6/__init__.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/imod/mf6/__init__.py b/imod/mf6/__init__.py index 51306e43c..32ec8eb40 100644 --- a/imod/mf6/__init__.py +++ b/imod/mf6/__init__.py @@ -2,7 +2,12 @@ Create a Modflow 6 model. """ -from imod.mf6.adv import AdvectionCentral, AdvectionTVD, AdvectionUpstream +from imod.mf6.adv import ( + AdvectionCentral, + AdvectionTVD, + AdvectionUpstream, + AdvectionUTVD, +) from imod.mf6.api_package import ApiPackage from imod.mf6.ats import AdaptiveTimeStepping from imod.mf6.buy import Buoyancy From e1ded2e0c06e7b771ff934885653c7d7e7808f64 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Thu, 8 Oct 2026 11:52:50 +0200 Subject: [PATCH 3/7] Add UTVD to API docs --- docs/api/mf6.rst | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/docs/api/mf6.rst b/docs/api/mf6.rst index dda794721..d800a1944 100644 --- a/docs/api/mf6.rst +++ b/docs/api/mf6.rst @@ -449,6 +449,13 @@ Transport Packages AdvectionTVD.copy AdvectionTVD.is_empty AdvectionTVD.get_regrid_methods + AdvectionUTVD + AdvectionUTVD.write + AdvectionUTVD.from_file + AdvectionUTVD.to_file + AdvectionUTVD.copy + AdvectionUTVD.is_empty + AdvectionUTVD.get_regrid_methods AdvectionUpstream AdvectionUpstream.write AdvectionUpstream.from_file From 4d5b00886cba95305df13ccf83e9a2c80c6ecd48 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Thu, 8 Oct 2026 11:55:43 +0200 Subject: [PATCH 4/7] Update changelog --- docs/api/changelog.rst | 3 +++ 1 file changed, 3 insertions(+) diff --git a/docs/api/changelog.rst b/docs/api/changelog.rst index 88cd2a0a2..1b4992f32 100644 --- a/docs/api/changelog.rst +++ b/docs/api/changelog.rst @@ -39,6 +39,9 @@ Added the final package has fully empty layers trimmed off afterwards, so this does not affect computed values. Set to False to keep the previous full-layer-coordinate behaviour. +- Added :class:`imod.mf6.AdvectionUTVD` to the transport packages. This is the + most accurate advection scheme in most cases, and is recommended for when + sharp fronts need to be preserved or when computing on an unstructured grid. Fixed ~~~~~ From e9743278d311b2138b3043bf1e71e324d7ab48cc Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Thu, 8 Oct 2026 11:57:45 +0200 Subject: [PATCH 5/7] Use UTVD in unstructured transport example --- examples/mf6/circle_transport.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/mf6/circle_transport.py b/examples/mf6/circle_transport.py index 7c0f8aeda..2a6bdd618 100644 --- a/examples/mf6/circle_transport.py +++ b/examples/mf6/circle_transport.py @@ -264,7 +264,7 @@ xt3d_off=False, xt3d_rhs=False, ) -transport_model["adv"] = imod.mf6.AdvectionTVD(ats_percel=0.95) +transport_model["adv"] = imod.mf6.AdvectionUTVD(ats_percel=0.95) transport_model["mst"] = imod.mf6.MobileStorageTransfer(porosity) # %% From 8184f3e36d68b7760ce07524a5075480861f39f3 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Thu, 8 Oct 2026 14:52:06 +0200 Subject: [PATCH 6/7] Fix docstring after review comment --- imod/mf6/adv.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/imod/mf6/adv.py b/imod/mf6/adv.py index b39c6da02..4f5a0f6b6 100644 --- a/imod/mf6/adv.py +++ b/imod/mf6/adv.py @@ -1,12 +1,13 @@ """ When simulating transport, MODFLOW6 needs to compute the concentration at a -cellface between 2 adjacent cells. It supports 3 ways of doing this. Each of +cellface between 2 adjacent cells. It supports 4 ways of doing this. Each of those has its own wrapper class. These numerical schemes differ in how much numerical dispersion they cause, how much oscillations, and how timestep and -grid size affect stability. Central-in-space weighting is not often used -because it can result in spurious oscillations in the simulated concentrations. -Upstream weighting is a fast alternative, and TVD is a more expensive and more -robust alternative. +grid size affect stability. Central-in-space weighting is not often used because +it can result in spurious oscillations in the simulated concentrations. Upstream +weighting is a fast alternative, TVD is a more expensive and more robust +alternative, and UTVD is the most accurate for sharp concentration fronts and +unstructured grids. """ from abc import ABC From 3cff2a54e1a0573176e5e1571b2abe8c71b3b7ab Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Mon, 12 Oct 2026 09:48:26 +0200 Subject: [PATCH 7/7] Bump MODFLOW version for CI to 6.8.1 --- .teamcity/Deploy/DeployProject.kt | 2 +- .teamcity/_Self/MainProject.kt | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/.teamcity/Deploy/DeployProject.kt b/.teamcity/Deploy/DeployProject.kt index a25f8c4f3..f30eda82f 100644 --- a/.teamcity/Deploy/DeployProject.kt +++ b/.teamcity/Deploy/DeployProject.kt @@ -106,7 +106,7 @@ object BuildPages : BuildType({ } } params { - param("reverse.dep.Modflow_Modflow6Release.MODFLOW6_Version", "6.6.3") + param("reverse.dep.Modflow_Modflow6Release.MODFLOW6_Version", "6.8.1") param("reverse.dep.Modflow_Modflow6Release.MODFLOW6_Platform", "win64") } dependencies { diff --git a/.teamcity/_Self/MainProject.kt b/.teamcity/_Self/MainProject.kt index 600e741c4..fa5c8e3d2 100644 --- a/.teamcity/_Self/MainProject.kt +++ b/.teamcity/_Self/MainProject.kt @@ -99,7 +99,7 @@ object UnitTests : BuildType({ name = "UnitTests" params { - param("reverse.dep.Modflow_Modflow6Release.MODFLOW6_Version", "6.6.3") + param("reverse.dep.Modflow_Modflow6Release.MODFLOW6_Version", "6.8.1") param("reverse.dep.Modflow_Modflow6Release.MODFLOW6_Platform", "win64") } templates(UnitTestsTemplate, GitHubIntegrationTemplate) @@ -128,7 +128,7 @@ object Examples : BuildType({ name = "Examples" params { - param("reverse.dep.Modflow_Modflow6Release.MODFLOW6_Version", "6.6.3") + param("reverse.dep.Modflow_Modflow6Release.MODFLOW6_Version", "6.8.1") param("reverse.dep.Modflow_Modflow6Release.MODFLOW6_Platform", "win64") }