Add CI workflows and implement unit tests for MPAS dynamical core

.gitattributes

@@ -0,0 +1 @@
+*.pf linguist-language=Fortran-Free-Form

.github/workflows/mpas_dynamical_core_ci.yml

@@ -0,0 +1,95 @@
+name: MPAS Dynamical Core CI
+
+permissions:
+  contents: read
+
+on:
+  pull_request:
+    paths:
+      # Only run this workflow if something changes in MPAS dynamical core. Skip otherwise.
+      - 'src/dynamics/mpas/**'
+    types:
+      - opened
+      - reopened
+      - synchronize
+  push:
+    branches:
+      - 'staging/**'
+      - development
+      - main
+    paths:
+      # Only run this workflow if something changes in MPAS dynamical core. Skip otherwise.
+      - 'src/dynamics/mpas/**'
+  workflow_dispatch:
+
+concurrency:
+  cancel-in-progress: true
+  group: ${{ github.workflow }} - ${{ github.ref }}
+
+jobs:
+  unit-tests:
+    name: Build and run unit tests (GCC ${{ matrix.gcc-version }})
+    timeout-minutes: 10
+    strategy:
+      fail-fast: false
+      matrix:
+        gcc-version:
+          - '12'
+          - '13'
+          - '14'
+    runs-on: ubuntu-24.04
+    env:
+      CC: gcc-${{ matrix.gcc-version }}
+      CXX: g++-${{ matrix.gcc-version }}
+      FC: gfortran-${{ matrix.gcc-version }}
+      PFUNIT_BUILD_TYPE: Release
+      PFUNIT_PATH: ${{ github.workspace }}/pFUnit
+      PFUNIT_REFERENCE: 'v4.12.0'
+      UNIT_TESTS_BUILD_TYPE: Debug
+      UNIT_TESTS_PATH: ${{ github.workspace }}/src/dynamics/mpas
+    steps:
+      - name: Checkout CAM-SIMA
+        uses: actions/checkout@v5
+
+      # pFUnit takes a while to build. Cache it if possible to save time on consecutive workflow runs.
+      - name: Cache pFUnit
+        id: cache-pfunit
+        uses: actions/cache@v4
+        with:
+          key: cache-pfunit-${{ env.PFUNIT_REFERENCE }}-gcc-v${{ matrix.gcc-version }}
+          path: ${{ env.PFUNIT_PATH }}/install
+
+      # If there is a cache hit, just use the cached pFUnit and skip this step.
+      - name: Checkout pFUnit
+        if: ${{ steps.cache-pfunit.outputs.cache-hit != 'true' }}
+        uses: actions/checkout@v5
+        with:
+          repository: Goddard-Fortran-Ecosystem/pFUnit
+          ref: ${{ env.PFUNIT_REFERENCE }}
+          path: pFUnit
+
+      # If there is a cache hit, just use the cached pFUnit and skip this step.
+      - name: Build pFUnit
+        if: ${{ steps.cache-pfunit.outputs.cache-hit != 'true' }}
+        run: |
+          cmake -D CMAKE_BUILD_TYPE="$PFUNIT_BUILD_TYPE" -D CMAKE_INSTALL_PREFIX="$PFUNIT_PATH/install" -B "$PFUNIT_PATH/build" -S "$PFUNIT_PATH"
+          cmake --build "$PFUNIT_PATH/build" --config "$PFUNIT_BUILD_TYPE"
+          cmake --install "$PFUNIT_PATH/build" --config "$PFUNIT_BUILD_TYPE" --prefix "$PFUNIT_PATH/install"
+
+      - name: Build unit tests
+        run: |
+          cmake -D CMAKE_BUILD_TYPE="$UNIT_TESTS_BUILD_TYPE" -D CMAKE_PREFIX_PATH="$PFUNIT_PATH/install" -B "$UNIT_TESTS_PATH/build" -S "$UNIT_TESTS_PATH"
+          cmake --build "$UNIT_TESTS_PATH/build" --config "$UNIT_TESTS_BUILD_TYPE"
+
+      - name: Run unit tests
+        run: |
+          ctest --build-config "$UNIT_TESTS_BUILD_TYPE" --output-log "$UNIT_TESTS_PATH/build/unit-tests.log" --output-on-failure --test-dir "$UNIT_TESTS_PATH/build" --verbose
+
+      - name: Upload unit tests log
+        if: ${{ always() }}
+        uses: actions/upload-artifact@v4
+        with:
+          if-no-files-found: ignore
+          name: unit-tests-log-gcc-v${{ matrix.gcc-version }}
+          path: ${{ env.UNIT_TESTS_PATH }}/build/unit-tests.log
+          retention-days: 7

