Skip to content

TR-BDF2 Timestepper #714

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
Thomas Bendall (tommbendall) wants to merge 40 commits into
base: main
Choose a base branch
from
Open

TR-BDF2 Timestepper #714

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
40 commits
Select commit Hold shift + click to select a range
4a7b1fc
Merge pull request #5 from MetOffice/stable
tommbendall Mar 9, 2026
92ac379
Merge pull request #13 from MetOffice/stable
tommbendall Jul 20, 2026
364168f
first pass at factoring constants out of SI operators
tommbendall Jul 23, 2026
3190af6
all tests passing
tommbendall Jul 24, 2026
306a844
update copyright statement
tommbendall Jul 24, 2026
2b97cf4
update KGOs
tommbendall Jul 24, 2026
db22e60
rose-meta tidy up and upgrade macro
tommbendall Aug 4, 2026
7376d41
fix upgrade macro
tommbendall Aug 4, 2026
c1cc154
fix upgrade macro again
tommbendall Aug 4, 2026
340712d
more rose-meta fixes
tommbendall Aug 4, 2026
5917577
TR-BDF2 compiling. Also add tasks
tommbendall Aug 7, 2026
bd7b8ab
fixes, which get a local example running
tommbendall Aug 7, 2026
a3d4e15
fix typo in tasks
tommbendall Aug 7, 2026
a20692e
fix task names
tommbendall Aug 7, 2026
2f2a264
tau_r_dt fix
tommbendall Aug 7, 2026
f817866
get tau_r_dt from stored value
tommbendall Aug 7, 2026
e151e8f
some fixes, and first attempt at TR-BDF2 diagnostics
tommbendall Aug 8, 2026
b807692
copy over KGOs
tommbendall Aug 8, 2026
189e6b8
Merge branch 'si_operator_dt' into tr_bdf2
tommbendall Aug 8, 2026
830205c
adopt new SI operator split
tommbendall Aug 8, 2026
1c946e7
fix silly error
tommbendall Aug 8, 2026
08396f3
fixes to test cases
tommbendall Aug 10, 2026
0a3cc39
linear and adjoint fixes
tommbendall Aug 11, 2026
35b1081
adjustments for RAL task
tommbendall Aug 11, 2026
9e08fec
revert diagnostic changes, and add documentation
tommbendall Aug 13, 2026
952d14d
new lines at end of file
tommbendall Aug 13, 2026
88292d1
fixes to shallow-water app and gungho unit-tests
tommbendall Aug 13, 2026
89859c6
actually fix shallow water
tommbendall Aug 13, 2026
9d90abf
copilot review suggestions
tommbendall Aug 17, 2026
d3cf600
fix linear integration test, improve config checks, correct upgrade
tommbendall Aug 17, 2026
1a556dc
Merge branch 'tr_bdf2' of github.com:tommbendall/lfric_apps into tr_bdf2
tommbendall Aug 17, 2026
ddaadcd
actually fix upgrade macro
tommbendall Aug 17, 2026
77c1268
urgh silly typo
tommbendall Aug 17, 2026
b08c1e1
ian comments
tommbendall Aug 19, 2026
4b3b8df
typo
tommbendall Aug 19, 2026
de850ab
and update memory
tommbendall Aug 19, 2026
b30b42a
change RAL settings
tommbendall Aug 19, 2026
b2a01f5
documentation improvements, and add solver flux to BDF tracer transport
tommbendall Aug 21, 2026
be6cc39
documentation improvements, and avoid unnecessary flux precomputations
tommbendall Aug 25, 2026
1d12ccf
fix precision issue
tommbendall Aug 25, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
The table of contents is too big for display.
Diff view
Diff view
  •  
  •  
  •  
