add fast mode to many inegrators

Author: szaghi

src/lib/foodie_integrator_lmm_ssp_vss.f90

@@ -53,16 +53,17 @@ module foodie_integrator_lmm_ssp_vss
                                                          trim(class_name_)//'_steps_4_order_3', &
                                                          trim(class_name_)//'_steps_5_order_3'] !< 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_lmm_ssp_vss
   !< FOODIE integrator: provide an explicit class of Linear Multi-step Methods (LLM) with Strong Stability Preserving property and
   !< variable stepsize (VSS), from 2nd to 3rd order accurate.
   !<
   !< @note The integrator must be created or initialized before used.
   private
-  integer(I_P), public                    :: steps=0                         !< Number of time steps.
-  procedure(integrate_interface), pointer :: integrate_ => integrate_order_2 !< Integrate integrand field.
+  integer(I_P), public                         :: steps=0                                   !< Number of time steps.
+  procedure(integrate_interface),      pointer :: integrate_ => integrate_order_2           !< Integrate integrand field.
+  procedure(integrate_fast_interface), pointer :: integrate_fast_ => integrate_order_2_fast !< Integrate integrand field, fast.
   contains
     ! deferred methods
     procedure, pass(self) :: class_name           !< Return the class name of schemes.
@@ -75,10 +76,13 @@ module foodie_integrator_lmm_ssp_vss
     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.
     ! private methods
-    procedure, pass(self), private :: integrate_order_2 !< Integrate integrand field by 2nd order formula.
-    procedure, pass(self), private :: integrate_order_3 !< Integrate integrand field by 3rd order formula.
+    procedure, pass(self), private :: integrate_order_2      !< Integrate integrand field by 2nd order formula.
+    procedure, pass(self), private :: integrate_order_3      !< Integrate integrand field by 3rd order formula.
+    procedure, pass(self), private :: integrate_order_2_fast !< Integrate integrand field by 2nd order formula, fast mode.
+    procedure, pass(self), private :: integrate_order_3_fast !< Integrate integrand field by 3rd order formula, fast mode.
 endtype integrator_lmm_ssp_vss
 
 abstract interface
@@ -93,6 +97,18 @@ subroutine integrate_interface(self, U, previous, Dt, t, autoupdate)
   real(R_P),                     intent(in)    :: t(:)         !< Times.
   logical, optional,             intent(in)    :: autoupdate   !< Perform cyclic autoupdate of previous time steps.
   endsubroutine integrate_interface
+
+  subroutine integrate_fast_interface(self, U, previous, buffer, Dt, t, autoupdate)
+  !< Integrate field with LMM-SSP class scheme.
+  import :: integrand_object, integrator_lmm_ssp_vss, R_P
+  class(integrator_lmm_ssp_vss), 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(inout) :: Dt(:)        !< Time steps.
+  real(R_P),                     intent(in)    :: t(:)         !< Times.
+  logical, optional,             intent(in)    :: autoupdate   !< Perform cyclic autoupdate of previous time steps.
+  endsubroutine integrate_fast_interface
 endinterface
 
 contains
@@ -195,18 +211,23 @@ subroutine initialize(self, scheme, stop_on_fail)
     case('lmm_ssp_vss_steps_2_order_2')
       self%steps = 2
       self%integrate_ => integrate_order_2
+      self%integrate_fast_ => integrate_order_2_fast
     case('lmm_ssp_vss_steps_3_order_2')
       self%steps = 3
       self%integrate_ => integrate_order_2
+      self%integrate_fast_ => integrate_order_2_fast
     case('lmm_ssp_vss_steps_3_order_3')
       self%steps = 3
       self%integrate_ => integrate_order_3
+      self%integrate_fast_ => integrate_order_3_fast
     case('lmm_ssp_vss_steps_4_order_3')
       self%steps = 4
       self%integrate_ => integrate_order_3
+      self%integrate_fast_ => integrate_order_3_fast
     case('lmm_ssp_vss_steps_5_order_3')
       self%steps = 5
       self%integrate_ => integrate_order_3
