Linear and adjoint boundary layer physics #182

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.

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

tom-j-h wants to merge 118 commits into MetOffice:main from tom-j-h:tlad_boundary_layer

applications/adjoint_tests/example/configuration.nml

            
                      Original file line number
                      Diff line number
                      Diff line change
                  
    @@ -68,8 +68,8 @@ stretching_method='smooth',
  
    ancil_directory='/data/users/lfricadmin/data/ancils/basic-gal/yak/C12',

    checkpoint_stem_name='',

    diag_stem_name='diagGungho',

    ls_directory='/data/users/lfricadmin/data/tangent-linear/Ticket354',

    ls_filename='final_ls',

    ls_directory='/data/users/tim.payne/lfric_apps/files',

    ls_filename='final_ls_with_land',

    orography_mean_ancil_path='orography/gmted_ramp2/qrparm.orog',

    start_dump_directory='/data/users/lfricadmin/data/tangent-linear/Ticket354',

    start_dump_filename='final_pert',

    @@ -221,12 +221,23 @@ write_minmax_tseries=.false.,
  
    /

    &linear

    fixed_ls=.true.,

    l_stabilise_bl=.true.,

    l_stabilise_bl=.false.,

    ls_read_w2h=.false.,

    max_bl_stabilisation=0.75,

    n_bl_levels_to_stabilise=15,

    pert_option='file',

    /

    &linear_physics

    blevs_m=15,

    e_folding_levs_m=10,

    l_0_m=80.0,

    l_boundary_layer=.true.,

    log_layer=2,

    u_land_m=0.4,

    u_sea_m=0.4,

    z_land_m=0.05,

    z_sea_m=0.0005,

    /

    &logging

    log_to_rank_zero_only=.false.,

    run_log_level='debug',

applications/adjoint_tests/example/iodef.xml

-Original file line number
+Diff line change
@@ Expand Up / @@ -50,6 +50,11 @@ @@
           <axis axis_ref="vert_axis_half_levels"/>
           <axis axis_ref="ls_axis"/>
         </grid>
+        <grid id="one_level_face_time_grid">
+          <domain domain_ref="face"/>
+          <axis axis_ref="1"/>
+          <axis axis_ref="ls_axis"/>
+        </grid>
        </grid_definition>
        <field_definition prec="8" >
@@ Expand Down Expand Up / @@ -318,6 +323,7 @@ @@
             <field id="ls_u_in_w3_data" name="u_in_w3" long_name="NS wind" unit="m s-1" grid_ref="half_level_face_time_grid" operation="once"/>
             <field id="ls_v_in_w3_data" name="v_in_w3" long_name="EW wind" unit="m s-1" grid_ref="half_level_face_time_grid" operation="once"/>
             <field id="ls_w_in_wth_data" name="w_in_wth" long_name="Upward wind" unit="m s-1" grid_ref="full_level_face_time_grid" operation="once"/>
+    	<field id="ls_land_fraction_data" name="land_area_fraction" long_name="land_area_fraction" unit="1" grid_ref="one_level_face_time_grid" operation="once" />
           </file>
         </file_definition>
@@ Expand Down @@

applications/adjoint_tests/source/algorithm/linear_physics/atlt_bdy_lyr_alg_mod.x90

