restore error computation into tester

Author: szaghi

src/tests/tester/foodie_test_integrand_ladvection.f90

@@ -59,10 +59,10 @@ module foodie_test_integrand_ladvection
       ! auxiliary methods
       procedure, pass(self), public :: destroy          !< Destroy field.
       procedure, pass(self), public :: dt => compute_dt !< Compute the current time step, by means of CFL condition.
-      procedure, pass(self), public :: exact_solution   !< Return exact solution.
       procedure, pass(self), public :: output           !< Extract integrand state field.
       ! integrand_tester_object deferred methods
       procedure, pass(self), public :: description    !< Return an informative description of the test.
+      procedure, pass(self), public :: exact_solution !< Return exact solution.
       procedure, pass(self), public :: export_tecplot !< Export integrand to Tecplot file.
       procedure, pass(self), public :: parse_cli      !< Initialize from command line interface.
       procedure, nopass,     public :: set_cli        !< Set command line interface.
@@ -125,15 +125,6 @@ pure function compute_dt(self, final_time, t) result(Dt)
    endassociate
    endfunction compute_dt
 
-   pure function exact_solution(self, u0, t) result(exact)
-   !< Return exact solution.
-   class(integrand_ladvection), intent(in) :: self     !< Integrand.
-   real(R_P),                   intent(in) :: u0(1:)   !< Initial state
-   real(R_P),                   intent(in) :: t        !< Time.
-   real(R_P), allocatable                  :: exact(:) !< Exact solution.
-
-   endfunction exact_solution
-
    pure function output(self) result(state)
    !< Output the advection field state.
    class(integrand_ladvection), intent(in) :: self     !< Advection field.
@@ -154,6 +145,15 @@ pure function description(self, prefix) result(desc)
    desc = prefix//'linear_advection-Ni_'//trim(strz(self%Ni, 10))
    endfunction description
 
+   pure function exact_solution(self, t, U0) result(exact)
+   !< Return exact solution.
+   class(integrand_ladvection), intent(in)           :: self     !< Integrand.
+   real(R_P),                   intent(in)           :: t        !< Time.
+   class(integrand_object),     intent(in), optional :: U0       !< Initial conditions.
+   real(R_P), allocatable                            :: exact(:) !< Exact solution.
+
+   endfunction exact_solution
+
    subroutine export_tecplot(self, file_name, t, scheme, close_file)
    !< Export integrand to Tecplot file.
    class(integrand_ladvection), intent(in)           :: self            !< Advection field.

src/tests/tester/foodie_test_integrand_oscillation.f90

@@ -50,11 +50,11 @@ module foodie_test_integrand_oscillation
    contains
       ! auxiliary methods
       procedure, pass(self), public :: amplitude_phase !< Return amplitude and phase of the oscillation.
-      procedure, pass(self), public :: exact_solution  !< Return exact solution.
       procedure, pass(self), public :: initialize      !< Initialize integrand.
       procedure, pass(self), public :: output          !< Extract integrand state field.
       ! integrand_tester_object deferred methods
       procedure, pass(self), public :: description    !< Return an informative description of the test.
+      procedure, pass(self), public :: exact_solution !< Return exact solution.
       procedure, pass(self), public :: export_tecplot !< Export integrand to Tecplot file.
       procedure, pass(self), public :: parse_cli      !< Initialize from command line interface.
       procedure, nopass,     public :: set_cli        !< Set command line interface.
@@ -99,16 +99,6 @@ function amplitude_phase(self) result(ap)
    ap(2) = atan(-self%U(1) / self%U(2))
    endfunction amplitude_phase
 
-   pure function exact_solution(self, t) result(exact)
-   !< Return exact solution.
-   class(integrand_oscillation), intent(in) :: self       !< Integrand.
-   real(R_P),                    intent(in) :: t          !< Time.
-   real(R_P)                                :: exact(1:2) !< Exact solution.
-
-   exact(1) = self%U0(1) * cos(self%f * t) - self%U0(2) * sin(self%f * t)
-   exact(2) = self%U0(1) * sin(self%f * t) + self%U0(2) * cos(self%f * t)
-   endfunction exact_solution
-
    pure subroutine initialize(self, U0, frequency)
    !< Initialize integrand.
    class(integrand_oscillation), intent(inout) :: self      !< Integrand.
