Add quadrature

doc/fortran-classes.md

@@ -10,6 +10,7 @@ classDiagram
 class tensor_1D_t
 class scalar_1D_t
 class vector_1D_t
+class gradient_1D_t
 class divergence_1D_t
 class laplacian_1D_t
 class gradient_operator_1D_t
@@ -19,25 +20,35 @@ 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
+
+tensor_1D_t <|-- weighted_product_1D_t
+tensor_1D_t <|-- vector_dot_gradient_1D_t
+tensor_1D_t <|-- scalar_x_divergence_1D_t
+
 divergence_1D_t <|-- laplacian_1D_t : is a
+vector_1D_t <|-- gradient_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
+    + operator(.grad.) gradient_1D_t
+    + operator(.laplacian.) laplacian_1D_t
 }
 
 class vector_1D_t{
     - divergence_operator_1D_ : divergence_operator_1D_t
     + operator(.div.) divergence_1D_t
 }
 
+class gradient_1D_t{
+   - weights : double precision[]
+}
+
 class mimetic_matrix_1D_t{
- - upper_ :: double precision
- - inner_ :: double precision
- - lower_ :: double precision
+ - upper_ :: double precision[]
+ - inner_ :: double precision[]
+ - lower_ :: double precision[]
 }
 
 class gradient_operator_1D_t{

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

@@ -22,25 +22,34 @@ double precision elemental function d2f_dx2(x)
 end module functions_m
 
 program div_grad_laplacian_1D
+  !! Compute the 2nd- and 4th-order mimetic approximations to the gradient and Laplacian of the
+  !! above function f(x) on a 1D uniform, staggered grid.
   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
+  use mole_m, only : vector_1D_t, laplacian_1D_t, gradient_1D_t
 #endif
   implicit none
 
   procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => f
 
   print *,new_line('')
-  print *,"            2nd-order approximations"
-  print *,"            ========================"
+  print *,"   Functions"
+  print *,"   ========================"
+  call execute_command_line("grep 'f =' example/div-grad-laplacian-1D.F90 | grep -v execute_command", wait=.true.)
+  call execute_command_line("grep 'df_dx =' example/div-grad-laplacian-1D.F90 | grep -v execute_command", wait=.true.)
+  call execute_command_line("grep 'd2f_dx2 =' example/div-grad-laplacian-1D.F90 | grep -v execute_command", wait=.true.)
+
+  print *,new_line('')
+  print *,"   2nd-order approximations"
+  print *,"   ========================"
 
   call output(order=2)
 
   print *,new_line('')
-  print *,"            4th-order approximations"
-  print *,"            ========================"
+  print *,"   4th-order approximations"
+  print *,"   ========================"
 
   call output(order=4)
 
@@ -51,22 +60,22 @@ program div_grad_laplacian_1D
   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(   s           => scalar_1D_t(scalar_1D_initializer, order=order, cells=20, 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"         ]) &
+                        string_t([character(len=22)::"x", "f(x) expected"         , "f(x) actual"         ]) &
                        ,s_grid, f(s_grid), s%values() &
                      ) &
                     ,grad_s_table      => tabulate( &
-                       string_t([character(len=18)::"x", ".grad. f exp"     , ".grad. f act"     ])  &
+                       string_t([character(len=22)::"x", ".grad. f expected"     , ".grad. f actual"     ])  &
                       ,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"])  &
+                       string_t([character(len=22)::"x", ".laplacian. f expected", ".laplacian. f actual"])  &
                       ,laplacian_s_grid, d2f_dx2(laplacian_s_grid), laplacian_s%values()) &
                      )
              call s_table%write_lines()
@@ -84,12 +93,12 @@ subroutine output(order)
     integer, intent(in) :: order
 
     type(scalar_1D_t) s
-    type(vector_1D_t) grad_s
+    type(gradient_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)
+    s = scalar_1D_t(scalar_1D_initializer, order=order, cells=20, x_min=0D0, x_max=20D0)
     grad_s = .grad. s
     laplacian_s = .laplacian. s
 
@@ -98,15 +107,15 @@ subroutine output(order)
     laplacian_s_grid = laplacian_s%grid()
 
     s_table           = tabulate( &
-       string_t([character(len=18)::"x", "f(x) exp."         , "f(x) act."         ]) &
+       string_t([character(len=22)::"x", "f(x) expected"         , "f(x) actual"         ]) &
       ,s_grid, f(s_grid), s%values() &
     )
     grad_s_table      = tabulate( &
-       string_t([character(len=18)::"x", ".grad. f exp."     , ".grad. f act."     ]) &
+       string_t([character(len=22)::"x", ".grad. f expected"     , ".grad. f actual"     ]) &
       ,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."]) &
+       string_t([character(len=22)::"x", ".laplacian. f expected", ".laplacian. f actual"]) &
       ,laplacian_s_grid, d2f_dx2(laplacian_s_grid), laplacian_s%values() &
     )
     call s_table%write_lines()
@@ -125,8 +134,8 @@ pure function tabulate(headings, abscissa, expected, actual) result(file)
     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))] &
