Features: 2D/3D types, operators, output, examples, & plotting scripts

.github/workflows/build.yml

@@ -224,7 +224,7 @@ jobs:
           fi
         elif test "$FC" = "lfortran" ; then
           echo "FPM_FC=lfortran" >> "$GITHUB_ENV"
-          echo "FFLAGS=--cpp --separate-compilation --realloc-lhs-arrays $FFLAGS" >> "$GITHUB_ENV"
+          echo "FFLAGS=--cpp --realloc-lhs-arrays $FFLAGS" >> "$GITHUB_ENV"
           echo "FPM_FLAGS=--profile debug --verbose" >> "$GITHUB_ENV" ; : fpm 0.13 workaround
         else
           echo "FPM_FC=gfortran-${COMPILER_VERSION}" >> "$GITHUB_ENV"

.gitignore

@@ -1,3 +1,9 @@
+# Scratch directory
+scratch
+
+# macOS folder display format
+.DS_Store
+
 # Prerequisites
 *.d
 

README.md

@@ -78,14 +78,7 @@ Building and testing
  LLVM     | `flang`     | 19                    | `fpm test --compiler flang --profile release --flag "-mmlir -allow-assumed-rank"`
 
 ### `fpm` versions before 0.13.0
-With LLVM 20-22, replace the above `flang` command with
-```
-  fpm test --compiler flang --flag "-O3"
-```
-With LLVM 19, replace the above `flang` command with
-```
- fpm test --compiler flang --flag "-O3 -mmlir -allow-assumed-rank"
-```
+With LLVM, replace the `flang` with `flang-new` and delete `--profile release` from the above commands. 
 
 Unsupported Compilers
 ---------------------

example/2D-sink.F90

@@ -0,0 +1,59 @@
+! Copyright (c) 2026, The Regents of the University of California
+! Terms of use are as specified in LICENSE.txt
+
+#include "julienne-assert-macros.h"
+
+module sink_2D_functions_m
+  use julienne_m, only : call_julienne_assert_
+  implicit none
+
+  integer, parameter :: space_dimension = 2
+  double precision, parameter :: pi = acos(-1D0)
+
+contains
+
+  pure function velocity(x,y) result(v)
+    double precision, intent(in) :: x(:), y(:)
+    double precision v(size(x),size(y),space_dimension), theta
+    double precision, parameter :: Q = 1D0
+    do concurrent(integer :: i=1:size(x), j=1:size(y))
+      associate(r => sqrt(x(i)**2 + y(j)**2))
+        call_julienne_assert(r /= 0D0)
+        v(i,j,:) = -(Q/(2*pi))*[x(i), y(j)]/(x(i)**2 + y(j)**2)
+      end associate
+    end do
+  end function
+
+  pure function divergence(x,y) result(div_v)
+    double precision, intent(in) :: x(:), y(:)
+    double precision div_v(size(x),size(y))
+    call_julienne_assert(.not. any(x == 0D0 .and. y == 0D0))
+    div_v = 0D0
+  end function
+
+end module
+
+program sink_2D
+  use julienne_m, only : file_t
+  use sink_2D_functions_m, only : velocity, divergence, pi
+  use formal_m, only : vector_2D_t, divergence_2D_t, divergence_2D_initializer_i, vector_2D_initializer_i
+  implicit none
+
+  integer, parameter ::  order = 4
+  procedure(vector_2D_initializer_i), pointer :: vector_2D_initializer
+  procedure(divergence_2D_initializer_i), pointer :: divergence_2D_initializer
+
+  divergence_2D_initializer => divergence
+  vector_2D_initializer => velocity
+
+  associate(v => vector_2D_t(vector_2D_initializer, order=order, cells=[11,11], x_min=[-1D0,-1D0], x_max=[1D0,1D0]))
+    associate(div_v => .div. v, expected_divergence => divergence_2D_t(divergence_2D_initializer, mold=v))
+      associate(v_file => v%to_file(),div_v_file => div_v%to_file(), expected_divergence_file => expected_divergence%to_file())
+        call v_file%write_lines("example/scripts/sink-velocity.csv") 
+        call div_v_file%write_lines("example/scripts/sink-divergence.csv") 
+        call expected_divergence_file%write_lines("example/scripts/expected-divergence.csv")
+      end associate
+    end associate
+  end associate
+
+end program

