! ! Compile with ! f95 -O4 -r8 -f77 -o burg burg.f90 -L/usr/local/lib -lfftw -lrfftw module constants complex, parameter :: zi = ( 0.0 , 1.0 ) real, parameter :: pi = 3.1415926535897931 end module constants module fft implicit none ! ! For scientific applications that require a fast Fourier transform, ! one wishes to use an optimized set of subroutines to obtain ! good performance. The best portable and well-documented FFT ! library that I am familiar with is called "FFTW", and is the ! one that we will use below. FFTW is freely available at ! http://www.fftw.org, where there is also extensive documentation. ! ! There is one important usage note, which I include for future ! reference: : If you wish to perform one-dimensional ! real-to-complex or complex-to-real transforms, I recommend ! using the rfftwnd set of calls, with the number of dimensions set ! to one. Otherwise, you will likely find that FFTW's ordering of ! array elements is inconvenient and confusing. ! ! The FFTW routines are not written in Fortran. To call them ! from Fortran programs requires a careful look at the documentation, ! or you can follow the examples below. The interface is reasonably ! intuitive. ! integer fftw_forward,fftw_backward parameter (fftw_forward=-1,fftw_backward=1) integer fftw_real_to_complex,fftw_complex_to_real parameter (fftw_real_to_complex=-1,fftw_complex_to_real=1) integer fftw_estimate,fftw_measure parameter (fftw_estimate=0,fftw_measure=1) integer fftw_out_of_place,fftw_in_place parameter (fftw_out_of_place=0) parameter (fftw_in_place=8) end module fft module nonlinearity interface nonlinear module procedure pseudospectral, spectral end interface contains subroutine pseudospectral (a, aprime, nl) use fft use constants complex, dimension(:), intent (in) :: a, aprime complex, dimension(:), intent (out) :: nl complex, dimension(:), allocatable, save :: a_tmp, aprime_tmp, nl_tmp real, dimension(:), allocatable, save :: u, uprime, nl_x ! ! kmax will be the maximum k, while n will be the number of ! grid points used to evaluate the FFT ! integer :: kmax, n, j, mult integer, save :: plan_x2k, plan_k2x logical, save :: first = .true. kmax = size(a) n = 2*(kmax-1) ! we might find we need to change this later ! mult = 1 n = (n*3)/2 if (first) then call rfftwnd_f77_create_plan(plan_k2x,1,n,fftw_backward,fftw_estimate) call rfftwnd_f77_create_plan(plan_x2k,1,n,fftw_forward,fftw_estimate) allocate (u(n), uprime(n), nl_x(n)) allocate (a_tmp((kmax*3)/2), aprime_tmp((kmax*3)/2), nl_tmp((kmax*3)/2)) first = .false. end if a_tmp = 0. aprime_tmp = 0. do j = 1, kmax a_tmp(j) = a(j) aprime_tmp(j) = aprime(j) end do call rfftwnd_f77_one_complex_to_real(plan_k2x,a_tmp,u) call rfftwnd_f77_one_complex_to_real(plan_k2x,aprime_tmp,uprime) nl_x = u*uprime call rfftwnd_f77_one_real_to_complex(plan_x2k,nl_x,nl_tmp) nl = nl_tmp(1:kmax)/n end subroutine pseudospectral subroutine spectral (a, nl) use constants complex, dimension(0:) :: a, nl complex :: aptmp, aqtmp integer :: p, q integer :: kmax kmax = size(a)-1 ! ! Calculate the nonlinear term from a given array of "a" values ! nl = 0. do p = -kmax, kmax do q = -kmax, kmax if (abs(p+q) > kmax) cycle if (p+q < 0) cycle if (p < 0) then aptmp = conjg(a(-p)) else aptmp = a(p) end if if (q < 0) then aqtmp = conjg(a(-q)) else aqtmp = a(q) end if nl(p+q) = nl(p+q) + zi * p * aptmp * aqtmp end do end do end subroutine spectral end module nonlinearity !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!! Main Program starts here !!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! program burg ! ! Reduce typographical errors with implicit none: ! implicit none ! ! For convenience, define i and pi ! complex, parameter :: zi = ( 0.0 , 1.0 ) real, parameter :: pi = 3.1415926535897931 ! ! Defer specification of problem size until runtime ! complex, dimension(:), allocatable :: a, astar, aold, nl real, dimension(:), allocatable :: gam ! ! Declare a few control variables ! real :: dt, nu, gamma, pert, time, fac integer :: nsteps, n, kmax, nout, kinj integer :: k logical :: spectral ! ! Read inputs, initialize arrays ! call init ! ! Main loop ! do n = 1, nsteps call step ! ! Every nout steps, do some output to make time-dependent graphs later. ! if (mod(n,nout) == 0) call output if (mod(n,nout) == 0) call final_output end do ! ! Write final info to a file ! call final_output ! ! This is the end of the main program block. ! contains ! ! Subroutines that are "contained" have access to all ! variables declared in the containing procedure, so ! common blocks are not needed. ! subroutine init ! ! Declare some local variables. ! real :: x, y ! ! Define the input namelist ! namelist /input/ dt, nu, kinj, gamma, kmax, nsteps, pert, nout, & fac, spectral ! ! set default values before reading the input file ! spectral = .true. ! ! Open the input file, read the namelist, and close it back up. ! open (unit=10, file='burg.in', status='old') read (10, input) close (10) ! ! Allocate storage for the calculation. Note that we do not need to keep ! information for k < 0 because the reality condition determines them from ! the relation a(-k) = conjg(a(k)). The total amount of information is ! the same as in the real-space representation, since the main arrays are ! complex. ! allocate (a (0:kmax)) allocate (aold (0:kmax)) allocate (astar(0:kmax)) allocate (nl (0:kmax)) allocate (gam (0:kmax)) ! ! Set allocated variables to zero initially. ! a = 0. aold = 0. astar = 0. nl = 0. gam = 0. ! ! Put in initial condition. ! if (pert > 0.) then do k = 1, kmax call random_number(x) call random_number(y) aold(k) = cmplx(x,y) * pert end do else aold(1) = -pert end if ! ! Set the growth rate and k for the injected energy. ! gam(kinj) = gamma ! ! Set time = 0. ! time = 0. ! ! Open unit 11 for later writes. ! open (unit=11, file='burg.out') end subroutine init subroutine step use nonlinearity, only: nonlinear complex, dimension(0:kmax) :: aoldprime, astarprime ! ! This is a second-order Runge-Kutta scheme in time. ! ! ! Take half-step first with aold on the RHS: ! ! Linear terms: ! do k = 0, kmax astar(k) = aold(k) * (1. - 0.5*nu*k**2*dt + 0.5*gam(k)*dt) end do ! ! Nonlinear term: ! if (spectral) then call nonlinear (aold, nl) else do k=0,kmax aoldprime(k) = zi*k*aold(k) end do call nonlinear (aold, aoldprime, nl) end if do k = 0, kmax astar(k) = astar(k) - 0.5*dt*nl(k)*fac end do ! ! Then take full step from aold, using astar on the RHS: ! do k = 0, kmax a(k) = aold(k) - nu*k**2*dt*astar(k) + gam(k)*dt*astar(k) end do if (spectral) then call nonlinear (astar, nl) else do k=0, kmax astarprime(k) = zi*k*astar(k) end do call nonlinear (astar, astarprime, nl) end if do k = 0, kmax a(k) = a(k) - dt*nl(k)*fac end do ! ! Update the "old" value ! aold = a ! ! Update the time for diagnostics. ! time = time + dt end subroutine step subroutine output ! ! Put quantities that you want time histories of here: ! write(11,*) time, real(a(1)), aimag(a(1)), real(sum(a*conjg(a))), & real(a(1)*conjg(a(1))) end subroutine output subroutine final_output ! ! These are output statements called at the end of the program. ! complex, allocatable, dimension(:) :: u integer :: j logical :: first = .true. if (first) then open (unit=12, file='burg.out2') open (unit=13, file='burg.x') write(13,*) '# x u(x)' first = .false. else write(12, *) write(13,*) end if do k = 0, kmax write(12,*) k, real(a(k)), aimag(a(k)), real(a(k)*conjg(a(k))) end do ! ! Construct real-space representation of u from spectral coefficients: ! allocate (u(-kmax:kmax)) u = 0. do j = -kmax, kmax ! ! These loops could be combined: ! do k = 1, kmax u(j) = u(j) + a(k)*exp(zi*k*pi*j/kmax) end do do k = 1, kmax u(j) = u(j) + conjg(a(k))*exp(-zi*k*pi*j/kmax) end do end do do j = -kmax, kmax write(13,*) pi*j/kmax, real(u(j))/2. end do deallocate (u) end subroutine final_output end program burg