-
Notifications
You must be signed in to change notification settings - Fork 0
/
1_5_12.jl
146 lines (108 loc) · 2.92 KB
/
1_5_12.jl
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
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
# solve 1.2.1 in the textbook
function driver()
xmin = 0
xmax = 1
tmax = .901
N = 11
# r = 0.5
delta_t = 0.02
nu = 1/6
ICFunc = ICSin
BCL = BCD2
BCR = BCD2
u, tmax_ret = solve(xmin, xmax, tmax, N, delta_t, nu, ICFunc, BCL, BCR)
# u = solve(xmin, xmax, tmax, N, delta_t, nu, ICFunc, BCL, BCR)
vals = [xmin, xmax, tmax_ret, nu, delta_t]
writedlm("counts.dat", vals)
writedlm("u.dat", u[2:end-1])
end
function solve(xmin, xmax, tmax, N, delta_t, nu, ICFunc::Function, BCL::Function, BCR::Function)
# xmin = minimum x coordinate
# xmax = maximum x coordinate
# tmax = maximum time value
# N = number of x points
# r = nu*delta_t/delta_x^2
# ICFunc = initial condition function with signature val = ICFunc(x)
# BCL = left boundary condition function with signature val = BCL(t)
# BCR = right boundary condition function
delta_x = (xmax - xmin)/(N-1)
#delta_t = (r*delta_x^2)/nu # nu*delta_t
r = nu*delta_t/(delta_x^2)
nStep = convert(Int, div(tmax, delta_t)) + 1
println("delta_x = ", delta_x)
println("delta_t = ", delta_t)
println("r = ", r)
println("nStep = ", nStep)
# allocate storage
n = N+2 # number of points in u
u_i = Array(Float64, n) # current timestep
u_i_1 = Array(Float64, n) # previous timestpe
# apply IC
# Not applying BC at initial condition
# hopefully IC and BC are consistent
for i=2:(n-1)
println("x = ", (i-2)*delta_x)
u_i_1[i] = ICFunc(xmin + (i-2)*delta_x)
end
flops = 0
time = @elapsed for tstep=2:nStep # loop over timesteps
# println("tstep = ", tstep)
# apply BC to u_i_1
BCL(u_i_1)
BCR(u_i_1, isleft=false)
println("u_i_1 =\n", u_i_1)
# calculate interior points
for i=3:(n-2)
u_k = u_i_1[i]
u_k_1 = u_i_1[i-1]
u_k_2 = u_i_1[i-2]
u_k_p1 = u_i_1[i+1]
u_k_p2 = u_i_1[i+2]
# u_i[i] = u_k + r*(u_k_p1 - 2*u_k + u_k_1)
u_i[i] = u_k + r*(-u_k_2/6 + 5*u_k_1/3 - 3*u_k + 5*u_k_p1/3 - u_k_p2/6)
end
flops += 12*(n-4)
# rebind names, using u_i as u_i_1, and reusing the array that was
# u_i_1 as u_i for the next timestep
tmp = u_i_1
u_i_1 = u_i
u_i = tmp
end
#apply BC to final time
BCL(u_i_1)
BCR(u_i_1, isleft=false)
# apply BC to fina
println("time = ", time)
flop_rate = 1e-6*flops/time # MFlops/seconds
println("flop rate = ", flop_rate, " MFlops/sec")
return u_i_1, delta_t*(nStep - 1) # return u_i_1
end
function BCZero(u; isleft=true)
# apply direchlet BC to boundary point, numerical BC =0 to ghost point
if isleft
u[2] = 0 # direchlet BC
u[1] = 0 # numerical BC
else
n = length(u)
u[n-1] = 0 # dirchlet BC
u[n] = 0 # numerical BC
end
return nothing
end
function BCD2(u; isleft=true)
# apply direchlet BC to boundary pointer, numerical u''=0 to ghost point
if isleft
u[2] = 0 # dirchlet BC
u[1] = 2*u[2] - u[3] # numerical BC
else
n = length(u)
u[n-1] = 0 # dirchlet BC
u[n] = 2*u[n-1] - u[n-2]
end
return nothing
end
function ICSin(x)
return sin(2*pi*x)
end
# run
driver()