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