Homec4science
Diffusion SPClibs (feature/refactoring) Browse

spclibs/examples/tmassmat.f90
03540012fedefeature/refactoring

Last Change
Raw File

tmassmat.f90
View Options

!>
!> @file tmassmat.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
!
! FT of mass matrix
!
USE bsplines
USE matrix
!
IMPLICIT NONE
INTEGER :: nx, nidbas
INTEGER :: ngauss, nrank, kl, ku
INTEGER :: i, k, kmin, kmax
DOUBLE PRECISION :: pi, xlenght, dx, arg0, arg
DOUBLE PRECISION, ALLOCATABLE :: xgrid(:), arow(:)
DOUBLE PRECISION, ALLOCATABLE :: fftmassm1(:), fftmassm2(:), fftmassm3(:)
DOUBLE PRECISION, ALLOCATABLE :: fftmass_shifted(:)
TYPE(spline1d) :: splx
TYPE(periodic_mat) :: massm
!
NAMELIST /newrun/ nx, nidbas, xlenght
!================================================================================
! 1.0 Prologue
!
pi = 4.0d0*ATAN(1.0d0)
!
nx = 8
nidbas = 3
xlenght = 2.0d0*pi
!
READ(*,newrun)
WRITE(*,newrun)
!
ngauss = nidbas+1 ! Exact integration for polynomials of degree 2*nidbas
!
ALLOCATE(xgrid(0:nx))
dx = xlenght/REAL(nx)
xgrid = (/ (i*dx,i=0,nx) /)
WRITE(*,'(a/(10f8.3))') 'xgrid', xgrid
!
CALL set_spline(nidbas, ngauss, xgrid, splx, .TRUE.)
!===========================================================================
! 2.0 Mass matrix
nrank = nx
kl = nidbas
ku = kl
CALL init(kl, ku, nrank, 1, massm)
CALL dismat(splx, massm)
!
ALLOCATE(arow(nrank))
!!$ WRITE(*,'(/a)') 'Mass matrix'
!!$ DO i=1,nrank
!!$ CALL getrow(massm, i, arow)
!!$ WRITE(*,'(10(1pe12.4))') arow
!!$ END DO
!===========================================================================
! 3.0 Fourier transform of Mass matrix
!
ALLOCATE(fftmassm1(0:nx-1))
ALLOCATE(fftmassm2(0:nx-1))
ALLOCATE(fftmassm3(0:nx-1))
IF(nidbas.LE.3) THEN
CALL analytic(nidbas, fftmassm1)
fftmassm1 = dx*fftmassm1
WRITE(*,'(/a/(10(1pe12.4)))') 'Analytic', fftmassm1
END IF
!
CALL calc_fftmass_old(splx, fftmassm2)
WRITE(*,'(/a/(10(1pe12.4)))') 'Old version', fftmassm2
IF(nidbas.LE.3) THEN
WRITE(*,'(a,1pe12.4)') 'error =', MAXVAL(ABS(fftmassm2-fftmassm1))
END IF
!
! Init DFT
kmin = -nx/2
kmax = nx/2-1
CALL init_dft(splx, kmin, kmax)
!
ALLOCATE(fftmass_shifted(kmin:kmax))
CALL calc_fftmass(splx, fftmass_shifted)
DO k=kmin, kmax
fftmassm3(MODULO(k+nx,nx))=fftmass_shifted(k)
END DO
WRITE(*,'(/a/(10(1pe12.4)))') 'New version', fftmassm3
WRITE(*,'(a,1pe12.4)') 'error =', MAXVAL(ABS(fftmassm3-fftmassm2))
!
! Check dftcoefs
WRITE(*,'(/a)') 'Check DFT of splines'
PRINT*, 'dims of dftcoefs', SHAPE(splx%dft%coefs)
DO i=0,nidbas
WRITE(*,'(a,i3,2(1pe12.4))') 'Sum of coefs for spline', i, &
& SUM(splx%dft%coefs(:,i))
END DO
!===========================================================================
! 9.0 Clean up
!
DEALLOCATE(xgrid)
DEALLOCATE(fftmassm1)
DEALLOCATE(fftmassm2)
DEALLOCATE(fftmassm3)
DEALLOCATE(fftmass_shifted)
DEALLOCATE(arow)
CALL destroy(massm)
CALL destroy_sp(splx)
CONTAINS
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE dismat(spl, mat)
TYPE(spline1d) :: spl
TYPE(periodic_mat) :: mat
INTEGER :: dim, nx, nidbas, ngauss
DOUBLE PRECISION, ALLOCATABLE :: fun(:,:), xgauss(:), wgauss(:)
INTEGER :: i, igauss, iw, jt, irow, jcol
DOUBLE PRECISION :: contrib
!
CALL get_dim(spl, dim, nx, nidbas)
ALLOCATE(fun(0:nidbas,1)) ! Spline
CALL get_gauss(spl, ngauss)
ALLOCATE(xgauss(ngauss), wgauss(ngauss))
!
DO i=1,nx
CALL get_gauss(spl, ngauss, i, xgauss, wgauss)
DO igauss=1,ngauss
CALL basfun(xgauss(igauss), spl, fun, i)
DO jt=0,nidbas
DO iw=0,nidbas
contrib = fun(jt,1) * fun(iw,1) * 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
!
DEALLOCATE(fun)
DEALLOCATE(xgauss, wgauss)
END SUBROUTINE dismat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE analytic(nidbas, mat)
!
! Analytic form for nidbas .le. 3
!
INTEGER, INTENT(in) :: nidbas
DOUBLE PRECISION, INTENT(out) :: mat(0:)
DOUBLE PRECISION :: arg0, arg, cosk
INTEGER :: n, k
!
n = SIZE(mat)
arg0 = 2.0d0*pi/REAL(n,8)
DO k=0,n-1
arg = k*arg0
cosk = COS(arg)
SELECT CASE (nidbas)
CASE (1)
mat(k) = (2.0d0 + cosk)/3.0d0
CASE (2)
mat(k) = (16.0d0 + cosk*(13.0d0+cosk))/30.0d0
CASE (3)
mat(k) = (272.0d0 + cosk*(297.0d0+cosk*(60.0d0+cosk)))/630.0d0
CASE default
WRITE(*,'(a,i4,a)') 'ANALYTIC: nidbas =', nidbas, ' is not implemented!'
STOP
END SELECT
END DO
END SUBROUTINE analytic
END PROGRAM main
c4science · Help