MPI-AMRVAC  2.2
The MPI - Adaptive Mesh Refinement - Versatile Advection Code
mod_limiter.t
Go to the documentation of this file.
1 !> Module with slope/flux limiters
2 module mod_limiter
3  use mod_ppm
4  use mod_mp5
5  use mod_weno
6  use mod_venk
7 
8  implicit none
9  public
10 
11  !> radius of the asymptotic region [0.001, 10], larger means more accurate in smooth
12  !> region but more overshooting at discontinuities
13  double precision :: cada3_radius
14  double precision :: schmid_rad^d
15  integer, parameter :: limiter_minmod = 1
16  integer, parameter :: limiter_woodward = 2
17  integer, parameter :: limiter_mcbeta = 3
18  integer, parameter :: limiter_superbee = 4
19  integer, parameter :: limiter_vanleer = 5
20  integer, parameter :: limiter_albada = 6
21  integer, parameter :: limiter_koren = 7
22  integer, parameter :: limiter_cada = 8
23  integer, parameter :: limiter_cada3 = 9
24  integer, parameter :: limiter_schmid = 10
25  integer, parameter :: limiter_venk = 11
26  ! Special cases
27  integer, parameter :: limiter_ppm = 12
28  integer, parameter :: limiter_mp5 = 13
29  integer, parameter :: limiter_weno3 = 14
30  integer, parameter :: limiter_wenoyc3 = 15
31  integer, parameter :: limiter_weno5 = 16
32  integer, parameter :: limiter_weno5nm = 17
33  integer, parameter :: limiter_wenoz5 = 18
34  integer, parameter :: limiter_wenoz5nm = 19
35  integer, parameter :: limiter_wenozp5 = 20
36  integer, parameter :: limiter_wenozp5nm = 21
37  integer, parameter :: limiter_weno5cu6 = 22
38  integer, parameter :: limiter_weno7 = 23
39  integer, parameter :: limiter_mpweno7 = 24
40  integer, parameter :: limiter_exeno7 = 25
41 
42 contains
43 
44  integer function limiter_type(namelim)
45  character(len=*), intent(in) :: namelim
46 
47  select case (namelim)
48  case ('minmod')
50  case ('woodward')
52  case ('mcbeta')
54  case ('superbee')
56  case ('vanleer')
58  case ('albada')
60  case ('koren')
62  case ('cada')
64  case ('cada3')
66  case ('schmid1')
68  case ('schmid2')
70  case('venk')
72  case ('ppm')
74  case ('mp5')
76  case ('weno3')
78  case ('wenoyc3')
80  case ('weno5')
82  case ('weno5nm')
84  case ('wenoz5')
86  case ('wenoz5nm')
88  case ('wenozp5')
90  case ('wenozp5nm')
92  case ('weno5cu6')
94  case ('weno7')
96  case ('mpweno7')
98  case ('exeno7')
100 
101  case default
102  limiter_type = -1
103  write(*,*) 'Unknown limiter: ', namelim
104  call mpistop("No such limiter")
105  end select
106  end function limiter_type
107 
108  pure logical function limiter_symmetric(typelim)
109  integer, intent(in) :: typelim
110 
111  select case (typelim)
113  limiter_symmetric = .false.
114  case default
115  limiter_symmetric = .true.
116  end select
117  end function limiter_symmetric
118 
119  !> Limit the centered dwC differences within ixC for iw in direction idim.
120  !> The limiter is chosen according to typelimiter.
121  !>
122  !> Note that this subroutine is called from upwindLR (hence from methods
123  !> like tvdlf, hancock, hll(c) etc) or directly from tvd.t,
124  !> but also from the gradientS and divvectorS subroutines in geometry.t
125  !> Accordingly, the typelimiter here corresponds to one of limiter
126  !> or one of gradient_limiter.
127  subroutine dwlimiter2(dwC,ixI^L,ixC^L,idims,typelim,ldw,rdw,a2max)
130 
131  integer, intent(in) :: ixI^L, ixC^L, idims
132  double precision, intent(in) :: dwC(ixi^s)
133  integer, intent(in) :: typelim
134  !> Result using left-limiter (same as right for symmetric)
135  double precision, intent(out), optional :: ldw(ixi^s)
136  !> Result using right-limiter (same as left for symmetric)
137  double precision, intent(out), optional :: rdw(ixi^s)
138  double precision, intent(in), optional :: a2max
139 
140  double precision :: tmp(ixi^s), tmp2(ixi^s)
141  integer :: ixO^L, hxO^L
142  double precision, parameter :: qsmall=1.d-12, qsmall2=2.d-12
143  double precision, parameter :: eps = sqrt(epsilon(1.0d0))
144 
145  ! mcbeta limiter parameter value
146  double precision, parameter :: c_mcbeta=1.4d0
147  ! cada limiter parameter values
148  double precision, parameter :: cadalfa=0.5d0, cadbeta=2.0d0, cadgamma=1.6d0
149  ! full third order cada limiter
150  double precision :: rdelinv
151  double precision :: ldwA(ixi^s),ldwB(ixi^s),tmpeta(ixi^s)
152  double precision, parameter :: cadepsilon=1.d-14, invcadepsilon=1.d14,cada3_radius=0.1d0
153  integer :: ix^D
154  !-----------------------------------------------------------------------------
155 
156  ! Contract indices in idim for output.
157  ixomin^d=ixcmin^d+kr(idims,^d); ixomax^d=ixcmax^d;
158  hxo^l=ixo^l-kr(idims,^d);
159 
160  ! About the notation: the conventional argument theta (the ratio of slopes)
161  ! would be given by dwC(ixO^S)/dwC(hxO^S). However, in the end one
162  ! multiplies phi(theta) by dwC(hxO^S), which is incorporated in the
163  ! equations below. The minmod limiter can for example be written as:
164  ! A:
165  ! max(0.0d0, min(1.0d0, dwC(ixO^S)/dwC(hxO^S))) * dwC(hxO^S)
166  ! B:
167  ! tmp(ixO^S)*max(0.0d0,min(abs(dwC(ixO^S)),tmp(ixO^S)*dwC(hxO^S)))
168  ! where tmp(ixO^S)=sign(1.0d0,dwC(ixO^S))
169 
170  select case (typelim)
171  case (limiter_minmod)
172  ! Minmod limiter eq(3.51e) and (eq.3.38e) with omega=1
173  tmp(ixo^s)=sign(one,dwc(ixo^s))
174  tmp(ixo^s)=tmp(ixo^s)* &
175  max(zero,min(abs(dwc(ixo^s)),tmp(ixo^s)*dwc(hxo^s)))
176  if (present(ldw)) ldw(ixo^s) = tmp(ixo^s)
177  if (present(rdw)) rdw(ixo^s) = tmp(ixo^s)
178  case (limiter_woodward)
179  ! Woodward and Collela limiter (eq.3.51h), a factor of 2 is pulled out
180  tmp(ixo^s)=sign(one,dwc(ixo^s))
181  tmp(ixo^s)=2*tmp(ixo^s)* &
182  max(zero,min(abs(dwc(ixo^s)),tmp(ixo^s)*dwc(hxo^s),&
183  tmp(ixo^s)*quarter*(dwc(hxo^s)+dwc(ixo^s))))
184  if (present(ldw)) ldw(ixo^s) = tmp(ixo^s)
185  if (present(rdw)) rdw(ixo^s) = tmp(ixo^s)
186  case (limiter_mcbeta)
187  ! Woodward and Collela limiter, with factor beta
188  tmp(ixo^s)=sign(one,dwc(ixo^s))
189  tmp(ixo^s)=tmp(ixo^s)* &
190  max(zero,min(c_mcbeta*abs(dwc(ixo^s)),c_mcbeta*tmp(ixo^s)*dwc(hxo^s),&
191  tmp(ixo^s)*half*(dwc(hxo^s)+dwc(ixo^s))))
192  if (present(ldw)) ldw(ixo^s) = tmp(ixo^s)
193  if (present(rdw)) rdw(ixo^s) = tmp(ixo^s)
194  case (limiter_superbee)
195  ! Roes superbee limiter (eq.3.51i)
196  tmp(ixo^s)=sign(one,dwc(ixo^s))
197  tmp(ixo^s)=tmp(ixo^s)* &
198  max(zero,min(2*abs(dwc(ixo^s)),tmp(ixo^s)*dwc(hxo^s)),&
199  min(abs(dwc(ixo^s)),2*tmp(ixo^s)*dwc(hxo^s)))
200  if (present(ldw)) ldw(ixo^s) = tmp(ixo^s)
201  if (present(rdw)) rdw(ixo^s) = tmp(ixo^s)
202  case (limiter_vanleer)
203  ! van Leer limiter (eq 3.51f), but a missing delta2=1.D-12 is added
204  tmp(ixo^s)=2*max(dwc(hxo^s)*dwc(ixo^s),zero) &
205  /(dwc(ixo^s)+dwc(hxo^s)+qsmall)
206  if (present(ldw)) ldw(ixo^s) = tmp(ixo^s)
207  if (present(rdw)) rdw(ixo^s) = tmp(ixo^s)
208  case (limiter_albada)
209  ! Albada limiter (eq.3.51g) with delta2=1D.-12
210  tmp(ixo^s)=(dwc(hxo^s)*(dwc(ixo^s)**2+qsmall)&
211  +dwc(ixo^s)*(dwc(hxo^s)**2+qsmall))&
212  /(dwc(ixo^s)**2+dwc(hxo^s)**2+qsmall2)
213  if (present(ldw)) ldw(ixo^s) = tmp(ixo^s)
214  if (present(rdw)) rdw(ixo^s) = tmp(ixo^s)
215  case (limiter_koren)
216  tmp(ixo^s)=sign(one,dwc(ixo^s))
217  tmp2(ixo^s)=min(2*abs(dwc(ixo^s)),2*tmp(ixo^s)*dwc(hxo^s))
218  if (present(ldw)) then
219  ldw(ixo^s)=tmp(ixo^s)* &
220  max(zero,min(tmp2(ixo^s),&
221  (dwc(hxo^s)*tmp(ixo^s)+2*abs(dwc(ixo^s)))*third))
222  end if
223  if (present(rdw)) then
224  rdw(ixo^s)=tmp(ixo^s)* &
225  max(zero,min(tmp2(ixo^s),&
226  (2*dwc(hxo^s)*tmp(ixo^s)+abs(dwc(ixo^s)))*third))
227  end if
228  case (limiter_cada)
229  ! This limiter has been rewritten in the usual form, and uses a division
230  ! of the gradients.
231  if (present(ldw)) then
232  ! Cada Left variant
233  ! Compute theta, but avoid division by zero
234  tmp(ixo^s)=dwc(hxo^s)/(dwc(ixo^s) + sign(eps, dwc(ixo^s)))
235  tmp2(ixo^s)=(2+tmp(ixo^s))*third
236  ldw(ixo^s)= max(zero,min(tmp2(ixo^s), &
237  max(-cadalfa*tmp(ixo^s), &
238  min(cadbeta*tmp(ixo^s), tmp2(ixo^s), &
239  cadgamma)))) * dwc(ixo^s)
240  end if
241 
242  if (present(rdw)) then
243  ! Cada Right variant
244  tmp(ixo^s)=dwc(ixo^s)/(dwc(hxo^s) + sign(eps, dwc(hxo^s)))
245  tmp2(ixo^s)=(2+tmp(ixo^s))*third
246  rdw(ixo^s)= max(zero,min(tmp2(ixo^s), &
247  max(-cadalfa*tmp(ixo^s), &
248  min(cadbeta*tmp(ixo^s), tmp2(ixo^s), &
249  cadgamma)))) * dwc(hxo^s)
250  end if
251  case (limiter_cada3)
252  rdelinv=one/(cada3_radius*dxlevel(idims))**2
253  tmpeta(ixo^s)=(dwc(ixo^s)**2+dwc(hxo^s)**2)*rdelinv
254  if (present(ldw)) then
255  tmp(ixo^s)=dwc(hxo^s)/(dwc(ixo^s) + sign(eps, dwc(ixo^s)))
256  ldwa(ixo^s)=(two+tmp(ixo^s))*third
257  where(tmpeta(ixo^s)<=one-cadepsilon)
258  ldw(ixo^s)=ldwa(ixo^s)
259  elsewhere(tmpeta(ixo^s)>=one+cadepsilon)
260  ldwb(ixo^s)= max(zero,min(ldwa(ixo^s), max(-cadalfa*tmp(ixo^s), &
261  min(cadbeta*tmp(ixo^s), ldwa(ixo^s), cadgamma))))
262  ldw(ixo^s)=ldwb(ixo^s)
263  elsewhere
264  ldwb(ixo^s)= max(zero,min(ldwa(ixo^s), max(-cadalfa*tmp(ixo^s), &
265  min(cadbeta*tmp(ixo^s), ldwa(ixo^s), cadgamma))))
266  tmp2(ixo^s)=(tmpeta(ixo^s)-one)*invcadepsilon
267  ldw(ixo^s)=half*( (one-tmp2(ixo^s))*ldwa(ixo^s) &
268  +(one+tmp2(ixo^s))*ldwb(ixo^s))
269  endwhere
270  ldw(ixo^s)=ldw(ixo^s) * dwc(ixo^s)
271  end if
272 
273  if (present(rdw)) then
274  tmp(ixo^s)=dwc(ixo^s)/(dwc(hxo^s) + sign(eps, dwc(hxo^s)))
275  ldwa(ixo^s)=(two+tmp(ixo^s))*third
276  where(tmpeta(ixo^s)<=one-cadepsilon)
277  rdw(ixo^s)=ldwa(ixo^s)
278  elsewhere(tmpeta(ixo^s)>=one+cadepsilon)
279  ldwb(ixo^s)= max(zero,min(ldwa(ixo^s), max(-cadalfa*tmp(ixo^s), &
280  min(cadbeta*tmp(ixo^s), ldwa(ixo^s), cadgamma))))
281  rdw(ixo^s)=ldwb(ixo^s)
282  elsewhere
283  ldwb(ixo^s)= max(zero,min(ldwa(ixo^s), max(-cadalfa*tmp(ixo^s), &
284  min(cadbeta*tmp(ixo^s), ldwa(ixo^s), cadgamma))))
285  tmp2(ixo^s)=(tmpeta(ixo^s)-one)*invcadepsilon
286  rdw(ixo^s)=half*( (one-tmp2(ixo^s))*ldwa(ixo^s) &
287  +(one+tmp2(ixo^s))*ldwb(ixo^s))
288  endwhere
289  rdw(ixo^s)=rdw(ixo^s) * dwc(hxo^s)
290  end if
291  case(limiter_schmid)
292  tmpeta(ixo^s)=(sqrt(0.4d0*(dwc(ixo^s)**2+dwc(hxo^s)**2)))&
293  /((a2max+cadepsilon)*dxlevel(idims)**2)
294  if(present(ldw)) then
295  tmp(ixo^s)=dwc(hxo^s)/(dwc(ixo^s)+sign(eps,dwc(ixo^s)))
296  ldwa(ixo^s)=(two+tmp(ixo^s))*third
297  where(tmpeta(ixo^s)<=one-cadepsilon)
298  ldw(ixo^s)=ldwa(ixo^s)
299  else where(tmpeta(ixo^s)>=one+cadepsilon)
300  ldwb(ixo^s)=max(zero,min(ldwa(ixo^s),max(-tmp(ixo^s),&
301  min(cadbeta*tmp(ixo^s),ldwa(ixo^s),1.5d0))))
302  ldw(ixo^s)=ldwb(ixo^s)
303  else where
304  ldwb(ixo^s)=max(zero,min(ldwa(ixo^s),max(-tmp(ixo^s),&
305  min(cadbeta*tmp(ixo^s),ldwa(ixo^s),1.5d0))))
306  tmp2(ixo^s)=(tmpeta(ixo^s)-one)*invcadepsilon
307  ldw(ixo^s)=half*((one-tmp2(ixo^s))*ldwa(ixo^s)&
308  +(one+tmp2(ixo^s))*ldwb(ixo^s))
309  end where
310  ldw(ixo^s)=ldw(ixo^s)*dwc(ixo^s)
311  end if
312  if(present(rdw)) then
313  tmp(ixo^s)=dwc(ixo^s)/(dwc(hxo^s)+sign(eps,dwc(hxo^s)))
314  ldwa(ixo^s)=(two+tmp(ixo^s))*third
315  where(tmpeta(ixo^s)<=one-cadepsilon)
316  rdw(ixo^s)=ldwa(ixo^s)
317  else where(tmpeta(ixo^s)>=one+cadepsilon)
318  ldwb(ixo^s)=max(zero,min(ldwa(ixo^s),max(-tmp(ixo^s),&
319  min(cadbeta*tmp(ixo^s),ldwa(ixo^s),1.5d0))))
320  rdw(ixo^s)=ldwb(ixo^s)
321  else where
322  ldwb(ixo^s)=max(zero,min(ldwa(ixo^s),max(-tmp(ixo^s),&
323  min(cadbeta*tmp(ixo^s), ldwa(ixo^s), 1.5d0))))
324  tmp2(ixo^s)=(tmpeta(ixo^s)-one)*invcadepsilon
325  rdw(ixo^s)=half*((one-tmp2(ixo^s))*ldwa(ixo^s)&
326  +(one+tmp2(ixo^s))*ldwb(ixo^s))
327  end where
328  rdw(ixo^s)=rdw(ixo^s)*dwc(hxo^s)
329  end if
330  case default
331  call mpistop("Error in dwLimiter: unknown limiter")
332  end select
333 
334  end subroutine dwlimiter2
335 
336 end module mod_limiter
This module contains definitions of global parameters and variables and some generic functions/subrou...
pure logical function limiter_symmetric(typelim)
Definition: mod_limiter.t:109
integer, parameter limiter_superbee
Definition: mod_limiter.t:18
integer, parameter limiter_cada3
Definition: mod_limiter.t:23
integer, parameter limiter_wenoz5
Definition: mod_limiter.t:33
integer, parameter limiter_wenozp5
Definition: mod_limiter.t:35
integer, parameter limiter_mp5
Definition: mod_limiter.t:28
integer, parameter limiter_albada
Definition: mod_limiter.t:20
integer, parameter limiter_wenoyc3
Definition: mod_limiter.t:30
integer, parameter limiter_mpweno7
Definition: mod_limiter.t:39
Module containing the MP5 (fifth order) flux scheme.
Definition: mod_mp5.t:2
integer, parameter limiter_woodward
Definition: mod_limiter.t:16
integer, parameter limiter_weno7
Definition: mod_limiter.t:38
integer, parameter limiter_ppm
Definition: mod_limiter.t:27
Module with slope/flux limiters.
Definition: mod_limiter.t:2
integer, parameter limiter_schmid
Definition: mod_limiter.t:24
integer, parameter limiter_wenozp5nm
Definition: mod_limiter.t:36
integer, parameter limiter_weno5nm
Definition: mod_limiter.t:32
integer function limiter_type(namelim)
Definition: mod_limiter.t:45
integer, parameter limiter_koren
Definition: mod_limiter.t:21
double precision cada3_radius
radius of the asymptotic region [0.001, 10], larger means more accurate in smooth region but more ove...
Definition: mod_limiter.t:13
integer, dimension(3, 3) kr
Kronecker delta tensor.
integer, dimension(:), allocatable, parameter d
integer, parameter limiter_weno5cu6
Definition: mod_limiter.t:37
integer, parameter limiter_mcbeta
Definition: mod_limiter.t:17
integer, parameter limiter_venk
Definition: mod_limiter.t:25
Definition: mod_ppm.t:1
subroutine dwlimiter2(dwC, ixIL, ixCL, idims, typelim, ldw, rdw, a2max)
Limit the centered dwC differences within ixC for iw in direction idim. The limiter is chosen accordi...
Definition: mod_limiter.t:128
integer, parameter limiter_exeno7
Definition: mod_limiter.t:40
subroutine mpistop(message)
Exit MPI-AMRVAC with an error message.
Definition: comm_lib.t:198
integer, parameter limiter_weno5
Definition: mod_limiter.t:31
double precision, dimension(ndim) dxlevel
double precision schmid_rad
Definition: mod_limiter.t:14
integer, parameter limiter_vanleer
Definition: mod_limiter.t:19
integer, parameter limiter_minmod
Definition: mod_limiter.t:15
integer, parameter limiter_cada
Definition: mod_limiter.t:22
integer, parameter limiter_weno3
Definition: mod_limiter.t:29
integer, parameter limiter_wenoz5nm
Definition: mod_limiter.t:34
double precision d
Definition: mod_limiter.t:14