From 24a9f6ac87c550a6f0fc16f9d1d4a7a6e91e926b Mon Sep 17 00:00:00 2001
From: Paul <ptb0890@gmail.com>
Date: Sun, 26 Feb 2023 11:23:45 +0000
Subject: [PATCH 01/22] Being implementing example program to compute
 derivatives

---
 CMakeLists.txt                   |  4 ++++
 examples/compute-derivatives.f90 | 17 +++++++++++++++++
 2 files changed, 21 insertions(+)
 create mode 100644 examples/compute-derivatives.f90

diff --git a/CMakeLists.txt b/CMakeLists.txt
index a11e24e..0dc60cf 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -78,6 +78,10 @@ add_test_libs(tridsolver system_00_anti-symm trid_test_utils)
 # Finite difference test suite
 define_test(fd-schemes verify_coeffs)
 
+## Examples
+add_executable(compute_derivatives examples/compute-derivatives.f90)
+target_link_libraries(compute_derivatives nucfd)
+
 ## Documentation
 add_custom_target(doc ford
   ${CMAKE_SOURCE_DIR}/NuCFD.md
diff --git a/examples/compute-derivatives.f90 b/examples/compute-derivatives.f90
new file mode 100644
index 0000000..6574cd1
--- /dev/null
+++ b/examples/compute-derivatives.f90
@@ -0,0 +1,17 @@
+!!!! examples/compute-derivatives.f90
+!!!
+!!!! Description
+!!!
+!!! A program that uses non-uniform compact finite difference schemes to compute derivatives on
+!!! (non-)uniform meshes.
+!!!
+!!!! LICENSE
+!!!
+!!! SPDX-License-Identifier: BSD-3-Clause
+!!
+
+program compute_derivatives
+
+  implicit none
+  
+end program compute_derivatives

From 5cfd84d584b57929d1d6dc1c76b7c74aa5b8b2d9 Mon Sep 17 00:00:00 2001
From: Paul <ptb0890@gmail.com>
Date: Sun, 26 Feb 2023 14:48:33 +0000
Subject: [PATCH 02/22] Update comments to work better with FORD

---
 NuCFD.md                         |  1 +
 examples/compute-derivatives.f90 | 16 ++++------------
 2 files changed, 5 insertions(+), 12 deletions(-)

diff --git a/NuCFD.md b/NuCFD.md
index 5af7c99..d992afc 100644
--- a/NuCFD.md
+++ b/NuCFD.md
@@ -4,6 +4,7 @@ author: Paul Bartholomew
 source: true
 src_dir: ./src
          ./tests
+		 ./examples
 ---
 
 {!README.md!}
diff --git a/examples/compute-derivatives.f90 b/examples/compute-derivatives.f90
index 6574cd1..4b7f768 100644
--- a/examples/compute-derivatives.f90
+++ b/examples/compute-derivatives.f90
@@ -1,16 +1,8 @@
-!!!! examples/compute-derivatives.f90
-!!!
-!!!! Description
-!!!
-!!! A program that uses non-uniform compact finite difference schemes to compute derivatives on
-!!! (non-)uniform meshes.
-!!!
-!!!! LICENSE
-!!!
-!!! SPDX-License-Identifier: BSD-3-Clause
-!!
-
 program compute_derivatives
+  !! A program that uses non-uniform compact finite difference schemes to compute derivatives on
+  !! (non-)uniform meshes.
+  !!
+  !! SPDX-License-Identifier: BSD-3-Clause
 
   implicit none
   

From e529de5ed391e96663261b8bee23bc07071c001f Mon Sep 17 00:00:00 2001
From: Paul <ptb0890@gmail.com>
Date: Sat, 4 Mar 2023 20:22:51 +0000
Subject: [PATCH 03/22] Add header to compute derivatives example

---
 examples/compute-derivatives.f90 | 8 ++++++--
 1 file changed, 6 insertions(+), 2 deletions(-)

diff --git a/examples/compute-derivatives.f90 b/examples/compute-derivatives.f90
index 4b7f768..e753b74 100644
--- a/examples/compute-derivatives.f90
+++ b/examples/compute-derivatives.f90
@@ -1,8 +1,12 @@
+! examples/compute-derivatives.f90
+!
+!! Example program to compute derivatives.
+!
+! SPDX-License-Identifier: BSD-3-Clause
+
 program compute_derivatives
   !! A program that uses non-uniform compact finite difference schemes to compute derivatives on
   !! (non-)uniform meshes.
-  !!
-  !! SPDX-License-Identifier: BSD-3-Clause
 
   implicit none
   

From d96377e295ffa044ed73cb34f5424adc37d2f523 Mon Sep 17 00:00:00 2001
From: Paul <ptb0890@gmail.com>
Date: Sat, 4 Mar 2023 23:57:58 +0000
Subject: [PATCH 04/22] Add a CMakeLists for examples

---
 CMakeLists.txt          |  3 +--
 examples/CMakeLists.txt | 13 +++++++++++++
 2 files changed, 14 insertions(+), 2 deletions(-)
 create mode 100644 examples/CMakeLists.txt

diff --git a/CMakeLists.txt b/CMakeLists.txt
index a13d768..749ee86 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -49,8 +49,7 @@ enable_testing()
 add_subdirectory(tests)
 
 ## Examples
-add_executable(compute_derivatives examples/compute-derivatives.f90)
-target_link_libraries(compute_derivatives nucfd)
+add_subdirectory(examples)
 
 ## Documentation
 add_custom_target(doc ford
diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
new file mode 100644
index 0000000..fdf0d43
--- /dev/null
+++ b/examples/CMakeLists.txt
@@ -0,0 +1,13 @@
+## examples/CMakeLists.txt
+#
+## Description
+#
+# CMakeLists file for building the NuCFD project examples.
+#
+## LICENSE
+#
+# SPDX-License-Identifier: BSD-3-Clause
+#
+
+add_executable(compute_derivatives compute-derivatives.f90)
+target_link_libraries(compute_derivatives nucfd)

From f94abc3743fd2c7711b07dc59bff37859d912224 Mon Sep 17 00:00:00 2001
From: Paul <ptb0890@gmail.com>
Date: Sun, 5 Mar 2023 11:16:44 +0000
Subject: [PATCH 05/22] Make examples build optional

---
 CMakeLists.txt          | 4 ++--
 examples/CMakeLists.txt | 5 +++++
 2 files changed, 7 insertions(+), 2 deletions(-)

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 749ee86..7d9a66f 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -48,8 +48,8 @@ add_subdirectory(src)
 enable_testing()
 add_subdirectory(tests)
 
-## Examples
-add_subdirectory(examples)
+## Examples (optional)
+add_subdirectory(examples EXCLUDE_FROM_ALL)
 
 ## Documentation
 add_custom_target(doc ford
diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
index fdf0d43..330c1f8 100644
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -9,5 +9,10 @@
 # SPDX-License-Identifier: BSD-3-Clause
 #
 
+project(examples)
+
 add_executable(compute_derivatives compute-derivatives.f90)
 target_link_libraries(compute_derivatives nucfd)
+
+add_custom_target(${PROJECT_NAME})
+add_dependencies(${PROJECT_NAME} compute_derivatives)

From 3f8e6428a59ed399ff6827049fb73e59b590f1b4 Mon Sep 17 00:00:00 2001
From: Paul <ptb0890@gmail.com>
Date: Sun, 5 Mar 2023 12:22:34 +0000
Subject: [PATCH 06/22] Make building tests optional

---
 CMakeLists.txt | 7 +++++--
 1 file changed, 5 insertions(+), 2 deletions(-)

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 7d9a66f..dc4d2f5 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -45,8 +45,11 @@ endif()
 add_subdirectory(src)
 
 ## Testing
-enable_testing()
-add_subdirectory(tests)
+# option(BUILD_TESTING "Build the testing tree" OFF) # Disables testing by default
+include(CTest)
+if(BUILD_TESTING)
+  add_subdirectory(tests)
+endif()
 
 ## Examples (optional)
 add_subdirectory(examples EXCLUDE_FROM_ALL)

From 7be116d52d6cdcb54fdee5fe534eb5ecc1ac058d Mon Sep 17 00:00:00 2001
From: Paul Bartholomew <p.bartholomew@epcc.ed.ac.uk>
Date: Sun, 5 Mar 2023 21:09:09 +0000
Subject: [PATCH 07/22] 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 1c9151a..d0d8dc8 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 ce7f9b06d650e1d25e3d65bac9ec9b8de1e47597 Mon Sep 17 00:00:00 2001
From: Paul Bartholomew <p.bartholomew@epcc.ed.ac.uk>
Date: Sun, 5 Mar 2023 21:09:33 +0000
Subject: [PATCH 08/22] Use object libraries for test utility libraries

---
 tests/CMakeLists.txt            |  5 +++--
 tests/tridsolver/CMakeLists.txt | 11 ++++++-----
 2 files changed, 9 insertions(+), 7 deletions(-)

diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 61dbdaf..fa8ce03 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -15,13 +15,14 @@ target_link_libraries(nucfd_test_flags INTERFACE
   nucfd_debug_flags
   nucfd_strict_flags)
 
-add_library(nucfd_test_framework nucfd_tests_mod.f90)
+add_library(nucfd_test_framework OBJECT
+  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} $<TARGET_OBJECTS:nucfd_test_framework>)
   target_link_libraries(${test_name} nucfd_test_flags)
   add_test(NAME ${suite}:${test_name} COMMAND ${test_name})
 endfunction()
diff --git a/tests/tridsolver/CMakeLists.txt b/tests/tridsolver/CMakeLists.txt
index 56cb4b5..2a9de12 100644
--- a/tests/tridsolver/CMakeLists.txt
+++ b/tests/tridsolver/CMakeLists.txt
@@ -9,15 +9,16 @@
 # SPDX-License-Identifier: BSD-3-Clause
 #
 
-add_library(trid_test_utils tridsol_test_utils_mod.f90)
+add_library(trid_test_utils OBJECT
+  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)
+add_test_libs(tridsolver system_11_symm $<TARGET_OBJECTS:trid_test_utils>)
+add_test_libs(tridsolver system_11_anti-symm $<TARGET_OBJECTS:trid_test_utils>)
 
 define_test(tridsolver system_00_symm nucfd_test_framework)
 define_test(tridsolver system_00_anti-symm)
-add_test_libs(tridsolver system_00_symm trid_test_utils)
-add_test_libs(tridsolver system_00_anti-symm trid_test_utils)
+add_test_libs(tridsolver system_00_symm $<TARGET_OBJECTS:trid_test_utils>)
+add_test_libs(tridsolver system_00_anti-symm $<TARGET_OBJECTS:trid_test_utils>)

From 10001869b84863faff544efe8553614dda33c719 Mon Sep 17 00:00:00 2001
From: Paul <ptb0890@gmail.com>
Date: Thu, 9 Mar 2023 21:10:15 +0000
Subject: [PATCH 09/22] Tidy up differentiation rules test suite

---
 .../differentiation_rules/CMakeLists.txt      |  5 +
 .../constant_function.f90                     | 32 ++-----
 .../diff_rules_utils.f90                      | 93 +++++++++++++++++++
 .../linear_falling_function.f90               | 34 ++-----
 .../linear_rising_function.f90                | 34 ++-----
 5 files changed, 120 insertions(+), 78 deletions(-)
 create mode 100644 tests/fd-schemes/differentiation_rules/diff_rules_utils.f90

diff --git a/tests/fd-schemes/differentiation_rules/CMakeLists.txt b/tests/fd-schemes/differentiation_rules/CMakeLists.txt
index 3ffc88a..e294d04 100644
--- a/tests/fd-schemes/differentiation_rules/CMakeLists.txt
+++ b/tests/fd-schemes/differentiation_rules/CMakeLists.txt
@@ -9,6 +9,11 @@
 # SPDX-License-Identifier: BSD-3-Clause
 #
 
+add_library(diff_rules_utils OBJECT diff_rules_utils.f90)
+
 define_test(fd-schemes constant_function)
+add_test_libs(fd-schemes constant_function diff_rules_utils)
 define_test(fd-schemes linear_rising_function)
+add_test_libs(fd-schemes linear_rising_function diff_rules_utils)
 define_test(fd-schemes linear_falling_function)
+add_test_libs(fd-schemes linear_falling_function diff_rules_utils)
diff --git a/tests/fd-schemes/differentiation_rules/constant_function.f90 b/tests/fd-schemes/differentiation_rules/constant_function.f90
index 06ea223..5b19e77 100644
--- a/tests/fd-schemes/differentiation_rules/constant_function.f90
+++ b/tests/fd-schemes/differentiation_rules/constant_function.f90
@@ -11,50 +11,30 @@ program constant_function
   use nucfd_deriv
   
   use nucfd_tests
+  use diff_rules_utils
   
   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
+  type(test_setup) :: ts
   
   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
+  ts = initialise_test(n, L)
+  f = allocate_test_array(ts)
+  x = create_mesh_array(ts)
+  stencil = build_stencil(ts)
 
   f(:) = 1.0
   call deriv_rhs(f, stencil, x, dfdx)
diff --git a/tests/fd-schemes/differentiation_rules/diff_rules_utils.f90 b/tests/fd-schemes/differentiation_rules/diff_rules_utils.f90
new file mode 100644
index 0000000..b4ce3bb
--- /dev/null
+++ b/tests/fd-schemes/differentiation_rules/diff_rules_utils.f90
@@ -0,0 +1,93 @@
+! tests/fd-schemes/differentiation_rules/diff_rules_utils.f90
+!
+!! Part of the fd-schemes test suite.
+!
+! SPDX-License-Identifier: BSD-3-Clause
+
+module diff_rules_utils
+  !! Utility module for differentiation rules test subsuite.
+
+  use nucfd_types
+  
+  implicit none
+
+  type test_setup
+     integer :: n
+     real :: L
+     real :: h
+  end type test_setup
+  
+  private
+  public :: test_setup
+  public :: initialise_test
+  public :: allocate_test_array
+  public :: create_mesh_array
+  public :: build_stencil
+  
+  integer, parameter :: width = 5
+  integer, parameter :: centre = 3
+  
+contains
+
+  type(test_setup) function initialise_test(n, L)
+
+    integer, intent(in) :: n
+    real, intent(in) :: L
+
+    real :: h
+
+    h = L / real(n - 1)
+    initialise_test = test_setup(n, L, h)
+
+  end function initialise_test
+
+  function allocate_test_array(ts) result(arr)
+
+    type(test_setup), intent(in) :: ts
+    real, dimension(:), allocatable :: arr
+
+    allocate(arr(ts%n))
+    
+  end function allocate_test_array
+
+  function create_mesh_array(ts) result(x)
+
+    type(test_setup), intent(in) :: ts
+    real, dimension(:), allocatable :: x
+
+    integer :: i
+    
+    x = allocate_test_array(ts)
+
+    associate(n => ts%n, h => ts%h)
+      do i = 1, n
+         x(i) = real(i - 1) * h
+      end do
+    end associate
+
+  end function create_mesh_array
+
+  type(nucfd_index_stencil) function build_stencil(ts)
+
+    type(test_setup), intent(in) :: ts
+
+    integer :: i
+    
+    print *, "+++ Initialising stencil +++"
+    call create_stencil(width, centre, build_stencil)
+    i = ts%n / 2
+    select type(indices => build_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
+
+  end function build_stencil
+  
+end module diff_rules_utils
diff --git a/tests/fd-schemes/differentiation_rules/linear_falling_function.f90 b/tests/fd-schemes/differentiation_rules/linear_falling_function.f90
index eab0827..fd16fc8 100644
--- a/tests/fd-schemes/differentiation_rules/linear_falling_function.f90
+++ b/tests/fd-schemes/differentiation_rules/linear_falling_function.f90
@@ -11,55 +11,37 @@ program linear_falling_function
   use nucfd_deriv
   
   use nucfd_tests
+  use diff_rules_utils
   
   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
+  type(test_setup) :: ts
+
   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
+  ts = initialise_test(n, L)
+  f = allocate_test_array(ts)
+  x = create_mesh_array(ts)
+  stencil = build_stencil(ts)
 
   f(1) = 0.0
   dfdx = 1.0
   do i = 2, n
-     f(i) = f(i - 1) - h * dfdx
+     f(i) = f(i - 1) - ts%h * dfdx
   end do
 
   call deriv_rhs(f, stencil, x, dfdx)
diff --git a/tests/fd-schemes/differentiation_rules/linear_rising_function.f90 b/tests/fd-schemes/differentiation_rules/linear_rising_function.f90
index 4400537..05d28b2 100644
--- a/tests/fd-schemes/differentiation_rules/linear_rising_function.f90
+++ b/tests/fd-schemes/differentiation_rules/linear_rising_function.f90
@@ -11,55 +11,37 @@ program linear_rising_function
   use nucfd_deriv
   
   use nucfd_tests
+  use diff_rules_utils
   
   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
+  type(test_setup) :: ts
+
   integer :: i
-  integer, parameter :: width = 5
-  integer, parameter :: centre = 3
   
   call initialise_suite("Linearly increasing 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
+  ts = initialise_test(n, L)
+  f = allocate_test_array(ts)
+  x = create_mesh_array(ts)
+  stencil = build_stencil(ts)
 
   f(1) = 0.0
   dfdx = 1.0
   do i = 2, n
-     f(i) = f(i - 1) + h * dfdx
+     f(i) = f(i - 1) + ts%h * dfdx
   end do
   call deriv_rhs(f, stencil, x, dfdx)
   call deriv_rhs(f + 1.0, stencil, x, dgdx)

From 2f65b6f8421186f411101af91305b2de34fcd377 Mon Sep 17 00:00:00 2001
From: Paul <ptb0890@gmail.com>
Date: Fri, 10 Mar 2023 11:08:53 +0000
Subject: [PATCH 10/22] Implement a grid generator for derivative example

---
 examples/compute-derivatives.f90 | 75 ++++++++++++++++++++++++++++++++
 1 file changed, 75 insertions(+)

diff --git a/examples/compute-derivatives.f90 b/examples/compute-derivatives.f90
index e753b74..2a2e8fa 100644
--- a/examples/compute-derivatives.f90
+++ b/examples/compute-derivatives.f90
@@ -9,5 +9,80 @@ program compute_derivatives
   !! (non-)uniform meshes.
 
   implicit none
+
+  real, dimension(:), allocatable :: x ! 1D "grid"
+  
+  call build_grid(11, 1.0, (/ 1, 2 /), x)
+  
+contains
+
+  subroutine build_grid(n, L, spacings, x)
+    !! Crude subroutine to build a 1-D grid.
+    !! The grid is described by the number of nodes, length and an array of relative grid spacings
+    !! for each block (constant within a block), returning a 1-D array of node centres.
+    !
+    !! @note The number of cells (n - 1) must be divisible by the number of blocks specified by the
+    !! spacings array otherwise an error will be raised.
+    
+    integer, intent(in) :: n                          !! The number of nodes in the mesh
+    real, intent(in) :: L                             !! The length of the domain
+    integer, dimension(:), intent(in) :: spacings     !! An array of the relative (integer) spacings
+                                                      !! for each block
+    real, dimension(:), allocatable, intent(out) :: x !! The grid node locations
+
+    integer :: ncells
+    integer :: nblocks
+    integer :: block_size
+    real :: block_length
+    real :: h
+
+    integer :: i
+    integer :: b
+    integer :: idx
+
+    allocate(x(n))
+
+    ncells = n - 1
+    nblocks = size(spacings, 1)
+    block_size = ncells / nblocks
+    if (block_size * nblocks + 1 /= n) then
+       print *, "Error: invalid grid specification!"
+       print *, " - Grid size n should be specified so that n - 1 is divisible by number of blocks"
+       print *, "   e.g. for 2 blocks a grid size of 11 is valid, 10 is not."
+       error stop
+    end if
+    block_length = L / real(sum(spacings))
+    h = block_length / real(block_size)
+
+    idx = 1
+    x(idx) = 0.0
+    idx = 2
+    do b = 1, nblocks
+       do i = 1, block_size
+          x(idx) = x(idx - 1) + real(spacings(b)) * h
+          idx = idx + 1
+       end do
+    end do
+
+    ! Self-check
+    if (any(x > L) .or. any(x < 0.0)) then
+       print *, "Error: grid exceeds specified length!"
+       print *, x
+       error stop
+    end if
+    if (abs(x(1)) > 0.0) then
+       print *, "Error: start point is wrong!"
+       print *, x(1)
+       error stop
+    end if
+    if (abs(x(n) - L) > (2 * epsilon(L) * L)) then
+       print *, "Error: end point is wrong!"
+       print *, "- Computed: ", x(n)
+       print *, "- Expected: ", L
+       print *, "- Relative error: ", abs(x(n) - L) / L
+       error stop
+    end if
+    
+  end subroutine build_grid
   
 end program compute_derivatives

From 0725ef992ea5a6b0bf37890120f93aef07055651 Mon Sep 17 00:00:00 2001
From: Paul <ptb0890@gmail.com>
Date: Fri, 10 Mar 2023 19:58:16 +0000
Subject: [PATCH 11/22] Derivative RHS calculation wasn't working with
 user-defined l/ubounds

---
 src/nucfd_deriv_mod.f90 | 16 ++++++++--------
 1 file changed, 8 insertions(+), 8 deletions(-)

diff --git a/src/nucfd_deriv_mod.f90 b/src/nucfd_deriv_mod.f90
index f9a946b..782aee9 100644
--- a/src/nucfd_deriv_mod.f90
+++ b/src/nucfd_deriv_mod.f90
@@ -28,6 +28,10 @@ subroutine deriv_rhs(f, stencil, x, dfdx)
     type(nucfd_stencil_points) :: stencil_coordinates
 
     real :: a, b, c, d, e
+
+    integer :: offset
+
+    offset = lbound(x, 1) - (-1)
     
     call create_stencil(5, 3, stencil_coordinates)
 
@@ -35,11 +39,7 @@ subroutine deriv_rhs(f, stencil, x, dfdx)
     type is(integer)
        select type(points => stencil_coordinates%stencil)
        type is(real)
-          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))
+          points(:) = x(indices(:) + offset)
        class default
           print *, "Error: Coordinate stencil is misallocated!"
           error stop
@@ -51,9 +51,9 @@ subroutine deriv_rhs(f, stencil, x, dfdx)
        d = coeff_d(stencil_coordinates)
        e = coeff_e(stencil_coordinates)
 
-       dfdx = a * (f(indices(+1)) + (b / a) * f(indices(-1))) &
-            + c * (f(indices(+2)) + (d / c) * f(indices(-2))) &
-            + e * f(indices(0))
+       dfdx = a * (f(indices(+1) + offset) + (b / a) * f(indices(-1) + offset)) &
+            + c * (f(indices(+2) + offset) + (d / c) * f(indices(-2) + offset)) &
+            + e * f(indices(0) + offset)
     class default
        print *, "Error: Index stencil is misallocated!"
        error stop

From 010989bab0a0b06fcb822019cd42be6a38e032fb Mon Sep 17 00:00:00 2001
From: Paul <ptb0890@gmail.com>
Date: Fri, 10 Mar 2023 19:59:04 +0000
Subject: [PATCH 12/22] Initial near working derivative computation
 implementation

---
 .gitignore                       |   6 +-
 examples/compute-derivatives.f90 | 165 +++++++++++++++++++++++++++----
 2 files changed, 153 insertions(+), 18 deletions(-)

diff --git a/.gitignore b/.gitignore
index 331e0b1..b9eb984 100644
--- a/.gitignore
+++ b/.gitignore
@@ -38,4 +38,8 @@ build/
 doc/
 
 # Emacs configuration
-.dir-locals.el
\ No newline at end of file
+.dir-locals.el
+
+# Output
+*.dat
+*.png
\ No newline at end of file
diff --git a/examples/compute-derivatives.f90 b/examples/compute-derivatives.f90
index 2a2e8fa..f22b77d 100644
--- a/examples/compute-derivatives.f90
+++ b/examples/compute-derivatives.f90
@@ -8,15 +8,91 @@ program compute_derivatives
   !! A program that uses non-uniform compact finite difference schemes to compute derivatives on
   !! (non-)uniform meshes.
 
+  use nucfd_types
+  use nucfd_deriv
+  use nucfd_trid_solver
+  
   implicit none
 
+  real, parameter :: pi = 4.0 * atan(1.0)
+  real, parameter :: alpha = 1.0 / 3.0
+  
+  integer :: n
+  real :: L
+  integer, dimension(:), allocatable :: block_spacings
+  logical :: periodic
   real, dimension(:), allocatable :: x ! 1D "grid"
+
+  type(nucfd_index_stencil) :: stencil
+  integer :: width, centre
+
+  integer :: idx, i
+
+  real, dimension(:), allocatable :: f
+  real, dimension(:), allocatable :: dfdx
+
+  real, dimension(:), allocatable :: a, b, c, r
+  
+  n = 66
+  L = 1.0
+  block_spacings = (/ 1, 2, 1 /)
+  periodic = .true.
+  call build_grid(n, L, block_spacings, periodic, x)
+
+  print *, "+++ Built grid +++"
+  print *, "- lbounds: ", lbound(x)
+  print *, "- ubounds: ", ubound(x)
+  print *, "- 1st/last node: ", x(1), x(n)
+  if (periodic) then
+     print *, "- periodic node: ", x(n + 1)
+  end if
+
+  allocate(f, mold=x)
+  f(:) = sin((x(:) / L) * (2.0 * pi))
   
-  call build_grid(11, 1.0, (/ 1, 2 /), x)
+  allocate(a(n))
+  allocate(b(n))
+  allocate(c(n))
+  allocate(r(n))
+  allocate(dfdx(n))
+  
+  width = 5
+  centre = 3
+  call create_stencil(width, centre, stencil)
+
+  print *, "+++ Computing RHS +++"
+  select type (indices => stencil%stencil)
+  type is(integer)
+     do idx = 1, n
+        ! Move stencil
+        do i = lbound(indices, 1), ubound(indices, 1)
+           indices(i) = idx + i 
+        end do
+        
+        call deriv_rhs(f, stencil, x, r(idx))
+     end do
+  end select
+
+  print *, "+++ Solving system +++"
+  a(:) = alpha
+  b(:) = 1.0
+  c(:) = alpha
+  call solve_cyclic(a, b, c, r, dfdx)
+
+  open (unit=10, file="dfdx.dat")
+  do i = 1, n
+     write(10, *) x(i), (2.0 * pi / L) * cos((x(i) / L) * (2.0 * pi)), dfdx(i)
+  end do
+  close(10)
+
+  deallocate(f, dfdx)
+  deallocate(x)
+
+  deallocate(a, b, c, r)
   
 contains
 
-  subroutine build_grid(n, L, spacings, x)
+  subroutine build_grid(n, L, spacings, periodic, x)
     !! Crude subroutine to build a 1-D grid.
     !! The grid is described by the number of nodes, length and an array of relative grid spacings
     !! for each block (constant within a block), returning a 1-D array of node centres.
@@ -28,6 +104,7 @@ subroutine build_grid(n, L, spacings, x)
     real, intent(in) :: L                             !! The length of the domain
     integer, dimension(:), intent(in) :: spacings     !! An array of the relative (integer) spacings
                                                       !! for each block
+    logical, intent(in) :: periodic                   !! Flag indicating periodicity of mesh
     real, dimension(:), allocatable, intent(out) :: x !! The grid node locations
 
     integer :: ncells
@@ -40,32 +117,80 @@ subroutine build_grid(n, L, spacings, x)
     integer :: b
     integer :: idx
 
-    allocate(x(n))
+    integer, parameter :: lghost = 2 ! Ghost node count at left boundary
+    integer, parameter :: rghost = 2 ! Ghost node count at right boundary
 
-    ncells = n - 1
+    allocate(x(1-lghost:n+rghost))
+
+    if (.not. periodic) then
+       ncells = n - 1
+    else
+       ncells = n
+    end if
     nblocks = size(spacings, 1)
-    block_size = ncells / nblocks
-    if (block_size * nblocks + 1 /= n) then
-       print *, "Error: invalid grid specification!"
-       print *, " - Grid size n should be specified so that n - 1 is divisible by number of blocks"
-       print *, "   e.g. for 2 blocks a grid size of 11 is valid, 10 is not."
-       error stop
+    if (.not. periodic) then
+       block_size = ncells / nblocks
+    else
+       block_size = ncells / nblocks
+    end if
+    if (.not. periodic) then
+       if (block_size * nblocks + 1 /= n) then
+          print *, "Error: invalid grid specification!"
+          print *, " - Grid size n should be specified so that n - 1 is divisible by number of blocks"
+          print *, "   e.g. for 2 blocks a grid size of 11 is valid, 10 is not."
+          error stop
+       end if
+    else
+       if (block_size * nblocks /= n) then
+          print *, "Error: invalid grid specification!"
+          error stop
+       end if
     end if
     block_length = L / real(sum(spacings))
     h = block_length / real(block_size)
-
+    
     idx = 1
     x(idx) = 0.0
-    idx = 2
     do b = 1, nblocks
+       if (periodic .and. (b == nblocks)) then
+          block_size = block_size - 1
+       end if
        do i = 1, block_size
-          x(idx) = x(idx - 1) + real(spacings(b)) * h
           idx = idx + 1
+          x(idx) = x(idx - 1) + real(spacings(b)) * h
        end do
     end do
+    if (idx /= n) then
+       print *, "Error: didn't fill all interior nodes!"
+       error stop
+    end if
+
+    ! Add ghost points
+    if (.not. periodic) then
+       do i = 1, lghost
+          idx = 1 - i
+          x(idx) = x(idx + 1) - real(spacings(1)) * h
+       end do
+       do i = 1, rghost
+          idx = n + i
+          x(idx) = x(idx - 1) + real(spacings(nblocks)) * h
+       end do
+    else
+       idx = n + 1
+       x(idx) = L
+       do i = 2, rghost
+          idx = n + i
+          x(idx) = x(idx - 1) + (x(i) - x(i - 1))
+       end do
+
+       do i = 1, lghost
+          idx = 1 - i
+          x(idx) = x(idx + 1) - (x((n + 1) - (i - 1)) - (x((n + 1) - (i - 1) - 1)))
+       end do
+    end if
 
     ! Self-check
-    if (any(x > L) .or. any(x < 0.0)) then
+    if (any(x(1:n) > L) .or. any(x(1:n) < 0.0)) then
        print *, "Error: grid exceeds specified length!"
        print *, x
        error stop
@@ -75,11 +200,17 @@ subroutine build_grid(n, L, spacings, x)
        print *, x(1)
        error stop
     end if
-    if (abs(x(n) - L) > (2 * epsilon(L) * L)) then
+
+    if (.not. periodic) then
+       idx = n
+    else
+       idx = n + 1
+    end if
+    if (abs(x(idx) - L) > (2 * epsilon(L) * L)) then
        print *, "Error: end point is wrong!"
-       print *, "- Computed: ", x(n)
+       print *, "- Computed: ", x(idx)
        print *, "- Expected: ", L
-       print *, "- Relative error: ", abs(x(n) - L) / L
+       print *, "- Relative error: ", abs(x(idx) - L) / L
        error stop
     end if
     

From dce6c3e5426a23a293073bdf2387e6959b243ff1 Mon Sep 17 00:00:00 2001
From: Paul <ptb0890@gmail.com>
Date: Sat, 11 Mar 2023 20:33:52 +0000
Subject: [PATCH 13/22] Add differentiation rules test for quadratic functions

---
 .../differentiation_rules/CMakeLists.txt      |  2 +
 .../quadratic_function.f90                    | 61 +++++++++++++++++++
 2 files changed, 63 insertions(+)
 create mode 100644 tests/fd-schemes/differentiation_rules/quadratic_function.f90

diff --git a/tests/fd-schemes/differentiation_rules/CMakeLists.txt b/tests/fd-schemes/differentiation_rules/CMakeLists.txt
index e294d04..920f892 100644
--- a/tests/fd-schemes/differentiation_rules/CMakeLists.txt
+++ b/tests/fd-schemes/differentiation_rules/CMakeLists.txt
@@ -17,3 +17,5 @@ define_test(fd-schemes linear_rising_function)
 add_test_libs(fd-schemes linear_rising_function diff_rules_utils)
 define_test(fd-schemes linear_falling_function)
 add_test_libs(fd-schemes linear_falling_function diff_rules_utils)
+define_test(fd-schemes quadratic_function)
+add_test_libs(fd-schemes quadratic_function diff_rules_utils)
diff --git a/tests/fd-schemes/differentiation_rules/quadratic_function.f90 b/tests/fd-schemes/differentiation_rules/quadratic_function.f90
new file mode 100644
index 0000000..218be94
--- /dev/null
+++ b/tests/fd-schemes/differentiation_rules/quadratic_function.f90
@@ -0,0 +1,61 @@
+! tests/fd-schemes/differentiation_rules/quadratic_function.f90
+!
+!! Part of the fd-schemes test suite.
+!
+! SPDX-License-Identifier: BSD-3-Clause
+
+program quadratic_function
+  !! Tests that rules of differentiation are respected for quadratic functions.
+
+  use nucfd_types
+  use nucfd_deriv
+  
+  use nucfd_tests
+  use diff_rules_utils
+  
+  implicit none
+
+  integer :: n ! Grid size
+  real :: L    ! Domain size
+
+  real, dimension(:), allocatable :: x ! Grid
+  real, dimension(:), allocatable :: f ! Function
+  real :: dfdx ! Derivative
+  real :: dgdx ! Derivative
+
+  type(nucfd_index_stencil) :: stencil
+  type(test_setup) :: ts
+
+  integer :: i
+  
+  call initialise_suite("Quadratic function differentiation rules")
+
+  n = 33
+  L = 1.0
+
+  ts = initialise_test(n, L)
+  f = allocate_test_array(ts)
+  x = create_mesh_array(ts)
+  stencil = build_stencil(ts)
+
+  do i = 1, n
+     f(i) = x(i)**2
+  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 quadratic_function

From d8d8f8f8ffbf6858055fbc130aeac9301db41ed3 Mon Sep 17 00:00:00 2001
From: Paul <ptb0890@gmail.com>
Date: Sat, 11 Mar 2023 20:34:17 +0000
Subject: [PATCH 14/22] Change differentiation rules to run on refined mesh

---
 tests/fd-schemes/differentiation_rules/diff_rules_utils.f90 | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

diff --git a/tests/fd-schemes/differentiation_rules/diff_rules_utils.f90 b/tests/fd-schemes/differentiation_rules/diff_rules_utils.f90
index b4ce3bb..94be47e 100644
--- a/tests/fd-schemes/differentiation_rules/diff_rules_utils.f90
+++ b/tests/fd-schemes/differentiation_rules/diff_rules_utils.f90
@@ -61,7 +61,11 @@ function create_mesh_array(ts) result(x)
 
     associate(n => ts%n, h => ts%h)
       do i = 1, n
-         x(i) = real(i - 1) * h
+         if (i < n / 2) then
+            x(i) = real(i - 1) * h
+         else
+            x(i) = x(i - 1) + h / 2.0
+         end if
       end do
     end associate
 

From 27d86db3574210a891a1b5ce2973a402f5557922 Mon Sep 17 00:00:00 2001
From: Paul <ptb0890@gmail.com>
Date: Sat, 11 Mar 2023 21:05:00 +0000
Subject: [PATCH 15/22] Add test of RHS on variable mesh (currently failing)

---
 .../differentiation_rules/CMakeLists.txt      |  2 +
 .../varmesh_constgrad.f90                     | 86 +++++++++++++++++++
 2 files changed, 88 insertions(+)
 create mode 100644 tests/fd-schemes/differentiation_rules/varmesh_constgrad.f90

diff --git a/tests/fd-schemes/differentiation_rules/CMakeLists.txt b/tests/fd-schemes/differentiation_rules/CMakeLists.txt
index 920f892..d3234ea 100644
--- a/tests/fd-schemes/differentiation_rules/CMakeLists.txt
+++ b/tests/fd-schemes/differentiation_rules/CMakeLists.txt
@@ -19,3 +19,5 @@ define_test(fd-schemes linear_falling_function)
 add_test_libs(fd-schemes linear_falling_function diff_rules_utils)
 define_test(fd-schemes quadratic_function)
 add_test_libs(fd-schemes quadratic_function diff_rules_utils)
+define_test(fd-schemes varmesh_constgrad)
+add_test_libs(fd-schemes varmesh_constgrad diff_rules_utils)
diff --git a/tests/fd-schemes/differentiation_rules/varmesh_constgrad.f90 b/tests/fd-schemes/differentiation_rules/varmesh_constgrad.f90
new file mode 100644
index 0000000..e4d571b
--- /dev/null
+++ b/tests/fd-schemes/differentiation_rules/varmesh_constgrad.f90
@@ -0,0 +1,86 @@
+! tests/fd-schemes/differentiation_rules/varmesh_constgrad.f90
+!
+!! Part of the fd-schemes test suite.
+!
+! SPDX-License-Identifier: BSD-3-Clause
+
+program linear_rising_function
+  !! Tests that the derivative of a linear function remains constant on a mesh with local
+  !! refinement.
+
+  use nucfd_types
+  use nucfd_deriv
+  
+  use nucfd_tests
+  use diff_rules_utils
+  
+  implicit none
+
+  integer :: n ! Grid size
+  real :: L    ! Domain size
+
+  real, dimension(:), allocatable :: x ! Grid
+  real, dimension(:), allocatable :: f ! Function
+  real :: dfdx ! Derivative
+  real :: dgdx ! Derivative
+
+  type(nucfd_index_stencil) :: stencil
+  type(test_setup) :: ts
+
+  integer :: i
+
+  logical :: passing
+  
+  call initialise_suite("Linearly increasing function differentiation rules")
+
+  n = 33
+  L = 1.0
+
+  ts = initialise_test(n, L)
+  f = allocate_test_array(ts)
+  x = create_mesh_array(ts)
+  stencil = build_stencil(ts)
+
+  f(1) = 0.0
+  dfdx = 1.0
+  do i = 2, n
+     f(i) = f(i - 1) + (x(i) - x(i - 1)) * dfdx
+  end do
+
+  i = 3
+  call centre_stencil(i, stencil)
+  call deriv_rhs(f, stencil, x, dfdx)
+
+  passing = .true.
+  do i = 4, n - 2
+     print *, i, x(i - 1) - x(i - 2), x(i) - x(i - 1), x(i + 1) - x(i), x(i + 2) - x(i + 1)
+     call centre_stencil(i, stencil)
+     call deriv_rhs(f, stencil, x, dgdx)
+     passing = passing .and. check_scalar(dgdx, dfdx)
+  end do
+  call test_report("Variable mesh, constant gradient", passing)
+  
+  deallocate(x)
+  deallocate(f)
+  
+  call finalise_suite()
+
+contains
+
+  subroutine centre_stencil(idx, stencil)
+    !! Move a stencil to be centred at specified index.
+
+    integer, intent(in) :: idx                          !! The index
+    type(nucfd_index_stencil), intent(inout) :: stencil !! The stencil to be moved
+
+    integer :: shift
+
+    select type(indices => stencil%stencil)
+    type is(integer)
+       shift = idx - indices(0)
+       indices(:) = indices(:) + shift
+    end select
+    
+  end subroutine centre_stencil
+ 
+end program linear_rising_function

From 666b7655f83b2e80b11bc93ca00c54cbe7441997 Mon Sep 17 00:00:00 2001
From: Paul <ptb0890@gmail.com>
Date: Sat, 11 Mar 2023 21:29:02 +0000
Subject: [PATCH 16/22] Test on variable meshes wasn't accounting for need to
 offset ghost nodes

---
 .../differentiation_rules/varmesh_constgrad.f90     | 13 +++++++------
 tests/fd-schemes/verify_coeffs/verify_coeff_a.f90   |  4 ++--
 2 files changed, 9 insertions(+), 8 deletions(-)

diff --git a/tests/fd-schemes/differentiation_rules/varmesh_constgrad.f90 b/tests/fd-schemes/differentiation_rules/varmesh_constgrad.f90
index e4d571b..aea4db2 100644
--- a/tests/fd-schemes/differentiation_rules/varmesh_constgrad.f90
+++ b/tests/fd-schemes/differentiation_rules/varmesh_constgrad.f90
@@ -27,11 +27,11 @@ program linear_rising_function
   type(nucfd_index_stencil) :: stencil
   type(test_setup) :: ts
 
-  integer :: i
+  integer :: i, idx
 
   logical :: passing
   
-  call initialise_suite("Linearly increasing function differentiation rules")
+  call initialise_suite("Constant gradient on locally refined mesh")
 
   n = 33
   L = 1.0
@@ -47,18 +47,19 @@ program linear_rising_function
      f(i) = f(i - 1) + (x(i) - x(i - 1)) * dfdx
   end do
 
-  i = 3
+  i = 1
   call centre_stencil(i, stencil)
   call deriv_rhs(f, stencil, x, dfdx)
 
   passing = .true.
-  do i = 4, n - 2
-     print *, i, x(i - 1) - x(i - 2), x(i) - x(i - 1), x(i + 1) - x(i), x(i + 2) - x(i + 1)
+  do i = 2, n - 4
+     idx = i + 2
+     print *, idx, x(idx - 1) - x(idx - 2), x(idx) - x(idx - 1), x(idx + 1) - x(idx), x(idx + 2) - x(idx + 1)
      call centre_stencil(i, stencil)
      call deriv_rhs(f, stencil, x, dgdx)
      passing = passing .and. check_scalar(dgdx, dfdx)
   end do
-  call test_report("Variable mesh, constant gradient", passing)
+  call test_report("Gradient remains constant", passing)
   
   deallocate(x)
   deallocate(f)
diff --git a/tests/fd-schemes/verify_coeffs/verify_coeff_a.f90 b/tests/fd-schemes/verify_coeffs/verify_coeff_a.f90
index 9c3e5b4..b29d062 100644
--- a/tests/fd-schemes/verify_coeffs/verify_coeff_a.f90
+++ b/tests/fd-schemes/verify_coeffs/verify_coeff_a.f90
@@ -4,7 +4,7 @@
 !
 ! SPDX-License-Identifier: BSD-3-Clause
 
-program verify_coeff_b
+program verify_coeff_a
   !! Tests the computation of finite difference coefficient for non-uniform grids acting on f_{i+1}.
 
   use nucfd_types
@@ -58,4 +58,4 @@ program verify_coeff_b
 
   call finalise_suite()
 
-end program verify_coeff_b
+end program verify_coeff_a

From a3f24d38de2faeacf042e3cb1b91293bde8a908a Mon Sep 17 00:00:00 2001
From: Paul <ptb0890@gmail.com>
Date: Sat, 11 Mar 2023 22:15:31 +0000
Subject: [PATCH 17/22] Update coefficient test suite

---
 tests/fd-schemes/verify_coeffs/CMakeLists.txt |  1 +
 .../verify_coeffs/verify_coeff_e.f90          |  2 +-
 .../verify_coeffs/verify_coeff_pairsums.f90   | 55 +++++++++++++++++++
 3 files changed, 57 insertions(+), 1 deletion(-)
 create mode 100644 tests/fd-schemes/verify_coeffs/verify_coeff_pairsums.f90

diff --git a/tests/fd-schemes/verify_coeffs/CMakeLists.txt b/tests/fd-schemes/verify_coeffs/CMakeLists.txt
index f95f68c..26e3d3f 100644
--- a/tests/fd-schemes/verify_coeffs/CMakeLists.txt
+++ b/tests/fd-schemes/verify_coeffs/CMakeLists.txt
@@ -14,3 +14,4 @@ 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)
+define_test(fd-schemes verify_coeff_pairsums)
diff --git a/tests/fd-schemes/verify_coeffs/verify_coeff_e.f90 b/tests/fd-schemes/verify_coeffs/verify_coeff_e.f90
index 9f5f75a..55409e0 100644
--- a/tests/fd-schemes/verify_coeffs/verify_coeff_e.f90
+++ b/tests/fd-schemes/verify_coeffs/verify_coeff_e.f90
@@ -43,7 +43,7 @@ program verify_coeff_e
   end select
 
   e = coeff_e(stencil)
-  call test_report("Coefficient E", check_scalar(e, eref / (4.0 * h)))
+  call test_report("Coefficient E", check_scalar(e, eref))
 
   call finalise_suite()
 
diff --git a/tests/fd-schemes/verify_coeffs/verify_coeff_pairsums.f90 b/tests/fd-schemes/verify_coeffs/verify_coeff_pairsums.f90
new file mode 100644
index 0000000..609c39a
--- /dev/null
+++ b/tests/fd-schemes/verify_coeffs/verify_coeff_pairsums.f90
@@ -0,0 +1,55 @@
+! tests/fd-schemes/verify_coeffs/verify_coeff_pairsums.f90
+!
+!! Part of the fd-schemes test suite.
+!
+! SPDX-License-Identifier: BSD-3-Clause
+
+program verify_coeff_pairsums
+  !! Tests the computation of finite difference coefficients for the pairs f_{i+/-1} and 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 :: a, b, c, d
+  real, parameter :: abref = 0.0, cdref = 0.0
+  
+  call initialise_suite("Verify coefficient pairsums")
+
+  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
+
+  a = coeff_a(stencil)
+  b = coeff_b(stencil)
+  call test_report("Coefficient A+B", check_scalar(a + b, abref))
+
+  c = coeff_c(stencil)
+  d = coeff_d(stencil)
+  call test_report("Coefficient A+B", check_scalar(c + d, cdref))
+
+  call finalise_suite()
+
+end program verify_coeff_pairsums

From b5a242ab7e13481a33f049166024e235c1d1129a Mon Sep 17 00:00:00 2001
From: Paul Bartholomew <p.bartholomew@epcc.ed.ac.uk>
Date: Mon, 20 Mar 2023 12:04:45 +0000
Subject: [PATCH 18/22] Update CHANGELOG

---
 CHANGELOG.md | 1 +
 1 file changed, 1 insertion(+)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index c180acf..fbe858c 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -6,6 +6,7 @@ This is the change log for `NuCFD`, it is based on the format suggested by [Keep
 
 ### Added
 
+- Added example code to compute derivative of a function [010989b].
 - Added `FORD` documentation system [c780fe1].
 - Added tridiagonal solver with support for periodic problems [f5d254f].
 - Added simple library to support writing tests [a5814c8].

From dc7fc46860d13ac76a5549f115ea80a468116422 Mon Sep 17 00:00:00 2001
From: Paul Bartholomew <p.bartholomew@epcc.ed.ac.uk>
Date: Mon, 20 Mar 2023 18:07:03 +0000
Subject: [PATCH 19/22] Add test that A coefficient is independent of scale

---
 .../verify_coeffs/verify_coeff_a.f90          | 28 +++++++++++++++++++
 1 file changed, 28 insertions(+)

diff --git a/tests/fd-schemes/verify_coeffs/verify_coeff_a.f90 b/tests/fd-schemes/verify_coeffs/verify_coeff_a.f90
index b29d062..cc41ce6 100644
--- a/tests/fd-schemes/verify_coeffs/verify_coeff_a.f90
+++ b/tests/fd-schemes/verify_coeffs/verify_coeff_a.f90
@@ -28,6 +28,9 @@ program verify_coeff_a
   real, parameter :: numerator_corr_f1ref = 0.0
   real, parameter :: denominator_f1ref = 3.0
   real, parameter :: divisor_f1ref = 2.0
+
+  real :: a_2h
+  real :: numerator_2h, numerator_corr_2h, denominator_2h, divisor_2h
   
   call initialise_suite("Verify coefficient A")
 
@@ -56,6 +59,31 @@ program verify_coeff_a
   a = coeff_a(stencil)
   call test_report("Coefficient A", check_scalar(a, aref / (2.0 * h)))
 
+  !! Check coefficient is independent of scale
+  select type(points => stencil%stencil)
+  type is(real)
+     points(:) = 0.0
+     points(-2) = -2.0 * (2 * h)
+     points(-1) = -1.0 * (2 * h)
+     points(0) = 0.0
+     points(1) = +1.0 * (2 * h)
+     points(2) = +2.0 * (2 * h)
+  class default
+     print *, "Shouldn't happen!"
+  end select
+  call coeff_a_components(points_to_deltas(stencil), numerator_2h, numerator_corr_2h, &
+       denominator_2h, divisor_2h)
+  call test_report("Coefficient A numerator (scale independence)", check_scalar(numerator / (h**3), &
+       numerator_2h / ((2 * h)**3)))
+  call test_report("Coefficient A numerator correction (scale independence)", check_scalar(numerator_corr / (h**3), &
+       numerator_corr_2h / ((2 * h)**3)))
+  call test_report("Coefficient A denominator (scale independence)", check_scalar(denominator / (h**3), &
+       denominator_2h / ((2 * h)**3)))
+  call test_report("Coefficient A divisor (scale independence)", check_scalar(divisor / h, &
+       divisor_2h / (2 * h)))
+  a_2h = coeff_a(stencil)
+  call test_report("Coefficient A (scale independence)", check_scalar(a * h, a_2h * (2 * h)))
+
   call finalise_suite()
 
 end program verify_coeff_a

From 282b7ce9e2cad542e159a82922797c280dd7b227 Mon Sep 17 00:00:00 2001
From: Paul <ptb0890@gmail.com>
Date: Sat, 15 Jul 2023 17:15:47 +0100
Subject: [PATCH 20/22] Correct typo

---
 src/coeffs/nucfd_coeffs_a.f90 | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/coeffs/nucfd_coeffs_a.f90 b/src/coeffs/nucfd_coeffs_a.f90
index 0d6ce03..46a0cfd 100644
--- a/src/coeffs/nucfd_coeffs_a.f90
+++ b/src/coeffs/nucfd_coeffs_a.f90
@@ -102,7 +102,7 @@ module subroutine coeff_a_components(h, numerator, numerator_corr, denominator,
     end select
 
     associate(beta => alpha) ! To match Gamet et al. (1999)
-      numerator = coeff_numerator(hm1, h0, hp2, hp2, beta)
+      numerator = coeff_numerator(hm1, h0, hp1, 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)

From 0ed2563e9a55d897a7eff3307a7e804eeea63d1a Mon Sep 17 00:00:00 2001
From: Paul <ptb0890@gmail.com>
Date: Sat, 15 Jul 2023 17:49:04 +0100
Subject: [PATCH 21/22] Simplify coeff A verification test

---
 .../verify_coeffs/verify_coeff_a.f90          | 137 +++++++++++-------
 1 file changed, 85 insertions(+), 52 deletions(-)

diff --git a/tests/fd-schemes/verify_coeffs/verify_coeff_a.f90 b/tests/fd-schemes/verify_coeffs/verify_coeff_a.f90
index cc41ce6..f12f589 100644
--- a/tests/fd-schemes/verify_coeffs/verify_coeff_a.f90
+++ b/tests/fd-schemes/verify_coeffs/verify_coeff_a.f90
@@ -21,16 +21,10 @@ program verify_coeff_a
   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
-
-  real :: a_2h
-  real :: numerator_2h, numerator_corr_2h, denominator_2h, divisor_2h
+  integer :: scale
+  real :: aref
+  real :: numerator_f1ref, numerator_corr_f1ref, denominator_f1ref, divisor_f1ref
   
   call initialise_suite("Verify coefficient A")
 
@@ -39,51 +33,90 @@ program verify_coeff_a
   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)))
 
