13 subroutine tvdlimit(method,qdt,ixI^L,ixO^L,idim^LIM,s,qt,snew,fC,dxs,x)
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)
24 integer :: idims, ixic^
l, jxic^
l
25 double precision,
dimension(ixI^S,nw) :: wr, wl
27 associate(w=>s%w,wnew=>snew%w)
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)
38 subroutine tvdlimit2(method,qdt,ixI^L,ixIC^L,ixO^L,idims,wL,wR,wnew,x,fC,dxs)
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)
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
63 hxo^
l=ixo^
l-
kr(idims,^
d);
64 ixcmax^
d=ixomax^
d; ixcmin^
d=hxomin^
d;
66 jxc^
l=ixc^
l+
kr(idims,^
d);
67 jxic^
l=ixic^
l+
kr(idims,^
d);
76 call phys_get_eigenjump(wl,wr,wroec,x,ixic^
l,il,idims,smallac,adtdxc,jumpc,workroe)
80 adtdxc(ixic^s)=adtdxc(ixic^s)*dxinv(idims)
82 smallac(ixic^s)=smallac(ixic^s)*dxinv(idims)
84 adtdxc(ixic^s)=adtdxc(ixic^s)*qdt*
block%surfaceC(ixic^s,idims)*&
85 2.0d0/(
block%dvolume(ixic^s)+
block%dvolume(jxic^s))
87 smallac(ixic^s)=smallac(ixic^s)*qdt*
block%surfaceC(ixic^s,idims)*&
88 2.0d0/(
block%dvolume(ixic^s)+
block%dvolume(jxic^s))
92 call getphi(method,jumpc,adtdxc,smallac,ixi^
l,ixic^
l,ixc^
l,il,idims,phic)
96 call phys_rtimes(phic,wroec,ixc^
l,iw,il,idims,rphic,workroe)
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)
103 rphic(ixc^s)=rphic(ixc^s)*quarter* &
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)
114 subroutine getphi(method,jumpC,adtdxC,smallaC,ixI^L,ixIC^L,ixC^L,il,idims,phiC)
122 integer,
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
126 double precision,
dimension(ixG^T) :: ljumpC, tmp
127 integer :: jxC^L, ix^L, hx^L, typelimiter
134 phic(ixc^s)=abs(adtdxc(ixc^s))*jumpc(ixc^s)
136 where(abs(adtdxc(ixc^s))>=smallac(ixc^s))
137 phic(ixc^s)=abs(adtdxc(ixc^s))*jumpc(ixc^s)
139 phic(ixc^s)=half*(smallac(ixc^s)+adtdxc(ixc^s)**2/smallac(ixc^s))&
150 phic(ixic^s)=abs(adtdxc(ixic^s))-adtdxc(ixic^s)**2
152 where(abs(adtdxc(ixic^s))>=smallac(ixic^s))
153 phic(ixic^s)=abs(adtdxc(ixic^s))-adtdxc(ixic^s)**2
155 phic(ixic^s)=half*smallac(ixic^s)+&
156 (half/smallac(ixic^s)-one)*adtdxc(ixic^s)**2
162 phic(ixic^s)=abs(adtdxc(ixic^s))
164 where(abs(adtdxc(ixic^s))>=smallac(ixic^s))
165 phic(ixic^s)=abs(adtdxc(ixic^s))
167 phic(ixic^s)=half*smallac(ixic^s)+&
168 half/smallac(ixic^s)*adtdxc(ixic^s)**2
173 jxc^l=ixc^l+
kr(idims,^
d);
174 hxmin^
d=ixicmin^
d; hxmax^
d=ixicmax^
d-
kr(idims,^
d);
175 ix^l=hx^l+
kr(idims,^
d);
178 call mpistop(
"TVD only supports symmetric limiters")
183 call dwlimiter2(jumpc,ixi^l,ixic^l,idims,typelimiter,ldw=ljumpc)
184 where(adtdxc(ixc^s)<=0)
185 phic(ixc^s)=phic(ixc^s)*(jumpc(ixc^s)-ljumpc(jxc^s))
187 phic(ixc^s)=phic(ixc^s)*(jumpc(ixc^s)-ljumpc(ixc^s))
190 if(method==
fs_tvd)phic(ixc^s)=phic(ixc^s)+adtdxc(ixc^s)**2*jumpc(ixc^s)
193 phic(ixic^s)=phic(ixic^s)*jumpc(ixic^s)
194 call dwlimiter2(phic,ixi^l,ixic^l,idims,typelimiter,ldw=ljumpc)
195 where(adtdxc(ixc^s)<=0)
196 phic(ixc^s)=phic(ixc^s)-ljumpc(jxc^s)
198 phic(ixc^s)=phic(ixc^s)-ljumpc(ixc^s)
201 if(method==
fs_tvd)phic(ixc^s)=phic(ixc^s)+adtdxc(ixc^s)**2*jumpc(ixc^s)
204 call dwlimiter2(jumpc,ixi^l,ixic^l,idims,typelimiter,ldw=ljumpc)
207 phic(ixc^s)=half*phic(ixc^s)
209 where(abs(jumpc(ixc^s))>smalldouble)
210 tmp(ixc^s)=adtdxc(ixc^s)+phic(ixc^s)*&
211 (ljumpc(jxc^s)-ljumpc(ixc^s))/jumpc(ixc^s)
213 tmp(ixc^s)=adtdxc(ixc^s)
218 phic(ixc^s)=-phic(ixc^s)*(ljumpc(jxc^s)+ljumpc(ixc^s))+&
219 abs(tmp(ixc^s))*jumpc(ixc^s)
221 where(abs(tmp(ixc^s))>=smallac(ixc^s))
222 phic(ixc^s)=-phic(ixc^s)*(ljumpc(jxc^s)+ljumpc(ixc^s))+&
223 abs(tmp(ixc^s))*jumpc(ixc^s)
225 phic(ixc^s)=-phic(ixc^s)*(ljumpc(jxc^s)+ljumpc(ixc^s))+&
226 (half*smallac(ixc^s)+half/smallac(ixc^s)*&
227 tmp(ixc^s)**2)*jumpc(ixc^s)
233 phic(ixic^s)=half*phic(ixic^s)*jumpc(ixic^s)
234 call dwlimiter2(phic,ixi^l,ixic^l,idims,typelimiter,ldw=ljumpc)
237 where(abs(jumpc(ixc^s))>smalldouble)
238 tmp(ixc^s)=adtdxc(ixc^s)+&
239 (ljumpc(jxc^s)-ljumpc(ixc^s))/jumpc(ixc^s)
241 tmp(ixc^s)=adtdxc(ixc^s)
245 phic(ixc^s)=-ljumpc(jxc^s)-ljumpc(ixc^s)+jumpc(ixc^s)*&
248 where(abs(tmp(ixc^s))>=smallac(ixc^s))
249 phic(ixc^s)=-ljumpc(jxc^s)-ljumpc(ixc^s)+jumpc(ixc^s)*&
252 phic(ixc^s)=-ljumpc(jxc^s)-ljumpc(ixc^s)+jumpc(ixc^s)*&
253 (half*smallac(ixc^s)+half/smallac(ixc^s)*tmp(ixc^s)**2)
258 call mpistop(
"Error in TVDLimit: Unknown TVD type")
270 integer,
intent(in) :: ix^
l, il
271 double precision,
dimension(ixG^T) :: al, ar, a, smalla
276 smalla(ix^s)=max(zero,a(ix^s)-al(ix^s),ar(ix^s)-a(ix^s))
278 smalla(ix^s)=max(zero,two*(ar(ix^s)-al(ix^s)))
287 call mpistop(
"No such type of entropy fix")
subroutine mpistop(message)
Exit MPI-AMRVAC with an error message.
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, parameter fs_tvd
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.
pure logical function limiter_symmetric(typelim)
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...
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.
subroutine, public tvdlimit2(method, qdt, ixIL, ixICL, ixOL, idims, wL, wR, wnew, x, fC, dxs)
subroutine getphi(method, jumpC, adtdxC, smallaC, ixIL, ixICL, ixCL, il, idims, phiC)
subroutine, public entropyfix(ixL, il, aL, aR, a, smalla)
subroutine, public tvdlimit(method, qdt, ixIL, ixOL, idimLIM, s, qt, snew, fC, dxs, x)