+      self%integrate_fast_ => integrate_order_3_fast
     endselect
   else
     call self%trigger_error(error=ERROR_UNSUPPORTED_SCHEME,                                   &
@@ -227,6 +248,19 @@ subroutine integrate(self, U, previous, Dt, t, autoupdate)
   call self%integrate_(U=U, previous=previous, Dt=Dt, t=t, autoupdate=autoupdate)
   endsubroutine integrate
 
+  subroutine integrate_fast(self, U, previous, buffer, Dt, t, autoupdate)
+  !< Integrate field with LMM-SSP class scheme.
+  class(integrator_lmm_ssp_vss), 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(inout) :: Dt(:)        !< Time steps.
+  real(R_P),                     intent(in)    :: t(:)         !< Times.
+  logical, optional,             intent(in)    :: autoupdate   !< Perform cyclic autoupdate of previous time steps.
+
+  call self%integrate_fast_(U=U, previous=previous, buffer=buffer, Dt=Dt, t=t, autoupdate=autoupdate)
+  endsubroutine integrate_fast
+
   subroutine update_previous(self, U, previous, Dt)
   !< Cyclic update previous time steps.
   class(integrator_lmm_ssp_vss), intent(in)    :: self         !< Integrator.
@@ -283,6 +317,54 @@ subroutine integrate_order_3(self, U, previous, Dt, t, autoupdate)
   if (autoupdate_) call self%update_previous(U=U, previous=previous, Dt=Dt)
   endsubroutine integrate_order_3
 
+  subroutine integrate_order_2_fast(self, U, previous, buffer, Dt, t, autoupdate)
+  !< Integrate field with LMM-SSP-VSS 2nd order class scheme.
+  class(integrator_lmm_ssp_vss), 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(inout) :: Dt(:)        !< Time steps.
+  real(R_P),                     intent(in)    :: t(:)         !< Times.
+  logical, optional,             intent(in)    :: autoupdate   !< Perform cyclic autoupdate of previous time steps.
+  logical                                      :: autoupdate_  !< Perform cyclic autoupdate of previous time steps, dummy var.
+  real(R_P)                                    :: omega_       !< Omega coefficient.
+  real(R_P)                                    :: omega_sq     !< Square of omega coefficient.
+
+  autoupdate_ = .true. ; if (present(autoupdate)) autoupdate_ = autoupdate
+  omega_ = omega(Dt=Dt, s=self%steps-1)
+  omega_sq = omega_ * omega_
+  call U%multiply_fast(lhs=previous(1), rhs=1._R_P / omega_sq)
+  call buffer%multiply_fast(lhs=previous(self%steps), rhs=(omega_sq - 1._R_P) / omega_sq)
+  call U%add_fast(lhs=U, rhs=buffer)
+  buffer = previous(self%steps)
+  call buffer%t_fast(t=t(self%steps))
+  call buffer%multiply_fast(lhs=buffer, rhs=Dt(self%steps) * (omega_ + 1._R_P) / omega_)
+  call U%add_fast(lhs=U, rhs=buffer)
+  if (autoupdate_) call self%update_previous(U=U, previous=previous, Dt=Dt)
+  endsubroutine integrate_order_2_fast
+
+  subroutine integrate_order_3_fast(self, U, previous, buffer, Dt, t, autoupdate)
+  !< Integrate field with LMM-SSP-VSS 3rd order class scheme, fast mode.
+  class(integrator_lmm_ssp_vss), 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(inout) :: Dt(:)        !< Time steps.
+  real(R_P),                     intent(in)    :: t(:)         !< Times.
+  logical, optional,             intent(in)    :: autoupdate   !< Perform cyclic autoupdate of previous time steps.
+  logical                                      :: autoupdate_  !< Perform cyclic autoupdate of previous time steps, dummy var.
+  real(R_P)                                    :: omega_       !< Omega coefficient.
+
+  autoupdate_ = .true. ; if (present(autoupdate)) autoupdate_ = autoupdate
+  omega_= omega(Dt=Dt, s=self%steps-1)
+  ! @TODO to be completed
+  U = (previous(1) * ((3._R_P * omega_ + 2._R_P) / omega_ ** 3 )) +                           &
+      (previous(self%steps) * (((omega_ + 1._R_P) ** 2) * (omega_ - 2._R_P) / omega_ ** 3)) + &
+      (previous(1)%t(t=t(1)) * (Dt(self%steps) * (omega_ + 1._R_P) / omega_ ** 2)) +          &
+      (previous(self%steps)%t(t=t(self%steps)) * (Dt(self%steps) * (omega_ + 1._R_P) ** 2 / omega_ ** 2))
+  if (autoupdate_) call self%update_previous(U=U, previous=previous, Dt=Dt)
+  endsubroutine integrate_order_3_fast
+
   ! private non TBP
   pure function dt_ratio(Dt, s) result(ratio)
   !< Return `Dt(n+s)/Dt(n+Ns)` ratio.

src/lib/foodie_integrator_ms_runge_kutta_ssp.f90

@@ -49,7 +49,7 @@ module foodie_integrator_ms_runge_kutta_ssp
                                                          trim(class_name_)//'_steps_4_stages_5_order_8'] !< 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_ms_runge_kutta_ssp
   !< FOODIE integrator: provide an explicit class of Multi-step Runge-Kutta Methods with Strong Stability Preserving property,
@@ -77,6 +77,7 @@ module foodie_integrator_ms_runge_kutta_ssp
     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_ms_runge_kutta_ssp
 
@@ -361,6 +362,77 @@ subroutine integrate(self, U, previous, stage, Dt, t, autoupdate)
   if (autoupdate_) call self%update_previous(U=U, previous=previous)
   endsubroutine integrate
 
+  subroutine integrate_fast(self, U, previous, stage, buffer, Dt, t, autoupdate)
+  !< Integrate field with LMM-SSP class scheme, fast mode.
+  class(integrator_ms_runge_kutta_ssp), 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)        :: stage(1:)    !< Runge-Kutta stages [1:stages].
+  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.
+  logical,                              intent(in), optional :: autoupdate   !< Perform cyclic autoupdate of previous steps.
+  logical                                                    :: autoupdate_  !< Perform cyclic autoupdate of previous steps,
+                                                                             !< local variable.
+  integer(I_P)                                               :: k, kk        !< Stages counters.
+  integer(I_P)                                               :: s            !< Steps counter.
+
+  autoupdate_ = .true. ; if (present(autoupdate)) autoupdate_ = autoupdate
+  ! computing stages
+  stage(1) = U
+  do k=2, self%stages
+    call stage(k)%multiply_fast(lhs=U, rhs=0._R_P)
+    do s=1, self%steps
+      if (self%D(k, s) /= 0._R_P) then
+         call buffer%multiply_fast(lhs=previous(s), rhs=self%D(k, s))
+         call stage(k)%add_fast(lhs=stage(k), rhs=buffer)
+      endif
+    enddo
+    do s=1, self%steps - 1
+      if (self%Ahat(k, s) /= 0._R_P) then
+         buffer = previous(s)
+         call buffer%t_fast(t=t(s))
+         call buffer%multiply_fast(lhs=buffer, rhs=Dt * self%Ahat(k, s))
+         call stage(k)%add_fast(lhs=stage(k), rhs=buffer)
+      endif
+    enddo
+    do kk=1, k - 1
+      if (self%A(k, kk) /= 0._R_P) then
+         buffer = stage(kk)
+         call buffer%t_fast(t=t(self%steps))
+         call buffer%multiply_fast(lhs=buffer, rhs=Dt * self%A(k, kk))
+         call stage(k)%add_fast(lhs=stage(k), rhs=buffer)
+      endif
+    enddo
+  enddo
+  ! computing new time step
+  U = U * 0._R_P
+  do s=1, self%steps
+    if (self%T(s) /= 0._R_P) then
+      call buffer%multiply_fast(lhs=previous(s), rhs=self%T(s))
+      call U%add_fast(lhs=U, rhs=buffer)
+    endif
+  enddo
+  do s=1, self%steps - 1
+    if (self%Bhat(s) /= 0._R_P) then
+       buffer = previous(s)
+       call buffer%t_fast(t=t(s))
+       call buffer%multiply_fast(lhs=buffer, rhs=Dt * self%Bhat(s))
+       call U%add_fast(lhs=U, rhs=buffer)
+    endif
+  enddo
+  do k=1, self%stages
+    if (self%B(k) /= 0._R_P) U = U + (stage(k)%t(t=t(self%steps)) * (Dt * self%B(k)))
+    if (self%B(k) /= 0._R_P) then
+       buffer = stage(k)
+       call buffer%t_fast(t=t(self%steps))
+       call buffer%multiply_fast(lhs=buffer, rhs=Dt * self%B(k))
+       call U%add_fast(lhs=U, rhs=buffer)
+    endif
+  enddo
+  if (autoupdate_) call self%update_previous(U=U, previous=previous)
+  endsubroutine integrate_fast
+
   subroutine update_previous(self, U, previous)
   !< Cyclic update previous time steps.
   class(integrator_ms_runge_kutta_ssp), intent(in)    :: self         !< Integrator.

