Homec4science
Diffusion SPClibs (feature/fft-split) Browse

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

Last Change
Raw File

adv.f90
View Options

!>
!> @file adv.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
!
! 1d Periodic Advection: F(x,t) = F(x-u*dt,t-dt)
! using module bsplines
!
USE bsplines
!
IMPLICIT NONE
INCLUDE 'fftw_f77.h'
TYPE(spline1d) :: spl
INTEGER, PARAMETER :: nhistmx=1000, ncomb=4
INTEGER :: nx, nidbas, dim
INTEGER :: nstep, nskipt, nhist, mhist
DOUBLE PRECISION :: a, b, dt, u, w, coefx(5)
DOUBLE PRECISION, DIMENSION(0:nhistmx) :: thist, tmass, tfmin, tfmax, ermass
DOUBLE PRECISION, ALLOCATABLE :: xgrid(:), fgrid(:), xshft(:)
DOUBLE PRECISION, ALLOCATABLE :: coefs(:), ferr(:), kx(:), ampl(:)
DOUBLE COMPLEX, ALLOCATABLE :: cfgrid(:), ffft(:)
INTEGER(8) :: forw
DOUBLE PRECISION :: time
INTEGER :: i
NAMELIST /newrun/ nx, nidbas, a, b, dt, u, w, coefx
!===========================================================================
! 1.0 Prologue
!
! Read in data specific to run
!
nx = 100 ! Number of intevals in x
a = 0.0 ! Left boundary of interval
b = 100.0 ! Right boundary of interval
dt = 0.1 ! Time step
u = 1.0 ! Velocity
w = 2.0 ! Shape of initial function
coefx(1:5) = (/1.0d0, 0.d0, 0.d0, 0.d0, 1.d0/) ! Mesh point distribution function
!
READ(*,newrun)
WRITE(*,newrun)
!===========================================================================
! 2.0 Define initial conditions
!
! Set up mesh
!
ALLOCATE(xgrid(0:nx), xshft(0:nx), fgrid(0:nx), ferr(0:nx))
xgrid(0) = a
xgrid(nx) = b
CALL meshdist(coefx, xgrid, nx)
!
! Set up the spline interpolation
!
CALL set_splcoef(nidbas, xgrid, spl, period=.TRUE.)
CALL get_dim(spl, dim)
WRITE(*,'(a,i6)') 'dimension of splines', dim
ALLOCATE(coefs(dim))
!
! Initial conditions
!
time = 0.0d0
nstep = 0
nskipt = 1
DO i=0,nx-1
fgrid(i) = finit(xgrid(i))
END DO
fgrid(nx) = fgrid(0)
ferr = 0.0
CALL get_splcoef(spl, fgrid, coefs)
WRITE(*,'(a/(10f8.3))') 'knots', spl%knots
WRITE(*,'(a/(10f8.3))') 'xgrid', xgrid
WRITE(*,'(a/(10f8.3))') 'fgrid', fgrid
WRITE(*,'(a/(10f8.3))') 'coefs', coefs
!
! Set up FFT
!
ALLOCATE(kx(-nx/2+1:nx/2), ampl(-nx/2+1:nx/2))
ALLOCATE(cfgrid(0:nx-1), ffft(0:nx-1))
DO i=-nx/2+1,nx/2
kx(i) = i
END DO
cfgrid(0:nx-1) = fgrid(0:nx-1)
CALL fftw_f77_create_plan(forw, nx, FFTW_FORWARD, FFTW_ESTIMATE)
CALL fftw_f77_one(forw, cfgrid, ffft)
ampl(0:nx/2) = ABS(ffft(0:nx/2))
ampl(-nx/2+1:-1) = ABS(ffft(nx/2+1:nx-1))
!
! Set up history arrays
!
nhist = 0
thist(nhist) = 0.0
!!$ tmass(nhist) = SUM(coefs(1:nx))
tmass(nhist) = SUM(fgrid(1:nx))
tfmin(nhist) = MINVAL(fgrid)
tfmax(nhist) = MAXVAL(fgrid)
ermass(nhist) = tmass(nhist)-tmass(0)
WRITE(*,'(a,(10f8.3))') 'Initial mass', tmass(nhist)
!
! Initialize Xgrafix
!
CALL xginit(3,'ADV','adv',' ',' ',' ',' ',time)
!
CALL xgset2d('linlin', 'X', 'F', 'open', 220, 60, 1.d0, 1.d0, &
& .FALSE., .FALSE., xgrid(0), xgrid(nx), -0.2d0, 1.2d0)
CALL xgcurve(xgrid, fgrid, nx+1, 1)
!
CALL xgset2d('linlin', 'X', 'FERR', 'open', 620, 60, 1.d0, 1.d0, &
& .FALSE., .FALSE., xgrid(0), xgrid(nx), -1.d0, 1.d0)
CALL xgcurve(xgrid, ferr, nx+1, 1)
!
CALL xgset2d('linlin', 'Time', 'Error Mass', 'open', 820, 400, 1.d0, 1.d0, &
& .TRUE., .TRUE., 0.d0, 1.d0, 0.d0, 1.d0)
CALL xgcurve(thist, ermass, nhist, 1)
!
CALL xgset2d('linlin', 'Time', 'Min/Max', 'open', 420, 400, 1.d0, 1.d0, &
& .TRUE., .TRUE., 0.d0, 1.d0, 0.d0, 1.d0)
CALL xgcurve(thist, tfmin, nhist, 1)
CALL xgcurve(thist, tfmax, nhist, 2)
!
CALL xgset2d('linlin', 'kx', 'Amplitude of F', 'open', 20, 400, 1.d0, 1.d0, &
& .FALSE., .FALSE., kx(-nx/2+1), kx(nx/2), &
& 0.0d0, MAXVAL(ampl))
CALL xgcurve(kx, ampl, nx, 1)
!
CALL xgupdate
!===========================================================================
! 3.0 Time loop
!
nskipt = 1
DO
nstep = nstep+1
time = time+dt
CALL xgevent
!
! Shift x
!
CALL get_splcoef(spl, fgrid, coefs)
xshft(0:nx) = xgrid(0:nx) - u*dt
CALL gridval(spl, xshft, fgrid, 0, coefs)
!
xshft(0:nx) = xgrid(0:nx) - u*time
DO i =0,nx
ferr(i) = fgrid(i) - finit(xshft(i))
END DO
!
cfgrid(0:nx-1) = fgrid(0:nx-1)
CALL fftw_f77_one(forw, cfgrid, ffft)
ampl(0:nx/2) = ABS(ffft(0:nx/2))
ampl(-nx/2+1:-1) = ABS(ffft(nx/2+1:nx-1))
!
! Diagnostics
!
IF( MOD(nstep,nskipt) .EQ. 0 ) THEN
nhist = nhist+1
IF( nhist .GT. nhistmx ) THEN
nskipt = ncomb*nskipt
mhist = nhist-1
CALL packarr(mhist, thist, ncomb, nhist)
CALL packarr(mhist, tmass, ncomb, nhist)
CALL packarr(mhist, tfmin, ncomb, nhist)
CALL packarr(mhist, tfmax, ncomb, nhist)
CALL packarr(mhist, ermass, ncomb, nhist)
END IF
thist(nhist) = time
tmass(nhist) = SUM(fgrid(1:nx))
!!$ tmass(nhist) = SUM(coefs(1:nx))
tfmin(nhist) = MINVAL(fgrid)
tfmax(nhist) = MAXVAL(fgrid)
ermass(nhist) = (tmass(nhist)-tmass(0))/tmass(0)
END IF
!
CALL xgupdate
END DO
!===========================================================================
! 9.0 Prologue
!
CALL fftw_f77_destroy_plan(forw)
DEALLOCATE(xgrid, fgrid, xshft, coefs)
CONTAINS
DOUBLE PRECISION FUNCTION finit(xx)
!
! A "box" function
!
DOUBLE PRECISION, INTENT(in) :: xx
DOUBLE PRECISION :: xl, xr, xl0, xr0, h, x, xlen
INTEGER :: kl, kr, klflag, krflag
!
xlen = b-a
x = a + MODULO(xx-a+xlen, xlen)
!
xl = 0.375*(b-a)
xr = 0.624*(b-a)
CALL interv(xgrid, nx+1, xl, kl, klflag)
CALL interv(xgrid, nx+1, xr, kr, krflag)
xl0 = xl + w*(xgrid(kl)-xgrid(kl-1))
xr0 = xr - w*(xgrid(kr)-xgrid(kr-1))
CALL interv(xgrid, nx+1, xl0, kl, klflag)
CALL interv(xgrid, nx+1, xr0, kr, krflag)
IF( x .LT. xl0 ) THEN
h = xgrid(kl)-xgrid(kl-1)
finit = EXP(-((x-xl0)/(w*h))**2)
ELSE IF( x .GT. xr0) THEN
h = xgrid(kr)-xgrid(kr-1)
finit = EXP(-((x-xr0)/(w*h))**2)
ELSE
finit = 1.0d0
END IF
END FUNCTION finit
END PROGRAM main
!+++
SUBROUTINE meshdist(c, x, nx)
!
! Construct an 1d non-equidistant mesh given a
! mesh distribution function.
!
IMPLICIT NONE
DOUBLE PRECISION, INTENT(in) :: c(5)
INTEGER, INTENT(iN) :: nx
DOUBLE PRECISION, INTENT(inout) :: x(0:nx)
INTEGER :: nintg
DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: xint, fint
DOUBLE PRECISION :: a, b, dx, f0, f1, scal
INTEGER :: i, k
!
a=x(0)
b=x(nx)
nintg = 10*nx
ALLOCATE(xint(0:nintg), fint(0:nintg))
!
! Mesh distribution
!
dx = (b-a)/REAL(nintg)
xint(0) = a
fint(0) = 0.0d0
f1 = fdist(xint(0))
DO i=1,nintg
f0 = f1
xint(i) = xint(i-1) + dx
f1 = fdist(xint(i))
fint(i) = fint(i-1) + 0.5*(f0+f1)
END DO
!
! Normalization
!
scal = REAL(nx) / fint(nintg)
fint(0:nintg) = fint(0:nintg) * scal
!!$ WRITE(*,'(a/(10f8.3))') 'FINT', fint
!
! Obtain mesh point by (inverse) interpolation
!
k = 1
DO i=1,nintg-1
IF( fint(i) .GE. REAL(k) ) THEN
x(k) = xint(i) + (xint(i+1)-xint(i))/(fint(i+1)-fint(i)) * &
& (k-fint(i))
k = k+1
END IF
END DO
!
DEALLOCATE(xint, fint)
CONTAINS
DOUBLE PRECISION FUNCTION fdist(x)
DOUBLE PRECISION, INTENT(in) :: x
fdist = c(1) + c(2)*x + c(3)*EXP(-((x-c(4))/c(5))**2)
END FUNCTION fdist
END SUBROUTINE meshdist
!+++
SUBROUTINE packarr(n, arr, skip, nhist)
IMPLICIT NONE
INTEGER :: n, skip, i, ii, nhist
DOUBLE PRECISION :: arr(0:n)
ii = 0
DO i=0,n,skip
arr(ii) = arr(i)
ii=ii+1
END DO
nhist = ii
END SUBROUTINE packarr
!+++
SUBROUTINE dump(filename, l)
!
! Is invoked when button "Dump" is pressed.
!
IMPLICIT NONE
CHARACTER(len=*) :: filename
INTEGER :: l
WRITE(*,'(a,a,a1)') 'Dumpfile = "', filename(1:l),'"'
END SUBROUTINE dump
SUBROUTINE quit()
!
! Is invoked when button "Quit" is pressed
!
IMPLICIT NONE
PRINT*, 'Program terminated ...'
END SUBROUTINE quit
c4science · Help