CliMA · glwagner · Sep 12, 2025 · Apr 11, 2025 · Apr 11, 2025 · Apr 11, 2025
diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml
@@ -71,7 +71,6 @@ steps:
           - "👻 abstract_operations"
           - "🕊 poisson_solvers_1"
           - "🦖 poisson_solvers_2"
-          - "🌷 matrix_poisson_solvers"
           - "🦤 general_solvers"
           - "🎣 turbulence_closures"
           - "🦀 time_stepping_1"

diff --git a/Project.toml b/Project.toml
@@ -15,7 +15,6 @@ GPUArrays = "0c68f7d7-f131-5f86-a1c3-88cf8149b2d7"
 GPUArraysCore = "46192b85-c4d5-4398-a991-12ede77f4527"
 Glob = "c27321d9-0574-5035-807b-f59d2c89b15c"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
-IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
 Krylov = "ba0b0d4f-ebba-5204-a429-3ac8c609bfb7"
@@ -75,7 +74,6 @@ GPUArrays = "10.3, 11.2.2"
 GPUArraysCore = "0.2"
 Glob = "1.3"
 InteractiveUtils = "1.9"
-IterativeSolvers = "0.9"
 JLD2 = "0.4, 0.5, 0.6"
 KernelAbstractions = "0.9.21"
 Krylov = "0.10.0"

diff --git a/benchmark/benchmark_implicit_surface_solvers.jl b/benchmark/benchmark_implicit_surface_solvers.jl
@@ -67,8 +67,6 @@ for N in 10:10:250
 
     implicit_free_surface_solvers = (:FastFourierTransform,
                                      :PreconditionedConjugateGradient,
-                                     :HeptadiagonalIterativeSolver,
-                                     :HeptadiagonalIterativeSolver_withMGpreconditioner,
                                      :PreconditionedConjugateGradient_withFFTpreconditioner,
                                     )
 
@@ -87,8 +85,6 @@ for N in 10:10:250
         if implicit_free_surface_solver == :PreconditionedConjugateGradient_withFFTpreconditioner
             fft_preconditioner = FFTImplicitFreeSurfaceSolver(grid)
             free_surface = ImplicitFreeSurface(solver_method=:PreconditionedConjugateGradient, preconditioner=fft_preconditioner, reltol=sqrt(eps(eltype(grid))), abstol=0)
-        elseif implicit_free_surface_solver == :HeptadiagonalIterativeSolver
-            free_surface = ImplicitFreeSurface(solver_method=implicit_free_surface_solver, tolerance=sqrt(eps(eltype(grid))))
         else
             free_surface = ImplicitFreeSurface(solver_method=implicit_free_surface_solver, reltol=sqrt(eps(eltype(grid))), abstol=0)
         end

diff --git a/benchmark/benchmark_spai_preconditioner.jl b/benchmark/benchmark_spai_preconditioner.jl
diff --git a/ext/OceananigansAMDGPUExt.jl b/ext/OceananigansAMDGPUExt.jl
@@ -67,15 +67,6 @@ end
 
 @inline AC.unsafe_free!(a::ROCArray) = AMDGPU.unsafe_free!(a)
 
