MPI-AMRVAC 3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
Loading...
Searching...
No Matches
mod_mp5.t
Go to the documentation of this file.
1!> Module containing the MP5 (fifth order) flux scheme
2module mod_mp5
3
4 implicit none
5 private
6
7 public :: mp5limiter
8 public :: mp5limitervar
9 public :: mp5limiterl
10 public :: mp5limiterr
11
12contains
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:nw_recon) = 1.0d0/60.0d0 * (&
55 2.0d0* w(ilmm^s,1:nw_recon) &
56 - 13.0d0* w(ilm^s,1:nw_recon) &
57 + 47.0d0* w(il^s,1:nw_recon) &
58 + 27.0d0* w(ilp^s,1:nw_recon) &
59 - 3.0d0* w(ilpp^s,1:nw_recon))
60
61 ! get fmp and ful:
62 do iw=1,nw_recon
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:nw_recon) = w(iep^s,1:nw_recon)-2.0d0*w(ie^s,1:nw_recon)+w(iem^s,1:nw_recon)
80
81 do iw=1,nw_recon
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:nw_recon) = (w(il^s,1:nw_recon)+w(ilp^s,1:nw_recon))/2.0d0-dm4(il^s,1:nw_recon)/2.0d0
94
95 !get flc:
96 flc(il^s,1:nw_recon) = half*(3.0d0*w(il^s,1:nw_recon) &
97 - w(ilm^s,1:nw_recon)) + 4.0d0/3.0d0*dm4(ilm^s,1:nw_recon)
98
99 fmin(il^s,1:nw_recon) = max(min(w(il^s,1:nw_recon),w(ilp^s,1:nw_recon),fmd(il^s,1:nw_recon)),&
100 min(w(il^s,1:nw_recon),ful(il^s,1:nw_recon),flc(il^s,1:nw_recon)))
101
102 fmax(il^s,1:nw_recon) = min(max(w(il^s,1:nw_recon),w(ilp^s,1:nw_recon),fmd(il^s,1:nw_recon)),&
103 max(w(il^s,1:nw_recon),ful(il^s,1:nw_recon),flc(il^s,1:nw_recon)))
104
105 do iw=1,nw_recon
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:nw_recon)-w(il^s,1:nw_recon))*(f(il^s,1:nw_recon)-fmp(il^s,1:nw_recon)) .le. eps)
115 wlctmp(il^s,1:nw_recon) = f(il^s,1:nw_recon)
116 elsewhere
117 wlctmp(il^s,1:nw_recon) = flim(il^s,1:nw_recon)
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:nw_recon) = 1.0d0/60.0d0 * (&
133 2.0d0* w(ilppp^s,1:nw_recon) &
134 - 13.0d0* w(ilpp^s,1:nw_recon) &
135 + 47.0d0* w(ilp^s,1:nw_recon) &
136 + 27.0d0* w(il^s,1:nw_recon) &
137 - 3.0d0* w(ilm^s,1:nw_recon))
138
139 ! get fmp and ful:
140 do iw=1,nw_recon
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:nw_recon) = w(ie^s,1:nw_recon)-2.0d0*w(iep^s,1:nw_recon)+w(iepp^s,1:nw_recon)
159
160 do iw=1,nw_recon
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:nw_recon) = (w(il^s,1:nw_recon)+w(ilp^s,1:nw_recon))/2.0d0-dm4(il^s,1:nw_recon)/2.0d0
173
174 !get flc:
175 flc(il^s,1:nw_recon) = half*(3.0d0*w(ilp^s,1:nw_recon) &
176 - w(ilpp^s,1:nw_recon)) + 4.0d0/3.0d0*dm4(ilp^s,1:nw_recon)
177
178 fmin(il^s,1:nw_recon) = max(min(w(ilp^s,1:nw_recon),w(il^s,1:nw_recon),fmd(il^s,1:nw_recon)),&
179 min(w(ilp^s,1:nw_recon),ful(il^s,1:nw_recon),flc(il^s,1:nw_recon)))
180
181 fmax(il^s,1:nw_recon) = min(max(w(ilp^s,1:nw_recon),w(il^s,1:nw_recon),fmd(il^s,1:nw_recon)),&
182 max(w(ilp^s,1:nw_recon),ful(il^s,1:nw_recon),flc(il^s,1:nw_recon)))
183
184 do iw=1,nw_recon
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:nw_recon)-w(ilp^s,1:nw_recon))*(f(il^s,1:nw_recon)-fmp(il^s,1:nw_recon)) .le. eps)
194 wrctmp(il^s,1:nw_recon) = f(il^s,1:nw_recon)
195 elsewhere
196 wrctmp(il^s,1:nw_recon) = flim(il^s,1:nw_recon)
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:nw_recon) = 1.0d0/60.0d0 * (&
231 2.0d0* w(ilmm^s,1:nw_recon) &
232 - 13.0d0* w(ilm^s,1:nw_recon) &
233 + 47.0d0* w(il^s,1:nw_recon) &
234 + 27.0d0* w(ilp^s,1:nw_recon) &
235 - 3.0d0* w(ilpp^s,1:nw_recon))
236
237 ! get fmp and ful:
238 do iw=1,nw_recon
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:nw_recon) = w(iep^s,1:nw_recon)-2.0d0*w(ie^s,1:nw_recon)+w(iem^s,1:nw_recon)
256
257 do iw=1,nw_recon
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:nw_recon) = (w(il^s,1:nw_recon)+w(ilp^s,1:nw_recon))/2.0d0-dm4(il^s,1:nw_recon)/2.0d0
270
271 !get flc:
272 flc(il^s,1:nw_recon) = half*(3.0d0*w(il^s,1:nw_recon) &
273 - w(ilm^s,1:nw_recon)) + 4.0d0/3.0d0*dm4(ilm^s,1:nw_recon)
274
275 fmin(il^s,1:nw_recon) = max(min(w(il^s,1:nw_recon),w(ilp^s,1:nw_recon),fmd(il^s,1:nw_recon)),&
276 min(w(il^s,1:nw_recon),ful(il^s,1:nw_recon),flc(il^s,1:nw_recon)))
277
278 fmax(il^s,1:nw_recon) = min(max(w(il^s,1:nw_recon),w(ilp^s,1:nw_recon),fmd(il^s,1:nw_recon)),&
279 max(w(il^s,1:nw_recon),ful(il^s,1:nw_recon),flc(il^s,1:nw_recon)))
280
281 do iw=1,nw_recon
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:nw_recon)-w(il^s,1:nw_recon))*(f(il^s,1:nw_recon)-fmp(il^s,1:nw_recon)) .le. eps)
291 wlc(il^s,1:nw_recon) = f(il^s,1:nw_recon)
292 elsewhere
293 wlc(il^s,1:nw_recon) = flim(il^s,1:nw_recon)
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:nw_recon) = 1.0d0/60.0d0 * (&
334 2.0d0* w(ilppp^s,1:nw_recon) &
335 - 13.0d0* w(ilpp^s,1:nw_recon) &
336 + 47.0d0* w(ilp^s,1:nw_recon) &
337 + 27.0d0* w(il^s,1:nw_recon) &
338 - 3.0d0* w(ilm^s,1:nw_recon))
339
340 ! get fmp and ful:
341 do iw=1,nw_recon
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:nw_recon) = w(ie^s,1:nw_recon)-2.0d0*w(iep^s,1:nw_recon)+w(iepp^s,1:nw_recon)
360
361 do iw=1,nw_recon
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:nw_recon) = (w(il^s,1:nw_recon)+w(ilp^s,1:nw_recon))/2.0d0-dm4(il^s,1:nw_recon)/2.0d0
374
375 !get flc:
376 flc(il^s,1:nw_recon) = half*(3.0d0*w(ilp^s,1:nw_recon) &
377 - w(ilpp^s,1:nw_recon)) + 4.0d0/3.0d0*dm4(ilp^s,1:nw_recon)
378
379 fmin(il^s,1:nw_recon) = max(min(w(ilp^s,1:nw_recon),w(il^s,1:nw_recon),fmd(il^s,1:nw_recon)),&
380 min(w(ilp^s,1:nw_recon),ful(il^s,1:nw_recon),flc(il^s,1:nw_recon)))
381
382 fmax(il^s,1:nw_recon) = min(max(w(ilp^s,1:nw_recon),w(il^s,1:nw_recon),fmd(il^s,1:nw_recon)),&
383 max(w(ilp^s,1:nw_recon),ful(il^s,1:nw_recon),flc(il^s,1:nw_recon)))
384
385 do iw=1,nw_recon
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:nw_recon)-w(ilp^s,1:nw_recon))*(f(il^s,1:nw_recon)-fmp(il^s,1:nw_recon)) .le. eps)
395 wrc(il^s,1:nw_recon) = f(il^s,1:nw_recon)
396 elsewhere
397 wrc(il^s,1:nw_recon) = flim(il^s,1:nw_recon)
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
600end 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::mp5limiterl
subroutine, public mp5limiterl(ixil, ill, idims, w, wlc)
Definition mod_mp5.t:204
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::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_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.8