Homec4science
Diffusion SPClibs (feature/fft-split) Browse

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

Last Change
Raw File

fit2d_cmpl.f90
View Options

!>
!> @file fit2d_cmpl.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
!
! Fit 2d grid value function to a 2d spline of any order
!
USE bsplines
!
IMPLICIT NONE
INTEGER :: nx, ny, nidbas(2), mbes, dims(2)
INTEGER, PARAMETER :: npt=10000
DOUBLE PRECISION :: pi, coefx(5), coefy(5)
DOUBLE PRECISION, ALLOCATABLE :: xgrid(:),ygrid(:)
DOUBLE COMPLEX, ALLOCATABLE :: fgrid(:,:), fgrid_calc(:,:), bcoefs(:,:)
DOUBLE PRECISION :: dx, dy, xpt(npt), ypt(npt), errs(npt)
DOUBLE COMPLEX, DIMENSION(npt) :: fcalc, fexact
TYPE(spline2d) :: splxy
INTEGER :: i, j
!
NAMELIST /newrun/ nx, ny, nidbas, mbes, coefx, coefy
!===========================================================================
! 1.0 Prologue
!
! Read in data specific to run
!
nx = 8 ! Number of intevals in x
ny = 8 ! Number of intevals in y
nidbas = (/3,3/) ! Degree of splines
mbes = 2
coefx(1:5) = (/1.0d0, 0.d0, 0.d0, 0.d0, 1.d0/) ! Mesh point distribution function
coefy(1:5) = (/1.0d0, 0.d0, 0.d0, 0.d0, 1.d0/) ! Mesh point distribution function
!
READ(*,newrun)
WRITE(*,newrun)
!
! Define grid on x axis
!
pi = 4.0d0*ATAN(1.0d0)
ALLOCATE(xgrid(0:nx), ygrid(0:ny))
ALLOCATE(fgrid(0:nx,0:ny), fgrid_calc(0:nx,0:ny))
xgrid(0) = 0.0d0; xgrid(nx) = 1.0d0
CALL meshdist(coefx, xgrid, nx)
ygrid(0) = 0.0d0; ygrid(ny) = 2.d0*pi
CALL meshdist(coefy, ygrid, ny)
fgrid = func2(xgrid,ygrid)
!===========================================================================
! 2.0 Spline interpolation
!
! Setup the spline interpolation
CALL set_splcoef(nidbas, xgrid, ygrid, splxy, (/.FALSE., .TRUE./))
!
! Compute spline interpolation coefficients
CALL get_dim(splxy, dims)
ALLOCATE(bcoefs(dims(1),dims(2)))
WRITE(*,'(a,2i4)') 'Dims of spline', dims
!
CALL get_splcoef(splxy, fgrid, bcoefs)
!===========================================================================
! 2.0 Check interpolation
!
CALL RANDOM_NUMBER(xpt)
CALL RANDOM_NUMBER(ypt)
ypt(:) = ypt(:)*2.0*pi
fexact = func1(xpt,ypt)
!
CALL gridval(splxy, xpt, ypt, fcalc, (/0,0/), bcoefs)
errs = ABS(fcalc-fexact)
WRITE(*,'(a,2(1pe12.3))') 'Max errors (on random points)', MAXVAL(errs)
!
CALL gridval(splxy, xgrid, ygrid, fgrid_calc, (/0,0/))
WRITE(*,'(a,2(1pe12.3))') 'Max errors (on grid points)', &
& MAXVAL(ABS(fgrid_calc-fgrid))
!
fgrid_calc = 0.0
DO j=0,ny
ypt(1:nx+1) = ygrid(j)
CALL gridval(splxy, xgrid, ypt(1:nx+1), fgrid_calc(:,j), (/0,0/))
END DO
WRITE(*,'(a,2(1pe12.3))') 'Max errors (on grid points)', &
& MAXVAL(ABS(fgrid_calc-fgrid))
!===========================================================================
! 9.0 Epilogue
!
DEALLOCATE(bcoefs)
DEALLOCATE(xgrid, ygrid, fgrid, fgrid_calc)
CALL destroy_sp(splxy)
!
CONTAINS
FUNCTION func2(x,y)
DOUBLE PRECISION, INTENT(in) :: x(:), y(:)
DOUBLE COMPLEX :: func2(SIZE(x), SIZE(y))
DOUBLE COMPLEX :: zy
INTEGER :: j
DO j=1,SIZE(y)
zy = -mbes * CMPLX(SIN(mbes*y(j)), COS(mbes*y(j)))
func2(:,j) =(1-x(:)**2) * x(:)**mbes * zy
END DO
END FUNCTION func2
FUNCTION func1(x,y)
DOUBLE PRECISION, INTENT(in) :: x(:), y(:)
DOUBLE COMPLEX :: func1(SIZE(x))
DOUBLE COMPLEX :: zy
INTEGER :: j
DO j=1,SIZE(x)
zy = -mbes * CMPLX(SIN(mbes*y(j)), COS(mbes*y(j)))
func1(j) =(1-x(j)**2) * x(j)**mbes * zy
END DO
END FUNCTION func1
END PROGRAM main
c4science · Help