improve tester

Author: szaghi

src/tests/tester/foodie_test_integrand_ladvection.f90

@@ -55,15 +55,19 @@ module foodie_test_integrand_ladvection
    real(R_P)                               :: a=0._R_P        !< Advection coefficient.
    real(R_P), allocatable                  :: u(:)            !< Integrand (state) variable.
    class(interpolator_object), allocatable :: interpolator    !< WENO interpolator.
+   character(99)                           :: initial_state   !< Initial state.
    contains
       ! 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 :: dt => compute_dt !< Compute the current time step by means of CFL condition.
+      procedure, pass(self), public ::       compute_dx !< Compute the space step by means of CFL condition.
       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 :: error          !< Return error.
       procedure, pass(self), public :: exact_solution !< Return exact solution.
       procedure, pass(self), public :: export_tecplot !< Export integrand to Tecplot file.
+      procedure, pass(self), public :: initialize     !< Initialize integrand.
       procedure, pass(self), public :: parse_cli      !< Initialize from command line interface.
       procedure, nopass,     public :: set_cli        !< Set command line interface.
       ! integrand_object deferred methods
@@ -117,14 +121,25 @@ pure function compute_dt(self, final_time, t) result(Dt)
    real(R_P),                   intent(in), optional :: t          !< Time.
    real(R_P)                                         :: Dt         !< Time step.
 
-   associate(a=>self%a, Ni=>self%Ni, Dx=>self%Dx, CFL=>self%CFL)
+   associate(a=>self%a, Dx=>self%Dx, CFL=>self%CFL)
       Dt = Dx * CFL / abs(a)
       if (present(t)) then
          if ((t + Dt) > final_time) Dt = final_time - t
       endif
    endassociate
    endfunction compute_dt
 
+   pure function compute_dx(self, Dt) result(Dx)
+   !< Compute the space step step by means of CFL condition.
+   class(integrand_ladvection), intent(in) :: self       !< Advection field.
+   real(R_P),                   intent(in) :: Dt         !< Time step.
+   real(R_P)                               :: Dx         !< Space step.
+
+   associate(a=>self%a, CFL=>self%CFL)
+      Dx = Dt / CFL * abs(a)
+   endassociate
+   endfunction compute_dx
+
    pure function output(self) result(state)
    !< Output the advection field state.
    class(integrand_ladvection), intent(in) :: self     !< Advection field.
@@ -145,6 +160,28 @@ pure function description(self, prefix) result(desc)
    desc = prefix//'linear_advection-Ni_'//trim(strz(self%Ni, 10))
    endfunction description
 
+   pure function error(self, t, U0)
+   !< Return error.
+   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                            :: error(:) !< Error.
+   integer(I_P)                                      :: i        !< Counter.
+
+   allocate(error(1:1))
+   error = 0._R_P
+   if (present(U0)) then
+      select type(U0)
+      type is(integrand_ladvection)
+         do i=1, self%Ni
+            ! error = error + (U0%u(i) - self%u(i)) ** 2
+            error = max(error, abs(U0%u(i) - self%u(i)))
+         enddo
+      endselect
+      ! error = sqrt(error)
+   endif
+   endfunction error
+
    pure function exact_solution(self, t, U0) result(exact)
    !< Return exact solution.
    class(integrand_ladvection), intent(in)           :: self     !< Integrand.
@@ -194,11 +231,24 @@ subroutine export_tecplot(self, file_name, t, scheme, close_file)
    endif
    endsubroutine export_tecplot
 
+   subroutine initialize(self, Dt)
+   !< Initialize integrand.
+   class(integrand_ladvection), intent(inout) :: self !< Integrand.
+   real(R_P),                   intent(in)    :: Dt   !< Time step.
+
+   self%Dx = self%compute_dx(Dt=Dt)
+   self%Ni = nint(self%length / self%Dx)
+
+   select case(trim(adjustl(self%initial_state)))
+   case('square_wave')
+      call self%set_square_wave_initial_state
+   endselect
+   endsubroutine initialize
+
    subroutine parse_cli(self, cli)
    !< Initialize from command line interface.
-   class(integrand_ladvection),  intent(inout) :: self          !< Advection field.
-   type(command_line_interface), intent(inout) :: cli           !< Command line interface handler.
-   character(99)                               :: initial_state !< Initial state.
+   class(integrand_ladvection),  intent(inout) :: self !< Advection field.
+   type(command_line_interface), intent(inout) :: cli  !< Command line interface handler.
 
    call self%destroy
 
@@ -209,16 +259,11 @@ subroutine parse_cli(self, cli)
    call cli%get(switch='-a', val=self%a, error=cli%error) ; if (cli%error/=0) stop
    call cli%get(switch='--length', val=self%length, error=cli%error) ; if (cli%error/=0) stop
    call cli%get(switch='--Ni', val=self%Ni, error=cli%error) ; if (cli%error/=0) stop
