FSILS RCS preconditioner alters small matrix diagonals and degrades convergence

Description

The FSILS row and column scaling (RCS) preconditioner changes small, unconstrained matrix diagonal entries while applying Dirichlet boundary conditions. This introduces rounding error into the operator before row and column scaling begins.

This is a shared linear-solver issue. The same failure has been reproduced with the existing heatS equation and the Darcy equation reviewed in PR #620. For an exactly representable linear solution on a small triangular mesh, both RCS runs have a maximum nodal error of approximately 0.383 after five outer iterations. Changing only the preconditioner to fsils gives zero error at the recorded output precision.

The problematic expression appears in every scalar and block-size branch of precond_rcs. Other physics using this FSILS preconditioner may therefore be affected when their assembled matrix diagonals are small.

Reproduction

Solver-level reproduction

Use a 2D unit square with a structured 5-by-3 node grid and 16 linear triangles:

• Prescribe u = 1 on the left edge and u = 0 on the right edge.
• Leave the top and bottom edges at their natural zero-flux boundary condition.
• Use a zero initial field and no source. The exact steady solution is u(x, y) = 1 - x, which this mesh can represent exactly.
• For heatS, set Density = 0 and Conductivity = 5e-16. The zero density deliberately removes the transient term for this steady diffusion test.
• For the equivalent Darcy check, set fluid density to 0.5, permeability to 1e-15, viscosity to 1, and compressibility to 0.
• Use one time step of size 0.1, spectral radius 0, and at most five outer iterations with tolerance 1e-10.
• Select the FSILS linear algebra backend with CG, relative tolerance 1e-12, absolute tolerance 1e-30, and at most 30 linear iterations. Compare rcs and fsils preconditioners.

Measured maximum nodal error, max_i |u_i - (1 - x_i)|:

Physics	Preconditioner	Maximum nodal error
heatS	rcs	3.8263694109545121e-01
heatS	fsils	0.0000000000000000e+00
darcy	rcs	3.8263694109545121e-01
darcy	fsils	0.0000000000000000e+00

The RCS runs reach the five-iteration outer limit without meeting the outer convergence tolerance, even though the inner linear solves report very small residual ratios. The resulting corrections produce slow outer convergence and inaccurate output at the iteration limit.

Tightening the linear absolute tolerance to 1e-30 does not eliminate the problem, as the table shows. In a separate Darcy run using 30 outer iterations and the default linear absolute tolerance, the maximum nodal error was still 3.9171435852137615e-02.

Expected behavior

Dirichlet handling should leave unconstrained matrix diagonals unchanged and set constrained diagonals to one before applying the intended row and column scaling. Preconditioning should preserve the equivalent linear system.

Uniformly reducing the diffusion coefficient should preserve the exact solution u = 1 - x. The test should converge to the requested accuracy without the large coefficient-dependent error shown above.

Additional context

Shared code path and numerical mechanism

The relevant code is precond_rcs in Code/Source/linear_solver/precond.cpp. After masking constrained rows and columns, it performs the following update for a scalar diagonal:

Val(0,d) = Wr(0,Ac) * (Val(0,d) - 1.0) + 1.0;

The same expression is repeated for block sizes 2, 3, 4, and the general block-size case.

For an unconstrained degree of freedom, Wr = 1, so this evaluates (a - 1.0) + 1.0 for a diagonal value a. Although algebraically equal to a, the subtraction rounds at the scale of one and can lose the small diagonal value. For example, ordinary double-precision evaluation gives:

Original diagonal a	After (a - 1.0) + 1.0
1e-11	approximately 1.000000082740371e-11
1e-16	approximately 1.1102230246251565e-16
1e-17	0

This operation modifies free diagonals as well as constrained ones. Subsequent scaling cannot recover the lost information.

fsils_solve dispatches to this shared routine, independently of the originating physics assembly. Existing CEP and Stokes test configurations also select FSILS RCS.

Version and environment

• Tested commit: 28f48ece9f79c29df7ae02bb1417e5f7ace69f38.
• precond.cpp is identical to the PR's base commit, 00f5c9204bb5a9991e4c8d18eddd3b1492f3b59a.
• Ubuntu 22.04.5 LTS; GCC 13.1.0; Release build (-O3 -DNDEBUG) with array index checking enabled.
• Open MPI 4.1.2; one MPI rank per solver run.
• Solver linked against VTK 9.1; reproducer used Python VTK 9.3.1 and ASCII mesh input files.

Suggested resolution and regression coverage

Preserve unconstrained diagonals directly and set only constrained diagonals to one, with consistent behavior in every block-size branch.

The standard FSILS diagonal preconditioner (<Preconditioner>fsils</Preconditioner>) is a verified workaround for the two scalar examples in this report.

Code of Conduct

• I agree to follow this project's Code of Conduct and Contributing Guidelines