diff options
Diffstat (limited to 'src/GRHydro_HLLEM.F90')
| -rw-r--r-- | src/GRHydro_HLLEM.F90 | 173 |
1 files changed, 114 insertions, 59 deletions
diff --git a/src/GRHydro_HLLEM.F90 b/src/GRHydro_HLLEM.F90 index 3c694ca..87df4d9 100644 --- a/src/GRHydro_HLLEM.F90 +++ b/src/GRHydro_HLLEM.F90 @@ -55,9 +55,10 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS) CCTK_REAL :: pressstarp,pressstarm,velxlowp,velxlowm,velylowp,velylowm,velzlowp,velzlowm CCTK_REAL :: psidcp, psidcm, psidcf, psidcdiff, psidcfp, psidcfm - CCTK_REAL :: charmax_dc, charmin_dc, charpm_dc - + CCTK_REAL :: charmax_dc, charmin_dc, charpm_dc + CCTK_INT :: type_bits, trivial, not_trivial + CCTK_REAL :: xtemp if (flux_direction == 1) then call SpaceMask_GetTypeBits(type_bits, "Hydro_RiemannProblemX") @@ -227,9 +228,9 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS) 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.0*sdet*uxxh*psidcp - cons_p(6)*avg_betax/avg_alp - f1(7)=f1(7) + 1.0*sdet*uxyh*psidcp - cons_p(6)*avg_betay/avg_alp - f1(8)=f1(8) + 1.0*sdet*uxzh*psidcp - cons_p(6)*avg_betaz/avg_alp + 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 @@ -240,9 +241,9 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS) 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.0*sdet*uxyh*psidcp - cons_p(7)*avg_betax/avg_alp - f1(7)=f1(7) + 1.0*sdet*uyyh*psidcp - cons_p(7)*avg_betay/avg_alp - f1(8)=f1(8) + 1.0*sdet*uyzh*psidcp - cons_p(7)*avg_betaz/avg_alp + 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 @@ -253,9 +254,9 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS) 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.0*sdet*uxzh*psidcp - cons_p(8)*avg_betax/avg_alp - f1(7)=f1(7) + 1.0*sdet*uyzh*psidcp - cons_p(8)*avg_betay/avg_alp - f1(8)=f1(8) + 1.0*sdet*uzzh*psidcp - cons_p(8)*avg_betaz/avg_alp + 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 @@ -293,16 +294,34 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS) !!$ Eigenvalues and fluxes either side of the cell interface 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), & - 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), & - prim_p(8),prim_p(9),prim_p(10),& - lamplus,gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,& - usendh,avg_alp,avg_beta) + 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(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(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+xoffset,j+yoffset,k+zoffset) + 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(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(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 + 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),& @@ -316,27 +335,45 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS) avg_det,avg_alp,avg_beta) if(clean_divergence.ne.0) then - fminus(6)=fminus(6) + 1.0*sdet*uxxh*psidcm - cons_m(6)*avg_betax/avg_alp - fminus(7)=fminus(7) + 1.0*sdet*uxyh*psidcm - cons_m(6)*avg_betay/avg_alp - fminus(8)=fminus(8) + 1.0*sdet*uxzh*psidcm - cons_m(6)*avg_betaz/avg_alp - fplus(6)=fplus(6) + 1.0*sdet*uxxh*psidcp - cons_p(6)*avg_betax/avg_alp - fplus(7)=fplus(7) + 1.0*sdet*uxyh*psidcp - cons_p(6)*avg_betay/avg_alp - fplus(8)=fplus(8) + 1.0*sdet*uxzh*psidcp - cons_p(6)*avg_betaz/avg_alp + 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 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), & - 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), & - prim_p(9),prim_p(10),prim_p(8),& - lamplus,gyyh,gyzh,gxyh,gzzh,gxzh,gxxh,& - usendh,avg_alp,avg_beta) + 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(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(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+xoffset,j+yoffset,k+zoffset) + 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(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(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),& @@ -350,27 +387,45 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS) avg_det,avg_alp,avg_beta) if(clean_divergence.ne.0) then - fminus(6)=fminus(6) + 1.0*sdet*uxyh*psidcm - cons_m(7)*avg_betax/avg_alp - fminus(7)=fminus(7) + 1.0*sdet*uyyh*psidcm - cons_m(7)*avg_betay/avg_alp - fminus(8)=fminus(8) + 1.0*sdet*uyzh*psidcm - cons_m(7)*avg_betaz/avg_alp - fplus(6)=fplus(6) + 1.0*sdet*uxyh*psidcp - cons_p(7)*avg_betax/avg_alp - fplus(7)=fplus(7) + 1.0*sdet*uyyh*psidcp - cons_p(7)*avg_betay/avg_alp - fplus(8)=fplus(8) + 1.0*sdet*uyzh*psidcp - cons_p(7)*avg_betaz/avg_alp + 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 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), & - 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), & - prim_p(10),prim_p(8),prim_p(9),& - lamplus,gzzh,gxzh,gyzh,gxxh,gxyh,gyyh,& - usendh,avg_alp,avg_beta) + 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(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(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+xoffset,j+yoffset,k+zoffset) + 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(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(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), & @@ -384,12 +439,12 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS) avg_det,avg_alp,avg_beta) if(clean_divergence.ne.0) then - fminus(6)=fminus(6) + 1.0*sdet*uxzh*psidcm - cons_m(8)*avg_betax/avg_alp - fminus(7)=fminus(7) + 1.0*sdet*uyzh*psidcm - cons_m(8)*avg_betay/avg_alp - fminus(8)=fminus(8) + 1.0*sdet*uzzh*psidcm - cons_m(8)*avg_betaz/avg_alp - fplus(6)=fplus(6) + 1.0*sdet*uxzh*psidcp - cons_p(8)*avg_betax/avg_alp - fplus(7)=fplus(7) + 1.0*sdet*uyzh*psidcp - cons_p(8)*avg_betay/avg_alp - fplus(8)=fplus(8) + 1.0*sdet*uzzh*psidcp - cons_p(8)*avg_betaz/avg_alp + 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 |