Homec4science
Diffusion SPClibs (feature/fft-split) Browse

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

Last Change
Raw File

test_stencilg.f90
View Options

!>
!> @file test_stencilg.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>
!>
MODULE mod
USE iso_fortran_env, ONLY : rkind => real64
IMPLICIT NONE
!
LOGICAL, PARAMETER :: nldebug=.FALSE.
REAL(rkind), PARAMETER :: PI=4.d0*ATAN(1.0d0)
CONTAINS
END MODULE mod
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
PROGRAM main
USE mpi
USE pputils2, ONLY : dist1d
USE gvector, ONLY : gvector_2d, ASSIGNMENT(=), OPERATOR(-)
USE parmg, ONLY : exchange, norm_vec
USE stencil, ONLY : stencil_2d, init, laplacian, putmat, OPERATOR(*)
USE mod
IMPLICIT NONE
!
INTEGER, PARAMETER :: ndims=2
!
INTEGER :: me, neighs(4), npes, ierr
INTEGER, DIMENSION(ndims) :: dims=[0,0]
INTEGER, DIMENSION(ndims) :: coords, comm1d
LOGICAL, DIMENSION(ndims) :: periods=[.FALSE.,.FALSE.]
LOGICAL :: reorder =.FALSE.
INTEGER :: comm_cart
!
INTEGER :: nx=4, ny=4 ! Number of intervals
INTEGER, DIMENSION(ndims) :: e, s, g, npt_glob, npt_loc
!
REAL(rkind), ALLOCATABLE :: xgrid(:), ygrid(:)
TYPE(gvector_2d) :: arr, fexact
TYPE(gvector_2d) :: barr1, barr2, barr3
REAL(rkind) :: dx, dy
REAL(rkind) :: err
INTEGER, DIMENSION(5,2) :: id ! 5-point stencil
INTEGER :: npoints
TYPE(stencil_2d) :: mat
INTEGER :: i, j
!
NAMELIST /in/ nx, ny
!================================================================================
! 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)
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)
!
! Read problem inputs
IF(me.EQ.0) THEN
READ(*,in)
WRITE(*,in)
END IF
CALL mpi_bcast(nx, 1, MPI_INTEGER, 0, comm_cart, ierr)
CALL mpi_bcast(ny, 1, MPI_INTEGER, 0, comm_cart, ierr)
!================================================================================
! 2.0 2d Grid construction
!
! Partition 2D grid
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
WRITE(*,'(a,i3.3,a,2(3i4,", "))') 'PE', me, ' coords, s,e:', &
& (coords(i),s(i),e(i),i=1,ndims)
!
! Global mesh
dx = 1.0d0/REAL(nx)
dy = 1.0d0/REAL(ny)
ALLOCATE(xgrid(0:nx))
ALLOCATE(ygrid(0:ny))
xgrid = [ (i*dx, i=0,nx) ]
ygrid = [ (i*dy, i=0,ny) ]
!================================================================================
! 3.0 FD Laplacian
!
id=RESHAPE([ 0, -1, 0, 1, 0, &
0, 0,-1, 0, 1], &
[5,2])
npoints = 5
CALL init(s, e, id, .FALSE., mat, comm_cart)
!
CALL laplacian(dx, dy, mat)
!================================================================================
! 4.0 Check matrice-vector product
!
! Local arrays with ghost cells
g = [1,1]
arr = gvector_2d(s, e, g)
barr1 = gvector_2d(s, e, g)
barr2 = gvector_2d(s, e, g)
barr3 = gvector_2d(s, e, g)
fexact = gvector_2d(s, e, g)
!
! Constant vector => Laplacian = 0
arr = 1.0d0
CALL exchange(comm_cart, arr)
barr1 = mat*arr
IF(s(1).EQ.0) barr1%val(0,:) = 0.0 ! discard boundary values
IF(e(1).EQ.nx) barr1%val(nx,:) = 0.0
IF(s(2).EQ.0) barr1%val(:,0) = 0.0
IF(e(2).EQ.ny) barr1%val(:,ny) = 0.0
err = norm_vec(barr1, comm_cart, root=0)
IF(me.EQ.0) THEN
WRITE(*,'(/a,1pe12.3)') 'Constant vector: ||B1|| =', err
END IF
!
! Bilinear vector => Laplacian = 0
DO j=s(2),e(2)
DO i=s(1),e(1)
arr%val(i,j) = xgrid(i)*ygrid(j)
END DO
END DO
CALL exchange(comm_cart, arr)
barr2 = mat*arr
IF(s(1).EQ.0) barr2%val(0,:) = 0.0 ! discard boundary values
IF(e(1).EQ.nx) barr2%val(nx,:) = 0.0
IF(s(2).EQ.0) barr2%val(:,0) = 0.0
IF(e(2).EQ.ny) barr2%val(:,ny) = 0.0
err = norm_vec(barr2, comm_cart, root=0)
IF(me.EQ.0) THEN
WRITE(*,'(a,1pe12.3)') 'Bilinear vector: ||B2|| =', err
END IF
!
! Biquadratic vector => Laplacian = fexact
DO j=s(2),e(2)
DO i=s(1),e(1)
arr%val(i,j) = (xgrid(i)*ygrid(j))**2/4.0d0
fexact%val(i,j) = (xgrid(i)**2 + ygrid(j)**2)/2.0d0
END DO
END DO
CALL exchange(comm_cart, arr)
CALL exchange(comm_cart, fexact)
barr3 = mat*arr - fexact
IF(s(1).EQ.0) barr3%val(0,:) = 0.0 ! discard boundary values
IF(e(1).EQ.nx) barr3%val(nx,:) = 0.0
IF(s(2).EQ.0) barr3%val(:,0) = 0.0
IF(e(2).EQ.ny) barr3%val(:,ny) = 0.0
err = norm_vec(barr3, comm_cart, root=0)
IF(me.EQ.0) THEN
WRITE(*,'(a,1pe12.3)') 'Biquadratic vector: ||B3|| =', err
END IF
!================================================================================
! 9.0 Epilogue
CALL h5file
CALL MPI_FINALIZE(ierr)
!
CONTAINS
SUBROUTINE h5file
USE futils
CHARACTER(len=128) :: file='test_stencilg.h5'
INTEGER :: fid
CALL creatf(file, fid, real_prec='d', mpicomm=comm_cart)
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, '/arr', arr%val,(/1,2/), garea=g)
CALL putarrnd(fid, '/barr1', barr1%val,(/1,2/), garea=g)
CALL putarrnd(fid, '/barr2', barr2%val,(/1,2/), garea=g)
CALL putarrnd(fid, '/barr3', barr3%val,(/1,2/), garea=g)
CALL putmat(fid, '/MAT', mat)
CALL closef(fid)
END SUBROUTINE h5file
!
END PROGRAM main
c4science · Help