-@inline AC.constructors(::AC.GPU{ROCBackend}, A::SparseMatrixCSC) = (ROCArray(A.colptr), ROCArray(A.rowval), ROCArray(A.nzval),  (A.m, A.n))
-@inline AC.constructors(::AC.CPU, A::ROCSparseMatrixCSC) = (A.dims[1], A.dims[2], Int64.(Array(A.colPtr)), Int64.(Array(A.rowVal)), Array(A.nzVal))
-@inline AC.constructors(::AC.GPU{ROCBackend}, A::ROCSparseMatrixCSC) = (A.colPtr, A.rowVal, A.nzVal,  A.dims)
-
-@inline AC.arch_sparse_matrix(::AC.GPU{ROCBackend}, constr::Tuple) = ROCSparseMatrixCSC(constr...)
-@inline AC.arch_sparse_matrix(::AC.CPU, A::ROCSparseMatrixCSC)   = SparseMatrixCSC(AC.constructors(AC.CPU(), A)...)
-@inline AC.arch_sparse_matrix(::AC.GPU{ROCBackend}, A::SparseMatrixCSC)     = ROCSparseMatrixCSC(AC.constructors(AC.GPU(), A)...)
-@inline AC.arch_sparse_matrix(::AC.GPU{ROCBackend}, A::ROCSparseMatrixCSC) = A
-
 @inline convert_to_device(::ROCGPU, args) = AMDGPU.rocconvert(args)
 @inline convert_to_device(::ROCGPU, args::Tuple) = map(AMDGPU.rocconvert, args)
 

diff --git a/ext/OceananigansCUDAExt.jl b/ext/OceananigansCUDAExt.jl
@@ -89,15 +89,6 @@ end
 
 @inline AC.unsafe_free!(a::CuArray) = CUDA.unsafe_free!(a)
 
-@inline AC.constructors(::AC.GPU{CUDABackend}, A::SparseMatrixCSC) = (CuArray(A.colptr), CuArray(A.rowval), CuArray(A.nzval),  (A.m, A.n))
-@inline AC.constructors(::AC.CPU, A::CuSparseMatrixCSC) = (A.dims[1], A.dims[2], Int64.(Array(A.colPtr)), Int64.(Array(A.rowVal)), Array(A.nzVal))
-@inline AC.constructors(::AC.GPU{CUDABackend}, A::CuSparseMatrixCSC) = (A.colPtr, A.rowVal, A.nzVal,  A.dims)
-
-@inline AC.arch_sparse_matrix(::AC.GPU{CUDABackend}, constr::Tuple) = CuSparseMatrixCSC(constr...)
-@inline AC.arch_sparse_matrix(::AC.CPU, A::CuSparseMatrixCSC)   = SparseMatrixCSC(AC.constructors(AC.CPU(), A)...)
-@inline AC.arch_sparse_matrix(::AC.GPU{CUDABackend}, A::SparseMatrixCSC)     = CuSparseMatrixCSC(AC.constructors(AC.GPU(), A)...)
-@inline AC.arch_sparse_matrix(::AC.GPU{CUDABackend}, A::CuSparseMatrixCSC) = A
-
 @inline AC.convert_to_device(::CUDAGPU, args) = CUDA.cudaconvert(args)
 @inline AC.convert_to_device(::CUDAGPU, args::Tuple) = map(CUDA.cudaconvert, args)
 

diff --git a/src/Architectures.jl b/src/Architectures.jl
@@ -4,8 +4,7 @@ export AbstractArchitecture, AbstractSerialArchitecture
 export CPU, GPU, ReactantState
 export device, device!, ndevices, synchronize, architecture, unified_array, device_copy_to!
 export array_type, on_architecture
-export constructors, unpack_constructors, copy_unpack_constructors
-export arch_sparse_matrix, child_architecture
+export child_architecture
 
 using Adapt
 using OffsetArrays
@@ -117,18 +116,4 @@ unified_array(::GPU, a) = a
 @inline convert_to_device(arch, args)  = args
 @inline convert_to_device(::CPU, args) = args
 
-# Utils for sparse matrix manipulation
-@inline constructors(::CPU, A::SparseMatrixCSC) = (A.m, A.n, A.colptr, A.rowval, A.nzval)
-@inline constructors(::CPU, m::Number, n::Number, constr::Tuple) = (m, n, constr...)
-@inline constructors(::GPU, m::Number, n::Number, constr::Tuple) = (constr..., (m, n))
-
-@inline unpack_constructors(::CPU, constr::Tuple) = (constr[3], constr[4], constr[5])
-@inline unpack_constructors(::GPU, constr::Tuple) = (constr[1], constr[2], constr[3])
-
-@inline copy_unpack_constructors(::CPU, constr::Tuple) = deepcopy((constr[3], constr[4], constr[5]))
-@inline copy_unpack_constructors(::GPU, constr::Tuple) = deepcopy((constr[1], constr[2], constr[3]))
-
-@inline arch_sparse_matrix(::CPU, constr::Tuple) = SparseMatrixCSC(constr...)
-@inline arch_sparse_matrix(::CPU, A::SparseMatrixCSC) = A
-
 end # module
