MPI-AMRVAC  3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
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,dxs,x)
15 
16  integer, intent(in) :: method
17  double precision, intent(in) :: qdt, qt, dxs(ndim)
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,dxs)
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,dxs)
39 
40  ! Limit the flow variables in wnew according to typetvd.
41  ! wroeC is based on wL and wR.
42  ! If method=fs_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  integer, intent(in) :: method
49  double precision, intent(in) :: qdt, dxs(ndim)
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  dxinv=qdt/dxs
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)
115 
116  ! Calculate the dissipative flux from jumpC=L*dw and adtdx=eigenvalue*dt/dx.
117  ! Add Lax-Wendroff type correction if method=fs_tvd.
118  ! Limit according to method and typetvd.
119  use mod_limiter
121  use mod_comm_lib, only: mpistop
122 
123  integer, intent(in) :: method
124  integer, intent(in) :: ixI^L, ixIC^L, ixC^L, il, idims
125  double precision, dimension(ixG^T) :: jumpC, adtdxC, smallaC, phiC
126 
127  double precision, dimension(ixG^T) :: ljumpC, tmp
128  integer :: jxC^L, ix^L, hx^L, typelimiter
129  !-----------------------------------------------------------------------------
130 
131  typelimiter=type_limiter(block%level)
132  if(method==fs_tvdmu)then
133  ! In the MUSCL scheme phi=|a|*jump, apply entropy fix to it
134  if(typeentropy(il)=='nul'.or.typeentropy(il)=='ratio')then
135  phic(ixc^s)=abs(adtdxc(ixc^s))*jumpc(ixc^s)
136  else
137  where(abs(adtdxc(ixc^s))>=smallac(ixc^s))
138  phic(ixc^s)=abs(adtdxc(ixc^s))*jumpc(ixc^s)
139  elsewhere
140  phic(ixc^s)=half*(smallac(ixc^s)+adtdxc(ixc^s)**2/smallac(ixc^s))&
141  *jumpc(ixc^s)
142  endwhere
143  endif
144  ! That's all for the MUSCL scheme
145  return
146  endif
147 
148  if(method==fs_tvd)then
149  !Entropy fix to |a|-a**2
150  if(typeentropy(il)=='nul'.or.typeentropy(il)=='ratio')then
151  phic(ixic^s)=abs(adtdxc(ixic^s))-adtdxc(ixic^s)**2
152  else
153  where(abs(adtdxc(ixic^s))>=smallac(ixic^s))
154  phic(ixic^s)=abs(adtdxc(ixic^s))-adtdxc(ixic^s)**2
155  elsewhere
156  phic(ixic^s)=half*smallac(ixic^s)+&
157  (half/smallac(ixic^s)-one)*adtdxc(ixic^s)**2
158  endwhere
159  endif
160  else
161  !Entropy fix to |a|
162  if(typeentropy(il)=='nul'.or.typeentropy(il)=='ratio')then
163  phic(ixic^s)=abs(adtdxc(ixic^s))
164  else
165  where(abs(adtdxc(ixic^s))>=smallac(ixic^s))
166  phic(ixic^s)=abs(adtdxc(ixic^s))
167  elsewhere
168  phic(ixic^s)=half*smallac(ixic^s)+&
169  half/smallac(ixic^s)*adtdxc(ixic^s)**2
170  endwhere
171  endif
172  endif
173 
174  jxc^l=ixc^l+kr(idims,^d);
175  hxmin^d=ixicmin^d; hxmax^d=ixicmax^d-kr(idims,^d);
176  ix^l=hx^l+kr(idims,^d);
177 
178  if (.not. limiter_symmetric(typelimiter)) then
179  call mpistop("TVD only supports symmetric limiters")
180  end if
181 
182  select case(typetvd)
183  case('roe')
184  call dwlimiter2(jumpc,ixi^l,ixic^l,idims,typelimiter,ldw=ljumpc)
185  where(adtdxc(ixc^s)<=0)
186  phic(ixc^s)=phic(ixc^s)*(jumpc(ixc^s)-ljumpc(jxc^s))
187  elsewhere
188  phic(ixc^s)=phic(ixc^s)*(jumpc(ixc^s)-ljumpc(ixc^s))
189  end where
190  !extra (a*lambda)**2*delta
191  if(method==fs_tvd)phic(ixc^s)=phic(ixc^s)+adtdxc(ixc^s)**2*jumpc(ixc^s)
192  case('sweby')
193  !Sweby eqs.4.11-4.15, but no 0.5 ?!
194  phic(ixic^s)=phic(ixic^s)*jumpc(ixic^s)
195  call dwlimiter2(phic,ixi^l,ixic^l,idims,typelimiter,ldw=ljumpc)
196  where(adtdxc(ixc^s)<=0)
197  phic(ixc^s)=phic(ixc^s)-ljumpc(jxc^s)
198  elsewhere
199  phic(ixc^s)=phic(ixc^s)-ljumpc(ixc^s)
200  end where
201  !extra (a*lambda)**2*delta
202  if(method==fs_tvd)phic(ixc^s)=phic(ixc^s)+adtdxc(ixc^s)**2*jumpc(ixc^s)
203  case('yee')
204  !eq.3.51 with correction
205  call dwlimiter2(jumpc,ixi^l,ixic^l,idims,typelimiter,ldw=ljumpc)
206 
207  !Use phiC as 0.5*(|nu|-nu**2) eq.3.45e for tvd otherwise 0.5*|nu|
208  phic(ixc^s)=half*phic(ixc^s)
209  !gamma*lambda eq.3.51d, use tmp to store agdtdxC
210  where(abs(jumpc(ixc^s))>smalldouble)
211  tmp(ixc^s)=adtdxc(ixc^s)+phic(ixc^s)*&
212  (ljumpc(jxc^s)-ljumpc(ixc^s))/jumpc(ixc^s)
213  elsewhere
214  tmp(ixc^s)=adtdxc(ixc^s)
215  end where
216 
217  !eq.3.51a
218  if(typeentropy(il)=='nul'.or.typeentropy(il)=='ratio')then
219  phic(ixc^s)=-phic(ixc^s)*(ljumpc(jxc^s)+ljumpc(ixc^s))+&
220  abs(tmp(ixc^s))*jumpc(ixc^s)
221  else
222  where(abs(tmp(ixc^s))>=smallac(ixc^s))
223  phic(ixc^s)=-phic(ixc^s)*(ljumpc(jxc^s)+ljumpc(ixc^s))+&
224  abs(tmp(ixc^s))*jumpc(ixc^s)
225  elsewhere
226  phic(ixc^s)=-phic(ixc^s)*(ljumpc(jxc^s)+ljumpc(ixc^s))+&
227  (half*smallac(ixc^s)+half/smallac(ixc^s)*&
228  tmp(ixc^s)**2)*jumpc(ixc^s)
229  endwhere
230  endif
231  case('harten')
232  !See Ryu, section 2.3
233  !Use phiC as 0.5*(|nu|-nu**2)*jumpC eq.3.45b,e
234  phic(ixic^s)=half*phic(ixic^s)*jumpc(ixic^s)
235  call dwlimiter2(phic,ixi^l,ixic^l,idims,typelimiter,ldw=ljumpc)
236 
237  !gamma*lambda eq.3.45d, use tmp as agdtdxC
238  where(abs(jumpc(ixc^s))>smalldouble)
239  tmp(ixc^s)=adtdxc(ixc^s)+&
240  (ljumpc(jxc^s)-ljumpc(ixc^s))/jumpc(ixc^s)
241  elsewhere
242  tmp(ixc^s)=adtdxc(ixc^s)
243  end where
244  !eq.3.45a with correction
245  if(typeentropy(il)=='nul'.or.typeentropy(il)=='ratio')then
246  phic(ixc^s)=-ljumpc(jxc^s)-ljumpc(ixc^s)+jumpc(ixc^s)*&
247  abs(tmp(ixc^s))
248  else
249  where(abs(tmp(ixc^s))>=smallac(ixc^s))
250  phic(ixc^s)=-ljumpc(jxc^s)-ljumpc(ixc^s)+jumpc(ixc^s)*&
251  abs(tmp(ixc^s))
252  elsewhere
253  phic(ixc^s)=-ljumpc(jxc^s)-ljumpc(ixc^s)+jumpc(ixc^s)*&
254  (half*smallac(ixc^s)+half/smallac(ixc^s)*tmp(ixc^s)**2)
255  endwhere
256  endif
257  !extra -(a*lambda)**2*delta
258  case default
259  call mpistop("Error in TVDLimit: Unknown TVD type")
260  end select
261 
262  end subroutine getphi
263 
264  subroutine entropyfix(ix^L,il,aL,aR,a,smalla)
265 
266  ! Apply entropyfix based on typeentropy(il),aL,aR, and a
267  ! Calculate "smalla" (Harten,Powell) or modify "a" (ratio)
268 
270  use mod_comm_lib, only: mpistop
271 
272  integer, intent(in) :: ix^l, il
273  double precision, dimension(ixG^T) :: al, ar, a, smalla
274  !-----------------------------------------------------------------------------
275 
276  select case(typeentropy(il))
277  case('harten')
278  smalla(ix^s)=max(zero,a(ix^s)-al(ix^s),ar(ix^s)-a(ix^s))
279  case('powell')
280  smalla(ix^s)=max(zero,two*(ar(ix^s)-al(ix^s)))
281  !!case('ratio')
282  !! where(aL(ix^S)<zero .and. aR(ix^S)>zero)&
283  !! a(ix^S)=a(ix^S)-2*aR(ix^S)*aL(ix^S)/(aR(ix^S)-aL(ix^S))
284  case('yee')
285  ! This has been done in geteigenjump already
286  case('nul')
287  ! No entropyfix is applied
288  case default
289  call mpistop("No such type of entropy fix")
290  end select
291 
292  end subroutine entropyfix
293 
294 end module mod_tvd
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...
character(len=std_len), dimension(:), allocatable typeentropy
Which type of entropy fix to use with Riemann-type solvers.
type(state), pointer block
Block pointer for using one block and its previous state.
integer, dimension(3, 3) kr
Kronecker delta tensor.
integer, parameter ndim
Number of spatial dimensions for grid variables.
integer, parameter fs_tvdmu
integer, dimension(:), allocatable, parameter d
integer, dimension(:), allocatable type_limiter
Type of slope limiter used for reconstructing variables on cell edges.
character(len=std_len) typetvd
Which type of TVD method to use.
logical slab_uniform
uniform Cartesian geometry or not (stretched Cartesian)
Module with slope/flux limiters.
Definition: mod_limiter.t:2
pure logical function limiter_symmetric(typelim)
Definition: mod_limiter.t:110
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
double precision d
Definition: mod_limiter.t:14
procedure(sub_rtimes), pointer phys_rtimes
procedure(sub_get_eigenjump), pointer phys_get_eigenjump
procedure(sub_average), pointer phys_average
Subroutines for TVD-MUSCL schemes.
Definition: mod_tvd.t:2
subroutine, public tvdlimit2(method, qdt, ixIL, ixICL, ixOL, idims, wL, wR, wnew, x, fC, dxs)
Definition: mod_tvd.t:39
subroutine getphi(method, jumpC, adtdxC, smallaC, ixIL, ixICL, ixCL, il, idims, phiC)
Definition: mod_tvd.t:115
subroutine, public entropyfix(ixL, il, aL, aR, a, smalla)
Definition: mod_tvd.t:265
subroutine, public tvdlimit(method, qdt, ixIL, ixOL, idimLIM, s, qt, snew, fC, dxs, x)
Definition: mod_tvd.t:14