-   call cli%get(switch='-is', val=initial_state, error=cli%error) ; if (cli%error/=0) stop
+   call cli%get(switch='-is', val=self%initial_state, error=cli%error) ; if (cli%error/=0) stop
 
    self%Ng = (self%weno_order + 1) / 2
    self%Dx = self%length / self%Ni
 
-   select case(trim(adjustl(initial_state)))
-   case('square_wave')
-      call self%set_square_wave_initial_state
-   endselect
-
    if (self%weno_order>1) call wenoof_create(interpolator_type=trim(adjustl(self%w_scheme)), &
                                              S=self%Ng,                                      &
                                              interpolator=self%interpolator,                 &
@@ -433,6 +478,7 @@ pure subroutine assign_integrand(lhs, rhs)
       lhs%CFL        = rhs%CFL
       lhs%Ni         = rhs%Ni
       lhs%Ng         = rhs%Ng
+      lhs%length     = rhs%length
       lhs%Dx         = rhs%Dx
       lhs%a          = rhs%a
       if (allocated(rhs%u)) then
@@ -447,6 +493,7 @@ pure subroutine assign_integrand(lhs, rhs)
       else
          if (allocated(lhs%interpolator)) deallocate(lhs%interpolator)
       endif
+      lhs%initial_state = rhs%initial_state
    endselect
    endsubroutine assign_integrand
 

src/tests/tester/foodie_test_integrand_oscillation.f90

@@ -50,12 +50,13 @@ 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 :: 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 :: error          !< Return error.
       procedure, pass(self), public :: exact_solution !< Return exact solution.
       procedure, pass(self), public :: export_tecplot !< Export integrand to Tecplot file.
+      procedure, pass(self), public :: initialize     !< Initialize field.
       procedure, pass(self), public :: parse_cli      !< Initialize from command line interface.
       procedure, nopass,     public :: set_cli        !< Set command line interface.
       ! integrand_object deferred methods
@@ -99,17 +100,6 @@ function amplitude_phase(self) result(ap)
    ap(2) = atan(-self%U(1) / self%U(2))
    endfunction amplitude_phase
 
-   pure subroutine initialize(self, U0, frequency)
-   !< Initialize integrand.
-   class(integrand_oscillation), intent(inout) :: self      !< Integrand.
-   real(R_P),                    intent(in)    :: U0(1:2)   !< Initial state.
-   real(R_P),                    intent(in)    :: frequency !< Frequency of oscillation.
-
-   self%U = U0
-   self%f = frequency
-   self%U0 = U0
-   endsubroutine initialize
-
    pure function output(self) result(state)
    !< Extract integrand state field.
    class(integrand_oscillation), intent(in) :: self       !< Integrand.
@@ -130,6 +120,17 @@ pure function description(self, prefix) result(desc)
    desc = prefix//'oscillation'
    endfunction description
 
+   pure function error(self, t, U0)
+   !< Return error.
+   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                             :: error(:) !< Error.
+
+   allocate(error(1:2))
+   error = abs(self%U - self%exact_solution(t=t))
+   endfunction error
+
    pure function exact_solution(self, t, U0) result(exact)
    !< Return exact solution.
    class(integrand_oscillation), intent(in)           :: self     !< Integrand.
@@ -173,6 +174,14 @@ subroutine export_tecplot(self, file_name, t, scheme, close_file)
    endif
    endsubroutine export_tecplot
 
+   pure subroutine initialize(self, Dt)
+   !< Initialize integrand.
+   !<
+   !< Intentionally empty, all is done in `parse_cli` method.
+   class(integrand_oscillation), intent(inout) :: self !< Integrand.
+   real(R_P),                    intent(in)    :: Dt   !< Time step.
+   endsubroutine initialize
+
    subroutine parse_cli(self, cli)
    !< Initialize from command line interface.
    class(integrand_oscillation), intent(inout) :: self !< Advection field.

src/tests/tester/foodie_test_integrand_tester_object.f90

@@ -17,8 +17,10 @@ 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(error_interface),          pass(self), deferred :: error          !< Return error.
       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(initialize_interface),     pass(self), deferred :: initialize     !< Initialize integrand.
       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.
 endtype integrand_tester_object
@@ -33,6 +35,15 @@ pure function description_interface(self, prefix) result(desc)
    character(len=:), allocatable                        :: desc   !< Description.
    endfunction description_interface
 
+   pure function error_interface(self, t, U0) result(error)
+   !< Return error.
+   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                               :: error(:) !< Error.
+   endfunction error_interface
+
    pure function exact_solution_interface(self, t, U0) result(exact)
    !< Return exact solution.
    import :: integrand_object, integrand_tester_object, R_P
@@ -45,17 +56,24 @@ pure function exact_solution_interface(self, t, U0) result(exact)
    subroutine export_tecplot_interface(self, file_name, t, scheme, close_file)
    !< Export integrand to Tecplot file.
    import :: integrand_tester_object, R_P
-   class(integrand_tester_object), intent(in)           :: self       !< Advection field.
+   class(integrand_tester_object), intent(in)           :: self       !< Integrand.
    character(*),                   intent(in), optional :: file_name  !< File name.
    real(R_P),                      intent(in), optional :: t          !< Time.
    character(*),                   intent(in), optional :: scheme     !< Scheme used to integrate integrand.
    logical,                        intent(in), optional :: close_file !< Flag for closing file.
    endsubroutine export_tecplot_interface
 
+   subroutine initialize_interface(self, Dt)
+   !< Initialize integrand.
+   import :: integrand_tester_object, R_P
+   class(integrand_tester_object), intent(inout) :: self !< Integrand.
+   real(R_P),                      intent(in)    :: Dt   !< Time step.
+   endsubroutine initialize_interface
+
    subroutine parse_cli_interface(self, cli)
    !< Initialize from command line interface.
    import :: command_line_interface, integrand_tester_object
-   class(integrand_tester_object), intent(inout) :: self !< Advection field.
+   class(integrand_tester_object), intent(inout) :: self !< Integrand.
    type(command_line_interface),   intent(inout) :: cli  !< Command line interface handler.
    endsubroutine parse_cli_interface
 

src/tests/tester/foodie_tester.f90

@@ -86,11 +86,12 @@ 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))
+   allocate(error(1:size(self%integrand_0%error(t=0._R_P), dim=1), 1:size(self%Dt, dim=1)))
+   if (size(self%Dt, dim=1) > 1) allocate(order(1:size(error, dim=1), 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 self%integrand_0%initialize(Dt=self%Dt(t))
          call integrate(scheme=trim(integrator_schemes(s)),          &
                         integrand_0=self%integrand_0,                &
                         Dt=self%Dt(t),                               &
@@ -102,10 +103,10 @@ subroutine execute(self)
                         output_base_name=trim(adjustl(self%output)), &
                         save_frequency=self%save_frequency,          &
                         error=error(:,t))
-         print*, 'Dt = ', self%Dt(t), ', error (max) = ', maxval(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,F6.2)', ' Observed order (min) =', minval(order(:,t-1))
+            print '(A,'//trim(str(size(order, dim=1), no_sign=.true.))//'F6.2)', ' Observed order =', order(:,t-1)
          endif
       enddo
    enddo
@@ -209,15 +210,15 @@ function is_dt_valid()
       is_dt_valid = .true.
       do t=1, size(self%Dt)
          is_dt_valid = ((self%final_time - int(self%final_time/self%Dt(t), I_P)*self%Dt(t))==0)
-         if (is_dt_valid) then
-            select type(integrand=>self%integrand_0)
-            type is(integrand_ladvection)
-               is_dt_valid = is_dt_valid .and. self%Dt(t) <= integrand%Dt(final_time=self%final_time)
-               if (.not.is_dt_valid) then
-                  write(stderr, '(A)') 'error: Dt violates CFL condition, Dt_max = '//str(integrand%Dt(final_time=self%final_time))
-               endif
-            endselect
-         endif
+         ! if (is_dt_valid) then
+         !    select type(integrand=>self%integrand_0)
+         !    type is(integrand_ladvection)
+         !       is_dt_valid = is_dt_valid .and. self%Dt(t) <= integrand%Dt(final_time=self%final_time)
+         !       if (.not.is_dt_valid) then
+         !          write(stderr, '(A)') 'error: Dt violates CFL condition, Dt_max = '//str(integrand%Dt(final_time=self%final_time))
+         !       endif
+         !    endselect
+         ! endif
          if (.not.is_dt_valid) exit
       enddo
       endfunction is_dt_valid
@@ -249,6 +250,7 @@ subroutine integrate(scheme, integrand_0, Dt, final_time, iterations, stages, is
    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.
+   real(R_P), allocatable                       :: 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.
@@ -330,7 +332,8 @@ subroutine integrate(scheme, integrand_0, Dt, final_time, iterations, stages, is
 
    call integrand_close_tecplot
 
-   call integrand_compute_error
+   error_ = integrand%error(t=time, U0=integrand_0)
+   error(:) = error_(:)
    contains
       subroutine integrand_close_tecplot
       !< Close current integrand tecplot file.
@@ -344,17 +347,6 @@ subroutine integrand_close_tecplot
       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)
-         error = abs(integrand%output() - integrand%exact_solution(t=time))
-      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.