User: tbode Date: 2011/04/28 11:47 AM
Modified: /branches/divclean/src/ GRHydro_CalcUpdate.F90, GRHydro_FluxM.F90, GRHydro_HLLEM.F90, GRHydro_Prim2ConM.F90, GRHydro_Reconstruct.F90, GRHydro_SourceM.F90
Log: Josh's fixes.
File Changes:
Directory: /branches/divclean/src/ ==================================
File [modified]: GRHydro_CalcUpdate.F90 Delta lines: +7 -0 =================================================================== --- branches/divclean/src/GRHydro_CalcUpdate.F90 2011-04-28 01:28:55 UTC (rev 242) +++ branches/divclean/src/GRHydro_CalcUpdate.F90 2011-04-28 16:47:14 UTC (rev 243) @@ -102,6 +102,13 @@ Bvec(i+xoffset,j+yoffset,k+zoffset,flux_direction)) divB(i,j,k) = divB(i,j,k) + ( alp_l * Bvec_l - alp_r * Bvec_r ) * idx endif + +!!$ if(i.eq.75.and.j.eq.5.and.k.eq.5)write(6,*)'update:', & +!!$ Bconsrhs(i,j,k,1),Bconsxflux(i,j,k), & +!!$ Bconsxflux(i-xoffset,j-yoffset,k-zoffset), & +!!$ Bconsrhs(i,j,k,3),Bconszflux(i,j,k), & +!!$ Bconszflux(i-xoffset,j-yoffset,k-zoffset) + endif
if (evolve_tracer .ne. 0) then
File [modified]: GRHydro_FluxM.F90 Delta lines: +7 -10 =================================================================== --- branches/divclean/src/GRHydro_FluxM.F90 2011-04-28 01:28:55 UTC (rev 242) +++ branches/divclean/src/GRHydro_FluxM.F90 2011-04-28 16:47:14 UTC (rev 243) @@ -23,13 +23,10 @@ CCTK_REAL :: vxt,vyt,vzt,bsubx,bsuby,bsubz,ab0,w CCTK_REAL :: det,alp,beta,pressstar CCTK_REAL :: velm - CCTK_REAL :: sdet,psipstar,psiBx,psiBy,psiBz + CCTK_REAL :: sdet,psipstar
sdet=sqrt(det) psipstar=pressstar*sdet - psiBx=Bx*sdet - psiBy=By*sdet - psiBz=Bz*sdet
!!$ We actually need all three values of vtilde = v^i - beta^i/alp, as well as velm=vxt+beta/alp @@ -38,19 +35,19 @@ !!$ In the notation of Anton et al.: [psi^6 D] vtilde^i densf = dens * vxt
- sxf = sx*vxt+psipstar-bsubx*psiBx/w + sxf = sx*vxt+psipstar-bsubx*Bx/w
- syf = sy*vxt-bsuby*psiBx/w + syf = sy*vxt-bsuby*Bx/w
- szf = sz*vxt-bsubz*psiBx/w + szf = sz*vxt-bsubz*Bx/w
!!$ [psi^6 tau] vtilde^i +p* v^i - alp b^0 B^i/w - tauf = tau*vxt + psipstar*velm - ab0*psiBx/w + tauf = tau*vxt + psipstar*velm - ab0*Bx/w
!!$ [psi^6 (B^k vtilde^i - B^i vtilde^k)] bxf = 0.0 - byf = psiBy * vxt - psiBx*vyt - bzf = psiBz * vxt - psiBx*vzt + byf = By * vxt - Bx*vyt + bzf = Bz * vxt - Bx*vzt
end subroutine num_x_fluxM
File [modified]: GRHydro_HLLEM.F90 Delta lines: +88 -78 =================================================================== --- branches/divclean/src/GRHydro_HLLEM.F90 2011-04-28 01:28:55 UTC (rev 242) +++ branches/divclean/src/GRHydro_HLLEM.F90 2011-04-28 16:47:14 UTC (rev 243) @@ -40,10 +40,9 @@ DECLARE_CCTK_PARAMETERS
integer :: i, j, k, m - CCTK_REAL, dimension(8) :: fplus,fminus,f1,qdiff - CCTK_REAL, dimension(7) :: prim_p, prim_m - CCTK_REAL, dimension(5) :: lamminus,lamplus,cons_p,cons_m - CCTK_REAL, dimension(3) :: mag_p,mag_m + CCTK_REAL, dimension(8) :: cons_p,cons_m,fplus,fminus,f1,qdiff + CCTK_REAL, dimension(10) :: prim_p, prim_m + CCTK_REAL, dimension(5) :: lamminus,lamplus CCTK_REAL :: charmin, charmax, charpm,avg_alp,avg_det, sdet CCTK_REAL :: gxxh, gxyh, gxzh, gyyh, gyzh, gzzh, uxxh, uxyh, & uxzh, uyyh, uyzh, uzzh, avg_beta, usendh, alp_l, alp_r, & @@ -83,8 +82,6 @@ lamplus = 0.d0 cons_p = 0.d0 cons_m = 0.d0 - mag_p=0.d0 - mag_m=0.d0 fplus = 0.d0 fminus = 0.d0 qdiff = 0.d0 @@ -101,20 +98,20 @@ cons_p(3) = syplus(i,j,k) cons_p(4) = szplus(i,j,k) cons_p(5) = tauplus(i,j,k) - - mag_p(1) = Bvecxplus(i,j,k) - mag_p(2) = Bvecyplus(i,j,k) - mag_p(3) = Bveczplus(i,j,k) + cons_p(6) = Bconsxplus(i,j,k) + cons_p(7) = Bconsyplus(i,j,k) + cons_p(8) = Bconszplus(i,j,k)
cons_m(1) = densminus(i+xoffset,j+yoffset,k+zoffset) cons_m(2) = sxminus(i+xoffset,j+yoffset,k+zoffset) cons_m(3) = syminus(i+xoffset,j+yoffset,k+zoffset) cons_m(4) = szminus(i+xoffset,j+yoffset,k+zoffset) cons_m(5) = tauminus(i+xoffset,j+yoffset,k+zoffset) + cons_m(6) = Bconsxminus(i+xoffset,j+yoffset,k+zoffset) + cons_m(7) = Bconsyminus(i+xoffset,j+yoffset,k+zoffset) + cons_m(8) = Bconszminus(i+xoffset,j+yoffset,k+zoffset)
- mag_m(1) = Bvecxminus(i+xoffset,j+yoffset,k+zoffset) - mag_m(2) = Bvecyminus(i+xoffset,j+yoffset,k+zoffset) - mag_m(3) = Bveczminus(i+xoffset,j+yoffset,k+zoffset) +!!$ if(i.eq.75.and.j.eq.5.and.k.eq.5)write(6,*)'HLLEM0:',cons_m(6),cons_p(6),cons_m(7),cons_p(7),cons_m(8),cons_p(8)
prim_p(1) = rhoplus(i,j,k) prim_p(2) = velxplus(i,j,k) @@ -123,6 +120,9 @@ prim_p(5) = epsplus(i,j,k) prim_p(6) = pressplus(i,j,k) prim_p(7) = w_lorentzplus(i,j,k) + prim_p(8) = Bvecxplus(i,j,k) + prim_p(9) = Bvecyplus(i,j,k) + prim_p(10) = Bveczplus(i,j,k)
prim_m(1) = rhominus(i+xoffset,j+yoffset,k+zoffset) prim_m(2) = velxminus(i+xoffset,j+yoffset,k+zoffset) @@ -131,6 +131,9 @@ prim_m(5) = epsminus(i+xoffset,j+yoffset,k+zoffset) prim_m(6) = pressminus(i+xoffset,j+yoffset,k+zoffset) prim_m(7) = w_lorentzminus(i+xoffset,j+yoffset,k+zoffset) + prim_m(8) = Bvecxminus(i+xoffset,j+yoffset,k+zoffset) + prim_m(9) = Bvecyminus(i+xoffset,j+yoffset,k+zoffset) + prim_m(10)= Bveczminus(i+xoffset,j+yoffset,k+zoffset)
if(clean_divergence.ne.0) then psidcp = psidcplus(i,j,k) @@ -191,11 +194,11 @@ vztm = prim_m(4)-avg_betaz/avg_alp
call calc_vlow_blow(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh, & - prim_p(2),prim_p(3),prim_p(4),mag_p(1),mag_p(2),mag_p(3), & + prim_p(2),prim_p(3),prim_p(4),prim_p(8),prim_p(9),prim_p(10), & velxlowp,velylowp,velzlowp,Bvecxlowp,Bvecylowp,Bveczlowp, & bdotvp,b2p,v2p,wp,bxlowp,bylowp,bzlowp) call calc_vlow_blow(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh, & - prim_m(2),prim_m(3),prim_m(4),mag_m(1),mag_m(2),mag_m(3), & + prim_m(2),prim_m(3),prim_m(4),prim_m(8),prim_m(9),prim_m(10), & velxlowm,velylowm,velzlowm,Bvecxlowm,Bvecylowm,Bveczlowm, & bdotvm,b2m,v2m,wm,bxlowm,bylowm,bzlowm)
@@ -224,41 +227,41 @@
if (flux_direction == 1) then call num_x_fluxM(cons_p(1),cons_p(2),cons_p(3),cons_p(4),cons_p(5),& - mag_p(1),mag_p(2),mag_p(3),& + cons_p(6),cons_p(7),cons_p(8),& f1(1),f1(2),f1(3),f1(4),f1(5),f1(6),f1(7),f1(8),& vxtp,vytp,vztp,pressstarp,bxlowp,bylowp,bzlowp,ab0p,wp, & avg_det,avg_alp,avg_beta) if(clean_divergence.ne.0) then - f1(6)=f1(6) + avg_alp*sdet*uxxh*psidcp - sdet*mag_p(1)*avg_betax - f1(7)=f1(7) + avg_alp*sdet*uxyh*psidcp - sdet*mag_p(1)*avg_betay - f1(8)=f1(8) + avg_alp*sdet*uxzh*psidcp - sdet*mag_p(1)*avg_betaz - psidcf = avg_alp*mag_p(1)-psidcp*avg_betax + f1(6)=f1(6) + avg_alp*sdet*uxxh*psidcp - cons_p(6)*avg_betax + f1(7)=f1(7) + avg_alp*sdet*uxyh*psidcp - cons_p(6)*avg_betay + f1(8)=f1(8) + avg_alp*sdet*uxzh*psidcp - cons_p(6)*avg_betaz + psidcf = avg_alp*cons_p(6)/sdet-psidcp*avg_betax endif
else if (flux_direction == 2) then call num_x_fluxM(cons_p(1),cons_p(3),cons_p(4),cons_p(2),cons_p(5),& - mag_p(2),mag_p(3),mag_p(1),& + cons_p(7),cons_p(8),cons_p(6),& f1(1),f1(3),f1(4),f1(2),f1(5),f1(7),f1(8),f1(6),& vytp,vztp,vxtp,pressstarp,bylowp,bzlowp,bxlowp,ab0p,wp, & avg_det,avg_alp,avg_beta) if(clean_divergence.ne.0) then - f1(6)=f1(6) + avg_alp*sdet*uxyh*psidcp - sdet*mag_p(2)*avg_betax - f1(7)=f1(7) + avg_alp*sdet*uyyh*psidcp - sdet*mag_p(2)*avg_betay - f1(8)=f1(8) + avg_alp*sdet*uyzh*psidcp - sdet*mag_p(2)*avg_betaz - psidcf = avg_alp*mag_p(2)-psidcp*avg_betay + f1(6)=f1(6) + avg_alp*sdet*uxyh*psidcp - cons_p(7)*avg_betax + f1(7)=f1(7) + avg_alp*sdet*uyyh*psidcp - cons_p(7)*avg_betay + f1(8)=f1(8) + avg_alp*sdet*uyzh*psidcp - cons_p(7)*avg_betaz + psidcf = avg_alp*cons_p(7)/sdet-psidcp*avg_betay endif
else if (flux_direction == 3) then call num_x_fluxM(cons_p(1),cons_p(4),cons_p(2),cons_p(3),cons_p(5),& - mag_p(3),mag_p(1),mag_p(2),& + cons_p(8),cons_p(6),cons_p(7),& f1(1),f1(4),f1(2),f1(3),f1(5),f1(8),f1(6),f1(7), & vztp,vxtp,vytp,pressstarp,bzlowp,bxlowp,bylowp,ab0p,wp, & avg_det,avg_alp,avg_beta) if(clean_divergence.ne.0) then - f1(6)=f1(6) + avg_alp*sdet*uxzh*psidcp - sdet*mag_p(3)*avg_betax - f1(7)=f1(7) + avg_alp*sdet*uyzh*psidcp - sdet*mag_p(3)*avg_betay - f1(8)=f1(8) + avg_alp*sdet*uzzh*psidcp - sdet*mag_p(3)*avg_betaz - psidcf = avg_alp*mag_p(3)-psidcp*avg_betaz + f1(6)=f1(6) + avg_alp*sdet*uxzh*psidcp - cons_p(8)*avg_betax + f1(7)=f1(7) + avg_alp*sdet*uyzh*psidcp - cons_p(8)*avg_betay + f1(8)=f1(8) + avg_alp*sdet*uzzh*psidcp - cons_p(8)*avg_betaz + psidcf = avg_alp*cons_p(8)/sdet-psidcp*avg_betaz endif
else @@ -284,11 +287,10 @@ qdiff(3) = cons_m(3) - cons_p(3) qdiff(4) = cons_m(4) - cons_p(4) qdiff(5) = cons_m(5) - cons_p(5) + qdiff(6) = cons_m(6) - cons_p(6) + qdiff(7) = cons_m(7) - cons_p(7) + qdiff(8) = cons_m(8) - cons_p(8)
- qdiff(6) = mag_m(1) - mag_p(1) - qdiff(7) = mag_m(2) - mag_p(2) - qdiff(8) = mag_m(3) - mag_p(3) - if (clean_divergence.ne.0) then psidcdiff = psidcm - psidcp endif @@ -298,107 +300,108 @@ if (flux_direction == 1) then call eigenvaluesM(GRHydro_eos_handle,& prim_m(1),prim_m(2),prim_m(3),prim_m(4),prim_m(5),prim_m(7), & - mag_m(1),mag_m(2),mag_m(3),& + prim_m(8),prim_m(9),prim_m(10),& lamminus,gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,& usendh,avg_alp,avg_beta) call eigenvaluesM(GRHydro_eos_handle, & prim_p(1),prim_p(2),prim_p(3),prim_p(4),prim_p(5),prim_p(7), & - mag_p(1),mag_p(2),mag_p(3),& + prim_p(8),prim_p(9),prim_p(10),& lamplus,gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,& usendh,avg_alp,avg_beta) call num_x_fluxM(cons_p(1),cons_p(2),cons_p(3),cons_p(4),cons_p(5),& - mag_p(1),mag_p(2),mag_p(3),& + cons_p(6),cons_p(7),cons_p(8),& fplus(1),fplus(2),fplus(3),fplus(4),fplus(5),fplus(6),fplus(7),fplus(8),& vxtp,vytp,vztp,pressstarp,bxlowp,bylowp,bzlowp,ab0p,wp, & avg_det,avg_alp,avg_beta) call num_x_fluxM(cons_m(1),cons_m(2),cons_m(3),cons_m(4),cons_m(5),& - mag_m(1),mag_m(2),mag_m(3),& + cons_m(6),cons_m(7),cons_m(8),& fminus(1),fminus(2),fminus(3),fminus(4),fminus(5),& fminus(6),fminus(7),fminus(8),& vxtm,vytm,vztm,pressstarm,bxlowm,bylowm,bzlowm,ab0m,wm, & avg_det,avg_alp,avg_beta)
if(clean_divergence.ne.0) then - fminus(6)=fminus(6) + avg_alp*sdet*uxxh*psidcm - sdet*mag_m(1)*avg_betax - fminus(7)=fminus(7) + avg_alp*sdet*uxyh*psidcm - sdet*mag_m(1)*avg_betay - fminus(8)=fminus(8) + avg_alp*sdet*uxzh*psidcm - sdet*mag_m(1)*avg_betaz - fplus(6)=fplus(6) + avg_alp*sdet*uxxh*psidcp - sdet*mag_p(1)*avg_betax - fplus(7)=fplus(7) + avg_alp*sdet*uxyh*psidcp - sdet*mag_p(1)*avg_betay - fplus(8)=fplus(8) + avg_alp*sdet*uxzh*psidcp - sdet*mag_p(1)*avg_betaz - psidcfp = avg_alp*mag_p(1)-avg_betax*psidcp - psidcfm = avg_alp*mag_m(1)-avg_betax*psidcm + fminus(6)=fminus(6) + avg_alp*sdet*uxxh*psidcm - cons_m(6)*avg_betax + fminus(7)=fminus(7) + avg_alp*sdet*uxyh*psidcm - cons_m(6)*avg_betay + fminus(8)=fminus(8) + avg_alp*sdet*uxzh*psidcm - cons_m(6)*avg_betaz + fplus(6)=fplus(6) + avg_alp*sdet*uxxh*psidcp - cons_p(6)*avg_betax + fplus(7)=fplus(7) + avg_alp*sdet*uxyh*psidcp - cons_p(6)*avg_betay + fplus(8)=fplus(8) + avg_alp*sdet*uxzh*psidcp - cons_p(6)*avg_betaz + psidcfp = avg_alp*cons_p(6)/sdet-avg_betax*psidcp + psidcfm = avg_alp*cons_m(6)/sdet-avg_betax*psidcm endif
else if (flux_direction == 2) then call eigenvaluesM(GRHydro_eos_handle,& prim_m(1),prim_m(3),prim_m(4),prim_m(2),prim_m(5),prim_m(7), & - mag_m(2),mag_m(3),mag_m(1),& + prim_m(9),prim_m(10),prim_m(8),& lamminus,gyyh,gyzh,gxyh,gzzh,gxzh,gxxh,& usendh,avg_alp,avg_beta) call eigenvaluesM(GRHydro_eos_handle, & prim_p(1),prim_p(3),prim_p(4),prim_p(2),prim_p(5),prim_p(7), & - mag_p(2),mag_p(3),mag_p(1),& + prim_p(9),prim_p(10),prim_p(8),& lamplus,gyyh,gyzh,gxyh,gzzh,gxzh,gxxh,& usendh,avg_alp,avg_beta) call num_x_fluxM(cons_p(1),cons_p(3),cons_p(4),cons_p(2),cons_p(5),& - mag_p(2),mag_p(3),mag_p(1),& + cons_p(7),cons_p(8),cons_p(6),& fplus(1),fplus(3),fplus(4),fplus(2),fplus(5),fplus(7),fplus(8),fplus(6),& vytp,vztp,vxtp,pressstarp,bylowp,bzlowp,bxlowp,ab0p,wp, & avg_det,avg_alp,avg_beta) call num_x_fluxM(cons_m(1),cons_m(3),cons_m(4),cons_m(2),cons_m(5),& - mag_m(2),mag_m(3),mag_m(1),& + cons_m(7),cons_m(8),cons_m(6),& fminus(1),fminus(3),fminus(4),fminus(2),fminus(5),& fminus(7),fminus(8),fminus(6),& vytm,vztm,vxtm,pressstarm,bylowm,bzlowm,bxlowm,ab0m,wm, & avg_det,avg_alp,avg_beta) if(clean_divergence.ne.0) then - fminus(6)=fminus(6) + avg_alp*sdet*uxyh*psidcm - sdet*mag_m(2)*avg_betax - fminus(7)=fminus(7) + avg_alp*sdet*uyyh*psidcm - sdet*mag_m(2)*avg_betay - fminus(8)=fminus(8) + avg_alp*sdet*uyzh*psidcm - sdet*mag_m(2)*avg_betaz - fplus(6)=fplus(6) + avg_alp*sdet*uxyh*psidcp - sdet*mag_p(2)*avg_betax - fplus(7)=fplus(7) + avg_alp*sdet*uyyh*psidcp - sdet*mag_p(2)*avg_betay - fplus(8)=fplus(8) + avg_alp*sdet*uyzh*psidcp - sdet*mag_p(2)*avg_betaz - psidcfp = avg_alp*mag_p(2)-avg_betay*psidcp - psidcfm = avg_alp*mag_m(2)-avg_betay*psidcm + fminus(6)=fminus(6) + avg_alp*sdet*uxyh*psidcm - cons_m(7)*avg_betax + fminus(7)=fminus(7) + avg_alp*sdet*uyyh*psidcm - cons_m(7)*avg_betay + fminus(8)=fminus(8) + avg_alp*sdet*uyzh*psidcm - cons_m(7)*avg_betaz + fplus(6)=fplus(6) + avg_alp*sdet*uxyh*psidcp - cons_p(7)*avg_betax + fplus(7)=fplus(7) + avg_alp*sdet*uyyh*psidcp - cons_p(7)*avg_betay + fplus(8)=fplus(8) + avg_alp*sdet*uyzh*psidcp - cons_p(7)*avg_betaz + psidcfp = avg_alp*cons_p(7)/sdet-avg_betay*psidcp + psidcfm = avg_alp*cons_m(7)/sdet-avg_betay*psidcm endif
else if (flux_direction == 3) then call eigenvaluesM(GRHydro_eos_handle,& prim_m(1),prim_m(4),prim_m(2),prim_m(3),prim_m(5),prim_m(7), & - mag_m(3),mag_m(1),mag_m(2),& + prim_m(10),prim_m(8),prim_m(9),& lamminus,gzzh,gxzh,gyzh,gxxh,gxyh,gyyh,& usendh,avg_alp,avg_beta) call eigenvaluesM(GRHydro_eos_handle,& prim_p(1),prim_p(4),prim_p(2),prim_p(3),prim_p(5),prim_p(7), & - mag_p(3),mag_p(1),mag_p(2),& + prim_p(10),prim_p(8),prim_p(9),& lamplus,gzzh,gxzh,gyzh,gxxh,gxyh,gyyh,& usendh,avg_alp,avg_beta) call num_x_fluxM(cons_p(1),cons_p(4),cons_p(2),cons_p(3),cons_p(5),& - mag_p(3),mag_p(1),mag_p(2),& + cons_p(8),cons_p(6),cons_p(7),& fplus(1),fplus(4),fplus(2),fplus(3),fplus(5),fplus(8),fplus(6),fplus(7), & vztp,vxtp,vytp,pressstarp,bzlowp,bxlowp,bylowp,ab0p,wp, & avg_det,avg_alp,avg_beta) call num_x_fluxM(cons_m(1),cons_m(4),cons_m(2),cons_m(3),cons_m(5),& - mag_m(3),mag_m(1),mag_m(2),& + cons_m(8),cons_m(6),cons_m(7),& fminus(1),fminus(4),fminus(2),fminus(3),fminus(5), & fminus(8),fminus(6),fminus(7), & vztm,vxtm,vytm,pressstarm,bzlowm,bxlowm,bylowm,ab0m,wm, & avg_det,avg_alp,avg_beta) if(clean_divergence.ne.0) then - fminus(6)=fminus(6) + avg_alp*sdet*uxzh*psidcm - sdet*mag_m(3)*avg_betax - fminus(7)=fminus(7) + avg_alp*sdet*uyzh*psidcm - sdet*mag_m(3)*avg_betay - fminus(8)=fminus(8) + avg_alp*sdet*uzzh*psidcm - sdet*mag_m(3)*avg_betaz - fplus(6)=fplus(6) + avg_alp*sdet*uxzh*psidcp - sdet*mag_p(3)*avg_betax - fplus(7)=fplus(7) + avg_alp*sdet*uyzh*psidcp - sdet*mag_p(3)*avg_betay - fplus(8)=fplus(8) + avg_alp*sdet*uzzh*psidcp - sdet*mag_p(3)*avg_betaz - psidcfp = avg_alp*mag_p(3)-avg_betaz*psidcp - psidcfm = avg_alp*mag_m(3)-avg_betaz*psidcm + fminus(6)=fminus(6) + avg_alp*sdet*uxzh*psidcm - cons_m(8)*avg_betax + fminus(7)=fminus(7) + avg_alp*sdet*uyzh*psidcm - cons_m(8)*avg_betay + fminus(8)=fminus(8) + avg_alp*sdet*uzzh*psidcm - cons_m(8)*avg_betaz + fplus(6)=fplus(6) + avg_alp*sdet*uxzh*psidcp - cons_p(8)*avg_betax + fplus(7)=fplus(7) + avg_alp*sdet*uyzh*psidcp - cons_p(8)*avg_betay + fplus(8)=fplus(8) + avg_alp*sdet*uzzh*psidcp - cons_p(8)*avg_betaz + psidcfp = avg_alp*cons_p(8)/sdet-avg_betaz*psidcp + psidcfm = avg_alp*cons_m(8)/sdet-avg_betaz*psidcm endif
else call CCTK_WARN(0, "Flux direction not x,y,z") end if - + + !!$ Find minimum and maximum wavespeeds
charmin = min(0.d0, lamplus(1), lamplus(2), lamplus(3), & @@ -410,18 +413,23 @@ lamminus(4),lamminus(5))
charpm = charmax - charmin + +!!$ if(i.eq.75.and.j.eq.5.and.k.eq.5)write(6,*)'HLLEM:',fplus(6),fminus(6),fplus(8),fminus(8),charmax,charmin,cons_p(6),cons_m(6)
!!$ Calculate flux by standard formula
- do m = 1,8 - - qdiff(m) = cons_m(m) - cons_p(m) - - f1(m) = (charmax * fplus(m) - charmin * fminus(m) + & - charmax * charmin * qdiff(m)) / charpm - + do m = 1,8 + + qdiff(m) = cons_m(m) - cons_p(m) + + f1(m) = (charmax * fplus(m) - charmin * fminus(m) + & + charmax * charmin * qdiff(m)) / charpm + end do
+!!$ if(i.eq.75.and.j.eq.5.and.k.eq.5)write(6,*)'HLLEM2:',fplus(6),fminus(6),fplus(8),fminus(8),charmax,charmin,qdiff(6),qdiff(8),f1(6),f1(8) + + if(clean_divergence.ne.0) then psidcdiff = psidcm - psidcp
@@ -448,6 +456,8 @@ Bconsyflux(i, j, k) = f1(7) Bconszflux(i, j, k) = f1(8)
+!!$ if(i.eq.75.and.j.eq.5.and.k.eq.5)write(6,*)'HLLEM2:',fplus(6),fminus(6),fplus(8),fminus(8),f1(6),f1(8) + if(clean_divergence.ne.0) then psidcflux(i,j,k) = psidcf endif
File [modified]: GRHydro_Prim2ConM.F90 Delta lines: +2 -2 =================================================================== --- branches/divclean/src/GRHydro_Prim2ConM.F90 2011-04-28 01:28:55 UTC (rev 242) +++ branches/divclean/src/GRHydro_Prim2ConM.F90 2011-04-28 16:47:14 UTC (rev 243) @@ -88,10 +88,10 @@ call prim2conM(GRHydro_eos_handle, gxxr,gxyr,gxzr,gyyr,gyzr,gzzr, & avg_detr, densplus(i,j,k),sxplus(i,j,k),& syplus(i,j,k),szplus(i,j ,k),tauplus(i,j,k),& - Bvecxplus(i,j,k),Bvecyplus(i,j,k),Bveczplus(i,j,k), & + Bconsxplus(i,j,k),Bconsyplus(i,j,k),Bconszplus(i,j,k), & rhoplus(i,j,k),velxplus(i,j,k),velyplus(i,j,k),& velzplus(i,j,k),epsplus(i,j,k),pressplus(i,j,k),& - Bconsxminus(i,j,k), Bconsyminus(i,j,k), Bconszminus(i,j,k), & + Bvecxminus(i,j,k), Bvecyminus(i,j,k), Bveczminus(i,j,k), & w_lorentzplus(i,j,k))
end do
File [modified]: GRHydro_Reconstruct.F90 Delta lines: +6 -0 =================================================================== --- branches/divclean/src/GRHydro_Reconstruct.F90 2011-04-28 01:28:55 UTC (rev 242) +++ branches/divclean/src/GRHydro_Reconstruct.F90 2011-04-28 16:47:14 UTC (rev 243) @@ -219,6 +219,9 @@ call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & Bvec(:,:,:,3), Bveczplus, Bveczminus, trivial_rp, hydro_excision_mask) endif + +!!$ write(6,*)'recon0:',Bvecxplus(75,5,5),Bvecxminus(75,5,5) + else if (CCTK_EQUALS(recon_vars,"conservative")) then call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & dens, densplus, densminus, trivial_rp, hydro_excision_mask) @@ -920,6 +923,9 @@ CCTK_EQUALS(recon_method,"ppm")) then if(evolve_mhd.ne.0) then call primitive2conservativeM(CCTK_PASS_FTOF) + +!!$ write(6,*)'recon1:',Bvecxplus(75,5,5),Bvecxminus(75,5,5),Bconsxplus(75,5,5),Bconsxminus(75,5,5) + else call primitive2conservative(CCTK_PASS_FTOF) endif
File [modified]: GRHydro_SourceM.F90 Delta lines: +7 -0 =================================================================== --- branches/divclean/src/GRHydro_SourceM.F90 2011-04-28 01:28:55 UTC (rev 242) +++ branches/divclean/src/GRHydro_SourceM.F90 2011-04-28 16:47:14 UTC (rev 243) @@ -442,6 +442,10 @@ srhs(i,j,k,2) = alp(i,j,k)*sqrtdet*sy_source srhs(i,j,k,3) = alp(i,j,k)*sqrtdet*sz_source taurhs(i,j,k) = alp(i,j,k)*sqrtdet*tau_source + +!!$ if(i.eq.75.and.j.eq.5.and.k.eq.5)write(6,*)'source:',i,j,k, & +!!$ Bconsrhs(i,j,k,1),Bconsrhs(i,j,k,3) + if(clean_divergence.ne.0) then psidcrhs(i,j,k) = -1.d0 * ( ch_dc*ch_dc/cp_dc/cp_dc*alp(i,j,k) + & dx_betax + dy_betay + dz_betaz ) * psidc(i,j,k) + & @@ -500,6 +504,9 @@ psidc(i,j,k)*alp(i,j,k)*sqrtdet*( uxz*gdg_x + uyz*gdg_y + & uzz*gdg_z )
+!!$ if(i.eq.5.and.j.eq.5.and.k.eq.5)write(6,*)'source:',i,j,k, & +!!$ bvcx_source,bvcz_source + Bconsrhsx(i,j,k) = bvcx_source Bconsrhsy(i,j,k) = bvcy_source Bconsrhsz(i,j,k) = bvcz_source