19 subroutine hancock(qdt,dtfactor,ixI^L,ixO^L,idims^LIM,qtC,sCT,qt,snew,dxs,x)
25 integer,
intent(in) :: ixi^
l, ixo^
l, idims^lim
26 double precision,
intent(in) :: qdt, dtfactor,qtc, qt, dxs(
ndim), x(ixi^s,1:
ndim)
27 type(state) :: sct, snew
29 double precision,
dimension(ixI^S,1:nw) :: wprim, wlc, wrc
31 double precision,
dimension(ixI^S,1:nw) :: wlp, wrp
32 double precision,
dimension(ixO^S) :: inv_volume
33 double precision :: flc(ixi^s, nwflux), frc(ixi^s, nwflux)
34 double precision :: dxinv(1:
ndim)
35 integer :: idims, iw, ix^
l, hxo^
l
38 associate(wct=>sct%w,wnew=>snew%w)
42 ix^
l=ix^
l^laddkr(idims,^
d);
44 if (ixi^
l^ltix^
l|.or.|.or.) &
45 call mpistop(
"Error in Hancock: Nonconforming input limits")
59 hxo^
l=ixo^
l-
kr(idims,^
d);
61 wrp(hxo^s,1:nwflux)=wprim(ixo^s,1:nwflux)
62 wlp(ixo^s,1:nwflux)=wprim(ixo^s,1:nwflux)
65 call reconstruct_lr(ixi^
l,ixo^
l,hxo^
l,idims,wprim,wlc,wrc,wlp,wrp,x,dxs(idims))
75 wnew(ixo^s,iw)=wnew(ixo^s,iw)-
block%dt(ixo^s)*dtfactor/dxs(idims)* &
76 (flc(ixo^s, iw)-frc(hxo^s, iw))
80 wnew(ixo^s,iw)=wnew(ixo^s,iw)+dxinv(idims)* &
81 (flc(ixo^s, iw)-frc(hxo^s, iw))
87 wnew(ixo^s,iw)=wnew(ixo^s,iw) -
block%dt(ixo^s)*dtfactor * inv_volume &
88 *(
block%surfaceC(ixo^s,idims)*flc(ixo^s, iw) &
89 -
block%surfaceC(hxo^s,idims)*frc(hxo^s, iw))
93 wnew(ixo^s,iw)=wnew(ixo^s,iw) - qdt * inv_volume &
94 *(
block%surfaceC(ixo^s,idims)*flc(ixo^s, iw) &
95 -
block%surfaceC(hxo^s,idims)*frc(hxo^s, iw))
105 dtfactor*dble(idimsmax-idimsmin+1)/dble(
ndim),&
106 ixi^
l,ixo^
l,1,nw,qtc,wct,wprim,qt,wnew,x,.false.,active)
117 qtC,sCT,qt,snew,fC,fE,dxs,x)
126 integer,
intent(in) :: method
127 double precision,
intent(in) :: qdt, dtfactor, qtc, qt, dxs(
ndim)
128 integer,
intent(in) :: ixi^
l, ixo^
l, idims^lim
129 double precision,
dimension(ixI^S,1:ndim),
intent(in) :: x
130 type(state) :: sct, snew
131 double precision,
dimension(ixI^S,1:nwflux,1:ndim) :: fc
132 double precision,
dimension(ixI^S,sdim:3) :: fe
136 integer :: idims, iw, ix^
d, hx^
d, ix^
l, hxo^
l, ixc^
l, ixcr^
l, kxc^
l, kxr^
l, ii
137 integer,
dimension(ixI^S) :: patchf
138 double precision,
dimension(ixI^S,1:nw) :: wprim
140 double precision,
dimension(ixI^S,1:nw) :: wlc, wrc
142 double precision,
dimension(ixI^S,1:nw) :: wlp, wrp
143 double precision,
dimension(ixI^S,1:nwflux) :: flc, frc
144 double precision,
dimension(ixI^S,1:number_species) :: cmaxc
145 double precision,
dimension(ixI^S,1:number_species) :: cminc
146 double precision,
dimension(ixI^S) :: hspeed
147 double precision,
dimension(ixO^S) :: inv_volume
148 double precision,
dimension(1:ndim) :: dxinv
149 type(ct_velocity) :: vcts
151 associate(wct=>sct%w, wnew=>snew%w)
163 ix^
l=ix^
l^ladd2*
kr(idims,^
d);
165 if (ixi^
l^ltix^
l|.or.|.or.) &
166 call mpistop(
"Error in fv : Nonconforming input limits")
175 hxo^
l=ixo^
l-
kr(idims,^
d);
177 kxcmin^
d=iximin^
d; kxcmax^
d=iximax^
d-
kr(idims,^
d);
178 kxr^
l=kxc^
l+
kr(idims,^
d);
186 ixcmax^
d=ixomax^
d; ixcmin^
d=hxomin^
d;
193 wrp(kxc^s,1:
nw)=wprim(kxr^s,1:
nw)
194 wlp(kxc^s,1:
nw)=wprim(kxc^s,1:
nw)
200 call reconstruct_lr(ixi^
l,ixcr^
l,ixcr^
l,idims,wprim,wlc,wrc,wlp,wrp,x,dxs(idims))
218 call phys_get_cbounds(wlc,wrc,wlp,wrp,x,ixi^
l,ixc^
l,idims,hspeed,cmaxc,cminc)
247 call mpistop(
'unkown Riemann flux in finite volume')
257 hxomin^
d=ixomin^
d-
kr(idims,^
d)\
259 {
do ix^db=hxomin^db,ixomax^db\}
260 fc(ix^
d,iw,idims)=
block%dt(ix^
d)*dxinv(idims)*fc(ix^
d,iw,idims)
262 {
do ix^db=ixomin^db,ixomax^db\}
263 hx^d=ix^d-kr(idims,^d)\
264 wnew(ix^d,iw)=wnew(ix^d,iw)+fc(ix^d,iw,idims)-fc(hx^d,iw,idims)
268 if(method==fs_tvdmu) &
269 call tvdlimit2(method,qdt,ixi^l,ixc^l,ixo^l,idims,wlc,wrc,wnew,x,fc,dxs)
272 dxinv(1:ndim)=-qdt/dxs(1:ndim)
274 hxomin^d=ixomin^d-kr(idims,^d)\
276 {
do ix^db=hxomin^db,ixomax^db\}
277 fc(ix^d,iw,idims)=dxinv(idims)*fc(ix^d,iw,idims)
279 {
do ix^db=ixomin^db,ixomax^db\}
280 hx^d=ix^d-kr(idims,^d)\
281 wnew(ix^d,iw)=wnew(ix^d,iw)+fc(ix^d,iw,idims)-fc(hx^d,iw,idims)
285 if(method==fs_tvdmu) &
286 call tvdlimit2(method,qdt,ixi^l,ixc^l,ixo^l,idims,wlc,wrc,wnew,x,fc,dxs)
290 inv_volume(ixo^s) = 1.d0/block%dvolume(ixo^s)
291 if(local_timestep)
then
293 hxomin^d=ixomin^d-kr(idims,^d)\
295 {
do ix^db=hxomin^db,ixomax^db\}
296 fc(ix^d,iw,idims)=-block%dt(ix^d)*dtfactor*fc(ix^d,iw,idims)*block%surfaceC(ix^d,idims)
298 {
do ix^db=ixomin^db,ixomax^db\}
299 hx^d=ix^d-kr(idims,^d)\
300 wnew(ix^d,iw)=wnew(ix^d,iw)+(fc(ix^d,iw,idims)-fc(hx^d,iw,idims))*inv_volume(ix^d)
304 if (method==fs_tvdmu) &
305 call tvdlimit2(method,qdt,ixi^l,ixc^l,ixo^l,idims,wlc,wrc,wnew,x,fc,dxs)
309 hxomin^d=ixomin^d-kr(idims,^d)\
311 {
do ix^db=hxomin^db,ixomax^db\}
312 fc(ix^d,iw,idims)=-qdt*fc(ix^d,iw,idims)*block%surfaceC(ix^d,idims)
314 {
do ix^db=ixomin^db,ixomax^db\}
315 hx^d=ix^d-kr(idims,^d)\
316 wnew(ix^d,iw)=wnew(ix^d,iw)+(fc(ix^d,iw,idims)-fc(hx^d,iw,idims))*inv_volume(ix^d)
323 if (.not.slab.and.idimsmin==1) &
324 call phys_add_source_geom(qdt,dtfactor,ixi^l,ixo^l,wct,wprim,wnew,x)
326 if(stagger_grid)
call phys_face_to_center(ixo^l,snew)
329 if(fix_small_values)
then
330 call phys_handle_small_values(.false.,wnew,x,ixi^l,ixo^l,
'multi-D finite_volume')
334 if(stagger_grid) snew%ws=sct%ws
338 call addsource2(qdt*dble(idimsmax-idimsmin+1)/dble(ndim),&
339 dtfactor*dble(idimsmax-idimsmin+1)/dble(ndim),&
340 ixi^l,ixo^l,1,nw,qtc,wct,wprim,qt,wnew,x,.false.,active)
347 fc(ixc^s,iw,idims)=half*(flc(ixc^s,iw)+frc(ixc^s,iw))
351 subroutine get_riemann_flux_tvdlf(iws,iwe)
352 integer,
intent(in) :: iws,iwe
355 double precision :: fac(ixC^S),phi
357 fac(ixc^s) = -0.5d0*tvdlfeps*cmaxc(ixc^s,ii)
359 fc(ixc^s,iw,idims)=0.5d0*(flc(ixc^s, iw)+frc(ixc^s, iw))
361 if(flux_type(idims, iw) /= flux_no_dissipation)
then
362 if(flux_adaptive_diffusion)
then
363 {
do ix^db=ixcmin^db,ixcmax^db\}
364 jx^d=ix^d+kr(idims,^d)\
367 if(((wrc(ix^d,iw)-wlc(ix^d,iw))*(sct%w(jx^d,iw)-sct%w(ix^d,iw))) .gt. 1.e-18)
then
368 phi=min((wrc(ix^d,iw)-wlc(ix^d,iw))**2/((sct%w(jx^d,iw)-sct%w(ix^d,iw))**2+1.e-18),one)
372 fc(ix^d,iw,idims)=fc(ix^d,iw,idims)+fac(ix^d)*(wrc(ix^d,iw)-wlc(ix^d,iw))*phi
375 {
do ix^db=ixcmin^db,ixcmax^db\}
376 fc(ix^d,iw,idims)=fc(ix^d,iw,idims)+fac(ix^d)*(wrc(ix^d,iw)-wlc(ix^d,iw))
382 end subroutine get_riemann_flux_tvdlf
384 subroutine get_riemann_flux_hll(iws,iwe)
385 integer,
intent(in) :: iws,iwe
387 double precision :: phi
390 if(flux_type(idims, iw) == flux_tvdlf)
then
392 if(stagger_grid) cycle
393 fc(ixc^s,iw,idims) = -tvdlfeps*half*max(cmaxc(ixc^s,ii),dabs(cminc(ixc^s,ii))) * &
394 (wrc(ixc^s,iw)-wlc(ixc^s,iw))
396 if(flux_adaptive_diffusion)
then
397 {
do ix^db=ixcmin^db,ixcmax^db\}
398 if(cminc(ix^d,ii) >= zero)
then
399 fc(ix^d,iw,idims)=flc(ix^d,iw)
400 else if(cmaxc(ix^d,ii) <= zero)
then
401 fc(ix^d,iw,idims)=frc(ix^d,iw)
404 phi=max(abs(cmaxc(ix^d,ii)),abs(cminc(ix^d,ii)))/(cmaxc(ix^d,ii)-cminc(ix^d,ii))
405 fc(ix^d,iw,idims)=(cmaxc(ix^d,ii)*flc(ix^d, iw)-cminc(ix^d,ii)*frc(ix^d,iw)&
406 +phi*cminc(ix^d,ii)*cmaxc(ix^d,ii)*(wrc(ix^d,iw)-wlc(ix^d,iw)))&
407 /(cmaxc(ix^d,ii)-cminc(ix^d,ii))
411 {
do ix^db=ixcmin^db,ixcmax^db\}
412 if(cminc(ix^d,ii) >= zero)
then
413 fc(ix^d,iw,idims)=flc(ix^d,iw)
414 else if(cmaxc(ix^d,ii) <= zero)
then
415 fc(ix^d,iw,idims)=frc(ix^d,iw)
417 fc(ix^d,iw,idims)=(cmaxc(ix^d,ii)*flc(ix^d, iw)-cminc(ix^d,ii)*frc(ix^d,iw)&
418 +cminc(ix^d,ii)*cmaxc(ix^d,ii)*(wrc(ix^d,iw)-wlc(ix^d,iw)))&
419 /(cmaxc(ix^d,ii)-cminc(ix^d,ii))
425 end subroutine get_riemann_flux_hll
427 subroutine get_riemann_flux_hllc(iws,iwe)
428 integer,
intent(in) :: iws, iwe
429 double precision,
dimension(ixI^S,1:nwflux) :: whll, Fhll, fCD
430 double precision,
dimension(ixI^S) :: lambdaCD
432 integer :: rho_, p_, e_, mom(1:ndir)
435 if (
allocated(iw_mom)) mom(:) = iw_mom(:)
438 if(
associated(phys_hllc_init_species))
then
439 call phys_hllc_init_species(ii, rho_, mom(:), e_)
445 where(cminc(ixc^s,1) >= zero)
447 elsewhere(cmaxc(ixc^s,1) <= zero)
451 if(method==fs_hllcd) &
452 call phys_diffuse_hllcd(ixi^l,ixc^l,idims,wlc,wrc,flc,frc,patchf)
455 if(any(patchf(ixc^s)==1)) &
456 call phys_get_lcd(wlc,wrc,flc,frc,cminc(ixi^s,ii),cmaxc(ixi^s,ii),idims,ixi^l,ixc^l, &
457 whll,fhll,lambdacd,patchf)
460 if(any(abs(patchf(ixc^s))== 1))
then
462 call phys_get_wcd(wlc,wrc,whll,frc,flc,fhll,patchf,lambdacd,&
463 cminc(ixi^s,ii),cmaxc(ixi^s,ii),ixi^l,ixc^l,idims,fcd)
467 if (flux_type(idims, iw) == flux_tvdlf)
then
468 flc(ixc^s,iw)=-tvdlfeps*half*max(cmaxc(ixc^s,ii),abs(cminc(ixc^s,ii))) * &
469 (wrc(ixc^s,iw) - wlc(ixc^s,iw))
471 where(patchf(ixc^s)==-2)
472 flc(ixc^s,iw)=flc(ixc^s,iw)
473 elsewhere(abs(patchf(ixc^s))==1)
474 flc(ixc^s,iw)=fcd(ixc^s,iw)
475 elsewhere(patchf(ixc^s)==2)
476 flc(ixc^s,iw)=frc(ixc^s,iw)
477 elsewhere(patchf(ixc^s)==3)
479 flc(ixc^s,iw)=fhll(ixc^s,iw)
480 elsewhere(patchf(ixc^s)==4)
482 flc(ixc^s,iw) = half*((flc(ixc^s,iw)+frc(ixc^s,iw)) &
483 -tvdlfeps * max(cmaxc(ixc^s,ii), dabs(cminc(ixc^s,ii))) * &
484 (wrc(ixc^s,iw)-wlc(ixc^s,iw)))
488 fc(ixc^s,iw,idims)=flc(ixc^s,iw)
491 end subroutine get_riemann_flux_hllc
494 subroutine get_riemann_flux_hlld(iws,iwe)
495 integer,
intent(in) :: iws, iwe
496 double precision,
dimension(ixI^S,1:nwflux) :: w1R,w1L,f1R,f1L,f2R,f2L
497 double precision,
dimension(ixI^S,1:nwflux) :: w2R,w2L
498 double precision,
dimension(ixI^S) :: sm,s1R,s1L,suR,suL,Bx
499 double precision,
dimension(ixI^S) :: pts,ptR,ptL,signBx,r1L,r1R,tmp
501 double precision,
dimension(ixI^S,ndir) :: vRC, vLC
503 double precision,
dimension(ixI^S,ndir) :: BR, BL
504 integer :: ip1,ip2,ip3,idir,ix^D
505 integer :: rho_, p_, e_, mom(1:ndir), mag(1:ndir)
507 associate(sr=>cmaxc,sl=>cminc)
525 vrc(ixc^s,:)=wrp(ixc^s,mom(:))
526 vlc(ixc^s,:)=wlp(ixc^s,mom(:))
528 br(ixc^s,:)=wrc(ixc^s,mag(:))+block%B0(ixc^s,:,ip1)
529 bl(ixc^s,:)=wlc(ixc^s,mag(:))+block%B0(ixc^s,:,ip1)
531 br(ixc^s,:)=wrc(ixc^s,mag(:))
532 bl(ixc^s,:)=wlc(ixc^s,mag(:))
534 if(stagger_grid)
then
535 bx(ixc^s)=block%ws(ixc^s,ip1)
539 bx(ixc^s)=(sr(ixc^s,ii)*br(ixc^s,ip1)-sl(ixc^s,ii)*bl(ixc^s,ip1))/(sr(ixc^s,ii)-sl(ixc^s,ii))
541 ptr(ixc^s)=wrp(ixc^s,p_)+0.5d0*sum(br(ixc^s,:)**2,dim=ndim+1)
542 ptl(ixc^s)=wlp(ixc^s,p_)+0.5d0*sum(bl(ixc^s,:)**2,dim=ndim+1)
543 if(iw_equi_rho>0)
then
544 sur(ixc^s) = wrc(ixc^s,rho_)+ block%equi_vars(ixc^s,iw_equi_rho,ip1)
546 sur(ixc^s) = wrc(ixc^s,rho_)
548 sur(ixc^s)=(sr(ixc^s,ii)-vrc(ixc^s,ip1))*sur(ixc^s)
549 if(iw_equi_rho>0)
then
550 sul(ixc^s) = wlc(ixc^s,rho_)+ block%equi_vars(ixc^s,iw_equi_rho,ip1)
552 sul(ixc^s) = wlc(ixc^s,rho_)
554 sul(ixc^s)=(sl(ixc^s,ii)-vlc(ixc^s,ip1))*sul(ixc^s)
556 sm(ixc^s)=(sur(ixc^s)*vrc(ixc^s,ip1)-sul(ixc^s)*vlc(ixc^s,ip1)-&
557 ptr(ixc^s)+ptl(ixc^s))/(sur(ixc^s)-sul(ixc^s))
559 w1r(ixc^s,mom(ip1))=sm(ixc^s)
560 w1l(ixc^s,mom(ip1))=sm(ixc^s)
561 w2r(ixc^s,mom(ip1))=sm(ixc^s)
562 w2l(ixc^s,mom(ip1))=sm(ixc^s)
564 w1r(ixc^s,mag(ip1))=bx(ixc^s)
565 w1l(ixc^s,mag(ip1))=bx(ixc^s)
567 ptr(ixc^s)=wrp(ixc^s,p_)+0.5d0*sum(wrc(ixc^s,mag(:))**2,dim=ndim+1)
568 ptl(ixc^s)=wlp(ixc^s,p_)+0.5d0*sum(wlc(ixc^s,mag(:))**2,dim=ndim+1)
572 w1r(ixc^s,rho_)=sur(ixc^s)/(sr(ixc^s,ii)-sm(ixc^s))
573 w1l(ixc^s,rho_)=sul(ixc^s)/(sl(ixc^s,ii)-sm(ixc^s))
577 r1r(ixc^s)=sur(ixc^s)*(sr(ixc^s,ii)-sm(ixc^s))-bx(ixc^s)**2
578 where(abs(r1r(ixc^s))>smalldouble)
579 r1r(ixc^s)=1.d0/r1r(ixc^s)
583 r1l(ixc^s)=sul(ixc^s)*(sl(ixc^s,ii)-sm(ixc^s))-bx(ixc^s)**2
584 where(abs(r1l(ixc^s))>smalldouble)
585 r1l(ixc^s)=1.d0/r1l(ixc^s)
590 w1r(ixc^s,mom(ip2))=vrc(ixc^s,ip2)-bx(ixc^s)*br(ixc^s,ip2)*&
591 (sm(ixc^s)-vrc(ixc^s,ip1))*r1r(ixc^s)
592 w1l(ixc^s,mom(ip2))=vlc(ixc^s,ip2)-bx(ixc^s)*bl(ixc^s,ip2)*&
593 (sm(ixc^s)-vlc(ixc^s,ip1))*r1l(ixc^s)
595 w1r(ixc^s,mag(ip2))=(sur(ixc^s)*(sr(ixc^s,ii)-vrc(ixc^s,ip1))-bx(ixc^s)**2)*r1r(ixc^s)
596 w1l(ixc^s,mag(ip2))=(sul(ixc^s)*(sl(ixc^s,ii)-vlc(ixc^s,ip1))-bx(ixc^s)**2)*r1l(ixc^s)
601 w1r(ixc^s,mom(ip3))=vrc(ixc^s,ip3)-bx(ixc^s)*br(ixc^s,ip3)*&
602 (sm(ixc^s)-vrc(ixc^s,ip1))*r1r(ixc^s)
603 w1l(ixc^s,mom(ip3))=vlc(ixc^s,ip3)-bx(ixc^s)*bl(ixc^s,ip3)*&
604 (sm(ixc^s)-vlc(ixc^s,ip1))*r1l(ixc^s)
606 w1r(ixc^s,mag(ip3))=br(ixc^s,ip3)*w1r(ixc^s,mag(ip2))
607 w1l(ixc^s,mag(ip3))=bl(ixc^s,ip3)*w1l(ixc^s,mag(ip2))
610 w1r(ixc^s,mag(ip2))=br(ixc^s,ip2)*w1r(ixc^s,mag(ip2))
611 w1l(ixc^s,mag(ip2))=bl(ixc^s,ip2)*w1l(ixc^s,mag(ip2))
614 w1r(ixc^s,mag(:))=w1r(ixc^s,mag(:))-block%B0(ixc^s,:,ip1)
615 w1l(ixc^s,mag(:))=w1l(ixc^s,mag(:))-block%B0(ixc^s,:,ip1)
620 w1r(ixc^s,p_)=sur(ixc^s)*(sm(ixc^s)-vrc(ixc^s,ip1))+ptr(ixc^s)
622 w1l(ixc^s,p_)=w1r(ixc^s,p_)
625 w1r(ixc^s,p_)=w1r(ixc^s,p_)+sum(block%B0(ixc^s,:,ip1)*(wrc(ixc^s,mag(:))-w1r(ixc^s,mag(:))),dim=ndim+1)
626 w1l(ixc^s,p_)=w1l(ixc^s,p_)+sum(block%B0(ixc^s,:,ip1)*(wlc(ixc^s,mag(:))-w1l(ixc^s,mag(:))),dim=ndim+1)
629 w1r(ixc^s,e_)=((sr(ixc^s,ii)-vrc(ixc^s,ip1))*wrc(ixc^s,e_)-ptr(ixc^s)*vrc(ixc^s,ip1)+&
630 w1r(ixc^s,p_)*sm(ixc^s)+bx(ixc^s)*(sum(vrc(ixc^s,:)*wrc(ixc^s,mag(:)),dim=ndim+1)-&
631 sum(w1r(ixc^s,mom(:))*w1r(ixc^s,mag(:)),dim=ndim+1)))/(sr(ixc^s,ii)-sm(ixc^s))
632 w1l(ixc^s,e_)=((sl(ixc^s,ii)-vlc(ixc^s,ip1))*wlc(ixc^s,e_)-ptl(ixc^s)*vlc(ixc^s,ip1)+&
633 w1l(ixc^s,p_)*sm(ixc^s)+bx(ixc^s)*(sum(vlc(ixc^s,:)*wlc(ixc^s,mag(:)),dim=ndim+1)-&
634 sum(w1l(ixc^s,mom(:))*w1l(ixc^s,mag(:)),dim=ndim+1)))/(sl(ixc^s,ii)-sm(ixc^s))
637 w1r(ixc^s,e_)=w1r(ixc^s,e_)+(sum(w1r(ixc^s,mag(:))*block%B0(ixc^s,:,ip1),dim=ndim+1)*sm(ixc^s)-&
638 sum(wrc(ixc^s,mag(:))*block%B0(ixc^s,:,ip1),dim=ndim+1)*vrc(ixc^s,ip1))/(sr(ixc^s,ii)-sm(ixc^s))
639 w1l(ixc^s,e_)=w1l(ixc^s,e_)+(sum(w1l(ixc^s,mag(:))*block%B0(ixc^s,:,ip1),dim=ndim+1)*sm(ixc^s)-&
640 sum(wlc(ixc^s,mag(:))*block%B0(ixc^s,:,ip1),dim=ndim+1)*vlc(ixc^s,ip1))/(sl(ixc^s,ii)-sm(ixc^s))
643 #if !defined(E_RM_W0) || E_RM_W0 == 1
644 w1r(ixc^s,e_)= w1r(ixc^s,e_) + 1d0/(phys_gamma - 1) * block%equi_vars(ixc^s,iw_equi_p,ip1) * &
645 (sm(ixc^s)-vrc(ixc^s,ip1))/(sr(ixc^s,ii)-sm(ixc^s))
646 w1l(ixc^s,e_)= w1l(ixc^s,e_) + 1d0/(phys_gamma - 1) * block%equi_vars(ixc^s,iw_equi_p,ip1) * &
647 (sm(ixc^s)-vlc(ixc^s,ip1))/(sl(ixc^s,ii)-sm(ixc^s))
649 w1r(ixc^s,e_)= w1r(ixc^s,e_) + phys_gamma /(phys_gamma - 1) * block%equi_vars(ixc^s,iw_equi_p,ip1) * &
650 (sm(ixc^s)-vrc(ixc^s,ip1))/(sr(ixc^s,ii)-sm(ixc^s))
651 w1l(ixc^s,e_)= w1l(ixc^s,e_) + phys_gamma /(phys_gamma - 1) * block%equi_vars(ixc^s,iw_equi_p,ip1) * &
652 (sm(ixc^s)-vlc(ixc^s,ip1))/(sl(ixc^s,ii)-sm(ixc^s))
658 w2r(ixc^s,rho_)=w1r(ixc^s,rho_)
659 w2l(ixc^s,rho_)=w1l(ixc^s,rho_)
660 w2r(ixc^s,mag(ip1))=w1r(ixc^s,mag(ip1))
661 w2l(ixc^s,mag(ip1))=w1l(ixc^s,mag(ip1))
663 r1r(ixc^s)=sqrt(w1r(ixc^s,rho_))
664 r1l(ixc^s)=sqrt(w1l(ixc^s,rho_))
665 tmp(ixc^s)=1.d0/(r1r(ixc^s)+r1l(ixc^s))
666 signbx(ixc^s)=sign(1.d0,bx(ixc^s))
668 s1r(ixc^s)=sm(ixc^s)+abs(bx(ixc^s))/r1r(ixc^s)
669 s1l(ixc^s)=sm(ixc^s)-abs(bx(ixc^s))/r1l(ixc^s)
671 w2r(ixc^s,mom(ip2))=(r1l(ixc^s)*w1l(ixc^s,mom(ip2))+r1r(ixc^s)*w1r(ixc^s,mom(ip2))+&
672 (w1r(ixc^s,mag(ip2))-w1l(ixc^s,mag(ip2)))*signbx(ixc^s))*tmp(ixc^s)
673 w2l(ixc^s,mom(ip2))=w2r(ixc^s,mom(ip2))
675 w2r(ixc^s,mag(ip2))=(r1l(ixc^s)*w1r(ixc^s,mag(ip2))+r1r(ixc^s)*w1l(ixc^s,mag(ip2))+&
676 r1l(ixc^s)*r1r(ixc^s)*(w1r(ixc^s,mom(ip2))-w1l(ixc^s,mom(ip2)))*signbx(ixc^s))*tmp(ixc^s)
677 w2l(ixc^s,mag(ip2))=w2r(ixc^s,mag(ip2))
680 w2r(ixc^s,mom(ip3))=(r1l(ixc^s)*w1l(ixc^s,mom(ip3))+r1r(ixc^s)*w1r(ixc^s,mom(ip3))+&
681 (w1r(ixc^s,mag(ip3))-w1l(ixc^s,mag(ip3)))*signbx(ixc^s))*tmp(ixc^s)
682 w2l(ixc^s,mom(ip3))=w2r(ixc^s,mom(ip3))
684 w2r(ixc^s,mag(ip3))=(r1l(ixc^s)*w1r(ixc^s,mag(ip3))+r1r(ixc^s)*w1l(ixc^s,mag(ip3))+&
685 r1l(ixc^s)*r1r(ixc^s)*(w1r(ixc^s,mom(ip3))-w1l(ixc^s,mom(ip3)))*signbx(ixc^s))*tmp(ixc^s)
686 w2l(ixc^s,mag(ip3))=w2r(ixc^s,mag(ip3))
690 w2r(ixc^s,e_)=w1r(ixc^s,e_)+r1r(ixc^s)*(sum(w1r(ixc^s,mom(:))*w1r(ixc^s,mag(:)),dim=ndim+1)-&
691 sum(w2r(ixc^s,mom(:))*w2r(ixc^s,mag(:)),dim=ndim+1))*signbx(ixc^s)
692 w2l(ixc^s,e_)=w1l(ixc^s,e_)-r1l(ixc^s)*(sum(w1l(ixc^s,mom(:))*w1l(ixc^s,mag(:)),dim=ndim+1)-&
693 sum(w2l(ixc^s,mom(:))*w2l(ixc^s,mag(:)),dim=ndim+1))*signbx(ixc^s)
698 w1r(ixc^s,mom(idir))=w1r(ixc^s,mom(idir))*w1r(ixc^s,rho_)
699 w1l(ixc^s,mom(idir))=w1l(ixc^s,mom(idir))*w1l(ixc^s,rho_)
700 w2r(ixc^s,mom(idir))=w2r(ixc^s,mom(idir))*w2r(ixc^s,rho_)
701 w2l(ixc^s,mom(idir))=w2l(ixc^s,mom(idir))*w2l(ixc^s,rho_)
703 if(iw_equi_rho>0)
then
704 w1r(ixc^s,rho_) = w1r(ixc^s,rho_) - block%equi_vars(ixc^s,iw_equi_rho,ip1)
705 w1l(ixc^s,rho_) = w1l(ixc^s,rho_) - block%equi_vars(ixc^s,iw_equi_rho,ip1)
706 w2r(ixc^s,rho_) = w2r(ixc^s,rho_) - block%equi_vars(ixc^s,iw_equi_rho,ip1)
707 w2l(ixc^s,rho_) = w2l(ixc^s,rho_) - block%equi_vars(ixc^s,iw_equi_rho,ip1)
711 if(flux_type(idims, iw) == flux_special)
then
713 f1l(ixc^s,iw)=flc(ixc^s,iw)
714 f1r(ixc^s,iw)=f1l(ixc^s,iw)
715 f2l(ixc^s,iw)=f1l(ixc^s,iw)
716 f2r(ixc^s,iw)=f1l(ixc^s,iw)
717 else if(flux_type(idims, iw) == flux_hll)
then
719 f1l(ixc^s,iw)=(sr(ixc^s,ii)*flc(ixc^s, iw)-sl(ixc^s,ii)*frc(ixc^s, iw) &
720 +sr(ixc^s,ii)*sl(ixc^s,ii)*(wrc(ixc^s,iw)-wlc(ixc^s,iw)))/(sr(ixc^s,ii)-sl(ixc^s,ii))
721 f1r(ixc^s,iw)=f1l(ixc^s,iw)
722 f2l(ixc^s,iw)=f1l(ixc^s,iw)
723 f2r(ixc^s,iw)=f1l(ixc^s,iw)
726 f1l(ixc^s,iw)=flc(ixc^s,iw)+sl(ixc^s,ii)*(w1l(ixc^s,iw)-wlc(ixc^s,iw))
727 f1r(ixc^s,iw)=frc(ixc^s,iw)+sr(ixc^s,ii)*(w1r(ixc^s,iw)-wrc(ixc^s,iw))
728 f2l(ixc^s,iw)=f1l(ixc^s,iw)+s1l(ixc^s)*(w2l(ixc^s,iw)-w1l(ixc^s,iw))
729 f2r(ixc^s,iw)=f1r(ixc^s,iw)+s1r(ixc^s)*(w2r(ixc^s,iw)-w1r(ixc^s,iw))
734 {
do ix^db=ixcmin^db,ixcmax^db\}
735 if(sl(ix^d,ii)>0.d0)
then
736 fc(ix^d,iws:iwe,ip1)=flc(ix^d,iws:iwe)
737 else if(s1l(ix^d)>=0.d0)
then
738 fc(ix^d,iws:iwe,ip1)=f1l(ix^d,iws:iwe)
739 else if(sm(ix^d)>=0.d0)
then
740 fc(ix^d,iws:iwe,ip1)=f2l(ix^d,iws:iwe)
741 else if(s1r(ix^d)>=0.d0)
then
742 fc(ix^d,iws:iwe,ip1)=f2r(ix^d,iws:iwe)
743 else if(sr(ix^d,ii)>=0.d0)
then
744 fc(ix^d,iws:iwe,ip1)=f1r(ix^d,iws:iwe)
745 else if(sr(ix^d,ii)<0.d0)
then
746 fc(ix^d,iws:iwe,ip1)=frc(ix^d,iws:iwe)
751 end subroutine get_riemann_flux_hlld
755 subroutine get_riemann_flux_hlld_mag2()
760 double precision,
dimension(ixI^S,1:nwflux) :: w1R,w1L,f1R,f1L,f2R,f2L
761 double precision,
dimension(ixI^S,1:nwflux) :: w2R,w2L
762 double precision,
dimension(ixI^S) :: sm,s1R,s1L,suR,suL,Bx
763 double precision,
dimension(ixI^S) :: pts,ptR,ptL,signBx,r1L,r1R,tmp
765 double precision,
dimension(ixI^S,ndir) :: vRC, vLC
767 double precision,
dimension(ixI^S,ndir) :: BR, BL
768 integer :: ip1,ip2,ip3,idir,ix^D
769 double precision :: phiPres, thetaSM, du, dv, dw
770 integer :: ixV^L, ixVb^L, ixVc^L, ixVd^L, ixVe^L, ixVf^L
771 integer :: rho_, p_, e_, mom(1:ndir), mag(1:ndir)
772 double precision,
parameter :: aParam = 4d0
780 associate(sr=>cmaxc,sl=>cminc)
787 vrc(ixc^s,:)=wrp(ixc^s,mom(:))
788 vlc(ixc^s,:)=wlp(ixc^s,mom(:))
794 phipres = min(1d0, maxval(max(s1l(ixo^s),s1r(ixo^s)))/maxval(cmaxc(ixo^s,1)))
795 phipres = phipres*(2d0 - phipres)
802 du = minval(wprim(ixv^s,mom(1))-wprim(ixo^s,mom(1)))
838 dv = minval(min(wprim(ixo^s,mom(2))-wprim(ixv^s,mom(2)),&
839 wprim(ixvb^s,mom(2))-wprim(ixo^s,mom(2)),&
840 wprim(ixvc^s,mom(2))-wprim(ixvd^s,mom(2)),&
841 wprim(ixve^s,mom(2))-wprim(ixvc^s,mom(2))&
875 dw = minval(min(wprim(ixo^s,mom(3))-wprim(ixv^s,mom(3)),&
876 wprim(ixvb^s,mom(3))-wprim(ixo^s,mom(3)),&
877 wprim(ixvc^s,mom(3))-wprim(ixvd^s,mom(3)),&
878 wprim(ixve^s,mom(3))-wprim(ixvc^s,mom(3))&
881 thetasm = maxval(cmaxc(ixo^s,1))
883 thetasm = (min(1d0, (thetasm-du)/(thetasm-min(dv,dw))))**aparam
887 br(ixc^s,:)=wrc(ixc^s,mag(:))+block%B0(ixc^s,:,ip1)
888 bl(ixc^s,:)=wlc(ixc^s,mag(:))+block%B0(ixc^s,:,ip1)
890 br(ixc^s,:)=wrc(ixc^s,mag(:))
891 bl(ixc^s,:)=wlc(ixc^s,mag(:))
896 ptr(ixc^s)=wrp(ixc^s,p_)+0.5d0*sum(br(ixc^s,:)**2,dim=ndim+1)
897 ptl(ixc^s)=wlp(ixc^s,p_)+0.5d0*sum(bl(ixc^s,:)**2,dim=ndim+1)
898 sur(ixc^s)=(sr(ixc^s,
index_v_mag)-vrc(ixc^s,ip1))*wrc(ixc^s,rho_)
899 sul(ixc^s)=(sl(ixc^s,
index_v_mag)-vlc(ixc^s,ip1))*wlc(ixc^s,rho_)
901 sm(ixc^s)=(sur(ixc^s)*vrc(ixc^s,ip1)-sul(ixc^s)*vlc(ixc^s,ip1)-&
902 thetasm*(ptr(ixc^s)-ptl(ixc^s)) )/(sur(ixc^s)-sul(ixc^s))
904 w1r(ixc^s,mom(ip1))=sm(ixc^s)
905 w1l(ixc^s,mom(ip1))=sm(ixc^s)
906 w2r(ixc^s,mom(ip1))=sm(ixc^s)
907 w2l(ixc^s,mom(ip1))=sm(ixc^s)
909 w1r(ixc^s,mag(ip1))=bx(ixc^s)
910 w1l(ixc^s,mag(ip1))=bx(ixc^s)
912 ptr(ixc^s)=wrp(ixc^s,p_)+0.5d0*sum(wrc(ixc^s,mag(:))**2,dim=ndim+1)
913 ptl(ixc^s)=wlp(ixc^s,p_)+0.5d0*sum(wlc(ixc^s,mag(:))**2,dim=ndim+1)
917 w1r(ixc^s,rho_)=sur(ixc^s)/(sr(ixc^s,
index_v_mag)-sm(ixc^s))
918 w1l(ixc^s,rho_)=sul(ixc^s)/(sl(ixc^s,
index_v_mag)-sm(ixc^s))
922 r1r(ixc^s)=sur(ixc^s)*(sr(ixc^s,
index_v_mag)-sm(ixc^s))-bx(ixc^s)**2
923 where(r1r(ixc^s)/=0.d0)
924 r1r(ixc^s)=1.d0/r1r(ixc^s)
926 r1l(ixc^s)=sul(ixc^s)*(sl(ixc^s,
index_v_mag)-sm(ixc^s))-bx(ixc^s)**2
927 where(r1l(ixc^s)/=0.d0)
928 r1l(ixc^s)=1.d0/r1l(ixc^s)
931 w1r(ixc^s,mom(ip2))=vrc(ixc^s,ip2)-bx(ixc^s)*br(ixc^s,ip2)*&
932 (sm(ixc^s)-vrc(ixc^s,ip1))*r1r(ixc^s)
933 w1l(ixc^s,mom(ip2))=vlc(ixc^s,ip2)-bx(ixc^s)*bl(ixc^s,ip2)*&
934 (sm(ixc^s)-vlc(ixc^s,ip1))*r1l(ixc^s)
936 w1r(ixc^s,mag(ip2))=(sur(ixc^s)*(sr(ixc^s,
index_v_mag)-vrc(ixc^s,ip1))-bx(ixc^s)**2)*r1r(ixc^s)
937 w1l(ixc^s,mag(ip2))=(sul(ixc^s)*(sl(ixc^s,
index_v_mag)-vlc(ixc^s,ip1))-bx(ixc^s)**2)*r1l(ixc^s)
942 w1r(ixc^s,mom(ip3))=vrc(ixc^s,ip3)-bx(ixc^s)*br(ixc^s,ip3)*&
943 (sm(ixc^s)-vrc(ixc^s,ip1))*r1r(ixc^s)
944 w1l(ixc^s,mom(ip3))=vlc(ixc^s,ip3)-bx(ixc^s)*bl(ixc^s,ip3)*&
945 (sm(ixc^s)-vlc(ixc^s,ip1))*r1l(ixc^s)
947 w1r(ixc^s,mag(ip3))=br(ixc^s,ip3)*w1r(ixc^s,mag(ip2))
948 w1l(ixc^s,mag(ip3))=bl(ixc^s,ip3)*w1l(ixc^s,mag(ip2))
951 w1r(ixc^s,mag(ip2))=br(ixc^s,ip2)*w1r(ixc^s,mag(ip2))
952 w1l(ixc^s,mag(ip2))=bl(ixc^s,ip2)*w1l(ixc^s,mag(ip2))
955 w1r(ixc^s,mag(:))=w1r(ixc^s,mag(:))-block%B0(ixc^s,:,ip1)
956 w1l(ixc^s,mag(:))=w1l(ixc^s,mag(:))-block%B0(ixc^s,:,ip1)
961 w1r(ixc^s,p_)=(sur(ixc^s)*ptl(ixc^s) - sul(ixc^s)*ptr(ixc^s) +&
962 phipres * sur(ixc^s)*sul(ixc^s)*(vrc(ixc^s,ip1)-vlc(ixc^s,ip1)))/&
963 (sur(ixc^s)-sul(ixc^s))
964 w1l(ixc^s,p_)=w1r(ixc^s,p_)
967 w1r(ixc^s,p_)=w1r(ixc^s,p_)+sum(block%B0(ixc^s,:,ip1)*(wrc(ixc^s,mag(:))-w1r(ixc^s,mag(:))),dim=ndim+1)
968 w1l(ixc^s,p_)=w1l(ixc^s,p_)+sum(block%B0(ixc^s,:,ip1)*(wlc(ixc^s,mag(:))-w1l(ixc^s,mag(:))),dim=ndim+1)
971 w1r(ixc^s,e_)=((sr(ixc^s,
index_v_mag)-vrc(ixc^s,ip1))*wrc(ixc^s,e_)-ptr(ixc^s)*vrc(ixc^s,ip1)+&
972 w1r(ixc^s,p_)*sm(ixc^s)+bx(ixc^s)*(sum(vrc(ixc^s,:)*wrc(ixc^s,mag(:)),dim=ndim+1)-&
973 sum(w1r(ixc^s,mom(:))*w1r(ixc^s,mag(:)),dim=ndim+1)))/(sr(ixc^s,
index_v_mag)-sm(ixc^s))
974 w1l(ixc^s,e_)=((sl(ixc^s,
index_v_mag)-vlc(ixc^s,ip1))*wlc(ixc^s,e_)-ptl(ixc^s)*vlc(ixc^s,ip1)+&
975 w1l(ixc^s,p_)*sm(ixc^s)+bx(ixc^s)*(sum(vlc(ixc^s,:)*wlc(ixc^s,mag(:)),dim=ndim+1)-&
976 sum(w1l(ixc^s,mom(:))*w1l(ixc^s,mag(:)),dim=ndim+1)))/(sl(ixc^s,
index_v_mag)-sm(ixc^s))
979 w1r(ixc^s,e_)=w1r(ixc^s,e_)+(sum(w1r(ixc^s,mag(:))*block%B0(ixc^s,:,ip1),dim=ndim+1)*sm(ixc^s)-&
980 sum(wrc(ixc^s,mag(:))*block%B0(ixc^s,:,ip1),dim=ndim+1)*vrc(ixc^s,ip1))/(sr(ixc^s,
index_v_mag)-sm(ixc^s))
981 w1l(ixc^s,e_)=w1l(ixc^s,e_)+(sum(w1l(ixc^s,mag(:))*block%B0(ixc^s,:,ip1),dim=ndim+1)*sm(ixc^s)-&
982 sum(wlc(ixc^s,mag(:))*block%B0(ixc^s,:,ip1),dim=ndim+1)*vlc(ixc^s,ip1))/(sl(ixc^s,
index_v_mag)-sm(ixc^s))
987 w2r(ixc^s,rho_)=w1r(ixc^s,rho_)
988 w2l(ixc^s,rho_)=w1l(ixc^s,rho_)
989 w2r(ixc^s,mag(ip1))=w1r(ixc^s,mag(ip1))
990 w2l(ixc^s,mag(ip1))=w1l(ixc^s,mag(ip1))
992 r1r(ixc^s)=sqrt(w1r(ixc^s,rho_))
993 r1l(ixc^s)=sqrt(w1l(ixc^s,rho_))
994 tmp(ixc^s)=1.d0/(r1r(ixc^s)+r1l(ixc^s))
995 signbx(ixc^s)=sign(1.d0,bx(ixc^s))
997 s1r(ixc^s)=sm(ixc^s)+abs(bx(ixc^s))/r1r(ixc^s)
998 s1l(ixc^s)=sm(ixc^s)-abs(bx(ixc^s))/r1l(ixc^s)
1000 w2r(ixc^s,mom(ip2))=(r1l(ixc^s)*w1l(ixc^s,mom(ip2))+r1r(ixc^s)*w1r(ixc^s,mom(ip2))+&
1001 (w1r(ixc^s,mag(ip2))-w1l(ixc^s,mag(ip2)))*signbx(ixc^s))*tmp(ixc^s)
1002 w2l(ixc^s,mom(ip2))=w2r(ixc^s,mom(ip2))
1004 w2r(ixc^s,mag(ip2))=(r1l(ixc^s)*w1r(ixc^s,mag(ip2))+r1r(ixc^s)*w1l(ixc^s,mag(ip2))+&
1005 r1l(ixc^s)*r1r(ixc^s)*(w1r(ixc^s,mom(ip2))-w1l(ixc^s,mom(ip2)))*signbx(ixc^s))*tmp(ixc^s)
1006 w2l(ixc^s,mag(ip2))=w2r(ixc^s,mag(ip2))
1009 w2r(ixc^s,mom(ip3))=(r1l(ixc^s)*w1l(ixc^s,mom(ip3))+r1r(ixc^s)*w1r(ixc^s,mom(ip3))+&
1010 (w1r(ixc^s,mag(ip3))-w1l(ixc^s,mag(ip3)))*signbx(ixc^s))*tmp(ixc^s)
1011 w2l(ixc^s,mom(ip3))=w2r(ixc^s,mom(ip3))
1013 w2r(ixc^s,mag(ip3))=(r1l(ixc^s)*w1r(ixc^s,mag(ip3))+r1r(ixc^s)*w1l(ixc^s,mag(ip3))+&
1014 r1l(ixc^s)*r1r(ixc^s)*(w1r(ixc^s,mom(ip3))-w1l(ixc^s,mom(ip3)))*signbx(ixc^s))*tmp(ixc^s)
1015 w2l(ixc^s,mag(ip3))=w2r(ixc^s,mag(ip3))
1019 w2r(ixc^s,e_)=w1r(ixc^s,e_)+r1r(ixc^s)*(sum(w1r(ixc^s,mom(:))*w1r(ixc^s,mag(:)),dim=ndim+1)-&
1020 sum(w2r(ixc^s,mom(:))*w2r(ixc^s,mag(:)),dim=ndim+1))*signbx(ixc^s)
1021 w2l(ixc^s,e_)=w1l(ixc^s,e_)-r1l(ixc^s)*(sum(w1l(ixc^s,mom(:))*w1l(ixc^s,mag(:)),dim=ndim+1)-&
1022 sum(w2l(ixc^s,mom(:))*w2l(ixc^s,mag(:)),dim=ndim+1))*signbx(ixc^s)
1027 w1r(ixc^s,mom(idir))=w1r(ixc^s,mom(idir))*w1r(ixc^s,rho_)
1028 w1l(ixc^s,mom(idir))=w1l(ixc^s,mom(idir))*w1l(ixc^s,rho_)
1029 w2r(ixc^s,mom(idir))=w2r(ixc^s,mom(idir))*w2r(ixc^s,rho_)
1030 w2l(ixc^s,mom(idir))=w2l(ixc^s,mom(idir))*w2l(ixc^s,rho_)
1039 f1l(ixc^s,iw)=flc(ixc^s,iw)
1040 f1r(ixc^s,iw)=f1l(ixc^s,iw)
1041 f2l(ixc^s,iw)=f1l(ixc^s,iw)
1042 f2r(ixc^s,iw)=f1l(ixc^s,iw)
1047 f1r(ixc^s,iw)=f1l(ixc^s,iw)
1048 f2l(ixc^s,iw)=f1l(ixc^s,iw)
1049 f2r(ixc^s,iw)=f1l(ixc^s,iw)
1051 f1l(ixc^s,iw)=flc(ixc^s,iw)+sl(ixc^s,
index_v_mag)*(w1l(ixc^s,iw)-wlc(ixc^s,iw))
1052 f1r(ixc^s,iw)=frc(ixc^s,iw)+sr(ixc^s,
index_v_mag)*(w1r(ixc^s,iw)-wrc(ixc^s,iw))
1053 f2l(ixc^s,iw)=f1l(ixc^s,iw)+s1l(ixc^s)*(w2l(ixc^s,iw)-w1l(ixc^s,iw))
1054 f2r(ixc^s,iw)=f1r(ixc^s,iw)+s1r(ixc^s)*(w2r(ixc^s,iw)-w1r(ixc^s,iw))
1059 {
do ix^db=ixcmin^db,ixcmax^db\}
1062 else if(s1l(ix^d)>=0.d0)
then
1064 else if(sm(ix^d)>=0.d0)
then
1066 else if(s1r(ix^d)>=0.d0)
then
1076 end subroutine get_riemann_flux_hlld_mag2
1082 integer,
intent(in) :: ixI^L, ixO^L
1083 double precision,
intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
1084 double precision,
intent(out):: csound(ixI^S)
1085 double precision :: cfast2(ixI^S), AvMinCs2(ixI^S), b2(ixI^S), kmax
1086 double precision :: inv_rho(ixO^S), gamma_A2(ixO^S)
1087 integer :: rho_, p_, e_, mom(1:ndir), mag(1:ndir)
1095 inv_rho=1.d0/w(ixo^s,rho_)
1100 b2(ixo^s) = sum((w(ixo^s, mag(:))+
block%B0(ixo^s,:,
b0i))**2, dim=ndim+1)
1102 b2(ixo^s) = sum(w(ixo^s, mag(:))**2, dim=ndim+1)
1107 avmincs2= w(ixo^s, mag(idims))+
block%B0(ixo^s,idims,
b0i)
1109 avmincs2= w(ixo^s, mag(idims))
1113 csound(ixo^s) = sum(w(ixo^s, mom(:))**2, dim=ndim+1)
1115 cfast2(ixo^s) = b2(ixo^s) * inv_rho+csound(ixo^s)
1116 avmincs2(ixo^s) = cfast2(ixo^s)**2-4.0d0*csound(ixo^s) &
1117 * avmincs2(ixo^s)**2 &
1120 where(avmincs2(ixo^s)<zero)
1121 avmincs2(ixo^s)=zero
1124 avmincs2(ixo^s)=sqrt(avmincs2(ixo^s))
1126 csound(ixo^s) = sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s)))
1134 subroutine reconstruct_lr(ixI^L,ixL^L,ixR^L,idims,w,wLC,wRC,wLp,wRp,x,dxdim)
1140 integer,
intent(in) :: ixi^
l, ixl^
l, ixr^
l, idims
1141 double precision,
intent(in) :: dxdim
1143 double precision,
dimension(ixI^S,1:nw) :: w
1145 double precision,
dimension(ixI^S,1:nw) :: wlc, wrc
1147 double precision,
dimension(ixI^S,1:nw) :: wlp, wrp
1148 double precision,
dimension(ixI^S,1:ndim) :: x
1150 integer :: jxr^
l, ixc^
l, jxc^
l, iw
1151 double precision :: ldw(ixi^s), rdw(ixi^s), dwc(ixi^s)
1152 double precision :: a2max
1156 call mp5limiter(ixi^
l,ixl^
l,idims,w,wlp,wrp)
1158 call weno3limiter(ixi^
l,ixl^
l,idims,dxdim,w,wlp,wrp,1)
1160 call weno3limiter(ixi^
l,ixl^
l,idims,dxdim,w,wlp,wrp,2)
1162 call weno5limiter(ixi^
l,ixl^
l,idims,dxdim,w,wlp,wrp,1)
1164 call weno5nmlimiter(ixi^
l,ixl^
l,idims,dxdim,w,wlp,wrp,1)
1166 call weno5limiter(ixi^
l,ixl^
l,idims,dxdim,w,wlp,wrp,2)
1168 call weno5nmlimiter(ixi^
l,ixl^
l,idims,dxdim,w,wlp,wrp,2)
1170 call weno5limiter(ixi^
l,ixl^
l,idims,dxdim,w,wlp,wrp,3)
1172 call weno5nmlimiter(ixi^
l,ixl^
l,idims,dxdim,w,wlp,wrp,3)
1174 call weno5cu6limiter(ixi^
l,ixl^
l,idims,w,wlp,wrp)
1176 call teno5adlimiter(ixi^
l,ixl^
l,idims,dxdim,w,wlp,wrp)
1178 call weno7limiter(ixi^
l,ixl^
l,idims,w,wlp,wrp,1)
1180 call weno7limiter(ixi^
l,ixl^
l,idims,w,wlp,wrp,2)
1182 call venklimiter(ixi^
l,ixl^
l,idims,dxdim,w,wlp,wrp)
1188 call ppmlimiter(ixi^
l,
ixm^
ll,idims,w,w,wlp,wrp)
1194 jxr^
l=ixr^
l+
kr(idims,^
d);
1195 ixcmax^
d=jxrmax^
d; ixcmin^
d=ixlmin^
d-
kr(idims,^
d);
1196 jxc^
l=ixc^
l+
kr(idims,^
d);
1199 w(ixcmin^
d:jxcmax^
d,iw)=dlog10(w(ixcmin^
d:jxcmax^
d,iw))
1200 wlp(ixl^s,iw)=dlog10(wlp(ixl^s,iw))
1201 wrp(ixr^s,iw)=dlog10(wrp(ixr^s,iw))
1204 dwc(ixc^s)=w(jxc^s,iw)-w(ixc^s,iw)
1220 call mpistop(
"idims is wrong in mod_limiter")
1226 wlp(ixl^s,iw)=wlp(ixl^s,iw)+half*ldw(ixl^s)
1227 wrp(ixr^s,iw)=wrp(ixr^s,iw)-half*rdw(jxr^s)
1230 w(ixcmin^
d:jxcmax^
d,iw)=10.0d0**w(ixcmin^
d:jxcmax^
d,iw)
1231 wlp(ixl^s,iw)=10.0d0**wlp(ixl^s,iw)
1232 wrp(ixr^s,iw)=10.0d0**wrp(ixr^s,iw)
1241 wlc(ixl^s,1:nwflux) = wlp(ixl^s,1:nwflux)
1242 wrc(ixr^s,1:nwflux) = wrp(ixr^s,1:nwflux)
1246 wlp(ixl^s,nwflux+1:nwflux+nwaux) = wlc(ixl^s,nwflux+1:nwflux+nwaux)
1247 wrp(ixr^s,nwflux+1:nwflux+nwaux) = wrc(ixr^s,nwflux+1:nwflux+nwaux)
subroutine get_hlld2_modif_c(w, x, ixIL, ixOL, csound)
Calculate fast magnetosonic wave speed.
subroutine get_riemann_flux_tvdmu()
subroutine, public mpistop(message)
Exit MPI-AMRVAC with an error message.
Module with finite volume methods for fluxes.
subroutine, public finite_volume(method, qdt, dtfactor, ixIL, ixOL, idimsLIM, qtC, sCT, qt, snew, fC, fE, dxs, x)
finite volume method
subroutine, public reconstruct_lr(ixIL, ixLL, ixRL, idims, w, wLC, wRC, wLp, wRp, x, dxdim)
Determine the upwinded wLC(ixL) and wRC(ixR) from w. the wCT is only used when PPM is exploited.
subroutine, public hancock(qdt, dtfactor, ixIL, ixOL, idimsLIM, qtC, sCT, qt, snew, dxs, x)
The non-conservative Hancock predictor for TVDLF.
This module contains definitions of global parameters and variables and some generic functions/subrou...
type(state), pointer block
Block pointer for using one block and its previous state.
logical h_correction
If true, do H-correction to fix the carbuncle problem at grid-aligned shocks.
integer, parameter fs_tvdlf
integer ixghi
Upper index of grid block arrays.
integer, parameter fs_hlld
integer, dimension(3, 3) kr
Kronecker delta tensor.
logical, dimension(:), allocatable loglimit
integer, parameter fs_hllcd
integer, parameter ndim
Number of spatial dimensions for grid variables.
logical stagger_grid
True for using stagger grid.
integer, parameter fs_tvdmu
integer b0i
background magnetic field location indicator
integer, dimension(:), allocatable, parameter d
logical local_timestep
each cell has its own timestep or not
logical need_global_a2max
global value for schmid scheme
integer ixm
the mesh range of a physical block without ghost cells
logical slab
Cartesian geometry or not.
integer, parameter fs_hll
flux schemes
logical b0field
split magnetic field as background B0 field
integer, dimension(:), allocatable type_limiter
Type of slope limiter used for reconstructing variables on cell edges.
integer nghostcells
Number of ghost cells surrounding a grid.
logical slab_uniform
uniform Cartesian geometry or not (stretched Cartesian)
integer, parameter fs_hllc
double precision, dimension(ndim) a2max_global
global largest a2 for schmid scheme
integer, parameter ixglo
Lower index of grid block arrays (always 1)
Module with slope/flux limiters.
integer, parameter limiter_weno5cu6
integer, parameter limiter_mpweno7
integer, parameter limiter_teno5ad
integer, parameter limiter_weno3
integer, parameter limiter_ppm
subroutine dwlimiter2(dwC, ixIL, ixCL, idims, typelim, ldw, rdw, a2max)
Limit the centered dwC differences within ixC for iw in direction idim. The limiter is chosen accordi...
integer, parameter limiter_wenozp5
integer, parameter limiter_weno5
integer, parameter limiter_wenoz5
integer, parameter limiter_wenoz5nm
integer, parameter limiter_weno7
integer, parameter limiter_wenozp5nm
integer, parameter limiter_mp5
integer, parameter limiter_venk
integer, parameter limiter_weno5nm
integer, parameter limiter_wenoyc3
This module defines the procedures of a physics module. It contains function pointers for the various...
procedure(sub_get_ct_velocity), pointer phys_get_ct_velocity
procedure(sub_convert), pointer phys_to_primitive
integer, parameter flux_tvdlf
Indicates the flux should be treated with tvdlf.
procedure(sub_small_values), pointer phys_handle_small_values
procedure(sub_get_flux), pointer phys_get_flux
procedure(sub_get_cbounds), pointer phys_get_cbounds
integer, parameter flux_hll
Indicates the flux should be treated with hll.
procedure(sub_add_source_geom), pointer phys_add_source_geom
integer, parameter flux_special
Indicates the flux should be specially treated.
integer, dimension(:, :), allocatable flux_type
Array per direction per variable, which can be used to specify that certain fluxes have to be treated...
integer phys_wider_stencil
To use wider stencils in flux calculations. A value of 1 will extend it by one cell in both direction...
procedure(sub_update_faces), pointer phys_update_faces
procedure(sub_convert), pointer phys_to_conserved
procedure(sub_modify_wlr), pointer phys_modify_wlr
procedure(sub_get_h_speed), pointer phys_get_h_speed
logical phys_energy
Solve energy equation or not.
Module for handling split source terms (split from the fluxes)
subroutine, public addsource2(qdt, dtfactor, ixIL, ixOL, iwLIM, qtC, wCT, wCTprim, qt, w, x, qsourcesplit, src_active)
Add source within ixO for iws: w=w+qdt*S[wCT].
Subroutines for TVD-MUSCL schemes.
subroutine, public tvdlimit2(method, qdt, ixIL, ixICL, ixOL, idims, wL, wR, wnew, x, fC, dxs)
Module with all the methods that users can customize in AMRVAC.
integer nw
Total number of variables.
integer number_species
number of species: each species has different characterictic speeds and should be used accordingly in...
integer, dimension(:), allocatable iw_mom
Indices of the momentum density.
integer, dimension(:), allocatable start_indices
the indices in 1:nwflux array are assumed consecutive for each species this array should be of size n...
integer, dimension(:), allocatable stop_indices
the indices in 1:nwflux array are assumed consecutive for each species this array should be of size n...
integer, dimension(:), allocatable, protected iw_mag
Indices of the magnetic field components.
integer iw_rho
Index of the (gas) density.
integer nwflux
Number of flux variables.
integer index_v_mag
index of the var whose velocity appears in the induction eq.
integer iw_e
Index of the energy density.