From dd5602c99dda20fd5aad32536e0fb69c3a5320ae Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Wed, 13 May 2026 18:56:45 -0700 Subject: [PATCH 01/34] 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 f54ee72038d7d2b38e517c19cb2c161c421c06bf Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Wed, 13 May 2026 19:00:29 -0700 Subject: [PATCH 02/34] 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 e11993f49a89054729dba9ecc469cede10af3df5 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Wed, 13 May 2026 19:03:51 -0700 Subject: [PATCH 03/34] 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 2c88e85994fe75fb497c333a3b6aa790db3c5048 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Wed, 13 May 2026 19:09:21 -0700 Subject: [PATCH 04/34] 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 35a67166a4d8341f9cec0906fdd6f864f8e7484f Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Wed, 13 May 2026 19:58:39 -0700 Subject: [PATCH 05/34] 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 eb92f783d02b989261901a3c8a367c8e4421ff55 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Wed, 13 May 2026 20:23:26 -0700 Subject: [PATCH 06/34] 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 4ea3c61470a47752c81da38a91bffb6a33fcaf7a Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Wed, 13 May 2026 20:33:19 -0700 Subject: [PATCH 07/34] 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 02124fbaefe48d00f76c276d0e2fb49738274725 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Sat, 16 May 2026 17:53:24 -0700 Subject: [PATCH 08/34] 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 cc37fb54988eda95475d0d0d678a5e625acbfd75 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Sun, 17 May 2026 01:05:21 -0700 Subject: [PATCH 09/34] 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 79c14a7b4399f06d3422201ee711ca0544825e7c Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Sun, 17 May 2026 01:10:03 -0700 Subject: [PATCH 10/34] 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 de712560ab3e544d7a2dce70642cb453a5310380 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Sun, 17 May 2026 09:53:42 -0700 Subject: [PATCH 11/34] 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 9a2e1c123c43bba62dc8f5fe1d3c58f0ebcf33f2 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Sun, 17 May 2026 12:19:05 -0700 Subject: [PATCH 12/34] 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 0c8743ccf0b2d41f771e5fc378ff1beafc078dde Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Sun, 17 May 2026 17:25:02 -0700 Subject: [PATCH 13/34] 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 0f090ec45542a9add05c9060e62eb73cb98da4dc Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Sun, 17 May 2026 19:52:54 -0700 Subject: [PATCH 14/34] 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 3cdc7c53f27fdc2dd8ee723c2e763603fe1fb24e Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Sun, 17 May 2026 20:07:35 -0700 Subject: [PATCH 15/34] 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 f2da3dceb69242f3b886bd37b6b3e6c9156ddbef Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Sun, 17 May 2026 23:13:48 -0700 Subject: [PATCH 16/34] 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 7885382cf219b45a8bc063ef81d6d28e1aede8a4 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Sun, 17 May 2026 23:42:05 -0700 Subject: [PATCH 17/34] 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 e0560af888591449922b7e909ba8e7d8db72b094 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Mon, 18 May 2026 22:18:35 -0700 Subject: [PATCH 18/34] 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 32d3a2a8dd075c0631985e9c1084da693905c406 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Tue, 19 May 2026 18:37:30 -0700 Subject: [PATCH 19/34] 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 f3cc871d762f564b6f6ebb2685c3d34c75a56406 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Tue, 19 May 2026 18:57:53 -0700 Subject: [PATCH 20/34] 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 7b5692e43853c438214e03f134bff2d7d52c35ba Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Fri, 22 May 2026 10:42:34 -0700 Subject: [PATCH 21/34] 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 78c923c6ce56360a638f2f4163e028ae9e99f6b8 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Fri, 22 May 2026 10:53:01 -0700 Subject: [PATCH 22/34] 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 74bcb1fd5be1cac04d1020893e527bef26ee811b Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Fri, 22 May 2026 15:02:30 -0700 Subject: [PATCH 23/34] 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 6ce911317ab4960635e252fcf933f818f783d2e2 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Fri, 22 May 2026 15:56:15 -0700 Subject: [PATCH 24/34] 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 33b7252f4e7ae47a217046d7aa479e5afc052ae9 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Sun, 24 May 2026 22:20:45 -0700 Subject: [PATCH 25/34] 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 6cb3cee44160168b527a09fbe099daa19d894937 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Sun, 24 May 2026 23:01:17 -0700 Subject: [PATCH 26/34] 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 cf4cbadc252b43b826c318fa4e134143deddcad0 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Mon, 25 May 2026 22:56:18 -0700 Subject: [PATCH 27/34] 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 8adb61f3921f9fd6a035b6383b703c290354ea40 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Mon, 25 May 2026 23:47:08 -0700 Subject: [PATCH 28/34] 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 5675810e9fdebad94db5a5dedd2a2b205f2b937d Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Tue, 26 May 2026 10:42:36 -0700 Subject: [PATCH 29/34] 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 d627c2bb0f74d3020a5eac23a6ed8f664e7bdb48 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Tue, 26 May 2026 10:49:38 -0700 Subject: [PATCH 30/34] 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 c815b5b8ed4f68ede8eb07746203f25c8730a097 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Tue, 26 May 2026 11:24:55 -0700 Subject: [PATCH 31/34] 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 2e2e540f4b80dd5296e9f71f6c930dcf2057c216 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Tue, 26 May 2026 11:27:47 -0700 Subject: [PATCH 32/34] 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 b5ef8c118811a233f38f9c698e389cb10261933a Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Fri, 29 May 2026 14:56:36 -0500 Subject: [PATCH 33/34] doc(README): update lfortran version to "latest" --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 938adbf..7d5341d 100644 --- a/README.md +++ b/README.md @@ -76,7 +76,7 @@ Building and testing GCC | `gfortran` | 14-15 | `fpm test --compiler gfortran --profile release` GCC | `gfortran` | 13 | `fpm test --compiler gfortran --profile release --flag "-ffree-line-length-none"` Intel | `ifx` | 2025.1.2 | `FOR_COARRAY_NUM_IMAGES=1 fpm test --compiler ifx --flag "-fpp -O3 -coarray" --profile release` - LFortran | `lfortran` | 0.60.0-421-ge2c448c79 | `fpm test --compiler lfortran --flag "--cpp --realloc-lhs-arrays"` + LFortran | `lfortran` | latest | `fpm test --compiler lfortran --flag "--cpp --realloc-lhs-arrays"` LLVM | `flang` | 20-21 | `fpm test --compiler flang --profile release LLVM | `flang` | 19 | `fpm test --compiler flang --profile release --flag "-mmlir -allow-assumed-rank"` NAG | `nagfor` | 7.2 Build 7242 | `fpm test --compiler nagfor --flag "-O3 -fpp"` From 3b5b98a055341d206cfbe70244dd4a7fc51bee90 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Fri, 29 May 2026 15:00:14 -0500 Subject: [PATCH 34/34] 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"