2 changes: 1 addition & 1 deletion applications/adjoint_tests/example/configuration.nml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -374,7 +374,7 @@ timestep_start='1',
&timestepping
alpha=0.55,
dt=1800,
inner_iterations=1,
inner_iterations_si=1,1
method='semi_implicit',
outer_iterations=2,
runge_kutta_method='forward_euler',
Expand Down
5 changes: 3 additions & 2 deletions applications/adjoint_tests/source/algorithm/core_dynamics/atlt_rhs_alg_mod.x90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,7 @@ module atlt_rhs_alg_mod
use adjoint_test_parameters_mod, only : ls_u_range, ls_theta_range, &
ls_rho_range, ls_exner_range, &
ls_md1_range, ls_md2_range, ls_md3_range
use dycore_constants_mod, only : stepper_siqn

implicit none

Expand Down Expand Up @@ -208,7 +209,7 @@ module atlt_rhs_alg_mod
! Tangent linear
call tl_rhs_alg( rhs, alpha_dt, base_state, state, moist_dyn, &
ls_state, ls_moist_dyn, compute_eos, &
dlayer_rhs, model_clock )
dlayer_rhs, stepper_siqn, model_clock )

! < Mx, Mx >
call invoke( x_innerproduct_x( base_state_inner_prod(igh_u), base_state(igh_u) ), &
Expand Down Expand Up @@ -312,7 +313,7 @@ module atlt_rhs_alg_mod
! Adjoint
call atl_rhs_alg( rhs, alpha_dt, base_state, state, moist_dyn, &
ls_state, ls_moist_dyn, compute_eos, &
dlayer_rhs, model_clock )
dlayer_rhs, stepper_siqn, model_clock )

