single_field_solve! / bycolumn bug

[juliasloan25]

Describe the bug

In ClimaLand's RichardsModel implicit solver, I encountered a bug in the MatrixFields solver upon upgrading to ClimaCore v0.14. See this build for a stacktrace.

This was introduced in #1653, perhaps in this line. It is not an issue when using the commit just before this PR (ClimaCore#033e5b0), and is an issue in all subsequent commits (starting with ClimaCore#31e69ef).

PR to add unit test for reproducer: #1750

To Reproduce

import ClimaComms
import ClimaCore: Spaces, MatrixFields, Fields, Domains, Meshes, Topologies, Geometry
import ClimaCore.MatrixFields: @name

FT = Float32
zmax = FT(0)
zmin = FT(-0.35)
nelems = 5

# Make FiniteDifferenceSpace
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)

tridiag_type = MatrixFields.TridiagonalMatrixRow{FT}
# Create a field containing a `TridiagonalMatrixRow` at each point
tridiag_field = Fields.Field(tridiag_type, space)

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)

# this passes with ClimaCore#033e5b0 and fails with ClimaCore#31e69ef
MatrixFields.field_matrix_solve!(solver, x, A, b)

Project

If not using the `examples` project: ``` paste your Project.toml here. ``` ``` paste your Manifest.toml here. ```

System details

Any relevant system information:

• Julia version
• operating system
• modules loaded on cluster (module list)

Related issues / PRs

Please add any relevant links.

[charleskawczynski]

Closed by #1750