-  !! Check coefficient is independent of scale
-  select type(points => stencil%stencil)
-  type is(real)
-     points(:) = 0.0
-     points(-2) = -2.0 * (2 * h)
-     points(-1) = -1.0 * (2 * h)
-     points(0) = 0.0
-     points(1) = +1.0 * (2 * h)
-     points(2) = +2.0 * (2 * h)
-  class default
-     print *, "Shouldn't happen!"
-  end select
-  call coeff_a_components(points_to_deltas(stencil), numerator_2h, numerator_corr_2h, &
-       denominator_2h, divisor_2h)
-  call test_report("Coefficient A numerator (scale independence)", check_scalar(numerator / (h**3), &
-       numerator_2h / ((2 * h)**3)))
-  call test_report("Coefficient A numerator correction (scale independence)", check_scalar(numerator_corr / (h**3), &
-       numerator_corr_2h / ((2 * h)**3)))
-  call test_report("Coefficient A denominator (scale independence)", check_scalar(denominator / (h**3), &
-       denominator_2h / ((2 * h)**3)))
-  call test_report("Coefficient A divisor (scale independence)", check_scalar(divisor / h, &
-       divisor_2h / (2 * h)))
-  a_2h = coeff_a(stencil)
-  call test_report("Coefficient A (scale independence)", check_scalar(a * h, a_2h * (2 * h)))
+  scale = 1
+  call setup_stencil(scale, scale, scale, scale)
+  call compute_uniform_ref_coeffs(scale)
+  call test_coefficient("Uniform", stencil, aref, numerator_f1ref, numerator_corr_f1ref, &
+       denominator_f1ref, divisor_f1ref)
 
