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