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  integer, parameter :: limiter_minmod = 1
15  integer, parameter :: limiter_woodward = 2
16  integer, parameter :: limiter_mcbeta = 3
17  integer, parameter :: limiter_superbee = 4
18  integer, parameter :: limiter_vanleer = 5
19  integer, parameter :: limiter_albada = 6
20  integer, parameter :: limiter_koren = 7
21  integer, parameter :: limiter_cada = 8
22  integer, parameter :: limiter_cada3 = 9
23  integer, parameter :: limiter_venk = 10
24  ! Special cases
25  integer, parameter :: limiter_ppm = 11
26  integer, parameter :: limiter_mp5 = 12
27  integer, parameter :: limiter_weno3 = 13
28  integer, parameter :: limiter_wenoyc3 = 14
29  integer, parameter :: limiter_weno5 = 15
30  integer, parameter :: limiter_weno5nm = 16
31  integer, parameter :: limiter_wenoz5 = 17
32  integer, parameter :: limiter_wenoz5nm = 18
33  integer, parameter :: limiter_wenozp5 = 19
34  integer, parameter :: limiter_wenozp5nm = 20
35  integer, parameter :: limiter_weno7 = 21
36  integer, parameter :: limiter_mpweno7 = 22
37  integer, parameter :: limiter_exeno7 = 23
38 
39 contains
40 
41  integer function limiter_type(namelim)
42  character(len=*), intent(in) :: namelim
43 
44  select case (namelim)
45  case ('minmod')
47  case ('woodward')
49  case ('mcbeta')
51  case ('superbee')
53  case ('vanleer')
55  case ('albada')
57  case ('koren')
59  case ('cada')
61  case ('cada3')
63  case('venk')
65  case ('ppm')
67  case ('mp5')
69  case ('weno3')
71  case ('wenoyc3')
73  case ('weno5')
75  case ('weno5nm')
77  case ('wenoz5')
79  case ('wenoz5nm')
81  case ('wenozp5')
83  case ('wenozp5nm')
85  case ('weno7')
87  case ('mpweno7')
89  case ('exeno7')
91 
92  case default
93  limiter_type = -1
94  write(*,*) 'Unknown limiter: ', namelim
95  call mpistop("No such limiter")
96  end select
97  end function limiter_type
98 
99  pure logical function limiter_symmetric(typelim)
100  integer, intent(in) :: typelim
101 
102  select case (typelim)
104  limiter_symmetric = .false.
105  case default
106  limiter_symmetric = .true.
107  end select
108  end function limiter_symmetric
109 
110  !> Limit the centered dwC differences within ixC for iw in direction idim.
111  !> The limiter is chosen according to typelimiter.
112  !>
113  !> Note that this subroutine is called from upwindLR (hence from methods
114  !> like tvdlf, hancock, hll(c) etc) or directly from tvd.t,
115  !> but also from the gradientS and divvectorS subroutines in geometry.t
116  !> Accordingly, the typelimiter here corresponds to one of limiter
117  !> or one of gradient_limiter.
118  subroutine dwlimiter2(dwC,ixI^L,ixC^L,idims,typelim,ldw,rdw)
121 
122  integer, intent(in) :: ixI^L, ixC^L, idims
123  double precision, intent(in) :: dwC(ixi^s)
124  integer, intent(in) :: typelim
125  !> Result using left-limiter (same as right for symmetric)
126  double precision, intent(out), optional :: ldw(ixi^s)
127  !> Result using right-limiter (same as left for symmetric)
128  double precision, intent(out), optional :: rdw(ixi^s)
129 
130  double precision :: tmp(ixi^s), tmp2(ixi^s)
131  integer :: ixO^L, hxO^L
132  double precision, parameter :: qsmall=1.d-12, qsmall2=2.d-12
133  double precision, parameter :: eps = sqrt(epsilon(1.0d0))
134 
135  ! mcbeta limiter parameter value
136  double precision, parameter :: c_mcbeta=1.4d0
137  ! cada limiter parameter values
138  double precision, parameter :: cadalfa=0.5d0, cadbeta=2.0d0, cadgamma=1.6d0
139  ! full third order cada limiter
140  double precision :: rdelinv
141  double precision :: ldwA(ixi^s),ldwB(ixi^s),tmpeta(ixi^s)
142  double precision, parameter :: cadepsilon=1.d-14, invcadepsilon=1.d14,cada3_radius=0.1d0
143  !-----------------------------------------------------------------------------
144 
145  ! Contract indices in idim for output.
146  ixomin^d=ixcmin^d+kr(idims,^d); ixomax^d=ixcmax^d;
147  hxo^l=ixo^l-kr(idims,^d);
148 
149  ! About the notation: the conventional argument theta (the ratio of slopes)
150  ! would be given by dwC(ixO^S)/dwC(hxO^S). However, in the end one
151  ! multiplies phi(theta) by dwC(hxO^S), which is incorporated in the
152  ! equations below. The minmod limiter can for example be written as:
153  ! A:
154  ! max(0.0d0, min(1.0d0, dwC(ixO^S)/dwC(hxO^S))) * dwC(hxO^S)
155  ! B:
156  ! tmp(ixO^S)*max(0.0d0,min(abs(dwC(ixO^S)),tmp(ixO^S)*dwC(hxO^S)))
157  ! where tmp(ixO^S)=sign(1.0d0,dwC(ixO^S))
158 
159  select case (typelim)
160  case (limiter_minmod)
161  ! Minmod limiter eq(3.51e) and (eq.3.38e) with omega=1
162  tmp(ixo^s)=sign(one,dwc(ixo^s))
163  tmp(ixo^s)=tmp(ixo^s)* &
164  max(zero,min(abs(dwc(ixo^s)),tmp(ixo^s)*dwc(hxo^s)))
165  if (present(ldw)) ldw(ixo^s) = tmp(ixo^s)
166  if (present(rdw)) rdw(ixo^s) = tmp(ixo^s)
167  case (limiter_woodward)
168  ! Woodward and Collela limiter (eq.3.51h), a factor of 2 is pulled out
169  tmp(ixo^s)=sign(one,dwc(ixo^s))
170  tmp(ixo^s)=2*tmp(ixo^s)* &
171  max(zero,min(abs(dwc(ixo^s)),tmp(ixo^s)*dwc(hxo^s),&
172  tmp(ixo^s)*quarter*(dwc(hxo^s)+dwc(ixo^s))))
173  if (present(ldw)) ldw(ixo^s) = tmp(ixo^s)
174  if (present(rdw)) rdw(ixo^s) = tmp(ixo^s)
175  case (limiter_mcbeta)
176  ! Woodward and Collela limiter, with factor beta
177  tmp(ixo^s)=sign(one,dwc(ixo^s))
178  tmp(ixo^s)=tmp(ixo^s)* &
179  max(zero,min(c_mcbeta*abs(dwc(ixo^s)),c_mcbeta*tmp(ixo^s)*dwc(hxo^s),&
180  tmp(ixo^s)*half*(dwc(hxo^s)+dwc(ixo^s))))
181  if (present(ldw)) ldw(ixo^s) = tmp(ixo^s)
182  if (present(rdw)) rdw(ixo^s) = tmp(ixo^s)
183  case (limiter_superbee)
184  ! Roes superbee limiter (eq.3.51i)
185  tmp(ixo^s)=sign(one,dwc(ixo^s))
186  tmp(ixo^s)=tmp(ixo^s)* &
187  max(zero,min(2*abs(dwc(ixo^s)),tmp(ixo^s)*dwc(hxo^s)),&
188  min(abs(dwc(ixo^s)),2*tmp(ixo^s)*dwc(hxo^s)))
189  if (present(ldw)) ldw(ixo^s) = tmp(ixo^s)
190  if (present(rdw)) rdw(ixo^s) = tmp(ixo^s)
191  case (limiter_vanleer)
192  ! van Leer limiter (eq 3.51f), but a missing delta2=1.D-12 is added
193  tmp(ixo^s)=2*max(dwc(hxo^s)*dwc(ixo^s),zero) &
194  /(dwc(ixo^s)+dwc(hxo^s)+qsmall)
195  if (present(ldw)) ldw(ixo^s) = tmp(ixo^s)
196  if (present(rdw)) rdw(ixo^s) = tmp(ixo^s)
197  case (limiter_albada)
198  ! Albada limiter (eq.3.51g) with delta2=1D.-12
199  tmp(ixo^s)=(dwc(hxo^s)*(dwc(ixo^s)**2+qsmall)&
200  +dwc(ixo^s)*(dwc(hxo^s)**2+qsmall))&
201  /(dwc(ixo^s)**2+dwc(hxo^s)**2+qsmall2)
202  if (present(ldw)) ldw(ixo^s) = tmp(ixo^s)
203  if (present(rdw)) rdw(ixo^s) = tmp(ixo^s)
204  case (limiter_koren)
205  tmp(ixo^s)=sign(one,dwc(ixo^s))
206  tmp2(ixo^s)=min(2*abs(dwc(ixo^s)),2*tmp(ixo^s)*dwc(hxo^s))
207  if (present(ldw)) then
208  ldw(ixo^s)=tmp(ixo^s)* &
209  max(zero,min(tmp2(ixo^s),&
210  (dwc(hxo^s)*tmp(ixo^s)+2*abs(dwc(ixo^s)))*third))
211  end if
212  if (present(rdw)) then
213  rdw(ixo^s)=tmp(ixo^s)* &
214  max(zero,min(tmp2(ixo^s),&
215  (2*dwc(hxo^s)*tmp(ixo^s)+abs(dwc(ixo^s)))*third))
216  end if
217  case (limiter_cada)
218  ! This limiter has been rewritten in the usual form, and uses a division
219  ! of the gradients.
220  if (present(ldw)) then
221  ! Cada Left variant
222  ! Compute theta, but avoid division by zero
223  tmp(ixo^s)=dwc(hxo^s)/(dwc(ixo^s) + sign(eps, dwc(ixo^s)))
224  tmp2(ixo^s)=(2+tmp(ixo^s))*third
225  ldw(ixo^s)= max(zero,min(tmp2(ixo^s), &
226  max(-cadalfa*tmp(ixo^s), &
227  min(cadbeta*tmp(ixo^s), tmp2(ixo^s), &
228  cadgamma)))) * dwc(ixo^s)
229  end if
230 
231  if (present(rdw)) then
232  ! Cada Right variant
233  tmp(ixo^s)=dwc(ixo^s)/(dwc(hxo^s) + sign(eps, dwc(hxo^s)))
234  tmp2(ixo^s)=(2+tmp(ixo^s))*third
235  rdw(ixo^s)= max(zero,min(tmp2(ixo^s), &
236  max(-cadalfa*tmp(ixo^s), &
237  min(cadbeta*tmp(ixo^s), tmp2(ixo^s), &
238  cadgamma)))) * dwc(hxo^s)
239  end if
240  case (limiter_cada3)
241  rdelinv=one/(cada3_radius*dxlevel(idims))**2
242  tmpeta(ixo^s)=(dwc(ixo^s)**2+dwc(hxo^s)**2)*rdelinv
243 
244  if (present(ldw)) then
245  tmp(ixo^s)=dwc(hxo^s)/(dwc(ixo^s) + sign(eps, dwc(ixo^s)))
246  ldwa(ixo^s)=(two+tmp(ixo^s))*third
247  ldwb(ixo^s)= max(zero,min(ldwa(ixo^s), max(-cadalfa*tmp(ixo^s), &
248  min(cadbeta*tmp(ixo^s), ldwa(ixo^s), cadgamma))))
249  where(tmpeta(ixo^s)<=one-cadepsilon)
250  ldw(ixo^s)=ldwa(ixo^s)
251  elsewhere(tmpeta(ixo^s)>=one+cadepsilon)
252  ldw(ixo^s)=ldwb(ixo^s)
253  elsewhere
254  tmp2(ixo^s)=(tmpeta(ixo^s)-one)*invcadepsilon
255  ldw(ixo^s)=half*( (one-tmp2(ixo^s))*ldwa(ixo^s) &
256  +(one+tmp2(ixo^s))*ldwb(ixo^s))
257  endwhere
258  ldw(ixo^s)=ldw(ixo^s) * dwc(ixo^s)
259  end if
260 
261  if (present(rdw)) then
262  tmp(ixo^s)=dwc(ixo^s)/(dwc(hxo^s) + sign(eps, dwc(hxo^s)))
263  ldwa(ixo^s)=(two+tmp(ixo^s))*third
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  where(tmpeta(ixo^s)<=one-cadepsilon)
267  rdw(ixo^s)=ldwa(ixo^s)
268  elsewhere(tmpeta(ixo^s)>=one+cadepsilon)
269  rdw(ixo^s)=ldwb(ixo^s)
270  elsewhere
271  tmp2(ixo^s)=(tmpeta(ixo^s)-one)*invcadepsilon
272  rdw(ixo^s)=half*( (one-tmp2(ixo^s))*ldwa(ixo^s) &
273  +(one+tmp2(ixo^s))*ldwb(ixo^s))
274  endwhere
275  rdw(ixo^s)=rdw(ixo^s) * dwc(hxo^s)
276  end if
277 
278  case default
279  call mpistop("Error in dwLimiter: unknown limiter")
280  end select
281 
282  end subroutine dwlimiter2
283 
284 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:100
integer, parameter limiter_superbee
Definition: mod_limiter.t:17
integer, parameter limiter_cada3
Definition: mod_limiter.t:22
integer, parameter limiter_wenoz5
Definition: mod_limiter.t:31
integer, parameter limiter_wenozp5
Definition: mod_limiter.t:33
integer, parameter limiter_mp5
Definition: mod_limiter.t:26
integer, parameter limiter_albada
Definition: mod_limiter.t:19
integer, parameter limiter_wenoyc3
Definition: mod_limiter.t:28
integer, parameter limiter_mpweno7
Definition: mod_limiter.t:36
Module containing the MP5 (fifth order) flux scheme.
Definition: mod_mp5.t:2
integer, parameter limiter_woodward
Definition: mod_limiter.t:15
integer, parameter limiter_weno7
Definition: mod_limiter.t:35
integer, parameter limiter_ppm
Definition: mod_limiter.t:25
Module with slope/flux limiters.
Definition: mod_limiter.t:2
integer, parameter limiter_wenozp5nm
Definition: mod_limiter.t:34
integer, parameter limiter_weno5nm
Definition: mod_limiter.t:30
integer function limiter_type(namelim)
Definition: mod_limiter.t:42
integer, parameter limiter_koren
Definition: mod_limiter.t:20
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_mcbeta
Definition: mod_limiter.t:16
integer, parameter limiter_venk
Definition: mod_limiter.t:23
Definition: mod_ppm.t:1
integer, parameter limiter_exeno7
Definition: mod_limiter.t:37
subroutine mpistop(message)
Exit MPI-AMRVAC with an error message.
Definition: comm_lib.t:198
integer, parameter limiter_weno5
Definition: mod_limiter.t:29
double precision, dimension(ndim) dxlevel
integer, parameter limiter_vanleer
Definition: mod_limiter.t:18
integer, parameter limiter_minmod
Definition: mod_limiter.t:14
integer, parameter limiter_cada
Definition: mod_limiter.t:21
subroutine dwlimiter2(dwC, ixIL, ixCL, idims, typelim, ldw, rdw)
Limit the centered dwC differences within ixC for iw in direction idim. The limiter is chosen accordi...
Definition: mod_limiter.t:119
integer, parameter limiter_weno3
Definition: mod_limiter.t:27
integer, parameter limiter_wenoz5nm
Definition: mod_limiter.t:32