+  !! Check coefficient is independent of scale
+  scale = 2
+  call setup_stencil(scale, scale, scale, scale)
+  call compute_uniform_ref_coeffs(scale)
+  call test_coefficient("Uniform (2h)", stencil, aref, numerator_f1ref, numerator_corr_f1ref, &
+       denominator_f1ref, divisor_f1ref)
+  
   call finalise_suite()
 
+contains
+
+  subroutine setup_stencil(sm2, sm1, sp1, sp2)
+
+    integer, intent(in) :: sm2
+    integer, intent(in) :: sm1
+    integer, intent(in) :: sp1
+    integer, intent(in) :: sp2
+
+    select type(points => stencil%stencil)
+    type is(real)
+       points(:) = 0.0
+       points(-2) = -(sm1 + sm2) * h
+       points(-1) = -sm1 * h
+       points(0) = 0.0
+       points(1) = +sp1 * h
+       points(2) = +(sp1 + sp2) * h
+    class default
+       print *, "Error: Coordinate stencil is misallocated!"
+    end select
+    
+  end subroutine setup_stencil
+
+  subroutine compute_uniform_ref_coeffs(scale)
+
+    integer, intent(in) :: scale
+    
+    aref = 14.0 / 9.0 / (2.0 * (scale * h))
+    numerator_f1ref = 14.0 / 3.0 * ((scale * h)**3)
+    numerator_corr_f1ref = 0.0 * ((scale * h)**3)
+    denominator_f1ref = 3.0 * ((scale * h)**3)
+    divisor_f1ref = 2.0 * (scale * h)
+  end subroutine compute_uniform_ref_coeffs
+  
+  subroutine test_coefficient(grid_desc, stencil, aref, numerator_f1ref, numerator_corr_f1ref, &
+       denominator_f1ref, divisor_f1ref)
+
+    character(len=*), intent(in) :: grid_desc
+    type(nucfd_stencil_points), intent(in) :: stencil
+    real, intent(in) :: aref
+    real, intent(in) :: numerator_f1ref
+    real, intent(in) :: numerator_corr_f1ref
+    real, intent(in) :: denominator_f1ref
+    real, intent(in) :: divisor_f1ref
+
+    character(len=:), allocatable :: prefix
+    real :: numerator, numerator_corr, denominator, divisor
+
+    prefix = "("//grid_desc//") "
+    
+    call coeff_a_components(points_to_deltas(stencil), numerator, numerator_corr, &
+         denominator, divisor)
+    call test_report(prefix//"Coefficient A numerator", &
+         check_scalar(numerator, numerator_f1ref))
+    call test_report(prefix//"Coefficient A numerator correction", &
+         check_scalar(numerator_corr, numerator_corr_f1ref))
+    call test_report(prefix//"Coefficient A denominator", &
+         check_scalar(denominator, denominator_f1ref))
+    call test_report(prefix//"Coefficient A divisor", &
+         check_scalar(divisor, divisor_f1ref))
+    call test_report(prefix//"Coefficient A den * div", &
+         check_scalar(denominator * divisor, &
+         (denominator_f1ref * divisor_f1ref)))
+    
+    a = coeff_a(stencil)
+    call test_report("Coefficient A", &
+         check_scalar(a, aref))
+
+  end subroutine test_coefficient
+  
 end program verify_coeff_a

