Issue #1928 add utvd support

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
 ~~~~~

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

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)
 
 # %%

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

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
@@ -135,8 +136,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 +161,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.
+    <https://modflow6-examples.readthedocs.io/en/latest/_examples/ex-gwt-adv-schemes.html>`_
+    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"

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")]