2D advection/diffusion equation solver

.gitignore

@@ -1,3 +1,13 @@
+# Imagery
+*.gif
+*.mp4
+*.png
+*.pdf
+
+# Data files
+*.csv
+*.dat
+
 # Scratch directory
 scratch
 

README.md

@@ -48,8 +48,26 @@ where the small residual of approximately $-.222 \times 10^{-15}$ evidences a hi
 **Future work:** Formal lays a foundation for defining DSL embedded in Fortran via template requirements,
 a feature of the forthcoming Fortran 2028 standard.
 
+
 Examples
 --------
+### Highlights
+As demonstrated in this repository's two-dimensional (2D) advection/diffusion partial differential equation
+(PDE) solver, Formal now supports two-dimensional (2D) operators that compute the gradient (`.grad.`) of a
+scalar field, the divergence (`.div.`) of vector field, and arithmetic operators sufficient to express the
+advection/diffusion PDE,
+
+$$ \partial s / \partial t = \nabla \cdot (D \nabla s) - \vec{v} \cdot \nabla s$$
+
+in code as
+```
+ds_dt = .div. (D * .grad. s) - (v .dot. .grad. s)
+```
+where `s` is the concentration of a passive scalar quantity, `D` is a diffusion coefficient, and
+`v` is a prescribed velocity field. The [2D-advection-diffusion.F90](./example/2D-advection-diffusion.F90) 
+program to advance the above equation over time using a Runge-Kutta scheme.
+
+### Other example programs
 See this repository's [example](./example) subdirectory for demonstrations of how
 to use Formal.  For usage information for each example, execute something like
 ```bash
@@ -71,11 +89,11 @@ compiling [fpm-0.12.0.F90] and placing the resulting executable file in your
 
 Building and testing
 --------------------
- Vendor   | Compiler    | Version(s) Tested     | Build/Test Command
-----------|-------------|-----------------------|-------------------
- LFortran | `lfortran`  | latest                | `fpm test --compiler lfortran --flag "--cpp --realloc-lhs-arrays"`
- LLVM     | `flang`     | 20-22                 | `fpm test --compiler flang --profile release`
- LLVM     | `flang`     | 19                    | `fpm test --compiler flang --profile release --flag "-mmlir -allow-assumed-rank"`
+ Vendor   | Compiler    | Version(s) Tested | Build/Test Command
+----------|-------------|-------------------|-------------------
+ LFortran | `lfortran`  | latest            | `fpm test --compiler lfortran --flag "--cpp --realloc-lhs-arrays"`
+ LLVM     | `flang`     | 20-22             | `fpm test --compiler flang --profile release`
+ LLVM     | `flang`     | 19                | `fpm test --compiler flang --profile release --flag "-mmlir -allow-assumed-rank"`
 
 ### `fpm` versions before 0.13.0
 With LLVM, replace the `flang` with `flang-new` and delete `--profile release` from the above commands. 

example/2D-advection-diffusion.F90

@@ -0,0 +1,79 @@
+! Copyright (c) 2026, The Regents of the University of California
+! Terms of use are as specified in LICENSE.txt
+
+module fields_m
+  implicit none
+
+  integer, parameter :: space_dimension = 2
+
+contains
+
+  pure function scalar_field(x,y) result(gaussian)
+    double precision, intent(in) :: x(:), y(:)
+    double precision gaussian(size(x),size(y))
+    double precision, parameter :: pi = acos(-1D0)
+    double precision, parameter :: x0 = 1D0, y0 = 1D0, sigma = pi/5
+    do concurrent(integer :: j=1:size(y)) default(none) shared(x,y,gaussian)
+      associate(r => sqrt((x-x0)**2 + (y(j)-y0)**2))
+        gaussian(:,j) = exp(-(r**2)/(2*sigma**2))
+      end associate
+    end do
+  end function
+
+  pure function stagnation_point_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,:) = [x(i), -y(j)]
+    end do
+  end function
+
+end module
+
+program advection_diffusion_2D
+  !! Solve the advection-diffusion equation for a passive scalar moving through a static
+  !! 2D velocity field.
+  use julienne_m, only : file_t
+  use fields_m, only : scalar_field, stagnation_point_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
+  procedure(scalar_2D_initializer_i), pointer :: scalar_2D_initializer
+  procedure(vector_2D_initializer_i), pointer :: velocity_2D_initializer
+  type(scalar_2D_t) s
+
+  scalar_2D_initializer => scalar_field
+  velocity_2D_initializer => stagnation_point_velocity
+
+  s = scalar_2D_t(scalar_2D_initializer, order=4, cells=[20,20], x_min=[-2D0,-2D0], x_max=[2D0,2D0])
+
+  associate(v => vector_2D_t(velocity_2D_initializer, mold=s))
+
+    advance_time: &
+    block
+      double precision :: dt = 1D-6
+      s = s + dt * d_dt(s, v)
+    end block advance_time
+
+    associate( &
+       scalar_file   => s%to_file("scalar") &
+      ,velocity_file => v%to_file("vector") &
+    )
+      call   scalar_file%write_lines("example/scripts/scalar.csv")
+      call velocity_file%write_lines("example/scripts/velocity.csv")
+    end associate
+
+  end associate
+
+contains
+
+  pure function d_dt(s, v) result(ds_dt)
+    type(scalar_2D_t), intent(in) :: s
+    type(vector_2D_t), intent(in) :: v
+    type(scalar_2D_t) ds_dt
+    double precision, parameter :: D = 1D0
+    ds_dt = .div. (D * .grad. s) - (v .dot. .grad. s)
+  end function
+
+end program

example/2D-sink.F90

@@ -48,7 +48,7 @@ program sink_2D
 
   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())
+      associate(v_file => v%to_file("v"),div_v_file => div_v%to_file(".div. v"), expected_divergence_file => expected_divergence%to_file("expected .div. v"))
         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")

example/2D-stagnation-point.F90

@@ -0,0 +1,69 @@
+! 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) = (x**2 - y(j)**2)/2
+    end do
+  end function
+
+  pure function potential_gradient(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,:) = [x(i), -y(j)]
+    end do
+  end function
+
+end module
+
+program stagnation_point_2D
+  use julienne_m, only : file_t
+  use velocity_potential_m, only : potential, potential_gradient
+  use formal_m, only : scalar_2D_t, vector_2D_t, scalar_2D_initializer_i, vector_2D_initializer_i
+  implicit none
+
+  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 => potential_gradient
+
+  associate( &
+    phi => scalar_2D_t(scalar_2D_initializer, order=4, cells=[20,20], x_min=[-2D0,-2D0], x_max=[2D0,2D0]) &
+  )
+    associate( &
+       v => .grad. phi &
+      ,v_exp => vector_2D_t(vector_2D_initializer, mold=phi) &
+    )
+      associate( &
+         velocity_potential_file => phi%to_file("phi") &
+        ,velocity_file           => v%to_file("velocity") &
+        ,velocity_expected_file  => v_exp%to_file("expected velocity") &
+      )
+        block 
+          character(len=*), parameter :: path = "example/scripts/"
+          call velocity_potential_file%write_lines(path // "velocity-potential.csv")
+          call velocity_file%write_lines(path // "velocity.csv")
+          call velocity_expected_file%write_lines(path // "expected-velocity.csv")
+          print *
+          print '(a)', "With gnuplot installed, plot the results by setting your present working" 
+          print '(a)', "directory to formal/example/scripts and executing the following commands:" // new_line('')
+          print '(a)', 'gnuplot -e "base_name=' // "'velocity'"           // '" 2D-vector-field.gnuplot'
+          print '(a)', 'gnuplot -e "base_name=' // "'expected-velocity'"  // '" 2D-vector-field.gnuplot'
+          print '(a)', 'gnuplot -e "base_name=' // "'velocity-potential'" // '" 2D-scalar-field.gnuplot' // new_line('')
+        end block
+      end associate
+    end associate
+  end associate
+
+end program

example/burgers-1D.F90

@@ -28,7 +28,7 @@ program burgers_1D
   !! * 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, csv, call_julienne_assert_, operator(.equalsExpected.), string_t
+  use julienne_m, only :  command_line_t, call_julienne_assert_, operator(.equalsExpected.), string_t
   use formal_m, only : scalar_1D_t, scalar_1D_initializer_i, d_dx, d2_dx2
   use iso_fortran_env, only : output_unit
   implicit none
@@ -67,7 +67,7 @@ program burgers_1D
   block
     procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer
     double precision, parameter :: pi = acos(-1D0), nu=1D0, t_final=0.6D0
-    double precision, allocatable :: u_surface(:,:), time(:)
+    double precision, allocatable :: u_surface(:,:)
     double precision dt
     type(scalar_1D_t) u
     integer step, n
@@ -122,7 +122,7 @@ program burgers_1D
 
         associate(x => u%grid())
           do n = 1, size(x)
-            write(file_unit,"(*(G13.6,:,'  '))"), x(n), u_surface(n,:) ! write space-separated values
+            write(file_unit,"(*(G13.6,:,'  '))") x(n), u_surface(n,:) ! write space-separated values
           end do
         end associate
 

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

@@ -99,6 +99,7 @@ subroutine output(order)
         associate( s_grid           => s%grid()      &
                   ,grad_s_grid      => grad_s%grid() &
                   ,laplacian_s_grid => laplacian_s%grid())
+#ifndef __GFORTRAN__
           associate( s_table           => tabulate( &
                         string_t([character(len=22)::"x", "f(x) expected"         , "f(x) actual"         ]) &
                        ,s_grid, f(s_grid), s%values() &
@@ -115,6 +116,20 @@ subroutine output(order)
              call grad_s_table%write_lines()
              call laplacian_s_table%write_lines()
           end associate
+#else
+          block
+            type(file_t) s_table, grad_s_table, laplacian_s_table
+            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 block
+#endif
         end associate
       end associate
     end associate

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

@@ -2,7 +2,7 @@
 #  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
+#  Usage:  gnuplot -e "base_name='velocity-potential'" 2D-scalar-field.gnuplot
 #  Default: base_name='velocity-potential'
 # ============================================================================