From aeb0739446353fca0f67f4ae79c46b587b50be0e Mon Sep 17 00:00:00 2001 From: thomasmelvin Date: Fri, 23 Jan 2026 16:28:46 +0000 Subject: [PATCH 1/2] Addition of code from SRS branch --- .../algorithm/sci_mapping_constants_mod.x90 | 65 ++++++++++++++++ .../sci_restrict_w2h_kernel_mod.F90 | 76 +++++++++---------- .../solver/sci_iterative_solver_mod.f90 | 40 +++++----- 3 files changed, 123 insertions(+), 58 deletions(-) diff --git a/components/science/source/algorithm/sci_mapping_constants_mod.x90 b/components/science/source/algorithm/sci_mapping_constants_mod.x90 index 7bcecf0f3..6512631b1 100644 --- a/components/science/source/algorithm/sci_mapping_constants_mod.x90 +++ b/components/science/source/algorithm/sci_mapping_constants_mod.x90 @@ -70,6 +70,7 @@ module sci_mapping_constants_mod ! Constants for mapping between meshes with different horizontal resolution type(inventory_by_mesh_type) :: intermesh_wghts_w2_inventory + type(inventory_by_mesh_type) :: intermesh_wghts_w2h_inventory type(inventory_by_mesh_type) :: intermesh_wghts_rdef_w3_inventory type(inventory_by_mesh_type) :: intermesh_wghts_rtran_w3_inventory @@ -107,6 +108,7 @@ module sci_mapping_constants_mod public :: get_intermesh_weights_w3_rdef public :: get_intermesh_weights_w3_rtran public :: get_intermesh_weights_w2 + public :: get_intermesh_weights_w2h public :: get_project_xdot_to_w2 public :: get_project_ydot_to_w2 public :: get_project_zdot_to_w2 @@ -631,6 +633,69 @@ contains end function get_intermesh_weights_w2 + !> @brief Returns a pointer to the weights for W2h prolongation + !> @param[in] fine_mesh The fine mesh for the transform + !> @param[in] coarse_mesh The coarse mesh for the transform + !> @return The field containing weights for W2h prolongation + function get_intermesh_weights_w2h(fine_mesh, coarse_mesh) result(weights) + + use sci_weights_prolong_w2h_kernel_mod, only: weights_prolong_w2h_kernel_type + + implicit none + + type(mesh_type), pointer, intent(in) :: coarse_mesh + type(mesh_type), pointer, intent(in) :: fine_mesh + integer(kind=i_def) :: intermesh_id + logical(kind=l_def) :: constant_exists + type(field_type), pointer :: weights + type(field_type) :: dummy_w2h_field + type(function_space_type), pointer :: fine_w2h_fs + type(function_space_type), pointer :: coarse_w2h_fs + + ! Check inventory is initialised + if (.not. intermesh_wghts_w2h_inventory%is_initialised()) then + call intermesh_wghts_w2h_inventory%initialise(name='intermesh_weights_w2h') + end if + + intermesh_id = intermesh_wghts_w2h_inventory%compute_intermesh_id( & + coarse_mesh, fine_mesh & + ) + constant_exists = intermesh_wghts_w2h_inventory%paired_object_exists( & + intermesh_id & + ) + + if (.not. constant_exists) then + + if (.not. coarse_mesh%query_mesh_map(fine_mesh)) then + write(log_scratch_space, '(A,I6,A,I6)') & + 'No mesh map exists for mapping between meshes ', & + coarse_mesh%get_id(), ' and ', fine_mesh%get_id() + call log_event(log_scratch_space, LOG_LEVEL_ERROR) + end if + + if ( subroutine_timers ) call timer('runtime_constants.mapping') + + coarse_w2h_fs => function_space_collection%get_fs(coarse_mesh, 0, 0, W2h) + fine_w2h_fs => function_space_collection%get_fs(fine_mesh, 0, 0, W2h) + + call dummy_w2h_field%initialise( coarse_w2h_fs ) + call invoke( setval_c(dummy_w2h_field, 0.0_r_def) ) + + call intermesh_wghts_w2h_inventory%add_field( & + weights, fine_w2h_fs, coarse_mesh, fine_mesh & + ) + + call invoke( setval_c(weights, 0.0_r_def), & + weights_prolong_w2h_kernel_type(weights, dummy_w2h_field) ) + + if ( subroutine_timers ) call timer('runtime_constants.mapping') + end if + + ! Get existing constant + call intermesh_wghts_w2h_inventory%get_field(coarse_mesh, fine_mesh, weights) + + end function get_intermesh_weights_w2h + !> @brief Returns a pointer to the weights for conservative W3 mapping !> @param[in] fine_mesh The fine mesh for the transform !> @param[in] coarse_mesh The coarse mesh for the transform diff --git a/components/science/source/kernel/inter_mesh/sci_restrict_w2h_kernel_mod.F90 b/components/science/source/kernel/inter_mesh/sci_restrict_w2h_kernel_mod.F90 index 32c727926..0d431e412 100644 --- a/components/science/source/kernel/inter_mesh/sci_restrict_w2h_kernel_mod.F90 +++ b/components/science/source/kernel/inter_mesh/sci_restrict_w2h_kernel_mod.F90 @@ -53,13 +53,13 @@ module sci_restrict_w2h_kernel_mod ! Contained functions/subroutines !------------------------------------------------------------------------------- -public :: restrict_w2h_code +public :: restrict_w2h_kernel_code ! Generic interface for real32 and real64 types - interface restrict_w2h_code + interface restrict_w2h_kernel_code module procedure & - restrict_w2h_code_r_single, & - restrict_w2h_code_r_double + restrict_w2h_kernel_code_r_single, & + restrict_w2h_kernel_code_r_double end interface contains @@ -93,22 +93,22 @@ module sci_restrict_w2h_kernel_mod ! R_SINGLE PRECISION ! ================== - subroutine restrict_w2h_code_r_single(nlayers, & - cell_map, & - ncell_fine_per_coarse_x, & - ncell_fine_per_coarse_y, & - ncell_fine, & - coarse_field, & - fine_field, & - face_selector_ew, & - face_selector_ns, & - undf_coarse, & - map_coarse, & - ndf, & - undf_fine, & - map_fine, & - undf_w3_2d, & - map_w3_2d ) + subroutine restrict_w2h_kernel_code_r_single(nlayers, & + cell_map, & + ncell_fine_per_coarse_x, & + ncell_fine_per_coarse_y, & + ncell_fine, & + coarse_field, & + fine_field, & + face_selector_ew, & + face_selector_ns, & + undf_coarse, & + map_coarse, & + ndf, & + undf_fine, & + map_fine, & + undf_w3_2d, & + map_w3_2d ) implicit none @@ -224,26 +224,26 @@ subroutine restrict_w2h_code_r_single(nlayers, & end do end do - end subroutine restrict_w2h_code_r_single + end subroutine restrict_w2h_kernel_code_r_single ! R_DOUBLE PRECISION ! ================== - subroutine restrict_w2h_code_r_double(nlayers, & - cell_map, & - ncell_fine_per_coarse_x, & - ncell_fine_per_coarse_y, & - ncell_fine, & - coarse_field, & - fine_field, & - face_selector_ew, & - face_selector_ns, & - undf_coarse, & - map_coarse, & - ndf, & - undf_fine, & - map_fine, & - undf_w3_2d, & - map_w3_2d ) + subroutine restrict_w2h_kernel_code_r_double(nlayers, & + cell_map, & + ncell_fine_per_coarse_x, & + ncell_fine_per_coarse_y, & + ncell_fine, & + coarse_field, & + fine_field, & + face_selector_ew, & + face_selector_ns, & + undf_coarse, & + map_coarse, & + ndf, & + undf_fine, & + map_fine, & + undf_w3_2d, & + map_w3_2d ) implicit none @@ -359,6 +359,6 @@ subroutine restrict_w2h_code_r_double(nlayers, & end do end do - end subroutine restrict_w2h_code_r_double + end subroutine restrict_w2h_kernel_code_r_double end module sci_restrict_w2h_kernel_mod diff --git a/components/science/source/solver/sci_iterative_solver_mod.f90 b/components/science/source/solver/sci_iterative_solver_mod.f90 index 48a1fde6a..e9f6a14ee 100644 --- a/components/science/source/solver/sci_iterative_solver_mod.f90 +++ b/components/science/source/solver/sci_iterative_solver_mod.f90 @@ -38,9 +38,9 @@ module sci_iterative_solver_mod type, public, abstract :: abstract_iterative_solver_type private !> Preconditioner - class(abstract_preconditioner_type), pointer :: prec => null() + class(abstract_preconditioner_type), allocatable :: prec !> Linear operator - class(abstract_linear_operator_type), pointer :: lin_op => null() + class(abstract_linear_operator_type), allocatable :: lin_op ! relative tolerance real(kind=r_def) :: r_tol @@ -475,8 +475,8 @@ module function cg_constructor(lin_op, prec, r_tol, a_tol, max_iter, & type(conjugate_gradient_type) :: self - self%lin_op => lin_op - self%prec => prec + allocate( self%lin_op, source = lin_op) + allocate( self%prec, source = prec) self%r_tol = r_tol self%a_tol = a_tol self%max_iter = max_iter @@ -627,8 +627,8 @@ module function bicgstab_constructor( lin_op, prec, r_tol, a_tol, max_iter, & write(log_scratch_space,'(A)') "bicgstab_constructor:" call log_event(log_scratch_space, LOG_LEVEL_INFO) - self%lin_op => lin_op - self%prec => prec + allocate(self%lin_op, source = lin_op ) + allocate(self%prec, source = prec ) self%r_tol = r_tol self%a_tol = a_tol self%max_iter = max_iter @@ -810,8 +810,8 @@ module function gmres_constructor( lin_op, prec, gcrk, r_tol, a_tol, max_iter, & logical(kind=l_def), intent(in) :: fail_on_non_converged type(gmres_type) :: self - self%lin_op => lin_op - self%prec => prec + allocate(self%lin_op, source = lin_op ) + allocate(self%prec, source = prec ) self%gcrk = gcrk self%r_tol = r_tol self%a_tol = a_tol @@ -1033,8 +1033,8 @@ module function fgmres_constructor( lin_op, prec, gcrk, r_tol, a_tol, max_iter, logical(kind=l_def), intent(in) :: fail_on_non_converged type(fgmres_type) :: self - self%lin_op => lin_op - self%prec => prec + allocate(self%lin_op, source = lin_op ) + allocate(self%prec, source = prec ) self%gcrk = gcrk self%r_tol = r_tol self%a_tol = a_tol @@ -1267,8 +1267,8 @@ module function gcr_constructor( lin_op, prec, gcrk, r_tol, a_tol, max_iter, & logical(kind=l_def), intent(in) :: fail_on_non_converged type(gcr_type) :: self - self%lin_op => lin_op - self%prec => prec + allocate(self%lin_op, source = lin_op ) + allocate(self%prec, source = prec ) self%gcrk = gcrk self%r_tol = r_tol self%a_tol = a_tol @@ -1438,8 +1438,8 @@ module function block_gcr_constructor( lin_op, prec, gcrk, r_tol, a_tol, max_ite logical(kind=l_def), intent(in) :: fail_on_non_converged type(block_gcr_type) :: self - self%lin_op => lin_op - self%prec => prec + allocate(self%lin_op, source = lin_op ) + allocate(self%prec, source = prec ) self%gcrk = gcrk self%r_tol = r_tol self%a_tol = a_tol @@ -1646,8 +1646,8 @@ module function precondition_only_constructor(lin_op, prec, monitor_convergence) logical(kind=l_def), intent(in) :: monitor_convergence type(precondition_only_type) :: self - self%lin_op => lin_op - self%prec => prec + allocate(self%lin_op, source = lin_op ) + allocate(self%prec, source = prec ) self%monitor_convergence = monitor_convergence if( .not. self%monitor_convergence ) then @@ -1726,8 +1726,8 @@ module function jacobi_constructor(lin_op, prec, r_tol, a_tol, max_iter, & real(kind=r_def), intent(in) :: rho_relax type(jacobi_type) :: self - self%lin_op => lin_op - self%prec => prec + allocate(self%lin_op, source = lin_op ) + allocate(self%prec, source = prec ) self%r_tol = r_tol self%a_tol = a_tol self%max_iter = max_iter @@ -1848,8 +1848,8 @@ module function chebyshev_constructor(lin_op, prec, r_tol, a_tol, max_iter, real(kind=r_def), intent(in) :: lmax type(chebyshev_type) :: self - self%lin_op => lin_op - self%prec => prec + allocate(self%lin_op, source = lin_op ) + allocate(self%prec, source = prec ) self%r_tol = r_tol self%a_tol = a_tol self%max_iter = max_iter From 33da686a4e934c7963e51ccb7b589dc320ad263a Mon Sep 17 00:00:00 2001 From: thomasmelvin Date: Mon, 26 Jan 2026 09:52:08 +0000 Subject: [PATCH 2/2] Add missing files --- .../inter_mesh/sci_prolong_w2h_kernel_mod.F90 | 390 ++++++++++++++++++ .../sci_weights_prolong_w2h_kernel_mod.F90 | 216 ++++++++++ 2 files changed, 606 insertions(+) create mode 100644 components/science/source/kernel/inter_mesh/sci_prolong_w2h_kernel_mod.F90 create mode 100644 components/science/source/kernel/inter_mesh/sci_weights_prolong_w2h_kernel_mod.F90 diff --git a/components/science/source/kernel/inter_mesh/sci_prolong_w2h_kernel_mod.F90 b/components/science/source/kernel/inter_mesh/sci_prolong_w2h_kernel_mod.F90 new file mode 100644 index 000000000..643a16056 --- /dev/null +++ b/components/science/source/kernel/inter_mesh/sci_prolong_w2h_kernel_mod.F90 @@ -0,0 +1,390 @@ +!----------------------------------------------------------------------------- +! (C) Crown copyright 2021 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 Perform the prolongation operation from a coarse grid W2h field to a +!! fine grid W2h field +!> @details Prolong the W2h coarse grid correction into all cells in a W2h field +!! on a fine grid. The fine grid cells are contained in that coarse +!! grid cell. Values are found by multiplying by weights. DoFs internal +!! to the coarse cell are treated differently to those DoFs on the +!! edges of the coarse cell: +!! - internal DoFs: values interpolated from coarse cell values +!! - edge DoFs: divide coarse cell flux by fine cell area fraction + +module sci_prolong_w2h_kernel_mod + +use argument_mod, only: arg_type, & + GH_FIELD, GH_REAL, & + GH_READ, GH_WRITE, & + GH_COARSE, GH_FINE, & + ANY_SPACE_2, CELL_COLUMN, & + ANY_DISCONTINUOUS_SPACE_2 +use constants_mod, only: i_def, r_def, r_single, r_double +use kernel_mod, only: kernel_type +use reference_element_mod, only: W, S, E, N + +implicit none + +private + +!------------------------------------------------------------------------------- +! Public types +!------------------------------------------------------------------------------- +!> The type declaration for the kernel. Contains the metadata needed by the +!> Psy layer. +!> + +type, public, extends(kernel_type) :: prolong_w2h_kernel_type + private + type(arg_type) :: meta_args(3) = (/ & + arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_2, & + mesh_arg=GH_FINE ), & + arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_SPACE_2, mesh_arg=GH_COARSE ), & + arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_2, & + mesh_arg=GH_FINE ) & + /) + integer :: operates_on = CELL_COLUMN +end type prolong_w2h_kernel_type + +!------------------------------------------------------------------------------- +! Contained functions/subroutines +!------------------------------------------------------------------------------- + +public :: prolong_w2h_kernel_code + + ! Generic interface for real32 and real64 types + interface prolong_w2h_kernel_code + module procedure & + prolong_w2h_kernel_code_r_single, & + prolong_w2h_kernel_code_r_double + end interface + +contains + + !> @brief Subroutine to perform the W2h prolongation operation + !> @param[in] nlayers Number of layers in a model column + !> @param[in] cell_map A 2D index map of which fine grid + !! cells lie in the coarse grid cell + !> @param[in] ncell_fine_per_coarse_x Number of fine cells per coarse + !! cell in the horizontal x-direction + !> @param[in] ncell_fine_per_coarse_x Number of fine cells per coarse + !! cell in the horizontal y-direction + !> @param[in] ncell_fine Number of cells in the partition + !! for the fine grid + !> @param[in,out] fine_field Fine grid W2h field to restrict + !> @param[in] coarse_field Coarse grid W2h field to compute + !> @param[in] weights A fine grid W2h field containing + !! weights for the mapping process + !> @param[in] ndf Num of DoFs per cell on both grids + !> @param[in] undf_fine Total num of DoFs on the fine grid + !! for this mesh partition + !> @param[in] map_fine DoFmap of cells on the fine grid + !> @param[in] undf_coarse Total num of DoFs on the coarse + !! grid for this mesh partition + !> @param[in] map_coarse DoFmap of cells on the coarse grid + + ! R_SINGLE PRECISION + ! ================== + subroutine prolong_w2h_kernel_code_r_single(nlayers, & + cell_map, & + ncell_fine_per_coarse_x, & + ncell_fine_per_coarse_y, & + ncell_fine, & + fine_field, & + coarse_field, & + weights, & + ndf, & + undf_fine, & + map_fine, & + undf_coarse, & + map_coarse ) + + implicit none + + integer(kind=i_def), intent(in) :: nlayers + integer(kind=i_def), intent(in) :: ncell_fine_per_coarse_x + integer(kind=i_def), intent(in) :: ncell_fine_per_coarse_y + integer(kind=i_def), intent(in) :: ncell_fine + integer(kind=i_def), intent(in) :: ndf + integer(kind=i_def), intent(in) :: undf_fine, undf_coarse + + ! Fields + real(kind=r_single), intent(inout) :: fine_field(undf_fine) + real(kind=r_single), intent(in) :: coarse_field(undf_coarse) + real(kind=r_def), intent(in) :: weights(undf_fine) + + ! Maps + integer(kind=i_def), intent(in) :: map_fine(ndf, ncell_fine) + integer(kind=i_def), intent(in) :: map_coarse(ndf) + integer(kind=i_def), intent(in) :: cell_map(ncell_fine_per_coarse_x, ncell_fine_per_coarse_y) + + ! Internal arguments + integer(kind=i_def) :: df, opp_df, k, x_idx, y_idx + real(kind=r_single) :: edge_weight, internal_weight + + !=========================================================================== + ! Horizontal components + !=========================================================================== + + do k = 0, nlayers-1 + + ! The rows and columns forming the cell map match the arrangment + ! of fine cells within the coarse cell + + ! These are aligned as follows with the LFRic directions: + ! N + ! |--------------| + ! | row 1 | + ! |c c| + ! |o o| + ! W |l l| E + ! | | + ! |1 nx| + ! | row ny | + ! |--------------| + ! S + + !------------------------------------------------------------------------- + ! Edges of coarse cell + !------------------------------------------------------------------------- + + ! The weights here describe the flux through the edge of the fine cell + ! This should be the fractional area of the coarse cell edge taken up + ! by the fine cell. If fine cells have equal area then this will be the + ! reciprocal of the number of fine cells along that edge + + do x_idx = 1, ncell_fine_per_coarse_x + ! N edge is first row of cell map + y_idx = 1 + df = N + fine_field(map_fine(df,cell_map(x_idx,y_idx))+k) = & + real(weights(map_fine(df,cell_map(x_idx,y_idx))+k), r_single) & + * coarse_field(map_coarse(df)+k) + + ! S edge is last row of cell map + y_idx = ncell_fine_per_coarse_y + df = S + fine_field(map_fine(df,cell_map(x_idx,y_idx))+k) = & + real(weights(map_fine(df,cell_map(x_idx,y_idx))+k), r_single) & + * coarse_field(map_coarse(df)+k) + end do + + do y_idx = 1, ncell_fine_per_coarse_y + ! W edge is first column of cell map + x_idx = 1 + df = W + fine_field(map_fine(df,cell_map(x_idx,y_idx))+k) = & + real(weights(map_fine(df,cell_map(x_idx,y_idx))+k), r_single) & + * coarse_field(map_coarse(df)+k) + + ! E edge is last column of cell map + x_idx = ncell_fine_per_coarse_x + df = E + fine_field(map_fine(df,cell_map(x_idx,y_idx))+k) = & + real(weights(map_fine(df,cell_map(x_idx,y_idx))+k), r_single) & + * coarse_field(map_coarse(df)+k) + end do + + !------------------------------------------------------------------------- + ! DoFs internal to coarse cell + !------------------------------------------------------------------------- + + ! For internal cells, we interpolate. The weights at these DoFs give the + ! relative weighting of the WEST and NORTH components. The contributions + ! from the east and south components are given by 1 minus those weights. + ! All must be combined with the area fraction weights to give the true + ! values. As the fine cells are internal to a coarse cell, they will all + ! have the same orientation as one another. + + ! Do West/East DoFs + ! These loops will only execute if there are internal DoFs + df = W + opp_df = E + do y_idx = 1, ncell_fine_per_coarse_y + edge_weight = real(weights(map_fine(df,cell_map(1,y_idx))+k), r_single) + + do x_idx = 2, ncell_fine_per_coarse_x + internal_weight = real(weights(map_fine(df,cell_map(x_idx,y_idx))+k), r_single) + + fine_field(map_fine(df,cell_map(x_idx,y_idx))+k) = & + edge_weight * ( internal_weight * coarse_field(map_coarse(df)+k) + & + (1.0_r_single - internal_weight) * coarse_field(map_coarse(opp_df)+k) ) + end do + end do + + ! Do North/South DoFs + ! These loops will only execute if there are internal DoFs + df = N + opp_df = S + do x_idx = 1, ncell_fine_per_coarse_x + edge_weight = real(weights(map_fine(df,cell_map(x_idx,1))+k), r_single) + + do y_idx = 2, ncell_fine_per_coarse_y + internal_weight = real(weights(map_fine(df,cell_map(x_idx,y_idx))+k), r_single) + + fine_field(map_fine(df,cell_map(x_idx,y_idx))+k) = & + edge_weight * ( internal_weight * coarse_field(map_coarse(df)+k) + & + (1.0_r_single - internal_weight) * coarse_field(map_coarse(opp_df)+k) ) + + end do + end do + + end do + + end subroutine prolong_w2h_kernel_code_r_single + + ! R_DOUBLE PRECISION + ! ================== + subroutine prolong_w2h_kernel_code_r_double(nlayers, & + cell_map, & + ncell_fine_per_coarse_x, & + ncell_fine_per_coarse_y, & + ncell_fine, & + fine_field, & + coarse_field, & + weights, & + ndf, & + undf_fine, & + map_fine, & + undf_coarse, & + map_coarse ) + + implicit none + + integer(kind=i_def), intent(in) :: nlayers + integer(kind=i_def), intent(in) :: ncell_fine_per_coarse_x + integer(kind=i_def), intent(in) :: ncell_fine_per_coarse_y + integer(kind=i_def), intent(in) :: ncell_fine + integer(kind=i_def), intent(in) :: ndf + integer(kind=i_def), intent(in) :: undf_fine, undf_coarse + + ! Fields + real(kind=r_double), intent(inout) :: fine_field(undf_fine) + real(kind=r_double), intent(in) :: coarse_field(undf_coarse) + real(kind=r_def), intent(in) :: weights(undf_fine) + + ! Maps + integer(kind=i_def), intent(in) :: map_fine(ndf, ncell_fine) + integer(kind=i_def), intent(in) :: map_coarse(ndf) + integer(kind=i_def), intent(in) :: cell_map(ncell_fine_per_coarse_x, ncell_fine_per_coarse_y) + + ! Internal arguments + integer(kind=i_def) :: df, opp_df, k, x_idx, y_idx + real(kind=r_double) :: edge_weight, internal_weight + + !=========================================================================== + ! Horizontal components + !=========================================================================== + + do k = 0, nlayers-1 + + ! The rows and columns forming the cell map match the arrangment + ! of fine cells within the coarse cell + + ! These are aligned as follows with the LFRic directions: + ! N + ! |--------------| + ! | row 1 | + ! |c c| + ! |o o| + ! W |l l| E + ! | | + ! |1 nx| + ! | row ny | + ! |--------------| + ! S + + !------------------------------------------------------------------------- + ! Edges of coarse cell + !------------------------------------------------------------------------- + + ! The weights here describe the flux through the edge of the fine cell + ! This should be the fractional area of the coarse cell edge taken up + ! by the fine cell. If fine cells have equal area then this will be the + ! reciprocal of the number of fine cells along that edge + + do x_idx = 1, ncell_fine_per_coarse_x + ! N edge is first row of cell map + y_idx = 1 + df = N + fine_field(map_fine(df,cell_map(x_idx,y_idx))+k) = & + real(weights(map_fine(df,cell_map(x_idx,y_idx))+k), r_double) & + * coarse_field(map_coarse(df)+k) + + ! S edge is last row of cell map + y_idx = ncell_fine_per_coarse_y + df = S + fine_field(map_fine(df,cell_map(x_idx,y_idx))+k) = & + real(weights(map_fine(df,cell_map(x_idx,y_idx))+k), r_double) & + * coarse_field(map_coarse(df)+k) + end do + + do y_idx = 1, ncell_fine_per_coarse_y + ! W edge is first column of cell map + x_idx = 1 + df = W + fine_field(map_fine(df,cell_map(x_idx,y_idx))+k) = & + real(weights(map_fine(df,cell_map(x_idx,y_idx))+k), r_double) & + * coarse_field(map_coarse(df)+k) + + ! E edge is last column of cell map + x_idx = ncell_fine_per_coarse_x + df = E + fine_field(map_fine(df,cell_map(x_idx,y_idx))+k) = & + real(weights(map_fine(df,cell_map(x_idx,y_idx))+k), r_double) & + * coarse_field(map_coarse(df)+k) + end do + + !------------------------------------------------------------------------- + ! DoFs internal to coarse cell + !------------------------------------------------------------------------- + + ! For internal cells, we interpolate. The weights at these DoFs give the + ! relative weighting of the WEST and NORTH components. The contributions + ! from the east and south components are given by 1 minus those weights. + ! All must be combined with the area fraction weights to give the true + ! values. As the fine cells are internal to a coarse cell, they will all + ! have the same orientation as one another. + + ! Do West/East DoFs + ! These loops will only execute if there are internal DoFs + df = W + opp_df = E + do y_idx = 1, ncell_fine_per_coarse_y + edge_weight = real(weights(map_fine(df,cell_map(1,y_idx))+k), r_double) + + do x_idx = 2, ncell_fine_per_coarse_x + internal_weight = real(weights(map_fine(df,cell_map(x_idx,y_idx))+k), r_double) + + fine_field(map_fine(df,cell_map(x_idx,y_idx))+k) = & + edge_weight * ( internal_weight * coarse_field(map_coarse(df)+k) + & + (1.0_r_double - internal_weight) * coarse_field(map_coarse(opp_df)+k) ) + end do + end do + + ! Do North/South DoFs + ! These loops will only execute if there are internal DoFs + df = N + opp_df = S + do x_idx = 1, ncell_fine_per_coarse_x + edge_weight = real(weights(map_fine(df,cell_map(x_idx,1))+k), r_double) + + do y_idx = 2, ncell_fine_per_coarse_y + internal_weight = real(weights(map_fine(df,cell_map(x_idx,y_idx))+k), r_double) + + fine_field(map_fine(df,cell_map(x_idx,y_idx))+k) = & + edge_weight * ( internal_weight * coarse_field(map_coarse(df)+k) + & + (1.0_r_double - internal_weight) * coarse_field(map_coarse(opp_df)+k) ) + + end do + end do + + end do + + end subroutine prolong_w2h_kernel_code_r_double + + +end module sci_prolong_w2h_kernel_mod diff --git a/components/science/source/kernel/inter_mesh/sci_weights_prolong_w2h_kernel_mod.F90 b/components/science/source/kernel/inter_mesh/sci_weights_prolong_w2h_kernel_mod.F90 new file mode 100644 index 000000000..eaf1f63ce --- /dev/null +++ b/components/science/source/kernel/inter_mesh/sci_weights_prolong_w2h_kernel_mod.F90 @@ -0,0 +1,216 @@ +!----------------------------------------------------------------------------- +! (C) Crown copyright 2021 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 weights for the prologation operation on a W2h field. +!> @details Computes weights for performing the prolongation operation to take a +!! W2h field from a coarse mesh to a fine mesh. The weights are stored +!! in a W2h field on the fine mesh. For horizontal W2h components, the +!! fine cell DoFs on the edges of coarse cells have weight given by +!! the reciprocal of the number of fine cells touching that edge of +!! the coarse cell. At fine cell DoFs on the interior of the coarse +!! cell, weights correspond to interpolation between the edges of the +!! coarse cell. +!! This is only designed to work with the lowest order W2h fields and +!! for fine cells that are exactly nested within a coarse cell. + +module sci_weights_prolong_w2h_kernel_mod + +use argument_mod, only: arg_type, & + GH_FIELD, GH_REAL, & + GH_READ, GH_WRITE, & + GH_COARSE, GH_FINE, & + ANY_DISCONTINUOUS_SPACE_1, & + ANY_DISCONTINUOUS_SPACE_2, & + CELL_COLUMN +use constants_mod, only: i_def, r_def +use kernel_mod, only: kernel_type +use reference_element_mod, only: W, S, E, N + +implicit none + +private + +!------------------------------------------------------------------------------- +! Public types +!------------------------------------------------------------------------------- +!> The type declaration for the kernel. Contains the metadata needed by the +!> Psy layer. +!> + +type, public, extends(kernel_type) :: weights_prolong_w2h_kernel_type + private + type(arg_type) :: meta_args(2) = (/ & + arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1, mesh_arg=GH_FINE ), & + arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_2, mesh_arg=GH_COARSE ) & + /) + integer :: operates_on = CELL_COLUMN +contains + procedure, nopass :: weights_prolong_w2h_kernel_code +end type weights_prolong_w2h_kernel_type + +!------------------------------------------------------------------------------- +! Contained functions/subroutines +!------------------------------------------------------------------------------- + +public :: weights_prolong_w2h_kernel_code + +contains + + !> @brief Computes weights for the prologation operation on a W2h field. + !> @param[in] nlayers Number of layers in a model column + !> @param[in] cell_map A 2D index map of which fine grid + !! cells lie in the coarse grid cell + !> @param[in] ncell_fine_per_coarse_x Number of fine cells per coarse + !! cell in the horizontal x-direction + !> @param[in] ncell_fine_per_coarse_x Number of fine cells per coarse + !! cell in the horizontal y-direction + !> @param[in] ncell_fine Number of cells in the partition + !! for the fine grid + !> @param[in,out] weights_fine The fine grid W2h field that contains + !! the weights to be computed here for + !! prolongation + !> @param[in] dummy_coarse An unused coarse W2h field supplied + !! to ensure that appropriate coarse + !! field properties are provided. + !> @param[in] ndf Num of DoFs per cell on both grids + !> @param[in] undf_fine Total num of DoFs on the fine grid + !! for this mesh partition + !> @param[in] map_fine DoFmap of cells on the fine grid + !> @param[in] undf_coarse Total num of DoFs on the coarse + !! grid for this mesh partition + !> @param[in] map_coarse DoFmap of cells on the coarse grid + subroutine weights_prolong_w2h_kernel_code(nlayers, & + cell_map, & + ncell_fine_per_coarse_x, & + ncell_fine_per_coarse_y, & + ncell_fine, & + weights_fine, & + dummy_coarse, & + ndf, & + undf_fine, & + map_fine, & + undf_coarse, & + map_coarse ) + + implicit none + + integer(kind=i_def), intent(in) :: nlayers + integer(kind=i_def), intent(in) :: ncell_fine_per_coarse_x + integer(kind=i_def), intent(in) :: ncell_fine_per_coarse_y + integer(kind=i_def), intent(in) :: ncell_fine + integer(kind=i_def), intent(in) :: ndf + integer(kind=i_def), intent(in) :: undf_fine, undf_coarse + + ! Fields + real(kind=r_def), intent(inout) :: weights_fine(undf_fine) + real(kind=r_def), intent(in) :: dummy_coarse(undf_coarse) + + ! Maps + integer(kind=i_def), intent(in) :: map_fine(ndf, ncell_fine) + integer(kind=i_def), intent(in) :: map_coarse(ndf) + integer(kind=i_def), intent(in) :: cell_map(ncell_fine_per_coarse_x, ncell_fine_per_coarse_y) + + ! Internal arguments + integer(kind=i_def) :: df, k, x_idx, y_idx + + !=========================================================================== + ! Horizontal components + !=========================================================================== + + do k = 0, nlayers-1 + + ! The rows and columns forming the cell map match the arrangment + ! of fine cells within the coarse cell + + ! These are aligned as follows with the LFRic directions: + ! N + ! |--------------| + ! | row 1 | + ! |c c| + ! |o o| + ! W |l l| E + ! | | + ! |1 nx| + ! | row ny | + ! |--------------| + ! S + + !------------------------------------------------------------------------- + ! Edges of coarse cell + !------------------------------------------------------------------------- + + ! The weights here describe the flux through the edge of the fine cell + ! This should be the fractional area of the coarse cell edge taken up + ! by the fine cell. If fine cells have equal area then this will be the + ! reciprocal of the number of fine cells along that edge + + do x_idx = 1, ncell_fine_per_coarse_x + ! N edge is first row of cell map + y_idx = 1 + df = N + weights_fine(map_fine(df,cell_map(x_idx,y_idx))+k) = & + 1.0_r_def / real(ncell_fine_per_coarse_x, r_def) + + ! S edge is last row of cell map + y_idx = ncell_fine_per_coarse_y + df = S + weights_fine(map_fine(df,cell_map(x_idx,y_idx))+k) = & + 1.0_r_def / real(ncell_fine_per_coarse_x, r_def) + end do + + do y_idx = 1, ncell_fine_per_coarse_y + ! W edge is first column of cell map + x_idx = 1 + df = W + weights_fine(map_fine(df,cell_map(x_idx,y_idx))+k) = & + 1.0_r_def / real(ncell_fine_per_coarse_y, r_def) + + ! E edge is last column of cell map + x_idx = ncell_fine_per_coarse_x + df = E + weights_fine(map_fine(df,cell_map(x_idx,y_idx))+k) = & + 1.0_r_def / real(ncell_fine_per_coarse_y, r_def) + end do + + !------------------------------------------------------------------------- + ! DoFs internal to coarse cell + !------------------------------------------------------------------------- + + ! For internal cells, we interpolate. The weights at these DoFs give the + ! relative weighting of the WEST and NORTH components. The contributions + ! from the east and south components are given by 1 minus those weights. + ! All must be combined with the area fraction weights to give the true + ! values. As the fine cells are internal to a coarse cell, they will all + ! have the same orientation as one another. + + ! Do West/East DoFs + ! These loops will only execute if there are internal DoFs + df = W + do y_idx = 1, ncell_fine_per_coarse_y + do x_idx = 2, ncell_fine_per_coarse_x + + weights_fine(map_fine(df,cell_map(x_idx,y_idx))+k) = & + real(ncell_fine_per_coarse_x + 1 - x_idx, r_def) / real(ncell_fine_per_coarse_x, r_def) + end do + end do + + ! Do North/South DoFs + ! These loops will only execute if there are internal DoFs + df = N + do x_idx = 1, ncell_fine_per_coarse_x + do y_idx = 2, ncell_fine_per_coarse_y + + weights_fine(map_fine(df,cell_map(x_idx,y_idx))+k) = & + real(ncell_fine_per_coarse_y + 1 - y_idx, r_def) / real(ncell_fine_per_coarse_y, r_def) + + end do + end do + + end do + + end subroutine weights_prolong_w2h_kernel_code + +end module sci_weights_prolong_w2h_kernel_mod