Sample Computer Programs
Three sample programs are provided in this appendix. The first set of programs is for solving the onedimensional convective wave equation. There are two programs. The first program uses the DRP scheme. The second program uses the standard 2nd, 4th and 6th order central difference scheme for spatial discretization and the RungeKutta scheme for time marching. The second set of programs is for solving the linearized twodimensional Euler equations. The DRP scheme is used. The third set of programs is for solving the nonlinear threedimensional Euler equations. The DRP scheme is used.
G.1 Sample Program for Solving the OneDimensional Convective Wave Equation
The numerical solution to the following problem in an infinite domain is considered. The convective equation in dimensionless form is
d u d u
+ = 0 91 dx
Two types of initial conditions are used in this example
1. A Gaussian function.
2. A box car function i. e.,
u(x, 0) = H(x + L) – H(x – L) where H (x) is the unit step function.
The Fortran 90 program below computes the solution using the DRP scheme.
program DRP_1D
This code solves the 1D convective wave equation on the real line using the DRP scheme
equation is nondimensionalized by mesh spacing and wave speed
equation is discretized with 7point spatial stencils and 4 level multistep time advancement
the solution is advanced from the initial conditions to the desired time and then output to file
Variables:
id – number of mesh points to the left and right of origin nmax – maximum number of time steps t – time
delt – time step u – solution values
k – time derivatives, dudt, at different time levels
kl – converts time level to index in k array a – spatial stencil coefficients
b – multistep advancement coefficients
implicit none
integer, parameter::rkind=8,id=800 integer::n, l,nmax, bh integer, dimension(3:0)::kl real(kind=rkind)::delt, h,t real(kind=rkind), dimension(id3:id+3)::u real(kind=rkind), dimension(3:0,id:id)::k real(kind=rkind), dimension(3:3)::a real(kind=rkind), dimension(3:0)::b integer::outdata=20 !
! Open output data file ‘out. dat’
!
open(unit=outdata, file="outDRP. dat", status="replace") !
! Optimized 7point spatial stencil
!
a(3)=0.02003774846106832_rkind
a(2)=0.163484327177606692_rkind
a(1)=0.766855408972008323_rkind
a(0)=0.0_rkind
a(1)=0.766855408972008323_rkind
a(2)=0.163484327177606692_rkind
a(3)=0.02003774846106832_rkind
!
! Optimized multistep
!
b(0)=2.3025580888_rkind
b(1)=2.4910075998_rkind
b(2)=1.5743409332_rkind
b(3)=0.3858914222_rkind
I
! Set simulation parameters I
delt=0.1_rkind! time step size (stability/accuracy limit is 0.2111) nmax=3000 ! maximum number of time steps I
! Set ‘boundary’ conditions (assume u=0 far from origin)
I
u(id3)=0.0_rkind
u(id2)=0.0_rkind
u(id1)=0.0_rkind
u(id+1)=0.0_rkind
u(id+1)=0.0_rkind
u(id+3)=0.0_rkind
!
! Set initial conditions!
t=0.0_rkind
!
! Gaussian function with halfwidth ‘bh’ and height ‘h’ !bh=3
!h=0.5_rkind! do l=id, id
! u(l)=h*exp(log(2.0_rkind)*(real(l, rkind)/real(bh, rkind))**2) !end do!
! Boxcar
!
u=0.0_rkind
u(50:50)=1.0_rkind
!
! Set solution to zero for t<0
!
k=0.0_rkind
kl(3)=3
kl(2)=2
kl(1)=1
kl(0)=0
!
!——————————————————————
!
! Main loop – march nmax steps with multilevel DRP
do n=1,nmax!
t=t+delt
write (*,*)"n: ",n," t: ",real(t,4)
!
! Update k
!
do l=id, id
k(kl(0),l)=a(3)*u(l3)+a(2)*u(l2)+a(1)*u(l1)+a(1)*u(l+1)+a(2)*u(l+2)+a(3)*u(l+3) end do!
! Update u to time level n
!
do l=id, id
u(l)=u(l)delt*(b(0)*k(kl(0),l)+b(1)*k(kl(1),l)+b(2)*k(kl(2),l)+b(3)*k(kl(3),l)) end do!
! Shift k array back kl =cshift(kl,1)
end do!
! Output u results !
50 format(i8,e16.7e3) do l=id, id
write (outdata,50)l, u(l) end do!
stop
end program DRP_1D
The Fortran 90 program below computes the solution to the convective wave equation by means of the standard 2nd, 4th and 6th order central difference scheme for spatial discretization and the RungeKutta scheme for time marching.
program RK4_1D
This code solves the 1D convective wave equation on the real line using the standard fourth order RungeKutta scheme – equation is nondimensionalized by mesh spacing and wave speed
equation can be discretized with either 6th,4th,2nd order spatial stencils
the solution is advanced from the initial conditions to the desired time and then output to file
Variables:
id – number of mesh points to the left and right of origin nmax – maximum number of time steps t time delt – time step
и, ut – solution values
к, k1,k2 – time derivatives, dudt
a2 
2nd 
order 
spatial 
stencil 
coefficients 
a4 
4th 
order 
spatial 
stencil 
coefficients 
a6 
6th 
order 
spatial 
stencil 
coefficients 
implicit none
integer, parameter::rkind=8,id=800 integer::n, l,nmax, bh, iorder real(kind=rkind)::delt, dth, dt6,h, t real(kind=rkind), dimension(id3:id+3)::u, ut real(kind=rkind), dimension(id:id)::k, k1,k2 real(kind=rkind), dimension(1:1)::a2 real(kind=rkind), dimension(2:2)::a4 real(kind=rkind), dimension(3:3)::a6 integer::outdata=20 !
! Open output data file ‘out. dat’
!
open(unit=outdata, file="outRK4.dat", status="replace") !
! 2nd order spatial stencil
!
a2(1)=0.5_rkind
a2(0)=0.0_rkind
a2(1)=0.5_rkind
!
! 4th order spatial stencil
!
a4(2)=1.0_rkind/12.0_rkind
a4(1)=2.0_rkind/3.0_rkind
a4(0)=0.0_rkind
a4(1)=2.0_rkind/3.0_rkind
a4(2)=1.0_rkind/12.0_rkind
6th order stencil
a6(3)=1.0_rkind/60.0_rkind
a6(2)=3.0_rkind/20.0_rkind
a6(1)=3.0_rkind/4.0_rkind
a6(0)=0.0_rkind
a6(1)=3.0_rkind/4.0_rkind
a6(2)=3.0_rkind/20.0_rkind
a6(3)=1.0_rkind/60.0_rkind
!
! Set simulation parameters!
delt=0.4_rkind! time step size
dth=delt*0.5_rkind
dt6=delt/6.0_rkind
nmax=750 ! maximum number of time steps iorder=6 ! order of the spatial stencil (2,4,6)
!
! Set ‘boundary’ conditions (assume u=0 far from origin) !
u(id3)=0.0_rkind
u(id2)=0.0_rkind
u(id1)=0.0_rkind
u(id+1)=0.0_rkind
u(id+1)=0.0_rkind
u(id+3)=0.0_rkind
!
ut=0.0_rkind
!
! Set initial conditions!
t=0.0_rkind
!
! Gaussian function with halfwidth ‘bh’ and height ‘h’ !bh=3
!h=0.5_rkind! do l=id, id
! u(l)=h*exp(log(2.0_rkind)*(real(l, rkind)/real(bh, rkind))**2) !end do!
! Boxcar!
u=0.0_rkind
u(50:50)=1.0_rkind
Initial du/dt
call dudt(u, k,iorder)
!
!………………………………………………………….
!
! Main loop – march nmax steps in time using RK4 integration!
do n=1,nmax!
t=t+delt
write (*,*)”n: ",n," t: M, real(t,4)
!
! Second step!
do l=id, id ut(l)=u(l)+dth*k(l) end do
call dudt(ut, k1,iorder)
!
! Third step
!
do l=id, id ut(l)=u(l)+dth*k1(l) end do
call dudt(ut, k2,iorder)
!
! Fourth step
!
do l=id, id ut(l)=u(l)+delt*k2(l) k2(l)=k2(l)+k1(l) end do
call dudt(ut, k1,iorder)
!
! Update u
!
do l=id, id
u(l)=u(l)+dt6*(k(l)+k1(l)+2.0_rkind*k2(l)) end do!
! Update dudt
!
call dudt(u, k,iorder) end do
! Output u results I
50 format(i8,e16.7e3) do l=id, id
write (outdata,50)l, u(l) end do!
stop
!
contains
!
subroutine dudt(f, kf, iflag) !
!
! Calculate the derivative du/dt=du/dx of the convectice wave eqn! – iflag is a flag for the spatial stencil (2nd,4th,6th order)
!
!
implicit none
integer, intent(in)::iflag
real(kind=rkind), intent(in), dimension(id3:id+3)::f real(kind=rkind), intent(out), dimension(id:id)::kf integer::l!
select case (iflag) case (2) do l =id, id
kf(l)=(a2(1)*f(l1)+a2(1)*f(l+1)) end do case (4) do l =id, id
kf(l)=(a4(2)*f(l2)+a4(1)*f(l1)+a4(1)*f(l+1)+a4(2)*f(l+2)) end do case (6) do l =id, id
kf(l)=(a6(3)*f(l3)+a6(2)*f(l2)+a6(1)*f(l1)+a6(1)*f(l+1)+a6(2)*f(l+2)+a6(3)*f(l+3)) end do end select!
return
end subroutine dudt !
end program RK4_1D