example/2D-vortex.F90

@@ -0,0 +1,58 @@
+! Copyright (c) 2026, The Regents of the University of California
+! Terms of use are as specified in LICENSE.txt
+
+module velocity_potential_m 
+  implicit none
+
+  integer, parameter :: space_dimension = 2
+
+contains
+
+  pure function potential(x,y) result(phi)
+    double precision, intent(in) :: x(:), y(:)
+    double precision phi(size(x),size(y))
+    do concurrent(integer :: j=1:size(y)) default(none) shared(x,y,phi)
+      phi(:,j) = atan(y(j)/x)
+    end do
+  end function
+
+  pure function velocity(x,y) result(grad_phi)
+    double precision, intent(in) :: x(:), y(:)
+    double precision grad_phi(size(x),size(y),space_dimension)
+    do concurrent(integer :: i=1:size(x), j=1:size(y))
+      grad_phi(i,j,:) = [-y(j)/(x(i)**2 + y(j)**2), x(i)/(x(i)**2 + y(j)**2)]
+    end do
+  end function
+
+end module
+
+program vortex_2D
+  use julienne_m, only : file_t
+  use velocity_potential_m, only : potential, velocity
+  use formal_m, only : scalar_2D_t, vector_2D_t, scalar_2D_initializer_i, vector_2D_initializer_i
+  implicit none
+
+  integer, parameter ::  order = 4
+  double precision, parameter :: pi = acos(-1D0)
+  procedure(scalar_2D_initializer_i), pointer :: scalar_2D_initializer
+  procedure(vector_2D_initializer_i), pointer :: vector_2D_initializer
+
+  scalar_2D_initializer => potential
+  vector_2D_initializer => velocity
+
+  associate(phi => scalar_2D_t(scalar_2D_initializer, order=order, cells=[11,11], x_min=[-pi,-pi], x_max=[pi,pi]))
+    associate(  velocity => .grad. phi & 
+      ,expected_velocity => vector_2D_t(vector_2D_initializer, mold=phi) &
+    ) 
+      associate(velocity_potential_file =>               phi%to_file() &
+                         ,velocity_file =>          velocity%to_file() &
+                ,expected_velocity_file => expected_velocity%to_file() &
+    )
+        call velocity_potential_file%write_lines("example/scripts/velocity-potential.csv") 
+        call velocity_file%write_lines("example/scripts/velocity.csv") 
+        call expected_velocity_file%write_lines("example/scripts/expected-velocity.csv")
+      end associate
+    end associate
+  end associate
+
+end program

example/burgers-1D.F90

@@ -33,9 +33,6 @@ program burgers_1D
   use iso_fortran_env, only : output_unit
   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
@@ -68,9 +65,7 @@ program burgers_1D
 
 
   block
-#ifndef __GFORTRAN__
     procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer
-#endif
     double precision, parameter :: pi = acos(-1D0), nu=1D0, t_final=0.6D0
     double precision, allocatable :: u_surface(:,:), time(:)
     double precision dt
@@ -149,10 +144,6 @@ program burgers_1D
     end associate
   end block
 
-#ifdef __GFORTRAN__
-    stop
-#endif
-
 contains
 
   pure function d_dt(u, nu) result(du_dt)

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

@@ -42,9 +42,6 @@ program div_grad_laplacian_1D
   use functions_m, only : f, df_dx, d2f_dx2
   use julienne_m, only :  file_t, string_t, operator(.separatedBy.), command_line_t
   use formal_m, only : scalar_1D_t, scalar_1D_initializer_i
-#ifdef __GFORTRAN__
-  use formal_m, only : vector_1D_t, laplacian_1D_t, gradient_1D_t
-#endif
   implicit none
 
   procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => f
