diff --git a/README.md b/README.md index 977e80b..938adbf 100644 --- a/README.md +++ b/README.md @@ -84,11 +84,11 @@ Building and testing ### `fpm` versions before 0.13.0 With LLVM 20-22, replace the above `flang` command with ``` - fpm test --compiler flang-new --flag "-O3" + fpm test --compiler flang --flag "-O3" ``` With LLVM 19, replace the above `flang` command with ``` - fpm test --compiler flang-new --flag "-O3 -mmlir -allow-assumed-rank" + fpm test --compiler flang --flag "-O3 -mmlir -allow-assumed-rank" ``` Documentation diff --git a/example/burgers-1D.F90 b/example/burgers-1D.F90 new file mode 100644 index 0000000..6002e1a --- /dev/null +++ b/example/burgers-1D.F90 @@ -0,0 +1,97 @@ +! Copyright (c) 2026, The Regents of the University of California +! Terms of use are as specified in LICENSE.txt + +module initial_condition_m + implicit none + +contains + + pure function initial_condition(x) + !! Initial solution to Burgers equation + double precision, intent(in) :: x(:) + double precision, allocatable :: initial_condition(:) + initial_condition = 10*sin(x) + ! To change this function, please edit only the right-hand-side (RHS) expression, + ! keeping the rest in place for proper display of the function at runtime. + end function + +end module + +program burgers_1D + !! This program demonstrates the use of Formal to solve the partial differential equation of + !! Burgers (1948) in conservative form using the 2nd- or 4th-order mimetic discretizations of + !! Corbino & Castillo (2020) and Dumett & Castillo (2022). + !! + !! * Burgers, J.M. (1948) https://doi.org/10.1016/S0065-2156(08)70100-5 + !! * Corbino & Castillo (2020) https://doi.org/10.1016/j.cam.2019.06.042. + !! * Dumett & Castillo (2022) https://doi.org/10.13140/RG.2.2.26630.14400 + use initial_condition_m, only : initial_condition + use julienne_m, only : command_line_t + use formal_m, only : scalar_1D_t, scalar_1D_initializer_i, d_dx, d2_dx2 + 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 + + if (command_line%argument_present([character(len=len("--help")) :: ("--help"), "-h"])) then + stop new_line('') // new_line('') & + // 'Usage:' // new_line('') & + // ' fpm run \' // new_line('') & + // ' --example burgers-1D \' // new_line('') & + // ' --compiler flang \' // new_line('') & + // ' --flag "-O3" \' // new_line('') & + // ' [--help|-h] | [--order ]' // new_line('') // new_line('') & + // 'where square brackets indicate optional arguments and angular brackets indicate user input values.' // new_line('') + end if + + print *, new_line('') + print *," Initial condition" + print *," =================" + + call execute_command_line("grep 'initial_condition =' example/burgers-1D.F90 | grep -v execute_command", wait=.true.) + + order_string = command_line%flag_value("--order") + + if (len(order_string)==0) then + order = 2 + else + read(order_string,"(i1)") order + end if + + print *, "order = ", order + + block +#ifndef __GFORTRAN__ + procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer +#endif + double precision, parameter :: dt = 1D0, pi = acos(-1D0) + type(scalar_1D_t) u + + scalar_1D_initializer => initial_condition + u = scalar_1D_t(scalar_1D_initializer, order, x_min=0D0, x_max=2*pi, cells=20) + + runge_kutta_2nd_order: & + associate(u_half => u + d_dt(u)*dt/2) ! first substep + u = u + d_dt(u_half)*dt ! second substep + end associate runge_kutta_2nd_order + + end block + +#ifdef __GFORTRAN__ + stop +#endif + +contains + + pure function d_dt(u) result(du_dt) + type(scalar_1D_t), intent(in) :: u + type(scalar_1D_t) du_dt + double precision, parameter :: nu=1D0 + du_dt = nu*d2_dx2(u) - d_dx((u**2)/2) + end function + +end program diff --git a/src/formal/mimetic_matrix_1D_s.F90 b/src/formal/differential_operator_matrix_1D_s.F90 similarity index 76% rename from src/formal/mimetic_matrix_1D_s.F90 rename to src/formal/differential_operator_matrix_1D_s.F90 index fd0b114..b40d573 100644 --- a/src/formal/mimetic_matrix_1D_s.F90 +++ b/src/formal/differential_operator_matrix_1D_s.F90 @@ -1,16 +1,16 @@ ! Copyright (c) 2026, The Regents of the University of California ! Terms of use are as specified in LICENSE.txt -submodule(mimetic_operators_1D_m) mimetic_matrix_1D_s +submodule(differential_operators_1D_m) differential_operator_matrix_1D_s use julienne_m, only : string_t, operator(.csv.) implicit none contains module procedure construct_matrix_operator - mimetic_matrix_1D%upper_ = upper - mimetic_matrix_1D%inner_ = inner - mimetic_matrix_1D%lower_ = lower + differential_operator_matrix_1D%upper_ = upper + differential_operator_matrix_1D%inner_ = inner + differential_operator_matrix_1D%lower_ = lower end procedure module procedure to_file_t @@ -33,4 +33,4 @@ end procedure -end submodule mimetic_matrix_1D_s \ No newline at end of file +end submodule differential_operator_matrix_1D_s \ No newline at end of file diff --git a/src/formal/mimetic_operators_1D_m.F90 b/src/formal/differential_operators_1D_m.F90 similarity index 91% rename from src/formal/mimetic_operators_1D_m.F90 rename to src/formal/differential_operators_1D_m.F90 index 22a4670..41f7bc2 100644 --- a/src/formal/mimetic_operators_1D_m.F90 +++ b/src/formal/differential_operators_1D_m.F90 @@ -3,7 +3,7 @@ #include "formal-language-support.F90" -module mimetic_operators_1D_m +module differential_operators_1D_m !! Define sparse matrix storage formats and operators tailored to the one-dimensional (1D) mimetic discretizations !! detaild by Corbino & Castillo (2020) https://doi.org/10.1016/j.cam.2019.06.042. use julienne_m, only : file_t @@ -13,8 +13,8 @@ module mimetic_operators_1D_m public :: gradient_operator_1D_t public :: divergence_operator_1D_t - type mimetic_matrix_1D_t - !! Encapsulate a mimetic matrix with a corresponding matrix-vector product operator + type differential_operator_matrix_1D_t + !! Encapsulate a mimetic matrix private double precision, allocatable :: upper_(:,:) !! A submatrix block (cf. Corbino & Castillo, 2020) double precision, allocatable :: inner_(:) !! M submatrix row (cf. Corbino & Castillo, 2020) @@ -23,20 +23,20 @@ module mimetic_operators_1D_m procedure, non_overridable :: to_file_t end type - interface mimetic_matrix_1D_t + interface differential_operator_matrix_1D_t - pure module function construct_matrix_operator(upper, inner, lower) result(mimetic_matrix_1D) + pure module function construct_matrix_operator(upper, inner, lower) result(differential_operator_matrix_1D) !! Construct discrete operator from matrix blocks implicit none double precision, intent(in) :: upper(:,:) !! A submatrix block (cf. Corbino & Castillo, 2020) double precision, intent(in) :: inner(:) !! M submatrix row (cf. Corbino & Castillo, 2020) double precision, intent(in) :: lower(:,:) !! A' submatrix block (cf. Corbino & Castillo, 2020) - type(mimetic_matrix_1D_t) mimetic_matrix_1D + type(differential_operator_matrix_1D_t) differential_operator_matrix_1D end function end interface - type, extends(mimetic_matrix_1D_t) :: gradient_operator_1D_t + type, extends(differential_operator_matrix_1D_t) :: gradient_operator_1D_t !! Encapsulate a 1D mimetic gradient operator private integer k_ !! order of accuracy @@ -62,7 +62,7 @@ pure module function construct_1D_gradient_operator(k, dx, cells) result(gradien end interface - type, extends(mimetic_matrix_1D_t) :: divergence_operator_1D_t + type, extends(differential_operator_matrix_1D_t) :: divergence_operator_1D_t !! Encapsulate kth-order mimetic divergence operator on m_ cells of width dx private integer k_, m_ @@ -129,7 +129,7 @@ pure module function divergence_matrix_multiply(self, vec) result(matvec_product pure module function to_file_t(self) result(file) implicit none - class(mimetic_matrix_1D_t), intent(in) :: self + class(differential_operator_matrix_1D_t), intent(in) :: self type(file_t) file end function @@ -164,4 +164,4 @@ pure function negate_and_flip(A) result(Ap) #endif -end module mimetic_operators_1D_m +end module differential_operators_1D_m diff --git a/src/formal/divergence_operator_1D_s.F90 b/src/formal/divergence_operator_1D_s.F90 index 3da38d0..68c23d4 100644 --- a/src/formal/divergence_operator_1D_s.F90 +++ b/src/formal/divergence_operator_1D_s.F90 @@ -4,7 +4,7 @@ #include "formal-language-support.F90" #include "julienne-assert-macros.h" -submodule(mimetic_operators_1D_m) divergence_operator_1D_s +submodule(differential_operators_1D_m) divergence_operator_1D_s use julienne_m, only : call_julienne_assert_, string_t #if ASSERTIONS use julienne_m, only : operator(.isAtLeast.), operator(.equalsExpected.) @@ -48,7 +48,7 @@ pure function negate_and_flip(A) result(Ap) else allocate(Ap, mold = A) end if - divergence_operator_1D%mimetic_matrix_1D_t = mimetic_matrix_1D_t(A, M(k, dx), Ap) + divergence_operator_1D%differential_operator_matrix_1D_t = differential_operator_matrix_1D_t(A, M(k, dx), Ap) divergence_operator_1D%k_ = k divergence_operator_1D%dx_ = dx divergence_operator_1D%m_ = cells diff --git a/src/formal/gradient_operator_1D_s.F90 b/src/formal/gradient_operator_1D_s.F90 index 9f6ebed..e758b01 100644 --- a/src/formal/gradient_operator_1D_s.F90 +++ b/src/formal/gradient_operator_1D_s.F90 @@ -4,7 +4,7 @@ #include "julienne-assert-macros.h" #include "formal-language-support.F90" -submodule(mimetic_operators_1D_m) gradient_operator_1D_s +submodule(differential_operators_1D_m) gradient_operator_1D_s use julienne_m, only : call_julienne_assert_, string_t #if ASSERTIONS use julienne_m, only : operator(.isAtLeast.) @@ -42,7 +42,7 @@ pure function negate_and_flip(A) result(Ap) call_julienne_assert(cells .isAtLeast. 2*k) associate(A => corbino_castillo_A(k, dx), M => corbino_castillo_M(k, dx)) - gradient_operator_1D%mimetic_matrix_1D_t = mimetic_matrix_1D_t(A, M, negate_and_flip(A)) + gradient_operator_1D%differential_operator_matrix_1D_t = differential_operator_matrix_1D_t(A, M, negate_and_flip(A)) gradient_operator_1D%k_ = k gradient_operator_1D%dx_ = dx gradient_operator_1D%m_ = cells diff --git a/src/formal/interpolator_1D_m.F90 b/src/formal/interpolator_1D_m.F90 new file mode 100644 index 0000000..537c6c4 --- /dev/null +++ b/src/formal/interpolator_1D_m.F90 @@ -0,0 +1,78 @@ +! Copyright (c) 2026, The Regents of the University of California +! Terms of use are as specified in LICENSE.txt + +module interpolator_1D_m + !! Define a sparse matrix storage format tailored to the staggered-grid interpolation matrix + !! operators detailed by Dumett & Castillo (2022) https://doi.org/10.13140/RG.2.2.26630.14400 + implicit none + + private + public :: centers_to_faces_1D_t + public :: faces_to_centers_1D_t + + type, abstract :: interpolator_1D_t + !! Encapsulate a staggered-grid interpolation matrix with a corresponding matrix-vector product operator + private + integer order_, cells_, dx_ + double precision first_ + double precision, allocatable :: upper_(:,:) + double precision, allocatable :: inner_(:) + double precision, allocatable :: lower_(:,:) + double precision final_ + end type + + type, extends(interpolator_1D_t) :: centers_to_faces_1D_t + contains + procedure, non_overridable :: face_values + end type + + type, extends(interpolator_1D_t) :: faces_to_centers_1D_t + contains + procedure, non_overridable :: center_values + end type + + interface centers_to_faces_1D_t + + pure module function c2f_constructor(order, cells, dx) result(centers_to_faces_1D) + !! Construct centers-to-faces interpolation operator + implicit none + integer, intent(in) :: order, cells + double precision, intent(in) :: dx + type(centers_to_faces_1D_t) centers_to_faces_1D + end function + + end interface + + interface faces_to_centers_1D_t + + pure module function f2c_constructor(order, cells, dx) result(faces_to_centers_1D) + !! Construct centers-to-faces interpolation operator + implicit none + integer, intent(in) :: order, cells + double precision, intent(in) :: dx + type(faces_to_centers_1D_t) faces_to_centers_1D + end function + + end interface + + interface + + pure module function face_values(self, centers_extended) result(faces) + !! Interpolate cell-centered values to face-centered values + implicit none + class(centers_to_faces_1D_t), intent(in) :: self + double precision, intent(in) :: centers_extended(:) + double precision, allocatable :: faces(:) + end function + + pure module function center_values(self, faces) result(centers_extended) + !! Interpolate face-centered values to cell-centered values + implicit none + class(faces_to_centers_1D_t), intent(in) :: self + double precision, intent(in) :: faces(:) + double precision, allocatable :: centers_extended(:) + end function + + end interface + +end module interpolator_1D_m \ No newline at end of file diff --git a/src/formal/interpolator_1D_s.F90 b/src/formal/interpolator_1D_s.F90 new file mode 100644 index 0000000..d7e5592 --- /dev/null +++ b/src/formal/interpolator_1D_s.F90 @@ -0,0 +1,108 @@ +! 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(interpolator_1D_m) interpolator_1D_s + use julienne_m, only : call_julienne_assert_, operator(.all.), operator(.equalsExpected.) + implicit none + +contains + + module procedure c2f_constructor + + centers_to_faces_1D%order_ = order + centers_to_faces_1D%cells_ = cells + centers_to_faces_1D%dx_ = dx + + select case(order) + case(2) + centers_to_faces_1D%first_ = (2D0 )/2 + centers_to_faces_1D%upper_ = (reshape([double precision::], [0,3]))/2 + centers_to_faces_1D%inner_ = ([1D0,1D0] )/2 + centers_to_faces_1D%lower_ = (reshape([double precision::], [0,3]))/2 + centers_to_faces_1D%final_ = (2D0 )/2 + case(4) + centers_to_faces_1D%first_ = 1D0 + centers_to_faces_1D%upper_ = reshape([-16, 70, 70, -14, 2], [1,5])/112D0 + centers_to_faces_1D%inner_ = [ -7, 63, 63, -7] /112D0 + centers_to_faces_1D%lower_ = reshape([ 2, -14, 70, 70, -16], [1,5])/112D0 + centers_to_faces_1D%final_ = 1D0 + case default + error stop "c2f_component_constructor: unsupported order" + end select + + call_julienne_assert(.all. (shape(centers_to_faces_1D%lower_) .equalsExpected. shape(centers_to_faces_1D%upper_))) + end procedure + + module procedure f2c_constructor + + faces_to_centers_1D%order_ = order + faces_to_centers_1D%cells_ = cells + faces_to_centers_1D%dx_ = dx + + select case(order) + case(2) + faces_to_centers_1D%first_ = (2D0 )/2 + faces_to_centers_1D%upper_ = reshape([double precision::], [0,3])/2 + faces_to_centers_1D%inner_ = [1D0,1D0]/2 + faces_to_centers_1D%lower_ = reshape([double precision::], [0,3])/2 + faces_to_centers_1D%final_ = (2D0 )/2 + case(4) + faces_to_centers_1D%first_ = 1D0 + faces_to_centers_1D%upper_ = reshape([35, 140, -70, 28, -5], [1,5])/128D0 + faces_to_centers_1D%inner_ = [-8, 72, 72, -8] /128D0 + faces_to_centers_1D%lower_ = reshape([-5, 28, -70, 140, 35], [1,5])/128D0 + faces_to_centers_1D%final_ = 1D0 + case default + error stop "f2c_component_constructor: unsupported order" + end select + + call_julienne_assert(.all. (shape(faces_to_centers_1D%lower_) .equalsExpected. shape(faces_to_centers_1D%upper_))) + end procedure + + module procedure face_values + integer row + integer, parameter :: end_point = 1 + associate( & + N => size(centers_extended) , inner_cols => size(self%inner_) & + ,upper_rows => size(self%upper_,1), upper_cols => size(self%upper_,2) & + ,lower_rows => size(self%lower_,1), lower_cols => size(self%lower_,2) & + ) + call_julienne_assert(N .equalsExpected. self%cells_ + 2*end_point) + associate(inner_rows => N - (2*end_point + upper_rows + lower_rows)) + faces = [ self%first_ * centers_extended(1) & + , matmul(self%upper_ , centers_extended(1:upper_cols)) & + ,[(dot_product(self%inner_ , centers_extended(row - upper_rows : row - upper_rows + inner_cols - 1)), & + row = end_point + upper_rows + 1, end_point + upper_rows + inner_rows - 1)] & + , matmul(self%lower_ , centers_extended(N-lower_cols+1:N)) & + , self%final_ * centers_extended(N) & + ] + end associate + call_julienne_assert(size(faces) .equalsExpected. self%cells_ + 1) + end associate + end procedure + + module procedure center_values + integer row + integer, parameter :: end_point = 1 + associate( & + N => size(faces) , inner_cols => size(self%inner_) & + ,upper_rows => size(self%upper_,1), upper_cols => size(self%upper_,2) & + ,lower_rows => size(self%lower_,1), lower_cols => size(self%lower_,2) & + ) + call_julienne_assert(N .equalsExpected. self%cells_ + 1) + associate(inner_rows => N + 1 - (2*end_point + upper_rows + lower_rows)) + centers_extended = & + [ self%first_ * faces(1) & + , matmul(self%upper_ , faces(1:upper_cols)) & + ,[(dot_product(self%inner_ , faces(row : row + inner_cols - 1)), row = 1, inner_rows )] & + , matmul(self%lower_ , faces(N-lower_cols+1:N)) & + , self%final_ * faces(N) & + ] + end associate + call_julienne_assert(size(centers_extended) .equalsExpected. self%cells_ + 2) + end associate + end procedure + +end submodule interpolator_1D_s \ No newline at end of file diff --git a/src/formal/scalar_1D_s.F90 b/src/formal/scalar_1D_s.F90 index 47fa585..645f0dd 100644 --- a/src/formal/scalar_1D_s.F90 +++ b/src/formal/scalar_1D_s.F90 @@ -16,6 +16,7 @@ ,operator(.greaterThan.) & ,operator(.within.) & ,string_t + use interpolator_1D_m, only : faces_to_centers_1d_t implicit none contains @@ -65,6 +66,66 @@ pure function cell_center_locations(x_min, x_max, cells) result(x) #endif + 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_ & + ) + ratio%gradient_operator_1D_ = gradient_operator_1D_t( & + k = self%order_, dx = (self%x_max_ - self%x_min_)/self%cells_, cells = self%cells_ & + ) + end procedure + + 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%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 = & + 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(lhs%cells_ .equalsExpected. rhs%cells_) + 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 = & + 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 premultiply_double + lhs_x_rhs%gradient_operator_1D_ = rhs%gradient_operator_1D_ + lhs_x_rhs%tensor_1D_t = & + tensor_1D_t(values = lhs*rhs%values_, x_min = rhs%x_min_, x_max = rhs%x_max_, cells = rhs%cells_, order = rhs%order_) + end procedure + + module procedure postmultiply_double + lhs_x_rhs%gradient_operator_1D_ = lhs%gradient_operator_1D_ + lhs_x_rhs%tensor_1D_t = & + tensor_1D_t(values = rhs*lhs%values_, x_min = lhs%x_min_, x_max = lhs%x_max_, cells = lhs%cells_, order = lhs%order_) + end procedure + + module procedure exponentiate + power%tensor_1D_t = tensor_1D_t( & + values = self%values_**exponent, x_min = self%x_min_, x_max = self%x_max_, cells = self%cells_, order = self%order_ & + ) + power%gradient_operator_1D_ = gradient_operator_1D_t( & + k = self%order_, dx = (self%x_max_ - self%x_min_)/self%cells_, cells = self%cells_ & + ) + end procedure + + module procedure construct_1D_scalar_from_parent + + call_julienne_assert(tensor_1D%is_cell_centers_extended()) + + scalar_1D%tensor_1D_t = tensor_1D_t( & + values = tensor_1D%values_, x_min = tensor_1D%x_min_, x_max = tensor_1D%x_max_, cells = tensor_1D%cells_, order = tensor_1D%order_ & + ) + scalar_1D%gradient_operator_1D_ = gradient_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 grad integer c @@ -82,6 +143,29 @@ pure function cell_center_locations(x_min, x_max, cells) result(x) end procedure + module procedure d_dx + + associate( & + dx => (self%x_max_ - self%x_min_)/self%cells_ & + ,interpolator => faces_to_centers_1D_t(order=self%order_, cells=self%cells_, dx=(self%x_max_ - self%x_min_)/self%cells_) & + ) + dself_dx%gradient_operator_1D_ = gradient_operator_1D_t(self%order_, dx, self%cells_) + associate(tensor_1D => & + tensor_1D_t(dself_dx%gradient_operator_1D_ .x. self%values_, self%x_min_, self%x_max_, cells=self%cells_, order=self%order_) & + ) + dself_dx%tensor_1D_t = & + tensor_1D_t(interpolator%center_values(tensor_1D%values_), self%x_min_, self%x_max_, cells=self%cells_, order=self%order_) + end associate + end associate + + call_julienne_assert(dself_dx%is_cell_centers_extended()) + + end procedure + + module procedure d2_dx2 + d2_self_dx2 = d_dx(d_dx(self)) + end procedure + module procedure laplacian #ifndef __GFORTRAN__ diff --git a/src/formal/tensor_1D_s.f90 b/src/formal/tensor_1D_s.f90 index c67558b..c085da3 100644 --- a/src/formal/tensor_1D_s.f90 +++ b/src/formal/tensor_1D_s.f90 @@ -1,10 +1,25 @@ ! 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_1D_m) tensor_1D_s + use julienne_m, only : call_julienne_assert_ implicit none contains + module procedure is_cell_centered + is_cell_centered = size(self%values_) == self%cells_ + end procedure + + module procedure is_face_centered + is_face_centered = size(self%values_) == self%cells_ + 1 + end procedure + + module procedure is_cell_centers_extended + is_cell_centers_extended = size(self%values_) == self%cells_ + 2 + end procedure + module procedure construct_1D_tensor_from_components tensor_1D%values_ = values tensor_1D%x_min_ = x_min diff --git a/src/formal/tensors_1D_m.F90 b/src/formal/tensors_1D_m.F90 index 2615a7c..eecd06f 100644 --- a/src/formal/tensors_1D_m.F90 +++ b/src/formal/tensors_1D_m.F90 @@ -8,7 +8,7 @@ module tensors_1D_m !! 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 mimetic_operators_1D_m, only : divergence_operator_1D_t, gradient_operator_1D_t + use differential_operators_1D_m, only : divergence_operator_1D_t, gradient_operator_1D_t implicit none @@ -21,6 +21,8 @@ module tensors_1D_m public :: divergence_1D_t public :: scalar_1D_initializer_i public :: vector_1D_initializer_i + public :: d_dx + public :: d2_dx2 abstract interface @@ -51,6 +53,9 @@ pure function vector_1D_initializer_i(x) result(v) integer order_ !! order of accuracy of mimetic discretization double precision, allocatable :: values_(:) !! tensor components at spatial locations contains + procedure, non_overridable, private :: is_cell_centered + procedure, non_overridable, private :: is_face_centered + procedure, non_overridable, private :: is_cell_centers_extended procedure, non_overridable, private :: gradient_1D_weights procedure, non_overridable, private :: divergence_1D_weights generic :: dV => dx @@ -77,14 +82,25 @@ pure module function construct_1D_tensor_from_components(values, x_min, x_max, c private type(gradient_operator_1D_t) gradient_operator_1D_ contains - generic :: operator(.grad.) => grad - generic :: operator(.laplacian.) => laplacian generic :: grid => scalar_1D_grid generic :: values => scalar_1D_values - procedure, non_overridable, private :: grad - procedure, non_overridable, private :: laplacian + generic :: operator(-) => subtract_scalar_1D + generic :: operator(+) => add_scalar_1D + generic :: operator(/) => divide_by_integer + generic :: operator(*) => premultiply_double, postmultiply_double + generic :: operator(**) => exponentiate + generic :: operator(.grad.) => grad + generic :: operator(.laplacian.) => laplacian + procedure, non_overridable, private :: add_scalar_1D + procedure, non_overridable, private :: subtract_scalar_1D 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, pass(rhs) :: premultiply_double + procedure, non_overridable, private :: postmultiply_double + procedure, non_overridable, private :: exponentiate + procedure, non_overridable, private :: grad + procedure, non_overridable, private :: laplacian end type interface scalar_1D_t @@ -100,6 +116,12 @@ 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_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 + type(scalar_1D_t) scalar_1D + end function + end interface type, extends(tensor_1D_t) :: vector_1D_t @@ -116,7 +138,7 @@ pure module function construct_1D_scalar_from_function(initializer, order, cells generic :: weights => gradient_1D_weights #endif procedure, non_overridable :: dA - procedure, non_overridable, pass(vector_1D) :: weighted_premultiply + procedure, non_overridable, private, pass(vector_1D) :: weighted_premultiply procedure, non_overridable, private :: div procedure, non_overridable, private :: dot_surface_normal procedure, non_overridable, private :: vector_1D_grid @@ -202,6 +224,24 @@ pure module function construct_from_components(tensor_1D, divergence_operator_1D interface + pure logical module function is_face_centered(self) + !! Result is .true. if the values are face-centered and .false. otherwise + implicit none + class(tensor_1D_t), intent(in) :: self + end function + + pure logical module function is_cell_centered(self) + !! Result is .true. if the values are cell-centered and .false. otherwise + implicit none + class(tensor_1D_t), intent(in) :: self + end function + + pure logical module function is_cell_centers_extended(self) + !! Result is .true. if the values are at cell centers + boundaries and .false. otherwise + implicit none + class(tensor_1D_t), intent(in) :: self + end function + pure module function dA(self) !! Result is the grid's discrete surface-area differential for use in surface integrals of the form !! .SS. (f .x. (v .dot. dA)) @@ -259,6 +299,20 @@ pure module function divergence_1D_values(self) result(cell_centered_values) double precision, allocatable :: cell_centered_values(:) end function + pure module function d_dx(self) result(dself_dx) + !! Result is a scalar_1D_t approximating the spatial derivative of "self" + implicit none + class(scalar_1D_t), intent(in) :: self + type(scalar_1D_t) dself_dx + end function + + pure module function d2_dx2(self) result(d2_self_dx2) + !! Result is a scalar_1D_t approximating the spatial second derivative of "self" + implicit none + class(scalar_1D_t), intent(in) :: self + type(scalar_1D_t) d2_self_dx2 + end function + pure module function grad(self) result(gradient_1D) !! Result is mimetic gradient of the scalar_1D_t "self" implicit none @@ -266,6 +320,54 @@ pure module function grad(self) result(gradient_1D) type(gradient_1D_t) gradient_1D end function + pure module function premultiply_double(lhs, rhs) result(lhs_x_rhs) + !! Result is the product of lhs and rhs + implicit none + double precision, intent(in) :: lhs + class(scalar_1D_t), intent(in) :: rhs + type(scalar_1D_t) lhs_x_rhs + end function + + pure module function postmultiply_double(lhs, rhs) result(lhs_x_rhs) + !! Result is the product of lhs and rhs + implicit none + class(scalar_1D_t), intent(in) :: lhs + double precision, intent(in) :: rhs + type(scalar_1D_t) lhs_x_rhs + end function + + pure module function add_scalar_1D(lhs, rhs) result(total) + !! Result is the sum 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) total + 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 + class(scalar_1D_t), intent(in) :: lhs + type(scalar_1D_t), intent(in) :: rhs + type(scalar_1D_t) difference + end function + + pure module function exponentiate(self, exponent) result(power) + !! Result is mimetic Laplacian of the scalar_1D_t "self" + implicit none + class(scalar_1D_t), intent(in) :: self + integer, intent(in) :: exponent + type(scalar_1D_t) power + end function + + pure module function divide_by_integer(self, denominator) result(ratio) + !! Result is scalar_1D_t "self" divided by the integer denominator + implicit none + class(scalar_1D_t), intent(in) :: self + integer, intent(in) :: denominator + type(scalar_1D_t) ratio + end function + pure module function laplacian(self) result(laplacian_1D) !! Result is mimetic Laplacian of the scalar_1D_t "self" implicit none diff --git a/src/formal/vector_1D_s.F90 b/src/formal/vector_1D_s.F90 index 5bd892e..8dc9e34 100644 --- a/src/formal/vector_1D_s.F90 +++ b/src/formal/vector_1D_s.F90 @@ -76,6 +76,7 @@ pure module function construct_1D_vector_from_function(initializer, order, cells #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_) +#if ASSERTIONS associate( & q => divergence_1D%weights() & ,dx => (self%x_max_ - self%x_min_)/self%cells_ & @@ -85,6 +86,7 @@ pure module function construct_1D_vector_from_function(initializer, order, cells 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 end associate end associate diff --git a/src/formal_m.f90 b/src/formal_m.f90 index 523cfa2..2fc9526 100644 --- a/src/formal_m.f90 +++ b/src/formal_m.f90 @@ -13,13 +13,19 @@ module formal_m ,gradient_1D_t & ! result of `.grad. s` for a scalar_1D_t s ,divergence_1D_t & ! result of `.div. v` for a vector_1D_t v ,laplacian_1D_t & ! result of `.laplacian. s` for a scalar_1D_t s - ,scalar_1D_initializer_i & ! scalar_1D_t initializer abstract interface - ,vector_1D_initializer_i ! vector_1D_t initializar abstract interface + ,scalar_1D_initializer_i & ! scalar_1D_t initializer abstract interface + ,vector_1D_initializer_i & ! vector_1D_t initializar abstract interface + ,d_dx & ! scalar_1D_t spatial derivative + ,d2_dx2 ! scalar_1D_t spatial derivative - use mimetic_operators_1D_m, only : & + 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 + use interpolator_1D_m, only : & + centers_to_faces_1D_t & ! 1D mimetic interpolator producing cell-centered values from face-centered values + ,faces_to_centers_1D_t ! 1D mimetic interpolator producing face-centered values from cell-centered values + implicit none end module formal_m diff --git a/test/driver.f90 b/test/driver.f90 index 5c12c1d..65809b8 100644 --- a/test/driver.f90 +++ b/test/driver.f90 @@ -2,11 +2,14 @@ ! Terms of use are as specified in LICENSE.txt program test_suite_driver + !! Run the test suite and report results use julienne_m, only : test_fixture_t, test_harness_t use gradient_operator_1D_test_m, only : gradient_operator_1D_test_t use divergence_operator_1D_test_m, only : divergence_operator_1D_test_t use laplacian_operator_1D_test_m, only : laplacian_operator_1D_test_t 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 implicit none associate(test_harness => test_harness_t([ & @@ -14,7 +17,9 @@ program test_suite_driver ,test_fixture_t(divergence_operator_1D_test_t()) & ,test_fixture_t(laplacian_operator_1D_test_t()) & ,test_fixture_t(integration_operators_1D_test_t()) & + ,test_fixture_t(interpolator_1D_test_t()) & + ,test_fixture_t(scalar_1D_test_t()) & ])) call test_harness%report_results end associate -end program test_suite_driver +end program test_suite_driver \ 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 e61710a..1f35875 100644 --- a/test/gradient_operator_1D_test_m.F90 +++ b/test/gradient_operator_1D_test_m.F90 @@ -145,7 +145,6 @@ pure function parabola(x) result(y) y = 7*x**2 + 3*x + 5 end function - function check_grad_parabola() result(test_diagnosis) type(test_diagnosis_t) test_diagnosis procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => parabola diff --git a/test/interpolator_1D_test_m.F90 b/test/interpolator_1D_test_m.F90 new file mode 100644 index 0000000..8fcda74 --- /dev/null +++ b/test/interpolator_1D_test_m.F90 @@ -0,0 +1,137 @@ +! Copyright (c) 2026, The Regents of the University of California +! Terms of use are as specified in LICENSE.txt + +module interpolator_1D_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 & + ,csv + use formal_m, only : & + centers_to_faces_1D_t, scalar_1D_t, scalar_1D_initializer_i & + ,faces_to_centers_1D_t, vector_1D_t, vector_1D_initializer_i + + implicit none + + type, extends(test_t) :: interpolator_1D_test_t + contains + procedure, nopass :: subject + procedure, nopass :: results + end type + + double precision, parameter :: tolerance = 1D-11 + +contains + + pure function subject() result(test_subject) + character(len=:), allocatable :: test_subject + test_subject = 'A 1D mimetic interpolator' + end function + + function results() result(test_results) + type(interpolator_1D_test_t) interpolator_1D_test + type(test_result_t), allocatable :: test_results(:) + + test_results = interpolator_1D_test%run([ & + test_description_t('estimating center values from face values with 2nd- & 4th-order interpolators', usher(check_centers_to_faces)) & + ,test_description_t('estimating face values from center values with 2nd- & 4th-order interpolators', usher(check_faces_to_centers)) & + ]) + end function + + pure function line(x) result(y) + double precision, intent(in) :: x(:) + double precision, allocatable :: y(:) + y = x + end function + + pure function cubic(x) result(y) + double precision, intent(in) :: x(:) + double precision, allocatable :: y(:) + y = 7*x**3 + 4*x**2 + x + 2 + end function + + function check_centers_to_faces() result(test_diagnosis) + type(test_diagnosis_t) test_diagnosis + procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => null() + procedure(vector_1D_initializer_i), pointer :: vector_1D_initializer => line + double precision, parameter :: x_min = 0D0, x_max = 20D0 + integer order, cells + + test_diagnosis = passing_test() + + do order = 2, 4, 2 + + select case(order) + case(2) + cells = 10 + scalar_1D_initializer => line + case(4) + cells = 20 + scalar_1D_initializer => cubic + case default + error stop "check_centers_to_faces (interpolator_1D_test_m) unsupported order of accuracy" + end select + + associate( & + scalar_1D => scalar_1D_t(scalar_1D_initializer, order=order, cells=cells, x_min=x_min, x_max=x_max) & + ,vector_1D => vector_1D_t(vector_1D_initializer, order=order, cells=cells, x_min=x_min, x_max=x_max) & + ,interpolator => centers_to_faces_1D_t(order=order, cells=cells, dx=(x_max - x_min)/cells) & + ) + associate(scalar_at_faces => interpolator%face_values(scalar_1D%values())) + test_diagnosis = test_diagnosis .also. .all. (scalar_at_faces .approximates. scalar_1D_initializer(vector_1D%grid()) .within. tolerance) & + // string_t(" for order ") // string_t(order) + end associate + end associate + + end do + end function + + function check_faces_to_centers() result(test_diagnosis) + type(test_diagnosis_t) test_diagnosis + procedure(vector_1D_initializer_i), pointer :: vector_1D_initializer => null() + procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => line + double precision, parameter :: x_min = 0D0, x_max = 20D0 + integer order, cells + + test_diagnosis = passing_test() + + do order = 2, 4, 2 + + select case(order) + case(2) + cells = 10 + vector_1D_initializer => line + case(4) + cells = 20 + vector_1D_initializer => cubic + case default + error stop "check_faces_to_centers (interpolator_1D_test_m) unsupported order of accuracy" + end select + + associate( & + vector_1D => vector_1D_t(vector_1D_initializer, order=order, cells=cells, x_min=x_min, x_max=x_max) & + ,scalar_1D => scalar_1D_t(scalar_1D_initializer, order=order, cells=cells, x_min=x_min, x_max=x_max) & + ,interpolator => faces_to_centers_1D_t(order=order, cells=cells, dx=(x_max - x_min)/cells) & + ) + associate(vector_at_centers_extended => interpolator%center_values(vector_1D%values())) + associate( centers_extended => scalar_1D%grid() ) + test_diagnosis = test_diagnosis .also. & + .all. (vector_at_centers_extended .approximates. vector_1D_initializer(centers_extended) .within. tolerance) & + // string_t(" for order ") // string_t(order) + end associate + end associate + end associate + + end do + end function + +end module interpolator_1D_test_m diff --git a/test/scalar_1D_test_m.F90 b/test/scalar_1D_test_m.F90 new file mode 100644 index 0000000..4b8edef --- /dev/null +++ b/test/scalar_1D_test_m.F90 @@ -0,0 +1,173 @@ +! Copyright (c) 2026, The Regents of the University of California +! Terms of use are as specified in LICENSE.txt + +module scalar_1D_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 & + ,csv + use formal_m, only : scalar_1D_t, scalar_1D_initializer_i, d_dx, d2_dx2 + + implicit none + + type, extends(test_t) :: scalar_1D_test_t + contains + procedure, nopass :: subject + procedure, nopass :: results + end type + + double precision, parameter :: tolerance = 1D-11 + +contains + + pure function subject() result(test_subject) + character(len=:), allocatable :: test_subject + test_subject = 'The scalar_1D_t derived type' + end function + + function results() result(test_results) + type(scalar_1D_test_t) scalar_1D_test + type(test_result_t), allocatable :: test_results(:) + + test_results = scalar_1D_test%run([ & + test_description_t('raising a 1D scalar field to a power', usher(check_exponentiation)) & + ,test_description_t('dividing a 1D scalar field by a constant', usher(check_divison_operator)) & + ,test_description_t('computing a 1D scalar field derivative at cell centers extended', usher(check_derivative)) & + ,test_description_t('computing a 1D scalar field second derivative at cell centers extended', usher(check_2nd_derivative)) & + ]) + end function + + pure function line(x) result(y) + double precision, intent(in) :: x(:) + double precision, allocatable :: y(:) + y = x + end function + + function check_exponentiation() result(test_diagnosis) + type(test_diagnosis_t) test_diagnosis + procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => null() + integer order + + scalar_1D_initializer => line + test_diagnosis = passing_test() + + 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 + + end function + + function check_divison_operator() result(test_diagnosis) + type(test_diagnosis_t) test_diagnosis + procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => null() + integer order + + scalar_1D_initializer => line + test_diagnosis = passing_test() + + 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 + + end function + + pure function parabola(x) result(y) + double precision, intent(in) :: x(:) + double precision, allocatable :: y(:) + y = 7*x**2 + 3*x + 5D0 + end function + + pure function d_parabola_dx(x) result(y) + double precision, intent(in) :: x(:) + double precision, allocatable :: y(:) + y = 14*x + 3D0 + end function + + pure function d2_parabola_dx2(x) result(y) + double precision, intent(in) :: x(:) + double precision, allocatable :: y(:) + allocate(y(size(x))) + y = 14D0 + end function + + function check_derivative() result(test_diagnosis) + type(test_diagnosis_t) test_diagnosis + procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => null() + integer order + + scalar_1D_initializer => parabola + test_diagnosis = passing_test() + + do order = 2, 4, 2 + associate(scalar_1D => scalar_1D_t(scalar_1D_initializer, order=order, cells=10, x_min=0D0, x_max=10D0) ) + associate( d_scalar_1D_dx => d_dx(scalar_1D)) + test_diagnosis = test_diagnosis .also. .all. & + (d_scalar_1D_dx%values() .approximates. d_parabola_dx(scalar_1D%grid()) .within. tolerance) & + // string_t(" for order ") // string_t(order) + end associate + end associate + end do + + end function + + function check_2nd_derivative() result(test_diagnosis) + type(test_diagnosis_t) test_diagnosis + procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => null() + integer order + + scalar_1D_initializer => parabola + test_diagnosis = passing_test() + + do order = 2, 4, 2 + associate(scalar_1D => scalar_1D_t(scalar_1D_initializer, order=order, cells=10, x_min=0D0, x_max=10D0) ) + associate(d2_scalar_1D_dx2 => d2_dx2(scalar_1D)) + test_diagnosis = test_diagnosis .also. .all. & + (d2_scalar_1D_dx2%values() .approximates. d2_parabola_dx2(scalar_1D%grid()) .within. tolerance) & + // string_t(" for order ") // string_t(order) + end associate + end associate + end do + + end function + +end module scalar_1D_test_m