Merge remote-tracking branch 'origin/master' into 2889_extract_non_al…

…loc_arrays
stfc · Feb 14, 2025 · a00b7aa · a00b7aa
2 parents cd7c3d3 + 6afc06a
commit a00b7aa
Show file tree

Hide file tree

Showing 55 changed files with 6,221 additions and 4,526 deletions.
diff --git a/changelog b/changelog
@@ -55,6 +55,12 @@
 
 	19) PR #2885. Update integration test dependencies.
 
+	20) PR #2799 for #2779. Implement split finite-element order for
+	function spaces in the LFRic API and update Adjoint tests to use it.
+
+	21) PR #2883 for #2875. Fix issues cannonicalising WHEREs with elemental
+	functions.
+
 release 3.0.0 6th of December 2024
 
 	1) PR #2477 for #2463. Add support for Fortran Namelist statements.

diff --git a/doc/user_guide/psyclone_kern.rst b/doc/user_guide/psyclone_kern.rst
@@ -32,7 +32,7 @@
 .. POSSIBILITY OF SUCH DAMAGE.
 .. -----------------------------------------------------------------------------
 .. Written by R. W. Ford and A. R. Porter, STFC Daresbury Lab
-.. Modified by I. Kavcic and L. Turner, Met Office
+.. Modified by I. Kavcic, L. Turner and J. Dendy, Met Office
 
 PSyclone Kernel Tools
 =====================
@@ -558,14 +558,15 @@ gives the following algorithm layer code:
       use mesh_mod, only : mesh_type
       use simple_mod, only : simple_type
       use constants_mod, only : i_def, r_def
-      integer(kind=i_def), parameter :: element_order = 1_i_def
+      integer(kind=i_def), parameter :: element_order_h = 1_i_def
+      integer(kind=i_def), parameter :: element_order_v = 1_i_def
       type(mesh_type), pointer, intent(in) :: mesh
       type(field_type), dimension(3), intent(in), optional :: chi
       type(field_type), intent(in), optional :: panel_id
       TYPE(function_space_type), POINTER :: vector_space_w1_ptr
       type(field_type) :: field_1
 
-      vector_space_w1_ptr => function_space_collection % get_fs(mesh, element_order, w1)
+      vector_space_w1_ptr => function_space_collection % get_fs(mesh, element_order_h, element_order_v, w1)
       call field_1 % initialise(vector_space=vector_space_w1_ptr, name='field_1')
       call invoke(setval_c(field_1, 1.0_r_def), simple_type(field_1))
 

diff --git a/doc/user_guide/transformations.rst b/doc/user_guide/transformations.rst
@@ -34,6 +34,7 @@
 .. Written by: R. W. Ford, A. R. Porter and S. Siso, STFC Daresbury Lab
 ..             A. B. G. Chalk and N. Nobre, STFC Daresbury Lab
 ..             I. Kavcic, Met Office
+..             J. Dendy, Met Office
 
 .. _transformations:
 
@@ -216,7 +217,7 @@ process, optional parameters for the transformation are also
 provided this way. A simple example::
 
     kctrans = Dynamo0p3KernelConstTrans()
-    kctrans.apply(kernel, {"element_order": 0, "quadrature": True})
+    kctrans.apply(kernel, {"element_order_h": 0, "element_order_v": 0, "quadrature": True})
 
 The same ``options`` dictionary will be used when calling ``validate``.
 
@@ -287,19 +288,19 @@ can be found in the API-specific sections).
 .. autoclass:: psyclone.psyir.transformations.AllArrayAccess2LoopTrans
     :members: apply
     :noindex:
-  
+
 ####
 
 .. autoclass:: psyclone.psyir.transformations.ArrayAccess2LoopTrans
     :members: apply
     :noindex:
-  
+
 ####
 
 .. autoclass:: psyclone.psyir.transformations.ArrayAssignment2LoopsTrans
     :members: apply
     :noindex:
-  
+
 ####
 
 .. autoclass:: psyclone.psyir.transformations.ChunkLoopTrans
@@ -802,7 +803,7 @@ PSyclone supports OpenMP Tasking, through the `OMPTaskloopTrans` and
 transformations can be applied to loops, whilst the `OMPTaskwaitTrans`
 operator is applied to an OpenMP Parallel Region, and computes the dependencies
 caused by Taskloops, and adds OpenMP Taskwait statements to satisfy those
-dependencies. An example of using OpenMP tasking is available in 
+dependencies. An example of using OpenMP tasking is available in
 `PSyclone/examples/nemo/eg1/openmp_taskloop_trans.py`.
 
 OpenCL

