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_teno5ad = 23
39  integer, parameter :: limiter_weno7 = 24
40  integer, parameter :: limiter_mpweno7 = 25
41  integer, parameter :: limiter_exeno7 = 26
42 
43 contains
44 
45  integer function limiter_type(namelim)
46  character(len=*), intent(in) :: namelim
47 
48  select case (namelim)
49  case ('minmod')
51  case ('woodward')
53  case ('mcbeta')
55  case ('superbee')
57  case ('vanleer')
59  case ('albada')
61  case ('koren')
63  case ('cada')
65  case ('cada3')
67  case ('schmid1')
69  case ('schmid2')
71  case('venk')
73  case ('ppm')
75  case ('mp5')
77  case ('weno3')
79  case ('wenoyc3')
81  case ('weno5')
83  case ('weno5nm')
85  case ('wenoz5')
87  case ('wenoz5nm')
89  case ('wenozp5')
91  case ('wenozp5nm')
93  case ('weno5cu6')
95  case ('teno5ad')
97  case ('weno7')
99  case ('mpweno7')
101  case ('exeno7')
103 
104  case default
105  limiter_type = -1
106  write(*,*) 'Unknown limiter: ', namelim
107  call mpistop("No such limiter")
108  end select
109  end function limiter_type
110 
111  pure logical function limiter_symmetric(typelim)
112  integer, intent(in) :: typelim
113 
114  select case (typelim)
116  limiter_symmetric = .false.
117  case default
118  limiter_symmetric = .true.
119  end select
120  end function limiter_symmetric
121 
122  !> Limit the centered dwC differences within ixC for iw in direction idim.
123  !> The limiter is chosen according to typelimiter.
124  !>
125  !> Note that this subroutine is called from upwindLR (hence from methods
126  !> like tvdlf, hancock, hll(c) etc) or directly from tvd.t,
127  !> but also from the gradientS and divvectorS subroutines in geometry.t
128  !> Accordingly, the typelimiter here corresponds to one of limiter
129  !> or one of gradient_limiter.
130  subroutine dwlimiter2(dwC,ixI^L,ixC^L,idims,typelim,ldw,rdw,a2max)
133 
134  integer, intent(in) :: ixI^L, ixC^L, idims
135  double precision, intent(in) :: dwC(ixi^s)
136  integer, intent(in) :: typelim
137  !> Result using left-limiter (same as right for symmetric)
138  double precision, intent(out), optional :: ldw(ixi^s)
139  !> Result using right-limiter (same as left for symmetric)
140  double precision, intent(out), optional :: rdw(ixi^s)
141  double precision, intent(in), optional :: a2max
142 
143  double precision :: tmp(ixi^s), tmp2(ixi^s)
144  integer :: ixO^L, hxO^L
145  double precision, parameter :: qsmall=1.d-12, qsmall2=2.d-12
146  double precision, parameter :: eps = sqrt(epsilon(1.0d0))
147 
148  ! mcbeta limiter parameter value
149  double precision, parameter :: c_mcbeta=1.4d0
150  ! cada limiter parameter values
151  double precision, parameter :: cadalfa=0.5d0, cadbeta=2.0d0, cadgamma=1.6d0
152  ! full third order cada limiter
153  double precision :: rdelinv
154  double precision :: ldwA(ixi^s),ldwB(ixi^s),tmpeta(ixi^s)
155  double precision, parameter :: cadepsilon=1.d-14, invcadepsilon=1.d14,cada3_radius=0.1d0
156  integer :: ix^D
157  !-----------------------------------------------------------------------------
158 
159  ! Contract indices in idim for output.
160  ixomin^d=ixcmin^d+kr(idims,^d); ixomax^d=ixcmax^d;
161  hxo^l=ixo^l-kr(idims,^d);
162 
163  ! About the notation: the conventional argument theta (the ratio of slopes)
164  ! would be given by dwC(ixO^S)/dwC(hxO^S). However, in the end one
165  ! multiplies phi(theta) by dwC(hxO^S), which is incorporated in the
166  ! equations below. The minmod limiter can for example be written as:
167  ! A:
168  ! max(0.0d0, min(1.0d0, dwC(ixO^S)/dwC(hxO^S))) * dwC(hxO^S)
169  ! B:
170  ! tmp(ixO^S)*max(0.0d0,min(abs(dwC(ixO^S)),tmp(ixO^S)*dwC(hxO^S)))
171  ! where tmp(ixO^S)=sign(1.0d0,dwC(ixO^S))
172 
173  select case (typelim)
174  case (limiter_minmod)
175  ! Minmod limiter eq(3.51e) and (eq.3.38e) with omega=1
176  tmp(ixo^s)=sign(one,dwc(ixo^s))
177  tmp(ixo^s)=tmp(ixo^s)* &
178  max(zero,min(abs(dwc(ixo^s)),tmp(ixo^s)*dwc(hxo^s)))
179  if (present(ldw)) ldw(ixo^s) = tmp(ixo^s)
180  if (present(rdw)) rdw(ixo^s) = tmp(ixo^s)
181  case (limiter_woodward)
182  ! Woodward and Collela limiter (eq.3.51h), a factor of 2 is pulled out
183  tmp(ixo^s)=sign(one,dwc(ixo^s))
184  tmp(ixo^s)=2*tmp(ixo^s)* &
185  max(zero,min(abs(dwc(ixo^s)),tmp(ixo^s)*dwc(hxo^s),&
186  tmp(ixo^s)*quarter*(dwc(hxo^s)+dwc(ixo^s))))
187  if (present(ldw)) ldw(ixo^s) = tmp(ixo^s)
188  if (present(rdw)) rdw(ixo^s) = tmp(ixo^s)
189  case (limiter_mcbeta)
190  ! Woodward and Collela limiter, with factor beta
191  tmp(ixo^s)=sign(one,dwc(ixo^s))
192  tmp(ixo^s)=tmp(ixo^s)* &
193  max(zero,min(c_mcbeta*abs(dwc(ixo^s)),c_mcbeta*tmp(ixo^s)*dwc(hxo^s),&
194  tmp(ixo^s)*half*(dwc(hxo^s)+dwc(ixo^s))))
195  if (present(ldw)) ldw(ixo^s) = tmp(ixo^s)
196  if (present(rdw)) rdw(ixo^s) = tmp(ixo^s)
197  case (limiter_superbee)
198  ! Roes superbee limiter (eq.3.51i)
199  tmp(ixo^s)=sign(one,dwc(ixo^s))
200  tmp(ixo^s)=tmp(ixo^s)* &
201  max(zero,min(2*abs(dwc(ixo^s)),tmp(ixo^s)*dwc(hxo^s)),&
202  min(abs(dwc(ixo^s)),2*tmp(ixo^s)*dwc(hxo^s)))
203  if (present(ldw)) ldw(ixo^s) = tmp(ixo^s)
204  if (present(rdw)) rdw(ixo^s) = tmp(ixo^s)
205  case (limiter_vanleer)
206  ! van Leer limiter (eq 3.51f), but a missing delta2=1.D-12 is added
207  tmp(ixo^s)=2*max(dwc(hxo^s)*dwc(ixo^s),zero) &
208  /(dwc(ixo^s)+dwc(hxo^s)+qsmall)
209  if (present(ldw)) ldw(ixo^s) = tmp(ixo^s)
210  if (present(rdw)) rdw(ixo^s) = tmp(ixo^s)
211  case (limiter_albada)
212  ! Albada limiter (eq.3.51g) with delta2=1D.-12
213  tmp(ixo^s)=(dwc(hxo^s)*(dwc(ixo^s)**2+qsmall)&
214  +dwc(ixo^s)*(dwc(hxo^s)**2+qsmall))&
215  /(dwc(ixo^s)**2+dwc(hxo^s)**2+qsmall2)
216  if (present(ldw)) ldw(ixo^s) = tmp(ixo^s)
217  if (present(rdw)) rdw(ixo^s) = tmp(ixo^s)
218  case (limiter_koren)
219  tmp(ixo^s)=sign(one,dwc(ixo^s))
220  tmp2(ixo^s)=min(2*abs(dwc(ixo^s)),2*tmp(ixo^s)*dwc(hxo^s))
221  if (present(ldw)) then
222  ldw(ixo^s)=tmp(ixo^s)* &
223  max(zero,min(tmp2(ixo^s),&
224  (dwc(hxo^s)*tmp(ixo^s)+2*abs(dwc(ixo^s)))*third))
225  end if
226  if (present(rdw)) then
227  rdw(ixo^s)=tmp(ixo^s)* &
228  max(zero,min(tmp2(ixo^s),&
229  (2*dwc(hxo^s)*tmp(ixo^s)+abs(dwc(ixo^s)))*third))
230  end if
231  case (limiter_cada)
232  ! This limiter has been rewritten in the usual form, and uses a division
233  ! of the gradients.
234  if (present(ldw)) then
235  ! Cada Left variant
236  ! Compute theta, but avoid division by zero
237  tmp(ixo^s)=dwc(hxo^s)/(dwc(ixo^s) + sign(eps, dwc(ixo^s)))
238  tmp2(ixo^s)=(2+tmp(ixo^s))*third
239  ldw(ixo^s)= max(zero,min(tmp2(ixo^s), &
240  max(-cadalfa*tmp(ixo^s), &
241  min(cadbeta*tmp(ixo^s), tmp2(ixo^s), &
242  cadgamma)))) * dwc(ixo^s)
243  end if
244 
245  if (present(rdw)) then
246  ! Cada Right variant
247  tmp(ixo^s)=dwc(ixo^s)/(dwc(hxo^s) + sign(eps, dwc(hxo^s)))
248  tmp2(ixo^s)=(2+tmp(ixo^s))*third
249  rdw(ixo^s)= max(zero,min(tmp2(ixo^s), &
250  max(-cadalfa*tmp(ixo^s), &
251  min(cadbeta*tmp(ixo^s), tmp2(ixo^s), &
252  cadgamma)))) * dwc(hxo^s)
253  end if
254  case (limiter_cada3)
255  rdelinv=one/(cada3_radius*dxlevel(idims))**2
256  tmpeta(ixo^s)=(dwc(ixo^s)**2+dwc(hxo^s)**2)*rdelinv
257  if (present(ldw)) then
258  tmp(ixo^s)=dwc(hxo^s)/(dwc(ixo^s) + sign(eps, dwc(ixo^s)))
259  ldwa(ixo^s)=(two+tmp(ixo^s))*third
260  where(tmpeta(ixo^s)<=one-cadepsilon)
261  ldw(ixo^s)=ldwa(ixo^s)
262  elsewhere(tmpeta(ixo^s)>=one+cadepsilon)
263  ldwb(ixo^s)= max(zero,min(ldwa(ixo^s), max(-cadalfa*tmp(ixo^s), &
264  min(cadbeta*tmp(ixo^s), ldwa(ixo^s), cadgamma))))
265  ldw(ixo^s)=ldwb(ixo^s)
266  elsewhere
267  ldwb(ixo^s)= max(zero,min(ldwa(ixo^s), max(-cadalfa*tmp(ixo^s), &
268  min(cadbeta*tmp(ixo^s), ldwa(ixo^s), cadgamma))))
269  tmp2(ixo^s)=(tmpeta(ixo^s)-one)*invcadepsilon
270  ldw(ixo^s)=half*( (one-tmp2(ixo^s))*ldwa(ixo^s) &
271  +(one+tmp2(ixo^s))*ldwb(ixo^s))
272  endwhere
273  ldw(ixo^s)=ldw(ixo^s) * dwc(ixo^s)
274  end if
275 
276  if (present(rdw)) then
277  tmp(ixo^s)=dwc(ixo^s)/(dwc(hxo^s) + sign(eps, dwc(hxo^s)))
278  ldwa(ixo^s)=(two+tmp(ixo^s))*third
279  where(tmpeta(ixo^s)<=one-cadepsilon)
280  rdw(ixo^s)=ldwa(ixo^s)
281  elsewhere(tmpeta(ixo^s)>=one+cadepsilon)
282  ldwb(ixo^s)= max(zero,min(ldwa(ixo^s), max(-cadalfa*tmp(ixo^s), &
283  min(cadbeta*tmp(ixo^s), ldwa(ixo^s), cadgamma))))
284  rdw(ixo^s)=ldwb(ixo^s)
285  elsewhere
286  ldwb(ixo^s)= max(zero,min(ldwa(ixo^s), max(-cadalfa*tmp(ixo^s), &
287  min(cadbeta*tmp(ixo^s), ldwa(ixo^s), cadgamma))))
288  tmp2(ixo^s)=(tmpeta(ixo^s)-one)*invcadepsilon
289  rdw(ixo^s)=half*( (one-tmp2(ixo^s))*ldwa(ixo^s) &
290  +(one+tmp2(ixo^s))*ldwb(ixo^s))
291  endwhere
292  rdw(ixo^s)=rdw(ixo^s) * dwc(hxo^s)
293  end if
294  case(limiter_schmid)
295  tmpeta(ixo^s)=(sqrt(0.4d0*(dwc(ixo^s)**2+dwc(hxo^s)**2)))&
296  /((a2max+cadepsilon)*dxlevel(idims)**2)
297  if(present(ldw)) then
298  tmp(ixo^s)=dwc(hxo^s)/(dwc(ixo^s)+sign(eps,dwc(ixo^s)))
299  ldwa(ixo^s)=(two+tmp(ixo^s))*third
300  where(tmpeta(ixo^s)<=one-cadepsilon)
301  ldw(ixo^s)=ldwa(ixo^s)
302  else where(tmpeta(ixo^s)>=one+cadepsilon)
303  ldwb(ixo^s)=max(zero,min(ldwa(ixo^s),max(-tmp(ixo^s),&
304  min(cadbeta*tmp(ixo^s),ldwa(ixo^s),1.5d0))))
305  ldw(ixo^s)=ldwb(ixo^s)
306  else where
307  ldwb(ixo^s)=max(zero,min(ldwa(ixo^s),max(-tmp(ixo^s),&
308  min(cadbeta*tmp(ixo^s),ldwa(ixo^s),1.5d0))))
309  tmp2(ixo^s)=(tmpeta(ixo^s)-one)*invcadepsilon
310  ldw(ixo^s)=half*((one-tmp2(ixo^s))*ldwa(ixo^s)&
311  +(one+tmp2(ixo^s))*ldwb(ixo^s))
312  end where
313  ldw(ixo^s)=ldw(ixo^s)*dwc(ixo^s)
314  end if
315  if(present(rdw)) then
316  tmp(ixo^s)=dwc(ixo^s)/(dwc(hxo^s)+sign(eps,dwc(hxo^s)))
317  ldwa(ixo^s)=(two+tmp(ixo^s))*third
318  where(tmpeta(ixo^s)<=one-cadepsilon)
319  rdw(ixo^s)=ldwa(ixo^s)
320  else where(tmpeta(ixo^s)>=one+cadepsilon)
321  ldwb(ixo^s)=max(zero,min(ldwa(ixo^s),max(-tmp(ixo^s),&
322  min(cadbeta*tmp(ixo^s),ldwa(ixo^s),1.5d0))))
323  rdw(ixo^s)=ldwb(ixo^s)
324  else where
325  ldwb(ixo^s)=max(zero,min(ldwa(ixo^s),max(-tmp(ixo^s),&
326  min(cadbeta*tmp(ixo^s), ldwa(ixo^s), 1.5d0))))
327  tmp2(ixo^s)=(tmpeta(ixo^s)-one)*invcadepsilon
328  rdw(ixo^s)=half*((one-tmp2(ixo^s))*ldwa(ixo^s)&
329  +(one+tmp2(ixo^s))*ldwb(ixo^s))
330  end where
331  rdw(ixo^s)=rdw(ixo^s)*dwc(hxo^s)
332  end if
333  case default
334  call mpistop("Error in dwLimiter: unknown limiter")
335  end select
336 
337  end subroutine dwlimiter2
338 
339 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:112
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:40
Module containing the MP5 (fifth order) flux scheme.
Definition: mod_mp5.t:2
integer, parameter limiter_teno5ad
Definition: mod_limiter.t:38
integer, parameter limiter_woodward
Definition: mod_limiter.t:16
integer, parameter limiter_weno7
Definition: mod_limiter.t:39
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:46
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:131
integer, parameter limiter_exeno7
Definition: mod_limiter.t:41
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