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