Homec4science
Diffusion SPClibs (feature/refactoring) Browse

spclibs/examples/fit2dbc_y.f90
03540012fedefeature/refactoring

Last Change
Raw File

fit2dbc_y.f90
View Options

!>
!> @file fit2dbc_y.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
! BC using derivatives at both ends, in the non-periodic direction.
!
! Testing BC on derivative along second direction
!
USE bsplines
USE futils
!
IMPLICIT NONE
CHARACTER(len=128) :: file='fit2d.h5'
INTEGER :: fid
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(npty,nptx) :: fcalc, fexact, errs
TYPE(spline2d) :: splxy
DOUBLE PRECISION :: mem
INTEGER :: i, j, ii
INTEGER :: ibc(2,10)
DOUBLE PRECISION, ALLOCATABLE:: fbc(:,:,:)
!
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:ny,0:nx))
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 = TRANSPOSE(func(xgrid,ygrid))
!
WRITE(*,'(a/(12f8.3))') 'XGRID', xgrid(0:nx)
WRITE(*,'(a/(12f8.3))') 'YGRID', ygrid(0:ny)
IF( nx.LE.10 .AND. ny.LE.10 ) THEN
WRITE(*,'(a)') 'FGRID(y, x)'
DO j=0,ny
WRITE(*,'(12f8.3)') fgrid(j,:)
END DO
WRITE(*,'(a,f8.3)') 'Memory used so far (MB) =', mem()
END IF
!
! Create hdf5 file
!
CALL creatf(file, fid, 'PDE2D Result File')
CALL attach(fid, '/', 'NX', nx)
CALL attach(fid, '/', 'NY', ny)
CALL attach(fid, '/', 'NIDBAS1', nidbas(1))
CALL attach(fid, '/', 'NIDBAS2', nidbas(2))
CALL attach(fid, '/', 'MBESS', mbes)
!===========================================================================
! 2.0 Spline interpolation
!
! Setup the spline interpolation
ii=1 ! Start with first derivative
DO i = 1, nidbas(1)/2
ibc(1,i) = ii+i-1
ibc(2,i) = ii+i-1
END DO
CALL set_splcoef((/nidbas(2), nidbas(1)/), ygrid, xgrid, splxy, (/.TRUE., .FALSE./),&
& ibc2=ibc)
WRITE(*,'(a,f8.3)') 'Memory used so far (MB) =', mem()
!
! Compute spline interpolation coefficients
CALL get_dim(splxy, dims)
ALLOCATE(bcoefs(dims(1),dims(2)))
WRITE(*,'(a,2i4)') 'Dims of spline', dims
!
ALLOCATE(fbc(2, nidbas(1)/2, 0:ny))
fbc=0.0d0
!
!!$! Exact first derivatives at boundaries
!!$ fbc(1,1:1,:) = func1(xgrid(0:0), ygrid(0:ny))
!!$ fbc(2,1:1,:) = func1(xgrid(nx:nx), ygrid(0:ny))
!
!!$! Derivatives at boundaries approximated with FD
!!$ DO j=0,ny
!!$ fbc(1,1,j+1) = fgrid(1,j)-fgrid(0,j)
!!$ fbc(2,1,j+1) = fgrid(nx,j)-fgrid(nx-1,j)
!!$ END DO
!!$ fbc(1,1,:) = fbc(1,1,:)/(xgrid(1)-xgrid(0))
!!$ fbc(2,1,:) = fbc(2,1,:)/(xgrid(nx)-xgrid(nx-1))
!
WRITE(*,'(a/(10f8.3))') 'fbc(1)', fbc(1,1,:)
WRITE(*,'(a/(10f8.3))') 'fbc(2)', fbc(2,1,:)
!
CALL get_splcoef(splxy, fgrid, bcoefs, fbc2=fbc)
!
DEALLOCATE(fbc)
!===========================================================================
! 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 = TRANSPOSE(func(xpt,ypt))
CALL gridval(splxy, ypt, xpt, fcalc, (/0,0/), bcoefs)
errs = fcalc-fexact
WRITE(*,'(a,2(1pe12.3))') 'Min max or errors', MINVAL(errs), MAXVAL(errs)
!
CALL putarr(fid, '/xpt', xpt, 'r')
CALL putarr(fid, '/ypt', ypt, '\theta')
CALL putarr(fid, '/fcalc', fcalc, 'Interpolated')
CALL putarr(fid, '/fexact', fexact, 'Exact')
CALL putarr(fid, '/errs', errs, 'Errors')
!===========================================================================
! 9.0 Epilogue
!
DEALLOCATE(bcoefs)
DEALLOCATE(xgrid, ygrid, fgrid)
CALL destroy_sp(splxy)
CALL closef(fid)
!
!===========================================================================
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), SIZE(y))
DOUBLE PRECISION :: zy
INTEGER :: j
DO j=1,SIZE(y)
zy = -mbes * SIN(mbes*y(j))
func1(:,j) = (mbes - (mbes+2.0d0)*x(:)**2) * x(:)**(mbes-1) * 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