Homec4science
Diffusion SPClibs (feature/fft-split) Browse

spclibs/examples/pde3d.f90
03540012fedefeature/fft-split

Last Change
Raw File

pde3d.f90
View Options

!>
!> @file pde3d.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 3d PDE using splines:
!
! -d/dx[x d/dx]f - 1/x[d/dy]^2 f = 4(m+1)x^(m+1)cos(my)cos(z)^n, with f(x=1,y) = 0
! exact solution: f(x,y) = (1-x^2) x^m cos(my)cos(z)^n
!
USE futils
USE fft
USE pde3d_mod
!
IMPLICIT NONE
INTEGER :: nx, ny, nz, nidbas(3), ngauss(3), mbess, npow, nterms
LOGICAL :: nlppform
INTEGER :: i, j, k, kk, ij, dimx, dimy, dimz, nrank, kl, ku
INTEGER :: jder(3), it
DOUBLE PRECISION :: pi, coefx(5)
DOUBLE PRECISION :: dy, dz
DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: xgrid, ygrid, zgrid
DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: fftmass
DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: rhs, sol
DOUBLE COMPLEX, DIMENSION(:,:), ALLOCATABLE :: crhs
!
TYPE(spline2d1d), TARGET :: splxyz
TYPE(spline2d), POINTER :: splxy
TYPE(gbmat) :: mat
!
CHARACTER(len=128) :: file='pde3d.h5'
INTEGER :: fid
DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: arr
DOUBLE PRECISION, DIMENSION(:,:,:), ALLOCATABLE :: bcoef, solcal, solana, errsol
DOUBLE PRECISION :: seconds, mem, dopla
DOUBLE PRECISION :: t0, tmat, tfact, tsolv, tgrid, gflops1, gflops2
INTEGER :: nits=500
!
INTEGER, PARAMETER :: npart=10
DOUBLE PRECISION, DIMENSION(npart) :: xp, yp, zp, fp_calc, fp_anal
!
INTEGER :: kmin, kmax
DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: fftmass_shifted
!
NAMELIST /newrun/ nx, ny, nz, nidbas, ngauss, mbess, npow, nlppform, coefx
!===========================================================================
! 1.0 Prologue
!
! Read in data specific to run
!
nx = 8 ! Number of intervals in x
ny = 8 ! Number of intervals in y
nz = 8 ! Number of intervals in z
nidbas = (/3,3,3/) ! Degree of splines
ngauss = (/4,4, 4/) ! Number of Gauss points/interval
mbess = 2 ! Exponent of differential problem
npow = 2 ! Exponent of differential problem
nterms = 2 ! Number of terms in weak form
nlppform = .TRUE. ! Use PPFORM for gridval or not
coefx(1:5) = (/1.0d0, 0.d0, 0.d0, 0.d0, 1.d0/) ! Mesh point distribution function
!
READ(*,newrun)
WRITE(*,newrun)
!
! Define grid on x (=radial) & y (=poloidal) axis
!
pi = 4.0d0*ATAN(1.0d0)
ALLOCATE(xgrid(0:nx), ygrid(0:ny), zgrid(0:nz))
!
xgrid(0) = 0.0d0; xgrid(nx) = 1.0d0
CALL meshdist(coefx, xgrid, nx)
!
dy = 2.d0*pi/REAL(ny,8) ! Equidistant in y
ygrid = (/ (j*dy, j=0,ny) /)
!
dz = 2.0d0*pi/REAL(nz,8) ! Equidistant in z
zgrid = (/ (k*dz, k=0,nz) /)
!
WRITE(*,'(a/(10f8.3))') 'XGRID', xgrid(0:nx)
WRITE(*,'(a/(10f8.3))') 'YGRID', ygrid(0:ny)
WRITE(*,'(a/(10f8.3))') 'ZGRID', zgrid(0:nz)
!
! Create hdf5 file
!
CALL creatf(file, fid, 'PDE2D Result File', real_prec='d')
CALL attach(fid, '/', 'NX', nx)
CALL attach(fid, '/', 'NY', ny)
CALL attach(fid, '/', 'NZ', nz)
CALL attach(fid, '/', 'NIDBAS1', nidbas(1))
CALL attach(fid, '/', 'NIDBAS2', nidbas(2))
CALL attach(fid, '/', 'NIDBAS3', nidbas(3))
CALL attach(fid, '/', 'NGAUSS1', ngauss(1))
CALL attach(fid, '/', 'NGAUSS2', ngauss(2))
CALL attach(fid, '/', 'NGAUSS2', ngauss(3))
CALL attach(fid, '/', 'MBESS', mbess)
!===========================================================================
! 2.0 Discretize the PDE
!
! Set up spline
!
t0 = seconds()
CALL set_spline(nidbas, ngauss, xgrid, ygrid, zgrid, splxyz, &
& (/.FALSE., .TRUE., .TRUE./), nlppform=nlppform)
splxy => splxyz%sp12
!
WRITE(*,'(/a,/(12(f8.3)))') 'KNOTS in X', splxy%sp1%knots
WRITE(*,'(/a,/(12(f8.3)))') 'KNOTS in Y', splxy%sp2%knots
WRITE(*,'(/a,/(12(f8.3)))') 'KNOTS in Z', splxyz%sp3%knots
!
! 2D FE matrix assembly (in plane x-y)
!
nrank = (nx+nidbas(1))*ny ! Rank of the FE matrix
kl = (nidbas(1)+1)*ny -1 ! Number of sub-diagnonals
ku = kl ! Number of super-diagnonals
WRITE(*,'(a,5i6)') 'nrank, kl, ku', nrank, kl, ku
!
CALL init(kl, ku, nrank, nterms, mat)
CALL dismat(splxy, mat)
!!$ CALL putmat(fid, '/MAT0', mat, 'Assembled GB matrice')
ALLOCATE(arr(nrank))
!
! BC on Matrix
!
CALL ibcmat(mat, ny)
tmat = seconds() - t0
!
! 3D RHS assembly
!
ALLOCATE(rhs(nrank,0:nz-1), sol(nrank,0:nz-1))
CALL disrhs3(mbess, npow, splxyz, rhs)
!
! FFT in z of RHS
!
ALLOCATE(crhs(nrank,0:nz-1))
crhs = rhs
CALL fourrow(crhs, 1)
crhs = crhs/REAL(nz,8)
!
! Apply Mass matrix to crhs
!
kmin =-nz/2
kmax = nz/2-1
CALL init_dft(splxyz%sp3, kmin, kmax)
ALLOCATE(fftmass_shifted(kmin:kmax))
ALLOCATE(fftmass(0:nz-1))
CALL calc_fftmass(splxyz%sp3, fftmass_shifted)
DO k=kmin,kmax
fftmass(MODULO(k+nz,nz)) = fftmass_shifted(k)
END DO
DO k=0,nz-1
crhs(:,k) = crhs(:,k)/fftmass(k)
END DO
!
! Fourier transform back crhs to real space in z
!
CALL fourrow(crhs, -1)
rhs(:,:) = REAL(crhs(:,:),8)
!
! BC on RHS
!
CALL ibcrhs3(rhs, ny)
!
CALL putarr(fid, '/RHS', rhs, 'RHS of linear system')
WRITE(*,'(a,1pe12.3)') 'Mem used so far (MB)', mem()
!===========================================================================
! 3.0 Solve the dicretized system
!
t0 = seconds()
CALL factor(mat)
tfact = seconds() - t0
gflops1 = dopla('DGBTRF',nrank,nrank,kl,ku,0) / tfact / 1.d9
t0 = seconds()
CALL bsolve(mat, rhs, sol)
!
! Backtransform of solution
!
DO k=0,nz-1
sol(1:ny-1,k) = sol(ny,k)
END DO
tsolv = seconds() - t0
gflops2 = dopla('DGBTRS',nrank,1,kl,ku,0) / tsolv / 1.d9
CALL putarr(fid, '/SOL', sol, 'Spline coefficients of solution')
!
! Spline coefficients, taking into account of periodicity in y and z
! Note: in SOL, y was numbered first.
!
dimx = splxy%sp1%dim
dimy = splxy%sp2%dim
dimz = splxyz%sp3%dim
WRITE(*,'(/a,3i6)') 'dimx, dimy, dimz =', dimx, dimy, dimz
ALLOCATE(bcoef(0:dimx-1, 0:dimy-1, 0:dimz-1))
DO j=0,dimy-1
DO i=0,dimx-1
ij = MODULO(j,ny) + i*ny + 1
DO k=0,dimz-1
kk = MODULO(k,nz)
bcoef(i,j,k) = sol(ij,kk)
END DO
END DO
END DO
CALL putarr(fid, '/BCOEF', bcoef, 'Spline coefficients of solution')
!===========================================================================
! 4.0 Check the solution
!
! Check function values computed with various method
!
WRITE(*,'(/a)') 'Check with analytical solutions ...'
CALL RANDOM_NUMBER(xp)
yp=0.0d0
zp=0.0d0
jder = (/0,0,0/)
CALL gridval(splxyz, xp, yp, zp, fp_calc, jder, bcoef)
!!$ WRITE(*,'(4a12)') 'X', 'CALC', 'ANAL', 'ERROR'
!!$ DO i=1,npart
!!$ fp_anal(i) = (1-xp(i)**2) * xp(i)**mbess &
!!$ & * COS(mbess*yp(i)) * COS(zp(i))**npow
!!$ WRITE(*,'(4(1pe12.3))') xp(i), fp_calc(i), fp_anal(i), fp_calc(i)-fp_anal(i)
!!$ END DO
!
ALLOCATE(solcal(0:nx,0:ny,0:nz))
ALLOCATE(solana(0:nx,0:ny,0:nz))
ALLOCATE(errsol(0:nx,0:ny,0:nz))
DO i=0,nx
DO j=0,ny
DO k=0,nz
solana(i,j,k) = (1-xgrid(i)**2) * xgrid(i)**mbess &
& * COS(mbess*ygrid(j)) * COS(zgrid(k))**npow
END DO
END DO
END DO
!
jder = (/0,0,0/)
CALL gridval(splxyz, xgrid, ygrid, zgrid, solcal, jder, bcoef)
t0 = seconds()
CALL gridval(splxyz, xgrid, ygrid, zgrid, solcal, jder)
tgrid = seconds()-t0
errsol = solana - solcal
WRITE(*,'(a,2(1pe12.3))') 'Relative discretization errors', &
& norm2(errsol) / norm2(solana)
!
CALL putarr(fid, '/xgrid', xgrid, 'r')
CALL putarr(fid, '/ygrid', ygrid, '\theta')
CALL putarr(fid, '/zgrid', zgrid, '\phi')
CALL putarr(fid, '/sol', solcal, 'Solutions')
CALL putarr(fid, '/solana', solana,'Exact solutions')
!
! Check derivative d/dx
!
DO i=0,nx
DO j=0,ny
DO k=0,nz
IF( mbess .EQ. 0 ) THEN
solana(i,j,k) = -2.0d0 * xgrid(i) * COS(zgrid(k))**npow
ELSE
solana(i,j,k) = (-(mbess+2)*xgrid(i)**2 + mbess) * &
& xgrid(i)**(mbess-1) * COS(mbess*ygrid(j)) * &
& COS(zgrid(k))**npow
END IF
END DO
END DO
END DO
!
jder = (/1,0,0/)
t0 = seconds()
CALL gridval(splxyz, xgrid, ygrid, zgrid, solcal, jder)
tgrid = tgrid + seconds()-t0
errsol = solana - solcal
CALL putarr(fid, '/derivx', solcal, 'd/dx of solutions')
CALL putarr(fid, '/derivx_exact', solana)
WRITE(*,'(a,2(1pe12.3))') 'Relative error in d/dx', norm2(errsol)/norm2(solana)
!
! Check derivative d/dy
!
DO i=0,nx
DO j=0,ny
DO k=0,nz
solana(i,j,k) = -mbess * (1-xgrid(i)**2) * xgrid(i)**mbess * &
& SIN(mbess*ygrid(j))* COS(zgrid(k))**npow
END DO
END DO
END DO
!
jder = (/0,1,0/)
t0 = seconds()
CALL gridval(splxyz, xgrid, ygrid, zgrid, solcal, jder)
tgrid = tgrid + seconds()-t0
CALL putarr(fid, '/derivy', solcal, 'd/dy of solutions')
CALL putarr(fid, '/derivy_exact', solana)
errsol = solana - solcal
WRITE(*,'(a,2(1pe12.3))') 'Relative error in d/dy', norm2(errsol)/norm2(solana)
!
! Check derivative d/dz
!
DO i=0,nx
DO j=0,ny
DO k=0,nz
solana(i,j,k) = -npow*(1-xgrid(i)**2) * xgrid(i)**mbess &
& * COS(mbess*ygrid(j)) * COS(zgrid(k))**(npow-1) &
& * SIN(zgrid(k))
END DO
END DO
END DO
!
jder = (/0,0,1/)
t0 = seconds()
IF(nlppform) THEN
CALL gridval(splxyz, xgrid, ygrid, zgrid, solcal, jder)
ELSE
CALL gridval(splxyz, xgrid, ygrid, zgrid, solcal, jder, bcoef)
END IF
tgrid = tgrid + seconds()-t0
CALL putarr(fid, '/derivz', solcal, 'd/dz of solutions')
CALL putarr(fid, '/derivz_exact', solana)
errsol = solana - solcal
WRITE(*,'(a,2(1pe12.3))') 'Relative error in d/dz', norm2(errsol)/norm2(solana)
!===========================================================================
! 9.0 Epilogue
!
WRITE(*,'(/a)') '---'
WRITE(*,'(a,1pe12.3)') 'Matrice construction time (s) ', tmat
WRITE(*,'(a,1pe12.3)') 'Matrice factorisation time (s)', tfact
WRITE(*,'(a,1pe12.3)') 'Backsolve time (s) ', tsolv
WRITE(*,'(a,1pe12.3)') 'gridval time (s) ', tgrid
WRITE(*,'(a,2f10.3)') 'Factor/solve Gflop/s', gflops1, gflops2
WRITE(*,'(a,1pe12.3)') 'Mem used so far (MB)', mem()
!
DEALLOCATE(xgrid, ygrid, zgrid)
DEALLOCATE(rhs, sol)
DEALLOCATE(crhs)
DEALLOCATE(fftmass)
DEALLOCATE(fftmass_shifted)
DEALLOCATE(bcoef)
DEALLOCATE(solcal, solana, errsol)
DEALLOCATE(arr)
CALL destroy_sp(splxyz)
CALL destroy(mat)
!
CALL closef(fid)
!===========================================================================
!
CONTAINS
FUNCTION norm2(x)
!
! Compute the 2-norm of array x
!
IMPLICIT NONE
DOUBLE PRECISION :: norm2
DOUBLE PRECISION, INTENT(in) :: x(:,:,:)
DOUBLE PRECISION :: sum2
INTEGER :: i, j
!
sum2 = 0.0d0
DO i=1,SIZE(x,1)
DO j=1,SIZE(x,2)
DO k=1,SIZE(x,3)
sum2 = sum2 + x(i,j,k)**2
END DO
END DO
END DO
norm2 = SQRT(sum2)
END FUNCTION norm2
SUBROUTINE prnmat(label, mat)
CHARACTER(len=*) :: label
DOUBLE PRECISION, DIMENSION(:,:), INTENT(in) :: mat
INTEGER :: i
WRITE(*,'(/a)') TRIM(label)
DO i=1,SIZE(mat,1)
WRITE(*,'(10(1pe12.3))') mat(i,:)
END DO
END SUBROUTINE prnmat
END PROGRAM main
!
!+++
c4science · Help