MPI-AMRVAC 3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
Loading...
Searching...
No Matches
mod_twofl_roe.t
Go to the documentation of this file.
1!> Subroutines for Roe-type Riemann solver for HD
2module mod_twofl_roe
3#include "amrvac.h"
4
5 use mod_twofl_phys
6 use mod_functions_bfield, only: mag
7 use mod_physics_roe
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
17contains
18
19 subroutine twofl_roe_init()
20 use mod_global_parameters
21 integer :: il
22
23 phys_average => twofl_average
24 phys_get_eigenjump => twofl_get_eigenjump
25 phys_rtimes => twofl_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 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)
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 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)
72 use mod_global_parameters
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)
179 use mod_global_parameters
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)
207 use mod_global_parameters
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)
348 use mod_global_parameters
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)
365 use mod_global_parameters
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
507end module mod_twofl_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_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_twofl_phys
Magneto-hydrodynamics module.
Definition mod_twofl_phys.t:2
mod_twofl_phys::twofl_get_csound2_c_from_conserved
subroutine, public twofl_get_csound2_c_from_conserved(w, x, ixil, ixol, csound2)
Definition mod_twofl_phys.t:3034
mod_twofl_phys::twofl_adiab
double precision, public twofl_adiab
The adiabatic constant.
Definition mod_twofl_phys.t:152
mod_twofl_phys::twofl_get_pthermal_c
subroutine, public twofl_get_pthermal_c(w, x, ixil, ixol, pth)
Definition mod_twofl_phys.t:3640
mod_twofl_phys::e_c_
integer, public e_c_
Index of the energy density (-1 if not present)
Definition mod_twofl_phys.t:97
mod_twofl_phys::mom_c
integer, dimension(:), allocatable, public mom_c
Indices of the momentum density.
Definition mod_twofl_phys.t:94
mod_twofl_phys::divbwave
logical, public divbwave
Add divB wave in Roe solver.
Definition mod_twofl_phys.t:188
mod_twofl_phys::rho_c_
integer, public rho_c_
Index of the density (in the w array)
Definition mod_twofl_phys.t:91
mod_twofl_phys::twofl_eq_energy
integer, public, protected twofl_eq_energy
Definition mod_twofl_phys.t:29
mod_twofl_phys::twofl_gamma
double precision, public twofl_gamma
The adiabatic index.
Definition mod_twofl_phys.t:149
mod_twofl_roe
Subroutines for Roe-type Riemann solver for HD.
Definition mod_twofl_roe.t:2
mod_twofl_roe::twofl_roe_init
subroutine, public twofl_roe_init()
Definition mod_twofl_roe.t:20
Generated by doxygen 1.9.8