add fast mode to all remaining solvers

Author: szaghi

src/lib/foodie_integrator_adams_bashforth_moulton.f90

@@ -105,7 +105,7 @@ module foodie_integrator_adams_bashforth_moulton
                                                           trim(class_name_)//'_15', &
                                                           trim(class_name_)//'_16'] !< List of supported schemes.
 
-logical, parameter :: has_fast_mode_=.false. !< Flag to check if integrator provides *fast mode* integrate.
+logical, parameter :: has_fast_mode_=.true. !< Flag to check if integrator provides *fast mode* integrate.
 
 type, extends(integrator_object) :: integrator_adams_bashforth_moulton
   !< FOODIE integrator: provide an explicit class of Adams-Bashforth-Moulton multi-step schemes, from 1st to 4rd order accurate.
@@ -124,10 +124,11 @@ module foodie_integrator_adams_bashforth_moulton
     procedure, pass(self) :: is_supported         !< Return .true. if the integrator class support the given scheme.
     procedure, pass(self) :: supported_schemes    !< Return the list of supported schemes.
     ! public methods
-    procedure, pass(self) :: destroy       !< Destroy the integrator.
-    procedure, pass(self) :: initialize    !< Initialize (create) the integrator.
-    procedure, pass(self) :: integrate     !< Integrate integrand field.
-    procedure, pass(self) :: scheme_number !< Return the scheme number in the list of supported schemes.
+    procedure, pass(self) :: destroy        !< Destroy the integrator.
+    procedure, pass(self) :: initialize     !< Initialize (create) the integrator.
+    procedure, pass(self) :: integrate      !< Integrate integrand field.
+    procedure, pass(self) :: integrate_fast !< Integrate integrand field.
+    procedure, pass(self) :: scheme_number  !< Return the scheme number in the list of supported schemes.
 endtype integrator_adams_bashforth_moulton
 
 contains
@@ -254,6 +255,22 @@ subroutine integrate(self, U, previous, Dt, t, iterations)
   call self%predictor%update_previous(U=U, previous=previous)
   endsubroutine integrate
 
+  subroutine integrate_fast(self, U, previous, buffer, Dt, t, iterations)
+  !< Integrate field with Adams-Bashforth-Moulton class scheme, fast mode.
+  class(integrator_adams_bashforth_moulton), intent(in)           :: self         !< Integrator.
+  class(integrand_object),                   intent(inout)        :: U            !< Field to be integrated.
+  class(integrand_object),                   intent(inout)        :: previous(1:) !< Previous time steps solutions of integrand.
+  class(integrand_object),                   intent(inout)        :: buffer       !< Temporary buffer for doing fast operation.
+  real(R_P),                                 intent(in)           :: Dt           !< Time steps.
+  real(R_P),                                 intent(in)           :: t(:)         !< Times.
+  integer(I_P),                              intent(in), optional :: iterations   !< Fixed point iterations of AM scheme.
+
+  call self%predictor%integrate_fast(U=U, previous=previous, buffer=buffer, Dt=Dt, t=t, autoupdate=.false.)
+  call self%corrector%integrate_fast(U=U, previous=previous(2:), buffer=buffer, Dt=Dt, t=t, iterations=iterations, &
+                                     autoupdate=.false.)
+  call self%predictor%update_previous(U=U, previous=previous)
+  endsubroutine integrate_fast
+
   elemental function scheme_number(self, scheme)
   !< Return the scheme number in the list of supported schemes.
   class(integrator_adams_bashforth_moulton), intent(in) :: self          !< Integrator.

src/lib/foodie_integrator_adams_moulton.f90

@@ -54,7 +54,7 @@ module foodie_integrator_adams_moulton
                                                           trim(class_name_)//'_14', &
                                                           trim(class_name_)//'_15'] !< List of supported schemes.
 
-logical, parameter :: has_fast_mode_=.false. !< Flag to check if integrator provides *fast mode* integrate.
+logical, parameter :: has_fast_mode_=.true. !< Flag to check if integrator provides *fast mode* integrate.
 
 type, extends(integrator_object) :: integrator_adams_moulton
   !< FOODIE integrator: provide an explicit class of Adams-Moulton multi-step schemes, from 1st to 16th order accurate.
