diff --git a/.codespell.ignore b/.codespell.ignore
index 30296cf990b..55ad98775f7 100644
--- a/.codespell.ignore
+++ b/.codespell.ignore
@@ -2,6 +2,7 @@ Alph
wel
nam
gage
+DisUtils
drob
delt
lke
diff --git a/autotest/test_gwt_adv01.py b/autotest/test_gwt_adv01.py
index 97069a7f29b..9150f6b9b5d 100644
--- a/autotest/test_gwt_adv01.py
+++ b/autotest/test_gwt_adv01.py
@@ -10,8 +10,8 @@
import pytest
from framework import TestFramework
-cases = ["adv01a", "adv01b", "adv01c"]
-scheme = ["upstream", "central", "tvd"]
+cases = ["adv01a", "adv01b", "adv01c", "adv01d"]
+scheme = ["upstream", "central", "tvd", "utvd"]
def build_models(idx, test):
@@ -569,7 +569,115 @@ def check_output(idx, test):
]
cres3 = np.array(cres3)
- creslist = [cres1, cres2, cres3]
+ cres4 = [
+ [
+ [
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 9.99999999e-01,
+ 9.99999997e-01,
+ 9.99999991e-01,
+ 9.99999975e-01,
+ 9.99999926e-01,
+ 9.99999789e-01,
+ 9.99999407e-01,
+ 9.99998374e-01,
+ 9.99995665e-01,
+ 9.99988785e-01,
+ 9.99971918e-01,
+ 9.99932078e-01,
+ 9.99841550e-01,
+ 9.99643930e-01,
+ 9.99229970e-01,
+ 9.98398720e-01,
+ 9.96800070e-01,
+ 9.93857995e-01,
+ 9.88681096e-01,
+ 9.79978744e-01,
+ 9.66015902e-01,
+ 9.44652308e-01,
+ 9.13514114e-01,
+ 8.70328697e-01,
+ 8.13410723e-01,
+ 7.42224213e-01,
+ 6.57879958e-01,
+ 5.63390872e-01,
+ 4.63530314e-01,
+ 3.64233330e-01,
+ 2.71628520e-01,
+ 1.90935414e-01,
+ 1.25541011e-01,
+ 7.65316278e-02,
+ 4.28052306e-02,
+ 2.16851952e-02,
+ 9.78980420e-03,
+ 3.85619551e-03,
+ 1.28879215e-03,
+ 3.52115803e-04,
+ 7.49350468e-05,
+ 1.17646292e-05,
+ 1.32776825e-06,
+ 9.64125296e-08,
+ -5.86919920e-08,
+ -7.40823239e-08,
+ -6.31800862e-08,
+ -5.27398214e-08,
+ -4.32274354e-08,
+ -3.48423105e-08,
+ -2.76485491e-08,
+ -2.16203169e-08,
+ -1.66732861e-08,
+ -1.26895729e-08,
+ -9.53670757e-09,
+ -7.08114746e-09,
+ -5.19714853e-09,
+ -3.77193598e-09,
+ -2.70809957e-09,
+ -1.92404050e-09,
+ -1.35315643e-09,
+ -9.42300766e-10,
+ -6.49908787e-10,
+ -4.44059825e-10,
+ -3.00645032e-10,
+ -2.01734707e-10,
+ -1.34185608e-10,
+ -8.84931133e-11,
+ -5.78716880e-11,
+ -3.75358997e-11,
+ -2.41500917e-11,
+ -1.54150906e-11,
+ -9.76314930e-12,
+ -6.13633357e-12,
+ -3.82789021e-12,
+ -2.37025821e-12,
+ -1.07644858e-12,
+ -2.60415045e-12,
+ ]
+ ]
+ ]
+ cres4 = np.array(cres4)
+
+ creslist = [cres1, cres2, cres3, cres4]
assert np.allclose(creslist[idx], conc), (
"simulated concentrations do not match with known solution."
diff --git a/autotest/test_gwt_adv01_fmi.py b/autotest/test_gwt_adv01_fmi.py
index 97a3919df19..2d7e6fe593a 100644
--- a/autotest/test_gwt_adv01_fmi.py
+++ b/autotest/test_gwt_adv01_fmi.py
@@ -12,8 +12,8 @@
from flopy.utils.gridutil import uniform_flow_field
from framework import TestFramework
-cases = ["adv01a_fmi", "adv01b_fmi", "adv01c_fmi"]
-scheme = ["upstream", "central", "tvd"]
+cases = ["adv01a_fmi", "adv01b_fmi", "adv01c_fmi", "adv01d_fmi"]
+scheme = ["upstream", "central", "tvd", "utvd"]
def build_models(idx, test):
@@ -539,7 +539,115 @@ def check_output(idx, test):
]
cres3 = np.array(cres3)
- creslist = [cres1, cres2, cres3]
+ cres4 = [
+ [
+ [
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 9.99999999e-01,
+ 9.99999997e-01,
+ 9.99999991e-01,
+ 9.99999975e-01,
+ 9.99999926e-01,
+ 9.99999789e-01,
+ 9.99999407e-01,
+ 9.99998374e-01,
+ 9.99995665e-01,
+ 9.99988785e-01,
+ 9.99971918e-01,
+ 9.99932078e-01,
+ 9.99841550e-01,
+ 9.99643930e-01,
+ 9.99229970e-01,
+ 9.98398720e-01,
+ 9.96800070e-01,
+ 9.93857995e-01,
+ 9.88681096e-01,
+ 9.79978744e-01,
+ 9.66015902e-01,
+ 9.44652308e-01,
+ 9.13514114e-01,
+ 8.70328697e-01,
+ 8.13410723e-01,
+ 7.42224213e-01,
+ 6.57879958e-01,
+ 5.63390872e-01,
+ 4.63530314e-01,
+ 3.64233330e-01,
+ 2.71628520e-01,
+ 1.90935414e-01,
+ 1.25541011e-01,
+ 7.65316278e-02,
+ 4.28052306e-02,
+ 2.16851952e-02,
+ 9.78980420e-03,
+ 3.85619551e-03,
+ 1.28879215e-03,
+ 3.52115803e-04,
+ 7.49350468e-05,
+ 1.17646292e-05,
+ 1.32776825e-06,
+ 9.64125296e-08,
+ -5.86919920e-08,
+ -7.40823239e-08,
+ -6.31800862e-08,
+ -5.27398214e-08,
+ -4.32274354e-08,
+ -3.48423105e-08,
+ -2.76485491e-08,
+ -2.16203169e-08,
+ -1.66732861e-08,
+ -1.26895729e-08,
+ -9.53670757e-09,
+ -7.08114746e-09,
+ -5.19714853e-09,
+ -3.77193598e-09,
+ -2.70809957e-09,
+ -1.92404050e-09,
+ -1.35315643e-09,
+ -9.42300766e-10,
+ -6.49908787e-10,
+ -4.44059825e-10,
+ -3.00645032e-10,
+ -2.01734707e-10,
+ -1.34185608e-10,
+ -8.84931133e-11,
+ -5.78716880e-11,
+ -3.75358997e-11,
+ -2.41500917e-11,
+ -1.54150906e-11,
+ -9.76314930e-12,
+ -6.13633357e-12,
+ -3.82789021e-12,
+ -2.37025821e-12,
+ -1.07644858e-12,
+ -2.60415045e-12,
+ ]
+ ]
+ ]
+ cres4 = np.array(cres4)
+
+ creslist = [cres1, cres2, cres3, cres4]
assert np.allclose(creslist[idx], conc), (
"simulated concentrations do not match with known solution."
diff --git a/autotest/test_gwt_adv01_gwtgwt.py b/autotest/test_gwt_adv01_gwtgwt.py
index 08cb6366aa6..d584e6bf6ab 100644
--- a/autotest/test_gwt_adv01_gwtgwt.py
+++ b/autotest/test_gwt_adv01_gwtgwt.py
@@ -10,8 +10,8 @@
import pytest
from framework import TestFramework
-cases = ["adv01a_gwtgwt", "adv01b_gwtgwt", "adv01c_gwtgwt"]
-scheme = ["upstream", "central", "tvd"]
+cases = ["adv01a_gwtgwt", "adv01b_gwtgwt", "adv01c_gwtgwt", "adv01d_gwtgwt"]
+scheme = ["upstream", "central", "tvd", "utvd"]
gdelr = 1.0
# solver settings
@@ -668,7 +668,115 @@ def check_output(idx, test):
]
cres3 = np.array(cres3)
- creslist = [cres1, cres2, cres3]
+ cres4 = [
+ [
+ [
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 9.99999999e-01,
+ 9.99999997e-01,
+ 9.99999991e-01,
+ 9.99999975e-01,
+ 9.99999926e-01,
+ 9.99999789e-01,
+ 9.99999407e-01,
+ 9.99998374e-01,
+ 9.99995665e-01,
+ 9.99988785e-01,
+ 9.99971918e-01,
+ 9.99932078e-01,
+ 9.99841550e-01,
+ 9.99643930e-01,
+ 9.99229970e-01,
+ 9.98398720e-01,
+ 9.96800070e-01,
+ 9.93857995e-01,
+ 9.88681096e-01,
+ 9.79978744e-01,
+ 9.66015902e-01,
+ 9.44652308e-01,
+ 9.13514114e-01,
+ 8.70328697e-01,
+ 8.13410723e-01,
+ 7.42224213e-01,
+ 6.57879958e-01,
+ 5.63390872e-01,
+ 4.63530314e-01,
+ 3.64233330e-01,
+ 2.71628520e-01,
+ 1.90935414e-01,
+ 1.25541011e-01,
+ 7.65316278e-02,
+ 4.28052306e-02,
+ 2.16851952e-02,
+ 9.78980420e-03,
+ 3.85619551e-03,
+ 1.28879215e-03,
+ 3.52115803e-04,
+ 7.49350468e-05,
+ 1.17646292e-05,
+ 1.32776825e-06,
+ 9.64125296e-08,
+ -5.86919920e-08,
+ -7.40823239e-08,
+ -6.31800862e-08,
+ -5.27398214e-08,
+ -4.32274354e-08,
+ -3.48423105e-08,
+ -2.76485491e-08,
+ -2.16203169e-08,
+ -1.66732861e-08,
+ -1.26895729e-08,
+ -9.53670757e-09,
+ -7.08114746e-09,
+ -5.19714853e-09,
+ -3.77193598e-09,
+ -2.70809957e-09,
+ -1.92404050e-09,
+ -1.35315643e-09,
+ -9.42300766e-10,
+ -6.49908787e-10,
+ -4.44059825e-10,
+ -3.00645032e-10,
+ -2.01734707e-10,
+ -1.34185608e-10,
+ -8.84931133e-11,
+ -5.78716880e-11,
+ -3.75358997e-11,
+ -2.41500917e-11,
+ -1.54150906e-11,
+ -9.76314930e-12,
+ -6.13633357e-12,
+ -3.82789021e-12,
+ -2.37025821e-12,
+ -1.07644858e-12,
+ -2.60415045e-12,
+ ]
+ ]
+ ]
+ cres4 = np.array(cres4)
+
+ creslist = [cres1, cres2, cres3, cres4]
assert np.allclose(creslist[idx], conc), (
"simulated concentrations do not match with known solution."
diff --git a/autotest/test_gwt_adv02.py b/autotest/test_gwt_adv02.py
index 07dbe19463d..bfd82b641fa 100644
--- a/autotest/test_gwt_adv02.py
+++ b/autotest/test_gwt_adv02.py
@@ -13,8 +13,8 @@
import pytest
from framework import TestFramework
-cases = ["adv02a", "adv02b", "adv02c"]
-scheme = ["upstream", "central", "tvd"]
+cases = ["adv02a", "adv02b", "adv02c", "adv02d"]
+scheme = ["upstream", "central", "tvd", "utvd"]
def grid_triangulator(itri, delr, delc):
@@ -922,7 +922,213 @@ def check_output(idx, test):
]
cres3 = np.array(cres3)
- creslist = [cres1, cres2, cres3]
+ cres4 = [
+ [
+ [
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 1.00000000e00,
+ 9.99999999e-01,
+ 9.99999999e-01,
+ 9.99999998e-01,
+ 9.99999997e-01,
+ 9.99999996e-01,
+ 9.99999993e-01,
+ 9.99999988e-01,
+ 9.99999980e-01,
+ 9.99999967e-01,
+ 9.99999945e-01,
+ 9.99999909e-01,
+ 9.99999849e-01,
+ 9.99999748e-01,
+ 9.99999580e-01,
+ 9.99999300e-01,
+ 9.99998835e-01,
+ 9.99998064e-01,
+ 9.99996788e-01,
+ 9.99994679e-01,
+ 9.99991207e-01,
+ 9.99985519e-01,
+ 9.99976245e-01,
+ 9.99961221e-01,
+ 9.99937056e-01,
+ 9.99898507e-01,
+ 9.99837576e-01,
+ 9.99742235e-01,
+ 9.99594683e-01,
+ 9.99369014e-01,
+ 9.99028179e-01,
+ 9.98520185e-01,
+ 9.97773486e-01,
+ 9.96691640e-01,
+ 9.95147449e-01,
+ 9.92976938e-01,
+ 9.89973781e-01,
+ 9.85884941e-01,
+ 9.80408526e-01,
+ 9.73194950e-01,
+ 9.63852496e-01,
+ 9.51958228e-01,
+ 9.37074848e-01,
+ 9.18773513e-01,
+ 8.96661955e-01,
+ 8.70416363e-01,
+ 8.39814727e-01,
+ 8.04768610e-01,
+ 7.65349931e-01,
+ 7.21809331e-01,
+ 6.74583139e-01,
+ 6.24286813e-01,
+ 5.71694587e-01,
+ 5.17708928e-01,
+ 4.63475646e-01,
+ 4.10621606e-01,
+ 3.60090882e-01,
+ 3.12561498e-01,
+ 2.68548170e-01,
+ 2.28398779e-01,
+ 1.92300754e-01,
+ 1.60295787e-01,
+ 1.32300488e-01,
+ 1.08130463e-01,
+ 8.75255237e-02,
+ 7.01741203e-02,
+ 5.57355716e-02,
+ 4.38591512e-02,
+ 3.41995388e-02,
+ 2.64285307e-02,
+ 2.02431998e-02,
+ 1.53708955e-02,
+ 1.15715937e-02,
+ 8.63815025e-03,
+ 6.39500342e-03,
+ 4.69581865e-03,
+ 3.42049930e-03,
+ 2.47190356e-03,
+ 1.77252555e-03,
+ 1.26132360e-03,
+ 8.90813017e-04,
+ 6.24488550e-04,
+ 4.34601338e-04,
+ 3.00286592e-04,
+ 2.06019188e-04,
+ 1.40363342e-04,
+ 9.49774268e-05,
+ 6.38341487e-05,
+ 4.26182859e-05,
+ 2.82678226e-05,
+ 1.86287813e-05,
+ 1.21986984e-05,
+ 7.93815930e-06,
+ 5.13384506e-06,
+ 3.30005212e-06,
+ 2.10858694e-06,
+ 1.33934239e-06,
+ 8.45779171e-07,
+ 5.31034450e-07,
+ 3.31530143e-07,
+ 2.05821848e-07,
+ 1.27075205e-07,
+ 7.80302654e-08,
+ 4.76573231e-08,
+ 2.89528335e-08,
+ 1.74975602e-08,
+ 1.05200557e-08,
+ 6.29276848e-09,
+ 3.74521020e-09,
+ 2.21793995e-09,
+ 1.30703695e-09,
+ 7.66512576e-10,
+ 4.47368886e-10,
+ 2.59871508e-10,
+ 1.50249706e-10,
+ 8.64692928e-11,
+ 4.95369461e-11,
+ 2.82485119e-11,
+ 1.60400425e-11,
+ 9.06356925e-12,
+ 5.10160108e-12,
+ 2.85609391e-12,
+ 1.59463729e-12,
+ 8.83001720e-13,
+ 4.88912295e-13,
+ 2.66402206e-13,
+ 1.47206301e-13,
+ 7.85891341e-14,
+ 4.65296803e-14,
+ 3.41047008e-14,
+ 2.48440107e-14,
+ 1.80183366e-14,
+ 1.30166484e-14,
+ 9.36861656e-15,
+ 6.71933416e-15,
+ 4.80311981e-15,
+ 3.42235516e-15,
+ 2.43094613e-15,
+ 1.72149946e-15,
+ 1.21547008e-15,
+ 8.55670685e-16,
+ 6.00634004e-16,
+ 4.20405055e-16,
+ 2.93422427e-16,
+ 2.04219766e-16,
+ 1.41739986e-16,
+ 9.81042000e-17,
+ 6.77163197e-17,
+ 4.66142906e-17,
+ 3.20017576e-17,
+ 2.19112061e-17,
+ 1.49625226e-17,
+ 7.71070649e-18,
+ ]
+ ]
+ ]
+ cres4 = np.array(cres4)
+
+ creslist = [cres1, cres2, cres3, cres4]
assert np.allclose(creslist[idx], conc), (
"simulated concentrations do not match with known solution."
diff --git a/autotest/test_gwt_adv02_gwtgwt.py b/autotest/test_gwt_adv02_gwtgwt.py
index 1bd064cf78e..fbcc6f73887 100644
--- a/autotest/test_gwt_adv02_gwtgwt.py
+++ b/autotest/test_gwt_adv02_gwtgwt.py
@@ -12,8 +12,8 @@
import pytest
from framework import TestFramework
-cases = ["adv02a_gwtgwt", "adv02b_gwtgwt", "adv02c_gwtgwt"]
-scheme = ["upstream", "central", "tvd"]
+cases = ["adv02a_gwtgwt", "adv02b_gwtgwt", "adv02c_gwtgwt", "adv02d_gwtgwt"]
+scheme = ["upstream", "central", "tvd", "utvd"]
gdelr = 1.0
# solver settings
diff --git a/autotest/test_gwt_adv03.py b/autotest/test_gwt_adv03.py
index bf660c4d515..312ba5d9371 100644
--- a/autotest/test_gwt_adv03.py
+++ b/autotest/test_gwt_adv03.py
@@ -13,8 +13,8 @@
import pytest
from framework import TestFramework
-cases = ["adv03a", "adv03b", "adv03c"]
-scheme = ["upstream", "central", "tvd"]
+cases = ["adv03a", "adv03b", "adv03c", "adv03d"]
+scheme = ["upstream", "central", "tvd", "utvd"]
def grid_triangulator(itri, delr, delc):
@@ -128,7 +128,8 @@ def build_models(idx, test):
sim.register_ims_package(imsgwf, [gwf.name])
itri = np.zeros((nrow, ncol), dtype=int)
- itri[:, 1 : ncol - 1] = 1
+ itri[int(nrow / 2) :, 1 : ncol - 1] = 1
+ itri[: int(nrow / 2), 1 : ncol - 1] = 2
verts, iverts = grid_triangulator(itri, delr, delc)
vertices, cell2d = cvfd_to_cell2d(verts, iverts)
ncpl = len(cell2d)
@@ -463,11 +464,53 @@ def check_output(idx, test):
]
cres3 = np.array(cres3)
+ cres4 = [
+ 9.75305992e-01,
+ 9.60923712e-01,
+ 9.40345041e-01,
+ 8.98518356e-01,
+ 8.43452227e-01,
+ 7.50966209e-01,
+ 6.55085404e-01,
+ 5.23794573e-01,
+ 4.15360236e-01,
+ 2.94697679e-01,
+ 2.10720136e-01,
+ 1.35962806e-01,
+ 8.93878174e-02,
+ 5.38698838e-02,
+ 3.32165717e-02,
+ 1.88821558e-02,
+ 1.10910530e-02,
+ 5.99079841e-03,
+ 3.39811390e-03,
+ 1.75392691e-03,
+ 9.71354219e-04,
+ 4.81791201e-04,
+ 2.62200432e-04,
+ 1.25921717e-04,
+ 6.72680080e-05,
+ 3.15465073e-05,
+ 1.64347475e-05,
+ 7.56625500e-06,
+ 3.83544217e-06,
+ 1.73771238e-06,
+ 8.59572442e-07,
+ 3.83960616e-07,
+ 1.85884337e-07,
+ 8.19811804e-08,
+ 3.89673178e-08,
+ 1.69141522e-08,
+ 8.17467605e-09,
+ 1.47129715e-09,
+ ]
+ cres4 = np.array(cres4)
+
# Compare the first row in the layer with the answer and compare the
# last row in the bottom layer with the answer. This will verify that
# the results are one-dimensional even though the model is three
# dimensional
- creslist = [cres1, cres2, cres3]
+ creslist = [cres1, cres2, cres3, cres4]
ncellsperrow = cres1.shape[0]
assert np.allclose(creslist[idx], conc[0, 0, 0:ncellsperrow]), (
"simulated concentrations do not match with known solution.",
diff --git a/autotest/test_gwt_adv04.py b/autotest/test_gwt_adv04.py
index 6c0b1b5702a..b79fab22b92 100644
--- a/autotest/test_gwt_adv04.py
+++ b/autotest/test_gwt_adv04.py
@@ -11,8 +11,8 @@
import pytest
from framework import TestFramework
-cases = ["adv04a", "adv04b", "adv04c"]
-scheme = ["upstream", "central", "tvd"]
+cases = ["adv04a", "adv04b", "adv04c", "adv04d"]
+scheme = ["upstream", "central", "tvd", "utvd"]
def build_models(idx, test):
diff --git a/autotest/test_gwt_disu01.py b/autotest/test_gwt_disu01.py
index 1addf3ec50a..bc2371b56ce 100644
--- a/autotest/test_gwt_disu01.py
+++ b/autotest/test_gwt_disu01.py
@@ -95,7 +95,9 @@ def get_nn(k, i, j):
sim.register_ims_package(imsgwf, [gwf.name])
# use utility to make a disu version of a regular grid
- disu_kwargs = get_disu_kwargs(nlay, nrow, ncol, delr, delc, top, botm)
+ disu_kwargs = get_disu_kwargs(
+ nlay, nrow, ncol, delr, delc, top, botm, return_vertices=True
+ )
disu = flopy.mf6.ModflowGwfdisu(gwf, **disu_kwargs)
# initial conditions
@@ -158,7 +160,9 @@ def get_nn(k, i, j):
sim.register_ims_package(imsgwt, [gwt.name])
# use utility to make a disu version of a regular grid
- disu_kwargs = get_disu_kwargs(nlay, nrow, ncol, delr, delc, top, botm)
+ disu_kwargs = get_disu_kwargs(
+ nlay, nrow, ncol, delr, delc, top, botm, return_vertices=True
+ )
disu = flopy.mf6.ModflowGwtdisu(gwt, **disu_kwargs)
# initial conditions
diff --git a/autotest/test_gwt_henry_gwtgwt.py b/autotest/test_gwt_henry_gwtgwt.py
index 7df20fbc6de..753fd59f9df 100644
--- a/autotest/test_gwt_henry_gwtgwt.py
+++ b/autotest/test_gwt_henry_gwtgwt.py
@@ -5,8 +5,13 @@
import pytest
from framework import TestFramework
-cases = ["henry01-gwtgwt-ups", "henry01-gwtgwt-cen", "henry01-gwtgwt-tvd"]
-advection_scheme = ["UPSTREAM", "CENTRAL", "TVD"]
+cases = [
+ "henry01-gwtgwt-ups",
+ "henry01-gwtgwt-cen",
+ "henry01-gwtgwt-tvd",
+ "henry01-gwtgwt-utvd",
+]
+advection_scheme = ["UPSTREAM", "CENTRAL", "TVD", "UTVD"]
lx = 2.0
lz = 1.0
diff --git a/autotest/test_par_gwt_henry.py b/autotest/test_par_gwt_henry.py
index 1c3f4c6cacd..445abe9123a 100644
--- a/autotest/test_par_gwt_henry.py
+++ b/autotest/test_par_gwt_henry.py
@@ -6,7 +6,7 @@
import pytest
from framework import TestFramework
-cases = ["par-henry-ups", "par-henry-cen", "par-henry-tvd"]
+cases = ["par-henry-ups", "par-henry-cen", "par-henry-tvd", "par-henry-utvd"]
def build_models(idx, test):
diff --git a/doc/mf6io/mf6ivar/dfn/gwt-adv.dfn b/doc/mf6io/mf6ivar/dfn/gwt-adv.dfn
index 4af532fe8de..519f0bc6d2b 100644
--- a/doc/mf6io/mf6ivar/dfn/gwt-adv.dfn
+++ b/doc/mf6io/mf6ivar/dfn/gwt-adv.dfn
@@ -3,11 +3,11 @@
block options
name scheme
type string
-valid central upstream tvd
+valid central upstream tvd utvd
reader urword
optional true
longname advective scheme
-description scheme used to solve the advection term. Can be upstream, central, or TVD. If not specified, upstream weighting is the default weighting scheme.
+description scheme used to solve the advection term. Can be upstream, central, TVD or UTVD. If not specified, upstream weighting is the default weighting scheme.
block options
name ats_percel
diff --git a/make/makefile b/make/makefile
index 7058838eb6e..0d89d9d465f 100644
--- a/make/makefile
+++ b/make/makefile
@@ -23,19 +23,19 @@ SOURCEDIR16=../src/Model/OverlandFlow
SOURCEDIR17=../src/Model/ParticleTracking
SOURCEDIR18=../src/Model/SurfaceWaterFlow
SOURCEDIR19=../src/Model/TransportModel
-SOURCEDIR20=../src/Model/TransportModel/InterpolationScheme
-SOURCEDIR21=../src/Solution
-SOURCEDIR22=../src/Solution/LinearMethods
-SOURCEDIR23=../src/Solution/PETSc
-SOURCEDIR24=../src/Solution/ParticleTracker
-SOURCEDIR25=../src/Solution/ParticleTracker/Domain
-SOURCEDIR26=../src/Solution/ParticleTracker/Method
-SOURCEDIR27=../src/Solution/ParticleTracker/Particle
-SOURCEDIR28=../src/Timing
-SOURCEDIR29=../src/Utilities
-SOURCEDIR30=../src/Utilities/ArrayRead
-SOURCEDIR31=../src/Utilities/Export
-SOURCEDIR32=../src/Utilities/Export/tmp
+SOURCEDIR20=../src/Model/TransportModel/Gradient
+SOURCEDIR21=../src/Model/TransportModel/InterpolationScheme
+SOURCEDIR22=../src/Solution
+SOURCEDIR23=../src/Solution/LinearMethods
+SOURCEDIR24=../src/Solution/PETSc
+SOURCEDIR25=../src/Solution/ParticleTracker
+SOURCEDIR26=../src/Solution/ParticleTracker/Domain
+SOURCEDIR27=../src/Solution/ParticleTracker/Method
+SOURCEDIR28=../src/Solution/ParticleTracker/Particle
+SOURCEDIR29=../src/Timing
+SOURCEDIR30=../src/Utilities
+SOURCEDIR31=../src/Utilities/ArrayRead
+SOURCEDIR32=../src/Utilities/Export
SOURCEDIR33=../src/Utilities/Idm
SOURCEDIR34=../src/Utilities/Idm/mf6blockfile
SOURCEDIR35=../src/Utilities/Idm/netcdf
@@ -314,6 +314,7 @@ $(OBJDIR)/SubcellRect.o \
$(OBJDIR)/MethodSubcellTernary.o \
$(OBJDIR)/MethodSubcellPollock.o \
$(OBJDIR)/TimeStepSelect.o \
+$(OBJDIR)/LinearAlgebraUtils.o \
$(OBJDIR)/SfrCrossSectionUtils.o \
$(OBJDIR)/MethodSubcellPool.o \
$(OBJDIR)/MethodCell.o \
@@ -322,6 +323,8 @@ $(OBJDIR)/CellRectQuad.o \
$(OBJDIR)/CellRect.o \
$(OBJDIR)/SwfCxsUtils.o \
$(OBJDIR)/PrintSaveManager.o \
+$(OBJDIR)/tsp-fmi.o \
+$(OBJDIR)/SingularValueDecomposition.o \
$(OBJDIR)/SfrCrossSectionManager.o \
$(OBJDIR)/dag_module.o \
$(OBJDIR)/BoundaryPackageExt.o \
@@ -333,8 +336,10 @@ $(OBJDIR)/VectorInterpolation.o \
$(OBJDIR)/swf-cxs.o \
$(OBJDIR)/OutputControlData.o \
$(OBJDIR)/gwf-ic.o \
-$(OBJDIR)/tsp-fmi.o \
$(OBJDIR)/InterpolationSchemeInterface.o \
+$(OBJDIR)/IGradient.o \
+$(OBJDIR)/DisUtils.o \
+$(OBJDIR)/PseudoInverse.o \
$(OBJDIR)/UzfEtUtil.o \
$(OBJDIR)/Xt3dAlgorithm.o \
$(OBJDIR)/TvBase.o \
@@ -355,9 +360,11 @@ $(OBJDIR)/BaseModel.o \
$(OBJDIR)/TspSpc.o \
$(OBJDIR)/GweInputData.o \
$(OBJDIR)/tsp-ic.o \
+$(OBJDIR)/UTVDScheme.o \
$(OBJDIR)/UpstreamScheme.o \
$(OBJDIR)/TVDScheme.o \
$(OBJDIR)/TspAdvOptions.o \
+$(OBJDIR)/LeastSquaresGradient.o \
$(OBJDIR)/CentralDifferenceScheme.o \
$(OBJDIR)/AdvSchemeEnum.o \
$(OBJDIR)/UzfCellGroup.o \
@@ -538,13 +545,10 @@ $(OBJDIR)/RunControlFactory.o \
$(OBJDIR)/IdmLoad.o \
$(OBJDIR)/ConnectionBuilder.o \
$(OBJDIR)/comarg.o \
-$(OBJDIR)/LinearAlgebraUtils.o \
$(OBJDIR)/mf6core.o \
-$(OBJDIR)/SingularValueDecomposition.o \
$(OBJDIR)/BaseGeometry.o \
$(OBJDIR)/mf6.o \
$(OBJDIR)/StringList.o \
-$(OBJDIR)/PseudoInverse.o \
$(OBJDIR)/MemorySetHandler.o \
$(OBJDIR)/ilut.o \
$(OBJDIR)/sparsekit.o \
diff --git a/msvs/mf6core.vfproj b/msvs/mf6core.vfproj
index 2617076250c..22117463b9c 100644
--- a/msvs/mf6core.vfproj
+++ b/msvs/mf6core.vfproj
@@ -261,6 +261,7 @@
+
@@ -378,12 +379,16 @@
+
+
+
-
+
+
diff --git a/src/Exchange/exg-gwegwe.f90 b/src/Exchange/exg-gwegwe.f90
index 9cff78d59b5..121929c29f7 100644
--- a/src/Exchange/exg-gwegwe.f90
+++ b/src/Exchange/exg-gwegwe.f90
@@ -29,6 +29,7 @@ module GweGweExchangeModule
use ObsModule, only: ObsType
use TableModule, only: TableType, table_cr
use MatrixBaseModule
+ use AdvSchemeEnumModule
implicit none
@@ -58,7 +59,7 @@ module GweGweExchangeModule
! -- GWT specific option block:
integer(I4B), pointer :: inewton => null() !< unneeded newton flag allows for mvt to be used here
integer(I4B), pointer :: iAdvScheme !< the advection scheme at the interface:
- !! 0 = upstream, 1 = central, 2 = TVD
+ !! 0 = upstream, 1 = central, 2 = TVD, 3 = UTVD
!
! -- Mover transport package
integer(I4B), pointer :: inmvt => null() !< unit number for mover transport (0 if off)
@@ -141,7 +142,7 @@ subroutine gweexchange_create(filename, name, id, m1_id, m2_id, input_mempath)
call exchange%allocate_scalars()
exchange%filename = filename
exchange%typename = 'GWE-GWE'
- exchange%iAdvScheme = 0
+ exchange%iAdvScheme = ADV_SCHEME_UPSTREAM
exchange%ixt3d = 1
!
! -- set gwemodel1
@@ -682,8 +683,8 @@ subroutine source_options(this, iout)
integer(I4B), intent(in) :: iout
! -- local
type(ExgGwegweParamFoundType) :: found
- character(len=LENVARNAME), dimension(3) :: adv_scheme = &
- &[character(len=LENVARNAME) :: 'UPSTREAM', 'CENTRAL', 'TVD']
+ character(len=LENVARNAME), dimension(4) :: adv_scheme = &
+ &[character(len=LENVARNAME) :: 'UPSTREAM', 'CENTRAL', 'TVD', 'UTVD']
character(len=LINELENGTH) :: mvt_fname
!
! -- update defaults with values sourced from input context
diff --git a/src/Exchange/exg-gwtgwt.f90 b/src/Exchange/exg-gwtgwt.f90
index 858ee653cd7..4f7877f397d 100644
--- a/src/Exchange/exg-gwtgwt.f90
+++ b/src/Exchange/exg-gwtgwt.f90
@@ -679,8 +679,8 @@ subroutine source_options(this, iout)
integer(I4B), intent(in) :: iout
! -- local
type(ExgGwtgwtParamFoundType) :: found
- character(len=LENVARNAME), dimension(3) :: adv_scheme = &
- &[character(len=LENVARNAME) :: 'UPSTREAM', 'CENTRAL', 'TVD']
+ character(len=LENVARNAME), dimension(4) :: adv_scheme = &
+ &[character(len=LENVARNAME) :: 'UPSTREAM', 'CENTRAL', 'TVD', 'UTVD']
character(len=LINELENGTH) :: mvt_fname
!
! -- update defaults with values sourced from input context
diff --git a/src/Model/Connection/GweGweConnection.f90 b/src/Model/Connection/GweGweConnection.f90
index d66aa5895ee..1c5ba428e32 100644
--- a/src/Model/Connection/GweGweConnection.f90
+++ b/src/Model/Connection/GweGweConnection.f90
@@ -264,9 +264,11 @@ subroutine setGridExtent(this)
hasCnd = this%gweModel%incnd > 0
if (hasAdv) then
- if (this%iIfaceAdvScheme == 2) then
+ if (this%iIfaceAdvScheme == ADV_SCHEME_TVD .or. &
+ this%iIfaceAdvScheme == ADV_SCHEME_UTVD) then
this%exg_stencil_depth = 2
- if (this%gweModel%adv%iadvwt == ADV_SCHEME_TVD) then
+ if (this%gweModel%adv%iadvwt == ADV_SCHEME_TVD .or. &
+ this%gweModel%adv%iadvwt == ADV_SCHEME_UTVD) then
this%int_stencil_depth = 2
end if
end if
diff --git a/src/Model/Connection/GwtGwtConnection.f90 b/src/Model/Connection/GwtGwtConnection.f90
index 9a4ee8fc2c9..cac6b44d14b 100644
--- a/src/Model/Connection/GwtGwtConnection.f90
+++ b/src/Model/Connection/GwtGwtConnection.f90
@@ -257,9 +257,11 @@ subroutine setGridExtent(this)
hasDsp = this%gwtModel%indsp > 0
if (hasAdv) then
- if (this%iIfaceAdvScheme == 2) then
+ if (this%iIfaceAdvScheme == ADV_SCHEME_TVD .or. &
+ this%iIfaceAdvScheme == ADV_SCHEME_UTVD) then
this%exg_stencil_depth = 2
- if (this%gwtModel%adv%iadvwt == ADV_SCHEME_TVD) then
+ if (this%gwtModel%adv%iadvwt == ADV_SCHEME_TVD .or. &
+ this%gwtModel%adv%iadvwt == ADV_SCHEME_UTVD) then
this%int_stencil_depth = 2
end if
end if
diff --git a/src/Model/Discretization/DisUtils.f90 b/src/Model/Discretization/DisUtils.f90
new file mode 100644
index 00000000000..cbf7ea7be7e
--- /dev/null
+++ b/src/Model/Discretization/DisUtils.f90
@@ -0,0 +1,99 @@
+module DisUtilsModule
+ use KindModule, only: DP, I4B
+ use ConstantsModule, only: DSAME, DONE
+ use BaseDisModule, only: DisBaseType
+
+ implicit none
+ private
+
+ public :: number_connected_faces
+ public :: cell_center
+ public :: node_distance
+
+contains
+
+ !> @brief Returns the number of connected faces for a given cell.
+ !!
+ !! This function computes the number of faces of cell `n` that are connected to neighboring cells
+ !! in the discretization. The value is determined from the connectivity information in the
+ !! connection arrays, and does not include boundary faces (faces not connected to another cell).
+ !<
+ function number_connected_faces(dis, n) result(connected_faces_count)
+ class(DisBaseType), intent(in) :: dis
+ integer(I4B), intent(in) :: n
+ integer(I4B) :: connected_faces_count
+
+ connected_faces_count = dis%con%ia(n + 1) - dis%con%ia(n) - 1
+ end function number_connected_faces
+
+ !> @brief Returns the vector distance from cell n to cell m.
+ !!
+ !! This function computes the vector from the center of cell `n` to the center of cell `m`
+ !! in the discretization. If the cells are directly connected, the vector is computed along
+ !! the connection direction, taking into account cell geometry and, if available, cell saturation.
+ !! If the cells are not directly connected (e.g., when using an extended stencil such as neighbors-of-neighbors),
+ !! the vector is simply the difference between their centroids: `d = centroid(m) - centroid(n)`.
+ !! The returned vector always points from cell `n` to cell `m`.
+ !<
+ function node_distance(dis, n, m) result(d)
+ !-- modules
+ use TspFmiModule, only: TspFmiType
+ ! -- return
+ real(DP), dimension(3) :: d
+ ! -- dummy
+ class(DisBaseType), intent(in) :: dis
+ integer(I4B), intent(in) :: n, m
+ ! -- local
+ real(DP) :: x_dir, y_dir, z_dir, length
+ integer(I4B) :: ipos, isympos, ihc
+ real(DP), dimension(3) :: xn, xm
+
+ ! -- Find the connection position (isympos) between cell n and cell m
+ isympos = -1
+ do ipos = dis%con%ia(n) + 1, dis%con%ia(n + 1) - 1
+ if (dis%con%ja(ipos) == m) then
+ isympos = dis%con%jas(ipos)
+ exit
+ end if
+ end do
+
+ ! -- if the connection is not found, then return the distance between the two nodes
+ ! -- This can happen when using an extended stencil (neighbours-of-neigbhours) to compute the gradients
+ if (isympos == -1) then
+ xn = cell_center(dis, n)
+ xm = cell_center(dis, m)
+
+ d = xm - xn
+ return
+ end if
+
+ ! -- Get the connection direction and length
+ ihc = dis%con%ihc(isympos)
+ call dis%connection_vector(n, m, .false., DONE, DONE, ihc, x_dir, &
+ y_dir, z_dir, length)
+
+ ! -- Compute the distance vector
+ d(1) = x_dir * length
+ d(2) = y_dir * length
+ d(3) = z_dir * length
+
+ end function node_distance
+
+ !> @brief Returns the center coordinates of a given cell.
+ !!
+ !! This function computes the center of cell `n` in the discretization.
+ !! The center is returned as a 3-element vector containing the x, y, and z coordinates.
+ !! The x and y coordinates are taken directly from the cell center arrays, while the z coordinate
+ !! is computed as the average of the cell's top and bottom elevations.
+ !<
+ function cell_center(dis, n) result(center)
+ ! -- dummy
+ class(DisBaseType), intent(in) :: dis
+ integer(I4B), intent(in) :: n
+ real(DP), dimension(3) :: center
+
+ center(1) = dis%xc(n)
+ center(2) = dis%yc(n)
+ center(3) = (dis%top(n) + dis%bot(n)) / 2.0_dp
+ end function cell_center
+end module DisUtilsModule
diff --git a/src/Model/TransportModel/Gradient/IGradient.f90 b/src/Model/TransportModel/Gradient/IGradient.f90
new file mode 100644
index 00000000000..6a3268a4485
--- /dev/null
+++ b/src/Model/TransportModel/Gradient/IGradient.f90
@@ -0,0 +1,38 @@
+module IGradient
+ use KindModule, only: DP, I4B
+
+ implicit none
+ private
+
+ public :: IGradientType
+
+ !> @brief Abstract interface for cell-based gradient computation.
+ !!
+ !! This module defines the abstract type `IGradientType`, which provides a deferred
+ !! interface for computing the gradient of a scalar field at a given cell.
+ !! Any concrete gradient implementation must extend this type and implement the `get` method.
+ !<
+ type, abstract :: IGradientType
+ contains
+ procedure(get_if), deferred :: get
+ end type IGradientType
+
+ abstract interface
+
+ function get_if(this, n, c) result(grad_c)
+ ! -- import
+ import IGradientType
+ import DP, I4B
+ ! -- dummy
+ class(IGradientType), target :: this
+ integer(I4B), intent(in) :: n
+ real(DP), dimension(:), intent(in) :: c
+ !-- return
+ real(DP), dimension(3) :: grad_c
+ end function
+
+ end interface
+
+contains
+
+end module IGradient
diff --git a/src/Model/TransportModel/Gradient/LeastSquaresGradient.f90 b/src/Model/TransportModel/Gradient/LeastSquaresGradient.f90
new file mode 100644
index 00000000000..5e78ec0cd41
--- /dev/null
+++ b/src/Model/TransportModel/Gradient/LeastSquaresGradient.f90
@@ -0,0 +1,156 @@
+module LeastSquaresGradientModule
+ use KindModule, only: DP, I4B
+ use ConstantsModule, only: DONE
+
+ Use IGradient
+ use BaseDisModule, only: DisBaseType
+ use PseudoInverseModule, only: pinv
+ use DisUtilsModule, only: number_connected_faces, node_distance
+
+ implicit none
+ private
+
+ public :: LeastSquaresGradientType
+
+ type Array2D
+ real(DP), dimension(:, :), allocatable :: data
+ end type Array2D
+
+ !> @brief Weighted least-squares gradient method for structured and unstructured grids.
+ !!
+ !! This class implements a least-squares gradient reconstruction for use on both structured and unstructured grids.
+ !! For each cell, it precomputes and caches a gradient reconstruction matrix using the Moore-Penrose pseudoinverse,
+ !! based on the geometry and connectivity of the mesh. The operator is created once during initialization
+ !! and can then be efficiently applied to any scalar field to compute the gradient in each cell.
+ !! The gradient can then be computed by multiplying the reconstruction matrix with the difference vector.
+ !! ∇ɸ = R * ∑(ɸ_i - ɸ_up), where i are the neighboring cells.
+ !!
+ !! - The gradient operator is constructed using normalized direction vectors between cell centers,
+ !! scaled by the inverse of the distance.
+ !! - The least-squares approach ensures robust gradients even for irregular or rank-deficient stencils.
+ !! - The operator is cached for each cell, so gradient computation is efficient for repeated queries.
+ !! - The class provides a `get` method to compute the gradient for any cell and scalar field.
+ !!
+ !! @note Boundary cells are not handled in a special manner. This may impact the quality of the gradient
+ !! near boundaries, especially if a cell does not have enough neighbors (fewer than three in 3D).
+ !<
+ type, extends(IGradientType) :: LeastSquaresGradientType
+ class(DisBaseType), pointer :: dis
+ type(Array2D), allocatable, dimension(:) :: R ! Gradient reconstruction matrix
+ contains
+ procedure :: get
+
+ procedure, private :: compute_cell_gradient
+ procedure, private :: create_gradient_reconstruction_matrix
+ end type LeastSquaresGradientType
+
+ interface LeastSquaresGradientType
+ module procedure Constructor
+ end interface LeastSquaresGradientType
+
+contains
+ function constructor(dis) Result(gradient)
+ ! --dummy
+ class(DisBaseType), pointer, intent(in) :: dis
+ !-- return
+ type(LeastSquaresGradientType) :: gradient
+ ! -- local
+ integer(I4B) :: n, nodes
+
+ gradient%dis => dis
+ nodes = dis%nodes
+
+ ! -- Compute the gradient rec
+ nodes = dis%nodes
+ allocate (gradient%R(dis%nodes))
+ do n = 1, nodes
+ gradient%R(n)%data = gradient%create_gradient_reconstruction_matrix(n)
+ end do
+ end function constructor
+
+ function create_gradient_reconstruction_matrix(this, n) result(R)
+ ! -- dummy
+ class(LeastSquaresGradientType) :: this
+ integer(I4B), intent(in) :: n ! Cell index for which to create the operator
+ real(DP), dimension(:, :), allocatable :: R ! The resulting gradient reconstruction matrix (3 x number_connections)
+ ! -- local
+ integer(I4B) :: number_connections ! Number of connected neighboring cells
+ integer(I4B) :: ipos, local_pos, m ! Loop indices and neighbor cell index
+ real(DP) :: length ! Distance between cell centers
+ real(DP), dimension(3) :: dnm ! Vector from cell n to neighbor m
+ real(DP), dimension(:, :), allocatable :: d ! Matrix of normalized direction vectors (number_connections x 3)
+ real(DP), dimension(:, :), allocatable :: grad_scale ! Diagonal scaling matrix (number_connections x number_connections),
+ ! where each diagonal entry is the inverse of the distance between
+
+ number_connections = number_connected_faces(this%dis, n)
+
+ allocate (d(number_connections, 3))
+ allocate (R(3, number_connections))
+ allocate (grad_scale(number_connections, number_connections))
+
+ grad_scale = 0
+ d = 0
+
+ ! Assemble the distance matrix
+ ! Handle the internal connections
+ local_pos = 1
+ do ipos = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
+ m = this%dis%con%ja(ipos)
+
+ dnm = node_distance(this%dis, n, m)
+ length = norm2(dnm)
+
+ d(local_pos, :) = dnm / length
+ grad_scale(local_pos, local_pos) = 1.0_dp / length
+
+ local_pos = local_pos + 1
+ end do
+
+ ! Compute the gradient reconstructions matrix
+ R = matmul(pinv(d), grad_scale)
+
+ end function create_gradient_reconstruction_matrix
+
+ function get(this, n, c) result(grad_c)
+ ! -- dummy
+ class(LeastSquaresGradientType), target :: this
+ integer(I4B), intent(in) :: n
+ real(DP), dimension(:), intent(in) :: c
+ !-- return
+ real(DP), dimension(3) :: grad_c
+
+ grad_c = this%compute_cell_gradient(n, c)
+ end function get
+
+ function compute_cell_gradient(this, n, phi_new) result(grad_c)
+ ! -- return
+ real(DP), dimension(3) :: grad_c
+ ! -- dummy
+ class(LeastSquaresGradientType), target :: this
+ integer(I4B), intent(in) :: n
+ real(DP), dimension(:), intent(in) :: phi_new
+ ! -- local
+ real(DP), dimension(:, :), pointer :: R
+ integer(I4B) :: ipos, local_pos
+ integer(I4B) :: number_connections
+
+ integer(I4B) :: m
+ real(DP), dimension(:), allocatable :: dc
+
+ ! Assemble the concentration difference vector
+ number_connections = number_connected_faces(this%dis, n)
+ allocate (dc(number_connections))
+ local_pos = 1
+ do ipos = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
+ m = this%dis%con%ja(ipos)
+ dc(local_pos) = phi_new(m) - phi_new(n)
+ local_pos = local_pos + 1
+ end do
+
+ ! Compute the cells gradient
+ R => this%R(n)%data
+ grad_c = matmul(R, dc)
+
+ end function compute_cell_gradient
+
+end module LeastSquaresGradientModule
diff --git a/src/Model/TransportModel/InterpolationScheme/AdvSchemeEnum.f90 b/src/Model/TransportModel/InterpolationScheme/AdvSchemeEnum.f90
index 3e6ee17b1b3..4efb6f3428b 100644
--- a/src/Model/TransportModel/InterpolationScheme/AdvSchemeEnum.f90
+++ b/src/Model/TransportModel/InterpolationScheme/AdvSchemeEnum.f90
@@ -7,5 +7,6 @@ module AdvSchemeEnumModule
integer(I4B), parameter :: ADV_SCHEME_UPSTREAM = 0
integer(I4B), parameter :: ADV_SCHEME_CENTRAL = 1
integer(I4B), parameter :: ADV_SCHEME_TVD = 2
+ integer(I4B), parameter :: ADV_SCHEME_UTVD = 3
-end module AdvSchemeEnumModule
\ No newline at end of file
+end module AdvSchemeEnumModule
diff --git a/src/Model/TransportModel/InterpolationScheme/UTVDScheme.f90 b/src/Model/TransportModel/InterpolationScheme/UTVDScheme.f90
new file mode 100644
index 00000000000..275dbd8ff08
--- /dev/null
+++ b/src/Model/TransportModel/InterpolationScheme/UTVDScheme.f90
@@ -0,0 +1,205 @@
+module UTVDSchemeModule
+ use KindModule, only: DP, I4B
+ use ConstantsModule, only: DONE, DZERO, DSAME, DHALF
+ use InterpolationSchemeInterfaceModule, only: InterpolationSchemeInterface, &
+ CoefficientsType
+ use BaseDisModule, only: DisBaseType
+ use TspFmiModule, only: TspFmiType
+ use IGradient, only: IGradientType
+ use DisUtilsModule, only: node_distance
+
+ implicit none
+ private
+
+ public :: UTVDSchemeType
+
+ !> @brief Total Variation Diminishing (TVD) interpolation scheme.
+ !!
+ !! This class implements a high-resolution, TVD interpolation scheme for use in transport modeling.
+ !! It extends a generic interpolation scheme interface and supports multiple TVD limiters (van Leer,
+ !! Koren, Superbee, van Albada, etc.) for controlling numerical diffusion and oscillations.
+ !! The default limiter is van Leer, but others can be selected by changing the `limiter_id` member.
+ !!
+ !! The scheme uses a combination of low-order upwind and high-order limited terms to compute
+ !! face concentrations. The high-order term is constructed using a gradient-based virtual node
+ !! value, following the approach described by Darwish et al. An additional TVD clamp is applied
+ !! to the virtual node value to enforce TVD compliance, especially on grids where the original
+ !! method may not guarantee monotonicity.
+ !!
+ !! - Supports both structured and unstructured grids via polymorphic discretization and gradient objects.
+ !! - The limiter can be selected via the `limiter_id` member (default is van Leer).
+ !! - The method `find_local_extrema` finds the minimum and maximum values among the current cell and its neighbors,
+ !! which is used to enforce the TVD condition.
+ !! - The `compute` method calculates the face coefficients for the transport equation.
+ !<
+ type, extends(InterpolationSchemeInterface) :: UTVDSchemeType
+ private
+ class(DisBaseType), pointer :: dis
+ type(TspFmiType), pointer :: fmi
+ class(IGradientType), pointer :: gradient
+ integer(I4B) :: limiter_id = 2 ! default to van Leer limiter
+ contains
+ procedure :: compute
+
+ procedure, private :: find_local_extrema
+ procedure, private :: limiter
+ end type UTVDSchemeType
+
+ interface UTVDSchemeType
+ module procedure constructor
+ end interface UTVDSchemeType
+
+contains
+ function constructor(dis, fmi, gradient) &
+ result(interpolation_scheme)
+ ! -- return
+ type(UTVDSchemeType) :: interpolation_scheme
+ ! --dummy
+ class(DisBaseType), pointer, intent(in) :: dis
+ type(TspFmiType), pointer, intent(in) :: fmi
+ class(IGradientType), allocatable, target, intent(in) :: gradient
+
+ interpolation_scheme%dis => dis
+ interpolation_scheme%fmi => fmi
+ interpolation_scheme%gradient => gradient
+
+ end function constructor
+
+ function compute(this, n, m, iposnm, phi) result(phi_face)
+ !-- return
+ type(CoefficientsType), target :: phi_face ! Output: coefficients for the face between cells n and m
+ ! -- dummy
+ class(UTVDSchemeType), target :: this
+ integer(I4B), intent(in) :: n ! Index of the first cell
+ integer(I4B), intent(in) :: m ! Index of the second cell
+ integer(I4B), intent(in) :: iposnm ! Position in the connectivity array for the n-m connection
+ real(DP), intent(in), dimension(:) :: phi ! Array of scalar values (e.g., concentrations) at all cells
+ ! -- local
+ integer(I4B) :: iup, idn, isympos ! Indices for upwind, downwind, and symmetric position in connectivity
+ real(DP) :: qnm ! Flow rate across the face between n and m
+ real(DP), pointer :: coef_up, coef_dn ! Pointers to upwind and downwind coefficients in phi_face
+ real(DP), dimension(3) :: grad_c ! gradient at upwind cell
+ real(DP), dimension(3) :: dnm ! vector from upwind to downwind cell
+ real(DP) :: smooth ! Smoothness indicator for limiter i.e. ratio of gradients
+ real(DP) :: alimiter ! Value of the TVD limiter
+ real(DP) :: cl1, cl2 ! Connection lengths from upwind and downwind cells to the face
+ real(DP) :: relative_distance ! Relative distance factor for high-order term
+ real(DP) :: c_virtual ! Virtual node concentration (Darwish method)
+ real(DP) :: min_phi, max_phi ! Local minimum and maximum among cell and neighbors
+
+ isympos = this%dis%con%jas(iposnm)
+ qnm = this%fmi%gwfflowja(iposnm)
+ !
+ ! -- Find upstream node
+ if (qnm > DZERO) then
+ ! -- positive flow into n means m is upstream
+ iup = m
+ idn = n
+
+ cl1 = this%dis%con%cl2(isympos)
+ cl2 = this%dis%con%cl1(isympos)
+
+ coef_up => phi_face%c_m
+ coef_dn => phi_face%c_n
+ else
+ iup = n
+ idn = m
+
+ cl1 = this%dis%con%cl1(isympos)
+ cl2 = this%dis%con%cl2(isympos)
+
+ coef_up => phi_face%c_n
+ coef_dn => phi_face%c_m
+ end if
+ !
+ ! -- Add low order terms
+ coef_up = DONE
+ !
+ ! -- Add high order terms
+ !
+ ! -- Return if straddled cells have same value
+ if (abs(phi(idn) - phi(iup)) < DSAME) return
+ !
+ ! -- Compute cell concentration gradient
+ grad_c = this%gradient%get(iup, phi)
+ !
+ ! Darwish's method to compute virtual node concentration
+ dnm = node_distance(this%dis, iup, idn)
+ c_virtual = phi(idn) - 2.0_dp * (dot_product(grad_c, dnm))
+ !
+ ! Enforce local TVD condition.
+ ! This is done by limiting the virtual concentration to the range of
+ ! the max and min concentration of the neighbouring cells.
+ call find_local_extrema(this, iup, phi, min_phi, max_phi)
+
+ if (c_virtual > max_phi) then
+ c_virtual = max_phi
+ end if
+
+ if (c_virtual < max(min_phi, DZERO)) then
+ c_virtual = max(min_phi, DZERO)
+ end if
+ !
+ ! -- Compute smoothness factor
+ smooth = (phi(iup) - c_virtual) / (phi(idn) - phi(iup))
+ !
+ ! -- Compute limiter
+ alimiter = this%limiter(smooth)
+
+ ! High order term is:
+ relative_distance = cl1 / (cl1 + cl2)
+ phi_face%rhs = -relative_distance * alimiter * (phi(idn) - phi(iup))
+
+ ! Alternative way of writing the high order term by adding it to the
+ ! coefficients matrix. The equation to be added is:
+ ! high_order = cl1 / (cl1 + cl2) * alimiter * qnm * (phi(idn) - phi(iup))
+ ! This is split into two parts:
+ ! coef_up = coef_up - relative_distance * alimiter
+ ! coef_dn = coef_dn + relative_distance * alimiter
+
+ end function compute
+
+ subroutine find_local_extrema(this, n, phi, min_phi, max_phi)
+ ! -- dummy
+ class(UTVDSchemeType) :: this
+ integer(I4B), intent(in) :: n
+ real(DP), intent(in), dimension(:) :: phi
+ real(DP), intent(out) :: min_phi, max_phi
+ ! -- local
+ integer(I4B) :: ipos, m
+
+ min_phi = phi(n)
+ max_phi = phi(n)
+ do ipos = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
+ m = this%dis%con%ja(ipos)
+ min_phi = min(min_phi, phi(m))
+ max_phi = max(max_phi, phi(m))
+ end do
+ end subroutine
+
+ function limiter(this, r) result(theta)
+ ! -- return
+ real(DP) :: theta ! limited slope
+ ! -- dummy
+ class(UTVDSchemeType) :: this
+ real(DP) :: r ! ratio of successive gradients
+
+ select case (this%limiter_id)
+ case (2) ! van Leer
+ theta = max(0.0_dp, min((r + dabs(r)) / (1.0_dp + dabs(r)), 2.0_dp))
+ case (3) ! Koren
+ theta = max(0.0_dp, min(2.0_dp * r, &
+ 1.0_dp / 3.0_dp + 2.0_dp / 3.0_dp * r, 2.0_dp))
+ case (4) ! Superbee
+ theta = max(0.0_dp, min(2.0_dp * r, 1.0_dp), min(r, 2.0_dp))
+ case (5) ! van Albada
+ theta = max(0.0_dp, (r * r + r) / (r * r + 1.0_dp))
+ case (6) ! Koren modified
+ theta = max(0.0_dp, min(4.0_dp * r * r + r, &
+ 1.0_dp / 3.0_dp + 2.0_dp / 3.0_dp * r, 2.0_dp))
+ case default
+ theta = DZERO
+ end select
+ end function
+
+end module UTVDSchemeModule
diff --git a/src/Model/TransportModel/tsp-adv.f90 b/src/Model/TransportModel/tsp-adv.f90
index 0184e46510f..e20262e7fc0 100644
--- a/src/Model/TransportModel/tsp-adv.f90
+++ b/src/Model/TransportModel/tsp-adv.f90
@@ -7,6 +7,9 @@ module TspAdvModule
use TspFmiModule, only: TspFmiType
use TspAdvOptionsModule, only: TspAdvOptionsType
use MatrixBaseModule, only: MatrixBaseType
+ ! -- Gradient schemes
+ use IGradient, only: IGradientType
+ use LeastSquaresGradientModule, only: LeastSquaresGradientType
! -- Interpolation schemes
use InterpolationSchemeInterfaceModule, only: InterpolationSchemeInterface, &
CoefficientsType
@@ -14,6 +17,7 @@ module TspAdvModule
use UpstreamSchemeModule, only: UpstreamSchemeType
use CentralDifferenceSchemeModule, only: CentralDifferenceSchemeType
use TVDSchemeModule, only: TVDSchemeType
+ use UTVDSchemeModule, only: UTVDSchemeType
implicit none
private
@@ -28,6 +32,7 @@ module TspAdvModule
real(DP), pointer :: eqnsclfac => null() !< governing equation scale factor; =1. for solute; =rhow*cpw for energy
class(InterpolationSchemeInterface), allocatable :: face_interpolation !< interpolation scheme for face values
+ class(IGradientType), allocatable :: gradient !< cell centered gradient
contains
procedure :: adv_df
@@ -133,6 +138,10 @@ subroutine adv_ar(this, dis, ibound)
case (ADV_SCHEME_TVD)
this%face_interpolation = &
TVDSchemeType(this%dis, this%fmi, this%ibound)
+ case (ADV_SCHEME_UTVD)
+ this%gradient = LeastSquaresGradientType(this%dis)
+ this%face_interpolation = &
+ UTVDSchemeType(this%dis, this%fmi, this%gradient)
case default
call store_error("Unknown advection scheme", terminate=.TRUE.)
end select
@@ -335,8 +344,8 @@ subroutine source_options(this)
! -- dummy
class(TspAdvType) :: this
! -- locals
- character(len=LENVARNAME), dimension(3) :: scheme = &
- &[character(len=LENVARNAME) :: 'UPSTREAM', 'CENTRAL', 'TVD']
+ character(len=LENVARNAME), dimension(4) :: scheme = &
+ &[character(len=LENVARNAME) :: 'UPSTREAM', 'CENTRAL', 'TVD', 'UTVD']
logical(LGP) :: found_scheme, found_atspercel
! -- formats
character(len=*), parameter :: fmtiadvwt = &
@@ -352,7 +361,7 @@ subroutine source_options(this)
! should currently be set to index of scheme names
if (this%iadvwt == 0) then
write (errmsg, '(a, a)') &
- 'Unknown scheme, must be "UPSTREAM", "CENTRAL" or "TVD"'
+ 'Unknown scheme, must be "UPSTREAM", "CENTRAL", "TVD" or "UTVD"'
call store_error(errmsg)
call store_error_filename(this%input_fname)
else
diff --git a/src/meson.build b/src/meson.build
index 72e9cc0390c..f6a15d9eb78 100644
--- a/src/meson.build
+++ b/src/meson.build
@@ -171,6 +171,7 @@ modflow_sources = files(
'Model' / 'Connection' / 'DistributedVariable.f90',
'Model' / 'Discretization' / 'Dis.f90',
'Model' / 'Discretization' / 'Dis2d.f90',
+ 'Model' / 'Discretization' / 'DisUtils.f90',
'Model' / 'Discretization' / 'Disu.f90',
'Model' / 'Discretization' / 'Disv.f90',
'Model' / 'Discretization' / 'Disv1d.f90',
@@ -277,11 +278,14 @@ modflow_sources = files(
'Model' / 'ModelUtilities' / 'VectorInterpolation.f90',
'Model' / 'ModelUtilities' / 'Xt3dAlgorithm.f90',
'Model' / 'ModelUtilities' / 'Xt3dInterface.f90',
+ 'Model' / 'TransportModel' / 'Gradient' / 'IGradient.f90',
+ 'Model' / 'TransportModel' / 'Gradient' / 'LeastSquaresGradient.f90',
'Model' / 'TransportModel' / 'InterpolationScheme' / 'AdvSchemeEnum.f90',
'Model' / 'TransportModel' / 'InterpolationScheme' / 'CentralDifferenceScheme.f90',
'Model' / 'TransportModel' / 'InterpolationScheme' / 'InterpolationSchemeInterface.f90',
'Model' / 'TransportModel' / 'InterpolationScheme' / 'TVDScheme.f90',
'Model' / 'TransportModel' / 'InterpolationScheme' / 'UpstreamScheme.f90',
+ 'Model' / 'TransportModel' / 'InterpolationScheme' / 'UTVDScheme.f90',
'Model' / 'TransportModel' / 'tsp.f90',
'Model' / 'TransportModel' / 'tsp-adv.f90',
'Model' / 'TransportModel' / 'tsp-apt.f90',