diff --git a/.gitmodules b/.gitmodules index 43f7e6b..038d4ca 100644 --- a/.gitmodules +++ b/.gitmodules @@ -6,3 +6,7 @@ path = src/third_party/PENF url = https://github.com/szaghi/PENF branch = master +[submodule "src/third_party/FLAP"] + path = src/third_party/FLAP + url = https://github.com/szaghi/FLAP + branch = master diff --git a/fobos b/fobos index 1d89474..14c35db 100644 --- a/fobos +++ b/fobos @@ -14,7 +14,7 @@ $DEBUG_GNU = -Og -g3 -Warray-bounds -Wcharacter-truncation -Wline-truncation - ; $DEBUG_GNU = -Og -g3 -Warray-bounds -Wcharacter-truncation -Wline-truncation -Wimplicit-interface -Wimplicit-procedure -Wunderflow -fcheck=all -fmodule-private -ffree-line-length-132 -fimplicit-none -fbacktrace -fdump-core -finit-real=nan -std=f2008 -fall-intrinsics $DEBUG_INT = -O0 -debug all -check all -warn all -extend-source 132 -traceback -gen-interfaces#-fpe-all=0 -fp-stack-check -fstack-protector-all -ftrapuv -no-ftz -std08 $OPTIMIZE = -O2 -$EXDIRS = PENF/src/tests/ pyplot-fortran/src/tests/ +$EXDIRS = FLAP/src/tests/ PENF/src/tests/ pyplot-fortran/src/tests/ # main modes # GNU diff --git a/src/lib/wenoof_interpolator_js.f90 b/src/lib/wenoof_interpolator_js.f90 index dc3adad..e6c99e2 100644 --- a/src/lib/wenoof_interpolator_js.f90 +++ b/src/lib/wenoof_interpolator_js.f90 @@ -138,7 +138,7 @@ pure subroutine interpolate(self, S, stencil, location, interpolation) enddo ! computing the convolution - interpolation = 0. + interpolation = 0._R_P do k = 0, S - 1 ! stencils loop do f = f1, f2 ! 1 => left interface (i-1/2), 2 => right interface (i+1/2) interpolation(f + ff) = interpolation(f + ff) + weights(f, k) * self%polynom%poly(f, k) diff --git a/src/tests/sin_reconstruction.f90 b/src/tests/sin_reconstruction.f90 index a7fa09e..23d60fb 100644 --- a/src/tests/sin_reconstruction.f90 +++ b/src/tests/sin_reconstruction.f90 @@ -1,65 +1,354 @@ !< WenOOF test: reconstruction of sin function. + +module test_module +!< Auxiliary module defining the test class. + +use flap, only : command_line_interface +use penf, only : I_P, R_P, FR_P, str, strz +use pyplot_module, only : pyplot +use wenoof, only : interpolator, wenoof_create + +implicit none +private +public :: test + +character(99), parameter :: interpolators(1:4) = ["all ", & + "JS ", & + "JS-Z", & + "JS-M"] !< List of available interpolators. +real(R_P), parameter :: pi = 4._R_P * atan(1._R_P) !< Pi greek. + +type :: solution_data + !< Class to handle solution data. + real(R_P), allocatable :: x_cell(:) !< Cell domain [1-S:points_number+S]. + real(R_P), allocatable :: fx_cell(:) !< Cell refecence values [1-S:points_number+S]. + real(R_P), allocatable :: x_face(:) !< Face domain [1:points_number]. + real(R_P), allocatable :: fx_face(:) !< Face reference values [1:points_number]. + real(R_P), allocatable :: interpolation(:) !< Interpolated values [1:points_number]. + real(R_P) :: error_L2 !< L2 norm of the numerical error. +endtype solution_data + +type :: test + !< Class to handle test(s). + !< + !< Test is driven by the Command Line Interface (CLI) options. + !< + !< Test has only 1 public method `execute`: it executes test(s) accordingly to cli options. + private + type(command_line_interface) :: cli !< Command line interface handler. + integer(I_P) :: error=0 !< Error handler. + character(99) :: interpolator_type='JS' !< Interpolator used. + character(99) :: output_bname='unset' !< Output files basename. + integer(I_P) :: pn_number !< Number of different points-number tested. + integer(I_P), allocatable :: points_number(:) !< Points number used to discretize the domain. + integer(I_P) :: S_number !< Number of different stencils tested. + integer(I_P), allocatable :: S(:) !< Stencils used. + type(solution_data), allocatable :: solution(:,:) !< Solution [1:pn_number, 1:S_number]. + real(R_P), allocatable :: accuracy(:,:) !< Accuracy (measured) [1:pn_number-1, 1:S_number]. + logical :: errors_analysis=.false. !< Flag for activating errors analysis. + logical :: plots=.false. !< Flag for activating plots saving. + logical :: results=.false. !< Flag for activating results saving. + contains + ! public methods + procedure, pass(self) :: execute !< Execute selected test(s). + ! private methods + procedure, pass(self), private :: allocate_solution_data !< Allocate solution data. + procedure, pass(self), private :: analize_errors !< Analize errors. + procedure, pass(self), private :: compute_reference_solution !< Compute reference solution. + procedure, pass(self), private :: deallocate_solution_data !< Deallocate solution data. + procedure, pass(self), private :: initialize !< Initialize test(s). + procedure, pass(self), private :: perform !< Perform test(s). + procedure, pass(self), private :: save_results_and_plots !< Save results and plots. +endtype test + +contains + ! public methods + subroutine execute(self) + !< Execute test(s). + class(test), intent(inout) :: self !< Test. + integer(I_P) :: s !< Counter. + + call self%initialize + if (trim(adjustl(self%interpolator_type))/='all') then + call self%perform + else + do s=2, size(interpolators, dim=1) + self%interpolator_type = trim(adjustl(interpolators(s))) + call self%perform + enddo + endif + endsubroutine execute + + ! private methods + subroutine allocate_solution_data(self) + !< Allocate solution data. + class(test), intent(inout) :: self !< Test. + integer(I_P) :: s !< Counter. + integer(I_P) :: pn !< Counter. + + call self%deallocate_solution_data + self%pn_number = size(self%points_number, dim=1) + self%S_number = size(self%S, dim=1) + allocate(self%solution(1:self%pn_number, 1:self%S_number)) + if (self%pn_number>1) then + allocate(self%accuracy(1:self%pn_number-1, 1:self%S_number)) + self%accuracy = 0._R_P + endif + do s=1, self%S_number + do pn=1, self%pn_number + allocate(self%solution(pn, s)%x_cell( 1-self%S(s):self%points_number(pn)+self%S(s))) + allocate(self%solution(pn, s)%fx_cell(1-self%S(s):self%points_number(pn)+self%S(s))) + allocate(self%solution(pn, s)%x_face( 1:self%points_number(pn) )) + allocate(self%solution(pn, s)%fx_face( 1:self%points_number(pn) )) + allocate(self%solution(pn, s)%interpolation( 1:self%points_number(pn) )) + self%solution(pn, s)%x_cell = 0._R_P + self%solution(pn, s)%fx_cell = 0._R_P + self%solution(pn, s)%x_face = 0._R_P + self%solution(pn, s)%fx_face = 0._R_P + self%solution(pn, s)%interpolation = 0._R_P + enddo + enddo + endsubroutine allocate_solution_data + + subroutine compute_reference_solution(self) + !< Allocate solution data. + class(test), intent(inout) :: self !< Test. + integer(I_P) :: s !< Counter. + integer(I_P) :: pn !< Counter. + integer(I_P) :: i !< Counter. + + call self%allocate_solution_data + do s=1, self%S_number + do pn=1, self%pn_number + ! compute the values used for the reconstruction of sin function: cell values + do i=1 - self%S(s), self%points_number(pn) + self%S(s) + self%solution(pn, s)%x_cell(i) = i * 2 * pi / self%points_number(pn) + self%solution(pn, s)%fx_cell(i) = sin(self%solution(pn, s)%x_cell(i)) + enddo + ! face values to which the reconstruction should tend + do i = 1, self%points_number(pn) + self%solution(pn, s)%x_face(i) = self%solution(pn, s)%x_cell(i) + pi / self%points_number(pn) + self%solution(pn, s)%fx_face(i) = sin(self%solution(pn, s)%x_face(i)) + enddo + enddo + enddo + endsubroutine compute_reference_solution + + subroutine deallocate_solution_data(self) + !< Deallocate solution data. + class(test), intent(inout) :: self !< Test. + + if (allocated(self%solution)) deallocate(self%solution) + if (allocated(self%accuracy)) deallocate(self%accuracy) + endsubroutine deallocate_solution_data + + subroutine initialize(self) + !< Initialize test: set Command Line Interface, parse it and check its validity. + class(test), intent(inout) :: self !< Test. + + call set_cli + call parse_cli + contains + subroutine set_cli() + !< Set Command Line Interface. + + associate(cli => self%cli) + call cli%init(progname = 'sin reconstruction', & + authors = 'Fortran-FOSS-Programmers', & + license = 'GNU GPLv3', & + description = 'Test WenOOF library on sin function reconstruction', & + examples = ["sin_reconstruction --interpolator JS --results", & + "sin_reconstruction --interpolator JS-Z -r ", & + "sin_reconstruction --interpolator JS-M ", & + "sin_reconstruction --interpolator all -p -r "]) + call cli%add(switch='--interpolator', switch_ab='-i', help='WENO interpolator type', required=.false., def='JS', act='store') + call cli%add(switch='--points_number', switch_ab='-pn', nargs='+', help='Number of points used to discretize the domain', & + required=.false., act='store', def='50') + call cli%add(switch='--stencils', switch_ab='-s', nargs='+', help='Stencils dimensions (and number)', & + required=.false., act='store', def='2', choices='2, 3, 4, 5, 6, 7, 8, 9') + call cli%add(switch='--results', switch_ab='-r', help='Save results', required=.false., act='store_true', def='.false.') + call cli%add(switch='--plots', switch_ab='-p', help='Save plots', required=.false., act='store_true', def='.false.') + call cli%add(switch='--output', help='Output files basename', required=.false., act='store', def='sin_reconstruction') + call cli%add(switch='--errors_analysis', help='Peform errors analysis', required=.false., act='store_true', def='.false.') + endassociate + endsubroutine set_cli + + subroutine parse_cli() + !< Parse Command Line Interface and check its validity. + character(len=:), allocatable :: valid_solvers_list !< Pretty printed list of available solvers. + + call self%cli%parse(error=self%error) ; if (self%error/=0) stop + call self%cli%get(switch='-i', val=self%interpolator_type, error=self%error) ; if (self%error/=0) stop + call self%cli%get_varying(switch='-pn', val=self%points_number, error=self%error) ; if (self%error/=0) stop + call self%cli%get_varying(switch='-s', val=self%S, error=self%error) ; if (self%error/=0) stop + call self%cli%get(switch='-r', val=self%results, error=self%error) ; if (self%error/=0) stop + call self%cli%get(switch='-p', val=self%plots, error=self%error) ; if (self%error/=0) stop + call self%cli%get(switch='--output', val=self%output_bname, error=self%error) ; if (self%error/=0) stop + call self%cli%get(switch='--errors_analysis', val=self%errors_analysis, error=self%error) ; if (self%error/=0) stop + + if (.not.is_interpolator_valid()) then + print "(A)", 'error: the interpolator type "'//trim(adjustl(self%interpolator_type))//'" is unknown!' + print "(A)", list_interpolators() + stop + endif + endsubroutine parse_cli + + function is_interpolator_valid() + !< Verify if the selected interpolator is valid. + logical :: is_interpolator_valid !< Return true is the selected interpolator is available. + integer(I_P) :: s !< Counter. + + is_interpolator_valid = .false. + do s=1, size(interpolators, dim=1) + is_interpolator_valid = (trim(adjustl(self%interpolator_type))==trim(adjustl(interpolators(s)))) + if (is_interpolator_valid) exit + enddo + endfunction is_interpolator_valid + + function list_interpolators() result(list) + !< List available solvers. + character(len=:), allocatable :: list !< Pretty printed list of available interpolators. + integer(I_P) :: s !< Counter. + + list = 'Valid interpolator names are:' // new_line('a') + do s=1, ubound(interpolators, dim=1) + list = list // ' + ' // trim(adjustl(interpolators(s))) // new_line('a') + enddo + endfunction list_interpolators + endsubroutine initialize + + subroutine perform(self) + !< Perform the test. + class(test), intent(inout) :: self !< Test. + real(R_P), allocatable :: error(:,:) !< Error (norm L2) with respect the exact solution. + real(R_P), allocatable :: order(:,:) !< Observed order based on subsequent refined solutions. + class(interpolator), allocatable :: weno_interpolator !< WENO interpolator. + integer(I_P) :: s !< Counter. + integer(I_P) :: pn !< Counter. + integer(I_P) :: i !< Counter. + + call self%compute_reference_solution + do s=1, self%S_number + call wenoof_create(interpolator_type=trim(adjustl(self%interpolator_type)), S=self%S(s), wenoof_interpolator=weno_interpolator) + do pn=1, self%pn_number + do i=1, self%points_number(pn) + call weno_interpolator%interpolate(S=self%S(s), & + stencil=reshape(source=self%solution(pn, s)%fx_cell(i+1-self%S(s):i-1+self%S(s)), & + shape=[1,2*self%S(s)-1]), & + location='right', & + interpolation=self%solution(pn, s)%interpolation(i:i)) + enddo + enddo + enddo + call self%analize_errors + call self%save_results_and_plots + endsubroutine perform + + subroutine save_results_and_plots(self) + !< Save results (and plots). + class(test), intent(inout) :: self !< Test. + type(pyplot) :: plt !< Plot handler. + character(len=:), allocatable :: title !< Plot title + integer(I_P) :: file_unit !< File unit. + integer(I_P) :: s !< Counter. + integer(I_P) :: pn !< Counter. + integer(I_P) :: i !< Counter. + + if (self%results) then + open(newunit=file_unit, file=trim(self%output_bname)//'.dat') + write(file_unit, "(A)") 'VARIABLES = "x" "sin(x)" "weno_interpolation"' + do s=1, self%S_number + do pn=1, self%pn_number + write(file_unit, "(A)") 'ZONE T = "'//'S_'//trim(str(self%S(s), .true.))//& + '-Np_'//trim(str(self%points_number(pn), .true.))//'"' + do i = 1, self%points_number(pn) + write(file_unit, "(A)") str(n=self%solution(pn, s)%x_face(i))//' '// & + str(n=self%solution(pn, s)%fx_face(i))//' '//& + str(n=self%solution(pn, s)%interpolation(i)) + enddo + enddo + enddo + close(file_unit) + + if (self%errors_analysis.and.self%pn_number>1) then + open(newunit=file_unit, file=trim(self%output_bname)//'-accuracy.dat') + write(file_unit, "(A)") 'VARIABLES = "S" "Np" "error (L2)" "observed order"' + do s=1, self%S_number + do pn=1, self%pn_number - 1 + write(file_unit, "(A)") trim(str(self%S(s), .true.))//' '//& + trim(str(self%points_number(pn+1), .true.))//' '//& + trim(str(self%solution(pn+1, s)%error_L2))//' '//& + trim(str(self%accuracy(pn, s))) + enddo + enddo + close(file_unit) + endif + endif + if (self%plots) then + do s=1, self%S_number + do pn=1, self%pn_number + title = 'WENO interpolation of $\sin(x)$; '//& + 'S='//trim(str(self%S(s), .true.))//'Np='//trim(str(self%points_number(pn), .true.)) + call plt%initialize(grid=.true., xlabel='angle (rad)', title=title, legend=.true.) + call plt%add_plot(x=self%solution(pn, s)%x_face(:), & + y=self%solution(pn, s)%fx_face(:), & + label='$\sin(x)$', & + linestyle='k-', & + linewidth=2, & + ylim=[-1.1_R_P, 1.1_R_P]) + call plt%add_plot(x=self%solution(pn, s)%x_face(:), & + y=self%solution(pn, s)%interpolation(:), & + label='WENO interpolation', & + linestyle='ro', & + markersize=6, & + ylim=[-1.1_R_P, 1.1_R_P]) + call plt%savefig('sin_reconstruction'//& + '-S_'//trim(str(self%S(s), .true.))//'-Np_'//trim(str(self%points_number(pn), .true.))//'.png') + enddo + enddo + endif + endsubroutine save_results_and_plots + + subroutine analize_errors(self) + !< Analize errors. + class(test), intent(inout) :: self !< Test. + integer(I_P) :: s !< Counter. + integer(I_P) :: pn !< Counter. + integer(I_P) :: i !< Counter. + + if (self%errors_analysis) then + do s=1, self%S_number + do pn=1, self%pn_number + associate(error_L2=>self%solution(pn, s)%error_L2, & + fx_face=>self%solution(pn, s)%fx_face, & + interpolation=>self%solution(pn, s)%interpolation) + error_L2 = 0._R_P + do i=1, self%points_number(pn) + error_L2 = error_L2 + (interpolation(i) - fx_face(i))**2 + enddo + error_L2 = sqrt(error_L2) + endassociate + enddo + enddo + if (self%pn_number>1) then + do s=1, self%S_number + do pn=1, self%pn_number-1 + self%accuracy(pn, s) = log(self%solution(pn, s)%error_L2 / self%solution(pn + 1, s)%error_L2) / & + log((1._R_P*self%points_number(pn)) / self%points_number(pn + 1)) + enddo + enddo + endif + endif + endsubroutine analize_errors +endmodule test_module + program sin_reconstruction !< WenOOF test: reconstruction of sin function. -use penf, only : I_P, R_P, str -use wenoof, only : interpolator, wenoof_create -use pyplot_module, only : pyplot +use test_module implicit none -class(interpolator), allocatable :: weno_interpolator !< WENO interpolator. -integer(I_P), parameter :: S = 6_I_P !< Stencils used. -integer(I_P), parameter :: Nv = 30_I_P !< Number of discretized values to be interpolated. -real(R_P), parameter :: pi = 4._R_P * atan(1._R_P) !< Extent of domain. -real(R_P) :: x(1-S:Nv+S) !< Whole domain. -real(R_P) :: fx(1-S:Nv+S) !< Discretized values to be interpolated. -real(R_P) :: xi(1:Nv) !< Domain of the interpolation. -real(R_P) :: fx_ref(1:Nv) !< Reference values. -real(R_P) :: interpolation(1:1, 1:Nv) !< Interpolated values. -type(pyplot) :: plt !< Plotter handler. -integer :: i, j, f !< Counters. - -! build the values used for the reconstruction of sin function: nodal values -x = 0. -do i = 1 - S, Nv + S - x(i) = i * 2 * pi / Nv - fx(i) = sin(x(i)) -enddo -! face values to which the reconstruction should tend -do i = 1, Nv - xi(i) = x(i) + pi / Nv - fx_ref(i) = sin(xi(i)) -enddo - -! create weno interpolator -call wenoof_create(interpolator_type='JS', S=S, wenoof_interpolator=weno_interpolator) - -! interpolate values -interpolation = 0. -do i = 1, Nv ! interpolated values loop - call weno_interpolator%interpolate(S=S, & - stencil=reshape(source=fx(i+1-S:i-1+S), shape=[1,2*S-1]), & - location='right', & - interpolation=interpolation(1:1, i)) -enddo - -! print results -print "(A)", '# x, sin(x), weno_interpolation(x)' -do i = 1, Nv - print "(A)", str(n=xi(i))//' '//str(n=fx_ref(i))//' '//str(n=interpolation(1, i)) -enddo - -! plotting graph to image file -call plt%initialize(grid=.true., xlabel='angle (rad)', title='WENO interpolation of $\sin(x)$', legend=.true.) -call plt%add_plot(x=xi(1:Nv), & - y=fx_ref(1:Nv), & - label='$\sin(x)$', & - linestyle='k-', & - linewidth=2) -call plt%add_plot(x=xi(1:Nv), & - y=interpolation(1, 1:Nv), & - label='WENO interpolation', & - linestyle='ro', & - markersize=6) -call plt%savefig('sin_reconstruction.png') +type(test) :: sin_test + +call sin_test%execute endprogram sin_reconstruction diff --git a/src/third_party/FLAP b/src/third_party/FLAP new file mode 160000 index 0000000..9e60188 --- /dev/null +++ b/src/third_party/FLAP @@ -0,0 +1 @@ +Subproject commit 9e6018852452fed4d1315fe83f21550d8fb4b0a9