diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml index 1d5461ef..789bfbd8 100644 --- a/.github/workflows/build.yml +++ b/.github/workflows/build.yml @@ -152,7 +152,7 @@ jobs: echo "FPM_FC=flang-new" >> "$GITHUB_ENV" ; \ elif [[ "$FC" =~ "ifx" ]] ; then \ echo "FPM_FC=ifx" >> "$GITHUB_ENV" ; \ - echo "FFLAGS=-fpp -coarray $FFLAGS" >> "$GITHUB_ENV" ; \ + echo "FFLAGS=-fpp -coarray -coarray-num-images=1 $FFLAGS" >> "$GITHUB_ENV" ; \ : ls -al /opt/intel/oneapi/compiler/2025.*/bin/ ; \ if type -p icpx ; then \ echo "FPM_CC=icx" >> "$GITHUB_ENV" ; \ diff --git a/doc/fortran-classes.md b/doc/fortran-classes.md new file mode 100644 index 00000000..76e4ea30 --- /dev/null +++ b/doc/fortran-classes.md @@ -0,0 +1,51 @@ +MOLE Fortran Class Diagram +-------------------------- + +```mermaid + +%%{init: { 'theme':'neo', "class" : {"hideEmptyMembersBox": true} } }%% + +classDiagram + +class tensor_1D_t +class scalar_1D_t +class vector_1D_t +class divergence_1D_t +class laplacian_1D_t +class gradient_operator_1D_t +class divergence_operator_1D_t +class mimetic_matrix_1D_t + +tensor_1D_t <|-- scalar_1D_t : is a +tensor_1D_t <|-- vector_1D_t : is a +tensor_1D_t <|-- divergence_1D_t : is a +divergence_1D_t <|-- laplacian_1D_t : is a +mimetic_matrix_1D_t <|-- gradient_operator_1D_t : is a +mimetic_matrix_1D_t <|-- divergence_operator_1D_t : is a + +class scalar_1D_t{ + - gradient_operator_1D_ : gradient_operator_1D_t + + operator(.grad.) vector_1D_t + + operator(.laplacian.) scalar_1D_t +} + +class vector_1D_t{ + - divergence_operator_1D_ : divergence_operator_1D_t + + operator(.div.) divergence_1D_t +} + +class mimetic_matrix_1D_t{ + - upper_ :: double precision + - inner_ :: double precision + - lower_ :: double precision +} + +class gradient_operator_1D_t{ + + operator(.x.) double precision[] + + assemble() double precision[] "2D array" +} + +class divergence_operator_1D_t{ + + operator(.x.) double precision[] + + assemble() double precision[] "2D array" +} diff --git a/example/div-grad-laplacian-1D.F90 b/example/div-grad-laplacian-1D.F90 new file mode 100644 index 00000000..eafc1b49 --- /dev/null +++ b/example/div-grad-laplacian-1D.F90 @@ -0,0 +1,133 @@ +module functions_m + implicit none + +contains + + pure function f(x) + double precision, intent(in) :: x(:) + double precision, allocatable :: f(:) + f = (x**3)/6 + (x**2)/2 + 1 + end function + + double precision elemental function df_dx(x) + double precision, intent(in) :: x + df_dx = (x**2)/2 + x + end function + + double precision elemental function d2f_dx2(x) + double precision, intent(in) :: x + d2f_dx2 = x + 1 + end function + +end module functions_m + +program div_grad_laplacian_1D + use functions_m + use julienne_m, only : file_t, string_t, operator(.separatedBy.) + use mole_m, only : scalar_1D_t, scalar_1D_initializer_i +#ifdef __GFORTRAN__ + use mole_m, only : vector_1D_t, laplacian_1D_t +#endif + implicit none + + procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => f + + print *,new_line('') + print *," 2nd-order approximations" + print *," ========================" + + call output(order=2) + + print *,new_line('') + print *," 4th-order approximations" + print *," ========================" + + call output(order=4) + +contains + +#ifndef __GFORTRAN__ + + subroutine output(order) + integer, intent(in) :: order + + associate( s => scalar_1D_t(scalar_1D_initializer, order=order, cells=10, x_min=0D0, x_max=20D0)) + associate( grad_s => .grad. s & + ,laplacian_s => .laplacian. s) + associate( s_grid => s%grid() & + ,grad_s_grid => grad_s%grid() & + ,laplacian_s_grid => laplacian_s%grid()) + associate( s_table => tabulate( & + string_t([character(len=18)::"x", "f(x) exp" , "f(x) act" ]) & + ,s_grid, f(s_grid), s%values() & + ) & + ,grad_s_table => tabulate( & + string_t([character(len=18)::"x", ".grad. f exp" , ".grad. f act" ]) & + ,grad_s_grid, df_dx(grad_s_grid), grad_s%values() & + ) & + ,laplacian_s_table => tabulate( & + string_t([character(len=18)::"x", ".laplacian. f exp", ".laplacian. f act"]) & + ,laplacian_s_grid, d2f_dx2(laplacian_s_grid), laplacian_s%values()) & + ) + call s_table%write_lines() + call grad_s_table%write_lines() + call laplacian_s_table%write_lines() + end associate + end associate + end associate + end associate + end subroutine + +#else + + subroutine output(order) + integer, intent(in) :: order + + type(scalar_1D_t) s + type(vector_1D_t) grad_s + type(laplacian_1D_t) laplacian_s + type(file_t) s_table, grad_s_table, laplacian_s_table + double precision, allocatable,dimension(:) :: s_grid, grad_s_grid, laplacian_s_grid + + s = scalar_1D_t(scalar_1D_initializer, order=order, cells=10, x_min=0D0, x_max=20D0) + grad_s = .grad. s + laplacian_s = .laplacian. s + + s_grid = s%grid() + grad_s_grid = grad_s%grid() + laplacian_s_grid = laplacian_s%grid() + + s_table = tabulate( & + string_t([character(len=18)::"x", "f(x) exp." , "f(x) act." ]) & + ,s_grid, f(s_grid), s%values() & + ) + grad_s_table = tabulate( & + string_t([character(len=18)::"x", ".grad. f exp." , ".grad. f act." ]) & + ,grad_s_grid, df_dx(grad_s_grid), grad_s%values() & + ) + laplacian_s_table = tabulate( & + string_t([character(len=18)::"x", ".laplacian. f exp.", ".laplacian. f act."]) & + ,laplacian_s_grid, d2f_dx2(laplacian_s_grid), laplacian_s%values() & + ) + call s_table%write_lines() + call grad_s_table%write_lines() + call laplacian_s_table%write_lines() + end subroutine + +#endif + + pure function tabulate(headings, abscissa, expected, actual) result(file) + double precision, intent(in), dimension(:) :: abscissa, expected, actual + type(string_t), intent(in) :: headings(:) + type(file_t) file + integer line + + file = file_t([ & + string_t("") & + ,headings .separatedBy. " " & + ,string_t("----------------------------------------------------------") & + ,[( string_t(abscissa(line)) // " " // string_t(expected(line)) // " " // string_t(actual(line)), line = 1, size(abscissa))] & + ]) + end function + +end program diff --git a/example/print-assembled-1D-operators.f90 b/example/print-assembled-1D-operators.f90 new file mode 100644 index 00000000..cf93b91c --- /dev/null +++ b/example/print-assembled-1D-operators.f90 @@ -0,0 +1,78 @@ +program print_assembled_1D_operators + !! Print fully assembled memetic 1D gradient, divergence, and Laplacian matrices, + !! including the zero elements. + use julienne_m, only : operator(.csv.), string_t, command_line_t + use mimetic_operators_1D_m, only : gradient_operator_1D_t, divergence_operator_1D_t + implicit none + + type(command_line_t) command_line + integer row + + command_line_settings: & + associate( & + gradient => command_line%argument_present(["--grad" ]) & + ,divergence => command_line%argument_present(["--div" ]) & + ,order => command_line%flag_value("--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 print-assembled-1D-operators \' // new_line('') & + // ' --compiler flang-new \' // new_line('') & + // ' --flag "-O3" \' // new_line('') & + // ' -- [--help|-h] | [--grad] [--div] [--order ]' // new_line('') // new_line('') & + // 'where square brackets indicate optional arguments and angular brackets indicate user input values.' // new_line('') + end if + + default_usage: & + associate(print_all => .not. any([gradient, divergence, len(order)/=0])) + + if (print_all .or. (gradient .and. len(order)==0) .or. (gradient .and. order=="2")) call print_gradient_operator( k=2, dx=1D0, m=5) + if (print_all .or. (divergence .and. len(order)==0) .or. (divergence .and. order=="2")) call print_divergence_operator(k=2, dx=1D0, m=5) + if (print_all .or. (gradient .and. len(order)==0) .or. (gradient .and. order=="4")) call print_gradient_operator( k=4, dx=1D0, m=9) + if (print_all .or. (divergence .and. len(order)==0) .or. (divergence .and. order=="4")) call print_divergence_operator(k=4, dx=1D0, m=9) + + end associate default_usage + end associate command_line_settings + +contains + + subroutine print_gradient_operator(k, dx, m) + integer, intent(in) :: k, m + double precision, intent(in) :: dx + + print *, new_line(""), "Gradient operator: order = ", k, " | cells = ", m, " | dx = ", dx + + associate(grad_op => gradient_operator_1D_t(k, dx, cells=m)) + associate(G => grad_op%assemble()) + do row = 1, size(G,1) + associate(csv_row => .csv. string_t(G(row,:))) + print '(a)', csv_row%string() + end associate + end do + end associate + end associate + + end subroutine + + subroutine print_divergence_operator(k, dx, m) + integer, intent(in) :: k, m + double precision, intent(in) :: dx + + print *, new_line(""), "Divergence operator: order = ", k, " | cells = ", m, " | dx = ", dx + + associate(div_op => divergence_operator_1D_t(k, dx, cells=m)) + associate(D => div_op%assemble()) + do row = 1, size(D,1) + associate(csv_row => .csv. string_t(D(row,:))) + print '(a)', csv_row%string() + end associate + end do + end associate + end associate + + end subroutine + +end program diff --git a/fpm.toml b/fpm.toml index 1a36b7b6..b79772c6 100644 --- a/fpm.toml +++ b/fpm.toml @@ -4,7 +4,7 @@ name = "MOLE" armadillo-code = {git = "https://gitlab.com/rouson/armadillo-code.git", tag = "fpm"} [dev-dependencies] -julienne = {git = "https://github.com/berkeleylab/julienne.git", tag = "3.1.5"} +julienne = {git = "https://github.com/berkeleylab/julienne.git", tag = "3.3.0"} [install] library = true diff --git a/src/fortran/divergence_1D_s.F90 b/src/fortran/divergence_1D_s.F90 new file mode 100644 index 00000000..254948b2 --- /dev/null +++ b/src/fortran/divergence_1D_s.F90 @@ -0,0 +1,33 @@ +submodule(tensors_1D_m) divergence_1D_s + implicit none + +contains + +#ifdef __GFORTRAN__ + + pure function cell_center_locations(x_min, x_max, cells) result(x) + double precision, intent(in) :: x_min, x_max + integer, intent(in) :: cells + double precision, allocatable:: x(:) + integer cell + + associate(dx => (x_max - x_min)/cells) + x = x_min + dx/2. + [((cell-1)*dx, cell = 1, cells)] + end associate + end function + +#endif + + module procedure construct_from_tensor + divergence_1D%tensor_1D_t = tensor_1D + end procedure + + module procedure divergence_1D_values + cell_centered_values = self%values_ + end procedure + + module procedure divergence_1D_grid + cell_centers = cell_center_locations(self%x_min_, self%x_max_, self%cells_) + end procedure + +end submodule divergence_1D_s \ No newline at end of file diff --git a/src/fortran/divergence_operator_1D_s.F90 b/src/fortran/divergence_operator_1D_s.F90 new file mode 100644 index 00000000..ba019d9e --- /dev/null +++ b/src/fortran/divergence_operator_1D_s.F90 @@ -0,0 +1,187 @@ +#include "mole-language-support.F90" +#include "julienne-assert-macros.h" + +submodule(mimetic_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.) +#endif + implicit none +contains + +#ifdef __GFORTRAN__ + + pure function negate_and_flip(A) result(Ap) + !! Transform a mimetic matrix upper block into a lower block + double precision, intent(in) :: A(:,:) + double precision, allocatable :: Ap(:,:) + integer row, column + + allocate(Ap, mold=A) + + reverse_elements_within_rows_and_flip_sign: & + do concurrent(row = 1:size(Ap,1)) + Ap(row,:) = -A(row,size(A,2):1:-1) + end do reverse_elements_within_rows_and_flip_sign + + reverse_elements_within_columns: & + do concurrent(column = 1 : size(Ap,2)) + Ap(:,column) = Ap(size(Ap,1):1:-1,column) + end do reverse_elements_within_columns + + end function + +#endif + + module procedure construct_1D_divergence_operator + + double precision, allocatable :: Ap(:,:) + + call_julienne_assert(cells .isAtLeast. 2*k+1) + + associate(A => A_block(k,dx)) + if (size(A) /= 0) then + Ap = negate_and_flip(A) + 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%k_ = k + divergence_operator_1D%dx_ = dx + divergence_operator_1D%m_ = cells + end associate + + contains + + pure function A_block(k, dx) result(matrix_block) + !! Compute the upper block submatrix "A" of the Corbino & Castillo (2020) mimetic divergence operator + integer, intent(in) :: k + double precision, intent(in) :: dx + double precision, allocatable :: matrix_block(:,:) + + order_of_accuracy: & + select case(k) + case(2) + matrix_block = reshape([ double precision :: & + ], shape=[0,0]) + case(4) + matrix_block = reshape([ & + -11/12D0, 17/24D0, 3/8D0, -5/24D0, 1/24D0 & + ], shape=[1,5], order=[2,1]) / dx + case default + associate(string_k => string_t(k)) + error stop "A (divergence_operator_1D_s): unsupported order of accuracy: " // string_k%string() + end associate + end select order_of_accuracy + + end function + + pure function M(k, dx) result(row) + !! Compute the middle block submatrix "M" of the Corbino & Castillo (2020) mimetic divergence operator + integer, intent(in) :: k + double precision, intent(in) :: dx + double precision, allocatable :: row(:) + + order_of_accuracy: & + select case(k) + case(2) + row = [-1D0, 1D0]/ dx + case(4) + row = [1D0/24D0, -9D0/8D0, 9D0/8D0, -1D0/24D0] / dx + case default + associate(string_k => string_t(k)) + error stop "M (divergence_operator_1D_s): unsupported order of accuracy: " // string_k%string() + end associate + end select order_of_accuracy + + end function + + end procedure construct_1D_divergence_operator + + module procedure submatrix_A_rows + call_julienne_assert(allocated(self%upper_)) + rows = size(self%upper_,1) + end procedure + + module procedure divergence_matrix_multiply + + double precision, allocatable :: product_inner(:) + + associate( & + upper_rows => size(self%upper_,1) & + ,lower_rows => size(self%lower_,1) & + ) + associate( & + inner_rows => self%m_ - (upper_rows + lower_rows) & ! rows(A) + rows(M) + rows(A') + 2 rows of zeros == m + 2 (Corbino & Castillo, 2020) + ,inner_columns => size(self%inner_) & + ) + call_julienne_assert((size(vec) .equalsExpected. self%m_ + 1)) + allocate(product_inner(inner_rows)) + +#if HAVE_DO_CONCURRENT_TYPE_SPEC_SUPPORT && HAVE_LOCALITY_SPECIFIER_SUPPORT + do concurrent(integer :: row = 1 : inner_rows) default(none) shared(product_inner, self, vec, inner_columns) + product_inner(row) = dot_product(self%inner_, vec(row : row + inner_columns - 1)) + end do +#else + block + integer row + do concurrent(row = 1 : inner_rows) + product_inner(row) = dot_product(self%inner_, vec(row : row + inner_columns - 1)) + end do + end block +#endif + + end associate + end associate + + associate( & + upper_columns => size(self%upper_,2) & + ,lower_columns => size(self%lower_,2) & + ) + matvec_product = [ & + 0D0 & + ,matmul(self%upper_, vec(1 : upper_columns )) & + ,product_inner & + ,matmul(self%lower_, vec(size(vec) - lower_columns + 1 : )) & + ,0D0 & + ] + call_julienne_assert(size(matvec_product) .equalsExpected. self%m_ + 2) + end associate + + end procedure + + module procedure assemble_divergence + + associate(rows => self%m_ + 2, cols => self%m_ + 1) + + allocate(D(rows, cols)) + +#if HAVE_DO_CONCURRENT_TYPE_SPEC_SUPPORT && HAVE_LOCALITY_SPECIFIER_SUPPORT + do concurrent(integer :: col=1:cols) default(none) shared(D, self, cols) + D(:,col) = self .x. e(dir=col, length=cols) + end do +#else + block + integer col + do concurrent(col=1:cols) + D(:,col) = divergence_matrix_multiply(self, e(dir=col, length=rows)) + ! work around gfortran 13-14 on Ubuntu (https://github.com/rouson/mole/actions/runs/19996603010/job/57345130923) + end do + end block +#endif + end associate + + contains + + pure function e(dir, length) result(unit_vector) + !! Result is the dir-th column of the len x len identity matrix + integer, intent(in) :: dir, length + double precision :: unit_vector(length) + unit_vector(1:dir-1) = 0D0 + unit_vector(dir) = 1D0 + unit_vector(dir+1:) = 0D0 + end function + + end procedure + +end submodule divergence_operator_1D_s \ No newline at end of file diff --git a/src/fortran/gradient_1D_m.f90 b/src/fortran/gradient_1D_m.f90 deleted file mode 100644 index 62ed6a80..00000000 --- a/src/fortran/gradient_1D_m.f90 +++ /dev/null @@ -1,49 +0,0 @@ -module gradient_1D_m - !! Define an abstraction for the collection of points used to compute gradidents: - !! cell centers plus oundaries. - implicit none - - private - public :: gradient_1D_t - - type gradient_1D_t - !! Encapsulate gradient_1D values produced only by .grad. (no other constructors) - private - double precision, allocatable :: vector_1D_(:) !! gradient_1D values at cell faces (nodes in 1D) - double precision x_min_ !! domain lower boundary - double precision x_max_ !! domain upper boundary - integer cells_ !! number of grid cells spanning the domain - contains - procedure values - procedure faces - end type - - interface gradient_1D_t - - pure module function construct_from_components(face_centered_values, x_min, x_max, cells) result(gradient_1D) - !! Result is an object storing gradient_1Ds at cell faces - implicit none - double precision, intent(in) :: face_centered_values(:), x_min, x_max - integer, intent(in) :: cells - type(gradient_1D_t) gradient_1D - end function - - end interface - - interface - - pure module function faces(self) result(x) - implicit none - class(gradient_1D_t), intent(in) :: self - double precision, allocatable :: x(:) - end function - - pure module function values(self) result(gradients) - implicit none - class(gradient_1D_t), intent(in) :: self - double precision, allocatable :: gradients(:) - end function - - end interface - -end module gradient_1D_m \ No newline at end of file diff --git a/src/fortran/gradient_1D_s.F90 b/src/fortran/gradient_1D_s.F90 deleted file mode 100644 index 9cff5f90..00000000 --- a/src/fortran/gradient_1D_s.F90 +++ /dev/null @@ -1,25 +0,0 @@ -submodule(gradient_1D_m) gradient_1D_s - implicit none - -contains - - module procedure construct_from_components - gradient_1D%vector_1D_ = face_centered_values - gradient_1D%x_min_ = x_min - gradient_1D%x_max_ = x_max - gradient_1D%cells_ = cells - end procedure - - module procedure values - gradients = self%vector_1D_ - end procedure - - module procedure faces - integer cell - x = [ self%x_min_ & - ,self%x_min_ + [(cell*(self%x_max_ - self%x_min_)/self%cells_, cell = 1, self%cells_-1)] & - ,self%x_max_ & - ] - end procedure - -end submodule gradient_1D_s \ No newline at end of file diff --git a/src/fortran/gradient_operator_1D_s.F90 b/src/fortran/gradient_operator_1D_s.F90 index f5f41ab2..fc1ce136 100644 --- a/src/fortran/gradient_operator_1D_s.F90 +++ b/src/fortran/gradient_operator_1D_s.F90 @@ -1,7 +1,7 @@ #include "julienne-assert-macros.h" #include "mole-language-support.F90" -submodule(scalar_1D_m) gradient_operator_1D_s +submodule(mimetic_operators_1D_m) gradient_operator_1D_s use julienne_m, only : call_julienne_assert_, string_t #if ASSERTIONS use julienne_m, only : operator(.isAtLeast.) @@ -10,102 +10,160 @@ contains - module procedure construct_from_parameters +#ifdef __GFORTRAN__ - call_julienne_assert(m .isAtLeast. 2*k) + pure function negate_and_flip(A) result(Ap) + !! Transform a mimetic matrix upper block into a lower block + double precision, intent(in) :: A(:,:) + double precision, allocatable :: Ap(:,:) + integer row, column - gradient_operator_1D%mimetic_matrix_1D_ = mimetic_matrix_1D_t( & - corbino_castillo_A( k, dx) & - ,corbino_castillo_M( k, dx) & - ,corbino_castillo_Ap(k, dx) & - ) - gradient_operator_1D%k_ = k - gradient_operator_1D%dx_ = dx - gradient_operator_1D%m_ = m - end procedure + allocate(Ap, mold=A) - pure function corbino_castillo_A(k, dx) result(matrix_block) - integer, intent(in) :: k - double precision, intent(in) :: dx - double precision, allocatable :: matrix_block(:,:) - - order_of_accuracy: & - select case(k) - case(2) - matrix_block = reshape([-8D0/3D0, 3D0, -1D0/3D0] , shape=[1,3]) / dx - case(4) - matrix_block = reshape([ & - -352D0/105D0, 35D0/ 8D0, -35D0/24D0, 21D0/40D0, -5D0/ 56D0 & - , 16D0/105D0, -31D0/24D0, 29D0/24D0, -3D0/40D0, 1D0/168D0 & - ], shape=[2,5], order=[2,1]) / dx - case default - associate(string_k => string_t(k)) - error stop "corbino_castillo_A: unsupported order of accuracy: " // string_k%string() - end associate - end select order_of_accuracy + reverse_elements_within_rows_and_flip_sign: & + do concurrent(row = 1:size(Ap,1)) + Ap(row,:) = -A(row,size(A,2):1:-1) + end do reverse_elements_within_rows_and_flip_sign - end function - - pure function corbino_castillo_M(k, dx) result(row) - integer, intent(in) :: k - double precision, intent(in) :: dx - double precision, allocatable :: row(:) - - order_of_accuracy: & - select case(k) - case(2) - row = [-1D0, 1D0]/ dx - case(4) - row = [1D0/24D0, -9D0/8D0, 9D0/8D0, -1D0/24D0] / dx - case default - associate(string_k => string_t(k)) - error stop "corbino_castillo_A: unsupported order of accuracy: " // string_k%string() - end associate - end select order_of_accuracy + reverse_elements_within_columns: & + do concurrent(column = 1 : size(Ap,2)) + Ap(:,column) = Ap(size(Ap,1):1:-1,column) + end do reverse_elements_within_columns end function -#if HAVE_DO_CONCURRENT_TYPE_SPEC_SUPPORT && HAVE_LOCALITY_SPECIFIER_SUPPORT +#endif + + module procedure construct_1D_gradient_operator + + call_julienne_assert(cells .isAtLeast. 2*k) - pure function corbino_castillo_Ap(k, dx) result(matrix_block) - integer, intent(in) :: k - double precision, intent(in) :: dx - double precision, allocatable :: matrix_block(:,:) - - associate(A => corbino_castillo_A(k, dx)) - allocate(matrix_block , mold=A) - reverse_elements_within_rows_and_flip_sign: & - do concurrent(integer :: row = 1:size(matrix_block,1)) default(none) shared(matrix_block, A) - matrix_block(row,:) = -A(row,size(A,2):1:-1) - end do reverse_elements_within_rows_and_flip_sign - reverse_elements_within_columns: & - do concurrent(integer :: column = 1 : size(matrix_block,2)) default(none) shared(matrix_block) - matrix_block(:,column) = matrix_block(size(matrix_block,1):1:-1,column) - end do reverse_elements_within_columns + 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%k_ = k + gradient_operator_1D%dx_ = dx + gradient_operator_1D%m_ = cells end associate - end function - + + contains + + pure function corbino_castillo_A(k, dx) result(matrix_block) + integer, intent(in) :: k + double precision, intent(in) :: dx + double precision, allocatable :: matrix_block(:,:) + + order_of_accuracy: & + select case(k) + case(2) + matrix_block = reshape([-8D0/3D0, 3D0, -1D0/3D0] , shape=[1,3]) / dx + case(4) + matrix_block = reshape([ & + -352D0/105D0, 35D0/ 8D0, -35D0/24D0, 21D0/40D0, -5D0/ 56D0 & + , 16D0/105D0, -31D0/24D0, 29D0/24D0, -3D0/40D0, 1D0/168D0 & + ], shape=[2,5], order=[2,1]) / dx + case default + associate(string_k => string_t(k)) + error stop "corbino_castillo_A: unsupported order of accuracy: " // string_k%string() + end associate + end select order_of_accuracy + + end function + + pure function corbino_castillo_M(k, dx) result(row) + integer, intent(in) :: k + double precision, intent(in) :: dx + double precision, allocatable :: row(:) + + order_of_accuracy: & + select case(k) + case(2) + row = [-1D0, 1D0]/ dx + case(4) + row = [1D0/24D0, -9D0/8D0, 9D0/8D0, -1D0/24D0] / dx + case default + associate(string_k => string_t(k)) + error stop "corbino_castillo_A: unsupported order of accuracy: " // string_k%string() + end associate + end select order_of_accuracy + + end function + + end procedure construct_1D_gradient_operator + + module procedure gradient_matrix_multiply + + double precision, allocatable :: product_inner(:) + + associate( & + upper_rows => size(self%upper_,1) & + ,lower_rows => size(self%lower_,1) & + ) + associate( & + inner_rows => size(vec) - (upper_rows + lower_rows + 1) & + ,inner_columns => size(self%inner_) & + ) + allocate(product_inner(inner_rows)) + +#if HAVE_DO_CONCURRENT_TYPE_SPEC_SUPPORT && HAVE_LOCALITY_SPECIFIER_SUPPORT + do concurrent(integer :: row = 1 : inner_rows) default(none) shared(product_inner, self, vec, inner_columns) + product_inner(row) = dot_product(self%inner_, vec(row + 1 : row + inner_columns)) + end do #else + block + integer row + do concurrent(row = 1 : inner_rows) + product_inner(row) = dot_product(self%inner_, vec(row + 1 : row + inner_columns)) + end do + end block +#endif - pure function corbino_castillo_Ap(k, dx) result(matrix_block) - integer, intent(in) :: k - double precision, intent(in) :: dx - double precision, allocatable :: matrix_block(:,:) - integer row, column + end associate + end associate - associate(A => corbino_castillo_A(k, dx)) - allocate(matrix_block , mold=A) - reverse_elements_within_rows_and_flip_sign: & - do concurrent(row = 1:size(matrix_block,1)) - matrix_block(row,:) = -A(row,size(A,2):1:-1) - end do reverse_elements_within_rows_and_flip_sign - reverse_elements_within_columns: & - do concurrent(column = 1 : size(matrix_block,2)) - matrix_block(:,column) = matrix_block(size(matrix_block,1):1:-1,column) - end do reverse_elements_within_columns + associate( & + upper_columns => size(self%upper_,2) & + ,lower_columns => size(self%lower_,2) & + ) + matvec_product = [ & + matmul(self%upper_, vec(1 : upper_columns)) & + ,product_inner & + ,matmul(self%lower_, vec(size(vec) - lower_columns + 1 : )) & + ] end associate - end function - + end procedure + + module procedure assemble_gradient + + associate(rows => self%m_ + 1, cols => self%m_ + 2) + + allocate(G(rows, cols), source = 0D0) + +#if HAVE_DO_CONCURRENT_TYPE_SPEC_SUPPORT && HAVE_LOCALITY_SPECIFIER_SUPPORT + do concurrent(integer :: col=1:cols) default(none) shared(G, self, cols) + G(:,col) = self .x. e(dir=col, length=cols) + end do +#else + block + integer col + do concurrent(col=1:cols) + G(:,col) = gradient_matrix_multiply(self, e(dir=col, length=cols)) + !! Work around gfortran 13-14 issue: https://github.com/rouson/mole/actions/runs/19996724722/job/57345438334 + end do + end block #endif + end associate + + contains + + pure function e(dir, length) result(unit_vector) + !! Result is the dir-th column of the len x len identity matrix + integer, intent(in) :: dir, length + double precision :: unit_vector(length) + unit_vector(1:dir-1) = 0D0 + unit_vector(dir) = 1D0 + unit_vector(dir+1:) = 0D0 + end function + + end procedure end submodule gradient_operator_1D_s \ No newline at end of file diff --git a/src/fortran/laplacian_s.f90 b/src/fortran/laplacian_s.f90 new file mode 100644 index 00000000..e882ced7 --- /dev/null +++ b/src/fortran/laplacian_s.f90 @@ -0,0 +1,9 @@ +submodule(tensors_1D_m) laplacian_s + implicit none +contains + + module procedure reduced_order_boundary_depth + num_nodes = self%boundary_depth_ + end procedure + +end submodule laplacian_s diff --git a/src/fortran/mimetic_matrix_1D_s.F90 b/src/fortran/mimetic_matrix_1D_s.F90 index 1138d8bc..d1881a80 100644 --- a/src/fortran/mimetic_matrix_1D_s.F90 +++ b/src/fortran/mimetic_matrix_1D_s.F90 @@ -1,69 +1,15 @@ -#include "mole-language-support.F90" -#include "julienne-assert-macros.h" - -submodule(scalar_1D_m) mimetic_matrix_1D_s - use julienne_m, only : call_julienne_assert_, string_t, operator(.equalsExpected.), operator(.csv.) +submodule(mimetic_operators_1D_m) mimetic_matrix_1D_s + use julienne_m, only : string_t, operator(.csv.) implicit none contains - module procedure construct_from_components + module procedure construct_matrix_operator mimetic_matrix_1D%upper_ = upper mimetic_matrix_1D%inner_ = inner mimetic_matrix_1D%lower_ = lower end procedure -#if HAVE_DO_CONCURRENT_TYPE_SPEC_SUPPORT && HAVE_LOCALITY_SPECIFIER_SUPPORT - - module procedure matvec - - double precision, allocatable :: product_inner(:) - - associate(upper => size(self%upper_,1), lower => size(self%lower_,1)) - associate(inner_rows => size(scalar_1D%scalar_1D_) - (upper + lower + 1)) - - allocate(product_inner(inner_rows)) - - do concurrent(integer :: row = 1 : inner_rows) default(none) shared(product_inner, self, scalar_1D) - product_inner(row) = dot_product(self%inner_, scalar_1D%scalar_1D_(row + 1 : row + size(self%inner_))) - end do - - matvec_product = [ & - matmul(self%upper_, scalar_1D%scalar_1D_(1 : size(self%upper_,2))) & - ,product_inner & - ,matmul(self%lower_, scalar_1D%scalar_1D_(size(scalar_1D%scalar_1D_) - size(self%lower_,2) + 1 : )) & - ] - end associate - end associate - end procedure - -#else - - module procedure matvec - - integer row - double precision, allocatable :: product_inner(:) - - associate(upper => size(self%upper_,1), lower => size(self%lower_,1)) - associate(inner_rows => size(scalar_1D%scalar_1D_) - (upper + lower + 1)) - - allocate(product_inner(inner_rows)) - - do concurrent(row = 1 : inner_rows) - product_inner(row) = dot_product(self%inner_, scalar_1D%scalar_1D_(row + 1 : row + size(self%inner_))) - end do - - matvec_product = [ & - matmul(self%upper_, scalar_1D%scalar_1D_(1 : size(self%upper_,2))) & - ,product_inner & - ,matmul(self%lower_, scalar_1D%scalar_1D_(size(scalar_1D%scalar_1D_) - size(self%lower_,2) + 1 : )) & - ] - end associate - end associate - end procedure - -#endif - module procedure to_file_t type(string_t), allocatable :: lines(:) integer, parameter :: inner_rows = 1 diff --git a/src/fortran/mimetic_operators_1D_m.F90 b/src/fortran/mimetic_operators_1D_m.F90 new file mode 100644 index 00000000..279be943 --- /dev/null +++ b/src/fortran/mimetic_operators_1D_m.F90 @@ -0,0 +1,163 @@ +#include "mole-language-support.F90" + +module mimetic_operators_1D_m + !! Define 1D scalar and vector abstractions and associated mimetic gradient and divergence operators. + use julienne_m, only : file_t + implicit none + + private + 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 + private + double precision, allocatable :: upper_(:,:) !! A submatrix block (cf. Corbino & Castillo, 2020) + double precision, allocatable :: inner_(:) !! M submatrix row (cf. Corbino & Castillo, 2020) + double precision, allocatable :: lower_(:,:) !! A' submatrix block (cf. Corbino & Castillo, 2020) + contains + procedure, non_overridable :: to_file_t + end type + + interface mimetic_matrix_1D_t + + pure module function construct_matrix_operator(upper, inner, lower) result(mimetic_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 + end function + + end interface + + type, extends(mimetic_matrix_1D_t) :: gradient_operator_1D_t + !! Encapsulate a 1D mimetic gradient operator + private + integer k_ !! order of accuracy + integer m_ !! number of cells + double precision dx_ !! cell width + contains + generic :: operator(.x.) => gradient_matrix_multiply + procedure, non_overridable, private :: gradient_matrix_multiply + generic :: assemble => assemble_gradient + procedure, non_overridable, private :: assemble_gradient + end type + + interface gradient_operator_1D_t + + pure module function construct_1D_gradient_operator(k, dx, cells) result(gradient_operator_1D) + !! Construct a mimetic gradient operator + implicit none + integer, intent(in) :: k !! order of accuracy + double precision, intent(in) :: dx !! step size + integer, intent(in) :: cells !! number of grid cells + type(gradient_operator_1D_t) gradient_operator_1D + end function + + end interface + + type, extends(mimetic_matrix_1D_t) :: divergence_operator_1D_t + !! Encapsulate kth-order mimetic divergence operator on m_ cells of width dx + private + integer k_, m_ + double precision dx_ + contains + generic :: operator(.x.) => divergence_matrix_multiply + procedure, non_overridable, private :: divergence_matrix_multiply + generic :: assemble => assemble_divergence + procedure, non_overridable, private :: assemble_divergence + procedure, non_overridable :: submatrix_A_rows + end type + + interface divergence_operator_1D_t + + pure module function construct_1D_divergence_operator(k, dx, cells) result(divergence_operator_1D) + !! Construct a mimetic gradient operator + implicit none + integer, intent(in) :: k !! order of accuracy + double precision, intent(in) :: dx !! step size + integer, intent(in) :: cells !! number of grid cells + type(divergence_operator_1D_t) divergence_operator_1D + end function + + end interface + + interface + + pure module function submatrix_A_rows(self) result(rows) + !! Result is number of rows in the A block of the mimetic divergence matrix operator + implicit none + class(divergence_operator_1D_t), intent(in) :: self + integer rows + end function + + pure module function gradient_matrix_multiply(self, vec) result(matvec_product) + !! Result is mimetic gradient vector + implicit none + class(gradient_operator_1D_t), intent(in) :: self + double precision, intent(in) :: vec(:) + double precision, allocatable :: matvec_product(:) + end function + + pure module function assemble_gradient(self) result(G) + !! Result is the assembled 1D mimetic gradient operator matrix + implicit none + class(gradient_operator_1D_t), intent(in) :: self + double precision, allocatable :: G(:,:) + end function + + pure module function assemble_divergence(self) result(D) + !! Result is the assembled 1D mimetic divergence operator matrix + implicit none + class(divergence_operator_1D_t), intent(in) :: self + double precision, allocatable :: D(:,:) + end function + + pure module function divergence_matrix_multiply(self, vec) result(matvec_product) + !! Result is mimetic divergence defined at cell centers + implicit none + class(divergence_operator_1D_t), intent(in) :: self + double precision, intent(in) :: vec(:) + double precision, allocatable :: matvec_product(:) + end function + + pure module function to_file_t(self) result(file) + implicit none + class(mimetic_matrix_1D_t), intent(in) :: self + type(file_t) file + end function + + end interface + +contains + +#if HAVE_DO_CONCURRENT_TYPE_SPEC_SUPPORT && HAVE_LOCALITY_SPECIFIER_SUPPORT + + pure function negate_and_flip(A) result(Ap) + !! Transform a mimetic matrix upper block into a lower block + double precision, intent(in) :: A(:,:) + double precision, allocatable :: Ap(:,:) + + allocate(Ap, mold=A) + + reverse_elements_within_rows_and_flip_sign: & + do concurrent(integer :: row = 1:size(Ap,1)) default(none) shared(Ap, A) + Ap(row,:) = -A(row,size(A,2):1:-1) + end do reverse_elements_within_rows_and_flip_sign + + reverse_elements_within_columns: & + do concurrent(integer :: column = 1 : size(Ap,2)) default(none) shared(Ap) + Ap(:,column) = Ap(size(Ap,1):1:-1,column) + end do reverse_elements_within_columns + + end function + +#else + +! see divergence_operator_1D_s and gradient_operator_1D_s + +#endif + +end module mimetic_operators_1D_m diff --git a/src/fortran/mole_m.f90 b/src/fortran/mole_m.f90 index 06fbda92..1eda112a 100644 --- a/src/fortran/mole_m.f90 +++ b/src/fortran/mole_m.f90 @@ -1,6 +1,14 @@ module mole_m - !! MOLE Fortran public entities - use scalar_1D_m, only : scalar_1D_t, scalar_1D_initializer_i - use gradient_1D_m, only : gradient_1D_t + !! All public MOLE Fortran entities: + + use tensors_1D_m, only : & + scalar_1D_t & ! discrete 1D scalar abstraction supporting mimetic gradient (.grad.) and Laplacian (.laplacian.) operators + ,vector_1D_t & ! discrete 1D vector abstraction supporting a mimetic divergence (.div.). operator + ,divergence_1D_t & ! result of applying the unary .div. operator to a vector_1D_t operand + ,laplacian_1D_t & ! result of applying the unary .laplacian. operator to a scalar_1D_t operand + ,scalar_1D_initializer_i & ! abstract interface for a scalar_1D_t initialization function + ,vector_1D_initializer_i ! abstract interface for a vector_1D_t initialization function + implicit none + end module mole_m diff --git a/src/fortran/scalar_1D_m.f90 b/src/fortran/scalar_1D_m.f90 deleted file mode 100644 index dc1e5c25..00000000 --- a/src/fortran/scalar_1D_m.f90 +++ /dev/null @@ -1,132 +0,0 @@ -module scalar_1D_m - !! Define an abstraction for the collection of points used to compute gradidents: - !! cell centers plus oundaries. - use julienne_m, only : file_t - use gradient_1D_m, only : gradient_1D_t - implicit none - - private - public :: scalar_1D_t - public :: gradient_operator_1D_t - public :: scalar_1D_initializer_i - - abstract interface - - pure function scalar_1D_initializer_i(x) result(f) - implicit none - double precision, intent(in) :: x(:) - double precision, allocatable :: f(:) - end function - - end interface - - type mimetic_matrix_1D_t - !! Encapsulate a mimetic matrix with a corresponding matrix-vector product operator - private - double precision, allocatable :: upper_(:,:), inner_(:), lower_(:,:) - contains - procedure to_file_t - end type - - type gradient_operator_1D_t - !! Encapsulate kth-order mimetic gradient operator on dx-sized cells - private - integer k_, m_ - double precision dx_ - type(mimetic_matrix_1D_t) mimetic_matrix_1D_ - end type - - type scalar_1D_t - !! Encapsulate information at cell centers and boundaries - private - double precision, allocatable :: scalar_1D_(:) - double precision x_min_, x_max_ - integer cells_ - type(gradient_operator_1D_t) gradient_operator_1D_ - contains - procedure grid - generic :: operator(.grad.) => grad - procedure, non_overridable, private :: grad - end type - - interface - - pure module function to_file_t(self) result(file) - implicit none - class(mimetic_matrix_1D_t), intent(in) :: self - type(file_t) file - end function - - end interface - - interface scalar_1D_t - - pure module function construct_from_function(initializer, order, cells, x_min, x_max) result(scalar_1D) - !! Result is a collection of cell-centered-extended values with a corresponding mimetic gradient operator - implicit none - procedure(scalar_1D_initializer_i), pointer :: initializer - integer, intent(in) :: order !! order of accuracy - integer, intent(in) :: cells !! number of grid cells spanning the domain - double precision, intent(in) :: x_min !! grid location minimum - double precision, intent(in) :: x_max !! grid location maximum - type(scalar_1D_t) scalar_1D - end function - - end interface - - interface - - pure module function grid(self) result(x) - !! Result is array of cell-centers-extended grid locations (cell centers + boundaries) - !! as described in Corbino & Castillo (2020) https://doi.org/10.1016/j.cam.2019.06.042 - implicit none - class(scalar_1D_t), intent(in) :: self - double precision, allocatable :: x(:) - end function - - pure module function grad(self) result(grad_f) - !! Result is mimetic gradient of f - implicit none - class(scalar_1D_t), intent(in) :: self - type(gradient_1D_t) grad_f !! discrete gradient approximation - end function - - end interface - - interface gradient_operator_1D_t - - pure module function construct_from_parameters(k, dx, m) result(gradient_operator_1D) - !! Construct a mimetic gradient operator - implicit none - integer, intent(in) :: k !! order of accuracy - double precision, intent(in) :: dx !! step siz - integer, intent(in) :: m !! number of grid cells - type(gradient_operator_1D_t) gradient_operator_1D - end function - - end interface - - interface mimetic_matrix_1D_t - - pure module function construct_from_components(upper, inner, lower) result(mimetic_matrix_1D) - !! Construct discrete operator from coefficient matrix - implicit none - double precision, intent(in) :: upper(:,:), inner(:), lower(:,:) - type(mimetic_matrix_1D_t) mimetic_matrix_1D - end function - - end interface - - interface - - pure module function matvec(self, scalar_1D) result(matvec_product) - !! Apply a mimetic matrix operator to a vector encapsulated in a scalar_1D_t object - implicit none - class(mimetic_matrix_1D_t), intent(in) :: self - type(scalar_1D_t), intent(in) :: scalar_1D - double precision, allocatable :: matvec_product(:) - end function - - end interface - -end module scalar_1D_m \ No newline at end of file diff --git a/src/fortran/scalar_1D_s.F90 b/src/fortran/scalar_1D_s.F90 index c7dde20f..9eae4be5 100644 --- a/src/fortran/scalar_1D_s.F90 +++ b/src/fortran/scalar_1D_s.F90 @@ -1,13 +1,27 @@ #include "julienne-assert-macros.h" -submodule(scalar_1D_m) scalar_1D_s - use julienne_m, only : call_julienne_assert_, operator(.equalsExpected.), operator(.greaterThan.), operator(.isAtLeast.) +submodule(tensors_1D_m) scalar_1D_s + use julienne_m, only : call_julienne_assert_, operator(.greaterThan.), operator(.isAtLeast.) implicit none contains - pure module function construct_from_function(initializer, order, cells, x_min, x_max) result(scalar_1D) - implicit none + +#ifndef __GFORTRAN__ + + module procedure construct_1D_scalar_from_function + call_julienne_assert(x_max .greaterThan. x_min) + call_julienne_assert(cells .isAtLeast. 2*order) + + associate(values => initializer(scalar_1D_grid_locations(x_min, x_max, cells))) + scalar_1D%tensor_1D_t = tensor_1D_t(values, x_min, x_max, cells, order) + end associate + scalar_1D%gradient_operator_1D_ = gradient_operator_1D_t(k=order, dx=(x_max - x_min)/cells, cells=cells) + end procedure + +#else + + pure module function construct_1D_scalar_from_function(initializer, order, cells, x_min, x_max) result(scalar_1D) procedure(scalar_1D_initializer_i), pointer :: initializer integer, intent(in) :: order !! order of accuracy integer, intent(in) :: cells !! number of grid cells spanning the domain @@ -18,30 +32,64 @@ pure module function construct_from_function(initializer, order, cells, x_min, x call_julienne_assert(x_max .greaterThan. x_min) call_julienne_assert(cells .isAtLeast. 2*order) - scalar_1D%x_min_ = x_min - scalar_1D%x_max_ = x_max - scalar_1D%cells_ = cells - scalar_1D%gradient_operator_1D_ = gradient_operator_1D_t(k=order, dx=(x_max - x_min)/cells, m=cells) - scalar_1D%scalar_1D_ = initializer(grid_(x_min, x_max, cells)) + associate(values => initializer(scalar_1D_grid_locations(x_min, x_max, cells))) + scalar_1D%tensor_1D_t = tensor_1D_t(values, x_min, x_max, cells, order) + end associate + scalar_1D%gradient_operator_1D_ = gradient_operator_1D_t(k=order, dx=(x_max - x_min)/cells, cells=cells) end function - pure function grid_(x_min, x_max, cells) result(x) + pure function cell_center_locations(x_min, x_max, cells) result(x) double precision, intent(in) :: x_min, x_max integer, intent(in) :: cells - double precision, allocatable :: x(:) + double precision, allocatable:: x(:) integer cell associate(dx => (x_max - x_min)/cells) - x = [x_min, x_min + dx/2. + [((cell-1)*dx, cell = 1, cells)], x_max] + x = x_min + dx/2. + [((cell-1)*dx, cell = 1, cells)] end associate end function - module procedure grid - x = grid_(self%x_min_, self%x_max_, self%cells_) - end procedure +#endif + module procedure grad - grad_f = gradient_1D_t(matvec(self%gradient_operator_1D_%mimetic_matrix_1D_, self), self%x_min_, self%x_max_, self%cells_) + gradient_1D = vector_1D_t( & + tensor_1D_t(self%gradient_operator_1D_ .x. self%values_, self%x_min_, self%x_max_, self%cells_, self%order_) & + ,divergence_operator_1D_t(self%order_, (self%x_max_-self%x_min_)/self%cells_, self%cells_) & + ) + end procedure + + module procedure laplacian + +#ifndef __GFORTRAN__ + laplacian_1D%divergence_1D_t = .div. (.grad. self) +#else + laplacian_1D%divergence_1D_t = div(grad(self)) +#endif + + associate(divergence_operator_1D => divergence_operator_1D_t(self%order_, (self%x_max_ - self%x_min_)/self%cells_, self%cells_)) + laplacian_1D%boundary_depth_ = divergence_operator_1D%submatrix_A_rows() + 1 + end associate + + end procedure + + module procedure scalar_1D_values + cell_centers_extended_values = self%values_ + end procedure + + pure function scalar_1D_grid_locations(x_min, x_max, cells) result(x) + double precision, intent(in) :: x_min, x_max + integer, intent(in) :: cells + double precision, allocatable:: x(:) + integer cell + + associate(dx => (x_max - x_min)/cells) + x = [x_min, cell_center_locations(x_min, x_max, cells), x_max] + end associate + end function + + module procedure scalar_1D_grid + cell_centers_extended = scalar_1D_grid_locations(self%x_min_, self%x_max_, self%cells_) end procedure end submodule scalar_1D_s \ No newline at end of file diff --git a/src/fortran/tensor_1D_s.f90 b/src/fortran/tensor_1D_s.f90 new file mode 100644 index 00000000..15630a1a --- /dev/null +++ b/src/fortran/tensor_1D_s.f90 @@ -0,0 +1,13 @@ +submodule(tensors_1D_m) tensor_1D_s + implicit none +contains + + module procedure construct_1D_tensor_from_components + tensor_1D%values_ = values + tensor_1D%x_min_ = x_min + tensor_1D%x_max_ = x_max + tensor_1D%cells_ = cells + tensor_1D%order_ = order + end procedure + +end submodule tensor_1D_s diff --git a/src/fortran/tensors_1D_m.F90 b/src/fortran/tensors_1D_m.F90 new file mode 100644 index 00000000..d9e76219 --- /dev/null +++ b/src/fortran/tensors_1D_m.F90 @@ -0,0 +1,249 @@ +#include "mole-language-support.F90" + +module tensors_1D_m + !! Define public 1D scalar and vector abstractions and associated mimetic gradient, + !! divergence, and Laplacian operators. + use julienne_m, only : file_t + use mimetic_operators_1D_m, only : divergence_operator_1D_t, gradient_operator_1D_t + + implicit none + + private + + public :: scalar_1D_t + public :: vector_1D_t + public :: laplacian_1D_t + public :: divergence_1D_t + public :: scalar_1D_initializer_i + public :: vector_1D_initializer_i + + abstract interface + + pure function scalar_1D_initializer_i(x) result(f) + !! Sampling function for initializing a scalar_1D_t object + implicit none + double precision, intent(in) :: x(:) + double precision, allocatable :: f(:) + end function + + pure function vector_1D_initializer_i(x) result(v) + !! Sampling function for initializing a vector_1D_t object + implicit none + double precision, intent(in) :: x(:) + double precision, allocatable :: v(:) + end function + + end interface + + type tensor_1D_t + !! Encapsulate the components that are common to all 1D tensors. + !! Child types define the operations supported by each child, including + !! gradient (.grad.) for scalars and divergence (.div.) for vectors. + private + double precision x_min_ !! domain lower boundary + double precision x_max_ !! domain upper boundary + integer cells_ !! number of grid cells spanning the domain + integer order_ !! order of accuracy of mimetic discretization + double precision, allocatable :: values_(:) !! tensor components at spatial locations + end type + + interface tensor_1D_t + + pure module function construct_1D_tensor_from_components(values, x_min, x_max, cells, order) result(tensor_1D) + !! User-defined constructor: result is a 1D tensor defined by assigning the dummy arguments to corresponding components + implicit none + double precision, intent(in) :: values(:) !! tensor components at grid locations define by child + double precision, intent(in) :: x_min !! grid location minimum + double precision, intent(in) :: x_max !! grid location maximum + integer, intent(in) :: cells !! number of grid cells spanning the domain + integer, intent(in) :: order !! order of accuracy + type(tensor_1D_t) tensor_1D + end function + + end interface + + type, extends(tensor_1D_t) :: scalar_1D_t + !! Encapsulate scalar values at cell centers and boundaries + private + type(gradient_operator_1D_t) gradient_operator_1D_ + contains + generic :: operator(.grad.) => grad + generic :: operator(.laplacian.) => laplacian + generic :: grid => scalar_1D_grid + generic :: values => scalar_1D_values + procedure, non_overridable, private :: grad + procedure, non_overridable, private :: laplacian + procedure, non_overridable, private :: scalar_1D_values + procedure, non_overridable, private :: scalar_1D_grid + end type + + interface scalar_1D_t + + pure module function construct_1D_scalar_from_function(initializer, order, cells, x_min, x_max) result(scalar_1D) + !! Result is a collection of cell-centered-extended values with a corresponding mimetic gradient operator + implicit none + procedure(scalar_1D_initializer_i), pointer :: initializer + integer, intent(in) :: order !! order of accuracy + integer, intent(in) :: cells !! number of grid cells spanning the domain + double precision, intent(in) :: x_min !! grid location minimum + double precision, intent(in) :: x_max !! grid location maximum + type(scalar_1D_t) scalar_1D + end function + + end interface + + type, extends(tensor_1D_t) :: vector_1D_t + !! Encapsulate 1D vector values at cell faces (nodes in 1D) and corresponding operators + private + type(divergence_operator_1D_t) divergence_operator_1D_ + contains + generic :: operator(.div.) => div + generic :: grid => vector_1D_grid + generic :: values => vector_1D_values + procedure, non_overridable, private :: div + procedure, non_overridable, private :: vector_1D_grid + procedure, non_overridable, private :: vector_1D_values + end type + + interface vector_1D_t + + pure module function construct_1D_vector_from_function(initializer, order, cells, x_min, x_max) result(vector_1D) + !! Result is a 1D vector with values initialized by the provided procedure pointer sampled on the specified + !! number of evenly spaced cells covering [x_min, x_max] + implicit none + procedure(vector_1D_initializer_i), pointer :: initializer + integer, intent(in) :: order !! order of accuracy + integer, intent(in) :: cells !! number of grid cells spanning the domain + double precision, intent(in) :: x_min !! grid location minimum + double precision, intent(in) :: x_max !! grid location maximum + type(vector_1D_t) vector_1D + end function + + pure module function construct_from_components(tensor_1D, divergence_operator_1D) result(vector_1D) + !! Result is a 1D vector with the provided parent component tensor_1D and the provided divergence operatror + type(tensor_1D_t), intent(in) :: tensor_1D + type(divergence_operator_1D_t), intent(in) :: divergence_operator_1D + type(vector_1D_t) vector_1D + end function + + end interface + + type, extends(tensor_1D_t) :: divergence_1D_t + !! Encapsulate divergences at cell centers + contains + generic :: grid => divergence_1D_grid + generic :: values => divergence_1D_values + procedure, non_overridable, private :: divergence_1D_values + procedure, non_overridable, private :: divergence_1D_grid + end type + + interface divergence_1D_t + + pure module function construct_from_tensor(tensor_1D) result(divergence_1D) + !! Result is a 1D divergence with the provided parent component + implicit none + type(tensor_1D_t), intent(in) :: tensor_1D + type(divergence_1D_t) divergence_1D + end function + + end interface + + type, extends(divergence_1D_t) :: laplacian_1D_t + private + integer boundary_depth_ + contains + procedure reduced_order_boundary_depth + end type + + interface + + pure module function scalar_1D_grid(self) result(cell_centers_extended) + !! Result is the array of locations at which 1D scalars are defined: cell centers agumented by spatial boundaries + implicit none + class(scalar_1D_t), intent(in) :: self + double precision, allocatable :: cell_centers_extended(:) + end function + + pure module function vector_1D_grid(self) result(cell_faces) + !! Result is the array of cell face locations (nodes in 1D) at which 1D vectors are defined + implicit none + class(vector_1D_t), intent(in) :: self + double precision, allocatable :: cell_faces(:) + end function + + pure module function divergence_1D_grid(self) result(cell_centers) + !! Result is the array of cell centers at which 1D divergences are defined + implicit none + class(divergence_1D_t), intent(in) :: self + double precision, allocatable :: cell_centers(:) + end function + + pure module function scalar_1D_values(self) result(cell_centers_extended_values) + !! Result is an array of 1D scalar values at boundaries and cell centers + implicit none + class(scalar_1D_t), intent(in) :: self + double precision, allocatable :: cell_centers_extended_values(:) + end function + + pure module function vector_1D_values(self) result(face_centered_values) + !! Result is an array of the 1D vector values at cell faces (nodes in 1D) + implicit none + class(vector_1D_t), intent(in) :: self + double precision, allocatable :: face_centered_values(:) + end function + + pure module function divergence_1D_values(self) result(cell_centered_values) + !! Result is an array of 1D divergences at cell centers + implicit none + class(divergence_1D_t), intent(in) :: self + double precision, allocatable :: cell_centered_values(:) + end function + + pure module function grad(self) result(gradient_1D) + !! Result is mimetic gradient of the scalar_1D_t "self" + implicit none + class(scalar_1D_t), intent(in) :: self + type(vector_1D_t) gradient_1D !! discrete gradient + end function + + pure module function laplacian(self) result(laplacian_1D) + !! Result is mimetic Laplacian of the scalar_1D_t "self" + implicit none + class(scalar_1D_t), intent(in) :: self + type(laplacian_1D_t) laplacian_1D !! discrete gradient + end function + + pure module function reduced_order_boundary_depth(self) result(num_nodes) + !! Result is number of nodes away from the boundary for which convergence rate is one degree lower + implicit none + class(laplacian_1D_t), intent(in) :: self + integer num_nodes + end function + + pure module function div(self) result(divergence_1D) + !! Result is mimetic divergence of the vector_1D_t "self" + implicit none + class(vector_1D_t), intent(in) :: self + type(divergence_1D_t) divergence_1D !! discrete divergence + end function + + end interface + +#ifndef __GFORTRAN__ + +contains + + pure function cell_center_locations(x_min, x_max, cells) result(x) + double precision, intent(in) :: x_min, x_max + integer, intent(in) :: cells + double precision, allocatable:: x(:) + integer cell + + associate(dx => (x_max - x_min)/cells) + x = x_min + dx/2. + [((cell-1)*dx, cell = 1, cells)] + end associate + end function + +#endif + +end module tensors_1D_m diff --git a/src/fortran/vector_1D_s.F90 b/src/fortran/vector_1D_s.F90 new file mode 100644 index 00000000..3099248f --- /dev/null +++ b/src/fortran/vector_1D_s.F90 @@ -0,0 +1,73 @@ +#include "julienne-assert-macros.h" + +submodule(tensors_1D_m) vector_1D_s + use julienne_m, only : call_julienne_assert_, operator(.greaterThan.), operator(.isAtLeast.), operator(.equalsExpected.) + implicit none + +contains + +#ifndef __GFORTRAN__ + + module procedure construct_1D_vector_from_function + call_julienne_assert(x_max .greaterThan. x_min) + call_julienne_assert(cells .isAtLeast. 2*order+1) + + associate(values => initializer(faces(x_min, x_max, cells))) + vector_1D%tensor_1D_t = tensor_1D_t(values, x_min, x_max, cells, order) + end associate + vector_1D%divergence_operator_1D_ = divergence_operator_1D_t(k=order, dx=(x_max - x_min)/cells, cells=cells) + end procedure + +#else + + pure module function construct_1D_vector_from_function(initializer, order, cells, x_min, x_max) result(vector_1D) + procedure(vector_1D_initializer_i), pointer :: initializer + integer, intent(in) :: order !! order of accuracy + integer, intent(in) :: cells !! number of grid cells spanning the domain + double precision, intent(in) :: x_min !! grid location minimum + double precision, intent(in) :: x_max !! grid location maximum + type(vector_1D_t) vector_1D + + call_julienne_assert(x_max .greaterThan. x_min) + call_julienne_assert(cells .isAtLeast. 2*order+1) + + associate(values => initializer(faces(x_min, x_max, cells))) + vector_1D%tensor_1D_t = tensor_1D_t(values, x_min, x_max, cells, order) + end associate + vector_1D%divergence_operator_1D_ = divergence_operator_1D_t(k=order, dx=(x_max - x_min)/cells, cells=cells) + end function + +#endif + + module procedure construct_from_components + vector_1D%tensor_1D_t = tensor_1D + vector_1D%divergence_operator_1D_ = divergence_operator_1D + end procedure + + module procedure div + associate(Dv => self%divergence_operator_1D_ .x. self%values_) + call_julienne_assert(size(Dv) .equalsExpected. self%cells_ + 2) + divergence_1D = divergence_1D_t( tensor_1D_t(Dv(2:size(Dv)-1), self%x_min_, self%x_max_, self%cells_, self%order_) ) + end associate + end procedure + + module procedure vector_1D_values + face_centered_values = self%values_ + end procedure + + pure function faces(x_min, x_max, cells) result(x) + double precision, intent(in) :: x_min, x_max + integer, intent(in) :: cells + double precision, allocatable:: x(:) + integer cell + + associate(dx => (x_max - x_min)/cells) + x = [x_min, x_min + [(cell*dx, cell = 1, cells-1)], x_max] + end associate + end function + + module procedure vector_1D_grid + cell_faces = faces(self%x_min_, self%x_max_, self%cells_) + end procedure + +end submodule vector_1D_s \ No newline at end of file diff --git a/test/divergence_operator_1D_test_m.F90 b/test/divergence_operator_1D_test_m.F90 new file mode 100644 index 00000000..8ed321bb --- /dev/null +++ b/test/divergence_operator_1D_test_m.F90 @@ -0,0 +1,197 @@ +#include "language-support.F90" + !! include Julienne preprocessor macros + +module divergence_operator_1D_test_m + use julienne_m, only : & + operator(//) & + ,operator(.all.) & + ,operator(.also.) & + ,operator(.approximates.) & + ,operator(.csv.) & + ,operator(.within.) & + ,string_t & + ,test_t & + ,test_description_t & + ,test_diagnosis_t & + ,test_result_t & + ,usher + use mole_m, only : vector_1D_t, vector_1D_initializer_i, scalar_1D_t, scalar_1D_initializer_i +#ifdef __GFORTRAN__ + use mole_m, only : divergence_1D_t +#endif + implicit none + + type, extends(test_t) :: divergence_operator_1D_test_t + contains + procedure, nopass :: subject + procedure, nopass :: results + end type + + double precision, parameter :: tight_tolerance = 1D-14, loose_tolerance = 1D-08, rough_tolerance = 1D-02, crude_tolerance = 2D-02 + +contains + + pure function subject() result(test_subject) + character(len=:), allocatable :: test_subject + test_subject = 'A 1D mimetic divergence operator' + end function + + function results() result(test_results) + type(divergence_operator_1D_test_t) divergence_operator_1D_test + type(test_result_t), allocatable :: test_results(:) + + test_results = divergence_operator_1D_test%run([ & + test_description_t( & + 'computing 2nd-order .div.(.grad. (x**2)/2) within ' // string_t(tight_tolerance) & + ,usher(check_2nd_order_div_grad_parabola)) & + ,test_description_t( & + 'computing 4th-order .div.(.grad. (x**2)/2) within ' // string_t(tight_tolerance) & + ,usher(check_4th_order_div_grad_parabola)) & + ,test_description_t( & + 'computing convergence rate of 2 for 2nd-order .div. [sin(x) + cos(x)] within ' // string_t(rough_tolerance) & + ,usher(check_2nd_order_div_sinusoid_convergence)) & + ,test_description_t( & + 'computing convergence rate of 4 for 4th-order .div. [sin(x) + cos(x)] within ' // string_t(crude_tolerance) & + ,usher(check_4th_order_div_sinusoid_convergence)) & + ]) + end function + + pure function parabola(x) result(y) + double precision, intent(in) :: x(:) + double precision, allocatable :: y(:) + y = (x**2)/2 + end function + + function check_2nd_order_div_grad_parabola() result(test_diagnosis) + type(test_diagnosis_t) test_diagnosis + procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => parabola + double precision, parameter :: expected_divergence = 1D0 +#ifdef __GFORTRAN__ + type(divergence_1D_t) div_grad_scalar + div_grad_scalar = .div. (.grad. scalar_1D_t(scalar_1D_initializer, order=2, cells=5, x_min=0D0, x_max=5D0)) +#else + associate(div_grad_scalar => .div. (.grad. scalar_1D_t(scalar_1D_initializer, order=2, cells=5, x_min=0D0, x_max=5D0))) +#endif + test_diagnosis = .all. (div_grad_scalar%values() .approximates. expected_divergence .within. tight_tolerance) & + // " (2nd-order .div. (.grad. (x**2)/2))" +#ifndef __GFORTRAN__ + end associate +#endif + end function + + function check_4th_order_div_grad_parabola() result(test_diagnosis) + type(test_diagnosis_t) test_diagnosis + procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => parabola + double precision, parameter :: expected_divergence = 1D0 +#ifdef __GFORTRAN__ + type(divergence_1D_t) div_grad_scalar + div_grad_scalar = .div. (.grad. scalar_1D_t(scalar_1D_initializer, order=4, cells=9, x_min=0D0, x_max=9D0)) +#else + associate(div_grad_scalar => .div. (.grad. scalar_1D_t(scalar_1D_initializer, order=4, cells=9, x_min=0D0, x_max=9D0))) +#endif + test_diagnosis = .all. (div_grad_scalar%values() .approximates. expected_divergence .within. tight_tolerance) & + // " (4th-order .div. (.grad. (x**2)/2))" +#ifndef __GFORTRAN__ + end associate +#endif + end function + + + pure function sinusoid(x) result(y) + double precision, intent(in) :: x(:) + double precision, allocatable :: y(:) + y = sin(x) + cos(x) + end function + + + function check_2nd_order_div_sinusoid_convergence() result(test_diagnosis) + type(test_diagnosis_t) test_diagnosis + procedure(vector_1D_initializer_i), pointer :: vector_1D_initializer => sinusoid + double precision, parameter :: pi = 3.141592653589793D0 + integer, parameter :: order_desired = 2, coarse_cells=100, fine_cells=200 +#ifdef __GFORTRAN__ + type(divergence_1D_t) div_coarse, div_fine + div_coarse = .div. vector_1D_t(vector_1D_initializer , order=order_desired, cells=coarse_cells, x_min=0D0, x_max=2*pi) + div_fine = .div. vector_1D_t(vector_1D_initializer , order=order_desired, cells=fine_cells , x_min=0D0, x_max=2*pi) +#else + associate( & + div_coarse => .div. vector_1D_t(vector_1D_initializer , order=order_desired, cells=coarse_cells, x_min=0D0, x_max=2*pi) & + ,div_fine => .div. vector_1D_t(vector_1D_initializer , order=order_desired, cells=fine_cells , x_min=0D0, x_max=2*pi) & + ) +#endif + associate( & + x_coarse => div_coarse%grid() & + ,x_fine => div_fine%grid()) + associate( & + grad_coarse => cos(x_coarse) - sin(x_coarse) & + ,grad_fine => cos(x_fine) - sin(x_fine) & + ,div_coarse_values => div_coarse%values() & + ,div_fine_values => div_fine%values() & + ) + test_diagnosis = .all. (div_coarse_values .approximates. grad_coarse .within. rough_tolerance) & + // " (coarse-grid 2nd-order .div. [sin(x) + cos(x)])" + test_diagnosis = test_diagnosis .also. (.all. (div_fine_values .approximates. grad_fine .within. rough_tolerance)) & + // " (fine-grid 2nd-order .div. [sin(x) + cos(x)])" + associate( & + error_coarse_max => maxval(abs(div_coarse_values - grad_coarse)) & + ,error_fine_max => maxval(abs(div_fine_values - grad_fine)) & + ) + associate(order_actual => log(error_coarse_max/error_fine_max)/log(dble(fine_cells)/coarse_cells)) + test_diagnosis = test_diagnosis .also. (order_actual .approximates. dble(order_desired) .within. rough_tolerance) & + // " (convergence rate for 2nd-order .div. [sin(x) + cos(x)])" + end associate + end associate + end associate + end associate +#ifndef __GFORTRAN__ + end associate +#endif + end function + + + function check_4th_order_div_sinusoid_convergence() result(test_diagnosis) + type(test_diagnosis_t) test_diagnosis + procedure(vector_1D_initializer_i), pointer :: vector_1D_initializer => sinusoid + double precision, parameter :: pi = 3.141592653589793D0 + integer, parameter :: order_desired = 4, coarse_cells=300, fine_cells=1500 +#ifdef __GFORTRAN__ + type(divergence_1D_t) div_coarse, div_fine + div_coarse = .div. vector_1D_t(vector_1D_initializer , order=order_desired, cells=coarse_cells, x_min=0D0, x_max=2*pi) + div_fine = .div. vector_1D_t(vector_1D_initializer , order=order_desired, cells=fine_cells , x_min=0D0, x_max=2*pi) +#else + associate( & + div_coarse => .div. vector_1D_t(vector_1D_initializer , order=order_desired, cells=coarse_cells, x_min=0D0, x_max=2*pi) & + ,div_fine => .div. vector_1D_t(vector_1D_initializer , order=order_desired, cells=fine_cells , x_min=0D0, x_max=2*pi) & + ) +#endif + associate( & + x_coarse => div_coarse%grid() & + ,x_fine => div_fine%grid() & + ) + associate( & + div_coarse_expected => cos(x_coarse) - sin(x_coarse) & + ,div_fine_expected => cos(x_fine) - sin(x_fine) & + ,div_coarse_values => div_coarse%values() & + ,div_fine_values => div_fine%values() & + ) + test_diagnosis = .all. (div_coarse_values .approximates. div_coarse_expected .within. loose_tolerance) & + // " (coarse-grid 4th-order .div. [sin(x) + cos(x)])" + test_diagnosis = test_diagnosis .also. (.all. (div_fine_values .approximates. div_fine_expected .within. loose_tolerance)) & + // " (fine-grid 4th-order .div. [sin(x) + cos(x)])" + associate( & + error_coarse_max => maxval(abs(div_coarse_values - div_coarse_expected)) & + ,error_fine_max => maxval(abs(div_fine_values - div_fine_expected)) & + ) + associate(order_actual => log(error_coarse_max/error_fine_max)/log(dble(fine_cells)/coarse_cells)) + test_diagnosis = test_diagnosis .also. (order_actual .approximates. dble(order_desired) .within. crude_tolerance) & + // " (convergence rate for 4th-order .div. [sin(x) + cos(x)])" + end associate + end associate + end associate + end associate +#ifndef __GFORTRAN__ + end associate +#endif + end function + +end module divergence_operator_1D_test_m \ No newline at end of file diff --git a/test/driver.f90 b/test/driver.f90 index 724b7fca..408507bd 100644 --- a/test/driver.f90 +++ b/test/driver.f90 @@ -1,10 +1,14 @@ program test_suite_driver 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 implicit none associate(test_harness => test_harness_t([ & test_fixture_t(gradient_operator_1D_test_t()) & + ,test_fixture_t(divergence_operator_1D_test_t()) & + ,test_fixture_t(laplacian_operator_1D_test_t()) & ])) call test_harness%report_results end associate diff --git a/test/gradient_operator_1D_test_m.F90 b/test/gradient_operator_1D_test_m.F90 index 1ef85f7b..fbab341f 100644 --- a/test/gradient_operator_1D_test_m.F90 +++ b/test/gradient_operator_1D_test_m.F90 @@ -3,22 +3,21 @@ module gradient_operator_1D_test_m use julienne_m, only : & - string_t & + operator(//) & + ,operator(.all.) & + ,operator(.also.) & + ,operator(.approximates.) & + ,operator(.within.) & + ,string_t & ,test_t & ,test_description_t & ,test_diagnosis_t & ,test_result_t & - ,operator(//) & - ,operator(.all.) & - ,operator(.also.) & - ,operator(.approximates.) & - ,operator(.csv.) & - ,operator(.within.) - use mole_m, only : scalar_1D_t, gradient_1D_t, scalar_1D_initializer_i -#if ! HAVE_PROCEDURE_ACTUAL_FOR_POINTER_DUMMY - use julienne_m, only : diagnosis_function_i + ,usher + use mole_m, only : scalar_1D_t, scalar_1D_initializer_i +#ifdef __GFORTRAN__ + use mole_m, only : vector_1D_t, vector_1D_initializer_i #endif - implicit none type, extends(test_t) :: gradient_operator_1D_test_t contains @@ -26,7 +25,7 @@ module gradient_operator_1D_test_m procedure, nopass :: results end type - double precision, parameter :: tight_tolerance = 1D-14, loose_tolerance = 1D-12, rough_tolerance = 1D-02, crude_tolerance = 5D-02 + double precision, parameter :: tight_tolerance = 1D-14, loose_tolerance = 1D-12, rough_tolerance = 2D-02 contains @@ -35,39 +34,24 @@ pure function subject() result(test_subject) test_subject = 'A 1D mimetic gradient operator' end function -#if HAVE_PROCEDURE_ACTUAL_FOR_POINTER_DUMMY - function results() result(test_results) type(gradient_operator_1D_test_t) gradient_operator_1D_test type(test_result_t), allocatable :: test_results(:) - test_results = gradient_operator_1D_test%run([ & - test_description_t('computing 2nd- & 4th-order 1D gradients of a constant within tolerance ' // string_t(tight_tolerance), check_grad_const) & - ,test_description_t('computing 2nd- & 4th-order 1D gradients of a line within tolerance ' // string_t(loose_tolerance), check_grad_line) & - ,test_description_t('computing 2nd- & 4th-order 1D gradients of a parabola within tolerance ' // string_t(loose_tolerance), check_grad_parabola) & - ,test_description_t('computing 2nd-order 1D gradients of a sinusoid with a convergence rate of 2 within tolerance ' // string_t(crude_tolerance), check_2nd_order_grad_sinusoid) & - ,test_description_t('computing 4th-order 1D gradients of a sinusoid with a convergence rate of 4 within tolerance ' // string_t(crude_tolerance), check_4th_order_grad_sinusoid) & - ]) - end function - -#else - function results() result(test_results) - type(gradient_operator_1D_test_t) gradient_operator_1D_test - type(test_result_t), allocatable :: test_results(:) - procedure(diagnosis_function_i), pointer :: & - check_grad_const_ptr => check_grad_const & - ,check_grad_line_ptr => check_grad_line & - ,check_grad_parabola_ptr => check_grad_parabola - - test_results = gradient_operator_1D_test%run([ & - test_description_t('computing 2nd & 4th-order 1D gradients of a constant within tolerance ' // string_t(tight_tolerance), check_grad_const_ptr) & - ,test_description_t('computing 2nd & 4th-order 1D gradients of a line within tolerance ' // string_t(loose_tolerance), check_grad_line_ptr) & - ,test_description_t('computing 2nd & 4th-order 1D gradients of a parabola within tolerance ' // string_t(loose_tolerance), check_grad_parabola_ptr) & + test_results = gradient_operator_1D_test%run([ & + test_description_t('computing 2nd- & 4th-order .grad. (5) within ' & + // string_t(tight_tolerance), usher(check_grad_const)) & + ,test_description_t('computing 2nd- & 4th-order .grad. (14*x + 3) within ' & + // string_t(loose_tolerance), usher(check_grad_line)) & + ,test_description_t('computing 2nd- & 4th-order .grad. (7*x**2 + 3*x + 5) within ' & + // string_t(loose_tolerance), usher(check_grad_parabola)) & + ,test_description_t('computing convergence rate of 2 for 2nd-order .grad. [sin(x) + cos(x)] within ' & + // string_t(rough_tolerance), usher(check_2nd_order_grad_convergence)) & + ,test_description_t('computing convergence rate of 4 for 4th-order .grad. [sin(x) + cos(x)] within ' & + // string_t(rough_tolerance), usher(check_4th_order_grad_convergence)) & ]) end function -#endif - pure function const(x) result(y) double precision, intent(in) :: x(:) double precision, allocatable :: y(:) @@ -75,17 +59,34 @@ pure function const(x) result(y) y = [(5D0, i=1,size(x))] end function + function check_grad_const() result(test_diagnosis) type(test_diagnosis_t) test_diagnosis - type(gradient_1D_t) grad_f - double precision, parameter :: df_dx = 0. + double precision, parameter :: grad_expected = 0. procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => const +#ifdef __GFORTRAN__ + type(vector_1D_t) grad - grad_f = .grad. scalar_1D_t(scalar_1D_initializer, order=2, cells=4, x_min=0D0, x_max=4D0) - test_diagnosis = .all. (grad_f%values() .approximates. df_dx .within. loose_tolerance) // " (2nd-order d(line)/dx)" + grad = .grad. scalar_1D_t(scalar_1D_initializer, order=2, cells=5, x_min=0D0, x_max=4D0) +#else + associate(grad => .grad. scalar_1D_t(scalar_1D_initializer, order=2, cells=5, x_min=0D0, x_max=4D0)) +#endif + test_diagnosis = .all. (grad%values() .approximates. grad_expected .within. loose_tolerance) & + // " (2nd-order .grad.(5))" +#ifndef __GFORTRAN__ + end associate +#endif - grad_f = .grad. scalar_1D_t(scalar_1D_initializer, order=4, cells=8, x_min=0D0, x_max=8D0) - test_diagnosis = test_diagnosis .also. (.all. (grad_f%values() .approximates. df_dx .within. loose_tolerance)) // " (4th-order d(line)/dx)" +#ifdef __GFORTRAN__ + grad = .grad. scalar_1D_t(scalar_1D_initializer, order=4, cells=9, x_min=0D0, x_max=8D0) +#else + associate(grad => .grad. scalar_1D_t(scalar_1D_initializer, order=4, cells=9, x_min=0D0, x_max=8D0)) +#endif + test_diagnosis = test_diagnosis .also. (.all. (grad%values() .approximates. grad_expected .within. loose_tolerance)) & + // " (4th-order .grad.(5))" +#ifndef __GFORTRAN__ + end associate +#endif end function pure function line(x) result(y) @@ -96,15 +97,32 @@ pure function line(x) result(y) function check_grad_line() result(test_diagnosis) type(test_diagnosis_t) test_diagnosis - type(gradient_1D_t) grad_f - double precision, parameter :: df_dx = 14D0 + double precision, parameter :: grad_expected = 14D0 procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => line +#ifdef __GFORTRAN__ + type(vector_1D_t) grad + + grad = .grad. scalar_1D_t(scalar_1D_initializer, order=2, cells=5, x_min=0D0, x_max=4D0) +#else + associate(grad => .grad. scalar_1D_t(scalar_1D_initializer, order=2, cells=5, x_min=0D0, x_max=4D0)) +#endif - grad_f = .grad. scalar_1D_t(scalar_1D_initializer, order=2, cells=4, x_min=0D0, x_max=4D0) - test_diagnosis = .all. (grad_f%values() .approximates. df_dx .within. loose_tolerance) // " (2nd-order d(line)/dx)" + test_diagnosis = .all. (grad%values() .approximates. grad_expected .within. loose_tolerance) & + // " (2nd-order .grad.(14*x + 3))" +#ifndef __GFORTRAN__ + end associate +#endif - grad_f = .grad. scalar_1D_t(scalar_1D_initializer, order=4, cells=8, x_min=0D0, x_max=8D0) - test_diagnosis = test_diagnosis .also. (.all. (grad_f%values() .approximates. df_dx .within. loose_tolerance)) // " (4th-order d(line)/dx)" +#ifdef __GFORTRAN__ + grad = .grad. scalar_1D_t(scalar_1D_initializer, order=4, cells=9, x_min=0D0, x_max=8D0) +#else + associate(grad => .grad. scalar_1D_t(scalar_1D_initializer, order=4, cells=9, x_min=0D0, x_max=8D0)) +#endif + test_diagnosis = test_diagnosis .also. (.all. (grad%values() .approximates. grad_expected .within. loose_tolerance)) & + // " (4th-order .grad.(14*x + 3))" +#ifndef __GFORTRAN__ + end associate +#endif end function @@ -114,29 +132,41 @@ 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 - type(scalar_1D_t) quadratic - type(gradient_1D_t) grad_f procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => parabola +#ifdef __GFORTRAN__ + type(vector_1D_t) grad - quadratic = scalar_1D_t(scalar_1D_initializer , order=2, cells=4, x_min=0D0, x_max=4D0) - grad_f = .grad. quadratic - - associate(x => grad_f%faces()) - associate(df_dx => 14*x + 3) - test_diagnosis = .all. (grad_f%values() .approximates. df_dx .within. loose_tolerance) // " (2nd-order d(parabola)/dx)" + grad = .grad. scalar_1D_t(scalar_1D_initializer , order=2, cells=5, x_min=0D0, x_max=4D0) +#else + associate(grad => .grad. scalar_1D_t(scalar_1D_initializer , order=2, cells=5, x_min=0D0, x_max=4D0)) +#endif + associate(x => grad%grid()) + associate(grad_expected => 14*x + 3) + test_diagnosis = .all. (grad%values() .approximates. grad_expected .within. loose_tolerance) & + // " (2nd-order .grad.(7*x**2 + 3*x + 5))" + end associate end associate +#ifndef __GFORTRAN__ end associate +#endif - quadratic = scalar_1D_t(scalar_1D_initializer , order=4, cells=8, x_min=0D0, x_max=8D0) - grad_f = .grad. quadratic - - associate(x => grad_f%faces()) - associate(df_dx => 14*x + 3) - test_diagnosis = test_diagnosis .also. (.all. (grad_f%values() .approximates. df_dx .within. loose_tolerance)) // " (4th-order d(parabola)/dx)" +#ifdef __GFORTRAN__ + grad = .grad. scalar_1D_t(scalar_1D_initializer , order=4, cells=9, x_min=0D0, x_max=8D0) +#else + associate(grad => .grad. scalar_1D_t(scalar_1D_initializer , order=4, cells=9, x_min=0D0, x_max=8D0)) +#endif + associate(x => grad%grid()) + associate(grad_expected => 14*x + 3) + test_diagnosis = test_diagnosis .also. (.all. (grad%values() .approximates. grad_expected .within. loose_tolerance)) & + // " (4th-order .grad.(7*x**2 + 3*x + 5))" + end associate end associate +#ifndef __GFORTRAN__ end associate +#endif end function pure function sinusoid(x) result(y) @@ -145,58 +175,97 @@ pure function sinusoid(x) result(y) y = sin(x) + cos(x) end function - function check_2nd_order_grad_sinusoid() result(test_diagnosis) + + function check_2nd_order_grad_convergence() result(test_diagnosis) type(test_diagnosis_t) test_diagnosis - type(scalar_1D_t) coarse, fine - type(gradient_1D_t) grad_coarse, grad_fine procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => sinusoid double precision, parameter :: pi = 3.141592653589793D0 - integer, parameter :: order_desired = 2, coarse_cells=100, fine_cells=1000 + integer, parameter :: order_desired = 2, coarse_cells=500, fine_cells=1500 +#ifdef __GFORTRAN__ + type(vector_1D_t) grad_coarse, grad_fine - coarse = scalar_1D_t(scalar_1D_initializer , order=order_desired, cells=coarse_cells, x_min=0D0, x_max=2*pi) - fine = scalar_1D_t(scalar_1D_initializer , order=order_desired, cells=fine_cells , x_min=0D0, x_max=2*pi) - - grad_coarse = .grad. coarse - grad_fine = .grad. fine - - associate(x_coarse => grad_coarse%faces(), x_fine => grad_fine%faces()) - associate(df_dx_coarse => cos(x_coarse) - sin(x_coarse), df_dx_fine => cos(x_fine) - sin(x_fine), grad_coarse_values => grad_coarse%values(), grad_fine_values => grad_fine%values()) - test_diagnosis = .all. (grad_coarse_values .approximates. df_dx_coarse .within. rough_tolerance) // " (2nd-order d(sinusoid)/dx point-wise errors)" - test_diagnosis = test_diagnosis .also. (.all. (grad_fine_values .approximates. df_dx_fine .within. rough_tolerance)) // " (2nd-order d(sinusoid)/dx point-wise)" - associate(error_coarse_max => maxval(abs(grad_coarse_values - df_dx_coarse)), error_fine_max => maxval(abs(grad_fine_values - df_dx_fine))) - associate(order_actual => log(error_coarse_max/error_fine_max)/log(dble(fine_cells)/coarse_cells)) - test_diagnosis = test_diagnosis .also. (order_actual .approximates. dble(order_desired) .within. crude_tolerance) // " (2nd-order d(sinusoid)/dx order of accuracy)" + grad_coarse = .grad. scalar_1D_t(scalar_1D_initializer , order=order_desired, cells=coarse_cells, x_min=0D0, x_max=2*pi) + grad_fine = .grad. scalar_1D_t(scalar_1D_initializer , order=order_desired, cells=fine_cells , x_min=0D0, x_max=2*pi) +#else + associate( & + grad_coarse => .grad. scalar_1D_t(scalar_1D_initializer , order=order_desired, cells=coarse_cells, x_min=0D0, x_max=2*pi) & + ,grad_fine => .grad. scalar_1D_t(scalar_1D_initializer , order=order_desired, cells=fine_cells , x_min=0D0, x_max=2*pi) & + ) +#endif + associate( & + x_coarse => grad_coarse%grid() & + ,x_fine => grad_fine%grid() & + ) + associate( & + grad_coarse_expected => cos(x_coarse) - sin(x_coarse) & + ,grad_fine_expected => cos(x_fine) - sin(x_fine) & + ,grad_coarse_values => grad_coarse%values() & + ,grad_fine_values => grad_fine%values() & + ) + test_diagnosis = .all. (grad_coarse_values .approximates. grad_coarse_expected .within. rough_tolerance) & + // " (coarse-grid 2nd-order .grad. [sin(x) + cos(x)])" + test_diagnosis = test_diagnosis .also. (.all. (grad_fine_values .approximates. grad_fine_expected .within. rough_tolerance)) & + // " (fine-grid 4th-order .grad. [sin(x) + cos(x)])" + associate( & + error_coarse_max => maxval(abs(grad_coarse_values - grad_coarse_expected)) & + ,error_fine_max => maxval(abs(grad_fine_values - grad_fine_expected)) & + ) + associate(order_actual => log(error_coarse_max/error_fine_max)/log(dble(fine_cells)/coarse_cells)) + test_diagnosis = test_diagnosis .also. (order_actual .approximates. dble(order_desired) .within. rough_tolerance) & + // " (2nd-order .grad. [sin(x) + cos(x)] order of accuracy)" + end associate end associate end associate end associate +#ifndef __GFORTRAN__ end associate +#endif end function - function check_4th_order_grad_sinusoid() result(test_diagnosis) + function check_4th_order_grad_convergence() result(test_diagnosis) type(test_diagnosis_t) test_diagnosis - type(scalar_1D_t) coarse, fine - type(gradient_1D_t) grad_coarse, grad_fine procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => sinusoid double precision, parameter :: pi = 3.141592653589793D0 - integer, parameter :: order_desired = 4, coarse_cells=100, fine_cells=1000 - - coarse = scalar_1D_t(scalar_1D_initializer , order=order_desired, cells=coarse_cells, x_min=0D0, x_max=2*pi) - fine = scalar_1D_t(scalar_1D_initializer , order=order_desired, cells=fine_cells , x_min=0D0, x_max=2*pi) - - grad_coarse = .grad. coarse - grad_fine = .grad. fine + integer, parameter :: order_desired = 4, coarse_cells=100, fine_cells=1600 +#ifdef __GFORTRAN__ + type(vector_1D_t) grad_coarse, grad_fine - associate(x_coarse => grad_coarse%faces(), x_fine => grad_fine%faces()) - associate(df_dx_coarse => cos(x_coarse) - sin(x_coarse), df_dx_fine => cos(x_fine) - sin(x_fine), grad_coarse_values => grad_coarse%values(), grad_fine_values => grad_fine%values()) - test_diagnosis = .all. (grad_coarse_values .approximates. df_dx_coarse .within. rough_tolerance) // " (4th-order d(sinusoid)/dx point-wise errors)" - test_diagnosis = test_diagnosis .also. (.all. (grad_fine_values .approximates. df_dx_fine .within. rough_tolerance)) // " (4th-order d(sinusoid)/dx point-wise)" - associate(error_coarse_max => maxval(abs(grad_coarse_values - df_dx_coarse)), error_fine_max => maxval(abs(grad_fine_values - df_dx_fine))) - associate(order_actual => log(error_coarse_max/error_fine_max)/log(dble(fine_cells)/coarse_cells)) - test_diagnosis = test_diagnosis .also. (order_actual .approximates. dble(order_desired) .within. crude_tolerance) // " (4th-order d(sinusoid)/dx order of accuracy)" + grad_coarse = .grad. scalar_1D_t(scalar_1D_initializer , order=order_desired, cells=coarse_cells, x_min=0D0, x_max=2*pi) + grad_fine = .grad. scalar_1D_t(scalar_1D_initializer , order=order_desired, cells=fine_cells , x_min=0D0, x_max=2*pi) +#else + associate( & + grad_coarse => .grad. scalar_1D_t(scalar_1D_initializer , order=order_desired, cells=coarse_cells, x_min=0D0, x_max=2*pi) & + ,grad_fine => .grad. scalar_1D_t(scalar_1D_initializer , order=order_desired, cells=fine_cells , x_min=0D0, x_max=2*pi) & + ) +#endif + associate( & + x_coarse => grad_coarse%grid() & + ,x_fine => grad_fine%grid() & + ) + associate( & + grad_coarse_expected => cos(x_coarse) - sin(x_coarse) & + ,grad_fine_expected => cos(x_fine) - sin(x_fine) & + ,grad_coarse_values => grad_coarse%values() & + ,grad_fine_values => grad_fine%values() & + ) + test_diagnosis = .all. (grad_coarse_values .approximates. grad_coarse_expected .within. rough_tolerance) & + // " (4th-order d(sinusoid)/dx point-wise errors)" + test_diagnosis = test_diagnosis .also. (.all. (grad_fine_values .approximates. grad_fine_expected .within. rough_tolerance)) & + // " (4th-order d(sinusoid)/dx point-wise)" + associate( & + error_coarse_max => maxval(abs(grad_coarse_values - grad_coarse_expected)) & + ,error_fine_max => maxval(abs(grad_fine_values - grad_fine_expected)) & + ) + associate(order_actual => log(error_coarse_max/error_fine_max)/log(dble(fine_cells)/coarse_cells)) + test_diagnosis = test_diagnosis .also. (order_actual .approximates. dble(order_desired) .within. rough_tolerance) & + // " (4th-order d(sinusoid)/dx order of accuracy)" + end associate end associate end associate end associate +#ifndef __GFORTRAN__ end associate +#endif end function end module \ No newline at end of file diff --git a/test/laplacian_operator_1D_test_m.F90 b/test/laplacian_operator_1D_test_m.F90 new file mode 100644 index 00000000..a6e96af7 --- /dev/null +++ b/test/laplacian_operator_1D_test_m.F90 @@ -0,0 +1,203 @@ +module laplacian_operator_1D_test_m + use julienne_m, only : & + file_t & + ,operator(//) & + ,operator(.all.) & + ,operator(.also.) & + ,operator(.approximates.) & + ,operator(.separatedBy.) & + ,operator(.within.) & + ,string_t & + ,test_t & + ,test_description_t & + ,test_diagnosis_t & + ,test_result_t & + ,usher + use mole_m, only : scalar_1D_t, scalar_1D_initializer_i +#ifdef __GFORTRAN__ + use mole_m, only : laplacian_1D_t +#endif + implicit none + + type, extends(test_t) :: laplacian_operator_1D_test_t + contains + procedure, nopass :: subject + procedure, nopass :: results + end type + + double precision, parameter :: tight_tolerance = 1D-14, loose_tolerance = 1D-09, crude_tolerance = 1D-02 + +contains + + pure function subject() result(test_subject) + character(len=:), allocatable :: test_subject + test_subject = 'A 1D mimetic laplacian operator' + end function + + function results() result(test_results) + type(laplacian_operator_1D_test_t) laplacian_operator_1D_test + type(test_result_t), allocatable :: test_results(:) + + test_results = laplacian_operator_1D_test%run([ & + test_description_t( & + 'computing 2nd-order .laplacian. [(x**2)/2] within ' // string_t(tight_tolerance) & + ,usher(check_2nd_order_laplacian_parabola)) & + ,test_description_t( & + 'computing 4th-order .laplacian. [(x**4)/12] within ' // string_t(loose_tolerance) & + ,usher(check_4th_order_laplacian_of_quartic)) & + ,test_description_t( & + 'converging as dx^2 internally and dx near boundary for 2nd-order .laplacian. sin(x) within ' // string_t(crude_tolerance) & + ,usher(check_2nd_order_laplacian_convergence)) & + ,test_description_t( & + 'converging as dx^4 internally and dx^3 near boundary for 4th-order .laplacian. sin(x) within ' // string_t(crude_tolerance) & + ,usher(check_4th_order_laplacian_convergence)) & + ]) + end function + + pure function parabola(x) result(y) + double precision, intent(in) :: x(:) + double precision, allocatable :: y(:) + y = (x**2)/2 + end function + + function check_2nd_order_laplacian_parabola() result(test_diagnosis) + type(test_diagnosis_t) test_diagnosis + procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => parabola + double precision, parameter :: expected_laplacian = 1D0 +#ifdef __GFORTRAN__ + type(laplacian_1D_t) laplacian_scalar + laplacian_scalar = .laplacian. scalar_1D_t(scalar_1D_initializer, order=2, cells=5, x_min=0D0, x_max=5D0) +#else + associate(laplacian_scalar => .laplacian. scalar_1D_t(scalar_1D_initializer, order=2, cells=5, x_min=0D0, x_max=5D0)) +#endif + test_diagnosis = .all. (laplacian_scalar%values() .approximates. expected_laplacian .within. tight_tolerance) & + // " (2nd-order .laplacian. [(x**2)/2]" +#ifndef __GFORTRAN__ + end associate +#endif + end function + + pure function quartic(x) result(y) + double precision, intent(in) :: x(:) + double precision, allocatable :: y(:) + y = (x**4)/12 + end function + + function check_4th_order_laplacian_of_quartic() result(test_diagnosis) + type(test_diagnosis_t) test_diagnosis + procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => quartic + +#ifndef __GFORTRAN__ + associate(laplacian_quartic => .laplacian. scalar_1D_t(scalar_1D_initializer, order=4, cells=20, x_min=0D0, x_max=40D0)) +#else + type(laplacian_1D_t) laplacian_quartic + laplacian_quartic = .laplacian. scalar_1D_t(scalar_1D_initializer, order=4, cells=20, x_min=0D0, x_max=40D0) +#endif + associate(x => laplacian_quartic%grid()) + associate(expected_laplacian => x**2, actual_laplacian => laplacian_quartic%values()) + test_diagnosis = .all. (actual_laplacian .approximates. expected_laplacian .within. loose_tolerance) & + // " (4th-order .laplacian. [(x**4)/24]" + end associate + end associate +#ifndef __GFORTRAN__ + end associate +#endif + end function + + pure function f(x) + double precision, intent(in) :: x(:) + double precision, allocatable :: f(:) + f = sin(x) + end function + + pure function d2f_dx2(x) + double precision, intent(in) :: x(:) + double precision, allocatable :: d2f_dx2(:) + d2f_dx2 = -sin(x) + end function + + function check_2nd_order_laplacian_convergence() result(test_diagnosis) + type(test_diagnosis_t) test_diagnosis + test_diagnosis = check_laplacian_convergence(order_desired=2, coarse_cells=500, fine_cells=1000) + end function + + function check_4th_order_laplacian_convergence() result(test_diagnosis) + type(test_diagnosis_t) test_diagnosis + test_diagnosis = check_laplacian_convergence(order_desired = 4, coarse_cells=100, fine_cells=200) + end function + + function check_laplacian_convergence(order_desired, coarse_cells, fine_cells) result(test_diagnosis) + type(test_diagnosis_t) test_diagnosis + procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => f + double precision, parameter :: pi = 3.141592653589793D0 + integer, intent(in) :: order_desired, coarse_cells, fine_cells + +#ifndef __GFORTRAN__ + associate( & + laplacian_coarse => .laplacian. scalar_1D_t(scalar_1D_initializer , order=order_desired, cells=coarse_cells, x_min=0D0, x_max=2*pi) & + ,laplacian_fine => .laplacian. scalar_1D_t(scalar_1D_initializer , order=order_desired, cells=fine_cells , x_min=0D0, x_max=2*pi) & + ) +#else + type(laplacian_1D_t) laplacian_coarse, laplacian_fine + laplacian_coarse = .laplacian. scalar_1D_t(scalar_1D_initializer , order=order_desired, cells=coarse_cells, x_min=0D0, x_max=2*pi) + laplacian_fine = .laplacian. scalar_1D_t(scalar_1D_initializer , order=order_desired, cells=fine_cells , x_min=0D0, x_max=2*pi) +#endif + grids: & + associate( & + x_coarse => laplacian_coarse%grid() & + ,x_fine => laplacian_fine%grid()) + + laplacian_values: & + associate( & + expected_coarse => d2f_dx2(x_coarse) & + ,expected_fine => d2f_dx2(x_fine) & + ,actual_coarse => laplacian_coarse%values() & + ,actual_fine => laplacian_fine%values() & + ,depth => laplacian_coarse%reduced_order_boundary_depth() & + ) + test_diagnosis = & + .all. (actual_coarse .approximates. expected_coarse .within. crude_tolerance) & + // " (coarse-grid 2nd-order .laplacian. sin(x))" + + test_diagnosis = test_diagnosis .also. & + (.all. (actual_fine .approximates. expected_fine .within. crude_tolerance)) & + // " (fine-grid 2nd-order .laplacian. sin(x))" + + check_internal_convergence_rate: & + associate( & + coarse_error_max => maxval( abs( & + actual_coarse(1+depth:size(actual_coarse)-depth) - expected_coarse(1+depth:size(expected_coarse)-depth) & + )) & + ,fine_error_max => maxval( abs( & + actual_fine(1+depth:size(actual_fine)-depth) - expected_fine(1+depth:size(expected_fine)-depth) & + ) )) + associate(order_actual => log(coarse_error_max/fine_error_max)/log(dble(fine_cells)/coarse_cells)) + test_diagnosis = test_diagnosis .also. (order_actual .approximates. dble(order_desired) .within. crude_tolerance) & + // " (boundary convergence rate as dx^" // string_t(order_desired) // " for .laplacian. sin(x))" + end associate + end associate check_internal_convergence_rate + + check_boundary_convergence_rate: & + associate( & + coarse_error_max => maxval( abs( & + [ actual_coarse(1:depth-1), actual_coarse(size(actual_coarse)-depth+1:)] & + -[expected_coarse(1:depth-1), expected_coarse(size(actual_coarse)-depth+1:)] & + )) & + ,fine_error_max => maxval( abs( & + [ actual_fine(1:depth-1), actual_fine(size(actual_fine)-depth+1:)] & + -[expected_fine(1:depth-1), expected_fine(size(actual_fine)-depth+1:)] & + ) )) + associate(order_actual => log(coarse_error_max/fine_error_max)/log(dble(fine_cells)/coarse_cells)) + test_diagnosis = test_diagnosis .also. (order_actual .approximates. dble(order_desired-1) .within. crude_tolerance) & + // " (boundary convergence rate as dx^" // string_t(order_desired-1) // " for .laplacian. sin(x))" + end associate + end associate check_boundary_convergence_rate + + end associate laplacian_values + end associate grids +#ifndef __GFORTRAN__ + end associate +#endif + end function + +end module \ No newline at end of file