!< Define [[integrand_oscillation]], the Oscillation test field that is a concrete extension of the abstract integrand type.

module foodie_test_integrand_oscillation
!< Define [[integrand_oscillation]], the Oscillation test field that is a concrete extension of the abstract integrand type.

use flap, only : command_line_interface
use foodie, only : integrand_object
use foodie_test_integrand_tester_object, only : integrand_tester_object
use penf, only : FR_P, R_P, I_P, str

implicit none
private
public :: integrand_oscillation

type, extends(integrand_tester_object) :: integrand_oscillation
   !< The oscillation equations field.
   !<
   !< It is a FOODIE integrand class concrete extension.
   !<
   !<### Oscillation ODEs system
   !<The (inertial) Oscillation equations system is a non linear system of pure ODEs and it can be written as:
   !<
   !<$$\begin{matrix}
   !< U_t = R(U)  \\
   !< U = \begin{bmatrix}
   !< v_1 \\
   !< v_2
   !< \end{bmatrix}\;\;\;
   !< R(U) = \begin{bmatrix}
   !< -f v_2 \\
   !< f v_1
   !< \end{bmatrix}
   !<\end{matrix}$$
   !<
   !<The frequency *f* is constant and it is here selected as *f=10^-4*.
   !<
   !<In the space *v1-v2* the path of the oscillation must be a circle, but for the frequency selected some ODE solvers are not
   !<stable leading to a wrong path, see [1].
   !<
   !<#### Bibliography
   !<
   !<[1] *Numerical Methods for Fluid Dynamics With Applications to Geophysics*, Dale R. Durran, Springer, 2010.
   !<
   !<#### State variables organization
   !< State variables are organized as an array (rank 1) of reals of *n=2* elements.
   private
   real(R_P) :: f=0._R_P                 !< Oscillation frequency (Hz).
   real(R_P) :: U0(1:2)=[0._R_P, 0._R_P] !< Initial state.
   real(R_P) :: U(1:2)=[0._R_P, 0._R_P]  !< Integrand (state) variables.
   contains
      ! auxiliary methods
      procedure, pass(self), public :: amplitude_phase !< Return amplitude and phase of the oscillation.
      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
      procedure, pass(self), public :: integrand_dimension !< Return integrand dimension.
      procedure, pass(self), public :: t => dU_dt          !< Time derivative, residuals.
      ! operators
      procedure, pass(lhs), public :: local_error !<`||integrand_oscillation - integrand_oscillation||` operator.
      ! +
      procedure, pass(lhs), public :: integrand_add_integrand !< `+` operator.
      procedure, pass(lhs), public :: integrand_add_real      !< `+ real` operator.
      procedure, pass(rhs), public :: real_add_integrand      !< `real +` operator.
      ! *
      procedure, pass(lhs), public :: integrand_multiply_integrand   !< `*` operator.
      procedure, pass(lhs), public :: integrand_multiply_real        !< `* real` operator.
      procedure, pass(rhs), public :: real_multiply_integrand        !< `real *` operator.
      procedure, pass(lhs), public :: integrand_multiply_real_scalar !< `* real_scalar` operator.
      procedure, pass(rhs), public :: real_scalar_multiply_integrand !< `real_scalar *` operator.
      ! -
      procedure, pass(lhs), public :: integrand_sub_integrand !< `-` operator.
      procedure, pass(lhs), public :: integrand_sub_real      !< `- real` operator.
      procedure, pass(rhs), public :: real_sub_integrand      !< `real -` operator.
      ! =
      procedure, pass(lhs), public :: assign_integrand !< `=` operator.
      procedure, pass(lhs), public :: assign_real      !< `= real` operator.
      ! override fast operators
      procedure, pass(self), public :: t_fast                              !< Time derivative, residuals, fast mode.
      procedure, pass(opr),  public :: integrand_add_integrand_fast        !< `+` fast operator.
      procedure, pass(opr),  public :: integrand_multiply_integrand_fast   !< `*` fast operator.
      procedure, pass(opr),  public :: integrand_multiply_real_scalar_fast !< `* real_scalar` fast operator.
      procedure, pass(opr),  public :: integrand_subtract_integrand_fast   !< `-` fast operator.
