!>
!> @file test_transf2d.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 .
!>
!> @author
!> (in alphabetical order)
!> @author Trach-Minh Tran
!>
PROGRAM main
!
! Test 2d multigrid
!
USE multigrid
USE csr
IMPLICIT NONE
INCLUDE 'mpif.h'
!
DOUBLE PRECISION, PARAMETER :: PI=4.0d0*ATAN(1.0d0)
INTEGER, DIMENSION(2) :: n, nidbas, ngauss, alpha
DOUBLE PRECISION :: kx=4.d0, ky=3.d0, sigma=10.0d0
INTEGER :: levels=1
DOUBLE PRECISION :: omega=2.0d0/3.0d0
!
DOUBLE PRECISION, ALLOCATABLE :: x(:), y(:)
DOUBLE PRECISION :: dx, dy
INTEGER :: ix, iy
!
DOUBLE PRECISION, ALLOCATABLE :: sol_direct_grid(:,:)
DOUBLE PRECISION, ALLOCATABLE :: sol_anal_grid(:,:)
DOUBLE PRECISION, ALLOCATABLE :: errdisc(:), resid(:)
!
DOUBLE PRECISION, ALLOCATABLE, TARGET :: fcoarse(:,:)
DOUBLE PRECISION, POINTER :: fcoarse_1d(:)
DOUBLE PRECISION, ALLOCATABLE, TARGET :: vfine(:,:)
DOUBLE PRECISION, POINTER :: vfine_1d(:)
DOUBLE PRECISION, ALLOCATABLE :: vfine_grid(:,:)
DOUBLE PRECISION, ALLOCATABLE :: err_restrict(:), err_prolong(:), &
& disc_err_prolong(:)
!
INTEGER :: ierr, me
INTEGER :: l, nterms
INTEGER :: its
!
TYPE(grid2d), ALLOCATABLE :: grids(:)
!
NAMELIST /newrun/ n, nidbas, ngauss, kx, ky, sigma, alpha, levels
!--------------------------------------------------------------------------------
! 1. Prologue
!
CALL mpi_init(ierr)
CALL mpi_comm_rank(MPI_COMM_WORLD, me, ierr)
!
! Inputs
!
n = (/8, 8/)
nidbas=(/3,3/)
ngauss=(/2,2/)
alpha = (/0,0/)
kx=4
ky=3
sigma=10.0d0
levels=2
!
IF(me.EQ.0) THEN
READ(*,newrun)
WRITE(*,newrun)
END IF
CALL mpi_bcast(n, 2, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
CALL mpi_bcast(nidbas, 2, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
CALL mpi_bcast(ngauss, 2, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
CALL mpi_bcast(alpha, 2, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
CALL mpi_bcast(kx, 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
CALL mpi_bcast(ky, 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
CALL mpi_bcast(sigma, 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
CALL mpi_bcast(levels, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
!
! Adjust number of levels
!
levels = MIN(levels, get_lmax(n(1)), get_lmax(n(2)))
!
! Create grids
!
dx = 1.0d0/REAL(n(1),8)
dy = 1.0d0/REAL(n(2),8)
ALLOCATE(x(0:n(1)), y(0:n(2)))
DO ix=0,n(1)
x(ix) = ix*dx
END DO
DO iy=0,n(2)
y(iy) = iy*dy
END DO
!
ALLOCATE(grids(levels))
CALL create_grid(x, y, nidbas, ngauss, alpha, grids)
WRITE(*,'(5a6)') 'l', 'nx', 'ny', 'rx', 'ry'
WRITE(*,'(5i6)') (l, grids(l)%n, grids(l)%rank, l=1,levels)
!
! Build FE matrices and set BC
!
nterms = 3
DO l=1,levels
CALL femat(grids(l)%spl, grids(l)%mata, coefeq, nterms)
CALL ibcmat(grids(l), grids(l)%mata)
CALL to_mat(grids(l)%mata)
END DO
!
! Set BC on grid transfer matrices
!
CALL ibc_transf(grids,1,3)
CALL ibc_transf(grids,2,3)
!--------------------------------------------------------------------------------
! 1. Direct solutions
!
WRITE(*,'(/a)') 'Direct solutions for all levels ...'
ALLOCATE(errdisc(levels))
ALLOCATE(resid(levels))
WRITE(*,'(3a5,2a12)') 'l', 'nx', 'ny', 'err', 'resid'
DO l=1,levels
CALL disrhs(grids(l)%spl, grids(l)%f, rhs)
CALL ibcrhs(grids(l), grids(l)%f)
grids(l)%v = grids(l)%f
CALL direct_solve(grids(l), grids(l)%v1d, debug=.FALSE.)
errdisc(l) = disc_err(grids(l)%spl, grids(l)%v, sol)
resid(l) = residue(grids(l)%mata, grids(l)%f1d, grids(l)%v1d)
WRITE(*,'(3i5,2(1pe12.3))') l, grids(l)%n, Errdisc(l), resid(l)
END DO
!
! Grid values of direct solutions at the finest levels
ALLOCATE(sol_direct_grid(0:n(1),0:n(2)))
ALLOCATE(sol_anal_grid(0:n(1),0:n(2)))
CALL gridval(grids(1)%spl, grids(1)%x, grids(1)%y, sol_direct_grid, &
& [0,0], grids(1)%v)
sol_anal_grid = sol(grids(1)%x, grids(1)%y)
!--------------------------------------------------------------------------------
! 2. Test restrict and prolong
!
WRITE(*,'(/a)') 'Testing restrict and prolong...'
WRITE(*,'(3a5,3a12)') 'l', 'nx', 'ny', 'rhs', 'sol', 'disc_err'
ALLOCATE(err_restrict(2:levels))
ALLOCATE(err_prolong(2:levels))
ALLOCATE(disc_err_prolong(2:levels))
ALLOCATE(vfine_grid(0:n(1),0:n(2)))
DO l=2,levels
ALLOCATE(fcoarse(SIZE(grids(l)%f,1),SIZE(grids(l)%f,2)))
fcoarse_1d(1:SIZE(grids(l)%f1d)) => fcoarse
ALLOCATE(vfine(SIZE(grids(l-1)%v,1),SIZE(grids(l-1)%v,2)))
vfine_1d(1:SIZE(grids(l-1)%v1d)) => vfine
!
fcoarse = restrict(grids(l)%matp, grids(l-1)%f)
err_restrict(l) = MAXVAL(ABS(fcoarse_1d-grids(l)%f1d))
!
CALL direct_solve(grids(l), fcoarse_1d)
vfine = prolong(grids(l)%matp, fcoarse)
disc_err_prolong(l) = disc_err(grids(l-1)%spl, vfine, sol)
err_prolong(l) = MAXVAL(ABS(vfine_1d-grids(l-1)%v1d))
!
IF(l.EQ.2) THEN ! Grid val on finest grid
CALL gridval(grids(1)%spl, grids(1)%x, grids(1)%y, vfine_grid, &
& [0,0], vfine)
END IF
!
WRITE(*,'(3i5,3(1pe12.3))') l, grids(l)%n, err_restrict(l), err_prolong(l), &
& disc_err_prolong(l)
DEALLOCATE(fcoarse)
DEALLOCATE(vfine)
END DO
!--------------------------------------------------------------------------------
! 9. Epilogue
!
! Creata HDF5 file
!
IF(me.EQ.0) CALL h5file
!
CALL mpi_finalize(ierr)
!--------------------------------------------------------------------------------
CONTAINS
!+++
FUNCTION rhs(x, y)
!
! Return problem RHS
!
DOUBLE PRECISION, INTENT(in) :: x, y
DOUBLE PRECISION :: rhs
rhs = SIN(PI*kx*x)*SIN(PI*ky*y)
END FUNCTION rhs
!+++
FUNCTION sol(x, y)
!
! Return exact problem solution
!
DOUBLE PRECISION, INTENT(in) :: x(:), y(:)
DOUBLE PRECISION :: sol(SIZE(x),SIZE(y))
DOUBLE PRECISION :: c
INTEGER :: j
DO j=1,SIZE(y)
c = SIN(PI*ky*y(j)) / (PI**2*(kx**2+ky**2) + sigma**2)
sol(:,j) = c * SIN(PI*kx*x(:))
END DO
END FUNCTION sol
!+++
SUBROUTINE coefeq(x, y, idt, idw, c)
!
! Weak form = Int( \nabla(w).\nabla(t) + \sigma.t.w) dV)
!
DOUBLE PRECISION, INTENT(in) :: x, y
INTEGER, INTENT(out) :: idt(:,:), idw(:,:)
DOUBLE PRECISION, INTENT(out) :: c(:)
!
c(1) = 1.0d0
idt(1,1) = 1
idt(1,2) = 0
idw(1,1) = 1
idw(1,2) = 0
!
c(2) = 1.0d0
idt(2,1) = 0
idt(2,2) = 1
idw(2,1) = 0
idw(2,2) = 1
!
c(3) = sigma
idt(3,1) = 0
idt(3,2) = 0
idw(3,1) = 0
idw(3,2) = 0
END SUBROUTINE coefeq
!+++
SUBROUTINE h5file
USE futils
CHARACTER(len=128) :: file='test_transf2d.h5'
INTEGER :: fid
INTEGER :: l
CHARACTER(len=64) :: dsname
CALL creatf(file, fid, real_prec='d')
CALL attach(fid, '/', 'NX', n(1))
CALL attach(fid, '/', 'NY', n(2))
CALL attach(fid, '/', 'NIDBAS1', nidbas(1))
CALL attach(fid, '/', 'NIDBAS2', nidbas(2))
CALL attach(fid, '/', 'KX', kx)
CALL attach(fid, '/', 'KY', ky)
CALL attach(fid, '/', 'SIGMA', sigma)
CALL attach(fid, '/', 'ALPHA1', alpha(1))
CALL attach(fid, '/', 'ALPHA2', alpha(2))
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', grids(l)%mata)
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)
CALL putarr(fid, TRIM(dsname)//'/f1d', grids(l)%f1d)
CALL putarr(fid, TRIM(dsname)//'/v1d', grids(l)%v1d)
END DO
!
! Solutions at finest grid
!
CALL creatg(fid, '/solutions')
CALL putarr(fid, '/solutions/xg', grids(1)%x)
CALL putarr(fid, '/solutions/yg', grids(1)%y)
CALL putarr(fid, '/solutions/direct', sol_direct_grid)
CALL putarr(fid, '/solutions/anal', sol_anal_grid)
CALL putarr(fid, '/solutions/vfine', vfine_grid)
!
! Some errors
!
CALL creatg(fid, '/errors')
CALL putarr(fid, '/errors/errdisc', errdisc)
CALL putarr(fid, '/errors/resid', resid)
CALL putarr(fid, '/errors/restrict', err_restrict)
CALL putarr(fid, '/errors/prolong', err_prolong)
CALL putarr(fid, '/errors/disc_err_prolong', disc_err_prolong)
!
CALL closef(fid)
END SUBROUTINE h5file
!+++
END PROGRAM