src/dynamics/mpas/CMakeLists.txt

@@ -0,0 +1,41 @@
+cmake_minimum_required(VERSION 3.21)
+
+# `mpas_dynamical_core_procedures` has not been integrated into the CMake build of any top level projects yet,
+# and this CMakeLists.txt file is currently for unit testing purposes only.
+# Therefore, making a change to this CMakeLists.txt file will not impact the build of a parent project at this time.
+project(mpas_dynamical_core_procedures
+    VERSION
+        0.1.0
+    DESCRIPTION
+        "Standardized computational and utility procedures in MPAS dynamical core"
+    LANGUAGES
+        Fortran
+)
+
+add_library(mpas_dynamical_core_procedures)
+target_sources(mpas_dynamical_core_procedures
+    PRIVATE
+        driver/dyn_mpas_procedures.F90
+        dyn_procedures.F90
+)
+target_include_directories(mpas_dynamical_core_procedures
+    INTERFACE
+        ${CMAKE_CURRENT_BINARY_DIR}
+)
+target_compile_options(mpas_dynamical_core_procedures
+    PRIVATE
+        $<$<AND:$<CONFIG:Debug>,$<Fortran_COMPILER_ID:GNU>>:-fbacktrace -fcheck=all -ffpe-trap=invalid,overflow,zero -fno-unsafe-math-optimizations -frounding-math -fsignaling-nans -std=f2018 -Wall -Wextra -Wpedantic>
+        $<$<AND:$<CONFIG:Debug>,$<Fortran_COMPILER_ID:Intel>>:-check all -fp-model=precise -stand f18 -traceback -warn all>
+        $<$<AND:$<CONFIG:Debug>,$<Fortran_COMPILER_ID:IntelLLVM>>:-check all -fp-model=precise -stand f18 -traceback -warn all>
+)
+
+# Enable testing only when this project is at the top level.
+# Otherwise, if this project is under a parent project, respect its choice of whether to enable testing.
+if(PROJECT_IS_TOP_LEVEL)
+    include(CTest)
+    enable_testing()
+endif(PROJECT_IS_TOP_LEVEL)
+
+if(BUILD_TESTING)
+    add_subdirectory(tests/unit)
+endif(BUILD_TESTING)

src/dynamics/mpas/driver/dyn_mpas_procedures.F90

