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