Skip to content

Allow substepping of just tracers #733

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 8 commits into
base: main
Choose a base branch
from
Open

Allow substepping of just tracers #733

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 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
0fb29a0
bring over changes to implement substepping for tracers
tommbendall Aug 24, 2026
8143cc1
only substep tracers if namelist says so!
tommbendall Aug 25, 2026
351c569
tighten up substepping criteria
tommbendall Aug 25, 2026
eba60b0
fixes: ensure advective transport kicks in correctly, clear shared
tommbendall Sep 24, 2026
ba71df5
fixes: use correct dt for share stencil extent, take substepping into
tommbendall Sep 24, 2026
5000f7f
protect total_ref_flux so that it is correct when set
tommbendall Sep 28, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
27 changes: 23 additions & 4 deletions science/gungho/source/algorithm/transport/common/flux_precomputations_mod.x90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -107,8 +107,9 @@ module flux_precomputations_alg_mod
integer(kind=i_def) :: max_stencil_extent
logical(kind=l_def) :: max_stencil_extent_computed = .false.

! Total reference flux (over whole transport step), needed for cheap update
type(r_tran_field_type) :: total_ref_flux
! Total reference flux (over whole transport step), needed when there is
! a separate tracer transport, e.g. cheap update, coarse tracers, TR-BDF2
type(r_tran_field_type), allocatable :: total_ref_flux

! Wind and dt for whole transport step, needed for logging
type(r_tran_field_type), pointer :: transporting_wind => null()
Expand Down Expand Up @@ -294,6 +295,7 @@ contains
self%to_initialise_by_step = present(flux)
self%reset_unity = .false.
if (present(reset_unity)) self%reset_unity = reset_unity
allocate(self%total_ref_flux)

! Determine settings for extended mesh remapping of density
if (panel_edge_treatment == panel_edge_treatment_extended_mesh) then
Expand Down Expand Up @@ -603,6 +605,7 @@ contains
nullify(self%transporting_wind)
self%is_initialised_flag = .false.
self%max_stencil_extent_computed = .false.
if ( allocated( self%total_ref_flux ) ) deallocate( self%total_ref_flux )
if ( allocated( self%mesh_ids ) ) deallocate( self%mesh_ids )
if ( allocated( self%ref_field ) ) deallocate( self%ref_field )
if ( allocated( self%ref_mass ) ) deallocate( self%ref_mass )
Expand Down Expand Up @@ -1758,7 +1761,8 @@ contains
!! reference field over the *whole* transport step. This is needed when
!! using the "cheap update" transport formulation, as the total flux
!! needs storing from previous outer loops for the transport of tracers
!! in the final outer loop.
!! in the final outer loop. It is also needed for coarse transport of
!! tracers, and the transport of tracers in TR-BDF2.
!> @return The total reference flux
function get_total_ref_flux(self) result(total_ref_flux)

Expand All @@ -1768,6 +1772,13 @@ contains
type(r_tran_field_type), pointer :: total_ref_flux
character(len=str_def) :: field_name

if (.not. allocated(self%total_ref_flux)) then
call log_event( &
'flux_precomp%get_total_ref_flux: self%total_ref_flux not allocated', &
LOG_LEVEL_ERROR &
)
end if

total_ref_flux => self%total_ref_flux

! Check field has been initialised
Expand All @@ -1780,7 +1791,8 @@ contains
!! field over the *whole* transport step. This is needed when using the
!! "cheap update" transport formulation, as the total flux needs
!! storing from previous outer loops for the transport of tracers in
!! the final outer loop.
!! the final outer loop. It is also needed for coarse transport of
!! tracers, and the transport of tracers in TR-BDF2.
!! As we may be substepping, the field is incremented rather than set
!> @param[in] total_ref_flux The total reference flux to set
subroutine set_total_ref_flux(self, total_ref_flux)
Expand All @@ -1790,6 +1802,13 @@ contains
class(flux_precomputations_type), target, intent(inout) :: self
type(r_tran_field_type), intent(in) :: total_ref_flux

