Feature: Interpolation

example/burgers-1D.F90

@@ -0,0 +1,75 @@
+! Copyright (c) 2026, The Regents of the University of California
+! Terms of use are as specified in LICENSE.txt
+
+module initial_condition_m
+  implicit none
+
+contains
+
+  pure function initial_condition(x)
+    !! Initial solution to Burgers equation
+    double precision, intent(in) :: x(:)
+    double precision, allocatable :: initial_condition(:)
+    double precision, parameter :: pi = acos(-1D0)
+    initial_condition = sin(2*pi*x)
+    ! To change this function, please edit only the right-hand-side (RHS) expression,
+    ! keeping the rest in place for proper display of the function at runtime.
+  end function
+
+end module
+
+program burgers_1D
+  !! Advance the 1D Burgers partial differential equation over time.
+  use initial_condition_m, only : initial_condition
+  use julienne_m, only :  command_line_t
+  use formal_m, only : vector_1D_t, vector_1D_initializer_i
+  implicit none
+
+  procedure(vector_1D_initializer_i), pointer :: vector_1D_initializer => initial_condition
+  character(len=:), allocatable :: order_string
+  type(command_line_t) command_line
+  type(vector_1D_t) :: u
+  integer order
+
+  if (command_line%argument_present([character(len=len("--help")) :: ("--help"), "-h"])) then
+    stop                  new_line('') // new_line('') &
+      // 'Usage:'                      // new_line('') &
+      // '  fpm run \'                 // new_line('') &
+      // '    --example burgers-1D \'  // new_line('') &
+      // '    --compiler flang \'      // new_line('') &
+      // '    --flag "-O3" \'          // new_line('') & 
+      // '    [--help|-h] | [--order <integer>]' // new_line('') // new_line('') &
+      // 'where square brackets indicate optional arguments and angular brackets indicate user input values.' // new_line('')
+  end if
+
+  print *, new_line('')
+  print *,"   Initial condition"
+  print *,"   ================="
+
+  call execute_command_line("grep 'initial_condition =' example/burgers-1D.F90 | grep -v execute_command", wait=.true.)
+
+  order_string = command_line%flag_value("--order")
+
+  if (len(order_string)==0) then 
+    order = 4 
+  else
+    read(order_string,"(i1)") order 
+  end if
+
+  print *, "order = ", order
+
+  u = vector_1D_t(vector_1D_initializer, order, x_min=0D0, x_max=1D0, cells=50)
+
+  block
+    double precision dt
+    dt = 1D0
+    !associate(dt_grad_u2_2 => dt * (.div. ((u*u))))
+    associate(dt_div_uu => .div. (u*u))
+    end associate
+  end block
+
+#ifdef __GFORTRAN__
+    stop
+#endif
+
+end program

src/formal/mimetic_matrix_1D_s.F90 → src/formal/differential_operator_matrix_1D_s.F90

@@ -1,16 +1,16 @@
 ! Copyright (c) 2026, The Regents of the University of California
 ! Terms of use are as specified in LICENSE.txt
 
-submodule(mimetic_operators_1D_m) mimetic_matrix_1D_s
+submodule(differential_operators_1D_m) differential_operator_matrix_1D_s
   use julienne_m, only : string_t, operator(.csv.)
   implicit none
 
 contains
 
   module procedure construct_matrix_operator
-    mimetic_matrix_1D%upper_ = upper
-    mimetic_matrix_1D%inner_ = inner
-    mimetic_matrix_1D%lower_ = lower
+    differential_operator_matrix_1D%upper_ = upper
+    differential_operator_matrix_1D%inner_ = inner
+    differential_operator_matrix_1D%lower_ = lower
   end procedure
 
   module procedure to_file_t
@@ -33,4 +33,4 @@
 
   end procedure
 
-end submodule mimetic_matrix_1D_s
+end submodule differential_operator_matrix_1D_s

src/formal/mimetic_operators_1D_m.F90 → src/formal/differential_operators_1D_m.F90

@@ -3,7 +3,7 @@
 
 #include "formal-language-support.F90"
 
-module mimetic_operators_1D_m
+module differential_operators_1D_m
   !! Define sparse matrix storage formats and operators tailored to the one-dimensional (1D) mimetic discretizations
   !! detaild by Corbino & Castillo (2020) https://doi.org/10.1016/j.cam.2019.06.042.
   use julienne_m, only : file_t
@@ -13,8 +13,8 @@ module mimetic_operators_1D_m
   public :: gradient_operator_1D_t
   public :: divergence_operator_1D_t
 
-  type mimetic_matrix_1D_t
-    !! Encapsulate a mimetic matrix with a corresponding matrix-vector product operator
+  type differential_operator_matrix_1D_t
+    !! Encapsulate a mimetic matrix
     private
     double precision, allocatable :: upper_(:,:) !! A  submatrix block (cf. Corbino & Castillo, 2020)
     double precision, allocatable :: inner_(:)   !! M  submatrix row   (cf. Corbino & Castillo, 2020)
