Factor timestep out of solver operators

[Thomas Bendall (tommbendall)]

PR Summary

Sci/Tech Reviewer: thomasmelvin
Code Reviewer: Lottie Turner (@mo-lottieturner)

This PR performs some refactoring of the si_operators used by Gungho's linear solver. Many of these operators currently include a scaling of the form tau*dt, and this PR pulls those scalings out so that they are applied in-place rather than included in the operators.

Motivation: the new TR-BDF2 timestepper being developed in #379 calls the solver multiple times with different values of tau*dt. Factoring tau*dt out of these operators means that the operators do not need duplicating multiple times for the TR-BDF2 timestepper.

However, some of the si_operators can't be straightforwardly refactored and will need duplicating. This PR anticipates this by storing these operators in an appropriate inventory.

Linked-To:

• LFRic-Core #435
• LFRic-JEDI 1323

Blocks: #379

Overview of Changes

Primary operators — tau*dt dependence

Operator	Factor of tau*dt
m3_rho_star	none
m3_exner_star	none
p3theta	none
rho_at_u	none
div_star	tau_u * dt * c_p
p2theta	tau_u * dt * c_p
ptheta2	tau_t * dt
ptheta2v	tau_t * dt
compound_div	tau_r * dt
eliminated_q22	tau_t * tau_u * (dt)^2 * c_p
eliminated_q32	tau_r * dt (from compound_div) + tau_t * dt
eliminated_q2t	tau_u * dt * c_p

Derived fields — dt dependence

Field	dt dependence
Hb_lumped_inv	Standard: inverse of $N_u(M_u + P_{2\theta}M_t^{-1}P_{\theta2}^v)\mathbf{1}$ — contains tau_u*tau_t*(dt)^2*c_p from the buoyancy term; M_u contributes a dt-independent offset, so the dependence is nonlinear. Analytic: inverse of $(N_u M_u - Q_{22})\mathbf{1}$ — same structure with tau_t*tau_u*(dt)^2*c_p from Q_22.
Helmholtz_operator	Assembled from Hb_lumped_inv (nonlinear in dt), div_star (tau_u*dt*c_p), compound_div or eliminated_q32 (tau_r*dt, tau_t*dt), and m3_exner_star (none). Dominant off-diagonal structure $\mathcal{A} = -N_u H_B^{-1} D^*$ is already $O(dt)$. Overall dependence nonlinear through $H_B^{-1}$.
Helm_diag	Pointwise inverse of the central stencil component of Helmholtz_operator; inherits the same nonlinear dt dependence. When normalisation is off, set to 1 (independent of dt).
tri_precon	Linear combination of Helmholtz_operator stencil components, optionally scaled by Helm_diag; inherits the same nonlinear dt dependence as Helmholtz_operator.

PR review walkthrough

Core refactor: SI operators are assembled unscaled, tau*dt applied later

• si_operators_alg_mod.x90: The underlying (div_star, p2theta, ptheta2, ptheta2v, compound_div, eliminated_q22, , eliminated_q2t) are now assembled without the tau*dt factor baked in. The scaling is applied later, at the point each operator is applied (matrix-vector kernels, mixed_operator/apply_mixed_operator, apply_elim_mixed_lp, opt_apply_variable_hx, etc.), via new scalar arguments (const_u, const_t, const_r, const_div, const_theta).
• eliminated_q32 has been split into eliminated_q32_rho and eliminated_q32_theta
• Derived fields with a non-linear dependence on tau*dt factors are stored in an inventory, so that there can
  be multiple copies of these.
• New getters are added (get_tau_u_dt_cp, get_tau_t_dt, get_tau_r_dt, get_tau_tu_dt2_cp). These take a stepper argument, since the constants depend on which timestepper scheme is active. Operator-assembly routines (get_tri_precon_*, get_helmholtz_operator_*, get_helm_diag_*, get_hb_lumped_inv_*) and create_si_operators were similarly threaded with stepper so the correct constants can be looked up from the inventory.
• dycore_constants_mod.x90: added the stepper_siqn enumerator, used to select which SI-scheme constants apply.

Gungho kernels/algorithms

• apply_mixed_operator_kernel_mod.F90: kernel metadata and argument list extended to take const_u, const_t, const_r; the scaling multiplications that used to be implicit in the operator are now explicit in the kernel body.
• mixed_operator_alg_mod.x90, mixed_schur_preconditioner_alg_mod.x90, pressure_operator_alg_mod.x90, pressure_precon_alg_mod.x90, semi_implicit_solver_alg_mod.x90: updated to source and pass through the new scaling constants and stepper argument at the call sites that build/apply these operators.
• apply_elim_mixed_lp_operator_kernel_mod and opt_apply_variable_hx_kernel_mod: same pattern — added scaling scalar arguments (const_r; const_div/const_theta) and renamed q32 → q32_rho to reflect the analytic-elimination naming.
• sample_eliminated_theta_q32_kernel_mod.F90 / project_eliminated_theta_q32_kernel_mod.F90: minor follow-on tweaks for the renamed/rescoped operator.

