From a645b9bb4211661dad81a0ce2e98521eb922a049 Mon Sep 17 00:00:00 2001 From: Sunny Titus Date: Mon, 23 Jun 2025 12:46:56 +0200 Subject: [PATCH 01/15] Refactor tsp-adv. Move different interpolation scheme to separate files --- make/makefile | 25 +- msvs/mf6core.vfproj | 6 + src/Model/Connection/GweGweConnection.f90 | 3 +- src/Model/Connection/GwtGwtConnection.f90 | 3 +- .../InterpolationScheme/AdvSchemeEnum.f90 | 11 + .../CentralDifferenceScheme.f90 | 68 +++++ .../IInterpolationScheme.f90 | 40 +++ .../InterpolationScheme/TVDScheme.f90 | 111 +++++++ .../InterpolationScheme/UpwindScheme.f90 | 60 ++++ src/Model/TransportModel/tsp-adv.f90 | 271 ++++-------------- src/meson.build | 5 + 11 files changed, 382 insertions(+), 221 deletions(-) create mode 100644 src/Model/TransportModel/InterpolationScheme/AdvSchemeEnum.f90 create mode 100644 src/Model/TransportModel/InterpolationScheme/CentralDifferenceScheme.f90 create mode 100644 src/Model/TransportModel/InterpolationScheme/IInterpolationScheme.f90 create mode 100644 src/Model/TransportModel/InterpolationScheme/TVDScheme.f90 create mode 100644 src/Model/TransportModel/InterpolationScheme/UpwindScheme.f90 diff --git a/make/makefile b/make/makefile index f6c71870589..d8d96b88d48 100644 --- a/make/makefile +++ b/make/makefile @@ -23,15 +23,15 @@ SOURCEDIR16=../src/Model/OverlandFlow SOURCEDIR17=../src/Model/ParticleTracking SOURCEDIR18=../src/Model/SurfaceWaterFlow SOURCEDIR19=../src/Model/TransportModel -SOURCEDIR20=../src/Solution -SOURCEDIR21=../src/Solution/LinearMethods -SOURCEDIR22=../src/Solution/PETSc -SOURCEDIR23=../src/Solution/ParticleTracker -SOURCEDIR24=../src/Timing -SOURCEDIR25=../src/Utilities -SOURCEDIR26=../src/Utilities/ArrayRead -SOURCEDIR27=../src/Utilities/Export -SOURCEDIR28=../src/Utilities/Export/tmp +SOURCEDIR20=../src/Model/TransportModel/InterpolationScheme +SOURCEDIR21=../src/Solution +SOURCEDIR22=../src/Solution/LinearMethods +SOURCEDIR23=../src/Solution/PETSc +SOURCEDIR24=../src/Solution/ParticleTracker +SOURCEDIR25=../src/Timing +SOURCEDIR26=../src/Utilities +SOURCEDIR27=../src/Utilities/ArrayRead +SOURCEDIR28=../src/Utilities/Export SOURCEDIR29=../src/Utilities/Idm SOURCEDIR30=../src/Utilities/Idm/mf6blockfile SOURCEDIR31=../src/Utilities/Idm/netcdf @@ -306,6 +306,8 @@ $(OBJDIR)/SimStages.o \ $(OBJDIR)/NumericalModel.o \ $(OBJDIR)/OutputControlData.o \ $(OBJDIR)/gwf-ic.o \ +$(OBJDIR)/tsp-fmi.o \ +$(OBJDIR)/IInterpolationScheme.o \ $(OBJDIR)/UzfEtUtil.o \ $(OBJDIR)/Xt3dAlgorithm.o \ $(OBJDIR)/TvBase.o \ @@ -328,11 +330,14 @@ $(OBJDIR)/SwfCxsUtils.o \ $(OBJDIR)/VirtualModel.o \ $(OBJDIR)/BaseExchange.o \ $(OBJDIR)/TspSpc.o \ -$(OBJDIR)/tsp-fmi.o \ $(OBJDIR)/GweInputData.o \ $(OBJDIR)/OutputControl.o \ $(OBJDIR)/tsp-ic.o \ +$(OBJDIR)/UpwindScheme.o \ +$(OBJDIR)/TVDScheme.o \ $(OBJDIR)/TspAdvOptions.o \ +$(OBJDIR)/CentralDifferenceScheme.o \ +$(OBJDIR)/AdvSchemeEnum.o \ $(OBJDIR)/UzfCellGroup.o \ $(OBJDIR)/Xt3dInterface.o \ $(OBJDIR)/gwf-tvk.o \ diff --git a/msvs/mf6core.vfproj b/msvs/mf6core.vfproj index 442825324d1..d6ed378136a 100644 --- a/msvs/mf6core.vfproj +++ b/msvs/mf6core.vfproj @@ -368,6 +368,12 @@ + + + + + + diff --git a/src/Model/Connection/GweGweConnection.f90 b/src/Model/Connection/GweGweConnection.f90 index b7cd571e282..d66aa5895ee 100644 --- a/src/Model/Connection/GweGweConnection.f90 +++ b/src/Model/Connection/GweGweConnection.f90 @@ -16,6 +16,7 @@ module GweGweConnectionModule use DistVariableModule use SimStagesModule use MatrixBaseModule + use AdvSchemeEnumModule implicit none private @@ -265,7 +266,7 @@ subroutine setGridExtent(this) if (hasAdv) then if (this%iIfaceAdvScheme == 2) then this%exg_stencil_depth = 2 - if (this%gweModel%adv%iadvwt == 2) then + if (this%gweModel%adv%iadvwt == ADV_SCHEME_TVD) 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 aa2a74037dd..9a4ee8fc2c9 100644 --- a/src/Model/Connection/GwtGwtConnection.f90 +++ b/src/Model/Connection/GwtGwtConnection.f90 @@ -16,6 +16,7 @@ module GwtGwtConnectionModule use DistVariableModule use SimStagesModule use MatrixBaseModule + use AdvSchemeEnumModule implicit none private @@ -258,7 +259,7 @@ subroutine setGridExtent(this) if (hasAdv) then if (this%iIfaceAdvScheme == 2) then this%exg_stencil_depth = 2 - if (this%gwtModel%adv%iadvwt == 2) then + if (this%gwtModel%adv%iadvwt == ADV_SCHEME_TVD) then this%int_stencil_depth = 2 end if end if diff --git a/src/Model/TransportModel/InterpolationScheme/AdvSchemeEnum.f90 b/src/Model/TransportModel/InterpolationScheme/AdvSchemeEnum.f90 new file mode 100644 index 00000000000..3e6ee17b1b3 --- /dev/null +++ b/src/Model/TransportModel/InterpolationScheme/AdvSchemeEnum.f90 @@ -0,0 +1,11 @@ +module AdvSchemeEnumModule + use KindModule, only: I4B + + implicit none + + ! Advection scheme codes + integer(I4B), parameter :: ADV_SCHEME_UPSTREAM = 0 + integer(I4B), parameter :: ADV_SCHEME_CENTRAL = 1 + integer(I4B), parameter :: ADV_SCHEME_TVD = 2 + +end module AdvSchemeEnumModule \ No newline at end of file diff --git a/src/Model/TransportModel/InterpolationScheme/CentralDifferenceScheme.f90 b/src/Model/TransportModel/InterpolationScheme/CentralDifferenceScheme.f90 new file mode 100644 index 00000000000..3a08d74743f --- /dev/null +++ b/src/Model/TransportModel/InterpolationScheme/CentralDifferenceScheme.f90 @@ -0,0 +1,68 @@ +module CentralDifferenceSchemeModule + use KindModule, only: DP, I4B + use ConstantsModule, only: DHALF, DONE + use IInterpolationSchemeModule, only: IInterpolationSchemeType, CoefficientsType + use BaseDisModule, only: DisBaseType + use TspFmiModule, only: TspFmiType + + implicit none + private + + public :: CentralDifferenceSchemeType + + type, extends(IInterpolationSchemeType) :: CentralDifferenceSchemeType + private + class(DisBaseType), pointer :: dis + type(TspFmiType), pointer :: fmi + contains + procedure :: compute + end type CentralDifferenceSchemeType + + interface CentralDifferenceSchemeType + module procedure constructor + end interface CentralDifferenceSchemeType + +contains + function constructor(dis, fmi) result(interpolation_scheme) + ! -- return + type(CentralDifferenceSchemeType) :: interpolation_scheme + ! --dummy + class(DisBaseType), pointer, intent(in) :: dis + type(TspFmiType), pointer, intent(in) :: fmi + + interpolation_scheme%dis => dis + interpolation_scheme%fmi => fmi + + end function constructor + + function compute(this, n, m, iposnm, phi) result(phi_face) + !-- return + type(CoefficientsType) :: phi_face + ! -- dummy + class(CentralDifferenceSchemeType), target :: this + integer(I4B), intent(in) :: n + integer(I4B), intent(in) :: m + integer(I4B), intent(in) :: iposnm + real(DP), intent(in), dimension(:) :: phi + ! -- local + real(DP) :: lnm, lmn, omega + + ! -- calculate weight based on distances between nodes and the shared + ! face of the connection + if (this%dis%con%ihc(this%dis%con%jas(iposnm)) == 0) then + ! -- vertical connection; assume cell is fully saturated + lnm = DHALF * (this%dis%top(n) - this%dis%bot(n)) + lmn = DHALF * (this%dis%top(m) - this%dis%bot(m)) + else + ! -- horizontal connection + lnm = this%dis%con%cl1(this%dis%con%jas(iposnm)) + lmn = this%dis%con%cl2(this%dis%con%jas(iposnm)) + end if + + omega = lmn / (lnm + lmn) + + phi_face%c_n = omega + phi_face%c_m = DONE - omega + + end function compute +end module CentralDifferenceSchemeModule diff --git a/src/Model/TransportModel/InterpolationScheme/IInterpolationScheme.f90 b/src/Model/TransportModel/InterpolationScheme/IInterpolationScheme.f90 new file mode 100644 index 00000000000..92dafbae5be --- /dev/null +++ b/src/Model/TransportModel/InterpolationScheme/IInterpolationScheme.f90 @@ -0,0 +1,40 @@ +module IInterpolationSchemeModule + use KindModule, only: DP, I4B + + implicit none + private + + public :: IInterpolationSchemeType + public :: CoefficientsType + + type :: CoefficientsType + real(DP) :: c_n = 0.0_dp + real(DP) :: c_m = 0.0_dp + real(DP) :: rhs = 0.0_dp + end type CoefficientsType + + type, abstract :: IInterpolationSchemeType + contains + procedure(compute_if), deferred :: compute + end type IInterpolationSchemeType + + abstract interface + + function compute_if(this, n, m, iposnm, phi) result(phi_face) + ! -- import + import DP, I4B + import IInterpolationSchemeType + import CoefficientsType + ! -- return + type(CoefficientsType) :: phi_face + ! -- dummy + class(IInterpolationSchemeType), target :: this + integer(I4B), intent(in) :: n + integer(I4B), intent(in) :: m + integer(I4B), intent(in) :: iposnm + real(DP), intent(in), dimension(:) :: phi + end function + + end interface + +end module IInterpolationSchemeModule diff --git a/src/Model/TransportModel/InterpolationScheme/TVDScheme.f90 b/src/Model/TransportModel/InterpolationScheme/TVDScheme.f90 new file mode 100644 index 00000000000..d0580534497 --- /dev/null +++ b/src/Model/TransportModel/InterpolationScheme/TVDScheme.f90 @@ -0,0 +1,111 @@ +module TVDSchemeModule + use KindModule, only: DP, I4B + use ConstantsModule, only: DONE, DZERO, DPREC, DHALF, DTWO + use IInterpolationSchemeModule, only: IInterpolationSchemeType, CoefficientsType + use BaseDisModule, only: DisBaseType + use TspFmiModule, only: TspFmiType + + implicit none + private + + public :: TVDSchemeType + + type, extends(IInterpolationSchemeType) :: TVDSchemeType + private + class(DisBaseType), pointer :: dis + type(TspFmiType), pointer :: fmi + integer(I4B), dimension(:), pointer, contiguous :: ibound => null() !< pointer to model ibound + contains + procedure :: compute + end type TVDSchemeType + + interface TVDSchemeType + module procedure constructor + end interface TVDSchemeType + +contains + function constructor(dis, fmi, ibound) result(interpolation_scheme) + ! -- return + type(TVDSchemeType) :: interpolation_scheme + ! --dummy + class(DisBaseType), pointer, intent(in) :: dis + type(TspFmiType), pointer, intent(in) :: fmi + integer(I4B), dimension(:), pointer, contiguous, intent(in) :: ibound + + interpolation_scheme%dis => dis + interpolation_scheme%fmi => fmi + interpolation_scheme%ibound => ibound + + end function constructor + + function compute(this, n, m, iposnm, phi) result(phi_face) + !-- return + type(CoefficientsType), target :: phi_face + ! -- dummy + class(TVDSchemeType), target :: this + integer(I4B), intent(in) :: n + integer(I4B), intent(in) :: m + integer(I4B), intent(in) :: iposnm + real(DP), intent(in), dimension(:) :: phi + ! -- local + integer(I4B) :: ipos, isympos, iup, idn, i2up, j + real(DP) :: qnm, qmax, qupj, elupdn, elup2up + real(DP) :: smooth, cdiff, alimiter + real(DP), pointer :: coef_up, coef_dn + ! + ! -- Find upstream node + isympos = this%dis%con%jas(iposnm) + qnm = this%fmi%gwfflowja(iposnm) + if (qnm > DZERO) then + ! -- positive flow into n means m is upstream + iup = m + idn = n + + coef_up => phi_face%c_m + coef_dn => phi_face%c_n + else + iup = n + idn = m + + coef_up => phi_face%c_n + coef_dn => phi_face%c_m + end if + elupdn = this%dis%con%cl1(isympos) + this%dis%con%cl2(isympos) + ! + ! -- Add low order terms + coef_up = DONE + ! + ! -- Add high order terms + ! + ! -- Find second node upstream to iup + i2up = 0 + qmax = DZERO + do ipos = this%dis%con%ia(iup) + 1, this%dis%con%ia(iup + 1) - 1 + j = this%dis%con%ja(ipos) + if (this%ibound(j) == 0) cycle + qupj = this%fmi%gwfflowja(ipos) + isympos = this%dis%con%jas(ipos) + if (qupj > qmax) then + qmax = qupj + i2up = j + elup2up = this%dis%con%cl1(isympos) + this%dis%con%cl2(isympos) + end if + end do + ! + ! -- Calculate flux limiting term + if (i2up > 0) then + smooth = DZERO + cdiff = ABS(phi(idn) - phi(iup)) + if (cdiff > DPREC) then + smooth = (phi(iup) - phi(i2up)) / elup2up * & + elupdn / (phi(idn) - phi(iup)) + end if + if (smooth > DZERO) then + alimiter = DTWO * smooth / (DONE + smooth) + phi_face%rhs = -DHALF * alimiter * (phi(idn) - phi(iup)) + end if + end if + + end function compute + +end module TVDSchemeModule diff --git a/src/Model/TransportModel/InterpolationScheme/UpwindScheme.f90 b/src/Model/TransportModel/InterpolationScheme/UpwindScheme.f90 new file mode 100644 index 00000000000..8e14b7caedb --- /dev/null +++ b/src/Model/TransportModel/InterpolationScheme/UpwindScheme.f90 @@ -0,0 +1,60 @@ +module UpwindSchemeModule + use KindModule, only: DP, I4B + use ConstantsModule, only: DZERO + use IInterpolationSchemeModule, only: IInterpolationSchemeType, CoefficientsType + use BaseDisModule, only: DisBaseType + use TspFmiModule, only: TspFmiType + + implicit none + private + + public :: UpwindSchemeType + + type, extends(IInterpolationSchemeType) :: UpwindSchemeType + private + class(DisBaseType), pointer :: dis + type(TspFmiType), pointer :: fmi + contains + procedure :: compute + end type UpwindSchemeType + + interface UpwindSchemeType + module procedure constructor + end interface UpwindSchemeType + +contains + function constructor(dis, fmi) result(interpolation_scheme) + ! -- return + type(UpwindSchemeType) :: interpolation_scheme + ! --dummy + class(DisBaseType), pointer, intent(in) :: dis + type(TspFmiType), pointer, intent(in) :: fmi + + interpolation_scheme%dis => dis + interpolation_scheme%fmi => fmi + + end function constructor + + function compute(this, n, m, iposnm, phi) result(phi_face) + !-- return + type(CoefficientsType) :: phi_face + ! -- dummy + class(UpwindSchemeType), target :: this + integer(I4B), intent(in) :: n + integer(I4B), intent(in) :: m + integer(I4B), intent(in) :: iposnm + real(DP), intent(in), dimension(:) :: phi + ! -- local + real(DP) :: qnm + + ! -- Compute the coefficients for the upwind scheme + qnm = this%fmi%gwfflowja(iposnm) + + if (qnm < DZERO) then + phi_face%c_n = 1.0_dp + else + phi_face%c_m = 1.0_dp + end if + + end function compute +end module UpwindSchemeModule diff --git a/src/Model/TransportModel/tsp-adv.f90 b/src/Model/TransportModel/tsp-adv.f90 index 8eb856e79fd..09412fe657a 100644 --- a/src/Model/TransportModel/tsp-adv.f90 +++ b/src/Model/TransportModel/tsp-adv.f90 @@ -1,13 +1,18 @@ module TspAdvModule use KindModule, only: DP, I4B - use ConstantsModule, only: DONE, DZERO, DHALF, DTWO, DNODATA, DPREC, & - LINELENGTH + use ConstantsModule, only: DONE, DZERO, DNODATA, DPREC, LINELENGTH use NumericalPackageModule, only: NumericalPackageType use BaseDisModule, only: DisBaseType use TspFmiModule, only: TspFmiType use TspAdvOptionsModule, only: TspAdvOptionsType - use MatrixBaseModule + use MatrixBaseModule, only: MatrixBaseType + ! -- Interpolation schemes + use IInterpolationSchemeModule, only: IInterpolationSchemeType, CoefficientsType + use AdvSchemeEnumModule + use UpwindSchemeModule, only: UpwindSchemeType + use CentralDifferenceSchemeModule, only: CentralDifferenceSchemeType + use TVDSchemeModule, only: TVDSchemeType implicit none private @@ -15,13 +20,13 @@ module TspAdvModule public :: adv_cr type, extends(NumericalPackageType) :: TspAdvType - - integer(I4B), pointer :: iadvwt => null() !< advection scheme (0 up, 1 central, 2 tvd) + integer(I4B), pointer :: iadvwt => null() !< advection scheme. See ADV_SCHEME_* constants real(DP), pointer :: ats_percel => null() !< user-specified fractional number of cells advection can move a particle during one time step integer(I4B), dimension(:), pointer, contiguous :: ibound => null() !< pointer to model ibound type(TspFmiType), pointer :: fmi => null() !< pointer to fmi object real(DP), pointer :: eqnsclfac => null() !< governing equation scale factor; =1. for solute; =rhow*cpw for energy + class(IInterpolationSchemeType), allocatable :: face_interpolation !< interpolation scheme for face values contains procedure :: adv_df @@ -33,10 +38,6 @@ module TspAdvModule procedure :: allocate_scalars procedure, private :: read_options - procedure, private :: advqtvd - procedure, private :: advtvd_bd - procedure :: adv_weight - procedure :: advtvd end type TspAdvType @@ -109,16 +110,33 @@ end subroutine adv_df !< subroutine adv_ar(this, dis, ibound) ! -- modules + use SimModule, only: store_error ! -- dummy class(TspAdvType) :: this class(DisBaseType), pointer, intent(in) :: dis integer(I4B), dimension(:), pointer, contiguous, intent(in) :: ibound ! -- local - ! -- formats + integer(I4B) :: iadvwt_value ! ! -- adv pointers to arguments that were passed in this%dis => dis this%ibound => ibound + ! + ! -- Create interpolation scheme + iadvwt_value = this%iadvwt ! Dereference iadvwt to work with case statement + select case (iadvwt_value) + case (ADV_SCHEME_UPSTREAM) + this%face_interpolation = & + UpwindSchemeType(this%dis, this%fmi) + case (ADV_SCHEME_CENTRAL) + this%face_interpolation = & + CentralDifferenceSchemeType(this%dis, this%fmi) + case (ADV_SCHEME_TVD) + this%face_interpolation = & + TVDSchemeType(this%dis, this%fmi, this%ibound) + case default + call store_error("Unknown advection scheme", terminate=.TRUE.) + end select end subroutine adv_ar !> @brief Calculate maximum time step length @@ -193,137 +211,36 @@ end subroutine adv_dt subroutine adv_fc(this, nodes, matrix_sln, idxglo, cnew, rhs) ! -- modules ! -- dummy - class(TspAdvType) :: this - integer, intent(in) :: nodes - class(MatrixBaseType), pointer :: matrix_sln - integer(I4B), intent(in), dimension(:) :: idxglo - real(DP), intent(in), dimension(:) :: cnew - real(DP), dimension(:), intent(inout) :: rhs + class(TspAdvType) :: this !< this instance + integer(I4B), intent(in) :: nodes !< number of nodes + class(MatrixBaseType), pointer :: matrix_sln !< pointer to solution matrix + integer(I4B), intent(in), dimension(:) :: idxglo !< global indices for matrix + real(DP), intent(in), dimension(:) :: cnew !< new concentration/temperature values + real(DP), dimension(:), intent(inout) :: rhs !< right-hand side vector ! -- local integer(I4B) :: n, m, idiag, ipos - real(DP) :: omega, qnm - ! - ! -- Calculate advection terms and add to solution rhs and hcof. qnm - ! is the volumetric flow rate and has dimensions of L^/T. + real(DP) :: qnm !< volumetric flow rate + type(CoefficientsType) :: coefficients + + ! Calculate internal domain fluxes and add to matrix_sln and rhs. do n = 1, nodes - if (this%ibound(n) == 0) cycle + if (this%ibound(n) == 0) cycle ! skip inactive nodes idiag = this%dis%con%ia(n) do ipos = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1 - if (this%dis%con%mask(ipos) == 0) cycle m = this%dis%con%ja(ipos) - if (this%ibound(m) == 0) cycle + if (this%ibound(m) == 0) cycle ! skip inactive nodes + if (this%dis%con%mask(ipos) == 0) cycle ! skip masked connections + qnm = this%fmi%gwfflowja(ipos) * this%eqnsclfac - omega = this%adv_weight(this%iadvwt, ipos, n, m, qnm) - call matrix_sln%add_value_pos(idxglo(ipos), qnm * (DONE - omega)) - call matrix_sln%add_value_pos(idxglo(idiag), qnm * omega) + coefficients = this%face_interpolation%compute(n, m, ipos, cnew) + + call matrix_sln%add_value_pos(idxglo(idiag), qnm * coefficients%c_n) + call matrix_sln%add_value_pos(idxglo(ipos), qnm * coefficients%c_m) + rhs(n) = rhs(n) + qnm * coefficients%rhs end do end do - ! - ! -- TVD - if (this%iadvwt == 2) then - do n = 1, nodes - if (this%ibound(n) == 0) cycle - call this%advtvd(n, cnew, rhs) - end do - end if end subroutine adv_fc - !> @brief Calculate TVD - !! - !! Use explicit scheme to calculate the advective component of transport. - !! TVD is an acronym for Total-Variation Diminishing - !< - subroutine advtvd(this, n, cnew, rhs) - ! -- modules - ! -- dummy - class(TspAdvType) :: this - integer(I4B), intent(in) :: n - real(DP), dimension(:), intent(in) :: cnew - real(DP), dimension(:), intent(inout) :: rhs - ! -- local - real(DP) :: qtvd - integer(I4B) :: m, ipos - ! - ! -- Loop through each n connection. This will - do ipos = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1 - if (this%dis%con%mask(ipos) == 0) cycle - m = this%dis%con%ja(ipos) - if (m > n .and. this%ibound(m) /= 0) then - qtvd = this%advqtvd(n, m, ipos, cnew) - rhs(n) = rhs(n) - qtvd - rhs(m) = rhs(m) + qtvd - end if - end do - end subroutine advtvd - - !> @brief Calculate TVD - !! - !! Use explicit scheme to calculate the advective component of transport. - !! TVD is an acronym for Total-Variation Diminishing - !< - function advqtvd(this, n, m, iposnm, cnew) result(qtvd) - ! -- modules - use ConstantsModule, only: DPREC - ! -- return - real(DP) :: qtvd - ! -- dummy - class(TspAdvType) :: this - integer(I4B), intent(in) :: n - integer(I4B), intent(in) :: m - integer(I4B), intent(in) :: iposnm - real(DP), dimension(:), intent(in) :: cnew - ! -- local - integer(I4B) :: ipos, isympos, iup, idn, i2up, j - real(DP) :: qnm, qmax, qupj, elupdn, elup2up - real(DP) :: smooth, cdiff, alimiter - ! - ! -- initialize - qtvd = DZERO - ! - ! -- Find upstream node - isympos = this%dis%con%jas(iposnm) - qnm = this%fmi%gwfflowja(iposnm) - if (qnm > DZERO) then - ! -- positive flow into n means m is upstream - iup = m - idn = n - else - iup = n - idn = m - end if - elupdn = this%dis%con%cl1(isympos) + this%dis%con%cl2(isympos) - ! - ! -- Find second node upstream to iup - i2up = 0 - qmax = DZERO - do ipos = this%dis%con%ia(iup) + 1, this%dis%con%ia(iup + 1) - 1 - j = this%dis%con%ja(ipos) - if (this%ibound(j) == 0) cycle - qupj = this%fmi%gwfflowja(ipos) - isympos = this%dis%con%jas(ipos) - if (qupj > qmax) then - qmax = qupj - i2up = j - elup2up = this%dis%con%cl1(isympos) + this%dis%con%cl2(isympos) - end if - end do - ! - ! -- Calculate flux limiting term - if (i2up > 0) then - smooth = DZERO - cdiff = ABS(cnew(idn) - cnew(iup)) - if (cdiff > DPREC) then - smooth = (cnew(iup) - cnew(i2up)) / elup2up * & - elupdn / (cnew(idn) - cnew(iup)) - end if - if (smooth > DZERO) then - alimiter = DTWO * smooth / (DONE + smooth) - qtvd = DHALF * alimiter * qnm * (cnew(idn) - cnew(iup)) - qtvd = qtvd * this%eqnsclfac - end if - end if - end function advqtvd - !> @brief Calculate advection contribution to flowja !< subroutine adv_cq(this, cnew, flowja) @@ -334,53 +251,30 @@ subroutine adv_cq(this, cnew, flowja) real(DP), intent(inout), dimension(:) :: flowja ! -- local integer(I4B) :: nodes - integer(I4B) :: n, m, idiag, ipos - real(DP) :: omega, qnm + integer(I4B) :: n, m, ipos + real(DP) :: qnm + type(CoefficientsType) :: coefficients ! ! -- Calculate advection and add to flowja. qnm is the volumetric flow ! rate and has dimensions of L^/T. nodes = this%dis%nodes + do n = 1, nodes if (this%ibound(n) == 0) cycle - idiag = this%dis%con%ia(n) do ipos = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1 m = this%dis%con%ja(ipos) if (this%ibound(m) == 0) cycle qnm = this%fmi%gwfflowja(ipos) * this%eqnsclfac - omega = this%adv_weight(this%iadvwt, ipos, n, m, qnm) - flowja(ipos) = flowja(ipos) + qnm * omega * cnew(n) + & - qnm * (DONE - omega) * cnew(m) - end do - end do - ! - ! -- TVD - if (this%iadvwt == 2) call this%advtvd_bd(cnew, flowja) - end subroutine adv_cq - !> @brief Add TVD contribution to flowja - subroutine advtvd_bd(this, cnew, flowja) - ! -- modules - ! -- dummy - class(TspAdvType) :: this - real(DP), dimension(:), intent(in) :: cnew - real(DP), dimension(:), intent(inout) :: flowja - ! -- local - real(DP) :: qtvd, qnm - integer(I4B) :: nodes, n, m, ipos - ! - nodes = this%dis%nodes - do n = 1, nodes - if (this%ibound(n) == 0) cycle - do ipos = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1 - m = this%dis%con%ja(ipos) - if (this%ibound(m) /= 0) then - qnm = this%fmi%gwfflowja(ipos) - qtvd = this%advqtvd(n, m, ipos, cnew) - flowja(ipos) = flowja(ipos) + qtvd - end if + coefficients = this%face_interpolation%compute(n, m, ipos, cnew) + flowja(ipos) = flowja(ipos) & + + qnm * coefficients%c_n * cnew(n) & + + qnm * coefficients%c_m * cnew(m) & + - qnm * coefficients%rhs end do end do - end subroutine advtvd_bd + + end subroutine adv_cq !> @brief Deallocate memory !< @@ -423,7 +317,7 @@ subroutine allocate_scalars(this) call mem_allocate(this%ats_percel, 'ATS_PERCEL', this%memoryPath) ! ! -- Initialize - this%iadvwt = 0 + this%iadvwt = ADV_SCHEME_UPSTREAM this%ats_percel = DNODATA ! ! -- Advection creates an asymmetric coefficient matrix @@ -464,13 +358,13 @@ subroutine read_options(this) call this%parser%GetStringCaps(keyword) select case (keyword) case ('UPSTREAM') - this%iadvwt = 0 + this%iadvwt = ADV_SCHEME_UPSTREAM write (this%iout, fmtiadvwt) 'UPSTREAM' case ('CENTRAL') - this%iadvwt = 1 + this%iadvwt = ADV_SCHEME_CENTRAL write (this%iout, fmtiadvwt) 'CENTRAL' case ('TVD') - this%iadvwt = 2 + this%iadvwt = ADV_SCHEME_TVD write (this%iout, fmtiadvwt) 'TVD' case default write (errmsg, '(a, a)') & @@ -497,45 +391,4 @@ subroutine read_options(this) end if end subroutine read_options - !> @ brief Advection weight - !! - !! Calculate the advection weight - !< - function adv_weight(this, iadvwt, ipos, n, m, qnm) result(omega) - ! -- return - real(DP) :: omega - ! -- dummy - class(TspAdvType) :: this - integer, intent(in) :: iadvwt - integer, intent(in) :: ipos - integer, intent(in) :: n - integer, intent(in) :: m - real(DP), intent(in) :: qnm - ! -- local - real(DP) :: lnm, lmn - - select case (iadvwt) - case (1) - ! -- calculate weight based on distances between nodes and the shared - ! face of the connection - if (this%dis%con%ihc(this%dis%con%jas(ipos)) == 0) then - ! -- vertical connection; assume cell is fully saturated - lnm = DHALF * (this%dis%top(n) - this%dis%bot(n)) - lmn = DHALF * (this%dis%top(m) - this%dis%bot(m)) - else - ! -- horizontal connection - lnm = this%dis%con%cl1(this%dis%con%jas(ipos)) - lmn = this%dis%con%cl2(this%dis%con%jas(ipos)) - end if - omega = lmn / (lnm + lmn) - case (0, 2) - ! -- use upstream weighting for upstream and tvd schemes - if (qnm > DZERO) then - omega = DZERO - else - omega = DONE - end if - end select - end function adv_weight - end module TspAdvModule diff --git a/src/meson.build b/src/meson.build index 55d80a403a6..b25f151c9fc 100644 --- a/src/meson.build +++ b/src/meson.build @@ -267,6 +267,11 @@ modflow_sources = files( 'Model' / 'ModelUtilities' / 'VectorInterpolation.f90', 'Model' / 'ModelUtilities' / 'Xt3dAlgorithm.f90', 'Model' / 'ModelUtilities' / 'Xt3dInterface.f90', + 'Model' / 'TransportModel' / 'InterpolationScheme' / 'AdvSchemeEnum.f90', + 'Model' / 'TransportModel' / 'InterpolationScheme' / 'CentralDifferenceScheme.f90', + 'Model' / 'TransportModel' / 'InterpolationScheme' / 'IInterpolationScheme.f90', + 'Model' / 'TransportModel' / 'InterpolationScheme' / 'TVDScheme.f90', + 'Model' / 'TransportModel' / 'InterpolationScheme' / 'UpwindScheme.f90', 'Model' / 'TransportModel' / 'tsp.f90', 'Model' / 'TransportModel' / 'tsp-adv.f90', 'Model' / 'TransportModel' / 'tsp-apt.f90', From e4f889885831a6663a12e0039aa3838d3e55a83c Mon Sep 17 00:00:00 2001 From: Sunny Titus Date: Fri, 27 Jun 2025 15:29:20 +0200 Subject: [PATCH 02/15] =?UTF-8?q?=EF=BB=BFAdd=20the=20SVD-algorithm=20and?= =?UTF-8?q?=20the=20pseudoinverse=20method?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- autotest/TestPseudoInverse.f90 | 182 +++++++++++ autotest/TestSVD.f90 | 460 ++++++++++++++++++++++++++ autotest/meson.build | 3 + autotest/tester.f90 | 4 + autotest/tester_utils.f90 | 184 +++++++++++ make/makefile | 62 ++-- msvs/mf6core.vfproj | 3 + src/Utilities/LinearAlgebraUtils.f90 | 51 +++ src/Utilities/PseudoInverse.f90 | 55 ++++ src/Utilities/SVD.f90 | 469 +++++++++++++++++++++++++++ src/meson.build | 3 + 11 files changed, 1443 insertions(+), 33 deletions(-) create mode 100644 autotest/TestPseudoInverse.f90 create mode 100644 autotest/TestSVD.f90 create mode 100644 autotest/tester_utils.f90 create mode 100644 src/Utilities/LinearAlgebraUtils.f90 create mode 100644 src/Utilities/PseudoInverse.f90 create mode 100644 src/Utilities/SVD.f90 diff --git a/autotest/TestPseudoInverse.f90 b/autotest/TestPseudoInverse.f90 new file mode 100644 index 00000000000..c1be5d62ffd --- /dev/null +++ b/autotest/TestPseudoInverse.f90 @@ -0,0 +1,182 @@ +module TestPseudoInverse + use KindModule, only: I4B, DP + use ConstantsModule, only: DSAME + use testdrive, only: error_type, unittest_type, new_unittest, check, test_failed + use PseudoInverseModule, only: pinv + use LinearAlgebraUtilsModule, only: eye + use TesterUtils, only: check_same_matrix + + implicit none + private + public :: collect_pinv + +contains + + subroutine collect_pinv(testsuite) + type(unittest_type), allocatable, intent(out) :: testsuite(:) + + testsuite = [ & + new_unittest("identity", test_identity), & + new_unittest("square", test_square), & + new_unittest("tall", test_tall), & + new_unittest("wide", test_wide), & + new_unittest("rank_def", test_rank_def), & + new_unittest("last_row_col_zero", test_last_row_col_zero), & + new_unittest("zero", test_zero) & + ] + end subroutine collect_pinv + + !> @brief Test the pseudoinverse of the identity matrix. + !! + !! This test verifies that the pseudoinverse of a 3x3 identity matrix is itself. + !< + subroutine test_identity(error) + type(error_type), allocatable, intent(out) :: error + real(DP), dimension(3, 3) :: A, B_expected, B_actual + + A = eye(3) ! Create a 3x3 identity matrix + B_expected = A + B_actual = pinv(A) + + call check_same_matrix(error, B_actual, B_expected, tolerance=DSAME) + if (allocated(error)) return + end subroutine + + !> @brief Test the pseudoinverse of a square invertible matrix. + !! + !! This test verifies that the pseudoinverse of a regular 2x2 matrix matches its + !! analytical inverse. + !< + subroutine test_square(error) + type(error_type), allocatable, intent(out) :: error + real(DP), dimension(2, 2) :: A, B_expected, B_actual + + A = reshape([ & + 4.0_DP, 7.0_DP, & + 2.0_DP, 6.0_DP & + ], [2, 2]) + B_expected = reshape([ & + 0.6_DP, -0.7_DP, & + -0.2_DP, 0.4_DP & + ], [2, 2]) + B_actual = pinv(A) + + call check_same_matrix(error, B_actual, B_expected, tolerance=DSAME) + if (allocated(error)) return + end subroutine + + !> @brief Test the pseudoinverse of a tall matrix (more rows than columns). + !! + !! This test checks that the pseudoinverse of a 3x2 matrix matches the expected result, + !! as computed analytically. + !< + subroutine test_tall(error) + type(error_type), allocatable, intent(out) :: error + real(DP), dimension(3, 2) :: A + real(DP), dimension(2, 3) :: B_expected, B_actual + + A = reshape([ & + 1.0_DP, 2.0_DP, 3.0_DP, & + 4.0_DP, 5.0_DP, 6.0_DP & + ], [3, 2]) + + B_expected = reshape([ & + -17.0_DP / 18.0_DP, 4.0_DP / 9.0_DP, & + -1.0_DP / 9.0_DP, 1.0_DP / 9.0_DP, & + 13.0_DP / 18.0_DP, -2.0_DP / 9.0_DP & + ], [2, 3]) + + B_actual = pinv(A) + + call check_same_matrix(error, B_actual, B_expected, tolerance=DSAME) + if (allocated(error)) return + end subroutine + + !> @brief Test the pseudoinverse of a wide matrix (more columns than rows). + !! + !! This test checks that the pseudoinverse of a 2x3 matrix matches the expected result, + !! as computed analytically. + !< + subroutine test_wide(error) + type(error_type), allocatable, intent(out) :: error + real(DP), dimension(2, 3) :: A + real(DP), dimension(3, 2) :: B_expected, B_actual + + A = reshape([ & + 1.0_DP, 2.0_DP, & + 3.0_DP, 4.0_DP, & + 5.0_DP, 6.0_DP & + ], [2, 3]) + B_expected = reshape([ & + -4.0_DP / 3.0_DP, -1.0_DP / 3.0_DP, 2.0_DP / 3.0_DP, & + 13.0_DP / 12.0_DP, 1.0_DP / 3.0_DP, -5.0_DP / 12.0_DP & + ], [3, 2]) + B_actual = pinv(A) + + call check_same_matrix(error, B_actual, B_expected, tolerance=DSAME) + if (allocated(error)) return + end subroutine + + !> @brief Test the pseudoinverse of a rank-deficient matrix. + !! + !! This test verifies that the pseudoinverse of a singular (rank-deficient) 2x2 matrix + !! matches the expected result, ensuring correct handling of singular values. + !! The test matrix is rank-deficient because its second row is a multiple of the first + !! (row 2 = 2 × row 1), so its rank is 1 instead of 2. This checks that the pseudoinverse + !! implementation correctly handles matrices that are not full rank. + !< + subroutine test_rank_def(error) + type(error_type), allocatable, intent(out) :: error + real(DP), dimension(2, 2) :: A, B_expected, B_actual + + A = reshape([ & + 1.0_DP, 2.0_DP, & + 2.0_DP, 4.0_DP & + ], [2, 2]) + B_expected = reshape([ & + 0.04_DP, 0.08_DP, & + 0.08_DP, 0.16_DP], [2, 2]) + B_actual = pinv(A) + + call check_same_matrix(error, B_actual, B_expected, tolerance=DSAME) + if (allocated(error)) return + end subroutine + + !> @brief Test the pseudoinverse of a 3x3 matrix with the last row and column as zeros. + !! + !! This test a common case encountered with the TVD limiter in which the last row and column + !! of the matrix are zero due to solving a 2D or 1D problem in a 3D context. + !< + subroutine test_last_row_col_zero(error) + type(error_type), allocatable, intent(out) :: error + real(DP), dimension(3, 3) :: A, B_expected, B_actual + + A = 0.0_DP + A(1:2, 1:2) = reshape([4.0_DP, 7.0_DP, 2.0_DP, 6.0_DP], [2, 2]) + + B_expected = 0.0_DP + B_expected(1:2, 1:2) = reshape([0.6_DP, -0.7_DP, -0.2_DP, 0.4_DP], [2, 2]) + + B_actual = pinv(A) + + call check_same_matrix(error, B_actual, B_expected, tolerance=DSAME) + if (allocated(error)) return + end subroutine test_last_row_col_zero + + !> @brief Test the pseudoinverse of a zero matrix. + !! + !! This test checks that the pseudoinverse of a 3x3 zero matrix is also a zero matrix. + !< + subroutine test_zero(error) + type(error_type), allocatable, intent(out) :: error + real(DP), dimension(3, 3) :: A, B_expected, B_actual + + A = 0.0_DP + B_expected = 0.0_DP + B_actual = pinv(A) + + call check_same_matrix(error, B_actual, B_expected, tolerance=DSAME) + if (allocated(error)) return + end subroutine + +end module TestPseudoInverse diff --git a/autotest/TestSVD.f90 b/autotest/TestSVD.f90 new file mode 100644 index 00000000000..acecf92913e --- /dev/null +++ b/autotest/TestSVD.f90 @@ -0,0 +1,460 @@ +module TestSVD + use KindModule, only: I4B, DP + use testdrive, only: error_type, unittest_type, new_unittest, check + use SVDModule, only: SVD, bidiagonal_decomposition, bidiagonal_qr_decomposition + use LinearAlgebraUtilsModule, only: eye + use TesterUtils, only: check_same_matrix, & + check_matrix_row_orthogonal, & + check_matrix_column_orthogonal, & + check_matrix_is_diagonal, & + check_matrix_bidiagonal + + implicit none + private + public :: collect_svd + +contains + + subroutine collect_svd(testsuite) + type(unittest_type), allocatable, intent(out) :: testsuite(:) + + testsuite = [ & + new_unittest("bidiagonal_decomposition", & + test_bidiagonal_decomposition), & + new_unittest("diagonalize_matrix", & + test_diagonalize_matrix), & + new_unittest("svd_input_bidiagonal", & + test_svd_input_bidiagonal), & + new_unittest("svd_input_square_matrix", & + test_svd_input_square_matrix), & + new_unittest("svd_input_tall_matrix", & + test_svd_input_tall_matrix), & + new_unittest("svd_input_wide_matrix", & + test_svd_input_wide_matrix), & + new_unittest("svd_zero_on_last_diagonal", & + test_svd_zero_on_last_diagonal), & + new_unittest("svd_zero_on_inner_diagonal", & + test_svd_zero_on_inner_diagonal) & + ] + end subroutine collect_svd + + !> @brief Test the bidiagonal decomposition (first phase of SVD) + !! + !! This test verifies the correctness of the **bidiagonal decomposition** step, + !! which is the first phase of the SVD algorithm. It checks that the input matrix + !! can be reduced to bidiagonal form using orthogonal matrices P and Qt, and that + !! these matrices are orthogonal. The test does **not** perform full diagonalization. + !! The focus is on ensuring that the reduction to bidiagonal form is correct and + !! reversible for various matrix shapes (square, tall, wide). + !< + subroutine test_bidiagonal_decomposition(error) + ! -- dummy + type(error_type), allocatable, intent(out) :: error + ! -- locals + real(DP), dimension(4, 3) :: A1, A1_reconstructed, A1_mod ! A1: more rows than columns (4x3) + real(DP), dimension(3, 4) :: A2, A2_reconstructed, A2_mod ! A2: more columns than rows (3x4) + real(DP), dimension(4, 4) :: A3, A3_reconstructed, A3_mod ! A3: square matrix (4x4) + + real(DP), dimension(:, :), allocatable :: P1, Qt1 + real(DP), dimension(:, :), allocatable :: P2, Qt2 + real(DP), dimension(:, :), allocatable :: P3, Qt3 + + ! - Arrange. + A1 = reshape( & + [1.0_DP, 0.0_DP, 1.0_DP, 1.0_DP, & + 0.0_DP, 1.0_DP, 1.0_DP, 1.0_DP, & + 0.0_DP, 1.0_DP, 1.0_DP, 0.0_DP & + ], [4, 3]) + A1_mod = A1 + + A2 = reshape( & + [1.0_DP, 0.0_DP, 0.0_DP, & + 0.0_DP, 1.0_DP, 1.0_DP, & + 1.0_DP, 1.0_DP, 1.0_DP, & + 1.0_DP, 1.0_DP, 0.0_DP & + ], [3, 4]) + A2_mod = A2 + + A3 = reshape( & + [1.0_DP, 0.0_DP, 0.0_DP, 1.0_DP, & + 0.0_DP, 1.0_DP, 1.0_DP, 1.0_DP, & + 1.0_DP, 1.0_DP, 1.0_DP, 1.0_DP, & + 1.0_DP, 1.0_DP, 0.0_DP, 0.0_DP & + ], [4, 4]) + A3_mod = A3 + + ! - Act. + call bidiagonal_decomposition(A1_mod, P1, Qt1) + call bidiagonal_decomposition(A2_mod, P2, Qt2) + call bidiagonal_decomposition(A3_mod, P3, Qt3) + + ! - Assert. + ! Test A1 + ! A1_reconstructed = P1 * A1_mod * Qt1 + A1_reconstructed = matmul(P1, matmul(A1_mod, Qt1)) + call check_same_matrix(error, A1_reconstructed, A1) + if (allocated(error)) return + + call check_matrix_row_orthogonal(error, P1) + call check_matrix_column_orthogonal(error, P1) + if (allocated(error)) return + + call check_matrix_row_orthogonal(error, Qt1) + call check_matrix_column_orthogonal(error, Qt1) + if (allocated(error)) return + + call check_matrix_bidiagonal(error, A1_mod) + if (allocated(error)) return + + ! Test A2 + ! A2_reconstructed = P2 * A2_mod * Qt2 + A2_reconstructed = matmul(P2, matmul(A2_mod, Qt2)) + call check_same_matrix(error, A2_reconstructed, A2) + if (allocated(error)) return + + call check_matrix_row_orthogonal(error, P2) + call check_matrix_column_orthogonal(error, P2) + if (allocated(error)) return + + call check_matrix_row_orthogonal(error, Qt2) + call check_matrix_column_orthogonal(error, Qt2) + if (allocated(error)) return + + call check_matrix_bidiagonal(error, A2_mod) + if (allocated(error)) return + + ! Test A3 + ! A3_reconstructed = P3 * A3_mod * Qt3 + A3_reconstructed = matmul(P3, matmul(A3_mod, Qt3)) + call check_same_matrix(error, A3_reconstructed, A3) + if (allocated(error)) return + + call check_matrix_row_orthogonal(error, P3) + call check_matrix_column_orthogonal(error, P3) + if (allocated(error)) return + + call check_matrix_row_orthogonal(error, Qt3) + call check_matrix_column_orthogonal(error, Qt3) + if (allocated(error)) return + + call check_matrix_bidiagonal(error, A3_mod) + if (allocated(error)) return + end subroutine test_bidiagonal_decomposition + + !> @brief Test the diagonalization of a bidiagonal matrix (second phase of SVD) + !! + !! This test verifies the correctness of the **diagonalization** step, + !! which is the second phase of the SVD algorithm. It assumes the input matrix + !! is already in bidiagonal form and applies the QR algorithm to diagonalize it. + !! The test checks that the resulting matrix is diagonal, that the orthogonal + !! matrices U and Vt remain orthogonal, and that the original matrix can be + !! reconstructed. This test ensures the full SVD process (from bidiagonal to diagonal) + !! is correct. + !< + subroutine test_diagonalize_matrix(error) + ! -- dummy + type(error_type), allocatable, intent(out) :: error + ! -- locals + real(DP), dimension(4, 4) :: A1, A1_reconstructed, A1_mod + real(DP), dimension(3, 4) :: A2, A2_reconstructed, A2_mod + real(DP), dimension(4, 3) :: A3, A3_reconstructed, A3_mod + real(DP), dimension(:, :), allocatable :: U1, Vt1 + real(DP), dimension(:, :), allocatable :: U2, Vt2 + real(DP), dimension(:, :), allocatable :: U3, Vt3 + + ! - Arrange. + A1 = reshape( & + [1.0_DP, 0.0_DP, 0.0_DP, 0.0_DP, & + 1.0_DP, 1.0_DP, 0.0_DP, 0.0_DP, & + 0.0_DP, 1.0_DP, 1.0_DP, 0.0_DP, & + 0.0_DP, 0.0_DP, 1.0_DP, 1.0_DP & + ], [4, 4]) + A1_mod = A1 + U1 = Eye(4) + Vt1 = Eye(4) + + A2 = reshape( & + [1.0_DP, 0.0_DP, 0.0_DP, & + 1.0_DP, 1.0_DP, 0.0_DP, & + 0.0_DP, 1.0_DP, 1.0_DP, & + 0.0_DP, 0.0_DP, 1.0_DP & + ], [3, 4]) + A2_mod = A2 + U2 = Eye(3) + Vt2 = Eye(4) + + A3 = reshape( & + [1.0_DP, 0.0_DP, 0.0_DP, 0.0_DP, & + 1.0_DP, 1.0_DP, 0.0_DP, 0.0_DP, & + 0.0_DP, 1.0_DP, 1.0_DP, 0.0_DP & + ], [4, 3]) + A3_mod = A3 + U3 = Eye(4) + Vt3 = Eye(3) + + ! - Act. + call bidiagonal_qr_decomposition(A1_mod, U1, Vt1) + call bidiagonal_qr_decomposition(A2_mod, U2, Vt2) + call bidiagonal_qr_decomposition(A3_mod, U3, Vt3) + + ! - Assert. + ! Test A1 + ! A1_reconstructed = U1 * A1_mod * Vt1 + A1_reconstructed = matmul(U1, matmul(A1_mod, Vt1)) + call check_same_matrix(error, A1_reconstructed, A1) + if (allocated(error)) return + + call check_matrix_row_orthogonal(error, U1) + call check_matrix_column_orthogonal(error, U1) + if (allocated(error)) return + + call check_matrix_row_orthogonal(error, Vt1) + call check_matrix_column_orthogonal(error, Vt1) + if (allocated(error)) return + + call check_matrix_bidiagonal(error, A1_mod) + if (allocated(error)) return + + ! Test A2 + ! A2_reconstructed = U2 * A2_mod * Vt2 + A2_reconstructed = matmul(U2, matmul(A2_mod, Vt2)) + call check_same_matrix(error, A2_reconstructed, A2) + if (allocated(error)) return + + call check_matrix_row_orthogonal(error, U2) + call check_matrix_column_orthogonal(error, U2) + if (allocated(error)) return + + call check_matrix_row_orthogonal(error, Vt2) + call check_matrix_column_orthogonal(error, Vt2) + if (allocated(error)) return + + call check_matrix_bidiagonal(error, A2_mod) + if (allocated(error)) return + + ! Test A3 + ! A3_reconstructed = U3 * A3_mod * Vt3 + A3_reconstructed = matmul(U3, matmul(A3_mod, Vt3)) + call check_same_matrix(error, A3_reconstructed, A3) + if (allocated(error)) return + + call check_matrix_row_orthogonal(error, U3) + call check_matrix_column_orthogonal(error, U3) + if (allocated(error)) return + + call check_matrix_row_orthogonal(error, Vt3) + call check_matrix_column_orthogonal(error, Vt3) + if (allocated(error)) return + + call check_matrix_bidiagonal(error, A3_mod) + if (allocated(error)) return + + end subroutine test_diagonalize_matrix + + !> @brief Test the SVD decomposition of a bidiagonal matrix + !! + !! This test checks that the SVD decomposition works correctly when the input matrix + !! is already in bidiagonal form. It verifies that the SVD reconstructs the original + !! matrix and that the resulting S matrix is diagonal. This ensures that the SVD + !! implementation can handle bidiagonal matrices as input and produce correct results. + !< + subroutine test_svd_input_bidiagonal(error) + ! -- dummy + type(error_type), allocatable, intent(out) :: error + ! -- locals + real(DP), dimension(4, 4) :: A, A_reconstructed + real(DP), DIMENSION(:, :), allocatable :: U, S, Vt + + ! - Arrange. + A = reshape( & + [1.0_DP, 0.0_DP, 0.0_DP, 0.0_DP, & + 1.0_DP, 2.0_DP, 0.0_DP, 0.0_DP, & + 0.0_DP, 1.0_DP, 3.0_DP, 0.0_DP, & + 0.0_DP, 0.0_DP, 1.0_DP, 4.0_DP & + ], [4, 4]) + + ! - Act. + call SVD(A, U, S, Vt) + + ! - Assert. + A_reconstructed = matmul(U, matmul(S, Vt)) + call check_same_matrix(error, A_reconstructed, A) + if (allocated(error)) return + + end subroutine test_svd_input_bidiagonal + + !> @brief Test the SVD decomposition of a square matrix + !! + !! This test checks that the SVD decomposition of a square matrix is correct. + !! The input matrix is chosen such that the returned S-matrix + !! has elements on the entire diagonal. + !< + subroutine test_svd_input_square_matrix(error) + ! -- dummy + type(error_type), allocatable, intent(out) :: error + ! -- locals + real(DP), dimension(4, 4) :: A, A_reconstructed + real(DP), DIMENSION(:, :), allocatable :: U, S, Vt + + ! - Arrange. + A = reshape( & + [1.0_DP, 0.0_DP, 0.0_DP, 1.0_DP, & + 0.0_DP, 1.0_DP, 0.0_DP, 0.0_DP, & + 1.0_DP, 0.0_DP, 1.0_DP, 0.0_DP, & + 1.0_DP, 1.0_DP, 1.0_DP, 1.0_DP & + ], [4, 4]) + + ! - Act. + call SVD(A, U, S, Vt) + + ! - Assert. + A_reconstructed = matmul(U, matmul(S, Vt)) + call check_same_matrix(error, A_reconstructed, A) + if (allocated(error)) return + + call check_matrix_is_diagonal(error, S) + if (allocated(error)) return + + end subroutine test_svd_input_square_matrix + + !> @brief Test the SVD decomposition of a tall matrix (m > n) + !! + !! This test checks that the SVD decomposition works correctly for a **tall matrix**, + !! i.e., a matrix with more rows than columns (m > n). It verifies that the SVD + !! reconstructs the original matrix, and that the resulting S matrix is diagonal. + !! This ensures the SVD implementation handles tall matrices as expected. + !< + subroutine test_svd_input_tall_matrix(error) + ! -- dummy + type(error_type), allocatable, intent(out) :: error + ! -- locals + real(DP), dimension(4, 3) :: A, A_reconstructed + real(DP), DIMENSION(:, :), allocatable :: U, S, Vt + + ! - Arrange. + A = reshape( & + [1.0_DP, 0.0_DP, 0.0_DP, 1.0_DP, & + 0.0_DP, 1.0_DP, 0.0_DP, 1.0_DP, & + 1.0_DP, 0.0_DP, 1.0_DP, 1.0_DP & + ], [4, 3]) + + ! - Act. + call SVD(A, U, S, Vt) + + ! - Assert. + A_reconstructed = matmul(U, matmul(S, Vt)) + call check_same_matrix(error, A_reconstructed, A) + if (allocated(error)) return + + call check_matrix_is_diagonal(error, S) + if (allocated(error)) return + + end subroutine test_svd_input_tall_matrix + + !> @brief Test the SVD decomposition of a wide matrix (n > m) + !! + !! This test checks that the SVD decomposition works correctly for a **wide matrix**, + !! i.e., a matrix with more columns than rows (n > m). It verifies that the SVD + !! reconstructs the original matrix, and that the resulting S matrix is diagonal. + !! This ensures the SVD implementation handles wide matrices as expected. + !< + subroutine test_svd_input_wide_matrix(error) + ! -- dummy + type(error_type), allocatable, intent(out) :: error + ! -- locals + real(DP), dimension(3, 5) :: A, A_reconstructed + real(DP), DIMENSION(:, :), allocatable :: U, S, Vt + + ! - Arrange. + A = reshape( & + [1.0_DP, 0.0_DP, 1.0_DP, & + 0.0_DP, 1.0_DP, 0.0_DP, & + 1.0_DP, 0.0_DP, 0.0_DP, & + 1.0_DP, 1.0_DP, 1.0_DP, & + 1.0_DP, 1.0_DP, 0.0_DP & + ], [3, 5]) + + ! - Act. + call SVD(A, U, S, Vt) + + ! - Assert. + A_reconstructed = matmul(U, matmul(S, Vt)) + call check_same_matrix(error, A_reconstructed, A) + if (allocated(error)) return + + call check_matrix_is_diagonal(error, S) + if (allocated(error)) return + end subroutine test_svd_input_wide_matrix + + !> @brief Test the SVD decomposition of a matrix with a zero on the last diagonal element + !! + !! This test checks that the SVD decomposition works correctly when the input matrix + !! has a zero on its last diagonal element. This situation can cause the matrix to be + !! split into smaller blocks during the SVD process. The test verifies that the SVD + !! reconstructs the original matrix and that the resulting S matrix is diagonal. + !! This ensures the implementation handles matrices with trailing zero diagonals robustly. + !< + subroutine test_svd_zero_on_last_diagonal(error) + ! -- dummy + type(error_type), allocatable, intent(out) :: error + ! -- locals + real(DP), dimension(4, 4) :: A, A_reconstructed + real(DP), DIMENSION(:, :), allocatable :: U, S, Vt + + A = reshape( & + [1.0_DP, 1.0_DP, 0.0_DP, 0.0_DP, & + 1.0_DP, 1.0_DP, 0.0_DP, 0.0_DP, & + 0.0_DP, 1.0_DP, 1.0_DP, 0.0_DP, & + 0.0_DP, 0.0_DP, 1.0_DP, 0.0_DP & + ], [4, 4]) + ! - Act. + call SVD(A, U, S, Vt) + + ! - Assert. + A_reconstructed = matmul(U, matmul(S, Vt)) + call check_same_matrix(error, A_reconstructed, A) + if (allocated(error)) return + + call check_matrix_is_diagonal(error, S) + if (allocated(error)) return + + end subroutine test_svd_zero_on_last_diagonal + + !> @brief Test the SVD decomposition of a matrix with a zero on an inner diagonal element + !! + !! This test checks that the SVD decomposition works correctly when the input matrix + !! has a zero on an inner (not last) diagonal element. This case is more challenging, + !! as it requires the SVD algorithm to correctly split and process the matrix into + !! independent blocks. The test verifies that the SVD reconstructs the original matrix + !! and that the resulting S matrix is diagonal. This ensures the implementation is robust + !! for matrices with inner zero diagonals. + !< + subroutine test_svd_zero_on_inner_diagonal(error) + ! -- dummy + type(error_type), allocatable, intent(out) :: error + ! -- locals + real(DP), dimension(6, 6) :: A, A_reconstructed + real(DP), DIMENSION(:, :), allocatable :: U, S, Vt + + A = reshape( & + [1.0_DP, 0.0_DP, 0.0_DP, 0.0_DP, 0.0_DP, 0.0_DP, & + 1.0_DP, 1.0_DP, 0.0_DP, 0.0_DP, 0.0_DP, 0.0_DP, & + 0.0_DP, 1.0_DP, 0.0_DP, 0.0_DP, 0.0_DP, 0.0_DP, & + 0.0_DP, 0.0_DP, 1.0_DP, 1.0_DP, 0.0_DP, 0.0_DP, & + 0.0_DP, 0.0_DP, 0.0_DP, 1.0_DP, 1.0_DP, 0.0_DP, & + 0.0_DP, 0.0_DP, 0.0_DP, 0.0_DP, 1.0_DP, 1.0_DP & + ], [6, 6]) + ! - Act. + call SVD(A, U, S, Vt) + + ! - Assert. + A_reconstructed = matmul(U, matmul(S, Vt)) + call check_same_matrix(error, A_reconstructed, A) + if (allocated(error)) return + + call check_matrix_is_diagonal(error, S) + if (allocated(error)) return + + end subroutine test_svd_zero_on_inner_diagonal + +end module TestSVD diff --git a/autotest/meson.build b/autotest/meson.build index d4b99b93ae2..681fc0eb6dc 100644 --- a/autotest/meson.build +++ b/autotest/meson.build @@ -13,8 +13,10 @@ if test_drive.found() and not fc_id.contains('intel') 'MemoryContainerIterator', 'MemoryStore', 'Message', + 'PseudoInverse', 'PtrHashTable', 'Sim', + 'SVD', 'SwfUtils', 'TimeSelect', 'TimeStepSelect', @@ -23,6 +25,7 @@ if test_drive.found() and not fc_id.contains('intel') test_srcs = files( 'tester.f90', + 'tester_utils.f90', ) foreach t : tests test_srcs += files('Test@0@.f90'.format(t.underscorify())) diff --git a/autotest/tester.f90 b/autotest/tester.f90 index d2ecba16014..52beff50d62 100644 --- a/autotest/tester.f90 +++ b/autotest/tester.f90 @@ -14,8 +14,10 @@ program tester use TestMemoryContainerIterator, only: collect_memorycontaineriterator use TestMemoryStore, only: collect_memorystore use TestMessage, only: collect_message + use TestPseudoInverse, only: collect_pinv use TestPtrHashTable, only: collect_ptrhashtable use TestSim, only: collect_sim + use TestSVD, only: collect_svd use TestSwfUtils, only: collect_swfutils use TestTimeSelect, only: collect_timeselect use TestTimeStepSelect, only: collect_timestepselect @@ -42,8 +44,10 @@ program tester collect_memorycontaineriterator), & new_testsuite("MemoryStore", collect_memorystore), & new_testsuite("Message", collect_message), & + new_testsuite("PseudoInverse", collect_pinv), & new_testsuite("PtrHashTable", collect_ptrhashtable), & new_testsuite("Sim", collect_sim), & + new_testsuite("SVD", collect_svd), & new_testsuite("SwfUtils", collect_swfutils), & new_testsuite("TimeSelect", collect_timeselect), & new_testsuite("TimeStepSelect", collect_timestepselect), & diff --git a/autotest/tester_utils.f90 b/autotest/tester_utils.f90 new file mode 100644 index 00000000000..d06633bb1e9 --- /dev/null +++ b/autotest/tester_utils.f90 @@ -0,0 +1,184 @@ +module TesterUtils + use testdrive, only: error_type, test_failed, check + use KindModule, only: DP, I4B + use ConstantsModule, only: DSAME + + implicit none + private + + public :: check_same_matrix + public :: check_matrix_row_orthogonal + public :: check_matrix_column_orthogonal + public :: check_matrix_is_diagonal + public :: check_matrix_bidiagonal + +contains + !> @brief Check if two matrices A and B are the same within a given tolerance. + !! + !! This subroutine compares two matrices element-wise and sets the + !! error flag if any corresponding elements differ by more than the + !! specified tolerance. + !< + subroutine check_same_matrix(error, A, B, tolerance) + ! -- dummy + type(error_type), allocatable, intent(out) :: error + real(DP), dimension(:, :), intent(in) :: A, B + real(DP), optional :: tolerance + ! -- locals + integer(I4B) :: i, j + real(DP) :: tol + logical :: is_different + + if (.not. present(tolerance)) then + tol = DSAME + else + tol = tolerance + end if + + if (size(A, 1) /= size(B, 1) .or. size(A, 2) /= size(B, 2)) then + call test_failed(error, "Matrices have different dimensions") + return + end if + + is_different = .false. + do i = 1, size(A, 1) + do j = 1, size(A, 2) + if (abs(A(i, j) - B(i, j)) > tol) then + call test_failed(error, "Matrices have different elements") + return + end if + end do + end do + + end subroutine check_same_matrix + + !> @brief Checks if the rows of a matrix are orthogonal within a given tolerance. + !! + !! This subroutine verifies whether the rows of the input matrix `A` are orthogonal + !! to each other within a specified tolerance. If the rows are not orthogonal, the + !! subroutine sets the `error` flag to true. + !< + subroutine check_matrix_row_orthogonal(error, A, tolerance) + ! -- dummy + type(error_type), allocatable, intent(out) :: error + real(DP), dimension(:, :), intent(in) :: A + real(DP), optional :: tolerance + ! -- locals + integer(I4B) :: i, j + real(DP) :: tol + real(DP) :: err + + if (.not. present(tolerance)) then + tol = DSAME + else + tol = tolerance + end if + + do i = 1, size(A, 1) - 1 + do j = i + 1, size(A, 1) + err = abs(dot_product(A(i, :), A(j, :))) + call check(error, err, 0.0_DP, thr=DSAME) + if (allocated(error)) return + end do + end do + end subroutine check_matrix_row_orthogonal + + !> @brief Checks if the columns of a matrix are orthogonal within a given tolerance. + !! + !! This subroutine verifies whether the columns of the input matrix `A` are orthogonal + !! to each other within a specified tolerance. If the columns are not orthogonal, the + !! subroutine sets the `error` flag to true. + !< + subroutine check_matrix_column_orthogonal(error, A, tolerance) + ! -- dummy + type(error_type), allocatable, intent(out) :: error + real(DP), dimension(:, :), intent(in) :: A + real(DP), optional :: tolerance + ! -- locals + integer(I4B) :: i, j + real(DP) :: tol + real(DP) :: err + + if (.not. present(tolerance)) then + tol = DSAME + else + tol = tolerance + end if + + do i = 1, size(A, 2) - 1 + do j = i + 1, size(A, 2) + err = abs(dot_product(A(:, i), A(:, j))) + call check(error, err, 0.0_DP, thr=DSAME) + if (allocated(error)) return + end do + end do + end subroutine check_matrix_column_orthogonal + + !> @brief Checks if the given matrix A is diagonal within a specified tolerance. + !! + !! This subroutine verifies whether the input matrix A is diagonal by + !! comparing its off-diagonal elements to zero within a given tolerance. + !< + subroutine check_matrix_is_diagonal(error, A, tolerance) + ! -- dummy + type(error_type), allocatable, intent(out) :: error + real(DP), dimension(:, :), intent(in) :: A + real(DP), optional :: tolerance + ! -- locals + integer(I4B) :: i, j + real(DP) :: tol + real(DP) :: err + + if (.not. present(tolerance)) then + tol = DSAME + else + tol = tolerance + end if + + do i = 1, size(A, 1) + do j = 1, size(A, 2) + if (i /= j) then + err = abs(A(i, j)) + call check(error, err, 0.0_DP, thr=DSAME) + if (allocated(error)) return + end if + end do + end do + + end subroutine check_matrix_is_diagonal + + !> @brief Checks if the given matrix A is bidiagonal within a specified tolerance. + !! + !! This subroutine verifies whether the input matrix A is bidiagonal by + !! comparing its off-diagonal and off-superdiagonal elements to zero within + !! a given tolerance. + !< + subroutine check_matrix_bidiagonal(error, A, tolerance) + ! -- dummy + type(error_type), allocatable, intent(out) :: error + real(DP), dimension(:, :), intent(in) :: A + real(DP), optional :: tolerance + ! -- locals + integer(I4B) :: i, j + real(DP) :: tol + real(DP) :: err + + if (.not. present(tolerance)) then + tol = DSAME + else + tol = tolerance + end if + + do i = 1, size(A, 1) + do j = 1, size(A, 2) + if (i /= j .and. i + 1 /= j) then + err = abs(A(i, j)) + call check(error, err, 0.0_DP, thr=DSAME) + if (allocated(error)) return + end if + end do + end do + + end subroutine check_matrix_bidiagonal + +end module TesterUtils diff --git a/make/makefile b/make/makefile index d8d96b88d48..48b64c4916b 100644 --- a/make/makefile +++ b/make/makefile @@ -23,31 +23,30 @@ 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/Timing -SOURCEDIR26=../src/Utilities -SOURCEDIR27=../src/Utilities/ArrayRead -SOURCEDIR28=../src/Utilities/Export -SOURCEDIR29=../src/Utilities/Idm -SOURCEDIR30=../src/Utilities/Idm/mf6blockfile -SOURCEDIR31=../src/Utilities/Idm/netcdf -SOURCEDIR32=../src/Utilities/Libraries -SOURCEDIR33=../src/Utilities/Libraries/blas -SOURCEDIR34=../src/Utilities/Libraries/daglib -SOURCEDIR35=../src/Utilities/Libraries/rcm -SOURCEDIR36=../src/Utilities/Libraries/sparsekit -SOURCEDIR37=../src/Utilities/Libraries/sparskit2 -SOURCEDIR38=../src/Utilities/Matrix -SOURCEDIR39=../src/Utilities/Memory -SOURCEDIR40=../src/Utilities/Observation -SOURCEDIR41=../src/Utilities/OutputControl -SOURCEDIR42=../src/Utilities/Performance -SOURCEDIR43=../src/Utilities/TimeSeries -SOURCEDIR44=../src/Utilities/Vector +SOURCEDIR20=../src/Solution +SOURCEDIR21=../src/Solution/LinearMethods +SOURCEDIR22=../src/Solution/PETSc +SOURCEDIR23=../src/Solution/ParticleTracker +SOURCEDIR24=../src/Timing +SOURCEDIR25=../src/Utilities +SOURCEDIR26=../src/Utilities/ArrayRead +SOURCEDIR27=../src/Utilities/Export +SOURCEDIR28=../src/Utilities/Idm +SOURCEDIR29=../src/Utilities/Idm/mf6blockfile +SOURCEDIR30=../src/Utilities/Idm/netcdf +SOURCEDIR31=../src/Utilities/Libraries +SOURCEDIR32=../src/Utilities/Libraries/blas +SOURCEDIR33=../src/Utilities/Libraries/daglib +SOURCEDIR34=../src/Utilities/Libraries/rcm +SOURCEDIR35=../src/Utilities/Libraries/sparsekit +SOURCEDIR36=../src/Utilities/Libraries/sparskit2 +SOURCEDIR37=../src/Utilities/Matrix +SOURCEDIR38=../src/Utilities/Memory +SOURCEDIR39=../src/Utilities/Observation +SOURCEDIR40=../src/Utilities/OutputControl +SOURCEDIR41=../src/Utilities/Performance +SOURCEDIR42=../src/Utilities/TimeSeries +SOURCEDIR43=../src/Utilities/Vector VPATH = \ ${SOURCEDIR1} \ @@ -92,8 +91,7 @@ ${SOURCEDIR39} \ ${SOURCEDIR40} \ ${SOURCEDIR41} \ ${SOURCEDIR42} \ -${SOURCEDIR43} \ -${SOURCEDIR44} +${SOURCEDIR43} .SUFFIXES: .f90 .F90 .o @@ -306,8 +304,6 @@ $(OBJDIR)/SimStages.o \ $(OBJDIR)/NumericalModel.o \ $(OBJDIR)/OutputControlData.o \ $(OBJDIR)/gwf-ic.o \ -$(OBJDIR)/tsp-fmi.o \ -$(OBJDIR)/IInterpolationScheme.o \ $(OBJDIR)/UzfEtUtil.o \ $(OBJDIR)/Xt3dAlgorithm.o \ $(OBJDIR)/TvBase.o \ @@ -330,14 +326,11 @@ $(OBJDIR)/SwfCxsUtils.o \ $(OBJDIR)/VirtualModel.o \ $(OBJDIR)/BaseExchange.o \ $(OBJDIR)/TspSpc.o \ +$(OBJDIR)/tsp-fmi.o \ $(OBJDIR)/GweInputData.o \ $(OBJDIR)/OutputControl.o \ $(OBJDIR)/tsp-ic.o \ -$(OBJDIR)/UpwindScheme.o \ -$(OBJDIR)/TVDScheme.o \ $(OBJDIR)/TspAdvOptions.o \ -$(OBJDIR)/CentralDifferenceScheme.o \ -$(OBJDIR)/AdvSchemeEnum.o \ $(OBJDIR)/UzfCellGroup.o \ $(OBJDIR)/Xt3dInterface.o \ $(OBJDIR)/gwf-tvk.o \ @@ -516,10 +509,13 @@ $(OBJDIR)/RunControlFactory.o \ $(OBJDIR)/IdmLoad.o \ $(OBJDIR)/ConnectionBuilder.o \ $(OBJDIR)/comarg.o \ +$(OBJDIR)/LinearAlgebraUtils.o \ $(OBJDIR)/mf6core.o \ +$(OBJDIR)/SVD.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 d6ed378136a..9e8db9778da 100644 --- a/msvs/mf6core.vfproj +++ b/msvs/mf6core.vfproj @@ -655,6 +655,7 @@ + @@ -664,6 +665,7 @@ + @@ -675,6 +677,7 @@ + diff --git a/src/Utilities/LinearAlgebraUtils.f90 b/src/Utilities/LinearAlgebraUtils.f90 new file mode 100644 index 00000000000..b62ea1e0d13 --- /dev/null +++ b/src/Utilities/LinearAlgebraUtils.f90 @@ -0,0 +1,51 @@ +module LinearAlgebraUtilsModule + use KindModule, only: DP, I4B + + implicit none + private + + public :: eye, outer_product, cross_product + +contains + + pure function eye(n) RESULT(A) + ! -- dummy + integer(I4B), intent(IN) :: n + real(DP), dimension(n, n) :: A + ! -- locals + integer(I4B) :: i + + A = 0.0_DP + do i = 1, n + A(i, i) = 1.0_DP + end do + + end function eye + + pure function outer_product(A, B) result(AB) + ! -- dummy + real(DP), intent(in) :: A(:), B(:) + real(DP) :: AB(size(A), size(B)) + ! -- locals + integer :: nA, nB + + nA = size(A) + nB = size(B) + + AB = spread(source=A, dim=2, ncopies=nB) * & + spread(source=B, dim=1, ncopies=nA) + end function outer_product + + function cross_product(a, b) result(c) + ! -- return + real(DP), dimension(3) :: c + ! -- dummy + real(DP), dimension(3), intent(in) :: a + real(DP), dimension(3), intent(in) :: b + + c(1) = a(2) * b(3) - a(3) * b(2) + c(2) = a(3) * b(1) - a(1) * b(3) + c(3) = a(1) * b(2) - a(2) * b(1) + end function cross_product + +end module LinearAlgebraUtilsModule diff --git a/src/Utilities/PseudoInverse.f90 b/src/Utilities/PseudoInverse.f90 new file mode 100644 index 00000000000..cfd0460fecc --- /dev/null +++ b/src/Utilities/PseudoInverse.f90 @@ -0,0 +1,55 @@ +module PseudoInverseModule + use KindModule, only: DP, I4B + use ConstantsModule, only: DPREC + + use SVDModule, only: SVD + + implicit none + private + + public :: pinv + +contains + + !!> @brief Computes the Moore-Penrose pseudoinverse of a matrix using SVD. + !! + !! This function computes the pseudoinverse \(A^+\) of the input matrix \(A\) using + !! the singular value decomposition (SVD). Small singular values (less than `DPREC`) + !! are set to zero to improve numerical stability. The pseudoinverse is useful for + !! solving least-squares problems and for handling rank-deficient or non-square matrices. + !! + !! The result satisfies: + !! B = V * Sigma^+ * U^T + !! where \(A = U \Sigma V^T\) is the SVD of \(A\), and \(\Sigma^+\) is the diagonal matrix + !! with reciprocals of the nonzero singular values. + !< + function pinv(A) result(B) + ! -- dummy + real(DP), intent(in) :: A(:, :) + real(DP) :: B(SIZE(A, DIM=2), SIZE(A, DIM=1)) !! The pseudoinverse of A (size n x m, where A is m x n) + ! -- locals + integer(I4B) :: m, n + integer(I4B) :: pos + real(DP), dimension(:, :), allocatable :: U + real(DP), dimension(:, :), allocatable :: Vt + real(DP), dimension(:, :), allocatable :: Sigma + + m = size(A, dim=1) + n = size(A, dim=2) + + CALL SVD(A, U, Sigma, Vt) + + ! Transform Sigma to its Sigma^+ + do pos = 1, min(m, n) + if (dabs(Sigma(pos, pos)) < DPREC) then + Sigma(pos, pos) = 0.0_DP + else + Sigma(pos, pos) = 1.0_DP / Sigma(pos, pos) + end if + end do + + B = matmul(transpose(Vt), matmul(transpose(Sigma), transpose(U))) + + end function pinv + +end module PseudoInverseModule diff --git a/src/Utilities/SVD.f90 b/src/Utilities/SVD.f90 new file mode 100644 index 00000000000..546cfb7bc92 --- /dev/null +++ b/src/Utilities/SVD.f90 @@ -0,0 +1,469 @@ +module SVDModule + use KindModule, only: DP, LGP, I4B + use ConstantsModule, only: DONE, DPREC, DSAME, DEM6 + use LinearAlgebraUtilsModule, only: eye, outer_product + + implicit none + private + + public :: SVD + public :: bidiagonal_decomposition + public :: bidiagonal_qr_decomposition + +contains + + !!> @brief Constructs a Householder reflection matrix for a given vector. + !! + !< + pure function house_holder_matrix(x) result(Q) + ! -- dummy + real(DP), intent(in) :: x(:) + real(DP), allocatable, dimension(:, :) :: Q + ! -- locals + real(DP) :: x_norm, y_norm + real(DP), allocatable :: v(:) + real(DP), allocatable :: w(:) + + x_norm = norm2(x) + y_norm = norm2(x(2:)) + + Q = Eye(size(x)) + + if (dabs(y_norm) < DSAME) then + return + end if + + v = x + if (x(1) > 0) then + ! Parlett (1971) suggested this modification to avoid cancellation + v(1) = -(y_norm**2) / (X(1) + x_norm) + else + v(1) = v(1) - x_norm + end if + + w = v / norm2(v) + + Q = Q - 2.0_DP * outer_product(w, w) + + end function house_holder_matrix + + !!> @brief Reduces a matrix to bidiagonal form using Householder transformations. + !! + !! This subroutine transforms the input matrix `A` into a bidiagonal matrix by applying a sequence of + !! Householder reflections from the left and right. The orthogonal transformations are accumulated in + !! the matrices `P` (left transformations) and `Qt` (right transformations). + !! + !! The result is: + !! P^T * A * Qt = B + !! where `B` is bidiagonal, `P` and `Qt` are orthogonal matrices. + !! + !! This step is the first phase of the Golub-Reinsch SVD algorithm and prepares the matrix for further + !! diagonalization using the QR algorithm. + !! + !< + pure subroutine bidiagonal_decomposition(A, P, Qt) + ! -- dummy + real(DP), intent(inout), dimension(:, :) :: A + real(DP), intent(out), dimension(:, :), allocatable :: P, Qt + ! -- locals + integer(I4B) :: M, N + integer(I4B) :: I + real(DP), allocatable, dimension(:, :) :: Qi, Pi + real(DP), allocatable, dimension(:) :: h + + M = size(A, dim=1) ! Number of rows + N = size(A, dim=2) ! Number of columns + + Qt = Eye(N) + P = Eye(M) + + do I = 1, min(M, N) + ! columns + h = A(I:M, I) + Pi = house_holder_matrix(h) + A(I:M, :) = matmul(Pi, A(I:M, :)) ! Apply householder transformation from left + P(:, I:M) = matmul(P(:, I:M), Pi) + + ! rows + if (I < N) then + h = A(I, I + 1:N) + Qi = transpose(house_holder_matrix(h)) + A(:, I + 1:N) = matmul(A(:, I + 1:N), Qi) ! Apply householder transformation from right + Qt(I + 1:N, :) = matmul(Qi, Qt(I + 1:N, :)) + end if + + end do + + end subroutine bidiagonal_decomposition + + !!> @brief Constructs a Givens rotation matrix. + !! + !< + pure function givens_rotation(a, b) result(G) + ! -- dummy + real(DP), intent(in) :: a, b + real(DP), dimension(2, 2) :: G + ! -- locals + real(DP) :: c, s, h, d + + if (abs(b) < DPREC) then + G = Eye(2) + return + end if + + h = hypot(a, b) + d = 1.0 / h + c = abs(a) * d + s = sign(d, a) * b + + G(1, 1) = c + G(1, 2) = s + G(2, 1) = -s + G(2, 2) = c + + end function givens_rotation + + !!> @brief Computes the Wilkinson shift for the lower-right 2x2 block of a bidiagonal matrix. + !! + !! This function calculates the Wilkinson shift, which is used to accelerate convergence in the implicit QR algorithm + !! for bidiagonal matrices during SVD. The shift is computed from the lower-right 2x2 block of the input matrix `A` + !! and is chosen to be the eigenvalue of that block which is closer to the bottom-right element. + !! + !! The shift is used to improve numerical stability and speed up the annihilation of off-diagonal elements. + !! + !< + pure function compute_shift(A) result(mu) + ! -- dummy + real(DP), intent(in), dimension(:, :) :: A + real(DP) :: mu + ! -- locals + integer(I4B) :: m, n, min_mn + real(DP) T11, T12, T21, T22 + real(DP) dm, fmmin, fm, dn + real(DP) :: mean, product, mean_product, mu1, mu2 + + m = size(A, dim=1) ! Number of rows + n = size(A, dim=2) ! Number of columns + + min_mn = min(m, n) + + dn = A(min_mn, min_mn) + dm = A(min_mn - 1, min_mn - 1) + fm = A(min_mn - 1, min_mn) + if (min_mn > 2) then + fmmin = A(min_mn - 2, min_mn - 1) + else + fmmin = 0.0_DP + end if + + T11 = dm**2 + fmmin**2 + T12 = dm * fm + T21 = T12 + T22 = fm**2 + dn**2 + + mean = (T11 + T22) / 2.0_DP + product = T11 * T22 - T12 * T21 + mean_product = mean**2 - product + if (abs(mean_product) < DEM6) then + mean_product = 0.0_DP + end if + mu1 = mean - sqrt(mean_product) + mu2 = mean + sqrt(mean_product) + if (abs(T22 - mu1) < abs(T22 - mu2)) then + mu = mu1 + else + mu = mu2 + end if + + end function compute_shift + + !!> @brief Performs QR iterations on a bidiagonal matrix as part of the SVD algorithm. + !! + !< + pure subroutine bidiagonal_qr_decomposition(A, U, VT) + ! -- dummy + real(DP), intent(inout), dimension(:, :) :: A + real(DP), intent(inout), dimension(:, :) :: U, Vt + ! -- locals + integer(I4B) :: m, n, I, J + real(DP), dimension(:, :), allocatable :: G + real(DP) :: t11, t12, mu + + m = size(A, dim=1) ! Number of rows + n = size(A, dim=2) ! Number of columns + + DO I = 1, n - 1 + J = I + 1 + if (I == 1) then + ! For the first iteration use the implicit shift + mu = compute_shift(A) + ! Apply the shift to the full tri-diagonal matrix. + t11 = A(1, 1)**2 - mu + t12 = A(1, 1) * A(1, 2) + + G = givens_rotation(t11, t12) + else + G = givens_rotation(A(I - 1, I), A(I - 1, J)) + end if + + A(:, I:J) = matmul(A(:, I:J), transpose(G)) + Vt(I:J, :) = matmul(G, Vt(I:J, :)) + + if (j > m) cycle + G = givens_rotation(A(I, I), A(J, I)) + A(I:J, :) = matmul(G, A(I:J, :)) + U(:, I:J) = matmul(U(:, I:J), transpose(G)) + END DO + + end subroutine bidiagonal_qr_decomposition + + !!> @brief Handles zero (or near-zero) diagonal elements in a bidiagonal matrix during SVD. + !! + !! This subroutine detects and processes zero (or near-zero) diagonal elements in the bidiagonal matrix `A`. + !! If a zero is found on the diagonal, it applies a sequence of Givens rotations to either zero out the + !! corresponding column (if the zero is at the end of the diagonal) or the corresponding row (otherwise). + !! The orthogonal transformations are accumulated in the `U` and `VT` matrices as appropriate. + !! This step is essential for splitting the matrix into smaller blocks and ensuring numerical stability + !! during the iterative QR step of the SVD algorithm. + !! + !< + pure subroutine handle_zero_diagonal(A, U, VT) + ! -- dummy + real(DP), intent(inout), dimension(:, :) :: A + real(DP), intent(inout), dimension(:, :) :: U, Vt + ! -- locals + integer(I4B) :: m, n, I + real(DP), dimension(:, :), allocatable :: G + integer(I4B) :: zero_index + + m = size(A, dim=1) ! Number of rows + n = size(A, dim=2) ! Number of columns + + do I = 1, min(m, n) + if (abs(A(I, I)) < DPREC) then + zero_index = I + exit + end if + end do + + if (zero_index == min(n, m)) then + ! If the zero index is the last element of the diagonal then zero out the column + do I = zero_index - 1, 1, -1 + G = givens_rotation(A(I, I), A(I, zero_index)) + A(:, [I, zero_index]) = matmul(A(:, [I, zero_index]), transpose(G)) + Vt([I, zero_index], :) = matmul(G, Vt([I, zero_index], :)) + end do + + else + ! Else zero out the row + do I = zero_index + 1, n + G = givens_rotation(A(I, I), A(zero_index, I)) + A([zero_index, I], :) = matmul(transpose(G), A([zero_index, I], :)) + U(:, [zero_index, I]) = matmul(U(:, [zero_index, I]), G) + end do + end if + end subroutine handle_zero_diagonal + + !!> @brief Computes the infinity norm of the superdiagonal elements of a matrix. + !! + !< + pure function superdiagonal_norm(A) result(norm) + ! -- dummy + real(DP), intent(in) :: A(:, :) + real(DP) :: norm + ! -- locals + integer(I4B) :: m, n, I + + m = size(A, DIM=1) ! Number of rows + n = size(A, DIM=2) ! Number of columns + + norm = 0.0_DP + do I = 1, min(m, n) - 1 + norm = max(norm, abs(A(I, I + 1))) + end do + + end function superdiagonal_norm + + !!> @brief Finds the range of nonzero superdiagonal elements in a bidiagonal matrix. + !! + !! This subroutine scans the superdiagonal of the matrix `A` and determines the indices `p` and `q` + !! such that all superdiagonal elements outside the range `A(p:q, p:q)` are (near-)zero, while + !! those within the range may be nonzero. This is useful for identifying active submatrices during + !! the iterative QR step of the SVD algorithm. + !! + !! The search starts from the lower-right corner and moves upward to find the last nonzero + !! superdiagonal (sets `q`), then moves upward again to find the first zero (sets `p`). + !< + pure subroutine find_nonzero_superdiagonal(A, p, q) + ! -- dummy + real(DP), intent(in), dimension(:, :) :: A + integer(I4B), intent(out) :: p ! Index of the first nonzero superdiagonal (start of active block) + integer(I4B), intent(out) :: q ! Index of the last nonzero superdiagonal (end of active block) + ! -- locals + integer(I4B) :: m, n, j, min_mn + + m = size(A, DIM=1) ! Number of rows + n = size(A, DIM=2) ! Number of columns + + min_mn = min(m, n) + p = 1 + q = min_mn + + do j = min_mn, 2, -1 + if (abs(A(j - 1, j)) > DPREC) then + q = j + exit + end if + end do + + do j = q - 1, 2, -1 + if (abs(A(j - 1, j)) < DPREC) then + p = j + exit + end if + end do + end subroutine find_nonzero_superdiagonal + + !!> @brief Checks if a matrix has a (near-)zero diagonal element. + !! + !< + pure function has_zero_diagonal(A) result(has_zero) + ! -- dummy + real(DP), intent(in) :: A(:, :) + logical(LGP) :: has_zero + ! -- locals + integer(I4B) :: m, n, I + + m = size(A, dim=1) ! Number of rows + n = size(A, dim=2) ! Number of columns + + has_zero = .FALSE. + do I = 1, min(m, n) + if (abs(A(I, I)) < DPREC) then + has_zero = .TRUE. + exit + end if + end do + + end function has_zero_diagonal + + !!> @brief Reduces a rectangular matrix to square form using Givens rotations. + !! + !! This subroutine transforms a rectangular matrix `A` into a square matrix by applying Givens rotations + !! to eliminate extra columns (when the number of columns exceeds the number of rows). + !! The corresponding orthogonal transformations are accumulated in the matrix `Qt`. + !< + pure subroutine make_matrix_square(A, Qt) + ! -- dummy + real(DP), intent(inout), dimension(:, :) :: A + real(DP), intent(inout), dimension(:, :), allocatable :: Qt + ! -- locals + real(DP), dimension(:, :), allocatable :: G + integer(I4B) :: m, n, I + + m = size(A, dim=1) ! Number of rows + n = size(A, dim=2) ! Number of columns + + do I = m, 1, -1 + G = givens_rotation(A(I, I), A(I, M + 1)) + A(:, [I, M + 1]) = matmul(A(:, [I, M + 1]), transpose(G)) + Qt([I, M + 1], :) = matmul(G, Qt([I, M + 1], :)) + end do + + end subroutine make_matrix_square + + !!> @brief Sets small superdiagonal elements to zero for numerical stability. + !! + !! This subroutine scans the superdiagonal elements of the matrix `A` and sets to zero any element + !! whose absolute value is less than or equal to a tolerance (relative to the neighboring diagonal elements). + !! This helps to maintain numerical stability and ensures that the matrix is properly treated as bidiagonal + !! in subsequent SVD steps. + !! + !! The criterion used is: + !! |A(j, j+1)| <= DPREC * (|A(j, j)| + |A(j+1, j+1)|) + !! In such cases, A(j, j+1) is set to zero. + !< + pure subroutine clean_superdiagonal(A) + ! -- dummy + real(DP), intent(inout), dimension(:, :) :: A + ! -- locals + integer(I4B) :: n, m, j + + m = size(A, DIM=1) ! Number of rows + n = size(A, DIM=2) ! Number of columns + + do j = 1, min(n, m) - 1 + if (abs(A(j, j + 1)) <= DPREC * (abs(A(j, j)) + abs(A(j + 1, j + 1)))) then + A(j, j + 1) = 0.0_DP + end if + end do + + end subroutine clean_superdiagonal + + !> @brief Singular Value Decomposition + !! + !! This method decomposes the matrix A into U, S and VT. + !! It follows the algorithm as described by Golub and Reinsch. + !! + !! The first step is to decompose the matrix A into a bidiagonal matrix. + !! This is done using Householder transformations. + !! Then second step is to decompose the bidiagonal matrix into U, S and VT + !! by repetitively applying the QR algorithm. + !! If there is a zero on the diagonal or superdiagonal the matrix can be split + !! into two smaller matrices and the QR algorithm can be applied to the smaller + !! matrices. + !! + !! The matrix U is the eigenvectors of A*A^T + !! The matrix VT is the eigenvectors of A^T*A + !! The matrix S is the square root of the eigenvalues of A*A^T or A^T*A + !! + !< + pure subroutine SVD(A, U, S, VT) + ! -- dummy + real(DP), intent(in), dimension(:, :) :: A + real(DP), intent(out), dimension(:, :), allocatable :: U + real(DP), intent(out), dimension(:, :), allocatable :: S + real(DP), intent(out), dimension(:, :), allocatable :: VT + ! -- locals + integer(I4B) :: i ! Loop counter for QR iterations + integer(I4B) :: max_itr ! Maximum number of QR iterations allowed for convergence + real(DP) :: error ! Convergence criterion (superdiagonal norm) + + integer(I4B) :: m ! Number of rows in input matrix A + integer(I4B) :: n ! Number of columns in input matrix A + + integer(I4B) :: r ! Start index of active submatrix (for QR step) + integer(I4B) :: q ! End index of active submatrix (for QR step) + + max_itr = 100 + + m = size(A, dim=1) ! Number of rows + n = size(A, dim=2) ! Number of columns + S = A + + call bidiagonal_decomposition(S, U, Vt) + + if (n > m) then + call make_matrix_square(S, Vt) + end if + + do i = 1, max_itr + call find_nonzero_superdiagonal(S, r, q) + + ! find zero diagonal + if (has_zero_diagonal(S(r:q, r:q))) then + ! write(*,*) 'Iteration: ', i, ' handle zero diagonal element' + call handle_zero_diagonal(S(r:q, r:q), U(:, r:q), Vt(r:q, :)) + else + call bidiagonal_qr_decomposition(S(r:q, r:q), U(:, r:q), Vt(r:q, :)) + call clean_superdiagonal(S(r:q, r:q)) + + error = superdiagonal_norm(S) + ! write(*,*) 'Iteration: ', i, ' Error: ', error + if (error < DPREC) exit + end if + end do + + end subroutine SVD + +end module diff --git a/src/meson.build b/src/meson.build index b25f151c9fc..4a214cc4662 100644 --- a/src/meson.build +++ b/src/meson.build @@ -412,6 +412,7 @@ modflow_sources = files( 'Utilities' / 'KeyValueList.f90', 'Utilities' / 'KeyValueListIterator.f90', 'Utilities' / 'KeyValueNode.f90', + 'Utilities' / 'LinearAlgebraUtils.f90', 'Utilities' / 'List.f90', 'Utilities' / 'ListIterator.f90', 'Utilities' / 'ListNode.f90', @@ -421,6 +422,7 @@ modflow_sources = files( 'Utilities' / 'Message.f90', 'Utilities' / 'OpenSpec.f90', 'Utilities' / 'PackageBudget.f90', + 'Utilities' / 'PseudoInverse.f90', 'Utilities' / 'PtrHashTable.f90', 'Utilities' / 'PtrHashTableIterator.f90', 'Utilities' / 'Sim.f90', @@ -432,6 +434,7 @@ modflow_sources = files( 'Utilities' / 'STLStackInt.f90', 'Utilities' / 'STLVecInt.f90', 'Utilities' / 'StringList.f90', + 'Utilities' / 'SVD.f90', 'Utilities' / 'Table.f90', 'Utilities' / 'TableTerm.f90', 'Utilities' / 'Timer.f90', From dfeb9324d885f596dba8b0235c37ee52088fe066 Mon Sep 17 00:00:00 2001 From: Sunny Titus Date: Mon, 30 Jun 2025 13:30:24 +0200 Subject: [PATCH 03/15] Add the new unstructured TVD limiter --- .codespell.ignore | 1 + autotest/test_gwt_disu01.py | 8 +- make/makefile | 71 +++--- msvs/mf6core.vfproj | 7 +- src/Model/Discretization/DisUtils.f90 | 110 ++++++++++ .../TransportModel/Gradient/IGradient.f90 | 38 ++++ .../Gradient/LeastSquaredGradient.f90 | 170 +++++++++++++++ .../InterpolationScheme/AdvSchemeEnum.f90 | 1 + .../InterpolationScheme/UTVDScheme.f90 | 204 ++++++++++++++++++ src/Model/TransportModel/tsp-adv.f90 | 14 ++ src/meson.build | 4 + 11 files changed, 596 insertions(+), 32 deletions(-) create mode 100644 src/Model/Discretization/DisUtils.f90 create mode 100644 src/Model/TransportModel/Gradient/IGradient.f90 create mode 100644 src/Model/TransportModel/Gradient/LeastSquaredGradient.f90 create mode 100644 src/Model/TransportModel/InterpolationScheme/UTVDScheme.f90 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_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/make/makefile b/make/makefile index 48b64c4916b..4172f7b088a 100644 --- a/make/makefile +++ b/make/makefile @@ -23,30 +23,32 @@ SOURCEDIR16=../src/Model/OverlandFlow SOURCEDIR17=../src/Model/ParticleTracking SOURCEDIR18=../src/Model/SurfaceWaterFlow SOURCEDIR19=../src/Model/TransportModel -SOURCEDIR20=../src/Solution -SOURCEDIR21=../src/Solution/LinearMethods -SOURCEDIR22=../src/Solution/PETSc -SOURCEDIR23=../src/Solution/ParticleTracker -SOURCEDIR24=../src/Timing -SOURCEDIR25=../src/Utilities -SOURCEDIR26=../src/Utilities/ArrayRead -SOURCEDIR27=../src/Utilities/Export -SOURCEDIR28=../src/Utilities/Idm -SOURCEDIR29=../src/Utilities/Idm/mf6blockfile -SOURCEDIR30=../src/Utilities/Idm/netcdf -SOURCEDIR31=../src/Utilities/Libraries -SOURCEDIR32=../src/Utilities/Libraries/blas -SOURCEDIR33=../src/Utilities/Libraries/daglib -SOURCEDIR34=../src/Utilities/Libraries/rcm -SOURCEDIR35=../src/Utilities/Libraries/sparsekit -SOURCEDIR36=../src/Utilities/Libraries/sparskit2 -SOURCEDIR37=../src/Utilities/Matrix -SOURCEDIR38=../src/Utilities/Memory -SOURCEDIR39=../src/Utilities/Observation -SOURCEDIR40=../src/Utilities/OutputControl -SOURCEDIR41=../src/Utilities/Performance -SOURCEDIR42=../src/Utilities/TimeSeries -SOURCEDIR43=../src/Utilities/Vector +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/Timing +SOURCEDIR27=../src/Utilities +SOURCEDIR28=../src/Utilities/ArrayRead +SOURCEDIR29=../src/Utilities/Export +SOURCEDIR30=../src/Utilities/Idm +SOURCEDIR31=../src/Utilities/Idm/mf6blockfile +SOURCEDIR32=../src/Utilities/Idm/netcdf +SOURCEDIR33=../src/Utilities/Libraries +SOURCEDIR34=../src/Utilities/Libraries/blas +SOURCEDIR35=../src/Utilities/Libraries/daglib +SOURCEDIR36=../src/Utilities/Libraries/rcm +SOURCEDIR37=../src/Utilities/Libraries/sparsekit +SOURCEDIR38=../src/Utilities/Libraries/sparskit2 +SOURCEDIR39=../src/Utilities/Matrix +SOURCEDIR40=../src/Utilities/Memory +SOURCEDIR41=../src/Utilities/Observation +SOURCEDIR42=../src/Utilities/OutputControl +SOURCEDIR43=../src/Utilities/Performance +SOURCEDIR44=../src/Utilities/TimeSeries +SOURCEDIR45=../src/Utilities/Vector VPATH = \ ${SOURCEDIR1} \ @@ -91,7 +93,9 @@ ${SOURCEDIR39} \ ${SOURCEDIR40} \ ${SOURCEDIR41} \ ${SOURCEDIR42} \ -${SOURCEDIR43} +${SOURCEDIR43} \ +${SOURCEDIR44} \ +${SOURCEDIR45} .SUFFIXES: .f90 .F90 .o @@ -284,6 +288,7 @@ $(OBJDIR)/TrackControl.o \ $(OBJDIR)/TimeSelect.o \ $(OBJDIR)/prt-fmi.o \ $(OBJDIR)/TimeStepSelect.o \ +$(OBJDIR)/LinearAlgebraUtils.o \ $(OBJDIR)/SfrCrossSectionUtils.o \ $(OBJDIR)/TernarySolveTrack.o \ $(OBJDIR)/SubcellTri.o \ @@ -293,6 +298,8 @@ $(OBJDIR)/VirtualBase.o \ $(OBJDIR)/STLVecInt.o \ $(OBJDIR)/BaseModel.o \ $(OBJDIR)/PrintSaveManager.o \ +$(OBJDIR)/tsp-fmi.o \ +$(OBJDIR)/SVD.o \ $(OBJDIR)/SfrCrossSectionManager.o \ $(OBJDIR)/dag_module.o \ $(OBJDIR)/BoundaryPackageExt.o \ @@ -304,6 +311,10 @@ $(OBJDIR)/SimStages.o \ $(OBJDIR)/NumericalModel.o \ $(OBJDIR)/OutputControlData.o \ $(OBJDIR)/gwf-ic.o \ +$(OBJDIR)/IInterpolationScheme.o \ +$(OBJDIR)/IGradient.o \ +$(OBJDIR)/DisUtils.o \ +$(OBJDIR)/PseudoInverse.o \ $(OBJDIR)/UzfEtUtil.o \ $(OBJDIR)/Xt3dAlgorithm.o \ $(OBJDIR)/TvBase.o \ @@ -326,11 +337,16 @@ $(OBJDIR)/SwfCxsUtils.o \ $(OBJDIR)/VirtualModel.o \ $(OBJDIR)/BaseExchange.o \ $(OBJDIR)/TspSpc.o \ -$(OBJDIR)/tsp-fmi.o \ $(OBJDIR)/GweInputData.o \ $(OBJDIR)/OutputControl.o \ $(OBJDIR)/tsp-ic.o \ +$(OBJDIR)/UTVDScheme.o \ +$(OBJDIR)/UpwindScheme.o \ +$(OBJDIR)/TVDScheme.o \ $(OBJDIR)/TspAdvOptions.o \ +$(OBJDIR)/LeastSquaredGradient.o \ +$(OBJDIR)/CentralDifferenceScheme.o \ +$(OBJDIR)/AdvSchemeEnum.o \ $(OBJDIR)/UzfCellGroup.o \ $(OBJDIR)/Xt3dInterface.o \ $(OBJDIR)/gwf-tvk.o \ @@ -509,13 +525,10 @@ $(OBJDIR)/RunControlFactory.o \ $(OBJDIR)/IdmLoad.o \ $(OBJDIR)/ConnectionBuilder.o \ $(OBJDIR)/comarg.o \ -$(OBJDIR)/LinearAlgebraUtils.o \ $(OBJDIR)/mf6core.o \ -$(OBJDIR)/SVD.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 9e8db9778da..1e58bac857b 100644 --- a/msvs/mf6core.vfproj +++ b/msvs/mf6core.vfproj @@ -248,6 +248,7 @@ + @@ -368,12 +369,16 @@ + + + - + + diff --git a/src/Model/Discretization/DisUtils.f90 b/src/Model/Discretization/DisUtils.f90 new file mode 100644 index 00000000000..a125b20f326 --- /dev/null +++ b/src/Model/Discretization/DisUtils.f90 @@ -0,0 +1,110 @@ +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_centroid + 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, fmi, n, m) result(d) + !-- modules + use TspFmiModule, only: TspFmiType + ! -- return + real(DP), dimension(3) :: d + ! -- dummy + class(DisBaseType), intent(in) :: dis + type(TspFmiType), pointer, intent(in) :: fmi + integer(I4B), intent(in) :: n, m + ! -- local + real(DP) :: x_dir, y_dir, z_dir, length + real(DP) :: satn, satm + 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 + ! -- Thi can happen when using an extended stencil (neighbours-of-neigbhours) to compute the gradients + if (isympos == -1) then + xn = cell_centroid(dis, n) + xm = cell_centroid(dis, m) + + d = xm - xn + return + end if + + ! -- Account for the saturation levels if available + ihc = dis%con%ihc(isympos) + if (associated(fmi%gwfsat)) then + satn = fmi%gwfsat(n) + satm = fmi%gwfsat(m) + else + satn = DONE + satm = DONE + end if + + ! -- Get the connection direction and length + call dis%connection_vector(n, m, .true., satn, satm, 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 centroid coordinates of a given cell. + !! + !! This function computes the centroid (geometric center) of cell `n` in the discretization. + !! The centroid 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_centroid(dis, n) result(centroid) + ! -- dummy + class(DisBaseType), intent(in) :: dis + integer(I4B), intent(in) :: n + real(DP), dimension(3) :: centroid + + centroid(1) = dis%xc(n) + centroid(2) = dis%yc(n) + centroid(3) = (dis%top(n) + dis%bot(n)) / 2.0_dp + end function cell_centroid +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/LeastSquaredGradient.f90 b/src/Model/TransportModel/Gradient/LeastSquaredGradient.f90 new file mode 100644 index 00000000000..9fb02d3c714 --- /dev/null +++ b/src/Model/TransportModel/Gradient/LeastSquaredGradient.f90 @@ -0,0 +1,170 @@ +module LeastSquaredGradientModule + use KindModule, only: DP, I4B + use ConstantsModule, only: DONE + + Use IGradient + use BaseDisModule, only: DisBaseType + use TspFmiModule, only: TspFmiType + use PseudoInverseModule, only: pinv + use DisUtilsModule, only: number_connected_faces, node_distance + + implicit none + private + + public :: LeastSquaredGradientType + + 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 operator 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 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) :: LeastSquaredGradientType + class(DisBaseType), pointer :: dis + type(TspFmiType), pointer :: fmi + type(Array2D), allocatable, dimension(:) :: grad_op + contains + procedure :: get + + procedure, private :: compute_cell_gradient + procedure, private :: create_grad_operator + end type LeastSquaredGradientType + + interface LeastSquaredGradientType + module procedure Constructor + end interface LeastSquaredGradientType + +contains + function constructor(dis, fmi) Result(gradient) + ! --dummy + class(DisBaseType), pointer, intent(in) :: dis + type(TspFmiType), pointer, intent(in) :: fmi + !-- return + type(LeastSquaredGradientType) :: gradient + ! -- local + integer(I4B) :: n, nodes + + gradient%dis => dis + gradient%fmi => fmi + + nodes = dis%nodes + + ! -- Compute the gradient operator + nodes = dis%nodes + allocate (gradient%grad_op(dis%nodes)) + do n = 1, nodes + gradient%grad_op(n)%data = gradient%create_grad_operator(n) + end do + end function constructor + + function create_grad_operator(this, n) result(grad_op) + ! -- dummy + class(LeastSquaredGradientType) :: this + integer(I4B), intent(in) :: n ! Cell index for which to create the operator + real(DP), dimension(:, :), allocatable :: grad_op ! The resulting gradient operator (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 :: d_trans ! Transpose of d (3 x number_connections) + 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 + ! the cell center and its neighbor. + real(DP), dimension(3, 3) :: g + real(DP), dimension(3, 3) :: g_inv + + number_connections = number_connected_faces(this%dis, n) + + allocate (d(number_connections, 3)) + allocate (d_trans(3, number_connections)) + allocate (grad_op(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, this%fmi, 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 + + d_trans = transpose(d) + + ! Compute the G and inverse G matrices + g = matmul(d_trans, d) + g_inv = pinv(g) + + ! Compute the gradient operator + grad_op = matmul(matmul(g_inv, d_trans), grad_scale) + + end function create_grad_operator + + function get(this, n, c) result(grad_c) + ! -- dummy + class(LeastSquaredGradientType), 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(LeastSquaredGradientType), target :: this + integer(I4B), intent(in) :: n + real(DP), dimension(:), intent(in) :: phi_new + ! -- local + real(DP), dimension(:, :), pointer :: grad_op + 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 + grad_op => this%grad_op(n)%data + grad_c = matmul(grad_op, dc) + + end function compute_cell_gradient + +end module LeastSquaredGradientModule diff --git a/src/Model/TransportModel/InterpolationScheme/AdvSchemeEnum.f90 b/src/Model/TransportModel/InterpolationScheme/AdvSchemeEnum.f90 index 3e6ee17b1b3..1289114d10c 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 diff --git a/src/Model/TransportModel/InterpolationScheme/UTVDScheme.f90 b/src/Model/TransportModel/InterpolationScheme/UTVDScheme.f90 new file mode 100644 index 00000000000..6067e6d78ac --- /dev/null +++ b/src/Model/TransportModel/InterpolationScheme/UTVDScheme.f90 @@ -0,0 +1,204 @@ +module UTVDSchemeModule + use KindModule, only: DP, I4B + use ConstantsModule, only: DONE, DZERO, DSAME, DHALF + use IInterpolationSchemeModule, only: IInterpolationSchemeType, 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(IInterpolationSchemeType) :: 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, this%fmi, 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), target :: 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 09412fe657a..f4bf933e17b 100644 --- a/src/Model/TransportModel/tsp-adv.f90 +++ b/src/Model/TransportModel/tsp-adv.f90 @@ -7,12 +7,16 @@ module TspAdvModule use TspFmiModule, only: TspFmiType use TspAdvOptionsModule, only: TspAdvOptionsType use MatrixBaseModule, only: MatrixBaseType + ! -- Gradient schemes + use IGradient, only: IGradientType + use LeastSquaredGradientModule, only: LeastSquaredGradientType ! -- Interpolation schemes use IInterpolationSchemeModule, only: IInterpolationSchemeType, CoefficientsType use AdvSchemeEnumModule use UpwindSchemeModule, only: UpwindSchemeType use CentralDifferenceSchemeModule, only: CentralDifferenceSchemeType use TVDSchemeModule, only: TVDSchemeType + use UTVDSchemeModule, only: UTVDSchemeType implicit none private @@ -27,6 +31,7 @@ module TspAdvModule real(DP), pointer :: eqnsclfac => null() !< governing equation scale factor; =1. for solute; =rhow*cpw for energy class(IInterpolationSchemeType), allocatable :: face_interpolation !< interpolation scheme for face values + class(IGradientType), allocatable :: gradient !< cell centered gradient contains procedure :: adv_df @@ -122,6 +127,9 @@ subroutine adv_ar(this, dis, ibound) this%dis => dis this%ibound => ibound ! + ! -- Compute the gradient operator + this%gradient = LeastSquaredGradientType(this%dis, this%fmi) + ! ! -- Create interpolation scheme iadvwt_value = this%iadvwt ! Dereference iadvwt to work with case statement select case (iadvwt_value) @@ -134,6 +142,9 @@ 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%face_interpolation = & + UTVDSchemeType(this%dis, this%fmi, this%gradient) case default call store_error("Unknown advection scheme", terminate=.TRUE.) end select @@ -366,6 +377,9 @@ subroutine read_options(this) case ('TVD') this%iadvwt = ADV_SCHEME_TVD write (this%iout, fmtiadvwt) 'TVD' + case ('UTVD') + this%iadvwt = ADV_SCHEME_UTVD + write (this%iout, fmtiadvwt) 'UTVD' case default write (errmsg, '(a, a)') & 'Unknown scheme: ', trim(keyword) diff --git a/src/meson.build b/src/meson.build index 4a214cc4662..fa28eb200f0 100644 --- a/src/meson.build +++ b/src/meson.build @@ -158,6 +158,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', @@ -267,11 +268,14 @@ modflow_sources = files( 'Model' / 'ModelUtilities' / 'VectorInterpolation.f90', 'Model' / 'ModelUtilities' / 'Xt3dAlgorithm.f90', 'Model' / 'ModelUtilities' / 'Xt3dInterface.f90', + 'Model' / 'TransportModel' / 'Gradient' / 'IGradient.f90', + 'Model' / 'TransportModel' / 'Gradient' / 'LeastSquaredGradient.f90', 'Model' / 'TransportModel' / 'InterpolationScheme' / 'AdvSchemeEnum.f90', 'Model' / 'TransportModel' / 'InterpolationScheme' / 'CentralDifferenceScheme.f90', 'Model' / 'TransportModel' / 'InterpolationScheme' / 'IInterpolationScheme.f90', 'Model' / 'TransportModel' / 'InterpolationScheme' / 'TVDScheme.f90', 'Model' / 'TransportModel' / 'InterpolationScheme' / 'UpwindScheme.f90', + 'Model' / 'TransportModel' / 'InterpolationScheme' / 'UTVDScheme.f90', 'Model' / 'TransportModel' / 'tsp.f90', 'Model' / 'TransportModel' / 'tsp-adv.f90', 'Model' / 'TransportModel' / 'tsp-apt.f90', From 1306626c635930ea80f834d42118bb476a437667 Mon Sep 17 00:00:00 2001 From: Sunny Titus Date: Tue, 1 Jul 2025 22:16:36 +0200 Subject: [PATCH 04/15] Fix merge --- make/makefile | 2 - msvs/mf6core.vfproj | 3 +- .../InterpolationScheme/UpwindScheme.f90 | 60 ------------------- 3 files changed, 2 insertions(+), 63 deletions(-) delete mode 100644 src/Model/TransportModel/InterpolationScheme/UpwindScheme.f90 diff --git a/make/makefile b/make/makefile index 22847fbf7c4..e65fcd8fcf7 100644 --- a/make/makefile +++ b/make/makefile @@ -533,7 +533,6 @@ $(OBJDIR)/IdmLoad.o \ $(OBJDIR)/ConnectionBuilder.o \ $(OBJDIR)/comarg.o \ $(OBJDIR)/mf6core.o \ -$(OBJDIR)/IInterpolationScheme.o \ $(OBJDIR)/BaseGeometry.o \ $(OBJDIR)/mf6.o \ $(OBJDIR)/StringList.o \ @@ -543,7 +542,6 @@ $(OBJDIR)/sparsekit.o \ $(OBJDIR)/rcm.o \ $(OBJDIR)/blas1_d.o \ $(OBJDIR)/Iunit.o \ -$(OBJDIR)/UpwindScheme.o \ $(OBJDIR)/GwfSfrCommon.o \ $(OBJDIR)/gwf-sfr-transient.o \ $(OBJDIR)/gwf-sfr-steady.o \ diff --git a/msvs/mf6core.vfproj b/msvs/mf6core.vfproj index 0aa809fba01..74cad34f4c7 100644 --- a/msvs/mf6core.vfproj +++ b/msvs/mf6core.vfproj @@ -377,7 +377,8 @@ - + + diff --git a/src/Model/TransportModel/InterpolationScheme/UpwindScheme.f90 b/src/Model/TransportModel/InterpolationScheme/UpwindScheme.f90 deleted file mode 100644 index 8e14b7caedb..00000000000 --- a/src/Model/TransportModel/InterpolationScheme/UpwindScheme.f90 +++ /dev/null @@ -1,60 +0,0 @@ -module UpwindSchemeModule - use KindModule, only: DP, I4B - use ConstantsModule, only: DZERO - use IInterpolationSchemeModule, only: IInterpolationSchemeType, CoefficientsType - use BaseDisModule, only: DisBaseType - use TspFmiModule, only: TspFmiType - - implicit none - private - - public :: UpwindSchemeType - - type, extends(IInterpolationSchemeType) :: UpwindSchemeType - private - class(DisBaseType), pointer :: dis - type(TspFmiType), pointer :: fmi - contains - procedure :: compute - end type UpwindSchemeType - - interface UpwindSchemeType - module procedure constructor - end interface UpwindSchemeType - -contains - function constructor(dis, fmi) result(interpolation_scheme) - ! -- return - type(UpwindSchemeType) :: interpolation_scheme - ! --dummy - class(DisBaseType), pointer, intent(in) :: dis - type(TspFmiType), pointer, intent(in) :: fmi - - interpolation_scheme%dis => dis - interpolation_scheme%fmi => fmi - - end function constructor - - function compute(this, n, m, iposnm, phi) result(phi_face) - !-- return - type(CoefficientsType) :: phi_face - ! -- dummy - class(UpwindSchemeType), target :: this - integer(I4B), intent(in) :: n - integer(I4B), intent(in) :: m - integer(I4B), intent(in) :: iposnm - real(DP), intent(in), dimension(:) :: phi - ! -- local - real(DP) :: qnm - - ! -- Compute the coefficients for the upwind scheme - qnm = this%fmi%gwfflowja(iposnm) - - if (qnm < DZERO) then - phi_face%c_n = 1.0_dp - else - phi_face%c_m = 1.0_dp - end if - - end function compute -end module UpwindSchemeModule From fc0aa31f1af7e998b78a8d2acdc925dbb1a1719d Mon Sep 17 00:00:00 2001 From: Sunny Titus Date: Wed, 2 Jul 2025 14:20:12 +0200 Subject: [PATCH 05/15] Add UTVD tests to existing unit tests. Fix correct stencil size on exchanges --- autotest/test_gwt_adv01.py | 114 +++++++++- autotest/test_gwt_adv01_fmi.py | 114 +++++++++- autotest/test_gwt_adv01_gwtgwt.py | 114 +++++++++- autotest/test_gwt_adv02.py | 212 +++++++++++++++++- autotest/test_gwt_adv02_gwtgwt.py | 4 +- autotest/test_gwt_adv03.py | 51 ++++- autotest/test_gwt_adv04.py | 4 +- autotest/test_gwt_henry_gwtgwt.py | 9 +- autotest/test_par_gwt_henry.py | 2 +- doc/mf6io/mf6ivar/dfn/gwt-adv.dfn | 4 +- src/Exchange/exg-gwegwe.f90 | 9 +- src/Exchange/exg-gwtgwt.f90 | 4 +- src/Model/Connection/GweGweConnection.f90 | 6 +- src/Model/Connection/GwtGwtConnection.f90 | 6 +- .../InterpolationScheme/AdvSchemeEnum.f90 | 2 +- 15 files changed, 619 insertions(+), 36 deletions(-) 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_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/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/TransportModel/InterpolationScheme/AdvSchemeEnum.f90 b/src/Model/TransportModel/InterpolationScheme/AdvSchemeEnum.f90 index 1289114d10c..4efb6f3428b 100644 --- a/src/Model/TransportModel/InterpolationScheme/AdvSchemeEnum.f90 +++ b/src/Model/TransportModel/InterpolationScheme/AdvSchemeEnum.f90 @@ -9,4 +9,4 @@ module AdvSchemeEnumModule 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 From bf2dcc7250442ebfb5928de001863a9f80279b52 Mon Sep 17 00:00:00 2001 From: Sunny Titus Date: Mon, 7 Jul 2025 14:23:56 +0200 Subject: [PATCH 06/15] Apply review comments --- make/makefile | 2 +- msvs/mf6core.vfproj | 2 +- src/Model/Discretization/DisUtils.f90 | 19 ++---- ...dGradient.f90 => LeastSquaresGradient.f90} | 65 +++++++++---------- .../InterpolationScheme/UTVDScheme.f90 | 2 +- src/Model/TransportModel/tsp-adv.f90 | 6 +- src/meson.build | 2 +- 7 files changed, 41 insertions(+), 57 deletions(-) rename src/Model/TransportModel/Gradient/{LeastSquaredGradient.f90 => LeastSquaresGradient.f90} (72%) diff --git a/make/makefile b/make/makefile index e65fcd8fcf7..9eefd285af8 100644 --- a/make/makefile +++ b/make/makefile @@ -351,7 +351,7 @@ $(OBJDIR)/UTVDScheme.o \ $(OBJDIR)/UpstreamScheme.o \ $(OBJDIR)/TVDScheme.o \ $(OBJDIR)/TspAdvOptions.o \ -$(OBJDIR)/LeastSquaredGradient.o \ +$(OBJDIR)/LeastSquaresGradient.o \ $(OBJDIR)/CentralDifferenceScheme.o \ $(OBJDIR)/AdvSchemeEnum.o \ $(OBJDIR)/UzfCellGroup.o \ diff --git a/msvs/mf6core.vfproj b/msvs/mf6core.vfproj index 74cad34f4c7..c1c3a8ec450 100644 --- a/msvs/mf6core.vfproj +++ b/msvs/mf6core.vfproj @@ -371,7 +371,7 @@ - + diff --git a/src/Model/Discretization/DisUtils.f90 b/src/Model/Discretization/DisUtils.f90 index a125b20f326..be9fa939e16 100644 --- a/src/Model/Discretization/DisUtils.f90 +++ b/src/Model/Discretization/DisUtils.f90 @@ -35,18 +35,16 @@ end function number_connected_faces !! 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, fmi, n, m) result(d) + 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 - type(TspFmiType), pointer, intent(in) :: fmi integer(I4B), intent(in) :: n, m ! -- local real(DP) :: x_dir, y_dir, z_dir, length - real(DP) :: satn, satm integer(I4B) :: ipos, isympos, ihc real(DP), dimension(3) :: xn, xm @@ -60,7 +58,7 @@ function node_distance(dis, fmi, n, m) result(d) end do ! -- if the connection is not found, then return the distance between the two nodes - ! -- Thi can happen when using an extended stencil (neighbours-of-neigbhours) to compute the gradients + ! -- This can happen when using an extended stencil (neighbours-of-neigbhours) to compute the gradients if (isympos == -1) then xn = cell_centroid(dis, n) xm = cell_centroid(dis, m) @@ -69,18 +67,9 @@ function node_distance(dis, fmi, n, m) result(d) return end if - ! -- Account for the saturation levels if available - ihc = dis%con%ihc(isympos) - if (associated(fmi%gwfsat)) then - satn = fmi%gwfsat(n) - satm = fmi%gwfsat(m) - else - satn = DONE - satm = DONE - end if - ! -- Get the connection direction and length - call dis%connection_vector(n, m, .true., satn, satm, ihc, x_dir, & + 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 diff --git a/src/Model/TransportModel/Gradient/LeastSquaredGradient.f90 b/src/Model/TransportModel/Gradient/LeastSquaresGradient.f90 similarity index 72% rename from src/Model/TransportModel/Gradient/LeastSquaredGradient.f90 rename to src/Model/TransportModel/Gradient/LeastSquaresGradient.f90 index 9fb02d3c714..e4694aa6b13 100644 --- a/src/Model/TransportModel/Gradient/LeastSquaredGradient.f90 +++ b/src/Model/TransportModel/Gradient/LeastSquaresGradient.f90 @@ -1,17 +1,16 @@ -module LeastSquaredGradientModule +module LeastSquaresGradientModule use KindModule, only: DP, I4B use ConstantsModule, only: DONE Use IGradient use BaseDisModule, only: DisBaseType - use TspFmiModule, only: TspFmiType use PseudoInverseModule, only: pinv use DisUtilsModule, only: number_connected_faces, node_distance implicit none private - public :: LeastSquaredGradientType + public :: LeastSquaresGradientType type Array2D real(DP), dimension(:, :), allocatable :: data @@ -20,9 +19,11 @@ module LeastSquaredGradientModule !> @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 operator using the Moore-Penrose pseudoinverse, + !! 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. @@ -33,49 +34,45 @@ module LeastSquaredGradientModule !! @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) :: LeastSquaredGradientType + type, extends(IGradientType) :: LeastSquaresGradientType class(DisBaseType), pointer :: dis - type(TspFmiType), pointer :: fmi - type(Array2D), allocatable, dimension(:) :: grad_op + type(Array2D), allocatable, dimension(:) :: R ! Gradient reconstruction matrix contains procedure :: get procedure, private :: compute_cell_gradient - procedure, private :: create_grad_operator - end type LeastSquaredGradientType + procedure, private :: create_gradient_reconstruction_matrix + end type LeastSquaresGradientType - interface LeastSquaredGradientType + interface LeastSquaresGradientType module procedure Constructor - end interface LeastSquaredGradientType + end interface LeastSquaresGradientType contains - function constructor(dis, fmi) Result(gradient) + function constructor(dis) Result(gradient) ! --dummy class(DisBaseType), pointer, intent(in) :: dis - type(TspFmiType), pointer, intent(in) :: fmi !-- return - type(LeastSquaredGradientType) :: gradient + type(LeastSquaresGradientType) :: gradient ! -- local integer(I4B) :: n, nodes gradient%dis => dis - gradient%fmi => fmi - nodes = dis%nodes - ! -- Compute the gradient operator + ! -- Compute the gradient rec nodes = dis%nodes - allocate (gradient%grad_op(dis%nodes)) + allocate (gradient%R(dis%nodes)) do n = 1, nodes - gradient%grad_op(n)%data = gradient%create_grad_operator(n) + gradient%R(n)%data = gradient%create_gradient_reconstruction_matrix(n) end do end function constructor - function create_grad_operator(this, n) result(grad_op) + function create_gradient_reconstruction_matrix(this, n) result(R) ! -- dummy - class(LeastSquaredGradientType) :: this + class(LeastSquaresGradientType) :: this integer(I4B), intent(in) :: n ! Cell index for which to create the operator - real(DP), dimension(:, :), allocatable :: grad_op ! The resulting gradient operator (3 x number_connections) + 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 @@ -93,7 +90,7 @@ function create_grad_operator(this, n) result(grad_op) allocate (d(number_connections, 3)) allocate (d_trans(3, number_connections)) - allocate (grad_op(3, number_connections)) + allocate (R(3, number_connections)) allocate (grad_scale(number_connections, number_connections)) grad_scale = 0 @@ -105,7 +102,7 @@ function create_grad_operator(this, n) result(grad_op) 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, this%fmi, n, m) + dnm = node_distance(this%dis, n, m) length = norm2(dnm) d(local_pos, :) = dnm / length @@ -116,18 +113,18 @@ function create_grad_operator(this, n) result(grad_op) d_trans = transpose(d) - ! Compute the G and inverse G matrices + ! Compute the G and pseudo-inverse G matrices g = matmul(d_trans, d) g_inv = pinv(g) - ! Compute the gradient operator - grad_op = matmul(matmul(g_inv, d_trans), grad_scale) + ! Compute the gradient recontructions matrix + R = matmul(matmul(g_inv, d_trans), grad_scale) - end function create_grad_operator + end function create_gradient_reconstruction_matrix function get(this, n, c) result(grad_c) ! -- dummy - class(LeastSquaredGradientType), target :: this + class(LeastSquaresGradientType), target :: this integer(I4B), intent(in) :: n real(DP), dimension(:), intent(in) :: c !-- return @@ -140,11 +137,11 @@ function compute_cell_gradient(this, n, phi_new) result(grad_c) ! -- return real(DP), dimension(3) :: grad_c ! -- dummy - class(LeastSquaredGradientType), target :: this + class(LeastSquaresGradientType), target :: this integer(I4B), intent(in) :: n real(DP), dimension(:), intent(in) :: phi_new ! -- local - real(DP), dimension(:, :), pointer :: grad_op + real(DP), dimension(:, :), pointer :: R integer(I4B) :: ipos, local_pos integer(I4B) :: number_connections @@ -162,9 +159,9 @@ function compute_cell_gradient(this, n, phi_new) result(grad_c) end do ! Compute the cells gradient - grad_op => this%grad_op(n)%data - grad_c = matmul(grad_op, dc) + R => this%R(n)%data + grad_c = matmul(R, dc) end function compute_cell_gradient -end module LeastSquaredGradientModule +end module LeastSquaresGradientModule diff --git a/src/Model/TransportModel/InterpolationScheme/UTVDScheme.f90 b/src/Model/TransportModel/InterpolationScheme/UTVDScheme.f90 index 736b6d96778..8c335df43c6 100644 --- a/src/Model/TransportModel/InterpolationScheme/UTVDScheme.f90 +++ b/src/Model/TransportModel/InterpolationScheme/UTVDScheme.f90 @@ -124,7 +124,7 @@ function compute(this, n, m, iposnm, phi) result(phi_face) grad_c = this%gradient%get(iup, phi) ! ! Darwish's method to compute virtual node concentration - dnm = node_distance(this%dis, this%fmi, iup, idn) + dnm = node_distance(this%dis, iup, idn) c_virtual = phi(idn) - 2.0_dp * (dot_product(grad_c, dnm)) ! ! Enforce local TVD condition. diff --git a/src/Model/TransportModel/tsp-adv.f90 b/src/Model/TransportModel/tsp-adv.f90 index 3248ff773da..61591d0da6b 100644 --- a/src/Model/TransportModel/tsp-adv.f90 +++ b/src/Model/TransportModel/tsp-adv.f90 @@ -9,7 +9,7 @@ module TspAdvModule use MatrixBaseModule, only: MatrixBaseType ! -- Gradient schemes use IGradient, only: IGradientType - use LeastSquaredGradientModule, only: LeastSquaredGradientType + use LeastSquaresGradientModule, only: LeastSquaresGradientType ! -- Interpolation schemes use InterpolationSchemeInterfaceModule, only: InterpolationSchemeInterface, & CoefficientsType @@ -128,9 +128,6 @@ subroutine adv_ar(this, dis, ibound) this%dis => dis this%ibound => ibound ! - ! -- Compute the gradient operator - this%gradient = LeastSquaredGradientType(this%dis, this%fmi) - ! ! -- Create interpolation scheme iadvwt_value = this%iadvwt ! Dereference iadvwt to work with case statement select case (iadvwt_value) @@ -144,6 +141,7 @@ subroutine adv_ar(this, dis, ibound) 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 diff --git a/src/meson.build b/src/meson.build index b5abe674346..674d3bc02a2 100644 --- a/src/meson.build +++ b/src/meson.build @@ -269,7 +269,7 @@ modflow_sources = files( 'Model' / 'ModelUtilities' / 'Xt3dAlgorithm.f90', 'Model' / 'ModelUtilities' / 'Xt3dInterface.f90', 'Model' / 'TransportModel' / 'Gradient' / 'IGradient.f90', - 'Model' / 'TransportModel' / 'Gradient' / 'LeastSquaredGradient.f90', + 'Model' / 'TransportModel' / 'Gradient' / 'LeastSquaresGradient.f90', 'Model' / 'TransportModel' / 'InterpolationScheme' / 'AdvSchemeEnum.f90', 'Model' / 'TransportModel' / 'InterpolationScheme' / 'CentralDifferenceScheme.f90', 'Model' / 'TransportModel' / 'InterpolationScheme' / 'InterpolationSchemeInterface.f90', From bdbffa9b13b3b58f28f3cd41236d7d18f8b9497d Mon Sep 17 00:00:00 2001 From: Sunny Titus Date: Mon, 7 Jul 2025 14:26:02 +0200 Subject: [PATCH 07/15] Fix lint error --- src/Model/TransportModel/Gradient/LeastSquaresGradient.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Model/TransportModel/Gradient/LeastSquaresGradient.f90 b/src/Model/TransportModel/Gradient/LeastSquaresGradient.f90 index e4694aa6b13..10b4f6c2ee1 100644 --- a/src/Model/TransportModel/Gradient/LeastSquaresGradient.f90 +++ b/src/Model/TransportModel/Gradient/LeastSquaresGradient.f90 @@ -117,7 +117,7 @@ function create_gradient_reconstruction_matrix(this, n) result(R) g = matmul(d_trans, d) g_inv = pinv(g) - ! Compute the gradient recontructions matrix + ! Compute the gradient reconstructions matrix R = matmul(matmul(g_inv, d_trans), grad_scale) end function create_gradient_reconstruction_matrix From 49612ab8154cb5760513bf2b61ad2b1810ccbb2e Mon Sep 17 00:00:00 2001 From: Sunny Titus Date: Tue, 8 Jul 2025 11:28:06 +0200 Subject: [PATCH 08/15] Apply review comments --- src/Model/Discretization/DisUtils.f90 | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/src/Model/Discretization/DisUtils.f90 b/src/Model/Discretization/DisUtils.f90 index be9fa939e16..cbf7ea7be7e 100644 --- a/src/Model/Discretization/DisUtils.f90 +++ b/src/Model/Discretization/DisUtils.f90 @@ -7,7 +7,7 @@ module DisUtilsModule private public :: number_connected_faces - public :: cell_centroid + public :: cell_center public :: node_distance contains @@ -60,8 +60,8 @@ function node_distance(dis, n, m) result(d) ! -- 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_centroid(dis, n) - xm = cell_centroid(dis, m) + xn = cell_center(dis, n) + xm = cell_center(dis, m) d = xm - xn return @@ -79,21 +79,21 @@ function node_distance(dis, n, m) result(d) end function node_distance - !> @brief Returns the centroid coordinates of a given cell. + !> @brief Returns the center coordinates of a given cell. !! - !! This function computes the centroid (geometric center) of cell `n` in the discretization. - !! The centroid is returned as a 3-element vector containing the x, y, and z coordinates. + !! 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_centroid(dis, n) result(centroid) + function cell_center(dis, n) result(center) ! -- dummy class(DisBaseType), intent(in) :: dis integer(I4B), intent(in) :: n - real(DP), dimension(3) :: centroid + real(DP), dimension(3) :: center - centroid(1) = dis%xc(n) - centroid(2) = dis%yc(n) - centroid(3) = (dis%top(n) + dis%bot(n)) / 2.0_dp - end function cell_centroid + 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 From aeb5f72d0e3160a4312bdeabd01f8ec58ae62442 Mon Sep 17 00:00:00 2001 From: Sunny Titus Date: Fri, 11 Jul 2025 20:43:14 +0200 Subject: [PATCH 09/15] Apply review comment. Use pinv directly on the distance matrix omstead of the G matrix --- .../Gradient/LeastSquaresGradient.f90 | 13 +------------ 1 file changed, 1 insertion(+), 12 deletions(-) diff --git a/src/Model/TransportModel/Gradient/LeastSquaresGradient.f90 b/src/Model/TransportModel/Gradient/LeastSquaresGradient.f90 index 10b4f6c2ee1..5e78ec0cd41 100644 --- a/src/Model/TransportModel/Gradient/LeastSquaresGradient.f90 +++ b/src/Model/TransportModel/Gradient/LeastSquaresGradient.f90 @@ -79,17 +79,12 @@ function create_gradient_reconstruction_matrix(this, n) result(R) 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 :: d_trans ! Transpose of d (3 x number_connections) 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 - ! the cell center and its neighbor. - real(DP), dimension(3, 3) :: g - real(DP), dimension(3, 3) :: g_inv number_connections = number_connected_faces(this%dis, n) allocate (d(number_connections, 3)) - allocate (d_trans(3, number_connections)) allocate (R(3, number_connections)) allocate (grad_scale(number_connections, number_connections)) @@ -111,14 +106,8 @@ function create_gradient_reconstruction_matrix(this, n) result(R) local_pos = local_pos + 1 end do - d_trans = transpose(d) - - ! Compute the G and pseudo-inverse G matrices - g = matmul(d_trans, d) - g_inv = pinv(g) - ! Compute the gradient reconstructions matrix - R = matmul(matmul(g_inv, d_trans), grad_scale) + R = matmul(pinv(d), grad_scale) end function create_gradient_reconstruction_matrix From 92028a6b1b2cf29d9ad2ce7d3a04c28dc13372be Mon Sep 17 00:00:00 2001 From: Sunny Titus Date: Wed, 30 Jul 2025 14:32:56 +0200 Subject: [PATCH 10/15] Optimize the limiter by caching several values --- .../CentralDifferenceScheme.f90 | 7 ++ .../InterpolationSchemeInterface.f90 | 8 ++ .../InterpolationScheme/TVDScheme.f90 | 6 ++ .../InterpolationScheme/UTVDScheme.f90 | 95 +++++++++++++++++-- .../InterpolationScheme/UpstreamScheme.f90 | 6 ++ src/Model/TransportModel/tsp-adv.f90 | 2 + 6 files changed, 118 insertions(+), 6 deletions(-) diff --git a/src/Model/TransportModel/InterpolationScheme/CentralDifferenceScheme.f90 b/src/Model/TransportModel/InterpolationScheme/CentralDifferenceScheme.f90 index 341a3f49e16..1ff19ff8166 100644 --- a/src/Model/TransportModel/InterpolationScheme/CentralDifferenceScheme.f90 +++ b/src/Model/TransportModel/InterpolationScheme/CentralDifferenceScheme.f90 @@ -17,6 +17,7 @@ module CentralDifferenceSchemeModule type(TspFmiType), pointer :: fmi contains procedure :: compute + procedure :: invalidate end type CentralDifferenceSchemeType interface CentralDifferenceSchemeType @@ -66,4 +67,10 @@ function compute(this, n, m, iposnm, phi) result(phi_face) phi_face%c_m = DONE - omega end function compute + + subroutine invalidate(this) + ! -- dummy + class(CentralDifferenceSchemeType), target :: this + end subroutine invalidate + end module CentralDifferenceSchemeModule diff --git a/src/Model/TransportModel/InterpolationScheme/InterpolationSchemeInterface.f90 b/src/Model/TransportModel/InterpolationScheme/InterpolationSchemeInterface.f90 index 0f4a853deb5..db7939d6502 100644 --- a/src/Model/TransportModel/InterpolationScheme/InterpolationSchemeInterface.f90 +++ b/src/Model/TransportModel/InterpolationScheme/InterpolationSchemeInterface.f90 @@ -16,6 +16,7 @@ module InterpolationSchemeInterfaceModule type, abstract :: InterpolationSchemeInterface contains procedure(compute_if), deferred :: compute + procedure(invalidate_if), deferred :: invalidate end type InterpolationSchemeInterface abstract interface @@ -35,6 +36,13 @@ function compute_if(this, n, m, iposnm, phi) result(phi_face) real(DP), intent(in), dimension(:) :: phi end function + subroutine invalidate_if(this) + ! -- import + import InterpolationSchemeInterface + ! -- dummy + class(InterpolationSchemeInterface), target :: this + end subroutine + end interface end module InterpolationSchemeInterfaceModule diff --git a/src/Model/TransportModel/InterpolationScheme/TVDScheme.f90 b/src/Model/TransportModel/InterpolationScheme/TVDScheme.f90 index 9a8bcfeb9c1..e51bb8bea57 100644 --- a/src/Model/TransportModel/InterpolationScheme/TVDScheme.f90 +++ b/src/Model/TransportModel/InterpolationScheme/TVDScheme.f90 @@ -18,6 +18,7 @@ module TVDSchemeModule integer(I4B), dimension(:), pointer, contiguous :: ibound => null() !< pointer to model ibound contains procedure :: compute + procedure :: invalidate end type TVDSchemeType interface TVDSchemeType @@ -109,4 +110,9 @@ function compute(this, n, m, iposnm, phi) result(phi_face) end function compute + subroutine invalidate(this) + ! -- dummy + class(TVDSchemeType), target :: this + end subroutine invalidate + end module TVDSchemeModule diff --git a/src/Model/TransportModel/InterpolationScheme/UTVDScheme.f90 b/src/Model/TransportModel/InterpolationScheme/UTVDScheme.f90 index 8c335df43c6..3289a18ad8d 100644 --- a/src/Model/TransportModel/InterpolationScheme/UTVDScheme.f90 +++ b/src/Model/TransportModel/InterpolationScheme/UTVDScheme.f90 @@ -38,11 +38,20 @@ module UTVDSchemeModule type(TspFmiType), pointer :: fmi class(IGradientType), pointer :: gradient integer(I4B) :: limiter_id = 2 ! default to van Leer limiter + logical :: cache_valid = .false. ! indicates if the cached gradients are valid + real(DP), dimension(:, :), allocatable :: cached_gradients + real(DP), dimension(:), allocatable :: cached_min_phi + real(DP), dimension(:), allocatable :: cached_max_phi + real(DP), dimension(:, :), allocatable :: cached_node_distance contains procedure :: compute + procedure :: invalidate procedure, private :: find_local_extrema procedure, private :: limiter + procedure, private :: compute_gradients + procedure, private :: compute_local_extrama + procedure, private :: compute_node_distance end type UTVDSchemeType interface UTVDSchemeType @@ -63,8 +72,69 @@ function constructor(dis, fmi, gradient) & interpolation_scheme%fmi => fmi interpolation_scheme%gradient => gradient + interpolation_scheme%cache_valid = .false. + allocate(interpolation_scheme%cached_gradients(dis%nodes, 3)) + allocate(interpolation_scheme%cached_min_phi(dis%nodes)) + allocate(interpolation_scheme%cached_max_phi(dis%nodes)) + + allocate(interpolation_scheme%cached_node_distance(dis%njas, 3)) + call compute_node_distance(interpolation_scheme) + end function constructor + subroutine invalidate(this) + ! -- dummy + class(UTVDSchemeType), target :: this + + this%cache_valid = .false. + end subroutine invalidate + + subroutine compute_gradients(this, phi) + ! -- dummy + class(UTVDSchemeType), target :: this + real(DP), intent(in), dimension(:) :: phi + ! -- local + integer(I4B) :: n + + do n = 1, this%dis%nodes + this%cached_gradients(n, :) = this%gradient%get(n, phi) + end do + end subroutine compute_gradients + + subroutine compute_local_extrama(this, phi) + ! -- dummy + class(UTVDSchemeType), target :: this + real(DP), intent(in), dimension(:) :: phi + ! -- local + integer(I4B) :: n + real(DP) :: min_phi, max_phi + + do n = 1, this%dis%nodes + call this%find_local_extrema(n, phi, min_phi, max_phi) + this%cached_min_phi(n) = min_phi + this%cached_max_phi(n) = max_phi + end do + end subroutine compute_local_extrama + + subroutine compute_node_distance(this) + ! -- dummy + class(UTVDSchemeType), target :: this + ! -- local + integer(I4B) :: n, m, ipos, isympos + + this%cached_node_distance = 0.0_dp + do n = 1, this%dis%nodes + do ipos = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1 + m = this%dis%con%ja(ipos) + if (m <= n) cycle + + isympos = this%dis%con%jas(ipos) + this%cached_node_distance(isympos, :) = node_distance(this%dis, n, m) + end do + end do + + end subroutine compute_node_distance + 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 @@ -78,14 +148,21 @@ function compute(this, n, m, iposnm, phi) result(phi_face) 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(:), pointer :: 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 + real(DP), pointer :: min_phi, max_phi ! Local minimum and maximum among cell and neighbors + + + if (.not. this%cache_valid) then + call this%compute_gradients(phi) + call this%compute_local_extrama(phi) + this%cache_valid = .true. + end if isympos = this%dis%con%jas(iposnm) qnm = this%fmi%gwfflowja(iposnm) @@ -107,10 +184,16 @@ function compute(this, n, m, iposnm, phi) result(phi_face) 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 + + if (iup > idn) then + dnm = -this%cached_node_distance(isympos, :) + else + dnm = this%cached_node_distance(isympos, :) + end if ! ! -- Add low order terms coef_up = DONE @@ -121,16 +204,16 @@ function compute(this, n, m, iposnm, phi) result(phi_face) if (abs(phi(idn) - phi(iup)) < DSAME) return ! ! -- Compute cell concentration gradient - grad_c = this%gradient%get(iup, phi) + grad_c => this%cached_gradients(iup, :) ! ! 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) + min_phi => this%cached_min_phi(iup) + max_phi => this%cached_max_phi(iup) if (c_virtual > max_phi) then c_virtual = max_phi diff --git a/src/Model/TransportModel/InterpolationScheme/UpstreamScheme.f90 b/src/Model/TransportModel/InterpolationScheme/UpstreamScheme.f90 index 1729d1c6a22..023ed4411bc 100644 --- a/src/Model/TransportModel/InterpolationScheme/UpstreamScheme.f90 +++ b/src/Model/TransportModel/InterpolationScheme/UpstreamScheme.f90 @@ -17,6 +17,7 @@ module UpstreamSchemeModule type(TspFmiType), pointer :: fmi contains procedure :: compute + procedure :: invalidate end type UpstreamSchemeType interface UpstreamSchemeType @@ -58,4 +59,9 @@ function compute(this, n, m, iposnm, phi) result(phi_face) end if end function compute + + subroutine invalidate(this) + ! -- dummy + class(UpstreamSchemeType), target :: this + end subroutine invalidate end module UpstreamSchemeModule diff --git a/src/Model/TransportModel/tsp-adv.f90 b/src/Model/TransportModel/tsp-adv.f90 index e20262e7fc0..205773cf133 100644 --- a/src/Model/TransportModel/tsp-adv.f90 +++ b/src/Model/TransportModel/tsp-adv.f90 @@ -231,6 +231,7 @@ subroutine adv_fc(this, nodes, matrix_sln, idxglo, cnew, rhs) type(CoefficientsType) :: coefficients ! Calculate internal domain fluxes and add to matrix_sln and rhs. + call this%face_interpolation%invalidate() do n = 1, nodes if (this%ibound(n) == 0) cycle ! skip inactive nodes idiag = this%dis%con%ia(n) @@ -267,6 +268,7 @@ subroutine adv_cq(this, cnew, flowja) ! rate and has dimensions of L^/T. nodes = this%dis%nodes + call this%face_interpolation%invalidate() do n = 1, nodes if (this%ibound(n) == 0) cycle do ipos = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1 From fa0bba042e5a3824aeee920de4ba4feffbaed4f0 Mon Sep 17 00:00:00 2001 From: Sunny Titus Date: Wed, 30 Jul 2025 16:18:06 +0200 Subject: [PATCH 11/15] Clean up --- .../CentralDifferenceScheme.f90 | 6 +++ .../InterpolationSchemeInterface.f90 | 2 +- .../InterpolationScheme/TVDScheme.f90 | 6 +++ .../InterpolationScheme/UTVDScheme.f90 | 52 +++++++++++++------ .../InterpolationScheme/UpstreamScheme.f90 | 7 +++ 5 files changed, 55 insertions(+), 18 deletions(-) diff --git a/src/Model/TransportModel/InterpolationScheme/CentralDifferenceScheme.f90 b/src/Model/TransportModel/InterpolationScheme/CentralDifferenceScheme.f90 index 1ff19ff8166..e4fccc05a53 100644 --- a/src/Model/TransportModel/InterpolationScheme/CentralDifferenceScheme.f90 +++ b/src/Model/TransportModel/InterpolationScheme/CentralDifferenceScheme.f90 @@ -68,9 +68,15 @@ function compute(this, n, m, iposnm, phi) result(phi_face) end function compute + !> @brief Invalidate cached data + !! + !! This method is required by the InterpolationSchemeInterface but does + !! nothing for this scheme since it does not cache anything. + !> subroutine invalidate(this) ! -- dummy class(CentralDifferenceSchemeType), target :: this + ! -- No cached data to invalidate end subroutine invalidate end module CentralDifferenceSchemeModule diff --git a/src/Model/TransportModel/InterpolationScheme/InterpolationSchemeInterface.f90 b/src/Model/TransportModel/InterpolationScheme/InterpolationSchemeInterface.f90 index db7939d6502..051e744a103 100644 --- a/src/Model/TransportModel/InterpolationScheme/InterpolationSchemeInterface.f90 +++ b/src/Model/TransportModel/InterpolationScheme/InterpolationSchemeInterface.f90 @@ -36,7 +36,7 @@ function compute_if(this, n, m, iposnm, phi) result(phi_face) real(DP), intent(in), dimension(:) :: phi end function - subroutine invalidate_if(this) + subroutine invalidate_if(this) ! -- import import InterpolationSchemeInterface ! -- dummy diff --git a/src/Model/TransportModel/InterpolationScheme/TVDScheme.f90 b/src/Model/TransportModel/InterpolationScheme/TVDScheme.f90 index e51bb8bea57..dc8e3c6b286 100644 --- a/src/Model/TransportModel/InterpolationScheme/TVDScheme.f90 +++ b/src/Model/TransportModel/InterpolationScheme/TVDScheme.f90 @@ -110,9 +110,15 @@ function compute(this, n, m, iposnm, phi) result(phi_face) end function compute + !> @brief Invalidate cached data + !! + !! This method is required by the InterpolationSchemeInterface but does + !! nothing for this scheme since it does not cache anything. + !> subroutine invalidate(this) ! -- dummy class(TVDSchemeType), target :: this + ! -- No cached data to invalidate end subroutine invalidate end module TVDSchemeModule diff --git a/src/Model/TransportModel/InterpolationScheme/UTVDScheme.f90 b/src/Model/TransportModel/InterpolationScheme/UTVDScheme.f90 index 3289a18ad8d..bf7b90fcdfd 100644 --- a/src/Model/TransportModel/InterpolationScheme/UTVDScheme.f90 +++ b/src/Model/TransportModel/InterpolationScheme/UTVDScheme.f90 @@ -39,18 +39,19 @@ module UTVDSchemeModule class(IGradientType), pointer :: gradient integer(I4B) :: limiter_id = 2 ! default to van Leer limiter logical :: cache_valid = .false. ! indicates if the cached gradients are valid - real(DP), dimension(:, :), allocatable :: cached_gradients - real(DP), dimension(:), allocatable :: cached_min_phi - real(DP), dimension(:), allocatable :: cached_max_phi - real(DP), dimension(:, :), allocatable :: cached_node_distance + real(DP), dimension(:, :), allocatable :: cached_gradients ! concentration gradients at nodes + real(DP), dimension(:), allocatable :: cached_min_phi ! minimum concentration among node and neighbors + real(DP), dimension(:), allocatable :: cached_max_phi ! maximum concentration among node and neighbors + real(DP), dimension(:, :), allocatable :: cached_node_distance ! distance vectors contains procedure :: compute procedure :: invalidate + procedure :: finalize procedure, private :: find_local_extrema procedure, private :: limiter procedure, private :: compute_gradients - procedure, private :: compute_local_extrama + procedure, private :: compute_local_extrema procedure, private :: compute_node_distance end type UTVDSchemeType @@ -63,7 +64,7 @@ function constructor(dis, fmi, gradient) & result(interpolation_scheme) ! -- return type(UTVDSchemeType) :: interpolation_scheme - ! --dummy + ! -- dummy class(DisBaseType), pointer, intent(in) :: dis type(TspFmiType), pointer, intent(in) :: fmi class(IGradientType), allocatable, target, intent(in) :: gradient @@ -73,15 +74,28 @@ function constructor(dis, fmi, gradient) & interpolation_scheme%gradient => gradient interpolation_scheme%cache_valid = .false. - allocate(interpolation_scheme%cached_gradients(dis%nodes, 3)) - allocate(interpolation_scheme%cached_min_phi(dis%nodes)) - allocate(interpolation_scheme%cached_max_phi(dis%nodes)) + allocate (interpolation_scheme%cached_gradients(dis%nodes, 3)) + allocate (interpolation_scheme%cached_min_phi(dis%nodes)) + allocate (interpolation_scheme%cached_max_phi(dis%nodes)) - allocate(interpolation_scheme%cached_node_distance(dis%njas, 3)) + allocate (interpolation_scheme%cached_node_distance(dis%njas, 3)) call compute_node_distance(interpolation_scheme) end function constructor + subroutine finalize(this) + ! -- dummy + class(UTVDSchemeType), intent(inout) :: this + + deallocate (this%cached_gradients) + deallocate (this%cached_min_phi) + deallocate (this%cached_max_phi) + deallocate (this%cached_node_distance) + end subroutine finalize + + !> @brief Invalidate cached data to force recomputation on next access + !! + !> subroutine invalidate(this) ! -- dummy class(UTVDSchemeType), target :: this @@ -97,11 +111,11 @@ subroutine compute_gradients(this, phi) integer(I4B) :: n do n = 1, this%dis%nodes - this%cached_gradients(n, :) = this%gradient%get(n, phi) + this%cached_gradients(n, :) = this%gradient%get(n, phi) end do end subroutine compute_gradients - subroutine compute_local_extrama(this, phi) + subroutine compute_local_extrema(this, phi) ! -- dummy class(UTVDSchemeType), target :: this real(DP), intent(in), dimension(:) :: phi @@ -114,7 +128,7 @@ subroutine compute_local_extrama(this, phi) this%cached_min_phi(n) = min_phi this%cached_max_phi(n) = max_phi end do - end subroutine compute_local_extrama + end subroutine compute_local_extrema subroutine compute_node_distance(this) ! -- dummy @@ -127,7 +141,7 @@ subroutine compute_node_distance(this) do ipos = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1 m = this%dis%con%ja(ipos) if (m <= n) cycle - + isympos = this%dis%con%jas(ipos) this%cached_node_distance(isympos, :) = node_distance(this%dis, n, m) end do @@ -157,10 +171,9 @@ function compute(this, n, m, iposnm, phi) result(phi_face) real(DP) :: c_virtual ! Virtual node concentration (Darwish method) real(DP), pointer :: min_phi, max_phi ! Local minimum and maximum among cell and neighbors - if (.not. this%cache_valid) then call this%compute_gradients(phi) - call this%compute_local_extrama(phi) + call this%compute_local_extrema(phi) this%cache_valid = .true. end if @@ -184,11 +197,16 @@ function compute(this, n, m, iposnm, phi) result(phi_face) 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 + ! Determine direction of distance vector from upwind to downwind cell + ! The cached_node_distance always stores vector from lower-numbered node to higher-numbered node. + ! Since we need dnm to point from upwind (iup) to downwind (idn), we must adjust the sign: + ! - If iup > idn: the cached vector points from idn to iup, so we negate it to get iup to idn + ! - If iup < idn: the cached vector already points from iup to idn, so use it as-is if (iup > idn) then dnm = -this%cached_node_distance(isympos, :) else diff --git a/src/Model/TransportModel/InterpolationScheme/UpstreamScheme.f90 b/src/Model/TransportModel/InterpolationScheme/UpstreamScheme.f90 index 023ed4411bc..05af7b140a0 100644 --- a/src/Model/TransportModel/InterpolationScheme/UpstreamScheme.f90 +++ b/src/Model/TransportModel/InterpolationScheme/UpstreamScheme.f90 @@ -60,8 +60,15 @@ function compute(this, n, m, iposnm, phi) result(phi_face) end function compute + !> @brief Invalidate cached data + !! + !! This method is required by the InterpolationSchemeInterface but does + !! nothing for this scheme since it does not cache anything. + !> subroutine invalidate(this) ! -- dummy class(UpstreamSchemeType), target :: this + ! -- No cached data to invalidate end subroutine invalidate + end module UpstreamSchemeModule From 85e752e74631c40d9157de8d449656ce7e224f53 Mon Sep 17 00:00:00 2001 From: Sunny Titus Date: Fri, 1 Aug 2025 11:23:45 +0200 Subject: [PATCH 12/15] Fix deconstructor call --- .../TransportModel/InterpolationScheme/UTVDScheme.f90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/Model/TransportModel/InterpolationScheme/UTVDScheme.f90 b/src/Model/TransportModel/InterpolationScheme/UTVDScheme.f90 index 3bd3c1901ce..42b9534dfb7 100644 --- a/src/Model/TransportModel/InterpolationScheme/UTVDScheme.f90 +++ b/src/Model/TransportModel/InterpolationScheme/UTVDScheme.f90 @@ -46,7 +46,7 @@ module UTVDSchemeModule contains procedure :: compute procedure :: invalidate - procedure :: finalize + final :: destructor procedure, private :: find_local_extrema procedure, private :: limiter @@ -83,15 +83,15 @@ function constructor(dis, fmi, gradient) & end function constructor - subroutine finalize(this) + subroutine destructor(this) ! -- dummy - class(UTVDSchemeType), intent(inout) :: this + type(UTVDSchemeType), intent(inout) :: this deallocate (this%cached_gradients) deallocate (this%cached_min_phi) deallocate (this%cached_max_phi) deallocate (this%cached_node_distance) - end subroutine finalize + end subroutine destructor !> @brief Invalidate cached data to force recomputation on next access !! From b3cd26cbefee69d9401d297dbad51f1255c67688 Mon Sep 17 00:00:00 2001 From: Sunny Titus Date: Fri, 1 Aug 2025 11:28:38 +0200 Subject: [PATCH 13/15] Move some function around --- .../InterpolationScheme/UTVDScheme.f90 | 92 +++++++++---------- 1 file changed, 46 insertions(+), 46 deletions(-) diff --git a/src/Model/TransportModel/InterpolationScheme/UTVDScheme.f90 b/src/Model/TransportModel/InterpolationScheme/UTVDScheme.f90 index 42b9534dfb7..c8195a64959 100644 --- a/src/Model/TransportModel/InterpolationScheme/UTVDScheme.f90 +++ b/src/Model/TransportModel/InterpolationScheme/UTVDScheme.f90 @@ -103,52 +103,6 @@ subroutine invalidate(this) this%cache_valid = .false. end subroutine invalidate - subroutine compute_gradients(this, phi) - ! -- dummy - class(UTVDSchemeType), target :: this - real(DP), intent(in), dimension(:) :: phi - ! -- local - integer(I4B) :: n - - do n = 1, this%dis%nodes - this%cached_gradients(n, :) = this%gradient%get(n, phi) - end do - end subroutine compute_gradients - - subroutine compute_local_extrema(this, phi) - ! -- dummy - class(UTVDSchemeType), target :: this - real(DP), intent(in), dimension(:) :: phi - ! -- local - integer(I4B) :: n - real(DP) :: min_phi, max_phi - - do n = 1, this%dis%nodes - call this%find_local_extrema(n, phi, min_phi, max_phi) - this%cached_min_phi(n) = min_phi - this%cached_max_phi(n) = max_phi - end do - end subroutine compute_local_extrema - - subroutine compute_node_distance(this) - ! -- dummy - class(UTVDSchemeType), target :: this - ! -- local - integer(I4B) :: n, m, ipos, isympos - - this%cached_node_distance = 0.0_dp - do n = 1, this%dis%nodes - do ipos = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1 - m = this%dis%con%ja(ipos) - if (m <= n) cycle - - isympos = this%dis%con%jas(ipos) - this%cached_node_distance(isympos, :) = node_distance(this%dis, n, m) - end do - end do - - end subroutine compute_node_distance - 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 @@ -303,4 +257,50 @@ function limiter(this, r) result(theta) end select end function + subroutine compute_gradients(this, phi) + ! -- dummy + class(UTVDSchemeType), target :: this + real(DP), intent(in), dimension(:) :: phi + ! -- local + integer(I4B) :: n + + do n = 1, this%dis%nodes + this%cached_gradients(n, :) = this%gradient%get(n, phi) + end do + end subroutine compute_gradients + + subroutine compute_local_extrema(this, phi) + ! -- dummy + class(UTVDSchemeType), target :: this + real(DP), intent(in), dimension(:) :: phi + ! -- local + integer(I4B) :: n + real(DP) :: min_phi, max_phi + + do n = 1, this%dis%nodes + call this%find_local_extrema(n, phi, min_phi, max_phi) + this%cached_min_phi(n) = min_phi + this%cached_max_phi(n) = max_phi + end do + end subroutine compute_local_extrema + + subroutine compute_node_distance(this) + ! -- dummy + class(UTVDSchemeType), target :: this + ! -- local + integer(I4B) :: n, m, ipos, isympos + + this%cached_node_distance = 0.0_dp + do n = 1, this%dis%nodes + do ipos = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1 + m = this%dis%con%ja(ipos) + if (m <= n) cycle + + isympos = this%dis%con%jas(ipos) + this%cached_node_distance(isympos, :) = node_distance(this%dis, n, m) + end do + end do + + end subroutine compute_node_distance + end module UTVDSchemeModule From 8f283007baba01483500d4f8589214b2c77f7a03 Mon Sep 17 00:00:00 2001 From: Sunny Titus Date: Thu, 7 Aug 2025 13:34:53 +0200 Subject: [PATCH 14/15] Rename variable to better match documentation nomenclature --- .../TransportModel/Gradient/LeastSquaresGradient.f90 | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/Model/TransportModel/Gradient/LeastSquaresGradient.f90 b/src/Model/TransportModel/Gradient/LeastSquaresGradient.f90 index 5e78ec0cd41..65521c84346 100644 --- a/src/Model/TransportModel/Gradient/LeastSquaresGradient.f90 +++ b/src/Model/TransportModel/Gradient/LeastSquaresGradient.f90 @@ -79,16 +79,16 @@ function create_gradient_reconstruction_matrix(this, n) result(R) 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), + real(DP), dimension(:, :), allocatable :: inverse_distance ! 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)) + allocate (inverse_distance(number_connections, number_connections)) - grad_scale = 0 + inverse_distance = 0 d = 0 ! Assemble the distance matrix @@ -101,13 +101,13 @@ function create_gradient_reconstruction_matrix(this, n) result(R) length = norm2(dnm) d(local_pos, :) = dnm / length - grad_scale(local_pos, local_pos) = 1.0_dp / length + inverse_distance(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) + R = matmul(pinv(d), inverse_distance) end function create_gradient_reconstruction_matrix From c2c0b1ab3426d36556ed364d035fd67a0853d3d0 Mon Sep 17 00:00:00 2001 From: Sunny Titus Date: Thu, 14 Aug 2025 21:27:23 +0200 Subject: [PATCH 15/15] Apply review comments. Make the flow more explicit --- make/makefile | 2 + msvs/mf6core.vfproj | 2 + .../Gradient/CachedGradient.f90 | 95 +++++++++++++ .../TransportModel/Gradient/IGradient.f90 | 13 +- .../Gradient/LeastSquaresGradient.f90 | 28 +++- .../CentralDifferenceScheme.f90 | 30 +++-- .../InterpolationSchemeInterface.f90 | 21 ++- .../InterpolationScheme/TVDScheme.f90 | 38 +++--- .../InterpolationScheme/UTVDScheme.f90 | 114 +++++----------- .../InterpolationScheme/UpstreamScheme.f90 | 30 +++-- src/Model/TransportModel/tsp-adv.f90 | 18 +-- src/Utilities/LocalCellExtrema.f90 | 127 ++++++++++++++++++ src/meson.build | 2 + 13 files changed, 370 insertions(+), 150 deletions(-) create mode 100644 src/Model/TransportModel/Gradient/CachedGradient.f90 create mode 100644 src/Utilities/LocalCellExtrema.f90 diff --git a/make/makefile b/make/makefile index 84253ddd2bc..985ba40037b 100644 --- a/make/makefile +++ b/make/makefile @@ -341,6 +341,7 @@ $(OBJDIR)/VectorInterpolation.o \ $(OBJDIR)/swf-cxs.o \ $(OBJDIR)/OutputControlData.o \ $(OBJDIR)/gwf-ic.o \ +$(OBJDIR)/LocalCellExtrema.o \ $(OBJDIR)/InterpolationSchemeInterface.o \ $(OBJDIR)/IGradient.o \ $(OBJDIR)/DisUtils.o \ @@ -372,6 +373,7 @@ $(OBJDIR)/TVDScheme.o \ $(OBJDIR)/TspAdvOptions.o \ $(OBJDIR)/LeastSquaresGradient.o \ $(OBJDIR)/CentralDifferenceScheme.o \ +$(OBJDIR)/CachedGradient.o \ $(OBJDIR)/AdvSchemeEnum.o \ $(OBJDIR)/UzfCellGroup.o \ $(OBJDIR)/Xt3dInterface.o \ diff --git a/msvs/mf6core.vfproj b/msvs/mf6core.vfproj index 66cff1389dd..a418134c97d 100644 --- a/msvs/mf6core.vfproj +++ b/msvs/mf6core.vfproj @@ -384,6 +384,7 @@ + @@ -685,6 +686,7 @@ + diff --git a/src/Model/TransportModel/Gradient/CachedGradient.f90 b/src/Model/TransportModel/Gradient/CachedGradient.f90 new file mode 100644 index 00000000000..244abf2f4bc --- /dev/null +++ b/src/Model/TransportModel/Gradient/CachedGradient.f90 @@ -0,0 +1,95 @@ +module CachedGradientModule + use KindModule, only: DP, I4B + use ConstantsModule, only: DONE + + Use IGradient + use BaseDisModule, only: DisBaseType + + implicit none + private + + public :: CachedGradientType + + !> @brief Decorator that adds caching to any gradient computation implementation + !! + !! This class wraps any IGradientType implementation and provides caching functionality + !! using the Decorator pattern. When set_field is called, it pre-computes gradients + !! for all nodes and stores them in memory for fast O(1) retrieval. This trades memory + !! for speed when gradients are accessed multiple times for the same scalar field. + !! + !! The class takes ownership of the wrapped gradient object via move semantics and + !! provides the same interface as any other gradient implementation. This allows it + !! to be used transparently in place of the original gradient object. + !! + !! Usage pattern: + !! 1. Create a base gradient implementation (e.g., LeastSquaresGradientType) + !! 2. Wrap it with CachedGradientType using the constructor + !! 3. Call set_field() once to pre-compute all gradients + !! 4. Call get() multiple times for fast cached lookups + !! + !! @note The wrapped gradient object is moved (not copied) during construction + !! for efficient memory management. + !< + type, extends(IGradientType) :: CachedGradientType + private + class(DisBaseType), pointer :: dis + class(IGradientType), allocatable :: gradient + real(DP), dimension(:, :), allocatable :: cached_gradients ! gradients at nodes + contains + procedure :: get + procedure :: set_field + final :: destructor + end type CachedGradientType + + interface CachedGradientType + module procedure Constructor + end interface CachedGradientType + +contains + + function constructor(gradient, dis) result(cached_gradient) + ! --dummy + class(IGradientType), allocatable, intent(inout) :: gradient + class(DisBaseType), pointer, intent(in) :: dis + !-- return + type(CachedGradientType) :: cached_gradient + + cached_gradient%dis => dis + + call move_alloc(gradient, cached_gradient%gradient) ! Take ownership + allocate (cached_gradient%cached_gradients(dis%nodes, 3)) + + end function constructor + + subroutine destructor(this) + ! -- dummy + type(CachedGradientType), intent(inout) :: this + + deallocate (this%cached_gradients) + end subroutine destructor + + function get(this, n) result(grad_c) + ! -- dummy + class(CachedGradientType), target :: this + integer(I4B), intent(in) :: n + !-- return + real(DP), dimension(3) :: grad_c + + grad_c = this%cached_gradients(n, :) + end function get + + subroutine set_field(this, phi) + ! -- dummy + class(CachedGradientType), target :: this + real(DP), dimension(:), pointer, intent(in) :: phi + ! -- local + integer(I4B) :: n + + call this%gradient%set_field(phi) + do n = 1, this%dis%nodes + this%cached_gradients(n, :) = this%gradient%get(n) + end do + + end subroutine set_field + +end module CachedGradientModule diff --git a/src/Model/TransportModel/Gradient/IGradient.f90 b/src/Model/TransportModel/Gradient/IGradient.f90 index 6a3268a4485..6ef20137772 100644 --- a/src/Model/TransportModel/Gradient/IGradient.f90 +++ b/src/Model/TransportModel/Gradient/IGradient.f90 @@ -15,22 +15,31 @@ module IGradient type, abstract :: IGradientType contains procedure(get_if), deferred :: get + procedure(set_field_if), deferred :: set_field end type IGradientType abstract interface - function get_if(this, n, c) result(grad_c) + function get_if(this, n) 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 + subroutine set_field_if(this, phi) + ! -- import + import IGradientType + import DP, I4B + ! -- dummy + class(IGradientType), target :: this + real(DP), dimension(:), pointer, intent(in) :: phi + end subroutine + end interface contains diff --git a/src/Model/TransportModel/Gradient/LeastSquaresGradient.f90 b/src/Model/TransportModel/Gradient/LeastSquaresGradient.f90 index 65521c84346..d89e1beb270 100644 --- a/src/Model/TransportModel/Gradient/LeastSquaresGradient.f90 +++ b/src/Model/TransportModel/Gradient/LeastSquaresGradient.f90 @@ -25,20 +25,28 @@ module LeastSquaresGradientModule !! The gradient can then be computed by multiplying the reconstruction matrix with the difference vector. !! ∇ɸ = R * ∑(ɸ_i - ɸ_up), where i are the neighboring cells. !! + !! Usage: + !! 1. Create the gradient object with the discretization + !! 2. Set the scalar field using `set_field(phi)` where phi is the field for which gradients are computed + !! 3. Retrieve gradients for any cell using the `get(n)` method + !! !! - 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. + !! - The `set_field` method establishes a pointer to the scalar field data. + !! - The `get` method computes the gradient for any cell using the current 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 + real(DP), dimension(:), pointer :: phi type(Array2D), allocatable, dimension(:) :: R ! Gradient reconstruction matrix contains procedure :: get + procedure :: set_field procedure, private :: compute_cell_gradient procedure, private :: create_gradient_reconstruction_matrix @@ -111,24 +119,30 @@ function create_gradient_reconstruction_matrix(this, n) result(R) end function create_gradient_reconstruction_matrix - function get(this, n, c) result(grad_c) + function get(this, n) 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) + grad_c = this%compute_cell_gradient(n) end function get - function compute_cell_gradient(this, n, phi_new) result(grad_c) + subroutine set_field(this, phi) + ! -- dummy + class(LeastSquaresGradientType), target :: this + real(DP), dimension(:), pointer, intent(in) :: phi + + this%phi => phi + end subroutine set_field + + function compute_cell_gradient(this, n) 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 @@ -143,7 +157,7 @@ function compute_cell_gradient(this, n, phi_new) result(grad_c) 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) + dc(local_pos) = this%phi(m) - this%phi(n) local_pos = local_pos + 1 end do diff --git a/src/Model/TransportModel/InterpolationScheme/CentralDifferenceScheme.f90 b/src/Model/TransportModel/InterpolationScheme/CentralDifferenceScheme.f90 index e4fccc05a53..a16246b91ec 100644 --- a/src/Model/TransportModel/InterpolationScheme/CentralDifferenceScheme.f90 +++ b/src/Model/TransportModel/InterpolationScheme/CentralDifferenceScheme.f90 @@ -15,9 +15,10 @@ module CentralDifferenceSchemeModule private class(DisBaseType), pointer :: dis type(TspFmiType), pointer :: fmi + real(DP), dimension(:), pointer :: phi contains procedure :: compute - procedure :: invalidate + procedure :: set_field end type CentralDifferenceSchemeType interface CentralDifferenceSchemeType @@ -37,7 +38,20 @@ function constructor(dis, fmi) result(interpolation_scheme) end function constructor - function compute(this, n, m, iposnm, phi) result(phi_face) + !> @brief Set the scalar field for which interpolation will be computed + !! + !! This method establishes a pointer to the scalar field data for + !! subsequent central difference interpolation computations. + !< + subroutine set_field(this, phi) + ! -- dummy + class(CentralDifferenceSchemeType), target :: this + real(DP), intent(in), dimension(:), pointer :: phi + + this%phi => phi + end subroutine set_field + + function compute(this, n, m, iposnm) result(phi_face) !-- return type(CoefficientsType) :: phi_face ! -- dummy @@ -45,7 +59,6 @@ function compute(this, n, m, iposnm, phi) result(phi_face) integer(I4B), intent(in) :: n integer(I4B), intent(in) :: m integer(I4B), intent(in) :: iposnm - real(DP), intent(in), dimension(:) :: phi ! -- local real(DP) :: lnm, lmn, omega @@ -68,15 +81,4 @@ function compute(this, n, m, iposnm, phi) result(phi_face) end function compute - !> @brief Invalidate cached data - !! - !! This method is required by the InterpolationSchemeInterface but does - !! nothing for this scheme since it does not cache anything. - !> - subroutine invalidate(this) - ! -- dummy - class(CentralDifferenceSchemeType), target :: this - ! -- No cached data to invalidate - end subroutine invalidate - end module CentralDifferenceSchemeModule diff --git a/src/Model/TransportModel/InterpolationScheme/InterpolationSchemeInterface.f90 b/src/Model/TransportModel/InterpolationScheme/InterpolationSchemeInterface.f90 index 051e744a103..320132b8115 100644 --- a/src/Model/TransportModel/InterpolationScheme/InterpolationSchemeInterface.f90 +++ b/src/Model/TransportModel/InterpolationScheme/InterpolationSchemeInterface.f90 @@ -16,12 +16,18 @@ module InterpolationSchemeInterfaceModule type, abstract :: InterpolationSchemeInterface contains procedure(compute_if), deferred :: compute - procedure(invalidate_if), deferred :: invalidate + procedure(set_field_if), deferred :: set_field end type InterpolationSchemeInterface abstract interface - function compute_if(this, n, m, iposnm, phi) result(phi_face) + !> @brief Compute interpolation coefficients for a face value + !! + !! This method computes the coefficients needed to interpolate a scalar field + !! value at the face between two connected cells. The method returns coefficients + !! that define how the face value depends on the cell-centered values. + !< + function compute_if(this, n, m, iposnm) result(phi_face) ! -- import import DP, I4B import InterpolationSchemeInterface @@ -33,14 +39,21 @@ function compute_if(this, n, m, iposnm, phi) result(phi_face) integer(I4B), intent(in) :: n integer(I4B), intent(in) :: m integer(I4B), intent(in) :: iposnm - real(DP), intent(in), dimension(:) :: phi end function - subroutine invalidate_if(this) + !> @brief Set the scalar field for which interpolation will be computed + !! + !! This method establishes a pointer to the scalar field data that will be + !! used for subsequent interpolation computations. Implementations may also + !! update any dependent cached data to ensure consistent results. + !< + subroutine set_field_if(this, phi) ! -- import + import DP import InterpolationSchemeInterface ! -- dummy class(InterpolationSchemeInterface), target :: this + real(DP), intent(in), dimension(:), pointer :: phi end subroutine end interface diff --git a/src/Model/TransportModel/InterpolationScheme/TVDScheme.f90 b/src/Model/TransportModel/InterpolationScheme/TVDScheme.f90 index dc8e3c6b286..5cd5ed2707c 100644 --- a/src/Model/TransportModel/InterpolationScheme/TVDScheme.f90 +++ b/src/Model/TransportModel/InterpolationScheme/TVDScheme.f90 @@ -15,10 +15,11 @@ module TVDSchemeModule private class(DisBaseType), pointer :: dis type(TspFmiType), pointer :: fmi + real(DP), dimension(:), pointer :: phi integer(I4B), dimension(:), pointer, contiguous :: ibound => null() !< pointer to model ibound contains procedure :: compute - procedure :: invalidate + procedure :: set_field end type TVDSchemeType interface TVDSchemeType @@ -40,7 +41,20 @@ function constructor(dis, fmi, ibound) result(interpolation_scheme) end function constructor - function compute(this, n, m, iposnm, phi) result(phi_face) + !> @brief Set the scalar field for which interpolation will be computed + !! + !! This method establishes a pointer to the scalar field data for + !! subsequent TVD interpolation computations. + !< + subroutine set_field(this, phi) + ! -- dummy + class(TVDSchemeType), target :: this + real(DP), intent(in), dimension(:), pointer :: phi + + this%phi => phi + end subroutine set_field + + function compute(this, n, m, iposnm) result(phi_face) !-- return type(CoefficientsType), target :: phi_face ! -- dummy @@ -48,7 +62,6 @@ function compute(this, n, m, iposnm, phi) result(phi_face) integer(I4B), intent(in) :: n integer(I4B), intent(in) :: m integer(I4B), intent(in) :: iposnm - real(DP), intent(in), dimension(:) :: phi ! -- local integer(I4B) :: ipos, isympos, iup, idn, i2up, j real(DP) :: qnm, qmax, qupj, elupdn, elup2up @@ -97,28 +110,17 @@ function compute(this, n, m, iposnm, phi) result(phi_face) ! -- Calculate flux limiting term if (i2up > 0) then smooth = DZERO - cdiff = ABS(phi(idn) - phi(iup)) + cdiff = ABS(this%phi(idn) - this%phi(iup)) if (cdiff > DPREC) then - smooth = (phi(iup) - phi(i2up)) / elup2up * & - elupdn / (phi(idn) - phi(iup)) + smooth = (this%phi(iup) - this%phi(i2up)) / elup2up * & + elupdn / (this%phi(idn) - this%phi(iup)) end if if (smooth > DZERO) then alimiter = DTWO * smooth / (DONE + smooth) - phi_face%rhs = -DHALF * alimiter * (phi(idn) - phi(iup)) + phi_face%rhs = -DHALF * alimiter * (this%phi(idn) - this%phi(iup)) end if end if end function compute - !> @brief Invalidate cached data - !! - !! This method is required by the InterpolationSchemeInterface but does - !! nothing for this scheme since it does not cache anything. - !> - subroutine invalidate(this) - ! -- dummy - class(TVDSchemeType), target :: this - ! -- No cached data to invalidate - end subroutine invalidate - end module TVDSchemeModule diff --git a/src/Model/TransportModel/InterpolationScheme/UTVDScheme.f90 b/src/Model/TransportModel/InterpolationScheme/UTVDScheme.f90 index c8195a64959..ed60f983139 100644 --- a/src/Model/TransportModel/InterpolationScheme/UTVDScheme.f90 +++ b/src/Model/TransportModel/InterpolationScheme/UTVDScheme.f90 @@ -6,7 +6,9 @@ module UTVDSchemeModule use BaseDisModule, only: DisBaseType use TspFmiModule, only: TspFmiType use IGradient, only: IGradientType + use DisUtilsModule, only: node_distance + use LocalCellExtremaModule, only: LocalCellExtremaType implicit none private @@ -37,21 +39,17 @@ module UTVDSchemeModule class(DisBaseType), pointer :: dis type(TspFmiType), pointer :: fmi class(IGradientType), pointer :: gradient + + real(DP), dimension(:), pointer :: phi + type(LocalCellExtremaType), allocatable :: min_max_phi ! local minimum values at nodes integer(I4B) :: limiter_id = 2 ! default to van Leer limiter - logical :: cache_valid = .false. ! indicates if the cached gradients are valid - real(DP), dimension(:, :), allocatable :: cached_gradients ! concentration gradients at nodes - real(DP), dimension(:), allocatable :: cached_min_phi ! minimum concentration among node and neighbors - real(DP), dimension(:), allocatable :: cached_max_phi ! maximum concentration among node and neighbors real(DP), dimension(:, :), allocatable :: cached_node_distance ! distance vectors contains procedure :: compute - procedure :: invalidate + procedure :: set_field final :: destructor - procedure, private :: find_local_extrema procedure, private :: limiter - procedure, private :: compute_gradients - procedure, private :: compute_local_extrema procedure, private :: compute_node_distance end type UTVDSchemeType @@ -60,6 +58,7 @@ module UTVDSchemeModule end interface UTVDSchemeType contains + function constructor(dis, fmi, gradient) & result(interpolation_scheme) ! -- return @@ -67,16 +66,13 @@ function constructor(dis, fmi, gradient) & ! -- dummy class(DisBaseType), pointer, intent(in) :: dis type(TspFmiType), pointer, intent(in) :: fmi - class(IGradientType), allocatable, target, intent(in) :: gradient + class(IGradientType), allocatable, intent(in), target :: gradient interpolation_scheme%dis => dis interpolation_scheme%fmi => fmi - interpolation_scheme%gradient => gradient - interpolation_scheme%cache_valid = .false. - allocate (interpolation_scheme%cached_gradients(dis%nodes, 3)) - allocate (interpolation_scheme%cached_min_phi(dis%nodes)) - allocate (interpolation_scheme%cached_max_phi(dis%nodes)) + interpolation_scheme%gradient => gradient + interpolation_scheme%min_max_phi = LocalCellExtremaType(dis) allocate (interpolation_scheme%cached_node_distance(dis%njas, 3)) call compute_node_distance(interpolation_scheme) @@ -87,23 +83,27 @@ subroutine destructor(this) ! -- dummy type(UTVDSchemeType), intent(inout) :: this - deallocate (this%cached_gradients) - deallocate (this%cached_min_phi) - deallocate (this%cached_max_phi) deallocate (this%cached_node_distance) end subroutine destructor - !> @brief Invalidate cached data to force recomputation on next access + !> @brief Set the scalar field for which interpolation will be computed + !! + !! This method establishes a pointer to the scalar field data and updates + !! any dependent cached data (gradients and local extrema) to ensure + !! subsequent interpolation computations use the current field values. !! - !> - subroutine invalidate(this) + !< + subroutine set_field(this, phi) ! -- dummy class(UTVDSchemeType), target :: this + real(DP), intent(in), dimension(:), pointer :: phi - this%cache_valid = .false. - end subroutine invalidate + this%phi => phi + call this%gradient%set_field(phi) + call this%min_max_phi%set_field(phi) + end subroutine set_field - function compute(this, n, m, iposnm, phi) result(phi_face) + function compute(this, n, m, iposnm) result(phi_face) !-- return type(CoefficientsType), target :: phi_face ! Output: coefficients for the face between cells n and m ! -- dummy @@ -111,12 +111,11 @@ function compute(this, n, m, iposnm, phi) result(phi_face) 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(:), pointer :: grad_c ! gradient at upwind cell + 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 @@ -125,12 +124,6 @@ function compute(this, n, m, iposnm, phi) result(phi_face) real(DP) :: c_virtual ! Virtual node concentration (Darwish method) real(DP), pointer :: min_phi, max_phi ! Local minimum and maximum among cell and neighbors - if (.not. this%cache_valid) then - call this%compute_gradients(phi) - call this%compute_local_extrema(phi) - this%cache_valid = .true. - end if - isympos = this%dis%con%jas(iposnm) qnm = this%fmi%gwfflowja(iposnm) ! @@ -173,19 +166,19 @@ function compute(this, n, m, iposnm, phi) result(phi_face) ! -- Add high order terms ! ! -- Return if straddled cells have same value - if (abs(phi(idn) - phi(iup)) < DSAME) return + if (abs(this%phi(idn) - this%phi(iup)) < DSAME) return ! ! -- Compute cell concentration gradient - grad_c => this%cached_gradients(iup, :) + grad_c = this%gradient%get(iup) ! ! Darwish's method to compute virtual node concentration - c_virtual = phi(idn) - 2.0_dp * (dot_product(grad_c, dnm)) + c_virtual = this%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. - min_phi => this%cached_min_phi(iup) - max_phi => this%cached_max_phi(iup) + min_phi => this%min_max_phi%get_min(iup) + max_phi => this%min_max_phi%get_max(iup) if (c_virtual > max_phi) then c_virtual = max_phi @@ -196,14 +189,14 @@ function compute(this, n, m, iposnm, phi) result(phi_face) end if ! ! -- Compute smoothness factor - smooth = (phi(iup) - c_virtual) / (phi(idn) - phi(iup)) + smooth = (this%phi(iup) - c_virtual) / (this%phi(idn) - this%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)) + phi_face%rhs = -relative_distance * alimiter * (this%phi(idn) - this%phi(iup)) ! Alternative way of writing the high order term by adding it to the ! coefficients matrix. The equation to be added is: @@ -214,24 +207,6 @@ function compute(this, n, m, iposnm, phi) result(phi_face) 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 @@ -257,33 +232,6 @@ function limiter(this, r) result(theta) end select end function - subroutine compute_gradients(this, phi) - ! -- dummy - class(UTVDSchemeType), target :: this - real(DP), intent(in), dimension(:) :: phi - ! -- local - integer(I4B) :: n - - do n = 1, this%dis%nodes - this%cached_gradients(n, :) = this%gradient%get(n, phi) - end do - end subroutine compute_gradients - - subroutine compute_local_extrema(this, phi) - ! -- dummy - class(UTVDSchemeType), target :: this - real(DP), intent(in), dimension(:) :: phi - ! -- local - integer(I4B) :: n - real(DP) :: min_phi, max_phi - - do n = 1, this%dis%nodes - call this%find_local_extrema(n, phi, min_phi, max_phi) - this%cached_min_phi(n) = min_phi - this%cached_max_phi(n) = max_phi - end do - end subroutine compute_local_extrema - subroutine compute_node_distance(this) ! -- dummy class(UTVDSchemeType), target :: this diff --git a/src/Model/TransportModel/InterpolationScheme/UpstreamScheme.f90 b/src/Model/TransportModel/InterpolationScheme/UpstreamScheme.f90 index 05af7b140a0..a12728e15e4 100644 --- a/src/Model/TransportModel/InterpolationScheme/UpstreamScheme.f90 +++ b/src/Model/TransportModel/InterpolationScheme/UpstreamScheme.f90 @@ -15,9 +15,10 @@ module UpstreamSchemeModule private class(DisBaseType), pointer :: dis type(TspFmiType), pointer :: fmi + real(DP), dimension(:), pointer :: phi contains procedure :: compute - procedure :: invalidate + procedure :: set_field end type UpstreamSchemeType interface UpstreamSchemeType @@ -37,7 +38,20 @@ function constructor(dis, fmi) result(interpolation_scheme) end function constructor - function compute(this, n, m, iposnm, phi) result(phi_face) + !> @brief Set the scalar field for which interpolation will be computed + !! + !! This method establishes a pointer to the scalar field data for + !! subsequent upstream interpolation computations. + !< + subroutine set_field(this, phi) + ! -- dummy + class(UpstreamSchemeType), target :: this + real(DP), intent(in), dimension(:), pointer :: phi + + this%phi => phi + end subroutine set_field + + function compute(this, n, m, iposnm) result(phi_face) !-- return type(CoefficientsType) :: phi_face ! -- dummy @@ -45,7 +59,6 @@ function compute(this, n, m, iposnm, phi) result(phi_face) integer(I4B), intent(in) :: n integer(I4B), intent(in) :: m integer(I4B), intent(in) :: iposnm - real(DP), intent(in), dimension(:) :: phi ! -- local real(DP) :: qnm @@ -60,15 +73,4 @@ function compute(this, n, m, iposnm, phi) result(phi_face) end function compute - !> @brief Invalidate cached data - !! - !! This method is required by the InterpolationSchemeInterface but does - !! nothing for this scheme since it does not cache anything. - !> - subroutine invalidate(this) - ! -- dummy - class(UpstreamSchemeType), target :: this - ! -- No cached data to invalidate - end subroutine invalidate - end module UpstreamSchemeModule diff --git a/src/Model/TransportModel/tsp-adv.f90 b/src/Model/TransportModel/tsp-adv.f90 index 8d91f8c1c31..43f6ce9307a 100644 --- a/src/Model/TransportModel/tsp-adv.f90 +++ b/src/Model/TransportModel/tsp-adv.f90 @@ -10,6 +10,7 @@ module TspAdvModule ! -- Gradient schemes use IGradient, only: IGradientType use LeastSquaresGradientModule, only: LeastSquaresGradientType + use CachedGradientModule, only: CachedGradientType ! -- Interpolation schemes use InterpolationSchemeInterfaceModule, only: InterpolationSchemeInterface, & CoefficientsType @@ -121,7 +122,7 @@ subroutine adv_ar(this, dis, ibound) integer(I4B), dimension(:), pointer, contiguous, intent(in) :: ibound ! -- local integer(I4B) :: iadvwt_value - ! + class(IGradientType), allocatable :: gradient ! -- adv pointers to arguments that were passed in this%dis => dis this%ibound => ibound @@ -139,7 +140,8 @@ subroutine adv_ar(this, dis, ibound) this%face_interpolation = & TVDSchemeType(this%dis, this%fmi, this%ibound) case (ADV_SCHEME_UTVD) - this%gradient = LeastSquaresGradientType(this%dis) + gradient = LeastSquaresGradientType(this%dis) + this%gradient = CachedGradientType(gradient, this%dis) this%face_interpolation = & UTVDSchemeType(this%dis, this%fmi, this%gradient) case default @@ -223,7 +225,7 @@ subroutine adv_fc(this, nodes, matrix_sln, idxglo, cnew, rhs) integer(I4B), intent(in) :: nodes !< number of nodes class(MatrixBaseType), pointer :: matrix_sln !< pointer to solution matrix integer(I4B), intent(in), dimension(:) :: idxglo !< global indices for matrix - real(DP), intent(in), dimension(:) :: cnew !< new concentration/temperature values + real(DP), intent(in), dimension(:), target :: cnew !< new concentration/temperature values real(DP), dimension(:), intent(inout) :: rhs !< right-hand side vector ! -- local integer(I4B) :: n, m, idiag, ipos @@ -231,7 +233,7 @@ subroutine adv_fc(this, nodes, matrix_sln, idxglo, cnew, rhs) type(CoefficientsType) :: coefficients ! Calculate internal domain fluxes and add to matrix_sln and rhs. - call this%face_interpolation%invalidate() + call this%face_interpolation%set_field(cnew) do n = 1, nodes if (this%ibound(n) == 0) cycle ! skip inactive nodes idiag = this%dis%con%ia(n) @@ -241,7 +243,7 @@ subroutine adv_fc(this, nodes, matrix_sln, idxglo, cnew, rhs) if (this%dis%con%mask(ipos) == 0) cycle ! skip masked connections qnm = this%fmi%gwfflowja(ipos) * this%eqnsclfac - coefficients = this%face_interpolation%compute(n, m, ipos, cnew) + coefficients = this%face_interpolation%compute(n, m, ipos) call matrix_sln%add_value_pos(idxglo(idiag), qnm * coefficients%c_n) call matrix_sln%add_value_pos(idxglo(ipos), qnm * coefficients%c_m) @@ -256,7 +258,7 @@ subroutine adv_cq(this, cnew, flowja) ! -- modules ! -- dummy class(TspAdvType) :: this - real(DP), intent(in), dimension(:) :: cnew + real(DP), intent(in), dimension(:), target :: cnew real(DP), intent(inout), dimension(:) :: flowja ! -- local integer(I4B) :: nodes @@ -268,7 +270,7 @@ subroutine adv_cq(this, cnew, flowja) ! rate and has dimensions of L^/T. nodes = this%dis%nodes - call this%face_interpolation%invalidate() + call this%face_interpolation%set_field(cnew) do n = 1, nodes if (this%ibound(n) == 0) cycle do ipos = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1 @@ -276,7 +278,7 @@ subroutine adv_cq(this, cnew, flowja) if (this%ibound(m) == 0) cycle qnm = this%fmi%gwfflowja(ipos) * this%eqnsclfac - coefficients = this%face_interpolation%compute(n, m, ipos, cnew) + coefficients = this%face_interpolation%compute(n, m, ipos) flowja(ipos) = flowja(ipos) & + qnm * coefficients%c_n * cnew(n) & + qnm * coefficients%c_m * cnew(m) & diff --git a/src/Utilities/LocalCellExtrema.f90 b/src/Utilities/LocalCellExtrema.f90 new file mode 100644 index 00000000000..cd46c7d1e2c --- /dev/null +++ b/src/Utilities/LocalCellExtrema.f90 @@ -0,0 +1,127 @@ +module LocalCellExtremaModule + use KindModule, only: DP, I4B + use BaseDisModule, only: DisBaseType + + implicit none + private + + public :: LocalCellExtremaType + + !> @brief Computes and caches local extrema for each cell and its connected neighbors + !! + !! This class computes the minimum and maximum values within the local + !! stencil of each cell (the cell itself plus all directly connected neighboring cells). + !! The extrema are computed once when the scalar field is set and then cached for + !! fast retrieval. This is particularly useful for TVD limiters and slope + !! limiting algorithms that need to enforce monotonicity constraints. + !! + !! The local extrema computation follows the connectivity pattern defined by the + !! discretization object, examining all cells that share a face with the target cell. + !! This creates a computational stencil that includes the central cell and its + !! immediate neighbors. + !! + !! @note The get_min() and get_max() methods return pointers for zero-copy performance + !! in performance-critical loops. Callers should treat returned values as read-only. + !< + type :: LocalCellExtremaType + private + class(DisBaseType), pointer :: dis + real(DP), dimension(:), allocatable :: min + real(DP), dimension(:), allocatable :: max + contains + procedure :: set_field + procedure :: get_min + procedure :: get_max + final :: destructor + + procedure, private :: compute_local_extrema + procedure, private :: find_local_extrema + end type LocalCellExtremaType + + interface LocalCellExtremaType + module procedure constructor + end interface LocalCellExtremaType + +contains + function constructor(dis) result(local_extrema) + ! -- return + type(LocalCellExtremaType) :: local_extrema + ! -- dummy + class(DisBaseType), pointer, intent(in) :: dis + + local_extrema%dis => dis + + allocate (local_extrema%min(dis%nodes)) + allocate (local_extrema%max(dis%nodes)) + end function constructor + + subroutine destructor(this) + ! -- dummy + type(LocalCellExtremaType), intent(inout) :: this + + deallocate (this%min) + deallocate (this%max) + end subroutine destructor + + subroutine set_field(this, phi) + ! -- dummy + class(LocalCellExtremaType), target :: this + real(DP), intent(in), dimension(:), pointer :: phi + + call this%compute_local_extrema(phi) + end subroutine set_field + + function get_min(this, n) result(min_val) + ! -- dummy + class(LocalCellExtremaType), target :: this + integer(I4B), intent(in) :: n ! Node index + !-- return + real(DP), pointer :: min_val + + min_val => this%min(n) + end function get_min + + function get_max(this, n) result(max_val) + ! -- dummy + class(LocalCellExtremaType), target :: this + integer(I4B), intent(in) :: n ! Node index + !-- return + real(DP), pointer :: max_val + + max_val => this%max(n) + end function get_max + + subroutine compute_local_extrema(this, phi) + ! -- dummy + class(LocalCellExtremaType), target :: this + real(DP), intent(in), dimension(:) :: phi + ! -- local + integer(I4B) :: n + real(DP) :: min_phi, max_phi + + do n = 1, this%dis%nodes + call this%find_local_extrema(n, phi, min_phi, max_phi) + this%min(n) = min_phi + this%max(n) = max_phi + end do + end subroutine compute_local_extrema + + subroutine find_local_extrema(this, n, phi, min_phi, max_phi) + ! -- dummy + class(LocalCellExtremaType), target :: 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 + +end module LocalCellExtremaModule diff --git a/src/meson.build b/src/meson.build index d8aba2dc023..670a9c1411e 100644 --- a/src/meson.build +++ b/src/meson.build @@ -282,6 +282,7 @@ modflow_sources = files( 'Model' / 'ModelUtilities' / 'VectorInterpolation.f90', 'Model' / 'ModelUtilities' / 'Xt3dAlgorithm.f90', 'Model' / 'ModelUtilities' / 'Xt3dInterface.f90', + 'Model' / 'TransportModel' / 'Gradient' / 'CachedGradient.f90', 'Model' / 'TransportModel' / 'Gradient' / 'IGradient.f90', 'Model' / 'TransportModel' / 'Gradient' / 'LeastSquaresGradient.f90', 'Model' / 'TransportModel' / 'InterpolationScheme' / 'AdvSchemeEnum.f90', @@ -441,6 +442,7 @@ modflow_sources = files( 'Utilities' / 'ListIterator.f90', 'Utilities' / 'ListNode.f90', 'Utilities' / 'ListReader.f90', + 'Utilities' / 'LocalCellExtrema.f90', 'Utilities' / 'LongLineReader.f90', 'Utilities' / 'MathUtil.f90', 'Utilities' / 'Message.f90',