! < AMx, x >
call invoke( x_innerproduct_y( base_state_base_state_input_inner_prod(igh_u), &
Expand Down
12 changes: 8 additions & 4 deletions ...s/adjoint_tests/source/algorithm/lookup/solver/adjt_apply_helmholtz_op_lookup_alg_mod.x90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,10 @@ module adjt_apply_helmholtz_op_lookup_alg_mod
use integer_field_mod, only: integer_field_type
use adj_lookup_table_mod, only: adj_lookup_table_type
use sci_r_solver_field_vector_mod, only: r_solver_field_vector_type
use si_operators_alg_mod, only: compute_si_operators, get_helmholtz_operator
use si_operators_alg_mod, only: compute_si_operators_common, &
compute_helmholtz_operators, &
get_helmholtz_operator
use dycore_constants_mod, only: stepper_siqn

implicit none

Expand Down Expand Up @@ -119,10 +122,11 @@ module adjt_apply_helmholtz_op_lookup_alg_mod
call invoke(setval_random(init_op_moist_dyn(idx)))
end do

call compute_si_operators( config, init_op_t, init_op_d, init_op_p, &
model_clock, init_op_moist_dyn )
call compute_si_operators_common( config, init_op_t, init_op_d, init_op_p, &
init_op_moist_dyn )
call compute_helmholtz_operators( model_clock, stepper_siqn )

Helmholtz_operator => get_helmholtz_operator(level)
Helmholtz_operator => get_helmholtz_operator(stepper_siqn, level)
lam_mesh = .false.

call vector_x%initialise(vector_space_w3_ptr)
Expand Down
13 changes: 8 additions & 5 deletions applications/adjoint_tests/source/algorithm/solver/adjt_mixed_operator_alg_mod.x90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -32,7 +32,9 @@ module adjt_mixed_operator_alg_mod
use function_space_collection_mod, only: function_space_collection
use moist_dyn_mod, only: num_moist_factors
use fs_continuity_mod, only: W2, W3, Wtheta
use si_operators_alg_mod, only: compute_si_operators
use si_operators_alg_mod, only: compute_si_operators_common, &
compute_helmholtz_operators
use dycore_constants_mod, only: stepper_siqn
use sci_enforce_bc_kernel_mod, only: enforce_bc_kernel_type
use model_clock_mod, only: model_clock_type
use function_space_mod, only: function_space_type
Expand Down Expand Up @@ -115,8 +117,9 @@ contains
call invoke(setval_random(moist_dyn(idx)))
end do

call compute_si_operators( config, rhs(igh_t), rhs(igh_d), rhs(igh_p), &
model_clock, moist_dyn )
call compute_si_operators_common( config, rhs(igh_t), rhs(igh_d), rhs(igh_p), &
moist_dyn )
call compute_helmholtz_operators( model_clock, stepper_siqn )

! Input fields are r_def fields so preliminary work uses field_types

Expand Down Expand Up @@ -144,7 +147,7 @@ contains
call construct_solver_state( vector_mx, rhs_rsol, import_fields=.false. )
call construct_solver_state( vector_amx, rhs_rsol, import_fields=.false. )

mixed_op = mixed_operator_type()
mixed_op = mixed_operator_type(stepper_siqn)

call vector_mx%set_scalar(0.0_r_def)
call vector_amx%set_scalar(0.0_r_def)
Expand Down Expand Up @@ -184,7 +187,7 @@ contains
! Adjoint code
! --------------------------------------------------------------------------

adj_mixed_op = adj_mixed_operator_type()
adj_mixed_op = adj_mixed_operator_type(stepper_siqn)

call adj_mixed_op%apply( vector_mx, vector_amx )

Expand Down
18 changes: 11 additions & 7 deletions ...cations/adjoint_tests/source/algorithm/solver/adjt_mixed_schur_preconditioner_alg_mod.x90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,9 @@ module adjt_mixed_schur_preconditioner_alg_mod
use function_space_collection_mod, only: function_space_collection
use driver_modeldb_mod, only: modeldb_type
use fs_continuity_mod, only: W2, W3, Wtheta
use si_operators_alg_mod, only: compute_si_operators
use si_operators_alg_mod, only: compute_si_operators_common, &
compute_helmholtz_operators
use dycore_constants_mod, only: stepper_siqn
use sci_enforce_bc_kernel_mod, only: enforce_bc_kernel_type
use model_clock_mod, only: model_clock_type
use function_space_mod, only: function_space_type
Expand Down Expand Up @@ -160,17 +162,19 @@ contains
call moisture_fields%get_field( "ls_moist_dyn", ls_moist_dyn_array )
ls_moist_dyn => ls_moist_dyn_array%bundle

call compute_si_operators( modeldb%config, ls_theta, ls_rho, ls_exner, &
model_clock, ls_moist_dyn )
call compute_si_operators_common( modeldb%config, ls_theta, ls_rho, ls_exner, &
ls_moist_dyn )
call compute_helmholtz_operators( model_clock, stepper_siqn )

call create_pressure_preconditioner( rhs, pressure_operator, pressure_preconditioner )
call create_pressure_preconditioner( rhs, pressure_operator, pressure_preconditioner, stepper_siqn )
call create_pressure_solver( pressure_operator, pressure_preconditioner, pressure_solver )
call create_mixed_preconditioner( rhs, pressure_solver, mixed_preconditioner )
call create_mixed_preconditioner( rhs, pressure_solver, mixed_preconditioner, stepper_siqn )

call create_adj_pressure_preconditioner( rhs, adj_lookup_table_cache, &
adj_pressure_operator, adj_pressure_preconditioner )
adj_pressure_operator, adj_pressure_preconditioner, &
stepper_siqn )
call create_adj_pressure_solver( adj_pressure_operator, adj_pressure_preconditioner, adj_pressure_solver )
call create_adj_mixed_preconditioner( rhs, adj_pressure_solver, adj_mixed_preconditioner )
call create_adj_mixed_preconditioner( rhs, adj_pressure_solver, adj_mixed_preconditioner, stepper_siqn )

call construct_solver_state( vector_x, rhs_rsol, import_fields=.true. )
call construct_solver_state( vector_mx, rhs_rsol, import_fields=.false. )
Expand Down
22 changes: 13 additions & 9 deletions applications/adjoint_tests/source/algorithm/solver/adjt_mixed_solver_alg_mod.x90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -34,7 +34,9 @@ module adjt_mixed_solver_alg_mod
use driver_modeldb_mod, only: modeldb_type
use solver_constants_mod, only: get_normalisation
use fs_continuity_mod, only: W2, W3, Wtheta
use si_operators_alg_mod, only: compute_si_operators
use si_operators_alg_mod, only: compute_si_operators_common, &
compute_helmholtz_operators
use dycore_constants_mod, only: stepper_siqn
use operator_mod, only: operator_type
use sci_enforce_bc_kernel_mod, only: enforce_bc_kernel_type
use model_clock_mod, only: model_clock_type
Expand Down Expand Up @@ -189,19 +191,21 @@ contains
call moisture_fields%get_field( "ls_moist_dyn", ls_moist_dyn_array )
ls_moist_dyn => ls_moist_dyn_array%bundle

call compute_si_operators( modeldb%config, ls_theta, ls_rho, ls_exner, &
model_clock, ls_moist_dyn )
call compute_si_operators_common( modeldb%config, ls_theta, ls_rho, ls_exner, &
ls_moist_dyn )
call compute_helmholtz_operators( model_clock, stepper_siqn )

call create_pressure_preconditioner( rhs, pressure_operator, pressure_preconditioner )
call create_pressure_preconditioner( rhs, pressure_operator, pressure_preconditioner, stepper_siqn )
call create_pressure_solver( pressure_operator, pressure_preconditioner, pressure_solver )
call create_mixed_preconditioner( rhs, pressure_solver, mixed_preconditioner )
call create_mixed_solver( mixed_preconditioner, mixed_operator, mixed_solver )
call create_mixed_preconditioner( rhs, pressure_solver, mixed_preconditioner, stepper_siqn )
call create_mixed_solver( mixed_preconditioner, mixed_operator, mixed_solver, stepper_siqn )

call create_adj_pressure_preconditioner( rhs, adj_lookup_table_cache, &
adj_pressure_operator, adj_pressure_preconditioner )
adj_pressure_operator, adj_pressure_preconditioner, &
stepper_siqn )
call create_adj_pressure_solver( adj_pressure_operator, adj_pressure_preconditioner, adj_pressure_solver )
call create_adj_mixed_preconditioner( rhs, adj_pressure_solver, adj_mixed_preconditioner )
call create_adj_mixed_solver( adj_mixed_preconditioner, adj_mixed_operator, adj_mixed_solver )
call create_adj_mixed_preconditioner( rhs, adj_pressure_solver, adj_mixed_preconditioner, stepper_siqn )
call create_adj_mixed_solver( adj_mixed_preconditioner, adj_mixed_operator, adj_mixed_solver, stepper_siqn )

