diff --git a/.gitignore b/.gitignore index c73fb8d..619556d 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,5 @@ +scratch + # Prerequisites *.d 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/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/velocity-potential.gnuplot b/example/scripts/velocity-potential.gnuplot new file mode 100644 index 0000000..82668b0 --- /dev/null +++ b/example/scripts/velocity-potential.gnuplot @@ -0,0 +1,34 @@ +# ============================================================ +# velocity-potential.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 +# ============================================================ + +base_name = "velocity-potential" +datafile = base_name . ".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 base_name . ".gif" + +splot datafile every ::1 using 1:2:3 with pm3d title "" + +set output # flush and close the file 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 diff --git a/src/formal/scalar_1D_s.F90 b/src/formal/scalar_1D_s.F90 index cabd16f..a7f0b06 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) @@ -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_ & @@ -143,14 +159,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 @@ -196,19 +210,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/scalar_2D_s.F90 b/src/formal/scalar_2D_s.F90 new file mode 100644 index 0000000..e6d23c6 --- /dev/null +++ b/src/formal/scalar_2D_s.F90 @@ -0,0 +1,110 @@ +! 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, scalar_1D_t + use julienne_m, only : string_t, operator(.csv.) + implicit none + +contains + + module procedure scalar_2D_values + scalar_values = self%values_(:,:,1,1,1,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(x => cell_centers_extended_1D(x_min(1), x_max(1), cells(1)), y => cell_centers_extended_1D(x_min(2), x_max(2), cells(2))) + scalar_2D%tensor_2D_t = tensor_2D_t( & + values = reshape(initializer(x,y), shape=[size(x),size(y),1,1,1,1]) & + ,cells = cells , x_min = x_min, x_max = x_max, order = order & + ) + scalar_2D%gradient_operator_1D_ = gradient_operator_1D_t(k=order, dx=(x_max - x_min)/cells, cells=cells) + end associate + end procedure + + 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 scalar_2D_gradient + + integer c, i, j + + gradient_2D%x_min_ = self%x_min_ + gradient_2D%x_max_ = self%x_max_ + gradient_2D%cells_ = self%cells_ + gradient_2D%order_ = self%order_ + + 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: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: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 + + 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 + + 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/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_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/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_1D_m.F90 b/src/formal/tensors_1D_m.F90 index 42b2823..59ceae9 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 @@ -23,6 +22,8 @@ module tensors_1D_m public :: vector_1D_initializer_i public :: d_dx public :: d2_dx2 + public :: cell_centers_extended_1D + public :: faces_1D abstract interface @@ -119,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 @@ -168,10 +180,20 @@ 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 + + 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 @@ -496,10 +518,29 @@ 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 + 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) + 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) double precision, intent(in) :: x_min, x_max integer, intent(in) :: cells diff --git a/src/formal/tensors_2D_m.F90 b/src/formal/tensors_2D_m.F90 new file mode 100644 index 0000000..3f1a83e --- /dev/null +++ b/src/formal/tensors_2D_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_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 + use julienne_m, only : file_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(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 + implicit none + double precision, intent(in) :: x(:), y(:) + double precision v(size(x),size(y),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, 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 + integer order_ !! order of accuracy of mimetic discretization + 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 + type(gradient_operator_1D_t) gradient_operator_1D_(space_dimension) + contains + generic :: operator(.grad.) => scalar_2D_gradient + generic :: values => scalar_2D_values + generic :: grid => scalar_2D_grid + 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 + + 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 + + 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 + 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 + + 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 + + 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 + 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 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_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 + double precision, allocatable :: vector_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 + 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 + + 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/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_1D_s.F90 b/src/formal/vector_1D_s.F90 index 8dc9e34..1467d38 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) @@ -60,34 +60,44 @@ 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 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 -#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 @@ -96,19 +106,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 diff --git a/src/formal/vector_2D_s.F90 b/src/formal/vector_2D_s.F90 new file mode 100644 index 0000000..73f9554 --- /dev/null +++ b/src/formal/vector_2D_s.F90 @@ -0,0 +1,117 @@ +! 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(.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_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(x => faces_1D(x_min(1), x_max(1), cells(1)), y => faces_1D(x_min(2), x_max(2), cells(2))) + 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 + + 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_)) + + 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))) + 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 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))) + 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_(:,:,:,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_)/space_dimension + 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 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 2fc9526..86d4319 100644 --- a/src/formal_m.f90 +++ b/src/formal_m.f90 @@ -18,6 +18,20 @@ 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 + ,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 65809b8..bffe1e5 100644 --- a/test/driver.f90 +++ b/test/driver.f90 @@ -10,6 +10,8 @@ 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 + use scalar_3D_test_m, only : scalar_3D_test_t implicit none associate(test_harness => test_harness_t([ & @@ -19,6 +21,8 @@ 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()) & + ,test_fixture_t(scalar_3D_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..d1a470f --- /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 = 5D-2 + 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 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)) 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 + + 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 => 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=[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) & + // 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 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