12 integer,
parameter :: fastRW_ = 3,fastlw_=4,slowrw_=5,slowlw_=6
13 integer,
parameter :: entroW_ = 8,diverw_=7,alfvrw_=1,alfvlw_=2
32 case(fastrw_,fastlw_,slowrw_,slowlw_)
49 subroutine twofl_average(wL,wR,x,ix^L,idim,wroe,workroe)
52 integer,
intent(in) :: ix^
l, idim
53 double precision,
intent(in) :: wl(ixg^t, nw), wr(ixg^t, nw)
54 double precision,
intent(inout) :: wroe(ixg^t, nw)
55 double precision,
intent(inout) :: workroe(ixg^t,
nworkroe)
56 double precision,
intent(in) :: x(ixg^t, 1:^nd)
58 call average2(wl,wr,x,ix^
l,idim,wroe,workroe(ixg^t,1),workroe(ixg^t,2), &
59 workroe(ixg^t,3),workroe(ixg^t,4),workroe(ixg^t,5),workroe(ixg^t,6), &
60 workroe(ixg^t,7),workroe(ixg^t,8))
62 end subroutine twofl_average
70 subroutine average2(wL,wR,x,ix^L,idim,wroe,cfast,cslow,afast,aslow,csound2,dp, &
73 integer :: ix^
l,idim,idir,jdir,iw
74 double precision,
dimension(ixG^T,nw) :: wl,wr,wroe
75 double precision,
intent(in) :: x(ixg^t,1:
ndim)
76 double precision,
dimension(ixG^T) :: cfast,cslow,afast,aslow,csound2,dp, &
79 if (
ndir==1)
call mpistop(
"twofl with d=11 is the same as HD")
85 wroe(ix^s,
mag(idir))=half*(wl(ix^s,
mag(idir))+wr(ix^s,
mag(idir)))
92 wroe(ix^s,
e_c_)=half*(afast(ix^s)+aslow(ix^s))
94 dp(ix^s)=aslow(ix^s)-afast(ix^s)
96 dp(ix^s)=aslow(ix^s)-afast(ix^s)
100 rhodv(ix^s)=wr(ix^s,
mom_c(idim))-wl(ix^s,
mom_c(idim))-&
113 cfast(ix^s)=sum(wroe(ix^s,
mag(:))**2,dim=
ndim+1)/wroe(ix^s,
rho_c_)+csound2(ix^s)
116 cslow(ix^s)=half*(cfast(ix^s)-dsqrt(cfast(ix^s)**2-&
117 4d0*csound2(ix^s)*wroe(ix^s,
mag(idim))**2/wroe(ix^s,
rho_c_)))
120 cfast(ix^s)=cfast(ix^s)-cslow(ix^s)
123 afast(ix^s)=(csound2(ix^s)-cslow(ix^s))/(cfast(ix^s)-cslow(ix^s))
124 afast(ix^s)=min(one,max(afast(ix^s),zero))
127 aslow(ix^s)=dsqrt(one-afast(ix^s))
130 afast(ix^s)=dsqrt(afast(ix^s))
133 cfast(ix^s)=dsqrt(cfast(ix^s))
136 cslow(ix^s)=dsqrt(cslow(ix^s))
143 where(dabs(wroe(ix^s,
mag(idim)))<smalldouble)&
144 wroe(ix^s,
mag(idim))=smalldouble
148 where(wroe(ix^s,
mag(idir))>=zero)
149 wroe(ix^s,
mag(idir))=one
151 wroe(ix^s,
mag(idir))=-one
156 tmp(ix^s)=dsqrt(wroe(ix^s,
mag(idir))**2+wroe(ix^s,
mag(jdir))**2)
157 where(tmp(ix^s)>smalldouble)
158 wroe(ix^s,
mag(idir))=wroe(ix^s,
mag(idir))/tmp(ix^s)
159 wroe(ix^s,
mag(jdir))=wroe(ix^s,
mag(jdir))/tmp(ix^s)
161 wroe(ix^s,
mag(idir))=dsqrt(half)
162 wroe(ix^s,
mag(jdir))=dsqrt(half)
166 end subroutine average2
178 subroutine twofl_get_eigenjump(wL,wR,wroe,x,ix^L,il,idim,smalla,a,jump,workroe)
181 integer,
intent(in) :: ix^
l,il,idim
182 double precision,
dimension(ixG^T,nw) :: wl,wr,wroe
183 double precision,
intent(in) :: x(ixg^t,1:
ndim)
184 double precision,
dimension(ixG^T) :: smalla,a,jump
185 double precision,
dimension(ixG^T,nworkroe) :: workroe
187 call geteigenjump2(wl,wr,wroe,x,ix^
l,il,idim,smalla,a,jump, &
188 workroe(ixg^t,1),workroe(ixg^t,2), &
189 workroe(ixg^t,3),workroe(ixg^t,4),workroe(ixg^t,5),workroe(ixg^t,6), &
190 workroe(ixg^t,7),workroe(ixg^t,8),workroe(ixg^t,9),workroe(ixg^t,10), &
191 workroe(ixg^t,11),workroe(ixg^t,12),workroe(ixg^t,13))
193 end subroutine twofl_get_eigenjump
205 subroutine geteigenjump2(wL,wR,wroe,x,ix^L,il,idim,smalla,a,jump, &
206 cfast,cslow,afast,aslow,csound2,dp,rhodv,bdv,bdb,cs2L,cs2R,cs2ca2L,cs2ca2R)
210 integer :: ix^
l,il,idim,idir,jdir
211 double precision,
dimension(ixG^T,nw) :: wl,wr,wroe
212 double precision,
intent(in) :: x(ixg^t,1:
ndim)
213 double precision,
dimension(ixG^T) :: smalla,a,jump
214 double precision,
dimension(ixG^T) :: cfast,cslow,afast,aslow,csound2,dp,rhodv
215 double precision,
dimension(ixG^T) :: bdv,bdb
216 double precision,
dimension(ixG^T) :: al,ar,cs2l,cs2r,cs2ca2l,cs2ca2r
224 bdv(ix^s)=wroe(ix^s,
mag(idir))* &
226 if(
ndir==3)bdv(ix^s)=bdv(ix^s)+wroe(ix^s,
mag(jdir))* &
228 bdv(ix^s)=bdv(ix^s)*sign(wroe(ix^s,
rho_c_)**2,wroe(ix^s,
mag(idim)))
230 bdb(ix^s)=wroe(ix^s,
mag(idir))*(wr(ix^s,
mag(idir))-wl(ix^s,
mag(idir)))
231 if(
ndir==3)bdb(ix^s)=bdb(ix^s)+&
232 wroe(ix^s,
mag(jdir))*(wr(ix^s,
mag(jdir))-wl(ix^s,
mag(jdir)))
233 bdb(ix^s)=bdb(ix^s)*dsqrt(csound2(ix^s))*wroe(ix^s,
rho_c_)
239 bdv(ix^s)=wroe(ix^s,
mag(jdir))* &
241 -wroe(ix^s,
mag(idir))* &
243 bdb(ix^s)=wroe(ix^s,
mag(jdir))*(wr(ix^s,
mag(idir))-wl(ix^s,
mag(idir))) &
244 -wroe(ix^s,
mag(idir))*(wr(ix^s,
mag(jdir))-wl(ix^s,
mag(jdir)))
245 bdv(ix^s)=bdv(ix^s)*half*wroe(ix^s,
rho_c_)**2
246 bdb(ix^s)=bdb(ix^s)*half*sign(wroe(ix^s,
rho_c_),wroe(ix^s,
mag(idim)))
251 a(ix^s)=wroe(ix^s,
mom_c(idim))+cfast(ix^s)
252 jump(ix^s)=half/csound2(ix^s)*(&
253 afast(ix^s)*(+cfast(ix^s)*rhodv(ix^s)+dp(ix^s))&
254 +aslow(ix^s)*(-cslow(ix^s)*bdv(ix^s)+bdb(ix^s)))
256 a(ix^s)=wroe(ix^s,
mom_c(idim))-cfast(ix^s)
257 jump(ix^s)=half/csound2(ix^s)*(&
258 afast(ix^s)*(-cfast(ix^s)*rhodv(ix^s)+dp(ix^s))&
259 +aslow(ix^s)*(+cslow(ix^s)*bdv(ix^s)+bdb(ix^s)))
261 a(ix^s)=wroe(ix^s,
mom_c(idim))+cslow(ix^s)
262 jump(ix^s)=half/csound2(ix^s)*(&
263 aslow(ix^s)*(+cslow(ix^s)*rhodv(ix^s)+dp(ix^s))&
264 +afast(ix^s)*(+cfast(ix^s)*bdv(ix^s)-bdb(ix^s)))
266 a(ix^s)=wroe(ix^s,
mom_c(idim))-cslow(ix^s)
267 jump(ix^s)=half/csound2(ix^s)*(&
268 aslow(ix^s)*(-cslow(ix^s)*rhodv(ix^s)+dp(ix^s))&
269 +afast(ix^s)*(-cfast(ix^s)*bdv(ix^s)-bdb(ix^s)))
271 a(ix^s)=wroe(ix^s,
mom_c(idim))
272 jump(ix^s)=wr(ix^s,
rho_c_)-wl(ix^s,
rho_c_)-dp(ix^s)/csound2(ix^s)
275 a(ix^s)=wroe(ix^s,
mom_c(idim))
276 jump(ix^s)=wr(ix^s,
mag(idim))-wl(ix^s,
mag(idim))
282 a(ix^s)=wroe(ix^s,
mom_c(idim))+dabs(wroe(ix^s,
mag(idim)))/wroe(ix^s,
rho_c_)
283 jump(ix^s)=+bdv(ix^s)-bdb(ix^s)
285 a(ix^s)=wroe(ix^s,
mom_c(idim))-dabs(wroe(ix^s,
mag(idim)))/wroe(ix^s,
rho_c_)
286 jump(ix^s)=-bdv(ix^s)-bdb(ix^s)
295 case(
'harten',
'powell',
'ratio')
307 cs2ca2l(ix^s)=cs2l(ix^s)+sum(wl(ix^s,
mag(:))**2,dim=
ndim+1)/wl(ix^s,
rho_c_)
308 cs2ca2r(ix^s)=cs2r(ix^s)+sum(wr(ix^s,
mag(:))**2,dim=
ndim+1)/wr(ix^s,
rho_c_)
311 dsqrt(cs2ca2l(ix^s)**2-4d0*cs2l(ix^s)*wl(ix^s,
mag(idim))**2/wl(ix^s,
rho_c_))
313 dsqrt(cs2ca2r(ix^s)**2-4d0*cs2r(ix^s)*wr(ix^s,
mag(idim))**2/wr(ix^s,
rho_c_))
316 al(ix^s)=al(ix^s) + dsqrt(half*(cs2ca2l(ix^s) + cs2l(ix^s)))
317 ar(ix^s)=ar(ix^s) + dsqrt(half*(cs2ca2r(ix^s) + cs2r(ix^s)))
319 al(ix^s)=al(ix^s) - dsqrt(half*(cs2ca2l(ix^s) + cs2l(ix^s)))
320 ar(ix^s)=ar(ix^s) - dsqrt(half*(cs2ca2r(ix^s) + cs2r(ix^s)))
322 al(ix^s)=al(ix^s) + dsqrt(half*(cs2ca2l(ix^s) - cs2l(ix^s)))
323 ar(ix^s)=ar(ix^s) + dsqrt(half*(cs2ca2r(ix^s) - cs2r(ix^s)))
325 al(ix^s)=al(ix^s) - dsqrt(half*(cs2ca2l(ix^s) - cs2l(ix^s)))
326 ar(ix^s)=ar(ix^s) - dsqrt(half*(cs2ca2r(ix^s) - cs2r(ix^s)))
327 case(entrow_,diverw_)
331 cs2ca2l(ix^s)=dabs(wl(ix^s,
mag(idim)))/dsqrt(wl(ix^s,
rho_c_))
332 cs2ca2r(ix^s)=dabs(wr(ix^s,
mag(idim)))/dsqrt(wr(ix^s,
rho_c_))
334 al(ix^s)=al(ix^s) + cs2ca2l(ix^s)
335 ar(ix^s)=ar(ix^s) + cs2ca2r(ix^s)
337 al(ix^s)=al(ix^s) - cs2ca2l(ix^s)
338 ar(ix^s)=ar(ix^s) - cs2ca2r(ix^s)
344 end subroutine geteigenjump2
347 subroutine twofl_rtimes(q,w,ix^L,iw,il,idim,rq,workroe)
350 integer,
intent(in) :: ix^
l, iw, il, idim
351 double precision,
intent(in) :: w(ixg^t, nw), q(ixg^t)
352 double precision,
intent(inout) :: rq(ixg^t)
353 double precision,
intent(inout) :: workroe(ixg^t,
nworkroe)
355 call rtimes2(q,w,ix^
l,iw,il,idim,rq,&
356 workroe(ixg^t,1),workroe(ixg^t,2), &
357 workroe(ixg^t,3),workroe(ixg^t,4),workroe(ixg^t,5),workroe(ixg^t,6), &
358 workroe(ixg^t,7),workroe(ixg^t,14),workroe(ixg^t,15))
360 end subroutine twofl_rtimes
363 subroutine rtimes2(q,wroe,ix^L,iw,il,idim,rq, &
364 cfast,cslow,afast,aslow,csound2,dp,rhodv,bv,v2a2)
367 integer :: ix^
l,iw,il,idim,idir,jdir
368 double precision :: wroe(ixg^t,nw)
369 double precision,
dimension(ixG^T) :: q,rq
370 double precision,
dimension(ixG^T) :: cfast,cslow,afast,aslow,csound2,dp,rhodv
371 double precision,
dimension(ixG^T) :: bv,v2a2
378 case(fastrw_,fastlw_)
379 rq(ix^s)=q(ix^s)*afast(ix^s)
380 case(slowrw_,slowlw_)
381 rq(ix^s)=q(ix^s)*aslow(ix^s)
384 case(diverw_,alfvrw_,alfvlw_)
387 else if(iw ==
e_c_)
then
390 v2a2(ix^s)=half*sum(wroe(ix^s,
mom_c(:))**2,dim=
ndim+1)+ &
393 bv(ix^s)=wroe(ix^s,
mag(idir))*wroe(ix^s,
mom_c(idir))
394 if(
ndir==3)bv(ix^s)=bv(ix^s)+wroe(ix^s,
mag(jdir))*wroe(ix^s,
mom_c(jdir))
395 bv(ix^s)=bv(ix^s)*sign(one,wroe(ix^s,
mag(idim)))
396 else if(il==alfvrw_)
then
398 bv(ix^s)=(wroe(ix^s,
mag(jdir))*wroe(ix^s,
mom_c(idir))-&
399 wroe(ix^s,
mag(idir))*wroe(ix^s,
mom_c(jdir)))
404 rq(ix^s)=q(ix^s)*(-aslow(ix^s)*cslow(ix^s)*bv(ix^s)+afast(ix^s)*&
405 (v2a2(ix^s)+cfast(ix^s)*(cfast(ix^s)+wroe(ix^s,
mom_c(idim)))))
407 rq(ix^s)=q(ix^s)*(+aslow(ix^s)*cslow(ix^s)*bv(ix^s)+afast(ix^s)*&
408 (v2a2(ix^s)+cfast(ix^s)*(cfast(ix^s)-wroe(ix^s,
mom_c(idim)))))
410 rq(ix^s)=q(ix^s)*(+afast(ix^s)*cfast(ix^s)*bv(ix^s)+aslow(ix^s)*&
411 (v2a2(ix^s)+cslow(ix^s)*(cslow(ix^s)+wroe(ix^s,
mom_c(idim)))))
413 rq(ix^s)=q(ix^s)*(-afast(ix^s)*cfast(ix^s)*bv(ix^s)+aslow(ix^s)*&
414 (v2a2(ix^s)+cslow(ix^s)*(cslow(ix^s)-wroe(ix^s,
mom_c(idim)))))
416 rq(ix^s)= q(ix^s)*half*sum(wroe(ix^s,
mom_c(:))**2,dim=
ndim+1)
419 rq(ix^s)= q(ix^s)*wroe(ix^s,
mag(idim))
424 rq(ix^s)=+q(ix^s)*bv(ix^s)
426 rq(ix^s)=-q(ix^s)*bv(ix^s)
428 else if(any(
mom_c(:)==iw))
then
429 if(iw==
mom_c(idim))
then
432 rq(ix^s)=q(ix^s)*afast(ix^s)*(wroe(ix^s,iw)+cfast(ix^s))
434 rq(ix^s)=q(ix^s)*afast(ix^s)*(wroe(ix^s,iw)-cfast(ix^s))
436 rq(ix^s)=q(ix^s)*aslow(ix^s)*(wroe(ix^s,iw)+cslow(ix^s))
438 rq(ix^s)=q(ix^s)*aslow(ix^s)*(wroe(ix^s,iw)-cslow(ix^s))
440 rq(ix^s)=q(ix^s)*wroe(ix^s,iw)
441 case(diverw_,alfvlw_,alfvrw_)
447 rq(ix^s)=q(ix^s)*(afast(ix^s)*wroe(ix^s,iw)-aslow(ix^s)*&
448 cslow(ix^s)*wroe(ix^s,
mag(1)-
mom_c(1)+iw)*sign(one,wroe(ix^s,
mag(idim))))
450 rq(ix^s)=q(ix^s)*(afast(ix^s)*wroe(ix^s,iw)+aslow(ix^s)*&
451 cslow(ix^s)*wroe(ix^s,
mag(1)-
mom_c(1)+iw)*sign(one,wroe(ix^s,
mag(idim))))
453 rq(ix^s)=q(ix^s)*(aslow(ix^s)*wroe(ix^s,iw)+afast(ix^s)*&
454 cfast(ix^s)*wroe(ix^s,
mag(1)-
mom_c(1)+iw)*sign(one,wroe(ix^s,
mag(idim))))
456 rq(ix^s)=q(ix^s)*(aslow(ix^s)*wroe(ix^s,iw)-afast(ix^s)*&
457 cfast(ix^s)*wroe(ix^s,
mag(1)-
mom_c(1)+iw)*sign(one,wroe(ix^s,
mag(idim))))
459 rq(ix^s)=q(ix^s)*wroe(ix^s,iw)
463 if(iw==
mom_c(idir))
then
464 rq(ix^s)=+q(ix^s)*wroe(ix^s,
mag(jdir))
466 rq(ix^s)=-q(ix^s)*wroe(ix^s,
mag(idir))
469 if(iw==
mom_c(idir))
then
470 rq(ix^s)=-q(ix^s)*wroe(ix^s,
mag(jdir))
472 rq(ix^s)=+q(ix^s)*wroe(ix^s,
mag(idir))
476 else if(any(
mag(:)==iw))
then
477 if(iw==
mag(idim))
then
485 case(fastrw_,fastlw_)
486 rq(ix^s)=+q(ix^s)*aslow(ix^s)*dsqrt(csound2(ix^s))*wroe(ix^s,iw)&
488 case(slowrw_,slowlw_)
489 rq(ix^s)=-q(ix^s)*afast(ix^s)*dsqrt(csound2(ix^s))*wroe(ix^s,iw)&
491 case(entrow_,diverw_)
493 case(alfvrw_,alfvlw_)
494 if(iw==
mag(idir))
then
495 rq(ix^s)=-q(ix^s)*wroe(ix^s,
mag(jdir))&
496 /sign(wroe(ix^s,
rho_c_),wroe(ix^s,
mag(idim)))
498 rq(ix^s)=+q(ix^s)*wroe(ix^s,
mag(idir))&
499 /sign(wroe(ix^s,
rho_c_),wroe(ix^s,
mag(idim)))
505 end subroutine rtimes2
integer, dimension(:), allocatable, public mag
Indices of the magnetic field.
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.
integer, parameter ndim
Number of spatial dimensions for grid variables.
integer ndir
Number of spatial dimensions (components) for vector variables.
double precision, dimension(:), allocatable entropycoef
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 entropyfix(ixl, il, al, ar, a, smalla)
Magneto-hydrodynamics module.
subroutine, public twofl_get_csound2_c_from_conserved(w, x, ixil, ixol, csound2)
double precision, public twofl_adiab
The adiabatic constant.
subroutine, public twofl_get_pthermal_c(w, x, ixil, ixol, pth)
integer, public e_c_
Index of the energy density (-1 if not present)
integer, dimension(:), allocatable, public mom_c
Indices of the momentum density.
logical, public divbwave
Add divB wave in Roe solver.
integer, public rho_c_
Index of the density (in the w array)
integer, public, protected twofl_eq_energy
double precision, public twofl_gamma
The adiabatic index.
Subroutines for Roe-type Riemann solver for HD.
subroutine, public twofl_roe_init()