From b912010a4dc3c61b186802528c8d70f3f1bf63e9 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Thu, 9 Jul 2026 21:50:44 -0700 Subject: [PATCH 1/4] refact(co_located_components): rename to_centers_extended --- src/formal/tensors_2D_m.F90 | 6 +++--- src/formal/vector_2D_s.F90 | 8 ++++---- test/scalar_2D_test_m.F90 | 2 +- 3 files changed, 8 insertions(+), 8 deletions(-) diff --git a/src/formal/tensors_2D_m.F90 b/src/formal/tensors_2D_m.F90 index ba48bf3..4936b0d 100644 --- a/src/formal/tensors_2D_m.F90 +++ b/src/formal/tensors_2D_m.F90 @@ -143,7 +143,7 @@ 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 :: to_file => vector_2D_to_file @@ -153,7 +153,7 @@ pure module function construct_2D_scalar_from_components(tensor_2D, gradient_ope 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 end type @@ -384,7 +384,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 diff --git a/src/formal/vector_2D_s.F90 b/src/formal/vector_2D_s.F90 index a225e0a..7817b6d 100644 --- a/src/formal/vector_2D_s.F90 +++ b/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()) @@ -184,8 +184,8 @@ pure function description(coordinate, component) result(point_cloud) 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 +217,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))) diff --git a/test/scalar_2D_test_m.F90 b/test/scalar_2D_test_m.F90 index 8da5816..ef9fec8 100644 --- a/test/scalar_2D_test_m.F90 +++ b/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 From 0f3b983a3cd827e53cc0413d8a9a56c7c7f3afc6 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Thu, 9 Jul 2026 22:47:18 -0700 Subject: [PATCH 2/4] feat(vector_2D): {pre,post}_multiply_scalar --- example/2D-advection-diffusion.F90 | 13 ++++++++---- src/formal/tensors_2D_m.F90 | 21 ++++++++++++++++++- src/formal/vector_2D_s.F90 | 33 ++++++++++++++++++++++++++++++ 3 files changed, 62 insertions(+), 5 deletions(-) diff --git a/example/2D-advection-diffusion.F90 b/example/2D-advection-diffusion.F90 index 1916a51..2c6d1dd 100644 --- a/example/2D-advection-diffusion.F90 +++ b/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,10 @@ 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) - (v .dot. .grad. s) + !ds_dt = .div. (D * .grad. s) - .div. (v * s) + associate(ds_dt => .div.(v * s)) + end associate end function end program diff --git a/src/formal/tensors_2D_m.F90 b/src/formal/tensors_2D_m.F90 index 4936b0d..bc083ae 100644 --- a/src/formal/tensors_2D_m.F90 +++ b/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 @@ -146,6 +146,7 @@ pure module function construct_2D_scalar_from_components(tensor_2D, gradient_ope 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 @@ -155,6 +156,8 @@ pure module function construct_2D_scalar_from_components(tensor_2D, gradient_ope procedure, non_overridable, private :: vector_2D_conformable_scalar 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 @@ -412,6 +415,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 diff --git a/src/formal/vector_2D_s.F90 b/src/formal/vector_2D_s.F90 index 7817b6d..240c360 100644 --- a/src/formal/vector_2D_s.F90 +++ b/src/formal/vector_2D_s.F90 @@ -179,6 +179,39 @@ 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)) From 3520e1c80c15cf4abef8a3e1c914ad5623eacde2 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Thu, 9 Jul 2026 23:04:21 -0700 Subject: [PATCH 3/4] feat(divergence_2D_t): def homogeneous operator(-) --- src/formal/divergence_2D_s.F90 | 14 ++++++++++++++ src/formal/tensors_2D_m.F90 | 10 +++++++++- 2 files changed, 23 insertions(+), 1 deletion(-) diff --git a/src/formal/divergence_2D_s.F90 b/src/formal/divergence_2D_s.F90 index d78f5b0..a8e756a 100644 --- a/src/formal/divergence_2D_s.F90 +++ b/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()) diff --git a/src/formal/tensors_2D_m.F90 b/src/formal/tensors_2D_m.F90 index bc083ae..afa5a28 100644 --- a/src/formal/tensors_2D_m.F90 +++ b/src/formal/tensors_2D_m.F90 @@ -230,7 +230,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 @@ -239,6 +239,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 @@ -470,6 +471,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 From 5e9fbae563fb30d13fbc7d9edaaa4611258b1f76 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Thu, 9 Jul 2026 23:32:25 -0700 Subject: [PATCH 4/4] feat(scalar_2D_t): assgin divergence_2D_t --- example/2D-advection-diffusion.F90 | 5 +---- src/formal/scalar_2D_s.F90 | 28 ++++++++++++++++++++++++++++ src/formal/tensors_2D_m.F90 | 9 +++++++++ 3 files changed, 38 insertions(+), 4 deletions(-) diff --git a/example/2D-advection-diffusion.F90 b/example/2D-advection-diffusion.F90 index 2c6d1dd..d197640 100644 --- a/example/2D-advection-diffusion.F90 +++ b/example/2D-advection-diffusion.F90 @@ -75,10 +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) - associate(ds_dt => .div.(v * s)) - end associate + ds_dt = .div. (D * .grad. s) - .div. (v * s) end function end program diff --git a/src/formal/scalar_2D_s.F90 b/src/formal/scalar_2D_s.F90 index 9e2c9ff..96115bd 100644 --- a/src/formal/scalar_2D_s.F90 +++ b/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 diff --git a/src/formal/tensors_2D_m.F90 b/src/formal/tensors_2D_m.F90 index afa5a28..efba343 100644 --- a/src/formal/tensors_2D_m.F90 +++ b/src/formal/tensors_2D_m.F90 @@ -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 @@ -282,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