MPI-AMRVAC  2.2
The MPI - Adaptive Mesh Refinement - Versatile Advection Code
mod_mhd_roe.t
Go to the documentation of this file.
1 !> Subroutines for Roe-type Riemann solver for HD
2 module mod_mhd_roe
3  use mod_mhd_phys
5 
6  implicit none
7  private
8 
9  integer, parameter :: fastrw_ = 3,fastlw_=4,slowrw_=5,slowlw_=6 ! Characteristic
10  integer, parameter :: entrow_ = 8,diverw_=7,alfvrw_=1,alfvlw_=2 ! waves
11 
12  public :: mhd_roe_init
13 
14 contains
15 
16  subroutine mhd_roe_init()
18  integer :: il
19 
23 
24  nworkroe=15
25  allocate(entropycoef(nw))
26 
27  do il = 1, nw
28  select case(il)
29  case(fastrw_,fastlw_,slowrw_,slowlw_)
30  entropycoef(il) = 0.2d0
31  case(alfvrw_,alfvlw_)
32  entropycoef(il) = 0.4d0
33  case default
34  entropycoef(il) = -1.0d0
35  end select
36  end do
37 
38  end subroutine mhd_roe_init
39 
40  ! Eight-wave MHD Riemann solver. See Powell, Notes on the eigensystem, Gombosi
41  ! Calculate the wroe average of primitive variables in wL and wR, assignment:
42  ! rho -> sqrho, m -> v, e -> p, B_idim -> B_idim, B_idir -> beta_idir
43  ! Calculate also alpha_f,alpha_s,c_f,c_s,csound2,dp,rhodv
44  !
45  ! wL,wR,wroe are all interface centered quantities
46  subroutine mhd_average(wL,wR,x,ix^L,idim,wroe,workroe)
48 
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)
54 
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))
58 
59  end subroutine mhd_average
60 
61  ! Eight-wave MHD Riemann solver. See Powell, Notes on the eigensystem, Gombosi
62  ! Calculate the wroe average of primitive variables in wL and wR, assignment:
63  ! rho -> sqrho, m -> v, e -> p, B_idim -> B_idim, B_idir -> beta_idir
64  ! Calculate also alpha_f,alpha_s,c_f,c_s,csound2,dp,rhodv
65  !
66  ! wL,wR,wroe are all interface centered quantities
67  subroutine average2(wL,wR,x,ix^L,idim,wroe,cfast,cslow,afast,aslow,csound2,dp, &
68  rhodv,tmp)
70 
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, &
75  rhodv,tmp
76 
77  if (ndir==1) call mpistop("MHD with d=11 is the same as HD")
78 
79  !Averaging primitive variables
80  wroe(ix^s,rho_)=half*(wl(ix^s,rho_)+wr(ix^s,rho_))
81  do idir=1,ndir
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)))
84  end do
85  ! Use afast and aslow for pressures pL and pR
86  call mhd_get_pthermal(wl,x,ixg^ll,ix^l,afast)
87  call mhd_get_pthermal(wr,x,ixg^ll,ix^l,aslow)
88 
89  if(mhd_energy) then
90  wroe(ix^s,e_)=half*(afast(ix^s)+aslow(ix^s))
91  ! dp=pR-pL
92  dp(ix^s)=aslow(ix^s)-afast(ix^s)
93  else
94  dp(ix^s)=aslow(ix^s)-afast(ix^s)
95  end if
96 
97  !CONSERVATIVE rho*dv_idim=dm_idim-v_idim*drho
98  rhodv(ix^s)=wr(ix^s,mom(idim))-wl(ix^s,mom(idim))-&
99  wroe(ix^s,mom(idim))*(wr(ix^s,rho_)-wl(ix^s,rho_))
100 
101  !Calculate csound2,cfast,cslow,alphafast and alphaslow
102 
103  ! get csound**2
104  if(mhd_energy) then
105  csound2(ix^s)=mhd_gamma*wroe(ix^s,p_)/wroe(ix^s,rho_)
106  else
107  csound2(ix^s)=mhd_gamma*mhd_adiab*wroe(ix^s,rho_)**(mhd_gamma-one)
108  end if
109 
110  ! aa=B**2/rho+a**2
111  cfast(ix^s)=sum(wroe(ix^s,mag(:))**2,dim=ndim+1)/wroe(ix^s,rho_)+csound2(ix^s)
112 
113  ! cs**2=0.5*(aa+dsqrt(aa**2-4*a**2*(b_i**2/rho)))
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_)))
116 
117  ! cf**2=aa-cs**2
118  cfast(ix^s)=cfast(ix^s)-cslow(ix^s)
119 
120  ! alpha_f**2=(a**2-cs**2)/(cf**2-cs**2)
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))
123 
124  ! alpha_s=dsqrt(1-alpha_f**2)
125  aslow(ix^s)=dsqrt(one-afast(ix^s))
126 
127  ! alpha_f=dsqrt(alpha_f**2)
128  afast(ix^s)=dsqrt(afast(ix^s))
129 
130  ! cf=dsqrt(cf**2)
131  cfast(ix^s)=dsqrt(cfast(ix^s))
132 
133  ! cs=dsqrt(cs**2)
134  cslow(ix^s)=dsqrt(cslow(ix^s))
135 
136  !Replace the primitive variables with more useful quantities:
137  ! rho -> dsqrt(rho)
138  wroe(ix^s,rho_)=dsqrt(wroe(ix^s,rho_))
139 
140  ! Avoid sgn(b_idim)==0
141  where(dabs(wroe(ix^s,mag(idim)))<smalldouble)&
142  wroe(ix^s,mag(idim))=smalldouble
143  ! B_idir,jdir -> beta_idir,jdir
144  idir=idim+1-ndir*(idim/ndir)
145  if(ndir==2)then
146  where(wroe(ix^s,mag(idir))>=zero)
147  wroe(ix^s,mag(idir))=one
148  elsewhere
149  wroe(ix^s,mag(idir))=-one
150  end where
151  else
152  !beta_j=B_j/dsqrt(B_i**2+B_j**2); beta_i=B_i/dsqrt(B_i**2+B_j**2)
153  jdir=idir+1-ndir*(idir/ndir)
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)
158  elsewhere
159  wroe(ix^s,mag(idir))=dsqrt(half)
160  wroe(ix^s,mag(jdir))=dsqrt(half)
161  end where
162  endif
163 
164  end subroutine average2
165 
166  ! Calculate the il-th characteristic speed and the jump in the il-th
167  ! characteristic variable in the idim direction within ixL.
168  ! The eigenvalues and the l=r**(-1) matrix is calculated from wroe.
169  ! jump(il)=Sum_il l(il,iw)*(wR(iw)-wL(iw)), where w are the conservative
170  ! variables. However part of the summation is done in advance and saved into
171  ! bdv,bdb,dp and dv variables. "smalla" contains a lower limit for "a" to be
172  ! used in the entropy fix.
173  !
174  ! All the variables are centered on the cell interface, thus the
175  ! "*C" notation is omitted for sake of brevity.
176  subroutine mhd_get_eigenjump(wL,wR,wroe,x,ix^L,il,idim,smalla,a,jump,workroe)
178 
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
184 
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))
190 
191  end subroutine mhd_get_eigenjump
192 
193  ! Calculate the il-th characteristic speed and the jump in the il-th
194  ! characteristic variable in the idim direction within ixL.
195  ! The eigenvalues and the l=r**(-1) matrix is calculated from wroe.
196  ! jump(il)=Sum_il l(il,iw)*(wR(iw)-wL(iw)), where w are the conservative
197  ! variables. However part of the summation is done in advance and saved into
198  ! bdv,bdb,dp and dv variables. "smalla" contains a lower limit for "a" to be
199  ! used in the entropy fix.
200  !
201  ! All the variables are centered on the cell interface, thus the
202  ! "*C" notation is omitted for sake of brevity.
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)
206  use mod_tvd
207 
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
215 
216  idir=idim+1-ndir*(idim/ndir)
217  jdir=idir+1-ndir*(idir/ndir)
218 
219  if(il==fastrw_)then
220  !Fast and slow waves use bdv=sqrho**2*sign(bx)*(betay*dvy+betaz*dvz)
221  ! bdb=sqrho*a* (betay*dBy+betaz*dBz)
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)))
227 
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_)
232  endif
233 
234  if(il==alfvrw_)then
235  !Alfven waves use bdv=0.5*sqrho**2* (betaz*dvy-betay*dvz)
236  ! bdb=0.5*sqrho*sign(bx)*(betaz*dBy-betay*dBz)
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)))
245  endif
246 
247  select case(il)
248  case(fastrw_)
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)))
253  case(fastlw_)
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)))
258  case(slowrw_)
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)))
263  case(slowlw_)
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)))
268  case(entrow_)
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)
271  case(diverw_)
272  if(divbwave)then
273  a(ix^s)=wroe(ix^s,mom(idim))
274  jump(ix^s)=wr(ix^s,mag(idim))-wl(ix^s,mag(idim))
275  else
276  a(ix^s)=zero
277  jump(ix^s)=zero
278  endif
279  case(alfvrw_)
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)
282  case(alfvlw_)
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)
285  end select
286 
287  ! Calculate "smalla" or modify "a" based on the "typeentropy" switch
288 
289  select case(typeentropy(il))
290  case('yee')
291  ! Based on Yee JCP 68,151 eq 3.23
292  smalla(ix^s)=entropycoef(il)
293  case('harten','powell', 'ratio')
294  ! Based on Harten & Hyman JCP 50, 235 and Zeeuw & Powell JCP 104,56
295  ! Initialize left and right eigenvalues by velocities
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_)
298  ! Calculate the final "aL" and "aR"
299  select case(il)
300  case(fastrw_)
301  ! These quantities will be used for all the fast and slow waves
302  ! Calculate soundspeed**2 and cs**2+ca**2.
303  call mhd_get_csound2(wl,x,ixg^ll,ix^l,cs2l)
304  call mhd_get_csound2(wr,x,ixg^ll,ix^l,cs2r)
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_)
307  ! Save the discriminants into cs2L and cs2R
308  cs2l(ix^s)=&
309  dsqrt(cs2ca2l(ix^s)**2-4d0*cs2l(ix^s)*wl(ix^s,mag(idim))**2/wl(ix^s,rho_))
310  cs2r(ix^s)=&
311  dsqrt(cs2ca2r(ix^s)**2-4d0*cs2r(ix^s)*wr(ix^s,mag(idim))**2/wr(ix^s,rho_))
312 
313  ! The left and right eigenvalues for the fast wave going to right
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)))
316  case(fastlw_)
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)))
319  case(slowrw_)
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)))
322  case(slowlw_)
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_)
326  ! These propagate by the velocity
327  case(alfvrw_)
328  ! Store the Alfven speeds into cs2ca2L and cs2ca2R
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_))
331 
332  al(ix^s)=al(ix^s) + cs2ca2l(ix^s)
333  ar(ix^s)=ar(ix^s) + cs2ca2r(ix^s)
334  case(alfvlw_)
335  al(ix^s)=al(ix^s) - cs2ca2l(ix^s)
336  ar(ix^s)=ar(ix^s) - cs2ca2r(ix^s)
337  end select
338  end select
339 
340  call entropyfix(ix^l,il,al,ar,a,smalla)
341 
342  end subroutine geteigenjump2
343 
344  ! Multiply q by R(il,iw), where R is the right eigenvalue matrix at wroe
345  subroutine mhd_rtimes(q,w,ix^L,iw,il,idim,rq,workroe)
347 
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)
352 
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))
357 
358  end subroutine mhd_rtimes
359 
360  ! Multiply q by R(il,iw), where R is the right eigenvalue matrix at wroe
361  subroutine rtimes2(q,wroe,ix^L,iw,il,idim,rq, &
362  cfast,cslow,afast,aslow,csound2,dp,rhodv,bv,v2a2)
364 
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
370 
371  idir=idim+1-ndir*(idim/ndir)
372  jdir=idir+1-ndir*(idir/ndir)
373 
374  if(iw == rho_) then
375  select case(il)
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)
380  case(entrow_)
381  rq(ix^s)=q(ix^s)
382  case(diverw_,alfvrw_,alfvlw_)
383  rq(ix^s)=zero
384  end select
385  else if(iw == e_) then
386  if(il==fastrw_)then
387  ! Store 0.5*v**2+(2-gamma)/(gamma-1)*a**2
388  v2a2(ix^s)=half*sum(wroe(ix^s,mom(:))**2,dim=ndim+1)+ &
389  (two-mhd_gamma)/(mhd_gamma-one)*csound2(ix^s)
390  ! Store sgn(bx)*(betay*vy+betaz*vz) in bv
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
395  !Store betaz*vy-betay*vz in bv
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)))
398  endif
399 
400  select case(il)
401  case(fastrw_)
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)))))
404  case(fastlw_)
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)))))
407  case(slowrw_)
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)))))
410  case(slowlw_)
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)))))
413  case(entrow_)
414  rq(ix^s)= q(ix^s)*half*sum(wroe(ix^s,mom(:))**2,dim=ndim+1)
415  case(diverw_)
416  if(divbwave)then
417  rq(ix^s)= q(ix^s)*wroe(ix^s,mag(idim))
418  else
419  rq(ix^s)= zero
420  endif
421  case(alfvrw_)
422  rq(ix^s)=+q(ix^s)*bv(ix^s)
423  case(alfvlw_)
424  rq(ix^s)=-q(ix^s)*bv(ix^s)
425  end select
426  else if(any(mom(:)==iw)) then
427  if(iw==mom(idim))then
428  select case(il)
429  case(fastrw_)
430  rq(ix^s)=q(ix^s)*afast(ix^s)*(wroe(ix^s,iw)+cfast(ix^s))
431  case(fastlw_)
432  rq(ix^s)=q(ix^s)*afast(ix^s)*(wroe(ix^s,iw)-cfast(ix^s))
433  case(slowrw_)
434  rq(ix^s)=q(ix^s)*aslow(ix^s)*(wroe(ix^s,iw)+cslow(ix^s))
435  case(slowlw_)
436  rq(ix^s)=q(ix^s)*aslow(ix^s)*(wroe(ix^s,iw)-cslow(ix^s))
437  case(entrow_)
438  rq(ix^s)=q(ix^s)*wroe(ix^s,iw)
439  case(diverw_,alfvlw_,alfvrw_)
440  rq(ix^s)=zero
441  end select
442  else
443  select case(il)
444  case(fastrw_)
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))))
447  case(fastlw_)
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))))
450  case(slowrw_)
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))))
453  case(slowlw_)
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))))
456  case(entrow_)
457  rq(ix^s)=q(ix^s)*wroe(ix^s,iw)
458  case(diverw_)
459  rq(ix^s)=zero
460  case(alfvrw_)
461  if(iw==mom(idir))then
462  rq(ix^s)=+q(ix^s)*wroe(ix^s,mag(jdir))
463  else
464  rq(ix^s)=-q(ix^s)*wroe(ix^s,mag(idir))
465  endif
466  case(alfvlw_)
467  if(iw==mom(idir))then
468  rq(ix^s)=-q(ix^s)*wroe(ix^s,mag(jdir))
469  else
470  rq(ix^s)=+q(ix^s)*wroe(ix^s,mag(idir))
471  endif
472  end select
473  end if ! iw=m_idir,m_jdir
474  else if(any(mag(:)==iw)) then
475  if(iw==mag(idim))then
476  if(il==diverw_ .and. divbwave)then
477  rq(ix^s)=q(ix^s)
478  else
479  rq(ix^s)=zero
480  endif
481  else
482  select case(il)
483  case(fastrw_,fastlw_)
484  rq(ix^s)=+q(ix^s)*aslow(ix^s)*dsqrt(csound2(ix^s))*wroe(ix^s,iw)&
485  /wroe(ix^s,rho_)
486  case(slowrw_,slowlw_)
487  rq(ix^s)=-q(ix^s)*afast(ix^s)*dsqrt(csound2(ix^s))*wroe(ix^s,iw)&
488  /wroe(ix^s,rho_)
489  case(entrow_,diverw_)
490  rq(ix^s)=zero
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)))
495  else
496  rq(ix^s)=+q(ix^s)*wroe(ix^s,mag(idir))&
497  /sign(wroe(ix^s,rho_),wroe(ix^s,mag(idim)))
498  end if
499  end select
500  end if ! iw=b_idir,b_jdir
501  end if
502 
503  end subroutine rtimes2
504 
505 end module mod_mhd_roe
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
integer, public, protected e_
Index of the energy density (-1 if not present)
Definition: mod_mhd_phys.t:66
integer, public, protected p_
Index of the gas pressure (-1 if not present) should equal e_.
Definition: mod_mhd_phys.t:69
procedure(sub_rtimes), pointer phys_rtimes
procedure(sub_average), pointer phys_average
Subroutines for Roe-type Riemann solver for HD.
Definition: mod_mhd_roe.t:2
integer ndir
Number of spatial dimensions (components) for vector variables.
subroutine mhd_rtimes(q, w, ixL, iw, il, idim, rq, workroe)
Definition: mod_mhd_roe.t:346
double precision, dimension(:), allocatable entropycoef
logical, public, protected mhd_energy
Whether an energy equation is used.
Definition: mod_mhd_phys.t:8
subroutine rtimes2(q, wroe, ixL, iw, il, idim, rq, cfast, cslow, afast, aslow, csound2, dp, rhodv, bv, v2a2)
Definition: mod_mhd_roe.t:363
double precision, public mhd_gamma
The adiabatic index.
Definition: mod_mhd_phys.t:81
subroutine mhd_average(wL, wR, x, ixL, idim, wroe, workroe)
Definition: mod_mhd_roe.t:47
integer, dimension(:), allocatable, public, protected mag
Indices of the magnetic field.
Definition: mod_mhd_phys.t:72
character(len=std_len), dimension(:), allocatable typeentropy
Which type of entropy fix to use with Riemann-type solvers.
integer, dimension(:), allocatable, public, protected mom
Indices of the momentum density.
Definition: mod_mhd_phys.t:63
subroutine, public mhd_roe_init()
Definition: mod_mhd_roe.t:17
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)
Definition: mod_mhd_roe.t:205
integer, public, protected rho_
Index of the density (in the w array)
Definition: mod_mhd_phys.t:60
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.
Definition: mod_mhd_phys.t:942
subroutine mpistop(message)
Exit MPI-AMRVAC with an error message.
Definition: comm_lib.t:198
subroutine average2(wL, wR, x, ixL, idim, wroe, cfast, cslow, afast, aslow, csound2, dp, rhodv, tmp)
Definition: mod_mhd_roe.t:69
procedure(sub_get_eigenjump), pointer phys_get_eigenjump
integer, parameter ndim
Number of spatial dimensions for grid variables.
subroutine, public mhd_get_pthermal(w, x, ixIL, ixOL, pth)
Calculate thermal pressure=(gamma-1)*(e-0.5*m**2/rho-b**2/2) within ixO^L.
Definition: mod_mhd_phys.t:919
logical, public divbwave
Add divB wave in Roe solver.
Definition: mod_mhd_phys.t:123
Subroutines for TVD-MUSCL schemes.
Definition: mod_tvd.t:2
Magneto-hydrodynamics module.
Definition: mod_mhd_phys.t:2
double precision, public mhd_adiab
The adiabatic constant.
Definition: mod_mhd_phys.t:84
subroutine mhd_get_eigenjump(wL, wR, wroe, x, ixL, il, idim, smalla, a, jump, workroe)
Definition: mod_mhd_roe.t:177