Homec4science
Diffusion SPClibs (feature/refactoring) Browse

spclibs/examples/optim3.f90
93a1e9712d04feature/refactoring

Last Change
Raw File

optim3.f90
View Options

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