-Original file line number
+Diff line change
@@ -0,0 +1,192 @@
+    !-----------------------------------------------------------------------------
+    ! (c) Crown copyright 2026 Met Office. All rights reserved.
+    ! The file LICENCE, distributed with this code, contains details of the terms
+    ! under which the code may be brief.
+    !-----------------------------------------------------------------------------
+    !> @brief Module containing adjoint test for atl_bdy_lyr_alg
+    module atlt_bdy_lyr_alg_mod
+      use sci_assign_field_random_range_alg_mod, &
+                                               only : assign_field_random_range
+      use field_mod,                           only : field_type
+      use function_space_mod,                  only : function_space_type
+      use constants_mod,                       only : r_def, i_def, l_def
+      use field_indices_mod,                   only : igh_u, igh_t, &
+                                                      igh_d, igh_p
+      use log_mod,                             only : log_event,         &
+                                                      log_scratch_space, &
+                                                      LOG_LEVEL_INFO,    &
+                                                      LOG_LEVEL_DEBUG,   &
+                                                      LOG_LEVEL_ERROR
+      use mesh_mod,                            only : mesh_type
+      use derived_config_mod,                  only : bundle_size
+      use driver_modeldb_mod,                  only : modeldb_type
+      use finite_element_config_mod,           only : element_order_h, element_order_v
+      use fs_continuity_mod,                   only : W2, W3, Wtheta
+      use function_space_collection_mod,       only : function_space_collection
+      use adjoint_test_parameters_mod,         only : ls_u_range, ls_theta_range,   &
+                                                      ls_rho_range, ls_exner_range
+      use timing_mod,                          only : start_timing, stop_timing, tik, LPROF
+      implicit none
+      public
+      contains
+      !=============================================================================
+      !> @brief   Adjoint test for atl_bdy_lyr_alg.
+      !> @details Passes if adjoint is transpose of tangent linear.
+      !>          Determined by testing the equality of inner products <Mx, Mx> and <AMx, x>,
+      !>          where M is the tangent linear and A is the adjoint.
+      !> @param[in,out] modeldb       Structure containing the model state
+      !> @param[in]     mesh          The model mesh
+      subroutine atlt_bdy_lyr_alg( modeldb, mesh )
+        use tl_bdy_lyr_alg_mod,                      only : tl_bdy_lyr_alg
+        use atl_bdy_lyr_alg_mod,                     only : atl_bdy_lyr_alg
+        implicit none
+        ! Arguments
+        type(modeldb_type), target,     intent(inout) :: modeldb
+        type(mesh_type), pointer,       intent(in)    :: mesh
+        ! Arguments for tl and atl calls
+        ! Form of state is [u,theta,rho,exner]
+        type(field_type)    :: state(bundle_size)
+        type(field_type)    :: u_bl_inc
+        type(field_type)    :: ls_state(bundle_size)
+        ! Copies of input fields used in inner products
+        type(field_type)    :: state_input(bundle_size)
+        type(field_type)    :: u_bl_inc_input
+        ! Pointers for initialising fields
+        type(function_space_type), pointer :: vector_space_wtheta_ptr
+        type(function_space_type), pointer :: vector_space_w2_ptr
+        type(function_space_type), pointer :: vector_space_w3_ptr
+        ! Inner products
+        real(kind=r_def)    :: inner1
+        real(kind=r_def)    :: inner2
+        real(kind=r_def)    :: ip1(2),ip2(2)
+        real(kind=r_def)    :: sf(2)
+        ! Test parameters and variables
+        real(kind=r_def), parameter :: overall_tolerance = 1500.0_r_def
+        real(kind=r_def) :: machine_tol
+        real(kind=r_def) :: relative_diff
+        real(kind=r_def), parameter :: eps = 1e-30_r_def
+        ! Misc
+        real(kind=r_def)    :: dt
+        integer(kind=tik) :: id
+        if ( LPROF ) call start_timing( id, 'atlt_bdy_lyr_alg' )
+        ! Determining time constants
+        dt  = real(modeldb%clock%get_seconds_per_step(), r_def)
+        if ( dt <= 0.0_r_def ) then
+          write( log_scratch_space, * ) "dt cannot be leq 0, found dt = ", dt
+          call log_event( log_scratch_space, log_level_error )
+        end if
+        ! Initialising fields
+        vector_space_wtheta_ptr => function_space_collection % get_fs( mesh, element_order_h, element_order_v, Wtheta )
+        vector_space_w2_ptr => function_space_collection % get_fs( mesh, element_order_h, element_order_v, W2 )
+        vector_space_w3_ptr => function_space_collection % get_fs( mesh, element_order_h, element_order_v, W3 )
+        call state(igh_u) % initialise( vector_space = vector_space_w2_ptr )
+        call state(igh_t) % initialise( vector_space = vector_space_wtheta_ptr )
+        call state(igh_d) % initialise( vector_space = vector_space_w3_ptr )
+        call state(igh_p) % initialise( vector_space = vector_space_w3_ptr )
+        call u_bl_inc     % initialise( vector_space = vector_space_w2_ptr )
+        call state(igh_u) % copy_field_properties( state_input(igh_u) )
+        call state(igh_t) % copy_field_properties( state_input(igh_t) )
+        call state(igh_d) % copy_field_properties( state_input(igh_d) )
+        call state(igh_p) % copy_field_properties( state_input(igh_p) )
+        call state(igh_u) % copy_field_properties( u_bl_inc_input )
+        call state(igh_u) % copy_field_properties( ls_state(igh_u) )
+        call state(igh_t) % copy_field_properties( ls_state(igh_t) )
+        call state(igh_d) % copy_field_properties( ls_state(igh_d) )
+        call state(igh_p) % copy_field_properties( ls_state(igh_p) )
+        ! Initialise values and call the tangent-linear alg.
+        call invoke( setval_random( u_bl_inc ),                     &
+                     setval_x(u_bl_inc_input, u_bl_inc ),           &
+                     setval_random( state(igh_u) ),                 &
+                     setval_x( state_input(igh_u), state(igh_u) ),  &
+                     setval_random( state(igh_t) ),                 &
+                     setval_x( state_input(igh_t), state(igh_t) ),  &
+                     setval_random( state(igh_d) ),                 &
+                     setval_x( state_input(igh_d), state(igh_d) ),  &
+                     setval_random( state(igh_p) ),                 &
+                     setval_x( state_input(igh_p), state(igh_p) ) )
+        ! LS init
+        call assign_field_random_range( ls_state(igh_u), ls_u_range(1), ls_u_range(2) )
+        call assign_field_random_range( ls_state(igh_t), ls_theta_range(1), ls_theta_range(2) )
+        call invoke( setval_random( ls_state(igh_d) ), setval_random( ls_state(igh_p) ) )
+        ! Tangent linear
+        call tl_bdy_lyr_alg(modeldb, u_bl_inc, state(igh_u), ls_state, dt )
+        ! < Mx, Mx >
+        call invoke( x_innerproduct_x( ip1(1), state(igh_u) ), &
+                     x_innerproduct_x( ip1(2), u_bl_inc ) )
+        sf(1) = 1.0_r_def / (ip1(1) + eps)
+        sf(2) = 1.0_r_def / (ip1(2) + eps)
+        inner1 = 0.0_r_def
+        inner1 = inner1 + ip1(1) * sf(1)
+        inner1 = inner1 + ip1(2) * sf(2)
+        ! Scaling fields
+        call invoke( inc_a_times_X( sf(1), state(igh_u) ), &
+                     inc_a_times_X( sf(2), u_bl_inc ) )
+        ! Adjoint alg call and inner products
+        ! Adjoint
+        call atl_bdy_lyr_alg(modeldb, u_bl_inc, state(igh_u), ls_state, dt )
+        ! < AMx, x >
+        call invoke( x_innerproduct_y( ip2(1),               &
+                                       state(igh_u),         &
+                                       state_input(igh_u) ), &
+                     x_innerproduct_y( ip2(2),               &
+                                       u_bl_inc,             &
+                                       u_bl_inc_input ) )
+        inner2 = 0.0_r_def
+        inner2 = inner2 + ip2(1)
+        inner2 = inner2 + ip2(2)
+        write( log_scratch_space, * ) 'atlt_bdy_lyr_alg: ad test for tl_bdy_lyr_alg: ip1=', ip1(1) * sf(1), ip1(2) * sf(2)
+        call log_event( log_scratch_space, LOG_LEVEL_DEBUG )
+        write( log_scratch_space, * ) 'atlt_bdy_lyr_alg: ad test for tl_bdy_lyr_alg: ip2=', ip2
+        call log_event( log_scratch_space, LOG_LEVEL_DEBUG )
+        write( log_scratch_space, * ) 'atlt_bdy_lyr_alg: ad test for tl_bdy_lyr_alg: s(ip1),s(ip2)=', inner1, inner2
+        call log_event( log_scratch_space, LOG_LEVEL_DEBUG )
+        ! Test the inner-product values for equality, allowing for the precision of the active variables
+        machine_tol = spacing( max( abs( inner1 ), abs( inner2 ) ) )
+        relative_diff = abs( inner1 - inner2 ) / machine_tol
+        if ( relative_diff < overall_tolerance ) then
+          write( log_scratch_space, * ) "PASSED tl_bdy_lyr_alg:", inner1, inner2, relative_diff
+          call log_event( log_scratch_space, LOG_LEVEL_INFO )
+        else
+          write( log_scratch_space, * ) "FAILED tl_bdy_lyr_alg:", inner1, inner2, relative_diff
+          call log_event( log_scratch_space, LOG_LEVEL_ERROR )
+        end if
+        if ( LPROF ) call stop_timing( id, 'atlt_bdy_lyr_alg' )
+      end subroutine atlt_bdy_lyr_alg
+    end module atlt_bdy_lyr_alg_mod

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Linear and adjoint boundary layer physics #182

Uh oh!

Diff view

Diff view

There are no files selected for viewing

Uh oh!

Uh oh!

Linear and adjoint boundary layer physics #182

Are you sure you want to change the base?

Uh oh!

Linear and adjoint boundary layer physics #182

Uh oh!

Uh oh!

Diff view

Diff view

There are no files selected for viewing

Uh oh!

Uh oh!