1 implement non uniform compact finite difference schemes

CHANGELOG.md

@@ -6,12 +6,13 @@ This is the change log for `NuCFD`, it is based on the format suggested by [Keep
 
 ### Added
 
-- Added Cray compiler support [42b7070].
-- Added Intel compiler support [2d3ce07].
 - Added `FORD` documentation system [c780fe1].
 - Added tridiagonal solver with support for periodic problems [f5d254f].
 - Added simple library to support writing tests [a5814c8].
 - Added build and test system based on `cmake` [8d56f30].
+  - Added Cray compiler support [42b7070].
+  - Added Intel compiler support [2d3ce07].
+  - Set minimum `cmake` version `v3.25.2` to support Fortran submodules [e8c4aba].
 
 ### Changed
 ### Deprecated

CMakeLists.txt

@@ -2,16 +2,23 @@
 #
 ## Description
 #
+# Main CMakeLists file for the NuCFD project.
+#
 ## LICENSE
 #
 # SPDX-License-Identifier: BSD-3-Clause
 #
 
 ## Basic configuration
-cmake_minimum_required(VERSION 3.16)
+cmake_minimum_required(VERSION 3.25.2)
 project(NuCFD)
 enable_language(Fortran)
 
+option(BUILD_SHARED_LIBS "Enable building as a shared library" OFF)
+if(BUILD_SHARED_LIBS)
+  set(CMAKE_POSITION_INDEPENDENT_CODE true)
+endif()
+
 include(GNUInstallDirs)
 set(LIBRARY_OUTPUT_PATH ${PROJECT_BINARY_DIR}/${CMAKE_INSTALL_LIBDIR})
 set(EXECUTABLE_OUTPUT_PATH ${PROJECT_BINARY_DIR}/${CMAKE_INSTALL_BINDIR})
@@ -23,9 +30,9 @@ set(icc_like_fc "$<COMPILE_LANG_AND_ID:Fortran,Intel>")
 set(cray_like_fc "$<COMPILE_LANG_AND_ID:Fortran,Cray>")
 add_library(nucfd_compiler_flags INTERFACE)
 target_compile_options(nucfd_compiler_flags INTERFACE
-  "$<${gcc_like_fc}:-std=f2008;-Wall;-Wextra;-Wpedantic;-Warray-bounds;-Wimplicit-interface;-Wimplicit-procedure;-fimplicit-none>"
-  "$<${icc_like_fc}:-std08;-warn;all>"
-  "$<${cray_like_fc}:-en;-emf;-eI>")
+  "$<${gcc_like_fc}:-cpp;-std=f2018;-Wall;-Wextra;-Wpedantic;-Warray-bounds;-Wimplicit-interface;-Wimplicit-procedure;-fimplicit-none>"
+  "$<${icc_like_fc}:-fpp;-std18;-warn;all>"
+  "$<${cray_like_fc}:-cpp;-hfp0;-en;-emf;-eI>") # Anything greater than fp0 fails tests!
 add_library(nucfd_debug_flags INTERFACE)
 target_compile_options(nucfd_debug_flags INTERFACE
   "$<${gcc_like_fc}:-g;-Og;-ffpe-trap=invalid,zero;-fcheck=bounds;-fbacktrace>"
@@ -35,48 +42,16 @@ add_library(nucfd_strict_flags INTERFACE)
 target_compile_options(nucfd_strict_flags INTERFACE
   "$<${gcc_like_fc}:-Werror;-Wno-error=integer-division>"
   "$<${cray_like_fc}:-eN>")
+if (CMAKE_BUILD_TYPE MATCHES "Debug")
+  target_link_libraries(nucfd_compiler_flags INTERFACE nucfd_debug_flags)
+endif()
 
 ## Build nucfd library
-add_library(nucfd
-  src/nucfd_trid_solver_mod.f90
-  src/nucfd_coeffs_mod.f90)
-target_link_libraries(nucfd nucfd_compiler_flags)
+add_subdirectory(src)
 
 ## Testing
 enable_testing()
-
-add_library(nucfd_test_flags INTERFACE)
-target_link_libraries(nucfd_test_flags INTERFACE
-  nucfd_compiler_flags
-  nucfd_debug_flags
-  nucfd_strict_flags)
-
-add_library(nucfd_test_framework tests/nucfd_tests_mod.f90)
-target_link_libraries(nucfd_test_framework nucfd_test_flags)
-
-function(define_test suite test_name)
-  add_executable(${test_name} tests/${suite}/${test_name}.f90)
-  target_link_libraries(${test_name} nucfd)
-  target_link_libraries(${test_name} nucfd_test_framework)
-  target_link_libraries(${test_name} nucfd_test_flags)
-  add_test(NAME ${suite}:${test_name} COMMAND ${test_name})
-endfunction()
-
-function(add_test_libs suite test_name lib)
-  target_link_libraries(${test_name} ${lib})
-endfunction()
-
-# Tidiagonal solver test suite
-add_library(trid_test_utils tests/tridsolver/tridsol_test_utils_mod.f90)
-target_link_libraries(trid_test_utils nucfd_test_flags)
-define_test(tridsolver system_11_symm nucfd_test_framework)
-define_test(tridsolver system_11_anti-symm)
-add_test_libs(tridsolver system_11_symm trid_test_utils)
-add_test_libs(tridsolver system_11_anti-symm trid_test_utils)
-define_test(tridsolver system_00_symm nucfd_test_framework)
-define_test(tridsolver system_00_anti-symm)
-add_test_libs(tridsolver system_00_symm trid_test_utils)
-add_test_libs(tridsolver system_00_anti-symm trid_test_utils)
+add_subdirectory(tests)
 
 ## Documentation
 add_custom_target(doc ford

README.md

@@ -4,36 +4,49 @@ difference schemes on non-uniform grids.
 ## Building
 
 The `NuCFD` build system is generated by `cmake`, a default configuration can be created by running
-``
+```
 cmake -B build
-``
+```
 from the root directory.
 Once the build system has been generated running
-``
+```
 make -C build
-``
+```
 will build the library.
 The build can be configured using the `ccmake` tool, however currently only `gfortran` is supported.
 
+Note that this project uses Fortran2008 submodules, support for these in CMake requires CMake v3.25.2
+or above.
+If you already have an older CMake it can be relatively easily updated for the local user by cloning
+the repository and checking out tag `v3.25.2` (or later), CMake iteslf can then be built using a
+fairly standard process:
+```
+cmake -B build -DCMAKE_INSTALL_PREFIX=/path/to/install
+# Any configuration required
+make -C build && make -C build install
+export PATH=/path/to/install/bin:${PATH}
+```
+you should then be able to use your new CMake to configure the `NuCFD` build.
+
 ## Testing
 
 After building the library it can be tested using `ctest`.
 From the root directory run
-``
+```
 make -C build test
-``
+```
 to launch the tests.
 If any test fail more detailed output can be shown by running `ctest` verbosely
-``
+```
 cd build
 ctest --verbose
-``
+```
 which will print any output from the tests, add `--rerun-failed` to only repeat failing tests.
 
 ## Documentation
 
 Documentation can be generated using the FORD tool by running
-``
+```
 make -C build/ doc
-``
+```
 the generated documentation can be viewed in a web browser by opening `./manual/index.html`.

src/CMakeLists.txt

@@ -0,0 +1,24 @@
+## src/CMakeLists.txt
+#
+## Description
+#
+# CMakeLists file for building the NuCFD library.
+#
+## LICENSE
+#
+# SPDX-License-Identifier: BSD-3-Clause
+#
+
+add_library(nucfd_types OBJECT
+  nucfd_types_mod.f90
+)
+target_link_libraries(nucfd_types PRIVATE nucfd_compiler_flags)
+
+add_subdirectory(coeffs)
+
+add_library(nucfd
+  nucfd_trid_solver_mod.f90
+  nucfd_deriv_mod.f90)
+target_link_libraries(nucfd PRIVATE $<TARGET_OBJECTS:nucfd_types>)
+target_link_libraries(nucfd PRIVATE $<TARGET_OBJECTS:nucfd_coeffs>)
+target_link_libraries(nucfd PUBLIC nucfd_compiler_flags)

src/coeffs/CMakeLists.txt

@@ -0,0 +1,21 @@
+## src/coeffs/CMakeLists.txt
+#
+## Description
+#
+# CMakeLists file for building the coefficients module of NuCFD.
+#
+## LICENSE
+#
+# SPDX-License-Identifier: BSD-3-Clause
+#
+
+add_library(nucfd_coeffs OBJECT
+  nucfd_coeffs_mod.f90
+  nucfd_coeffs_a.f90
+  nucfd_coeffs_b.f90
+  nucfd_coeffs_c.f90
+  nucfd_coeffs_d.f90
+  nucfd_coeffs_e.f90
+)
+target_link_libraries(nucfd_coeffs PRIVATE $<TARGET_OBJECTS:nucfd_types>)
+target_link_libraries(nucfd_coeffs PRIVATE nucfd_compiler_flags)

src/coeffs/nucfd_coeffs_a.f90

@@ -0,0 +1,173 @@
+! src/coeffs/nucfd_coeff_a.f90
+!
+!! Implements coefficient A of the nucfd_coeffs module.
+!
+! SPDX-License-Identifier: BSD-3-Clause
+
+submodule (nucfd_coeffs) nucfd_coeffs_a
+  !! Submodule defining the coefficient acting on f_{i+1} for compact finite difference schemes on
+  !! non-uniform grids.
+
+  implicit none
+
+contains
+
+  module real function coeff_a_points(x)
+    !! Compute the coefficient acting on f_{i+1} of the finite diference given a stencil of points.
+
+    type(nucfd_stencil_points), intent(in) :: x !! Stencil of points for the finite
+                                                !! difference.
+
+    type(nucfd_stencil_deltas) :: h  ! Stencil of grid spacings for the finite difference.
+
+#ifndef NDEBUG
+    if (size(x%stencil) /= 5) then
+       print *, "Error@coeff_a_points: expecting width 5 stencil, received width = ", &
+            size(x%stencil)
+       error stop
+    end if
+    if (1 - lbound(x%stencil, 1) /= 3) then
+       print *, "Error@coeff_a_points: expecting centre 3 stencil, received centre = ", &
+            1 - lbound(x%stencil, 1)
+       print *, size(x%stencil), lbound(x%stencil), ubound(x%stencil)
+       error stop
+    end if
+#endif
+
+    h = points_to_deltas(x)
+    coeff_a_points = coeff_a(h)
+
+  end function coeff_a_points
+
+  module real function coeff_a_deltas(h)
+    !! Compute the coefficient acting on f_{i+1} of the finite diference given a stencil of grid
+    !! spacings.
+
+    type(nucfd_stencil_deltas), intent(in) :: h !! Stencil of grid spacings for the finite difference.
+
+    real :: numerator, numerator_corr, denominator, divisor
+
+#ifndef NDEBUG
+    if (size(h%stencil) /= 4) then
+       print *, "Error@coeff_a_deltas: expecting width 4 stencil, received width = ", &
+            size(h%stencil)
+       error stop
+    end if
+    if (1 - lbound(h%stencil, 1) /= 2) then
+       print *, "Error@coeff_a_deltas: expecting centre 2 stencil, received centre = ", &
+            1 - lbound(h%stencil, 1)
+       print *, size(h%stencil), lbound(h%stencil), ubound(h%stencil)
+       error stop
+    end if
+#endif
+
+    call coeff_a_components(h, numerator, numerator_corr, denominator, divisor)
+    coeff_a_deltas = ((numerator + numerator_corr) / denominator) / divisor
+
+  end function coeff_a_deltas
+
+  module subroutine coeff_a_components(h, numerator, numerator_corr, denominator, divisor)
+
+    type(nucfd_stencil_deltas), intent(in) :: h !! Stencil of grid spacings for the finite difference.
+
+    real, intent(out) :: numerator
+    real, intent(out) :: numerator_corr
+    real, intent(out) :: denominator
+    real, intent(out) :: divisor
+
+    real :: hm1, h0, hp1, hp2 ! Grid deltas at i -2, -1, 0, +1, +2
+
+#ifndef NDEBUG
+    if (size(h%stencil) /= 4) then
+       print *, "Error@coeff_a_components: expecting width 4 stencil, received width = ", &
+            size(h%stencil)
+       error stop
+    end if
+    if (1 - lbound(h%stencil, 1) /= 2) then
+       print *, "Error@coeff_a_components: expecting centre 2 stencil, received centre = ", &
+            1 - lbound(h%stencil, 1)
+       print *, size(h%stencil), lbound(h%stencil), ubound(h%stencil)
+       error stop
+    end if
+#endif
+
+    select type(deltas => h%stencil)
+    type is(real)
+       hm1 = deltas(-1)
+       h0 = deltas(0)
+       hp1 = deltas(1)
+       hp2 = deltas(2)
+    class default
+       error stop
+    end select
+
+    associate(beta => alpha) ! To match Gamet et al. (1999)
+      numerator = coeff_numerator(hm1, h0, hp2, hp2, beta)
+      numerator_corr = coeff_numerator_corr(hm1, h0, hp1, hp2, alpha, beta)
+      denominator = coeff_denominator(h0, hp1, hp2)
+      divisor = coeff_divisor(hm1, h0, hp1)
+    end associate
+
+  end subroutine coeff_a_components
+
+  pure real function coeff_numerator(hm1, h0, hp1, hp2, beta)
+    !! Computes the numerator of the coefficient acting on f_{i+1}.
+    !!
+    !! Reduces to (14/3) h^3 when h=const and beta=1/3.
+
+    real, intent(in) :: hm1 
+    real, intent(in) :: h0 
+    real, intent(in) :: hp1 
+    real, intent(in) :: hp2
+    real, intent(in) :: beta
+
+    coeff_numerator = h0 * ((hm1 + h0) * (hp1 + hp2) + 2.0 * hp1 * hp2 * beta) ! = (14/3) h^3
+  end function coeff_numerator
+
+  pure real function coeff_numerator_corr(hm1, h0, hp1, hp2, alpha, beta)
+    !! Computes the non-uniform correction to the numerator acting on f_{i+1}.
+    !!
+    !! Reduces to zero when h=const and alpha=beta.
+
+    real, intent(in) :: hm1 
+    real, intent(in) :: h0 
+    real, intent(in) :: hp1 
+    real, intent(in) :: hp2
+    real, intent(in) :: alpha
+    real, intent(in) :: beta
+
+    coeff_numerator_corr =  hm1 * h0 * hp2 * (beta - alpha) & ! Should cancel for case alpha = beta
+         + (h0**2) * (hp2 - hp1) * beta                     & ! Should cancel for constant h
+         + hm1 * hp1 * (2.0 * hp2 - h0 - hp1) * beta        & ! Should cancel for constant h
+         + (hp1**2) * (3.0 * hp2 - 2.0 * h0 - hp1) * beta   & ! Should cancel for constant h
+         + h0 * (2.0 * hp1 * hp2 * beta - hm1 * (h0 + hp1) * alpha) ! Should cancel for constant h,
+                                                                    ! alpha=beta
+  end function coeff_numerator_corr
+
+  pure real function coeff_denominator(h0, hp1, hp2)
+    !! Computes the denominator of the coefficient acting on f_{i+1}.
+    !!
+    !! Reduces to 3 h^3 when h=const. Dividing the numerator by this term yields the coefficient
+    !! 14/9 when h=const, alpha=beta=1/3.
+
+    real, intent(in) :: h0 
+    real, intent(in) :: hp1 
+    real, intent(in) :: hp2 
+
+    coeff_denominator = (3.0 * hp1 * ((h0 + hp1) / 2.0) * hp2)
+  end function coeff_denominator
+
+  pure real function coeff_divisor(hm1, h0, hp1)
+    !! Computes the non-uniform equivalent to 2h divisor of the coefficient acting on f_{i+1}.
+    !!
+    !! Dividing the coefficient by this term should reduce to (14/9)/(2h) when h=const,
+    !! alpha=beta=1/3.
+
+    real, intent(in) :: hm1 
+    real, intent(in) :: h0 
+    real, intent(in) :: hp1 
+
+    coeff_divisor = (2.0 * ((hm1 + h0 + hp1) / 3.0))
+  end function coeff_divisor
+
+end submodule nucfd_coeffs_a