diff --git a/example/burgers-1D.F90 b/example/burgers-1D.F90 new file mode 100644 index 0000000..66ac30c --- /dev/null +++ b/example/burgers-1D.F90 @@ -0,0 +1,75 @@ +! 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(:) + double precision, parameter :: pi = acos(-1D0) + initial_condition = sin(2*pi*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 + !! Advance the 1D Burgers partial differential equation over time. + use initial_condition_m, only : initial_condition + use julienne_m, only : command_line_t + use formal_m, only : vector_1D_t, vector_1D_initializer_i + implicit none + + procedure(vector_1D_initializer_i), pointer :: vector_1D_initializer => initial_condition + character(len=:), allocatable :: order_string + type(command_line_t) command_line + type(vector_1D_t) :: u + 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 = 4 + else + read(order_string,"(i1)") order + end if + + print *, "order = ", order + + u = vector_1D_t(vector_1D_initializer, order, x_min=0D0, x_max=1D0, cells=50) + + block + double precision dt + dt = 1D0 + !associate(dt_grad_u2_2 => dt * (.div. ((u*u)))) + associate(dt_div_uu => .div. (u*u)) + end associate + end block + +#ifdef __GFORTRAN__ + stop +#endif + +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/dyad_1D_s.F90 b/src/formal/dyad_1D_s.F90 new file mode 100644 index 0000000..50c3ed5 --- /dev/null +++ b/src/formal/dyad_1D_s.F90 @@ -0,0 +1,66 @@ +! 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) dyad_1D_s + use julienne_m, only : & + call_julienne_assert_ & + ,operator(.all.) & + ,operator(.approximates.) & + ,operator(.equalsExpected.) & + ,operator(.within.) + use interpolator_1D_m, only : centers_to_faces_1D_t + implicit none + + double precision, parameter :: double_equivalence = 2D-4 + +contains + + module procedure dyad_over_integer + ratio%tensor_1D_t = tensor_1D_t(self%values_/numerator, self%x_min_, self%x_max_, self%cells_, order = self%order_) + ratio%divergence_operator_1D_ = self%divergence_operator_1D_ + end procedure + + module procedure construct_1D_dyad_from_components + dyad_1D%tensor_1D_t = tensor_1D + dyad_1D%divergence_operator_1D_ = divergence_operator_1D + end procedure + + module procedure div_dyad + + integer center + +#ifdef NAGFOR + associate(D => self%divergence_operator_1D_) +#else + associate(D => (self%divergence_operator_1D_)) +#endif + associate( & + Dv => D .x. self%values_ & + ,dx => (self%x_max_ - self%x_min_)/self%cells_ & + ) + associate(interpolator => centers_to_faces_1D_t(order=self%order_, cells=self%cells_, dx=dx)) + vector_1D = vector_1D_t( & + tensor_1D_t(interpolator%face_values(Dv(2:size(Dv)-1)), self%x_min_, self%x_max_, self%cells_, self%order_) & + ,divergence_operator_1D_t(k=self%order_, dx=dx, cells=self%cells_) & + ) + end associate +#if ASSERTIONS + associate(divergence_1D => divergence_1D_t(tensor_1D_t(Dv(2:size(Dv)-1), self%x_min_, self%x_max_, self%cells_, self%order_))) + associate( & + q => divergence_1D%weights() & + ,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 + end associate +#endif + end associate + end associate + + end procedure + +end submodule dyad_1D_s 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/tensors_1D_m.F90 b/src/formal/tensors_1D_m.F90 index 2615a7c..39ccc0e 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 @@ -16,6 +16,7 @@ module tensors_1D_m public :: scalar_1D_t public :: vector_1D_t + public :: dyad_1D_t public :: gradient_1D_t public :: laplacian_1D_t public :: divergence_1D_t @@ -110,19 +111,41 @@ pure module function construct_1D_scalar_from_function(initializer, order, cells generic :: operator(.x.) => weighted_premultiply generic :: operator(.div.) => div generic :: operator(.dot.) => dot_surface_normal + generic :: operator(*) => outer_product generic :: grid => vector_1D_grid generic :: values => vector_1D_values #ifdef __INTEL_COMPILER 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 :: outer_product procedure, non_overridable, private :: div procedure, non_overridable, private :: dot_surface_normal procedure, non_overridable, private :: vector_1D_grid procedure, non_overridable, private :: vector_1D_values end type + type, extends(tensor_1D_t) :: dyad_1D_t + type(divergence_operator_1D_t) divergence_operator_1D_ + contains + generic :: operator(.div.) => div_dyad + generic :: operator(/) => dyad_over_integer + procedure, non_overridable, private :: dyad_over_integer + procedure, non_overridable, private :: div_dyad + end type + + interface dyad_1D_t + + pure module function construct_1D_dyad_from_components(tensor_1D, divergence_operator_1D) result(dyad_1D) + !! Result is a 1D dyadic tensor 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(dyad_1D_t) dyad_1D + end function + + end interface + type, extends(tensor_1D_t) :: weighted_product_1D_t contains generic :: operator(.SS.) => surface_integrate_vector_x_scalar_1D @@ -202,6 +225,14 @@ pure module function construct_from_components(tensor_1D, divergence_operator_1D interface + pure module function dyad_over_integer(self, numerator) result(ratio) + !! Result is the dyad divided by the numerator + implicit none + class(dyad_1D_t), intent(in) :: self + integer, intent(in) :: numerator + type(dyad_1D_t) ratio + 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)) @@ -280,6 +311,13 @@ pure module function reduced_order_boundary_depth(self) result(num_nodes) integer num_nodes end function + pure module function div_dyad(self) result(vector_1D) + !! Result is the vector divergence of the dyad_1D_t "self" + implicit none + class(dyad_1D_t), intent(in) :: self + type(vector_1D_t) vector_1D + end function + pure module function div(self) result(divergence_1D) !! Result is mimetic divergence of the vector_1D_t "self" implicit none @@ -334,6 +372,13 @@ pure module function weighted_premultiply(scalar_1D, vector_1D) result(weighted_ type(weighted_product_1D_t) weighted_product_1D end function + pure module function outer_product(lhs, rhs) result(dyad_1D) + !! Result is the dyadic product of two vector operands + implicit none + class(vector_1D_t), intent(in) :: lhs, rhs + type(dyad_1D_t) dyad_1D + end function + pure module function gradient_1D_weights(self) result(weights) !! Result is an array of quadrature coefficients that can be used to compute a weighted !! inner product of a vector_1D_t object and a gradient_1D_t object. diff --git a/src/formal/vector_1D_s.F90 b/src/formal/vector_1D_s.F90 index 5bd892e..f4f7159 100644 --- a/src/formal/vector_1D_s.F90 +++ b/src/formal/vector_1D_s.F90 @@ -22,6 +22,16 @@ contains + module procedure outer_product + + call_julienne_assert(.all. ([lhs%x_min_, lhs%x_max_] .approximates. [rhs%x_min_, rhs%x_max_] .within. 0D0)) + call_julienne_assert(.all. ([lhs%cells_, lhs%order_] .equalsExpected. [rhs%cells_, rhs%order_])) + + dyad_1D%tensor_1D_t = tensor_1D_t(lhs%values_ * rhs%values_, lhs%x_min_, lhs%x_max_, lhs%cells_, lhs%order_ ) + dyad_1D%divergence_operator_1D_ = divergence_operator_1D_t(k=lhs%order_, dx=(lhs%x_max_ - lhs%x_min_)/lhs%cells_, cells=lhs%cells_) + + end procedure + module procedure dot_surface_normal v_dot_dS%tensor_1D_t = tensor_1D_t(vector_1D%values_*dS, vector_1D%x_min_, vector_1D%x_max_, vector_1D%cells_, vector_1D%order_) v_dot_dS%divergence_operator_1D_ = vector_1D%divergence_operator_1D_ @@ -76,6 +86,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 +96,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..7f507d9 100644 --- a/src/formal_m.f90 +++ b/src/formal_m.f90 @@ -16,10 +16,14 @@ module formal_m ,scalar_1D_initializer_i & ! scalar_1D_t initializer abstract interface ,vector_1D_initializer_i ! vector_1D_t initializar abstract interface - 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..72454f5 100644 --- a/test/driver.f90 +++ b/test/driver.f90 @@ -2,11 +2,13 @@ ! 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 implicit none associate(test_harness => test_harness_t([ & @@ -14,7 +16,8 @@ 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()) & ])) 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/interpolator_1D_test_m.F90 b/test/interpolator_1D_test_m.F90 new file mode 100644 index 0000000..9262c9f --- /dev/null +++ b/test/interpolator_1D_test_m.F90 @@ -0,0 +1,143 @@ +! 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()) & + ,face_locations => vector_1D%grid() & + ) + test_diagnosis = test_diagnosis .also. .all. (scalar_at_faces .approximates. scalar_1D_initializer(face_locations) .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 => interpolator%center_values(vector_1D%values()) & + ,faces => vector_1D%grid() & + ) + associate( centers_extended => scalar_1D%grid() ) + test_diagnosis = test_diagnosis .also. & + .all. (vector_at_centers .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