19 subroutine hancock(qdt,ixI^L,ixO^L,idims^LIM,qtC,sCT,qt,snew,dxs,x)
24 integer,
intent(in) :: ixi^
l, ixo^
l, idims^lim
25 double precision,
intent(in) :: qdt, qtc, qt, dxs(
ndim), x(ixi^s,1:
ndim)
26 type(state) :: sct, snew
28 double precision,
dimension(ixI^S,1:nw) :: wprim, wlc, wrc
30 double precision,
dimension(ixI^S,1:nw) :: wlp, wrp
31 double precision,
dimension(ixO^S) :: inv_volume
32 double precision :: flc(ixi^s, nwflux), frc(ixi^s, nwflux)
33 double precision :: dxinv(1:
ndim)
34 integer :: idims, iw, ix^
l, hxo^
l
37 associate(wct=>sct%w,wnew=>snew%w)
41 ix^
l=ix^
l^laddkr(idims,^
d);
43 if (ixi^
l^ltix^
l|.or.|.or.) &
44 call mpistop(
"Error in Hancock: Nonconforming input limits")
56 hxo^
l=ixo^
l-
kr(idims,^
d);
58 wrp(hxo^s,1:nwflux)=wprim(ixo^s,1:nwflux)
59 wlp(ixo^s,1:nwflux)=wprim(ixo^s,1:nwflux)
62 call reconstruct_lr(ixi^
l,ixo^
l,hxo^
l,idims,wprim,wlc,wrc,wlp,wrp,x,dxs(idims))
71 wnew(ixo^s,iw)=wnew(ixo^s,iw)+dxinv(idims)* &
72 (flc(ixo^s, iw)-frc(hxo^s, iw))
76 wnew(ixo^s,iw)=wnew(ixo^s,iw) - qdt * inv_volume &
77 *(
block%surfaceC(ixo^s,idims)*flc(ixo^s, iw) &
78 -
block%surfaceC(hxo^s,idims)*frc(hxo^s, iw))
87 ixi^
l,ixo^
l,1,nw,qtc,wct,qt,wnew,x,.false.,active,wprim)
98 qtC,sCT,qt,snew,sold,fC,fE,dxs,x)
106 integer,
intent(in) :: method
107 double precision,
intent(in) :: qdt, qtc, qt, dxs(
ndim)
108 integer,
intent(in) :: ixi^
l, ixo^
l, idims^lim
109 double precision,
dimension(ixI^S,1:ndim),
intent(in) :: x
110 type(state) :: sct, snew, sold
111 double precision,
dimension(ixI^S,1:nwflux,1:ndim) :: fc
112 double precision,
dimension(ixI^S,7-2*ndim:3) :: fe
115 double precision,
dimension(ixI^S,1:nw) :: wprim
117 double precision,
dimension(ixI^S,1:nw) :: wlc, wrc
119 double precision,
dimension(ixI^S,1:nw) :: wlp, wrp
120 double precision,
dimension(ixI^S,1:nwflux) :: flc, frc
121 double precision,
dimension(ixI^S,1:number_species) :: cmaxc
122 double precision,
dimension(ixI^S,1:number_species) :: cminc
123 double precision,
dimension(ixI^S) :: hspeed
124 double precision,
dimension(ixO^S) :: inv_volume
125 double precision,
dimension(1:ndim) :: dxinv
126 integer,
dimension(ixI^S) :: patchf
127 integer :: idims, iw, ix^
l, hxo^
l, ixc^
l, ixcr^
l, kxc^
l, kxr^
l, ii
129 type(ct_velocity) :: vcts
131 associate(wct=>sct%w, wnew=>snew%w, wold=>sold%w)
141 ix^
l=ix^
l^ladd2*
kr(idims,^
d);
143 if (ixi^
l^ltix^
l|.or.|.or.) &
144 call mpistop(
"Error in fv : Nonconforming input limits")
153 hxo^
l=ixo^
l-
kr(idims,^
d);
155 kxcmin^
d=iximin^
d; kxcmax^
d=iximax^
d-
kr(idims,^
d);
156 kxr^
l=kxc^
l+
kr(idims,^
d);
164 ixcmax^
d=ixomax^
d; ixcmin^
d=hxomin^
d;
171 wrp(kxc^s,1:
nw)=wprim(kxr^s,1:
nw)
172 wlp(kxc^s,1:
nw)=wprim(kxc^s,1:
nw)
179 call reconstruct_lr(ixi^
l,ixcr^
l,ixcr^
l,idims,wprim,wlc,wrc,wlp,wrp,x,dxs(idims))
197 call phys_get_cbounds(wlc,wrc,wlp,wrp,x,ixi^
l,ixc^
l,idims,hspeed,cmaxc,cminc)
226 call mpistop(
'unkown Riemann flux in finite volume')
237 hxo^
l=ixo^
l-
kr(idims,^
d);
240 fc(ixi^s,1:
nwflux,idims)=dxinv(idims)*fc(ixi^s,1:
nwflux,idims)
247 call tvdlimit2(method,qdt,ixi^
l,ixc^
l,ixo^
l,idims,wlc,wrc,wnew,x,fc,dxs)
251 inv_volume = 1.d0/
block%dvolume(ixo^s)
253 hxo^
l=ixo^
l-
kr(idims,^
d);
257 fc(ixi^s,iw,idims)=-qdt*fc(ixi^s,iw,idims)*
block%surfaceC(ixi^s,idims)
258 wnew(ixo^s,iw)=wnew(ixo^s,iw) + (fc(ixo^s,iw,idims)-fc(hxo^s,iw,idims)) * &
269 call tvdlimit2(method,qdt,ixi^
l,ixc^
l,ixo^
l,idims,wlc,wrc,wnew,x,fc,dxs)
274 if (.not.
slab.and.idimsmin==1) &
285 ixi^
l,ixo^
l,1,
nw,qtc,wct,qt,wnew,x,.false.,active,wprim)
298 flc(ixc^s, iw)=half*(flc(ixc^s, iw)+frc(ixc^s, iw))
299 fc(ixc^s,iw,idims)=flc(ixc^s, iw)
303 subroutine get_riemann_flux_tvdlf(iws,iwe)
304 integer,
intent(in) :: iws,iwe
305 double precision :: fac(ixC^S)
307 fac(ixc^s) = -0.5d0*tvdlfeps*cmaxc(ixc^s,ii)
311 flc(ixc^s, iw)=0.5d0*(flc(ixc^s, iw)+frc(ixc^s, iw))
313 if (flux_type(idims, iw) /= flux_no_dissipation)
then
314 flc(ixc^s, iw)=flc(ixc^s, iw) + fac(ixc^s)*(wrc(ixc^s,iw)-wlc(ixc^s,iw))
316 fc(ixc^s,iw,idims)=flc(ixc^s, iw)
319 end subroutine get_riemann_flux_tvdlf
321 subroutine get_riemann_flux_hll(iws,iwe)
322 integer,
intent(in) :: iws,iwe
326 if(flux_type(idims, iw) == flux_tvdlf)
then
328 if(stagger_grid) cycle
329 fc(ixc^s,iw,idims) = -tvdlfeps*half*max(cmaxc(ixc^s,ii),dabs(cminc(ixc^s,ii))) * &
330 (wrc(ixc^s,iw)-wlc(ixc^s,iw))
332 {
do ix^db=ixcmin^db,ixcmax^db\}
333 if(cminc(ix^d,ii) >= zero)
then
334 fc(ix^d,iw,idims)=flc(ix^d,iw)
335 else if(cmaxc(ix^d,ii) <= zero)
then
336 fc(ix^d,iw,idims)=frc(ix^d,iw)
339 fc(ix^d,iw,idims)=(cmaxc(ix^d,ii)*flc(ix^d, iw)-cminc(ix^d,ii)*frc(ix^d,iw)&
340 +cminc(ix^d,ii)*cmaxc(ix^d,ii)*(wrc(ix^d,iw)-wlc(ix^d,iw)))&
341 /(cmaxc(ix^d,ii)-cminc(ix^d,ii))
347 end subroutine get_riemann_flux_hll
349 subroutine get_riemann_flux_hllc(iws,iwe)
350 integer,
intent(in) :: iws, iwe
351 double precision,
dimension(ixI^S,1:nwflux) :: whll, Fhll, fCD
352 double precision,
dimension(ixI^S) :: lambdaCD
354 integer :: rho_, p_, e_, eaux_, mom(1:ndir)
357 if (
allocated(iw_mom)) mom(:) = iw_mom(:)
361 if(
associated(phys_hllc_init_species))
then
362 call phys_hllc_init_species(ii, rho_, mom(:), e_, eaux_)
368 where(cminc(ixc^s,1) >= zero)
370 elsewhere(cmaxc(ixc^s,1) <= zero)
374 if(method==fs_hllcd) &
375 call phys_diffuse_hllcd(ixi^l,ixc^l,idims,wlc,wrc,flc,frc,patchf)
378 if(any(patchf(ixc^s)==1)) &
379 call phys_get_lcd(wlc,wrc,flc,frc,cminc(ixi^s,ii),cmaxc(ixi^s,ii),idims,ixi^l,ixc^l, &
380 whll,fhll,lambdacd,patchf)
383 if(any(abs(patchf(ixc^s))== 1))
then
385 call phys_get_wcd(wlc,wrc,whll,frc,flc,fhll,patchf,lambdacd,&
386 cminc(ixi^s,ii),cmaxc(ixi^s,ii),ixi^l,ixc^l,idims,fcd)
390 if(phys_energy.and.phys_solve_eaux .and. eaux_>0)
then
392 fcd(ixc^s, iw) = (cmaxc(ixc^s,ii)*flc(ixc^s, iw)-cminc(ixc^s,ii) * frc(ixc^s, iw) &
393 +cminc(ixc^s,ii)*cmaxc(ixc^s,ii)*(wrc(ixc^s,iw)-wlc(ixc^s,iw)))/(cmaxc(ixc^s,ii)-cminc(ixc^s,ii))
397 if (flux_type(idims, iw) == flux_tvdlf)
then
398 flc(ixc^s,iw)=-tvdlfeps*half*max(cmaxc(ixc^s,ii),abs(cminc(ixc^s,ii))) * &
399 (wrc(ixc^s,iw) - wlc(ixc^s,iw))
401 where(patchf(ixc^s)==-2)
402 flc(ixc^s,iw)=flc(ixc^s,iw)
403 elsewhere(abs(patchf(ixc^s))==1)
404 flc(ixc^s,iw)=fcd(ixc^s,iw)
405 elsewhere(patchf(ixc^s)==2)
406 flc(ixc^s,iw)=frc(ixc^s,iw)
407 elsewhere(patchf(ixc^s)==3)
409 flc(ixc^s,iw)=fhll(ixc^s,iw)
410 elsewhere(patchf(ixc^s)==4)
412 flc(ixc^s,iw) = half*((flc(ixc^s,iw)+frc(ixc^s,iw)) &
413 -tvdlfeps * max(cmaxc(ixc^s,ii), dabs(cminc(ixc^s,ii))) * &
414 (wrc(ixc^s,iw)-wlc(ixc^s,iw)))
418 fc(ixc^s,iw,idims)=flc(ixc^s,iw)
421 end subroutine get_riemann_flux_hllc
424 subroutine get_riemann_flux_hlld(iws,iwe)
425 integer,
intent(in) :: iws, iwe
426 double precision,
dimension(ixI^S,1:nwflux) :: w1R,w1L,f1R,f1L,f2R,f2L
427 double precision,
dimension(ixI^S,1:nwflux) :: w2R,w2L
428 double precision,
dimension(ixI^S) :: sm,s1R,s1L,suR,suL,Bx
429 double precision,
dimension(ixI^S) :: pts,ptR,ptL,signBx,r1L,r1R,tmp
431 double precision,
dimension(ixI^S,ndir) :: vRC, vLC
433 double precision,
dimension(ixI^S,ndir) :: BR, BL
434 integer :: ip1,ip2,ip3,idir,ix^D
435 integer :: rho_, p_, e_, eaux_, mom(1:ndir), mag(1:ndir)
437 associate(sr=>cmaxc,sl=>cminc)
457 vrc(ixc^s,:)=wrp(ixc^s,mom(:))
458 vlc(ixc^s,:)=wlp(ixc^s,mom(:))
460 br(ixc^s,:)=wrc(ixc^s,mag(:))+block%B0(ixc^s,:,ip1)
461 bl(ixc^s,:)=wlc(ixc^s,mag(:))+block%B0(ixc^s,:,ip1)
463 br(ixc^s,:)=wrc(ixc^s,mag(:))
464 bl(ixc^s,:)=wlc(ixc^s,mag(:))
466 if(stagger_grid)
then
467 bx(ixc^s)=block%ws(ixc^s,ip1)
471 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))
473 ptr(ixc^s)=wrp(ixc^s,p_)+0.5d0*sum(br(ixc^s,:)**2,dim=ndim+1)
474 ptl(ixc^s)=wlp(ixc^s,p_)+0.5d0*sum(bl(ixc^s,:)**2,dim=ndim+1)
475 if(iw_equi_rho>0)
then
476 sur(ixc^s) = wrc(ixc^s,rho_)+ block%equi_vars(ixc^s,iw_equi_rho,ip1)
478 sur(ixc^s) = wrc(ixc^s,rho_)
480 sur(ixc^s)=(sr(ixc^s,ii)-vrc(ixc^s,ip1))*sur(ixc^s)
481 if(iw_equi_rho>0)
then
482 sul(ixc^s) = wlc(ixc^s,rho_)+ block%equi_vars(ixc^s,iw_equi_rho,ip1)
484 sul(ixc^s) = wlc(ixc^s,rho_)
486 sul(ixc^s)=(sl(ixc^s,ii)-vlc(ixc^s,ip1))*sul(ixc^s)
488 sm(ixc^s)=(sur(ixc^s)*vrc(ixc^s,ip1)-sul(ixc^s)*vlc(ixc^s,ip1)-&
489 ptr(ixc^s)+ptl(ixc^s))/(sur(ixc^s)-sul(ixc^s))
491 w1r(ixc^s,mom(ip1))=sm(ixc^s)
492 w1l(ixc^s,mom(ip1))=sm(ixc^s)
493 w2r(ixc^s,mom(ip1))=sm(ixc^s)
494 w2l(ixc^s,mom(ip1))=sm(ixc^s)
496 w1r(ixc^s,mag(ip1))=bx(ixc^s)
497 w1l(ixc^s,mag(ip1))=bx(ixc^s)
499 ptr(ixc^s)=wrp(ixc^s,p_)+0.5d0*sum(wrc(ixc^s,mag(:))**2,dim=ndim+1)
500 ptl(ixc^s)=wlp(ixc^s,p_)+0.5d0*sum(wlc(ixc^s,mag(:))**2,dim=ndim+1)
504 w1r(ixc^s,rho_)=sur(ixc^s)/(sr(ixc^s,ii)-sm(ixc^s))
505 w1l(ixc^s,rho_)=sul(ixc^s)/(sl(ixc^s,ii)-sm(ixc^s))
509 r1r(ixc^s)=sur(ixc^s)*(sr(ixc^s,ii)-sm(ixc^s))-bx(ixc^s)**2
510 where(abs(r1r(ixc^s))>smalldouble)
511 r1r(ixc^s)=1.d0/r1r(ixc^s)
515 r1l(ixc^s)=sul(ixc^s)*(sl(ixc^s,ii)-sm(ixc^s))-bx(ixc^s)**2
516 where(abs(r1l(ixc^s))>smalldouble)
517 r1l(ixc^s)=1.d0/r1l(ixc^s)
522 w1r(ixc^s,mom(ip2))=vrc(ixc^s,ip2)-bx(ixc^s)*br(ixc^s,ip2)*&
523 (sm(ixc^s)-vrc(ixc^s,ip1))*r1r(ixc^s)
524 w1l(ixc^s,mom(ip2))=vlc(ixc^s,ip2)-bx(ixc^s)*bl(ixc^s,ip2)*&
525 (sm(ixc^s)-vlc(ixc^s,ip1))*r1l(ixc^s)
527 w1r(ixc^s,mag(ip2))=(sur(ixc^s)*(sr(ixc^s,ii)-vrc(ixc^s,ip1))-bx(ixc^s)**2)*r1r(ixc^s)
528 w1l(ixc^s,mag(ip2))=(sul(ixc^s)*(sl(ixc^s,ii)-vlc(ixc^s,ip1))-bx(ixc^s)**2)*r1l(ixc^s)
533 w1r(ixc^s,mom(ip3))=vrc(ixc^s,ip3)-bx(ixc^s)*br(ixc^s,ip3)*&
534 (sm(ixc^s)-vrc(ixc^s,ip1))*r1r(ixc^s)
535 w1l(ixc^s,mom(ip3))=vlc(ixc^s,ip3)-bx(ixc^s)*bl(ixc^s,ip3)*&
536 (sm(ixc^s)-vlc(ixc^s,ip1))*r1l(ixc^s)
538 w1r(ixc^s,mag(ip3))=br(ixc^s,ip3)*w1r(ixc^s,mag(ip2))
539 w1l(ixc^s,mag(ip3))=bl(ixc^s,ip3)*w1l(ixc^s,mag(ip2))
542 w1r(ixc^s,mag(ip2))=br(ixc^s,ip2)*w1r(ixc^s,mag(ip2))
543 w1l(ixc^s,mag(ip2))=bl(ixc^s,ip2)*w1l(ixc^s,mag(ip2))
546 w1r(ixc^s,mag(:))=w1r(ixc^s,mag(:))-block%B0(ixc^s,:,ip1)
547 w1l(ixc^s,mag(:))=w1l(ixc^s,mag(:))-block%B0(ixc^s,:,ip1)
552 w1r(ixc^s,p_)=sur(ixc^s)*(sm(ixc^s)-vrc(ixc^s,ip1))+ptr(ixc^s)
554 w1l(ixc^s,p_)=w1r(ixc^s,p_)
557 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)
558 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)
561 w1r(ixc^s,e_)=((sr(ixc^s,ii)-vrc(ixc^s,ip1))*wrc(ixc^s,e_)-ptr(ixc^s)*vrc(ixc^s,ip1)+&
562 w1r(ixc^s,p_)*sm(ixc^s)+bx(ixc^s)*(sum(vrc(ixc^s,:)*wrc(ixc^s,mag(:)),dim=ndim+1)-&
563 sum(w1r(ixc^s,mom(:))*w1r(ixc^s,mag(:)),dim=ndim+1)))/(sr(ixc^s,ii)-sm(ixc^s))
564 w1l(ixc^s,e_)=((sl(ixc^s,ii)-vlc(ixc^s,ip1))*wlc(ixc^s,e_)-ptl(ixc^s)*vlc(ixc^s,ip1)+&
565 w1l(ixc^s,p_)*sm(ixc^s)+bx(ixc^s)*(sum(vlc(ixc^s,:)*wlc(ixc^s,mag(:)),dim=ndim+1)-&
566 sum(w1l(ixc^s,mom(:))*w1l(ixc^s,mag(:)),dim=ndim+1)))/(sl(ixc^s,ii)-sm(ixc^s))
569 w1r(ixc^s,e_)=w1r(ixc^s,e_)+(sum(w1r(ixc^s,mag(:))*block%B0(ixc^s,:,ip1),dim=ndim+1)*sm(ixc^s)-&
570 sum(wrc(ixc^s,mag(:))*block%B0(ixc^s,:,ip1),dim=ndim+1)*vrc(ixc^s,ip1))/(sr(ixc^s,ii)-sm(ixc^s))
571 w1l(ixc^s,e_)=w1l(ixc^s,e_)+(sum(w1l(ixc^s,mag(:))*block%B0(ixc^s,:,ip1),dim=ndim+1)*sm(ixc^s)-&
572 sum(wlc(ixc^s,mag(:))*block%B0(ixc^s,:,ip1),dim=ndim+1)*vlc(ixc^s,ip1))/(sl(ixc^s,ii)-sm(ixc^s))
575 #if !defined(E_RM_W0) || E_RM_W0 == 1
576 w1r(ixc^s,e_)= w1r(ixc^s,e_) + 1d0/(phys_gamma - 1) * block%equi_vars(ixc^s,iw_equi_p,ip1) * &
577 (sm(ixc^s)-vrc(ixc^s,ip1))/(sr(ixc^s,ii)-sm(ixc^s))
578 w1l(ixc^s,e_)= w1l(ixc^s,e_) + 1d0/(phys_gamma - 1) * block%equi_vars(ixc^s,iw_equi_p,ip1) * &
579 (sm(ixc^s)-vlc(ixc^s,ip1))/(sl(ixc^s,ii)-sm(ixc^s))
581 w1r(ixc^s,e_)= w1r(ixc^s,e_) + phys_gamma /(phys_gamma - 1) * block%equi_vars(ixc^s,iw_equi_p,ip1) * &
582 (sm(ixc^s)-vrc(ixc^s,ip1))/(sr(ixc^s,ii)-sm(ixc^s))
583 w1l(ixc^s,e_)= w1l(ixc^s,e_) + phys_gamma /(phys_gamma - 1) * block%equi_vars(ixc^s,iw_equi_p,ip1) * &
584 (sm(ixc^s)-vlc(ixc^s,ip1))/(sl(ixc^s,ii)-sm(ixc^s))
590 w2r(ixc^s,rho_)=w1r(ixc^s,rho_)
591 w2l(ixc^s,rho_)=w1l(ixc^s,rho_)
592 w2r(ixc^s,mag(ip1))=w1r(ixc^s,mag(ip1))
593 w2l(ixc^s,mag(ip1))=w1l(ixc^s,mag(ip1))
595 r1r(ixc^s)=sqrt(w1r(ixc^s,rho_))
596 r1l(ixc^s)=sqrt(w1l(ixc^s,rho_))
597 tmp(ixc^s)=1.d0/(r1r(ixc^s)+r1l(ixc^s))
598 signbx(ixc^s)=sign(1.d0,bx(ixc^s))
600 s1r(ixc^s)=sm(ixc^s)+abs(bx(ixc^s))/r1r(ixc^s)
601 s1l(ixc^s)=sm(ixc^s)-abs(bx(ixc^s))/r1l(ixc^s)
603 w2r(ixc^s,mom(ip2))=(r1l(ixc^s)*w1l(ixc^s,mom(ip2))+r1r(ixc^s)*w1r(ixc^s,mom(ip2))+&
604 (w1r(ixc^s,mag(ip2))-w1l(ixc^s,mag(ip2)))*signbx(ixc^s))*tmp(ixc^s)
605 w2l(ixc^s,mom(ip2))=w2r(ixc^s,mom(ip2))
607 w2r(ixc^s,mag(ip2))=(r1l(ixc^s)*w1r(ixc^s,mag(ip2))+r1r(ixc^s)*w1l(ixc^s,mag(ip2))+&
608 r1l(ixc^s)*r1r(ixc^s)*(w1r(ixc^s,mom(ip2))-w1l(ixc^s,mom(ip2)))*signbx(ixc^s))*tmp(ixc^s)
609 w2l(ixc^s,mag(ip2))=w2r(ixc^s,mag(ip2))
612 w2r(ixc^s,mom(ip3))=(r1l(ixc^s)*w1l(ixc^s,mom(ip3))+r1r(ixc^s)*w1r(ixc^s,mom(ip3))+&
613 (w1r(ixc^s,mag(ip3))-w1l(ixc^s,mag(ip3)))*signbx(ixc^s))*tmp(ixc^s)
614 w2l(ixc^s,mom(ip3))=w2r(ixc^s,mom(ip3))
616 w2r(ixc^s,mag(ip3))=(r1l(ixc^s)*w1r(ixc^s,mag(ip3))+r1r(ixc^s)*w1l(ixc^s,mag(ip3))+&
617 r1l(ixc^s)*r1r(ixc^s)*(w1r(ixc^s,mom(ip3))-w1l(ixc^s,mom(ip3)))*signbx(ixc^s))*tmp(ixc^s)
618 w2l(ixc^s,mag(ip3))=w2r(ixc^s,mag(ip3))
622 w2r(ixc^s,e_)=w1r(ixc^s,e_)+r1r(ixc^s)*(sum(w1r(ixc^s,mom(:))*w1r(ixc^s,mag(:)),dim=ndim+1)-&
623 sum(w2r(ixc^s,mom(:))*w2r(ixc^s,mag(:)),dim=ndim+1))*signbx(ixc^s)
624 w2l(ixc^s,e_)=w1l(ixc^s,e_)-r1l(ixc^s)*(sum(w1l(ixc^s,mom(:))*w1l(ixc^s,mag(:)),dim=ndim+1)-&
625 sum(w2l(ixc^s,mom(:))*w2l(ixc^s,mag(:)),dim=ndim+1))*signbx(ixc^s)
630 w1r(ixc^s,mom(idir))=w1r(ixc^s,mom(idir))*w1r(ixc^s,rho_)
631 w1l(ixc^s,mom(idir))=w1l(ixc^s,mom(idir))*w1l(ixc^s,rho_)
632 w2r(ixc^s,mom(idir))=w2r(ixc^s,mom(idir))*w2r(ixc^s,rho_)
633 w2l(ixc^s,mom(idir))=w2l(ixc^s,mom(idir))*w2l(ixc^s,rho_)
635 if(iw_equi_rho>0)
then
636 w1r(ixc^s,rho_) = w1r(ixc^s,rho_) - block%equi_vars(ixc^s,iw_equi_rho,ip1)
637 w1l(ixc^s,rho_) = w1l(ixc^s,rho_) - block%equi_vars(ixc^s,iw_equi_rho,ip1)
638 w2r(ixc^s,rho_) = w2r(ixc^s,rho_) - block%equi_vars(ixc^s,iw_equi_rho,ip1)
639 w2l(ixc^s,rho_) = w2l(ixc^s,rho_) - block%equi_vars(ixc^s,iw_equi_rho,ip1)
643 if(flux_type(idims, iw) == flux_special)
then
645 f1l(ixc^s,iw)=flc(ixc^s,iw)
646 f1r(ixc^s,iw)=f1l(ixc^s,iw)
647 f2l(ixc^s,iw)=f1l(ixc^s,iw)
648 f2r(ixc^s,iw)=f1l(ixc^s,iw)
649 else if(flux_type(idims, iw) == flux_hll)
then
651 f1l(ixc^s,iw)=(sr(ixc^s,ii)*flc(ixc^s, iw)-sl(ixc^s,ii)*frc(ixc^s, iw) &
652 +sr(ixc^s,ii)*sl(ixc^s,ii)*(wrc(ixc^s,iw)-wlc(ixc^s,iw)))/(sr(ixc^s,ii)-sl(ixc^s,ii))
653 f1r(ixc^s,iw)=f1l(ixc^s,iw)
654 f2l(ixc^s,iw)=f1l(ixc^s,iw)
655 f2r(ixc^s,iw)=f1l(ixc^s,iw)
658 f1l(ixc^s,iw)=flc(ixc^s,iw)+sl(ixc^s,ii)*(w1l(ixc^s,iw)-wlc(ixc^s,iw))
659 f1r(ixc^s,iw)=frc(ixc^s,iw)+sr(ixc^s,ii)*(w1r(ixc^s,iw)-wrc(ixc^s,iw))
660 f2l(ixc^s,iw)=f1l(ixc^s,iw)+s1l(ixc^s)*(w2l(ixc^s,iw)-w1l(ixc^s,iw))
661 f2r(ixc^s,iw)=f1r(ixc^s,iw)+s1r(ixc^s)*(w2r(ixc^s,iw)-w1r(ixc^s,iw))
666 {
do ix^db=ixcmin^db,ixcmax^db\}
667 if(sl(ix^d,ii)>0.d0)
then
668 fc(ix^d,iws:iwe,ip1)=flc(ix^d,iws:iwe)
669 else if(s1l(ix^d)>=0.d0)
then
670 fc(ix^d,iws:iwe,ip1)=f1l(ix^d,iws:iwe)
671 else if(sm(ix^d)>=0.d0)
then
672 fc(ix^d,iws:iwe,ip1)=f2l(ix^d,iws:iwe)
673 else if(s1r(ix^d)>=0.d0)
then
674 fc(ix^d,iws:iwe,ip1)=f2r(ix^d,iws:iwe)
675 else if(sr(ix^d,ii)>=0.d0)
then
676 fc(ix^d,iws:iwe,ip1)=f1r(ix^d,iws:iwe)
677 else if(sr(ix^d,ii)<0.d0)
then
678 fc(ix^d,iws:iwe,ip1)=frc(ix^d,iws:iwe)
683 end subroutine get_riemann_flux_hlld
687 subroutine get_riemann_flux_hlld_mag2()
692 double precision,
dimension(ixI^S,1:nwflux) :: w1R,w1L,f1R,f1L,f2R,f2L
693 double precision,
dimension(ixI^S,1:nwflux) :: w2R,w2L
694 double precision,
dimension(ixI^S) :: sm,s1R,s1L,suR,suL,Bx
695 double precision,
dimension(ixI^S) :: pts,ptR,ptL,signBx,r1L,r1R,tmp
697 double precision,
dimension(ixI^S,ndir) :: vRC, vLC
699 double precision,
dimension(ixI^S,ndir) :: BR, BL
700 integer :: ip1,ip2,ip3,idir,ix^D
701 double precision :: phiPres, thetaSM, du, dv, dw
702 integer :: ixV^L, ixVb^L, ixVc^L, ixVd^L, ixVe^L, ixVf^L
703 integer :: rho_, p_, e_, eaux_, mom(1:ndir), mag(1:ndir)
704 double precision,
parameter :: aParam = 4d0
713 associate(sr=>cmaxc,sl=>cminc)
720 vrc(ixc^s,:)=wrp(ixc^s,mom(:))
721 vlc(ixc^s,:)=wlp(ixc^s,mom(:))
727 phipres = min(1d0, maxval(max(s1l(ixo^s),s1r(ixo^s)))/maxval(cmaxc(ixo^s,1)))
728 phipres = phipres*(2d0 - phipres)
735 du = minval(wprim(ixv^s,mom(1))-wprim(ixo^s,mom(1)))
771 dv = minval(min(wprim(ixo^s,mom(2))-wprim(ixv^s,mom(2)),&
772 wprim(ixvb^s,mom(2))-wprim(ixo^s,mom(2)),&
773 wprim(ixvc^s,mom(2))-wprim(ixvd^s,mom(2)),&
774 wprim(ixve^s,mom(2))-wprim(ixvc^s,mom(2))&
808 dw = minval(min(wprim(ixo^s,mom(3))-wprim(ixv^s,mom(3)),&
809 wprim(ixvb^s,mom(3))-wprim(ixo^s,mom(3)),&
810 wprim(ixvc^s,mom(3))-wprim(ixvd^s,mom(3)),&
811 wprim(ixve^s,mom(3))-wprim(ixvc^s,mom(3))&
814 thetasm = maxval(cmaxc(ixo^s,1))
816 thetasm = (min(1d0, (thetasm-du)/(thetasm-min(dv,dw))))**aparam
820 br(ixc^s,:)=wrc(ixc^s,mag(:))+block%B0(ixc^s,:,ip1)
821 bl(ixc^s,:)=wlc(ixc^s,mag(:))+block%B0(ixc^s,:,ip1)
823 br(ixc^s,:)=wrc(ixc^s,mag(:))
824 bl(ixc^s,:)=wlc(ixc^s,mag(:))
829 ptr(ixc^s)=wrp(ixc^s,p_)+0.5d0*sum(br(ixc^s,:)**2,dim=ndim+1)
830 ptl(ixc^s)=wlp(ixc^s,p_)+0.5d0*sum(bl(ixc^s,:)**2,dim=ndim+1)
831 sur(ixc^s)=(sr(ixc^s,
index_v_mag)-vrc(ixc^s,ip1))*wrc(ixc^s,rho_)
832 sul(ixc^s)=(sl(ixc^s,
index_v_mag)-vlc(ixc^s,ip1))*wlc(ixc^s,rho_)
834 sm(ixc^s)=(sur(ixc^s)*vrc(ixc^s,ip1)-sul(ixc^s)*vlc(ixc^s,ip1)-&
835 thetasm*(ptr(ixc^s)-ptl(ixc^s)) )/(sur(ixc^s)-sul(ixc^s))
837 w1r(ixc^s,mom(ip1))=sm(ixc^s)
838 w1l(ixc^s,mom(ip1))=sm(ixc^s)
839 w2r(ixc^s,mom(ip1))=sm(ixc^s)
840 w2l(ixc^s,mom(ip1))=sm(ixc^s)
842 w1r(ixc^s,mag(ip1))=bx(ixc^s)
843 w1l(ixc^s,mag(ip1))=bx(ixc^s)
845 ptr(ixc^s)=wrp(ixc^s,p_)+0.5d0*sum(wrc(ixc^s,mag(:))**2,dim=ndim+1)
846 ptl(ixc^s)=wlp(ixc^s,p_)+0.5d0*sum(wlc(ixc^s,mag(:))**2,dim=ndim+1)
850 w1r(ixc^s,rho_)=sur(ixc^s)/(sr(ixc^s,
index_v_mag)-sm(ixc^s))
851 w1l(ixc^s,rho_)=sul(ixc^s)/(sl(ixc^s,
index_v_mag)-sm(ixc^s))
855 r1r(ixc^s)=sur(ixc^s)*(sr(ixc^s,
index_v_mag)-sm(ixc^s))-bx(ixc^s)**2
856 where(r1r(ixc^s)/=0.d0)
857 r1r(ixc^s)=1.d0/r1r(ixc^s)
859 r1l(ixc^s)=sul(ixc^s)*(sl(ixc^s,
index_v_mag)-sm(ixc^s))-bx(ixc^s)**2
860 where(r1l(ixc^s)/=0.d0)
861 r1l(ixc^s)=1.d0/r1l(ixc^s)
864 w1r(ixc^s,mom(ip2))=vrc(ixc^s,ip2)-bx(ixc^s)*br(ixc^s,ip2)*&
865 (sm(ixc^s)-vrc(ixc^s,ip1))*r1r(ixc^s)
866 w1l(ixc^s,mom(ip2))=vlc(ixc^s,ip2)-bx(ixc^s)*bl(ixc^s,ip2)*&
867 (sm(ixc^s)-vlc(ixc^s,ip1))*r1l(ixc^s)
869 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)
870 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)
875 w1r(ixc^s,mom(ip3))=vrc(ixc^s,ip3)-bx(ixc^s)*br(ixc^s,ip3)*&
876 (sm(ixc^s)-vrc(ixc^s,ip1))*r1r(ixc^s)
877 w1l(ixc^s,mom(ip3))=vlc(ixc^s,ip3)-bx(ixc^s)*bl(ixc^s,ip3)*&
878 (sm(ixc^s)-vlc(ixc^s,ip1))*r1l(ixc^s)
880 w1r(ixc^s,mag(ip3))=br(ixc^s,ip3)*w1r(ixc^s,mag(ip2))
881 w1l(ixc^s,mag(ip3))=bl(ixc^s,ip3)*w1l(ixc^s,mag(ip2))
884 w1r(ixc^s,mag(ip2))=br(ixc^s,ip2)*w1r(ixc^s,mag(ip2))
885 w1l(ixc^s,mag(ip2))=bl(ixc^s,ip2)*w1l(ixc^s,mag(ip2))
888 w1r(ixc^s,mag(:))=w1r(ixc^s,mag(:))-block%B0(ixc^s,:,ip1)
889 w1l(ixc^s,mag(:))=w1l(ixc^s,mag(:))-block%B0(ixc^s,:,ip1)
894 w1r(ixc^s,p_)=(sur(ixc^s)*ptl(ixc^s) - sul(ixc^s)*ptr(ixc^s) +&
895 phipres * sur(ixc^s)*sul(ixc^s)*(vrc(ixc^s,ip1)-vlc(ixc^s,ip1)))/&
896 (sur(ixc^s)-sul(ixc^s))
897 w1l(ixc^s,p_)=w1r(ixc^s,p_)
904 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)
905 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)
908 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)+&
909 w1r(ixc^s,p_)*sm(ixc^s)+bx(ixc^s)*(sum(vrc(ixc^s,:)*wrc(ixc^s,mag(:)),dim=ndim+1)-&
910 sum(w1r(ixc^s,mom(:))*w1r(ixc^s,mag(:)),dim=ndim+1)))/(sr(ixc^s,
index_v_mag)-sm(ixc^s))
911 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)+&
912 w1l(ixc^s,p_)*sm(ixc^s)+bx(ixc^s)*(sum(vlc(ixc^s,:)*wlc(ixc^s,mag(:)),dim=ndim+1)-&
913 sum(w1l(ixc^s,mom(:))*w1l(ixc^s,mag(:)),dim=ndim+1)))/(sl(ixc^s,
index_v_mag)-sm(ixc^s))
916 w1r(ixc^s,e_)=w1r(ixc^s,e_)+(sum(w1r(ixc^s,mag(:))*block%B0(ixc^s,:,ip1),dim=ndim+1)*sm(ixc^s)-&
917 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))
918 w1l(ixc^s,e_)=w1l(ixc^s,e_)+(sum(w1l(ixc^s,mag(:))*block%B0(ixc^s,:,ip1),dim=ndim+1)*sm(ixc^s)-&
919 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))
924 w2r(ixc^s,rho_)=w1r(ixc^s,rho_)
925 w2l(ixc^s,rho_)=w1l(ixc^s,rho_)
926 w2r(ixc^s,mag(ip1))=w1r(ixc^s,mag(ip1))
927 w2l(ixc^s,mag(ip1))=w1l(ixc^s,mag(ip1))
929 r1r(ixc^s)=sqrt(w1r(ixc^s,rho_))
930 r1l(ixc^s)=sqrt(w1l(ixc^s,rho_))
931 tmp(ixc^s)=1.d0/(r1r(ixc^s)+r1l(ixc^s))
932 signbx(ixc^s)=sign(1.d0,bx(ixc^s))
934 s1r(ixc^s)=sm(ixc^s)+abs(bx(ixc^s))/r1r(ixc^s)
935 s1l(ixc^s)=sm(ixc^s)-abs(bx(ixc^s))/r1l(ixc^s)
937 w2r(ixc^s,mom(ip2))=(r1l(ixc^s)*w1l(ixc^s,mom(ip2))+r1r(ixc^s)*w1r(ixc^s,mom(ip2))+&
938 (w1r(ixc^s,mag(ip2))-w1l(ixc^s,mag(ip2)))*signbx(ixc^s))*tmp(ixc^s)
939 w2l(ixc^s,mom(ip2))=w2r(ixc^s,mom(ip2))
941 w2r(ixc^s,mag(ip2))=(r1l(ixc^s)*w1r(ixc^s,mag(ip2))+r1r(ixc^s)*w1l(ixc^s,mag(ip2))+&
942 r1l(ixc^s)*r1r(ixc^s)*(w1r(ixc^s,mom(ip2))-w1l(ixc^s,mom(ip2)))*signbx(ixc^s))*tmp(ixc^s)
943 w2l(ixc^s,mag(ip2))=w2r(ixc^s,mag(ip2))
946 w2r(ixc^s,mom(ip3))=(r1l(ixc^s)*w1l(ixc^s,mom(ip3))+r1r(ixc^s)*w1r(ixc^s,mom(ip3))+&
947 (w1r(ixc^s,mag(ip3))-w1l(ixc^s,mag(ip3)))*signbx(ixc^s))*tmp(ixc^s)
948 w2l(ixc^s,mom(ip3))=w2r(ixc^s,mom(ip3))
950 w2r(ixc^s,mag(ip3))=(r1l(ixc^s)*w1r(ixc^s,mag(ip3))+r1r(ixc^s)*w1l(ixc^s,mag(ip3))+&
951 r1l(ixc^s)*r1r(ixc^s)*(w1r(ixc^s,mom(ip3))-w1l(ixc^s,mom(ip3)))*signbx(ixc^s))*tmp(ixc^s)
952 w2l(ixc^s,mag(ip3))=w2r(ixc^s,mag(ip3))
956 w2r(ixc^s,e_)=w1r(ixc^s,e_)+r1r(ixc^s)*(sum(w1r(ixc^s,mom(:))*w1r(ixc^s,mag(:)),dim=ndim+1)-&
957 sum(w2r(ixc^s,mom(:))*w2r(ixc^s,mag(:)),dim=ndim+1))*signbx(ixc^s)
958 w2l(ixc^s,e_)=w1l(ixc^s,e_)-r1l(ixc^s)*(sum(w1l(ixc^s,mom(:))*w1l(ixc^s,mag(:)),dim=ndim+1)-&
959 sum(w2l(ixc^s,mom(:))*w2l(ixc^s,mag(:)),dim=ndim+1))*signbx(ixc^s)
964 w1r(ixc^s,mom(idir))=w1r(ixc^s,mom(idir))*w1r(ixc^s,rho_)
965 w1l(ixc^s,mom(idir))=w1l(ixc^s,mom(idir))*w1l(ixc^s,rho_)
966 w2r(ixc^s,mom(idir))=w2r(ixc^s,mom(idir))*w2r(ixc^s,rho_)
967 w2l(ixc^s,mom(idir))=w2l(ixc^s,mom(idir))*w2l(ixc^s,rho_)
976 f1l(ixc^s,iw)=flc(ixc^s,iw)
977 f1r(ixc^s,iw)=f1l(ixc^s,iw)
978 f2l(ixc^s,iw)=f1l(ixc^s,iw)
979 f2r(ixc^s,iw)=f1l(ixc^s,iw)
984 f1r(ixc^s,iw)=f1l(ixc^s,iw)
985 f2l(ixc^s,iw)=f1l(ixc^s,iw)
986 f2r(ixc^s,iw)=f1l(ixc^s,iw)
988 f1l(ixc^s,iw)=flc(ixc^s,iw)+sl(ixc^s,
index_v_mag)*(w1l(ixc^s,iw)-wlc(ixc^s,iw))
989 f1r(ixc^s,iw)=frc(ixc^s,iw)+sr(ixc^s,
index_v_mag)*(w1r(ixc^s,iw)-wrc(ixc^s,iw))
990 f2l(ixc^s,iw)=f1l(ixc^s,iw)+s1l(ixc^s)*(w2l(ixc^s,iw)-w1l(ixc^s,iw))
991 f2r(ixc^s,iw)=f1r(ixc^s,iw)+s1r(ixc^s)*(w2r(ixc^s,iw)-w1r(ixc^s,iw))
996 {
do ix^db=ixcmin^db,ixcmax^db\}
999 else if(s1l(ix^d)>=0.d0)
then
1001 else if(sm(ix^d)>=0.d0)
then
1003 else if(s1r(ix^d)>=0.d0)
then
1013 end subroutine get_riemann_flux_hlld_mag2
1019 integer,
intent(in) :: ixI^L, ixO^L
1020 double precision,
intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
1021 double precision,
intent(out):: csound(ixI^S)
1022 double precision :: cfast2(ixI^S), AvMinCs2(ixI^S), b2(ixI^S), kmax
1023 double precision :: inv_rho(ixO^S), gamma_A2(ixO^S)
1024 integer :: rho_, p_, e_, eaux_, mom(1:ndir), mag(1:ndir)
1033 inv_rho=1.d0/w(ixo^s,rho_)
1038 b2(ixo^s) = sum((w(ixo^s, mag(:))+
block%B0(ixo^s,:,
b0i))**2, dim=ndim+1)
1040 b2(ixo^s) = sum(w(ixo^s, mag(:))**2, dim=ndim+1)
1045 avmincs2= w(ixo^s, mag(idims))+
block%B0(ixo^s,idims,
b0i)
1047 avmincs2= w(ixo^s, mag(idims))
1051 csound(ixo^s) = sum(w(ixo^s, mom(:))**2, dim=ndim+1)
1053 cfast2(ixo^s) = b2(ixo^s) * inv_rho+csound(ixo^s)
1054 avmincs2(ixo^s) = cfast2(ixo^s)**2-4.0d0*csound(ixo^s) &
1055 * avmincs2(ixo^s)**2 &
1058 where(avmincs2(ixo^s)<zero)
1059 avmincs2(ixo^s)=zero
1062 avmincs2(ixo^s)=sqrt(avmincs2(ixo^s))
1064 csound(ixo^s) = sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s)))
1072 subroutine reconstruct_lr(ixI^L,ixL^L,ixR^L,idims,w,wLC,wRC,wLp,wRp,x,dxdim)
1077 integer,
intent(in) :: ixi^
l, ixl^
l, ixr^
l, idims
1078 double precision,
intent(in) :: dxdim
1079 double precision,
dimension(ixI^S,1:nw) :: w
1081 double precision,
dimension(ixI^S,1:nw) :: wlc, wrc
1083 double precision,
dimension(ixI^S,1:nw) :: wlp, wrp
1084 double precision,
dimension(ixI^S,1:ndim) :: x
1086 integer :: jxr^
l, ixc^
l, jxc^
l, iw
1087 double precision :: ldw(ixi^s), rdw(ixi^s), dwc(ixi^s)
1088 double precision :: a2max
1092 call mp5limiter(ixi^
l,ixl^
l,idims,w,wlp,wrp)
1094 call weno3limiter(ixi^
l,ixl^
l,idims,dxdim,w,wlp,wrp,1)
1096 call weno3limiter(ixi^
l,ixl^
l,idims,dxdim,w,wlp,wrp,2)
1098 call weno5limiter(ixi^
l,ixl^
l,idims,dxdim,w,wlp,wrp,1)
1100 call weno5nmlimiter(ixi^
l,ixl^
l,idims,dxdim,w,wlp,wrp,1)
1102 call weno5limiter(ixi^
l,ixl^
l,idims,dxdim,w,wlp,wrp,2)
1104 call weno5nmlimiter(ixi^
l,ixl^
l,idims,dxdim,w,wlp,wrp,2)
1106 call weno5limiter(ixi^
l,ixl^
l,idims,dxdim,w,wlp,wrp,3)
1108 call weno5nmlimiter(ixi^
l,ixl^
l,idims,dxdim,w,wlp,wrp,3)
1110 call weno5cu6limiter(ixi^
l,ixl^
l,idims,w,wlp,wrp)
1112 call teno5adlimiter(ixi^
l,ixl^
l,idims,dxdim,w,wlp,wrp)
1114 call weno7limiter(ixi^
l,ixl^
l,idims,w,wlp,wrp,1)
1116 call weno7limiter(ixi^
l,ixl^
l,idims,w,wlp,wrp,2)
1118 call venklimiter(ixi^
l,ixl^
l,idims,dxdim,w,wlp,wrp)
1124 call ppmlimiter(ixi^
l,
ixm^
ll,idims,w,w,wlp,wrp)
1130 jxr^
l=ixr^
l+
kr(idims,^
d);
1131 ixcmax^
d=jxrmax^
d; ixcmin^
d=ixlmin^
d-
kr(idims,^
d);
1132 jxc^
l=ixc^
l+
kr(idims,^
d);
1135 w(ixcmin^
d:jxcmax^
d,iw)=dlog10(w(ixcmin^
d:jxcmax^
d,iw))
1136 wlp(ixl^s,iw)=dlog10(wlp(ixl^s,iw))
1137 wrp(ixr^s,iw)=dlog10(wrp(ixr^s,iw))
1140 dwc(ixc^s)=w(jxc^s,iw)-w(ixc^s,iw)
1156 call mpistop(
"idims is wrong in mod_limiter")
1162 wlp(ixl^s,iw)=wlp(ixl^s,iw)+half*ldw(ixl^s)
1163 wrp(ixr^s,iw)=wrp(ixr^s,iw)-half*rdw(jxr^s)
1166 w(ixcmin^
d:jxcmax^
d,iw)=10.0d0**w(ixcmin^
d:jxcmax^
d,iw)
1167 wlp(ixl^s,iw)=10.0d0**wlp(ixl^s,iw)
1168 wrp(ixr^s,iw)=10.0d0**wrp(ixr^s,iw)
1177 wlc(ixl^s,1:nw)=wlp(ixl^s,1:nw)
1178 wrc(ixr^s,1:nw)=wrp(ixr^s,1:nw)
subroutine mpistop(message)
Exit MPI-AMRVAC with an error message.
subroutine get_hlld2_modif_c(w, x, ixIL, ixOL, csound)
Calculate fast magnetosonic wave speed.
subroutine get_riemann_flux_tvdmu()
Module with finite volume methods for fluxes.
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 finite_volume(method, qdt, ixIL, ixOL, idimsLIM, qtC, sCT, qt, snew, sold, fC, fE, dxs, x)
finite volume method
subroutine, public hancock(qdt, 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 angmomfix
Enable to strictly conserve the angular momentum (works both in cylindrical and spherical coordinates...
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 need_global_a2max
global value for schmid scheme
integer ixm
the mesh range (within a block with 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
logical phys_solve_eaux
Solve internal energy and total energy equations.
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_energy_synchro), pointer phys_energy_synchro
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
procedure(sub_angmomfix), pointer phys_angmomfix
logical phys_energy
Solve energy equation or not.
Module for handling split source terms (split from the fluxes)
subroutine, public addsource2(qdt, ixIL, ixOL, iwLIM, qtC, wCT, qt, w, x, qsourcesplit, src_active, wCTprim)
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 iw_eaux
Index of the internal energy density.
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.