if (.not. allocated(self%total_ref_flux)) then
call log_event( &
'flux_precomp%set_total_ref_flux: self%total_ref_flux not allocated', &
LOG_LEVEL_ERROR &
)
end if

if (.not. self%total_ref_flux%is_initialised()) then
call self%total_ref_flux%initialise(total_ref_flux%get_function_space())
call invoke( setval_c(self%total_ref_flux, 0.0_r_tran) )
Expand Down
3 changes: 2 additions & 1 deletion science/gungho/source/algorithm/transport/control/gungho_transport_control_alg_mod.X90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -244,6 +244,7 @@ contains
! Internal variables
logical(kind=l_def) :: do_moisture_diagnostics
logical(kind=l_def) :: cheap_update_step
logical(kind=l_def) :: safe_ref

! Temporary fields or pointers
type(field_type) :: fields_np1(bundle_size)
Expand Down Expand Up @@ -346,7 +347,7 @@ contains

! Check negative reference fields at this point
if (check_any_eqn_consistent()) then
call transport_controller%check_reference_field()
safe_ref = transport_controller%check_reference_field(fail_on_neg=.true.)
end if

! ------------------------------------------------------------------------ !
Expand Down
98 changes: 84 additions & 14 deletions science/gungho/source/algorithm/transport/control/transport_controller_mod.x90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -396,7 +396,8 @@ contains
watkins_fail_field, &
wind_flux, &
dt_one, &
detj_at_w3), &
detj_at_w3, &
1), &
! Determine the adjusted wind for the final step
X_minus_Y(step_3_wind_flux, vert_wind_flux, step_1_wind_flux) )

Expand Down Expand Up @@ -557,7 +558,10 @@ contains
!> @brief Calculates the reference fields to be used in tracer transport, and
!! if any of these are going to be negative, may perform the Watkins
!! algorithm to adjust how the fluxes are distributed between steps.
subroutine check_reference_field(self)
!> @param[in] fail_on_neg If true, the model will stop if a reference field
!! is negative. If false, the model will continue.
!> @return True if the ref field is safe to use, false if not
function check_reference_field(self, fail_on_neg) result(ref_field_safe)

use transport_config_mod, only: ffsl_unity_3d, &
adjust_vhv_wind, &
Expand All @@ -572,6 +576,7 @@ contains
implicit none

class(transport_controller_type), target, intent(inout) :: self
logical(kind=l_def), intent(in) :: fail_on_neg

type(flux_precomputations_type), pointer :: flux_precomputations

Expand Down Expand Up @@ -668,7 +673,8 @@ contains
watkins_fail_field, &
ref_flux_one_step, &
dt_one, &
ref_mass), &
ref_mass, &
self%num_substeps), &
! Determine the adjusted wind for the final step
X_minus_Y(step_3_flux, vert_flux, step_1_flux) )

Expand All @@ -688,22 +694,37 @@ contains
call self%flux_precomputations%initialise_step(step=(3+3*(i-1)), ref_flux=step_3_flux)
end do

ref_field_safe = .true.
else
! We're about to fail... log the Lipschitz numbers
ref_field_safe = .false.

call lipschitz_diag_alg( &
self%config, &
self%reference_splitting, self%transporting_wind, &
self%dt_substep, LOG_LEVEL_WARNING &
)

if (adjust_tracer_equation) then
if (.not. fail_on_neg &
.and. substep_transport == substep_transport_adaptive &
.and. self%num_substeps < 3) then
! Time to attempt substepping...
call log_event( &
'The Watkins algorithm failed to keep all reference fields ' // &
'positive for the transport of tracers. Now try substepping...', &
LOG_LEVEL_WARNING &
)

else if (adjust_tracer_equation) then
! Use advective form for tracer transport instead of consistent form
call log_event( &
'Reference field for tracer transport is negative, so ' // &
'making all tracers advective for this step. This will ' // &
'make tracer transport non-conservative. ', LOG_LEVEL_INFO &
'make tracer transport non-conservative. ', LOG_LEVEL_WARNING &
)
self%make_tracers_advective = .true.
ref_field_safe = .true.

