MPI-AMRVAC  2.2
The MPI - Adaptive Mesh Refinement - Versatile Advection Code
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 by;
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 
15  implicit none
16  private
17 
18  public :: weno3limiter
19  public :: weno5limiter
20  public :: weno5nmlimiter
21  public :: weno5limiterl
22  public :: weno5nmlimiterl
23  public :: weno5limiterr
24  public :: weno5nmlimiterr
25  public :: weno7limiter
26  public :: exeno7limiter
27 
28 contains
29 
30  subroutine weno3limiter(ixI^L,iL^L,idims,dxdim,w,wLC,wRC,var)
32 
33  integer, intent(in) :: ixI^L, iL^L, idims, var
34  double precision, intent(in) :: dxdim
35  double precision, intent(in) :: w(ixi^s,1:nw)
36  double precision, intent(inout) :: wRC(ixi^s,1:nw),wLC(ixi^s,1:nw)
37  !> local
38  integer :: iLm^L, iLp^L, iLpp^L
39  double precision :: f_array(ixi^s,1:nw,2), d_array(2)
40  double precision :: beta(ixi^s,1:nw,2),tau(ixi^s,1:nw)
41  double precision :: u1_coeff(2), u2_coeff(2)
42  double precision :: alpha_array(ixi^s,1:nw,2), alpha_sum(ixi^s,1:nw), flux(ixi^s,1:nw)
43  integer :: i, iw
44  double precision, parameter :: weno_eps_machine = 1.0d-18
45 
46  ! iL^L holds the indices of interfaces to reconstruct to. Convention is that a center index holds the _right-side_ interface.
47  ilm^l=il^l-kr(idims,^d);
48  ilp^l=il^l+kr(idims,^d);
49  ilpp^l=ilp^l+kr(idims,^d);
50  d_array(1:2) = (/ 1.0d0/3.0d0, 2.0d0/3.0d0 /)
51  u1_coeff(1:2) = (/ -1.d0/2.d0, 3.d0/2.d0 /)
52  u2_coeff(1:2) = (/ 1.d0/2.d0, 1.d0/2.d0 /)
53 
54  !> left side
55  f_array(il^s,1:nwflux,1) = u1_coeff(1) * w(ilm^s,1:nwflux) + u1_coeff(2) * w(il^s,1:nwflux)
56  f_array(il^s,1:nwflux,2) = u2_coeff(1) * w(il^s,1:nwflux) + u2_coeff(2) * w(ilp^s,1:nwflux)
57 
58  beta(il^s,1:nwflux,1) = (w(il^s,1:nwflux) - w(ilm^s,1:nwflux))**2
59  beta(il^s,1:nwflux,2) = (w(ilp^s,1:nwflux) - w(il^s,1:nwflux))**2
60 
61  alpha_sum(il^s,1:nwflux) = 0.0d0
62  select case(var)
63  case(1)
64  do i = 1,2
65  alpha_array(il^s,1:nwflux,i) = d_array(i)/(beta(il^s,1:nwflux,i) + weno_eps_machine)**2
66  alpha_sum(il^s,1:nwflux) = alpha_sum(il^s,1:nwflux) + alpha_array(il^s,1:nwflux,i)
67  end do
68  case(2)
69  tau(il^s,1:nwflux) = abs(beta(il^s,1:nwflux,2) - beta(il^s,1:nwflux,1))
70  do i = 1,2
71  alpha_array(il^s,1:nwflux,i) = d_array(i) * (1.d0 + (tau(il^s,1:nwflux) / &
72  (beta(il^s,1:nwflux,i) + dxdim**2)))
73  alpha_sum(il^s,1:nwflux) = alpha_sum(il^s,1:nwflux) + alpha_array(il^s,1:nwflux,i)
74  end do
75  end select
76 
77  flux(il^s,1:nwflux) = 0.0d0
78  do i = 1,2
79  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))
80  end do
81 
82  !> left value at right interface
83  wlc(il^s,1:nwflux) = flux(il^s,1:nwflux)
84 
85  !> right side
86  f_array(il^s,1:nwflux,1) = u1_coeff(1) * w(ilpp^s,1:nwflux) + u1_coeff(2) * w(ilp^s,1:nwflux)
87  f_array(il^s,1:nwflux,2) = u2_coeff(1) * w(ilp^s,1:nwflux) + u2_coeff(2) * w(il^s,1:nwflux)
88 
89  beta(il^s,1:nwflux,1) = (w(ilpp^s,1:nwflux) - w(ilp^s,1:nwflux))**2
90  beta(il^s,1:nwflux,2) = (w(ilp^s,1:nwflux) - w(il^s,1:nwflux))**2
91 
92  alpha_sum(il^s,1:nwflux) = 0.0d0
93  select case(var)
94  case(1)
95  do i = 1,2
96  alpha_array(il^s,1:nwflux,i) = d_array(i)/(beta(il^s,1:nwflux,i) + weno_eps_machine)**2
97  alpha_sum(il^s,1:nwflux) = alpha_sum(il^s,1:nwflux) + alpha_array(il^s,1:nwflux,i)
98  end do
99  case(2)
100  tau(il^s,1:nwflux) = abs(beta(il^s,1:nwflux,2) - beta(il^s,1:nwflux,1))
101  do i = 1,2
102  alpha_array(il^s,1:nwflux,i) = d_array(i) * (1.d0 + (tau(il^s,1:nwflux) / &
103  (beta(il^s,1:nwflux,i) + dxdim**2)))
104  alpha_sum(il^s,1:nwflux) = alpha_sum(il^s,1:nwflux) + alpha_array(il^s,1:nwflux,i)
105  end do
106  end select
107 
108  flux(il^s,1:nwflux) = 0.0d0
109  do i = 1,2
110  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))
111  end do
112 
113  !> right value at right interface
114  wrc(il^s,1:nwflux) = flux(il^s,1:nwflux)
115 
116  end subroutine weno3limiter
117 
118  subroutine weno5limiter(ixI^L,iL^L,idims,dxdim,w,wLC,wRC,var)
120 
121  integer, intent(in) :: ixI^L, iL^L, idims, var
122  double precision, intent(in) :: dxdim
123  double precision, intent(in) :: w(ixi^s,1:nw)
124  double precision, intent(inout) :: wRC(ixi^s,1:nw),wLC(ixi^s,1:nw)
125  !> local
126  integer :: iLm^L, iLmm^L, iLp^L, iLpp^L, iLppp^L
127  double precision :: f_array(ixi^s,1:nw,3), d_array(3)
128  double precision :: beta(ixi^s,1:nw,3), beta_coeff(2)
129  double precision :: tau(ixi^s,1:nw), tmp(ixi^s,1:nw)
130  double precision :: u1_coeff(3), u2_coeff(3), u3_coeff(3)
131  double precision :: alpha_array(ixi^s,1:nw,3), alpha_sum(ixi^s,1:nw), flux(ixi^s,1:nw)
132  integer :: i
133  double precision, parameter :: weno_eps_machine = 1.0d-18
134  double precision :: lambda
135  double precision, parameter :: weno_dx_exp = 2.0d0/3.0d0
136 
137  ilm^l=il^l-kr(idims,^d);
138  ilmm^l=ilm^l-kr(idims,^d);
139  ilp^l=il^l+kr(idims,^d);
140  ilpp^l=ilp^l+kr(idims,^d);
141  ilppp^l=ilpp^l+kr(idims,^d);
142  lambda = dxdim**weno_dx_exp
143  beta_coeff(1:2) = (/ 1.0833333333333333d0, 0.25d0/)
144 ! reconstruction variation
145  d_array(1:3) = (/ 1.0d0/10.0d0, 3.0d0/5.0d0, 3.0d0/10.0d0 /)
146  u1_coeff(1:3) = (/ 1.d0/3.d0, -7.d0/6.d0, 11.d0/6.d0 /)
147  u2_coeff(1:3) = (/ -1.d0/6.d0, 5.d0/6.d0, 1.d0/3.d0 /)
148  u3_coeff(1:3) = (/ 1.d0/3.d0, 5.d0/6.d0, -1.d0/6.d0 /)
149 ! interpolation variation
150 ! d_array(1:3) = (/ 1.0d0/16.0d0, 10.0d0/16.0d0, 5.0d0/16.0d0 /)
151 ! u1_coeff(1:3) = (/ 3.d0/8.d0, -10.d0/8.d0, 15.d0/8.d0 /)
152 ! u2_coeff(1:3) = (/ -1.d0/8.d0, 6.d0/8.d0, 3.d0/8.d0 /)
153 ! u3_coeff(1:3) = (/ 3.d0/8.d0, 6.d0/8.d0, -1.d0/8.d0 /)
154 
155  !> left side
156  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)
157  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)
158  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)
159 
160  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 &
161  + beta_coeff(2) * (w(ilmm^s,1:nwflux) - 4.0d0 * w(ilm^s,1:nwflux) + 3.0d0*w(il^s,1:nwflux))**2
162  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 &
163  + beta_coeff(2) * (w(ilm^s,1:nwflux) - w(ilp^s,1:nwflux))**2
164  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 &
165  + beta_coeff(2) * (3.0d0 * w(il^s, 1:nwflux) - 4.0d0 * w(ilp^s,1:nwflux) + w(ilpp^s,1:nwflux))**2
166 
167  alpha_sum(il^s,1:nwflux) = 0.0d0
168  select case(var)
169  ! case1 for wenojs, case2 for wenoz, case3 for wenoz+
170  case(1)
171  do i = 1,3
172  alpha_array(il^s,1:nwflux,i) = d_array(i)/(beta(il^s,1:nwflux,i) + weno_eps_machine)**2
173  alpha_sum(il^s,1:nwflux) = alpha_sum(il^s,1:nwflux) + alpha_array(il^s,1:nwflux,i)
174  end do
175  case(2)
176  tau(il^s,1:nwflux) = abs(beta(il^s,1:nwflux,1) - beta(il^s,1:nwflux,3))
177  do i = 1,3
178  alpha_array(il^s,1:nwflux,i) = d_array(i) * (1.d0 + (tau(il^s,1:nwflux) / &
179  (beta(il^s,1:nwflux,i) + weno_eps_machine))**2)
180  alpha_sum(il^s,1:nwflux) = alpha_sum(il^s,1:nwflux) + alpha_array(il^s,1:nwflux,i)
181  end do
182  case(3)
183  tau(il^s,1:nwflux) = abs(beta(il^s,1:nwflux,1) - beta(il^s,1:nwflux,3))
184  do i = 1,3
185  tmp(il^s,1:nwflux) = (tau(il^s,1:nwflux) + weno_eps_machine) / (beta(il^s,1:nwflux,i) + weno_eps_machine)
186  alpha_array(il^s,1:nwflux,i) = d_array(i) * (1.0d0 + tmp(il^s,1:nwflux)**2 + lambda/tmp(il^s,1:nwflux))
187  alpha_sum(il^s,1:nwflux) = alpha_sum(il^s,1:nwflux) + alpha_array(il^s,1:nwflux,i)
188  end do
189  end select
190 
191  flux(il^s,1:nwflux) = 0.0d0
192  do i = 1,3
193  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))
194  end do
195 
196  !> left value at right interface
197  wlc(il^s,1:nwflux) = flux(il^s,1:nwflux)
198 
199  !> right side
200  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)
201  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)
202  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)
203 
204  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 &
205  + beta_coeff(2) * (w(ilppp^s,1:nwflux) - 4.0d0 * w(ilpp^s,1:nwflux) + 3.0d0*w(ilp^s,1:nwflux))**2
206  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 &
207  + beta_coeff(2) * (w(ilpp^s,1:nwflux) - w(il^s,1:nwflux))**2
208  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 &
209  + beta_coeff(2) * (3.0d0 * w(ilp^s, 1:nwflux) - 4.0d0 * w(il^s,1:nwflux) + w(ilm^s,1:nwflux))**2
210 
211  alpha_sum(il^s,1:nwflux) = 0.0d0
212  select case(var)
213  case(1)
214  do i = 1,3
215  alpha_array(il^s,1:nwflux,i) = d_array(i)/(beta(il^s,1:nwflux,i) + weno_eps_machine)**2
216  alpha_sum(il^s,1:nwflux) = alpha_sum(il^s,1:nwflux) + alpha_array(il^s,1:nwflux,i)
217  end do
218  case(2)
219  tau(il^s,1:nwflux) = abs(beta(il^s,1:nwflux,1) - beta(il^s,1:nwflux,3))
220  do i = 1,3
221  alpha_array(il^s,1:nwflux,i) = d_array(i) * (1.d0 + (tau(il^s,1:nwflux) / &
222  (beta(il^s,1:nwflux,i) + weno_eps_machine))**2)
223  alpha_sum(il^s,1:nwflux) = alpha_sum(il^s,1:nwflux) + alpha_array(il^s,1:nwflux,i)
224  end do
225  case(3)
226  tau(il^s,1:nwflux) = abs(beta(il^s,1:nwflux,1) - beta(il^s,1:nwflux,3))
227  do i = 1,3
228  tmp(il^s,1:nwflux) = (tau(il^s,1:nwflux) + weno_eps_machine) / (beta(il^s,1:nwflux,i) + weno_eps_machine)
229  alpha_array(il^s,1:nwflux,i) = d_array(i) * (1.0d0 + tmp(il^s,1:nwflux)**2 + lambda/tmp(il^s,1:nwflux))
230  alpha_sum(il^s,1:nwflux) = alpha_sum(il^s,1:nwflux) + alpha_array(il^s,1:nwflux,i)
231  end do
232  end select
233  flux(il^s,1:nwflux) = 0.0d0
234  do i = 1,3
235  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))
236  end do
237 
238  !> right value at right interface
239  wrc(il^s,1:nwflux) = flux(il^s,1:nwflux)
240 
241  end subroutine weno5limiter
242 
243  subroutine weno5nmlimiter(ixI^L,iL^L,idims,dxdim,w,wLC,wRC,var)
245 
246  integer, intent(in) :: ixI^L,iL^L,idims,var
247  double precision, intent(in) :: dxdim
248  double precision, intent(in) :: w(ixi^s,1:nw)
249  double precision, intent(inout) :: wRC(ixi^s,1:nw),wLC(ixi^s,1:nw)
250  !> local
251  integer :: iLm^L, iLmm^L, iLp^L, iLpp^L, iLppp^L
252  double precision :: f_array(ixi^s,1:nw,3), d_array(3)
253  double precision :: beta(ixi^s,1:nw,3), beta_coeff(2)
254  double precision :: tau(ixi^s,1:nw), tmp(ixi^s,1:nw)
255  double precision :: u1_coeff(3), u2_coeff(3), u3_coeff(3)
256  double precision :: alpha_array(ixi^s,1:nw,3), alpha_sum(ixi^s,1:nw), flux(ixi^s,1:nw)
257  double precision :: wc(ixi^s,1:nw), wd(ixi^s,1:nw), we(ixi^s,1:nw)
258  integer :: i,j
259  double precision, parameter :: weno_eps_machine = 1.0d-18
260  double precision, parameter :: weno_dx_exp = 2.0d0/3.0d0
261  double precision :: lambda(ixi^s)
262 
263  ilm^l=il^l-kr(idims,^d);
264  ilmm^l=ilm^l-kr(idims,^d);
265  ilp^l=il^l+kr(idims,^d);
266  ilpp^l=ilp^l+kr(idims,^d);
267  ilppp^l=ilpp^l+kr(idims,^d);
268  lambda=block%dx(il^s,idims)**weno_dx_exp
269  beta_coeff(1:2) = (/ 1.0833333333333333d0, 0.25d0/)
270  d_array(1:3) = (/ 1.0d0/10.0d0, 3.0d0/5.0d0, 3.0d0/10.0d0 /)
271  u1_coeff(1:3) = (/ -2.d0/3.d0, -1.d0/3.d0, 2.d0 /)
272  u2_coeff(1:3) = (/ -1.d0/3.d0, 2.d0/3.d0, 2.d0/3.d0 /)
273  u3_coeff(1:3) = (/ 2.d0/3.d0, 2.d0/3.d0, -1.d0/3.d0 /)
274  do i = 1, nwflux
275  wc(il^s,i) = (block%dx(ilp^s,idims) * w(il^s,i) + block%dx(il^s,idims) * w(ilp^s,i)) / &
276  (block%dx(ilp^s,idims) + block%dx(il^s,idims))
277  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)) / &
278  (block%dx(ilmm^s,idims) + block%dx(ilm^s,idims))
279  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)) / &
280  (block%dx(ilppp^s,idims) + block%dx(ilpp^s,idims))
281  enddo
282  !> left side
283  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)
284  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)
285  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)
286 
287  beta(il^s,1:nwflux,1) = beta_coeff(1) * (wc(ilm^s,1:nwflux) - wd(il^s,1:nwflux))**2 &
288  + beta_coeff(2) * (2.d0 * w(il^s,1:nwflux) - wc(ilm^s,1:nwflux) - wd(il^s,1:nwflux))**2
289  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 &
290  + beta_coeff(2) * (wc(ilm^s,1:nwflux) - wc(il^s,1:nwflux))**2
291  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 &
292  + beta_coeff(2) * (3.0d0 * wc(il^s, 1:nwflux) - 4.0d0 * w(ilp^s,1:nwflux) + wc(ilp^s,1:nwflux))**2
293  alpha_sum(il^s,1:nwflux) = 0.0d0
294  select case(var)
295  ! case1 for wenojs, case2 for wenoz, case3 for wenoz+
296  case(1)
297  do i = 1,3
298  alpha_array(il^s,1:nwflux,i) = d_array(i)/(4.d0 * beta(il^s,1:nwflux,i) + weno_eps_machine)**2
299  alpha_sum(il^s,1:nwflux) = alpha_sum(il^s,1:nwflux) + alpha_array(il^s,1:nwflux,i)
300  end do
301  case(2)
302  tau(il^s,1:nwflux) = abs(beta(il^s,1:nwflux,1) - beta(il^s,1:nwflux,3)) * 4.d0
303  do i = 1,3
304  alpha_array(il^s,1:nwflux,i) = d_array(i) * (1.d0 + (tau(il^s,1:nwflux) / &
305  (4.d0 * beta(il^s,1:nwflux,i) + weno_eps_machine))**2)
306  alpha_sum(il^s,1:nwflux) = alpha_sum(il^s,1:nwflux) + alpha_array(il^s,1:nwflux,i)
307  end do
308  case(3)
309  tau(il^s,1:nwflux) = abs(beta(il^s,1:nwflux,1) - beta(il^s,1:nwflux,3)) * 4.d0
310  do i=1,3
311  do j=1,nwflux
312  tmp(il^s,j) = (tau(il^s,j) + weno_eps_machine) / (4.d0 * beta(il^s,j,i) + weno_eps_machine)
313  alpha_array(il^s,j,i) = d_array(i) * (1.0d0 + tmp(il^s,j)**2 + lambda(il^s)/tmp(il^s,j))
314  alpha_sum(il^s,j) = alpha_sum(il^s,j) + alpha_array(il^s,j,i)
315  end do
316  end do
317  end select
318  flux(il^s,1:nwflux) = 0.0d0
319  do i = 1,3
320  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))
321  end do
322  !> left value at right interface
323  wlc(il^s,1:nwflux) = flux(il^s,1:nwflux)
324  !> right side
325  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)
326  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)
327  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)
328  beta(il^s,1:nwflux,1) = beta_coeff(1) * (wc(ilp^s,1:nwflux) - we(il^s,1:nwflux))**2 &
329  + beta_coeff(2) * (2.d0 * w(ilp^s,1:nwflux) - wc(ilp^s,1:nwflux) - we(il^s,1:nwflux))**2
330  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 &
331  + beta_coeff(2) * (wc(ilp^s,1:nwflux) - wc(il^s,1:nwflux))**2
332  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 &
333  + beta_coeff(2) * (3.0d0 * wc(il^s, 1:nwflux) - 4.0d0 * w(il^s,1:nwflux) + wc(ilm^s,1:nwflux))**2
334  alpha_sum(il^s,1:nwflux) = 0.0d0
335  select case(var)
336  case(1)
337  do i = 1,3
338  alpha_array(il^s,1:nwflux,i) = d_array(i)/(4.d0 * beta(il^s,1:nwflux,i) + weno_eps_machine)**2
339  alpha_sum(il^s,1:nwflux) = alpha_sum(il^s,1:nwflux) + alpha_array(il^s,1:nwflux,i)
340  end do
341  case(2)
342  tau(il^s,1:nwflux) = abs(beta(il^s,1:nwflux,1) - beta(il^s,1:nwflux,3)) * 4.d0
343  do i = 1,3
344  alpha_array(il^s,1:nwflux,i) = d_array(i) * (1.d0 + (tau(il^s,1:nwflux) / &
345  (4.d0 * beta(il^s,1:nwflux,i) + weno_eps_machine))**2)
346  alpha_sum(il^s,1:nwflux) = alpha_sum(il^s,1:nwflux) + alpha_array(il^s,1:nwflux,i)
347  end do
348  case(3)
349  tau(il^s,1:nwflux) = abs(beta(il^s,1:nwflux,1) - beta(il^s,1:nwflux,3)) * 4.d0
350  do i = 1,3
351  do j = 1,nwflux
352  tmp(il^s,j) = (tau(il^s,j) + weno_eps_machine) / (4.d0 * beta(il^s,j,i) + weno_eps_machine)
353  alpha_array(il^s,j,i) = d_array(i) * (1.0d0 + tmp(il^s,j)**2 + lambda(il^s)/tmp(il^s,j))
354  alpha_sum(il^s,j) = alpha_sum(il^s,j) + alpha_array(il^s,j,i)
355  end do
356  end do
357  end select
358  flux(il^s,1:nwflux) = 0.0d0
359  do i = 1,3
360  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))
361  end do
362  !> right value at right interface
363  wrc(il^s,1:nwflux) = flux(il^s,1:nwflux)
364 
365  end subroutine weno5nmlimiter
366 
367  subroutine weno5limiterl(ixI^L,iL^L,idims,w,wLC,var)
369 
370  integer, intent(in) :: ixI^L, iL^L, idims
371  integer, intent(in) :: var
372  double precision, intent(in) :: w(ixi^s,1:nw)
373  double precision, intent(inout) :: wLC(ixi^s,1:nw)
374  !> local
375  integer :: iLm^L, iLmm^L, iLp^L, iLpp^L, iLppp^L
376  double precision :: f_array(ixi^s,1:nw,3), d_array(3)
377  double precision :: beta(ixi^s,1:nw,3), beta_coeff(2)
378  double precision :: tau(ixi^s,1:nw), tmp(ixi^s,1:nw)
379  double precision :: u1_coeff(3), u2_coeff(3), u3_coeff(3)
380  double precision :: alpha_array(ixi^s,1:nw,3), alpha_sum(ixi^s,1:nw), flux(ixi^s,1:nw)
381  integer :: i,j
382  double precision :: lambda(ixi^s)
383  double precision, parameter :: weno_dx_exp = 2.0d0/3.0d0
384  double precision, parameter :: weno_eps_machine = 1.0d-18
385 
386  ilm^l=il^l-kr(idims,^d);
387  ilmm^l=ilm^l-kr(idims,^d);
388  ilp^l=il^l+kr(idims,^d);
389  ilpp^l=ilp^l+kr(idims,^d);
390  lambda=block%dx(il^s,idims)**weno_dx_exp
391  beta_coeff(1:2) = (/ 1.0833333333333333d0, 0.25d0/)
392 ! reconstruction variation
393  d_array(1:3) = (/ 1.0d0/10.0d0, 3.0d0/5.0d0, 3.0d0/10.0d0 /)
394  u1_coeff(1:3) = (/ 1.d0/3.d0, -7.d0/6.d0, 11.d0/6.d0 /)
395  u2_coeff(1:3) = (/ -1.d0/6.d0, 5.d0/6.d0, 1.d0/3.d0 /)
396  u3_coeff(1:3) = (/ 1.d0/3.d0, 5.d0/6.d0, -1.d0/6.d0 /)
397 ! interpolation variation
398 ! d_array(1:3) = (/ 1.0d0/16.0d0, 10.0d0/16.0d0, 5.0d0/16.0d0 /)
399 ! u1_coeff(1:3) = (/ 3.d0/8.d0, -10.d0/8.d0, 15.d0/8.d0 /)
400 ! u2_coeff(1:3) = (/ -1.d0/8.d0, 6.d0/8.d0, 3.d0/8.d0 /)
401 ! u3_coeff(1:3) = (/ 3.d0/8.d0, 6.d0/8.d0, -1.d0/8.d0 /)
402 
403  !> left side
404  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)
405  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)
406  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)
407 
408  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 &
409  + beta_coeff(2) * (w(ilmm^s,1:nwflux) - 4.0d0 * w(ilm^s,1:nwflux) + 3.0d0*w(il^s,1:nwflux))**2
410  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 &
411  + beta_coeff(2) * (w(ilm^s,1:nwflux) - w(ilp^s,1:nwflux))**2
412  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 &
413  + beta_coeff(2) * (3.0d0 * w(il^s, 1:nwflux) - 4.0d0 * w(ilp^s,1:nwflux) + w(ilpp^s,1:nwflux))**2
414 
415  alpha_sum(il^s,1:nwflux) = 0.0d0
416  select case(var)
417  ! case1 for wenojs, case2 for wenoz
418  case(1)
419  do i = 1,3
420  alpha_array(il^s,1:nwflux,i) = d_array(i)/(beta(il^s,1:nwflux,i) + weno_eps_machine)**2
421  alpha_sum(il^s,1:nwflux) = alpha_sum(il^s,1:nwflux) + alpha_array(il^s,1:nwflux,i)
422  end do
423  case(2)
424  tau(il^s,1:nwflux) = abs(beta(il^s,1:nwflux,1) - beta(il^s,1:nwflux,3))
425  do i = 1,3
426  alpha_array(il^s,1:nwflux,i) = d_array(i) * (1.d0 + (tau(il^s,1:nwflux) / &
427  (beta(il^s,1:nwflux,i) + weno_eps_machine))**2)
428  alpha_sum(il^s,1:nwflux) = alpha_sum(il^s,1:nwflux) + alpha_array(il^s,1:nwflux,i)
429  end do
430  case(3)
431  tau(il^s,1:nwflux) = abs(beta(il^s,1:nwflux,1) - beta(il^s,1:nwflux,3)) * 4.d0
432  do i=1,3
433  do j=1,nwflux
434  tmp(il^s,j) = (tau(il^s,j) + weno_eps_machine) / (4.d0 * beta(il^s,j,i) + weno_eps_machine)
435  alpha_array(il^s,j,i) = d_array(i) * (1.0d0 + tmp(il^s,j)**2 + lambda(il^s)/tmp(il^s,j))
436  alpha_sum(il^s,j) = alpha_sum(il^s,j) + alpha_array(il^s,j,i)
437  end do
438  end do
439  end select
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) * alpha_array(il^s,1:nwflux,i)/(alpha_sum(il^s,1:nwflux))
443  end do
444 
445  !> left value at right interface
446  wlc(il^s,1:nwflux) = flux(il^s,1:nwflux)
447 
448  end subroutine weno5limiterl
449 
450  subroutine weno5limiterr(ixI^L,iL^L,idims,w,wRC,var)
452 
453  integer, intent(in) :: ixI^L, iL^L, idims
454  integer, intent(in) :: var
455  double precision, intent(in) :: w(ixi^s,1:nw)
456  double precision, intent(inout) :: wRC(ixi^s,1:nw)
457  !> local
458  integer :: iLm^L, iLmm^L, iLp^L, iLpp^L, iLppp^L
459  double precision :: f_array(ixi^s,1:nw,3), d_array(3)
460  double precision :: beta(ixi^s,1:nw,3), beta_coeff(2)
461  double precision :: tau(ixi^s,1:nw), tmp(ixi^s,1:nw)
462  double precision :: u1_coeff(3), u2_coeff(3), u3_coeff(3)
463  double precision :: alpha_array(ixi^s,1:nw,3), alpha_sum(ixi^s,1:nw), flux(ixi^s,1:nw)
464  integer :: i,j
465  double precision :: lambda(ixi^s)
466  double precision, parameter :: weno_dx_exp = 2.0d0/3.0d0
467  double precision, parameter :: weno_eps_machine = 1.0d-18
468 
469  ilm^l=il^l-kr(idims,^d);
470  ilp^l=il^l+kr(idims,^d);
471  ilpp^l=ilp^l+kr(idims,^d);
472  ilppp^l=ilpp^l+kr(idims,^d);
473  lambda=block%dx(il^s,idims)**weno_dx_exp
474  beta_coeff(1:2) = (/ 1.0833333333333333d0, 0.25d0/)
475 ! reconstruction variation
476  d_array(1:3) = (/ 1.0d0/10.0d0, 3.0d0/5.0d0, 3.0d0/10.0d0 /)
477  u1_coeff(1:3) = (/ 1.d0/3.d0, -7.d0/6.d0, 11.d0/6.d0 /)
478  u2_coeff(1:3) = (/ -1.d0/6.d0, 5.d0/6.d0, 1.d0/3.d0 /)
479  u3_coeff(1:3) = (/ 1.d0/3.d0, 5.d0/6.d0, -1.d0/6.d0 /)
480 ! interpolation variation
481 ! d_array(1:3) = (/ 1.0d0/16.0d0, 10.0d0/16.0d0, 5.0d0/16.0d0 /)
482 ! u1_coeff(1:3) = (/ 3.d0/8.d0, -10.d0/8.d0, 15.d0/8.d0 /)
483 ! u2_coeff(1:3) = (/ -1.d0/8.d0, 6.d0/8.d0, 3.d0/8.d0 /)
484 ! u3_coeff(1:3) = (/ 3.d0/8.d0, 6.d0/8.d0, -1.d0/8.d0 /)
485 
486  !> right side
487  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)
488  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)
489  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)
490 
491  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 &
492  + beta_coeff(2) * (w(ilppp^s,1:nwflux) - 4.0d0 * w(ilpp^s,1:nwflux) + 3.0d0*w(ilp^s,1:nwflux))**2
493  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 &
494  + beta_coeff(2) * (w(ilpp^s,1:nwflux) - w(il^s,1:nwflux))**2
495  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 &
496  + beta_coeff(2) * (3.0d0 * w(ilp^s, 1:nwflux) - 4.0d0 * w(il^s,1:nwflux) + w(ilm^s,1:nwflux))**2
497 
498  alpha_sum(il^s,1:nwflux) = 0.0d0
499  select case(var)
500  case(1)
501  do i = 1,3
502  alpha_array(il^s,1:nwflux,i) = d_array(i)/(beta(il^s,1:nwflux,i) + weno_eps_machine)**2
503  alpha_sum(il^s,1:nwflux) = alpha_sum(il^s,1:nwflux) + alpha_array(il^s,1:nwflux,i)
504  end do
505  case(2)
506  tau(il^s,1:nwflux) = abs(beta(il^s,1:nwflux,1) - beta(il^s,1:nwflux,3))
507  do i = 1,3
508  alpha_array(il^s,1:nwflux,i) = d_array(i) * (1.d0 + (tau(il^s,1:nwflux) / &
509  (beta(il^s,1:nwflux,i) + weno_eps_machine))**2)
510  alpha_sum(il^s,1:nwflux) = alpha_sum(il^s,1:nwflux) + alpha_array(il^s,1:nwflux,i)
511  end do
512  case(3)
513  tau(il^s,1:nwflux) = abs(beta(il^s,1:nwflux,1) - beta(il^s,1:nwflux,3)) * 4.d0
514  do i = 1,3
515  do j = 1,nwflux
516  tmp(il^s,j) = (tau(il^s,j) + weno_eps_machine) / (4.d0 * beta(il^s,j,i) + weno_eps_machine)
517  alpha_array(il^s,j,i) = d_array(i) * (1.0d0 + tmp(il^s,j)**2 + lambda(il^s)/tmp(il^s,j))
518  alpha_sum(il^s,j) = alpha_sum(il^s,j) + alpha_array(il^s,j,i)
519  end do
520  end do
521  end select
522  flux(il^s,1:nwflux) = 0.0d0
523  do i = 1,3
524  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))
525  end do
526 
527  !> right value at right interface
528  wrc(il^s,1:nwflux) = flux(il^s,1:nwflux)
529 
530  end subroutine weno5limiterr
531 
532  subroutine weno5nmlimiterl(ixI^L,iL^L,idims,w,wLC,var)
534 
535  integer, intent(in) :: ixI^L,iL^L,idims,var
536  double precision, intent(in) :: w(ixi^s,1:nw)
537  double precision, intent(inout) :: wLC(ixi^s,1:nw)
538  !> local
539  integer :: iLm^L, iLmm^L, iLp^L, iLpp^L
540  double precision :: f_array(ixi^s,1:nw,3), d_array(3)
541  double precision :: beta(ixi^s,1:nw,3), beta_coeff(2)
542  double precision :: tau(ixi^s,1:nw), tmp(ixi^s,1:nw)
543  double precision :: u1_coeff(3), u2_coeff(3), u3_coeff(3)
544  double precision :: alpha_array(ixi^s,1:nw,3), alpha_sum(ixi^s,1:nw), flux(ixi^s,1:nw)
545  double precision :: wc(ixi^s,1:nw), wd(ixi^s,1:nw)
546  integer :: i,j
547  double precision, parameter :: weno_eps_machine = 1.0d-18
548  double precision, parameter :: weno_dx_exp = 2.0d0/3.0d0
549  double precision :: lambda(ixi^s)
550 
551  ilm^l=il^l-kr(idims,^d);
552  ilmm^l=ilm^l-kr(idims,^d);
553  ilp^l=il^l+kr(idims,^d);
554  ilpp^l=ilp^l+kr(idims,^d);
555  lambda=block%dx(il^s,idims)**weno_dx_exp
556  beta_coeff(1:2) = (/ 1.0833333333333333d0, 0.25d0/)
557  d_array(1:3) = (/ 1.0d0/10.0d0, 3.0d0/5.0d0, 3.0d0/10.0d0 /)
558  u1_coeff(1:3) = (/ -2.d0/3.d0, -1.d0/3.d0, 2.d0 /)
559  u2_coeff(1:3) = (/ -1.d0/3.d0, 2.d0/3.d0, 2.d0/3.d0 /)
560  u3_coeff(1:3) = (/ 2.d0/3.d0, 2.d0/3.d0, -1.d0/3.d0 /)
561  do i = 1, nwflux
562  wc(il^s,i) = (block%dx(ilp^s,idims) * w(il^s,i) + block%dx(il^s,idims) * w(ilp^s,i)) / &
563  (block%dx(ilp^s,idims) + block%dx(il^s,idims))
564  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)) / &
565  (block%dx(ilmm^s,idims) + block%dx(ilm^s,idims))
566  enddo
567  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)
568  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)
569  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)
570  beta(il^s,1:nwflux,1) = beta_coeff(1) * (wc(ilm^s,1:nwflux) - wd(il^s,1:nwflux))**2 &
571  + beta_coeff(2) * (2.d0 * w(il^s,1:nwflux) - wc(ilm^s,1:nwflux) - wd(il^s,1:nwflux))**2
572  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 &
573  + beta_coeff(2) * (wc(ilm^s,1:nwflux) - wc(il^s,1:nwflux))**2
574  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 &
575  + beta_coeff(2) * (3.0d0 * wc(il^s, 1:nwflux) - 4.0d0 * w(ilp^s,1:nwflux) + wc(ilp^s,1:nwflux))**2
576  alpha_sum(il^s,1:nwflux) = 0.0d0
577  select case(var)
578  ! case1 for wenojs, case2 for wenoz, case3 for wenoz+
579  case(1)
580  do i = 1,3
581  alpha_array(il^s,1:nwflux,i) = d_array(i)/(4.d0 * beta(il^s,1:nwflux,i) + weno_eps_machine)**2
582  alpha_sum(il^s,1:nwflux) = alpha_sum(il^s,1:nwflux) + alpha_array(il^s,1:nwflux,i)
583  end do
584  case(2)
585  tau(il^s,1:nwflux) = abs(beta(il^s,1:nwflux,1) - beta(il^s,1:nwflux,3)) * 4.d0
586  do i = 1,3
587  alpha_array(il^s,1:nwflux,i) = d_array(i) * (1.d0 + (tau(il^s,1:nwflux) / &
588  (4.d0 * beta(il^s,1:nwflux,i) + weno_eps_machine))**2)
589  alpha_sum(il^s,1:nwflux) = alpha_sum(il^s,1:nwflux) + alpha_array(il^s,1:nwflux,i)
590  end do
591  case(3)
592  tau(il^s,1:nwflux) = abs(beta(il^s,1:nwflux,1) - beta(il^s,1:nwflux,3)) * 4.d0
593  do i=1,3
594  do j=1,nwflux
595  tmp(il^s,j) = (tau(il^s,j) + weno_eps_machine) / (4.d0 * beta(il^s,j,i) + weno_eps_machine)
596  alpha_array(il^s,j,i) = d_array(i) * (1.0d0 + tmp(il^s,j)**2 + lambda(il^s)/tmp(il^s,j))
597  alpha_sum(il^s,j) = alpha_sum(il^s,j) + alpha_array(il^s,j,i)
598  end do
599  end do
600  end select
601  flux(il^s,1:nwflux) = 0.0d0
602  do i = 1,3
603  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))
604  end do
605  !> left value at right interface
606  wlc(il^s,1:nwflux) = flux(il^s,1:nwflux)
607 
608  end subroutine weno5nmlimiterl
609 
610  subroutine weno5nmlimiterr(ixI^L,iL^L,idims,w,wRC,var)
612 
613  integer, intent(in) :: ixI^L,iL^L,idims,var
614  double precision, intent(in) :: w(ixi^s,1:nw)
615  double precision, intent(inout) :: wRC(ixi^s,1:nw)
616  !> local
617  integer :: iLm^L,iLp^L,iLpp^L,iLppp^L
618  double precision :: f_array(ixi^s,1:nw,3), d_array(3)
619  double precision :: beta(ixi^s,1:nw,3), beta_coeff(2)
620  double precision :: tau(ixi^s,1:nw), tmp(ixi^s,1:nw)
621  double precision :: u1_coeff(3), u2_coeff(3), u3_coeff(3)
622  double precision :: alpha_array(ixi^s,1:nw,3), alpha_sum(ixi^s,1:nw), flux(ixi^s,1:nw)
623  double precision :: wc(ixi^s,1:nw), we(ixi^s,1:nw)
624  integer :: i,j
625  double precision, parameter :: weno_eps_machine = 1.0d-18
626  double precision, parameter :: weno_dx_exp = 2.0d0/3.0d0
627  double precision :: lambda(ixi^s)
628 
629  ilm^l=il^l-kr(idims,^d);
630  ilp^l=il^l+kr(idims,^d);
631  ilpp^l=ilp^l+kr(idims,^d);
632  ilppp^l=ilpp^l+kr(idims,^d);
633  lambda=block%dx(il^s,idims)**weno_dx_exp
634  beta_coeff(1:2) = (/ 1.0833333333333333d0, 0.25d0/)
635  d_array(1:3) = (/ 1.0d0/10.0d0, 3.0d0/5.0d0, 3.0d0/10.0d0 /)
636  u1_coeff(1:3) = (/ -2.d0/3.d0, -1.d0/3.d0, 2.d0 /)
637  u2_coeff(1:3) = (/ -1.d0/3.d0, 2.d0/3.d0, 2.d0/3.d0 /)
638  u3_coeff(1:3) = (/ 2.d0/3.d0, 2.d0/3.d0, -1.d0/3.d0 /)
639  do i = 1, nwflux
640  wc(il^s,i) = (block%dx(ilp^s,idims) * w(il^s,i) + block%dx(il^s,idims) * w(ilp^s,i)) / &
641  (block%dx(ilp^s,idims) + block%dx(il^s,idims))
642  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)) / &
643  (block%dx(ilppp^s,idims) + block%dx(ilpp^s,idims))
644  enddo
645  !> right side
646  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)
647  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)
648  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)
649  beta(il^s,1:nwflux,1) = beta_coeff(1) * (wc(ilp^s,1:nwflux) - we(il^s,1:nwflux))**2 &
650  + beta_coeff(2) * (2.d0 * w(ilp^s,1:nwflux) - wc(ilp^s,1:nwflux) - we(il^s,1:nwflux))**2
651  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 &
652  + beta_coeff(2) * (wc(ilp^s,1:nwflux) - wc(il^s,1:nwflux))**2
653  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 &
654  + beta_coeff(2) * (3.0d0 * wc(il^s, 1:nwflux) - 4.0d0 * w(il^s,1:nwflux) + wc(ilm^s,1:nwflux))**2
655  alpha_sum(il^s,1:nwflux) = 0.0d0
656  select case(var)
657  case(1)
658  do i = 1,3
659  alpha_array(il^s,1:nwflux,i) = d_array(i)/(4.d0 * beta(il^s,1:nwflux,i) + weno_eps_machine)**2
660  alpha_sum(il^s,1:nwflux) = alpha_sum(il^s,1:nwflux) + alpha_array(il^s,1:nwflux,i)
661  end do
662  case(2)
663  tau(il^s,1:nwflux) = abs(beta(il^s,1:nwflux,1) - beta(il^s,1:nwflux,3)) * 4.d0
664  do i = 1,3
665  alpha_array(il^s,1:nwflux,i) = d_array(i) * (1.d0 + (tau(il^s,1:nwflux) / &
666  (4.d0 * beta(il^s,1:nwflux,i) + weno_eps_machine))**2)
667  alpha_sum(il^s,1:nwflux) = alpha_sum(il^s,1:nwflux) + alpha_array(il^s,1:nwflux,i)
668  end do
669  case(3)
670  tau(il^s,1:nwflux) = abs(beta(il^s,1:nwflux,1) - beta(il^s,1:nwflux,3)) * 4.d0
671  do i = 1,3
672  do j = 1,nwflux
673  tmp(il^s,j) = (tau(il^s,j) + weno_eps_machine) / (4.d0 * beta(il^s,j,i) + weno_eps_machine)
674  alpha_array(il^s,j,i) = d_array(i) * (1.0d0 + tmp(il^s,j)**2 + lambda(il^s)/tmp(il^s,j))
675  alpha_sum(il^s,j) = alpha_sum(il^s,j) + alpha_array(il^s,j,i)
676  end do
677  end do
678  end select
679  flux(il^s,1:nwflux) = 0.0d0
680  do i = 1,3
681  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))
682  end do
683  !> right value at right interface
684  wrc(il^s,1:nwflux) = flux(il^s,1:nwflux)
685 
686  end subroutine weno5nmlimiterr
687 
688  subroutine weno7limiter(ixI^L,iL^L,idims,w,wLC,wRC,var)
690 
691  integer, intent(in) :: ixI^L, iL^L, idims, var
692  double precision, intent(in) :: w(ixi^s,1:nw)
693  double precision, intent(inout) :: wRC(ixi^s,1:nw),wLC(ixi^s,1:nw)
694  !> local
695  integer :: iLm^L, iLmm^L, iLmmm^L
696  integer :: iLp^L, iLpp^L, iLppp^L, iLpppp^L
697  integer :: id^L, idp^L, idpp^L, idm^L, ie^L, iem^L, iep^L, iepp^L
698 
699  double precision, dimension(4) :: d_array, u1_coeff, u2_coeff, u3_coeff, u4_coeff
700  double precision, dimension(ixI^S,1:nw,4) :: f_array, beta, alpha_array
701  double precision, dimension(ixI^S) :: a, b, c, tmp, tmp2, tmp3
702  double precision, dimension(ixI^S,1:nw) :: alpha_sum, d, dm4
703  double precision, dimension(ixI^S,1:nw) :: flux, flux_min, flux_max, flux_ul, flux_md, flux_lc
704  integer :: i,iw
705  double precision, parameter :: mpalpha = 2.d0, mpbeta = 4.d0
706  double precision, parameter :: weno_eps_machine = 1.0d-18
707 
708  ilm^l=il^l-kr(idims,^d);
709  ilmm^l=ilm^l-kr(idims,^d);
710  ilmmm^l=ilmm^l-kr(idims,^d);
711  ilp^l=il^l+kr(idims,^d);
712  ilpp^l=ilp^l+kr(idims,^d);
713  ilppp^l=ilpp^l+kr(idims,^d);
714  ilpppp^l=ilppp^l+kr(idims,^d);
715 
716  d_array(1:4) = (/ 1.d0/35.d0, 12.d0/35.d0, 18.d0/35.d0, 4.d0/35.d0 /)
717  u1_coeff(1:4) = (/ -1.d0/4.d0, 13.d0/12.d0, -23.d0/12.d0, 25.d0/12.d0 /)
718  u2_coeff(1:4) = (/ 1.d0/12.d0, -5.d0/12.d0, 13.d0/12.d0, 1.d0/4.d0 /)
719  u3_coeff(1:4) = (/ -1.d0/12.d0, 7.d0/12.d0, 7.d0/12.d0, -1.d0/12.d0 /)
720  u4_coeff(1:4) = (/ 1.d0/4.d0, 13.d0/12.d0, -5.d0/12.d0, 1.d0/12.d0 /)
721 
722  !> left side
723  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) &
724  + u1_coeff(4) * w(il^s,1:nwflux)
725  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) &
726  + u2_coeff(4) * w(ilp^s,1:nwflux)
727  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) &
728  + u3_coeff(4) * w(ilpp^s,1:nwflux)
729  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) &
730  + u4_coeff(4) * w(ilppp^s,1:nwflux)
731 
732  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) &
733  - 1854.d0 * w(il^s,1:nwflux)) &
734  + 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)) &
735  + 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
736  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) &
737  - 494.d0 * w(ilp^s,1:nwflux)) &
738  + 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)) &
739  + 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
740  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) &
741  - 494.d0 * w(ilpp^s,1:nwflux)) &
742  + 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)) &
743  + 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
744  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) &
745  - 1854.d0 * w(ilppp^s,1:nwflux)) &
746  + 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)) &
747  + w(ilpp^s,1:nwflux) * (7043.d0 * w(ilpp^s,1:nwflux) - 3882.d0 * w(ilppp^s,1:nwflux)) &
748  + 547.d0 * w(ilppp^s,1:nwflux) ** 2
749 
750  alpha_sum(il^s,1:nwflux) = 0.0d0
751  do i = 1,4
752  alpha_array(il^s,1:nwflux,i) = d_array(i)/(beta(il^s,1:nwflux,i) + weno_eps_machine)**2
753  alpha_sum(il^s,1:nwflux) = alpha_sum(il^s,1:nwflux) + alpha_array(il^s,1:nwflux,i)
754  end do
755  flux(il^s,1:nwflux) = 0.0d0
756  do i = 1,4
757  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))
758  end do
759 
760  select case(var)
761  ! case1 for wenojs, case2 for mpweno
762  case(1)
763  wlc(il^s,1:nwflux) = flux(il^s,1:nwflux)
764  case(2)
765  idmax^d=ilmax^d; idmin^d=ilmin^d-kr(idims,^d);
766  idm^l=id^l-kr(idims,^d);
767  idp^l=id^l+kr(idims,^d);
768 
769  iemax^d=idmax^d+kr(idims,^d); iemin^d=idmin^d;
770  iem^l=ie^l-kr(idims,^d);
771  iep^l=ie^l+kr(idims,^d);
772 
773  d(ie^s,1:nwflux) = w(iep^s,1:nwflux)-2.0d0*w(ie^s,1:nwflux)+w(iem^s,1:nwflux)
774 
775  do iw=1,nwflux
776  a(id^s) = 4.0d0*d(id^s,iw)-d(idp^s,iw)
777  b(id^s) = 4.0d0*d(idp^s,iw)-d(id^s,iw)
778  call minmod(ixi^l,id^l,a,b,tmp)
779  a(id^s) = d(id^s,iw)
780  b(id^s) = d(idp^s,iw)
781  call minmod(ixi^l,id^l,a,b,tmp2)
782  call minmod(ixi^l,id^l,tmp,tmp2,tmp3)
783  dm4(id^s,iw) = tmp3(id^s)
784  end do
785 
786  flux_ul(il^s,1:nwflux) = w(il^s,1:nwflux) + mpalpha * (w(il^s,1:nwflux) - w(ilm^s,1:nwflux))
787  flux_md(il^s,1:nwflux) = half * (w(il^s,1:nwflux) + w(ilp^s,1:nwflux) - dm4(il^s,1:nwflux))
788  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)
789 
790  flux_min(il^s,1:nwflux) = max(min(w(il^s,1:nwflux), w(ilp^s,1:nwflux), flux_md(il^s,1:nwflux)), &
791  min(w(il^s,1:nwflux), flux_ul(il^s,1:nwflux),flux_lc(il^s,1:nwflux)))
792 
793  flux_max(il^s,1:nwflux) = min(max(w(il^s,1:nwflux), w(ilp^s,1:nwflux), flux_md(il^s,1:nwflux)), &
794  max(w(il^s,1:nwflux), flux_ul(il^s,1:nwflux),flux_lc(il^s,1:nwflux)))
795  do iw=1,nwflux
796  a(il^s) = flux(il^s,iw)
797  b(il^s) = flux_min(il^s,iw)
798  c(il^s) = flux_max(il^s,iw)
799  call median(ixi^l, il^l, a, b, c, tmp)
800  wlc(il^s,iw) = tmp(il^s)
801  end do
802  end select
803 
804  !> right side
805  !>> mmm -> pppp
806  !>> mm -> ppp
807  !>> m -> pp
808  !>> 0 -> p
809  !>> p -> 0
810  !>> pp -> m
811  !>> ppp -> mm
812  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) &
813  + u1_coeff(4) * w(ilp^s,1:nwflux)
814  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) &
815  + u2_coeff(4) * w(il^s,1:nwflux)
816  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) &
817  + u3_coeff(4) * w(ilm^s,1:nwflux)
818  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) &
819  + u4_coeff(4) * w(ilmm^s,1:nwflux)
820 
821  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) &
822  - 1854.d0 * w(ilp^s,1:nwflux)) &
823  + 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)) &
824  + 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
825  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) &
826  - 494.d0 * w(il^s,1:nwflux)) &
827  + 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)) &
828  + 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
829  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) &
830  - 494.d0 * w(ilm^s,1:nwflux)) &
831  + 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)) &
832  + 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
833  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) &
834  - 1854.d0 * w(ilmm^s,1:nwflux)) &
835  + 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)) &
836  + 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
837 
838  alpha_sum(il^s,1:nwflux) = 0.0d0
839  do i = 1,4
840  alpha_array(il^s,1:nwflux,i) = d_array(i)/(beta(il^s,1:nwflux,i) + weno_eps_machine)**2
841  alpha_sum(il^s,1:nwflux) = alpha_sum(il^s,1:nwflux) + alpha_array(il^s,1:nwflux,i)
842  end do
843  flux(il^s,1:nwflux) = 0.0d0
844  do i = 1,4
845  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))
846  end do
847 
848  select case(var)
849  case(1)
850  wrc(il^s,1:nwflux) = flux(il^s,1:nwflux)
851  case(2)
852  idmax^d=ilmax^d+kr(idims,^d); idmin^d=ilmin^d;
853  idm^l=id^l-kr(idims,^d);
854  idp^l=id^l+kr(idims,^d);
855 
856  iemax^d=idmax^d; iemin^d=idmin^d-kr(idims,^d);
857  iem^l=ie^l-kr(idims,^d);
858  iep^l=ie^l+kr(idims,^d);
859  iepp^l=iep^l+kr(idims,^d);
860 
861  d(ie^s,1:nwflux) = w(ie^s,1:nwflux)-2.0d0*w(iep^s,1:nwflux)+w(iepp^s,1:nwflux)
862 
863  do iw = 1,nwflux
864  a(id^s) = 4.0d0*d(id^s,iw)-d(idm^s,iw)
865  b(id^s) = 4.0d0*d(idm^s,iw)-d(id^s,iw)
866  call minmod(ixi^l,id^l,a,b,tmp)
867  a(id^s) = d(id^s,iw)
868  b(id^s) = d(idm^s,iw)
869  call minmod(ixi^l,id^l,a,b,tmp2)
870  call minmod(ixi^l,id^l,tmp,tmp2,tmp3)
871  dm4(id^s,iw) = tmp3(id^s)
872  end do
873 
874  flux_ul(il^s,1:nwflux) = w(ilp^s,1:nwflux) + mpalpha * (w(ilp^s,1:nwflux) - w(ilpp^s,1:nwflux))
875  flux_md(il^s,1:nwflux) = half * (w(il^s,1:nwflux) + w(ilp^s,1:nwflux) - dm4(il^s,1:nwflux))
876  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)
877 
878  flux_min(il^s,1:nwflux) = max(min(w(ilp^s,1:nwflux), w(il^s,1:nwflux), flux_md(il^s,1:nwflux)), &
879  min(w(ilp^s,1:nwflux), flux_ul(il^s,1:nwflux),flux_lc(il^s,1:nwflux)))
880 
881  flux_max(il^s,1:nwflux) = min(max(w(ilp^s,1:nwflux), w(il^s,1:nwflux), flux_md(il^s,1:nwflux)), &
882  max(w(ilp^s,1:nwflux), flux_ul(il^s,1:nwflux),flux_lc(il^s,1:nwflux)))
883  do iw=1,nwflux
884  a(il^s) = flux(il^s,iw)
885  b(il^s) = flux_min(il^s,iw)
886  c(il^s) = flux_max(il^s,iw)
887  call median(ixi^l, il^l, a, b, c, tmp)
888  wrc(il^s,iw) = tmp(il^s)
889  end do
890  end select
891  end subroutine weno7limiter
892 
893  subroutine exeno7limiter(ixI^L,iL^L,idims,w,wLC,wRC)
895 
896  integer, intent(in) :: ixI^L, iL^L, idims
897  double precision, intent(in) :: w(ixi^s,1:nw)
898  double precision, intent(inout) :: wRC(ixi^s,1:nw),wLC(ixi^s,1:nw)
899  !> local
900  integer :: i, iw
901  integer :: iLm^L, iLmm^L, iLmmm^L
902  integer :: iLp^L, iLpp^L, iLppp^L, iLpppp^L
903  integer :: iM^L, iMm^L, iMmm^L
904  integer :: iMp^L, iMpp^L, iMppp^L
905  integer :: id^L, idp^L, idpp^L, idm^L, ie^L, iem^L, iep^L, iepp^L
906  integer, dimension(ixI^S,1:nw) :: delta_sum
907  integer, dimension(ixI^S,1:nw,3):: delta
908  double precision, dimension(2) :: beta_coeff
909  double precision, dimension(3) :: d_array
910  double precision, dimension(ixI^S,1:nw) :: gamma_sum, flux
911  double precision, dimension(ixI^S,1:nw,3) :: beta, gamma_array, kai_array
912  double precision, parameter :: exeno_ct = 1.d-1
913  double precision, parameter :: weno_eps_machine = 1.d-18
914 
915  ilm^l=il^l-kr(idims,^d);
916  ilmm^l=ilm^l-kr(idims,^d);
917  ilmmm^l=ilmm^l-kr(idims,^d);
918  ilp^l=il^l+kr(idims,^d);
919  ilpp^l=ilp^l+kr(idims,^d);
920  ilppp^l=ilpp^l+kr(idims,^d);
921  ilpppp^l=ilppp^l+kr(idims,^d);
922 
923  immin^d=ilmin^d-kr(idims,^d);
924  immax^d=ilmax^d+kr(idims,^d);
925  imm^l=im^l-kr(idims,^d);
926  immm^l=imm^l-kr(idims,^d);
927  imp^l=im^l+kr(idims,^d);
928  impp^l=imp^l+kr(idims,^d);
929  imppp^l=impp^l+kr(idims,^d);
930 
931  beta_coeff(1:2) = (/ 1.0833333333333333d0, 0.25d0/)
932  d_array(1:3) = (/ 1.0d0/10.0d0, 3.0d0/5.0d0, 3.0d0/10.0d0 /)
933 
934  !>> left side
935  beta(im^s,1:nwflux,1) = beta_coeff(1) * (w(immm^s,1:nwflux) + w(im^s,1:nwflux) - 2.0d0 * w(imm^s,1:nwflux))**2 &
936  + beta_coeff(2) * (w(immm^s,1:nwflux) - 4.0d0 * w(imm^s,1:nwflux) + 3.0d0 * w(im^s,1:nwflux))**2
937  beta(im^s,1:nwflux,2) = beta_coeff(1) * (w(imm^s,1:nwflux) + w(imp^s,1:nwflux) - 2.0d0 * w(im^s,1:nwflux))**2 &
938  + beta_coeff(2) * (w(imm^s,1:nwflux) - w(imp^s,1:nwflux))**2
939  beta(im^s,1:nwflux,3) = beta_coeff(1) * (w(im^s,1:nwflux) + w(impp^s,1:nwflux) - 2.0d0 * w(imp^s,1:nwflux))**2 &
940  + beta_coeff(2) * (3.0d0 * w(im^s, 1:nwflux) - 4.0d0 * w(imp^s,1:nwflux) + w(impp^s,1:nwflux))**2
941 
942  gamma_sum(im^s,1:nwflux) = 0.0d0
943  do i = 1,3
944  gamma_array(im^s,1:nwflux,i) = d_array(i) / (beta(im^s,1:nwflux,i) + weno_eps_machine)**2
945  gamma_sum(im^s,1:nwflux) = gamma_sum(im^s,1:nwflux) + gamma_array(im^s,1:nwflux,i)
946  end do
947  do i = 1,3
948  kai_array(im^s,1:nwflux,i) = gamma_array(im^s,1:nwflux,i) / gamma_sum(im^s,1:nwflux)
949  where(kai_array(im^s,1:nwflux,i) .lt. exeno_ct)
950  delta(im^s,1:nwflux,i) = 0
951  elsewhere
952  delta(im^s,1:nwflux,i) = 1
953  endwhere
954  end do
955 
956  delta_sum(il^s,1:nwflux) = delta(ilm^s,1:nwflux,1) * 1 + delta(ilp^s,1:nwflux,3) * 2 + delta(il^s,1:nwflux,1) * 4 &
957  + delta(il^s,1:nwflux,2) * 8 + delta(il^s,1:nwflux,3) * 16
958 
959  !> f3
960  where(delta_sum(il^s,1:nwflux) .eq. 31)
961  flux(il^s,1:nwflux) = (- 3.d0 * w(ilmmm^s,1:nwflux) + 25.d0 * w(ilmm^s,1:nwflux) - 101.d0 * w(ilm^s,1:nwflux) + 319.d0 * w(il^s,1:nwflux) &
962  + 214.d0 * w(ilp^s,1:nwflux) - 38.d0 * w(ilpp^s,1:nwflux) + 4.d0 * w(ilppp^s,1:nwflux)) / 420.d0
963  !> f4
964  elsewhere(delta_sum(il^s,1:nwflux) .eq. 30)
965  flux(il^s,1:nwflux) = (w(ilmm^s,1:nwflux) - 8.d0 * w(ilm^s,1:nwflux) + 37.d0 * w(il^s,1:nwflux) + 37.d0 * w(ilp^s,1:nwflux) &
966  - 8.d0 * w(ilpp^s,1:nwflux) + w(ilppp^s,1:nwflux)) / 60.d0
967  !> f5
968  elsewhere(delta_sum(il^s,1:nwflux) .eq. 29)
969  flux(il^s,1:nwflux) = (- w(ilmmm^s,1:nwflux) + 7.d0 * w(ilmm^s,1:nwflux) - 23.d0 * w(ilm^s,1:nwflux) + 57.d0 * w(il^s,1:nwflux) &
970  + 22.d0 * w(ilp^s,1:nwflux) - 2.d0 * w(ilpp^s,1:nwflux)) / 60.d0
971  !> f6
972  elsewhere(delta_sum(il^s,1:nwflux) .eq. 28)
973  flux(il^s,1:nwflux) = (2.d0 * w(ilmm^s,1:nwflux) - 13.d0 * w(ilm^s,1:nwflux) + 47.d0 * w(il^s,1:nwflux) + 27.d0 * w(ilp^s,1:nwflux) &
974  - 3.d0 * w(ilpp^s,1:nwflux)) / 60.d0
975  !> f7
976  elsewhere(delta_sum(il^s,1:nwflux) .ge. 24)
977  flux(il^s,1:nwflux) = (- w(ilm^s,1:nwflux) + 7.d0 * w(il^s,1:nwflux) + 7.d0 * w(ilp^s,1:nwflux) - w(ilpp^s,1:nwflux)) / 12.d0
978  !> f9
979  elsewhere(delta_sum(il^s,1:nwflux) .ge. 16)
980  flux(il^s,1:nwflux) = (2.d0 * w(il^s,1:nwflux) + 5.d0 * w(ilp^s,1:nwflux) - w(ilpp^s,1:nwflux)) / 6.d0
981  !> f8
982  elsewhere(delta_sum(il^s,1:nwflux) .ge. 12)
983  flux(il^s,1:nwflux) = (w(ilmm^s,1:nwflux) - 5.d0 * w(ilm^s,1:nwflux) + 13.d0 * w(il^s,1:nwflux) + 3.d0 * w(ilp^s,1:nwflux)) / 12.d0
984  !> f10
985  elsewhere(delta_sum(il^s,1:nwflux) .ge. 8)
986  flux(il^s,1:nwflux) = (- w(ilm^s,1:nwflux) + 5.d0 * w(il^s,1:nwflux) + 2.d0 * w(ilp^s,1:nwflux)) / 6.d0
987  !> f11
988  elsewhere
989  flux(il^s,1:nwflux) = (2.d0 * w(ilmm^s,1:nwflux) - 7.d0 * w(ilm^s,1:nwflux) + 11.d0 * w(il^s,1:nwflux)) / 6.d0
990  endwhere
991 
992  wlc(il^s,1:nwflux) = flux(il^s,1:nwflux)
993 
994  !> right side
995  beta(im^s,1:nwflux,1) = beta_coeff(1) * (w(imppp^s,1:nwflux) + w(imp^s,1:nwflux) - 2.0d0 * w(impp^s,1:nwflux))**2 &
996  + beta_coeff(2) * (w(imppp^s,1:nwflux) - 4.0d0 * w(impp^s,1:nwflux) + 3.0d0 * w(imp^s,1:nwflux))**2
997  beta(im^s,1:nwflux,2) = beta_coeff(1) * (w(impp^s,1:nwflux) + w(im^s,1:nwflux) - 2.0d0 * w(imp^s,1:nwflux))**2 &
998  + beta_coeff(2) * (w(impp^s,1:nwflux) - w(im^s,1:nwflux))**2
999  beta(im^s,1:nwflux,3) = beta_coeff(1) * (w(imp^s,1:nwflux) + w(imm^s,1:nwflux) - 2.0d0 * w(im^s,1:nwflux))**2 &
1000  + beta_coeff(2) * (3.0d0 * w(imp^s, 1:nwflux) - 4.0d0 * w(im^s,1:nwflux) + w(imm^s,1:nwflux))**2
1001 
1002  gamma_sum(im^s,1:nwflux) = 0.0d0
1003  do i = 1,3
1004  gamma_array(im^s,1:nwflux,i) = d_array(i) / (beta(im^s,1:nwflux,i) + weno_eps_machine)**2
1005  gamma_sum(im^s,1:nwflux) = gamma_sum(im^s,1:nwflux) + gamma_array(im^s,1:nwflux,i)
1006  end do
1007  do i = 1,3
1008  kai_array(im^s,1:nwflux,i) = gamma_array(im^s,1:nwflux,i) / gamma_sum(im^s,1:nwflux)
1009  where(kai_array(im^s,1:nwflux,i) .lt. exeno_ct)
1010  delta(im^s,1:nwflux,i) = 0
1011  elsewhere
1012  delta(im^s,1:nwflux,i) = 1
1013  endwhere
1014  end do
1015 
1016  delta_sum(il^s,1:nwflux) = delta(ilp^s,1:nwflux,1) * 1 + delta(ilm^s,1:nwflux,3) * 2 + delta(il^s,1:nwflux,1) * 4 &
1017  + delta(il^s,1:nwflux,2) * 8 + delta(il^s,1:nwflux,3) * 16
1018 
1019  where(delta_sum(il^s,1:nwflux) .eq. 31)
1020  flux(il^s,1:nwflux) = (- 3.d0 * w(ilpppp^s,1:nwflux) + 25.d0 * w(ilppp^s,1:nwflux) - 101.d0 * w(ilpp^s,1:nwflux) &
1021  + 319.d0 * w(ilp^s,1:nwflux) + 214.d0 * w(il^s,1:nwflux) - 38.d0 * w(ilm^s,1:nwflux) &
1022  + 4.d0 * w(ilmm^s,1:nwflux)) / 420.d0
1023  elsewhere(delta_sum(il^s,1:nwflux) .eq. 30)
1024  flux(il^s,1:nwflux) = (w(ilppp^s,1:nwflux) - 8.d0 * w(ilpp^s,1:nwflux) + 37.d0 * w(ilp^s,1:nwflux) + 37.d0 * w(il^s,1:nwflux) &
1025  - 8.d0 * w(ilm^s,1:nwflux) + w(ilmm^s,1:nwflux)) / 60.d0
1026  elsewhere(delta_sum(il^s,1:nwflux) .eq. 29)
1027  flux(il^s,1:nwflux) = (- w(ilpppp^s,1:nwflux) + 7.d0 * w(ilppp^s,1:nwflux) - 23.d0 * w(ilpp^s,1:nwflux) + 57.d0 * w(ilp^s,1:nwflux) &
1028  + 22.d0 * w(il^s,1:nwflux) - 2.d0 * w(ilm^s,1:nwflux)) / 60.d0
1029  elsewhere(delta_sum(il^s,1:nwflux) .eq. 28)
1030  flux(il^s,1:nwflux) = (2.d0 * w(ilppp^s,1:nwflux) - 13.d0 * w(ilpp^s,1:nwflux) + 47.d0 * w(ilp^s,1:nwflux) + 27.d0 * w(il^s,1:nwflux) &
1031  - 3.d0 * w(ilm^s,1:nwflux)) / 60.d0
1032  elsewhere(delta_sum(il^s,1:nwflux) .ge. 24)
1033  flux(il^s,1:nwflux) = (- w(ilpp^s,1:nwflux) + 7.d0 * w(ilp^s,1:nwflux) + 7.d0 * w(il^s,1:nwflux) - w(ilm^s,1:nwflux)) / 12.d0
1034  elsewhere(delta_sum(il^s,1:nwflux) .ge. 16)
1035  flux(il^s,1:nwflux) = (2.d0 * w(ilp^s,1:nwflux) + 5.d0 * w(il^s,1:nwflux) - w(ilm^s,1:nwflux)) / 6.d0
1036  elsewhere(delta_sum(il^s,1:nwflux) .ge. 12)
1037  flux(il^s,1:nwflux) = (w(ilppp^s,1:nwflux) - 5.d0 * w(ilpp^s,1:nwflux) + 13.d0 * w(ilp^s,1:nwflux) + 3.d0 * w(il^s,1:nwflux)) / 12.d0
1038  elsewhere(delta_sum(il^s,1:nwflux) .ge. 8)
1039  flux(il^s,1:nwflux) = (- w(ilpp^s,1:nwflux) + 5.d0 * w(ilp^s,1:nwflux) + 2.d0 * w(il^s,1:nwflux)) / 6.d0
1040  elsewhere
1041  flux(il^s,1:nwflux) = (2.d0 * w(ilppp^s,1:nwflux) - 7.d0 * w(ilpp^s,1:nwflux) + 11.d0 * w(ilp^s,1:nwflux)) / 6.d0
1042  endwhere
1043 
1044  wrc(il^s,1:nwflux) = flux(il^s,1:nwflux)
1045 
1046  end subroutine exeno7limiter
1047 
1048  subroutine minmod(ixI^L,ixO^L,a,b,minm)
1051 
1052  integer, intent(in) :: ixI^L, ixO^L
1053  double precision, intent(in) :: a(ixi^s), b(ixi^s)
1054  double precision, intent(out):: minm(ixi^s)
1055 
1056  minm(ixo^s) = (sign(one,a(ixo^s))+sign(one,b(ixo^s)))/2.0d0 * &
1057  min(abs(a(ixo^s)),abs(b(ixo^s)))
1058 
1059  end subroutine minmod
1060 
1061  subroutine median(ixI^L,ixO^L,a,b,c,med)
1064 
1065  integer, intent(in) :: ixI^L, ixO^L
1066  double precision, intent(in) :: a(ixi^s), b(ixi^s), c(ixi^s)
1067  double precision, intent(out):: med(ixi^s)
1068 
1069  double precision :: tmp1(ixi^s),tmp2(ixi^s)
1070 
1071  tmp1(ixo^s) = b(ixo^s) - a(ixo^s); tmp2(ixo^s) = c(ixo^s) - a(ixo^s)
1072 
1073  med(ixo^s) = a(ixo^s) + (sign(one,tmp1(ixo^s))+sign(one,tmp2(ixo^s)))/2.0d0 * &
1074  min(abs(tmp1(ixo^s)),abs(tmp2(ixo^s)))
1075 
1076  end subroutine median
1077 end module mod_weno
This module contains definitions of global parameters and variables and some generic functions/subrou...
subroutine, public weno3limiter(ixIL, iLL, idims, dxdim, w, wLC, wRC, var)
Definition: mod_weno.t:31
subroutine, public weno5nmlimiter(ixIL, iLL, idims, dxdim, w, wLC, wRC, var)
Definition: mod_weno.t:244
subroutine, public exeno7limiter(ixIL, iLL, idims, w, wLC, wRC)
Definition: mod_weno.t:894
subroutine, public weno5limiterl(ixIL, iLL, idims, w, wLC, var)
Definition: mod_weno.t:368
subroutine median(ixIL, ixOL, a, b, c, med)
Definition: mod_weno.t:1062
subroutine minmod(ixIL, ixOL, a, b, minm)
Definition: mod_weno.t:1049
subroutine, public weno5limiterr(ixIL, iLL, idims, w, wRC, var)
Definition: mod_weno.t:451
integer, dimension(3, 3) kr
Kronecker delta tensor.
integer, dimension(:), allocatable, parameter d
subroutine, public weno7limiter(ixIL, iLL, idims, w, wLC, wRC, var)
Definition: mod_weno.t:689
subroutine, public weno5nmlimiterl(ixIL, iLL, idims, w, wLC, var)
Definition: mod_weno.t:533
subroutine, public weno5limiter(ixIL, iLL, idims, dxdim, w, wLC, wRC, var)
Definition: mod_weno.t:119
subroutine, public weno5nmlimiterr(ixIL, iLL, idims, w, wRC, var)
Definition: mod_weno.t:611