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