@@ -75,6 +75,7 @@ module foodie_integrator_adams_moulton
     procedure, pass(self) :: destroy         !< Destroy the integrator.
     procedure, pass(self) :: initialize      !< Initialize (create) the integrator.
     procedure, pass(self) :: integrate       !< Integrate integrand field.
+    procedure, pass(self) :: integrate_fast  !< Integrate integrand field, fast mode.
     procedure, pass(self) :: update_previous !< Cyclic update previous time steps.
 endtype integrator_adams_moulton
 
@@ -383,6 +384,58 @@ subroutine integrate(self, U, previous, Dt, t, iterations, autoupdate)
   endif
   endsubroutine integrate
 
+  subroutine integrate_fast(self, U, previous, buffer, Dt, t, iterations, autoupdate)
+  !< Integrate field with Adams-Moulton class scheme, fast mode.
+  class(integrator_adams_moulton), intent(in)           :: self         !< Integrator.
+  class(integrand_object),         intent(inout)        :: U            !< Field to be integrated.
+  class(integrand_object),         intent(inout)        :: previous(1:) !< Previous time steps solutions of integrand field.
+  class(integrand_object),         intent(inout)        :: buffer       !< Temporary buffer for doing fast operation.
+  real(R_P),                       intent(in)           :: Dt           !< Time steps.
+  real(R_P),                       intent(in)           :: t(:)         !< Times.
+  integer(I_P),                    intent(in), optional :: iterations   !< Fixed point iterations.
+  logical,                         intent(in), optional :: autoupdate   !< Cyclic autoupdate of previous time steps flag.
+  logical                                               :: autoupdate_  !< Cyclic autoupdate of previous time steps flag, dummy var.
+  class(integrand_object), allocatable                  :: delta        !< Delta RHS for fixed point iterations.
+  integer(I_P)                                          :: s            !< Steps counter.
+
+  autoupdate_ = .true. ; if (present(autoupdate)) autoupdate_ = autoupdate
+  if (self%steps>0) then
+    if (present(iterations)) then ! perform fixed point iterations
+      allocate(delta, mold=U)
+      delta = previous(self%steps)
+      do s=0, self%steps - 1
+        buffer = previous(s+1)
+        call buffer%t_fast(t=t(s+1))
+        call buffer%multiply_fast(lhs=buffer, rhs=Dt * self%b(s))
+        call delta%add_fast(lhs=delta, rhs=buffer)
+      enddo
+      do s=1, iterations
+        buffer = U
+        call buffer%t_fast(t=t(self%steps) + Dt)
+        call buffer%multiply_fast(lhs=buffer, rhs=Dt * self%b(self%steps))
+        call U%add_fast(lhs=delta, rhs=buffer)
+      enddo
+    else
+      buffer = U
+      call buffer%t_fast(t=t(self%steps) + Dt)
+      call buffer%multiply_fast(lhs=buffer, rhs=Dt * self%b(self%steps))
+      call U%add_fast(lhs=previous(self%steps), rhs=buffer)
+      do s=0, self%steps - 1
+        buffer = previous(s+1)
+        call buffer%t_fast(t=t(s+1))
+        call buffer%multiply_fast(lhs=buffer, rhs=Dt * self%b(s))
+        call U%add_fast(lhs=U, rhs=buffer)
+      enddo
+    endif
+    if (autoupdate_) call self%update_previous(U=U, previous=previous)
+  else
+    buffer = U
+    call buffer%t_fast(t=t(1))
+    call buffer%multiply_fast(lhs=buffer, rhs=Dt * self%b(0))
+    call U%add_fast(lhs=U, rhs=buffer)
+  endif
+  endsubroutine integrate_fast
+
   subroutine update_previous(self, U, previous)
   !< Cyclic update previous time steps.
   class(integrator_adams_moulton), intent(in)    :: self         !< Integrator.

src/tests/accuracy/oscillation/oscillation.f90

@@ -284,7 +284,7 @@ subroutine initialize_test
     if     (index(trim(adjustl(scheme)), trim(int_adams_bashforth_moulton%class_name())) > 0) then
        if (self%is_fast) then
           call check_scheme_has_fast_mode(scheme=trim(adjustl(scheme)), integr=int_adams_bashforth_moulton)
-          ! integrate => integrate_adams_bashforth_moulton_fast
+          integrate => integrate_adams_bashforth_moulton_fast
        else
           integrate => integrate_adams_bashforth_moulton
        endif
