diff --git a/NEWS.md b/NEWS.md
index 2b877ef64c..85c3f4d808 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -4,6 +4,16 @@ ClimaCore.jl Release Notes
 main
 -------
 
+v0.14.5
+-------
+
+- ![][badge-🐛bugfix] `run_field_matrix_solver!` was fixed for column spaces, and tests were added to ensure it doesn't break in the future.
+  PR [#1750](https://github.com/CliMA/ClimaCore.jl/pull/1750)
+- ![][badge-🚀performance] We're now using local memory (MArrays) in the `band_matrix_solve!`, which has improved performance. PR [#1735](https://github.com/CliMA/ClimaCore.jl/pull/1735).
+- ![][badge-🚀performance] We've specialized some cases in `run_field_matrix_solver!`, which results in more efficient kernels being launched. PR [#1732](https://github.com/CliMA/ClimaCore.jl/pull/1732).
+- ![][badge-🚀performance] We've reduced memory reads in the `band_matrix_solve!` for tridiagonal systems, improving its performance. PR [#1731](https://github.com/CliMA/ClimaCore.jl/pull/1731).
+- ![][badge-🚀performance] We've added NVTX annotations in ClimaCore functions, so that we have a more granular trace of performance. PRs [#1726](https://github.com/CliMA/ClimaCore.jl/pull/1726), [#1723](https://github.com/CliMA/ClimaCore.jl/pull/1723).
+
 v0.14.0
 -------
 
diff --git a/Project.toml b/Project.toml
index 55793f7af0..71ba5458da 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "ClimaCore"
 uuid = "d414da3d-4745-48bb-8d80-42e94e092884"
 authors = ["CliMA Contributors <clima-software@caltech.edu>"]
-version = "0.14.4"
+version = "0.14.5"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
diff --git a/src/MatrixFields/single_field_solver.jl b/src/MatrixFields/single_field_solver.jl
index 77435213a6..bbef957895 100644
--- a/src/MatrixFields/single_field_solver.jl
+++ b/src/MatrixFields/single_field_solver.jl
@@ -70,16 +70,28 @@ single_field_solve!(::ClimaComms.AbstractCPUDevice, cache, x, A, b) =
     _single_field_solve!(ClimaComms.device(axes(A)), cache, x, A, b)
 
 # CPU (GPU has already called Spaces.column on arg)
-_single_field_solve!(device::ClimaComms.AbstractCPUDevice, cache, x, A, b) =
-    Fields.bycolumn(axes(A)) do colidx
-        _single_field_solve_col!(
-            ClimaComms.device(axes(A)),
-            cache[colidx],
-            x[colidx],
-            A[colidx],
-            b[colidx],
-        )
+function _single_field_solve!(
+    device::ClimaComms.AbstractCPUDevice,
+    cache,
+    x,
+    A,
+    b,
+)
+    space = axes(x)
+    if space isa Spaces.FiniteDifferenceSpace
+        _single_field_solve_col!(device, cache, x, A, b)
+    else
+        Fields.bycolumn(space) do colidx
+            _single_field_solve_col!(
+                device,
+                cache[colidx],
+                x[colidx],
+                A[colidx],
+                b[colidx],
+            )
+        end
     end
+end
 
 function _single_field_solve_col!(
     ::ClimaComms.AbstractCPUDevice,
diff --git a/test/MatrixFields/field_matrix_solvers.jl b/test/MatrixFields/field_matrix_solvers.jl
index 0b7026546c..1a75db1b3f 100644
--- a/test/MatrixFields/field_matrix_solvers.jl
+++ b/test/MatrixFields/field_matrix_solvers.jl
@@ -9,6 +9,8 @@ import ClimaComms
 import ClimaCore.Utilities: half
 import ClimaCore.RecursiveApply: ⊠
 import ClimaCore.MatrixFields: @name
+import ClimaCore:
+    Spaces, MatrixFields, Fields, Domains, Meshes, Topologies, Geometry
 
 include("matrix_field_test_utils.jl")
 
@@ -79,462 +81,495 @@ function test_field_matrix_solver(; test_name, alg, A, b, use_rel_error = false)
     end
 end
 
-@testset "FieldMatrixSolver Unit Tests" begin
-    FT = Float64
-    center_space, face_space = test_spaces(FT)
-    surface_space = Spaces.level(face_space, half)
-
-    seed!(1) # ensures reproducibility
-
-    ᶜvec = random_field(FT, center_space)
-    ᶠvec = random_field(FT, face_space)
-    sfc_vec = random_field(FT, surface_space)
-
-    # Make each random square matrix diagonally dominant in order to avoid large
-    # large roundoff errors when computing its inverse. Scale the non-square
-    # matrices by the same amount as the square matrices.
-    λ = 10 # scale factor
-    ᶜᶜmat1 = random_field(DiagonalMatrixRow{FT}, center_space) ./ λ .+ (I,)
-    ᶠᶠmat1 = random_field(DiagonalMatrixRow{FT}, face_space) ./ λ .+ (I,)
-    ᶜᶠmat2 = random_field(BidiagonalMatrixRow{FT}, center_space) ./ λ
-    ᶠᶜmat2 = random_field(BidiagonalMatrixRow{FT}, face_space) ./ λ
-    ᶜᶜmat3 = random_field(TridiagonalMatrixRow{FT}, center_space) ./ λ .+ (I,)
-    ᶠᶠmat3 = random_field(TridiagonalMatrixRow{FT}, face_space) ./ λ .+ (I,)
-    ᶜᶠmat4 = random_field(QuaddiagonalMatrixRow{FT}, center_space) ./ λ
-    ᶠᶜmat4 = random_field(QuaddiagonalMatrixRow{FT}, face_space) ./ λ
-    ᶜᶜmat5 = random_field(PentadiagonalMatrixRow{FT}, center_space) ./ λ .+ (I,)
-    ᶠᶠmat5 = random_field(PentadiagonalMatrixRow{FT}, face_space) ./ λ .+ (I,)
-
-    for (vector, matrix, string1, string2) in (
-        (sfc_vec, I, "UniformScaling", "a single level"),
-        (ᶜvec, I, "UniformScaling", "cell centers"),
-        (ᶠvec, I, "UniformScaling", "cell faces"),
-        (ᶜvec, ᶜᶜmat1, "diagonal matrix", "cell centers"),
-        (ᶠvec, ᶠᶠmat1, "diagonal matrix", "cell faces"),
-        (ᶜvec, ᶜᶜmat3, "tri-diagonal matrix", "cell centers"),
-        (ᶠvec, ᶠᶠmat3, "tri-diagonal matrix", "cell faces"),
-        (ᶜvec, ᶜᶜmat5, "penta-diagonal matrix", "cell centers"),
-        (ᶠvec, ᶠᶠmat5, "penta-diagonal matrix", "cell faces"),
-    )
-        test_field_matrix_solver(;
-            test_name = "$string1 solve on $string2",
-            alg = MatrixFields.BlockDiagonalSolve(),
-            A = MatrixFields.FieldMatrix((@name(_), @name(_)) => matrix),
-            b = Fields.FieldVector(; _ = vector),
-        )
-    end
-
-    # TODO: Add a simple test where typeof(x) != typeof(b).
-
-    # Note: The round-off error of StationaryIterativeSolve can be much larger
-    # on GPUs, so n_iters often has to be increased when using_cuda is true.
-
-    for alg in (
-        MatrixFields.BlockDiagonalSolve(),
-        MatrixFields.BlockLowerTriangularSolve(@name(c)),
-        MatrixFields.BlockArrowheadSolve(@name(c)),
-        MatrixFields.ApproximateBlockArrowheadIterativeSolve(@name(c)),
-        MatrixFields.StationaryIterativeSolve(; n_iters = using_cuda ? 28 : 18),
-    )
-        test_field_matrix_solver(;
-            test_name = "$(typeof(alg).name.name) for a block diagonal matrix \
-                         with diagonal and penta-diagonal blocks",
-            alg,
-            A = MatrixFields.FieldMatrix(
-                (@name(c), @name(c)) => ᶜᶜmat1,
-                (@name(f), @name(f)) => ᶠᶠmat5,
-            ),
-            b = Fields.FieldVector(; c = ᶜvec, f = ᶠvec),
-        )
-    end
-
-    test_field_matrix_solver(;
-        test_name = "BlockDiagonalSolve for a block diagonal matrix with \
-                     tri-diagonal and penta-diagonal blocks",
-        alg = MatrixFields.BlockDiagonalSolve(),
-        A = MatrixFields.FieldMatrix(
-            (@name(c), @name(c)) => ᶜᶜmat3,
-            (@name(f), @name(f)) => ᶠᶠmat5,
-        ),
-        b = Fields.FieldVector(; c = ᶜvec, f = ᶠvec),
-    )
-
-    test_field_matrix_solver(;
-        test_name = "BlockLowerTriangularSolve for a block lower triangular \
-                     matrix with tri-diagonal, bi-diagonal, and penta-diagonal \
-                     blocks",
-        alg = MatrixFields.BlockLowerTriangularSolve(@name(c)),
-        A = MatrixFields.FieldMatrix(
-            (@name(c), @name(c)) => ᶜᶜmat3,
-            (@name(f), @name(c)) => ᶠᶜmat2,
-            (@name(f), @name(f)) => ᶠᶠmat5,
-        ),
-        b = Fields.FieldVector(; c = ᶜvec, f = ᶠvec),
-    )
-
-    test_field_matrix_solver(;
-        test_name = "BlockArrowheadSolve for a block matrix with diagonal, \
-                     quad-diagonal, bi-diagonal, and penta-diagonal blocks",
-        alg = MatrixFields.BlockArrowheadSolve(@name(c)),
-        A = MatrixFields.FieldMatrix(
-            (@name(c), @name(c)) => ᶜᶜmat1,
-            (@name(c), @name(f)) => ᶜᶠmat4,
-            (@name(f), @name(c)) => ᶠᶜmat2,
-            (@name(f), @name(f)) => ᶠᶠmat5,
-        ),
-        b = Fields.FieldVector(; c = ᶜvec, f = ᶠvec),
-    )
-
-    # Since test_field_matrix_solver runs the solver many times with the same
-    # values of x, A, and b for benchmarking, setting correlated_solves to true
-    # is equivalent to setting n_iters to some very large number.
-    test_field_matrix_solver(;
-        test_name = "StationaryIterativeSolve with correlated_solves for a \
-                     block matrix with tri-diagonal, quad-diagonal, \
-                     bi-diagonal, and penta-diagonal blocks",
-        alg = MatrixFields.StationaryIterativeSolve(; correlated_solves = true),
-        A = MatrixFields.FieldMatrix(
-            (@name(c), @name(c)) => ᶜᶜmat3,
-            (@name(c), @name(f)) => ᶜᶠmat4,
-            (@name(f), @name(c)) => ᶠᶜmat2,
-            (@name(f), @name(f)) => ᶠᶠmat5,
-        ),
-        b = Fields.FieldVector(; c = ᶜvec, f = ᶠvec),
-    )
-
-    # Each of the scaled identity matrices below was chosen to minimize the
-    # value of ρ(I - P⁻¹ * A), which was found by setting print_radius to true.
-    # Each value of n_iters below was then chosen to be the smallest value for
-    # which the relative error was less than 1e-6.
-    scaled_identity_matrix(scalar) =
-        MatrixFields.FieldMatrix((@name(), @name()) => scalar * I)
-    for (P_name, alg) in (
-        (
-            "no (identity matrix)",
-            MatrixFields.StationaryIterativeSolve(;
-                n_iters = using_cuda ? 10 : 7,
-            ),
-        ), # ρ(I - P⁻¹ * A) ≈ 0.3777
-        (
-            "Richardson (damped identity matrix)",
-            MatrixFields.StationaryIterativeSolve(;
-                P_alg = MatrixFields.CustomPreconditioner(
-                    scaled_identity_matrix(FT(1.12)),
-                ),
-                n_iters = using_cuda ? 8 : 7,
-            ),
-        ), # ρ(I - P⁻¹ * A) ≈ 0.2294
-        (
-            "Jacobi (diagonal)",
-            MatrixFields.StationaryIterativeSolve(;
-                P_alg = MatrixFields.MainDiagonalPreconditioner(),
-                n_iters = using_cuda ? 8 : 6,
-            ),
-        ), # ρ(I - P⁻¹ * A) ≈ 0.3241
-        (
-            "damped Jacobi (diagonal)",
-            MatrixFields.StationaryIterativeSolve(;
-                P_alg = MatrixFields.WeightedPreconditioner(
-                    scaled_identity_matrix(FT(1.08)),
-                    MatrixFields.MainDiagonalPreconditioner(),
-                ),
-                n_iters = using_cuda ? 8 : 7,
-            ),
-        ), # ρ(I - P⁻¹ * A) ≈ 0.2249
-        (
-            "block Jacobi (diagonal)",
-            MatrixFields.StationaryIterativeSolve(;
-                P_alg = MatrixFields.BlockDiagonalPreconditioner(),
-                n_iters = 7,
-            ),
-        ), # ρ(I - P⁻¹ * A) ≈ 0.1450
-        (
-            "damped block Jacobi (diagonal)",
-            MatrixFields.StationaryIterativeSolve(;
-                P_alg = MatrixFields.WeightedPreconditioner(
-                    scaled_identity_matrix(FT(1.002)),
-                    MatrixFields.BlockDiagonalPreconditioner(),
-                ),
-                n_iters = 7,
-            ),
-        ), # ρ(I - P⁻¹ * A) ≈ 0.1427
-        (
-            "block arrowhead",
-            MatrixFields.StationaryIterativeSolve(;
-                P_alg = MatrixFields.BlockArrowheadPreconditioner(@name(c)),
-                n_iters = 6,
-            ),
-        ), # ρ(I - P⁻¹ * A) ≈ 0.1356
-        (
-            "damped block arrowhead",
-            MatrixFields.StationaryIterativeSolve(;
-                P_alg = MatrixFields.BlockArrowheadPreconditioner(
-                    @name(c);
-                    P_alg₁ = MatrixFields.WeightedPreconditioner(
-                        scaled_identity_matrix(FT(1.0001)),
-                        MatrixFields.MainDiagonalPreconditioner(),
-                    ),
-                ),
-                n_iters = 6,
-            ),
-        ), # ρ(I - P⁻¹ * A) ≈ 0.1355
-        (
-            "block arrowhead Schur complement",
-            MatrixFields.ApproximateBlockArrowheadIterativeSolve(
-                @name(c);
-                n_iters = 3,
-            ),
-        ), # ρ(I - P⁻¹ * A) ≈ 0.0009
-        (
-            "damped block arrowhead Schur complement",
-            MatrixFields.ApproximateBlockArrowheadIterativeSolve(
-                @name(c);
-                P_alg₁ = MatrixFields.WeightedPreconditioner(
-                    scaled_identity_matrix(FT(1.09)),
-                    MatrixFields.MainDiagonalPreconditioner(),
-                ),
-                n_iters = 2,
-            ),
-        ), # ρ(I - P⁻¹ * A) ≈ 0.000006
-    )
-        test_field_matrix_solver(;
-            test_name = "approximate iterative solve with $P_name \
-                         preconditioning for a block matrix with tri-diagonal, \
-                         quad-diagonal, bi-diagonal, and penta-diagonal blocks",
-            alg,
-            A = MatrixFields.FieldMatrix(
-                (@name(c), @name(c)) => ᶜᶜmat3,
-                (@name(c), @name(f)) => ᶜᶠmat4,
-                (@name(f), @name(c)) => ᶠᶜmat2,
-                (@name(f), @name(f)) => ᶠᶠmat5,
-            ),
-            b = Fields.FieldVector(; c = ᶜvec, f = ᶠvec),
-            use_rel_error = true,
-        )
-    end
-
-    @testset "approximate iterative solve with debugging" begin
-        Logging.with_logger(Logging.SimpleLogger(stderr, Logging.Debug)) do
-            # Recreate the setup from the previous unit test.
-            alg = MatrixFields.ApproximateBlockArrowheadIterativeSolve(
-                @name(c);
-                P_alg₁ = MatrixFields.WeightedPreconditioner(
-                    scaled_identity_matrix(FT(1.09)),
-                    MatrixFields.MainDiagonalPreconditioner(),
-                ),
-                n_iters = 2,
-            )
-            A = MatrixFields.FieldMatrix(
-                (@name(c), @name(c)) => ᶜᶜmat3,
-                (@name(c), @name(f)) => ᶜᶠmat4,
-                (@name(f), @name(c)) => ᶠᶜmat2,
-                (@name(f), @name(f)) => ᶠᶠmat5,
-            )
-            b = Fields.FieldVector(; c = ᶜvec, f = ᶠvec)
-
-            x = similar(b)
-            solver = FieldMatrixSolver(alg, A, b)
-            args = (solver, x, A, b)
-
-            # Compare the debugging logs to RegEx strings. Note that debugging the
-            # spectral radius is currently not possible on GPUs.
-            spectral_radius_logs =
-                using_cuda ? () : ((:debug, r"ρ\(I \- inv\(P\) \* A\) ≈"),)
-            error_norm_logs = (
-                (:debug, r"||x[0] - x'||₂ ≈"),
-                (:debug, r"||x[1] - x'||₂ ≈"),
-                (:debug, r"||x[2] - x'||₂ ≈"),
-            )
-            logs = (spectral_radius_logs..., error_norm_logs...)
-            @test_logs logs... min_level = Logging.Debug field_matrix_solve!(
-                args...,
-            )
-        end
-    end
-end
-
-@testset "FieldMatrixSolver ClimaAtmos-Based Tests" begin
-    FT = Float64
-    center_space, face_space = test_spaces(FT)
-    surface_space = Spaces.level(face_space, half)
-
-    seed!(1) # ensures reproducibility
-
-    ᶜvec = random_field(FT, center_space)
-    ᶠvec = random_field(FT, face_space)
-    sfc_vec = random_field(FT, surface_space)
-
-    # Make each random square matrix diagonally dominant in order to avoid large
-    # large roundoff errors when computing its inverse. Scale the non-square
-    # matrices by the same amount as the square matrices.
-    λ = 10 # scale factor
-    ᶜᶜmat1 = random_field(DiagonalMatrixRow{FT}, center_space) ./ λ .+ (I,)
-    ᶜᶠmat2 = random_field(BidiagonalMatrixRow{FT}, center_space) ./ λ
-    ᶠᶜmat2 = random_field(BidiagonalMatrixRow{FT}, face_space) ./ λ
-    ᶜᶜmat3 = random_field(TridiagonalMatrixRow{FT}, center_space) ./ λ .+ (I,)
-    ᶠᶠmat3 = random_field(TridiagonalMatrixRow{FT}, face_space) ./ λ .+ (I,)
-
-    e¹² = Geometry.Covariant12Vector(1, 1)
-    e³ = Geometry.Covariant3Vector(1)
-    e₃ = Geometry.Contravariant3Vector(1)
-
-    ρχ_unit = (; ρq_tot = 1, ρq_liq = 1, ρq_ice = 1, ρq_rai = 1, ρq_sno = 1)
-    ρaχ_unit =
-        (; ρaq_tot = 1, ρaq_liq = 1, ρaq_ice = 1, ρaq_rai = 1, ρaq_sno = 1)
-
-    dry_center_gs_unit = (; ρ = 1, ρe_tot = 1, uₕ = e¹²)
-    center_gs_unit = (; dry_center_gs_unit..., ρatke = 1, ρχ = ρχ_unit)
-    center_sgsʲ_unit = (; ρa = 1, ρae_tot = 1, ρaχ = ρaχ_unit)
-
-    ᶜᶜmat3_uₕ_scalar = ᶜᶜmat3 .* (e¹²,)
-    ᶠᶜmat2_u₃_scalar = ᶠᶜmat2 .* (e³,)
-    ᶜᶠmat2_scalar_u₃ = ᶜᶠmat2 .* (e₃',)
-    ᶜᶠmat2_uₕ_u₃ = ᶜᶠmat2 .* (e¹² * e₃',)
-    ᶠᶠmat3_u₃_u₃ = ᶠᶠmat3 .* (e³ * e₃',)
-    ᶜᶜmat3_ρχ_scalar = map(Base.Fix1(map, Base.Fix2(⊠, ρχ_unit)), ᶜᶜmat3)
-    ᶜᶜmat3_ρaχ_scalar = map(Base.Fix1(map, Base.Fix2(⊠, ρaχ_unit)), ᶜᶜmat3)
-    ᶜᶠmat2_ρχ_u₃ = map(Base.Fix1(map, Base.Fix2(⊠, ρχ_unit ⊠ e₃')), ᶜᶠmat2)
-    ᶜᶠmat2_ρaχ_u₃ = map(Base.Fix1(map, Base.Fix2(⊠, ρaχ_unit ⊠ e₃')), ᶜᶠmat2)
-    # We need to use Fix1 and Fix2 instead of defining anonymous functions in
-    # order for the result of map to be inferrable.
-
-    b_dry_dycore = Fields.FieldVector(;
-        c = ᶜvec .* (dry_center_gs_unit,),
-        f = ᶠvec .* ((; u₃ = e³),),
-    )
-
-    b_moist_dycore_diagnostic_edmf = Fields.FieldVector(;
-        c = ᶜvec .* (center_gs_unit,),
-        f = ᶠvec .* ((; u₃ = e³),),
-    )
-
-    b_moist_dycore_prognostic_edmf_prognostic_surface = Fields.FieldVector(;
-        sfc = sfc_vec .* ((; T = 1),),
-        c = ᶜvec .* ((; center_gs_unit..., sgsʲs = (center_sgsʲ_unit,)),),
-        f = ᶠvec .* ((; u₃ = e³, sgsʲs = ((; u₃ = e³),)),),
-    )
-
-    test_field_matrix_solver(;
-        test_name = "similar solve to ClimaAtmos's dry dycore with implicit \
-                     acoustic waves",
-        alg = MatrixFields.BlockArrowheadSolve(@name(c)),
-        A = MatrixFields.FieldMatrix(
-            (@name(c.ρ), @name(c.ρ)) => I,
-            (@name(c.ρe_tot), @name(c.ρe_tot)) => I,
-            (@name(c.uₕ), @name(c.uₕ)) => I,
-            (@name(c.ρ), @name(f.u₃)) => ᶜᶠmat2_scalar_u₃,
-            (@name(c.ρe_tot), @name(f.u₃)) => ᶜᶠmat2_scalar_u₃,
-            (@name(f.u₃), @name(c.ρ)) => ᶠᶜmat2_u₃_scalar,
-            (@name(f.u₃), @name(c.ρe_tot)) => ᶠᶜmat2_u₃_scalar,
-            (@name(f.u₃), @name(f.u₃)) => ᶠᶠmat3_u₃_u₃,
-        ),
-        b = b_dry_dycore,
-    )
-
-    test_field_matrix_solver(;
-        test_name = "similar solve to ClimaAtmos's dry dycore with implicit \
-                     acoustic waves and diffusion",
-        alg = MatrixFields.ApproximateBlockArrowheadIterativeSolve(
-            @name(c);
-            n_iters = 6,
-        ),
-        A = MatrixFields.FieldMatrix(
-            (@name(c.ρ), @name(c.ρ)) => I,
-            (@name(c.ρe_tot), @name(c.ρe_tot)) => ᶜᶜmat3,
-            (@name(c.uₕ), @name(c.uₕ)) => ᶜᶜmat3,
-            (@name(c.ρ), @name(f.u₃)) => ᶜᶠmat2_scalar_u₃,
-            (@name(c.ρe_tot), @name(f.u₃)) => ᶜᶠmat2_scalar_u₃,
-            (@name(f.u₃), @name(c.ρ)) => ᶠᶜmat2_u₃_scalar,
-            (@name(f.u₃), @name(c.ρe_tot)) => ᶠᶜmat2_u₃_scalar,
-            (@name(f.u₃), @name(f.u₃)) => ᶠᶠmat3_u₃_u₃,
-        ),
-        b = b_dry_dycore,
-    )
-
-    test_field_matrix_solver(;
-        test_name = "similar solve to ClimaAtmos's moist dycore + diagnostic \
-                     EDMF with implicit acoustic waves and SGS fluxes",
-        alg = MatrixFields.ApproximateBlockArrowheadIterativeSolve(
-            @name(c);
-            n_iters = 6,
-        ),
-        A = MatrixFields.FieldMatrix(
-            (@name(c.ρ), @name(c.ρ)) => I,
-            (@name(c.ρe_tot), @name(c.ρe_tot)) => ᶜᶜmat3,
-            (@name(c.ρatke), @name(c.ρatke)) => ᶜᶜmat3,
-            (@name(c.ρχ), @name(c.ρχ)) => ᶜᶜmat3,
-            (@name(c.uₕ), @name(c.uₕ)) => ᶜᶜmat3,
-            (@name(c.ρ), @name(f.u₃)) => ᶜᶠmat2_scalar_u₃,
-            (@name(c.ρe_tot), @name(f.u₃)) => ᶜᶠmat2_scalar_u₃,
-            (@name(c.ρatke), @name(f.u₃)) => ᶜᶠmat2_scalar_u₃,
-            (@name(c.ρχ), @name(f.u₃)) => ᶜᶠmat2_ρχ_u₃,
-            (@name(f.u₃), @name(c.ρ)) => ᶠᶜmat2_u₃_scalar,
-            (@name(f.u₃), @name(c.ρe_tot)) => ᶠᶜmat2_u₃_scalar,
-            (@name(f.u₃), @name(f.u₃)) => ᶠᶠmat3_u₃_u₃,
-        ),
-        b = b_moist_dycore_diagnostic_edmf,
-    )
-
-    test_field_matrix_solver(;
-        test_name = "similar solve to ClimaAtmos's moist dycore + prognostic \
-                     EDMF + prognostic surface temperature with implicit \
-                     acoustic waves and SGS fluxes",
-        alg = MatrixFields.BlockLowerTriangularSolve(
-            @name(c.sgsʲs),
-            @name(f.sgsʲs);
-            alg₁ = MatrixFields.BlockArrowheadSolve(@name(c)),
-            alg₂ = MatrixFields.ApproximateBlockArrowheadIterativeSolve(
-                @name(c);
-                n_iters = 6,
-            ),
-        ),
-        A = MatrixFields.FieldMatrix(
-            # GS-GS blocks:
-            (@name(sfc), @name(sfc)) => I,
-            (@name(c.ρ), @name(c.ρ)) => I,
-            (@name(c.ρe_tot), @name(c.ρe_tot)) => ᶜᶜmat3,
-            (@name(c.ρatke), @name(c.ρatke)) => ᶜᶜmat3,
-            (@name(c.ρχ), @name(c.ρχ)) => ᶜᶜmat3,
-            (@name(c.uₕ), @name(c.uₕ)) => ᶜᶜmat3,
-            (@name(c.ρ), @name(f.u₃)) => ᶜᶠmat2_scalar_u₃,
-            (@name(c.ρe_tot), @name(f.u₃)) => ᶜᶠmat2_scalar_u₃,
-            (@name(c.ρatke), @name(f.u₃)) => ᶜᶠmat2_scalar_u₃,
-            (@name(c.ρχ), @name(f.u₃)) => ᶜᶠmat2_ρχ_u₃,
-            (@name(f.u₃), @name(c.ρ)) => ᶠᶜmat2_u₃_scalar,
-            (@name(f.u₃), @name(c.ρe_tot)) => ᶠᶜmat2_u₃_scalar,
-            (@name(f.u₃), @name(f.u₃)) => ᶠᶠmat3_u₃_u₃,
-            # GS-SGS blocks:
-            (@name(c.ρe_tot), @name(c.sgsʲs.:(1).ρae_tot)) => ᶜᶜmat3,
-            (@name(c.ρχ.ρq_tot), @name(c.sgsʲs.:(1).ρaχ.ρaq_tot)) => ᶜᶜmat3,
-            (@name(c.ρχ.ρq_liq), @name(c.sgsʲs.:(1).ρaχ.ρaq_liq)) => ᶜᶜmat3,
-            (@name(c.ρχ.ρq_ice), @name(c.sgsʲs.:(1).ρaχ.ρaq_ice)) => ᶜᶜmat3,
-            (@name(c.ρχ.ρq_rai), @name(c.sgsʲs.:(1).ρaχ.ρaq_rai)) => ᶜᶜmat3,
-            (@name(c.ρχ.ρq_sno), @name(c.sgsʲs.:(1).ρaχ.ρaq_sno)) => ᶜᶜmat3,
-            (@name(c.ρe_tot), @name(c.sgsʲs.:(1).ρa)) => ᶜᶜmat3,
-            (@name(c.ρatke), @name(c.sgsʲs.:(1).ρa)) => ᶜᶜmat3,
-            (@name(c.ρχ), @name(c.sgsʲs.:(1).ρa)) => ᶜᶜmat3_ρχ_scalar,
-            (@name(c.uₕ), @name(c.sgsʲs.:(1).ρa)) => ᶜᶜmat3_uₕ_scalar,
-            (@name(c.ρe_tot), @name(f.sgsʲs.:(1).u₃)) => ᶜᶠmat2_scalar_u₃,
-            (@name(c.ρatke), @name(f.sgsʲs.:(1).u₃)) => ᶜᶠmat2_scalar_u₃,
-            (@name(c.ρχ), @name(f.sgsʲs.:(1).u₃)) => ᶜᶠmat2_ρχ_u₃,
-            (@name(c.uₕ), @name(f.sgsʲs.:(1).u₃)) => ᶜᶠmat2_uₕ_u₃,
-            (@name(f.u₃), @name(c.sgsʲs.:(1).ρa)) => ᶠᶜmat2_u₃_scalar,
-            (@name(f.u₃), @name(f.sgsʲs.:(1).u₃)) => ᶠᶠmat3_u₃_u₃,
-            # SGS-SGS blocks:
-            (@name(c.sgsʲs.:(1).ρa), @name(c.sgsʲs.:(1).ρa)) => I,
-            (@name(c.sgsʲs.:(1).ρae_tot), @name(c.sgsʲs.:(1).ρae_tot)) => I,
-            (@name(c.sgsʲs.:(1).ρaχ), @name(c.sgsʲs.:(1).ρaχ)) => I,
-            (@name(c.sgsʲs.:(1).ρa), @name(f.sgsʲs.:(1).u₃)) =>
-                ᶜᶠmat2_scalar_u₃,
-            (@name(c.sgsʲs.:(1).ρae_tot), @name(f.sgsʲs.:(1).u₃)) =>
-                ᶜᶠmat2_scalar_u₃,
-            (@name(c.sgsʲs.:(1).ρaχ), @name(f.sgsʲs.:(1).u₃)) => ᶜᶠmat2_ρaχ_u₃,
-            (@name(f.sgsʲs.:(1).u₃), @name(c.sgsʲs.:(1).ρa)) =>
-                ᶠᶜmat2_u₃_scalar,
-            (@name(f.sgsʲs.:(1).u₃), @name(c.sgsʲs.:(1).ρae_tot)) =>
-                ᶠᶜmat2_u₃_scalar,
-            (@name(f.sgsʲs.:(1).u₃), @name(f.sgsʲs.:(1).u₃)) => ᶠᶠmat3_u₃_u₃,
-        ),
-        b = b_moist_dycore_prognostic_edmf_prognostic_surface,
+# @testset "FieldMatrixSolver Unit Tests" begin
+#     FT = Float64
+#     center_space, face_space = test_spaces(FT)
+#     surface_space = Spaces.level(face_space, half)
+
+#     seed!(1) # ensures reproducibility
+
+#     ᶜvec = random_field(FT, center_space)
+#     ᶠvec = random_field(FT, face_space)
+#     sfc_vec = random_field(FT, surface_space)
+
+#     # Make each random square matrix diagonally dominant in order to avoid large
+#     # large roundoff errors when computing its inverse. Scale the non-square
+#     # matrices by the same amount as the square matrices.
+#     λ = 10 # scale factor
+#     ᶜᶜmat1 = random_field(DiagonalMatrixRow{FT}, center_space) ./ λ .+ (I,)
+#     ᶠᶠmat1 = random_field(DiagonalMatrixRow{FT}, face_space) ./ λ .+ (I,)
+#     ᶜᶠmat2 = random_field(BidiagonalMatrixRow{FT}, center_space) ./ λ
+#     ᶠᶜmat2 = random_field(BidiagonalMatrixRow{FT}, face_space) ./ λ
+#     ᶜᶜmat3 = random_field(TridiagonalMatrixRow{FT}, center_space) ./ λ .+ (I,)
+#     ᶠᶠmat3 = random_field(TridiagonalMatrixRow{FT}, face_space) ./ λ .+ (I,)
+#     ᶜᶠmat4 = random_field(QuaddiagonalMatrixRow{FT}, center_space) ./ λ
+#     ᶠᶜmat4 = random_field(QuaddiagonalMatrixRow{FT}, face_space) ./ λ
+#     ᶜᶜmat5 = random_field(PentadiagonalMatrixRow{FT}, center_space) ./ λ .+ (I,)
+#     ᶠᶠmat5 = random_field(PentadiagonalMatrixRow{FT}, face_space) ./ λ .+ (I,)
+
+#     for (vector, matrix, string1, string2) in (
+#         (sfc_vec, I, "UniformScaling", "a single level"),
+#         (ᶜvec, I, "UniformScaling", "cell centers"),
+#         (ᶠvec, I, "UniformScaling", "cell faces"),
+#         (ᶜvec, ᶜᶜmat1, "diagonal matrix", "cell centers"),
+#         (ᶠvec, ᶠᶠmat1, "diagonal matrix", "cell faces"),
+#         (ᶜvec, ᶜᶜmat3, "tri-diagonal matrix", "cell centers"),
+#         (ᶠvec, ᶠᶠmat3, "tri-diagonal matrix", "cell faces"),
+#         (ᶜvec, ᶜᶜmat5, "penta-diagonal matrix", "cell centers"),
+#         (ᶠvec, ᶠᶠmat5, "penta-diagonal matrix", "cell faces"),
+#     )
+#         test_field_matrix_solver(;
+#             test_name = "$string1 solve on $string2",
+#             alg = MatrixFields.BlockDiagonalSolve(),
+#             A = MatrixFields.FieldMatrix((@name(_), @name(_)) => matrix),
+#             b = Fields.FieldVector(; _ = vector),
+#         )
+#     end
+
+#     # TODO: Add a simple test where typeof(x) != typeof(b).
+
+#     # Note: The round-off error of StationaryIterativeSolve can be much larger
+#     # on GPUs, so n_iters often has to be increased when using_cuda is true.
+
+#     for alg in (
+#         MatrixFields.BlockDiagonalSolve(),
+#         MatrixFields.BlockLowerTriangularSolve(@name(c)),
+#         MatrixFields.BlockArrowheadSolve(@name(c)),
+#         MatrixFields.ApproximateBlockArrowheadIterativeSolve(@name(c)),
+#         MatrixFields.StationaryIterativeSolve(; n_iters = using_cuda ? 28 : 18),
+#     )
+#         test_field_matrix_solver(;
+#             test_name = "$(typeof(alg).name.name) for a block diagonal matrix \
+#                          with diagonal and penta-diagonal blocks",
+#             alg,
+#             A = MatrixFields.FieldMatrix(
+#                 (@name(c), @name(c)) => ᶜᶜmat1,
+#                 (@name(f), @name(f)) => ᶠᶠmat5,
+#             ),
+#             b = Fields.FieldVector(; c = ᶜvec, f = ᶠvec),
+#         )
+#     end
+
+#     test_field_matrix_solver(;
+#         test_name = "BlockDiagonalSolve for a block diagonal matrix with \
+#                      tri-diagonal and penta-diagonal blocks",
+#         alg = MatrixFields.BlockDiagonalSolve(),
+#         A = MatrixFields.FieldMatrix(
+#             (@name(c), @name(c)) => ᶜᶜmat3,
+#             (@name(f), @name(f)) => ᶠᶠmat5,
+#         ),
+#         b = Fields.FieldVector(; c = ᶜvec, f = ᶠvec),
+#     )
+
+#     test_field_matrix_solver(;
+#         test_name = "BlockLowerTriangularSolve for a block lower triangular \
+#                      matrix with tri-diagonal, bi-diagonal, and penta-diagonal \
+#                      blocks",
+#         alg = MatrixFields.BlockLowerTriangularSolve(@name(c)),
+#         A = MatrixFields.FieldMatrix(
+#             (@name(c), @name(c)) => ᶜᶜmat3,
+#             (@name(f), @name(c)) => ᶠᶜmat2,
+#             (@name(f), @name(f)) => ᶠᶠmat5,
+#         ),
+#         b = Fields.FieldVector(; c = ᶜvec, f = ᶠvec),
+#     )
+
+#     test_field_matrix_solver(;
+#         test_name = "BlockArrowheadSolve for a block matrix with diagonal, \
+#                      quad-diagonal, bi-diagonal, and penta-diagonal blocks",
+#         alg = MatrixFields.BlockArrowheadSolve(@name(c)),
+#         A = MatrixFields.FieldMatrix(
+#             (@name(c), @name(c)) => ᶜᶜmat1,
+#             (@name(c), @name(f)) => ᶜᶠmat4,
+#             (@name(f), @name(c)) => ᶠᶜmat2,
+#             (@name(f), @name(f)) => ᶠᶠmat5,
+#         ),
+#         b = Fields.FieldVector(; c = ᶜvec, f = ᶠvec),
+#     )
+
+#     # Since test_field_matrix_solver runs the solver many times with the same
+#     # values of x, A, and b for benchmarking, setting correlated_solves to true
+#     # is equivalent to setting n_iters to some very large number.
+#     test_field_matrix_solver(;
+#         test_name = "StationaryIterativeSolve with correlated_solves for a \
+#                      block matrix with tri-diagonal, quad-diagonal, \
+#                      bi-diagonal, and penta-diagonal blocks",
+#         alg = MatrixFields.StationaryIterativeSolve(; correlated_solves = true),
+#         A = MatrixFields.FieldMatrix(
+#             (@name(c), @name(c)) => ᶜᶜmat3,
+#             (@name(c), @name(f)) => ᶜᶠmat4,
+#             (@name(f), @name(c)) => ᶠᶜmat2,
+#             (@name(f), @name(f)) => ᶠᶠmat5,
+#         ),
+#         b = Fields.FieldVector(; c = ᶜvec, f = ᶠvec),
+#     )
+
+#     # Each of the scaled identity matrices below was chosen to minimize the
+#     # value of ρ(I - P⁻¹ * A), which was found by setting print_radius to true.
+#     # Each value of n_iters below was then chosen to be the smallest value for
+#     # which the relative error was less than 1e-6.
+#     scaled_identity_matrix(scalar) =
+#         MatrixFields.FieldMatrix((@name(), @name()) => scalar * I)
+#     for (P_name, alg) in (
+#         (
+#             "no (identity matrix)",
+#             MatrixFields.StationaryIterativeSolve(;
+#                 n_iters = using_cuda ? 10 : 7,
+#             ),
+#         ), # ρ(I - P⁻¹ * A) ≈ 0.3777
+#         (
+#             "Richardson (damped identity matrix)",
+#             MatrixFields.StationaryIterativeSolve(;
+#                 P_alg = MatrixFields.CustomPreconditioner(
+#                     scaled_identity_matrix(FT(1.12)),
+#                 ),
+#                 n_iters = using_cuda ? 8 : 7,
+#             ),
+#         ), # ρ(I - P⁻¹ * A) ≈ 0.2294
+#         (
+#             "Jacobi (diagonal)",
+#             MatrixFields.StationaryIterativeSolve(;
+#                 P_alg = MatrixFields.MainDiagonalPreconditioner(),
+#                 n_iters = using_cuda ? 8 : 6,
+#             ),
+#         ), # ρ(I - P⁻¹ * A) ≈ 0.3241
+#         (
+#             "damped Jacobi (diagonal)",
+#             MatrixFields.StationaryIterativeSolve(;
+#                 P_alg = MatrixFields.WeightedPreconditioner(
+#                     scaled_identity_matrix(FT(1.08)),
+#                     MatrixFields.MainDiagonalPreconditioner(),
+#                 ),
+#                 n_iters = using_cuda ? 8 : 7,
+#             ),
+#         ), # ρ(I - P⁻¹ * A) ≈ 0.2249
+#         (
+#             "block Jacobi (diagonal)",
+#             MatrixFields.StationaryIterativeSolve(;
+#                 P_alg = MatrixFields.BlockDiagonalPreconditioner(),
+#                 n_iters = 7,
+#             ),
+#         ), # ρ(I - P⁻¹ * A) ≈ 0.1450
+#         (
+#             "damped block Jacobi (diagonal)",
+#             MatrixFields.StationaryIterativeSolve(;
+#                 P_alg = MatrixFields.WeightedPreconditioner(
+#                     scaled_identity_matrix(FT(1.002)),
+#                     MatrixFields.BlockDiagonalPreconditioner(),
+#                 ),
+#                 n_iters = 7,
+#             ),
+#         ), # ρ(I - P⁻¹ * A) ≈ 0.1427
+#         (
+#             "block arrowhead",
+#             MatrixFields.StationaryIterativeSolve(;
+#                 P_alg = MatrixFields.BlockArrowheadPreconditioner(@name(c)),
+#                 n_iters = 6,
+#             ),
+#         ), # ρ(I - P⁻¹ * A) ≈ 0.1356
+#         (
+#             "damped block arrowhead",
+#             MatrixFields.StationaryIterativeSolve(;
+#                 P_alg = MatrixFields.BlockArrowheadPreconditioner(
+#                     @name(c);
+#                     P_alg₁ = MatrixFields.WeightedPreconditioner(
+#                         scaled_identity_matrix(FT(1.0001)),
+#                         MatrixFields.MainDiagonalPreconditioner(),
+#                     ),
+#                 ),
+#                 n_iters = 6,
+#             ),
+#         ), # ρ(I - P⁻¹ * A) ≈ 0.1355
+#         (
+#             "block arrowhead Schur complement",
+#             MatrixFields.ApproximateBlockArrowheadIterativeSolve(
+#                 @name(c);
+#                 n_iters = 3,
+#             ),
+#         ), # ρ(I - P⁻¹ * A) ≈ 0.0009
+#         (
+#             "damped block arrowhead Schur complement",
+#             MatrixFields.ApproximateBlockArrowheadIterativeSolve(
+#                 @name(c);
+#                 P_alg₁ = MatrixFields.WeightedPreconditioner(
+#                     scaled_identity_matrix(FT(1.09)),
+#                     MatrixFields.MainDiagonalPreconditioner(),
+#                 ),
+#                 n_iters = 2,
+#             ),
+#         ), # ρ(I - P⁻¹ * A) ≈ 0.000006
+#     )
+#         test_field_matrix_solver(;
+#             test_name = "approximate iterative solve with $P_name \
+#                          preconditioning for a block matrix with tri-diagonal, \
+#                          quad-diagonal, bi-diagonal, and penta-diagonal blocks",
+#             alg,
+#             A = MatrixFields.FieldMatrix(
+#                 (@name(c), @name(c)) => ᶜᶜmat3,
+#                 (@name(c), @name(f)) => ᶜᶠmat4,
+#                 (@name(f), @name(c)) => ᶠᶜmat2,
+#                 (@name(f), @name(f)) => ᶠᶠmat5,
+#             ),
+#             b = Fields.FieldVector(; c = ᶜvec, f = ᶠvec),
+#             use_rel_error = true,
+#         )
+#     end
+
+#     @testset "approximate iterative solve with debugging" begin
+#         Logging.with_logger(Logging.SimpleLogger(stderr, Logging.Debug)) do
+#             # Recreate the setup from the previous unit test.
+#             alg = MatrixFields.ApproximateBlockArrowheadIterativeSolve(
+#                 @name(c);
+#                 P_alg₁ = MatrixFields.WeightedPreconditioner(
+#                     scaled_identity_matrix(FT(1.09)),
+#                     MatrixFields.MainDiagonalPreconditioner(),
+#                 ),
+#                 n_iters = 2,
+#             )
+#             A = MatrixFields.FieldMatrix(
+#                 (@name(c), @name(c)) => ᶜᶜmat3,
+#                 (@name(c), @name(f)) => ᶜᶠmat4,
+#                 (@name(f), @name(c)) => ᶠᶜmat2,
+#                 (@name(f), @name(f)) => ᶠᶠmat5,
+#             )
+#             b = Fields.FieldVector(; c = ᶜvec, f = ᶠvec)
+
+#             x = similar(b)
+#             solver = FieldMatrixSolver(alg, A, b)
+#             args = (solver, x, A, b)
+
+#             # Compare the debugging logs to RegEx strings. Note that debugging the
+#             # spectral radius is currently not possible on GPUs.
+#             spectral_radius_logs =
+#                 using_cuda ? () : ((:debug, r"ρ\(I \- inv\(P\) \* A\) ≈"),)
+#             error_norm_logs = (
+#                 (:debug, r"||x[0] - x'||₂ ≈"),
+#                 (:debug, r"||x[1] - x'||₂ ≈"),
+#                 (:debug, r"||x[2] - x'||₂ ≈"),
+#             )
+#             logs = (spectral_radius_logs..., error_norm_logs...)
+#             @test_logs logs... min_level = Logging.Debug field_matrix_solve!(
+#                 args...,
+#             )
+#         end
+#     end
+# end
+
+# @testset "FieldMatrixSolver ClimaAtmos-Based Tests" begin
+#     FT = Float64
+#     center_space, face_space = test_spaces(FT)
+#     surface_space = Spaces.level(face_space, half)
+
+#     seed!(1) # ensures reproducibility
+
+#     ᶜvec = random_field(FT, center_space)
+#     ᶠvec = random_field(FT, face_space)
+#     sfc_vec = random_field(FT, surface_space)
+
+#     # Make each random square matrix diagonally dominant in order to avoid large
+#     # large roundoff errors when computing its inverse. Scale the non-square
+#     # matrices by the same amount as the square matrices.
+#     λ = 10 # scale factor
+#     ᶜᶜmat1 = random_field(DiagonalMatrixRow{FT}, center_space) ./ λ .+ (I,)
+#     ᶜᶠmat2 = random_field(BidiagonalMatrixRow{FT}, center_space) ./ λ
+#     ᶠᶜmat2 = random_field(BidiagonalMatrixRow{FT}, face_space) ./ λ
+#     ᶜᶜmat3 = random_field(TridiagonalMatrixRow{FT}, center_space) ./ λ .+ (I,)
+#     ᶠᶠmat3 = random_field(TridiagonalMatrixRow{FT}, face_space) ./ λ .+ (I,)
+
+#     e¹² = Geometry.Covariant12Vector(1, 1)
+#     e³ = Geometry.Covariant3Vector(1)
+#     e₃ = Geometry.Contravariant3Vector(1)
+
+#     ρχ_unit = (; ρq_tot = 1, ρq_liq = 1, ρq_ice = 1, ρq_rai = 1, ρq_sno = 1)
+#     ρaχ_unit =
+#         (; ρaq_tot = 1, ρaq_liq = 1, ρaq_ice = 1, ρaq_rai = 1, ρaq_sno = 1)
+
+#     dry_center_gs_unit = (; ρ = 1, ρe_tot = 1, uₕ = e¹²)
+#     center_gs_unit = (; dry_center_gs_unit..., ρatke = 1, ρχ = ρχ_unit)
+#     center_sgsʲ_unit = (; ρa = 1, ρae_tot = 1, ρaχ = ρaχ_unit)
+
+#     ᶜᶜmat3_uₕ_scalar = ᶜᶜmat3 .* (e¹²,)
+#     ᶠᶜmat2_u₃_scalar = ᶠᶜmat2 .* (e³,)
+#     ᶜᶠmat2_scalar_u₃ = ᶜᶠmat2 .* (e₃',)
+#     ᶜᶠmat2_uₕ_u₃ = ᶜᶠmat2 .* (e¹² * e₃',)
+#     ᶠᶠmat3_u₃_u₃ = ᶠᶠmat3 .* (e³ * e₃',)
+#     ᶜᶜmat3_ρχ_scalar = map(Base.Fix1(map, Base.Fix2(⊠, ρχ_unit)), ᶜᶜmat3)
+#     ᶜᶜmat3_ρaχ_scalar = map(Base.Fix1(map, Base.Fix2(⊠, ρaχ_unit)), ᶜᶜmat3)
+#     ᶜᶠmat2_ρχ_u₃ = map(Base.Fix1(map, Base.Fix2(⊠, ρχ_unit ⊠ e₃')), ᶜᶠmat2)
+#     ᶜᶠmat2_ρaχ_u₃ = map(Base.Fix1(map, Base.Fix2(⊠, ρaχ_unit ⊠ e₃')), ᶜᶠmat2)
+#     # We need to use Fix1 and Fix2 instead of defining anonymous functions in
+#     # order for the result of map to be inferrable.
+
+#     b_dry_dycore = Fields.FieldVector(;
+#         c = ᶜvec .* (dry_center_gs_unit,),
+#         f = ᶠvec .* ((; u₃ = e³),),
+#     )
+
+#     b_moist_dycore_diagnostic_edmf = Fields.FieldVector(;
+#         c = ᶜvec .* (center_gs_unit,),
+#         f = ᶠvec .* ((; u₃ = e³),),
+#     )
+
+#     b_moist_dycore_prognostic_edmf_prognostic_surface = Fields.FieldVector(;
+#         sfc = sfc_vec .* ((; T = 1),),
+#         c = ᶜvec .* ((; center_gs_unit..., sgsʲs = (center_sgsʲ_unit,)),),
+#         f = ᶠvec .* ((; u₃ = e³, sgsʲs = ((; u₃ = e³),)),),
+#     )
+
+#     test_field_matrix_solver(;
+#         test_name = "similar solve to ClimaAtmos's dry dycore with implicit \
+#                      acoustic waves",
+#         alg = MatrixFields.BlockArrowheadSolve(@name(c)),
+#         A = MatrixFields.FieldMatrix(
+#             (@name(c.ρ), @name(c.ρ)) => I,
+#             (@name(c.ρe_tot), @name(c.ρe_tot)) => I,
+#             (@name(c.uₕ), @name(c.uₕ)) => I,
+#             (@name(c.ρ), @name(f.u₃)) => ᶜᶠmat2_scalar_u₃,
+#             (@name(c.ρe_tot), @name(f.u₃)) => ᶜᶠmat2_scalar_u₃,
+#             (@name(f.u₃), @name(c.ρ)) => ᶠᶜmat2_u₃_scalar,
+#             (@name(f.u₃), @name(c.ρe_tot)) => ᶠᶜmat2_u₃_scalar,
+#             (@name(f.u₃), @name(f.u₃)) => ᶠᶠmat3_u₃_u₃,
+#         ),
+#         b = b_dry_dycore,
+#     )
+
+#     test_field_matrix_solver(;
+#         test_name = "similar solve to ClimaAtmos's dry dycore with implicit \
+#                      acoustic waves and diffusion",
+#         alg = MatrixFields.ApproximateBlockArrowheadIterativeSolve(
+#             @name(c);
+#             n_iters = 6,
+#         ),
+#         A = MatrixFields.FieldMatrix(
+#             (@name(c.ρ), @name(c.ρ)) => I,
+#             (@name(c.ρe_tot), @name(c.ρe_tot)) => ᶜᶜmat3,
+#             (@name(c.uₕ), @name(c.uₕ)) => ᶜᶜmat3,
+#             (@name(c.ρ), @name(f.u₃)) => ᶜᶠmat2_scalar_u₃,
+#             (@name(c.ρe_tot), @name(f.u₃)) => ᶜᶠmat2_scalar_u₃,
+#             (@name(f.u₃), @name(c.ρ)) => ᶠᶜmat2_u₃_scalar,
+#             (@name(f.u₃), @name(c.ρe_tot)) => ᶠᶜmat2_u₃_scalar,
+#             (@name(f.u₃), @name(f.u₃)) => ᶠᶠmat3_u₃_u₃,
+#         ),
+#         b = b_dry_dycore,
+#     )
+
+#     test_field_matrix_solver(;
+#         test_name = "similar solve to ClimaAtmos's moist dycore + diagnostic \
+#                      EDMF with implicit acoustic waves and SGS fluxes",
+#         alg = MatrixFields.ApproximateBlockArrowheadIterativeSolve(
+#             @name(c);
+#             n_iters = 6,
+#         ),
+#         A = MatrixFields.FieldMatrix(
+#             (@name(c.ρ), @name(c.ρ)) => I,
+#             (@name(c.ρe_tot), @name(c.ρe_tot)) => ᶜᶜmat3,
+#             (@name(c.ρatke), @name(c.ρatke)) => ᶜᶜmat3,
+#             (@name(c.ρχ), @name(c.ρχ)) => ᶜᶜmat3,
+#             (@name(c.uₕ), @name(c.uₕ)) => ᶜᶜmat3,
+#             (@name(c.ρ), @name(f.u₃)) => ᶜᶠmat2_scalar_u₃,
+#             (@name(c.ρe_tot), @name(f.u₃)) => ᶜᶠmat2_scalar_u₃,
+#             (@name(c.ρatke), @name(f.u₃)) => ᶜᶠmat2_scalar_u₃,
+#             (@name(c.ρχ), @name(f.u₃)) => ᶜᶠmat2_ρχ_u₃,
+#             (@name(f.u₃), @name(c.ρ)) => ᶠᶜmat2_u₃_scalar,
+#             (@name(f.u₃), @name(c.ρe_tot)) => ᶠᶜmat2_u₃_scalar,
+#             (@name(f.u₃), @name(f.u₃)) => ᶠᶠmat3_u₃_u₃,
+#         ),
+#         b = b_moist_dycore_diagnostic_edmf,
+#     )
+
+#     test_field_matrix_solver(;
+#         test_name = "similar solve to ClimaAtmos's moist dycore + prognostic \
+#                      EDMF + prognostic surface temperature with implicit \
+#                      acoustic waves and SGS fluxes",
+#         alg = MatrixFields.BlockLowerTriangularSolve(
+#             @name(c.sgsʲs),
+#             @name(f.sgsʲs);
+#             alg₁ = MatrixFields.BlockArrowheadSolve(@name(c)),
+#             alg₂ = MatrixFields.ApproximateBlockArrowheadIterativeSolve(
+#                 @name(c);
+#                 n_iters = 6,
+#             ),
+#         ),
+#         A = MatrixFields.FieldMatrix(
+#             # GS-GS blocks:
+#             (@name(sfc), @name(sfc)) => I,
+#             (@name(c.ρ), @name(c.ρ)) => I,
+#             (@name(c.ρe_tot), @name(c.ρe_tot)) => ᶜᶜmat3,
+#             (@name(c.ρatke), @name(c.ρatke)) => ᶜᶜmat3,
+#             (@name(c.ρχ), @name(c.ρχ)) => ᶜᶜmat3,
+#             (@name(c.uₕ), @name(c.uₕ)) => ᶜᶜmat3,
+#             (@name(c.ρ), @name(f.u₃)) => ᶜᶠmat2_scalar_u₃,
+#             (@name(c.ρe_tot), @name(f.u₃)) => ᶜᶠmat2_scalar_u₃,
+#             (@name(c.ρatke), @name(f.u₃)) => ᶜᶠmat2_scalar_u₃,
+#             (@name(c.ρχ), @name(f.u₃)) => ᶜᶠmat2_ρχ_u₃,
+#             (@name(f.u₃), @name(c.ρ)) => ᶠᶜmat2_u₃_scalar,
+#             (@name(f.u₃), @name(c.ρe_tot)) => ᶠᶜmat2_u₃_scalar,
+#             (@name(f.u₃), @name(f.u₃)) => ᶠᶠmat3_u₃_u₃,
+#             # GS-SGS blocks:
+#             (@name(c.ρe_tot), @name(c.sgsʲs.:(1).ρae_tot)) => ᶜᶜmat3,
+#             (@name(c.ρχ.ρq_tot), @name(c.sgsʲs.:(1).ρaχ.ρaq_tot)) => ᶜᶜmat3,
+#             (@name(c.ρχ.ρq_liq), @name(c.sgsʲs.:(1).ρaχ.ρaq_liq)) => ᶜᶜmat3,
+#             (@name(c.ρχ.ρq_ice), @name(c.sgsʲs.:(1).ρaχ.ρaq_ice)) => ᶜᶜmat3,
+#             (@name(c.ρχ.ρq_rai), @name(c.sgsʲs.:(1).ρaχ.ρaq_rai)) => ᶜᶜmat3,
+#             (@name(c.ρχ.ρq_sno), @name(c.sgsʲs.:(1).ρaχ.ρaq_sno)) => ᶜᶜmat3,
+#             (@name(c.ρe_tot), @name(c.sgsʲs.:(1).ρa)) => ᶜᶜmat3,
+#             (@name(c.ρatke), @name(c.sgsʲs.:(1).ρa)) => ᶜᶜmat3,
+#             (@name(c.ρχ), @name(c.sgsʲs.:(1).ρa)) => ᶜᶜmat3_ρχ_scalar,
+#             (@name(c.uₕ), @name(c.sgsʲs.:(1).ρa)) => ᶜᶜmat3_uₕ_scalar,
+#             (@name(c.ρe_tot), @name(f.sgsʲs.:(1).u₃)) => ᶜᶠmat2_scalar_u₃,
+#             (@name(c.ρatke), @name(f.sgsʲs.:(1).u₃)) => ᶜᶠmat2_scalar_u₃,
+#             (@name(c.ρχ), @name(f.sgsʲs.:(1).u₃)) => ᶜᶠmat2_ρχ_u₃,
+#             (@name(c.uₕ), @name(f.sgsʲs.:(1).u₃)) => ᶜᶠmat2_uₕ_u₃,
+#             (@name(f.u₃), @name(c.sgsʲs.:(1).ρa)) => ᶠᶜmat2_u₃_scalar,
+#             (@name(f.u₃), @name(f.sgsʲs.:(1).u₃)) => ᶠᶠmat3_u₃_u₃,
+#             # SGS-SGS blocks:
+#             (@name(c.sgsʲs.:(1).ρa), @name(c.sgsʲs.:(1).ρa)) => I,
+#             (@name(c.sgsʲs.:(1).ρae_tot), @name(c.sgsʲs.:(1).ρae_tot)) => I,
+#             (@name(c.sgsʲs.:(1).ρaχ), @name(c.sgsʲs.:(1).ρaχ)) => I,
+#             (@name(c.sgsʲs.:(1).ρa), @name(f.sgsʲs.:(1).u₃)) =>
+#                 ᶜᶠmat2_scalar_u₃,
+#             (@name(c.sgsʲs.:(1).ρae_tot), @name(f.sgsʲs.:(1).u₃)) =>
+#                 ᶜᶠmat2_scalar_u₃,
+#             (@name(c.sgsʲs.:(1).ρaχ), @name(f.sgsʲs.:(1).u₃)) => ᶜᶠmat2_ρaχ_u₃,
+#             (@name(f.sgsʲs.:(1).u₃), @name(c.sgsʲs.:(1).ρa)) =>
+#                 ᶠᶜmat2_u₃_scalar,
+#             (@name(f.sgsʲs.:(1).u₃), @name(c.sgsʲs.:(1).ρae_tot)) =>
+#                 ᶠᶜmat2_u₃_scalar,
+#             (@name(f.sgsʲs.:(1).u₃), @name(f.sgsʲs.:(1).u₃)) => ᶠᶠmat3_u₃_u₃,
+#         ),
+#         b = b_moist_dycore_prognostic_edmf_prognostic_surface,
+#     )
+# end
+
+@testset "FieldMatrixSolver with CenterFiniteDifferenceSpace" begin
+    # Set up FiniteDifferenceSpace
+    FT = Float32
+    zmax = FT(0)
+    zmin = FT(-0.35)
+    nelems = 5
+
+    context = ClimaComms.context()
+    z_domain = Domains.IntervalDomain(
+        Geometry.ZPoint(zmin),
+        Geometry.ZPoint(zmax);
+        boundary_names = (:bottom, :top),
     )
+    z_mesh = Meshes.IntervalMesh(z_domain, nelems = nelems)
+    z_topology = Topologies.IntervalTopology(context, z_mesh)
+    space = Spaces.CenterFiniteDifferenceSpace(z_topology)
+
+    # Create a field containing a `TridiagonalMatrixRow` at each point
+    tridiag_type = MatrixFields.TridiagonalMatrixRow{FT}
+    tridiag_field = Fields.Field(tridiag_type, space)
+
+    # Set up objects for matrix solve
+    A = MatrixFields.FieldMatrix((@name(_), @name(_)) => tridiag_field)
+    field = Fields.ones(space)
+    b = Fields.FieldVector(; _ = field)
+    x = similar(b)
+    solver =
+        MatrixFields.FieldMatrixSolver(MatrixFields.BlockDiagonalSolve(), A, b)
+
+    # Run matrix solve
+    MatrixFields.field_matrix_solve!(solver, x, A, b)
 end