@@ -93,8 +90,6 @@ program div_grad_laplacian_1D
 
 contains
 
-#ifndef __GFORTRAN__
-
   subroutine output(order)
     integer, intent(in) :: order
 
@@ -125,44 +120,6 @@ subroutine output(order)
     end associate
   end subroutine
 
-#else
-
-  subroutine output(order)
-    integer, intent(in) :: order
-
-    type(scalar_1D_t) 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=20, 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=22)::"x", "f(x) expected"         , "f(x) actual"         ]) &
-      ,s_grid, f(s_grid), s%values() &
-    )
-    grad_s_table      = tabulate( &
-       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=22)::"x", ".laplacian. f expected", ".laplacian. f actual"]) &
-      ,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(:)

example/extended-gauss-divergence.F90

@@ -61,15 +61,8 @@ program extended_gauss_divergence
   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" ] ) &
@@ -114,39 +107,24 @@ program extended_gauss_divergence
         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

@@ -40,9 +40,6 @@ program print_assembled_1D_operators
       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

example/scripts/2D-scalar-field.gnuplot

@@ -0,0 +1,37 @@
+# ============================================================================
+#  2D-scalar-field.gnuplot  --  surface plot CSV
+#  Line 1: column labels
+#  Lines 2+: x, y, z data with blank lines between x-slices
+#  Usage:  gnuplot -d "base_name='velocity-potential'" 2D-scalar-field.gnuplot
+#  Default: base_name='velocity-potential'
+# ============================================================================
+
+if (!exists("base_name")) base_name = "velocity-potential"
+
+datafile = base_name . ".csv"
+
+set datafile separator ","
+
+# --- 1. Read column headers from line 1 ---
+xlabel = "" ; ylabel = "" ; zlabel = ""
+set table $Dummy
+    plot datafile every ::0::0 \
+        using (xlabel=strcol(1), ylabel=strcol(2), zlabel=strcol(3), 0):(0) \
+        with table
+unset table
+
+# --- 2. Plot ---
+set title  zlabel . "(" . xlabel . ", " . ylabel . ")"
+set xlabel xlabel ; set ylabel ylabel
+set zlabel zlabel offset 3,0 ; set cblabel zlabel
+set hidden3d
+set pm3d depthorder
+set palette rgbformulae 33,13,10
+set ticslevel 0 ; set key off
+
+set terminal gif size 800,600
+set output base_name . ".gif"
+
+splot datafile every ::1 using 1:2:3 with pm3d title ""
+
+set output    # flush and close the file

example/scripts/2D-vector-field.gnuplot

@@ -0,0 +1,35 @@
+# ===============================================================
+#  vector-field.gnuplot  --  2D vector/quiver plot from a CSV
+#  Line 1: column labels
+#  Lines 2+: x, y, velocity_x, velocity_y data
+#  Usage:  gnuplot -e "base_name='velocity'" vector-field.gnuplot
+#  Default: base_name='velocity'
+# ===============================================================
+
+if (!exists("base_name")) base_name = "velocity"
+
+datafile = base_name . ".csv"
+set datafile separator ","
+
+# --- 1. Read column headers from line 1 ---
+xlabel = "" ; ylabel = "" ; dxlabel = "" ; dylabel = ""
+set table $Dummy
+    plot datafile every ::0::0 \
+        using (xlabel=strcol(1), ylabel=strcol(2), \
+               dxlabel=strcol(3), dylabel=strcol(4), 0):(0) \
+        with table
+unset table
+
+# --- 2. Plot ---
+set title  dxlabel . "," . dylabel . " at each " . xlabel . "," . ylabel
+set xlabel xlabel
+set ylabel ylabel
+set key off
+set cblabel "magnitude"
+
+set terminal gif size 800,600
+set output base_name . ".gif"
+
+plot datafile every ::1 \
+    using ($1-$3/2):($2-$4/2):3:4:(sqrt($3**2+$4**2)) \
+    with vectors head filled size screen 0.02,15 lw 1.5 lc palette z title ""