diff --git a/components/science/source/algorithm/sci_geometric_constants_mod.x90 b/components/science/source/algorithm/sci_geometric_constants_mod.x90 index 744700203..85c9430fd 100644 --- a/components/science/source/algorithm/sci_geometric_constants_mod.x90 +++ b/components/science/source/algorithm/sci_geometric_constants_mod.x90 @@ -1162,8 +1162,12 @@ contains !> @return A height field function get_height_fe(space_id, mesh_id) result(height) - use sci_get_height_kernel_mod, only: get_height_kernel_type - use planet_config_mod, only: scaled_radius + use sci_height_continuous_kernel_mod, only: height_continuous_kernel_type + use sci_height_discontinuous_kernel_mod, & + only: height_discontinuous_kernel_type + use base_mesh_config_mod, only: geometry + use finite_element_config_mod, only: coord_system + use planet_config_mod, only: scaled_radius implicit none @@ -1175,6 +1179,9 @@ contains type(function_space_type), pointer :: space type(field_type), pointer :: chi(:) type(field_type), pointer :: height + type(field_type) :: rmultiplicity + type(field_type) :: nodal_multiplicity + type(field_type) :: ones character(len=str_def) :: inventory_name ! If running at lowest order, use finite volume @@ -1219,11 +1226,37 @@ contains if ( subroutine_timers ) call timer('runtime_constants.geometric') - space => function_space_collection%get_fs(mesh, element_order_h, & - element_order_v, space_id) + space => function_space_collection%get_fs( & + mesh, element_order_h, element_order_v, space_id & + ) call inventory%add_field(height, space, mesh) - call invoke( get_height_kernel_type(height, chi, scaled_radius) ) + select case (space_id) + case (W3, Wtheta) + call invoke( & + height_discontinuous_kernel_type( & + height, chi, geometry, coord_system, scaled_radius & + ) & + ) + + case default + ! Can't import multiplicity, so must calculate it + call ones%initialise( space ) + call nodal_multiplicity%initialise( space ) + call rmultiplicity%initialise( space ) + + call invoke( & + setval_c(ones, 1.0_r_def), & + setval_c(nodal_multiplicity, 0.0_r_def), & + multiplicity_kernel_type(nodal_multiplicity), & + X_divideby_Y(rmultiplicity, ones, nodal_multiplicity), & + setval_c(height, 0.0_r_def), & + height_continuous_kernel_type( & + height, chi, rmultiplicity, & + geometry, coord_system, scaled_radius & + ) & + ) + end select if ( subroutine_timers ) call timer('runtime_constants.geometric') else @@ -1238,8 +1271,12 @@ contains !> @return A height field function get_height_fv(space_id, mesh_id) result(height) - use sci_get_height_kernel_mod, only: get_height_kernel_type - use planet_config_mod, only: scaled_radius + use sci_height_continuous_kernel_mod, only: height_continuous_kernel_type + use sci_height_discontinuous_kernel_mod, & + only: height_discontinuous_kernel_type + use base_mesh_config_mod, only: geometry + use finite_element_config_mod, only: coord_system + use planet_config_mod, only: scaled_radius implicit none @@ -1251,6 +1288,9 @@ contains type(function_space_type), pointer :: space type(field_type), pointer :: chi(:) type(field_type), pointer :: height + type(field_type) :: rmultiplicity + type(field_type) :: nodal_multiplicity + type(field_type) :: ones character(len=str_def) :: inventory_name ! Determine inventory based on space @@ -1292,7 +1332,32 @@ contains space => function_space_collection%get_fs(mesh, 0, 0, space_id) call inventory%add_field(height, space, mesh) - call invoke( get_height_kernel_type(height, chi, scaled_radius) ) + select case (space_id) + case (W3, Wtheta) + call invoke( & + height_discontinuous_kernel_type( & + height, chi, geometry, coord_system, scaled_radius & + ) & + ) + + case default + ! Can't import multiplicity, so must calculate it + call ones%initialise( space ) + call nodal_multiplicity%initialise( space ) + call rmultiplicity%initialise( space ) + + call invoke( & + setval_c(ones, 1.0_r_def), & + setval_c(nodal_multiplicity, 0.0_r_def), & + multiplicity_kernel_type(nodal_multiplicity), & + X_divideby_Y(rmultiplicity, ones, nodal_multiplicity), & + setval_c(height, 0.0_r_def), & + height_continuous_kernel_type( & + height, chi, rmultiplicity, & + geometry, coord_system, scaled_radius & + ) & + ) + end select if ( subroutine_timers ) call timer('runtime_constants.geometric') else diff --git a/components/science/source/kernel/geometry/sci_get_height_kernel_mod.F90 b/components/science/source/kernel/geometry/sci_get_height_kernel_mod.F90 deleted file mode 100644 index 9778fdc6f..000000000 --- a/components/science/source/kernel/geometry/sci_get_height_kernel_mod.F90 +++ /dev/null @@ -1,144 +0,0 @@ -!----------------------------------------------------------------------------- -! Copyright (c) 2017, Met Office, on behalf of HMSO and Queen's Printer -! For further details please refer to the file LICENCE which you -! should have received as part of this distribution. -!----------------------------------------------------------------------------- -!> @brief Returns a height field (r or z) from the chi array. -!> -!> @details Returns a height field (r or z) from the chi array -!> -module sci_get_height_kernel_mod - - use argument_mod, only: arg_type, func_type, & - GH_FIELD, GH_REAL, & - GH_SCALAR, & - GH_WRITE, GH_READ, GH_INC, & - ANY_DISCONTINUOUS_SPACE_1, & - ANY_SPACE_9, GH_BASIS, & - CELL_COLUMN, GH_EVALUATOR - use base_mesh_config_mod, only: geometry, & - geometry_spherical - use constants_mod, only: r_def, i_def - use finite_element_config_mod, only: coord_system, & - coord_system_xyz - use kernel_mod, only: kernel_type - - implicit none - private - - !--------------------------------------------------------------------------- - ! Public types - !--------------------------------------------------------------------------- - !> The type declaration for the kernel. Contains the metadata needed by the - !> Psy layer. - !> - type, public, extends(kernel_type) :: get_height_kernel_type - private - type(arg_type) :: meta_args(3) = (/ & - arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1), & - arg_type(GH_FIELD*3, GH_REAL, GH_READ, ANY_SPACE_9), & - arg_type(GH_SCALAR, GH_REAL, GH_READ) & - /) - type(func_type) :: meta_funcs(1) = (/ & - func_type(ANY_SPACE_9, GH_BASIS) & - /) - integer :: operates_on = CELL_COLUMN - integer :: gh_shape = GH_EVALUATOR - contains - procedure, nopass :: get_height_code - end type - - !--------------------------------------------------------------------------- - ! Contained functions/subroutines - !--------------------------------------------------------------------------- - public :: get_height_code - -contains - -!> @brief Returns a height field (r or z) from the chi array -!> Will only work at lowest order for now -!! @param[in] nlayers Number of layers -!! @param[in,out] height The height field -!! @param[in] chi_1 1st component of the coordinate -!! @param[in] chi_2 2nd component of the coordinate -!! @param[in] chi_3 3rd component of the coordinate -!! @param[in] planet_radius The planet radius -!! @param[in] ndf_x Number of degrees of freedom per cell for height -!! @param[in] undf_x Number of unique degrees of freedom for height -!! @param[in] map_x Dofmap for the cell at the base of the column for height -!! @param[in] ndf_chi The number of degrees of freedom per cell for chi -!! @param[in] undf_chi The number of unique degrees of freedom for chi -!! @param[in] map_chi Dofmap for the cell at the base of the column for chi -!! @param[in] basis_chi Basis functions evaluated at nodal points for height -subroutine get_height_code(nlayers, & - height, & - chi_1, chi_2, chi_3, & - planet_radius, & - ndf_x, undf_x, map_x, & - ndf_chi, undf_chi, map_chi, & - basis_chi & - ) - implicit none - - ! Arguments - integer(kind=i_def), intent(in) :: nlayers - - integer(kind=i_def), intent(in) :: ndf_x, undf_x - integer(kind=i_def), intent(in) :: ndf_chi, undf_chi - real(kind=r_def), dimension(undf_x), intent(inout) :: height - real(kind=r_def), dimension(undf_chi), intent(in) :: chi_1, chi_2, chi_3 - real(kind=r_def), intent(in) :: planet_radius - - integer(kind=i_def), dimension(ndf_x), intent(in) :: map_x - integer(kind=i_def), dimension(ndf_chi), intent(in) :: map_chi - real(kind=r_def), dimension(1,ndf_chi,ndf_x), intent(in) :: basis_chi - - ! Internal variables - integer(kind=i_def) :: df_chi, df_x, k - real(kind=r_def) :: coord(3), coord_radius, height_at_dof - - if ( (geometry == geometry_spherical) .and. & - (coord_system == coord_system_xyz) ) then - ! NB This will result in the height above - ! the spherical representation of the planet - ! but not necessarily the height above the bottom - ! of the mesh - ! This should be reviewed with ticket #562 - - do k = 0, nlayers-1 - do df_x = 1, ndf_x - coord(:) = 0.0_r_def - do df_chi = 1, ndf_chi - coord(1) = coord(1) + chi_1(map_chi(df_chi)+k)*basis_chi(1,df_chi,df_x) - coord(2) = coord(2) + chi_2(map_chi(df_chi)+k)*basis_chi(1,df_chi,df_x) - coord(3) = coord(3) + chi_3(map_chi(df_chi)+k)*basis_chi(1,df_chi,df_x) - end do - - coord_radius = sqrt(coord(1)**2 + coord(2)**2 + coord(3)**2) - - height( map_x(df_x) + k ) = coord_radius - planet_radius - - end do - end do - - else - - ! Either the domain is Cartesian or we are using a non-Cartesian spherical - ! coordinate system. In both these cases, chi_3 will be the height. - do k = 0, nlayers-1 - do df_x = 1, ndf_x - height_at_dof = 0.0_r_def - do df_chi = 1, ndf_chi - height_at_dof = height_at_dof + & - chi_3(map_chi(df_chi)+k)*basis_chi(1,df_chi,df_x) - end do - - height( map_x(df_x) + k ) = height_at_dof - - end do - end do - end if - -end subroutine get_height_code - -end module sci_get_height_kernel_mod diff --git a/components/science/source/kernel/geometry/sci_height_continuous_kernel_mod.F90 b/components/science/source/kernel/geometry/sci_height_continuous_kernel_mod.F90 new file mode 100644 index 000000000..4964cb6aa --- /dev/null +++ b/components/science/source/kernel/geometry/sci_height_continuous_kernel_mod.F90 @@ -0,0 +1,179 @@ +!------------------------------------------------------------------------------- +! (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 Computes the height above the Earth's mean radius at the degrees of +!! freedom for a continuous function space. +!> @details Uses the chi coordinate fields to determine the height of nodes +!! above the Earth's mean radius, for continuous function spaces. +module sci_height_continuous_kernel_mod + + use argument_mod, only: arg_type, func_type, & + GH_FIELD, GH_SCALAR, & + GH_REAL, GH_INTEGER, & + GH_READ, GH_INC, & + ANY_SPACE_1, ANY_SPACE_9, & + CELL_COLUMN, GH_BASIS, GH_EVALUATOR + use base_mesh_config_mod, only: geometry_spherical + use constants_mod, only: r_def, i_def, l_def + use finite_element_config_mod, only: coord_system_xyz + use kernel_mod, only: kernel_type + + implicit none + private + + !----------------------------------------------------------------------------- + ! Public types + !----------------------------------------------------------------------------- + !> The type declaration for the kernel. Contains the metadata needed by the + !! Psy layer. + type, public, extends(kernel_type) :: height_continuous_kernel_type + private + type(arg_type) :: meta_args(6) = (/ & + arg_type(GH_FIELD, GH_REAL, GH_INC, ANY_SPACE_1), & + arg_type(GH_FIELD*3, GH_REAL, GH_READ, ANY_SPACE_9), & + arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_SPACE_1), & + arg_type(GH_SCALAR, GH_INTEGER, GH_READ), & + arg_type(GH_SCALAR, GH_INTEGER, GH_READ), & + arg_type(GH_SCALAR, GH_REAL, GH_READ) & + /) + type(func_type) :: meta_funcs(1) = (/ & + func_type(ANY_SPACE_9, GH_BASIS) & + /) + integer :: operates_on = CELL_COLUMN + integer :: gh_shape = GH_EVALUATOR + contains + procedure, nopass :: height_continuous_code + end type + + !----------------------------------------------------------------------------- + ! Contained functions/subroutines + !----------------------------------------------------------------------------- + public :: height_continuous_code + +contains + +!> @brief Computes the height above the Earth's mean radius at the degrees of +!! freedom for a continuous function space. +!> @param[in] nlayers Number of layers in the mesh +!> @param[in,out] height The height field to compute +!> @param[in] chi_1 1st component of the coordinate fields +!> @param[in] chi_2 2nd component of the coordinate fields +!> @param[in] chi_3 3rd component of the coordinate fields +!> @param[in] rmultiplicity Reciprocal of DoF multiplicity for height field +!> @param[in] geometry The geometry of the domain +!> @param[in] coord_system The coordinate system of the domain +!> @param[in] planet_radius The planet radius +!> @param[in] ndf_h Num DoFs per cell for height field +!> @param[in] undf_h Num DoFs in this partition for height field +!> @param[in] map_h DoF index map for height field +!> @param[in] ndf_chi Num DoFs per cell for chi fields +!> @param[in] undf_chi Num DoFs in this partition for chi fields +!> @param[in] map_chi DoF index map for chi fields +!> @param[in] basis_chi Chi basis functions evaluated at height DoFs +subroutine height_continuous_code( & + nlayers, & + height, & + chi_1, chi_2, chi_3, & + rmultiplicity, & + geometry, coord_system, planet_radius, & + ndf_h, undf_h, map_h, & + ndf_chi, undf_chi, map_chi, & + basis_chi & +) + implicit none + + ! Arguments + integer(kind=i_def), intent(in) :: nlayers + integer(kind=i_def), intent(in) :: ndf_h, undf_h + integer(kind=i_def), intent(in) :: ndf_chi, undf_chi + integer(kind=i_def), intent(in) :: geometry, coord_system + integer(kind=i_def), intent(in) :: map_h(ndf_h) + integer(kind=i_def), intent(in) :: map_chi(ndf_chi) + real(kind=r_def), intent(in) :: planet_radius + real(kind=r_def), intent(in) :: basis_chi(1,ndf_chi,ndf_h) + real(kind=r_def), intent(inout) :: height(undf_h) + real(kind=r_def), intent(in) :: rmultiplicity(undf_h) + real(kind=r_def), intent(in) :: chi_1(undf_chi) + real(kind=r_def), intent(in) :: chi_2(undf_chi) + real(kind=r_def), intent(in) :: chi_3(undf_chi) + + ! Internal variables + logical(kind=l_def) :: is_spherical_xyz + integer(kind=i_def) :: df_chi, df_h + integer(kind=i_def) :: h_b_idx, h_t_idx, chi_b_idx, chi_t_idx + real(kind=r_def) :: coord_1(nlayers), coord_2(nlayers), coord_3(nlayers) + real(kind=r_def) :: coord_radius(nlayers), basis_val + + ! Two cases: if the geometry is spherical with geocentric Cartesian coords, + ! all three chi components are needed to compute the height. Otherwise, only + ! the last chi component is needed, as this is already the vertical coordinate + + is_spherical_xyz = ( & + geometry == geometry_spherical .and. coord_system == coord_system_xyz & + ) + + ! -------------------------------------------------------------------------- ! + ! Cartesian system and spherical + ! -------------------------------------------------------------------------- ! + if (is_spherical_xyz) then + + do df_h = 1, ndf_h + ! Indices for bottom and top height DoFs in this column + h_b_idx = map_h(df_h) + h_t_idx = map_h(df_h) + nlayers - 1 + + ! Determine coordinates at this DoF using chi basis functions + coord_1(:) = 0.0_r_def + coord_2(:) = 0.0_r_def + coord_3(:) = 0.0_r_def + do df_chi = 1, ndf_chi + ! Indices for bottom and top chi DoFs in this column + chi_b_idx = map_chi(df_chi) + chi_t_idx = map_chi(df_chi) + nlayers - 1 + basis_val = basis_chi(1, df_chi, df_h) + + coord_1(:) = coord_1(:) + chi_1(chi_b_idx:chi_t_idx)*basis_val + coord_2(:) = coord_2(:) + chi_2(chi_b_idx:chi_t_idx)*basis_val + coord_3(:) = coord_3(:) + chi_3(chi_b_idx:chi_t_idx)*basis_val + end do + + ! Convert to radial coordinate + coord_radius(:) = sqrt(coord_1(:)**2 + coord_2(:)**2 + coord_3(:)**2) + + ! Increment height field at this DoF + height(h_b_idx:h_t_idx) = ( & + height(h_b_idx:h_t_idx) & + + rmultiplicity(h_b_idx:h_t_idx)*(coord_radius(:) - planet_radius) & + ) + end do + + ! -------------------------------------------------------------------------- ! + ! Native coordinates, or planar domain + ! -------------------------------------------------------------------------- ! + else + ! Third coordinate already gives the height above the Earth's mean radius + do df_h = 1, ndf_h + ! Indices for bottom and top height DoFs in this column + h_b_idx = map_h(df_h) + h_t_idx = map_h(df_h) + nlayers - 1 + + do df_chi = 1, ndf_chi + ! Indices for bottom and top chi DoFs in this column + chi_b_idx = map_chi(df_chi) + chi_t_idx = map_chi(df_chi) + nlayers - 1 + basis_val = basis_chi(1, df_chi, df_h) + + ! Increment height field at this DoF + height(h_b_idx:h_t_idx) = ( & + height(h_b_idx:h_t_idx) + rmultiplicity(h_b_idx:h_t_idx) & + * chi_3(chi_b_idx:chi_t_idx)*basis_val & + ) + end do + end do + end if + +end subroutine height_continuous_code + +end module sci_height_continuous_kernel_mod diff --git a/components/science/source/kernel/geometry/sci_height_discontinuous_kernel_mod.F90 b/components/science/source/kernel/geometry/sci_height_discontinuous_kernel_mod.F90 new file mode 100644 index 000000000..101e8aa4b --- /dev/null +++ b/components/science/source/kernel/geometry/sci_height_discontinuous_kernel_mod.F90 @@ -0,0 +1,173 @@ +!------------------------------------------------------------------------------- +! (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 Computes the height above the Earth's mean radius at the degrees of +!! freedom for a discontinuous function space. +!> @details Uses the chi coordinate fields to determine the height of nodes +!! above the Earth's mean radius, for discontinuous function spaces. + +module sci_height_discontinuous_kernel_mod + + use argument_mod, only: arg_type, func_type, & + GH_FIELD, GH_SCALAR, & + GH_REAL, GH_INTEGER, & + GH_READ, GH_WRITE, & + ANY_DISCONTINUOUS_SPACE_1, ANY_SPACE_9, & + CELL_COLUMN, GH_BASIS, GH_EVALUATOR + use base_mesh_config_mod, only: geometry_spherical + use constants_mod, only: r_def, i_def, l_def + use finite_element_config_mod, only: coord_system_xyz + use kernel_mod, only: kernel_type + + implicit none + private + + !----------------------------------------------------------------------------- + ! Public types + !----------------------------------------------------------------------------- + !> The type declaration for the kernel. Contains the metadata needed by the + !! Psy layer. + type, public, extends(kernel_type) :: height_discontinuous_kernel_type + private + type(arg_type) :: meta_args(5) = (/ & + arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1), & + arg_type(GH_FIELD*3, GH_REAL, GH_READ, ANY_SPACE_9), & + arg_type(GH_SCALAR, GH_INTEGER, GH_READ), & + arg_type(GH_SCALAR, GH_INTEGER, GH_READ), & + arg_type(GH_SCALAR, GH_REAL, GH_READ) & + /) + type(func_type) :: meta_funcs(1) = (/ & + func_type(ANY_SPACE_9, GH_BASIS) & + /) + integer :: operates_on = CELL_COLUMN + integer :: gh_shape = GH_EVALUATOR + contains + procedure, nopass :: height_discontinuous_code + end type + + !----------------------------------------------------------------------------- + ! Contained functions/subroutines + !----------------------------------------------------------------------------- + public :: height_discontinuous_code + +contains + +!> @brief Computes the height above the Earth's mean radius at the degrees of +!! freedom for a discontinuous function space. +!> @param[in] nlayers Number of layers in the mesh +!> @param[in,out] height The height field to compute +!> @param[in] chi_1 1st component of the coordinate fields +!> @param[in] chi_2 2nd component of the coordinate fields +!> @param[in] chi_3 3rd component of the coordinate fields +!> @param[in] geometry The geometry of the domain +!> @param[in] coord_system The coordinate system of the domain +!> @param[in] planet_radius The planet radius +!> @param[in] ndf_h Num DoFs per cell for height field +!> @param[in] undf_h Num DoFs in this partition for height field +!> @param[in] map_h DoF index map for height field +!> @param[in] ndf_chi Num DoFs per cell for chi fields +!> @param[in] undf_chi Num DoFs in this partition for chi fields +!> @param[in] map_chi DoF index map for chi fields +!> @param[in] basis_chi Chi basis functions evaluated at height DoFs +subroutine height_discontinuous_code( & + nlayers, & + height, & + chi_1, chi_2, chi_3, & + geometry, coord_system, planet_radius, & + ndf_h, undf_h, map_h, & + ndf_chi, undf_chi, map_chi, & + basis_chi & +) + implicit none + + ! Arguments + integer(kind=i_def), intent(in) :: nlayers + integer(kind=i_def), intent(in) :: ndf_h, undf_h + integer(kind=i_def), intent(in) :: ndf_chi, undf_chi + integer(kind=i_def), intent(in) :: geometry, coord_system + integer(kind=i_def), intent(in) :: map_h(ndf_h) + integer(kind=i_def), intent(in) :: map_chi(ndf_chi) + real(kind=r_def), intent(in) :: planet_radius + real(kind=r_def), intent(in) :: basis_chi(1,ndf_chi,ndf_h) + real(kind=r_def), intent(inout) :: height(undf_h) + real(kind=r_def), intent(in) :: chi_1(undf_chi) + real(kind=r_def), intent(in) :: chi_2(undf_chi) + real(kind=r_def), intent(in) :: chi_3(undf_chi) + + ! Internal variables + logical(kind=l_def) :: is_spherical_xyz + integer(kind=i_def) :: df_chi, df_h + integer(kind=i_def) :: h_b_idx, h_t_idx, chi_b_idx, chi_t_idx + real(kind=r_def) :: coord_1(nlayers), coord_2(nlayers), coord_3(nlayers) + real(kind=r_def) :: coord_radius(nlayers), basis_val + + ! Two cases: if the geometry is spherical with geocentric Cartesian coords, + ! all three chi components are needed to compute the height. Otherwise, only + ! the last chi component is needed, as this is already the vertical coordinate + + is_spherical_xyz = ( & + geometry == geometry_spherical .and. coord_system == coord_system_xyz & + ) + + ! -------------------------------------------------------------------------- ! + ! Cartesian system and spherical + ! -------------------------------------------------------------------------- ! + if (is_spherical_xyz) then + + do df_h = 1, ndf_h + ! Indices for bottom and top height DoFs in this column + h_b_idx = map_h(df_h) + h_t_idx = map_h(df_h) + nlayers - 1 + + ! Determine coordinates at this DoF using chi basis functions + coord_1(:) = 0.0_r_def + coord_2(:) = 0.0_r_def + coord_3(:) = 0.0_r_def + do df_chi = 1, ndf_chi + ! Indices for bottom and top chi DoFs in this column + chi_b_idx = map_chi(df_chi) + chi_t_idx = map_chi(df_chi) + nlayers - 1 + basis_val = basis_chi(1, df_chi, df_h) + + coord_1(:) = coord_1(:) + chi_1(chi_b_idx:chi_t_idx)*basis_val + coord_2(:) = coord_2(:) + chi_2(chi_b_idx:chi_t_idx)*basis_val + coord_3(:) = coord_3(:) + chi_3(chi_b_idx:chi_t_idx)*basis_val + end do + + ! Convert to radial coordinate + coord_radius(:) = sqrt(coord_1(:)**2 + coord_2(:)**2 + coord_3(:)**2) + + ! Increment height field at this DoF + height(h_b_idx:h_t_idx) = coord_radius(:) - planet_radius + end do + + ! -------------------------------------------------------------------------- ! + ! Native coordinates, or planar domain + ! -------------------------------------------------------------------------- ! + else + ! Third coordinate already gives the height above the Earth's mean radius + do df_h = 1, ndf_h + ! Indices for bottom and top height DoFs in this column + h_b_idx = map_h(df_h) + h_t_idx = map_h(df_h) + nlayers - 1 + height(h_b_idx:h_t_idx) = 0.0_r_def + + do df_chi = 1, ndf_chi + ! Indices for bottom and top chi DoFs in this column + chi_b_idx = map_chi(df_chi) + chi_t_idx = map_chi(df_chi) + nlayers - 1 + basis_val = basis_chi(1, df_chi, df_h) + + ! Increment height field at this DoF + height(h_b_idx:h_t_idx) = ( & + height(h_b_idx:h_t_idx) + chi_3(chi_b_idx:chi_t_idx)*basis_val & + ) + end do + end do + end if + +end subroutine height_discontinuous_code + +end module sci_height_discontinuous_kernel_mod diff --git a/components/science/unit-test/kernel/geometry/height_continuous_kernel_mod_test.pf b/components/science/unit-test/kernel/geometry/height_continuous_kernel_mod_test.pf new file mode 100644 index 000000000..b3e3d31e6 --- /dev/null +++ b/components/science/unit-test/kernel/geometry/height_continuous_kernel_mod_test.pf @@ -0,0 +1,267 @@ +!----------------------------------------------------------------------------- +! (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. +!----------------------------------------------------------------------------- +!> Test the height_continuous kernel +!! This test is parametrised to test both planar and spherical geometries. +module height_continuous_kernel_mod_test + + use base_mesh_config_mod, only : geometry_planar, & + geometry_spherical + use constants_mod, only : i_def, r_def, str_long + use funit + use get_unit_test_w2nodal_basis_mod, only : get_w0_w2nodal_basis + use get_unit_test_m3x3_q3x3x3_sizes_mod, only : get_w0_m3x3_q3x3x3_size, & + get_w2_m3x3_q3x3x3_size + use get_unit_test_m3x3_dofmap_mod, only : get_w0_m3x3_dofmap, & + get_w2_m3x3_dofmap + + implicit none + + private + public :: test_all, get_height_parameters, test_height_constructor, & + height_continuous_kernel_test_type + + @testParameter + type, public, extends(AbstractTestParameter) :: height_parameters_type + integer(kind=i_def) :: geometry + contains + procedure :: toString + end type height_parameters_type + + @TestCase(testParameters={get_height_parameters()}, constructor=test_height_constructor) + type, extends(ParameterizedTestCase) :: height_continuous_kernel_test_type + private + integer(kind=i_def) :: geometry + contains + procedure test_all + end type height_continuous_kernel_test_type + +contains + + function test_height_constructor(test_parameter) result (new_test) + + implicit none + + type(height_parameters_type), intent(in) :: test_parameter + type(height_continuous_kernel_test_type) :: new_test + + new_test%geometry = test_parameter%geometry + + end function test_height_constructor + + !> @brief A routine to convert test parameters to a string + function toString(this) result(output_string) + + implicit none + + class( height_parameters_type ), intent(in) :: this + character(:), allocatable :: output_string + character(str_long) :: geometry_string + + select case (this%geometry) + case (geometry_spherical ) + write(geometry_string, '(A)') 'spherical' + case (geometry_planar) + write(geometry_string, '(A)') 'planar' + end select + + output_string = trim(geometry_string) + + end function toString + + !> @brief Sets the parameters to be tested + function get_height_parameters() result (height_parameters) + + implicit none + + type(height_parameters_type) :: height_parameters(2) + + height_parameters = [ & + height_parameters_type(geometry_planar), & + height_parameters_type(geometry_spherical) & + ] + + end function get_height_parameters + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + @Test + subroutine test_all(this) + + use sci_height_continuous_kernel_mod, only : height_continuous_code + use finite_element_config_mod, only : coord_system_xyz + + implicit none + + class(height_continuous_kernel_test_type), intent(inout) :: this + + real(kind=r_def), parameter :: tol = 1.0e-3_r_def + + integer(kind=i_def) :: i, j, k, idx + integer(kind=i_def) :: nlayers, cell + integer(kind=i_def) :: ndf_w0, undf_w0 + integer(kind=i_def) :: ndf_w2, undf_w2 + integer(kind=i_def) :: ncells, dim_space, dim_space_diff + integer(kind=i_def) :: nqp_h, nqp_v + real(kind=r_def) :: dx, dy, dz, radius + real(kind=r_def) :: lon, lat, dlon, dlat, lon0, lat0, r + + integer(kind=i_def), allocatable :: map_w0(:,:) + integer(kind=i_def), allocatable :: map_w2(:,:) + real(kind=r_def), allocatable :: basis_w0_on_w2(:,:,:) + real(kind=r_def), allocatable :: chi_1(:), chi_2(:), chi_3(:) + real(kind=r_def), allocatable :: rmultiplicity(:) + real(kind=r_def), allocatable :: height_w2(:) + real(kind=r_def), allocatable :: height_w2_answer(:) + + ! ------------------------------------------------------------------------ ! + ! Set up domain + ! ------------------------------------------------------------------------ ! + ! Set number of layers + nlayers = 3 + + ! Set up degrees of freedom information + call get_w0_m3x3_q3x3x3_size( & + ndf_w0, undf_w0, ncells, dim_space, dim_space_diff, & + nqp_h, nqp_v, nlayers & + ) + call get_w2_m3x3_q3x3x3_size( & + ndf_w2, undf_w2, ncells, dim_space, dim_space_diff, & + nqp_h, nqp_v, nlayers & + ) + + ! Set up dof maps + call get_w0_m3x3_dofmap(map_w0) + call get_w2_m3x3_dofmap(map_w2) + + ! chi basis + call get_w0_w2nodal_basis(basis_w0_on_w2) + + ! Allocate height arrays + allocate(height_w2(undf_w2)) + allocate(height_w2_answer(undf_w2)) + + ! ------------------------------------------------------------------------ ! + ! Set up coordinates + ! ------------------------------------------------------------------------ ! + + ! Set up chi coefficients with quadratic extrusion + allocate(chi_1(undf_w0)) + allocate(chi_2(undf_w0)) + allocate(chi_3(undf_w0)) + + dz = 2.0_r_def + + if (this%geometry == geometry_planar) then + radius = 0.0_r_def + dx = 4.0_r_def + dy = 3.0_r_def + + ! Final coordinate field is height + cell = 1 + do j = 1, 3 + do i = 1, 3 + do k = 0, nlayers + chi_1(map_w0(1,cell)+k) = real(i-1)*dx + chi_2(map_w0(1,cell)+k) = real(j-1)*dy + chi_3(map_w0(1,cell)+k) = (real(k)**2)*dz + end do + cell = cell + 1 + end do + end do + + else + radius = 6.0_r_def + lon0 = 0.05_r_def + lat0 = 0.04_r_def + dlon = 0.01_r_def + dlat = 0.01_r_def + + ! Convert from spherical coordinates + cell = 1 + do j = 1, 3 + lat = lat0 + real(j-1)*dlat + do i = 1, 3 + lon = lon0 + real(i-1)*dlon + do k = 0, nlayers + r = radius + (real(k)**2)*dz + ! Spherical polar coordinates + chi_1(map_w0(1,cell)+k) = r*cos(lat)*cos(lon) + chi_2(map_w0(1,cell)+k) = r*cos(lat)*sin(lon) + chi_3(map_w0(1,cell)+k) = r*sin(lat) + end do + cell = cell + 1 + end do + end do + end if + + ! ------------------------------------------------------------------------ ! + ! Set rmultiplicity + ! ------------------------------------------------------------------------ ! + + allocate(rmultiplicity(undf_w2)) + rmultiplicity(:) = 0.5_r_def ! default value + ! Top and bottom DoFs need setting to 1 + do cell = 1, 9 + rmultiplicity(map_w2(5,cell)) = 1.0_r_def + rmultiplicity(map_w2(5,cell)+nlayers) = 1.0_r_def + end do + + ! ------------------------------------------------------------------------ ! + ! Set answer + ! ------------------------------------------------------------------------ ! + + ! Set answer for central cell: horizontal DoFs + cell = 5 + do i = 1, 4 + do k = 0, nlayers-1 + height_w2_answer(map_w2(i,cell)+k) = & + 0.5_r_def*dz*(real(k)**2 + real(k+1)**2) + end do + end do + + ! Set answer for central cell: vertical DoFs + cell = 5 + i = 5 + do k = 0, nlayers + height_w2_answer(map_w2(i,cell)+k) = (real(k)**2)*dz + end do + + ! ------------------------------------------------------------------------ ! + ! Test kernel + ! ------------------------------------------------------------------------ ! + + height_w2(:) = 0.0_r_def + + do cell = 1, ncells + call height_continuous_code( & + nlayers, & + height_w2, & + chi_1, chi_2, chi_3, & + rmultiplicity, & + this%geometry, coord_system_xyz, radius, & + ndf_w2, undf_w2, map_w2(:,cell), & + ndf_w0, undf_w0, map_w0(:,cell), & + basis_w0_on_w2(:,:,:) & + ) + end do + + cell = 5 + do i = 1, ndf_w2 + do k = 0, nlayers-1 + idx = map_w2(i,cell)+k + @assertEqual(height_w2_answer(idx), height_w2(idx), tol) + end do + end do + + deallocate(chi_1, chi_2, chi_3) + deallocate(rmultiplicity) + deallocate(height_w2, height_w2_answer) + deallocate(map_w0) + deallocate(map_w2) + deallocate(basis_w0_on_w2) + + end subroutine test_all + +end module height_continuous_kernel_mod_test diff --git a/components/science/unit-test/kernel/geometry/get_height_kernel_mod_test.pf b/components/science/unit-test/kernel/geometry/height_discontinuous_kernel_mod_test.pf similarity index 73% rename from components/science/unit-test/kernel/geometry/get_height_kernel_mod_test.pf rename to components/science/unit-test/kernel/geometry/height_discontinuous_kernel_mod_test.pf index 066d58f8d..400f083eb 100644 --- a/components/science/unit-test/kernel/geometry/get_height_kernel_mod_test.pf +++ b/components/science/unit-test/kernel/geometry/height_discontinuous_kernel_mod_test.pf @@ -7,7 +7,7 @@ !> Test the get height kernel. Sets up a planar domain with quadratic height !> extrusion, and tests answer in W3 and Wtheta spaces !> -module get_height_kernel_mod_test +module height_discontinuous_kernel_mod_test use constants_mod, only : i_def, r_def use funit @@ -26,69 +26,25 @@ module get_height_kernel_mod_test public :: test_all @TestCase - type, extends(TestCase), public :: get_height_kernel_test_type + type, extends(TestCase), public :: height_discontinuous_kernel_test_type private contains - procedure setUp - procedure tearDown procedure test_all - end type get_height_kernel_test_type + end type height_discontinuous_kernel_test_type contains - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - subroutine setUp( this ) - - use base_mesh_config_mod, only : geometry_planar, & - topology_fully_periodic - use finite_element_config_mod, only : cellshape_quadrilateral, & - coord_system_xyz - use feign_config_mod, only : feign_finite_element_config, & - feign_base_mesh_config - - implicit none - - class(get_height_kernel_test_type), intent(inout) :: this - - call feign_base_mesh_config( file_prefix='foo', & - prime_mesh_name='unit_test', & - geometry=geometry_planar, & - prepartitioned=.false., & - topology=topology_fully_periodic, & - fplane=.false., f_lat_deg=45.0_r_def ) - - call feign_finite_element_config( & - cellshape=cellshape_quadrilateral, & - coord_order=0_i_def, & - coord_system=coord_system_xyz, & - element_order_h=0_i_def, & - element_order_v=0_i_def, & - rehabilitate=.true. ) - - end subroutine setUp - - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - subroutine tearDown( this ) - - use configuration_mod, only: final_configuration - - implicit none - - class(get_height_kernel_test_type), intent(inout) :: this - - call final_configuration() - - end subroutine tearDown - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @Test subroutine test_all( this ) - use sci_get_height_kernel_mod, only : get_height_code + use base_mesh_config_mod, only : geometry_planar + use finite_element_config_mod, only : coord_system_xyz + use sci_height_discontinuous_kernel_mod, only : height_discontinuous_code implicit none - class(get_height_kernel_test_type), intent(inout) :: this + class(height_discontinuous_kernel_test_type), intent(inout) :: this real(r_def), parameter :: dx = 4.0_r_def, & dy = 3.0_r_def, & @@ -172,17 +128,19 @@ contains height_w3(:) = 0.0_r_def height_wtheta(:) = 0.0_r_def - call get_height_code(nlayers, & + call height_discontinuous_code(nlayers, & height_w3, & chi_1, chi_2, chi_3, & + geometry_planar, coord_system_xyz, & radius, & ndf_w3, undf_w3, map_w3(:,cell), & ndf_w0, undf_w0, map_w0(:,cell), & basis_w0_on_w3(:,:,:) ) - call get_height_code(nlayers, & + call height_discontinuous_code(nlayers, & height_wtheta, & chi_1, chi_2, chi_3, & + geometry_planar, coord_system_xyz, & radius, & ndf_wtheta, undf_wtheta, map_wtheta(:,cell), & ndf_w0, undf_w0, map_w0(:,cell), & @@ -216,4 +174,4 @@ contains end subroutine test_all -end module get_height_kernel_mod_test +end module height_discontinuous_kernel_mod_test