[Commits] [svn:einsteintoolkit] GRHydro/branches/divclean/src/ (Rev. 243)
tanja.bode at physics.gatech.edu
tanja.bode at physics.gatech.edu
Thu Apr 28 11:47:15 CDT 2011
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
More information about the Commits
mailing list