Homec4science
Diffusion SPClibs (feature/refactoring) Browse

spclibs/examples/gyro.f90
03540012fedefeature/refactoring

Last Change
Raw File

gyro.f90
View Options

!>
!> @file gyro.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
!
! Gyro-average using splines.
! F(x)=cos(x) => \bar F(x,rho) = J_0(rho) cos(x)
!
USE bsplines
USE futils
!
IMPLICIT NONE
INTEGER, PARAMETER :: nx=10, nidbas=3, dim=nx+nidbas, npt=100, &
& nrho=21, nnq=5
DOUBLE PRECISION :: xgrid(0:nx), fgrid(0:nx), coefs(dim)
DOUBLE PRECISION :: xpt(npt), fcalc(npt), fexact(npt)
DOUBLE PRECISION :: averf(0:nx), averfexact(0:nx)
DOUBLE PRECISION, POINTER :: splines(:,:)
TYPE(spline1d) :: spl
DOUBLE PRECISION :: pi, twopi, dx, lperiod, dth
DOUBLE PRECISION :: drho, rho(nrho), erraver(nrho,nnq)
INTEGER :: i, j, nq(nnq)
DOUBLE PRECISION :: dbesj0
!
CHARACTER(len=128) :: file='gyro.h5'
INTEGER :: fid
!
INTERFACE
SUBROUTINE gyro(spl, xgrid, coefs, rho, nq, averf)
USE bsplines
TYPE(spline1d) :: spl
DOUBLE PRECISION, INTENT(in) :: xgrid(0:), coefs(:), rho
INTEGER, INTENT(in) :: nq
DOUBLE PRECISION, INTENT(out) :: averf(0:)
END SUBROUTINE gyro
END INTERFACE
!
pi = 4.0d0*ATAN(1.0d0)
twopi = 2.0d0*pi
!
! Create hdf5 file
CALL creatf(file, fid, 'gyro Result File')
CALL attach(fid, '/', 'NX', nx)
CALL attach(fid, '/', 'NIDBAS', nidbas)
!
! Grid and function values
dx = twopi/nx
xgrid(0) = 0.0d0
DO i=1,nx
xgrid(i) = xgrid(0) + i*dx
END DO
lperiod = xgrid(nx)-xgrid(0)
fgrid = func(xgrid)
!
! Spline interpolation
CALL set_splcoef(nidbas, xgrid, spl, period=.TRUE.)
CALL get_splcoef(spl, fgrid, coefs)
WRITE(*,'(a)') 'Spline coefficients'
WRITE(*,'(i5,f12.4)') (i-1, coefs(i), i=1,dim)
!
! Error of interpolation
CALL RANDOM_NUMBER(xpt)
xpt = twopi*xpt
CALL gridval(spl, xpt, fcalc, 0, coefs)
fexact = func(xpt)
!!$ WRITE(*,'(a)') 'Interpolated and exact f'
!!$ WRITE(*,'(3(1pe12.3))') (xpt(i), fcalc(i), fexact(i), i=1,npt)
WRITE(*,'(a,1pe12.3)') 'Interpolation error', norm2(fcalc-fexact)
CALL putarr(fid, '/X', xpt)
CALL putarr(fid, '/FEXACT', fexact, 'Exact values')
CALL putarr(fid, '/FCALC', fcalc, 'Interpolated values')
!
! Gyro-averaged of F at grid points xgrid(0:nx-1)
drho = 5.0d0/nrho
DO i=1,nrho
rho(i) = i*drho
DO j=1,nnq
nq(j) = 2**(j+1)
CALL gyro(spl, xgrid, coefs, rho(i), nq(j), averf)
averfexact = dbesj0(rho(i))*COS(xgrid)
erraver(i,j) = norm2(averfexact-averf)
END DO
END DO
!
WRITE(*,'(a,f8.3,i5)') 'averaged f at rho, nq =', rho(nrho), nq(nnq)
WRITE(*,'(3(1pe12.3))') (xgrid(i),averf(i),averfexact(i), i=0,nx)
CALL putarr(fid, '/XGRID', xgrid)
CALL putarr(fid, '/AVERF', averf, 'Averaged F')
CALL putarr(fid, '/AVERFEXACT', averfexact, 'Averaged F exact')
!
CALL putarr(fid, '/RHO', rho)
CALL putarr(fid, '/NQ', nq)
CALL putarr(fid, '/ERRAVER', erraver)
!
! Clean up
CALL destroy_sp(spl)
CALL closef(fid)
CONTAINS
FUNCTION func(x)
DOUBLE PRECISION, INTENT(in) :: x(:)
DOUBLE PRECISION :: func(SIZE(x))
func = COS(x)
END FUNCTION func
FUNCTION norm2(x)
DOUBLE PRECISION :: norm2
DOUBLE PRECISION, INTENT(in) :: x(:)
norm2 = SQRT(DOT_PRODUCT(x,x))
END FUNCTION norm2
END PROGRAM main
!
SUBROUTINE gyro(spl, xgrid, coefs, rho, nq, averf)
!
! Gyro-average using spline SPL and NQ point quadratire for
! theta-integration.
!
USE bsplines
IMPLICIT NONE
TYPE(spline1d) :: spl
DOUBLE PRECISION, INTENT(in) :: xgrid(0:), coefs(:), rho
INTEGER, INTENT(in) :: nq
DOUBLE PRECISION, INTENT(out) :: averf(0:)
!
DOUBLE PRECISION :: th(nq), wth(nq), xq(nq)
DOUBLE PRECISION, ALLOCATABLE :: avermat(:,:)
DOUBLE PRECISION :: pi, twopi, lperiod, dth
DOUBLE PRECISION, POINTER :: splines(:,:)
INTEGER :: i, j, iq, nx, dim
!
pi = 4.0d0*ATAN(1.0d0)
twopi = 2.0d0*pi
!
! Quadrature in theta
dth = twopi/nq
th(1) = -pi + dth/2.0d0
DO iq=2,nq
th(iq) = th(iq-1)+dth
END DO
wth = dth
!
! Gyro-averaging matrix
CALL get_dim(spl, dim, nx)
lperiod = xgrid(nx)-xgrid(0)
ALLOCATE(avermat(0:nx,dim))
DO i=0,nx
xq = xgrid(i) + rho*COS(th)
xq = xgrid(0) + MODULO(xq-xgrid(0), lperiod)
CALL allsplines(spl, xq, splines)
DO j=1,dim
avermat(i,j) = DOT_PRODUCT(wth, splines(:,j))/twopi
END DO
END DO
!
! Gyro-averaged of F at grid points xgrid(0:nx)
averf = MATMUL(avermat, coefs)
!
DEALLOCATE(avermat)
END SUBROUTINE gyro
c4science · Help