!>
!> @file test_intergrid0.f90
!>
!> @brief
!>
!> @copyright
!> Copyright (©) 2021 EPFL (Ecole Polytechnique Fédérale de Lausanne)
!> SPC (Swiss Plasma Center)
!>
!> SPClibs is free software: you can redistribute it and/or modify it under
!> the terms of the GNU Lesser General Public License as published by the Free
!> Software Foundation, either version 3 of the License, or (at your option)
!> any later version.
!>
!> SPClibs is distributed in the hope that it will be useful, but WITHOUT ANY
!> WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
!> FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
!>
!> You should have received a copy of the GNU Lesser General Public License
!> along with this program. If not, see <https://www.gnu.org/licenses/>.
!>
!> @authors
!> (in alphabetical order)
!> @author Trach-Minh Tran <trach-minh.tran@epfl.ch>
!>
PROGRAM main
!
!   intergrid transfer using *serial* multigrid module:
!     - restriction of rhs
!     - prolongation of sol
!
  USE multigrid, ONLY : grid2d, mg_info, &
       &                get_lmax, create_grid_fd, ibc_transf, &
       &                prolong, restrict
  IMPLICIT NONE
  DOUBLE PRECISION, PARAMETER :: pi=4.0d0*ATAN(1.0d0)
  DOUBLE PRECISION :: Lx, Ly, kx, ky, icrosst, beta, miome
  INTEGER :: nx, ny, levels
  CHARACTER(len=4) :: prb
  LOGICAL :: nldebug
!
  DOUBLE PRECISION :: dx, dy
  DOUBLE PRECISION, ALLOCATABLE :: x(:),y(:)
!
  TYPE(mg_info) :: info  ! info for MG
  TYPE(grid2d), ALLOCATABLE :: grids(:)
!
  INTEGER :: i, l
!
  NAMELIST /parameters/ prb, nx, ny, levels, Lx, Ly, kx, ky, icrosst, beta, &
       &                miome, nldebug
!--------------------------------------------------------------------------------
!
!   Default inputs
!
  nx=32
  ny=32
  levels = 2
  kx=1
  ky=1
  icrosst=1.0d0
  Lx = 1.0D0
  Ly = 1.0D0
  miome = 200d0
  beta = 0d0
  prb = 'dddd'
  nldebug = .FALSE.
!
  READ(*,parameters)
  WRITE(*,parameters)
!
!   Fine grid
!
  dx = lx/REAL(nx,8)
  dy = ly/REAL(ny,8)
  ALLOCATE(x(0:nx), y(0:ny))
  x = dx * [(i,i=0,nx)]
  y = dy * [(i,i=0,ny)]
  WRITE(*,'(a/10(1pe12.3))') 'x =', x
  WRITE(*,'(a/10(1pe12.3))') 'y =', y
!
!   Create array of grids
!
  levels = MIN(levels, get_lmax(nx), get_lmax(ny))
  WRITE(*,'(a,i4)') 'Number of levels', levels
  ALLOCATE(grids(levels))
  info%nu1 = 1
  info%nu2 = 1
  info%mu = 1
  info%nu0 = 1
  info%levels = levels
  info%relax = 'jac'
  info%omega = 1
  CALL create_grid_fd(x, y, grids, info, mat_type='cds', debug=nldebug)
!
!   Set BC on grid transfer matrices
!
  IF(prb.EQ.'dddd') CALL ibc_transf(grids,1,3) ! Direction X
  CALL ibc_transf(grids,2,3)                   ! Direction Y
!
!   Define RHS at l=1, compute RHS at l=2,...,levels by "restriction".
!
  grids(1)%f(:,:) = frhs(grids(1)%x,grids(1)%y)
  DO l=2,levels
     grids(l)%f = restrict(grids(l)%matp, grids(l-1)%f)
     grids(l)%f = 0.25d0*grids(l)%f  ! Scaling for FD
  END DO
!
!   Define SOL at l=levels, compute SOL at l=levels-1,..,1 by "prolongation"
!
  grids(levels)%v(:,:) = fsol(grids(levels)%x,grids(levels)%y)
  DO l=levels-1,1,-1
     grids(l)%v = prolong(grids(l+1)%matp, grids(l+1)%v)
  END DO
!
  IF(nldebug) THEN
     DO l=1,levels
        WRITE(*,'(a,i3)') '==== Level', l
        WRITE(*,'(a)') 'f ='
        DO i=0,grids(l)%n(1)
           WRITE(*,'(10f8.3)') grids(l)%f(i,:)
        END DO
        WRITE(*,'(a)') 'v ='
        DO i=0,grids(l)%n(1)
           WRITE(*,'(10f8.3)') grids(l)%v(i,:)
        END DO
     END DO
  END IF
!
!   Epilogue
!
  CALL h5file
!--------------------------------------------------------------------------------
CONTAINS
!+++
  FUNCTION fdense(x)