From 090cebb704d24dec6e1cac8cc64454438fde7e85 Mon Sep 17 00:00:00 2001
From: Paul <ptb0890@gmail.com>
Date: Sun, 16 Jul 2023 21:51:00 +0100
Subject: [PATCH 22/22] Rearrange denominator/divisor expressions

The way these had been written wasn't wise, too many FP operations
with potential for numerical errors. Also, the divisors should be
2h / 4h for the i+/-1 and i+/-2 cases, respectively which is the
case for the new arrangement
---
 src/coeffs/nucfd_coeffs_a.f90 | 13 +++++++------
 src/coeffs/nucfd_coeffs_b.f90 | 14 +++++++-------
 src/coeffs/nucfd_coeffs_c.f90 | 16 ++++++++--------
 src/coeffs/nucfd_coeffs_d.f90 | 16 ++++++++--------
 4 files changed, 30 insertions(+), 29 deletions(-)

diff --git a/src/coeffs/nucfd_coeffs_a.f90 b/src/coeffs/nucfd_coeffs_a.f90
index 46a0cfd..98780ef 100644
--- a/src/coeffs/nucfd_coeffs_a.f90
+++ b/src/coeffs/nucfd_coeffs_a.f90
@@ -104,8 +104,8 @@ module subroutine coeff_a_components(h, numerator, numerator_corr, denominator,
     associate(beta => alpha) ! To match Gamet et al. (1999)
       numerator = coeff_numerator(hm1, h0, hp1, 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)
+      denominator = coeff_denominator(hm1, h0, hp1, hp2)
+      divisor = coeff_divisor(h0, hp1)
     end associate
 
   end subroutine coeff_a_components
@@ -145,29 +145,30 @@ pure real function coeff_numerator_corr(hm1, h0, hp1, hp2, alpha, beta)
   end function coeff_numerator_corr
 
   pure real function coeff_denominator(h0, hp1, hp2)
+  pure real function coeff_denominator(hm1, h0, hp1, hp2)
     !! Computes the denominator of the coefficient acting on f_{i+1}.
     !!
     !! Reduces to 3 h^3 when h=const. Dividing the numerator by this term yields the coefficient
     !! 14/9 when h=const, alpha=beta=1/3.
 
+    real, intent(in) :: hm1 
     real, intent(in) :: h0 
     real, intent(in) :: hp1 
     real, intent(in) :: hp2 
 
-    coeff_denominator = (3.0 * hp1 * ((h0 + hp1) / 2.0) * hp2)
+    coeff_denominator = (hp1 * (hm1 + h0 + hp1) * hp2)
   end function coeff_denominator
 
-  pure real function coeff_divisor(hm1, h0, hp1)
+  pure real function coeff_divisor(h0, hp1)
     !! Computes the non-uniform equivalent to 2h divisor of the coefficient acting on f_{i+1}.
     !!
     !! Dividing the coefficient by this term should reduce to (14/9)/(2h) when h=const,
     !! alpha=beta=1/3.
     
-    real, intent(in) :: hm1 
     real, intent(in) :: h0 
     real, intent(in) :: hp1 
 
-    coeff_divisor = (2.0 * ((hm1 + h0 + hp1) / 3.0))
+    coeff_divisor = (h0 + hp1)
   end function coeff_divisor
   
 end submodule nucfd_coeffs_a
diff --git a/src/coeffs/nucfd_coeffs_b.f90 b/src/coeffs/nucfd_coeffs_b.f90
index 2e104cc..c7b092b 100644
--- a/src/coeffs/nucfd_coeffs_b.f90
+++ b/src/coeffs/nucfd_coeffs_b.f90
@@ -103,8 +103,8 @@ module subroutine coeff_b_components(h, numerator, numerator_corr, denominator,
     associate(beta => alpha) ! To match Gamet et al. (1999)
       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)
+      denominator = coeff_denominator(hm1, h0, hp1, hp2)
+      divisor = coeff_divisor(h0, hp1)
     end associate
 
   end subroutine coeff_b_components
@@ -143,7 +143,7 @@ pure real function coeff_numerator_corr(hm1, h0, hp1, hp2, alpha, beta)
                                                                     ! alpha = beta
   end function coeff_numerator_corr
 
-  pure real function coeff_denominator(hm1, h0, hp1)
+  pure real function coeff_denominator(hm1, 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
@@ -152,11 +152,12 @@ pure real function coeff_denominator(hm1, h0, hp1)
     real, intent(in) :: hm1
     real, intent(in) :: h0 
     real, intent(in) :: hp1
+    real, intent(in) :: hp2
     
-    coeff_denominator = 3.0 * hm1 * h0 * ((h0 + hp1) / 2.0)
+    coeff_denominator = hm1 * h0 * (h0 + hp1 + hp2)
   end function coeff_denominator
 
-  pure real function coeff_divisor(h0, hp1, hp2)
+  pure real function coeff_divisor(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,
@@ -164,9 +165,8 @@ pure real function coeff_divisor(h0, hp1, hp2)
 
     real, intent(in) :: h0
     real, intent(in) :: hp1
-    real, intent(in) :: hp2
 
-    coeff_divisor = 2.0 * (h0 + hp1 + hp2) / 3.0
+    coeff_divisor = (h0 + hp1) 
   end function coeff_divisor
   
 end submodule nucfd_coeffs_b
diff --git a/src/coeffs/nucfd_coeffs_c.f90 b/src/coeffs/nucfd_coeffs_c.f90
index 256a716..e90416a 100644
--- a/src/coeffs/nucfd_coeffs_c.f90
+++ b/src/coeffs/nucfd_coeffs_c.f90
@@ -106,8 +106,8 @@ module subroutine coeff_c_components(h, numerator, denominator, divisor)
 
     associate(beta => alpha) ! To match Gamet et al. (1999)
       numerator = coeff_numerator(hm1, h0, hp1, alpha, beta)
-      denominator = coeff_denominator(hm1, h0, hp1, hp2)
-      divisor = coeff_divisor(h0, hp1, hp2)
+      denominator = coeff_denominator(h0, hp1, hp2)
+      divisor = coeff_divisor(hm1, h0, hp1, hp2)
     end associate
     
   end subroutine coeff_c_components
@@ -129,33 +129,33 @@ pure real function coeff_numerator(hm1, h0, hp1, alpha, beta)
     
   end function coeff_numerator
   
-  pure real function coeff_denominator(hm1, h0, hp1, hp2)
+  pure real function coeff_denominator(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)
+    coeff_denominator = hp2 * (hp1 + hp2) &
+         * (h0 + hp1 + hp2)
     
   end function coeff_denominator
 
-  pure real function coeff_divisor(h0, hp1, hp2)
+  pure real function coeff_divisor(hm1, 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) :: hm1
     real, intent(in) :: h0
     real, intent(in) :: hp1
     real, intent(in) :: hp2
 
-    coeff_divisor = 4.0 * ((h0 + hp1 + hp2) / 3.0) 
+    coeff_divisor = (hm1 + h0 + hp1 + hp2) 
   end function coeff_divisor
   
 end submodule nucfd_coeffs_c
diff --git a/src/coeffs/nucfd_coeffs_d.f90 b/src/coeffs/nucfd_coeffs_d.f90
index 9dd0f80..b2d347d 100644
--- a/src/coeffs/nucfd_coeffs_d.f90
+++ b/src/coeffs/nucfd_coeffs_d.f90
@@ -106,8 +106,8 @@ module subroutine coeff_d_components(h, numerator, denominator, divisor)
 
     associate(beta => alpha) ! To match Gamet et al. (1999)
       numerator = coeff_numerator(h0, hp1, hp2, alpha, beta)
-      denominator = coeff_denominator(hm1, h0, hp1, hp2)
-      divisor = coeff_divisor(hm1, h0, hp1)
+      denominator = coeff_denominator(hm1, h0, hp1)
+      divisor = coeff_divisor(hm1, h0, hp1, hp2)
     end associate
     
   end subroutine coeff_d_components
@@ -129,7 +129,7 @@ pure real function coeff_numerator(h0, hp1, hp2, alpha, beta)
     
   end function coeff_numerator
   
-  pure real function coeff_denominator(hm1, h0, hp1, hp2)
+  pure real function coeff_denominator(hm1, h0, hp1)
     !! 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
@@ -138,14 +138,13 @@ pure real function coeff_denominator(hm1, h0, hp1, hp2)
     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)
+    coeff_denominator = hm1 * (hm1 + h0) &
+         * (hm1 + h0 + hp1)
     
   end function coeff_denominator
 
-  pure real function coeff_divisor(hm1, h0, hp1)
+  pure real function coeff_divisor(hm1, 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,
@@ -154,8 +153,9 @@ pure real function coeff_divisor(hm1, h0, hp1)
     real, intent(in) :: hm1
     real, intent(in) :: h0
     real, intent(in) :: hp1
+    real, intent(in) :: hp2
 
-    coeff_divisor = 4.0 * ((hm1 + h0 + hp1) / 3.0) 
+    coeff_divisor = (hm1 + h0 + hp1 + hp2)
   end function coeff_divisor
 
 end submodule nucfd_coeffs_d