Homec4science
Diffusion SPClibs (feature/fft-split) Browse

spclibs/examples/fit2d1d.f90
93a1e9712d04feature/fft-split

Last Change
Raw File

fit2d1d.f90
View Options

!>
!> @file fit2d1d.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
! Interpolating on an grid (x_i,y_j) or a set of particle
! positions (x_p,y_p).
!
USE bsplines
!
IMPLICIT NONE
INTEGER :: nx, ny, nidbas(2), mbes, dims(2)
INTEGER, PARAMETER :: nptx=100, npty=100
DOUBLE PRECISION :: pi, coefx(5), coefy(5)
DOUBLE PRECISION, ALLOCATABLE :: xgrid(:),ygrid(:),fgrid(:,:),bcoefs(:,:)
DOUBLE PRECISION :: dx, dy, xpt(nptx), ypt(npty)
DOUBLE PRECISION, DIMENSION(nptx,npty) :: fcalc, fexact, errs
DOUBLE PRECISION, DIMENSION(nptx*npty) :: xp, yp, fcalcp, fexactp, errsp
TYPE(spline2d) :: splxy
DOUBLE PRECISION :: mem
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), fgrid(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 = func(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
!
dx=(xgrid(nx)-xgrid(0))/REAL(nptx-1)
dy=(ygrid(ny)-ygrid(0))/REAL(npty-1)
DO i=1,nptx
xpt(i) = xgrid(0) + (i-1)*dx
END DO
DO i=1,npty
ypt(i) = ygrid(0) + (i-1)*dy
END DO
fexact = func(xpt,ypt)
CALL gridval(splxy, xpt, ypt, fcalc, (/0,0/), bcoefs)
errs = fcalc-fexact
WRITE(*,*) 'Using the GRIDVAL2D2D'
WRITE(*,'(a,2(1pe12.3))') 'Min max or errors', MINVAL(errs), MAXVAL(errs)
!
! The 2d1d version
WRITE(*,*) 'Using the GRIDVAL2D1D'
CALL RANDOM_NUMBER(xp)
CALL RANDOM_NUMBER(yp)
yp=2.0*pi*yp
fexactp = func1(xp,yp)
!!$ CALL gridval(splxy, xp, yp, fcalcp, (/0,0/), bcoefs)
CALL gridval(splxy, xp, yp, fcalcp, (/0,0/))
errsp = fcalcp-fexactp
WRITE(*,'(a,2(1pe12.3))') 'Min max or errors', MINVAL(errsp), MAXVAL(errsp)
!===========================================================================
! 9.0 Epilogue
!
DEALLOCATE(bcoefs)
DEALLOCATE(xgrid, ygrid, fgrid)
CALL destroy_sp(splxy)
!
CONTAINS
FUNCTION func(x,y)
DOUBLE PRECISION, INTENT(in) :: x(:), y(:)
DOUBLE PRECISION :: func(SIZE(x), SIZE(y))
DOUBLE PRECISION :: zy
INTEGER :: j
DO j=1,SIZE(y)
zy = -mbes * SIN(mbes*y(j))
func(:,j) =(1-x(:)**2) * x(:)**mbes * zy
END DO
END FUNCTION func
FUNCTION func1(x,y)
DOUBLE PRECISION, INTENT(in) :: x(:), y(:)
DOUBLE PRECISION :: func1(SIZE(x))
DOUBLE PRECISION :: zy
INTEGER :: j
DO j=1,SIZE(x)
zy = -mbes * SIN(mbes*y(j))
func1(j) =(1-x(j)**2) * x(j)**mbes * zy
END DO
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
c4science · Help