!
!   Return density
!
    DOUBLE PRECISION, INTENT(in) :: x(:)
    DOUBLE PRECISION :: fdense(SIZE(x))
    fdense(:) = 0.5d0*beta*miome*EXP( -(x(:)-Lx/3d0)**2d0/(Lx/2d0)**2d0 );
  END FUNCTION fdense
!+++
  FUNCTION frhs(x,y)
!
!   Return RHS
!
    DOUBLE PRECISION, INTENT(in) :: x(:), y(:)
    DOUBLE PRECISION :: frhs(SIZE(x),SIZE(y))
    DOUBLE PRECISION :: c, s, d(SIZE(x))
    DOUBLE PRECISION :: corr
    INTEGER :: j
    corr = 1.d0+icrosst**2/4.0d0
    d(:) = fdense(x(:))
    IF(prb.EQ.'dddd') THEN
       DO j=1,SIZE(y)
          c = COS(2.0d0*pi*ky*y(j)/Ly)
          s = SIN(2.0d0*pi*ky*y(j)/Ly)
          frhs(:,j) =  -4.0d0*pi*pi*( (kx*kx/Lx/Lx+corr*ky*ky/Ly/Ly) * SIN(2.0d0*pi*kx*x(:)/Lx)*s &
               &     -icrosst*kx*ky/Lx/Ly * COS(2.0d0*pi*kx*x(:)/Lx)*c ) &
               &     + d(:) * SIN(2.0d0*pi*kx*x(:)/Lx)*s
       END DO
    ELSE IF (prb.EQ.'nndd') THEN
       DO j=1,SIZE(y)
          c = COS(2.0d0*pi*ky*y(j)/Ly)
          s = SIN(2.0d0*pi*ky*y(j)/Ly)
          frhs(:,j) =  -4.0d0*pi*pi*( (kx*kx/Lx/Lx+corr*ky*ky/Ly/Ly) * COS(2.0d0*pi*kx*x(:)/Lx)*s &
               &     +icrosst*kx*ky/Lx/Ly * SIN(2.0d0*pi*kx*x(:)/Lx)*c ) &
               &     + d(:) * COS(2.0d0*pi*kx*x(:)/Lx)*s
       END DO
    END IF
  END FUNCTION frhs
!+++
  FUNCTION fsol(x,y)
!
!   Return analytical solution
!
    DOUBLE PRECISION, INTENT(in) :: x(:), y(:)
    DOUBLE PRECISION :: fsol(SIZE(x),SIZE(y))
    DOUBLE PRECISION :: c
    INTEGER :: j
    IF(prb.EQ.'dddd') THEN
       DO j=1,SIZE(y)
          c = SIN(2.0d0*pi*ky*y(j)/Ly)
          fsol(:,j) = SIN(2.0d0*pi*kx*x(:)/Lx)*c
       END DO
    ELSE IF (prb.EQ.'nndd') THEN
       DO j=1,SIZE(y)
          c = SIN(2.0d0*pi*ky*y(j)/Ly)
          fsol(:,j) = COS(2.0d0*pi*kx*x(:)/Lx)*c
       END DO
    END IF
  END FUNCTION fsol
!+++
  SUBROUTINE h5file
    USE futils
    USE csr, ONLY : putmat
    CHARACTER(len=128) :: file='test_intergrid0.h5'
    INTEGER :: fid
    INTEGER :: l
    CHARACTER(len=64) :: dsname
    CALL creatf(file, fid, real_prec='d')
    CALL attach(fid, '/', 'NX', nx)
    CALL attach(fid, '/', 'NY', ny)
    CALL attach(fid, '/', 'KX', kx)
    CALL attach(fid, '/', 'KY', ky)
    CALL attach(fid, '/', 'LX', Lx)
    CALL attach(fid, '/', 'LY', Ly)
    CALL attach(fid, '/', 'BETA', beta)
    CALL attach(fid, '/', 'LEVELS', levels)
    CALL attach(fid, '/', 'PRB', prb)
    CALL attach(fid, '/', 'NLDEBUG', nldebug)
    CALL creatg(fid, '/mglevels')
    DO l=1,levels
       WRITE(dsname,'("/mglevels/level.",i2.2)') l
       CALL creatg(fid, TRIM(dsname))
       IF(l.GT.1) THEN
          CALL putmat(fid, TRIM(dsname)//'/matpx', grids(l)%matp(1))
          CALL putmat(fid, TRIM(dsname)//'/matpy', grids(l)%matp(2))
       END IF
       CALL putarr(fid, TRIM(dsname)//'/x', grids(l)%x)
       CALL putarr(fid, TRIM(dsname)//'/y', grids(l)%y)
       CALL putarr(fid, TRIM(dsname)//'/f', grids(l)%f)
       CALL putarr(fid, TRIM(dsname)//'/v', grids(l)%v)
    END DO
    CALL closef(fid)
  END SUBROUTINE h5file
!+++
END PROGRAM main