MPI-AMRVAC  3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
mod_twofl_roe.t
Go to the documentation of this file.
1 !> Subroutines for Roe-type Riemann solver for HD
3 #include "amrvac.h"
4 
6  use mod_functions_bfield, only: mag
8 
9  implicit none
10  private
11 
12  integer, parameter :: fastRW_ = 3,fastlw_=4,slowrw_=5,slowlw_=6 ! Characteristic
13  integer, parameter :: entroW_ = 8,diverw_=7,alfvrw_=1,alfvlw_=2 ! waves
14 
15  public :: twofl_roe_init
16 
17 contains
18 
19  subroutine twofl_roe_init()
21  integer :: il
22 
26 
27  nworkroe=15
28  allocate(entropycoef(nw))
29 
30  do il = 1, nw
31  select case(il)
32  case(fastrw_,fastlw_,slowrw_,slowlw_)
33  entropycoef(il) = 0.2d0
34  case(alfvrw_,alfvlw_)
35  entropycoef(il) = 0.4d0
36  case default
37  entropycoef(il) = -1.0d0
38  end select
39  end do
40 
41  end subroutine twofl_roe_init
42 
43  ! Eight-wave MHD Riemann solver. See Powell, Notes on the eigensystem, Gombosi
44  ! Calculate the wroe average of primitive variables in wL and wR, assignment:
45  ! rho -> sqrho, m -> v, e -> p, B_idim -> B_idim, B_idir -> beta_idir
46  ! Calculate also alpha_f,alpha_s,c_f,c_s,csound2,dp,rhodv
47  !
48  ! wL,wR,wroe are all interface centered quantities
49  subroutine twofl_average(wL,wR,x,ix^L,idim,wroe,workroe)
51 
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)
57 
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))
61 
62  end subroutine twofl_average
63 
64  ! Eight-wave MHD Riemann solver. See Powell, Notes on the eigensystem, Gombosi
65  ! Calculate the wroe average of primitive variables in wL and wR, assignment:
66  ! rho -> sqrho, m -> v, e -> p, B_idim -> B_idim, B_idir -> beta_idir
67  ! Calculate also alpha_f,alpha_s,c_f,c_s,csound2,dp,rhodv
68  !
69  ! wL,wR,wroe are all interface centered quantities
70  subroutine average2(wL,wR,x,ix^L,idim,wroe,cfast,cslow,afast,aslow,csound2,dp, &
71  rhodv,tmp)
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, &
77  rhodv,tmp
78 
79  if (ndir==1) call mpistop("twofl with d=11 is the same as HD")
80 
81  !Averaging primitive variables
82  wroe(ix^s,rho_c_)=half*(wl(ix^s,rho_c_)+wr(ix^s,rho_c_))
83  do idir=1,ndir
84  wroe(ix^s,mom_c(idir))=half*(wl(ix^s,mom_c(idir))/wl(ix^s,rho_c_)+wr(ix^s,mom_c(idir))/wr(ix^s,rho_c_))
85  wroe(ix^s,mag(idir))=half*(wl(ix^s,mag(idir))+wr(ix^s,mag(idir)))
86  end do
87  ! Use afast and aslow for pressures pL and pR
88  call twofl_get_pthermal_c(wl,x,ixg^ll,ix^l,afast)
89  call twofl_get_pthermal_c(wr,x,ixg^ll,ix^l,aslow)
90 
91  if(twofl_eq_energy/=0) then
92  wroe(ix^s,e_c_)=half*(afast(ix^s)+aslow(ix^s))
93  ! dp=pR-pL
94  dp(ix^s)=aslow(ix^s)-afast(ix^s)
95  else
96  dp(ix^s)=aslow(ix^s)-afast(ix^s)
97  end if
98 
99  !CONSERVATIVE rho*dv_idim=dm_idim-v_idim*drho
100  rhodv(ix^s)=wr(ix^s,mom_c(idim))-wl(ix^s,mom_c(idim))-&
101  wroe(ix^s,mom_c(idim))*(wr(ix^s,rho_c_)-wl(ix^s,rho_c_))
102 
103  !Calculate csound2,cfast,cslow,alphafast and alphaslow
104 
105  ! get csound**2
106  if(twofl_eq_energy/=0) then
107  csound2(ix^s)=twofl_gamma*wroe(ix^s,e_c_)/wroe(ix^s,rho_c_)
108  else
109  csound2(ix^s)=twofl_gamma*twofl_adiab*wroe(ix^s,rho_c_)**(twofl_gamma-one)
110  end if
111 
112  ! aa=B**2/rho+a**2
113  cfast(ix^s)=sum(wroe(ix^s,mag(:))**2,dim=ndim+1)/wroe(ix^s,rho_c_)+csound2(ix^s)
114 
115  ! cs**2=0.5*(aa+dsqrt(aa**2-4*a**2*(b_i**2/rho)))
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_)))
118 
119  ! cf**2=aa-cs**2
120  cfast(ix^s)=cfast(ix^s)-cslow(ix^s)
121 
122  ! alpha_f**2=(a**2-cs**2)/(cf**2-cs**2)
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))
125 
126  ! alpha_s=dsqrt(1-alpha_f**2)
127  aslow(ix^s)=dsqrt(one-afast(ix^s))
128 
129  ! alpha_f=dsqrt(alpha_f**2)
130  afast(ix^s)=dsqrt(afast(ix^s))
131 
132  ! cf=dsqrt(cf**2)
133  cfast(ix^s)=dsqrt(cfast(ix^s))
134 
135  ! cs=dsqrt(cs**2)
136  cslow(ix^s)=dsqrt(cslow(ix^s))
137 
138  !Replace the primitive variables with more useful quantities:
139  ! rho -> dsqrt(rho)
140  wroe(ix^s,rho_c_)=dsqrt(wroe(ix^s,rho_c_))
141 
142  ! Avoid sgn(b_idim)==0
143  where(dabs(wroe(ix^s,mag(idim)))<smalldouble)&
144  wroe(ix^s,mag(idim))=smalldouble
145  ! B_idir,jdir -> beta_idir,jdir
146  idir=idim+1-ndir*(idim/ndir)
147  if(ndir==2)then
148  where(wroe(ix^s,mag(idir))>=zero)
149  wroe(ix^s,mag(idir))=one
150  elsewhere
151  wroe(ix^s,mag(idir))=-one
152  end where
153  else
154  !beta_j=B_j/dsqrt(B_i**2+B_j**2); beta_i=B_i/dsqrt(B_i**2+B_j**2)
155  jdir=idir+1-ndir*(idir/ndir)
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)
160  elsewhere
161  wroe(ix^s,mag(idir))=dsqrt(half)
162  wroe(ix^s,mag(jdir))=dsqrt(half)
163  end where
164  endif
165 
166  end subroutine average2
167 
168  ! Calculate the il-th characteristic speed and the jump in the il-th
169  ! characteristic variable in the idim direction within ixL.
170  ! The eigenvalues and the l=r**(-1) matrix is calculated from wroe.
171  ! jump(il)=Sum_il l(il,iw)*(wR(iw)-wL(iw)), where w are the conservative
172  ! variables. However part of the summation is done in advance and saved into
173  ! bdv,bdb,dp and dv variables. "smalla" contains a lower limit for "a" to be
174  ! used in the entropy fix.
175  !
176  ! All the variables are centered on the cell interface, thus the
177  ! "*C" notation is omitted for sake of brevity.
178  subroutine twofl_get_eigenjump(wL,wR,wroe,x,ix^L,il,idim,smalla,a,jump,workroe)
180 
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
186 
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))
192 
193  end subroutine twofl_get_eigenjump
194 
195  ! Calculate the il-th characteristic speed and the jump in the il-th
196  ! characteristic variable in the idim direction within ixL.
197  ! The eigenvalues and the l=r**(-1) matrix is calculated from wroe.
198  ! jump(il)=Sum_il l(il,iw)*(wR(iw)-wL(iw)), where w are the conservative
199  ! variables. However part of the summation is done in advance and saved into
200  ! bdv,bdb,dp and dv variables. "smalla" contains a lower limit for "a" to be
201  ! used in the entropy fix.
202  !
203  ! All the variables are centered on the cell interface, thus the
204  ! "*C" notation is omitted for sake of brevity.
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)
208  use mod_tvd
209 
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
217 
218  idir=idim+1-ndir*(idim/ndir)
219  jdir=idir+1-ndir*(idir/ndir)
220 
221  if(il==fastrw_)then
222  !Fast and slow waves use bdv=sqrho**2*sign(bx)*(betay*dvy+betaz*dvz)
223  ! bdb=sqrho*a* (betay*dBy+betaz*dBz)
224  bdv(ix^s)=wroe(ix^s,mag(idir))* &
225  (wr(ix^s,mom_c(idir))/wr(ix^s,rho_c_)-wl(ix^s,mom_c(idir))/wl(ix^s,rho_c_))
226  if(ndir==3)bdv(ix^s)=bdv(ix^s)+wroe(ix^s,mag(jdir))* &
227  (wr(ix^s,mom_c(jdir))/wr(ix^s,rho_c_)-wl(ix^s,mom_c(jdir))/wl(ix^s,rho_c_))
228  bdv(ix^s)=bdv(ix^s)*sign(wroe(ix^s,rho_c_)**2,wroe(ix^s,mag(idim)))
229 
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_)
234  endif
235 
236  if(il==alfvrw_)then
237  !Alfven waves use bdv=0.5*sqrho**2* (betaz*dvy-betay*dvz)
238  ! bdb=0.5*sqrho*sign(bx)*(betaz*dBy-betay*dBz)
239  bdv(ix^s)=wroe(ix^s,mag(jdir))* &
240  (wr(ix^s,mom_c(idir))/wr(ix^s,rho_c_)-wl(ix^s,mom_c(idir))/wl(ix^s,rho_c_)) &
241  -wroe(ix^s,mag(idir))* &
242  (wr(ix^s,mom_c(jdir))/wr(ix^s,rho_c_)-wl(ix^s,mom_c(jdir))/wl(ix^s,rho_c_))
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)))
247  endif
248 
249  select case(il)
250  case(fastrw_)
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)))
255  case(fastlw_)
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)))
260  case(slowrw_)
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)))
265  case(slowlw_)
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)))
270  case(entrow_)
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)
273  case(diverw_)
274  if(divbwave)then
275  a(ix^s)=wroe(ix^s,mom_c(idim))
276  jump(ix^s)=wr(ix^s,mag(idim))-wl(ix^s,mag(idim))
277  else
278  a(ix^s)=zero
279  jump(ix^s)=zero
280  endif
281  case(alfvrw_)
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)
284  case(alfvlw_)
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)
287  end select
288 
289  ! Calculate "smalla" or modify "a" based on the "typeentropy" switch
290 
291  select case(typeentropy(il))
292  case('yee')
293  ! Based on Yee JCP 68,151 eq 3.23
294  smalla(ix^s)=entropycoef(il)
295  case('harten','powell', 'ratio')
296  ! Based on Harten & Hyman JCP 50, 235 and Zeeuw & Powell JCP 104,56
297  ! Initialize left and right eigenvalues by velocities
298  al(ix^s)= wl(ix^s,mom_c(idim))/wl(ix^s,rho_c_)
299  ar(ix^s)= wr(ix^s,mom_c(idim))/wr(ix^s,rho_c_)
300  ! Calculate the final "aL" and "aR"
301  select case(il)
302  case(fastrw_)
303  ! These quantities will be used for all the fast and slow waves
304  ! Calculate soundspeed**2 and cs**2+ca**2.
305  call twofl_get_csound2_c_from_conserved(wl,x,ixg^ll,ix^l,cs2l)
306  call twofl_get_csound2_c_from_conserved(wr,x,ixg^ll,ix^l,cs2r)
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_)
309  ! Save the discriminants into cs2L and cs2R
310  cs2l(ix^s)=&
311  dsqrt(cs2ca2l(ix^s)**2-4d0*cs2l(ix^s)*wl(ix^s,mag(idim))**2/wl(ix^s,rho_c_))
312  cs2r(ix^s)=&
313  dsqrt(cs2ca2r(ix^s)**2-4d0*cs2r(ix^s)*wr(ix^s,mag(idim))**2/wr(ix^s,rho_c_))
314 
315  ! The left and right eigenvalues for the fast wave going to right
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)))
318  case(fastlw_)
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)))
321  case(slowrw_)
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)))
324  case(slowlw_)
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_)
328  ! These propagate by the velocity
329  case(alfvrw_)
330  ! Store the Alfven speeds into cs2ca2L and cs2ca2R
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_))
333 
334  al(ix^s)=al(ix^s) + cs2ca2l(ix^s)
335  ar(ix^s)=ar(ix^s) + cs2ca2r(ix^s)
336  case(alfvlw_)
337  al(ix^s)=al(ix^s) - cs2ca2l(ix^s)
338  ar(ix^s)=ar(ix^s) - cs2ca2r(ix^s)
339  end select
340  end select
341 
342  call entropyfix(ix^l,il,al,ar,a,smalla)
343 
344  end subroutine geteigenjump2
345 
346  ! Multiply q by R(il,iw), where R is the right eigenvalue matrix at wroe
347  subroutine twofl_rtimes(q,w,ix^L,iw,il,idim,rq,workroe)
349 
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)
354 
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))
359 
360  end subroutine twofl_rtimes
361 
362  ! Multiply q by R(il,iw), where R is the right eigenvalue matrix at wroe
363  subroutine rtimes2(q,wroe,ix^L,iw,il,idim,rq, &
364  cfast,cslow,afast,aslow,csound2,dp,rhodv,bv,v2a2)
366 
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
372 
373  idir=idim+1-ndir*(idim/ndir)
374  jdir=idir+1-ndir*(idir/ndir)
375 
376  if(iw == rho_c_) then
377  select case(il)
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)
382  case(entrow_)
383  rq(ix^s)=q(ix^s)
384  case(diverw_,alfvrw_,alfvlw_)
385  rq(ix^s)=zero
386  end select
387  else if(iw == e_c_) then
388  if(il==fastrw_)then
389  ! Store 0.5*v**2+(2-gamma)/(gamma-1)*a**2
390  v2a2(ix^s)=half*sum(wroe(ix^s,mom_c(:))**2,dim=ndim+1)+ &
391  (two-twofl_gamma)/(twofl_gamma-one)*csound2(ix^s)
392  ! Store sgn(bx)*(betay*vy+betaz*vz) in bv
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
397  !Store betaz*vy-betay*vz in bv
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)))
400  endif
401 
402  select case(il)
403  case(fastrw_)
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)))))
406  case(fastlw_)
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)))))
409  case(slowrw_)
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)))))
412  case(slowlw_)
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)))))
415  case(entrow_)
416  rq(ix^s)= q(ix^s)*half*sum(wroe(ix^s,mom_c(:))**2,dim=ndim+1)
417  case(diverw_)
418  if(divbwave)then
419  rq(ix^s)= q(ix^s)*wroe(ix^s,mag(idim))
420  else
421  rq(ix^s)= zero
422  endif
423  case(alfvrw_)
424  rq(ix^s)=+q(ix^s)*bv(ix^s)
425  case(alfvlw_)
426  rq(ix^s)=-q(ix^s)*bv(ix^s)
427  end select
428  else if(any(mom_c(:)==iw)) then
429  if(iw==mom_c(idim))then
430  select case(il)
431  case(fastrw_)
432  rq(ix^s)=q(ix^s)*afast(ix^s)*(wroe(ix^s,iw)+cfast(ix^s))
433  case(fastlw_)
434  rq(ix^s)=q(ix^s)*afast(ix^s)*(wroe(ix^s,iw)-cfast(ix^s))
435  case(slowrw_)
436  rq(ix^s)=q(ix^s)*aslow(ix^s)*(wroe(ix^s,iw)+cslow(ix^s))
437  case(slowlw_)
438  rq(ix^s)=q(ix^s)*aslow(ix^s)*(wroe(ix^s,iw)-cslow(ix^s))
439  case(entrow_)
440  rq(ix^s)=q(ix^s)*wroe(ix^s,iw)
441  case(diverw_,alfvlw_,alfvrw_)
442  rq(ix^s)=zero
443  end select
444  else
445  select case(il)
446  case(fastrw_)
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))))
449  case(fastlw_)
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))))
452  case(slowrw_)
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))))
455  case(slowlw_)
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))))
458  case(entrow_)
459  rq(ix^s)=q(ix^s)*wroe(ix^s,iw)
460  case(diverw_)
461  rq(ix^s)=zero
462  case(alfvrw_)
463  if(iw==mom_c(idir))then
464  rq(ix^s)=+q(ix^s)*wroe(ix^s,mag(jdir))
465  else
466  rq(ix^s)=-q(ix^s)*wroe(ix^s,mag(idir))
467  endif
468  case(alfvlw_)
469  if(iw==mom_c(idir))then
470  rq(ix^s)=-q(ix^s)*wroe(ix^s,mag(jdir))
471  else
472  rq(ix^s)=+q(ix^s)*wroe(ix^s,mag(idir))
473  endif
474  end select
475  end if ! iw=m_idir,m_jdir
476  else if(any(mag(:)==iw)) then
477  if(iw==mag(idim))then
478  if(il==diverw_ .and. divbwave)then
479  rq(ix^s)=q(ix^s)
480  else
481  rq(ix^s)=zero
482  endif
483  else
484  select case(il)
485  case(fastrw_,fastlw_)
486  rq(ix^s)=+q(ix^s)*aslow(ix^s)*dsqrt(csound2(ix^s))*wroe(ix^s,iw)&
487  /wroe(ix^s,rho_c_)
488  case(slowrw_,slowlw_)
489  rq(ix^s)=-q(ix^s)*afast(ix^s)*dsqrt(csound2(ix^s))*wroe(ix^s,iw)&
490  /wroe(ix^s,rho_c_)
491  case(entrow_,diverw_)
492  rq(ix^s)=zero
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)))
497  else
498  rq(ix^s)=+q(ix^s)*wroe(ix^s,mag(idir))&
499  /sign(wroe(ix^s,rho_c_),wroe(ix^s,mag(idim)))
500  end if
501  end select
502  end if ! iw=b_idir,b_jdir
503  end if
504 
505  end subroutine rtimes2
506 
507 end module mod_twofl_roe
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.
Definition: mod_tvd.t:2
subroutine, public entropyfix(ixL, il, aL, aR, a, smalla)
Definition: mod_tvd.t:265
Magneto-hydrodynamics module.
Definition: mod_twofl_phys.t:2
double precision, public twofl_adiab
The adiabatic constant.
subroutine, public twofl_get_csound2_c_from_conserved(w, x, ixIL, ixOL, csound2)
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.
Definition: mod_twofl_roe.t:2
subroutine average2(wL, wR, x, ixL, idim, wroe, cfast, cslow, afast, aslow, csound2, dp, rhodv, tmp)
Definition: mod_twofl_roe.t:72
subroutine twofl_get_eigenjump(wL, wR, wroe, x, ixL, il, idim, smalla, a, jump, workroe)
subroutine twofl_rtimes(q, w, ixL, iw, il, idim, rq, workroe)
subroutine twofl_average(wL, wR, x, ixL, idim, wroe, workroe)
Definition: mod_twofl_roe.t:50
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 rtimes2(q, wroe, ixL, iw, il, idim, rq, cfast, cslow, afast, aslow, csound2, dp, rhodv, bv, v2a2)
subroutine, public twofl_roe_init()
Definition: mod_twofl_roe.t:20