@@ -140,6 +130,17 @@ pure function description(self, prefix) result(desc)
    desc = prefix//'oscillation'
    endfunction description
 
+   pure function exact_solution(self, t, U0) result(exact)
+   !< Return exact solution.
+   class(integrand_oscillation), intent(in)           :: self     !< Integrand.
+   real(R_P),                    intent(in)           :: t        !< Time.
+   class(integrand_object),      intent(in), optional :: U0       !< Initial conditions.
+   real(R_P), allocatable                             :: exact(:) !< Exact solution.
+
+   exact = [self%U0(1) * cos(self%f * t) - self%U0(2) * sin(self%f * t), &
+            self%U0(1) * sin(self%f * t) + self%U0(2) * cos(self%f * t)]
+   endfunction exact_solution
+
    subroutine export_tecplot(self, file_name, t, scheme, close_file)
    !< Export integrand to Tecplot file.
    class(integrand_oscillation), intent(in)           :: self            !< Advection field.

src/tests/tester/foodie_test_integrand_tester_object.f90

@@ -17,6 +17,7 @@ module foodie_test_integrand_tester_object
    !< This abstract provided some auxiliary methods useful for the tester machinery.
    contains
       procedure(description_interface),    pass(self), deferred :: description    !< Return an informative description of the test.
+      procedure(exact_solution_interface), pass(self), deferred :: exact_solution !< Return exact solution.
       procedure(export_tecplot_interface), pass(self), deferred :: export_tecplot !< Export integrand to Tecplot file.
       procedure(parse_cli_interface),      pass(self), deferred :: parse_cli      !< Initialize from command line interface.
       procedure(set_cli_interface),        nopass,     deferred :: set_cli        !< Set command line interface.
@@ -32,6 +33,15 @@ pure function description_interface(self, prefix) result(desc)
    character(len=:), allocatable                        :: desc   !< Description.
    endfunction description_interface
 
+   pure function exact_solution_interface(self, t, U0) result(exact)
+   !< Return exact solution.
+   import :: integrand_object, integrand_tester_object, R_P
+   class(integrand_tester_object), intent(in)           :: self     !< Integrand.
+   real(R_P),                      intent(in)           :: t        !< Time.
+   class(integrand_object),        intent(in), optional :: U0       !< Initial conditions.
+   real(R_P), allocatable                               :: exact(:) !< Exact solution.
+   endfunction exact_solution_interface
+
    subroutine export_tecplot_interface(self, file_name, t, scheme, close_file)
    !< Export integrand to Tecplot file.
    import :: integrand_tester_object, R_P

src/tests/tester/foodie_tester.f90

@@ -55,7 +55,6 @@ module foodie_test_object
    logical                                     :: save_results        !< Flag for activating results saving.
    character(99)                               :: output              !< Output files basename.
    integer(I_P)                                :: save_frequency      !< Save frequency.
-   logical                                     :: verbose             !< Flag for activating verbose output.
    type(integrand_ladvection)                  :: ladvection_0        !< Initial conditions for linear advection test.
    type(integrand_oscillation)                 :: oscillation_0       !< Initial conditions for oscillation test.
    class(integrand_tester_object), allocatable :: integrand_0         !< Initial conditions.
@@ -72,6 +71,8 @@ subroutine execute(self)
    !< Execute test(s).
    class(test_object), intent(inout) :: self                  !< Test.
    character(99), allocatable        :: integrator_schemes(:) !< Name of FOODIE integrator schemes.
+   real(R_P), allocatable            :: error(:,:)            !< Error of integrand integration.
+   real(R_P), allocatable            :: order(:,:)            !< Observed order of accuravy.
    integer(I_P)                      :: s                     !< Counter.
    integer(I_P)                      :: t                     !< Counter.
 