endtype integrand_oscillation

contains
   ! auxiliary methods
   function amplitude_phase(self) result(ap)
   !< Compute amplitude and phase of the oscillation.
   class(integrand_oscillation), intent(in) :: self    !< Advection field.
   real(R_P)                                :: ap(1:2) !< Amplitude and phase.

   ap(1) = sqrt(self%U(1)**2 + self%U(2)**2)
   ap(2) = atan(-self%U(1) / self%U(2))
   endfunction amplitude_phase

   pure function output(self) result(state)
   !< Extract integrand state field.
   class(integrand_oscillation), intent(in) :: self       !< Integrand.
   real(R_P)                                :: state(1:2) !< State.

   state = self%U
   endfunction output

   ! integrand_tester_object deferred methods
   pure function description(self, prefix) result(desc)
   !< Return informative integrator description.
   class(integrand_oscillation), intent(in)           :: self    !< Integrand.
   character(*),                 intent(in), optional :: prefix  !< Prefixing string.
   character(len=:), allocatable                      :: desc    !< Description.
   character(len=:), allocatable                      :: prefix_ !< Prefixing string, local variable.

   prefix_ = '' ; if (present(prefix)) prefix_ = prefix
   desc = prefix//'oscillation'
   endfunction description

   pure function error(self, t, t0, U0)
   !< Return error.
   class(integrand_oscillation), intent(in)           :: self     !< Integrand.
   real(R_P),                    intent(in)           :: t        !< Time.
   real(R_P),                    intent(in), optional :: t0       !< Initial 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, t0, U0) result(exact)
   !< Return exact solution.
   class(integrand_oscillation), intent(in)           :: self     !< Integrand.
   real(R_P),                    intent(in)           :: t        !< Time.
   real(R_P),                    intent(in), optional :: t0       !< Initial 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, with_exact_solution, U0)
   !< Export integrand to Tecplot file.
   class(integrand_oscillation), intent(in)           :: self                 !< Advection field.
   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.
   logical,                      intent(in), optional :: with_exact_solution  !< Flag for export also exact solution.
   class(integrand_object),      intent(in), optional :: U0                   !< Initial conditions.
   logical                                            :: with_exact_solution_ !< Flag for export also exact solution, local variable.
   logical, save                                      :: is_open=.false.      !< Flag for checking if file is open.
   integer(I_P), save                                 :: file_unit            !< File unit.

   if (present(close_file)) then
      if (close_file .and. is_open) then
         close(unit=file_unit)
         is_open = .false.
      endif
   else
      with_exact_solution_ = .false. ; if (present(with_exact_solution)) with_exact_solution_ = with_exact_solution
      if (present(file_name)) then
         if (is_open) close(unit=file_unit)
         open(newunit=file_unit, file=trim(adjustl(file_name)))
         is_open = .true.
         if (with_exact_solution_) then
            write(unit=file_unit, fmt='(A)') 'VARIABLES="t" "x" "y" "amplitude" "phase" "x_e" "y_e"'
         else
            write(unit=file_unit, fmt='(A)') 'VARIABLES="t" "x" "y" "amplitude" "phase"'
         endif
      endif
      if (present(t) .and. present(scheme) .and. is_open) then
         write(unit=file_unit, fmt='(A)') 'ZONE T="'//trim(adjustl(scheme))//'"'
         if (with_exact_solution_) then
            write(unit=file_unit, fmt='(7('//FR_P//',1X))') t, self%U, self%amplitude_phase(), self%exact_solution(t=t)
         else
            write(unit=file_unit, fmt='(5('//FR_P//',1X))') t, self%U, self%amplitude_phase()
         endif
      elseif (present(t) .and. is_open) then
         if (with_exact_solution_) then
            write(unit=file_unit, fmt='(7('//FR_P//',1X))') t, self%U, self%amplitude_phase(), self%exact_solution(t=t)
         else
            write(unit=file_unit, fmt='(5('//FR_P//',1X))') t, self%U, self%amplitude_phase()
         endif
      endif
   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.
   type(command_line_interface), intent(inout) :: cli  !< Command line interface handler.

   call cli%get(group='oscillation', switch='-f', val=self%f, error=cli%error) ; if (cli%error/=0) stop
   call cli%get(group='oscillation', switch='-U0', val=self%U0, error=cli%error) ; if (cli%error/=0) stop
   self%U = self%U0
   endsubroutine parse_cli

   subroutine set_cli(cli)
   !< Set command line interface.
   type(command_line_interface), intent(inout) :: cli !< Command line interface handler.

   call cli%add_group(description='oscillation test settings', group='oscillation')
   call cli%add(group='oscillation', switch='--frequency', switch_ab='-f', help='frequency', required=.false., def='1e-4', &
                act='store')
   call cli%add(group='oscillation', switch='--U0', switch_ab='-U0', nargs='2', help='initial state', required=.false., &
                def='0.0 1.0', act='store')
   endsubroutine set_cli

   ! integrand_object deferred methods
   pure function integrand_dimension(self)
   !< return integrand dimension.
   class(integrand_oscillation), intent(in) :: self                !< integrand.
   integer(I_P)                             :: integrand_dimension !< integrand dimension.

   integrand_dimension = size(self%U, dim=1)
   endfunction integrand_dimension

   pure function dU_dt(self, t) result(dState_dt)
   !< Time derivative of integrand_oscillation field.
   class(integrand_oscillation), intent(in)           :: self         !< Integrand.
   real(R_P),                    intent(in), optional :: t            !< Time.
   real(R_P), allocatable                             :: dState_dt(:) !< Integrand time derivative.

   dState_dt = [-self%f * self%U(2), &
                 self%f * self%U(1)]
   endfunction dU_dt

   pure function local_error(lhs, rhs) result(error)
   !< Estimate local truncation error between 2 oscillation approximations.
   !<
   !< The estimation is done by norm L2 of U:
   !<
   !< $$ error = \sqrt{ \sum_i{ \frac{(lhs\%U_i - rhs\%U_i)^2}{lhs\%U_i^2} }} $$
   class(integrand_oscillation), intent(in) :: lhs   !< Left hand side.
   class(integrand_object),      intent(in) :: rhs   !< Right hand side.
   real(R_P)                                :: error !< Error estimation.
   integer(I_P)                             :: i     !< Space counter.

   select type(rhs)
   class is(integrand_oscillation)
      error = 0._R_P
      do i=1, size(lhs%U, dim=1)
         if (lhs%U(i) /= 0._R_P) then
            error = error + (lhs%U(i) - rhs%U(i)) ** 2 / lhs%U(i) ** 2
         elseif (rhs%U(i) /= 0._R_P) then
            error = error + (lhs%U(i) - rhs%U(i)) ** 2 / rhs%U(i) ** 2
         endif
      enddo
      error = sqrt(error)
   endselect
   endfunction local_error

   ! +
   pure function integrand_add_integrand(lhs, rhs) result(opr)
   !< `+` operator.
   class(integrand_oscillation), intent(in) :: lhs    !< Left hand side.
   class(integrand_object),      intent(in) :: rhs    !< Right hand side.
   real(R_P), allocatable                   :: opr(:) !< Operator result.

   select type(rhs)
   class is(integrand_oscillation)
      opr = lhs%U + rhs%U
   endselect
   endfunction integrand_add_integrand

   pure function integrand_add_real(lhs, rhs) result(opr)
   !< `+ real` operator.
   class(integrand_oscillation), intent(in) :: lhs     !< Left hand side.
   real(R_P),                    intent(in) :: rhs(1:) !< Right hand side.
   real(R_P), allocatable                   :: opr(:)  !< Operator result.

   opr = lhs%U + rhs
   endfunction integrand_add_real

   pure function real_add_integrand(lhs, rhs) result(opr)
   !< `real +` operator.
   real(R_P),                    intent(in) :: lhs(1:) !< Left hand side.
   class(integrand_oscillation), intent(in) :: rhs     !< Left hand side.
   real(R_P), allocatable                   :: opr(:)  !< Operator result.

   opr = lhs + rhs%U
   endfunction real_add_integrand

   ! *
   pure function integrand_multiply_integrand(lhs, rhs) result(opr)
   !< `*` operator.
   class(integrand_oscillation), intent(in) :: lhs    !< Left hand side.
   class(integrand_object),      intent(in) :: rhs    !< Right hand side.
   real(R_P), allocatable                   :: opr(:) !< Operator result.

   select type(rhs)
   class is(integrand_oscillation)
      opr = lhs%U * rhs%U
   endselect
   endfunction integrand_multiply_integrand

   pure function integrand_multiply_real(lhs, rhs) result(opr)
   !< `* real_scalar` operator.
   class(integrand_oscillation), intent(in) :: lhs     !< Left hand side.
   real(R_P),                    intent(in) :: rhs(1:) !< Right hand side.
   real(R_P), allocatable                   :: opr(:)  !< Operator result.

   opr = lhs%U * rhs
   endfunction integrand_multiply_real

   pure function real_multiply_integrand(lhs, rhs) result(opr)
   !< `real_scalar *` operator.
   class(integrand_oscillation), intent(in) :: rhs     !< Right hand side.
   real(R_P),                    intent(in) :: lhs(1:) !< Left hand side.
   real(R_P), allocatable                   :: opr(:)  !< Operator result.

   opr = lhs * rhs%U
   endfunction real_multiply_integrand

   pure function integrand_multiply_real_scalar(lhs, rhs) result(opr)
   !< `* real_scalar` operator.
   class(integrand_oscillation), intent(in) :: lhs    !< Left hand side.
   real(R_P),                    intent(in) :: rhs    !< Right hand side.
   real(R_P), allocatable                   :: opr(:) !< Operator result.

   opr = lhs%U * rhs
   endfunction integrand_multiply_real_scalar

   pure function real_scalar_multiply_integrand(lhs, rhs) result(opr)
   !< `real_scalar *` operator.
   real(R_P),                    intent(in) :: lhs    !< Left hand side.
   class(integrand_oscillation), intent(in) :: rhs    !< Right hand side.
   real(R_P), allocatable                   :: opr(:) !< Operator result.

   opr = lhs * rhs%U
   endfunction real_scalar_multiply_integrand

   ! -
   pure function integrand_sub_integrand(lhs, rhs) result(opr)
   !< `-` operator.
   class(integrand_oscillation), intent(in) :: lhs    !< Left hand side.
   class(integrand_object),      intent(in) :: rhs    !< Right hand side.
   real(R_P), allocatable                   :: opr(:) !< Operator result.

   select type(rhs)
   class is(integrand_oscillation)
      opr = lhs%U - rhs%U
   endselect
   endfunction integrand_sub_integrand

   pure function integrand_sub_real(lhs, rhs) result(opr)
   !< `- real` operator.
   class(integrand_oscillation), intent(in) :: lhs     !< Left hand side.
   real(R_P),                    intent(in) :: rhs(1:) !< Right hand side.
   real(R_P), allocatable                   :: opr(:)  !< Operator result.

   opr = lhs%U - rhs
   endfunction integrand_sub_real

   pure function real_sub_integrand(lhs, rhs) result(opr)
   !< `real -` operator.
   class(integrand_oscillation), intent(in) :: rhs     !< Left hand side.
   real(R_P),                    intent(in) :: lhs(1:) !< Left hand side.
   real(R_P), allocatable                   :: opr(:)  !< Operator result.

   opr = lhs - rhs%U
   endfunction real_sub_integrand

   ! =
   pure subroutine assign_integrand(lhs, rhs)
   !< `=` operator.
   class(integrand_oscillation), intent(inout) :: lhs !< Left hand side.
   class(integrand_object),      intent(in)    :: rhs !< Right hand side.

   select type(rhs)
   class is(integrand_oscillation)
      lhs%f = rhs%f
      lhs%U0 = rhs%U0
      lhs%U = rhs%U
   endselect
   endsubroutine assign_integrand

   pure subroutine assign_real(lhs, rhs)
   !< `= real` operator.
   class(integrand_oscillation), intent(inout) :: lhs     !< Left hand side.
   real(R_P),                    intent(in)    :: rhs(1:) !< Right hand side.

   lhs%U = rhs
   endsubroutine assign_real

   ! fast operators
   ! time derivative
   subroutine t_fast(self, t)
   !< Time derivative function of integrand class, i.e. the residuals function. Fast mode acting directly on self.
   class(integrand_oscillation), intent(inout)        :: self   !< Oscillation field.
   real(R_P),                    intent(in), optional :: t      !< Time.
   real(R_P)                                          :: buffer !< Temporary buffer.

   buffer = self%U(1)
   self%U(1) = - self%f * self%U(2)
   self%U(2) =   self%f * buffer
   endsubroutine t_fast

   ! +
   pure subroutine integrand_add_integrand_fast(opr, lhs, rhs)
   !< `+` fast operator.
   class(integrand_oscillation), intent(inout) :: opr !< Operator result.
   class(integrand_object),      intent(in)    :: lhs !< Left hand side.
   class(integrand_object),      intent(in)    :: rhs !< Right hand side.

   select type(lhs)
   class is(integrand_oscillation)
      select type(rhs)
      class is(integrand_oscillation)
         opr%U = lhs%U + rhs%U
      endselect
   endselect
   endsubroutine integrand_add_integrand_fast

   ! *
   pure subroutine integrand_multiply_integrand_fast(opr, lhs, rhs)
   !< `*` fast operator.
   class(integrand_oscillation), intent(inout) :: opr !< Operator result.
   class(integrand_object),      intent(in)    :: lhs !< Left hand side.
   class(integrand_object),      intent(in)    :: rhs !< Right hand side.

   select type(lhs)
   class is(integrand_oscillation)
      select type(rhs)
      class is(integrand_oscillation)
         opr%U = lhs%U * rhs%U
      endselect
   endselect
   endsubroutine integrand_multiply_integrand_fast

   pure subroutine integrand_multiply_real_scalar_fast(opr, lhs, rhs)
   !< `* real_scalar` fast operator.
   class(integrand_oscillation), intent(inout) :: opr !< Operator result.
   class(integrand_object),      intent(in)    :: lhs !< Left hand side.
   real(R_P),                    intent(in)    :: rhs !< Right hand side.

   select type(lhs)
   class is(integrand_oscillation)
      opr%U = lhs%U * rhs
   endselect
   endsubroutine integrand_multiply_real_scalar_fast

   ! -
   pure subroutine integrand_subtract_integrand_fast(opr, lhs, rhs)
   !< `-` fast operator.
   class(integrand_oscillation), intent(inout) :: opr !< Operator result.
   class(integrand_object),      intent(in)    :: lhs !< Left hand side.
   class(integrand_object),      intent(in)    :: rhs !< Right hand side.

   select type(lhs)
   class is(integrand_oscillation)
      select type(rhs)
      class is(integrand_oscillation)
         opr%U = lhs%U - rhs%U
      endselect
   endselect
   endsubroutine integrand_subtract_integrand_fast
endmodule foodie_test_integrand_oscillation
