From c015ef1062bbc2940bab32180ae19f4d6ae80315 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Tue, 9 Jun 2026 09:56:20 -0700 Subject: [PATCH 1/5] operators: wire in plane slope operator build targets --- components/mpas-framework/src/operators/Makefile | 2 ++ components/mpas-framework/src/operators/operators.cmake | 1 + 2 files changed, 3 insertions(+) diff --git a/components/mpas-framework/src/operators/Makefile b/components/mpas-framework/src/operators/Makefile index b4d7085e90f7..7f253275269d 100644 --- a/components/mpas-framework/src/operators/Makefile +++ b/components/mpas-framework/src/operators/Makefile @@ -5,6 +5,7 @@ OBJS = mpas_vector_operations.o \ mpas_tensor_operations.o \ mpas_rbf_interpolation.o \ mpas_vector_reconstruction.o \ + mpas_plane_slope_operators.o \ mpas_spline_interpolation.o \ mpas_tracer_advection_helpers.o \ mpas_tracer_advection_mono.o \ @@ -23,6 +24,7 @@ mpas_matrix_operations.o: $(DEPS) mpas_tensor_operations.o: mpas_vector_operations.o mpas_matrix_operations.o $(DEPS) mpas_rbf_interpolation.o: mpas_vector_operations.o mpas_vector_reconstruction.o: mpas_rbf_interpolation.o +mpas_plane_slope_operators.o: mpas_spline_interpolation: mpas_tracer_advection_helpers.o: mpas_geometry_utils.o $(DEPS) mpas_tracer_advection_mono.o: mpas_tracer_advection_helpers.o diff --git a/components/mpas-framework/src/operators/operators.cmake b/components/mpas-framework/src/operators/operators.cmake index d65c7c661eb3..79e670a46690 100644 --- a/components/mpas-framework/src/operators/operators.cmake +++ b/components/mpas-framework/src/operators/operators.cmake @@ -5,6 +5,7 @@ list(APPEND COMMON_RAW_SOURCES operators/mpas_tensor_operations.F operators/mpas_rbf_interpolation.F operators/mpas_vector_reconstruction.F + operators/mpas_plane_slope_operators.F operators/mpas_spline_interpolation.F operators/mpas_tracer_advection_helpers.F operators/mpas_tracer_advection_mono.F From 0f91c5329ffd94d8a3af6929cc75abf5020dc7c8 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Tue, 9 Jun 2026 09:57:48 -0700 Subject: [PATCH 2/5] operators: add ring-based planar max-slope plane-fit operator --- .../operators/mpas_plane_slope_operators.F | 337 ++++++++++++++++++ 1 file changed, 337 insertions(+) create mode 100644 components/mpas-framework/src/operators/mpas_plane_slope_operators.F diff --git a/components/mpas-framework/src/operators/mpas_plane_slope_operators.F b/components/mpas-framework/src/operators/mpas_plane_slope_operators.F new file mode 100644 index 000000000000..961ef0505d82 --- /dev/null +++ b/components/mpas-framework/src/operators/mpas_plane_slope_operators.F @@ -0,0 +1,337 @@ +! Copyright (c) 2013, Los Alamos National Security, LLC (LANS) +! and the University Corporation for Atmospheric Research (UCAR). +! +! Unless noted otherwise source code is licensed under the BSD license. +! Additional copyright and license information can be found in the LICENSE file +! distributed with this code, or at http://mpas-dev.github.io/license.html +! +!*********************************************************************** +! +! mpas_plane_slope_operators +! +!> \brief MPAS plane-fit slope operators +!> \author MPAS Development Team +!> \date 06/09/26 +!> \details +!> This module computes per-cell unsigned slope magnitude from a linear +!> plane fit over a contiguous ring stencil on planar meshes. +! +!----------------------------------------------------------------------- +module mpas_plane_slope_operators + + use mpas_kind_types + use mpas_derived_types + use mpas_pool_routines + use mpas_log + + implicit none + private + save + + public :: mpas_max_plane_slope_ring + +contains + +!*********************************************************************** +! +! routine mpas_max_plane_slope_ring +! +!> \brief Compute unsigned plane-fit slope on configurable ring stencils +!> \author MPAS Development Team +!> \date 06/09/26 +!> \details +!> For each owned cell, this routine fits a linear plane to neighboring +!> cell-center values over a contiguous ring stencil and returns the +!> slope magnitude sqrt(a^2+b^2) where z = a*x + b*y + c in local +!> cell-centered coordinates. The routine is planar-mesh only. +! +!----------------------------------------------------------------------- + subroutine mpas_max_plane_slope_ring(fieldCell, meshPool, configPool, ringCount, maxSlope, validMask)!{{{ + + real(kind=RKIND), dimension(:), intent(in) :: fieldCell + type(mpas_pool_type), intent(in) :: meshPool + type(mpas_pool_type), intent(in) :: configPool + integer, intent(in) :: ringCount + real(kind=RKIND), dimension(:), intent(out) :: maxSlope + logical, dimension(:), intent(in), optional :: validMask + + integer, pointer :: nCells, nCellsSolve + integer, dimension(:), pointer :: nEdgesOnCell + integer, dimension(:,:), pointer :: cellsOnCell + real(kind=RKIND), dimension(:), pointer :: xCell, yCell + logical, pointer :: on_a_sphere + integer, pointer :: config_num_halos + + integer, allocatable, dimension(:) :: visited, ringByCell, frontier, nextFrontier, stencil + integer :: iCell, i, iRing, iFront, iNbr, iCellNbr + integer :: frontierCount, nextFrontierCount, stencilCount + integer :: maxBuiltRing, effectiveRing + integer :: visitTag + real(kind=RKIND) :: aCoeff, bCoeff + logical :: fitOK + + call mpas_pool_get_dimension(meshPool, 'nCells', nCells) + call mpas_pool_get_dimension(meshPool, 'nCellsSolve', nCellsSolve) + + call mpas_pool_get_array(meshPool, 'nEdgesOnCell', nEdgesOnCell) + call mpas_pool_get_array(meshPool, 'cellsOnCell', cellsOnCell) + call mpas_pool_get_array(meshPool, 'xCell', xCell) + call mpas_pool_get_array(meshPool, 'yCell', yCell) + call mpas_pool_get_config(meshPool, 'on_a_sphere', on_a_sphere) + call mpas_pool_get_config(configPool, 'config_num_halos', config_num_halos) + + if (ringCount < 1) then + call mpas_log_write('mpas_max_plane_slope_ring: ringCount must be >= 1.', MPAS_LOG_CRIT) + end if + + if (ringCount > config_num_halos) then + call mpas_log_write('mpas_max_plane_slope_ring: requested ring count $i exceeds config_num_halos $i.', & + MPAS_LOG_CRIT, intArgs=(/ringCount, config_num_halos/)) + end if + + if (on_a_sphere) then + call mpas_log_write('mpas_max_plane_slope_ring currently supports planar meshes only.', MPAS_LOG_CRIT) + end if + + if (size(maxSlope) < nCellsSolve) then + call mpas_log_write('mpas_max_plane_slope_ring: output array maxSlope is too small.', MPAS_LOG_CRIT) + end if + + if (size(fieldCell) < nCells) then + call mpas_log_write('mpas_max_plane_slope_ring: input array fieldCell is too small.', MPAS_LOG_CRIT) + end if + + if (present(validMask)) then + if (size(validMask) < nCells) then + call mpas_log_write('mpas_max_plane_slope_ring: optional validMask is too small.', MPAS_LOG_CRIT) + end if + end if + + maxSlope = 0.0_RKIND + + allocate(visited(nCells)) + allocate(ringByCell(nCells)) + allocate(frontier(nCells)) + allocate(nextFrontier(nCells)) + allocate(stencil(nCells)) + + visited = 0 + ringByCell = -1 + visitTag = 0 + + do iCell = 1, nCellsSolve + + if (present(validMask)) then + if (.not. validMask(iCell)) then + maxSlope(iCell) = 0.0_RKIND + cycle + end if + end if + + visitTag = visitTag + 1 + if (visitTag == huge(visitTag)) then + visited = 0 + visitTag = 1 + end if + + frontierCount = 1 + frontier(1) = iCell + visited(iCell) = visitTag + ringByCell(iCell) = 0 + + stencilCount = 1 + stencil(1) = iCell + maxBuiltRing = 0 + + do iRing = 1, ringCount + nextFrontierCount = 0 + + do iFront = 1, frontierCount + do iNbr = 1, nEdgesOnCell(frontier(iFront)) + iCellNbr = cellsOnCell(iNbr, frontier(iFront)) + + if (iCellNbr < 1 .or. iCellNbr > nCells) cycle + if (visited(iCellNbr) == visitTag) cycle + + if (present(validMask)) then + if (.not. validMask(iCellNbr)) cycle + end if + + visited(iCellNbr) = visitTag + ringByCell(iCellNbr) = iRing + + nextFrontierCount = nextFrontierCount + 1 + nextFrontier(nextFrontierCount) = iCellNbr + + stencilCount = stencilCount + 1 + stencil(stencilCount) = iCellNbr + end do + end do + + if (nextFrontierCount == 0) exit + + frontier(1:nextFrontierCount) = nextFrontier(1:nextFrontierCount) + frontierCount = nextFrontierCount + maxBuiltRing = iRing + end do + + effectiveRing = maxBuiltRing + fitOK = .false. + + do while (effectiveRing >= 1 .and. .not. fitOK) + call mpas_try_plane_fit(iCell, effectiveRing, stencil, stencilCount, ringByCell, & + xCell, yCell, fieldCell, aCoeff, bCoeff, fitOK) + if (.not. fitOK) effectiveRing = effectiveRing - 1 + end do + + if (fitOK) then + maxSlope(iCell) = sqrt(aCoeff*aCoeff + bCoeff*bCoeff) + else + maxSlope(iCell) = 0.0_RKIND + end if + + end do + + deallocate(visited) + deallocate(ringByCell) + deallocate(frontier) + deallocate(nextFrontier) + deallocate(stencil) + + end subroutine mpas_max_plane_slope_ring!}}} + + +!*********************************************************************** +! +! routine mpas_try_plane_fit +! +!> \brief Attempt linear plane fit over a reduced ring stencil +!> \author MPAS Development Team +!> \date 06/09/26 +!> \details +!> Fits dz = a*dx + b*dy + c and returns (a,b) when matrix is well posed. +! +!----------------------------------------------------------------------- + subroutine mpas_try_plane_fit(iCenter, maxRing, stencil, stencilCount, ringByCell, xCell, yCell, fieldCell, aCoeff, bCoeff, fitOK)!{{{ + + integer, intent(in) :: iCenter, maxRing, stencilCount + integer, dimension(:), intent(in) :: stencil, ringByCell + real(kind=RKIND), dimension(:), intent(in) :: xCell, yCell, fieldCell + real(kind=RKIND), intent(out) :: aCoeff, bCoeff + logical, intent(out) :: fitOK + + integer :: i, iCell + integer :: nPoints + real(kind=RKIND) :: x0, y0, z0, dx, dy, dz + real(kind=RKIND), dimension(3,3) :: mat + real(kind=RKIND), dimension(3) :: rhs, soln + real(kind=RKIND) :: detScale + + fitOK = .false. + aCoeff = 0.0_RKIND + bCoeff = 0.0_RKIND + + x0 = xCell(iCenter) + y0 = yCell(iCenter) + z0 = fieldCell(iCenter) + + mat = 0.0_RKIND + rhs = 0.0_RKIND + nPoints = 0 + + do i = 1, stencilCount + iCell = stencil(i) + if (ringByCell(iCell) > maxRing) cycle + + dx = xCell(iCell) - x0 + dy = yCell(iCell) - y0 + dz = fieldCell(iCell) - z0 + + mat(1,1) = mat(1,1) + dx*dx + mat(1,2) = mat(1,2) + dx*dy + mat(1,3) = mat(1,3) + dx + + mat(2,2) = mat(2,2) + dy*dy + mat(2,3) = mat(2,3) + dy + + mat(3,3) = mat(3,3) + 1.0_RKIND + + rhs(1) = rhs(1) + dx*dz + rhs(2) = rhs(2) + dy*dz + rhs(3) = rhs(3) + dz + + nPoints = nPoints + 1 + end do + + if (nPoints < 3) return + + mat(2,1) = mat(1,2) + mat(3,1) = mat(1,3) + mat(3,2) = mat(2,3) + + detScale = max(1.0_RKIND, maxval(abs(mat))) + call mpas_solve_3x3(mat, rhs, detScale*1.0e-12_RKIND, soln, fitOK) + if (.not. fitOK) return + + aCoeff = soln(1) + bCoeff = soln(2) + + end subroutine mpas_try_plane_fit!}}} + + +!*********************************************************************** +! +! routine mpas_solve_3x3 +! +!> \brief Solve a 3x3 linear system with Cramer's rule +!> \author MPAS Development Team +!> \date 06/09/26 +! +!----------------------------------------------------------------------- + subroutine mpas_solve_3x3(a, b, detTol, x, ok)!{{{ + + real(kind=RKIND), dimension(3,3), intent(in) :: a + real(kind=RKIND), dimension(3), intent(in) :: b + real(kind=RKIND), intent(in) :: detTol + real(kind=RKIND), dimension(3), intent(out) :: x + logical, intent(out) :: ok + + real(kind=RKIND) :: detA, detX1, detX2, detX3 + real(kind=RKIND), dimension(3,3) :: work + + detA = a(1,1)*(a(2,2)*a(3,3) - a(2,3)*a(3,2)) - & + a(1,2)*(a(2,1)*a(3,3) - a(2,3)*a(3,1)) + & + a(1,3)*(a(2,1)*a(3,2) - a(2,2)*a(3,1)) + + if (abs(detA) <= detTol) then + ok = .false. + x = 0.0_RKIND + return + end if + + work = a + work(:,1) = b + detX1 = work(1,1)*(work(2,2)*work(3,3) - work(2,3)*work(3,2)) - & + work(1,2)*(work(2,1)*work(3,3) - work(2,3)*work(3,1)) + & + work(1,3)*(work(2,1)*work(3,2) - work(2,2)*work(3,1)) + + work = a + work(:,2) = b + detX2 = work(1,1)*(work(2,2)*work(3,3) - work(2,3)*work(3,2)) - & + work(1,2)*(work(2,1)*work(3,3) - work(2,3)*work(3,1)) + & + work(1,3)*(work(2,1)*work(3,2) - work(2,2)*work(3,1)) + + work = a + work(:,3) = b + detX3 = work(1,1)*(work(2,2)*work(3,3) - work(2,3)*work(3,2)) - & + work(1,2)*(work(2,1)*work(3,3) - work(2,3)*work(3,1)) + & + work(1,3)*(work(2,1)*work(3,2) - work(2,2)*work(3,1)) + + x(1) = detX1 / detA + x(2) = detX2 / detA + x(3) = detX3 / detA + ok = .true. + + end subroutine mpas_solve_3x3!}}} + +end module mpas_plane_slope_operators From 765351b684f9cbd4ac6af60d395f9cb49fcd28b4 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Tue, 9 Jun 2026 11:27:46 -0700 Subject: [PATCH 3/5] Add comments --- .../operators/mpas_plane_slope_operators.F | 41 ++++++++++++------- 1 file changed, 26 insertions(+), 15 deletions(-) diff --git a/components/mpas-framework/src/operators/mpas_plane_slope_operators.F b/components/mpas-framework/src/operators/mpas_plane_slope_operators.F index 961ef0505d82..5d1df08ed7ce 100644 --- a/components/mpas-framework/src/operators/mpas_plane_slope_operators.F +++ b/components/mpas-framework/src/operators/mpas_plane_slope_operators.F @@ -48,12 +48,12 @@ module mpas_plane_slope_operators !----------------------------------------------------------------------- subroutine mpas_max_plane_slope_ring(fieldCell, meshPool, configPool, ringCount, maxSlope, validMask)!{{{ - real(kind=RKIND), dimension(:), intent(in) :: fieldCell - type(mpas_pool_type), intent(in) :: meshPool - type(mpas_pool_type), intent(in) :: configPool - integer, intent(in) :: ringCount - real(kind=RKIND), dimension(:), intent(out) :: maxSlope - logical, dimension(:), intent(in), optional :: validMask + real(kind=RKIND), dimension(:), intent(in) :: fieldCell ! Cell-centered scalar field to fit + type(mpas_pool_type), intent(in) :: meshPool ! Mesh geometry/connectivity pool + type(mpas_pool_type), intent(in) :: configPool ! Config pool (used for halo-depth guard) + integer, intent(in) :: ringCount ! Requested contiguous ring depth + real(kind=RKIND), dimension(:), intent(out) :: maxSlope ! Output unsigned slope magnitude per solve cell + logical, dimension(:), intent(in), optional :: validMask ! Optional mask of cells eligible for stencils integer, pointer :: nCells, nCellsSolve integer, dimension(:), pointer :: nEdgesOnCell @@ -70,6 +70,7 @@ subroutine mpas_max_plane_slope_ring(fieldCell, meshPool, configPool, ringCount, real(kind=RKIND) :: aCoeff, bCoeff logical :: fitOK + ! Gather mesh/config metadata and connectivity needed to build local ring stencils. call mpas_pool_get_dimension(meshPool, 'nCells', nCells) call mpas_pool_get_dimension(meshPool, 'nCellsSolve', nCellsSolve) @@ -80,6 +81,7 @@ subroutine mpas_max_plane_slope_ring(fieldCell, meshPool, configPool, ringCount, call mpas_pool_get_config(meshPool, 'on_a_sphere', on_a_sphere) call mpas_pool_get_config(configPool, 'config_num_halos', config_num_halos) + ! Input and configuration validation. if (ringCount < 1) then call mpas_log_write('mpas_max_plane_slope_ring: ringCount must be >= 1.', MPAS_LOG_CRIT) end if @@ -109,6 +111,7 @@ subroutine mpas_max_plane_slope_ring(fieldCell, meshPool, configPool, ringCount, maxSlope = 0.0_RKIND + ! Work arrays are sized once to avoid repeated allocation in the per-cell loop. allocate(visited(nCells)) allocate(ringByCell(nCells)) allocate(frontier(nCells)) @@ -121,6 +124,7 @@ subroutine mpas_max_plane_slope_ring(fieldCell, meshPool, configPool, ringCount, do iCell = 1, nCellsSolve + ! If the center cell is masked out, return the requested fallback immediately. if (present(validMask)) then if (.not. validMask(iCell)) then maxSlope(iCell) = 0.0_RKIND @@ -143,6 +147,7 @@ subroutine mpas_max_plane_slope_ring(fieldCell, meshPool, configPool, ringCount, stencil(1) = iCell maxBuiltRing = 0 + ! Breadth-first expansion over cellsOnCell builds a contiguous stencil up to ringCount. do iRing = 1, ringCount nextFrontierCount = 0 @@ -175,6 +180,7 @@ subroutine mpas_max_plane_slope_ring(fieldCell, meshPool, configPool, ringCount, maxBuiltRing = iRing end do + ! Reduce outer rings until a stable linear fit is available. effectiveRing = maxBuiltRing fitOK = .false. @@ -192,6 +198,7 @@ subroutine mpas_max_plane_slope_ring(fieldCell, meshPool, configPool, ringCount, end do + ! Explicit cleanup of workspace arrays. deallocate(visited) deallocate(ringByCell) deallocate(frontier) @@ -214,11 +221,11 @@ end subroutine mpas_max_plane_slope_ring!}}} !----------------------------------------------------------------------- subroutine mpas_try_plane_fit(iCenter, maxRing, stencil, stencilCount, ringByCell, xCell, yCell, fieldCell, aCoeff, bCoeff, fitOK)!{{{ - integer, intent(in) :: iCenter, maxRing, stencilCount - integer, dimension(:), intent(in) :: stencil, ringByCell + integer, intent(in) :: iCenter, maxRing, stencilCount ! Center cell index and reduced-ring depth + integer, dimension(:), intent(in) :: stencil, ringByCell ! Candidate stencil cells and ring assignment real(kind=RKIND), dimension(:), intent(in) :: xCell, yCell, fieldCell - real(kind=RKIND), intent(out) :: aCoeff, bCoeff - logical, intent(out) :: fitOK + real(kind=RKIND), intent(out) :: aCoeff, bCoeff ! Plane slopes dz/dx and dz/dy + logical, intent(out) :: fitOK ! .true. when solve is well-posed integer :: i, iCell integer :: nPoints @@ -239,6 +246,7 @@ subroutine mpas_try_plane_fit(iCenter, maxRing, stencil, stencilCount, ringByCel rhs = 0.0_RKIND nPoints = 0 + ! Build normal equations for dz = a*dx + b*dy + c in center-relative coordinates. do i = 1, stencilCount iCell = stencil(i) if (ringByCell(iCell) > maxRing) cycle @@ -263,12 +271,14 @@ subroutine mpas_try_plane_fit(iCenter, maxRing, stencil, stencilCount, ringByCel nPoints = nPoints + 1 end do + ! Fewer than 3 points cannot constrain a linear plane robustly. if (nPoints < 3) return mat(2,1) = mat(1,2) mat(3,1) = mat(1,3) mat(3,2) = mat(2,3) + ! Use matrix scale to set a relative singularity tolerance. detScale = max(1.0_RKIND, maxval(abs(mat))) call mpas_solve_3x3(mat, rhs, detScale*1.0e-12_RKIND, soln, fitOK) if (.not. fitOK) return @@ -290,15 +300,16 @@ end subroutine mpas_try_plane_fit!}}} !----------------------------------------------------------------------- subroutine mpas_solve_3x3(a, b, detTol, x, ok)!{{{ - real(kind=RKIND), dimension(3,3), intent(in) :: a - real(kind=RKIND), dimension(3), intent(in) :: b - real(kind=RKIND), intent(in) :: detTol - real(kind=RKIND), dimension(3), intent(out) :: x - logical, intent(out) :: ok + real(kind=RKIND), dimension(3,3), intent(in) :: a ! Coefficient matrix + real(kind=RKIND), dimension(3), intent(in) :: b ! Right-hand side + real(kind=RKIND), intent(in) :: detTol ! Determinant tolerance for singularity + real(kind=RKIND), dimension(3), intent(out) :: x ! Solution vector + logical, intent(out) :: ok ! .true. when matrix is invertible real(kind=RKIND) :: detA, detX1, detX2, detX3 real(kind=RKIND), dimension(3,3) :: work + ! Cramer's rule is sufficient here because the system is fixed-size (3x3). detA = a(1,1)*(a(2,2)*a(3,3) - a(2,3)*a(3,2)) - & a(1,2)*(a(2,1)*a(3,3) - a(2,3)*a(3,1)) + & a(1,3)*(a(2,1)*a(3,2) - a(2,2)*a(3,1)) From 1cb9f19a5adb2897eb1e0b594ee69641ebd4e26d Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Tue, 9 Jun 2026 11:52:41 -0700 Subject: [PATCH 4/5] Refactor plane slope operator internals for shared stencil and fit helpers --- .../operators/mpas_plane_slope_operators.F | 194 ++++++++++++------ 1 file changed, 136 insertions(+), 58 deletions(-) diff --git a/components/mpas-framework/src/operators/mpas_plane_slope_operators.F b/components/mpas-framework/src/operators/mpas_plane_slope_operators.F index 5d1df08ed7ce..a69729f490ff 100644 --- a/components/mpas-framework/src/operators/mpas_plane_slope_operators.F +++ b/components/mpas-framework/src/operators/mpas_plane_slope_operators.F @@ -63,10 +63,10 @@ subroutine mpas_max_plane_slope_ring(fieldCell, meshPool, configPool, ringCount, integer, pointer :: config_num_halos integer, allocatable, dimension(:) :: visited, ringByCell, frontier, nextFrontier, stencil - integer :: iCell, i, iRing, iFront, iNbr, iCellNbr - integer :: frontierCount, nextFrontierCount, stencilCount + integer :: iCell, stencilCount integer :: maxBuiltRing, effectiveRing integer :: visitTag + integer :: seedCells(1) real(kind=RKIND) :: aCoeff, bCoeff logical :: fitOK @@ -132,53 +132,12 @@ subroutine mpas_max_plane_slope_ring(fieldCell, meshPool, configPool, ringCount, end if end if - visitTag = visitTag + 1 - if (visitTag == huge(visitTag)) then - visited = 0 - visitTag = 1 - end if - - frontierCount = 1 - frontier(1) = iCell - visited(iCell) = visitTag - ringByCell(iCell) = 0 - - stencilCount = 1 - stencil(1) = iCell - maxBuiltRing = 0 - - ! Breadth-first expansion over cellsOnCell builds a contiguous stencil up to ringCount. - do iRing = 1, ringCount - nextFrontierCount = 0 - - do iFront = 1, frontierCount - do iNbr = 1, nEdgesOnCell(frontier(iFront)) - iCellNbr = cellsOnCell(iNbr, frontier(iFront)) - - if (iCellNbr < 1 .or. iCellNbr > nCells) cycle - if (visited(iCellNbr) == visitTag) cycle - - if (present(validMask)) then - if (.not. validMask(iCellNbr)) cycle - end if - - visited(iCellNbr) = visitTag - ringByCell(iCellNbr) = iRing + seedCells(1) = iCell + call mpas_build_contiguous_cell_stencil(seedCells, 1, ringCount, nCells, nEdgesOnCell, cellsOnCell, & + visited, ringByCell, frontier, nextFrontier, stencil, visitTag, & + stencilCount, maxBuiltRing, validMask) - nextFrontierCount = nextFrontierCount + 1 - nextFrontier(nextFrontierCount) = iCellNbr - - stencilCount = stencilCount + 1 - stencil(stencilCount) = iCellNbr - end do - end do - - if (nextFrontierCount == 0) exit - - frontier(1:nextFrontierCount) = nextFrontier(1:nextFrontierCount) - frontierCount = nextFrontierCount - maxBuiltRing = iRing - end do + if (stencilCount <= 0) cycle ! Reduce outer rings until a stable linear fit is available. effectiveRing = maxBuiltRing @@ -208,6 +167,98 @@ subroutine mpas_max_plane_slope_ring(fieldCell, meshPool, configPool, ringCount, end subroutine mpas_max_plane_slope_ring!}}} +!*********************************************************************** +! +! routine mpas_build_contiguous_cell_stencil +! +!> \brief Build a mask-aware contiguous ring stencil from one or more seed cells +!> \author MPAS Development Team +!> \date 06/09/26 +! +!----------------------------------------------------------------------- + subroutine mpas_build_contiguous_cell_stencil(seedCells, nSeed, ringCount, nCells, nEdgesOnCell, cellsOnCell, & + visited, ringByCell, frontier, nextFrontier, stencil, visitTag, & + stencilCount, maxBuiltRing, validMask)!{{{ + + integer, dimension(:), intent(in) :: seedCells + integer, intent(in) :: nSeed, ringCount, nCells + integer, dimension(:), intent(in) :: nEdgesOnCell + integer, dimension(:,:), intent(in) :: cellsOnCell + integer, dimension(:), intent(inout) :: visited, ringByCell, frontier, nextFrontier, stencil + integer, intent(inout) :: visitTag + integer, intent(out) :: stencilCount, maxBuiltRing + logical, dimension(:), intent(in), optional :: validMask + + integer :: iSeed, iRing, iFront, iNbr, iCell, iCellNbr + integer :: frontierCount, nextFrontierCount + + visitTag = visitTag + 1 + if (visitTag == huge(visitTag)) then + visited = 0 + visitTag = 1 + end if + + frontierCount = 0 + stencilCount = 0 + maxBuiltRing = 0 + + do iSeed = 1, nSeed + iCell = seedCells(iSeed) + + if (iCell < 1 .or. iCell > nCells) cycle + if (visited(iCell) == visitTag) cycle + + if (present(validMask)) then + if (.not. validMask(iCell)) cycle + end if + + visited(iCell) = visitTag + ringByCell(iCell) = 0 + + frontierCount = frontierCount + 1 + frontier(frontierCount) = iCell + + stencilCount = stencilCount + 1 + stencil(stencilCount) = iCell + end do + + if (frontierCount == 0) return + + do iRing = 1, ringCount + nextFrontierCount = 0 + + do iFront = 1, frontierCount + do iNbr = 1, nEdgesOnCell(frontier(iFront)) + iCellNbr = cellsOnCell(iNbr, frontier(iFront)) + + if (iCellNbr < 1 .or. iCellNbr > nCells) cycle + if (visited(iCellNbr) == visitTag) cycle + + if (present(validMask)) then + if (.not. validMask(iCellNbr)) cycle + end if + + visited(iCellNbr) = visitTag + ringByCell(iCellNbr) = iRing + + nextFrontierCount = nextFrontierCount + 1 + nextFrontier(nextFrontierCount) = iCellNbr + + stencilCount = stencilCount + 1 + stencil(stencilCount) = iCellNbr + end do + end do + + if (nextFrontierCount == 0) exit + + frontier(1:nextFrontierCount) = nextFrontier(1:nextFrontierCount) + frontierCount = nextFrontierCount + maxBuiltRing = iRing + end do + + end subroutine mpas_build_contiguous_cell_stencil!}}} + + !*********************************************************************** ! ! routine mpas_try_plane_fit @@ -227,9 +278,40 @@ subroutine mpas_try_plane_fit(iCenter, maxRing, stencil, stencilCount, ringByCel real(kind=RKIND), intent(out) :: aCoeff, bCoeff ! Plane slopes dz/dx and dz/dy logical, intent(out) :: fitOK ! .true. when solve is well-posed + real(kind=RKIND) :: x0, y0, z0 + + x0 = xCell(iCenter) + y0 = yCell(iCenter) + z0 = fieldCell(iCenter) + + call mpas_try_plane_fit_reference(x0, y0, z0, maxRing, stencil, stencilCount, ringByCell, & + xCell, yCell, fieldCell, aCoeff, bCoeff, fitOK) + + end subroutine mpas_try_plane_fit!}}} + + +!*********************************************************************** +! +! routine mpas_try_plane_fit_reference +! +!> \brief Attempt linear plane fit over a reduced ring stencil around a reference point +!> \author MPAS Development Team +!> \date 06/09/26 +! +!----------------------------------------------------------------------- + subroutine mpas_try_plane_fit_reference(xRef, yRef, zRef, maxRing, stencil, stencilCount, ringByCell, & + xCell, yCell, fieldCell, aCoeff, bCoeff, fitOK)!{{{ + + real(kind=RKIND), intent(in) :: xRef, yRef, zRef + integer, intent(in) :: maxRing, stencilCount + integer, dimension(:), intent(in) :: stencil, ringByCell + real(kind=RKIND), dimension(:), intent(in) :: xCell, yCell, fieldCell + real(kind=RKIND), intent(out) :: aCoeff, bCoeff + logical, intent(out) :: fitOK + integer :: i, iCell integer :: nPoints - real(kind=RKIND) :: x0, y0, z0, dx, dy, dz + real(kind=RKIND) :: dx, dy, dz real(kind=RKIND), dimension(3,3) :: mat real(kind=RKIND), dimension(3) :: rhs, soln real(kind=RKIND) :: detScale @@ -238,22 +320,18 @@ subroutine mpas_try_plane_fit(iCenter, maxRing, stencil, stencilCount, ringByCel aCoeff = 0.0_RKIND bCoeff = 0.0_RKIND - x0 = xCell(iCenter) - y0 = yCell(iCenter) - z0 = fieldCell(iCenter) - mat = 0.0_RKIND rhs = 0.0_RKIND nPoints = 0 - ! Build normal equations for dz = a*dx + b*dy + c in center-relative coordinates. + ! Build normal equations for dz = a*dx + b*dy + c in reference-relative coordinates. do i = 1, stencilCount iCell = stencil(i) if (ringByCell(iCell) > maxRing) cycle - dx = xCell(iCell) - x0 - dy = yCell(iCell) - y0 - dz = fieldCell(iCell) - z0 + dx = xCell(iCell) - xRef + dy = yCell(iCell) - yRef + dz = fieldCell(iCell) - zRef mat(1,1) = mat(1,1) + dx*dx mat(1,2) = mat(1,2) + dx*dy @@ -286,7 +364,7 @@ subroutine mpas_try_plane_fit(iCenter, maxRing, stencil, stencilCount, ringByCel aCoeff = soln(1) bCoeff = soln(2) - end subroutine mpas_try_plane_fit!}}} + end subroutine mpas_try_plane_fit_reference!}}} !*********************************************************************** From 92c03a0514d86367bd27544ab6aaf4af68056957 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Tue, 9 Jun 2026 11:53:23 -0700 Subject: [PATCH 5/5] Add edge-based plane-fit normal and tangent slope operator --- .../operators/mpas_plane_slope_operators.F | 147 ++++++++++++++++++ 1 file changed, 147 insertions(+) diff --git a/components/mpas-framework/src/operators/mpas_plane_slope_operators.F b/components/mpas-framework/src/operators/mpas_plane_slope_operators.F index a69729f490ff..e77c891dd44a 100644 --- a/components/mpas-framework/src/operators/mpas_plane_slope_operators.F +++ b/components/mpas-framework/src/operators/mpas_plane_slope_operators.F @@ -29,6 +29,7 @@ module mpas_plane_slope_operators save public :: mpas_max_plane_slope_ring + public :: mpas_edge_plane_slopes_ring contains @@ -167,6 +168,152 @@ subroutine mpas_max_plane_slope_ring(fieldCell, meshPool, configPool, ringCount, end subroutine mpas_max_plane_slope_ring!}}} +!*********************************************************************** +! +! routine mpas_edge_plane_slopes_ring +! +!> \brief Compute signed normal and tangent plane-fit slopes on edges +!> \author MPAS Development Team +!> \date 06/09/26 +!> \details +!> For each owned edge, this routine fits a linear plane to neighboring +!> cell-center values over a contiguous ring stencil around the edge and +!> returns the projected slope components along edge-normal and edge- +!> tangent directions. The routine is planar-mesh only. +! +!----------------------------------------------------------------------- + subroutine mpas_edge_plane_slopes_ring(fieldCell, meshPool, configPool, ringCount, normalSlopeEdge, & + tangentSlopeEdge, validMask)!{{{ + + real(kind=RKIND), dimension(:), intent(in) :: fieldCell + type(mpas_pool_type), intent(in) :: meshPool + type(mpas_pool_type), intent(in) :: configPool + integer, intent(in) :: ringCount + real(kind=RKIND), dimension(:), intent(out) :: normalSlopeEdge + real(kind=RKIND), dimension(:), intent(out) :: tangentSlopeEdge + logical, dimension(:), intent(in), optional :: validMask + + integer, pointer :: nCells, nEdgesSolve + integer, dimension(:), pointer :: nEdgesOnCell + integer, dimension(:,:), pointer :: cellsOnCell, cellsOnEdge + real(kind=RKIND), dimension(:), pointer :: xCell, yCell, xEdge, yEdge, angleEdge + logical, pointer :: on_a_sphere + integer, pointer :: config_num_halos + + integer, allocatable, dimension(:) :: visited, ringByCell, frontier, nextFrontier, stencil + integer :: iEdge, iCell1, iCell2, stencilCount + integer :: maxBuiltRing, effectiveRing + integer :: visitTag + integer :: seedCells(2) + real(kind=RKIND) :: xRef, yRef, zRef + real(kind=RKIND) :: aCoeff, bCoeff + real(kind=RKIND) :: edgeCos, edgeSin + logical :: fitOK + + call mpas_pool_get_dimension(meshPool, 'nCells', nCells) + call mpas_pool_get_dimension(meshPool, 'nEdgesSolve', nEdgesSolve) + + call mpas_pool_get_array(meshPool, 'nEdgesOnCell', nEdgesOnCell) + call mpas_pool_get_array(meshPool, 'cellsOnCell', cellsOnCell) + call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge) + call mpas_pool_get_array(meshPool, 'xCell', xCell) + call mpas_pool_get_array(meshPool, 'yCell', yCell) + call mpas_pool_get_array(meshPool, 'xEdge', xEdge) + call mpas_pool_get_array(meshPool, 'yEdge', yEdge) + call mpas_pool_get_array(meshPool, 'angleEdge', angleEdge) + call mpas_pool_get_config(meshPool, 'on_a_sphere', on_a_sphere) + call mpas_pool_get_config(configPool, 'config_num_halos', config_num_halos) + + if (ringCount < 1) then + call mpas_log_write('mpas_edge_plane_slopes_ring: ringCount must be >= 1.', MPAS_LOG_CRIT) + end if + + if (ringCount > config_num_halos) then + call mpas_log_write('mpas_edge_plane_slopes_ring: requested ring count $i exceeds config_num_halos $i.', & + MPAS_LOG_CRIT, intArgs=(/ringCount, config_num_halos/)) + end if + + if (on_a_sphere) then + call mpas_log_write('mpas_edge_plane_slopes_ring currently supports planar meshes only.', MPAS_LOG_CRIT) + end if + + if (size(normalSlopeEdge) < nEdgesSolve) then + call mpas_log_write('mpas_edge_plane_slopes_ring: output array normalSlopeEdge is too small.', MPAS_LOG_CRIT) + end if + + if (size(tangentSlopeEdge) < nEdgesSolve) then + call mpas_log_write('mpas_edge_plane_slopes_ring: output array tangentSlopeEdge is too small.', MPAS_LOG_CRIT) + end if + + if (size(fieldCell) < nCells) then + call mpas_log_write('mpas_edge_plane_slopes_ring: input array fieldCell is too small.', MPAS_LOG_CRIT) + end if + + if (present(validMask)) then + if (size(validMask) < nCells) then + call mpas_log_write('mpas_edge_plane_slopes_ring: optional validMask is too small.', MPAS_LOG_CRIT) + end if + end if + + normalSlopeEdge = 0.0_RKIND + tangentSlopeEdge = 0.0_RKIND + + allocate(visited(nCells)) + allocate(ringByCell(nCells)) + allocate(frontier(nCells)) + allocate(nextFrontier(nCells)) + allocate(stencil(nCells)) + + visited = 0 + ringByCell = -1 + visitTag = 0 + + do iEdge = 1, nEdgesSolve + + iCell1 = cellsOnEdge(1, iEdge) + iCell2 = cellsOnEdge(2, iEdge) + + seedCells(1) = iCell1 + seedCells(2) = iCell2 + + call mpas_build_contiguous_cell_stencil(seedCells, 2, ringCount, nCells, nEdgesOnCell, cellsOnCell, & + visited, ringByCell, frontier, nextFrontier, stencil, visitTag, & + stencilCount, maxBuiltRing, validMask) + + if (stencilCount <= 0) cycle + + xRef = xEdge(iEdge) + yRef = yEdge(iEdge) + zRef = fieldCell(stencil(1)) + + effectiveRing = maxBuiltRing + fitOK = .false. + + do while (effectiveRing >= 1 .and. .not. fitOK) + call mpas_try_plane_fit_reference(xRef, yRef, zRef, effectiveRing, stencil, stencilCount, ringByCell, & + xCell, yCell, fieldCell, aCoeff, bCoeff, fitOK) + if (.not. fitOK) effectiveRing = effectiveRing - 1 + end do + + if (.not. fitOK) cycle + + edgeCos = cos(angleEdge(iEdge)) + edgeSin = sin(angleEdge(iEdge)) + + normalSlopeEdge(iEdge) = aCoeff*edgeCos + bCoeff*edgeSin + tangentSlopeEdge(iEdge) = -aCoeff*edgeSin + bCoeff*edgeCos + + end do + + deallocate(visited) + deallocate(ringByCell) + deallocate(frontier) + deallocate(nextFrontier) + deallocate(stencil) + + end subroutine mpas_edge_plane_slopes_ring!}}} + + !*********************************************************************** ! ! routine mpas_build_contiguous_cell_stencil