Add 1D divergence and laplacian operators

.github/workflows/build.yml

@@ -152,7 +152,7 @@ jobs:
           echo "FPM_FC=flang-new" >> "$GITHUB_ENV" ; \
         elif [[ "$FC" =~ "ifx" ]] ; then \
           echo "FPM_FC=ifx" >> "$GITHUB_ENV" ; \
-          echo "FFLAGS=-fpp -coarray $FFLAGS" >> "$GITHUB_ENV" ; \
+          echo "FFLAGS=-fpp -coarray -coarray-num-images=1 $FFLAGS" >> "$GITHUB_ENV" ; \
           : ls -al /opt/intel/oneapi/compiler/2025.*/bin/ ; \
           if type -p icpx ; then \
             echo "FPM_CC=icx" >> "$GITHUB_ENV" ; \

doc/fortran-classes.md

@@ -0,0 +1,51 @@
+MOLE Fortran Class Diagram
+--------------------------
+
+```mermaid
+
+%%{init: { 'theme':'neo',  "class" : {"hideEmptyMembersBox": true} } }%%
+
+classDiagram
+
+class tensor_1D_t
+class scalar_1D_t
+class vector_1D_t
+class divergence_1D_t
+class laplacian_1D_t
+class gradient_operator_1D_t
+class divergence_operator_1D_t
+class mimetic_matrix_1D_t
+
+tensor_1D_t <|-- scalar_1D_t : is a
+tensor_1D_t <|-- vector_1D_t : is a
+tensor_1D_t <|-- divergence_1D_t : is a
+divergence_1D_t <|-- laplacian_1D_t : is a
+mimetic_matrix_1D_t <|-- gradient_operator_1D_t : is a
+mimetic_matrix_1D_t <|-- divergence_operator_1D_t : is a
+
+class scalar_1D_t{
+    - gradient_operator_1D_ : gradient_operator_1D_t
+    + operator(.grad.) vector_1D_t
+    + operator(.laplacian.) scalar_1D_t
+}
+
+class vector_1D_t{
+    - divergence_operator_1D_ : divergence_operator_1D_t
+    + operator(.div.) divergence_1D_t
+}
+
+class mimetic_matrix_1D_t{
+ - upper_ :: double precision
+ - inner_ :: double precision
+ - lower_ :: double precision
+}
+
+class gradient_operator_1D_t{
+ + operator(.x.) double precision[]
+ + assemble() double precision[] "2D array"
+}
+
+class divergence_operator_1D_t{
+ + operator(.x.) double precision[]
+ + assemble() double precision[] "2D array"
+}

example/div-grad-laplacian-1D.F90

@@ -0,0 +1,133 @@
+module functions_m
+  implicit none
+
+contains
+
+  pure function f(x)
+    double precision, intent(in) :: x(:)
+    double precision, allocatable :: f(:)
+    f = (x**3)/6 + (x**2)/2 + 1
+  end function
+
+  double precision elemental function df_dx(x)
+    double precision, intent(in) :: x
+    df_dx = (x**2)/2 + x
+  end function
+
+  double precision elemental function d2f_dx2(x)
+    double precision, intent(in) :: x
+    d2f_dx2 = x + 1
+  end function
+
+end module functions_m
+
+program div_grad_laplacian_1D
+  use functions_m
+  use julienne_m, only :  file_t, string_t, operator(.separatedBy.)
+  use mole_m, only : scalar_1D_t, scalar_1D_initializer_i
+#ifdef __GFORTRAN__
+  use mole_m, only : vector_1D_t, laplacian_1D_t
+#endif
+  implicit none
+
+  procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => f
+
+  print *,new_line('')
+  print *,"            2nd-order approximations"
+  print *,"            ========================"
+
+  call output(order=2)
+
+  print *,new_line('')
+  print *,"            4th-order approximations"
+  print *,"            ========================"
+
+  call output(order=4)
+
+contains
+
+#ifndef __GFORTRAN__
+
+  subroutine output(order)
+    integer, intent(in) :: order
+
+    associate(   s           => scalar_1D_t(scalar_1D_initializer, order=order, cells=10, x_min=0D0, x_max=20D0))
+      associate( grad_s      => .grad. s &
+                ,laplacian_s => .laplacian. s)
+        associate( s_grid           => s%grid()      &
+                  ,grad_s_grid      => grad_s%grid() &
+                  ,laplacian_s_grid => laplacian_s%grid())
+          associate( s_table           => tabulate( &
+                        string_t([character(len=18)::"x", "f(x) exp"         , "f(x) act"         ]) &
+                       ,s_grid, f(s_grid), s%values() &
+                     ) &
+                    ,grad_s_table      => tabulate( &
+                       string_t([character(len=18)::"x", ".grad. f exp"     , ".grad. f act"     ])  &
+                      ,grad_s_grid, df_dx(grad_s_grid), grad_s%values() &
+                     ) &
+                    ,laplacian_s_table => tabulate( &
+                       string_t([character(len=18)::"x", ".laplacian. f exp", ".laplacian. f act"])  &
+                      ,laplacian_s_grid, d2f_dx2(laplacian_s_grid), laplacian_s%values()) &
+                     )
+             call s_table%write_lines()
+             call grad_s_table%write_lines()
+             call laplacian_s_table%write_lines()
+          end associate
+        end associate
+      end associate
+    end associate
+  end subroutine
+
+#else
+
+  subroutine output(order)
+    integer, intent(in) :: order
+
+    type(scalar_1D_t) s
+    type(vector_1D_t) grad_s
+    type(laplacian_1D_t) laplacian_s
+    type(file_t) s_table, grad_s_table, laplacian_s_table
+    double precision, allocatable,dimension(:) :: s_grid, grad_s_grid, laplacian_s_grid
+
+    s = scalar_1D_t(scalar_1D_initializer, order=order, cells=10, x_min=0D0, x_max=20D0)
+    grad_s = .grad. s
+    laplacian_s = .laplacian. s
+
+    s_grid = s%grid()
+    grad_s_grid = grad_s%grid()
+    laplacian_s_grid = laplacian_s%grid()
+
+    s_table           = tabulate( &
+       string_t([character(len=18)::"x", "f(x) exp."         , "f(x) act."         ]) &
+      ,s_grid, f(s_grid), s%values() &
+    )
+    grad_s_table      = tabulate( &
+       string_t([character(len=18)::"x", ".grad. f exp."     , ".grad. f act."     ]) &
+      ,grad_s_grid, df_dx(grad_s_grid), grad_s%values() &
+    )
+    laplacian_s_table = tabulate( &
+       string_t([character(len=18)::"x", ".laplacian. f exp.", ".laplacian. f act."]) &
+      ,laplacian_s_grid, d2f_dx2(laplacian_s_grid), laplacian_s%values() &
+    )
+    call s_table%write_lines()
+    call grad_s_table%write_lines()
+    call laplacian_s_table%write_lines()
+  end subroutine
+
+#endif
+
+  pure function tabulate(headings, abscissa, expected, actual) result(file)
+    double precision, intent(in), dimension(:) :: abscissa, expected, actual
+    type(string_t), intent(in) :: headings(:)
+    type(file_t) file
+    integer line
+
+    file = file_t([ &
+       string_t("") &
+      ,headings .separatedBy. "  " &
+      ,string_t("----------------------------------------------------------") &
+      ,[( string_t(abscissa(line)) // "      " // string_t(expected(line)) // "        " // string_t(actual(line)), line = 1, size(abscissa))] &
+    ])
+  end function
+
+end program

example/print-assembled-1D-operators.f90

@@ -0,0 +1,78 @@
+program print_assembled_1D_operators
+  !! Print fully assembled memetic 1D gradient, divergence, and Laplacian matrices,
+  !! including the zero elements.
+  use julienne_m, only : operator(.csv.), string_t, command_line_t
+  use mimetic_operators_1D_m, only : gradient_operator_1D_t, divergence_operator_1D_t
+  implicit none
+
+  type(command_line_t) command_line
+  integer row
+
+  command_line_settings: &
+  associate( &
+     gradient   => command_line%argument_present(["--grad" ]) &
+    ,divergence => command_line%argument_present(["--div"  ]) &
+    ,order      => command_line%flag_value("--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 print-assembled-1D-operators \' // new_line('') &
+        // '  --compiler flang-new \'                   // new_line('') &
+        // '  --flag "-O3" \'                           // new_line('') & 
+        // '  -- [--help|-h] | [--grad]  [--div] [--order <integer>]' // new_line('') // new_line('') &
+        // 'where square brackets indicate optional arguments and angular brackets indicate user input values.' // new_line('')
+    end if
+
+    default_usage: &
+    associate(print_all => .not. any([gradient, divergence, len(order)/=0]))
+
+      if (print_all .or. (gradient   .and. len(order)==0) .or. (gradient   .and. order=="2")) call print_gradient_operator(  k=2, dx=1D0, m=5)
+      if (print_all .or. (divergence .and. len(order)==0) .or. (divergence .and. order=="2")) call print_divergence_operator(k=2, dx=1D0, m=5)
+      if (print_all .or. (gradient   .and. len(order)==0) .or. (gradient   .and. order=="4")) call print_gradient_operator(  k=4, dx=1D0, m=9)
+      if (print_all .or. (divergence .and. len(order)==0) .or. (divergence .and. order=="4")) call print_divergence_operator(k=4, dx=1D0, m=9)
+
+    end associate default_usage
+  end associate command_line_settings
+
+contains
+
+  subroutine print_gradient_operator(k, dx, m)
+    integer, intent(in) :: k, m 
+    double precision, intent(in) :: dx
+
+    print *, new_line(""), "Gradient operator: order = ", k, " | cells = ", m, " | dx = ", dx
+
+    associate(grad_op => gradient_operator_1D_t(k, dx, cells=m))
+      associate(G => grad_op%assemble())
+        do row = 1, size(G,1)
+          associate(csv_row => .csv. string_t(G(row,:)))
+            print '(a)', csv_row%string()
+          end associate
+        end do
+      end associate
+    end associate
+
+  end subroutine
+
+  subroutine print_divergence_operator(k, dx, m)
+    integer, intent(in) :: k, m 
+    double precision, intent(in) :: dx
+
+    print *, new_line(""), "Divergence operator: order = ", k, " | cells = ", m, " | dx = ", dx
+
+    associate(div_op => divergence_operator_1D_t(k, dx, cells=m))
+      associate(D => div_op%assemble())
+        do row = 1, size(D,1)
+          associate(csv_row => .csv. string_t(D(row,:)))
+            print '(a)', csv_row%string()
+          end associate
+        end do
+      end associate
+    end associate
+
+  end subroutine
+
+end program

fpm.toml

@@ -4,7 +4,7 @@ name = "MOLE"
 armadillo-code = {git = "https://gitlab.com/rouson/armadillo-code.git", tag = "fpm"}
 
 [dev-dependencies]
-julienne = {git = "https://github.com/berkeleylab/julienne.git", tag = "3.1.5"}
+julienne = {git = "https://github.com/berkeleylab/julienne.git", tag = "3.3.0"}
 
 [install]
 library = true

src/fortran/divergence_1D_s.F90

@@ -0,0 +1,33 @@
+submodule(tensors_1D_m) divergence_1D_s
+  implicit none
+
+contains
+
+#ifdef __GFORTRAN__
+
+  pure function cell_center_locations(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 => (x_max - x_min)/cells)
+      x = x_min + dx/2. + [((cell-1)*dx, cell = 1, cells)]
+    end associate
+  end function
+
+#endif
+
+  module procedure construct_from_tensor
+    divergence_1D%tensor_1D_t = tensor_1D
+  end procedure
+
+  module procedure divergence_1D_values
+    cell_centered_values = self%values_
+  end procedure
+
+  module procedure divergence_1D_grid
+    cell_centers = cell_center_locations(self%x_min_, self%x_max_, self%cells_)
+  end procedure
+
+end submodule divergence_1D_s