MPI-AMRVAC 3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
Loading...
Searching...
No Matches
mod_weno.t
Go to the documentation of this file.
1module mod_weno
2 ! All kinds of (W)ENO schemes
3 !
4 ! 2019.9.19 WENO(-JS)5 transplant from the BHAC code;
5 ! 2019.9.20 WENO3;
6 ! 2019.9.21 WENO-Z5;
7 ! 2019.9.22 WENO-Z+5 transplant from the BHAC code;
8 ! 2019.10.30 WENO(-JS)7;
9 ! 2019.10.31 MPWENO7;
10 ! 2019.11.1 exENO7;
11 ! 2019.11.7 clean up the code, comment out the interpolation variation;
12 ! 2019.12.9 WENO-YC3;
13 ! 2020.1.15 new WENO5 limiter WENO5NM for stretched grid.
14 ! 2020.4.15 WENO5-CU6: hybrid sixth-order linear & WENO5
15 ! 2021.6.12 generic treatment for fixing unphysical reconstructions
16 ! 2022.10.21 remove exENO7
17
18 implicit none
19 private
20
21 double precision, parameter :: weno_eps_machine = 1.0d-18
22 double precision, parameter :: weno_dx_exp = 2.0d0/3.0d0
23
24 public :: weno3limiter
25 public :: weno5limiter
26 public :: weno5nmlimiter
27 public :: weno5limiterl
28 public :: weno5nmlimiterl
29 public :: weno5limiterr
30 public :: weno5nmlimiterr
31 public :: teno5adlimiter
32 public :: weno5cu6limiter
33 public :: weno7limiter
34 public :: fix_limiter
35 public :: fix_limiter1
36 public :: fix_onelimiter
37 public :: fix_onelimiter1
38
39contains
40
41 subroutine fix_onelimiter(ixI^L,iL^L,wCin,wCout)
42 use mod_global_parameters
43 use mod_physics, only: phys_check_w
44
45 integer, intent(in) :: ixi^l, il^l
46 double precision, intent(in) :: wcin(ixi^s,1:nw)
47 double precision, intent(inout) :: wcout(ixi^s,1:nw)
48
49 integer :: iw
50 logical :: flagc(ixi^s,1:nw), flag(ixi^s)
51
52 ! When limiter not TVD, negative pressures or densities could result.
53 ! Fall back to flat interpolation
54 ! flagC(*,iw) indicates failed state (T when failed)
55 ! assumes wCin contains primitive variables
56 call phys_check_w(.true.,ixi^l,il^l,wcin,flagc)
57
58 ! collect all failures
59 flag(il^s)=.false.
60 do iw=1,nw_recon
61 flag(il^s)=flag(il^s).or.(flagc(il^s,iw))
62 end do
63 ! only use WENO reconstructions when no field failed
64 ! in other places: do not modify the initial state in wCout
65 do iw=1,nw_recon
66 where (flag(il^s) .eqv. .false.)
67 wcout(il^s,iw)=wcin(il^s,iw)
68 end where
69 enddo
70
71 end subroutine fix_onelimiter
72
73 subroutine fix_onelimiter1(ixI^L,iL^L,wCin,wCout)
74 use mod_global_parameters
75 use mod_physics, only: phys_check_w
76
77 integer, intent(in) :: ixi^l, il^l
78 double precision, intent(in) :: wcin(ixi^s,1:nw)
79 double precision, intent(inout) :: wcout(ixi^s,1:nw)
80
81 integer :: iw
82 logical :: flagc(ixi^s,1:nw)
83
84 ! When limiter not TVD, negative pressures or densities could result.
85 ! Fall back to flat interpolation
86 ! flagC(*,iw) indicates failed state (T when failed)
87 ! assumes wCin contains primitive variables
88 call phys_check_w(.true.,ixi^l,il^l,wcin,flagc)
89
90 ! only use WENO reconstructions when no field failed
91 ! in other places: do not modify the initial state in wCout
92 do iw=1,nw_recon
93 where (flagc(il^s,iw) .eqv. .false.)
94 wcout(il^s,iw)=wcin(il^s,iw)
95 end where
96 enddo
97
98 end subroutine fix_onelimiter1
99
100 subroutine fix_limiter(ixI^L,iL^L,wLCin,wRCin,wLCout,wRCout)
101 use mod_global_parameters
102 use mod_physics, only: phys_check_w
103
104 integer, intent(in) :: ixi^l, il^l
105 double precision, intent(in) :: wrcin(ixi^s,1:nw),wlcin(ixi^s,1:nw)
106 double precision, intent(inout) :: wrcout(ixi^s,1:nw),wlcout(ixi^s,1:nw)
107
108 integer :: iw
109 logical :: flagl(ixi^s,1:nw), flagr(ixi^s,1:nw), flag(ixi^s)
110
111 ! When limiter not TVD, negative pressures or densities could result.
112 ! Fall back to flat interpolation
113 ! flagL(*,iw) indicates failed L state (T when failed)
114 ! flagR(*,iw) indicates failed R state (T when failed)
115 ! assumes wLCin and wRCin contain primitive variables
116 call phys_check_w(.true.,ixi^l,il^l,wlcin,flagl)
117 call phys_check_w(.true.,ixi^l,il^l,wrcin,flagr)
118
119 ! collect all failures
120 flag(il^s)=.false.
121 do iw=1,nw_recon
122 flag(il^s)=flag(il^s).or.(flagl(il^s,iw).or.flagr(il^s,iw))
123 end do
124 ! only use WENO reconstructions L and R when no neighbour field failed
125 ! in other places: do not modify the initial state in wLCout/wRCout
126 do iw=1,nw_recon
127 where (flag(il^s) .eqv. .false.)
128 wlcout(il^s,iw)=wlcin(il^s,iw)
129 wrcout(il^s,iw)=wrcin(il^s,iw)
130 end where
131 enddo
132
133 end subroutine fix_limiter
134
135 subroutine fix_limiter1(ixI^L,iL^L,wLCin,wRCin,wLCout,wRCout)
136 use mod_global_parameters
137 use mod_physics, only: phys_check_w
138
139 integer, intent(in) :: ixi^l, il^l
140 double precision, intent(in) :: wrcin(ixi^s,1:nw),wlcin(ixi^s,1:nw)
141 double precision, intent(inout) :: wrcout(ixi^s,1:nw),wlcout(ixi^s,1:nw)
142
143 integer :: iw
144 logical :: flagl(ixi^s,1:nw), flagr(ixi^s,1:nw)
145
146 ! When limiter not TVD, negative pressures or densities could result.
147 ! Fall back to flat interpolation
148 ! flagL(*,iw) indicates failed L state (T when failed)
149 ! flagR(*,iw) indicates failed R state (T when failed)
150 ! assumes wLCin and wRCin contain primitive variables
151 call phys_check_w(.true.,ixi^l,il^l,wlcin,flagl)
152 call phys_check_w(.true.,ixi^l,il^l,wrcin,flagr)
153
154 ! only use WENO reconstructions L and R when no neighbour field failed
155 ! in other places: do not modify the initial state in wLCout/wRCout
156 do iw=1,nw_recon
157 where ((flagl(il^s,iw) .eqv. .false.) .and. (flagr(il^s,iw) .eqv. .false.))
158 wlcout(il^s,iw)=wlcin(il^s,iw)
159 wrcout(il^s,iw)=wrcin(il^s,iw)
160 end where
161 enddo
162
163 end subroutine fix_limiter1
164
165 subroutine weno3limiter(ixI^L,iL^L,idims,dxdim,w,wLC,wRC,var)
166 use mod_global_parameters
167
168 integer, intent(in) :: ixi^l, il^l, idims, var
169 double precision, intent(in) :: dxdim
170 double precision, intent(in) :: w(ixi^s,1:nw)
171 double precision, intent(inout) :: wrc(ixi^s,1:nw),wlc(ixi^s,1:nw)
172 !> local
173 double precision :: f_array(ixi^s,1:nw,2), d_array(2)
174 double precision :: beta(ixi^s,1:nw,2),tau(ixi^s,1:nw)
175 double precision :: u1_coeff(2), u2_coeff(2)
176 double precision :: alpha_array(ixi^s,1:nw,2), alpha_sum(ixi^s,1:nw), flux(ixi^s,1:nw)
177 double precision, dimension(ixI^S,1:nw) :: wrctmp, wlctmp
178 integer :: ilm^l, ilp^l, ilpp^l, i
179
180 ! iL^L holds the indices of interfaces to reconstruct to. Convention is that a center index holds the _right-side_ interface.
181 ilm^l=il^l-kr(idims,^d);
182 ilp^l=il^l+kr(idims,^d);
183 ilpp^l=ilp^l+kr(idims,^d);
184 d_array(1:2) = (/ 1.0d0/3.0d0, 2.0d0/3.0d0 /)
185 u1_coeff(1:2) = (/ -1.d0/2.d0, 3.d0/2.d0 /)
186 u2_coeff(1:2) = (/ 1.d0/2.d0, 1.d0/2.d0 /)
187
188 !> left side
189 f_array(il^s,1:nw_recon,1) = u1_coeff(1) * w(ilm^s,1:nw_recon) + u1_coeff(2) * w(il^s,1:nw_recon)
190 f_array(il^s,1:nw_recon,2) = u2_coeff(1) * w(il^s,1:nw_recon) + u2_coeff(2) * w(ilp^s,1:nw_recon)
191
192 beta(il^s,1:nw_recon,1) = (w(il^s,1:nw_recon) - w(ilm^s,1:nw_recon))**2
193 beta(il^s,1:nw_recon,2) = (w(ilp^s,1:nw_recon) - w(il^s,1:nw_recon))**2
194
195 alpha_sum(il^s,1:nw_recon) = 0.0d0
196 select case(var)
197 case(1)
198 do i = 1,2
199 alpha_array(il^s,1:nw_recon,i) = d_array(i)/(beta(il^s,1:nw_recon,i) + weno_eps_machine)**2
200 alpha_sum(il^s,1:nw_recon) = alpha_sum(il^s,1:nw_recon) + alpha_array(il^s,1:nw_recon,i)
201 end do
202 case(2)
203 tau(il^s,1:nw_recon) = abs(beta(il^s,1:nw_recon,2) - beta(il^s,1:nw_recon,1))
204 do i = 1,2
205 alpha_array(il^s,1:nw_recon,i) = d_array(i) * (1.d0 + (tau(il^s,1:nw_recon) / &
206 (beta(il^s,1:nw_recon,i) + dxdim**2)))
207 alpha_sum(il^s,1:nw_recon) = alpha_sum(il^s,1:nw_recon) + alpha_array(il^s,1:nw_recon,i)
208 end do
209 end select
210
211 flux(il^s,1:nw_recon) = 0.0d0
212 do i = 1,2
213 flux(il^s,1:nw_recon) = flux(il^s,1:nw_recon) + f_array(il^s,1:nw_recon,i) * alpha_array(il^s,1:nw_recon,i)/(alpha_sum(il^s,1:nw_recon))
214 end do
215
216 !> left value at right interface
217 wlctmp(il^s,1:nw_recon) = flux(il^s,1:nw_recon)
218
219 !> right side
220 f_array(il^s,1:nw_recon,1) = u1_coeff(1) * w(ilpp^s,1:nw_recon) + u1_coeff(2) * w(ilp^s,1:nw_recon)
221 f_array(il^s,1:nw_recon,2) = u2_coeff(1) * w(ilp^s,1:nw_recon) + u2_coeff(2) * w(il^s,1:nw_recon)
222
223 beta(il^s,1:nw_recon,1) = (w(ilpp^s,1:nw_recon) - w(ilp^s,1:nw_recon))**2
224 beta(il^s,1:nw_recon,2) = (w(ilp^s,1:nw_recon) - w(il^s,1:nw_recon))**2
225
226 alpha_sum(il^s,1:nw_recon) = 0.0d0
227 select case(var)
228 case(1)
229 do i = 1,2
230 alpha_array(il^s,1:nw_recon,i) = d_array(i)/(beta(il^s,1:nw_recon,i) + weno_eps_machine)**2
231 alpha_sum(il^s,1:nw_recon) = alpha_sum(il^s,1:nw_recon) + alpha_array(il^s,1:nw_recon,i)
232 end do
233 case(2)
234 tau(il^s,1:nw_recon) = abs(beta(il^s,1:nw_recon,2) - beta(il^s,1:nw_recon,1))
235 do i = 1,2
236 alpha_array(il^s,1:nw_recon,i) = d_array(i) * (1.d0 + (tau(il^s,1:nw_recon) / &
237 (beta(il^s,1:nw_recon,i) + dxdim**2)))
238 alpha_sum(il^s,1:nw_recon) = alpha_sum(il^s,1:nw_recon) + alpha_array(il^s,1:nw_recon,i)
239 end do
240 end select
241
242 flux(il^s,1:nw_recon) = 0.0d0
243 do i = 1,2
244 flux(il^s,1:nw_recon) = flux(il^s,1:nw_recon) + f_array(il^s,1:nw_recon,i) * alpha_array(il^s,1:nw_recon,i)/(alpha_sum(il^s,1:nw_recon))
245 end do
246
247 !> right value at right interface
248 wrctmp(il^s,1:nw_recon) = flux(il^s,1:nw_recon)
249
250 call fix_limiter(ixi^l,il^l,wlctmp,wrctmp,wlc,wrc)
251
252 end subroutine weno3limiter
253
254 subroutine weno5limiter(ixI^L,iL^L,idims,dxdim,w,wLC,wRC,var)
255 use mod_global_parameters
256
257 integer, intent(in) :: ixi^l, il^l, idims, var
258 double precision, intent(in) :: dxdim
259 double precision, intent(in) :: w(ixi^s,1:nw)
260 double precision, intent(inout) :: wrc(ixi^s,1:nw),wlc(ixi^s,1:nw)
261 !> local
262 double precision :: f_array(ixi^s,1:nw,3), d_array(3)
263 double precision :: beta(ixi^s,1:nw,3), beta_coeff(2)
264 double precision :: tau(ixi^s,1:nw), tmp(ixi^s,1:nw)
265 double precision :: u1_coeff(3), u2_coeff(3), u3_coeff(3)
266 double precision :: alpha_array(ixi^s,1:nw,3), alpha_sum(ixi^s,1:nw), flux(ixi^s,1:nw)
267 double precision :: lambda
268 double precision, dimension(ixI^S,1:nw) :: wrctmp, wlctmp
269 integer :: ilm^l, ilmm^l, ilp^l, ilpp^l, ilppp^l, i
270
271 ilm^l=il^l-kr(idims,^d);
272 ilmm^l=ilm^l-kr(idims,^d);
273 ilp^l=il^l+kr(idims,^d);
274 ilpp^l=ilp^l+kr(idims,^d);
275 ilppp^l=ilpp^l+kr(idims,^d);
276 lambda = dxdim**weno_dx_exp
277 beta_coeff(1:2) = (/ 1.0833333333333333d0, 0.25d0/)
278! reconstruction variation
279 d_array(1:3) = (/ 1.0d0/10.0d0, 3.0d0/5.0d0, 3.0d0/10.0d0 /)
280 u1_coeff(1:3) = (/ 1.d0/3.d0, -7.d0/6.d0, 11.d0/6.d0 /)
281 u2_coeff(1:3) = (/ -1.d0/6.d0, 5.d0/6.d0, 1.d0/3.d0 /)
282 u3_coeff(1:3) = (/ 1.d0/3.d0, 5.d0/6.d0, -1.d0/6.d0 /)
283! interpolation variation
284! d_array(1:3) = (/ 1.0d0/16.0d0, 10.0d0/16.0d0, 5.0d0/16.0d0 /)
285! u1_coeff(1:3) = (/ 3.d0/8.d0, -10.d0/8.d0, 15.d0/8.d0 /)
286! u2_coeff(1:3) = (/ -1.d0/8.d0, 6.d0/8.d0, 3.d0/8.d0 /)
287! u3_coeff(1:3) = (/ 3.d0/8.d0, 6.d0/8.d0, -1.d0/8.d0 /)
288
289 !> left side
290 f_array(il^s,1:nw_recon,1) = u1_coeff(1) * w(ilmm^s,1:nw_recon) + u1_coeff(2) * w(ilm^s,1:nw_recon) + u1_coeff(3) * w(il^s,1:nw_recon)
291 f_array(il^s,1:nw_recon,2) = u2_coeff(1) * w(ilm^s,1:nw_recon) + u2_coeff(2) * w(il^s,1:nw_recon) + u2_coeff(3) * w(ilp^s,1:nw_recon)
292 f_array(il^s,1:nw_recon,3) = u3_coeff(1) * w(il^s,1:nw_recon) + u3_coeff(2) * w(ilp^s,1:nw_recon) + u3_coeff(3) * w(ilpp^s,1:nw_recon)
293
294 beta(il^s,1:nw_recon,1) = beta_coeff(1) * (w(ilmm^s,1:nw_recon) + w(il^s,1:nw_recon) - 2.0d0*w(ilm^s,1:nw_recon))**2 &
295 + beta_coeff(2) * (w(ilmm^s,1:nw_recon) - 4.0d0 * w(ilm^s,1:nw_recon) + 3.0d0*w(il^s,1:nw_recon))**2
296 beta(il^s,1:nw_recon,2) = beta_coeff(1) * (w(ilm^s,1:nw_recon) + w(ilp^s,1:nw_recon) - 2.0d0 * w(il^s,1:nw_recon))**2 &
297 + beta_coeff(2) * (w(ilm^s,1:nw_recon) - w(ilp^s,1:nw_recon))**2
298 beta(il^s,1:nw_recon,3) = beta_coeff(1) * (w(il^s,1:nw_recon) + w(ilpp^s,1:nw_recon) - 2.0d0 * w(ilp^s,1:nw_recon))**2 &
299 + beta_coeff(2) * (3.0d0 * w(il^s, 1:nw_recon) - 4.0d0 * w(ilp^s,1:nw_recon) + w(ilpp^s,1:nw_recon))**2
300
301 select case(var)
302 ! case1 for wenojs, case2 for wenoz, case3 for wenoz+
303 case(1)
304 do i = 1,3
305 alpha_array(il^s,1:nw_recon,i) = d_array(i)/(beta(il^s,1:nw_recon,i) + weno_eps_machine)**2
306 end do
307 case(2)
308 tau(il^s,1:nw_recon) = abs(beta(il^s,1:nw_recon,1) - beta(il^s,1:nw_recon,3))
309 do i = 1,3
310 alpha_array(il^s,1:nw_recon,i) = d_array(i) * (1.d0 + (tau(il^s,1:nw_recon) / &
311 (beta(il^s,1:nw_recon,i) + weno_eps_machine))**2)
312 end do
313 case(3)
314 tau(il^s,1:nw_recon) = abs(beta(il^s,1:nw_recon,1) - beta(il^s,1:nw_recon,3))
315 do i = 1,3
316 tmp(il^s,1:nw_recon) = (tau(il^s,1:nw_recon) + weno_eps_machine) / (beta(il^s,1:nw_recon,i) + weno_eps_machine)
317 alpha_array(il^s,1:nw_recon,i) = d_array(i) * (1.0d0 + tmp(il^s,1:nw_recon)**2 + lambda/tmp(il^s,1:nw_recon))
318 end do
319 end select
320
321 alpha_sum(il^s,1:nw_recon) = 0.0d0
322 do i = 1,3
323 alpha_sum(il^s,1:nw_recon) = alpha_sum(il^s,1:nw_recon) + alpha_array(il^s,1:nw_recon,i)
324 end do
325 flux(il^s,1:nw_recon) = 0.0d0
326 do i = 1,3
327 flux(il^s,1:nw_recon) = flux(il^s,1:nw_recon) + f_array(il^s,1:nw_recon,i) &
328 *alpha_array(il^s,1:nw_recon,i)/(alpha_sum(il^s,1:nw_recon))
329 end do
330 wlctmp(il^s,1:nw_recon) = flux(il^s,1:nw_recon)
331
332 !> right side
333 f_array(il^s,1:nw_recon,1) = u1_coeff(1) * w(ilppp^s,1:nw_recon) + u1_coeff(2) * w(ilpp^s,1:nw_recon) + u1_coeff(3) * w(ilp^s,1:nw_recon)
334 f_array(il^s,1:nw_recon,2) = u2_coeff(1) * w(ilpp^s,1:nw_recon) + u2_coeff(2) * w(ilp^s,1:nw_recon) + u2_coeff(3) * w(il^s,1:nw_recon)
335 f_array(il^s,1:nw_recon,3) = u3_coeff(1) * w(ilp^s,1:nw_recon) + u3_coeff(2) * w(il^s,1:nw_recon) + u3_coeff(3) * w(ilm^s,1:nw_recon)
336
337 beta(il^s,1:nw_recon,1) = beta_coeff(1) * (w(ilppp^s,1:nw_recon) + w(ilp^s,1:nw_recon) - 2.0d0*w(ilpp^s,1:nw_recon))**2 &
338 + beta_coeff(2) * (w(ilppp^s,1:nw_recon) - 4.0d0 * w(ilpp^s,1:nw_recon) + 3.0d0*w(ilp^s,1:nw_recon))**2
339 beta(il^s,1:nw_recon,2) = beta_coeff(1) * (w(ilpp^s,1:nw_recon) + w(il^s,1:nw_recon) - 2.0d0 * w(ilp^s,1:nw_recon))**2 &
340 + beta_coeff(2) * (w(ilpp^s,1:nw_recon) - w(il^s,1:nw_recon))**2
341 beta(il^s,1:nw_recon,3) = beta_coeff(1) * (w(ilp^s,1:nw_recon) + w(ilm^s,1:nw_recon) - 2.0d0 * w(il^s,1:nw_recon))**2 &
342 + beta_coeff(2) * (3.0d0 * w(ilp^s, 1:nw_recon) - 4.0d0 * w(il^s,1:nw_recon) + w(ilm^s,1:nw_recon))**2
343
344 select case(var)
345 case(1)
346 do i = 1,3
347 alpha_array(il^s,1:nw_recon,i) = d_array(i)/(beta(il^s,1:nw_recon,i) + weno_eps_machine)**2
348 end do
349 case(2)
350 tau(il^s,1:nw_recon) = abs(beta(il^s,1:nw_recon,1) - beta(il^s,1:nw_recon,3))
351 do i = 1,3
352 alpha_array(il^s,1:nw_recon,i) = d_array(i) * (1.d0 + (tau(il^s,1:nw_recon) / &
353 (beta(il^s,1:nw_recon,i) + weno_eps_machine))**2)
354 end do
355 case(3)
356 tau(il^s,1:nw_recon) = abs(beta(il^s,1:nw_recon,1) - beta(il^s,1:nw_recon,3))
357 do i = 1,3
358 tmp(il^s,1:nw_recon) = (tau(il^s,1:nw_recon) + weno_eps_machine) / (beta(il^s,1:nw_recon,i) + weno_eps_machine)
359 alpha_array(il^s,1:nw_recon,i) = d_array(i) * (1.0d0 + tmp(il^s,1:nw_recon)**2 + lambda/tmp(il^s,1:nw_recon))
360 end do
361 end select
362
363 alpha_sum(il^s,1:nw_recon) = 0.0d0
364 ! do nothing, normal case
365 do i = 1,3
366 alpha_sum(il^s,1:nw_recon) = alpha_sum(il^s,1:nw_recon) + alpha_array(il^s,1:nw_recon,i)
367 end do
368 flux(il^s,1:nw_recon) = 0.0d0
369 do i = 1,3
370 flux(il^s,1:nw_recon) = flux(il^s,1:nw_recon) + f_array(il^s,1:nw_recon,i) &
371 *alpha_array(il^s,1:nw_recon,i)/(alpha_sum(il^s,1:nw_recon))
372 end do
373 wrctmp(il^s,1:nw_recon) = flux(il^s,1:nw_recon)
374
375 call fix_limiter(ixi^l,il^l,wlctmp,wrctmp,wlc,wrc)
376
377 end subroutine weno5limiter
378
379 subroutine teno5adlimiter(ixI^L,iL^L,idims,dxdim,w,wLC,wRC)
380 use mod_global_parameters
381 integer, intent(in) :: ixi^l, il^l, idims
382 double precision, intent(in) :: dxdim
383 double precision, intent(in) :: w(ixi^s,1:nw)
384 double precision, intent(inout) :: wrc(ixi^s,1:nw),wlc(ixi^s,1:nw)
385 !> local
386 double precision :: f_array(ixi^s,1:nw,3), d_array(3)
387 double precision :: beta(ixi^s,1:nw,3), beta_coeff(2)
388 double precision :: u1_coeff(3), u2_coeff(3), u3_coeff(3)
389 double precision :: gam_sum(ixi^s,1:nw),tau(ixi^s,1:nw),delta_sum(ixi^s,1:nw)
390 double precision :: gam(ixi^s,1:nw,3), kai(ixi^s,1:nw,3), delta(ixi^s,1:nw,3)
391 double precision :: flux(ixi^s,1:nw), kai1(ixi^s,1:nw,3), theta(ixi^s,1:nw)
392 double precision, dimension(ixI^S,1:nw) :: wrctmp, wlctmp
393 integer :: marray(ixi^s,1:nw)
394 integer :: ilm^l, ilmm^l, ilp^l, ilpp^l, ilppp^l, i
395
396 ilm^l=il^l-kr(idims,^d);
397 ilmm^l=ilm^l-kr(idims,^d);
398 ilp^l=il^l+kr(idims,^d);
399 ilpp^l=ilp^l+kr(idims,^d);
400 ilppp^l=ilpp^l+kr(idims,^d);
401 beta_coeff(1:2) = (/ 1.0833333333333333d0, 0.25d0/)
402 d_array(1:3) = (/ 1.0d0/10.0d0, 3.0d0/5.0d0, 3.0d0/10.0d0 /)
403 u1_coeff(1:3) = (/ 1.d0/3.d0, -7.d0/6.d0, 11.d0/6.d0 /)
404 u2_coeff(1:3) = (/ -1.d0/6.d0, 5.d0/6.d0, 1.d0/3.d0 /)
405 u3_coeff(1:3) = (/ 1.d0/3.d0, 5.d0/6.d0, -1.d0/6.d0 /)
406
407 !> left side
408 f_array(il^s,1:nw_recon,1) = u1_coeff(1) * w(ilmm^s,1:nw_recon) + u1_coeff(2) * w(ilm^s,1:nw_recon) + u1_coeff(3) * w(il^s,1:nw_recon)
409 f_array(il^s,1:nw_recon,2) = u2_coeff(1) * w(ilm^s,1:nw_recon) + u2_coeff(2) * w(il^s,1:nw_recon) + u2_coeff(3) * w(ilp^s,1:nw_recon)
410 f_array(il^s,1:nw_recon,3) = u3_coeff(1) * w(il^s,1:nw_recon) + u3_coeff(2) * w(ilp^s,1:nw_recon) + u3_coeff(3) * w(ilpp^s,1:nw_recon)
411
412 beta(il^s,1:nw_recon,1) = beta_coeff(1) * (w(ilmm^s,1:nw_recon) + w(il^s,1:nw_recon) - 2.0d0*w(ilm^s,1:nw_recon))**2 &
413 + beta_coeff(2) * (w(ilmm^s,1:nw_recon) - 4.0d0 * w(ilm^s,1:nw_recon) + 3.0d0*w(il^s,1:nw_recon))**2
414 beta(il^s,1:nw_recon,2) = beta_coeff(1) * (w(ilm^s,1:nw_recon) + w(ilp^s,1:nw_recon) - 2.0d0 * w(il^s,1:nw_recon))**2 &
415 + beta_coeff(2) * (w(ilm^s,1:nw_recon) - w(ilp^s,1:nw_recon))**2
416 beta(il^s,1:nw_recon,3) = beta_coeff(1) * (w(il^s,1:nw_recon) + w(ilpp^s,1:nw_recon) - 2.0d0 * w(ilp^s,1:nw_recon))**2 &
417 + beta_coeff(2) * (3.0d0 * w(il^s, 1:nw_recon) - 4.0d0 * w(ilp^s,1:nw_recon) + w(ilpp^s,1:nw_recon))**2
418
419 gam_sum(il^s,1:nw_recon) = 0.0d0
420 tau(il^s,1:nw_recon) = abs(beta(il^s,1:nw_recon,1) - beta(il^s,1:nw_recon,3))
421 do i = 1,3
422 kai1(il^s,1:nw_recon,i) = (tau(il^s,1:nw_recon) / (beta(il^s,1:nw_recon,i) + weno_eps_machine))
423 gam(il^s,1:nw_recon,i) = (1.d0 + kai1(il^s,1:nw_recon,i))**2
424 gam_sum(il^s,1:nw_recon) = gam_sum(il^s,1:nw_recon) + gam(il^s,1:nw_recon,i)
425 end do
426 theta(il^s,1:nw_recon) = one / (one + maxval(kai1(il^s,1:nw_recon,1:3)/10.d0, dim=ndim+2))
427 marray(il^s,1:nw_recon)=-floor(4.d0 + theta(il^s,1:nw_recon)*6.d0)
428 do i = 1,3
429 kai(il^s,1:nw_recon,i) = gam(il^s,1:nw_recon,i) / gam_sum(il^s,1:nw_recon)
430 where(kai(il^s,1:nw_recon,i) .lt. 10**marray(il^s,1:nw_recon))
431 delta(il^s,1:nw_recon,i)=zero
432 else where
433 delta(il^s,1:nw_recon,i)=one
434 end where
435 end do
436 delta_sum=zero
437 do i = 1,3
438 delta_sum(il^s,1:nw_recon)=delta_sum(il^s,1:nw_recon)+delta(il^s,1:nw_recon,i)*d_array(i)
439 end do
440 flux(il^s,1:nw_recon)=0.0d0
441 do i = 1,3
442 flux(il^s,1:nw_recon)=flux(il^s,1:nw_recon)+f_array(il^s,1:nw_recon,i)*delta(il^s,1:nw_recon,i)*d_array(i)/(delta_sum(il^s,1:nw_recon))
443 end do
444
445 !> left value at right interface
446 wlctmp(il^s,1:nw_recon) = flux(il^s,1:nw_recon)
447
448 !> right side
449 f_array(il^s,1:nw_recon,1) = u1_coeff(1) * w(ilppp^s,1:nw_recon) + u1_coeff(2) * w(ilpp^s,1:nw_recon) + u1_coeff(3) * w(ilp^s,1:nw_recon)
450 f_array(il^s,1:nw_recon,2) = u2_coeff(1) * w(ilpp^s,1:nw_recon) + u2_coeff(2) * w(ilp^s,1:nw_recon) + u2_coeff(3) * w(il^s,1:nw_recon)
451 f_array(il^s,1:nw_recon,3) = u3_coeff(1) * w(ilp^s,1:nw_recon) + u3_coeff(2) * w(il^s,1:nw_recon) + u3_coeff(3) * w(ilm^s,1:nw_recon)
452
453 beta(il^s,1:nw_recon,1) = beta_coeff(1) * (w(ilppp^s,1:nw_recon) + w(ilp^s,1:nw_recon) - 2.0d0*w(ilpp^s,1:nw_recon))**2 &
454 + beta_coeff(2) * (w(ilppp^s,1:nw_recon) - 4.0d0 * w(ilpp^s,1:nw_recon) + 3.0d0*w(ilp^s,1:nw_recon))**2
455 beta(il^s,1:nw_recon,2) = beta_coeff(1) * (w(ilpp^s,1:nw_recon) + w(il^s,1:nw_recon) - 2.0d0 * w(ilp^s,1:nw_recon))**2 &
456 + beta_coeff(2) * (w(ilpp^s,1:nw_recon) - w(il^s,1:nw_recon))**2
457 beta(il^s,1:nw_recon,3) = beta_coeff(1) * (w(ilp^s,1:nw_recon) + w(ilm^s,1:nw_recon) - 2.0d0 * w(il^s,1:nw_recon))**2 &
458 + beta_coeff(2) * (3.0d0 * w(ilp^s, 1:nw_recon) - 4.0d0 * w(il^s,1:nw_recon) + w(ilm^s,1:nw_recon))**2
459
460
461 gam_sum(il^s,1:nw_recon)=0.0d0
462 tau(il^s,1:nw_recon)=abs(beta(il^s,1:nw_recon,1)-beta(il^s,1:nw_recon,3))
463 do i=1,3
464 kai1(il^s,1:nw_recon,i)=(tau(il^s,1:nw_recon)/(beta(il^s,1:nw_recon,i)+weno_eps_machine))
465 gam(il^s,1:nw_recon,i)=(1.d0+kai1(il^s,1:nw_recon,i))**6
466 gam_sum(il^s,1:nw_recon)=gam_sum(il^s,1:nw_recon)+gam(il^s,1:nw_recon,i)
467 end do
468 theta(il^s,1:nw_recon)=one/(one+maxval(kai1(il^s,1:nw_recon,1:3)/10.d0,dim=ndim+2))
469 marray(il^s,1:nw_recon)=-floor(4.d0+theta(il^s,1:nw_recon)*6.d0)
470 do i=1,3
471 kai(il^s,1:nw_recon,i) = gam(il^s,1:nw_recon,i)/gam_sum(il^s,1:nw_recon)
472 where(kai(il^s,1:nw_recon,i) .lt. 10**marray(il^s,1:nw_recon))
473 delta(il^s,1:nw_recon,i)=zero
474 else where
475 delta(il^s,1:nw_recon,i)=one
476 end where
477 end do
478 delta_sum=zero
479 do i = 1,3
480 delta_sum(il^s,1:nw_recon)=delta_sum(il^s,1:nw_recon)+delta(il^s,1:nw_recon,i)*d_array(i)
481 end do
482 flux(il^s,1:nw_recon)=0.0d0
483 do i = 1,3
484 flux(il^s,1:nw_recon)=flux(il^s,1:nw_recon)+f_array(il^s,1:nw_recon,i)*delta(il^s,1:nw_recon,i)*d_array(i)/(delta_sum(il^s,1:nw_recon))
485 end do
486
487 !> right value at right interface
488 wrctmp(il^s,1:nw_recon)=flux(il^s,1:nw_recon)
489
490 call fix_limiter(ixi^l,il^l,wlctmp,wrctmp,wlc,wrc)
491
492 end subroutine teno5adlimiter
493
494 subroutine weno5nmlimiter(ixI^L,iL^L,idims,dxdim,w,wLC,wRC,var)
495 use mod_global_parameters
496
497 integer, intent(in) :: ixi^l,il^l,idims,var
498 double precision, intent(in) :: dxdim
499 double precision, intent(in) :: w(ixi^s,1:nw)
500 double precision, intent(inout) :: wrc(ixi^s,1:nw),wlc(ixi^s,1:nw)
501 !> local
502 double precision :: f_array(ixi^s,1:nw,3), d_array(3)
503 double precision :: beta(ixi^s,1:nw,3), beta_coeff(2)
504 double precision :: tau(ixi^s,1:nw), tmp(ixi^s,1:nw)
505 double precision :: u1_coeff(3), u2_coeff(3), u3_coeff(3)
506 double precision :: alpha_array(ixi^s,1:nw,3), alpha_sum(ixi^s,1:nw), flux(ixi^s,1:nw)
507 double precision :: wc(ixi^s,1:nw), wd(ixi^s,1:nw), we(ixi^s,1:nw)
508 double precision :: lambda(ixi^s)
509 double precision, dimension(ixI^S,1:nw) :: wrctmp, wlctmp
510 integer :: ilm^l, ilmm^l, ilp^l, ilpp^l, ilppp^l
511 integer :: im^l, imp^l
512 integer :: i,j
513
514 ilm^l=il^l-kr(idims,^d);
515 ilmm^l=ilm^l-kr(idims,^d);
516 ilp^l=il^l+kr(idims,^d);
517 ilpp^l=ilp^l+kr(idims,^d);
518 ilppp^l=ilpp^l+kr(idims,^d);
519 immin^d=ilmmin^d;
520 immax^d=ilpmax^d;
521 imp^l=im^l+kr(idims,^d);
522 lambda(il^s)=block%dx(il^s,idims)**weno_dx_exp
523 beta_coeff(1:2) = (/ 1.0833333333333333d0, 0.25d0/)
524 d_array(1:3) = (/ 1.0d0/10.0d0, 3.0d0/5.0d0, 3.0d0/10.0d0 /)
525 u1_coeff(1:3) = (/ -2.d0/3.d0, -1.d0/3.d0, 2.d0 /)
526 u2_coeff(1:3) = (/ -1.d0/3.d0, 2.d0/3.d0, 2.d0/3.d0 /)
527 u3_coeff(1:3) = (/ 2.d0/3.d0, 2.d0/3.d0, -1.d0/3.d0 /)
528 do i = 1, nw_recon
529 wc(im^s,i) = (block%dx(imp^s,idims) * w(im^s,i) + block%dx(im^s,idims) * w(imp^s,i)) / &
530 (block%dx(imp^s,idims) + block%dx(im^s,idims))
531 wd(il^s,i) = ((2.d0 * block%dx(ilm^s,idims) + block%dx(ilmm^s,idims)) * w(ilm^s,i) - block%dx(ilm^s,idims) * w(ilmm^s,i)) / &
532 (block%dx(ilmm^s,idims) + block%dx(ilm^s,idims))
533 we(il^s,i) = ((2.d0 * block%dx(ilpp^s,idims) + block%dx(ilppp^s,idims)) * w(ilpp^s,i) - block%dx(ilpp^s,idims) * w(ilppp^s,i)) / &
534 (block%dx(ilppp^s,idims) + block%dx(ilpp^s,idims))
535 enddo
536 !> left side
537 f_array(il^s,1:nw_recon,1) = u1_coeff(1) * wd(il^s,1:nw_recon) + u1_coeff(2) * wc(ilm^s,1:nw_recon)+ u1_coeff(3) * w(il^s,1:nw_recon)
538 f_array(il^s,1:nw_recon,2) = u2_coeff(1) * wc(ilm^s,1:nw_recon) + u2_coeff(2) * w(il^s,1:nw_recon) + u2_coeff(3) * wc(il^s,1:nw_recon)
539 f_array(il^s,1:nw_recon,3) = u3_coeff(1) * wc(il^s,1:nw_recon) + u3_coeff(2) * w(ilp^s,1:nw_recon) + u3_coeff(3) * wc(ilp^s,1:nw_recon)
540
541 beta(il^s,1:nw_recon,1) = beta_coeff(1) * (wc(ilm^s,1:nw_recon) - wd(il^s,1:nw_recon))**2 &
542 + beta_coeff(2) * (2.d0 * w(il^s,1:nw_recon) - wc(ilm^s,1:nw_recon) - wd(il^s,1:nw_recon))**2
543 beta(il^s,1:nw_recon,2) = beta_coeff(1) * (wc(ilm^s,1:nw_recon) + wc(il^s,1:nw_recon) - 2.0d0 * w(il^s,1:nw_recon))**2 &
544 + beta_coeff(2) * (wc(ilm^s,1:nw_recon) - wc(il^s,1:nw_recon))**2
545 beta(il^s,1:nw_recon,3) = beta_coeff(1) * (wc(il^s,1:nw_recon) + wc(ilp^s,1:nw_recon) - 2.0d0 * w(ilp^s,1:nw_recon))**2 &
546 + beta_coeff(2) * (3.0d0 * wc(il^s, 1:nw_recon) - 4.0d0 * w(ilp^s,1:nw_recon) + wc(ilp^s,1:nw_recon))**2
547 alpha_sum(il^s,1:nw_recon) = 0.0d0
548 select case(var)
549 ! case1 for wenojs, case2 for wenoz, case3 for wenoz+
550 case(1)
551 do i = 1,3
552 alpha_array(il^s,1:nw_recon,i) = d_array(i)/(4.d0 * beta(il^s,1:nw_recon,i) + weno_eps_machine)**2
553 alpha_sum(il^s,1:nw_recon) = alpha_sum(il^s,1:nw_recon) + alpha_array(il^s,1:nw_recon,i)
554 end do
555 case(2)
556 tau(il^s,1:nw_recon) = abs(beta(il^s,1:nw_recon,1) - beta(il^s,1:nw_recon,3)) * 4.d0
557 do i = 1,3
558 alpha_array(il^s,1:nw_recon,i) = d_array(i) * (1.d0 + (tau(il^s,1:nw_recon) / &
559 (4.d0 * beta(il^s,1:nw_recon,i) + weno_eps_machine))**2)
560 alpha_sum(il^s,1:nw_recon) = alpha_sum(il^s,1:nw_recon) + alpha_array(il^s,1:nw_recon,i)
561 end do
562 case(3)
563 tau(il^s,1:nw_recon) = abs(beta(il^s,1:nw_recon,1) - beta(il^s,1:nw_recon,3)) * 4.d0
564 do i=1,3
565 do j=1,nw_recon
566 tmp(il^s,j) = (tau(il^s,j) + weno_eps_machine) / (4.d0 * beta(il^s,j,i) + weno_eps_machine)
567 alpha_array(il^s,j,i) = d_array(i) * (1.0d0 + tmp(il^s,j)**2 + lambda(il^s)/tmp(il^s,j))
568 alpha_sum(il^s,j) = alpha_sum(il^s,j) + alpha_array(il^s,j,i)
569 end do
570 end do
571 end select
572 flux(il^s,1:nw_recon) = 0.0d0
573 do i = 1,3
574 flux(il^s,1:nw_recon) = flux(il^s,1:nw_recon) + f_array(il^s,1:nw_recon,i) * alpha_array(il^s,1:nw_recon,i)/(alpha_sum(il^s,1:nw_recon))
575 end do
576 !> left value at right interface
577 wlctmp(il^s,1:nw_recon) = flux(il^s,1:nw_recon)
578
579 !> right side
580 f_array(il^s,1:nw_recon,1) = u1_coeff(1) * we(il^s,1:nw_recon) + u1_coeff(2) * wc(ilp^s,1:nw_recon) + u1_coeff(3) * w(ilp^s,1:nw_recon)
581 f_array(il^s,1:nw_recon,2) = u2_coeff(1) * wc(ilp^s,1:nw_recon) + u2_coeff(2) * w(ilp^s,1:nw_recon) + u2_coeff(3) * wc(il^s,1:nw_recon)
582 f_array(il^s,1:nw_recon,3) = u3_coeff(1) * wc(il^s,1:nw_recon) + u3_coeff(2) * w(il^s,1:nw_recon) + u3_coeff(3) * wc(ilm^s,1:nw_recon)
583 beta(il^s,1:nw_recon,1) = beta_coeff(1) * (wc(ilp^s,1:nw_recon) - we(il^s,1:nw_recon))**2 &
584 + beta_coeff(2) * (2.d0 * w(ilp^s,1:nw_recon) - wc(ilp^s,1:nw_recon) - we(il^s,1:nw_recon))**2
585 beta(il^s,1:nw_recon,2) = beta_coeff(1) * (wc(ilp^s,1:nw_recon) + wc(il^s,1:nw_recon) - 2.0d0 * w(ilp^s,1:nw_recon))**2 &
586 + beta_coeff(2) * (wc(ilp^s,1:nw_recon) - wc(il^s,1:nw_recon))**2
587 beta(il^s,1:nw_recon,3) = beta_coeff(1) * (wc(il^s,1:nw_recon) + wc(ilm^s,1:nw_recon) - 2.0d0 * w(il^s,1:nw_recon))**2 &
588 + beta_coeff(2) * (3.0d0 * wc(il^s, 1:nw_recon) - 4.0d0 * w(il^s,1:nw_recon) + wc(ilm^s,1:nw_recon))**2
589 alpha_sum(il^s,1:nw_recon) = 0.0d0
590 select case(var)
591 case(1)
592 do i = 1,3
593 alpha_array(il^s,1:nw_recon,i) = d_array(i)/(4.d0 * beta(il^s,1:nw_recon,i) + weno_eps_machine)**2
594 alpha_sum(il^s,1:nw_recon) = alpha_sum(il^s,1:nw_recon) + alpha_array(il^s,1:nw_recon,i)
595 end do
596 case(2)
597 tau(il^s,1:nw_recon) = abs(beta(il^s,1:nw_recon,1) - beta(il^s,1:nw_recon,3)) * 4.d0
598 do i = 1,3
599 alpha_array(il^s,1:nw_recon,i) = d_array(i) * (1.d0 + (tau(il^s,1:nw_recon) / &
600 (4.d0 * beta(il^s,1:nw_recon,i) + weno_eps_machine))**2)
601 alpha_sum(il^s,1:nw_recon) = alpha_sum(il^s,1:nw_recon) + alpha_array(il^s,1:nw_recon,i)
602 end do
603 case(3)
604 tau(il^s,1:nw_recon) = abs(beta(il^s,1:nw_recon,1) - beta(il^s,1:nw_recon,3)) * 4.d0
605 do i = 1,3
606 do j = 1,nw_recon
607 tmp(il^s,j) = (tau(il^s,j) + weno_eps_machine) / (4.d0 * beta(il^s,j,i) + weno_eps_machine)
608 alpha_array(il^s,j,i) = d_array(i) * (1.0d0 + tmp(il^s,j)**2 + lambda(il^s)/tmp(il^s,j))
609 alpha_sum(il^s,j) = alpha_sum(il^s,j) + alpha_array(il^s,j,i)
610 end do
611 end do
612 end select
613 flux(il^s,1:nw_recon) = 0.0d0
614 do i = 1,3
615 flux(il^s,1:nw_recon) = flux(il^s,1:nw_recon) + f_array(il^s,1:nw_recon,i) * alpha_array(il^s,1:nw_recon,i)/(alpha_sum(il^s,1:nw_recon))
616 end do
617 !> right value at right interface
618 wrctmp(il^s,1:nw_recon) = flux(il^s,1:nw_recon)
619
620 call fix_limiter(ixi^l,il^l,wlctmp,wrctmp,wlc,wrc)
621
622 end subroutine weno5nmlimiter
623
624 subroutine weno5limiterl(ixI^L,iL^L,idims,w,wLC,var)
625 use mod_global_parameters
626
627 integer, intent(in) :: ixi^l, il^l, idims
628 integer, intent(in) :: var
629 double precision, intent(in) :: w(ixi^s,1:nw)
630 double precision, intent(inout) :: wlc(ixi^s,1:nw)
631 !> local
632 double precision :: f_array(ixi^s,1:nw,3), d_array(3)
633 double precision :: beta(ixi^s,1:nw,3), beta_coeff(2)
634 double precision :: tau(ixi^s,1:nw), tmp(ixi^s,1:nw)
635 double precision :: u1_coeff(3), u2_coeff(3), u3_coeff(3)
636 double precision :: alpha_array(ixi^s,1:nw,3), alpha_sum(ixi^s,1:nw), flux(ixi^s,1:nw)
637 double precision :: lambda(ixi^s)
638 double precision, dimension(ixI^S,1:nw) :: wlctmp
639 integer :: ilm^l, ilmm^l, ilp^l, ilpp^l, ilppp^l
640 integer :: i,j
641
642 ilm^l=il^l-kr(idims,^d);
643 ilmm^l=ilm^l-kr(idims,^d);
644 ilp^l=il^l+kr(idims,^d);
645 ilpp^l=ilp^l+kr(idims,^d);
646 lambda=block%dx(il^s,idims)**weno_dx_exp
647 beta_coeff(1:2) = (/ 1.0833333333333333d0, 0.25d0/)
648! reconstruction variation
649 d_array(1:3) = (/ 1.0d0/10.0d0, 3.0d0/5.0d0, 3.0d0/10.0d0 /)
650 u1_coeff(1:3) = (/ 1.d0/3.d0, -7.d0/6.d0, 11.d0/6.d0 /)
651 u2_coeff(1:3) = (/ -1.d0/6.d0, 5.d0/6.d0, 1.d0/3.d0 /)
652 u3_coeff(1:3) = (/ 1.d0/3.d0, 5.d0/6.d0, -1.d0/6.d0 /)
653! interpolation variation
654! d_array(1:3) = (/ 1.0d0/16.0d0, 10.0d0/16.0d0, 5.0d0/16.0d0 /)
655! u1_coeff(1:3) = (/ 3.d0/8.d0, -10.d0/8.d0, 15.d0/8.d0 /)
656! u2_coeff(1:3) = (/ -1.d0/8.d0, 6.d0/8.d0, 3.d0/8.d0 /)
657! u3_coeff(1:3) = (/ 3.d0/8.d0, 6.d0/8.d0, -1.d0/8.d0 /)
658
659 !> left side
660 f_array(il^s,1:nw_recon,1) = u1_coeff(1) * w(ilmm^s,1:nw_recon) + u1_coeff(2) * w(ilm^s,1:nw_recon) + u1_coeff(3) * w(il^s,1:nw_recon)
661 f_array(il^s,1:nw_recon,2) = u2_coeff(1) * w(ilm^s,1:nw_recon) + u2_coeff(2) * w(il^s,1:nw_recon) + u2_coeff(3) * w(ilp^s,1:nw_recon)
662 f_array(il^s,1:nw_recon,3) = u3_coeff(1) * w(il^s,1:nw_recon) + u3_coeff(2) * w(ilp^s,1:nw_recon) + u3_coeff(3) * w(ilpp^s,1:nw_recon)
663
664 beta(il^s,1:nw_recon,1) = beta_coeff(1) * (w(ilmm^s,1:nw_recon) + w(il^s,1:nw_recon) - 2.0d0*w(ilm^s,1:nw_recon))**2 &
665 + beta_coeff(2) * (w(ilmm^s,1:nw_recon) - 4.0d0 * w(ilm^s,1:nw_recon) + 3.0d0*w(il^s,1:nw_recon))**2
666 beta(il^s,1:nw_recon,2) = beta_coeff(1) * (w(ilm^s,1:nw_recon) + w(ilp^s,1:nw_recon) - 2.0d0 * w(il^s,1:nw_recon))**2 &
667 + beta_coeff(2) * (w(ilm^s,1:nw_recon) - w(ilp^s,1:nw_recon))**2
668 beta(il^s,1:nw_recon,3) = beta_coeff(1) * (w(il^s,1:nw_recon) + w(ilpp^s,1:nw_recon) - 2.0d0 * w(ilp^s,1:nw_recon))**2 &
669 + beta_coeff(2) * (3.0d0 * w(il^s, 1:nw_recon) - 4.0d0 * w(ilp^s,1:nw_recon) + w(ilpp^s,1:nw_recon))**2
670
671 alpha_sum(il^s,1:nw_recon) = 0.0d0
672 select case(var)
673 ! case1 for wenojs, case2 for wenoz
674 case(1)
675 do i = 1,3
676 alpha_array(il^s,1:nw_recon,i) = d_array(i)/(beta(il^s,1:nw_recon,i) + weno_eps_machine)**2
677 alpha_sum(il^s,1:nw_recon) = alpha_sum(il^s,1:nw_recon) + alpha_array(il^s,1:nw_recon,i)
678 end do
679 case(2)
680 tau(il^s,1:nw_recon) = abs(beta(il^s,1:nw_recon,1) - beta(il^s,1:nw_recon,3))
681 do i = 1,3
682 alpha_array(il^s,1:nw_recon,i) = d_array(i) * (1.d0 + (tau(il^s,1:nw_recon) / &
683 (beta(il^s,1:nw_recon,i) + weno_eps_machine))**2)
684 alpha_sum(il^s,1:nw_recon) = alpha_sum(il^s,1:nw_recon) + alpha_array(il^s,1:nw_recon,i)
685 end do
686 case(3)
687 tau(il^s,1:nw_recon) = abs(beta(il^s,1:nw_recon,1) - beta(il^s,1:nw_recon,3)) * 4.d0
688 do i=1,3
689 do j=1,nw_recon
690 tmp(il^s,j) = (tau(il^s,j) + weno_eps_machine) / (4.d0 * beta(il^s,j,i) + weno_eps_machine)
691 alpha_array(il^s,j,i) = d_array(i) * (1.0d0 + tmp(il^s,j)**2 + lambda(il^s)/tmp(il^s,j))
692 alpha_sum(il^s,j) = alpha_sum(il^s,j) + alpha_array(il^s,j,i)
693 end do
694 end do
695 end select
696 flux(il^s,1:nw_recon) = 0.0d0
697 do i = 1,3
698 flux(il^s,1:nw_recon) = flux(il^s,1:nw_recon) + f_array(il^s,1:nw_recon,i) * alpha_array(il^s,1:nw_recon,i)/(alpha_sum(il^s,1:nw_recon))
699 end do
700
701 !> left value at right interface
702 wlctmp(il^s,1:nw_recon) = flux(il^s,1:nw_recon)
703
704 call fix_onelimiter(ixi^l,il^l,wlctmp,wlc)
705
706 end subroutine weno5limiterl
707
708 subroutine weno5limiterr(ixI^L,iL^L,idims,w,wRC,var)
709 use mod_global_parameters
710
711 integer, intent(in) :: ixi^l, il^l, idims
712 integer, intent(in) :: var
713 double precision, intent(in) :: w(ixi^s,1:nw)
714 double precision, intent(inout) :: wrc(ixi^s,1:nw)
715 !> local
716 double precision :: f_array(ixi^s,1:nw,3), d_array(3)
717 double precision :: beta(ixi^s,1:nw,3), beta_coeff(2)
718 double precision :: tau(ixi^s,1:nw), tmp(ixi^s,1:nw)
719 double precision :: u1_coeff(3), u2_coeff(3), u3_coeff(3)
720 double precision :: alpha_array(ixi^s,1:nw,3), alpha_sum(ixi^s,1:nw), flux(ixi^s,1:nw)
721 double precision :: lambda(ixi^s)
722 double precision, dimension(ixI^S,1:nw) :: wrctmp
723 integer :: ilm^l, ilmm^l, ilp^l, ilpp^l, ilppp^l
724 integer :: i,j
725
726 ilm^l=il^l-kr(idims,^d);
727 ilp^l=il^l+kr(idims,^d);
728 ilpp^l=ilp^l+kr(idims,^d);
729 ilppp^l=ilpp^l+kr(idims,^d);
730 lambda=block%dx(il^s,idims)**weno_dx_exp
731 beta_coeff(1:2) = (/ 1.0833333333333333d0, 0.25d0/)
732! reconstruction variation
733 d_array(1:3) = (/ 1.0d0/10.0d0, 3.0d0/5.0d0, 3.0d0/10.0d0 /)
734 u1_coeff(1:3) = (/ 1.d0/3.d0, -7.d0/6.d0, 11.d0/6.d0 /)
735 u2_coeff(1:3) = (/ -1.d0/6.d0, 5.d0/6.d0, 1.d0/3.d0 /)
736 u3_coeff(1:3) = (/ 1.d0/3.d0, 5.d0/6.d0, -1.d0/6.d0 /)
737! interpolation variation
738! d_array(1:3) = (/ 1.0d0/16.0d0, 10.0d0/16.0d0, 5.0d0/16.0d0 /)
739! u1_coeff(1:3) = (/ 3.d0/8.d0, -10.d0/8.d0, 15.d0/8.d0 /)
740! u2_coeff(1:3) = (/ -1.d0/8.d0, 6.d0/8.d0, 3.d0/8.d0 /)
741! u3_coeff(1:3) = (/ 3.d0/8.d0, 6.d0/8.d0, -1.d0/8.d0 /)
742
743 !> right side
744 f_array(il^s,1:nw_recon,1) = u1_coeff(1) * w(ilppp^s,1:nw_recon) + u1_coeff(2) * w(ilpp^s,1:nw_recon) + u1_coeff(3) * w(ilp^s,1:nw_recon)
745 f_array(il^s,1:nw_recon,2) = u2_coeff(1) * w(ilpp^s,1:nw_recon) + u2_coeff(2) * w(ilp^s,1:nw_recon) + u2_coeff(3) * w(il^s,1:nw_recon)
746 f_array(il^s,1:nw_recon,3) = u3_coeff(1) * w(ilp^s,1:nw_recon) + u3_coeff(2) * w(il^s,1:nw_recon) + u3_coeff(3) * w(ilm^s,1:nw_recon)
747
748 beta(il^s,1:nw_recon,1) = beta_coeff(1) * (w(ilppp^s,1:nw_recon) + w(ilp^s,1:nw_recon) - 2.0d0*w(ilpp^s,1:nw_recon))**2 &
749 + beta_coeff(2) * (w(ilppp^s,1:nw_recon) - 4.0d0 * w(ilpp^s,1:nw_recon) + 3.0d0*w(ilp^s,1:nw_recon))**2
750 beta(il^s,1:nw_recon,2) = beta_coeff(1) * (w(ilpp^s,1:nw_recon) + w(il^s,1:nw_recon) - 2.0d0 * w(ilp^s,1:nw_recon))**2 &
751 + beta_coeff(2) * (w(ilpp^s,1:nw_recon) - w(il^s,1:nw_recon))**2
752 beta(il^s,1:nw_recon,3) = beta_coeff(1) * (w(ilp^s,1:nw_recon) + w(ilm^s,1:nw_recon) - 2.0d0 * w(il^s,1:nw_recon))**2 &
753 + beta_coeff(2) * (3.0d0 * w(ilp^s, 1:nw_recon) - 4.0d0 * w(il^s,1:nw_recon) + w(ilm^s,1:nw_recon))**2
754
755 alpha_sum(il^s,1:nw_recon) = 0.0d0
756 select case(var)
757 case(1)
758 do i = 1,3
759 alpha_array(il^s,1:nw_recon,i) = d_array(i)/(beta(il^s,1:nw_recon,i) + weno_eps_machine)**2
760 alpha_sum(il^s,1:nw_recon) = alpha_sum(il^s,1:nw_recon) + alpha_array(il^s,1:nw_recon,i)
761 end do
762 case(2)
763 tau(il^s,1:nw_recon) = abs(beta(il^s,1:nw_recon,1) - beta(il^s,1:nw_recon,3))
764 do i = 1,3
765 alpha_array(il^s,1:nw_recon,i) = d_array(i) * (1.d0 + (tau(il^s,1:nw_recon) / &
766 (beta(il^s,1:nw_recon,i) + weno_eps_machine))**2)
767 alpha_sum(il^s,1:nw_recon) = alpha_sum(il^s,1:nw_recon) + alpha_array(il^s,1:nw_recon,i)
768 end do
769 case(3)
770 tau(il^s,1:nw_recon) = abs(beta(il^s,1:nw_recon,1) - beta(il^s,1:nw_recon,3)) * 4.d0
771 do i = 1,3
772 do j = 1,nw_recon
773 tmp(il^s,j) = (tau(il^s,j) + weno_eps_machine) / (4.d0 * beta(il^s,j,i) + weno_eps_machine)
774 alpha_array(il^s,j,i) = d_array(i) * (1.0d0 + tmp(il^s,j)**2 + lambda(il^s)/tmp(il^s,j))
775 alpha_sum(il^s,j) = alpha_sum(il^s,j) + alpha_array(il^s,j,i)
776 end do
777 end do
778 end select
779 flux(il^s,1:nw_recon) = 0.0d0
780 do i = 1,3
781 flux(il^s,1:nw_recon) = flux(il^s,1:nw_recon) + f_array(il^s,1:nw_recon,i) * alpha_array(il^s,1:nw_recon,i)/(alpha_sum(il^s,1:nw_recon))
782 end do
783
784 !> right value at right interface
785 wrctmp(il^s,1:nw_recon) = flux(il^s,1:nw_recon)
786
787 call fix_onelimiter(ixi^l,il^l,wrctmp,wrc)
788
789 end subroutine weno5limiterr
790
791 subroutine weno5nmlimiterl(ixI^L,iL^L,idims,w,wLC,var)
792 use mod_global_parameters
793
794 integer, intent(in) :: ixi^l,il^l,idims,var
795 double precision, intent(in) :: w(ixi^s,1:nw)
796 double precision, intent(inout) :: wlc(ixi^s,1:nw)
797 !> local
798 double precision :: f_array(ixi^s,1:nw,3), d_array(3)
799 double precision :: beta(ixi^s,1:nw,3), beta_coeff(2)
800 double precision :: tau(ixi^s,1:nw), tmp(ixi^s,1:nw)
801 double precision :: u1_coeff(3), u2_coeff(3), u3_coeff(3)
802 double precision :: alpha_array(ixi^s,1:nw,3), alpha_sum(ixi^s,1:nw), flux(ixi^s,1:nw)
803 double precision :: wc(ixi^s,1:nw), wd(ixi^s,1:nw)
804 double precision :: lambda(ixi^s)
805 double precision, dimension(ixI^S,1:nw) :: wlctmp
806 integer :: ilm^l, ilmm^l, ilp^l, ilpp^l
807 integer :: im^l, imp^l
808 integer :: i,j
809
810 ilm^l=il^l-kr(idims,^d);
811 ilmm^l=ilm^l-kr(idims,^d);
812 ilp^l=il^l+kr(idims,^d);
813 ilpp^l=ilp^l+kr(idims,^d);
814 immin^d=ilmmin^d;
815 immax^d=ilpmax^d;
816 imp^l=im^l+kr(idims,^d);
817 lambda(il^s)=block%dx(il^s,idims)**weno_dx_exp
818 beta_coeff(1:2) = (/ 1.0833333333333333d0, 0.25d0/)
819 d_array(1:3) = (/ 1.0d0/10.0d0, 3.0d0/5.0d0, 3.0d0/10.0d0 /)
820 u1_coeff(1:3) = (/ -2.d0/3.d0, -1.d0/3.d0, 2.d0 /)
821 u2_coeff(1:3) = (/ -1.d0/3.d0, 2.d0/3.d0, 2.d0/3.d0 /)
822 u3_coeff(1:3) = (/ 2.d0/3.d0, 2.d0/3.d0, -1.d0/3.d0 /)
823 do i = 1, nw_recon
824 wc(im^s,i) = (block%dx(imp^s,idims) * w(im^s,i) + block%dx(im^s,idims) * w(imp^s,i)) / &
825 (block%dx(imp^s,idims) + block%dx(im^s,idims))
826 wd(il^s,i) = ((2.d0 * block%dx(ilm^s,idims) + block%dx(ilmm^s,idims)) * w(ilm^s,i) - block%dx(ilm^s,idims) * w(ilmm^s,i)) / &
827 (block%dx(ilmm^s,idims) + block%dx(ilm^s,idims))
828 enddo
829 f_array(il^s,1:nw_recon,1) = u1_coeff(1) * wd(il^s,1:nw_recon) + u1_coeff(2) * wc(ilm^s,1:nw_recon)+ u1_coeff(3) * w(il^s,1:nw_recon)
830 f_array(il^s,1:nw_recon,2) = u2_coeff(1) * wc(ilm^s,1:nw_recon) + u2_coeff(2) * w(il^s,1:nw_recon) + u2_coeff(3) * wc(il^s,1:nw_recon)
831 f_array(il^s,1:nw_recon,3) = u3_coeff(1) * wc(il^s,1:nw_recon) + u3_coeff(2) * w(ilp^s,1:nw_recon) + u3_coeff(3) * wc(ilp^s,1:nw_recon)
832 beta(il^s,1:nw_recon,1) = beta_coeff(1) * (wc(ilm^s,1:nw_recon) - wd(il^s,1:nw_recon))**2 &
833 + beta_coeff(2) * (2.d0 * w(il^s,1:nw_recon) - wc(ilm^s,1:nw_recon) - wd(il^s,1:nw_recon))**2
834 beta(il^s,1:nw_recon,2) = beta_coeff(1) * (wc(ilm^s,1:nw_recon) + wc(il^s,1:nw_recon) - 2.0d0 * w(il^s,1:nw_recon))**2 &
835 + beta_coeff(2) * (wc(ilm^s,1:nw_recon) - wc(il^s,1:nw_recon))**2
836 beta(il^s,1:nw_recon,3) = beta_coeff(1) * (wc(il^s,1:nw_recon) + wc(ilp^s,1:nw_recon) - 2.0d0 * w(ilp^s,1:nw_recon))**2 &
837 + beta_coeff(2) * (3.0d0 * wc(il^s, 1:nw_recon) - 4.0d0 * w(ilp^s,1:nw_recon) + wc(ilp^s,1:nw_recon))**2
838 alpha_sum(il^s,1:nw_recon) = 0.0d0
839 select case(var)
840 ! case1 for wenojs, case2 for wenoz, case3 for wenoz+
841 case(1)
842 do i = 1,3
843 alpha_array(il^s,1:nw_recon,i) = d_array(i)/(4.d0 * beta(il^s,1:nw_recon,i) + weno_eps_machine)**2
844 alpha_sum(il^s,1:nw_recon) = alpha_sum(il^s,1:nw_recon) + alpha_array(il^s,1:nw_recon,i)
845 end do
846 case(2)
847 tau(il^s,1:nw_recon) = abs(beta(il^s,1:nw_recon,1) - beta(il^s,1:nw_recon,3)) * 4.d0
848 do i = 1,3
849 alpha_array(il^s,1:nw_recon,i) = d_array(i) * (1.d0 + (tau(il^s,1:nw_recon) / &
850 (4.d0 * beta(il^s,1:nw_recon,i) + weno_eps_machine))**2)
851 alpha_sum(il^s,1:nw_recon) = alpha_sum(il^s,1:nw_recon) + alpha_array(il^s,1:nw_recon,i)
852 end do
853 case(3)
854 tau(il^s,1:nw_recon) = abs(beta(il^s,1:nw_recon,1) - beta(il^s,1:nw_recon,3)) * 4.d0
855 do i=1,3
856 do j=1,nw_recon
857 tmp(il^s,j) = (tau(il^s,j) + weno_eps_machine) / (4.d0 * beta(il^s,j,i) + weno_eps_machine)
858 alpha_array(il^s,j,i) = d_array(i) * (1.0d0 + tmp(il^s,j)**2 + lambda(il^s)/tmp(il^s,j))
859 alpha_sum(il^s,j) = alpha_sum(il^s,j) + alpha_array(il^s,j,i)
860 end do
861 end do
862 end select
863 flux(il^s,1:nw_recon) = 0.0d0
864 do i = 1,3
865 flux(il^s,1:nw_recon) = flux(il^s,1:nw_recon) + f_array(il^s,1:nw_recon,i) * alpha_array(il^s,1:nw_recon,i)/(alpha_sum(il^s,1:nw_recon))
866 end do
867 !> left value at right interface
868 wlctmp(il^s,1:nw_recon) = flux(il^s,1:nw_recon)
869
870 call fix_onelimiter(ixi^l,il^l,wlctmp,wlc)
871
872 end subroutine weno5nmlimiterl
873
874 subroutine weno5nmlimiterr(ixI^L,iL^L,idims,w,wRC,var)
875 use mod_global_parameters
876
877 integer, intent(in) :: ixi^l,il^l,idims,var
878 double precision, intent(in) :: w(ixi^s,1:nw)
879 double precision, intent(inout) :: wrc(ixi^s,1:nw)
880 !> local
881 double precision :: f_array(ixi^s,1:nw,3), d_array(3)
882 double precision :: beta(ixi^s,1:nw,3), beta_coeff(2)
883 double precision :: tau(ixi^s,1:nw), tmp(ixi^s,1:nw)
884 double precision :: u1_coeff(3), u2_coeff(3), u3_coeff(3)
885 double precision :: alpha_array(ixi^s,1:nw,3), alpha_sum(ixi^s,1:nw), flux(ixi^s,1:nw)
886 double precision :: wc(ixi^s,1:nw), we(ixi^s,1:nw)
887 double precision :: lambda(ixi^s)
888 double precision, dimension(ixI^S,1:nw) :: wrctmp
889 integer :: ilm^l,ilp^l,ilpp^l,ilppp^l
890 integer :: im^l, imp^l
891 integer :: i,j
892
893 ilm^l=il^l-kr(idims,^d);
894 ilp^l=il^l+kr(idims,^d);
895 ilpp^l=ilp^l+kr(idims,^d);
896 ilppp^l=ilpp^l+kr(idims,^d);
897 immin^d=ilmmin^d;
898 immax^d=ilpmax^d;
899 imp^l=im^l+kr(idims,^d);
900 lambda(il^s)=block%dx(il^s,idims)**weno_dx_exp
901 beta_coeff(1:2) = (/ 1.0833333333333333d0, 0.25d0/)
902 d_array(1:3) = (/ 1.0d0/10.0d0, 3.0d0/5.0d0, 3.0d0/10.0d0 /)
903 u1_coeff(1:3) = (/ -2.d0/3.d0, -1.d0/3.d0, 2.d0 /)
904 u2_coeff(1:3) = (/ -1.d0/3.d0, 2.d0/3.d0, 2.d0/3.d0 /)
905 u3_coeff(1:3) = (/ 2.d0/3.d0, 2.d0/3.d0, -1.d0/3.d0 /)
906 do i = 1, nw_recon
907 wc(im^s,i) = (block%dx(imp^s,idims) * w(im^s,i) + block%dx(im^s,idims) * w(imp^s,i)) / &
908 (block%dx(imp^s,idims) + block%dx(im^s,idims))
909 we(il^s,i) = ((2.d0 * block%dx(ilpp^s,idims) + block%dx(ilppp^s,idims)) * w(ilpp^s,i) - block%dx(ilpp^s,idims) * w(ilppp^s,i)) / &
910 (block%dx(ilppp^s,idims) + block%dx(ilpp^s,idims))
911 enddo
912 !> right side
913 f_array(il^s,1:nw_recon,1) = u1_coeff(1) * we(il^s,1:nw_recon) + u1_coeff(2) * wc(ilp^s,1:nw_recon) + u1_coeff(3) * w(ilp^s,1:nw_recon)
914 f_array(il^s,1:nw_recon,2) = u2_coeff(1) * wc(ilp^s,1:nw_recon) + u2_coeff(2) * w(ilp^s,1:nw_recon) + u2_coeff(3) * wc(il^s,1:nw_recon)
915 f_array(il^s,1:nw_recon,3) = u3_coeff(1) * wc(il^s,1:nw_recon) + u3_coeff(2) * w(il^s,1:nw_recon) + u3_coeff(3) * wc(ilm^s,1:nw_recon)
916 beta(il^s,1:nw_recon,1) = beta_coeff(1) * (wc(ilp^s,1:nw_recon) - we(il^s,1:nw_recon))**2 &
917 + beta_coeff(2) * (2.d0 * w(ilp^s,1:nw_recon) - wc(ilp^s,1:nw_recon) - we(il^s,1:nw_recon))**2
918 beta(il^s,1:nw_recon,2) = beta_coeff(1) * (wc(ilp^s,1:nw_recon) + wc(il^s,1:nw_recon) - 2.0d0 * w(ilp^s,1:nw_recon))**2 &
919 + beta_coeff(2) * (wc(ilp^s,1:nw_recon) - wc(il^s,1:nw_recon))**2
920 beta(il^s,1:nw_recon,3) = beta_coeff(1) * (wc(il^s,1:nw_recon) + wc(ilm^s,1:nw_recon) - 2.0d0 * w(il^s,1:nw_recon))**2 &
921 + beta_coeff(2) * (3.0d0 * wc(il^s, 1:nw_recon) - 4.0d0 * w(il^s,1:nw_recon) + wc(ilm^s,1:nw_recon))**2
922 alpha_sum(il^s,1:nw_recon) = 0.0d0
923 select case(var)
924 case(1)
925 do i = 1,3
926 alpha_array(il^s,1:nw_recon,i) = d_array(i)/(4.d0 * beta(il^s,1:nw_recon,i) + weno_eps_machine)**2
927 alpha_sum(il^s,1:nw_recon) = alpha_sum(il^s,1:nw_recon) + alpha_array(il^s,1:nw_recon,i)
928 end do
929 case(2)
930 tau(il^s,1:nw_recon) = abs(beta(il^s,1:nw_recon,1) - beta(il^s,1:nw_recon,3)) * 4.d0
931 do i = 1,3
932 alpha_array(il^s,1:nw_recon,i) = d_array(i) * (1.d0 + (tau(il^s,1:nw_recon) / &
933 (4.d0 * beta(il^s,1:nw_recon,i) + weno_eps_machine))**2)
934 alpha_sum(il^s,1:nw_recon) = alpha_sum(il^s,1:nw_recon) + alpha_array(il^s,1:nw_recon,i)
935 end do
936 case(3)
937 tau(il^s,1:nw_recon) = abs(beta(il^s,1:nw_recon,1) - beta(il^s,1:nw_recon,3)) * 4.d0
938 do i = 1,3
939 do j = 1,nw_recon
940 tmp(il^s,j) = (tau(il^s,j) + weno_eps_machine) / (4.d0 * beta(il^s,j,i) + weno_eps_machine)
941 alpha_array(il^s,j,i) = d_array(i) * (1.0d0 + tmp(il^s,j)**2 + lambda(il^s)/tmp(il^s,j))
942 alpha_sum(il^s,j) = alpha_sum(il^s,j) + alpha_array(il^s,j,i)
943 end do
944 end do
945 end select
946 flux(il^s,1:nw_recon) = 0.0d0
947 do i = 1,3
948 flux(il^s,1:nw_recon) = flux(il^s,1:nw_recon) + f_array(il^s,1:nw_recon,i) * alpha_array(il^s,1:nw_recon,i)/(alpha_sum(il^s,1:nw_recon))
949 end do
950 !> right value at right interface
951 wrctmp(il^s,1:nw_recon) = flux(il^s,1:nw_recon)
952
953 call fix_onelimiter(ixi^l,il^l,wrctmp,wrc)
954
955 end subroutine weno5nmlimiterr
956
957 subroutine weno5cu6limiter(ixI^L,iL^L,idims,w,wLC,wRC)
958 use mod_global_parameters
959
960 integer, intent(in) :: ixi^l, il^l, idims
961 double precision, intent(in) :: w(ixi^s,1:nw)
962 double precision, intent(inout) :: wrc(ixi^s,1:nw),wlc(ixi^s,1:nw)
963 !> local
964 double precision :: f_array(ixi^s,1:nw,3), d_array(3)
965 double precision :: beta(ixi^s,1:nw,3), beta_coeff(2)
966 double precision :: u1_coeff(3), u2_coeff(3), u3_coeff(3)
967 double precision :: alpha_array(ixi^s,1:nw,3), alpha_sum(ixi^s,1:nw)
968 double precision :: theta2(ixi^s,1:nw)
969 double precision, parameter :: theta_limit=0.7d0
970 double precision, dimension(ixI^S,1:nw) :: wrctmp, wlctmp
971 integer :: ilm^l, ilmm^l, ilp^l, ilpp^l, ilppp^l
972 integer :: i
973
974 ilm^l=il^l-kr(idims,^d);
975 ilmm^l=ilm^l-kr(idims,^d);
976 ilp^l=il^l+kr(idims,^d);
977 ilpp^l=ilp^l+kr(idims,^d);
978 ilppp^l=ilpp^l+kr(idims,^d);
979
980 beta_coeff(1:2) = (/ 1.0833333333333333d0, 0.25d0/)
981 d_array(1:3) = (/ 1.0d0/10.0d0, 3.0d0/5.0d0, 3.0d0/10.0d0 /)
982 u1_coeff(1:3) = (/ 1.d0/3.d0, -7.d0/6.d0, 11.d0/6.d0 /)
983 u2_coeff(1:3) = (/ -1.d0/6.d0, 5.d0/6.d0, 1.d0/3.d0 /)
984 u3_coeff(1:3) = (/ 1.d0/3.d0, 5.d0/6.d0, -1.d0/6.d0 /)
985
986 !> left side
987 beta(il^s,1:nw_recon,1)=beta_coeff(1)*(w(ilmm^s,1:nw_recon)+w(il^s,1:nw_recon)-2.0d0*w(ilm^s,1:nw_recon))**2&
988 +beta_coeff(2)*(w(ilmm^s,1:nw_recon)-4.0d0*w(ilm^s,1:nw_recon)+3.0d0*w(il^s,1:nw_recon))**2
989 beta(il^s,1:nw_recon,2)=beta_coeff(1)*(w(ilm^s,1:nw_recon)+w(ilp^s,1:nw_recon)-2.0d0*w(il^s,1:nw_recon))**2&
990 +beta_coeff(2)*(w(ilm^s,1:nw_recon)-w(ilp^s,1:nw_recon))**2
991 beta(il^s,1:nw_recon,3)=beta_coeff(1)*(w(il^s,1:nw_recon)+w(ilpp^s,1:nw_recon)-2.0d0*w(ilp^s,1:nw_recon))**2&
992 +beta_coeff(2)*(3.0d0*w(il^s,1:nw_recon)-4.0d0*w(ilp^s,1:nw_recon)+w(ilpp^s,1:nw_recon))**2
993 alpha_sum(il^s,1:nw_recon)=zero
994 do i=1,3
995 alpha_array(il^s,1:nw_recon,i)=d_array(i)/(beta(il^s,1:nw_recon,i)+weno_eps_machine)**2
996 alpha_sum(il^s,1:nw_recon)=alpha_sum(il^s,1:nw_recon)+alpha_array(il^s,1:nw_recon,i)
997 end do
998 do i=1,3
999 alpha_array(il^s,1:nw_recon,i)=alpha_array(il^s,1:nw_recon,i)/alpha_sum(il^s,1:nw_recon)
1000 end do
1001 theta2(il^s,1:nw_recon)=((alpha_array(il^s,1:nw_recon,1)/d_array(1)-1.d0)**2&
1002 +(alpha_array(il^s,1:nw_recon,2)/d_array(2)-1.d0)**2&
1003 +(alpha_array(il^s,1:nw_recon,3)/d_array(3)-1.d0)**2)/83.d0
1004 where(theta2(il^s,1:nw_recon) .gt. theta_limit)
1005 f_array(il^s,1:nw_recon,1)=u1_coeff(1)*w(ilmm^s,1:nw_recon)+u1_coeff(2)*w(ilm^s,1:nw_recon)+u1_coeff(3)*w(il^s,1:nw_recon)
1006 f_array(il^s,1:nw_recon,2)=u2_coeff(1)*w(ilm^s,1:nw_recon)+u2_coeff(2)*w(il^s,1:nw_recon)+u2_coeff(3)*w(ilp^s,1:nw_recon)
1007 f_array(il^s,1:nw_recon,3)=u3_coeff(1)*w(il^s,1:nw_recon)+u3_coeff(2)*w(ilp^s,1:nw_recon)+u3_coeff(3)*w(ilpp^s,1:nw_recon)
1008 wlctmp(il^s,1:nw_recon)=f_array(il^s,1:nw_recon,1)*alpha_array(il^s,1:nw_recon,1)&
1009 +f_array(il^s,1:nw_recon,2)*alpha_array(il^s,1:nw_recon,2)&
1010 +f_array(il^s,1:nw_recon,3)*alpha_array(il^s,1:nw_recon,3)
1011 else where
1012 wlctmp(il^s,1:nw_recon)=1.d0/60.d0*(w(ilmm^s,1:nw_recon)-8.d0*w(ilm^s,1:nw_recon)+37.d0*w(il^s,1:nw_recon)&
1013 +37.d0*w(ilp^s,1:nw_recon)-8.d0*w(ilpp^s,1:nw_recon)+w(ilppp^s,1:nw_recon))
1014 end where
1015
1016 !> right side
1017 beta(il^s,1:nw_recon,1)=beta_coeff(1)*(w(ilppp^s,1:nw_recon)+w(ilp^s,1:nw_recon)-2.0d0*w(ilpp^s,1:nw_recon))**2&
1018 +beta_coeff(2)*(w(ilppp^s,1:nw_recon)-4.0d0*w(ilpp^s,1:nw_recon)+3.0d0*w(ilp^s,1:nw_recon))**2
1019 beta(il^s,1:nw_recon,2)=beta_coeff(1)*(w(ilpp^s,1:nw_recon)+w(il^s,1:nw_recon)-2.0d0*w(ilp^s,1:nw_recon))**2&
1020 +beta_coeff(2)*(w(ilpp^s,1:nw_recon)-w(il^s,1:nw_recon))**2
1021 beta(il^s,1:nw_recon,3)=beta_coeff(1)*(w(ilp^s,1:nw_recon)+w(ilm^s,1:nw_recon)-2.0d0*w(il^s,1:nw_recon))**2&
1022 +beta_coeff(2)*(3.0d0*w(ilp^s,1:nw_recon)-4.0d0*w(il^s,1:nw_recon)+w(ilm^s,1:nw_recon))**2
1023 alpha_sum(il^s,1:nw_recon)=zero
1024 do i=1,3
1025 alpha_array(il^s,1:nw_recon,i)=d_array(i)/(beta(il^s,1:nw_recon,i)+weno_eps_machine)**2
1026 alpha_sum(il^s,1:nw_recon)=alpha_sum(il^s,1:nw_recon)+alpha_array(il^s,1:nw_recon,i)
1027 end do
1028 do i=1,3
1029 alpha_array(il^s,1:nw_recon,i)=alpha_array(il^s,1:nw_recon,i)/alpha_sum(il^s,1:nw_recon)
1030 end do
1031 theta2(il^s,1:nw_recon)=((alpha_array(il^s,1:nw_recon,1)/d_array(1)-1.d0)**2&
1032 +(alpha_array(il^s,1:nw_recon,2)/d_array(2)-1.d0)**2&
1033 +(alpha_array(il^s,1:nw_recon,3)/d_array(3)-1.d0)**2)/83.d0
1034 where(theta2(il^s,1:nw_recon) .gt. theta_limit)
1035 f_array(il^s,1:nw_recon,1)=u1_coeff(1)*w(ilppp^s,1:nw_recon)+u1_coeff(2)*w(ilpp^s,1:nw_recon)+u1_coeff(3)*w(ilp^s,1:nw_recon)
1036 f_array(il^s,1:nw_recon,2)=u2_coeff(1)*w(ilpp^s,1:nw_recon)+u2_coeff(2)*w(ilp^s,1:nw_recon)+u2_coeff(3)*w(il^s,1:nw_recon)
1037 f_array(il^s,1:nw_recon,3)=u3_coeff(1)*w(ilp^s,1:nw_recon)+u3_coeff(2)*w(il^s,1:nw_recon)+u3_coeff(3)*w(ilm^s,1:nw_recon)
1038 wrctmp(il^s,1:nw_recon)=f_array(il^s,1:nw_recon,1)*alpha_array(il^s,1:nw_recon,1)&
1039 +f_array(il^s,1:nw_recon,2)*alpha_array(il^s,1:nw_recon,2)&
1040 +f_array(il^s,1:nw_recon,3)*alpha_array(il^s,1:nw_recon,3)
1041 else where
1042 wrctmp(il^s,1:nw_recon)=1.d0/60.d0*(w(ilppp^s,1:nw_recon)-8.d0*w(ilpp^s,1:nw_recon)+37.d0*w(ilp^s,1:nw_recon)&
1043 +37.d0*w(il^s,1:nw_recon)-8.d0*w(ilm^s,1:nw_recon)+w(ilmm^s,1:nw_recon))
1044 end where
1045
1046 call fix_limiter(ixi^l,il^l,wlctmp,wrctmp,wlc,wrc)
1047
1048 end subroutine weno5cu6limiter
1049
1050 subroutine weno7limiter(ixI^L,iL^L,idims,w,wLC,wRC,var)
1051 use mod_global_parameters
1052
1053 integer, intent(in) :: ixi^l, il^l, idims, var
1054 double precision, intent(in) :: w(ixi^s,1:nw)
1055 double precision, intent(inout) :: wrc(ixi^s,1:nw),wlc(ixi^s,1:nw)
1056 !> local
1057 double precision, dimension(4) :: d_array, u1_coeff, u2_coeff, u3_coeff, u4_coeff
1058 double precision, dimension(ixI^S,1:nw,4) :: f_array, beta, alpha_array
1059 double precision, dimension(ixI^S) :: a, b, c, tmp, tmp2, tmp3
1060 double precision, dimension(ixI^S,1:nw) :: alpha_sum, d, dm4
1061 double precision, dimension(ixI^S,1:nw) :: flux, flux_min, flux_max, flux_ul, flux_md, flux_lc
1062 double precision, parameter :: mpalpha = 2.d0, mpbeta = 4.d0
1063 double precision, dimension(ixI^S,1:nw) :: wrctmp, wlctmp
1064 integer :: ilm^l, ilmm^l, ilmmm^l
1065 integer :: ilp^l, ilpp^l, ilppp^l, ilpppp^l
1066 integer :: id^l, idp^l, idpp^l, idm^l, ie^l, iem^l, iep^l, iepp^l
1067 integer :: i,iw
1068
1069 ilm^l=il^l-kr(idims,^d);
1070 ilmm^l=ilm^l-kr(idims,^d);
1071 ilmmm^l=ilmm^l-kr(idims,^d);
1072 ilp^l=il^l+kr(idims,^d);
1073 ilpp^l=ilp^l+kr(idims,^d);
1074 ilppp^l=ilpp^l+kr(idims,^d);
1075 ilpppp^l=ilppp^l+kr(idims,^d);
1076
1077 d_array(1:4) = (/ 1.d0/35.d0, 12.d0/35.d0, 18.d0/35.d0, 4.d0/35.d0 /)
1078 u1_coeff(1:4) = (/ -1.d0/4.d0, 13.d0/12.d0, -23.d0/12.d0, 25.d0/12.d0 /)
1079 u2_coeff(1:4) = (/ 1.d0/12.d0, -5.d0/12.d0, 13.d0/12.d0, 1.d0/4.d0 /)
1080 u3_coeff(1:4) = (/ -1.d0/12.d0, 7.d0/12.d0, 7.d0/12.d0, -1.d0/12.d0 /)
1081 u4_coeff(1:4) = (/ 1.d0/4.d0, 13.d0/12.d0, -5.d0/12.d0, 1.d0/12.d0 /)
1082
1083 !> left side
1084 f_array(il^s,1:nw_recon,1) = u1_coeff(1) * w(ilmmm^s,1:nw_recon) + u1_coeff(2) * w(ilmm^s,1:nw_recon) + u1_coeff(3) * w(ilm^s,1:nw_recon) &
1085 + u1_coeff(4) * w(il^s,1:nw_recon)
1086 f_array(il^s,1:nw_recon,2) = u2_coeff(1) * w(ilmm^s,1:nw_recon) + u2_coeff(2) * w(ilm^s,1:nw_recon) + u2_coeff(3) * w(il^s,1:nw_recon) &
1087 + u2_coeff(4) * w(ilp^s,1:nw_recon)
1088 f_array(il^s,1:nw_recon,3) = u3_coeff(1) * w(ilm^s,1:nw_recon) + u3_coeff(2) * w(il^s,1:nw_recon) + u3_coeff(3) * w(ilp^s,1:nw_recon) &
1089 + u3_coeff(4) * w(ilpp^s,1:nw_recon)
1090 f_array(il^s,1:nw_recon,4) = u4_coeff(1) * w(il^s,1:nw_recon) + u4_coeff(2) * w(ilp^s,1:nw_recon) + u4_coeff(3) * w(ilpp^s,1:nw_recon) &
1091 + u4_coeff(4) * w(ilppp^s,1:nw_recon)
1092
1093 beta(il^s,1:nw_recon,1) = w(ilmmm^s,1:nw_recon) * (547.d0 * w(ilmmm^s,1:nw_recon) - 3882.d0 * w(ilmm^s,1:nw_recon) + 4642.d0 * w(ilm^s,1:nw_recon) &
1094 - 1854.d0 * w(il^s,1:nw_recon)) &
1095 + w(ilmm^s,1:nw_recon) * (7043.d0 * w(ilmm^s,1:nw_recon) - 17246.d0 * w(ilm^s,1:nw_recon) + 7042.d0 * w(il^s,1:nw_recon)) &
1096 + w(ilm^s,1:nw_recon) * (11003.d0 * w(ilm^s,1:nw_recon) - 9402.d0 * w(il^s,1:nw_recon)) + 2107.d0 * w(il^s,1:nw_recon)**2
1097 beta(il^s,1:nw_recon,2) = w(ilmm^s,1:nw_recon) * (267.d0 * w(ilmm^s,1:nw_recon) - 1642.d0 * w(ilm^s,1:nw_recon) + 1602.d0 * w(il^s,1:nw_recon) &
1098 - 494.d0 * w(ilp^s,1:nw_recon)) &
1099 + w(ilm^s,1:nw_recon) * (2843.d0 * w(ilm^s,1:nw_recon) - 5966.d0 * w(il^s,1:nw_recon) + 1922.d0 * w(ilp^s,1:nw_recon)) &
1100 + w(il^s,1:nw_recon) * (3443.d0 * w(il^s,1:nw_recon) - 2522.d0 * w(ilp^s,1:nw_recon)) + 547.d0 * w(ilp^s,1:nw_recon) ** 2
1101 beta(il^s,1:nw_recon,3) = w(ilm^s,1:nw_recon) * (547.d0 * w(ilm^s,1:nw_recon) - 2522.d0 * w(il^s,1:nw_recon) + 1922.d0 * w(ilp^s,1:nw_recon) &
1102 - 494.d0 * w(ilpp^s,1:nw_recon)) &
1103 + w(il^s,1:nw_recon) * (3443.d0 * w(il^s,1:nw_recon) - 5966.d0 * w(ilp^s,1:nw_recon) + 1602.d0 * w(ilpp^s,1:nw_recon)) &
1104 + w(ilp^s,1:nw_recon) * (2843.d0 * w(ilp^s,1:nw_recon) - 1642.d0 * w(ilpp^s,1:nw_recon)) + 267.d0 * w(ilpp^s,1:nw_recon) ** 2
1105 beta(il^s,1:nw_recon,4) = w(il^s,1:nw_recon) * (2107.d0 * w(il^s,1:nw_recon) - 9402.d0 * w(ilp^s,1:nw_recon) + 7042.d0 * w(ilpp^s,1:nw_recon) &
1106 - 1854.d0 * w(ilppp^s,1:nw_recon)) &
1107 + w(ilp^s,1:nw_recon) * (11003.d0 * w(ilp^s,1:nw_recon) - 17246.d0 * w(ilpp^s,1:nw_recon) + 4642.d0 * w(ilppp^s,1:nw_recon)) &
1108 + w(ilpp^s,1:nw_recon) * (7043.d0 * w(ilpp^s,1:nw_recon) - 3882.d0 * w(ilppp^s,1:nw_recon)) &
1109 + 547.d0 * w(ilppp^s,1:nw_recon) ** 2
1110
1111 alpha_sum(il^s,1:nw_recon) = 0.0d0
1112 do i = 1,4
1113 alpha_array(il^s,1:nw_recon,i) = d_array(i)/(beta(il^s,1:nw_recon,i) + weno_eps_machine)**2
1114 alpha_sum(il^s,1:nw_recon) = alpha_sum(il^s,1:nw_recon) + alpha_array(il^s,1:nw_recon,i)
1115 end do
1116 flux(il^s,1:nw_recon) = 0.0d0
1117 do i = 1,4
1118 flux(il^s,1:nw_recon) = flux(il^s,1:nw_recon) + f_array(il^s,1:nw_recon,i) * alpha_array(il^s,1:nw_recon,i)/(alpha_sum(il^s,1:nw_recon))
1119 end do
1120
1121 select case(var)
1122 ! case1 for wenojs, case2 for mpweno
1123 case(1)
1124 wlctmp(il^s,1:nw_recon) = flux(il^s,1:nw_recon)
1125 case(2)
1126 idmax^d=ilmax^d; idmin^d=ilmin^d-kr(idims,^d);
1127 idm^l=id^l-kr(idims,^d);
1128 idp^l=id^l+kr(idims,^d);
1129
1130 iemax^d=idmax^d+kr(idims,^d); iemin^d=idmin^d;
1131 iem^l=ie^l-kr(idims,^d);
1132 iep^l=ie^l+kr(idims,^d);
1133
1134 d(ie^s,1:nw_recon) = w(iep^s,1:nw_recon)-2.0d0*w(ie^s,1:nw_recon)+w(iem^s,1:nw_recon)
1135
1136 do iw=1,nw_recon
1137 a(id^s) = 4.0d0*d(id^s,iw)-d(idp^s,iw)
1138 b(id^s) = 4.0d0*d(idp^s,iw)-d(id^s,iw)
1139 call minmod(ixi^l,id^l,a,b,tmp)
1140 a(id^s) = d(id^s,iw)
1141 b(id^s) = d(idp^s,iw)
1142 call minmod(ixi^l,id^l,a,b,tmp2)
1143 call minmod(ixi^l,id^l,tmp,tmp2,tmp3)
1144 dm4(id^s,iw) = tmp3(id^s)
1145 end do
1146
1147 flux_ul(il^s,1:nw_recon) = w(il^s,1:nw_recon) + mpalpha * (w(il^s,1:nw_recon) - w(ilm^s,1:nw_recon))
1148 flux_md(il^s,1:nw_recon) = half * (w(il^s,1:nw_recon) + w(ilp^s,1:nw_recon) - dm4(il^s,1:nw_recon))
1149 flux_lc(il^s,1:nw_recon) = half * (3.d0 * w(il^s,1:nw_recon) - w(ilm^s,1:nw_recon)) + mpbeta / 3.d0 * dm4(ilm^s,1:nw_recon)
1150
1151 flux_min(il^s,1:nw_recon) = max(min(w(il^s,1:nw_recon), w(ilp^s,1:nw_recon), flux_md(il^s,1:nw_recon)), &
1152 min(w(il^s,1:nw_recon), flux_ul(il^s,1:nw_recon),flux_lc(il^s,1:nw_recon)))
1153
1154 flux_max(il^s,1:nw_recon) = min(max(w(il^s,1:nw_recon), w(ilp^s,1:nw_recon), flux_md(il^s,1:nw_recon)), &
1155 max(w(il^s,1:nw_recon), flux_ul(il^s,1:nw_recon),flux_lc(il^s,1:nw_recon)))
1156 do iw=1,nw_recon
1157 a(il^s) = flux(il^s,iw)
1158 b(il^s) = flux_min(il^s,iw)
1159 c(il^s) = flux_max(il^s,iw)
1160 call median(ixi^l, il^l, a, b, c, tmp)
1161 wlctmp(il^s,iw) = tmp(il^s)
1162 end do
1163 end select
1164
1165 !> right side
1166 !>> mmm -> pppp
1167 !>> mm -> ppp
1168 !>> m -> pp
1169 !>> 0 -> p
1170 !>> p -> 0
1171 !>> pp -> m
1172 !>> ppp -> mm
1173 f_array(il^s,1:nw_recon,1) = u1_coeff(1) * w(ilpppp^s,1:nw_recon) + u1_coeff(2) * w(ilppp^s,1:nw_recon) + u1_coeff(3) * w(ilpp^s,1:nw_recon) &
1174 + u1_coeff(4) * w(ilp^s,1:nw_recon)
1175 f_array(il^s,1:nw_recon,2) = u2_coeff(1) * w(ilppp^s,1:nw_recon) + u2_coeff(2) * w(ilpp^s,1:nw_recon) + u2_coeff(3) * w(ilp^s,1:nw_recon) &
1176 + u2_coeff(4) * w(il^s,1:nw_recon)
1177 f_array(il^s,1:nw_recon,3) = u3_coeff(1) * w(ilpp^s,1:nw_recon) + u3_coeff(2) * w(ilp^s,1:nw_recon) + u3_coeff(3) * w(il^s,1:nw_recon) &
1178 + u3_coeff(4) * w(ilm^s,1:nw_recon)
1179 f_array(il^s,1:nw_recon,4) = u4_coeff(1) * w(ilp^s,1:nw_recon) + u4_coeff(2) * w(il^s,1:nw_recon) + u4_coeff(3) * w(ilm^s,1:nw_recon) &
1180 + u4_coeff(4) * w(ilmm^s,1:nw_recon)
1181
1182 beta(il^s,1:nw_recon,1) = w(ilpppp^s,1:nw_recon) * (547.d0 * w(ilpppp^s,1:nw_recon) - 3882.d0 * w(ilppp^s,1:nw_recon) + 4642.d0 * w(ilpp^s,1:nw_recon) &
1183 - 1854.d0 * w(ilp^s,1:nw_recon)) &
1184 + w(ilppp^s,1:nw_recon) * (7043.d0 * w(ilppp^s,1:nw_recon) - 17246.d0 * w(ilpp^s,1:nw_recon) + 7042.d0 * w(ilp^s,1:nw_recon)) &
1185 + w(ilpp^s,1:nw_recon) * (11003.d0 * w(ilpp^s,1:nw_recon) - 9402.d0 * w(ilp^s,1:nw_recon)) + 2107.d0 * w(ilp^s,1:nw_recon)**2
1186 beta(il^s,1:nw_recon,2) = w(ilppp^s,1:nw_recon) * (267.d0 * w(ilppp^s,1:nw_recon) - 1642.d0 * w(ilpp^s,1:nw_recon) + 1602.d0 * w(ilp^s,1:nw_recon) &
1187 - 494.d0 * w(il^s,1:nw_recon)) &
1188 + w(ilpp^s,1:nw_recon) * (2843.d0 * w(ilpp^s,1:nw_recon) - 5966.d0 * w(ilp^s,1:nw_recon) + 1922.d0 * w(il^s,1:nw_recon)) &
1189 + w(ilp^s,1:nw_recon) * (3443.d0 * w(ilp^s,1:nw_recon) - 2522.d0 * w(il^s,1:nw_recon)) + 547.d0 * w(il^s,1:nw_recon) ** 2
1190 beta(il^s,1:nw_recon,3) = w(ilpp^s,1:nw_recon) * (547.d0 * w(ilpp^s,1:nw_recon) - 2522.d0 * w(ilp^s,1:nw_recon) + 1922.d0 * w(il^s,1:nw_recon) &
1191 - 494.d0 * w(ilm^s,1:nw_recon)) &
1192 + w(ilp^s,1:nw_recon) * (3443.d0 * w(ilp^s,1:nw_recon) - 5966.d0 * w(il^s,1:nw_recon) + 1602.d0 * w(ilm^s,1:nw_recon)) &
1193 + w(il^s,1:nw_recon) * (2843.d0 * w(il^s,1:nw_recon) - 1642.d0 * w(ilm^s,1:nw_recon)) + 267.d0 * w(ilm^s,1:nw_recon) ** 2
1194 beta(il^s,1:nw_recon,4) = w(ilp^s,1:nw_recon) * (2107.d0 * w(ilp^s,1:nw_recon) - 9402.d0 * w(il^s,1:nw_recon) + 7042.d0 * w(ilm^s,1:nw_recon) &
1195 - 1854.d0 * w(ilmm^s,1:nw_recon)) &
1196 + w(il^s,1:nw_recon) * (11003.d0 * w(il^s,1:nw_recon) - 17246.d0 * w(ilm^s,1:nw_recon) + 4642.d0 * w(ilmm^s,1:nw_recon)) &
1197 + w(ilm^s,1:nw_recon) * (7043.d0 * w(ilm^s,1:nw_recon) - 3882.d0 * w(ilmm^s,1:nw_recon)) + 547.d0 * w(ilmm^s,1:nw_recon) ** 2
1198
1199 alpha_sum(il^s,1:nw_recon) = 0.0d0
1200 do i = 1,4
1201 alpha_array(il^s,1:nw_recon,i) = d_array(i)/(beta(il^s,1:nw_recon,i) + weno_eps_machine)**2
1202 alpha_sum(il^s,1:nw_recon) = alpha_sum(il^s,1:nw_recon) + alpha_array(il^s,1:nw_recon,i)
1203 end do
1204 flux(il^s,1:nw_recon) = 0.0d0
1205 do i = 1,4
1206 flux(il^s,1:nw_recon) = flux(il^s,1:nw_recon) + f_array(il^s,1:nw_recon,i) * alpha_array(il^s,1:nw_recon,i)/(alpha_sum(il^s,1:nw_recon))
1207 end do
1208
1209 select case(var)
1210 case(1)
1211 wrctmp(il^s,1:nw_recon) = flux(il^s,1:nw_recon)
1212 case(2)
1213 idmax^d=ilmax^d+kr(idims,^d); idmin^d=ilmin^d;
1214 idm^l=id^l-kr(idims,^d);
1215 idp^l=id^l+kr(idims,^d);
1216
1217 iemax^d=idmax^d; iemin^d=idmin^d-kr(idims,^d);
1218 iem^l=ie^l-kr(idims,^d);
1219 iep^l=ie^l+kr(idims,^d);
1220 iepp^l=iep^l+kr(idims,^d);
1221
1222 d(ie^s,1:nw_recon) = w(ie^s,1:nw_recon)-2.0d0*w(iep^s,1:nw_recon)+w(iepp^s,1:nw_recon)
1223
1224 do iw = 1,nw_recon
1225 a(id^s) = 4.0d0*d(id^s,iw)-d(idm^s,iw)
1226 b(id^s) = 4.0d0*d(idm^s,iw)-d(id^s,iw)
1227 call minmod(ixi^l,id^l,a,b,tmp)
1228 a(id^s) = d(id^s,iw)
1229 b(id^s) = d(idm^s,iw)
1230 call minmod(ixi^l,id^l,a,b,tmp2)
1231 call minmod(ixi^l,id^l,tmp,tmp2,tmp3)
1232 dm4(id^s,iw) = tmp3(id^s)
1233 end do
1234
1235 flux_ul(il^s,1:nw_recon) = w(ilp^s,1:nw_recon) + mpalpha * (w(ilp^s,1:nw_recon) - w(ilpp^s,1:nw_recon))
1236 flux_md(il^s,1:nw_recon) = half * (w(il^s,1:nw_recon) + w(ilp^s,1:nw_recon) - dm4(il^s,1:nw_recon))
1237 flux_lc(il^s,1:nw_recon) = half * (3.d0 * w(ilp^s,1:nw_recon) - w(ilpp^s,1:nw_recon)) + mpbeta / 3.d0 * dm4(ilp^s,1:nw_recon)
1238
1239 flux_min(il^s,1:nw_recon) = max(min(w(ilp^s,1:nw_recon), w(il^s,1:nw_recon), flux_md(il^s,1:nw_recon)), &
1240 min(w(ilp^s,1:nw_recon), flux_ul(il^s,1:nw_recon),flux_lc(il^s,1:nw_recon)))
1241
1242 flux_max(il^s,1:nw_recon) = min(max(w(ilp^s,1:nw_recon), w(il^s,1:nw_recon), flux_md(il^s,1:nw_recon)), &
1243 max(w(ilp^s,1:nw_recon), flux_ul(il^s,1:nw_recon),flux_lc(il^s,1:nw_recon)))
1244 do iw=1,nw_recon
1245 a(il^s) = flux(il^s,iw)
1246 b(il^s) = flux_min(il^s,iw)
1247 c(il^s) = flux_max(il^s,iw)
1248 call median(ixi^l, il^l, a, b, c, tmp)
1249 wrctmp(il^s,iw) = tmp(il^s)
1250 end do
1251 end select
1252
1253 call fix_limiter(ixi^l,il^l,wlctmp,wrctmp,wlc,wrc)
1254
1255 end subroutine weno7limiter
1256
1257 subroutine minmod(ixI^L,ixO^L,a,b,minm)
1258
1259 use mod_global_parameters
1260
1261 integer, intent(in) :: ixi^l, ixo^l
1262 double precision, intent(in) :: a(ixi^s), b(ixi^s)
1263 double precision, intent(out):: minm(ixi^s)
1264
1265 minm(ixo^s) = (sign(one,a(ixo^s))+sign(one,b(ixo^s)))/2.0d0 * &
1266 min(abs(a(ixo^s)),abs(b(ixo^s)))
1267
1268 end subroutine minmod
1269
1270 subroutine median(ixI^L,ixO^L,a,b,c,med)
1271
1272 use mod_global_parameters
1273
1274 integer, intent(in) :: ixi^l, ixo^l
1275 double precision, intent(in) :: a(ixi^s), b(ixi^s), c(ixi^s)
1276 double precision, intent(out):: med(ixi^s)
1277
1278 double precision :: tmp1(ixi^s),tmp2(ixi^s)
1279
1280 tmp1(ixo^s) = b(ixo^s) - a(ixo^s); tmp2(ixo^s) = c(ixo^s) - a(ixo^s)
1281
1282 med(ixo^s) = a(ixo^s) + (sign(one,tmp1(ixo^s))+sign(one,tmp2(ixo^s)))/2.0d0 * &
1283 min(abs(tmp1(ixo^s)),abs(tmp2(ixo^s)))
1284
1285 end subroutine median
1286end module mod_weno
mod_global_parameters
This module contains definitions of global parameters and variables and some generic functions/subrou...
Definition mod_global_parameters.t:4
mod_global_parameters::block
type(state), pointer block
Block pointer for using one block and its previous state.
Definition mod_global_parameters.t:771
mod_global_parameters::kr
integer, dimension(3, 3) kr
Kronecker delta tensor.
Definition mod_global_parameters.t:426
mod_global_parameters::ndim
integer, parameter ndim
Number of spatial dimensions for grid variables.
Definition mod_global_parameters.t:217
mod_global_parameters::d
double precision, dimension(:), allocatable, parameter d
Definition mod_global_parameters.t:23
mod_global_parameters::l
double precision l
Definition mod_global_parameters.t:16
mod_physics
This module defines the procedures of a physics module. It contains function pointers for the various...
Definition mod_physics.t:4
mod_physics::phys_check_w
procedure(sub_check_w), pointer phys_check_w
Definition mod_physics.t:72
mod_weno
Definition mod_weno.t:1
mod_weno::weno5limiterl
subroutine, public weno5limiterl(ixil, ill, idims, w, wlc, var)
Definition mod_weno.t:625
mod_weno::weno5limiterr
subroutine, public weno5limiterr(ixil, ill, idims, w, wrc, var)
Definition mod_weno.t:709
mod_weno::weno7limiter
subroutine, public weno7limiter(ixil, ill, idims, w, wlc, wrc, var)
Definition mod_weno.t:1051
mod_weno::fix_onelimiter1
subroutine, public fix_onelimiter1(ixil, ill, wcin, wcout)
Definition mod_weno.t:74
mod_weno::weno3limiter
subroutine, public weno3limiter(ixil, ill, idims, dxdim, w, wlc, wrc, var)
Definition mod_weno.t:166
mod_weno::weno5limiter
subroutine, public weno5limiter(ixil, ill, idims, dxdim, w, wlc, wrc, var)
Definition mod_weno.t:255
mod_weno::teno5adlimiter
subroutine, public teno5adlimiter(ixil, ill, idims, dxdim, w, wlc, wrc)
Definition mod_weno.t:380
mod_weno::weno5nmlimiterl
subroutine, public weno5nmlimiterl(ixil, ill, idims, w, wlc, var)
Definition mod_weno.t:792
mod_weno::fix_limiter
subroutine, public fix_limiter(ixil, ill, wlcin, wrcin, wlcout, wrcout)
Definition mod_weno.t:101
mod_weno::weno5nmlimiterr
subroutine, public weno5nmlimiterr(ixil, ill, idims, w, wrc, var)
Definition mod_weno.t:875
mod_weno::weno5cu6limiter
subroutine, public weno5cu6limiter(ixil, ill, idims, w, wlc, wrc)
Definition mod_weno.t:958
mod_weno::fix_onelimiter
subroutine, public fix_onelimiter(ixil, ill, wcin, wcout)
Definition mod_weno.t:42
mod_weno::fix_limiter1
subroutine, public fix_limiter1(ixil, ill, wlcin, wrcin, wlcout, wrcout)
Definition mod_weno.t:136
mod_weno::weno5nmlimiter
subroutine, public weno5nmlimiter(ixil, ill, idims, dxdim, w, wlc, wrc, var)
Definition mod_weno.t:495
Generated by doxygen 1.9.8