Homec4science
Diffusion SPClibs (feature/fft-split) Browse

spclibs/multigrid/src/ppoisson_fd.f90
03540012fedefeature/fft-split

Last Change
Raw File

ppoisson_fd.f90
View Options

!>
!> @file ppoisson_fd.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>
!>
!
! Test 2D parallel multigrid V-cycle
!
MODULE mod
USE iso_fortran_env, ONLY : rkind => real64
IMPLICIT NONE
!
REAL(rkind), PARAMETER :: PI=4.d0*ATAN(1.0d0)
CONTAINS
END MODULE mod
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
PROGRAM main
USE mpi
USE fdmat_mod, ONLY : fdmat, ibc_fdmat, ibc_rhs
USE pputils2, ONLY : dist1d, timera, hostlist
USE gvector, ONLY : gvector_2d, norm2, ASSIGNMENT(=), OPERATOR(-)
USE parmg, ONLY : grid2_type, mg_info, create_grid, mg, exchange, &
& get_resids, disp, norm_vec, norm_mat
USE stencil, ONLY : stencil_2d, putmat
USE mod
IMPLICIT NONE
!
INTEGER, PARAMETER :: ndims=2
!
INTEGER :: me, neighs(4), npes, ierr
INTEGER, DIMENSION(ndims) :: coords, comm1d
LOGICAL, DIMENSION(ndims) :: periods=[.FALSE.,.FALSE.]
LOGICAL :: reorder =.FALSE.
INTEGER :: comm_cart
INTEGER, DIMENSION(ndims) :: e0, s0, e, s, npt_glob, npt_loc
!
REAL(rkind), ALLOCATABLE :: xgrid(:), ygrid(:)
INTEGER, ALLOCATABLE :: id(:,:)
REAL(rkind) :: dx, dy
INTEGER :: npoints ! Number of points in FD stencil
REAL(rkind) :: t_mg, t_mg0, t_mg_min, t_mg_max
!
TYPE(gvector_2d) :: v_exact, resids, errs
REAL(rkind) :: norma, normb, normv
REAL(rkind), ALLOCATABLE :: resid_it(:), err_it(:), rresid(:)
REAL(rkind) :: ratio_err, ratio_resid, ratio_rresid
INTEGER, DIMENSION(ndims) :: g, npt_loc_min
INTEGER :: l, i, it
CHARACTER(len=64) :: str
!
TYPE(grid2_type), ALLOCATABLE :: grids(:)
TYPE(mg_info) :: info ! info for MG
!
! Input quantities
!
LOGICAL :: nldebug=.FALSE.
CHARACTER(len=64) :: filein = 'ppoisson_fd.in'
INTEGER :: dims(2)=[0,0]
CHARACTER(len=4) :: prb='dddd'
CHARACTER(len=4) :: relax='jac'
INTEGER :: nx=4, ny=4 ! Number of intervals
INTEGER :: kx=1, ky=1
REAL(rkind) :: Lx=1.0, Ly=1.0
REAL(rkind) :: icrosst=1.0, beta=0.0, miome=200.0
REAL(rkind) :: omega=1.0d0
INTEGER :: nits=100, direct_solve_nits=5
INTEGER :: levels=2, nu1=3, nu2=3, mu=1, nu0=1
REAL(rkind) :: rtol=1.e-8, atol=1.e-8, errtol=1.e-3
!
NAMELIST /in/ nldebug, dims, prb, nx, ny, kx, ky, Lx, Ly, icrosst, beta, &
& miome, omega, nits, levels, relax, nu1, nu2, mu, nu0, &
& direct_solve_nits, rtol, atol, errtol
!================================================================================
! 1.0 Prologue
!
! 2D process grid
!
CALL mpi_init(ierr)
CALL mpi_comm_rank(MPI_COMM_WORLD, me, ierr)
CALL mpi_comm_size(MPI_COMM_WORLD, npes, ierr)
!
! Read input filename from commmand line argument
!
IF( command_argument_count() > 0 ) THEN
CALL get_command_argument(1, filein)
END IF
IF(me.EQ.0) WRITE(*,'(a,a)') 'filein = ', TRIM(filein)
!
! Read problem inputs
!
OPEN(unit=99, file=filein, form='formatted')
READ(99,in)
CLOSE(99)
!
CALL mpi_dims_create(npes, ndims, dims, ierr)
CALL mpi_cart_create(MPI_COMM_WORLD, ndims, dims, periods, reorder, comm_cart,&
& ierr)
CALL mpi_comm_rank(comm_cart, me, ierr)
CALL mpi_cart_coords(comm_cart, me, ndims, coords, ierr)
CALL mpi_cart_shift(comm_cart, 0, 1, neighs(1), neighs(2), ierr)
CALL mpi_cart_shift(comm_cart, 1, 1, neighs(3), neighs(4), ierr)
!
CALL mpi_cart_sub(comm_cart, [.TRUE.,.FALSE.], comm1d(1), ierr)
CALL mpi_cart_sub(comm_cart, [.FALSE.,.TRUE.], comm1d(2), ierr)
!
info%comm = comm_cart
info%nu1 = nu1
info%nu2 = nu2
info%mu = mu
info%nu0 = nu0
info%levels = levels
info%direct_solve_nits = direct_solve_nits
info%relax = relax
info%omega = omega
!
IF(me.EQ.0) THEN
WRITE(*, in)
END IF
IF(nldebug) THEN
CALL hostlist(comm_cart)
END IF
!================================================================================
! 2.0 2d Grid construction
!
! Partition 2D grid
!
CALL timera(0, 'Grid_construction')
npt_glob(1) = nx+1
npt_glob(2) = ny+1
DO i=1,ndims
CALL dist1d(comm1d(i), 0, npt_glob(i), s(i), npt_loc(i))
e(i) = s(i) + npt_loc(i) - 1
END DO
IF(nldebug) THEN
WRITE(*,'(a,i3.3,a,2(3i8,", "))') 'PE', me, ' coords, s,e:', &
& (coords(i),s(i),e(i),i=1,ndims)
END IF
!
! Global mesh
!
dx = Lx/REAL(nx)
dy = Ly/REAL(ny)
ALLOCATE(xgrid(0:nx))
ALLOCATE(ygrid(0:ny))
xgrid = [ (i*dx, i=0,nx) ]
ygrid = [ (i*dy, i=0,ny) ]
!
! Create grid structure
!
ALLOCATE(grids(levels))
npoints = 9 ! Size of FD stencil
ALLOCATE(id(npoints,2))
id=RESHAPE([ 0, -1, 0, 1, -1, 1, -1, 0, 1, &
0, -1, -1, -1, 0, 0, 1, 1, 1], &
[npoints,2])
CALL create_grid(xgrid, ygrid, s, e, id, prb, grids, comm_cart)
!
IF(nldebug) THEN
DO l=1,levels
WRITE(str,'(a,i0)') 'Number of local points at level ', l
CALL disp(TRIM(str), grids(l)%npt_loc, comm_cart)
END DO
END IF
CALL mpi_reduce(grids(levels)%npt_loc, npt_loc_min, 2, MPI_INTEGER, MPI_MIN, &
& 0, comm_cart, ierr)
IF(me.EQ.0) THEN
WRITE(*,'(a,2i4)') 'Minimum local npt at coarsest grid:', npt_loc_min
END IF
!
CALL timera(1, 'Grid_construction')
!================================================================================
! 3.0 FD Operator
!
CALL timera(0, 'FD Operator')
!
DO l=1,levels
CALL fdmat(grids(l), fdense, icrosst, grids(l)%fdmat)
CALL ibc_fdmat(grids(l)%fdmat, prb)
END DO
!
CALL timera(1, 'FD Operator')
!================================================================================
! 4.0 RHS and exact solution at the finest grid (l=1)
!
! Allocate memory
!
CALL timera(0, 'RHS and exact sol')
!
s0 = grids(1)%s0; e0 = grids(1)%e0
s = grids(1)%s; e = grids(1)%e
g = [1,1]
v_exact = gvector_2d(s, e, g) ! Exact solutions
errs = gvector_2d(s, e, g) ! Disc. errors
resids = gvector_2d(s, e, g) ! Residuals
ALLOCATE(resid_it(0:nits))
ALLOCATE(rresid(0:nits))
ALLOCATE(err_it(0:nits))
!
! Set RHS at the finest grid and impose Dirichlet/Neuman BC.
!
grids(1)%f = frhs(xgrid(s(1):e(1)),ygrid(s(2):e(2)))
CALL ibc_rhs(grids(1)%f, s0, e0, prb)
!
! Exact solutions
!
v_exact = fexact(xgrid(s(1):e(1)),ygrid(s(2):e(2)))
!
CALL timera(1, 'RHS and exact sol')
!================================================================================
! 5.0 MG V-cycle iteration loop
!
! Norm of A and b
!
norma = norm_mat(grids(1)%fdmat, comm_cart)
normb = norm_vec(grids(1)%f, comm_cart)
IF(me.EQ.0) THEN
WRITE(*,'(a,2(1pe12.3))') 'Norm A and RHS', norma, normb
END IF
!
grids(1)%v = 0.0d0
resids = get_resids(comm_cart, grids(1)%fdmat, grids(1)%v, grids(1)%f)
errs = grids(1)%v - v_exact
err_it(0) = norm_vec(errs, comm_cart, root=0)
resid_it(0) = norm_vec(resids, comm_cart)
normv = norm_vec(grids(1)%v, comm_cart)
rresid(0) = resid_it(0) / ( norma*normv + normb )
!
IF(me.EQ.0) THEN
WRITE(*,'(a4,3(a12,a8))') 'its', 'residue', 'ratio', 'disc. err', 'ratio', &
& 'rel. resid', 'ratio'
WRITE(*,'(i4,3(1pe12.3,8X))') 0, resid_it(0), err_it(0), rresid(0)
END IF
!
CALL timera(0, 'MG V-cycle loop')
t_mg = 0.0d0
DO it=1,nits
t_mg0 = mpi_wtime()
CALL mg(grids, info, 1)
t_mg = t_mg + (mpi_wtime()-t_mg0)
resids = get_resids(comm_cart, grids(1)%fdmat, grids(1)%v, grids(1)%f)
errs = grids(1)%v - v_exact
err_it(it) = norm_vec(errs, comm_cart)
resid_it(it) = norm_vec(resids, comm_cart)
normv = norm_vec(grids(1)%v, comm_cart)
rresid(it) = resid_it(it) / ( norma*normv + normb )
ratio_err = err_it(it)/err_it(it-1)
ratio_resid = resid_it(it)/resid_it(it-1)
ratio_rresid= rresid(it)/ rresid(it-1)
IF(me.EQ.0) THEN
WRITE(*,'(i4,3(1pe12.3,0pf8.2))') it, &
& resid_it(it), ratio_resid,&
& err_it(it), ratio_err, &
& rresid(it), ratio_rresid
END IF
IF(resid_it(it) .LE. atol .OR. rresid(it) .LE. rtol .OR. &
& ABS(ratio_err-1._rkind).LT.errtol) THEN
nits = it
EXIT
END IF
END DO
!
CALL timera(1, 'MG V-cycle loop')
CALL mpi_reduce(t_mg, t_mg_max, 1, MPI_DOUBLE_PRECISION, MPI_MAX, 0, comm_cart, ierr)
CALL mpi_reduce(t_mg, t_mg_min, 1, MPI_DOUBLE_PRECISION, MPI_MIN, 0, comm_cart, ierr)
IF(me.EQ.0) THEN
WRITE(*,'(a,2(1pe10.3)/)') 'Minmax of MG (only) time (s):', t_mg_min, t_mg_max
END IF
!================================================================================
! 9.0 Epilogue
!
IF(nldebug) THEN
CALL h5file
END IF
!
CALL timera(9, '')
CALL MPI_FINALIZE(ierr)
CONTAINS
!
!+++
FUNCTION fdense(x)
!
! Return density
!
REAL(rkind), INTENT(in) :: x(:)
REAL(rkind) :: fdense(SIZE(x))
fdense(:) = 0.5d0*beta*miome*EXP( -(x(:)-Lx/3d0)**2d0/(Lx/2d0)**2d0 );
END FUNCTION fdense
!+++
FUNCTION fexact(x,y)
!
! Return analytical solution
!
REAL(rkind), INTENT(in) :: x(:), y(:)
REAL(rkind) :: fexact(SIZE(x),SIZE(y))
REAL(rkind) :: c
INTEGER :: j
IF(prb.EQ.'dddd') THEN
DO j=1,SIZE(y)
c = SIN(2.0d0*pi*ky*y(j)/Ly)
fexact(:,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)
fexact(:,j) = COS(2.0d0*pi*kx*x(:)/Lx)*c
END DO
END IF
END FUNCTION fexact
!+++
FUNCTION frhs(x,y)
!
! Return RHS
!
REAL(rkind), INTENT(in) :: x(:), y(:)
REAL(rkind) :: frhs(SIZE(x),SIZE(y))
REAL(rkind) :: c, s, d(SIZE(x))
REAL(rkind) :: 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 outerprod(x, y) RESULT(r)
!
! outer product
!
REAL(rkind), INTENT(in) :: x(:), y(:)
REAL(rkind) :: r(SIZE(x),SIZE(y))
INTEGER :: i, j
DO j=1,SIZE(y)
DO i=1,SIZE(x)
r(i,j) = x(i)*y(j)
END DO
END DO
END FUNCTION outerprod
!+++
SUBROUTINE h5file
!
! Result hdf5 file
!
USE futils
CHARACTER(len=128) :: file='ppoisson_fd.h5'
INTEGER :: fid
CALL creatf(file, fid, real_prec='d', mpicomm=comm_cart)
CALL attach(fid, '/', 'PRB', prb)
CALL attach(fid, '/', 'RELAX', relax)
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, '/', 'ICROSST', icrosst)
CALL attach(fid, '/', 'BETA', beta)
CALL attach(fid, '/', 'MIOME', miome)
CALL attach(fid, '/', 'OMEGA', omega)
CALL attach(fid, '/', 'NITS', nits)
CALL attach(fid, '/', 'DIRECT_SOLVE_NITS', direct_solve_nits)
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 putarr(fid, '/xgrid', xgrid, ionode=0) ! only rank 0 does IO
CALL putarr(fid, '/ygrid', ygrid, ionode=0) ! only rank 0 does IO
!
CALL putarrnd(fid, '/f', grids(1)%f%val, (/1,2/), garea=g)
CALL putarrnd(fid, '/v', v_exact%val, (/1,2/), garea=g)
CALL putarrnd(fid, '/u', grids(1)%v%val, (/1,2/), garea=g)
CALL putarrnd(fid, '/errs', errs%val, (/1,2/), garea=(/1,1/))
CALL putarrnd(fid, '/resids', resids%val,(/1,2/), garea=(/1,1/))
!
CALL putarr(fid, '/resid', resid_it(0:nits), ionode=0)
CALL putarr(fid, '/error', err_it(0:nits), ionode=0)
!
CALL putmat(fid, '/MAT', grids(1)%fdmat)
CALL closef(fid)
END SUBROUTINE h5file
!+++
END PROGRAM main
c4science · Help