diff --git a/examples/lfric/code/gw_mixed_schur_preconditioner_alg_mod.x90 b/examples/lfric/code/gw_mixed_schur_preconditioner_alg_mod.x90
@@ -38,6 +38,7 @@
 ! -----------------------------------------------------------------------------
 ! Modified by I. Kavcic, Met Office
 ! Modified by R. W. Ford, STFC Daresbury Lab
+! Modified by J. Dendy, Met Office
 
 !>@brief Preconditioner for the gravity-wave system
 module gw_mixed_schur_preconditioner_alg_mod
@@ -125,7 +126,6 @@ contains
 						     result(self)
 
     use function_space_mod,            only: function_space_type
-    use finite_element_config_mod,     only: element_order
     use quadrature_xyoz_mod,           only: quadrature_xyoz_type
     use quadrature_rule_gaussian_mod,  only: quadrature_rule_gaussian_type
     use matrix_vector_kernel_mod,      only: matrix_vector_kernel_type

diff --git a/examples/lfric/eg13/kernel_constants.py b/examples/lfric/eg13/kernel_constants.py
@@ -32,6 +32,7 @@
 # POSSIBILITY OF SUCH DAMAGE.
 # -----------------------------------------------------------------------------
 # Authors: R. W. Ford, N. Nobre and S. Siso, STFC Daresbury Lab