src/lib/foodie_integrator_runge_kutta_embedded.f90

@@ -281,7 +281,7 @@ module foodie_integrator_runge_kutta_emd
                                                          trim(class_name_)//'_stages_9_order_6  ', &
                                                          trim(class_name_)//'_stages_17_order_10'] !< 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_runge_kutta_emd
   !< FOODIE integrator: provide an explicit class of embedded Runge-Kutta schemes, from 2nd to 10th order accurate.
@@ -303,9 +303,10 @@ module foodie_integrator_runge_kutta_emd
     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) :: 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.
     ! private methods
     procedure, pass(self), private :: new_Dt !< Compute new estimation of the time step Dt.
 endtype integrator_runge_kutta_emd
@@ -677,8 +678,6 @@ subroutine integrate(self, U, stage, Dt, t)
   !< Integrate field with explicit embedded Runge-Kutta scheme.
   !<
   !< The time steps is adaptively resized using the local truncation error estimation by means of the embedded pairs of RK schemes.
-  !<
-  !< @note This method can be used **after** the integrator is created (i.e. the RK coefficients are initialized).
   class(integrator_runge_kutta_emd), intent(in)    :: self      !< Integrator.
   class(integrand_object),           intent(inout) :: U         !< Field to be integrated.
   class(integrand_object),           intent(inout) :: stage(1:) !< Runge-Kutta stages [1:stages].