Results

By refactoring code, this PR changes the KGOs of all tests that use the linear solver. However these changes are only at the bit-level and do not impact science.

These PDFs contain the plots from the test-suite, compared with those from vn3.2, demonstrating that this change is indeed science-neutral:

• gungho:
  results_si_operator_dt_gungho.pdf
• lfric_atm:
  results_si_operator_dt_lfric_atm.pdf

Performance

By moving the tau*dt factors out of the operators, these need applying during the application of the solver. The following plots show timings for lfric_atm and the solver, from jobs in the lfric_atm weekly suite. Each was run three times to get a sense for variability. Any increase in solver cost seems to be very small, and within the variability of the total runtime.

Code Quality Checklist

• I have performed a self-review of my own code
• My code follows the project's style guidelines
• Comments have been included that aid understanding and enhance the readability of the code
• My changes generate no new warnings
• All automated checks in the CI pipeline have completed successfully

Testing

• I have tested this change locally, using the LFRic Apps rose-stem suite
• If any tests fail (rose-stem or CI) the reason is understood and acceptable (e.g. kgo changes)
• I have added tests to cover new functionality as appropriate (e.g. system tests, unit tests, etc.)
• Any new tests have been assigned an appropriate amount of compute resource and have been allocated to an appropriate testing group (i.e. the developer tests are for jobs which use a small amount of compute resource and complete in a matter of minutes)

trac.log

Test Suite Results - lfric_apps - si_operator_dt/run1

Suite Information

Item	Value
Suite Name	si_operator_dt/run1
Suite User	thomas.bendall
Workflow Start	2026-08-10T09:05:05
Groups Run	all

Dependency	Reference	Main Like
casim	MetOffice/casim@2026.07.1	True
jules	MetOffice/jules@2026.07.1	True
lfric_apps	tommbendall/lfric_apps@si_operator_dt	False
lfric_core	tommbendall/lfric_core@932de27	True
moci	MetOffice/moci@2026.07.1	True
SimSys_Scripts	MetOffice/SimSys_Scripts@2026.07.1	True
socrates	MetOffice/socrates@2026.07.1	True
socrates-spectral	MetOffice/socrates-spectral@2026.07.1	True
ukca	MetOffice/ukca@2026.07.1	True

Task Information

✅ succeeded tasks - 1595

Security Considerations

• I have reviewed my changes for potential security issues
• Sensitive data is properly handled (if applicable)
• Authentication and authorisation are properly implemented (if applicable)

Performance Impact

• Performance of the code has been considered and, if applicable, suitable performance measurements have been conducted

AI Assistance and Attribution

• Some of the content of this change has been produced with the assistance of Generative AI tool name (e.g., Met Office Github Copilot Enterprise, Github Copilot Personal, ChatGPT GPT-4, etc) and I have followed the Simulation Systems AI policy (including attribution labels)

Copilot was used under human design and governance to perform this refactor and to fix the adjoint patches.

Documentation

• Where appropriate I have updated documentation related to this change and confirmed that it builds correctly

PSyclone Approval

• If you have edited any PSyclone-related code (e.g. PSyKAl-lite, Kernel interface, optimisation scripts, LFRic data structure code) then please contact the TCD Team

Sci/Tech Review

• I understand this area of code and the changes being added
• The proposed changes correspond to the pull request description
• Documentation is sufficient (do documentation papers need updating)
• Sufficient testing has been completed

(Please alert the code reviewer via a tag when you have approved the SR)

Code Review

• All dependencies have been resolved
• Related Issues have been properly linked and addressed
• CLA compliance has been confirmed
• Code quality standards have been met
• Tests are adequate and have passed
• Documentation is complete and accurate
• Security considerations have been addressed
• Performance impact is acceptable

[cjohnson-pi]

Question: should/does this have an impact on gungho/source/algorithm/solver/lam_rhs_alg_mod.x90 ? Thanks.

[Thomas Bendall (tommbendall)]

Question: should/does this have an impact on gungho/source/algorithm/solver/lam_rhs_alg_mod.x90 ? Thanks.

No it doesn't impact lam_rhs_alg_mod, as I have not factored the timestep out of the various science-adjusted W2 mass matrices

[DanStoneMO]

Looks like this will need a PR in JEDI, I will link it here once ready

[DanStoneMO]

Linked JEDI PR is now live at: https://github.com/MetOffice/lfric-jedi/pull/1323

[Steven Sandbach (ss421)]

I see that you have referenced the change of results:

"By refactoring code, this PR changes the KGOs of all tests that use the linear solver. However these changes are only at the bit-level and do not impact science."

I assume that the change to results is due to an order of operations change? Thanks.

