Homec4science
Diffusion SPClibs (feature/refactoring) Browse

spclibs/examples/fit1d.f90
3324689e437cfeature/refactoring

Last Change
Raw File

fit1d.f90
View Options

!>
!> @file fit1d.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
!
! Fit grid value function to a spline of any order
!
USE bsplines
USE futils
!
IMPLICIT NONE
INTEGER :: nx, nidbas
DOUBLE PRECISION :: a, b, coefx(5)
!!$ DOUBLE PRECISION, ALLOCATABLE :: xgrid(:), fgrid(:), coefs(:)
DOUBLE PRECISION, ALLOCATABLE :: xgrid(:), fgrid(:,:), coefs(:,:)
INTEGER :: i, dim, left
TYPE(spline1d) :: spl
DOUBLE PRECISION :: dx
INTEGER :: npts
DOUBLE PRECISION, ALLOCATABLE :: xpt(:), fcalc(:), fexact(:), err(:)
DOUBLE PRECISION, ALLOCATABLE :: fun(:,:)
DOUBLE PRECISION, POINTER :: splines(:,:) => null()
!
CHARACTER(len=128) :: file='fit1d.h5'
INTEGER :: fid
!
NAMELIST /newrun/ nx, nidbas, a, b, coefx
!===========================================================================
! 1.0 Prologue
!
! Read in data specific to run
!
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 and function values
!
ALLOCATE(xgrid(0:nx), fgrid(0:nx,1))
xgrid(0) = a
xgrid(nx) = b
CALL meshdist(coefx, xgrid, nx)
fgrid(:,1) = func(xgrid)
WRITE(*,'(a/(10f8.3))') 'xgrid', xgrid
WRITE(*,'(a/(10f8.3))') 'fgrid', fgrid
!
! Create hdf5 file
!
CALL creatf(file, fid, 'FIT1D Result File')
CALL attach(fid, '/', 'NX', nx)
CALL attach(fid, '/', 'NIDBAS', nidbas)
!===========================================================================
! 2.0 Spline interpolation
!
! Set up the spline interpolation
!!$ CALL splcoef_setup(nidbas, xgrid, spl)
CALL set_splcoef(nidbas, xgrid, spl)
PRINT*, 'nlequid =', spl%nlequid
!
! Compute spline values and derivatives at Boundaries
ALLOCATE(fun(nidbas+1,0:nidbas))
WRITE(*,'(/a)') 'spline at the left boundary'
CALL locintv(spl, a, left)
CALL basfun(a, spl, fun, left+1)
DO i=0,nidbas
WRITE(*,'(8(1pe12.4))') fun(:,i)
END DO
!
WRITE(*,'(/a)') 'spline at the right boundary'
CALL locintv(spl, b, left)
CALL basfun(b, spl, fun, left+1)
DO i=0,nidbas
WRITE(*,'(8(1pe12.4))') fun(:,i)
END DO
DEALLOCATE(fun)
!
CALL get_dim(spl, dim)
ALLOCATE(coefs(dim,1))
WRITE(*,'(2(a,i5, 2x))') 'NX =', nx, 'DIM =', dim
WRITE(*,'(/a,i6,a,i6/(10(f8.3)))') 'KNOTS of splines', LBOUND(spl%knots), &
& ':',UBOUND(spl%knots), spl%knots
!
! From given grid values fgrid, compute the spline coefs
CALL get_splcoef(spl, fgrid, coefs)
WRITE(*,'(a/(10f8.3))') 'coefs', coefs
!
! Plot all splines
npts = 100
ALLOCATE(xpt(npts))
dx = (b-a)/REAL(npts-1)
DO i=1,npts
xpt(i) = a + (i-1)*dx
END DO
CALL allsplines(spl, xpt, splines)
CALL putarr(fid, '/X', xpt)
CALL putarr(fid,'/SPLINES', splines)
!
! Check interpolation
ALLOCATE(fcalc(npts), fexact(npts), err(npts))
fexact = func(xpt)
!
! Function values
CALL gridval(spl, xpt, fcalc, 0, coefs(:,1))
err = fexact - fcalc
CALL putarr(fid, '/FEXACT', fexact, 'Exact values')
CALL putarr(fid, '/FCALC', fcalc, 'Interpolated values')
CALL putarr(fid, '/ERROR', err, 'Interpolation errors')
WRITE(*,'(/a,1pe12.3)') 'Norm of error', norm2(err)/norm2(fexact)
!
! Derivatives values
CALL gridval(spl, xpt, fcalc, 1)
fexact = func1(xpt)
err = fexact - fcalc
CALL putarr(fid, '/FEXACT1', fexact, 'Exact values of first derivative')
CALL putarr(fid, '/FCALC1', fcalc, 'Interpolated first derivative')
CALL putarr(fid, '/ERROR1', err, 'Interpolation errors on first derivative')
WRITE(*,'(/a,1pe12.3)') 'Norm of error', norm2(err)/norm2(fexact)
!
!===========================================================================
! 9.0 Epilogue
!
DEALLOCATE(xpt, splines, fexact, fcalc, err)
DEALLOCATE(xgrid, fgrid)
CALL destroy_sp(spl)
CALL closef(fid)
!===========================================================================
!
CONTAINS
FUNCTION func(x)
DOUBLE PRECISION, INTENT(in) :: x(:)
DOUBLE PRECISION :: func(SIZE(x))
!!$ INTEGER :: n
!!$ n = SIZE(x)
!!$ func = 1.d0+x*(1.d0+x*(1.d0+x))
!!$ func(1:n/2) = 1.0d0
!!$ func(n/2+1:n) = 0.5d0
func = EXP(-8.*x*x)
END FUNCTION func
!
FUNCTION func1(x)
DOUBLE PRECISION, INTENT(in) :: x(:)
DOUBLE PRECISION :: func1(SIZE(x))
!!$ INTEGER :: n
!!$ n = SIZE(x)
!!$ func = 1.d0+x*(1.d0+x*(1.d0+x))
!!$ func(1:n/2) = 1.0d0
!!$ func(n/2+1:n) = 0.5d0
func1 = -16.d0*x*EXP(-8.*x*x)
END FUNCTION func1
!
FUNCTION norm2(x)
!
! Compute the 2-norm of array x
!
DOUBLE PRECISION :: norm2
DOUBLE PRECISION, INTENT(in) :: x(:)
DOUBLE PRECISION :: sum2
INTEGER :: i
!
sum2 = 0.0d0
DO i=1,SIZE(x,1)
sum2 = sum2 + x(i)**2
END DO
norm2 = SQRT(sum2)
END FUNCTION norm2
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