diff --git a/applications/adjoint_tests/source/algorithm/algebra/adjt_dg_inc_matrix_vector_alg_mod.x90 b/applications/adjoint_tests/source/algorithm/algebra/adjt_dg_inc_matrix_vector_alg_mod.x90 new file mode 100644 index 000000000..4bf0b46a8 --- /dev/null +++ b/applications/adjoint_tests/source/algorithm/algebra/adjt_dg_inc_matrix_vector_alg_mod.x90 @@ -0,0 +1,145 @@ +!----------------------------------------------------------------------------- +! (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 used. +!----------------------------------------------------------------------------- +!> @brief Module containing adjoint test for adj_dg_inc_matrix_vector_kernel_type +module adjt_dg_inc_matrix_vector_alg_mod + + use field_mod, only : field_type + use function_space_mod, only : function_space_type + use mesh_mod, only : mesh_type + use function_space_collection_mod, only : function_space_collection + use finite_element_config_mod, only : element_order_h, element_order_v + use fs_continuity_mod, only : W0, W3 + use operator_mod, only : operator_type + use constants_mod, only : i_def, r_def, EPS + use setop_random_alg_mod, only : setop_random_alg + use log_mod, only : log_event, & + log_scratch_space, & + LOG_LEVEL_INFO, & + LOG_LEVEL_DEBUG, & + LOG_LEVEL_ERROR + + implicit none + + private + + public :: adjt_dg_inc_matrix_vector_alg + + contains + + !============================================================================= + !> @brief Adjoint test for adj_dg_inc_matrix_vector_kernel_type. + !> @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] mesh The model mesh + subroutine adjt_dg_inc_matrix_vector_alg(mesh) + + use dg_inc_matrix_vector_kernel_mod, only : dg_inc_matrix_vector_kernel_type + use adj_dg_inc_matrix_vector_kernel_mod, only : adj_dg_inc_matrix_vector_kernel_type + + implicit none + + ! Arguments + type(mesh_type), pointer, intent(in) :: mesh + + ! Arguments for tl and adj calls + type(field_type) :: lhs + type(field_type) :: x + type(operator_type) :: matrix + + ! Copies of input fields used in inner products + type(field_type) :: lhs_input + type(field_type) :: x_input + + ! Variables for initialising fields + type(function_space_type), pointer :: vector_space_w0_ptr + type(function_space_type), pointer :: vector_space_w3_ptr + + ! Inner products + real(kind=r_def) :: lhs_inner_prod + real(kind=r_def) :: x_inner_prod + real(kind=r_def) :: lhs_sf + real(kind=r_def) :: x_sf + real(kind=r_def) :: inner1 + real(kind=r_def) :: lhs_lhs_input_inner_prod + real(kind=r_def) :: x_x_input_inner_prod + real(kind=r_def) :: inner2 + + ! 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 + + vector_space_w0_ptr => function_space_collection%get_fs(mesh,element_order_h,element_order_v,w0) + vector_space_w3_ptr => function_space_collection%get_fs(mesh,element_order_h,element_order_v,w3) + + call lhs%initialise(vector_space = vector_space_w3_ptr, name = 'lhs') + call x%initialise(vector_space = vector_space_w0_ptr, name = 'x') + call matrix%initialise(vector_space_w3_ptr, vector_space_w0_ptr) + call lhs_input%initialise(vector_space = vector_space_w3_ptr, name = 'lhs_input') + call x_input%initialise(vector_space = vector_space_w0_ptr, name = 'x_input') + + lhs_inner_prod = 0.0_r_def + x_inner_prod = 0.0_r_def + + ! Initialise arguments and call the tangent-linear kernel. + call invoke( setval_random(lhs), & + setval_x(lhs_input, lhs), & + setval_random(x), & + setval_x(x_input, x) ) + call setop_random_alg( matrix ) + + ! Tangent linear and + call invoke( dg_inc_matrix_vector_kernel_type(lhs, x, matrix), & + x_innerproduct_x(lhs_inner_prod, lhs), & + x_innerproduct_x(x_inner_prod, x) ) + + ! Determining scale factors + lhs_sf = 1.0_r_def / (lhs_inner_prod + eps) + x_sf = 1.0_r_def / (x_inner_prod + eps) + + inner1 = 0.0_r_def + inner1 = inner1 + lhs_inner_prod * lhs_sf + inner1 = inner1 + x_inner_prod * x_sf + + write( log_scratch_space, * ) "adjt_dg_matrix_vector_alg inner products:" + call log_event( log_scratch_space, log_level_debug ) + write( log_scratch_space, * ) "lhs inner product = ", lhs_inner_prod + call log_event( log_scratch_space, log_level_debug ) + write( log_scratch_space, * ) "x inner product = ", x_inner_prod + call log_event( log_scratch_space, log_level_debug ) + + ! Scaling fields + call invoke( inc_a_times_X( lhs_sf, lhs ), & + inc_a_times_X( x_sf, x ) ) + + ! Adjoint and + lhs_lhs_input_inner_prod = 0.0_r_def + x_x_input_inner_prod = 0.0_r_def + + call invoke( adj_dg_inc_matrix_vector_kernel_type(lhs, x, matrix), & + x_innerproduct_y(lhs_lhs_input_inner_prod, lhs, lhs_input), & + x_innerproduct_y(x_x_input_inner_prod, x, x_input) ) + + inner2 = 0.0_r_def + inner2 = inner2 + lhs_lhs_input_inner_prod + inner2 = inner2 + x_x_input_inner_prod + + ! 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 dg_inc_matrix_vector_kernel_type:", inner1, inner2, relative_diff + call log_event( log_scratch_space, LOG_LEVEL_INFO ) + else + write(log_scratch_space, *) "FAILED dg_inc_matrix_vector_kernel_type:", inner1, inner2, relative_diff + call log_event( log_scratch_space, LOG_LEVEL_ERROR ) + end if + + end subroutine adjt_dg_inc_matrix_vector_alg + +end module adjt_dg_inc_matrix_vector_alg_mod diff --git a/applications/adjoint_tests/source/algorithm/algebra/adjt_dg_matrix_vector_alg_mod.x90 b/applications/adjoint_tests/source/algorithm/algebra/adjt_dg_matrix_vector_alg_mod.x90 new file mode 100644 index 000000000..7b4d662f0 --- /dev/null +++ b/applications/adjoint_tests/source/algorithm/algebra/adjt_dg_matrix_vector_alg_mod.x90 @@ -0,0 +1,145 @@ +!----------------------------------------------------------------------------- +! (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 used. +!----------------------------------------------------------------------------- +!> @brief Module containing adjoint test for adj_dg_matrix_vector_kernel_type +module adjt_dg_matrix_vector_alg_mod + + use field_mod, only : field_type + use function_space_mod, only : function_space_type + use mesh_mod, only : mesh_type + use function_space_collection_mod, only : function_space_collection + use finite_element_config_mod, only : element_order_h, element_order_v + use fs_continuity_mod, only : W0, W3 + use operator_mod, only : operator_type + use constants_mod, only : i_def, r_def, EPS + use setop_random_alg_mod, only : setop_random_alg + use log_mod, only : log_event, & + log_scratch_space, & + LOG_LEVEL_INFO, & + LOG_LEVEL_DEBUG, & + LOG_LEVEL_ERROR + + implicit none + + private + + public :: adjt_dg_matrix_vector_alg + + contains + + !============================================================================= + !> @brief Adjoint test for adj_dg_matrix_vector_kernel_type. + !> @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] mesh The model mesh + subroutine adjt_dg_matrix_vector_alg(mesh) + + use dg_matrix_vector_kernel_mod, only : dg_matrix_vector_kernel_type + use adj_dg_matrix_vector_kernel_mod, only : adj_dg_matrix_vector_kernel_type + + implicit none + + ! Arguments + type(mesh_type), pointer, intent(in) :: mesh + + ! Arguments for tl and adj calls + type(field_type) :: lhs + type(field_type) :: x + type(operator_type) :: matrix + + ! Copies of input fields used in inner products + type(field_type) :: lhs_input + type(field_type) :: x_input + + ! Variables for initialising fields + type(function_space_type), pointer :: vector_space_w0_ptr + type(function_space_type), pointer :: vector_space_w3_ptr + + ! Inner products + real(kind=r_def) :: lhs_inner_prod + real(kind=r_def) :: x_inner_prod + real(kind=r_def) :: lhs_sf + real(kind=r_def) :: x_sf + real(kind=r_def) :: inner1 + real(kind=r_def) :: lhs_lhs_input_inner_prod + real(kind=r_def) :: x_x_input_inner_prod + real(kind=r_def) :: inner2 + + ! 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 + + vector_space_w0_ptr => function_space_collection%get_fs(mesh,element_order_h,element_order_v,w0) + vector_space_w3_ptr => function_space_collection%get_fs(mesh,element_order_h,element_order_v,w3) + + call lhs%initialise(vector_space = vector_space_w3_ptr, name = 'lhs') + call x%initialise(vector_space = vector_space_w0_ptr, name = 'x') + call matrix%initialise(vector_space_w3_ptr, vector_space_w0_ptr) + call lhs_input%initialise(vector_space = vector_space_w3_ptr, name = 'lhs_input') + call x_input%initialise(vector_space = vector_space_w0_ptr, name = 'x_input') + + lhs_inner_prod = 0.0_r_def + x_inner_prod = 0.0_r_def + + ! Initialise arguments and call the tangent-linear kernel. + call invoke( setval_random(lhs), & + setval_x(lhs_input, lhs), & + setval_random(x), & + setval_x(x_input, x) ) + call setop_random_alg( matrix ) + + ! Tangent linear and + call invoke( dg_matrix_vector_kernel_type(lhs, x, matrix), & + x_innerproduct_x(lhs_inner_prod, lhs), & + x_innerproduct_x(x_inner_prod, x) ) + + ! Determining scale factors + lhs_sf = 1.0_r_def / (lhs_inner_prod + eps) + x_sf = 1.0_r_def / (x_inner_prod + eps) + + inner1 = 0.0_r_def + inner1 = inner1 + lhs_inner_prod * lhs_sf + inner1 = inner1 + x_inner_prod * x_sf + + write( log_scratch_space, * ) "adjt_dg_matrix_vector_alg inner products:" + call log_event( log_scratch_space, log_level_debug ) + write( log_scratch_space, * ) "lhs inner product = ", lhs_inner_prod + call log_event( log_scratch_space, log_level_debug ) + write( log_scratch_space, * ) "x inner product = ", x_inner_prod + call log_event( log_scratch_space, log_level_debug ) + + ! Scaling fields + call invoke( inc_a_times_X( lhs_sf, lhs ), & + inc_a_times_X( x_sf, x ) ) + + ! Adjoint and + lhs_lhs_input_inner_prod = 0.0_r_def + x_x_input_inner_prod = 0.0_r_def + + call invoke( adj_dg_matrix_vector_kernel_type(lhs, x, matrix), & + x_innerproduct_y(lhs_lhs_input_inner_prod, lhs, lhs_input), & + x_innerproduct_y(x_x_input_inner_prod, x, x_input) ) + + inner2 = 0.0_r_def + inner2 = inner2 + lhs_lhs_input_inner_prod + inner2 = inner2 + x_x_input_inner_prod + + ! 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 dg_matrix_vector_kernel_type:", inner1, inner2, relative_diff + call log_event( log_scratch_space, LOG_LEVEL_INFO ) + else + write(log_scratch_space, *) "FAILED dg_matrix_vector_kernel_type:", inner1, inner2, relative_diff + call log_event( log_scratch_space, LOG_LEVEL_ERROR ) + end if + + end subroutine adjt_dg_matrix_vector_alg + +end module adjt_dg_matrix_vector_alg_mod diff --git a/applications/adjoint_tests/source/algorithm/algebra/adjt_matrix_vector_alg_mod.x90 b/applications/adjoint_tests/source/algorithm/algebra/adjt_matrix_vector_alg_mod.x90 new file mode 100644 index 000000000..e19cc6748 --- /dev/null +++ b/applications/adjoint_tests/source/algorithm/algebra/adjt_matrix_vector_alg_mod.x90 @@ -0,0 +1,142 @@ +!----------------------------------------------------------------------------- +! (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 used. +!----------------------------------------------------------------------------- +!> @brief Module containing adjoint test for adj_matrix_vector_kernel_type +module adjt_matrix_vector_alg_mod + + use field_mod, only : field_type + use function_space_mod, only : function_space_type + use mesh_mod, only : mesh_type + use function_space_collection_mod, only : function_space_collection + use finite_element_config_mod, only : element_order_h, element_order_v + use fs_continuity_mod, only : W0 + use operator_mod, only : operator_type + use constants_mod, only : i_def, r_def, EPS + use setop_random_alg_mod, only : setop_random_alg + use log_mod, only : log_event, & + log_scratch_space, & + LOG_LEVEL_INFO, & + LOG_LEVEL_DEBUG, & + LOG_LEVEL_ERROR + implicit none + + private + + public :: adjt_matrix_vector_alg + + contains + + !============================================================================= + !> @brief Adjoint test for adj_matrix_vector_kernel_type. + !> @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] mesh The model mesh + subroutine adjt_matrix_vector_alg(mesh) + + use matrix_vector_kernel_mod, only : matrix_vector_kernel_type + use adj_matrix_vector_kernel_mod, only : adj_matrix_vector_kernel_type + + implicit none + + ! Arguments + type(mesh_type), pointer, intent(in) :: mesh + + ! Arguments for tl and adj calls + type(field_type) :: lhs + type(field_type) :: x + type(operator_type) :: matrix + + ! Copies of input fields used in inner products + type(field_type) :: lhs_input + type(field_type) :: x_input + + ! Variables for initialising fields + type(function_space_type), pointer :: vector_space_w0_ptr + + ! Inner products + real(kind=r_def) :: lhs_inner_prod + real(kind=r_def) :: x_inner_prod + real(kind=r_def) :: lhs_sf + real(kind=r_def) :: x_sf + real(kind=r_def) :: inner1 + real(kind=r_def) :: lhs_lhs_input_inner_prod + real(kind=r_def) :: x_x_input_inner_prod + real(kind=r_def) :: inner2 + + ! 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 + + vector_space_w0_ptr => function_space_collection%get_fs(mesh, element_order_h, & + element_order_v, W0) + call lhs%initialise(vector_space = vector_space_w0_ptr, name = 'lhs') + call x%initialise(vector_space = vector_space_w0_ptr, name = 'x') + call matrix%initialise(vector_space_w0_ptr, vector_space_w0_ptr) + call lhs_input%initialise(vector_space = vector_space_w0_ptr, name = 'lhs_input') + call x_input%initialise(vector_space = vector_space_w0_ptr, name = 'x_input') + + lhs_inner_prod = 0.0_r_def + x_inner_prod = 0.0_r_def + + ! Initialise arguments and call the tangent-linear kernel. + call invoke( setval_random(lhs), & + setval_x(lhs_input, lhs), & + setval_random(x), & + setval_x(x_input, x) ) + call setop_random_alg( matrix ) + + ! Tangent linear and + call invoke( matrix_vector_kernel_type(lhs, x, matrix), & + x_innerproduct_x(lhs_inner_prod, lhs), & + x_innerproduct_x(x_inner_prod, x) ) + + ! Determining scale factors + lhs_sf = 1.0_r_def / (lhs_inner_prod + EPS) + x_sf = 1.0_r_def / (x_inner_prod + EPS) + + inner1 = 0.0_r_def + inner1 = inner1 + lhs_inner_prod * lhs_sf + inner1 = inner1 + x_inner_prod * x_sf + + write( log_scratch_space, * ) "adjt_matrix_vector_alg inner products:" + call log_event( log_scratch_space, log_level_debug ) + write( log_scratch_space, * ) "lhs inner product = ", lhs_inner_prod + call log_event( log_scratch_space, log_level_debug ) + write( log_scratch_space, * ) "x inner product = ", x_inner_prod + call log_event( log_scratch_space, log_level_debug ) + + ! Scaling fields + call invoke( inc_a_times_X( lhs_sf, lhs ), & + inc_a_times_X( x_sf, x ) ) + + ! Adjoint and + lhs_lhs_input_inner_prod = 0.0_r_def + x_x_input_inner_prod = 0.0_r_def + + call invoke( adj_matrix_vector_kernel_type(lhs, x, matrix), & + x_innerproduct_y(lhs_lhs_input_inner_prod, lhs, lhs_input), & + x_innerproduct_y(x_x_input_inner_prod, x, x_input)) + + inner2 = 0.0_r_def + inner2 = inner2 + lhs_lhs_input_inner_prod + inner2 = inner2 + x_x_input_inner_prod + + ! 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 matrix_vector_kernel_type:", inner1, inner2, relative_diff + call log_event( log_scratch_space, LOG_LEVEL_INFO ) + else + write(log_scratch_space, *) "FAILED matrix_vector_kernel_type:", inner1, inner2, relative_diff + call log_event( log_scratch_space, LOG_LEVEL_ERROR ) + end if + + end subroutine adjt_matrix_vector_alg + +end module adjt_matrix_vector_alg_mod diff --git a/applications/adjoint_tests/source/algorithm/algebra/adjt_transpose_matrix_vector_alg_mod.x90 b/applications/adjoint_tests/source/algorithm/algebra/adjt_transpose_matrix_vector_alg_mod.x90 new file mode 100644 index 000000000..f6b654236 --- /dev/null +++ b/applications/adjoint_tests/source/algorithm/algebra/adjt_transpose_matrix_vector_alg_mod.x90 @@ -0,0 +1,141 @@ +!----------------------------------------------------------------------------- +! (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 used. +!----------------------------------------------------------------------------- +!> @brief Module containing adjoint test for adj_transpose_matrix_vector_kernel_type +module adjt_transpose_matrix_vector_alg_mod + use field_mod, only : field_type + use function_space_mod, only : function_space_type + use mesh_mod, only : mesh_type + use function_space_collection_mod, only : function_space_collection + use finite_element_config_mod, only : element_order_h, element_order_v + use fs_continuity_mod, only : W0 + use operator_mod, only : operator_type + use constants_mod, only : i_def, r_def, EPS + use setop_random_alg_mod, only : setop_random_alg + use log_mod, only : log_event, & + log_scratch_space, & + LOG_LEVEL_INFO, & + LOG_LEVEL_DEBUG, & + LOG_LEVEL_ERROR + implicit none + + private + + public :: adjt_transpose_matrix_vector_alg + + contains + + !============================================================================= + !> @brief Adjoint test for adj_transpose_matrix_vector_kernel_type. + !> @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] mesh The model mesh + subroutine adjt_transpose_matrix_vector_alg(mesh) + + use transpose_matrix_vector_kernel_mod, only : transpose_matrix_vector_kernel_type + use adj_transpose_matrix_vector_kernel_mod, only : adj_transpose_matrix_vector_kernel_type + +implicit none + + ! Arguments + type(mesh_type), pointer, intent(in) :: mesh + + ! Arguments for tl and adj calls + type(field_type) :: lhs + type(field_type) :: x + type(operator_type) :: matrix + + ! Copies of input fields used in inner products + type(field_type) :: lhs_input + type(field_type) :: x_input + + ! Variables for initialising fields + type(function_space_type), pointer :: vector_space_w0_ptr + + ! Inner products + real(kind=r_def) :: lhs_inner_prod + real(kind=r_def) :: x_inner_prod + real(kind=r_def) :: lhs_sf + real(kind=r_def) :: x_sf + real(kind=r_def) :: inner1 + real(kind=r_def) :: lhs_lhs_input_inner_prod + real(kind=r_def) :: x_x_input_inner_prod + real(kind=r_def) :: inner2 + + ! 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 + + vector_space_w0_ptr => function_space_collection%get_fs(mesh, element_order_h, & + element_order_v, W0) + call lhs%initialise(vector_space = vector_space_w0_ptr, name = 'lhs') + call x%initialise(vector_space = vector_space_w0_ptr, name = 'x') + call matrix%initialise(vector_space_w0_ptr, vector_space_w0_ptr) + call lhs_input%initialise(vector_space = vector_space_w0_ptr, name = 'lhs_input') + call x_input%initialise(vector_space = vector_space_w0_ptr, name = 'x_input') + + lhs_inner_prod = 0.0_r_def + x_inner_prod = 0.0_r_def + + ! Initialise arguments and call the tangent-linear kernel. + call invoke( setval_random(lhs), & + setval_x(lhs_input, lhs), & + setval_random(x), & + setval_x(x_input, x) ) + call setop_random_alg( matrix ) + + ! Tangent linear and + call invoke( transpose_matrix_vector_kernel_type(lhs, x, matrix), & + x_innerproduct_x(lhs_inner_prod, lhs), & + x_innerproduct_x(x_inner_prod, x) ) + + ! Determining scale factors + lhs_sf = 1.0_r_def / (lhs_inner_prod + EPS) + x_sf = 1.0_r_def / (x_inner_prod + EPS) + + inner1 = 0.0_r_def + inner1 = inner1 + lhs_inner_prod * lhs_sf + inner1 = inner1 + x_inner_prod * x_sf + + write( log_scratch_space, * ) "adjt_matrix_vector_alg inner products:" + call log_event( log_scratch_space, log_level_debug ) + write( log_scratch_space, * ) "lhs inner product = ", lhs_inner_prod + call log_event( log_scratch_space, log_level_debug ) + write( log_scratch_space, * ) "x inner product = ", x_inner_prod + call log_event( log_scratch_space, log_level_debug ) + + ! Scaling fields + call invoke( inc_a_times_X( lhs_sf, lhs ), & + inc_a_times_X( x_sf, x ) ) + + ! Adjoint and + lhs_lhs_input_inner_prod = 0.0_r_def + x_x_input_inner_prod = 0.0_r_def + + call invoke( adj_transpose_matrix_vector_kernel_type(lhs, x, matrix), & + x_innerproduct_y(lhs_lhs_input_inner_prod, lhs, lhs_input), & + x_innerproduct_y(x_x_input_inner_prod, x, x_input)) + + inner2 = 0.0_r_def + inner2 = inner2 + lhs_lhs_input_inner_prod + inner2 = inner2 + x_x_input_inner_prod + + ! 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 transpose_matrix_vector_kernel_type:", inner1, inner2, relative_diff + call log_event( log_scratch_space, LOG_LEVEL_INFO ) + else + write(log_scratch_space, *) "FAILED transpose_matrix_vector_kernel_type:", inner1, inner2, relative_diff + call log_event( log_scratch_space, LOG_LEVEL_ERROR ) + end if + + end subroutine adjt_transpose_matrix_vector_alg + +end module adjt_transpose_matrix_vector_alg_mod diff --git a/applications/adjoint_tests/source/driver/adjoint_test_driver_mod.f90 b/applications/adjoint_tests/source/driver/adjoint_test_driver_mod.f90 index 9f426f45d..737d81dbb 100644 --- a/applications/adjoint_tests/source/driver/adjoint_test_driver_mod.f90 +++ b/applications/adjoint_tests/source/driver/adjoint_test_driver_mod.f90 @@ -43,6 +43,12 @@ subroutine run( modeldb ) ! ./ use adjt_sci_psykal_builtin_alg_mod, only : adjt_convert_cart2sphere_vector_alg + ! ./algebra + use adjt_matrix_vector_alg_mod, only : adjt_matrix_vector_alg + use adjt_dg_matrix_vector_alg_mod, only : adjt_dg_matrix_vector_alg + use adjt_dg_inc_matrix_vector_alg_mod, only : adjt_dg_inc_matrix_vector_alg + use adjt_transpose_matrix_vector_alg_mod, only : adjt_transpose_matrix_vector_alg + ! ./inter_function_space use adjt_sci_convert_hdiv_field_alg_mod, only : adjt_sci_convert_hdiv_field_alg @@ -159,6 +165,12 @@ subroutine run( modeldb ) ! ./inter_function_space call adjt_sci_convert_hdiv_field_alg( mesh, chi, panel_id ) + ! ./algebra + call adjt_matrix_vector_alg( mesh ) + call adjt_dg_matrix_vector_alg( mesh ) + call adjt_dg_inc_matrix_vector_alg( mesh ) + call adjt_transpose_matrix_vector_alg( mesh ) + call log_event( "TESTING misc adjoints", LOG_LEVEL_INFO ) ! ./ call adjt_convert_cart2sphere_vector_alg( mesh ) diff --git a/science/adjoint/build/psyad_files_list_core.txt b/science/adjoint/build/psyad_files_list_core.txt index 2f1e378ec..410e9fcdf 100644 --- a/science/adjoint/build/psyad_files_list_core.txt +++ b/science/adjoint/build/psyad_files_list_core.txt @@ -7,7 +7,3 @@ components/science/source/kernel/inter_function_space/sample_field_kernel_mod.F9 components/science/source/kernel/inter_function_space/sample_flux_kernel_mod.F90 components/science/source/kernel/inter_function_space/split_w2_field_kernel_mod.F90 components/science/source/kernel/fem/strong_curl_kernel_mod.F90 -components/science/source/kernel/algebra/matrix_vector_kernel_mod.F90 -components/science/source/kernel/algebra/dg_inc_matrix_vector_kernel_mod.F90 -components/science/source/kernel/algebra/transpose_matrix_vector_kernel_mod.F90 -components/science/source/kernel/algebra/dg_matrix_vector_kernel_mod.F90 diff --git a/science/adjoint/build/psyad_vars.mk b/science/adjoint/build/psyad_vars.mk index 15900420a..c6b526ca9 100644 --- a/science/adjoint/build/psyad_vars.mk +++ b/science/adjoint/build/psyad_vars.mk @@ -52,10 +52,6 @@ all: export ACTIVE_sample_field_kernel_mod := field_1 field_2 f_at_n all: export ACTIVE_sample_flux_kernel_mod := flux u all: export ACTIVE_split_w2_field_kernel_mod := uvw w uv all: export ACTIVE_strong_curl_kernel_mod := xi res_dot_product curl_u u -all: export ACTIVE_matrix_vector_kernel_mod := lhs x -all: export ACTIVE_dg_inc_matrix_vector_kernel_mod := lhs x -all: export ACTIVE_transpose_matrix_vector_kernel_mod := lhs lhs_e x x_e -all: export ACTIVE_dg_matrix_vector_kernel_mod := lhs x all: export ACTIVE_sci_average_w2b_to_w2_kernel_mod := field_w2 field_w2_broken all: export ACTIVE_sci_extract_w_kernel_mod := velocity_w2v u_physics all: export ACTIVE_sci_combine_multidata_field_kernel_mod := field1_in field2_in field_out diff --git a/science/adjoint/patches/algorithm/adjt_dg_inc_matrix_vector_alg_mod.patch b/science/adjoint/patches/algorithm/adjt_dg_inc_matrix_vector_alg_mod.patch deleted file mode 100644 index 0b76557a5..000000000 --- a/science/adjoint/patches/algorithm/adjt_dg_inc_matrix_vector_alg_mod.patch +++ /dev/null @@ -1,65 +0,0 @@ -@@ -14,8 +14,8 @@ - use fs_continuity_mod, only : w0, w3 - use operator_mod, only : operator_type - use constants_mod, only : i_def, r_def -- use setop_random_kernel_mod, only : setop_random_kernel_type -- use log_mod, only : log_event, log_level_error, log_level_info, log_scratch_space -+ use setop_random_alg_mod, only : setop_random_alg -+ use log_mod, only : log_event, log_level_error, log_level_info, log_scratch_space, log_level_debug - real(kind=r_def), parameter :: overall_tolerance = 1500.0_r_def - type(mesh_type), pointer, intent(in) :: mesh - type(field_type), dimension(3), intent(in), optional :: chi -@@ -29,12 +29,15 @@ - type(field_type) :: x_input - real(kind=r_def) :: lhs_inner_prod - real(kind=r_def) :: x_inner_prod -+ real(kind=r_def) :: lhs_sf -+ real(kind=r_def) :: x_sf - real(kind=r_def) :: inner1 - real(kind=r_def) :: lhs_lhs_input_inner_prod - real(kind=r_def) :: x_x_input_inner_prod - real(kind=r_def) :: inner2 - real(kind=r_def) :: MachineTol - real(kind=r_def) :: relative_diff -+ real(kind=r_def), parameter :: eps = 1.0e-30_r_def - - vector_space_w0_ptr => function_space_collection%get_fs(mesh,element_order_h,element_order_v,w0) - vector_space_w3_ptr => function_space_collection%get_fs(mesh,element_order_h,element_order_v,w3) -@@ -45,14 +48,24 @@ - call x_input%initialise(vector_space=vector_space_w0_ptr, name='x_input') - lhs_inner_prod = 0.0_r_def - x_inner_prod = 0.0_r_def -- - ! Initialise arguments and call the tangent-linear kernel. -- call invoke(setval_random(lhs), setval_x(lhs_input, lhs), setval_random(x), setval_x(x_input, x), & --&setop_random_kernel_type(matrix), dg_inc_matrix_vector_kernel_type(lhs, x, matrix), x_innerproduct_x(lhs_inner_prod, lhs), & -+ call invoke(setval_random(lhs), setval_x(lhs_input, lhs), setval_random(x), setval_x(x_input, x)) -+ call setop_random_alg(matrix) -+ call invoke(dg_inc_matrix_vector_kernel_type(lhs, x, matrix), x_innerproduct_x(lhs_inner_prod, lhs), & - &x_innerproduct_x(x_inner_prod, x)) -+ write( log_scratch_space, * ) "adjt_dg_inc_matrix_vector_alg inner products:" -+ call log_event( log_scratch_space, log_level_debug ) -+ write( log_scratch_space, * ) "lhs inner product = ", lhs_inner_prod -+ call log_event( log_scratch_space, log_level_debug ) -+ write( log_scratch_space, * ) "x inner product = ", x_inner_prod -+ call log_event( log_scratch_space, log_level_debug ) -+ lhs_sf = 1.0_r_def / (lhs_inner_prod + eps) -+ x_sf = 1.0_r_def / (x_inner_prod + eps) - inner1 = 0.0_r_def -- inner1 = inner1 + lhs_inner_prod -- inner1 = inner1 + x_inner_prod -+ inner1 = inner1 + lhs_inner_prod * lhs_sf -+ inner1 = inner1 + x_inner_prod * x_sf -+ call invoke( inc_a_times_X( lhs_sf, lhs ), & -+ inc_a_times_X( x_sf, x ) ) - lhs_lhs_input_inner_prod = 0.0_r_def - x_x_input_inner_prod = 0.0_r_def - call invoke(adj_dg_inc_matrix_vector_kernel_type(lhs, x, matrix), x_innerproduct_y(lhs_lhs_input_inner_prod, lhs, lhs_input), & -@@ -60,7 +73,6 @@ - inner2 = 0.0_r_def - inner2 = inner2 + lhs_lhs_input_inner_prod - inner2 = inner2 + x_x_input_inner_prod -- - ! Test the inner-product values for equality, allowing for the precision of the active variables - MachineTol = SPACING(MAX(ABS(inner1), ABS(inner2))) - relative_diff = ABS(inner1 - inner2) / MachineTol diff --git a/science/adjoint/patches/algorithm/adjt_dg_matrix_vector_alg_mod.patch b/science/adjoint/patches/algorithm/adjt_dg_matrix_vector_alg_mod.patch deleted file mode 100644 index 64ff46d28..000000000 --- a/science/adjoint/patches/algorithm/adjt_dg_matrix_vector_alg_mod.patch +++ /dev/null @@ -1,65 +0,0 @@ -@@ -14,8 +14,8 @@ - use fs_continuity_mod, only : w0, w3 - use operator_mod, only : operator_type - use constants_mod, only : i_def, r_def -- use setop_random_kernel_mod, only : setop_random_kernel_type -- use log_mod, only : log_event, log_level_error, log_level_info, log_scratch_space -+ use setop_random_alg_mod, only : setop_random_alg -+ use log_mod, only : log_event, log_level_error, log_level_info, log_scratch_space, log_level_debug - real(kind=r_def), parameter :: overall_tolerance = 1500.0_r_def - type(mesh_type), pointer, intent(in) :: mesh - type(field_type), dimension(3), intent(in), optional :: chi -@@ -29,12 +29,15 @@ - type(field_type) :: x_input - real(kind=r_def) :: lhs_inner_prod - real(kind=r_def) :: x_inner_prod -+ real(kind=r_def) :: lhs_sf -+ real(kind=r_def) :: x_sf - real(kind=r_def) :: inner1 - real(kind=r_def) :: lhs_lhs_input_inner_prod - real(kind=r_def) :: x_x_input_inner_prod - real(kind=r_def) :: inner2 - real(kind=r_def) :: MachineTol - real(kind=r_def) :: relative_diff -+ real(kind=r_def), parameter :: eps = 1.0e-30_r_def - - vector_space_w0_ptr => function_space_collection%get_fs(mesh,element_order_h,element_order_v,w0) - vector_space_w3_ptr => function_space_collection%get_fs(mesh,element_order_h,element_order_v,w3) -@@ -45,14 +48,24 @@ - call x_input%initialise(vector_space=vector_space_w0_ptr, name='x_input') - lhs_inner_prod = 0.0_r_def - x_inner_prod = 0.0_r_def -- - ! Initialise arguments and call the tangent-linear kernel. -- call invoke(setval_random(lhs), setval_x(lhs_input, lhs), setval_random(x), setval_x(x_input, x), & --&setop_random_kernel_type(matrix), dg_matrix_vector_kernel_type(lhs, x, matrix), x_innerproduct_x(lhs_inner_prod, lhs), & -+ call invoke(setval_random(lhs), setval_x(lhs_input, lhs), setval_random(x), setval_x(x_input, x)) -+ call setop_random_alg(matrix) -+ call invoke(dg_matrix_vector_kernel_type(lhs, x, matrix), x_innerproduct_x(lhs_inner_prod, lhs), & - &x_innerproduct_x(x_inner_prod, x)) -+ write( log_scratch_space, * ) "adjt_dg_matrix_vector_alg inner products:" -+ call log_event( log_scratch_space, log_level_debug ) -+ write( log_scratch_space, * ) "lhs inner product = ", lhs_inner_prod -+ call log_event( log_scratch_space, log_level_debug ) -+ write( log_scratch_space, * ) "x inner product = ", x_inner_prod -+ call log_event( log_scratch_space, log_level_debug ) -+ lhs_sf = 1.0_r_def / (lhs_inner_prod + eps) -+ x_sf = 1.0_r_def / (x_inner_prod + eps) - inner1 = 0.0_r_def -- inner1 = inner1 + lhs_inner_prod -- inner1 = inner1 + x_inner_prod -+ inner1 = inner1 + lhs_inner_prod * lhs_sf -+ inner1 = inner1 + x_inner_prod * x_sf -+ call invoke( inc_a_times_X( lhs_sf, lhs ), & -+ inc_a_times_X( x_sf, x ) ) - lhs_lhs_input_inner_prod = 0.0_r_def - x_x_input_inner_prod = 0.0_r_def - call invoke(adj_dg_matrix_vector_kernel_type(lhs, x, matrix), x_innerproduct_y(lhs_lhs_input_inner_prod, lhs, lhs_input), & -@@ -60,7 +73,6 @@ - inner2 = 0.0_r_def - inner2 = inner2 + lhs_lhs_input_inner_prod - inner2 = inner2 + x_x_input_inner_prod -- - ! Test the inner-product values for equality, allowing for the precision of the active variables - MachineTol = SPACING(MAX(ABS(inner1), ABS(inner2))) - relative_diff = ABS(inner1 - inner2) / MachineTol diff --git a/science/adjoint/patches/algorithm/adjt_matrix_vector_alg_mod.patch b/science/adjoint/patches/algorithm/adjt_matrix_vector_alg_mod.patch deleted file mode 100644 index ba1e86ecf..000000000 --- a/science/adjoint/patches/algorithm/adjt_matrix_vector_alg_mod.patch +++ /dev/null @@ -1,65 +0,0 @@ -@@ -14,8 +14,8 @@ - use fs_continuity_mod, only : w0 - use operator_mod, only : operator_type - use constants_mod, only : i_def, r_def -- use setop_random_kernel_mod, only : setop_random_kernel_type -- use log_mod, only : log_event, log_level_error, log_level_info, log_scratch_space -+ use setop_random_alg_mod, only : setop_random_alg -+ use log_mod, only : log_event, log_level_error, log_level_info, log_scratch_space, log_level_debug - real(kind=r_def), parameter :: overall_tolerance = 1500.0_r_def - type(mesh_type), pointer, intent(in) :: mesh - type(field_type), dimension(3), intent(in), optional :: chi -@@ -28,12 +28,15 @@ - type(field_type) :: x_input - real(kind=r_def) :: lhs_inner_prod - real(kind=r_def) :: x_inner_prod -+ real(kind=r_def) :: lhs_sf -+ real(kind=r_def) :: x_sf - real(kind=r_def) :: inner1 - real(kind=r_def) :: lhs_lhs_input_inner_prod - real(kind=r_def) :: x_x_input_inner_prod - real(kind=r_def) :: inner2 - real(kind=r_def) :: MachineTol - real(kind=r_def) :: relative_diff -+ real(kind=r_def), parameter :: eps = 1.0e-30_r_def - - vector_space_w0_ptr => function_space_collection%get_fs(mesh,element_order_h,element_order_v,w0) - call lhs%initialise(vector_space=vector_space_w0_ptr, name='lhs') -@@ -43,14 +46,24 @@ - call x_input%initialise(vector_space=vector_space_w0_ptr, name='x_input') - lhs_inner_prod = 0.0_r_def - x_inner_prod = 0.0_r_def -- - ! Initialise arguments and call the tangent-linear kernel. -- call invoke(setval_random(lhs), setval_x(lhs_input, lhs), setval_random(x), setval_x(x_input, x), & --&setop_random_kernel_type(matrix), matrix_vector_kernel_type(lhs, x, matrix), x_innerproduct_x(lhs_inner_prod, lhs), & -+ call invoke(setval_random(lhs), setval_x(lhs_input, lhs), setval_random(x), setval_x(x_input, x)) -+ call setop_random_alg(matrix) -+ call invoke(matrix_vector_kernel_type(lhs, x, matrix), x_innerproduct_x(lhs_inner_prod, lhs), & - &x_innerproduct_x(x_inner_prod, x)) -+ write( log_scratch_space, * ) "adjt_matrix_vector_alg inner products:" -+ call log_event( log_scratch_space, log_level_debug ) -+ write( log_scratch_space, * ) "lhs inner product = ", lhs_inner_prod -+ call log_event( log_scratch_space, log_level_debug ) -+ write( log_scratch_space, * ) "x inner product = ", x_inner_prod -+ call log_event( log_scratch_space, log_level_debug ) -+ lhs_sf = 1.0_r_def / (lhs_inner_prod + eps) -+ x_sf = 1.0_r_def / (x_inner_prod + eps) - inner1 = 0.0_r_def -- inner1 = inner1 + lhs_inner_prod -- inner1 = inner1 + x_inner_prod -+ inner1 = inner1 + lhs_inner_prod * lhs_sf -+ inner1 = inner1 + x_inner_prod * x_sf -+ call invoke( inc_a_times_X( lhs_sf, lhs ), & -+ inc_a_times_X( x_sf, x ) ) - lhs_lhs_input_inner_prod = 0.0_r_def - x_x_input_inner_prod = 0.0_r_def - call invoke(adj_matrix_vector_kernel_type(lhs, x, matrix), x_innerproduct_y(lhs_lhs_input_inner_prod, lhs, lhs_input), & -@@ -58,7 +71,6 @@ - inner2 = 0.0_r_def - inner2 = inner2 + lhs_lhs_input_inner_prod - inner2 = inner2 + x_x_input_inner_prod -- - ! Test the inner-product values for equality, allowing for the precision of the active variables - MachineTol = SPACING(MAX(ABS(inner1), ABS(inner2))) - relative_diff = ABS(inner1 - inner2) / MachineTol diff --git a/science/adjoint/patches/algorithm/adjt_transpose_matrix_vector_alg_mod.patch b/science/adjoint/patches/algorithm/adjt_transpose_matrix_vector_alg_mod.patch deleted file mode 100644 index 858e597e1..000000000 --- a/science/adjoint/patches/algorithm/adjt_transpose_matrix_vector_alg_mod.patch +++ /dev/null @@ -1,65 +0,0 @@ -@@ -14,8 +14,8 @@ - use fs_continuity_mod, only : w0 - use operator_mod, only : operator_type - use constants_mod, only : i_def, r_def -- use setop_random_kernel_mod, only : setop_random_kernel_type -- use log_mod, only : log_event, log_level_error, log_level_info, log_scratch_space -+ use setop_random_alg_mod, only : setop_random_alg -+ use log_mod, only : log_event, log_level_error, log_level_info, log_scratch_space, log_level_debug - real(kind=r_def), parameter :: overall_tolerance = 1500.0_r_def - type(mesh_type), pointer, intent(in) :: mesh - type(field_type), dimension(3), intent(in), optional :: chi -@@ -28,12 +28,15 @@ - type(field_type) :: x_input - real(kind=r_def) :: lhs_inner_prod - real(kind=r_def) :: x_inner_prod -+ real(kind=r_def) :: lhs_sf -+ real(kind=r_def) :: x_sf - real(kind=r_def) :: inner1 - real(kind=r_def) :: lhs_lhs_input_inner_prod - real(kind=r_def) :: x_x_input_inner_prod - real(kind=r_def) :: inner2 - real(kind=r_def) :: MachineTol - real(kind=r_def) :: relative_diff -+ real(kind=r_def), parameter :: eps = 1.0e-30_r_def - - vector_space_w0_ptr => function_space_collection%get_fs(mesh,element_order_h,element_order_v,w0) - call lhs%initialise(vector_space=vector_space_w0_ptr, name='lhs') -@@ -43,14 +46,24 @@ - call x_input%initialise(vector_space=vector_space_w0_ptr, name='x_input') - lhs_inner_prod = 0.0_r_def - x_inner_prod = 0.0_r_def -- - ! Initialise arguments and call the tangent-linear kernel. -- call invoke(setval_random(lhs), setval_x(lhs_input, lhs), setval_random(x), setval_x(x_input, x), & --&setop_random_kernel_type(matrix), transpose_matrix_vector_kernel_type(lhs, x, matrix), x_innerproduct_x(lhs_inner_prod, lhs), & -+ call invoke(setval_random(lhs), setval_x(lhs_input, lhs), setval_random(x), setval_x(x_input, x)) -+ call setop_random_alg(matrix) -+ call invoke(transpose_matrix_vector_kernel_type(lhs, x, matrix), x_innerproduct_x(lhs_inner_prod, lhs), & - &x_innerproduct_x(x_inner_prod, x)) -+ write( log_scratch_space, * ) "adjt_transpose_matrix_vector_alg inner products:" -+ call log_event( log_scratch_space, log_level_debug ) -+ write( log_scratch_space, * ) "lhs inner product = ", lhs_inner_prod -+ call log_event( log_scratch_space, log_level_debug ) -+ write( log_scratch_space, * ) "x inner product = ", x_inner_prod -+ call log_event( log_scratch_space, log_level_debug ) -+ lhs_sf = 1.0_r_def / (lhs_inner_prod + eps) -+ x_sf = 1.0_r_def / (x_inner_prod + eps) - inner1 = 0.0_r_def -- inner1 = inner1 + lhs_inner_prod -- inner1 = inner1 + x_inner_prod -+ inner1 = inner1 + lhs_inner_prod * lhs_sf -+ inner1 = inner1 + x_inner_prod * x_sf -+ call invoke( inc_a_times_X( lhs_sf, lhs ), & -+ inc_a_times_X( x_sf, x ) ) - lhs_lhs_input_inner_prod = 0.0_r_def - x_x_input_inner_prod = 0.0_r_def - call invoke(adj_transpose_matrix_vector_kernel_type(lhs, x, matrix), x_innerproduct_y(lhs_lhs_input_inner_prod, lhs, & -@@ -58,7 +71,6 @@ - inner2 = 0.0_r_def - inner2 = inner2 + lhs_lhs_input_inner_prod - inner2 = inner2 + x_x_input_inner_prod -- - ! Test the inner-product values for equality, allowing for the precision of the active variables - MachineTol = SPACING(MAX(ABS(inner1), ABS(inner2))) - relative_diff = ABS(inner1 - inner2) / MachineTol diff --git a/science/adjoint/patches/kernel/adj_dg_inc_matrix_vector_kernel_mod.patch b/science/adjoint/patches/kernel/adj_dg_inc_matrix_vector_kernel_mod.patch deleted file mode 100644 index ac07611db..000000000 --- a/science/adjoint/patches/kernel/adj_dg_inc_matrix_vector_kernel_mod.patch +++ /dev/null @@ -1,30 +0,0 @@ -@@ -1,23 +1,23 @@ - module adj_dg_inc_matrix_vector_kernel_mod - use argument_mod, only : any_discontinuous_space_1, any_space_1, arg_type, cell_column, gh_field, gh_operator, gh_read, & --&gh_readwrite, gh_real -+&gh_real, gh_inc - use constants_mod, only : i_def, r_double, r_single - use kernel_mod, only : kernel_type - implicit none -- interface dg_inc_matrix_vector_code -+ interface adj_dg_inc_matrix_vector_code - module procedure :: adj_dg_inc_matrix_vector_code_r_single, adj_dg_inc_matrix_vector_code_r_double -- end interface dg_inc_matrix_vector_code -+ end interface adj_dg_inc_matrix_vector_code - type, public, extends(kernel_type) :: adj_dg_inc_matrix_vector_kernel_type - type(ARG_TYPE) :: META_ARGS(3) = (/ & -- arg_type(gh_field, gh_real, gh_readwrite, any_discontinuous_space_1), & -- arg_type(gh_field, gh_real, gh_read, any_space_1), & -+ arg_type(gh_field, gh_real, gh_read, any_discontinuous_space_1), & -+ arg_type(gh_field, gh_real, gh_inc, any_space_1), & - arg_type(gh_operator, gh_real, gh_read, any_discontinuous_space_1, any_space_1)/) - INTEGER :: OPERATES_ON = cell_column - END TYPE adj_dg_inc_matrix_vector_kernel_type - - private - -- public :: dg_inc_matrix_vector_code -+ public :: adj_dg_inc_matrix_vector_code - - contains - subroutine adj_dg_inc_matrix_vector_code_r_single(cell, nlayers, lhs, x, ncell_3d, matrix, ndf1, undf1, map1, ndf2, undf2, map2) diff --git a/science/adjoint/patches/kernel/adj_dg_matrix_vector_kernel_mod.patch b/science/adjoint/patches/kernel/adj_dg_matrix_vector_kernel_mod.patch deleted file mode 100644 index ad1aa671a..000000000 --- a/science/adjoint/patches/kernel/adj_dg_matrix_vector_kernel_mod.patch +++ /dev/null @@ -1,47 +0,0 @@ -@@ -2,22 +2,22 @@ - use constants_mod, only : i_def, r_double, r_single - use kernel_mod, only : kernel_type - use argument_mod, only : any_discontinuous_space_1, any_space_1, arg_type, cell_column, gh_field, gh_operator, gh_read, gh_real, & --&gh_write -+&gh_readwrite, gh_inc - implicit none -- interface dg_matrix_vector_code -+ interface adj_dg_matrix_vector_code - module procedure :: adj_dg_matrix_vector_code_r_single, adj_dg_matrix_vector_code_r_double -- end interface dg_matrix_vector_code -+ end interface adj_dg_matrix_vector_code - type, public, extends(kernel_type) :: adj_dg_matrix_vector_kernel_type - type(ARG_TYPE) :: META_ARGS(3) = (/ & -- arg_type(gh_field, gh_real, gh_write, any_discontinuous_space_1), & -- arg_type(gh_field, gh_real, gh_read, any_space_1), & -+ arg_type(gh_field, gh_real, gh_readwrite, any_discontinuous_space_1), & -+ arg_type(gh_field, gh_real, gh_inc, any_space_1), & - arg_type(gh_operator, gh_real, gh_read, any_discontinuous_space_1, any_space_1)/) - INTEGER :: OPERATES_ON = cell_column - END TYPE adj_dg_matrix_vector_kernel_type - - private - -- public :: dg_matrix_vector_code -+ public :: adj_dg_matrix_vector_code - - contains - subroutine adj_dg_matrix_vector_code_r_single(cell, nlayers, lhs, x, ncell_3d, matrix, ndf1, undf1, map1, ndf2, undf2, map2) -@@ -56,7 +56,7 @@ - do df1 = ndf1, 1, -1 - i1 = map1(df1) - do idx = i1 + nl, i1, -1 -- lhs(idx) = 0.0 -+ lhs(idx) = 0.0_r_single - enddo - enddo - -@@ -97,7 +97,7 @@ - do df1 = ndf1, 1, -1 - i1 = map1(df1) - do idx = i1 + nl, i1, -1 -- lhs(idx) = 0.0 -+ lhs(idx) = 0.0_r_double - enddo - enddo - diff --git a/science/adjoint/patches/kernel/adj_matrix_vector_kernel_mod.patch b/science/adjoint/patches/kernel/adj_matrix_vector_kernel_mod.patch deleted file mode 100644 index d86ef04fc..000000000 --- a/science/adjoint/patches/kernel/adj_matrix_vector_kernel_mod.patch +++ /dev/null @@ -1,26 +0,0 @@ -@@ -3,20 +3,20 @@ - use constants_mod, only : i_def, r_double, r_single - use kernel_mod, only : kernel_type - implicit none -- interface matrix_vector_code -+ interface adj_matrix_vector_code - module procedure :: adj_matrix_vector_code_r_single, adj_matrix_vector_code_r_double -- end interface matrix_vector_code -+ end interface adj_matrix_vector_code - type, public, extends(kernel_type) :: adj_matrix_vector_kernel_type - type(ARG_TYPE) :: META_ARGS(3) = (/ & -- arg_type(gh_field, gh_real, gh_inc, any_space_1), & -- arg_type(gh_field, gh_real, gh_read, any_space_2), & -+ arg_type(gh_field, gh_real, gh_read, any_space_1), & -+ arg_type(gh_field, gh_real, gh_inc, any_space_2), & - arg_type(gh_operator, gh_real, gh_read, any_space_1, any_space_2)/) - INTEGER :: OPERATES_ON = cell_column - END TYPE adj_matrix_vector_kernel_type - - private - -- public :: matrix_vector_code -+ public :: adj_matrix_vector_code - - contains - subroutine adj_matrix_vector_code_r_single(cell, nlayers, lhs, x, ncell_3d, matrix, ndf1, undf1, map1, ndf2, undf2, map2) diff --git a/science/adjoint/patches/kernel/adj_transpose_matrix_vector_kernel_mod.patch b/science/adjoint/patches/kernel/adj_transpose_matrix_vector_kernel_mod.patch deleted file mode 100644 index 46e87a21e..000000000 --- a/science/adjoint/patches/kernel/adj_transpose_matrix_vector_kernel_mod.patch +++ /dev/null @@ -1,26 +0,0 @@ -@@ -3,20 +3,20 @@ - use constants_mod, only : i_def, r_double, r_single - use kernel_mod, only : kernel_type - implicit none -- interface transpose_matrix_vector_code -+ interface adj_transpose_matrix_vector_code - module procedure :: adj_transpose_matrix_vector_code_r_single, adj_transpose_matrix_vector_code_r_double -- end interface transpose_matrix_vector_code -+ end interface adj_transpose_matrix_vector_code - type, public, extends(kernel_type) :: adj_transpose_matrix_vector_kernel_type - type(ARG_TYPE) :: META_ARGS(3) = (/ & -- arg_type(gh_field, gh_real, gh_inc, any_space_1), & -- arg_type(gh_field, gh_real, gh_read, any_space_2), & -+ arg_type(gh_field, gh_real, gh_read, any_space_1), & -+ arg_type(gh_field, gh_real, gh_inc, any_space_2), & - arg_type(gh_operator, gh_real, gh_read, any_space_2, any_space_1)/) - INTEGER :: OPERATES_ON = cell_column - END TYPE adj_transpose_matrix_vector_kernel_type - - private - -- public :: transpose_matrix_vector_code -+ public :: adj_transpose_matrix_vector_code - - contains - subroutine adj_transpose_matrix_vector_code_r_single(cell, nlayers, lhs, x, ncell_3d, matrix, ndf1, undf1, map1, ndf2, undf2, & diff --git a/science/adjoint/patches/kernel/transpose_matrix_vector_kernel_mod.patch b/science/adjoint/patches/kernel/transpose_matrix_vector_kernel_mod.patch deleted file mode 100644 index 9a1fcbf0e..000000000 --- a/science/adjoint/patches/kernel/transpose_matrix_vector_kernel_mod.patch +++ /dev/null @@ -1,34 +0,0 @@ -@@ -88,13 +88,15 @@ - integer(kind=i_def) :: df, k, ik - real(kind=r_single), dimension(ndf2) :: x_e - real(kind=r_single), dimension(ndf1) :: lhs_e -+ real(kind=r_single), dimension(ndf1,ndf2) :: transposed_matrix - - do k = 0, nlayers-1 - do df = 1, ndf2 - x_e(df) = x(map2(df)+k) - end do - ik = (cell-1)*nlayers + k + 1 -- lhs_e = matmul(transpose(matrix(ik,:,:)),x_e) -+ transposed_matrix(:,:) = transpose(matrix(ik,:,:)) -+ lhs_e = matmul(transposed_matrix,x_e) - do df = 1,ndf1 - lhs(map1(df)+k) = lhs(map1(df)+k) + lhs_e(df) - end do -@@ -129,13 +131,15 @@ - integer(kind=i_def) :: df, k, ik - real(kind=r_double), dimension(ndf2) :: x_e - real(kind=r_double), dimension(ndf1) :: lhs_e -+ real(kind=r_double), dimension(ndf1,ndf2) :: transposed_matrix - - do k = 0, nlayers-1 - do df = 1, ndf2 - x_e(df) = x(map2(df)+k) - end do - ik = (cell-1)*nlayers + k + 1 -- lhs_e = matmul(transpose(matrix(ik,:,:)),x_e) -+ transposed_matrix(:,:) = transpose(matrix(ik,:,:)) -+ lhs_e = matmul(transposed_matrix,x_e) - do df = 1,ndf1 - lhs(map1(df)+k) = lhs(map1(df)+k) + lhs_e(df) - end do diff --git a/science/adjoint/source/kernel/algebra/adj_dg_inc_matrix_vector_kernel_mod.F90 b/science/adjoint/source/kernel/algebra/adj_dg_inc_matrix_vector_kernel_mod.F90 new file mode 100644 index 000000000..e2f91cc41 --- /dev/null +++ b/science/adjoint/source/kernel/algebra/adj_dg_inc_matrix_vector_kernel_mod.F90 @@ -0,0 +1,142 @@ +!----------------------------------------------------------------------------- +! (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 used. +!----------------------------------------------------------------------------- +!> @brief Adjoint of dg_inc_matrix_vector_kernel, for use on discontinous spaces +module adj_dg_inc_matrix_vector_kernel_mod + use argument_mod, only : arg_type, & + GH_FIELD, GH_OPERATOR, & + GH_REAL, GH_READ, & + GH_INC, ANY_SPACE_1, & + ANY_DISCONTINUOUS_SPACE_1, & + CELL_COLUMN + use constants_mod, only : i_def, r_single, r_double + use kernel_mod, only : kernel_type + + implicit none + + private + + !------------------------------------------------------------------------------- + ! Public types + !------------------------------------------------------------------------------- + type, public, extends(kernel_type) :: adj_dg_inc_matrix_vector_kernel_type + type(arg_type) :: meta_args(3) = (/ & + arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_1), & ! lhs + arg_type(GH_FIELD, GH_REAL, GH_INC, ANY_SPACE_1), & ! x + arg_type(GH_OPERATOR, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_1, & ! matrix + ANY_SPACE_1) & + /) + integer :: operates_on = CELL_COLUMN + end type adj_dg_inc_matrix_vector_kernel_type + + !------------------------------------------------------------------------------- + ! Contained functions/subroutines + !------------------------------------------------------------------------------- + public :: adj_dg_inc_matrix_vector_code + + ! Generic interface for real32 and real64 types + interface adj_dg_inc_matrix_vector_code + module procedure & + adj_dg_inc_matrix_vector_code_r_single, & + adj_dg_inc_matrix_vector_code_r_double + end interface adj_dg_inc_matrix_vector_code + + contains + + !> @brief Computes adjoint of lhs = lhs + matrix*x for discontinuous function spaces + !! real32 and real64 variants + !> @param[in] cell Horizontal cell index + !> @param[in] nlayers Number of layers + !> @param[inout] lhs Left hand side field of lhs = matrix*x + !> @param[in] x Right hand side field of lhs = matrix*x + !> @param[in] ncell_3d Total number of cells + !> @param[in] matrix Local matrix assembly form of the operator A + !> @param[in] ndf1 Number of degrees of freedom per cell for the output field + !> @param[in] undf1 Unique number of degrees of freedom for the output field + !> @param[in] map1 Dofmap for the cell at the base of the column for the + !! output field + !> @param[in] ndf2 Number of degrees of freedom per cell for the input field + !> @param[in] undf2 Unique number of degrees of freedom for the input field + !> @param[in] map2 Dofmap for the cell at the base of the column for the input + !! field + + ! R_SINGLE PRECISION + ! ================== + subroutine adj_dg_inc_matrix_vector_code_r_single(cell, & + nlayers, & + lhs, x, & + ncell_3d, & + matrix, & + ndf1, undf1, map1, & + ndf2, undf2, map2) + + implicit none + + ! Arguments + integer(kind=i_def), intent(in) :: cell, nlayers, ncell_3d + integer(kind=i_def), intent(in) :: undf1, ndf1 + integer(kind=i_def), intent(in) :: undf2, ndf2 + integer(kind=i_def), dimension(ndf1), intent(in) :: map1 + integer(kind=i_def), dimension(ndf2), intent(in) :: map2 + + real(kind=r_single), dimension(undf1), intent(in) :: lhs + real(kind=r_single), dimension(undf2), intent(inout) :: x + real(kind=r_single), dimension(ncell_3d,ndf1,ndf2), intent(in) :: matrix + + ! Internal variables + integer(kind=i_def) :: df, df2, ik, i1, i2, nl + + nl = nlayers - 1 + ik = cell * nlayers - nlayers + 1 + do df2 = ndf2, 1, -1 + i2 = map2(df2) + do df = ndf1, 1, -1 + i1 = map1(df) + x(i2:i2+nl) = x(i2:i2+nl) + matrix(ik:ik+nl,df,df2)*lhs(i1:i1+nl) + end do + end do + + end subroutine adj_dg_inc_matrix_vector_code_r_single + + ! R_DOUBLE PRECISION + ! ================== + subroutine adj_dg_inc_matrix_vector_code_r_double(cell, & + nlayers, & + lhs, x, & + ncell_3d, & + matrix, & + ndf1, undf1, map1, & + ndf2, undf2, map2) + + implicit none + + ! Arguments + integer(kind=i_def), intent(in) :: cell, nlayers, ncell_3d + integer(kind=i_def), intent(in) :: undf1, ndf1 + integer(kind=i_def), intent(in) :: undf2, ndf2 + integer(kind=i_def), dimension(ndf1), intent(in) :: map1 + integer(kind=i_def), dimension(ndf2), intent(in) :: map2 + + + real(kind=r_double), dimension(undf1), intent(in) :: lhs + real(kind=r_double), dimension(undf2), intent(inout) :: x + real(kind=r_double), dimension(ncell_3d,ndf1,ndf2), intent(in) :: matrix + + ! Internal variables + integer(kind=i_def) :: df, df2, ik, i1, i2, nl + + nl = nlayers - 1 + ik = cell * nlayers - nlayers + 1 + do df2 = ndf2, 1, -1 + i2 = map2(df2) + do df = ndf1, 1, -1 + i1 = map1(df) + x(i2:i2+nl) = x(i2:i2+nl) + matrix(ik:ik+nl,df,df2)*lhs(i1:i1+nl) + end do + end do + + end subroutine adj_dg_inc_matrix_vector_code_r_double + +end module adj_dg_inc_matrix_vector_kernel_mod diff --git a/science/adjoint/source/kernel/algebra/adj_dg_matrix_vector_kernel_mod.F90 b/science/adjoint/source/kernel/algebra/adj_dg_matrix_vector_kernel_mod.F90 new file mode 100644 index 000000000..cb8714f78 --- /dev/null +++ b/science/adjoint/source/kernel/algebra/adj_dg_matrix_vector_kernel_mod.F90 @@ -0,0 +1,156 @@ +!----------------------------------------------------------------------------- +! (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 used. +!----------------------------------------------------------------------------- +!> @brief Adjoint of dg_matrix_vector_kernel. This version is for use on +!> discontinous spaces and over-writes the output field entirely, +!> hence can only be used for W3 spaces +module adj_dg_matrix_vector_kernel_mod + + use constants_mod, only : i_def, r_single, r_double + use kernel_mod, only : kernel_type + use argument_mod, only : arg_type, & + GH_FIELD, GH_OPERATOR, & + GH_READ, & + GH_INC, GH_READWRITE, & + GH_REAL, ANY_SPACE_1, & + ANY_DISCONTINUOUS_SPACE_1, & + CELL_COLUMN + + implicit none + + private + + !--------------------------------------------------------------------------- + ! Public types + !--------------------------------------------------------------------------- + + type, public, extends(kernel_type) :: adj_dg_matrix_vector_kernel_type + private + type(arg_type) :: meta_args(3) = (/ & + arg_type(GH_FIELD, GH_REAL, GH_READWRITE, ANY_DISCONTINUOUS_SPACE_1), & ! lhs + arg_type(GH_FIELD, GH_REAL, GH_INC, ANY_SPACE_1), & ! x + arg_type(GH_OPERATOR, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_1, & + ANY_SPACE_1) & ! matrix + /) + integer :: operates_on = CELL_COLUMN + end type adj_dg_matrix_vector_kernel_type + + !--------------------------------------------------------------------------- + ! Contained functions/subroutines + !--------------------------------------------------------------------------- + public :: adj_dg_matrix_vector_code + + ! Generic interface for real32 and real64 types + interface adj_dg_matrix_vector_code + module procedure & + adj_dg_matrix_vector_code_r_single, & + adj_dg_matrix_vector_code_r_double + end interface + + contains + + !> @brief Computes adjoint of lhs = matrix*x for discontinuous function spaces + !> @param[in] cell Horizontal cell index + !> @param[in] nlayers Number of layers + !> @param[in,out] lhs Left hand side field of lhs = matrix*x + !> @param[in,out] x Right hand side field of lhs = matrix*x + !> @param[in] ncell_3d Total number of cells + !> @param[in] matrix Matrix values in LMA form + !> @param[in] ndf1 Number of degrees of freedom per cell for the output field + !> @param[in] undf1 Unique number of degrees of freedom for the output field + !> @param[in] map1 Dofmap for the cell at the base of the column for the output field + !> @param[in] ndf2 Number of degrees of freedom per cell for the input field + !> @param[in] undf2 Unique number of degrees of freedom for the input field + !> @param[in] map2 Dofmap for the cell at the base of the column for the input field + + ! R_SINGLE PRECISION + ! ================== + subroutine adj_dg_matrix_vector_code_r_single(cell, & + nlayers, & + lhs, x, & + ncell_3d, & + matrix, & + ndf1, undf1, map1, & + ndf2, undf2, map2) + + implicit none + + ! Arguments + integer(kind=i_def), intent(in) :: cell, nlayers, ncell_3d + integer(kind=i_def), intent(in) :: undf1, ndf1 + integer(kind=i_def), intent(in) :: undf2, ndf2 + integer(kind=i_def), dimension(ndf1), intent(in) :: map1 + integer(kind=i_def), dimension(ndf2), intent(in) :: map2 + + real(kind=r_single), dimension(undf2), intent(inout) :: x + real(kind=r_single), dimension(undf1), intent(inout) :: lhs + real(kind=r_single), dimension(ncell_3d,ndf1,ndf2), intent(in) :: matrix + + ! Internal variables + integer(kind=i_def) :: df1, df2, ik, i1, i2, nl + + nl = nlayers - 1 + + ik = cell * nlayers - nlayers + 1 + do df2 = ndf2, 1, -1 + i2 = map2(df2) + do df1 = ndf1, 1, -1 + i1 = map1(df1) + x(i2:i2+nl) = x(i2:i2+nl) + matrix(ik:ik+nl,df1,df2)*lhs(i1:i1+nl) + end do + end do + + do df1 = ndf1, 1, -1 + i1 = map1(df1) + lhs(i1:i1+nl) = 0.0_r_single + end do + + end subroutine adj_dg_matrix_vector_code_r_single + + ! R_DOUBLE PRECISION + ! ================== + subroutine adj_dg_matrix_vector_code_r_double(cell, & + nlayers, & + lhs, x, & + ncell_3d, & + matrix, & + ndf1, undf1, map1, & + ndf2, undf2, map2) + + implicit none + + ! Arguments + integer(kind=i_def), intent(in) :: cell, nlayers, ncell_3d + integer(kind=i_def), intent(in) :: undf1, ndf1 + integer(kind=i_def), intent(in) :: undf2, ndf2 + integer(kind=i_def), dimension(ndf1), intent(in) :: map1 + integer(kind=i_def), dimension(ndf2), intent(in) :: map2 + + real(kind=r_double), dimension(undf2), intent(inout) :: x + real(kind=r_double), dimension(undf1), intent(inout) :: lhs + real(kind=r_double), dimension(ncell_3d,ndf1,ndf2), intent(in) :: matrix + + ! Internal variables + integer(kind=i_def) :: df1, df2, ik, i1, i2, nl + + nl = nlayers - 1 + + ik = cell * nlayers - nlayers + 1 + do df2 = ndf2, 1, -1 + i2 = map2(df2) + do df1 = ndf1, 1, -1 + i1 = map1(df1) + x(i2:i2+nl) = x(i2:i2+nl) + matrix(ik:ik+nl,df1,df2)*lhs(i1:i1+nl) + end do + end do + + do df1 = ndf1, 1, -1 + i1 = map1(df1) + lhs(i1:i1+nl) = 0.0_r_double + end do + + end subroutine adj_dg_matrix_vector_code_r_double + +end module adj_dg_matrix_vector_kernel_mod diff --git a/science/adjoint/source/kernel/algebra/adj_matrix_vector_kernel_mod.F90 b/science/adjoint/source/kernel/algebra/adj_matrix_vector_kernel_mod.F90 new file mode 100644 index 000000000..8af235360 --- /dev/null +++ b/science/adjoint/source/kernel/algebra/adj_matrix_vector_kernel_mod.F90 @@ -0,0 +1,141 @@ +!----------------------------------------------------------------------------- +! (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 used. +!----------------------------------------------------------------------------- +!> @brief Adjoint kernel of matrix vector operation for use on continuous spaces +module adj_matrix_vector_kernel_mod + use argument_mod, only : arg_type, & + GH_FIELD, GH_REAL, & + GH_OPERATOR, & + GH_INC, GH_READ, & + ANY_SPACE_1, ANY_SPACE_2, & + CELL_COLUMN + use constants_mod, only : i_def, r_double, r_single + use kernel_mod, only : kernel_type + + implicit none + + private + + !------------------------------------------------------------------------------- + ! Public types + !------------------------------------------------------------------------------- + type, public, extends(kernel_type) :: adj_matrix_vector_kernel_type + private + type(arg_type) :: meta_args(3) = (/ & + arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_SPACE_1), & ! lhs + arg_type(GH_FIELD, GH_REAL, GH_INC, ANY_SPACE_2), & ! x + arg_type(GH_OPERATOR, GH_REAL, GH_READ, ANY_SPACE_1, ANY_SPACE_2) & ! matrix + /) + integer :: operates_on = CELL_COLUMN + end type adj_matrix_vector_kernel_type + + !------------------------------------------------------------------------------- + ! Contained functions/subroutines + !------------------------------------------------------------------------------- + public :: adj_matrix_vector_code + + ! Generic interface for real32 and real64 types + interface adj_matrix_vector_code + module procedure & + adj_matrix_vector_code_r_single, & + adj_matrix_vector_code_r_double + end interface adj_matrix_vector_code + +contains + + !> @brief Computes adjoint of lhs = lhs + matrix*x + !! (x = x + matrix*lhs) + !! with real32 and real64 variants + !> @param[in] cell Horizontal cell index + !> @param[in] nlayers Number of layers + !> @param[in] lhs Left hand side field of lhs = matrix*x + !> @param[inout] x Right hand side field of lhs = matrix*x + !> @param[in] ncell_3d Total number of cells + !> @param[in] matrix Local matrix assembly form of the operator A + !> @param[in] ndf1 Number of degrees of freedom per cell for the output field + !> @param[in] undf1 Unique number of degrees of freedom for the output field + !> @param[in] map1 Dofmap for the cell at the base of the column for the + !! output field + !> @param[in] ndf2 Number of degrees of freedom per cell for the input field + !> @param[in] undf2 Unique number of degrees of freedom for the input field + !> @param[in] map2 Dofmap for the cell at the base of the column for the input + !! field + + ! R_SINGLE PRECISION + ! ================== + subroutine adj_matrix_vector_code_r_single(cell, & + nlayers, & + lhs, & + x, & + ncell_3d, & + matrix, & + ndf1, undf1, map1, & + ndf2, undf2, map2) + implicit none + + ! Arguments + integer(kind=i_def), intent(in) :: cell, nlayers, ncell_3d + integer(kind=i_def), intent(in) :: undf1, ndf1 + integer(kind=i_def), intent(in) :: undf2, ndf2 + integer(kind=i_def), dimension(ndf1), intent(in) :: map1 + integer(kind=i_def), dimension(ndf2), intent(in) :: map2 + + real(kind=r_single), dimension(undf2), intent(inout) :: x + real(kind=r_single), dimension(undf1), intent(in) :: lhs + real(kind=r_single), dimension(ncell_3d,ndf1,ndf2), intent(in) :: matrix + + ! Internal variables + integer(kind=i_def) :: df, ik, df2, i1, i2, nl + + nl = nlayers - 1 + ik = cell * nlayers - nlayers + 1 + do df2 = ndf2, 1, -1 + i2 = map2(df2) + do df = ndf1, 1, -1 + i1 = map1(df) + x(i2:i2+nl) = x(i2:i2+nl) + matrix(ik:ik+nl,df,df2)*lhs(i1:i1+nl) + end do + end do + + end subroutine adj_matrix_vector_code_r_single + + ! R_DOUBLE PRECISION + ! ================== + subroutine adj_matrix_vector_code_r_double(cell, & + nlayers, & + lhs, x, & + ncell_3d, & + matrix, & + ndf1, undf1, map1, & + ndf2, undf2, map2) + + implicit none + + ! Arguments + integer(kind=i_def), intent(in) :: cell, nlayers, ncell_3d + integer(kind=i_def), intent(in) :: undf1, ndf1 + integer(kind=i_def), intent(in) :: undf2, ndf2 + integer(kind=i_def), dimension(ndf1), intent(in) :: map1 + integer(kind=i_def), dimension(ndf2), intent(in) :: map2 + + real(kind=r_double), dimension(undf2), intent(inout) :: x + real(kind=r_double), dimension(undf1), intent(in) :: lhs + real(kind=r_double), dimension(ncell_3d,ndf1,ndf2), intent(in) :: matrix + + integer(kind=i_def) :: df, ik, df2, i1, i2, nl + + nl = nlayers-1 + ik = (cell-1)*nlayers +1 + do df2 = 1, ndf2 + i2 = map2(df2) + do df = 1, ndf1 + i1 = map1(df) + x(i2:i2+nl) = x(i2:i2+nl) + matrix(ik:ik+nl,df,df2)*lhs(i1:i1+nl) + end do + end do + + end subroutine adj_matrix_vector_code_r_double + +end module adj_matrix_vector_kernel_mod diff --git a/science/adjoint/source/kernel/algebra/adj_transpose_matrix_vector_kernel_mod.F90 b/science/adjoint/source/kernel/algebra/adj_transpose_matrix_vector_kernel_mod.F90 new file mode 100644 index 000000000..96696da3e --- /dev/null +++ b/science/adjoint/source/kernel/algebra/adj_transpose_matrix_vector_kernel_mod.F90 @@ -0,0 +1,144 @@ +!----------------------------------------------------------------------------- +! (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 used. +!----------------------------------------------------------------------------- +!> @brief Adjoint of transpose_matrix_vector_kernel +module adj_transpose_matrix_vector_kernel_mod + use argument_mod, only : arg_type, & + GH_FIELD, GH_OPERATOR, & + GH_REAL, GH_READ, & + GH_INC, & + ANY_SPACE_1, ANY_SPACE_2, & + CELL_COLUMN + use constants_mod, only : i_def, r_single, r_double + use kernel_mod, only : kernel_type + + implicit none + + private + + !------------------------------------------------------------------------------- + ! Public types + !------------------------------------------------------------------------------- + type, public, extends(kernel_type) :: adj_transpose_matrix_vector_kernel_type + private + type(arg_type) :: meta_args(3) = (/ & + arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_SPACE_1), & + arg_type(GH_FIELD, GH_REAL, GH_INC, ANY_SPACE_2), & + arg_type(GH_OPERATOR, GH_REAL, GH_READ, ANY_SPACE_2, ANY_SPACE_1) & + /) + integer :: operates_on = CELL_COLUMN + end type + + !------------------------------------------------------------------------------- + ! Contained functions/subroutines + !------------------------------------------------------------------------------- + public :: adj_transpose_matrix_vector_code + + ! Generic interface for real32 and real64 types + interface adj_transpose_matrix_vector_code + module procedure & + adj_transpose_matrix_vector_code_r_single, & + adj_transpose_matrix_vector_code_r_double + end interface + + contains + + !> @brief Computes adjoint of lhs = matrix^T*x where matrix^T is the transpose of the matrix + !! @param[in] cell Horizontal cell index + !! @param[in] nlayers Number of layers + !! @param[in] lhs Left hand side field of lhs = matrix^T*x + !! @param[in,out] x Right hand side field of lhs = matrix^T*x + !! @param[in] ncell_3d Total number of cells + !! @param[in] matrix Matrix values + !! @param[in] ndf1 Number of degrees of freedom per cell for the output field + !! @param[in] undf1 Unique number of degrees of freedom for the output field + !! @param[in] map1 Dofmap for the cell at the base of the column for the output field + !! @param[in] ndf2 Number of degrees of freedom per cell for the input field + !! @param[in] undf2 Unique number of degrees of freedom for the input field + !! @param[in] map2 Dofmap for the cell at the base of the column for the input field + + ! R_SINGLE PRECISION + ! ================== + subroutine adj_transpose_matrix_vector_code_r_single(cell, & + nlayers, & + lhs, x, & + ncell_3d, & + matrix, & + ndf1, undf1, map1, & + ndf2, undf2, map2) + + implicit none + + ! Arguments + integer(kind=i_def), intent(in) :: cell, nlayers, ncell_3d + integer(kind=i_def), intent(in) :: undf1, ndf1 + integer(kind=i_def), intent(in) :: undf2, ndf2 + integer(kind=i_def), dimension(ndf1), intent(in) :: map1 + integer(kind=i_def), dimension(ndf2), intent(in) :: map2 + + real(kind=r_single), dimension(undf1), intent(in) :: lhs + real(kind=r_single), dimension(undf2), intent(inout) :: x + real(kind=r_single), dimension(ncell_3d,ndf2,ndf1), intent(in) :: matrix + + ! Internal variables + integer(kind=i_def) :: df, k, ik + real(kind=r_single), dimension(ndf2) :: x_e + real(kind=r_single), dimension(ndf1) :: lhs_e + + do k = nlayers - 1, 0, -1 + ik = cell * nlayers + k - nlayers + 1 + do df = ndf1, 1, -1 + lhs_e(df) = lhs(map1(df) + k) + end do + x_e = matmul(matrix(ik,:,:), lhs_e) + do df = ndf2, 1, -1 + x(map2(df) + k) = x(map2(df) + k) + x_e(df) + end do + end do + + end subroutine adj_transpose_matrix_vector_code_r_single + + ! R_DOUBLE PRECISION + ! ================== + subroutine adj_transpose_matrix_vector_code_r_double(cell, & + nlayers, & + lhs, x, & + ncell_3d, & + matrix, & + ndf1, undf1, map1, & + ndf2, undf2, map2) + + implicit none + + ! Arguments + integer(kind=i_def), intent(in) :: cell, nlayers, ncell_3d + integer(kind=i_def), intent(in) :: undf1, ndf1 + integer(kind=i_def), intent(in) :: undf2, ndf2 + integer(kind=i_def), dimension(ndf1), intent(in) :: map1 + integer(kind=i_def), dimension(ndf2), intent(in) :: map2 + + real(kind=r_double), dimension(undf1), intent(in) :: lhs + real(kind=r_double), dimension(undf2), intent(inout) :: x + real(kind=r_double), dimension(ncell_3d,ndf2,ndf1), intent(in) :: matrix + + ! Internal variables + integer(kind=i_def) :: df, k, ik + real(kind=r_double), dimension(ndf2) :: x_e + real(kind=r_double), dimension(ndf1) :: lhs_e + + do k = nlayers - 1, 0, -1 + ik = cell * nlayers + k - nlayers + 1 + do df = ndf1, 1, -1 + lhs_e(df) = lhs(map1(df) + k) + end do + x_e = matmul(matrix(ik,:,:), lhs_e) + do df = ndf2, 1, -1 + x(map2(df) + k) = x(map2(df) + k) + x_e(df) + end do + end do + + end subroutine adj_transpose_matrix_vector_code_r_double + +end module adj_transpose_matrix_vector_kernel_mod