+# Modified by: J. Dendy, Met Office
 
 
 '''An example PSyclone transformation script which makes ndofs, nqp*
@@ -62,10 +63,11 @@
 # associated kernel value constant (rather than passing it in by
 # argument).
 NUMBER_OF_LAYERS = 20
-# The element order to use when modifying a kernel to make the
+# The element orders to use when modifying a kernel to make the
 # associated degrees of freedom values constant (rather than passing
 # them in by argument).
-ELEMENT_ORDER = 0
+ELEMENT_ORDER_H = 0
+ELEMENT_ORDER_V = 0
 # Whether or not to make the number of quadrature points constant in a
 # kernel (rather than passing them in by argument).
 CONSTANT_QUADRATURE = True
@@ -86,7 +88,8 @@ def trans(psyir):
         try:
             const_trans.apply(kernel,
                               {"number_of_layers": NUMBER_OF_LAYERS,
-                               "element_order": ELEMENT_ORDER,
+                               "element_order_h": ELEMENT_ORDER_H,
+                               "element_order_v": ELEMENT_ORDER_V,
                                "quadrature": CONSTANT_QUADRATURE})
         except TransformationError:
             print(f"    Failed to modify kernel '{kernel.name}'")
diff --git a/examples/lfric/eg14/main.x90 b/examples/lfric/eg14/main.x90
@@ -30,6 +30,7 @@
 ! -----------------------------------------------------------------------------
 ! Author: J. Henrichs, Bureau of Meteorology
 ! Modifications: A. R. Porter, STFC Daresbury Laboratory
+!                J. Dendy, Met Office
 
 program main
 
@@ -57,7 +58,8 @@ program main
     procedure (partitioner_interface), pointer :: partitioner_ptr
     type(field_type)                           :: field1, field2
     integer(kind=i_def)                        :: lfric_fs = W0   ! W0
-    integer(kind=i_def)                        :: element_order = 1
+    integer(kind=i_def)                        :: element_order_h = 1
+    integer(kind=i_def)                        :: element_order_v = 1
     integer(kind=i_def)                        :: ndata_sz, istp
     real(kind=r_def) :: chksm
 
@@ -84,9 +86,10 @@ program main
     mesh = mesh_type(global_mesh_ptr, partition, extrusion_ptr)
     write (*,*) "Mesh has", mesh%get_nlayers(), "layers."
     ndata_sz = 1
-    vector_space = function_space_type( mesh,          &
-                                        element_order, &
-                                        lfric_fs,      &
+    vector_space = function_space_type( mesh,            &
+                                        element_order_h, &
+                                        element_order_v, &
+                                        lfric_fs,        &
                                         ndata_sz)
     vector_space_ptr => vector_space
     call field1%initialise( vector_space = vector_space_ptr, name="field1" )

diff --git a/examples/lfric/eg17/full_example/README.md b/examples/lfric/eg17/full_example/README.md
@@ -9,7 +9,7 @@ The following steps are required for this (using simplified code examples):
     ```fortran
     global_mesh = global_mesh_type()
     ```
-   
+
 2) A 1x1 planar partition for one process is created:
     ```fortran
     partitioner_ptr => partitioner_planar
@@ -34,9 +34,10 @@ The following steps are required for this (using simplified code examples):
 
 5) Create a function/vector space:
     ```fortran
-    vector_space = function_space_type( mesh,          &
-                                        element_order, &
-                                        lfric_fs,      &
+    vector_space = function_space_type( mesh,            &
+                                        element_order_h, &
+                                        element_order_v, &
+                                        lfric_fs,        &
                                         ndata_sz)
     ```
 

diff --git a/examples/lfric/eg17/full_example/main.x90 b/examples/lfric/eg17/full_example/main.x90
@@ -30,6 +30,7 @@
 ! -----------------------------------------------------------------------------
 ! Author: J. Henrichs, Bureau of Meteorology
 ! Modifications: A. R. Porter, STFC Daresbury Laboratory
+!                J. Dendy, Met Office
 
 program main
 
@@ -57,7 +58,8 @@ program main
     procedure (partitioner_interface), pointer :: partitioner_ptr
     type(field_type)                           :: field1, field2
     integer(kind=i_def)                        :: lfric_fs = W0   ! W0
-    integer(kind=i_def)                        :: element_order = 1
+    integer(kind=i_def)                        :: element_order_h = 1
+    integer(kind=i_def)                        :: element_order_v = 1
     integer(kind=i_def)                        :: ndata_sz
 
     ! Use the unit-testing constructor:
@@ -81,9 +83,10 @@ program main
     mesh = mesh_type(global_mesh_ptr, partition, extrusion_ptr)
     write (*,*) "Mesh has", mesh%get_nlayers(), "layers."
     ndata_sz = 1
-    vector_space = function_space_type( mesh,          &
-                                        element_order, &
-                                        lfric_fs,      &
+    vector_space = function_space_type( mesh,            &
+                                        element_order_h, &
+                                        element_order_v, &
+                                        lfric_fs,        &
                                         ndata_sz)
     vector_space_ptr => vector_space
     call field1%initialise( vector_space = vector_space_ptr, name="field1" )

diff --git a/examples/lfric/eg17/full_example_extract/main.X90 b/examples/lfric/eg17/full_example_extract/main.X90
@@ -30,6 +30,7 @@
 ! -----------------------------------------------------------------------------
 ! Author: J. Henrichs, Bureau of Meteorology
 ! Modifications: A. R. Porter, STFC Daresbury Laboratory
+!                J. Dendy, Met Office
 
 program main
     !> This program is a simple LFRic program that contains two invokes,
@@ -70,7 +71,8 @@ program main
     type(field_type)                           :: field1, field2
     type(field_type)                           :: chi(3)
     integer(kind=i_def)                        :: lfric_fs = W0   ! W0
-    integer(kind=i_def)                        :: element_order = 1
+    integer(kind=i_def)                        :: element_order_h = 1
+    integer(kind=i_def)                        :: element_order_v = 1
     integer(kind=i_def)                        :: ndata_sz
     real(kind=r_def)                           :: one
     logical(kind=l_def)                        :: some_logical
@@ -101,9 +103,10 @@ program main
     mesh = mesh_type(global_mesh_ptr, partition, extrusion_ptr)
     write (*,*) "Mesh has", mesh%get_nlayers(), "layers."
     ndata_sz = 1
-    vector_space = function_space_type( mesh,          &
-                                        element_order, &
-                                        lfric_fs,      &
+    vector_space = function_space_type( mesh,            &
+                                        element_order_h, &
+                                        element_order_v, &
+                                        lfric_fs,        &
                                         ndata_sz)
     vector_space_ptr => vector_space
     do i=1, 3

diff --git a/examples/lfric/eg17/full_example_netcdf/README.md b/examples/lfric/eg17/full_example_netcdf/README.md
@@ -9,7 +9,7 @@ steps are required for this (using simplified code examples):
     ```fortran
     global_mesh = global_mesh_type("mesh_BiP128x16-400x100.nc", "dynamics")
     ```
-   
+
 2) A 1x1 planar partition for one process is created:
     ```fortran
     partitioner_ptr => partitioner_planar
@@ -34,9 +34,10 @@ steps are required for this (using simplified code examples):
 
 5) Create a function/vector space:
     ```fortran
-    vector_space = function_space_type( mesh,          &
-                                        element_order, &
-                                        lfric_fs,      &
+    vector_space = function_space_type( mesh,            &
+                                        element_order_h, &
+                                        element_order_v, &
+                                        lfric_fs,        &
                                         ndata_sz)
     ```
 
@@ -61,7 +62,7 @@ steps are required for this (using simplified code examples):
 
 ## Compilation
 
-A simple makefile is provided to compile the example. It needs 
+A simple makefile is provided to compile the example. It needs
 a full installation of NetCDF, since it is using ``nf-config`` to
 query the required compiler and linker flags, and the
 infrastructure library ``liblfric_netcdf.a`` provided in

diff --git a/examples/lfric/eg17/full_example_netcdf/main.x90 b/examples/lfric/eg17/full_example_netcdf/main.x90
@@ -30,6 +30,7 @@
 ! -----------------------------------------------------------------------------
 ! Author: J. Henrichs, Bureau of Meteorology
 ! Modifications: A. R. Porter, STFC Daresbury Laboratory
+!                J. Dendy, Met Office
 
 program main
 
@@ -58,7 +59,8 @@ program main
     procedure (partitioner_interface), pointer :: partitioner_ptr
     type(field_type)                           :: field1, field2
     integer(kind=i_def)                        :: lfric_fs = W0   ! W0
-    integer(kind=i_def)                        :: element_order = 1
+    integer(kind=i_def)                        :: element_order_h = 1
+    integer(kind=i_def)                        :: element_order_v = 1
     integer(kind=i_def)                        :: ndata_sz
 
     ! Use the unit-testing constructor:
@@ -82,9 +84,10 @@ program main
     mesh = mesh_type(global_mesh_ptr, partition, extrusion_ptr)
     write (*,*) "Mesh has", mesh%get_nlayers(), "layers."
     ndata_sz = 1
-    vector_space = function_space_type( mesh,          &
-                                        element_order, &
-                                        lfric_fs,      &
+    vector_space = function_space_type( mesh,            &
+                                        element_order_h, &
+                                        element_order_v, &
+                                        lfric_fs,        &
                                         ndata_sz)
     vector_space_ptr => vector_space
     call field1%initialise( vector_space = vector_space_ptr, name="field1" )

diff --git a/examples/lfric/eg4/solver_mod.x90 b/examples/lfric/eg4/solver_mod.x90
@@ -38,6 +38,7 @@
 ! -----------------------------------------------------------------------------
 ! Modified by A. Porter, STFC Daresbury Lab
 ! Modified by I. Kavcic, Met Office
+!             J. Dendy, Met Office
 !
 !> @brief Contains methods and algorithms for solving a system A.x = b for known
 !! input field b and matrix A and returns field x
@@ -51,7 +52,7 @@ module solver_mod
   use constants_mod,                 only : r_def, str_def, i_def
   use function_space_collection_mod, only : function_space_collection
   use field_mod,                     only : field_type
-  use finite_element_config_mod,     only : element_order
+  use finite_element_config_mod,     only : element_order_h, element_order_v
   use fs_continuity_mod,             only : W3
   use function_space_mod,            only : function_space_type
   use log_mod,                   only : log_event,         &
@@ -121,8 +122,9 @@ subroutine jacobi_solver_algorithm(lhs, rhs, mm, mesh, n_iter)
 
   flag_obj%rflag = 0.0_r_def
 
-  rhs_fs => function_space_collection%get_fs( mesh, &
-                                              element_order, &
+  rhs_fs => function_space_collection%get_fs( mesh,            &
+                                              element_order_h, &
+                                              element_order_v, &
                                               rhs%which_function_space() )
 
   diagonal = field_type( vector_space = rhs_fs )

diff --git a/examples/lfric/eg5/alg.f90 b/examples/lfric/eg5/alg.f90
@@ -39,6 +39,7 @@
 ! POSSIBILITY OF SUCH DAMAGE.
 !------------------------------------------------------------------------------
 ! Modified by I. Kavcic, Met Office
+!             J. Dendy, Met Office
 
 module oned_conservative_flux_alg_mod
 
@@ -48,7 +49,8 @@ module oned_conservative_flux_alg_mod
   use function_space_collection_mod,     only: function_space_collection
   use quadrature_mod,                    only: quadrature_type, GAUSSIAN
   use fs_continuity_mod,                 only: W0, W3
-  use finite_element_config_mod,         only: element_order
+  use finite_element_config_mod,         only: element_order_h,                &
+                                               element_order_v
   use subgrid_config_mod,                only: transport_stencil_length,       &
                                                rho_stencil_length
 
@@ -81,8 +83,9 @@ subroutine oned_conservative_flux_alg( direction,    &
 
     type(function_space_type), pointer :: rho_fs   => null()
 
-    rho_fs => function_space_collection%get_fs( mesh_id, element_order,      &
-                                              rho_in%which_function_space() )
+    rho_fs => function_space_collection%get_fs( mesh_id, element_order_h, &
+                                                element_order_v,          &
+                                                rho_in%which_function_space() )
 
     a0 = field_type( vector_space = rho_fs )
     a1 = field_type( vector_space = rho_fs )