Homec4science
Diffusion SPClibs (feature/refactoring) Browse

spclibs/multigrid/src/poisson_fd.f90
03540012fedefeature/refactoring

Last Change
Raw File

poisson_fd.f90
View Options

!>
!> @file poisson_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>
!>
PROGRAM main
!
! Solving the following 2d PDE using finite differences:
!
! d2/dx2 f(x,y) + tau*d2/dxdy f(x,y) + d2/dy2 f(x,y) = w(x,y), x,y in [0:Lx][0:Ly]
! w(x,y) = -4*pi^2 *[(kx^2/Lx^2+ky^2/Ly^2)*cos(2*kx*pi*x/Lx)*sin(2*ky*pi*y/Ly)
! -tau*(kx*ky)/(Lx*Ly)*sin(2*kx*pi*x/Lx)*cos(2*ky*pi*y/Lx)]
!
! West, East boundaries: Neumann
! South, North boundaries: Dirichlet
!
! Analytic solution : f(x,y) = cos(2*kx*pi*x/Lx)*sin(2*kx*pi*y/Ly)
!
USE multigrid
USE fdmat_mod
IMPLICIT NONE
INCLUDE 'mpif.h'
!
INTEGER, PARAMETER :: nnumx=32
!
INTEGER :: ierr, np, me
DOUBLE PRECISION, PARAMETER :: PI=4.0d0*ATAN(1.0d0)
DOUBLE PRECISION :: Lx, Ly, icrosst, beta, miome
INTEGER :: n, nx,ny,nz,kx,ky
CHARACTER(len=4) :: prb
INTEGER :: nits
DOUBLE PRECISION :: atol, rtol
LOGICAL :: nldirect, nldebug
!
TYPE(mg_info) :: info ! info for MG
INTEGER :: levels, nnu, mu, nu0
!
INTEGER :: inu, nu1(nnumx), nu2(nnumx), niter(nnumx)
DOUBLE PRECISION :: titer(nnumx)
!
LOGICAL :: nlfixed
DOUBLE PRECISION :: omega
CHARACTER(len=4) :: mat_type, relax
!
DOUBLE PRECISION :: dx, dy
INTEGER :: ix, iy, l, its
DOUBLE PRECISION, ALLOCATABLE :: x(:), y(:), dense(:)
TYPE(grid2d), ALLOCATABLE :: grids(:)
!
DOUBLE PRECISION, ALLOCATABLE, TARGET :: sol_ana2d(:,:), sol_direct2d(:,:)
DOUBLE PRECISION, POINTER :: sol_ana(:), sol_direct(:)
DOUBLE PRECISION :: err_direct, resid_direct
DOUBLE PRECISION :: norma, normb
DOUBLE PRECISION, ALLOCATABLE :: rresid(:), resid(:), errdisc(:)
DOUBLE PRECISION :: t0, tsetup, tmat(2), tdirect, tbsolve
DOUBLE PRECISION, EXTERNAL :: mem
!
NAMELIST /parameters/ prb, mat_type,nx, ny, nz, kx, ky, Lx, Ly, icrosst, beta, &
& miome, nldebug, nlfixed, levels, nnu, nu1, nu2, mu, nu0, &
& relax,omega, nldirect, nits, atol, rtol
!--------------------------------------------------------------------------------
! 1.0 Prologue
!
CALL MPI_INIT(ierr)
CALL MPI_COMM_RANK(MPI_COMM_WORLD,me,ierr)
CALL MPI_COMM_SIZE(MPI_COMM_WORLD,np,ierr)
!
! Default inputs
!
nx=32
ny=32
nz=1
kx=1
ky=1
icrosst=1.0d0
Lx = 1.0D0
Ly = 1.0D0
beta = 0d0
miome = 200d0
nldebug = .FALSE.
prb = 'dddd'
mat_type = 'cds'
nldirect = .TRUE.
!
nlfixed = .FALSE.
levels = 2
nnu = 1
nu1 = 1
nu2 = 1
mu = 1
nu0 = 1
nits = 10
atol = 1.e-8; rtol = 1.e-8
relax = 'jac'
omega = 0.6667
!
IF(me==0) THEN
READ(*,parameters)
WRITE(*,parameters)
END IF
!
! Send input parameters to other processors
!
CALL MPI_BCAST(nx, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
CALL MPI_BCAST(ny, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
CALL MPI_BCAST(nz, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
CALL MPI_BCAST(kx, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
CALL MPI_BCAST(ky, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
CALL MPI_BCAST(icrosst, 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
CALL MPI_BCAST(Lx,1,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD, ierr)
CALL MPI_BCAST(Ly,1,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD, ierr)
CALL MPI_BCAST(beta,1,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD, ierr)
CALL MPI_BCAST(miome, 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
!
CALL mpi_bcast(nldebug, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr)
CALL mpi_bcast(nlfixed, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr)
CALL mpi_bcast(nldirect, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr)
CALL mpi_bcast(nnu, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
CALL mpi_bcast(nu1, nnumx, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
CALL mpi_bcast(nu2, nnumx, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
CALL mpi_bcast(nu0, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
CALL mpi_bcast(mu, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
CALL mpi_bcast(nits, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
CALL mpi_bcast(relax, 4, MPI_CHARACTER, 0, MPI_COMM_WORLD, ierr)
CALL mpi_bcast(mat_type, 4, MPI_CHARACTER, 0, MPI_COMM_WORLD, ierr)
CALL mpi_bcast(omega, 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
CALL mpi_bcast(atol, 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
CALL mpi_bcast(rtol, 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
CALL mpi_bcast(prb, LEN(prb), MPI_CHARACTER, 0, MPI_COMM_WORLD, ierr)
!
IF(nnu.GT.nnumx) THEN
IF(me.EQ.0) THEN
PRINT*, 'Value of nnu larger than', nnumx
END IF
CALL mpi_finalize(ierr)
STOP
END IF
!
! Adjust number of levels and fill mg info.
!
levels = MIN(levels, get_lmax(nx), get_lmax(ny))
info%nu1 = nu1(1)
info%nu2 = nu2(1)
info%mu = mu
info%nu0 = nu0
info%levels = levels
info%relax = relax
info%omega = omega
!--------------------------------------------------------------------------------
! 2.0 Setup grids
!
! Grid on the finest level
!
dx = lx/REAL(nx,8)
dy = ly/REAL(ny,8)
ALLOCATE(x(0:nx), y(0:ny))
DO ix=0,nx
x(ix) = ix*dx
END DO
DO iy=0,ny
y(iy) = iy*dy
END DO
WRITE(*,'(a,3(1pe12.3))') 'dx, dy, dx/dy =', dx, dy, dx/dy
!
ALLOCATE(dense(0:nx))
dense = fdense(x)
!
! Set up grids
!
t0 = mpi_wtime()
ALLOCATE(grids(levels))
CALL create_grid_fd(x, y, grids, info, mat_type=mat_type, debug=nldebug)
WRITE(*,'(5a6)') 'l', 'nx', 'ny', 'rx', 'ry'
WRITE(*,'(5i6)') (l, grids(l)%n, grids(l)%rank, l=1,levels)
IF(nldebug) THEN
CALL printmat('** Prolongation matrix in 1st dim.**', grids(2)%transf(1))
CALL printmat('** Prolongation matrix in 2nd dim.**', grids(2)%transf(2))
END IF
!
! Set BC on grid transfer matrices
!
IF(prb.EQ.'dddd') CALL ibc_transf(grids,1,3) ! Direction X
CALL ibc_transf(grids,2,3) ! Direction Y
tsetup = mpi_wtime()-t0
!--------------------------------------------------------------------------------
! 3.0 Problem discretization
!
! Construct FD matrix and impose BC on all grids
!
t0=mpi_wtime()
DO l=1,levels
CALL fdmat(grids(l), fdense, icrosst)
IF(mat_type.EQ.'csr') CALL to_mat(grids(l)%mata)
END DO
tmat(1) = mpi_wtime()-t0
!
t0=mpi_wtime()
DO l=1,levels
CALL ibc_fdmat(grids(l), prb)
END DO
tmat(2) = mpi_wtime()-t0
!
! Set RHS and impose BC on the fiest grid
!
grids(1)%f(:,:) = frhs(x,y)
!
IF(prb.EQ.'dddd') THEN
grids(1)%f(0,:) = 0.0d0 ! Dirichlet on west and east
grids(1)%f(nx,:) = 0.0d0
ELSE IF(prb.EQ.'nndd') THEN ! Neumann on west and east
grids(1)%f(0,:) = 0.5d0*grids(1)%f(0,:)
grids(1)%f(nx,:) = 0.5d0*grids(1)%f(nx,:)
END IF
grids(1)%f(:,0) = 0.0d0 ! Dirichlet on south and north
grids(1)%f(:,ny) = 0.0d0
!
!--------------------------------------------------------------------------------
! 4.0 Analytical solutions and RHS at the finest grid (l=1)
!
n = (nx+1)*(ny+1) ! Number of unknowns
ALLOCATE(sol_ana2d(0:nx,0:ny))
sol_ana(1:n) => sol_ana2d
sol_ana2d(:,:) = fsol(x,y)
!--------------------------------------------------------------------------------
! 5.0 Direct solution at the finest grid (l=1)
!
IF(nldirect) THEN
WRITE(*,'(/a)') 'Direct solution for the finest grid problem ...'
ALLOCATE(sol_direct2d(0:nx,0:ny))
sol_direct(1:n) => sol_direct2d
!
t0 = mpi_wtime()
sol_direct = grids(1)%f1d
CALL direct_solve(grids(1), sol_direct, debug=nldebug)
tdirect = mpi_wtime()-t0
!
t0 = mpi_wtime()
sol_direct = grids(1)%f1d
CALL direct_solve(grids(1), sol_direct, debug=nldebug)
tbsolve = mpi_wtime()-t0
!
! Max norm and residual
!
err_direct = MAXVAL(ABS(sol_direct-sol_ana))
resid_direct = residue(grids(1), grids(1)%f1d, sol_direct, 'inf')
WRITE(*,'(a,2(1pe12.3))') 'Max norm of error and residual norm', &
& err_direct, resid_direct
END IF
!--------------------------------------------------------------------------------
! 5.0 Iterative solution using MG V-cycle
!
WRITE(*,'(/a)') 'Multigrid MG V-cycles ...'
ALLOCATE(errdisc(0:nits))
ALLOCATE(resid(0:nits))
ALLOCATE(rresid(0:nits))
!
! Norm of A and b
!
IF(mat_type.EQ.'csr') THEN
norma = matnorm(grids(1)%mata, 'inf')
ELSE
norma = matnorm(grids(1)%mata_cds, 'inf')
END IF
normb = MAXVAL(ABS(grids(1)%f1d))
WRITE(*,'(a,2(1pe12.3))') 'Norm A and RHS', norma, normb
!
! Initial guess
!
DO inu=1,nnu
info%nu1 = nu1(inu)
info%nu2 = nu2(inu)
WRITE(*,'(/2(a5,i3,2x))') 'nu1 =', nu1(inu), 'nu2 =', nu2(inu)
IF(nlfixed .AND. nldirect) THEN
grids(1)%v = sol_direct2d
ELSE
grids(1)%v = 0.0d0
END IF
!
errdisc(0) = MAXVAL(ABS(grids(1)%v1d-sol_ana))
resid(0) = residue(grids(1), grids(1)%f1d, grids(1)%v1d, 'inf')
rresid(0) = resid(0) / ( norma*MAXVAL(ABS(grids(1)%v1d)) + normb )
WRITE(*,'(a4,3(a12,a8))') 'its', 'residue', 'ratio', 'disc. err', 'ratio', &
& 'rel. resid', 'ratio'
WRITE(*,'(i4,3(1pe12.3,8x))') 0, resid(0), errdisc(0), rresid(0)
!
! Iterations
!
t0 = mpi_wtime()
DO its=1,nits
CALL mg(grids, info, 1)
errdisc(its) = MAXVAL(ABS(grids(1)%v1d-sol_ana))
resid(its) = residue(grids(1), grids(1)%f1d, grids(1)%v1d, 'inf')
rresid(its) = resid(its) / ( norma*MAXVAL(ABS(grids(1)%v1d)) + normb )
WRITE(*,'((i4,3(1pe12.3,0pf8.2)))') its, &
& resid(its), resid(its)/resid(its-1), &
& errdisc(its), errdisc(its)/errdisc(its-1), &
& rresid(its), rresid(its)/rresid(its-1)
IF(resid(its) .LE. atol .or. rresid(its) .le. rtol) EXIT
END DO
niter(inu) = MIN(nits,its)
titer(inu) = mpi_wtime() - t0
END DO
!--------------------------------------------------------------------------------
! 9.0 Epilogue
!
! Display timing
!
WRITE(*,'(/a)') 'Timing ...'
WRITE(*,'(a,1pe12.3,i5)') 'Setup time (s) ', tsetup
WRITE(*,'(a,2(1pe12.3))') 'Matrix construction time(s)', tmat
WRITE(*,'(a,2(1pe12.3))') 'Direct and bsolve time (s) ', tdirect, tbsolve
WRITE(*,'(/3a6,a15)') 'nu1', 'nu2', 'niter', 'Iter time(s)'
DO inu=1,nnu
WRITE(*,'(3i6,3x,1pe12.3)') nu1(inu), nu2(inu), niter(inu), titer(inu)
END DO
!
WRITE(*,'(/a,f12.3)') 'Mem used so far (MB)', mem()
!
! Creata HDF5 file
!
IF(me.EQ.0) CALL h5file
!
! Clean up
!
DEALLOCATE(x)
DEALLOCATE(y)
DEALLOCATE(dense)
DEALLOCATE(grids)
DEALLOCATE(sol_ana2d)
IF(nldirect) DEALLOCATE(sol_direct2d)
DEALLOCATE(errdisc)
DEALLOCATE(resid)
DEALLOCATE(rresid)
!
CALL MPI_FINALIZE(ierr)
!--------------------------------------------------------------------------------
CONTAINS
!+++
FUNCTION fdense(x)
!
! Return density
!
DOUBLE PRECISION, INTENT(in) :: x(:)
DOUBLE PRECISION :: fdense(SIZE(x))
fdense(:) = 0.5d0*beta*miome*EXP( -(x(:)-Lx/3d0)**2d0/(Lx/2d0)**2d0 );
END FUNCTION fdense
!+++
FUNCTION frhs(x,y)
!
! Return RHS
!
DOUBLE PRECISION, INTENT(in) :: x(:), y(:)
DOUBLE PRECISION :: frhs(SIZE(x),SIZE(y))
DOUBLE PRECISION :: c, s, d(SIZE(x))
DOUBLE PRECISION :: 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
!!$ frhs = -frhs
END FUNCTION frhs
!+++
FUNCTION fsol(x,y)
!
! Return analytical solution
!
DOUBLE PRECISION, INTENT(in) :: x(:), y(:)
DOUBLE PRECISION :: fsol(SIZE(x),SIZE(y))
DOUBLE PRECISION :: c
INTEGER :: j
IF(prb.EQ.'dddd') THEN
DO j=1,SIZE(y)
c = SIN(2.0d0*pi*ky*y(j)/Ly)
fsol(:,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)
fsol(:,j) = COS(2.0d0*pi*kx*x(:)/Lx)*c
END DO
END IF
END FUNCTION fsol
!+++
SUBROUTINE h5file
USE futils
CHARACTER(len=128) :: file='poisson_mg.h5'
INTEGER :: fid
INTEGER :: l
CHARACTER(len=64) :: dsname
CALL creatf(file, fid, real_prec='d')
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, '/', 'BETA', beta)
CALL attach(fid, '/', 'OMEGA', omega)
CALL attach(fid, '/', 'RELAX', relax)
CALL attach(fid, '/', 'MAT_TYPE', mat_type)
CALL attach(fid, '/', 'NITS', nits)
CALL attach(fid, '/', 'LEVELS', levels)
CALL attach(fid, '/', 'NNU', nnu)
CALL attach(fid, '/', 'NU0', nu0)
CALL attach(fid, '/', 'MU', mu)
!
CALL putarr(fid, '/nu1', nu1(1:nnu))
CALL putarr(fid, '/nu2', nu2(1:nnu))
CALL putarr(fid, '/dense', dense)
CALL creatg(fid, '/mglevels')
DO l=1,levels
WRITE(dsname,'("/mglevels/level.",i2.2)') l
CALL creatg(fid, TRIM(dsname))
IF(mat_type.EQ.'csr') THEN
CALL putmat(fid, TRIM(dsname)//'/mata', grids(l)%mata)
ELSE
CALL putmat(fid, TRIM(dsname)//'/mata', grids(l)%mata_cds)
END IF
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/calc', grids(1)%v)
CALL putarr(fid, '/solutions/anal', sol_ana2d)
IF(nldirect) CALL putarr(fid, '/solutions/direct', sol_direct2d)
!
nits=niter(nnu)
CALL creatg(fid, '/Iterations')
CALL putarr(fid, '/Iterations/residues', resid(0:nits))
CALL putarr(fid, '/Iterations/disc_errors', errdisc(0:nits))
!
CALL closef(fid)
END SUBROUTINE h5file
!+++
END PROGRAM main
c4science · Help