9 integer,
parameter :: fastRW_ = 3,fastlw_=4,slowrw_=5,slowlw_=6
10 integer,
parameter :: entroW_ = 8,diverw_=7,alfvrw_=1,alfvlw_=2
29 case(fastrw_,fastlw_,slowrw_,slowlw_)
49 integer,
intent(in) :: ix^L, idim
50 double precision,
intent(in) :: wL(ixG^T, nw), wR(ixG^T, nw)
51 double precision,
intent(inout) :: wroe(ixG^T, nw)
52 double precision,
intent(inout) :: workroe(ixG^T, nworkroe)
53 double precision,
intent(in) :: x(ixG^T, 1:^ND)
55 call average2(wl,wr,x,ix^l,idim,wroe,workroe(ixg^t,1),workroe(ixg^t,2), &
56 workroe(ixg^t,3),workroe(ixg^t,4),workroe(ixg^t,5),workroe(ixg^t,6), &
57 workroe(ixg^t,7),workroe(ixg^t,8))
67 subroutine average2(wL,wR,x,ix^L,idim,wroe,cfast,cslow,afast,aslow,csound2,dp, &
71 integer :: ix^L,idim,idir,jdir,iw
72 double precision,
dimension(ixG^T,nw) :: wL,wR,wroe
73 double precision,
intent(in) :: x(ixG^T,1:ndim)
74 double precision,
dimension(ixG^T) :: cfast,cslow,afast,aslow,csound2,dp, &
77 if (
ndir==1)
call mpistop(
"MHD with d=11 is the same as HD")
82 wroe(ix^s,
mom(idir))=half*(wl(ix^s,
mom(idir))/wl(ix^s,
rho_)+wr(ix^s,
mom(idir))/wr(ix^s,
rho_))
83 wroe(ix^s,
mag(idir))=half*(wl(ix^s,
mag(idir))+wr(ix^s,
mag(idir)))
90 wroe(ix^s,
e_)=half*(afast(ix^s)+aslow(ix^s))
92 dp(ix^s)=aslow(ix^s)-afast(ix^s)
94 dp(ix^s)=aslow(ix^s)-afast(ix^s)
98 rhodv(ix^s)=wr(ix^s,
mom(idim))-wl(ix^s,
mom(idim))-&
111 cfast(ix^s)=sum(wroe(ix^s,
mag(:))**2,dim=ndim+1)/wroe(ix^s,
rho_)+csound2(ix^s)
114 cslow(ix^s)=half*(cfast(ix^s)-dsqrt(cfast(ix^s)**2-&
115 4d0*csound2(ix^s)*wroe(ix^s,
mag(idim))**2/wroe(ix^s,
rho_)))
118 cfast(ix^s)=cfast(ix^s)-cslow(ix^s)
121 afast(ix^s)=(csound2(ix^s)-cslow(ix^s))/(cfast(ix^s)-cslow(ix^s))
122 afast(ix^s)=min(one,max(afast(ix^s),zero))
125 aslow(ix^s)=dsqrt(one-afast(ix^s))
128 afast(ix^s)=dsqrt(afast(ix^s))
131 cfast(ix^s)=dsqrt(cfast(ix^s))
134 cslow(ix^s)=dsqrt(cslow(ix^s))
138 wroe(ix^s,
rho_)=dsqrt(wroe(ix^s,
rho_))
141 where(dabs(wroe(ix^s,
mag(idim)))<smalldouble)&
142 wroe(ix^s,
mag(idim))=smalldouble
146 where(wroe(ix^s,
mag(idir))>=zero)
147 wroe(ix^s,
mag(idir))=one
149 wroe(ix^s,
mag(idir))=-one
154 tmp(ix^s)=dsqrt(wroe(ix^s,
mag(idir))**2+wroe(ix^s,
mag(jdir))**2)
155 where(tmp(ix^s)>smalldouble)
156 wroe(ix^s,
mag(idir))=wroe(ix^s,
mag(idir))/tmp(ix^s)
157 wroe(ix^s,
mag(jdir))=wroe(ix^s,
mag(jdir))/tmp(ix^s)
159 wroe(ix^s,
mag(idir))=dsqrt(half)
160 wroe(ix^s,
mag(jdir))=dsqrt(half)
176 subroutine mhd_get_eigenjump(wL,wR,wroe,x,ix^L,il,idim,smalla,a,jump,workroe)
179 integer,
intent(in) :: ix^L,il,idim
180 double precision,
dimension(ixG^T,nw) :: wL,wR,wroe
181 double precision,
intent(in) :: x(ixG^T,1:ndim)
182 double precision,
dimension(ixG^T) :: smalla,a,jump
183 double precision,
dimension(ixG^T,nworkroe) :: workroe
185 call geteigenjump2(wl,wr,wroe,x,ix^l,il,idim,smalla,a,jump, &
186 workroe(ixg^t,1),workroe(ixg^t,2), &
187 workroe(ixg^t,3),workroe(ixg^t,4),workroe(ixg^t,5),workroe(ixg^t,6), &
188 workroe(ixg^t,7),workroe(ixg^t,8),workroe(ixg^t,9),workroe(ixg^t,10), &
189 workroe(ixg^t,11),workroe(ixg^t,12),workroe(ixg^t,13))
203 subroutine geteigenjump2(wL,wR,wroe,x,ix^L,il,idim,smalla,a,jump, &
204 cfast,cslow,afast,aslow,csound2,dp,rhodv,bdv,bdb,cs2L,cs2R,cs2ca2L,cs2ca2R)
208 integer :: ix^L,il,idim,idir,jdir
209 double precision,
dimension(ixG^T,nw) :: wL,wR,wroe
210 double precision,
intent(in) :: x(ixG^T,1:ndim)
211 double precision,
dimension(ixG^T) :: smalla,a,jump
212 double precision,
dimension(ixG^T) :: cfast,cslow,afast,aslow,csound2,dp,rhodv
213 double precision,
dimension(ixG^T) :: bdv,bdb
214 double precision,
dimension(ixG^T) :: aL,aR,cs2L,cs2R,cs2ca2L,cs2ca2R
222 bdv(ix^s)=wroe(ix^s,
mag(idir))* &
223 (wr(ix^s,
mom(idir))/wr(ix^s,
rho_)-wl(ix^s,
mom(idir))/wl(ix^s,
rho_))
224 if(
ndir==3)bdv(ix^s)=bdv(ix^s)+wroe(ix^s,
mag(jdir))* &
225 (wr(ix^s,
mom(jdir))/wr(ix^s,
rho_)-wl(ix^s,
mom(jdir))/wl(ix^s,
rho_))
226 bdv(ix^s)=bdv(ix^s)*sign(wroe(ix^s,
rho_)**2,wroe(ix^s,
mag(idim)))
228 bdb(ix^s)=wroe(ix^s,
mag(idir))*(wr(ix^s,
mag(idir))-wl(ix^s,
mag(idir)))
229 if(
ndir==3)bdb(ix^s)=bdb(ix^s)+&
230 wroe(ix^s,
mag(jdir))*(wr(ix^s,
mag(jdir))-wl(ix^s,
mag(jdir)))
231 bdb(ix^s)=bdb(ix^s)*dsqrt(csound2(ix^s))*wroe(ix^s,
rho_)
237 bdv(ix^s)=wroe(ix^s,
mag(jdir))* &
238 (wr(ix^s,
mom(idir))/wr(ix^s,
rho_)-wl(ix^s,
mom(idir))/wl(ix^s,
rho_)) &
239 -wroe(ix^s,
mag(idir))* &
240 (wr(ix^s,
mom(jdir))/wr(ix^s,
rho_)-wl(ix^s,
mom(jdir))/wl(ix^s,
rho_))
241 bdb(ix^s)=wroe(ix^s,
mag(jdir))*(wr(ix^s,
mag(idir))-wl(ix^s,
mag(idir))) &
242 -wroe(ix^s,
mag(idir))*(wr(ix^s,
mag(jdir))-wl(ix^s,
mag(jdir)))
243 bdv(ix^s)=bdv(ix^s)*half*wroe(ix^s,
rho_)**2
244 bdb(ix^s)=bdb(ix^s)*half*sign(wroe(ix^s,
rho_),wroe(ix^s,
mag(idim)))
249 a(ix^s)=wroe(ix^s,
mom(idim))+cfast(ix^s)
250 jump(ix^s)=half/csound2(ix^s)*(&
251 afast(ix^s)*(+cfast(ix^s)*rhodv(ix^s)+dp(ix^s))&
252 +aslow(ix^s)*(-cslow(ix^s)*bdv(ix^s)+bdb(ix^s)))
254 a(ix^s)=wroe(ix^s,
mom(idim))-cfast(ix^s)
255 jump(ix^s)=half/csound2(ix^s)*(&
256 afast(ix^s)*(-cfast(ix^s)*rhodv(ix^s)+dp(ix^s))&
257 +aslow(ix^s)*(+cslow(ix^s)*bdv(ix^s)+bdb(ix^s)))
259 a(ix^s)=wroe(ix^s,
mom(idim))+cslow(ix^s)
260 jump(ix^s)=half/csound2(ix^s)*(&
261 aslow(ix^s)*(+cslow(ix^s)*rhodv(ix^s)+dp(ix^s))&
262 +afast(ix^s)*(+cfast(ix^s)*bdv(ix^s)-bdb(ix^s)))
264 a(ix^s)=wroe(ix^s,
mom(idim))-cslow(ix^s)
265 jump(ix^s)=half/csound2(ix^s)*(&
266 aslow(ix^s)*(-cslow(ix^s)*rhodv(ix^s)+dp(ix^s))&
267 +afast(ix^s)*(-cfast(ix^s)*bdv(ix^s)-bdb(ix^s)))
269 a(ix^s)=wroe(ix^s,
mom(idim))
270 jump(ix^s)=wr(ix^s,
rho_)-wl(ix^s,
rho_)-dp(ix^s)/csound2(ix^s)
273 a(ix^s)=wroe(ix^s,
mom(idim))
274 jump(ix^s)=wr(ix^s,
mag(idim))-wl(ix^s,
mag(idim))
280 a(ix^s)=wroe(ix^s,
mom(idim))+dabs(wroe(ix^s,
mag(idim)))/wroe(ix^s,
rho_)
281 jump(ix^s)=+bdv(ix^s)-bdb(ix^s)
283 a(ix^s)=wroe(ix^s,
mom(idim))-dabs(wroe(ix^s,
mag(idim)))/wroe(ix^s,
rho_)
284 jump(ix^s)=-bdv(ix^s)-bdb(ix^s)
293 case(
'harten',
'powell',
'ratio')
296 al(ix^s)= wl(ix^s,
mom(idim))/wl(ix^s,
rho_)
297 ar(ix^s)= wr(ix^s,
mom(idim))/wr(ix^s,
rho_)
305 cs2ca2l(ix^s)=cs2l(ix^s)+sum(wl(ix^s,
mag(:))**2,dim=ndim+1)/wl(ix^s,
rho_)
306 cs2ca2r(ix^s)=cs2r(ix^s)+sum(wr(ix^s,
mag(:))**2,dim=ndim+1)/wr(ix^s,
rho_)
309 dsqrt(cs2ca2l(ix^s)**2-4d0*cs2l(ix^s)*wl(ix^s,
mag(idim))**2/wl(ix^s,
rho_))
311 dsqrt(cs2ca2r(ix^s)**2-4d0*cs2r(ix^s)*wr(ix^s,
mag(idim))**2/wr(ix^s,
rho_))
314 al(ix^s)=al(ix^s) + dsqrt(half*(cs2ca2l(ix^s) + cs2l(ix^s)))
315 ar(ix^s)=ar(ix^s) + dsqrt(half*(cs2ca2r(ix^s) + cs2r(ix^s)))
317 al(ix^s)=al(ix^s) - dsqrt(half*(cs2ca2l(ix^s) + cs2l(ix^s)))
318 ar(ix^s)=ar(ix^s) - dsqrt(half*(cs2ca2r(ix^s) + cs2r(ix^s)))
320 al(ix^s)=al(ix^s) + dsqrt(half*(cs2ca2l(ix^s) - cs2l(ix^s)))
321 ar(ix^s)=ar(ix^s) + dsqrt(half*(cs2ca2r(ix^s) - cs2r(ix^s)))
323 al(ix^s)=al(ix^s) - dsqrt(half*(cs2ca2l(ix^s) - cs2l(ix^s)))
324 ar(ix^s)=ar(ix^s) - dsqrt(half*(cs2ca2r(ix^s) - cs2r(ix^s)))
325 case(entrow_,diverw_)
329 cs2ca2l(ix^s)=dabs(wl(ix^s,
mag(idim)))/dsqrt(wl(ix^s,
rho_))
330 cs2ca2r(ix^s)=dabs(wr(ix^s,
mag(idim)))/dsqrt(wr(ix^s,
rho_))
332 al(ix^s)=al(ix^s) + cs2ca2l(ix^s)
333 ar(ix^s)=ar(ix^s) + cs2ca2r(ix^s)
335 al(ix^s)=al(ix^s) - cs2ca2l(ix^s)
336 ar(ix^s)=ar(ix^s) - cs2ca2r(ix^s)
348 integer,
intent(in) :: ix^L, iw, il, idim
349 double precision,
intent(in) :: w(ixG^T, nw), q(ixG^T)
350 double precision,
intent(inout) :: rq(ixG^T)
351 double precision,
intent(inout) :: workroe(ixG^T, nworkroe)
353 call rtimes2(q,w,ix^l,iw,il,idim,rq,&
354 workroe(ixg^t,1),workroe(ixg^t,2), &
355 workroe(ixg^t,3),workroe(ixg^t,4),workroe(ixg^t,5),workroe(ixg^t,6), &
356 workroe(ixg^t,7),workroe(ixg^t,14),workroe(ixg^t,15))
361 subroutine rtimes2(q,wroe,ix^L,iw,il,idim,rq, &
362 cfast,cslow,afast,aslow,csound2,dp,rhodv,bv,v2a2)
365 integer :: ix^L,iw,il,idim,idir,jdir
366 double precision :: wroe(ixG^T,nw)
367 double precision,
dimension(ixG^T) :: q,rq
368 double precision,
dimension(ixG^T) :: cfast,cslow,afast,aslow,csound2,dp,rhodv
369 double precision,
dimension(ixG^T) :: bv,v2a2
376 case(fastrw_,fastlw_)
377 rq(ix^s)=q(ix^s)*afast(ix^s)
378 case(slowrw_,slowlw_)
379 rq(ix^s)=q(ix^s)*aslow(ix^s)
382 case(diverw_,alfvrw_,alfvlw_)
385 else if(iw ==
e_)
then
388 v2a2(ix^s)=half*sum(wroe(ix^s,
mom(:))**2,dim=
ndim+1)+ &
391 bv(ix^s)=wroe(ix^s,
mag(idir))*wroe(ix^s,
mom(idir))
392 if(
ndir==3)bv(ix^s)=bv(ix^s)+wroe(ix^s,
mag(jdir))*wroe(ix^s,
mom(jdir))
393 bv(ix^s)=bv(ix^s)*sign(one,wroe(ix^s,
mag(idim)))
394 else if(il==alfvrw_)
then
396 bv(ix^s)=(wroe(ix^s,
mag(jdir))*wroe(ix^s,
mom(idir))-&
397 wroe(ix^s,
mag(idir))*wroe(ix^s,
mom(jdir)))
402 rq(ix^s)=q(ix^s)*(-aslow(ix^s)*cslow(ix^s)*bv(ix^s)+afast(ix^s)*&
403 (v2a2(ix^s)+cfast(ix^s)*(cfast(ix^s)+wroe(ix^s,
mom(idim)))))
405 rq(ix^s)=q(ix^s)*(+aslow(ix^s)*cslow(ix^s)*bv(ix^s)+afast(ix^s)*&
406 (v2a2(ix^s)+cfast(ix^s)*(cfast(ix^s)-wroe(ix^s,
mom(idim)))))
408 rq(ix^s)=q(ix^s)*(+afast(ix^s)*cfast(ix^s)*bv(ix^s)+aslow(ix^s)*&
409 (v2a2(ix^s)+cslow(ix^s)*(cslow(ix^s)+wroe(ix^s,
mom(idim)))))
411 rq(ix^s)=q(ix^s)*(-afast(ix^s)*cfast(ix^s)*bv(ix^s)+aslow(ix^s)*&
412 (v2a2(ix^s)+cslow(ix^s)*(cslow(ix^s)-wroe(ix^s,
mom(idim)))))
414 rq(ix^s)= q(ix^s)*half*sum(wroe(ix^s,
mom(:))**2,dim=
ndim+1)
417 rq(ix^s)= q(ix^s)*wroe(ix^s,
mag(idim))
422 rq(ix^s)=+q(ix^s)*bv(ix^s)
424 rq(ix^s)=-q(ix^s)*bv(ix^s)
426 else if(any(
mom(:)==iw))
then
427 if(iw==
mom(idim))
then
430 rq(ix^s)=q(ix^s)*afast(ix^s)*(wroe(ix^s,iw)+cfast(ix^s))
432 rq(ix^s)=q(ix^s)*afast(ix^s)*(wroe(ix^s,iw)-cfast(ix^s))
434 rq(ix^s)=q(ix^s)*aslow(ix^s)*(wroe(ix^s,iw)+cslow(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)*wroe(ix^s,iw)
439 case(diverw_,alfvlw_,alfvrw_)
445 rq(ix^s)=q(ix^s)*(afast(ix^s)*wroe(ix^s,iw)-aslow(ix^s)*&
446 cslow(ix^s)*wroe(ix^s,
mag(1)-
mom(1)+iw)*sign(one,wroe(ix^s,
mag(idim))))
448 rq(ix^s)=q(ix^s)*(afast(ix^s)*wroe(ix^s,iw)+aslow(ix^s)*&
449 cslow(ix^s)*wroe(ix^s,
mag(1)-
mom(1)+iw)*sign(one,wroe(ix^s,
mag(idim))))
451 rq(ix^s)=q(ix^s)*(aslow(ix^s)*wroe(ix^s,iw)+afast(ix^s)*&
452 cfast(ix^s)*wroe(ix^s,
mag(1)-
mom(1)+iw)*sign(one,wroe(ix^s,
mag(idim))))
454 rq(ix^s)=q(ix^s)*(aslow(ix^s)*wroe(ix^s,iw)-afast(ix^s)*&
455 cfast(ix^s)*wroe(ix^s,
mag(1)-
mom(1)+iw)*sign(one,wroe(ix^s,
mag(idim))))
457 rq(ix^s)=q(ix^s)*wroe(ix^s,iw)
461 if(iw==
mom(idir))
then
462 rq(ix^s)=+q(ix^s)*wroe(ix^s,
mag(jdir))
464 rq(ix^s)=-q(ix^s)*wroe(ix^s,
mag(idir))
467 if(iw==
mom(idir))
then
468 rq(ix^s)=-q(ix^s)*wroe(ix^s,
mag(jdir))
470 rq(ix^s)=+q(ix^s)*wroe(ix^s,
mag(idir))
474 else if(any(
mag(:)==iw))
then
475 if(iw==
mag(idim))
then
483 case(fastrw_,fastlw_)
484 rq(ix^s)=+q(ix^s)*aslow(ix^s)*dsqrt(csound2(ix^s))*wroe(ix^s,iw)&
486 case(slowrw_,slowlw_)
487 rq(ix^s)=-q(ix^s)*afast(ix^s)*dsqrt(csound2(ix^s))*wroe(ix^s,iw)&
489 case(entrow_,diverw_)
491 case(alfvrw_,alfvlw_)
492 if(iw==
mag(idir))
then
493 rq(ix^s)=-q(ix^s)*wroe(ix^s,
mag(jdir))&
494 /sign(wroe(ix^s,
rho_),wroe(ix^s,
mag(idim)))
496 rq(ix^s)=+q(ix^s)*wroe(ix^s,
mag(idir))&
497 /sign(wroe(ix^s,
rho_),wroe(ix^s,
mag(idim)))
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.
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
Magneto-hydrodynamics module.
integer, dimension(:), allocatable, public, protected mom
Indices of the momentum density.
double precision, public mhd_adiab
The adiabatic constant.
double precision, public mhd_gamma
The adiabatic index.
integer, dimension(:), allocatable, public, protected mag
Indices of the magnetic field.
logical, public, protected mhd_energy
Whether an energy equation is used.
logical, public divbwave
Add divB wave in Roe solver.
subroutine, public mhd_get_csound2(w, x, ixIL, ixOL, csound2)
Calculate the square of the thermal sound speed csound2 within ixO^L. csound2=gamma*p/rho.
integer, public, protected p_
Index of the gas pressure (-1 if not present) should equal e_.
integer, public, protected rho_
Index of the density (in the w array)
procedure(sub_get_pthermal), pointer, public mhd_get_pthermal
integer, public, protected e_
Index of the energy density (-1 if not present)
Subroutines for Roe-type Riemann solver for MHD.
subroutine average2(wL, wR, x, ixL, idim, wroe, cfast, cslow, afast, aslow, csound2, dp, rhodv, tmp)
subroutine mhd_get_eigenjump(wL, wR, wroe, x, ixL, il, idim, smalla, a, jump, workroe)
subroutine, public mhd_roe_init()
subroutine rtimes2(q, wroe, ixL, iw, il, idim, rq, cfast, cslow, afast, aslow, csound2, dp, rhodv, bv, v2a2)
subroutine mhd_average(wL, wR, x, ixL, idim, wroe, workroe)
subroutine geteigenjump2(wL, wR, wroe, x, ixL, il, idim, smalla, a, jump, cfast, cslow, afast, aslow, csound2, dp, rhodv, bdv, bdb, cs2L, cs2R, cs2ca2L, cs2ca2R)
subroutine mhd_rtimes(q, w, ixL, iw, il, idim, rq, workroe)
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)