diff --git a/src/Models/HydrostaticFreeSurfaceModels/HydrostaticFreeSurfaceModels.jl b/src/Models/HydrostaticFreeSurfaceModels/HydrostaticFreeSurfaceModels.jl
@@ -64,7 +64,6 @@ include("implicit_free_surface_utils.jl")
 include("compute_vertically_integrated_variables.jl")
 include("fft_based_implicit_free_surface_solver.jl")
 include("pcg_implicit_free_surface_solver.jl")
-include("matrix_implicit_free_surface_solver.jl")
 include("implicit_free_surface.jl")
 include("hydrostatic_free_surface_field_tuples.jl")
 

diff --git a/src/Models/HydrostaticFreeSurfaceModels/implicit_free_surface.jl b/src/Models/HydrostaticFreeSurfaceModels/implicit_free_surface.jl
@@ -68,11 +68,10 @@ surface can be obtained using the [`FFTBasedPoissonSolver`](@ref).
 
 `solver_method` can be either of:
 * `:FastFourierTransform` for [`FFTBasedPoissonSolver`](@ref)
-* `:HeptadiagonalIterativeSolver`  for [`HeptadiagonalIterativeSolver`](@ref)
 * `:PreconditionedConjugateGradient` for [`ConjugateGradientSolver`](@ref)
 
 By default, if the grid has regular spacing in the horizontal directions then the `:FastFourierTransform` is chosen,
-otherwise the `:HeptadiagonalIterativeSolver`.
+otherwise the `:PreconditionedConjugateGradient`.
 """
 ImplicitFreeSurface(; solver_method=:Default, gravitational_acceleration=g_Earth, solver_settings...) =
     ImplicitFreeSurface(nothing, gravitational_acceleration, nothing, nothing, solver_method, solver_settings)
@@ -115,7 +114,7 @@ build_implicit_step_solver(::Val{:Default}, grid::XYRegularRG, settings, gravita
     build_implicit_step_solver(Val(:FastFourierTransform), grid, settings, gravitational_acceleration)
 
 build_implicit_step_solver(::Val{:Default}, grid, settings, gravitational_acceleration) =
-    build_implicit_step_solver(Val(:HeptadiagonalIterativeSolver), grid, settings, gravitational_acceleration)
+    build_implicit_step_solver(Val(:PreconditionedConjugateGradient), grid, settings, gravitational_acceleration)
 
 @inline explicit_barotropic_pressure_x_gradient(i, j, k, grid, ::ImplicitFreeSurface) = 0
 @inline explicit_barotropic_pressure_y_gradient(i, j, k, grid, ::ImplicitFreeSurface) = 0

diff --git a/src/Models/HydrostaticFreeSurfaceModels/matrix_implicit_free_surface_solver.jl b/src/Models/HydrostaticFreeSurfaceModels/matrix_implicit_free_surface_solver.jl
diff --git a/src/MultiRegion/MultiRegion.jl b/src/MultiRegion/MultiRegion.jl
@@ -75,7 +75,6 @@ include("multi_region_field.jl")
 include("multi_region_abstract_operations.jl")
 include("multi_region_boundary_conditions.jl")
 include("multi_region_reductions.jl")
-include("unified_implicit_free_surface_solver.jl")
 include("multi_region_split_explicit_free_surface.jl")
 include("multi_region_models.jl")
 include("multi_region_output_writers.jl")