Add Burgers Equation solver example

README.md

@@ -84,11 +84,11 @@ Building and testing
 ### `fpm` versions before 0.13.0
 With LLVM 20-22, replace the above `flang` command with
 ```
-  fpm test --compiler flang-new --flag "-O3"
+  fpm test --compiler flang --flag "-O3"
 ```
 With LLVM 19, replace the above `flang` command with
 ```
- fpm test --compiler flang-new --flag "-O3 -mmlir -allow-assumed-rank"
+ fpm test --compiler flang --flag "-O3 -mmlir -allow-assumed-rank"
 ```
 
 Documentation

example/burgers-1D.F90

@@ -0,0 +1,97 @@
+! 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(:)
+    initial_condition = 10*sin(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
+  !! This program demonstrates the use of Formal to solve the partial differential equation of
+  !! Burgers (1948) in conservative form using the 2nd- or 4th-order mimetic discretizations of
+  !! Corbino & Castillo (2020) and Dumett & Castillo (2022).
+  !!
+  !! * Burgers, J.M.      (1948) https://doi.org/10.1016/S0065-2156(08)70100-5
+  !! * Corbino & Castillo (2020) https://doi.org/10.1016/j.cam.2019.06.042.
+  !! * Dumett & Castillo  (2022) https://doi.org/10.13140/RG.2.2.26630.14400
+  use initial_condition_m, only : initial_condition
+  use julienne_m, only :  command_line_t
+  use formal_m, only : scalar_1D_t, scalar_1D_initializer_i, d_dx, d2_dx2
+  implicit none
+
+#ifdef __GFORTRAN__
+  procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer
+#endif
+  character(len=:), allocatable :: order_string
+  type(command_line_t) command_line
+  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 = 2
+  else
+    read(order_string,"(i1)") order 
+  end if
+
+  print *, "order = ", order
+
+  block
+#ifndef __GFORTRAN__
+    procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer
+#endif
+    double precision, parameter :: dt = 1D0, pi = acos(-1D0)
+    type(scalar_1D_t) u
+
+    scalar_1D_initializer => initial_condition
+    u = scalar_1D_t(scalar_1D_initializer, order, x_min=0D0, x_max=2*pi, cells=20)
+
+    runge_kutta_2nd_order: &
+    associate(u_half => u + d_dt(u)*dt/2) ! first substep
+      u = u + d_dt(u_half)*dt ! second substep
+    end associate runge_kutta_2nd_order
+
+  end block
+
+#ifdef __GFORTRAN__
+    stop
+#endif
+
+contains
+
+  pure function d_dt(u) result(du_dt)
+    type(scalar_1D_t), intent(in) :: u
+    type(scalar_1D_t) du_dt
+    double precision, parameter :: nu=1D0
+    du_dt = nu*d2_dx2(u) - d_dx((u**2)/2)
+  end function
+
+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/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