call construct_solver_state( vector_x, rhs_rsol, import_fields=.true. )
call construct_solver_state( vector_mx, rhs_rsol, import_fields=.false. )
Expand Down
20 changes: 12 additions & 8 deletions applications/adjoint_tests/source/algorithm/solver/adjt_pressure_precon_alg_mod.x90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -27,7 +27,9 @@ module adjt_pressure_precon_alg_mod
use function_space_collection_mod, only: function_space_collection
use driver_modeldb_mod, only: modeldb_type
use fs_continuity_mod, only: W2, W3, Wtheta
use si_operators_alg_mod, only: compute_si_operators
use si_operators_alg_mod, only: compute_si_operators_common, &
compute_helmholtz_operators
use dycore_constants_mod, only: stepper_siqn
use sci_enforce_bc_kernel_mod, only: enforce_bc_kernel_type
use model_clock_mod, only: model_clock_type
use function_space_mod, only: function_space_type
Expand Down Expand Up @@ -149,15 +151,17 @@ contains
call moisture_fields%get_field( "ls_moist_dyn", ls_moist_dyn_array )
ls_moist_dyn => ls_moist_dyn_array%bundle

call compute_si_operators( modeldb%config, ls_theta, ls_rho, ls_exner, &
model_clock, ls_moist_dyn )
call compute_si_operators_common( modeldb%config, ls_theta, ls_rho, ls_exner, &
ls_moist_dyn )
call compute_helmholtz_operators( model_clock, stepper_siqn )

