diff --git a/make/makefile b/make/makefile index 3d9a270ff64..60d161ef20b 100644 --- a/make/makefile +++ b/make/makefile @@ -357,6 +357,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 \ @@ -388,6 +389,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 1dd10fb8aa0..4406f5d5cac 100644 --- a/msvs/mf6core.vfproj +++ b/msvs/mf6core.vfproj @@ -390,6 +390,7 @@ + @@ -700,6 +701,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 5e78ec0cd41..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 @@ -79,16 +87,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,34 +109,40 @@ 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 - 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 341a3f49e16..a16246b91ec 100644 --- a/src/Model/TransportModel/InterpolationScheme/CentralDifferenceScheme.f90 +++ b/src/Model/TransportModel/InterpolationScheme/CentralDifferenceScheme.f90 @@ -15,8 +15,10 @@ module CentralDifferenceSchemeModule private class(DisBaseType), pointer :: dis type(TspFmiType), pointer :: fmi + real(DP), dimension(:), pointer :: phi contains procedure :: compute + procedure :: set_field end type CentralDifferenceSchemeType interface CentralDifferenceSchemeType @@ -36,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 @@ -44,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 @@ -66,4 +80,5 @@ function compute(this, n, m, iposnm, phi) result(phi_face) phi_face%c_m = DONE - omega end function compute + end module CentralDifferenceSchemeModule diff --git a/src/Model/TransportModel/InterpolationScheme/InterpolationSchemeInterface.f90 b/src/Model/TransportModel/InterpolationScheme/InterpolationSchemeInterface.f90 index 0f4a853deb5..320132b8115 100644 --- a/src/Model/TransportModel/InterpolationScheme/InterpolationSchemeInterface.f90 +++ b/src/Model/TransportModel/InterpolationScheme/InterpolationSchemeInterface.f90 @@ -16,11 +16,18 @@ module InterpolationSchemeInterfaceModule type, abstract :: InterpolationSchemeInterface contains procedure(compute_if), deferred :: compute + 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 @@ -32,9 +39,23 @@ 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 + !> @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 end module InterpolationSchemeInterfaceModule diff --git a/src/Model/TransportModel/InterpolationScheme/TVDScheme.f90 b/src/Model/TransportModel/InterpolationScheme/TVDScheme.f90 index 9a8bcfeb9c1..5cd5ed2707c 100644 --- a/src/Model/TransportModel/InterpolationScheme/TVDScheme.f90 +++ b/src/Model/TransportModel/InterpolationScheme/TVDScheme.f90 @@ -15,9 +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 :: set_field end type TVDSchemeType interface TVDSchemeType @@ -39,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 @@ -47,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 @@ -96,14 +110,14 @@ 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 diff --git a/src/Model/TransportModel/InterpolationScheme/UTVDScheme.f90 b/src/Model/TransportModel/InterpolationScheme/UTVDScheme.f90 index 275dbd8ff08..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,12 +39,18 @@ 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 + real(DP), dimension(:, :), allocatable :: cached_node_distance ! distance vectors contains procedure :: compute + procedure :: set_field + final :: destructor - procedure, private :: find_local_extrema procedure, private :: limiter + procedure, private :: compute_node_distance end type UTVDSchemeType interface UTVDSchemeType @@ -50,22 +58,52 @@ module UTVDSchemeModule end interface UTVDSchemeType contains + 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 + class(IGradientType), allocatable, intent(in), target :: gradient interpolation_scheme%dis => dis interpolation_scheme%fmi => fmi + 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) end function constructor - function compute(this, n, m, iposnm, phi) result(phi_face) + subroutine destructor(this) + ! -- dummy + type(UTVDSchemeType), intent(inout) :: this + + deallocate (this%cached_node_distance) + end subroutine destructor + + !> @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 set_field(this, phi) + ! -- dummy + class(UTVDSchemeType), target :: this + real(DP), intent(in), dimension(:), pointer :: phi + + 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) result(phi_face) !-- return type(CoefficientsType), target :: phi_face ! Output: coefficients for the face between cells n and m ! -- dummy @@ -73,7 +111,6 @@ 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 @@ -85,7 +122,7 @@ function compute(this, n, m, iposnm, phi) result(phi_face) 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 isympos = this%dis%con%jas(iposnm) qnm = this%fmi%gwfflowja(iposnm) @@ -112,25 +149,36 @@ function compute(this, n, m, iposnm, phi) result(phi_face) 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 + dnm = this%cached_node_distance(isympos, :) + 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 + if (abs(this%phi(idn) - this%phi(iup)) < DSAME) return ! ! -- Compute cell concentration gradient - grad_c = this%gradient%get(iup, phi) + grad_c = this%gradient%get(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)) + 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. - call find_local_extrema(this, iup, phi, min_phi, max_phi) + 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 @@ -141,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: @@ -159,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 @@ -202,4 +232,23 @@ function limiter(this, r) result(theta) end select end function + 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 diff --git a/src/Model/TransportModel/InterpolationScheme/UpstreamScheme.f90 b/src/Model/TransportModel/InterpolationScheme/UpstreamScheme.f90 index 1729d1c6a22..a12728e15e4 100644 --- a/src/Model/TransportModel/InterpolationScheme/UpstreamScheme.f90 +++ b/src/Model/TransportModel/InterpolationScheme/UpstreamScheme.f90 @@ -15,8 +15,10 @@ module UpstreamSchemeModule private class(DisBaseType), pointer :: dis type(TspFmiType), pointer :: fmi + real(DP), dimension(:), pointer :: phi contains procedure :: compute + procedure :: set_field end type UpstreamSchemeType interface UpstreamSchemeType @@ -36,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 @@ -44,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 @@ -58,4 +72,5 @@ function compute(this, n, m, iposnm, phi) result(phi_face) end if end function compute + end module UpstreamSchemeModule diff --git a/src/Model/TransportModel/tsp-adv.f90 b/src/Model/TransportModel/tsp-adv.f90 index 9e6608fbee8..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,6 +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%set_field(cnew) do n = 1, nodes if (this%ibound(n) == 0) cycle ! skip inactive nodes idiag = this%dis%con%ia(n) @@ -240,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) @@ -255,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 @@ -267,6 +270,7 @@ subroutine adv_cq(this, cnew, flowja) ! rate and has dimensions of L^/T. nodes = this%dis%nodes + 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 @@ -274,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 1412f4ac141..7b043fe981d 100644 --- a/src/meson.build +++ b/src/meson.build @@ -288,6 +288,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', @@ -457,6 +458,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',