Homec4science
Diffusion SPClibs (feature/refactoring) Browse

spclibs/examples/pde1dp_cmpl_dft.f90
03540012fedefeature/refactoring

Last Change
Raw File

pde1dp_cmpl_dft.f90
View Options

!>
!> @file pde1dp_cmpl_dft.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>
!>
MODULE pde1dp_cmpl_dft_mod
USE bsplines
IMPLICIT NONE
!
CONTAINS
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE dismat(spl, mat, alpha, beta)
!
USE bsplines
TYPE(spline1d) :: spl
DOUBLE COMPLEX :: mat(:)
DOUBLE COMPLEX, INTENT(in) :: alpha, beta
!
INTEGER :: dim, nx, nidbas, ngauss, intv, igauss
DOUBLE COMPLEX, ALLOCATABLE :: ft_fun(:,:)
DOUBLE PRECISION, ALLOCATABLE :: xgauss(:), wgauss(:)
!===========================================================================
! 1.0 Prologue
!
! Properties of spline space
!
CALL get_dim(spl, dim, nx, nidbas)
CALL get_gauss(spl, ngauss)
ALLOCATE(xgauss(ngauss), wgauss(ngauss))
ALLOCATE(ft_fun(0:nx-1,2)) ! Up to first derivative
!
! Weak form: integration on first interval
!
intv = 1
CALL get_gauss(spl, ngauss, intv, xgauss, wgauss)
mat = 0.0d0
DO igauss=1,ngauss
CALL ft_basfun(xgauss(igauss), spl, ft_fun, intv)
mat(:) = mat(:) + wgauss(igauss) * ( &
& alpha*ft_fun(:,2)*CONJG(ft_fun(:,2)) &
& - beta*ft_fun(:,1)*CONJG(ft_fun(:,1)) &
& )
END DO
!===========================================================================
! 9.0 Epilogue
!
DEALLOCATE(ft_fun)
DEALLOCATE(xgauss, wgauss)
!
END SUBROUTINE dismat
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
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
END MODULE pde1dp_cmpl_dft_mod
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
PROGRAM main
!
! 1D complex PDE with periodic BC, using DFT
!
USE pde1dp_cmpl_dft_mod
USE bsplines
USE matrix
USE futils
USE fft
!
IMPLICIT NONE
TYPE(spline1d) :: splx
DOUBLE COMPLEX, ALLOCATABLE :: mat(:)
INTEGER ::nrank
!
CHARACTER(len=128) :: file='pde1dp_cmpl_dft.h5'
INTEGER :: fid
INTEGER :: nx, nidbas, ngauss, mmode, npt, dim
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 :: sol_shifted, rhs_shifted
DOUBLE COMPLEX, DIMENSION(:), ALLOCATABLE :: solcal, solana
DOUBLE PRECISION, ALLOCATABLE :: xgauss(:), wgauss(:)
DOUBLE PRECISION :: err_norm
INTEGER :: i
INTEGER :: k, kmin, kmax
!
NAMELIST /newrun/ nx, nidbas, ngauss, alpha, beta, mmode, npt
!===========================================================================
! 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
mmode = 1 ! Fourier mode
alpha = (1.0, 1.0) ! Complex "diffusion"
beta = 1.0
npt = 100
!
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) /)
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.)
nrank = nx ! Rank of the FE matrix
CALL get_dim(splx, dim)
WRITE(*,'(/a,4i6)') 'nrank, dim', nrank, dim
!
! Init DFT
kmin = -nx/2
kmax = nx/2-1
CALL init_dft(splx, kmin, kmax)
!
! FE matrix assembly in Fourier space
!
ALLOCATE(mat(0:nx-1))
CALL dismat(splx, mat, alpha, beta)
CALL putarr(fid, '/mat', mat)
!
! RHS assembly in real space
!
ALLOCATE(rhs(nrank), sol(nrank))
CALL disrhs(splx, rhs, mmode, alpha, beta)
!===========================================================================
! 3.0 Solve the dicretized system
!
! Fourier solve
!
CALL putarr(fid, '/rhs', rhs)
!
CALL fourcol(rhs, 1)
!
ALLOCATE(rhs_shifted(kmin:kmax))
ALLOCATE(sol_shifted(kmin:kmax))
DO k=kmin,kmax
rhs_shifted(k) = rhs(MODULO(k+nx,nx)+1)/REAL(nx,8)
END DO
sol_shifted = rhs_shifted / mat
DO k=kmin,kmax
sol(MODULO(k+nx,nx)+1) = sol_shifted(k)
END DO
!
CALL putarr(fid, '/rhs_fft', rhs)
CALL putarr(fid, '/sol_fft', sol)
!
! Solution in real space
!
CALL fourcol(sol,-1)
CALL putarr(fid, '/sol', sol)
!===========================================================================
! 4.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
!!
!===========================================================================
! 9.0 Clean up
!
DEALLOCATE(xgrid)
DEALLOCATE(mat)
CALL destroy_sp(splx)
CALL closef(fid)
END PROGRAM main
c4science · Help