@@ -95,6 +95,22 @@ pure elemental function almost_equal_real32(a, b, absolute_tolerance, relative_t
 
         real(real32) :: error_tolerance
 
+        if (a /= a .or. b /= b) then
+            ! NaN does not equal to anything, including itself.
+            almost_equal = .false.
+
+            return
+        end if
+
+        if (abs(a) > huge(a) .or. abs(b) > huge(b)) then
+            ! Infinity of the same sign equals to each other.
+            ! Infinity of different signs does not equal to each other.
+            ! Infinity does not equal to anything that is finite.
+            almost_equal = (a == b)
+
+            return
+        end if
+
         if (present(relative_tolerance)) then
             error_tolerance = relative_tolerance * max(abs(a), abs(b))
         else
@@ -103,15 +119,11 @@ pure elemental function almost_equal_real32(a, b, absolute_tolerance, relative_t
 
         if (present(absolute_tolerance)) then
             error_tolerance = max(absolute_tolerance, error_tolerance)
+        else
+            error_tolerance = max(epsilon(1.0_real32), error_tolerance)
         end if
 
-        if (abs(a - b) <= error_tolerance) then
-            almost_equal = .true.
-
-            return
-        end if
-
-        almost_equal = .false.
+        almost_equal = (abs(a - b) <= error_tolerance)
     end function almost_equal_real32
 
     !> Test `a` and `b` for approximate equality, where `a` and `b` are both reals.
@@ -125,6 +137,22 @@ pure elemental function almost_equal_real64(a, b, absolute_tolerance, relative_t
 
         real(real64) :: error_tolerance
 
+        if (a /= a .or. b /= b) then
+            ! NaN does not equal to anything, including itself.
+            almost_equal = .false.
+
+            return
+        end if
+
+        if (abs(a) > huge(a) .or. abs(b) > huge(b)) then
+            ! Infinity of the same sign equals to each other.
+            ! Infinity of different signs does not equal to each other.
+            ! Infinity does not equal to anything that is finite.
+            almost_equal = (a == b)
+
+            return
+        end if
+
         if (present(relative_tolerance)) then
             error_tolerance = relative_tolerance * max(abs(a), abs(b))
         else
@@ -133,15 +161,11 @@ pure elemental function almost_equal_real64(a, b, absolute_tolerance, relative_t
 
         if (present(absolute_tolerance)) then
             error_tolerance = max(absolute_tolerance, error_tolerance)
+        else
+            error_tolerance = max(epsilon(1.0_real64), error_tolerance)
         end if
 
-        if (abs(a - b) <= error_tolerance) then
-            almost_equal = .true.
-
-            return
-        end if
-
-        almost_equal = .false.
+        almost_equal = (abs(a - b) <= error_tolerance)
     end function almost_equal_real64
 
     !> Clamp/Limit the value of `x` to the range of [`xmin`, `xmax`], where `x`, `xmin`, and `xmax` are all integers.

src/dynamics/mpas/dyn_comp_impl.F90

@@ -330,6 +330,8 @@ subroutine check_topography_data(pio_file)
         use shr_kind_mod, only: kind_r8 => shr_kind_r8
         ! Module(s) from external libraries.
         use pio, only: file_desc_t, pio_file_is_open
+        ! Module(s) from MPAS.
+        use dyn_mpas_procedures, only: almost_equal
 
         type(file_desc_t), pointer, intent(in) :: pio_file
 
@@ -373,7 +375,9 @@ subroutine check_topography_data(pio_file)
             surface_geometric_height(:) = surface_geopotential(:) / constant_g
 
             ! Surface geometric height in MPAS should match the values in topography file.
-            if (any(abs(real(zgrid(1, 1:ncells_solve), kind_r8) - surface_geometric_height(:)) > error_tolerance)) then
+            if (.not. all(almost_equal( &
+                real(zgrid(1, 1:ncells_solve), kind_r8), surface_geometric_height, &
+                relative_tolerance=error_tolerance))) then
                 call endrun('Surface geometric height in MPAS is not consistent with topography data', subname, __LINE__)
             end if
 
@@ -383,7 +387,9 @@ subroutine check_topography_data(pio_file)
             call dyn_debug_print(debugout_verbose, 'Topography file is not used for consistency check')
 
             ! Surface geometric height in MPAS should be zero.
-            if (any(abs(real(zgrid(1, 1:ncells_solve), kind_r8)) > error_tolerance)) then
+            if (.not. all(almost_equal( &
+                real(zgrid(1, 1:ncells_solve), kind_r8), 0.0_kind_r8, &
+                relative_tolerance=error_tolerance))) then
                 call endrun('Surface geometric height in MPAS is not zero', subname, __LINE__)
             end if
         end if
@@ -587,7 +593,7 @@ subroutine set_mpas_state_scalars()
             use cam_constituents, only: num_advected
             use cam_logfile, only: debugout_verbose
             use dyn_grid, only: ncells_solve
-            use dyn_procedures, only: reverse
+            use dyn_procedures, only: qv_of_sh, reverse
             use dyn_tests_utils, only: vc_height
             use inic_analytic, only: dyn_set_inic_col
             use vert_coord, only: pver
@@ -641,8 +647,8 @@ subroutine set_mpas_state_scalars()
                 call dyn_debug_print(debugout_verbose, 'Conversion is needed and applied for water vapor mixing ratio')
 
                 ! Convert specific humidity to water vapor mixing ratio.
-                scalars(index_qv, :, 1:ncells_solve) = &
-                    scalars(index_qv, :, 1:ncells_solve) / (1.0_kind_dyn_mpas - scalars(index_qv, :, 1:ncells_solve))
+                scalars(index_qv, :, 1:ncells_solve) = real(qv_of_sh( &
+                    real(scalars(index_qv, :, 1:ncells_solve), kind_r8)), kind_dyn_mpas)
             end if
 
             deallocate(buffer_3d_real)

src/dynamics/mpas/dyn_coupling_impl.F90

@@ -361,10 +361,11 @@ end subroutine final_shared_variables
         !> (KCW, 2024-07-30)
         subroutine update_shared_variables(i)
             ! Module(s) from CAM-SIMA.
-            use dyn_mpas_procedures, only: clamp
             use dyn_procedures, only: dp_by_hydrostatic_equation, omega_of_w_rho, p_by_equation_of_state, t_of_tm_qv
             use dynconst, only: constant_g => gravit, constant_rd => rair, constant_rv => rh2o
             use vert_coord, only: pver, pverp
+            ! Module(s) from MPAS.
+            use dyn_mpas_procedures, only: clamp
 
             integer, intent(in) :: i
 
@@ -489,6 +490,7 @@ subroutine set_physics_state_external()
             use cam_thermo, only: cam_thermo_dry_air_update, cam_thermo_water_update
             use cam_thermo_formula, only: energy_formula_dycore_mpas
             use dyn_comp, only: mpas_dynamical_core
+            use dyn_procedures, only: exner_function
             use dynconst, only: constant_g => gravit
             use physics_types, only: cappav, cp_or_cv_dycore, cpairv, lagrangian_vertical, phys_state, rairv, zvirv
             use runtime_obj, only: cam_runtime_opts
@@ -553,7 +555,7 @@ subroutine set_physics_state_external()
             ! the paragraph below equation 1.5.1c in doi:10.1007/978-94-009-3027-8.
             ! Also note that `cappav` is updated externally by `cam_thermo_dry_air_update`.
             do i = 1, ncells_solve
-                phys_state % exner(i, :) = (phys_state % ps(i) / phys_state % pmid(i, :)) ** cappav(i, :)
+                phys_state % exner(i, :) = 1.0_kind_r8 / exner_function(cappav(i, :), phys_state % ps(i), phys_state % pmid(i, :))
             end do
 
             ! Note that constituents become moist after this.

src/dynamics/mpas/dyn_procedures.F90

@@ -28,6 +28,11 @@ module dyn_procedures
     ! Utility procedures.
     public :: reverse
     public :: sec_to_hour_min_sec
+
+    interface exner_function
+        module procedure exner_function_of_cpd_p0_rd_p
+        module procedure exner_function_of_kappa_p0_p
+    end interface exner_function
 contains
     !> Compute the pressure `p` from the density `rho` and the temperature `t` by equation of state. Essentially,
     !> \( P = \rho R T \). Equation of state may take other forms, such as \( P_d = \rho_d R_d T \), \( P = \rho R_d T_v \),
@@ -69,14 +74,24 @@ pure elemental function t_by_equation_of_state(constant_r, p, rho) result(t)
     end function t_by_equation_of_state
 
     !> Compute the Exner function `pi` from the pressure `p`. Essentially, \( \Pi = (\frac{P}{P_0})^{\frac{R_d}{C_{pd}}} \).
-    pure elemental function exner_function(constant_cpd, constant_p0, constant_rd, p) result(pi)
+    pure elemental function exner_function_of_cpd_p0_rd_p(constant_cpd, constant_p0, constant_rd, p) result(pi)
         use, intrinsic :: iso_fortran_env, only: real64
 
         real(real64), intent(in) :: constant_cpd, constant_p0, constant_rd, p
         real(real64) :: pi
 
         pi = (p / constant_p0) ** (constant_rd / constant_cpd)
-    end function exner_function
+    end function exner_function_of_cpd_p0_rd_p
+
+    !> Compute the Exner function `pi` from the pressure `p`. Essentially, \( \Pi = (\frac{P}{P_0})^{\kappa} \).
+    pure elemental function exner_function_of_kappa_p0_p(constant_kappa, constant_p0, p) result(pi)
+        use, intrinsic :: iso_fortran_env, only: real64
+
+        real(real64), intent(in) :: constant_kappa, constant_p0, p
+        real(real64) :: pi
+
+        pi = (p / constant_p0) ** constant_kappa
+    end function exner_function_of_kappa_p0_p
 
     !> Compute the pressure difference `dp` from the density `rho` and the height difference `dz` by hydrostatic equation.
     !> Essentially, \( \mathrm{d} P = -\rho g \mathrm{d} z \).
@@ -348,12 +363,14 @@ end function reverse
     !> Convert second(s) to hour(s), minute(s), and second(s).
     !> (KCW, 2024-02-07)
     pure function sec_to_hour_min_sec(sec) result(hour_min_sec)
-        integer, intent(in) :: sec
-        integer :: hour_min_sec(3)
+        use, intrinsic :: iso_fortran_env, only: int32
+
+        integer(int32), intent(in) :: sec
+        integer(int32) :: hour_min_sec(3)
 
-        ! These are all intended to be integer arithmetics.
-        hour_min_sec(1) = sec / 3600
-        hour_min_sec(2) = sec / 60 - hour_min_sec(1) * 60
-        hour_min_sec(3) = sec - hour_min_sec(1) * 3600 - hour_min_sec(2) * 60
+        ! These are all intended to be integer arithmetic.
+        hour_min_sec(1) = sec / 3600_int32
+        hour_min_sec(2) = sec / 60_int32 - hour_min_sec(1) * 60_int32
+        hour_min_sec(3) = sec - hour_min_sec(1) * 3600_int32 - hour_min_sec(2) * 60_int32
     end function sec_to_hour_min_sec
 end module dyn_procedures

src/dynamics/mpas/tests/unit/CMakeLists.txt

@@ -0,0 +1,27 @@
+find_package(PFUNIT REQUIRED)
+
+add_pfunit_ctest(test_dyn_mpas_procedures
+    TEST_SOURCES
+        test_dyn_mpas_procedures.pf
+    LINK_LIBRARIES
+        mpas_dynamical_core_procedures
+)
+target_compile_options(test_dyn_mpas_procedures
+    PRIVATE
+        $<$<AND:$<CONFIG:Debug>,$<Fortran_COMPILER_ID:GNU>>:-fbacktrace -fcheck=all -ffpe-trap=invalid,overflow,zero -fno-unsafe-math-optimizations -frounding-math -fsignaling-nans -std=f2018 -Wall -Wextra -Wpedantic>
+        $<$<AND:$<CONFIG:Debug>,$<Fortran_COMPILER_ID:Intel>>:-check all -fp-model=precise -stand f18 -traceback -warn all>
+        $<$<AND:$<CONFIG:Debug>,$<Fortran_COMPILER_ID:IntelLLVM>>:-check all -fp-model=precise -stand f18 -traceback -warn all>
+)
+
+add_pfunit_ctest(test_dyn_procedures
+    TEST_SOURCES
+        test_dyn_procedures.pf
+    LINK_LIBRARIES
+        mpas_dynamical_core_procedures
+)
+target_compile_options(test_dyn_procedures
+    PRIVATE
+        $<$<AND:$<CONFIG:Debug>,$<Fortran_COMPILER_ID:GNU>>:-fbacktrace -fcheck=all -ffpe-trap=invalid,overflow,zero -fno-unsafe-math-optimizations -frounding-math -fsignaling-nans -std=f2018 -Wall -Wextra -Wpedantic>
+        $<$<AND:$<CONFIG:Debug>,$<Fortran_COMPILER_ID:Intel>>:-check all -fp-model=precise -stand f18 -traceback -warn all>
+        $<$<AND:$<CONFIG:Debug>,$<Fortran_COMPILER_ID:IntelLLVM>>:-check all -fp-model=precise -stand f18 -traceback -warn all>
+)