@@ -85,7 +86,10 @@ subroutine execute(self)
    else
       integrator_schemes = foodie_integrator_schemes()
    endif
+   allocate(error(1:self%integrand_0%integrand_dimension(), 1:size(self%Dt, dim=1)))
+   if (size(self%Dt, dim=1) > 1) allocate(order(1:self%integrand_0%integrand_dimension(), 1:size(error, dim=2)-1))
    do s=1, size(integrator_schemes, dim=1)
+      print '(A)', trim(integrator_schemes(s))
       do t=1, size(self%Dt)
          call integrate(scheme=trim(integrator_schemes(s)),          &
                         integrand_0=self%integrand_0,                &
@@ -96,7 +100,13 @@ subroutine execute(self)
                         is_fast=self%is_fast,                        &
                         save_results=self%save_results,              &
                         output_base_name=trim(adjustl(self%output)), &
-                        save_frequency=self%save_frequency)
+                        save_frequency=self%save_frequency,          &
+                        error=error(:,t))
+         print*, 'Dt = ', self%Dt(t), ', error = ', error(:,t)
+         if (t > 1) then
+            order(:, t-1) = observed_order(error=error(:, t-1:t), Dt=self%Dt(t-1:t))
+            print '(A,'//trim(str(size(order, dim=1), no_sign=.true.))//'F6.2)', ' Observed order =', order(:,t-1)
+         endif
       enddo
    enddo
    endsubroutine execute
@@ -130,7 +140,6 @@ subroutine set_cli()
          call cli%add(switch='--save_results', switch_ab='-r',help='save result', required=.false., act='store_true', def='.false.')
          call cli%add(switch='--output', help='output file basename', required=.false., act='store', def='foodie_test')
          call cli%add(switch='--save_frequency', help='save frequency', required=.false., act='store', def='1')
-         call cli%add(switch='--verbose', help='Verbose output', required=.false., act='store_true', def='.false.')
       endassociate
       call self%ladvection_0%set_cli(cli=self%cli)
       call self%oscillation_0%set_cli(cli=self%cli)
@@ -154,7 +163,6 @@ subroutine parse_cli()
       call self%cli%get(switch='-r', val=self%save_results, error=self%error) ; if (self%error/=0) stop
       call self%cli%get(switch='--output', val=self%output, error=self%error) ; if (self%error/=0) stop
       call self%cli%get(switch='--save_frequency', val=self%save_frequency, error=self%error) ; if (self%error/=0) stop
-      call self%cli%get(switch='--verbose', val=self%verbose, error=self%error) ; if (self%error/=0) stop
       call self%ladvection_0%parse_cli(cli=self%cli)
       call self%oscillation_0%parse_cli(cli=self%cli)
 
@@ -228,23 +236,24 @@ subroutine check_scheme_has_fast_mode(scheme, integrator)
    endsubroutine check_scheme_has_fast_mode
 
    subroutine integrate(scheme, integrand_0, Dt, final_time, iterations, stages, is_fast, save_results, output_base_name, &
-                        save_frequency)
+                        save_frequency, error)
    !< Integrate integrand by means of the given scheme.
