(feat) Utvd optimization (#2427)

* Refactor tsp-adv. Move different interpolation scheme to separate files

* Add the SVD-algorithm and the pseudoinverse method

* Add the new unstructured TVD limiter

* Fix merge

* Add UTVD tests to existing unit tests. Fix correct stencil size on exchanges

* Apply review comments

* Fix lint error

* Apply review comments

* Apply review comment. Use pinv directly on the distance matrix omstead of the G matrix

* Optimize the limiter by caching several values

* Clean up

* Fix deconstructor call

* Move some function around

* Rename variable to better match documentation nomenclature

* Apply review comments.  Make the flow more explicit

Author: Manangka

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 \

msvs/mf6core.vfproj

@@ -390,6 +390,7 @@
 		<File RelativePath="..\src\Model\SurfaceWaterFlow\swf.f90"/></Filter>
 		<Filter Name="TransportModel">
 		<Filter Name="Gradient">
+		<File RelativePath="..\src\Model\TransportModel\Gradient\CachedGradient.f90"/>
 		<File RelativePath="..\src\Model\TransportModel\Gradient\IGradient.f90"/>
 		<File RelativePath="..\src\Model\TransportModel\Gradient\LeastSquaresGradient.f90"/></Filter>
 		<Filter Name="InterpolationScheme">
@@ -700,6 +701,7 @@
 		<File RelativePath="..\src\Utilities\ListIterator.f90"/>
 		<File RelativePath="..\src\Utilities\ListNode.f90"/>
 		<File RelativePath="..\src\Utilities\ListReader.f90"/>
+		<File RelativePath="..\src\Utilities\LocalCellExtrema.f90"/>
 		<File RelativePath="..\src\Utilities\LongLineReader.f90"/>
 		<File RelativePath="..\src\Utilities\MathUtil.f90"/>
 		<File RelativePath="..\src\Utilities\Message.f90"/>

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

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

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
 

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,15 +38,27 @@ 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
     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
 
@@ -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

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