[Thomas Bendall (tommbendall)]

I see that you have referenced the change of results:

"By refactoring code, this PR changes the KGOs of all tests that use the linear solver. However these changes are only at the bit-level and do not impact science."

I assume that the change to results is due to an order of operations change? Thanks.

Yes that's right, the point at which fields are multiplied by tau*dt changes

[thomasmelvin]

                 inc_X_times_Y(rhs_tmp, m2_diag),                                &
                 inc_a_times_X(const_u, rhs_tmp),                                  &

could become
inc_aX_times_Y(rhs_rtmp, const_u, m2_diag), &

[thomasmelvin]

                 inc_a_times_X( const_u, rhs_tmp ),                                  &
                 inc_X_plus_Y( rhs_rsol(igh_u), rhs_tmp ),                           &

could become
inc_X_plus_bY(rhs_rsol(igh_u), const_u, rhs_tmp), &

[thomasmelvin]

I think the multiplication by const_u could be taken inside the scaled_matrix_vector kernel to avoid the extra built in

[thomasmelvin]

Generally looks good, I've made a few suggestions where I think the code could be optimised a little more (generally removing a built-in of extra temporary field)

[thomasmelvin]

I think the multiplication by const_u could be taken inside the scaled_matrix_vector kernel to avoid the extra built in

[thomasmelvin]

I think the multiplication by const_u could be taken inside the scaled_matrix_vector kernel to avoid the extra built in

[thomasmelvin]

               X_times_Y(uvw_norm, u_normalisation, Hb_lumped_inv),       &
               inc_a_times_X(const_u, uvw_norm),  

could become
aX_times_Y(uvw_norm, const_u, u_normalisation, Hb_lumped_inv), &

[Thomas Bendall (tommbendall)]

So there actually isn't an aX_times_Y built-in! I think this is still the simplest combination. Hopefully a clever compiler will fuse the loops and there will be no impact...

[thomasmelvin]

I'm a bit concerned by the replacement of a single kernel with 5 kernels in terms of performance. I think it would be more efficient to create a new kernel to compute r_p = const_rq32_rho_opself%rhs_u + const_tq32_theta_opself%rhs_u
or more generally to compute
y = aMx + bNx for scalars a, b, matrices M, N and vectors x, y
this would also remove the need to create the r_p_theta vector

[thomasmelvin]

Do all these copy_field_array mean we are creating multiple copies of the fields? i.e. one in the inventory and a local copy? if so it feels it would be better to expose the field from the inventory to the algorithm, something like tri_precon => tri_precon_inventory%get_field_array(stepper)

[Thomas Bendall (tommbendall)]

I've now done this. One complication is that the elements of the field array are on different meshes. I couldn't think of a nice interface in the inventory to handle that! So instead I'm reinitialising the fields once they've been created.

I think this is not too bad since this routine is only done once at model start-up. I could make this cleaner if I stored the fields in the inventory by mesh and by stepper, but this would require changed interfaces everywhere these fields are used... I think that would be neater but is probably best left to a future PR

[thomasmelvin]

You could replace the used of scaled ones by
matrix_vector_kernel_type(m2_u, ones, eliminated_q22(level)), &
inc_X_minus_bY(Hb_lumped_inc(level), tau_tu_dt2_cp, m2_u), &

[thomasmelvin]

Again, you could replace the used of scaled ones by
dg_matrix_vector_kernel_type(theta_adv_term, ones, ptheta2v(level)), &
inc_X_times_bY(theta_adv_term, tau_t_dt, t_normalisation), &
and remove the extra invoke

[Thomas Bendall (tommbendall)]

Generally looks good, I've made a few suggestions where I think the code could be optimised a little more (generally removing a built-in of extra temporary field)

Thanks very much for the review and the suggested improvements. I have:

• added an extra argument to scaled_matrix_vector_kernel_mod.F90 to avoid needing a following
• added a new kernel
• made your other suggested changes of combining built-ins

I also realised I had forgotten to do one important thing: split the compute_si_operators into a generic part which computes those operators shared between steppers, and a part that computers the fields that different between steppers. I have now done that too.

The changes were all made in this commit: b807692

I have also given specific replies about two changes:

• there actually isn't an aX_times_Y built-in so I have left the suggestion to do that
• I have changed the copy_field_array calls for the inventory to add_field_array, but flagged that this means that all fields in the array are initialised on the fine mesh. However different fields in the array are needed on different meshes, so I reinitialise them.

[DrTVockerodtMO]

Thanks Tom, adjoint code changes look good. I've suggested only one minor change but in general I approve.

[github-actions]

⚠️ Hello Thomas Bendall (@tommbendall)!

Your CLA signature was found on the base branch, but you appear to have modified the CONTRIBUTORS.md file in this PR.

Please do not edit the CONTRIBUTORS.md file. If you have already signed the CLA, revert changes to the file and your signature will be picked up.