MPI-AMRVAC  3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
mod_venk.t
Go to the documentation of this file.
1 module mod_venk
2  ! Venkatakrishnan limiter
3  !
4  ! 2019.10.11 coded up by nanami;
5  !
6  ! see Venkatakrishnan 1993 for this limiter;
7 
8 
9  implicit none
10  private
11 
12  public :: venklimiter
13 
14 contains
15 
16  subroutine venklimiter(ixI^L,iL^L,idims,dxdim,w,wLC,wRC)
18 
19  integer, intent(in) :: ixi^l, il^l, idims
20  double precision, intent(in) :: dxdim
21  double precision, intent(in) :: w(ixi^s,1:nw)
22  double precision, intent(inout) :: wrc(ixi^s,1:nw),wlc(ixi^s,1:nw)
23  !> local
24  integer :: im^l, imm^l, imp^l
25  integer :: ilm^l, ilp^l, ilpp^l
26  double precision :: wmax(ixi^s,1:nw),wmin(ixi^s,1:nw)
27  double precision :: westp(ixi^s,1:nw),westm(ixi^s,1:nw)
28  double precision :: phi1(ixi^s,1:nw),phi2(ixi^s,1:nw)
29  double precision :: phi3(ixi^s,1:nw),phi4(ixi^s,1:nw)
30  double precision :: phim(ixi^s,1:nw),phip(ixi^s,1:nw),phi(ixi^s,1:nw)
31  double precision :: deltap(ixi^s,1:nw),deltam(ixi^s,1:nw)
32  double precision :: eps2
33  double precision, parameter :: venk_omega = 1.d-12
34  double precision, parameter :: venk_k = 0.3d0
35 
36  im^l=il^l;
37  immax^d=ilmax^d+kr(idims,^d);
38  ilm^l=il^l-kr(idims,^d);
39  imm^l=im^l-kr(idims,^d);
40  ilp^l=il^l+kr(idims,^d);
41  imp^l=im^l+kr(idims,^d);
42  ilpp^l=ilp^l+kr(idims,^d);
43 
44  eps2 = (venk_k * dxdim) ** 3
45 
46  wmax(im^s,1:nwflux) = max(w(imm^s,1:nwflux), w(im^s,1:nwflux), w(imp^s,1:nwflux))
47  wmin(im^s,1:nwflux) = min(w(imm^s,1:nwflux), w(im^s,1:nwflux), w(imp^s,1:nwflux))
48  !> use central difference approximation as (eq.5) and take phi = 1
49  westp(im^s,1:nwflux) = w(im^s,1:nwflux) + (w(imp^s,1:nwflux) - w(imm^s,1:nwflux)) * 0.25d0
50  westm(im^s,1:nwflux) = w(im^s,1:nwflux) - (w(imp^s,1:nwflux) - w(imm^s,1:nwflux)) * 0.25d0
51 
52  !> (eq.30) & (eq.31)
53  deltap = 0
54  deltam = 0
55  phi1 = 0
56  phi2 = 0
57  phip = 1
58  where(westp(im^s,1:nwflux) .gt. w(im^s,1:nwflux))
59  deltap(im^s,1:nwflux) = wmax(im^s,1:nwflux) - w(im^s,1:nwflux)
60  deltam(im^s,1:nwflux) = westp(im^s,1:nwflux) - w(im^s,1:nwflux)
61  deltam(im^s,1:nwflux) = sign(dabs(deltam(im^s,1:nwflux)) + venk_omega, deltam(im^s,1:nwflux))
62  phi1(im^s,1:nwflux) = (deltap(im^s,1:nwflux)**2 + eps2) + 2.d0 * deltap(im^s,1:nwflux) * deltam(im^s,1:nwflux)
63  phi2(im^s,1:nwflux) = deltap(im^s,1:nwflux)**2 + 2.d0 * deltam(im^s,1:nwflux)**2 + deltap(im^s,1:nwflux) * deltam(im^s,1:nwflux) + eps2
64  phip(im^s,1:nwflux) = phi1(im^s,1:nwflux) / phi2(im^s, 1:nwflux)
65  elsewhere(westp(im^s,1:nwflux) .lt. w(im^s,1:nwflux))
66  deltap(im^s,1:nwflux) = wmin(im^s,1:nwflux) - w(im^s,1:nwflux)
67  deltam(im^s,1:nwflux) = westp(im^s,1:nwflux) - w(im^s,1:nwflux)
68  deltam(im^s,1:nwflux) = sign(dabs(deltam(im^s,1:nwflux)) + venk_omega, deltam(im^s,1:nwflux))
69  phi1(im^s,1:nwflux) = (deltap(im^s,1:nwflux)**2 + eps2) + 2.d0 * deltap(im^s,1:nwflux) * deltam(im^s,1:nwflux)
70  phi2(im^s,1:nwflux) = deltap(im^s,1:nwflux)**2 + 2.d0 * deltam(im^s,1:nwflux)**2 + deltap(im^s,1:nwflux) * deltam(im^s,1:nwflux) + eps2
71  phip(im^s,1:nwflux) = phi1(im^s,1:nwflux) / phi2(im^s, 1:nwflux)
72  elsewhere
73  phip(im^s,1:nwflux) = one
74  endwhere
75 
76  deltap = 0
77  deltam = 0
78  phi3 = 0
79  phi4 = 0
80  phim = 0
81  where(westm(im^s,1:nwflux) .lt. w(im^s,1:nwflux))
82  deltap(im^s,1:nwflux) = - (wmax(im^s,1:nwflux) - w(im^s,1:nwflux))
83  deltam(im^s,1:nwflux) = westm(im^s,1:nwflux) - w(im^s,1:nwflux)
84  deltam(im^s,1:nwflux) = sign(dabs(deltam(im^s,1:nwflux)) + venk_omega, deltam(im^s,1:nwflux))
85  phi3(im^s,1:nwflux) = (deltap(im^s,1:nwflux)**2 + eps2) + 2.d0 * deltap(im^s,1:nwflux) * deltam(im^s,1:nwflux)
86  phi4(im^s,1:nwflux) = deltap(im^s,1:nwflux)**2 + 2.d0 * deltam(im^s,1:nwflux)**2 + deltap(im^s,1:nwflux) * deltam(im^s,1:nwflux) + eps2
87  phim(im^s,1:nwflux) = phi3(im^s,1:nwflux) / phi4(im^s, 1:nwflux)
88  elsewhere(westm(im^s,1:nwflux) .gt. w(im^s,1:nwflux))
89  deltap(im^s,1:nwflux) = - (wmin(im^s,1:nwflux) - w(im^s,1:nwflux))
90  deltam(im^s,1:nwflux) = westm(im^s,1:nwflux) - w(im^s,1:nwflux)
91  deltam(im^s,1:nwflux) = sign(dabs(deltam(im^s,1:nwflux)) + venk_omega, deltam(im^s,1:nwflux))
92  phi3(im^s,1:nwflux) = (deltap(im^s,1:nwflux)**2 + eps2) + 2.d0 * deltap(im^s,1:nwflux) * deltam(im^s,1:nwflux)
93  phi4(im^s,1:nwflux) = deltap(im^s,1:nwflux)**2 + 2.d0 * deltam(im^s,1:nwflux)**2 + deltap(im^s,1:nwflux) * deltam(im^s,1:nwflux) + eps2
94  phim(im^s,1:nwflux) = phi3(im^s,1:nwflux) / phi4(im^s, 1:nwflux)
95  elsewhere
96  phim(im^s,1:nwflux) = one
97  endwhere
98  !> (eq.3)
99  phi(im^s,1:nwflux) = min(phim(im^s,1:nwflux),phip(im^s,1:nwflux))
100  !> (eq.5)
101  wlc(il^s,1:nwflux) = w(il^s,1:nwflux) + 0.25d0 * (w(ilp^s,1:nwflux)-w(ilm^s,1:nwflux)) * phi(il^s,1:nwflux)
102  wrc(il^s,1:nwflux) = w(ilp^s,1:nwflux) - 0.25d0 * (w(ilpp^s,1:nwflux)-w(il^s,1:nwflux)) * phi(ilp^s,1:nwflux)
103 
104  end subroutine venklimiter
105 
106 end module mod_venk
This module contains definitions of global parameters and variables and some generic functions/subrou...
integer, dimension(3, 3) kr
Kronecker delta tensor.
integer, dimension(:), allocatable, parameter d
subroutine, public venklimiter(ixIL, iLL, idims, dxdim, w, wLC, wRC)
Definition: mod_venk.t:17