diff --git a/autotest/test_gwt_mst05.py b/autotest/test_gwt_mst05.py index ca5fa12a9dc..8294f1dcc50 100644 --- a/autotest/test_gwt_mst05.py +++ b/autotest/test_gwt_mst05.py @@ -276,9 +276,9 @@ def check_output(idx, test): assert False, f'could not load data from "{fpth}"' cnorm = obs["X008"] / 0.05 - cnorm_max = [0.32842034, 0.875391418] + cnorm_max = [0.324119806, 0.873678873] msg = f"{cnorm_max[idx]} /= {cnorm.max()}" - assert np.allclose(cnorm_max[idx], cnorm.max(), atol=0.001), msg + assert np.allclose(cnorm_max[idx], cnorm.max(), atol=1e-6), msg savefig = False if savefig: diff --git a/autotest/test_gwt_mst08.py b/autotest/test_gwt_mst08.py new file mode 100644 index 00000000000..046dafaa0af --- /dev/null +++ b/autotest/test_gwt_mst08.py @@ -0,0 +1,222 @@ +import os + +import flopy +import numpy as np +import pytest +from framework import TestFramework +from scipy import optimize + +cases = ["mst08a", "mst08b", "mst08c"] +sorption_type = ["linear", "freundlich", "langmuir"] + +# Test description +# +# This test builds and runs a single-cell with Sorption and a Source. +# The MST package is tested with several equilibrium isotherm. +# +# We inject a known total mass over a finite time and validate the +# end-of-period concentration against the **analytical mass-balance** +# solution: $$M_{total} = \theta V C - f_m \rho_b V \overline{C}$$. +# +# This solution doesn't have an exact closed form for all isotherms, +# so we solve it algebraically using a root-finding method. + +# Parameters +## Time parameters +nper = 1 +perlen = [10.0] +nstp = [10] +tsmult = [1.0] + +## Solver parameters +nouter, ninner = 100, 300 +hclose, rclose, relax = 1e-6, 1e-6, 1.0 + +## Model domain and grid definition +nlay, nrow, ncol = 1, 1, 1 +delr, delc = 1.0, 1.0 +top = 1.0 +botm = 0.0 + +## Mass source parameters +Msrc = 0.05 # mass source term (M/T) + +## MST parameters +porosity = 0.25 # porosity (dimensionless) +fm = 1 # fraction of mass in the solid phase (dimensionless) +rho_bulk = 1700.0 # bulk density of the solid phase (M/L^3) + +## Sorption parameters +### Linear sorption parameters +Kd = 5.0e-4 # linear sorption coefficient (L^3/M) +# Freundlich sorption parameters +Kf = 1.0e-4 # Freundlich sorption coefficient (L^3M^(-1) )**n +n_exp = 0.8 # Freundlich exponent (dimensionless) +### Langmuir sorption parameters +Kl = 5.0 # Langmuir sorption coefficient (L^3/M) +Sbar = 1e-4 # Langmuir maximum solid phase concentration + + +# Algebraic solution for equilibrium concentration +# solved using a root-finding method +def algebraic_solution(sorption_type): + Mtot = Msrc * perlen[0] # total mass added over the period (M) + V = delr * delc * (top - botm) # volume of the cell (L^3) + + C0 = Mtot / (porosity * V) + + # Find exact concentration at equilibrium for the listed isotherm + def C_equilibrium_linear(): + def f(C): + return porosity * V * C + fm * rho_bulk * Kd * C * V - Mtot + + def fprime(C): + return porosity * V + fm * rho_bulk * Kd * V + + c_exact = optimize.newton(f, C0, fprime=fprime) + return c_exact + + # Find exact concentration at equilibrium for the Freundlich isotherm + def C_equilibrium_freundlich(): + def f(C): + return porosity * V * C + fm * rho_bulk * Kf * C**n_exp * V - Mtot + + def fprime(C): + return porosity * V + fm * rho_bulk * Kf * n_exp * C ** (n_exp - 1) * V + + c_exact = optimize.newton(f, C0, fprime=fprime) + return c_exact + + # Find exact concentration at equilibrium for the Langmuir isotherm + def C_equilibrium_langmuir(): + def f(C): + return ( + porosity * V * C + + fm * rho_bulk * (Kl * Sbar * C) / (1 + Kl * C) * V + - Mtot + ) + + def fprime(C): + return porosity * V + fm * rho_bulk * (Kl * Sbar) / ((1 + Kl * C) ** 2) * V + + c_exact = optimize.newton(f, C0, fprime=fprime) + return c_exact + + if sorption_type == "linear": + c_exact = C_equilibrium_linear() + elif sorption_type == "freundlich": + c_exact = C_equilibrium_freundlich() + elif sorption_type == "langmuir": + c_exact = C_equilibrium_langmuir() + else: + raise ValueError(f"Unknown sorption type: {sorption_type}") + + return c_exact + + +def build_models(idx, test): + name = cases[idx] + sorption = sorption_type[idx] + + # build MODFLOW 6 files + ws = test.workspace + sim = flopy.mf6.MFSimulation(sim_name=name, sim_ws=ws, exe_name="mf6") + + # create tdis package + tdis_rc = [perlen[0], nstp[0], tsmult[0]] + tdis = flopy.mf6.ModflowTdis(sim, time_units="DAYS", nper=nper, perioddata=tdis_rc) + + # create gwt model + gwt = flopy.mf6.ModflowGwt(sim, modelname=name) + + # create iterative model solution + ims = flopy.mf6.ModflowIms( + sim, + print_option="SUMMARY", + outer_dvclose=hclose, + outer_maximum=nouter, + inner_maximum=ninner, + inner_dvclose=hclose, + rcloserecord=rclose, + linear_acceleration="BICGSTAB", + ) + + # create discretization package + dis = flopy.mf6.ModflowGwtdis( + gwt, nlay=nlay, nrow=nrow, ncol=ncol, delr=delr, delc=delc, top=top, botm=botm + ) + + # initial conditions + ic = flopy.mf6.ModflowGwtic(gwt, strt=0.0) + + # source term + spd = {0: [((0, 0, 0), Msrc)]} + src = flopy.mf6.ModflowGwtsrc(gwt, stress_period_data=spd) + + # create mst package + if sorption == "linear": + mst = flopy.mf6.ModflowGwtmst( + gwt, + sorption="linear", + porosity=porosity, + bulk_density=rho_bulk, + distcoef=Kd, + ) + elif sorption == "freundlich": + mst = flopy.mf6.ModflowGwtmst( + gwt, + sorption="freundlich", + porosity=porosity, + bulk_density=rho_bulk, + distcoef=Kf, + sp2=n_exp, + ) + elif sorption == "langmuir": + mst = flopy.mf6.ModflowGwtmst( + gwt, + sorption="langmuir", + porosity=porosity, + bulk_density=rho_bulk, + distcoef=Kl, + sp2=Sbar, + ) + else: + raise ValueError(f"Unknown sorption type: {sorption}") + + # output control + oc = flopy.mf6.ModflowGwtoc( + gwt, + budget_filerecord=f"{name}.cbc", + concentration_filerecord=f"{name}.ucn", + saverecord=[("CONCENTRATION", "ALL")], + printrecord=[("CONCENTRATION", "LAST"), ("BUDGET", "LAST")], + ) + + return sim, None + + +# Algebraic solution for equilibrium concentration + + +def check_output(idx, test): + sorption = sorption_type[idx] + c_exact = algebraic_solution(sorption) + + name = cases[idx] + fpth = os.path.join(test.workspace, f"{name}.ucn") + cobj = flopy.utils.HeadFile(fpth, precision="double", text="CONCENTRATION") + conc = cobj.get_data(idx=-1) + + assert np.isclose(conc[0, 0, 0], c_exact, rtol=1e-6), "Concentrations do not match!" + + +@pytest.mark.parametrize("idx, name", enumerate(cases)) +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() diff --git a/make/makefile b/make/makefile index 94f3a0eaeef..bce23aec7ee 100644 --- a/make/makefile +++ b/make/makefile @@ -41,19 +41,20 @@ SOURCEDIR34=../src/Utilities/Export SOURCEDIR35=../src/Utilities/Idm SOURCEDIR36=../src/Utilities/Idm/mf6blockfile SOURCEDIR37=../src/Utilities/Idm/netcdf -SOURCEDIR38=../src/Utilities/Libraries -SOURCEDIR39=../src/Utilities/Libraries/blas -SOURCEDIR40=../src/Utilities/Libraries/daglib -SOURCEDIR41=../src/Utilities/Libraries/rcm -SOURCEDIR42=../src/Utilities/Libraries/sparsekit -SOURCEDIR43=../src/Utilities/Libraries/sparskit2 -SOURCEDIR44=../src/Utilities/Matrix -SOURCEDIR45=../src/Utilities/Memory -SOURCEDIR46=../src/Utilities/Observation -SOURCEDIR47=../src/Utilities/OutputControl -SOURCEDIR48=../src/Utilities/Performance -SOURCEDIR49=../src/Utilities/TimeSeries -SOURCEDIR50=../src/Utilities/Vector +SOURCEDIR38=../src/Utilities/InterpolationScheme +SOURCEDIR39=../src/Utilities/Libraries +SOURCEDIR40=../src/Utilities/Libraries/blas +SOURCEDIR41=../src/Utilities/Libraries/daglib +SOURCEDIR42=../src/Utilities/Libraries/rcm +SOURCEDIR43=../src/Utilities/Libraries/sparsekit +SOURCEDIR44=../src/Utilities/Libraries/sparskit2 +SOURCEDIR45=../src/Utilities/Matrix +SOURCEDIR46=../src/Utilities/Memory +SOURCEDIR47=../src/Utilities/Observation +SOURCEDIR48=../src/Utilities/OutputControl +SOURCEDIR49=../src/Utilities/Performance +SOURCEDIR50=../src/Utilities/TimeSeries +SOURCEDIR51=../src/Utilities/Vector VPATH = \ ${SOURCEDIR1} \ @@ -105,7 +106,8 @@ ${SOURCEDIR46} \ ${SOURCEDIR47} \ ${SOURCEDIR48} \ ${SOURCEDIR49} \ -${SOURCEDIR50} +${SOURCEDIR50} \ +${SOURCEDIR51} .SUFFIXES: .f90 .F90 .o diff --git a/src/Model/GroundWaterTransport/gwt-mst.f90 b/src/Model/GroundWaterTransport/gwt-mst.f90 index fa39da39001..651d4f1eaad 100644 --- a/src/Model/GroundWaterTransport/gwt-mst.f90 +++ b/src/Model/GroundWaterTransport/gwt-mst.f90 @@ -340,14 +340,10 @@ subroutine mst_fc_srb(this, nodes, cold, nja, matrix_sln, idxglo, rhs, & real(DP) :: vcell real(DP) :: volfracm real(DP) :: rhobm - real(DP), dimension(nodes) :: c_half - real(DP) :: cbar_derv_half - real(DP) :: cbar_new, cbar_half, cbar_old - real(DP) :: sat_new, sat_half, sat_old + real(DP) :: sat_new, sat_old ! ! -- set variables tled = DONE / delt - c_half = 0.5_DP * (cold + cnew) ! ! -- loop through and calculate sorption contribution to hcof and rhs do n = 1, this%dis%nodes @@ -362,23 +358,24 @@ subroutine mst_fc_srb(this, nodes, cold, nja, matrix_sln, idxglo, rhs, & sat_new = this%fmi%gwfsat(n) sat_old = this%fmi%gwfsatold(n, delt) - ! -- Midpoint formulation using average values - cbar_new = this%isotherm%value(cnew, n) - cbar_old = this%isotherm%value(cold, n) - cbar_half = 0.5 * (cbar_new + cbar_old) - - sat_half = 0.5 * (sat_new + sat_old) - cbar_derv_half = this%isotherm%derivative(c_half, n) - - hhcof = -volfracm * rhobm * cbar_derv_half * sat_half * Vcell * tled + ! -- Matrix contribution for sorption term + hhcof = -volfracm * rhobm * sat_new * this%isotherm%derivative(cnew, n) & + * Vcell * tled idiag = this%dis%con%ia(n) call matrix_sln%add_value_pos(idxglo(idiag), hhcof) - rrhs = -volfracm * rhobm * cbar_derv_half * sat_half * cold(n) & - * Vcell * tled + ! -- Right-hand side contribution due to linearization + rrhs = -volfracm * rhobm * sat_new * this%isotherm%derivative(cnew, n) & + * cnew(n) * Vcell * tled rhs(n) = rhs(n) + rrhs - rrhs = volfracm * rhobm * cbar_half * (sat_new - sat_old) * Vcell * tled + rrhs = volfracm * rhobm * sat_new * this%isotherm%value(cnew, n) * & + Vcell * tled + rhs(n) = rhs(n) + rrhs + + ! -- Right-hand side contribution from previous time step + rrhs = -volfracm * rhobm * sat_old * this%isotherm%value(cold, n) * & + Vcell * tled rhs(n) = rhs(n) + rrhs end do @@ -639,14 +636,11 @@ subroutine mst_cq_srb(this, nodes, cnew, cold, flowja) real(DP) :: vcell real(DP) :: volfracm real(DP) :: rhobm - real(DP), dimension(nodes) :: c_half - real(DP) :: cbar_derv_half - real(DP) :: cbar_new, cbar_half, cbar_old - real(DP) :: sat_new, sat_half, sat_old + real(DP) :: sat_new, sat_old + real(DP) :: contribution ! ! -- initialize tled = DONE / delt - c_half = 0.5_DP * (cold + cnew) ! ! -- Calculate sorption change do n = 1, nodes @@ -666,20 +660,25 @@ subroutine mst_cq_srb(this, nodes, cnew, cold, flowja) sat_new = this%fmi%gwfsat(n) sat_old = this%fmi%gwfsatold(n, delt) - ! -- Midpoint formulation using average values - cbar_new = this%isotherm%value(cnew, n) - cbar_old = this%isotherm%value(cold, n) - cbar_half = 0.5 * (cbar_new + cbar_old) + ! -- Matrix contribution for sorption term + contribution = -volfracm * rhobm * sat_new & + * this%isotherm%derivative(cnew, n) * cnew(n) * Vcell * tled + rate = rate + contribution + + ! -- Right-hand side contribution due to linearization + ! -- Note: this contribution should cancel with the matrix contribution when the outer loop is converged + contribution = -volfracm * rhobm * sat_new * & + this%isotherm%derivative(cnew, n) * cnew(n) * Vcell * tled + rate = rate - contribution - sat_half = 0.5 * (sat_new + sat_old) - cbar_derv_half = this%isotherm%derivative(c_half, n) + contribution = volfracm * rhobm * sat_new * this%isotherm%value(cnew, n) & + * Vcell * tled + rate = rate - contribution - rate = -volfracm * rhobm * cbar_derv_half * sat_half * cnew(n) & - * Vcell * tled - rate = rate + volfracm * rhobm * cbar_derv_half * sat_half * cold(n) & - * Vcell * tled - rate = rate - volfracm * rhobm * cbar_half * (sat_new - sat_old) & - * Vcell * tled + ! -- Right-hand side contribution from previous time step + contribution = -volfracm * rhobm * sat_old * this%isotherm%value(cold, n) & + * Vcell * tled + rate = rate - contribution this%ratesrb(n) = rate idiag = this%dis%con%ia(n)