Boxcar, t=300, DRP Scheme x 
Boxcar, t=300, RK42nd 

1.1. The governing equation and boundary conditions for a rectangular vibrating membrane of width W and breadth B are
92U 2 / d2Ы d2U
dt2 dx2^~ dy2
At x = 0, x = W, y = 0, y = B, u = 0,
where U is the displacement of the membrane. It is easy to show that the normal mode frequencies of the vibrating membrane are
[2]/2
, where j and к are integers.
Now divide the width into N subdivisions and the breadth into M subdivisions so that Ax = W/N and Ay = B/M. Suppose a secondorder central difference is used to approximate the x and y derivatives. This converts the partial differential equation into a partial difference equation, as follows:
[3]
The conditions for E to be a minimum are
[4] For convenience, D(a Ax) may be normalized so that
D (n) = 1. (7.18)
The right side of Eq. (7.14) is a truncated Fourier cosine series. To ensure that D (a Ax) is small when a Ax is small but becomes large when a Ax ^ n, one may use a Gaussian function centered at a Ax = n with halfwidth a as a template for choosing the Fourier coefficients d. This is done by determining dj such that the integral
[5] There should be no damping for long waves. That is, D (a Ax) ^ 0 as aAx^0. This requires
3
d0 + 2 dj = 0. (7.17)
j=1
7.1. Let us reconsider the initial value problem
9 u d u
+ = 0 9t dx
t = 0, u = H (x + 50) – H (x – 50), where H (x) is the unit step function.
Compute the solution using the 7point DRP scheme with artificial selective damping. Use a At that satisfies the numerical stability requirement.
[7] Use RA = 0.3, where Ra is the artificial mesh Reynolds number.
[8] Use RA = 0.3/At.
Plot u(x, t) at t = 300 for each case.
7.2. Real fluids have molecular viscosity that tends to damp out long and short wave disturbances. This problem focuses on how best to compute viscous terms to simulate real viscous effects. At the same time, it provides an example illustrating why
[9] Use RA = 3.
[10] 3
— aj ФJ.
Ді j j
j=3
[11]+dMA)U+A і+(1 – m2 ) 2 ВI=<938)
where u(x, y, t) = U(x, y)e1 ®f.
The second step is to introduce damping into the equation through a socalled complex change of variables. Let о x(x) > 0, о (y) > 0 be the absorption coefficients. The complex change of variables is to let
9.1. Many aeroacoustic problems involve the scattering of vorticity or acoustic waves by an object. To simulate this class of problems in the time domain, the incoming waves are generated by the boundary conditions imposed on the artificial external boundaries of the computation domain.
Dimensionless variables with respect to length scale L (L will be chosen later), velocity scale aTO (ambient sound speed), time scale L/aTO, density scale pTO (ambient density), and pressure scale p^a^ are used. Suppose there is a uniform mean flow
11.1. The optimized extrapolation method can be applied even when the stencil points are not spaced uniformly apart. Let the coordinates of the stencil points be x0, x_1, x_2,…, x__7 (a 7point stencil). Suppose the extrapolation point is at x = x0 + n. The extrapolation formula is
[14]
f (x0 + n) = Y^ Sjf(x j )•
j=0
[15]
[16] (R, в, ф, ю)
[17] 0 < r < 0.8