Conservative advection/diffusion equation

example/2D-advection-diffusion.F90

@@ -53,15 +53,17 @@ program advection_diffusion_2D
     advance_time: &
     block
       double precision :: dt = 1D-6
-      s = s + dt * d_dt(s, v)
+      associate(s_half => s + (dt/2) * d_dt(s, v))
+        s = s + dt * d_dt(s_half, v)
+      end associate
     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")
+      call   scalar_file%write_lines("example/scripts/scalar-adv-dif.csv")
+      call velocity_file%write_lines("example/scripts/velocity-adv-dif.csv")
     end associate
 
   end associate
@@ -73,7 +75,7 @@ pure function d_dt(s, v) result(ds_dt)
     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)
+    ds_dt = .div. (D * .grad. s) - .div. (v * s)
   end function
 
 end program

src/formal/divergence_2D_s.F90

@@ -68,6 +68,20 @@
 
   end procedure
 
+  module procedure divergence_2D_minus_divergence
+    call_julienne_assert(lhs%conformable(rhs))
+
+    difference%tensor_2D_t =  tensor_2D_t( &
+       points = reshape([points_2D_t(lhs%points_(1,1,1,1)%values_ - rhs%points_(1,1,1,1)%values_)], shape = [1,1,1,1]) &
+      ,cells  = lhs%cells_ &
+      ,x_min  = lhs%x_min_ &
+      ,x_max  = lhs%x_max_ &
+      ,order  = lhs%order_ &
+    )
+    call_julienne_assert(difference%consistent())
+
+  end procedure
+
   module procedure construct_2D_divergence_from_vector_mold
 
     call_julienne_assert(mold%consistent())

src/formal/scalar_2D_s.F90

@@ -168,6 +168,34 @@
 
   end procedure
 
+  module procedure scalar_2D_assign_divergence
+
+     call_julienne_assert(lhs%conformable(rhs))
+
+     allocate(lhs%points_(1,1,1,1))
+     allocate(lhs%points_(1,1,1,1)%values_(rhs%cells_(x_dir)+2, rhs%cells_(y_dir)+2))
+
+     associate( &
+        x_last => size(rhs%points_(1,1,1,1)%values_,x_dir) - 1 &
+       ,y_last => size(rhs%points_(1,1,1,1)%values_,y_dir) - 1 &
+     )
+       lhs%points_(1,1,1,1)%values_(2:x_last-1, 2:y_last-1) = rhs%points_(1,1,1,1)%values_(2:x_last-1, 2:y_last-1) ! internal points
+       lhs%points_(1,1,1,1)%values_(1         ,  :        ) = 0D0 ! x_min boundary
+       lhs%points_(1,1,1,1)%values_(  x_last  ,  :        ) = 0D0 ! x_max boundary
+       lhs%points_(1,1,1,1)%values_( :        , 1         ) = 0D0 ! y_min boundary
+       lhs%points_(1,1,1,1)%values_( :        ,   y_last  ) = 0D0 ! y_max boundary
+     end associate
+
+     lhs%cells_ = rhs%cells_
+     lhs%x_min_ = rhs%x_min_
+     lhs%x_max_ = rhs%x_max_
+     lhs%order_ = rhs%order_
+
+     call_julienne_assert(lhs%consistent())
+
+  end procedure
+
+
   module procedure scalar_2D_to_file
     type(string_t), allocatable :: lines(:)
     integer i, j, l, m, n, p, q

src/formal/tensors_2D_m.F90

@@ -9,7 +9,7 @@ module tensors_2D_m
   use julienne_m, only : file_t 
   implicit none
 
-  private
+  !private
   public :: scalar_2D_t
   public :: vector_2D_t
   public :: gradient_2D_t
