!>
!> @file test_mgp.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
!
!   Test multigrid V-cycle for periodic problems
!
  USE multigrid
  IMPLICIT NONE
!
  INTEGER          :: nx=8, nidbas=1, ngauss=2, nits=40
  INTEGER          :: levels=2, nu1=1, nu2=1, mu=1, nu0=1
  CHARACTER(len=4) :: relax='jac '
  LOGICAL          :: nlfixed = .FALSE.
  DOUBLE PRECISION :: sigma=1.0d0, kmode=10.0, pi=4.0d0*ATAN(1.0d0)
  DOUBLE PRECISION :: omega=2.0d0/3.0d0
  INTEGER          :: l, nrank, dim, its
  DOUBLE PRECISION :: errdisc_dir
  DOUBLE PRECISION, ALLOCATABLE :: u_direct(:), u_exact(:), u_calc(:)
  DOUBLE PRECISION, ALLOCATABLE :: sol_direct(:), sol_calc(:), sol_grid(:)
  DOUBLE PRECISION, ALLOCATABLE :: err(:), resid(:), errdisc(:)
  DOUBLE PRECISION, ALLOCATABLE :: errdisc_fmg(:)
!
  TYPE(grid1d), ALLOCATABLE :: gridx(:)
  TYPE(mg_info)             :: info
!
  NAMELIST /newrun/ nx, nidbas, ngauss, sigma, kmode, &
       &            relax, nits, nlfixed, levels, nu1, nu2, mu, nu0
!--------------------------------------------------------------------------------
!                   1.  Prologue
!   Inputs
!
  READ(*,newrun)
  WRITE(*,newrun)
!
  levels = MIN(levels, get_lmax(nx))
!
  info%nu1 = nu1
  info%nu2 = nu2
  info%mu = mu
  info%nu0 = nu0
  info%levels = levels
  info%relax = relax
  info%omega = omega
!
!   Create grids
!
  ALLOCATE(gridx(levels))
  CALL create_grid(nx, nidbas, ngauss, 0, gridx, period=.TRUE.)
  WRITE(*,'(a/(20i6))') 'Number of intervals in grids', (gridx(l)%n, l=1,levels)
!
!   Create FE matrice and set BC u(0)=u(1)=0
!
  DO l=1,levels
     CALL femat(gridx(l)%spl, gridx(l)%matap, coefeq)
  END DO
!
!   Construct RHS only on the finest grid
!
  nrank = gridx(1)%rank    !  Rank of the system (number of unknowns)
  dim = nrank+nidbas       !  Dimension of Splines space
  CALL disrhs(gridx(1)%spl, gridx(1)%f, rhs)
!--------------------------------------------------------------------------------
!                   2.  Direct solution
!
  WRITE(*,'(//a)') 'Direct solution for the finest grid problem'
  ALLOCATE(u_direct(0:nx))
  ALLOCATE(sol_direct(nrank))
  ALLOCATE(sol_grid(dim))      ! Required by GRIDVAL
!
  CALL direct_solve(gridx(1), sol_direct)
  sol_grid(1:nrank) = sol_direct(1:nrank)
  sol_grid(nrank+1:dim) = sol_direct(1:nidbas)
  CALL gridval(gridx(1)%spl, gridx(1)%x, u_direct, 0, sol_grid)
!
  errdisc_dir = disc_err(gridx(1)%spl, sol_grid, sol)
  WRITE(*,'(a,1pe12.3)') 'Discretization error', errdisc_dir
!--------------------------------------------------------------------------------
!                   3.  Solution from MG V-cycles
!
  WRITE(*,'(//a)') 'Multigrid MG V-cycles'
  ALLOCATE(sol_calc(nrank))
  ALLOCATE(err(0:nits))
  ALLOCATE(errdisc(0:nits))
  ALLOCATE(resid(0:nits))
!
!   Initial guess
!
  sol_calc(:) = 0.0d0
  sol_grid(:) = 0.0d0
  IF(nlfixed) THEN
     sol_calc(:) = sol_direct(:)
     sol_grid(1:nrank) = sol_calc(1:nrank)
     sol_grid(nrank+1:dim) = sol_calc(1:nidbas)
  END IF
  gridx(1)%v(:) = sol_calc(:)
  err(0) = normf(gridx(1)%matmp, sol_calc-sol_direct)
  errdisc(0) = disc_err(gridx(1)%spl, sol_grid, sol)
  resid(0) = residue(gridx(1)%matap, gridx(1)%f, sol_calc)
!
!   Iterations
!
  DO its=1,nits
!
     CALL mg(gridx, info, 1)
     sol_calc(:) = gridx(1)%v(:)
     sol_grid(1:nrank) = sol_calc(1:nrank)
     sol_grid(nrank+1:dim) = sol_calc(1:nidbas)
!
     err(its) = normf(gridx(1)%matmp, sol_calc-sol_direct)
     errdisc(its) = disc_err(gridx(1)%spl, sol_grid, sol)  ! will call GRIDVAL
     resid(its) = residue(gridx(1)%matap, gridx(1)%f, sol_calc)
  END DO
!
  WRITE(*,'(a4,3(a12,a8))') 'its', 'error', 'ratio', 'residue', 'ratio', &
       &                'disc. err', 'ratio'
  WRITE(*,'(i4,3(1pe12.3,8x))') 0, err(0), resid(0), errdisc(0)
  WRITE(*,'((i4,3(1pe12.3,0pf8.2)))') (its, err(its), err(its)/err(its-1), &
       &   resid(its), resid(its)/resid(its-1), &
       &   errdisc(its), errdisc(its)/errdisc(its-1), its=1,nits)
