From 0a3c601f9ff2c8e96ec9dcae3127ddddafd95495 Mon Sep 17 00:00:00 2001 From: Paul Date: Sat, 25 Feb 2023 23:26:56 +0000 Subject: [PATCH 01/58] Start implementing first coefficient of non-uniform FD scheme This coefficient acts on f_{i+1} - coefficient A following Gamet et al. (1999). --- CMakeLists.txt | 3 +++ src/nucfd_coeffs_mod.f90 | 11 ++++++++- tests/fd-schemes/verify_coeffs.f90 | 39 ++++++++++++++++++++++++++++++ tests/nucfd_tests_mod.f90 | 23 ++++++++++++++++++ 4 files changed, 75 insertions(+), 1 deletion(-) create mode 100644 tests/fd-schemes/verify_coeffs.f90 diff --git a/CMakeLists.txt b/CMakeLists.txt index ab89ffe..374ba62 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -65,3 +65,6 @@ 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) + +# Finite difference test suite +define_test(fd-schemes verify_coeffs) diff --git a/src/nucfd_coeffs_mod.f90 b/src/nucfd_coeffs_mod.f90 index 626eb9b..54485aa 100644 --- a/src/nucfd_coeffs_mod.f90 +++ b/src/nucfd_coeffs_mod.f90 @@ -16,8 +16,17 @@ module nucfd_coeffs implicit none private - + public :: coeff_a + real, parameter, public :: alpha = 1.0 / 3.0 ! Off-diagonal coefficient for first-derivative ! system. + +contains + + pure real function coeff_a() + + coeff_a = 0.0 + + end function coeff_a end module nucfd_coeffs diff --git a/tests/fd-schemes/verify_coeffs.f90 b/tests/fd-schemes/verify_coeffs.f90 new file mode 100644 index 0000000..0a95bc9 --- /dev/null +++ b/tests/fd-schemes/verify_coeffs.f90 @@ -0,0 +1,39 @@ +!!!! tests/fd-schemes/verify_coeffs.f90 +!!! +!!!! Description +!!! +!!! Part of the fd-schemes test suite. +!!! Tests the computation of finite difference coefficients for non-uniform grids. +!!! +!!!! LICENSE +!!! +!!! SPDX-License-Identifier: BSD-3-Clause +!!! + +program verify_coeffs + + use nucfd_coeffs + + use nucfd_tests + + implicit none + + real :: n + real :: L + real :: h + + real :: a + real, parameter :: aref = 14.0 / 9.0 + + call initialise_suite("Verify coefficients") + + n = 128 + L = 1.0 + h = L / real(n - 1) + + a = coeff_a() + call test_report("Coefficient A", check_scalar(a, aref / (2.0 * h))) + + call finalise_suite() + +end program verify_coeffs diff --git a/tests/nucfd_tests_mod.f90 b/tests/nucfd_tests_mod.f90 index dee8edf..7c38db6 100644 --- a/tests/nucfd_tests_mod.f90 +++ b/tests/nucfd_tests_mod.f90 @@ -19,6 +19,7 @@ module nucfd_tests public :: initialise_suite, finalise_suite public :: test_report public :: check_rms + public :: check_scalar character(len=:), allocatable :: suite_name logical, save :: passing @@ -73,6 +74,28 @@ subroutine test_report(test_name, test_status) end subroutine test_report + logical function check_scalar(test, ref) + ! Compute error and report errors. + + real, intent(in) :: test ! The test data + real, intent(in) :: ref ! The reference data + + real :: err + logical :: test_passing + + err = abs(test - ref) + if (err > (2 * epsilon(ref))) then + test_passing = .false. + + print *, "Delta = ", err, " exceeds tolerance: ", 2 * epsilon(ref) + else + test_passing = .true. + end if + + check_scalar = test_passing + + end function check_scalar + logical function check_rms(test, ref) ! Compute RMS of error and report errors. From 29aaf14c01fd59036e8a61c5a52a3336b3b5fb22 Mon Sep 17 00:00:00 2001 From: Paul Date: Sat, 25 Feb 2023 23:41:28 +0000 Subject: [PATCH 02/58] Describe paramters of verify coefficients test. --- tests/fd-schemes/verify_coeffs.f90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/tests/fd-schemes/verify_coeffs.f90 b/tests/fd-schemes/verify_coeffs.f90 index 0a95bc9..7ccdff5 100644 --- a/tests/fd-schemes/verify_coeffs.f90 +++ b/tests/fd-schemes/verify_coeffs.f90 @@ -18,9 +18,9 @@ program verify_coeffs implicit none - real :: n - real :: L - real :: h + real :: n ! Mesh size + real :: L ! Domain size + real :: h ! Average grid spacing. real :: a real, parameter :: aref = 14.0 / 9.0 From 243183d55bcbe3f8172c3426bbb482298875add8 Mon Sep 17 00:00:00 2001 From: Paul Date: Sun, 26 Feb 2023 01:45:45 +0000 Subject: [PATCH 03/58] Fleshing out the first coefficient implementation --- CMakeLists.txt | 1 + src/nucfd_coeffs_mod.f90 | 31 +++++++++++++++++--- src/nucfd_types_mod.f90 | 46 ++++++++++++++++++++++++++++++ tests/fd-schemes/verify_coeffs.f90 | 5 +++- 4 files changed, 78 insertions(+), 5 deletions(-) create mode 100644 src/nucfd_types_mod.f90 diff --git a/CMakeLists.txt b/CMakeLists.txt index f3df556..5ee7d6b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -27,6 +27,7 @@ target_compile_options(nucfd_strict_flags INTERFACE ## Build nucfd library add_library(nucfd src/nucfd_trid_solver_mod.f90 + src/nucfd_types_mod.f90 src/nucfd_coeffs_mod.f90) target_link_libraries(nucfd nucfd_compiler_flags) diff --git a/src/nucfd_coeffs_mod.f90 b/src/nucfd_coeffs_mod.f90 index 47728a2..05c374c 100644 --- a/src/nucfd_coeffs_mod.f90 +++ b/src/nucfd_coeffs_mod.f90 @@ -13,22 +13,45 @@ module nucfd_coeffs !! Module defining the coefficients for compact finite differenceschemes on non-uniform grids. + + use nucfd_types implicit none private public :: coeff_a + interface coeff_a + module procedure coeff_a_points + module procedure coeff_a_deltas + end interface coeff_a + real, parameter, public :: alpha = 1.0 / 3.0 !! Off-diagonal coefficient for first derivative !! system. contains - pure real function coeff_a() - !! Compute the coefficient acting on f_{i+1} of the finite diference. + pure 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(6, 4)), intent(in) :: x !! Stencil of points for the finite + !! difference. + + type(nucfd_stencil_deltas(5, 3)) :: h !! Stencil of grid spacings for the finite difference. + + h = points_to_deltas(x) + coeff_a_points = coeff_a(h) - coeff_a = 0.0 + end function coeff_a_points + + pure 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(5, 3)), intent(in) :: h !! Stencil of grid spacings for the finite difference. + + coeff_a_deltas = 0.0 - end function coeff_a + end function coeff_a_deltas end module nucfd_coeffs diff --git a/src/nucfd_types_mod.f90 b/src/nucfd_types_mod.f90 new file mode 100644 index 0000000..7a23460 --- /dev/null +++ b/src/nucfd_types_mod.f90 @@ -0,0 +1,46 @@ +!!!! src/nucfd_types_mod.f90 +!!! +!!!! Description +!!! +!!! Defines the custom types used by NuCFD. +!!! +!!! Provides the nucfd_types module. +!!! +!!!! LICENSE +!!! +!!! SPDX-License-Identifier: BSD-3-Clause +!!! + +module nucfd_types + + implicit none + + type nucfd_stencil(width, centre) + !! Type representing a stencil. + integer, kind :: width = 1 !! Stencil width. + integer, kind :: centre = 1 !! Stencil centre. + real, dimension(1-centre:(1-centre)+(width-1)) :: stencil !! Stencil data. The lbound, ubound + !! are set to allow stencil access by + !! offsets. + end type nucfd_stencil + + type, extends(nucfd_stencil) :: nucfd_stencil_points + !! Type representing the point coordinates of a stencil. + end type nucfd_stencil_points + + type, extends(nucfd_stencil) :: nucfd_stencil_deltas + !! Type representing the grid spacings of a stencil. Following Gamet et al. (1999) these are + !! defined as h_i = x_i - x_{i-1}. + end type nucfd_stencil_deltas + +contains + + pure function points_to_deltas(x) result(h) + !! Helper function converting a stencil's point representation to a grid spacing representation. + + type(nucfd_stencil_points(6, 4)), intent(in) :: x !! The stencil's point representation. + type(nucfd_stencil_deltas(5, 3)) :: h !! The stencil's grid spacing representation. + + end function points_to_deltas + +end module nucfd_types diff --git a/tests/fd-schemes/verify_coeffs.f90 b/tests/fd-schemes/verify_coeffs.f90 index 7ccdff5..80d05bf 100644 --- a/tests/fd-schemes/verify_coeffs.f90 +++ b/tests/fd-schemes/verify_coeffs.f90 @@ -12,12 +12,15 @@ program verify_coeffs + use nucfd_types use nucfd_coeffs use nucfd_tests implicit none + type(nucfd_stencil_points(6, 4)) :: stencil ! Stencil of grid points. + real :: n ! Mesh size real :: L ! Domain size real :: h ! Average grid spacing. @@ -31,7 +34,7 @@ program verify_coeffs L = 1.0 h = L / real(n - 1) - a = coeff_a() + a = coeff_a(stencil) call test_report("Coefficient A", check_scalar(a, aref / (2.0 * h))) call finalise_suite() From 9d4e9c1dd17db727ef62a2322f2da61c2dbf03c7 Mon Sep 17 00:00:00 2001 From: Paul Date: Sun, 26 Feb 2023 10:04:53 +0000 Subject: [PATCH 04/58] Add debugging flags for Debug builds --- CMakeLists.txt | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index aa8a505..a11e24e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -24,7 +24,7 @@ target_compile_options(nucfd_compiler_flags INTERFACE "$<${gcc_like_fc}:-std=f2008;-Wall;-Wextra;-Warray-bounds;-Wimplicit-interface;-Wimplicit-procedure;-fimplicit-none>") add_library(nucfd_debug_flags INTERFACE) target_compile_options(nucfd_debug_flags INTERFACE - "$<${gcc_like_fc}:-ffpe-trap=invalid,zero;-fcheck=bounds>") + "$<${gcc_like_fc}:-ffpe-trap=invalid,zero;-fcheck=bounds;-fbacktrace>") add_library(nucfd_strict_flags INTERFACE) target_compile_options(nucfd_strict_flags INTERFACE "$<${gcc_like_fc}:-Werror>") @@ -35,6 +35,9 @@ add_library(nucfd src/nucfd_types_mod.f90 src/nucfd_coeffs_mod.f90) target_link_libraries(nucfd nucfd_compiler_flags) +if (CMAKE_BUILD_TYPE MATCHES "Debug") + target_link_libraries(nucfd nucfd_debug_flags) +endif() ## Testing enable_testing() From edc4a6efd74626f6e2aab8c0099083e27e6a6a80 Mon Sep 17 00:00:00 2001 From: Paul Date: Sun, 26 Feb 2023 10:24:58 +0000 Subject: [PATCH 05/58] gfortran doesn't support parameterised types Atleast not as used here. --- src/nucfd_coeffs_mod.f90 | 30 ++++++++++---- src/nucfd_types_mod.f90 | 64 +++++++++++++++++++++++++----- tests/fd-schemes/verify_coeffs.f90 | 12 +++++- 3 files changed, 88 insertions(+), 18 deletions(-) diff --git a/src/nucfd_coeffs_mod.f90 b/src/nucfd_coeffs_mod.f90 index 05c374c..b5343de 100644 --- a/src/nucfd_coeffs_mod.f90 +++ b/src/nucfd_coeffs_mod.f90 @@ -31,14 +31,14 @@ module nucfd_coeffs contains - pure real function coeff_a_points(x) + 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(6, 4)), intent(in) :: x !! Stencil of points for the finite + type(nucfd_stencil_points), intent(in) :: x !! Stencil of points for the finite !! difference. - type(nucfd_stencil_deltas(5, 3)) :: h !! Stencil of grid spacings for the finite difference. - + type(nucfd_stencil_deltas) :: h ! Stencil of grid spacings for the finite difference. + h = points_to_deltas(x) coeff_a_points = coeff_a(h) @@ -48,10 +48,26 @@ pure 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(5, 3)), intent(in) :: h !! Stencil of grid spacings for the finite difference. + type(nucfd_stencil_deltas), intent(in) :: h !! Stencil of grid spacings for the finite difference. - coeff_a_deltas = 0.0 - + real :: hm2, hm1, h0, hp1, hp2 ! Grid deltas at i -2, -1, 0, +1, +2 + + hm2 = h%stencil(-2) + hm1 = h%stencil(-1) + h0 = h%stencil(0) + hp1 = h%stencil(1) + hp2 = h%stencil(2) + + associate(beta => alpha) ! To match Gamet et al. (1999) + coeff_a_deltas = hm1 * h0 * hp1 + h0**2 * hp1 + hm1 * h0 * hp2 + h0**2 * hp2 & + - hm1 * h0**2 * alpha - hm1 * h0 * hp1 * alpha & + - hm1 * h0 * hp2 * alpha - hm1 * h0 * hp1 * beta - h0**2 * hp1 * beta & + - hm1 * hp1**2 * beta - 2.0 * h0 * hp1**2 * beta - hp1**3 * beta & + + hm1 * h0 * hp2 * beta + h0**2 * hp2 * beta + 2.0 * hm1 * hp1 * hp2 * beta & + + 4.0 * h0 * hp1 * hp2 * beta + 3.0 * hp1**2 * hp2 * beta + coeff_a_deltas = coeff_a_deltas & + / (hp1 * (h0 + hp1) * (hm1 + h0 + hp1) * hp2) + end associate end function coeff_a_deltas end module nucfd_coeffs diff --git a/src/nucfd_types_mod.f90 b/src/nucfd_types_mod.f90 index 7a23460..6108aec 100644 --- a/src/nucfd_types_mod.f90 +++ b/src/nucfd_types_mod.f90 @@ -15,13 +15,11 @@ module nucfd_types implicit none - type nucfd_stencil(width, centre) + type nucfd_stencil !! Type representing a stencil. - integer, kind :: width = 1 !! Stencil width. - integer, kind :: centre = 1 !! Stencil centre. - real, dimension(1-centre:(1-centre)+(width-1)) :: stencil !! Stencil data. The lbound, ubound - !! are set to allow stencil access by - !! offsets. + real, dimension(:), allocatable :: stencil !! Stencil data + contains + final :: free_stencil end type nucfd_stencil type, extends(nucfd_stencil) :: nucfd_stencil_points @@ -35,12 +33,58 @@ module nucfd_types contains - pure function points_to_deltas(x) result(h) + subroutine create_stencil(width, centre, stencil) + !! Subroutine to create a stencil, specifying the width and centre. The underlying array is + !! allocated with bounds (1-centre):(1-centre)+(width-1) so that accesses are by offset from the + !! central position. + + integer, intent(in) :: width !! Specifies the width of the stencil + integer, intent(in) :: centre !! Specifies the index the stencil should be centred about + class(nucfd_stencil), intent(out) :: stencil !! The stencil object. + + integer :: lb, ub + + lb = 1 - centre + ub = (1 - centre) + (width - 1) + allocate(stencil%stencil(lb:ub)) + + end subroutine create_stencil + + function points_to_deltas(x) result(h) !! Helper function converting a stencil's point representation to a grid spacing representation. - type(nucfd_stencil_points(6, 4)), intent(in) :: x !! The stencil's point representation. - type(nucfd_stencil_deltas(5, 3)) :: h !! The stencil's grid spacing representation. + type(nucfd_stencil_points), intent(in) :: x !! The stencil's point representation. + type(nucfd_stencil_deltas) :: h !! The stencil's grid spacing representation. + + integer :: i + real :: xm1, x0 ! Start and end points of grid spacing at i. + integer :: ub, lb, width, centre ! Upper bound, lower bound, width and centre of grid stencil. + + lb = lbound(x%stencil, 1) + 1 + ub = ubound(x%stencil, 1) + width = (ub - lb) + 1 + centre = 1 - lb + call create_stencil(width, centre, h) + + do i = lbound(h%stencil, 1) + 1, ubound(h%stencil, 1) + xm1 = x%stencil(i - 1) + x0 = x%stencil(i) + + h%stencil(i) = x0 - xm1 + end do + end function points_to_deltas - + + subroutine free_stencil(stencil) + !! Frees a stencil object + + type(nucfd_stencil) :: stencil !! The stencil object to be freed. + + if (allocated(stencil%stencil)) then + deallocate(stencil%stencil) + end if + + end subroutine free_stencil + end module nucfd_types diff --git a/tests/fd-schemes/verify_coeffs.f90 b/tests/fd-schemes/verify_coeffs.f90 index 80d05bf..202388a 100644 --- a/tests/fd-schemes/verify_coeffs.f90 +++ b/tests/fd-schemes/verify_coeffs.f90 @@ -19,7 +19,7 @@ program verify_coeffs implicit none - type(nucfd_stencil_points(6, 4)) :: stencil ! Stencil of grid points. + type(nucfd_stencil_points) :: stencil ! Stencil of grid points. real :: n ! Mesh size real :: L ! Domain size @@ -33,6 +33,16 @@ program verify_coeffs n = 128 L = 1.0 h = L / real(n - 1) + + call create_stencil(6, 4, stencil) + + stencil%stencil(:) = 0.0 + stencil%stencil(-3) = -3.0 * h + stencil%stencil(-2) = -2.0 * h + stencil%stencil(-1) = -1.0 * h + stencil%stencil(0) = 0.0 + stencil%stencil(1) = +1.0 * h + stencil%stencil(2) = +2.0 * h a = coeff_a(stencil) call test_report("Coefficient A", check_scalar(a, aref / (2.0 * h))) From 94e09ef912fdbe285453b1c4d9d1cc2280161cd8 Mon Sep 17 00:00:00 2001 From: Paul Date: Sun, 26 Feb 2023 10:56:06 +0000 Subject: [PATCH 06/58] Implement B coefficient --- src/nucfd_coeffs_mod.f90 | 47 +++++++++++++++++++++++++++++- tests/fd-schemes/verify_coeffs.f90 | 13 +++++++-- 2 files changed, 57 insertions(+), 3 deletions(-) diff --git a/src/nucfd_coeffs_mod.f90 b/src/nucfd_coeffs_mod.f90 index b5343de..31894eb 100644 --- a/src/nucfd_coeffs_mod.f90 +++ b/src/nucfd_coeffs_mod.f90 @@ -20,11 +20,17 @@ module nucfd_coeffs private public :: coeff_a + public :: coeff_b interface coeff_a module procedure coeff_a_points module procedure coeff_a_deltas end interface coeff_a + + interface coeff_b + module procedure coeff_b_points + module procedure coeff_b_deltas + end interface coeff_b real, parameter, public :: alpha = 1.0 / 3.0 !! Off-diagonal coefficient for first derivative !! system. @@ -35,7 +41,7 @@ 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. + !! difference. type(nucfd_stencil_deltas) :: h ! Stencil of grid spacings for the finite difference. @@ -69,5 +75,44 @@ pure real function coeff_a_deltas(h) / (hp1 * (h0 + hp1) * (hm1 + h0 + hp1) * hp2) end associate end function coeff_a_deltas + + real function coeff_b_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. + + h = points_to_deltas(x) + coeff_b_points = coeff_b_deltas(h) + + end function coeff_b_points + + pure real function coeff_b_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 :: hm2, hm1, h0, hp1, hp2 ! Grid deltas at i -2, -1, 0, +1, +2 + + hm2 = h%stencil(-2) + hm1 = h%stencil(-1) + h0 = h%stencil(0) + hp1 = h%stencil(1) + hp2 = h%stencil(2) + + associate(beta => alpha) ! To match Gamet et al. (1999) + coeff_b_deltas = -hm1 * hp1**2 - h0 * hp1**2 - hm1 * hp1 * hp2 - h0 * hp1 * hp2 & + - 3.0 * hm1 * h0**2 * alpha - 4.0 * hm1 * h0 * hp1 * alpha & ! End line 1 + + h0**3 * alpha + 2.0 * h0**2 * hp1 * alpha - hm1 * hp1**2 * alpha & + + h0 * hp1**2 * alpha - 2.0 * hm1 * h0 * hp2 * alpha - hm1 * hp1 * hp2 * alpha & ! End line 2 + + h0**2 * hp2 * alpha + h0 * hp1 * hp2 * alpha + hm1 * hp1 * hp2 * beta & + + h0 * hp1 * hp2 * beta + hp1**2 * hp2 * beta ! End line 3 + coeff_b_deltas = coeff_b_deltas & + / (hm1 * h0 * (h0 + hp1) * (h0 + hp1 + hp2)) + end associate + end function coeff_b_deltas end module nucfd_coeffs diff --git a/tests/fd-schemes/verify_coeffs.f90 b/tests/fd-schemes/verify_coeffs.f90 index 202388a..fb694c8 100644 --- a/tests/fd-schemes/verify_coeffs.f90 +++ b/tests/fd-schemes/verify_coeffs.f90 @@ -25,8 +25,9 @@ program verify_coeffs real :: L ! Domain size real :: h ! Average grid spacing. - real :: a + real :: a, b real, parameter :: aref = 14.0 / 9.0 + real, parameter :: bref = -aref call initialise_suite("Verify coefficients") @@ -46,7 +47,15 @@ program verify_coeffs a = coeff_a(stencil) call test_report("Coefficient A", check_scalar(a, aref / (2.0 * h))) - + b = coeff_b(stencil) + call test_report("Coefficient B", check_scalar(b, bref / (2.0 * h))) + + stencil%stencil(-3) = -5.0 * h + stencil%stencil(2) = +3.0 * h + a = coeff_a(stencil) + b = coeff_b(stencil) + call test_report("Coefficient A /= B", .not. check_scalar(a, b)) + call finalise_suite() end program verify_coeffs From cac6613351e21ff8e46310c6495762928726b815 Mon Sep 17 00:00:00 2001 From: Paul Date: Sun, 26 Feb 2023 11:01:16 +0000 Subject: [PATCH 07/58] Implement coefficient C --- src/nucfd_coeffs_mod.f90 | 42 ++++++++++++++++++++++++++++++ tests/fd-schemes/verify_coeffs.f90 | 7 +++-- 2 files changed, 47 insertions(+), 2 deletions(-) diff --git a/src/nucfd_coeffs_mod.f90 b/src/nucfd_coeffs_mod.f90 index 31894eb..5f3d31b 100644 --- a/src/nucfd_coeffs_mod.f90 +++ b/src/nucfd_coeffs_mod.f90 @@ -21,6 +21,7 @@ module nucfd_coeffs private public :: coeff_a public :: coeff_b + public :: coeff_c interface coeff_a module procedure coeff_a_points @@ -31,6 +32,11 @@ module nucfd_coeffs module procedure coeff_b_points module procedure coeff_b_deltas end interface coeff_b + + interface coeff_c + module procedure coeff_c_points + module procedure coeff_c_deltas + end interface coeff_c real, parameter, public :: alpha = 1.0 / 3.0 !! Off-diagonal coefficient for first derivative !! system. @@ -114,5 +120,41 @@ pure real function coeff_b_deltas(h) / (hm1 * h0 * (h0 + hp1) * (h0 + hp1 + hp2)) end associate end function coeff_b_deltas + + real function coeff_c_points(x) + !! Compute the coefficient acting on f_{i+2} 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. + + h = points_to_deltas(x) + coeff_c_points = coeff_c_deltas(h) + + end function coeff_c_points + + pure real function coeff_c_deltas(h) + !! Compute the coefficient acting on f_{i+21} 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 :: hm2, hm1, h0, hp1, hp2 ! Grid deltas at i -2, -1, 0, +1, +2 + + hm2 = h%stencil(-2) + hm1 = h%stencil(-1) + h0 = h%stencil(0) + hp1 = h%stencil(1) + hp2 = h%stencil(2) + + associate(beta => alpha) ! To match Gamet et al. (1999) + coeff_c_deltas = -hm1 * h0 * hp1 - h0**2 * hp1 + hm1 * h0**2 * alpha & + + hm1 * h0 * hp1 * alpha + hm1 * h0 * hp1 * beta + h0**2 * hp1 * beta & ! End line 1 + + hm1 * hp1**2 * beta + 2.0 * h0 * hp1**2 * beta + hp1**3 * beta ! End line 2 + coeff_c_deltas = coeff_c_deltas & + / (hp2 * (hp1 + hp2) * (h0 + hp1 + hp2) * (hm1 + h0 + hp1 + hp2)) + end associate + end function coeff_c_deltas end module nucfd_coeffs diff --git a/tests/fd-schemes/verify_coeffs.f90 b/tests/fd-schemes/verify_coeffs.f90 index fb694c8..328da1d 100644 --- a/tests/fd-schemes/verify_coeffs.f90 +++ b/tests/fd-schemes/verify_coeffs.f90 @@ -25,9 +25,10 @@ program verify_coeffs real :: L ! Domain size real :: h ! Average grid spacing. - real :: a, b + real :: a, b, c real, parameter :: aref = 14.0 / 9.0 real, parameter :: bref = -aref + real, parameter :: cref = 1.0 / 9.0 call initialise_suite("Verify coefficients") @@ -49,7 +50,9 @@ program verify_coeffs call test_report("Coefficient A", check_scalar(a, aref / (2.0 * h))) b = coeff_b(stencil) call test_report("Coefficient B", check_scalar(b, bref / (2.0 * h))) - + c = coeff_c(stencil) + call test_report("Coefficient C", check_scalar(c, cref / (4.0 * h))) + stencil%stencil(-3) = -5.0 * h stencil%stencil(2) = +3.0 * h a = coeff_a(stencil) From 8034c7d511eb32eab3e824adc061daa71f3650d3 Mon Sep 17 00:00:00 2001 From: Paul Date: Sun, 26 Feb 2023 11:06:24 +0000 Subject: [PATCH 08/58] Implement coefficient D --- src/nucfd_coeffs_mod.f90 | 42 ++++++++++++++++++++++++++++++ tests/fd-schemes/verify_coeffs.f90 | 9 ++++++- 2 files changed, 50 insertions(+), 1 deletion(-) diff --git a/src/nucfd_coeffs_mod.f90 b/src/nucfd_coeffs_mod.f90 index 5f3d31b..2fda432 100644 --- a/src/nucfd_coeffs_mod.f90 +++ b/src/nucfd_coeffs_mod.f90 @@ -22,6 +22,7 @@ module nucfd_coeffs public :: coeff_a public :: coeff_b public :: coeff_c + public :: coeff_d interface coeff_a module procedure coeff_a_points @@ -37,6 +38,11 @@ module nucfd_coeffs module procedure coeff_c_points module procedure coeff_c_deltas end interface coeff_c + + interface coeff_d + module procedure coeff_d_points + module procedure coeff_d_deltas + end interface coeff_d real, parameter, public :: alpha = 1.0 / 3.0 !! Off-diagonal coefficient for first derivative !! system. @@ -156,5 +162,41 @@ pure real function coeff_c_deltas(h) / (hp2 * (hp1 + hp2) * (h0 + hp1 + hp2) * (hm1 + h0 + hp1 + hp2)) end associate end function coeff_c_deltas + + real function coeff_d_points(x) + !! Compute the coefficient acting on f_{i+2} 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. + + h = points_to_deltas(x) + coeff_d_points = coeff_d_deltas(h) + + end function coeff_d_points + + pure real function coeff_d_deltas(h) + !! Compute the coefficient acting on f_{i+21} 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 :: hm2, hm1, h0, hp1, hp2 ! Grid deltas at i -2, -1, 0, +1, +2 + + hm2 = h%stencil(-2) + hm1 = h%stencil(-1) + h0 = h%stencil(0) + hp1 = h%stencil(1) + hp2 = h%stencil(2) + + associate(beta => alpha) ! To match Gamet et al. (1999) + coeff_d_deltas = h0 * hp1**2 + h0 * hp1 * hp2 - h0**3 * alpha - 2.0 * h0**2 * hp1 * alpha & + - h0 * hp1**2 * alpha - h0**2 * hp2 * alpha - h0 * hp1 * hp2 * alpha & ! End line 1 + - h0 * hp1 * hp2 * beta - hp1**2 * hp2 * beta ! End line 2 + coeff_d_deltas = coeff_d_deltas & + / (hm1 * (hm1 + h0) * (hm1 + h0 + hp1) * (hm1 + h0 + hp1 + hp2)) + end associate + end function coeff_d_deltas end module nucfd_coeffs diff --git a/tests/fd-schemes/verify_coeffs.f90 b/tests/fd-schemes/verify_coeffs.f90 index 328da1d..e53aa30 100644 --- a/tests/fd-schemes/verify_coeffs.f90 +++ b/tests/fd-schemes/verify_coeffs.f90 @@ -25,10 +25,12 @@ program verify_coeffs real :: L ! Domain size real :: h ! Average grid spacing. - real :: a, b, c + real :: a, b, c, d, e real, parameter :: aref = 14.0 / 9.0 real, parameter :: bref = -aref real, parameter :: cref = 1.0 / 9.0 + real, parameter :: dref = -cref + real, parameter :: eref = 0.0 call initialise_suite("Verify coefficients") @@ -52,12 +54,17 @@ program verify_coeffs call test_report("Coefficient B", check_scalar(b, bref / (2.0 * h))) c = coeff_c(stencil) call test_report("Coefficient C", check_scalar(c, cref / (4.0 * h))) + d = coeff_d(stencil) + call test_report("Coefficient D", check_scalar(d, dref / (4.0 * h))) stencil%stencil(-3) = -5.0 * h stencil%stencil(2) = +3.0 * h a = coeff_a(stencil) b = coeff_b(stencil) call test_report("Coefficient A /= B", .not. check_scalar(a, b)) + c = coeff_c(stencil) + d = coeff_d(stencil) + call test_report("Coefficient C /= D", .not. check_scalar(c, d)) call finalise_suite() From 9852f7e8cf537c03abd719e6a8d49bf9ebd2b21a Mon Sep 17 00:00:00 2001 From: Paul Date: Sun, 26 Feb 2023 11:09:36 +0000 Subject: [PATCH 09/58] Implement coefficient E --- src/nucfd_coeffs_mod.f90 | 25 +++++++++++++++++++++---- tests/fd-schemes/verify_coeffs.f90 | 2 ++ 2 files changed, 23 insertions(+), 4 deletions(-) diff --git a/src/nucfd_coeffs_mod.f90 b/src/nucfd_coeffs_mod.f90 index 2fda432..540c674 100644 --- a/src/nucfd_coeffs_mod.f90 +++ b/src/nucfd_coeffs_mod.f90 @@ -23,6 +23,7 @@ module nucfd_coeffs public :: coeff_b public :: coeff_c public :: coeff_d + public :: coeff_e interface coeff_a module procedure coeff_a_points @@ -141,7 +142,7 @@ real function coeff_c_points(x) end function coeff_c_points pure real function coeff_c_deltas(h) - !! Compute the coefficient acting on f_{i+21} of the finite diference given a stencil of grid + !! Compute the coefficient acting on f_{i+2} 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. @@ -164,7 +165,7 @@ pure real function coeff_c_deltas(h) end function coeff_c_deltas real function coeff_d_points(x) - !! Compute the coefficient acting on f_{i+2} of the finite diference given a stencil of points. + !! Compute the coefficient acting on f_{i-2} of the finite diference given a stencil of points. type(nucfd_stencil_points), intent(in) :: x !! Stencil of points for the finite !! difference. @@ -177,7 +178,7 @@ real function coeff_d_points(x) end function coeff_d_points pure real function coeff_d_deltas(h) - !! Compute the coefficient acting on f_{i+21} of the finite diference given a stencil of grid + !! Compute the coefficient acting on f_{i-2} 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. @@ -198,5 +199,21 @@ pure real function coeff_d_deltas(h) / (hm1 * (hm1 + h0) * (hm1 + h0 + hp1) * (hm1 + h0 + hp1 + hp2)) end associate end function coeff_d_deltas - + + real function coeff_e(x) + !! Compute the coefficient acting on f_{i} of the finite diference given a stencil of points. + + type(nucfd_stencil_points), intent(in) :: x !! Stencil of points for the finite + !! difference. + + real :: a, b, c, d ! The neighbouring coefficients + + a = coeff_a(x) + b = coeff_b(x) + c = coeff_c(x) + d = coeff_d(x) + + coeff_e = -(a + b + c + d) + + end function coeff_e end module nucfd_coeffs diff --git a/tests/fd-schemes/verify_coeffs.f90 b/tests/fd-schemes/verify_coeffs.f90 index e53aa30..3727005 100644 --- a/tests/fd-schemes/verify_coeffs.f90 +++ b/tests/fd-schemes/verify_coeffs.f90 @@ -56,6 +56,8 @@ program verify_coeffs call test_report("Coefficient C", check_scalar(c, cref / (4.0 * h))) d = coeff_d(stencil) call test_report("Coefficient D", check_scalar(d, dref / (4.0 * h))) + e = coeff_e(stencil) + call test_report("Coefficient E", check_scalar(e, eref)) stencil%stencil(-3) = -5.0 * h stencil%stencil(2) = +3.0 * h From 3a12b0fe91917a0e256eb1e33d8ed3754d0148c5 Mon Sep 17 00:00:00 2001 From: Paul Date: Sun, 26 Feb 2023 14:46:20 +0000 Subject: [PATCH 10/58] Update comments to improve FORD documentation --- src/nucfd_coeffs_mod.f90 | 4 ++++ src/nucfd_types_mod.f90 | 16 +++------------- tests/fd-schemes/verify_coeffs.f90 | 17 +++++------------ tests/nucfd_tests_mod.f90 | 10 +++++----- 4 files changed, 17 insertions(+), 30 deletions(-) diff --git a/src/nucfd_coeffs_mod.f90 b/src/nucfd_coeffs_mod.f90 index a8e54c4..735547f 100644 --- a/src/nucfd_coeffs_mod.f90 +++ b/src/nucfd_coeffs_mod.f90 @@ -15,21 +15,25 @@ module nucfd_coeffs public :: coeff_e interface coeff_a + !! Compute the coefficient acting on f_{i+1} of the finite diference stencil. module procedure coeff_a_points module procedure coeff_a_deltas end interface coeff_a interface coeff_b + !! Compute the coefficient acting on f_{i-1} of the finite diference stencil. module procedure coeff_b_points module procedure coeff_b_deltas end interface coeff_b interface coeff_c + !! Compute the coefficient acting on f_{i+2} of the finite diference stencil. module procedure coeff_c_points module procedure coeff_c_deltas end interface coeff_c interface coeff_d + !! Compute the coefficient acting on f_{i-2} of the finite diference stencil. module procedure coeff_d_points module procedure coeff_d_deltas end interface coeff_d diff --git a/src/nucfd_types_mod.f90 b/src/nucfd_types_mod.f90 index 6108aec..bd88a38 100644 --- a/src/nucfd_types_mod.f90 +++ b/src/nucfd_types_mod.f90 @@ -1,17 +1,7 @@ -!!!! src/nucfd_types_mod.f90 -!!! -!!!! Description -!!! -!!! Defines the custom types used by NuCFD. -!!! -!!! Provides the nucfd_types module. -!!! -!!!! LICENSE -!!! -!!! SPDX-License-Identifier: BSD-3-Clause -!!! - module nucfd_types + !! Module defining the custom types used by NuCFD. + !! + !! SPDX-License-Identifier: BSD-3-Clause implicit none diff --git a/tests/fd-schemes/verify_coeffs.f90 b/tests/fd-schemes/verify_coeffs.f90 index 3727005..11343dd 100644 --- a/tests/fd-schemes/verify_coeffs.f90 +++ b/tests/fd-schemes/verify_coeffs.f90 @@ -1,16 +1,9 @@ -!!!! tests/fd-schemes/verify_coeffs.f90 -!!! -!!!! Description -!!! -!!! Part of the fd-schemes test suite. -!!! Tests the computation of finite difference coefficients for non-uniform grids. -!!! -!!!! LICENSE -!!! -!!! SPDX-License-Identifier: BSD-3-Clause -!!! - program verify_coeffs + !! Tests the computation of finite difference coefficients for non-uniform grids. + !! + !! Part of the fd-schemes test suite. + !! + !! SPDX-License-Identifier: BSD-3-Clause use nucfd_types use nucfd_coeffs diff --git a/tests/nucfd_tests_mod.f90 b/tests/nucfd_tests_mod.f90 index 2823e83..85ba249 100644 --- a/tests/nucfd_tests_mod.f90 +++ b/tests/nucfd_tests_mod.f90 @@ -69,10 +69,10 @@ subroutine test_report(test_name, test_status) end subroutine test_report logical function check_scalar(test, ref) - ! Compute error and report errors. + !! Compute error and report for scalar values. - real, intent(in) :: test ! The test data - real, intent(in) :: ref ! The reference data + real, intent(in) :: test !! The test data + real, intent(in) :: ref !! The reference data real :: err logical :: test_passing @@ -93,8 +93,8 @@ end function check_scalar logical function check_rms(test, ref) !! Compute RMS of error and report errors. - real, dimension(:), intent(in) :: test ! The test data - real, dimension(:), intent(in) :: ref ! The reference data + real, dimension(:), intent(in) :: test !! The test data + real, dimension(:), intent(in) :: ref !! The reference data integer :: n integer :: i From bec3ac8415b014a0a5fa8fda9eac8b095cdb038d Mon Sep 17 00:00:00 2001 From: Paul Date: Sun, 26 Feb 2023 17:05:11 +0000 Subject: [PATCH 11/58] Add index stencil type --- src/nucfd_coeffs_mod.f90 | 60 ++++++++++++++++++++---------- src/nucfd_types_mod.f90 | 44 ++++++++++++++++++---- tests/fd-schemes/verify_coeffs.f90 | 44 ++++++++++++++-------- 3 files changed, 104 insertions(+), 44 deletions(-) diff --git a/src/nucfd_coeffs_mod.f90 b/src/nucfd_coeffs_mod.f90 index 735547f..3c8602d 100644 --- a/src/nucfd_coeffs_mod.f90 +++ b/src/nucfd_coeffs_mod.f90 @@ -64,11 +64,16 @@ pure real function coeff_a_deltas(h) real :: hm2, hm1, h0, hp1, hp2 ! Grid deltas at i -2, -1, 0, +1, +2 - hm2 = h%stencil(-2) - hm1 = h%stencil(-1) - h0 = h%stencil(0) - hp1 = h%stencil(1) - hp2 = h%stencil(2) + select type(deltas => h%stencil) + type is(real) + hm2 = deltas(-2) + 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) coeff_a_deltas = hm1 * h0 * hp1 + h0**2 * hp1 + hm1 * h0 * hp2 + h0**2 * hp2 & @@ -103,11 +108,16 @@ pure real function coeff_b_deltas(h) real :: hm2, hm1, h0, hp1, hp2 ! Grid deltas at i -2, -1, 0, +1, +2 - hm2 = h%stencil(-2) - hm1 = h%stencil(-1) - h0 = h%stencil(0) - hp1 = h%stencil(1) - hp2 = h%stencil(2) + select type(deltas => h%stencil) + type is(real) + hm2 = deltas(-2) + 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) coeff_b_deltas = -hm1 * hp1**2 - h0 * hp1**2 - hm1 * hp1 * hp2 - h0 * hp1 * hp2 & @@ -142,11 +152,16 @@ pure real function coeff_c_deltas(h) real :: hm2, hm1, h0, hp1, hp2 ! Grid deltas at i -2, -1, 0, +1, +2 - hm2 = h%stencil(-2) - hm1 = h%stencil(-1) - h0 = h%stencil(0) - hp1 = h%stencil(1) - hp2 = h%stencil(2) + select type(deltas => h%stencil) + type is(real) + hm2 = deltas(-2) + 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) coeff_c_deltas = -hm1 * h0 * hp1 - h0**2 * hp1 + hm1 * h0**2 * alpha & @@ -178,11 +193,16 @@ pure real function coeff_d_deltas(h) real :: hm2, hm1, h0, hp1, hp2 ! Grid deltas at i -2, -1, 0, +1, +2 - hm2 = h%stencil(-2) - hm1 = h%stencil(-1) - h0 = h%stencil(0) - hp1 = h%stencil(1) - hp2 = h%stencil(2) + select type(deltas => h%stencil) + type is(real) + hm2 = deltas(-2) + 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) coeff_d_deltas = h0 * hp1**2 + h0 * hp1 * hp2 - h0**3 * alpha - 2.0 * h0**2 * hp1 * alpha & diff --git a/src/nucfd_types_mod.f90 b/src/nucfd_types_mod.f90 index bd88a38..6e3f17a 100644 --- a/src/nucfd_types_mod.f90 +++ b/src/nucfd_types_mod.f90 @@ -7,7 +7,7 @@ module nucfd_types type nucfd_stencil !! Type representing a stencil. - real, dimension(:), allocatable :: stencil !! Stencil data + class(*), dimension(:), allocatable :: stencil !! Stencil data contains final :: free_stencil end type nucfd_stencil @@ -21,6 +21,10 @@ module nucfd_types !! defined as h_i = x_i - x_{i-1}. end type nucfd_stencil_deltas + type, extends(nucfd_stencil) :: nucfd_index_stencil + !! Type representing a stencil's indices. + end type nucfd_index_stencil + contains subroutine create_stencil(width, centre, stencil) @@ -36,7 +40,17 @@ subroutine create_stencil(width, centre, stencil) lb = 1 - centre ub = (1 - centre) + (width - 1) - allocate(stencil%stencil(lb:ub)) + select type(stencil) + type is(nucfd_stencil_points) + allocate(real::stencil%stencil(lb:ub)) + type is(nucfd_stencil_deltas) + allocate(real::stencil%stencil(lb:ub)) + type is(nucfd_index_stencil) + allocate(integer::stencil%stencil(lb:ub)) + class default + print *, "Error: trying to create unknown stencil type" + error stop + end select end subroutine create_stencil @@ -56,13 +70,27 @@ function points_to_deltas(x) result(h) width = (ub - lb) + 1 centre = 1 - lb call create_stencil(width, centre, h) - - do i = lbound(h%stencil, 1) + 1, ubound(h%stencil, 1) - xm1 = x%stencil(i - 1) - x0 = x%stencil(i) - h%stencil(i) = x0 - xm1 - end do + associate(points => x%stencil, & + deltas => h%stencil) + select type(points) + type is(real) + select type(deltas) + type is(real) + do i = lbound(h%stencil, 1) + 1, ubound(h%stencil, 1) + xm1 = points(i - 1) + x0 = points(i) + + deltas(i) = x0 - xm1 + end do + class default + print *, "Error: Difference stencil is misallocated!" + end select + class default + print *, "Error: Coordinate stencil is misallocated!" + error stop + end select + end associate end function points_to_deltas diff --git a/tests/fd-schemes/verify_coeffs.f90 b/tests/fd-schemes/verify_coeffs.f90 index 11343dd..9ccb8cc 100644 --- a/tests/fd-schemes/verify_coeffs.f90 +++ b/tests/fd-schemes/verify_coeffs.f90 @@ -1,9 +1,11 @@ +! tests/fd-schemes/verify_coeffs.f90 +! +!! Part of the fd-schemes test suite. +! +! SPDX-License-Identifier: BSD-3-Clause + program verify_coeffs !! Tests the computation of finite difference coefficients for non-uniform grids. - !! - !! Part of the fd-schemes test suite. - !! - !! SPDX-License-Identifier: BSD-3-Clause use nucfd_types use nucfd_coeffs @@ -11,12 +13,12 @@ program verify_coeffs use nucfd_tests implicit none - - type(nucfd_stencil_points) :: stencil ! Stencil of grid points. real :: n ! Mesh size real :: L ! Domain size real :: h ! Average grid spacing. + + type(nucfd_stencil_points) :: stencil ! Stencil of grid points. real :: a, b, c, d, e real, parameter :: aref = 14.0 / 9.0 @@ -32,14 +34,18 @@ program verify_coeffs h = L / real(n - 1) call create_stencil(6, 4, stencil) - - stencil%stencil(:) = 0.0 - stencil%stencil(-3) = -3.0 * h - stencil%stencil(-2) = -2.0 * h - stencil%stencil(-1) = -1.0 * h - stencil%stencil(0) = 0.0 - stencil%stencil(1) = +1.0 * h - stencil%stencil(2) = +2.0 * h + select type(points => stencil%stencil) + type is(real) + points(:) = 0.0 + points(-3) = -3.0 * h + points(-2) = -2.0 * h + points(-1) = -1.0 * h + points(0) = 0.0 + points(1) = +1.0 * h + points(2) = +2.0 * h + class default + print *, "Error: Coordinate stencil is misallocated!" + end select a = coeff_a(stencil) call test_report("Coefficient A", check_scalar(a, aref / (2.0 * h))) @@ -52,8 +58,14 @@ program verify_coeffs e = coeff_e(stencil) call test_report("Coefficient E", check_scalar(e, eref)) - stencil%stencil(-3) = -5.0 * h - stencil%stencil(2) = +3.0 * h + select type(points => stencil%stencil) + type is(real) + points(-3) = -5.0 * h + points(2) = +3.0 * h + class default + print *, "Error: Coordinate stencil is misallocated!" + end select + a = coeff_a(stencil) b = coeff_b(stencil) call test_report("Coefficient A /= B", .not. check_scalar(a, b)) From 3b2e1e9511e97ee87771abaed925e6ee1227bf7f Mon Sep 17 00:00:00 2001 From: Paul Date: Sun, 26 Feb 2023 17:42:55 +0000 Subject: [PATCH 12/58] Initial implementation of RHS derivative function Adds test that derivative RHS obeys law of derivatives when function is shifted by a constant --- CMakeLists.txt | 6 +- src/nucfd_deriv_mod.f90 | 63 ++++++++++++++++++++ tests/fd-schemes/differentiation_rules.f90 | 69 ++++++++++++++++++++++ 3 files changed, 136 insertions(+), 2 deletions(-) create mode 100644 src/nucfd_deriv_mod.f90 create mode 100644 tests/fd-schemes/differentiation_rules.f90 diff --git a/CMakeLists.txt b/CMakeLists.txt index 9b89256..e1bb2f1 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -30,13 +30,14 @@ target_compile_options(nucfd_debug_flags INTERFACE "$<${icc_like_fc}:-g;-O0;-check;all;-fpe0;-traceback;-init:snan,arrays>") add_library(nucfd_strict_flags INTERFACE) target_compile_options(nucfd_strict_flags INTERFACE - "$<${gcc_like_fc}:-Werror>") + "$<${gcc_like_fc}:-Werror;-Wno-error=integer-division>") ## Build nucfd library add_library(nucfd src/nucfd_trid_solver_mod.f90 src/nucfd_types_mod.f90 - src/nucfd_coeffs_mod.f90) + src/nucfd_coeffs_mod.f90 + src/nucfd_deriv_mod.f90) target_link_libraries(nucfd nucfd_compiler_flags) if (CMAKE_BUILD_TYPE MATCHES "Debug") target_link_libraries(nucfd nucfd_debug_flags) @@ -80,6 +81,7 @@ add_test_libs(tridsolver system_00_anti-symm trid_test_utils) # Finite difference test suite define_test(fd-schemes verify_coeffs) +define_test(fd-schemes differentiation_rules) ## Documentation add_custom_target(doc ford diff --git a/src/nucfd_deriv_mod.f90 b/src/nucfd_deriv_mod.f90 new file mode 100644 index 0000000..b5b11af --- /dev/null +++ b/src/nucfd_deriv_mod.f90 @@ -0,0 +1,63 @@ +! src/nucfd_deriv_mod.f90 +! +!! Defines the nucfd_deriv module. +! +! SPDX-License-Identifier: BSD-3-Clause + +module nucfd_deriv + !! Module defining the coefficients for compact finite difference schemes on non-uniform grids. + + use nucfd_types + use nucfd_coeffs + + implicit none + + private + public :: deriv_rhs + +contains + + subroutine deriv_rhs(f, stencil, x, dfdx) + + real, dimension(:), intent(in) :: f + type(nucfd_index_stencil), intent(in) :: stencil + real, dimension(:), intent(in) :: x + real, intent(out) :: dfdx + + type(nucfd_stencil_points) :: stencil_coordinates + + real :: a, b, c, d, e + + call create_stencil(6, 4, stencil_coordinates) + + select type(indices => stencil%stencil) + type is(integer) + select type(points => stencil_coordinates%stencil) + type is(real) + points(-3) = x(indices(-2) - 1) + points(-2) = x(indices(-2)) + points(-1) = x(indices(-1)) + points(+0) = x(indices(+0)) + points(+1) = x(indices(+1)) + points(+2) = x(indices(+2)) + class default + print *, "Error: Coordinate stencil is misallocated!" + error stop + end select + + a = coeff_a(stencil_coordinates) + b = coeff_b(stencil_coordinates) + c = coeff_c(stencil_coordinates) + d = coeff_d(stencil_coordinates) + e = coeff_e(stencil_coordinates) + dfdx = (a * f(indices(+1)) + b * f(indices(-1))) & + + (c * f(indices(+2)) + d * f(indices(-2))) & + + e * f(indices(0)) + class default + print *, "Error: Index stencil is misallocated!" + error stop + end select + + end subroutine deriv_rhs + +end module nucfd_deriv diff --git a/tests/fd-schemes/differentiation_rules.f90 b/tests/fd-schemes/differentiation_rules.f90 new file mode 100644 index 0000000..5cac99b --- /dev/null +++ b/tests/fd-schemes/differentiation_rules.f90 @@ -0,0 +1,69 @@ +! tests/fd-schemes/differentiation_rules.f90 +! +!! Part of the fd-schemes test suite. +! +! SPDX-License-Identifier: BSD-3-Clause + +program differentiation_rules + !! Tests that rules of differentiation are respected. + + use nucfd_types + use nucfd_deriv + + use nucfd_tests + + implicit none + + integer :: n ! Grid size + real :: L ! Domain size + real :: h ! Grid spacing + + real, dimension(:), allocatable :: x ! Grid + real, dimension(:), allocatable :: f ! Function + real :: dfdx ! Derivative + real :: dgdx ! Derivative + + ! Stencil + type(nucfd_index_stencil) :: stencil + integer :: i + integer, parameter :: width = 5 + integer, parameter :: centre = 3 + + call initialise_suite("Differentiation rules") + + n = 33 + L = 1.0 + h = L / real(n - 1) + + allocate(x(n)) + allocate(f(n)) + do i = 1, n + x(i) = real(i - 1) * h + f(i) = sin(x(i)) + end do + + call create_stencil(width, centre, stencil) + i = 33 / 2 + select type(indices => stencil%stencil) + type is(integer) + indices(-2) = i - 2 + indices(-1) = i - 1 + indices(+0) = i + 0 + indices(+1) = i + 1 + indices(+2) = i + 2 + class default + print *, "Error: Index stencil is misallocated!" + error stop + end select + + call deriv_rhs(f, stencil, x, dfdx) + call deriv_rhs(f + 1.0, stencil, x, dgdx) + + call test_report("Shifted derivative", check_scalar(dgdx, dfdx)) + + deallocate(x) + deallocate(f) + + call finalise_suite() + +end program differentiation_rules From 9e6ce2f79909c1823151627e95fa68e2b291b37e Mon Sep 17 00:00:00 2001 From: Paul Bartholomew Date: Sun, 26 Feb 2023 19:07:47 +0000 Subject: [PATCH 13/58] Move to f2018 standard Calling error stop from PURE functions is only supported in f2018. --- CMakeLists.txt | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index e1bb2f1..00d704b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -22,8 +22,8 @@ set(gcc_like_fc "$") set(icc_like_fc "$") 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>") + "$<${gcc_like_fc}:-std=f2018;-Wall;-Wextra;-Wpedantic;-Warray-bounds;-Wimplicit-interface;-Wimplicit-procedure;-fimplicit-none>" + "$<${icc_like_fc}:-std18;-warn;all>") 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>" From bfe44c059fa0b4fca3bbddbee4e521d94981998e Mon Sep 17 00:00:00 2001 From: Paul Bartholomew Date: Sun, 26 Feb 2023 21:02:21 +0000 Subject: [PATCH 14/58] Need to disable FP optimisation in Cray compiler Otherwise it fails the tests --- CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 391af6d..ed9c2ee 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -25,7 +25,7 @@ add_library(nucfd_compiler_flags INTERFACE) target_compile_options(nucfd_compiler_flags INTERFACE "$<${gcc_like_fc}:-std=f2018;-Wall;-Wextra;-Wpedantic;-Warray-bounds;-Wimplicit-interface;-Wimplicit-procedure;-fimplicit-none>" "$<${icc_like_fc}:-std18;-warn;all>" - "$<${cray_like_fc}:-en;-emf;-eI>") + "$<${cray_like_fc}:-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>" From db9a91ab6087b647237d0654fe2322ded66c8095 Mon Sep 17 00:00:00 2001 From: Paul Date: Mon, 27 Feb 2023 23:16:37 +0000 Subject: [PATCH 15/58] WIP Implement coeff_a_alt Currently missing terms that cancel, but passing coefficient verification with optimisation --- src/nucfd_coeffs_mod.f90 | 51 +++++++++++++++++++++++++++++++++++++++- 1 file changed, 50 insertions(+), 1 deletion(-) diff --git a/src/nucfd_coeffs_mod.f90 b/src/nucfd_coeffs_mod.f90 index 3c8602d..4c7096d 100644 --- a/src/nucfd_coeffs_mod.f90 +++ b/src/nucfd_coeffs_mod.f90 @@ -8,7 +8,7 @@ module nucfd_coeffs implicit none private - public :: coeff_a + public :: coeff_a, coeff_a_alt public :: coeff_b public :: coeff_c public :: coeff_d @@ -20,6 +20,12 @@ module nucfd_coeffs module procedure coeff_a_deltas end interface coeff_a + interface coeff_a_alt + !! Compute the coefficient acting on f_{i+1} of the finite diference stencil. + module procedure coeff_a_points_alt + module procedure coeff_a_deltas_alt + end interface coeff_a_alt + interface coeff_b !! Compute the coefficient acting on f_{i-1} of the finite diference stencil. module procedure coeff_b_points @@ -87,6 +93,49 @@ pure real function coeff_a_deltas(h) end associate end function coeff_a_deltas + real function coeff_a_points_alt(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. + + h = points_to_deltas(x) + coeff_a_points_alt = coeff_a_alt(h) + + end function coeff_a_points_alt + + real function coeff_a_deltas_alt(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 :: hm2, hm1, h0, hp1, hp2 ! Grid deltas at i -2, -1, 0, +1, +2 + + select type(deltas => h%stencil) + type is(real) + hm2 = deltas(-2) + 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) + coeff_a_deltas_alt = h0 * ((hm1 + h0) * (hp1 + hp2) + 2.0 * hp1 * hp2 * beta) ! = (14/3) h^3 + print *, "+++", coeff_a_deltas_alt / (h0**3), 14.0 / 3.0, "+++" + coeff_a_deltas_alt = coeff_a_deltas_alt & + / (3.0 * hp1 * ((h0 + hp1) / 2.0) * hp2) ! => 14.0 / 9.0 when h = const + print *, "+++", coeff_a_deltas_alt, 14.0 / 9.0, "+++" + coeff_a_deltas_alt = coeff_a_deltas_alt & + / (2.0 * ((hm1 + h0 + hp1) / 3.0)) ! => (14.0 / 9.0) / (2h) when h = const + end associate + end function coeff_a_deltas_alt + real function coeff_b_points(x) !! Compute the coefficient acting on f_{i-1} of the finite diference given a stencil of points. From 01f9e1ecc25bba8e12757546368a339601afb77d Mon Sep 17 00:00:00 2001 From: Paul Date: Mon, 27 Feb 2023 23:20:03 +0000 Subject: [PATCH 16/58] WIP coeff_a_alt adding additional terms --- src/nucfd_coeffs_mod.f90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/nucfd_coeffs_mod.f90 b/src/nucfd_coeffs_mod.f90 index 4c7096d..3370d60 100644 --- a/src/nucfd_coeffs_mod.f90 +++ b/src/nucfd_coeffs_mod.f90 @@ -126,7 +126,8 @@ real function coeff_a_deltas_alt(h) end select associate(beta => alpha) ! To match Gamet et al. (1999) - coeff_a_deltas_alt = h0 * ((hm1 + h0) * (hp1 + hp2) + 2.0 * hp1 * hp2 * beta) ! = (14/3) h^3 + coeff_a_deltas_alt = h0 * ((hm1 + h0) * (hp1 + hp2) + 2.0 * hp1 * hp2 * beta) & ! = (14/3) h^3 + + hm1 * h0 * hp2 * (beta - alpha) ! Should cancel for case alpha = beta print *, "+++", coeff_a_deltas_alt / (h0**3), 14.0 / 3.0, "+++" coeff_a_deltas_alt = coeff_a_deltas_alt & / (3.0 * hp1 * ((h0 + hp1) / 2.0) * hp2) ! => 14.0 / 9.0 when h = const From d9e7c6b6ad8ff30bb5feab7c846cfa64a6985640 Mon Sep 17 00:00:00 2001 From: Paul Date: Mon, 27 Feb 2023 23:21:28 +0000 Subject: [PATCH 17/58] WIP coeff_a_alt add additional terms --- src/nucfd_coeffs_mod.f90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/nucfd_coeffs_mod.f90 b/src/nucfd_coeffs_mod.f90 index 3370d60..0828b8d 100644 --- a/src/nucfd_coeffs_mod.f90 +++ b/src/nucfd_coeffs_mod.f90 @@ -127,7 +127,8 @@ real function coeff_a_deltas_alt(h) associate(beta => alpha) ! To match Gamet et al. (1999) coeff_a_deltas_alt = h0 * ((hm1 + h0) * (hp1 + hp2) + 2.0 * hp1 * hp2 * beta) & ! = (14/3) h^3 - + hm1 * h0 * hp2 * (beta - alpha) ! Should cancel for case alpha = beta + + hm1 * h0 * hp2 * (beta - alpha) & ! Should cancel for case alpha = beta + + (h0**2) * (hp2 - hp1) * beta ! Should cancel for constant h print *, "+++", coeff_a_deltas_alt / (h0**3), 14.0 / 3.0, "+++" coeff_a_deltas_alt = coeff_a_deltas_alt & / (3.0 * hp1 * ((h0 + hp1) / 2.0) * hp2) ! => 14.0 / 9.0 when h = const From f5a16d68c505abdb8f740ae8915b08be6549c748 Mon Sep 17 00:00:00 2001 From: Paul Date: Mon, 27 Feb 2023 23:24:32 +0000 Subject: [PATCH 18/58] WIP coeff_a_alt add additional coefficients --- src/nucfd_coeffs_mod.f90 | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/nucfd_coeffs_mod.f90 b/src/nucfd_coeffs_mod.f90 index 0828b8d..4ff4317 100644 --- a/src/nucfd_coeffs_mod.f90 +++ b/src/nucfd_coeffs_mod.f90 @@ -127,8 +127,10 @@ real function coeff_a_deltas_alt(h) associate(beta => alpha) ! To match Gamet et al. (1999) coeff_a_deltas_alt = h0 * ((hm1 + h0) * (hp1 + hp2) + 2.0 * hp1 * hp2 * beta) & ! = (14/3) h^3 - + hm1 * h0 * hp2 * (beta - alpha) & ! Should cancel for case alpha = beta - + (h0**2) * (hp2 - hp1) * beta ! Should cancel for constant h + + hm1 * h0 * hp2 * (beta - alpha) & ! (G) Should cancel for case alpha = beta + + (h0**2) * (hp2 - hp1) * beta & ! (R) Should cancel for constant h + + hm1 * hp1 * (2.0 * hp2 - h0 - hp1) * beta ! (Y) Should cancel for constant h + print *, "+++", coeff_a_deltas_alt / (h0**3), 14.0 / 3.0, "+++" coeff_a_deltas_alt = coeff_a_deltas_alt & / (3.0 * hp1 * ((h0 + hp1) / 2.0) * hp2) ! => 14.0 / 9.0 when h = const From 36b1d1aa7b2fe47b9bd0c1d343a6216702cfecab Mon Sep 17 00:00:00 2001 From: Paul Date: Mon, 27 Feb 2023 23:26:49 +0000 Subject: [PATCH 19/58] WIP coeff_a_alt add additional terms --- src/nucfd_coeffs_mod.f90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/nucfd_coeffs_mod.f90 b/src/nucfd_coeffs_mod.f90 index 4ff4317..43e6134 100644 --- a/src/nucfd_coeffs_mod.f90 +++ b/src/nucfd_coeffs_mod.f90 @@ -129,7 +129,8 @@ real function coeff_a_deltas_alt(h) coeff_a_deltas_alt = h0 * ((hm1 + h0) * (hp1 + hp2) + 2.0 * hp1 * hp2 * beta) & ! = (14/3) h^3 + hm1 * h0 * hp2 * (beta - alpha) & ! (G) Should cancel for case alpha = beta + (h0**2) * (hp2 - hp1) * beta & ! (R) Should cancel for constant h - + hm1 * hp1 * (2.0 * hp2 - h0 - hp1) * beta ! (Y) Should cancel for constant h + + hm1 * hp1 * (2.0 * hp2 - h0 - hp1) * beta & ! (Y) Should cancel for constant h + + (hp1**3) * (3.0 * hp2 - 2.0 * h0 - hp1) * beta ! (B) Should cancel for constant h print *, "+++", coeff_a_deltas_alt / (h0**3), 14.0 / 3.0, "+++" coeff_a_deltas_alt = coeff_a_deltas_alt & From 0a2ce4c45cc9f295bba8118e36be16852e2aa26e Mon Sep 17 00:00:00 2001 From: Paul Date: Mon, 27 Feb 2023 23:29:58 +0000 Subject: [PATCH 20/58] WIP coeff_a_alt add additional terms --- src/nucfd_coeffs_mod.f90 | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/nucfd_coeffs_mod.f90 b/src/nucfd_coeffs_mod.f90 index 43e6134..82709ea 100644 --- a/src/nucfd_coeffs_mod.f90 +++ b/src/nucfd_coeffs_mod.f90 @@ -129,8 +129,9 @@ real function coeff_a_deltas_alt(h) coeff_a_deltas_alt = h0 * ((hm1 + h0) * (hp1 + hp2) + 2.0 * hp1 * hp2 * beta) & ! = (14/3) h^3 + hm1 * h0 * hp2 * (beta - alpha) & ! (G) Should cancel for case alpha = beta + (h0**2) * (hp2 - hp1) * beta & ! (R) Should cancel for constant h - + hm1 * hp1 * (2.0 * hp2 - h0 - hp1) * beta & ! (Y) Should cancel for constant h - + (hp1**3) * (3.0 * hp2 - 2.0 * h0 - hp1) * beta ! (B) Should cancel for constant h + + hm1 * hp1 * (2.0 * hp2 - h0 - hp1) * beta & ! (Y) Should cancel for constant h + + (hp1**3) * (3.0 * hp2 - 2.0 * h0 - hp1) * beta & ! (B) Should cancel for constant h + + h0 * (2.0 * hp1 * hp2 * beta - hm1 * (h0 + hp1) * alpha) ! (P) Should cancel for constant h print *, "+++", coeff_a_deltas_alt / (h0**3), 14.0 / 3.0, "+++" coeff_a_deltas_alt = coeff_a_deltas_alt & From fe2ffe0e5c195d66ab3c0399ae4905da4404c21c Mon Sep 17 00:00:00 2001 From: Paul Date: Mon, 27 Feb 2023 23:37:34 +0000 Subject: [PATCH 21/58] Remove colour codes --- src/nucfd_coeffs_mod.f90 | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/nucfd_coeffs_mod.f90 b/src/nucfd_coeffs_mod.f90 index 82709ea..32ee4cc 100644 --- a/src/nucfd_coeffs_mod.f90 +++ b/src/nucfd_coeffs_mod.f90 @@ -127,11 +127,11 @@ real function coeff_a_deltas_alt(h) associate(beta => alpha) ! To match Gamet et al. (1999) coeff_a_deltas_alt = h0 * ((hm1 + h0) * (hp1 + hp2) + 2.0 * hp1 * hp2 * beta) & ! = (14/3) h^3 - + hm1 * h0 * hp2 * (beta - alpha) & ! (G) Should cancel for case alpha = beta - + (h0**2) * (hp2 - hp1) * beta & ! (R) Should cancel for constant h - + hm1 * hp1 * (2.0 * hp2 - h0 - hp1) * beta & ! (Y) Should cancel for constant h - + (hp1**3) * (3.0 * hp2 - 2.0 * h0 - hp1) * beta & ! (B) Should cancel for constant h - + h0 * (2.0 * hp1 * hp2 * beta - hm1 * (h0 + hp1) * alpha) ! (P) Should cancel for constant h + + 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**3) * (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 print *, "+++", coeff_a_deltas_alt / (h0**3), 14.0 / 3.0, "+++" coeff_a_deltas_alt = coeff_a_deltas_alt & From 5848e59c35e7b2ab91a043768cc2773464be563b Mon Sep 17 00:00:00 2001 From: Paul Date: Mon, 27 Feb 2023 23:42:12 +0000 Subject: [PATCH 22/58] Add the test for coeff_a_alt --- tests/fd-schemes/verify_coeffs.f90 | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/tests/fd-schemes/verify_coeffs.f90 b/tests/fd-schemes/verify_coeffs.f90 index 9ccb8cc..c7f58e4 100644 --- a/tests/fd-schemes/verify_coeffs.f90 +++ b/tests/fd-schemes/verify_coeffs.f90 @@ -57,7 +57,10 @@ program verify_coeffs call test_report("Coefficient D", check_scalar(d, dref / (4.0 * h))) e = coeff_e(stencil) call test_report("Coefficient E", check_scalar(e, eref)) - + + a = coeff_a_alt(stencil) + call test_report("Coefficient A (alt)", check_scalar(a, aref / (2.0 * h))) + select type(points => stencil%stencil) type is(real) points(-3) = -5.0 * h From 37530bfe5a7e62a888eb1389da25813de38092f4 Mon Sep 17 00:00:00 2001 From: Paul Date: Mon, 27 Feb 2023 23:47:26 +0000 Subject: [PATCH 23/58] Replace coeff_a with coeff_a_alt Rearranges coefficient computation to evaluate terms that should cancel together. --- src/nucfd_coeffs_mod.f90 | 63 +++--------------------------- tests/fd-schemes/verify_coeffs.f90 | 3 -- 2 files changed, 6 insertions(+), 60 deletions(-) diff --git a/src/nucfd_coeffs_mod.f90 b/src/nucfd_coeffs_mod.f90 index 32ee4cc..a80d8cd 100644 --- a/src/nucfd_coeffs_mod.f90 +++ b/src/nucfd_coeffs_mod.f90 @@ -8,7 +8,7 @@ module nucfd_coeffs implicit none private - public :: coeff_a, coeff_a_alt + public :: coeff_a public :: coeff_b public :: coeff_c public :: coeff_d @@ -20,12 +20,6 @@ module nucfd_coeffs module procedure coeff_a_deltas end interface coeff_a - interface coeff_a_alt - !! Compute the coefficient acting on f_{i+1} of the finite diference stencil. - module procedure coeff_a_points_alt - module procedure coeff_a_deltas_alt - end interface coeff_a_alt - interface coeff_b !! Compute the coefficient acting on f_{i-1} of the finite diference stencil. module procedure coeff_b_points @@ -82,65 +76,20 @@ pure real function coeff_a_deltas(h) end select associate(beta => alpha) ! To match Gamet et al. (1999) - coeff_a_deltas = hm1 * h0 * hp1 + h0**2 * hp1 + hm1 * h0 * hp2 + h0**2 * hp2 & - - hm1 * h0**2 * alpha - hm1 * h0 * hp1 * alpha & - - hm1 * h0 * hp2 * alpha - hm1 * h0 * hp1 * beta - h0**2 * hp1 * beta & - - hm1 * hp1**2 * beta - 2.0 * h0 * hp1**2 * beta - hp1**3 * beta & - + hm1 * h0 * hp2 * beta + h0**2 * hp2 * beta + 2.0 * hm1 * hp1 * hp2 * beta & - + 4.0 * h0 * hp1 * hp2 * beta + 3.0 * hp1**2 * hp2 * beta - coeff_a_deltas = coeff_a_deltas & - / (hp1 * (h0 + hp1) * (hm1 + h0 + hp1) * hp2) - end associate - end function coeff_a_deltas - - real function coeff_a_points_alt(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. - - h = points_to_deltas(x) - coeff_a_points_alt = coeff_a_alt(h) - - end function coeff_a_points_alt - - real function coeff_a_deltas_alt(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 :: hm2, hm1, h0, hp1, hp2 ! Grid deltas at i -2, -1, 0, +1, +2 - - select type(deltas => h%stencil) - type is(real) - hm2 = deltas(-2) - 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) - coeff_a_deltas_alt = h0 * ((hm1 + h0) * (hp1 + hp2) + 2.0 * hp1 * hp2 * beta) & ! = (14/3) h^3 + coeff_a_deltas = h0 * ((hm1 + h0) * (hp1 + hp2) + 2.0 * hp1 * hp2 * beta) & ! = (14/3) h^3 + 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**3) * (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 - print *, "+++", coeff_a_deltas_alt / (h0**3), 14.0 / 3.0, "+++" - coeff_a_deltas_alt = coeff_a_deltas_alt & + coeff_a_deltas = coeff_a_deltas & / (3.0 * hp1 * ((h0 + hp1) / 2.0) * hp2) ! => 14.0 / 9.0 when h = const - print *, "+++", coeff_a_deltas_alt, 14.0 / 9.0, "+++" - coeff_a_deltas_alt = coeff_a_deltas_alt & + + coeff_a_deltas = coeff_a_deltas & / (2.0 * ((hm1 + h0 + hp1) / 3.0)) ! => (14.0 / 9.0) / (2h) when h = const end associate - end function coeff_a_deltas_alt + end function coeff_a_deltas real function coeff_b_points(x) !! Compute the coefficient acting on f_{i-1} of the finite diference given a stencil of points. diff --git a/tests/fd-schemes/verify_coeffs.f90 b/tests/fd-schemes/verify_coeffs.f90 index c7f58e4..ef63300 100644 --- a/tests/fd-schemes/verify_coeffs.f90 +++ b/tests/fd-schemes/verify_coeffs.f90 @@ -57,9 +57,6 @@ program verify_coeffs call test_report("Coefficient D", check_scalar(d, dref / (4.0 * h))) e = coeff_e(stencil) call test_report("Coefficient E", check_scalar(e, eref)) - - a = coeff_a_alt(stencil) - call test_report("Coefficient A (alt)", check_scalar(a, aref / (2.0 * h))) select type(points => stencil%stencil) type is(real) From be787e712aba96acb0abeaf085b419315d72570c Mon Sep 17 00:00:00 2001 From: Paul Bartholomew Date: Tue, 28 Feb 2023 13:21:51 +0000 Subject: [PATCH 24/58] Initial alternative implementation of coeff_b --- src/nucfd_coeffs_mod.f90 | 52 +++++++++++++++++++++++++++++- tests/fd-schemes/verify_coeffs.f90 | 3 ++ 2 files changed, 54 insertions(+), 1 deletion(-) diff --git a/src/nucfd_coeffs_mod.f90 b/src/nucfd_coeffs_mod.f90 index a80d8cd..d721ca4 100644 --- a/src/nucfd_coeffs_mod.f90 +++ b/src/nucfd_coeffs_mod.f90 @@ -9,7 +9,7 @@ module nucfd_coeffs private public :: coeff_a - public :: coeff_b + public :: coeff_b, coeff_b_alt public :: coeff_c public :: coeff_d public :: coeff_e @@ -25,6 +25,11 @@ module nucfd_coeffs module procedure coeff_b_points module procedure coeff_b_deltas end interface coeff_b + interface coeff_b_alt + !! Compute the coefficient acting on f_{i-1} of the finite diference stencil. + module procedure coeff_b_alt_points + module procedure coeff_b_alt_deltas + end interface coeff_b_alt interface coeff_c !! Compute the coefficient acting on f_{i+2} of the finite diference stencil. @@ -135,6 +140,51 @@ pure real function coeff_b_deltas(h) end associate end function coeff_b_deltas + real function coeff_b_alt_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. + + h = points_to_deltas(x) + coeff_b_alt_points = coeff_b_alt_deltas(h) + + end function coeff_b_alt_points + + real function coeff_b_alt_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 :: hm2, hm1, h0, hp1, hp2 ! Grid deltas at i -2, -1, 0, +1, +2 + + select type(deltas => h%stencil) + type is(real) + hm2 = deltas(-2) + 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) + coeff_b_alt_deltas = -hp1 * ((hm1 + h0) * (hp1 + hp2) + 2.0 * hm1 * h0 * alpha) ! = -(14/3) h^3 + print *, "+++", coeff_b_alt_deltas, -(14.0 / 3.0) * (h0**3), "+++" + + coeff_b_alt_deltas = coeff_b_alt_deltas & + / (3.0 * hm1 * h0 * ((h0 + hp1) / 2.0)) ! => -14.0 / 9.0 when h = const + print *, "+++", coeff_b_alt_deltas, -(14.0 / 9.0), "+++" + + coeff_b_alt_deltas = coeff_b_alt_deltas & + / (2.0 * (h0 + hp1 + hp2) / 3.0) ! => -(14.0 / 9.0) / (2h) when h = const + end associate + end function coeff_b_alt_deltas + real function coeff_c_points(x) !! Compute the coefficient acting on f_{i+2} of the finite diference given a stencil of points. diff --git a/tests/fd-schemes/verify_coeffs.f90 b/tests/fd-schemes/verify_coeffs.f90 index ef63300..454f211 100644 --- a/tests/fd-schemes/verify_coeffs.f90 +++ b/tests/fd-schemes/verify_coeffs.f90 @@ -57,6 +57,9 @@ program verify_coeffs call test_report("Coefficient D", check_scalar(d, dref / (4.0 * h))) e = coeff_e(stencil) call test_report("Coefficient E", check_scalar(e, eref)) + + b = coeff_b_alt(stencil) + call test_report("Coefficient B (ALT)", check_scalar(b, bref / (2.0 * h))) select type(points => stencil%stencil) type is(real) From 34c39c6495c56bebe5d022a806308983c3172025 Mon Sep 17 00:00:00 2001 From: Paul Bartholomew Date: Tue, 28 Feb 2023 13:25:37 +0000 Subject: [PATCH 25/58] WIP add additional terms to coeff_b_alt These should cancel for constant h --- src/nucfd_coeffs_mod.f90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/nucfd_coeffs_mod.f90 b/src/nucfd_coeffs_mod.f90 index d721ca4..7098789 100644 --- a/src/nucfd_coeffs_mod.f90 +++ b/src/nucfd_coeffs_mod.f90 @@ -173,7 +173,8 @@ real function coeff_b_alt_deltas(h) end select associate(beta => alpha) ! To match Gamet et al. (1999) - coeff_b_alt_deltas = -hp1 * ((hm1 + h0) * (hp1 + hp2) + 2.0 * hm1 * h0 * alpha) ! = -(14/3) h^3 + coeff_b_alt_deltas = -hp1 * ((hm1 + h0) * (hp1 + hp2) + 2.0 * hm1 * h0 * alpha) & ! = -(14/3) h^3 + - hm1 * hp1 * hp2 * (alpha - beta) ! Should cancel for case alpha = beta print *, "+++", coeff_b_alt_deltas, -(14.0 / 3.0) * (h0**3), "+++" coeff_b_alt_deltas = coeff_b_alt_deltas & From 4d5f63fad70f888cf3fd176d7149d1bbfe41a59e Mon Sep 17 00:00:00 2001 From: Paul Bartholomew Date: Tue, 28 Feb 2023 13:28:13 +0000 Subject: [PATCH 26/58] Adding additional terms to coeff_b_alt These cancel for constant h --- src/nucfd_coeffs_mod.f90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/nucfd_coeffs_mod.f90 b/src/nucfd_coeffs_mod.f90 index 7098789..90c5e83 100644 --- a/src/nucfd_coeffs_mod.f90 +++ b/src/nucfd_coeffs_mod.f90 @@ -174,7 +174,8 @@ real function coeff_b_alt_deltas(h) associate(beta => alpha) ! To match Gamet et al. (1999) coeff_b_alt_deltas = -hp1 * ((hm1 + h0) * (hp1 + hp2) + 2.0 * hm1 * h0 * alpha) & ! = -(14/3) h^3 - - hm1 * hp1 * hp2 * (alpha - beta) ! Should cancel for case alpha = beta + - hm1 * hp1 * hp2 * (alpha - beta) & ! (G) Should cancel for case alpha = beta + - (hp1**2) * (hm1 - h0) * alpha ! (R) Should cancel for case h = const print *, "+++", coeff_b_alt_deltas, -(14.0 / 3.0) * (h0**3), "+++" coeff_b_alt_deltas = coeff_b_alt_deltas & From 052935b1d817ef0a76dcf3523753cba28b36d265 Mon Sep 17 00:00:00 2001 From: Paul Bartholomew Date: Tue, 28 Feb 2023 13:30:41 +0000 Subject: [PATCH 27/58] Adding additional terms to coeff_b_alt These should cancel for h = const --- src/nucfd_coeffs_mod.f90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/nucfd_coeffs_mod.f90 b/src/nucfd_coeffs_mod.f90 index 90c5e83..b59bca2 100644 --- a/src/nucfd_coeffs_mod.f90 +++ b/src/nucfd_coeffs_mod.f90 @@ -175,7 +175,8 @@ real function coeff_b_alt_deltas(h) associate(beta => alpha) ! To match Gamet et al. (1999) coeff_b_alt_deltas = -hp1 * ((hm1 + h0) * (hp1 + hp2) + 2.0 * hm1 * h0 * alpha) & ! = -(14/3) h^3 - hm1 * hp1 * hp2 * (alpha - beta) & ! (G) Should cancel for case alpha = beta - - (hp1**2) * (hm1 - h0) * alpha ! (R) Should cancel for case h = const + - (hp1**2) * (hm1 - h0) * alpha & ! (R) Should cancel for case h = const + - h0 * hp2 * (2.0 * hm1 - h0 - hp1) * alpha !(Y) Should cancel for case h = const print *, "+++", coeff_b_alt_deltas, -(14.0 / 3.0) * (h0**3), "+++" coeff_b_alt_deltas = coeff_b_alt_deltas & From 6df5f75640326e6fe5afeafce41424836f11f761 Mon Sep 17 00:00:00 2001 From: Paul Bartholomew Date: Tue, 28 Feb 2023 13:32:36 +0000 Subject: [PATCH 28/58] Adding additional terms to coeff_b_alt These cancel under h= const --- src/nucfd_coeffs_mod.f90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/nucfd_coeffs_mod.f90 b/src/nucfd_coeffs_mod.f90 index b59bca2..d358f75 100644 --- a/src/nucfd_coeffs_mod.f90 +++ b/src/nucfd_coeffs_mod.f90 @@ -176,7 +176,8 @@ real function coeff_b_alt_deltas(h) coeff_b_alt_deltas = -hp1 * ((hm1 + h0) * (hp1 + hp2) + 2.0 * hm1 * h0 * alpha) & ! = -(14/3) h^3 - hm1 * hp1 * hp2 * (alpha - beta) & ! (G) Should cancel for case alpha = beta - (hp1**2) * (hm1 - h0) * alpha & ! (R) Should cancel for case h = const - - h0 * hp2 * (2.0 * hm1 - h0 - hp1) * alpha !(Y) Should cancel for case h = const + - h0 * hp2 * (2.0 * hm1 - h0 - hp1) * alpha & !(Y) Should cancel for case h = const + - (h0**2) * (3.0 * hm1 - h0 - 2.0 * hp1) * alpha print *, "+++", coeff_b_alt_deltas, -(14.0 / 3.0) * (h0**3), "+++" coeff_b_alt_deltas = coeff_b_alt_deltas & From 75eae924d5b5f4fd7f92c76da20a47ce5e93fc3e Mon Sep 17 00:00:00 2001 From: Paul Bartholomew Date: Tue, 28 Feb 2023 13:33:28 +0000 Subject: [PATCH 29/58] Spotted error in coeff_A calculation --- src/nucfd_coeffs_mod.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/nucfd_coeffs_mod.f90 b/src/nucfd_coeffs_mod.f90 index d358f75..665c482 100644 --- a/src/nucfd_coeffs_mod.f90 +++ b/src/nucfd_coeffs_mod.f90 @@ -85,7 +85,7 @@ pure real function coeff_a_deltas(h) + 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**3) * (3.0 * hp2 - 2.0 * 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 coeff_a_deltas = coeff_a_deltas & @@ -177,7 +177,7 @@ real function coeff_b_alt_deltas(h) - hm1 * hp1 * hp2 * (alpha - beta) & ! (G) Should cancel for case alpha = beta - (hp1**2) * (hm1 - h0) * alpha & ! (R) Should cancel for case h = const - h0 * hp2 * (2.0 * hm1 - h0 - hp1) * alpha & !(Y) Should cancel for case h = const - - (h0**2) * (3.0 * hm1 - h0 - 2.0 * hp1) * alpha + - (h0**2) * (3.0 * hm1 - h0 - 2.0 * hp1) * alpha ! (B) Should cancel for case h = const print *, "+++", coeff_b_alt_deltas, -(14.0 / 3.0) * (h0**3), "+++" coeff_b_alt_deltas = coeff_b_alt_deltas & From eeb047753d5a0f776ca49d07c3a448b2e5c07609 Mon Sep 17 00:00:00 2001 From: Paul Bartholomew Date: Tue, 28 Feb 2023 13:36:19 +0000 Subject: [PATCH 30/58] Implement last term in coeff_b_alt --- src/nucfd_coeffs_mod.f90 | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/nucfd_coeffs_mod.f90 b/src/nucfd_coeffs_mod.f90 index 665c482..a045296 100644 --- a/src/nucfd_coeffs_mod.f90 +++ b/src/nucfd_coeffs_mod.f90 @@ -176,8 +176,9 @@ real function coeff_b_alt_deltas(h) coeff_b_alt_deltas = -hp1 * ((hm1 + h0) * (hp1 + hp2) + 2.0 * hm1 * h0 * alpha) & ! = -(14/3) h^3 - hm1 * hp1 * hp2 * (alpha - beta) & ! (G) Should cancel for case alpha = beta - (hp1**2) * (hm1 - h0) * alpha & ! (R) Should cancel for case h = const - - h0 * hp2 * (2.0 * hm1 - h0 - hp1) * alpha & !(Y) Should cancel for case h = const - - (h0**2) * (3.0 * hm1 - h0 - 2.0 * hp1) * alpha ! (B) Should cancel for case h = const + - h0 * hp2 * (2.0 * hm1 - h0 - hp1) * alpha & ! (Y) Should cancel for case h = const + - (h0**2) * (3.0 * hm1 - h0 - 2.0 * hp1) * alpha & ! (B) Should cancel for case h = const + - hp1 * (2.0 * hm1 * h0 * alpha - (h0 + hp1) * hp2 * beta) ! (P) Should cancel for constant h & alpha = beta print *, "+++", coeff_b_alt_deltas, -(14.0 / 3.0) * (h0**3), "+++" coeff_b_alt_deltas = coeff_b_alt_deltas & From bfc178b63bbc0b0fcf44a5a17c29153f047ab50d Mon Sep 17 00:00:00 2001 From: Paul Bartholomew Date: Tue, 28 Feb 2023 13:36:49 +0000 Subject: [PATCH 31/58] Remove colour codes --- src/nucfd_coeffs_mod.f90 | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/nucfd_coeffs_mod.f90 b/src/nucfd_coeffs_mod.f90 index a045296..c48811b 100644 --- a/src/nucfd_coeffs_mod.f90 +++ b/src/nucfd_coeffs_mod.f90 @@ -174,11 +174,11 @@ real function coeff_b_alt_deltas(h) associate(beta => alpha) ! To match Gamet et al. (1999) coeff_b_alt_deltas = -hp1 * ((hm1 + h0) * (hp1 + hp2) + 2.0 * hm1 * h0 * alpha) & ! = -(14/3) h^3 - - hm1 * hp1 * hp2 * (alpha - beta) & ! (G) Should cancel for case alpha = beta - - (hp1**2) * (hm1 - h0) * alpha & ! (R) Should cancel for case h = const - - h0 * hp2 * (2.0 * hm1 - h0 - hp1) * alpha & ! (Y) Should cancel for case h = const - - (h0**2) * (3.0 * hm1 - h0 - 2.0 * hp1) * alpha & ! (B) Should cancel for case h = const - - hp1 * (2.0 * hm1 * h0 * alpha - (h0 + hp1) * hp2 * beta) ! (P) Should cancel for constant h & alpha = beta + - hm1 * hp1 * hp2 * (alpha - beta) & ! Should cancel for case alpha = beta + - (hp1**2) * (hm1 - h0) * alpha & ! Should cancel for case h = const + - h0 * hp2 * (2.0 * hm1 - h0 - hp1) * alpha & ! Should cancel for case h = const + - (h0**2) * (3.0 * hm1 - h0 - 2.0 * hp1) * alpha & ! Should cancel for case h = const + - hp1 * (2.0 * hm1 * h0 * alpha - (h0 + hp1) * hp2 * beta) ! Should cancel for constant h & alpha = beta print *, "+++", coeff_b_alt_deltas, -(14.0 / 3.0) * (h0**3), "+++" coeff_b_alt_deltas = coeff_b_alt_deltas & From 7fea668a8e9fa18929ae74d6286f454a455177a8 Mon Sep 17 00:00:00 2001 From: Paul Bartholomew Date: Tue, 28 Feb 2023 13:39:34 +0000 Subject: [PATCH 32/58] Replace coeff_b with coeff_b_alt implementation This should be more resilient to round-off error. --- src/nucfd_coeffs_mod.f90 | 65 ++++-------------------------- tests/fd-schemes/verify_coeffs.f90 | 3 -- 2 files changed, 7 insertions(+), 61 deletions(-) diff --git a/src/nucfd_coeffs_mod.f90 b/src/nucfd_coeffs_mod.f90 index c48811b..e196f4e 100644 --- a/src/nucfd_coeffs_mod.f90 +++ b/src/nucfd_coeffs_mod.f90 @@ -9,7 +9,7 @@ module nucfd_coeffs private public :: coeff_a - public :: coeff_b, coeff_b_alt + public :: coeff_b public :: coeff_c public :: coeff_d public :: coeff_e @@ -19,17 +19,12 @@ module nucfd_coeffs module procedure coeff_a_points module procedure coeff_a_deltas end interface coeff_a - + interface coeff_b !! Compute the coefficient acting on f_{i-1} of the finite diference stencil. module procedure coeff_b_points module procedure coeff_b_deltas end interface coeff_b - interface coeff_b_alt - !! Compute the coefficient acting on f_{i-1} of the finite diference stencil. - module procedure coeff_b_alt_points - module procedure coeff_b_alt_deltas - end interface coeff_b_alt interface coeff_c !! Compute the coefficient acting on f_{i+2} of the finite diference stencil. @@ -109,7 +104,7 @@ real function coeff_b_points(x) end function coeff_b_points - pure real function coeff_b_deltas(h) + real function coeff_b_deltas(h) !! Compute the coefficient acting on f_{i-1} of the finite diference given a stencil of grid !! spacings. @@ -129,66 +124,20 @@ pure real function coeff_b_deltas(h) end select associate(beta => alpha) ! To match Gamet et al. (1999) - coeff_b_deltas = -hm1 * hp1**2 - h0 * hp1**2 - hm1 * hp1 * hp2 - h0 * hp1 * hp2 & - - 3.0 * hm1 * h0**2 * alpha - 4.0 * hm1 * h0 * hp1 * alpha & ! End line 1 - + h0**3 * alpha + 2.0 * h0**2 * hp1 * alpha - hm1 * hp1**2 * alpha & - + h0 * hp1**2 * alpha - 2.0 * hm1 * h0 * hp2 * alpha - hm1 * hp1 * hp2 * alpha & ! End line 2 - + h0**2 * hp2 * alpha + h0 * hp1 * hp2 * alpha + hm1 * hp1 * hp2 * beta & - + h0 * hp1 * hp2 * beta + hp1**2 * hp2 * beta ! End line 3 - coeff_b_deltas = coeff_b_deltas & - / (hm1 * h0 * (h0 + hp1) * (h0 + hp1 + hp2)) - end associate - end function coeff_b_deltas - - real function coeff_b_alt_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. - - h = points_to_deltas(x) - coeff_b_alt_points = coeff_b_alt_deltas(h) - - end function coeff_b_alt_points - - real function coeff_b_alt_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 :: hm2, hm1, h0, hp1, hp2 ! Grid deltas at i -2, -1, 0, +1, +2 - - select type(deltas => h%stencil) - type is(real) - hm2 = deltas(-2) - 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) - coeff_b_alt_deltas = -hp1 * ((hm1 + h0) * (hp1 + hp2) + 2.0 * hm1 * h0 * alpha) & ! = -(14/3) h^3 + coeff_b_deltas = -hp1 * ((hm1 + h0) * (hp1 + hp2) + 2.0 * hm1 * h0 * alpha) & ! = -(14/3) h^3 - hm1 * hp1 * hp2 * (alpha - beta) & ! Should cancel for case alpha = beta - (hp1**2) * (hm1 - h0) * alpha & ! Should cancel for case h = const - h0 * hp2 * (2.0 * hm1 - h0 - hp1) * alpha & ! Should cancel for case h = const - (h0**2) * (3.0 * hm1 - h0 - 2.0 * hp1) * alpha & ! Should cancel for case h = const - hp1 * (2.0 * hm1 * h0 * alpha - (h0 + hp1) * hp2 * beta) ! Should cancel for constant h & alpha = beta - print *, "+++", coeff_b_alt_deltas, -(14.0 / 3.0) * (h0**3), "+++" - coeff_b_alt_deltas = coeff_b_alt_deltas & + coeff_b_deltas = coeff_b_deltas & / (3.0 * hm1 * h0 * ((h0 + hp1) / 2.0)) ! => -14.0 / 9.0 when h = const - print *, "+++", coeff_b_alt_deltas, -(14.0 / 9.0), "+++" - coeff_b_alt_deltas = coeff_b_alt_deltas & + coeff_b_deltas = coeff_b_deltas & / (2.0 * (h0 + hp1 + hp2) / 3.0) ! => -(14.0 / 9.0) / (2h) when h = const end associate - end function coeff_b_alt_deltas + end function coeff_b_deltas real function coeff_c_points(x) !! Compute the coefficient acting on f_{i+2} of the finite diference given a stencil of points. diff --git a/tests/fd-schemes/verify_coeffs.f90 b/tests/fd-schemes/verify_coeffs.f90 index 454f211..ef63300 100644 --- a/tests/fd-schemes/verify_coeffs.f90 +++ b/tests/fd-schemes/verify_coeffs.f90 @@ -57,9 +57,6 @@ program verify_coeffs call test_report("Coefficient D", check_scalar(d, dref / (4.0 * h))) e = coeff_e(stencil) call test_report("Coefficient E", check_scalar(e, eref)) - - b = coeff_b_alt(stencil) - call test_report("Coefficient B (ALT)", check_scalar(b, bref / (2.0 * h))) select type(points => stencil%stencil) type is(real) From 0d70de236e3003926cb5983e4ad63ec8a0896376 Mon Sep 17 00:00:00 2001 From: Paul Bartholomew Date: Tue, 28 Feb 2023 15:23:46 +0000 Subject: [PATCH 33/58] Move coeffs module to subdirectory In preparation for submodules --- CMakeLists.txt | 2 +- src/{ => coeffs}/nucfd_coeffs_mod.f90 | 0 2 files changed, 1 insertion(+), 1 deletion(-) rename src/{ => coeffs}/nucfd_coeffs_mod.f90 (100%) diff --git a/CMakeLists.txt b/CMakeLists.txt index ed9c2ee..498c764 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -40,7 +40,7 @@ target_compile_options(nucfd_strict_flags INTERFACE add_library(nucfd src/nucfd_trid_solver_mod.f90 src/nucfd_types_mod.f90 - src/nucfd_coeffs_mod.f90 + src/coeffs/nucfd_coeffs_mod.f90 src/nucfd_deriv_mod.f90) target_link_libraries(nucfd nucfd_compiler_flags) if (CMAKE_BUILD_TYPE MATCHES "Debug") diff --git a/src/nucfd_coeffs_mod.f90 b/src/coeffs/nucfd_coeffs_mod.f90 similarity index 100% rename from src/nucfd_coeffs_mod.f90 rename to src/coeffs/nucfd_coeffs_mod.f90 From 155cd1be0e410f591af8b57914654707fccd1e34 Mon Sep 17 00:00:00 2001 From: Paul Bartholomew Date: Tue, 28 Feb 2023 15:45:16 +0000 Subject: [PATCH 34/58] Move coefficient evaluations to submodules --- CMakeLists.txt | 2 +- src/coeffs/nucfd_coeffs_a.f90 | 55 ++++++++ src/coeffs/nucfd_coeffs_b.f90 | 54 ++++++++ src/coeffs/nucfd_coeffs_c.f90 | 48 +++++++ src/coeffs/nucfd_coeffs_d.f90 | 47 +++++++ src/coeffs/nucfd_coeffs_e.f90 | 24 ++++ src/coeffs/nucfd_coeffs_mod.f90 | 235 ++++++-------------------------- 7 files changed, 268 insertions(+), 197 deletions(-) create mode 100644 src/coeffs/nucfd_coeffs_a.f90 create mode 100644 src/coeffs/nucfd_coeffs_b.f90 create mode 100644 src/coeffs/nucfd_coeffs_c.f90 create mode 100644 src/coeffs/nucfd_coeffs_d.f90 create mode 100644 src/coeffs/nucfd_coeffs_e.f90 diff --git a/CMakeLists.txt b/CMakeLists.txt index 498c764..762f528 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -40,7 +40,7 @@ target_compile_options(nucfd_strict_flags INTERFACE add_library(nucfd src/nucfd_trid_solver_mod.f90 src/nucfd_types_mod.f90 - src/coeffs/nucfd_coeffs_mod.f90 + src/coeffs/nucfd_coeffs_mod.f90 src/coeffs/nucfd_coeffs_a.f90 src/coeffs/nucfd_coeffs_b.f90 src/coeffs/nucfd_coeffs_c.f90 src/coeffs/nucfd_coeffs_d.f90 src/coeffs/nucfd_coeffs_e.f90 src/nucfd_deriv_mod.f90) target_link_libraries(nucfd nucfd_compiler_flags) if (CMAKE_BUILD_TYPE MATCHES "Debug") diff --git a/src/coeffs/nucfd_coeffs_a.f90 b/src/coeffs/nucfd_coeffs_a.f90 new file mode 100644 index 0000000..add0723 --- /dev/null +++ b/src/coeffs/nucfd_coeffs_a.f90 @@ -0,0 +1,55 @@ +submodule (nucfd_coeffs) nucfd_coeffs_a + + 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. + + h = points_to_deltas(x) + coeff_a_points = coeff_a(h) + + end function coeff_a_points + + module pure 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 :: hm2, hm1, h0, hp1, hp2 ! Grid deltas at i -2, -1, 0, +1, +2 + + select type(deltas => h%stencil) + type is(real) + hm2 = deltas(-2) + 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) + coeff_a_deltas = h0 * ((hm1 + h0) * (hp1 + hp2) + 2.0 * hp1 * hp2 * beta) & ! = (14/3) h^3 + + 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 + + coeff_a_deltas = coeff_a_deltas & + / (3.0 * hp1 * ((h0 + hp1) / 2.0) * hp2) ! => 14.0 / 9.0 when h = const + + coeff_a_deltas = coeff_a_deltas & + / (2.0 * ((hm1 + h0 + hp1) / 3.0)) ! => (14.0 / 9.0) / (2h) when h = const + end associate + end function coeff_a_deltas + +end submodule nucfd_coeffs_a diff --git a/src/coeffs/nucfd_coeffs_b.f90 b/src/coeffs/nucfd_coeffs_b.f90 new file mode 100644 index 0000000..fab288b --- /dev/null +++ b/src/coeffs/nucfd_coeffs_b.f90 @@ -0,0 +1,54 @@ +submodule (nucfd_coeffs) nucfd_coeffs_b + + implicit none + +contains + + module real function coeff_b_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. + + h = points_to_deltas(x) + coeff_b_points = coeff_b_deltas(h) + + end function coeff_b_points + + module pure real function coeff_b_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 :: hm2, hm1, h0, hp1, hp2 ! Grid deltas at i -2, -1, 0, +1, +2 + + select type(deltas => h%stencil) + type is(real) + hm2 = deltas(-2) + 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) + coeff_b_deltas = -hp1 * ((hm1 + h0) * (hp1 + hp2) + 2.0 * hm1 * h0 * alpha) & ! = -(14/3) h^3 + - hm1 * hp1 * hp2 * (alpha - beta) & ! Should cancel for case alpha = beta + - (hp1**2) * (hm1 - h0) * alpha & ! Should cancel for case h = const + - h0 * hp2 * (2.0 * hm1 - h0 - hp1) * alpha & ! Should cancel for case h = const + - (h0**2) * (3.0 * hm1 - h0 - 2.0 * hp1) * alpha & ! Should cancel for case h = const + - hp1 * (2.0 * hm1 * h0 * alpha - (h0 + hp1) * hp2 * beta) ! Should cancel for constant h & alpha = beta + + coeff_b_deltas = coeff_b_deltas & + / (3.0 * hm1 * h0 * ((h0 + hp1) / 2.0)) ! => -14.0 / 9.0 when h = const + + coeff_b_deltas = coeff_b_deltas & + / (2.0 * (h0 + hp1 + hp2) / 3.0) ! => -(14.0 / 9.0) / (2h) when h = const + end associate + end function coeff_b_deltas +end submodule nucfd_coeffs_b diff --git a/src/coeffs/nucfd_coeffs_c.f90 b/src/coeffs/nucfd_coeffs_c.f90 new file mode 100644 index 0000000..a1616a9 --- /dev/null +++ b/src/coeffs/nucfd_coeffs_c.f90 @@ -0,0 +1,48 @@ +submodule (nucfd_coeffs) nucfd_coeffs_c + + implicit none + +contains + + module real function coeff_c_points(x) + !! Compute the coefficient acting on f_{i+2} 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. + + h = points_to_deltas(x) + coeff_c_points = coeff_c_deltas(h) + + end function coeff_c_points + + module pure real function coeff_c_deltas(h) + !! Compute the coefficient acting on f_{i+2} 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 :: hm2, hm1, h0, hp1, hp2 ! Grid deltas at i -2, -1, 0, +1, +2 + + select type(deltas => h%stencil) + type is(real) + hm2 = deltas(-2) + 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) + coeff_c_deltas = -hm1 * h0 * hp1 - h0**2 * hp1 + hm1 * h0**2 * alpha & + + hm1 * h0 * hp1 * alpha + hm1 * h0 * hp1 * beta + h0**2 * hp1 * beta & ! End line 1 + + hm1 * hp1**2 * beta + 2.0 * h0 * hp1**2 * beta + hp1**3 * beta ! End line 2 + coeff_c_deltas = coeff_c_deltas & + / (hp2 * (hp1 + hp2) * (h0 + hp1 + hp2) * (hm1 + h0 + hp1 + hp2)) + end associate + end function coeff_c_deltas + +end submodule nucfd_coeffs_c diff --git a/src/coeffs/nucfd_coeffs_d.f90 b/src/coeffs/nucfd_coeffs_d.f90 new file mode 100644 index 0000000..72c4cfe --- /dev/null +++ b/src/coeffs/nucfd_coeffs_d.f90 @@ -0,0 +1,47 @@ +submodule (nucfd_coeffs) nucfd_coeffs_d + + implicit none + +contains + + module real function coeff_d_points(x) + !! Compute the coefficient acting on f_{i-2} 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. + + h = points_to_deltas(x) + coeff_d_points = coeff_d_deltas(h) + + end function coeff_d_points + + module pure real function coeff_d_deltas(h) + !! Compute the coefficient acting on f_{i-2} 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 :: hm2, hm1, h0, hp1, hp2 ! Grid deltas at i -2, -1, 0, +1, +2 + + select type(deltas => h%stencil) + type is(real) + hm2 = deltas(-2) + 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) + coeff_d_deltas = h0 * hp1**2 + h0 * hp1 * hp2 - h0**3 * alpha - 2.0 * h0**2 * hp1 * alpha & + - h0 * hp1**2 * alpha - h0**2 * hp2 * alpha - h0 * hp1 * hp2 * alpha & ! End line 1 + - h0 * hp1 * hp2 * beta - hp1**2 * hp2 * beta ! End line 2 + coeff_d_deltas = coeff_d_deltas & + / (hm1 * (hm1 + h0) * (hm1 + h0 + hp1) * (hm1 + h0 + hp1 + hp2)) + end associate + end function coeff_d_deltas +end submodule nucfd_coeffs_d diff --git a/src/coeffs/nucfd_coeffs_e.f90 b/src/coeffs/nucfd_coeffs_e.f90 new file mode 100644 index 0000000..450beb5 --- /dev/null +++ b/src/coeffs/nucfd_coeffs_e.f90 @@ -0,0 +1,24 @@ +submodule (nucfd_coeffs) nucfd_coeffs_e + + implicit none + +contains + + module real function coeff_e(x) + !! Compute the coefficient acting on f_{i} of the finite diference given a stencil of points. + + type(nucfd_stencil_points), intent(in) :: x !! Stencil of points for the finite + !! difference. + + real :: a, b, c, d ! The neighbouring coefficients + + a = coeff_a(x) + b = coeff_b(x) + c = coeff_c(x) + d = coeff_d(x) + + coeff_e = -(a + b + c + d) + + end function coeff_e + +end submodule nucfd_coeffs_e diff --git a/src/coeffs/nucfd_coeffs_mod.f90 b/src/coeffs/nucfd_coeffs_mod.f90 index e196f4e..2a76792 100644 --- a/src/coeffs/nucfd_coeffs_mod.f90 +++ b/src/coeffs/nucfd_coeffs_mod.f90 @@ -37,204 +37,47 @@ module nucfd_coeffs module procedure coeff_d_points module procedure coeff_d_deltas end interface coeff_d + + interface + module real function coeff_a_points(x) + type(nucfd_stencil_points), intent(in) :: x !! Stencil of points for the finite + !! difference. + end function + module pure real function coeff_a_deltas(h) + type(nucfd_stencil_deltas), intent(in) :: h !! Stencil of grid spacings for the finite + !! difference. + end function + module real function coeff_b_points(x) + type(nucfd_stencil_points), intent(in) :: x !! Stencil of points for the finite + !! difference. + end function + module pure real function coeff_b_deltas(h) + type(nucfd_stencil_deltas), intent(in) :: h !! Stencil of grid spacings for the finite + !! difference. + end function + module real function coeff_c_points(x) + type(nucfd_stencil_points), intent(in) :: x !! Stencil of points for the finite + !! difference. + end function + module pure real function coeff_c_deltas(h) + type(nucfd_stencil_deltas), intent(in) :: h !! Stencil of grid spacings for the finite + !! difference. + end function + module real function coeff_d_points(x) + type(nucfd_stencil_points), intent(in) :: x !! Stencil of points for the finite + !! difference. + end function + module pure real function coeff_d_deltas(h) + type(nucfd_stencil_deltas), intent(in) :: h !! Stencil of grid spacings for the finite + !! difference. + end function + module real function coeff_e(x) + type(nucfd_stencil_points), intent(in) :: x !! Stencil of points for the finite + !! difference. + end function + end interface real, parameter, public :: alpha = 1.0 / 3.0 !! Off-diagonal coefficient for first derivative !! system. - -contains - - 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. - - h = points_to_deltas(x) - coeff_a_points = coeff_a(h) - - end function coeff_a_points - - pure 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 :: hm2, hm1, h0, hp1, hp2 ! Grid deltas at i -2, -1, 0, +1, +2 - - select type(deltas => h%stencil) - type is(real) - hm2 = deltas(-2) - 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) - coeff_a_deltas = h0 * ((hm1 + h0) * (hp1 + hp2) + 2.0 * hp1 * hp2 * beta) & ! = (14/3) h^3 - + 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 - - coeff_a_deltas = coeff_a_deltas & - / (3.0 * hp1 * ((h0 + hp1) / 2.0) * hp2) ! => 14.0 / 9.0 when h = const - - coeff_a_deltas = coeff_a_deltas & - / (2.0 * ((hm1 + h0 + hp1) / 3.0)) ! => (14.0 / 9.0) / (2h) when h = const - end associate - end function coeff_a_deltas - - real function coeff_b_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. - - h = points_to_deltas(x) - coeff_b_points = coeff_b_deltas(h) - - end function coeff_b_points - - real function coeff_b_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 :: hm2, hm1, h0, hp1, hp2 ! Grid deltas at i -2, -1, 0, +1, +2 - - select type(deltas => h%stencil) - type is(real) - hm2 = deltas(-2) - 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) - coeff_b_deltas = -hp1 * ((hm1 + h0) * (hp1 + hp2) + 2.0 * hm1 * h0 * alpha) & ! = -(14/3) h^3 - - hm1 * hp1 * hp2 * (alpha - beta) & ! Should cancel for case alpha = beta - - (hp1**2) * (hm1 - h0) * alpha & ! Should cancel for case h = const - - h0 * hp2 * (2.0 * hm1 - h0 - hp1) * alpha & ! Should cancel for case h = const - - (h0**2) * (3.0 * hm1 - h0 - 2.0 * hp1) * alpha & ! Should cancel for case h = const - - hp1 * (2.0 * hm1 * h0 * alpha - (h0 + hp1) * hp2 * beta) ! Should cancel for constant h & alpha = beta - - coeff_b_deltas = coeff_b_deltas & - / (3.0 * hm1 * h0 * ((h0 + hp1) / 2.0)) ! => -14.0 / 9.0 when h = const - - coeff_b_deltas = coeff_b_deltas & - / (2.0 * (h0 + hp1 + hp2) / 3.0) ! => -(14.0 / 9.0) / (2h) when h = const - end associate - end function coeff_b_deltas - - real function coeff_c_points(x) - !! Compute the coefficient acting on f_{i+2} 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. - - h = points_to_deltas(x) - coeff_c_points = coeff_c_deltas(h) - - end function coeff_c_points - - pure real function coeff_c_deltas(h) - !! Compute the coefficient acting on f_{i+2} 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 :: hm2, hm1, h0, hp1, hp2 ! Grid deltas at i -2, -1, 0, +1, +2 - - select type(deltas => h%stencil) - type is(real) - hm2 = deltas(-2) - 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) - coeff_c_deltas = -hm1 * h0 * hp1 - h0**2 * hp1 + hm1 * h0**2 * alpha & - + hm1 * h0 * hp1 * alpha + hm1 * h0 * hp1 * beta + h0**2 * hp1 * beta & ! End line 1 - + hm1 * hp1**2 * beta + 2.0 * h0 * hp1**2 * beta + hp1**3 * beta ! End line 2 - coeff_c_deltas = coeff_c_deltas & - / (hp2 * (hp1 + hp2) * (h0 + hp1 + hp2) * (hm1 + h0 + hp1 + hp2)) - end associate - end function coeff_c_deltas - - real function coeff_d_points(x) - !! Compute the coefficient acting on f_{i-2} 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. - - h = points_to_deltas(x) - coeff_d_points = coeff_d_deltas(h) - - end function coeff_d_points - - pure real function coeff_d_deltas(h) - !! Compute the coefficient acting on f_{i-2} 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 :: hm2, hm1, h0, hp1, hp2 ! Grid deltas at i -2, -1, 0, +1, +2 - - select type(deltas => h%stencil) - type is(real) - hm2 = deltas(-2) - 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) - coeff_d_deltas = h0 * hp1**2 + h0 * hp1 * hp2 - h0**3 * alpha - 2.0 * h0**2 * hp1 * alpha & - - h0 * hp1**2 * alpha - h0**2 * hp2 * alpha - h0 * hp1 * hp2 * alpha & ! End line 1 - - h0 * hp1 * hp2 * beta - hp1**2 * hp2 * beta ! End line 2 - coeff_d_deltas = coeff_d_deltas & - / (hm1 * (hm1 + h0) * (hm1 + h0 + hp1) * (hm1 + h0 + hp1 + hp2)) - end associate - end function coeff_d_deltas - - real function coeff_e(x) - !! Compute the coefficient acting on f_{i} of the finite diference given a stencil of points. - - type(nucfd_stencil_points), intent(in) :: x !! Stencil of points for the finite - !! difference. - - real :: a, b, c, d ! The neighbouring coefficients - - a = coeff_a(x) - b = coeff_b(x) - c = coeff_c(x) - d = coeff_d(x) - coeff_e = -(a + b + c + d) - - end function coeff_e end module nucfd_coeffs From 2652a2c2d0271da5c352301267e4b52fca1c0cb5 Mon Sep 17 00:00:00 2001 From: Paul Bartholomew Date: Tue, 28 Feb 2023 16:47:17 +0000 Subject: [PATCH 35/58] Split coeff_a into PURE subfunctions This should make it easier to debug as return values at intermediate stages can be checked, whereas intermediate values of a larger pure function cannot (atleast using print statements). --- src/coeffs/nucfd_coeffs_a.f90 | 85 ++++++++++++++++++++++++++++----- src/coeffs/nucfd_coeffs_mod.f90 | 2 +- 2 files changed, 75 insertions(+), 12 deletions(-) diff --git a/src/coeffs/nucfd_coeffs_a.f90 b/src/coeffs/nucfd_coeffs_a.f90 index add0723..924a920 100644 --- a/src/coeffs/nucfd_coeffs_a.f90 +++ b/src/coeffs/nucfd_coeffs_a.f90 @@ -1,4 +1,8 @@ submodule (nucfd_coeffs) nucfd_coeffs_a + !! Submodule defining the coefficient acting on f_{i+1} for compact finite difference schemes on + !! non-uniform grids. + !! + !! SPDX-License-Identifier: BSD-3-Clause implicit none @@ -17,7 +21,7 @@ module real function coeff_a_points(x) end function coeff_a_points - module pure real function coeff_a_deltas(h) + 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. @@ -37,19 +41,78 @@ module pure real function coeff_a_deltas(h) end select associate(beta => alpha) ! To match Gamet et al. (1999) - coeff_a_deltas = h0 * ((hm1 + h0) * (hp1 + hp2) + 2.0 * hp1 * hp2 * beta) & ! = (14/3) h^3 - + 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 - - coeff_a_deltas = coeff_a_deltas & - / (3.0 * hp1 * ((h0 + hp1) / 2.0) * hp2) ! => 14.0 / 9.0 when h = const + coeff_a_deltas = coeff_numerator(hm2, hm1, h0, hp1, hp2, alpha, beta) ! => (14/3) h^3, + ! h=const, alpha=beta + coeff_a_deltas = coeff_a_deltas & ! => Zero correction, + + coeff_numerator_corr(hm2, hm1, h0, hp1, hp2, alpha, beta) ! h=const, alpha=beta + coeff_a_deltas = coeff_a_deltas & + / coeff_denominator(h0, hp1, hp2) ! => 14/9, h=const coeff_a_deltas = coeff_a_deltas & - / (2.0 * ((hm1 + h0 + hp1) / 3.0)) ! => (14.0 / 9.0) / (2h) when h = const + / coeff_denominator_2h(hm1, h0, hp1) ! => (14/9)/(2h), h=const end associate end function coeff_a_deltas + + pure real function coeff_numerator(hm2, hm1, h0, hp1, hp2, alpha, beta) + !! Computes the numerator of the coefficient acting on f_{i+1}. + !! + !! Reduces to (14/3) h^3 when h=const and alpha=beta. + + real, intent(in) :: hm2 + 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 = h0 * ((hm1 + h0) * (hp1 + hp2) + 2.0 * hp1 * hp2 * beta) ! = (14/3) h^3 + end function coeff_numerator + + pure real function coeff_numerator_corr(hm2, 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) :: hm2 + 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. + + 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_denominator_2h(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. + + real, intent(in) :: hm1 + real, intent(in) :: h0 + real, intent(in) :: hp1 + + coeff_denominator_2h = (2.0 * ((hm1 + h0 + hp1) / 3.0)) + end function coeff_denominator_2h end submodule nucfd_coeffs_a diff --git a/src/coeffs/nucfd_coeffs_mod.f90 b/src/coeffs/nucfd_coeffs_mod.f90 index 2a76792..98b2ea8 100644 --- a/src/coeffs/nucfd_coeffs_mod.f90 +++ b/src/coeffs/nucfd_coeffs_mod.f90 @@ -43,7 +43,7 @@ module real function coeff_a_points(x) type(nucfd_stencil_points), intent(in) :: x !! Stencil of points for the finite !! difference. end function - module pure real function coeff_a_deltas(h) + module real function coeff_a_deltas(h) type(nucfd_stencil_deltas), intent(in) :: h !! Stencil of grid spacings for the finite !! difference. end function From e53685b965b34b88a3bbb28d9d66241ed6db23b8 Mon Sep 17 00:00:00 2001 From: Paul Bartholomew Date: Tue, 28 Feb 2023 17:27:52 +0000 Subject: [PATCH 36/58] Tidy up coeff_a submodule --- src/coeffs/nucfd_coeffs_a.f90 | 26 ++++++++++++-------------- 1 file changed, 12 insertions(+), 14 deletions(-) diff --git a/src/coeffs/nucfd_coeffs_a.f90 b/src/coeffs/nucfd_coeffs_a.f90 index 924a920..cbc9fee 100644 --- a/src/coeffs/nucfd_coeffs_a.f90 +++ b/src/coeffs/nucfd_coeffs_a.f90 @@ -42,39 +42,36 @@ module real function coeff_a_deltas(h) associate(beta => alpha) ! To match Gamet et al. (1999) - coeff_a_deltas = coeff_numerator(hm2, hm1, h0, hp1, hp2, alpha, beta) ! => (14/3) h^3, - ! h=const, alpha=beta - coeff_a_deltas = coeff_a_deltas & ! => Zero correction, - + coeff_numerator_corr(hm2, hm1, h0, hp1, hp2, alpha, beta) ! h=const, alpha=beta + coeff_a_deltas = coeff_numerator(hm1, h0, hp1, hp2, beta) ! => (14/3) h^3, + ! h=const, beta=1/3 + coeff_a_deltas = coeff_a_deltas & ! => Zero correction, + + coeff_numerator_corr(hm1, h0, hp1, hp2, alpha, beta) ! h=const, alpha=beta coeff_a_deltas = coeff_a_deltas & - / coeff_denominator(h0, hp1, hp2) ! => 14/9, h=const + / coeff_denominator(h0, hp1, hp2) ! => 14/9, h=const coeff_a_deltas = coeff_a_deltas & - / coeff_denominator_2h(hm1, h0, hp1) ! => (14/9)/(2h), h=const + / coeff_denominator_2h(hm1, h0, hp1) ! => (14/9)/(2h), h=const end associate end function coeff_a_deltas - pure real function coeff_numerator(hm2, hm1, h0, hp1, hp2, alpha, beta) + 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 alpha=beta. + !! Reduces to (14/3) h^3 when h=const and beta=1/3. - real, intent(in) :: hm2 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 = h0 * ((hm1 + h0) * (hp1 + hp2) + 2.0 * hp1 * hp2 * beta) ! = (14/3) h^3 end function coeff_numerator - pure real function coeff_numerator_corr(hm2, hm1, h0, hp1, hp2, alpha, beta) + 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) :: hm2 real, intent(in) :: hm1 real, intent(in) :: h0 real, intent(in) :: hp1 @@ -94,7 +91,7 @@ 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. + !! 14/9 when h=const, alpha=beta=1/3. real, intent(in) :: h0 real, intent(in) :: hp1 @@ -106,7 +103,8 @@ end function coeff_denominator pure real function coeff_denominator_2h(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. + !! 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 From a06ba88690342dc7d917c9d486cdf2b5db141246 Mon Sep 17 00:00:00 2001 From: Paul Bartholomew Date: Tue, 28 Feb 2023 17:28:15 +0000 Subject: [PATCH 37/58] Split coeff_b evaluation into smaller PURE subfunctions The larger calling function can then be debugged more easily --- src/coeffs/nucfd_coeffs_b.f90 | 83 ++++++++++++++++++++++++++++----- src/coeffs/nucfd_coeffs_mod.f90 | 2 +- 2 files changed, 73 insertions(+), 12 deletions(-) diff --git a/src/coeffs/nucfd_coeffs_b.f90 b/src/coeffs/nucfd_coeffs_b.f90 index fab288b..2a3b506 100644 --- a/src/coeffs/nucfd_coeffs_b.f90 +++ b/src/coeffs/nucfd_coeffs_b.f90 @@ -1,4 +1,8 @@ submodule (nucfd_coeffs) nucfd_coeffs_b + !! Submodule defining the coefficient acting on f_{i+1} for compact finite difference schemes on + !! non-uniform grids. + !! + !! SPDX-License-Identifier: BSD-3-Clause implicit none @@ -17,7 +21,7 @@ module real function coeff_b_points(x) end function coeff_b_points - module pure real function coeff_b_deltas(h) + module real function coeff_b_deltas(h) !! Compute the coefficient acting on f_{i-1} of the finite diference given a stencil of grid !! spacings. @@ -37,18 +41,75 @@ module pure real function coeff_b_deltas(h) end select associate(beta => alpha) ! To match Gamet et al. (1999) - coeff_b_deltas = -hp1 * ((hm1 + h0) * (hp1 + hp2) + 2.0 * hm1 * h0 * alpha) & ! = -(14/3) h^3 - - hm1 * hp1 * hp2 * (alpha - beta) & ! Should cancel for case alpha = beta - - (hp1**2) * (hm1 - h0) * alpha & ! Should cancel for case h = const - - h0 * hp2 * (2.0 * hm1 - h0 - hp1) * alpha & ! Should cancel for case h = const - - (h0**2) * (3.0 * hm1 - h0 - 2.0 * hp1) * alpha & ! Should cancel for case h = const - - hp1 * (2.0 * hm1 * h0 * alpha - (h0 + hp1) * hp2 * beta) ! Should cancel for constant h & alpha = beta - + coeff_b_deltas = coeff_numerator(hm1, h0, hp1, hp2, alpha) ! => (14/3) h^3, h=const, + ! alpha=1/3 + coeff_b_deltas = coeff_b_deltas & ! => Zero correction, + + coeff_numerator_corr(hm1, h0, hp1, hp2, alpha, beta) ! h=const, alpha=beta coeff_b_deltas = coeff_b_deltas & - / (3.0 * hm1 * h0 * ((h0 + hp1) / 2.0)) ! => -14.0 / 9.0 when h = const - + / coeff_denominator(hm1, h0, hp1) ! => -14/9, h=const coeff_b_deltas = coeff_b_deltas & - / (2.0 * (h0 + hp1 + hp2) / 3.0) ! => -(14.0 / 9.0) / (2h) when h = const + / coeff_denominator_2h(h0, hp1, hp2) ! => -(14/9)/(2h), h = const end associate end function coeff_b_deltas + + pure real function coeff_numerator(hm1, h0, hp1, hp2, alpha) + !! Computes the numerator of the coefficient acting on f_{i-1}. + !! + !! Reduces to -(14/3) h^3 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 + + coeff_numerator = -hp1 * ((hm1 + h0) * (hp1 + hp2) + 2.0 * hm1 * h0 * alpha) ! = -(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 * hp1 * hp2 * (alpha - beta) & ! Should cancel for case alpha = beta + - (hp1**2) * (hm1 - h0) * alpha & ! Should cancel for case h = const + - h0 * hp2 * (2.0 * hm1 - h0 - hp1) * alpha & ! Should cancel for case h = const + - (h0**2) * (3.0 * hm1 - h0 - 2.0 * hp1) * alpha & ! Should cancel for case h = const + - hp1 * (2.0 * hm1 * h0 * alpha - (h0 + hp1) * hp2 * beta) ! Should cancel for constant h & + ! alpha = beta + end function coeff_numerator_corr + + pure real function coeff_denominator(hm1, h0, hp1) + !! 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) :: hm1 + real, intent(in) :: h0 + real, intent(in) :: hp1 + + coeff_denominator = 3.0 * hm1 * h0 * ((h0 + hp1) / 2.0) + end function coeff_denominator + + pure real function coeff_denominator_2h(h0, hp1, hp2) + !! 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) :: h0 + real, intent(in) :: hp1 + real, intent(in) :: hp2 + + coeff_denominator_2h = 2.0 * (h0 + hp1 + hp2) / 3.0 + end function coeff_denominator_2h + end submodule nucfd_coeffs_b diff --git a/src/coeffs/nucfd_coeffs_mod.f90 b/src/coeffs/nucfd_coeffs_mod.f90 index 98b2ea8..7d4be4f 100644 --- a/src/coeffs/nucfd_coeffs_mod.f90 +++ b/src/coeffs/nucfd_coeffs_mod.f90 @@ -51,7 +51,7 @@ module real function coeff_b_points(x) type(nucfd_stencil_points), intent(in) :: x !! Stencil of points for the finite !! difference. end function - module pure real function coeff_b_deltas(h) + module real function coeff_b_deltas(h) type(nucfd_stencil_deltas), intent(in) :: h !! Stencil of grid spacings for the finite !! difference. end function From 46fdf9fa7cf7ebf022c5ccbe985f074c3566f3a9 Mon Sep 17 00:00:00 2001 From: Paul Bartholomew Date: Tue, 28 Feb 2023 19:32:18 +0000 Subject: [PATCH 38/58] Split coeff_c evaluation into smaller PURE functions Each of these functions can be tested independently --- src/coeffs/nucfd_coeffs_c.f90 | 85 ++++++++++++++++++++++++++---- src/coeffs/nucfd_coeffs_mod.f90 | 17 +++++- tests/fd-schemes/verify_coeffs.f90 | 11 ++++ 3 files changed, 102 insertions(+), 11 deletions(-) diff --git a/src/coeffs/nucfd_coeffs_c.f90 b/src/coeffs/nucfd_coeffs_c.f90 index a1616a9..87fb1ca 100644 --- a/src/coeffs/nucfd_coeffs_c.f90 +++ b/src/coeffs/nucfd_coeffs_c.f90 @@ -1,4 +1,8 @@ submodule (nucfd_coeffs) nucfd_coeffs_c + !! Submodule defining the coefficient acting on f_{i+2} for compact finite difference schemes on + !! non-uniform grids. + !! + !! SPDX-License-Identifier: BSD-3-Clause implicit none @@ -17,17 +21,35 @@ module real function coeff_c_points(x) end function coeff_c_points - module pure real function coeff_c_deltas(h) + module real function coeff_c_deltas(h) !! Compute the coefficient acting on f_{i+2} 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 :: hm2, hm1, h0, hp1, hp2 ! Grid deltas at i -2, -1, 0, +1, +2 + real :: numerator, denominator, divisor + + call coeff_c_components(h, numerator, denominator, divisor) + coeff_c_deltas = (numerator / denominator) / divisor + + end function coeff_c_deltas + + module subroutine coeff_c_components(h, numerator, denominator, divisor) + !! Compute the components of the coefficient acting on f_{i+2} of the finite diference given a + !! stencil of grid spacings. + !! + !! For uniform grids the numerator should reduce to (2/3) h^3, the denominator to 6h^3 (for a + !! coefficient value of (1/9) h^3) and the finite difference divisor to 4h. + + type(nucfd_stencil_deltas), intent(in) :: h !! Stencil of grid spacings for the finite difference. + real, intent(out) :: numerator !! The numerator of the coefficient + real, intent(out) :: denominator !! The denominator of the coefficient + real, intent(out) :: divisor !! The finite-difference divisor of the coefficient + + real :: hm1, h0, hp1, hp2 ! Grid deltas at i -1, 0, +1, +2 select type(deltas => h%stencil) type is(real) - hm2 = deltas(-2) hm1 = deltas(-1) h0 = deltas(0) hp1 = deltas(1) @@ -37,12 +59,57 @@ module pure real function coeff_c_deltas(h) end select associate(beta => alpha) ! To match Gamet et al. (1999) - coeff_c_deltas = -hm1 * h0 * hp1 - h0**2 * hp1 + hm1 * h0**2 * alpha & - + hm1 * h0 * hp1 * alpha + hm1 * h0 * hp1 * beta + h0**2 * hp1 * beta & ! End line 1 - + hm1 * hp1**2 * beta + 2.0 * h0 * hp1**2 * beta + hp1**3 * beta ! End line 2 - coeff_c_deltas = coeff_c_deltas & - / (hp2 * (hp1 + hp2) * (h0 + hp1 + hp2) * (hm1 + h0 + hp1 + hp2)) + numerator = coeff_numerator(hm1, h0, hp1, alpha, beta) + denominator = coeff_denominator(hm1, h0, hp1, hp2) + divisor = coeff_divisor(h0, hp1, hp2) end associate - end function coeff_c_deltas + + end subroutine coeff_c_components + + pure real function coeff_numerator(hm1, h0, hp1, alpha, beta) + !! Computes the numerator of the coefficient acting on f_{i+2}. + !! + !! Reduces to (2/3) h^3 when h=const and alpha=beta=1/3. + + real, intent(in) :: hm1 + real, intent(in) :: h0 + real, intent(in) :: hp1 + real, intent(in) :: alpha + real, intent(in) :: beta + + coeff_numerator = -h0 * hp1 * (hm1 + h0) & + + hm1 * h0 * (h0 + hp1) * alpha & + + hp1 * (h0 * (hm1 + h0) + hp1 * (hm1 + 2.0 * h0 + hp1)) * beta + + end function coeff_numerator + + pure real function coeff_denominator(hm1, h0, hp1, hp2) + !! Computes the denominator of the coefficient acting on f_{i+2}. + !! + !! Reduces to 6 h^3 when h=const. Dividing the numerator by this term yields the coefficient + !! 1/9 when h=const, alpha=beta=1/3. + + real, intent(in) :: hm1 + real, intent(in) :: h0 + real, intent(in) :: hp1 + real, intent(in) :: hp2 + + coeff_denominator = 3.0 * hp2 * (hp1 + hp2) & + * ((hm1 + h0 + hp1 + hp2) / 4.0) + + end function coeff_denominator + + pure real function coeff_divisor(h0, hp1, hp2) + !! Computes the non-uniform equivalent to 4h divisor of the coefficient acting on f_{i+2}. + !! + !! Dividing the coefficient by this term should reduce to (1/9)/(4h) when h=const, + !! alpha=beta=1/3. + + real, intent(in) :: h0 + real, intent(in) :: hp1 + real, intent(in) :: hp2 + + coeff_divisor = 4.0 * ((h0 + hp1 + hp2) / 3.0) + end function coeff_divisor end submodule nucfd_coeffs_c diff --git a/src/coeffs/nucfd_coeffs_mod.f90 b/src/coeffs/nucfd_coeffs_mod.f90 index 7d4be4f..5827828 100644 --- a/src/coeffs/nucfd_coeffs_mod.f90 +++ b/src/coeffs/nucfd_coeffs_mod.f90 @@ -10,7 +10,7 @@ module nucfd_coeffs private public :: coeff_a public :: coeff_b - public :: coeff_c + public :: coeff_c, coeff_c_components public :: coeff_d public :: coeff_e @@ -55,14 +55,27 @@ module real function coeff_b_deltas(h) type(nucfd_stencil_deltas), intent(in) :: h !! Stencil of grid spacings for the finite !! difference. end function + module real function coeff_c_points(x) type(nucfd_stencil_points), intent(in) :: x !! Stencil of points for the finite !! difference. end function - module pure real function coeff_c_deltas(h) + module real function coeff_c_deltas(h) type(nucfd_stencil_deltas), intent(in) :: h !! Stencil of grid spacings for the finite !! difference. end function + module subroutine coeff_c_components(h, numerator, denominator, divisor) + !! Compute the components of the coefficient acting on f_{i+2} of the finite diference given a + !! stencil of grid spacings. + !! + !! For uniform grids the numerator should reduce to (2/3) h^3, the denominator to 6h^3 (for a + !! coefficient value of (1/9) h^3) and the finite difference divisor to 4h. + type(nucfd_stencil_deltas), intent(in) :: h !! Stencil of grid spacings for the finite difference. + real, intent(out) :: numerator !! The numerator of the coefficient + real, intent(out) :: denominator !! The denominator of the coefficient + real, intent(out) :: divisor !! The finite-difference divisor of the coefficient + end subroutine coeff_c_components + module real function coeff_d_points(x) type(nucfd_stencil_points), intent(in) :: x !! Stencil of points for the finite !! difference. diff --git a/tests/fd-schemes/verify_coeffs.f90 b/tests/fd-schemes/verify_coeffs.f90 index ef63300..4f6fe4f 100644 --- a/tests/fd-schemes/verify_coeffs.f90 +++ b/tests/fd-schemes/verify_coeffs.f90 @@ -26,6 +26,11 @@ program verify_coeffs real, parameter :: cref = 1.0 / 9.0 real, parameter :: dref = -cref real, parameter :: eref = 0.0 + + real :: numerator, denominator, divisor + real, parameter :: numerator_f2ref = 2.0 / 3.0 + real, parameter :: denominator_f2ref = 6.0 + real, parameter :: divisor_f2ref = 4.0 call initialise_suite("Verify coefficients") @@ -51,8 +56,14 @@ program verify_coeffs call test_report("Coefficient A", check_scalar(a, aref / (2.0 * h))) b = coeff_b(stencil) call test_report("Coefficient B", check_scalar(b, bref / (2.0 * h))) + + call coeff_c_components(points_to_deltas(stencil), numerator, denominator, divisor) + call test_report("Coefficient C numerator", check_scalar(numerator, numerator_f2ref * (h**3))) + call test_report("Coefficient C denominator", check_scalar(denominator, denominator_f2ref * (h**3))) + call test_report("Coefficient C divisor", check_scalar(divisor, divisor_f2ref * h)) c = coeff_c(stencil) call test_report("Coefficient C", check_scalar(c, cref / (4.0 * h))) + d = coeff_d(stencil) call test_report("Coefficient D", check_scalar(d, dref / (4.0 * h))) e = coeff_e(stencil) From 7e5f51f98b5d0b7e8b383719486011ceb9652b70 Mon Sep 17 00:00:00 2001 From: Paul Bartholomew Date: Tue, 28 Feb 2023 20:02:40 +0000 Subject: [PATCH 39/58] Split coeff_d computation into smaller PURE functions --- src/coeffs/nucfd_coeffs_c.f90 | 2 +- src/coeffs/nucfd_coeffs_d.f90 | 86 ++++++++++++++++++++++++++---- src/coeffs/nucfd_coeffs_mod.f90 | 16 +++++- tests/fd-schemes/verify_coeffs.f90 | 5 ++ 4 files changed, 97 insertions(+), 12 deletions(-) diff --git a/src/coeffs/nucfd_coeffs_c.f90 b/src/coeffs/nucfd_coeffs_c.f90 index 87fb1ca..c726532 100644 --- a/src/coeffs/nucfd_coeffs_c.f90 +++ b/src/coeffs/nucfd_coeffs_c.f90 @@ -39,7 +39,7 @@ module subroutine coeff_c_components(h, numerator, denominator, divisor) !! stencil of grid spacings. !! !! For uniform grids the numerator should reduce to (2/3) h^3, the denominator to 6h^3 (for a - !! coefficient value of (1/9) h^3) and the finite difference divisor to 4h. + !! coefficient value of 1/9) and the finite difference divisor to 4h. type(nucfd_stencil_deltas), intent(in) :: h !! Stencil of grid spacings for the finite difference. real, intent(out) :: numerator !! The numerator of the coefficient diff --git a/src/coeffs/nucfd_coeffs_d.f90 b/src/coeffs/nucfd_coeffs_d.f90 index 72c4cfe..4b06556 100644 --- a/src/coeffs/nucfd_coeffs_d.f90 +++ b/src/coeffs/nucfd_coeffs_d.f90 @@ -1,4 +1,8 @@ submodule (nucfd_coeffs) nucfd_coeffs_d + !! Submodule defining the coefficient acting on f_{i+2} for compact finite difference schemes on + !! non-uniform grids. + !! + !! SPDX-License-Identifier: BSD-3-Clause implicit none @@ -17,17 +21,35 @@ module real function coeff_d_points(x) end function coeff_d_points - module pure real function coeff_d_deltas(h) + module real function coeff_d_deltas(h) !! Compute the coefficient acting on f_{i-2} 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 :: hm2, hm1, h0, hp1, hp2 ! Grid deltas at i -2, -1, 0, +1, +2 + real :: numerator, denominator, divisor + + call coeff_d_components(h, numerator, denominator, divisor) + coeff_d_deltas = (numerator / denominator) / divisor + + end function coeff_d_deltas + + module subroutine coeff_d_components(h, numerator, denominator, divisor) + !! Compute the components of the coefficient acting on f_{i-2} of the finite diference given a + !! stencil of grid spacings. + !! + !! For uniform grids the numerator should reduce to -(2/3) h^3, the denominator to 6h^3 (for a + !! coefficient value of -1/9) and the finite difference divisor to 4h. + + type(nucfd_stencil_deltas), intent(in) :: h !! Stencil of grid spacings for the finite difference. + real, intent(out) :: numerator !! The numerator of the coefficient + real, intent(out) :: denominator !! The denominator of the coefficient + real, intent(out) :: divisor !! The finite-difference divisor of the coefficient + + real :: hm1, h0, hp1, hp2 ! Grid deltas at i -1, 0, +1, +2 select type(deltas => h%stencil) type is(real) - hm2 = deltas(-2) hm1 = deltas(-1) h0 = deltas(0) hp1 = deltas(1) @@ -37,11 +59,57 @@ module pure real function coeff_d_deltas(h) end select associate(beta => alpha) ! To match Gamet et al. (1999) - coeff_d_deltas = h0 * hp1**2 + h0 * hp1 * hp2 - h0**3 * alpha - 2.0 * h0**2 * hp1 * alpha & - - h0 * hp1**2 * alpha - h0**2 * hp2 * alpha - h0 * hp1 * hp2 * alpha & ! End line 1 - - h0 * hp1 * hp2 * beta - hp1**2 * hp2 * beta ! End line 2 - coeff_d_deltas = coeff_d_deltas & - / (hm1 * (hm1 + h0) * (hm1 + h0 + hp1) * (hm1 + h0 + hp1 + hp2)) + numerator = coeff_numerator(h0, hp1, hp2, alpha, beta) + denominator = coeff_denominator(hm1, h0, hp1, hp2) + divisor = coeff_divisor(hm1, h0, hp1) end associate - end function coeff_d_deltas + + end subroutine coeff_d_components + + pure real function coeff_numerator(h0, hp1, hp2, alpha, beta) + !! Computes the numerator of the coefficient acting on f_{i-2}. + !! + !! Reduces to -(2/3) h^3 when h=const and alpha=beta=1/3. + + real, intent(in) :: h0 + real, intent(in) :: hp1 + real, intent(in) :: hp2 + real, intent(in) :: alpha + real, intent(in) :: beta + + coeff_numerator = h0 * hp1 * (hp1 + hp2) & + - h0 * (h0 * (h0 + 2.0 * hp1 + hp2) + hp1 * (hp1 + hp2)) * alpha & + - hp1 * hp2 * (h0 + hp1) * beta + + end function coeff_numerator + + pure real function coeff_denominator(hm1, h0, hp1, hp2) + !! Computes the denominator of the coefficient acting on f_{f-2}. + !! + !! Reduces to 6 h^3 when h=const. Dividing the numerator by this term yields the coefficient + !! -1/9 when h=const, alpha=beta=1/3. + + real, intent(in) :: hm1 + real, intent(in) :: h0 + real, intent(in) :: hp1 + real, intent(in) :: hp2 + + coeff_denominator = 3.0 * hm1 * (hp1 + hp2) & + * ((hm1 + h0 + hp1 + hp2) / 4.0) + + end function coeff_denominator + + pure real function coeff_divisor(hm1, h0, hp1) + !! Computes the non-uniform equivalent to 4h divisor of the coefficient acting on f_{i-2}. + !! + !! Dividing the coefficient by this term should reduce to -(1/9)/(4h) when h=const, + !! alpha=beta=1/3. + + real, intent(in) :: hm1 + real, intent(in) :: h0 + real, intent(in) :: hp1 + + coeff_divisor = 4.0 * ((hm1 + h0 + hp1) / 3.0) + end function coeff_divisor + end submodule nucfd_coeffs_d diff --git a/src/coeffs/nucfd_coeffs_mod.f90 b/src/coeffs/nucfd_coeffs_mod.f90 index 5827828..198e766 100644 --- a/src/coeffs/nucfd_coeffs_mod.f90 +++ b/src/coeffs/nucfd_coeffs_mod.f90 @@ -11,7 +11,7 @@ module nucfd_coeffs public :: coeff_a public :: coeff_b public :: coeff_c, coeff_c_components - public :: coeff_d + public :: coeff_d, coeff_d_components public :: coeff_e interface coeff_a @@ -80,10 +80,22 @@ module real function coeff_d_points(x) type(nucfd_stencil_points), intent(in) :: x !! Stencil of points for the finite !! difference. end function - module pure real function coeff_d_deltas(h) + module real function coeff_d_deltas(h) type(nucfd_stencil_deltas), intent(in) :: h !! Stencil of grid spacings for the finite !! difference. end function + module subroutine coeff_d_components(h, numerator, denominator, divisor) + !! Compute the components of the coefficient acting on f_{i-2} of the finite diference given a + !! stencil of grid spacings. + !! + !! For uniform grids the numerator should reduce to -(2/3) h^3, the denominator to 6h^3 (for a + !! coefficient value of -1/9) and the finite difference divisor to 4h. + type(nucfd_stencil_deltas), intent(in) :: h !! Stencil of grid spacings for the finite difference. + real, intent(out) :: numerator !! The numerator of the coefficient + real, intent(out) :: denominator !! The denominator of the coefficient + real, intent(out) :: divisor !! The finite-difference divisor of the coefficient + end subroutine coeff_d_components + module real function coeff_e(x) type(nucfd_stencil_points), intent(in) :: x !! Stencil of points for the finite !! difference. diff --git a/tests/fd-schemes/verify_coeffs.f90 b/tests/fd-schemes/verify_coeffs.f90 index 4f6fe4f..25ad49c 100644 --- a/tests/fd-schemes/verify_coeffs.f90 +++ b/tests/fd-schemes/verify_coeffs.f90 @@ -64,8 +64,13 @@ program verify_coeffs c = coeff_c(stencil) call test_report("Coefficient C", check_scalar(c, cref / (4.0 * h))) + call coeff_d_components(points_to_deltas(stencil), numerator, denominator, divisor) + call test_report("Coefficient C numerator", check_scalar(numerator, -numerator_f2ref * (h**3))) + call test_report("Coefficient C denominator", check_scalar(denominator, denominator_f2ref * (h**3))) + call test_report("Coefficient C divisor", check_scalar(divisor, divisor_f2ref * h)) d = coeff_d(stencil) call test_report("Coefficient D", check_scalar(d, dref / (4.0 * h))) + e = coeff_e(stencil) call test_report("Coefficient E", check_scalar(e, eref)) From 5e6770265010ca51b4feef5b155be90cdbc9bf5a Mon Sep 17 00:00:00 2001 From: Paul Bartholomew Date: Tue, 28 Feb 2023 20:17:02 +0000 Subject: [PATCH 40/58] Enable testing coefficient A by component --- src/coeffs/nucfd_coeffs_a.f90 | 38 ++++++++++++++++++++---------- src/coeffs/nucfd_coeffs_mod.f90 | 10 +++++++- tests/fd-schemes/verify_coeffs.f90 | 14 +++++++++-- 3 files changed, 46 insertions(+), 16 deletions(-) diff --git a/src/coeffs/nucfd_coeffs_a.f90 b/src/coeffs/nucfd_coeffs_a.f90 index cbc9fee..db62e48 100644 --- a/src/coeffs/nucfd_coeffs_a.f90 +++ b/src/coeffs/nucfd_coeffs_a.f90 @@ -27,6 +27,22 @@ module real function coeff_a_deltas(h) type(nucfd_stencil_deltas), intent(in) :: h !! Stencil of grid spacings for the finite difference. + real :: numerator, numerator_corr, denominator, divisor + + 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 :: hm2, hm1, h0, hp1, hp2 ! Grid deltas at i -2, -1, 0, +1, +2 select type(deltas => h%stencil) @@ -41,18 +57,14 @@ module real function coeff_a_deltas(h) end select associate(beta => alpha) ! To match Gamet et al. (1999) - - coeff_a_deltas = coeff_numerator(hm1, h0, hp1, hp2, beta) ! => (14/3) h^3, - ! h=const, beta=1/3 - coeff_a_deltas = coeff_a_deltas & ! => Zero correction, - + coeff_numerator_corr(hm1, h0, hp1, hp2, alpha, beta) ! h=const, alpha=beta - coeff_a_deltas = coeff_a_deltas & - / coeff_denominator(h0, hp1, hp2) ! => 14/9, h=const - coeff_a_deltas = coeff_a_deltas & - / coeff_denominator_2h(hm1, h0, hp1) ! => (14/9)/(2h), h=const + 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 function coeff_a_deltas + 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}. !! @@ -100,7 +112,7 @@ pure real function coeff_denominator(h0, hp1, hp2) coeff_denominator = (3.0 * hp1 * ((h0 + hp1) / 2.0) * hp2) end function coeff_denominator - pure real function coeff_denominator_2h(hm1, h0, hp1) + 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, @@ -110,7 +122,7 @@ pure real function coeff_denominator_2h(hm1, h0, hp1) real, intent(in) :: h0 real, intent(in) :: hp1 - coeff_denominator_2h = (2.0 * ((hm1 + h0 + hp1) / 3.0)) - end function coeff_denominator_2h + coeff_divisor = (2.0 * ((hm1 + h0 + hp1) / 3.0)) + end function coeff_divisor end submodule nucfd_coeffs_a diff --git a/src/coeffs/nucfd_coeffs_mod.f90 b/src/coeffs/nucfd_coeffs_mod.f90 index 198e766..5cd8493 100644 --- a/src/coeffs/nucfd_coeffs_mod.f90 +++ b/src/coeffs/nucfd_coeffs_mod.f90 @@ -8,7 +8,7 @@ module nucfd_coeffs implicit none private - public :: coeff_a + public :: coeff_a, coeff_a_components public :: coeff_b public :: coeff_c, coeff_c_components public :: coeff_d, coeff_d_components @@ -47,6 +47,14 @@ module real function coeff_a_deltas(h) type(nucfd_stencil_deltas), intent(in) :: h !! Stencil of grid spacings for the finite !! difference. end function + 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 + end subroutine coeff_a_components + module real function coeff_b_points(x) type(nucfd_stencil_points), intent(in) :: x !! Stencil of points for the finite !! difference. diff --git a/tests/fd-schemes/verify_coeffs.f90 b/tests/fd-schemes/verify_coeffs.f90 index 25ad49c..f36ad8f 100644 --- a/tests/fd-schemes/verify_coeffs.f90 +++ b/tests/fd-schemes/verify_coeffs.f90 @@ -27,7 +27,11 @@ program verify_coeffs real, parameter :: dref = -cref real, parameter :: eref = 0.0 - real :: numerator, denominator, divisor + real :: numerator, numerator_corr, denominator, divisor + real, parameter :: numerator_f1ref = 14.0 / 3.0 + real, parameter :: numerator_corr_f1ref = 0.0 + real, parameter :: denominator_f1ref = 3.0 + real, parameter :: divisor_f1ref = 2.0 real, parameter :: numerator_f2ref = 2.0 / 3.0 real, parameter :: denominator_f2ref = 6.0 real, parameter :: divisor_f2ref = 4.0 @@ -51,9 +55,15 @@ program verify_coeffs class default print *, "Error: Coordinate stencil is misallocated!" end select - + + call coeff_a_components(points_to_deltas(stencil), numerator, numerator_corr, denominator, divisor) + call test_report("Coefficient A numerator", check_scalar(numerator, numerator_f1ref * (h**3))) + call test_report("Coefficient A numerator correction", check_scalar(numerator_corr, numerator_corr_f1ref * (h**3))) + call test_report("Coefficient A denominator", check_scalar(denominator, denominator_f1ref * (h**3))) + call test_report("Coefficient A divisor", check_scalar(divisor, divisor_f1ref * h)) a = coeff_a(stencil) call test_report("Coefficient A", check_scalar(a, aref / (2.0 * h))) + b = coeff_b(stencil) call test_report("Coefficient B", check_scalar(b, bref / (2.0 * h))) From b3a4cf10624910e88ad941dc5c0804de0eb9f01a Mon Sep 17 00:00:00 2001 From: Paul Bartholomew Date: Tue, 28 Feb 2023 20:28:25 +0000 Subject: [PATCH 41/58] Enable testing coefficient B by component --- src/coeffs/nucfd_coeffs_b.f90 | 37 ++++++++++++++++++++---------- src/coeffs/nucfd_coeffs_mod.f90 | 9 +++++++- tests/fd-schemes/verify_coeffs.f90 | 5 ++++ 3 files changed, 38 insertions(+), 13 deletions(-) diff --git a/src/coeffs/nucfd_coeffs_b.f90 b/src/coeffs/nucfd_coeffs_b.f90 index 2a3b506..bcaa6fe 100644 --- a/src/coeffs/nucfd_coeffs_b.f90 +++ b/src/coeffs/nucfd_coeffs_b.f90 @@ -27,6 +27,22 @@ module real function coeff_b_deltas(h) type(nucfd_stencil_deltas), intent(in) :: h !! Stencil of grid spacings for the finite difference. + real :: numerator, numerator_corr, denominator, divisor + + call coeff_b_components(h, numerator, numerator_corr, denominator, divisor) + + coeff_b_deltas = ((numerator + numerator_corr) / denominator) / divisor + + end function coeff_b_deltas + + module subroutine coeff_b_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 :: hm2, hm1, h0, hp1, hp2 ! Grid deltas at i -2, -1, 0, +1, +2 select type(deltas => h%stencil) @@ -41,16 +57,13 @@ module real function coeff_b_deltas(h) end select associate(beta => alpha) ! To match Gamet et al. (1999) - coeff_b_deltas = coeff_numerator(hm1, h0, hp1, hp2, alpha) ! => (14/3) h^3, h=const, - ! alpha=1/3 - coeff_b_deltas = coeff_b_deltas & ! => Zero correction, - + coeff_numerator_corr(hm1, h0, hp1, hp2, alpha, beta) ! h=const, alpha=beta - coeff_b_deltas = coeff_b_deltas & - / coeff_denominator(hm1, h0, hp1) ! => -14/9, h=const - coeff_b_deltas = coeff_b_deltas & - / coeff_denominator_2h(h0, hp1, hp2) ! => -(14/9)/(2h), h = const + numerator = coeff_numerator(hm1, h0, hp1, hp2, alpha) + numerator_corr = coeff_numerator_corr(hm1, h0, hp1, hp2, alpha, beta) + denominator = coeff_denominator(hm1, h0, hp1) + divisor = coeff_divisor(h0, hp1, hp2) end associate - end function coeff_b_deltas + + end subroutine coeff_b_components pure real function coeff_numerator(hm1, h0, hp1, hp2, alpha) !! Computes the numerator of the coefficient acting on f_{i-1}. @@ -99,7 +112,7 @@ pure real function coeff_denominator(hm1, h0, hp1) coeff_denominator = 3.0 * hm1 * h0 * ((h0 + hp1) / 2.0) end function coeff_denominator - pure real function coeff_denominator_2h(h0, hp1, hp2) + pure real function coeff_divisor(h0, hp1, hp2) !! 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, @@ -109,7 +122,7 @@ pure real function coeff_denominator_2h(h0, hp1, hp2) real, intent(in) :: hp1 real, intent(in) :: hp2 - coeff_denominator_2h = 2.0 * (h0 + hp1 + hp2) / 3.0 - end function coeff_denominator_2h + coeff_divisor = 2.0 * (h0 + hp1 + hp2) / 3.0 + end function coeff_divisor end submodule nucfd_coeffs_b diff --git a/src/coeffs/nucfd_coeffs_mod.f90 b/src/coeffs/nucfd_coeffs_mod.f90 index 5cd8493..766bd84 100644 --- a/src/coeffs/nucfd_coeffs_mod.f90 +++ b/src/coeffs/nucfd_coeffs_mod.f90 @@ -9,7 +9,7 @@ module nucfd_coeffs private public :: coeff_a, coeff_a_components - public :: coeff_b + public :: coeff_b, coeff_b_components public :: coeff_c, coeff_c_components public :: coeff_d, coeff_d_components public :: coeff_e @@ -63,6 +63,13 @@ module real function coeff_b_deltas(h) type(nucfd_stencil_deltas), intent(in) :: h !! Stencil of grid spacings for the finite !! difference. end function + module subroutine coeff_b_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 + end subroutine coeff_b_components module real function coeff_c_points(x) type(nucfd_stencil_points), intent(in) :: x !! Stencil of points for the finite diff --git a/tests/fd-schemes/verify_coeffs.f90 b/tests/fd-schemes/verify_coeffs.f90 index f36ad8f..f128387 100644 --- a/tests/fd-schemes/verify_coeffs.f90 +++ b/tests/fd-schemes/verify_coeffs.f90 @@ -64,6 +64,11 @@ program verify_coeffs a = coeff_a(stencil) call test_report("Coefficient A", check_scalar(a, aref / (2.0 * h))) + call coeff_b_components(points_to_deltas(stencil), numerator, numerator_corr, denominator, divisor) + call test_report("Coefficient B numerator", check_scalar(numerator, -numerator_f1ref * (h**3))) + call test_report("Coefficient B numerator correction", check_scalar(numerator_corr, -numerator_corr_f1ref * (h**3))) + call test_report("Coefficient B denominator", check_scalar(denominator, denominator_f1ref * (h**3))) + call test_report("Coefficient B divisor", check_scalar(divisor, divisor_f1ref * h)) b = coeff_b(stencil) call test_report("Coefficient B", check_scalar(b, bref / (2.0 * h))) From e8c4aba6cc3d39630a929b4f5c26ad3e2851ced7 Mon Sep 17 00:00:00 2001 From: Paul Bartholomew Date: Tue, 28 Feb 2023 21:19:23 +0000 Subject: [PATCH 42/58] Fortran submodules require use of CMake v3.25.2 --- CMakeLists.txt | 2 +- README.md | 13 +++++++++++++ 2 files changed, 14 insertions(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 762f528..9912ca9 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -8,7 +8,7 @@ # ## Basic configuration -cmake_minimum_required(VERSION 3.16) +cmake_minimum_required(VERSION 3.25.2) project(NuCFD) enable_language(Fortran) diff --git a/README.md b/README.md index a59a942..2b2684c 100644 --- a/README.md +++ b/README.md @@ -15,6 +15,19 @@ 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 configuratoin 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 2da2353eb9d215a92e99a1d544c9e1d8c5f62aa1 Mon Sep 17 00:00:00 2001 From: Paul Bartholomew Date: Tue, 28 Feb 2023 21:58:28 +0000 Subject: [PATCH 43/58] Update derivative rules to test range of properties + functions --- src/nucfd_deriv_mod.f90 | 6 ++- tests/fd-schemes/differentiation_rules.f90 | 52 +++++++++++++++++++--- 2 files changed, 50 insertions(+), 8 deletions(-) diff --git a/src/nucfd_deriv_mod.f90 b/src/nucfd_deriv_mod.f90 index b5b11af..6a23df5 100644 --- a/src/nucfd_deriv_mod.f90 +++ b/src/nucfd_deriv_mod.f90 @@ -30,6 +30,7 @@ subroutine deriv_rhs(f, stencil, x, dfdx) call create_stencil(6, 4, stencil_coordinates) + print *, "*** Creating coordinate stencil ***" select type(indices => stencil%stencil) type is(integer) select type(points => stencil_coordinates%stencil) @@ -44,12 +45,15 @@ subroutine deriv_rhs(f, stencil, x, dfdx) print *, "Error: Coordinate stencil is misallocated!" error stop end select - + + print *, "*** Computing coefficients ***" a = coeff_a(stencil_coordinates) b = coeff_b(stencil_coordinates) c = coeff_c(stencil_coordinates) d = coeff_d(stencil_coordinates) e = coeff_e(stencil_coordinates) + + print *, "*** Assembling RHS ***" dfdx = (a * f(indices(+1)) + b * f(indices(-1))) & + (c * f(indices(+2)) + d * f(indices(-2))) & + e * f(indices(0)) diff --git a/tests/fd-schemes/differentiation_rules.f90 b/tests/fd-schemes/differentiation_rules.f90 index 5cac99b..795a896 100644 --- a/tests/fd-schemes/differentiation_rules.f90 +++ b/tests/fd-schemes/differentiation_rules.f90 @@ -37,11 +37,8 @@ program differentiation_rules allocate(x(n)) allocate(f(n)) - do i = 1, n - x(i) = real(i - 1) * h - f(i) = sin(x(i)) - end do - + + print *, "+++ Initialising stencil +++" call create_stencil(width, centre, stencil) i = 33 / 2 select type(indices => stencil%stencil) @@ -55,11 +52,52 @@ program differentiation_rules print *, "Error: Index stencil is misallocated!" error stop end select - + + print *, "+++ Testing derivative of constant function +++" + f(:) = 1.0 call deriv_rhs(f, stencil, x, dfdx) call deriv_rhs(f + 1.0, stencil, x, dgdx) + call test_report("Shifted derivative +, const f", check_scalar(dgdx, dfdx)) + call deriv_rhs(f - 1.0, stencil, x, dgdx) + call test_report("Shifted derivative -, const f", check_scalar(dgdx, dfdx)) + call deriv_rhs(f, stencil, x, dgdx) + call test_report("Shifted derivative 0, const f", check_scalar(dgdx, dfdx)) + call deriv_rhs(2.0 * f, stencil, x, dgdx) + call test_report("Scaled derivative 2x, const f", check_scalar(dgdx, 2.0 * dfdx)) + call deriv_rhs(-f, stencil, x, dgdx) + call test_report("Scaled derivative -1, const f", check_scalar(dgdx, -dfdx)) - call test_report("Shifted derivative", check_scalar(dgdx, dfdx)) + print *, "+++ Testing derivative of linearly increasing function +++" + f(1) = 0.0 + dfdx = 1.0 + do i = 2, n + f(i) = f(i - 1) + h * dfdx + end do + call deriv_rhs(f, stencil, x, dfdx) + call deriv_rhs(f + 1.0, stencil, x, dgdx) + call test_report("Shifted derivative +, linear+ f", check_scalar(dgdx, dfdx)) + call deriv_rhs(f - 1.0, stencil, x, dgdx) + call test_report("Shifted derivative -, linear+ f", check_scalar(dgdx, dfdx)) + call deriv_rhs(f, stencil, x, dgdx) + call test_report("Shifted derivative 0, linear+ f", check_scalar(dgdx, dfdx)) + call deriv_rhs(2.0 * f, stencil, x, dgdx) + call test_report("Scaled derivative 2x, linear+ f", check_scalar(dgdx, 2.0 * dfdx)) + call deriv_rhs(-f, stencil, x, dgdx) + call test_report("Scaled derivative -1, linear+ f", check_scalar(dgdx, -dfdx)) + + print *, "+++ Testing derivative of linearly decreasing function +++" + f(:) = -f(:) + call deriv_rhs(f, stencil, x, dfdx) + call deriv_rhs(f + 1.0, stencil, x, dgdx) + call test_report("Shifted derivative +, linear- f", check_scalar(dgdx, dfdx)) + call deriv_rhs(f - 1.0, stencil, x, dgdx) + call test_report("Shifted derivative -, linear- f", check_scalar(dgdx, dfdx)) + call deriv_rhs(f, stencil, x, dgdx) + call test_report("Shifted derivative 0, linear- f", check_scalar(dgdx, dfdx)) + call deriv_rhs(2.0 * f, stencil, x, dgdx) + call test_report("Scaled derivative 2x, linear- f", check_scalar(dgdx, 2.0 * dfdx)) + call deriv_rhs(-f, stencil, x, dgdx) + call test_report("Scaled derivative -1, linear- f", check_scalar(dgdx, -dfdx)) deallocate(x) deallocate(f) From c97b1380228b50b2805f9cf4d737d688eb97f44f Mon Sep 17 00:00:00 2001 From: Paul Bartholomew Date: Tue, 28 Feb 2023 22:31:31 +0000 Subject: [PATCH 44/58] Delta stencils were wider than necessary The grid stencil is symmetric (i-2 i-1 i i+1 i+2) - the delta stencil is not (i-1 i i+1 i+2) where h_k = x_k - x_k-1 (Gammet et al 1999) Added checking to coefficeints that stencil is correctly sized --- CMakeLists.txt | 6 ++-- src/coeffs/nucfd_coeffs_a.f90 | 45 +++++++++++++++++++++++++++-- src/coeffs/nucfd_coeffs_b.f90 | 46 ++++++++++++++++++++++++++++-- src/coeffs/nucfd_coeffs_c.f90 | 42 +++++++++++++++++++++++++++ src/coeffs/nucfd_coeffs_d.f90 | 42 +++++++++++++++++++++++++++ src/nucfd_deriv_mod.f90 | 3 +- src/nucfd_types_mod.f90 | 2 +- tests/fd-schemes/verify_coeffs.f90 | 5 ++-- 8 files changed, 177 insertions(+), 14 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 9912ca9..28f1244 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -23,9 +23,9 @@ set(icc_like_fc "$") set(cray_like_fc "$") add_library(nucfd_compiler_flags INTERFACE) target_compile_options(nucfd_compiler_flags INTERFACE - "$<${gcc_like_fc}:-std=f2018;-Wall;-Wextra;-Wpedantic;-Warray-bounds;-Wimplicit-interface;-Wimplicit-procedure;-fimplicit-none>" - "$<${icc_like_fc}:-std18;-warn;all>" - "$<${cray_like_fc}:-hfp0;-en;-emf;-eI>") # Anything greater than fp0 fails tests! + "$<${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>" diff --git a/src/coeffs/nucfd_coeffs_a.f90 b/src/coeffs/nucfd_coeffs_a.f90 index db62e48..877ff01 100644 --- a/src/coeffs/nucfd_coeffs_a.f90 +++ b/src/coeffs/nucfd_coeffs_a.f90 @@ -15,6 +15,20 @@ module real function coeff_a_points(x) !! 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(h%stencil, 1) + print *, size(h%stencil), lbound(h%stencil), ubound(h%stencil) + error stop + end if +#endif h = points_to_deltas(x) coeff_a_points = coeff_a(h) @@ -29,6 +43,20 @@ module real function coeff_a_deltas(h) 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 @@ -43,11 +71,24 @@ module subroutine coeff_a_components(h, numerator, numerator_corr, denominator, real, intent(out) :: denominator real, intent(out) :: divisor - real :: hm2, hm1, h0, hp1, hp2 ! Grid deltas at i -2, -1, 0, +1, +2 + 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) - hm2 = deltas(-2) hm1 = deltas(-1) h0 = deltas(0) hp1 = deltas(1) diff --git a/src/coeffs/nucfd_coeffs_b.f90 b/src/coeffs/nucfd_coeffs_b.f90 index bcaa6fe..308ed3c 100644 --- a/src/coeffs/nucfd_coeffs_b.f90 +++ b/src/coeffs/nucfd_coeffs_b.f90 @@ -15,6 +15,20 @@ module real function coeff_b_points(x) !! 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_b_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_b_points: expecting centre 3 stencil, received centre = ", & + 1 - lbound(h%stencil, 1) + print *, size(h%stencil), lbound(h%stencil), ubound(h%stencil) + error stop + end if +#endif h = points_to_deltas(x) coeff_b_points = coeff_b_deltas(h) @@ -29,8 +43,21 @@ module real function coeff_b_deltas(h) real :: numerator, numerator_corr, denominator, divisor - call coeff_b_components(h, numerator, numerator_corr, denominator, divisor) +#ifndef NDEBUG + if (size(h%stencil) /= 4) then + print *, "Error@coeff_b_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_b_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_b_components(h, numerator, numerator_corr, denominator, divisor) coeff_b_deltas = ((numerator + numerator_corr) / denominator) / divisor end function coeff_b_deltas @@ -43,11 +70,24 @@ module subroutine coeff_b_components(h, numerator, numerator_corr, denominator, real, intent(out) :: denominator real, intent(out) :: divisor - real :: hm2, hm1, h0, hp1, hp2 ! Grid deltas at i -2, -1, 0, +1, +2 + 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_b_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_b_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) - hm2 = deltas(-2) hm1 = deltas(-1) h0 = deltas(0) hp1 = deltas(1) diff --git a/src/coeffs/nucfd_coeffs_c.f90 b/src/coeffs/nucfd_coeffs_c.f90 index c726532..6e31ca2 100644 --- a/src/coeffs/nucfd_coeffs_c.f90 +++ b/src/coeffs/nucfd_coeffs_c.f90 @@ -15,6 +15,20 @@ module real function coeff_c_points(x) !! 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_c_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_c_points: expecting centre 3 stencil, received centre = ", & + 1 - lbound(h%stencil, 1) + print *, size(h%stencil), lbound(h%stencil), ubound(h%stencil) + error stop + end if +#endif h = points_to_deltas(x) coeff_c_points = coeff_c_deltas(h) @@ -28,6 +42,20 @@ module real function coeff_c_deltas(h) type(nucfd_stencil_deltas), intent(in) :: h !! Stencil of grid spacings for the finite difference. real :: numerator, denominator, divisor + +#ifndef NDEBUG + if (size(h%stencil) /= 4) then + print *, "Error@coeff_c_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_c_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_c_components(h, numerator, denominator, divisor) coeff_c_deltas = (numerator / denominator) / divisor @@ -48,6 +76,20 @@ module subroutine coeff_c_components(h, numerator, denominator, divisor) real :: hm1, h0, hp1, hp2 ! Grid deltas at i -1, 0, +1, +2 +#ifndef NDEBUG + if (size(h%stencil) /= 4) then + print *, "Error@coeff_c_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_c_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) diff --git a/src/coeffs/nucfd_coeffs_d.f90 b/src/coeffs/nucfd_coeffs_d.f90 index 4b06556..c5da23b 100644 --- a/src/coeffs/nucfd_coeffs_d.f90 +++ b/src/coeffs/nucfd_coeffs_d.f90 @@ -15,6 +15,20 @@ module real function coeff_d_points(x) !! 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_d_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_d_points: expecting centre 3 stencil, received centre = ", & + 1 - lbound(h%stencil, 1) + print *, size(h%stencil), lbound(h%stencil), ubound(h%stencil) + error stop + end if +#endif h = points_to_deltas(x) coeff_d_points = coeff_d_deltas(h) @@ -28,6 +42,20 @@ module real function coeff_d_deltas(h) type(nucfd_stencil_deltas), intent(in) :: h !! Stencil of grid spacings for the finite difference. real :: numerator, denominator, divisor + +#ifndef NDEBUG + if (size(h%stencil) /= 4) then + print *, "Error@coeff_d_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_d_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_d_components(h, numerator, denominator, divisor) coeff_d_deltas = (numerator / denominator) / divisor @@ -48,6 +76,20 @@ module subroutine coeff_d_components(h, numerator, denominator, divisor) real :: hm1, h0, hp1, hp2 ! Grid deltas at i -1, 0, +1, +2 +#ifndef NDEBUG + if (size(h%stencil) /= 4) then + print *, "Error@coeff_d_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_d_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) diff --git a/src/nucfd_deriv_mod.f90 b/src/nucfd_deriv_mod.f90 index 6a23df5..d3adcb2 100644 --- a/src/nucfd_deriv_mod.f90 +++ b/src/nucfd_deriv_mod.f90 @@ -28,14 +28,13 @@ subroutine deriv_rhs(f, stencil, x, dfdx) real :: a, b, c, d, e - call create_stencil(6, 4, stencil_coordinates) + call create_stencil(5, 3, stencil_coordinates) print *, "*** Creating coordinate stencil ***" select type(indices => stencil%stencil) type is(integer) select type(points => stencil_coordinates%stencil) type is(real) - points(-3) = x(indices(-2) - 1) points(-2) = x(indices(-2)) points(-1) = x(indices(-1)) points(+0) = x(indices(+0)) diff --git a/src/nucfd_types_mod.f90 b/src/nucfd_types_mod.f90 index 6e3f17a..ac09ee1 100644 --- a/src/nucfd_types_mod.f90 +++ b/src/nucfd_types_mod.f90 @@ -77,7 +77,7 @@ function points_to_deltas(x) result(h) type is(real) select type(deltas) type is(real) - do i = lbound(h%stencil, 1) + 1, ubound(h%stencil, 1) + do i = lb, ub xm1 = points(i - 1) x0 = points(i) diff --git a/tests/fd-schemes/verify_coeffs.f90 b/tests/fd-schemes/verify_coeffs.f90 index f128387..e27c298 100644 --- a/tests/fd-schemes/verify_coeffs.f90 +++ b/tests/fd-schemes/verify_coeffs.f90 @@ -42,11 +42,10 @@ program verify_coeffs L = 1.0 h = L / real(n - 1) - call create_stencil(6, 4, stencil) + call create_stencil(5, 3, stencil) select type(points => stencil%stencil) type is(real) points(:) = 0.0 - points(-3) = -3.0 * h points(-2) = -2.0 * h points(-1) = -1.0 * h points(0) = 0.0 @@ -91,7 +90,7 @@ program verify_coeffs select type(points => stencil%stencil) type is(real) - points(-3) = -5.0 * h + points(-2) = -4.0 * h points(2) = +3.0 * h class default print *, "Error: Coordinate stencil is misallocated!" From fc47ba87a01f84780d40a061f97fd527c9603b99 Mon Sep 17 00:00:00 2001 From: Paul Bartholomew Date: Tue, 28 Feb 2023 22:49:13 +0000 Subject: [PATCH 45/58] Need to actually set the grid in the test... --- src/nucfd_deriv_mod.f90 | 3 --- tests/fd-schemes/differentiation_rules.f90 | 3 +++ tests/nucfd_tests_mod.f90 | 1 + 3 files changed, 4 insertions(+), 3 deletions(-) diff --git a/src/nucfd_deriv_mod.f90 b/src/nucfd_deriv_mod.f90 index d3adcb2..df3cdd3 100644 --- a/src/nucfd_deriv_mod.f90 +++ b/src/nucfd_deriv_mod.f90 @@ -30,7 +30,6 @@ subroutine deriv_rhs(f, stencil, x, dfdx) call create_stencil(5, 3, stencil_coordinates) - print *, "*** Creating coordinate stencil ***" select type(indices => stencil%stencil) type is(integer) select type(points => stencil_coordinates%stencil) @@ -45,14 +44,12 @@ subroutine deriv_rhs(f, stencil, x, dfdx) error stop end select - print *, "*** Computing coefficients ***" a = coeff_a(stencil_coordinates) b = coeff_b(stencil_coordinates) c = coeff_c(stencil_coordinates) d = coeff_d(stencil_coordinates) e = coeff_e(stencil_coordinates) - print *, "*** Assembling RHS ***" dfdx = (a * f(indices(+1)) + b * f(indices(-1))) & + (c * f(indices(+2)) + d * f(indices(-2))) & + e * f(indices(0)) diff --git a/tests/fd-schemes/differentiation_rules.f90 b/tests/fd-schemes/differentiation_rules.f90 index 795a896..f5d507f 100644 --- a/tests/fd-schemes/differentiation_rules.f90 +++ b/tests/fd-schemes/differentiation_rules.f90 @@ -37,6 +37,9 @@ program differentiation_rules allocate(x(n)) allocate(f(n)) + do i = 1, n + x(i) = real(i - 1) * h + end do print *, "+++ Initialising stencil +++" call create_stencil(width, centre, stencil) diff --git a/tests/nucfd_tests_mod.f90 b/tests/nucfd_tests_mod.f90 index 85ba249..a40122e 100644 --- a/tests/nucfd_tests_mod.f90 +++ b/tests/nucfd_tests_mod.f90 @@ -82,6 +82,7 @@ logical function check_scalar(test, ref) test_passing = .false. print *, "Delta = ", err, " exceeds tolerance: ", 2 * epsilon(ref) + print *, "Value = ", test, " expected: ", ref else test_passing = .true. end if From 6c8b1c5e85793015b134d1144213fea44a52eb56 Mon Sep 17 00:00:00 2001 From: Paul Bartholomew Date: Tue, 28 Feb 2023 22:56:16 +0000 Subject: [PATCH 46/58] Compute finite difference in floating-point robust way --- src/nucfd_deriv_mod.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/nucfd_deriv_mod.f90 b/src/nucfd_deriv_mod.f90 index df3cdd3..aa68e7a 100644 --- a/src/nucfd_deriv_mod.f90 +++ b/src/nucfd_deriv_mod.f90 @@ -50,8 +50,8 @@ subroutine deriv_rhs(f, stencil, x, dfdx) d = coeff_d(stencil_coordinates) e = coeff_e(stencil_coordinates) - dfdx = (a * f(indices(+1)) + b * f(indices(-1))) & - + (c * f(indices(+2)) + d * f(indices(-2))) & + dfdx = a * (f(indices(+1)) + (b / a) * f(indices(-1))) & + + c * (f(indices(+2)) + (d / c) * f(indices(-2))) & + e * f(indices(0)) class default print *, "Error: Index stencil is misallocated!" From 1832cb5b1f36dc0e5ac2f0596a7d01a6788d9bc0 Mon Sep 17 00:00:00 2001 From: Paul Bartholomew Date: Fri, 3 Mar 2023 18:34:12 +0000 Subject: [PATCH 47/58] Document the coeff_e module --- src/coeffs/nucfd_coeffs_e.f90 | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/src/coeffs/nucfd_coeffs_e.f90 b/src/coeffs/nucfd_coeffs_e.f90 index 450beb5..f1e33d0 100644 --- a/src/coeffs/nucfd_coeffs_e.f90 +++ b/src/coeffs/nucfd_coeffs_e.f90 @@ -1,4 +1,8 @@ submodule (nucfd_coeffs) nucfd_coeffs_e + !! Submodule defining the coefficient acting on f_{i} for compact finite difference schemes on + !! non-uniform grids. + !! + !! SPDX-License-Identifier: BSD-3-Clause implicit none @@ -12,6 +16,20 @@ module real function coeff_e(x) real :: a, b, c, d ! The neighbouring coefficients +#ifndef NDEBUG + if (size(x%stencil) /= 5) then + print *, "Error@coeff_e: expecting width 5 stencil, received width = ", & + size(x%stencil) + error stop + end if + if (1 - lbound(x%stencil, 1) /= 3) then + print *, "Error@coeff_e: 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 + a = coeff_a(x) b = coeff_b(x) c = coeff_c(x) From 296abaf72e0cae3cc8045311824116a8570465e5 Mon Sep 17 00:00:00 2001 From: Paul Bartholomew Date: Fri, 3 Mar 2023 18:35:54 +0000 Subject: [PATCH 48/58] Correct copy-paste error --- src/coeffs/nucfd_coeffs_a.f90 | 4 ++-- src/coeffs/nucfd_coeffs_b.f90 | 4 ++-- src/coeffs/nucfd_coeffs_c.f90 | 4 ++-- src/coeffs/nucfd_coeffs_d.f90 | 4 ++-- 4 files changed, 8 insertions(+), 8 deletions(-) diff --git a/src/coeffs/nucfd_coeffs_a.f90 b/src/coeffs/nucfd_coeffs_a.f90 index 877ff01..78ee436 100644 --- a/src/coeffs/nucfd_coeffs_a.f90 +++ b/src/coeffs/nucfd_coeffs_a.f90 @@ -24,8 +24,8 @@ module real function coeff_a_points(x) end if if (1 - lbound(x%stencil, 1) /= 3) then print *, "Error@coeff_a_points: expecting centre 3 stencil, received centre = ", & - 1 - lbound(h%stencil, 1) - print *, size(h%stencil), lbound(h%stencil), ubound(h%stencil) + 1 - lbound(x%stencil, 1) + print *, size(x%stencil), lbound(x%stencil), ubound(x%stencil) error stop end if #endif diff --git a/src/coeffs/nucfd_coeffs_b.f90 b/src/coeffs/nucfd_coeffs_b.f90 index 308ed3c..9249419 100644 --- a/src/coeffs/nucfd_coeffs_b.f90 +++ b/src/coeffs/nucfd_coeffs_b.f90 @@ -24,8 +24,8 @@ module real function coeff_b_points(x) end if if (1 - lbound(x%stencil, 1) /= 3) then print *, "Error@coeff_b_points: expecting centre 3 stencil, received centre = ", & - 1 - lbound(h%stencil, 1) - print *, size(h%stencil), lbound(h%stencil), ubound(h%stencil) + 1 - lbound(x%stencil, 1) + print *, size(x%stencil), lbound(x%stencil), ubound(x%stencil) error stop end if #endif diff --git a/src/coeffs/nucfd_coeffs_c.f90 b/src/coeffs/nucfd_coeffs_c.f90 index 6e31ca2..30e7012 100644 --- a/src/coeffs/nucfd_coeffs_c.f90 +++ b/src/coeffs/nucfd_coeffs_c.f90 @@ -24,8 +24,8 @@ module real function coeff_c_points(x) end if if (1 - lbound(x%stencil, 1) /= 3) then print *, "Error@coeff_c_points: expecting centre 3 stencil, received centre = ", & - 1 - lbound(h%stencil, 1) - print *, size(h%stencil), lbound(h%stencil), ubound(h%stencil) + 1 - lbound(x%stencil, 1) + print *, size(x%stencil), lbound(x%stencil), ubound(x%stencil) error stop end if #endif diff --git a/src/coeffs/nucfd_coeffs_d.f90 b/src/coeffs/nucfd_coeffs_d.f90 index c5da23b..77013a0 100644 --- a/src/coeffs/nucfd_coeffs_d.f90 +++ b/src/coeffs/nucfd_coeffs_d.f90 @@ -24,8 +24,8 @@ module real function coeff_d_points(x) end if if (1 - lbound(x%stencil, 1) /= 3) then print *, "Error@coeff_d_points: expecting centre 3 stencil, received centre = ", & - 1 - lbound(h%stencil, 1) - print *, size(h%stencil), lbound(h%stencil), ubound(h%stencil) + 1 - lbound(x%stencil, 1) + print *, size(x%stencil), lbound(x%stencil), ubound(x%stencil) error stop end if #endif From cdfb9fcd7f9ac9033dfad350f82b1c992cb02ddd Mon Sep 17 00:00:00 2001 From: Paul Date: Sat, 4 Mar 2023 19:55:36 +0000 Subject: [PATCH 49/58] Fix code blocks in README --- README.md | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/README.md b/README.md index 2b2684c..ec573b7 100644 --- a/README.md +++ b/README.md @@ -4,14 +4,14 @@ 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. @@ -20,32 +20,32 @@ 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 configuratoin required +# 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 -`` +``` From 30a26a6018cd80756cc4f08061b37fd58b184844 Mon Sep 17 00:00:00 2001 From: Paul Date: Sat, 4 Mar 2023 20:10:10 +0000 Subject: [PATCH 50/58] Ensuring each source file has a documentation header --- src/coeffs/nucfd_coeffs_a.f90 | 8 ++++++-- src/coeffs/nucfd_coeffs_b.f90 | 8 ++++++-- src/coeffs/nucfd_coeffs_c.f90 | 8 ++++++-- src/coeffs/nucfd_coeffs_d.f90 | 8 ++++++-- src/coeffs/nucfd_coeffs_e.f90 | 8 ++++++++ src/coeffs/nucfd_coeffs_mod.f90 | 8 ++++++-- src/nucfd_trid_solver_mod.f90 | 8 ++++++-- src/nucfd_types_mod.f90 | 8 ++++++-- tests/nucfd_tests_mod.f90 | 8 ++++++-- tests/tridsolver/system_00_anti-symm.f90 | 12 ++++++++---- tests/tridsolver/system_00_symm.f90 | 12 ++++++++---- tests/tridsolver/system_11_anti-symm.f90 | 12 ++++++++---- tests/tridsolver/system_11_symm.f90 | 12 ++++++++---- tests/tridsolver/tridsol_test_utils_mod.f90 | 12 ++++++++---- 14 files changed, 96 insertions(+), 36 deletions(-) diff --git a/src/coeffs/nucfd_coeffs_a.f90 b/src/coeffs/nucfd_coeffs_a.f90 index 877ff01..98111e9 100644 --- a/src/coeffs/nucfd_coeffs_a.f90 +++ b/src/coeffs/nucfd_coeffs_a.f90 @@ -1,8 +1,12 @@ +! 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. - !! - !! SPDX-License-Identifier: BSD-3-Clause implicit none diff --git a/src/coeffs/nucfd_coeffs_b.f90 b/src/coeffs/nucfd_coeffs_b.f90 index 308ed3c..0a1cee1 100644 --- a/src/coeffs/nucfd_coeffs_b.f90 +++ b/src/coeffs/nucfd_coeffs_b.f90 @@ -1,8 +1,12 @@ +! src/coeffs/nucfd_coeff_b.f90 +! +!! Implements coefficient B of the nucfd_coeffs module. +! +! SPDX-License-Identifier: BSD-3-Clause + submodule (nucfd_coeffs) nucfd_coeffs_b !! Submodule defining the coefficient acting on f_{i+1} for compact finite difference schemes on !! non-uniform grids. - !! - !! SPDX-License-Identifier: BSD-3-Clause implicit none diff --git a/src/coeffs/nucfd_coeffs_c.f90 b/src/coeffs/nucfd_coeffs_c.f90 index 6e31ca2..3843cef 100644 --- a/src/coeffs/nucfd_coeffs_c.f90 +++ b/src/coeffs/nucfd_coeffs_c.f90 @@ -1,8 +1,12 @@ +! src/coeffs/nucfd_coeff_c.f90 +! +!! Implements coefficient C of the nucfd_coeffs module. +! +! SPDX-License-Identifier: BSD-3-Clause + submodule (nucfd_coeffs) nucfd_coeffs_c !! Submodule defining the coefficient acting on f_{i+2} for compact finite difference schemes on !! non-uniform grids. - !! - !! SPDX-License-Identifier: BSD-3-Clause implicit none diff --git a/src/coeffs/nucfd_coeffs_d.f90 b/src/coeffs/nucfd_coeffs_d.f90 index c5da23b..37e673e 100644 --- a/src/coeffs/nucfd_coeffs_d.f90 +++ b/src/coeffs/nucfd_coeffs_d.f90 @@ -1,8 +1,12 @@ +! src/coeffs/nucfd_coeff_d.f90 +! +!! Implements coefficient D of the nucfd_coeffs module. +! +! SPDX-License-Identifier: BSD-3-Clause + submodule (nucfd_coeffs) nucfd_coeffs_d !! Submodule defining the coefficient acting on f_{i+2} for compact finite difference schemes on !! non-uniform grids. - !! - !! SPDX-License-Identifier: BSD-3-Clause implicit none diff --git a/src/coeffs/nucfd_coeffs_e.f90 b/src/coeffs/nucfd_coeffs_e.f90 index 450beb5..356180e 100644 --- a/src/coeffs/nucfd_coeffs_e.f90 +++ b/src/coeffs/nucfd_coeffs_e.f90 @@ -1,4 +1,12 @@ +! src/coeffs/nucfd_coeff_e.f90 +! +!! Implements coefficient E of the nucfd_coeffs module. +! +! SPDX-License-Identifier: BSD-3-Clause + submodule (nucfd_coeffs) nucfd_coeffs_e + !! Submodule defining the coefficient acting on f_i for compact finite difference schemes on + !! non-uniform grids. implicit none diff --git a/src/coeffs/nucfd_coeffs_mod.f90 b/src/coeffs/nucfd_coeffs_mod.f90 index 766bd84..6c57d55 100644 --- a/src/coeffs/nucfd_coeffs_mod.f90 +++ b/src/coeffs/nucfd_coeffs_mod.f90 @@ -1,7 +1,11 @@ +! src/coeffs/nucfd_coeff_mod.f90 +! +!! Defines the nucfd_coeffs module. +! +! SPDX-License-Identifier: BSD-3-Clause + module nucfd_coeffs !! Module defining the coefficients for compact finite difference schemes on non-uniform grids. - !! - !! SPDX-License-Identifier: BSD-3-Clause use nucfd_types diff --git a/src/nucfd_trid_solver_mod.f90 b/src/nucfd_trid_solver_mod.f90 index be2f8e7..eda61dc 100644 --- a/src/nucfd_trid_solver_mod.f90 +++ b/src/nucfd_trid_solver_mod.f90 @@ -1,9 +1,13 @@ +! src/nucfd_trid_solver_mod.f90 +! +!! Defines the nucfd_trid_solver module. +! +! SPDX-License-Identifier: BSD-3-Clause + module nucfd_trid_solver !! Module implementing a simple tridiagonal solver based on the algorithm given at !! https://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm including support for cyclic !! systems. - !! - !! SPDX-License-Identifier: BSD-3-Clause implicit none diff --git a/src/nucfd_types_mod.f90 b/src/nucfd_types_mod.f90 index ac09ee1..1c2bace 100644 --- a/src/nucfd_types_mod.f90 +++ b/src/nucfd_types_mod.f90 @@ -1,7 +1,11 @@ +! src/nucfd_types_mod.f90 +! +!! Defines the nucfd_types module. +! +! SPDX-License-Identifier: BSD-3-Clause + module nucfd_types !! Module defining the custom types used by NuCFD. - !! - !! SPDX-License-Identifier: BSD-3-Clause implicit none diff --git a/tests/nucfd_tests_mod.f90 b/tests/nucfd_tests_mod.f90 index a40122e..cc783fa 100644 --- a/tests/nucfd_tests_mod.f90 +++ b/tests/nucfd_tests_mod.f90 @@ -1,8 +1,12 @@ +! tests/nucfd_tests_mod.f90 +! +!! Defines the nucfd_tests module. +! +! SPDX-License-Identifier: BSD-3-Clause + module nucfd_tests !! NuCFD test module, enables defining test suites, reporting results of tests and overall status !! of the suite, and provides utilities for checking floating-point numbers. - !! - !! SPDX-License-Identifier: BSD-3-Clause implicit none diff --git a/tests/tridsolver/system_00_anti-symm.f90 b/tests/tridsolver/system_00_anti-symm.f90 index 808e3a7..c1642ff 100644 --- a/tests/tridsolver/system_00_anti-symm.f90 +++ b/tests/tridsolver/system_00_anti-symm.f90 @@ -1,10 +1,14 @@ +! tests/tridsolver/system_00_anti-symm.f90 +! +!! Defines the test for solving tridiagonal systems of anti-symmetric functions on periodic domains. +!! +!! Part of the tridsolver test suite. +! +! SPDX-License-Identifier: BSD-3-Clause + program test_system_00_anti_symm !! Tests the solution of tridiagonal systems arising from compact finite differences of !! anti-symmetric functions on periodic domains. - !! - !! Part of the tridsolver test suite. - !! - !! SPDX-License-Identifier: BSD-3-Clause use nucfd_coeffs use nucfd_trid_solver diff --git a/tests/tridsolver/system_00_symm.f90 b/tests/tridsolver/system_00_symm.f90 index c89abbc..85aa342 100644 --- a/tests/tridsolver/system_00_symm.f90 +++ b/tests/tridsolver/system_00_symm.f90 @@ -1,10 +1,14 @@ +! tests/tridsolver/system_00_symm.f90 +! +!! Defines the test for solving tridiagonal systems of symmetric functions on periodic domains. +!! +!! Part of the tridsolver test suite. +! +! SPDX-License-Identifier: BSD-3-Clause + program test_system_00_symm !! Tests the solution of tridiagonal systems arising from compact finite differences of symmetric !! functions on periodic domains. - !! - !! Part of the tridsolver test suite. - !! - !! SPDX-License-Identifier: BSD-3-Clause use nucfd_coeffs use nucfd_trid_solver diff --git a/tests/tridsolver/system_11_anti-symm.f90 b/tests/tridsolver/system_11_anti-symm.f90 index 80e2563..e8a3804 100644 --- a/tests/tridsolver/system_11_anti-symm.f90 +++ b/tests/tridsolver/system_11_anti-symm.f90 @@ -1,10 +1,14 @@ +! tests/tridsolver/system_11_anti-symm.f90 +! +!! Defines the test for solving tridiagonal systems of anti-symmetric functions. +!! +!! Part of the tridsolver test suite. +! +! SPDX-License-Identifier: BSD-3-Clause + program test_system_11_anti_symm !! Tests the solution of tridiagonal systems arising from compact finite differences of !! anti-symmetric functions. - !! - !! Part of the tridsolver test suite. - !! - !! SPDX-License-Identifier: BSD-3-Clause use nucfd_coeffs use nucfd_trid_solver diff --git a/tests/tridsolver/system_11_symm.f90 b/tests/tridsolver/system_11_symm.f90 index 77b2367..603a110 100644 --- a/tests/tridsolver/system_11_symm.f90 +++ b/tests/tridsolver/system_11_symm.f90 @@ -1,10 +1,14 @@ +! tests/tridsolver/system_11_symm.f90 +! +!! Defines the test for solving tridiagonal systems of symmetric functions. +!! +!! Part of the tridsolver test suite. +! +! SPDX-License-Identifier: BSD-3-Clause + program test_system_11_symm !! Tests the solution of tridiagonal systems arising from compact finite differences of symmetric !! functions. - !! - !! Part of the tridsolver test suite. - !! - !! SPDX-License-Identifier: BSD-3-Clause use nucfd_coeffs use nucfd_trid_solver diff --git a/tests/tridsolver/tridsol_test_utils_mod.f90 b/tests/tridsolver/tridsol_test_utils_mod.f90 index 636fb8d..8b58ed1 100644 --- a/tests/tridsolver/tridsol_test_utils_mod.f90 +++ b/tests/tridsolver/tridsol_test_utils_mod.f90 @@ -1,9 +1,13 @@ +! tests/tridsolver/tridsol_test_utils_mod.f90 +! +!! Defines the tridsol_test_utils module. +!! +!! Part of the tridsolver test suite. +! +! SPDX-License-Identifier: BSD-3-Clause + module tridsol_test_utils !! Module defining utility functions for the tridiagonal solver test suite. - !! - !! Part of the tridsolver test suite. - !! - !! SPDX-License-Identifier: BSD-3-Clause implicit none From b89a55f0eddacb31bbe7a8b687b129a679e20070 Mon Sep 17 00:00:00 2001 From: Paul Date: Sat, 4 Mar 2023 20:17:05 +0000 Subject: [PATCH 51/58] Documenting procedures --- src/nucfd_deriv_mod.f90 | 1 + 1 file changed, 1 insertion(+) diff --git a/src/nucfd_deriv_mod.f90 b/src/nucfd_deriv_mod.f90 index aa68e7a..f9a946b 100644 --- a/src/nucfd_deriv_mod.f90 +++ b/src/nucfd_deriv_mod.f90 @@ -18,6 +18,7 @@ module nucfd_deriv contains subroutine deriv_rhs(f, stencil, x, dfdx) + !! Compute the RHS of the derivative of a function f based on a stencil of indices. real, dimension(:), intent(in) :: f type(nucfd_index_stencil), intent(in) :: stencil From f99e95a1dd983fb964dbf142782ee41bd47d2e04 Mon Sep 17 00:00:00 2001 From: Paul Date: Sat, 4 Mar 2023 23:15:37 +0000 Subject: [PATCH 52/58] Tidy up CMakeLists.txt --- CMakeLists.txt | 53 +++++---------------------------- src/CMakeLists.txt | 24 +++++++++++++++ src/coeffs/CMakeLists.txt | 21 +++++++++++++ tests/CMakeLists.txt | 37 +++++++++++++++++++++++ tests/fd-schemes/CMakeLists.txt | 13 ++++++++ tests/tridsolver/CMakeLists.txt | 23 ++++++++++++++ 6 files changed, 125 insertions(+), 46 deletions(-) create mode 100644 src/CMakeLists.txt create mode 100644 src/coeffs/CMakeLists.txt create mode 100644 tests/CMakeLists.txt create mode 100644 tests/fd-schemes/CMakeLists.txt create mode 100644 tests/tridsolver/CMakeLists.txt diff --git a/CMakeLists.txt b/CMakeLists.txt index 28f1244..82a2cde 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -2,6 +2,8 @@ # ## Description # +# Main CMakeLists file for the NuCFD project. +# ## LICENSE # # SPDX-License-Identifier: BSD-3-Clause @@ -35,57 +37,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>") - -## Build nucfd library -add_library(nucfd - src/nucfd_trid_solver_mod.f90 - src/nucfd_types_mod.f90 - src/coeffs/nucfd_coeffs_mod.f90 src/coeffs/nucfd_coeffs_a.f90 src/coeffs/nucfd_coeffs_b.f90 src/coeffs/nucfd_coeffs_c.f90 src/coeffs/nucfd_coeffs_d.f90 src/coeffs/nucfd_coeffs_e.f90 - src/nucfd_deriv_mod.f90) -target_link_libraries(nucfd nucfd_compiler_flags) if (CMAKE_BUILD_TYPE MATCHES "Debug") - target_link_libraries(nucfd nucfd_debug_flags) + target_link_libraries(nucfd_compiler_flags nucfd_debug_flags) endif() +## Build nucfd library +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) - -# Finite difference test suite -define_test(fd-schemes verify_coeffs) -define_test(fd-schemes differentiation_rules) +add_subdirectory(tests) ## Documentation add_custom_target(doc ford diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt new file mode 100644 index 0000000..9efc08b --- /dev/null +++ b/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_link_libraries(nucfd PRIVATE $) +target_link_libraries(nucfd PUBLIC nucfd_compiler_flags) diff --git a/src/coeffs/CMakeLists.txt b/src/coeffs/CMakeLists.txt new file mode 100644 index 0000000..d0b4672 --- /dev/null +++ b/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_link_libraries(nucfd_coeffs PRIVATE nucfd_compiler_flags) diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt new file mode 100644 index 0000000..61dbdaf --- /dev/null +++ b/tests/CMakeLists.txt @@ -0,0 +1,37 @@ +## tests/CMakeLists.txt +# +## Description +# +# CMakeLists file controlling the tests for the NuCFD project. +# +## LICENSE +# +# SPDX-License-Identifier: BSD-3-Clause +# + +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 nucfd_tests_mod.f90) +target_link_libraries(nucfd_test_framework nucfd_test_flags) + +function(define_test suite test_name) + add_executable(${test_name} ${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_subdirectory(tridsolver) + +# Finite difference test suite +add_subdirectory(fd-schemes) diff --git a/tests/fd-schemes/CMakeLists.txt b/tests/fd-schemes/CMakeLists.txt new file mode 100644 index 0000000..f32d3dc --- /dev/null +++ b/tests/fd-schemes/CMakeLists.txt @@ -0,0 +1,13 @@ +## tests/fd-schemes/CMakeLists.txt +# +## Description +# +# CMakeLists file for the finite difference schemes test suite. +# +## LICENSE +# +# SPDX-License-Identifier: BSD-3-Clause +# + +define_test(fd-schemes verify_coeffs) +define_test(fd-schemes differentiation_rules) diff --git a/tests/tridsolver/CMakeLists.txt b/tests/tridsolver/CMakeLists.txt new file mode 100644 index 0000000..56cb4b5 --- /dev/null +++ b/tests/tridsolver/CMakeLists.txt @@ -0,0 +1,23 @@ +## tests/tridsolver/CMakeLists.txt +# +## Description +# +# CMakeLists file for the tri-diagonal solver test suite. +# +## LICENSE +# +# SPDX-License-Identifier: BSD-3-Clause +# + +add_library(trid_test_utils 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) From 7001531f8a5085f4dfd4d752f5604cd23b4b66de Mon Sep 17 00:00:00 2001 From: Paul Date: Sun, 5 Mar 2023 12:45:53 +0000 Subject: [PATCH 53/58] Add option to build shared libraries OFF by default --- CMakeLists.txt | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/CMakeLists.txt b/CMakeLists.txt index 82a2cde..783a8d9 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -14,6 +14,11 @@ 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}) From ae103985e28c9e78edfa190365ea2523f0b04596 Mon Sep 17 00:00:00 2001 From: Paul Bartholomew Date: Sun, 5 Mar 2023 20:39:57 +0000 Subject: [PATCH 54/58] Split FD schemes test suite for separate function types --- tests/fd-schemes/CMakeLists.txt | 3 +- .../differentiation_rules/CMakeLists.txt | 14 ++++ .../constant_function.f90 | 77 +++++++++++++++++ .../linear_falling_function.f90 | 82 +++++++++++++++++++ .../linear_rising_function.f90} | 39 ++------- 5 files changed, 180 insertions(+), 35 deletions(-) create mode 100644 tests/fd-schemes/differentiation_rules/CMakeLists.txt create mode 100644 tests/fd-schemes/differentiation_rules/constant_function.f90 create mode 100644 tests/fd-schemes/differentiation_rules/linear_falling_function.f90 rename tests/fd-schemes/{differentiation_rules.f90 => differentiation_rules/linear_rising_function.f90} (51%) diff --git a/tests/fd-schemes/CMakeLists.txt b/tests/fd-schemes/CMakeLists.txt index f32d3dc..ed204b6 100644 --- a/tests/fd-schemes/CMakeLists.txt +++ b/tests/fd-schemes/CMakeLists.txt @@ -10,4 +10,5 @@ # define_test(fd-schemes verify_coeffs) -define_test(fd-schemes differentiation_rules) + +add_subdirectory(differentiation_rules) diff --git a/tests/fd-schemes/differentiation_rules/CMakeLists.txt b/tests/fd-schemes/differentiation_rules/CMakeLists.txt new file mode 100644 index 0000000..3ffc88a --- /dev/null +++ b/tests/fd-schemes/differentiation_rules/CMakeLists.txt @@ -0,0 +1,14 @@ +## tests/fd-schemes/differentiation_rules/CMakeLists.txt +# +## Description +# +# CMakeLists file for the finite difference schemes differentiation rules test sub-suite. +# +## LICENSE +# +# SPDX-License-Identifier: BSD-3-Clause +# + +define_test(fd-schemes constant_function) +define_test(fd-schemes linear_rising_function) +define_test(fd-schemes linear_falling_function) diff --git a/tests/fd-schemes/differentiation_rules/constant_function.f90 b/tests/fd-schemes/differentiation_rules/constant_function.f90 new file mode 100644 index 0000000..06ea223 --- /dev/null +++ b/tests/fd-schemes/differentiation_rules/constant_function.f90 @@ -0,0 +1,77 @@ +! tests/fd-schemes/differentiation_rules/constant_function.f90 +! +!! Part of the fd-schemes test suite. +! +! SPDX-License-Identifier: BSD-3-Clause + +program constant_function + !! Tests that rules of differentiation are respected for constant functions. + + use nucfd_types + use nucfd_deriv + + use nucfd_tests + + implicit none + + integer :: n ! Grid size + real :: L ! Domain size + real :: h ! Grid spacing + + real, dimension(:), allocatable :: x ! Grid + real, dimension(:), allocatable :: f ! Function + real :: dfdx ! Derivative + real :: dgdx ! Derivative + + ! Stencil + type(nucfd_index_stencil) :: stencil + integer :: i + integer, parameter :: width = 5 + integer, parameter :: centre = 3 + + call initialise_suite("Constant function differentiation rules") + + n = 33 + L = 1.0 + h = L / real(n - 1) + + allocate(x(n)) + allocate(f(n)) + do i = 1, n + x(i) = real(i - 1) * h + end do + + print *, "+++ Initialising stencil +++" + call create_stencil(width, centre, stencil) + i = 33 / 2 + select type(indices => stencil%stencil) + type is(integer) + indices(-2) = i - 2 + indices(-1) = i - 1 + indices(+0) = i + 0 + indices(+1) = i + 1 + indices(+2) = i + 2 + class default + print *, "Error: Index stencil is misallocated!" + error stop + end select + + f(:) = 1.0 + call deriv_rhs(f, stencil, x, dfdx) + call deriv_rhs(f + 1.0, stencil, x, dgdx) + call test_report("Shifted derivative +, const f", check_scalar(dgdx, dfdx)) + call deriv_rhs(f - 1.0, stencil, x, dgdx) + call test_report("Shifted derivative -, const f", check_scalar(dgdx, dfdx)) + call deriv_rhs(f, stencil, x, dgdx) + call test_report("Shifted derivative 0, const f", check_scalar(dgdx, dfdx)) + call deriv_rhs(2.0 * f, stencil, x, dgdx) + call test_report("Scaled derivative 2x, const f", check_scalar(dgdx, 2.0 * dfdx)) + call deriv_rhs(-f, stencil, x, dgdx) + call test_report("Scaled derivative -1, const f", check_scalar(dgdx, -dfdx)) + + deallocate(x) + deallocate(f) + + call finalise_suite() + +end program constant_function diff --git a/tests/fd-schemes/differentiation_rules/linear_falling_function.f90 b/tests/fd-schemes/differentiation_rules/linear_falling_function.f90 new file mode 100644 index 0000000..eab0827 --- /dev/null +++ b/tests/fd-schemes/differentiation_rules/linear_falling_function.f90 @@ -0,0 +1,82 @@ +! tests/fd-schemes/differentiation_rules/linear_falling_function.f90 +! +!! Part of the fd-schemes test suite. +! +! SPDX-License-Identifier: BSD-3-Clause + +program linear_falling_function + !! Tests that rules of differentiation are respected for linearly decreasing functions. + + use nucfd_types + use nucfd_deriv + + use nucfd_tests + + implicit none + + integer :: n ! Grid size + real :: L ! Domain size + real :: h ! Grid spacing + + real, dimension(:), allocatable :: x ! Grid + real, dimension(:), allocatable :: f ! Function + real :: dfdx ! Derivative + real :: dgdx ! Derivative + + ! Stencil + type(nucfd_index_stencil) :: stencil + integer :: i + integer, parameter :: width = 5 + integer, parameter :: centre = 3 + + call initialise_suite("Linearly decreasing function differentiation rules") + + n = 33 + L = 1.0 + h = L / real(n - 1) + + allocate(x(n)) + allocate(f(n)) + do i = 1, n + x(i) = real(i - 1) * h + end do + + print *, "+++ Initialising stencil +++" + call create_stencil(width, centre, stencil) + i = 33 / 2 + select type(indices => stencil%stencil) + type is(integer) + indices(-2) = i - 2 + indices(-1) = i - 1 + indices(+0) = i + 0 + indices(+1) = i + 1 + indices(+2) = i + 2 + class default + print *, "Error: Index stencil is misallocated!" + error stop + end select + + f(1) = 0.0 + dfdx = 1.0 + do i = 2, n + f(i) = f(i - 1) - h * dfdx + end do + + call deriv_rhs(f, stencil, x, dfdx) + call deriv_rhs(f + 1.0, stencil, x, dgdx) + call test_report("Shifted derivative +, linear- f", check_scalar(dgdx, dfdx)) + call deriv_rhs(f - 1.0, stencil, x, dgdx) + call test_report("Shifted derivative -, linear- f", check_scalar(dgdx, dfdx)) + call deriv_rhs(f, stencil, x, dgdx) + call test_report("Shifted derivative 0, linear- f", check_scalar(dgdx, dfdx)) + call deriv_rhs(2.0 * f, stencil, x, dgdx) + call test_report("Scaled derivative 2x, linear- f", check_scalar(dgdx, 2.0 * dfdx)) + call deriv_rhs(-f, stencil, x, dgdx) + call test_report("Scaled derivative -1, linear- f", check_scalar(dgdx, -dfdx)) + + deallocate(x) + deallocate(f) + + call finalise_suite() + +end program linear_falling_function diff --git a/tests/fd-schemes/differentiation_rules.f90 b/tests/fd-schemes/differentiation_rules/linear_rising_function.f90 similarity index 51% rename from tests/fd-schemes/differentiation_rules.f90 rename to tests/fd-schemes/differentiation_rules/linear_rising_function.f90 index f5d507f..4400537 100644 --- a/tests/fd-schemes/differentiation_rules.f90 +++ b/tests/fd-schemes/differentiation_rules/linear_rising_function.f90 @@ -1,11 +1,11 @@ -! tests/fd-schemes/differentiation_rules.f90 +! tests/fd-schemes/differentiation_rules/linear_rising_function.f90 ! !! Part of the fd-schemes test suite. ! ! SPDX-License-Identifier: BSD-3-Clause -program differentiation_rules - !! Tests that rules of differentiation are respected. +program linear_rising_function + !! Tests that rules of differentiation are respected for linearly increasing functions. use nucfd_types use nucfd_deriv @@ -29,7 +29,7 @@ program differentiation_rules integer, parameter :: width = 5 integer, parameter :: centre = 3 - call initialise_suite("Differentiation rules") + call initialise_suite("Linearly increasing function differentiation rules") n = 33 L = 1.0 @@ -56,21 +56,6 @@ program differentiation_rules error stop end select - print *, "+++ Testing derivative of constant function +++" - f(:) = 1.0 - call deriv_rhs(f, stencil, x, dfdx) - call deriv_rhs(f + 1.0, stencil, x, dgdx) - call test_report("Shifted derivative +, const f", check_scalar(dgdx, dfdx)) - call deriv_rhs(f - 1.0, stencil, x, dgdx) - call test_report("Shifted derivative -, const f", check_scalar(dgdx, dfdx)) - call deriv_rhs(f, stencil, x, dgdx) - call test_report("Shifted derivative 0, const f", check_scalar(dgdx, dfdx)) - call deriv_rhs(2.0 * f, stencil, x, dgdx) - call test_report("Scaled derivative 2x, const f", check_scalar(dgdx, 2.0 * dfdx)) - call deriv_rhs(-f, stencil, x, dgdx) - call test_report("Scaled derivative -1, const f", check_scalar(dgdx, -dfdx)) - - print *, "+++ Testing derivative of linearly increasing function +++" f(1) = 0.0 dfdx = 1.0 do i = 2, n @@ -87,24 +72,10 @@ program differentiation_rules call test_report("Scaled derivative 2x, linear+ f", check_scalar(dgdx, 2.0 * dfdx)) call deriv_rhs(-f, stencil, x, dgdx) call test_report("Scaled derivative -1, linear+ f", check_scalar(dgdx, -dfdx)) - - print *, "+++ Testing derivative of linearly decreasing function +++" - f(:) = -f(:) - call deriv_rhs(f, stencil, x, dfdx) - call deriv_rhs(f + 1.0, stencil, x, dgdx) - call test_report("Shifted derivative +, linear- f", check_scalar(dgdx, dfdx)) - call deriv_rhs(f - 1.0, stencil, x, dgdx) - call test_report("Shifted derivative -, linear- f", check_scalar(dgdx, dfdx)) - call deriv_rhs(f, stencil, x, dgdx) - call test_report("Shifted derivative 0, linear- f", check_scalar(dgdx, dfdx)) - call deriv_rhs(2.0 * f, stencil, x, dgdx) - call test_report("Scaled derivative 2x, linear- f", check_scalar(dgdx, 2.0 * dfdx)) - call deriv_rhs(-f, stencil, x, dgdx) - call test_report("Scaled derivative -1, linear- f", check_scalar(dgdx, -dfdx)) deallocate(x) deallocate(f) call finalise_suite() -end program differentiation_rules +end program linear_rising_function From 25fd40361dd5fe33b2303a8e7b3aa11e35e00997 Mon Sep 17 00:00:00 2001 From: Paul Bartholomew Date: Sun, 5 Mar 2023 20:57:07 +0000 Subject: [PATCH 55/58] Split FD schemes coefficient test suite --- tests/fd-schemes/CMakeLists.txt | 3 +- tests/fd-schemes/verify_coeffs.f90 | 108 ------------------ tests/fd-schemes/verify_coeffs/CMakeLists.txt | 16 +++ .../verify_coeffs/verify_coeff_a.f90 | 61 ++++++++++ .../verify_coeffs/verify_coeff_b.f90 | 61 ++++++++++ .../verify_coeffs/verify_coeff_c.f90 | 59 ++++++++++ .../verify_coeffs/verify_coeff_d.f90 | 59 ++++++++++ .../verify_coeffs/verify_coeff_e.f90 | 50 ++++++++ 8 files changed, 307 insertions(+), 110 deletions(-) delete mode 100644 tests/fd-schemes/verify_coeffs.f90 create mode 100644 tests/fd-schemes/verify_coeffs/CMakeLists.txt create mode 100644 tests/fd-schemes/verify_coeffs/verify_coeff_a.f90 create mode 100644 tests/fd-schemes/verify_coeffs/verify_coeff_b.f90 create mode 100644 tests/fd-schemes/verify_coeffs/verify_coeff_c.f90 create mode 100644 tests/fd-schemes/verify_coeffs/verify_coeff_d.f90 create mode 100644 tests/fd-schemes/verify_coeffs/verify_coeff_e.f90 diff --git a/tests/fd-schemes/CMakeLists.txt b/tests/fd-schemes/CMakeLists.txt index ed204b6..27bea18 100644 --- a/tests/fd-schemes/CMakeLists.txt +++ b/tests/fd-schemes/CMakeLists.txt @@ -9,6 +9,5 @@ # SPDX-License-Identifier: BSD-3-Clause # -define_test(fd-schemes verify_coeffs) - +add_subdirectory(verify_coeffs) add_subdirectory(differentiation_rules) diff --git a/tests/fd-schemes/verify_coeffs.f90 b/tests/fd-schemes/verify_coeffs.f90 deleted file mode 100644 index e27c298..0000000 --- a/tests/fd-schemes/verify_coeffs.f90 +++ /dev/null @@ -1,108 +0,0 @@ -! tests/fd-schemes/verify_coeffs.f90 -! -!! Part of the fd-schemes test suite. -! -! SPDX-License-Identifier: BSD-3-Clause - -program verify_coeffs - !! Tests the computation of finite difference coefficients for non-uniform grids. - - use nucfd_types - use nucfd_coeffs - - use nucfd_tests - - implicit none - - real :: n ! Mesh size - real :: L ! Domain size - real :: h ! Average grid spacing. - - type(nucfd_stencil_points) :: stencil ! Stencil of grid points. - - real :: a, b, c, d, e - real, parameter :: aref = 14.0 / 9.0 - real, parameter :: bref = -aref - real, parameter :: cref = 1.0 / 9.0 - real, parameter :: dref = -cref - real, parameter :: eref = 0.0 - - real :: numerator, numerator_corr, denominator, divisor - real, parameter :: numerator_f1ref = 14.0 / 3.0 - real, parameter :: numerator_corr_f1ref = 0.0 - real, parameter :: denominator_f1ref = 3.0 - real, parameter :: divisor_f1ref = 2.0 - real, parameter :: numerator_f2ref = 2.0 / 3.0 - real, parameter :: denominator_f2ref = 6.0 - real, parameter :: divisor_f2ref = 4.0 - - call initialise_suite("Verify coefficients") - - n = 128 - L = 1.0 - h = L / real(n - 1) - - call create_stencil(5, 3, stencil) - select type(points => stencil%stencil) - type is(real) - points(:) = 0.0 - points(-2) = -2.0 * h - points(-1) = -1.0 * h - points(0) = 0.0 - points(1) = +1.0 * h - points(2) = +2.0 * h - class default - print *, "Error: Coordinate stencil is misallocated!" - end select - - call coeff_a_components(points_to_deltas(stencil), numerator, numerator_corr, denominator, divisor) - call test_report("Coefficient A numerator", check_scalar(numerator, numerator_f1ref * (h**3))) - call test_report("Coefficient A numerator correction", check_scalar(numerator_corr, numerator_corr_f1ref * (h**3))) - call test_report("Coefficient A denominator", check_scalar(denominator, denominator_f1ref * (h**3))) - call test_report("Coefficient A divisor", check_scalar(divisor, divisor_f1ref * h)) - a = coeff_a(stencil) - call test_report("Coefficient A", check_scalar(a, aref / (2.0 * h))) - - call coeff_b_components(points_to_deltas(stencil), numerator, numerator_corr, denominator, divisor) - call test_report("Coefficient B numerator", check_scalar(numerator, -numerator_f1ref * (h**3))) - call test_report("Coefficient B numerator correction", check_scalar(numerator_corr, -numerator_corr_f1ref * (h**3))) - call test_report("Coefficient B denominator", check_scalar(denominator, denominator_f1ref * (h**3))) - call test_report("Coefficient B divisor", check_scalar(divisor, divisor_f1ref * h)) - b = coeff_b(stencil) - call test_report("Coefficient B", check_scalar(b, bref / (2.0 * h))) - - call coeff_c_components(points_to_deltas(stencil), numerator, denominator, divisor) - call test_report("Coefficient C numerator", check_scalar(numerator, numerator_f2ref * (h**3))) - call test_report("Coefficient C denominator", check_scalar(denominator, denominator_f2ref * (h**3))) - call test_report("Coefficient C divisor", check_scalar(divisor, divisor_f2ref * h)) - c = coeff_c(stencil) - call test_report("Coefficient C", check_scalar(c, cref / (4.0 * h))) - - call coeff_d_components(points_to_deltas(stencil), numerator, denominator, divisor) - call test_report("Coefficient C numerator", check_scalar(numerator, -numerator_f2ref * (h**3))) - call test_report("Coefficient C denominator", check_scalar(denominator, denominator_f2ref * (h**3))) - call test_report("Coefficient C divisor", check_scalar(divisor, divisor_f2ref * h)) - d = coeff_d(stencil) - call test_report("Coefficient D", check_scalar(d, dref / (4.0 * h))) - - e = coeff_e(stencil) - call test_report("Coefficient E", check_scalar(e, eref)) - - select type(points => stencil%stencil) - type is(real) - points(-2) = -4.0 * h - points(2) = +3.0 * h - class default - print *, "Error: Coordinate stencil is misallocated!" - end select - - a = coeff_a(stencil) - b = coeff_b(stencil) - call test_report("Coefficient A /= B", .not. check_scalar(a, b)) - c = coeff_c(stencil) - d = coeff_d(stencil) - call test_report("Coefficient C /= D", .not. check_scalar(c, d)) - - call finalise_suite() - -end program verify_coeffs diff --git a/tests/fd-schemes/verify_coeffs/CMakeLists.txt b/tests/fd-schemes/verify_coeffs/CMakeLists.txt new file mode 100644 index 0000000..f95f68c --- /dev/null +++ b/tests/fd-schemes/verify_coeffs/CMakeLists.txt @@ -0,0 +1,16 @@ +## tests/fd-schemes/verify_coeffs/CMakeLists.txt +# +## Description +# +# CMakeLists file for the finite difference schemes differentiation rules test sub-suite. +# +## LICENSE +# +# SPDX-License-Identifier: BSD-3-Clause +# + +define_test(fd-schemes verify_coeff_a) +define_test(fd-schemes verify_coeff_b) +define_test(fd-schemes verify_coeff_c) +define_test(fd-schemes verify_coeff_d) +define_test(fd-schemes verify_coeff_e) diff --git a/tests/fd-schemes/verify_coeffs/verify_coeff_a.f90 b/tests/fd-schemes/verify_coeffs/verify_coeff_a.f90 new file mode 100644 index 0000000..9c3e5b4 --- /dev/null +++ b/tests/fd-schemes/verify_coeffs/verify_coeff_a.f90 @@ -0,0 +1,61 @@ +! tests/fd-schemes/verify_coeffs/verify_coeff_a.f90 +! +!! Part of the fd-schemes test suite. +! +! SPDX-License-Identifier: BSD-3-Clause + +program verify_coeff_b + !! Tests the computation of finite difference coefficient for non-uniform grids acting on f_{i+1}. + + use nucfd_types + use nucfd_coeffs + + use nucfd_tests + + implicit none + + real :: n ! Mesh size + real :: L ! Domain size + real :: h ! Average grid spacing. + + type(nucfd_stencil_points) :: stencil ! Stencil of grid points. + + real :: a + real, parameter :: aref = 14.0 / 9.0 + + real :: numerator, numerator_corr, denominator, divisor + real, parameter :: numerator_f1ref = 14.0 / 3.0 + real, parameter :: numerator_corr_f1ref = 0.0 + real, parameter :: denominator_f1ref = 3.0 + real, parameter :: divisor_f1ref = 2.0 + + call initialise_suite("Verify coefficient A") + + n = 128 + L = 1.0 + h = L / real(n - 1) + + call create_stencil(5, 3, stencil) + select type(points => stencil%stencil) + type is(real) + points(:) = 0.0 + points(-2) = -2.0 * h + points(-1) = -1.0 * h + points(0) = 0.0 + points(1) = +1.0 * h + points(2) = +2.0 * h + class default + print *, "Error: Coordinate stencil is misallocated!" + end select + + call coeff_a_components(points_to_deltas(stencil), numerator, numerator_corr, denominator, divisor) + call test_report("Coefficient A numerator", check_scalar(numerator, numerator_f1ref * (h**3))) + call test_report("Coefficient A numerator correction", check_scalar(numerator_corr, numerator_corr_f1ref * (h**3))) + call test_report("Coefficient A denominator", check_scalar(denominator, denominator_f1ref * (h**3))) + call test_report("Coefficient A divisor", check_scalar(divisor, divisor_f1ref * h)) + a = coeff_a(stencil) + call test_report("Coefficient A", check_scalar(a, aref / (2.0 * h))) + + call finalise_suite() + +end program verify_coeff_b diff --git a/tests/fd-schemes/verify_coeffs/verify_coeff_b.f90 b/tests/fd-schemes/verify_coeffs/verify_coeff_b.f90 new file mode 100644 index 0000000..dcd013f --- /dev/null +++ b/tests/fd-schemes/verify_coeffs/verify_coeff_b.f90 @@ -0,0 +1,61 @@ +! tests/fd-schemes/verify_coeffs/verify_coeff_b.f90 +! +!! Part of the fd-schemes test suite. +! +! SPDX-License-Identifier: BSD-3-Clause + +program verify_coeff_b + !! Tests the computation of finite difference coefficient for non-uniform grids acting on f_{i-1}. + + use nucfd_types + use nucfd_coeffs + + use nucfd_tests + + implicit none + + real :: n ! Mesh size + real :: L ! Domain size + real :: h ! Average grid spacing. + + type(nucfd_stencil_points) :: stencil ! Stencil of grid points. + + real :: b + real, parameter :: bref = -14.0 / 9.0 + + real :: numerator, numerator_corr, denominator, divisor + real, parameter :: numerator_f1ref = -14.0 / 3.0 + real, parameter :: numerator_corr_f1ref = -0.0 + real, parameter :: denominator_f1ref = 3.0 + real, parameter :: divisor_f1ref = 2.0 + + call initialise_suite("Verify coefficient B") + + n = 128 + L = 1.0 + h = L / real(n - 1) + + call create_stencil(5, 3, stencil) + select type(points => stencil%stencil) + type is(real) + points(:) = 0.0 + points(-2) = -2.0 * h + points(-1) = -1.0 * h + points(0) = 0.0 + points(1) = +1.0 * h + points(2) = +2.0 * h + class default + print *, "Error: Coordinate stencil is misallocated!" + end select + + call coeff_b_components(points_to_deltas(stencil), numerator, numerator_corr, denominator, divisor) + call test_report("Coefficient B numerator", check_scalar(numerator, numerator_f1ref * (h**3))) + call test_report("Coefficient B numerator correction", check_scalar(numerator_corr, numerator_corr_f1ref * (h**3))) + call test_report("Coefficient B denominator", check_scalar(denominator, denominator_f1ref * (h**3))) + call test_report("Coefficient B divisor", check_scalar(divisor, divisor_f1ref * h)) + b = coeff_b(stencil) + call test_report("Coefficient B", check_scalar(b, bref / (2.0 * h))) + + call finalise_suite() + +end program verify_coeff_b diff --git a/tests/fd-schemes/verify_coeffs/verify_coeff_c.f90 b/tests/fd-schemes/verify_coeffs/verify_coeff_c.f90 new file mode 100644 index 0000000..16a8787 --- /dev/null +++ b/tests/fd-schemes/verify_coeffs/verify_coeff_c.f90 @@ -0,0 +1,59 @@ +! tests/fd-schemes/verify_coeffs/verify_coeff_c.f90 +! +!! Part of the fd-schemes test suite. +! +! SPDX-License-Identifier: BSD-3-Clause + +program verify_coeff_c + !! Tests the computation of finite difference coefficient for non-uniform grids acting on f_{i+2}. + + use nucfd_types + use nucfd_coeffs + + use nucfd_tests + + implicit none + + real :: n ! Mesh size + real :: L ! Domain size + real :: h ! Average grid spacing. + + type(nucfd_stencil_points) :: stencil ! Stencil of grid points. + + real :: c + real, parameter :: cref = 1.0 / 9.0 + + real :: numerator, denominator, divisor + real, parameter :: numerator_f2ref = 2.0 / 3.0 + real, parameter :: denominator_f2ref = 6.0 + real, parameter :: divisor_f2ref = 4.0 + + call initialise_suite("Verify coefficient C") + + n = 128 + L = 1.0 + h = L / real(n - 1) + + call create_stencil(5, 3, stencil) + select type(points => stencil%stencil) + type is(real) + points(:) = 0.0 + points(-2) = -2.0 * h + points(-1) = -1.0 * h + points(0) = 0.0 + points(1) = +1.0 * h + points(2) = +2.0 * h + class default + print *, "Error: Coordinate stencil is misallocated!" + end select + + call coeff_c_components(points_to_deltas(stencil), numerator, denominator, divisor) + call test_report("Coefficient C numerator", check_scalar(numerator, numerator_f2ref * (h**3))) + call test_report("Coefficient C denominator", check_scalar(denominator, denominator_f2ref * (h**3))) + call test_report("Coefficient C divisor", check_scalar(divisor, divisor_f2ref * h)) + c = coeff_c(stencil) + call test_report("Coefficient C", check_scalar(c, cref / (4.0 * h))) + + call finalise_suite() + +end program verify_coeff_c diff --git a/tests/fd-schemes/verify_coeffs/verify_coeff_d.f90 b/tests/fd-schemes/verify_coeffs/verify_coeff_d.f90 new file mode 100644 index 0000000..fcc22ca --- /dev/null +++ b/tests/fd-schemes/verify_coeffs/verify_coeff_d.f90 @@ -0,0 +1,59 @@ +! tests/fd-schemes/verify_coeffs/verify_coeff_d.f90 +! +!! Part of the fd-schemes test suite. +! +! SPDX-License-Identifier: BSD-3-Clause + +program verify_coeff_d + !! Tests the computation of finite difference coefficient for non-uniform grids acting on f_{i-2}. + + use nucfd_types + use nucfd_coeffs + + use nucfd_tests + + implicit none + + real :: n ! Mesh size + real :: L ! Domain size + real :: h ! Average grid spacing. + + type(nucfd_stencil_points) :: stencil ! Stencil of grid points. + + real :: d + real, parameter :: dref = -1.0 / 9.0 + + real :: numerator, denominator, divisor + real, parameter :: numerator_f2ref =- 2.0 / 3.0 + real, parameter :: denominator_f2ref = 6.0 + real, parameter :: divisor_f2ref = 4.0 + + call initialise_suite("Verify coefficient D") + + n = 128 + L = 1.0 + h = L / real(n - 1) + + call create_stencil(5, 3, stencil) + select type(points => stencil%stencil) + type is(real) + points(:) = 0.0 + points(-2) = -2.0 * h + points(-1) = -1.0 * h + points(0) = 0.0 + points(1) = +1.0 * h + points(2) = +2.0 * h + class default + print *, "Error: Coordinate stencil is misallocated!" + end select + + call coeff_d_components(points_to_deltas(stencil), numerator, denominator, divisor) + call test_report("Coefficient D numerator", check_scalar(numerator, numerator_f2ref * (h**3))) + call test_report("Coefficient D denominator", check_scalar(denominator, denominator_f2ref * (h**3))) + call test_report("Coefficient D divisor", check_scalar(divisor, divisor_f2ref * h)) + d = coeff_d(stencil) + call test_report("Coefficient D", check_scalar(d, dref / (4.0 * h))) + + call finalise_suite() + +end program verify_coeff_d diff --git a/tests/fd-schemes/verify_coeffs/verify_coeff_e.f90 b/tests/fd-schemes/verify_coeffs/verify_coeff_e.f90 new file mode 100644 index 0000000..9f5f75a --- /dev/null +++ b/tests/fd-schemes/verify_coeffs/verify_coeff_e.f90 @@ -0,0 +1,50 @@ +! tests/fd-schemes/verify_coeffs/verify_coeff_e.f90 +! +!! Part of the fd-schemes test suite. +! +! SPDX-License-Identifier: BSD-3-Clause + +program verify_coeff_e + !! Tests the computation of finite difference coefficient for non-uniform grids acting on f_{i}. + + use nucfd_types + use nucfd_coeffs + + use nucfd_tests + + implicit none + + real :: n ! Mesh size + real :: L ! Domain size + real :: h ! Average grid spacing. + + type(nucfd_stencil_points) :: stencil ! Stencil of grid points. + + real :: e + real, parameter :: eref = 0.0 + + call initialise_suite("Verify coefficient E") + + n = 128 + L = 1.0 + h = L / real(n - 1) + + call create_stencil(5, 3, stencil) + select type(points => stencil%stencil) + type is(real) + points(:) = 0.0 + points(-2) = -2.0 * h + points(-1) = -1.0 * h + points(0) = 0.0 + points(1) = +1.0 * h + points(2) = +2.0 * h + class default + print *, "Error: Coordinate stencil is misallocated!" + end select + + e = coeff_e(stencil) + call test_report("Coefficient E", check_scalar(e, eref / (4.0 * h))) + + call finalise_suite() + +end program verify_coeff_e From 4a217accba42e1ee37475f34a329528a4f83abfc Mon Sep 17 00:00:00 2001 From: Paul Bartholomew Date: Sun, 5 Mar 2023 21:09:09 +0000 Subject: [PATCH 56/58] Fix error in setting Debug flags --- CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 783a8d9..2f62c41 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -43,7 +43,7 @@ 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 nucfd_debug_flags) + target_link_libraries(nucfd_compiler_flags INTERFACE nucfd_debug_flags) endif() ## Build nucfd library From 1c01908a65cc1852618840a8ae6cbf0a8ed1b5d1 Mon Sep 17 00:00:00 2001 From: Paul Bartholomew Date: Sat, 18 Mar 2023 23:18:01 +0000 Subject: [PATCH 57/58] Fix typo --- tests/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 61dbdaf..d2342c1 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -30,7 +30,7 @@ function(add_test_libs suite test_name lib) target_link_libraries(${test_name} ${lib}) endfunction() -# Tidiagonal solver test suite +# Tridiagonal solver test suite add_subdirectory(tridsolver) # Finite difference test suite From 9352fc5fccf6cffc4c9c8e092f4647e73e98ec22 Mon Sep 17 00:00:00 2001 From: Paul Bartholomew Date: Mon, 20 Mar 2023 12:02:08 +0000 Subject: [PATCH 58/58] Update CHANGELOG --- CHANGELOG.md | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 2eff767..c180acf 100644 --- a/CHANGELOG.md +++ b/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