-
Notifications
You must be signed in to change notification settings - Fork 2
/
example2_pbe_2d_fv.f90
256 lines (208 loc) · 8.02 KB
/
example2_pbe_2d_fv.f90
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
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
program example_pbe_2d_fv
!! This program illustrates the application of modules 'hrweno' and 'mstvd' for solving a 2D
!! population balance equation (PBE):
!!```
!! d/dt u(x,t) + d/dx1 f1(u(x,t)) + d/dx2 f2(u(x,t)) = 0
!! with f1(u,t) = u
!! with f2(u,t) = u
!! and u(x,t=0) = ic(x)
!!```
!! The internal coordinates 'x=(x1,x2)' are discretized according to a finite-volume approach
!! and the time variable 't' is left continuous (method of lines). See example1 for notation.
!! In this particular example, we use the 3rd order 'mstvd' ode solver. The reconstruction
!! is done with the 5th order WENO scheme; to try other orders, we can change the parameter 'k'
!! below.
use, intrinsic :: iso_fortran_env, only: stderr => error_unit, stdout => output_unit
use hrweno_kinds, only: rk
use hrweno_tvdode, only: mstvd
use hrweno_weno, only: weno
use hrweno_fluxes, only: godunov
use hrweno_grids, only: grid1
use stdlib_strings, only: to_string
use omp_lib
implicit none
integer, parameter :: nc(2) = [250, 250], k = 3
real(rk) :: u(product(nc))
type(grid1) :: gx(2)
type(weno) :: myweno(2)
type(mstvd) :: ode
real(rk) :: dt, time, time_out, time_start, time_end
integer :: num_time_points, ii, jj
! OMP settings
! call omp_set_num_threads(min(2, omp_get_num_procs()))
! Define grids for x1 and x2
! In this example, we use linear grids, but any smooth grid can be used
call gx(1)%linear(xmin=0._rk, xmax=10._rk, ncells=nc(1))
call gx(2)%linear(xmin=0._rk, xmax=10._rk, ncells=nc(2))
! Init weno objects
myweno(1) = weno(ncells=nc(1), k=k, eps=1e-6_rk)
myweno(2) = weno(ncells=nc(2), k=k, eps=1e-6_rk)
! Open file where results will be stored
call output(1)
! Initial condition u(x,t=0)
do concurrent(ii=1:nc(1), jj=1:nc(2))
u((jj - 1)*nc(1) + ii) = ic([gx(1)%center(ii), gx(2)%center(jj)])
end do
! Call ODE time solver
ode = mstvd(rhs, size(u))
time_start = 0._rk
time_end = 5._rk
dt = 5e-3_rk
time = time_start
num_time_points = 100
do ii = 0, num_time_points
time_out = time_end*ii/num_time_points
call ode%integrate(u, time, time_out, dt)
call output(2)
end do
! End of simulation
call output(3)
contains
subroutine rhs(t, v, vdot)
!! This subroutine computes the *numerical approximation* to the right hand side of:
!!```
!! du(i,j,t)/dt = -1/dx1(i)*( f1(u(x1(i+1/2),x2(j),t)) - f1(u(x1(i-1/2),x2(j),t)) )
!! -1/dx2(j)*( f2(u(x1(i),x2(j+1/2),t)) - f2(u(x1(i),x2(j-1/2),t)) )
!!```
!! There are two main steps. First, we use the WENO scheme to obtain the reconstructed
!! values of 'u' at the left and right cell boundaries along each dimension. Second, we use
!! a suitable flux method (e.g., Godunov, Lax-Friedrichs) to compute the flux from the
!! reconstructed 'u' values.
real(rk), intent(in) :: t
!! time variable
real(rk), intent(in) :: v(:)
!! vector(N) with v(z,t) values
real(rk), intent(out) :: vdot(:)
!! vector(N) with v'(z,t) values
real(rk) :: vl1(nc(1)), vl2(nc(2)), vr1(nc(1)), vr2(nc(2))
real(rk) :: fedges1(0:nc(1), 0:nc(2)), fedges2(0:nc(2), 0:nc(1))
integer :: i, j
! Fluxes along x1 and x2 at interior cell boundaries
!$omp parallel
!$omp do private(vl1, vr1)
do j = 1, nc(2)
call myweno(1)%reconstruct(v(1+(j-1)*nc(1):j*nc(1)), vl1, vr1)
do i = 1, nc(1) - 1
fedges1(i, j) = godunov(flux1, vr1(i), vl1(i + 1), &
[gx(1)%right(i), gx(2)%center(j)], t)
end do
end do
!$omp end do nowait
!$omp do private(vl2, vr2)
do i = 1, nc(1)
call myweno(2)%reconstruct(v(i:i+(nc(2)-1)*nc(1):nc(1)), vl2, vr2)
do j = 1, nc(2) - 1
fedges2(j, i) = godunov(flux2, vr2(j), vl2(j + 1), &
[gx(1)%center(i), gx(2)%right(j)], t)
end do
end do
!$omp end do
!$omp end parallel
! Apply problem-specific flux constraints at domain boundaries
fedges1(0, :) = 0
fedges1(nc(1), :) = 0
fedges2(0, :) = 0
fedges2(nc(2), :) = 0
! Evaluate du/dt
do concurrent(i=1:nc(1), j=1:nc(2))
vdot((j - 1)*nc(1) + i) = &
- (fedges1(i, j) - fedges1(i - 1, j))/gx(1)%width(i) &
- (fedges2(j, i) - fedges2(j - 1, i))/gx(2)%width(j)
end do
end subroutine rhs
pure real(rk) function flux1(v, x, t)
!! Flux function along x1.
real(rk), intent(in) :: v
!! function v(z,t)
real(rk), intent(in) :: x(:)
!! vector of internal coordinates
real(rk), intent(in) :: t
!! time
flux1 = v !*x(1)**2
end function flux1
pure real(rk) function flux2(v, x, t)
!! Flux function along x2.
real(rk), intent(in) :: v
!! function v(z,t)
real(rk), intent(in) :: x(:)
!! vector of internal coordinates
real(rk), intent(in) :: t
!! time
flux2 = v !*x(1)*x(2)
end function flux2
pure real(rk) function ic(x)
!! Initial condition. Here we used rectangular pulse in both coordinates.
real(rk), intent(in) :: x(2)
!! vector of internal coordinates
if ((x(1) >= 1._rk .and. x(1) <= 3._rk) .and. (x(2) >= 1._rk .and. x(2) <= 3._rk)) then
ic = 1._rk
else
ic = 0._rk
end if
end function ic
subroutine output(action)
!! Auxiliary routine to save results to file.
integer, intent(in) :: action
!! parameter to select output action
character(*), parameter :: folder = ".\output\example2\"
real(rk) :: time_elapsed
integer, save :: funit_x(size(nc)) = 0, funit_u = 0
integer :: i, j
select case (action)
! Open files and write headers and grid
case (1)
write (stdout, '(1x, a)') "Running example2 ..."
write (stdout, '(1x, a, i3)') "Max # threads: ", omp_get_max_threads()
! Write grid
do concurrent(j=1:size(nc))
open (newunit=funit_x(j), file=folder//"x"//to_string(j)//".txt", &
status="replace", action="write", position="rewind")
write (funit_x(j), '(a5, 2(1x, a15))') "i", "x(i)", "dx(i)"
do i = 1, nc(j)
write (funit_x(j), '(i5, 2(1x, es15.5))') i, gx(j)%center(i), gx(j)%width(i)
end do
end do
call timer()
! Write header u
open (newunit=funit_u, file=folder//"u.txt", status="replace", &
action="write", position="rewind")
write (funit_u, '(a16)', advance="no") "t"
do i = 1, size(u)
write (funit_u, '(1x, a16)', advance="no") "u("//to_string(i)//")"
end do
write (funit_u, *) ""
! Write values u(x,t)
case (2)
write (funit_u, '(es16.5e3)', advance="no") time
do i = 1, size(u)
write (funit_u, '(1x, es16.5e3)', advance="no") u(i)
end do
write (funit_u, *) ""
! Close files
case (3)
do concurrent(j=1:size(nc))
close (funit_x(j))
end do
close (funit_u)
call timer(time_elapsed)
write (stdout, '(1x, a, 1x, f5.2)') "Elapsed time (s) :", time_elapsed
end select
end subroutine output
subroutine timer(res)
!! Quick and dirty timer
real(rk), intent(out), optional :: res
real(rk) :: res_
logical, save :: first_call = .true.
integer, save :: values_previous(8)
integer :: values_now(8), delta(8)
call date_and_time(values=values_now)
if (first_call) then
values_previous = values_now
first_call = .false.
end if
delta = values_now - values_previous
res_ = 1._rk*delta(3)*24*60**2 + 1._rk*delta(5)*60**2 + 1._rk*delta(6)*60 &
+ 1._rk*delta(7) + 1._rk*delta(8)/1000
if (present(res)) res = res_
end subroutine
end program example_pbe_2d_fv