DISPATCH
riemann_mod.f90
1 module riemann_mod
2  use amr_parameters
3  use hydro_parameters
4  use const
5  USE io_mod
6  USE trace_mod
7  USE mpi_mod
8  implicit none
9  private
10  type riemann_t
11  contains
12  procedure, nopass:: approx
13  procedure, nopass:: acoustic
14  procedure, nopass:: llf
15  procedure, nopass:: hll
16  procedure, nopass:: hllc
17  procedure, nopass:: hllc_a
18  procedure, nopass:: hllc_v
19  end type
20  type(riemann_t), public:: rieman
21 contains
22 
23 !===============================================================================
24 !===============================================================================
25 subroutine approx(qleft,qright,fgdnv,ngrid,nvar)
26  integer::ngrid,nvar
27  real(dp),dimension(1:ngrid,1:nvar)::qleft,qright,qgdnv,fgdnv
28 
29  ! local arrays
30  real(dp),dimension(1:ngrid)::rl ,ul ,pl ,cl
31  real(dp),dimension(1:ngrid)::rr ,ur ,pr ,cr
32  real(dp),dimension(1:ngrid)::ro ,uo ,po ,co
33  real(dp),dimension(1:ngrid)::rstar,ustar,pstar,cstar
34  real(dp),dimension(1:ngrid)::wl ,wr ,wo
35  real(dp),dimension(1:ngrid)::sgnm ,spin ,spout,ushock
36  real(dp),dimension(1:ngrid)::frac ,delp ,pold
37  integer ,dimension(1:ngrid)::ind ,ind2
38 
39  ! local variables
40  real(dp)::smallp, gamma6, ql, qr, usr, usl, wwl, wwr, smallpp, entho, etot
41  integer ::i, j, n, iter, n_new
42 
43  ! constants
44  smallp = smallc**2/gamma
45  smallpp = smallr*smallp
46  gamma6 = (gamma+one)/(two*gamma)
47  entho = one/(gamma-one)
48 
49  ! Pressure, density and velocity
50  do i=3,ngrid-1
51  rl(i)=max(qleft(i,1),smallr)
52  ul(i)= qleft(i,2)
53  pl(i)=max(qleft(i,3),rl(i)*smallp)
54  rr(i)=max(qright(i,1),smallr)
55  ur(i)= qright(i,2)
56  pr(i)=max(qright(i,3),rr(i)*smallp)
57  end do
58 
59  ! Lagrangian sound speed
60  do i=3,ngrid-1
61  cl(i) = gamma*pl(i)*rl(i)
62  cr(i) = gamma*pr(i)*rr(i)
63  end do
64 
65  ! First guess
66  do i=3,ngrid-1
67  wl(i)= sqrt(cl(i)); wr(i) = sqrt(cr(i))
68  pstar(i) = ((wr(i)*pl(i)+wl(i)*pr(i))+wl(i)*wr(i)*(ul(i)-ur(i)))/(wl(i)+wr(i))
69  pstar(i) = max(pstar(i),0.0_dp)
70  pold(i)= pstar(i)
71  end do
72  n = ngrid
73  do i = 1,n
74  ind(i)=i
75  end do
76 
77  ! Newton-Raphson iterations to find pstar at the required accuracy
78  ! for a two-shock Riemann problem
79  do iter = 1,niter_riemann
80  do i=1,n
81  wwl=sqrt(cl(ind(i))*(one+gamma6*(pold(i)-pl(ind(i)))/pl(ind(i))))
82  wwr=sqrt(cr(ind(i))*(one+gamma6*(pold(i)-pr(ind(i)))/pr(ind(i))))
83  ql=two*wwl**3/(wwl**2+cl(ind(i)))
84  qr=two*wwr**3/(wwr**2+cr(ind(i)))
85  usl=ul(ind(i))-(pold(i)-pl(ind(i)))/wwl
86  usr=ur(ind(i))+(pold(i)-pr(ind(i)))/wwr
87  delp(i)=max(qr*ql/(qr+ql)*(usl-usr),-pold(i))
88  end do
89  do i=1,n
90  pold(i)=pold(i)+delp(i)
91  end do
92  ! Convergence indicator
93  do i=1,n
94  uo(i)=abs(delp(i)/(pold(i)+smallpp))
95  end do
96  n_new=0
97  do i=1,n
98  if(uo(i)>1.d-06)then
99  n_new=n_new+1
100  ind2(n_new)=ind(i)
101  po(n_new)=pold(i)
102  end if
103  end do
104  j=n_new
105  do i=1,n
106  if(uo(i)<=1.d-06)then
107  n_new=n_new+1
108  ind2(n_new)=ind(i)
109  po(n_new)=pold(i)
110  end if
111  end do
112  ind(1:n)=ind2(1:n)
113  pold(1:n)=po(1:n)
114  n=j
115  end do
116 
117  ! Star region pressure
118  ! for a two-shock Riemann problem
119  do i=3,ngrid-1
120  pstar(ind(i))=pold(i)
121  end do
122  do i=3,ngrid-1
123  wl(i) = sqrt(cl(i)*(one+gamma6*(pstar(i)-pl(i))/pl(i)))
124  wr(i) = sqrt(cr(i)*(one+gamma6*(pstar(i)-pr(i))/pr(i)))
125  end do
126 
127  ! Star region velocity
128  ! for a two shock Riemann problem
129  do i=3,ngrid-1
130  ustar(i) = half*(ul(i) + (pl(i)-pstar(i))/wl(i) + &
131  & ur(i) - (pr(i)-pstar(i))/wr(i) )
132  end do
133 
134  ! Left going or right going contact wave
135  do i=3,ngrid-1
136  sgnm(i) = sign(one,ustar(i))
137  end do
138 
139  ! Left or right unperturbed state
140  do i=3,ngrid-1
141  if(sgnm(i)==one)then
142  ro(i) = rl(i)
143  uo(i) = ul(i)
144  po(i) = pl(i)
145  wo(i) = wl(i)
146  else
147  ro(i) = rr(i)
148  uo(i) = ur(i)
149  po(i) = pr(i)
150  wo(i) = wr(i)
151  end if
152  end do
153  do i=3,ngrid-1
154  co(i) = max(smallc,sqrt(abs(gamma*po(i)/ro(i))))
155  end do
156 
157  ! Star region density
158  do i=3,ngrid-1
159  if(pstar(i)>= po(i))then
160  ! Shock
161  rstar(i) = ro(i)/(one+ro(i)*(po(i)-pstar(i))/wo(i)**2)
162  else
163  ! Rarefaction
164  rstar(i) = ro(i)*(pstar(i)/po(i))**(one/gamma)
165  end if
166  end do
167 
168  do i=3,ngrid-1
169  ! Prevent vacuum formation in star region
170  rstar(i) = max(rstar(i),smallr)
171  ! Star region sound speed
172  cstar(i) = sqrt(abs(gamma*pstar(i)/rstar(i)))
173  cstar(i) = max(cstar(i),smallc)
174  ! Compute rarefaction head and tail speed
175  spout(i) = co(i)-sgnm(i)*uo(i)
176  spin(i) = cstar(i)-sgnm(i)*ustar(i)
177  ! Compute shock speed
178  ushock(i) = wo(i)/ro(i)-sgnm(i)*uo(i)
179  end do
180 
181  do i=3,ngrid-1
182  if(pstar(i)>=po(i))then
183  spout(i)=ushock(i)
184  spin(i)=spout(i)
185  end if
186  end do
187 
188  ! Sample the solution at x/t=0
189  do i=3,ngrid-1
190  if(spout(i)<=zero)then
191  qgdnv(i,1) = ro(i)
192  qgdnv(i,2) = uo(i)
193  qgdnv(i,3) = po(i)
194  else if(spin(i)>=zero)then
195  qgdnv(i,1) = rstar(i)
196  qgdnv(i,2) = ustar(i)
197  qgdnv(i,3) = pstar(i)
198  else
199  frac(i)=spout(i)/(spout(i)-spin(i))
200  qgdnv(i,2) = frac(i)*ustar(i) + (one - frac(i))*uo(i)
201  qgdnv(i,3) = frac(i)*pstar(i) + (one - frac(i))*po(i)
202  qgdnv(i,1) = ro(i)*(qgdnv(i,3)/po(i))**(one/gamma)
203  end if
204  end do
205 
206  ! Passive scalars
207  do n = 4,nvar
208  do i=3,ngrid-1
209  if(sgnm(i)==one)then
210  qgdnv(i,n) = qleft(i,n)
211  else
212  qgdnv(i,n) = qright(i,n)
213  end if
214  end do
215  end do
216 
217  ! Compute fluxes
218  do i = 3,ngrid-1
219  fgdnv(i,1) = qgdnv(i,1)*qgdnv(i,2) ! Mass density
220  fgdnv(i,2) = qgdnv(i,3)+qgdnv(i,1)*qgdnv(i,2)**2 ! Normal momentum
221  etot = qgdnv(i,3)*entho + half*qgdnv(i,1)*qgdnv(i,2)**2
222  etot = etot + half*qgdnv(i,1)*qgdnv(i,4)**2
223  etot = etot + half*qgdnv(i,1)*qgdnv(i,5)**2
224  fgdnv(i,3) = qgdnv(i,2)*(etot+qgdnv(i,3)) ! Total energy
225  end do
226  ! Other advected quantities
227  do n = 4, nvar
228  do i = 3,ngrid-1
229  fgdnv(i,n) = fgdnv(i,1)*qgdnv(i,n)
230  end do
231  end do
232 end subroutine approx
233 
234 !===============================================================================
235 !===============================================================================
236 subroutine acoustic(qleft,qright,fgdnv,ngrid,nvar)
237  integer::ngrid,nvar
238  real(dp),dimension(1:ngrid,1:nvar)::qleft,qright,qgdnv,fgdnv
239 
240  ! local variables
241  integer::i,n
242  real(dp)::smallp, entho, etot
243 
244  ! local arrays
245  real(dp),dimension(1:ngrid)::rl ,ul ,pl ,cl
246  real(dp),dimension(1:ngrid)::rr ,ur ,pr ,cr
247  real(dp),dimension(1:ngrid)::ro ,uo ,po ,co
248  real(dp),dimension(1:ngrid)::rstar,ustar,pstar,cstar
249  real(dp),dimension(1:ngrid)::wl ,wr ,wo
250  real(dp),dimension(1:ngrid)::sgnm ,spin ,spout,ushock
251  real(dp),dimension(1:ngrid)::frac
252 
253  ! constants
254  smallp = smallc**2/gamma
255  entho = one/(gamma-one)
256 
257  ! Initial states pressure, density and velocity
258  do i=3,ngrid-1
259  rl(i)=max(qleft(i,1),smallr)
260  ul(i)= qleft(i,2)
261  pl(i)=max(qleft(i,3),rl(i)*smallp)
262  rr(i)=max(qright(i,1),smallr)
263  ur(i)= qright(i,2)
264  pr(i)=max(qright(i,3),rr(i)*smallp)
265  end do
266 
267  ! Acoustic star state
268  do i=3,ngrid-1
269  cl(i) = sqrt(gamma*pl(i)/rl(i))
270  cr(i) = sqrt(gamma*pr(i)/rr(i))
271  wl(i) = cl(i)*rl(i)
272  wr(i) = cr(i)*rr(i)
273  pstar(i) = ((wr(i)*pl(i)+wl(i)*pr(i))+wl(i)*wr(i)*(ul(i)-ur(i))) &
274  & / (wl(i)+wr(i))
275  ustar(i) = ((wr(i)*ur(i)+wl(i)*ul(i))+(pl(i)-pr(i))) &
276  & / (wl(i)+wr(i))
277 !!$ pstar(i) = MAX(pstar(i),zero)
278  end do
279 
280  ! Left going or right going contact wave
281  do i=3,ngrid-1
282  sgnm(i) = sign(one,ustar(i))
283  end do
284 
285  ! Left or right unperturbed state
286  do i=3,ngrid-1
287  if(sgnm(i)==one)then
288  ro(i) = rl(i)
289  uo(i) = ul(i)
290  po(i) = pl(i)
291  wo(i) = wl(i)
292  co(i) = cl(i)
293  else
294  ro(i) = rr(i)
295  uo(i) = ur(i)
296  po(i) = pr(i)
297  wo(i) = wr(i)
298  co(i) = cr(i)
299  end if
300  end do
301 
302  ! Star region density and sound speed
303  do i=3,ngrid-1
304  rstar(i) = ro(i)+(pstar(i)-po(i))/co(i)**2
305  rstar(i) = max(rstar(i),smallr)
306  cstar(i) = sqrt(abs(gamma*pstar(i)/rstar(i)))
307  cstar(i) = max(cstar(i),smallc)
308  end do
309 
310  ! Head and tail speed of rarefaction
311  do i=3,ngrid-1
312  spout(i) = co(i)-sgnm(i)*uo(i)
313  spin(i) = cstar(i)-sgnm(i)*ustar(i)
314  end do
315 
316  ! Shock speed
317  do i=3,ngrid-1
318  ushock(i) = half*(spin(i)+spout(i))
319  ushock(i) = max(ushock(i),-sgnm(i)*ustar(i))
320  end do
321  do i=3,ngrid-1
322  if(pstar(i)>=po(i))then
323  spout(i)=ushock(i)
324  spin(i)=spout(i)
325  end if
326  end do
327 
328  ! Sample the solution at x/t=0
329  do i=3,ngrid-1
330  if(spout(i)<zero)then ! Initial state
331  qgdnv(i,1) = ro(i)
332  qgdnv(i,2) = uo(i)
333  qgdnv(i,3) = po(i)
334  else if(spin(i)>=zero)then ! Star region
335  qgdnv(i,1) = rstar(i)
336  qgdnv(i,2) = ustar(i)
337  qgdnv(i,3) = pstar(i)
338  else ! Rarefaction
339  frac(i) = spout(i)/(spout(i)-spin(i))
340  qgdnv(i,1) = frac(i)*rstar(i) + (one - frac(i))*ro(i)
341  qgdnv(i,2) = frac(i)*ustar(i) + (one - frac(i))*uo(i)
342  qgdnv(i,3) = frac(i)*pstar(i) + (one - frac(i))*po(i)
343  end if
344  end do
345 
346  ! Passive scalars
347  do n = 4,nvar
348  do i=3,ngrid-1
349  if(sgnm(i)==one)then
350  qgdnv(i,n) = qleft(i,n)
351  else
352  qgdnv(i,n) = qright(i,n)
353  end if
354  end do
355  end do
356 
357  ! Compute fluxes
358  do i = 3,ngrid-1
359  fgdnv(i,1) = qgdnv(i,1)*qgdnv(i,2) ! Mass density
360  fgdnv(i,2) = qgdnv(i,3)+qgdnv(i,1)*qgdnv(i,2)**2 ! Normal momentum
361  etot = qgdnv(i,3)*entho + half*qgdnv(i,1)*qgdnv(i,2)**2
362  etot = etot + half*qgdnv(i,1)*qgdnv(i,4)**2
363  etot = etot + half*qgdnv(i,1)*qgdnv(i,5)**2
364  fgdnv(i,3) = qgdnv(i,2)*(etot+qgdnv(i,3)) ! Total energy
365  end do
366  ! Other advected quantities
367  do n = 4, nvar
368  do i = 3,ngrid-1
369  fgdnv(i,n) = fgdnv(i,1)*qgdnv(i,n)
370  end do
371  end do
372 end subroutine acoustic
373 
374 !===============================================================================
375 !===============================================================================
376 subroutine llf(qleft,qright,fgdnv,ngrid,nvar)
377  integer::ngrid,nvar
378  real(dp),dimension(1:ngrid,1:nvar)::qleft,qright,fgdnv
379 
380  ! local arrays
381  real(dp),dimension(1:ngrid,1:nvar)::fleft,fright,uleft,uright
382  real(dp),dimension(1:ngrid)::cmax
383 
384  ! local variables
385  integer::i,n
386  real(dp)::smallp, entho
387  real(dp)::rl ,ul ,pl ,cl
388  real(dp)::rr ,ur ,pr ,cr
389 
390  ! constants
391  smallp = smallc**2/gamma
392  entho = one/(gamma-one)
393 
394  ! Maximum wave speed
395  do i=3,ngrid-1
396  rl=max(qleft(i,1),smallr)
397  ul= qleft(i,2)
398  pl=max(qleft(i,3),rl*smallp)
399  rr=max(qright(i,1),smallr)
400  ur= qright(i,2)
401  pr=max(qright(i,3),rr*smallp)
402  cl= sqrt(gamma*pl/rl)
403  cr= sqrt(gamma*pr/rr)
404  cmax(i)=max(abs(ul)+cl,abs(ur)+cr)
405  end do
406 
407  ! Compute conservative variables
408  do i = 3,ngrid-1
409  ! Mass density
410  uleft(i,1) = qleft(i,1)
411  uright(i,1) = qright(i,1)
412  ! Normal momentum
413  uleft(i,2) = qleft(i,1)*qleft(i,2)
414  uright(i,2) = qright(i,1)*qright(i,2)
415  ! Total energy
416  uleft(i,3) = qleft(i,3)*entho + half*qleft(i,1)*qleft(i,2)**2
417  uright(i,3) = qright(i,3)*entho + half*qright(i,1)*qright(i,2)**2
418  uleft(i,3) = uleft(i,3) + half*qleft(i,1)*qleft(i,4)**2
419  uright(i,3) = uright(i,3) + half*qright(i,1)*qright(i,4)**2
420  uleft(i,3) = uleft(i,3) + half*qleft(i,1)*qleft(i,5)**2
421  uright(i,3) = uright(i,3) + half*qright(i,1)*qright(i,5)**2
422  end do
423  ! Other advected quantities
424  do n = 4, nvar
425  do i = 3,ngrid-1
426  uleft(i,n) = qleft(i,1)*qleft(i,n)
427  uright(i,n) = qright(i,1)*qright(i,n)
428  end do
429  end do
430 
431  ! Compute left and right fluxes
432  do i = 3,ngrid-1
433  ! Mass density
434  fleft(i,1) = uleft(i,2)
435  fright(i,1) = uright(i,2)
436  ! Normal momentum
437  fleft(i,2) = qleft(i,3)+uleft(i,2)*qleft(i,2)
438  fright(i,2) = qright(i,3)+uright(i,2)*qright(i,2)
439  ! Total energy
440  fleft(i,3) = qleft(i,2)*(uleft(i,3)+qleft(i,3))
441  fright(i,3) = qright(i,2)*(uright(i,3)+qright(i,3))
442  end do
443  ! Other advected quantities
444  do n = 4, nvar
445  do i = 3,ngrid-1
446  fleft(i,n) = fleft(i,1)*qleft(i,n)
447  fright(i,n) = fright(i,1)*qright(i,n)
448  end do
449  end do
450 
451  ! Compute Lax-Friedrich fluxes
452  do n = 1, nvar
453  do i = 3,ngrid-1
454  fgdnv(i,n) = half*(fleft(i,n)+fright(i,n)-cmax(i)*(uright(i,n)-uleft(i,n)))
455  end do
456  end do
457 end subroutine llf
458 
459 !===============================================================================
460 !===============================================================================
461 subroutine hllc (qleft,qright,fgdnv,ngrid,nvar,isothermal,error,detailed_timer)
462  integer::ngrid,nvar
463  real(dp),dimension(1:ngrid,1:nvar)::qleft,qright,fgdnv
464  REAL(dp)::sl,sr
465  REAL(dp)::entho
466  REAL(dp)::rl,pl,ul,ecinl,etotl,ptotl
467  REAL(dp)::rr,pr,ur,ecinr,etotr,ptotr
468  REAL(dp)::cfastl,rcl,rstarl
469  REAL(dp)::cfastr,rcr,rstarr
470  REAL(dp)::etotstarl,etotstarr
471  REAL(dp)::ustar,ptotstar
472  REAL(dp)::ro,uo,ptoto,etoto
473  REAL(dp)::smallp
474  INTEGER ::ivar, i
475  LOGICAL :: isothermal,error,detailed_timer
476  INTEGER, SAVE:: itimer=0
477  !
478  if (detailed_timer) call trace_begin ('riemann::hllc', itimer=itimer)
479  ! constants
480  smallp = smallc**2/gamma
481  if (isothermal .or. gamma==1d0) then
482  entho=one
483  else
484  entho=one/(gamma-one)
485  end if
486  error = .false.
487  !
488  if (io%verbose>1) then
489  do i=3,ngrid-1
490  !
491  ! Left variables
492  rl=max(qleft(i,1),smallr)
493  pl=max(qleft(i,3),rl*smallp)
494  ul= qleft(i,2)
495  !
496  ecinl = half*rl*(ul*ul+qleft(i,4)**2+qleft(i,5)**2)
497  etotl = pl*entho+ecinl
498  ptotl = pl
499  !
500  ! Right variables
501  rr=max(qright(i,1),smallr)
502  pr=max(qright(i,3),rr*smallp)
503  ur= qright(i,2)
504  !
505  ecinr = half*rr*(ur*ur+qright(i,4)**2+qright(i,5)**2)
506  etotr = pr*entho+ecinr
507  ptotr = pr
508  !
509  ! Find the largest eigenvalues in the normal direction to the interface
510  cfastl=sqrt(max(gamma*pl/rl,smallc**2))
511  cfastr=sqrt(max(gamma*pr/rr,smallc**2))
512  !
513  ! Compute HLL wave speed
514  sl=min(ul,ur)-max(cfastl,cfastr)
515  sr=max(ul,ur)+max(cfastl,cfastr)
516  !
517  ! Compute lagrangian sound speed
518  rcl=rl*(ul-sl)
519  rcr=rr*(sr-ur)
520  if ((rcl+rcl)==0.0) then
521  print 1,'HLLC: rl,Pl,ul,ecinl,etotl,Ptotl,cfastl,SL =', &
522  i,rl,pl,ul,ecinl,etotl,ptotl,cfastl,sl,sl-ul
523  print 1,'HLLC: rr,Pr,ur,ecinr,etotr,Ptotr,cfastr,SR =', &
524  i,rr,pr,ur,ecinr,etotr,ptotr,cfastr,sr,sr-ur
525  1 format(a,i4,1p,10e12.3)
526  error = .true.
527  return
528  end if
529  !
530  ! Compute acoustic star state
531  ustar =(rcr*ur +rcl*ul + (ptotl-ptotr))/(rcr+rcl)
532  ptotstar=(rcr*ptotl+rcl*ptotr+rcl*rcr*(ul-ur))/(rcr+rcl)
533  if ((sl-ustar)==0.0 .or. (sr-ustar)==0.0) then
534  print 1,'HLLC: rl,Pl,ul,ecinl,etotl,Ptotl,cfastl,SL =', &
535  i,rl,pl,ul,ecinl,etotl,ptotl,cfastl,sl,sl-ul
536  print 1,'HLLC: rr,Pr,ur,ecinr,etotr,Ptotr,cfastr,SR =', &
537  i,rr,pr,ur,ecinr,etotr,ptotr,cfastr,sr,sr-ur
538  print 1, 'HLLC: ustar,Ptotstar', i,ustar,ptotstar
539  call mpi%abort ('HLLC')
540  error = .true.
541  return
542  end if
543  !
544  ! Left star region variables
545  rstarl=rl*(sl-ul)/(sl-ustar)
546  etotstarl=((sl-ul)*etotl-ptotl*ul+ptotstar*ustar)/(sl-ustar)
547  !
548  ! Right star region variables
549  rstarr=rr*(sr-ur)/(sr-ustar)
550  etotstarr=((sr-ur)*etotr-ptotr*ur+ptotstar*ustar)/(sr-ustar)
551  !
552  ! Sample the solution at x/t=0
553  if(sl>0d0)then
554  ro=rl
555  uo=ul
556  ptoto=ptotl
557  etoto=etotl
558  else if(ustar>0d0)then
559  ro=rstarl
560  uo=ustar
561  ptoto=ptotstar
562  etoto=etotstarl
563  else if (sr>0d0)then
564  ro=rstarr
565  uo=ustar
566  ptoto=ptotstar
567  etoto=etotstarr
568  else
569  ro=rr
570  uo=ur
571  ptoto=ptotr
572  etoto=etotr
573  end if
574  !
575  ! Compute the Godunov flux
576  fgdnv(i,1) = ro*uo
577  fgdnv(i,2) = ro*uo*uo+ptoto
578  fgdnv(i,3) = (etoto+ptoto)*uo
579  end do
580  else
581  do i=3,ngrid-1
582  !
583  ! Left variables
584  rl=max(qleft(i,1),smallr)
585  pl=max(qleft(i,3),rl*smallp)
586  ul= qleft(i,2)
587  !
588  ecinl = half*rl*(ul*ul+qleft(i,4)**2+qleft(i,5)**2)
589  etotl = pl*entho+ecinl
590  ptotl = pl
591  !
592  ! Right variables
593  rr=max(qright(i,1),smallr)
594  pr=max(qright(i,3),rr*smallp)
595  ur= qright(i,2)
596  !
597  ecinr = half*rr*(ur*ur+qright(i,4)**2+qright(i,5)**2)
598  etotr = pr*entho+ecinr
599  ptotr = pr
600  !
601  ! Find the largest eigenvalues in the normal direction to the interface
602  cfastl=sqrt(max(gamma*pl/rl,smallc**2))
603  cfastr=sqrt(max(gamma*pr/rr,smallc**2))
604  !
605  ! Compute HLL wave speed
606  sl=min(ul,ur)-max(cfastl,cfastr)
607  sr=max(ul,ur)+max(cfastl,cfastr)
608  !
609  ! Compute lagrangian sound speed
610  rcl=rl*(ul-sl)
611  rcr=rr*(sr-ur)
612  !
613  ! Compute acoustic star state
614  ustar =(rcr*ur +rcl*ul + (ptotl-ptotr))/(rcr+rcl)
615  ptotstar=(rcr*ptotl+rcl*ptotr+rcl*rcr*(ul-ur))/(rcr+rcl)
616  !
617  ! Left star region variables
618  rstarl=rl*(sl-ul)/(sl-ustar)
619  etotstarl=((sl-ul)*etotl-ptotl*ul+ptotstar*ustar)/(sl-ustar)
620  !
621  ! Right star region variables
622  rstarr=rr*(sr-ur)/(sr-ustar)
623  etotstarr=((sr-ur)*etotr-ptotr*ur+ptotstar*ustar)/(sr-ustar)
624  !
625  ! Sample the solution at x/t=0
626 !define MERGE
627 #ifdef MERGE
628  ro=rr
629  uo=ur
630  ptoto=ptotr
631  etoto=etotr
632  ro =merge(rstarr ,ro ,sr>0d0)
633  uo =merge(ustar ,uo ,sr>0d0)
634  ptoto=merge(ptotstar ,ptoto ,sr>0d0)
635  etoto=merge(etotstarr,etoto ,sr>0d0)
636  ro =merge(rstarr ,ro ,ustar>0d0)
637  uo =merge(ustar ,uo ,ustar>0d0)
638  ptoto=merge(ptotstar ,ptoto ,ustar>0d0)
639  etoto=merge(etotstarl,etoto ,ustar>0d0)
640  ro =merge(rl ,ro ,sl>0d0)
641  uo =merge(ul ,uo ,sl>0d0)
642  ptoto=merge(ptotl ,ptoto ,sl>0d0)
643  etoto=merge(etotl ,etoto ,sl>0d0)
644 #else
645  if(sl>0d0)then
646  ro=rl
647  uo=ul
648  ptoto=ptotl
649  etoto=etotl
650  else if(ustar>0d0)then
651  ro=rstarl
652  uo=ustar
653  ptoto=ptotstar
654  etoto=etotstarl
655  else if (sr>0d0)then
656  ro=rstarr
657  uo=ustar
658  ptoto=ptotstar
659  etoto=etotstarr
660  else
661  ro=rr
662  uo=ur
663  ptoto=ptotr
664  etoto=etotr
665  end if
666 #endif
667  !
668  ! Compute the Godunov flux
669  fgdnv(i,1) = ro*uo
670  fgdnv(i,2) = ro*uo*uo+ptoto
671  fgdnv(i,3) = (etoto+ptoto)*uo
672  end do
673  end if
674  !
675  do ivar=4,nvar
676  do i=3,ngrid-1
677  fgdnv(i,ivar) = fgdnv(i,1)*merge(qleft(i,ivar),qright(i,ivar),fgdnv(i,1)>0)
678  end do
679  end do
680  if (detailed_timer) call trace_end (itimer)
681 end subroutine hllc
682 
683 !===============================================================================
684 !===============================================================================
685 subroutine hllc_a (qleft,qright,fgdnv,ngrid,nvar)
686  integer::ngrid,nvar
687  real(dp),dimension(1:ngrid,1:nvar)::qleft,qright,fgdnv
688  REAL(dp)::entho
689  REAL(dp)::pl,ecinl
690  REAL(dp)::pr,ecinr
691  REAL(dp)::cfastl,rcl
692  REAL(dp)::cfastr,rcr
693  REAL(dp)::smallp
694  REAL(dp),dimension(ngrid)::sl,sr
695  REAL(dp),dimension(ngrid)::ro,uo,ptoto,etoto
696  REAL(dp),dimension(ngrid)::rl,ul,ptotl,etotl
697  REAL(dp),dimension(ngrid)::rr,ur,ptotr,etotr
698  REAL(dp),dimension(ngrid)::rstarl,etotstarl
699  REAL(dp),dimension(ngrid)::rstarr,etotstarr
700  REAL(dp),dimension(ngrid)::ustar,ptotstar
701  INTEGER ::ivar, i
702  !
703  ! constants
704  smallp = smallc**2/gamma
705  entho = one/(gamma-one)
706  !
707  do i=3,ngrid-1
708  !
709  ! Left variables
710  rl(i)=max(qleft(i,1),smallr)
711  pl=max(qleft(i,3),rl(i)*smallp)
712  ul(i)= qleft(i,2)
713  !
714  ecinl = half*rl(i)*(ul(i)*ul(i)+qleft(i,4)**2+qleft(i,5)**2)
715  etotl(i) = pl*entho+ecinl
716  ptotl(i) = pl
717  !
718  ! Right variables
719  rr(i)=max(qright(i,1),smallr)
720  pr=max(qright(i,3),rr(i)*smallp)
721  ur(i)= qright(i,2)
722  !
723  ecinr = half*rr(i)*(ur(i)*ur(i)+qright(i,4)**2+qright(i,5)**2)
724  etotr(i) = pr*entho+ecinr
725  ptotr(i) = pr
726  !
727  ! Find the largest eigenvalues in the normal direction to the interface
728  cfastl=sqrt(max(gamma*pl/rl(i),smallc**2))
729  cfastr=sqrt(max(gamma*pr/rr(i),smallc**2))
730  !
731  ! Compute HLL wave speed
732  sl(i)=min(ul(i),ur(i))-max(cfastl,cfastr)
733  sr(i)=max(ul(i),ur(i))+max(cfastl,cfastr)
734  !
735  ! Compute lagrangian sound speed
736  rcl=rl(i)*(ul(i)-sl(i))
737  rcr=rr(i)*(sr(i)-ur(i))
738  !
739  ! Compute acoustic star state
740  ustar(i) =(rcr*ur(i) +rcl*ul(i) + (ptotl(i)-ptotr(i)))/(rcr+rcl)
741  ptotstar(i)=(rcr*ptotl(i)+rcl*ptotr(i)+rcl*rcr*(ul(i)-ur(i)))/(rcr+rcl)
742  end do
743  do i=3,ngrid-1
744  !
745  ! Left star region variables
746  rstarl(i)=rl(i)*(sl(i)-ul(i))/(sl(i)-ustar(i))
747  etotstarl(i)=((sl(i)-ul(i))*etotl(i)-ptotl(i)*ul(i)+ptotstar(i)*ustar(i))/(sl(i)-ustar(i))
748  !
749  ! Right star region variables
750  rstarr(i)=rr(i)*(sr(i)-ur(i))/(sr(i)-ustar(i))
751  etotstarr(i)=((sr(i)-ur(i))*etotr(i)-ptotr(i)*ur(i)+ptotstar(i)*ustar(i))/(sr(i)-ustar(i))
752  end do
753  !
754  ! Sample the solution at x/t=0
755  do i=3,ngrid-1
756  if(sl(i)>0d0)then
757  ro(i)=rl(i)
758  uo(i)=ul(i)
759  ptoto(i)=ptotl(i)
760  etoto(i)=etotl(i)
761  else if(ustar(i)>0d0)then
762  ro(i)=rstarl(i)
763  uo(i)=ustar(i)
764  ptoto(i)=ptotstar(i)
765  etoto(i)=etotstarl(i)
766  else if (sr(i)>0d0)then
767  ro(i)=rstarr(i)
768  uo(i)=ustar(i)
769  ptoto(i)=ptotstar(i)
770  etoto(i)=etotstarr(i)
771  else
772  ro(i)=rr(i)
773  uo(i)=ur(i)
774  ptoto(i)=ptotr(i)
775  etoto(i)=etotr(i)
776  end if
777  end do
778  !
779  ! Compute the Godunov flux
780  do i=3,ngrid-1
781  fgdnv(i,1) = ro(i)*uo(i)
782  fgdnv(i,2) = ro(i)*uo(i)*uo(i)+ptoto(i)
783  fgdnv(i,3) = (etoto(i)+ptoto(i))*uo(i)
784  end do
785  !
786  do ivar=4,nvar
787  do i=3,ngrid-1
788  fgdnv(i,ivar) = merge(fgdnv(i,1)*qleft(i,ivar), &
789  fgdnv(i,1)*qright(i,ivar), fgdnv(i,1)>0)
790  end do
791  end do
792 end subroutine hllc_a
793 
794 !===============================================================================
795 !===============================================================================
796 subroutine hllc_v (qleft,qright,fgdnv,ngrid,nvar)
797  integer::ngrid,nvar
798  real(dp),dimension(1:ngrid,1:nvar)::qleft,qright,fgdnv
799  REAL(dp),dimension(ngrid)::sl,sr
800  REAL(dp),dimension(ngrid)::entho
801  REAL(dp),dimension(ngrid)::rl,pl,ul,ecinl,etotl,ptotl
802  REAL(dp),dimension(ngrid)::rr,pr,ur,ecinr,etotr,ptotr
803  REAL(dp),dimension(ngrid)::cfastl,rcl,rstarl
804  REAL(dp),dimension(ngrid)::cfastr,rcr,rstarr
805  REAL(dp),dimension(ngrid)::etotstarl,etotstarr
806  REAL(dp),dimension(ngrid)::ustar,ptotstar
807  REAL(dp),dimension(ngrid)::ro,uo,ptoto,etoto
808  REAL(dp),dimension(ngrid)::smallp
809  INTEGER ::ivar, i
810  !
811  ! constants
812  smallp = smallc**2/gamma
813  entho = one/(gamma-one)
814  !
815  !
816  ! Left variables
817  rl=max(qleft(:,1),smallr)
818  pl=max(qleft(:,3),rl*smallp)
819  ul= qleft(:,2)
820  !
821  ecinl = half*rl*(ul*ul+qleft(:,4)**2+qleft(:,5)**2)
822  etotl = pl*entho+ecinl
823  ptotl = pl
824  !
825  ! Right variables
826  rr=max(qright(:,1),smallr)
827  pr=max(qright(:,3),rr*smallp)
828  ur= qright(:,2)
829  !
830  ecinr = half*rr*(ur*ur+qright(:,4)**2+qright(:,5)**2)
831  etotr = pr*entho+ecinr
832  ptotr = pr
833  !
834  ! Find the largest eigenvalues in the normal direction to the interface
835  cfastl=sqrt(max(gamma*pl/rl,smallc**2))
836  cfastr=sqrt(max(gamma*pr/rr,smallc**2))
837  !
838  ! Compute HLL wave speed
839  sl=min(ul,ur)-max(cfastl,cfastr)
840  sr=max(ul,ur)+max(cfastl,cfastr)
841  !
842  ! Compute lagrangian sound speed
843  rcl=rl*(ul-sl)
844  rcr=rr*(sr-ur)
845  !
846  ! Compute acoustic star state
847  ustar =(rcr*ur +rcl*ul + (ptotl-ptotr))/(rcr+rcl)
848  ptotstar=(rcr*ptotl+rcl*ptotr+rcl*rcr*(ul-ur))/(rcr+rcl)
849  !
850  ! Left star region variables
851  rstarl=rl*(sl-ul)/(sl-ustar)
852  etotstarl=((sl-ul)*etotl-ptotl*ul+ptotstar*ustar)/(sl-ustar)
853  !
854  ! Right star region variables
855  rstarr=rr*(sr-ur)/(sr-ustar)
856  etotstarr=((sr-ur)*etotr-ptotr*ur+ptotstar*ustar)/(sr-ustar)
857  !
858  ! Sample the solution at x/t=0
859  where (sl>0d0)
860  ro=rl
861  uo=ul
862  ptoto=ptotl
863  etoto=etotl
864  else where (ustar>0d0)
865  ro=rstarl
866  uo=ustar
867  ptoto=ptotstar
868  etoto=etotstarl
869  else where (sr>0d0)
870  ro=rstarr
871  uo=ustar
872  ptoto=ptotstar
873  etoto=etotstarr
874  else where
875  ro=rr
876  uo=ur
877  ptoto=ptotr
878  etoto=etotr
879  end where
880  !
881  ! Compute the Godunov flux
882  fgdnv(:,1) = ro*uo
883  fgdnv(:,2) = ro*uo*uo+ptoto
884  fgdnv(:,3) = (etoto+ptoto)*uo
885  do ivar = 4,nvar
886  where (fgdnv(:,1)>0)
887  fgdnv(:,ivar) = fgdnv(:,1)*qleft(:,ivar)
888  else where
889  fgdnv(:,ivar) = fgdnv(:,1)*qright(:,ivar)
890  end where
891  end do
892 end subroutine hllc_v
893 
894 !===============================================================================
895 !===============================================================================
896 subroutine hll(qleft,qright,fgdnv,ngrid,nvar)
897  integer::ngrid,nvar
898  real(dp),dimension(1:ngrid,1:nvar)::qleft,qright,fgdnv
899  real(dp),dimension(1:ngrid,1:nvar)::fleft,fright,uleft,uright
900  real(dp),dimension(1:ngrid)::sl,sr
901  integer::i,n
902  real(dp)::smallp, entho
903  real(dp)::rl ,ul ,pl ,cl
904  real(dp)::rr ,ur ,pr ,cr
905 
906  ! constants
907  smallp = smallc**2/gamma
908  entho = one/(gamma-one)
909 
910  ! Maximum wave speed
911  do i=3,ngrid-1
912  rl=max(qleft(i,1),smallr)
913  ul= qleft(i,2)
914  pl=max(qleft(i,3),rl*smallp)
915  rr=max(qright(i,1),smallr)
916  ur= qright(i,2)
917  pr=max(qright(i,3),rr*smallp)
918  cl= sqrt(gamma*pl/rl)
919  cr= sqrt(gamma*pr/rr)
920  sl(i)=min(min(ul,ur)-max(cl,cr),zero)
921  sr(i)=max(max(ul,ur)+max(cl,cr),zero)
922  end do
923 
924  ! Compute conservative variables
925  do i = 3,ngrid-1
926  uleft(i,1) = qleft(i,1)
927  uright(i,1) = qright(i,1)
928  uleft(i,2) = qleft(i,1)*qleft(i,2)
929  uright(i,2) = qright(i,1)*qright(i,2)
930  uleft(i,3) = qleft(i,3)*entho + half*qleft(i,1)*qleft(i,2)**2
931  uright(i,3) = qright(i,3)*entho + half*qright(i,1)*qright(i,2)**2
932  uleft(i,3) = uleft(i,3) + half*qleft(i,1)*qleft(i,4)**2
933  uright(i,3) = uright(i,3) + half*qright(i,1)*qright(i,4)**2
934  uleft(i,3) = uleft(i,3) + half*qleft(i,1)*qleft(i,5)**2
935  uright(i,3) = uright(i,3) + half*qright(i,1)*qright(i,5)**2
936  end do
937  ! Other advected quantities
938  do n = 4, nvar
939  do i = 3,ngrid-1
940  uleft(i,n) = qleft(i,1)*qleft(i,n)
941  uright(i,n) = qright(i,1)*qright(i,n)
942  end do
943  end do
944 
945  ! Compute left and right fluxes
946  do i = 3,ngrid-1
947  fleft(i,1) = uleft(i,2)
948  fright(i,1) = uright(i,2)
949  fleft(i,2) = qleft(i,3)+uleft(i,2)*qleft(i,2)
950  fright(i,2) = qright(i,3)+uright(i,2)*qright(i,2)
951  fleft(i,3) = qleft(i,2)*(uleft(i,3)+qleft(i,3))
952  fright(i,3) = qright(i,2)*(uright(i,3)+qright(i,3))
953  end do
954  ! Other advected quantities
955  do n = 4, nvar
956  do i = 3,ngrid-1
957  fleft(i,n) = fleft(i,1)*qleft(i,n)
958  fright(i,n) = fright(i,1)*qright(i,n)
959  end do
960  end do
961 
962  ! Compute HLL fluxes
963  do n = 1, nvar
964  do i = 3,ngrid-1
965  fgdnv(i,n) = (sr(i)*fleft(i,n)-sl(i)*fright(i,n) &
966  & + sr(i)*sl(i)*(uright(i,n)-uleft(i,n)))/(sr(i)-sl(i))
967  end do
968  end do
969 end subroutine hll
970 
971 end module riemann_mod
Definition: io_mod.f90:4