From ac990b60203df9bb949f2d27688cc5ad4e9c1932 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Wed, 13 May 2026 18:56:45 -0700 Subject: [PATCH 01/36] chore(tensors_1D_m): rm superfluous use statement --- src/formal/tensors_1D_m.F90 | 1 - 1 file changed, 1 deletion(-) diff --git a/src/formal/tensors_1D_m.F90 b/src/formal/tensors_1D_m.F90 index 42b2823..7ae8f45 100644 --- a/src/formal/tensors_1D_m.F90 +++ b/src/formal/tensors_1D_m.F90 @@ -7,7 +7,6 @@ module tensors_1D_m !! Define public 1D scalar and vector abstractions and associated mimetic gradient, !! divergence, and Laplacian operators as detailed by Corbino & Castillo (2020) !! https://doi.org/10.1016/j.cam.2019.06.042. - use julienne_m, only : file_t use differential_operators_1D_m, only : divergence_operator_1D_t, gradient_operator_1D_t implicit none From b764dfe44189c4b91003fa2e2e5a0eee7606e6ee Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Wed, 13 May 2026 19:00:29 -0700 Subject: [PATCH 02/36] chore(grad): rm redundant operator construction --- src/formal/scalar_1D_s.F90 | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/src/formal/scalar_1D_s.F90 b/src/formal/scalar_1D_s.F90 index cabd16f..92615e9 100644 --- a/src/formal/scalar_1D_s.F90 +++ b/src/formal/scalar_1D_s.F90 @@ -143,14 +143,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 From 98567ad0d5cfbd49482cf67231f8ce76078074db Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Wed, 13 May 2026 19:03:51 -0700 Subject: [PATCH 03/36] chore(div): rm redundant operator construction --- src/formal/vector_1D_s.F90 | 28 +++++++++++----------------- 1 file changed, 11 insertions(+), 17 deletions(-) diff --git a/src/formal/vector_1D_s.F90 b/src/formal/vector_1D_s.F90 index 8dc9e34..fd19586 100644 --- a/src/formal/vector_1D_s.F90 +++ b/src/formal/vector_1D_s.F90 @@ -69,25 +69,19 @@ pure module function construct_1D_vector_from_function(initializer, order, cells integer center -#ifdef NAGFOR - associate(D => self%divergence_operator_1D_) -#else - associate(D => (self%divergence_operator_1D_)) -#endif - associate(Dv => D .x. self%values_) - divergence_1D%tensor_1D_t = tensor_1D_t(Dv(2:size(Dv)-1), self%x_min_, self%x_max_, self%cells_, self%order_) + associate(Dv => self%divergence_operator_1D_ .x. self%values_) + divergence_1D%tensor_1D_t = tensor_1D_t(Dv(2:size(Dv)-1), self%x_min_, self%x_max_, self%cells_, self%order_) #if ASSERTIONS - associate( & - q => divergence_1D%weights() & - ,dx => (self%x_max_ - self%x_min_)/self%cells_ & - ,b => [-1D0, [(0D0, center = 1, self%cells_-1)], 1D0] & - ) - call_julienne_assert(.all. ([size(Dv), size(q)] .equalsExpected. self%cells_+2)) - call_julienne_assert((.all. (matmul(transpose(D%assemble()), q) .approximates. b/dx .within. double_equivalence))) - ! Check D^T * a = b_{m+1}, Eq. (19), Corbino & Castillo (2020) - end associate -#endif + associate( & + q => divergence_1D%weights() & + ,dx => (self%x_max_ - self%x_min_)/self%cells_ & + ,b => [-1D0, [(0D0, center = 1, self%cells_-1)], 1D0] & + ) + call_julienne_assert(.all. ([size(Dv), size(q)] .equalsExpected. self%cells_+2)) + call_julienne_assert((.all. (matmul(transpose(self%divergence_operator_1D_%assemble()), q) .approximates. b/dx .within. double_equivalence))) + ! Check D^T * a = b_{m+1}, Eq. (19), Corbino & Castillo (2020) end associate +#endif end associate end procedure From da440186ecc60954b4aebea8bb311ee5a41b36fe Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Wed, 13 May 2026 19:09:21 -0700 Subject: [PATCH 04/36] feat(differential_ops): mk constructors elemental --- src/formal/differential_operators_1D_m.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/formal/differential_operators_1D_m.F90 b/src/formal/differential_operators_1D_m.F90 index 41f7bc2..0b77038 100644 --- a/src/formal/differential_operators_1D_m.F90 +++ b/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 From 682843c51f895af56dd337cd29054168ee79f235 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Wed, 13 May 2026 19:58:39 -0700 Subject: [PATCH 05/36] refac(scalar_1D): rename & expose scalar grid calc This commit hoists the 1D scalar grid locations calculator up from scalar_1D_s to tensors_1D_m, renames it from "scalar_1D_grid_locations" to "cell_centers_extended_1D", and makes it public in anticipation of wider use in multidimensional calculations. --- src/formal/scalar_1D_s.F90 | 17 +++-------------- src/formal/tensors_1D_m.F90 | 16 ++++++++++++++-- 2 files changed, 17 insertions(+), 16 deletions(-) diff --git a/src/formal/scalar_1D_s.F90 b/src/formal/scalar_1D_s.F90 index 92615e9..109e5f3 100644 --- a/src/formal/scalar_1D_s.F90 +++ b/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) @@ -194,19 +194,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 \ No newline at end of file diff --git a/src/formal/tensors_1D_m.F90 b/src/formal/tensors_1D_m.F90 index 7ae8f45..49198d2 100644 --- a/src/formal/tensors_1D_m.F90 +++ b/src/formal/tensors_1D_m.F90 @@ -22,6 +22,7 @@ module tensors_1D_m public :: vector_1D_initializer_i public :: d_dx public :: d2_dx2 + public :: cell_centers_extended_1D abstract interface @@ -495,10 +496,21 @@ pure module function postmultiply_scalar_1D(divergence_1D, scalar_1D) result(sca end interface -#ifndef __GFORTRAN__ - contains + pure function cell_centers_extended_1D(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 + +#ifndef __GFORTRAN__ + pure function cell_center_locations(x_min, x_max, cells) result(x) double precision, intent(in) :: x_min, x_max integer, intent(in) :: cells From e97f28c93993281a740b2cceb06d3e22e96d26e9 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Wed, 13 May 2026 20:23:26 -0700 Subject: [PATCH 06/36] refac(vector_1D): rename & expose vector grid calc This commit hoists the 1D vector grid locations calculator up from vector_1D_s to tensors_1D_m, renames it from "faces" to "faces_1D", and makes it public in anticipation of wider use in multidimensional calculations. --- src/formal/tensors_1D_m.F90 | 12 ++++++++++++ src/formal/vector_1D_s.F90 | 17 +++-------------- 2 files changed, 15 insertions(+), 14 deletions(-) diff --git a/src/formal/tensors_1D_m.F90 b/src/formal/tensors_1D_m.F90 index 49198d2..bc2d740 100644 --- a/src/formal/tensors_1D_m.F90 +++ b/src/formal/tensors_1D_m.F90 @@ -23,6 +23,7 @@ module tensors_1D_m public :: d_dx public :: d2_dx2 public :: cell_centers_extended_1D + public :: faces_1D abstract interface @@ -509,6 +510,17 @@ pure function cell_centers_extended_1D(x_min, x_max, cells) result(x) end associate end function + pure function faces_1D(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, x_min + [(cell*dx, cell = 1, cells-1)], x_max] + end associate + end function + #ifndef __GFORTRAN__ pure function cell_center_locations(x_min, x_max, cells) result(x) diff --git a/src/formal/vector_1D_s.F90 b/src/formal/vector_1D_s.F90 index fd19586..1be73da 100644 --- a/src/formal/vector_1D_s.F90 +++ b/src/formal/vector_1D_s.F90 @@ -33,7 +33,7 @@ call_julienne_assert(x_max .greaterThan. x_min) call_julienne_assert(cells .isAtLeast. 2*order+1) - associate(values => initializer(faces(x_min, x_max, cells))) + associate(values => initializer(faces_1D(x_min, x_max, cells))) vector_1D%tensor_1D_t = tensor_1D_t(values, x_min, x_max, cells, order) end associate vector_1D%divergence_operator_1D_ = divergence_operator_1D_t(k=order, dx=(x_max - x_min)/cells, cells=cells) @@ -52,7 +52,7 @@ pure module function construct_1D_vector_from_function(initializer, order, cells call_julienne_assert(x_max .greaterThan. x_min) call_julienne_assert(cells .isAtLeast. 2*order+1) - associate(values => initializer(faces(x_min, x_max, cells))) + associate(values => initializer(faces_1D(x_min, x_max, cells))) vector_1D%tensor_1D_t = tensor_1D_t(values, x_min, x_max, cells, order) end associate vector_1D%divergence_operator_1D_ = divergence_operator_1D_t(k=order, dx=(x_max - x_min)/cells, cells=cells) @@ -90,19 +90,8 @@ pure module function construct_1D_vector_from_function(initializer, order, cells face_centered_values = self%values_ end procedure - pure function faces(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, x_min + [(cell*dx, cell = 1, cells-1)], x_max] - end associate - end function - module procedure vector_1D_grid - cell_faces = faces(self%x_min_, self%x_max_, self%cells_) + cell_faces = faces_1D(self%x_min_, self%x_max_, self%cells_) end procedure module procedure weighted_premultiply From bf1f66a322178863c65d62a489f9281cd75c71e2 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Wed, 13 May 2026 20:33:19 -0700 Subject: [PATCH 07/36] feat: add {tensor,scalar,vector}_2D_t types --- src/formal/scalar_2D_s.F90 | 60 +++++++++++++++++ src/formal/tensors_2D_m.F90 | 124 ++++++++++++++++++++++++++++++++++++ src/formal/vector_2D_s.F90 | 41 ++++++++++++ src/formal_m.f90 | 5 ++ 4 files changed, 230 insertions(+) create mode 100644 src/formal/scalar_2D_s.F90 create mode 100644 src/formal/tensors_2D_m.F90 create mode 100644 src/formal/vector_2D_s.F90 diff --git a/src/formal/scalar_2D_s.F90 b/src/formal/scalar_2D_s.F90 new file mode 100644 index 0000000..4b30cb5 --- /dev/null +++ b/src/formal/scalar_2D_s.F90 @@ -0,0 +1,60 @@ +! 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 + implicit none + +contains + + 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( & + x1 => cell_centers_extended_1D(x_min(1), x_max(1), cells(1)) & + ,x2 => cell_centers_extended_1D(x_min(2), x_max(2), cells(2)) & + ) + allocate(scalar_2D%values_(cells(1)+2, cells(2)+2,1)) + + do concurrent(integer :: i=1:cells(1)+2, j=1:cells(2)+2) default(none) shared(scalar_2D, x1, x2) + scalar_2D%values_(i,j,1) = initializer(x1(i), x2(j)) + end do + end associate + + scalar_2D%order_ = order + scalar_2D%x_min_ = x_min + scalar_2D%x_max_ = x_max + scalar_2D%cells_ = cells + scalar_2D%gradient_operator_1D_ = gradient_operator_1D_t(k=order, dx=(x_max - x_min)/cells, cells=cells) + end procedure + + module procedure grad + + integer c + + associate(dx => (self%x_max_ - self%x_min_)/self%cells_) + +! gradient_2D%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_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 + +end submodule scalar_2D_s diff --git a/src/formal/tensors_2D_m.F90 b/src/formal/tensors_2D_m.F90 new file mode 100644 index 0000000..dcd0a40 --- /dev/null +++ b/src/formal/tensors_2D_m.F90 @@ -0,0 +1,124 @@ +! Copyright (c) 2026, The Regents of the University of California +! Terms of use are as specified in LICENSE.txt + +#include "formal-language-support.F90" + +module tensors_2D_m + !! Define public 2D scalar and vector abstractions and associated mimetic gradient, + !! divergence, and Laplacian operators as detailed by Corbino & Castillo (2020) + !! https://doi.org/10.1016/j.cam.2019.06.042. + use differential_operators_1D_m, only : gradient_operator_1D_t, divergence_operator_1D_t + + implicit none + + private + + public :: scalar_2D_t + public :: vector_2D_t + public :: gradient_2D_t + public :: scalar_2D_initializer_i + public :: vector_2D_initializer_i + + integer, parameter :: space_dimension = 2 + + abstract interface + + pure function scalar_2D_initializer_i(x1, x2) result(f) + !! Sampling function for initializing a scalar_2D_t object + implicit none + double precision, intent(in) :: x1, x2 + double precision, allocatable :: f + end function + + pure function vector_2D_initializer_i(x1, x2 ) result(v) + !! Sampling function for initializing a vector_2D_t object + import space_dimension + implicit none + double precision, intent(in) :: x1, x2 + double precision v(space_dimension) + end function + + end interface + + type tensor_2D_t + !! Encapsulate the components that are common to all 2D tensors. + !! Child types define the operations supported by each child, including + !! gradient (.grad.) for scalars and divergence (.div.) for vectors. + private + double precision x_min_(space_dimension) !! domain lower boundary + double precision x_max_(space_dimension) !! domain upper boundary + integer cells_(space_dimension) !! number of grid cells spanning the domain + integer order_ !! order of accuracy of mimetic discretization + double precision, allocatable :: values_(:,:,:) !! tensor components at spatial locations + end type + + type, extends(tensor_2D_t) :: scalar_2D_t + !! Encapsulate scalar values at cell centers and boundaries + private + type(gradient_operator_1D_t) gradient_operator_1D_(space_dimension) + contains + generic :: operator(.grad.) => grad + generic :: values => scalar_2D_values + procedure, non_overridable, private :: grad + procedure, non_overridable, private :: scalar_2D_values + end type + + interface scalar_2D_t + + pure module function construct_2D_scalar_from_function(initializer, order, cells, x_min, x_max) result(scalar_2D) + !! Result is a collection of cell-centered-extended values with a corresponding mimetic gradient operator + implicit none + procedure(scalar_2D_initializer_i), pointer :: initializer + integer, intent(in) :: order !! order of accuracy + integer, intent(in) :: cells(:) !! number of grid cells spanning each spatial direction + double precision, intent(in) :: x_min(:) !! grid location minima + double precision, intent(in) :: x_max(:) !! grid location maxima + type(scalar_2D_t) scalar_2D + end function + + end interface + + type, extends(tensor_2D_t) :: vector_2D_t + !! Encapsulate 2D vector values at cell faces (of unit area for 2D) and corresponding operators + private + end type + + interface vector_2D_t + + pure module function construct_2D_vector_from_function(initializer, order, cells, x_min, x_max) result(vector_2D) + !! Result is a 2D vector with values initialized by the provided procedure pointer sampled on the specified + !! number of evenly spaced cells covering [x_min, x_max] + implicit none + procedure(vector_2D_initializer_i), pointer :: initializer + integer, intent(in) :: order !! order of accuracy + integer, intent(in) :: cells(:) !! number of grid cells spanning each spatial direction + double precision, intent(in) :: x_min(:) !! grid location minima + double precision, intent(in) :: x_max(:) !! grid location maxima + type(vector_2D_t) vector_2D + end function + + end interface + + type, extends(vector_2D_t) :: gradient_2D_t + !! A 2D mimetic gradient vector field abstraction with a public method that produces corresponding numerical quadrature weights + type(divergence_operator_1D_t) divergence_operator_1D_(space_dimension) + end type + + interface + + pure module function scalar_2D_values(self) result(scalar_values) + !! Scalar values getter + class(scalar_2D_t), intent(in) :: self + double precision, allocatable :: scalar_values(:,:) + end function + + pure module function grad(self) result(gradient_2D) + !! Result is mimetic gradient of the scalar_2D_t "self" + implicit none + class(scalar_2D_t), intent(in) :: self + type(gradient_2D_t) gradient_2D + end function + + end interface + +end module tensors_2D_m diff --git a/src/formal/vector_2D_s.F90 b/src/formal/vector_2D_s.F90 new file mode 100644 index 0000000..4f66e2a --- /dev/null +++ b/src/formal/vector_2D_s.F90 @@ -0,0 +1,41 @@ +! 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) vector_2D_s + use julienne_m, only : & + call_julienne_assert_ & + ,operator(.all.) & + ,operator(.equalsExpected.) & + ,operator(.greaterThan.) & + ,operator(.isAtLeast.) + use tensors_1D_m, only : faces_1D + implicit none + +contains + + module procedure construct_2D_vector_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( & + x1 => faces_1D(x_min(1), x_max(1), cells(1)) & + ,x2 => faces_1D(x_min(2), x_max(2), cells(2)) & + ) + allocate(vector_2D%values_(cells(1)+1, cells(2)+1, space_dimension)) + + do concurrent(integer :: i=1:cells(1)+1, j=1:cells(2)+1, dir=1:space_dimension) default(none) shared(vector_2D, x1, x2) + vector_2D%values_(i,j,:) = initializer(x1(i), x2(j)) + end do + end associate + + vector_2D%order_ = order + vector_2D%x_min_ = x_min + vector_2D%x_max_ = x_max + vector_2D%cells_ = cells + end procedure + +end submodule vector_2D_s diff --git a/src/formal_m.f90 b/src/formal_m.f90 index 2fc9526..d03f6bd 100644 --- a/src/formal_m.f90 +++ b/src/formal_m.f90 @@ -18,6 +18,11 @@ module formal_m ,d_dx & ! scalar_1D_t spatial derivative ,d2_dx2 ! scalar_1D_t spatial derivative + use tensors_2D_m !, only : & + ! scalar_2D_t & ! discrete 2D scalar field derived type + ! ,vector_2D_t & ! discrete 2D vector field derived type + ! ,gradient_2D_t & ! result of `.grad. s` for a scalar_2D_t s + use differential_operators_1D_m, only : & gradient_operator_1D_t & ! matrix operator defining a 1D mimetic gradient ,divergence_operator_1D_t ! matrix operator defining a 1D mimetic divergence From 08e119e579da1aa78f0c8afa939880b028d12080 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Sat, 16 May 2026 17:53:24 -0700 Subject: [PATCH 08/36] feat(tensor_2D): .grad.scalar_2D_t->gradient_2D_t The commit contains the first passing test of a 2D differential operator: .grad. correctly computes a gradient_2D_t when given a scalar_2D_t operand. --- src/formal/scalar_1D_s.F90 | 16 +++++++ src/formal/scalar_2D_s.F90 | 70 +++++++++++++++++++++---------- src/formal/tensor_2D_s.F90 | 14 +++++++ src/formal/tensors_1D_m.F90 | 11 +++++ src/formal/tensors_2D_m.F90 | 76 ++++++++++++++++++++++++++++----- src/formal/vector_2D_s.F90 | 48 +++++++++++++++------ test/driver.f90 | 2 + test/scalar_2D_test_m.F90 | 84 +++++++++++++++++++++++++++++++++++++ 8 files changed, 275 insertions(+), 46 deletions(-) create mode 100644 src/formal/tensor_2D_s.F90 create mode 100644 test/scalar_2D_test_m.F90 diff --git a/src/formal/scalar_1D_s.F90 b/src/formal/scalar_1D_s.F90 index 109e5f3..a7f0b06 100644 --- a/src/formal/scalar_1D_s.F90 +++ b/src/formal/scalar_1D_s.F90 @@ -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_ & diff --git a/src/formal/scalar_2D_s.F90 b/src/formal/scalar_2D_s.F90 index 4b30cb5..b42262d 100644 --- a/src/formal/scalar_2D_s.F90 +++ b/src/formal/scalar_2D_s.F90 @@ -10,49 +10,73 @@ ,operator(.equalsExpected.) & ,operator(.greaterThan.) & ,operator(.isAtLeast.) - use tensors_1D_m, only : cell_centers_extended_1D + use tensors_1D_m, only : cell_centers_extended_1D, scalar_1D_t implicit none contains + module procedure scalar_2D_values + scalar_values = self%values_(:,:,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( & - x1 => cell_centers_extended_1D(x_min(1), x_max(1), cells(1)) & - ,x2 => cell_centers_extended_1D(x_min(2), x_max(2), cells(2)) & - ) - allocate(scalar_2D%values_(cells(1)+2, cells(2)+2,1)) - - do concurrent(integer :: i=1:cells(1)+2, j=1:cells(2)+2) default(none) shared(scalar_2D, x1, x2) - scalar_2D%values_(i,j,1) = initializer(x1(i), x2(j)) - end do + 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]) & + ,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 - scalar_2D%order_ = order - scalar_2D%x_min_ = x_min - scalar_2D%x_max_ = x_max - scalar_2D%cells_ = cells - scalar_2D%gradient_operator_1D_ = gradient_operator_1D_t(k=order, dx=(x_max - x_min)/cells, cells=cells) + 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 grad - integer c + integer c, i, j - associate(dx => (self%x_max_ - self%x_min_)/self%cells_) + gradient_2D%x_min_ = self%x_min_ + gradient_2D%x_max_ = self%x_max_ + gradient_2D%cells_ = self%cells_ + gradient_2D%order_ = self%order_ -! gradient_2D%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_) + allocate(gradient_2D%values_(self%cells_(1)+1, self%cells_(2)+1, space_dimension)) - gradient_2D%divergence_operator_1D_ = divergence_operator_1D_t(self%order_, dx, self%cells_) + gradient_x_component: & + do concurrent(integer :: j=1:self%cells_(2)+2) default(none) shared(gradient_2D, self) + gradient_2D%values_(:,j,1) = self%gradient_operator_1D_(1) .x. self%values_(:,j,1) + end do gradient_x_component - !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 + gradient_y_component: & + do concurrent(integer :: i=1:self%cells_(1)+2) default(none) shared(gradient_2D, self) + gradient_2D%values_(i,:,2) = self%gradient_operator_1D_(2) .x. self%values_(i,:,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 diff --git a/src/formal/tensor_2D_s.F90 b/src/formal/tensor_2D_s.F90 new file mode 100644 index 0000000..cf2a2df --- /dev/null +++ b/src/formal/tensor_2D_s.F90 @@ -0,0 +1,14 @@ +submodule(tensors_2D_m) tensor_2D_s + implicit none + +contains + + module procedure construct_2D_tensor_from_components + tensor_2D%values_ = values + tensor_2D%cells_ = cells + tensor_2D%x_min_ = x_min + tensor_2D%x_max_ = x_max + tensor_2D%order_ = order + end procedure + +end submodule \ No newline at end of file diff --git a/src/formal/tensors_1D_m.F90 b/src/formal/tensors_1D_m.F90 index bc2d740..b0f4dd6 100644 --- a/src/formal/tensors_1D_m.F90 +++ b/src/formal/tensors_1D_m.F90 @@ -120,6 +120,17 @@ pure module function construct_1D_scalar_from_function(initializer, order, cells type(scalar_1D_t) scalar_1D end function + pure module function construct_1D_scalar_constant(constant, order, cells, x_min, x_max) result(scalar_1D) + !! Result is a collection of cell-centered-extended values with a corresponding mimetic gradient operator + implicit none + double precision, intent(in) :: constant !! scalar value + integer, intent(in) :: order !! order of accuracy + integer, intent(in) :: cells !! number of grid cells spanning the domain + double precision, intent(in) :: x_min !! grid location minimum + double precision, intent(in) :: x_max !! grid location maximum + type(scalar_1D_t) scalar_1D + end function + pure module function construct_1D_scalar_from_parent(tensor_1D) result(scalar_1D) !! Result is a 1D vector with the provided parent component tensor_1D and the provided divergence operatror type(tensor_1D_t), intent(in) :: tensor_1D diff --git a/src/formal/tensors_2D_m.F90 b/src/formal/tensors_2D_m.F90 index dcd0a40..6473bfa 100644 --- a/src/formal/tensors_2D_m.F90 +++ b/src/formal/tensors_2D_m.F90 @@ -1,8 +1,6 @@ ! Copyright (c) 2026, The Regents of the University of California ! Terms of use are as specified in LICENSE.txt -#include "formal-language-support.F90" - module tensors_2D_m !! Define public 2D scalar and vector abstractions and associated mimetic gradient, !! divergence, and Laplacian operators as detailed by Corbino & Castillo (2020) @@ -23,19 +21,19 @@ module tensors_2D_m abstract interface - pure function scalar_2D_initializer_i(x1, x2) result(f) + pure function scalar_2D_initializer_i(x,y) result(f) !! Sampling function for initializing a scalar_2D_t object implicit none - double precision, intent(in) :: x1, x2 - double precision, allocatable :: f + double precision, intent(in) :: x(:), y(:) + double precision f(size(x),size(y)) end function - pure function vector_2D_initializer_i(x1, x2 ) result(v) + pure function vector_2D_initializer_i(x,y) result(v) !! Sampling function for initializing a vector_2D_t object import space_dimension implicit none - double precision, intent(in) :: x1, x2 - double precision v(space_dimension) + double precision, intent(in) :: x(:), y(:) + double precision v(size(x),size(y),space_dimension) end function end interface @@ -45,13 +43,27 @@ pure function vector_2D_initializer_i(x1, x2 ) result(v) !! Child types define the operations supported by each child, including !! gradient (.grad.) for scalars and divergence (.div.) for vectors. private + double precision, allocatable :: values_(:,:,:) !! tensor components at 2D spatial locations double precision x_min_(space_dimension) !! domain lower boundary double precision x_max_(space_dimension) !! domain upper boundary integer cells_(space_dimension) !! number of grid cells spanning the domain integer order_ !! order of accuracy of mimetic discretization - double precision, allocatable :: values_(:,:,:) !! tensor components at spatial locations end type + interface tensor_2D_t + + pure module function construct_2D_tensor_from_components(values, cells, x_min, x_max, order) result(tensor_2D) + implicit none + double precision, intent(in) :: values(:,:,:) !! tensor components at 2D spatial locations + double precision, intent(in) :: x_min(:) !! domain lower boundary + double precision, intent(in) :: x_max(:) !! domain upper boundary + integer, intent(in) :: cells(:) !! number of grid cells spanning the domain + integer, intent(in) :: order !! order of accuracy of mimetic discretization + type(tensor_2D_t) tensor_2D + end function + + end interface + type, extends(tensor_2D_t) :: scalar_2D_t !! Encapsulate scalar values at cell centers and boundaries private @@ -59,8 +71,10 @@ pure function vector_2D_initializer_i(x1, x2 ) result(v) contains generic :: operator(.grad.) => grad generic :: values => scalar_2D_values + generic :: grid => scalar_2D_grid procedure, non_overridable, private :: grad procedure, non_overridable, private :: scalar_2D_values + procedure, non_overridable, private :: scalar_2D_grid end type interface scalar_2D_t @@ -76,11 +90,23 @@ pure module function construct_2D_scalar_from_function(initializer, order, cells type(scalar_2D_t) scalar_2D end function + pure module function construct_2D_scalar_from_mold(initializer, mold) result(scalar_2D) + !! Result is a 2D scalar field using a mold for all components other than the field values + implicit none + procedure(scalar_2D_initializer_i), pointer :: initializer + type(scalar_2D_t), intent(in) :: mold + type(scalar_2D_t) scalar_2D + end function + end interface type, extends(tensor_2D_t) :: vector_2D_t !! Encapsulate 2D vector values at cell faces (of unit area for 2D) and corresponding operators private + type(divergence_operator_1D_t) divergence_operator_1D_(space_dimension) + contains + generic :: values => vector_2D_values + procedure, non_overridable, private :: vector_2D_values end type interface vector_2D_t @@ -96,12 +122,29 @@ pure module function construct_2D_vector_from_function(initializer, order, cells double precision, intent(in) :: x_max(:) !! grid location maxima type(vector_2D_t) vector_2D end function + + pure module function construct_2D_vector_from_vector_mold(initializer, mold) result(vector_2D) + !! Result is a 2D vector with values initialized by the provided procedure pointer sampled on the + !! same grid as the mold + implicit none + procedure(vector_2D_initializer_i), pointer :: initializer + type(vector_2D_t), intent(in) :: mold + type(vector_2D_t) vector_2D + end function + pure module function construct_2D_vector_from_scalar_mold(initializer, mold) result(vector_2D) + !! Result is a 2D vector with values initialized by the provided procedure pointer sampled on the + !! face-centered grid corresponding to the cell-centered grid of the mold + implicit none + procedure(vector_2D_initializer_i), pointer :: initializer + type(scalar_2D_t), intent(in) :: mold + type(vector_2D_t) vector_2D + end function + end interface type, extends(vector_2D_t) :: gradient_2D_t !! A 2D mimetic gradient vector field abstraction with a public method that produces corresponding numerical quadrature weights - type(divergence_operator_1D_t) divergence_operator_1D_(space_dimension) end type interface @@ -112,6 +155,19 @@ pure module function scalar_2D_values(self) result(scalar_values) double precision, allocatable :: scalar_values(:,:) end function + pure module function scalar_2D_grid(self, direction) result(scalar_grid_1D) + !! Result array contains scalar grid locations along the requested spatial direction + class(scalar_2D_t), intent(in) :: self + integer, intent(in) :: direction + double precision, allocatable :: scalar_grid_1D(:) + end function + + pure module function vector_2D_values(self) result(vector_values) + !! Vector values getter + class(vector_2D_t), intent(in) :: self + double precision, allocatable :: vector_values(:,:,:) + end function + pure module function grad(self) result(gradient_2D) !! Result is mimetic gradient of the scalar_2D_t "self" implicit none diff --git a/src/formal/vector_2D_s.F90 b/src/formal/vector_2D_s.F90 index 4f66e2a..3f8aff6 100644 --- a/src/formal/vector_2D_s.F90 +++ b/src/formal/vector_2D_s.F90 @@ -11,31 +11,53 @@ ,operator(.greaterThan.) & ,operator(.isAtLeast.) use tensors_1D_m, only : faces_1D + use differential_operators_1D_m, only : divergence_operator_1D_t implicit none contains module procedure construct_2D_vector_from_function + integer dir + 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( & - x1 => faces_1D(x_min(1), x_max(1), cells(1)) & - ,x2 => faces_1D(x_min(2), x_max(2), cells(2)) & - ) - allocate(vector_2D%values_(cells(1)+1, cells(2)+1, space_dimension)) + associate(x => faces_1D(x_min(1), x_max(1), cells(1)), y => faces_1D(x_min(2), x_max(2), cells(2))) + vector_2D%tensor_2D_t = tensor_2D_t(values = initializer(x,y), cells = cells, x_min = x_min, x_max = x_max, order = order) + vector_2D%divergence_operator_1D_ = [(divergence_operator_1D_t(k=order, dx=((x_max(dir)-x_min(dir))/cells(dir)), cells=cells(dir)), dir=1,space_dimension)] + end associate + end procedure + + module procedure construct_2D_vector_from_vector_mold + integer dir + + call_julienne_assert(.all. ([size(mold%cells_), size(mold%x_min_), size(mold%x_max_)] .equalsExpected. space_dimension)) + call_julienne_assert(.all. (mold%x_max_ .greaterThan. mold%x_min_)) + call_julienne_assert(.all. (mold%cells_ .isAtLeast. 2*mold%order_)) - do concurrent(integer :: i=1:cells(1)+1, j=1:cells(2)+1, dir=1:space_dimension) default(none) shared(vector_2D, x1, x2) - vector_2D%values_(i,j,:) = initializer(x1(i), x2(j)) - end do + associate(x => faces_1D(mold%x_min_(1), mold%x_max_(1), mold%cells_(1)), y => faces_1D(mold%x_min_(2), mold%x_max_(2), mold%cells_(2))) + vector_2D%tensor_2D_t = tensor_2D_t(values = initializer(x,y), cells = mold%cells_, x_min = mold%x_min_, x_max = mold%x_max_, order = mold%order_) + vector_2D%divergence_operator_1D_ = [(divergence_operator_1D_t(k=mold%order_, dx=((mold%x_max_(dir)-mold%x_min_(dir))/mold%cells_(dir)), cells=mold%cells_(dir)), dir=1,space_dimension)] end associate + end procedure + + module procedure construct_2D_vector_from_scalar_mold + integer dir + + call_julienne_assert(.all. ([size(mold%cells_), size(mold%x_min_), size(mold%x_max_)] .equalsExpected. space_dimension)) + call_julienne_assert(.all. (mold%x_max_ .greaterThan. mold%x_min_)) + call_julienne_assert(.all. (mold%cells_ .isAtLeast. 2*mold%order_)) + + associate(x => faces_1D(mold%x_min_(1), mold%x_max_(1), mold%cells_(1)), y => faces_1D(mold%x_min_(2), mold%x_max_(2), mold%cells_(2))) + vector_2D%tensor_2D_t = tensor_2D_t(values = initializer(x,y), cells = mold%cells_, x_min = mold%x_min_, x_max = mold%x_max_, order = mold%order_) + vector_2D%divergence_operator_1D_ = [(divergence_operator_1D_t(k=mold%order_, dx=((mold%x_max_(dir)-mold%x_min_(dir))/mold%cells_(dir)), cells=mold%cells_(dir)), dir=1,space_dimension)] + end associate + end procedure - vector_2D%order_ = order - vector_2D%x_min_ = x_min - vector_2D%x_max_ = x_max - vector_2D%cells_ = cells + module procedure vector_2D_values + vector_values = self%values_ end procedure -end submodule vector_2D_s +end submodule vector_2D_s \ No newline at end of file diff --git a/test/driver.f90 b/test/driver.f90 index 65809b8..1464cc7 100644 --- a/test/driver.f90 +++ b/test/driver.f90 @@ -10,6 +10,7 @@ program test_suite_driver use integration_operators_1D_test_m, only : integration_operators_1D_test_t 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 implicit none associate(test_harness => test_harness_t([ & @@ -19,6 +20,7 @@ program test_suite_driver ,test_fixture_t(integration_operators_1D_test_t()) & ,test_fixture_t(interpolator_1D_test_t()) & ,test_fixture_t(scalar_1D_test_t()) & + ,test_fixture_t(scalar_2D_test_t()) & ])) call test_harness%report_results end associate diff --git a/test/scalar_2D_test_m.F90 b/test/scalar_2D_test_m.F90 new file mode 100644 index 0000000..45c1d4f --- /dev/null +++ b/test/scalar_2D_test_m.F90 @@ -0,0 +1,84 @@ +! Copyright (c) 2026, The Regents of the University of California +! Terms of use are as specified in LICENSE.txt + +module scalar_2D_test_m + use julienne_m, only : & + operator(//) & + ,operator(.all.) & + ,operator(.also.) & + ,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, scalar_2D_initializer_i, vector_2D_initializer_i + + implicit none + + type, extends(test_t) :: scalar_2D_test_t + contains + procedure, nopass :: subject + procedure, nopass :: results + end type + + double precision, parameter :: tolerance = 2D-11 + integer, parameter :: space_dimension = 2 + +contains + + pure function subject() result(test_subject) + character(len=:), allocatable :: test_subject + test_subject = 'The scalar_2D_t derived type' + end function + + function results() result(test_results) + type(scalar_2D_test_t) scalar_2D_test + type(test_result_t), allocatable :: test_results(:) + + test_results = scalar_2D_test%run([ & + test_description_t('computing the gradient of a scalar field', usher(check_gradient)) & + ]) + end function + + pure function bilinear(x,y) result(z) + double precision, intent(in) :: x(:), y(:) + double precision z(size(x),size(y)) + do concurrent(integer :: j=1:size(y)) default(none) shared(x,y,z) + z(:,j) = 2*x + 3*y(j) + end do + end function + + pure function bilinear_gradient(x,y) result(gradient) + double precision, intent(in) :: x(:), y(:) + double precision gradient(size(x),size(y),space_dimension) + do concurrent(integer :: i=1:size(x), j=1:size(y)) + gradient(i,j,:) = [2,3] + end do + end function + + function check_gradient() result(test_diagnosis) + type(test_diagnosis_t) test_diagnosis + procedure(scalar_2D_initializer_i), pointer :: scalar_2D_initializer + procedure(vector_2D_initializer_i), pointer :: expected_gradient_initializer + integer order + + scalar_2D_initializer => bilinear + expected_gradient_initializer => bilinear_gradient + test_diagnosis = passing_test() + + do order = 2, 4, 2 + associate(scalar_2D => scalar_2D_t(scalar_2D_initializer, order=order, cells=[10,10], x_min=[0D0,0D0], x_max=[10D0,10D0])) + associate(grad_scalar => .grad. scalar_2D, expected_gradient => vector_2D_t(expected_gradient_initializer, mold=scalar_2D)) + test_diagnosis = test_diagnosis .also. & + .all. (grad_scalar%values() .approximates. expected_gradient%values() .within. tolerance) & + // string_t(" for order ") // string_t(order) + end associate + end associate + end do + end function + +end module scalar_2D_test_m \ No newline at end of file From 70adca5f4b8c48919a715e7d9998a3a16494fbec Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Sun, 17 May 2026 01:05:21 -0700 Subject: [PATCH 09/36] feat(tensor_2D): mk values_ component 6D This commit increases the rank of the tensor_2D_t values_ component to facilitate storing tensors of rank up to and including rank 4. --- src/formal/scalar_2D_s.F90 | 10 +++++----- src/formal/tensors_2D_m.F90 | 4 ++-- src/formal/vector_2D_s.F90 | 23 +++++++++++++++++++---- 3 files changed, 26 insertions(+), 11 deletions(-) diff --git a/src/formal/scalar_2D_s.F90 b/src/formal/scalar_2D_s.F90 index b42262d..0e9b44a 100644 --- a/src/formal/scalar_2D_s.F90 +++ b/src/formal/scalar_2D_s.F90 @@ -16,7 +16,7 @@ contains module procedure scalar_2D_values - scalar_values = self%values_(:,:,1) + scalar_values = self%values_(:,:,1,1,1,1) end procedure module procedure scalar_2D_grid @@ -39,7 +39,7 @@ 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]) & + 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) @@ -59,16 +59,16 @@ gradient_2D%cells_ = self%cells_ gradient_2D%order_ = self%order_ - allocate(gradient_2D%values_(self%cells_(1)+1, self%cells_(2)+1, space_dimension)) + 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:self%cells_(2)+2) default(none) shared(gradient_2D, self) - gradient_2D%values_(:,j,1) = self%gradient_operator_1D_(1) .x. self%values_(:,j,1) + 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:self%cells_(1)+2) default(none) shared(gradient_2D, self) - gradient_2D%values_(i,:,2) = self%gradient_operator_1D_(2) .x. self%values_(i,:,1) + 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_) diff --git a/src/formal/tensors_2D_m.F90 b/src/formal/tensors_2D_m.F90 index 6473bfa..e0ed9de 100644 --- a/src/formal/tensors_2D_m.F90 +++ b/src/formal/tensors_2D_m.F90 @@ -43,7 +43,7 @@ pure function vector_2D_initializer_i(x,y) result(v) !! Child types define the operations supported by each child, including !! gradient (.grad.) for scalars and divergence (.div.) for vectors. private - double precision, allocatable :: values_(:,:,:) !! tensor components at 2D spatial locations + double precision, allocatable :: values_(:,:,:,:,:,:) !! tensor components for rank<=4 at 2D locations double precision x_min_(space_dimension) !! domain lower boundary double precision x_max_(space_dimension) !! domain upper boundary integer cells_(space_dimension) !! number of grid cells spanning the domain @@ -54,7 +54,7 @@ pure function vector_2D_initializer_i(x,y) result(v) pure module function construct_2D_tensor_from_components(values, cells, x_min, x_max, order) result(tensor_2D) implicit none - double precision, intent(in) :: values(:,:,:) !! tensor components at 2D spatial locations + double precision, intent(in) :: values(:,:,:,:,:,:) !! tensor components at 2D spatial locations double precision, intent(in) :: x_min(:) !! domain lower boundary double precision, intent(in) :: x_max(:) !! domain upper boundary integer, intent(in) :: cells(:) !! number of grid cells spanning the domain diff --git a/src/formal/vector_2D_s.F90 b/src/formal/vector_2D_s.F90 index 3f8aff6..09a5e5b 100644 --- a/src/formal/vector_2D_s.F90 +++ b/src/formal/vector_2D_s.F90 @@ -25,7 +25,12 @@ call_julienne_assert(.all. (cells .isAtLeast. 2*order)) associate(x => faces_1D(x_min(1), x_max(1), cells(1)), y => faces_1D(x_min(2), x_max(2), cells(2))) - vector_2D%tensor_2D_t = tensor_2D_t(values = initializer(x,y), cells = cells, x_min = x_min, x_max = x_max, order = order) + associate(vector_values => initializer(x,y)) + vector_2D%tensor_2D_t = tensor_2D_t( & + values = reshape(vector_values, shape=[shape(vector_values),1,1,1]) & + ,cells = cells, x_min = x_min, x_max = x_max, order = order & + ) + end associate vector_2D%divergence_operator_1D_ = [(divergence_operator_1D_t(k=order, dx=((x_max(dir)-x_min(dir))/cells(dir)), cells=cells(dir)), dir=1,space_dimension)] end associate end procedure @@ -38,7 +43,12 @@ call_julienne_assert(.all. (mold%cells_ .isAtLeast. 2*mold%order_)) associate(x => faces_1D(mold%x_min_(1), mold%x_max_(1), mold%cells_(1)), y => faces_1D(mold%x_min_(2), mold%x_max_(2), mold%cells_(2))) - vector_2D%tensor_2D_t = tensor_2D_t(values = initializer(x,y), cells = mold%cells_, x_min = mold%x_min_, x_max = mold%x_max_, order = mold%order_) + associate(vector_values => initializer(x,y)) + vector_2D%tensor_2D_t = tensor_2D_t( & + values = reshape(vector_values, shape=[shape(vector_values),1,1,1]) & + ,cells = mold%cells_, x_min = mold%x_min_, x_max = mold%x_max_, order = mold%order_ & + ) + end associate vector_2D%divergence_operator_1D_ = [(divergence_operator_1D_t(k=mold%order_, dx=((mold%x_max_(dir)-mold%x_min_(dir))/mold%cells_(dir)), cells=mold%cells_(dir)), dir=1,space_dimension)] end associate end procedure @@ -51,13 +61,18 @@ call_julienne_assert(.all. (mold%cells_ .isAtLeast. 2*mold%order_)) associate(x => faces_1D(mold%x_min_(1), mold%x_max_(1), mold%cells_(1)), y => faces_1D(mold%x_min_(2), mold%x_max_(2), mold%cells_(2))) - vector_2D%tensor_2D_t = tensor_2D_t(values = initializer(x,y), cells = mold%cells_, x_min = mold%x_min_, x_max = mold%x_max_, order = mold%order_) + associate(vector_values => initializer(x,y)) + vector_2D%tensor_2D_t = tensor_2D_t( & + values = reshape(vector_values, shape=[shape(vector_values),1,1,1]) & + ,cells = mold%cells_, x_min = mold%x_min_, x_max = mold%x_max_, order = mold%order_ & + ) + end associate vector_2D%divergence_operator_1D_ = [(divergence_operator_1D_t(k=mold%order_, dx=((mold%x_max_(dir)-mold%x_min_(dir))/mold%cells_(dir)), cells=mold%cells_(dir)), dir=1,space_dimension)] end associate end procedure module procedure vector_2D_values - vector_values = self%values_ + vector_values = self%values_(:,:,:,1,1,1) end procedure end submodule vector_2D_s \ No newline at end of file From 662eec327668254910f443b4e908b755384cd540 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Sun, 17 May 2026 01:10:03 -0700 Subject: [PATCH 10/36] chore: reintroduce only clause in formal_m --- src/formal_m.f90 | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/formal_m.f90 b/src/formal_m.f90 index d03f6bd..c51649e 100644 --- a/src/formal_m.f90 +++ b/src/formal_m.f90 @@ -18,10 +18,12 @@ module formal_m ,d_dx & ! scalar_1D_t spatial derivative ,d2_dx2 ! scalar_1D_t spatial derivative - use tensors_2D_m !, only : & - ! scalar_2D_t & ! discrete 2D scalar field derived type - ! ,vector_2D_t & ! discrete 2D vector field derived type - ! ,gradient_2D_t & ! result of `.grad. s` for a scalar_2D_t s + use tensors_2D_m, only : & + scalar_2D_t & ! discrete 2D scalar field derived type + ,vector_2D_t & ! discrete 2D vector field derived type + ,gradient_2D_t & ! result of `.grad. s` for a scalar_2D_t s + ,scalar_2D_initializer_i & ! scalar_2D_t initializer abstract interface + ,vector_2D_initializer_i ! vector_2D_t initializar abstract interface use differential_operators_1D_m, only : & gradient_operator_1D_t & ! matrix operator defining a 1D mimetic gradient From c6b6a0788168b97f78b49e23acec18a5c31f96f1 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Sun, 17 May 2026 09:53:42 -0700 Subject: [PATCH 11/36] test(scalar_2D): more stringent test This commit tests the .grad. operator with a scalar_2D_t defined as the biquadratic function z = 1 - 2*x + 3*x**2 - x*y/5 + 3*y**2 - 2*y which has the gradient g = [-2 + 6*x - y/5, -x/5 + 6*y - 2] on the domain cells=[30,20], x_min=[-1D0,1D0], x_max=[9D0,4D0]. --- src/formal/scalar_2D_s.F90 | 4 ++-- test/scalar_2D_test_m.F90 | 16 ++++++++-------- 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/src/formal/scalar_2D_s.F90 b/src/formal/scalar_2D_s.F90 index 0e9b44a..4cd98f1 100644 --- a/src/formal/scalar_2D_s.F90 +++ b/src/formal/scalar_2D_s.F90 @@ -62,12 +62,12 @@ 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:self%cells_(2)+2) default(none) shared(gradient_2D, self) + 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:self%cells_(1)+2) default(none) shared(gradient_2D, self) + 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 diff --git a/test/scalar_2D_test_m.F90 b/test/scalar_2D_test_m.F90 index 45c1d4f..2be8f16 100644 --- a/test/scalar_2D_test_m.F90 +++ b/test/scalar_2D_test_m.F90 @@ -25,7 +25,7 @@ module scalar_2D_test_m procedure, nopass :: results end type - double precision, parameter :: tolerance = 2D-11 + double precision, parameter :: tolerance = 5D-2 integer, parameter :: space_dimension = 2 contains @@ -44,19 +44,19 @@ function results() result(test_results) ]) end function - pure function bilinear(x,y) result(z) + pure function biquadratic(x,y) result(z) double precision, intent(in) :: x(:), y(:) double precision z(size(x),size(y)) do concurrent(integer :: j=1:size(y)) default(none) shared(x,y,z) - z(:,j) = 2*x + 3*y(j) + z(:,j) = 1 - 2*x + 3*x**2 - x*y(j)/5 + 3*y(j)**2 - 2*y(j) end do end function - pure function bilinear_gradient(x,y) result(gradient) + pure function biquadratic_gradient(x,y) result(gradient) double precision, intent(in) :: x(:), y(:) double precision gradient(size(x),size(y),space_dimension) do concurrent(integer :: i=1:size(x), j=1:size(y)) - gradient(i,j,:) = [2,3] + gradient(i,j,:) = [-2 + 6*x(i) - y(j)/5, -x(i)/5 + 6*y(j) - 2] end do end function @@ -66,12 +66,12 @@ function check_gradient() result(test_diagnosis) procedure(vector_2D_initializer_i), pointer :: expected_gradient_initializer integer order - scalar_2D_initializer => bilinear - expected_gradient_initializer => bilinear_gradient + scalar_2D_initializer => biquadratic + expected_gradient_initializer => biquadratic_gradient test_diagnosis = passing_test() do order = 2, 4, 2 - associate(scalar_2D => scalar_2D_t(scalar_2D_initializer, order=order, cells=[10,10], x_min=[0D0,0D0], x_max=[10D0,10D0])) + associate(scalar_2D => scalar_2D_t(scalar_2D_initializer, order=order, cells=[30,20], x_min=[-1D0,1D0], x_max=[9D0,4D0])) associate(grad_scalar => .grad. scalar_2D, expected_gradient => vector_2D_t(expected_gradient_initializer, mold=scalar_2D)) test_diagnosis = test_diagnosis .also. & .all. (grad_scalar%values() .approximates. expected_gradient%values() .within. tolerance) & From 0d65be47a7fce8bcdae7e91b6c79b7021480c9bf Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Sun, 17 May 2026 12:19:05 -0700 Subject: [PATCH 12/36] feat(example): save & plot scalar surface This commit adds 1. A scalar_2D_t "to_file" type-bound procedure that creates a Juliennne file_t object containing points for a surface plot, 2. An scalar-surface example that creates a scalar and saves it to example/scripts/scalar-surface.csv, and 3. A scalar-surface.gnuplot script that plots the surface and saves it to scalar-surface.gif. --- example/scalar-surface.F90 | 50 ++++++++++++++++++++++++++ example/scripts/scalar-surface.gnuplot | 33 +++++++++++++++++ src/formal/scalar_2D_s.F90 | 30 ++++++++++++++-- src/formal/tensors_1D_m.F90 | 3 +- src/formal/tensors_2D_m.F90 | 16 +++++++-- src/formal/vector_1D_s.F90 | 4 +-- test/scalar_2D_test_m.F90 | 2 +- 7 files changed, 128 insertions(+), 10 deletions(-) create mode 100644 example/scalar-surface.F90 create mode 100644 example/scripts/scalar-surface.gnuplot diff --git a/example/scalar-surface.F90 b/example/scalar-surface.F90 new file mode 100644 index 0000000..8d5aab2 --- /dev/null +++ b/example/scalar-surface.F90 @@ -0,0 +1,50 @@ +! Copyright (c) 2026, The Regents of the University of California +! Terms of use are as specified in LICENSE.txt + +module scalar_2D_functions_m + implicit none + + integer, parameter :: space_dimension = 2 + +contains + + pure function biquadratic(x,y) result(z) + double precision, intent(in) :: x(:), y(:) + double precision z(size(x),size(y)) + do concurrent(integer :: j=1:size(y)) default(none) shared(x,y,z) + z(:,j) = 1 - 2*x + 3*x**2 - x*y(j)/5 + 3*y(j)**2 - 2*y(j) + end do + end function + + pure function biquadratic_gradient(x,y) result(gradient) + double precision, intent(in) :: x(:), y(:) + double precision gradient(size(x),size(y),space_dimension) + do concurrent(integer :: i=1:size(x), j=1:size(y)) + gradient(i,j,:) = [-2 + 6*x(i) - y(j)/5, -x(i)/5 + 6*y(j) - 2] + end do + end function + +end module scalar_2D_functions_m + +program scalar_surface + use julienne_m, only : file_t + use scalar_2D_functions_m, only : biquadratic, biquadratic_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 :: expected_gradient_initializer + integer, parameter :: order = 4 + + scalar_2D_initializer => biquadratic + expected_gradient_initializer => biquadratic_gradient + + associate(scalar_2D => scalar_2D_t(scalar_2D_initializer, order=order, cells=[30,20], x_min=[-1D0,1D0], x_max=[9D0,4D0])) + associate(grad_scalar => .grad. scalar_2D, expected_gradient => vector_2D_t(expected_gradient_initializer, mold=scalar_2D)) + associate(scalar_2D_file => scalar_2D%to_file()) + call scalar_2D_file%write_lines("example/scripts/scalar-surface.csv") + end associate + end associate + end associate + +end program scalar_surface \ No newline at end of file diff --git a/example/scripts/scalar-surface.gnuplot b/example/scripts/scalar-surface.gnuplot new file mode 100644 index 0000000..d7418c8 --- /dev/null +++ b/example/scripts/scalar-surface.gnuplot @@ -0,0 +1,33 @@ +# ============================================================ +# surface.gnuplot -- surface plot from a pre-gridded CSV +# Line 1: column labels +# Lines 2+: x, y, z data with blank lines between x-slices +# Usage: gnuplot scalar-surface.gnuplot +# ============================================================ + +datafile = "scalar-surface.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 "scalar-surface.gif" + +splot datafile every ::1 using 1:2:3 with pm3d title "" + +set output # flush and close the file diff --git a/src/formal/scalar_2D_s.F90 b/src/formal/scalar_2D_s.F90 index 4cd98f1..e6d23c6 100644 --- a/src/formal/scalar_2D_s.F90 +++ b/src/formal/scalar_2D_s.F90 @@ -11,6 +11,7 @@ ,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 @@ -50,7 +51,7 @@ 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 grad + module procedure scalar_2D_gradient integer c, i, j @@ -81,4 +82,29 @@ end procedure -end submodule scalar_2D_s + 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 \ No newline at end of file diff --git a/src/formal/tensors_1D_m.F90 b/src/formal/tensors_1D_m.F90 index b0f4dd6..7089a94 100644 --- a/src/formal/tensors_1D_m.F90 +++ b/src/formal/tensors_1D_m.F90 @@ -180,10 +180,9 @@ pure module function construct_1D_vector_from_function(initializer, order, cells type(vector_1D_t) vector_1D end function - pure module function construct_from_components(tensor_1D, divergence_operator_1D) result(vector_1D) + pure module function construct_1D_vector_from_parent(tensor_1D) result(vector_1D) !! Result is a 1D vector with the provided parent component tensor_1D and the provided divergence operatror type(tensor_1D_t), intent(in) :: tensor_1D - type(divergence_operator_1D_t), intent(in) :: divergence_operator_1D type(vector_1D_t) vector_1D end function diff --git a/src/formal/tensors_2D_m.F90 b/src/formal/tensors_2D_m.F90 index e0ed9de..8a2b0c0 100644 --- a/src/formal/tensors_2D_m.F90 +++ b/src/formal/tensors_2D_m.F90 @@ -6,6 +6,7 @@ module tensors_2D_m !! divergence, and Laplacian operators as detailed by Corbino & Castillo (2020) !! https://doi.org/10.1016/j.cam.2019.06.042. use differential_operators_1D_m, only : gradient_operator_1D_t, divergence_operator_1D_t + use julienne_m, only : file_t implicit none @@ -69,10 +70,12 @@ pure module function construct_2D_tensor_from_components(values, cells, x_min, x private type(gradient_operator_1D_t) gradient_operator_1D_(space_dimension) contains - generic :: operator(.grad.) => grad + generic :: operator(.grad.) => scalar_2D_gradient generic :: values => scalar_2D_values generic :: grid => scalar_2D_grid - procedure, non_overridable, private :: grad + generic :: to_file => scalar_2D_to_file + procedure, non_overridable, private :: scalar_2D_to_file + procedure, non_overridable, private :: scalar_2D_gradient procedure, non_overridable, private :: scalar_2D_values procedure, non_overridable, private :: scalar_2D_grid end type @@ -168,13 +171,20 @@ pure module function vector_2D_values(self) result(vector_values) double precision, allocatable :: vector_values(:,:,:) end function - pure module function grad(self) result(gradient_2D) + pure module function scalar_2D_gradient(self) result(gradient_2D) !! Result is mimetic gradient of the scalar_2D_t "self" implicit none class(scalar_2D_t), intent(in) :: self type(gradient_2D_t) gradient_2D end function + pure module function scalar_2D_to_file(self) result(file) + !! Result is a file_t object containing the grid points and the corresponding scalar values + implicit none + class(scalar_2D_t), intent(in) :: self + type(file_t) file + end function + end interface end module tensors_2D_m diff --git a/src/formal/vector_1D_s.F90 b/src/formal/vector_1D_s.F90 index 1be73da..3bd8f00 100644 --- a/src/formal/vector_1D_s.F90 +++ b/src/formal/vector_1D_s.F90 @@ -60,9 +60,9 @@ pure module function construct_1D_vector_from_function(initializer, order, cells #endif - module procedure construct_from_components + module procedure construct_1D_vector_from_parent vector_1D%tensor_1D_t = tensor_1D - vector_1D%divergence_operator_1D_ = divergence_operator_1D + vector_1D%divergence_operator_1D_ = divergence_operator_1D_t(k=tensor_1D%order_, dx=(tensor_1D%x_max_ - tensor_1D%x_min_)/tensor_1D%cells_, cells=tensor_1D%cells_) end procedure module procedure div diff --git a/test/scalar_2D_test_m.F90 b/test/scalar_2D_test_m.F90 index 2be8f16..d1a470f 100644 --- a/test/scalar_2D_test_m.F90 +++ b/test/scalar_2D_test_m.F90 @@ -55,7 +55,7 @@ pure function biquadratic(x,y) result(z) pure function biquadratic_gradient(x,y) result(gradient) double precision, intent(in) :: x(:), y(:) double precision gradient(size(x),size(y),space_dimension) - do concurrent(integer :: i=1:size(x), j=1:size(y)) + do concurrent(integer :: i=1:size(x), j=1:size(y)) default(none) shared(gradient,x,y) gradient(i,j,:) = [-2 + 6*x(i) - y(j)/5, -x(i)/5 + 6*y(j) - 2] end do end function From c4f42f9a5e9ce11fddf263e394c16033205bb7e0 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Sun, 17 May 2026 17:25:02 -0700 Subject: [PATCH 13/36] feat: vector_1D const construct, vector_2D grid This commit adds 1. A vector_1D constructor for constant vector fields and 2. A vector_2D grid calculator. --- src/formal/tensors_1D_m.F90 | 11 ++++++++++ src/formal/tensors_2D_m.F90 | 18 ++++++++++++++++ src/formal/vector_1D_s.F90 | 16 ++++++++++++++ src/formal/vector_2D_s.F90 | 43 +++++++++++++++++++++++++++++++++++-- 4 files changed, 86 insertions(+), 2 deletions(-) diff --git a/src/formal/tensors_1D_m.F90 b/src/formal/tensors_1D_m.F90 index 7089a94..4555c6f 100644 --- a/src/formal/tensors_1D_m.F90 +++ b/src/formal/tensors_1D_m.F90 @@ -186,6 +186,17 @@ pure module function construct_1D_vector_from_parent(tensor_1D) result(vector_1D type(vector_1D_t) vector_1D end function + pure module function construct_1D_vector_constant(constant, order, cells, x_min, x_max) result(vector_1D) + !! Result is a collection of cell-centered-extended values with a corresponding mimetic gradient operator + implicit none + double precision, intent(in) :: constant !! scalar value + integer, intent(in) :: order !! order of accuracy + integer, intent(in) :: cells !! number of grid cells spanning the domain + double precision, intent(in) :: x_min !! grid location minimum + double precision, intent(in) :: x_max !! grid location maximum + type(vector_1D_t) vector_1D + end function + end interface type, extends(vector_1D_t) :: gradient_1D_t diff --git a/src/formal/tensors_2D_m.F90 b/src/formal/tensors_2D_m.F90 index 8a2b0c0..d298966 100644 --- a/src/formal/tensors_2D_m.F90 +++ b/src/formal/tensors_2D_m.F90 @@ -109,7 +109,11 @@ pure module function construct_2D_scalar_from_mold(initializer, mold) result(sca type(divergence_operator_1D_t) divergence_operator_1D_(space_dimension) contains generic :: values => vector_2D_values + generic :: to_file => vector_2D_to_file + generic :: grid => vector_2D_grid procedure, non_overridable, private :: vector_2D_values + procedure, non_overridable, private :: vector_2D_to_file + procedure, non_overridable, private :: vector_2D_grid end type interface vector_2D_t @@ -165,6 +169,13 @@ pure module function scalar_2D_grid(self, direction) result(scalar_grid_1D) double precision, allocatable :: scalar_grid_1D(:) end function + pure module function vector_2D_grid(self, direction) result(vector_grid_1D) + !! Result array contains scalar grid locations along the requested spatial direction + class(vector_2D_t), intent(in) :: self + integer, intent(in) :: direction + double precision, allocatable :: vector_grid_1D(:) !! grid points along one the requested coordinate direction + end function + pure module function vector_2D_values(self) result(vector_values) !! Vector values getter class(vector_2D_t), intent(in) :: self @@ -185,6 +196,13 @@ pure module function scalar_2D_to_file(self) result(file) type(file_t) file end function + pure module function vector_2D_to_file(self) result(file) + !! Result is a file_t object containing the grid points and the corresponding vector components + implicit none + class(vector_2D_t), intent(in) :: self + type(file_t) file + end function + end interface end module tensors_2D_m diff --git a/src/formal/vector_1D_s.F90 b/src/formal/vector_1D_s.F90 index 3bd8f00..1467d38 100644 --- a/src/formal/vector_1D_s.F90 +++ b/src/formal/vector_1D_s.F90 @@ -65,6 +65,22 @@ pure module function construct_1D_vector_from_function(initializer, order, cells vector_1D%divergence_operator_1D_ = divergence_operator_1D_t(k=tensor_1D%order_, dx=(tensor_1D%x_max_ - tensor_1D%x_min_)/tensor_1D%cells_, cells=tensor_1D%cells_) end procedure + module procedure construct_1D_vector_constant + + integer i + + call_julienne_assert(x_max .greaterThan. x_min) + call_julienne_assert(cells .isAtLeast. 2*order) + + vector_1D = vector_1D_t( tensor_1D_t( & + values = [(constant, i = 1, size(faces_1D(x_min, x_max, cells)))] & + ,x_min = x_min & + ,x_max = x_max & + ,cells = cells & + ,order = order & + ) ) + end procedure + module procedure div integer center diff --git a/src/formal/vector_2D_s.F90 b/src/formal/vector_2D_s.F90 index 09a5e5b..b04f3c5 100644 --- a/src/formal/vector_2D_s.F90 +++ b/src/formal/vector_2D_s.F90 @@ -7,10 +7,12 @@ use julienne_m, only : & call_julienne_assert_ & ,operator(.all.) & + ,operator(.csv.) & ,operator(.equalsExpected.) & ,operator(.greaterThan.) & - ,operator(.isAtLeast.) - use tensors_1D_m, only : faces_1D + ,operator(.isAtLeast.) & + ,string_t + use tensors_1D_m, only : faces_1D, vector_1D_t use differential_operators_1D_m, only : divergence_operator_1D_t implicit none @@ -75,4 +77,41 @@ vector_values = self%values_(:,:,:,1,1,1) end procedure + module procedure vector_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,vector_x,vector_y")]) + 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:space_dimension,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 + + module procedure vector_2D_grid + associate(vector_1D => vector_1D_t( & + constant = 0D0 & + ,cells = self%cells_(direction) & + ,x_min = self%x_min_(direction) & + ,x_max = self%x_max_(direction) & + ,order = self%order_ & + )) + vector_grid_1D = vector_1D%grid() + end associate + end procedure + end submodule vector_2D_s \ No newline at end of file From 00efe687503e8661af0034487bcb080f617b13c0 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Sun, 17 May 2026 19:52:54 -0700 Subject: [PATCH 14/36] feat(3D): scalar, vector, & gradient This commit 1. Adds public scalar_3D_t and vector_3D_t types and a supporting private tensors_3D_t type and 2. A passing unit test for the gradient of a 3D scalar field. --- src/formal/scalar_3D_s.F90 | 126 ++++++++++++++++++++++ src/formal/tensor_3D_s.F90 | 14 +++ src/formal/tensors_2D_m.F90 | 4 +- src/formal/tensors_3D_m.F90 | 208 ++++++++++++++++++++++++++++++++++++ src/formal/vector_3D_s.F90 | 134 +++++++++++++++++++++++ src/formal_m.f90 | 7 ++ test/driver.f90 | 2 + test/scalar_3D_test_m.F90 | 84 +++++++++++++++ 8 files changed, 577 insertions(+), 2 deletions(-) create mode 100644 src/formal/scalar_3D_s.F90 create mode 100644 src/formal/tensor_3D_s.F90 create mode 100644 src/formal/tensors_3D_m.F90 create mode 100644 src/formal/vector_3D_s.F90 create mode 100644 test/scalar_3D_test_m.F90 diff --git a/src/formal/scalar_3D_s.F90 b/src/formal/scalar_3D_s.F90 new file mode 100644 index 0000000..28e4fc8 --- /dev/null +++ b/src/formal/scalar_3D_s.F90 @@ -0,0 +1,126 @@ +! 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_3D_m) scalar_3D_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_3D_values + scalar_values = self%values_(:,:,:,1,1,1,1) + end procedure + + module procedure scalar_3D_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_3D_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)) & + ,z => cell_centers_extended_1D(x_min(3), x_max(3), cells(3)) & + ) + scalar_3D%tensor_3D_t = tensor_3D_t( & + values = reshape(initializer(x,y,z), shape=[size(x),size(y),size(z),1,1,1,1]) & + ,cells = cells , x_min = x_min, x_max = x_max, order = order & + ) + scalar_3D%gradient_operator_1D_ = gradient_operator_1D_t(k=order, dx=(x_max - x_min)/cells, cells=cells) + end associate + end procedure + + module procedure construct_3D_scalar_from_mold + scalar_3D = scalar_3D_t(initializer, cells = mold%cells_, x_min = mold%x_min_, x_max = mold%x_max_, order = mold%order_) + end procedure + + module procedure scalar_3D_gradient + + integer c, i, j + + gradient_3D%x_min_ = self%x_min_ + gradient_3D%x_max_ = self%x_max_ + gradient_3D%cells_ = self%cells_ + gradient_3D%order_ = self%order_ + + allocate(gradient_3D%values_(self%cells_(1)+1, self%cells_(2)+1, self%cells_(3)+1, space_dimension, 1, 1, 1)) + + gradient_x_component: & + do concurrent(integer :: j=1:size(gradient_3D%values_,2), k=1:size(gradient_3D%values_,3)) default(none) shared(gradient_3D, self) + gradient_3D%values_(:,j,k,1,1,1,1) = self%gradient_operator_1D_(1) .x. self%values_(:,j,k,1,1,1,1) + end do gradient_x_component + + gradient_y_component: & + do concurrent(integer :: i=1:size(gradient_3D%values_,1), k=1:size(gradient_3D%values_,3)) default(none) shared(gradient_3D, self) + gradient_3D%values_(i,:,k,2,1,1,1) = self%gradient_operator_1D_(2) .x. self%values_(i,:,k,1,1,1,1) + end do gradient_y_component + + gradient_z_component: & + do concurrent(integer :: i=1:size(gradient_3D%values_,1), j=1:size(gradient_3D%values_,2)) default(none) shared(gradient_3D, self) + gradient_3D%values_(i,j,:,3,1,1,1) = self%gradient_operator_1D_(3) .x. self%values_(i,j,:,1,1,1,1) + end do gradient_z_component + + associate(dx => (self%x_max_ - self%x_min_)/self%cells_) + gradient_3D%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. 3D-3))) + !end associate check_corbino_castillo_eq_17 + end associate + + end procedure + + module procedure scalar_3D_to_file + type(string_t), allocatable :: lines(:) + integer i, j, k, l + + associate( & + x => self%grid(1) & + ,y => self%grid(2) & + ,z => self%grid(3) & + ,header => [string_t("x,y,z,scalar")] & + ) + associate(num_blank_lines => size(y)*size(z) - 1) + allocate(lines(size(header) + size(self%values_) + num_blank_lines)) + end associate + lines(1:size(header)) = header + l = size(header) + do k = 1, size(z) + do j = 1, size(y) + do i = 1, size(x) + l = l + 1 + lines(l) = .csv. string_t([x(i), y(j), z(k), self%values_(i,j,k,1,1,1,1)]) + end do + if (j/=size(y) .or. k/=size(z)) then + l = l + 1 + lines(l) = "" + end if + end do + end do + end associate + + file = file_t(lines) + end procedure + +end submodule scalar_3D_s \ No newline at end of file diff --git a/src/formal/tensor_3D_s.F90 b/src/formal/tensor_3D_s.F90 new file mode 100644 index 0000000..e26a1c5 --- /dev/null +++ b/src/formal/tensor_3D_s.F90 @@ -0,0 +1,14 @@ +submodule(tensors_3D_m) tensor_3D_s + implicit none + +contains + + module procedure construct_3D_tensor_from_components + tensor_3D%values_ = values + tensor_3D%cells_ = cells + tensor_3D%x_min_ = x_min + tensor_3D%x_max_ = x_max + tensor_3D%order_ = order + end procedure + +end submodule \ No newline at end of file diff --git a/src/formal/tensors_2D_m.F90 b/src/formal/tensors_2D_m.F90 index d298966..3f1a83e 100644 --- a/src/formal/tensors_2D_m.F90 +++ b/src/formal/tensors_2D_m.F90 @@ -44,7 +44,7 @@ pure function vector_2D_initializer_i(x,y) result(v) !! Child types define the operations supported by each child, including !! gradient (.grad.) for scalars and divergence (.div.) for vectors. private - double precision, allocatable :: values_(:,:,:,:,:,:) !! tensor components for rank<=4 at 2D locations + double precision, allocatable :: values_(:,:, :,:,:,:) !! tensor components for rank<=4 at 2D locations double precision x_min_(space_dimension) !! domain lower boundary double precision x_max_(space_dimension) !! domain upper boundary integer cells_(space_dimension) !! number of grid cells spanning the domain @@ -55,7 +55,7 @@ pure function vector_2D_initializer_i(x,y) result(v) pure module function construct_2D_tensor_from_components(values, cells, x_min, x_max, order) result(tensor_2D) implicit none - double precision, intent(in) :: values(:,:,:,:,:,:) !! tensor components at 2D spatial locations + double precision, intent(in) :: values(:,:, :,:,:,:) !! tensor components at 2D spatial locations double precision, intent(in) :: x_min(:) !! domain lower boundary double precision, intent(in) :: x_max(:) !! domain upper boundary integer, intent(in) :: cells(:) !! number of grid cells spanning the domain diff --git a/src/formal/tensors_3D_m.F90 b/src/formal/tensors_3D_m.F90 new file mode 100644 index 0000000..7b919a5 --- /dev/null +++ b/src/formal/tensors_3D_m.F90 @@ -0,0 +1,208 @@ +! Copyright (c) 2026, The Regents of the University of California +! Terms of use are as specified in LICENSE.txt + +module tensors_3D_m + !! Define public 3D scalar and vector abstractions and associated mimetic gradient, + !! divergence, and Laplacian operators as detailed by Corbino & Castillo (2020) + !! https://doi.org/10.1016/j.cam.2019.06.042. + use differential_operators_1D_m, only : gradient_operator_1D_t, divergence_operator_1D_t + use julienne_m, only : file_t + + implicit none + + private + + public :: scalar_3D_t + public :: vector_3D_t + public :: gradient_3D_t + public :: scalar_3D_initializer_i + public :: vector_3D_initializer_i + + integer, parameter :: space_dimension = 3 + + abstract interface + + pure function scalar_3D_initializer_i(x,y,z) result(f) + !! Sampling function for initializing a scalar_3D_t object + implicit none + double precision, intent(in) :: x(:), y(:), z(:) + double precision f(size(x),size(y),size(z)) + end function + + pure function vector_3D_initializer_i(x,y,z) result(v) + !! Sampling function for initializing a vector_3D_t object + import space_dimension + implicit none + double precision, intent(in) :: x(:), y(:), z(:) + double precision v(size(x),size(y),size(z),space_dimension) + end function + + end interface + + type tensor_3D_t + !! Encapsulate the components that are common to all 3D tensors. + !! Child types define the operations supported by each child, including + !! gradient (.grad.) for scalars and divergence (.div.) for vectors. + private + double precision, allocatable :: values_(:,:,:, :,:,:,:) !! tensor components for rank<=4 at 3D locations + double precision x_min_(space_dimension) !! domain lower boundary + double precision x_max_(space_dimension) !! domain upper boundary + integer cells_(space_dimension) !! number of grid cells spanning the domain + integer order_ !! order of accuracy of mimetic discretization + end type + + interface tensor_3D_t + + pure module function construct_3D_tensor_from_components(values, cells, x_min, x_max, order) result(tensor_3D) + implicit none + double precision, intent(in) :: values(:,:,:, :,:,:,:) !! tensor components for rank<=4 at 3D locations + double precision, intent(in) :: x_min(:) !! domain lower boundary + double precision, intent(in) :: x_max(:) !! domain upper boundary + integer, intent(in) :: cells(:) !! number of grid cells spanning the domain + integer, intent(in) :: order !! order of accuracy of mimetic discretization + type(tensor_3D_t) tensor_3D + end function + + end interface + + type, extends(tensor_3D_t) :: scalar_3D_t + !! Encapsulate scalar values at cell centers and boundaries + private + type(gradient_operator_1D_t) gradient_operator_1D_(space_dimension) + contains + generic :: operator(.grad.) => scalar_3D_gradient + generic :: values => scalar_3D_values + generic :: grid => scalar_3D_grid + generic :: to_file => scalar_3D_to_file + procedure, non_overridable, private :: scalar_3D_to_file + procedure, non_overridable, private :: scalar_3D_gradient + procedure, non_overridable, private :: scalar_3D_values + procedure, non_overridable, private :: scalar_3D_grid + end type + + interface scalar_3D_t + + pure module function construct_3D_scalar_from_function(initializer, order, cells, x_min, x_max) result(scalar_3D) + !! Result is a collection of cell-centered-extended values with a corresponding mimetic gradient operator + implicit none + procedure(scalar_3D_initializer_i), pointer :: initializer + integer, intent(in) :: order !! order of accuracy + integer, intent(in) :: cells(:) !! number of grid cells spanning each spatial direction + double precision, intent(in) :: x_min(:) !! grid location minima + double precision, intent(in) :: x_max(:) !! grid location maxima + type(scalar_3D_t) scalar_3D + end function + + pure module function construct_3D_scalar_from_mold(initializer, mold) result(scalar_3D) + !! Result is a 3D scalar field using a mold for all components other than the field values + implicit none + procedure(scalar_3D_initializer_i), pointer :: initializer + type(scalar_3D_t), intent(in) :: mold + type(scalar_3D_t) scalar_3D + end function + + end interface + + type, extends(tensor_3D_t) :: vector_3D_t + !! Encapsulate 3D vector values at cell faces (of unit area for 3D) and corresponding operators + private + type(divergence_operator_1D_t) divergence_operator_1D_(space_dimension) + contains + generic :: values => vector_3D_values + generic :: to_file => vector_3D_to_file + generic :: grid => vector_3D_grid + procedure, non_overridable, private :: vector_3D_values + procedure, non_overridable, private :: vector_3D_to_file + procedure, non_overridable, private :: vector_3D_grid + end type + + interface vector_3D_t + + pure module function construct_3D_vector_from_function(initializer, order, cells, x_min, x_max) result(vector_3D) + !! Result is a 3D vector with values initialized by the provided procedure pointer sampled on the faces of + !! the specified number of evenly spaced cells covering [x_min, x_max] + implicit none + procedure(vector_3D_initializer_i), pointer :: initializer + integer, intent(in) :: order !! order of accuracy + integer, intent(in) :: cells(:) !! number of grid cells spanning each spatial direction + double precision, intent(in) :: x_min(:) !! grid location minima + double precision, intent(in) :: x_max(:) !! grid location maxima + type(vector_3D_t) vector_3D + end function + + pure module function construct_3D_vector_from_vector_mold(initializer, mold) result(vector_3D) + !! Result is a 3D vector with values initialized by the provided procedure pointer sampled on the + !! same grid as the mold + implicit none + procedure(vector_3D_initializer_i), pointer :: initializer + type(vector_3D_t), intent(in) :: mold + type(vector_3D_t) vector_3D + end function + + pure module function construct_3D_vector_from_scalar_mold(initializer, mold) result(vector_3D) + !! Result is a 3D vector with values initialized by the provided procedure pointer sampled on the + !! face-centered grid corresponding to the cell-centered grid of the mold + implicit none + procedure(vector_3D_initializer_i), pointer :: initializer + type(scalar_3D_t), intent(in) :: mold + type(vector_3D_t) vector_3D + end function + + end interface + + type, extends(vector_3D_t) :: gradient_3D_t + !! A 3D mimetic gradient vector field abstraction with a public method that produces corresponding numerical quadrature weights + end type + + interface + + pure module function scalar_3D_values(self) result(scalar_values) + !! Scalar values getter + class(scalar_3D_t), intent(in) :: self + double precision, allocatable :: scalar_values(:,:,:) + end function + + pure module function scalar_3D_grid(self, direction) result(scalar_grid_1D) + !! Result contains scalar grid locations along the requested spatial direction + class(scalar_3D_t), intent(in) :: self + integer, intent(in) :: direction + double precision, allocatable :: scalar_grid_1D(:) + end function + + pure module function vector_3D_grid(self, direction) result(vector_grid_1D) + !! Result contains scalar grid locations along the requested spatial direction + class(vector_3D_t), intent(in) :: self + integer, intent(in) :: direction + double precision, allocatable :: vector_grid_1D(:) !! grid points along one the requested coordinate direction + end function + + pure module function vector_3D_values(self) result(vector_values) + !! Vector values getter + class(vector_3D_t), intent(in) :: self + double precision, allocatable :: vector_values(:,:,:,:) + end function + + pure module function scalar_3D_gradient(self) result(gradient_3D) + !! Result is the mimetic gradient of the scalar_3D_t "self" + implicit none + class(scalar_3D_t), intent(in) :: self + type(gradient_3D_t) gradient_3D + end function + + pure module function scalar_3D_to_file(self) result(file) + !! Result is a file_t object containing the grid points and the corresponding scalar values + implicit none + class(scalar_3D_t), intent(in) :: self + type(file_t) file + end function + + pure module function vector_3D_to_file(self) result(file) + !! Result is a file_t object containing the grid points and the corresponding vector components + implicit none + class(vector_3D_t), intent(in) :: self + type(file_t) file + end function + + end interface + +end module tensors_3D_m diff --git a/src/formal/vector_3D_s.F90 b/src/formal/vector_3D_s.F90 new file mode 100644 index 0000000..8191bbe --- /dev/null +++ b/src/formal/vector_3D_s.F90 @@ -0,0 +1,134 @@ +! 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_3D_m) vector_3D_s + use julienne_m, only : & + call_julienne_assert_ & + ,operator(.all.) & + ,operator(.csv.) & + ,operator(.equalsExpected.) & + ,operator(.greaterThan.) & + ,operator(.isAtLeast.) & + ,string_t + use tensors_1D_m, only : faces_1D, vector_1D_t + use differential_operators_1D_m, only : divergence_operator_1D_t + implicit none + +contains + + module procedure construct_3D_vector_from_function + + integer dir + + 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 => faces_1D(x_min(1), x_max(1), cells(1)) & + ,y => faces_1D(x_min(2), x_max(2), cells(2)) & + ,z => faces_1D(x_min(3), x_max(3), cells(3)) & + ) + associate(vector_values => initializer(x,y,z)) + vector_3D%tensor_3D_t = tensor_3D_t( & + values = reshape(vector_values, shape=[shape(vector_values),1,1,1]) & + ,cells = cells, x_min = x_min, x_max = x_max, order = order & + ) + end associate + vector_3D%divergence_operator_1D_ = & + [(divergence_operator_1D_t(k=order, dx=((x_max(dir)-x_min(dir))/cells(dir)), cells=cells(dir)), dir=1,space_dimension)] + end associate + end procedure + + module procedure construct_3D_vector_from_vector_mold + integer dir + + call_julienne_assert(.all. ([size(mold%cells_), size(mold%x_min_), size(mold%x_max_)] .equalsExpected. space_dimension)) + call_julienne_assert(.all. (mold%x_max_ .greaterThan. mold%x_min_)) + call_julienne_assert(.all. (mold%cells_ .isAtLeast. 2*mold%order_)) + + associate( & + x => faces_1D(mold%x_min_(1), mold%x_max_(1), mold%cells_(1)) & + ,y => faces_1D(mold%x_min_(2), mold%x_max_(2), mold%cells_(2)) & + ,z => faces_1D(mold%x_min_(3), mold%x_max_(3), mold%cells_(3)) & + ) + associate(vector_values => initializer(x,y,z)) + vector_3D%tensor_3D_t = tensor_3D_t( & + values = reshape(vector_values, shape=[shape(vector_values),1,1,1]) & + ,cells = mold%cells_, x_min = mold%x_min_, x_max = mold%x_max_, order = mold%order_ & + ) + end associate + vector_3D%divergence_operator_1D_ = & + [(divergence_operator_1D_t(k=mold%order_, dx=((mold%x_max_(dir)-mold%x_min_(dir))/mold%cells_(dir)), cells=mold%cells_(dir)), dir=1,space_dimension)] + end associate + end procedure + + module procedure construct_3D_vector_from_scalar_mold + integer dir + + call_julienne_assert(.all. ([size(mold%cells_), size(mold%x_min_), size(mold%x_max_)] .equalsExpected. space_dimension)) + call_julienne_assert(.all. (mold%x_max_ .greaterThan. mold%x_min_)) + call_julienne_assert(.all. (mold%cells_ .isAtLeast. 2*mold%order_)) + + associate( & + x => faces_1D(mold%x_min_(1), mold%x_max_(1), mold%cells_(1)) & + ,y => faces_1D(mold%x_min_(2), mold%x_max_(2), mold%cells_(2)) & + ,z => faces_1D(mold%x_min_(3), mold%x_max_(3), mold%cells_(3)) & + ) + associate(vector_values => initializer(x,y,z)) + vector_3D%tensor_3D_t = tensor_3D_t( & + values = reshape(vector_values, shape=[shape(vector_values),1,1,1]) & + ,cells = mold%cells_, x_min = mold%x_min_, x_max = mold%x_max_, order = mold%order_ & + ) + end associate + vector_3D%divergence_operator_1D_ = & + [(divergence_operator_1D_t(k=mold%order_, dx=((mold%x_max_(dir)-mold%x_min_(dir))/mold%cells_(dir)), cells=mold%cells_(dir)), dir=1,space_dimension)] + end associate + end procedure + + module procedure vector_3D_values + vector_values = self%values_(:,:,:,:,1,1,1) + end procedure + + module procedure vector_3D_to_file + type(string_t), allocatable :: lines(:) + integer i, j, k, l + + associate(x => self%grid(1), y => self%grid(2), z => self%grid(3), header => [string_t("x,y,vector_x,vector_y")]) + associate(num_blank_lines => size(y)*size(z) - 1) + allocate(lines(size(header) + size(self%values_) + num_blank_lines)) + end associate + lines(1:size(header)) = header + l = size(header) + do k = 1, size(z) + do j = 1, size(y) + do i = 1, size(x) + l = l + 1 + lines(l) = .csv. string_t([x(i), y(j), z(k), self%values_(i,j,k,1:space_dimension,1,1,1)]) + end do + if (j/=size(y) .or. k/=size(z)) then + l = l + 1 + lines(l) = "" + end if + end do + end do + end associate + + file = file_t(lines) + end procedure + + module procedure vector_3D_grid + associate(vector_1D => vector_1D_t( & + constant = 0D0 & + ,cells = self%cells_(direction) & + ,x_min = self%x_min_(direction) & + ,x_max = self%x_max_(direction) & + ,order = self%order_ & + )) + vector_grid_1D = vector_1D%grid() + end associate + end procedure + +end submodule vector_3D_s \ No newline at end of file diff --git a/src/formal_m.f90 b/src/formal_m.f90 index c51649e..86d4319 100644 --- a/src/formal_m.f90 +++ b/src/formal_m.f90 @@ -25,6 +25,13 @@ module formal_m ,scalar_2D_initializer_i & ! scalar_2D_t initializer abstract interface ,vector_2D_initializer_i ! vector_2D_t initializar abstract interface + use tensors_3D_m, only : & + scalar_3D_t & ! discrete 3D scalar field derived type + ,vector_3D_t & ! discrete 3D vector field derived type + ,gradient_3D_t & ! result of `.grad. s` for a scalar_3D_t s + ,scalar_3D_initializer_i & ! scalar_3D_t initializer abstract interface + ,vector_3D_initializer_i ! vector_3D_t initializar abstract interface + use differential_operators_1D_m, only : & gradient_operator_1D_t & ! matrix operator defining a 1D mimetic gradient ,divergence_operator_1D_t ! matrix operator defining a 1D mimetic divergence diff --git a/test/driver.f90 b/test/driver.f90 index 1464cc7..bffe1e5 100644 --- a/test/driver.f90 +++ b/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 scalar_3D_test_m, only : scalar_3D_test_t implicit none associate(test_harness => test_harness_t([ & @@ -21,6 +22,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(scalar_3D_test_t()) & ])) call test_harness%report_results end associate diff --git a/test/scalar_3D_test_m.F90 b/test/scalar_3D_test_m.F90 new file mode 100644 index 0000000..6252230 --- /dev/null +++ b/test/scalar_3D_test_m.F90 @@ -0,0 +1,84 @@ +! Copyright (c) 2026, The Regents of the University of California +! Terms of use are as specified in LICENSE.txt + +module scalar_3D_test_m + use julienne_m, only : & + operator(//) & + ,operator(.all.) & + ,operator(.also.) & + ,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_3D_t, vector_3D_t, scalar_3D_initializer_i, vector_3D_initializer_i + + implicit none + + type, extends(test_t) :: scalar_3D_test_t + contains + procedure, nopass :: subject + procedure, nopass :: results + end type + + double precision, parameter :: tolerance = 1D-2 + integer, parameter :: space_dimension = 3 + +contains + + pure function subject() result(test_subject) + character(len=:), allocatable :: test_subject + test_subject = 'The scalar_3D_t derived type' + end function + + function results() result(test_results) + type(scalar_3D_test_t) scalar_3D_test + type(test_result_t), allocatable :: test_results(:) + + test_results = scalar_3D_test%run([ & + test_description_t('computing the gradient of a scalar field', usher(check_gradient)) & + ]) + end function + + pure function triquadratic(x,y,z) result(f) + double precision, intent(in) :: x(:), y(:), z(:) + double precision f(size(x),size(y),size(z)) + do concurrent(integer :: j=1:size(y), k=1:size(z)) default(none) shared(x,y,z,f) + f(:,j,k) = 1 - 2*x + 3*x**2 - x*y(j)/5 + 3*y(j)**2 - 2*y(j) - 2*z(k) + end do + end function + + pure function triquadratic_gradient(x,y,z) result(gradient) + double precision, intent(in) :: x(:), y(:), z(:) + double precision gradient(size(x),size(y),size(z),space_dimension) + do concurrent(integer :: i=1:size(x), j=1:size(y), k=1:size(z)) default(none) shared(gradient,x,y,z) + gradient(i,j,k,:) = [-2 + 6*x(i) - y(j)/5, -x(i)/5 + 6*y(j) - 2, -2D0] + end do + end function + + function check_gradient() result(test_diagnosis) + type(test_diagnosis_t) test_diagnosis + procedure(scalar_3D_initializer_i), pointer :: scalar_3D_initializer + procedure(vector_3D_initializer_i), pointer :: expected_gradient_initializer + integer order + + scalar_3D_initializer => triquadratic + expected_gradient_initializer => triquadratic_gradient + test_diagnosis = passing_test() + + do order = 2, 4, 2 + associate(scalar_3D => scalar_3D_t(scalar_3D_initializer, order=order, cells=[30,20,10], x_min=[0D0,0D0,0D0], x_max=[1D0,1D0,1D0])) + associate(grad_scalar => .grad. scalar_3D, expected_gradient => vector_3D_t(expected_gradient_initializer, mold=scalar_3D)) + test_diagnosis = test_diagnosis .also. & + .all. (grad_scalar%values() .approximates. expected_gradient%values() .within. tolerance) & + // string_t(" for order ") // string_t(order) + end associate + end associate + end do + end function + +end module scalar_3D_test_m \ No newline at end of file From e86f82a7745ecfc8110bd684f8b53b124c3827d4 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Sun, 17 May 2026 20:07:35 -0700 Subject: [PATCH 15/36] fix(vector_2D%to_file()): trim trailing blanks --- example/{scalar-surface.F90 => plot-2D-scalar-gradient.F90} | 3 ++- src/formal/vector_2D_s.F90 | 2 +- 2 files changed, 3 insertions(+), 2 deletions(-) rename example/{scalar-surface.F90 => plot-2D-scalar-gradient.F90} (90%) diff --git a/example/scalar-surface.F90 b/example/plot-2D-scalar-gradient.F90 similarity index 90% rename from example/scalar-surface.F90 rename to example/plot-2D-scalar-gradient.F90 index 8d5aab2..f23541b 100644 --- a/example/scalar-surface.F90 +++ b/example/plot-2D-scalar-gradient.F90 @@ -41,8 +41,9 @@ program scalar_surface associate(scalar_2D => scalar_2D_t(scalar_2D_initializer, order=order, cells=[30,20], x_min=[-1D0,1D0], x_max=[9D0,4D0])) associate(grad_scalar => .grad. scalar_2D, expected_gradient => vector_2D_t(expected_gradient_initializer, mold=scalar_2D)) - associate(scalar_2D_file => scalar_2D%to_file()) + associate(scalar_2D_file => scalar_2D%to_file(), grad_scalar_file => grad_scalar%to_file()) call scalar_2D_file%write_lines("example/scripts/scalar-surface.csv") + call grad_scalar_file%write_lines("example/scripts/gradient-field.csv") end associate end associate end associate diff --git a/src/formal/vector_2D_s.F90 b/src/formal/vector_2D_s.F90 index b04f3c5..73f9554 100644 --- a/src/formal/vector_2D_s.F90 +++ b/src/formal/vector_2D_s.F90 @@ -83,7 +83,7 @@ associate(x => self%grid(1), y => self%grid(2), header => [string_t("x,y,vector_x,vector_y")]) associate(num_blank_lines => size(y)-1) - allocate(lines(size(header) + size(self%values_) + num_blank_lines)) + allocate(lines(size(header) + size(self%values_)/space_dimension + num_blank_lines)) end associate lines(1:size(header)) = header l = size(header) From 842fabe39d5f74279a8250c7962636b5b38451de Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Sun, 17 May 2026 23:13:48 -0700 Subject: [PATCH 16/36] feat(example): refactor/enhance scalar/vector plot This commit 1. Updates the scalar-surface plot example to also output files that can be used to plot the gradient of the surface, 2. Redefines the surface so that it corresponds to a velocity potential defining an irrotational vortex, 3. Adds a script that plots the resulting velocity vield and the expected velocity field. --- example/2D-vortex.F90 | 58 +++++++++++++++++++ example/plot-2D-scalar-gradient.F90 | 51 ---------------- example/scripts/velocities.gnuplot | 34 +++++++++++ ...ace.gnuplot => velocity-potential.gnuplot} | 9 +-- 4 files changed, 97 insertions(+), 55 deletions(-) create mode 100644 example/2D-vortex.F90 delete mode 100644 example/plot-2D-scalar-gradient.F90 create mode 100644 example/scripts/velocities.gnuplot rename example/scripts/{scalar-surface.gnuplot => velocity-potential.gnuplot} (82%) diff --git a/example/2D-vortex.F90 b/example/2D-vortex.F90 new file mode 100644 index 0000000..edfe4be --- /dev/null +++ b/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 \ No newline at end of file diff --git a/example/plot-2D-scalar-gradient.F90 b/example/plot-2D-scalar-gradient.F90 deleted file mode 100644 index f23541b..0000000 --- a/example/plot-2D-scalar-gradient.F90 +++ /dev/null @@ -1,51 +0,0 @@ -! Copyright (c) 2026, The Regents of the University of California -! Terms of use are as specified in LICENSE.txt - -module scalar_2D_functions_m - implicit none - - integer, parameter :: space_dimension = 2 - -contains - - pure function biquadratic(x,y) result(z) - double precision, intent(in) :: x(:), y(:) - double precision z(size(x),size(y)) - do concurrent(integer :: j=1:size(y)) default(none) shared(x,y,z) - z(:,j) = 1 - 2*x + 3*x**2 - x*y(j)/5 + 3*y(j)**2 - 2*y(j) - end do - end function - - pure function biquadratic_gradient(x,y) result(gradient) - double precision, intent(in) :: x(:), y(:) - double precision gradient(size(x),size(y),space_dimension) - do concurrent(integer :: i=1:size(x), j=1:size(y)) - gradient(i,j,:) = [-2 + 6*x(i) - y(j)/5, -x(i)/5 + 6*y(j) - 2] - end do - end function - -end module scalar_2D_functions_m - -program scalar_surface - use julienne_m, only : file_t - use scalar_2D_functions_m, only : biquadratic, biquadratic_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 :: expected_gradient_initializer - integer, parameter :: order = 4 - - scalar_2D_initializer => biquadratic - expected_gradient_initializer => biquadratic_gradient - - associate(scalar_2D => scalar_2D_t(scalar_2D_initializer, order=order, cells=[30,20], x_min=[-1D0,1D0], x_max=[9D0,4D0])) - associate(grad_scalar => .grad. scalar_2D, expected_gradient => vector_2D_t(expected_gradient_initializer, mold=scalar_2D)) - associate(scalar_2D_file => scalar_2D%to_file(), grad_scalar_file => grad_scalar%to_file()) - call scalar_2D_file%write_lines("example/scripts/scalar-surface.csv") - call grad_scalar_file%write_lines("example/scripts/gradient-field.csv") - end associate - end associate - end associate - -end program scalar_surface \ No newline at end of file diff --git a/example/scripts/velocities.gnuplot b/example/scripts/velocities.gnuplot new file mode 100644 index 0000000..fb062c1 --- /dev/null +++ b/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 "" diff --git a/example/scripts/scalar-surface.gnuplot b/example/scripts/velocity-potential.gnuplot similarity index 82% rename from example/scripts/scalar-surface.gnuplot rename to example/scripts/velocity-potential.gnuplot index d7418c8..82668b0 100644 --- a/example/scripts/scalar-surface.gnuplot +++ b/example/scripts/velocity-potential.gnuplot @@ -1,11 +1,12 @@ # ============================================================ -# surface.gnuplot -- surface plot from a pre-gridded CSV +# velocity-potential.gnuplot -- surface plot CSV # Line 1: column labels # Lines 2+: x, y, z data with blank lines between x-slices -# Usage: gnuplot scalar-surface.gnuplot +# Usage: gnuplot velocity-potential.gnuplot # ============================================================ -datafile = "scalar-surface.csv" +base_name = "velocity-potential" +datafile = base_name . ".csv" set datafile separator "," # --- 1. Read column headers from line 1 --- @@ -26,7 +27,7 @@ set palette rgbformulae 33,13,10 set ticslevel 0 ; set key off set terminal gif size 800,600 -set output "scalar-surface.gif" +set output base_name . ".gif" splot datafile every ::1 using 1:2:3 with pm3d title "" From 1308b21499a3e8f4000343dcca58162658738b6c Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Sun, 17 May 2026 23:42:05 -0700 Subject: [PATCH 17/36] chore(tensors_1D_m): rm unused associate name --- src/formal/tensors_1D_m.F90 | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/src/formal/tensors_1D_m.F90 b/src/formal/tensors_1D_m.F90 index 4555c6f..59ceae9 100644 --- a/src/formal/tensors_1D_m.F90 +++ b/src/formal/tensors_1D_m.F90 @@ -525,10 +525,7 @@ pure function cell_centers_extended_1D(x_min, x_max, cells) result(x) 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 + x = [x_min, cell_center_locations(x_min, x_max, cells), x_max] end function pure function faces_1D(x_min, x_max, cells) result(x) From e3b26bffaf6c064085db40095163c991d5b7acd9 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Mon, 18 May 2026 22:18:35 -0700 Subject: [PATCH 18/36] fix(lfortran): work around integer-type-spec --- src/formal/differential_operators_1D_m.F90 | 21 ++++++++++++++++++- src/formal/divergence_operator_1D_s.F90 | 24 ---------------------- src/formal/scalar_2D_s.F90 | 17 ++++++++++++++- src/formal/scalar_3D_s.F90 | 20 ++++++++++++++++++ 4 files changed, 56 insertions(+), 26 deletions(-) diff --git a/src/formal/differential_operators_1D_m.F90 b/src/formal/differential_operators_1D_m.F90 index 0b77038..24268d2 100644 --- a/src/formal/differential_operators_1D_m.F90 +++ b/src/formal/differential_operators_1D_m.F90 @@ -160,7 +160,26 @@ pure function negate_and_flip(A) result(Ap) #else -! see divergence_operator_1D_s and gradient_operator_1D_s + pure function negate_and_flip(A) result(Ap) + !! Transform a mimetic matrix upper block into a lower block + double precision, intent(in) :: A(:,:) + double precision, allocatable :: Ap(:,:) + integer row, column + + allocate(Ap, mold=A) + + reverse_elements_within_rows_and_flip_sign: & + do concurrent(row = 1:size(Ap,1)) + Ap(row,:) = -A(row,size(A,2):1:-1) + end do reverse_elements_within_rows_and_flip_sign + + reverse_elements_within_columns: & + do concurrent(column = 1 : size(Ap,2)) + Ap(:,column) = Ap(size(Ap,1):1:-1,column) + end do reverse_elements_within_columns + + end function + #endif diff --git a/src/formal/divergence_operator_1D_s.F90 b/src/formal/divergence_operator_1D_s.F90 index 68c23d4..b3cb100 100644 --- a/src/formal/divergence_operator_1D_s.F90 +++ b/src/formal/divergence_operator_1D_s.F90 @@ -12,30 +12,6 @@ implicit none contains -#if !(HAVE_DO_CONCURRENT_TYPE_SPEC_SUPPORT && HAVE_LOCALITY_SPECIFIER_SUPPORT) - - pure function negate_and_flip(A) result(Ap) - !! Transform a mimetic matrix upper block into a lower block - double precision, intent(in) :: A(:,:) - double precision, allocatable :: Ap(:,:) - integer row, column - - allocate(Ap, mold=A) - - reverse_elements_within_rows_and_flip_sign: & - do concurrent(row = 1:size(Ap,1)) - Ap(row,:) = -A(row,size(A,2):1:-1) - end do reverse_elements_within_rows_and_flip_sign - - reverse_elements_within_columns: & - do concurrent(column = 1 : size(Ap,2)) - Ap(:,column) = Ap(size(Ap,1):1:-1,column) - end do reverse_elements_within_columns - - end function - -#endif - module procedure construct_1D_divergence_operator double precision, allocatable :: Ap(:,:) diff --git a/src/formal/scalar_2D_s.F90 b/src/formal/scalar_2D_s.F90 index e6d23c6..b41f301 100644 --- a/src/formal/scalar_2D_s.F90 +++ b/src/formal/scalar_2D_s.F90 @@ -53,7 +53,7 @@ module procedure scalar_2D_gradient - integer c, i, j + integer c gradient_2D%x_min_ = self%x_min_ gradient_2D%x_max_ = self%x_max_ @@ -62,6 +62,7 @@ allocate(gradient_2D%values_(self%cells_(1)+1, self%cells_(2)+1, space_dimension, 1, 1, 1)) +#if HAVE_DO_CONCURRENT_TYPE_SPEC_SUPPORT && HAVE_LOCALITY_SPECIFIER_SUPPORT 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) @@ -71,6 +72,20 @@ 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 +#else + block + integer i, j + gradient_x_component: & + do concurrent(j=1:size(gradient_2D%values_,2)) + 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(i=1:size(gradient_2D%values_,1)) + 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 + end block +#endif associate(dx => (self%x_max_ - self%x_min_)/self%cells_) gradient_2D%divergence_operator_1D_ = divergence_operator_1D_t(self%order_, dx, self%cells_) diff --git a/src/formal/scalar_3D_s.F90 b/src/formal/scalar_3D_s.F90 index 28e4fc8..de13768 100644 --- a/src/formal/scalar_3D_s.F90 +++ b/src/formal/scalar_3D_s.F90 @@ -66,6 +66,7 @@ allocate(gradient_3D%values_(self%cells_(1)+1, self%cells_(2)+1, self%cells_(3)+1, space_dimension, 1, 1, 1)) +#if HAVE_DO_CONCURRENT_TYPE_SPEC_SUPPORT && HAVE_LOCALITY_SPECIFIER_SUPPORT gradient_x_component: & do concurrent(integer :: j=1:size(gradient_3D%values_,2), k=1:size(gradient_3D%values_,3)) default(none) shared(gradient_3D, self) gradient_3D%values_(:,j,k,1,1,1,1) = self%gradient_operator_1D_(1) .x. self%values_(:,j,k,1,1,1,1) @@ -80,6 +81,25 @@ do concurrent(integer :: i=1:size(gradient_3D%values_,1), j=1:size(gradient_3D%values_,2)) default(none) shared(gradient_3D, self) gradient_3D%values_(i,j,:,3,1,1,1) = self%gradient_operator_1D_(3) .x. self%values_(i,j,:,1,1,1,1) end do gradient_z_component +#else + block + integer i,j,k + gradient_x_component: & + do concurrent(j=1:size(gradient_3D%values_,2), k=1:size(gradient_3D%values_,3)) + gradient_3D%values_(:,j,k,1,1,1,1) = self%gradient_operator_1D_(1) .x. self%values_(:,j,k,1,1,1,1) + end do gradient_x_component + + gradient_y_component: & + do concurrent(i=1:size(gradient_3D%values_,1), k=1:size(gradient_3D%values_,3)) + gradient_3D%values_(i,:,k,2,1,1,1) = self%gradient_operator_1D_(2) .x. self%values_(i,:,k,1,1,1,1) + end do gradient_y_component + + gradient_z_component: & + do concurrent(i=1:size(gradient_3D%values_,1), j=1:size(gradient_3D%values_,2)) + gradient_3D%values_(i,j,:,3,1,1,1) = self%gradient_operator_1D_(3) .x. self%values_(i,j,:,1,1,1,1) + end do gradient_z_component + end block +#endif associate(dx => (self%x_max_ - self%x_min_)/self%cells_) gradient_3D%divergence_operator_1D_ = divergence_operator_1D_t(self%order_, dx, self%cells_) From 89aa9efafd03de40f5b2c94c3442539be5180a45 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Tue, 19 May 2026 18:37:30 -0700 Subject: [PATCH 19/36] Revert "fix(lfortran): work around integer-type-spec" This reverts commit 3b86aa2939686883c9e35307b4dbd3cc5aa207cb. --- src/formal/differential_operators_1D_m.F90 | 21 +------------------ src/formal/divergence_operator_1D_s.F90 | 24 ++++++++++++++++++++++ src/formal/scalar_2D_s.F90 | 17 +-------------- src/formal/scalar_3D_s.F90 | 20 ------------------ 4 files changed, 26 insertions(+), 56 deletions(-) diff --git a/src/formal/differential_operators_1D_m.F90 b/src/formal/differential_operators_1D_m.F90 index 24268d2..0b77038 100644 --- a/src/formal/differential_operators_1D_m.F90 +++ b/src/formal/differential_operators_1D_m.F90 @@ -160,26 +160,7 @@ pure function negate_and_flip(A) result(Ap) #else - pure function negate_and_flip(A) result(Ap) - !! Transform a mimetic matrix upper block into a lower block - double precision, intent(in) :: A(:,:) - double precision, allocatable :: Ap(:,:) - integer row, column - - allocate(Ap, mold=A) - - reverse_elements_within_rows_and_flip_sign: & - do concurrent(row = 1:size(Ap,1)) - Ap(row,:) = -A(row,size(A,2):1:-1) - end do reverse_elements_within_rows_and_flip_sign - - reverse_elements_within_columns: & - do concurrent(column = 1 : size(Ap,2)) - Ap(:,column) = Ap(size(Ap,1):1:-1,column) - end do reverse_elements_within_columns - - end function - +! see divergence_operator_1D_s and gradient_operator_1D_s #endif diff --git a/src/formal/divergence_operator_1D_s.F90 b/src/formal/divergence_operator_1D_s.F90 index b3cb100..68c23d4 100644 --- a/src/formal/divergence_operator_1D_s.F90 +++ b/src/formal/divergence_operator_1D_s.F90 @@ -12,6 +12,30 @@ implicit none contains +#if !(HAVE_DO_CONCURRENT_TYPE_SPEC_SUPPORT && HAVE_LOCALITY_SPECIFIER_SUPPORT) + + pure function negate_and_flip(A) result(Ap) + !! Transform a mimetic matrix upper block into a lower block + double precision, intent(in) :: A(:,:) + double precision, allocatable :: Ap(:,:) + integer row, column + + allocate(Ap, mold=A) + + reverse_elements_within_rows_and_flip_sign: & + do concurrent(row = 1:size(Ap,1)) + Ap(row,:) = -A(row,size(A,2):1:-1) + end do reverse_elements_within_rows_and_flip_sign + + reverse_elements_within_columns: & + do concurrent(column = 1 : size(Ap,2)) + Ap(:,column) = Ap(size(Ap,1):1:-1,column) + end do reverse_elements_within_columns + + end function + +#endif + module procedure construct_1D_divergence_operator double precision, allocatable :: Ap(:,:) diff --git a/src/formal/scalar_2D_s.F90 b/src/formal/scalar_2D_s.F90 index b41f301..e6d23c6 100644 --- a/src/formal/scalar_2D_s.F90 +++ b/src/formal/scalar_2D_s.F90 @@ -53,7 +53,7 @@ module procedure scalar_2D_gradient - integer c + integer c, i, j gradient_2D%x_min_ = self%x_min_ gradient_2D%x_max_ = self%x_max_ @@ -62,7 +62,6 @@ allocate(gradient_2D%values_(self%cells_(1)+1, self%cells_(2)+1, space_dimension, 1, 1, 1)) -#if HAVE_DO_CONCURRENT_TYPE_SPEC_SUPPORT && HAVE_LOCALITY_SPECIFIER_SUPPORT 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) @@ -72,20 +71,6 @@ 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 -#else - block - integer i, j - gradient_x_component: & - do concurrent(j=1:size(gradient_2D%values_,2)) - 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(i=1:size(gradient_2D%values_,1)) - 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 - end block -#endif associate(dx => (self%x_max_ - self%x_min_)/self%cells_) gradient_2D%divergence_operator_1D_ = divergence_operator_1D_t(self%order_, dx, self%cells_) diff --git a/src/formal/scalar_3D_s.F90 b/src/formal/scalar_3D_s.F90 index de13768..28e4fc8 100644 --- a/src/formal/scalar_3D_s.F90 +++ b/src/formal/scalar_3D_s.F90 @@ -66,7 +66,6 @@ allocate(gradient_3D%values_(self%cells_(1)+1, self%cells_(2)+1, self%cells_(3)+1, space_dimension, 1, 1, 1)) -#if HAVE_DO_CONCURRENT_TYPE_SPEC_SUPPORT && HAVE_LOCALITY_SPECIFIER_SUPPORT gradient_x_component: & do concurrent(integer :: j=1:size(gradient_3D%values_,2), k=1:size(gradient_3D%values_,3)) default(none) shared(gradient_3D, self) gradient_3D%values_(:,j,k,1,1,1,1) = self%gradient_operator_1D_(1) .x. self%values_(:,j,k,1,1,1,1) @@ -81,25 +80,6 @@ do concurrent(integer :: i=1:size(gradient_3D%values_,1), j=1:size(gradient_3D%values_,2)) default(none) shared(gradient_3D, self) gradient_3D%values_(i,j,:,3,1,1,1) = self%gradient_operator_1D_(3) .x. self%values_(i,j,:,1,1,1,1) end do gradient_z_component -#else - block - integer i,j,k - gradient_x_component: & - do concurrent(j=1:size(gradient_3D%values_,2), k=1:size(gradient_3D%values_,3)) - gradient_3D%values_(:,j,k,1,1,1,1) = self%gradient_operator_1D_(1) .x. self%values_(:,j,k,1,1,1,1) - end do gradient_x_component - - gradient_y_component: & - do concurrent(i=1:size(gradient_3D%values_,1), k=1:size(gradient_3D%values_,3)) - gradient_3D%values_(i,:,k,2,1,1,1) = self%gradient_operator_1D_(2) .x. self%values_(i,:,k,1,1,1,1) - end do gradient_y_component - - gradient_z_component: & - do concurrent(i=1:size(gradient_3D%values_,1), j=1:size(gradient_3D%values_,2)) - gradient_3D%values_(i,j,:,3,1,1,1) = self%gradient_operator_1D_(3) .x. self%values_(i,j,:,1,1,1,1) - end do gradient_z_component - end block -#endif associate(dx => (self%x_max_ - self%x_min_)/self%cells_) gradient_3D%divergence_operator_1D_ = divergence_operator_1D_t(self%order_, dx, self%cells_) From 12ed6ea12af02514e685614c6c22892bf976ea74 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Tue, 19 May 2026 18:57:53 -0700 Subject: [PATCH 20/36] chore(.gitignore): add scratch dir --- .gitignore | 2 ++ 1 file changed, 2 insertions(+) diff --git a/.gitignore b/.gitignore index c73fb8d..619556d 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,5 @@ +scratch + # Prerequisites *.d From e07c0f2d4c79cae9f835c5fff5c57f4122c60baf Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Fri, 22 May 2026 10:42:34 -0700 Subject: [PATCH 21/36] test(scalar_1D): assert same-order operands --- src/formal/scalar_1D_s.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/formal/scalar_1D_s.F90 b/src/formal/scalar_1D_s.F90 index a7f0b06..84669e8 100644 --- a/src/formal/scalar_1D_s.F90 +++ b/src/formal/scalar_1D_s.F90 @@ -93,7 +93,7 @@ pure function cell_center_locations(x_min, x_max, cells) result(x) module procedure subtract_scalar_1D call_julienne_assert(size(lhs%values_) .equalsExpected. size(rhs%values_)) - call_julienne_assert(lhs%cells_ .equalsExpected. rhs%cells_) + call_julienne_assert(.all.([lhs%cells_,lhs%order_] .equalsExpected. [rhs%cells_,rhs%order_])) call_julienne_assert(.all.([lhs%x_min_,lhs%x_max_] .approximates. [rhs%x_min_,rhs%x_max_] .within. 1D-08)) difference%gradient_operator_1D_ = lhs%gradient_operator_1D_ difference%tensor_1D_t = & @@ -102,7 +102,7 @@ pure function cell_center_locations(x_min, x_max, cells) result(x) module procedure add_scalar_1D call_julienne_assert(size(lhs%values_) .equalsExpected. size(rhs%values_)) - call_julienne_assert(lhs%cells_ .equalsExpected. rhs%cells_) + call_julienne_assert(.all.([lhs%cells_,lhs%order_] .equalsExpected. [rhs%cells_,rhs%order_])) call_julienne_assert(.all.([lhs%x_min_,lhs%x_max_] .approximates. [rhs%x_min_,rhs%x_max_] .within. 1D-08)) total%gradient_operator_1D_ = lhs%gradient_operator_1D_ total%tensor_1D_t = & From ef7181bf0684fbfd033cdd6f8aed7406b1fbf747 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Fri, 22 May 2026 10:53:01 -0700 Subject: [PATCH 22/36] chore(scalar_1D): operand conformabality assertion This commit gathers a pattern of scalar_1D_t operand conformability assertions into one function and calls that function where applicable. --- src/formal/scalar_1D_s.F90 | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/src/formal/scalar_1D_s.F90 b/src/formal/scalar_1D_s.F90 index 84669e8..007b5d0 100644 --- a/src/formal/scalar_1D_s.F90 +++ b/src/formal/scalar_1D_s.F90 @@ -91,19 +91,23 @@ pure function cell_center_locations(x_min, x_max, cells) result(x) ) end procedure - module procedure subtract_scalar_1D + pure logical function conformable(lhs, rhs) + type(scalar_1D_t), intent(in) :: lhs, rhs call_julienne_assert(size(lhs%values_) .equalsExpected. size(rhs%values_)) call_julienne_assert(.all.([lhs%cells_,lhs%order_] .equalsExpected. [rhs%cells_,rhs%order_])) call_julienne_assert(.all.([lhs%x_min_,lhs%x_max_] .approximates. [rhs%x_min_,rhs%x_max_] .within. 1D-08)) + conformable = .true. + end function + + module procedure subtract_scalar_1D + call_julienne_assert(conformable(lhs,rhs)) difference%gradient_operator_1D_ = lhs%gradient_operator_1D_ difference%tensor_1D_t = & tensor_1D_t(values = lhs%values_ - rhs%values_, x_min = rhs%x_min_, x_max = rhs%x_max_, cells = rhs%cells_, order = rhs%order_) end procedure module procedure add_scalar_1D - call_julienne_assert(size(lhs%values_) .equalsExpected. size(rhs%values_)) - call_julienne_assert(.all.([lhs%cells_,lhs%order_] .equalsExpected. [rhs%cells_,rhs%order_])) - call_julienne_assert(.all.([lhs%x_min_,lhs%x_max_] .approximates. [rhs%x_min_,rhs%x_max_] .within. 1D-08)) + call_julienne_assert(conformable(lhs,rhs)) total%gradient_operator_1D_ = lhs%gradient_operator_1D_ total%tensor_1D_t = & tensor_1D_t(values = lhs%values_ + rhs%values_, x_min = rhs%x_min_, x_max = rhs%x_max_, cells = rhs%cells_, order = rhs%order_) From bbc7c01c92991fcbb7d5dbed0c37142a46765773 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Fri, 22 May 2026 15:02:30 -0700 Subject: [PATCH 23/36] feat(scalar_1D): add scalar-scalar multiplication operator Co-Authored-By: Claude Opus 4.6 --- src/formal/scalar_1D_s.F90 | 7 +++++++ src/formal/tensors_1D_m.F90 | 11 ++++++++++- 2 files changed, 17 insertions(+), 1 deletion(-) diff --git a/src/formal/scalar_1D_s.F90 b/src/formal/scalar_1D_s.F90 index 007b5d0..c898797 100644 --- a/src/formal/scalar_1D_s.F90 +++ b/src/formal/scalar_1D_s.F90 @@ -106,6 +106,13 @@ pure logical function conformable(lhs, rhs) tensor_1D_t(values = lhs%values_ - rhs%values_, x_min = rhs%x_min_, x_max = rhs%x_max_, cells = rhs%cells_, order = rhs%order_) end procedure + module procedure multiply_1D_scalars + call_julienne_assert(conformable(lhs,rhs)) + lhs_x_rhs%gradient_operator_1D_ = lhs%gradient_operator_1D_ + lhs_x_rhs%tensor_1D_t = & + tensor_1D_t(values = lhs%values_ * rhs%values_, x_min = rhs%x_min_, x_max = rhs%x_max_, cells = rhs%cells_, order = rhs%order_) + end procedure + module procedure add_scalar_1D call_julienne_assert(conformable(lhs,rhs)) total%gradient_operator_1D_ = lhs%gradient_operator_1D_ diff --git a/src/formal/tensors_1D_m.F90 b/src/formal/tensors_1D_m.F90 index 59ceae9..db904c7 100644 --- a/src/formal/tensors_1D_m.F90 +++ b/src/formal/tensors_1D_m.F90 @@ -89,7 +89,7 @@ pure module function construct_1D_tensor_from_components(values, x_min, x_max, c generic :: operator(-) => subtract_scalar_1D generic :: operator(+) => add_scalar_1D generic :: operator(/) => divide_by_integer - generic :: operator(*) => premultiply_double, postmultiply_double, premultiply_integer, postmultiply_integer + generic :: operator(*) => premultiply_double, postmultiply_double, premultiply_integer, postmultiply_integer, multiply_1D_scalars generic :: operator(**) => exponentiate generic :: operator(.grad.) => grad generic :: operator(.laplacian.) => laplacian @@ -98,6 +98,7 @@ pure module function construct_1D_tensor_from_components(values, x_min, x_max, c procedure, non_overridable, private :: scalar_1D_values procedure, non_overridable, private :: scalar_1D_grid procedure, non_overridable, private :: divide_by_integer + procedure, non_overridable, private :: multiply_1D_scalars procedure, non_overridable, private, pass(rhs) :: premultiply_double procedure, non_overridable, private :: postmultiply_double procedure, non_overridable, private, pass(rhs) :: premultiply_integer @@ -392,6 +393,14 @@ pure module function add_scalar_1D(lhs, rhs) result(total) type(scalar_1D_t) total end function + pure module function multiply_1D_scalars(lhs, rhs) result(lhs_x_rhs) + !! Result is the product of scalar_1D_t lhs and rhs + implicit none + class(scalar_1D_t), intent(in) :: lhs + type(scalar_1D_t), intent(in) :: rhs + type(scalar_1D_t) lhs_x_rhs + end function + pure module function subtract_scalar_1D(lhs, rhs) result(difference) !! Result is the difference between scalar_1D_t lhs and rhs implicit none From 96e6d7d2aac1e65f911ce10dab32db6b4792afab Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Fri, 22 May 2026 15:56:15 -0700 Subject: [PATCH 24/36] chore(.gitignore): add macOS .DS_Store file --- .gitignore | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/.gitignore b/.gitignore index 619556d..780a618 100644 --- a/.gitignore +++ b/.gitignore @@ -1,5 +1,9 @@ +# Scratch directory scratch +# macOS folder display format +.DS_Store + # Prerequisites *.d From 4e23d25df5e99024d8472a5350a8ed4b58612b97 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Sun, 24 May 2026 22:20:45 -0700 Subject: [PATCH 25/36] feat(2D .div.): add/test 2D divergence type/ops --- src/formal/divergence_1D_s.F90 | 20 ++---- src/formal/divergence_2D_s.F90 | 78 ++++++++++++++++++++++ src/formal/scalar_1D_s.F90 | 35 ---------- src/formal/tensors_1D_m.F90 | 24 ++++--- src/formal/tensors_2D_m.F90 | 78 ++++++++++++++++++++-- src/formal/vector_2D_s.F90 | 25 ++++++- src/formal_m.f90 | 4 +- test/divergence_operator_1D_test_m.F90 | 7 +- test/driver.f90 | 2 + test/vector_2D_test_m.F90 | 90 ++++++++++++++++++++++++++ 10 files changed, 298 insertions(+), 65 deletions(-) create mode 100644 src/formal/divergence_2D_s.F90 create mode 100644 test/vector_2D_test_m.F90 diff --git a/src/formal/divergence_1D_s.F90 b/src/formal/divergence_1D_s.F90 index b7aa0d4..608ebbd 100644 --- a/src/formal/divergence_1D_s.F90 +++ b/src/formal/divergence_1D_s.F90 @@ -12,20 +12,10 @@ contains -#ifdef __GFORTRAN__ - - pure function cell_center_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 + dx/2. + [((cell-1)*dx, cell = 1, cells)] - end associate - end function - -#endif + module procedure construct_1D_divergence_constant + integer i + divergence_1D%tensor_1D_t = tensor_1D_t([(constant, i=1,cells)], x_min, x_max, cells, order) + end procedure module procedure premultiply_scalar_1D call_julienne_assert(size(scalar_1D%values_) .equalsExpected. size(divergence_1D%values_) + 2) @@ -48,7 +38,7 @@ pure function cell_center_locations(x_min, x_max, cells) result(x) end procedure module procedure divergence_1D_grid - cell_centers = cell_center_locations(self%x_min_, self%x_max_, self%cells_) + cell_centers = cell_centers_1D(self%x_min_, self%x_max_, self%cells_) end procedure module procedure divergence_1D_weights diff --git a/src/formal/divergence_2D_s.F90 b/src/formal/divergence_2D_s.F90 new file mode 100644 index 0000000..d5aafce --- /dev/null +++ b/src/formal/divergence_2D_s.F90 @@ -0,0 +1,78 @@ +! 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) divergence_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, divergence_1D_t, cell_centers_1D + use julienne_m, only : string_t, operator(.csv.) + implicit none + +contains + + module procedure divergence_2D_values + divergence_values = self%values_(:,:,1,1,1,1) + end procedure + + module procedure divergence_2D_grid + associate(divergence_1D => divergence_1D_t( & + constant = 0D0 & + ,cells = self%cells_(direction) & + ,x_min = self%x_min_(direction) & + ,x_max = self%x_max_(direction) & + ,order = self%order_ & + )) + divergence_grid_1D = divergence_1D%grid() + end associate + end procedure + + module procedure construct_2D_divergence_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_1D(x_min(1), x_max(1), cells(1)), y => cell_centers_1D(x_min(2), x_max(2), cells(2))) + divergence_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 & + ) + end associate + end procedure + + module procedure construct_2D_divergence_from_vector_mold + divergence_2D = divergence_2D_t(initializer, cells = mold%cells_, x_min = mold%x_min_, x_max = mold%x_max_, order = mold%order_) + end procedure + + module procedure divergence_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,divergence")]) + 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 divergence_2D_s \ No newline at end of file diff --git a/src/formal/scalar_1D_s.F90 b/src/formal/scalar_1D_s.F90 index c898797..d799756 100644 --- a/src/formal/scalar_1D_s.F90 +++ b/src/formal/scalar_1D_s.F90 @@ -21,9 +21,6 @@ contains - -#ifndef __GFORTRAN__ - module procedure construct_1D_scalar_from_function call_julienne_assert(x_max .greaterThan. x_min) call_julienne_assert(cells .isAtLeast. 2*order) @@ -34,38 +31,6 @@ scalar_1D%gradient_operator_1D_ = gradient_operator_1D_t(k=order, dx=(x_max - x_min)/cells, cells=cells) end procedure -#else - - pure module function construct_1D_scalar_from_function(initializer, order, cells, x_min, x_max) result(scalar_1D) - procedure(scalar_1D_initializer_i), pointer :: initializer - integer, intent(in) :: order !! order of accuracy - integer, intent(in) :: cells !! number of grid cells spanning the domain - double precision, intent(in) :: x_min !! grid location minimum - double precision, intent(in) :: x_max !! grid location maximum - type(scalar_1D_t) scalar_1D - - call_julienne_assert(x_max .greaterThan. x_min) - call_julienne_assert(cells .isAtLeast. 2*order) - - 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) - end function - - pure function cell_center_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 + dx/2. + [((cell-1)*dx, cell = 1, cells)] - end associate - end function - -#endif - module procedure construct_1D_scalar_constant integer i diff --git a/src/formal/tensors_1D_m.F90 b/src/formal/tensors_1D_m.F90 index db904c7..67a5ec9 100644 --- a/src/formal/tensors_1D_m.F90 +++ b/src/formal/tensors_1D_m.F90 @@ -232,6 +232,20 @@ pure module function construct_1D_vector_constant(constant, order, cells, x_min, procedure, non_overridable, private :: divergence_1D_grid end type + interface divergence_1D_t + + pure module function construct_1D_divergence_constant(constant, order, cells, x_min, x_max) result(divergence_1D) + implicit none + double precision, intent(in) :: constant !! scalar value + integer, intent(in) :: order !! order of accuracy + integer, intent(in) :: cells !! number of grid cells spanning the domain + double precision, intent(in) :: x_min !! grid location minimum + double precision, intent(in) :: x_max !! grid location maximum + type(divergence_1D_t) divergence_1D + end function + + end interface + type, extends(tensor_1D_t) :: scalar_x_divergence_1D_t !! product of a 1D scalar field and a 1D divergence field private @@ -534,7 +548,7 @@ pure function cell_centers_extended_1D(x_min, x_max, cells) result(x) integer, intent(in) :: cells double precision, allocatable:: x(:) integer cell - x = [x_min, cell_center_locations(x_min, x_max, cells), x_max] + x = [x_min, cell_centers_1D(x_min, x_max, cells), x_max] end function pure function faces_1D(x_min, x_max, cells) result(x) @@ -542,25 +556,19 @@ pure function faces_1D(x_min, x_max, cells) result(x) integer, intent(in) :: cells double precision, allocatable:: x(:) integer cell - associate(dx => (x_max - x_min)/cells) x = [x_min, x_min + [(cell*dx, cell = 1, cells-1)], x_max] end associate end function -#ifndef __GFORTRAN__ - - pure function cell_center_locations(x_min, x_max, cells) result(x) + pure function cell_centers_1D(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 + dx/2. + [((cell-1)*dx, cell = 1, cells)] end associate end function -#endif - end module tensors_1D_m diff --git a/src/formal/tensors_2D_m.F90 b/src/formal/tensors_2D_m.F90 index 3f1a83e..f5c5647 100644 --- a/src/formal/tensors_2D_m.F90 +++ b/src/formal/tensors_2D_m.F90 @@ -7,16 +7,16 @@ module tensors_2D_m !! https://doi.org/10.1016/j.cam.2019.06.042. use differential_operators_1D_m, only : gradient_operator_1D_t, divergence_operator_1D_t use julienne_m, only : file_t - implicit none private - public :: scalar_2D_t public :: vector_2D_t public :: gradient_2D_t + public :: divergence_2D_t public :: scalar_2D_initializer_i public :: vector_2D_initializer_i + public :: divergence_2D_initializer_i integer, parameter :: space_dimension = 2 @@ -29,6 +29,13 @@ pure function scalar_2D_initializer_i(x,y) result(f) double precision f(size(x),size(y)) end function + pure function divergence_2D_initializer_i(x,y) result(f) + !! Sampling function for initializing a scalar_2D_t object + implicit none + double precision, intent(in) :: x(:), y(:) + double precision f(size(x),size(y)) + end function + pure function vector_2D_initializer_i(x,y) result(v) !! Sampling function for initializing a vector_2D_t object import space_dimension @@ -111,9 +118,11 @@ pure module function construct_2D_scalar_from_mold(initializer, mold) result(sca generic :: values => vector_2D_values generic :: to_file => vector_2D_to_file generic :: grid => vector_2D_grid + generic :: operator(.div.) => vector_2D_divergence procedure, non_overridable, private :: vector_2D_values procedure, non_overridable, private :: vector_2D_to_file procedure, non_overridable, private :: vector_2D_grid + procedure, non_overridable, private :: vector_2D_divergence end type interface vector_2D_t @@ -129,7 +138,7 @@ pure module function construct_2D_vector_from_function(initializer, order, cells double precision, intent(in) :: x_max(:) !! grid location maxima type(vector_2D_t) vector_2D end function - + pure module function construct_2D_vector_from_vector_mold(initializer, mold) result(vector_2D) !! Result is a 2D vector with values initialized by the provided procedure pointer sampled on the !! same grid as the mold @@ -154,6 +163,40 @@ pure module function construct_2D_vector_from_scalar_mold(initializer, mold) res !! A 2D mimetic gradient vector field abstraction with a public method that produces corresponding numerical quadrature weights end type + type, extends(tensor_2D_t) :: divergence_2D_t + !! A 2D mimetic divergence field abstraction with a public method that produces corresponding numerical quadrature weights + contains + generic :: values => divergence_2D_values + generic :: grid => divergence_2D_grid + generic :: to_file => divergence_2D_to_file + procedure, private, non_overridable :: divergence_2D_values + procedure, private, non_overridable :: divergence_2D_grid + procedure, private, non_overridable :: divergence_2D_to_file + end type + + interface divergence_2D_t + + pure module function construct_2D_divergence_from_function(initializer, order, cells, x_min, x_max) result(divergence_2D) + !! Result is a 2D divergence initialized by sampling the initializer at cell centers defined by the other arguments + implicit none + procedure(scalar_2D_initializer_i), pointer, intent(in) :: initializer + integer, intent(in) :: order !! order of accuracy + integer, intent(in) :: cells(:) !! number of grid cells spanning each spatial direction + double precision, intent(in) :: x_min(:) !! grid location minima + double precision, intent(in) :: x_max(:) !! grid location maxima + type(divergence_2D_t) divergence_2D + end function + + pure module function construct_2D_divergence_from_vector_mold(initializer, mold) result(divergence_2D) + !! Result is a 2D divergence initialized by sampling the initializer on cell centers defined by the mold + implicit none + procedure(divergence_2D_initializer_i), pointer, intent(in) :: initializer + type(vector_2D_t), intent(in) :: mold + type(divergence_2D_t) divergence_2D + end function + + end interface + interface pure module function scalar_2D_values(self) result(scalar_values) @@ -173,7 +216,14 @@ pure module function vector_2D_grid(self, direction) result(vector_grid_1D) !! Result array contains scalar grid locations along the requested spatial direction class(vector_2D_t), intent(in) :: self integer, intent(in) :: direction - double precision, allocatable :: vector_grid_1D(:) !! grid points along one the requested coordinate direction + double precision, allocatable :: vector_grid_1D(:) !! grid points along the requested coordinate direction + end function + + pure module function divergence_2D_grid(self, direction) result(divergence_grid_1D) + !! Result array contains divergence grid locations along the requested spatial direction + class(divergence_2D_t), intent(in) :: self + integer, intent(in) :: direction + double precision, allocatable :: divergence_grid_1D(:) !! grid points along the requested coordinate direction end function pure module function vector_2D_values(self) result(vector_values) @@ -182,6 +232,12 @@ pure module function vector_2D_values(self) result(vector_values) double precision, allocatable :: vector_values(:,:,:) end function + pure module function divergence_2D_values(self) result(divergence_values) + !! Vector values getter + class(divergence_2D_t), intent(in) :: self + double precision, allocatable :: divergence_values(:,:) + end function + pure module function scalar_2D_gradient(self) result(gradient_2D) !! Result is mimetic gradient of the scalar_2D_t "self" implicit none @@ -189,6 +245,13 @@ pure module function scalar_2D_gradient(self) result(gradient_2D) type(gradient_2D_t) gradient_2D end function + pure module function vector_2D_divergence(self) result(divergence_2D) + !! Result is mimetic divergence of the scalar_2D_t "self" + implicit none + class(vector_2D_t), intent(in) :: self + type(divergence_2D_t) divergence_2D + end function + pure module function scalar_2D_to_file(self) result(file) !! Result is a file_t object containing the grid points and the corresponding scalar values implicit none @@ -196,6 +259,13 @@ pure module function scalar_2D_to_file(self) result(file) type(file_t) file end function + pure module function divergence_2D_to_file(self) result(file) + !! Result is a file_t object containing the grid points and the corresponding divergence values + implicit none + class(divergence_2D_t), intent(in) :: self + type(file_t) file + end function + pure module function vector_2D_to_file(self) result(file) !! Result is a file_t object containing the grid points and the corresponding vector components implicit none diff --git a/src/formal/vector_2D_s.F90 b/src/formal/vector_2D_s.F90 index 73f9554..f4ad4da 100644 --- a/src/formal/vector_2D_s.F90 +++ b/src/formal/vector_2D_s.F90 @@ -13,7 +13,6 @@ ,operator(.isAtLeast.) & ,string_t use tensors_1D_m, only : faces_1D, vector_1D_t - use differential_operators_1D_m, only : divergence_operator_1D_t implicit none contains @@ -114,4 +113,28 @@ end associate end procedure + module procedure vector_2D_divergence + + integer c, i, j + + divergence_2D%x_min_ = self%x_min_ + divergence_2D%x_max_ = self%x_max_ + divergence_2D%cells_ = self%cells_ + divergence_2D%order_ = self%order_ + + allocate(divergence_2D%values_(self%cells_(1)+2, self%cells_(2)+2, 1, 1, 1, 1)) + + divergence_x_term: & + do concurrent(integer :: j=1:size(divergence_2D%values_,2)) default(none) shared(divergence_2D, self) + divergence_2D%values_(:,j,1,1,1,1) = self%divergence_operator_1D_(1) .x. self%values_(:,j,1,1,1,1) + end do divergence_x_term + + add_y_term: & + do concurrent(integer :: i=1:size(divergence_2D%values_,1)) default(none) shared(divergence_2D, self) + divergence_2D%values_(i,:,2,1,1,1) = divergence_2D%values_(i,:,1,1,1,1) + & + (self%divergence_operator_1D_(2) .x. self%values_(i,:,1,1,1,1)) + end do add_y_term + + end procedure + end submodule vector_2D_s \ No newline at end of file diff --git a/src/formal_m.f90 b/src/formal_m.f90 index 86d4319..c6a07f3 100644 --- a/src/formal_m.f90 +++ b/src/formal_m.f90 @@ -21,9 +21,11 @@ module formal_m use tensors_2D_m, only : & scalar_2D_t & ! discrete 2D scalar field derived type ,vector_2D_t & ! discrete 2D vector field derived type + ,divergence_2D_t & ! discrete 2D divergence field derived type ,gradient_2D_t & ! result of `.grad. s` for a scalar_2D_t s ,scalar_2D_initializer_i & ! scalar_2D_t initializer abstract interface - ,vector_2D_initializer_i ! vector_2D_t initializar abstract interface + ,vector_2D_initializer_i & ! vector_2D_t initializar abstract interface + ,divergence_2D_initializer_i ! divergence_2D_t initializar abstract interface use tensors_3D_m, only : & scalar_3D_t & ! discrete 3D scalar field derived type diff --git a/test/divergence_operator_1D_test_m.F90 b/test/divergence_operator_1D_test_m.F90 index 87adb5a..9a37857 100644 --- a/test/divergence_operator_1D_test_m.F90 +++ b/test/divergence_operator_1D_test_m.F90 @@ -10,6 +10,7 @@ module divergence_operator_1D_test_m ,operator(.all.) & ,operator(.also.) & ,operator(.approximates.) & + ,operator(.equalsExpected.) & ,operator(.within.) & ,passing_test & ,string_t & @@ -126,6 +127,9 @@ function check_2nd_order_div_sinusoid_convergence() result(test_diagnosis) ,div_fine => .div. vector_1D_t(vector_1D_initializer , order=order_desired, cells=fine_cells , x_min=0D0, x_max=2*pi) & ) #endif + test_diagnosis = passing_test() + test_diagnosis = test_diagnosis .also. (size(div_coarse%values()) .equalsExpected. coarse_cells) + test_diagnosis = test_diagnosis .also. (size(div_fine%values()) .equalsExpected. fine_cells) associate( & x_coarse => div_coarse%grid() & ,x_fine => div_fine%grid()) @@ -135,7 +139,6 @@ function check_2nd_order_div_sinusoid_convergence() result(test_diagnosis) ,div_coarse_values => div_coarse%values() & ,div_fine_values => div_fine%values() & ) - test_diagnosis = passing_test() test_diagnosis = test_diagnosis .also. (.all. (div_coarse_values .approximates. grad_coarse .within. rough_tolerance)) & // " (coarse-grid 2nd-order .div. [sin(x) + cos(x)])" test_diagnosis = test_diagnosis .also. (.all. (div_fine_values .approximates. grad_fine .within. rough_tolerance)) & @@ -183,6 +186,8 @@ function check_4th_order_div_sinusoid_convergence() result(test_diagnosis) ) test_diagnosis = passing_test() + test_diagnosis = test_diagnosis .also. (size(div_coarse_values) .equalsExpected. coarse_cells) + test_diagnosis = test_diagnosis .also. (size(div_fine_values) .equalsExpected. fine_cells) test_diagnosis = test_diagnosis .also. (.all. (div_coarse_values .approximates. div_coarse_expected .within. loose_tolerance)) & // " (coarse-grid 4th-order .div. [sin(x) + cos(x)])" test_diagnosis = test_diagnosis .also. (.all. (div_fine_values .approximates. div_fine_expected .within. loose_tolerance)) & diff --git a/test/driver.f90 b/test/driver.f90 index bffe1e5..8cc6f35 100644 --- a/test/driver.f90 +++ b/test/driver.f90 @@ -12,6 +12,7 @@ program test_suite_driver use scalar_1D_test_m, only : scalar_1D_test_t use scalar_2D_test_m, only : scalar_2D_test_t use scalar_3D_test_m, only : scalar_3D_test_t + use vector_2D_test_m, only : vector_2D_test_t implicit none associate(test_harness => test_harness_t([ & @@ -23,6 +24,7 @@ program test_suite_driver ,test_fixture_t(scalar_1D_test_t()) & ,test_fixture_t(scalar_2D_test_t()) & ,test_fixture_t(scalar_3D_test_t()) & + ,test_fixture_t(vector_2D_test_t()) & ])) call test_harness%report_results end associate diff --git a/test/vector_2D_test_m.F90 b/test/vector_2D_test_m.F90 new file mode 100644 index 0000000..0131a0e --- /dev/null +++ b/test/vector_2D_test_m.F90 @@ -0,0 +1,90 @@ +! Copyright (c) 2026, The Regents of the University of California +! Terms of use are as specified in LICENSE.txt + +module vector_2D_test_m + use julienne_m, only : & + operator(//) & + ,operator(.all.) & + ,operator(.also.) & + ,operator(.approximates.) & + ,operator(.within.) & + ,passing_test & + ,string_t & + ,test_description_t & + ,test_diagnosis_t & + ,test_result_t & + ,test_t & + ,usher + use formal_m, only : & + vector_2D_t & + ,vector_2D_initializer_i & + ,divergence_2D_t & + ,divergence_2D_initializer_i + + implicit none + + type, extends(test_t) :: vector_2D_test_t + contains + procedure, nopass :: subject + procedure, nopass :: results + end type + + double precision, parameter :: tolerance = 5D-2 + integer, parameter :: space_dimension = 2 + +contains + + pure function subject() result(test_subject) + character(len=:), allocatable :: test_subject + test_subject = 'The vector_2D_t derived type' + end function + + function results() result(test_results) + type(vector_2D_test_t) vector_2D_test + type(test_result_t), allocatable :: test_results(:) + + test_results = vector_2D_test%run([ & + test_description_t('computing the divergence of a vector field', usher(check_divergence)) & + ]) + end function + + pure function v(x,y) result(z) + double precision, intent(in) :: x(:), y(:) + double precision z(size(x),size(y),space_dimension) + do concurrent(integer :: i=1:size(x), j=1:size(y)) default(none) shared(x,y,z) + z(i,j,:) = [0D0,0D0] + end do + end function + + pure function div_v(x,y) result(divergence) + double precision, intent(in) :: x(:), y(:) + double precision divergence(size(x),size(y)) + do concurrent(integer :: i=1:size(x), j=1:size(y)) default(none) shared(divergence,x,y) + divergence(i,j) = 0D0 + end do + end function + + function check_divergence() result(test_diagnosis) + type(test_diagnosis_t) test_diagnosis + procedure(vector_2D_initializer_i), pointer :: vector_2D_initializer + procedure(divergence_2D_initializer_i), pointer :: expected_divergence_initializer + integer order + + vector_2D_initializer => v + expected_divergence_initializer => div_v + test_diagnosis = passing_test() + + do order = 2, 4, 2 + associate(vector_2D => vector_2D_t(vector_2D_initializer, order=order, cells=[30,20], x_min=[-1D0,1D0], x_max=[9D0,4D0])) + associate(div_vector => .div. vector_2D) + associate(expected_divergence => divergence_2D_t(expected_divergence_initializer, mold=vector_2D)) + test_diagnosis = test_diagnosis .also. & + (.all. (div_vector%values() .approximates. expected_divergence%values() .within. tolerance)) & + // string_t(" for order ") // string_t(order) + end associate + end associate + end associate + end do + end function + +end module vector_2D_test_m \ No newline at end of file From 6e94a7e10afaa268451eec162d6dd9ae7dd0c8a6 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Sun, 24 May 2026 23:01:17 -0700 Subject: [PATCH 26/36] chore(gfortran): rm workarounds --- example/burgers-1D.F90 | 9 --- example/div-grad-laplacian-1D.F90 | 43 -------------- example/extended-gauss-divergence.F90 | 22 ------- example/print-assembled-1D-operators.F90 | 3 - src/formal/divergence_1D_s.F90 | 4 -- src/formal/gradient_1D_s.F90 | 4 -- src/formal/scalar_1D_s.F90 | 4 -- src/formal/vector_1D_s.F90 | 23 -------- test/divergence_operator_1D_test_m.F90 | 37 ------------ test/gradient_operator_1D_test_m.F90 | 75 ------------------------ test/laplacian_operator_1D_test_m.F90 | 28 +-------- test/scalar_1D_test_m.F90 | 19 ------ 12 files changed, 1 insertion(+), 270 deletions(-) diff --git a/example/burgers-1D.F90 b/example/burgers-1D.F90 index a62b191..9ae6076 100644 --- a/example/burgers-1D.F90 +++ b/example/burgers-1D.F90 @@ -33,9 +33,6 @@ program burgers_1D use iso_fortran_env, only : output_unit implicit none -#ifdef __GFORTRAN__ - procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer -#endif character(len=:), allocatable :: order_string type(command_line_t) command_line integer order @@ -68,9 +65,7 @@ program burgers_1D block -#ifndef __GFORTRAN__ procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer -#endif double precision, parameter :: pi = acos(-1D0), nu=1D0, t_final=0.6D0 double precision, allocatable :: u_surface(:,:), time(:) double precision dt @@ -149,10 +144,6 @@ program burgers_1D end associate end block -#ifdef __GFORTRAN__ - stop -#endif - contains pure function d_dt(u, nu) result(du_dt) diff --git a/example/div-grad-laplacian-1D.F90 b/example/div-grad-laplacian-1D.F90 index 00caf76..6a8f749 100644 --- a/example/div-grad-laplacian-1D.F90 +++ b/example/div-grad-laplacian-1D.F90 @@ -42,9 +42,6 @@ program div_grad_laplacian_1D use functions_m, only : f, df_dx, d2f_dx2 use julienne_m, only : file_t, string_t, operator(.separatedBy.), command_line_t use formal_m, only : scalar_1D_t, scalar_1D_initializer_i -#ifdef __GFORTRAN__ - use formal_m, only : vector_1D_t, laplacian_1D_t, gradient_1D_t -#endif implicit none procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => f @@ -93,8 +90,6 @@ program div_grad_laplacian_1D contains -#ifndef __GFORTRAN__ - subroutine output(order) integer, intent(in) :: order @@ -125,44 +120,6 @@ subroutine output(order) end associate end subroutine -#else - - subroutine output(order) - integer, intent(in) :: order - - type(scalar_1D_t) s - type(gradient_1D_t) grad_s - type(laplacian_1D_t) laplacian_s - type(file_t) s_table, grad_s_table, laplacian_s_table - double precision, allocatable,dimension(:) :: s_grid, grad_s_grid, laplacian_s_grid - - s = scalar_1D_t(scalar_1D_initializer, order=order, cells=20, x_min=0D0, x_max=20D0) - grad_s = .grad. s - laplacian_s = .laplacian. s - - s_grid = s%grid() - grad_s_grid = grad_s%grid() - laplacian_s_grid = laplacian_s%grid() - - 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 subroutine - -#endif - pure function tabulate(headings, abscissa, expected, actual) result(file) double precision, intent(in), dimension(:) :: abscissa, expected, actual type(string_t), intent(in) :: headings(:) diff --git a/example/extended-gauss-divergence.F90 b/example/extended-gauss-divergence.F90 index 89fc8fb..5a57c54 100644 --- a/example/extended-gauss-divergence.F90 +++ b/example/extended-gauss-divergence.F90 @@ -61,15 +61,8 @@ program extended_gauss_divergence call execute_command_line("grep '<-- scalar' example/extended-gauss-divergence.F90 | grep -v execute_command", wait=.true.) call execute_command_line("grep '<-- vector' example/extended-gauss-divergence.F90 | grep -v execute_command", wait=.true.) -#ifdef __GFORTRAN__ - command_line_arguments: & - block - type(numerical_arguments_t) args - args = get_numerical_arguments() -#else command_line_arguments: & associate(args => get_numerical_arguments()) -#endif text_flags: & associate(flags => text_flags_t( & div_ = command_line%argument_present( [ character(len=len("--div" )) :: "--div" , "-d" ] ) & @@ -114,39 +107,24 @@ program extended_gauss_divergence end associate integrand_factors end associate print_all end associate text_flags -#ifndef __GFORTRAN__ end associate command_line_arguments -#else - end block command_line_arguments -#endif contains function get_numerical_arguments() result(numerical_arguments) type(numerical_arguments_t) numerical_arguments -#ifdef __GFORTRAN__ - character(len=:), allocatable :: cells_string, order_string, x_min_string, x_max_string - cells_string = command_line%flag_value("--cells") - order_string = command_line%flag_value("--order") - x_min_string = command_line%flag_value("--x_min") - x_max_string = command_line%flag_value("--x_max") -#else associate( & cells_string => command_line%flag_value("--cells") & ,order_string => command_line%flag_value("--order") & ,x_min_string => command_line%flag_value("--x_min") & ,x_max_string => command_line%flag_value("--x_max") & ) -#endif if (len(cells_string)/=0) read(cells_string,*) numerical_arguments%cells_ if (len(order_string)/=0) read(order_string,*) numerical_arguments%order_ if (len(x_min_string)/=0) read(x_min_string,*) numerical_arguments%x_min_ if (len(x_max_string)/=0) read(x_max_string,*) numerical_arguments%x_max_ -#ifndef __GFORTRAN__ end associate -#endif - end function end program diff --git a/example/print-assembled-1D-operators.F90 b/example/print-assembled-1D-operators.F90 index dad3844..29b1831 100644 --- a/example/print-assembled-1D-operators.F90 +++ b/example/print-assembled-1D-operators.F90 @@ -40,9 +40,6 @@ program print_assembled_1D_operators if (print_all .or. (divergence .and. len(order)==0) .or. (divergence .and. order=="4")) call print_divergence_operator(k=4, dx=1D0, m=16) end associate default_usage -#ifdef __GFORTRAN__ - stop -#endif end associate command_line_settings contains diff --git a/src/formal/divergence_1D_s.F90 b/src/formal/divergence_1D_s.F90 index 608ebbd..bb52f6f 100644 --- a/src/formal/divergence_1D_s.F90 +++ b/src/formal/divergence_1D_s.F90 @@ -21,11 +21,7 @@ call_julienne_assert(size(scalar_1D%values_) .equalsExpected. size(divergence_1D%values_) + 2) scalar_x_divergence_1D%tensor_1D_t = & tensor_1D_t(scalar_1D%values_(2:size(scalar_1D%values_)-1) * divergence_1D%values_, scalar_1D%x_min_, scalar_1D%x_max_, scalar_1D%cells_, scalar_1D%order_) -#ifndef __GFORTRAN__ scalar_x_divergence_1D%weights_ = divergence_1D%weights() -#else - scalar_x_divergence_1D%weights_ = divergence_1D%divergence_1D_weights() -#endif call_julienne_assert(size(scalar_x_divergence_1D%weights_) .equalsExpected. size(divergence_1D%values_)+2) end procedure diff --git a/src/formal/gradient_1D_s.F90 b/src/formal/gradient_1D_s.F90 index e3589f1..010a0ef 100644 --- a/src/formal/gradient_1D_s.F90 +++ b/src/formal/gradient_1D_s.F90 @@ -53,11 +53,7 @@ ,cells = gradient_1D%cells_ & ,order = gradient_1D%order_ & ) -#ifndef __GFORTRAN__ vector_dot_gradient_1D%weights_ = gradient_1D%weights() -#else - vector_dot_gradient_1D%weights_ = gradient_1D%gradient_1D_weights() -#endif end procedure end submodule gradient_1D_s diff --git a/src/formal/scalar_1D_s.F90 b/src/formal/scalar_1D_s.F90 index d799756..a785d08 100644 --- a/src/formal/scalar_1D_s.F90 +++ b/src/formal/scalar_1D_s.F90 @@ -170,11 +170,7 @@ pure logical function conformable(lhs, rhs) module procedure laplacian -#ifndef __GFORTRAN__ laplacian_1D%divergence_1D_t = .div. (.grad. self) -#else - laplacian_1D%divergence_1D_t = div(grad(self)) -#endif associate(divergence_operator_1D => divergence_operator_1D_t(self%order_, (self%x_max_ - self%x_min_)/self%cells_, self%cells_)) laplacian_1D%boundary_depth_ = divergence_operator_1D%submatrix_A_rows() + 1 diff --git a/src/formal/vector_1D_s.F90 b/src/formal/vector_1D_s.F90 index 1467d38..c8cf0db 100644 --- a/src/formal/vector_1D_s.F90 +++ b/src/formal/vector_1D_s.F90 @@ -27,8 +27,6 @@ v_dot_dS%divergence_operator_1D_ = vector_1D%divergence_operator_1D_ end procedure -#ifndef __GFORTRAN__ - module procedure construct_1D_vector_from_function call_julienne_assert(x_max .greaterThan. x_min) call_julienne_assert(cells .isAtLeast. 2*order+1) @@ -39,27 +37,6 @@ vector_1D%divergence_operator_1D_ = divergence_operator_1D_t(k=order, dx=(x_max - x_min)/cells, cells=cells) end procedure -#else - - pure module function construct_1D_vector_from_function(initializer, order, cells, x_min, x_max) result(vector_1D) - procedure(vector_1D_initializer_i), pointer :: initializer - integer, intent(in) :: order !! order of accuracy - integer, intent(in) :: cells !! number of grid cells spanning the domain - double precision, intent(in) :: x_min !! grid location minimum - double precision, intent(in) :: x_max !! grid location maximum - type(vector_1D_t) vector_1D - - call_julienne_assert(x_max .greaterThan. x_min) - call_julienne_assert(cells .isAtLeast. 2*order+1) - - associate(values => initializer(faces_1D(x_min, x_max, cells))) - vector_1D%tensor_1D_t = tensor_1D_t(values, x_min, x_max, cells, order) - end associate - vector_1D%divergence_operator_1D_ = divergence_operator_1D_t(k=order, dx=(x_max - x_min)/cells, cells=cells) - end function - -#endif - module procedure construct_1D_vector_from_parent vector_1D%tensor_1D_t = tensor_1D vector_1D%divergence_operator_1D_ = divergence_operator_1D_t(k=tensor_1D%order_, dx=(tensor_1D%x_max_ - tensor_1D%x_min_)/tensor_1D%cells_, cells=tensor_1D%cells_) diff --git a/test/divergence_operator_1D_test_m.F90 b/test/divergence_operator_1D_test_m.F90 index 9a37857..5ecd49b 100644 --- a/test/divergence_operator_1D_test_m.F90 +++ b/test/divergence_operator_1D_test_m.F90 @@ -20,9 +20,6 @@ module divergence_operator_1D_test_m ,test_result_t & ,usher use formal_m, only : vector_1D_t, vector_1D_initializer_i, scalar_1D_t, scalar_1D_initializer_i -#ifdef __GFORTRAN__ - use formal_m, only : divergence_1D_t -#endif implicit none type, extends(test_t) :: divergence_operator_1D_test_t @@ -70,40 +67,22 @@ function check_2nd_order_div_grad_parabola() result(test_diagnosis) type(test_diagnosis_t) test_diagnosis procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => parabola double precision, parameter :: expected_divergence = 1D0 -#ifdef __GFORTRAN__ - type(divergence_1D_t) div_grad_scalar - div_grad_scalar = .div. (.grad. scalar_1D_t(scalar_1D_initializer, order=2, cells=16, x_min=0D0, x_max=5D0)) -#else associate(div_grad_scalar => .div. (.grad. scalar_1D_t(scalar_1D_initializer, order=2, cells=16, x_min=0D0, x_max=5D0))) -#endif - test_diagnosis = passing_test() test_diagnosis = test_diagnosis .also. (.all. (div_grad_scalar%values() .approximates. expected_divergence .within. tight_tolerance)) & // " (2nd-order .div. (.grad. (x**2)/2))" - -#ifndef __GFORTRAN__ end associate -#endif end function function check_4th_order_div_grad_parabola() result(test_diagnosis) type(test_diagnosis_t) test_diagnosis procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => parabola double precision, parameter :: expected_divergence = 1D0 -#ifdef __GFORTRAN__ - type(divergence_1D_t) div_grad_scalar - div_grad_scalar = .div. (.grad. scalar_1D_t(scalar_1D_initializer, order=4, cells=16, x_min=0D0, x_max=9D0)) -#else associate(div_grad_scalar => .div. (.grad. scalar_1D_t(scalar_1D_initializer, order=4, cells=16, x_min=0D0, x_max=9D0))) -#endif - test_diagnosis = passing_test() test_diagnosis = test_diagnosis .also. (.all. (div_grad_scalar%values() .approximates. expected_divergence .within. tight_tolerance)) & // " (4th-order .div. (.grad. (x**2)/2))" - -#ifndef __GFORTRAN__ end associate -#endif end function pure function sinusoid(x) result(y) @@ -117,16 +96,10 @@ function check_2nd_order_div_sinusoid_convergence() result(test_diagnosis) procedure(vector_1D_initializer_i), pointer :: vector_1D_initializer => sinusoid double precision, parameter :: pi = 3.141592653589793D0 integer, parameter :: order_desired = 2, coarse_cells=100, fine_cells=coarse_cells+1 -#ifdef __GFORTRAN__ - type(divergence_1D_t) div_coarse, div_fine - div_coarse = .div. vector_1D_t(vector_1D_initializer , order=order_desired, cells=coarse_cells, x_min=0D0, x_max=2*pi) - div_fine = .div. vector_1D_t(vector_1D_initializer , order=order_desired, cells=fine_cells , x_min=0D0, x_max=2*pi) -#else associate( & div_coarse => .div. vector_1D_t(vector_1D_initializer , order=order_desired, cells=coarse_cells, x_min=0D0, x_max=2*pi) & ,div_fine => .div. vector_1D_t(vector_1D_initializer , order=order_desired, cells=fine_cells , x_min=0D0, x_max=2*pi) & ) -#endif test_diagnosis = passing_test() test_diagnosis = test_diagnosis .also. (size(div_coarse%values()) .equalsExpected. coarse_cells) test_diagnosis = test_diagnosis .also. (size(div_fine%values()) .equalsExpected. fine_cells) @@ -154,9 +127,7 @@ function check_2nd_order_div_sinusoid_convergence() result(test_diagnosis) end associate end associate end associate -#ifndef __GFORTRAN__ end associate -#endif end function function check_4th_order_div_sinusoid_convergence() result(test_diagnosis) @@ -164,16 +135,10 @@ function check_4th_order_div_sinusoid_convergence() result(test_diagnosis) procedure(vector_1D_initializer_i), pointer :: vector_1D_initializer => sinusoid double precision, parameter :: pi = 3.141592653589793D0 integer, parameter :: order_desired = 4, coarse_cells=500, fine_cells=coarse_cells+1 -#ifdef __GFORTRAN__ - type(divergence_1D_t) div_coarse, div_fine - div_coarse = .div. vector_1D_t(vector_1D_initializer , order=order_desired, cells=coarse_cells, x_min=0D0, x_max=2*pi) - div_fine = .div. vector_1D_t(vector_1D_initializer , order=order_desired, cells=fine_cells , x_min=0D0, x_max=2*pi) -#else associate( & div_coarse => .div. vector_1D_t(vector_1D_initializer , order=order_desired, cells=coarse_cells, x_min=0D0, x_max=2*pi) & ,div_fine => .div. vector_1D_t(vector_1D_initializer , order=order_desired, cells=fine_cells , x_min=0D0, x_max=2*pi) & ) -#endif associate( & x_coarse => div_coarse%grid() & ,x_fine => div_fine%grid() & @@ -204,9 +169,7 @@ function check_4th_order_div_sinusoid_convergence() result(test_diagnosis) end associate end associate end associate -#ifndef __GFORTRAN__ end associate -#endif end function end module divergence_operator_1D_test_m \ No newline at end of file diff --git a/test/gradient_operator_1D_test_m.F90 b/test/gradient_operator_1D_test_m.F90 index 1f35875..e9564ad 100644 --- a/test/gradient_operator_1D_test_m.F90 +++ b/test/gradient_operator_1D_test_m.F90 @@ -19,9 +19,6 @@ module gradient_operator_1D_test_m ,test_result_t & ,usher use formal_m, only : scalar_1D_t, scalar_1D_initializer_i -#ifdef __GFORTRAN__ - use formal_m, only : vector_1D_t, vector_1D_initializer_i, gradient_1D_t -#endif type, extends(test_t) :: gradient_operator_1D_test_t contains @@ -69,34 +66,16 @@ function check_grad_const() result(test_diagnosis) double precision, parameter :: grad_expected = 0. procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => const -#ifdef __GFORTRAN__ - type(gradient_1D_t) grad - - grad = .grad. scalar_1D_t(scalar_1D_initializer, order=2, cells=16, x_min=0D0, x_max=4D0) -#else associate(grad => .grad. scalar_1D_t(scalar_1D_initializer, order=2, cells=16, x_min=0D0, x_max=4D0)) -#endif - test_diagnosis = passing_test() test_diagnosis = test_diagnosis .also. (.all. (grad%values() .approximates. grad_expected .within. loose_tolerance)) & // " (2nd-order .grad.(5))" - -#ifndef __GFORTRAN__ end associate -#endif -#ifdef __GFORTRAN__ - grad = .grad. scalar_1D_t(scalar_1D_initializer, order=4, cells=16, x_min=0D0, x_max=8D0) -#else associate(grad => .grad. scalar_1D_t(scalar_1D_initializer, order=4, cells=16, x_min=0D0, x_max=8D0)) -#endif - test_diagnosis = test_diagnosis .also. (.all. (grad%values() .approximates. grad_expected .within. loose_tolerance)) & // " (4th-order .grad.(5))" - -#ifndef __GFORTRAN__ end associate -#endif end function pure function line(x) result(y) @@ -109,34 +88,17 @@ function check_grad_line() result(test_diagnosis) type(test_diagnosis_t) test_diagnosis double precision, parameter :: grad_expected = 14D0 procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => line -#ifdef __GFORTRAN__ - type(gradient_1D_t) grad - grad = .grad. scalar_1D_t(scalar_1D_initializer, order=2, cells=16, x_min=0D0, x_max=4D0) -#else associate(grad => .grad. scalar_1D_t(scalar_1D_initializer, order=2, cells=16, x_min=0D0, x_max=4D0)) -#endif - test_diagnosis = passing_test() test_diagnosis = test_diagnosis .also. (.all. (grad%values() .approximates. grad_expected .within. loose_tolerance)) & // " (2nd-order .grad.(14*x + 3))" - -#ifndef __GFORTRAN__ end associate -#endif -#ifdef __GFORTRAN__ - grad = .grad. scalar_1D_t(scalar_1D_initializer, order=4, cells=16, x_min=0D0, x_max=8D0) -#else associate(grad => .grad. scalar_1D_t(scalar_1D_initializer, order=4, cells=16, x_min=0D0, x_max=8D0)) -#endif - test_diagnosis = test_diagnosis .also. (.all. (grad%values() .approximates. grad_expected .within. loose_tolerance)) & // " (4th-order .grad.(14*x + 3))" - -#ifndef __GFORTRAN__ end associate -#endif end function pure function parabola(x) result(y) @@ -148,43 +110,25 @@ pure function parabola(x) result(y) function check_grad_parabola() result(test_diagnosis) type(test_diagnosis_t) test_diagnosis procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => parabola -#ifdef __GFORTRAN__ - type(gradient_1D_t) grad - grad = .grad. scalar_1D_t(scalar_1D_initializer , order=2, cells=16, x_min=0D0, x_max=4D0) -#else associate(grad => .grad. scalar_1D_t(scalar_1D_initializer , order=2, cells=16, x_min=0D0, x_max=4D0)) -#endif - test_diagnosis = passing_test() - associate(x => grad%grid()) associate(grad_expected => 14*x + 3) test_diagnosis = test_diagnosis .also. (.all. (grad%values() .approximates. grad_expected .within. loose_tolerance)) & // " (2nd-order .grad.(7*x**2 + 3*x + 5))" end associate end associate - -#ifndef __GFORTRAN__ end associate -#endif -#ifdef __GFORTRAN__ - grad = .grad. scalar_1D_t(scalar_1D_initializer , order=4, cells=16, x_min=0D0, x_max=8D0) -#else associate(grad => .grad. scalar_1D_t(scalar_1D_initializer , order=4, cells=16, x_min=0D0, x_max=8D0)) -#endif - associate(x => grad%grid()) associate(grad_expected => 14*x + 3) test_diagnosis = test_diagnosis .also. (.all. (grad%values() .approximates. grad_expected .within. loose_tolerance)) & // " (4th-order .grad.(7*x**2 + 3*x + 5))" end associate end associate - -#ifndef __GFORTRAN__ end associate -#endif end function pure function sinusoid(x) result(y) @@ -199,17 +143,10 @@ function check_2nd_order_grad_convergence() result(test_diagnosis) procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => sinusoid double precision, parameter :: pi = 3.141592653589793D0 integer, parameter :: order_desired = 2, coarse_cells=200, fine_cells=coarse_cells+1 -#ifdef __GFORTRAN__ - type(gradient_1D_t) grad_coarse, grad_fine - - grad_coarse = .grad. scalar_1D_t(scalar_1D_initializer , order=order_desired, cells=coarse_cells, x_min=0D0, x_max=2*pi) - grad_fine = .grad. scalar_1D_t(scalar_1D_initializer , order=order_desired, cells=fine_cells , x_min=0D0, x_max=2*pi) -#else associate( & grad_coarse => .grad. scalar_1D_t(scalar_1D_initializer , order=order_desired, cells=coarse_cells, x_min=0D0, x_max=2*pi) & ,grad_fine => .grad. scalar_1D_t(scalar_1D_initializer , order=order_desired, cells=fine_cells , x_min=0D0, x_max=2*pi) & ) -#endif associate( & x_coarse => grad_coarse%grid() & ,x_fine => grad_fine%grid() & @@ -236,28 +173,18 @@ function check_2nd_order_grad_convergence() result(test_diagnosis) end associate end associate end associate -#ifndef __GFORTRAN__ end associate -#endif end function function check_4th_order_grad_convergence() result(test_diagnosis) type(test_diagnosis_t) test_diagnosis procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => sinusoid double precision, parameter :: pi = 3.141592653589793D0 -#ifdef __GFORTRAN__ - integer, parameter :: order_desired = 4, coarse_cells=300, fine_cells=coarse_cells+1 - type(gradient_1D_t) grad_coarse, grad_fine - - grad_coarse = .grad. scalar_1D_t(scalar_1D_initializer , order=order_desired, cells=coarse_cells, x_min=0D0, x_max=2*pi) - grad_fine = .grad. scalar_1D_t(scalar_1D_initializer , order=order_desired, cells=fine_cells , x_min=0D0, x_max=2*pi) -#else integer, parameter :: order_desired = 4, coarse_cells=400, fine_cells=coarse_cells+1 associate( & grad_coarse => .grad. scalar_1D_t(scalar_1D_initializer , order=order_desired, cells=coarse_cells, x_min=0D0, x_max=2*pi) & ,grad_fine => .grad. scalar_1D_t(scalar_1D_initializer , order=order_desired, cells=fine_cells , x_min=0D0, x_max=2*pi) & ) -#endif associate( & x_coarse => grad_coarse%grid() & ,x_fine => grad_fine%grid() & @@ -284,9 +211,7 @@ function check_4th_order_grad_convergence() result(test_diagnosis) end associate end associate end associate -#ifndef __GFORTRAN__ end associate -#endif end function end module \ No newline at end of file diff --git a/test/laplacian_operator_1D_test_m.F90 b/test/laplacian_operator_1D_test_m.F90 index ccd5486..578065a 100644 --- a/test/laplacian_operator_1D_test_m.F90 +++ b/test/laplacian_operator_1D_test_m.F90 @@ -18,9 +18,6 @@ module laplacian_operator_1D_test_m ,test_result_t & ,usher use formal_m, only : scalar_1D_t, scalar_1D_initializer_i -#ifdef __GFORTRAN__ - use formal_m, only : laplacian_1D_t -#endif implicit none type, extends(test_t) :: laplacian_operator_1D_test_t @@ -68,20 +65,12 @@ function check_2nd_order_laplacian_parabola() result(test_diagnosis) type(test_diagnosis_t) test_diagnosis procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => parabola double precision, parameter :: expected_laplacian = 1D0 -#ifdef __GFORTRAN__ - type(laplacian_1D_t) laplacian_scalar - laplacian_scalar = .laplacian. scalar_1D_t(scalar_1D_initializer, order=2, cells=16, x_min=0D0, x_max=5D0) -#else - associate(laplacian_scalar => .laplacian. scalar_1D_t(scalar_1D_initializer, order=2, cells=16, x_min=0D0, x_max=5D0)) -#endif + associate(laplacian_scalar => .laplacian. scalar_1D_t(scalar_1D_initializer, order=2, cells=16, x_min=0D0, x_max=5D0)) test_diagnosis = passing_test() test_diagnosis = test_diagnosis .also. (.all. (laplacian_scalar%values() .approximates. expected_laplacian .within. tight_tolerance)) & // " (2nd-order .laplacian. [(x**2)/2]" - -#ifndef __GFORTRAN__ end associate -#endif end function pure function quartic(x) result(y) @@ -94,12 +83,7 @@ function check_4th_order_laplacian_of_quartic() result(test_diagnosis) type(test_diagnosis_t) test_diagnosis procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => quartic -#ifndef __GFORTRAN__ associate(laplacian_quartic => .laplacian. scalar_1D_t(scalar_1D_initializer, order=4, cells=16, x_min=0D0, x_max=40D0)) -#else - type(laplacian_1D_t) laplacian_quartic - laplacian_quartic = .laplacian. scalar_1D_t(scalar_1D_initializer, order=4, cells=16, x_min=0D0, x_max=40D0) -#endif associate(x => laplacian_quartic%grid()) associate(expected_laplacian => x**2, actual_laplacian => laplacian_quartic%values()) test_diagnosis = passing_test() @@ -107,9 +91,7 @@ function check_4th_order_laplacian_of_quartic() result(test_diagnosis) // " (4th-order .laplacian. [(x**4)/24]" end associate end associate -#ifndef __GFORTRAN__ end associate -#endif end function pure function f(x) @@ -140,16 +122,10 @@ function check_laplacian_convergence(order_desired, coarse_cells, fine_cells) re double precision, parameter :: pi = 3.141592653589793D0 integer, intent(in) :: order_desired, coarse_cells, fine_cells -#ifndef __GFORTRAN__ associate( & laplacian_coarse => .laplacian. scalar_1D_t(scalar_1D_initializer , order=order_desired, cells=coarse_cells, x_min=0D0, x_max=2*pi) & ,laplacian_fine => .laplacian. scalar_1D_t(scalar_1D_initializer , order=order_desired, cells=fine_cells , x_min=0D0, x_max=2*pi) & ) -#else - type(laplacian_1D_t) laplacian_coarse, laplacian_fine - laplacian_coarse = .laplacian. scalar_1D_t(scalar_1D_initializer , order=order_desired, cells=coarse_cells, x_min=0D0, x_max=2*pi) - laplacian_fine = .laplacian. scalar_1D_t(scalar_1D_initializer , order=order_desired, cells=fine_cells , x_min=0D0, x_max=2*pi) -#endif grids: & associate( & x_coarse => laplacian_coarse%grid() & @@ -204,9 +180,7 @@ function check_laplacian_convergence(order_desired, coarse_cells, fine_cells) re end associate laplacian_values end associate grids -#ifndef __GFORTRAN__ end associate -#endif end function end module \ No newline at end of file diff --git a/test/scalar_1D_test_m.F90 b/test/scalar_1D_test_m.F90 index 4b8edef..42c6393 100644 --- a/test/scalar_1D_test_m.F90 +++ b/test/scalar_1D_test_m.F90 @@ -63,21 +63,11 @@ function check_exponentiation() result(test_diagnosis) do order = 2, 4, 2 associate(scalar_1D => scalar_1D_t(scalar_1D_initializer, order=order, cells=10, x_min=0D0, x_max=10D0) ) -#ifndef __GFORTRAN__ associate(cube => scalar_1D**3 ) test_diagnosis = test_diagnosis .also. .all. & (cube%values() .approximates. scalar_1D%values()**3 .within. tolerance) & // string_t(" for order ") // string_t(order) end associate -#else - block - type(scalar_1D_t) cube - cube = scalar_1D**3 - test_diagnosis = test_diagnosis .also. .all. & - (cube%values() .approximates. scalar_1D%values()**3 .within. tolerance) & - // string_t(" for order ") // string_t(order) - end block -#endif end associate end do @@ -93,19 +83,10 @@ function check_divison_operator() result(test_diagnosis) do order = 2, 4, 2 associate(scalar_1D => scalar_1D_t(scalar_1D_initializer, order=order, cells=10, x_min=0D0, x_max=10D0) ) -#ifndef __GFORTRAN__ associate( half => scalar_1D/2 ) test_diagnosis = test_diagnosis .also. .all. (half%values() .approximates. scalar_1D%values()/2 .within. tolerance) & // string_t(" for order ") // string_t(order) end associate -#else - block - type(scalar_1D_t) half - half = scalar_1D/2 - test_diagnosis = test_diagnosis .also. .all. (half%values() .approximates. scalar_1D%values()/2 .within. tolerance) & - // string_t(" for order ") // string_t(order) - end block -#endif end associate end do From 40e4da07465cdc5aa4171583ba74dbc23229db91 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Mon, 25 May 2026 22:56:18 -0700 Subject: [PATCH 27/36] feat(vector_2D_t): add and test operator(.div.) --- src/formal/vector_2D_s.F90 | 24 +++++++++++------ test/vector_2D_test_m.F90 | 53 ++++++++++++++++++++++++++++++-------- 2 files changed, 58 insertions(+), 19 deletions(-) diff --git a/src/formal/vector_2D_s.F90 b/src/formal/vector_2D_s.F90 index f4ad4da..2c4fc02 100644 --- a/src/formal/vector_2D_s.F90 +++ b/src/formal/vector_2D_s.F90 @@ -15,6 +15,8 @@ use tensors_1D_m, only : faces_1D, vector_1D_t implicit none + integer, parameter :: x_dir=1, y_dir=2 + contains module procedure construct_2D_vector_from_function @@ -73,6 +75,7 @@ end procedure module procedure vector_2D_values + call_julienne_assert(allocated(self%values_)) vector_values = self%values_(:,:,:,1,1,1) end procedure @@ -80,7 +83,9 @@ type(string_t), allocatable :: lines(:) integer i, j, l - associate(x => self%grid(1), y => self%grid(2), header => [string_t("x,y,vector_x,vector_y")]) + call_julienne_assert(allocated(self%values_)) + + associate(x => self%grid(x_dir), y => self%grid(y_dir), header => [string_t("x,y,vector_x,vector_y")]) associate(num_blank_lines => size(y)-1) allocate(lines(size(header) + size(self%values_)/space_dimension + num_blank_lines)) end associate @@ -115,24 +120,27 @@ module procedure vector_2D_divergence - integer c, i, j + call_julienne_assert(allocated(self%values_)) divergence_2D%x_min_ = self%x_min_ divergence_2D%x_max_ = self%x_max_ divergence_2D%cells_ = self%cells_ divergence_2D%order_ = self%order_ - allocate(divergence_2D%values_(self%cells_(1)+2, self%cells_(2)+2, 1, 1, 1, 1)) + allocate(divergence_2D%values_(self%cells_(x_dir), self%cells_(y_dir), 1, 1, 1, 1)) divergence_x_term: & - do concurrent(integer :: j=1:size(divergence_2D%values_,2)) default(none) shared(divergence_2D, self) - divergence_2D%values_(:,j,1,1,1,1) = self%divergence_operator_1D_(1) .x. self%values_(:,j,1,1,1,1) + do concurrent(integer :: j=1:size(divergence_2D%values_,y_dir)) default(none) shared(divergence_2D, self) + associate(padded_divergence => self%divergence_operator_1D_(x_dir) .x. self%values_(:,j,x_dir,1,1,1)) + divergence_2D%values_(:,j,1,1,1,1) = padded_divergence(2:size(padded_divergence)-1) + end associate end do divergence_x_term add_y_term: & - do concurrent(integer :: i=1:size(divergence_2D%values_,1)) default(none) shared(divergence_2D, self) - divergence_2D%values_(i,:,2,1,1,1) = divergence_2D%values_(i,:,1,1,1,1) + & - (self%divergence_operator_1D_(2) .x. self%values_(i,:,1,1,1,1)) + do concurrent(integer :: i=1:size(divergence_2D%values_,x_dir)) default(none) shared(divergence_2D, self) + associate(padded_divergence => self%divergence_operator_1D_(y_dir) .x. self%values_(i,:,y_dir,1,1,1)) + divergence_2D%values_(i,:,1,1,1,1) = divergence_2D%values_(i,:,1,1,1,1) + padded_divergence(2:size(padded_divergence)-1) + end associate end do add_y_term end procedure diff --git a/test/vector_2D_test_m.F90 b/test/vector_2D_test_m.F90 index 0131a0e..ed145f9 100644 --- a/test/vector_2D_test_m.F90 +++ b/test/vector_2D_test_m.F90 @@ -8,6 +8,7 @@ module vector_2D_test_m ,operator(.also.) & ,operator(.approximates.) & ,operator(.within.) & + ,operator(.withinPercentage.) & ,passing_test & ,string_t & ,test_description_t & @@ -29,8 +30,8 @@ module vector_2D_test_m procedure, nopass :: results end type - double precision, parameter :: tolerance = 5D-2 integer, parameter :: space_dimension = 2 + double precision, parameter :: tolerance = 1D-2 contains @@ -48,19 +49,41 @@ function results() result(test_results) ]) end function - pure function v(x,y) result(z) + pure function biquadratic(x,y) result(z) double precision, intent(in) :: x(:), y(:) double precision z(size(x),size(y),space_dimension) do concurrent(integer :: i=1:size(x), j=1:size(y)) default(none) shared(x,y,z) - z(i,j,:) = [0D0,0D0] + z(i,j,:) = [ & + 1 - 2*x(i) + 3*x(i)**2 - x(i)*y(j)/5 + 3*y(j)**2 - 2*y(j) & + ,1 - 2*x(i) + 3*x(i)**2 - x(i)*y(j)/5 + 3*y(j)**2 - 2*y(j) & + ] end do end function - pure function div_v(x,y) result(divergence) + pure function biquadratic_divergence(x,y) result(divergence) double precision, intent(in) :: x(:), y(:) double precision divergence(size(x),size(y)) do concurrent(integer :: i=1:size(x), j=1:size(y)) default(none) shared(divergence,x,y) - divergence(i,j) = 0D0 + divergence(i,j) = (-2 + 6*x(i) - y(j)/5) + (-2 + 6*y(j) - x(i)/5) + end do + end function + + pure function cubic(x,y) result(z) + double precision, intent(in) :: x(:), y(:) + double precision z(size(x),size(y),space_dimension) + do concurrent(integer :: i=1:size(x), j=1:size(y)) default(none) shared(x,y,z) + z(i,j,:) = [ & + 1 - 2*x(i) + 3*x(i)**3 - x(i)*y(j)/5 + 3*y(j)**3 - 2*y(j) & + ,1 - 2*x(i) + 3*x(i)**3 - x(i)*y(j)/5 + 3*y(j)**3 - 2*y(j) & + ] + end do + end function + + pure function cubic_divergence(x,y) result(divergence) + double precision, intent(in) :: x(:), y(:) + double precision divergence(size(x),size(y)) + do concurrent(integer :: i=1:size(x), j=1:size(y)) default(none) shared(divergence,x,y) + divergence(i,j) = (-2 + 9*x(i)**2 - y(j)/5) + (-2 + 9*y(j)**2 - x(i)/5) end do end function @@ -70,17 +93,25 @@ function check_divergence() result(test_diagnosis) procedure(divergence_2D_initializer_i), pointer :: expected_divergence_initializer integer order - vector_2D_initializer => v - expected_divergence_initializer => div_v test_diagnosis = passing_test() do order = 2, 4, 2 - associate(vector_2D => vector_2D_t(vector_2D_initializer, order=order, cells=[30,20], x_min=[-1D0,1D0], x_max=[9D0,4D0])) + select case(order) + case(2) + vector_2D_initializer => biquadratic + expected_divergence_initializer => biquadratic_divergence + case(4) + vector_2D_initializer => cubic + expected_divergence_initializer => cubic_divergence + case default + error stop "check_divergence(vector_2D_test_m): unsupported order" + end select + associate(vector_2D => vector_2D_t(vector_2D_initializer, order=order, cells=[40,30], x_min=[0D0,0D0], x_max=[2D0,1D0])) associate(div_vector => .div. vector_2D) associate(expected_divergence => divergence_2D_t(expected_divergence_initializer, mold=vector_2D)) - test_diagnosis = test_diagnosis .also. & - (.all. (div_vector%values() .approximates. expected_divergence%values() .within. tolerance)) & - // string_t(" for order ") // string_t(order) + test_diagnosis = test_diagnosis .also. & + (.all. (div_vector%values() .approximates. expected_divergence%values() .within. tolerance)) & + // string_t(" for order ") // string_t(order) end associate end associate end associate From 9127b8e50f5a1990800e269c01ed7e46b63c5ef9 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Mon, 25 May 2026 23:47:08 -0700 Subject: [PATCH 28/36] feat(vector_3D_t): add and test operator(.div.) --- src/formal/divergence_3D_s.F90 | 89 +++++++++++++++++++++++ src/formal/tensors_3D_m.F90 | 72 ++++++++++++++++++ src/formal/vector_3D_s.F90 | 39 ++++++++++ src/formal_m.f90 | 4 +- test/driver.f90 | 2 + test/vector_3D_test_m.F90 | 129 +++++++++++++++++++++++++++++++++ 6 files changed, 334 insertions(+), 1 deletion(-) create mode 100644 src/formal/divergence_3D_s.F90 create mode 100644 test/vector_3D_test_m.F90 diff --git a/src/formal/divergence_3D_s.F90 b/src/formal/divergence_3D_s.F90 new file mode 100644 index 0000000..ed47b32 --- /dev/null +++ b/src/formal/divergence_3D_s.F90 @@ -0,0 +1,89 @@ +! 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_3D_m) divergence_3D_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, divergence_1D_t, cell_centers_1D + use julienne_m, only : string_t, operator(.csv.) + implicit none + +contains + + module procedure divergence_3D_values + divergence_values = self%values_(:,:,:,1,1,1,1) + end procedure + + module procedure divergence_3D_grid + associate(divergence_1D => divergence_1D_t( & + constant = 0D0 & + ,cells = self%cells_(direction) & + ,x_min = self%x_min_(direction) & + ,x_max = self%x_max_(direction) & + ,order = self%order_ & + )) + divergence_grid_1D = divergence_1D%grid() + end associate + end procedure + + module procedure construct_3D_divergence_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_1D(x_min(1), x_max(1), cells(1)) & + ,y => cell_centers_1D(x_min(2), x_max(2), cells(2)) & + ,z => cell_centers_1D(x_min(3), x_max(3), cells(3)) & + ) + divergence_3D%tensor_3D_t = tensor_3D_t( & + values = reshape(initializer(x,y,z), shape=[size(x),size(y),size(z),1,1,1,1]) & + ,cells = cells , x_min = x_min, x_max = x_max, order = order & + ) + end associate + end procedure + + module procedure construct_3D_divergence_from_vector_mold + divergence_3D = divergence_3D_t(initializer, cells = mold%cells_, x_min = mold%x_min_, x_max = mold%x_max_, order = mold%order_) + end procedure + + module procedure divergence_3D_to_file + type(string_t), allocatable :: lines(:) + integer i, j, k, l + + associate( & + x => self%grid(1) & + ,y => self%grid(2) & + ,z => self%grid(3) & + ,header => [string_t("x,y,z,divergence")] & + ) + associate(num_blank_lines => size(y)*size(z)-1) + allocate(lines(size(header) + size(self%values_) + num_blank_lines)) + end associate + lines(1:size(header)) = header + l = size(header) + do k = 1, size(z) + do j = 1, size(y) + do i = 1, size(x) + l = l + 1 + lines(l) = .csv. string_t([x(i), y(j), z(k), self%values_(i,j,k,1,1,1,1)]) + end do + if (j/=size(y) .or. k/=size(z)) then + l = l + 1 + lines(l) = "" + end if + end do + end do + end associate + + file = file_t(lines) + end procedure + +end submodule divergence_3D_s \ No newline at end of file diff --git a/src/formal/tensors_3D_m.F90 b/src/formal/tensors_3D_m.F90 index 7b919a5..62e9779 100644 --- a/src/formal/tensors_3D_m.F90 +++ b/src/formal/tensors_3D_m.F90 @@ -15,8 +15,10 @@ module tensors_3D_m public :: scalar_3D_t public :: vector_3D_t public :: gradient_3D_t + public :: divergence_3D_t public :: scalar_3D_initializer_i public :: vector_3D_initializer_i + public :: divergence_3D_initializer_i integer, parameter :: space_dimension = 3 @@ -29,6 +31,13 @@ pure function scalar_3D_initializer_i(x,y,z) result(f) double precision f(size(x),size(y),size(z)) end function + pure function divergence_3D_initializer_i(x,y,z) result(f) + !! Sampling function for initializing a divergence_3D_t object + implicit none + double precision, intent(in) :: x(:), y(:), z(:) + double precision f(size(x),size(y),size(z)) + end function + pure function vector_3D_initializer_i(x,y,z) result(v) !! Sampling function for initializing a vector_3D_t object import space_dimension @@ -111,9 +120,11 @@ pure module function construct_3D_scalar_from_mold(initializer, mold) result(sca generic :: values => vector_3D_values generic :: to_file => vector_3D_to_file generic :: grid => vector_3D_grid + generic :: operator(.div.) => vector_3D_divergence procedure, non_overridable, private :: vector_3D_values procedure, non_overridable, private :: vector_3D_to_file procedure, non_overridable, private :: vector_3D_grid + procedure, non_overridable, private :: vector_3D_divergence end type interface vector_3D_t @@ -154,6 +165,40 @@ pure module function construct_3D_vector_from_scalar_mold(initializer, mold) res !! A 3D mimetic gradient vector field abstraction with a public method that produces corresponding numerical quadrature weights end type + type, extends(tensor_3D_t) :: divergence_3D_t + !! A 3D mimetic divergence field abstraction with a public method that produces corresponding numerical quadrature weights + contains + generic :: values => divergence_3D_values + generic :: grid => divergence_3D_grid + generic :: to_file => divergence_3D_to_file + procedure, private, non_overridable :: divergence_3D_values + procedure, private, non_overridable :: divergence_3D_grid + procedure, private, non_overridable :: divergence_3D_to_file + end type + + interface divergence_3D_t + + pure module function construct_3D_divergence_from_function(initializer, order, cells, x_min, x_max) result(divergence_3D) + !! Result is a 3D divergence initialized by sampling the initializer at cell centers defined by the other arguments + implicit none + procedure(scalar_3D_initializer_i), pointer, intent(in) :: initializer + integer, intent(in) :: order !! order of accuracy + integer, intent(in) :: cells(:) !! number of grid cells spanning each spatial direction + double precision, intent(in) :: x_min(:) !! grid location minima + double precision, intent(in) :: x_max(:) !! grid location maxima + type(divergence_3D_t) divergence_3D + end function + + pure module function construct_3D_divergence_from_vector_mold(initializer, mold) result(divergence_3D) + !! Result is a 3D divergence initialized by sampling the initializer on cell centers defined by the mold + implicit none + procedure(divergence_3D_initializer_i), pointer, intent(in) :: initializer + type(vector_3D_t), intent(in) :: mold + type(divergence_3D_t) divergence_3D + end function + + end interface + interface pure module function scalar_3D_values(self) result(scalar_values) @@ -189,6 +234,13 @@ pure module function scalar_3D_gradient(self) result(gradient_3D) type(gradient_3D_t) gradient_3D end function + pure module function vector_3D_divergence(self) result(divergence_3D) + !! Result is mimetic divergence of the scalar_3D_t "self" + implicit none + class(vector_3D_t), intent(in) :: self + type(divergence_3D_t) divergence_3D + end function + pure module function scalar_3D_to_file(self) result(file) !! Result is a file_t object containing the grid points and the corresponding scalar values implicit none @@ -203,6 +255,26 @@ pure module function vector_3D_to_file(self) result(file) type(file_t) file end function + pure module function divergence_3D_grid(self, direction) result(divergence_grid_1D) + !! Result array contains divergence grid locations along the requested spatial direction + class(divergence_3D_t), intent(in) :: self + integer, intent(in) :: direction + double precision, allocatable :: divergence_grid_1D(:) !! grid points along the requested coordinate direction + end function + + pure module function divergence_3D_values(self) result(divergence_values) + !! Vector values getter + class(divergence_3D_t), intent(in) :: self + double precision, allocatable :: divergence_values(:,:,:) + end function + + pure module function divergence_3D_to_file(self) result(file) + !! Result is a file_t object containing the grid points and the corresponding divergence values + implicit none + class(divergence_3D_t), intent(in) :: self + type(file_t) file + end function + end interface end module tensors_3D_m diff --git a/src/formal/vector_3D_s.F90 b/src/formal/vector_3D_s.F90 index 8191bbe..0369f3e 100644 --- a/src/formal/vector_3D_s.F90 +++ b/src/formal/vector_3D_s.F90 @@ -16,6 +16,8 @@ use differential_operators_1D_m, only : divergence_operator_1D_t implicit none + integer, parameter :: x_dir=1, y_dir=2, z_dir=3 + contains module procedure construct_3D_vector_from_function @@ -131,4 +133,41 @@ end associate end procedure + module procedure vector_3D_divergence + + call_julienne_assert(allocated(self%values_)) + + divergence_3D%x_min_ = self%x_min_ + divergence_3D%x_max_ = self%x_max_ + divergence_3D%cells_ = self%cells_ + divergence_3D%order_ = self%order_ + + allocate(divergence_3D%values_(self%cells_(x_dir), self%cells_(y_dir), self%cells_(z_dir), 1, 1, 1, 1)) + + divergence_x_term: & + do concurrent(integer :: j=1:size(divergence_3D%values_,y_dir),k=1:size(divergence_3D%values_,z_dir)) & + default(none) shared(divergence_3D, self) + associate(padded_divergence => self%divergence_operator_1D_(x_dir) .x. self%values_(:,j,k,x_dir,1,1,1)) + divergence_3D%values_(:,j,k,1,1,1,1) = padded_divergence(2:size(padded_divergence)-1) + end associate + end do divergence_x_term + + add_y_term: & + do concurrent(integer :: i=1:size(divergence_3D%values_,x_dir), k=1:size(divergence_3D%values_,z_dir)) & + default(none) shared(divergence_3D, self) + associate(padded_divergence => self%divergence_operator_1D_(y_dir) .x. self%values_(i,:,k,y_dir,1,1,1)) + divergence_3D%values_(i,:,k,1,1,1,1) = divergence_3D%values_(i,:,k,1,1,1,1) + padded_divergence(2:size(padded_divergence)-1) + end associate + end do add_y_term + + add_z_term: & + do concurrent(integer :: i=1:size(divergence_3D%values_,x_dir), j=1:size(divergence_3D%values_,y_dir)) & + default(none) shared(divergence_3D, self) + associate(padded_divergence => self%divergence_operator_1D_(z_dir) .x. self%values_(i,j,:,z_dir,1,1,1)) + divergence_3D%values_(i,j,:,1,1,1,1) = divergence_3D%values_(i,j,:,1,1,1,1) + padded_divergence(2:size(padded_divergence)-1) + end associate + end do add_z_term + + end procedure + end submodule vector_3D_s \ No newline at end of file diff --git a/src/formal_m.f90 b/src/formal_m.f90 index c6a07f3..0377f88 100644 --- a/src/formal_m.f90 +++ b/src/formal_m.f90 @@ -30,9 +30,11 @@ module formal_m use tensors_3D_m, only : & scalar_3D_t & ! discrete 3D scalar field derived type ,vector_3D_t & ! discrete 3D vector field derived type + ,divergence_3D_t & ! discrete 3D divergence field derived type ,gradient_3D_t & ! result of `.grad. s` for a scalar_3D_t s ,scalar_3D_initializer_i & ! scalar_3D_t initializer abstract interface - ,vector_3D_initializer_i ! vector_3D_t initializar abstract interface + ,vector_3D_initializer_i & ! vector_3D_t initializar abstract interface + ,divergence_3D_initializer_i ! divergence_2D_t initializar abstract interface use differential_operators_1D_m, only : & gradient_operator_1D_t & ! matrix operator defining a 1D mimetic gradient diff --git a/test/driver.f90 b/test/driver.f90 index 8cc6f35..35d7202 100644 --- a/test/driver.f90 +++ b/test/driver.f90 @@ -13,6 +13,7 @@ program test_suite_driver use scalar_2D_test_m, only : scalar_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 implicit none associate(test_harness => test_harness_t([ & @@ -25,6 +26,7 @@ program test_suite_driver ,test_fixture_t(scalar_2D_test_t()) & ,test_fixture_t(scalar_3D_test_t()) & ,test_fixture_t(vector_2D_test_t()) & + ,test_fixture_t(vector_3D_test_t()) & ])) call test_harness%report_results end associate diff --git a/test/vector_3D_test_m.F90 b/test/vector_3D_test_m.F90 new file mode 100644 index 0000000..231925a --- /dev/null +++ b/test/vector_3D_test_m.F90 @@ -0,0 +1,129 @@ +! Copyright (c) 2026, The Regents of the University of California +! Terms of use are as specified in LICENSE.txt + +module vector_3D_test_m + use julienne_m, only : & + operator(//) & + ,operator(.all.) & + ,operator(.also.) & + ,operator(.approximates.) & + ,operator(.within.) & + ,operator(.withinPercentage.) & + ,passing_test & + ,string_t & + ,test_description_t & + ,test_diagnosis_t & + ,test_result_t & + ,test_t & + ,usher + use formal_m, only : & + vector_3D_t & + ,vector_3D_initializer_i & + ,divergence_3D_t & + ,divergence_3D_initializer_i + + implicit none + + type, extends(test_t) :: vector_3D_test_t + contains + procedure, nopass :: subject + procedure, nopass :: results + end type + + integer, parameter :: space_dimension = 3 + double precision, parameter :: tolerance = 1D-2 + +contains + + pure function subject() result(test_subject) + character(len=:), allocatable :: test_subject + test_subject = 'The vector_3D_t derived type' + end function + + function results() result(test_results) + type(vector_3D_test_t) vector_3D_test + type(test_result_t), allocatable :: test_results(:) + + test_results = vector_3D_test%run([ & + test_description_t('computing the divergence of a vector field', usher(check_divergence)) & + ]) + end function + + pure function biquadratic(x,y,z) result(vector_field) + double precision, intent(in) :: x(:), y(:), z(:) + double precision vector_field(size(x),size(y),size(z),space_dimension) + do concurrent(integer :: i=1:size(x), j=1:size(y), k=1:size(z)) default(none) shared(x,y,z,vector_field) + vector_field(i,j,k,:) = [ & + 1 - 2*x(i) + 3*x(i)**2 - x(i)*y(j)/5 + 3*y(j)**2 - 2*y(j) & + ,1 - 2*x(i) + 3*x(i)**2 - x(i)*y(j)/5 + 3*y(j)**2 - 2*y(j) & + ,z(k) & + ] + end do + end function + + pure function biquadratic_divergence(x,y,z) result(divergence) + double precision, intent(in) :: x(:), y(:), z(:) + double precision divergence(size(x),size(y),size(z)) + do concurrent(integer :: i=1:size(x), j=1:size(y), k=1:size(z)) default(none) shared(x,y,z,divergence) + divergence(i,j,k) = & + (-2 + 6*x(i) - y(j)/5) & ! du/dx + + (-2 + 6*y(j) - x(i)/5) & ! dv/dy + + (1D0) ! dw/dz + end do + end function + + pure function cubic(x,y,z) result(vector_field) + double precision, intent(in) :: x(:), y(:), z(:) + double precision vector_field(size(x),size(y),size(z),space_dimension) + do concurrent(integer :: i=1:size(x), j=1:size(y), k=1:size(z)) default(none) shared(x,y,z,vector_field) + vector_field(i,j,k,:) = [ & + 1 - 2*x(i) + 3*x(i)**3 - x(i)*y(j)/5 + 3*y(j)**3 - 2*y(j) & + ,1 - 2*x(i) + 3*x(i)**3 - x(i)*y(j)/5 + 3*y(j)**3 - 2*y(j) & + ,-z(k) & + ] + end do + end function + + pure function cubic_divergence(x,y,z) result(divergence) + double precision, intent(in) :: x(:), y(:), z(:) + double precision divergence(size(x),size(y),size(z)) + do concurrent(integer :: i=1:size(x), j=1:size(y), k=1:size(z)) default(none) shared(x,y,z,divergence) + divergence(i,j,k) = & + (-2 + 9*x(i)**2 - y(j)/5) & ! du/dx + + (-2 + 9*y(j)**2 - x(i)/5) & ! dv/dx + + (-1D0) ! dw/dz + end do + end function + + function check_divergence() result(test_diagnosis) + type(test_diagnosis_t) test_diagnosis + procedure(vector_3D_initializer_i), pointer :: vector_3D_initializer + procedure(divergence_3D_initializer_i), pointer :: expected_divergence_initializer + integer order + + test_diagnosis = passing_test() + + do order = 2, 4, 2 + select case(order) + case(2) + vector_3D_initializer => biquadratic + expected_divergence_initializer => biquadratic_divergence + case(4) + vector_3D_initializer => cubic + expected_divergence_initializer => cubic_divergence + case default + error stop "check_divergence(vector_3D_test_m): unsupported order" + end select + associate(vector_3D => vector_3D_t(vector_3D_initializer, order=order, cells=[40,20,30], x_min=[0D0,-.5D0,0D0], x_max=[2D0,0D0,3D0])) + associate(div_vector => .div. vector_3D) + associate(expected_divergence => divergence_3D_t(expected_divergence_initializer, mold=vector_3D)) + test_diagnosis = test_diagnosis .also. & + (.all. (div_vector%values() .approximates. expected_divergence%values() .within. tolerance)) & + // string_t(" for order ") // string_t(order) + end associate + end associate + end associate + end do + end function + +end module vector_3D_test_m \ No newline at end of file From 00a77ae73df0319f9b021735944b6f1a64822d11 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Tue, 26 May 2026 10:42:36 -0700 Subject: [PATCH 29/36] feat(2D-sink): add example divergence calculation --- example/2D-sink.F90 | 59 +++++++++++++++++++ ...ential.gnuplot => 2D-scalar-field.gnuplot} | 13 ++-- ...cities.gnuplot => 2D-vector-field.gnuplot} | 9 +-- 3 files changed, 72 insertions(+), 9 deletions(-) create mode 100644 example/2D-sink.F90 rename example/scripts/{velocity-potential.gnuplot => 2D-scalar-field.gnuplot} (67%) rename example/scripts/{velocities.gnuplot => 2D-vector-field.gnuplot} (74%) diff --git a/example/2D-sink.F90 b/example/2D-sink.F90 new file mode 100644 index 0000000..f6bbfcb --- /dev/null +++ b/example/2D-sink.F90 @@ -0,0 +1,59 @@ +! Copyright (c) 2026, The Regents of the University of California +! Terms of use are as specified in LICENSE.txt + +#include "julienne-assert-macros.h" + +module sink_2D_functions_m + use julienne_m, only : call_julienne_assert_ + implicit none + + integer, parameter :: space_dimension = 2 + double precision, parameter :: pi = acos(-1D0) + +contains + + pure function velocity(x,y) result(v) + double precision, intent(in) :: x(:), y(:) + double precision v(size(x),size(y),space_dimension), theta + double precision, parameter :: Q = 1D0 + do concurrent(integer :: i=1:size(x), j=1:size(y)) + associate(r => sqrt(x(i)**2 + y(j)**2)) + call_julienne_assert(r /= 0D0) + v(i,j,:) = -(Q/(2*pi))*[x(i), y(j)]/(x(i)**2 + y(j)**2) + end associate + end do + end function + + pure function divergence(x,y) result(div_v) + double precision, intent(in) :: x(:), y(:) + double precision div_v(size(x),size(y)) + call_julienne_assert(.not. any(x == 0D0 .and. y == 0D0)) + div_v = 0D0 + end function + +end module + +program sink_2D + use julienne_m, only : file_t + use sink_2D_functions_m, only : velocity, divergence, pi + use formal_m, only : vector_2D_t, divergence_2D_t, divergence_2D_initializer_i, vector_2D_initializer_i + implicit none + + integer, parameter :: order = 4 + procedure(vector_2D_initializer_i), pointer :: vector_2D_initializer + procedure(divergence_2D_initializer_i), pointer :: divergence_2D_initializer + + divergence_2D_initializer => divergence + vector_2D_initializer => velocity + + 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()) + 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") + end associate + end associate + end associate + +end program \ No newline at end of file diff --git a/example/scripts/velocity-potential.gnuplot b/example/scripts/2D-scalar-field.gnuplot similarity index 67% rename from example/scripts/velocity-potential.gnuplot rename to example/scripts/2D-scalar-field.gnuplot index 82668b0..b2e69b5 100644 --- a/example/scripts/velocity-potential.gnuplot +++ b/example/scripts/2D-scalar-field.gnuplot @@ -1,12 +1,15 @@ -# ============================================================ -# velocity-potential.gnuplot -- surface plot CSV +# ============================================================================ +# 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 velocity-potential.gnuplot -# ============================================================ +# Usage: gnuplot -d "base_name='velocity-potential'" 2D-scalar-field.gnuplot +# Default: base_name='velocity-potential' +# ============================================================================ + +if (!exists("base_name")) base_name = "velocity-potential" -base_name = "velocity-potential" datafile = base_name . ".csv" + set datafile separator "," # --- 1. Read column headers from line 1 --- diff --git a/example/scripts/velocities.gnuplot b/example/scripts/2D-vector-field.gnuplot similarity index 74% rename from example/scripts/velocities.gnuplot rename to example/scripts/2D-vector-field.gnuplot index fb062c1..1b09746 100644 --- a/example/scripts/velocities.gnuplot +++ b/example/scripts/2D-vector-field.gnuplot @@ -1,9 +1,10 @@ -# ============================================================= -# velocities.gnuplot -- 2D vector/quiver plot from a CSV +# =============================================================== +# vector-field.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 -# ============================================================= +# Usage: gnuplot -e "base_name='velocity'" vector-field.gnuplot +# Default: base_name='velocity' +# =============================================================== if (!exists("base_name")) base_name = "velocity" From 7d72f1c67674f312b2dc253de9bb568a9ee5cd80 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Tue, 26 May 2026 10:49:38 -0700 Subject: [PATCH 30/36] chore(2D-vortex): center the example output --- example/2D-vortex.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/example/2D-vortex.F90 b/example/2D-vortex.F90 index edfe4be..6aae958 100644 --- a/example/2D-vortex.F90 +++ b/example/2D-vortex.F90 @@ -40,7 +40,7 @@ program vortex_2D 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(phi => scalar_2D_t(scalar_2D_initializer, order=order, cells=[11,11], x_min=[-pi,-pi], x_max=[pi,pi])) associate( velocity => .grad. phi & ,expected_velocity => vector_2D_t(vector_2D_initializer, mold=phi) & ) From c9f44a87f16d171454f2f8230d2b2b6d28538c4d Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Tue, 26 May 2026 11:24:55 -0700 Subject: [PATCH 31/36] chore(div_{2D,3D}): rm unused use-association --- src/formal/divergence_2D_s.F90 | 2 +- src/formal/divergence_3D_s.F90 | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/formal/divergence_2D_s.F90 b/src/formal/divergence_2D_s.F90 index d5aafce..e1d5ac3 100644 --- a/src/formal/divergence_2D_s.F90 +++ b/src/formal/divergence_2D_s.F90 @@ -10,7 +10,7 @@ ,operator(.equalsExpected.) & ,operator(.greaterThan.) & ,operator(.isAtLeast.) - use tensors_1D_m, only : cell_centers_extended_1D, divergence_1D_t, cell_centers_1D + use tensors_1D_m, only : divergence_1D_t, cell_centers_1D use julienne_m, only : string_t, operator(.csv.) implicit none diff --git a/src/formal/divergence_3D_s.F90 b/src/formal/divergence_3D_s.F90 index ed47b32..fa03ff1 100644 --- a/src/formal/divergence_3D_s.F90 +++ b/src/formal/divergence_3D_s.F90 @@ -10,7 +10,7 @@ ,operator(.equalsExpected.) & ,operator(.greaterThan.) & ,operator(.isAtLeast.) - use tensors_1D_m, only : cell_centers_extended_1D, divergence_1D_t, cell_centers_1D + use tensors_1D_m, only : divergence_1D_t, cell_centers_1D use julienne_m, only : string_t, operator(.csv.) implicit none From 9df220c8774f94bfabefab1e6807fbc34fae47a9 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Tue, 26 May 2026 11:27:47 -0700 Subject: [PATCH 32/36] fix(tensors_1D_m): mk cell_centers_1D public --- src/formal/tensors_1D_m.F90 | 1 + 1 file changed, 1 insertion(+) diff --git a/src/formal/tensors_1D_m.F90 b/src/formal/tensors_1D_m.F90 index 67a5ec9..e7fd6aa 100644 --- a/src/formal/tensors_1D_m.F90 +++ b/src/formal/tensors_1D_m.F90 @@ -23,6 +23,7 @@ module tensors_1D_m public :: d_dx public :: d2_dx2 public :: cell_centers_extended_1D + public :: cell_centers_1D public :: faces_1D abstract interface From fefe657ebdb59485db1c40df6c8fc3f0f7513aa7 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Sun, 31 May 2026 21:06:33 -0700 Subject: [PATCH 33/36] doc(README): update lfortran version to "latest" --- README.md | 9 +-------- 1 file changed, 1 insertion(+), 8 deletions(-) diff --git a/README.md b/README.md index 59c30fd..1d55458 100644 --- a/README.md +++ b/README.md @@ -78,14 +78,7 @@ Building and testing LLVM | `flang` | 19 | `fpm test --compiler flang --profile release --flag "-mmlir -allow-assumed-rank"` ### `fpm` versions before 0.13.0 -With LLVM 20-22, replace the above `flang` command with -``` - fpm test --compiler flang --flag "-O3" -``` -With LLVM 19, replace the above `flang` command with -``` - fpm test --compiler flang --flag "-O3 -mmlir -allow-assumed-rank" -``` +With LLVM, replace the `flang` with `flang-new` and delete `--profile release` from the above commands. Unsupported Compilers --------------------- From 58556a83f6b67386e73cae84ea1f37cf63292fde Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Fri, 29 May 2026 15:00:14 -0500 Subject: [PATCH 34/36] chore(CI): rm --separate-compilation with lfortran --- .github/workflows/build.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml index 715152e..88c9b42 100644 --- a/.github/workflows/build.yml +++ b/.github/workflows/build.yml @@ -224,7 +224,7 @@ jobs: fi elif test "$FC" = "lfortran" ; then echo "FPM_FC=lfortran" >> "$GITHUB_ENV" - echo "FFLAGS=--cpp --separate-compilation --realloc-lhs-arrays $FFLAGS" >> "$GITHUB_ENV" + echo "FFLAGS=--cpp --realloc-lhs-arrays $FFLAGS" >> "$GITHUB_ENV" echo "FPM_FLAGS=--profile debug --verbose" >> "$GITHUB_ENV" ; : fpm 0.13 workaround else echo "FPM_FC=gfortran-${COMPILER_VERSION}" >> "$GITHUB_ENV" From 8ab3ddbd86dc6a07c78d3f8dfa98df847f4b9f3a Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Sat, 30 May 2026 18:36:16 -0500 Subject: [PATCH 35/36] feat(2D): gather consistency/conformability checks This commit adds type-bound consistency- and conformability- checking functions to the following types: - tensor_2D_t - scalar_2D_t - vector_2D_t - gradient_2D_t - divergence_2D_t where the tensor_2D_t functions are private and are inherited and used by the other listed types and where each of the other types has `consistent` and `conformable` generic bindings. Consistency functions assert that the values_ component is allocated and that the remaining components are consistent in size with the spatial dimension and consistent in value with each other, e.g., that the number of cells is sufficient for the order of accuracy. Calls to the new functions are distributed throughout the {tensor,scalar,vector,gradient,divergence}_2D_s submodules as effective pre- and/or post-conditions for each module procedure. --- src/formal/divergence_2D_s.F90 | 34 +++++++--- src/formal/scalar_2D_s.F90 | 40 ++++++++++-- src/formal/tensor_2D_s.F90 | 35 ++++++++++- src/formal/tensors_2D_m.F90 | 70 +++++++++++++++++++++ src/formal/vector_2D_s.F90 | 112 +++++++++++++++++++++++++-------- test/scalar_2D_test_m.F90 | 5 +- 6 files changed, 254 insertions(+), 42 deletions(-) diff --git a/src/formal/divergence_2D_s.F90 b/src/formal/divergence_2D_s.F90 index e1d5ac3..dfc1e83 100644 --- a/src/formal/divergence_2D_s.F90 +++ b/src/formal/divergence_2D_s.F90 @@ -17,10 +17,18 @@ contains module procedure divergence_2D_values + + call_julienne_assert(self%consistent()) + divergence_values = self%values_(:,:,1,1,1,1) + end procedure module procedure divergence_2D_grid + + call_julienne_assert(self%consistent()) + + construct_prototype: & associate(divergence_1D => divergence_1D_t( & constant = 0D0 & ,cells = self%cells_(direction) & @@ -29,31 +37,43 @@ ,order = self%order_ & )) divergence_grid_1D = divergence_1D%grid() - end associate + end associate construct_prototype end procedure module procedure construct_2D_divergence_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_1D(x_min(1), x_max(1), cells(1)), y => cell_centers_1D(x_min(2), x_max(2), cells(2))) + define_grid: & + associate( & + x => cell_centers_1D(x_min(1), x_max(1), cells(1)) & + ,y => cell_centers_1D(x_min(2), x_max(2), cells(2)) & + ) divergence_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 & ) - end associate + end associate define_grid + + call_julienne_assert(divergence_2D%consistent()) + end procedure module procedure construct_2D_divergence_from_vector_mold + + call_julienne_assert(mold%consistent()) + divergence_2D = divergence_2D_t(initializer, cells = mold%cells_, x_min = mold%x_min_, x_max = mold%x_max_, order = mold%order_) + + call_julienne_assert(divergence_2D%consistent()) + call_julienne_assert(divergence_2D%conformable(mold)) + end procedure module procedure divergence_2D_to_file type(string_t), allocatable :: lines(:) integer i, j, l + call_julienne_assert(self%consistent()) + associate(x => self%grid(1), y => self%grid(2), header => [string_t("x,y,divergence")]) associate(num_blank_lines => size(y)-1) allocate(lines(size(header) + size(self%values_) + num_blank_lines)) diff --git a/src/formal/scalar_2D_s.F90 b/src/formal/scalar_2D_s.F90 index e6d23c6..c958711 100644 --- a/src/formal/scalar_2D_s.F90 +++ b/src/formal/scalar_2D_s.F90 @@ -16,11 +16,27 @@ contains + module procedure scalar_2D_consistent + call_julienne_assert(self%tensor_2D_consistent()) + call_julienne_assert(.all. (self%cells_ .isAtLeast. 2*self%order_)) + call_julienne_assert(size(self%gradient_operator_1D_) .equalsExpected. space_dimension) + self_consistent = .true. + end procedure + + module procedure scalar_2D_conformable_scalar + call_julienne_assert(scalar_2D_consistent(self)) + call_julienne_assert(scalar_2D_consistent(scalar_2D)) + call_julienne_assert(self%tensor_2D_conformable(scalar_2D)) + conformable = .true. + end procedure + module procedure scalar_2D_values + call_julienne_assert(self%consistent()) scalar_values = self%values_(:,:,1,1,1,1) end procedure module procedure scalar_2D_grid + call_julienne_assert(self%consistent()) associate(scalar_1D => scalar_1D_t( & constant = 0D0 & ,cells = self%cells_(direction) & @@ -34,27 +50,35 @@ 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))) + define_grid: & + 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 associate define_grid + + call_julienne_assert(scalar_2D%consistent()) + end procedure module procedure construct_2D_scalar_from_mold + call_julienne_assert(mold%consistent()) scalar_2D = scalar_2D_t(initializer, cells = mold%cells_, x_min = mold%x_min_, x_max = mold%x_max_, order = mold%order_) + call_julienne_assert(scalar_2D%consistent()) end procedure module procedure scalar_2D_gradient integer c, i, j + call_julienne_assert(self%consistent()) + gradient_2D%x_min_ = self%x_min_ gradient_2D%x_max_ = self%x_max_ gradient_2D%cells_ = self%cells_ @@ -80,12 +104,16 @@ !end associate check_corbino_castillo_eq_17 end associate + call_julienne_assert(gradient_2D%consistent()) + end procedure module procedure scalar_2D_to_file type(string_t), allocatable :: lines(:) integer i, j, l + call_julienne_assert(self%consistent()) + 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)) diff --git a/src/formal/tensor_2D_s.F90 b/src/formal/tensor_2D_s.F90 index cf2a2df..dc65b40 100644 --- a/src/formal/tensor_2D_s.F90 +++ b/src/formal/tensor_2D_s.F90 @@ -1,4 +1,16 @@ +! 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) tensor_2D_s + use julienne_m, only : & + call_julienne_assert_ & + ,operator(.all.) & + ,operator(.approximates.) & + ,operator(.equalsExpected.) & + ,operator(.greaterThan.) & + ,operator(.within.) implicit none contains @@ -11,4 +23,25 @@ tensor_2D%order_ = order end procedure -end submodule \ No newline at end of file + module procedure tensor_2D_consistent + call_julienne_assert(allocated(self%values_)) + call_julienne_assert(.all. ([size(self%cells_), size(self%x_min_), size(self%x_max_)] .equalsExpected. space_dimension)) + call_julienne_assert(.all. (self%x_max_ .greaterThan. self%x_min_)) + self_consistent = .true. + end procedure + + module procedure tensor_2D_conformable + call_julienne_assert( tensor_2D_consistent(self) ) + call_julienne_assert( tensor_2D_consistent(tensor_2D) ) + call_julienne_assert(.all. (shape(self%cells_) .equalsExpected. shape(tensor_2D%cells_))) + call_julienne_assert(.all. (shape(self%order_) .equalsExpected. shape(tensor_2D%order_))) + call_julienne_assert(.all. (shape(self%x_min_) .equalsExpected. shape(tensor_2D%x_min_))) + call_julienne_assert(.all. (shape(self%x_max_) .equalsExpected. shape(tensor_2D%x_max_))) + call_julienne_assert(.all. (self%cells_ .equalsExpected. tensor_2D%cells_)) + call_julienne_assert(.all. (self%order_ .equalsExpected. tensor_2D%order_)) + call_julienne_assert(.all. (self%x_min_ .approximates. tensor_2D%x_min_ .within. 0D0)) + call_julienne_assert(.all. (self%x_max_ .approximates. tensor_2D%x_max_ .within. 0D0)) + conformable = .true. + end procedure + +end submodule diff --git a/src/formal/tensors_2D_m.F90 b/src/formal/tensors_2D_m.F90 index f5c5647..1c235be 100644 --- a/src/formal/tensors_2D_m.F90 +++ b/src/formal/tensors_2D_m.F90 @@ -56,6 +56,9 @@ pure function vector_2D_initializer_i(x,y) result(v) double precision x_max_(space_dimension) !! domain upper boundary integer cells_(space_dimension) !! number of grid cells spanning the domain integer order_ !! order of accuracy of mimetic discretization + contains + procedure, non_overridable, private :: tensor_2D_consistent + procedure, non_overridable, private :: tensor_2D_conformable end type interface tensor_2D_t @@ -81,10 +84,14 @@ pure module function construct_2D_tensor_from_components(values, cells, x_min, x generic :: values => scalar_2D_values generic :: grid => scalar_2D_grid generic :: to_file => scalar_2D_to_file + generic :: consistent => scalar_2D_consistent + generic :: conformable => scalar_2D_conformable_scalar procedure, non_overridable, private :: scalar_2D_to_file procedure, non_overridable, private :: scalar_2D_gradient procedure, non_overridable, private :: scalar_2D_values procedure, non_overridable, private :: scalar_2D_grid + procedure, non_overridable, private :: scalar_2D_consistent + procedure, non_overridable, private :: scalar_2D_conformable_scalar end type interface scalar_2D_t @@ -118,11 +125,16 @@ pure module function construct_2D_scalar_from_mold(initializer, mold) result(sca generic :: values => vector_2D_values generic :: to_file => vector_2D_to_file generic :: grid => vector_2D_grid + generic :: consistent => vector_2D_consistent + generic :: conformable => vector_2D_conformable_vector, vector_2D_conformable_scalar generic :: operator(.div.) => vector_2D_divergence procedure, non_overridable, private :: vector_2D_values 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 end type interface vector_2D_t @@ -169,6 +181,8 @@ pure module function construct_2D_vector_from_scalar_mold(initializer, mold) res generic :: values => divergence_2D_values generic :: grid => divergence_2D_grid generic :: to_file => divergence_2D_to_file + generic :: consistent => tensor_2D_consistent + generic :: conformable => tensor_2D_conformable procedure, private, non_overridable :: divergence_2D_values procedure, private, non_overridable :: divergence_2D_grid procedure, private, non_overridable :: divergence_2D_to_file @@ -199,14 +213,66 @@ pure module function construct_2D_divergence_from_vector_mold(initializer, mold) interface + pure module function tensor_2D_consistent(self) result(self_consistent) + !! Assert self-consistent tensor component shapes and values except for the value of the values_ component + implicit none + class(tensor_2D_t), intent(in) :: self + logical self_consistent + end function + + pure module function tensor_2D_conformable(self, tensor_2D) result(conformable) + !! Assert same-shaped self & vector_2D components + implicit none + class(tensor_2D_t), intent(in) :: self, tensor_2D + logical conformable + end function + + pure module function scalar_2D_consistent(self) result(self_consistent) + !! Assert components allocated and self-consistent, including sufficient accuracy for gradient operator + implicit none + class(scalar_2D_t), intent(in) :: self + logical self_consistent + end function + + pure module function scalar_2D_conformable_scalar(self, scalar_2D) result(conformable) + !! Assert the arguments' components are conformable, self-consistent, and consistent with each other + implicit none + class(scalar_2D_t), intent(in) :: self, scalar_2D + logical conformable + end function + + pure module function vector_2D_consistent(self) result(self_consistent) + !! Assert components allocated and self-consistent, including sufficient accuracy for divergence operator + implicit none + class(vector_2D_t), intent(in) :: self + logical self_consistent + end function + + pure module function vector_2D_conformable_vector(self, vector_2D) result(conformable) + !! Assert the arguments' components are conformable, self-consistent, and consistent with each other + implicit none + class(vector_2D_t), intent(in) :: self, vector_2D + logical conformable + end function + + pure module function vector_2D_conformable_scalar(self, scalar_2D) result(conformable) + !! Assert components allocated and self-consistent + implicit none + class(vector_2D_t), intent(in) :: self + class(scalar_2D_t), intent(in) :: scalar_2D + logical conformable + end function + pure module function scalar_2D_values(self) result(scalar_values) !! Scalar values getter + implicit none class(scalar_2D_t), intent(in) :: self double precision, allocatable :: scalar_values(:,:) end function pure module function scalar_2D_grid(self, direction) result(scalar_grid_1D) !! Result array contains scalar grid locations along the requested spatial direction + implicit none class(scalar_2D_t), intent(in) :: self integer, intent(in) :: direction double precision, allocatable :: scalar_grid_1D(:) @@ -214,6 +280,7 @@ pure module function scalar_2D_grid(self, direction) result(scalar_grid_1D) pure module function vector_2D_grid(self, direction) result(vector_grid_1D) !! Result array contains scalar grid locations along the requested spatial direction + implicit none class(vector_2D_t), intent(in) :: self integer, intent(in) :: direction double precision, allocatable :: vector_grid_1D(:) !! grid points along the requested coordinate direction @@ -221,6 +288,7 @@ pure module function vector_2D_grid(self, direction) result(vector_grid_1D) pure module function divergence_2D_grid(self, direction) result(divergence_grid_1D) !! Result array contains divergence grid locations along the requested spatial direction + implicit none class(divergence_2D_t), intent(in) :: self integer, intent(in) :: direction double precision, allocatable :: divergence_grid_1D(:) !! grid points along the requested coordinate direction @@ -228,12 +296,14 @@ pure module function divergence_2D_grid(self, direction) result(divergence_grid_ pure module function vector_2D_values(self) result(vector_values) !! Vector values getter + implicit none class(vector_2D_t), intent(in) :: self double precision, allocatable :: vector_values(:,:,:) end function pure module function divergence_2D_values(self) result(divergence_values) !! Vector values getter + implicit none class(divergence_2D_t), intent(in) :: self double precision, allocatable :: divergence_values(:,:) end function diff --git a/src/formal/vector_2D_s.F90 b/src/formal/vector_2D_s.F90 index 2c4fc02..02bf621 100644 --- a/src/formal/vector_2D_s.F90 +++ b/src/formal/vector_2D_s.F90 @@ -7,10 +7,12 @@ use julienne_m, only : & call_julienne_assert_ & ,operator(.all.) & + ,operator(.approximates.) & ,operator(.csv.) & ,operator(.equalsExpected.) & ,operator(.greaterThan.) & ,operator(.isAtLeast.) & + ,operator(.within.) & ,string_t use tensors_1D_m, only : faces_1D, vector_1D_t implicit none @@ -19,63 +21,117 @@ contains - module procedure construct_2D_vector_from_function + module procedure vector_2D_consistent + call_julienne_assert(self%tensor_2D_consistent()) + call_julienne_assert(.all. (self%cells_ .isAtLeast. 2*self%order_ + 1)) + call_julienne_assert(size(self%divergence_operator_1D_) .equalsExpected. space_dimension) + self_consistent = .true. + end procedure - integer dir + module procedure vector_2D_conformable_vector + call_julienne_assert(vector_2D_consistent(self)) + call_julienne_assert(vector_2D_consistent(vector_2D)) + call_julienne_assert(self%tensor_2D_conformable(vector_2D)) + conformable = .true. + end procedure - 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)) + module procedure vector_2D_conformable_scalar + call_julienne_assert(vector_2D_consistent(self)) + call_julienne_assert(scalar_2D_consistent(scalar_2D)) + call_julienne_assert(self%tensor_2D_conformable(scalar_2D)) + conformable = .true. + end procedure + + module procedure construct_2D_vector_from_function - associate(x => faces_1D(x_min(1), x_max(1), cells(1)), y => faces_1D(x_min(2), x_max(2), cells(2))) + define_grid: & + associate( & + x => faces_1D(x_min(1), x_max(1), cells(1)) & + ,y => faces_1D(x_min(2), x_max(2), cells(2)) & + ) + define_parent_tensor: & associate(vector_values => initializer(x,y)) vector_2D%tensor_2D_t = tensor_2D_t( & values = reshape(vector_values, shape=[shape(vector_values),1,1,1]) & ,cells = cells, x_min = x_min, x_max = x_max, order = order & ) - end associate - vector_2D%divergence_operator_1D_ = [(divergence_operator_1D_t(k=order, dx=((x_max(dir)-x_min(dir))/cells(dir)), cells=cells(dir)), dir=1,space_dimension)] - end associate + end associate define_parent_tensor + + define_divergence_operators: & + block + integer dir + vector_2D%divergence_operator_1D_ = & + [(divergence_operator_1D_t(k=order, dx=((x_max(dir)-x_min(dir))/cells(dir)), cells=cells(dir)), dir=1,space_dimension)] + end block define_divergence_operators + + end associate define_grid + + call_julienne_assert( vector_2D%consistent() ) + end procedure module procedure construct_2D_vector_from_vector_mold - integer dir - call_julienne_assert(.all. ([size(mold%cells_), size(mold%x_min_), size(mold%x_max_)] .equalsExpected. space_dimension)) - call_julienne_assert(.all. (mold%x_max_ .greaterThan. mold%x_min_)) - call_julienne_assert(.all. (mold%cells_ .isAtLeast. 2*mold%order_)) + call_julienne_assert( mold%consistent() ) - associate(x => faces_1D(mold%x_min_(1), mold%x_max_(1), mold%cells_(1)), y => faces_1D(mold%x_min_(2), mold%x_max_(2), mold%cells_(2))) + define_grid: & + associate( & + x => faces_1D(mold%x_min_(1), mold%x_max_(1), mold%cells_(1)) & + ,y => faces_1D(mold%x_min_(2), mold%x_max_(2), mold%cells_(2)) & + ) + define_parent_tensor: & associate(vector_values => initializer(x,y)) vector_2D%tensor_2D_t = tensor_2D_t( & values = reshape(vector_values, shape=[shape(vector_values),1,1,1]) & ,cells = mold%cells_, x_min = mold%x_min_, x_max = mold%x_max_, order = mold%order_ & ) - end associate - vector_2D%divergence_operator_1D_ = [(divergence_operator_1D_t(k=mold%order_, dx=((mold%x_max_(dir)-mold%x_min_(dir))/mold%cells_(dir)), cells=mold%cells_(dir)), dir=1,space_dimension)] - end associate + end associate define_parent_tensor + + define_divergence_operators: & + block + integer dir + vector_2D%divergence_operator_1D_ = [( & + divergence_operator_1D_t(k=mold%order_, dx=((mold%x_max_(dir)-mold%x_min_(dir))/mold%cells_(dir)), cells=mold%cells_(dir)) & + ,dir = 1, space_dimension & + )] + end block define_divergence_operators + + end associate define_grid + + call_julienne_assert( vector_2D%conformable(mold) ) + end procedure module procedure construct_2D_vector_from_scalar_mold integer dir - call_julienne_assert(.all. ([size(mold%cells_), size(mold%x_min_), size(mold%x_max_)] .equalsExpected. space_dimension)) - call_julienne_assert(.all. (mold%x_max_ .greaterThan. mold%x_min_)) - call_julienne_assert(.all. (mold%cells_ .isAtLeast. 2*mold%order_)) + call_julienne_assert( mold%consistent() ) - associate(x => faces_1D(mold%x_min_(1), mold%x_max_(1), mold%cells_(1)), y => faces_1D(mold%x_min_(2), mold%x_max_(2), mold%cells_(2))) + define_grid: & + associate( & + x => faces_1D(mold%x_min_(1), mold%x_max_(1), mold%cells_(1)) & + ,y => faces_1D(mold%x_min_(2), mold%x_max_(2), mold%cells_(2)) & + ) + define_parent_tensor: & associate(vector_values => initializer(x,y)) vector_2D%tensor_2D_t = tensor_2D_t( & values = reshape(vector_values, shape=[shape(vector_values),1,1,1]) & ,cells = mold%cells_, x_min = mold%x_min_, x_max = mold%x_max_, order = mold%order_ & ) - end associate - vector_2D%divergence_operator_1D_ = [(divergence_operator_1D_t(k=mold%order_, dx=((mold%x_max_(dir)-mold%x_min_(dir))/mold%cells_(dir)), cells=mold%cells_(dir)), dir=1,space_dimension)] - end associate + end associate define_parent_tensor + + vector_2D%divergence_operator_1D_ = [( & + divergence_operator_1D_t(k=mold%order_, dx=((mold%x_max_(dir)-mold%x_min_(dir))/mold%cells_(dir)), cells=mold%cells_(dir)) & + ,dir = 1, space_dimension & + )] + end associate define_grid + + call_julienne_assert( vector_2D%conformable(mold) ) + end procedure module procedure vector_2D_values - call_julienne_assert(allocated(self%values_)) + call_julienne_assert(self%consistent()) vector_values = self%values_(:,:,:,1,1,1) end procedure @@ -83,7 +139,7 @@ type(string_t), allocatable :: lines(:) integer i, j, l - call_julienne_assert(allocated(self%values_)) + call_julienne_assert(self%consistent()) associate(x => self%grid(x_dir), y => self%grid(y_dir), header => [string_t("x,y,vector_x,vector_y")]) associate(num_blank_lines => size(y)-1) @@ -120,7 +176,7 @@ module procedure vector_2D_divergence - call_julienne_assert(allocated(self%values_)) + call_julienne_assert(self%consistent()) divergence_2D%x_min_ = self%x_min_ divergence_2D%x_max_ = self%x_max_ @@ -143,6 +199,8 @@ end associate end do add_y_term + call_julienne_assert(divergence_2D%conformable(self)) + end procedure end submodule vector_2D_s \ No newline at end of file diff --git a/test/scalar_2D_test_m.F90 b/test/scalar_2D_test_m.F90 index d1a470f..11b729f 100644 --- a/test/scalar_2D_test_m.F90 +++ b/test/scalar_2D_test_m.F90 @@ -72,7 +72,10 @@ function check_gradient() result(test_diagnosis) do order = 2, 4, 2 associate(scalar_2D => scalar_2D_t(scalar_2D_initializer, order=order, cells=[30,20], x_min=[-1D0,1D0], x_max=[9D0,4D0])) - associate(grad_scalar => .grad. scalar_2D, expected_gradient => vector_2D_t(expected_gradient_initializer, mold=scalar_2D)) + associate( & + grad_scalar => .grad. scalar_2D & + ,expected_gradient => vector_2D_t(expected_gradient_initializer, mold=scalar_2D) & + ) test_diagnosis = test_diagnosis .also. & .all. (grad_scalar%values() .approximates. expected_gradient%values() .within. tolerance) & // string_t(" for order ") // string_t(order) From 5e7e3a73370b3abf2be82caa036b286ac6730a51 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Sat, 30 May 2026 20:24:13 -0500 Subject: [PATCH 36/36] feat(divergence_2D): operator(*) for const operand --- src/formal/divergence_2D_s.F90 | 16 ++++++++++++++++ src/formal/tensors_2D_m.F90 | 21 ++++++++++++++++++++- src/formal/vector_2D_s.F90 | 16 ++++++++++++---- 3 files changed, 48 insertions(+), 5 deletions(-) diff --git a/src/formal/divergence_2D_s.F90 b/src/formal/divergence_2D_s.F90 index dfc1e83..24f4bea 100644 --- a/src/formal/divergence_2D_s.F90 +++ b/src/formal/divergence_2D_s.F90 @@ -95,4 +95,20 @@ file = file_t(lines) end procedure + module procedure divergence_2D_postmultiply_constant + + call_julienne_assert(lhs%consistent()) + + product%tensor_2D_t = & + tensor_2D_t(values = lhs%values_ * rhs, cells = lhs%cells_, x_min = lhs%x_min_, x_max = lhs%x_max_, order = lhs%order_) + + call_julienne_assert(product%consistent()) + + end procedure + + module procedure divergence_2D_premultiply_constant + product = rhs * lhs + end procedure + + end submodule divergence_2D_s \ No newline at end of file diff --git a/src/formal/tensors_2D_m.F90 b/src/formal/tensors_2D_m.F90 index 1c235be..6927855 100644 --- a/src/formal/tensors_2D_m.F90 +++ b/src/formal/tensors_2D_m.F90 @@ -183,9 +183,12 @@ pure module function construct_2D_vector_from_scalar_mold(initializer, mold) res generic :: to_file => divergence_2D_to_file generic :: consistent => tensor_2D_consistent generic :: conformable => tensor_2D_conformable + generic :: operator(*) => divergence_2D_premultiply_constant, divergence_2D_postmultiply_constant procedure, private, non_overridable :: divergence_2D_values procedure, private, non_overridable :: divergence_2D_grid procedure, private, non_overridable :: divergence_2D_to_file + procedure, private, non_overridable, private :: divergence_2D_postmultiply_constant + procedure, private, non_overridable, pass(rhs) :: divergence_2D_premultiply_constant end type interface divergence_2D_t @@ -316,12 +319,28 @@ pure module function scalar_2D_gradient(self) result(gradient_2D) end function pure module function vector_2D_divergence(self) result(divergence_2D) - !! Result is mimetic divergence of the scalar_2D_t "self" + !! Result is mimetic divergence of the vector_2D_t "self" implicit none class(vector_2D_t), intent(in) :: self type(divergence_2D_t) divergence_2D end function + pure module function divergence_2D_postmultiply_constant(lhs, rhs) result(product) + !! Result is product of the divergence_2D_t lhs and the constant rhs + implicit none + class(divergence_2D_t), intent(in) :: lhs + double precision, intent(in) :: rhs + type(divergence_2D_t) product + end function + + pure module function divergence_2D_premultiply_constant(lhs, rhs) result(product) + !! Result is product of the divergence_2D_t rhs and the constant lhs + implicit none + class(divergence_2D_t), intent(in) :: rhs + double precision, intent(in) :: lhs + type(divergence_2D_t) product + end function + pure module function scalar_2D_to_file(self) result(file) !! Result is a file_t object containing the grid points and the corresponding scalar values implicit none diff --git a/src/formal/vector_2D_s.F90 b/src/formal/vector_2D_s.F90 index 02bf621..5a5e02d 100644 --- a/src/formal/vector_2D_s.F90 +++ b/src/formal/vector_2D_s.F90 @@ -120,10 +120,15 @@ ) end associate define_parent_tensor - vector_2D%divergence_operator_1D_ = [( & - divergence_operator_1D_t(k=mold%order_, dx=((mold%x_max_(dir)-mold%x_min_(dir))/mold%cells_(dir)), cells=mold%cells_(dir)) & - ,dir = 1, space_dimension & - )] + define_divergence_operators: & + block + integer dir + vector_2D%divergence_operator_1D_ = [( & + divergence_operator_1D_t(k=mold%order_, dx=((mold%x_max_(dir)-mold%x_min_(dir))/mold%cells_(dir)), cells=mold%cells_(dir)) & + ,dir = 1, space_dimension & + )] + end block define_divergence_operators + end associate define_grid call_julienne_assert( vector_2D%conformable(mold) ) @@ -131,8 +136,11 @@ end procedure module procedure vector_2D_values + call_julienne_assert(self%consistent()) + vector_values = self%values_(:,:,:,1,1,1) + end procedure module procedure vector_2D_to_file