Homec4science
Diffusion SPClibs (feature/fft-split) Browse

spclibs/examples/tbasfun.f90
03540012fedefeature/fft-split

Last Change
Raw File

tbasfun.f90
View Options

!>
!> @file tbasfun.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
!
! Test scalar and vector versions of def_basfun
!
USE bsplines
IMPLICIT NONE
INTEGER :: nx, nidbas, nrank, npt, jdermx
DOUBLE PRECISION :: dx
DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: xgrid
DOUBLE PRECISION, ALLOCATABLE :: xpt(:), fun(:, :), vfun(:,:,:)
DOUBLE PRECISION, ALLOCATABLE :: fun1(:, :), vfun1(:,:,:)
DOUBLE PRECISION :: errfun
INTEGER :: left, i, nerrs, k
INTEGER, ALLOCATABLE :: vleft(:)
LOGICAL :: nlper=.FALSE.
TYPE(spline1d) :: splx
!
NAMELIST /newrun/ nx, nidbas, npt, jdermx, nlper
!
!===============================================================================
!
! 1D grid
!
nx = 10
nidbas = 3
npt = 1000000
jdermx = 0
READ(*,newrun)
WRITE(*,newrun)
ALLOCATE(xgrid(0:nx))
dx = 1.0d0/REAL(nx)
xgrid = (/ (i*dx,i=0,nx) /)
WRITE(*,'(a/(10f8.3))') 'xgrid', xgrid
!
! Set up spline
!
CALL set_spline(nidbas, 4, xgrid, splx, period=nlper)
nrank = splx%dim
WRITE(*,'(a, i5)') 'nrank =', nrank
WRITE(*,'(a/(10f8.3))') 'knots', splx%knots
!
IF(nx.LE.2) THEN
WRITE(*,'(/a)') 'VAL0'
DO i=1,nx
WRITE(*,'(a,i3)') 'Interval', i
DO k=1,nidbas+1 ! Spline number
WRITE(*,'(10f12.4)') splx%val0(:,k,i)
END DO
END DO
IF(nlper) THEN
WRITE(*,'(/a)') 'VALC'
DO k=1,nidbas+1 ! Spline number
WRITE(*,'(10f12.4)') splx%valc(:,k)
END DO
END IF
END IF
!
ALLOCATE(xpt(npt))
ALLOCATE(vleft(npt))
ALLOCATE(fun(0:nidbas,0:jdermx)) ! Values and derivatives of all Splines
ALLOCATE(vfun(0:nidbas,0:jdermx,npt))
ALLOCATE(fun1(0:nidbas,0:jdermx)) ! Values and derivatives of all Splines
ALLOCATE(vfun1(0:nidbas,0:jdermx,npt))
CALL RANDOM_NUMBER(xpt)
!===============================================================================
!
! Check def_basfun
!
CALL def_basfun(xpt, splx, vfun, vleft)
!
WRITE(*,'(/a)') 'vector def_basfun versus scalar def_basfun'
WRITE(*,'(a6,a12, 2a6, a12)') 'i', 'x', 'left', 'vleft', 'Max. err'
DO i=1,npt
CALL def_basfun(xpt(i), splx, fun, left)
errfun= MAXVAL(ABS(fun(:,:)-vfun(:,:,i)))
WRITE(*,'(i6,1pe12.4,2i6,1pe12.4)') i, xpt(i), left, vleft(i), errfun
END DO
!
IF(npt.LE.10) THEN
WRITE(*,'(/a)') 'Scalar/vector basfun'
DO i=1,npt
CALL basfun(xpt(i), splx, fun, vleft(i)+1)
WRITE(*,'(a,1pe12.4/10(1pe12.4))') 'x = ', xpt(i), fun(:,:)
WRITE(*,'(10(1pe12.4))') vfun(:,:,i)
END DO
END IF
!===============================================================================
!
! Check basfun
!
CALL basfun(xpt, splx, vfun1, vleft+1)
WRITE(*,'(/a,1pe12.4)') 'vector basfun versus vector def_basun: Max err', &
& MAXVAL(ABS(vfun-vfun1))
!
WRITE(*,'(/a)') 'vector basfun versus scalar basfun'
WRITE(*,'(a6,a12,a12)') 'i', 'x', 'Max. err'
DO i=1,npt
CALL basfun(xpt(i), splx, fun1, vleft(i)+1)
errfun= MAXVAL(ABS(fun1(:,:)-vfun1(:,:,i)))
WRITE(*,'(i6,1pe12.4,1pe12.4)') i, xpt(i), errfun
END DO
!!===============================================================================
!
! Clean up
!
CALL destroy_sp(splx)
DEALLOCATE(xgrid)
DEALLOCATE(xpt)
DEALLOCATE(vleft)
DEALLOCATE(fun)
END PROGRAM main
c4science · Help