From 18515675862a9d00a43ec92cdb0eadbae74c41df Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Thu, 15 Jan 2026 11:46:07 -0800 Subject: [PATCH 01/33] fix(vector_1D_t): mk weighted_premultiply private --- src/formal/tensors_1D_m.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/formal/tensors_1D_m.F90 b/src/formal/tensors_1D_m.F90 index 2615a7c..aa379cf 100644 --- a/src/formal/tensors_1D_m.F90 +++ b/src/formal/tensors_1D_m.F90 @@ -116,7 +116,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 From 146cc772e69c3f6224daaf49beb3656d2a0d2141 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Wed, 28 Jan 2026 12:35:47 -0800 Subject: [PATCH 02/33] feat(burgers-1D): define initial condition, order --- example/burgers-1D.F90 | 64 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 64 insertions(+) create mode 100644 example/burgers-1D.F90 diff --git a/example/burgers-1D.F90 b/example/burgers-1D.F90 new file mode 100644 index 0000000..b923d56 --- /dev/null +++ b/example/burgers-1D.F90 @@ -0,0 +1,64 @@ +! 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 \' // 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 + + order_string = command_line%flag_value("--order") + + print *,new_line('') + print *," Initial condition" + print *," =================" + call execute_command_line("grep 'initial_condition =' example/burgers-1D.F90 | grep -v execute_command", wait=.true.) + + if (len(order_string)==0) then + order = 4 + else + read(order_string,"(i)") order + end if + + print *,"order = ",order + +#ifdef __GFORTRAN__ + stop +#endif + +end program From bcb9fbef1bba916e99da3afac21092aec2978066 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Thu, 29 Jan 2026 23:12:11 -0800 Subject: [PATCH 03/33] feat(dyad_1D_t): type, operator(/), operator(*) This commit 1. Adds a dyad_1D_t type,, 2. Adds a vector_1D_t operator(/) with a RHS integer operand, 3. Adds operator(*) for vector_1D_t operands & a dyad_1D_t result, 4. Updates the burgers-1D example to compute u*u/2. --- example/burgers-1D.F90 | 12 +++++++++--- src/formal/dyad_1D_s.F90 | 12 ++++++++++++ src/formal/tensors_1D_m.F90 | 24 ++++++++++++++++++++++++ src/formal/vector_1D_s.F90 | 8 ++++++++ 4 files changed, 53 insertions(+), 3 deletions(-) create mode 100644 src/formal/dyad_1D_s.F90 diff --git a/example/burgers-1D.F90 b/example/burgers-1D.F90 index b923d56..8914f59 100644 --- a/example/burgers-1D.F90 +++ b/example/burgers-1D.F90 @@ -42,20 +42,26 @@ program burgers_1D // 'where square brackets indicate optional arguments and angular brackets indicate user input values.' // new_line('') end if - order_string = command_line%flag_value("--order") - print *,new_line('') + 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,"(i)") order end if - print *,"order = ",order + print *, "order = ", order + + u = vector_1D_t(vector_1D_initializer, order, x_min=0D0, x_max=1D0, cells=50) + associate(uu_2 => u*u/2) + end associate #ifdef __GFORTRAN__ stop diff --git a/src/formal/dyad_1D_s.F90 b/src/formal/dyad_1D_s.F90 new file mode 100644 index 0000000..73718bf --- /dev/null +++ b/src/formal/dyad_1D_s.F90 @@ -0,0 +1,12 @@ +! Copyright (c) 2026, The Regents of the University of California +! Terms of use are as specified in LICENSE.txt + +submodule(tensors_1D_m) dyad_1D_s + implicit none +contains + + module procedure dyad_over_integer + ratio%tensor_1D_t = tensor_1D_t(self%tensor_1D_t%values_/numerator, self%x_min_, self%x_max_, self%cells_, order = self%order_) + end procedure + +end submodule dyad_1D_s diff --git a/src/formal/tensors_1D_m.F90 b/src/formal/tensors_1D_m.F90 index aa379cf..b6e635d 100644 --- a/src/formal/tensors_1D_m.F90 +++ b/src/formal/tensors_1D_m.F90 @@ -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,6 +111,7 @@ 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 @@ -117,12 +119,19 @@ pure module function construct_1D_scalar_from_function(initializer, order, cells #endif procedure, non_overridable :: dA 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 + contains + generic :: operator(/) => dyad_over_integer + procedure, non_overridable, private :: dyad_over_integer + end type + type, extends(tensor_1D_t) :: weighted_product_1D_t contains generic :: operator(.SS.) => surface_integrate_vector_x_scalar_1D @@ -202,6 +211,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)) @@ -334,6 +351,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..3c06737 100644 --- a/src/formal/vector_1D_s.F90 +++ b/src/formal/vector_1D_s.F90 @@ -22,6 +22,14 @@ 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_ ) + 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_ From 82c1e4dd20a2c427bdbd918adddb653ebecc59fa Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Fri, 30 Jan 2026 00:09:45 -0800 Subject: [PATCH 04/33] feat(dyad_1D_t): add operator(.div.) This commit adds support expressions of the form .div. uu, where uu is a 1D dyadic tensor and the result is stored at cell centers. To Do: 1. Embed an interpolation operator into the divergence operator to interpolate the result back out to the cell faces. 2. Use the interpolated values to construct the vector_1D_t result. --- example/burgers-1D.F90 | 2 +- src/formal/dyad_1D_s.F90 | 38 +++++++++++++++++++++++++++++++++++++ src/formal/tensors_1D_m.F90 | 23 +++++++++++++++++++++- src/formal/vector_1D_s.F90 | 10 +++++++--- 4 files changed, 68 insertions(+), 5 deletions(-) diff --git a/example/burgers-1D.F90 b/example/burgers-1D.F90 index 8914f59..1e2bde9 100644 --- a/example/burgers-1D.F90 +++ b/example/burgers-1D.F90 @@ -60,7 +60,7 @@ program burgers_1D print *, "order = ", order u = vector_1D_t(vector_1D_initializer, order, x_min=0D0, x_max=1D0, cells=50) - associate(uu_2 => u*u/2) + associate(div_uu_2 => .div. (u*u/2)) ! result is at cell centers; To Do: interpolate to cell faces end associate #ifdef __GFORTRAN__ diff --git a/src/formal/dyad_1D_s.F90 b/src/formal/dyad_1D_s.F90 index 73718bf..7687ab4 100644 --- a/src/formal/dyad_1D_s.F90 +++ b/src/formal/dyad_1D_s.F90 @@ -1,12 +1,50 @@ ! 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_ implicit none contains module procedure dyad_over_integer ratio%tensor_1D_t = tensor_1D_t(self%tensor_1D_t%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 + call_julienne_assert(size(self%values_)) + 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_) + error stop "To Do: interploate the dyad's divergence to the cell faces and construct a vector_1D_t" +#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 + end associate + end associate + end procedure end submodule dyad_1D_s diff --git a/src/formal/tensors_1D_m.F90 b/src/formal/tensors_1D_m.F90 index b6e635d..27cec5c 100644 --- a/src/formal/tensors_1D_m.F90 +++ b/src/formal/tensors_1D_m.F90 @@ -127,11 +127,25 @@ pure module function construct_1D_scalar_from_function(initializer, order, cells end type type, extends(tensor_1D_t) :: dyad_1D_t + type(divergence_operator_1D_t) divergence_operator_1D_ contains - generic :: operator(/) => dyad_over_integer + 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 @@ -297,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(divergence_1D) + !! Result is mimetic divergence of the dyad_1D_t "self" + implicit none + class(dyad_1D_t), intent(in) :: self + type(divergence_1D_t) divergence_1D !! discrete divergence + end function + pure module function div(self) result(divergence_1D) !! Result is mimetic divergence of the vector_1D_t "self" implicit none diff --git a/src/formal/vector_1D_s.F90 b/src/formal/vector_1D_s.F90 index 3c06737..f4f7159 100644 --- a/src/formal/vector_1D_s.F90 +++ b/src/formal/vector_1D_s.F90 @@ -24,10 +24,12 @@ 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_])) + 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_) - dyad_1D%tensor_1D_t = tensor_1D_t(lhs%values_ * rhs%values_, lhs%x_min_, lhs%x_max_, lhs%cells_, lhs%order_ ) end procedure module procedure dot_surface_normal @@ -84,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_ & @@ -93,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 From b36d47a57f5e100640fa925ea5f36714736da2b2 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Thu, 19 Feb 2026 21:55:43 -0800 Subject: [PATCH 05/33] refac: rename mimetic_* -> differential_operator_* This commit prepares for disambiguation of two types of mimetic matrix abstractions: one for differential operators and one for interpolation operators. --- ... => differential_operator_matrix_1D_s.F90} | 10 +++++----- ..._m.F90 => differential_operators_1D_m.F90} | 20 +++++++++---------- src/formal/divergence_operator_1D_s.F90 | 4 ++-- src/formal/gradient_operator_1D_s.F90 | 4 ++-- src/formal/tensors_1D_m.F90 | 2 +- src/formal_m.f90 | 2 +- 6 files changed, 21 insertions(+), 21 deletions(-) rename src/formal/{mimetic_matrix_1D_s.F90 => differential_operator_matrix_1D_s.F90} (76%) rename src/formal/{mimetic_operators_1D_m.F90 => differential_operators_1D_m.F90} (91%) 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/tensors_1D_m.F90 b/src/formal/tensors_1D_m.F90 index 27cec5c..0a3527f 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 diff --git a/src/formal_m.f90 b/src/formal_m.f90 index 523cfa2..8128923 100644 --- a/src/formal_m.f90 +++ b/src/formal_m.f90 @@ -16,7 +16,7 @@ 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 From 6c351be6510d63bcbc1331ac550c886f35c726de Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Thu, 19 Feb 2026 23:27:26 -0800 Subject: [PATCH 06/33] feat(interpolation): add operators This commit adds interpolation operator matrix abstractions to go from cell centers to faces and vice versa based on the definitions in Dumett & Castillo (2022). https://doi.org/10.13140/RG.2.2.26630.14400 --- src/formal/interpolation_operators_1D_m.F90 | 54 ++++++++++++++++++++ src/formal/interpolation_operators_1D_s.F90 | 55 +++++++++++++++++++++ 2 files changed, 109 insertions(+) create mode 100644 src/formal/interpolation_operators_1D_m.F90 create mode 100644 src/formal/interpolation_operators_1D_s.F90 diff --git a/src/formal/interpolation_operators_1D_m.F90 b/src/formal/interpolation_operators_1D_m.F90 new file mode 100644 index 0000000..df8b156 --- /dev/null +++ b/src/formal/interpolation_operators_1D_m.F90 @@ -0,0 +1,54 @@ +! Copyright (c) 2026, The Regents of the University of California +! Terms of use are as specified in LICENSE.txt + +module interpolation_operators_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 :: cells_to_faces_t + public :: faces_to_cells_t + + type interpolation_operator_1D_t + !! Encapsulate a staggered-grid interpolation matrix with a corresponding matrix-vector product operator + private + integer order + double precision first_ + double precision, allocatable :: upper_(:,:) + double precision, allocatable :: inner_(:) + double precision, allocatable :: lower_(:,:) + double precision final_ + end type + + type, extends(interpolation_operator_1D_t) :: centers_to_faces_1D_t + end type + + type, extends(interpolation_operator_1D_t) :: faces_to_centers_1D_t + end type + + interface centers_to_faces_1D_t + + pure module function c2f_component_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_component_constructor(order, cells, dx) result(faces_to_centers_1D) + !! Construct faces-to-centers 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 + +end module interpolation_operators_1D_m diff --git a/src/formal/interpolation_operators_1D_s.F90 b/src/formal/interpolation_operators_1D_s.F90 new file mode 100644 index 0000000..05549ab --- /dev/null +++ b/src/formal/interpolation_operators_1D_s.F90 @@ -0,0 +1,55 @@ +! Copyright (c) 2026, The Regents of the University of California +! Terms of use are as specified in LICENSE.txt + +submodule(interpolation_operators_1D_m) interpolation_operators_1D_s + implicit none + +contains + + module procedure c2f_component_constructor + + centers_to_faces_1D%cells_ = cells + + select case(order) + case(2) + centers_to_faces_1D%first_ = 1D0 + centers_to_faces_1D%upper_ = reshape([double precision::], [0,3]) + centers_to_faces_1D%inner_ = [1D0,1D0]/2D0 + centers_to_faces_1D%lower_ = reshape([double precision::], [0,3]) + centers_to_faces_1D%final_ = 1D0 + 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 + + end procedure + + module procedure f2c_component_constructor + + faces_to_centers_1D%cells_ = cells + + select case(order) + case(2) + faces_to_centers_1D%first_ = 1D0 + faces_to_centers_1D%upper_ = reshape([double precision::], [0,3]) + faces_to_centers_1D%inner_ = [1D0,1D0]/2D0 + faces_to_centers_1D%lower_ = reshape([double precision::], [0,3]) + faces_to_centers_1D%final_ = 1D0 + 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 + + end procedure + +end submodule interpolation_operators_1D_s From 44ad05430576db544003fa6d99aa7cf8d1d4a581 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Mon, 23 Feb 2026 14:47:50 -0800 Subject: [PATCH 07/33] WIP: start defining divergence of dyad --- src/formal/tensors_1D_m.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/formal/tensors_1D_m.F90 b/src/formal/tensors_1D_m.F90 index 0a3527f..0249276 100644 --- a/src/formal/tensors_1D_m.F90 +++ b/src/formal/tensors_1D_m.F90 @@ -315,7 +315,7 @@ pure module function div_dyad(self) result(divergence_1D) !! Result is mimetic divergence of the dyad_1D_t "self" implicit none class(dyad_1D_t), intent(in) :: self - type(divergence_1D_t) divergence_1D !! discrete divergence + type(vector_1D_t) vector_1D end function pure module function div(self) result(divergence_1D) From bb9e304fa5b853f191aa1e3f98d6ab928e629471 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Sat, 14 Mar 2026 10:12:38 -0700 Subject: [PATCH 08/33] doc(README): merge upstream changes --- README.md | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/README.md b/README.md index 977e80b..fd10780 100644 --- a/README.md +++ b/README.md @@ -77,16 +77,16 @@ Building and testing GCC | `gfortran` | 13 | `fpm test --compiler gfortran --profile release --flag "-ffree-line-length-none"` Intel | `ifx` | 2025.1.2 | `FOR_COARRAY_NUM_IMAGES=1 fpm test --compiler ifx --flag "-fpp -O3 -coarray" --profile release` LFortran | `lfortran` | 0.60.0-421-ge2c448c79 | `fpm test --compiler lfortran --flag "--cpp --realloc-lhs-arrays"` - LLVM | `flang` | 20-21 | `fpm test --compiler flang --profile release + LLVM | `flang` | 20-22 | `fpm test --compiler flang --profile release` LLVM | `flang` | 19 | `fpm test --compiler flang --profile release --flag "-mmlir -allow-assumed-rank"` NAG | `nagfor` | 7.2 Build 7242 | `fpm test --compiler nagfor --flag "-O3 -fpp"` ### `fpm` versions before 0.13.0 -With LLVM 20-22, replace the above `flang` command with +With LLVM 20-22, replace the corresponding command above with ``` fpm test --compiler flang-new --flag "-O3" ``` -With LLVM 19, replace the above `flang` command with +With LLVM 19, replace the corresponding command above with ``` fpm test --compiler flang-new --flag "-O3 -mmlir -allow-assumed-rank" ``` From d3cb45313c71b196a78a531beba2e35cdc4f244d Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Sat, 14 Mar 2026 10:14:31 -0700 Subject: [PATCH 09/33] doc(README): resolve git-stash conflict --- README.md | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/README.md b/README.md index fd10780..977e80b 100644 --- a/README.md +++ b/README.md @@ -77,16 +77,16 @@ Building and testing GCC | `gfortran` | 13 | `fpm test --compiler gfortran --profile release --flag "-ffree-line-length-none"` Intel | `ifx` | 2025.1.2 | `FOR_COARRAY_NUM_IMAGES=1 fpm test --compiler ifx --flag "-fpp -O3 -coarray" --profile release` LFortran | `lfortran` | 0.60.0-421-ge2c448c79 | `fpm test --compiler lfortran --flag "--cpp --realloc-lhs-arrays"` - LLVM | `flang` | 20-22 | `fpm test --compiler flang --profile release` + LLVM | `flang` | 20-21 | `fpm test --compiler flang --profile release LLVM | `flang` | 19 | `fpm test --compiler flang --profile release --flag "-mmlir -allow-assumed-rank"` NAG | `nagfor` | 7.2 Build 7242 | `fpm test --compiler nagfor --flag "-O3 -fpp"` ### `fpm` versions before 0.13.0 -With LLVM 20-22, replace the corresponding command above with +With LLVM 20-22, replace the above `flang` command with ``` fpm test --compiler flang-new --flag "-O3" ``` -With LLVM 19, replace the corresponding command above with +With LLVM 19, replace the above `flang` command with ``` fpm test --compiler flang-new --flag "-O3 -mmlir -allow-assumed-rank" ``` From 2d17f8477307f957b44f14c52f1d2c5a556ea995 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Mon, 16 Mar 2026 01:47:12 -0700 Subject: [PATCH 10/33] feat(interpolate): .div.dyad maps centers to faces This commit adds a centers-to-faces interpolator and uses that interpolator when computing the divergence of a dyad. --- example/burgers-1D.F90 | 19 ++++++---- src/formal/dyad_1D_s.F90 | 31 +++++++++------ ...erators_1D_m.F90 => interpolator_1D_m.F90} | 38 +++++++++++++------ ...erators_1D_s.F90 => interpolator_1D_s.F90} | 33 ++++++++++++++-- src/formal/tensors_1D_m.F90 | 4 +- 5 files changed, 89 insertions(+), 36 deletions(-) rename src/formal/{interpolation_operators_1D_m.F90 => interpolator_1D_m.F90} (52%) rename src/formal/{interpolation_operators_1D_s.F90 => interpolator_1D_s.F90} (57%) diff --git a/example/burgers-1D.F90 b/example/burgers-1D.F90 index 1e2bde9..d7d7085 100644 --- a/example/burgers-1D.F90 +++ b/example/burgers-1D.F90 @@ -35,14 +35,13 @@ program burgers_1D stop new_line('') // new_line('') & // 'Usage:' // new_line('') & // ' fpm run \' // new_line('') & - // ' --example burgers-1D \' // new_line('') & - // ' --compiler flang-new \' // new_line('') & - // ' --flag "-O3" \' // new_line('') & - // ' -- [--help|-h] | [--order ]' // new_line('') // 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 *," =================" @@ -60,8 +59,14 @@ program burgers_1D print *, "order = ", order u = vector_1D_t(vector_1D_initializer, order, x_min=0D0, x_max=1D0, cells=50) - associate(div_uu_2 => .div. (u*u/2)) ! result is at cell centers; To Do: interpolate to cell faces - end associate + + 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 diff --git a/src/formal/dyad_1D_s.F90 b/src/formal/dyad_1D_s.F90 index 7687ab4..3b685e4 100644 --- a/src/formal/dyad_1D_s.F90 +++ b/src/formal/dyad_1D_s.F90 @@ -5,6 +5,7 @@ submodule(tensors_1D_m) dyad_1D_s use julienne_m, only : call_julienne_assert_ + use interpolator_1D_m, only : centers_to_faces_1D_t implicit none contains @@ -28,18 +29,26 @@ associate(D => (self%divergence_operator_1D_)) #endif call_julienne_assert(size(self%values_)) - 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_) - error stop "To Do: interploate the dyad's divergence to the cell faces and construct a vector_1D_t" + 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 .center. 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( & - 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) + 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 diff --git a/src/formal/interpolation_operators_1D_m.F90 b/src/formal/interpolator_1D_m.F90 similarity index 52% rename from src/formal/interpolation_operators_1D_m.F90 rename to src/formal/interpolator_1D_m.F90 index df8b156..38d85c0 100644 --- a/src/formal/interpolation_operators_1D_m.F90 +++ b/src/formal/interpolator_1D_m.F90 @@ -1,19 +1,18 @@ ! Copyright (c) 2026, The Regents of the University of California ! Terms of use are as specified in LICENSE.txt -module interpolation_operators_1D_m +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 :: cells_to_faces_t - public :: faces_to_cells_t + public :: centers_to_faces_1D_t - type interpolation_operator_1D_t + type, abstract :: interpolator_1D_t !! Encapsulate a staggered-grid interpolation matrix with a corresponding matrix-vector product operator private - integer order + integer order_, cells_, dx_ double precision first_ double precision, allocatable :: upper_(:,:) double precision, allocatable :: inner_(:) @@ -21,15 +20,18 @@ module interpolation_operators_1D_m double precision final_ end type - type, extends(interpolation_operator_1D_t) :: centers_to_faces_1D_t + type, extends(interpolator_1D_t) :: centers_to_faces_1D_t + contains + generic :: operator(.center.) => interpolate_to_faces + procedure, non_overridable, private :: interpolate_to_faces end type - type, extends(interpolation_operator_1D_t) :: faces_to_centers_1D_t + type, extends(interpolator_1D_t) :: faces_to_centers_1D_t end type interface centers_to_faces_1D_t - pure module function c2f_component_constructor(order, cells, dx) result(centers_to_faces_1D) + 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 @@ -41,14 +43,26 @@ pure module function c2f_component_constructor(order, cells, dx) result(centers_ interface faces_to_centers_1D_t - pure module function f2c_component_constructor(order, cells, dx) result(faces_to_centers_1D) - !! Construct faces-to-centers interpolation operator + 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 + type(centers_to_faces_1D_t) faces_to_centers_1D + end function + + end interface + + interface + + pure module function interpolate_to_faces(self, centers) 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(:) + double precision, allocatable :: faces(:) end function end interface -end module interpolation_operators_1D_m +end module interpolator_1D_m \ No newline at end of file diff --git a/src/formal/interpolation_operators_1D_s.F90 b/src/formal/interpolator_1D_s.F90 similarity index 57% rename from src/formal/interpolation_operators_1D_s.F90 rename to src/formal/interpolator_1D_s.F90 index 05549ab..ce5fe2a 100644 --- a/src/formal/interpolation_operators_1D_s.F90 +++ b/src/formal/interpolator_1D_s.F90 @@ -1,14 +1,19 @@ ! Copyright (c) 2026, The Regents of the University of California ! Terms of use are as specified in LICENSE.txt -submodule(interpolation_operators_1D_m) interpolation_operators_1D_s +#include "julienne-assert-macros.h" + +submodule(interpolator_1D_m) interpolator_1D_s + use julienne_m, only : call_julienne_assert_, operator(.equalsExpected.) implicit none contains - module procedure c2f_component_constructor + 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) @@ -27,11 +32,14 @@ error stop "c2f_component_constructor: unsupported order" end select + call_julienne_assert(shape(self%lower_) .equalsExpected. shape(self%upper_)) end procedure - module procedure f2c_component_constructor + 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) @@ -50,6 +58,23 @@ error stop "f2c_component_constructor: unsupported order" end select + call_julienne_assert(shape(self%lower_) .equalsExpected. shape(self%upper_)) + end procedure + + module procedure interpolate_to_faces + integer row + associate( & + N => size(centers) , 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) & + ) + faces = [ self%first_ * centers(1) & + , matmul(self%upper_ , centers(1:upper_cols)) & + ,[(dot_product(self%inner_ , centers(row-upper_rows+1 : row-upper_rows+inner_cols)), row = upper_rows+1, N-lower_rows-1)] & + , matmul(self%lower_ , centers(N-lower_cols+1:N)) & + , self%final_ * centers(N) & + ] + end associate end procedure -end submodule interpolation_operators_1D_s +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 0249276..39ccc0e 100644 --- a/src/formal/tensors_1D_m.F90 +++ b/src/formal/tensors_1D_m.F90 @@ -311,8 +311,8 @@ pure module function reduced_order_boundary_depth(self) result(num_nodes) integer num_nodes end function - pure module function div_dyad(self) result(divergence_1D) - !! Result is mimetic divergence of the dyad_1D_t "self" + 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 From 67c748bb2ab396b21c669c60000e8fadbe2deda4 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Thu, 2 Apr 2026 23:18:06 -0700 Subject: [PATCH 11/33] feat(interpolation): centers-extended to faces This commit adds the interpolate_1D_t type with a type-bound procedure that interpolates using values at cell centers plus boundaries to estimate values at faces using the 2nd- and 4th-order schemes of Dumett & Castillo (2022). --- src/formal/dyad_1D_s.F90 | 13 ++++- src/formal/interpolator_1D_m.F90 | 7 +-- src/formal/interpolator_1D_s.F90 | 38 +++++++------ src/formal_m.f90 | 3 + test/driver.f90 | 5 +- test/interpolator_1D_test_m.F90 | 97 ++++++++++++++++++++++++++++++++ 6 files changed, 139 insertions(+), 24 deletions(-) create mode 100644 test/interpolator_1D_test_m.F90 diff --git a/src/formal/dyad_1D_s.F90 b/src/formal/dyad_1D_s.F90 index 3b685e4..260ef3b 100644 --- a/src/formal/dyad_1D_s.F90 +++ b/src/formal/dyad_1D_s.F90 @@ -4,9 +4,17 @@ #include "julienne-assert-macros.h" submodule(tensors_1D_m) dyad_1D_s - use julienne_m, only : call_julienne_assert_ + 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 @@ -28,14 +36,13 @@ #else associate(D => (self%divergence_operator_1D_)) #endif - call_julienne_assert(size(self%values_)) 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 .center. Dv(2:size(Dv)-1), self%x_min_, self%x_max_, self%cells_, self%order_) & + 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 diff --git a/src/formal/interpolator_1D_m.F90 b/src/formal/interpolator_1D_m.F90 index 38d85c0..e8a2fb0 100644 --- a/src/formal/interpolator_1D_m.F90 +++ b/src/formal/interpolator_1D_m.F90 @@ -22,8 +22,7 @@ module interpolator_1D_m type, extends(interpolator_1D_t) :: centers_to_faces_1D_t contains - generic :: operator(.center.) => interpolate_to_faces - procedure, non_overridable, private :: interpolate_to_faces + procedure, non_overridable :: face_values end type type, extends(interpolator_1D_t) :: faces_to_centers_1D_t @@ -55,11 +54,11 @@ pure module function f2c_constructor(order, cells, dx) result(faces_to_centers_1 interface - pure module function interpolate_to_faces(self, centers) result(faces) + 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(:) + double precision, intent(in) :: centers_extended(:) double precision, allocatable :: faces(:) end function diff --git a/src/formal/interpolator_1D_s.F90 b/src/formal/interpolator_1D_s.F90 index ce5fe2a..2cc7559 100644 --- a/src/formal/interpolator_1D_s.F90 +++ b/src/formal/interpolator_1D_s.F90 @@ -4,7 +4,7 @@ #include "julienne-assert-macros.h" submodule(interpolator_1D_m) interpolator_1D_s - use julienne_m, only : call_julienne_assert_, operator(.equalsExpected.) + use julienne_m, only : call_julienne_assert_, operator(.all.), operator(.equalsExpected.) implicit none contains @@ -17,11 +17,11 @@ select case(order) case(2) - centers_to_faces_1D%first_ = 1D0 - centers_to_faces_1D%upper_ = reshape([double precision::], [0,3]) - centers_to_faces_1D%inner_ = [1D0,1D0]/2D0 - centers_to_faces_1D%lower_ = reshape([double precision::], [0,3]) - centers_to_faces_1D%final_ = 1D0 + 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 @@ -32,7 +32,7 @@ error stop "c2f_component_constructor: unsupported order" end select - call_julienne_assert(shape(self%lower_) .equalsExpected. shape(self%upper_)) + call_julienne_assert(.all. (shape(centers_to_faces_1D%lower_) .equalsExpected. shape(centers_to_faces_1D%upper_))) end procedure module procedure f2c_constructor @@ -58,22 +58,28 @@ error stop "f2c_component_constructor: unsupported order" end select - call_julienne_assert(shape(self%lower_) .equalsExpected. shape(self%upper_)) + call_julienne_assert(.all. (shape(faces_to_centers_1D%lower_) .equalsExpected. shape(faces_to_centers_1D%upper_))) end procedure - module procedure interpolate_to_faces + module procedure face_values integer row + integer, parameter :: end_point = 1 associate( & - N => size(centers) , inner_cols => size(self%inner_) & + 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) & ) - faces = [ self%first_ * centers(1) & - , matmul(self%upper_ , centers(1:upper_cols)) & - ,[(dot_product(self%inner_ , centers(row-upper_rows+1 : row-upper_rows+inner_cols)), row = upper_rows+1, N-lower_rows-1)] & - , matmul(self%lower_ , centers(N-lower_cols+1:N)) & - , self%final_ * centers(N) & - ] + 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 diff --git a/src/formal_m.f90 b/src/formal_m.f90 index 8128923..14e481f 100644 --- a/src/formal_m.f90 +++ b/src/formal_m.f90 @@ -20,6 +20,9 @@ module formal_m 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 estimator producing cell-centered values given facial 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..bc1360f --- /dev/null +++ b/test/interpolator_1D_test_m.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 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 + use formal_m, only : centers_to_faces_1D_t, scalar_1D_t, scalar_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 values at cell centers given face values with 2nd- & 4th-order interpolators', usher(check_centers_to_faces)) & + ]) + end function + + pure function line(x) result(y) + double precision, intent(in) :: x(:) + double precision, allocatable :: y(:) + y = 3*x + 5 + 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() + 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) & + ,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()) & + ,gradient => .grad. scalar_1D & + ) + associate(face_locations => gradient%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 associate + + end do + end function + +end module interpolator_1D_test_m \ No newline at end of file From 0c839bd6ece0d79c1f1ea0dfb3c74bc50bdd487f Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Sat, 4 Apr 2026 14:02:36 -0700 Subject: [PATCH 12/33] feat(interpolation): faces to centers-extended This commit adds the faces_to_centers_t type with a type-bound procedure that interpolates using values at the faces to estimate to estimate values at cell centers plus boundaries using the 2nd- and 4th-order schemes of Dumett & Castillo (2022). --- src/formal/interpolator_1D_m.F90 | 13 +++++++- src/formal/interpolator_1D_s.F90 | 32 +++++++++++++++--- src/formal_m.f90 | 3 +- test/interpolator_1D_test_m.F90 | 56 ++++++++++++++++++++++++++++---- 4 files changed, 91 insertions(+), 13 deletions(-) diff --git a/src/formal/interpolator_1D_m.F90 b/src/formal/interpolator_1D_m.F90 index e8a2fb0..537c6c4 100644 --- a/src/formal/interpolator_1D_m.F90 +++ b/src/formal/interpolator_1D_m.F90 @@ -8,6 +8,7 @@ module interpolator_1D_m 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 @@ -26,6 +27,8 @@ module interpolator_1D_m 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 @@ -47,7 +50,7 @@ pure module function f2c_constructor(order, cells, dx) result(faces_to_centers_1 implicit none integer, intent(in) :: order, cells double precision, intent(in) :: dx - type(centers_to_faces_1D_t) faces_to_centers_1D + type(faces_to_centers_1D_t) faces_to_centers_1D end function end interface @@ -62,6 +65,14 @@ pure module function face_values(self, centers_extended) result(faces) 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 index 2cc7559..d7e5592 100644 --- a/src/formal/interpolator_1D_s.F90 +++ b/src/formal/interpolator_1D_s.F90 @@ -43,11 +43,11 @@ select case(order) case(2) - faces_to_centers_1D%first_ = 1D0 - faces_to_centers_1D%upper_ = reshape([double precision::], [0,3]) - faces_to_centers_1D%inner_ = [1D0,1D0]/2D0 - faces_to_centers_1D%lower_ = reshape([double precision::], [0,3]) - faces_to_centers_1D%final_ = 1D0 + 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 @@ -83,4 +83,26 @@ 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_m.f90 b/src/formal_m.f90 index 14e481f..7f507d9 100644 --- a/src/formal_m.f90 +++ b/src/formal_m.f90 @@ -21,7 +21,8 @@ module formal_m ,divergence_operator_1D_t ! matrix operator defining a 1D mimetic divergence use interpolator_1D_m, only : & - centers_to_faces_1D_t ! 1D mimetic estimator producing cell-centered values given facial values + 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 diff --git a/test/interpolator_1D_test_m.F90 b/test/interpolator_1D_test_m.F90 index bc1360f..4cd32bd 100644 --- a/test/interpolator_1D_test_m.F90 +++ b/test/interpolator_1D_test_m.F90 @@ -14,8 +14,11 @@ module interpolator_1D_test_m ,test_diagnosis_t & ,test_result_t & ,test_t & - ,usher - use formal_m, only : centers_to_faces_1D_t, scalar_1D_t, scalar_1D_initializer_i + ,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 @@ -39,14 +42,15 @@ function results() result(test_results) type(test_result_t), allocatable :: test_results(:) test_results = interpolator_1D_test%run([ & - test_description_t('estimating values at cell centers given face values with 2nd- & 4th-order interpolators', usher(check_centers_to_faces)) & + 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 = 3*x + 5 + y = x end function pure function cubic(x) result(y) @@ -54,7 +58,6 @@ pure function cubic(x) result(y) 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() @@ -94,4 +97,45 @@ function check_centers_to_faces() result(test_diagnosis) end do end function -end module interpolator_1D_test_m \ No newline at end of file + function check_faces_to_centers() result(test_diagnosis) + type(test_diagnosis_t) test_diagnosis + procedure(vector_1D_initializer_i), pointer :: vector_1D_initializer => null() + 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) & + ,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()) & + ,divergence => .div. vector_1D & + ,faces => vector_1D%grid() & + ) + associate( centers_extended => [faces(lbound(faces,1)), divergence%grid(), faces(ubound(faces,1))] ) + 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 From 5d8110186ed936dc1cf281d38312661b193ca416 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Sun, 5 Apr 2026 16:14:57 -0600 Subject: [PATCH 13/33] fix(burgers-1D): work around gfortran issue --- example/burgers-1D.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/example/burgers-1D.F90 b/example/burgers-1D.F90 index d7d7085..66ac30c 100644 --- a/example/burgers-1D.F90 +++ b/example/burgers-1D.F90 @@ -53,7 +53,7 @@ program burgers_1D if (len(order_string)==0) then order = 4 else - read(order_string,"(i)") order + read(order_string,"(i1)") order end if print *, "order = ", order From 420c46f16eab35b105259a8c4841c4adb6ff9d4c Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Sun, 5 Apr 2026 16:29:55 -0600 Subject: [PATCH 14/33] fix(interpolator_test): work around gfortran issue --- test/interpolator_1D_test_m.F90 | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/test/interpolator_1D_test_m.F90 b/test/interpolator_1D_test_m.F90 index 4cd32bd..9262c9f 100644 --- a/test/interpolator_1D_test_m.F90 +++ b/test/interpolator_1D_test_m.F90 @@ -58,9 +58,11 @@ pure function cubic(x) result(y) 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 @@ -81,16 +83,15 @@ function check_centers_to_faces() result(test_diagnosis) 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()) & - ,gradient => .grad. scalar_1D & + ,face_locations => vector_1D%grid() & ) - associate(face_locations => gradient%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 associate @@ -100,6 +101,7 @@ function check_centers_to_faces() result(test_diagnosis) 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 @@ -120,14 +122,14 @@ function check_faces_to_centers() result(test_diagnosis) 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()) & - ,divergence => .div. vector_1D & ,faces => vector_1D%grid() & ) - associate( centers_extended => [faces(lbound(faces,1)), divergence%grid(), faces(ubound(faces,1))] ) + 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) From 7bb94ce52c0109ad4ca5c59ffa34d3847b9dabba Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Sun, 5 Apr 2026 17:09:41 -0600 Subject: [PATCH 15/33] fix(dyad_1D_s): work around lfortran issue --- src/formal/dyad_1D_s.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/formal/dyad_1D_s.F90 b/src/formal/dyad_1D_s.F90 index 260ef3b..50c3ed5 100644 --- a/src/formal/dyad_1D_s.F90 +++ b/src/formal/dyad_1D_s.F90 @@ -18,7 +18,7 @@ contains module procedure dyad_over_integer - ratio%tensor_1D_t = tensor_1D_t(self%tensor_1D_t%values_/numerator, self%x_min_, self%x_max_, self%cells_, order = self%order_) + 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 From 6f1531f00c4d842ac2aab7effc7263ecd36e38c9 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Fri, 13 Mar 2026 17:05:43 -0700 Subject: [PATCH 16/33] doc(README): fix test command for fpm version<0.13 --- README.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) 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 From 03886cc9a6f32a2779ec18190c8fb1c83038235a Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Sun, 5 Apr 2026 21:28:03 -0600 Subject: [PATCH 17/33] test(dyad): assert .div. dyad_1D_t has N+1 points The divergence of a dyad is a vector and therefore should be stored at cell faces. However, the mimetic divergence operator computes values only at cell centers. Therefore, the dyad_1D_t .div. operator result is interpolated to the facs from the centers. To Do ----- Fix an erroneous zero value at the right boundary in the burgers-1D example (and possibly fix a left-boundary zero that might be only coincidentally correct). --- example/burgers-1D.F90 | 11 ++++++----- src/formal/dyad_1D_s.F90 | 11 +++++++---- 2 files changed, 13 insertions(+), 9 deletions(-) diff --git a/example/burgers-1D.F90 b/example/burgers-1D.F90 index 66ac30c..40227d3 100644 --- a/example/burgers-1D.F90 +++ b/example/burgers-1D.F90 @@ -11,7 +11,7 @@ pure function initial_condition(x) double precision, intent(in) :: x(:) double precision, allocatable :: initial_condition(:) double precision, parameter :: pi = acos(-1D0) - initial_condition = sin(2*pi*x) + initial_condition = x/sqrt(2D0) ! 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 @@ -51,20 +51,21 @@ program burgers_1D order_string = command_line%flag_value("--order") if (len(order_string)==0) then - order = 4 + order = 2 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) + u = vector_1D_t(vector_1D_initializer, order, x_min=0D0, x_max=20D0, cells=10) block double precision dt dt = 1D0 - !associate(dt_grad_u2_2 => dt * (.div. ((u*u)))) - associate(dt_div_uu => .div. (u*u)) + !associate(dt_div_uu=> dt * (.div. ((u*u)))) + associate(div_uu => .div. (u*u)) + print *,div_uu%values() end associate end block diff --git a/src/formal/dyad_1D_s.F90 b/src/formal/dyad_1D_s.F90 index 50c3ed5..02bd068 100644 --- a/src/formal/dyad_1D_s.F90 +++ b/src/formal/dyad_1D_s.F90 @@ -41,10 +41,13 @@ ,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_) & - ) + associate(face_values => interpolator%face_values(Dv)) + vector_1D = vector_1D_t( & + tensor_1D_t(face_values, self%x_min_, self%x_max_, self%cells_, self%order_) & + ,divergence_operator_1D_t(k=self%order_, dx=dx, cells=self%cells_) & + ) + call_julienne_assert(size(face_values) .equalsExpected. self%cells_ + 1) + end associate 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_))) From 91e6ef0e6c53587300e91d085f87e65016674f6e Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Mon, 4 May 2026 14:31:11 -0700 Subject: [PATCH 18/33] feat(scalar): power op, new construct, grid check This commit enhances the scalar_1D_t derived type by adding 1. A component constructor with a parent tensor_1D_t dummy arg; 2. A type-bound expoentiation operator; and 3. Type-bound logical functions: a. is_cell_centered, b. is_face_centered, and c. is_cell_centers_extended; 4. A unit test exercising 1-2 and 3c. --- example/burgers-1D.F90 | 13 +++---- src/formal/scalar_1D_s.F90 | 21 +++++++++++ src/formal/tensor_1D_s.f90 | 15 ++++++++ src/formal/tensors_1D_m.F90 | 37 +++++++++++++++++++ test/driver.f90 | 2 + test/scalar_1D_test_m.F90 | 73 +++++++++++++++++++++++++++++++++++++ 6 files changed, 154 insertions(+), 7 deletions(-) create mode 100644 test/scalar_1D_test_m.F90 diff --git a/example/burgers-1D.F90 b/example/burgers-1D.F90 index 40227d3..c1bacbe 100644 --- a/example/burgers-1D.F90 +++ b/example/burgers-1D.F90 @@ -22,13 +22,13 @@ 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 + use formal_m, only : scalar_1D_t, scalar_1D_initializer_i implicit none - procedure(vector_1D_initializer_i), pointer :: vector_1D_initializer => initial_condition + procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => initial_condition character(len=:), allocatable :: order_string type(command_line_t) command_line - type(vector_1D_t) :: u + type(scalar_1D_t) :: u integer order if (command_line%argument_present([character(len=len("--help")) :: ("--help"), "-h"])) then @@ -58,14 +58,13 @@ program burgers_1D print *, "order = ", order - u = vector_1D_t(vector_1D_initializer, order, x_min=0D0, x_max=20D0, cells=10) + u = scalar_1D_t(scalar_1D_initializer, order, x_min=0D0, x_max=20D0, cells=10) block double precision dt dt = 1D0 - !associate(dt_div_uu=> dt * (.div. ((u*u)))) - associate(div_uu => .div. (u*u)) - print *,div_uu%values() + !associate(dt_du22_dx => dt * .d_dx. ((u**2)/2)) + associate(dt_du22_dx => u**2) end associate end block diff --git a/src/formal/scalar_1D_s.F90 b/src/formal/scalar_1D_s.F90 index 47fa585..701d9c1 100644 --- a/src/formal/scalar_1D_s.F90 +++ b/src/formal/scalar_1D_s.F90 @@ -65,6 +65,27 @@ pure function cell_center_locations(x_min, x_max, cells) result(x) #endif + module procedure exponentiate + power%tensor_1D_t = tensor_1D_t( & + values = self%values_**2, 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 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 39ccc0e..85261f4 100644 --- a/src/formal/tensors_1D_m.F90 +++ b/src/formal/tensors_1D_m.F90 @@ -52,6 +52,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 @@ -78,10 +81,12 @@ 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(**) => exponentiate generic :: operator(.grad.) => grad generic :: operator(.laplacian.) => laplacian generic :: grid => scalar_1D_grid generic :: values => scalar_1D_values + procedure, non_overridable, private :: exponentiate procedure, non_overridable, private :: grad procedure, non_overridable, private :: laplacian procedure, non_overridable, private :: scalar_1D_values @@ -101,6 +106,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 @@ -225,6 +236,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 dyad_over_integer(self, numerator) result(ratio) !! Result is the dyad divided by the numerator implicit none @@ -297,6 +326,14 @@ pure module function grad(self) result(gradient_1D) type(gradient_1D_t) gradient_1D 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 laplacian(self) result(laplacian_1D) !! Result is mimetic Laplacian of the scalar_1D_t "self" implicit none diff --git a/test/driver.f90 b/test/driver.f90 index 72454f5..65809b8 100644 --- a/test/driver.f90 +++ b/test/driver.f90 @@ -9,6 +9,7 @@ program test_suite_driver 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([ & @@ -17,6 +18,7 @@ program test_suite_driver ,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 diff --git a/test/scalar_1D_test_m.F90 b/test/scalar_1D_test_m.F90 new file mode 100644 index 0000000..5e39eaa --- /dev/null +++ b/test/scalar_1D_test_m.F90 @@ -0,0 +1,73 @@ +! 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 + + 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)) & + ]) + 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) ) + associate( scalar_1D_squared => scalar_1D**2 ) + test_diagnosis = test_diagnosis .also. .all. & + (scalar_1D_squared%values() .approximates. scalar_1D%values()**2 .within. tolerance) & + // string_t(" for order ") // string_t(order) + end associate + end associate + end do + + end function + +end module scalar_1D_test_m From dc493f5e10acd583eabc75e63b35c627fb443ff5 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Mon, 4 May 2026 14:49:57 -0700 Subject: [PATCH 19/33] feat(scalar_1D_t): type-bound division operator This commit adds 1. operator(/) for an integer denominator and 2. a corresponding unit test. --- example/burgers-1D.F90 | 2 +- src/formal/scalar_1D_s.F90 | 9 +++++++++ src/formal/tensors_1D_m.F90 | 18 ++++++++++++++---- test/scalar_1D_test_m.F90 | 21 +++++++++++++++++++++ 4 files changed, 45 insertions(+), 5 deletions(-) diff --git a/example/burgers-1D.F90 b/example/burgers-1D.F90 index c1bacbe..9ae9934 100644 --- a/example/burgers-1D.F90 +++ b/example/burgers-1D.F90 @@ -64,7 +64,7 @@ program burgers_1D double precision dt dt = 1D0 !associate(dt_du22_dx => dt * .d_dx. ((u**2)/2)) - associate(dt_du22_dx => u**2) + associate(dt_du22_dx => (u**2)/2) end associate end block diff --git a/src/formal/scalar_1D_s.F90 b/src/formal/scalar_1D_s.F90 index 701d9c1..937a777 100644 --- a/src/formal/scalar_1D_s.F90 +++ b/src/formal/scalar_1D_s.F90 @@ -65,6 +65,15 @@ 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_/2, 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 exponentiate power%tensor_1D_t = tensor_1D_t( & values = self%values_**2, x_min = self%x_min_, x_max = self%x_max_, cells = self%cells_, order = self%order_ & diff --git a/src/formal/tensors_1D_m.F90 b/src/formal/tensors_1D_m.F90 index 85261f4..91a2b32 100644 --- a/src/formal/tensors_1D_m.F90 +++ b/src/formal/tensors_1D_m.F90 @@ -81,16 +81,18 @@ 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 :: grid => scalar_1D_grid + generic :: values => scalar_1D_values + generic :: operator(/) => divide_by_integer generic :: operator(**) => exponentiate generic :: operator(.grad.) => grad generic :: operator(.laplacian.) => laplacian - generic :: grid => scalar_1D_grid - generic :: values => scalar_1D_values + 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 :: exponentiate procedure, non_overridable, private :: grad procedure, non_overridable, private :: laplacian - procedure, non_overridable, private :: scalar_1D_values - procedure, non_overridable, private :: scalar_1D_grid end type interface scalar_1D_t @@ -334,6 +336,14 @@ pure module function exponentiate(self, exponent) result(power) 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/test/scalar_1D_test_m.F90 b/test/scalar_1D_test_m.F90 index 5e39eaa..51e6987 100644 --- a/test/scalar_1D_test_m.F90 +++ b/test/scalar_1D_test_m.F90 @@ -41,6 +41,7 @@ function results() result(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)) & ]) end function @@ -70,4 +71,24 @@ function check_exponentiation() result(test_diagnosis) 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) ) + associate( scalar_1D_squared => scalar_1D/2 ) + test_diagnosis = test_diagnosis .also. .all. & + (scalar_1D_squared%values() .approximates. scalar_1D%values()/2 .within. tolerance) & + // string_t(" for order ") // string_t(order) + end associate + end associate + end do + + end function + end module scalar_1D_test_m From eae9489fad774f164cc49c1420faa9da18d71e58 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Mon, 4 May 2026 16:37:52 -0700 Subject: [PATCH 20/33] refac(interpolator test): rm intermediate vars --- test/interpolator_1D_test_m.F90 | 16 +++++----------- 1 file changed, 5 insertions(+), 11 deletions(-) diff --git a/test/interpolator_1D_test_m.F90 b/test/interpolator_1D_test_m.F90 index 9262c9f..8fcda74 100644 --- a/test/interpolator_1D_test_m.F90 +++ b/test/interpolator_1D_test_m.F90 @@ -86,12 +86,9 @@ function check_centers_to_faces() result(test_diagnosis) ,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) + 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 @@ -125,13 +122,10 @@ function check_faces_to_centers() result(test_diagnosis) ,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(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 .approximates. vector_1D_initializer(centers_extended) .within. tolerance) & + .all. (vector_at_centers_extended .approximates. vector_1D_initializer(centers_extended) .within. tolerance) & // string_t(" for order ") // string_t(order) end associate end associate From 41f867eb63f2c7368a1f4cfdb3ed94d9b430f21f Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Mon, 4 May 2026 17:10:53 -0700 Subject: [PATCH 21/33] feat(d_dx): spatial partial deriv of 1D scalar --- example/burgers-1D.F90 | 5 +++-- src/formal/scalar_1D_s.F90 | 20 +++++++++++++++++ src/formal/tensors_1D_m.F90 | 9 ++++++++ test/gradient_operator_1D_test_m.F90 | 1 - test/scalar_1D_test_m.F90 | 33 ++++++++++++++++++++++++++++ 5 files changed, 65 insertions(+), 3 deletions(-) diff --git a/example/burgers-1D.F90 b/example/burgers-1D.F90 index 9ae9934..47be739 100644 --- a/example/burgers-1D.F90 +++ b/example/burgers-1D.F90 @@ -63,8 +63,9 @@ program burgers_1D block double precision dt dt = 1D0 - !associate(dt_du22_dx => dt * .d_dx. ((u**2)/2)) - associate(dt_du22_dx => (u**2)/2) + ! u_next = u + dt * .d_dt. u + ! u_next = u + dt * (nu * .d2_dx2. u - .ddx. (u**2)/2) + associate( u_next => .ddx. ((u**2)/2) ) end associate end block diff --git a/src/formal/scalar_1D_s.F90 b/src/formal/scalar_1D_s.F90 index 937a777..713fb66 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 @@ -112,6 +113,25 @@ 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 laplacian #ifndef __GFORTRAN__ diff --git a/src/formal/tensors_1D_m.F90 b/src/formal/tensors_1D_m.F90 index 91a2b32..d66abc8 100644 --- a/src/formal/tensors_1D_m.F90 +++ b/src/formal/tensors_1D_m.F90 @@ -85,12 +85,14 @@ pure module function construct_1D_tensor_from_components(values, x_min, x_max, c generic :: values => scalar_1D_values generic :: operator(/) => divide_by_integer generic :: operator(**) => exponentiate + generic :: operator(.ddx.) => d_dx generic :: operator(.grad.) => grad generic :: operator(.laplacian.) => laplacian 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 :: exponentiate + procedure, non_overridable, private :: d_dx procedure, non_overridable, private :: grad procedure, non_overridable, private :: laplacian end type @@ -321,6 +323,13 @@ 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 the "self" + implicit none + class(scalar_1D_t), intent(in) :: self + type(scalar_1D_t) dself_dx + end function + pure module function grad(self) result(gradient_1D) !! Result is mimetic gradient of the scalar_1D_t "self" implicit none 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/scalar_1D_test_m.F90 b/test/scalar_1D_test_m.F90 index 51e6987..49a9843 100644 --- a/test/scalar_1D_test_m.F90 +++ b/test/scalar_1D_test_m.F90 @@ -42,6 +42,7 @@ function results() result(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)) & ]) end function @@ -91,4 +92,36 @@ function check_divison_operator() result(test_diagnosis) end function + pure function parabola(x) result(y) + double precision, intent(in) :: x(:) + double precision, allocatable :: y(:) + y = 7*x**2 + 3*x + 5 + end function + + pure function d_parabola_dx(x) result(y) + double precision, intent(in) :: x(:) + double precision, allocatable :: y(:) + y = 14*x + 3 + 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 => .ddx. 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 + end module scalar_1D_test_m From 394b1b2cc763b437c606906dcf3f8d12de8a2ff0 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Mon, 4 May 2026 18:02:39 -0700 Subject: [PATCH 22/33] feat(d2_dx2): scalar_1D spatial 2nd partial deriv --- example/burgers-1D.F90 | 10 ++++++---- src/formal/scalar_1D_s.F90 | 4 ++++ src/formal/tensors_1D_m.F90 | 13 ++++++++++--- src/formal_m.f90 | 6 ++++-- test/scalar_1D_test_m.F90 | 36 ++++++++++++++++++++++++++++++++---- 5 files changed, 56 insertions(+), 13 deletions(-) diff --git a/example/burgers-1D.F90 b/example/burgers-1D.F90 index 47be739..26d9448 100644 --- a/example/burgers-1D.F90 +++ b/example/burgers-1D.F90 @@ -22,7 +22,7 @@ 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 : scalar_1D_t, scalar_1D_initializer_i + use formal_m, only : scalar_1D_t, scalar_1D_initializer_i, d_dx, d2_dx2 implicit none procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => initial_condition @@ -63,9 +63,11 @@ program burgers_1D block double precision dt dt = 1D0 - ! u_next = u + dt * .d_dt. u - ! u_next = u + dt * (nu * .d2_dx2. u - .ddx. (u**2)/2) - associate( u_next => .ddx. ((u**2)/2) ) + ! u_next = u + dt * d_dt(u) + ! u_next = u + dt * (nu * d2_dx2(u) - d_dx((u**2)/2)) + associate(du_dx => d_dx((u**2)/2)) + associate(d2u_dx2 => d2_dx2(u)) + end associate end associate end block diff --git a/src/formal/scalar_1D_s.F90 b/src/formal/scalar_1D_s.F90 index 713fb66..e54a7b2 100644 --- a/src/formal/scalar_1D_s.F90 +++ b/src/formal/scalar_1D_s.F90 @@ -132,6 +132,10 @@ pure function cell_center_locations(x_min, x_max, cells) result(x) 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/tensors_1D_m.F90 b/src/formal/tensors_1D_m.F90 index d66abc8..fedeb10 100644 --- a/src/formal/tensors_1D_m.F90 +++ b/src/formal/tensors_1D_m.F90 @@ -22,6 +22,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 @@ -85,14 +87,12 @@ pure module function construct_1D_tensor_from_components(values, x_min, x_max, c generic :: values => scalar_1D_values generic :: operator(/) => divide_by_integer generic :: operator(**) => exponentiate - generic :: operator(.ddx.) => d_dx generic :: operator(.grad.) => grad generic :: operator(.laplacian.) => laplacian 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 :: exponentiate - procedure, non_overridable, private :: d_dx procedure, non_overridable, private :: grad procedure, non_overridable, private :: laplacian end type @@ -324,12 +324,19 @@ pure module function divergence_1D_values(self) result(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 the "self" + !! 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 diff --git a/src/formal_m.f90 b/src/formal_m.f90 index 7f507d9..2fc9526 100644 --- a/src/formal_m.f90 +++ b/src/formal_m.f90 @@ -13,8 +13,10 @@ 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 differential_operators_1D_m, only : & gradient_operator_1D_t & ! matrix operator defining a 1D mimetic gradient diff --git a/test/scalar_1D_test_m.F90 b/test/scalar_1D_test_m.F90 index 49a9843..44e444b 100644 --- a/test/scalar_1D_test_m.F90 +++ b/test/scalar_1D_test_m.F90 @@ -16,7 +16,7 @@ module scalar_1D_test_m ,test_t & ,usher & ,csv - use formal_m, only : scalar_1D_t, scalar_1D_initializer_i + use formal_m, only : scalar_1D_t, scalar_1D_initializer_i, d_dx, d2_dx2 implicit none @@ -43,6 +43,7 @@ function results() result(test_results) 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 @@ -95,13 +96,20 @@ function check_divison_operator() result(test_diagnosis) pure function parabola(x) result(y) double precision, intent(in) :: x(:) double precision, allocatable :: y(:) - y = 7*x**2 + 3*x + 5 + 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 + 3 + 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) @@ -114,7 +122,7 @@ function check_derivative() result(test_diagnosis) do order = 2, 4, 2 associate(scalar_1D => scalar_1D_t(scalar_1D_initializer, order=order, cells=10, x_min=0D0, x_max=10D0) ) - associate( d_scalar_1D_dx => .ddx. scalar_1D ) + 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) @@ -124,4 +132,24 @@ function check_derivative() result(test_diagnosis) 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 From 62fa003e3ef8b4ec504803d2cb8191dbaf818ae7 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Mon, 4 May 2026 18:13:38 -0700 Subject: [PATCH 23/33] feat(scalar_1D): operator(*) for premult by double --- example/burgers-1D.F90 | 5 +++-- src/formal/scalar_1D_s.F90 | 6 ++++++ src/formal/tensors_1D_m.F90 | 10 ++++++++++ 3 files changed, 19 insertions(+), 2 deletions(-) diff --git a/example/burgers-1D.F90 b/example/burgers-1D.F90 index 26d9448..5bf7e65 100644 --- a/example/burgers-1D.F90 +++ b/example/burgers-1D.F90 @@ -61,12 +61,13 @@ program burgers_1D u = scalar_1D_t(scalar_1D_initializer, order, x_min=0D0, x_max=20D0, cells=10) block - double precision dt + double precision dt, nu dt = 1D0 + nu = 1D0 ! u_next = u + dt * d_dt(u) ! u_next = u + dt * (nu * d2_dx2(u) - d_dx((u**2)/2)) associate(du_dx => d_dx((u**2)/2)) - associate(d2u_dx2 => d2_dx2(u)) + associate(nu_d2u_dx2 => nu*d2_dx2(u)) end associate end associate end block diff --git a/src/formal/scalar_1D_s.F90 b/src/formal/scalar_1D_s.F90 index e54a7b2..c50d60f 100644 --- a/src/formal/scalar_1D_s.F90 +++ b/src/formal/scalar_1D_s.F90 @@ -75,6 +75,12 @@ pure function cell_center_locations(x_min, x_max, cells) result(x) ) 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 exponentiate power%tensor_1D_t = tensor_1D_t( & values = self%values_**2, x_min = self%x_min_, x_max = self%x_max_, cells = self%cells_, order = self%order_ & diff --git a/src/formal/tensors_1D_m.F90 b/src/formal/tensors_1D_m.F90 index fedeb10..0f79492 100644 --- a/src/formal/tensors_1D_m.F90 +++ b/src/formal/tensors_1D_m.F90 @@ -86,12 +86,14 @@ pure module function construct_1D_tensor_from_components(values, x_min, x_max, c generic :: grid => scalar_1D_grid generic :: values => scalar_1D_values generic :: operator(/) => divide_by_integer + generic :: operator(*) => premultiply_double generic :: operator(**) => exponentiate generic :: operator(.grad.) => grad generic :: operator(.laplacian.) => laplacian 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 :: exponentiate procedure, non_overridable, private :: grad procedure, non_overridable, private :: laplacian @@ -344,6 +346,14 @@ 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 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 exponentiate(self, exponent) result(power) !! Result is mimetic Laplacian of the scalar_1D_t "self" implicit none From b86d9b4f489edb891d26ba3c3ec5c68ad8cfcc14 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Mon, 4 May 2026 18:22:15 -0700 Subject: [PATCH 24/33] feat(scalar_1D): operator(-) for scalar_1D lhs/rhs --- example/burgers-1D.F90 | 4 +--- src/formal/scalar_1D_s.F90 | 6 ++++++ src/formal/tensors_1D_m.F90 | 12 +++++++++++- 3 files changed, 18 insertions(+), 4 deletions(-) diff --git a/example/burgers-1D.F90 b/example/burgers-1D.F90 index 5bf7e65..a6f8f2f 100644 --- a/example/burgers-1D.F90 +++ b/example/burgers-1D.F90 @@ -66,9 +66,7 @@ program burgers_1D nu = 1D0 ! u_next = u + dt * d_dt(u) ! u_next = u + dt * (nu * d2_dx2(u) - d_dx((u**2)/2)) - associate(du_dx => d_dx((u**2)/2)) - associate(nu_d2u_dx2 => nu*d2_dx2(u)) - end associate + associate(du_dt => nu*d2_dx2(u) - d_dx((u**2)/2)) end associate end block diff --git a/src/formal/scalar_1D_s.F90 b/src/formal/scalar_1D_s.F90 index c50d60f..8384ef7 100644 --- a/src/formal/scalar_1D_s.F90 +++ b/src/formal/scalar_1D_s.F90 @@ -75,6 +75,12 @@ pure function cell_center_locations(x_min, x_max, cells) result(x) ) end procedure + module procedure subtract_scalar_1D + 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 premultiply_double lhs_x_rhs%gradient_operator_1D_ = rhs%gradient_operator_1D_ lhs_x_rhs%tensor_1D_t = & diff --git a/src/formal/tensors_1D_m.F90 b/src/formal/tensors_1D_m.F90 index 0f79492..2538722 100644 --- a/src/formal/tensors_1D_m.F90 +++ b/src/formal/tensors_1D_m.F90 @@ -85,11 +85,13 @@ pure module function construct_1D_tensor_from_components(values, x_min, x_max, c contains generic :: grid => scalar_1D_grid generic :: values => scalar_1D_values + generic :: operator(-) => subtract_scalar_1D generic :: operator(/) => divide_by_integer generic :: operator(*) => premultiply_double generic :: operator(**) => exponentiate generic :: operator(.grad.) => grad generic :: operator(.laplacian.) => laplacian + 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 @@ -347,13 +349,21 @@ pure module function grad(self) result(gradient_1D) end function pure module function premultiply_double(lhs, rhs) result(lhs_x_rhs) - !! Result is product of lhs and 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 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 From 3dac0642a43bb2f3313e4bc36f06a3bfc53008b9 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Mon, 4 May 2026 18:53:06 -0700 Subject: [PATCH 25/33] feat(scalar_1D): operator(+), operator(*) This commit implements 1. operator(+) for lhs/rhs scalar_1D_t 2. operator(*) to premultiply scalar_1D_t scalar by real double --- src/formal/scalar_1D_s.F90 | 12 ++++++++++++ src/formal/tensors_1D_m.F90 | 21 ++++++++++++++++++++- 2 files changed, 32 insertions(+), 1 deletion(-) diff --git a/src/formal/scalar_1D_s.F90 b/src/formal/scalar_1D_s.F90 index 8384ef7..ce2d084 100644 --- a/src/formal/scalar_1D_s.F90 +++ b/src/formal/scalar_1D_s.F90 @@ -81,12 +81,24 @@ pure function cell_center_locations(x_min, x_max, cells) result(x) 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 + 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_**2, x_min = self%x_min_, x_max = self%x_max_, cells = self%cells_, order = self%order_ & diff --git a/src/formal/tensors_1D_m.F90 b/src/formal/tensors_1D_m.F90 index 2538722..5129dd5 100644 --- a/src/formal/tensors_1D_m.F90 +++ b/src/formal/tensors_1D_m.F90 @@ -86,16 +86,19 @@ pure module function construct_1D_tensor_from_components(values, x_min, x_max, c generic :: grid => scalar_1D_grid generic :: values => scalar_1D_values generic :: operator(-) => subtract_scalar_1D + generic :: operator(+) => add_scalar_1D generic :: operator(/) => divide_by_integer - generic :: operator(*) => premultiply_double + 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 @@ -356,6 +359,22 @@ pure module function premultiply_double(lhs, rhs) result(lhs_x_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 From 4e2f5c76c5ad404de7cebe175ff485e83aac6904 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Mon, 4 May 2026 19:01:29 -0700 Subject: [PATCH 26/33] feat(example): 1st-draft Burgers eq. solver This commit completes one full times step of the 1D Burgers equation solver, using 2nd- or 4th-order mimetic discretization and 2nd-order Runge-Kutta time advancement. To Do ----- 1. Set time step based on stability criterion 3. Loop over more time steps 2. Compare to exact solution in Fig. 6 of Rouson, Xia & Xu (2011) --- example/burgers-1D.F90 | 37 ++++++++++++++++++++++++------------- 1 file changed, 24 insertions(+), 13 deletions(-) diff --git a/example/burgers-1D.F90 b/example/burgers-1D.F90 index a6f8f2f..7b7972e 100644 --- a/example/burgers-1D.F90 +++ b/example/burgers-1D.F90 @@ -10,8 +10,7 @@ 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 = x/sqrt(2D0) + 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 @@ -25,10 +24,8 @@ program burgers_1D use formal_m, only : scalar_1D_t, scalar_1D_initializer_i, d_dx, d2_dx2 implicit none - procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => initial_condition character(len=:), allocatable :: order_string type(command_line_t) command_line - type(scalar_1D_t) :: u integer order if (command_line%argument_present([character(len=len("--help")) :: ("--help"), "-h"])) then @@ -58,20 +55,34 @@ program burgers_1D print *, "order = ", order - u = scalar_1D_t(scalar_1D_initializer, order, x_min=0D0, x_max=20D0, cells=10) - block - double precision dt, nu - dt = 1D0 - nu = 1D0 - ! u_next = u + dt * d_dt(u) - ! u_next = u + dt * (nu * d2_dx2(u) - d_dx((u**2)/2)) - associate(du_dt => nu*d2_dx2(u) - d_dx((u**2)/2)) - end associate + procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => initial_condition + double precision, parameter :: dt = 1D0, pi = acos(-1D0) + type(scalar_1D_t) u + + u = scalar_1D_t(scalar_1D_initializer, order, x_min=0D0, x_max=2*pi, cells=20) + !dt = u%runge_kutta_2nd_step(nu ,grid_resolution) + + !do while (t<=t_final) ! 2nd-order Runge-Kutta: + associate(u_half => u + d_dt(u)*dt/2) ! first substep + u = u + d_dt(u_half)*dt ! second substep + end associate + !t = t + dt + !end do + 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 From e64aee1689a0435691c023885894d1673629e910 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Tue, 5 May 2026 13:08:07 -0700 Subject: [PATCH 27/33] chore(dyad_1D_t): remove unused type This commit removes the dyad_1D_t type, which existed to support a different approach to solving the Burgers equation than the approach ultimately adopted in example/burgers-1D.F90. --- src/formal/dyad_1D_s.F90 | 69 ------------------------------------- src/formal/tensors_1D_m.F90 | 45 ------------------------ src/formal/vector_1D_s.F90 | 10 ------ 3 files changed, 124 deletions(-) delete mode 100644 src/formal/dyad_1D_s.F90 diff --git a/src/formal/dyad_1D_s.F90 b/src/formal/dyad_1D_s.F90 deleted file mode 100644 index 02bd068..0000000 --- a/src/formal/dyad_1D_s.F90 +++ /dev/null @@ -1,69 +0,0 @@ -! 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)) - associate(face_values => interpolator%face_values(Dv)) - vector_1D = vector_1D_t( & - tensor_1D_t(face_values, self%x_min_, self%x_max_, self%cells_, self%order_) & - ,divergence_operator_1D_t(k=self%order_, dx=dx, cells=self%cells_) & - ) - call_julienne_assert(size(face_values) .equalsExpected. self%cells_ + 1) - end associate - 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/tensors_1D_m.F90 b/src/formal/tensors_1D_m.F90 index 5129dd5..eecd06f 100644 --- a/src/formal/tensors_1D_m.F90 +++ b/src/formal/tensors_1D_m.F90 @@ -16,7 +16,6 @@ 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 @@ -133,7 +132,6 @@ pure module function construct_1D_scalar_from_parent(tensor_1D) result(scalar_1D 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 @@ -141,33 +139,12 @@ pure module function construct_1D_scalar_from_parent(tensor_1D) result(scalar_1D #endif procedure, non_overridable :: dA 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 @@ -265,14 +242,6 @@ pure logical module function is_cell_centers_extended(self) class(tensor_1D_t), intent(in) :: self end function - 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)) @@ -413,13 +382,6 @@ 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 @@ -474,13 +436,6 @@ 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 f4f7159..8dc9e34 100644 --- a/src/formal/vector_1D_s.F90 +++ b/src/formal/vector_1D_s.F90 @@ -22,16 +22,6 @@ 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_ From 61f0da55e02b050b6491aa0c769714f4777f66f4 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Tue, 5 May 2026 13:24:08 -0700 Subject: [PATCH 28/33] fix(scalar_1D): operator(/) for integer denominator --- src/formal/scalar_1D_s.F90 | 2 +- test/scalar_1D_test_m.F90 | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/formal/scalar_1D_s.F90 b/src/formal/scalar_1D_s.F90 index ce2d084..660c936 100644 --- a/src/formal/scalar_1D_s.F90 +++ b/src/formal/scalar_1D_s.F90 @@ -68,7 +68,7 @@ pure function cell_center_locations(x_min, x_max, cells) result(x) module procedure divide_by_integer ratio%tensor_1D_t = tensor_1D_t( & - values = self%values_/2, x_min = self%x_min_, x_max = self%x_max_, cells = self%cells_, order = self%order_ & + 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_ & diff --git a/test/scalar_1D_test_m.F90 b/test/scalar_1D_test_m.F90 index 44e444b..4b1ed35 100644 --- a/test/scalar_1D_test_m.F90 +++ b/test/scalar_1D_test_m.F90 @@ -83,9 +83,9 @@ function check_divison_operator() result(test_diagnosis) do order = 2, 4, 2 associate(scalar_1D => scalar_1D_t(scalar_1D_initializer, order=order, cells=10, x_min=0D0, x_max=10D0) ) - associate( scalar_1D_squared => scalar_1D/2 ) + associate( ratio => scalar_1D/2 ) test_diagnosis = test_diagnosis .also. .all. & - (scalar_1D_squared%values() .approximates. scalar_1D%values()/2 .within. tolerance) & + (ratio%values() .approximates. scalar_1D%values()/2 .within. tolerance) & // string_t(" for order ") // string_t(order) end associate end associate From 35a90031beeb335c2e4be420c242e543b3d6dff3 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Tue, 5 May 2026 18:30:54 -0700 Subject: [PATCH 29/33] chore(scalar_1): assert preconditions in operators This commit verifies that the scalar_1D_t subtraction and addition type-bound operators or on the same grid by comparing the related components. --- src/formal/scalar_1D_s.F90 | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/formal/scalar_1D_s.F90 b/src/formal/scalar_1D_s.F90 index 660c936..2903115 100644 --- a/src/formal/scalar_1D_s.F90 +++ b/src/formal/scalar_1D_s.F90 @@ -76,12 +76,18 @@ pure function cell_center_locations(x_min, x_max, cells) result(x) 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_) From 6a25084ddab87ebbae46a2c4abfba6a235653d60 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Tue, 5 May 2026 18:45:51 -0700 Subject: [PATCH 30/33] fix(scalar_1D_t): exponentiate & better test --- src/formal/scalar_1D_s.F90 | 2 +- test/scalar_1D_test_m.F90 | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/formal/scalar_1D_s.F90 b/src/formal/scalar_1D_s.F90 index 2903115..645f0dd 100644 --- a/src/formal/scalar_1D_s.F90 +++ b/src/formal/scalar_1D_s.F90 @@ -107,7 +107,7 @@ pure function cell_center_locations(x_min, x_max, cells) result(x) module procedure exponentiate power%tensor_1D_t = tensor_1D_t( & - values = self%values_**2, x_min = self%x_min_, x_max = self%x_max_, cells = self%cells_, order = self%order_ & + 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_ & diff --git a/test/scalar_1D_test_m.F90 b/test/scalar_1D_test_m.F90 index 4b1ed35..dbc436b 100644 --- a/test/scalar_1D_test_m.F90 +++ b/test/scalar_1D_test_m.F90 @@ -63,9 +63,9 @@ function check_exponentiation() result(test_diagnosis) do order = 2, 4, 2 associate(scalar_1D => scalar_1D_t(scalar_1D_initializer, order=order, cells=10, x_min=0D0, x_max=10D0) ) - associate( scalar_1D_squared => scalar_1D**2 ) + associate(cube => scalar_1D**3 ) test_diagnosis = test_diagnosis .also. .all. & - (scalar_1D_squared%values() .approximates. scalar_1D%values()**2 .within. tolerance) & + (cube%values() .approximates. scalar_1D%values()**3 .within. tolerance) & // string_t(" for order ") // string_t(order) end associate end associate From 158ba0a81931eaac68a789e78fa100a2b67d3dd4 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Tue, 5 May 2026 19:27:57 -0700 Subject: [PATCH 31/33] fix(burgers-1D): work around gfortran issues --- example/burgers-1D.F90 | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/example/burgers-1D.F90 b/example/burgers-1D.F90 index 7b7972e..251bf20 100644 --- a/example/burgers-1D.F90 +++ b/example/burgers-1D.F90 @@ -24,6 +24,9 @@ program burgers_1D 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 @@ -56,10 +59,13 @@ program burgers_1D print *, "order = ", order block - procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => initial_condition +#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) !dt = u%runge_kutta_2nd_step(nu ,grid_resolution) @@ -67,7 +73,7 @@ program burgers_1D associate(u_half => u + d_dt(u)*dt/2) ! first substep u = u + d_dt(u_half)*dt ! second substep end associate - !t = t + dt + !t = t + dt !end do end block From 365a7f8751752497534fc8ae00dd5f495b741533 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Tue, 5 May 2026 19:35:17 -0700 Subject: [PATCH 32/33] fix(scalar_1D_test): work around gfortran issues --- test/scalar_1D_test_m.F90 | 26 ++++++++++++++++++++++---- 1 file changed, 22 insertions(+), 4 deletions(-) diff --git a/test/scalar_1D_test_m.F90 b/test/scalar_1D_test_m.F90 index dbc436b..4b8edef 100644 --- a/test/scalar_1D_test_m.F90 +++ b/test/scalar_1D_test_m.F90 @@ -63,11 +63,21 @@ function check_exponentiation() result(test_diagnosis) do order = 2, 4, 2 associate(scalar_1D => scalar_1D_t(scalar_1D_initializer, order=order, cells=10, x_min=0D0, x_max=10D0) ) +#ifndef __GFORTRAN__ associate(cube => scalar_1D**3 ) test_diagnosis = test_diagnosis .also. .all. & (cube%values() .approximates. scalar_1D%values()**3 .within. tolerance) & // string_t(" for order ") // string_t(order) end associate +#else + block + type(scalar_1D_t) cube + cube = scalar_1D**3 + test_diagnosis = test_diagnosis .also. .all. & + (cube%values() .approximates. scalar_1D%values()**3 .within. tolerance) & + // string_t(" for order ") // string_t(order) + end block +#endif end associate end do @@ -83,11 +93,19 @@ function check_divison_operator() result(test_diagnosis) do order = 2, 4, 2 associate(scalar_1D => scalar_1D_t(scalar_1D_initializer, order=order, cells=10, x_min=0D0, x_max=10D0) ) - associate( ratio => scalar_1D/2 ) - test_diagnosis = test_diagnosis .also. .all. & - (ratio%values() .approximates. scalar_1D%values()/2 .within. tolerance) & - // string_t(" for order ") // string_t(order) +#ifndef __GFORTRAN__ + associate( half => scalar_1D/2 ) + test_diagnosis = test_diagnosis .also. .all. (half%values() .approximates. scalar_1D%values()/2 .within. tolerance) & + // string_t(" for order ") // string_t(order) end associate +#else + block + type(scalar_1D_t) half + half = scalar_1D/2 + test_diagnosis = test_diagnosis .also. .all. (half%values() .approximates. scalar_1D%values()/2 .within. tolerance) & + // string_t(" for order ") // string_t(order) + end block +#endif end associate end do From ae59e8d261a1327b88da6aed6cd0b4e89a1fe3ff Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Wed, 6 May 2026 17:44:03 -0700 Subject: [PATCH 33/33] chore(burgers): tidy up This commit better documents the main program of example/burgers-1D.F90 and removes commented code. --- example/burgers-1D.F90 | 21 ++++++++++++--------- 1 file changed, 12 insertions(+), 9 deletions(-) diff --git a/example/burgers-1D.F90 b/example/burgers-1D.F90 index 251bf20..6002e1a 100644 --- a/example/burgers-1D.F90 +++ b/example/burgers-1D.F90 @@ -18,7 +18,13 @@ pure function initial_condition(x) end module program burgers_1D - !! Advance the 1D Burgers partial differential equation over time. + !! 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 @@ -67,14 +73,11 @@ program burgers_1D scalar_1D_initializer => initial_condition u = scalar_1D_t(scalar_1D_initializer, order, x_min=0D0, x_max=2*pi, cells=20) - !dt = u%runge_kutta_2nd_step(nu ,grid_resolution) - - !do while (t<=t_final) ! 2nd-order Runge-Kutta: - associate(u_half => u + d_dt(u)*dt/2) ! first substep - u = u + d_dt(u_half)*dt ! second substep - end associate - !t = t + dt - !end do + + 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