!--------------------------------------------------------------------------------
!                   4.  Solution from FMG
!
  WRITE(*,'(//a)') 'Full Multigrid'

  ALLOCATE(errdisc_fmg(nits))
  DO its=1,nits
     info%nu0 = its
!
     CALL fmg(gridx, info, 1)
     sol_calc(:) = gridx(1)%v(:)
     sol_grid(1:nrank) = sol_calc(1:nrank)
     sol_grid(nrank+1:dim) = sol_calc(1:nidbas)
!
     errdisc_fmg(its) = disc_err(gridx(1)%spl, sol_grid, sol)  ! will call GRIDVAL
     resid(its) = residue(gridx(1)%matap, gridx(1)%f, sol_calc)
  END DO
  WRITE(*,'(a4,2(a12,a8))') 'nu0', 'residue', 'ratio','disc. err', 'ratio'
  WRITE(*,'((i4,2(1pe12.3,0pf8.3)))') (its, resid(its), resid(its)/resid(its-1), &
       &     errdisc_fmg(its),errdisc_fmg(its)/errdisc_dir, its=1,nits)
!
!   Grid values at final iteration
!
  ALLOCATE(u_exact(0:nx))
  ALLOCATE(u_calc(0:nx))
  u_exact = sol(gridx(1)%x)
  CALL gridval(gridx(1)%spl, gridx(1)%x, u_calc, 0, sol_grid)
!--------------------------------------------------------------------------------
!                   9.  Epilogue
!
!   Creata HDF5 file
!
  CALL h5file
!--------------------------------------------------------------------------------
CONTAINS
  FUNCTION rhs(x)
    DOUBLE PRECISION, INTENT(in) :: x
    DOUBLE PRECISION :: rhs
    rhs = SIN(pi*kmode*x)
  END FUNCTION rhs
  FUNCTION sol(x)
    DOUBLE PRECISION, INTENT(in) :: x(:)
    DOUBLE PRECISION :: sol(SIZE(x))
    sol(:) = 1.0d0/((pi*kmode)**2+sigma)*SIN(pi*kmode*x(:))
  END FUNCTION sol
  SUBROUTINE coefeq(x, idt, idw, c)
    DOUBLE PRECISION, INTENT(in) :: x
    INTEGER, INTENT(out) :: idt(:), idw(:)
    DOUBLE PRECISION, INTENT(out) :: c(:)
       c(1) = 1.0d0
       idt(1) = 1
       idw(1) = 1
       c(2) = sigma
       idt(2) = 0
       idw(2) = 0
  END SUBROUTINE coefeq
  SUBROUTINE h5file
    USE futils
    CHARACTER(len=128) :: file='test_mgp.h5'
    INTEGER :: fid
    INTEGER :: l
    CHARACTER(len=64) :: dsname
    CALL creatf(file, fid, real_prec='d')
    CALL attach(fid, '/', 'NX', nx)
    CALL attach(fid, '/', 'NIDBAS', nidbas)
    CALL attach(fid, '/', 'SIGMA', sigma)
    CALL attach(fid, '/', 'KMODE', kmode)
    CALL attach(fid, '/', 'RELAX', relax)
    CALL attach(fid, '/', 'NITS', nits)
    CALL attach(fid, '/', 'NLFIXED', nlfixed)
    CALL attach(fid, '/', 'LEVELS', levels)
    CALL attach(fid, '/', 'NU1', nu1)
    CALL attach(fid, '/', 'NU2', nu2)
    CALL attach(fid, '/', 'NU0', nu0)
    CALL attach(fid, '/', 'MU', mu)
    CALL creatg(fid, '/mglevels')
    DO l=1,levels
       WRITE(dsname,'("/mglevels/level.",i2.2)') l
       CALL creatg(fid, TRIM(dsname))
       CALL putarr(fid, TRIM(dsname)//'/mata', gridx(l)%matap%val)
       IF(l.GT.1) THEN
          CALL putarr(fid, TRIM(dsname)//'/matp', gridx(l)%transf%val)
          CALL attach(fid, TRIM(dsname)//'/matp', 'M', gridx(l)%transf%mrows)
          CALL attach(fid, TRIM(dsname)//'/matp', 'N', gridx(l)%transf%ncols)
       END IF
       CALL putarr(fid, TRIM(dsname)//'/f', gridx(l)%f)
       CALL putarr(fid, TRIM(dsname)//'/v', gridx(l)%v)
    END DO
    CALL creatg(fid, '/Iterations')
    CALL putarr(fid, '/Iterations/errors', err)
    CALL putarr(fid, '/Iterations/residues', resid)
    CALL putarr(fid, '/Iterations/disc_errors', errdisc)
    CALL putarr(fid, '/Iterations/disc_errors_fmg', errdisc_fmg)
    CALL putarr(fid, '/Iterations/xgrid', gridx(1)%x)
    CALL putarr(fid, '/Iterations/u_direct', u_direct)
    CALL putarr(fid, '/Iterations/u_exact', u_exact)
    CALL putarr(fid, '/Iterations/u_calc', u_calc)
    CALL closef(fid)
  END SUBROUTINE h5file
END PROGRAM main