@@ -23,20 +23,20 @@ module mimetic_operators_1D_m
     procedure, non_overridable :: to_file_t
   end type
 
-  interface mimetic_matrix_1D_t
+  interface differential_operator_matrix_1D_t
 
-    pure module function construct_matrix_operator(upper, inner, lower) result(mimetic_matrix_1D)
+    pure module function construct_matrix_operator(upper, inner, lower) result(differential_operator_matrix_1D)
       !! Construct discrete operator from matrix blocks
       implicit none
       double precision, intent(in) :: upper(:,:) !! A  submatrix block (cf. Corbino & Castillo, 2020)
       double precision, intent(in) :: inner(:)   !! M  submatrix row   (cf. Corbino & Castillo, 2020)
       double precision, intent(in) :: lower(:,:) !! A' submatrix block (cf. Corbino & Castillo, 2020)
-      type(mimetic_matrix_1D_t) mimetic_matrix_1D
+      type(differential_operator_matrix_1D_t) differential_operator_matrix_1D
     end function
 
   end interface
 
-  type, extends(mimetic_matrix_1D_t) :: gradient_operator_1D_t
+  type, extends(differential_operator_matrix_1D_t) :: gradient_operator_1D_t
     !! Encapsulate a 1D mimetic gradient operator
     private
     integer k_ !! order of accuracy
