Two Dimensional Wave Equation with inflow in x and periodic in y.
10 time frames of a circle moving up from the lower left corner and exit
to the right boundary.
Two Dimensional Wave Equation with inflow in x and periodic in y.
10 time frames of a circle moving up from the lower left corner and exit
to the right boundary.
Note : Portions of Code that colored
Red are essential for used with PseudoPack.
List of Routines
! | Back to the Top |
!------
! Demonstration Program for the use of PseudoPack.
!
! 2-D Scalar Wave Equation
!
! Chebyshev/Fourier Implementation
!
! 1: U_t+a*U_x+b*U_y=0
! 2: U_t+y*U_x-x*U_y=0
!
!
! Initial Condition : U(x,y,0) = Smoothed Gaussian
! Boundary Condition : Periodic in y direction
! Inflow/Outflow in x direction
!
!------
c------
c Must include this include statement for used with Macros in pseudopack.h
c------
#include "pseudopack.h"
!
!============================================================
!
MODULE Initial_Parameter
Integer :: Flux_Choice
REALTYPE :: Speed_x, Speed_y
REALTYPE :: C_x, C_y, Radius_Ratio, Gaussian_Order
REALTYPE :: Left_x, Right_x, Bottom_y, Top_y
END MODULE Initial_Parameter
!
!============================================================
!
Program Wave_2D
USE PseudoPack
USE Initial_Parameter
implicit NONE
integer, parameter :: N_Frame = 10
Integer :: Index_x, Method_x, Point_Type_x, Max_Diff_x, &
Algorithm_x, Symmetry_x, &
N_x, BlockSize_x, Map_x, Manual_x, Map_F_x, &
Filter_Choice_D_x, Mode_CutOff_D_x, &
Smooth_1_D_x, Smooth_2_D_x, &
Filter_Choice_S_x, Mode_CutOff_S_x
REALTYPE :: alpha_x, Angle_x, &
Omega_D_x, Order_D_x, Omega_S_x, Order_S_x
Integer :: Index_y, Method_y, Point_Type_y, Max_Diff_y, &
Algorithm_y, Symmetry_y, &
N_y, BlockSize_y, Map_y, Manual_y, Map_F_y, &
Filter_Choice_D_y, Mode_CutOff_D_y, &
Smooth_1_D_y, Smooth_2_D_y, &
Filter_Choice_S_y, Mode_CutOff_S_y
REALTYPE :: alpha_y, Angle_y, &
Omega_D_y, Order_D_y, Omega_S_y, Order_S_y
Integer :: N, M, LDY, M_D, M_S
Integer :: Step, N_Time_Step, N_Inc, Status
logical :: BlowUp
REALTYPE :: Time, Final_Time, dt, CFL
REALTYPE, dimension(:,:), ALLOCATABLE :: u
! | Back to the Top |
REALTYPE, dimension(:) , ALLOCATABLE :: x, y
REALTYPE, dimension(:) , ALLOCATABLE :: D_x, D_y, S_x, S_y
TYPE (PS_Property ) :: Property
TYPE (PS_Grid_Index) :: Grid_Index
TYPE (PS_Domain ) :: Domain
TYPE (PS_Mapping ) :: Mapping
TYPE (PS_Filtering ) :: Filtering_D, Filtering_S
!--------------------------------------------------------------
! First Executable Statement start here.....
!--------------------------------------------------------------
!------
! Input Data
!------
call Input (7)
N = N_x-1
M = N_y-1
Allocate (x(0:N), y(0:M), u(0:N,0:M))
! | Back to the Top |
!------
! Specify the leading dimension of the data array
!------
LDY = SIZE(u, DIM=1)
!------
! Setup in x direction
!------
call PS_Setup ( Property=Property , Index=Index_x, Method=Method_x, &
Algorithm=Algorithm_x, &
Grid_Index=Grid_Index , N=N_x, M=N_y, LDY=LDY, &
Mapping=Mapping , Map=Map_x)
call PS_Setup_Filtering (Filtering_D , Filter_Choice_D_x)
call PS_Setup_Filtering (Filtering_S , Filter_Choice_S_x)
call PS_Get_Operator_Size ('D', N_x, M_D, Property)
Allocate (D_x(M_D))
call PS_Setup_Operator (N_x, D_x, x, Property, Grid_Index, &
Mapping=Mapping, Filtering=Filtering_D)
call PS_Get_Operator_Size ('S', N_x, M_S, Filtering_S, Property)
Allocate (S_x(M_S))
call PS_Setup_Operator (N_x, S_x, Property, Grid_Index, Filtering_S)
!------
! Setup in y direction
!------
call PS_Setup ( Property=Property , Index=Index_y, Method=Method_y, &
Algorithm=Algorithm_y, &
Grid_Index=Grid_Index , N=N_y, M=N_x, LDY=LDY, &
Domain=Domain , x0=Bottom_y, x1=Top_y, &
Map_F=Map_F_y, &
Mapping=Mapping , Map=Map_y)
call PS_Setup_Filtering (Filtering_D , Filter_Choice_D_y)
call PS_Setup_Filtering (Filtering_S , Filter_Choice_S_y)
call PS_Get_Operator_Size ('D', N_y, M_D, Property)
Allocate (D_y(M_D))
call PS_Setup_Operator (N_y, D_y, y, Property, Grid_Index, &
Domain, Filtering=Filtering_D)
call PS_Get_Operator_Size ('S', N_y, M_S, Filtering_S, Property)
Allocate (S_y(M_S))
call PS_Setup_Operator (N_y, S_y, Property, Grid_Index, Filtering_S)
!------
! Finished Setup
!------
call Echo_Code_Information
!------
! Compute Initial Condition
!------
call Initial_Conditions (N, M, x, y, u, LDY)
!------
! Compute Time Step dt
!------
call Time_Step (N, M, x, y, u, LDY, CFL, dt)
!------
! Output Initial Condition
!------
OPEN (UNIT=10, FILE='wave_Movie.dat', STATUS="UNKNOWN") ; REWIND (10)
OPEN (UNIT=11, FILE='wave.dat' , STATUS="UNKNOWN") ; REWIND (11)
OPEN (UNIT=12, FILE='wave_Movie.m' , STATUS="UNKNOWN") ; REWIND (12)
OPEN (UNIT=13, FILE='wave.m' , STATUS="UNKNOWN") ; REWIND (13)
Status = 0
call TecPlot (N, M, Step, Time, x, y, u, LDY, Status, 10)
call Matlab (N, M, Step, Time, x, y, u, LDY, Status, 12)
Status = 1
N_Time_Step = INT(Final_Time/dt)
N_Inc = MAX(N_Time_Step/(N_Frame-2),1)
call Echo_Time_Step_Information
!------
! Start the Runge Kutta Time Stepping
!------
Step = 0
Time = 0
100 continue
Step = Step + 1
Time = Time + dt
call RK_TVD (N, M, Time, dt, x, y, u, LDY, D_x, D_y, S_x, S_y)
if (Step/N_Inc*N_Inc == Step) then
write (6,*) 'The Current Step, Time is ',Step, Time
call TecPlot (N, M, Step, Time, x, y, u, LDY, Status, 10)
call Matlab (N, M, Step, Time, x, y, u, LDY, Status, 12)
endif
BlowUp = MAXVAL(ABS(u)) > TWO ! Do not fool around with this variable
if (BlowUp) then
Final_Time = Time
write (6,*)
write (6,*) '--------------------------------------------------'
write (6,*) 'Stopped! Unstable Numerical Solution ....'
write (6,*) '--------------------------------------------------'
write (6,*)
endif
if (Time < Final_Time) goto 100
!------
! Output Data
!------
call TecPlot (N, M, Step, Time, x, y, u, LDY, Status, 10)
call Matlab (N, M, Step, Time, x, y, u, LDY, Status, 12)
Status = 0
call TecPlot (N, M, Step, Time, x, y, u, LDY, Status, 11)
call Matlab (N, M, Step, Time, x, y, u, LDY, Status, 13)
call Echo_Final_Message
CONTAINS
! | Back to the Top |
Subroutine Input (LID)
integer :: LID, IOS
OPEN (UNIT=LID, FILE='Wave.input', STATUS='OLD', IOSTAT=IOS)
call PS_Input (Index_x, Method_x, Point_Type_x, Max_Diff_x, &
Algorithm_x, Symmetry_x, &
N_x, BlockSize_x, Left_x, Right_x, Map_F_x, &
Map_x, Manual_x, alpha_x, Angle_x, &
Filter_Choice_D_x, Mode_CutOff_D_x, &
Omega_D_x, Order_D_x, &
Smooth_1_D_x, Smooth_2_D_x, &
Filter_Choice_S_x, Mode_CutOff_S_x, &
Omega_S_x, Order_S_x, &
LID)
call PS_Input (Index_y, Method_y, Point_Type_y, Max_Diff_y, &
Algorithm_y, Symmetry_y, &
N_y, BlockSize_y, Bottom_y, Top_y, Map_F_y, &
Map_y, Manual_y, alpha_y, Angle_y, &
Filter_Choice_D_y, Mode_CutOff_D_y, &
Omega_D_y, Order_D_y, &
Smooth_1_D_y, Smooth_2_D_y, &
Filter_Choice_S_y, Mode_CutOff_S_y, &
Omega_S_y, Order_S_y, &
LID)
read (LID,'(//)')
read (LID,*) Final_Time
read (LID,*) CFL
read (LID,'(//)')
read (LID,*) Flux_Choice
read (LID,*)
read (LID,*) Speed_x, Speed_y
read (LID,*) C_x, C_y
read (LID,*) Radius_Ratio
read (LID,*) Gaussian_Order
CLOSE (LID)
if (Flux_Choice == 1) then
Speed_x = ONE ; Speed_y = ONE
endif
END Subroutine Input
! ==================================================================
Subroutine Echo_Code_Information
write (6,*)
write (6,*) '--------------------------------------------------'
write (6,*)
if (OPERATOR_METHOD(D_x) == 0) &
write (6,*) 'Method in x is Fourier'
if (OPERATOR_METHOD(D_x) == 1) &
write (6,*) 'Method in x is Chebyshev'
if (OPERATOR_METHOD(D_x) == 2) &
write (6,*) 'Method in x is Legendre '
if (OPERATOR_POINT_TYPE(D_x) == 1) &
write (6,*) 'Point_type in x is Gauss-Lobatto'
if (OPERATOR_POINT_TYPE(D_x) == 2) &
write (6,*) 'Point_type in x is Gauss_Radau'
if (OPERATOR_POINT_TYPE(D_x) == 3) &
write (6,*) 'Point_type in x is Gauss'
if (OPERATOR_ALGORITHM(D_x) == 0) &
write (6,*) 'Algorithm in x is MXM'
if (OPERATOR_ALGORITHM(D_x) == 1) &
write (6,*) 'Algorithm in x is EOD'
if (OPERATOR_ALGORITHM(D_x) == 2) &
write (6,*) 'Algorithm in x is Transform'
if (OPERATOR_METHOD(D_y) == 0) &
write (6,*) 'Method in y is Fourier'
if (OPERATOR_METHOD(D_y) == 1) &
write (6,*) 'Method in y is Chebyshev'
if (OPERATOR_METHOD(D_y) == 2) &
write (6,*) 'Method in y is Legendre '
if (OPERATOR_POINT_TYPE(D_y) == 1) &
write (6,*) 'Point_type in y is Gauss-Lobatto'
if (OPERATOR_POINT_TYPE(D_y) == 2) &
write (6,*) 'Point_type in y is Gauss_Radau'
if (OPERATOR_POINT_TYPE(D_y) == 3) &
write (6,*) 'Point_type in y is Gauss'
if (OPERATOR_ALGORITHM(D_y) == 0) &
write (6,*) 'Algorithm in y is MXM'
if (OPERATOR_ALGORITHM(D_y) == 1) &
write (6,*) 'Algorithm in y is EOD'
if (OPERATOR_ALGORITHM(D_y) == 2) &
write (6,*) 'Algorithm in y is Transform'
write (6,*)
write (6,*) '--------------------------------------------------'
write (6,*)
END Subroutine Echo_Code_Information
! ==================================================================
Subroutine Echo_Time_Step_Information
write (6,*)
write (6,*) '--------------------------------------------------'
write (6,*)
write (6,*) ' Time Step is ',dt
write (6,*) ' Final Time is ',Final_Time
write (6,*) 'Total Number of Time Step is ',int(Final_Time/dt)
write (6,*)
write (6,*) 'Frame Increment Time Step is ',N_Inc
write (6,*) ' Number of Frame is ',N_Frame
write (6,*)
write (6,*) '--------------------------------------------------'
write (6,*)
END Subroutine Echo_Time_Step_Information
! ==================================================================
Subroutine Echo_Final_Message
write (6,*)
write (6,*) '--------------------------------------------------'
write (6,*) 'Files wave_Movie.dat and wave.dat are output in', &
' TecPlot .dat format.'
write (6,*)
write (6,*) ' wave_Movie.dat --- contains a sequence of solution'
write (6,*) ' wave.dat --- contains Solution at Time ',Time
write (6,*)
write (6,*) 'Files wave_Movie.m and wave.m are output in', &
' Matlab .m format.'
write (6,*)
write (6,*) ' wave_Movie.m --- contains a sequence of solution'
write (6,*) ' wave.m --- contains Solution at Time ',Time
write (6,*)
write (6,*) 'If user does not have TecPlot or Matlab, Please replace'
write (6,*) " the plotting Subroutine TecPlot with user's own"
write (6,*) ' prefered plotting subroutine!'
write (6,*)
write (6,*) 'Thank you for trying this package! -- WSDON'
write (6,*)
END Subroutine Echo_Final_Message
END Program
!
! ==================================================================
! | Back to the Top |
Subroutine Initial_Conditions (N, M, x, y, u, LDY)
USE Initial_Parameter, Order=>Gaussian_Order
implicit NONE
Integer :: N, M, LDY, i, j
REALTYPE :: alpha, Radius, R, L_x, L_y, xx, yy
REALTYPE, dimension(0:N) :: x
REALTYPE, dimension(0:M) :: y
REALTYPE, dimension(0:LDY-1,0:M) :: u
alpha = -LOG(EPSILON(ONE))
L_x = Right_x-Left_x
L_y = Top_y-Bottom_y
Radius = MIN(L_y,L_x)*Radius_Ratio
if (Flux_Choice == 1) then
C_x = Left_x+HALF*L_x
C_y = Bottom_y+HALF*L_y
else
C_x = Left_x+Radius
C_y = Bottom_y+Radius
endif
do j = 0,M
do i = 0,N
xx = x(i)-C_x
yy = y(j)-C_y
R = SQRT(xx**2+yy**2)
if (R <= Radius) then
u(i,j) = EXP(-alpha*(R/Radius)**Order)
else
u(i,j) = ZERO
endif
enddo
enddo
END Subroutine Initial_Conditions
!
! ==================================================================
! | Back to the Top |
Subroutine Time_Step (N, M, x, y, u, LDY, CFL, dt)
USE Initial_Parameter, ONLY: Speed_x, Speed_y
implicit NONE
Integer :: N, M, LDY, i, j
REALTYPE :: CFL, dt, ddt, dx, dy
REALTYPE, dimension(0:N) :: x
REALTYPE, dimension(0:M) :: y
REALTYPE, dimension(0:LDY-1,0:M) :: u
ddt = ZERO
do j = 1,M
do i = 1,N
dx = Speed_x/ABS(x(i)-x(i-1))
dy = Speed_y/ABS(y(j)-y(j-1))*PI
ddt = MAX(dx+dy, ddt)
enddo
enddo
dt = CFL/ddt
END Subroutine Time_Step
!
! ==================================================================
! | Back to the Top |
Subroutine RK_TVD (N, M, Time, dt, x, y, u, LDY, D_x, D_y, S_x, S_y)
USE PseudoPack
implicit NONE
integer :: N, M, LDY
REALTYPE :: Time, dt, Time_n, Time_Now
REALTYPE, dimension(0:N) :: x
REALTYPE, dimension(0:M) :: y
REALTYPE, dimension(0:LDY-1,0:M) :: u, u1, du
REALTYPE, dimension(*) :: D_x, D_y, S_x, S_y
Time_n = Time-dt
Time_Now = Time_n
!
! Stage 1 :
!
call Flux (N, M, x, y, u , du, LDY, D_x, D_y)
u1(0:N,:) = u(0:N,:) - dt*du(0:N,:)
Time_Now = Time_n+dt
call Boundary_Condition (N, M, Time_Now, x, y, u1, LDY)
!
! Stage 2 :
!
call Flux (N, M, x, y, u1, du, LDY, D_x, D_y)
u1(0:N,:) = (THREE*u(0:N,:) + u1(0:N,:) - dt*du(0:N,:))/FOUR
Time_Now = Time_n+dt/2
call Boundary_Condition (N, M, Time_Now, x, y, u1, LDY)
!
! Stage 3 :
!
call Flux (N, M, x, y, u1, du, LDY, D_x, D_y)
u(0:N,:) = (u(0:N,:)+TWO*u1(0:N,:)-TWO*dt*du(0:N,:))/THREE
Time_Now = Time_n+dt
!------
! Smooth the solution u along the First index
!------
call PS_Smooth (LDY, S_x, u, M+1)
!------
! Smooth the solution u along the Second index
!------
call PS_Smooth (LDY, S_y, u, N+1)
!------
! Done Smoothing
!------
call Boundary_Condition (N, M, Time_Now, x, y, u, LDY)
END Subroutine RK_TVD
!
! ==================================================================
! | Back to the Top |
Subroutine Flux (N, M, x, y, u, du, LDY, D_x, D_y)
USE PseudoPack
USE Initial_Parameter, ONLY: Flux_Choice, Speed_x, Speed_y, C_x, C_y
integer :: N, M, LDY, i, j
REALTYPE :: xx, yy
REALTYPE, dimension(0:N) :: x
REALTYPE, dimension(0:M) :: y
REALTYPE, dimension(0:LDY-1,0:M) :: u, du, df
REALTYPE, dimension(*) :: D_x, D_y
!------
! Compute the du/dx along the First index
!------
call PS_Diff (LDY, D_x, u, df, M+1)
!------
! Compute the du/dy along the Second index
!------
call PS_Diff (LDY, D_y, u, du, N+1)
!------
! Compute the Total Flux
!------
SELECT CASE (Flux_Choice)
CASE (0)
du(0:N,:) = Speed_x*df(0:N,:)+Speed_y*du(0:N,:)
CASE (1)
do j = 0,M
do i = 0,N
xx = x(i)-C_x
yy = y(j)-C_y
du(i,j) = yy*df(i,j)-xx*du(i,j)
enddo
enddo
END SELECT
END Subroutine Flux
!
! ==================================================================
! | Back to the Top |
Subroutine Boundary_Condition (N, M, Time, x, y, u, LDY)
implicit NONE
Integer :: N, M, LDY
REALTYPE :: Time
REALTYPE, dimension(0:N) :: x
REALTYPE, dimension(0:M) :: y
REALTYPE, dimension(0:LDY-1,0:M) :: u
u(N,:) = ZERO
END Subroutine Boundary_Condition
!
! ==================================================================
! | Back to the Top |
Subroutine TecPlot (N, M, Step, Time, x, y, u, LDY, Status, lid)
implicit NONE
integer :: N, M, LDY, lid, Step, Status, i, j
REALTYPE :: Time
REALTYPE, dimension(0:N) :: x
REALTYPE, dimension(0:M) :: y
REALTYPE, dimension(0:LDY-1,0:M) :: u
if (Status == 0) then
write (lid,1000) Step, Time
write (lid,1001)
write (lid,1002) N+1, M+1
write (lid,1100) ((x(i) , i=0,N),j=0,M)
write (lid,1100) ((y(j) , i=0,N),j=0,M)
else
write (lid,1003) N+1, M+1
endif
write (lid,1100) ((u(i,j) , i=0,N),j=0,M)
1000 format (1x,'Title="2-D Linear Wave Equation: Step=',i5,' Time=',f10.7,'"')
1001 format (1x,'Variables=x,y, u')
1002 format (1x,'Zone I=',i4,' J=',i4,' F=Block')
1003 format (1x,'Zone I=',i4,' J=',i4,' F=Block, D=(1,2)')
1100 format ( 8(1x,e18.11))
END Subroutine TecPlot
!
! ==================================================================
!
Subroutine Matlab (N, M, Step, Time, x, y, u, LDY, Status, lid)
implicit NONE
integer :: N, M, LDY, lid, Step, Status, i, j
REALTYPE :: Time
REALTYPE, dimension(0:N) :: x
REALTYPE, dimension(0:M) :: y
REALTYPE, dimension(0:LDY-1,0:M) :: u
if (Status == 0) then
write (lid,1000) 'x=transpose(['
write (lid,1001) (x(i) , i=0,N)
write (lid,1000) ']);'
write (lid,1000) 'y=transpose(['
write (lid,1001) (y(j) , j=0,M)
write (lid,1000) ']);'
endif
write (lid,1000) 'z=reshape(['
write (lid,1001) ((u(i,j) , j=0,M),i=0,N)
write (lid,1000) '],size(y,2),size(x,2));'
write (lid,1003)
1000 format (1x,a40)
1001 format (e13.4)
1003 format (1x,'surf(x,y,z); shading interp; pause(2)')
END Subroutine Matlab
! | Back to the Top |