MPI-AMRVAC 3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
Loading...
Searching...
No Matches
mod_limiter.t
Go to the documentation of this file.
1!> Module with slope/flux limiters
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
42contains
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 gradientL 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 double precision, parameter :: qsmall=1.d-12, qsmall2=2.d-12
144 double precision, parameter :: eps = sqrt(epsilon(1.0d0))
145
146 ! mcbeta limiter parameter value
147 double precision, parameter :: c_mcbeta=1.4d0
148 ! cada limiter parameter values
149 double precision, parameter :: cadalfa=0.5d0, cadbeta=2.0d0, cadgamma=1.6d0
150 ! full third order cada limiter
151 double precision :: rdelinv
152 double precision :: ldwA(ixI^S),ldwB(ixI^S),tmpeta(ixI^S)
153 double precision, parameter :: cadepsilon=1.d-14, invcadepsilon=1.d14,cada3_radius=0.1d0
154 integer :: ixO^L, hxO^L, ix^D, hx^D
155 !-----------------------------------------------------------------------------
156
157 ! Contract indices in idim for output.
158 ixomin^d=ixcmin^d+kr(idims,^d); ixomax^d=ixcmax^d;
159 hxo^l=ixo^l-kr(idims,^d);
160 hx^d=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 {!DEC$ VECTOR ALWAYS
256 do ix^db=ixomin^db,ixomax^db\}
257 tmpeta(ix^d)=(dwc(ix^d)**2+dwc(ix^d-hx^d)**2)*rdelinv
258 if (present(ldw)) then
259 tmp(ix^d)=dwc(ix^d-hx^d)/(dwc(ix^d)+sign(eps,dwc(ix^d)))
260 ldwa(ix^d)=(two+tmp(ix^d))*third
261 if(tmpeta(ix^d)<=one-cadepsilon) then
262 ldw(ix^d)=ldwa(ix^d)
263 else if(tmpeta(ix^d)>=one+cadepsilon) then
264 ldw(ix^d)=max(zero,min(ldwa(ix^d),max(-cadalfa*tmp(ix^d),&
265 min(cadbeta*tmp(ix^d),ldwa(ix^d),cadgamma))))
266 else
267 tmp2(ix^d)=(tmpeta(ix^d)-one)*invcadepsilon
268 ldw(ix^d)=half*((one-tmp2(ix^d))*ldwa(ix^d)+(one+tmp2(ix^d))*&
269 max(zero,min(ldwa(ix^d),max(-cadalfa*tmp(ix^d),&
270 min(cadbeta*tmp(ix^d),ldwa(ix^d),cadgamma)))))
271 end if
272 ldw(ix^d)=ldw(ix^d)*dwc(ix^d)
273 end if
274 if (present(rdw)) then
275 tmp(ix^d)=dwc(ix^d)/(dwc(ix^d-hx^d)+sign(eps,dwc(ix^d-hx^d)))
276 ldwa(ix^d)=(two+tmp(ix^d))*third
277 if(tmpeta(ix^d)<=one-cadepsilon) then
278 rdw(ix^d)=ldwa(ix^d)
279 else if(tmpeta(ix^d)>=one+cadepsilon) then
280 rdw(ix^d)=max(zero,min(ldwa(ix^d),max(-cadalfa*tmp(ix^d),&
281 min(cadbeta*tmp(ix^d),ldwa(ix^d),cadgamma))))
282 else
283 tmp2(ix^d)=(tmpeta(ix^d)-one)*invcadepsilon
284 rdw(ix^d)=half*((one-tmp2(ix^d))*ldwa(ix^d)+(one+tmp2(ix^d))*&
285 max(zero,min(ldwa(ix^d),max(-cadalfa*tmp(ix^d),&
286 min(cadbeta*tmp(ix^d),ldwa(ix^d),cadgamma)))))
287 end if
288 rdw(ix^d)=rdw(ix^d)*dwc(ix^d-hx^d)
289 end if
290 {end do\}
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
336end module mod_limiter
subroutine, public mpistop(message)
Exit MPI-AMRVAC with an error message.
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(^nd) dxlevel
store unstretched cell size of current level
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)
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
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
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...
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