MPI-AMRVAC 3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
Loading...
Searching...
No Matches
mod_hd_roe.t
Go to the documentation of this file.
1!> Module with Roe-type Riemann solver for hydrodynamics
2module mod_hd_roe
3 use mod_hd_phys
4 use mod_physics_roe
5
6 implicit none
7 private
8
9 integer :: soundRW_ = -1
10 integer :: soundLW_ = -1
11 integer :: entropW_ = -1
12 integer :: shearW0_ = -1
13
14 public :: hd_roe_init
15
16contains
17
18 subroutine hd_roe_init()
19 use mod_global_parameters, only: entropycoef, nw
20
21 integer :: iw
22
23 if (hd_energy) then
24 ! Characteristic waves
25 soundrw_ = 1
26 soundlw_ = 2
27 entropw_ = 3
28 shearw0_ = 3
29 nworkroe = 3
30
31 phys_average => hd_average
32 phys_get_eigenjump => hd_get_eigenjump
33 phys_rtimes => hd_rtimes
34 else
35 ! Characteristic waves
36 soundrw_ = 1
37 soundlw_ = 2
38 shearw0_ = 2
39 nworkroe = 1
40
41 phys_average => hd_average_iso
42 phys_get_eigenjump => hd_get_eigenjump_iso
43 phys_rtimes => hd_rtimes_iso
44 end if
45
46 allocate(entropycoef(nw))
47
48 do iw = 1, nw
49 if (iw == soundrw_ .or. iw == soundlw_) then
50 ! TODO: Jannis: what's this?
51 entropycoef(iw) = 0.2d0
52 else
53 entropycoef(iw) = -1.0d0
54 end if
55 end do
56
57 end subroutine hd_roe_init
58
59 !> Calculate the Roe average of w, assignment of variables:
60 !> rho -> rho, m -> v, e -> h
61 subroutine hd_average(wL,wR,x,ix^L,idim,wroe,workroe)
62 use mod_global_parameters
63 integer, intent(in) :: ix^l, idim
64 double precision, intent(in) :: wl(ixg^t, nw), wr(ixg^t, nw)
65 double precision, intent(inout) :: wroe(ixg^t, nw)
66 double precision, intent(inout) :: workroe(ixg^t, nworkroe)
67 double precision, intent(in) :: x(ixg^t, 1:^nd)
68 integer :: idir
69
70 ! call average2(wL,wR,x,ix^L,idim,wroe,workroe(ixG^T,1),workroe(ixG^T,2))
71 workroe(ix^s, 1) = sqrt(wl(ix^s,rho_))
72 workroe(ix^s, 2) = sqrt(wr(ix^s,rho_))
73
74 ! The averaged density is sqrt(rhoL*rhoR)
75 wroe(ix^s,rho_) = workroe(ix^s, 1)*workroe(ix^s, 2)
76
77 ! Now the ratio sqrt(rhoL/rhoR) is put into workroe(ix^S, 1)
78 workroe(ix^s, 1) = workroe(ix^s, 1)/workroe(ix^s, 2)
79
80 ! Roe-average velocities
81 do idir = 1, ndir
82 wroe(ix^s,mom(idir)) = (wl(ix^s,mom(idir))/wl(ix^s,rho_) * workroe(ix^s, 1)+&
83 wr(ix^s,mom(idir))/wr(ix^s,rho_))/(one+workroe(ix^s, 1))
84 end do
85
86 ! Calculate enthalpyL, then enthalpyR, then Roe-average. Use tmp2 for pressure.
87 call hd_get_pthermal(wl,x,ixg^ll,ix^l, workroe(ixg^t, 2))
88
89 wroe(ix^s,e_) = (workroe(ix^s, 2)+wl(ix^s,e_))/wl(ix^s,rho_)
90
91 call hd_get_pthermal(wr,x,ixg^ll,ix^l, workroe(ixg^t, 2))
92
93 workroe(ix^s, 2) = (workroe(ix^s, 2)+wr(ix^s,e_))/wr(ix^s,rho_)
94 wroe(ix^s,e_) = (wroe(ix^s,e_)*workroe(ix^s, 1) + workroe(ix^s, 2))/(one+workroe(ix^s, 1))
95 end subroutine hd_average
96
97 subroutine average2(wL,wR,x,ix^L,idim,wroe,tmp,tmp2)
98
99 ! Calculate the Roe average of w, assignment of variables:
100 ! rho -> rho, m -> v, e -> h
101
102 use mod_global_parameters
103
104 integer :: ix^l,idim,idir
105 double precision, dimension(ixG^T,nw) :: wl,wr,wroe
106 double precision, dimension(ixG^T,ndim), intent(in) :: x
107 double precision, dimension(ixG^T) :: tmp,tmp2
108 !-----------------------------------------------------------------------------
109
110
111 end subroutine average2
112
113 !> Calculate the il-th characteristic speed and the jump in the il-th
114 !> characteristic variable in the idim direction within ixL.
115 !> The eigenvalues and the L=R**(-1) matrix is calculated from wroe.
116 !> jump(il)=Sum_il L(il,iw)*(wR(iw)-wL(iw))
117 subroutine hd_get_eigenjump(wL,wR,wroe,x,ix^L,il,idim,smalla,a,jump,workroe)
118 use mod_global_parameters
119
120 integer, intent(in) :: ix^l,il,idim
121 double precision, dimension(ixG^T,nw):: wl,wr,wroe
122 double precision, dimension(ixG^T,ndim), intent(in) :: x
123 double precision, dimension(ixG^T) :: smalla,a,jump
124 double precision, dimension(ixG^T,nworkroe) :: workroe
125
126 call geteigenjump2(wl,wr,wroe,x,ix^l,il,idim,smalla,a,jump, &
127 workroe(ixg^t,1),workroe(ixg^t,2),workroe(ixg^t,3))
128
129 end subroutine hd_get_eigenjump
130
131 subroutine geteigenjump2(wL,wR,wroe,x,ix^L,il,idim,smalla,a,jump,&
132 csound,dpperc2,dvperc)
133
134 ! Calculate the il-th characteristic speed and the jump in the il-th
135 ! characteristic variable in the idim direction within ixL.
136 ! The eigenvalues and the L=R**(-1) matrix is calculated from wroe.
137 ! jump(il)=Sum_il L(il,iw)*(wR(iw)-wL(iw))
138
139 use mod_global_parameters
140 use mod_tvd
141
142 integer :: ix^l,il,idim,idir
143 double precision, dimension(ixG^T,nw) :: wl,wr,wroe
144 double precision, dimension(ixG^T,ndim), intent(in) :: x
145 double precision, dimension(ixG^T) :: smalla,a,jump,tmp,tmp2
146 double precision, dimension(ixG^T) :: csound,dpperc2,dvperc
147 double precision :: kin_en(ixg^t)
148
149 if(il==1)then
150 !First calculate the square of the sound speed: c**2=(gamma-1)*(h-0.5*v**2)
151 kin_en(ix^s) = 0.5d0 * sum(wroe(ix^s, mom(:))**2, dim=^nd+1)
152 csound(ix^s)=(hd_gamma-one)*(wroe(ix^s,e_) - kin_en(ix^s))
153 ! Make sure that csound**2 is positive
154 csound(ix^s)=max(hd_gamma*smalldouble/wroe(ix^s,rho_),csound(ix^s))
155
156 ! Calculate (pR-pL)/c**2
157 ! To save memory we use tmp amnd tmp2 for pL and pR (hd_get_pthermal is OK)
158 call hd_get_pthermal(wl,x,ixg^ll,ix^l,tmp)
159 call hd_get_pthermal(wr,x,ixg^ll,ix^l,tmp2)
160 dpperc2(ix^s)=(tmp2(ix^s)-tmp(ix^s))/csound(ix^s)
161
162 !Now get the correct sound speed
163 csound(ix^s)=sqrt(csound(ix^s))
164
165 ! Calculate (vR_idim-vL_idim)/c
166 dvperc(ix^s)=(wr(ix^s,mom(idim))/wr(ix^s,rho_)-&
167 wl(ix^s,mom(idim))/wl(ix^s,rho_))/csound(ix^s)
168
169 endif
170
171 if (il == soundrw_) then
172 a(ix^s)=wroe(ix^s,mom(idim))+csound(ix^s)
173 jump(ix^s)=half*(dpperc2(ix^s)+wroe(ix^s,rho_)*dvperc(ix^s))
174 else if (il == soundlw_) then
175 a(ix^s)=wroe(ix^s,mom(idim))-csound(ix^s)
176 jump(ix^s)=half*(dpperc2(ix^s)-wroe(ix^s,rho_)*dvperc(ix^s))
177 else if (il == entropw_) then
178 a(ix^s)=wroe(ix^s,mom(idim))
179 jump(ix^s)=-dpperc2(ix^s)+wr(ix^s,rho_)-wl(ix^s,rho_)
180 else
181 !Determine the direction of the shear wave
182 idir=il-shearw0_; if(idir>=idim)idir=idir+1
183 a(ix^s)=wroe(ix^s,mom(idim))
184 jump(ix^s)=wroe(ix^s,rho_)*&
185 (wr(ix^s,mom(idir))/wr(ix^s,rho_)-wl(ix^s,mom(idir))/wl(ix^s,rho_))
186 end if
187
188 ! Calculate "smalla" or modify "a" based on the "typeentropy" switch
189 ! Put left and right eigenvalues, if needed, into tmp and tmp2
190 ! OK, since subroutines hd_get_pthermal and entropyfix do not use tmp and tmp2
191
192 select case(typeentropy(il))
193 case('yee')
194 ! Based on Yee JCP 68,151 eq 3.23
195 smalla(ix^s)=entropycoef(il)
196 case('harten','powell')
197 ! Based on Harten & Hyman JCP 50, 235 and Zeeuw & Powell JCP 104,56
198 if (il == soundrw_) then
199 call hd_get_pthermal(wl,x,ixg^ll,ix^l,tmp)
200 tmp(ix^s)=wl(ix^s,mom(idim))/wl(ix^s,rho_)&
201 + sqrt(hd_gamma*tmp(ix^s)/wl(ix^s,rho_))
202 call hd_get_pthermal(wr,x,ixg^ll,ix^l,tmp2)
203 tmp2(ix^s)=wr(ix^s,mom(idim))/wr(ix^s,rho_)&
204 + sqrt(hd_gamma*tmp2(ix^s)/wr(ix^s,rho_))
205 else if (il == soundlw_) then
206 call hd_get_pthermal(wl,x,ixg^ll,ix^l,tmp)
207 tmp(ix^s)=wl(ix^s,mom(idim))/wl(ix^s,rho_)&
208 - sqrt(hd_gamma*tmp(ix^s)/wl(ix^s,rho_))
209 call hd_get_pthermal(wr,x,ixg^ll,ix^l,tmp2)
210 tmp2(ix^s)=wr(ix^s,mom(idim))/wr(ix^s,rho_)&
211 - sqrt(hd_gamma*tmp2(ix^s)/wr(ix^s,rho_))
212 else
213 tmp(ix^s) =wl(ix^s,mom(idim))/wl(ix^s,rho_)
214 tmp2(ix^s)=wr(ix^s,mom(idim))/wr(ix^s,rho_)
215 end if
216 end select
217
218 call entropyfix(ix^l,il,tmp,tmp2,a,smalla)
219
220 end subroutine geteigenjump2
221
222 !> Multiply q by R(il,iw), where R is the right eigenvalue matrix at wroe
223 subroutine hd_rtimes(q,wroe,ix^L,iw,il,idim,rq,workroe)
224 use mod_global_parameters
225
226 integer, intent(in) :: ix^l,iw,il,idim
227 double precision, intent(in) :: wroe(ixg^t,nw)
228 double precision, intent(in) :: q(ixg^t)
229 double precision, intent(inout) :: rq(ixg^t)
230 double precision, intent(inout) :: workroe(ixg^t,nworkroe)
231 !-----------------------------------------------------------------------------
232 call rtimes2(q,wroe,ix^l,iw,il,idim,rq,workroe(ixg^t,1))
233 end subroutine hd_rtimes
234
235 subroutine rtimes2(q,wroe,ix^L,iw,il,idim,rq,csound)
236
237 ! Multiply q by R(il,iw), where R is the right eigenvalue matrix at wroe
238
239 use mod_global_parameters
240
241 integer :: ix^l,iw,il,idim,idir
242 double precision :: wroe(ixg^t,nw)
243 double precision, dimension(ixG^T) :: q,rq,csound
244 logical :: shearwave
245
246 shearwave=il>shearw0_
247 idir=idim
248 if(shearwave)then
249 ! Direction of shearwave increases with il plus idir==idim is jumped over
250 idir=il-shearw0_; if(idir>=idim)idir=idir+1
251 endif
252
253 if (iw == rho_) then
254 if(shearwave)then
255 rq(ix^s)=zero
256 else
257 rq(ix^s)=q(ix^s)
258 endif
259 else if (iw == e_) then
260 if (il == soundrw_) then
261 rq(ix^s)=q(ix^s)*(wroe(ix^s,e_)+wroe(ix^s,mom(idim))*csound(ix^s))
262 else if (il == soundlw_) then
263 rq(ix^s)=q(ix^s)*(wroe(ix^s,e_)-wroe(ix^s,mom(idim))*csound(ix^s))
264 else if (il == entropw_) then
265 rq(ix^s)=q(ix^s) * 0.5d0 * sum(wroe(ix^s, mom(:))**2, dim=^nd+1)
266 else
267 rq(ix^s)=q(ix^s)*wroe(ix^s,mom(idir))
268 end if
269 else
270 if(iw==mom(idim))then
271 if (il == soundrw_) then
272 rq(ix^s)=q(ix^s)*(wroe(ix^s,mom(idim))+csound(ix^s))
273 else if (il == soundlw_) then
274 rq(ix^s)=q(ix^s)*(wroe(ix^s,mom(idim))-csound(ix^s))
275 else if (il == entropw_) then
276 rq(ix^s)=q(ix^s)*wroe(ix^s,mom(idim))
277 else
278 rq(ix^s)=zero
279 end if
280 else
281 if(shearwave)then
282 if(iw==mom(idir))then
283 rq(ix^s)=q(ix^s)
284 else
285 rq(ix^s)=zero
286 endif
287 else
288 rq(ix^s)=q(ix^s)*wroe(ix^s,iw)
289 endif
290 endif
291 end if
292
293 end subroutine rtimes2
294
295 subroutine hd_average_iso(wL,wR,x,ix^L,idim,wroe,workroe)
296
297 ! Calculate the Roe average of w, assignment of variables:
298 ! rho -> rho, m -> v
299 use mod_global_parameters
300
301 integer, intent(in) :: ix^l, idim
302 double precision, intent(in) :: wl(ixg^t, nw), wr(ixg^t, nw)
303 double precision, intent(inout) :: wroe(ixg^t, nw)
304 double precision, intent(inout) :: workroe(ixg^t, nworkroe)
305 double precision, intent(in) :: x(ixg^t, 1:^nd)
306
307 call average2_iso(wl,wr,x,ix^l,idim,wroe,workroe(ixg^t,1))
308
309 end subroutine hd_average_iso
310
311 subroutine average2_iso(wL,wR,x,ix^L,idim,wroe,tmp)
312
313 ! Calculate the Roe average of w, assignment of variables:
314 ! rho -> rho, m -> v
315
316 use mod_global_parameters
317
318 integer:: ix^l,idim,idir
319 double precision, dimension(ixG^T,nw):: wl,wr,wroe
320 double precision, intent(in) :: x(ixg^t,1:ndim)
321 double precision, dimension(ixG^T):: tmp
322 !-----------------------------------------------------------------------------
323
324 select case (typeaverage)
325 case ('arithmetic')
326 ! This is the simple arithmetic average
327 wroe(ix^s,rho_)=half*(wl(ix^s,rho_)+wr(ix^s,rho_))
328 do idir = 1, ndir
329 wroe(ix^s, mom(idir)) = half * (wl(ix^s,mom(idir))/wl(ix^s,rho_) + &
330 wr(ix^s,mom(idir))/wr(ix^s,rho_))
331 end do
332 case ('roe','default')
333 ! Calculate the Roe-average
334 wroe(ix^s,rho_)=sqrt(wl(ix^s,rho_)*wr(ix^s,rho_))
335 ! Roe-average velocities
336 tmp(ix^s)=sqrt(wl(ix^s,rho_)/wr(ix^s,rho_))
337 do idir=1,ndir
338 wroe(ix^s,mom(idir))=(wl(ix^s,mom(idir))/wl(ix^s,rho_)*tmp(ix^s)+&
339 wr(ix^s,mom(idir))/wr(ix^s,rho_))/(one+tmp(ix^s))
340 end do
341 end select
342
343 end subroutine average2_iso
344
345 subroutine hd_get_eigenjump_iso(wL,wR,wroe,x,ix^L,il,idim,smalla,a,jump,workroe)
346
347 ! Calculate the il-th characteristic speed and the jump in the il-th
348 ! characteristic variable in the idim direction within ixL.
349 ! The eigenvalues and the L=R**(-1) matrix is calculated from wroe.
350 ! jump(il)=Sum_il L(il,iw)*(wR(iw)-wL(iw))
351
352 use mod_global_parameters
353
354 integer, intent(in) :: ix^l,il,idim
355 double precision, dimension(ixG^T,nw):: wl,wr,wroe
356 double precision, intent(in) :: x(ixg^t,1:ndim)
357 double precision, dimension(ixG^T) :: smalla,a,jump
358 double precision, dimension(ixG^T,nworkroe) :: workroe
359 !-----------------------------------------------------------------------------
360 call geteigenjump2_iso(wl,wr,wroe,x,ix^l,il,idim,smalla,a,jump,workroe(ixg^t,1))
361
362 end subroutine hd_get_eigenjump_iso
363
364 subroutine geteigenjump2_iso(wL,wR,wroe,x,ix^L,il,idim,smalla,a,jump,csound)
365
366 ! Calculate the il-th characteristic speed and the jump in the il-th
367 ! characteristic variable in the idim direction within ixL.
368 ! The eigenvalues and the L=R**(-1) matrix is calculated from wroe.
369 ! jump(il)=Sum_il L(il,iw)*(wR(iw)-wL(iw))
370
371 use mod_global_parameters
372 use mod_tvd
373
374 integer:: ix^l,il,idim,idir
375 double precision, dimension(ixG^T,nw):: wl,wr,wroe
376 double precision, intent(in) :: x(ixg^t,1:ndim)
377 double precision, dimension(ixG^T) :: smalla,a,jump,tmp,tmp2
378 double precision, dimension(ixG^T) :: csound
379 DOUBLE PRECISION,PARAMETER:: qsmall=1.d-6
380 !-----------------------------------------------------------------------------
381
382 select case (typeaverage)
383 case ('arithmetic')
384 call hd_get_pthermal(wroe,x,ixg^ll,ix^l,csound)
385 csound(ix^s) = sqrt(hd_gamma*csound(ix^s)/wroe(ix^s,rho_))
386 ! This is the original simple Roe-solver
387 if (il == soundrw_) then
388 a(ix^s)=wroe(ix^s,mom(idim))+csound(ix^s)
389 jump(ix^s)=half*((one-wroe(ix^s,mom(idim))/csound(ix^s))*&
390 (wr(ix^s,rho_)-wl(ix^s,rho_))&
391 +(wr(ix^s,mom(idim))-wl(ix^s,mom(idim)))/csound(ix^s))
392 else if (il == soundlw_) then
393 a(ix^s)=wroe(ix^s,mom(idim))-csound(ix^s)
394 jump(ix^s)=half*((one+wroe(ix^s,mom(idim))/csound(ix^s))*&
395 (wr(ix^s,rho_)-wl(ix^s,rho_))&
396 -(wr(ix^s,mom(idim))-wl(ix^s,mom(idim)))/csound(ix^s))
397 else
398 ! Determine direction of shear wave
399 idir=il-shearw0_; if(idir>=idim)idir=idir+1
400 a(ix^s)=wroe(ix^s,mom(idim))
401 jump(ix^s)=-wroe(ix^s,mom(idir))*(wr(ix^s,rho_)-wl(ix^s,rho_))&
402 +(wr(ix^s,mom(idir))-wl(ix^s,mom(idir)))
403 end if
404 case ('roe','default')
405 call hd_get_pthermal(wroe,x,ixg^ll,ix^l,csound)
406 call hd_get_pthermal(wl,x,ixg^ll,ix^l,tmp)
407 call hd_get_pthermal(wr,x,ixg^ll,ix^l,tmp2)
408 where(abs(wl(ix^s,rho_)-wr(ix^s,rho_))<=qsmall*(wl(ix^s,rho_)+wr(ix^s,rho_)))
409 csound(ix^s) = sqrt(hd_gamma*csound(ix^s)/wroe(ix^s,rho_))
410 elsewhere
411 csound(ix^s) = sqrt(hd_gamma*(tmp2(ix^s)-tmp(ix^s))/&
412 (wr(ix^s,rho_)-wl(ix^s,rho_)))
413 end where
414 ! This is the Roe solver by Glaister
415 ! based on P. Glaister JCP 93, 477-480 (1991)
416 if (il == soundrw_) then
417 a(ix^s)=wroe(ix^s,mom(idim))+csound(ix^s)
418 jump(ix^s)=half*((wr(ix^s,rho_)-wl(ix^s,rho_))+&
419 wroe(ix^s,rho_)/csound(ix^s)*(wr(ix^s,mom(idim))/wr(ix^s,rho_)-&
420 wl(ix^s,mom(idim))/wl(ix^s,rho_)))
421 else if (il == soundlw_) then
422 a(ix^s)=wroe(ix^s,mom(idim))-csound(ix^s)
423 jump(ix^s)=half*((wr(ix^s,rho_)-wl(ix^s,rho_))-&
424 wroe(ix^s,rho_)/csound(ix^s)*(wr(ix^s,mom(idim))/wr(ix^s,rho_)-&
425 wl(ix^s,mom(idim))/wl(ix^s,rho_)))
426 else
427 ! Determine direction of shear wave
428 idir=il-shearw0_; if(idir>=idim)idir=idir+1
429 a(ix^s)=wroe(ix^s,mom(idim))
430 jump(ix^s)=wroe(ix^s,rho_)*(wr(ix^s,mom(idir))/wr(ix^s,rho_)-&
431 wl(ix^s,mom(idir))/wl(ix^s,rho_))
432 end if
433 end select
434
435 ! Calculate "smalla" or modify "a" based on the "typeentropy" switch
436 ! Use tmp and tmp2 for the left and right eigenvalues if needed
437 select case(typeentropy(il))
438 case('yee')
439 ! Based on Yee JCP 68,151 eq 3.23
440 smalla(ix^s)=entropycoef(il)
441 case('harten','powell')
442 ! Based on Harten & Hyman JCP 50, 235 and Zeeuw & Powell JCP 104,56
443 if (il == soundrw_) then
444 call hd_get_pthermal(wl,x,ixg^ll,ix^l,tmp)
445 tmp(ix^s) =wl(ix^s,mom(idim))/wl(ix^s,rho_)&
446 + sqrt(hd_gamma*tmp(ix^s)/wl(ix^s,rho_))
447 call hd_get_pthermal(wr,x,ixg^ll,ix^l,tmp2)
448 tmp2(ix^s)=wr(ix^s,mom(idim))/wr(ix^s,rho_)&
449 + sqrt(hd_gamma*tmp2(ix^s)/wr(ix^s,rho_))
450 else if (il == soundlw_) then
451 call hd_get_pthermal(wl,x,ixg^ll,ix^l,tmp)
452 tmp(ix^s) =wl(ix^s,mom(idim))/wl(ix^s,rho_)&
453 - sqrt(hd_gamma*tmp(ix^s)/wl(ix^s,rho_))
454 call hd_get_pthermal(wr,x,ixg^ll,ix^l,tmp2)
455 tmp2(ix^s)=wr(ix^s,mom(idim))/wr(ix^s,rho_)&
456 - sqrt(hd_gamma*tmp2(ix^s)/wr(ix^s,rho_))
457 else
458 tmp(ix^s) =wl(ix^s,mom(idim))/wl(ix^s,rho_)
459 tmp2(ix^s)=wr(ix^s,mom(idim))/wr(ix^s,rho_)
460 end if
461 end select
462
463 call entropyfix(ix^l,il,tmp,tmp2,a,smalla)
464
465 end subroutine geteigenjump2_iso
466
467 subroutine hd_rtimes_iso(q,wroe,ix^L,iw,il,idim,rq,workroe)
468
469 ! Multiply q by R(il,iw), where R is the right eigenvalue matrix at wroe
470
471 use mod_global_parameters
472
473 integer, intent(in) :: ix^l,iw,il,idim
474 double precision, intent(in) :: wroe(ixg^t,nw)
475 double precision, intent(in) :: q(ixg^t)
476 double precision, intent(inout) :: rq(ixg^t)
477 double precision, intent(inout) :: workroe(ixg^t,nworkroe)
478 !-----------------------------------------------------------------------------
479
480 call rtimes2_iso(q,wroe,ix^l,iw,il,idim,rq,workroe(ixg^t,1))
481
482 end subroutine hd_rtimes_iso
483
484 subroutine rtimes2_iso(q,wroe,ix^L,iw,il,idim,rq,csound)
485
486 ! Multiply q by R(il,iw), where R is the right eigenvalue matrix at wroe
487
488 use mod_global_parameters
489
490 integer:: ix^l,iw,il,idim,idir
491 double precision:: wroe(ixg^t,nw)
492 double precision, dimension(ixG^T):: q,rq,csound
493
494 if(iw==rho_)then
495 if (il == soundrw_ .or. il == soundlw_) then
496 rq(ix^s)=q(ix^s)
497 else
498 rq(ix^s)=zero
499 end if
500 else if(iw==mom(idim))then
501 if (il == soundrw_) then
502 rq(ix^s)=q(ix^s)*(wroe(ix^s,mom(idim))+csound(ix^s))
503 else if (il == soundlw_) then
504 rq(ix^s)=q(ix^s)*(wroe(ix^s,mom(idim))-csound(ix^s))
505 else
506 rq(ix^s)=zero
507 end if
508 else
509 if (il == soundrw_ .or. il == soundlw_) then
510 rq(ix^s)=q(ix^s)*wroe(ix^s,iw)
511 else
512 !Determine direction of shear wave
513 idir=il-shearw0_; if(idir>=idim)idir=idir+1
514 if(iw==mom(idir)) then
515 rq(ix^s)=q(ix^s)
516 else
517 rq(ix^s)=zero
518 endif
519 end if
520 endif
521
522 end subroutine rtimes2_iso
523
524end module mod_hd_roe
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::d
double precision, dimension(:), allocatable, parameter d
Definition mod_global_parameters.t:23
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_global_parameters::typeaverage
character(len=std_len) typeaverage
Definition mod_global_parameters.t:746
mod_hd_phys
Hydrodynamics physics module.
Definition mod_hd_phys.t:2
mod_hd_phys::hd_energy
logical, public, protected hd_energy
Whether an energy equation is used.
Definition mod_hd_phys.t:12
mod_hd_phys::e_
integer, public, protected e_
Index of the energy density (-1 if not present)
Definition mod_hd_phys.t:57
mod_hd_phys::hd_gamma
double precision, public hd_gamma
The adiabatic index.
Definition mod_hd_phys.t:69
mod_hd_phys::mom
integer, dimension(:), allocatable, public, protected mom
Indices of the momentum density.
Definition mod_hd_phys.t:51
mod_hd_phys::rho_
integer, public, protected rho_
Index of the density (in the w array)
Definition mod_hd_phys.t:48
mod_hd_phys::hd_get_pthermal
subroutine, public hd_get_pthermal(w, x, ixil, ixol, pth)
Calculate thermal pressure=(gamma-1)*(e-0.5*m**2/rho) within ixO^L.
Definition mod_hd_phys.t:1011
mod_hd_roe
Module with Roe-type Riemann solver for hydrodynamics.
Definition mod_hd_roe.t:2
mod_hd_roe::hd_roe_init
subroutine, public hd_roe_init()
Definition mod_hd_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
Generated by doxygen 1.9.8