Fix scalar assignment from a 2D divergence

src/formal/scalar_2D_s.F90

@@ -15,6 +15,7 @@
     ,operator(.isAtLeast.) &
     ,string_t
   use tensors_1D_m, only : cell_centers_extended_1D, scalar_1D_t
+  use differential_operators_1D_m, only : gradient_operator_1D_t
   use interpolator_1D_m, only : centers_to_faces_1D_t
   implicit none
 
@@ -195,10 +196,10 @@
      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 &
+        x_last => size(rhs%points_(1,1,1,1)%values_,x_dir) + 2 &
+       ,y_last => size(rhs%points_(1,1,1,1)%values_,y_dir) + 2 &
      )
-       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_(2:x_last-1, 2:y_last-1) = rhs%points_(1,1,1,1)%values_ ! internal points
        lhs%points_(1,1,1,1)%values_(1         ,  :        ) = 0E0 ! x_min boundary
        lhs%points_(1,1,1,1)%values_(  x_last  ,  :        ) = 0E0 ! x_max boundary
        lhs%points_(1,1,1,1)%values_( :        , 1         ) = 0E0 ! y_min boundary
@@ -209,6 +210,7 @@
      lhs%x_min_ = rhs%x_min_
      lhs%x_max_ = rhs%x_max_
      lhs%order_ = rhs%order_
+     lhs%gradient_operator_1D_ = gradient_operator_1D_t(k=lhs%order_, dx=(lhs%x_max_ - lhs%x_min_)/lhs%cells_, cells=lhs%cells_)
 
      call_julienne_assert(lhs%consistent())
      call_julienne_assert(lhs%conformable(rhs))

test/divergence_assignment_2D_test_m.F90

