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