call create_pressure_preconditioner( rhs, pressure_operator, pressure_preconditioner )
call create_pressure_preconditioner( rhs, pressure_operator, pressure_preconditioner, stepper_siqn )

call create_adj_pressure_preconditioner( rhs, &
adj_lookup_table_cache, &
adj_pressure_operator, &
adj_pressure_preconditioner )
call create_adj_pressure_preconditioner( rhs, &
adj_lookup_table_cache, &
adj_pressure_operator, &
adj_pressure_preconditioner, &
stepper_siqn )

! Size = 1 as only P field needed
vector_x = r_solver_field_vector_type(1_i_def)
Expand Down
25 changes: 19 additions & 6 deletions applications/adjoint_tests/source/algorithm/solver/adjt_scaled_matrix_vector_alg_mod.x90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,9 @@
! (C) Crown copyright 2025 Met Office. All rights reserved.
! The file LICENCE, distributed with this code, contains details of the terms
! under which the code may be used.
!
! Some of the content of this file has been produced with the assistance of
! Met Office GitHub Copilot Enterprise.
!-----------------------------------------------------------------------------

! @todo Causes failure of later test (mixed_schur_preconditioner); not calling until resolved (#616).
Expand All @@ -23,7 +26,11 @@ module adjt_scaled_matrix_vector_alg_mod
use operator_mod, only: r_solver_operator_type
use r_solver_field_mod, only: r_solver_field_type
use sci_r_solver_field_vector_mod, only: r_solver_field_vector_type
use si_operators_alg_mod, only: compute_si_operators, get_div_star, get_Hb_lumped_inv
use si_operators_alg_mod, only: compute_si_operators_common, &
compute_helmholtz_operators, &
get_div_star, get_Hb_lumped_inv, &
get_tau_u_dt_cp
use dycore_constants_mod, only: stepper_siqn
use solver_constants_mod, only: get_normalisation_r_solver

! Object types
Expand Down Expand Up @@ -83,6 +90,7 @@ module adjt_scaled_matrix_vector_alg_mod
real(kind=r_def) :: machine_tolerance
real(kind=r_def) :: relative_diff
real(kind=r_def), parameter :: overall_tolerance = 1500.0_r_def
real(kind=r_solver) :: const_u

! Setup

Expand All @@ -106,24 +114,28 @@ module adjt_scaled_matrix_vector_alg_mod
call invoke(setval_random(init_op_moist_dyn(idx)))
end do

call compute_si_operators( config, init_op_t, init_op_d, init_op_p, &
model_clock, init_op_moist_dyn )
call compute_si_operators_common( config, init_op_t, init_op_d, init_op_p, &
init_op_moist_dyn )
call compute_helmholtz_operators( model_clock, stepper_siqn )

call vector_x%initialise(vector_space_w3_ptr)
call vector_mx%initialise(vector_space_w2_ptr)
call vector_x%copy_field_properties(vector_amx)

u_normalisation => get_normalisation_r_solver( W2, vector_mx%get_mesh_id() )
div_star => get_div_star()
Hb_lumped_inv => get_Hb_lumped_inv()
Hb_lumped_inv => get_Hb_lumped_inv(stepper_siqn)

const_u = get_tau_u_dt_cp(stepper_siqn)

call invoke( setval_random(vector_x), &
setval_c( vector_mx, 0.0_r_solver ), &
setval_c( vector_amx, 0.0_r_solver ), &

! Tangent linear (Mx)

scaled_matrix_vector_kernel_type( vector_mx, vector_x, div_star, u_normalisation, Hb_lumped_inv ), &
scaled_matrix_vector_kernel_type( vector_mx, vector_x, div_star, u_normalisation, Hb_lumped_inv, &
const_u ), &

! <Mx,Mx>

Expand All @@ -143,7 +155,8 @@ module adjt_scaled_matrix_vector_alg_mod

! Adjoint (AMx)

adj_scaled_matrix_vector_kernel_type( vector_mx, vector_amx, div_star, u_normalisation, Hb_lumped_inv ), &
adj_scaled_matrix_vector_kernel_type( vector_mx, vector_amx, div_star, u_normalisation, &
Hb_lumped_inv, const_u ), &

! <AMx,x>

Expand Down
13 changes: 8 additions & 5 deletions ...ications/adjoint_tests/source/algorithm/solver/adjt_semi_implicit_solver_step_alg_mod.x90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -28,7 +28,9 @@ module adjt_semi_implicit_solver_step_alg_mod
use driver_modeldb_mod, only: modeldb_type
use fs_continuity_mod, only: W2, W3, Wtheta
use mr_indices_mod, only: nummr
use si_operators_alg_mod, only: compute_si_operators
use si_operators_alg_mod, only: compute_si_operators_common, &
compute_helmholtz_operators
use dycore_constants_mod, only: stepper_siqn
use operator_mod, only: operator_type
use sci_enforce_bc_kernel_mod, only: enforce_bc_kernel_type
use model_clock_mod, only: model_clock_type
Expand Down Expand Up @@ -163,9 +165,10 @@ contains
call moisture_fields%get_field( "ls_moist_dyn", ls_moist_dyn_array )
ls_moist_dyn => ls_moist_dyn_array%bundle

call compute_si_operators( modeldb%config, ls_theta, ls_rho, ls_exner, &
model_clock, ls_moist_dyn )
call adj_semi_implicit_solver%initialise(rhs, adj_lookup_table_cache)
call compute_si_operators_common( modeldb%config, ls_theta, ls_rho, ls_exner, &
ls_moist_dyn )
call compute_helmholtz_operators( model_clock, stepper_siqn )
call adj_semi_implicit_solver%initialise(rhs, adj_lookup_table_cache, stepper_siqn)

! Set up moist_dyn_gas_law & mr
call moist_dyn_gas_law%initialise(vector_space=vector_space_wtheta_ptr)
Expand Down Expand Up @@ -288,7 +291,7 @@ contains
! --------------------------------------------------------------------------

call semi_implicit_solver_alg_step( state, rhs, dry_flux_solver, &
moist_dyn_gas_law, mr, first_iteration )
moist_dyn_gas_law, mr, first_iteration, stepper_siqn )

! --------------------------------------------------------------------------
! Compute <Mx,Mx>
Expand Down
11 changes: 7 additions & 4 deletions applications/adjoint_tests/source/algorithm/timestepping/atlt_si_timestep_alg_mod.x90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -35,7 +35,9 @@ module atlt_si_timestep_alg_mod
tl_semi_implicit_alg_final
use atl_si_timestep_alg_mod, only: atl_si_timestep_type
use driver_modeldb_mod, only: modeldb_type
use si_operators_alg_mod, only: compute_si_operators
use si_operators_alg_mod, only: compute_si_operators_common, &
compute_helmholtz_operators
use dycore_constants_mod, only: stepper_siqn
use sci_enforce_bc_kernel_mod, only: enforce_bc_kernel_type
use field_array_mod, only: field_array_type

Expand Down Expand Up @@ -185,10 +187,11 @@ contains

call tl_semi_implicit_alg_init( mesh, u, rho, theta, exner, mr, &
ls_u, ls_rho, ls_theta, ls_exner, &
ls_mr, ls_moist_dyn, (num_steps == 1) )
ls_mr, ls_moist_dyn, modeldb, (num_steps == 1) )

call compute_si_operators( modeldb%config, ls_theta, ls_rho, ls_exner, &
modeldb%clock, ls_moist_dyn )
call compute_si_operators_common( modeldb%config, ls_theta, ls_rho, ls_exner, &
ls_moist_dyn )
call compute_helmholtz_operators( modeldb%clock, stepper_siqn )

do i = 1, num_steps
call tl_semi_implicit_alg_step( modeldb, u, rho, theta, exner, mr, moist_dyn, &
Expand Down
Loading
Loading