@@ -0,0 +1,172 @@
+! Copyright (c) 2026, The Regents of the University of California
+! Terms of use are as specified in LICENSE.txt
+
+module divergence_assignment_2D_test_m
+  use julienne_m, only : &
+     operator(//) &
+    ,operator(.all.) &
+    ,operator(.also.) &
+    ,operator(.equalsExpected.) &
+    ,operator(.approximates.) &
+    ,operator(.within.) &
+    ,passing_test &
+    ,string_t &
+    ,test_description_t &
+    ,test_diagnosis_t &
+    ,test_result_t &
+    ,test_t &
+    ,usher
+  use formal_m, only : scalar_2D_t, vector_2D_t, gradient_2D_t, divergence_2D_t, scalar_2D_initializer_i
+  implicit none
+
+  integer, parameter :: x_cells_1 = 9, y_cells_1 = 10
+  integer, parameter :: x_cells_2 = 11, y_cells_2 = 12
+  real, parameter :: tolerance = 1E-6
+
+  type, extends(test_t) :: divergence_assignment_2D_test_t
+  contains
+    procedure, nopass :: subject
+    procedure, nopass :: results
+  end type
+
+contains
+
+  pure function subject() result(test_subject)
+    character(len=:), allocatable :: test_subject
+    test_subject = 'The scalar_2D_t divergence assignment'
+  end function
+
+  function results() result(test_results)
+    type(divergence_assignment_2D_test_t) test_case
+    type(test_result_t), allocatable :: test_results(:)
+
+    test_results = test_case%run([ &
+       test_description_t('copying a nonsquare divergence field into a scalar, preserving boundaries, repeated assignment, and .grad.', usher(check_divergence_assignment)) &
+    ])
+  end function
+
+  pure function pattern_a(x, y) result(values)
+    real, intent(in) :: x(:), y(:)
+    real values(size(x), size(y))
+    do concurrent(integer :: i = 1:size(x), j = 1:size(y)) default(none) shared(x, y, values)
+      values(i,j) = 1000E0 + 10E0*x(i) + y(j)
+    end do
+  end function
+
+  pure function pattern_b(x, y) result(values)
+    real, intent(in) :: x(:), y(:)
+    real values(size(x), size(y))
+    do concurrent(integer :: i = 1:size(x), j = 1:size(y)) default(none) shared(x, y, values)
+      values(i,j) = -500E0 + x(i) - 3E0*y(j)
+    end do
+  end function
+
+  pure function expected_pattern_a(x, y) result(values)
+    real, intent(in) :: x(:), y(:)
+    real values(size(x), size(y))
+    values = 0E0
+    values(2:size(x)-1, 2:size(y)-1) = pattern_a(x(2:size(x)-1), y(2:size(y)-1))
+  end function
+
+  pure function expected_pattern_b(x, y) result(values)
+    real, intent(in) :: x(:), y(:)
+    real values(size(x), size(y))
+    values = 0E0
+    values(2:size(x)-1, 2:size(y)-1) = pattern_b(x(2:size(x)-1), y(2:size(y)-1))
+  end function
+
+  function check_divergence_assignment() result(test_diagnosis)
+    type(test_diagnosis_t) test_diagnosis
+    test_diagnosis = passing_test()
+    test_diagnosis = test_diagnosis .also. check_case(4, [x_cells_1, y_cells_1], [-2E0,-1E0], [1E0,3E0])
+    test_diagnosis = test_diagnosis .also. check_case(2, [x_cells_2, y_cells_2], [-1E0,-2E0], [2E0,4E0])
+  end function
+
+  function check_case(order, cells, x_min, x_max) result(test_diagnosis)
+    integer, intent(in) :: order, cells(2)
+    real, intent(in) :: x_min(2), x_max(2)
+    type(test_diagnosis_t) test_diagnosis
+    procedure(scalar_2D_initializer_i), pointer :: initializer_a, initializer_b
+    procedure(scalar_2D_initializer_i), pointer :: expected_initializer_a, expected_initializer_b
+    type(divergence_2D_t) rhs_a, rhs_b
+    type(scalar_2D_t) lhs, expected_a, expected_b
+    type(gradient_2D_t) grad_lhs, grad_expected
+    real, allocatable :: lhs_values(:,:), rhs_values(:,:)
+
+    test_diagnosis = passing_test()
+    initializer_a => pattern_a
+    initializer_b => pattern_b
+    expected_initializer_a => expected_pattern_a
+    expected_initializer_b => expected_pattern_b
+
+    rhs_a = divergence_2D_t(initializer_a, order=order, cells=cells, x_min=x_min, x_max=x_max)
+    rhs_b = divergence_2D_t(initializer_b, order=order, cells=cells, x_min=x_min, x_max=x_max)
+    expected_a = scalar_2D_t(expected_initializer_a, order=order, cells=cells, x_min=x_min, x_max=x_max)
+    expected_b = scalar_2D_t(expected_initializer_b, order=order, cells=cells, x_min=x_min, x_max=x_max)
+
+    lhs = rhs_a
+    lhs_values = lhs%values()
+    rhs_values = rhs_a%values()
+    test_diagnosis = test_diagnosis .also. (check_shape(lhs_values)) .also. (check_contents(lhs_values, rhs_values)) .also. (check_gradient(lhs, expected_a))
+
+    lhs = rhs_b
+    lhs_values = lhs%values()
+    rhs_values = rhs_b%values()
+    test_diagnosis = test_diagnosis .also. (check_shape(lhs_values)) .also. (check_contents(lhs_values, rhs_values)) .also. (check_gradient(lhs, expected_b))
+
+    lhs = rhs_a
+    lhs_values = lhs%values()
+    rhs_values = rhs_a%values()
+    test_diagnosis = test_diagnosis .also. (check_shape(lhs_values)) .also. (check_contents(lhs_values, rhs_values)) .also. (check_gradient(lhs, expected_a))
+
+  contains
+
+    function check_shape(values) result(diagnosis)
+      real, intent(in) :: values(:,:)
+      type(test_diagnosis_t) diagnosis
+      diagnosis = passing_test()
+      diagnosis = diagnosis .also. (.all. (shape(values) .equalsExpected. [cells(1) + 2, cells(2) + 2]))
+    end function
+
+    function check_contents(values, rhs_values) result(diagnosis)
+      real, intent(in) :: values(:,:), rhs_values(:,:)
+      type(test_diagnosis_t) diagnosis
+      diagnosis = passing_test()
+      diagnosis = diagnosis .also. (boundary_zero(values)) .also. (interior_copy(values, rhs_values))
+    end function
+
+    function boundary_zero(values) result(diagnosis)
+      real, intent(in) :: values(:,:)
+      type(test_diagnosis_t) diagnosis
+      integer :: nx, ny
+      nx = size(values,1)
+      ny = size(values,2)
+      diagnosis = passing_test()
+      diagnosis = diagnosis .also. (maxval(abs(values(1,:))) .approximates. 0E0 .within. tolerance)
+      diagnosis = diagnosis .also. (maxval(abs(values(nx,:))) .approximates. 0E0 .within. tolerance)
+      diagnosis = diagnosis .also. (maxval(abs(values(:,1))) .approximates. 0E0 .within. tolerance)
+      diagnosis = diagnosis .also. (maxval(abs(values(:,ny))) .approximates. 0E0 .within. tolerance)
+    end function
+
+    function interior_copy(values, rhs_values) result(diagnosis)
+      real, intent(in) :: values(:,:), rhs_values(:,:)
+      type(test_diagnosis_t) diagnosis
+      diagnosis = passing_test()
+      diagnosis = diagnosis .also. (.all. (values(2:size(values,1)-1, 2:size(values,2)-1) .approximates. rhs_values .within. tolerance))
+    end function
+
+    function check_gradient(scalar_field, expected_scalar) result(diagnosis)
+      type(scalar_2D_t), intent(in) :: scalar_field
+      type(scalar_2D_t), intent(in) :: expected_scalar
+      type(test_diagnosis_t) diagnosis
+      type(gradient_2D_t) grad_expected
+      diagnosis = passing_test()
+      grad_lhs = .grad. scalar_field
+      grad_expected = .grad. expected_scalar
+      diagnosis = diagnosis .also. (.all. (shape(grad_lhs%to_centers_extended()) .equalsExpected. shape(grad_expected%to_centers_extended())))
+      diagnosis = diagnosis .also. (.all. (grad_lhs%to_centers_extended() .approximates. grad_expected%to_centers_extended() .within. tolerance))
+    end function
+
+  end function check_case
+
+end module divergence_assignment_2D_test_m

test/driver.f90

@@ -11,6 +11,7 @@ program test_suite_driver
   use interpolator_1D_test_m, only : interpolator_1D_test_t
   use scalar_1D_test_m, only : scalar_1D_test_t
   use scalar_2D_test_m, only : scalar_2D_test_t
+  use divergence_assignment_2D_test_m, only : divergence_assignment_2D_test_t
   use scalar_3D_test_m, only : scalar_3D_test_t
   use vector_2D_test_m, only : vector_2D_test_t
   use vector_3D_test_m, only : vector_3D_test_t
@@ -24,6 +25,7 @@ program test_suite_driver
     ,test_fixture_t(interpolator_1D_test_t()) &
     ,test_fixture_t(scalar_1D_test_t()) &
     ,test_fixture_t(scalar_2D_test_t()) &
+    ,test_fixture_t(divergence_assignment_2D_test_t()) &
     ,test_fixture_t(scalar_3D_test_t()) &
     ,test_fixture_t(vector_2D_test_t()) &
     ,test_fixture_t(vector_3D_test_t()) &