!>
!> @file two_grid.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/>.
!>
!> @author
!> (in alphabetical order)
!> @author Trach-Minh Tran <trach-minh.tran@epfl.ch>
!>
PROGRAM main
!
!   Check some properties of grid transfer
!
  USE multigrid
  USE math_util, ONLY : root_bessj
  IMPLICIT NONE
!
  INTEGER          :: nx=8, nidbas=1, ngauss=2, alpha=0
  INTEGER          :: modem=22, modep=10
  INTEGER          :: levels=2
  INTEGER          :: l, nrank
  DOUBLE PRECISION :: sigma=1.0d0, kmode=10.0, pi=4.0d0*ATAN(1.0d0)
  DOUBLE PRECISION, ALLOCATABLE :: v_prolong(:)
!
  TYPE(grid1d)  :: gridx(2)
!
  NAMELIST /newrun/ nx, nidbas, ngauss, sigma, kmode, modem, modep, alpha
!--------------------------------------------------------------------------------
!                   1.  Prologue
!   Inputs
!
  READ(*,newrun)
  WRITE(*,newrun)
!
!   Create grids
!
  CALL create_grid(nx, nidbas, ngauss, alpha, gridx)
  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)%mata, coefeq)
!
!     Left Dirichlet BC (only for Cartesian geometry)
     IF(alpha .EQ. 0) THEN
        CALL ibcmat(1, gridx(l)%mata)
     END IF
!
!     Right Dirichlet BC
     CALL ibcmat(gridx(l)%mata%rank, gridx(l)%mata)
!
!     BC on grid transfer operator
     IF(l.GT.1) THEN
        WHERE( ABS(gridx(l)%transf%val) < 1.d-8) gridx(l)%transf%val=0.0d0
        IF(alpha .EQ. 0) gridx(l)%transf%val(2:,1)=0.0d0
        gridx(l)%transf%val(1:gridx(l-1)%rank-1,gridx(l)%rank)=0.0d0
     END IF
  END DO
!
!   Construct RHS and set BC only on the finest grid
!
  nrank = gridx(1)%rank
  CALL disrhs(gridx(1)%spl, gridx(1)%f, rhs)
!
!     Left Dirichlet BC (only for Cartesian geometry)
  IF(alpha .EQ. 0) THEN
     gridx(1)%f(1) = 0.0d0
  END IF
!
!     Right Dirichlet BC
  gridx(1)%f(nrank) = 0.0d0
!
!   RHS on coarse grid by restriction
!
  gridx(2)%f = restrict(gridx(2)%transf,gridx(1)%f)
!--------------------------------------------------------------------------------
!                   2.  Direct solutions
!
  DO l=1,levels
     CALL direct_solve(gridx(l), gridx(l)%v)
     WRITE(*,'(a,i3/(10(1pe12.3)))') 'Sol at level', l, gridx(l)%v
  END DO
!
!   Prolongation of coarse solution  
!
  ALLOCATE(v_prolong(SIZE(gridx(1)%v)))
!
  v_prolong = prolong(gridx(2)%transf, gridx(2)%v)
  WRITE(*,'(a,i3/(10(1pe12.3)))') 'Prolong. sol.', l, v_prolong
  WRITE(*,'(a,1pe12.3)') 'Error ||V_prolong-V_fine||', normf(gridx(1)%matm, v_prolong-gridx(1)%v)
!--------------------------------------------------------------------------------
!                   9.  Epilogue
!
!   Creata HDF5 file
!
  CALL h5file
!--------------------------------------------------------------------------------
CONTAINS
  SUBROUTINE h5file
    USE futils
    CHARACTER(len=128) :: file='two_grid.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, '/', 'ALPHA', alpha)
    CALL attach(fid, '/', 'LEVELS', levels)
    CALL creatg(fid, '/mglevels')
    DO l=1,levels
       WRITE(dsname,'("/mglevels/level.",i2.2)') l
       CALL creatg(fid, TRIM(dsname))
       CALL putmat(fid, TRIM(dsname)//'/mata', gridx(l)%mata)
       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 putarr(fid, '/v_prolong', v_prolong)
  END SUBROUTINE h5file
  FUNCTION rhs(x)
    DOUBLE PRECISION, INTENT(in) :: x
    DOUBLE PRECISION :: rhs
    DOUBLE PRECISION :: nump
    SELECT CASE (alpha)
    CASE(0)                  ! Cartesian geometry
       rhs = SIN(pi*kmode*x)
    CASE(1)                  ! Cylindrical
       nump = root_bessj(modem, modep)
       rhs = x * nump**2 * bessel_jn(modem, nump*x)
    END SELECT
  END FUNCTION rhs
  FUNCTION sol(x)
    DOUBLE PRECISION, INTENT(in) :: x(:)
    DOUBLE PRECISION :: sol(SIZE(x))
    DOUBLE PRECISION :: nump
    SELECT CASE (alpha)
    CASE(0)                  ! Cartesian geometry
       sol(:) = 1.0d0/((pi*kmode)**2+sigma)*SIN(pi*kmode*x(:))
    CASE(1)                  ! Cylindrical
       nump = root_bessj(modem, modep)
       sol(:) = bessel_jn(modem, nump*x(:))
    END SELECT
  END FUNCTION sol
  SUBROUTINE coefeq(x, idt, idw, c)
    DOUBLE PRECISION, INTENT(in) :: x
    INTEGER, INTENT(out) :: idt(:), idw(:)
    DOUBLE PRECISION, INTENT(out) :: c(:)
    SELECT CASE (alpha)
    CASE(0)           ! Cartesian geometry
       c(1) = 1.0d0
       idt(1) = 1
       idw(1) = 1
       c(2) = sigma
       idt(2) = 0
       idw(2) = 0
    CASE(1)           ! Cylindrical
       c(1) = x
       idt(1) = 1
       idw(1) = 1
       c(2) = REAL(modem,8)**2/x
       idt(2) = 0
       idw(2) = 0
    END SELECT
  END SUBROUTINE coefeq
END PROGRAM main