@@ -298,7 +298,7 @@ subroutine initialize_test
     elseif (index(trim(adjustl(scheme)), trim(int_adams_moulton%class_name())) > 0) then
        if (self%is_fast) then
           call check_scheme_has_fast_mode(scheme=trim(adjustl(scheme)), integr=int_adams_moulton)
-          ! integrate => integrate_adams_moulton_fast
+          integrate => integrate_adams_moulton_fast
        else
           integrate => integrate_adams_moulton
        endif
@@ -567,6 +567,63 @@ subroutine integrate_adams_bashforth_moulton(scheme, frequency, final_time, solu
   error = error_L2(frequency=frequency, solution=solution(:, 0:last_step))
   endsubroutine integrate_adams_bashforth_moulton
 
+  subroutine integrate_adams_bashforth_moulton_fast(scheme, frequency, final_time, solution, error, last_step, iterations, Dt, &
+                                                    tolerance, stages)
+  !< Integrate domain by means of the Adams-Bashforth-Moulton scheme.
+  character(*),           intent(in)           :: scheme        !< Selected scheme.
+  real(R_P),              intent(in)           :: frequency     !< Oscillation frequency.
+  real(R_P),              intent(in)           :: final_time    !< Final integration time.
+  real(R_P), allocatable, intent(out)          :: solution(:,:) !< Solution at each time step, X-Y.
+  real(R_P),              intent(out)          :: error(1:)     !< Error (norm L2) with respect the exact solution.
+  integer(I_P),           intent(out)          :: last_step     !< Last time step computed.
+  integer(I_P),           intent(in), optional :: iterations    !< Number of fixed point iterations.
+  real(R_P),              intent(in), optional :: Dt            !< Time step.
+  real(R_P),              intent(in), optional :: tolerance     !< Local error tolerance.
+  integer(I_P),           intent(in), optional :: stages        !< Number of stages.
+  type(integrator_adams_bashforth_moulton)     :: integrator    !< The integrator.
+  type(integrator_runge_kutta_ssp)             :: integrator_rk !< RK integrator for starting non self-starting integrators.
+  type(oscillator)                             :: domain        !< Oscillation field.
+  type(oscillator), allocatable                :: rk_stage(:)   !< Runge-Kutta stages.
+  type(oscillator), allocatable                :: previous(:)   !< Previous time steps solutions.
+  type(oscillator)                             :: buffer        !< Buffer oscillation field.
+  integer                                      :: step          !< Time steps counter.
+
+  call domain%init(initial_state=initial_state, frequency=frequency)
+
+  if (allocated(solution)) deallocate(solution) ; allocate(solution(0:space_dimension, 0:int(final_time/Dt)))
+  solution = 0.0_R_P
+  solution(1:, 0) = domain%output()
+
+  call integrator%initialize(scheme=scheme)
+  if (allocated(previous)) deallocate(previous) ; allocate(previous(1:integrator%steps))
+
+  call integrator_rk%initialize(scheme='runge_kutta_ssp_stages_5_order_4')
+  if (allocated(rk_stage)) deallocate(rk_stage) ; allocate(rk_stage(1:integrator_rk%stages))
+
+  step = 0
+  do while(solution(0, step) < final_time .and. step < ubound(solution, dim=2))
+    step = step + 1
+
+    if (integrator%steps >= step) then
+      call integrator_rk%integrate_fast(U=domain, stage=rk_stage, buffer=buffer, Dt=Dt, t=solution(0, step))
+      previous(step) = domain
+    else
+      call integrator%integrate_fast(U=domain,          &
+                                     previous=previous, &
+                                     buffer=buffer,     &
+                                     Dt=Dt,             &
+                                     t=solution(0, step-integrator%steps:step-1))
+    endif
+
+    solution(0, step) = step * Dt
+
+    solution(1:, step) = domain%output()
+  enddo
+  last_step = step
+
+  error = error_L2(frequency=frequency, solution=solution(:, 0:last_step))
+  endsubroutine integrate_adams_bashforth_moulton_fast
+
   subroutine integrate_adams_moulton(scheme, frequency, final_time, solution, error, last_step, iterations, Dt, tolerance, stages)
   !< Integrate domain by means of the Adams-Moulton scheme.
   character(*),           intent(in)           :: scheme        !< Selected scheme.
@@ -635,6 +692,78 @@ subroutine integrate_adams_moulton(scheme, frequency, final_time, solution, erro
   error = error_L2(frequency=frequency, solution=solution(:, 0:last_step))
   endsubroutine integrate_adams_moulton
 
+  subroutine integrate_adams_moulton_fast(scheme, frequency, final_time, solution, error, last_step, iterations, Dt, &
+                                          tolerance, stages)
+  !< Integrate domain by means of the Adams-Moulton scheme, fast mode.
+  character(*),           intent(in)           :: scheme        !< Selected scheme.
+  real(R_P),              intent(in)           :: frequency     !< Oscillation frequency.
+  real(R_P),              intent(in)           :: final_time    !< Final integration time.
+  real(R_P), allocatable, intent(out)          :: solution(:,:) !< Solution at each time step, X-Y.
+  real(R_P),              intent(out)          :: error(1:)     !< Error (norm L2) with respect the exact solution.
+  integer(I_P),           intent(out)          :: last_step     !< Last time step computed.
+  integer(I_P),           intent(in), optional :: iterations    !< Number of fixed point iterations.
+  real(R_P),              intent(in), optional :: Dt            !< Time step.
+  real(R_P),              intent(in), optional :: tolerance     !< Local error tolerance.
+  integer(I_P),           intent(in), optional :: stages        !< Number of stages.
+  type(integrator_adams_moulton)               :: integrator    !< The integrator.
+  type(integrator_runge_kutta_ssp)             :: integrator_rk !< RK integrator for starting non self-starting integrators.
+  type(oscillator)                             :: domain        !< Oscillation field.
+  type(oscillator), allocatable                :: rk_stage(:)   !< Runge-Kutta stages.
+  type(oscillator), allocatable                :: previous(:)   !< Previous time steps solutions.
+  type(oscillator)                             :: buffer        !< Buffer oscillation field.
+  integer                                      :: step          !< Time steps counter.
+  integer                                      :: step_offset   !< Time steps counter offset for slicing previous data array.
+
+  call domain%init(initial_state=initial_state, frequency=frequency)
+
+  if (allocated(solution)) deallocate(solution) ; allocate(solution(0:space_dimension, 0:int(final_time/Dt)))
+  solution = 0.0_R_P
+  solution(1:, 0) = domain%output()
+
+  call integrator%initialize(scheme=scheme)
+  if (allocated(previous)) deallocate(previous) ; allocate(previous(1:integrator%steps+1))
+  if (integrator%steps==0) then
+    step_offset = 1                ! for 0 step-(a convention)-solver offset is 1
+  else
+    step_offset = integrator%steps ! for >0 step-solver offset is steps
+  endif
+
+  call integrator_rk%initialize(scheme='runge_kutta_ssp_stages_5_order_4')
+  if (allocated(rk_stage)) deallocate(rk_stage) ; allocate(rk_stage(1:integrator_rk%stages))
+
+  step = 0
+  do while(solution(0, step) < final_time .and. step < ubound(solution, dim=2))
+    step = step + 1
+
+    if (integrator%steps >= step) then
+      call integrator_rk%integrate_fast(U=domain, stage=rk_stage, buffer=buffer, Dt=Dt, t=solution(0, step))
+      previous(step) = domain
+    else
+      if (iterations>1) then
+        call integrator%integrate_fast(U=domain,                              &
+                                       previous=previous,                     &
+                                       buffer=buffer,                         &
+                                       Dt=Dt,                                 &
+                                       t=solution(0,step-step_offset:step-1), &
+                                       iterations=iterations)
+      else
+        call integrator%integrate_fast(U=domain,          &
+                                       previous=previous, &
+                                       buffer=buffer,     &
+                                       Dt=Dt,             &
+                                       t=solution(0,step-step_offset:step-1))
+      endif
+    endif
+
+    solution(0, step) = step * Dt
+
+    solution(1:, step) = domain%output()
+  enddo
+  last_step = step
+
+  error = error_L2(frequency=frequency, solution=solution(:, 0:last_step))
+  endsubroutine integrate_adams_moulton_fast
+
   subroutine integrate_back_df(scheme, frequency, final_time, solution, error, last_step, iterations, Dt, tolerance, stages)
   !< Integrate domain by means of the back differentiation formula scheme.
   character(*),           intent(in)           :: scheme        !< Selected scheme.