+      ,string_t("------------------------------------------------------------------") &
+      ,[( string_t(abscissa(line)) // "          " // string_t(expected(line)) // "          " // string_t(actual(line)), line = 1, size(abscissa))] &
     ])
   end function
 

example/extended-gauss-divergence.F90

@@ -0,0 +1,149 @@
+module integrand_operands_m
+  implicit none
+contains
+
+  pure function scalar(x) result(f)
+    double precision, intent(in) :: x(:)
+    double precision, allocatable :: f(:)
+    f = (x**2)/2 ! <-- scalar function
+  end function
+
+  pure function vector(x) result(v)
+    double precision, intent(in) :: x(:)
+    double precision, allocatable :: v(:)
+    v = x        ! <-- vector function
+  end function
+
+end module
+
+program extended_gauss_divergence
+  !! Print each term in the following residual formed from the extended Gauss-divergence
+  !! theorem using one-dimensional (1D) 4th-(default) or 2nd-order mimetic discretizations:
+  !! `residual = .SSS. (v .dot. .grad. f) * dV +.SSS. (f * .div. v) * dV - .SS. (f .x. (v .dot. dA))`
+  !! where `.SSS.` and `.SS.` are the 1D equivalents of a volume integral over the whole
+  !! domain and a surface integral over a domain boundary of unit area, respectively.
+  use julienne_m, only : command_line_t
+  use integrand_operands_m, only : scalar, vector
+  use mole_m, only : scalar_1D_t, scalar_1D_initializer_i, vector_1D_t, vector_1D_initializer_i
+  implicit none
+  procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => scalar
+  procedure(vector_1D_initializer_i), pointer :: vector_1D_initializer => vector
+
+  type numerical_arguments_t
+    !! Define default initializations that can be overridden with the command-line arguments
+    !! detailed by the usage information below
+    integer :: cells_=200, order_=4
+    double precision :: x_min_=0D0, x_max_=1D0
+  end type
+
+  type text_flags_t
+    logical div_, grad_, vf_
+  end type
+
+  type(command_line_t) command_line
+  double precision SSS_v_dot_grad_f_dV, SSS_f_div_v_dV, SS_f_v_dot_dA
+
+  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 extended-gauss-divergence \' // new_line('') &
+      // '  --compiler flang-new \'                // new_line('') &
+      // '  --flag "-O3" \'                        // new_line('') &
+      // '  -- [--help|-h] | [[--cells <integer>] [--order <integer>] [--xmin <double precision>] [--xmax <double precision>] [--div|d] [--grad|g] [--vf|f]]' &
+      //                              new_line('') // new_line('') &
+      // 'where pipes (|) separate square-bracketed optional arguments and angular brackets indicate user input values.' // new_line('')
+  end if
+
+  call execute_command_line("grep '<-- scalar' example/extended-gauss-divergence.F90 | grep -v execute_command", wait=.true.)
+  call execute_command_line("grep '<-- vector' example/extended-gauss-divergence.F90 | grep -v execute_command", wait=.true.)
+
+#ifdef __GFORTRAN__
+  command_line_arguments: &
+  block
+    type(numerical_arguments_t) args
+    args  = get_numerical_arguments()
+#else
+  command_line_arguments: &
+  associate(args => get_numerical_arguments())
+#endif
+    text_flags: &
+    associate(flags => text_flags_t( &
+       div_  = command_line%argument_present( [ character(len=len("--div" )) :: "--div" , "-d" ] ) &
+      ,grad_ = command_line%argument_present( [ character(len=len("--grad")) :: "--grad", "-g" ] ) &
+      ,vf_   = command_line%argument_present( [ character(len=len("--vf"  )) :: "--vf"  , "-f" ] ) &
+    ))
+      print_all: &
+      associate(all_terms => merge(.true., .false., all([flags%div_, flags%grad_, flags%vf_]) .or. .not. any([flags%div_, flags%grad_, flags%vf_])))
+        integrand_factors: &
+        associate( &
+           f => scalar_1D_t(scalar_1D_initializer, args%order_, args%cells_, args%x_min_, args%x_max_) &
+          ,v => vector_1D_t(vector_1D_initializer, args%order_, args%cells_, args%x_min_, args%x_max_) &
+        )
+          differential_volume: &
+          associate(dV => f%dV())
+
+            if (flags%grad_ .or. all_terms) then
+              SSS_v_dot_grad_f_dV = .SSS. (v .dot. .grad. f) * dV
+              print '(a,g0)', ".SSS. (v .dot. .grad. f) * dV = ", SSS_v_dot_grad_f_dV
+            end if
+
+            if (flags%div_ .or. all_terms) then
+              SSS_f_div_v_dV      = .SSS. (f * .div. v) * dV
+              print '(a,g0)', ".SSS. (     f * .div. v) * dV = ", SSS_f_div_v_dV
+            end if
+
+          end associate differential_volume
+
+          differential_area: &
+          associate(dA => v%dA())
+            if (flags%vf_ .or. all_terms) then
+              SS_f_v_dot_dA       =  .SS. (f .x. (v .dot. dA))
+              print '(a,g0)', "     -.SS. (f .x. (v .dot. dA)) = ", -SS_f_v_dot_dA
+            end if
+
+            if (all_terms) then
+              print '(a)'   , "----------------------------------------------------"
+              print '(26x,a,g0,a)',"sum = ", SSS_v_dot_grad_f_dV  +  SSS_f_div_v_dV - SS_f_v_dot_dA, " (residual)"
+            end if
+
+          end associate differential_area
+        end associate integrand_factors
+      end associate print_all
+    end associate text_flags
+#ifndef __GFORTRAN__
+  end associate command_line_arguments
+#else
+  end block command_line_arguments
+#endif
+
+contains
+
+  function get_numerical_arguments() result(numerical_arguments)
+    type(numerical_arguments_t) numerical_arguments
+
+#ifdef __GFORTRAN__
+    character(len=:), allocatable :: cells_string, order_string, x_min_string, x_max_string
+    cells_string = command_line%flag_value("--cells")
+    order_string = command_line%flag_value("--order")
+    x_min_string = command_line%flag_value("--x_min")
+    x_max_string = command_line%flag_value("--x_max")
+#else
+    associate( &
+       cells_string => command_line%flag_value("--cells") &
+      ,order_string => command_line%flag_value("--order") &
+      ,x_min_string => command_line%flag_value("--x_min") &
+      ,x_max_string => command_line%flag_value("--x_max") &
+    )
+#endif
+      if (len(cells_string)/=0) read(cells_string,*) numerical_arguments%cells_
+      if (len(order_string)/=0) read(order_string,*) numerical_arguments%order_
+      if (len(x_min_string)/=0) read(x_min_string,*) numerical_arguments%x_min_
+      if (len(x_max_string)/=0) read(x_max_string,*) numerical_arguments%x_max_
+#ifndef __GFORTRAN__
+    end associate
+#endif
+
+  end function
+
+end program

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