else
! Fail gracefully
call log_event( &
Expand Down Expand Up @@ -734,7 +755,7 @@ contains
end if
end if

end subroutine check_reference_field
end function check_reference_field

! ============================================================================ !
! SETTERS
Expand Down Expand Up @@ -864,14 +885,16 @@ contains
type(flux_precomputations_type), pointer :: ffsl_precomputations
type(mesh_type), pointer :: mesh
type(r_tran_field_type) :: wind_flux
real(kind=r_tran) :: dt_whole_step

! If object is not created, then create it here
if (.not. self%ffsl_precomputations%is_initialised()) then
mesh => self%transporting_wind%get_mesh()

! Wind needs multiplying by dt
dt_whole_step = self%dt_substep * real(self%num_substeps, r_tran)
call self%transporting_wind%copy_field_properties(wind_flux)
call invoke( a_times_X(wind_flux, self%dt_substep, self%transporting_wind) )
call invoke( a_times_X(wind_flux, dt_whole_step, self%transporting_wind) )

! Since the object has not already been created, we initialise it with the
! 3D wind, rather than a wind flux for each split step
Expand Down Expand Up @@ -1021,12 +1044,17 @@ contains
!> @param[in] ref_flux The reference flux field
subroutine initialise_flux_precomputations(self, ref_flux)

use transport_config_mod, only: substep_transport, &
substep_transport_adaptive

implicit none

class(transport_controller_type), target, intent(inout) :: self
type(r_tran_field_type), intent(in) :: ref_flux

type(mesh_type), pointer :: mesh
logical(kind=l_def) :: safe_ref
integer(kind=i_def) :: i, start_substeps, end_substeps

if (self%flux_precomputations%is_initialised()) then
call log_event( &
Expand All @@ -1043,13 +1071,55 @@ contains
! coarser mesh. In either case, the fluxes calculated during the transport
! of the density are not the same as those that will be used to transport
! the tracers. We set the correct flux here.
call self%flux_precomputations%initialise( &
self%config, &
mesh, self%reference_splitting, self%num_substeps, &
self%transporting_wind, self%dt_substep, &
flux=ref_flux, ref_field=self%ref_field_rtran &
)
call self%check_reference_field()

! Loop through substeps, as we may need to increase the number to keep
! the reference fields positive
safe_ref = .false.
start_substeps = self%num_substeps
if (substep_transport == substep_transport_adaptive) then
end_substeps = 3 ! Force max substeps of 3
else
end_substeps = 1 ! No substepping
end if

do i = start_substeps, end_substeps
call self%flux_precomputations%initialise( &
self%config, &
mesh, self%reference_splitting, self%num_substeps, &
self%transporting_wind, self%dt_substep, &
flux=ref_flux, ref_field=self%ref_field_rtran &
)
if (i == end_substeps) then
! Fail if we have reached max substeps
safe_ref = self%check_reference_field(fail_on_neg=.true.)
else
safe_ref = self%check_reference_field(fail_on_neg=.false.)
end if

if (.not. safe_ref) then
! Watkins algorithm failed to keep reference field positive
! This is likely to be BDF tracer transport near start of run
! Try enforcing substepping for tracer transport
self%dt_substep = ( &
self%dt_substep * self%num_substeps &
/ real(self%num_substeps + 1, r_tran) &
)
self%num_substeps = self%num_substeps + 1

! Clear all precomputations in case they were set
call self%ffsl_precomputations%finalise()
call self%flux_precomputations%finalise()
call self%wind_precomputations%finalise()
self%dep_stencil_extent_computed = .false.

write(log_scratch_space, '(A,I4)') &
'Transport: Reference field for tracer transport cannot stay ' // &
'positive. Increasing number of substeps to ', self%num_substeps
call log_event(log_scratch_space, LOG_LEVEL_WARNING)
else
EXIT
end if
end do

end subroutine initialise_flux_precomputations

Expand Down
Loading
Loading