Homec4science
Diffusion SPClibs (master) Browse

spclibs/examples/optim2.f90
93a1e9712d04master

Last Change
Raw File

optim2.f90
View Options

!>
!> @file optim2.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
!
! Test and compare performance of using "spline" and
! "pp" forms. 2D case
!
USE bsplines
!
IMPLICIT NONE
INTEGER :: nx, ny, ngauss(2), nidbas(2)
INTEGER :: npt, d1, d2
INTEGER :: i, j
DOUBLE PRECISION :: pi, dx, dy
DOUBLE PRECISION :: seconds, t0, t1
DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: xgrid, ygrid
DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: coefs
DOUBLE PRECISION, ALLOCATABLE :: xpt(:), ypt(:), fpt00(:), fpt01(:), fpt10(:)
TYPE(spline2d) :: splxy
!
NAMELIST /newrun/ nx, ny, nidbas, npt
!===============================================================================
!
! 2D grid
!
nx = 8
ny = 8
nidbas = (/ 3, 3 /)
npt = 1000000
READ(*,newrun)
WRITE(*,newrun)
!
pi = 4.0d0*ATAN(1.0d0)
ALLOCATE(xgrid(0:nx), ygrid(0:ny))
dx = 1.0d0/REAL(nx)
xgrid = (/ (i*dx,i=0,nx) /)
dy = 2.0d0*pi/REAL(ny)
ygrid = (/ (j*dy,j=0,ny) /)
!
WRITE(*,'(a/(10f8.3))') 'XGRID', xgrid(0:nx)
WRITE(*,'(a/(10f8.3))') 'YGRID', ygrid(0:ny)
!
! Set up spline
!
ngauss = 4
CALL set_spline(nidbas, ngauss, xgrid, ygrid, splxy, (/.FALSE., .TRUE./))
d1 = splxy%sp1%dim
d2 = splxy%sp2%dim
WRITE(*,'(a,3i4)') 'd1, d2 =', d1, d2
WRITE(*,'(/a,/(12(f8.3)))') 'KNOTS in X', splxy%sp1%knots
WRITE(*,'(/a,/(12(f8.3)))') 'KNOTS in Y', splxy%sp2%knots
!===============================================================================
!
! Check and performance of GRIDVAL using PPFORM and SPLINE expansion
!
ALLOCATE(xpt(npt), ypt(npt))
CALL RANDOM_NUMBER(xpt)
CALL RANDOM_NUMBER(ypt); ypt = 2.0d0*pi*ypt
!
ALLOCATE(coefs(d1,d2))
ALLOCATE(fpt00(npt), fpt01(npt), fpt10(npt))
!
coefs = 1.0d0
!
splxy%sp1%nlppform = .TRUE.
splxy%sp2%nlppform = .TRUE.
CALL gridval(splxy, xpt, ypt, fpt00, (/0,0/), coefs)
!
t0 = seconds()
CALL gridval(splxy, xpt, ypt, fpt00, (/0,0/))
CALL gridval(splxy, xpt, ypt, fpt01, (/0,1/))
CALL gridval(splxy, xpt, ypt, fpt10, (/1,0/))
t1 = seconds()-t0
WRITE(*,'(/a,3(1pe12.3))') 'PPFORM: Max errors', &
& MAXVAL(ABS(fpt00-1.0d0)), MAXVAL(ABS(fpt01)), MAXVAL(ABS(fpt10))
WRITE(*,'(a,3(1pe12.3))') 'time(s)', t1/REAL(npt)
!
splxy%sp1%nlppform = .FALSE.
splxy%sp2%nlppform = .FALSE.
CALL gridval(splxy, xpt, ypt, fpt00, (/0,0/), coefs)
t0 = seconds()
CALL gridval(splxy, xpt, ypt, fpt00, (/0,0/))
CALL gridval(splxy, xpt, ypt, fpt01, (/0,1/))
CALL gridval(splxy, xpt, ypt, fpt10, (/1,0/))
t1 = seconds()-t0
WRITE(*,'(/a,3(1pe12.3))') 'SPLINES: Max errors', &
& MAXVAL(ABS(fpt00-1.0d0)), MAXVAL(ABS(fpt01)), MAXVAL(ABS(fpt10))
WRITE(*,'(a,3(1pe12.3))') 'time(s)', t1/REAL(npt)
!===============================================================================
!
! Clean up
!
CALL destroy_sp(splxy)
DEALLOCATE(xgrid, ygrid, coefs)
DEALLOCATE(xpt, ypt, fpt00, fpt01, fpt10)
END PROGRAM main
c4science · Help