MPI-AMRVAC  3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
mod_mp5.t
Go to the documentation of this file.
1 !> Module containing the MP5 (fifth order) flux scheme
2 module mod_mp5
3 
4  implicit none
5  private
6 
7  public :: mp5limiter
8  public :: mp5limitervar
9  public :: mp5limiterl
10  public :: mp5limiterr
11 
12 contains
13 
14  !> MP5 limiter from Suresh & Huynh 1997 Following the convention of Mignone et
15  !> al. 2010. Needs at least three ghost cells
16  subroutine mp5limiter(ixI^L,iL^L,idims,w,wLC,wRC)
17  use mod_global_parameters
18  use mod_weno, only: fix_limiter1
19 
20  integer, intent(in) :: ixi^l, il^l, idims
21  double precision, intent(in) :: w(ixi^s,1:nw)
22 
23  double precision, intent(inout) :: wrc(ixi^s,1:nw),wlc(ixi^s,1:nw)
24  ! .. local ..
25  double precision, dimension(ixI^S,1:nw) :: f, fmp, fmin, fmax, ful, dm4, d, fmd, flc, flim
26  double precision, dimension(ixI^S,1:nw) :: wrctmp, wlctmp
27  double precision, dimension(ixI^S) :: tmp, tmp2, tmp3, a, b, c
28  double precision, parameter :: eps=0.d0, alpha=4.0d0
29  !double precision :: alpha
30  integer :: ilm^l, ilmm^l, ilp^l, ilpp^l, ilppp^l
31  integer :: id^l, idp^l, idpp^l, idm^l, ie^l, iem^l, iep^l, iepp^l
32  integer :: iw
33  !----------------------------------------------------------------------------
34 
35  ! Variable alpha:
36  !alpha = float(nstep)/courantpar - one
37 
38  ! Left side:
39  ! range to process:
40  !iLmin^D=ixmin^D-kr(idims,^D);iLmax^D=ixmax^D;
41 
42  ! HALL
43  ! For Hall, we need one more reconstructed layer since currents are computed in getflux:
44  ! also add one ghost zone!
45  ! {iL^L=iL^L^LADD1;}
46 
47  ! iL^L holds the indices of interfaces to reconstruct to. Convention is that a center index holds the _right-side_ interface.
48 
49  ilm^l=il^l-kr(idims,^d);
50  ilmm^l=ilm^l-kr(idims,^d);
51  ilp^l=il^l+kr(idims,^d);
52  ilpp^l=ilp^l+kr(idims,^d);
53 
54  f(il^s,1:nwflux) = 1.0d0/60.0d0 * (&
55  2.0d0* w(ilmm^s,1:nwflux) &
56  - 13.0d0* w(ilm^s,1:nwflux) &
57  + 47.0d0* w(il^s,1:nwflux) &
58  + 27.0d0* w(ilp^s,1:nwflux) &
59  - 3.0d0* w(ilpp^s,1:nwflux))
60 
61  ! get fmp and ful:
62  do iw=1,nwflux
63  a(il^s) = w(ilp^s,iw)-w(il^s,iw)
64  b(il^s) = alpha*(w(il^s,iw)-w(ilm^s,iw))
65  call minmod(ixi^l,il^l,a,b,tmp)
66  fmp(il^s,iw) = w(il^s,iw) + tmp(il^s)
67  ful(il^s,iw) = w(il^s,iw) + b(il^s)
68  end do ! iw loop
69 
70  ! get dm4:
71  idmax^d=ilmax^d; idmin^d=ilmin^d-kr(idims,^d);
72  idm^l=id^l-kr(idims,^d);
73  idp^l=id^l+kr(idims,^d);
74 
75  iemax^d=idmax^d+kr(idims,^d); iemin^d=idmin^d;
76  iem^l=ie^l-kr(idims,^d);
77  iep^l=ie^l+kr(idims,^d);
78 
79  d(ie^s,1:nwflux) = w(iep^s,1:nwflux)-2.0d0*w(ie^s,1:nwflux)+w(iem^s,1:nwflux)
80 
81  do iw=1,nwflux
82  a(id^s) = 4.0d0*d(id^s,iw)-d(idp^s,iw)
83  b(id^s) = 4.0d0*d(idp^s,iw)-d(id^s,iw)
84  call minmod(ixi^l,id^l,a,b,tmp)
85  a(id^s) = d(id^s,iw)
86  b(id^s) = d(idp^s,iw)
87  call minmod(ixi^l,id^l,a,b,tmp2)
88  call minmod(ixi^l,id^l,tmp,tmp2,tmp3)
89  dm4(id^s,iw) = tmp3(id^s)
90  end do
91 
92  ! get fmd:
93  fmd(il^s,1:nwflux) = (w(il^s,1:nwflux)+w(ilp^s,1:nwflux))/2.0d0-dm4(il^s,1:nwflux)/2.0d0
94 
95  !get flc:
96  flc(il^s,1:nwflux) = half*(3.0d0*w(il^s,1:nwflux) &
97  - w(ilm^s,1:nwflux)) + 4.0d0/3.0d0*dm4(ilm^s,1:nwflux)
98 
99  fmin(il^s,1:nwflux) = max(min(w(il^s,1:nwflux),w(ilp^s,1:nwflux),fmd(il^s,1:nwflux)),&
100  min(w(il^s,1:nwflux),ful(il^s,1:nwflux),flc(il^s,1:nwflux)))
101 
102  fmax(il^s,1:nwflux) = min(max(w(il^s,1:nwflux),w(ilp^s,1:nwflux),fmd(il^s,1:nwflux)),&
103  max(w(il^s,1:nwflux),ful(il^s,1:nwflux),flc(il^s,1:nwflux)))
104 
105  do iw=1,nwflux
106  a(il^s) = fmin(il^s,iw)
107  b(il^s) = f(il^s,iw)
108  c(il^s) = fmax(il^s,iw)
109  call median(ixi^l,il^l,a,b,c,tmp)
110  flim(il^s,iw) = tmp(il^s)
111  end do
112 
113  ! check case
114  where ((f(il^s,1:nwflux)-w(il^s,1:nwflux))*(f(il^s,1:nwflux)-fmp(il^s,1:nwflux)) .le. eps)
115  wlctmp(il^s,1:nwflux) = f(il^s,1:nwflux)
116  elsewhere
117  wlctmp(il^s,1:nwflux) = flim(il^s,1:nwflux)
118  end where
119 
120  ! Right side:
121  ! the interpolation from the right is obtained when the left-hand formula is applied to
122  ! data mirrored about the interface.
123  ! thus substitute:
124  ! i-2 -> i+3
125  ! i-1 -> i+2
126  ! i -> i+1
127  ! i+1 -> i
128  ! i+2 -> i-1
129 
130  ilppp^l=ilpp^l+kr(idims,^d);
131 
132  f(il^s,1:nwflux) = 1.0d0/60.0d0 * (&
133  2.0d0* w(ilppp^s,1:nwflux) &
134  - 13.0d0* w(ilpp^s,1:nwflux) &
135  + 47.0d0* w(ilp^s,1:nwflux) &
136  + 27.0d0* w(il^s,1:nwflux) &
137  - 3.0d0* w(ilm^s,1:nwflux))
138 
139  ! get fmp and ful:
140  do iw=1,nwflux
141  a(il^s) = w(il^s,iw)-w(ilp^s,iw)
142  b(il^s) = alpha*(w(ilp^s,iw)-w(ilpp^s,iw))
143  call minmod(ixi^l,il^l,a,b,tmp)
144  fmp(il^s,iw) = w(ilp^s,iw) + tmp(il^s)
145  ful(il^s,iw) = w(ilp^s,iw) + b(il^s)
146  end do ! iw loop
147 
148  ! get dm4:
149  idmax^d=ilmax^d+kr(idims,^d); idmin^d=ilmin^d;
150  idm^l=id^l-kr(idims,^d);
151  idp^l=id^l+kr(idims,^d);
152 
153  iemax^d=idmax^d; iemin^d=idmin^d-kr(idims,^d);
154  iem^l=ie^l-kr(idims,^d);
155  iep^l=ie^l+kr(idims,^d);
156  iepp^l=iep^l+kr(idims,^d);
157 
158  d(ie^s,1:nwflux) = w(ie^s,1:nwflux)-2.0d0*w(iep^s,1:nwflux)+w(iepp^s,1:nwflux)
159 
160  do iw=1,nwflux
161  a(id^s) = 4.0d0*d(id^s,iw)-d(idm^s,iw)
162  b(id^s) = 4.0d0*d(idm^s,iw)-d(id^s,iw)
163  call minmod(ixi^l,id^l,a,b,tmp)
164  a(id^s) = d(id^s,iw)
165  b(id^s) = d(idm^s,iw)
166  call minmod(ixi^l,id^l,a,b,tmp2)
167  call minmod(ixi^l,id^l,tmp,tmp2,tmp3)
168  dm4(id^s,iw) = tmp3(id^s)
169  end do
170 
171  ! get fmd:
172  fmd(il^s,1:nwflux) = (w(il^s,1:nwflux)+w(ilp^s,1:nwflux))/2.0d0-dm4(il^s,1:nwflux)/2.0d0
173 
174  !get flc:
175  flc(il^s,1:nwflux) = half*(3.0d0*w(ilp^s,1:nwflux) &
176  - w(ilpp^s,1:nwflux)) + 4.0d0/3.0d0*dm4(ilp^s,1:nwflux)
177 
178  fmin(il^s,1:nwflux) = max(min(w(ilp^s,1:nwflux),w(il^s,1:nwflux),fmd(il^s,1:nwflux)),&
179  min(w(ilp^s,1:nwflux),ful(il^s,1:nwflux),flc(il^s,1:nwflux)))
180 
181  fmax(il^s,1:nwflux) = min(max(w(ilp^s,1:nwflux),w(il^s,1:nwflux),fmd(il^s,1:nwflux)),&
182  max(w(ilp^s,1:nwflux),ful(il^s,1:nwflux),flc(il^s,1:nwflux)))
183 
184  do iw=1,nwflux
185  a(il^s) = fmin(il^s,iw)
186  b(il^s) = f(il^s,iw)
187  c(il^s) = fmax(il^s,iw)
188  call median(ixi^l,il^l,a,b,c,tmp)
189  flim(il^s,iw) = tmp(il^s)
190  end do
191 
192  ! check case
193  where ((f(il^s,1:nwflux)-w(ilp^s,1:nwflux))*(f(il^s,1:nwflux)-fmp(il^s,1:nwflux)) .le. eps)
194  wrctmp(il^s,1:nwflux) = f(il^s,1:nwflux)
195  elsewhere
196  wrctmp(il^s,1:nwflux) = flim(il^s,1:nwflux)
197  end where
198 
199  call fix_limiter1(ixi^l,il^l,wlctmp,wrctmp,wlc,wrc)
200 
201  end subroutine mp5limiter
202 
203  subroutine mp5limiterl(ixI^L,iL^L,idims,w,wLC)
204  use mod_global_parameters
205  !use mod_weno, only: fix_onelimiter1
206 
207  integer, intent(in) :: ixi^l, il^l, idims
208  double precision, intent(in) :: w(ixi^s,1:nw)
209 
210  double precision, intent(inout) :: wlc(ixi^s,1:nw)
211  ! .. local ..
212  double precision, dimension(ixI^S,1:nw) :: f, fmp, fmin, fmax, ful, dm4, d, fmd, flc, flim
213  double precision, dimension(ixI^S) :: tmp, tmp2, tmp3, a, b, c
214  double precision, parameter :: eps=0.d0, alpha=4.0d0
215  !double precision :: alpha
216  !double precision, dimension(ixI^S,1:nw) :: wLCtmp
217  integer :: ilm^l, ilmm^l, ilp^l, ilpp^l
218  integer :: id^l, idp^l, idpp^l, idm^l, ie^l, iem^l, iep^l, iepp^l
219  integer :: iw
220 
221  ! Variable alpha:
222  !alpha = float(nstep)/courantpar - one
223 
224  ! Left side:
225  ilm^l=il^l-kr(idims,^d);
226  ilmm^l=ilm^l-kr(idims,^d);
227  ilp^l=il^l+kr(idims,^d);
228  ilpp^l=ilp^l+kr(idims,^d);
229 
230  f(il^s,1:nwflux) = 1.0d0/60.0d0 * (&
231  2.0d0* w(ilmm^s,1:nwflux) &
232  - 13.0d0* w(ilm^s,1:nwflux) &
233  + 47.0d0* w(il^s,1:nwflux) &
234  + 27.0d0* w(ilp^s,1:nwflux) &
235  - 3.0d0* w(ilpp^s,1:nwflux))
236 
237  ! get fmp and ful:
238  do iw=1,nwflux
239  a(il^s) = w(ilp^s,iw)-w(il^s,iw)
240  b(il^s) = alpha*(w(il^s,iw)-w(ilm^s,iw))
241  call minmod(ixi^l,il^l,a,b,tmp)
242  fmp(il^s,iw) = w(il^s,iw) + tmp(il^s)
243  ful(il^s,iw) = w(il^s,iw) + b(il^s)
244  end do ! iw loop
245 
246  ! get dm4:
247  idmax^d=ilmax^d; idmin^d=ilmin^d-kr(idims,^d);
248  idm^l=id^l-kr(idims,^d);
249  idp^l=id^l+kr(idims,^d);
250 
251  iemax^d=idmax^d+kr(idims,^d); iemin^d=idmin^d;
252  iem^l=ie^l-kr(idims,^d);
253  iep^l=ie^l+kr(idims,^d);
254 
255  d(ie^s,1:nwflux) = w(iep^s,1:nwflux)-2.0d0*w(ie^s,1:nwflux)+w(iem^s,1:nwflux)
256 
257  do iw=1,nwflux
258  a(id^s) = 4.0d0*d(id^s,iw)-d(idp^s,iw)
259  b(id^s) = 4.0d0*d(idp^s,iw)-d(id^s,iw)
260  call minmod(ixi^l,id^l,a,b,tmp)
261  a(id^s) = d(id^s,iw)
262  b(id^s) = d(idp^s,iw)
263  call minmod(ixi^l,id^l,a,b,tmp2)
264  call minmod(ixi^l,id^l,tmp,tmp2,tmp3)
265  dm4(id^s,iw) = tmp3(id^s)
266  end do
267 
268  ! get fmd:
269  fmd(il^s,1:nwflux) = (w(il^s,1:nwflux)+w(ilp^s,1:nwflux))/2.0d0-dm4(il^s,1:nwflux)/2.0d0
270 
271  !get flc:
272  flc(il^s,1:nwflux) = half*(3.0d0*w(il^s,1:nwflux) &
273  - w(ilm^s,1:nwflux)) + 4.0d0/3.0d0*dm4(ilm^s,1:nwflux)
274 
275  fmin(il^s,1:nwflux) = max(min(w(il^s,1:nwflux),w(ilp^s,1:nwflux),fmd(il^s,1:nwflux)),&
276  min(w(il^s,1:nwflux),ful(il^s,1:nwflux),flc(il^s,1:nwflux)))
277 
278  fmax(il^s,1:nwflux) = min(max(w(il^s,1:nwflux),w(ilp^s,1:nwflux),fmd(il^s,1:nwflux)),&
279  max(w(il^s,1:nwflux),ful(il^s,1:nwflux),flc(il^s,1:nwflux)))
280 
281  do iw=1,nwflux
282  a(il^s) = fmin(il^s,iw)
283  b(il^s) = f(il^s,iw)
284  c(il^s) = fmax(il^s,iw)
285  call median(ixi^l,il^l,a,b,c,tmp)
286  flim(il^s,iw) = tmp(il^s)
287  end do
288 
289  ! check case
290  where ((f(il^s,1:nwflux)-w(il^s,1:nwflux))*(f(il^s,1:nwflux)-fmp(il^s,1:nwflux)) .le. eps)
291  wlc(il^s,1:nwflux) = f(il^s,1:nwflux)
292  elsewhere
293  wlc(il^s,1:nwflux) = flim(il^s,1:nwflux)
294  end where
295 
296  !call fix_onelimiter1(ixI^L,iL^L,wLCtmp,wLC)
297 
298  end subroutine mp5limiterl
299 
300  subroutine mp5limiterr(ixI^L,iL^L,idims,w,wRC)
301  use mod_global_parameters
302  !use mod_weno, only: fix_onelimiter1
303 
304  integer, intent(in) :: ixi^l, il^l, idims
305  double precision, intent(in) :: w(ixi^s,1:nw)
306 
307  double precision, intent(inout) :: wrc(ixi^s,1:nw)
308  ! .. local ..
309  double precision, dimension(ixI^S,1:nw) :: f, fmp, fmin, fmax, ful, dm4, d, fmd, flc, flim
310  double precision, dimension(ixI^S) :: tmp, tmp2, tmp3, a, b, c
311  double precision, parameter :: eps=0.d0, alpha=4.0d0
312  integer :: ilm^l, ilp^l, ilpp^l, ilppp^l
313  integer :: id^l, idp^l, idpp^l, idm^l, ie^l, iem^l, iep^l, iepp^l
314  integer :: iw
315  !double precision :: alpha
316  !double precision, dimension(ixI^S,1:nw) :: wRCtmp
317 
318  ! Right side:
319  ! the interpolation from the right is obtained when the left-hand formula is applied to
320  ! data mirrored about the interface.
321  ! thus substitute:
322  ! i-2 -> i+3
323  ! i-1 -> i+2
324  ! i -> i+1
325  ! i+1 -> i
326  ! i+2 -> i-1
327 
328  ilm^l=il^l-kr(idims,^d);
329  ilp^l=il^l+kr(idims,^d);
330  ilpp^l=ilp^l+kr(idims,^d);
331  ilppp^l=ilpp^l+kr(idims,^d);
332 
333  f(il^s,1:nwflux) = 1.0d0/60.0d0 * (&
334  2.0d0* w(ilppp^s,1:nwflux) &
335  - 13.0d0* w(ilpp^s,1:nwflux) &
336  + 47.0d0* w(ilp^s,1:nwflux) &
337  + 27.0d0* w(il^s,1:nwflux) &
338  - 3.0d0* w(ilm^s,1:nwflux))
339 
340  ! get fmp and ful:
341  do iw=1,nwflux
342  a(il^s) = w(il^s,iw)-w(ilp^s,iw)
343  b(il^s) = alpha*(w(ilp^s,iw)-w(ilpp^s,iw))
344  call minmod(ixi^l,il^l,a,b,tmp)
345  fmp(il^s,iw) = w(ilp^s,iw) + tmp(il^s)
346  ful(il^s,iw) = w(ilp^s,iw) + b(il^s)
347  end do ! iw loop
348 
349  ! get dm4:
350  idmax^d=ilmax^d+kr(idims,^d); idmin^d=ilmin^d;
351  idm^l=id^l-kr(idims,^d);
352  idp^l=id^l+kr(idims,^d);
353 
354  iemax^d=idmax^d; iemin^d=idmin^d-kr(idims,^d);
355  iem^l=ie^l-kr(idims,^d);
356  iep^l=ie^l+kr(idims,^d);
357  iepp^l=iep^l+kr(idims,^d);
358 
359  d(ie^s,1:nwflux) = w(ie^s,1:nwflux)-2.0d0*w(iep^s,1:nwflux)+w(iepp^s,1:nwflux)
360 
361  do iw=1,nwflux
362  a(id^s) = 4.0d0*d(id^s,iw)-d(idm^s,iw)
363  b(id^s) = 4.0d0*d(idm^s,iw)-d(id^s,iw)
364  call minmod(ixi^l,id^l,a,b,tmp)
365  a(id^s) = d(id^s,iw)
366  b(id^s) = d(idm^s,iw)
367  call minmod(ixi^l,id^l,a,b,tmp2)
368  call minmod(ixi^l,id^l,tmp,tmp2,tmp3)
369  dm4(id^s,iw) = tmp3(id^s)
370  end do
371 
372  ! get fmd:
373  fmd(il^s,1:nwflux) = (w(il^s,1:nwflux)+w(ilp^s,1:nwflux))/2.0d0-dm4(il^s,1:nwflux)/2.0d0
374 
375  !get flc:
376  flc(il^s,1:nwflux) = half*(3.0d0*w(ilp^s,1:nwflux) &
377  - w(ilpp^s,1:nwflux)) + 4.0d0/3.0d0*dm4(ilp^s,1:nwflux)
378 
379  fmin(il^s,1:nwflux) = max(min(w(ilp^s,1:nwflux),w(il^s,1:nwflux),fmd(il^s,1:nwflux)),&
380  min(w(ilp^s,1:nwflux),ful(il^s,1:nwflux),flc(il^s,1:nwflux)))
381 
382  fmax(il^s,1:nwflux) = min(max(w(ilp^s,1:nwflux),w(il^s,1:nwflux),fmd(il^s,1:nwflux)),&
383  max(w(ilp^s,1:nwflux),ful(il^s,1:nwflux),flc(il^s,1:nwflux)))
384 
385  do iw=1,nwflux
386  a(il^s) = fmin(il^s,iw)
387  b(il^s) = f(il^s,iw)
388  c(il^s) = fmax(il^s,iw)
389  call median(ixi^l,il^l,a,b,c,tmp)
390  flim(il^s,iw) = tmp(il^s)
391  end do
392 
393  ! check case
394  where ((f(il^s,1:nwflux)-w(ilp^s,1:nwflux))*(f(il^s,1:nwflux)-fmp(il^s,1:nwflux)) .le. eps)
395  wrc(il^s,1:nwflux) = f(il^s,1:nwflux)
396  elsewhere
397  wrc(il^s,1:nwflux) = flim(il^s,1:nwflux)
398  end where
399 
400  !call fix_onelimiter1(ixI^L,iL^L,wRCtmp,wRC)
401 
402  end subroutine mp5limiterr
403 
404  subroutine minmod(ixI^L,ixO^L,a,b,minm)
405 
406  use mod_global_parameters
407 
408  integer, intent(in) :: ixI^L, ixO^L
409  double precision, intent(in) :: a(ixI^S), b(ixI^S)
410  double precision, intent(out):: minm(ixI^S)
411 
412  minm(ixo^s) = (sign(one,a(ixo^s))+sign(one,b(ixo^s)))/2.0d0 * &
413  min(abs(a(ixo^s)),abs(b(ixo^s)))
414 
415  end subroutine minmod
416 
417  subroutine median(ixI^L,ixO^L,a,b,c,med)
418 
419  use mod_global_parameters
420 
421  integer, intent(in) :: ixI^L, ixO^L
422  double precision, intent(in) :: a(ixI^S), b(ixI^S), c(ixI^S)
423  double precision, intent(out):: med(ixI^S)
424 
425  double precision :: tmp1(ixI^S),tmp2(ixI^S)
426 
427  tmp1(ixo^s) = b(ixo^s) - a(ixo^s); tmp2(ixo^s) = c(ixo^s) - a(ixo^s)
428 
429  med(ixo^s) = a(ixo^s) + (sign(one,tmp1(ixo^s))+sign(one,tmp2(ixo^s)))/2.0d0 * &
430  min(abs(tmp1(ixo^s)),abs(tmp2(ixo^s)))
431 
432  end subroutine median
433 
434  !> MP5 limiter from Suresh & Huynh 1997
435  !> Following the convention of Mignone et al. 2010.
436  !> Needs at least three ghost cells. Set nghostcells=3.
437  ! for one variable only: no fixing applied
438  subroutine mp5limitervar(ixI^L,iL^L,idims,w,wLC,wRC)
439  use mod_global_parameters
440 
441  integer, intent(in) :: ixi^l, il^l, idims
442  double precision, intent(in) :: w(ixi^s)
443  double precision, intent(inout) :: wrc(ixi^s),wlc(ixi^s)
444 
445  double precision, dimension(ixI^S) :: f, fmp, fmin, fmax, ful, dm4, d, fmd, flc, flim
446  double precision, dimension(ixI^S) :: wrctmp, wlctmp
447  double precision, dimension(ixI^S) :: tmp, tmp2, tmp3, a, b, c
448  double precision, parameter :: eps=0.0d0, alpha=4.0d0
449  !double precision :: alpha
450  integer :: ilm^l, ilmm^l, ilp^l, ilpp^l, ilppp^l
451  integer :: id^l, idp^l, idpp^l, idm^l, ie^l, iem^l, iep^l, iepp^l
452 
453  ! Variable alpha:
454  !alpha = float(nstep)/courantpar - one
455 
456  ! Left side:
457  ! range to process:
458  !iLmin^D=ixmin^D-kr(idims,^D);iLmax^D=ixmax^D;
459 
460  ! HALL
461  ! For Hall, we need one more reconstructed layer since currents are computed in getflux:
462  ! also add one ghost zone!
463  ! {iL^L=iL^L^LADD1;}
464 
465  ! iL^L holds the indices of interfaces to reconstruct to. Convention is that a center index holds the _right-side_ interface.
466 
467  ilm^l=il^l-kr(idims,^d);
468  ilmm^l=ilm^l-kr(idims,^d);
469  ilp^l=il^l+kr(idims,^d);
470  ilpp^l=ilp^l+kr(idims,^d);
471 
472  f(il^s) = 1.0d0/60.0d0 * (&
473  2.0d0* w(ilmm^s) &
474  - 13.0d0* w(ilm^s) &
475  + 47.0d0* w(il^s) &
476  + 27.0d0* w(ilp^s) &
477  - 3.0d0* w(ilpp^s))
478 
479  ! get fmp and ful:
480  a(il^s) = w(ilp^s)-w(il^s)
481  b(il^s) = alpha*(w(il^s)-w(ilm^s))
482  call minmod(ixi^l,il^l,a,b,tmp)
483  fmp(il^s) = w(il^s) + tmp(il^s)
484  ful(il^s) = w(il^s) + b(il^s)
485 
486  ! get dm4:
487  idmax^d=ilmax^d; idmin^d=ilmin^d-kr(idims,^d);
488  idm^l=id^l-kr(idims,^d);
489  idp^l=id^l+kr(idims,^d);
490 
491  iemax^d=idmax^d+kr(idims,^d); iemin^d=idmin^d;
492  iem^l=ie^l-kr(idims,^d);
493  iep^l=ie^l+kr(idims,^d);
494 
495  d(ie^s) = w(iep^s)-2.0d0*w(ie^s)+w(iem^s)
496 
497  a(id^s) = 4.0d0*d(id^s)-d(idp^s)
498  b(id^s) = 4.0d0*d(idp^s)-d(id^s)
499  call minmod(ixi^l,id^l,a,b,tmp)
500  a(id^s) = d(id^s)
501  b(id^s) = d(idp^s)
502  call minmod(ixi^l,id^l,a,b,tmp2)
503  call minmod(ixi^l,id^l,tmp,tmp2,tmp3)
504  dm4(id^s) = tmp3(id^s)
505 
506  ! get fmd:
507  fmd(il^s) = (w(il^s)+w(ilp^s))/2.0d0-dm4(il^s)/2.0d0
508 
509  !get flc:
510  flc(il^s) = half*(3.0d0*w(il^s) &
511  - w(ilm^s)) + 4.0d0/3.0d0*dm4(ilm^s)
512 
513  fmin(il^s) = max(min(w(il^s),w(ilp^s),fmd(il^s)),&
514  min(w(il^s),ful(il^s),flc(il^s)))
515 
516  fmax(il^s) = min(max(w(il^s),w(ilp^s),fmd(il^s)),&
517  max(w(il^s),ful(il^s),flc(il^s)))
518 
519  call median(ixi^l,il^l,fmin,f,fmax,tmp)
520  flim(il^s) = tmp(il^s)
521 
522  ! check case
523  where ((f(il^s)-w(il^s))*(f(il^s)-fmp(il^s)) .le. eps)
524  wlc(il^s) = f(il^s)
525  elsewhere
526  wlc(il^s) = flim(il^s)
527  end where
528 
529  ! Right side:
530  ! the interpolation from the right is obtained when the left-hand formula is applied to
531  ! data mirrored about the interface.
532  ! thus substitute:
533  ! i-2 -> i+3
534  ! i-1 -> i+2
535  ! i -> i+1
536  ! i+1 -> i
537  ! i+2 -> i-1
538 
539  ilppp^l=ilpp^l+kr(idims,^d);
540 
541  f(il^s) = 1.0d0/60.0d0 * (&
542  2.0d0* w(ilppp^s) &
543  - 13.0d0* w(ilpp^s) &
544  + 47.0d0* w(ilp^s) &
545  + 27.0d0* w(il^s) &
546  - 3.0d0* w(ilm^s))
547 
548  ! get fmp and ful:
549  a(il^s) = w(il^s)-w(ilp^s)
550  b(il^s) = alpha*(w(ilp^s)-w(ilpp^s))
551  call minmod(ixi^l,il^l,a,b,tmp)
552  fmp(il^s) = w(ilp^s) + tmp(il^s)
553  ful(il^s) = w(ilp^s) + b(il^s)
554 
555  ! get dm4:
556  idmax^d=ilmax^d+kr(idims,^d); idmin^d=ilmin^d;
557  idm^l=id^l-kr(idims,^d);
558  idp^l=id^l+kr(idims,^d);
559 
560  iemax^d=idmax^d; iemin^d=idmin^d-kr(idims,^d);
561  iem^l=ie^l-kr(idims,^d);
562  iep^l=ie^l+kr(idims,^d);
563  iepp^l=iep^l+kr(idims,^d);
564 
565  d(ie^s) = w(ie^s)-2.0d0*w(iep^s)+w(iepp^s)
566 
567  a(id^s) = 4.0d0*d(id^s)-d(idm^s)
568  b(id^s) = 4.0d0*d(idm^s)-d(id^s)
569  call minmod(ixi^l,id^l,a,b,tmp)
570  a(id^s) = d(id^s)
571  b(id^s) = d(idm^s)
572  call minmod(ixi^l,id^l,a,b,tmp2)
573  call minmod(ixi^l,id^l,tmp,tmp2,tmp3)
574  dm4(id^s) = tmp3(id^s)
575 
576  ! get fmd:
577  fmd(il^s) = (w(il^s)+w(ilp^s))/2.0d0-dm4(il^s)/2.0d0
578 
579  !get flc:
580  flc(il^s) = half*(3.0d0*w(ilp^s) &
581  - w(ilpp^s)) + 4.0d0/3.0d0*dm4(ilp^s)
582 
583  fmin(il^s) = max(min(w(ilp^s),w(il^s),fmd(il^s)),&
584  min(w(ilp^s),ful(il^s),flc(il^s)))
585 
586  fmax(il^s) = min(max(w(ilp^s),w(il^s),fmd(il^s)),&
587  max(w(ilp^s),ful(il^s),flc(il^s)))
588 
589  call median(ixi^l,il^l,fmin,f,fmax,flim)
590 
591  ! check case
592  where ((f(il^s)-w(ilp^s))*(f(il^s)-fmp(il^s)) .le. eps)
593  wrc(il^s) = f(il^s)
594  elsewhere
595  wrc(il^s) = flim(il^s)
596  end where
597 
598  end subroutine mp5limitervar
599 
600 end module mod_mp5
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::kr
integer, dimension(3, 3) kr
Kronecker delta tensor.
Definition: mod_global_parameters.t:426
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_mp5
Module containing the MP5 (fifth order) flux scheme.
Definition: mod_mp5.t:2
mod_mp5::mp5limiter
subroutine, public mp5limiter(ixIL, iLL, idims, w, wLC, wRC)
MP5 limiter from Suresh & Huynh 1997 Following the convention of Mignone et al. 2010....
Definition: mod_mp5.t:17
mod_mp5::median
subroutine median(ixIL, ixOL, a, b, c, med)
Definition: mod_mp5.t:418
mod_mp5::minmod
subroutine minmod(ixIL, ixOL, a, b, minm)
Definition: mod_mp5.t:405
mod_mp5::mp5limitervar
subroutine, public mp5limitervar(ixIL, iLL, idims, w, wLC, wRC)
MP5 limiter from Suresh & Huynh 1997 Following the convention of Mignone et al. 2010....
Definition: mod_mp5.t:439
mod_mp5::mp5limiterr
subroutine, public mp5limiterr(ixIL, iLL, idims, w, wRC)
Definition: mod_mp5.t:301
mod_mp5::mp5limiterl
subroutine, public mp5limiterl(ixIL, iLL, idims, w, wLC)
Definition: mod_mp5.t:204
mod_weno
Definition: mod_weno.t:1
mod_weno::fix_limiter1
subroutine, public fix_limiter1(ixIL, iLL, wLCin, wRCin, wLCout, wRCout)
Definition: mod_weno.t:136
Generated by doxygen 1.9.1