User: rhaas
Date: 2014/04/15 02:48 PM
Modified:
/trunk/src/
GRHydro_HLLEM.F90
Log:
GRHydro: rewrite HLLEM without if(flux_direction) in inner loop
this comes at the cost of turning these into array indices, so we trade
if statements in the inner loop for array indices that are unknown at
compile time.
File Changes:
Directory: /trunk/src/
======================
File [modified]: GRHydro_HLLEM.F90
Delta lines: +139 -296
===================================================================
--- trunk/src/GRHydro_HLLEM.F90 2014-01-14 04:04:28 UTC (rev 598)
+++ trunk/src/GRHydro_HLLEM.F90 2014-04-15 19:48:05 UTC (rev 599)
@@ -45,13 +45,15 @@
CCTK_REAL, dimension(10) :: prim_p, prim_m
CCTK_REAL, dimension(5) :: lamminus,lamplus
CCTK_REAL :: charmin, charmax, charpm,chartop,avg_alp,avg_det, sdet
- CCTK_REAL :: gxxh, gxyh, gxzh, gyyh, gyzh, gzzh, uxxh, uxyh, &
- uxzh, uyyh, uyzh, uzzh, avg_beta, usendh
- CCTK_REAL :: rhoenth_p, rhoenth_m, avg_betax, avg_betay, avg_betaz
- CCTK_REAL :: vxtp,vytp,vztp,vxtm,vytm,vztm,ab0p,ab0m,b2p,b2m,bdotvp,bdotvm
- CCTK_REAL :: wp,wm,v2p,v2m,bxlowp,bxlowm,bylowp,bylowm,bzlowp,bzlowm,vA2m,vA2p
- CCTK_REAL :: Bvecxlowp,Bvecxlowm,Bvecylowp,Bvecylowm,Bveczlowp,Bveczlowm
- CCTK_REAL :: pressstarp,pressstarm,velxlowp,velxlowm,velylowp,velylowm,velzlowp,velzlowm
+ CCTK_REAL, dimension(3,3) :: gh, uh
+ CCTK_REAL :: usendh
+ CCTK_REAL :: rhoenth_p, rhoenth_m
+ CCTK_REAL, dimension(3) :: avg_beta
+ CCTK_REAL, dimension(3) :: vtp,vtm,blowp,blowm,Bveclowp,Bveclowm
+ CCTK_REAL, dimension(3) :: vellowp,vellowm
+ CCTK_REAL :: ab0p,ab0m,b2p,b2m,bdotvp,bdotvm
+ CCTK_REAL :: wp,wm,v2p,v2m,vA2m,vA2p
+ CCTK_REAL :: pressstarp,pressstarm
CCTK_REAL :: entropyconsp,entropyconsm,entropyp,entropym,entropyf,entropydiff,entropyfp,entropyfm
@@ -61,6 +63,8 @@
CCTK_INT :: type_bits, trivial
CCTK_REAL :: xtemp
+ integer :: ix,iy,iz
+
! save memory when MP is not used
CCTK_INT :: GRHydro_UseGeneralCoordinates
CCTK_REAL, DIMENSION(cctk_ash1,cctk_ash2,cctk_ash3) :: g11, g12, g13, g22, g23, g33
@@ -105,14 +109,17 @@
call SpaceMask_GetTypeBits(type_bits, "Hydro_RiemannProblemX")
call SpaceMask_GetStateBits(trivial, "Hydro_RiemannProblemX", &
&"trivial")
+ ix = 1 ; iy = 2 ; iz = 3
else if (flux_direction == 2) then
call SpaceMask_GetTypeBits(type_bits, "Hydro_RiemannProblemY")
call SpaceMask_GetStateBits(trivial, "Hydro_RiemannProblemY", &
&"trivial")
+ ix = 2 ; iy = 3 ; iz = 1
else if (flux_direction == 3) then
call SpaceMask_GetTypeBits(type_bits, "Hydro_RiemannProblemZ")
call SpaceMask_GetStateBits(trivial, "Hydro_RiemannProblemZ", &
&"trivial")
+ ix = 3 ; iy = 1 ; iz = 2
else
call CCTK_ERROR("Flux direction not x,y,z")
end if
@@ -120,19 +127,29 @@
! constraint transport needs to be able to average fluxes in the directions
! other that flux_direction
- !$OMP PARALLEL DO PRIVATE(k,j,i,f1,lamminus,lamplus,cons_p,cons_m,fplus,fminus,qdiff,psidcf,psidcp,psidcm,prim_p,prim_m,&
- !$OMP avg_betax,avg_betay,avg_betaz,avg_beta,avg_alp,&
- !$OMP gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,avg_det,sdet,uxxh,uxyh,uxzh,uyyh,uyzh,uzzh,&
- !$OMP vxtp,vxtm,vytp,vytm,vztp,vztm,&
- !$OMP velxlowp,velxlowm,Bvecxlowp,Bvecxlowm,&
- !$OMP velylowp,velylowm,Bvecylowp,Bvecylowm,&
- !$OMP velzlowp,velzlowm,Bveczlowp,Bveczlowm,&
+ !$OMP PARALLEL PRIVATE(k,j,i,f1,lamminus,lamplus,cons_p,cons_m,fplus,fminus,qdiff,psidcf,psidcp,psidcm,prim_p,prim_m,&
+ !$OMP avg_beta,avg_alp,&
+ !$OMP gh,avg_det,sdet,uh,&
+ !$OMP vtp,vtm,vellowp,vellowm,Bveclowp,Bveclowm,&
!$OMP bdotvp,bdotvm,b2p,b2m,v2p,v2m,wp,wm,&
- !$OMP bxlowp,bxlowm,bylowp,bylowm,bzlowp,bzlowm,&
+ !$OMP blowp,blowm,&
!$OMP rhoenth_p,rhoenth_m,ab0p,ab0m,vA2p,vA2m,pressstarp,pressstarm,&
!$OMP usendh,psidcdiff,psidcfp,psidcfm,charmin,charmax,chartop,charpm,&
!$OMP charmin_dc,charmax_dc,charpm_dc,m,xtemp,&
!$OMP entropyconsp,entropyconsm,entropyp,entropym,entropyf,entropydiff,entropyfp,entropyfm)
+ ! avoid compiler warnings
+ psidcf = 0d0
+ psidcp = 0.d0
+ psidcm = 0d0
+ psidcfp = 0d0
+ psidcfm = 0d0
+
+ entropyf = 0d0
+ entropyfp = 0d0
+ entropyfm = 0d0
+ entropyconsm = 0d0
+ entropyconsp = 0d0
+ !$OMP DO
do k = GRHydro_stencil, cctk_lsh(3) - GRHydro_stencil + transport_constraints*(1-zoffset)
do j = GRHydro_stencil, cctk_lsh(2) - GRHydro_stencil + transport_constraints*(1-yoffset)
do i = GRHydro_stencil, cctk_lsh(1) - GRHydro_stencil + transport_constraints*(1-xoffset)
@@ -146,11 +163,6 @@
fminus = 0.d0
qdiff = 0.d0
- if(clean_divergence.ne.0) then
- psidcp = 0.d0
- psidcm = 0.d0
- endif
-
!!$ Set the left (p for plus) and right (m_i for minus, i+1) states
cons_p(1) = densplus(i,j,k)
@@ -211,54 +223,47 @@
!!$
!!$ In MHD, we need all three shift components regardless of the flux direction
- avg_betax = 0.5d0 * (beta1(i+xoffset,j+yoffset,k+zoffset) + &
+ avg_beta(1) = 0.5d0 * (beta1(i+xoffset,j+yoffset,k+zoffset) + &
beta1(i,j,k))
- avg_betay = 0.5d0 * (beta2(i+xoffset,j+yoffset,k+zoffset) + &
+ avg_beta(2) = 0.5d0 * (beta2(i+xoffset,j+yoffset,k+zoffset) + &
beta2(i,j,k))
- avg_betaz = 0.5d0 * (beta3(i+xoffset,j+yoffset,k+zoffset) + &
+ avg_beta(3) = 0.5d0 * (beta3(i+xoffset,j+yoffset,k+zoffset) + &
beta3(i,j,k))
- if (flux_direction == 1) then
- avg_beta=avg_betax
- else if (flux_direction == 2) then
- avg_beta=avg_betay
- else if (flux_direction == 3) then
- avg_beta=avg_betaz
- else
- call CCTK_ERROR("Flux direction not x,y,z")
- end if
avg_alp = 0.5 * (alp(i,j,k) + alp(i+xoffset,j+yoffset,k+zoffset))
- gxxh = 0.5d0 * (g11(i+xoffset,j+yoffset,k+zoffset) + g11(i,j,k))
- gxyh = 0.5d0 * (g12(i+xoffset,j+yoffset,k+zoffset) + g12(i,j,k))
- gxzh = 0.5d0 * (g13(i+xoffset,j+yoffset,k+zoffset) + g13(i,j,k))
- gyyh = 0.5d0 * (g22(i+xoffset,j+yoffset,k+zoffset) + g22(i,j,k))
- gyzh = 0.5d0 * (g23(i+xoffset,j+yoffset,k+zoffset) + g23(i,j,k))
- gzzh = 0.5d0 * (g33(i+xoffset,j+yoffset,k+zoffset) + g33(i,j,k))
+ gh(1,1) = 0.5d0 * (g11(i+xoffset,j+yoffset,k+zoffset) + g11(i,j,k))
+ gh(1,2) = 0.5d0 * (g12(i+xoffset,j+yoffset,k+zoffset) + g12(i,j,k))
+ gh(1,3) = 0.5d0 * (g13(i+xoffset,j+yoffset,k+zoffset) + g13(i,j,k))
+ gh(2,2) = 0.5d0 * (g22(i+xoffset,j+yoffset,k+zoffset) + g22(i,j,k))
+ gh(2,3) = 0.5d0 * (g23(i+xoffset,j+yoffset,k+zoffset) + g23(i,j,k))
+ gh(3,3) = 0.5d0 * (g33(i+xoffset,j+yoffset,k+zoffset) + g33(i,j,k))
+ gh(2,1) = gh(1,2) ; gh(3,1) = gh(1,3) ; gh(3,2) = gh(2,3)
- avg_det = SPATIAL_DETERMINANT(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh)
+ avg_det = SPATIAL_DETERMINANT(gh(1,1),gh(1,2),gh(1,3),gh(2,2),gh(2,3),gh(3,3))
sdet = sqrt(avg_det)
- call UpperMetric(uxxh, uxyh, uxzh, uyyh, uyzh, uzzh, &
- avg_det,gxxh, gxyh, gxzh, &
- gyyh, gyzh, gzzh)
+ call UpperMetric(uh(1,1), uh(1,2), uh(1,3), uh(2,2), uh(2,3), uh(3,3), &
+ avg_det,gh(1,1), gh(1,2), gh(1,3), &
+ gh(2,2), gh(2,3), gh(3,3))
+ uh(2,1) = uh(1,2) ; uh(3,1) = uh(1,3) ; uh(3,2) = uh(2,3)
- vxtp = prim_p(2)-avg_betax/avg_alp
- vytp = prim_p(3)-avg_betay/avg_alp
- vztp = prim_p(4)-avg_betaz/avg_alp
- vxtm = prim_m(2)-avg_betax/avg_alp
- vytm = prim_m(3)-avg_betay/avg_alp
- vztm = prim_m(4)-avg_betaz/avg_alp
+ vtp(1) = prim_p(2)-avg_beta(1)/avg_alp
+ vtp(2) = prim_p(3)-avg_beta(2)/avg_alp
+ vtp(3) = prim_p(4)-avg_beta(3)/avg_alp
+ vtm(1) = prim_m(2)-avg_beta(1)/avg_alp
+ vtm(2) = prim_m(3)-avg_beta(2)/avg_alp
+ vtm(3) = prim_m(4)-avg_beta(3)/avg_alp
- call calc_vlow_blow(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh, &
+ call calc_vlow_blow(gh(1,1),gh(1,2),gh(1,3),gh(2,2),gh(2,3),gh(3,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, &
+ vellowp(1),vellowp(2),vellowp(3),Bveclowp(1),Bveclowp(2),Bveclowp(3), &
+ bdotvp,b2p,v2p,wp,blowp(1),blowp(2),blowp(3))
+ call calc_vlow_blow(gh(1,1),gh(1,2),gh(1,3),gh(2,2),gh(2,3),gh(3,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)
+ vellowm(1),vellowm(2),vellowm(3),Bveclowm(1),Bveclowm(2),Bveclowm(3), &
+ bdotvm,b2m,v2m,wm,blowm(1),blowm(2),blowm(3))
rhoenth_p = prim_p(1)*(1.d0+prim_p(5))+prim_p(6)
rhoenth_m = prim_m(1)*(1.d0+prim_m(5))+prim_m(6)
@@ -283,80 +288,30 @@
!!$ pressstar instead of P, b_i, alp b^0, w, metric determinant,
!!$ alp, and beta in the flux dir
- if (flux_direction == 1) then
- call num_x_fluxM(cons_p(1),cons_p(2),cons_p(3),cons_p(4),cons_p(5),&
- 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) + 1.0d0*sdet*uxxh*psidcp - cons_p(6)*avg_betax/avg_alp
- f1(7)=f1(7) + 1.0d0*sdet*uxyh*psidcp - cons_p(6)*avg_betay/avg_alp
- f1(8)=f1(8) + 1.0d0*sdet*uxzh*psidcp - cons_p(6)*avg_betaz/avg_alp
- psidcf = cons_p(6)/sdet-psidcp*avg_betax/avg_alp
- endif
- if(evolve_entropy.ne.0) then
- entropyf = entropyconsp*vxtp
- 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),&
- 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) + 1.0d0*sdet*uxyh*psidcp - cons_p(7)*avg_betax/avg_alp
- f1(7)=f1(7) + 1.0d0*sdet*uyyh*psidcp - cons_p(7)*avg_betay/avg_alp
- f1(8)=f1(8) + 1.0d0*sdet*uyzh*psidcp - cons_p(7)*avg_betaz/avg_alp
- psidcf = cons_p(7)/sdet-psidcp*avg_betay/avg_alp
- endif
- if(evolve_entropy.ne.0) then
- entropyf = entropyconsp*vytp
- 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),&
- 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) + 1.0d0*sdet*uxzh*psidcp - cons_p(8)*avg_betax/avg_alp
- f1(7)=f1(7) + 1.0d0*sdet*uyzh*psidcp - cons_p(8)*avg_betay/avg_alp
- f1(8)=f1(8) + 1.0d0*sdet*uzzh*psidcp - cons_p(8)*avg_betaz/avg_alp
- psidcf = cons_p(8)/sdet-psidcp*avg_betaz/avg_alp
- endif
- if(evolve_entropy.ne.0) then
- entropyf = entropyconsp*vztp
- endif
-
- else
- call CCTK_ERROR("Flux direction not x,y,z")
- end if
+ call num_x_fluxM(cons_p(1),cons_p(1+ix),cons_p(1+iy),cons_p(1+iz),&
+ cons_p(5),cons_p(5+ix),cons_p(5+iy),cons_p(5+iz),&
+ f1(1),f1(1+ix),f1(1+iy),f1(1+iz),f1(5),f1(5+ix),f1(5+iy),&
+ f1(5+iz),&
+ vtp(ix),vtp(iy),vtp(iz),pressstarp,blowp(ix),blowp(iy),&
+ blowp(iz),ab0p,wp, &
+ avg_det,avg_alp,avg_beta(flux_direction))
+ if(clean_divergence.ne.0) then
+ f1(6)=f1(6) + 1.0d0*sdet*uh(flux_direction,1)*psidcp - cons_p(5+flux_direction)*avg_beta(1)/avg_alp
+ f1(7)=f1(7) + 1.0d0*sdet*uh(flux_direction,2)*psidcp - cons_p(5+flux_direction)*avg_beta(2)/avg_alp
+ f1(8)=f1(8) + 1.0d0*sdet*uh(flux_direction,3)*psidcp - cons_p(5+flux_direction)*avg_beta(3)/avg_alp
+ psidcf = cons_p(5+flux_direction)/sdet-psidcp*avg_beta(flux_direction)/avg_alp
+ endif
+ if(evolve_entropy.ne.0) then
+ entropyf = entropyconsp*vtp(flux_direction)
+ endif
else !!! The end of this branch is right at the bottom of the routine
- if (flux_direction == 1) then
- usendh = uxxh
- else if (flux_direction == 2) then
- usendh = uyyh
- else if (flux_direction == 3) then
- usendh = uzzh
- else
- call CCTK_ERROR("Flux direction not x,y,z")
- end if
+ usendh = uh(flux_direction,flux_direction)
!!$ Calculate the jumps in the conserved variables
- qdiff(1) = cons_m(1) - cons_p(1)
- qdiff(2) = cons_m(2) - cons_p(2)
- 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 = cons_m - cons_p
if (clean_divergence.ne.0) then
psidcdiff = psidcm - psidcp
@@ -365,181 +320,68 @@
entropydiff = entropyconsm - entropyconsp
endif
-!!$ Eigenvalues and fluxes either side of the cell interface
+!!$ Eigenvalues and fluxes either side of the cell interface
- if (flux_direction == 1) then
- if(evolve_temper.ne.1) then
- call eigenvaluesM(GRHydro_eos_handle,&
- prim_m(1),prim_m(2),prim_m(3),prim_m(4),prim_m(5),prim_m(6),prim_m(7), &
- 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(6),prim_p(7), &
- prim_p(8),prim_p(9),prim_p(10),&
- lamplus,gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,&
- usendh,avg_alp,avg_beta)
- else
- xtemp = temperature(i,j,k)
- call eigenvaluesM_hot(GRHydro_eos_handle,i,j,k,&
- prim_m(1),prim_m(2),prim_m(3),prim_m(4),prim_m(5),prim_m(6),prim_m(7), &
- prim_m(8),prim_m(9),prim_m(10),&
- xtemp,y_e_minus(i+xoffset,j+yoffset,k+zoffset),&
- lamminus,gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,&
- usendh,avg_alp,avg_beta)
- xtemp = temperature(i,j,k)
- call eigenvaluesM_hot(GRHydro_eos_handle,i,j,k,&
- prim_p(1),prim_p(2),prim_p(3),prim_p(4),prim_p(5),prim_p(6),prim_p(7), &
- prim_p(8),prim_p(9),prim_p(10),&
- xtemp,y_e_plus(i,j,k),&
- lamplus,gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,&
- usendh,avg_alp,avg_beta)
- endif
+ if(evolve_temper.ne.1) then
+ call eigenvaluesM(GRHydro_eos_handle,&
+ prim_m(1),prim_m(1+ix),prim_m(1+iy),prim_m(1+iz),prim_m(5),prim_m(6),prim_m(7), &
+ prim_m(7+ix),prim_m(7+iy),prim_m(7+iz),&
+ lamminus,gh(1,1),gh(1,2),gh(1,3),gh(2,2),gh(2,3),gh(3,3),&
+ usendh,avg_alp,avg_beta(flux_direction))
+ call eigenvaluesM(GRHydro_eos_handle,&
+ prim_p(1),prim_p(1+ix),prim_p(1+iy),prim_p(1+iz),prim_p(5),prim_p(6),prim_p(7), &
+ prim_p(7+ix),prim_p(7+iy),prim_p(7+iz),&
+ lamplus,gh(1,1),gh(1,2),gh(1,3),gh(2,2),gh(2,3),gh(3,3),&
+ usendh,avg_alp,avg_beta(flux_direction))
+ else
+ xtemp = temperature(i,j,k)
+ call eigenvaluesM_hot(GRHydro_eos_handle,i,j,k,&
+ prim_m(1),prim_m(1+ix),prim_m(1+iy),prim_m(1+iz),prim_m(5),prim_m(6),prim_m(7), &
+ prim_m(7+ix),prim_m(7+iy),prim_m(7+iz),&
+ xtemp,y_e_minus(i+xoffset,j+yoffset,k+zoffset),&
+ lamminus,gh(1,1),gh(1,2),gh(1,3),gh(2,2),gh(2,3),gh(3,3),&
+ usendh,avg_alp,avg_beta(flux_direction))
+ xtemp = temperature(i,j,k)
+ call eigenvaluesM_hot(GRHydro_eos_handle,i,j,k,&
+ prim_p(1),prim_p(1+ix),prim_p(1+iy),prim_p(1+iz),prim_p(5),prim_p(6),prim_p(7), &
+ prim_p(7+ix),prim_p(7+iy),prim_p(7+iz),&
+ xtemp,y_e_plus(i,j,k),&
+ lamplus,gh(1,1),gh(1,2),gh(1,3),gh(2,2),gh(2,3),gh(3,3),&
+ usendh,avg_alp,avg_beta(flux_direction))
+ endif
- call num_x_fluxM(cons_p(1),cons_p(2),cons_p(3),cons_p(4),cons_p(5),&
- 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),&
- 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)
+ call num_x_fluxM(cons_p(1),cons_p(1+ix),cons_p(1+iy),cons_p(1+iz),&
+ cons_p(5),&
+ cons_p(5+ix),cons_p(5+iy),cons_p(5+iz),&
+ fplus(1),fplus(1+ix),fplus(1+iy),fplus(1+iz),fplus(5),&
+ fplus(5+ix),fplus(5+iy),fplus(5+iz),&
+ vtp(ix),vtp(iy),vtp(iz),pressstarp,blowp(ix),blowp(iy),&
+ blowp(iz),ab0p,wp, &
+ avg_det,avg_alp,avg_beta(flux_direction))
+ call num_x_fluxM(cons_m(1),cons_m(1+ix),cons_m(1+iy),cons_m(1+iz),&
+ cons_m(5),&
+ cons_m(5+ix),cons_m(5+iy),cons_m(5+iz),&
+ fminus(1),fminus(1+ix),fminus(1+iy),fminus(1+iz),fminus(5),&
+ fminus(5+ix),fminus(5+iy),fminus(5+iz),&
+ vtm(ix),vtm(iy),vtm(iz),pressstarm,blowm(ix),blowm(iy),&
+ blowm(iz),ab0m,wm, &
+ avg_det,avg_alp,avg_beta(flux_direction))
- if(clean_divergence.ne.0) then
- fminus(6)=fminus(6) + 1.0d0*sdet*uxxh*psidcm - cons_m(6)*avg_betax/avg_alp
- fminus(7)=fminus(7) + 1.0d0*sdet*uxyh*psidcm - cons_m(6)*avg_betay/avg_alp
- fminus(8)=fminus(8) + 1.0d0*sdet*uxzh*psidcm - cons_m(6)*avg_betaz/avg_alp
- fplus(6)=fplus(6) + 1.0d0*sdet*uxxh*psidcp - cons_p(6)*avg_betax/avg_alp
- fplus(7)=fplus(7) + 1.0d0*sdet*uxyh*psidcp - cons_p(6)*avg_betay/avg_alp
- fplus(8)=fplus(8) + 1.0d0*sdet*uxzh*psidcp - cons_p(6)*avg_betaz/avg_alp
- psidcfp = cons_p(6)/sdet-avg_betax*psidcp/avg_alp
- psidcfm = cons_m(6)/sdet-avg_betax*psidcm/avg_alp
- endif
- if(evolve_entropy.ne.0) then
- entropyfp = entropyconsp*vxtp
- entropyfm = entropyconsm*vxtm
- endif
+ if(clean_divergence.ne.0) then
+ fminus(6)=fminus(6) + 1.0d0*sdet*uh(flux_direction,1)*psidcm - cons_m(5+flux_direction)*avg_beta(1)/avg_alp
+ fminus(7)=fminus(7) + 1.0d0*sdet*uh(flux_direction,2)*psidcm - cons_m(5+flux_direction)*avg_beta(2)/avg_alp
+ fminus(8)=fminus(8) + 1.0d0*sdet*uh(flux_direction,3)*psidcm - cons_m(5+flux_direction)*avg_beta(3)/avg_alp
+ fplus(6)=fplus(6) + 1.0d0*sdet*uh(flux_direction,1)*psidcp - cons_p(5+flux_direction)*avg_beta(1)/avg_alp
+ fplus(7)=fplus(7) + 1.0d0*sdet*uh(flux_direction,2)*psidcp - cons_p(5+flux_direction)*avg_beta(2)/avg_alp
+ fplus(8)=fplus(8) + 1.0d0*sdet*uh(flux_direction,3)*psidcp - cons_p(5+flux_direction)*avg_beta(3)/avg_alp
+ psidcfp = cons_p(5+flux_direction)/sdet-avg_beta(flux_direction)*psidcp/avg_alp
+ psidcfm = cons_m(5+flux_direction)/sdet-avg_beta(flux_direction)*psidcm/avg_alp
+ endif
+ if(evolve_entropy.ne.0) then
+ entropyfp = entropyconsp*vtp(flux_direction)
+ entropyfm = entropyconsm*vtm(flux_direction)
+ endif
- else if (flux_direction == 2) then
- if(evolve_temper.ne.1) then
- call eigenvaluesM(GRHydro_eos_handle,&
- prim_m(1),prim_m(3),prim_m(4),prim_m(2),prim_m(5),prim_m(6),prim_m(7), &
- 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(6),prim_p(7), &
- prim_p(9),prim_p(10),prim_p(8),&
- lamplus,gyyh,gyzh,gxyh,gzzh,gxzh,gxxh,&
- usendh,avg_alp,avg_beta)
- else
- xtemp = temperature(i,j,k)
- call eigenvaluesM_hot(GRHydro_eos_handle,i,j,k,&
- prim_m(1),prim_m(3),prim_m(4),prim_m(2),prim_m(5),prim_m(6),prim_m(7), &
- prim_m(9),prim_m(10),prim_m(8),&
- xtemp,y_e_minus(i+xoffset,j+yoffset,k+zoffset),&
- lamminus,gyyh,gyzh,gxyh,gzzh,gxzh,gxxh,&
- usendh,avg_alp,avg_beta)
- xtemp = temperature(i,j,k)
- call eigenvaluesM_hot(GRHydro_eos_handle,i,j,k,&
- prim_p(1),prim_p(3),prim_p(4),prim_p(2),prim_p(5),prim_p(6),prim_p(7), &
- prim_p(9),prim_p(10),prim_p(8),&
- xtemp,y_e_plus(i,j,k),&
- lamplus,gyyh,gyzh,gxyh,gzzh,gxzh,gxxh,&
- usendh,avg_alp,avg_beta)
- endif
-
- call num_x_fluxM(cons_p(1),cons_p(3),cons_p(4),cons_p(2),cons_p(5),&
- 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),&
- 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) + 1.0d0*sdet*uxyh*psidcm - cons_m(7)*avg_betax/avg_alp
- fminus(7)=fminus(7) + 1.0d0*sdet*uyyh*psidcm - cons_m(7)*avg_betay/avg_alp
- fminus(8)=fminus(8) + 1.0d0*sdet*uyzh*psidcm - cons_m(7)*avg_betaz/avg_alp
- fplus(6)=fplus(6) + 1.0d0*sdet*uxyh*psidcp - cons_p(7)*avg_betax/avg_alp
- fplus(7)=fplus(7) + 1.0d0*sdet*uyyh*psidcp - cons_p(7)*avg_betay/avg_alp
- fplus(8)=fplus(8) + 1.0d0*sdet*uyzh*psidcp - cons_p(7)*avg_betaz/avg_alp
- psidcfp = cons_p(7)/sdet-avg_betay*psidcp/avg_alp
- psidcfm = cons_m(7)/sdet-avg_betay*psidcm/avg_alp
- endif
- if(evolve_entropy.ne.0) then
- entropyfp = entropyconsp*vytp
- entropyfm = entropyconsm*vytm
- endif
-
- else if (flux_direction == 3) then
- if(evolve_temper.ne.1) then
- call eigenvaluesM(GRHydro_eos_handle,&
- prim_m(1),prim_m(4),prim_m(2),prim_m(3),prim_m(5),prim_m(6),prim_m(7), &
- 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(6),prim_p(7), &
- prim_p(10),prim_p(8),prim_p(9),&
- lamplus,gzzh,gxzh,gyzh,gxxh,gxyh,gyyh,&
- usendh,avg_alp,avg_beta)
- else
- xtemp = temperature(i,j,k)
- call eigenvaluesM_hot(GRHydro_eos_handle,i,j,k,&
- prim_m(1),prim_m(4),prim_m(2),prim_m(3),prim_m(5),prim_m(6),prim_m(7), &
- prim_m(10),prim_m(8),prim_m(9),&
- xtemp,y_e_minus(i+xoffset,j+yoffset,k+zoffset),&
- lamminus,gzzh,gxzh,gyzh,gxxh,gxyh,gyyh,&
- usendh,avg_alp,avg_beta)
- xtemp = temperature(i,j,k)
- call eigenvaluesM_hot(GRHydro_eos_handle,i,j,k,&
- prim_p(1),prim_p(4),prim_p(2),prim_p(3),prim_p(5),prim_p(6),prim_p(7), &
- prim_p(10),prim_p(8),prim_p(9),&
- xtemp,y_e_plus(i,j,k),&
- lamplus,gzzh,gxzh,gyzh,gxxh,gxyh,gyyh,&
- usendh,avg_alp,avg_beta)
- endif
-
- call num_x_fluxM(cons_p(1),cons_p(4),cons_p(2),cons_p(3),cons_p(5),&
- 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),&
- 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) + 1.0d0*sdet*uxzh*psidcm - cons_m(8)*avg_betax/avg_alp
- fminus(7)=fminus(7) + 1.0d0*sdet*uyzh*psidcm - cons_m(8)*avg_betay/avg_alp
- fminus(8)=fminus(8) + 1.0d0*sdet*uzzh*psidcm - cons_m(8)*avg_betaz/avg_alp
- fplus(6)=fplus(6) + 1.0d0*sdet*uxzh*psidcp - cons_p(8)*avg_betax/avg_alp
- fplus(7)=fplus(7) + 1.0d0*sdet*uyzh*psidcp - cons_p(8)*avg_betay/avg_alp
- fplus(8)=fplus(8) + 1.0d0*sdet*uzzh*psidcp - cons_p(8)*avg_betaz/avg_alp
- psidcfp = cons_p(8)/sdet-avg_betaz*psidcp/avg_alp
- psidcfm = cons_m(8)/sdet-avg_betaz*psidcm/avg_alp
- endif
- if(evolve_entropy.ne.0) then
- entropyfp = entropyconsp*vztp
- entropyfm = entropyconsm*vztm
- endif
-
- else
- call CCTK_ERROR("Flux direction not x,y,z")
- end if
-
-
!!$ Find minimum and maximum wavespeeds
charmin = min(0.d0, lamplus(1), lamplus(2), lamplus(3), &
@@ -589,8 +431,8 @@
!!$ that the normal vector to the characteristic hypersurface
!!$ be spacelike (Eq. 60 of Anton et al.)
- charmax_dc = sqrt(usendh) - avg_beta/avg_alp
- charmin_dc = -1.d0*sqrt(usendh) - avg_beta/avg_alp
+ charmax_dc = sqrt(usendh) - avg_beta(flux_direction)/avg_alp
+ charmin_dc = -1.d0*sqrt(usendh) - avg_beta(flux_direction)/avg_alp
charpm_dc = charmax_dc - charmin_dc
psidcf = (charmax_dc * psidcfp - charmin_dc * psidcfm + &
@@ -664,7 +506,8 @@
end do
end do
end do
- !$OMP END PARALLEL DO
+ !$OMP END DO
+ !$OMP END PARALLEL
#undef faulty_gxx
#undef faulty_gxy
#undef faulty_gxz