@@ -1,6 +1,7 @@
 program print_assembled_1D_operators
-  !! Print fully assembled memetic 1D gradient, divergence, and Laplacian matrices,
-  !! including the zero elements.
+  !! Print the fully assembled 2nd- and 4th-order mimetic gradient, divergence, and Laplacian 
+  !! operator matrices as comma-separated rows for 1D grids with the minimum number of cells (16)
+  !! required for computing 4th-order gradient quadrature weights (Q).
   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
@@ -29,12 +30,15 @@ program print_assembled_1D_operators
     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)
+      if (print_all .or. (gradient   .and. len(order)==0) .or. (gradient   .and. order=="2")) call print_gradient_operator(  k=2, dx=1D0, m=16)
+      if (print_all .or. (divergence .and. len(order)==0) .or. (divergence .and. order=="2")) call print_divergence_operator(k=2, dx=1D0, m=16)
+      if (print_all .or. (gradient   .and. len(order)==0) .or. (gradient   .and. order=="4")) call print_gradient_operator(  k=4, dx=1D0, m=16)
+      if (print_all .or. (divergence .and. len(order)==0) .or. (divergence .and. order=="4")) call print_divergence_operator(k=4, dx=1D0, m=16)
 
     end associate default_usage
+#ifdef __GFORTRAN__
+    stop
+#endif
   end associate command_line_settings
 
 contains

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.3.0"}
+julienne = {git = "https://github.com/berkeleylab/julienne.git", tag = "3.6.0"}
 
 [install]
 library = true

src/fortran/README.md

@@ -15,7 +15,7 @@ and placing the resulting executable program in your `PATH` suffices.
 | Vendor | Compiler    | Version  | Build/Test Command                                    |
 |--------|-------------|----------|-------------------------------------------------------|
 | GCC    | `gfortran`  | 13       | fpm test --compiler gfortran --profile release        |
-| Intel  | `ifx`       | 2025.1.2 | fpm test --compiler ifx --profile release --flag -fpp |
+| Intel  | `ifx`       | 2025.1.2 | FOR_COARRAY_NUM_IMAGES=1 fpm test --compiler ifx --flag "-fpp -O3 -coarray" --profile release |
 | LLVM   | `flang-new` | 19       | fpm test --compiler flang-new --flag "-O3"            |
 | NAG    | `nagfor`    | 7.2      | fpm test --compiler nagfor --flag "-O3 -fpp"          |