diff --git a/science/gungho/source/algorithm/transport/common/flux_precomputations_mod.x90 b/science/gungho/source/algorithm/transport/common/flux_precomputations_mod.x90 index eada1ca879..5de81674e7 100644 --- a/science/gungho/source/algorithm/transport/common/flux_precomputations_mod.x90 +++ b/science/gungho/source/algorithm/transport/common/flux_precomputations_mod.x90 @@ -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() @@ -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 @@ -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 ) @@ -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) @@ -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 @@ -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) @@ -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) ) diff --git a/science/gungho/source/algorithm/transport/control/gungho_transport_control_alg_mod.X90 b/science/gungho/source/algorithm/transport/control/gungho_transport_control_alg_mod.X90 index f799cef47f..4656203b1c 100644 --- a/science/gungho/source/algorithm/transport/control/gungho_transport_control_alg_mod.X90 +++ b/science/gungho/source/algorithm/transport/control/gungho_transport_control_alg_mod.X90 @@ -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) @@ -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 ! ------------------------------------------------------------------------ ! diff --git a/science/gungho/source/algorithm/transport/control/transport_controller_mod.x90 b/science/gungho/source/algorithm/transport/control/transport_controller_mod.x90 index fe2eaeeb55..f0d6f6b7ec 100644 --- a/science/gungho/source/algorithm/transport/control/transport_controller_mod.x90 +++ b/science/gungho/source/algorithm/transport/control/transport_controller_mod.x90 @@ -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) ) @@ -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, & @@ -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 @@ -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) ) @@ -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( & @@ -734,7 +755,7 @@ contains end if end if - end subroutine check_reference_field + end function check_reference_field ! ============================================================================ ! ! SETTERS @@ -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 @@ -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( & @@ -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 diff --git a/science/gungho/source/kernel/transport/common/watkins_kernel_mod.F90 b/science/gungho/source/kernel/transport/common/watkins_kernel_mod.F90 index 2aacc73f3c..a2c95a9494 100644 --- a/science/gungho/source/kernel/transport/common/watkins_kernel_mod.F90 +++ b/science/gungho/source/kernel/transport/common/watkins_kernel_mod.F90 @@ -31,13 +31,14 @@ module watkins_kernel_mod !> The type declaration for the kernel. Contains the metadata needed by the Psy layer type, public, extends(kernel_type) :: watkins_kernel_type private - type(arg_type) :: meta_args(6) = (/ & + type(arg_type) :: meta_args(7) = (/ & arg_type(GH_FIELD, GH_REAL, GH_WRITE, W2v), & ! first_v_wind arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_2), & ! lipschitz_max_field arg_type(GH_FIELD, GH_INTEGER, GH_WRITE, ANY_DISCONTINUOUS_SPACE_2), & ! watkins_failures arg_type(GH_FIELD, GH_REAL, GH_READ, W2), & ! wind arg_type(GH_SCALAR, GH_REAL, GH_READ), & ! dt - arg_type(GH_FIELD, GH_REAL, GH_READ, W3) & ! detj + arg_type(GH_FIELD, GH_REAL, GH_READ, W3), & ! detj + arg_type(GH_SCALAR, GH_INTEGER, GH_READ) & ! num_substeps /) integer :: operates_on = CELL_COLUMN contains @@ -60,6 +61,7 @@ module watkins_kernel_mod !> @param[in] wind 3D wind !> @param[in] dt Transport time step !> @param[in] detj Det(J) at W3: the volume of cells +!> @param[in] num_substeps Maximum number of transport substeps !> @param[in] ndf_w2v Num of DoFs per cell for W2V !> @param[in] undf_w2v Num of W2V DoFs for this partition !> @param[in] map_w2v Map of lowest-cell W2V DoFs @@ -79,6 +81,7 @@ subroutine watkins_code( nlayers, & wind, & dt, & detj, & + num_substeps, & ndf_w2v, & undf_w2v, & map_w2v, & @@ -102,6 +105,7 @@ subroutine watkins_code( nlayers, & integer(kind=i_def), intent(in) :: map_w2v(ndf_w2v) integer(kind=i_def), intent(in) :: map_w3(ndf_w3) integer(kind=i_def), intent(in) :: map_w3_2d(ndf_w3_2d) + integer(kind=i_def), intent(in) :: num_substeps real(kind=r_tran), intent(in) :: dt real(kind=r_tran), intent(in) :: detj(undf_w3) real(kind=r_tran), intent(in) :: wind(undf_w2) @@ -115,12 +119,14 @@ subroutine watkins_code( nlayers, & real(kind=r_tran) :: lip_x_kp1, lip_y_kp1, lip_hori_kp1, lip_z_half_kp1 real(kind=r_tran) :: lip_lim_km1, lip_lim, lip_lim_kp1, slack_km1, remainder real(kind=r_tran) :: max_lip_init, max_lip_new + real(kind=r_tran) :: threshold, threshold_kp1 + real(kind=r_tran) :: lip_3d_kp1, lip_3d, lip_3d_km1 ! The threshold sets the maximum allowed Lipschitz number under this scheme, ! while the tolerance is used for checking whether the algorithm has failed, ! and this is slacker than machine precision to avoid machine precision errors ! being detected as failures. - real(kind=r_tran), parameter :: threshold = 0.9_r_tran + real(kind=r_tran), parameter :: max_threshold = 0.9_r_tran real(kind=r_tran), parameter :: tolerance = 1.0E-3_r_tran ! Set failures to zero @@ -154,8 +160,10 @@ subroutine watkins_code( nlayers, & lip_z_half_kp1 = dt*(first_v_wind(map_w2v(2)) - first_v_wind(map_w2v(1)))/detj(map_w3(1)) lip_hori_kp1 = MAX(lip_x_kp1, lip_y_kp1, lip_x_kp1 + lip_y_kp1) lip_lim_kp1 = MAX(lip_z_half_kp1, lip_z_half_kp1 + lip_hori_kp1) - max_lip_init = lip_lim_kp1 + lip_3d_kp1 = MAX(0.0_r_tran, lip_x_kp1 + lip_y_kp1 + 2.0_r_tran*lip_z_half_kp1) + max_lip_init = lip_lim_kp1 + REAL(num_substeps - 1, r_tran)*lip_3d_kp1 max_lip_new = 0.0_r_tran + threshold_kp1 = max_threshold - REAL(num_substeps - 1, r_tran)*lip_3d_kp1 do k = 0, nlayers - 1 @@ -167,6 +175,8 @@ subroutine watkins_code( nlayers, & if (k > 0) then lip_hori_km1 = lip_hori + lip_3d_km1 = lip_3d + threshold = max_threshold - REAL(num_substeps - 1, r_tran)*lip_3d_km1 lip_z_half_km1 = dt*(first_v_wind(map_w2v(2)+k-1) - first_v_wind(map_w2v(1)+k-1))/detj(map_w3(1)+k-1) lip_lim_km1 = MAX(lip_z_half_km1, lip_z_half_km1 + lip_hori_km1) @@ -176,8 +186,12 @@ subroutine watkins_code( nlayers, & ! Lipschitz numbers for this cell ------------------------------------------ lip_hori = lip_hori_kp1 + lip_3d = lip_3d_kp1 lip_z_half = dt*(first_v_wind(map_w2v(2)+k) - first_v_wind(map_w2v(1)+k))/detj(map_w3(1)+k) lip_lim = MAX(lip_z_half, lip_z_half + lip_hori) + ! Anticipate subsequent substeps, so build 3D Lipschitz number into threshold + ! Each substep will add lip_3d to the Lipschitz number + threshold = max_threshold - REAL(num_substeps - 1, r_tran)*lip_3d ! Lipschitz numbers for cell above ----------------------------------------- if (k < nlayers - 1) then @@ -186,7 +200,14 @@ subroutine watkins_code( nlayers, & lip_z_half_kp1 = dt*(first_v_wind(map_w2v(2)+k+1) - first_v_wind(map_w2v(1)+k+1))/detj(map_w3(1)+k+1) lip_hori_kp1 = MAX(lip_x_kp1, lip_y_kp1, lip_x_kp1 + lip_y_kp1) lip_lim_kp1 = MAX(lip_z_half_kp1, lip_z_half_kp1 + lip_hori_kp1) - max_lip_init = MAX(max_lip_init, lip_lim_kp1) + lip_3d_kp1 = MAX(0.0_r_tran, & + lip_x_kp1 + lip_y_kp1 & + + dt*(wind(map_w2(T)+k+1) - wind(map_w2(B)+k+1))/detj(map_w3(1)+k+1) & + ) + ! For tracking if max Lipschitz number has increased + max_lip_init = MAX(max_lip_init, & + lip_lim_kp1 + REAL(num_substeps-1, r_tran)*lip_3d_kp1) + threshold_kp1 = max_threshold - REAL(num_substeps - 1, r_tran)*lip_3d_kp1 end if ! ------------------------------------------------------------------------ ! @@ -201,7 +222,7 @@ subroutine watkins_code( nlayers, & if (slack_km1 > (lip_lim - threshold) * detj(map_w3(1)+k) / dt) then ! Check whether cell above might need some of this slack too - if (lip_lim_kp1 > threshold .AND. k < nlayers - 1) then + if (lip_lim_kp1 > threshold_kp1 .AND. k < nlayers - 1) then ! TODO: don't need to use all slack here! Arbitrarily using half first_v_wind(map_w2v(1)+k) = first_v_wind(map_w2v(1)+k) & + (lip_lim - threshold) * detj(map_w3(1)+k) / dt & @@ -229,7 +250,8 @@ subroutine watkins_code( nlayers, & ! Calculate Lipschitz number from cell below (which now won't change) if (k > 0) then lip_z_half_km1 = dt*(first_v_wind(map_w2v(2)+k-1) - first_v_wind(map_w2v(1)+k-1))/detj(map_w3(1)+k-1) - lip_lim_km1 = MAX(lip_z_half_km1, lip_z_half_km1 + lip_hori_km1) + lip_lim_km1 = MAX(lip_z_half_km1, lip_z_half_km1 + lip_hori_km1) & + + REAL(num_substeps - 1, r_tran)*lip_3d_km1 max_lip_new = MAX(max_lip_new, lip_lim_km1) end if @@ -248,14 +270,14 @@ subroutine watkins_code( nlayers, & k = nlayers - 1 lip_z_half = dt*(first_v_wind(map_w2v(2)+k) - first_v_wind(map_w2v(1)+k))/detj(map_w3(1)+k) lip_lim = MAX(lip_z_half, lip_z_half + lip_hori) - max_lip_new = MAX(max_lip_new, lip_lim) + max_lip_new = MAX(max_lip_new, lip_lim + REAL(num_substeps - 1, r_tran)*lip_3d) ! -------------------------------------------------------------------------- ! ! Check if algorithm has failed, and if so take original wind ! -------------------------------------------------------------------------- ! if (max_lip_new > max_lip_init + tolerance .AND. & - max_lip_new > threshold + tolerance) then + max_lip_new > max_threshold + tolerance) then watkins_failures(map_w3_2d(1)) = 1_i_def ! Set the bottom value @@ -269,7 +291,7 @@ subroutine watkins_code( nlayers, & ! Set the top values first_v_wind(map_w2v(1)+nlayers) = 0.0_r_tran - else if (max_lip_new > threshold + tolerance) then + else if (max_lip_new > max_threshold + tolerance) then watkins_failures(map_w3_2d(1)) = 1_i_def end if diff --git a/science/gungho/unit-test/kernel/transport/common/watkins_kernel_mod_test.pf b/science/gungho/unit-test/kernel/transport/common/watkins_kernel_mod_test.pf index 88956d7366..296fc5f928 100644 --- a/science/gungho/unit-test/kernel/transport/common/watkins_kernel_mod_test.pf +++ b/science/gungho/unit-test/kernel/transport/common/watkins_kernel_mod_test.pf @@ -160,6 +160,7 @@ contains wind, & dt, & detj_at_w3, & + 1, & ndf_w2v, & undf_w2v, & map_w2v, & @@ -261,6 +262,7 @@ contains wind, & dt, & detj_at_w3, & + 1, & ndf_w2v, & undf_w2v, & map_w2v, &