@@ -715,6 +714,50 @@ subroutine integrate(self, U, stage, Dt, t)
   U = U1
   endsubroutine integrate
 
+  subroutine integrate_fast(self, U, stage, buffer, Dt, t)
+  !< Integrate field with explicit embedded Runge-Kutta scheme, fast mode.
+  !<
+  !< The time steps is adaptively resized using the local truncation error estimation by means of the embedded pairs of RK schemes.
+  class(integrator_runge_kutta_emd), intent(in)    :: self      !< Integrator.
+  class(integrand_object),           intent(inout) :: U         !< Field to be integrated.
+  class(integrand_object),           intent(inout) :: stage(1:) !< Runge-Kutta stages [1:stages].
+  class(integrand_object),           intent(inout) :: buffer    !< Temporary buffer for doing fast operation.
+  real(R_P),                         intent(inout) :: Dt        !< Time step.
+  real(R_P),                         intent(in)    :: t         !< Time.
+  class(integrand_object), allocatable             :: U1        !< First U evaluation.
+  class(integrand_object), allocatable             :: U2        !< Second U evaluation.
+  real(R_P)                                        :: error     !< Local truncation error estimation.
+  integer(I_P)                                     :: s         !< First stages counter.
+  integer(I_P)                                     :: ss        !< Second stages counter.
+
+  allocate(U1, mold=U) ; U1 = U
+  allocate(U2, mold=U) ; U2 = U
+  error = 1e6
+  do while(error>self%tolerance)
+    ! compute stages
+    do s=1, self%stages
+      stage(s) = U
+      do ss=1, s - 1
+        call buffer%multiply_fast(lhs=stage(ss), rhs=Dt * self%alph(s, ss))
+        call stage(s)%add_fast(lhs=stage(s), rhs=buffer)
+      enddo
+      call stage(s)%t_fast(t=t + self%gamm(s) * Dt)
+    enddo
+    ! compute new time step
+    U1 = U
+    U2 = U
+    do s=1, self%stages
+      call buffer%multiply_fast(lhs=stage(s), rhs=Dt * self%beta(s, 1))
+      call U1%add_fast(lhs=U1, rhs=buffer)
+      call buffer%multiply_fast(lhs=stage(s), rhs=Dt * self%beta(s, 2))
+      call U2%add_fast(lhs=U2, rhs=buffer)
+    enddo
+    error = U2.lterror.U1
+    call self%new_Dt(error=error, Dt=Dt)
+  enddo
+  U = U1
+  endsubroutine integrate_fast
+
   ! private methods
   elemental subroutine new_Dt(self, error, Dt)
   !< Compute new estimation of the time step Dt.