-   character(*),                   intent(in)  :: scheme           !< Selected scheme.
-   class(integrand_tester_object), intent(in)  :: integrand_0      !< Initial conditions.
-   real(R_P),                      intent(in)  :: Dt               !< Time step.
-   real(R_P),                      intent(in)  :: final_time       !< Final integration time.
-   integer(I_P),                   intent(in)  :: iterations       !< Number of fixed point iterations.
-   integer(I_P),                   intent(in)  :: stages           !< Number of stages.
-   logical,                        intent(in)  :: is_fast          !< Activate fast mode integration.
-   logical,                        intent(in)  :: save_results     !< Save results.
-   character(*),                   intent(in)  :: output_base_name !< Base name of output results file.
-   integer(I_P),                   intent(in)  :: save_frequency   !< Save frequency.
-   class(integrator_object), allocatable       :: integrator       !< The integrator.
-   type(integrator_runge_kutta_ssp)            :: integrator_start !< The (auto) start integrator.
-   class(integrand_tester_object), allocatable :: integrand        !< Integrand.
-   real(R_P)                                   :: time             !< Time.
-   integer(I_P)                                :: step             !< Time steps counter.
+   character(*),                   intent(in)   :: scheme           !< Selected scheme.
+   class(integrand_tester_object), intent(in)   :: integrand_0      !< Initial conditions.
+   real(R_P),                      intent(in)   :: Dt               !< Time step.
+   real(R_P),                      intent(in)   :: final_time       !< Final integration time.
+   integer(I_P),                   intent(in)   :: iterations       !< Number of fixed point iterations.
+   integer(I_P),                   intent(in)   :: stages           !< Number of stages.
+   logical,                        intent(in)   :: is_fast          !< Activate fast mode integration.
+   logical,                        intent(in)   :: save_results     !< Save results.
+   character(*),                   intent(in)   :: output_base_name !< Base name of output results file.
+   integer(I_P),                   intent(in)   :: save_frequency   !< Save frequency.
+   real(R_P),                      intent(out)  :: error(:)         !< Error of integrand integration.
+   class(integrand_tester_object) , allocatable :: integrand        !< Integrand.
+   class(integrator_object), allocatable        :: integrator       !< The integrator.
+   type(integrator_runge_kutta_ssp)             :: integrator_start !< The (auto) start integrator.
+   real(R_P)                                    :: time             !< Time.
+   integer(I_P)                                 :: step             !< Time steps counter.
 
    allocate(integrand, mold=integrand_0) ; integrand = integrand_0
 
@@ -319,14 +328,32 @@ subroutine integrate(scheme, integrand_0, Dt, final_time, iterations, stages, is
       enddo
    endselect
 
-   select type(integrand)
-   type is(integrand_ladvection)
-      if (save_results .and. mod(step, save_frequency)==0) call integrand%export_tecplot(t=time, scheme=scheme)
-   type is(integrand_oscillation)
-      if (save_results .and. mod(step, save_frequency)==0) call integrand%export_tecplot(t=time)
-   endselect
-   if (save_results) call integrand%export_tecplot(close_file=.true.)
+   call integrand_close_tecplot
+
+   call integrand_compute_error
    contains
+      subroutine integrand_close_tecplot
+      !< Close current integrand tecplot file.
+
+      select type(integrand)
+      type is(integrand_ladvection)
+         if (save_results .and. mod(step, save_frequency)==0) call integrand%export_tecplot(t=time, scheme=scheme)
+      type is(integrand_oscillation)
+         if (save_results .and. mod(step, save_frequency)==0) call integrand%export_tecplot(t=time)
+      endselect
+      if (save_results) call integrand%export_tecplot(close_file=.true.)
+      endsubroutine integrand_close_tecplot
+
+      subroutine integrand_compute_error
+      !< Close current integrand tecplot file.
+
+      select type(integrand)
+      type is(integrand_ladvection)
+      type is(integrand_oscillation)
+         error = abs(integrand%output() - integrand%exact_solution(t=time))
+      endselect
+      endsubroutine integrand_compute_error
+
       subroutine integrand_export_tecplot
       !< Export current integrand solution to tecplot file.
 
@@ -338,6 +365,18 @@ subroutine integrand_export_tecplot
       endselect
       endsubroutine integrand_export_tecplot
    endsubroutine integrate
+
+   pure function observed_order(error, Dt)
+   !< Estimate the order of accuracy using 2 subsequent refined numerical solutions.
+   real(R_P), intent(in) :: error(1:, 1:)                        !< Computed errors.
+   real(R_P), intent(in) :: Dt(1:)                               !< Time steps used.
+   real(R_P)             :: observed_order(1:size(error, dim=1)) !< Estimation of the order of accuracy.
+   integer(I_P)          :: v                                    !< Variables counter.
+
+   do v=1, size(error, dim=1)
+     observed_order(v) = log(error(v, 1) / error(v, 2)) / log(Dt(1) / Dt(2))
+   enddo
+   endfunction observed_order
 endmodule foodie_test_object
 
 program foodie_tester