@@ -62,7 +62,7 @@ pure module function construct_1D_gradient_operator(k, dx, cells) result(gradien
 
   end interface
 
-  type, extends(mimetic_matrix_1D_t) :: divergence_operator_1D_t
+  type, extends(differential_operator_matrix_1D_t) :: divergence_operator_1D_t
     !! Encapsulate kth-order mimetic divergence operator on m_ cells of width dx
     private
     integer k_, m_
@@ -129,7 +129,7 @@ pure module function divergence_matrix_multiply(self, vec) result(matvec_product
 
      pure module function to_file_t(self) result(file)
        implicit none
-       class(mimetic_matrix_1D_t), intent(in) :: self
+       class(differential_operator_matrix_1D_t), intent(in) :: self
        type(file_t) file
      end function
 
@@ -164,4 +164,4 @@ pure function negate_and_flip(A) result(Ap)
 
 #endif
 
-end module mimetic_operators_1D_m
+end module differential_operators_1D_m

src/formal/divergence_operator_1D_s.F90

@@ -4,7 +4,7 @@
 #include "formal-language-support.F90"
 #include "julienne-assert-macros.h"
 
-submodule(mimetic_operators_1D_m) divergence_operator_1D_s
+submodule(differential_operators_1D_m) divergence_operator_1D_s
   use julienne_m, only : call_julienne_assert_, string_t
 #if ASSERTIONS
   use julienne_m, only : operator(.isAtLeast.), operator(.equalsExpected.)
@@ -48,7 +48,7 @@ pure function negate_and_flip(A) result(Ap)
       else
         allocate(Ap, mold = A)
       end if
-      divergence_operator_1D%mimetic_matrix_1D_t = mimetic_matrix_1D_t(A, M(k, dx), Ap)
+      divergence_operator_1D%differential_operator_matrix_1D_t = differential_operator_matrix_1D_t(A, M(k, dx), Ap)
       divergence_operator_1D%k_  = k
       divergence_operator_1D%dx_ = dx
       divergence_operator_1D%m_  = cells

src/formal/dyad_1D_s.F90

@@ -0,0 +1,66 @@
+! Copyright (c) 2026, The Regents of the University of California
+! Terms of use are as specified in LICENSE.txt
+
+#include "julienne-assert-macros.h"
+
+submodule(tensors_1D_m) dyad_1D_s
+  use julienne_m, only : &
+     call_julienne_assert_ &
+    ,operator(.all.) &
+    ,operator(.approximates.) &
+    ,operator(.equalsExpected.) &
+    ,operator(.within.)
+  use interpolator_1D_m, only : centers_to_faces_1D_t
+  implicit none
+
+  double precision, parameter :: double_equivalence = 2D-4
+
+contains
+
+  module procedure dyad_over_integer
+    ratio%tensor_1D_t = tensor_1D_t(self%values_/numerator, self%x_min_, self%x_max_, self%cells_, order = self%order_)
+    ratio%divergence_operator_1D_ = self%divergence_operator_1D_
+  end procedure
+
+  module procedure construct_1D_dyad_from_components
+    dyad_1D%tensor_1D_t = tensor_1D
+    dyad_1D%divergence_operator_1D_ = divergence_operator_1D
+  end procedure
+
+  module procedure div_dyad
+
+    integer center
+
+#ifdef NAGFOR
+    associate(D => self%divergence_operator_1D_)
+#else
+    associate(D => (self%divergence_operator_1D_))
+#endif
+      associate( &
+         Dv => D .x. self%values_ &
+        ,dx => (self%x_max_ - self%x_min_)/self%cells_ &
+      )
+        associate(interpolator => centers_to_faces_1D_t(order=self%order_, cells=self%cells_, dx=dx))
+          vector_1D = vector_1D_t( &
+              tensor_1D_t(interpolator%face_values(Dv(2:size(Dv)-1)), self%x_min_, self%x_max_, self%cells_, self%order_) &
+            ,divergence_operator_1D_t(k=self%order_, dx=dx, cells=self%cells_) &
+          )
+        end associate
+#if ASSERTIONS
+        associate(divergence_1D => divergence_1D_t(tensor_1D_t(Dv(2:size(Dv)-1), self%x_min_, self%x_max_, self%cells_, self%order_)))
+          associate( &
+             q => divergence_1D%weights() &
+            ,b => [-1D0, [(0D0, center = 1, self%cells_-1)], 1D0] &
+          )
+            call_julienne_assert(.all. ([size(Dv), size(q)] .equalsExpected. self%cells_+2))
+            call_julienne_assert((.all. (matmul(transpose(D%assemble()), q) .approximates. b/dx .within. double_equivalence)))
+              ! Check D^T * a = b_{m+1},  Eq. (19), Corbino & Castillo (2020)
+          end associate
+        end associate
+#endif
+      end associate
+    end associate
+
+  end procedure
+
+end submodule dyad_1D_s

src/formal/gradient_operator_1D_s.F90

@@ -4,7 +4,7 @@
 #include "julienne-assert-macros.h"
 #include "formal-language-support.F90"
 
-submodule(mimetic_operators_1D_m) gradient_operator_1D_s
+submodule(differential_operators_1D_m) gradient_operator_1D_s
   use julienne_m, only : call_julienne_assert_, string_t
 #if ASSERTIONS
   use julienne_m, only : operator(.isAtLeast.)
@@ -42,7 +42,7 @@ pure function negate_and_flip(A) result(Ap)
     call_julienne_assert(cells .isAtLeast. 2*k)
 
     associate(A => corbino_castillo_A(k, dx), M => corbino_castillo_M(k, dx))
-      gradient_operator_1D%mimetic_matrix_1D_t = mimetic_matrix_1D_t(A, M, negate_and_flip(A))
+      gradient_operator_1D%differential_operator_matrix_1D_t = differential_operator_matrix_1D_t(A, M, negate_and_flip(A))
       gradient_operator_1D%k_  = k
       gradient_operator_1D%dx_ = dx
       gradient_operator_1D%m_  = cells

src/formal/interpolator_1D_m.F90

@@ -0,0 +1,78 @@
+! Copyright (c) 2026, The Regents of the University of California
+! Terms of use are as specified in LICENSE.txt
+
+module interpolator_1D_m
+  !! Define a sparse matrix storage format tailored to the staggered-grid interpolation matrix
+  !! operators detailed by Dumett & Castillo (2022) https://doi.org/10.13140/RG.2.2.26630.14400
+  implicit none
+
+  private
+  public :: centers_to_faces_1D_t
+  public :: faces_to_centers_1D_t
+
+  type, abstract :: interpolator_1D_t
+    !! Encapsulate a staggered-grid interpolation matrix with a corresponding matrix-vector product operator
+    private
+    integer order_, cells_, dx_
+    double precision                 first_
+    double precision, allocatable :: upper_(:,:) 
+    double precision, allocatable :: inner_(:)
+    double precision, allocatable :: lower_(:,:)
+    double precision                 final_
+  end type
+
+  type, extends(interpolator_1D_t) :: centers_to_faces_1D_t
+  contains
+    procedure, non_overridable :: face_values
+  end type
+
+  type, extends(interpolator_1D_t) :: faces_to_centers_1D_t
+  contains
+    procedure, non_overridable :: center_values
+  end type
+
+  interface centers_to_faces_1D_t
+
+    pure module function c2f_constructor(order, cells, dx) result(centers_to_faces_1D)
+      !! Construct centers-to-faces interpolation operator
+      implicit none
+      integer, intent(in) :: order, cells
+      double precision, intent(in) :: dx
+      type(centers_to_faces_1D_t) centers_to_faces_1D
+    end function
+
+  end interface
+
+  interface faces_to_centers_1D_t
+
+    pure module function f2c_constructor(order, cells, dx) result(faces_to_centers_1D)
+      !! Construct centers-to-faces interpolation operator
+      implicit none
+      integer, intent(in) :: order, cells
+      double precision, intent(in) :: dx
+      type(faces_to_centers_1D_t) faces_to_centers_1D
+    end function
+
+  end interface
+
+  interface
+
+    pure module function face_values(self, centers_extended) result(faces)
+      !! Interpolate cell-centered values to face-centered values
+      implicit none
+      class(centers_to_faces_1D_t), intent(in) :: self
+      double precision, intent(in) :: centers_extended(:)
+      double precision, allocatable :: faces(:)
+    end function
+
+    pure module function center_values(self, faces) result(centers_extended)
+      !! Interpolate face-centered values to cell-centered values
+      implicit none
+      class(faces_to_centers_1D_t), intent(in) :: self
+      double precision, intent(in) :: faces(:)
+      double precision, allocatable :: centers_extended(:)
+    end function
+
+  end interface
+
+end module interpolator_1D_m