2D and 3D gradients and potential vortex example

.gitignore

@@ -1,3 +1,5 @@
+scratch
+
 # Prerequisites
 *.d
 

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=[15,20], x_min=[pi/2,-pi], x_max=[2*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/scripts/velocities.gnuplot

@@ -0,0 +1,34 @@
+# =============================================================
+#  velocities.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'" velocities.gnuplot
+# =============================================================
+
+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 ""

example/scripts/velocity-potential.gnuplot

@@ -0,0 +1,34 @@
+# ============================================================
+#  velocity-potential.gnuplot  --  surface plot CSV
+#  Line 1: column labels
+#  Lines 2+: x, y, z data with blank lines between x-slices
+#  Usage:  gnuplot velocity-potential.gnuplot
+# ============================================================
+
+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

src/formal/differential_operators_1D_m.F90

@@ -51,7 +51,7 @@ pure module function construct_matrix_operator(upper, inner, lower) result(diffe
 
   interface gradient_operator_1D_t
 
-    pure module function construct_1D_gradient_operator(k, dx, cells) result(gradient_operator_1D)
+    elemental module function construct_1D_gradient_operator(k, dx, cells) result(gradient_operator_1D)
       !! Construct a mimetic gradient operator
       implicit none
       integer, intent(in) :: k !! order of accuracy
@@ -77,7 +77,7 @@ pure module function construct_1D_gradient_operator(k, dx, cells) result(gradien
 
   interface divergence_operator_1D_t
 
-    pure module function construct_1D_divergence_operator(k, dx, cells) result(divergence_operator_1D)
+    elemental module function construct_1D_divergence_operator(k, dx, cells) result(divergence_operator_1D)
       !! Construct a mimetic gradient operator
       implicit none
       integer, intent(in) :: k !! order of accuracy

src/formal/scalar_1D_s.F90

@@ -28,7 +28,7 @@
     call_julienne_assert(x_max .greaterThan. x_min)
     call_julienne_assert(cells .isAtLeast. 2*order)
 
-    associate(values => initializer(scalar_1D_grid_locations(x_min, x_max, cells)))
+    associate(values => initializer(cell_centers_extended_1D(x_min, x_max, cells)))
       scalar_1D%tensor_1D_t = tensor_1D_t(values, x_min, x_max, cells, order)
     end associate
     scalar_1D%gradient_operator_1D_ = gradient_operator_1D_t(k=order, dx=(x_max - x_min)/cells, cells=cells)
@@ -47,7 +47,7 @@ pure module function construct_1D_scalar_from_function(initializer, order, cells
     call_julienne_assert(x_max .greaterThan. x_min)
     call_julienne_assert(cells .isAtLeast. 2*order)
 
-    associate(values => initializer(scalar_1D_grid_locations(x_min, x_max, cells)))
+    associate(values => initializer(cell_centers_extended_1D(x_min, x_max, cells)))
       scalar_1D%tensor_1D_t = tensor_1D_t(values, x_min, x_max, cells, order)
     end associate
     scalar_1D%gradient_operator_1D_ = gradient_operator_1D_t(k=order, dx=(x_max - x_min)/cells, cells=cells)
@@ -66,6 +66,22 @@ pure function cell_center_locations(x_min, x_max, cells) result(x)
 
 #endif
 
+  module procedure construct_1D_scalar_constant
+
+    integer i
+
+    call_julienne_assert(x_max .greaterThan. x_min)
+    call_julienne_assert(cells .isAtLeast. 2*order)
+
+    scalar_1D = scalar_1D_t( tensor_1D_t( &
+         values = [(constant, i = 1, size(cell_centers_extended_1D(x_min, x_max, cells)))] &
+        ,x_min = x_min &
+        ,x_max = x_max &
+        ,cells = cells &
+        ,order = order &
+    )   )
+  end procedure
+
   module procedure divide_by_integer
     ratio%tensor_1D_t = tensor_1D_t( &
       values = self%values_/denominator, x_min = self%x_min_, x_max = self%x_max_, cells = self%cells_, order = self%order_ &
@@ -143,14 +159,12 @@ pure function cell_center_locations(x_min, x_max, cells) result(x)
     integer c
 
     associate(dx => (self%x_max_ - self%x_min_)/self%cells_)
-      associate(G => gradient_operator_1D_t(self%order_, dx, self%cells_))
-        gradient_1D%tensor_1D_t = tensor_1D_t(G .x. self%values_, self%x_min_, self%x_max_, cells=self%cells_, order=self%order_)
-        gradient_1D%divergence_operator_1D_ = divergence_operator_1D_t(self%order_, dx, self%cells_)
-        check_corbino_castillo_eq_17: &
-        associate(p => gradient_1D%weights(), b => [-1D0, [(0D0, c = 1, self%cells_)], 1D0])
-          call_julienne_assert((.all. (matmul(transpose(G%assemble()), p) .approximates. b/dx .within. 2D-3)))
-        end associate check_corbino_castillo_eq_17
-      end associate
+      gradient_1D%tensor_1D_t = tensor_1D_t(self%gradient_operator_1D_ .x. self%values_, self%x_min_, self%x_max_, cells=self%cells_, order=self%order_)
+      gradient_1D%divergence_operator_1D_ = divergence_operator_1D_t(self%order_, dx, self%cells_)
+      check_corbino_castillo_eq_17: &
+      associate(p => gradient_1D%weights(), b => [-1D0, [(0D0, c = 1, self%cells_)], 1D0])
+        call_julienne_assert((.all. (matmul(transpose(self%gradient_operator_1D_%assemble()), p) .approximates. b/dx .within. 2D-3)))
+      end associate check_corbino_castillo_eq_17
     end associate
 
   end procedure
@@ -196,19 +210,8 @@ pure function cell_center_locations(x_min, x_max, cells) result(x)
     cell_centers_extended_values = self%values_
   end procedure
 
-  pure function scalar_1D_grid_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, cell_center_locations(x_min, x_max, cells), x_max]
-    end associate
-  end function
-
   module procedure scalar_1D_grid
-    cell_centers_extended  = scalar_1D_grid_locations(self%x_min_, self%x_max_, self%cells_)
+    cell_centers_extended  = cell_centers_extended_1D(self%x_min_, self%x_max_, self%cells_)
   end procedure
 
 end submodule scalar_1D_s

src/formal/scalar_2D_s.F90

@@ -0,0 +1,110 @@
+! 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_2D_m) scalar_2D_s
+  use julienne_m, only : &
+     call_julienne_assert_ &
+    ,operator(.all.) &
+    ,operator(.equalsExpected.) &
+    ,operator(.greaterThan.) &
+    ,operator(.isAtLeast.)
+  use tensors_1D_m, only : cell_centers_extended_1D, scalar_1D_t
+  use julienne_m, only : string_t, operator(.csv.)
+  implicit none
+
+contains
+
+  module procedure scalar_2D_values
+    scalar_values = self%values_(:,:,1,1,1,1)
+  end procedure
+
+  module procedure scalar_2D_grid
+    associate(scalar_1D => scalar_1D_t( &
+       constant = 0D0 &
+      ,cells = self%cells_(direction) &
+      ,x_min = self%x_min_(direction) &
+      ,x_max = self%x_max_(direction) &
+      ,order = self%order_ &
+    ))
+      scalar_grid_1D = scalar_1D%grid()
+    end associate
+  end procedure
+
+  module procedure construct_2D_scalar_from_function
+
+    call_julienne_assert(.all. ([size(cells), size(x_min), size(x_max)] .equalsExpected. space_dimension))
+    call_julienne_assert(.all. (x_max .greaterThan. x_min))
+    call_julienne_assert(.all. (cells .isAtLeast. 2*order))
+
+    associate(x => cell_centers_extended_1D(x_min(1), x_max(1), cells(1)), y => cell_centers_extended_1D(x_min(2), x_max(2), cells(2)))
+      scalar_2D%tensor_2D_t = tensor_2D_t( &
+         values = reshape(initializer(x,y), shape=[size(x),size(y),1,1,1,1]) &
+        ,cells = cells , x_min = x_min, x_max = x_max, order = order &
+      )
+      scalar_2D%gradient_operator_1D_ = gradient_operator_1D_t(k=order, dx=(x_max - x_min)/cells, cells=cells)
+    end associate
+  end procedure
+
+  module procedure construct_2D_scalar_from_mold
+    scalar_2D = scalar_2D_t(initializer, cells = mold%cells_, x_min = mold%x_min_, x_max = mold%x_max_, order = mold%order_)
+  end procedure
+
+  module procedure scalar_2D_gradient
+
+    integer c, i, j
+
+    gradient_2D%x_min_ = self%x_min_
+    gradient_2D%x_max_ = self%x_max_
+    gradient_2D%cells_ = self%cells_
+    gradient_2D%order_ = self%order_
+
+    allocate(gradient_2D%values_(self%cells_(1)+1, self%cells_(2)+1, space_dimension, 1, 1, 1))
+
+    gradient_x_component: &
+    do concurrent(integer :: j=1:size(gradient_2D%values_,2)) default(none) shared(gradient_2D, self)
+      gradient_2D%values_(:,j,1,1,1,1) = self%gradient_operator_1D_(1) .x. self%values_(:,j,1,1,1,1)
+    end do gradient_x_component
+
+    gradient_y_component: &
+    do concurrent(integer :: i=1:size(gradient_2D%values_,1)) default(none) shared(gradient_2D, self)
+      gradient_2D%values_(i,:,2,1,1,1) = self%gradient_operator_1D_(2) .x. self%values_(i,:,1,1,1,1)
+    end do gradient_y_component
+
+    associate(dx => (self%x_max_ - self%x_min_)/self%cells_)
+      gradient_2D%divergence_operator_1D_ = divergence_operator_1D_t(self%order_, dx, self%cells_)
+     !check_corbino_castillo_eq_17: &
+     !associate(p => gradient_1D%weights(), b => [-1D0, [(0D0, c = 1, self%cells_)], 1D0])
+     !  call_julienne_assert((.all. (matmul(transpose(self%gradient_operator_1D_%assemble()), p) .approximates. b/dx .within. 2D-3)))
+     !end associate check_corbino_castillo_eq_17
+    end associate
+
+  end procedure
+
+  module procedure scalar_2D_to_file
+    type(string_t), allocatable :: lines(:)
+    integer i, j, l
+
+    associate(x => self%grid(1), y => self%grid(2), header => [string_t("x,y,scalar")])
+      associate(num_blank_lines => size(y)-1)
+        allocate(lines(size(header) + size(self%values_) + num_blank_lines))
+      end associate
+      lines(1:size(header)) = header
+      l = size(header)
+      do j = 1, size(y)
+        do i = 1, size(x)
+          l = l + 1
+          lines(l) = .csv. string_t([x(i), y(j), self%values_(i,j,1,1,1,1)])
+        end do
+        if (j/=size(y)) then
+          l = l + 1
+          lines(l) = ""
+        end if
+      end do
+    end associate
+
+    file = file_t(lines)
+  end procedure
+
+end submodule scalar_2D_s