From 3a8cfe9c6bc8b0b3ddfdc7954a51e484142c3e4d Mon Sep 17 00:00:00 2001
From: Damian Rouson <rouson@lbl.gov>
Date: Fri, 3 Oct 2025 17:06:03 -0400
Subject: [PATCH 01/40] refac: rename face_centered_vals -> cell_centered

---
 ...alues_m.f90 => cell_centers_extended_m.f90} | 18 +++++++++---------
 src/fortran/cell_centers_extended_s.F90        | 15 +++++++++++++++
 src/fortran/face_values_s.F90                  | 15 ---------------
 src/fortran/gradient_operator_s.F90            | 16 +++++++---------
 src/fortran/gradient_s.F90                     |  2 +-
 src/fortran/mimetic_matrix_s.F90               |  2 +-
 src/fortran/mole_m.f90                         |  2 +-
 7 files changed, 34 insertions(+), 36 deletions(-)
 rename src/fortran/{face_values_m.f90 => cell_centers_extended_m.f90} (87%)
 create mode 100644 src/fortran/cell_centers_extended_s.F90
 delete mode 100644 src/fortran/face_values_s.F90

diff --git a/src/fortran/face_values_m.f90 b/src/fortran/cell_centers_extended_m.f90
similarity index 87%
rename from src/fortran/face_values_m.f90
rename to src/fortran/cell_centers_extended_m.f90
index 22999f89..0ca4e97c 100644
--- a/src/fortran/face_values_m.f90
+++ b/src/fortran/cell_centers_extended_m.f90
@@ -1,9 +1,9 @@
-module face_values_m
+module cell_centers_extended_m
   !! Define an abstraction for face-centered values with a corresonding mimetic gradient operator
   implicit none
 
   private
-  public :: face_values_t
+  public :: cell_centers_extended_t
   public :: gradient_t
   public :: gradient_operator_t
 
@@ -42,7 +42,7 @@ pure module function values(self) result(gradients)
 
   end interface
 
-  type face_values_t
+  type cell_centers_extended_t
     !! Face-centered values 
     private
     double precision, allocatable :: f_(:)
@@ -52,15 +52,15 @@ pure module function values(self) result(gradients)
     procedure, non_overridable, private :: grad
   end type
 
-  interface face_values_t
+  interface cell_centers_extended_t
 
-    pure module function construct(f, k, dx) result(face_values)
+    pure module function construct(f, k, dx) result(cell_centers_extended)
       !! Result is a collection of face-centered values with a mimetic gradient operator
       implicit none
       double precision, intent(in) :: f(:) !! face-centered values
       double precision, intent(in) :: dx !! face spacing (cell width)
       integer, intent(in) :: k !! order of accuracy
-      type(face_values_t) face_values
+      type(cell_centers_extended_t) cell_centers_extended
     end function
 
   end interface
@@ -70,7 +70,7 @@ pure module function construct(f, k, dx) result(face_values)
     pure module function grad(self) result(grad_f)
       !! Result is mimetic gradient of f
       implicit none
-      class(face_values_t), intent(in) :: self
+      class(cell_centers_extended_t), intent(in) :: self
       type(gradient_t) grad_f !! discrete gradient approximation
     end function
 
@@ -105,10 +105,10 @@ pure module function matvec(self, vector) result(gradient)
       !! Apply a matrix operator to a vector
       implicit none
       class(mimetic_matrix_t), intent(in) :: self
-      type(face_values_t), intent(in) :: vector
+      type(cell_centers_extended_t), intent(in) :: vector
       type(gradient_t) gradient
     end function
 
   end interface
 
-end module face_values_m
\ No newline at end of file
+end module cell_centers_extended_m
\ No newline at end of file
diff --git a/src/fortran/cell_centers_extended_s.F90 b/src/fortran/cell_centers_extended_s.F90
new file mode 100644
index 00000000..b679bf87
--- /dev/null
+++ b/src/fortran/cell_centers_extended_s.F90
@@ -0,0 +1,15 @@
+submodule(cell_centers_extended_m) cell_centers_extended_s
+  implicit none
+
+contains
+
+  module procedure construct
+    cell_centers_extended%f_ = f
+    cell_centers_extended%gradient_operator_ = gradient_operator_t(k, dx)
+  end procedure
+
+  module procedure grad
+    grad_f = self%gradient_operator_%mimetic_matrix_ .x. self
+  end procedure
+
+end submodule cell_centers_extended_s
\ No newline at end of file
diff --git a/src/fortran/face_values_s.F90 b/src/fortran/face_values_s.F90
deleted file mode 100644
index cc4a0a59..00000000
--- a/src/fortran/face_values_s.F90
+++ /dev/null
@@ -1,15 +0,0 @@
-submodule(face_values_m) face_values_s
-  implicit none
-
-contains
-
-  module procedure construct
-    face_values%f_ = f
-    face_values%gradient_operator_ = gradient_operator_t(k, dx)
-  end procedure
-
-  module procedure grad
-    grad_f = self%gradient_operator_%mimetic_matrix_ .x. self
-  end procedure
-
-end submodule face_values_s
\ No newline at end of file
diff --git a/src/fortran/gradient_operator_s.F90 b/src/fortran/gradient_operator_s.F90
index 36053a35..bc3c74ce 100644
--- a/src/fortran/gradient_operator_s.F90
+++ b/src/fortran/gradient_operator_s.F90
@@ -1,7 +1,7 @@
 #include "julienne-assert-macros.h"
 #include "mole-language-support.F90"
 
-submodule(face_values_m) gradient_operator_s
+submodule(cell_centers_extended_m) gradient_operator_s
   use julienne_m, only : call_julienne_assert_, string_t
 #if ASSERTIONS
   use julienne_m, only : operator(.isAtLeast.)
@@ -66,6 +66,7 @@ pure function corbino_castillo_M(k, dx) result(row)
 
 #if HAVE_DO_CONCURRENT_TYPE_SPEC_SUPPORT
 
+
   pure function corbino_castillo_Ap(k, dx) result(rows)
     integer, intent(in) :: k
     double precision, intent(in) :: dx
@@ -79,23 +80,20 @@ pure function corbino_castillo_Ap(k, dx) result(rows)
       end do reverse_and_flip_sign
     end associate
   end function
-
 #else
 
   pure function corbino_castillo_Ap(k, dx) result(rows)
     integer, intent(in) :: k
     double precision, intent(in) :: dx
     double precision, allocatable :: rows(:,:)
+    integer row
 
     associate(A => corbino_castillo_A(k, dx))
       allocate(rows , mold=A)
-      block
-        integer row
-        reverse_and_flip_sign: &
-        do concurrent(row = 1:size(rows,1)) default(none) shared(rows, A) 
-          rows(row,:) = -A(row,size(A,2):1)
-        end do reverse_and_flip_sign
-      end block
+      reverse_and_flip_sign: &
+      do concurrent(row = 1:size(rows,1)) default(none) shared(rows, A) 
+        rows(row,:) = -A(row,size(A,2):1:-1)
+      end do reverse_and_flip_sign
     end associate
   end function
 
diff --git a/src/fortran/gradient_s.F90 b/src/fortran/gradient_s.F90
index fc289f04..93c7b410 100644
--- a/src/fortran/gradient_s.F90
+++ b/src/fortran/gradient_s.F90
@@ -1,4 +1,4 @@
-submodule(face_values_m) gradient_s
+submodule(cell_centers_extended_m) gradient_s
   implicit none
 
 contains
diff --git a/src/fortran/mimetic_matrix_s.F90 b/src/fortran/mimetic_matrix_s.F90
index ebe8165e..cdbb8842 100644
--- a/src/fortran/mimetic_matrix_s.F90
+++ b/src/fortran/mimetic_matrix_s.F90
@@ -1,7 +1,7 @@
 #include "mole-language-support.F90"
 #include "julienne-assert-macros.h"
 
-submodule(face_values_m) mimetic_matrix_s
+submodule(cell_centers_extended_m) mimetic_matrix_s
   use julienne_m, only : call_julienne_assert_
   implicit none
 
diff --git a/src/fortran/mole_m.f90 b/src/fortran/mole_m.f90
index e81b61c1..3c49bbad 100644
--- a/src/fortran/mole_m.f90
+++ b/src/fortran/mole_m.f90
@@ -1,5 +1,5 @@
 module mole_m 
   !! MOLE Fortran public entities
-  use face_values_m, only : face_values_t, gradient_t
+  use cell_centers_extended_m, only : cell_centers_extended_t, gradient_t
   implicit none
 end module mole_m

From c8f1cd88998056b92d26ef69d5dcc532ef84117d Mon Sep 17 00:00:00 2001
From: Damian Rouson <rouson@lbl.gov>
Date: Fri, 3 Oct 2025 17:44:17 -0400
Subject: [PATCH 02/40] fix(gfortran): work around gfortran issues

This commit
1. Synchronizes a conditionally compiled
   gfortran workaround to match the updated alternative block.
2. Updates the unit test to work around a gfortran issue with
   `associate`.
---
 src/fortran/cell_centers_extended_m.f90 |  9 +++++----
 src/fortran/gradient_operator_s.F90     |  4 ++--
 src/fortran/mimetic_matrix_s.F90        | 27 +++++++++++++++----------
 test/gradient_operator_test_m.F90       | 23 +++++++--------------
 4 files changed, 30 insertions(+), 33 deletions(-)

diff --git a/src/fortran/cell_centers_extended_m.f90 b/src/fortran/cell_centers_extended_m.f90
index 0ca4e97c..949eca40 100644
--- a/src/fortran/cell_centers_extended_m.f90
+++ b/src/fortran/cell_centers_extended_m.f90
@@ -1,5 +1,6 @@
 module cell_centers_extended_m
-  !! Define an abstraction for face-centered values with a corresonding mimetic gradient operator
+  !! Define an abstraction for the collection of points used to compute gradidents:
+  !! cell centers plus oundaries.
   implicit none
 
   private
@@ -25,7 +26,7 @@ module cell_centers_extended_m
   end type
 
   type gradient_t
-    !! Encapsulate gradient values produced only by .grad. (private data, no constructors)
+    !! Encapsulate gradient values produced only by .grad. (no other constructors)
     private
     double precision, allocatable :: g_(:)
   contains
@@ -43,7 +44,7 @@ pure module function values(self) result(gradients)
   end interface
 
   type cell_centers_extended_t
-    !! Face-centered values 
+    !! Encapsulate information at cell centers and boundaries
     private
     double precision, allocatable :: f_(:)
     type(gradient_operator_t) gradient_operator_
@@ -55,7 +56,7 @@ pure module function values(self) result(gradients)
   interface cell_centers_extended_t
 
     pure module function construct(f, k, dx) result(cell_centers_extended)
-      !! Result is a collection of face-centered values with a mimetic gradient operator
+      !! Result is a collection of cell-centered-extended values with a mimetic gradient operator
       implicit none
       double precision, intent(in) :: f(:) !! face-centered values
       double precision, intent(in) :: dx !! face spacing (cell width)
diff --git a/src/fortran/gradient_operator_s.F90 b/src/fortran/gradient_operator_s.F90
index bc3c74ce..e02f2eaf 100644
--- a/src/fortran/gradient_operator_s.F90
+++ b/src/fortran/gradient_operator_s.F90
@@ -75,7 +75,7 @@ pure function corbino_castillo_Ap(k, dx) result(rows)
     associate(A => corbino_castillo_A(k, dx))
       allocate(rows , mold=A)
       reverse_and_flip_sign: &
-      do concurrent(integer :: row = 1:size(rows,1)) default(none) shared(rows, A) 
+      do concurrent(integer :: row = 1:size(rows,1)) default(none) shared(rows, A)
         rows(row,:) = -A(row,size(A,2):1:-1)
       end do reverse_and_flip_sign
     end associate
@@ -91,7 +91,7 @@ pure function corbino_castillo_Ap(k, dx) result(rows)
     associate(A => corbino_castillo_A(k, dx))
       allocate(rows , mold=A)
       reverse_and_flip_sign: &
-      do concurrent(row = 1:size(rows,1)) default(none) shared(rows, A) 
+      do concurrent(row = 1:size(rows,1)) default(none) shared(rows, A)
         rows(row,:) = -A(row,size(A,2):1:-1)
       end do reverse_and_flip_sign
     end associate
diff --git a/src/fortran/mimetic_matrix_s.F90 b/src/fortran/mimetic_matrix_s.F90
index cdbb8842..20156484 100644
--- a/src/fortran/mimetic_matrix_s.F90
+++ b/src/fortran/mimetic_matrix_s.F90
@@ -48,23 +48,28 @@
 
     double precision, allocatable :: product_inner(:)
 
-    associate(upper_rows  => size(self%upper_,1), inner_bandwidth => size(self%inner_))
+    associate(upper_rows  => size(self%upper_,1))
       associate(inner_rows => (size(vector%f_) - 1) - 2*upper_rows) ! inner_rows = matrix rows - (upper rows + lower rows)
 
         allocate(product_inner(inner_rows))
 
-        block
+        associate(inner_bandwidth => size(self%inner_))
+          block
           integer row
-          do concurrent(row = 1 : inner_rows) default(none) shared(product_inner, self, vector, upper_rows, inner_bandwidth)
-            product_inner(row) = dot_product(self%inner_, vector%f_(upper_rows + row : upper_rows + inner_bandwidth))
-          end do
-        end block
+            do concurrent(row = 1 : inner_rows) default(none) &
+              shared(product_inner, self, vector, upper_rows, inner_bandwidth)
+              product_inner(row) = dot_product(self%inner_, vector%f_(upper_rows + row : upper_rows + row + inner_bandwidth - 1))
+            end do
+          end block
+        end associate
 
-        gradient%g_ = [ &
-           matmul(self%upper_, vector%f_(1 : upper_rows)) &
-          ,product_inner &
-          ,matmul(self%lower_, vector%f_(upper_rows + inner_rows + 1 : )) &
-        ]
+        associate(upper_bandwidth => size(self%upper_,2), lower_bandwidth => size(self%lower_,2))
+          gradient%g_ = [ &
+             matmul(self%upper_, vector%f_(1 : upper_bandwidth)) &
+            ,product_inner &
+            ,matmul(self%lower_, vector%f_(size(vector%f_) - lower_bandwidth + 1 : )) &
+          ]
+        end associate
       end associate
     end associate
   end procedure
diff --git a/test/gradient_operator_test_m.F90 b/test/gradient_operator_test_m.F90
index 2df4dd6a..4dd7f4f1 100644
--- a/test/gradient_operator_test_m.F90
+++ b/test/gradient_operator_test_m.F90
@@ -7,7 +7,7 @@
 module gradient_operator_test_m
   use julienne_m, only : &
      test_t, test_description_t, test_diagnosis_t, test_result_t, operator(.all.), operator(.approximates.), operator(.within.)
-  use mole_m, only : face_values_t
+  use mole_m, only : cell_centers_extended_t, gradient_t
 #if ! HAVE_PROCEDURE_ACTUAL_FOR_POINTER_DUMMY
   use julienne_m, only : diagnosis_function_i
 #endif
@@ -52,21 +52,12 @@ function results() result(test_results)
 
   function check_gradient_of_line() result(test_diagnosis)
     type(test_diagnosis_t) test_diagnosis
-    double precision, parameter :: dx = 1D0, df_dx_exact = 1D0
-    double precision, parameter :: x(*) = [0D0,.5D0, 1.5D0, 2.5D0, 3.5D0, 4D0]*dx
-    double precision, parameter :: tolerance = 1D-15
-
-    associate(f => face_values_t(linear(x), k=2, dx=dx))
-      associate(grad_f => .grad. f)
-        test_diagnosis = .all. (grad_f%values() .approximates. df_dx_exact .within. tolerance)
-      end associate
-    end associate
-
-  contains
-    double precision elemental function linear(x)
-      double precision, intent(in) :: x
-      linear = x
-    end function
+    double precision, parameter :: dx = 1D0, tolerance = 1D-15
+    double precision, parameter :: x(*) = [0D0,.5D0, 1.5D0, 2.5D0, 3.5D0, 4D0]*dx !! grid from Corbino & Castillo (2020) Fig. 6
+    double precision, parameter :: f(*) = x, df_dx_exact = 1D0 !! f(x) = x, df_dx = 1
+    type(gradient_t) grad_f
+    grad_f = .grad. cell_centers_extended_t(f, k=2, dx=dx) ! gfortran blocks use of association
+    test_diagnosis = .all. (grad_f%values() .approximates. df_dx_exact .within. tolerance)
   end function
 
 end module

From c1ee1219e615bf97a88203d2a2b6b2bb693f00d1 Mon Sep 17 00:00:00 2001
From: Damian Rouson <rouson@lbl.gov>
Date: Sat, 4 Oct 2025 21:12:15 -0400
Subject: [PATCH 03/40] test(gradient): add 2nd line differentiation case

---
 test/gradient_operator_test_m.F90 | 38 +++++++++++++++++++++++--------
 1 file changed, 28 insertions(+), 10 deletions(-)

diff --git a/test/gradient_operator_test_m.F90 b/test/gradient_operator_test_m.F90
index 4dd7f4f1..7a8f4ef2 100644
--- a/test/gradient_operator_test_m.F90
+++ b/test/gradient_operator_test_m.F90
@@ -6,7 +6,15 @@
 
 module gradient_operator_test_m
   use julienne_m, only : &
