Homec4science
Diffusion SPClibs (feature/refactoring) Browse

spclibs/examples/driv2.f90
93a1e9712d04feature/refactoring

Last Change
Raw File

driv2.f90
View Options

!>
!> @file driv2.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
!
! Integration of splines
!
USE bsplines
!
IMPLICIT NONE
TYPE(spline1d) :: spx
INTEGER :: nx, nidbas, ngauss
DOUBLE PRECISION :: a, b, coefx(5)
DOUBLE PRECISION, ALLOCATABLE :: xgrid(:), fun(:,:), finteg(:)
DOUBLE PRECISION, ALLOCATABLE :: xg(:), wg(:)
DOUBLE PRECISION :: support, res, err
INTEGER :: i, ig, j, jj, left
INTEGER :: dim, ng
LOGICAL :: periodic
NAMELIST /newrun/ periodic, nx, nidbas, ngauss, a, b, coefx
!===========================================================================
! 1.0 Set up grids
!
! Read in data specific to run
!
periodic = .FALSE.
nidbas = 3
ngauss = 4
nx = 10 ! Number of intevals in x
a = 0.0d0 ! Left boundary of interval
b = 1.0d0 ! Right boundary of interval
coefx(1:5) = (/1.0d0, 0.d0, 0.d0, 0.d0, 1.d0/) ! Mesh point distribution function
READ(*,newrun)
WRITE(*,newrun)
!
! Define grid/knots
!
ALLOCATE(xgrid(0:nx))
xgrid(0) = a
xgrid(nx) = b
CALL meshdist(coefx, xgrid, nx)
WRITE(*,'(a/(10f8.3))') 'XGRID', xgrid(0:nx)
WRITE(*,'(a/(10f8.3))') 'Diff XGRID', (xgrid(i)-xgrid(i-1), i=1,nx)
WRITE(*,'(a/(10f8.3))') 'Diff XGRID**2', (xgrid(i)**2-xgrid(i-1)**2, i=1,nx)
!===========================================================================
! 2.0 Set up splines
!
CALL set_spline(nidbas, ngauss, xgrid, spx, periodic)
WRITE(*,'(/a,i6,a,i6/(10(f8.3)))') 'KNOTS of x-splines', LBOUND(spx%knots), &
& ':',UBOUND(spx%knots), spx%knots
WRITE(*,'(2(a,i5, 2x))') 'NX =', nx, 'DIM =', spx%dim
!===========================================================================
! 3.0 Integrate all splines
!
CALL get_dim(spx, dim)
ALLOCATE( finteg(0:dim-1), xg(ngauss), wg(ngauss), fun(0:nidbas,1) )
finteg = 0.0
fun = 0.0
DO i=1,nx ! Loop thru the intervals
CALL get_gauss(spx, ng, i, xg, wg)
DO ig=1,ng ! Loop thru Gauss points
CALL basfun(xg(ig), spx, fun, i)
left = i-1
DO j=0,nidbas ! Loop thru the splines in this interval
jj = left+j
IF( periodic ) jj = MODULO(left+j, nx)
finteg(jj) = finteg(jj) + wg(ig)*fun(j,1)
END DO
END DO
END DO
!!$ IF( periodic ) THEN
!!$ DO i=nx,dim-1
!!$ finteg(i) = finteg(i-nx)
!!$ END DO
!!$ END IF
!
WRITE(*,'(a/(10f10.5))') 'Integrals of splines', finteg
PRINT*, 'Sum of finteg', SUM(finteg)
!!$ IF( periodic ) THEN
!!$ PRINT*, 'Sum of finteg', SUM(finteg(0:nx-1))
!!$ ELSE
!!$ PRINT*, 'Sum of finteg', SUM(finteg)
!!$ END IF
!
WRITE(*,'(a/(10f10.5))') 'Integrals of splines from module', spx%intspl
PRINT*, 'Sum of finteg', SUM(spx%intspl)
WRITE(*,'(a5,4a12)') '#', 'I', 'S', '(p+1)I/S', '(p+1)I/S-1'
DO i=0,spx%dim-1
support = spx%knots(i+1)-spx%knots(i-nidbas)
res = spx%intspl(i)/support*(nidbas+1)
err = res - 1.0d0
WRITE(*,'(i5,4(1pe12.4))') i, spx%intspl(i), support, res, err
END DO
!===========================================================================
! 9.0 Epilogue
!
DEALLOCATE( finteg, xg, wg, fun )
DEALLOCATE(xgrid)
CALL destroy_sp(spx)
END PROGRAM main
!+++
SUBROUTINE meshdist(c, x, nx)
!
! Construct an 1d non-equidistant mesh given a
! mesh distribution function.
!
IMPLICIT NONE
DOUBLE PRECISION, INTENT(in) :: c(5)
INTEGER, INTENT(iN) :: nx
DOUBLE PRECISION, INTENT(inout) :: x(0:nx)
INTEGER :: nintg
DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: xint, fint
DOUBLE PRECISION :: a, b, dx, f0, f1, scal
INTEGER :: i, k
!
a=x(0)
b=x(nx)
nintg = 10*nx
ALLOCATE(xint(0:nintg), fint(0:nintg))
!
! Mesh distribution
!
dx = (b-a)/REAL(nintg)
xint(0) = a
fint(0) = 0.0d0
f1 = fdist(xint(0))
DO i=1,nintg
f0 = f1
xint(i) = xint(i-1) + dx
f1 = fdist(xint(i))
fint(i) = fint(i-1) + 0.5*(f0+f1)
END DO
!
! Normalization
!
scal = REAL(nx) / fint(nintg)
fint(0:nintg) = fint(0:nintg) * scal
!!$ WRITE(*,'(a/(10f8.3))') 'FINT', fint
!
! Obtain mesh point by (inverse) interpolation
!
k = 1
DO i=1,nintg-1
IF( fint(i) .GE. REAL(k) ) THEN
x(k) = xint(i) + (xint(i+1)-xint(i))/(fint(i+1)-fint(i)) * &
& (k-fint(i))
k = k+1
END IF
END DO
!
DEALLOCATE(xint, fint)
CONTAINS
DOUBLE PRECISION FUNCTION fdist(x)
DOUBLE PRECISION, INTENT(in) :: x
fdist = c(1) + c(2)*x + c(3)*EXP(-((x-c(4))/c(5))**2)
END FUNCTION fdist
END SUBROUTINE meshdist
!+++
c4science · Help