1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
|
SUBROUTINE acoustic3d_f90(N, niters, check)
INTEGER, INTENT( IN ) :: N, niters
REAL, INTENT( OUT ) :: check
REAL, DIMENSION (N,N,N) :: P1, P2, P3, c
INTEGER iter
CALL acoustic3d_f90_setup(P1, P2, P3, c, N)
DO iter=1, niters
P3(2:N-1,2:N-1,2:N-1) = (2-6*c(2:N-1,2:N-1,2:N-1)) &
* P2(2:N-1,2:N-1,2:N-1) &
+ c(2:N-1,2:N-1,2:N-1)*(P2(1:N-2,2:N-1,2:N-1) + P2(3:N,2:N-1,2:N-1) &
+ P2(2:N-1,1:N-2,2:N-1)+P2(2:N-1,3:N,2:N-1) &
+ P2(2:N-1,2:N-1,1:N-2)+P2(2:N-1,2:N-1,3:N)) - P1(2:N-1,2:N-1,2:N-1)
P1 = P2
P2 = P3
END DO
check = P1(N/2,N/2,N/2)
RETURN
END
SUBROUTINE acoustic3d_f90_setup(P1, P2, P3, c, N)
INTEGER, INTENT( IN ) :: N
REAL, DIMENSION (N,N,N), INTENT( INOUT ) :: P1, P2, P3, c
INTEGER cavityLeft, cavityRight, cavityFront, cavityBack, cavityTop, &
cavityBottom, cavityTop2, cavityBottom2
REAL ci, cj, ck, s2
! Set the velocity field
c(1:N/2,:,:) = 0.05
c(N/2+1:N,:,:) = 0.3;
cavityLeft = 3 * N / 7.0
cavityRight = 4 * N / 7.0
cavityFront = 3 * N / 7.0
cavityBack = 4 * N / 7.0
cavityTop = 5 * N / 7.0
cavityBottom = 6 * N / 7.0
cavityTop2 = 1 * N / 7.0
cavityBottom2 = 2 * N / 7.0
c(cavityTop:cavityBottom,cavityLeft:cavityRight, &
cavityFront:cavityBack) = 0.02;
c(cavityTop2:cavityBottom2,cavityLeft:cavityRight, &
cavityFront:cavityBack) = 0.001;
! Initial pressure distribution
P1 = 0.0
P3 = 0.0
ci = N/2.0
cj = N/2.0
ck = N/2.0
s2 = 64.0 * 9.0 / ((N/2.0) ** 2)
DO k=1,N
DO j=1,N
DO i=1,N
P2(i,j,k) = exp(-((i-ci)**2+(j-cj)**2+(k-ck)**2) * s2)
END DO
END DO
END DO
CALL acoustic3d_f90_setup_check(P2, N)
CALL acoustic3d_f90_setup_check(c, N)
RETURN
END
SUBROUTINE acoustic3d_f90_setup_check(A, N)
INTEGER, INTENT( IN ) :: N
REAL, INTENT( IN ), DIMENSION(N,N,N) :: A
INTEGER i, j, k
REAL check
check = 0.0
DO k=1,N
DO j=1,N
DO i=1,N
check = check + A(i,j,k) * (i+N*j+N*N*k)
END DO
END DO
END DO
PRINT *, 'Array check: ', check
RETURN
END
|
