MPI-AMRVAC  2.2
The MPI - Adaptive Mesh Refinement - Versatile Advection Code
mod_tvd.t
Go to the documentation of this file.
1 !> Subroutines for TVD-MUSCL schemes
2 module mod_tvd
3 
4  implicit none
5  private
6 
7  public :: tvdlimit
8  public :: tvdlimit2
9  public :: entropyfix
10 
11 contains
12 
13  subroutine tvdlimit(method,qdt,ixI^L,ixO^L,idim^LIM,s,qt,snew,fC,dx^D,x)
15 
16  character(len=*), intent(in) :: method
17  double precision, intent(in) :: qdt, qt, dx^D
18  integer, intent(in) :: ixI^L, ixO^L, idim^LIM
19  double precision, dimension(ixI^S,nw) :: w, wnew
20  type(state) :: s, snew
21  double precision, intent(in) :: x(ixi^s,1:ndim)
22  double precision :: fC(ixi^s,1:nwflux,1:ndim)
23 
24  integer :: idims, ixIC^L, jxIC^L
25  double precision, dimension(ixI^S,nw) :: wR, wL
26 
27  associate(w=>s%w,wnew=>snew%w)
28  do idims= idim^lim
29  ixicmax^d=ixomax^d+kr(idims,^d); ixicmin^d=ixomin^d-2*kr(idims,^d);
30  wl(ixic^s,1:nw)=w(ixic^s,1:nw)
31  jxic^l=ixic^l+kr(idims,^d);
32  wr(ixic^s,1:nw)=w(jxic^s,1:nw)
33  call tvdlimit2(method,qdt,ixi^l,ixic^l,ixo^l,idims,wl,wr,wnew,x,fc,dx^d)
34  end do
35  end associate
36  end subroutine tvdlimit
37 
38  subroutine tvdlimit2(method,qdt,ixI^L,ixIC^L,ixO^L,idims,wL,wR,wnew,x,fC,dx^D)
39 
40  ! Limit the flow variables in wnew according to typetvd.
41  ! wroeC is based on wL and wR.
42  ! If method=='tvd' an extra adtdx**2*jumpC is added to phiC for 2nd order
43  ! accuracy in time.
44 
46  use mod_physics_roe
47 
48  character(len=*), intent(in) :: method
49  double precision, intent(in) :: qdt, dx^D
50  integer, intent(in) :: ixI^L, ixIC^L, ixO^L, idims
51  double precision, dimension(ixG^T,nw) :: wL, wR
52  double precision, intent(in) :: x(ixi^s,1:ndim)
53  double precision :: wnew(ixi^s,1:nw)
54  double precision :: fC(ixi^s,1:nwflux,1:ndim)
55 
56  double precision:: workroe(ixg^t,1:nworkroe)
57  double precision, dimension(ixG^T,nw) :: wroeC
58  double precision, dimension(ixG^T) :: phiC, rphiC, jumpC, adtdxC, smallaC
59  double precision :: dxinv(1:ndim)
60  integer :: hxO^L, ixC^L, jxC^L, jxIC^L, iw, il
61  !-----------------------------------------------------------------------------
62 
63  hxo^l=ixo^l-kr(idims,^d);
64  ixcmax^d=ixomax^d; ixcmin^d=hxomin^d;
65 
66  jxc^l=ixc^l+kr(idims,^d);
67  jxic^l=ixic^l+kr(idims,^d);
68 
69  call phys_average(wl,wr,x,ixic^l,idims,wroec,workroe)
70 
71  ^d&dxinv(^d)=qdt/dx^d;
72 
73  ! A loop on characteristic variables to calculate the dissipative flux phiC.
74  do il=1,nwflux
75  !Calculate the jump in the il-th characteristic variable: L(wroe)*dw
76  call phys_get_eigenjump(wl,wr,wroec,x,ixic^l,il,idims,smallac,adtdxc,jumpc,workroe)
77 
78  ! Normalize the eigenvalue "a" (and its limit "smalla" if needed):
79  if (slab_uniform) then
80  adtdxc(ixic^s)=adtdxc(ixic^s)*dxinv(idims)
81  if (typeentropy(il)=='harten' .or. typeentropy(il)=='powell')&
82  smallac(ixic^s)=smallac(ixic^s)*dxinv(idims)
83  else
84  adtdxc(ixic^s)=adtdxc(ixic^s)*qdt*block%surfaceC(ixic^s,idims)*&
85  2.0d0/(block%dvolume(ixic^s)+block%dvolume(jxic^s))
86  if (typeentropy(il)=='harten' .or. typeentropy(il)=='powell')&
87  smallac(ixic^s)=smallac(ixic^s)*qdt*block%surfaceC(ixic^s,idims)*&
88  2.0d0/(block%dvolume(ixic^s)+block%dvolume(jxic^s))
89  endif
90 
91  ! Calculate the flux limiter function phi
92  call getphi(method,jumpc,adtdxc,smallac,ixi^l,ixic^l,ixc^l,il,idims,phic)
93 
94  !Add R(iw,il)*phiC(il) to each variable iw in wnew
95  do iw=1,nwflux
96  call phys_rtimes(phic,wroec,ixc^l,iw,il,idims,rphic,workroe)
97 
98  if (slab_uniform) then
99  rphic(ixc^s)=rphic(ixc^s)*half
100  fc(ixc^s,iw,idims)=fc(ixc^s,iw,idims)+rphic(ixc^s)
101  wnew(ixo^s,iw)=wnew(ixo^s,iw)+rphic(ixo^s)-rphic(hxo^s)
102  else
103  rphic(ixc^s)=rphic(ixc^s)*quarter* &
104  (block%dvolume(ixc^s)+block%dvolume(jxc^s))
105  fc(ixc^s,iw,idims)=fc(ixc^s,iw,idims)+rphic(ixc^s)
106  wnew(ixo^s,iw)=wnew(ixo^s,iw)+(rphic(ixo^s)-rphic(hxo^s)) &
107  /block%dvolume(ixo^s)
108  endif
109  end do !iw
110  end do !il
111 
112  end subroutine tvdlimit2
113 
114  subroutine getphi(method,jumpC,adtdxC,smallaC,ixI^L,ixIC^L,ixC^L,il,idims,phiC)
116  ! Calculate the dissipative flux from jumpC=L*dw and adtdx=eigenvalue*dt/dx.
117  ! Add Lax-Wendroff type correction if method=='tvd'.
118  ! Limit according to method and typetvd.
119  use mod_limiter
121 
122  character(len=*), intent(in) :: method
123  integer, intent(in) :: ixI^L, ixIC^L, ixC^L, il, idims
124  double precision, dimension(ixG^T) :: jumpC, adtdxC, smallaC, phiC
125 
126  double precision, dimension(ixG^T) :: ljumpC, tmp
127  integer :: jxC^L, ix^L, hx^L
128  !-----------------------------------------------------------------------------
129 
130  if(method=='tvdmu'.or.method=='tvdmu1')then
131  ! In the MUSCL scheme phi=|a|*jump, apply entropy fix to it
132  if(typeentropy(il)=='nul'.or.typeentropy(il)=='ratio')then
133  phic(ixc^s)=abs(adtdxc(ixc^s))*jumpc(ixc^s)
134  else
135  where(abs(adtdxc(ixc^s))>=smallac(ixc^s))
136  phic(ixc^s)=abs(adtdxc(ixc^s))*jumpc(ixc^s)
137  elsewhere
138  phic(ixc^s)=half*(smallac(ixc^s)+adtdxc(ixc^s)**2/smallac(ixc^s))&
139  *jumpc(ixc^s)
140  endwhere
141  endif
142  ! That's all for the MUSCL scheme
143  return
144  endif
145 
146  if(method=='tvd')then
147  !Entropy fix to |a|-a**2
148  if(typeentropy(il)=='nul'.or.typeentropy(il)=='ratio')then
149  phic(ixic^s)=abs(adtdxc(ixic^s))-adtdxc(ixic^s)**2
150  else
151  where(abs(adtdxc(ixic^s))>=smallac(ixic^s))
152  phic(ixic^s)=abs(adtdxc(ixic^s))-adtdxc(ixic^s)**2
153  elsewhere
154  phic(ixic^s)=half*smallac(ixic^s)+&
155  (half/smallac(ixic^s)-one)*adtdxc(ixic^s)**2
156  endwhere
157  endif
158  else
159  !Entropy fix to |a|
160  if(typeentropy(il)=='nul'.or.typeentropy(il)=='ratio')then
161  phic(ixic^s)=abs(adtdxc(ixic^s))
162  else
163  where(abs(adtdxc(ixic^s))>=smallac(ixic^s))
164  phic(ixic^s)=abs(adtdxc(ixic^s))
165  elsewhere
166  phic(ixic^s)=half*smallac(ixic^s)+&
167  half/smallac(ixic^s)*adtdxc(ixic^s)**2
168  endwhere
169  endif
170  endif
171 
172  jxc^l=ixc^l+kr(idims,^d);
173  hxmin^d=ixicmin^d; hxmax^d=ixicmax^d-kr(idims,^d);
174  ix^l=hx^l+kr(idims,^d);
175 
176  if (.not. limiter_symmetric(typelimiter)) then
177  call mpistop("TVD only supports symmetric limiters")
178  end if
179 
180  select case(typetvd)
181  case('roe')
182  call dwlimiter2(jumpc,ixi^l,ixic^l,idims,typelimiter,ldw=ljumpc)
183  where(adtdxc(ixc^s)<=0)
184  phic(ixc^s)=phic(ixc^s)*(jumpc(ixc^s)-ljumpc(jxc^s))
185  elsewhere
186  phic(ixc^s)=phic(ixc^s)*(jumpc(ixc^s)-ljumpc(ixc^s))
187  end where
188  !extra (a*lambda)**2*delta
189  if(method=='tvd')phic(ixc^s)=phic(ixc^s)+adtdxc(ixc^s)**2*jumpc(ixc^s)
190  case('sweby')
191  !Sweby eqs.4.11-4.15, but no 0.5 ?!
192  phic(ixic^s)=phic(ixic^s)*jumpc(ixic^s)
193  call dwlimiter2(phic,ixi^l,ixic^l,idims,typelimiter,ldw=ljumpc)
194  where(adtdxc(ixc^s)<=0)
195  phic(ixc^s)=phic(ixc^s)-ljumpc(jxc^s)
196  elsewhere
197  phic(ixc^s)=phic(ixc^s)-ljumpc(ixc^s)
198  end where
199  !extra (a*lambda)**2*delta
200  if(method=='tvd')phic(ixc^s)=phic(ixc^s)+adtdxc(ixc^s)**2*jumpc(ixc^s)
201  case('yee')
202  !eq.3.51 with correction
203  call dwlimiter2(jumpc,ixi^l,ixic^l,idims,typelimiter,ldw=ljumpc)
204 
205  !Use phiC as 0.5*(|nu|-nu**2) eq.3.45e for tvd otherwise 0.5*|nu|
206  phic(ixc^s)=half*phic(ixc^s)
207  !gamma*lambda eq.3.51d, use tmp to store agdtdxC
208  where(abs(jumpc(ixc^s))>smalldouble)
209  tmp(ixc^s)=adtdxc(ixc^s)+phic(ixc^s)*&
210  (ljumpc(jxc^s)-ljumpc(ixc^s))/jumpc(ixc^s)
211  elsewhere
212  tmp(ixc^s)=adtdxc(ixc^s)
213  end where
214 
215  !eq.3.51a
216  if(typeentropy(il)=='nul'.or.typeentropy(il)=='ratio')then
217  phic(ixc^s)=-phic(ixc^s)*(ljumpc(jxc^s)+ljumpc(ixc^s))+&
218  abs(tmp(ixc^s))*jumpc(ixc^s)
219  else
220  where(abs(tmp(ixc^s))>=smallac(ixc^s))
221  phic(ixc^s)=-phic(ixc^s)*(ljumpc(jxc^s)+ljumpc(ixc^s))+&
222  abs(tmp(ixc^s))*jumpc(ixc^s)
223  elsewhere
224  phic(ixc^s)=-phic(ixc^s)*(ljumpc(jxc^s)+ljumpc(ixc^s))+&
225  (half*smallac(ixc^s)+half/smallac(ixc^s)*&
226  tmp(ixc^s)**2)*jumpc(ixc^s)
227  endwhere
228  endif
229  case('harten')
230  !See Ryu, section 2.3
231  !Use phiC as 0.5*(|nu|-nu**2)*jumpC eq.3.45b,e
232  phic(ixic^s)=half*phic(ixic^s)*jumpc(ixic^s)
233  call dwlimiter2(phic,ixi^l,ixic^l,idims,typelimiter,ldw=ljumpc)
234 
235  !gamma*lambda eq.3.45d, use tmp as agdtdxC
236  where(abs(jumpc(ixc^s))>smalldouble)
237  tmp(ixc^s)=adtdxc(ixc^s)+&
238  (ljumpc(jxc^s)-ljumpc(ixc^s))/jumpc(ixc^s)
239  elsewhere
240  tmp(ixc^s)=adtdxc(ixc^s)
241  end where
242  !eq.3.45a with correction
243  if(typeentropy(il)=='nul'.or.typeentropy(il)=='ratio')then
244  phic(ixc^s)=-ljumpc(jxc^s)-ljumpc(ixc^s)+jumpc(ixc^s)*&
245  abs(tmp(ixc^s))
246  else
247  where(abs(tmp(ixc^s))>=smallac(ixc^s))
248  phic(ixc^s)=-ljumpc(jxc^s)-ljumpc(ixc^s)+jumpc(ixc^s)*&
249  abs(tmp(ixc^s))
250  elsewhere
251  phic(ixc^s)=-ljumpc(jxc^s)-ljumpc(ixc^s)+jumpc(ixc^s)*&
252  (half*smallac(ixc^s)+half/smallac(ixc^s)*tmp(ixc^s)**2)
253  endwhere
254  endif
255  !extra -(a*lambda)**2*delta
256  case default
257  call mpistop("Error in TVDLimit: Unknown TVD type")
258  end select
259 
260  end subroutine getphi
261 
262  subroutine entropyfix(ix^L,il,aL,aR,a,smalla)
264  ! Apply entropyfix based on typeentropy(il),aL,aR, and a
265  ! Calculate "smalla" (Harten,Powell) or modify "a" (ratio)
266 
268 
269  integer, intent(in) :: ix^L, il
270  double precision, dimension(ixG^T) :: aL, aR, a, smalla
271  !-----------------------------------------------------------------------------
272 
273  select case(typeentropy(il))
274  case('harten')
275  smalla(ix^s)=max(zero,a(ix^s)-al(ix^s),ar(ix^s)-a(ix^s))
276  case('powell')
277  smalla(ix^s)=max(zero,two*(ar(ix^s)-al(ix^s)))
278  !!case('ratio')
279  !! where(aL(ix^S)<zero .and. aR(ix^S)>zero)&
280  !! a(ix^S)=a(ix^S)-2*aR(ix^S)*aL(ix^S)/(aR(ix^S)-aL(ix^S))
281  case('yee')
282  ! This has been done in geteigenjump already
283  case('nul')
284  ! No entropyfix is applied
285  case default
286  call mpistop("No such type of entropy fix")
287  end select
288 
289  end subroutine entropyfix
290 
291 end module mod_tvd
This module contains definitions of global parameters and variables and some generic functions/subrou...
subroutine, public entropyfix(ixL, il, aL, aR, a, smalla)
Definition: mod_tvd.t:263
pure logical function limiter_symmetric(typelim)
Definition: mod_limiter.t:109
procedure(sub_rtimes), pointer phys_rtimes
character(len=std_len) typetvd
Which type of TVD method to use.
procedure(sub_average), pointer phys_average
logical slab_uniform
uniform Cartesian geometry or not (stretched Cartesian)
Module with slope/flux limiters.
Definition: mod_limiter.t:2
character(len=std_len), dimension(:), allocatable typeentropy
Which type of entropy fix to use with Riemann-type solvers.
integer, dimension(3, 3) kr
Kronecker delta tensor.
integer, dimension(:), allocatable, parameter d
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:128
subroutine mpistop(message)
Exit MPI-AMRVAC with an error message.
Definition: comm_lib.t:198
procedure(sub_get_eigenjump), pointer phys_get_eigenjump
integer, parameter ndim
Number of spatial dimensions for grid variables.
subroutine getphi(method, jumpC, adtdxC, smallaC, ixIL, ixICL, ixCL, il, idims, phiC)
Definition: mod_tvd.t:115
subroutine, public tvdlimit2(method, qdt, ixIL, ixICL, ixOL, idims, wL, wR, wnew, x, fC, dxD)
Definition: mod_tvd.t:39
Subroutines for TVD-MUSCL schemes.
Definition: mod_tvd.t:2
subroutine, public tvdlimit(method, qdt, ixIL, ixOL, idimLIM, s, qt, snew, fC, dxD, x)
Definition: mod_tvd.t:14