MPI-AMRVAC  3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
mod_hd_roe.t
Go to the documentation of this file.
1 !> Module with Roe-type Riemann solver for hydrodynamics
2 module 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 
16 contains
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 
524 end 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::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::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::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:1022
mod_hd_phys::rho_
integer, public, protected rho_
Index of the density (in the w array)
Definition: mod_hd_phys.t:48
mod_hd_roe
Module with Roe-type Riemann solver for hydrodynamics.
Definition: mod_hd_roe.t:2
mod_hd_roe::hd_rtimes
subroutine hd_rtimes(q, wroe, ixL, iw, il, idim, rq, workroe)
Multiply q by R(il,iw), where R is the right eigenvalue matrix at wroe.
Definition: mod_hd_roe.t:224
mod_hd_roe::hd_average
subroutine hd_average(wL, wR, x, ixL, idim, wroe, workroe)
Calculate the Roe average of w, assignment of variables: rho -> rho, m -> v, e -> h.
Definition: mod_hd_roe.t:62
mod_hd_roe::hd_roe_init
subroutine, public hd_roe_init()
Definition: mod_hd_roe.t:19
mod_hd_roe::geteigenjump2
subroutine geteigenjump2(wL, wR, wroe, x, ixL, il, idim, smalla, a, jump, csound, dpperc2, dvperc)
Definition: mod_hd_roe.t:133
mod_hd_roe::rtimes2_iso
subroutine rtimes2_iso(q, wroe, ixL, iw, il, idim, rq, csound)
Definition: mod_hd_roe.t:485
mod_hd_roe::rtimes2
subroutine rtimes2(q, wroe, ixL, iw, il, idim, rq, csound)
Definition: mod_hd_roe.t:236
mod_hd_roe::geteigenjump2_iso
subroutine geteigenjump2_iso(wL, wR, wroe, x, ixL, il, idim, smalla, a, jump, csound)
Definition: mod_hd_roe.t:365
mod_hd_roe::hd_average_iso
subroutine hd_average_iso(wL, wR, x, ixL, idim, wroe, workroe)
Definition: mod_hd_roe.t:296
mod_hd_roe::average2_iso
subroutine average2_iso(wL, wR, x, ixL, idim, wroe, tmp)
Definition: mod_hd_roe.t:312
mod_hd_roe::hd_get_eigenjump
subroutine hd_get_eigenjump(wL, wR, wroe, x, ixL, il, idim, smalla, a, jump, workroe)
Calculate the il-th characteristic speed and the jump in the il-th characteristic variable in the idi...
Definition: mod_hd_roe.t:118
mod_hd_roe::hd_rtimes_iso
subroutine hd_rtimes_iso(q, wroe, ixL, iw, il, idim, rq, workroe)
Definition: mod_hd_roe.t:468
mod_hd_roe::hd_get_eigenjump_iso
subroutine hd_get_eigenjump_iso(wL, wR, wroe, x, ixL, il, idim, smalla, a, jump, workroe)
Definition: mod_hd_roe.t:346
mod_hd_roe::average2
subroutine average2(wL, wR, x, ixL, idim, wroe, tmp, tmp2)
Definition: mod_hd_roe.t:98
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.1