Homec4science
Diffusion SPClibs (feature/refactoring) Browse

spclibs/examples/pde1d_eig_zmumps.f90
03540012fedefeature/refactoring

Last Change
Raw File

pde1d_eig_zmumps.f90
View Options

!>
!> @file pde1d_eig_zmumps.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 1d differential eqation using splines:
!
! Solve the standard eigenvalue:
! A*x = \lambda *x or inv(A)*x = 1/\lambda * x using Arpack and MUMPS.
! where A is obtained from discretozation of
! -d/dr[r d/dr] f = k^2 r^(k-1)
!
USE bsplines
USE mumps_bsplines
USE futils
IMPLICIT NONE
INTEGER :: nx, nidbas, ngauss, kdiff
INTEGER :: i, nrank
LOGICAL :: nlppform
DOUBLE PRECISION :: coefs(5)
DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: xgrid
TYPE(spline1d) :: splx
TYPE(zmumps_mat) :: mat
!
INTEGER :: ierr
INTEGER :: fid
CHARACTER(len=32) :: str
DOUBLE PRECISION :: seconds, t0, tmat, tfac, tarpack
!
! Arpack: Solve the standard eigenvalue problem
!
INTEGER :: nev = 10, ncv = 10
LOGICAL :: nlinv = .FALSE. ! Solve inv(A) = 1/\lambda * x if nlinv=.TRUE.
CHARACTER(len=2) :: which='SM'
INTEGER :: info=0 ! Use random vector to start the Arnoldi iterations
INTEGER :: ido=0 ! Reverse communications
LOGICAL :: rvec
LOGICAL, ALLOCATABLE :: select(:)
INTEGER :: iparam(11), ipntr(14), nconv
DOUBLE PRECISION :: tol=0.0d0
CHARACTER(len=1) :: bmat='I'
!
INTEGER :: lworkl
DOUBLE COMPLEX, ALLOCATABLE :: workl(:), workd(:), workev(:)
DOUBLE COMPLEX, ALLOCATABLE :: eig_vals(:), eig_vecs(:,:), resid(:)
DOUBLE COMPLEX :: sigma
DOUBLE PRECISION, ALLOCATABLE :: rwork(:)
!
INTERFACE
SUBROUTINE dismat(spl, mat)
USE bsplines
USE mumps_bsplines
TYPE(spline1d), INTENT(in) :: spl
TYPE(zmumps_mat), INTENT(inout) :: mat
END SUBROUTINE dismat
SUBROUTINE disrhs(kdiff, spl, rhs)
USE bsplines
INTEGER, INTENT(in) :: kdiff
TYPE(spline1d), INTENT(in) :: spl
DOUBLE PRECISION, INTENT(out) :: rhs(:)
END SUBROUTINE disrhs
SUBROUTINE meshdist(coefs, x, nx)
DOUBLE PRECISION, INTENT(in) :: coefs(5)
INTEGER, INTENT(iN) :: nx
DOUBLE PRECISION, INTENT(inout) :: x(0:nx)
END SUBROUTINE meshdist
END INTERFACE
!
NAMELIST /newrun/ nx, nidbas, ngauss, kdiff, nlppform, coefs, &
& nev, ncv, nlinv, which, tol
!===========================================================================
! 1.0 Prologue
!
CALL mpi_init(ierr)
!
! 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
kdiff = 2 ! Exponent of differential problem
nlppform = .TRUE. ! Use PPFORM for gridval or not
coefs(1:5) = (/0.2d0, 1.d0, 0.d0, 0.d0, 1.d0/) ! Mesh point distribution function
! see function FDIST in MESHDIST
!
READ(*,newrun)
WRITE(*,newrun)
!
! Define grid on x axis
!
ALLOCATE(xgrid(0:nx))
xgrid(0) = 0.0d0
xgrid(nx) = 1.0d0
CALL meshdist(coefs, xgrid, nx)
WRITE(*,'(a/(10f8.3))') 'XGRID', xgrid(0:nx)
WRITE(*,'(a/(10f8.3))') 'Diff XGRID', (xgrid(i)-xgrid(i-1), i=1,nx)
WRITE(*,'(a/(10f8.3))') 'Diff XGRID**2', (xgrid(i)**2-xgrid(i-1)**2, i=1,nx)
!===========================================================================
! 2.0 Discretize the PDE
!
! Set up spline
!
t0 = seconds()
CALL set_spline(nidbas, ngauss, xgrid, splx, nlppform=nlppform)
CALL get_dim(splx, nrank) ! Rank of the FE matrix
WRITE(*,'(a,l6)') 'splx%nlequid', splx%nlequid
!
! FE matrix assembly
!
WRITE(*,'(a,i6)') ' nrank', nrank
CALL init(nrank, 1, mat)
CALL dismat(splx, mat)
CALL to_mat(mat)
tmat = seconds() - t0
!
mat%mumps_par%IRN => mat%mumps_par%IRN_loc ! Work around for single proc.
PRINT*, 'nnz_loc', mat%nnz_loc
PRINT*, 'mat%mumps_par%N', mat%mumps_par%N
PRINT*, 'mat%mumps_par%NZ_loc', mat%mumps_par%NZ_loc
PRINT*, 'size of mat%mumps_par%IRN', SIZE(mat%mumps_par%IRN)
PRINT*, 'mat%istart,mat%iend', mat%istart,mat%iend
!
CALL creatf('pde1d_eig_zmumps.h5', fid, 'PDE1D Result File', real_prec='d')
PRINT*, 'rank', mat%rank
PRINT*, 'nnz', mat%nnz
!
CALL putmat(fid,'/MAT',mat,'FE matrix')
!
IF(nlinv) THEN
t0 = seconds()
CALL factor(mat)
tfac = seconds()-t0
END IF
!===========================================================================
! 3.0 Solve the standard eigenvalue problem
!
lworkl = 3*ncv**2 + 5*ncv
ALLOCATE(workl(lworkl))
ALLOCATE(workd(3*nrank))
ALLOCATE(workev(2*ncv))
ALLOCATE(eig_vals(ncv), eig_vecs(nrank,ncv))
ALLOCATE(resid(nrank))
ALLOCATE(rwork(ncv))
!
iparam(1) = 1 ! shfts
iparam(3) = 300 ! Max. number of iterations
iparam(7) = 1 ! Regular mode
!
! The reverse communication loop
!
t0 = seconds()
DO
CALL znaupd (ido, bmat, nrank, which, nev, tol, resid, ncv, &
& eig_vecs, nrank, iparam, ipntr, workd, workl, lworkl, &
& rwork, info )
!
IF(ido .EQ. -1 .OR. ido .EQ. 1) THEN ! Compute A*v
CALL av(nrank, workd(ipntr(1)), workd(ipntr(2)))
CYCLE
END IF
!
IF(info .LT. 0) THEN ! Error
PRINT*, 'Error in _naupd with info =', info
EXIT
ELSE
rvec = .TRUE.
ALLOCATE(select(ncv))
CALL zneupd (rvec, 'A', select, eig_vals, eig_vecs, nrank, &
& sigma, workev, bmat, nrank, which, nev, tol, resid, &
& ncv, eig_vecs, nrank, iparam, ipntr, workd, workl, lworkl,&
& rwork, ierr)
IF(ierr .NE. 0) THEN
PRINT*, 'Error in _neupd with ierr =', ierr
EXIT
ELSE
nconv = iparam(5)
PRINT*,'Number of converged eigenvalues', nconv
IF(nlinv) THEN
eig_vals(1:nconv) = (1.d0,0.0d0) / eig_vals(1:nconv)
END IF
WRITE(*,'(2(1pe12.3))') eig_vals(1:nconv)
CALL putarr(fid, '/eig_vals', eig_vals(1:nconv))
CALL putarr(fid, '/eig_vecs', eig_vecs(1:nrank,1:nconv))
DO i=1,nconv
!!$ WRITE(*,'(/a,2(pe20.6))') '*** eigen value =', eig_vals(i)
!!$ WRITE(*,'(a/(10(1pe12.4)))') 'Real of eigen vector', &
!!$ & REAL(eig_vecs(1:nrank,i))
!!$ WRITE(*,'(a/(10(1pe12.4)))') 'Imag of eigen vector', &
!!$ & aimag(eig_vecs(1:nrank,i))
WRITE(str,'(a,i3.3)') '/eig_vecs_',i
CALL putarr(fid, TRIM(str), eig_vecs(1:nrank,i))
END DO
EXIT
END IF
END IF
END DO ! End of reverse commuinication loop
IF(info .EQ. 1) THEN
PRINT*, 'Maximum number of iterations reached!'
PRINT*, 'IPARAM(5) =', iparam(5)
END IF
PRINT*, 'Number of Arnoldi iterations', iparam(3)
tarpack = seconds() - t0
!===========================================================================
! 9.0 Epilogue
!
WRITE(*,'(a)') '---'
WRITE(*,'(a,1pe12.3)') 'Matrice construction time (s) ', tmat
IF(nlinv) THEN
WRITE(*,'(a,1pe12.3)') 'Matrice factorization time (s) ', tfac
END IF
WRITE(*,'(a,1pe12.3)') 'Arpack time (s) ', tarpack
!
DEALLOCATE(xgrid)
CALL destroy(mat)
CALL destroy_sp(splx)
CALL closef(fid)
CALL mpi_finalize(ierr)
!
CONTAINS
SUBROUTINE av(n,v,w)
!
INTEGER, INTENT(in) :: n
DOUBLE COMPLEX, INTENT(in) :: v(*)
DOUBLE COMPLEX, INTENT(out) :: w(*)
!
IF(nlinv) THEN
w(1:n) = v(1:n)
CALL bsolve(mat,w(1:n)) ! Solve A*w = v or w=inv(A)*v
ELSE
w(1:n) = vmx(mat,v(1:n)) ! A*v matrix product
END IF
END SUBROUTINE av
!
SUBROUTINE coefeq(x, idt, idw, c)
DOUBLE PRECISION, INTENT(in) :: x
INTEGER, INTENT(out) :: idt(:), idw(:)
DOUBLE PRECISION, INTENT(out) :: c(:)
!
! Weak form = Int(x*dw/dx*dt/dx)dx
!
c(1) = x
idt(1) = 1
idw(1) = 1
END SUBROUTINE coefeq
END PROGRAM main
!
!+++
SUBROUTINE meshdist(c, x, nx)
!
! Construct an 1d non-equidistant mesh given a
! mesh distribution function defined in FDIST
!
IMPLICIT NONE
DOUBLE PRECISION, INTENT(in) :: c(5)
INTEGER, INTENT(iN) :: nx
DOUBLE PRECISION, INTENT(inout) :: x(0:nx)
INTEGER :: nintg
DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: xint, fint
DOUBLE PRECISION :: a, b, dx, f0, f1, scal
INTEGER :: i, k
!
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
!!$ WRITE(*,'(a/(10f8.3))') 'FINT', fint
!
! 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 = c(1) + c(2)*x + c(3)*EXP(-((x-c(4))/c(5))**2)
END FUNCTION fdist
END SUBROUTINE meshdist
!+++
SUBROUTINE dismat(spl, mat)
!
! Assembly FE matrix mat using spline spl
!
USE bsplines
USE mumps_bsplines
IMPLICIT NONE
TYPE(spline1d), INTENT(in) :: spl
TYPE(zmumps_mat), INTENT(inout) :: mat
INTEGER :: nrank, 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 PRECISION, ALLOCATABLE :: coefs(:)
!===========================================================================
! 1.0 Prologue
!
! Properties of spline space
!
CALL get_dim(spl, nrank, 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=i+iw; jcol=i+jt
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 PRECISION, INTENT(out) :: c(SIZE(idt))
!
! Weak form = Int(x*dw/dx*dt/dx)dx
!
c(1) = x
idt(1) = 1
idw(1) = 1
END SUBROUTINE coefeq
END SUBROUTINE dismat
!+++
SUBROUTINE disrhs(kdiff, spl, rhs)
!
! Assenbly the RHS using spline spl
!
USE bsplines
IMPLICIT NONE
INTEGER, INTENT(in) :: kdiff
TYPE(spline1d), INTENT(in) :: spl
DOUBLE PRECISION, INTENT(out) :: rhs(:)
INTEGER :: nrank, nx, nidbas, ngauss
INTEGER :: i, igauss, left
DOUBLE PRECISION, ALLOCATABLE :: fun(:,:)
DOUBLE PRECISION, ALLOCATABLE :: xgauss(:), wgauss(:)
DOUBLE PRECISION:: contrib
!===========================================================================
! 1.0 Prologue
!
! Properties of spline space
!
CALL get_dim(spl, nrank, nx, nidbas)
!!$ WRITE(*,'(/a, 5i3)') 'nrank, nx, nidbas =', nrank, nx, nidbas
!
ALLOCATE(fun(nidbas+1,1)) ! needs only basis functions (no derivatives)
!
! Gauss quadature
!
CALL get_gauss(spl, ngauss)
!!$ WRITE(*,'(/a, i3)') 'Gauss points and weights, ngauss =', ngauss
ALLOCATE(xgauss(ngauss), wgauss(ngauss))
!===========================================================================
! 2.0 Assembly loop
!
rhs(1:nrank) = 0.0d0
!
DO i=1,nx
CALL get_gauss(spl, ngauss, i, xgauss, wgauss)
!!$ WRITE(*,'(a,i3,(8f10.4))') 'i=',i, xgauss, wgauss
DO igauss=1,ngauss
CALL basfun(xgauss(igauss), spl, fun, i)
left = i-1
!!$ WRITE(*,'(a,5i4)') '>>>igauss, left', igauss, left
contrib = wgauss(igauss)*rhseq(xgauss(igauss), kdiff)
rhs(i:i+nidbas) = rhs(i:i+nidbas) + contrib*fun(:,1)
END DO
END DO
!===========================================================================
! 9.0 Epilogue
!
DEALLOCATE(fun)
DEALLOCATE(xgauss, wgauss)
!
CONTAINS
DOUBLE PRECISION FUNCTION rhseq(x,k)
DOUBLE PRECISION, INTENT(in) :: x
INTEGER, INTENT(in) :: k
rhseq = k*k*x**(k-1)
END FUNCTION rhseq
END SUBROUTINE disrhs
c4science · Help