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 :: smallw(1:nw)
35  !double precision :: alpha
36  !----------------------------------------------------------------------------
37 
38  ! Variable alpha:
39  !alpha = float(nstep)/courantpar - one
40 
41  ! Left side:
42  ! range to process:
43  !iLmin^D=ixmin^D-kr(idims,^D);iLmax^D=ixmax^D;
44 
45  ! HALL
46  ! For Hall, we need one more reconstructed layer since currents are computed in getflux:
47  ! also add one ghost zone!
48  ! {iL^L=iL^L^LADD1;}
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,smallw)
205  call phys_check_w(.true.,ixg^ll,il^l,wrctmp,flagr,smallw)
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  ! 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  ! iL^L holds the indices of interfaces to reconstruct to. Convention is that a center index holds the _right-side_ interface.
475 
476  ilm^l=il^l-kr(idims,^d);
477  ilmm^l=ilm^l-kr(idims,^d);
478  ilp^l=il^l+kr(idims,^d);
479  ilpp^l=ilp^l+kr(idims,^d);
480 
481  f(il^s) = 1.0d0/60.0d0 * (&
482  2.0d0* w(ilmm^s) &
483  - 13.0d0* w(ilm^s) &
484  + 47.0d0* w(il^s) &
485  + 27.0d0* w(ilp^s) &
486  - 3.0d0* w(ilpp^s))
487 
488  ! get fmp and ful:
489  a(il^s) = w(ilp^s)-w(il^s)
490  b(il^s) = alpha*(w(il^s)-w(ilm^s))
491  call minmod(ixi^l,il^l,a,b,tmp)
492  fmp(il^s) = w(il^s) + tmp(il^s)
493  ful(il^s) = w(il^s) + b(il^s)
494 
495  ! get dm4:
496  idmax^d=ilmax^d; idmin^d=ilmin^d-kr(idims,^d);
497  idm^l=id^l-kr(idims,^d);
498  idp^l=id^l+kr(idims,^d);
499 
500  iemax^d=idmax^d+kr(idims,^d); iemin^d=idmin^d;
501  iem^l=ie^l-kr(idims,^d);
502  iep^l=ie^l+kr(idims,^d);
503 
504  d(ie^s) = w(iep^s)-2.0d0*w(ie^s)+w(iem^s)
505 
506  a(id^s) = 4.0d0*d(id^s)-d(idp^s)
507  b(id^s) = 4.0d0*d(idp^s)-d(id^s)
508  call minmod(ixi^l,id^l,a,b,tmp)
509  a(id^s) = d(id^s)
510  b(id^s) = d(idp^s)
511  call minmod(ixi^l,id^l,a,b,tmp2)
512  call minmod(ixi^l,id^l,tmp,tmp2,tmp3)
513  dm4(id^s) = tmp3(id^s)
514 
515  ! get fmd:
516  fmd(il^s) = (w(il^s)+w(ilp^s))/2.0d0-dm4(il^s)/2.0d0
517 
518  !get flc:
519  flc(il^s) = half*(3.0d0*w(il^s) &
520  - w(ilm^s)) + 4.0d0/3.0d0*dm4(ilm^s)
521 
522  fmin(il^s) = max(min(w(il^s),w(ilp^s),fmd(il^s)),&
523  min(w(il^s),ful(il^s),flc(il^s)))
524 
525  fmax(il^s) = min(max(w(il^s),w(ilp^s),fmd(il^s)),&
526  max(w(il^s),ful(il^s),flc(il^s)))
527 
528  call median(ixi^l,il^l,fmin,f,fmax,tmp)
529  flim(il^s) = tmp(il^s)
530 
531  ! check case
532  where ((f(il^s)-w(il^s))*(f(il^s)-fmp(il^s)) .le. eps)
533  wlc(il^s) = f(il^s)
534  elsewhere
535  wlc(il^s) = flim(il^s)
536  end where
537 
538  ! Right side:
539  ! the interpolation from the right is obtained when the left-hand formula is applied to
540  ! data mirrored about the interface.
541  ! thus substitute:
542  ! i-2 -> i+3
543  ! i-1 -> i+2
544  ! i -> i+1
545  ! i+1 -> i
546  ! i+2 -> i-1
547 
548  ilppp^l=ilpp^l+kr(idims,^d);
549 
550  f(il^s) = 1.0d0/60.0d0 * (&
551  2.0d0* w(ilppp^s) &
552  - 13.0d0* w(ilpp^s) &
553  + 47.0d0* w(ilp^s) &
554  + 27.0d0* w(il^s) &
555  - 3.0d0* w(ilm^s))
556 
557  ! get fmp and ful:
558  a(il^s) = w(il^s)-w(ilp^s)
559  b(il^s) = alpha*(w(ilp^s)-w(ilpp^s))
560  call minmod(ixi^l,il^l,a,b,tmp)
561  fmp(il^s) = w(ilp^s) + tmp(il^s)
562  ful(il^s) = w(ilp^s) + b(il^s)
563 
564  ! get dm4:
565  idmax^d=ilmax^d+kr(idims,^d); idmin^d=ilmin^d;
566  idm^l=id^l-kr(idims,^d);
567  idp^l=id^l+kr(idims,^d);
568 
569  iemax^d=idmax^d; iemin^d=idmin^d-kr(idims,^d);
570  iem^l=ie^l-kr(idims,^d);
571  iep^l=ie^l+kr(idims,^d);
572  iepp^l=iep^l+kr(idims,^d);
573 
574  d(ie^s) = w(ie^s)-2.0d0*w(iep^s)+w(iepp^s)
575 
576  a(id^s) = 4.0d0*d(id^s)-d(idm^s)
577  b(id^s) = 4.0d0*d(idm^s)-d(id^s)
578  call minmod(ixi^l,id^l,a,b,tmp)
579  a(id^s) = d(id^s)
580  b(id^s) = d(idm^s)
581  call minmod(ixi^l,id^l,a,b,tmp2)
582  call minmod(ixi^l,id^l,tmp,tmp2,tmp3)
583  dm4(id^s) = tmp3(id^s)
584 
585  ! get fmd:
586  fmd(il^s) = (w(il^s)+w(ilp^s))/2.0d0-dm4(il^s)/2.0d0
587 
588  !get flc:
589  flc(il^s) = half*(3.0d0*w(ilp^s) &
590  - w(ilpp^s)) + 4.0d0/3.0d0*dm4(ilp^s)
591 
592  fmin(il^s) = max(min(w(ilp^s),w(il^s),fmd(il^s)),&
593  min(w(ilp^s),ful(il^s),flc(il^s)))
594 
595  fmax(il^s) = min(max(w(ilp^s),w(il^s),fmd(il^s)),&
596  max(w(ilp^s),ful(il^s),flc(il^s)))
597 
598  call median(ixi^l,il^l,fmin,f,fmax,flim)
599 
600  ! check case
601  where ((f(il^s)-w(ilp^s))*(f(il^s)-fmp(il^s)) .le. eps)
602  wrc(il^s) = f(il^s)
603  elsewhere
604  wrc(il^s) = flim(il^s)
605  end where
606 
607  end subroutine mp5limitervar
608 
609 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:78
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