-     test_t, test_description_t, test_diagnosis_t, test_result_t, operator(.all.), operator(.approximates.), operator(.within.)
+     string_t &
+    ,test_t &
+    ,test_description_t &
+    ,test_diagnosis_t &
+    ,test_result_t &
+    ,operator(//) &
+    ,operator(.all.) &
+    ,operator(.approximates.) &
+    ,operator(.within.)
   use mole_m, only : cell_centers_extended_t, gradient_t
 #if ! HAVE_PROCEDURE_ACTUAL_FOR_POINTER_DUMMY
   use julienne_m, only : diagnosis_function_i
@@ -19,6 +27,8 @@ module gradient_operator_test_m
     procedure, nopass :: results
   end type
 
+  double precision, parameter :: line_tolerance = 1D-14
+
 contains
 
   pure function subject() result(test_subject)
@@ -32,7 +42,7 @@ function results() result(test_results)
     type(gradient_operator_test_t) gradient_operator_test
     type(test_result_t), allocatable :: test_results(:)
     test_results = gradient_operator_test%run([ & 
-       test_description_t('computing the gradient of a linear function', check_gradient_of_line) &
+       test_description_t('computing gradients of lines within a tolerance of ' // string_t(line_tolerance), check_line_slope) &
     ])
   end function
 
@@ -42,22 +52,30 @@ function results() result(test_results)
     type(gradient_operator_test_t) gradient_operator_test
     type(test_result_t), allocatable :: test_results(:)
     procedure(diagnosis_function_i), pointer :: &
-      check_gradient_of_line_ptr => check_gradient_of_line
+      check_line_slope_ptr => check_line_slope
     test_results = gradient_operator_test%run([ &
-       test_description_t('computing the gradient of a linear function', check_gradient_of_line_ptr) &
+       test_description_t('computing gradients of lines within a tolerance of ' // string_t(line_tolerance), check_line_slope_ptr) &
     ])
   end function
 
 #endif
 
-  function check_gradient_of_line() result(test_diagnosis)
+  function check_line_slope() result(test_diagnosis)
     type(test_diagnosis_t) test_diagnosis
-    double precision, parameter :: dx = 1D0, tolerance = 1D-15
-    double precision, parameter :: x(*) = [0D0,.5D0, 1.5D0, 2.5D0, 3.5D0, 4D0]*dx !! grid from Corbino & Castillo (2020) Fig. 6
-    double precision, parameter :: f(*) = x, df_dx_exact = 1D0 !! f(x) = x, df_dx = 1
-    type(gradient_t) grad_f
+    integer i
+    integer, parameter :: cells = 10
+    double precision, parameter :: x_min = 0., x_max = 10., dx = (x_max - x_min)/dble(cells)
+    double precision, parameter :: x(*) = [x_min, x_min + dx/2. + [(dble(i-1)*dx, i = 1, cells)], x_max]*dx
+      !! boundaries + grid cell centers -- see Corbino & Castillo (2020) https://doi.org/10.1016/j.cam.2019.06.042
+    double precision, parameter :: m = 2D0, b = 3D0, n = 5D0, c = 7D0
+    double precision, parameter :: f(*) = m*x + b, df_dx = m
+    double precision, parameter :: g(*) = n*x + c, dg_dx = n
+    type(gradient_t) grad_f, grad_g
     grad_f = .grad. cell_centers_extended_t(f, k=2, dx=dx) ! gfortran blocks use of association
-    test_diagnosis = .all. (grad_f%values() .approximates. df_dx_exact .within. tolerance)
+    test_diagnosis = .all. (grad_f%values() .approximates. df_dx .within. line_tolerance) // " (df_dx)"
+
+    grad_g = .grad. cell_centers_extended_t(g, k=2, dx=dx) ! gfortran blocks use of association
+    test_diagnosis = .all. (grad_g%values() .approximates. dg_dx .within. line_tolerance) // " (dg_dx)"
   end function
 
 end module

From 14ef8f5fa2a9587423a4b7c9d75a275f5a39b261 Mon Sep 17 00:00:00 2001
From: Damian Rouson <rouson@lbl.gov>
Date: Tue, 7 Oct 2025 08:45:29 -0400
Subject: [PATCH 04/40] fix(matvec): use size() to set loop limits

This commit fixes issues with the matrix-vector multiplication
function "matvec" to ensure infers the loop limits and array
section bounds from the dimensions of the arrays accessed by the
loop or array section.
---
 src/fortran/cell_centers_extended_m.f90 | 21 +++++----
 src/fortran/cell_centers_extended_s.F90 | 20 ++++++++-
 src/fortran/gradient_operator_s.F90     |  2 +
 src/fortran/initializers_m.f90          | 28 ++++++++++++
 src/fortran/mimetic_matrix_s.F90        | 59 +++++--------------------
 src/fortran/mole_m.f90                  |  1 +
 test/gradient_operator_test_m.F90       | 35 ++++++++-------
 7 files changed, 92 insertions(+), 74 deletions(-)
 create mode 100644 src/fortran/initializers_m.f90

diff --git a/src/fortran/cell_centers_extended_m.f90 b/src/fortran/cell_centers_extended_m.f90
index 949eca40..e1021568 100644
--- a/src/fortran/cell_centers_extended_m.f90
+++ b/src/fortran/cell_centers_extended_m.f90
@@ -1,6 +1,7 @@
 module cell_centers_extended_m
   !! Define an abstraction for the collection of points used to compute gradidents:
   !! cell centers plus oundaries.
+  use initializers_m, only : scalar_1D_initializer_t
   implicit none
 
   private
@@ -20,7 +21,7 @@ module cell_centers_extended_m
   type gradient_operator_t
     !! Encapsulate kth-order mimetic gradient operator on dx-sized cells
     private
-    integer k_
+    integer k_, m_
     double precision dx_
     type(mimetic_matrix_t) mimetic_matrix_
   end type
@@ -46,7 +47,7 @@ pure module function values(self) result(gradients)
   type cell_centers_extended_t
     !! Encapsulate information at cell centers and boundaries
     private
-    double precision, allocatable :: f_(:)
+    double precision, allocatable :: scalar_1D_(:), domain_(:)
     type(gradient_operator_t) gradient_operator_
   contains
     generic :: operator(.grad.) => grad
@@ -55,12 +56,13 @@ pure module function values(self) result(gradients)
 
   interface cell_centers_extended_t
 
-    pure module function construct(f, k, dx) result(cell_centers_extended)
-      !! Result is a collection of cell-centered-extended values with a mimetic gradient operator
+    pure module function construct(scalar_1D_initializer, order, cells, domain) result(cell_centers_extended)
+      !! Result is a collection of cell-centered-extended values with a corresponding mimetic gradient operator
       implicit none
-      double precision, intent(in) :: f(:) !! face-centered values
-      double precision, intent(in) :: dx !! face spacing (cell width)
-      integer, intent(in) :: k !! order of accuracy
+      class(scalar_1D_initializer_t), intent(in) :: scalar_1D_initializer !! elemental initialization function hook
+      integer, intent(in) :: order !! order of accuracy
+      integer, intent(in) :: cells !! number of grid cells spanning domain
+      double precision, intent(in) :: domain(:) !! [grid minimum, grid maximum]
       type(cell_centers_extended_t) cell_centers_extended
     end function
 
@@ -79,11 +81,12 @@ pure module function grad(self) result(grad_f)
 
   interface gradient_operator_t
 
-    pure module function construct_from_parameters(k, dx) result(gradient_operator)
+    pure module function construct_from_parameters(k, dx, m) result(gradient_operator)
       !! Construct a mimetic gradient operator
       implicit none
       integer, intent(in) :: k !! order of accuracy
-      double precision, intent(in) :: dx !! step size
+      double precision, intent(in) :: dx !! step siz
+      integer, intent(in) :: m !! number of grid cells
       type(gradient_operator_t) gradient_operator
     end function
 
diff --git a/src/fortran/cell_centers_extended_s.F90 b/src/fortran/cell_centers_extended_s.F90
index b679bf87..58c0bd9b 100644
--- a/src/fortran/cell_centers_extended_s.F90
+++ b/src/fortran/cell_centers_extended_s.F90
@@ -1,11 +1,27 @@
+#include "julienne-assert-macros.h"
+
 submodule(cell_centers_extended_m) cell_centers_extended_s
+  use julienne_m, only : call_julienne_assert_, operator(.equalsExpected.)
   implicit none
 
 contains
 
   module procedure construct
-    cell_centers_extended%f_ = f
-    cell_centers_extended%gradient_operator_ = gradient_operator_t(k, dx)
+
+    integer cell
+
+    call_julienne_assert(size(domain) .equalsExpected. 2)
+
+    associate(x_min => domain(1), x_max => domain(2))
+      associate(dx => dble(domain(2) - domain(1))/dble(cells))
+        associate(x => [x_min, x_min + dx/2. + [(dble(cell-1)*dx, cell = 1, cells)], x_max]*dx) !! boundaries + cell centers
+          cell_centers_extended%scalar_1D_ = scalar_1D_initializer%f(x)                         !! Corbino & Castillo (2020)
+        end associate                                                                           !! https://doi.org/10.1016/j.cam.2019.06.042
+      end associate
+    end associate
+
+    cell_centers_extended%domain_ = domain
+    cell_centers_extended%gradient_operator_ = gradient_operator_t(k=order, dx=(domain(2)-domain(1))/dble(cells), m=cells)
   end procedure
 
   module procedure grad
diff --git a/src/fortran/gradient_operator_s.F90 b/src/fortran/gradient_operator_s.F90
index e02f2eaf..8514e15e 100644
--- a/src/fortran/gradient_operator_s.F90
+++ b/src/fortran/gradient_operator_s.F90
@@ -21,6 +21,7 @@
     )
     gradient_operator%k_  = k
     gradient_operator%dx_ = dx
+    gradient_operator%m_ = m
   end procedure
 
   pure function corbino_castillo_A(k, dx) result(rows)
@@ -80,6 +81,7 @@ pure function corbino_castillo_Ap(k, dx) result(rows)
       end do reverse_and_flip_sign
     end associate
   end function
+ 
 #else
 
   pure function corbino_castillo_Ap(k, dx) result(rows)
diff --git a/src/fortran/initializers_m.f90 b/src/fortran/initializers_m.f90
new file mode 100644
index 00000000..88fbe1a1
--- /dev/null
+++ b/src/fortran/initializers_m.f90
@@ -0,0 +1,28 @@
+module initializers_m
+  !! Implement a workaround for the Fortran standard's prohibition against
+  !! elemental procedures as dummy arguments. Users can extend the abstract
+  !! type(s) in this module and define the deferred binding as an elemental
+  !! function for use in initializing variables at grid locations without
+  !! requiring loops.
+  implicit none
+
+  private
+  public :: scalar_1D_initializer_t
+
+  abstract interface
+
+    elemental function scalar_1D_initializer_i(x) result(f)
+      implicit none
+      double precision, intent(in) :: x 
+      double precision f
+    end function
+
+  end interface
+
+  type, abstract :: scalar_1D_initializer_t
+    !! Define a hook on which to hang elemental grid-variable initializers
+  contains
+    procedure(scalar_1D_initializer_i), deferred, nopass :: f
+  end type
+
+end module initializers_m
diff --git a/src/fortran/mimetic_matrix_s.F90 b/src/fortran/mimetic_matrix_s.F90
index 20156484..cb89443a 100644
--- a/src/fortran/mimetic_matrix_s.F90
+++ b/src/fortran/mimetic_matrix_s.F90
@@ -2,7 +2,7 @@
 #include "julienne-assert-macros.h"
 
 submodule(cell_centers_extended_m) mimetic_matrix_s
-  use julienne_m, only : call_julienne_assert_
+  use julienne_m, only : call_julienne_assert_, string_t
   implicit none
 
 contains
@@ -19,61 +19,26 @@
 
     double precision, allocatable :: product_inner(:)
 
-    associate(upper_rows  => size(self%upper_,1))
-      associate(inner_rows => (size(vector%f_) - 1) - 2*upper_rows) ! inner_rows = matrix rows - (upper rows + lower rows)
+    associate(upper => merge(0, 1, size(self%upper_)==0), lower => merge(0, 1, size(self%lower_)==0))
+      associate(inner_rows => size(vector%scalar_1D_) - (upper + lower + 1))
 
         allocate(product_inner(inner_rows))
 
-        associate(inner_bandwidth => size(self%inner_))
-          do concurrent(integer :: row = 1 : inner_rows) default(none) &
-            shared(product_inner, self, vector, upper_rows, inner_bandwidth)
-            product_inner(row) = dot_product(self%inner_, vector%f_(upper_rows + row : upper_rows + row + inner_bandwidth - 1))
-          end do
-        end associate
+        do concurrent(integer :: row = 1 : size(product_inner)) default(none) shared(product_inner, self, vector)
+          product_inner(row) = dot_product(self%inner_, vector%scalar_1D_(row + 1 : row + size(self%inner_)))
+        end do
 
-        associate(upper_bandwidth => size(self%upper_,2), lower_bandwidth => size(self%lower_,2))
-          gradient%g_ = [ &
-             matmul(self%upper_, vector%f_(1 : upper_bandwidth)) &
-            ,product_inner &
-            ,matmul(self%lower_, vector%f_(size(vector%f_) - lower_bandwidth + 1 : )) &
-          ]
-        end associate
+        gradient%g_ = [ &
+           matmul(self%upper_, vector%scalar_1D_(1 : size(self%upper_,2))) &
+          ,product_inner &
+          ,matmul(self%lower_, vector%scalar_1D_(size(vector%scalar_1D_) - size(self%lower_,2) + 1 : )) &
+        ]
       end associate
     end associate
   end procedure
 
 #else
 
-  module procedure matvec
-
-    double precision, allocatable :: product_inner(:)
-
-    associate(upper_rows  => size(self%upper_,1))
-      associate(inner_rows => (size(vector%f_) - 1) - 2*upper_rows) ! inner_rows = matrix rows - (upper rows + lower rows)
-
-        allocate(product_inner(inner_rows))
-
-        associate(inner_bandwidth => size(self%inner_))
-          block
-          integer row
-            do concurrent(row = 1 : inner_rows) default(none) &
-              shared(product_inner, self, vector, upper_rows, inner_bandwidth)
-              product_inner(row) = dot_product(self%inner_, vector%f_(upper_rows + row : upper_rows + row + inner_bandwidth - 1))
-            end do
-          end block
-        end associate
-
-        associate(upper_bandwidth => size(self%upper_,2), lower_bandwidth => size(self%lower_,2))
-          gradient%g_ = [ &
-             matmul(self%upper_, vector%f_(1 : upper_bandwidth)) &
-            ,product_inner &
-            ,matmul(self%lower_, vector%f_(size(vector%f_) - lower_bandwidth + 1 : )) &
-          ]
-        end associate
-      end associate
-    end associate
-  end procedure
-
 #endif
 
-end submodule mimetic_matrix_s
+end submodule mimetic_matrix_s
\ No newline at end of file
diff --git a/src/fortran/mole_m.f90 b/src/fortran/mole_m.f90
index 3c49bbad..a97aaee9 100644
--- a/src/fortran/mole_m.f90
+++ b/src/fortran/mole_m.f90
@@ -1,5 +1,6 @@
 module mole_m 
   !! MOLE Fortran public entities
   use cell_centers_extended_m, only : cell_centers_extended_t, gradient_t
+  use initializers_m, only : scalar_1D_initializer_t
   implicit none
 end module mole_m
diff --git a/test/gradient_operator_test_m.F90 b/test/gradient_operator_test_m.F90
index 7a8f4ef2..51b399de 100644
--- a/test/gradient_operator_test_m.F90
+++ b/test/gradient_operator_test_m.F90
@@ -15,7 +15,7 @@ module gradient_operator_test_m
     ,operator(.all.) &
     ,operator(.approximates.) &
     ,operator(.within.)
-  use mole_m, only : cell_centers_extended_t, gradient_t
+  use mole_m, only : cell_centers_extended_t, gradient_t, scalar_1D_initializer_t
 #if ! HAVE_PROCEDURE_ACTUAL_FOR_POINTER_DUMMY
   use julienne_m, only : diagnosis_function_i
 #endif
@@ -62,20 +62,23 @@ function results() result(test_results)
 
   function check_line_slope() result(test_diagnosis)
     type(test_diagnosis_t) test_diagnosis
-    integer i
-    integer, parameter :: cells = 10
-    double precision, parameter :: x_min = 0., x_max = 10., dx = (x_max - x_min)/dble(cells)
-    double precision, parameter :: x(*) = [x_min, x_min + dx/2. + [(dble(i-1)*dx, i = 1, cells)], x_max]*dx
-      !! boundaries + grid cell centers -- see Corbino & Castillo (2020) https://doi.org/10.1016/j.cam.2019.06.042
-    double precision, parameter :: m = 2D0, b = 3D0, n = 5D0, c = 7D0
-    double precision, parameter :: f(*) = m*x + b, df_dx = m
-    double precision, parameter :: g(*) = n*x + c, dg_dx = n
-    type(gradient_t) grad_f, grad_g
-    grad_f = .grad. cell_centers_extended_t(f, k=2, dx=dx) ! gfortran blocks use of association
-    test_diagnosis = .all. (grad_f%values() .approximates. df_dx .within. line_tolerance) // " (df_dx)"
-
-    grad_g = .grad. cell_centers_extended_t(g, k=2, dx=dx) ! gfortran blocks use of association
-    test_diagnosis = .all. (grad_g%values() .approximates. dg_dx .within. line_tolerance) // " (dg_dx)"
+    type(gradient_t) grad_f
+    type, extends(scalar_1D_initializer_t) :: const_initializer_1D_t
+    contains
+      procedure, nopass, non_overridable :: f => const
+    end type
+    type(const_initializer_1D_t) :: const_initializer_1D
+
+    double precision, parameter :: df_dx = 0D0
+
+    grad_f = .grad. cell_centers_extended_t(const_initializer_1D, order=2, cells=4, domain=[0D0,1D0]) ! gfortran blocks use of association
+    test_diagnosis = .all. (grad_f%values() .approximates. df_dx .within. line_tolerance) // " (d(const)/dx)"
+  end function
+
+  elemental function const(x) result(c)
+    double precision, intent(in) :: x
+    double precision c
+    c = 2D0
   end function
 
-end module
+end module
\ No newline at end of file

From f602339dce2f0533f492d1876389cb7702a2ef4f Mon Sep 17 00:00:00 2001
From: Damian Rouson <rouson@lbl.gov>
Date: Tue, 7 Oct 2025 14:31:16 -0400
Subject: [PATCH 05/40] refac(test): mv code, rename function

---
 test/gradient_operator_test_m.F90 | 20 ++++++++++----------
 1 file changed, 10 insertions(+), 10 deletions(-)

diff --git a/test/gradient_operator_test_m.F90 b/test/gradient_operator_test_m.F90
index 51b399de..fd8f7a18 100644
--- a/test/gradient_operator_test_m.F90
+++ b/test/gradient_operator_test_m.F90
@@ -42,7 +42,7 @@ function results() result(test_results)
     type(gradient_operator_test_t) gradient_operator_test
     type(test_result_t), allocatable :: test_results(:)
     test_results = gradient_operator_test%run([ & 
-       test_description_t('computing gradients of lines within a tolerance of ' // string_t(line_tolerance), check_line_slope) &
+       test_description_t('computing gradients of lines within a tolerance of ' // string_t(line_tolerance), check_grad_const) &
     ])
   end function
 
@@ -52,15 +52,21 @@ function results() result(test_results)
     type(gradient_operator_test_t) gradient_operator_test
     type(test_result_t), allocatable :: test_results(:)
     procedure(diagnosis_function_i), pointer :: &
-      check_line_slope_ptr => check_line_slope
+      check_grad_const_ptr => check_grad_const
     test_results = gradient_operator_test%run([ &
-       test_description_t('computing gradients of lines within a tolerance of ' // string_t(line_tolerance), check_line_slope_ptr) &
+       test_description_t('computing gradients of lines within a tolerance of ' // string_t(line_tolerance), check_grad_const_ptr) &
     ])
   end function
 
 #endif
 
-  function check_line_slope() result(test_diagnosis)
+  elemental function const(x) result(c)
+    double precision, intent(in) :: x
+    double precision c
+    c = 2D0
+  end function
+
+  function check_grad_const() result(test_diagnosis)
     type(test_diagnosis_t) test_diagnosis
     type(gradient_t) grad_f
     type, extends(scalar_1D_initializer_t) :: const_initializer_1D_t
@@ -75,10 +81,4 @@ function check_line_slope() result(test_diagnosis)
     test_diagnosis = .all. (grad_f%values() .approximates. df_dx .within. line_tolerance) // " (d(const)/dx)"
   end function
 
-  elemental function const(x) result(c)
-    double precision, intent(in) :: x
-    double precision c
-    c = 2D0
-  end function
-
 end module
\ No newline at end of file

From 2ed78fc0ef10f004673abe508497a5454bdb1973 Mon Sep 17 00:00:00 2001
From: Damian Rouson <rouson@lbl.gov>
Date: Tue, 7 Oct 2025 14:44:17 -0400
Subject: [PATCH 06/40] fix(grid): rm extra multiplicative factor of dx

---
 src/fortran/cell_centers_extended_s.F90 | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/fortran/cell_centers_extended_s.F90 b/src/fortran/cell_centers_extended_s.F90
index 58c0bd9b..02115528 100644
--- a/src/fortran/cell_centers_extended_s.F90
+++ b/src/fortran/cell_centers_extended_s.F90
@@ -14,7 +14,7 @@
 
     associate(x_min => domain(1), x_max => domain(2))
       associate(dx => dble(domain(2) - domain(1))/dble(cells))
-        associate(x => [x_min, x_min + dx/2. + [(dble(cell-1)*dx, cell = 1, cells)], x_max]*dx) !! boundaries + cell centers
+        associate(x => [x_min, x_min + dx/2. + [(dble(cell-1)*dx, cell = 1, cells)], x_max]) !! boundaries + cell centers
           cell_centers_extended%scalar_1D_ = scalar_1D_initializer%f(x)                         !! Corbino & Castillo (2020)
         end associate                                                                           !! https://doi.org/10.1016/j.cam.2019.06.042
       end associate

From fdca554b1547e3489ab33c73dd27c0eb8e12b881 Mon Sep 17 00:00:00 2001
From: Damian Rouson <rouson@lbl.gov>
Date: Tue, 7 Oct 2025 15:57:41 -0400
Subject: [PATCH 07/40] test(grad): add passing test differentiataing line

---
 src/fortran/cell_centers_extended_m.f90 |  4 +-
 src/fortran/cell_centers_extended_s.F90 | 12 ++---
 test/gradient_operator_test_m.F90       | 62 ++++++++++++++++++++++---
 3 files changed, 61 insertions(+), 17 deletions(-)

diff --git a/src/fortran/cell_centers_extended_m.f90 b/src/fortran/cell_centers_extended_m.f90
index e1021568..33e5cc08 100644
--- a/src/fortran/cell_centers_extended_m.f90
+++ b/src/fortran/cell_centers_extended_m.f90
@@ -46,8 +46,8 @@ pure module function values(self) result(gradients)
 
   type cell_centers_extended_t
     !! Encapsulate information at cell centers and boundaries
-    private
-    double precision, allocatable :: scalar_1D_(:), domain_(:)
+    !private
+    double precision, allocatable :: scalar_1D_(:), domain_(:), grid_(:)
     type(gradient_operator_t) gradient_operator_
   contains
     generic :: operator(.grad.) => grad
diff --git a/src/fortran/cell_centers_extended_s.F90 b/src/fortran/cell_centers_extended_s.F90
index 02115528..98ab9eb0 100644
--- a/src/fortran/cell_centers_extended_s.F90
+++ b/src/fortran/cell_centers_extended_s.F90
@@ -12,16 +12,12 @@
 
     call_julienne_assert(size(domain) .equalsExpected. 2)
 
-    associate(x_min => domain(1), x_max => domain(2))
-      associate(dx => dble(domain(2) - domain(1))/dble(cells))
-        associate(x => [x_min, x_min + dx/2. + [(dble(cell-1)*dx, cell = 1, cells)], x_max]) !! boundaries + cell centers
-          cell_centers_extended%scalar_1D_ = scalar_1D_initializer%f(x)                         !! Corbino & Castillo (2020)
-        end associate                                                                           !! https://doi.org/10.1016/j.cam.2019.06.042
-      end associate
+    associate(x_min => domain(1), x_max => domain(2), dx => dble(domain(2) - domain(1))/dble(cells))
+      cell_centers_extended%grid_ = [x_min, x_min + dx/2. + [((cell-1)*dx, cell = 1, cells)], x_max] ! boundaries + cell centers as described in
+      cell_centers_extended%scalar_1D_ = scalar_1D_initializer%f(cell_centers_extended%grid_)        ! Corbino & Castillo (2020)
+      cell_centers_extended%gradient_operator_ = gradient_operator_t(k=order, dx=dx, m=cells)        ! https://doi.org/10.1016/j.cam.2019.06.042
     end associate
 
-    cell_centers_extended%domain_ = domain
-    cell_centers_extended%gradient_operator_ = gradient_operator_t(k=order, dx=(domain(2)-domain(1))/dble(cells), m=cells)
   end procedure
 
   module procedure grad
diff --git a/test/gradient_operator_test_m.F90 b/test/gradient_operator_test_m.F90
index fd8f7a18..76ae2c32 100644
--- a/test/gradient_operator_test_m.F90
+++ b/test/gradient_operator_test_m.F90
@@ -27,7 +27,7 @@ module gradient_operator_test_m
     procedure, nopass :: results
   end type
 
-  double precision, parameter :: line_tolerance = 1D-14
+  double precision, parameter :: tight_tolerance = 1D-14, loose_tolerance = 1D-10
 
 contains
 
@@ -42,7 +42,8 @@ function results() result(test_results)
     type(gradient_operator_test_t) gradient_operator_test
     type(test_result_t), allocatable :: test_results(:)
     test_results = gradient_operator_test%run([ & 
-       test_description_t('computing gradients of lines within a tolerance of ' // string_t(line_tolerance), check_grad_const) &
+       test_description_t('computing the gradient of a constant function within a tolerance of ' // string_t(tight_tolerance), check_grad_const) &
+      ,test_description_t('computing the gradient of a line within a tolerance of ' // string_t(loose_tolerance), check_grad_line) &
     ])
   end function
 
@@ -52,9 +53,11 @@ function results() result(test_results)
     type(gradient_operator_test_t) gradient_operator_test
     type(test_result_t), allocatable :: test_results(:)
     procedure(diagnosis_function_i), pointer :: &
-      check_grad_const_ptr => check_grad_const
+       check_grad_const_ptr => check_grad_const &
+      ,check_grad_line_ptr => check_grad_line
     test_results = gradient_operator_test%run([ &
-       test_description_t('computing gradients of lines within a tolerance of ' // string_t(line_tolerance), check_grad_const_ptr) &
+       test_description_t('computing the gradient of a constant function within a tolerance of ' // string_t(tight_tolerance), check_grad_const_ptr) &
+      ,test_description_t('computing the gradient of a line within a tolerance of ' // string_t(loose_tolerance), check_grad_line_ptr) &
     ])
   end function
 
@@ -63,7 +66,7 @@ function results() result(test_results)
   elemental function const(x) result(c)
     double precision, intent(in) :: x
     double precision c
-    c = 2D0
+    c = 3D0
   end function
 
   function check_grad_const() result(test_diagnosis)
@@ -74,11 +77,56 @@ function check_grad_const() result(test_diagnosis)
       procedure, nopass, non_overridable :: f => const
     end type
     type(const_initializer_1D_t) :: const_initializer_1D
-
     double precision, parameter :: df_dx = 0D0
 
     grad_f = .grad. cell_centers_extended_t(const_initializer_1D, order=2, cells=4, domain=[0D0,1D0]) ! gfortran blocks use of association
-    test_diagnosis = .all. (grad_f%values() .approximates. df_dx .within. line_tolerance) // " (d(const)/dx)"
+    test_diagnosis = .all. (grad_f%values() .approximates. df_dx .within. tight_tolerance) // " (d(const)/dx)"
+  end function
+
+  elemental function line(x) result(y)
+    double precision, intent(in) :: x
+    double precision y
+    y = 14*x + 3
+  end function
+
+  function check_grad_line() result(test_diagnosis)
+    type(test_diagnosis_t) test_diagnosis
+    type(gradient_t) grad_f
+    type, extends(scalar_1D_initializer_t) :: line_initializer_1D_t
+    contains
+      procedure, nopass, non_overridable :: f => line
+    end type
+    type(line_initializer_1D_t) :: line_initializer_1D
+    double precision, parameter :: df_dx = 14D0
+
+    grad_f = .grad. cell_centers_extended_t(line_initializer_1D, order=2, cells=4, domain=[0D0,1D0]) ! gfortran blocks use of association
+    test_diagnosis = .all. (grad_f%values() .approximates. df_dx .within. loose_tolerance) // " (d(line)/dx)"
+
+  end function
+
+  elemental function parabola(x) result(y)
+    double precision, intent(in) :: x
+    double precision y
+    y = 7*x**2 + 3*x + 5
+  end function
+
+  function check_grad_parabola() result(test_diagnosis)
+    type(test_diagnosis_t) test_diagnosis
+    type(gradient_t) grad_f
+    type(cell_centers_extended_t) quadratic
+    type, extends(scalar_1D_initializer_t) :: parabola_initializer_1D_t
+    contains
+      procedure, nopass, non_overridable :: f => line
+    end type
+    type(parabola_initializer_1D_t) parabola_initializer_1D
+
+    quadratic = cell_centers_extended_t(parabola_initializer_1D, order=2, cells=4, domain=[0D0,1D0]) ! gfortran blocks use of association
+    grad_f = .grad. quadratic ! gfortran blocks use of association
+    print *, "grad_f = ", grad_f%values()
+    print *, "df_dx = ",parabola(quadratic%grid_)
+
+    test_diagnosis = test_diagnosis_t(.true., "") !.all. (grad_f%values() .approximates. df_dx .within. loose_tolerance) // " (d(line)/dx)"
+
   end function
 
 end module
\ No newline at end of file

From 88e3b4bde979db9213b3dc260f6308cfa36d5044 Mon Sep 17 00:00:00 2001
From: Damian Rouson <rouson@lbl.gov>
Date: Tue, 7 Oct 2025 17:58:30 -0400
Subject: [PATCH 08/40] refac(cell_centers_ex): domain(:) -> x_{min,max}

---
 src/fortran/cell_centers_extended_m.f90 | 10 ++++++----
 src/fortran/cell_centers_extended_s.F90 |  4 ++--
 test/gradient_operator_test_m.F90       |  6 +++---
 3 files changed, 11 insertions(+), 9 deletions(-)

diff --git a/src/fortran/cell_centers_extended_m.f90 b/src/fortran/cell_centers_extended_m.f90
index 33e5cc08..25662572 100644
--- a/src/fortran/cell_centers_extended_m.f90
+++ b/src/fortran/cell_centers_extended_m.f90
@@ -47,7 +47,8 @@ pure module function values(self) result(gradients)
   type cell_centers_extended_t
     !! Encapsulate information at cell centers and boundaries
     !private
-    double precision, allocatable :: scalar_1D_(:), domain_(:), grid_(:)
+    double precision, allocatable :: scalar_1D_(:), grid_(:)
+    double precision x_min_, x_max_
     type(gradient_operator_t) gradient_operator_
   contains
     generic :: operator(.grad.) => grad
@@ -56,13 +57,14 @@ pure module function values(self) result(gradients)
 
   interface cell_centers_extended_t
 
-    pure module function construct(scalar_1D_initializer, order, cells, domain) result(cell_centers_extended)
+    pure module function construct(scalar_1D_initializer, order, cells, x_min, x_max) result(cell_centers_extended)
       !! Result is a collection of cell-centered-extended values with a corresponding mimetic gradient operator
       implicit none
       class(scalar_1D_initializer_t), intent(in) :: scalar_1D_initializer !! elemental initialization function hook
       integer, intent(in) :: order !! order of accuracy
-      integer, intent(in) :: cells !! number of grid cells spanning domain
-      double precision, intent(in) :: domain(:) !! [grid minimum, grid maximum]
+      integer, intent(in) :: cells !! number of grid cells spanning the domain
+      double precision, intent(in) :: x_min !! grid location minimum
+      double precision, intent(in) :: x_max !! grid location maximum
       type(cell_centers_extended_t) cell_centers_extended
     end function
 
diff --git a/src/fortran/cell_centers_extended_s.F90 b/src/fortran/cell_centers_extended_s.F90
index 98ab9eb0..1d7626d8 100644
--- a/src/fortran/cell_centers_extended_s.F90
+++ b/src/fortran/cell_centers_extended_s.F90
@@ -10,9 +10,9 @@
 
     integer cell
 
-    call_julienne_assert(size(domain) .equalsExpected. 2)
+    call_julienne_assert(x_max .isGreaterThan. x_min)
 
-    associate(x_min => domain(1), x_max => domain(2), dx => dble(domain(2) - domain(1))/dble(cells))
+    associate(dx => dble(x_max - x_min)/dble(cells))
       cell_centers_extended%grid_ = [x_min, x_min + dx/2. + [((cell-1)*dx, cell = 1, cells)], x_max] ! boundaries + cell centers as described in
       cell_centers_extended%scalar_1D_ = scalar_1D_initializer%f(cell_centers_extended%grid_)        ! Corbino & Castillo (2020)
       cell_centers_extended%gradient_operator_ = gradient_operator_t(k=order, dx=dx, m=cells)        ! https://doi.org/10.1016/j.cam.2019.06.042
diff --git a/test/gradient_operator_test_m.F90 b/test/gradient_operator_test_m.F90
index 76ae2c32..32f3defe 100644
--- a/test/gradient_operator_test_m.F90
+++ b/test/gradient_operator_test_m.F90
@@ -79,7 +79,7 @@ function check_grad_const() result(test_diagnosis)
     type(const_initializer_1D_t) :: const_initializer_1D
     double precision, parameter :: df_dx = 0D0
 
-    grad_f = .grad. cell_centers_extended_t(const_initializer_1D, order=2, cells=4, domain=[0D0,1D0]) ! gfortran blocks use of association
+    grad_f = .grad. cell_centers_extended_t(const_initializer_1D, order=2, cells=4, x_min=0D0, x_max=1D0) ! gfortran blocks use of association
     test_diagnosis = .all. (grad_f%values() .approximates. df_dx .within. tight_tolerance) // " (d(const)/dx)"
   end function
 
@@ -99,7 +99,7 @@ function check_grad_line() result(test_diagnosis)
     type(line_initializer_1D_t) :: line_initializer_1D
     double precision, parameter :: df_dx = 14D0
 
-    grad_f = .grad. cell_centers_extended_t(line_initializer_1D, order=2, cells=4, domain=[0D0,1D0]) ! gfortran blocks use of association
+    grad_f = .grad. cell_centers_extended_t(line_initializer_1D, order=2, cells=4, x_min=0D0, x_max=1D0) ! gfortran blocks use of association
     test_diagnosis = .all. (grad_f%values() .approximates. df_dx .within. loose_tolerance) // " (d(line)/dx)"
 
   end function
@@ -120,7 +120,7 @@ function check_grad_parabola() result(test_diagnosis)
     end type
     type(parabola_initializer_1D_t) parabola_initializer_1D
 
-    quadratic = cell_centers_extended_t(parabola_initializer_1D, order=2, cells=4, domain=[0D0,1D0]) ! gfortran blocks use of association
+    quadratic = cell_centers_extended_t(parabola_initializer_1D, order=2, cells=4, x_min=0D0, x_max=1D0) ! gfortran blocks use of association
     grad_f = .grad. quadratic ! gfortran blocks use of association
     print *, "grad_f = ", grad_f%values()
     print *, "df_dx = ",parabola(quadratic%grid_)

From 2c09eb8ef92dc90985786e864a8de6efa26f4d33 Mon Sep 17 00:00:00 2001
From: Damian Rouson <rouson@lbl.gov>
Date: Tue, 7 Oct 2025 18:01:26 -0400
Subject: [PATCH 09/40] feat(cell_centers_ex): add cells_ component

---
 src/fortran/cell_centers_extended_m.f90 | 1 +
 src/fortran/cell_centers_extended_s.F90 | 1 +
 2 files changed, 2 insertions(+)

diff --git a/src/fortran/cell_centers_extended_m.f90 b/src/fortran/cell_centers_extended_m.f90
index 25662572..1b45fd36 100644
--- a/src/fortran/cell_centers_extended_m.f90
+++ b/src/fortran/cell_centers_extended_m.f90
@@ -49,6 +49,7 @@ pure module function values(self) result(gradients)
     !private
     double precision, allocatable :: scalar_1D_(:), grid_(:)
     double precision x_min_, x_max_
+    integer cells_
     type(gradient_operator_t) gradient_operator_
   contains
     generic :: operator(.grad.) => grad
diff --git a/src/fortran/cell_centers_extended_s.F90 b/src/fortran/cell_centers_extended_s.F90
index 1d7626d8..b2387848 100644
--- a/src/fortran/cell_centers_extended_s.F90
+++ b/src/fortran/cell_centers_extended_s.F90
@@ -14,6 +14,7 @@
 
     associate(dx => dble(x_max - x_min)/dble(cells))
       cell_centers_extended%grid_ = [x_min, x_min + dx/2. + [((cell-1)*dx, cell = 1, cells)], x_max] ! boundaries + cell centers as described in
+      cell_centers_extended%cells_ = cells
       cell_centers_extended%scalar_1D_ = scalar_1D_initializer%f(cell_centers_extended%grid_)        ! Corbino & Castillo (2020)
       cell_centers_extended%gradient_operator_ = gradient_operator_t(k=order, dx=dx, m=cells)        ! https://doi.org/10.1016/j.cam.2019.06.042
     end associate

From 478f8b182206191a77ada46113b87213e9b1b9f6 Mon Sep 17 00:00:00 2001
From: Damian Rouson <rouson@lbl.gov>
Date: Tue, 7 Oct 2025 18:14:34 -0400
Subject: [PATCH 10/40] refac(cell_centers_extended): rm grid_ component

---
 src/fortran/cell_centers_extended_m.f90 |  9 +++++++++
 src/fortran/cell_centers_extended_s.F90 | 20 +++++++++++++++-----
 test/gradient_operator_test_m.F90       |  2 +-
 3 files changed, 25 insertions(+), 6 deletions(-)

diff --git a/src/fortran/cell_centers_extended_m.f90 b/src/fortran/cell_centers_extended_m.f90
index 1b45fd36..d813c43c 100644
--- a/src/fortran/cell_centers_extended_m.f90
+++ b/src/fortran/cell_centers_extended_m.f90
@@ -52,6 +52,7 @@ pure module function values(self) result(gradients)
     integer cells_
     type(gradient_operator_t) gradient_operator_
   contains
+    procedure grid
     generic :: operator(.grad.) => grad
     procedure, non_overridable, private :: grad
   end type
@@ -73,6 +74,14 @@ pure module function construct(scalar_1D_initializer, order, cells, x_min, x_max
 
   interface
 
+    pure module function grid(self) result(x)
+      !! Result is array of cell-centers-extended grid locations (cell centers + boundaries) 
+      !! as described in Corbino & Castillo (2020) https://doi.org/10.1016/j.cam.2019.06.042
+      implicit none
+      class(cell_centers_extended_t), intent(in) :: self
+      double precision, allocatable :: x(:)
+    end function
+
     pure module function grad(self) result(grad_f)
       !! Result is mimetic gradient of f
       implicit none
diff --git a/src/fortran/cell_centers_extended_s.F90 b/src/fortran/cell_centers_extended_s.F90
index b2387848..da51b778 100644
--- a/src/fortran/cell_centers_extended_s.F90
+++ b/src/fortran/cell_centers_extended_s.F90
@@ -12,13 +12,23 @@
 
     call_julienne_assert(x_max .isGreaterThan. x_min)
 
-    associate(dx => dble(x_max - x_min)/dble(cells))
-      cell_centers_extended%grid_ = [x_min, x_min + dx/2. + [((cell-1)*dx, cell = 1, cells)], x_max] ! boundaries + cell centers as described in
-      cell_centers_extended%cells_ = cells
-      cell_centers_extended%scalar_1D_ = scalar_1D_initializer%f(cell_centers_extended%grid_)        ! Corbino & Castillo (2020)
-      cell_centers_extended%gradient_operator_ = gradient_operator_t(k=order, dx=dx, m=cells)        ! https://doi.org/10.1016/j.cam.2019.06.042
+    associate(dx => (x_max - x_min)/cells)
+      cell_centers_extended%gradient_operator_ = gradient_operator_t(k=order, dx=dx, m=cells)
     end associate
 
+    cell_centers_extended%x_min_ = x_min
+    cell_centers_extended%x_max_ = x_max
+    cell_centers_extended%cells_ = cells
+    cell_centers_extended%scalar_1D_ = scalar_1D_initializer%f(cell_centers_extended%grid())
+
+  end procedure
+
+  module procedure grid
+    integer cell
+
+    associate(dx => (self%x_max_ - self%x_min_)/self%cells_)
+      x = [self%x_min_, self%x_min_ + dx/2. + [((cell-1)*dx, cell = 1, self%cells_)], self%x_max_]
+    end associate
   end procedure
 
   module procedure grad
diff --git a/test/gradient_operator_test_m.F90 b/test/gradient_operator_test_m.F90
index 32f3defe..3ad07160 100644
--- a/test/gradient_operator_test_m.F90
+++ b/test/gradient_operator_test_m.F90
@@ -123,7 +123,7 @@ function check_grad_parabola() result(test_diagnosis)
     quadratic = cell_centers_extended_t(parabola_initializer_1D, order=2, cells=4, x_min=0D0, x_max=1D0) ! gfortran blocks use of association
     grad_f = .grad. quadratic ! gfortran blocks use of association
     print *, "grad_f = ", grad_f%values()
-    print *, "df_dx = ",parabola(quadratic%grid_)
+    print *, "df_dx = ",parabola(quadratic%grid())
 
     test_diagnosis = test_diagnosis_t(.true., "") !.all. (grad_f%values() .approximates. df_dx .within. loose_tolerance) // " (d(line)/dx)"
 

From a5b551fe87d29c1950d5795cdea03814e308152c Mon Sep 17 00:00:00 2001
From: Damian Rouson <rouson@lbl.gov>
Date: Wed, 8 Oct 2025 10:46:54 -0400
Subject: [PATCH 11/40] test(gradient): unit test of d(parabola)/dx passes

---
 src/fortran/cell_centers_extended_m.f90 | 35 +++++++++++++++++++------
 src/fortran/cell_centers_extended_s.F90 | 27 ++++++++++---------
 src/fortran/gradient_s.F90              | 17 +++++++++++-
 src/fortran/mimetic_matrix_s.F90        |  2 +-
 test/gradient_operator_test_m.F90       | 30 ++++++++++++---------
 5 files changed, 76 insertions(+), 35 deletions(-)

diff --git a/src/fortran/cell_centers_extended_m.f90 b/src/fortran/cell_centers_extended_m.f90
index d813c43c..a4c2d03e 100644
--- a/src/fortran/cell_centers_extended_m.f90
+++ b/src/fortran/cell_centers_extended_m.f90
@@ -13,9 +13,6 @@ module cell_centers_extended_m
     !! Encapsulate a mimetic matrix with a corresponding matrix-vector product operator
     private
     double precision, allocatable :: upper_(:,:), inner_(:), lower_(:,:)
-  contains
-    generic :: operator(.x.) => matvec
-    procedure, private, non_overridable :: matvec
   end type
 
   type gradient_operator_t
@@ -29,13 +26,35 @@ module cell_centers_extended_m
   type gradient_t
     !! Encapsulate gradient values produced only by .grad. (no other constructors)
     private
-    double precision, allocatable :: g_(:)
+    double precision, allocatable :: vector_1D_(:) !! gradient values at cell faces (nodes in 1D)
+    double precision x_min_ !! domain lower boundary
+    double precision x_max_ !! domain upper boundary
+    integer cells_ !! number of grid cells spanning the domain
   contains
     procedure values
+    procedure faces
   end type
 
+  interface gradient_t
+
+      pure module function construct_gradient(face_centered_values, x_min, x_max, cells) result(gradient)
+        !! Result is an object storing gradients at cell faces
+        implicit none
+        double precision, intent(in) :: face_centered_values(:), x_min, x_max
+        integer, intent(in) :: cells
+        type(gradient_t) gradient
+      end function
+
+  end interface
+
   interface 
 
+    pure module function faces(self) result(x)
+      implicit none
+      class(gradient_t), intent(in) :: self
+      double precision, allocatable :: x(:)
+    end function
+
      pure module function values(self) result(gradients)
        implicit none
        class(gradient_t), intent(in) :: self
@@ -46,8 +65,8 @@ pure module function values(self) result(gradients)
 
   type cell_centers_extended_t
     !! Encapsulate information at cell centers and boundaries
-    !private
-    double precision, allocatable :: scalar_1D_(:), grid_(:)
+    private
+    double precision, allocatable :: scalar_1D_(:)
     double precision x_min_, x_max_
     integer cells_
     type(gradient_operator_t) gradient_operator_
@@ -117,12 +136,12 @@ pure module function construct_from_components(upper, inner, lower) result(mimet
 
   interface
 
-    pure module function matvec(self, vector) result(gradient)
+    pure module function matvec(self, vector) result(matvec_product)
       !! Apply a matrix operator to a vector
       implicit none
       class(mimetic_matrix_t), intent(in) :: self
       type(cell_centers_extended_t), intent(in) :: vector
-      type(gradient_t) gradient
+      double precision, allocatable :: matvec_product(:)
     end function
 
   end interface
diff --git a/src/fortran/cell_centers_extended_s.F90 b/src/fortran/cell_centers_extended_s.F90
index da51b778..9c3367d3 100644
--- a/src/fortran/cell_centers_extended_s.F90
+++ b/src/fortran/cell_centers_extended_s.F90
@@ -7,32 +7,33 @@
 contains
 
   module procedure construct
-
-    integer cell
-
     call_julienne_assert(x_max .isGreaterThan. x_min)
-
-    associate(dx => (x_max - x_min)/cells)
-      cell_centers_extended%gradient_operator_ = gradient_operator_t(k=order, dx=dx, m=cells)
-    end associate
+    call_julienne_assert(cells .isAtLeast. 2*order)
 
     cell_centers_extended%x_min_ = x_min
     cell_centers_extended%x_max_ = x_max
     cell_centers_extended%cells_ = cells
-    cell_centers_extended%scalar_1D_ = scalar_1D_initializer%f(cell_centers_extended%grid())
-
+    cell_centers_extended%gradient_operator_ = gradient_operator_t(k=order, dx=(x_max - x_min)/cells, m=cells)
+    cell_centers_extended%scalar_1D_ = scalar_1D_initializer%f(grid_(x_min, x_max, cells))
   end procedure
 
-  module procedure grid
+  pure function grid_(x_min, x_max, cells) result(x)
+    double precision, intent(in) :: x_min, x_max
+    integer, intent(in) :: cells
+    double precision, allocatable :: x(:)
     integer cell
 
-    associate(dx => (self%x_max_ - self%x_min_)/self%cells_)
-      x = [self%x_min_, self%x_min_ + dx/2. + [((cell-1)*dx, cell = 1, self%cells_)], self%x_max_]
+    associate(dx => (x_max - x_min)/cells)
+      x = [x_min, x_min + dx/2. + [((cell-1)*dx, cell = 1, cells)], x_max]
     end associate
+  end function
+
+  module procedure grid
+      x = grid_(self%x_min_, self%x_max_, self%cells_)
   end procedure
 
   module procedure grad
-    grad_f = self%gradient_operator_%mimetic_matrix_ .x. self
+    grad_f = gradient_t(matvec(self%gradient_operator_%mimetic_matrix_, self), self%x_min_, self%x_max_, self%cells_)
   end procedure
 
 end submodule cell_centers_extended_s
\ No newline at end of file
diff --git a/src/fortran/gradient_s.F90 b/src/fortran/gradient_s.F90
index 93c7b410..0c380af0 100644
--- a/src/fortran/gradient_s.F90
+++ b/src/fortran/gradient_s.F90
@@ -3,8 +3,23 @@
 
 contains
 
+  module procedure construct_gradient
+    gradient%vector_1D_ = face_centered_values
+    gradient%x_min_ = x_min
+    gradient%x_max_ = x_max
+    gradient%cells_ = cells
+  end procedure
+
   module procedure values
-    gradients =  self%g_
+    gradients =  self%vector_1D_
+  end procedure
+
+  module procedure faces
+    integer cell
+    x = [ self%x_min_ &
+         ,self%x_min_ + [(cell*(self%x_max_ - self%x_min_)/self%cells_, cell = 1, self%cells_-1)] &
+         ,self%x_max_ &
+        ]
   end procedure
 
 end submodule gradient_s
\ No newline at end of file
diff --git a/src/fortran/mimetic_matrix_s.F90 b/src/fortran/mimetic_matrix_s.F90
index cb89443a..e791e558 100644
--- a/src/fortran/mimetic_matrix_s.F90
+++ b/src/fortran/mimetic_matrix_s.F90
@@ -28,7 +28,7 @@
           product_inner(row) = dot_product(self%inner_, vector%scalar_1D_(row + 1 : row + size(self%inner_)))
         end do
 
-        gradient%g_ = [ &
+        matvec_product = [ &
            matmul(self%upper_, vector%scalar_1D_(1 : size(self%upper_,2))) &
           ,product_inner &
           ,matmul(self%lower_, vector%scalar_1D_(size(vector%scalar_1D_) - size(self%lower_,2) + 1 : )) &
diff --git a/test/gradient_operator_test_m.F90 b/test/gradient_operator_test_m.F90
index 3ad07160..050b6531 100644
--- a/test/gradient_operator_test_m.F90
+++ b/test/gradient_operator_test_m.F90
@@ -27,7 +27,7 @@ module gradient_operator_test_m
     procedure, nopass :: results
   end type
 
-  double precision, parameter :: tight_tolerance = 1D-14, loose_tolerance = 1D-10
+  double precision, parameter :: tight_tolerance = 1D-14, loose_tolerance = 1D-12
 
 contains
 
@@ -42,8 +42,9 @@ function results() result(test_results)
     type(gradient_operator_test_t) gradient_operator_test
     type(test_result_t), allocatable :: test_results(:)
     test_results = gradient_operator_test%run([ & 
-       test_description_t('computing the gradient of a constant function within a tolerance of ' // string_t(tight_tolerance), check_grad_const) &
+       test_description_t('computing the gradient of a constant within a tolerance of ' // string_t(tight_tolerance), check_grad_const) &
       ,test_description_t('computing the gradient of a line within a tolerance of ' // string_t(loose_tolerance), check_grad_line) &
+      ,test_description_t('computing the gradient of a parabola within a tolerance of ' // string_t(loose_tolerance), check_grad_parabola) &
     ])
   end function
 
@@ -55,9 +56,12 @@ function results() result(test_results)
     procedure(diagnosis_function_i), pointer :: &
        check_grad_const_ptr => check_grad_const &
       ,check_grad_line_ptr => check_grad_line
+      ,check_grad_parabola_ptr => check_grad_parabola
+
     test_results = gradient_operator_test%run([ &
-       test_description_t('computing the gradient of a constant function within a tolerance of ' // string_t(tight_tolerance), check_grad_const_ptr) &
+       test_description_t('computing the gradient of a constant within a tolerance of ' // string_t(tight_tolerance), check_grad_const_ptr) &
       ,test_description_t('computing the gradient of a line within a tolerance of ' // string_t(loose_tolerance), check_grad_line_ptr) &
+      ,test_description_t('computing the gradient of a parabola within a tolerance of ' // string_t(loose_tolerance), check_grad_parabola_ptr) &
     ])
   end function
 
@@ -99,7 +103,7 @@ function check_grad_line() result(test_diagnosis)
     type(line_initializer_1D_t) :: line_initializer_1D
     double precision, parameter :: df_dx = 14D0
 
-    grad_f = .grad. cell_centers_extended_t(line_initializer_1D, order=2, cells=4, x_min=0D0, x_max=1D0) ! gfortran blocks use of association
+    grad_f = .grad. cell_centers_extended_t(line_initializer_1D, order=2, cells=4, x_min=0D0, x_max=1D0)
     test_diagnosis = .all. (grad_f%values() .approximates. df_dx .within. loose_tolerance) // " (d(line)/dx)"
 
   end function
@@ -112,20 +116,22 @@ elemental function parabola(x) result(y)
 
   function check_grad_parabola() result(test_diagnosis)
     type(test_diagnosis_t) test_diagnosis
-    type(gradient_t) grad_f
-    type(cell_centers_extended_t) quadratic
     type, extends(scalar_1D_initializer_t) :: parabola_initializer_1D_t
     contains
-      procedure, nopass, non_overridable :: f => line
+      procedure, nopass, non_overridable :: f => parabola
     end type
     type(parabola_initializer_1D_t) parabola_initializer_1D
+    type(cell_centers_extended_t) quadratic
+    type(gradient_t) grad_f
 
-    quadratic = cell_centers_extended_t(parabola_initializer_1D, order=2, cells=4, x_min=0D0, x_max=1D0) ! gfortran blocks use of association
-    grad_f = .grad. quadratic ! gfortran blocks use of association
-    print *, "grad_f = ", grad_f%values()
-    print *, "df_dx = ",parabola(quadratic%grid())
+    quadratic = cell_centers_extended_t(parabola_initializer_1D, order=2, cells=4, x_min=0D0, x_max=1D0)
+    grad_f = .grad. quadratic
 
-    test_diagnosis = test_diagnosis_t(.true., "") !.all. (grad_f%values() .approximates. df_dx .within. loose_tolerance) // " (d(line)/dx)"
+    associate(x => grad_f%faces())
+      associate(df_dx => 14*x + 3)
+        test_diagnosis = .all. (grad_f%values() .approximates. df_dx .within. loose_tolerance) // " (d(parabola)/dx)"
+      end associate
+    end associate
 
   end function
 

From bd0bb76e20f122ef24a6442ec53520e4dd323e7d Mon Sep 17 00:00:00 2001
From: Damian Rouson <rouson@lbl.gov>
Date: Wed, 8 Oct 2025 16:44:11 -0400
Subject: [PATCH 12/40] refac(initializers): rm abstract initializer type

---
 src/fortran/cell_centers_extended_m.f90 | 16 ++++++--
 src/fortran/cell_centers_extended_s.F90 | 11 +++++
 src/fortran/initializers_m.f90          | 28 -------------
 src/fortran/mole_m.f90                  |  3 +-
 test/gradient_operator_test_m.F90       | 54 ++++++++++---------------
 5 files changed, 46 insertions(+), 66 deletions(-)
 delete mode 100644 src/fortran/initializers_m.f90

diff --git a/src/fortran/cell_centers_extended_m.f90 b/src/fortran/cell_centers_extended_m.f90
index a4c2d03e..159f8845 100644
--- a/src/fortran/cell_centers_extended_m.f90
+++ b/src/fortran/cell_centers_extended_m.f90
@@ -1,13 +1,23 @@
 module cell_centers_extended_m
   !! Define an abstraction for the collection of points used to compute gradidents:
   !! cell centers plus oundaries.
-  use initializers_m, only : scalar_1D_initializer_t
   implicit none
 
   private
   public :: cell_centers_extended_t
   public :: gradient_t
   public :: gradient_operator_t
+  public :: scalar_1D_initializer_i
+
+  abstract interface
+
+    pure function scalar_1D_initializer_i(x) result(f)
+      implicit none
+      double precision, intent(in) :: x(:)
+      double precision, allocatable :: f(:)
+    end function
+
+  end interface
 
   type mimetic_matrix_t
     !! Encapsulate a mimetic matrix with a corresponding matrix-vector product operator
@@ -78,10 +88,10 @@ pure module function values(self) result(gradients)
 
   interface cell_centers_extended_t
 
-    pure module function construct(scalar_1D_initializer, order, cells, x_min, x_max) result(cell_centers_extended)
+    pure module function construct_from_function(initializer, order, cells, x_min, x_max) result(cell_centers_extended)
       !! Result is a collection of cell-centered-extended values with a corresponding mimetic gradient operator
       implicit none
-      class(scalar_1D_initializer_t), intent(in) :: scalar_1D_initializer !! elemental initialization function hook
+      procedure(scalar_1D_initializer_i), pointer :: initializer 
       integer, intent(in) :: order !! order of accuracy
       integer, intent(in) :: cells !! number of grid cells spanning the domain
       double precision, intent(in) :: x_min !! grid location minimum
diff --git a/src/fortran/cell_centers_extended_s.F90 b/src/fortran/cell_centers_extended_s.F90
index 9c3367d3..18565814 100644
--- a/src/fortran/cell_centers_extended_s.F90
+++ b/src/fortran/cell_centers_extended_s.F90
@@ -17,6 +17,17 @@
     cell_centers_extended%scalar_1D_ = scalar_1D_initializer%f(grid_(x_min, x_max, cells))
   end procedure
 
+  module procedure construct_from_function
+    call_julienne_assert(x_max .isGreaterThan. x_min)
+    call_julienne_assert(cells .isAtLeast. 2*order)
+
+    cell_centers_extended%x_min_ = x_min
+    cell_centers_extended%x_max_ = x_max
+    cell_centers_extended%cells_ = cells
+    cell_centers_extended%gradient_operator_ = gradient_operator_t(k=order, dx=(x_max - x_min)/cells, m=cells)
+    cell_centers_extended%scalar_1D_ = initializer(grid_(x_min, x_max, cells))
+  end procedure
+
   pure function grid_(x_min, x_max, cells) result(x)
     double precision, intent(in) :: x_min, x_max
     integer, intent(in) :: cells
diff --git a/src/fortran/initializers_m.f90 b/src/fortran/initializers_m.f90
deleted file mode 100644
index 88fbe1a1..00000000
--- a/src/fortran/initializers_m.f90
+++ /dev/null
@@ -1,28 +0,0 @@
-module initializers_m
-  !! Implement a workaround for the Fortran standard's prohibition against
-  !! elemental procedures as dummy arguments. Users can extend the abstract
-  !! type(s) in this module and define the deferred binding as an elemental
-  !! function for use in initializing variables at grid locations without
-  !! requiring loops.
-  implicit none
-
-  private
-  public :: scalar_1D_initializer_t
-
-  abstract interface
-
-    elemental function scalar_1D_initializer_i(x) result(f)
-      implicit none
-      double precision, intent(in) :: x 
-      double precision f
-    end function
-
-  end interface
-
-  type, abstract :: scalar_1D_initializer_t
-    !! Define a hook on which to hang elemental grid-variable initializers
-  contains
-    procedure(scalar_1D_initializer_i), deferred, nopass :: f
-  end type
-
-end module initializers_m
diff --git a/src/fortran/mole_m.f90 b/src/fortran/mole_m.f90
index a97aaee9..e651a00a 100644
--- a/src/fortran/mole_m.f90
+++ b/src/fortran/mole_m.f90
@@ -1,6 +1,5 @@
 module mole_m 
   !! MOLE Fortran public entities
-  use cell_centers_extended_m, only : cell_centers_extended_t, gradient_t
-  use initializers_m, only : scalar_1D_initializer_t
+  use cell_centers_extended_m, only : cell_centers_extended_t, gradient_t, scalar_1D_initializer_i
   implicit none
 end module mole_m
diff --git a/test/gradient_operator_test_m.F90 b/test/gradient_operator_test_m.F90
index 050b6531..d161a264 100644
--- a/test/gradient_operator_test_m.F90
+++ b/test/gradient_operator_test_m.F90
@@ -15,7 +15,7 @@ module gradient_operator_test_m
     ,operator(.all.) &
     ,operator(.approximates.) &
     ,operator(.within.)
-  use mole_m, only : cell_centers_extended_t, gradient_t, scalar_1D_initializer_t
+  use mole_m, only : cell_centers_extended_t, gradient_t, scalar_1D_initializer_i
 #if ! HAVE_PROCEDURE_ACTUAL_FOR_POINTER_DUMMY
   use julienne_m, only : diagnosis_function_i
 #endif
@@ -67,64 +67,52 @@ function results() result(test_results)
 
 #endif
 
-  elemental function const(x) result(c)
-    double precision, intent(in) :: x
-    double precision c
-    c = 3D0
+  pure function const(x) result(y)
+    double precision, intent(in) :: x(:)
+    double precision, allocatable :: y(:)
+    integer i
+    y = [(5D0, i=1,size(x))]
   end function
 
   function check_grad_const() result(test_diagnosis)
     type(test_diagnosis_t) test_diagnosis
     type(gradient_t) grad_f
-    type, extends(scalar_1D_initializer_t) :: const_initializer_1D_t
-    contains
-      procedure, nopass, non_overridable :: f => const
-    end type
-    type(const_initializer_1D_t) :: const_initializer_1D
-    double precision, parameter :: df_dx = 0D0
-
-    grad_f = .grad. cell_centers_extended_t(const_initializer_1D, order=2, cells=4, x_min=0D0, x_max=1D0) ! gfortran blocks use of association
-    test_diagnosis = .all. (grad_f%values() .approximates. df_dx .within. tight_tolerance) // " (d(const)/dx)"
+    double precision, parameter :: df_dx = 0.
+    procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => const
+
+    grad_f = .grad. cell_centers_extended_t(scalar_1D_initializer, order=2, cells=4, x_min=0D0, x_max=1D0)
+    test_diagnosis = .all. (grad_f%values() .approximates. df_dx .within. loose_tolerance) // " (d(line)/dx)"
   end function
 
-  elemental function line(x) result(y)
-    double precision, intent(in) :: x
-    double precision y
+  pure function line(x) result(y)
+    double precision, intent(in) :: x(:)
+    double precision, allocatable :: y(:)
     y = 14*x + 3
   end function
 
   function check_grad_line() result(test_diagnosis)
     type(test_diagnosis_t) test_diagnosis
     type(gradient_t) grad_f
-    type, extends(scalar_1D_initializer_t) :: line_initializer_1D_t
-    contains
-      procedure, nopass, non_overridable :: f => line
-    end type
-    type(line_initializer_1D_t) :: line_initializer_1D
     double precision, parameter :: df_dx = 14D0
+    procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => line
 
-    grad_f = .grad. cell_centers_extended_t(line_initializer_1D, order=2, cells=4, x_min=0D0, x_max=1D0)
+    grad_f = .grad. cell_centers_extended_t(scalar_1D_initializer, order=2, cells=4, x_min=0D0, x_max=1D0)
     test_diagnosis = .all. (grad_f%values() .approximates. df_dx .within. loose_tolerance) // " (d(line)/dx)"
-
   end function
 
-  elemental function parabola(x) result(y)
-    double precision, intent(in) :: x
-    double precision y
+  pure function parabola(x) result(y)
+    double precision, intent(in) :: x(:)
+    double precision, allocatable :: y(:)
     y = 7*x**2 + 3*x + 5
   end function
 
   function check_grad_parabola() result(test_diagnosis)
     type(test_diagnosis_t) test_diagnosis
-    type, extends(scalar_1D_initializer_t) :: parabola_initializer_1D_t
-    contains
-      procedure, nopass, non_overridable :: f => parabola
-    end type
-    type(parabola_initializer_1D_t) parabola_initializer_1D
     type(cell_centers_extended_t) quadratic
     type(gradient_t) grad_f
+    procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => parabola
 
-    quadratic = cell_centers_extended_t(parabola_initializer_1D, order=2, cells=4, x_min=0D0, x_max=1D0)
+    quadratic = cell_centers_extended_t(scalar_1D_initializer , order=2, cells=4, x_min=0D0, x_max=1D0)
     grad_f = .grad. quadratic
 
     associate(x => grad_f%faces())

From 6653b8b86b96a90ed611b2a890a0513f804f2fff Mon Sep 17 00:00:00 2001
From: Damian Rouson <rouson@lbl.gov>
Date: Wed, 8 Oct 2025 19:16:50 -0400
Subject: [PATCH 13/40] fix(cell_centers_extended_s): rm dead code

---
 src/fortran/cell_centers_extended_s.F90 | 11 -----------
 1 file changed, 11 deletions(-)

diff --git a/src/fortran/cell_centers_extended_s.F90 b/src/fortran/cell_centers_extended_s.F90
index 18565814..6fd589b3 100644
--- a/src/fortran/cell_centers_extended_s.F90
+++ b/src/fortran/cell_centers_extended_s.F90
@@ -6,17 +6,6 @@
 
 contains
 
-  module procedure construct
-    call_julienne_assert(x_max .isGreaterThan. x_min)
-    call_julienne_assert(cells .isAtLeast. 2*order)
-
-    cell_centers_extended%x_min_ = x_min
-    cell_centers_extended%x_max_ = x_max
-    cell_centers_extended%cells_ = cells
-    cell_centers_extended%gradient_operator_ = gradient_operator_t(k=order, dx=(x_max - x_min)/cells, m=cells)
-    cell_centers_extended%scalar_1D_ = scalar_1D_initializer%f(grid_(x_min, x_max, cells))
-  end procedure
-
   module procedure construct_from_function
     call_julienne_assert(x_max .isGreaterThan. x_min)
     call_julienne_assert(cells .isAtLeast. 2*order)

From aa036a831efb7976045f0fe4460f67c120696b42 Mon Sep 17 00:00:00 2001
From: Damian Rouson <rouson@lbl.gov>
Date: Thu, 9 Oct 2025 00:25:18 -0400
Subject: [PATCH 14/40] refac(gradient_t): mv to separate module/submodule

---
 src/fortran/cell_centers_extended_m.f90 | 42 +--------------------
 src/fortran/gradient_m.f90              | 49 +++++++++++++++++++++++++
 src/fortran/gradient_s.F90              |  4 +-
 src/fortran/mole_m.f90                  |  3 +-
 4 files changed, 54 insertions(+), 44 deletions(-)
 create mode 100644 src/fortran/gradient_m.f90

diff --git a/src/fortran/cell_centers_extended_m.f90 b/src/fortran/cell_centers_extended_m.f90
index 159f8845..a67827e8 100644
--- a/src/fortran/cell_centers_extended_m.f90
+++ b/src/fortran/cell_centers_extended_m.f90
@@ -1,11 +1,11 @@
 module cell_centers_extended_m
   !! Define an abstraction for the collection of points used to compute gradidents:
   !! cell centers plus oundaries.
+  use gradient_m, only : gradient_t
   implicit none
 
   private
   public :: cell_centers_extended_t
-  public :: gradient_t
   public :: gradient_operator_t
   public :: scalar_1D_initializer_i
 
@@ -33,46 +33,6 @@ pure function scalar_1D_initializer_i(x) result(f)
     type(mimetic_matrix_t) mimetic_matrix_
   end type
 
-  type gradient_t
-    !! Encapsulate gradient values produced only by .grad. (no other constructors)
-    private
-    double precision, allocatable :: vector_1D_(:) !! gradient values at cell faces (nodes in 1D)
-    double precision x_min_ !! domain lower boundary
-    double precision x_max_ !! domain upper boundary
-    integer cells_ !! number of grid cells spanning the domain
-  contains
-    procedure values
-    procedure faces
-  end type
-
-  interface gradient_t
-
-      pure module function construct_gradient(face_centered_values, x_min, x_max, cells) result(gradient)
-        !! Result is an object storing gradients at cell faces
-        implicit none
-        double precision, intent(in) :: face_centered_values(:), x_min, x_max
-        integer, intent(in) :: cells
-        type(gradient_t) gradient
-      end function
-
-  end interface
-
-  interface 
-
-    pure module function faces(self) result(x)
-      implicit none
-      class(gradient_t), intent(in) :: self
-      double precision, allocatable :: x(:)
-    end function
-
-     pure module function values(self) result(gradients)
-       implicit none
-       class(gradient_t), intent(in) :: self
-       double precision, allocatable :: gradients(:)
-     end function
-
-  end interface
-
   type cell_centers_extended_t
     !! Encapsulate information at cell centers and boundaries
     private
diff --git a/src/fortran/gradient_m.f90 b/src/fortran/gradient_m.f90
new file mode 100644
index 00000000..3c669a9b
--- /dev/null
+++ b/src/fortran/gradient_m.f90
@@ -0,0 +1,49 @@
+module gradient_m
+  !! Define an abstraction for the collection of points used to compute gradidents:
+  !! cell centers plus oundaries.
+  implicit none
+
+  private
+  public :: gradient_t
+
+  type gradient_t
+    !! Encapsulate gradient values produced only by .grad. (no other constructors)
+    private
+    double precision, allocatable :: vector_1D_(:) !! gradient values at cell faces (nodes in 1D)
+    double precision x_min_ !! domain lower boundary
+    double precision x_max_ !! domain upper boundary
+    integer cells_ !! number of grid cells spanning the domain
+  contains
+    procedure values
+    procedure faces
+  end type
+
+  interface gradient_t
+
+      pure module function construct_from_components(face_centered_values, x_min, x_max, cells) result(gradient)
+        !! Result is an object storing gradients at cell faces
+        implicit none
+        double precision, intent(in) :: face_centered_values(:), x_min, x_max
+        integer, intent(in) :: cells
+        type(gradient_t) gradient
+      end function
+
+  end interface
+
+  interface
+
+    pure module function faces(self) result(x)
+      implicit none
+      class(gradient_t), intent(in) :: self
+      double precision, allocatable :: x(:)
+    end function
+
+     pure module function values(self) result(gradients)
+       implicit none
+       class(gradient_t), intent(in) :: self
+       double precision, allocatable :: gradients(:)
+     end function
+
+  end interface
+
+end module gradient_m
\ No newline at end of file
diff --git a/src/fortran/gradient_s.F90 b/src/fortran/gradient_s.F90
index 0c380af0..125262c9 100644
--- a/src/fortran/gradient_s.F90
+++ b/src/fortran/gradient_s.F90
@@ -1,9 +1,9 @@
-submodule(cell_centers_extended_m) gradient_s
+submodule(gradient_m) gradient_s
   implicit none
 
 contains
 
-  module procedure construct_gradient
+  module procedure construct_from_components
     gradient%vector_1D_ = face_centered_values
     gradient%x_min_ = x_min
     gradient%x_max_ = x_max
diff --git a/src/fortran/mole_m.f90 b/src/fortran/mole_m.f90
index e651a00a..2862f65c 100644
--- a/src/fortran/mole_m.f90
+++ b/src/fortran/mole_m.f90
@@ -1,5 +1,6 @@
 module mole_m 
   !! MOLE Fortran public entities
-  use cell_centers_extended_m, only : cell_centers_extended_t, gradient_t, scalar_1D_initializer_i
+  use cell_centers_extended_m, only : cell_centers_extended_t, scalar_1D_initializer_i
+  use gradient_m, only : gradient_t
   implicit none
 end module mole_m

From cbb27853a5e36aaaa59d0a78ce497c534b4cd7cc Mon Sep 17 00:00:00 2001
From: Damian Rouson <rouson@lbl.gov>
Date: Thu, 9 Oct 2025 00:39:21 -0400
Subject: [PATCH 15/40] refac(cell_centers_extended_t): rename scalar_1D_t

---
 src/fortran/gradient_operator_s.F90           |  2 +-
 src/fortran/mimetic_matrix_s.F90              |  2 +-
 src/fortran/mole_m.f90                        |  2 +-
 ...centers_extended_m.f90 => scalar_1D_m.f90} | 20 ++++++++---------
 ...centers_extended_s.F90 => scalar_1D_s.F90} | 14 ++++++------
 test/gradient_operator_test_m.F90             | 22 +++++++++----------
 6 files changed, 31 insertions(+), 31 deletions(-)
 rename src/fortran/{cell_centers_extended_m.f90 => scalar_1D_m.f90} (88%)
 rename src/fortran/{cell_centers_extended_s.F90 => scalar_1D_s.F90} (68%)

diff --git a/src/fortran/gradient_operator_s.F90 b/src/fortran/gradient_operator_s.F90
index 8514e15e..e649e1ff 100644
--- a/src/fortran/gradient_operator_s.F90
+++ b/src/fortran/gradient_operator_s.F90
@@ -1,7 +1,7 @@
 #include "julienne-assert-macros.h"
 #include "mole-language-support.F90"
 
-submodule(cell_centers_extended_m) gradient_operator_s
+submodule(scalar_1D_m) gradient_operator_s
   use julienne_m, only : call_julienne_assert_, string_t
 #if ASSERTIONS
   use julienne_m, only : operator(.isAtLeast.)
diff --git a/src/fortran/mimetic_matrix_s.F90 b/src/fortran/mimetic_matrix_s.F90
index e791e558..7a6ec067 100644
--- a/src/fortran/mimetic_matrix_s.F90
+++ b/src/fortran/mimetic_matrix_s.F90
@@ -1,7 +1,7 @@
 #include "mole-language-support.F90"
 #include "julienne-assert-macros.h"
 
-submodule(cell_centers_extended_m) mimetic_matrix_s
+submodule(scalar_1D_m) mimetic_matrix_s
   use julienne_m, only : call_julienne_assert_, string_t
   implicit none
 
diff --git a/src/fortran/mole_m.f90 b/src/fortran/mole_m.f90
index 2862f65c..672a2a94 100644
--- a/src/fortran/mole_m.f90
+++ b/src/fortran/mole_m.f90
@@ -1,6 +1,6 @@
 module mole_m 
   !! MOLE Fortran public entities
-  use cell_centers_extended_m, only : cell_centers_extended_t, scalar_1D_initializer_i
+  use scalar_1D_m, only : scalar_1D_t, scalar_1D_initializer_i
   use gradient_m, only : gradient_t
   implicit none
 end module mole_m
diff --git a/src/fortran/cell_centers_extended_m.f90 b/src/fortran/scalar_1D_m.f90
similarity index 88%
rename from src/fortran/cell_centers_extended_m.f90
rename to src/fortran/scalar_1D_m.f90
index a67827e8..21f4d0b5 100644
--- a/src/fortran/cell_centers_extended_m.f90
+++ b/src/fortran/scalar_1D_m.f90
@@ -1,11 +1,11 @@
-module cell_centers_extended_m
+module scalar_1D_m
   !! Define an abstraction for the collection of points used to compute gradidents:
   !! cell centers plus oundaries.
   use gradient_m, only : gradient_t
   implicit none
 
   private
-  public :: cell_centers_extended_t
+  public :: scalar_1D_t
   public :: gradient_operator_t
   public :: scalar_1D_initializer_i
 
@@ -33,7 +33,7 @@ pure function scalar_1D_initializer_i(x) result(f)
     type(mimetic_matrix_t) mimetic_matrix_
   end type
 
-  type cell_centers_extended_t
+  type scalar_1D_t
     !! Encapsulate information at cell centers and boundaries
     private
     double precision, allocatable :: scalar_1D_(:)
@@ -46,9 +46,9 @@ pure function scalar_1D_initializer_i(x) result(f)
     procedure, non_overridable, private :: grad
   end type
 
-  interface cell_centers_extended_t
+  interface scalar_1D_t
 
-    pure module function construct_from_function(initializer, order, cells, x_min, x_max) result(cell_centers_extended)
+    pure module function construct_from_function(initializer, order, cells, x_min, x_max) result(scalar_1D)
       !! Result is a collection of cell-centered-extended values with a corresponding mimetic gradient operator
       implicit none
       procedure(scalar_1D_initializer_i), pointer :: initializer 
@@ -56,7 +56,7 @@ pure module function construct_from_function(initializer, order, cells, x_min, x
       integer, intent(in) :: cells !! number of grid cells spanning the domain
       double precision, intent(in) :: x_min !! grid location minimum
       double precision, intent(in) :: x_max !! grid location maximum
-      type(cell_centers_extended_t) cell_centers_extended
+      type(scalar_1D_t) scalar_1D
     end function
 
   end interface
@@ -67,14 +67,14 @@ pure module function grid(self) result(x)
       !! Result is array of cell-centers-extended grid locations (cell centers + boundaries) 
       !! as described in Corbino & Castillo (2020) https://doi.org/10.1016/j.cam.2019.06.042
       implicit none
-      class(cell_centers_extended_t), intent(in) :: self
+      class(scalar_1D_t), intent(in) :: self
       double precision, allocatable :: x(:)
     end function
 
     pure module function grad(self) result(grad_f)
       !! Result is mimetic gradient of f
       implicit none
-      class(cell_centers_extended_t), intent(in) :: self
+      class(scalar_1D_t), intent(in) :: self
       type(gradient_t) grad_f !! discrete gradient approximation
     end function
 
@@ -110,10 +110,10 @@ pure module function matvec(self, vector) result(matvec_product)
       !! Apply a matrix operator to a vector
       implicit none
       class(mimetic_matrix_t), intent(in) :: self
-      type(cell_centers_extended_t), intent(in) :: vector
+      type(scalar_1D_t), intent(in) :: vector
       double precision, allocatable :: matvec_product(:)
     end function
 
   end interface
 
-end module cell_centers_extended_m
\ No newline at end of file
+end module scalar_1D_m
\ No newline at end of file
diff --git a/src/fortran/cell_centers_extended_s.F90 b/src/fortran/scalar_1D_s.F90
similarity index 68%
rename from src/fortran/cell_centers_extended_s.F90
rename to src/fortran/scalar_1D_s.F90
index 6fd589b3..859feed6 100644
--- a/src/fortran/cell_centers_extended_s.F90
+++ b/src/fortran/scalar_1D_s.F90
@@ -1,6 +1,6 @@
 #include "julienne-assert-macros.h"
 
-submodule(cell_centers_extended_m) cell_centers_extended_s
+submodule(scalar_1D_m) scalar_1D_s
   use julienne_m, only : call_julienne_assert_, operator(.equalsExpected.)
   implicit none
 
@@ -10,11 +10,11 @@
     call_julienne_assert(x_max .isGreaterThan. x_min)
     call_julienne_assert(cells .isAtLeast. 2*order)
 
-    cell_centers_extended%x_min_ = x_min
-    cell_centers_extended%x_max_ = x_max
-    cell_centers_extended%cells_ = cells
-    cell_centers_extended%gradient_operator_ = gradient_operator_t(k=order, dx=(x_max - x_min)/cells, m=cells)
-    cell_centers_extended%scalar_1D_ = initializer(grid_(x_min, x_max, cells))
+    scalar_1D%x_min_ = x_min
+    scalar_1D%x_max_ = x_max
+    scalar_1D%cells_ = cells
+    scalar_1D%gradient_operator_ = gradient_operator_t(k=order, dx=(x_max - x_min)/cells, m=cells)
+    scalar_1D%scalar_1D_ = initializer(grid_(x_min, x_max, cells))
   end procedure
 
   pure function grid_(x_min, x_max, cells) result(x)
@@ -36,4 +36,4 @@ pure function grid_(x_min, x_max, cells) result(x)
     grad_f = gradient_t(matvec(self%gradient_operator_%mimetic_matrix_, self), self%x_min_, self%x_max_, self%cells_)
   end procedure
 
-end submodule cell_centers_extended_s
\ No newline at end of file
+end submodule scalar_1D_s
\ No newline at end of file
diff --git a/test/gradient_operator_test_m.F90 b/test/gradient_operator_test_m.F90
index d161a264..f0674a84 100644
--- a/test/gradient_operator_test_m.F90
+++ b/test/gradient_operator_test_m.F90
@@ -15,7 +15,7 @@ module gradient_operator_test_m
     ,operator(.all.) &
     ,operator(.approximates.) &
     ,operator(.within.)
-  use mole_m, only : cell_centers_extended_t, gradient_t, scalar_1D_initializer_i
+  use mole_m, only : scalar_1D_t, gradient_t, scalar_1D_initializer_i
 #if ! HAVE_PROCEDURE_ACTUAL_FOR_POINTER_DUMMY
   use julienne_m, only : diagnosis_function_i
 #endif
@@ -42,9 +42,9 @@ function results() result(test_results)
     type(gradient_operator_test_t) gradient_operator_test
     type(test_result_t), allocatable :: test_results(:)
     test_results = gradient_operator_test%run([ & 
-       test_description_t('computing the gradient of a constant within a tolerance of ' // string_t(tight_tolerance), check_grad_const) &
-      ,test_description_t('computing the gradient of a line within a tolerance of ' // string_t(loose_tolerance), check_grad_line) &
-      ,test_description_t('computing the gradient of a parabola within a tolerance of ' // string_t(loose_tolerance), check_grad_parabola) &
+       test_description_t('computing the 2nd-order 1D gradient of a constant within tolerance ' // string_t(tight_tolerance), check_grad_const) &
+      ,test_description_t('computing the 2nd-order 1D gradient of a line within tolerance ' // string_t(loose_tolerance), check_grad_line) &
+      ,test_description_t('computing the 2nd-order 1D gradient of a parabola within tolerance ' // string_t(loose_tolerance), check_grad_parabola) &
     ])
   end function
 
@@ -59,9 +59,9 @@ function results() result(test_results)
       ,check_grad_parabola_ptr => check_grad_parabola
 
     test_results = gradient_operator_test%run([ &
-       test_description_t('computing the gradient of a constant within a tolerance of ' // string_t(tight_tolerance), check_grad_const_ptr) &
-      ,test_description_t('computing the gradient of a line within a tolerance of ' // string_t(loose_tolerance), check_grad_line_ptr) &
-      ,test_description_t('computing the gradient of a parabola within a tolerance of ' // string_t(loose_tolerance), check_grad_parabola_ptr) &
+       test_description_t('computing the 2nd-order 1D gradient of a constant within tolerance ' // string_t(tight_tolerance), check_grad_const_ptr) &
+      ,test_description_t('computing the 2nd-order 1D gradient of a line within tolerance ' // string_t(loose_tolerance), check_grad_line_ptr) &
+      ,test_description_t('computing the 2nd-order 1D gradient of a parabola within tolerance ' // string_t(loose_tolerance), check_grad_parabola_ptr) &
     ])
   end function
 
@@ -80,7 +80,7 @@ function check_grad_const() result(test_diagnosis)
     double precision, parameter :: df_dx = 0.
     procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => const
 
-    grad_f = .grad. cell_centers_extended_t(scalar_1D_initializer, order=2, cells=4, x_min=0D0, x_max=1D0)
+    grad_f = .grad. scalar_1D_t(scalar_1D_initializer, order=2, cells=4, x_min=0D0, x_max=1D0)
     test_diagnosis = .all. (grad_f%values() .approximates. df_dx .within. loose_tolerance) // " (d(line)/dx)"
   end function
 
@@ -96,7 +96,7 @@ function check_grad_line() result(test_diagnosis)
     double precision, parameter :: df_dx = 14D0
     procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => line
 
-    grad_f = .grad. cell_centers_extended_t(scalar_1D_initializer, order=2, cells=4, x_min=0D0, x_max=1D0)
+    grad_f = .grad. scalar_1D_t(scalar_1D_initializer, order=2, cells=4, x_min=0D0, x_max=1D0)
     test_diagnosis = .all. (grad_f%values() .approximates. df_dx .within. loose_tolerance) // " (d(line)/dx)"
   end function
 
@@ -108,11 +108,11 @@ pure function parabola(x) result(y)
 
   function check_grad_parabola() result(test_diagnosis)
     type(test_diagnosis_t) test_diagnosis
-    type(cell_centers_extended_t) quadratic
+    type(scalar_1D_t) quadratic
     type(gradient_t) grad_f
     procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => parabola
 
-    quadratic = cell_centers_extended_t(scalar_1D_initializer , order=2, cells=4, x_min=0D0, x_max=1D0)
+    quadratic = scalar_1D_t(scalar_1D_initializer , order=2, cells=4, x_min=0D0, x_max=1D0)
     grad_f = .grad. quadratic
 
     associate(x => grad_f%faces())

From 95054b4d65825cfdc081ad64cfcd7f937a1d2fe9 Mon Sep 17 00:00:00 2001
From: Damian Rouson <rouson@lbl.gov>
Date: Thu, 9 Oct 2025 00:53:43 -0400
Subject: [PATCH 16/40] refac(gradient_{t,m,s}):rename gradient_1D_{t,m,s}

---
 .../{gradient_m.f90 => gradient_1D_m.f90}     | 24 +++++++++----------
 .../{gradient_s.F90 => gradient_1D_s.F90}     | 12 +++++-----
 src/fortran/mole_m.f90                        |  2 +-
 src/fortran/scalar_1D_m.f90                   |  4 ++--
 src/fortran/scalar_1D_s.F90                   |  2 +-
 test/gradient_operator_test_m.F90             | 22 ++++++++---------
 6 files changed, 33 insertions(+), 33 deletions(-)
 rename src/fortran/{gradient_m.f90 => gradient_1D_m.f90} (62%)
 rename src/fortran/{gradient_s.F90 => gradient_1D_s.F90} (64%)

diff --git a/src/fortran/gradient_m.f90 b/src/fortran/gradient_1D_m.f90
similarity index 62%
rename from src/fortran/gradient_m.f90
rename to src/fortran/gradient_1D_m.f90
index 3c669a9b..62ed6a80 100644
--- a/src/fortran/gradient_m.f90
+++ b/src/fortran/gradient_1D_m.f90
@@ -1,15 +1,15 @@
-module gradient_m
+module gradient_1D_m
   !! Define an abstraction for the collection of points used to compute gradidents:
   !! cell centers plus oundaries.
   implicit none
 
   private
-  public :: gradient_t
+  public :: gradient_1D_t
 
-  type gradient_t
-    !! Encapsulate gradient values produced only by .grad. (no other constructors)
+  type gradient_1D_t
+    !! Encapsulate gradient_1D values produced only by .grad. (no other constructors)
     private
-    double precision, allocatable :: vector_1D_(:) !! gradient values at cell faces (nodes in 1D)
+    double precision, allocatable :: vector_1D_(:) !! gradient_1D values at cell faces (nodes in 1D)
     double precision x_min_ !! domain lower boundary
     double precision x_max_ !! domain upper boundary
     integer cells_ !! number of grid cells spanning the domain
@@ -18,14 +18,14 @@ module gradient_m
     procedure faces
   end type
 
-  interface gradient_t
+  interface gradient_1D_t
 
-      pure module function construct_from_components(face_centered_values, x_min, x_max, cells) result(gradient)
-        !! Result is an object storing gradients at cell faces
+      pure module function construct_from_components(face_centered_values, x_min, x_max, cells) result(gradient_1D)
+        !! Result is an object storing gradient_1Ds at cell faces
         implicit none
         double precision, intent(in) :: face_centered_values(:), x_min, x_max
         integer, intent(in) :: cells
-        type(gradient_t) gradient
+        type(gradient_1D_t) gradient_1D
       end function
 
   end interface
@@ -34,16 +34,16 @@ pure module function construct_from_components(face_centered_values, x_min, x_ma
 
     pure module function faces(self) result(x)
       implicit none
-      class(gradient_t), intent(in) :: self
+      class(gradient_1D_t), intent(in) :: self
       double precision, allocatable :: x(:)
     end function
 
      pure module function values(self) result(gradients)
        implicit none
-       class(gradient_t), intent(in) :: self
+       class(gradient_1D_t), intent(in) :: self
        double precision, allocatable :: gradients(:)
      end function
 
   end interface
 
-end module gradient_m
\ No newline at end of file
+end module gradient_1D_m
\ No newline at end of file
diff --git a/src/fortran/gradient_s.F90 b/src/fortran/gradient_1D_s.F90
similarity index 64%
rename from src/fortran/gradient_s.F90
rename to src/fortran/gradient_1D_s.F90
index 125262c9..9cff5f90 100644
--- a/src/fortran/gradient_s.F90
+++ b/src/fortran/gradient_1D_s.F90
@@ -1,13 +1,13 @@
-submodule(gradient_m) gradient_s
+submodule(gradient_1D_m) gradient_1D_s
   implicit none
 
 contains
 
   module procedure construct_from_components
-    gradient%vector_1D_ = face_centered_values
-    gradient%x_min_ = x_min
-    gradient%x_max_ = x_max
-    gradient%cells_ = cells
+    gradient_1D%vector_1D_ = face_centered_values
+    gradient_1D%x_min_ = x_min
+    gradient_1D%x_max_ = x_max
+    gradient_1D%cells_ = cells
   end procedure
 
   module procedure values
@@ -22,4 +22,4 @@
         ]
   end procedure
 
-end submodule gradient_s
\ No newline at end of file
+end submodule gradient_1D_s
\ No newline at end of file
diff --git a/src/fortran/mole_m.f90 b/src/fortran/mole_m.f90
index 672a2a94..06fbda92 100644
--- a/src/fortran/mole_m.f90
+++ b/src/fortran/mole_m.f90
@@ -1,6 +1,6 @@
 module mole_m 
   !! MOLE Fortran public entities
   use scalar_1D_m, only : scalar_1D_t, scalar_1D_initializer_i
-  use gradient_m, only : gradient_t
+  use gradient_1D_m, only : gradient_1D_t
   implicit none
 end module mole_m
diff --git a/src/fortran/scalar_1D_m.f90 b/src/fortran/scalar_1D_m.f90
index 21f4d0b5..a3c037e7 100644
--- a/src/fortran/scalar_1D_m.f90
+++ b/src/fortran/scalar_1D_m.f90
@@ -1,7 +1,7 @@
 module scalar_1D_m
   !! Define an abstraction for the collection of points used to compute gradidents:
   !! cell centers plus oundaries.
-  use gradient_m, only : gradient_t
+  use gradient_1D_m, only : gradient_1D_t
   implicit none
 
   private
@@ -75,7 +75,7 @@ pure module function grad(self) result(grad_f)
       !! Result is mimetic gradient of f
       implicit none
       class(scalar_1D_t), intent(in) :: self
-      type(gradient_t) grad_f !! discrete gradient approximation
+      type(gradient_1D_t) grad_f !! discrete gradient approximation
     end function
 
   end interface
diff --git a/src/fortran/scalar_1D_s.F90 b/src/fortran/scalar_1D_s.F90
index 859feed6..66c3ce20 100644
--- a/src/fortran/scalar_1D_s.F90
+++ b/src/fortran/scalar_1D_s.F90
@@ -33,7 +33,7 @@ pure function grid_(x_min, x_max, cells) result(x)
   end procedure
 
   module procedure grad
-    grad_f = gradient_t(matvec(self%gradient_operator_%mimetic_matrix_, self), self%x_min_, self%x_max_, self%cells_)
+    grad_f = gradient_1D_t(matvec(self%gradient_operator_%mimetic_matrix_, self), self%x_min_, self%x_max_, self%cells_)
   end procedure
 
 end submodule scalar_1D_s
\ No newline at end of file
diff --git a/test/gradient_operator_test_m.F90 b/test/gradient_operator_test_m.F90
index f0674a84..7a44dc72 100644
--- a/test/gradient_operator_test_m.F90
+++ b/test/gradient_operator_test_m.F90
@@ -15,7 +15,7 @@ module gradient_operator_test_m
     ,operator(.all.) &
     ,operator(.approximates.) &
     ,operator(.within.)
-  use mole_m, only : scalar_1D_t, gradient_t, scalar_1D_initializer_i
+  use mole_m, only : scalar_1D_t, gradient_1D_t, scalar_1D_initializer_i
 #if ! HAVE_PROCEDURE_ACTUAL_FOR_POINTER_DUMMY
   use julienne_m, only : diagnosis_function_i
 #endif
@@ -42,9 +42,9 @@ function results() result(test_results)
     type(gradient_operator_test_t) gradient_operator_test
     type(test_result_t), allocatable :: test_results(:)
     test_results = gradient_operator_test%run([ & 
-       test_description_t('computing the 2nd-order 1D gradient of a constant within tolerance ' // string_t(tight_tolerance), check_grad_const) &
-      ,test_description_t('computing the 2nd-order 1D gradient of a line within tolerance ' // string_t(loose_tolerance), check_grad_line) &
-      ,test_description_t('computing the 2nd-order 1D gradient of a parabola within tolerance ' // string_t(loose_tolerance), check_grad_parabola) &
+       test_description_t('computing the 2nd-order 1D gradient of a constant within a tolerance of ' // string_t(tight_tolerance), check_grad_const) &
+      ,test_description_t('computing the 2nd-order 1D gradient of a line within a tolerance of ' // string_t(loose_tolerance), check_grad_line) &
+      ,test_description_t('computing the 2nd-order 1D gradient of a parabola within a tolerance of ' // string_t(loose_tolerance), check_grad_parabola) &
     ])
   end function
 
@@ -55,13 +55,13 @@ function results() result(test_results)
     type(test_result_t), allocatable :: test_results(:)
     procedure(diagnosis_function_i), pointer :: &
        check_grad_const_ptr => check_grad_const &
-      ,check_grad_line_ptr => check_grad_line
+      ,check_grad_line_ptr => check_grad_line &
       ,check_grad_parabola_ptr => check_grad_parabola
 
     test_results = gradient_operator_test%run([ &
-       test_description_t('computing the 2nd-order 1D gradient of a constant within tolerance ' // string_t(tight_tolerance), check_grad_const_ptr) &
-      ,test_description_t('computing the 2nd-order 1D gradient of a line within tolerance ' // string_t(loose_tolerance), check_grad_line_ptr) &
-      ,test_description_t('computing the 2nd-order 1D gradient of a parabola within tolerance ' // string_t(loose_tolerance), check_grad_parabola_ptr) &
+       test_description_t('computing the 2nd-order 1D gradient of a constant within a tolerance of ' // string_t(tight_tolerance), check_grad_const_ptr) &
+      ,test_description_t('computing the 2nd-order 1D gradient of a line within a tolerance of ' // string_t(loose_tolerance), check_grad_line_ptr) &
+      ,test_description_t('computing the 2nd-order 1D gradient of a parabola within a tolerance of ' // string_t(loose_tolerance), check_grad_parabola_ptr) &
     ])
   end function
 
@@ -76,7 +76,7 @@ pure function const(x) result(y)
 
   function check_grad_const() result(test_diagnosis)
     type(test_diagnosis_t) test_diagnosis
-    type(gradient_t) grad_f
+    type(gradient_1D_t) grad_f
     double precision, parameter :: df_dx = 0.
     procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => const
 
@@ -92,7 +92,7 @@ pure function line(x) result(y)
 
   function check_grad_line() result(test_diagnosis)
     type(test_diagnosis_t) test_diagnosis
-    type(gradient_t) grad_f
+    type(gradient_1D_t) grad_f
     double precision, parameter :: df_dx = 14D0
     procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => line
 
@@ -109,7 +109,7 @@ pure function parabola(x) result(y)
   function check_grad_parabola() result(test_diagnosis)
     type(test_diagnosis_t) test_diagnosis
     type(scalar_1D_t) quadratic
-    type(gradient_t) grad_f
+    type(gradient_1D_t) grad_f
     procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => parabola
 
     quadratic = scalar_1D_t(scalar_1D_initializer , order=2, cells=4, x_min=0D0, x_max=1D0)

From 80451cfda706d5d3a489612d95bc0a118209d951 Mon Sep 17 00:00:00 2001
From: Damian Rouson <rouson@lbl.gov>
Date: Thu, 9 Oct 2025 01:13:37 -0400
Subject: [PATCH 17/40] refac(gradient_operator):name gradient_operator_1D

---
 ..._operator_s.F90 => gradient_operator_1D_s.F90} | 12 ++++++------
 src/fortran/scalar_1D_m.f90                       | 12 ++++++------
 src/fortran/scalar_1D_s.F90                       |  4 ++--
 test/driver.f90                                   |  6 ++----
 ...test_m.F90 => gradient_operator_1D_test_m.F90} | 15 ++++++---------
 5 files changed, 22 insertions(+), 27 deletions(-)
 rename src/fortran/{gradient_operator_s.F90 => gradient_operator_1D_s.F90} (91%)
 rename test/{gradient_operator_test_m.F90 => gradient_operator_1D_test_m.F90} (89%)

diff --git a/src/fortran/gradient_operator_s.F90 b/src/fortran/gradient_operator_1D_s.F90
similarity index 91%
rename from src/fortran/gradient_operator_s.F90
rename to src/fortran/gradient_operator_1D_s.F90
index e649e1ff..c101fde6 100644
--- a/src/fortran/gradient_operator_s.F90
+++ b/src/fortran/gradient_operator_1D_s.F90
@@ -1,7 +1,7 @@
 #include "julienne-assert-macros.h"
 #include "mole-language-support.F90"
 
-submodule(scalar_1D_m) gradient_operator_s
+submodule(scalar_1D_m) gradient_operator_1D_s
   use julienne_m, only : call_julienne_assert_, string_t
 #if ASSERTIONS
   use julienne_m, only : operator(.isAtLeast.)
@@ -14,14 +14,14 @@
 
     call_julienne_assert(m .isAtLeast. 2*k)
 
-    gradient_operator%mimetic_matrix_ = mimetic_matrix_t( &
+    gradient_operator_1D%mimetic_matrix_ = mimetic_matrix_t( &
        corbino_castillo_A( k, dx) & 
       ,corbino_castillo_M( k, dx) & 
       ,corbino_castillo_Ap(k, dx) &
     )
-    gradient_operator%k_  = k
-    gradient_operator%dx_ = dx
-    gradient_operator%m_ = m
+    gradient_operator_1D%k_  = k
+    gradient_operator_1D%dx_ = dx
+    gradient_operator_1D%m_ = m
   end procedure
 
   pure function corbino_castillo_A(k, dx) result(rows)
@@ -101,4 +101,4 @@ pure function corbino_castillo_Ap(k, dx) result(rows)
 
 #endif
 
-end submodule gradient_operator_s
\ No newline at end of file
+end submodule gradient_operator_1D_s
\ No newline at end of file
diff --git a/src/fortran/scalar_1D_m.f90 b/src/fortran/scalar_1D_m.f90
index a3c037e7..40a80104 100644
--- a/src/fortran/scalar_1D_m.f90
+++ b/src/fortran/scalar_1D_m.f90
@@ -6,7 +6,7 @@ module scalar_1D_m
 
   private
   public :: scalar_1D_t
-  public :: gradient_operator_t
+  public :: gradient_operator_1D_t
   public :: scalar_1D_initializer_i
 
   abstract interface
@@ -25,7 +25,7 @@ pure function scalar_1D_initializer_i(x) result(f)
     double precision, allocatable :: upper_(:,:), inner_(:), lower_(:,:)
   end type
 
-  type gradient_operator_t
+  type gradient_operator_1D_t
     !! Encapsulate kth-order mimetic gradient operator on dx-sized cells
     private
     integer k_, m_
@@ -39,7 +39,7 @@ pure function scalar_1D_initializer_i(x) result(f)
     double precision, allocatable :: scalar_1D_(:)
     double precision x_min_, x_max_
     integer cells_
-    type(gradient_operator_t) gradient_operator_
+    type(gradient_operator_1D_t) gradient_operator_1D_
   contains
     procedure grid
     generic :: operator(.grad.) => grad
@@ -80,15 +80,15 @@ pure module function grad(self) result(grad_f)
 
   end interface
 
-  interface gradient_operator_t
+  interface gradient_operator_1D_t
 
-    pure module function construct_from_parameters(k, dx, m) result(gradient_operator)
+    pure module function construct_from_parameters(k, dx, m) result(gradient_operator_1D)
       !! Construct a mimetic gradient operator
       implicit none
       integer, intent(in) :: k !! order of accuracy
       double precision, intent(in) :: dx !! step siz
       integer, intent(in) :: m !! number of grid cells
-      type(gradient_operator_t) gradient_operator
+      type(gradient_operator_1D_t) gradient_operator_1D
     end function
 
   end interface
diff --git a/src/fortran/scalar_1D_s.F90 b/src/fortran/scalar_1D_s.F90
index 66c3ce20..633dcfdf 100644
--- a/src/fortran/scalar_1D_s.F90
+++ b/src/fortran/scalar_1D_s.F90
@@ -13,7 +13,7 @@
     scalar_1D%x_min_ = x_min
     scalar_1D%x_max_ = x_max
     scalar_1D%cells_ = cells
-    scalar_1D%gradient_operator_ = gradient_operator_t(k=order, dx=(x_max - x_min)/cells, m=cells)
+    scalar_1D%gradient_operator_1D_ = gradient_operator_1D_t(k=order, dx=(x_max - x_min)/cells, m=cells)
     scalar_1D%scalar_1D_ = initializer(grid_(x_min, x_max, cells))
   end procedure
 
@@ -33,7 +33,7 @@ pure function grid_(x_min, x_max, cells) result(x)
   end procedure
 
   module procedure grad
-    grad_f = gradient_1D_t(matvec(self%gradient_operator_%mimetic_matrix_, self), self%x_min_, self%x_max_, self%cells_)
+    grad_f = gradient_1D_t(matvec(self%gradient_operator_1D_%mimetic_matrix_, self), self%x_min_, self%x_max_, self%cells_)
   end procedure
 
 end submodule scalar_1D_s
\ No newline at end of file
diff --git a/test/driver.f90 b/test/driver.f90
index 703733f4..724b7fca 100644
--- a/test/driver.f90
+++ b/test/driver.f90
@@ -1,12 +1,10 @@
-! Copyright (c) 2024-2025, The Regents of the University of California and Sourcery Institute
-! Terms of use are as specified in https://github.com/BerkeleyLab/julienne/blob/3.1.2/LICENSE.txt
 program test_suite_driver
   use julienne_m, only : test_fixture_t, test_harness_t
-  use gradient_operator_test_m, only : gradient_operator_test_t
+  use gradient_operator_1D_test_m, only : gradient_operator_1D_test_t
   implicit none
 
   associate(test_harness => test_harness_t([ &
-     test_fixture_t(gradient_operator_test_t()) &
+     test_fixture_t(gradient_operator_1D_test_t()) &
   ]))
     call test_harness%report_results
   end associate
diff --git a/test/gradient_operator_test_m.F90 b/test/gradient_operator_1D_test_m.F90
similarity index 89%
rename from test/gradient_operator_test_m.F90
rename to test/gradient_operator_1D_test_m.F90
index 7a44dc72..64cba1e8 100644
--- a/test/gradient_operator_test_m.F90
+++ b/test/gradient_operator_1D_test_m.F90
@@ -1,10 +1,7 @@
-! Copyright (c) 2024-2025, The Regents of the University of California and Sourcery Institute
-! Terms of use are as specified in https://github.com/BerkeleyLab/julienne/blob/3.1.2/LICENSE.txt
-
 #include "language-support.F90"
   !! include Julienne preprocessor  macros
 
-module gradient_operator_test_m
+module gradient_operator_1D_test_m
   use julienne_m, only : &
      string_t &
     ,test_t &
@@ -21,7 +18,7 @@ module gradient_operator_test_m
 #endif
   implicit none
 
-  type, extends(test_t) :: gradient_operator_test_t
+  type, extends(test_t) :: gradient_operator_1D_test_t
   contains
     procedure, nopass :: subject
     procedure, nopass :: results
@@ -39,9 +36,9 @@ pure function subject() result(test_subject)
 #if HAVE_PROCEDURE_ACTUAL_FOR_POINTER_DUMMY
 
   function results() result(test_results)
-    type(gradient_operator_test_t) gradient_operator_test
+    type(gradient_operator_1D_test_t) gradient_operator_1D_test
     type(test_result_t), allocatable :: test_results(:)
-    test_results = gradient_operator_test%run([ & 
+    test_results = gradient_operator_1D_test%run([ & 
        test_description_t('computing the 2nd-order 1D gradient of a constant within a tolerance of ' // string_t(tight_tolerance), check_grad_const) &
       ,test_description_t('computing the 2nd-order 1D gradient of a line within a tolerance of ' // string_t(loose_tolerance), check_grad_line) &
       ,test_description_t('computing the 2nd-order 1D gradient of a parabola within a tolerance of ' // string_t(loose_tolerance), check_grad_parabola) &
@@ -51,14 +48,14 @@ function results() result(test_results)
 #else
 
   function results() result(test_results)
-    type(gradient_operator_test_t) gradient_operator_test
+    type(gradient_operator_1D_test_t) gradient_operator_1D_test
     type(test_result_t), allocatable :: test_results(:)
     procedure(diagnosis_function_i), pointer :: &
        check_grad_const_ptr => check_grad_const &
       ,check_grad_line_ptr => check_grad_line &
       ,check_grad_parabola_ptr => check_grad_parabola
 
-    test_results = gradient_operator_test%run([ &
+    test_results = gradient_operator_1D_test%run([ &
        test_description_t('computing the 2nd-order 1D gradient of a constant within a tolerance of ' // string_t(tight_tolerance), check_grad_const_ptr) &
       ,test_description_t('computing the 2nd-order 1D gradient of a line within a tolerance of ' // string_t(loose_tolerance), check_grad_line_ptr) &
       ,test_description_t('computing the 2nd-order 1D gradient of a parabola within a tolerance of ' // string_t(loose_tolerance), check_grad_parabola_ptr) &

From 46e074d133ff1023c2778fc59ddd57f7d6e3af04 Mon Sep 17 00:00:00 2001
From: Damian Rouson <rouson@lbl.gov>
Date: Thu, 9 Oct 2025 01:17:44 -0400
Subject: [PATCH 18/40] refac(mimetic_matrix}: append "_1D"

---
 src/fortran/gradient_operator_1D_s.F90               |  2 +-
 ...{mimetic_matrix_s.F90 => mimetic_matrix_1D_s.F90} | 10 +++++-----
 src/fortran/scalar_1D_m.f90                          | 12 ++++++------
 src/fortran/scalar_1D_s.F90                          |  2 +-
 4 files changed, 13 insertions(+), 13 deletions(-)
 rename src/fortran/{mimetic_matrix_s.F90 => mimetic_matrix_1D_s.F90} (85%)

diff --git a/src/fortran/gradient_operator_1D_s.F90 b/src/fortran/gradient_operator_1D_s.F90
index c101fde6..d74a9d18 100644
--- a/src/fortran/gradient_operator_1D_s.F90
+++ b/src/fortran/gradient_operator_1D_s.F90
@@ -14,7 +14,7 @@
 
     call_julienne_assert(m .isAtLeast. 2*k)
 
-    gradient_operator_1D%mimetic_matrix_ = mimetic_matrix_t( &
+    gradient_operator_1D%mimetic_matrix_1D_ = mimetic_matrix_1D_t( &
        corbino_castillo_A( k, dx) & 
       ,corbino_castillo_M( k, dx) & 
       ,corbino_castillo_Ap(k, dx) &
diff --git a/src/fortran/mimetic_matrix_s.F90 b/src/fortran/mimetic_matrix_1D_s.F90
similarity index 85%
rename from src/fortran/mimetic_matrix_s.F90
rename to src/fortran/mimetic_matrix_1D_s.F90
index 7a6ec067..aeb0997a 100644
--- a/src/fortran/mimetic_matrix_s.F90
+++ b/src/fortran/mimetic_matrix_1D_s.F90
@@ -1,16 +1,16 @@
 #include "mole-language-support.F90"
 #include "julienne-assert-macros.h"
 
-submodule(scalar_1D_m) mimetic_matrix_s
+submodule(scalar_1D_m) mimetic_matrix_1D_s
   use julienne_m, only : call_julienne_assert_, string_t
   implicit none
 
 contains
 
   module procedure construct_from_components
-    mimetic_matrix%upper_ = upper
-    mimetic_matrix%inner_ = inner
-    mimetic_matrix%lower_ = lower
+    mimetic_matrix_1D%upper_ = upper
+    mimetic_matrix_1D%inner_ = inner
+    mimetic_matrix_1D%lower_ = lower
   end procedure
 
 #if HAVE_DO_CONCURRENT_TYPE_SPEC_SUPPORT
@@ -41,4 +41,4 @@
 
 #endif
 
-end submodule mimetic_matrix_s
\ No newline at end of file
+end submodule mimetic_matrix_1D_s
\ No newline at end of file
diff --git a/src/fortran/scalar_1D_m.f90 b/src/fortran/scalar_1D_m.f90
index 40a80104..2e41979d 100644
--- a/src/fortran/scalar_1D_m.f90
+++ b/src/fortran/scalar_1D_m.f90
@@ -19,7 +19,7 @@ pure function scalar_1D_initializer_i(x) result(f)
 
   end interface
 
-  type mimetic_matrix_t
+  type mimetic_matrix_1D_t
     !! Encapsulate a mimetic matrix with a corresponding matrix-vector product operator
     private
     double precision, allocatable :: upper_(:,:), inner_(:), lower_(:,:)
@@ -30,7 +30,7 @@ pure function scalar_1D_initializer_i(x) result(f)
     private
     integer k_, m_
     double precision dx_
-    type(mimetic_matrix_t) mimetic_matrix_
+    type(mimetic_matrix_1D_t) mimetic_matrix_1D_
   end type
 
   type scalar_1D_t
@@ -93,13 +93,13 @@ pure module function construct_from_parameters(k, dx, m) result(gradient_operato
 
   end interface
 
-  interface mimetic_matrix_t
+  interface mimetic_matrix_1D_t
 
-    pure module function construct_from_components(upper, inner, lower) result(mimetic_matrix)
+    pure module function construct_from_components(upper, inner, lower) result(mimetic_matrix_1D)
       !! Construct discrete operator from coefficient matrix
       implicit none
       double precision, intent(in) :: upper(:,:), inner(:), lower(:,:)
-      type(mimetic_matrix_t) mimetic_matrix
+      type(mimetic_matrix_1D_t) mimetic_matrix_1D
     end function
 
   end interface
@@ -109,7 +109,7 @@ pure module function construct_from_components(upper, inner, lower) result(mimet
     pure module function matvec(self, vector) result(matvec_product)
       !! Apply a matrix operator to a vector
       implicit none
-      class(mimetic_matrix_t), intent(in) :: self
+      class(mimetic_matrix_1D_t), intent(in) :: self
       type(scalar_1D_t), intent(in) :: vector
       double precision, allocatable :: matvec_product(:)
     end function
diff --git a/src/fortran/scalar_1D_s.F90 b/src/fortran/scalar_1D_s.F90
index 633dcfdf..df1672ab 100644
--- a/src/fortran/scalar_1D_s.F90
+++ b/src/fortran/scalar_1D_s.F90
@@ -33,7 +33,7 @@ pure function grid_(x_min, x_max, cells) result(x)
   end procedure
 
   module procedure grad
-    grad_f = gradient_1D_t(matvec(self%gradient_operator_1D_%mimetic_matrix_, self), self%x_min_, self%x_max_, self%cells_)
+    grad_f = gradient_1D_t(matvec(self%gradient_operator_1D_%mimetic_matrix_1D_, self), self%x_min_, self%x_max_, self%cells_)
   end procedure
 
 end submodule scalar_1D_s
\ No newline at end of file

From 17261500fee6e51a9b159d13fdcdd189a2908a55 Mon Sep 17 00:00:00 2001
From: Damian Rouson <rouson@lbl.gov>
Date: Thu, 9 Oct 2025 01:31:42 -0400
Subject: [PATCH 19/40] fix(scalar_1D_s): import operators

---
 src/fortran/scalar_1D_s.F90 | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/fortran/scalar_1D_s.F90 b/src/fortran/scalar_1D_s.F90
index df1672ab..57425e9a 100644
--- a/src/fortran/scalar_1D_s.F90
+++ b/src/fortran/scalar_1D_s.F90
@@ -1,13 +1,13 @@
 #include "julienne-assert-macros.h"
 
 submodule(scalar_1D_m) scalar_1D_s
-  use julienne_m, only : call_julienne_assert_, operator(.equalsExpected.)
+  use julienne_m, only : call_julienne_assert_, operator(.equalsExpected.), operator(.greaterThan.), operator(.isAtLeast.)
   implicit none
 
 contains
 
   module procedure construct_from_function
-    call_julienne_assert(x_max .isGreaterThan. x_min)
+    call_julienne_assert(x_max .greaterThan. x_min)
     call_julienne_assert(cells .isAtLeast. 2*order)
 
     scalar_1D%x_min_ = x_min

From b9c69d6dbc2a4257b834ae2429ed47af1f97b4a0 Mon Sep 17 00:00:00 2001
From: Damian Rouson <rouson@lbl.gov>
Date: Thu, 9 Oct 2025 01:36:06 -0400
Subject: [PATCH 20/40] fix(mimetic_matrix_s): count upper/lower rows

---
 src/fortran/mimetic_matrix_1D_s.F90 | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/fortran/mimetic_matrix_1D_s.F90 b/src/fortran/mimetic_matrix_1D_s.F90
index aeb0997a..b76db466 100644
--- a/src/fortran/mimetic_matrix_1D_s.F90
+++ b/src/fortran/mimetic_matrix_1D_s.F90
@@ -2,7 +2,7 @@
 #include "julienne-assert-macros.h"
 
 submodule(scalar_1D_m) mimetic_matrix_1D_s
-  use julienne_m, only : call_julienne_assert_, string_t
+  use julienne_m, only : call_julienne_assert_, string_t, operator(.equalsExpected.)
   implicit none
 
 contains
@@ -19,7 +19,7 @@
 
     double precision, allocatable :: product_inner(:)
 
-    associate(upper => merge(0, 1, size(self%upper_)==0), lower => merge(0, 1, size(self%lower_)==0))
+    associate(upper => size(self%upper_,1), lower => size(self%lower_,1))
       associate(inner_rows => size(vector%scalar_1D_) - (upper + lower + 1))
 
         allocate(product_inner(inner_rows))

From ac815c25893a53af3a5acddc5df57791f5ec5f58 Mon Sep 17 00:00:00 2001
From: Damian Rouson <rouson@lbl.gov>
Date: Thu, 9 Oct 2025 01:38:50 -0400
Subject: [PATCH 21/40] refac(mimetic_matrix_s): simpler loop bound

---
 src/fortran/mimetic_matrix_1D_s.F90 | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/fortran/mimetic_matrix_1D_s.F90 b/src/fortran/mimetic_matrix_1D_s.F90
index b76db466..ae9888fc 100644
--- a/src/fortran/mimetic_matrix_1D_s.F90
+++ b/src/fortran/mimetic_matrix_1D_s.F90
@@ -24,7 +24,7 @@
 
         allocate(product_inner(inner_rows))
 
-        do concurrent(integer :: row = 1 : size(product_inner)) default(none) shared(product_inner, self, vector)
+        do concurrent(integer :: row = 1 : inner_rows) default(none) shared(product_inner, self, vector)
           product_inner(row) = dot_product(self%inner_, vector%scalar_1D_(row + 1 : row + size(self%inner_)))
         end do
 

From 33bcfb80dadfa568e79d7ae3e25932f8985ed384 Mon Sep 17 00:00:00 2001
From: Damian Rouson <rouson@lbl.gov>
Date: Thu, 9 Oct 2025 01:44:58 -0400
Subject: [PATCH 22/40] fix(4th-order): flip some coefficient signs

---
 src/fortran/gradient_operator_1D_s.F90 | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/fortran/gradient_operator_1D_s.F90 b/src/fortran/gradient_operator_1D_s.F90
index d74a9d18..c9886506 100644
--- a/src/fortran/gradient_operator_1D_s.F90
+++ b/src/fortran/gradient_operator_1D_s.F90
@@ -35,8 +35,8 @@ pure function corbino_castillo_A(k, dx) result(rows)
       rows = reshape([-8D0/3D0, 3D0, -1D0/3D0] , shape=[1,3]) / dx
     case(4)
       rows = transpose(reshape([ & 
-         [-352D0/105D0,  35D0/ 8D0, -35D0/24D0, -21D0/40D0, -5D0/ 56D0] &
-        ,[ -16D0/105D0, -31D0/24D0,  29D0/24D0,  -3D0/40D0,  1D0/168D0] &
+         [-352D0/105D0,  35D0/ 8D0, -35D0/24D0, 21D0/40D0, -5D0/ 56D0] &
+        ,[  16D0/105D0, -31D0/24D0,  29D0/24D0, -3D0/40D0,  1D0/168D0] &
       ], shape=[5,2]))
     case default
       associate(string_k => string_t(k))

From 623342906be01921afb7881105e60e99c17993df Mon Sep 17 00:00:00 2001
From: Damian Rouson <rouson@lbl.gov>
Date: Thu, 9 Oct 2025 17:14:26 -0700
Subject: [PATCH 23/40] feat(mimetic_matrix_1D): add file_t constructor

This commit adds the ability to convert a mimetic matrix to a
Julienne file_t object so that the matrix elements can be printed
with code like the following

  type(mimetic_matrix_t) mimetic_matrix
  ! define matrix here
  associate(file => mimetic_matrix%to_file_t())
    call file%write_lines()
  end associate
---
 src/fortran/mimetic_matrix_1D_s.F90 | 22 +++++++++++++++++++++-
 1 file changed, 21 insertions(+), 1 deletion(-)

diff --git a/src/fortran/mimetic_matrix_1D_s.F90 b/src/fortran/mimetic_matrix_1D_s.F90
index ae9888fc..aa900d26 100644
--- a/src/fortran/mimetic_matrix_1D_s.F90
+++ b/src/fortran/mimetic_matrix_1D_s.F90
@@ -2,7 +2,7 @@
 #include "julienne-assert-macros.h"
 
 submodule(scalar_1D_m) mimetic_matrix_1D_s
-  use julienne_m, only : call_julienne_assert_, string_t, operator(.equalsExpected.)
+  use julienne_m, only : call_julienne_assert_, string_t, operator(.equalsExpected.), operator(.csv.)
   implicit none
 
 contains
@@ -41,4 +41,24 @@
 
 #endif
 
+  module procedure to_file_t
+    type(string_t), allocatable :: lines(:)
+    integer, parameter :: inner_rows = 1
+    integer row
+
+    associate(upper_rows => size(self%upper_,1), lower_rows => size(self%lower_,1))
+      allocate(lines(upper_rows + inner_rows + lower_rows))
+      do row = 1, upper_rows
+        lines(row) = .csv. string_t(self%upper_(row,:))
+      end do
+      lines(upper_rows + inner_rows) = .csv. string_t(self%inner_)
+      do row = 1, lower_rows
+        lines(upper_rows + inner_rows + row) = .csv. string_t(self%lower_(row,:))
+      end do
+    end associate
+
+    file = file_t(lines)
+
+  end procedure
+
 end submodule mimetic_matrix_1D_s
\ No newline at end of file

From c5ad4283d2123e7c739fa2010915f5da7b565a4e Mon Sep 17 00:00:00 2001
From: Damian Rouson <rouson@lbl.gov>
Date: Thu, 9 Oct 2025 17:18:51 -0700
Subject: [PATCH 24/40] fix(mimetic_matrix_t): reshape upper block (A)

Now all tests pass for computing 1D gradients of 0th, 1st, &
2nd-order polynomials, each with mimetic discretizations of
2nd- and 4th-order accuracy.
---
 src/fortran/gradient_operator_1D_s.F90 | 47 ++++++++++----------------
 src/fortran/scalar_1D_m.f90            | 13 +++++++
 test/gradient_operator_1D_test_m.F90   | 41 +++++++++++++++-------
 3 files changed, 60 insertions(+), 41 deletions(-)

diff --git a/src/fortran/gradient_operator_1D_s.F90 b/src/fortran/gradient_operator_1D_s.F90
index c9886506..ef9eb275 100644
--- a/src/fortran/gradient_operator_1D_s.F90
+++ b/src/fortran/gradient_operator_1D_s.F90
@@ -24,20 +24,20 @@
     gradient_operator_1D%m_ = m
   end procedure
 
-  pure function corbino_castillo_A(k, dx) result(rows)
+  pure function corbino_castillo_A(k, dx) result(matrix_block)
     integer, intent(in) :: k
     double precision, intent(in) :: dx
-    double precision, allocatable :: rows(:,:)
+    double precision, allocatable :: matrix_block(:,:)
 
     order_of_accuracy: &
     select case(k)
     case(2)
-      rows = reshape([-8D0/3D0, 3D0, -1D0/3D0] , shape=[1,3]) / dx
+      matrix_block = reshape([-8D0/3D0, 3D0, -1D0/3D0] , shape=[1,3]) / dx
     case(4)
-      rows = transpose(reshape([ & 
-         [-352D0/105D0,  35D0/ 8D0, -35D0/24D0, 21D0/40D0, -5D0/ 56D0] &
-        ,[  16D0/105D0, -31D0/24D0,  29D0/24D0, -3D0/40D0,  1D0/168D0] &
-      ], shape=[5,2]))
+      matrix_block = reshape([ &
+         -352D0/105D0,  35D0/ 8D0, -35D0/24D0, 21D0/40D0, -5D0/ 56D0 &
+        ,  16D0/105D0, -31D0/24D0,  29D0/24D0, -3D0/40D0,  1D0/168D0 &
+      ], shape=[2,5], order=[2,1]) / dx
     case default
       associate(string_k => string_t(k))
         error stop "corbino_castillo_A: unsupported order of accuracy: " // string_k%string()
@@ -68,37 +68,26 @@ pure function corbino_castillo_M(k, dx) result(row)
 #if HAVE_DO_CONCURRENT_TYPE_SPEC_SUPPORT
 
 
-  pure function corbino_castillo_Ap(k, dx) result(rows)
+  pure function corbino_castillo_Ap(k, dx) result(matrix_block)
     integer, intent(in) :: k
     double precision, intent(in) :: dx
-    double precision, allocatable :: rows(:,:)
+    double precision, allocatable :: matrix_block(:,:)
 
     associate(A => corbino_castillo_A(k, dx))
-      allocate(rows , mold=A)
-      reverse_and_flip_sign: &
-      do concurrent(integer :: row = 1:size(rows,1)) default(none) shared(rows, A)
-        rows(row,:) = -A(row,size(A,2):1:-1)
-      end do reverse_and_flip_sign
+      allocate(matrix_block , mold=A)
+      reverse_elements_within_matrix_block_and_flip_sign: &
+      do concurrent(integer :: row = 1:size(matrix_block,1)) default(none) shared(matrix_block, A)
+        matrix_block(row,:) = -A(row,size(A,2):1:-1)
+      end do reverse_elements_within_matrix_block_and_flip_sign
+      reverse_elements_within_columns: &
+      do concurrent(integer :: column = 1 : size(matrix_block,2)) default(none) shared(matrix_block)
+        matrix_block(:,column) = matrix_block(size(matrix_block,1):1:-1,column)
+      end do reverse_elements_within_columns
     end associate
   end function
  
 #else
 
-  pure function corbino_castillo_Ap(k, dx) result(rows)
-    integer, intent(in) :: k
-    double precision, intent(in) :: dx
-    double precision, allocatable :: rows(:,:)
-    integer row
-
-    associate(A => corbino_castillo_A(k, dx))
-      allocate(rows , mold=A)
-      reverse_and_flip_sign: &
-      do concurrent(row = 1:size(rows,1)) default(none) shared(rows, A)
-        rows(row,:) = -A(row,size(A,2):1:-1)
-      end do reverse_and_flip_sign
-    end associate
-  end function
-
 #endif
 
 end submodule gradient_operator_1D_s
\ No newline at end of file
diff --git a/src/fortran/scalar_1D_m.f90 b/src/fortran/scalar_1D_m.f90
index 2e41979d..fca7b831 100644
--- a/src/fortran/scalar_1D_m.f90
+++ b/src/fortran/scalar_1D_m.f90
@@ -1,6 +1,7 @@
 module scalar_1D_m
   !! Define an abstraction for the collection of points used to compute gradidents:
   !! cell centers plus oundaries.
+  use julienne_m, only : file_t
   use gradient_1D_m, only : gradient_1D_t
   implicit none
 
@@ -23,6 +24,8 @@ pure function scalar_1D_initializer_i(x) result(f)
     !! Encapsulate a mimetic matrix with a corresponding matrix-vector product operator
     private
     double precision, allocatable :: upper_(:,:), inner_(:), lower_(:,:)
+  contains
+    procedure to_file_t
   end type
 
   type gradient_operator_1D_t
@@ -46,6 +49,16 @@ pure function scalar_1D_initializer_i(x) result(f)
     procedure, non_overridable, private :: grad
   end type
 
+  interface
+
+     pure module function to_file_t(self) result(file)
+       implicit none
+       class(mimetic_matrix_1D_t), intent(in) :: self
+       type(file_t) file
+     end function
+
+  end interface
+
   interface scalar_1D_t
 
     pure module function construct_from_function(initializer, order, cells, x_min, x_max) result(scalar_1D)
diff --git a/test/gradient_operator_1D_test_m.F90 b/test/gradient_operator_1D_test_m.F90
index 64cba1e8..ce11e948 100644
--- a/test/gradient_operator_1D_test_m.F90
+++ b/test/gradient_operator_1D_test_m.F90
@@ -10,7 +10,9 @@ module gradient_operator_1D_test_m
     ,test_result_t &
     ,operator(//) &
     ,operator(.all.) &
+    ,operator(.also.) &
     ,operator(.approximates.) &
+    ,operator(.csv.) &
     ,operator(.within.)
   use mole_m, only : scalar_1D_t, gradient_1D_t, scalar_1D_initializer_i
 #if ! HAVE_PROCEDURE_ACTUAL_FOR_POINTER_DUMMY
@@ -39,9 +41,9 @@ function results() result(test_results)
     type(gradient_operator_1D_test_t) gradient_operator_1D_test
     type(test_result_t), allocatable :: test_results(:)
     test_results = gradient_operator_1D_test%run([ & 
-       test_description_t('computing the 2nd-order 1D gradient of a constant within a tolerance of ' // string_t(tight_tolerance), check_grad_const) &
-      ,test_description_t('computing the 2nd-order 1D gradient of a line within a tolerance of ' // string_t(loose_tolerance), check_grad_line) &
-      ,test_description_t('computing the 2nd-order 1D gradient of a parabola within a tolerance of ' // string_t(loose_tolerance), check_grad_parabola) &
+       test_description_t('computing 2nd & 4th-order 1D gradients of a constant within tolerance ' // string_t(tight_tolerance), check_grad_const) &
+      ,test_description_t('computing 2nd & 4th-order 1D gradients of a line within tolerance ' // string_t(loose_tolerance), check_grad_line) &
+      ,test_description_t('computing 2nd & 4th-order 1D gradients of a parabola within tolerance ' // string_t(loose_tolerance), check_grad_parabola) &
     ])
   end function
 
@@ -56,9 +58,9 @@ function results() result(test_results)
       ,check_grad_parabola_ptr => check_grad_parabola
 
     test_results = gradient_operator_1D_test%run([ &
-       test_description_t('computing the 2nd-order 1D gradient of a constant within a tolerance of ' // string_t(tight_tolerance), check_grad_const_ptr) &
-      ,test_description_t('computing the 2nd-order 1D gradient of a line within a tolerance of ' // string_t(loose_tolerance), check_grad_line_ptr) &
-      ,test_description_t('computing the 2nd-order 1D gradient of a parabola within a tolerance of ' // string_t(loose_tolerance), check_grad_parabola_ptr) &
+       test_description_t('computing 2nd & 4th-order 1D gradients of a constant within tolerance ' // string_t(tight_tolerance), check_grad_const_ptr) &
+      ,test_description_t('computing 2nd & 4th-order 1D gradients of a line within tolerance ' // string_t(loose_tolerance), check_grad_line_ptr) &
+      ,test_description_t('computing 2nd & 4th-order 1D gradients of a parabola within tolerance ' // string_t(loose_tolerance), check_grad_parabola_ptr) &
     ])
   end function
 
@@ -77,8 +79,11 @@ function check_grad_const() result(test_diagnosis)
     double precision, parameter :: df_dx = 0.
     procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => const
 
-    grad_f = .grad. scalar_1D_t(scalar_1D_initializer, order=2, cells=4, x_min=0D0, x_max=1D0)
-    test_diagnosis = .all. (grad_f%values() .approximates. df_dx .within. loose_tolerance) // " (d(line)/dx)"
+    grad_f = .grad. scalar_1D_t(scalar_1D_initializer, order=2, cells=4, x_min=0D0, x_max=4D0)
+    test_diagnosis = .all. (grad_f%values() .approximates. df_dx .within. loose_tolerance) // " (2nd-order d(line)/dx)"
+
+    grad_f = .grad. scalar_1D_t(scalar_1D_initializer, order=4, cells=8, x_min=0D0, x_max=8D0)
+    test_diagnosis = test_diagnosis .also. (.all. (grad_f%values() .approximates. df_dx .within. loose_tolerance)) // " (4th-order d(line)/dx)"
   end function
 
   pure function line(x) result(y)
@@ -93,8 +98,12 @@ function check_grad_line() result(test_diagnosis)
     double precision, parameter :: df_dx = 14D0
     procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => line
 
-    grad_f = .grad. scalar_1D_t(scalar_1D_initializer, order=2, cells=4, x_min=0D0, x_max=1D0)
-    test_diagnosis = .all. (grad_f%values() .approximates. df_dx .within. loose_tolerance) // " (d(line)/dx)"
+    grad_f = .grad. scalar_1D_t(scalar_1D_initializer, order=2, cells=4, x_min=0D0, x_max=4D0)
+    test_diagnosis = .all. (grad_f%values() .approximates. df_dx .within. loose_tolerance) // " (2nd-order d(line)/dx)"
+
+    grad_f = .grad. scalar_1D_t(scalar_1D_initializer, order=4, cells=8, x_min=0D0, x_max=8D0)
+    test_diagnosis = test_diagnosis .also. (.all. (grad_f%values() .approximates. df_dx .within. loose_tolerance)) // " (4th-order d(line)/dx)"
+
   end function
 
   pure function parabola(x) result(y)
@@ -109,15 +118,23 @@ function check_grad_parabola() result(test_diagnosis)
     type(gradient_1D_t) grad_f
     procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => parabola
 
-    quadratic = scalar_1D_t(scalar_1D_initializer , order=2, cells=4, x_min=0D0, x_max=1D0)
+    quadratic = scalar_1D_t(scalar_1D_initializer , order=2, cells=4, x_min=0D0, x_max=4D0)
     grad_f = .grad. quadratic
 
     associate(x => grad_f%faces())
       associate(df_dx => 14*x + 3)
-        test_diagnosis = .all. (grad_f%values() .approximates. df_dx .within. loose_tolerance) // " (d(parabola)/dx)"
+        test_diagnosis = .all. (grad_f%values() .approximates. df_dx .within. loose_tolerance) // " (2nd-order d(parabola)/dx)"
       end associate
     end associate
 
+    quadratic = scalar_1D_t(scalar_1D_initializer , order=4, cells=8, x_min=0D0, x_max=8D0)
+    grad_f = .grad. quadratic
+
+    associate(x => grad_f%faces())
+      associate(df_dx => 14*x + 3)
+        test_diagnosis = test_diagnosis .also. (.all. (grad_f%values() .approximates. df_dx .within. loose_tolerance)) // " (4th-order d(parabola)/dx)"
+      end associate
+    end associate
   end function
 
 end module
\ No newline at end of file

From ac62a7fd99cfe3edfcc722b54e0562d6f409a261 Mon Sep 17 00:00:00 2001
From: Damian Rouson <rouson@lbl.gov>
Date: Thu, 9 Oct 2025 22:18:07 -0700
Subject: [PATCH 25/40] refac(vector): rename matrix-vector RHS scalar_1D

---
 src/fortran/mimetic_matrix_1D_s.F90 | 10 +++++-----
 src/fortran/scalar_1D_m.f90         |  6 +++---
 2 files changed, 8 insertions(+), 8 deletions(-)

diff --git a/src/fortran/mimetic_matrix_1D_s.F90 b/src/fortran/mimetic_matrix_1D_s.F90
index aa900d26..607c00fe 100644
--- a/src/fortran/mimetic_matrix_1D_s.F90
+++ b/src/fortran/mimetic_matrix_1D_s.F90
@@ -20,18 +20,18 @@
     double precision, allocatable :: product_inner(:)
 
     associate(upper => size(self%upper_,1), lower => size(self%lower_,1))
-      associate(inner_rows => size(vector%scalar_1D_) - (upper + lower + 1))
+      associate(inner_rows => size(scalar_1D%scalar_1D_) - (upper + lower + 1))
 
         allocate(product_inner(inner_rows))
 
-        do concurrent(integer :: row = 1 : inner_rows) default(none) shared(product_inner, self, vector)
-          product_inner(row) = dot_product(self%inner_, vector%scalar_1D_(row + 1 : row + size(self%inner_)))
+        do concurrent(integer :: row = 1 : inner_rows) default(none) shared(product_inner, self, scalar_1D)
+          product_inner(row) = dot_product(self%inner_, scalar_1D%scalar_1D_(row + 1 : row + size(self%inner_)))
         end do
 
         matvec_product = [ &
-           matmul(self%upper_, vector%scalar_1D_(1 : size(self%upper_,2))) &
+           matmul(self%upper_, scalar_1D%scalar_1D_(1 : size(self%upper_,2))) &
           ,product_inner &
-          ,matmul(self%lower_, vector%scalar_1D_(size(vector%scalar_1D_) - size(self%lower_,2) + 1 : )) &
+          ,matmul(self%lower_, scalar_1D%scalar_1D_(size(scalar_1D%scalar_1D_) - size(self%lower_,2) + 1 : )) &
         ]
       end associate
     end associate
diff --git a/src/fortran/scalar_1D_m.f90 b/src/fortran/scalar_1D_m.f90
index fca7b831..dc1e5c25 100644
--- a/src/fortran/scalar_1D_m.f90
+++ b/src/fortran/scalar_1D_m.f90
@@ -119,11 +119,11 @@ pure module function construct_from_components(upper, inner, lower) result(mimet
 
   interface
 
-    pure module function matvec(self, vector) result(matvec_product)
-      !! Apply a matrix operator to a vector
+    pure module function matvec(self, scalar_1D) result(matvec_product)
+      !! Apply a mimetic matrix operator to a vector encapsulated in a scalar_1D_t object
       implicit none
       class(mimetic_matrix_1D_t), intent(in) :: self
-      type(scalar_1D_t), intent(in) :: vector
+      type(scalar_1D_t), intent(in) :: scalar_1D
       double precision, allocatable :: matvec_product(:)
     end function
 

From 79f42ba2179aa8b5bd745eeceac23b0e7fb90c25 Mon Sep 17 00:00:00 2001
From: Damian Rouson <rouson@lbl.gov>
Date: Thu, 9 Oct 2025 22:31:27 -0700
Subject: [PATCH 26/40] doc(gradient_operator_1D): clarify statement label

---
 src/fortran/gradient_operator_1D_s.F90 | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/fortran/gradient_operator_1D_s.F90 b/src/fortran/gradient_operator_1D_s.F90
index ef9eb275..873978a7 100644
--- a/src/fortran/gradient_operator_1D_s.F90
+++ b/src/fortran/gradient_operator_1D_s.F90
@@ -75,10 +75,10 @@ pure function corbino_castillo_Ap(k, dx) result(matrix_block)
 
     associate(A => corbino_castillo_A(k, dx))
       allocate(matrix_block , mold=A)
-      reverse_elements_within_matrix_block_and_flip_sign: &
+      reverse_elements_within_rows_and_flip_sign: &
       do concurrent(integer :: row = 1:size(matrix_block,1)) default(none) shared(matrix_block, A)
         matrix_block(row,:) = -A(row,size(A,2):1:-1)
-      end do reverse_elements_within_matrix_block_and_flip_sign
+      end do reverse_elements_within_rows_and_flip_sign
       reverse_elements_within_columns: &
       do concurrent(integer :: column = 1 : size(matrix_block,2)) default(none) shared(matrix_block)
         matrix_block(:,column) = matrix_block(size(matrix_block,1):1:-1,column)

From 88e8c5e6ce63e613818605691ce02b8d411b5eb0 Mon Sep 17 00:00:00 2001
From: Damian Rouson <rouson@lbl.gov>
Date: Fri, 10 Oct 2025 14:18:51 -0700
Subject: [PATCH 27/40] test(order): 2nd-/4th-order accuracy tests pass

This commit adds tests of the ratio of the log of the maximum
absolute error for gradients of a sinusoidal function.  The tests
verify that 2nd- and 4th-order mimetic discretizations converge
to the expected gradient values at a rate proportional to dx raised
to the 2nd and 4th powers of the error within a tolerance of 5%.
---
 src/fortran/gradient_operator_1D_s.F90 | 20 +++++++-
 src/fortran/scalar_1D_s.F90            | 12 ++++-
 test/gradient_operator_1D_test_m.F90   | 70 ++++++++++++++++++++++++--
 3 files changed, 95 insertions(+), 7 deletions(-)

diff --git a/src/fortran/gradient_operator_1D_s.F90 b/src/fortran/gradient_operator_1D_s.F90
index 873978a7..4d06263f 100644
--- a/src/fortran/gradient_operator_1D_s.F90
+++ b/src/fortran/gradient_operator_1D_s.F90
@@ -67,7 +67,6 @@ pure function corbino_castillo_M(k, dx) result(row)
 
 #if HAVE_DO_CONCURRENT_TYPE_SPEC_SUPPORT
 
-
   pure function corbino_castillo_Ap(k, dx) result(matrix_block)
     integer, intent(in) :: k
     double precision, intent(in) :: dx
@@ -88,6 +87,25 @@ pure function corbino_castillo_Ap(k, dx) result(matrix_block)
  
 #else
 
+  pure function corbino_castillo_Ap(k, dx) result(matrix_block)
+    integer, intent(in) :: k
+    double precision, intent(in) :: dx
+    double precision, allocatable :: matrix_block(:,:)
+    integer row, column
+
+    associate(A => corbino_castillo_A(k, dx))
+      allocate(matrix_block , mold=A)
+      reverse_elements_within_rows_and_flip_sign: &
+      do concurrent(row = 1:size(matrix_block,1)) default(none) shared(matrix_block, A)
+        matrix_block(row,:) = -A(row,size(A,2):1:-1)
+      end do reverse_elements_within_rows_and_flip_sign
+      reverse_elements_within_columns: &
+      do concurrent(column = 1 : size(matrix_block,2)) default(none) shared(matrix_block)
+        matrix_block(:,column) = matrix_block(size(matrix_block,1):1:-1,column)
+      end do reverse_elements_within_columns
+    end associate
+  end function
+ 
 #endif
 
 end submodule gradient_operator_1D_s
\ No newline at end of file
diff --git a/src/fortran/scalar_1D_s.F90 b/src/fortran/scalar_1D_s.F90
index 57425e9a..c7dde20f 100644
--- a/src/fortran/scalar_1D_s.F90
+++ b/src/fortran/scalar_1D_s.F90
@@ -6,7 +6,15 @@
 
 contains
 
-  module procedure construct_from_function
+  pure module function construct_from_function(initializer, order, cells, x_min, x_max) result(scalar_1D)
+    implicit none
+    procedure(scalar_1D_initializer_i), pointer :: initializer
+    integer, intent(in) :: order !! order of accuracy
+    integer, intent(in) :: cells !! number of grid cells spanning the domain
+    double precision, intent(in) :: x_min !! grid location minimum
+    double precision, intent(in) :: x_max !! grid location maximum
+    type(scalar_1D_t) scalar_1D
+
     call_julienne_assert(x_max .greaterThan. x_min)
     call_julienne_assert(cells .isAtLeast. 2*order)
 
@@ -15,7 +23,7 @@
     scalar_1D%cells_ = cells
     scalar_1D%gradient_operator_1D_ = gradient_operator_1D_t(k=order, dx=(x_max - x_min)/cells, m=cells)
     scalar_1D%scalar_1D_ = initializer(grid_(x_min, x_max, cells))
-  end procedure
+  end function
 
   pure function grid_(x_min, x_max, cells) result(x)
     double precision, intent(in) :: x_min, x_max
diff --git a/test/gradient_operator_1D_test_m.F90 b/test/gradient_operator_1D_test_m.F90
index ce11e948..1ef85f7b 100644
--- a/test/gradient_operator_1D_test_m.F90
+++ b/test/gradient_operator_1D_test_m.F90
@@ -26,7 +26,7 @@ module gradient_operator_1D_test_m
     procedure, nopass :: results
   end type
 
-  double precision, parameter :: tight_tolerance = 1D-14, loose_tolerance = 1D-12
+  double precision, parameter :: tight_tolerance = 1D-14, loose_tolerance = 1D-12, rough_tolerance = 1D-02, crude_tolerance = 5D-02
 
 contains
 
@@ -41,9 +41,11 @@ function results() result(test_results)
     type(gradient_operator_1D_test_t) gradient_operator_1D_test
     type(test_result_t), allocatable :: test_results(:)
     test_results = gradient_operator_1D_test%run([ & 
-       test_description_t('computing 2nd & 4th-order 1D gradients of a constant within tolerance ' // string_t(tight_tolerance), check_grad_const) &
-      ,test_description_t('computing 2nd & 4th-order 1D gradients of a line within tolerance ' // string_t(loose_tolerance), check_grad_line) &
-      ,test_description_t('computing 2nd & 4th-order 1D gradients of a parabola within tolerance ' // string_t(loose_tolerance), check_grad_parabola) &
+       test_description_t('computing 2nd- & 4th-order 1D gradients of a constant within tolerance ' // string_t(tight_tolerance), check_grad_const) &
+      ,test_description_t('computing 2nd- & 4th-order 1D gradients of a line within tolerance ' // string_t(loose_tolerance), check_grad_line) &
+      ,test_description_t('computing 2nd- & 4th-order 1D gradients of a parabola within tolerance ' // string_t(loose_tolerance), check_grad_parabola) &
+      ,test_description_t('computing 2nd-order 1D gradients of a sinusoid with a convergence rate of 2 within tolerance ' // string_t(crude_tolerance), check_2nd_order_grad_sinusoid) &
+      ,test_description_t('computing 4th-order 1D gradients of a sinusoid with a convergence rate of 4 within tolerance ' // string_t(crude_tolerance), check_4th_order_grad_sinusoid) &
     ])
   end function
 
@@ -137,4 +139,64 @@ function check_grad_parabola() result(test_diagnosis)
     end associate
   end function
 
+  pure function sinusoid(x) result(y)
+    double precision, intent(in) :: x(:)
+    double precision, allocatable :: y(:)
+    y = sin(x) + cos(x)
+  end function
+
+  function check_2nd_order_grad_sinusoid() result(test_diagnosis)
+    type(test_diagnosis_t) test_diagnosis
+    type(scalar_1D_t) coarse, fine
+    type(gradient_1D_t) grad_coarse, grad_fine
+    procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => sinusoid
+    double precision, parameter :: pi = 3.141592653589793D0
+    integer, parameter :: order_desired = 2, coarse_cells=100, fine_cells=1000
+
+    coarse = scalar_1D_t(scalar_1D_initializer , order=order_desired, cells=coarse_cells, x_min=0D0, x_max=2*pi)
+    fine   = scalar_1D_t(scalar_1D_initializer , order=order_desired, cells=fine_cells  , x_min=0D0, x_max=2*pi)
+
+    grad_coarse = .grad. coarse
+    grad_fine   = .grad. fine
+
+    associate(x_coarse => grad_coarse%faces(), x_fine => grad_fine%faces())
+      associate(df_dx_coarse => cos(x_coarse) - sin(x_coarse), df_dx_fine => cos(x_fine) - sin(x_fine), grad_coarse_values => grad_coarse%values(), grad_fine_values => grad_fine%values())
+        test_diagnosis = .all. (grad_coarse_values .approximates. df_dx_coarse .within. rough_tolerance) // " (2nd-order d(sinusoid)/dx point-wise errors)"
+        test_diagnosis = test_diagnosis .also. (.all. (grad_fine_values .approximates. df_dx_fine .within. rough_tolerance)) // " (2nd-order d(sinusoid)/dx point-wise)"
+        associate(error_coarse_max => maxval(abs(grad_coarse_values - df_dx_coarse)), error_fine_max => maxval(abs(grad_fine_values - df_dx_fine)))
+          associate(order_actual => log(error_coarse_max/error_fine_max)/log(dble(fine_cells)/coarse_cells))
+            test_diagnosis = test_diagnosis .also. (order_actual .approximates. dble(order_desired) .within. crude_tolerance)  // " (2nd-order d(sinusoid)/dx order of accuracy)"
+          end associate
+        end associate
+      end associate
+    end associate
+  end function
+
+  function check_4th_order_grad_sinusoid() result(test_diagnosis)
+    type(test_diagnosis_t) test_diagnosis
+    type(scalar_1D_t) coarse, fine
+    type(gradient_1D_t) grad_coarse, grad_fine
+    procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => sinusoid
+    double precision, parameter :: pi = 3.141592653589793D0
+    integer, parameter :: order_desired = 4, coarse_cells=100, fine_cells=1000
+
+    coarse = scalar_1D_t(scalar_1D_initializer , order=order_desired, cells=coarse_cells, x_min=0D0, x_max=2*pi)
+    fine   = scalar_1D_t(scalar_1D_initializer , order=order_desired, cells=fine_cells  , x_min=0D0, x_max=2*pi)
+
+    grad_coarse = .grad. coarse
+    grad_fine   = .grad. fine
+
+    associate(x_coarse => grad_coarse%faces(), x_fine => grad_fine%faces())
+      associate(df_dx_coarse => cos(x_coarse) - sin(x_coarse), df_dx_fine => cos(x_fine) - sin(x_fine), grad_coarse_values => grad_coarse%values(), grad_fine_values => grad_fine%values())
+        test_diagnosis = .all. (grad_coarse_values .approximates. df_dx_coarse .within. rough_tolerance) // " (4th-order d(sinusoid)/dx point-wise errors)"
+        test_diagnosis = test_diagnosis .also. (.all. (grad_fine_values .approximates. df_dx_fine .within. rough_tolerance)) // " (4th-order d(sinusoid)/dx point-wise)"
+        associate(error_coarse_max => maxval(abs(grad_coarse_values - df_dx_coarse)), error_fine_max => maxval(abs(grad_fine_values - df_dx_fine)))
+          associate(order_actual => log(error_coarse_max/error_fine_max)/log(dble(fine_cells)/coarse_cells))
+            test_diagnosis = test_diagnosis .also. (order_actual .approximates. dble(order_desired) .within. crude_tolerance)  // " (4th-order d(sinusoid)/dx order of accuracy)"
+          end associate
+        end associate
+      end associate
+    end associate
+  end function
+
 end module
\ No newline at end of file

From bf84505a00a7f80413db8ac95fc1e1c5637be4bb Mon Sep 17 00:00:00 2001
From: Damian Rouson <rouson@lbl.gov>
Date: Sat, 11 Oct 2025 12:14:51 -0700
Subject: [PATCH 28/40] build(gfortran): work around concurrent type-spec

---
 src/fortran/mimetic_matrix_1D_s.F90 | 23 +++++++++++++++++++++++
 1 file changed, 23 insertions(+)

diff --git a/src/fortran/mimetic_matrix_1D_s.F90 b/src/fortran/mimetic_matrix_1D_s.F90
index 607c00fe..c7b0739e 100644
--- a/src/fortran/mimetic_matrix_1D_s.F90
+++ b/src/fortran/mimetic_matrix_1D_s.F90
@@ -39,6 +39,29 @@
 
 #else
 
+  module procedure matvec
+
+    integer row
+    double precision, allocatable :: product_inner(:)
+
+    associate(upper => size(self%upper_,1), lower => size(self%lower_,1))
+      associate(inner_rows => size(scalar_1D%scalar_1D_) - (upper + lower + 1))
+
+        allocate(product_inner(inner_rows))
+
+        do concurrent(row = 1 : inner_rows) default(none) shared(product_inner, self, scalar_1D)
+          product_inner(row) = dot_product(self%inner_, scalar_1D%scalar_1D_(row + 1 : row + size(self%inner_)))
+        end do
+
+        matvec_product = [ &
+           matmul(self%upper_, scalar_1D%scalar_1D_(1 : size(self%upper_,2))) &
+          ,product_inner &
+          ,matmul(self%lower_, scalar_1D%scalar_1D_(size(scalar_1D%scalar_1D_) - size(self%lower_,2) + 1 : )) &
+        ]
+      end associate
+    end associate
+  end procedure
+
 #endif
 
   module procedure to_file_t

From 1e0baaf48902a7ca7a6d509f6be33dac429e7b27 Mon Sep 17 00:00:00 2001
From: Damian Rouson <rouson@lbl.gov>
Date: Sun, 12 Oct 2025 11:44:48 -0700
Subject: [PATCH 29/40] build(gfortran-14): add locality specifier macro

To work around gfortran-14 not supporting `do concurrent` locality
specifiers, this commit removes specifiers conditionally via a new
macro: HAVE_LOCALITY_SPECIFIER_SUPPORT defined in
include/mole-language-support.F90.
---
 include/mole-language-support.F90      | 13 ++++++++-----
 src/fortran/gradient_operator_1D_s.F90 |  6 +++---
 src/fortran/mimetic_matrix_1D_s.F90    |  4 ++--
 3 files changed, 13 insertions(+), 10 deletions(-)

diff --git a/include/mole-language-support.F90 b/include/mole-language-support.F90
index 50171f15..f8896b27 100644
--- a/include/mole-language-support.F90
+++ b/include/mole-language-support.F90
@@ -4,17 +4,12 @@
 #ifndef MOLE_LANGUAGE_SUPPORT
 #define MOLE_LANGUAGE_SUPPORT
 
-
 #ifdef __GNUC__
 #  define GCC_VERSION (__GNUC__ * 10000 + __GNUC_MINOR__ * 100 + __GNUC_PATCHLEVEL__)
 #else
 #  define GCC_VERSION 0
 #endif
 
-#if __GNUC__ && ( __GNUC__ < 14 || (__GNUC__ == 14 && __GNUC_MINOR__ < 3) )
-#define GCC_GE_MINIMUM
-#endif
-
 #ifndef HAVE_DO_CONCURRENT_TYPE_SPEC_SUPPORT
 #  if defined(_CRAYFTN) || defined(__INTEL_COMPILER) || defined(NAGFOR) || defined(__flang__)
 #    define HAVE_DO_CONCURRENT_TYPE_SPEC_SUPPORT 1
@@ -23,4 +18,12 @@
 #  endif
 #endif
 
+#ifndef HAVE_LOCALITY_SPECIFIER_SUPPORT
+#  if defined(NAGFOR) || defined(__flang__) || defined(__INTEL_COMPILER) || defined(_CRAYFTN)
+#    define HAVE_LOCALITY_SPECIFIER_SUPPORT 1
+#  else
+#    define HAVE_LOCALITY_SPECIFIER_SUPPORT 0
+#  endif
+#endif
+
 #endif
diff --git a/src/fortran/gradient_operator_1D_s.F90 b/src/fortran/gradient_operator_1D_s.F90
index 4d06263f..f5f41ab2 100644
--- a/src/fortran/gradient_operator_1D_s.F90
+++ b/src/fortran/gradient_operator_1D_s.F90
@@ -65,7 +65,7 @@ pure function corbino_castillo_M(k, dx) result(row)
 
   end function
 
-#if HAVE_DO_CONCURRENT_TYPE_SPEC_SUPPORT
+#if HAVE_DO_CONCURRENT_TYPE_SPEC_SUPPORT && HAVE_LOCALITY_SPECIFIER_SUPPORT
 
   pure function corbino_castillo_Ap(k, dx) result(matrix_block)
     integer, intent(in) :: k
@@ -96,11 +96,11 @@ pure function corbino_castillo_Ap(k, dx) result(matrix_block)
     associate(A => corbino_castillo_A(k, dx))
       allocate(matrix_block , mold=A)
       reverse_elements_within_rows_and_flip_sign: &
-      do concurrent(row = 1:size(matrix_block,1)) default(none) shared(matrix_block, A)
+      do concurrent(row = 1:size(matrix_block,1))
         matrix_block(row,:) = -A(row,size(A,2):1:-1)
       end do reverse_elements_within_rows_and_flip_sign
       reverse_elements_within_columns: &
-      do concurrent(column = 1 : size(matrix_block,2)) default(none) shared(matrix_block)
+      do concurrent(column = 1 : size(matrix_block,2))
         matrix_block(:,column) = matrix_block(size(matrix_block,1):1:-1,column)
       end do reverse_elements_within_columns
     end associate
diff --git a/src/fortran/mimetic_matrix_1D_s.F90 b/src/fortran/mimetic_matrix_1D_s.F90
index c7b0739e..1138d8bc 100644
--- a/src/fortran/mimetic_matrix_1D_s.F90
+++ b/src/fortran/mimetic_matrix_1D_s.F90
@@ -13,7 +13,7 @@
     mimetic_matrix_1D%lower_ = lower
   end procedure
 
-#if HAVE_DO_CONCURRENT_TYPE_SPEC_SUPPORT
+#if HAVE_DO_CONCURRENT_TYPE_SPEC_SUPPORT && HAVE_LOCALITY_SPECIFIER_SUPPORT
 
   module procedure matvec
 
@@ -49,7 +49,7 @@
 
         allocate(product_inner(inner_rows))
 
-        do concurrent(row = 1 : inner_rows) default(none) shared(product_inner, self, scalar_1D)
+        do concurrent(row = 1 : inner_rows)
           product_inner(row) = dot_product(self%inner_, scalar_1D%scalar_1D_(row + 1 : row + size(self%inner_)))
         end do
 

From 9abbde7f0765a2d5e5dc419c333c88644181d310 Mon Sep 17 00:00:00 2001
From: Damian Rouson <rouson@lbl.gov>
Date: Sun, 12 Oct 2025 12:09:14 -0700
Subject: [PATCH 30/40] build(gfortran): use GCC ver to define local spec

---
 include/mole-language-support.F90 | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/include/mole-language-support.F90 b/include/mole-language-support.F90
index f8896b27..9b332563 100644
--- a/include/mole-language-support.F90
+++ b/include/mole-language-support.F90
@@ -19,7 +19,7 @@
 #endif
 
 #ifndef HAVE_LOCALITY_SPECIFIER_SUPPORT
-#  if defined(NAGFOR) || defined(__flang__) || defined(__INTEL_COMPILER) || defined(_CRAYFTN)
+#  if defined(NAGFOR) || defined(__flang__) || defined(__INTEL_COMPILER) || defined(_CRAYFTN) || (GCC_VERSION >= 150100)
 #    define HAVE_LOCALITY_SPECIFIER_SUPPORT 1
 #  else
 #    define HAVE_LOCALITY_SPECIFIER_SUPPORT 0

From 3cc1e8135380e5f22dd3dd81a20bfeaba2fe407c Mon Sep 17 00:00:00 2001
From: Damian Rouson <rouson@lbl.gov>
Date: Sun, 12 Oct 2025 23:18:35 -0700
Subject: [PATCH 31/40] test(CI): build/test gradient-operator branch

---
 .github/workflows/build.yml | 182 ++++++++++++++++++++++++++++++++++++
 1 file changed, 182 insertions(+)
 create mode 100644 .github/workflows/build.yml

diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml
new file mode 100644
index 00000000..7b83ba88
--- /dev/null
+++ b/.github/workflows/build.yml
@@ -0,0 +1,182 @@
+name: Build
+
+on: [push]
+
+
+jobs:
+  build:
+    name: ${{ matrix.compiler }}-${{ matrix.version }} (${{ matrix.os }})
+    runs-on: ${{ matrix.os }}
+    strategy:
+      fail-fast: false
+      matrix:
+        os: [macos-13, macos-14, macos-15, ubuntu-22.04, ubuntu-24.04]
+        compiler: [ gfortran ]
+        version: [ 13, 14 ]
+        extra_flags: [ -g -O3 ]
+
+        exclude:
+          - os: ubuntu-22.04
+            version: 13 # no package available
+          - os: ubuntu-22.04
+            version: 14 # no package available
+
+        include:
+          # --- LLVM flang coverage ---
+
+          - os: macos-14
+            compiler: flang
+            version: 20
+          - os: macos-15
+            compiler: flang
+            version: 20
+
+          # https://hub.docker.com/r/snowstep/llvm/tags
+          - os: ubuntu-24.04
+            compiler: flang
+            version: latest
+            container: snowstep/llvm:noble
+          - os: ubuntu-22.04
+            compiler: flang
+            version: latest
+            container: snowstep/llvm:jammy
+
+          # https://hub.docker.com/r/phhargrove/llvm-flang/tags
+          - os: ubuntu-24.04
+            compiler: flang
+            version: 20
+            container: phhargrove/llvm-flang:20.1.0-1
+          - os: ubuntu-24.04
+            compiler: flang
+            version: 19
+            extra_flags: -g -mmlir -allow-assumed-rank -O3
+            container: phhargrove/llvm-flang:19.1.1-1
+
+          # --- Intel coverage ---
+
+# https://hub.docker.com/r/intel/fortran-essentials/tags
+          - os: ubuntu-24.04
+            compiler: ifx
+            version: 2025.2.0
+            error_stop_code: 128
+            container: intel/fortran-essentials:2025.2.0-0-devel-ubuntu24.04
+
+          - os: ubuntu-24.04
+            compiler: ifx
+            version: 2025.2.2
+            error_stop_code: 128
+            container: intel/fortran-essentials:2025.2.2-0-devel-ubuntu24.04
+
+          - os: ubuntu-24.04
+            compiler: ifx
+            version: latest
+            error_stop_code: 128
+            container: intel/fortran-essentials:latest
+
+          # --- LFortran coverage ---
+
+          # https://hub.docker.com/r/phhargrove/lfortran/tags
+          #- os: ubuntu-24.04
+          #  compiler: lfortran
+          #  version: 0.54.0
+          #  container: phhargrove/lfortran:0.54.0-1
+
+    container:
+      image: ${{ matrix.container }}
+
+    env:
+      COMPILER_VERSION: ${{ matrix.version }}
+      FC: ${{ matrix.compiler }}
+      FFLAGS: ${{ matrix.extra_flags }}
+      FPM_FLAGS:  --profile release --verbose
+
+    steps:
+    - name: Checkout code
+      uses: actions/checkout@v5
+      with:
+        ref: gradient-operator
+
+    - name: Install Dependencies Ubuntu
+      if: ${{ contains(matrix.os, 'ubuntu') && matrix.compiler == 'gfortran' && 0 }}
+      run: |
+        sudo apt-get update
+        sudo apt list -a 'gfortran-*'
+        sudo apt install -y gfortran-${COMPILER_VERSION} build-essential
+        sudo update-alternatives --install /usr/bin/gcc gcc /usr/bin/gcc-${COMPILER_VERSION} 100 \
+        --slave /usr/bin/gfortran gfortran /usr/bin/gfortran-${COMPILER_VERSION}
+
+
+    - name: Install Ubuntu Container Dependencies
+      if: ${{ contains(matrix.os, 'ubuntu') && matrix.container != '' && !contains(matrix.container, 'phhargrove') }}
+      run: |
+        set -x
+        apt update
+        apt install -y build-essential pkg-config make git curl
+
+    - name: Install macOS Dependencies
+      if: contains(matrix.os, 'macos') && matrix.compiler == 'flang'
+      run: |
+        set -x
+        brew update
+        brew install llvm@${COMPILER_VERSION} flang
+        # workaround issue #228: clang cannot find homebrew flang's C header
+        for p in /opt/homebrew /usr/local $(brew --prefix) ; do find $p/Cellar/flang -name ISO_Fortran_binding.h 2>/dev/null || true ; done
+        echo "CFLAGS=-I$(dirname $(find $(brew --prefix)/Cellar/flang -name ISO_Fortran_binding.h | head -1)) ${CFLAGS}" >> "$GITHUB_ENV"
+        # Prepend homebrew clang to PATH:
+        echo "PATH=$(brew --prefix)/opt/llvm/bin:${PATH}" >> "$GITHUB_ENV"
+
+    - name: Setup Compilers
+      run: |
+        set -x
+        if test "$FC" = "flang" ; then \
+          echo "FPM_FC=flang-new" >> "$GITHUB_ENV" ; \
+        elif test "$FC" = "ifx" ; then \
+          echo "FPM_FC=ifx" >> "$GITHUB_ENV" ; \
+          echo "FFLAGS=-fpp $FFLAGS" >> "$GITHUB_ENV" ; \
+          echo "FFLAGS=-coarray $FFLAGS" >> "$GITHUB_ENV" ; \
+        elif test "$FC" = "lfortran" ; then \
+          echo "FPM_FC=lfortran" >> "$GITHUB_ENV" ; \
+          echo "FFLAGS=--cpp $FFLAGS" >> "$GITHUB_ENV" ; \
+        else \
+          echo "FPM_FC=gfortran-${COMPILER_VERSION}" >> "$GITHUB_ENV" ; \
+          echo "FFLAGS=-ffree-line-length-0 $FFLAGS" >> "$GITHUB_ENV" ; \
+        fi
+        if test -n "${{ matrix.error_stop_code }}" ; then \
+          echo "ERROR_STOP_CODE=${{ matrix.error_stop_code }}" >> "$GITHUB_ENV" ; \
+        else \
+          echo "ERROR_STOP_CODE=1" >> "$GITHUB_ENV" ; \
+        fi
+
+    - name: Setup FPM
+      uses: fortran-lang/setup-fpm@main
+      if: ${{ !contains(matrix.os, 'macos') || matrix.os == 'macos-13' }}
+      with:
+        github-token: ${{ secrets.GITHUB_TOKEN }}
+        fpm-version: latest
+
+    - name: Build FPM
+      # no macos-arm64 fpm distro, build from source
+      if: ${{ contains(matrix.os, 'macos') && matrix.os != 'macos-13' }}
+      run: |
+        set -x
+        curl --retry 5 -LOsS https://github.com/fortran-lang/fpm/releases/download/v0.11.0/fpm-0.11.0.F90
+        mkdir fpm-temp
+        gfortran-14 -o fpm-temp/fpm fpm-0.11.0.F90
+        echo "PATH=`pwd`/fpm-temp:${PATH}" >> "$GITHUB_ENV"
+
+    - name: Version info
+      run: |
+        set -x
+        echo == TOOL VERSIONS ==
+        uname -a
+        if test -r /etc/os-release ; then cat /etc/os-release ; fi
+        ${FPM_FC} --version
+        fpm --version
+
+    - name: Build and Test (Assertions OFF)
+      run: |
+        set -x
+        echo "FPM_FLAGS=${FPM_FLAGS}"
+        echo "FFLAGS=$FFLAGS"
+        ls -lt
+        fpm test ${FPM_FLAGS} --flag "$FFLAGS"

From 5a6020b436bceca13d7af5209c1ae10239e3a7f7 Mon Sep 17 00:00:00 2001
From: Damian Rouson <rouson@lbl.gov>
Date: Sun, 12 Oct 2025 23:28:18 -0700
Subject: [PATCH 32/40] test(CI): reduce test matrix

---
 .github/workflows/build.yml | 33 ++-------------------------------
 1 file changed, 2 insertions(+), 31 deletions(-)

diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml
index 7b83ba88..12f0d0f7 100644
--- a/.github/workflows/build.yml
+++ b/.github/workflows/build.yml
@@ -1,6 +1,6 @@
 name: Build
 
-on: [push]
+on: [push, pull_request]
 
 
 jobs:
@@ -10,7 +10,7 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        os: [macos-13, macos-14, macos-15, ubuntu-22.04, ubuntu-24.04]
+        os: [macos-15, ubuntu-22.04, ubuntu-24.04]
         compiler: [ gfortran ]
         version: [ 13, 14 ]
         extra_flags: [ -g -O3 ]
@@ -52,35 +52,6 @@ jobs:
             extra_flags: -g -mmlir -allow-assumed-rank -O3
             container: phhargrove/llvm-flang:19.1.1-1
 
-          # --- Intel coverage ---
-
-# https://hub.docker.com/r/intel/fortran-essentials/tags
-          - os: ubuntu-24.04
-            compiler: ifx
-            version: 2025.2.0
-            error_stop_code: 128
-            container: intel/fortran-essentials:2025.2.0-0-devel-ubuntu24.04
-
-          - os: ubuntu-24.04
-            compiler: ifx
-            version: 2025.2.2
-            error_stop_code: 128
-            container: intel/fortran-essentials:2025.2.2-0-devel-ubuntu24.04
-
-          - os: ubuntu-24.04
-            compiler: ifx
-            version: latest
-            error_stop_code: 128
-            container: intel/fortran-essentials:latest
-
-          # --- LFortran coverage ---
-
-          # https://hub.docker.com/r/phhargrove/lfortran/tags
-          #- os: ubuntu-24.04
-          #  compiler: lfortran
-          #  version: 0.54.0
-          #  container: phhargrove/lfortran:0.54.0-1
-
     container:
       image: ${{ matrix.container }}
 

From 04a1307934b13e6e66bba3a7d49c20985ab831b5 Mon Sep 17 00:00:00 2001
From: Dan Bonachea <dobonachea@lbl.gov>
Date: Thu, 16 Oct 2025 14:55:54 -0700
Subject: [PATCH 33/40] Remove hard-coded checkout branch, expand version
 status step

---
 .github/workflows/build.yml | 17 ++++++++++-------
 1 file changed, 10 insertions(+), 7 deletions(-)

diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml
index 12f0d0f7..ca75d65c 100644
--- a/.github/workflows/build.yml
+++ b/.github/workflows/build.yml
@@ -64,8 +64,6 @@ jobs:
     steps:
     - name: Checkout code
       uses: actions/checkout@v5
-      with:
-        ref: gradient-operator
 
     - name: Install Dependencies Ubuntu
       if: ${{ contains(matrix.os, 'ubuntu') && matrix.compiler == 'gfortran' && 0 }}
@@ -137,12 +135,17 @@ jobs:
 
     - name: Version info
       run: |
-        set -x
-        echo == TOOL VERSIONS ==
+        ( type -p git > /dev/null && (echo == Commit Info == ; git log -n 1 ; echo ; echo ) ; exit 0 )
+        echo == Platform version info == 
         uname -a
-        if test -r /etc/os-release ; then cat /etc/os-release ; fi
-        ${FPM_FC} --version
-        fpm --version
+        if test -r /etc/os-release ; then grep -e NAME -e VERSION /etc/os-release  ; fi
+        if test -x /usr/bin/sw_vers ; then /usr/bin/sw_vers ; fi
+        echo
+        echo PATH="$PATH"
+        for tool in ${FPM_FC} fpm ; do
+          ( echo ; set -x ; w=$(which $tool) ; ls -al $w ; ls -alhL $w ; $tool --version )
+        done
+
 
     - name: Build and Test (Assertions OFF)
       run: |

From 80429ae02e40c95d6ca3bf05d325a8355193b704 Mon Sep 17 00:00:00 2001
From: Dan Bonachea <dobonachea@lbl.gov>
Date: Thu, 16 Oct 2025 15:36:08 -0700
Subject: [PATCH 34/40] Expand set of OS and flang/gfortran

Also, always use setup-fpm
---
 .github/workflows/build.yml | 37 +++++++++++++++++++------------------
 1 file changed, 19 insertions(+), 18 deletions(-)

diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml
index ca75d65c..1572af4b 100644
--- a/.github/workflows/build.yml
+++ b/.github/workflows/build.yml
@@ -10,26 +10,25 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        os: [macos-15, ubuntu-22.04, ubuntu-24.04]
+        os: [ macos-14, macos-15, macos-15-intel, ubuntu-24.04]
         compiler: [ gfortran ]
-        version: [ 13, 14 ]
+        version: [ 13, 14, 15 ]
         extra_flags: [ -g -O3 ]
 
         exclude:
-          - os: ubuntu-22.04
-            version: 13 # no package available
-          - os: ubuntu-22.04
-            version: 14 # no package available
+          - os: ubuntu-24.04
+            version: 15 # no package available (yet?)
 
         include:
+
           # --- LLVM flang coverage ---
 
           - os: macos-14
             compiler: flang
-            version: 20
+            version: 21
           - os: macos-15
             compiler: flang
-            version: 20
+            version: 21
 
           # https://hub.docker.com/r/snowstep/llvm/tags
           - os: ubuntu-24.04
@@ -42,15 +41,20 @@ jobs:
             container: snowstep/llvm:jammy
 
           # https://hub.docker.com/r/phhargrove/llvm-flang/tags
+          - os: ubuntu-24.04
+            compiler: flang
+            version: 21
+            network: smp
+            container: phhargrove/llvm-flang:21.1.0-latest
           - os: ubuntu-24.04
             compiler: flang
             version: 20
-            container: phhargrove/llvm-flang:20.1.0-1
+            container: phhargrove/llvm-flang:20.1.0-latest
           - os: ubuntu-24.04
             compiler: flang
             version: 19
             extra_flags: -g -mmlir -allow-assumed-rank -O3
-            container: phhargrove/llvm-flang:19.1.1-1
+            container: phhargrove/llvm-flang:19.1.1-latest
 
     container:
       image: ${{ matrix.container }}
@@ -60,20 +64,19 @@ jobs:
       FC: ${{ matrix.compiler }}
       FFLAGS: ${{ matrix.extra_flags }}
       FPM_FLAGS:  --profile release --verbose
+      FPM_CXXFLAGS: -std=c++14 # needed for Armadillo on macos-14
 
     steps:
     - name: Checkout code
       uses: actions/checkout@v5
 
-    - name: Install Dependencies Ubuntu
-      if: ${{ contains(matrix.os, 'ubuntu') && matrix.compiler == 'gfortran' && 0 }}
+    - name: Install Ubuntu Native Dependencies
+      if: ${{ 0 && contains(matrix.os, 'ubuntu') && matrix.container == '' && matrix.compiler == 'gfortran' }}
       run: |
+        set -x
         sudo apt-get update
         sudo apt list -a 'gfortran-*'
         sudo apt install -y gfortran-${COMPILER_VERSION} build-essential
-        sudo update-alternatives --install /usr/bin/gcc gcc /usr/bin/gcc-${COMPILER_VERSION} 100 \
-        --slave /usr/bin/gfortran gfortran /usr/bin/gfortran-${COMPILER_VERSION}
-
 
     - name: Install Ubuntu Container Dependencies
       if: ${{ contains(matrix.os, 'ubuntu') && matrix.container != '' && !contains(matrix.container, 'phhargrove') }}
@@ -118,14 +121,12 @@ jobs:
 
     - name: Setup FPM
       uses: fortran-lang/setup-fpm@main
-      if: ${{ !contains(matrix.os, 'macos') || matrix.os == 'macos-13' }}
       with:
         github-token: ${{ secrets.GITHUB_TOKEN }}
         fpm-version: latest
 
     - name: Build FPM
-      # no macos-arm64 fpm distro, build from source
-      if: ${{ contains(matrix.os, 'macos') && matrix.os != 'macos-13' }}
+      if: false
       run: |
         set -x
         curl --retry 5 -LOsS https://github.com/fortran-lang/fpm/releases/download/v0.11.0/fpm-0.11.0.F90

From 4eec9b230861f0910335f1e512b9a68893b68e91 Mon Sep 17 00:00:00 2001
From: Dan Bonachea <dobonachea@lbl.gov>
Date: Thu, 16 Oct 2025 19:33:50 -0700
Subject: [PATCH 35/40] Misc cleanups

---
 .github/workflows/build.yml | 22 ++++++++++++----------
 1 file changed, 12 insertions(+), 10 deletions(-)

diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml
index 1572af4b..15a402d0 100644
--- a/.github/workflows/build.yml
+++ b/.github/workflows/build.yml
@@ -2,6 +2,9 @@ name: Build
 
 on: [push, pull_request]
 
+defaults:
+  run:
+    shell: bash
 
 jobs:
   build:
@@ -12,7 +15,7 @@ jobs:
       matrix:
         os: [ macos-14, macos-15, macos-15-intel, ubuntu-24.04]
         compiler: [ gfortran ]
-        version: [ 13, 14, 15 ]
+        version: [ 13, 14, 15 ] # Julienne supports GFortran 13+
         extra_flags: [ -g -O3 ]
 
         exclude:
@@ -64,7 +67,6 @@ jobs:
       FC: ${{ matrix.compiler }}
       FFLAGS: ${{ matrix.extra_flags }}
       FPM_FLAGS:  --profile release --verbose
-      FPM_CXXFLAGS: -std=c++14 # needed for Armadillo on macos-14
 
     steps:
     - name: Checkout code
@@ -112,11 +114,7 @@ jobs:
         else \
           echo "FPM_FC=gfortran-${COMPILER_VERSION}" >> "$GITHUB_ENV" ; \
           echo "FFLAGS=-ffree-line-length-0 $FFLAGS" >> "$GITHUB_ENV" ; \
-        fi
-        if test -n "${{ matrix.error_stop_code }}" ; then \
-          echo "ERROR_STOP_CODE=${{ matrix.error_stop_code }}" >> "$GITHUB_ENV" ; \
-        else \
-          echo "ERROR_STOP_CODE=1" >> "$GITHUB_ENV" ; \
+          echo "FPM_CXXFLAGS=-std=c++14" >> "$GITHUB_ENV" ; : needed for Armadillo on macos-14 ; \
         fi
 
     - name: Setup FPM
@@ -136,19 +134,23 @@ jobs:
 
     - name: Version info
       run: |
-        ( type -p git > /dev/null && (echo == Commit Info == ; git log -n 1 ; echo ; echo ) ; exit 0 )
+        if test -d .git && type -p git > /dev/null ; then \
+          git config --global --add safe.directory `pwd` ; \
+          echo == Commit Info == ; \
+          git log -n 1 ; echo ; echo ; \
+        fi
         echo == Platform version info == 
         uname -a
         if test -r /etc/os-release ; then grep -e NAME -e VERSION /etc/os-release  ; fi
         if test -x /usr/bin/sw_vers ; then /usr/bin/sw_vers ; fi
         echo
         echo PATH="$PATH"
-        for tool in ${FPM_FC} fpm ; do
+        for tool in ${FPM_FC} ${FPM_CC} ${FPM_CXX} fpm ; do
           ( echo ; set -x ; w=$(which $tool) ; ls -al $w ; ls -alhL $w ; $tool --version )
         done
 
 
-    - name: Build and Test (Assertions OFF)
+    - name: Build and Test
       run: |
         set -x
         echo "FPM_FLAGS=${FPM_FLAGS}"

From bca86221e801a0ff9eed7f40b66dabbc0899e539 Mon Sep 17 00:00:00 2001
From: Dan Bonachea <dobonachea@lbl.gov>
Date: Thu, 16 Oct 2025 19:36:00 -0700
Subject: [PATCH 36/40] Add Intel Fortran coverage

This is complicated by the need for a compatible C++ compiler.

Intel provides two docker packages with Intel Fortran:

1.  The first is very compact, but lacks Intel C++.
    This can be used with GNU C/C++. This *probably* works in
    most cases, but might lead to subtle issues if there are
    ever ABI differences between the two tool-chains.

2.  The second includes Intel C++ and tons of other stuff,
    hence is much larger and can lead to slower container load.
    However it provides a "pure" Intel toolchain.

The CI scripts are adjusted to allow for either possibility.

Add a few configs of each type.
---
 .github/workflows/build.yml | 39 ++++++++++++++++++++++++++++++++++---
 1 file changed, 36 insertions(+), 3 deletions(-)

diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml
index 15a402d0..836ccc14 100644
--- a/.github/workflows/build.yml
+++ b/.github/workflows/build.yml
@@ -59,6 +59,32 @@ jobs:
             extra_flags: -g -mmlir -allow-assumed-rank -O3
             container: phhargrove/llvm-flang:19.1.1-latest
 
+          # --- Intel coverage ---
+
+          # Intel provides two docker packages with Intel Fortran.
+          # The first is very compact, but lacks Intel C++:
+          #   https://hub.docker.com/r/intel/fortran-essentials/tags
+          #   This can be used with GNU C/C++
+          # The second includes Intel C++ and tons of other stuff,
+          # hence is much larger and can lead to slower container load:
+          #   https://hub.docker.com/r/intel/oneapi-hpckit/tags
+
+          # Julienne requires 2025.2.0 or later
+          - os: ubuntu-24.04
+            compiler: ifx-gcc
+            version: latest
+            container: intel/fortran-essentials:latest
+
+          - os: ubuntu-24.04
+            compiler: ifx
+            version: 2025.2.2
+            container: intel/oneapi-hpckit:2025.2.2-0-devel-ubuntu24.04
+
+          - os: ubuntu-24.04
+            compiler: ifx-gcc
+            version: 2025.2.0
+            container: intel/fortran-essentials:2025.2.0-0-devel-ubuntu24.04
+
     container:
       image: ${{ matrix.container }}
 
@@ -104,10 +130,17 @@ jobs:
         set -x
         if test "$FC" = "flang" ; then \
           echo "FPM_FC=flang-new" >> "$GITHUB_ENV" ; \
-        elif test "$FC" = "ifx" ; then \
+        elif [[ "$FC" =~ "ifx" ]] ; then \
           echo "FPM_FC=ifx" >> "$GITHUB_ENV" ; \
-          echo "FFLAGS=-fpp $FFLAGS" >> "$GITHUB_ENV" ; \
-          echo "FFLAGS=-coarray $FFLAGS" >> "$GITHUB_ENV" ; \
+          echo "FFLAGS=-fpp -coarray $FFLAGS" >> "$GITHUB_ENV" ; \
+          : ls -al /opt/intel/oneapi/compiler/2025.*/bin/ ; \
+          if type -p icpx ; then \
+            echo "FPM_CC=icx" >> "$GITHUB_ENV" ; \
+            echo "FPM_CXX=icpx" >> "$GITHUB_ENV" ; \
+          else \
+            echo "FPM_CC=gcc" >> "$GITHUB_ENV" ; \
+            echo "FPM_CXX=g++" >> "$GITHUB_ENV" ; \
+          fi ; \
         elif test "$FC" = "lfortran" ; then \
           echo "FPM_FC=lfortran" >> "$GITHUB_ENV" ; \
           echo "FFLAGS=--cpp $FFLAGS" >> "$GITHUB_ENV" ; \

From d732f341fa7dfc8690f52c6cea65eb07bc88d117 Mon Sep 17 00:00:00 2001
From: Dan Bonachea <dobonachea@lbl.gov>
Date: Fri, 14 Nov 2025 14:29:34 -0800
Subject: [PATCH 37/40] CI: Fix FPM failures on macos-15-intel

---
 .github/workflows/build.yml | 9 ++++++++-
 1 file changed, 8 insertions(+), 1 deletion(-)

diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml
index 836ccc14..64e6785d 100644
--- a/.github/workflows/build.yml
+++ b/.github/workflows/build.yml
@@ -114,10 +114,17 @@ jobs:
         apt install -y build-essential pkg-config make git curl
 
     - name: Install macOS Dependencies
-      if: contains(matrix.os, 'macos') && matrix.compiler == 'flang'
+      if: contains(matrix.os, 'macos')
       run: |
         set -x
         brew update
+        # fpm binary distribution for macOS requires gfortran shared libraries from gcc@12
+        brew install gcc@12
+
+    - name: Install LLVM flang on macOS
+      if: contains(matrix.os, 'macos') && matrix.compiler == 'flang'
+      run: |
+        set -x
         brew install llvm@${COMPILER_VERSION} flang
         # workaround issue #228: clang cannot find homebrew flang's C header
         for p in /opt/homebrew /usr/local $(brew --prefix) ; do find $p/Cellar/flang -name ISO_Fortran_binding.h 2>/dev/null || true ; done

From 4940165a10902c4627b6ca0e9ae02e69c32336fb Mon Sep 17 00:00:00 2001
From: Dan Bonachea <dobonachea@lbl.gov>
Date: Fri, 14 Nov 2025 14:35:13 -0800
Subject: [PATCH 38/40] CI: Use Homebrew to install gfortran 15 on Ubuntu

---
 .github/workflows/build.yml | 27 ++++++++++++++++++++-------
 1 file changed, 20 insertions(+), 7 deletions(-)

diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml
index 64e6785d..12cf89ae 100644
--- a/.github/workflows/build.yml
+++ b/.github/workflows/build.yml
@@ -13,15 +13,11 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        os: [ macos-14, macos-15, macos-15-intel, ubuntu-24.04]
+        os: [ macos-14, macos-15, macos-15-intel, ubuntu-24.04 ]
         compiler: [ gfortran ]
         version: [ 13, 14, 15 ] # Julienne supports GFortran 13+
         extra_flags: [ -g -O3 ]
 
-        exclude:
-          - os: ubuntu-24.04
-            version: 15 # no package available (yet?)
-
         include:
 
           # --- LLVM flang coverage ---
@@ -99,12 +95,29 @@ jobs:
       uses: actions/checkout@v5
 
     - name: Install Ubuntu Native Dependencies
-      if: ${{ 0 && contains(matrix.os, 'ubuntu') && matrix.container == '' && matrix.compiler == 'gfortran' }}
+      if: ${{ contains(matrix.os, 'ubuntu') && matrix.container == '' && matrix.compiler == 'gfortran' }}
       run: |
         set -x
         sudo apt-get update
         sudo apt list -a 'gfortran-*'
-        sudo apt install -y gfortran-${COMPILER_VERSION} build-essential
+        sudo apt install -y build-essential
+        if [[ ${COMPILER_VERSION} < 15 ]] ; then \
+          sudo apt install -y gfortran-${COMPILER_VERSION} ; \
+        else \
+          curl -L https://raw.githubusercontent.com/Homebrew/install/HEAD/install.sh -o install-homebrew.sh ; \
+          chmod +x install-homebrew.sh ; \
+          env CI=1 ./install-homebrew.sh ; \
+          HOMEBREW_PREFIX="/home/linuxbrew/.linuxbrew" ; \
+          ${HOMEBREW_PREFIX}/bin/brew install -v gcc@${COMPILER_VERSION} binutils ; \
+          ls -al ${HOMEBREW_PREFIX}/bin ; \
+          echo "PATH=${HOMEBREW_PREFIX}/bin:${PATH}" >> "$GITHUB_ENV" ; \
+          : Homebrew GCC@15 needs binutils 2.44+ ; \
+          HOMEBREW_BINUTILS=$(ls -d ${HOMEBREW_PREFIX}/Cellar/binutils/2.*/bin ) ; \
+          ls -al ${HOMEBREW_BINUTILS} ; \
+          echo "FFLAGS=$FFLAGS -B ${HOMEBREW_BINUTILS}" >> "$GITHUB_ENV" ; \
+          echo "CFLAGS=$CFLAGS -B ${HOMEBREW_BINUTILS}" >> "$GITHUB_ENV" ; \
+          echo "CXXFLAGS=$CXXFLAGS -B ${HOMEBREW_BINUTILS}" >> "$GITHUB_ENV" ; \
+        fi
 
     - name: Install Ubuntu Container Dependencies
       if: ${{ contains(matrix.os, 'ubuntu') && matrix.container != '' && !contains(matrix.container, 'phhargrove') }}

From 96f3ef12e030e437bef1da2b76e2e9e42002ab9f Mon Sep 17 00:00:00 2001
From: Dan Bonachea <dobonachea@lbl.gov>
Date: Fri, 14 Nov 2025 18:12:19 -0800
Subject: [PATCH 39/40] CI: Fix eigen3 dependency on macOS that was breaking
 build-MOLE-macOSX

---
 .github/workflows/ci.yml | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 9e4d58a3..cb015c0c 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -75,7 +75,7 @@ jobs:
       - name: Install dependencies
         run: |
           brew update && (brew list cmake || brew install cmake)
-          brew install openblas eigen libomp lapack
+          brew install openblas eigen@3 libomp lapack
           sudo ln -sf /opt/homebrew/bin/gfortran-14 /opt/homebrew/bin/gfortran
         
       - name: Export Environment Variables
@@ -94,6 +94,8 @@ jobs:
           export CPPFLAGS+=" -I/opt/homebrew/opt/lapack/include"
           export PKG_CONFIG_PATH+=" /opt/homebrew/opt/lapack/lib/pkgconfig"
 
+          echo "Eigen3_DIR=/opt/homebrew/opt/eigen@3/share/eigen3/cmake" >> "$GITHUB_ENV"
+
       - name: Create build directory
         run: mkdir -p build
 

From 98892925c8a9c5a5c0efe0403af40e58fdf0c5ee Mon Sep 17 00:00:00 2001
From: Dan Bonachea <dobonachea@lbl.gov>
Date: Wed, 19 Nov 2025 15:37:44 -0800
Subject: [PATCH 40/40] CI: Incorporate recent platform changes

* Add the new macos-26 runner
* Add a workaround needed for recent snowstep/llvm containers
---
 .github/workflows/build.yml | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml
index 12cf89ae..1d5461ef 100644
--- a/.github/workflows/build.yml
+++ b/.github/workflows/build.yml
@@ -13,7 +13,7 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        os: [ macos-14, macos-15, macos-15-intel, ubuntu-24.04 ]
+        os: [ macos-14, macos-15, macos-15-intel, macos-26, ubuntu-24.04 ]
         compiler: [ gfortran ]
         version: [ 13, 14, 15 ] # Julienne supports GFortran 13+
         extra_flags: [ -g -O3 ]
@@ -169,6 +169,9 @@ jobs:
           echo "FFLAGS=-ffree-line-length-0 $FFLAGS" >> "$GITHUB_ENV" ; \
           echo "FPM_CXXFLAGS=-std=c++14" >> "$GITHUB_ENV" ; : needed for Armadillo on macos-14 ; \
         fi
+        if [[ "${{ matrix.container }}" =~ "snowstep/llvm" ]] ; then \
+          echo "LD_LIBRARY_PATH=/usr/lib/llvm-22/lib:$LD_LIBRARY_PATH" >> "$GITHUB_ENV" ; \
+        fi
 
     - name: Setup FPM
       uses: fortran-lang/setup-fpm@main