From 76a816cf418ccb622c90a92258de4d1b9e0d9239 Mon Sep 17 00:00:00 2001 From: Sunny Titus Date: Mon, 22 Dec 2025 17:19:52 +0100 Subject: [PATCH 01/12] Refactor Isotherm functionality. Split them of in separate classes. Use it in the sorption method --- make/makefile | 68 ++++--- msvs/mf6core.vfproj | 7 + src/Model/GroundWaterTransport/gwt-mst.f90 | 168 ++++++------------ .../Isotherm/FreundlichIsotherm.f90 | 61 +++++++ .../TransportModel/Isotherm/IsothermEnum.f90 | 11 ++ .../Isotherm/IsothermFactory.f90 | 40 +++++ .../Isotherm/IsothermInterface.f90 | 32 ++++ .../Isotherm/LangmuirIsotherm.f90 | 64 +++++++ .../Isotherm/LinearIsotherm.f90 | 49 +++++ src/meson.build | 6 + 10 files changed, 367 insertions(+), 139 deletions(-) create mode 100644 src/Model/TransportModel/Isotherm/FreundlichIsotherm.f90 create mode 100644 src/Model/TransportModel/Isotherm/IsothermEnum.f90 create mode 100644 src/Model/TransportModel/Isotherm/IsothermFactory.f90 create mode 100644 src/Model/TransportModel/Isotherm/IsothermInterface.f90 create mode 100644 src/Model/TransportModel/Isotherm/LangmuirIsotherm.f90 create mode 100644 src/Model/TransportModel/Isotherm/LinearIsotherm.f90 diff --git a/make/makefile b/make/makefile index 38ff4ad25c0..098d33da2d3 100644 --- a/make/makefile +++ b/make/makefile @@ -25,34 +25,36 @@ SOURCEDIR18=../src/Model/SurfaceWaterFlow SOURCEDIR19=../src/Model/TransportModel 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/Solution/ParticleTracker/Particle/Events -SOURCEDIR30=../src/Timing -SOURCEDIR31=../src/Utilities -SOURCEDIR32=../src/Utilities/ArrayRead -SOURCEDIR33=../src/Utilities/Export -SOURCEDIR34=../src/Utilities/Idm -SOURCEDIR35=../src/Utilities/Idm/mf6blockfile -SOURCEDIR36=../src/Utilities/Idm/netcdf -SOURCEDIR37=../src/Utilities/Libraries -SOURCEDIR38=../src/Utilities/Libraries/blas -SOURCEDIR39=../src/Utilities/Libraries/daglib -SOURCEDIR40=../src/Utilities/Libraries/rcm -SOURCEDIR41=../src/Utilities/Libraries/sparsekit -SOURCEDIR42=../src/Utilities/Libraries/sparskit2 -SOURCEDIR43=../src/Utilities/Matrix -SOURCEDIR44=../src/Utilities/Memory -SOURCEDIR45=../src/Utilities/Observation -SOURCEDIR46=../src/Utilities/OutputControl -SOURCEDIR47=../src/Utilities/Performance -SOURCEDIR48=../src/Utilities/TimeSeries -SOURCEDIR49=../src/Utilities/Vector +SOURCEDIR22=../src/Model/TransportModel/Isotherm +SOURCEDIR23=../src/Solution +SOURCEDIR24=../src/Solution/LinearMethods +SOURCEDIR25=../src/Solution/PETSc +SOURCEDIR26=../src/Solution/ParticleTracker +SOURCEDIR27=../src/Solution/ParticleTracker/Domain +SOURCEDIR28=../src/Solution/ParticleTracker/Method +SOURCEDIR29=../src/Solution/ParticleTracker/Particle +SOURCEDIR30=../src/Solution/ParticleTracker/Particle/Events +SOURCEDIR31=../src/Timing +SOURCEDIR32=../src/Utilities +SOURCEDIR33=../src/Utilities/ArrayRead +SOURCEDIR34=../src/Utilities/Export +SOURCEDIR35=../src/Utilities/Idm +SOURCEDIR36=../src/Utilities/Idm/mf6blockfile +SOURCEDIR37=../src/Utilities/Idm/netcdf +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} \ @@ -103,7 +105,9 @@ ${SOURCEDIR45} \ ${SOURCEDIR46} \ ${SOURCEDIR47} \ ${SOURCEDIR48} \ -${SOURCEDIR49} +${SOURCEDIR49} \ +${SOURCEDIR50} \ +${SOURCEDIR51} .SUFFIXES: .f90 .F90 .o @@ -356,6 +360,7 @@ $(OBJDIR)/SwfCxsUtils.o \ $(OBJDIR)/PrintSaveManager.o \ $(OBJDIR)/tsp-fmi.o \ $(OBJDIR)/SingularValueDecomposition.o \ +$(OBJDIR)/IsothermInterface.o \ $(OBJDIR)/SfrCrossSectionManager.o \ $(OBJDIR)/dag_module.o \ $(OBJDIR)/BoundaryPackageExt.o \ @@ -373,6 +378,10 @@ $(OBJDIR)/IGradient.o \ $(OBJDIR)/DisUtils.o \ $(OBJDIR)/PseudoInverse.o \ $(OBJDIR)/UzfEtUtil.o \ +$(OBJDIR)/LinearIsotherm.o \ +$(OBJDIR)/LangmuirIsotherm.o \ +$(OBJDIR)/IsothermEnum.o \ +$(OBJDIR)/FreundlichIsotherm.o \ $(OBJDIR)/Xt3dAlgorithm.o \ $(OBJDIR)/TvBase.o \ $(OBJDIR)/gwf-sfr.o \ @@ -402,6 +411,7 @@ $(OBJDIR)/CentralDifferenceScheme.o \ $(OBJDIR)/CachedGradient.o \ $(OBJDIR)/AdvSchemeEnum.o \ $(OBJDIR)/UzfCellGroup.o \ +$(OBJDIR)/IsothermFactory.o \ $(OBJDIR)/Xt3dInterface.o \ $(OBJDIR)/gwf-tvk.o \ $(OBJDIR)/SpdisWorkArray.o \ diff --git a/msvs/mf6core.vfproj b/msvs/mf6core.vfproj index c40e654f7d7..3049551b6c8 100644 --- a/msvs/mf6core.vfproj +++ b/msvs/mf6core.vfproj @@ -415,6 +415,13 @@ + + + + + + + diff --git a/src/Model/GroundWaterTransport/gwt-mst.f90 b/src/Model/GroundWaterTransport/gwt-mst.f90 index 0a5be24924b..fd02ed824e2 100644 --- a/src/Model/GroundWaterTransport/gwt-mst.f90 +++ b/src/Model/GroundWaterTransport/gwt-mst.f90 @@ -19,6 +19,9 @@ module GwtMstModule use NumericalPackageModule, only: NumericalPackageType use BaseDisModule, only: DisBaseType use TspFmiModule, only: TspFmiType + use IsothermInterfaceModule, only: IsothermType + use IsothermFactoryModule, only: create_isotherm + use IsothermEnumModule implicit none public :: GwtMstType @@ -35,10 +38,6 @@ module GwtMstModule ENUMERATOR :: DECAY_OFF = 0 !< Decay (or production) of mass inactive (default) ENUMERATOR :: DECAY_FIRST_ORDER = 1 !< First-order decay ENUMERATOR :: DECAY_ZERO_ORDER = 2 !< Zeroth-order decay - ENUMERATOR :: SORPTION_OFF = 0 !< Sorption is inactive (default) - ENUMERATOR :: SORPTION_LINEAR = 1 !< Linear sorption between aqueous and solid phases - ENUMERATOR :: SORPTION_FREUND = 2 !< Freundlich sorption between aqueous and solid phases - ENUMERATOR :: SORPTION_LANG = 3 !< Langmuir sorption between aqueous and solid phases END ENUM !> @ brief Mobile storage and transfer @@ -72,7 +71,7 @@ module GwtMstModule real(DP), dimension(:), pointer, contiguous :: ratesrb => null() !< rate of sorption real(DP), dimension(:), pointer, contiguous :: ratedcys => null() !< rate of sorbed mass decay real(DP), dimension(:), pointer, contiguous :: csrb => null() !< sorbate concentration - ! + class(IsothermType), pointer :: isotherm => null() !< pointer to isotherm object ! -- misc integer(I4B), dimension(:), pointer, contiguous :: ibound => null() !< pointer to model ibound type(TspFmiType), pointer :: fmi => null() !< pointer to fmi object @@ -167,6 +166,9 @@ subroutine mst_ar(this, dis, ibound) ! ! -- source the data block call this%source_data() + ! + ! -- Create isotherm object if sorption is active + this%isotherm => create_isotherm(this%isrb, this%distcoef, this%sp2) end subroutine mst_ar !> @ brief Fill coefficient method for package @@ -335,15 +337,17 @@ subroutine mst_fc_srb(this, nodes, cold, nja, matrix_sln, idxglo, rhs, & integer(I4B) :: n, idiag real(DP) :: tled real(DP) :: hhcof, rrhs - real(DP) :: swt, swtpdt real(DP) :: vcell - real(DP) :: const1 - real(DP) :: const2 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 ! ! -- 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 @@ -353,100 +357,33 @@ subroutine mst_fc_srb(this, nodes, cold, nja, matrix_sln, idxglo, rhs, & ! ! -- assign variables vcell = this%dis%area(n) * (this%dis%top(n) - this%dis%bot(n)) - swtpdt = this%fmi%gwfsat(n) - swt = this%fmi%gwfsatold(n, delt) - idiag = this%dis%con%ia(n) - const1 = this%distcoef(n) - const2 = 0. - if (this%isrb == SORPTION_FREUND .or. this%isrb == SORPTION_LANG) then - const2 = this%sp2(n) - end if volfracm = this%get_volfracm(n) rhobm = this%bulk_density(n) - call mst_srb_term(this%isrb, volfracm, rhobm, vcell, tled, cnew(n), & - cold(n), swtpdt, swt, const1, const2, & - hcofval=hhcof, rhsval=rrhs) - ! - ! -- Add hhcof to diagonal and rrhs to right-hand side + 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 + 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 rhs(n) = rhs(n) + rrhs - ! + + rrhs = volfracm * rhobm * cbar_half * (sat_new - sat_old) * Vcell * tled + rhs(n) = rhs(n) + rrhs + end do end subroutine mst_fc_srb - !> @ brief Calculate sorption terms - !! - !! Subroutine to calculate sorption terms - !< - subroutine mst_srb_term(isrb, volfracm, rhobm, vcell, tled, cnew, cold, & - swnew, swold, const1, const2, rate, hcofval, rhsval) - ! -- dummy - integer(I4B), intent(in) :: isrb !< sorption flag 1, 2, 3 are linear, freundlich, and langmuir - real(DP), intent(in) :: volfracm !< volume fraction of mobile domain (fhat_m) - real(DP), intent(in) :: rhobm !< bulk density of mobile domain (rhob_m) - real(DP), intent(in) :: vcell !< volume of cell - real(DP), intent(in) :: tled !< one over time step length - real(DP), intent(in) :: cnew !< concentration at end of this time step - real(DP), intent(in) :: cold !< concentration at end of last time step - real(DP), intent(in) :: swnew !< cell saturation at end of this time step - real(DP), intent(in) :: swold !< cell saturation at end of last time step - real(DP), intent(in) :: const1 !< distribution coefficient or freundlich or langmuir constant - real(DP), intent(in) :: const2 !< zero, freundlich exponent, or langmuir sorption sites - real(DP), intent(out), optional :: rate !< calculated sorption rate - real(DP), intent(out), optional :: hcofval !< diagonal contribution to solution coefficient matrix - real(DP), intent(out), optional :: rhsval !< contribution to solution right-hand-side - ! -- local - real(DP) :: term - real(DP) :: derv - real(DP) :: cbarnew - real(DP) :: cbarold - real(DP) :: cavg - real(DP) :: cbaravg - real(DP) :: swavg - ! - ! -- Calculate based on type of sorption - if (isrb == SORPTION_LINEAR) then - ! -- linear - term = -volfracm * rhobm * vcell * tled * const1 - if (present(hcofval)) hcofval = term * swnew - if (present(rhsval)) rhsval = term * swold * cold - if (present(rate)) rate = term * swnew * cnew - term * swold * cold - else - ! - ! -- calculate average aqueous concentration - cavg = DHALF * (cold + cnew) - ! - ! -- set values based on isotherm - select case (isrb) - case (SORPTION_FREUND) - ! -- freundlich - cbarnew = get_freundlich_conc(cnew, const1, const2) - cbarold = get_freundlich_conc(cold, const1, const2) - derv = get_freundlich_derivative(cavg, const1, const2) - case (SORPTION_LANG) - ! -- langmuir - cbarnew = get_langmuir_conc(cnew, const1, const2) - cbarold = get_langmuir_conc(cold, const1, const2) - derv = get_langmuir_derivative(cavg, const1, const2) - end select - ! - ! -- calculate hcof, rhs, and rate for freundlich and langmuir - term = -volfracm * rhobm * vcell * tled - cbaravg = (cbarold + cbarnew) * DHALF - swavg = (swnew + swold) * DHALF - if (present(hcofval)) then - hcofval = term * derv * swavg - end if - if (present(rhsval)) then - rhsval = term * derv * swavg * cold - term * cbaravg * (swnew - swold) - end if - if (present(rate)) then - rate = term * derv * swavg * (cnew - cold) & - + term * cbaravg * (swnew - swold) - end if - end if - end subroutine mst_srb_term - !> @ brief Fill sorption-decay coefficient method for package !! !! Method to calculate and fill sorption-decay coefficients for the package. @@ -708,12 +645,13 @@ subroutine mst_cq_srb(this, nodes, cnew, cold, flowja) integer(I4B) :: idiag real(DP) :: rate real(DP) :: tled - real(DP) :: swt, swtpdt real(DP) :: vcell real(DP) :: volfracm real(DP) :: rhobm - real(DP) :: const1 - real(DP) :: const2 + 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 ! ! -- initialize tled = DONE / delt @@ -723,24 +661,34 @@ subroutine mst_cq_srb(this, nodes, cnew, cold, flowja) ! ! -- initialize rates this%ratesrb(n) = DZERO + rate = 0.0_dp ! ! -- skip if transport inactive if (this%ibound(n) <= 0) cycle ! ! -- assign variables - vcell = this%dis%area(n) * (this%dis%top(n) - this%dis%bot(n)) - swtpdt = this%fmi%gwfsat(n) - swt = this%fmi%gwfsatold(n, delt) - volfracm = this%get_volfracm(n) + Vcell = this%dis%area(n) * (this%dis%top(n) - this%dis%bot(n)) + rhobm = this%bulk_density(n) - const1 = this%distcoef(n) - const2 = 0. - if (this%isrb == SORPTION_FREUND .or. this%isrb == SORPTION_LANG) then - const2 = this%sp2(n) - end if - call mst_srb_term(this%isrb, volfracm, rhobm, vcell, tled, cnew(n), & - cold(n), swtpdt, swt, const1, const2, & - rate=rate) + volfracm = this%get_volfracm(n) + 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) + + 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 + this%ratesrb(n) = rate idiag = this%dis%con%ia(n) flowja(idiag) = flowja(idiag) + rate diff --git a/src/Model/TransportModel/Isotherm/FreundlichIsotherm.f90 b/src/Model/TransportModel/Isotherm/FreundlichIsotherm.f90 new file mode 100644 index 00000000000..1d4bcb5a235 --- /dev/null +++ b/src/Model/TransportModel/Isotherm/FreundlichIsotherm.f90 @@ -0,0 +1,61 @@ +Module FreundlichIsothermModule + + use KindModule, only: DP, I4B + use IsothermInterfaceModule, only: IsothermType + + Implicit None + Private + Public :: FreundlichIsothermType + + type, extends(IsothermType) :: FreundlichIsothermType + real(DP), pointer, dimension(:) :: Kf => null() !< Freundlich constant + real(DP), pointer, dimension(:) :: a => null() !< Freundlich exponent + contains + procedure :: value + procedure :: derivative + end type FreundlichIsothermType + + interface FreundlichIsothermType + module procedure constructor + end interface FreundlichIsothermType + +contains + + function constructor(Kf, a) Result(isotherm) + type(FreundlichIsothermType) :: isotherm + ! -- dummy + real(DP), pointer, dimension(:), intent(in) :: Kf + real(DP), pointer, dimension(:), intent(in) :: a + ! -- local + isotherm%Kf => Kf + isotherm%a => a + + end function constructor + + function value(this, c, n) result(val) + class(FreundlichIsothermType), intent(in) :: this + real(DP), dimension(:), intent(in) :: c + integer(I4B), intent(in) :: n + real(DP) :: val + + if (c(n) > 0.0_DP) then + val = this%Kf(n) * c(n)**this%a(n) + else + val = 0.0_DP + end if + end function value + + function derivative(this, c, n) result(derv) + class(FreundlichIsothermType), intent(in) :: this + real(DP), dimension(:), intent(in) :: c + integer(I4B), intent(in) :: n + real(DP) :: derv + + if (c(n) > 0.0_DP) then + derv = this%a(n) * this%Kf(n) * c(n)**(this%a(n) - 1.0_DP) + else + derv = 0.0_DP + end if + end function derivative + +end module FreundlichIsothermModule diff --git a/src/Model/TransportModel/Isotherm/IsothermEnum.f90 b/src/Model/TransportModel/Isotherm/IsothermEnum.f90 new file mode 100644 index 00000000000..5df41802daa --- /dev/null +++ b/src/Model/TransportModel/Isotherm/IsothermEnum.f90 @@ -0,0 +1,11 @@ +module IsothermEnumModule + use KindModule, only: I4B + + implicit none + + integer(I4B), parameter :: SORPTION_OFF = 0 !< Sorption is inactive (default) + integer(I4B), parameter :: SORPTION_LINEAR = 1 !< Linear sorption between aqueous and solid phases + integer(I4B), parameter :: SORPTION_FREUND = 2 !< Freundlich sorption between aqueous and solid phases + integer(I4B), parameter :: SORPTION_LANG = 3 !< Langmuir sorption between aqueous and solid phases + +end module IsothermEnumModule diff --git a/src/Model/TransportModel/Isotherm/IsothermFactory.f90 b/src/Model/TransportModel/Isotherm/IsothermFactory.f90 new file mode 100644 index 00000000000..5406cf72895 --- /dev/null +++ b/src/Model/TransportModel/Isotherm/IsothermFactory.f90 @@ -0,0 +1,40 @@ +module IsothermFactoryModule + + use KindModule, only: I4B, DP + use SimModule, only: store_error + + use IsothermEnumModule + use IsothermInterfaceModule, only: IsothermType + use LinearIsothermsModule, only: LinearIsothermType + use FreundlichIsothermModule, only: FreundlichIsothermType + use LangmuirIsothermModule, only: LangmuirIsothermType + + implicit none + private + public :: create_isotherm + +contains + + function create_isotherm(isotherm_type, distcoef, sp2) result(isotherm) + ! -- result + class(IsothermType), pointer :: isotherm + ! -- dummy + integer(I4B), intent(in) :: isotherm_type + real(DP), dimension(:), pointer, contiguous :: distcoef + real(DP), dimension(:), pointer, contiguous :: sp2 + + select case (isotherm_type) + case (SORPTION_OFF) + nullify (isotherm) + case (SORPTION_LINEAR) + allocate (isotherm, source=LinearIsothermType(distcoef)) + case (SORPTION_FREUND) + allocate (isotherm, source=FreundlichIsothermType(distcoef, sp2)) + case (SORPTION_LANG) + allocate (isotherm, source=LangmuirIsothermType(distcoef, sp2)) + case default + call store_error('Sorption type not implemented yet.') + end select + + end function create_isotherm +end module IsothermFactoryModule diff --git a/src/Model/TransportModel/Isotherm/IsothermInterface.f90 b/src/Model/TransportModel/Isotherm/IsothermInterface.f90 new file mode 100644 index 00000000000..3a8bf9dd39e --- /dev/null +++ b/src/Model/TransportModel/Isotherm/IsothermInterface.f90 @@ -0,0 +1,32 @@ +module IsothermInterfaceModule + use KindModule, only: DP, I4B + + implicit none + private + public :: IsothermType + + type, abstract :: IsothermType + contains + procedure(value_interface), deferred :: value + procedure(derivative_interface), deferred :: derivative + end type IsothermType + + abstract interface + function value_interface(this, c, n) result(val) + import :: IsothermType, DP, I4B + class(IsothermType), intent(in) :: this + real(DP), dimension(:), intent(in) :: c + integer(I4B), intent(in) :: n + real(DP) :: val + end function value_interface + + function derivative_interface(this, c, n) result(derv) + import :: IsothermType, DP, I4B + class(IsothermType), intent(in) :: this + real(DP), dimension(:), intent(in) :: c + integer(I4B), intent(in) :: n + real(DP) :: derv + end function derivative_interface + end interface + +end module IsothermInterfaceModule diff --git a/src/Model/TransportModel/Isotherm/LangmuirIsotherm.f90 b/src/Model/TransportModel/Isotherm/LangmuirIsotherm.f90 new file mode 100644 index 00000000000..e00ffdbef75 --- /dev/null +++ b/src/Model/TransportModel/Isotherm/LangmuirIsotherm.f90 @@ -0,0 +1,64 @@ +module LangmuirIsothermModule + + use KindModule, only: DP, I4B + use IsothermInterfaceModule, only: IsothermType + + implicit none + private + public :: LangmuirIsothermType + + type, extends(IsothermType) :: LangmuirIsothermType + real(DP), pointer, dimension(:) :: Kl => null() !< Langmuir constant + real(DP), pointer, dimension(:) :: Sbar => null() !< Total concentration of sorption sites + contains + procedure :: value + procedure :: derivative + end type LangmuirIsothermType + + interface LangmuirIsothermType + module procedure constructor + end interface LangmuirIsothermType + +contains + function constructor(Kl, Sbar) Result(isotherm) + type(LangmuirIsothermType) :: isotherm + ! -- dummy + real(DP), pointer, dimension(:), intent(in) :: Kl + real(DP), pointer, dimension(:), intent(in) :: Sbar + ! -- local + isotherm%Kl => Kl + isotherm%Sbar => Sbar + + end function constructor + + function value(this, c, n) result(val) + class(LangmuirIsothermType), intent(in) :: this + real(DP), dimension(:), intent(in) :: c + integer(I4B), intent(in) :: n + real(DP) :: val + real(DP) :: denom + + if (c(n) > 0.0_DP) then + denom = 1.0_DP + this%Kl(n) * c(n) + val = (this%Sbar(n) * this%Kl(n) * c(n)) / denom + else + val = 0.0_DP + end if + end function value + + function derivative(this, c, n) result(derv) + class(LangmuirIsothermType), intent(in) :: this + real(DP), dimension(:), intent(in) :: c + integer(I4B), intent(in) :: n + real(DP) :: derv + real(DP) :: denom + + if (c(n) > 0.0_DP) then + denom = (1.0_DP + this%Kl(n) * c(n))**2.0_dp + derv = (this%Sbar(n) * this%Kl(n)) / denom + else + derv = 0.0_DP + end if + end function derivative + +end module LangmuirIsothermModule diff --git a/src/Model/TransportModel/Isotherm/LinearIsotherm.f90 b/src/Model/TransportModel/Isotherm/LinearIsotherm.f90 new file mode 100644 index 00000000000..770adbd59e4 --- /dev/null +++ b/src/Model/TransportModel/Isotherm/LinearIsotherm.f90 @@ -0,0 +1,49 @@ +module LinearIsothermsModule + + use KindModule, only: DP, I4B + use IsothermInterfaceModule, only: IsothermType + + implicit none + private + public :: LinearIsothermType + + type, extends(IsothermType) :: LinearIsothermType + real(DP), pointer, dimension(:) :: Kd => null() !< distribution coefficient + contains + procedure :: value + procedure :: derivative + end type LinearIsothermType + + interface LinearIsothermType + module procedure constructor + end interface LinearIsothermType + +contains + function constructor(Kd) Result(isotherm) + type(LinearIsothermType) :: isotherm + ! -- dummy + real(DP), pointer, dimension(:), intent(in) :: Kd + ! -- local + isotherm%Kd => Kd + + end function constructor + + function value(this, c, n) result(val) + class(LinearIsothermType), intent(in) :: this + real(DP), dimension(:), intent(in) :: c + integer(I4B), intent(in) :: n + real(DP) :: val + + val = this%Kd(n) * c(n) + end function value + + function derivative(this, c, n) result(derv) + class(LinearIsothermType), intent(in) :: this + real(DP), dimension(:), intent(in) :: c + integer(I4B), intent(in) :: n + real(DP) :: derv + + derv = this%Kd(n) + end function derivative + +end module LinearIsothermsModule diff --git a/src/meson.build b/src/meson.build index 4949b4e28bf..7b152e0c5ba 100644 --- a/src/meson.build +++ b/src/meson.build @@ -300,6 +300,12 @@ modflow_sources = files( 'Model' / 'TransportModel' / 'InterpolationScheme' / 'TVDScheme.f90', 'Model' / 'TransportModel' / 'InterpolationScheme' / 'UpstreamScheme.f90', 'Model' / 'TransportModel' / 'InterpolationScheme' / 'UTVDScheme.f90', + 'Model' / 'TransportModel' / 'Isotherm' / 'FreundlichIsotherm.f90', + 'Model' / 'TransportModel' / 'Isotherm' / 'IsothermEnum.f90', + 'Model' / 'TransportModel' / 'Isotherm' / 'IsothermFactory.f90', + 'Model' / 'TransportModel' / 'Isotherm' / 'IsothermInterface.f90', + 'Model' / 'TransportModel' / 'Isotherm' / 'LangmuirIsotherm.f90', + 'Model' / 'TransportModel' / 'Isotherm' / 'LinearIsotherm.f90', 'Model' / 'TransportModel' / 'tsp.f90', 'Model' / 'TransportModel' / 'tsp-adv.f90', 'Model' / 'TransportModel' / 'tsp-apt.f90', From bae906ee16062b18f12d2850dfb7d2ec0dbca480 Mon Sep 17 00:00:00 2001 From: Sunny Titus Date: Mon, 22 Dec 2025 17:49:05 +0100 Subject: [PATCH 02/12] Fix unitialized variable issue --- src/Model/GroundWaterTransport/gwt-mst.f90 | 1 + 1 file changed, 1 insertion(+) diff --git a/src/Model/GroundWaterTransport/gwt-mst.f90 b/src/Model/GroundWaterTransport/gwt-mst.f90 index fd02ed824e2..37dbb7dd753 100644 --- a/src/Model/GroundWaterTransport/gwt-mst.f90 +++ b/src/Model/GroundWaterTransport/gwt-mst.f90 @@ -655,6 +655,7 @@ subroutine mst_cq_srb(this, nodes, cnew, cold, flowja) ! ! -- initialize tled = DONE / delt + c_half = 0.5_DP * (cold + cnew) ! ! -- Calculate sorption change do n = 1, nodes From a9c0aa55e10fee0f34c85c1af3b373d588ae4687 Mon Sep 17 00:00:00 2001 From: Sunny Titus Date: Mon, 22 Dec 2025 19:13:51 +0100 Subject: [PATCH 03/12] Remove unused method. Move methd that is only used by IST to that module --- src/Model/GroundWaterTransport/gwt-ist.f90 | 39 +++++++++++- src/Model/GroundWaterTransport/gwt-mst.f90 | 72 ---------------------- 2 files changed, 36 insertions(+), 75 deletions(-) diff --git a/src/Model/GroundWaterTransport/gwt-ist.f90 b/src/Model/GroundWaterTransport/gwt-ist.f90 index dd07e3a5288..0d26b55bfcc 100644 --- a/src/Model/GroundWaterTransport/gwt-ist.f90 +++ b/src/Model/GroundWaterTransport/gwt-ist.f90 @@ -21,9 +21,8 @@ module GwtIstModule use BudgetModule, only: BudgetType use TspFmiModule, only: TspFmiType use GwtMstModule, only: GwtMstType, get_zero_order_decay, & - get_freundlich_conc, get_freundlich_derivative, & - get_langmuir_conc, get_langmuir_derivative, & - get_freundlich_kd, get_langmuir_kd + get_freundlich_conc, & + get_langmuir_conc use OutputControlDataModule, only: OutputControlDataType use MatrixBaseModule ! @@ -1586,4 +1585,38 @@ subroutine accumulate_budterm(budterm, ddterm, cimnew, cimold, cnew, idcy) ! end subroutine accumulate_budterm + !> @ brief Get effective Freundlich distribution coefficient + !< + function get_freundlich_kd(conc, kf, a) result(kd) + ! -- dummy + real(DP), intent(in) :: conc !< solute concentration + real(DP), intent(in) :: kf !< freundlich constant + real(DP), intent(in) :: a !< freundlich exponent + ! -- return + real(DP) :: kd !< effective distribution coefficient + ! + if (conc > DZERO) then + kd = kf * conc**(a - DONE) + else + kd = DZERO + end if + end function get_freundlich_kd + + !> @ brief Get effective Langmuir distribution coefficient + !< + function get_langmuir_kd(conc, kl, sbar) result(kd) + ! -- dummy + real(DP), intent(in) :: conc !< solute concentration + real(DP), intent(in) :: kl !< langmuir constant + real(DP), intent(in) :: sbar !< langmuir sorption sites + ! -- return + real(DP) :: kd !< effective distribution coefficient + ! + if (conc > DZERO) then + kd = (kl * sbar) / (DONE + kl * conc) + else + kd = DZERO + end if + end function get_langmuir_kd + end module GwtIstModule diff --git a/src/Model/GroundWaterTransport/gwt-mst.f90 b/src/Model/GroundWaterTransport/gwt-mst.f90 index 37dbb7dd753..9a2a7c8cdfe 100644 --- a/src/Model/GroundWaterTransport/gwt-mst.f90 +++ b/src/Model/GroundWaterTransport/gwt-mst.f90 @@ -1496,78 +1496,6 @@ function get_langmuir_conc(conc, kl, sbar) result(cbar) end if end function - !> @ brief Calculate sorption derivative using Freundlich - !! - !! Function to calculate sorption derivative using Freundlich - !< - function get_freundlich_derivative(conc, kf, a) result(derv) - ! -- dummy - real(DP), intent(in) :: conc !< solute concentration - real(DP), intent(in) :: kf !< freundlich constant - real(DP), intent(in) :: a !< freundlich exponent - ! -- return - real(DP) :: derv - ! - if (conc > DZERO) then - derv = kf * a * conc**(a - DONE) - else - derv = DZERO - end if - end function - - !> @ brief Calculate sorption derivative using Langmuir - !! - !! Function to calculate sorption derivative using Langmuir - !< - function get_langmuir_derivative(conc, kl, sbar) result(derv) - ! -- dummy - real(DP), intent(in) :: conc !< solute concentration - real(DP), intent(in) :: kl !< langmuir constant - real(DP), intent(in) :: sbar !< langmuir sorption sites - ! -- return - real(DP) :: derv - ! - if (conc > DZERO) then - derv = (kl * sbar) / (DONE + kl * conc)**DTWO - else - derv = DZERO - end if - end function - - !> @ brief Get effective Freundlich distribution coefficient - !< - function get_freundlich_kd(conc, kf, a) result(kd) - ! -- dummy - real(DP), intent(in) :: conc !< solute concentration - real(DP), intent(in) :: kf !< freundlich constant - real(DP), intent(in) :: a !< freundlich exponent - ! -- return - real(DP) :: kd !< effective distribution coefficient - ! - if (conc > DZERO) then - kd = kf * conc**(a - DONE) - else - kd = DZERO - end if - end function get_freundlich_kd - - !> @ brief Get effective Langmuir distribution coefficient - !< - function get_langmuir_kd(conc, kl, sbar) result(kd) - ! -- dummy - real(DP), intent(in) :: conc !< solute concentration - real(DP), intent(in) :: kl !< langmuir constant - real(DP), intent(in) :: sbar !< langmuir sorption sites - ! -- return - real(DP) :: kd !< effective distribution coefficient - ! - if (conc > DZERO) then - kd = (kl * sbar) / (DONE + kl * conc) - else - kd = DZERO - end if - end function get_langmuir_kd - !> @ brief Calculate zero-order decay rate and constrain if necessary !! !! Function to calculate the zero-order decay rate from the user specified From 3197e20940556a6b26a1c29c00bc9c5f98a342f4 Mon Sep 17 00:00:00 2001 From: Sunny Titus Date: Tue, 23 Dec 2025 13:28:22 +0100 Subject: [PATCH 04/12] Replace all occurences of the isotherm in in the MST package by the new classes. Refactor the IST package to use the new isotherm classes --- src/Model/GroundWaterTransport/gwt-ist.f90 | 69 ++++++------------ src/Model/GroundWaterTransport/gwt-mst.f90 | 84 +++------------------- 2 files changed, 32 insertions(+), 121 deletions(-) diff --git a/src/Model/GroundWaterTransport/gwt-ist.f90 b/src/Model/GroundWaterTransport/gwt-ist.f90 index 0d26b55bfcc..4e320cab31c 100644 --- a/src/Model/GroundWaterTransport/gwt-ist.f90 +++ b/src/Model/GroundWaterTransport/gwt-ist.f90 @@ -20,11 +20,12 @@ module GwtIstModule use BndModule, only: BndType use BudgetModule, only: BudgetType use TspFmiModule, only: TspFmiType - use GwtMstModule, only: GwtMstType, get_zero_order_decay, & - get_freundlich_conc, & - get_langmuir_conc + use GwtMstModule, only: GwtMstType, get_zero_order_decay use OutputControlDataModule, only: OutputControlDataType use MatrixBaseModule + use IsothermInterfaceModule, only: IsothermType + use IsothermFactoryModule, only: create_isotherm + use IsothermEnumModule ! implicit none ! @@ -79,6 +80,7 @@ module GwtIstModule real(DP), pointer, contiguous :: decay_sorbed(:) => null() !< first or zero order rate constant for sorbed mass real(DP), pointer, contiguous :: strg(:) => null() !< mass transfer rate real(DP) :: budterm(2, NBDITEMS) !< immmobile domain mass summaries + class(IsothermType), pointer :: isotherm => null() !< pointer to isotherm object contains @@ -201,6 +203,9 @@ subroutine ist_ar(this) call this%budget%budget_df(NBDITEMS, 'MASS', 'M', bdzone=this%packName) call this%budget%set_ibudcsv(this%ibudcsv) ! + ! -- Create isotherm object if sorption is active + this%isotherm => create_isotherm(this%isrb, this%distcoef, this%sp2) + ! ! -- Perform a check to ensure that sorption and decay are set ! consistently between the MST and IST packages. if (this%idcy /= this%mst%idcy) then @@ -347,33 +352,22 @@ subroutine ist_fc(this, rhs, ia, idxglo, matrix_sln) rhobim = this%bulk_density(n) ! set isotherm dependent sorption variables + cimsrbnew = this%isotherm%value(this%cimnew, n) + cimsrbold = this%isotherm%value(this%cimold, n) select case (this%isrb) - case (1) - ! linear + case (SORPTION_LINEAR) kdnew = this%distcoef(n) kdold = this%distcoef(n) - cimsrbnew = this%cimnew(n) * kdnew - cimsrbold = this%cimold(n) * kdold - case (2) - ! freundlich + case (SORPTION_FREUND) kdnew = get_freundlich_kd(this%cimnew(n), this%distcoef(n), & this%sp2(n)) kdold = get_freundlich_kd(this%cimold(n), this%distcoef(n), & this%sp2(n)) - cimsrbnew = get_freundlich_conc(this%cimnew(n), this%distcoef(n), & - this%sp2(n)) - cimsrbold = get_freundlich_conc(this%cimold(n), this%distcoef(n), & - this%sp2(n)) - case (3) - ! langmuir + case (SORPTION_LANG) kdnew = get_langmuir_kd(this%cimnew(n), this%distcoef(n), & this%sp2(n)) kdold = get_langmuir_kd(this%cimold(n), this%distcoef(n), & this%sp2(n)) - cimsrbnew = get_langmuir_conc(this%cimnew(n), this%distcoef(n), & - this%sp2(n)) - cimsrbold = get_langmuir_conc(this%cimold(n), this%distcoef(n), & - this%sp2(n)) end select ! set decay of sorbed solute parameters @@ -485,33 +479,22 @@ subroutine ist_cq(this, x, flowja, iadv) rhobim = this%bulk_density(n) ! set isotherm dependent sorption variables + cimsrbnew = this%isotherm%value(this%cimnew, n) + cimsrbold = this%isotherm%value(this%cimold, n) select case (this%isrb) - case (1) - ! linear + case (SORPTION_LINEAR) kdnew = this%distcoef(n) kdold = this%distcoef(n) - cimsrbnew = this%cimnew(n) * kdnew - cimsrbold = this%cimold(n) * kdold - case (2) - ! freundlich + case (SORPTION_FREUND) kdnew = get_freundlich_kd(this%cimnew(n), this%distcoef(n), & this%sp2(n)) kdold = get_freundlich_kd(this%cimold(n), this%distcoef(n), & this%sp2(n)) - cimsrbnew = get_freundlich_conc(this%cimnew(n), this%distcoef(n), & - this%sp2(n)) - cimsrbold = get_freundlich_conc(this%cimold(n), this%distcoef(n), & - this%sp2(n)) - case (3) - ! langmuir + case (SORPTION_LANG) kdnew = get_langmuir_kd(this%cimnew(n), this%distcoef(n), & this%sp2(n)) kdold = get_langmuir_kd(this%cimold(n), this%distcoef(n), & this%sp2(n)) - cimsrbnew = get_langmuir_conc(this%cimnew(n), this%distcoef(n), & - this%sp2(n)) - cimsrbold = get_langmuir_conc(this%cimold(n), this%distcoef(n), & - this%sp2(n)) end select ! set decay of sorbed solute parameters @@ -568,21 +551,13 @@ subroutine ist_calc_csrb(this, cim) real(DP), intent(in), dimension(:) :: cim !< immobile domain aqueous concentration at end of this time step ! -- local integer(I4B) :: n - real(DP) :: distcoef real(DP) :: csrb ! Calculate sorbed concentration do n = 1, size(cim) csrb = DZERO if (this%ibound(n) > 0) then - distcoef = this%distcoef(n) - if (this%isrb == 1) then - csrb = cim(n) * distcoef - else if (this%isrb == 2) then - csrb = get_freundlich_conc(cim(n), distcoef, this%sp2(n)) - else if (this%isrb == 3) then - csrb = get_langmuir_conc(cim(n), distcoef, this%sp2(n)) - end if + csrb = this%isotherm%value(cim, n) end if this%cimsrb(n) = csrb end do @@ -1117,11 +1092,11 @@ subroutine log_options(this, found, cim6_fname, budget_fname, & end if if (found%sorption) then select case (this%isrb) - case (1) + case (SORPTION_LINEAR) write (this%iout, fmtlinear) - case (2) + case (SORPTION_FREUND) write (this%iout, fmtfreundlich) - case (3) + case (SORPTION_LANG) write (this%iout, fmtlangmuir) end select end if diff --git a/src/Model/GroundWaterTransport/gwt-mst.f90 b/src/Model/GroundWaterTransport/gwt-mst.f90 index 9a2a7c8cdfe..111d890fb53 100644 --- a/src/Model/GroundWaterTransport/gwt-mst.f90 +++ b/src/Model/GroundWaterTransport/gwt-mst.f90 @@ -447,12 +447,12 @@ subroutine mst_fc_dcy_srb(this, nodes, cold, nja, matrix_sln, idxglo, & case (SORPTION_FREUND) ! ! -- nonlinear Freundlich sorption, so add to RHS - csrb = get_freundlich_conc(cnew(n), distcoef, this%sp2(n)) + csrb = this%isotherm%value(cnew, n) rrhs = term * csrb case (SORPTION_LANG) ! ! -- nonlinear Lanmuir sorption, so add to RHS - csrb = get_langmuir_conc(cnew(n), distcoef, this%sp2(n)) + csrb = this%isotherm%value(cnew, n) rrhs = term * csrb end select case (DECAY_ZERO_ORDER) @@ -460,17 +460,8 @@ subroutine mst_fc_dcy_srb(this, nodes, cold, nja, matrix_sln, idxglo, & ! -- call function to get zero-order decay rate, which may be changed ! from the user-specified rate to prevent negative concentrations if (distcoef > DZERO) then - select case (this%isrb) - case (SORPTION_LINEAR) - csrbold = cold(n) * distcoef - csrbnew = cnew(n) * distcoef - case (SORPTION_FREUND) - csrbold = get_freundlich_conc(cold(n), distcoef, this%sp2(n)) - csrbnew = get_freundlich_conc(cnew(n), distcoef, this%sp2(n)) - case (SORPTION_LANG) - csrbold = get_langmuir_conc(cold(n), distcoef, this%sp2(n)) - csrbnew = get_langmuir_conc(cnew(n), distcoef, this%sp2(n)) - end select + csrbold = this%isotherm%value(cold, n) + csrbnew = this%isotherm%value(cnew, n) ! decay_rate = get_zero_order_decay(this%decay_sorbed(n), & this%decayslast(n), & @@ -762,12 +753,12 @@ subroutine mst_cq_dcy_srb(this, nodes, cnew, cold, flowja) case (SORPTION_FREUND) ! ! -- nonlinear Freundlich sorption, so add to RHS - csrb = get_freundlich_conc(cnew(n), distcoef, this%sp2(n)) + csrb = this%isotherm%value(cnew, n) rrhs = term * csrb case (SORPTION_LANG) ! ! -- nonlinear Lanmuir sorption, so add to RHS - csrb = get_langmuir_conc(cnew(n), distcoef, this%sp2(n)) + csrb = this%isotherm%value(cnew, n) rrhs = term * csrb end select case (DECAY_ZERO_ORDER) @@ -775,17 +766,9 @@ subroutine mst_cq_dcy_srb(this, nodes, cnew, cold, flowja) ! -- Call function to get zero-order decay rate, which may be changed ! from the user-specified rate to prevent negative concentrations if (distcoef > DZERO) then - select case (this%isrb) - case (SORPTION_LINEAR) - csrbold = cold(n) * distcoef - csrbnew = cnew(n) * distcoef - case (SORPTION_FREUND) - csrbold = get_freundlich_conc(cold(n), distcoef, this%sp2(n)) - csrbnew = get_freundlich_conc(cnew(n), distcoef, this%sp2(n)) - case (SORPTION_LANG) - csrbold = get_langmuir_conc(cold(n), distcoef, this%sp2(n)) - csrbnew = get_langmuir_conc(cnew(n), distcoef, this%sp2(n)) - end select + csrbold = this%isotherm%value(cold, n) + csrbnew = this%isotherm%value(cnew, n) + decay_rate = get_zero_order_decay(this%decay_sorbed(n), & this%decayslast(n), & 0, csrbold, csrbnew, delt) @@ -810,22 +793,13 @@ subroutine mst_calc_csrb(this, cnew) real(DP), intent(in), dimension(:) :: cnew !< concentration at end of this time step ! -- local integer(I4B) :: n - real(DP) :: distcoef real(DP) :: csrb ! Calculate sorbed concentration do n = 1, size(cnew) csrb = DZERO if (this%ibound(n) > 0) then - distcoef = this%distcoef(n) - select case (this%isrb) - case (SORPTION_LINEAR) - csrb = cnew(n) * distcoef - case (SORPTION_FREUND) - csrb = get_freundlich_conc(cnew(n), distcoef, this%sp2(n)) - case (SORPTION_LANG) - csrb = get_langmuir_conc(cnew(n), distcoef, this%sp2(n)) - end select + csrb = this%isotherm%value(cnew, n) end if this%csrb(n) = csrb end do @@ -1458,44 +1432,6 @@ function get_volfracm(this, node) result(volfracm) volfracm = DONE - this%volfracim(node) end function get_volfracm - !> @ brief Calculate sorption concentration using Freundlich - !! - !! Function to calculate sorption concentration using Freundlich - !< - function get_freundlich_conc(conc, kf, a) result(cbar) - ! -- dummy - real(DP), intent(in) :: conc !< solute concentration - real(DP), intent(in) :: kf !< freundlich constant - real(DP), intent(in) :: a !< freundlich exponent - ! -- return - real(DP) :: cbar - ! - if (conc > DZERO) then - cbar = kf * conc**a - else - cbar = DZERO - end if - end function - - !> @ brief Calculate sorption concentration using Langmuir - !! - !! Function to calculate sorption concentration using Langmuir - !< - function get_langmuir_conc(conc, kl, sbar) result(cbar) - ! -- dummy - real(DP), intent(in) :: conc !< solute concentration - real(DP), intent(in) :: kl !< langmuir constant - real(DP), intent(in) :: sbar !< langmuir sorption sites - ! -- return - real(DP) :: cbar - ! - if (conc > DZERO) then - cbar = (kl * sbar * conc) / (DONE + kl * conc) - else - cbar = DZERO - end if - end function - !> @ brief Calculate zero-order decay rate and constrain if necessary !! !! Function to calculate the zero-order decay rate from the user specified From 09c23a772de04d4ca538c6b777ae3544ef69086b Mon Sep 17 00:00:00 2001 From: Sunny Titus Date: Tue, 23 Dec 2025 14:36:53 +0100 Subject: [PATCH 05/12] Add documentation to the isotherm files. Deallocate the isothem on finalization --- src/Model/GroundWaterTransport/gwt-ist.f90 | 4 ++ src/Model/GroundWaterTransport/gwt-mst.f90 | 6 +++ .../Isotherm/FreundlichIsotherm.f90 | 26 +++++++++--- .../Isotherm/IsothermFactory.f90 | 16 +++++--- .../Isotherm/IsothermInterface.f90 | 40 ++++++++++++------- .../Isotherm/LangmuirIsotherm.f90 | 34 ++++++++++++---- .../Isotherm/LinearIsotherm.f90 | 31 ++++++++++---- 7 files changed, 117 insertions(+), 40 deletions(-) diff --git a/src/Model/GroundWaterTransport/gwt-ist.f90 b/src/Model/GroundWaterTransport/gwt-ist.f90 index 4e320cab31c..5f3f3d93b8f 100644 --- a/src/Model/GroundWaterTransport/gwt-ist.f90 +++ b/src/Model/GroundWaterTransport/gwt-ist.f90 @@ -815,6 +815,10 @@ subroutine ist_da(this) deallocate (this%budget) call this%ocd%ocd_da() deallocate (this%ocd) + if (associated(this%isotherm)) then + deallocate (this%isotherm) + nullify (this%isotherm) + end if ! ! -- deallocate parent call this%BndType%bnd_da() diff --git a/src/Model/GroundWaterTransport/gwt-mst.f90 b/src/Model/GroundWaterTransport/gwt-mst.f90 index 111d890fb53..4e1650c0aeb 100644 --- a/src/Model/GroundWaterTransport/gwt-mst.f90 +++ b/src/Model/GroundWaterTransport/gwt-mst.f90 @@ -974,6 +974,12 @@ subroutine mst_da(this) ! ! -- Scalars ! + ! -- Objects + if (associated(this%isotherm)) then + deallocate (this%isotherm) + nullify (this%isotherm) + end if + ! ! -- deallocate parent call this%NumericalPackageType%da() end subroutine mst_da diff --git a/src/Model/TransportModel/Isotherm/FreundlichIsotherm.f90 b/src/Model/TransportModel/Isotherm/FreundlichIsotherm.f90 index 1d4bcb5a235..80d48d3f348 100644 --- a/src/Model/TransportModel/Isotherm/FreundlichIsotherm.f90 +++ b/src/Model/TransportModel/Isotherm/FreundlichIsotherm.f90 @@ -7,6 +7,10 @@ Module FreundlichIsothermModule Private Public :: FreundlichIsothermType + !> @brief Freundlich isotherm implementation of `IsothermType`. + !> + !> Sorbed concentration is cs = Kf*c^a. + !< type, extends(IsothermType) :: FreundlichIsothermType real(DP), pointer, dimension(:) :: Kf => null() !< Freundlich constant real(DP), pointer, dimension(:) :: a => null() !< Freundlich exponent @@ -21,6 +25,8 @@ Module FreundlichIsothermModule contains + !> @brief Constructor for Freundlich isotherm + !< function constructor(Kf, a) Result(isotherm) type(FreundlichIsothermType) :: isotherm ! -- dummy @@ -32,11 +38,15 @@ function constructor(Kf, a) Result(isotherm) end function constructor + !> @brief Evaluate the isotherm at a given node + !< function value(this, c, n) result(val) + ! -- return + real(DP) :: val !< isotherm value + ! -- dummy class(FreundlichIsothermType), intent(in) :: this - real(DP), dimension(:), intent(in) :: c - integer(I4B), intent(in) :: n - real(DP) :: val + real(DP), dimension(:), intent(in) :: c !< concentration array + integer(I4B), intent(in) :: n !< node index if (c(n) > 0.0_DP) then val = this%Kf(n) * c(n)**this%a(n) @@ -45,11 +55,15 @@ function value(this, c, n) result(val) end if end function value + !> @brief Evaluate derivative of the isotherm at a given node + !< function derivative(this, c, n) result(derv) + ! -- return + real(DP) :: derv !< derivative d(value)/dc evaluated at c + ! -- dummy class(FreundlichIsothermType), intent(in) :: this - real(DP), dimension(:), intent(in) :: c - integer(I4B), intent(in) :: n - real(DP) :: derv + real(DP), dimension(:), intent(in) :: c !< concentration array + integer(I4B), intent(in) :: n !< node index if (c(n) > 0.0_DP) then derv = this%a(n) * this%Kf(n) * c(n)**(this%a(n) - 1.0_DP) diff --git a/src/Model/TransportModel/Isotherm/IsothermFactory.f90 b/src/Model/TransportModel/Isotherm/IsothermFactory.f90 index 5406cf72895..b758537cb20 100644 --- a/src/Model/TransportModel/Isotherm/IsothermFactory.f90 +++ b/src/Model/TransportModel/Isotherm/IsothermFactory.f90 @@ -15,13 +15,19 @@ module IsothermFactoryModule contains + !> @brief Create an isotherm object based on type and parameters. + !> + !> Returns a pointer to a concrete `IsothermType` (or `null()` if sorption + !> is off). Uses `isotherm_type` to select Linear, Freundlich, or Langmuir, + !> passing `distcoef` and `sp2` as required by the chosen model. + !< function create_isotherm(isotherm_type, distcoef, sp2) result(isotherm) ! -- result - class(IsothermType), pointer :: isotherm + class(IsothermType), pointer :: isotherm !< allocated concrete isotherm or null() ! -- dummy - integer(I4B), intent(in) :: isotherm_type - real(DP), dimension(:), pointer, contiguous :: distcoef - real(DP), dimension(:), pointer, contiguous :: sp2 + integer(I4B), intent(in) :: isotherm_type !< enumerator from `IsothermEnumModule` + real(DP), dimension(:), pointer, contiguous :: distcoef !< primary coefficient (Kd, Kf, or Kl) + real(DP), dimension(:), pointer, contiguous :: sp2 !< secondary parameter (a for Freundlich, Sbar for Langmuir) select case (isotherm_type) case (SORPTION_OFF) @@ -33,7 +39,7 @@ function create_isotherm(isotherm_type, distcoef, sp2) result(isotherm) case (SORPTION_LANG) allocate (isotherm, source=LangmuirIsothermType(distcoef, sp2)) case default - call store_error('Sorption type not implemented yet.') + call store_error('Sorption type not implemented.') end select end function create_isotherm diff --git a/src/Model/TransportModel/Isotherm/IsothermInterface.f90 b/src/Model/TransportModel/Isotherm/IsothermInterface.f90 index 3a8bf9dd39e..2ab84ad3f56 100644 --- a/src/Model/TransportModel/Isotherm/IsothermInterface.f90 +++ b/src/Model/TransportModel/Isotherm/IsothermInterface.f90 @@ -5,28 +5,40 @@ module IsothermInterfaceModule private public :: IsothermType + !> @brief Interface for isotherms. + !< type, abstract :: IsothermType contains - procedure(value_interface), deferred :: value - procedure(derivative_interface), deferred :: derivative + procedure(value_if), deferred :: value !< Evaluate isotherm value at node n + procedure(derivative_if), deferred :: derivative !< Evaluate derivative d(value)/dc at node n end type IsothermType abstract interface - function value_interface(this, c, n) result(val) + !> @brief Evaluate the isotherm at a given node + !< + function value_if(this, c, n) result(val) + ! -- import import :: IsothermType, DP, I4B - class(IsothermType), intent(in) :: this - real(DP), dimension(:), intent(in) :: c - integer(I4B), intent(in) :: n - real(DP) :: val - end function value_interface + ! -- return + real(DP) :: val !< isotherm value + ! -- dummy + class(IsothermType), intent(in) :: this !< isotherm object + real(DP), dimension(:), intent(in) :: c !< concentration array + integer(I4B), intent(in) :: n !< node index + end function value_if - function derivative_interface(this, c, n) result(derv) + !> @brief Evaluate derivative of the isotherm at a given node + !< + function derivative_if(this, c, n) result(derv) + ! -- import import :: IsothermType, DP, I4B - class(IsothermType), intent(in) :: this - real(DP), dimension(:), intent(in) :: c - integer(I4B), intent(in) :: n - real(DP) :: derv - end function derivative_interface + ! -- return + real(DP) :: derv !< derivative d(value)/dc evaluated at c + ! -- dummy + class(IsothermType), intent(in) :: this !< isotherm object + real(DP), dimension(:), intent(in) :: c !< concentration array + integer(I4B), intent(in) :: n !< node index + end function derivative_if end interface end module IsothermInterfaceModule diff --git a/src/Model/TransportModel/Isotherm/LangmuirIsotherm.f90 b/src/Model/TransportModel/Isotherm/LangmuirIsotherm.f90 index e00ffdbef75..4b92c88a721 100644 --- a/src/Model/TransportModel/Isotherm/LangmuirIsotherm.f90 +++ b/src/Model/TransportModel/Isotherm/LangmuirIsotherm.f90 @@ -7,6 +7,10 @@ module LangmuirIsothermModule private public :: LangmuirIsothermType + !> @brief Langmuir isotherm implementation of `IsothermType`. + !! + !! Sorbed concentration is cs = (Sbar*Kl*c)/(1 + Kl*c). + !< type, extends(IsothermType) :: LangmuirIsothermType real(DP), pointer, dimension(:) :: Kl => null() !< Langmuir constant real(DP), pointer, dimension(:) :: Sbar => null() !< Total concentration of sorption sites @@ -20,22 +24,31 @@ module LangmuirIsothermModule end interface LangmuirIsothermType contains + !> @brief Constructor for Langmuir isotherm + !< function constructor(Kl, Sbar) Result(isotherm) + ! -- return type(LangmuirIsothermType) :: isotherm ! -- dummy - real(DP), pointer, dimension(:), intent(in) :: Kl - real(DP), pointer, dimension(:), intent(in) :: Sbar + real(DP), pointer, dimension(:), intent(in) :: Kl !< Langmuir constant + real(DP), pointer, dimension(:), intent(in) :: Sbar !< Total concentration of sorption sites ! -- local + isotherm%Kl => Kl isotherm%Sbar => Sbar end function constructor + !> @brief Evaluate the isotherm at a given node + !< function value(this, c, n) result(val) + ! -- return + real(DP) :: val !< isotherm value + ! -- dummy class(LangmuirIsothermType), intent(in) :: this - real(DP), dimension(:), intent(in) :: c - integer(I4B), intent(in) :: n - real(DP) :: val + real(DP), dimension(:), intent(in) :: c !< concentration array + integer(I4B), intent(in) :: n !< node index + ! -- local real(DP) :: denom if (c(n) > 0.0_DP) then @@ -46,11 +59,16 @@ function value(this, c, n) result(val) end if end function value + !> @brief Evaluate derivative of the isotherm at a given node + !< function derivative(this, c, n) result(derv) + ! -- return + real(DP) :: derv !< derivative d(value)/dc evaluated at c + ! -- dummy class(LangmuirIsothermType), intent(in) :: this - real(DP), dimension(:), intent(in) :: c - integer(I4B), intent(in) :: n - real(DP) :: derv + real(DP), dimension(:), intent(in) :: c !< concentration array + integer(I4B), intent(in) :: n !< node index + ! -- local real(DP) :: denom if (c(n) > 0.0_DP) then diff --git a/src/Model/TransportModel/Isotherm/LinearIsotherm.f90 b/src/Model/TransportModel/Isotherm/LinearIsotherm.f90 index 770adbd59e4..adc77301101 100644 --- a/src/Model/TransportModel/Isotherm/LinearIsotherm.f90 +++ b/src/Model/TransportModel/Isotherm/LinearIsotherm.f90 @@ -7,6 +7,10 @@ module LinearIsothermsModule private public :: LinearIsothermType + !> @brief Linear (Kd) isotherm implementation of `IsothermType`. + !> + !> Sorbed concentration is computed as cs = Kd*c. + !< type, extends(IsothermType) :: LinearIsothermType real(DP), pointer, dimension(:) :: Kd => null() !< distribution coefficient contains @@ -19,29 +23,42 @@ module LinearIsothermsModule end interface LinearIsothermType contains + !> @brief Constructor for Linear isotherm + !< function constructor(Kd) Result(isotherm) + ! -- return type(LinearIsothermType) :: isotherm ! -- dummy - real(DP), pointer, dimension(:), intent(in) :: Kd + real(DP), pointer, dimension(:), intent(in) :: Kd !< distribution coefficient ! -- local + isotherm%Kd => Kd end function constructor + !> @brief Evaluate the isotherm at a given node + !< function value(this, c, n) result(val) + ! -- return + real(DP) :: val !< isotherm value + ! -- dummy class(LinearIsothermType), intent(in) :: this - real(DP), dimension(:), intent(in) :: c - integer(I4B), intent(in) :: n - real(DP) :: val + real(DP), dimension(:), intent(in) :: c !< concentration array + integer(I4B), intent(in) :: n !< node index val = this%Kd(n) * c(n) end function value + !> @brief Evaluate derivative of the isotherm at a given node + !< function derivative(this, c, n) result(derv) + ! -- return + real(DP) :: derv !< derivative d(value)/dc evaluated at c + ! -- dummy class(LinearIsothermType), intent(in) :: this - real(DP), dimension(:), intent(in) :: c - integer(I4B), intent(in) :: n - real(DP) :: derv + real(DP), dimension(:), intent(in) :: c !< concentration array + integer(I4B), intent(in) :: n !< node index + ! -- local derv = this%Kd(n) end function derivative From 52fa2f57f10cbd350d32410f9c484822edb94ea6 Mon Sep 17 00:00:00 2001 From: Sunny Titus Date: Tue, 23 Dec 2025 14:51:49 +0100 Subject: [PATCH 06/12] Extend if statement to make sure the refactor doesnt change functional behavior --- src/Model/GroundWaterTransport/gwt-ist.f90 | 2 +- src/Model/GroundWaterTransport/gwt-mst.f90 | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Model/GroundWaterTransport/gwt-ist.f90 b/src/Model/GroundWaterTransport/gwt-ist.f90 index 5f3f3d93b8f..5d4a30cd3cf 100644 --- a/src/Model/GroundWaterTransport/gwt-ist.f90 +++ b/src/Model/GroundWaterTransport/gwt-ist.f90 @@ -556,7 +556,7 @@ subroutine ist_calc_csrb(this, cim) ! Calculate sorbed concentration do n = 1, size(cim) csrb = DZERO - if (this%ibound(n) > 0) then + if (this%ibound(n) > 0 .and. this%isrb /= SORPTION_OFF) then csrb = this%isotherm%value(cim, n) end if this%cimsrb(n) = csrb diff --git a/src/Model/GroundWaterTransport/gwt-mst.f90 b/src/Model/GroundWaterTransport/gwt-mst.f90 index 4e1650c0aeb..fa39da39001 100644 --- a/src/Model/GroundWaterTransport/gwt-mst.f90 +++ b/src/Model/GroundWaterTransport/gwt-mst.f90 @@ -798,7 +798,7 @@ subroutine mst_calc_csrb(this, cnew) ! Calculate sorbed concentration do n = 1, size(cnew) csrb = DZERO - if (this%ibound(n) > 0) then + if (this%ibound(n) > 0 .and. this%isrb /= SORPTION_OFF) then csrb = this%isotherm%value(cnew, n) end if this%csrb(n) = csrb From aa38b41dfcdebfb0fdbbc2df020d1d67c7ee64a7 Mon Sep 17 00:00:00 2001 From: Sunny Titus Date: Wed, 24 Dec 2025 15:17:48 +0100 Subject: [PATCH 07/12] Add unit test that captures mass inbalance with the current sorption implementation --- autotest/test_gwt_mst08.py | 222 +++++++++++++++++++++++++++++++++++++ 1 file changed, 222 insertions(+) create mode 100644 autotest/test_gwt_mst08.py 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() From 87ef0547bcefda5d7b7c75bca958843fc07a2812 Mon Sep 17 00:00:00 2001 From: Sunny Titus Date: Wed, 24 Dec 2025 15:32:43 +0100 Subject: [PATCH 08/12] Add new sorption formulation --- src/Model/GroundWaterTransport/gwt-mst.f90 | 67 +++++++++++----------- 1 file changed, 33 insertions(+), 34 deletions(-) diff --git a/src/Model/GroundWaterTransport/gwt-mst.f90 b/src/Model/GroundWaterTransport/gwt-mst.f90 index fa39da39001..dd66c083a46 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 contrubtion 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) From 6ed97c6b34611faf520ad2690824c6937c737555 Mon Sep 17 00:00:00 2001 From: Sunny Titus Date: Fri, 2 Jan 2026 14:07:27 +0100 Subject: [PATCH 09/12] Update unit tests --- autotest/test_gwt_mst05.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/autotest/test_gwt_mst05.py b/autotest/test_gwt_mst05.py index ca5fa12a9dc..e1a245c0dda 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.32355767, 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: From f4a3e56ddad67439636812594c8594dbc359473e Mon Sep 17 00:00:00 2001 From: Sunny Titus Date: Mon, 5 Jan 2026 12:37:55 +0100 Subject: [PATCH 10/12] Fix convergence issues witht he Freundlich Isotherm due to its infinite derivative at C=0 --- .../Isotherm/FreundlichIsotherm.f90 | 25 ++++++++++++++++--- 1 file changed, 21 insertions(+), 4 deletions(-) diff --git a/src/Model/TransportModel/Isotherm/FreundlichIsotherm.f90 b/src/Model/TransportModel/Isotherm/FreundlichIsotherm.f90 index 80d48d3f348..4208fd43ea5 100644 --- a/src/Model/TransportModel/Isotherm/FreundlichIsotherm.f90 +++ b/src/Model/TransportModel/Isotherm/FreundlichIsotherm.f90 @@ -2,6 +2,7 @@ Module FreundlichIsothermModule use KindModule, only: DP, I4B use IsothermInterfaceModule, only: IsothermType + use ConstantsModule, only: DZERO Implicit None Private @@ -10,6 +11,14 @@ Module FreundlichIsothermModule !> @brief Freundlich isotherm implementation of `IsothermType`. !> !> Sorbed concentration is cs = Kf*c^a. + !> + !> However, this expression has a singularity at c = 0 when a < 1, + !> leading to infinite derivative. To avoid this, the Freundlich isotherm + !> is modified as follows. + !> + !> Modified Sorbed concentration is cs = Kf*(c + eps)^a - Kf*eps^a + !> where eps = (K/(a*Kf))^(1/(a-1)) and K is a large constant (default 10). + !> This ensures that the derivative at c = 0 is below K. !< type, extends(IsothermType) :: FreundlichIsothermType real(DP), pointer, dimension(:) :: Kf => null() !< Freundlich constant @@ -47,9 +56,13 @@ function value(this, c, n) result(val) class(FreundlichIsothermType), intent(in) :: this real(DP), dimension(:), intent(in) :: c !< concentration array integer(I4B), intent(in) :: n !< node index + real(DP), parameter :: K = 10.0_dp !< constant to limit derivative at c=0 + real(DP) :: eps !< small concentration offset - if (c(n) > 0.0_DP) then - val = this%Kf(n) * c(n)**this%a(n) + eps = (K / (this%a(n) * this%Kf(n)))**(1.0_dp / (this%a(n) - 1.0_dp)) + + if (c(n) > DZERO) then + val = this%Kf(n) * (c(n) + eps)**this%a(n) - this%Kf(n) * eps**this%a(n) else val = 0.0_DP end if @@ -64,9 +77,13 @@ function derivative(this, c, n) result(derv) class(FreundlichIsothermType), intent(in) :: this real(DP), dimension(:), intent(in) :: c !< concentration array integer(I4B), intent(in) :: n !< node index + real(DP), parameter :: K = 10.0_dp !< constant to limit derivative at c=0 + real(DP) :: eps !< small concentration offset + + eps = (K / (this%a(n) * this%Kf(n)))**(1.0_dp / (this%a(n) - 1.0_dp)) - if (c(n) > 0.0_DP) then - derv = this%a(n) * this%Kf(n) * c(n)**(this%a(n) - 1.0_DP) + if (c(n) > DZERO) then + derv = this%a(n) * this%Kf(n) * ((c(n) + eps)**(this%a(n) - 1.0_DP)) else derv = 0.0_DP end if From 24241fb2da5701e161b1072329901fc261a8a987 Mon Sep 17 00:00:00 2001 From: Sunny Titus Date: Mon, 5 Jan 2026 16:56:23 +0100 Subject: [PATCH 11/12] Fix failing unit tests --- autotest/test_gwt_ist03.py | 11 ++--------- autotest/test_gwt_mst05.py | 2 +- src/Model/GroundWaterTransport/gwt-ist.f90 | 5 ++++- 3 files changed, 7 insertions(+), 11 deletions(-) diff --git a/autotest/test_gwt_ist03.py b/autotest/test_gwt_ist03.py index e97127b2411..7c67a4f7845 100644 --- a/autotest/test_gwt_ist03.py +++ b/autotest/test_gwt_ist03.py @@ -69,14 +69,10 @@ def build_models(idx, test): print_option="SUMMARY", outer_dvclose=hclose, outer_maximum=nouter, - under_relaxation="DBD", - under_relaxation_theta=0.7, inner_maximum=ninner, inner_dvclose=hclose, rcloserecord=rclose, linear_acceleration="BICGSTAB", - scaling_method="NONE", - reordering_method="NONE", relaxation_factor=relax, filename=f"{gwfname}.ims", ) @@ -132,13 +128,10 @@ def build_models(idx, test): print_option="SUMMARY", outer_dvclose=hclose, outer_maximum=nouter, - under_relaxation="NONE", inner_maximum=ninner, inner_dvclose=hclose, rcloserecord=rclose, linear_acceleration="BICGSTAB", - scaling_method="NONE", - reordering_method="NONE", relaxation_factor=relax, filename=f"{gwtname}.ims", ) @@ -339,7 +332,7 @@ def check_output(idx, test): csrb_answer = np.where( conc > 0, distcoef * sp2 * conc / (1 + distcoef * conc), 0 ) - if not np.allclose(csrb[:, 1:], csrb_answer[:, 1:]): + if not np.allclose(csrb[:, 1:], csrb_answer[:, 1:], atol=1e-6): diff = csrb - csrb_answer print("min and max difference") print(diff) @@ -355,7 +348,7 @@ def check_output(idx, test): cimsrb_answer = np.where( cim > 0, distcoef * sp2 * cim / (1 + distcoef * cim), 0 ) - if not np.allclose(cimsrb[:, 1:], cimsrb_answer[:, 1:]): + if not np.allclose(cimsrb[:, 1:], cimsrb_answer[:, 1:], atol=1e-6): diff = cimsrb - cimsrb_answer print("min and max difference") print(diff.min(), diff.max()) diff --git a/autotest/test_gwt_mst05.py b/autotest/test_gwt_mst05.py index e1a245c0dda..8294f1dcc50 100644 --- a/autotest/test_gwt_mst05.py +++ b/autotest/test_gwt_mst05.py @@ -276,7 +276,7 @@ def check_output(idx, test): assert False, f'could not load data from "{fpth}"' cnorm = obs["X008"] / 0.05 - cnorm_max = [0.32355767, 0.873678873] + cnorm_max = [0.324119806, 0.873678873] msg = f"{cnorm_max[idx]} /= {cnorm.max()}" assert np.allclose(cnorm_max[idx], cnorm.max(), atol=1e-6), msg diff --git a/src/Model/GroundWaterTransport/gwt-ist.f90 b/src/Model/GroundWaterTransport/gwt-ist.f90 index 5d4a30cd3cf..ebaa941fb04 100644 --- a/src/Model/GroundWaterTransport/gwt-ist.f90 +++ b/src/Model/GroundWaterTransport/gwt-ist.f90 @@ -1571,11 +1571,14 @@ function get_freundlich_kd(conc, kf, a) result(kd) real(DP), intent(in) :: conc !< solute concentration real(DP), intent(in) :: kf !< freundlich constant real(DP), intent(in) :: a !< freundlich exponent + real(DP), parameter :: K = 10.0_dp !< constant to limit derivative at c=0 + real(DP) :: eps !< small concentration offset ! -- return real(DP) :: kd !< effective distribution coefficient ! + eps = (K / (a * kf))**(1.0_dp / (a - 1.0_dp)) if (conc > DZERO) then - kd = kf * conc**(a - DONE) + kd = kf * (conc + eps)**(a - DONE) - kf * eps**(a - DONE) else kd = DZERO end if From 61b7da811bb568e62506333a46415b9e509ad287 Mon Sep 17 00:00:00 2001 From: Sunny Titus Date: Mon, 30 Mar 2026 11:40:00 +0200 Subject: [PATCH 12/12] Apply review comments. Split off Freundlich modification to a separate PR. Fix typo --- src/Model/GroundWaterTransport/gwt-ist.f90 | 5 +--- src/Model/GroundWaterTransport/gwt-mst.f90 | 2 +- .../Isotherm/FreundlichIsotherm.f90 | 25 +++---------------- 3 files changed, 6 insertions(+), 26 deletions(-) diff --git a/src/Model/GroundWaterTransport/gwt-ist.f90 b/src/Model/GroundWaterTransport/gwt-ist.f90 index ebaa941fb04..5d4a30cd3cf 100644 --- a/src/Model/GroundWaterTransport/gwt-ist.f90 +++ b/src/Model/GroundWaterTransport/gwt-ist.f90 @@ -1571,14 +1571,11 @@ function get_freundlich_kd(conc, kf, a) result(kd) real(DP), intent(in) :: conc !< solute concentration real(DP), intent(in) :: kf !< freundlich constant real(DP), intent(in) :: a !< freundlich exponent - real(DP), parameter :: K = 10.0_dp !< constant to limit derivative at c=0 - real(DP) :: eps !< small concentration offset ! -- return real(DP) :: kd !< effective distribution coefficient ! - eps = (K / (a * kf))**(1.0_dp / (a - 1.0_dp)) if (conc > DZERO) then - kd = kf * (conc + eps)**(a - DONE) - kf * eps**(a - DONE) + kd = kf * conc**(a - DONE) else kd = DZERO end if diff --git a/src/Model/GroundWaterTransport/gwt-mst.f90 b/src/Model/GroundWaterTransport/gwt-mst.f90 index dd66c083a46..651d4f1eaad 100644 --- a/src/Model/GroundWaterTransport/gwt-mst.f90 +++ b/src/Model/GroundWaterTransport/gwt-mst.f90 @@ -666,7 +666,7 @@ subroutine mst_cq_srb(this, nodes, cnew, cold, flowja) rate = rate + contribution ! -- Right-hand side contribution due to linearization - ! -- Note: this contrubtion should cancel with the matrix contribution when the outer loop is converged + ! -- 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 diff --git a/src/Model/TransportModel/Isotherm/FreundlichIsotherm.f90 b/src/Model/TransportModel/Isotherm/FreundlichIsotherm.f90 index 4208fd43ea5..80d48d3f348 100644 --- a/src/Model/TransportModel/Isotherm/FreundlichIsotherm.f90 +++ b/src/Model/TransportModel/Isotherm/FreundlichIsotherm.f90 @@ -2,7 +2,6 @@ Module FreundlichIsothermModule use KindModule, only: DP, I4B use IsothermInterfaceModule, only: IsothermType - use ConstantsModule, only: DZERO Implicit None Private @@ -11,14 +10,6 @@ Module FreundlichIsothermModule !> @brief Freundlich isotherm implementation of `IsothermType`. !> !> Sorbed concentration is cs = Kf*c^a. - !> - !> However, this expression has a singularity at c = 0 when a < 1, - !> leading to infinite derivative. To avoid this, the Freundlich isotherm - !> is modified as follows. - !> - !> Modified Sorbed concentration is cs = Kf*(c + eps)^a - Kf*eps^a - !> where eps = (K/(a*Kf))^(1/(a-1)) and K is a large constant (default 10). - !> This ensures that the derivative at c = 0 is below K. !< type, extends(IsothermType) :: FreundlichIsothermType real(DP), pointer, dimension(:) :: Kf => null() !< Freundlich constant @@ -56,13 +47,9 @@ function value(this, c, n) result(val) class(FreundlichIsothermType), intent(in) :: this real(DP), dimension(:), intent(in) :: c !< concentration array integer(I4B), intent(in) :: n !< node index - real(DP), parameter :: K = 10.0_dp !< constant to limit derivative at c=0 - real(DP) :: eps !< small concentration offset - eps = (K / (this%a(n) * this%Kf(n)))**(1.0_dp / (this%a(n) - 1.0_dp)) - - if (c(n) > DZERO) then - val = this%Kf(n) * (c(n) + eps)**this%a(n) - this%Kf(n) * eps**this%a(n) + if (c(n) > 0.0_DP) then + val = this%Kf(n) * c(n)**this%a(n) else val = 0.0_DP end if @@ -77,13 +64,9 @@ function derivative(this, c, n) result(derv) class(FreundlichIsothermType), intent(in) :: this real(DP), dimension(:), intent(in) :: c !< concentration array integer(I4B), intent(in) :: n !< node index - real(DP), parameter :: K = 10.0_dp !< constant to limit derivative at c=0 - real(DP) :: eps !< small concentration offset - - eps = (K / (this%a(n) * this%Kf(n)))**(1.0_dp / (this%a(n) - 1.0_dp)) - if (c(n) > DZERO) then - derv = this%a(n) * this%Kf(n) * ((c(n) + eps)**(this%a(n) - 1.0_DP)) + if (c(n) > 0.0_DP) then + derv = this%a(n) * this%Kf(n) * c(n)**(this%a(n) - 1.0_DP) else derv = 0.0_DP end if