Homec4science
Diffusion SPClibs (master) Browse

spclibs/examples/tcdsmat.f90
253e659ccb68master

Last Change
Raw File

tcdsmat.f90
View Options

!>
!> @file tcdsmat.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/>.
!>
!> @author
!> (in alphabetical order)
!> @author Trach-Minh Tran <trach-minh.tran@epfl.ch>
!>
PROGRAM main
!
! Solving the following 2d PDE using splines and iterative method
!
! -d/dx[x d/dx]f - 1/x[d/dy]^2 f = 4(m+1)x^(m+1)cos(my), with f(x=1,y) = 0
! exact solution: f(x,y) = (1-x^2) x^m cos(my)
!
USE tcdsmat_mod
USE bsplines
USE cds
USE futils
!
IMPLICIT NONE
INCLUDE 'mpif.h'
!
INTEGER :: nints(2), nidbas(2), ngauss(2), mbess, nterms
DOUBLE PRECISION :: pi, coefx(5), coefy(5)
DOUBLE PRECISION, ALLOCATABLE :: xgrid(:), ygrid(:)
DOUBLE PRECISION, ALLOCATABLE :: rhs(:), sol(:)
INTEGER :: mrank, bw0
TYPE(spline2d) :: splxy
TYPE(cds_mat) :: mat
LOGICAL :: readmat, verbose
CHARACTER(len=128) :: file='tcdsmat.h5'
CHARACTER(len=128) :: filein
INTEGER :: fid, fidin
DOUBLE PRECISION :: mem, seconds
DOUBLE PRECISION :: t0, tmat, tbal, tsolv, tgrid, tmumps(2)
INTEGER :: nitmx=100, niter, nssor
DOUBLE PRECISION :: rtolmx=1.0d-6, omega=0.0d0, resid
!
INTEGER :: i, j, ij, dimx, dimy
DOUBLE PRECISION, ALLOCATABLE :: bcoef(:,:)
DOUBLE PRECISION, ALLOCATABLE :: solcal(:,:), solana(:,:), errsol(:,:)
INTEGER, ALLOCATABLE :: dists(:)
!
INTEGER :: ierr, me
TYPE(mumps_mat) :: mat_mumps
!
NAMELIST /newrun/ nints, nidbas, ngauss, mbess, coefx, coefy, &
& nitmx, rtolmx, omega, nssor, readmat, verbose, filein
!===========================================================================
! 1.0 Prologue
!
CALL mpi_init(ierr)
CALL mpi_comm_rank(MPI_COMM_WORLD, me, ierr)
!
! Read in data specific to run
!
WRITE(*,'(/a)') 'Prologue ...'
readmat = .FALSE. ! Read matrix and rhs from file
filein = 'mat.h5'
nints = (/8,8/) ! Number of intervals in x, y
nidbas = (/3,3/) ! Degree of splines
ngauss = (/4,4/) ! Number of Gauss points/interval
mbess = 2 ! Exponent of differential problem
nterms = 2 ! Number of terms in weak form
nitmx = 1000 ! Max number of iterations
rtolmx = 1.e-12 ! Max relative tolerance
nssor = 1 ! Number of SSOR precond steps
verbose = .FALSE. ! Output residue at each iteration
coefx(1:5) = (/1.0d0, 0.d0, 0.d0, 0.d0, 1.d0/) ! Mesh point distribution function
coefy(1:5) = (/1.0d0, 0.d0, 0.d0, 0.d0, 1.d0/) ! Mesh point distribution function
!
READ(*,newrun)
!
! Overwrite some input by reading in "filein" if required
!
IF(readmat) THEN
CALL openf(filein, fidin, mode='r')
CALL getatt(fidin, '/', 'NX', nints(1))
CALL getatt(fidin, '/', 'NY', nints(2))
CALL getatt(fidin, '/', 'NIDBAS1', nidbas(1))
CALL getatt(fidin, '/', 'NIDBAS2', nidbas(2))
CALL getatt(fidin, '/', 'NGAUSS1', ngauss(1))
CALL getatt(fidin, '/', 'NGAUSS2', ngauss(2))
CALL getatt(fidin, '/', 'MBESS', mbess)
END IF
WRITE(*,newrun)
!
! Define grid on x (=radial) & y (=poloidal) axis
!
pi = 4.0d0*ATAN(1.0d0)
ALLOCATE(xgrid(0:nints(1)), ygrid(0:nints(2)))
xgrid(0) = 0.0d0; xgrid(nints(1)) = 1.0d0
CALL meshdist(coefx, xgrid, nints(1))
ygrid(0) = 0.0d0; ygrid(nints(2)) = 2.d0*pi
CALL meshdist(coefy, ygrid, nints(2))
!
! Create hdf5 file
!
CALL creatf(file, fid, 'TCDSMAT Result File', real_prec='d')
CALL attach(fid, '/', 'NX', nints(1))
CALL attach(fid, '/', 'NY', nints(2))
CALL attach(fid, '/', 'NIDBAS1', nidbas(1))
CALL attach(fid, '/', 'NIDBAS2', nidbas(2))
CALL attach(fid, '/', 'NGAUSS1', ngauss(1))
CALL attach(fid, '/', 'NGAUSS2', ngauss(2))
CALL attach(fid, '/', 'MBESS', mbess)
!===========================================================================
! 2.0 Discretize the PDE
!
! Set up spline
!
WRITE(*,'(/a)') 'Discretize the PDE ...'
t0 = seconds()
CALL set_spline(nidbas, ngauss, &
& xgrid, ygrid, splxy, (/.FALSE., .TRUE./) )
!
! Set up CDS matrix for solver
!
CALL mstruct(nidbas, nints, mrank, dists)
bw0=(nidbas(1)+1)*nints(2) ! Half band including all zero diagonals
CALL init(mrank, dists, nterms, mat, bw0=bw0)
WRITE(*,'(a,4i8)') 'rank, kl, ku, bw0 = ', mat%rank, mat%kl, mat%ku, bw0
WRITE(*,'(i4,a/(10i8))') mat%ndiags, ' diagonals:', mat%dists
!
! FE matrix assembly and apply BC
!
IF(readmat) THEN
CALL getmat(fidin, '/MAT1', mat)
ELSE
CALL dismat(splxy, mat)
CALL ibcmat(mat, nints(2), nidbas(1))
END IF
!
WRITE(*,'(a,1pe12.3)') 'Mem used so far (MB)', mem()
!
! Assembly RHS and apply BC
!
ALLOCATE(rhs(mrank), sol(mrank))
!
IF(readmat) THEN
CALL getarr(fidin, '/RHS', rhs)
ELSE
CALL disrhs(mbess, splxy, rhs)
CALL ibcrhs(rhs, nints(2))
CALL putarr(fid, '/RHS', rhs, 'RHS of linear system')
END IF
tmat = seconds() - t0
!===========================================================================
! 3.0 Diagonal balance of matrix
!
WRITE(*,'(/a)') 'Diagonal balance of matrix ...'
tbal = seconds()
IF( .NOT. readmat ) THEN
CALL diagbal(mat)
CALL putmat(fid,'/MAT1', mat, 'CDS matrix with BC')
END IF
rhs = mat%bal * rhs
tbal = seconds()-tbal
!===========================================================================
! 4.0 Analytical solutions
!
ALLOCATE(solcal(0:nints(1),0:nints(2)))
ALLOCATE(solana(0:nints(1),0:nints(2)))
ALLOCATE(errsol(0:nints(1),0:nints(2)))
DO i=0,nints(1)
DO j=0,nints(2)
solana(i,j) = (1-xgrid(i)**2) * xgrid(i)**mbess * COS(mbess*ygrid(j))
END DO
END DO
!===========================================================================
! 5.0 Direct solve with MUMPS
!
!
WRITE(*,'(/a)') 'Solve the linear system using MUMPS ...'
!
tmumps(1) = seconds()
CALL cds2mumps(mat, mat_mumps)
CALL factor(mat_mumps, debug=.FALSE.)
tmumps(2) = seconds()
sol = rhs
CALL bsolve(mat_mumps, sol, debug=.FALSE.)
sol = mat%bal * sol
sol(1:nints(2)-1) = sol(nints(2)) ! Unicity
tmumps(1:2) = seconds()-tmumps(1:2)
!
tgrid = seconds()
dimx = splxy%sp1%dim
dimy = splxy%sp2%dim
ALLOCATE(bcoef(0:dimx-1, 0:dimy-1))
!
DO j=0,dimy-1
DO i=0,dimx-1
ij = MODULO(j,nints(2)) + i*nints(2) + 1
bcoef(i,j) = sol(ij)
END DO
END DO
CALL gridval(splxy, xgrid, ygrid, solcal, (/0,0/), bcoef)
errsol = solana - solcal
tgrid = seconds() - tgrid
!
PRINT*, 'Relative discretization errors', norm2(errsol) / norm2(solana)
!!$ WRITE(*,'(a,2(1pe15.6))') 'Relative discretization errors', &
!!$ & norm2(errsol) / norm2(solana)
!===========================================================================
! 5.0 Solve the linear system using CG
!
WRITE(*,'(/a)') 'Solve the linear system using CG ...'
!
tsolv = seconds()
sol(:) = 0.0d0 ! Initial guest for solution
IF( nssor .EQ. 0 ) THEN
CALL cg(mat, rhs, omega, nitmx, rtolmx, sol, resid, niter, &
& verbose=verbose)
ELSE
CALL cg(mat, rhs, omega, nitmx, rtolmx, sol, resid, niter, &
& verbose=verbose, nssor=nssor)
END IF
sol = mat%bal * sol
sol(1:nints(2)-1) = sol(nints(2)) ! Unicity
tsolv = seconds()-tsolv
!
DO j=0,dimy-1
DO i=0,dimx-1
ij = MODULO(j,nints(2)) + i*nints(2) + 1
bcoef(i,j) = sol(ij)
END DO
END DO
CALL gridval(splxy, xgrid, ygrid, solcal, (/0,0/), bcoef)
errsol = solana - solcal
!
PRINT*, 'Relative discretization errors', norm2(errsol) / norm2(solana)
!===========================================================================
! 9.0 Epilogue
!
CALL putarr(fid, '/SOL', sol)
CALL putarr(fid, '/xgrid', xgrid, 'r')
CALL putarr(fid, '/ygrid', ygrid, '\theta')
CALL putarr(fid, '/sol', solcal, 'Solutions')
CALL putarr(fid, '/solana', solana,'Exact solutions')
CALL putarr(fid, '/errors', errsol, 'Errors')
!
WRITE(*,'(/a)') '---'
WRITE(*,'(a,1pe12.3)') 'Matrice construction time (s) ', tmat
WRITE(*,'(a,1pe12.3)') 'Matrice balancing time (s) ', tbal
WRITE(*,'(a,2(1pe12.3))') 'MUMPS solver time (s) ', tmumps
WRITE(*,'(a,1pe12.3)') 'Solution at grid time (s) ', tgrid
WRITE(*,'(a,i8,2(1pe12.3))') 'nits, resid, t(s)', niter, resid, tsolv
!
DEALLOCATE(xgrid,ygrid)
DEALLOCATE(rhs, sol)
DEALLOCATE(dists)
DEALLOCATE(bcoef)
DEALLOCATE(solcal,solana,errsol)
CALL destroy_sp(splxy)
CALL destroy(mat)
CALL closef(fid)
IF(readmat) THEN
CALL closef(fidin)
END IF
!
CALL mpi_finalize(ierr)
!===========================================================================
END PROGRAM main
c4science · Help