MPI-AMRVAC  3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
mod_weno.t
Go to the documentation of this file.
1 module 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 
39 contains
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,nwflux
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,nwflux
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,nwflux
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,nwflux
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,nwflux
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,nwflux
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:nwflux,1) = u1_coeff(1) * w(ilm^s,1:nwflux) + u1_coeff(2) * w(il^s,1:nwflux)
190  f_array(il^s,1:nwflux,2) = u2_coeff(1) * w(il^s,1:nwflux) + u2_coeff(2) * w(ilp^s,1:nwflux)
191 
192  beta(il^s,1:nwflux,1) = (w(il^s,1:nwflux) - w(ilm^s,1:nwflux))**2
193  beta(il^s,1:nwflux,2) = (w(ilp^s,1:nwflux) - w(il^s,1:nwflux))**2
194 
195  alpha_sum(il^s,1:nwflux) = 0.0d0
196  select case(var)
197  case(1)
198  do i = 1,2
199  alpha_array(il^s,1:nwflux,i) = d_array(i)/(beta(il^s,1:nwflux,i) + weno_eps_machine)**2
200  alpha_sum(il^s,1:nwflux) = alpha_sum(il^s,1:nwflux) + alpha_array(il^s,1:nwflux,i)
201  end do
202  case(2)
203  tau(il^s,1:nwflux) = abs(beta(il^s,1:nwflux,2) - beta(il^s,1:nwflux,1))
204  do i = 1,2
205  alpha_array(il^s,1:nwflux,i) = d_array(i) * (1.d0 + (tau(il^s,1:nwflux) / &
206  (beta(il^s,1:nwflux,i) + dxdim**2)))
207  alpha_sum(il^s,1:nwflux) = alpha_sum(il^s,1:nwflux) + alpha_array(il^s,1:nwflux,i)
208  end do
209  end select
210 
211  flux(il^s,1:nwflux) = 0.0d0
212  do i = 1,2
213  flux(il^s,1:nwflux) = flux(il^s,1:nwflux) + f_array(il^s,1:nwflux,i) * alpha_array(il^s,1:nwflux,i)/(alpha_sum(il^s,1:nwflux))
214  end do
215 
216  !> left value at right interface
217  wlctmp(il^s,1:nwflux) = flux(il^s,1:nwflux)
218 
219  !> right side
220  f_array(il^s,1:nwflux,1) = u1_coeff(1) * w(ilpp^s,1:nwflux) + u1_coeff(2) * w(ilp^s,1:nwflux)
221  f_array(il^s,1:nwflux,2) = u2_coeff(1) * w(ilp^s,1:nwflux) + u2_coeff(2) * w(il^s,1:nwflux)
222 
223  beta(il^s,1:nwflux,1) = (w(ilpp^s,1:nwflux) - w(ilp^s,1:nwflux))**2
224  beta(il^s,1:nwflux,2) = (w(ilp^s,1:nwflux) - w(il^s,1:nwflux))**2
225 
226  alpha_sum(il^s,1:nwflux) = 0.0d0
227  select case(var)
228  case(1)
229  do i = 1,2
230  alpha_array(il^s,1:nwflux,i) = d_array(i)/(beta(il^s,1:nwflux,i) + weno_eps_machine)**2
231  alpha_sum(il^s,1:nwflux) = alpha_sum(il^s,1:nwflux) + alpha_array(il^s,1:nwflux,i)
232  end do
233  case(2)
234  tau(il^s,1:nwflux) = abs(beta(il^s,1:nwflux,2) - beta(il^s,1:nwflux,1))
235  do i = 1,2
236  alpha_array(il^s,1:nwflux,i) = d_array(i) * (1.d0 + (tau(il^s,1:nwflux) / &
237  (beta(il^s,1:nwflux,i) + dxdim**2)))
238  alpha_sum(il^s,1:nwflux) = alpha_sum(il^s,1:nwflux) + alpha_array(il^s,1:nwflux,i)
239  end do
240  end select
241 
242  flux(il^s,1:nwflux) = 0.0d0
243  do i = 1,2
244  flux(il^s,1:nwflux) = flux(il^s,1:nwflux) + f_array(il^s,1:nwflux,i) * alpha_array(il^s,1:nwflux,i)/(alpha_sum(il^s,1:nwflux))
245  end do
246 
247  !> right value at right interface
248  wrctmp(il^s,1:nwflux) = flux(il^s,1:nwflux)
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:nwflux,1) = u1_coeff(1) * w(ilmm^s,1:nwflux) + u1_coeff(2) * w(ilm^s,1:nwflux) + u1_coeff(3) * w(il^s,1:nwflux)
291  f_array(il^s,1:nwflux,2) = u2_coeff(1) * w(ilm^s,1:nwflux) + u2_coeff(2) * w(il^s,1:nwflux) + u2_coeff(3) * w(ilp^s,1:nwflux)
292  f_array(il^s,1:nwflux,3) = u3_coeff(1) * w(il^s,1:nwflux) + u3_coeff(2) * w(ilp^s,1:nwflux) + u3_coeff(3) * w(ilpp^s,1:nwflux)
293 
294  beta(il^s,1:nwflux,1) = beta_coeff(1) * (w(ilmm^s,1:nwflux) + w(il^s,1:nwflux) - 2.0d0*w(ilm^s,1:nwflux))**2 &
295  + beta_coeff(2) * (w(ilmm^s,1:nwflux) - 4.0d0 * w(ilm^s,1:nwflux) + 3.0d0*w(il^s,1:nwflux))**2
296  beta(il^s,1:nwflux,2) = beta_coeff(1) * (w(ilm^s,1:nwflux) + w(ilp^s,1:nwflux) - 2.0d0 * w(il^s,1:nwflux))**2 &
297  + beta_coeff(2) * (w(ilm^s,1:nwflux) - w(ilp^s,1:nwflux))**2
298  beta(il^s,1:nwflux,3) = beta_coeff(1) * (w(il^s,1:nwflux) + w(ilpp^s,1:nwflux) - 2.0d0 * w(ilp^s,1:nwflux))**2 &
299  + beta_coeff(2) * (3.0d0 * w(il^s, 1:nwflux) - 4.0d0 * w(ilp^s,1:nwflux) + w(ilpp^s,1:nwflux))**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:nwflux,i) = d_array(i)/(beta(il^s,1:nwflux,i) + weno_eps_machine)**2
306  end do
307  case(2)
308  tau(il^s,1:nwflux) = abs(beta(il^s,1:nwflux,1) - beta(il^s,1:nwflux,3))
309  do i = 1,3
310  alpha_array(il^s,1:nwflux,i) = d_array(i) * (1.d0 + (tau(il^s,1:nwflux) / &
311  (beta(il^s,1:nwflux,i) + weno_eps_machine))**2)
312  end do
313  case(3)
314  tau(il^s,1:nwflux) = abs(beta(il^s,1:nwflux,1) - beta(il^s,1:nwflux,3))
315  do i = 1,3
316  tmp(il^s,1:nwflux) = (tau(il^s,1:nwflux) + weno_eps_machine) / (beta(il^s,1:nwflux,i) + weno_eps_machine)
317  alpha_array(il^s,1:nwflux,i) = d_array(i) * (1.0d0 + tmp(il^s,1:nwflux)**2 + lambda/tmp(il^s,1:nwflux))
318  end do
319  end select
320 
321  alpha_sum(il^s,1:nwflux) = 0.0d0
322  do i = 1,3
323  alpha_sum(il^s,1:nwflux) = alpha_sum(il^s,1:nwflux) + alpha_array(il^s,1:nwflux,i)
324  end do
325  flux(il^s,1:nwflux) = 0.0d0
326  do i = 1,3
327  flux(il^s,1:nwflux) = flux(il^s,1:nwflux) + f_array(il^s,1:nwflux,i) &
328  *alpha_array(il^s,1:nwflux,i)/(alpha_sum(il^s,1:nwflux))
329  end do
330  wlctmp(il^s,1:nwflux) = flux(il^s,1:nwflux)
331 
332  !> right side
333  f_array(il^s,1:nwflux,1) = u1_coeff(1) * w(ilppp^s,1:nwflux) + u1_coeff(2) * w(ilpp^s,1:nwflux) + u1_coeff(3) * w(ilp^s,1:nwflux)
334  f_array(il^s,1:nwflux,2) = u2_coeff(1) * w(ilpp^s,1:nwflux) + u2_coeff(2) * w(ilp^s,1:nwflux) + u2_coeff(3) * w(il^s,1:nwflux)
335  f_array(il^s,1:nwflux,3) = u3_coeff(1) * w(ilp^s,1:nwflux) + u3_coeff(2) * w(il^s,1:nwflux) + u3_coeff(3) * w(ilm^s,1:nwflux)
336 
337  beta(il^s,1:nwflux,1) = beta_coeff(1) * (w(ilppp^s,1:nwflux) + w(ilp^s,1:nwflux) - 2.0d0*w(ilpp^s,1:nwflux))**2 &
338  + beta_coeff(2) * (w(ilppp^s,1:nwflux) - 4.0d0 * w(ilpp^s,1:nwflux) + 3.0d0*w(ilp^s,1:nwflux))**2
339  beta(il^s,1:nwflux,2) = beta_coeff(1) * (w(ilpp^s,1:nwflux) + w(il^s,1:nwflux) - 2.0d0 * w(ilp^s,1:nwflux))**2 &
340  + beta_coeff(2) * (w(ilpp^s,1:nwflux) - w(il^s,1:nwflux))**2
341  beta(il^s,1:nwflux,3) = beta_coeff(1) * (w(ilp^s,1:nwflux) + w(ilm^s,1:nwflux) - 2.0d0 * w(il^s,1:nwflux))**2 &
342  + beta_coeff(2) * (3.0d0 * w(ilp^s, 1:nwflux) - 4.0d0 * w(il^s,1:nwflux) + w(ilm^s,1:nwflux))**2
343 
344  select case(var)
345  case(1)
346  do i = 1,3
347  alpha_array(il^s,1:nwflux,i) = d_array(i)/(beta(il^s,1:nwflux,i) + weno_eps_machine)**2
348  end do
349  case(2)
350  tau(il^s,1:nwflux) = abs(beta(il^s,1:nwflux,1) - beta(il^s,1:nwflux,3))
351  do i = 1,3
352  alpha_array(il^s,1:nwflux,i) = d_array(i) * (1.d0 + (tau(il^s,1:nwflux) / &
353  (beta(il^s,1:nwflux,i) + weno_eps_machine))**2)
354  end do
355  case(3)
356  tau(il^s,1:nwflux) = abs(beta(il^s,1:nwflux,1) - beta(il^s,1:nwflux,3))
357  do i = 1,3
358  tmp(il^s,1:nwflux) = (tau(il^s,1:nwflux) + weno_eps_machine) / (beta(il^s,1:nwflux,i) + weno_eps_machine)
359  alpha_array(il^s,1:nwflux,i) = d_array(i) * (1.0d0 + tmp(il^s,1:nwflux)**2 + lambda/tmp(il^s,1:nwflux))
360  end do
361  end select
362 
363  alpha_sum(il^s,1:nwflux) = 0.0d0
364  ! do nothing, normal case
365  do i = 1,3
366  alpha_sum(il^s,1:nwflux) = alpha_sum(il^s,1:nwflux) + alpha_array(il^s,1:nwflux,i)
367  end do
368  flux(il^s,1:nwflux) = 0.0d0
369  do i = 1,3
370  flux(il^s,1:nwflux) = flux(il^s,1:nwflux) + f_array(il^s,1:nwflux,i) &
371  *alpha_array(il^s,1:nwflux,i)/(alpha_sum(il^s,1:nwflux))
372  end do
373  wrctmp(il^s,1:nwflux) = flux(il^s,1:nwflux)
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:nwflux,1) = u1_coeff(1) * w(ilmm^s,1:nwflux) + u1_coeff(2) * w(ilm^s,1:nwflux) + u1_coeff(3) * w(il^s,1:nwflux)
409  f_array(il^s,1:nwflux,2) = u2_coeff(1) * w(ilm^s,1:nwflux) + u2_coeff(2) * w(il^s,1:nwflux) + u2_coeff(3) * w(ilp^s,1:nwflux)
410  f_array(il^s,1:nwflux,3) = u3_coeff(1) * w(il^s,1:nwflux) + u3_coeff(2) * w(ilp^s,1:nwflux) + u3_coeff(3) * w(ilpp^s,1:nwflux)
411 
412  beta(il^s,1:nwflux,1) = beta_coeff(1) * (w(ilmm^s,1:nwflux) + w(il^s,1:nwflux) - 2.0d0*w(ilm^s,1:nwflux))**2 &
413  + beta_coeff(2) * (w(ilmm^s,1:nwflux) - 4.0d0 * w(ilm^s,1:nwflux) + 3.0d0*w(il^s,1:nwflux))**2
414  beta(il^s,1:nwflux,2) = beta_coeff(1) * (w(ilm^s,1:nwflux) + w(ilp^s,1:nwflux) - 2.0d0 * w(il^s,1:nwflux))**2 &
415  + beta_coeff(2) * (w(ilm^s,1:nwflux) - w(ilp^s,1:nwflux))**2
416  beta(il^s,1:nwflux,3) = beta_coeff(1) * (w(il^s,1:nwflux) + w(ilpp^s,1:nwflux) - 2.0d0 * w(ilp^s,1:nwflux))**2 &
417  + beta_coeff(2) * (3.0d0 * w(il^s, 1:nwflux) - 4.0d0 * w(ilp^s,1:nwflux) + w(ilpp^s,1:nwflux))**2
418 
419  gam_sum(il^s,1:nwflux) = 0.0d0
420  tau(il^s,1:nwflux) = abs(beta(il^s,1:nwflux,1) - beta(il^s,1:nwflux,3))
421  do i = 1,3
422  kai1(il^s,1:nwflux,i) = (tau(il^s,1:nwflux) / (beta(il^s,1:nwflux,i) + weno_eps_machine))
423  gam(il^s,1:nwflux,i) = (1.d0 + kai1(il^s,1:nwflux,i))**2
424  gam_sum(il^s,1:nwflux) = gam_sum(il^s,1:nwflux) + gam(il^s,1:nwflux,i)
425  end do
426  theta(il^s,1:nwflux) = one / (one + maxval(kai1(il^s,1:nwflux,1:3)/10.d0, dim=ndim+2))
427  marray(il^s,1:nwflux)=-floor(4.d0 + theta(il^s,1:nwflux)*6.d0)
428  do i = 1,3
429  kai(il^s,1:nwflux,i) = gam(il^s,1:nwflux,i) / gam_sum(il^s,1:nwflux)
430  where(kai(il^s,1:nwflux,i) .lt. 10**marray(il^s,1:nwflux))
431  delta(il^s,1:nwflux,i)=zero
432  else where
433  delta(il^s,1:nwflux,i)=one
434  end where
435  end do
436  delta_sum=zero
437  do i = 1,3
438  delta_sum(il^s,1:nwflux)=delta_sum(il^s,1:nwflux)+delta(il^s,1:nwflux,i)*d_array(i)
439  end do
440  flux(il^s,1:nwflux)=0.0d0
441  do i = 1,3
442  flux(il^s,1:nwflux)=flux(il^s,1:nwflux)+f_array(il^s,1:nwflux,i)*delta(il^s,1:nwflux,i)*d_array(i)/(delta_sum(il^s,1:nwflux))
443  end do
444 
445  !> left value at right interface
446  wlctmp(il^s,1:nwflux) = flux(il^s,1:nwflux)
447 
448  !> right side
449  f_array(il^s,1:nwflux,1) = u1_coeff(1) * w(ilppp^s,1:nwflux) + u1_coeff(2) * w(ilpp^s,1:nwflux) + u1_coeff(3) * w(ilp^s,1:nwflux)
450  f_array(il^s,1:nwflux,2) = u2_coeff(1) * w(ilpp^s,1:nwflux) + u2_coeff(2) * w(ilp^s,1:nwflux) + u2_coeff(3) * w(il^s,1:nwflux)
451  f_array(il^s,1:nwflux,3) = u3_coeff(1) * w(ilp^s,1:nwflux) + u3_coeff(2) * w(il^s,1:nwflux) + u3_coeff(3) * w(ilm^s,1:nwflux)
452 
453  beta(il^s,1:nwflux,1) = beta_coeff(1) * (w(ilppp^s,1:nwflux) + w(ilp^s,1:nwflux) - 2.0d0*w(ilpp^s,1:nwflux))**2 &
454  + beta_coeff(2) * (w(ilppp^s,1:nwflux) - 4.0d0 * w(ilpp^s,1:nwflux) + 3.0d0*w(ilp^s,1:nwflux))**2
455  beta(il^s,1:nwflux,2) = beta_coeff(1) * (w(ilpp^s,1:nwflux) + w(il^s,1:nwflux) - 2.0d0 * w(ilp^s,1:nwflux))**2 &
456  + beta_coeff(2) * (w(ilpp^s,1:nwflux) - w(il^s,1:nwflux))**2
457  beta(il^s,1:nwflux,3) = beta_coeff(1) * (w(ilp^s,1:nwflux) + w(ilm^s,1:nwflux) - 2.0d0 * w(il^s,1:nwflux))**2 &
458  + beta_coeff(2) * (3.0d0 * w(ilp^s, 1:nwflux) - 4.0d0 * w(il^s,1:nwflux) + w(ilm^s,1:nwflux))**2
459 
460 
461  gam_sum(il^s,1:nwflux)=0.0d0
462  tau(il^s,1:nwflux)=abs(beta(il^s,1:nwflux,1)-beta(il^s,1:nwflux,3))
463  do i=1,3
464  kai1(il^s,1:nwflux,i)=(tau(il^s,1:nwflux)/(beta(il^s,1:nwflux,i)+weno_eps_machine))
465  gam(il^s,1:nwflux,i)=(1.d0+kai1(il^s,1:nwflux,i))**6
466  gam_sum(il^s,1:nwflux)=gam_sum(il^s,1:nwflux)+gam(il^s,1:nwflux,i)
467  end do
468  theta(il^s,1:nwflux)=one/(one+maxval(kai1(il^s,1:nwflux,1:3)/10.d0,dim=ndim+2))
469  marray(il^s,1:nwflux)=-floor(4.d0+theta(il^s,1:nwflux)*6.d0)
470  do i=1,3
471  kai(il^s,1:nwflux,i) = gam(il^s,1:nwflux,i)/gam_sum(il^s,1:nwflux)
472  where(kai(il^s,1:nwflux,i) .lt. 10**marray(il^s,1:nwflux))
473  delta(il^s,1:nwflux,i)=zero
474  else where
475  delta(il^s,1:nwflux,i)=one
476  end where
477  end do
478  delta_sum=zero
479  do i = 1,3
480  delta_sum(il^s,1:nwflux)=delta_sum(il^s,1:nwflux)+delta(il^s,1:nwflux,i)*d_array(i)
481  end do
482  flux(il^s,1:nwflux)=0.0d0
483  do i = 1,3
484  flux(il^s,1:nwflux)=flux(il^s,1:nwflux)+f_array(il^s,1:nwflux,i)*delta(il^s,1:nwflux,i)*d_array(i)/(delta_sum(il^s,1:nwflux))
485  end do
486 
487  !> right value at right interface
488  wrctmp(il^s,1:nwflux)=flux(il^s,1:nwflux)
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, nwflux
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:nwflux,1) = u1_coeff(1) * wd(il^s,1:nwflux) + u1_coeff(2) * wc(ilm^s,1:nwflux)+ u1_coeff(3) * w(il^s,1:nwflux)
538  f_array(il^s,1:nwflux,2) = u2_coeff(1) * wc(ilm^s,1:nwflux) + u2_coeff(2) * w(il^s,1:nwflux) + u2_coeff(3) * wc(il^s,1:nwflux)
539  f_array(il^s,1:nwflux,3) = u3_coeff(1) * wc(il^s,1:nwflux) + u3_coeff(2) * w(ilp^s,1:nwflux) + u3_coeff(3) * wc(ilp^s,1:nwflux)
540 
541  beta(il^s,1:nwflux,1) = beta_coeff(1) * (wc(ilm^s,1:nwflux) - wd(il^s,1:nwflux))**2 &
542  + beta_coeff(2) * (2.d0 * w(il^s,1:nwflux) - wc(ilm^s,1:nwflux) - wd(il^s,1:nwflux))**2
543  beta(il^s,1:nwflux,2) = beta_coeff(1) * (wc(ilm^s,1:nwflux) + wc(il^s,1:nwflux) - 2.0d0 * w(il^s,1:nwflux))**2 &
544  + beta_coeff(2) * (wc(ilm^s,1:nwflux) - wc(il^s,1:nwflux))**2
545  beta(il^s,1:nwflux,3) = beta_coeff(1) * (wc(il^s,1:nwflux) + wc(ilp^s,1:nwflux) - 2.0d0 * w(ilp^s,1:nwflux))**2 &
546  + beta_coeff(2) * (3.0d0 * wc(il^s, 1:nwflux) - 4.0d0 * w(ilp^s,1:nwflux) + wc(ilp^s,1:nwflux))**2
547  alpha_sum(il^s,1:nwflux) = 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:nwflux,i) = d_array(i)/(4.d0 * beta(il^s,1:nwflux,i) + weno_eps_machine)**2
553  alpha_sum(il^s,1:nwflux) = alpha_sum(il^s,1:nwflux) + alpha_array(il^s,1:nwflux,i)
554  end do
555  case(2)
556  tau(il^s,1:nwflux) = abs(beta(il^s,1:nwflux,1) - beta(il^s,1:nwflux,3)) * 4.d0
557  do i = 1,3
558  alpha_array(il^s,1:nwflux,i) = d_array(i) * (1.d0 + (tau(il^s,1:nwflux) / &
559  (4.d0 * beta(il^s,1:nwflux,i) + weno_eps_machine))**2)
560  alpha_sum(il^s,1:nwflux) = alpha_sum(il^s,1:nwflux) + alpha_array(il^s,1:nwflux,i)
561  end do
562  case(3)
563  tau(il^s,1:nwflux) = abs(beta(il^s,1:nwflux,1) - beta(il^s,1:nwflux,3)) * 4.d0
564  do i=1,3
565  do j=1,nwflux
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:nwflux) = 0.0d0
573  do i = 1,3
574  flux(il^s,1:nwflux) = flux(il^s,1:nwflux) + f_array(il^s,1:nwflux,i) * alpha_array(il^s,1:nwflux,i)/(alpha_sum(il^s,1:nwflux))
575  end do
576  !> left value at right interface
577  wlctmp(il^s,1:nwflux) = flux(il^s,1:nwflux)
578 
579  !> right side
580  f_array(il^s,1:nwflux,1) = u1_coeff(1) * we(il^s,1:nwflux) + u1_coeff(2) * wc(ilp^s,1:nwflux) + u1_coeff(3) * w(ilp^s,1:nwflux)
581  f_array(il^s,1:nwflux,2) = u2_coeff(1) * wc(ilp^s,1:nwflux) + u2_coeff(2) * w(ilp^s,1:nwflux) + u2_coeff(3) * wc(il^s,1:nwflux)
582  f_array(il^s,1:nwflux,3) = u3_coeff(1) * wc(il^s,1:nwflux) + u3_coeff(2) * w(il^s,1:nwflux) + u3_coeff(3) * wc(ilm^s,1:nwflux)
583  beta(il^s,1:nwflux,1) = beta_coeff(1) * (wc(ilp^s,1:nwflux) - we(il^s,1:nwflux))**2 &
584  + beta_coeff(2) * (2.d0 * w(ilp^s,1:nwflux) - wc(ilp^s,1:nwflux) - we(il^s,1:nwflux))**2
585  beta(il^s,1:nwflux,2) = beta_coeff(1) * (wc(ilp^s,1:nwflux) + wc(il^s,1:nwflux) - 2.0d0 * w(ilp^s,1:nwflux))**2 &
586  + beta_coeff(2) * (wc(ilp^s,1:nwflux) - wc(il^s,1:nwflux))**2
587  beta(il^s,1:nwflux,3) = beta_coeff(1) * (wc(il^s,1:nwflux) + wc(ilm^s,1:nwflux) - 2.0d0 * w(il^s,1:nwflux))**2 &
588  + beta_coeff(2) * (3.0d0 * wc(il^s, 1:nwflux) - 4.0d0 * w(il^s,1:nwflux) + wc(ilm^s,1:nwflux))**2
589  alpha_sum(il^s,1:nwflux) = 0.0d0
590  select case(var)
591  case(1)
592  do i = 1,3
593  alpha_array(il^s,1:nwflux,i) = d_array(i)/(4.d0 * beta(il^s,1:nwflux,i) + weno_eps_machine)**2
594  alpha_sum(il^s,1:nwflux) = alpha_sum(il^s,1:nwflux) + alpha_array(il^s,1:nwflux,i)
595  end do
596  case(2)
597  tau(il^s,1:nwflux) = abs(beta(il^s,1:nwflux,1) - beta(il^s,1:nwflux,3)) * 4.d0
598  do i = 1,3
599  alpha_array(il^s,1:nwflux,i) = d_array(i) * (1.d0 + (tau(il^s,1:nwflux) / &
600  (4.d0 * beta(il^s,1:nwflux,i) + weno_eps_machine))**2)
601  alpha_sum(il^s,1:nwflux) = alpha_sum(il^s,1:nwflux) + alpha_array(il^s,1:nwflux,i)
602  end do
603  case(3)
604  tau(il^s,1:nwflux) = abs(beta(il^s,1:nwflux,1) - beta(il^s,1:nwflux,3)) * 4.d0
605  do i = 1,3
606  do j = 1,nwflux
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:nwflux) = 0.0d0
614  do i = 1,3
615  flux(il^s,1:nwflux) = flux(il^s,1:nwflux) + f_array(il^s,1:nwflux,i) * alpha_array(il^s,1:nwflux,i)/(alpha_sum(il^s,1:nwflux))
616  end do
617  !> right value at right interface
618  wrctmp(il^s,1:nwflux) = flux(il^s,1:nwflux)
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:nwflux,1) = u1_coeff(1) * w(ilmm^s,1:nwflux) + u1_coeff(2) * w(ilm^s,1:nwflux) + u1_coeff(3) * w(il^s,1:nwflux)
661  f_array(il^s,1:nwflux,2) = u2_coeff(1) * w(ilm^s,1:nwflux) + u2_coeff(2) * w(il^s,1:nwflux) + u2_coeff(3) * w(ilp^s,1:nwflux)
662  f_array(il^s,1:nwflux,3) = u3_coeff(1) * w(il^s,1:nwflux) + u3_coeff(2) * w(ilp^s,1:nwflux) + u3_coeff(3) * w(ilpp^s,1:nwflux)
663 
664  beta(il^s,1:nwflux,1) = beta_coeff(1) * (w(ilmm^s,1:nwflux) + w(il^s,1:nwflux) - 2.0d0*w(ilm^s,1:nwflux))**2 &
665  + beta_coeff(2) * (w(ilmm^s,1:nwflux) - 4.0d0 * w(ilm^s,1:nwflux) + 3.0d0*w(il^s,1:nwflux))**2
666  beta(il^s,1:nwflux,2) = beta_coeff(1) * (w(ilm^s,1:nwflux) + w(ilp^s,1:nwflux) - 2.0d0 * w(il^s,1:nwflux))**2 &
667  + beta_coeff(2) * (w(ilm^s,1:nwflux) - w(ilp^s,1:nwflux))**2
668  beta(il^s,1:nwflux,3) = beta_coeff(1) * (w(il^s,1:nwflux) + w(ilpp^s,1:nwflux) - 2.0d0 * w(ilp^s,1:nwflux))**2 &
669  + beta_coeff(2) * (3.0d0 * w(il^s, 1:nwflux) - 4.0d0 * w(ilp^s,1:nwflux) + w(ilpp^s,1:nwflux))**2
670 
671  alpha_sum(il^s,1:nwflux) = 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:nwflux,i) = d_array(i)/(beta(il^s,1:nwflux,i) + weno_eps_machine)**2
677  alpha_sum(il^s,1:nwflux) = alpha_sum(il^s,1:nwflux) + alpha_array(il^s,1:nwflux,i)
678  end do
679  case(2)
680  tau(il^s,1:nwflux) = abs(beta(il^s,1:nwflux,1) - beta(il^s,1:nwflux,3))
681  do i = 1,3
682  alpha_array(il^s,1:nwflux,i) = d_array(i) * (1.d0 + (tau(il^s,1:nwflux) / &
683  (beta(il^s,1:nwflux,i) + weno_eps_machine))**2)
684  alpha_sum(il^s,1:nwflux) = alpha_sum(il^s,1:nwflux) + alpha_array(il^s,1:nwflux,i)
685  end do
686  case(3)
687  tau(il^s,1:nwflux) = abs(beta(il^s,1:nwflux,1) - beta(il^s,1:nwflux,3)) * 4.d0
688  do i=1,3
689  do j=1,nwflux
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:nwflux) = 0.0d0
697  do i = 1,3
698  flux(il^s,1:nwflux) = flux(il^s,1:nwflux) + f_array(il^s,1:nwflux,i) * alpha_array(il^s,1:nwflux,i)/(alpha_sum(il^s,1:nwflux))
699  end do
700 
701  !> left value at right interface
702  wlctmp(il^s,1:nwflux) = flux(il^s,1:nwflux)
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:nwflux,1) = u1_coeff(1) * w(ilppp^s,1:nwflux) + u1_coeff(2) * w(ilpp^s,1:nwflux) + u1_coeff(3) * w(ilp^s,1:nwflux)
745  f_array(il^s,1:nwflux,2) = u2_coeff(1) * w(ilpp^s,1:nwflux) + u2_coeff(2) * w(ilp^s,1:nwflux) + u2_coeff(3) * w(il^s,1:nwflux)
746  f_array(il^s,1:nwflux,3) = u3_coeff(1) * w(ilp^s,1:nwflux) + u3_coeff(2) * w(il^s,1:nwflux) + u3_coeff(3) * w(ilm^s,1:nwflux)
747 
748  beta(il^s,1:nwflux,1) = beta_coeff(1) * (w(ilppp^s,1:nwflux) + w(ilp^s,1:nwflux) - 2.0d0*w(ilpp^s,1:nwflux))**2 &
749  + beta_coeff(2) * (w(ilppp^s,1:nwflux) - 4.0d0 * w(ilpp^s,1:nwflux) + 3.0d0*w(ilp^s,1:nwflux))**2
750  beta(il^s,1:nwflux,2) = beta_coeff(1) * (w(ilpp^s,1:nwflux) + w(il^s,1:nwflux) - 2.0d0 * w(ilp^s,1:nwflux))**2 &
751  + beta_coeff(2) * (w(ilpp^s,1:nwflux) - w(il^s,1:nwflux))**2
752  beta(il^s,1:nwflux,3) = beta_coeff(1) * (w(ilp^s,1:nwflux) + w(ilm^s,1:nwflux) - 2.0d0 * w(il^s,1:nwflux))**2 &
753  + beta_coeff(2) * (3.0d0 * w(ilp^s, 1:nwflux) - 4.0d0 * w(il^s,1:nwflux) + w(ilm^s,1:nwflux))**2
754 
755  alpha_sum(il^s,1:nwflux) = 0.0d0
756  select case(var)
757  case(1)
758  do i = 1,3
759  alpha_array(il^s,1:nwflux,i) = d_array(i)/(beta(il^s,1:nwflux,i) + weno_eps_machine)**2
760  alpha_sum(il^s,1:nwflux) = alpha_sum(il^s,1:nwflux) + alpha_array(il^s,1:nwflux,i)
761  end do
762  case(2)
763  tau(il^s,1:nwflux) = abs(beta(il^s,1:nwflux,1) - beta(il^s,1:nwflux,3))
764  do i = 1,3
765  alpha_array(il^s,1:nwflux,i) = d_array(i) * (1.d0 + (tau(il^s,1:nwflux) / &
766  (beta(il^s,1:nwflux,i) + weno_eps_machine))**2)
767  alpha_sum(il^s,1:nwflux) = alpha_sum(il^s,1:nwflux) + alpha_array(il^s,1:nwflux,i)
768  end do
769  case(3)
770  tau(il^s,1:nwflux) = abs(beta(il^s,1:nwflux,1) - beta(il^s,1:nwflux,3)) * 4.d0
771  do i = 1,3
772  do j = 1,nwflux
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:nwflux) = 0.0d0
780  do i = 1,3
781  flux(il^s,1:nwflux) = flux(il^s,1:nwflux) + f_array(il^s,1:nwflux,i) * alpha_array(il^s,1:nwflux,i)/(alpha_sum(il^s,1:nwflux))
782  end do
783 
784  !> right value at right interface
785  wrctmp(il^s,1:nwflux) = flux(il^s,1:nwflux)
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, nwflux
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:nwflux,1) = u1_coeff(1) * wd(il^s,1:nwflux) + u1_coeff(2) * wc(ilm^s,1:nwflux)+ u1_coeff(3) * w(il^s,1:nwflux)
830  f_array(il^s,1:nwflux,2) = u2_coeff(1) * wc(ilm^s,1:nwflux) + u2_coeff(2) * w(il^s,1:nwflux) + u2_coeff(3) * wc(il^s,1:nwflux)
831  f_array(il^s,1:nwflux,3) = u3_coeff(1) * wc(il^s,1:nwflux) + u3_coeff(2) * w(ilp^s,1:nwflux) + u3_coeff(3) * wc(ilp^s,1:nwflux)
832  beta(il^s,1:nwflux,1) = beta_coeff(1) * (wc(ilm^s,1:nwflux) - wd(il^s,1:nwflux))**2 &
833  + beta_coeff(2) * (2.d0 * w(il^s,1:nwflux) - wc(ilm^s,1:nwflux) - wd(il^s,1:nwflux))**2
834  beta(il^s,1:nwflux,2) = beta_coeff(1) * (wc(ilm^s,1:nwflux) + wc(il^s,1:nwflux) - 2.0d0 * w(il^s,1:nwflux))**2 &
835  + beta_coeff(2) * (wc(ilm^s,1:nwflux) - wc(il^s,1:nwflux))**2
836  beta(il^s,1:nwflux,3) = beta_coeff(1) * (wc(il^s,1:nwflux) + wc(ilp^s,1:nwflux) - 2.0d0 * w(ilp^s,1:nwflux))**2 &
837  + beta_coeff(2) * (3.0d0 * wc(il^s, 1:nwflux) - 4.0d0 * w(ilp^s,1:nwflux) + wc(ilp^s,1:nwflux))**2
838  alpha_sum(il^s,1:nwflux) = 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:nwflux,i) = d_array(i)/(4.d0 * beta(il^s,1:nwflux,i) + weno_eps_machine)**2
844  alpha_sum(il^s,1:nwflux) = alpha_sum(il^s,1:nwflux) + alpha_array(il^s,1:nwflux,i)
845  end do
846  case(2)
847  tau(il^s,1:nwflux) = abs(beta(il^s,1:nwflux,1) - beta(il^s,1:nwflux,3)) * 4.d0
848  do i = 1,3
849  alpha_array(il^s,1:nwflux,i) = d_array(i) * (1.d0 + (tau(il^s,1:nwflux) / &
850  (4.d0 * beta(il^s,1:nwflux,i) + weno_eps_machine))**2)
851  alpha_sum(il^s,1:nwflux) = alpha_sum(il^s,1:nwflux) + alpha_array(il^s,1:nwflux,i)
852  end do
853  case(3)
854  tau(il^s,1:nwflux) = abs(beta(il^s,1:nwflux,1) - beta(il^s,1:nwflux,3)) * 4.d0
855  do i=1,3
856  do j=1,nwflux
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:nwflux) = 0.0d0
864  do i = 1,3
865  flux(il^s,1:nwflux) = flux(il^s,1:nwflux) + f_array(il^s,1:nwflux,i) * alpha_array(il^s,1:nwflux,i)/(alpha_sum(il^s,1:nwflux))
866  end do
867  !> left value at right interface
868  wlctmp(il^s,1:nwflux) = flux(il^s,1:nwflux)
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, nwflux
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:nwflux,1) = u1_coeff(1) * we(il^s,1:nwflux) + u1_coeff(2) * wc(ilp^s,1:nwflux) + u1_coeff(3) * w(ilp^s,1:nwflux)
914  f_array(il^s,1:nwflux,2) = u2_coeff(1) * wc(ilp^s,1:nwflux) + u2_coeff(2) * w(ilp^s,1:nwflux) + u2_coeff(3) * wc(il^s,1:nwflux)
915  f_array(il^s,1:nwflux,3) = u3_coeff(1) * wc(il^s,1:nwflux) + u3_coeff(2) * w(il^s,1:nwflux) + u3_coeff(3) * wc(ilm^s,1:nwflux)
916  beta(il^s,1:nwflux,1) = beta_coeff(1) * (wc(ilp^s,1:nwflux) - we(il^s,1:nwflux))**2 &
917  + beta_coeff(2) * (2.d0 * w(ilp^s,1:nwflux) - wc(ilp^s,1:nwflux) - we(il^s,1:nwflux))**2
918  beta(il^s,1:nwflux,2) = beta_coeff(1) * (wc(ilp^s,1:nwflux) + wc(il^s,1:nwflux) - 2.0d0 * w(ilp^s,1:nwflux))**2 &
919  + beta_coeff(2) * (wc(ilp^s,1:nwflux) - wc(il^s,1:nwflux))**2
920  beta(il^s,1:nwflux,3) = beta_coeff(1) * (wc(il^s,1:nwflux) + wc(ilm^s,1:nwflux) - 2.0d0 * w(il^s,1:nwflux))**2 &
921  + beta_coeff(2) * (3.0d0 * wc(il^s, 1:nwflux) - 4.0d0 * w(il^s,1:nwflux) + wc(ilm^s,1:nwflux))**2
922  alpha_sum(il^s,1:nwflux) = 0.0d0
923  select case(var)
924  case(1)
925  do i = 1,3
926  alpha_array(il^s,1:nwflux,i) = d_array(i)/(4.d0 * beta(il^s,1:nwflux,i) + weno_eps_machine)**2
927  alpha_sum(il^s,1:nwflux) = alpha_sum(il^s,1:nwflux) + alpha_array(il^s,1:nwflux,i)
928  end do
929  case(2)
930  tau(il^s,1:nwflux) = abs(beta(il^s,1:nwflux,1) - beta(il^s,1:nwflux,3)) * 4.d0
931  do i = 1,3
932  alpha_array(il^s,1:nwflux,i) = d_array(i) * (1.d0 + (tau(il^s,1:nwflux) / &
933  (4.d0 * beta(il^s,1:nwflux,i) + weno_eps_machine))**2)
934  alpha_sum(il^s,1:nwflux) = alpha_sum(il^s,1:nwflux) + alpha_array(il^s,1:nwflux,i)
935  end do
936  case(3)
937  tau(il^s,1:nwflux) = abs(beta(il^s,1:nwflux,1) - beta(il^s,1:nwflux,3)) * 4.d0
938  do i = 1,3
939  do j = 1,nwflux
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:nwflux) = 0.0d0
947  do i = 1,3
948  flux(il^s,1:nwflux) = flux(il^s,1:nwflux) + f_array(il^s,1:nwflux,i) * alpha_array(il^s,1:nwflux,i)/(alpha_sum(il^s,1:nwflux))
949  end do
950  !> right value at right interface
951  wrctmp(il^s,1:nwflux) = flux(il^s,1:nwflux)
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:nwflux,1)=beta_coeff(1)*(w(ilmm^s,1:nwflux)+w(il^s,1:nwflux)-2.0d0*w(ilm^s,1:nwflux))**2&
988  +beta_coeff(2)*(w(ilmm^s,1:nwflux)-4.0d0*w(ilm^s,1:nwflux)+3.0d0*w(il^s,1:nwflux))**2
989  beta(il^s,1:nwflux,2)=beta_coeff(1)*(w(ilm^s,1:nwflux)+w(ilp^s,1:nwflux)-2.0d0*w(il^s,1:nwflux))**2&
990  +beta_coeff(2)*(w(ilm^s,1:nwflux)-w(ilp^s,1:nwflux))**2
991  beta(il^s,1:nwflux,3)=beta_coeff(1)*(w(il^s,1:nwflux)+w(ilpp^s,1:nwflux)-2.0d0*w(ilp^s,1:nwflux))**2&
992  +beta_coeff(2)*(3.0d0*w(il^s,1:nwflux)-4.0d0*w(ilp^s,1:nwflux)+w(ilpp^s,1:nwflux))**2
993  alpha_sum(il^s,1:nwflux)=zero
994  do i=1,3
995  alpha_array(il^s,1:nwflux,i)=d_array(i)/(beta(il^s,1:nwflux,i)+weno_eps_machine)**2
996  alpha_sum(il^s,1:nwflux)=alpha_sum(il^s,1:nwflux)+alpha_array(il^s,1:nwflux,i)
997  end do
998  do i=1,3
999  alpha_array(il^s,1:nwflux,i)=alpha_array(il^s,1:nwflux,i)/alpha_sum(il^s,1:nwflux)
1000  end do
1001  theta2(il^s,1:nwflux)=((alpha_array(il^s,1:nwflux,1)/d_array(1)-1.d0)**2&
1002  +(alpha_array(il^s,1:nwflux,2)/d_array(2)-1.d0)**2&
1003  +(alpha_array(il^s,1:nwflux,3)/d_array(3)-1.d0)**2)/83.d0
1004  where(theta2(il^s,1:nwflux) .gt. theta_limit)
1005  f_array(il^s,1:nwflux,1)=u1_coeff(1)*w(ilmm^s,1:nwflux)+u1_coeff(2)*w(ilm^s,1:nwflux)+u1_coeff(3)*w(il^s,1:nwflux)
1006  f_array(il^s,1:nwflux,2)=u2_coeff(1)*w(ilm^s,1:nwflux)+u2_coeff(2)*w(il^s,1:nwflux)+u2_coeff(3)*w(ilp^s,1:nwflux)
1007  f_array(il^s,1:nwflux,3)=u3_coeff(1)*w(il^s,1:nwflux)+u3_coeff(2)*w(ilp^s,1:nwflux)+u3_coeff(3)*w(ilpp^s,1:nwflux)
1008  wlctmp(il^s,1:nwflux)=f_array(il^s,1:nwflux,1)*alpha_array(il^s,1:nwflux,1)&
1009  +f_array(il^s,1:nwflux,2)*alpha_array(il^s,1:nwflux,2)&
1010  +f_array(il^s,1:nwflux,3)*alpha_array(il^s,1:nwflux,3)
1011  else where
1012  wlctmp(il^s,1:nwflux)=1.d0/60.d0*(w(ilmm^s,1:nwflux)-8.d0*w(ilm^s,1:nwflux)+37.d0*w(il^s,1:nwflux)&
1013  +37.d0*w(ilp^s,1:nwflux)-8.d0*w(ilpp^s,1:nwflux)+w(ilppp^s,1:nwflux))
1014  end where
1015 
1016  !> right side
1017  beta(il^s,1:nwflux,1)=beta_coeff(1)*(w(ilppp^s,1:nwflux)+w(ilp^s,1:nwflux)-2.0d0*w(ilpp^s,1:nwflux))**2&
1018  +beta_coeff(2)*(w(ilppp^s,1:nwflux)-4.0d0*w(ilpp^s,1:nwflux)+3.0d0*w(ilp^s,1:nwflux))**2
1019  beta(il^s,1:nwflux,2)=beta_coeff(1)*(w(ilpp^s,1:nwflux)+w(il^s,1:nwflux)-2.0d0*w(ilp^s,1:nwflux))**2&
1020  +beta_coeff(2)*(w(ilpp^s,1:nwflux)-w(il^s,1:nwflux))**2
1021  beta(il^s,1:nwflux,3)=beta_coeff(1)*(w(ilp^s,1:nwflux)+w(ilm^s,1:nwflux)-2.0d0*w(il^s,1:nwflux))**2&
1022  +beta_coeff(2)*(3.0d0*w(ilp^s,1:nwflux)-4.0d0*w(il^s,1:nwflux)+w(ilm^s,1:nwflux))**2
1023  alpha_sum(il^s,1:nwflux)=zero
1024  do i=1,3
1025  alpha_array(il^s,1:nwflux,i)=d_array(i)/(beta(il^s,1:nwflux,i)+weno_eps_machine)**2
1026  alpha_sum(il^s,1:nwflux)=alpha_sum(il^s,1:nwflux)+alpha_array(il^s,1:nwflux,i)
1027  end do
1028  do i=1,3
1029  alpha_array(il^s,1:nwflux,i)=alpha_array(il^s,1:nwflux,i)/alpha_sum(il^s,1:nwflux)
1030  end do
1031  theta2(il^s,1:nwflux)=((alpha_array(il^s,1:nwflux,1)/d_array(1)-1.d0)**2&
1032  +(alpha_array(il^s,1:nwflux,2)/d_array(2)-1.d0)**2&
1033  +(alpha_array(il^s,1:nwflux,3)/d_array(3)-1.d0)**2)/83.d0
1034  where(theta2(il^s,1:nwflux) .gt. theta_limit)
1035  f_array(il^s,1:nwflux,1)=u1_coeff(1)*w(ilppp^s,1:nwflux)+u1_coeff(2)*w(ilpp^s,1:nwflux)+u1_coeff(3)*w(ilp^s,1:nwflux)
1036  f_array(il^s,1:nwflux,2)=u2_coeff(1)*w(ilpp^s,1:nwflux)+u2_coeff(2)*w(ilp^s,1:nwflux)+u2_coeff(3)*w(il^s,1:nwflux)
1037  f_array(il^s,1:nwflux,3)=u3_coeff(1)*w(ilp^s,1:nwflux)+u3_coeff(2)*w(il^s,1:nwflux)+u3_coeff(3)*w(ilm^s,1:nwflux)
1038  wrctmp(il^s,1:nwflux)=f_array(il^s,1:nwflux,1)*alpha_array(il^s,1:nwflux,1)&
1039  +f_array(il^s,1:nwflux,2)*alpha_array(il^s,1:nwflux,2)&
1040  +f_array(il^s,1:nwflux,3)*alpha_array(il^s,1:nwflux,3)
1041  else where
1042  wrctmp(il^s,1:nwflux)=1.d0/60.d0*(w(ilppp^s,1:nwflux)-8.d0*w(ilpp^s,1:nwflux)+37.d0*w(ilp^s,1:nwflux)&
1043  +37.d0*w(il^s,1:nwflux)-8.d0*w(ilm^s,1:nwflux)+w(ilmm^s,1:nwflux))
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:nwflux,1) = u1_coeff(1) * w(ilmmm^s,1:nwflux) + u1_coeff(2) * w(ilmm^s,1:nwflux) + u1_coeff(3) * w(ilm^s,1:nwflux) &
1085  + u1_coeff(4) * w(il^s,1:nwflux)
1086  f_array(il^s,1:nwflux,2) = u2_coeff(1) * w(ilmm^s,1:nwflux) + u2_coeff(2) * w(ilm^s,1:nwflux) + u2_coeff(3) * w(il^s,1:nwflux) &
1087  + u2_coeff(4) * w(ilp^s,1:nwflux)
1088  f_array(il^s,1:nwflux,3) = u3_coeff(1) * w(ilm^s,1:nwflux) + u3_coeff(2) * w(il^s,1:nwflux) + u3_coeff(3) * w(ilp^s,1:nwflux) &
1089  + u3_coeff(4) * w(ilpp^s,1:nwflux)
1090  f_array(il^s,1:nwflux,4) = u4_coeff(1) * w(il^s,1:nwflux) + u4_coeff(2) * w(ilp^s,1:nwflux) + u4_coeff(3) * w(ilpp^s,1:nwflux) &
1091  + u4_coeff(4) * w(ilppp^s,1:nwflux)
1092 
1093  beta(il^s,1:nwflux,1) = w(ilmmm^s,1:nwflux) * (547.d0 * w(ilmmm^s,1:nwflux) - 3882.d0 * w(ilmm^s,1:nwflux) + 4642.d0 * w(ilm^s,1:nwflux) &
1094  - 1854.d0 * w(il^s,1:nwflux)) &
1095  + w(ilmm^s,1:nwflux) * (7043.d0 * w(ilmm^s,1:nwflux) - 17246.d0 * w(ilm^s,1:nwflux) + 7042.d0 * w(il^s,1:nwflux)) &
1096  + w(ilm^s,1:nwflux) * (11003.d0 * w(ilm^s,1:nwflux) - 9402.d0 * w(il^s,1:nwflux)) + 2107.d0 * w(il^s,1:nwflux)**2
1097  beta(il^s,1:nwflux,2) = w(ilmm^s,1:nwflux) * (267.d0 * w(ilmm^s,1:nwflux) - 1642.d0 * w(ilm^s,1:nwflux) + 1602.d0 * w(il^s,1:nwflux) &
1098  - 494.d0 * w(ilp^s,1:nwflux)) &
1099  + w(ilm^s,1:nwflux) * (2843.d0 * w(ilm^s,1:nwflux) - 5966.d0 * w(il^s,1:nwflux) + 1922.d0 * w(ilp^s,1:nwflux)) &
1100  + w(il^s,1:nwflux) * (3443.d0 * w(il^s,1:nwflux) - 2522.d0 * w(ilp^s,1:nwflux)) + 547.d0 * w(ilp^s,1:nwflux) ** 2
1101  beta(il^s,1:nwflux,3) = w(ilm^s,1:nwflux) * (547.d0 * w(ilm^s,1:nwflux) - 2522.d0 * w(il^s,1:nwflux) + 1922.d0 * w(ilp^s,1:nwflux) &
1102  - 494.d0 * w(ilpp^s,1:nwflux)) &
1103  + w(il^s,1:nwflux) * (3443.d0 * w(il^s,1:nwflux) - 5966.d0 * w(ilp^s,1:nwflux) + 1602.d0 * w(ilpp^s,1:nwflux)) &
1104  + w(ilp^s,1:nwflux) * (2843.d0 * w(ilp^s,1:nwflux) - 1642.d0 * w(ilpp^s,1:nwflux)) + 267.d0 * w(ilpp^s,1:nwflux) ** 2
1105  beta(il^s,1:nwflux,4) = w(il^s,1:nwflux) * (2107.d0 * w(il^s,1:nwflux) - 9402.d0 * w(ilp^s,1:nwflux) + 7042.d0 * w(ilpp^s,1:nwflux) &
1106  - 1854.d0 * w(ilppp^s,1:nwflux)) &
1107  + w(ilp^s,1:nwflux) * (11003.d0 * w(ilp^s,1:nwflux) - 17246.d0 * w(ilpp^s,1:nwflux) + 4642.d0 * w(ilppp^s,1:nwflux)) &
1108  + w(ilpp^s,1:nwflux) * (7043.d0 * w(ilpp^s,1:nwflux) - 3882.d0 * w(ilppp^s,1:nwflux)) &
1109  + 547.d0 * w(ilppp^s,1:nwflux) ** 2
1110 
1111  alpha_sum(il^s,1:nwflux) = 0.0d0
1112  do i = 1,4
1113  alpha_array(il^s,1:nwflux,i) = d_array(i)/(beta(il^s,1:nwflux,i) + weno_eps_machine)**2
1114  alpha_sum(il^s,1:nwflux) = alpha_sum(il^s,1:nwflux) + alpha_array(il^s,1:nwflux,i)
1115  end do
1116  flux(il^s,1:nwflux) = 0.0d0
1117  do i = 1,4
1118  flux(il^s,1:nwflux) = flux(il^s,1:nwflux) + f_array(il^s,1:nwflux,i) * alpha_array(il^s,1:nwflux,i)/(alpha_sum(il^s,1:nwflux))
1119  end do
1120 
1121  select case(var)
1122  ! case1 for wenojs, case2 for mpweno
1123  case(1)
1124  wlctmp(il^s,1:nwflux) = flux(il^s,1:nwflux)
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:nwflux) = w(iep^s,1:nwflux)-2.0d0*w(ie^s,1:nwflux)+w(iem^s,1:nwflux)
1135 
1136  do iw=1,nwflux
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:nwflux) = w(il^s,1:nwflux) + mpalpha * (w(il^s,1:nwflux) - w(ilm^s,1:nwflux))
1148  flux_md(il^s,1:nwflux) = half * (w(il^s,1:nwflux) + w(ilp^s,1:nwflux) - dm4(il^s,1:nwflux))
1149  flux_lc(il^s,1:nwflux) = half * (3.d0 * w(il^s,1:nwflux) - w(ilm^s,1:nwflux)) + mpbeta / 3.d0 * dm4(ilm^s,1:nwflux)
1150 
1151  flux_min(il^s,1:nwflux) = max(min(w(il^s,1:nwflux), w(ilp^s,1:nwflux), flux_md(il^s,1:nwflux)), &
1152  min(w(il^s,1:nwflux), flux_ul(il^s,1:nwflux),flux_lc(il^s,1:nwflux)))
1153 
1154  flux_max(il^s,1:nwflux) = min(max(w(il^s,1:nwflux), w(ilp^s,1:nwflux), flux_md(il^s,1:nwflux)), &
1155  max(w(il^s,1:nwflux), flux_ul(il^s,1:nwflux),flux_lc(il^s,1:nwflux)))
1156  do iw=1,nwflux
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:nwflux,1) = u1_coeff(1) * w(ilpppp^s,1:nwflux) + u1_coeff(2) * w(ilppp^s,1:nwflux) + u1_coeff(3) * w(ilpp^s,1:nwflux) &
1174  + u1_coeff(4) * w(ilp^s,1:nwflux)
1175  f_array(il^s,1:nwflux,2) = u2_coeff(1) * w(ilppp^s,1:nwflux) + u2_coeff(2) * w(ilpp^s,1:nwflux) + u2_coeff(3) * w(ilp^s,1:nwflux) &
1176  + u2_coeff(4) * w(il^s,1:nwflux)
1177  f_array(il^s,1:nwflux,3) = u3_coeff(1) * w(ilpp^s,1:nwflux) + u3_coeff(2) * w(ilp^s,1:nwflux) + u3_coeff(3) * w(il^s,1:nwflux) &
1178  + u3_coeff(4) * w(ilm^s,1:nwflux)
1179  f_array(il^s,1:nwflux,4) = u4_coeff(1) * w(ilp^s,1:nwflux) + u4_coeff(2) * w(il^s,1:nwflux) + u4_coeff(3) * w(ilm^s,1:nwflux) &
1180  + u4_coeff(4) * w(ilmm^s,1:nwflux)
1181 
1182  beta(il^s,1:nwflux,1) = w(ilpppp^s,1:nwflux) * (547.d0 * w(ilpppp^s,1:nwflux) - 3882.d0 * w(ilppp^s,1:nwflux) + 4642.d0 * w(ilpp^s,1:nwflux) &
1183  - 1854.d0 * w(ilp^s,1:nwflux)) &
1184  + w(ilppp^s,1:nwflux) * (7043.d0 * w(ilppp^s,1:nwflux) - 17246.d0 * w(ilpp^s,1:nwflux) + 7042.d0 * w(ilp^s,1:nwflux)) &
1185  + w(ilpp^s,1:nwflux) * (11003.d0 * w(ilpp^s,1:nwflux) - 9402.d0 * w(ilp^s,1:nwflux)) + 2107.d0 * w(ilp^s,1:nwflux)**2
1186  beta(il^s,1:nwflux,2) = w(ilppp^s,1:nwflux) * (267.d0 * w(ilppp^s,1:nwflux) - 1642.d0 * w(ilpp^s,1:nwflux) + 1602.d0 * w(ilp^s,1:nwflux) &
1187  - 494.d0 * w(il^s,1:nwflux)) &
1188  + w(ilpp^s,1:nwflux) * (2843.d0 * w(ilpp^s,1:nwflux) - 5966.d0 * w(ilp^s,1:nwflux) + 1922.d0 * w(il^s,1:nwflux)) &
1189  + w(ilp^s,1:nwflux) * (3443.d0 * w(ilp^s,1:nwflux) - 2522.d0 * w(il^s,1:nwflux)) + 547.d0 * w(il^s,1:nwflux) ** 2
1190  beta(il^s,1:nwflux,3) = w(ilpp^s,1:nwflux) * (547.d0 * w(ilpp^s,1:nwflux) - 2522.d0 * w(ilp^s,1:nwflux) + 1922.d0 * w(il^s,1:nwflux) &
1191  - 494.d0 * w(ilm^s,1:nwflux)) &
1192  + w(ilp^s,1:nwflux) * (3443.d0 * w(ilp^s,1:nwflux) - 5966.d0 * w(il^s,1:nwflux) + 1602.d0 * w(ilm^s,1:nwflux)) &
1193  + w(il^s,1:nwflux) * (2843.d0 * w(il^s,1:nwflux) - 1642.d0 * w(ilm^s,1:nwflux)) + 267.d0 * w(ilm^s,1:nwflux) ** 2
1194  beta(il^s,1:nwflux,4) = w(ilp^s,1:nwflux) * (2107.d0 * w(ilp^s,1:nwflux) - 9402.d0 * w(il^s,1:nwflux) + 7042.d0 * w(ilm^s,1:nwflux) &
1195  - 1854.d0 * w(ilmm^s,1:nwflux)) &
1196  + w(il^s,1:nwflux) * (11003.d0 * w(il^s,1:nwflux) - 17246.d0 * w(ilm^s,1:nwflux) + 4642.d0 * w(ilmm^s,1:nwflux)) &
1197  + w(ilm^s,1:nwflux) * (7043.d0 * w(ilm^s,1:nwflux) - 3882.d0 * w(ilmm^s,1:nwflux)) + 547.d0 * w(ilmm^s,1:nwflux) ** 2
1198 
1199  alpha_sum(il^s,1:nwflux) = 0.0d0
1200  do i = 1,4
1201  alpha_array(il^s,1:nwflux,i) = d_array(i)/(beta(il^s,1:nwflux,i) + weno_eps_machine)**2
1202  alpha_sum(il^s,1:nwflux) = alpha_sum(il^s,1:nwflux) + alpha_array(il^s,1:nwflux,i)
1203  end do
1204  flux(il^s,1:nwflux) = 0.0d0
1205  do i = 1,4
1206  flux(il^s,1:nwflux) = flux(il^s,1:nwflux) + f_array(il^s,1:nwflux,i) * alpha_array(il^s,1:nwflux,i)/(alpha_sum(il^s,1:nwflux))
1207  end do
1208 
1209  select case(var)
1210  case(1)
1211  wrctmp(il^s,1:nwflux) = flux(il^s,1:nwflux)
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:nwflux) = w(ie^s,1:nwflux)-2.0d0*w(iep^s,1:nwflux)+w(iepp^s,1:nwflux)
1223 
1224  do iw = 1,nwflux
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:nwflux) = w(ilp^s,1:nwflux) + mpalpha * (w(ilp^s,1:nwflux) - w(ilpp^s,1:nwflux))
1236  flux_md(il^s,1:nwflux) = half * (w(il^s,1:nwflux) + w(ilp^s,1:nwflux) - dm4(il^s,1:nwflux))
1237  flux_lc(il^s,1:nwflux) = half * (3.d0 * w(ilp^s,1:nwflux) - w(ilpp^s,1:nwflux)) + mpbeta / 3.d0 * dm4(ilp^s,1:nwflux)
1238 
1239  flux_min(il^s,1:nwflux) = max(min(w(ilp^s,1:nwflux), w(il^s,1:nwflux), flux_md(il^s,1:nwflux)), &
1240  min(w(ilp^s,1:nwflux), flux_ul(il^s,1:nwflux),flux_lc(il^s,1:nwflux)))
1241 
1242  flux_max(il^s,1:nwflux) = min(max(w(ilp^s,1:nwflux), w(il^s,1:nwflux), flux_md(il^s,1:nwflux)), &
1243  max(w(ilp^s,1:nwflux), flux_ul(il^s,1:nwflux),flux_lc(il^s,1:nwflux)))
1244  do iw=1,nwflux
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
1286 end 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:76
mod_weno
Definition: mod_weno.t:1
mod_weno::fix_onelimiter
subroutine, public fix_onelimiter(ixIL, iLL, wCin, wCout)
Definition: mod_weno.t:42
mod_weno::weno5limiterl
subroutine, public weno5limiterl(ixIL, iLL, idims, w, wLC, var)
Definition: mod_weno.t:625
mod_weno::fix_onelimiter1
subroutine, public fix_onelimiter1(ixIL, iLL, wCin, wCout)
Definition: mod_weno.t:74
mod_weno::weno5nmlimiter
subroutine, public weno5nmlimiter(ixIL, iLL, idims, dxdim, w, wLC, wRC, var)
Definition: mod_weno.t:495
mod_weno::teno5adlimiter
subroutine, public teno5adlimiter(ixIL, iLL, idims, dxdim, w, wLC, wRC)
Definition: mod_weno.t:380
mod_weno::weno5limiterr
subroutine, public weno5limiterr(ixIL, iLL, idims, w, wRC, var)
Definition: mod_weno.t:709
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::weno5cu6limiter
subroutine, public weno5cu6limiter(ixIL, iLL, idims, w, wLC, wRC)
Definition: mod_weno.t:958
mod_weno::weno7limiter
subroutine, public weno7limiter(ixIL, iLL, idims, w, wLC, wRC, var)
Definition: mod_weno.t:1051
mod_weno::fix_limiter1
subroutine, public fix_limiter1(ixIL, iLL, wLCin, wRCin, wLCout, wRCout)
Definition: mod_weno.t:136
mod_weno::weno3limiter
subroutine, public weno3limiter(ixIL, iLL, idims, dxdim, w, wLC, wRC, var)
Definition: mod_weno.t:166
mod_weno::minmod
subroutine minmod(ixIL, ixOL, a, b, minm)
Definition: mod_weno.t:1258
mod_weno::weno5limiter
subroutine, public weno5limiter(ixIL, iLL, idims, dxdim, w, wLC, wRC, var)
Definition: mod_weno.t:255
mod_weno::weno5nmlimiterr
subroutine, public weno5nmlimiterr(ixIL, iLL, idims, w, wRC, var)
Definition: mod_weno.t:875
mod_weno::median
subroutine median(ixIL, ixOL, a, b, c, med)
Definition: mod_weno.t:1271
Generated by doxygen 1.9.1