11 integer(8):: ndiff=0, nsolu=0
15 procedure,
nopass:: hlld_2d
16 procedure,
nopass:: llf
17 procedure,
nopass:: llf_2d
19 type(riemann_t),
public:: rieman
24 subroutine hlld (self, qleft, qright, fgdnv, ngrid, detailed_timer)
25 class(riemann_t):: self
27 real(dp),
dimension(:,:)::qleft,qright,fgdnv
28 LOGICAL,
OPTIONAL:: detailed_timer
31 REAL(dp)::rl,pl,ul,ecinl,etotl,ptotl,bl,cl,wstarl,vtotbl,emagl,sal,sqrrstarl,vl,wl
32 REAL(dp)::rr,pr,ur,ecinr,etotr,ptotr,br,cr,wstarr,vtotbr,emagr,sar,sqrrstarr,vr,wr
33 REAL(dp)::cfastl,rcl,rstarl,el,vdotbstarl,vstarl,etotstarstarl,etotstarl,vdotbl
34 REAL(dp)::cfastr,rcr,rstarr,er,vdotbstarr,vstarr,etotstarstarr,etotstarr,vdotbr
35 REAL(dp)::ustar,ptotstar,estar,vstarstar,wstarstar
36 REAL(dp)::ro,uo,ptoto,etoto,vo,wo
38 REAL(dp)::bstarl,bstarr,cstarl,cstarr,calfvenl,calfvenr
39 REAL(dp)::bo,co,bstarstar,cstarstar,vdotbstarstar,vdotbo,emago
42 REAL(dp)::d,p,u,v,w,b,c,b2,c2,d2,cf
44 INTEGER ::ivar,i,icase, ndiff
46 integer,
save:: itimer=0
54 entho = one/max(gamma-one,1e-6)
63 a=half*(qleft(i,4)+qright(i,4))
65 qleft(i,4)=a; qright(i,4)=a
68 rl=qleft(i,1); pl=qleft(i,2); ul=qleft(i,3)
69 vl=qleft(i,5); bl=qleft(i,6); wl=qleft(i,7); cl=qleft(i,8)
70 ecinl = half*(ul*ul+vl*vl+wl*wl)*rl
71 emagl = half*(a*a+bl*bl+cl*cl)
72 etotl = pl*entho+ecinl+emagl
74 vdotbl= ul*a+vl*bl+wl*cl
77 rr=qright(i,1); pr=qright(i,2); ur=qright(i,3)
78 vr=qright(i,5); br=qright(i,6); wr=qright(i,7); cr=qright(i,8)
79 ecinr = half*(ur*ur+vr*vr+wr*wr)*rr
80 emagr = half*(a*a+br*br+cr*cr)
81 etotr = pr*entho+ecinr+emagr
83 vdotbr= ur*a+vr*br+wr*cr
87 d=qleft(i,1); p=qleft(i,2); u=qleft(i,3); a=qleft(i,4)
88 v=qleft(i,5); b=qleft(i,6); w=qleft(i,7); c=qleft(i,8)
92 cf = sqrt( d2 + sqrt(max(d2*d2-c2*a*a/d,zero)) )
94 d=qright(i,1); p=qright(i,2); u=qright(i,3); a=qright(i,4)
95 v=qright(i,5); b=qright(i,6); w=qright(i,7); c=qright(i,8)
99 cf = sqrt( d2 + sqrt(max(d2*d2-c2*a*a/d,zero)) )
102 CALL find_speed_fast(qleft(i,:),cfastl)
103 CALL find_speed_fast(qright(i,:),cfastr)
107 sl=min(ul,ur)-max(cfastl,cfastr)
108 sr=max(ul,ur)+max(cfastl,cfastr)
117 ustar =(rcr*ur +rcl*ul + (ptotl-ptotr))/(rcr+rcl)
118 ptotstar=(rcr*ptotl+rcl*ptotr+rcl*rcr*(ul-ur))/(rcr+rcl)
121 rstarl=rl*(sl-ul)/(sl-ustar)
122 estar =rl*(sl-ul)*(sl-ustar)-a**2
123 el =rl*(sl-ul)*(sl-ul )-a**2
125 if(abs(estar) < 1d-4*a**2)
then 131 vstarl=vl-a*bl*(ustar-ul)/estar
133 wstarl=wl-a*cl*(ustar-ul)/estar
136 vdotbstarl=ustar*a+vstarl*bstarl+wstarl*cstarl
137 etotstarl=((sl-ul)*etotl-ptotl*ul+ptotstar*ustar+a*(vdotbl-vdotbstarl))/(sl-ustar)
138 sqrrstarl=sqrt(rstarl)
139 calfvenl=abs(a)/sqrrstarl
143 rstarr=rr*(sr-ur)/(sr-ustar)
144 estar =rr*(sr-ur)*(sr-ustar)-a**2
145 er =rr*(sr-ur)*(sr-ur )-a**2
146 if(abs(estar) < 1d-4*a**2)
then 152 vstarr=vr-a*br*(ustar-ur)/estar
154 wstarr=wr-a*cr*(ustar-ur)/estar
158 vdotbstarr=ustar*a+vstarr*bstarr+wstarr*cstarr
159 etotstarr=((sr-ur)*etotr-ptotr*ur+ptotstar*ustar+a*(vdotbr-vdotbstarr))/(sr-ustar)
160 sqrrstarr=sqrt(rstarr)
161 calfvenr=abs(a)/sqrrstarr
165 vstarstar=(sqrrstarl*vstarl+sqrrstarr*vstarr+sgnm*(bstarr-bstarl))/(sqrrstarl+sqrrstarr)
166 wstarstar=(sqrrstarl*wstarl+sqrrstarr*wstarr+sgnm*(cstarr-cstarl))/(sqrrstarl+sqrrstarr)
167 bstarstar=(sqrrstarl*bstarr+sqrrstarr*bstarl+sgnm*sqrrstarl*sqrrstarr*(vstarr-vstarl))/(sqrrstarl+sqrrstarr)
168 cstarstar=(sqrrstarl*cstarr+sqrrstarr*cstarl+sgnm*sqrrstarl*sqrrstarr*(wstarr-wstarl))/(sqrrstarl+sqrrstarr)
169 vdotbstarstar=ustar*a+vstarstar*bstarstar+wstarstar*cstarstar
170 etotstarstarl=etotstarl-sgnm*sqrrstarl*(vdotbstarl-vdotbstarstar)
171 etotstarstarr=etotstarr+sgnm*sqrrstarr*(vdotbstarr-vdotbstarstar)
197 emago=half*(a*a+bo*bo+co*co)
198 else if(ustar>0d0)
then 209 emago=half*(a*a+bo*bo+co*co)
221 emago=half*(a*a+bo*bo+co*co)
233 emago=half*(a*a+bo*bo+co*co)
250 ustar = alog(rl/rr)*max(0.0, alog10((10.*smallr)/min(rl,rr)))
251 ustar = 0.25*max(cfastr,cfastl)*min(1.0, max(-1.0, ustar))
253 ndiff = ndiff + merge(1,0,abs(ustar)>0.0)
257 fgdnv(i,2) = (etoto+ptoto)*uo-a*vdotbo
258 fgdnv(i,3) = ro*uo*uo+ptoto+emago*(lim-1.0_dp)-a*a*lim
260 fgdnv(i,5) = ro*uo*vo-a*bo*lim
261 fgdnv(i,6) = bo*uo-a*vo
262 fgdnv(i,7) = ro*uo*wo-a*co*lim
263 fgdnv(i,8) = co*uo-a*wo
266 self%ndiff = self%ndiff + ndiff
268 self%nsolu = self%nsolu + ngrid
275 SUBROUTINE hlld_2d (qLL,qRL,qLR,qRR,emf,j,k,l0,l1,nvar,xdim,detailed_timer)
276 INTEGER,
INTENT(IN) :: j,k,l0,l1,nvar,xdim
277 REAL(dp),
DIMENSION(:,:),
INTENT(IN)::qll,qrl,qlr,qrr
278 REAL(dp),
DIMENSION(:,:,:):: emf
279 LOGICAL,
OPTIONAL:: detailed_timer
282 REAL(dp) :: e,ell,erl,elr,err,s,sl,sr,sb,st,sal,sar,sat,sab
283 REAL(dp) :: cfastllx,cfastrlx,cfastlrx,cfastrrx,cfastlly,cfastrly,cfastlry,cfastrry
284 REAL(dp) :: calfvenr,calfvenl,calfvent,calfvenb
285 REAL(dp) :: vllx,vrlx,vlrx,vrrx,vlly,vrly,vlry,vrry
286 REAL(dp) :: ptotll,ptotlr,ptotrl,ptotrr,rcllx,rclrx,rcrlx,rcrrx,rclly,rclry,rcrly,rcrry
287 REAL(dp) :: ustar,vstar,rstarllx,rstarlrx,rstarrlx,rstarrrx,rstarlly,rstarlry,rstarrly,rstarrry
288 REAL(dp) :: rstarll,rstarlr,rstarrl,rstarrr,astarll,astarlr,astarrl,astarrr,bstarll,bstarlr,bstarrl,bstarrr
289 REAL(dp) :: estarllx,estarlrx,estarrlx,estarrrx,estarlly,estarlry,estarrly,estarrry,estarll,estarlr,estarrl,estarrr
290 REAL(dp) :: astart,astarb,bstarr,bstarl
292 REAL(dp) :: rll,rlr,rrl,rrr,pll,plr,prl,prr,ull,ulr,url,urr,vll,vlr,vrl,vrr, &
293 all,alr,arl,arr,bll,blr,brl,brr,cll,clr,crl,crr
294 REAL(dp),
DIMENSION(1:nvar) :: qtmp
295 REAL(dp) :: b2, c2, d2, d1
297 integer,
save:: itimer=0
302 rll=qll(l,1); pll=qll(l,2); ull=qll(l,3); vll=qll(l,4); all=qll(l,6); bll=qll(l,7) ; cll=qll(l,8)
303 rlr=qlr(l,1); plr=qlr(l,2); ulr=qlr(l,3); vlr=qlr(l,4); alr=qlr(l,6); blr=qlr(l,7) ; clr=qlr(l,8)
304 rrl=qrl(l,1); prl=qrl(l,2); url=qrl(l,3); vrl=qrl(l,4); arl=qrl(l,6); brl=qrl(l,7) ; crl=qrl(l,8)
305 rrr=qrr(l,1); prr=qrr(l,2); urr=qrr(l,3); vrr=qrr(l,4); arr=qrr(l,6); brr=qrr(l,7) ; crr=qrr(l,8)
308 qtmp(1)=qll(l,1); qtmp(2)=qll(l,2); qtmp(7)=qll(l,5); qtmp(8)=qll(l,8)
309 qtmp(3)=qll(l,3); qtmp(4)=qll(l,6); qtmp(5)=qll(l,4); qtmp(6)=qll(l,7)
312 b2 = qtmp(4)**2+qtmp(6)**2+qtmp(8)**2
314 c2 = gamma*qtmp(2)*d1
316 cfastllx=sqrt(d2 + sqrt(max(d2*d2-c2*qtmp(4)**2*d1,zero)))
318 qtmp(1)=qlr(l,1); qtmp(2)=qlr(l,2); qtmp(7)=qlr(l,5); qtmp(8)=qlr(l,8)
319 qtmp(3)=qlr(l,3); qtmp(4)=qlr(l,6); qtmp(5)=qlr(l,4); qtmp(6)=qlr(l,7)
322 b2 = qtmp(4)**2+qtmp(6)**2+qtmp(8)**2
324 c2 = gamma*qtmp(2)*d1
326 cfastlrx=sqrt(d2 + sqrt(max(d2*d2-c2*qtmp(4)**2*d1,zero)))
328 qtmp(1)=qrl(l,1); qtmp(2)=qrl(l,2); qtmp(7)=qrl(l,5); qtmp(8)=qrl(l,8)
329 qtmp(3)=qrl(l,3); qtmp(4)=qrl(l,6); qtmp(5)=qrl(l,4); qtmp(6)=qrl(l,7)
332 b2 = qtmp(4)**2+qtmp(6)**2+qtmp(8)**2
334 c2 = gamma*qtmp(2)*d1
337 cfastrlx=sqrt(d2 + sqrt(max(d2*d2-c2*qtmp(4)**2*d1,zero)))
339 qtmp(1)=qrr(l,1); qtmp(2)=qrr(l,2); qtmp(7)=qrr(l,5); qtmp(8)=qrr(l,8)
340 qtmp(3)=qrr(l,3); qtmp(4)=qrr(l,6); qtmp(5)=qrr(l,4); qtmp(6)=qrr(l,7)
341 qtmp(1)=max(qtmp(1),smallr)
344 b2 = qtmp(4)**2+qtmp(6)**2+qtmp(8)**2
346 c2 = gamma*qtmp(2)*d1
348 cfastrrx=sqrt(d2 + sqrt(max(d2*d2-c2*qtmp(4)**2*d1,zero)))
351 qtmp(1)=qll(l,1); qtmp(2)=qll(l,2); qtmp(7)=qll(l,5); qtmp(8)=qll(l,8)
352 qtmp(3)=qll(l,4); qtmp(4)=qll(l,7); qtmp(5)=qll(l,3); qtmp(6)=qll(l,6)
355 b2 = qtmp(4)**2+qtmp(6)**2+qtmp(8)**2
357 c2 = gamma*qtmp(2)*d1
359 cfastlly=sqrt(d2 + sqrt(max(d2*d2-c2*qtmp(4)**2*d1,zero)))
361 qtmp(1)=qlr(l,1); qtmp(2)=qlr(l,2); qtmp(7)=qlr(l,5); qtmp(8)=qlr(l,8)
362 qtmp(3)=qlr(l,4); qtmp(4)=qlr(l,7); qtmp(5)=qlr(l,3); qtmp(6)=qlr(l,6)
365 b2 = qtmp(4)**2+qtmp(6)**2+qtmp(8)**2
367 c2 = gamma*qtmp(2)*d1
369 cfastlry=sqrt(d2 + sqrt(max(d2*d2-c2*qtmp(4)**2*d1,zero)))
371 qtmp(1)=qrl(l,1); qtmp(2)=qrl(l,2); qtmp(7)=qrl(l,5); qtmp(8)=qrl(l,8)
372 qtmp(3)=qrl(l,4); qtmp(4)=qrl(l,7); qtmp(5)=qrl(l,3); qtmp(6)=qrl(l,6)
375 b2 = qtmp(4)**2+qtmp(6)**2+qtmp(8)**2
377 c2 = gamma*qtmp(2)*d1
379 cfastrly=sqrt(d2 + sqrt(max(d2*d2-c2*qtmp(4)**2*d1,zero)))
381 qtmp(1)=qrr(l,1); qtmp(2)=qrr(l,2); qtmp(7)=qrr(l,5); qtmp(8)=qrr(l,8)
382 qtmp(3)=qrr(l,4); qtmp(4)=qrr(l,7); qtmp(5)=qrr(l,3); qtmp(6)=qrr(l,6)
385 b2 = qtmp(4)**2+qtmp(6)**2+qtmp(8)**2
387 c2 = gamma*qtmp(2)*d1
389 cfastrry=sqrt(d2 + sqrt(max(d2*d2-c2*qtmp(4)**2*d1,zero)))
391 sl=min(ull,ulr,url,urr)-max(cfastllx,cfastlrx,cfastrlx,cfastrrx)
392 sr=max(ull,ulr,url,urr)+max(cfastllx,cfastlrx,cfastrlx,cfastrrx)
393 sb=min(vll,vlr,vrl,vrr)-max(cfastlly,cfastlry,cfastrly,cfastrry)
394 st=max(vll,vlr,vrl,vrr)+max(cfastlly,cfastlry,cfastrly,cfastrry)
412 ptotll=pll+half*(all*all+bll*bll+cll*cll)
413 ptotlr=plr+half*(alr*alr+blr*blr+clr*clr)
414 ptotrl=prl+half*(arl*arl+brl*brl+crl*crl)
415 ptotrr=prr+half*(arr*arr+brr*brr+crr*crr)
417 rcllx=rll*(ull-sl); rcrlx=rrl*(sr-url)
418 rclrx=rlr*(ulr-sl); rcrrx=rrr*(sr-urr)
419 rclly=rll*(vll-sb); rclry=rlr*(st-vlr)
420 rcrly=rrl*(vrl-sb); rcrry=rrr*(st-vrr)
422 ustar=(rcllx*ull+rclrx*ulr+rcrlx*url+rcrrx*urr+(ptotll-ptotrl+ptotlr-ptotrr))/(rcllx+rclrx+rcrlx+rcrrx)
423 vstar=(rclly*vll+rclry*vlr+rcrly*vrl+rcrry*vrr+(ptotll-ptotlr+ptotrl-ptotrr))/(rclly+rclry+rcrly+rcrry)
425 rstarllx=rll*(sl-ull)/(sl-ustar)
426 rstarll =rll*(sl-ull)/(sl-ustar)*(sb-vll)/(sb-vstar)
427 astarll=all*(sb-vll)/(sb-vstar)
429 rstarlrx=rlr*(sl-ulr)/(sl-ustar)
430 rstarlr =rlr*(sl-ulr)/(sl-ustar)*(st-vlr)/(st-vstar)
431 astarlr=alr*(st-vlr)/(st-vstar)
433 rstarrlx=rrl*(sr-url)/(sr-ustar)
434 rstarrl =rrl*(sr-url)/(sr-ustar)*(sb-vrl)/(sb-vstar)
435 astarrl=arl*(sb-vrl)/(sb-vstar)
437 rstarrrx=rrr*(sr-urr)/(sr-ustar)
438 rstarrr =rrr*(sr-urr)/(sr-ustar)*(st-vrr)/(st-vstar)
439 astarrr=arr*(st-vrr)/(st-vstar)
441 calfvenl=max(abs(alr)/sqrt(rstarlrx),abs(astarlr)/sqrt(rstarlr), &
442 abs(all)/sqrt(rstarllx),abs(astarll)/sqrt(rstarll),smallc)
443 calfvenr=max(abs(arr)/sqrt(rstarrrx),abs(astarrr)/sqrt(rstarrr), &
444 abs(arl)/sqrt(rstarrlx),abs(astarrl)/sqrt(rstarrl),smallc)
446 sal=min(ustar-calfvenl,zero)
447 sar=max(ustar+calfvenr,zero)
449 bstarll=bll*(sl-ull)/(sl-ustar)
450 bstarrl=brl*(sr-url)/(sr-ustar)
452 estarllx=ustar*bstarll-vll*all
453 estarrlx=ustar*bstarrl-vrl*arl
455 e=(sar*estarllx-sal*estarrlx+sar*sal*(brl-bll))/(sar-sal)
473 ptotll=pll+half*(all*all+bll*bll+cll*cll)
474 ptotlr=plr+half*(alr*alr+blr*blr+clr*clr)
475 ptotrl=prl+half*(arl*arl+brl*brl+crl*crl)
476 ptotrr=prr+half*(arr*arr+brr*brr+crr*crr)
478 rcllx=rll*(ull-sl); rcrlx=rrl*(sr-url)
479 rclrx=rlr*(ulr-sl); rcrrx=rrr*(sr-urr)
480 rclly=rll*(vll-sb); rclry=rlr*(st-vlr)
481 rcrly=rrl*(vrl-sb); rcrry=rrr*(st-vrr)
483 ustar=(rcllx*ull+rclrx*ulr+rcrlx*url+rcrrx*urr+(ptotll-ptotrl+ptotlr-ptotrr))/(rcllx+rclrx+rcrlx+rcrrx)
484 vstar=(rclly*vll+rclry*vlr+rcrly*vrl+rcrry*vrr+(ptotll-ptotlr+ptotrl-ptotrr))/(rclly+rclry+rcrly+rcrry)
486 rstarllx=rll*(sl-ull)/(sl-ustar)
487 rstarll =rll*(sl-ull)/(sl-ustar)*(sb-vll)/(sb-vstar)
488 astarll=all*(sb-vll)/(sb-vstar)
490 rstarlrx=rlr*(sl-ulr)/(sl-ustar)
491 rstarlr =rlr*(sl-ulr)/(sl-ustar)*(st-vlr)/(st-vstar)
492 astarlr=alr*(st-vlr)/(st-vstar)
494 rstarrlx=rrl*(sr-url)/(sr-ustar)
495 rstarrl =rrl*(sr-url)/(sr-ustar)*(sb-vrl)/(sb-vstar)
496 astarrl=arl*(sb-vrl)/(sb-vstar)
498 rstarrrx=rrr*(sr-urr)/(sr-ustar)
499 rstarrr =rrr*(sr-urr)/(sr-ustar)*(st-vrr)/(st-vstar)
500 astarrr=arr*(st-vrr)/(st-vstar)
502 calfvenl=max(abs(alr)/sqrt(rstarlrx),abs(astarlr)/sqrt(rstarlr), &
503 abs(all)/sqrt(rstarllx),abs(astarll)/sqrt(rstarll),smallc)
504 calfvenr=max(abs(arr)/sqrt(rstarrrx),abs(astarrr)/sqrt(rstarrr), &
505 abs(arl)/sqrt(rstarrlx),abs(astarrl)/sqrt(rstarrl),smallc)
507 sal=min(ustar-calfvenl,zero)
508 sar=max(ustar+calfvenr,zero)
510 bstarlr=blr*(sl-ulr)/(sl-ustar)
511 bstarrr=brr*(sr-urr)/(sr-ustar)
513 estarlrx=ustar*bstarlr-vlr*alr
514 estarrrx=ustar*bstarrr-vrr*arr
516 e=(sar*estarlrx-sal*estarrrx+sar*sal*(brr-blr))/(sar-sal)
521 ptotll=pll+half*(all*all+bll*bll+cll*cll)
522 ptotlr=plr+half*(alr*alr+blr*blr+clr*clr)
523 ptotrl=prl+half*(arl*arl+brl*brl+crl*crl)
524 ptotrr=prr+half*(arr*arr+brr*brr+crr*crr)
526 rcllx=rll*(ull-sl); rcrlx=rrl*(sr-url)
527 rclrx=rlr*(ulr-sl); rcrrx=rrr*(sr-urr)
528 rclly=rll*(vll-sb); rclry=rlr*(st-vlr)
529 rcrly=rrl*(vrl-sb); rcrry=rrr*(st-vrr)
531 ustar=(rcllx*ull+rclrx*ulr+rcrlx*url+rcrrx*urr+(ptotll-ptotrl+ptotlr-ptotrr))/(rcllx+rclrx+rcrlx+rcrrx)
532 vstar=(rclly*vll+rclry*vlr+rcrly*vrl+rcrry*vrr+(ptotll-ptotlr+ptotrl-ptotrr))/(rclly+rclry+rcrly+rcrry)
534 rstarlly=rll*(sb-vll)/(sb-vstar)
535 rstarll =rll*(sl-ull)/(sl-ustar)*(sb-vll)/(sb-vstar)
536 bstarll=bll*(sl-ull)/(sl-ustar)
538 rstarlry=rlr*(st-vlr)/(st-vstar)
539 rstarlr =rlr*(sl-ulr)/(sl-ustar)*(st-vlr)/(st-vstar)
540 bstarlr=blr*(sl-ulr)/(sl-ustar)
542 rstarrly=rrl*(sb-vrl)/(sb-vstar)
543 rstarrl =rrl*(sr-url)/(sr-ustar)*(sb-vrl)/(sb-vstar)
544 bstarrl=brl*(sr-url)/(sr-ustar)
546 rstarrry=rrr*(st-vrr)/(st-vstar)
547 rstarrr =rrr*(sr-urr)/(sr-ustar)*(st-vrr)/(st-vstar)
548 bstarrr=brr*(sr-urr)/(sr-ustar)
550 calfvenb=max(abs(bll)/sqrt(rstarlly),abs(bstarll)/sqrt(rstarll), &
551 abs(brl)/sqrt(rstarrly),abs(bstarrl)/sqrt(rstarrl),smallc)
552 calfvent=max(abs(blr)/sqrt(rstarlry),abs(bstarlr)/sqrt(rstarlr), &
553 abs(brr)/sqrt(rstarrry),abs(bstarrr)/sqrt(rstarrr),smallc)
555 sab=min(vstar-calfvenb,zero)
556 sat=max(vstar+calfvent,zero)
558 astarll=all*(sb-vll)/(sb-vstar)
559 astarlr=alr*(st-vlr)/(st-vstar)
561 estarlly=ull*bll-vstar*astarll
562 estarlry=ulr*blr-vstar*astarlr
564 e=(sat*estarlly-sab*estarlry-sat*sab*(alr-all))/(sat-sab)
568 ptotll=pll+half*(all*all+bll*bll+cll*cll)
569 ptotlr=plr+half*(alr*alr+blr*blr+clr*clr)
570 ptotrl=prl+half*(arl*arl+brl*brl+crl*crl)
571 ptotrr=prr+half*(arr*arr+brr*brr+crr*crr)
573 rcllx=rll*(ull-sl); rcrlx=rrl*(sr-url)
574 rclrx=rlr*(ulr-sl); rcrrx=rrr*(sr-urr)
575 rclly=rll*(vll-sb); rclry=rlr*(st-vlr)
576 rcrly=rrl*(vrl-sb); rcrry=rrr*(st-vrr)
578 ustar=(rcllx*ull+rclrx*ulr+rcrlx*url+rcrrx*urr+(ptotll-ptotrl+ptotlr-ptotrr))/(rcllx+rclrx+rcrlx+rcrrx)
579 vstar=(rclly*vll+rclry*vlr+rcrly*vrl+rcrry*vrr+(ptotll-ptotlr+ptotrl-ptotrr))/(rclly+rclry+rcrly+rcrry)
581 rstarlly=rll*(sb-vll)/(sb-vstar)
582 rstarll =rll*(sl-ull)/(sl-ustar)*(sb-vll)/(sb-vstar)
583 bstarll=bll*(sl-ull)/(sl-ustar)
585 rstarlry=rlr*(st-vlr)/(st-vstar)
586 rstarlr =rlr*(sl-ulr)/(sl-ustar)*(st-vlr)/(st-vstar)
587 bstarlr=blr*(sl-ulr)/(sl-ustar)
589 rstarrly=rrl*(sb-vrl)/(sb-vstar)
590 rstarrl =rrl*(sr-url)/(sr-ustar)*(sb-vrl)/(sb-vstar)
591 bstarrl=brl*(sr-url)/(sr-ustar)
593 rstarrry=rrr*(st-vrr)/(st-vstar)
594 rstarrr =rrr*(sr-urr)/(sr-ustar)*(st-vrr)/(st-vstar)
595 bstarrr=brr*(sr-urr)/(sr-ustar)
597 calfvenb=max(abs(bll)/sqrt(rstarlly),abs(bstarll)/sqrt(rstarll), &
598 abs(brl)/sqrt(rstarrly),abs(bstarrl)/sqrt(rstarrl),smallc)
599 calfvent=max(abs(blr)/sqrt(rstarlry),abs(bstarlr)/sqrt(rstarlr), &
600 abs(brr)/sqrt(rstarrry),abs(bstarrr)/sqrt(rstarrr),smallc)
602 sab=min(vstar-calfvenb,zero)
603 sat=max(vstar+calfvent,zero)
605 astarrl=arl*(sb-vrl)/(sb-vstar)
606 astarrr=arr*(st-vrr)/(st-vstar)
608 estarrly=url*brl-vstar*astarrl
609 estarrry=urr*brr-vstar*astarrr
611 e=(sat*estarrly-sab*estarrry-sat*sab*(arr-arl))/(sat-sab)
615 ptotll=pll+half*(all*all+bll*bll+cll*cll)
616 ptotlr=plr+half*(alr*alr+blr*blr+clr*clr)
617 ptotrl=prl+half*(arl*arl+brl*brl+crl*crl)
618 ptotrr=prr+half*(arr*arr+brr*brr+crr*crr)
620 rcllx=rll*(ull-sl); rcrlx=rrl*(sr-url)
621 rclrx=rlr*(ulr-sl); rcrrx=rrr*(sr-urr)
622 rclly=rll*(vll-sb); rclry=rlr*(st-vlr)
623 rcrly=rrl*(vrl-sb); rcrry=rrr*(st-vrr)
625 ustar=(rcllx*ull+rclrx*ulr+rcrlx*url+rcrrx*urr+(ptotll-ptotrl+ptotlr-ptotrr))/(rcllx+rclrx+rcrlx+rcrrx)
626 vstar=(rclly*vll+rclry*vlr+rcrly*vrl+rcrry*vrr+(ptotll-ptotlr+ptotrl-ptotrr))/(rclly+rclry+rcrly+rcrry)
628 rstarllx=rll*(sl-ull)/(sl-ustar); bstarll=bll*(sl-ull)/(sl-ustar)
629 rstarlly=rll*(sb-vll)/(sb-vstar); astarll=all*(sb-vll)/(sb-vstar)
630 rstarll =rll*(sl-ull)/(sl-ustar)*(sb-vll)/(sb-vstar)
632 rstarlrx=rlr*(sl-ulr)/(sl-ustar); bstarlr=blr*(sl-ulr)/(sl-ustar)
633 rstarlry=rlr*(st-vlr)/(st-vstar); astarlr=alr*(st-vlr)/(st-vstar)
634 rstarlr =rlr*(sl-ulr)/(sl-ustar)*(st-vlr)/(st-vstar)
636 rstarrlx=rrl*(sr-url)/(sr-ustar); bstarrl=brl*(sr-url)/(sr-ustar)
637 rstarrly=rrl*(sb-vrl)/(sb-vstar); astarrl=arl*(sb-vrl)/(sb-vstar)
638 rstarrl =rrl*(sr-url)/(sr-ustar)*(sb-vrl)/(sb-vstar)
640 rstarrrx=rrr*(sr-urr)/(sr-ustar); bstarrr=brr*(sr-urr)/(sr-ustar)
641 rstarrry=rrr*(st-vrr)/(st-vstar); astarrr=arr*(st-vrr)/(st-vstar)
642 rstarrr =rrr*(sr-urr)/(sr-ustar)*(st-vrr)/(st-vstar)
644 calfvenl=max(abs(alr)/sqrt(rstarlrx),abs(astarlr)/sqrt(rstarlr), &
645 abs(all)/sqrt(rstarllx),abs(astarll)/sqrt(rstarll),smallc)
646 calfvenr=max(abs(arr)/sqrt(rstarrrx),abs(astarrr)/sqrt(rstarrr), &
647 abs(arl)/sqrt(rstarrlx),abs(astarrl)/sqrt(rstarrl),smallc)
648 calfvenb=max(abs(bll)/sqrt(rstarlly),abs(bstarll)/sqrt(rstarll), &
649 abs(brl)/sqrt(rstarrly),abs(bstarrl)/sqrt(rstarrl),smallc)
650 calfvent=max(abs(blr)/sqrt(rstarlry),abs(bstarlr)/sqrt(rstarlr), &
651 abs(brr)/sqrt(rstarrry),abs(bstarrr)/sqrt(rstarrr),smallc)
652 sal=min(ustar-calfvenl,zero)
653 sar=max(ustar+calfvenr,zero)
654 sab=min(vstar-calfvenb,zero)
655 sat=max(vstar+calfvent,zero)
656 estarll =ustar*bstarll-vstar*astarll
657 estarrl =ustar*bstarrl-vstar*astarrl
658 estarlr =ustar*bstarlr-vstar*astarlr
659 estarrr =ustar*bstarrr-vstar*astarrr
661 astart=(sar*astarrr-sal*astarlr)/(sar-sal); astarb=(sar*astarrl-sal*astarll)/(sar-sal)
662 bstarr=(sat*bstarrr-sab*bstarrl)/(sat-sab); bstarl=(sat*bstarlr-sab*bstarll)/(sat-sab)
664 e=(sal*sab*estarrr-sal*sat*estarrl-sar*sab*estarlr+sar*sat*estarll)/(sar-sal)/(sat-sab) &
665 & -sat*sab/(sat-sab)*(astart-astarb)+sar*sal/(sar-sal)*(bstarr-bstarl)
670 else if (xdim == 2)
then 672 else if (xdim == 3)
then 677 END SUBROUTINE hlld_2d
682 SUBROUTINE llf_2d (qLL,qRL,qLR,qRR,emf,j,k,l0,l1,nvar,xdim,detailed_timer)
683 INTEGER,
INTENT(IN) :: j,k,l0,l1,nvar,xdim
684 REAL(dp),
DIMENSION(:,:),
INTENT(IN)::qll,qrl,qlr,qrr
685 REAL(dp),
DIMENSION(:,:,:):: emf
686 LOGICAL,
OPTIONAL:: detailed_timer
689 REAL(dp) :: e, ev, ell, erl, elr, err, qleft(8), qright(8), fmean_x(8), fmean_y(8)
690 integer,
save:: itimer=0
696 ell = qll(l,3)*qll(l,7) - qll(l,4)*qll(l,6)
697 erl = qrl(l,3)*qrl(l,7) - qrl(l,4)*qrl(l,6)
698 elr = qlr(l,3)*qlr(l,7) - qlr(l,4)*qlr(l,6)
699 err = qrr(l,3)*qrr(l,7) - qrr(l,4)*qrr(l,6)
702 ev = forth*(ell+erl+elr+err)
706 qleft(1) = half*(qll(l,1)+qlr(l,1))
707 qright(1) = half*(qrr(l,1)+qrl(l,1))
710 qleft(2) = half*(qll(l,2)+qlr(l,2))
711 qright(2) = half*(qrr(l,2)+qrl(l,2))
714 qleft(3) = half*(qll(l,3)+qlr(l,3))
715 qright(3) = half*(qrr(l,3)+qrl(l,3))
718 qleft(4) = half*(qll(l,6)+qlr(l,6))
719 qright(4) = half*(qrr(l,6)+qrl(l,6))
722 qleft(5) = half*(qll(l,4)+qlr(l,4))
723 qright(5) = half*(qrr(l,4)+qrl(l,4))
726 qleft(6) = half*(qll(l,7)+qlr(l,7))
727 qright(6) = half*(qrr(l,7)+qrl(l,7))
730 qleft(7) = half*(qll(l,5)+qlr(l,5))
731 qright(7) = half*(qrr(l,5)+qrl(l,5))
734 qleft(8) = half*(qll(l,8)+qlr(l,8))
735 qright(8) = half*(qrr(l,8)+qrl(l,8))
737 CALL lax_friedrich(qleft,qright,fmean_x)
741 qleft(1) = half*(qll(l,1)+qrl(l,1))
742 qright(1) = half*(qrr(l,1)+qlr(l,1))
745 qleft(2) = half*(qll(l,2)+qrl(l,2))
746 qright(2) = half*(qrr(l,2)+qlr(l,2))
749 qleft(3) = half*(qll(l,4)+qrl(l,4))
750 qright(3) = half*(qrr(l,4)+qlr(l,4))
753 qleft(4) = half*(qll(l,7)+qrl(l,7))
754 qright(4) = half*(qrr(l,7)+qlr(l,7))
757 qleft(5) = half*(qll(l,3)+qrl(l,3))
758 qright(5) = half*(qrr(l,3)+qlr(l,3))
761 qleft(6) = half*(qll(l,6)+qrl(l,6))
762 qright(6) = half*(qrr(l,6)+qlr(l,6))
765 qleft(7) = half*(qll(l,5)+qrl(l,5))
766 qright(7) = half*(qrr(l,5)+qlr(l,5))
769 qleft(8) = half*(qll(l,8)+qrl(l,8))
770 qright(8) = half*(qrr(l,8)+qlr(l,8))
772 CALL lax_friedrich(qleft,qright,fmean_y)
774 e = ev + (fmean_x(6) - fmean_y(6))
777 else if (xdim == 2)
then 779 else if (xdim == 3)
then 784 END SUBROUTINE llf_2d
789 SUBROUTINE llf (qleft, qright, fgdnv, ngrid, detailed_timer)
790 real(DP),
dimension(:,:):: qleft, qright, fgdnv
792 logical,
optional:: detailed_timer
795 call lax_friedrich (qleft(l,:), qright(l,:), fgdnv(l,:))
802 SUBROUTINE lax_friedrich (qleft, qright, fgdnv)
803 REAL(dp),
DIMENSION(:):: qleft, qright, fgdnv
804 REAL(dp),
DIMENSION(size(qleft)):: fleft, fright, fmean
805 REAL(dp),
DIMENSION(size(qleft)):: uleft ,uright, udiff
806 REAL(dp):: vleft, vright, bx_mean
809 bx_mean = half*(qleft(4)+qright(4))
813 CALL find_mhd_flux (qleft ,uleft ,fleft )
814 CALL find_mhd_flux (qright,uright,fright)
817 fmean = half * ( fright + fleft )
820 CALL find_speed_info (qleft ,vleft )
821 CALL find_speed_info (qright,vright)
824 udiff = half * ( uright - uleft )
827 fgdnv = fmean - max(vleft,vright) * udiff
828 END SUBROUTINE lax_friedrich
832 SUBROUTINE find_mhd_flux (qvar, cvar, ff)
837 INTEGER :: i , ivar, ib
838 REAL(dp),
DIMENSION(:):: qvar, cvar, ff
839 REAL(dp):: ecin,emag,etot,d,u,v,w,a,b,c,p,ptot,entho
842 entho = one/(gamma-one)
846 entho = one/max(gamma-one,1e-6)
848 d=qvar(1); p=qvar(2); u=qvar(3); a=qvar(4)
849 v=qvar(5); b=qvar(6); w=qvar(7); c=qvar(8)
850 ecin = half*(u*u+v*v+w*w)*d
851 emag = half*(a*a+b*b+c*c)
852 etot = p*entho+ecin+emag
867 ff(2) = (etot+ptot)*u-a*(a*u+b*v+c*w)
868 ff(3) = d*u*u+ptot-a*a
874 END SUBROUTINE find_mhd_flux
878 SUBROUTINE find_speed_fast(qvar,vel_info)
882 REAL(dp),
DIMENSION(:):: qvar
884 REAL(dp) :: d,p,u,v,w,a,b,c,b2,c2,d2,cf
885 integer,
save:: itimer=0
887 d=qvar(1); p=qvar(2); u=qvar(3); a=qvar(4)
888 v=qvar(5); b=qvar(6); w=qvar(7); c=qvar(8)
890 p=max(p,smallr*smallc*smallc)
894 cf = sqrt( d2 + sqrt(max(d2*d2-c2*a*a/d,zero)) )
896 END SUBROUTINE find_speed_fast
900 SUBROUTINE find_speed_info (qvar, vel_info)
906 REAL(dp),
DIMENSION(:):: qvar
908 REAL(dp):: d, p, u, v, w, a, b, c, b2, c2, d2, cf
910 d=qvar(1); p=qvar(2); u=qvar(3); a=qvar(4)
911 v=qvar(5); b=qvar(6); w=qvar(7); c=qvar(8)
915 cf = sqrt( d2 + sqrt(d2**2-c2*a*a/d) )
917 END SUBROUTINE find_speed_info
919 END MODULE riemann_mod