Homec4science
Diffusion SPClibs (feature/refactoring) Browse

spclibs/examples/pde1dp_cmpl_pardiso.f90
c723281d03b6feature/refactoring

Last Change
Raw File

pde1dp_cmpl_pardiso.f90
View Options

!>
!> @file pde1dp_cmpl_pardiso.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>
!>
MODULE pde1dp_cmpl_pardiso_mod
USE bsplines
USE pardiso_bsplines
IMPLICIT NONE
!
CONTAINS
SUBROUTINE disrhs(spl, rhs, mmode, alpha, beta)
!
TYPE(spline1d) :: spl
INTEGER, INTENT(in) :: mmode
DOUBLE COMPLEX, INTENT(in) :: alpha, beta
DOUBLE COMPLEX, INTENT(out) :: rhs(:)
!
INTEGER :: dim, nrank, nx, nidbas, ngauss
INTEGER :: i, igauss, it, irow
DOUBLE PRECISION, ALLOCATABLE :: fun(:,:)
DOUBLE PRECISION, ALLOCATABLE :: xgauss(:), wgauss(:)
DOUBLE COMPLEX :: contrib
!===========================================================================
! 1.0 Prologue
!
! Properties of spline space
!
nrank = SIZE(rhs)
CALL get_dim(spl, dim, nx, nidbas)
!
ALLOCATE(fun(0:nidbas,1)) ! needs only basis functions (no derivatives)
!
! Gauss quadature
!
CALL get_gauss(spl, ngauss)
ALLOCATE(xgauss(ngauss), wgauss(ngauss))
!===========================================================================
! 2.0 Assembly loop
!
rhs(:) = 0.0d0
!
DO i=1,nx
CALL get_gauss(spl, ngauss, i, xgauss, wgauss)
DO igauss=1,ngauss
CALL basfun(xgauss(igauss), spl, fun, i)
contrib = wgauss(igauss) * rhseq(xgauss(igauss))
DO it=0,nidbas
irow=MODULO(i+it-1,nx) + 1 ! Periodic BC
rhs(irow) = rhs(irow) + contrib*fun(it,1)
END DO
END DO
END DO
!===========================================================================
! 9.0 Epilogue
!
DEALLOCATE(fun)
DEALLOCATE(xgauss, wgauss)
CONTAINS
DOUBLE COMPLEX FUNCTION rhseq(x)
DOUBLE PRECISION, INTENT(in) :: x
DOUBLE PRECISION :: arg
arg = mmode*x
rhseq = (mmode**2*alpha-beta)*COS(arg)
END FUNCTION rhseq
END SUBROUTINE disrhs
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE dismat(spl, mat, alpha, beta)
!
TYPE(spline1d) :: spl
TYPE(zpardiso_mat) :: mat
DOUBLE COMPLEX, INTENT(in) :: alpha, beta
!
INTEGER :: dim, nx, nidbas, ngauss
INTEGER :: i, igauss, iterm, iw, jt, irow, jcol
DOUBLE PRECISION, ALLOCATABLE :: fun(:,:)
DOUBLE PRECISION, ALLOCATABLE :: xgauss(:), wgauss(:)
DOUBLE COMPLEX :: contrib
!
INTEGER :: kterms ! Number of terms in weak form
INTEGER, ALLOCATABLE :: idert(:), iderw(:) ! Derivative order
DOUBLE COMPLEX, ALLOCATABLE :: coefs(:)
!===========================================================================
! 1.0 Prologue
!
! Properties of spline space
!
CALL get_dim(spl, dim, nx, nidbas)
ALLOCATE(fun(0:nidbas,0:1)) ! Spline and its 1st derivative
!
! Weak form
!
kterms = mat%nterms
ALLOCATE(idert(kterms), iderw(kterms), coefs(kterms))
!
! Gauss quadature
!
CALL get_gauss(spl, ngauss)
ALLOCATE(xgauss(ngauss), wgauss(ngauss))
!===========================================================================
! 2.0 Assembly loop
!
DO i=1,nx
CALL get_gauss(spl, ngauss, i, xgauss, wgauss)
DO igauss=1,ngauss
CALL basfun(xgauss(igauss), spl, fun, i)
CALL coefeq(xgauss(igauss), idert, iderw, coefs)
DO iterm=1,kterms
DO jt=0,nidbas
DO iw=0,nidbas
contrib = fun(jt,idert(iterm)) * coefs(iterm) * &
& fun(iw,iderw(iterm)) * wgauss(igauss)
irow=MODULO(i+iw-1,nx) + 1 ! Periodic BC
jcol=MODULO(i+jt-1,nx) + 1
CALL updtmat(mat, irow, jcol, contrib)
END DO
END DO
END DO
END DO
END DO
!===========================================================================
! 9.0 Epilogue
!
DEALLOCATE(fun)
DEALLOCATE(xgauss, wgauss)
DEALLOCATE(iderw, idert, coefs)
!
CONTAINS
SUBROUTINE coefeq(x, idt, idw, c)
DOUBLE PRECISION, INTENT(in) :: x
INTEGER, INTENT(out) :: idt(:), idw(:)
DOUBLE COMPLEX, INTENT(out) :: c(SIZE(idt))
!
c(1) = alpha
idt(1) = 1
idw(1) = 1
!
c(2) = -beta
idt(2) = 0
idw(2) = 0
END SUBROUTINE coefeq
END SUBROUTINE dismat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
FUNCTION norm2(x)
!
! Compute the 2-norm of complex array x
!
IMPLICIT NONE
DOUBLE PRECISION :: norm2
DOUBLE COMPLEX, INTENT(in) :: x(:)
DOUBLE PRECISION :: sum2
!
sum2 = DOT_PRODUCT(x,x)
norm2 = SQRT(sum2)
END FUNCTION norm2
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE meshdist(mmode, x)
!
! Construct a 1d non-equidistant mesh given a
! mesh distribution function.
!
INTEGER, INTENT(in) :: mmode
DOUBLE PRECISION, INTENT(inout) :: x(0:)
INTEGER :: nx, nintg
DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: xint, fint
DOUBLE PRECISION :: a, b, dx, f0, f1, scal
INTEGER :: i, k
!
nx = SIZE(x)-1
a=x(0)
b=x(nx)
nintg = 10*nx
ALLOCATE(xint(0:nintg), fint(0:nintg))
!
! Mesh distribution
!
dx = (b-a)/REAL(nintg)
xint(0) = a
fint(0) = 0.0d0
f1 = fdist(xint(0))
DO i=1,nintg
f0 = f1
xint(i) = xint(i-1) + dx
f1 = fdist(xint(i))
fint(i) = fint(i-1) + 0.5*(f0+f1)
END DO
!
! Normalization
!
scal = REAL(nx) / fint(nintg)
fint(0:nintg) = fint(0:nintg) * scal
!
! Obtain mesh point by (inverse) interpolation
!
k = 1
DO i=1,nintg-1
IF( fint(i) .GE. REAL(k) ) THEN
x(k) = xint(i) + (xint(i+1)-xint(i))/(fint(i+1)-fint(i)) * &
& (k-fint(i))
k = k+1
END IF
END DO
!
DEALLOCATE(xint, fint)
CONTAINS
DOUBLE PRECISION FUNCTION fdist(x)
DOUBLE PRECISION, INTENT(in) :: x
fdist = 2.0 + COS(mmode*x)
END FUNCTION fdist
END SUBROUTINE meshdist
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
END MODULE pde1dp_cmpl_pardiso_mod
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
PROGRAM main
!
! 1D complex PDE with periodic BC
!
USE pde1dp_cmpl_pardiso_mod
USE futils
USE conmat_mod
!
IMPLICIT NONE
TYPE(spline1d) :: splx
TYPE(zpardiso_mat) :: mat
TYPE(zpardiso_mat) :: newmat
INTEGER :: kl, ku, nrank
!
CHARACTER(len=128) :: file='pde1dp_cmpl_pardiso.h5'
INTEGER :: fid
INTEGER :: nx, nidbas, ngauss, mmode, npt, dim
LOGICAL :: nlequid
LOGICAL :: nlsym, nlherm, nlpos
DOUBLE PRECISION, PARAMETER :: pi=3.1415926535897931d0
DOUBLE PRECISION :: dx
DOUBLE COMPLEX :: alpha, beta
DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: xgrid, x, err
DOUBLE COMPLEX, DIMENSION(:), ALLOCATABLE :: sol, rhs, bcoef
DOUBLE COMPLEX, DIMENSION(:), ALLOCATABLE :: newsol, arow
DOUBLE COMPLEX, DIMENSION(:), ALLOCATABLE :: solcal, solana
DOUBLE PRECISION, ALLOCATABLE :: xgauss(:), wgauss(:)
DOUBLE PRECISION :: err_norm
INTEGER :: i
!
NAMELIST /newrun/ nx, nidbas, ngauss, nlequid, alpha, beta, mmode, npt, &
& nlsym, nlherm, nlpos
!===========================================================================
! 1.0 Prologue
!
! Read in data specific to run
!
nx = 8 ! Number oh intevals in x
nidbas = 3 ! Degree of splines
ngauss = 4 ! Number of Gauss points/interval
nlequid = .TRUE. ! Use exact sol. as mesh dist. function if .FALSE.
mmode = 1 ! Fourier mode
alpha = (1.0, 1.0) ! Complex "diffusion"
beta = 1.0
npt = 100
nlsym = .TRUE. ! Is matrice symmetric
nlherm = .FALSE. ! Is matrice hermitian
nlpos = .TRUE. ! and positive definite ?
!
READ(*,newrun)
WRITE(*,newrun)
!
! Define grid on x axis
!
ALLOCATE(xgrid(0:nx))
dx = 2.d0*pi/REAL(nx,8)
xgrid = (/ (i*dx,i=0,nx) /)
IF( .NOT. nlequid ) THEN
CALL meshdist(mmode, xgrid)
END IF
WRITE(*,'(a/(10f8.3))') 'xgrid', xgrid
!
! Create hdf5 file
!
CALL creatf(file, fid, 'PDE1DP Result File', real_prec='d')
CALL attach(fid, '/', 'NX', nx)
CALL attach(fid, '/', 'NIDBAS', nidbas)
CALL putarr(fid, '/xgrid', xgrid)
!===========================================================================
! 2.0 Discretize the PDE
!
! Set up periodic spline
!
CALL set_spline(nidbas, ngauss, xgrid, splx, period=.TRUE.)
WRITE(*,'(a,l6)') 'nlequid =', nlequid
nrank = nx ! Rank of the FE matrix
!
! FE matrix assembly
!
CALL init(nrank, 2, mat, nlsym=nlsym, nlherm=nlherm, nlpos=nlpos)
CALL get_dim(splx, dim)
WRITE(*,'(/a,4i6)') 'nrank, dim', nrank, dim
!!$ CALL dismat(splx, mat, alpha, beta)
CALL conmat(splx, mat, coefeq)
!
! RHS assembly
!
ALLOCATE(rhs(nrank), sol(nrank))
CALL disrhs(splx, rhs, mmode, alpha, beta)
!
CALL putarr(fid, '/rhs', rhs)
!
! Factor and solve
!
WRITE(*,'(a/(10i6))') 'iparm', mat%p%iparm
CALL factor(mat)
CALL putmat(fid,'/MAT', mat)
CALL bsolve(mat, rhs, sol)
CALL putarr(fid, '/sol', sol)
WRITE(*,'(/a,i8)') 'Number of non-zeros of matrix = ',get_count(mat)
WRITE(*,'(a,i8)') 'Number of nonzeros in factors of A = ',mat%p%iparm(18)
!
! Compute residue
!
WRITE(*,'(a,1pe12.4)') 'Residue =', norm2(vmx(mat,sol)-rhs)
!===========================================================================
! 3.0 Check solution
!
! Exact solution
ALLOCATE(x(0:npt), solcal(0:npt), solana(0:npt), err(0:npt))
dx=2.0d0*pi/REAL(npt,8)
x = (/ (i*dx, i=0,npt) /)
solana = COS(mmode*x)
!
! Prolongate solution using periodicity
!
ALLOCATE(bcoef(dim))
bcoef(1:nrank) = sol(1:nrank)
DO i=nrank+1,dim
bcoef(i) = bcoef(MODULO(i-1,nrank)+1)
END DO
!
! Interpolate field
!
CALL gridval(splx, x, solcal, 0, bcoef)
!
err = ABS(solcal-solana)
CALL putarr(fid, '/x', x)
CALL putarr(fid, '/solana', solana)
CALL putarr(fid, '/solcal', solcal)
CALL putarr(fid, '/err', err)
!
! Compute discretization error norm by Gauss integration
!
err_norm=0.0
ALLOCATE(xgauss(ngauss), wgauss(ngauss))
DO i=1,nx
CALL get_gauss(splx, ngauss, i, xgauss, wgauss)
CALL gridval(splx, xgauss, solcal(1:ngauss), 0)
solana(1:ngauss) = COS(mmode*xgauss)
err(1:ngauss) = DOT_PRODUCT(solana(1:ngauss)-solcal(1:ngauss), &
& solana(1:ngauss)-solcal(1:ngauss))
err_norm = err_norm + SUM(wgauss*err(1:ngauss))
END DO
err_norm = SQRT(err_norm)
WRITE(*,'(a,1pe12.3)') 'Discretization error ', err_norm
!
!===========================================================================
! 4.0 Test of getrow/putrow, getcol/putcol and mcopy
!
CALL init(nrank, 2, newmat, nlsym=nlsym, nlherm=nlherm, nlpos=nlpos)
ALLOCATE(arow(nrank), newsol(nrank))
!
DO i=1,nrank
CALL getrow(mat, i, arow)
CALL putrow(newmat, i, arow)
END DO
CALL factor(newmat)
CALL bsolve(newmat, rhs, newsol)
WRITE(*,'(/a)') 'putrow/getrow ...'
WRITE(*,'(a,1pe12.4)') 'Residue =', norm2(vmx(newmat,newsol)-rhs)
WRITE(*,'(a,1pe12.3)') 'Error ', norm2(sol-newsol)
!
DO i=1,nrank
CALL getcol(mat, i, arow)
CALL putcol(newmat, i, arow)
END DO
CALL factor(newmat)
CALL bsolve(newmat, rhs, newsol)
WRITE(*,'(/a)') 'putcol/getcol ...'
WRITE(*,'(a,1pe12.4)') 'Residue =', norm2(vmx(newmat,newsol)-rhs)
WRITE(*,'(a,1pe12.3)') 'Error ', norm2(sol-newsol)
!
CALL clear_mat(newmat)
CALL mcopy(mat, newmat)
WRITE(*,'(/a)') 'mcopy ...'
newmat%val = (1000.0d0,0.0d0)*newmat%val
CALL factor(newmat)
CALL bsolve(newmat, rhs, newsol)
WRITE(*,'(a)') 'Backsolve the new system'
WRITE(*,'(a,1pe12.4)') 'Norm =', norm2(newsol)
!
WRITE(*,'(a)') 'Destroy NEWMAT ...'
CALL destroy(newmat)
!
CALL bsolve(mat, rhs, sol)
WRITE(*,'(/a)') 'Backsolve the old system'
WRITE(*,'(a,1pe12.4)') 'Norm =', norm2(sol)
!
WRITE(*,'(a)') 'Destroy MAT ...'
CALL destroy(mat)
!!$!
!!$ WRITE(*,'(/a)') 'Should crash since NEWMAT is gone!'
!!$ CALL bsolve(newmat, rhs, newsol)
!!$ WRITE(*,'(a)') 'Backsolve the new system'
!!$ WRITE(*,'(a,1pe12.4)') 'Norm =', norm2(newsol)
!===========================================================================
! 9.0 Clean up
!
DEALLOCATE(x, solcal, solana, err)
DEALLOCATE(xgauss, wgauss)
DEALLOCATE(bcoef)
DEALLOCATE(xgrid)
DEALLOCATE(rhs, sol)
DEALLOCATE(arow, newsol)
CALL destroy_sp(splx)
CALL closef(fid)
CONTAINS
SUBROUTINE coefeq(x, idt, idw, c)
DOUBLE PRECISION, INTENT(in) :: x
INTEGER, INTENT(out) :: idt(:), idw(:)
DOUBLE COMPLEX, INTENT(out) :: c(:)
!
c(1) = alpha
idt(1) = 1
idw(1) = 1
!
c(2) = -beta
idt(2) = 0
idw(2) = 0
END SUBROUTINE coefeq
END PROGRAM main
c4science · Help