@@ -85,6 +85,7 @@ pure module function construct_2D_tensor_from_components(points, cells, x_min, x
     private
     type(gradient_operator_1D_t) gradient_operator_1D_(space_dimension)
   contains
+    generic :: assignment(=) => scalar_2D_assign_divergence
     generic :: operator(.grad.) => scalar_2D_gradient
     generic :: operator(*) => scalar_2D_postmultiply_double, scalar_2D_premultiply_double
     generic :: operator(+) => scalar_2D_plus_scalar
@@ -93,6 +94,7 @@ pure module function construct_2D_tensor_from_components(points, cells, x_min, x
     generic :: consistent => scalar_2D_consistent
     generic :: conformable => scalar_2D_conformable_scalar
     generic :: to_file => scalar_2D_to_file
+    procedure, non_overridable, private :: scalar_2D_assign_divergence
     procedure, non_overridable, private :: scalar_2D_to_file
     procedure, non_overridable, private :: scalar_2D_gradient
     procedure, non_overridable, private :: scalar_2D_values
@@ -143,18 +145,21 @@ pure module function construct_2D_scalar_from_components(tensor_2D, gradient_ope
     generic :: grid => vector_2D_grid
     generic :: consistent => vector_2D_consistent
     generic :: conformable => vector_2D_conformable_vector, vector_2D_conformable_scalar
-    generic :: co_located_components => vector_2D_co_located_components
+    generic :: to_centers_extended => vector_2D_to_centers_extended
     generic :: operator(.div.) => vector_2D_divergence
     generic :: operator(.dot.) => vector_2D_dot_vector
+    generic :: operator(*) => vector_2D_postmultiply_scalar, vector_2D_premultiply_scalar
     generic :: to_file => vector_2D_to_file
     procedure, non_overridable, private :: vector_2D_to_file
     procedure, non_overridable, private :: vector_2D_grid
     procedure, non_overridable, private :: vector_2D_divergence
     procedure, non_overridable, private :: vector_2D_consistent
     procedure, non_overridable, private :: vector_2D_conformable_vector
     procedure, non_overridable, private :: vector_2D_conformable_scalar
-    procedure, non_overridable, private :: vector_2D_co_located_components
+    procedure, non_overridable, private :: vector_2D_to_centers_extended
     procedure, non_overridable, private :: vector_2D_dot_vector
+    procedure, non_overridable, private :: vector_2D_postmultiply_scalar
+    procedure, non_overridable, private, pass(vector_2D) :: vector_2D_premultiply_scalar
   end type
 
   interface vector_2D_t
@@ -227,7 +232,7 @@ pure module function construct_2D_gradient_from_components(tensor_2D, divergence
     generic :: consistent => tensor_2D_consistent
     generic :: conformable => divergence_2D_conformable_scalar, divergence_2D_conformable_vector
     generic :: operator(*) => divergence_2D_premultiply_constant, divergence_2D_postmultiply_constant
-    generic :: operator(-) => divergence_2D_minus_scalar
+    generic :: operator(-) => divergence_2D_minus_scalar, divergence_2D_minus_divergence
     generic :: to_file => divergence_2D_to_file
     procedure, non_overridable, private :: divergence_2D_to_file
     procedure, private, non_overridable :: divergence_2D_values
@@ -236,6 +241,7 @@ pure module function construct_2D_gradient_from_components(tensor_2D, divergence
     procedure, private, non_overridable :: divergence_2D_conformable_vector
     procedure, private, non_overridable :: divergence_2D_conformable_scalar
     procedure, private, non_overridable :: divergence_2D_postmultiply_constant
+    procedure, private, non_overridable :: divergence_2D_minus_divergence
     procedure, private, non_overridable, pass(rhs) :: divergence_2D_premultiply_constant
   end type
 
@@ -278,6 +284,13 @@ pure module function tensor_2D_conformable(self, tensor_2D) result(conformable)
       logical conformable
     end function
 
+    pure module subroutine scalar_2D_assign_divergence(lhs, rhs)
+      !! Assign 2D divergence to 2D scalar at internal points
+      implicit none
+      class(scalar_2D_t), intent(inout) :: lhs
+      type(divergence_2D_t), intent(in) :: rhs
+    end subroutine
+
     pure module function scalar_2D_consistent(self) result(self_consistent)
       !! Assert components allocated and self-consistent, including sufficient accuracy for gradient operator
       implicit none
@@ -384,7 +397,7 @@ pure module function divergence_2D_grid(self, direction) result(divergence_grid_
       double precision, allocatable :: divergence_grid_1D(:) !! grid points along the requested coordinate direction
     end function
 
-    pure module function vector_2D_co_located_components(self) result(vectors)
+    pure module function vector_2D_to_centers_extended(self) result(vectors)
       !! Vector values getter
       implicit none
       class(vector_2D_t), intent(in) :: self
@@ -412,6 +425,22 @@ pure module function vector_2D_divergence(self) result(divergence_2D)
       type(divergence_2D_t) divergence_2D
     end function
 
+    pure module function vector_2D_postmultiply_scalar(vector_2D, scalar_2D) result(vector_x_scalar)
+      !! Result is product of the 2D vector and scalar arguments
+      implicit none
+      class(vector_2D_t), intent(in) :: vector_2D
+      type(scalar_2D_t), intent(in) :: scalar_2D
+      type(vector_2D_t) vector_x_scalar
+    end function
+
+    pure module function vector_2D_premultiply_scalar(scalar_2D, vector_2D) result(scalar_x_vector)
+      !! Result is product of the 2D vector and scalar arguments
+      implicit none
+      class(vector_2D_t), intent(in) :: vector_2D
+      type(scalar_2D_t), intent(in) :: scalar_2D
+      type(vector_2D_t) scalar_x_vector
+    end function
+
     pure module function vector_2D_dot_vector(lhs, rhs) result(scalar_2D)
       !! Result is scalar product of the 2D-vector arguments
       implicit none
@@ -451,6 +480,13 @@ pure module function divergence_2D_premultiply_constant(lhs, rhs) result(lhs_x_r
       type(divergence_2D_t) lhs_x_rhs
     end function
 
+    pure module function divergence_2D_minus_divergence(lhs, rhs) result(difference)
+      !! Result is the pointwise difference between the lhs and rhs
+      implicit none
+      class(divergence_2D_t), intent(in) :: lhs, rhs
+      type(divergence_2D_t) difference
+    end function
+
     pure module function divergence_2D_minus_scalar(lhs, rhs) result(difference)
       !! Result is the pointwise difference between the lhs and rhs
       implicit none

src/formal/vector_2D_s.F90

@@ -156,7 +156,7 @@ pure function description(coordinate, component) result(point_cloud)
 
   end procedure
 
-  module procedure vector_2D_co_located_components
+  module procedure vector_2D_to_centers_extended
 
     call_julienne_assert(self%consistent())
 
@@ -179,13 +179,46 @@ pure function description(coordinate, component) result(point_cloud)
 
   end procedure
 
+  module procedure vector_2D_postmultiply_scalar
+
+    call_julienne_assert(vector_2D%conformable(scalar_2D))
+
+    associate(vector => vector_2D%to_centers_extended())
+
+      call_julienne_assert(.all. ([size(vector,x_dir), size(vector,y_dir)] .equalsExpected. shape(scalar_2D%points(1,1,1,1)%values_)))
+
+      vector_x_scalar = construct_2D_vector_from_components( &
+         tensor_2D_t( &
+            points = reshape( &
+               source = [ points_2D_t(vector(:,:,x_dir) * scalar_2D%points_(1,1,1,1)%values_) &
+                         ,points_2D_t(vector(:,:,y_dir) * scalar_2D%points_(1,1,1,1)%values_)] &
+              ,shape  = [space_dimension,1,1,1] &
+            ) &
+           ,cells = vector_2D%cells_ &
+           ,x_min = vector_2D%x_min_ &
+           ,x_max = vector_2D%x_max_ &
+           ,order = vector_2D%order_ &
+         ) &
+        ,divergence_operator_1D_t( &
+           k = vector_2D%order_ &
+          ,dx = (vector_2D%x_max_ - vector_2D%x_min_)/vector_2D%cells_ &
+          ,cells = vector_2D%cells_ &
+      )  )
+    end associate
+
+  end procedure
+
+  module procedure vector_2D_premultiply_scalar
+    scalar_x_vector = vector_2D * scalar_2D
+  end procedure
+
   module procedure vector_2D_dot_vector
 
     call_julienne_assert(lhs%conformable(rhs))
 
     associate( &
-      lhs_ => lhs%co_located_components() &
-     ,rhs_ => rhs%co_located_components() &
+      lhs_ => lhs%to_centers_extended() &
+     ,rhs_ => rhs%to_centers_extended() &
     )
       call_julienne_assert(.all. (shape(lhs_) .equalsExpected. shape(rhs_)))
 
@@ -217,7 +250,7 @@ pure function description(coordinate, component) result(point_cloud)
        header => [string_t("x, y, " // name)] &
       ,x => cell_centers_extended_1D(self%x_min_(x_dir), self%x_max_(x_dir), self%cells_(x_dir)) &
       ,y => cell_centers_extended_1D(self%x_min_(y_dir), self%x_max_(y_dir), self%cells_(y_dir)) &
-      ,vectors => self%co_located_components() &
+      ,vectors => self%to_centers_extended() &
     )
       allocate(lines(size(header) +  size(x)*size(y)))
 

test/scalar_2D_test_m.F90

@@ -77,7 +77,7 @@ function check_gradient() result(test_diagnosis)
          ,expected_gradient => vector_2D_t(expected_gradient_initializer, mold=scalar_2D) &
         )
           test_diagnosis = test_diagnosis .also. &
-            .all. (grad_scalar%co_located_components() .approximates. expected_gradient%co_located_components() .within. tolerance) &
+            .all. (grad_scalar%to_centers_extended() .approximates. expected_gradient%to_centers_extended() .within. tolerance) &
             // string_t(" for order ") // string_t(order)
         end associate
       end associate