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