* A. Salih, Dept. of Mechanical Engg., NIT - Trichy, India. *
*************************************************************************
* THIS PROGRAM SOLVES TWO-DIMENSIONAL LINEAR WAVE EQUATION IN A *
* DOUBLY PERIODIC DOMAIN USING FINITE DIFFERENCE METHOD *
* -Schemes: 1. Upwind Scheme *
* 2. MacCormack Method *
* 3. ENO-2 Scheme *
*************************************************************************
*
PROGRAM LINEAR_WAVE_2D
implicit doubleprecision (a-h,o-z)
include 'scalars/integers.inc'
include 'scalars/parameters.inc'
open(unit=11,file='input.dat',status='unknown')
open(unit=21,file='output.dat',status='unknown')
open(unit=31,file='outputs/x.dat',status='unknown')
open(unit=32,file='outputs/y.dat',status='unknown')
open(unit=33,file='outputs/phi.dat',status='unknown')
open(unit=34,file='outputs/time.dat',status='unknown')
open(unit=35,file='outputs/time2.dat',status='unknown')
open(unit=36,file='outputs/phi_trans.dat',status='unknown')
open(unit=37,file='outputs/phi_cent.dat',status='unknown')
c
pi = 4*datan(1d0)
c
CALL READ_IN ! read-in data
close (unit=11)
CALL GRID ! setting up the grid points in the domain
CALL INITIAL_COND ! setting up initial condition
CALL SOLVE ! setting up system of equations
CALL PRINTOUT ! printing out the computational results
CALL WRITE_OUT ! write-out the basic data used for computation
c
stop
end
*
***SUBROUTINE: READING IN THE BASIC DATA *******************************
*
subroutine READ_IN
implicit doubleprecision (a-h,o-z)
include 'array_dimension.inc'
include 'scalars/integers.inc'
include 'scalars/reals.inc'
c
read(11,*) alx ! actual length of domain in x-direction
read(11,*) aly ! actual length of domain in y-direction
read(11,*) ic_type ! initial U profile (1-U_INITIAL1)
read(11,*) u ! value of advection velocity in x-direction
read(11,*) v ! value of advection velocity in y-direction
read(11,*) Co ! Courant number for calculating time step
read(11,*) tfinal ! time at which the solution is desired
read(11,*) scheme ! 1-Upwind, 2-MacCormack, 3-ENO-2
read(11,*) print_freq !
c
c calculation of time-step dt based on the stability condition
dx = alx /(m-1)
dy = aly /(n-1)
dt = Co /(dabs(u)/dx + dabs(v)/dy)
maxntimestp = tfinal /dt
dt = tfinal /maxntimestp ! revised time-step
Co = dt*(dabs(u)/dx + dabs(v)/dy) ! revised Courant number
if (tfinal < dt) then
print*, 'Warning: final time is less than time-step!'
stop
endif
c
return
end
*
***SUBROUTINE: WRITING OUT THE BASIC DATA USED FOR COMPUTATION**********
*
subroutine WRITE_OUT
implicit doubleprecision (a-h,o-z)
include 'array_dimension.inc'
include 'scalars/integers.inc'
include 'scalars/reals.inc'
c
write(21,*) 'Linear, 1-dimensional wave equation with periodic bou
$ndary conditions:'
c
if (scheme == 1) then
write(21,*) 'Upwind Scheme'
elseif (scheme == 2) then
write(21,*) 'MacCormack Method'
elseif (scheme == 3) then
write(21,*) 'Essentially NonOscillatory Scheme'
endif
c
write(21,*) 'Advection velocity in x-direction, u =', u
write(21,*) 'Advection velocity in y-direction, v =', v
write(21,*) 'Courant number, Co =', Co
write(21,*) 'Length of domain in x-direction', alx
write(21,*) 'Length of domain in x-direction', aly
write(21,*) 'Number of grids (including boundaries):', mx, '
$X', ny
write(21,*) 'dx =', dx
write(21,*) 'dy =', dy
write(21,*) 'dt =', dt
write(21,*) 'No. of time step = ', ntimestp
write(21,*) 'Final time = ', time
c
return
end
*
****SUBROUTINE: GRID****************************************************
*
SUBROUTINE GRID
implicit doubleprecision (a-h,o-z)
include 'array_dimension.inc'
include 'scalars/integers.inc'
include 'scalars/reals.inc'
include 'arrays/grid.inc'
c
c m = no. of grid points (including boundary nodes) in x-direction
c n = no. of grid points (including boundary nodes) in y-direction
dx = alx /(m-1)
x(1) = 0
do i = 2, m
x(i) = x(i-1) + dx
enddo
c
dy = aly /(n-1)
y(1) = 0
do j = 2, n
y(j) = y(j-1) + dy
enddo
c
return
end
*
***SUBROUTINE: INITIAL CONDITIONS FOR U ********************************
*
SUBROUTINE INITIAL_COND
implicit doubleprecision (a-h,o-z)
include 'scalars/integers.inc'
c
c defining the initial conditions
if (ic_type == 1) then
CALL PHI_INITIAL1
endif
c
return
end
*
****SUBROUTINE: INITIAL CONDITION
*
SUBROUTINE PHI_INITIAL1
implicit doubleprecision (a-h,o-z)
include 'array_dimension.inc'
include 'scalars/integers.inc'
include 'scalars/reals.inc'
include 'scalars/parameters.inc'
include 'arrays/grid.inc'
include 'arrays/phi.inc'
c
c parameters
half = 1d0/2d0
onet = 1d0/3d0
c
do j = 1, n
do i = 1, m
term = dsqrt( (x(i)-half)**2 + (y(j)-half)**2 )
if (term < onet) then
PHI(i,j) = dcos(2*pi*x(i)) + dcos(2*pi*y(j)) + 2d0
else
PHI(i,j) = dcos(2*pi*x(i)) + dcos(2*pi*y(j)) - 2d0
endif
enddo
enddo
c
return
end
*
***SUBROUTINE: NUMERICAL SOLUTION USING DIFFERENT SCHEMES***************
*
subroutine SOLVE
implicit doubleprecision (a-h,o-z)
include 'array_dimension.inc'
include 'scalars/integers.inc'
include 'scalars/reals.inc'
include 'scalars/parameters.inc'
include 'arrays/grid.inc'
include 'arrays/phi.inc'
doubleprecision PHIN(-1:m1+1,-1:n1+1), PHIPRED(-1:m1+1,-1:n1+1)
doubleprecision minmod_xm,minmod_xp,minmod_ym,minmod_yp
c
c storing the values of u at nth time level
PHIN = PHI
c
ntimestp = 0
maxdif = 1d9
time = 0d0
write(34,101) time ! continuous time
write(35,101) time ! intermittent time
monitorx = m/4
monitory = n/4
write(36,*) PHI(monitorx,monitory)
c
10 if (ntimestp < maxntimestp) then
c periodic boundary conditions
c west
do j = 1,n
do i = 0,-1,-1
PHIN(i,j) = PHIN(m+i-1,j)
enddo
enddo
c east
do j = 1,n
do i = 1,2
PHIN(m+i,j) = PHIN(i+1,j)
enddo
enddo
c south
do i = 1,m
do j = 0,-1,-1
PHIN(i,j) = PHIN(i,n+j-1)
enddo
enddo
c north
do i = 1,m
do j = 1,2
PHIN(i,n+j) = PHIN(i,j+1)
enddo
enddo
c
if (scheme == 1) then ! Upwind
do j = 1,n
do i = 1,m
dPhidx_m = (PHIN(i,j) - PHIN(i-1,j)) /dx
dPhidx_p = (PHIN(i+1,j) - PHIN(i,j)) /dx
dPhidy_m = (PHIN(i,j) - PHIN(i,j-1)) /dy
dPhidy_p = (PHIN(i,j+1) - PHIN(i,j)) /dy
VdPhids = max(u,0d0)*dPhidx_m + min(u,0d0)*dPhidx_p +
$ max(v,0d0)*dPhidy_m + min(v,0d0)*dPhidy_p
PHI(i,j) = PHIN(i,j) - VdPhids*dt
enddo
enddo
c
elseif (scheme == 2) then ! MacCormack
c predictor step
do j = 1,n
do i = 1,m
udPhi = u*(PHIN(i+1,j) - PHIN(i,j))
vdPhi = v*(PHIN(i,j+1) - PHIN(i,j))
PHIPRED(i,j) = PHIN(i,j) - (udPhi/dx + vdPhi/dy)*dt
enddo
enddo
c
do j = 1,n
do i = 0,-1,-1
PHIPRED(i,j) = PHIPRED(m+i-1,j)
enddo
enddo
do j = 1,n
do i = 1,2
PHIPRED(m+i,j) = PHIPRED(i+1,j)
enddo
enddo
do i = 1,m
do j = 0,-1,-1
PHIPRED(i,j) = PHIPRED(i,n+j-1)
enddo
enddo
do i = 1,m
do j = 1,2
PHIPRED(i,n+j) = PHIPRED(i,j+1)
enddo
enddo
c
c corrector step
do j = 1,n
do i = 1,m
udPhi = u*(PHIPRED(i,j) - PHIPRED(i-1,j))
vdPhi = v*(PHIPRED(i,j) - PHIPRED(i,j-1))
PHI(i,j) = (PHIN(i,j) + PHIPRED(i,j))/2d0 -
& (udPhi/dx + vdPhi/dy)*dt/2d0
enddo
enddo
c
elseif (scheme == 3) then ! ENO-2
c second-order Essentially NonOscillatory scheme
do j = 1,n
do i = 1,m
dPhidxxm =(PHIN(i,j) - 2*PHIN(i-1,j) + PHIN(i-2,j))/(dx**2)
dPhidxx =(PHIN(i+1,j) - 2*PHIN(i,j) + PHIN(i-1,j))/(dx**2)
dPhidxxp =(PHIN(i+2,j) - 2*PHIN(i+1,j) + PHIN(i,j))/(dx**2)
dPhidyym =(PHIN(i,j) - 2*PHIN(i,j-1) + PHIN(i,j-2))/(dy**2)
dPhidyy =(PHIN(i,j+1) - 2*PHIN(i,j) + PHIN(i,j-1))/(dy**2)
dPhidyyp =(PHIN(i,j+2) - 2*PHIN(i,j+1) + PHIN(i,j))/(dy**2)
c
if ( dPhidxx * dPhidxxm > 0) then
minmod_xm = dPhidxx/(dabs(dPhidxx)+1d-30) *
$ min( dabs(dPhidxx), dabs(dPhidxxm) )
else
minmod_xm = 0
endif
if ( dPhidxx * dPhidxxp > 0) then
minmod_xp = dPhidxx/(dabs(dPhidxx)+1d-30) *
$ min( dabs(dPhidxx), dabs(dPhidxxp) )
else
minmod_xp = 0
endif
c
if ( dPhidyy * dPhidyym > 0) then
minmod_ym = dPhidyy/(dabs(dPhidyy)+1d-30) *
$ min( dabs(dPhidyy), dabs(dPhidyym) )
else
minmod_ym = 0
endif
if ( dPhidyy * dPhidyyp > 0) then
minmod_yp = dPhidyy/(dabs(dPhidyy)+1d-30) *
$ min( dabs(dPhidyy), dabs(dPhidyyp) )
else
minmod_yp = 0
endif
c
dPhidx_m = (PHIN(i,j) - PHIN(i-1,j))/dx + (dx/2)*minmod_xm
dPhidx_p = (PHIN(i+1,j) - PHIN(i,j))/dx - (dx/2)*minmod_xp
dPhidy_m = (PHIN(i,j) - PHIN(i,j-1))/dy + (dy/2)*minmod_ym
dPhidy_p = (PHIN(i,j+1) - PHIN(i,j))/dy - (dy/2)*minmod_yp
c
VdPhids = max(u,0d0)*dPhidx_m + min(u,0d0)*dPhidx_p +
$ max(v,0d0)*dPhidy_m + min(v,0d0)*dPhidy_p
PHI(i,j) = PHIN(i,j) - VdPhids*dt
enddo
enddo
endif
c
time = time + dt
ntimestp = ntimestp + 1
c
c updating phi values at fictitious nodes
do j = 1,n
do i = 0,-1,-1
PHIN(i,j) = PHIN(m+i-1,j)
enddo
enddo
do j = 1,n
do i = 1,2
PHIN(m+i,j) = PHIN(i+1,j)
enddo
enddo
do i = 1,m
do j = 0,-1,-1
PHIN(i,j) = PHIN(i,n+j-1)
enddo
enddo
do i = 1,m
do j = 1,2
PHIN(i,n+j) = PHIN(i,j+1)
enddo
enddo
c
PHIN = PHI
c
write(34,101) time ! continuous time
write(6,102) ntimestp
c
c printing out transient phi at selected point of domain
j = jcounter + 1
jcounter = mod(j,print_freq)
if (jcounter == 0) then
write(35,101) (time + dt) ! intermittent time
write(36,*) PHI(monitorx,monitory)
endif
go to 10
endif
c
101 format(e14.7)
102 format(i7,5x,e14.7)
*
return
end
*
***SUBROUTINE: PRINTING-OUT THE OUTPUT DATA*****************************
*
SUBROUTINE PRINTOUT
implicit doubleprecision (a-h,o-z)
include 'array_dimension.inc'
include 'scalars/integers.inc'
include 'scalars/reals.inc'
include 'arrays/phi.inc'
include 'arrays/grid.inc'
c
ncent = n/2
write(31,101) (x(i), i=1,m)
write(32,101) (y(i), i=1,n)
write(33,102) ((PHI(i,j), i=1,m), j=1,n)
write(37,101) ((PHI(i,ncent)+ PHI(i,ncent+1)/2), i=1,m) !PHI center
101 format(e13.6)
102 format(128(e13.6,1x))
c
return
end
|
