-
Notifications
You must be signed in to change notification settings - Fork 2
/
example1_burgers_1d_fv.f90
190 lines (154 loc) · 6.18 KB
/
example1_burgers_1d_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
program example1_burgers_1d_fv
!! This program illustrates the application of modules 'hrweno' and 'tvdode' for solving a 1D
!! hyperbolic equation (Burger's equation):
!!```
!! d/dt u(x,t) = - d/dx f(u(x,t))
!! with f(u,t) = (u**2)/2
!! and u(x,t=0) = ic(x)
!!```
!! The spatial variable 'x' is discretized according to a finite-volume approach and the time
!! variable 't' is left continuous (method of lines), leading to:
!!```
!! du(i,t)/dt = -1/dx(i)*( f(u(x(i+1/2),t)) - f(u(x(i-1/2),t)) )
!!
!! ul=u(i-1/2)^+ ur=u(i+1/2)^-
!! --|-----(i-1)------|---------(i)----------|------(i+1)-----|--
!! x(i-1/2) x(i+1/2)
!! |<--- dx(i) --->|
!!```
!! In this particular example, we use the 3rd order 'rktvd' ode solver (we could equally well
!! employ the 'mstvd' solver). The reconstruction is done with the 5th order WENO scheme; to
!! try other orders, we can change the parameter 'k'.
use, intrinsic :: iso_fortran_env, only: stderr => error_unit, stdout => output_unit
use hrweno_kinds, only: rk
use hrweno_tvdode, only: rktvd
use hrweno_weno, only: weno
use hrweno_fluxes, only: godunov, lax_friedrichs
use hrweno_grids, only: grid1
use stdlib_strings, only: to_string
implicit none
integer, parameter :: nc = 100
real(rk) :: u(nc)
real(rk) :: dt, time, time_out, time_start, time_end
integer :: num_time_points, ii
type(grid1) :: gx
type(weno) :: myweno
type(rktvd) :: ode
! Define the spatial grid
! In this example, we use a linear grid, but any smooth grid can be used
call gx%linear(xmin=-5._rk, xmax=5._rk, ncells=nc)
! Init weno object
myweno = weno(ncells=nc, k=3, eps=1e-6_rk)
! Open file where results will be stored
call output(1)
! Initial condition u(x,t=0)
u = ic(gx%center)
! Call ODE time solver
ode = rktvd(rhs, nc, order=3)
time_start = 0._rk
time_end = 12._rk
dt = 1e-2_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
pure subroutine rhs(t, v, vdot)
!! This subroutine computes the *numerical approximation* to the right hand side of:
!!```
!! du(i,t)/dt = -1/dx(i)*( f(u(x(i+1/2),t)) - f(u(x(i-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. Note that, in general, because of
!! discontinuities, \( u_{i+1/2}^+ \neq u_{(i+1)+1/2}^- \). 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(x,t) values
real(rk), intent(out) :: vdot(:)
!! vector(N) with v'(x,t) values
real(rk) :: fedges(0:nc), vl(nc), vr(nc)
integer :: i
! Get reconstructed values at cell boundaries
call myweno%reconstruct(v, vl, vr)
! Fluxes at interior cell boundaries
! One can use the Lax-Friedrichs or the Godunov method
do concurrent(i=1:nc - 1)
!fedges(i) = lax_friedrichs(flux, vr(i), vl(i+1), gx%r(i), t, alpha = 1._rk)
fedges(i) = godunov(flux, vr(i), vl(i + 1), [gx%right(i)], t)
end do
! Apply problem-specific flux constraints at domain boundaries
fedges(0) = fedges(1)
fedges(nc) = fedges(nc - 1)
! Evaluate du/dt
vdot = -(fedges(1:) - fedges(:nc - 1))/gx%width
end subroutine rhs
pure real(rk) function flux(v, x, t)
!! Flux function. Here we define the flux corresponding to Burger's equation.
real(rk), intent(in) :: v
!! function v(x,t)
real(rk), intent(in) :: x(:)
!! spatial variable
real(rk), intent(in) :: t
!! time variable
flux = (v**2)/2
end function flux
elemental real(rk) function ic(x)
!! Initial condition. Here we used a limited linear profile.
real(rk), intent(in) :: x
!! spatial variable
real(rk), parameter :: xa = -4._rk, xb = 2._rk, va = 1._rk, vb = -0.5_rk
ic = va + (vb - va)/(xb - xa)*(x - xa)
ic = max(min(ic, va), vb)
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\example1\"
real(rk), save :: cpu_start = 0._rk, cpu_end = 0._rk
integer, save :: funit_x = 0, funit_u = 0
integer :: i
select case (action)
! Open files and write headers and grid
case (1)
write (stdout, '(1x, a)') "Running example1 ..."
call cpu_time(cpu_start)
! Write grid
open (newunit=funit_x, file=folder//"x.txt", status="replace", &
action="write", position="rewind")
write (funit_x, '(a5, 2(1x, a15))') "i", "x(i)", "dx(i)"
do i = 1, nc
write (funit_x, '(i5, 2(1x, es15.5))') i, gx%center(i), gx%width(i)
end do
! 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, nc
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, nc
write (funit_u, '(1x, es16.5e3)', advance="no") u(i)
end do
write (funit_u, *) ""
! Close files
case (3)
close (funit_x)
close (funit_u)
call cpu_time(cpu_end)
write (stdout, '(1x, a, 1x, f6.1)') "Elapsed time (ms) :", 1e3_rk*(cpu_end - cpu_start)
end select
end subroutine output
end program example1_burgers_1d_fv