Homec4science
Diffusion SPClibs (feature/refactoring) Browse

spclibs/examples/pde1dp.f90
3324689e437cfeature/refactoring

Last Change
Raw File

pde1dp.f90
View Options

!>
!> @file pde1dp.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
!
! 1D PDE with priodic BC
!
USE pde1dp_mod
USE bsplines
USE matrix
USE futils
USE conmat_mod
!
IMPLICIT NONE
CHARACTER(len=128) :: file='pde1dp.h5'
INTEGER :: fid
INTEGER :: nx, nidbas, ngauss, ibcoef
INTEGER :: nrank, kl, ku, dim
DOUBLE PRECISION :: coefs(5)
DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: xgrid, sol, rhs
INTEGER, PARAMETER :: npts=100
DOUBLE PRECISION, DIMENSION(0:npts-1) :: xpts, frhs
DOUBLE PRECISION :: dx, errmx
DOUBLE PRECISION, ALLOCATABLE :: arr(:,:)
TYPE(periodic_mat) :: mat
INTEGER :: i, j
!
NAMELIST /newrun/ nx, nidbas, ngauss, ibcoef, coefs
!===========================================================================
! 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
ibcoef = 1 ! Index of non-zero spline coef for RHS
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)
!
! Create hdf5 file
!
CALL creatf(file, fid, 'PDE1DP Result File')
CALL attach(fid, '/', 'NX', nx)
CALL attach(fid, '/', 'NIDBAS', nidbas)
CALL attach(fid, '/', 'NGAUSS', ngauss)
!===========================================================================
! 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
!
! Mass matrix assembly
!
kl = nidbas
ku = kl
CALL init(kl, ku, nrank, 1, mat)
CALL get_dim(splx, dim)
WRITE(*,'(/a,4i6)') 'kl, ku, nrank, dim', kl, ku, nrank, dim
!!$ CALL dismat(splx, mat)
CALL conmat(splx, mat, coefeq_mass)
!
! Store matrix in hdf5 file
!
ALLOCATE(arr(nrank,nrank))
DO j=1,nrank
CALL getcol(mat, j, arr(:,j))
END DO
CALL putarr(fid, '/mata', arr)
DEALLOCATE(arr)
!
! Check RHS constructed using input spline coefs.
!
ALLOCATE(bcoef(0:dim-1))
bcoef = 0.0d0; bcoef(ibcoef-1) = 1.0d0
!
DO i=nrank,dim-1 ! Periodicity to fill array of spline coefs
bcoef(i) = bcoef(MODULO(i,nrank))
END DO
WRITE(*,'(/a/(10f8.3))') 'bcoef from input', bcoef
dx = (xgrid(nx)-xgrid(0))/npts
DO i=0,npts-1
xpts(i) = xgrid(0) + i*dx
frhs(i) = rhseq(xpts(i))
END DO
CALL creatg(fid, '/rhs')
CALL putarr(fid,'/rhs/x', xpts)
CALL putarr(fid,'/rhs/f', frhs)
!
! Assembly RHS and check A*x = f, using method vmx
!
ALLOCATE(rhs(nrank), sol(nrank))
CALL disrhs(splx, rhs)
sol = vmx(mat, bcoef(0:nrank-1))
WRITE(*,'(/6x,3a12)') 'A*x', 'rhs', 'Err'
errmx = 0.0d0
DO i=1,nrank
WRITE(*,'(i6,3(1pe12.3))') i, sol(i), rhs(i), sol(i)-rhs(i)
errmx=MAX(errmx,ABS(sol(i)-rhs(i)))
END DO
WRITE(*,'(a,1pe12.3)') 'Max. error =', errmx
!
! Factor and solve
!
CALL factor(mat)
CALL bsolve(mat, rhs, sol)
WRITE(*,'(/6x,3a12)') 'Computed', 'Exact', 'Err'
errmx = 0.0d0
DO i=1,nrank
WRITE(*,'(i6,3(1pe12.3))') i, sol(i), bcoef(i-1), sol(i)-bcoef(i-1)
errmx=MAX(errmx,ABS(sol(i)-bcoef(i-1)))
END DO
WRITE(*,'(a,1pe12.3)') 'Max. error =', errmx
!===========================================================================
! 9.0 Clean up
!
DEALLOCATE(xgrid)
DEALLOCATE(bcoef)
DEALLOCATE(rhs, sol)
CALL destroy(mat)
CALL destroy_sp(splx)
CALL closef(fid)
CONTAINS
SUBROUTINE coefeq_mass(x, idt, idw, c)
DOUBLE PRECISION, INTENT(in) :: x
INTEGER, INTENT(out) :: idt(:), idw(:)
DOUBLE PRECISION, INTENT(out) :: c(:)
!
! Mass matrix
!
c(1) = 1.0d0
idt(1) = 0
idw(1) = 0
END SUBROUTINE coefeq_mass
END PROGRAM main
c4science · Help