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
135 double precision,
dimension(ixI^S,1:nw) :: wprim, wcons
137 double precision,
dimension(ixI^S,1:nw) :: wlc, wrc
139 double precision,
dimension(ixI^S,1:nw) :: wlp, wrp
140 double precision,
dimension(ixI^S,1:nwflux) :: flc, frc
141 double precision,
dimension(ixI^S,1:number_species) :: cmaxc
142 double precision,
dimension(ixI^S,1:number_species) :: cminc
143 double precision,
dimension(ixI^S) :: hspeed
144 double precision,
dimension(ixO^S) :: inv_volume
145 double precision,
dimension(1:ndim) :: dxinv
146 integer,
dimension(ixI^S) :: patchf
147 integer :: idims, iw, ix^
l, hxo^
l, ixc^
l, jxc^
l,ixcr^
l, kxc^
l, kxr^
l, ii
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")
176 hxo^
l=ixo^
l-
kr(idims,^
d);
178 kxcmin^
d=iximin^
d; kxcmax^
d=iximax^
d-
kr(idims,^
d);
179 kxr^
l=kxc^
l+
kr(idims,^
d);
187 ixcmax^
d=ixomax^
d; ixcmin^
d=hxomin^
d;
194 wrp(kxc^s,1:
nw)=wprim(kxr^s,1:
nw)
195 wlp(kxc^s,1:
nw)=wprim(kxc^s,1:
nw)
201 call reconstruct_lr(ixi^
l,ixcr^
l,ixcr^
l,idims,wprim,wlc,wrc,wlp,wrp,x,dxs(idims))
219 call phys_get_cbounds(wlc,wrc,wlp,wrp,x,ixi^
l,ixc^
l,idims,hspeed,cmaxc,cminc)
248 call mpistop(
'unkown Riemann flux in finite volume')
257 hxo^
l=ixo^
l-
kr(idims,^
d);
262 fc(ixi^s,iw,idims)=-
block%dt(ixi^s)*dtfactor/dxs(idims)*fc(ixi^s,iw,idims)
266 fc(ixi^s,1:
nwflux,idims)=dxinv(idims)*fc(ixi^s,1:
nwflux,idims)
274 call tvdlimit2(method,qdt,ixi^
l,ixc^
l,ixo^
l,idims,wlc,wrc,wnew,x,fc,dxs)
278 inv_volume = 1.d0/
block%dvolume(ixo^s)
280 hxo^
l=ixo^
l-
kr(idims,^
d);
284 fc(ixi^s,iw,idims)=-
block%dt(ixi^s)*dtfactor*fc(ixi^s,iw,idims)*
block%surfaceC(ixi^s,idims)
285 wnew(ixo^s,iw)=wnew(ixo^s,iw) + (fc(ixo^s,iw,idims)-fc(hxo^s,iw,idims)) * &
290 fc(ixi^s,iw,idims)=-qdt*fc(ixi^s,iw,idims)*
block%surfaceC(ixi^s,idims)
291 wnew(ixo^s,iw)=wnew(ixo^s,iw) + (fc(ixo^s,iw,idims)-fc(hxo^s,iw,idims)) * &
297 call tvdlimit2(method,qdt,ixi^
l,ixc^
l,ixo^
l,idims,wlc,wrc,wnew,x,fc,dxs)
302 if (.not.
slab.and.idimsmin==1) &
318 dtfactor*dble(idimsmax-idimsmin+1)/dble(
ndim),&
319 ixi^
l,ixo^
l,1,
nw,qtc,wct,wprim,qt,wnew,x,.false.,active)
327 flc(ixc^s, iw)=half*(flc(ixc^s, iw)+frc(ixc^s, iw))
328 fc(ixc^s,iw,idims)=flc(ixc^s, iw)
332 subroutine get_riemann_flux_tvdlf(iws,iwe)
333 integer,
intent(in) :: iws,iwe
334 double precision :: fac(ixC^S),phi(ixC^S)
336 jxc^l=ixc^l+kr(idims,^d);
337 fac(ixc^s) = -0.5d0*tvdlfeps*cmaxc(ixc^s,ii)
341 flc(ixc^s, iw)=0.5d0*(flc(ixc^s, iw)+frc(ixc^s, iw))
343 if(flux_type(idims, iw) /= flux_no_dissipation)
then
344 if(flux_adaptive_diffusion)
then
348 where(((wrc(ixc^s,iw)-wlc(ixc^s,iw))*(wcons(jxc^s,iw)-wcons(ixc^s,iw))) .gt. 1.e-18)
349 phi(ixc^s)=min((wrc(ixc^s,iw)-wlc(ixc^s,iw))**2/((wcons(jxc^s,iw)-wcons(ixc^s,iw))**2+1.e-18),one)
354 flc(ixc^s,iw)=flc(ixc^s,iw) + fac(ixc^s)*(wrc(ixc^s,iw)-wlc(ixc^s,iw))*phi(ixc^s)
356 fc(ixc^s,iw,idims)=flc(ixc^s, iw)
358 end subroutine get_riemann_flux_tvdlf
360 subroutine get_riemann_flux_hll(iws,iwe)
361 integer,
intent(in) :: iws,iwe
363 double precision :: phi(ixI^S)
366 if(flux_type(idims, iw) == flux_tvdlf)
then
368 if(stagger_grid) cycle
369 fc(ixc^s,iw,idims) = -tvdlfeps*half*max(cmaxc(ixc^s,ii),dabs(cminc(ixc^s,ii))) * &
370 (wrc(ixc^s,iw)-wlc(ixc^s,iw))
372 {
do ix^db=ixcmin^db,ixcmax^db\}
373 if(cminc(ix^d,ii) >= zero)
then
374 fc(ix^d,iw,idims)=flc(ix^d,iw)
375 else if(cmaxc(ix^d,ii) <= zero)
then
376 fc(ix^d,iw,idims)=frc(ix^d,iw)
379 if(flux_adaptive_diffusion)
then
381 phi(ix^d)=max(abs(cmaxc(ix^d,ii)),abs(cminc(ix^d,ii)))/(cmaxc(ix^d,ii)-cminc(ix^d,ii))
382 fc(ix^d,iw,idims)=(cmaxc(ix^d,ii)*flc(ix^d, iw)-cminc(ix^d,ii)*frc(ix^d,iw)&
383 +phi(ix^d)*cminc(ix^d,ii)*cmaxc(ix^d,ii)*(wrc(ix^d,iw)-wlc(ix^d,iw)))&
384 /(cmaxc(ix^d,ii)-cminc(ix^d,ii))
386 fc(ix^d,iw,idims)=(cmaxc(ix^d,ii)*flc(ix^d, iw)-cminc(ix^d,ii)*frc(ix^d,iw)&
387 +cminc(ix^d,ii)*cmaxc(ix^d,ii)*(wrc(ix^d,iw)-wlc(ix^d,iw)))&
388 /(cmaxc(ix^d,ii)-cminc(ix^d,ii))
394 end subroutine get_riemann_flux_hll
396 subroutine get_riemann_flux_hllc(iws,iwe)
397 integer,
intent(in) :: iws, iwe
398 double precision,
dimension(ixI^S,1:nwflux) :: whll, Fhll, fCD
399 double precision,
dimension(ixI^S) :: lambdaCD
401 integer :: rho_, p_, e_, mom(1:ndir)
404 if (
allocated(iw_mom)) mom(:) = iw_mom(:)
407 if(
associated(phys_hllc_init_species))
then
408 call phys_hllc_init_species(ii, rho_, mom(:), e_)
414 where(cminc(ixc^s,1) >= zero)
416 elsewhere(cmaxc(ixc^s,1) <= zero)
420 if(method==fs_hllcd) &
421 call phys_diffuse_hllcd(ixi^l,ixc^l,idims,wlc,wrc,flc,frc,patchf)
424 if(any(patchf(ixc^s)==1)) &
425 call phys_get_lcd(wlc,wrc,flc,frc,cminc(ixi^s,ii),cmaxc(ixi^s,ii),idims,ixi^l,ixc^l, &
426 whll,fhll,lambdacd,patchf)
429 if(any(abs(patchf(ixc^s))== 1))
then
431 call phys_get_wcd(wlc,wrc,whll,frc,flc,fhll,patchf,lambdacd,&
432 cminc(ixi^s,ii),cmaxc(ixi^s,ii),ixi^l,ixc^l,idims,fcd)
436 if (flux_type(idims, iw) == flux_tvdlf)
then
437 flc(ixc^s,iw)=-tvdlfeps*half*max(cmaxc(ixc^s,ii),abs(cminc(ixc^s,ii))) * &
438 (wrc(ixc^s,iw) - wlc(ixc^s,iw))
440 where(patchf(ixc^s)==-2)
441 flc(ixc^s,iw)=flc(ixc^s,iw)
442 elsewhere(abs(patchf(ixc^s))==1)
443 flc(ixc^s,iw)=fcd(ixc^s,iw)
444 elsewhere(patchf(ixc^s)==2)
445 flc(ixc^s,iw)=frc(ixc^s,iw)
446 elsewhere(patchf(ixc^s)==3)
448 flc(ixc^s,iw)=fhll(ixc^s,iw)
449 elsewhere(patchf(ixc^s)==4)
451 flc(ixc^s,iw) = half*((flc(ixc^s,iw)+frc(ixc^s,iw)) &
452 -tvdlfeps * max(cmaxc(ixc^s,ii), dabs(cminc(ixc^s,ii))) * &
453 (wrc(ixc^s,iw)-wlc(ixc^s,iw)))
457 fc(ixc^s,iw,idims)=flc(ixc^s,iw)
460 end subroutine get_riemann_flux_hllc
463 subroutine get_riemann_flux_hlld(iws,iwe)
464 integer,
intent(in) :: iws, iwe
465 double precision,
dimension(ixI^S,1:nwflux) :: w1R,w1L,f1R,f1L,f2R,f2L
466 double precision,
dimension(ixI^S,1:nwflux) :: w2R,w2L
467 double precision,
dimension(ixI^S) :: sm,s1R,s1L,suR,suL,Bx
468 double precision,
dimension(ixI^S) :: pts,ptR,ptL,signBx,r1L,r1R,tmp
470 double precision,
dimension(ixI^S,ndir) :: vRC, vLC
472 double precision,
dimension(ixI^S,ndir) :: BR, BL
473 integer :: ip1,ip2,ip3,idir,ix^D
474 integer :: rho_, p_, e_, mom(1:ndir), mag(1:ndir)
476 associate(sr=>cmaxc,sl=>cminc)
494 vrc(ixc^s,:)=wrp(ixc^s,mom(:))
495 vlc(ixc^s,:)=wlp(ixc^s,mom(:))
497 br(ixc^s,:)=wrc(ixc^s,mag(:))+block%B0(ixc^s,:,ip1)
498 bl(ixc^s,:)=wlc(ixc^s,mag(:))+block%B0(ixc^s,:,ip1)
500 br(ixc^s,:)=wrc(ixc^s,mag(:))
501 bl(ixc^s,:)=wlc(ixc^s,mag(:))
503 if(stagger_grid)
then
504 bx(ixc^s)=block%ws(ixc^s,ip1)
508 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))
510 ptr(ixc^s)=wrp(ixc^s,p_)+0.5d0*sum(br(ixc^s,:)**2,dim=ndim+1)
511 ptl(ixc^s)=wlp(ixc^s,p_)+0.5d0*sum(bl(ixc^s,:)**2,dim=ndim+1)
512 if(iw_equi_rho>0)
then
513 sur(ixc^s) = wrc(ixc^s,rho_)+ block%equi_vars(ixc^s,iw_equi_rho,ip1)
515 sur(ixc^s) = wrc(ixc^s,rho_)
517 sur(ixc^s)=(sr(ixc^s,ii)-vrc(ixc^s,ip1))*sur(ixc^s)
518 if(iw_equi_rho>0)
then
519 sul(ixc^s) = wlc(ixc^s,rho_)+ block%equi_vars(ixc^s,iw_equi_rho,ip1)
521 sul(ixc^s) = wlc(ixc^s,rho_)
523 sul(ixc^s)=(sl(ixc^s,ii)-vlc(ixc^s,ip1))*sul(ixc^s)
525 sm(ixc^s)=(sur(ixc^s)*vrc(ixc^s,ip1)-sul(ixc^s)*vlc(ixc^s,ip1)-&
526 ptr(ixc^s)+ptl(ixc^s))/(sur(ixc^s)-sul(ixc^s))
528 w1r(ixc^s,mom(ip1))=sm(ixc^s)
529 w1l(ixc^s,mom(ip1))=sm(ixc^s)
530 w2r(ixc^s,mom(ip1))=sm(ixc^s)
531 w2l(ixc^s,mom(ip1))=sm(ixc^s)
533 w1r(ixc^s,mag(ip1))=bx(ixc^s)
534 w1l(ixc^s,mag(ip1))=bx(ixc^s)
536 ptr(ixc^s)=wrp(ixc^s,p_)+0.5d0*sum(wrc(ixc^s,mag(:))**2,dim=ndim+1)
537 ptl(ixc^s)=wlp(ixc^s,p_)+0.5d0*sum(wlc(ixc^s,mag(:))**2,dim=ndim+1)
541 w1r(ixc^s,rho_)=sur(ixc^s)/(sr(ixc^s,ii)-sm(ixc^s))
542 w1l(ixc^s,rho_)=sul(ixc^s)/(sl(ixc^s,ii)-sm(ixc^s))
546 r1r(ixc^s)=sur(ixc^s)*(sr(ixc^s,ii)-sm(ixc^s))-bx(ixc^s)**2
547 where(abs(r1r(ixc^s))>smalldouble)
548 r1r(ixc^s)=1.d0/r1r(ixc^s)
552 r1l(ixc^s)=sul(ixc^s)*(sl(ixc^s,ii)-sm(ixc^s))-bx(ixc^s)**2
553 where(abs(r1l(ixc^s))>smalldouble)
554 r1l(ixc^s)=1.d0/r1l(ixc^s)
559 w1r(ixc^s,mom(ip2))=vrc(ixc^s,ip2)-bx(ixc^s)*br(ixc^s,ip2)*&
560 (sm(ixc^s)-vrc(ixc^s,ip1))*r1r(ixc^s)
561 w1l(ixc^s,mom(ip2))=vlc(ixc^s,ip2)-bx(ixc^s)*bl(ixc^s,ip2)*&
562 (sm(ixc^s)-vlc(ixc^s,ip1))*r1l(ixc^s)
564 w1r(ixc^s,mag(ip2))=(sur(ixc^s)*(sr(ixc^s,ii)-vrc(ixc^s,ip1))-bx(ixc^s)**2)*r1r(ixc^s)
565 w1l(ixc^s,mag(ip2))=(sul(ixc^s)*(sl(ixc^s,ii)-vlc(ixc^s,ip1))-bx(ixc^s)**2)*r1l(ixc^s)
570 w1r(ixc^s,mom(ip3))=vrc(ixc^s,ip3)-bx(ixc^s)*br(ixc^s,ip3)*&
571 (sm(ixc^s)-vrc(ixc^s,ip1))*r1r(ixc^s)
572 w1l(ixc^s,mom(ip3))=vlc(ixc^s,ip3)-bx(ixc^s)*bl(ixc^s,ip3)*&
573 (sm(ixc^s)-vlc(ixc^s,ip1))*r1l(ixc^s)
575 w1r(ixc^s,mag(ip3))=br(ixc^s,ip3)*w1r(ixc^s,mag(ip2))
576 w1l(ixc^s,mag(ip3))=bl(ixc^s,ip3)*w1l(ixc^s,mag(ip2))
579 w1r(ixc^s,mag(ip2))=br(ixc^s,ip2)*w1r(ixc^s,mag(ip2))
580 w1l(ixc^s,mag(ip2))=bl(ixc^s,ip2)*w1l(ixc^s,mag(ip2))
583 w1r(ixc^s,mag(:))=w1r(ixc^s,mag(:))-block%B0(ixc^s,:,ip1)
584 w1l(ixc^s,mag(:))=w1l(ixc^s,mag(:))-block%B0(ixc^s,:,ip1)
589 w1r(ixc^s,p_)=sur(ixc^s)*(sm(ixc^s)-vrc(ixc^s,ip1))+ptr(ixc^s)
591 w1l(ixc^s,p_)=w1r(ixc^s,p_)
594 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)
595 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)
598 w1r(ixc^s,e_)=((sr(ixc^s,ii)-vrc(ixc^s,ip1))*wrc(ixc^s,e_)-ptr(ixc^s)*vrc(ixc^s,ip1)+&
599 w1r(ixc^s,p_)*sm(ixc^s)+bx(ixc^s)*(sum(vrc(ixc^s,:)*wrc(ixc^s,mag(:)),dim=ndim+1)-&
600 sum(w1r(ixc^s,mom(:))*w1r(ixc^s,mag(:)),dim=ndim+1)))/(sr(ixc^s,ii)-sm(ixc^s))
601 w1l(ixc^s,e_)=((sl(ixc^s,ii)-vlc(ixc^s,ip1))*wlc(ixc^s,e_)-ptl(ixc^s)*vlc(ixc^s,ip1)+&
602 w1l(ixc^s,p_)*sm(ixc^s)+bx(ixc^s)*(sum(vlc(ixc^s,:)*wlc(ixc^s,mag(:)),dim=ndim+1)-&
603 sum(w1l(ixc^s,mom(:))*w1l(ixc^s,mag(:)),dim=ndim+1)))/(sl(ixc^s,ii)-sm(ixc^s))
606 w1r(ixc^s,e_)=w1r(ixc^s,e_)+(sum(w1r(ixc^s,mag(:))*block%B0(ixc^s,:,ip1),dim=ndim+1)*sm(ixc^s)-&
607 sum(wrc(ixc^s,mag(:))*block%B0(ixc^s,:,ip1),dim=ndim+1)*vrc(ixc^s,ip1))/(sr(ixc^s,ii)-sm(ixc^s))
608 w1l(ixc^s,e_)=w1l(ixc^s,e_)+(sum(w1l(ixc^s,mag(:))*block%B0(ixc^s,:,ip1),dim=ndim+1)*sm(ixc^s)-&
609 sum(wlc(ixc^s,mag(:))*block%B0(ixc^s,:,ip1),dim=ndim+1)*vlc(ixc^s,ip1))/(sl(ixc^s,ii)-sm(ixc^s))
612 #if !defined(E_RM_W0) || E_RM_W0 == 1
613 w1r(ixc^s,e_)= w1r(ixc^s,e_) + 1d0/(phys_gamma - 1) * block%equi_vars(ixc^s,iw_equi_p,ip1) * &
614 (sm(ixc^s)-vrc(ixc^s,ip1))/(sr(ixc^s,ii)-sm(ixc^s))
615 w1l(ixc^s,e_)= w1l(ixc^s,e_) + 1d0/(phys_gamma - 1) * block%equi_vars(ixc^s,iw_equi_p,ip1) * &
616 (sm(ixc^s)-vlc(ixc^s,ip1))/(sl(ixc^s,ii)-sm(ixc^s))
618 w1r(ixc^s,e_)= w1r(ixc^s,e_) + phys_gamma /(phys_gamma - 1) * block%equi_vars(ixc^s,iw_equi_p,ip1) * &
619 (sm(ixc^s)-vrc(ixc^s,ip1))/(sr(ixc^s,ii)-sm(ixc^s))
620 w1l(ixc^s,e_)= w1l(ixc^s,e_) + phys_gamma /(phys_gamma - 1) * block%equi_vars(ixc^s,iw_equi_p,ip1) * &
621 (sm(ixc^s)-vlc(ixc^s,ip1))/(sl(ixc^s,ii)-sm(ixc^s))
627 w2r(ixc^s,rho_)=w1r(ixc^s,rho_)
628 w2l(ixc^s,rho_)=w1l(ixc^s,rho_)
629 w2r(ixc^s,mag(ip1))=w1r(ixc^s,mag(ip1))
630 w2l(ixc^s,mag(ip1))=w1l(ixc^s,mag(ip1))
632 r1r(ixc^s)=sqrt(w1r(ixc^s,rho_))
633 r1l(ixc^s)=sqrt(w1l(ixc^s,rho_))
634 tmp(ixc^s)=1.d0/(r1r(ixc^s)+r1l(ixc^s))
635 signbx(ixc^s)=sign(1.d0,bx(ixc^s))
637 s1r(ixc^s)=sm(ixc^s)+abs(bx(ixc^s))/r1r(ixc^s)
638 s1l(ixc^s)=sm(ixc^s)-abs(bx(ixc^s))/r1l(ixc^s)
640 w2r(ixc^s,mom(ip2))=(r1l(ixc^s)*w1l(ixc^s,mom(ip2))+r1r(ixc^s)*w1r(ixc^s,mom(ip2))+&
641 (w1r(ixc^s,mag(ip2))-w1l(ixc^s,mag(ip2)))*signbx(ixc^s))*tmp(ixc^s)
642 w2l(ixc^s,mom(ip2))=w2r(ixc^s,mom(ip2))
644 w2r(ixc^s,mag(ip2))=(r1l(ixc^s)*w1r(ixc^s,mag(ip2))+r1r(ixc^s)*w1l(ixc^s,mag(ip2))+&
645 r1l(ixc^s)*r1r(ixc^s)*(w1r(ixc^s,mom(ip2))-w1l(ixc^s,mom(ip2)))*signbx(ixc^s))*tmp(ixc^s)
646 w2l(ixc^s,mag(ip2))=w2r(ixc^s,mag(ip2))
649 w2r(ixc^s,mom(ip3))=(r1l(ixc^s)*w1l(ixc^s,mom(ip3))+r1r(ixc^s)*w1r(ixc^s,mom(ip3))+&
650 (w1r(ixc^s,mag(ip3))-w1l(ixc^s,mag(ip3)))*signbx(ixc^s))*tmp(ixc^s)
651 w2l(ixc^s,mom(ip3))=w2r(ixc^s,mom(ip3))
653 w2r(ixc^s,mag(ip3))=(r1l(ixc^s)*w1r(ixc^s,mag(ip3))+r1r(ixc^s)*w1l(ixc^s,mag(ip3))+&
654 r1l(ixc^s)*r1r(ixc^s)*(w1r(ixc^s,mom(ip3))-w1l(ixc^s,mom(ip3)))*signbx(ixc^s))*tmp(ixc^s)
655 w2l(ixc^s,mag(ip3))=w2r(ixc^s,mag(ip3))
659 w2r(ixc^s,e_)=w1r(ixc^s,e_)+r1r(ixc^s)*(sum(w1r(ixc^s,mom(:))*w1r(ixc^s,mag(:)),dim=ndim+1)-&
660 sum(w2r(ixc^s,mom(:))*w2r(ixc^s,mag(:)),dim=ndim+1))*signbx(ixc^s)
661 w2l(ixc^s,e_)=w1l(ixc^s,e_)-r1l(ixc^s)*(sum(w1l(ixc^s,mom(:))*w1l(ixc^s,mag(:)),dim=ndim+1)-&
662 sum(w2l(ixc^s,mom(:))*w2l(ixc^s,mag(:)),dim=ndim+1))*signbx(ixc^s)
667 w1r(ixc^s,mom(idir))=w1r(ixc^s,mom(idir))*w1r(ixc^s,rho_)
668 w1l(ixc^s,mom(idir))=w1l(ixc^s,mom(idir))*w1l(ixc^s,rho_)
669 w2r(ixc^s,mom(idir))=w2r(ixc^s,mom(idir))*w2r(ixc^s,rho_)
670 w2l(ixc^s,mom(idir))=w2l(ixc^s,mom(idir))*w2l(ixc^s,rho_)
672 if(iw_equi_rho>0)
then
673 w1r(ixc^s,rho_) = w1r(ixc^s,rho_) - block%equi_vars(ixc^s,iw_equi_rho,ip1)
674 w1l(ixc^s,rho_) = w1l(ixc^s,rho_) - block%equi_vars(ixc^s,iw_equi_rho,ip1)
675 w2r(ixc^s,rho_) = w2r(ixc^s,rho_) - block%equi_vars(ixc^s,iw_equi_rho,ip1)
676 w2l(ixc^s,rho_) = w2l(ixc^s,rho_) - block%equi_vars(ixc^s,iw_equi_rho,ip1)
680 if(flux_type(idims, iw) == flux_special)
then
682 f1l(ixc^s,iw)=flc(ixc^s,iw)
683 f1r(ixc^s,iw)=f1l(ixc^s,iw)
684 f2l(ixc^s,iw)=f1l(ixc^s,iw)
685 f2r(ixc^s,iw)=f1l(ixc^s,iw)
686 else if(flux_type(idims, iw) == flux_hll)
then
688 f1l(ixc^s,iw)=(sr(ixc^s,ii)*flc(ixc^s, iw)-sl(ixc^s,ii)*frc(ixc^s, iw) &
689 +sr(ixc^s,ii)*sl(ixc^s,ii)*(wrc(ixc^s,iw)-wlc(ixc^s,iw)))/(sr(ixc^s,ii)-sl(ixc^s,ii))
690 f1r(ixc^s,iw)=f1l(ixc^s,iw)
691 f2l(ixc^s,iw)=f1l(ixc^s,iw)
692 f2r(ixc^s,iw)=f1l(ixc^s,iw)
695 f1l(ixc^s,iw)=flc(ixc^s,iw)+sl(ixc^s,ii)*(w1l(ixc^s,iw)-wlc(ixc^s,iw))
696 f1r(ixc^s,iw)=frc(ixc^s,iw)+sr(ixc^s,ii)*(w1r(ixc^s,iw)-wrc(ixc^s,iw))
697 f2l(ixc^s,iw)=f1l(ixc^s,iw)+s1l(ixc^s)*(w2l(ixc^s,iw)-w1l(ixc^s,iw))
698 f2r(ixc^s,iw)=f1r(ixc^s,iw)+s1r(ixc^s)*(w2r(ixc^s,iw)-w1r(ixc^s,iw))
703 {
do ix^db=ixcmin^db,ixcmax^db\}
704 if(sl(ix^d,ii)>0.d0)
then
705 fc(ix^d,iws:iwe,ip1)=flc(ix^d,iws:iwe)
706 else if(s1l(ix^d)>=0.d0)
then
707 fc(ix^d,iws:iwe,ip1)=f1l(ix^d,iws:iwe)
708 else if(sm(ix^d)>=0.d0)
then
709 fc(ix^d,iws:iwe,ip1)=f2l(ix^d,iws:iwe)
710 else if(s1r(ix^d)>=0.d0)
then
711 fc(ix^d,iws:iwe,ip1)=f2r(ix^d,iws:iwe)
712 else if(sr(ix^d,ii)>=0.d0)
then
713 fc(ix^d,iws:iwe,ip1)=f1r(ix^d,iws:iwe)
714 else if(sr(ix^d,ii)<0.d0)
then
715 fc(ix^d,iws:iwe,ip1)=frc(ix^d,iws:iwe)
720 end subroutine get_riemann_flux_hlld
724 subroutine get_riemann_flux_hlld_mag2()
729 double precision,
dimension(ixI^S,1:nwflux) :: w1R,w1L,f1R,f1L,f2R,f2L
730 double precision,
dimension(ixI^S,1:nwflux) :: w2R,w2L
731 double precision,
dimension(ixI^S) :: sm,s1R,s1L,suR,suL,Bx
732 double precision,
dimension(ixI^S) :: pts,ptR,ptL,signBx,r1L,r1R,tmp
734 double precision,
dimension(ixI^S,ndir) :: vRC, vLC
736 double precision,
dimension(ixI^S,ndir) :: BR, BL
737 integer :: ip1,ip2,ip3,idir,ix^D
738 double precision :: phiPres, thetaSM, du, dv, dw
739 integer :: ixV^L, ixVb^L, ixVc^L, ixVd^L, ixVe^L, ixVf^L
740 integer :: rho_, p_, e_, mom(1:ndir), mag(1:ndir)
741 double precision,
parameter :: aParam = 4d0
749 associate(sr=>cmaxc,sl=>cminc)
756 vrc(ixc^s,:)=wrp(ixc^s,mom(:))
757 vlc(ixc^s,:)=wlp(ixc^s,mom(:))
763 phipres = min(1d0, maxval(max(s1l(ixo^s),s1r(ixo^s)))/maxval(cmaxc(ixo^s,1)))
764 phipres = phipres*(2d0 - phipres)
771 du = minval(wprim(ixv^s,mom(1))-wprim(ixo^s,mom(1)))
807 dv = minval(min(wprim(ixo^s,mom(2))-wprim(ixv^s,mom(2)),&
808 wprim(ixvb^s,mom(2))-wprim(ixo^s,mom(2)),&
809 wprim(ixvc^s,mom(2))-wprim(ixvd^s,mom(2)),&
810 wprim(ixve^s,mom(2))-wprim(ixvc^s,mom(2))&
844 dw = minval(min(wprim(ixo^s,mom(3))-wprim(ixv^s,mom(3)),&
845 wprim(ixvb^s,mom(3))-wprim(ixo^s,mom(3)),&
846 wprim(ixvc^s,mom(3))-wprim(ixvd^s,mom(3)),&
847 wprim(ixve^s,mom(3))-wprim(ixvc^s,mom(3))&
850 thetasm = maxval(cmaxc(ixo^s,1))
852 thetasm = (min(1d0, (thetasm-du)/(thetasm-min(dv,dw))))**aparam
856 br(ixc^s,:)=wrc(ixc^s,mag(:))+block%B0(ixc^s,:,ip1)
857 bl(ixc^s,:)=wlc(ixc^s,mag(:))+block%B0(ixc^s,:,ip1)
859 br(ixc^s,:)=wrc(ixc^s,mag(:))
860 bl(ixc^s,:)=wlc(ixc^s,mag(:))
865 ptr(ixc^s)=wrp(ixc^s,p_)+0.5d0*sum(br(ixc^s,:)**2,dim=ndim+1)
866 ptl(ixc^s)=wlp(ixc^s,p_)+0.5d0*sum(bl(ixc^s,:)**2,dim=ndim+1)
867 sur(ixc^s)=(sr(ixc^s,
index_v_mag)-vrc(ixc^s,ip1))*wrc(ixc^s,rho_)
868 sul(ixc^s)=(sl(ixc^s,
index_v_mag)-vlc(ixc^s,ip1))*wlc(ixc^s,rho_)
870 sm(ixc^s)=(sur(ixc^s)*vrc(ixc^s,ip1)-sul(ixc^s)*vlc(ixc^s,ip1)-&
871 thetasm*(ptr(ixc^s)-ptl(ixc^s)) )/(sur(ixc^s)-sul(ixc^s))
873 w1r(ixc^s,mom(ip1))=sm(ixc^s)
874 w1l(ixc^s,mom(ip1))=sm(ixc^s)
875 w2r(ixc^s,mom(ip1))=sm(ixc^s)
876 w2l(ixc^s,mom(ip1))=sm(ixc^s)
878 w1r(ixc^s,mag(ip1))=bx(ixc^s)
879 w1l(ixc^s,mag(ip1))=bx(ixc^s)
881 ptr(ixc^s)=wrp(ixc^s,p_)+0.5d0*sum(wrc(ixc^s,mag(:))**2,dim=ndim+1)
882 ptl(ixc^s)=wlp(ixc^s,p_)+0.5d0*sum(wlc(ixc^s,mag(:))**2,dim=ndim+1)
886 w1r(ixc^s,rho_)=sur(ixc^s)/(sr(ixc^s,
index_v_mag)-sm(ixc^s))
887 w1l(ixc^s,rho_)=sul(ixc^s)/(sl(ixc^s,
index_v_mag)-sm(ixc^s))
891 r1r(ixc^s)=sur(ixc^s)*(sr(ixc^s,
index_v_mag)-sm(ixc^s))-bx(ixc^s)**2
892 where(r1r(ixc^s)/=0.d0)
893 r1r(ixc^s)=1.d0/r1r(ixc^s)
895 r1l(ixc^s)=sul(ixc^s)*(sl(ixc^s,
index_v_mag)-sm(ixc^s))-bx(ixc^s)**2
896 where(r1l(ixc^s)/=0.d0)
897 r1l(ixc^s)=1.d0/r1l(ixc^s)
900 w1r(ixc^s,mom(ip2))=vrc(ixc^s,ip2)-bx(ixc^s)*br(ixc^s,ip2)*&
901 (sm(ixc^s)-vrc(ixc^s,ip1))*r1r(ixc^s)
902 w1l(ixc^s,mom(ip2))=vlc(ixc^s,ip2)-bx(ixc^s)*bl(ixc^s,ip2)*&
903 (sm(ixc^s)-vlc(ixc^s,ip1))*r1l(ixc^s)
905 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)
906 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)
911 w1r(ixc^s,mom(ip3))=vrc(ixc^s,ip3)-bx(ixc^s)*br(ixc^s,ip3)*&
912 (sm(ixc^s)-vrc(ixc^s,ip1))*r1r(ixc^s)
913 w1l(ixc^s,mom(ip3))=vlc(ixc^s,ip3)-bx(ixc^s)*bl(ixc^s,ip3)*&
914 (sm(ixc^s)-vlc(ixc^s,ip1))*r1l(ixc^s)
916 w1r(ixc^s,mag(ip3))=br(ixc^s,ip3)*w1r(ixc^s,mag(ip2))
917 w1l(ixc^s,mag(ip3))=bl(ixc^s,ip3)*w1l(ixc^s,mag(ip2))
920 w1r(ixc^s,mag(ip2))=br(ixc^s,ip2)*w1r(ixc^s,mag(ip2))
921 w1l(ixc^s,mag(ip2))=bl(ixc^s,ip2)*w1l(ixc^s,mag(ip2))
924 w1r(ixc^s,mag(:))=w1r(ixc^s,mag(:))-block%B0(ixc^s,:,ip1)
925 w1l(ixc^s,mag(:))=w1l(ixc^s,mag(:))-block%B0(ixc^s,:,ip1)
930 w1r(ixc^s,p_)=(sur(ixc^s)*ptl(ixc^s) - sul(ixc^s)*ptr(ixc^s) +&
931 phipres * sur(ixc^s)*sul(ixc^s)*(vrc(ixc^s,ip1)-vlc(ixc^s,ip1)))/&
932 (sur(ixc^s)-sul(ixc^s))
933 w1l(ixc^s,p_)=w1r(ixc^s,p_)
936 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)
937 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)
940 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)+&
941 w1r(ixc^s,p_)*sm(ixc^s)+bx(ixc^s)*(sum(vrc(ixc^s,:)*wrc(ixc^s,mag(:)),dim=ndim+1)-&
942 sum(w1r(ixc^s,mom(:))*w1r(ixc^s,mag(:)),dim=ndim+1)))/(sr(ixc^s,
index_v_mag)-sm(ixc^s))
943 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)+&
944 w1l(ixc^s,p_)*sm(ixc^s)+bx(ixc^s)*(sum(vlc(ixc^s,:)*wlc(ixc^s,mag(:)),dim=ndim+1)-&
945 sum(w1l(ixc^s,mom(:))*w1l(ixc^s,mag(:)),dim=ndim+1)))/(sl(ixc^s,
index_v_mag)-sm(ixc^s))
948 w1r(ixc^s,e_)=w1r(ixc^s,e_)+(sum(w1r(ixc^s,mag(:))*block%B0(ixc^s,:,ip1),dim=ndim+1)*sm(ixc^s)-&
949 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))
950 w1l(ixc^s,e_)=w1l(ixc^s,e_)+(sum(w1l(ixc^s,mag(:))*block%B0(ixc^s,:,ip1),dim=ndim+1)*sm(ixc^s)-&
951 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))
956 w2r(ixc^s,rho_)=w1r(ixc^s,rho_)
957 w2l(ixc^s,rho_)=w1l(ixc^s,rho_)
958 w2r(ixc^s,mag(ip1))=w1r(ixc^s,mag(ip1))
959 w2l(ixc^s,mag(ip1))=w1l(ixc^s,mag(ip1))
961 r1r(ixc^s)=sqrt(w1r(ixc^s,rho_))
962 r1l(ixc^s)=sqrt(w1l(ixc^s,rho_))
963 tmp(ixc^s)=1.d0/(r1r(ixc^s)+r1l(ixc^s))
964 signbx(ixc^s)=sign(1.d0,bx(ixc^s))
966 s1r(ixc^s)=sm(ixc^s)+abs(bx(ixc^s))/r1r(ixc^s)
967 s1l(ixc^s)=sm(ixc^s)-abs(bx(ixc^s))/r1l(ixc^s)
969 w2r(ixc^s,mom(ip2))=(r1l(ixc^s)*w1l(ixc^s,mom(ip2))+r1r(ixc^s)*w1r(ixc^s,mom(ip2))+&
970 (w1r(ixc^s,mag(ip2))-w1l(ixc^s,mag(ip2)))*signbx(ixc^s))*tmp(ixc^s)
971 w2l(ixc^s,mom(ip2))=w2r(ixc^s,mom(ip2))
973 w2r(ixc^s,mag(ip2))=(r1l(ixc^s)*w1r(ixc^s,mag(ip2))+r1r(ixc^s)*w1l(ixc^s,mag(ip2))+&
974 r1l(ixc^s)*r1r(ixc^s)*(w1r(ixc^s,mom(ip2))-w1l(ixc^s,mom(ip2)))*signbx(ixc^s))*tmp(ixc^s)
975 w2l(ixc^s,mag(ip2))=w2r(ixc^s,mag(ip2))
978 w2r(ixc^s,mom(ip3))=(r1l(ixc^s)*w1l(ixc^s,mom(ip3))+r1r(ixc^s)*w1r(ixc^s,mom(ip3))+&
979 (w1r(ixc^s,mag(ip3))-w1l(ixc^s,mag(ip3)))*signbx(ixc^s))*tmp(ixc^s)
980 w2l(ixc^s,mom(ip3))=w2r(ixc^s,mom(ip3))
982 w2r(ixc^s,mag(ip3))=(r1l(ixc^s)*w1r(ixc^s,mag(ip3))+r1r(ixc^s)*w1l(ixc^s,mag(ip3))+&
983 r1l(ixc^s)*r1r(ixc^s)*(w1r(ixc^s,mom(ip3))-w1l(ixc^s,mom(ip3)))*signbx(ixc^s))*tmp(ixc^s)
984 w2l(ixc^s,mag(ip3))=w2r(ixc^s,mag(ip3))
988 w2r(ixc^s,e_)=w1r(ixc^s,e_)+r1r(ixc^s)*(sum(w1r(ixc^s,mom(:))*w1r(ixc^s,mag(:)),dim=ndim+1)-&
989 sum(w2r(ixc^s,mom(:))*w2r(ixc^s,mag(:)),dim=ndim+1))*signbx(ixc^s)
990 w2l(ixc^s,e_)=w1l(ixc^s,e_)-r1l(ixc^s)*(sum(w1l(ixc^s,mom(:))*w1l(ixc^s,mag(:)),dim=ndim+1)-&
991 sum(w2l(ixc^s,mom(:))*w2l(ixc^s,mag(:)),dim=ndim+1))*signbx(ixc^s)
996 w1r(ixc^s,mom(idir))=w1r(ixc^s,mom(idir))*w1r(ixc^s,rho_)
997 w1l(ixc^s,mom(idir))=w1l(ixc^s,mom(idir))*w1l(ixc^s,rho_)
998 w2r(ixc^s,mom(idir))=w2r(ixc^s,mom(idir))*w2r(ixc^s,rho_)
999 w2l(ixc^s,mom(idir))=w2l(ixc^s,mom(idir))*w2l(ixc^s,rho_)
1008 f1l(ixc^s,iw)=flc(ixc^s,iw)
1009 f1r(ixc^s,iw)=f1l(ixc^s,iw)
1010 f2l(ixc^s,iw)=f1l(ixc^s,iw)
1011 f2r(ixc^s,iw)=f1l(ixc^s,iw)
1016 f1r(ixc^s,iw)=f1l(ixc^s,iw)
1017 f2l(ixc^s,iw)=f1l(ixc^s,iw)
1018 f2r(ixc^s,iw)=f1l(ixc^s,iw)
1020 f1l(ixc^s,iw)=flc(ixc^s,iw)+sl(ixc^s,
index_v_mag)*(w1l(ixc^s,iw)-wlc(ixc^s,iw))
1021 f1r(ixc^s,iw)=frc(ixc^s,iw)+sr(ixc^s,
index_v_mag)*(w1r(ixc^s,iw)-wrc(ixc^s,iw))
1022 f2l(ixc^s,iw)=f1l(ixc^s,iw)+s1l(ixc^s)*(w2l(ixc^s,iw)-w1l(ixc^s,iw))
1023 f2r(ixc^s,iw)=f1r(ixc^s,iw)+s1r(ixc^s)*(w2r(ixc^s,iw)-w1r(ixc^s,iw))
1028 {
do ix^db=ixcmin^db,ixcmax^db\}
1031 else if(s1l(ix^d)>=0.d0)
then
1033 else if(sm(ix^d)>=0.d0)
then
1035 else if(s1r(ix^d)>=0.d0)
then
1045 end subroutine get_riemann_flux_hlld_mag2
1051 integer,
intent(in) :: ixI^L, ixO^L
1052 double precision,
intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
1053 double precision,
intent(out):: csound(ixI^S)
1054 double precision :: cfast2(ixI^S), AvMinCs2(ixI^S), b2(ixI^S), kmax
1055 double precision :: inv_rho(ixO^S), gamma_A2(ixO^S)
1056 integer :: rho_, p_, e_, mom(1:ndir), mag(1:ndir)
1064 inv_rho=1.d0/w(ixo^s,rho_)
1069 b2(ixo^s) = sum((w(ixo^s, mag(:))+
block%B0(ixo^s,:,
b0i))**2, dim=ndim+1)
1071 b2(ixo^s) = sum(w(ixo^s, mag(:))**2, dim=ndim+1)
1076 avmincs2= w(ixo^s, mag(idims))+
block%B0(ixo^s,idims,
b0i)
1078 avmincs2= w(ixo^s, mag(idims))
1082 csound(ixo^s) = sum(w(ixo^s, mom(:))**2, dim=ndim+1)
1084 cfast2(ixo^s) = b2(ixo^s) * inv_rho+csound(ixo^s)
1085 avmincs2(ixo^s) = cfast2(ixo^s)**2-4.0d0*csound(ixo^s) &
1086 * avmincs2(ixo^s)**2 &
1089 where(avmincs2(ixo^s)<zero)
1090 avmincs2(ixo^s)=zero
1093 avmincs2(ixo^s)=sqrt(avmincs2(ixo^s))
1095 csound(ixo^s) = sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s)))
1103 subroutine reconstruct_lr(ixI^L,ixL^L,ixR^L,idims,w,wLC,wRC,wLp,wRp,x,dxdim)
1109 integer,
intent(in) :: ixi^
l, ixl^
l, ixr^
l, idims
1110 double precision,
intent(in) :: dxdim
1112 double precision,
dimension(ixI^S,1:nw) :: w
1114 double precision,
dimension(ixI^S,1:nw) :: wlc, wrc
1116 double precision,
dimension(ixI^S,1:nw) :: wlp, wrp
1117 double precision,
dimension(ixI^S,1:ndim) :: x
1119 integer :: jxr^
l, ixc^
l, jxc^
l, iw
1120 double precision :: ldw(ixi^s), rdw(ixi^s), dwc(ixi^s)
1121 double precision :: a2max
1125 call mp5limiter(ixi^
l,ixl^
l,idims,w,wlp,wrp)
1127 call weno3limiter(ixi^
l,ixl^
l,idims,dxdim,w,wlp,wrp,1)
1129 call weno3limiter(ixi^
l,ixl^
l,idims,dxdim,w,wlp,wrp,2)
1131 call weno5limiter(ixi^
l,ixl^
l,idims,dxdim,w,wlp,wrp,1)
1133 call weno5nmlimiter(ixi^
l,ixl^
l,idims,dxdim,w,wlp,wrp,1)
1135 call weno5limiter(ixi^
l,ixl^
l,idims,dxdim,w,wlp,wrp,2)
1137 call weno5nmlimiter(ixi^
l,ixl^
l,idims,dxdim,w,wlp,wrp,2)
1139 call weno5limiter(ixi^
l,ixl^
l,idims,dxdim,w,wlp,wrp,3)
1141 call weno5nmlimiter(ixi^
l,ixl^
l,idims,dxdim,w,wlp,wrp,3)
1143 call weno5cu6limiter(ixi^
l,ixl^
l,idims,w,wlp,wrp)
1145 call teno5adlimiter(ixi^
l,ixl^
l,idims,dxdim,w,wlp,wrp)
1147 call weno7limiter(ixi^
l,ixl^
l,idims,w,wlp,wrp,1)
1149 call weno7limiter(ixi^
l,ixl^
l,idims,w,wlp,wrp,2)
1151 call venklimiter(ixi^
l,ixl^
l,idims,dxdim,w,wlp,wrp)
1157 call ppmlimiter(ixi^
l,
ixm^
ll,idims,w,w,wlp,wrp)
1163 jxr^
l=ixr^
l+
kr(idims,^
d);
1164 ixcmax^
d=jxrmax^
d; ixcmin^
d=ixlmin^
d-
kr(idims,^
d);
1165 jxc^
l=ixc^
l+
kr(idims,^
d);
1168 w(ixcmin^
d:jxcmax^
d,iw)=dlog10(w(ixcmin^
d:jxcmax^
d,iw))
1169 wlp(ixl^s,iw)=dlog10(wlp(ixl^s,iw))
1170 wrp(ixr^s,iw)=dlog10(wrp(ixr^s,iw))
1173 dwc(ixc^s)=w(jxc^s,iw)-w(ixc^s,iw)
1189 call mpistop(
"idims is wrong in mod_limiter")
1195 wlp(ixl^s,iw)=wlp(ixl^s,iw)+half*ldw(ixl^s)
1196 wrp(ixr^s,iw)=wrp(ixr^s,iw)-half*rdw(jxr^s)
1199 w(ixcmin^
d:jxcmax^
d,iw)=10.0d0**w(ixcmin^
d:jxcmax^
d,iw)
1200 wlp(ixl^s,iw)=10.0d0**wlp(ixl^s,iw)
1201 wrp(ixr^s,iw)=10.0d0**wrp(ixr^s,iw)
1210 wlc(ixl^s,1:nwflux) = wlp(ixl^s,1:nwflux)
1211 wrc(ixr^s,1:nwflux) = wrp(ixr^s,1:nwflux)
1215 wlp(ixl^s,nwflux+1:nwflux+nwaux) = wlc(ixl^s,nwflux+1:nwflux+nwaux)
1216 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 crash
Save a snapshot before crash a run met unphysical values.
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_face_to_center), pointer phys_face_to_center
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.