!> Anderson acceleration, it should accelerate the convergence of the Newton-like method
module anderson_mod
   implicit none

      ! data for Anderson Acceleration used for finding a better damping
   type, public :: AndersonAcceleration_t

      real, dimension(:,:,:), allocatable ::  AA_x
      real, dimension(:), allocatable ::  AA_update
      real, dimension(:), pointer :: alpha   ! for Anderson acceleration
      real, dimension(:,:), pointer :: dd   !!, dd1
      ! AA_max === > re-compilation necessary
      ! AA_max defined in  subroutine NewtonSolve,  newton.f90
      integer :: AA_max    ! Anderson acceleration, maximal number of considered vectors, ~ 4
      integer :: AA_F = 1     ! update d^{k+1} = w^{k+1} - w^{k}
      integer :: AA_w = 2    ! result of update w^{k+1}
      integer :: AA_m = 2    ! the maximal value AA_F, AA_w, ...

      integer :: AA_dF = -1  ! delta F (=update)  NOT NEVESSARY, ONLY in setAndersonInNewton_OLD

      contains

      procedure :: init => initAndersonData
      procedure :: dealloc => deallocateAndersonData

   end type AndersonAcceleration_t

   contains


!!!!!!!!!!!!!!!! ANDERSON !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   !> allocate data into the Newton solver needed for the Anderson acceleration
   !> nsize = state%nsize
   subroutine initAndersonData( this, nsize )
      class( AndersonAcceleration_t ) :: this
      integer, intent(in) :: nsize

!      this%AA_F = 1  ! update d^{k+1} = w^{k+1} - w^{k}
!      this%AA_w = 2  ! result of update w^{k+1}
!      this%AA_m = 2   ! the maximal value AA_F, ...

      ! Anderson acceleration
      ! data should be allocated before!
      !this%AA_max = 1  ! NO acceleration
      !this%AA_max = AA_max

      allocate( this%AA_x(1:this%AA_m, 0:this%AA_max, 1:nsize), source = 0.0 )
      allocate( this%AA_update( 1:nsize) )
      allocate( this%dd(1:this%AA_max, 1:this%AA_max ), source = 0.0 )
      !allocate( this%dd1(1:this%AA_max, 1:this%AA_max ) )
      allocate( this%alpha (0:this%AA_max ) )

   end subroutine initAndersonData

   subroutine deallocateAndersonData( this )
      class( AndersonAcceleration_t ) :: this

      deallocate( this%AA_x )
      deallocate( this%AA_update )
      deallocate( this%dd )
      !!deallocate( this%dd1 )
      deallocate( this%alpha  )

   end subroutine deallocateAndersonData


end module anderson_mod


!> module for the pseudo-time stepping, it should ensure the convergence for, e.g., porous media
module pseudotime_mod
   implicit none

      ! data for Anderson Acceleration used for finding a better damping
   type, public :: PseudoTimeStepping_t
      integer :: nPTS  ! maximal number of the pseudotime steppings
      real :: PTS_rho
      real, dimension(:,:), allocatable ::  PTS_x

      integer :: PTS_max    ! Anderson acceleration, maximal number of considered vectors, ~ 4
      integer :: PTS_xold = 1  
      integer :: PTS_actual = 2

      contains

      procedure :: init => initPseudoTimeData
      procedure :: dealloc => deallocatePseudoTimeData

   end type PseudoTimeStepping_t

   contains

   !> allocate data into the Newton solver needed for the Pseudotime acceleration
   !> nsize = state%nsize
   subroutine initPseudoTimeData( this, nsize )
      class( PseudoTimeStepping_t ) :: this
      integer, intent(in) :: nsize

!      this%AA_F = 1  ! update d^{k+1} = w^{k+1} - w^{k}
!      this%AA_w = 2  ! result of update w^{k+1}
!      this%AA_m = 2   ! the maximal value AA_F, ...

      ! Pseudo-time steping
      ! data should be allocated before!
      this%PTS_max = 2

      
      allocate( this%PTS_x(1:this%PTS_max, 1:nsize), source = 0.0 )

   end subroutine initPseudotimeData

   subroutine deallocatePseudotimeData( this )
      class( PseudoTimeStepping_t ) :: this

      deallocate( this%PTS_x )

   end subroutine deallocatePseudotimeData


 end module pseudotime_mod
