diff options
Diffstat (limited to 'src/GRHydro_Con2PrimAM.F90')
| -rw-r--r-- | src/GRHydro_Con2PrimAM.F90 | 99 |
1 files changed, 49 insertions, 50 deletions
diff --git a/src/GRHydro_Con2PrimAM.F90 b/src/GRHydro_Con2PrimAM.F90 index 7d2103e..374d794 100644 --- a/src/GRHydro_Con2PrimAM.F90 +++ b/src/GRHydro_Con2PrimAM.F90 @@ -183,10 +183,9 @@ subroutine Conservative2PrimitiveAM(CCTK_ARGUMENTS) epsnegative = 0 - det = SPATIAL_DETERMINANT(g11(i,j,k),g12(i,j,k),g13(i,j,k),g22(i,j,k),g23(i,j,k),g33(i,j,k)) - sdet = sqrt(det) + sdet = sdetg(i,j,k) - call UpperMetric(uxx,uxy,uxz,uyy,uyz,uzz,det,& + call UpperMetric(uxx,uxy,uxz,uyy,uyz,uzz,sdet*sdet,& g11(i,j,k),g12(i,j,k),g13(i,j,k),g22(i,j,k),& g23(i,j,k),g33(i,j,k)) @@ -356,7 +355,7 @@ subroutine Conservative2PrimitiveAM(CCTK_ARGUMENTS) velx_tmp,vely_tmp,velz_tmp,eps_tmp,press_tmp, & Bvecx_tmp,Bvecy_tmp,Bvecz_tmp,b2, w_lorentz_tmp,& g11(i,j,k),g12(i,j,k),g13(i,j,k),g22(i,j,k),g23(i,j,k),g33(i,j,k), & - uxx,uxy,uxz,uyy,uyz,uzz,det, & + uxx,uxy,uxz,uyy,uyz,uzz,sdet, & epsnegative,GRHydro_C2P_failed(i,j,k)) else @@ -367,7 +366,7 @@ subroutine Conservative2PrimitiveAM(CCTK_ARGUMENTS) eps_tmp,press_tmp,Bvecx_tmp,Bvecy_tmp,Bvecz_tmp,b2,& w_lorentz_tmp,g11(i,j,k),g12(i,j,k),g13(i,j,k),& g22(i,j,k),g23(i,j,k),g33(i,j,k), & - uxx,uxy,uxz,uyy,uyz,uzz,det, & + uxx,uxy,uxz,uyy,uyz,uzz,sdet, & epsnegative,GRHydro_C2P_failed(i,j,k)) endif @@ -477,7 +476,7 @@ subroutine Conservative2PrimitiveAM(CCTK_ARGUMENTS) eps_tmp,press_tmp,Bvecx_tmp,Bvecy_tmp,Bvecz_tmp,b2,& w_lorentz_tmp,g11(i,j,k),g12(i,j,k),g13(i,j,k),& g22(i,j,k),g23(i,j,k),g33(i,j,k), & - uxx,uxy,uxz,uyy,uyz,uzz,det, & + uxx,uxy,uxz,uyy,uyz,uzz,sdet, & epsnegative,GRHydro_C2P_failed(i,j,k)) else call GRHydro_Con2PrimM_Polytype_pt(GRHydro_eos_handle, local_pgam, & @@ -486,7 +485,7 @@ subroutine Conservative2PrimitiveAM(CCTK_ARGUMENTS) velx_tmp,vely_tmp,velz_tmp,eps_tmp,press_tmp,& Bvecx_tmp,Bvecy_tmp,Bvecz_tmp,b2,w_lorentz_tmp,& g11(i,j,k),g12(i,j,k),g13(i,j,k),g22(i,j,k),g23(i,j,k),g33(i,j,k), & - uxx,uxy,uxz,uyy,uyz,uzz,det, & + uxx,uxy,uxz,uyy,uyz,uzz,sdet, & epsnegative,GRHydro_C2P_failed(i,j,k)) end if @@ -529,7 +528,7 @@ subroutine Conservative2PrimitiveAM(CCTK_ARGUMENTS) eps_tmp,press_tmp,Bvecx_tmp,Bvecy_tmp,Bvecz_tmp,b2,& w_lorentz_tmp,g11(i,j,k),g12(i,j,k),g13(i,j,k),& g22(i,j,k),g23(i,j,k),g33(i,j,k), & - uxx,uxy,uxz,uyy,uyz,uzz,det, & + uxx,uxy,uxz,uyy,uyz,uzz,sdet, & epsnegative,GRHydro_C2P_failed(i,j,k)) if(GRHydro_C2P_failed(i,j,k).ne.0) then ! this means c2p did not converge. @@ -558,7 +557,7 @@ subroutine Conservative2PrimitiveAM(CCTK_ARGUMENTS) eps_tmp,press_tmp,Bvecx_tmp,Bvecy_tmp,Bvecz_tmp,b2,& w_lorentz_tmp,g11(i,j,k),g12(i,j,k),g13(i,j,k),& g22(i,j,k),g23(i,j,k),g33(i,j,k), & - uxx,uxy,uxz,uyy,uyz,uzz,det, & + uxx,uxy,uxz,uyy,uyz,uzz,sdet, & epsnegative,GRHydro_C2P_failed(i,j,k)) if(GRHydro_C2P_failed(i,j,k).ne.0) then !$OMP CRITICAL @@ -603,7 +602,7 @@ subroutine Conservative2PrimitiveAM(CCTK_ARGUMENTS) !$OMP END CRITICAL ! for safety, let's set the point to atmosphere - dens(i,j,k) = sqrt(det)*GRHydro_rho_min !/(1.d0+GRHydro_atmo_tolerance) + dens(i,j,k) = sdet*GRHydro_rho_min !/(1.d0+GRHydro_atmo_tolerance) rho(i,j,k) = GRHydro_rho_min scon(i,j,k,:) = 0.d0 vup(i,j,k,:) = 0.d0 @@ -621,7 +620,7 @@ subroutine Conservative2PrimitiveAM(CCTK_ARGUMENTS) ! w_lorentz=1, so the expression for tau reduces to [see above]: - tau(i,j,k) = sqrt(det) * (rho(i,j,k)*(1.0+eps(i,j,k)+b2/2.0)) - dens(i,j,k) + tau(i,j,k) = sdet * (rho(i,j,k)*(1.0+eps(i,j,k)+b2/2.0)) - dens(i,j,k) #else ! cott 2010/03/27: ! Honestly, this should never happen. We need to flag the point where @@ -662,7 +661,7 @@ subroutine Conservative2PrimitiveAM(CCTK_ARGUMENTS) velx_tmp,vely_tmp,velz_tmp,eps_tmp,press_tmp,& Bvecx_tmp,Bvecy_tmp,Bvecz_tmp,b2,w_lorentz_tmp,& g11(i,j,k),g12(i,j,k),g13(i,j,k),g22(i,j,k),g23(i,j,k),g33(i,j,k), & - uxx,uxy,uxz,uyy,uyz,uzz,det, & + uxx,uxy,uxz,uyy,uyz,uzz,sdet, & epsnegative,GRHydro_C2P_failed(i,j,k)) #endif @@ -778,8 +777,8 @@ subroutine Conservative2PrimitiveBoundsAM(CCTK_ARGUMENTS) integer :: i, j, k, itracer, nx, ny, nz CCTK_REAL :: uxxl, uxyl, uxzl, uyyl, uyzl, uzzl,& uxxr, uxyr, uxzr, uyyr, uyzr, uzzr, pmin(1), epsmin(1) - CCTK_REAL :: gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_detl,& - gxxr,gxyr,gxzr,gyyr,gyzr,gzzr,avg_detr + CCTK_REAL :: gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_sdetl,& + gxxr,gxyr,gxzr,gyyr,gyzr,gzzr,avg_sdetr CCTK_REAL :: b2minus, b2plus, local_gam(1), local_pgam,local_K,scminus,scplus CCTK_INT :: epsnegative CCTK_REAL :: Bconsx_tmp, Bconsy_tmp, Bconsz_tmp @@ -873,11 +872,11 @@ subroutine Conservative2PrimitiveBoundsAM(CCTK_ARGUMENTS) epsnegative = 0 - avg_detl = SPATIAL_DETERMINANT(gxxl,gxyl,gxzl,gyyl, gyzl,gzzl) - avg_detr = SPATIAL_DETERMINANT(gxxr,gxyr,gxzr,gyyr, gyzr,gzzr) - call UpperMetric(uxxl,uxyl,uxzl,uyyl,uyzl,uzzl,avg_detl,& + avg_sdetl = sqrt(SPATIAL_DETERMINANT(gxxl,gxyl,gxzl,gyyl, gyzl,gzzl)) + avg_sdetr = sqrt(SPATIAL_DETERMINANT(gxxr,gxyr,gxzr,gyyr, gyzr,gzzr)) + call UpperMetric(uxxl,uxyl,uxzl,uyyl,uyzl,uzzl,avg_sdetl*avg_sdetl,& gxxl,gxyl,gxzl,gyyl,gyzl,gzzl) - call UpperMetric(uxxr,uxyr,uxzr,uyyr,uyzr,uzzr,avg_detr,& + call UpperMetric(uxxr,uxyr,uxzr,uyyr,uyzr,uzzr,avg_sdetr*avg_sdetr,& gxxr,gxyr,gxzr,gyyr,gyzr,gzzr) !!$ Tracers get no update for MHD! @@ -899,16 +898,16 @@ subroutine Conservative2PrimitiveBoundsAM(CCTK_ARGUMENTS) Y_e(i,j,k) = Y_e_con(i,j,k) / dens(i,j,k) endif - Bconsx_tmp = sqrt(avg_detl)*Bvecxminus(i,j,k) - Bconsy_tmp = sqrt(avg_detl)*Bvecyminus(i,j,k) - Bconsz_tmp = sqrt(avg_detl)*Bveczminus(i,j,k) + Bconsx_tmp = avg_sdetl*Bvecxminus(i,j,k) + Bconsy_tmp = avg_sdetl*Bvecyminus(i,j,k) + Bconsz_tmp = avg_sdetl*Bveczminus(i,j,k) call GRHydro_Con2PrimM_pt2(GRHydro_eos_handle, keytemp, GRHydro_eos_rf_prec, GRHydro_perc_ptol, local_gam(1), densminus(i,j,k), & sxminus(i,j,k),syminus(i,j,k),szminus(i,j,k), tauminus(i,j,k), & Bconsx_tmp, Bconsy_tmp, Bconsz_tmp, xye(1), xtemp(1), rhominus(i,j,k),& velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k),epsminus(i,j,k),pressminus(i,j,k),& Bvecxminus(i,j,k), Bvecyminus(i,j,k), Bveczminus(i,j,k),b2minus, w_lorentzminus(i,j,k),& gxxl,gxyl,gxzl,gyyl,gyzl,gzzl, & - uxxl,uxyl,uxzl,uyyl,uyzl,uzzl,avg_detl, & + uxxl,uxyl,uxzl,uyyl,uyzl,uzzl,avg_sdetl, & epsnegative,GRHydro_C2P_failed(i,j,k)) if (epsnegative .ne. 0) then @@ -932,7 +931,7 @@ subroutine Conservative2PrimitiveBoundsAM(CCTK_ARGUMENTS) velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k),epsminus(i,j,k),pressminus(i,j,k),& Bvecxminus(i,j,k), Bvecyminus(i,j,k), Bveczminus(i,j,k),b2minus, w_lorentzminus(i,j,k),& gxxl,gxyl,gxzl,gyyl,gyzl,gzzl, & - uxxl,uxyl,uxzl,uyyl,uyzl,uzzl,avg_detl, & + uxxl,uxyl,uxzl,uyyl,uyzl,uzzl,avg_sdetl, & epsnegative,GRHydro_C2P_failed(i,j,k)) end if @@ -960,16 +959,16 @@ subroutine Conservative2PrimitiveBoundsAM(CCTK_ARGUMENTS) epsnegative = 0 - Bconsx_tmp = sqrt(avg_detr)*Bvecxplus(i,j,k) - Bconsy_tmp = sqrt(avg_detr)*Bvecyplus(i,j,k) - Bconsz_tmp = sqrt(avg_detr)*Bveczplus(i,j,k) + Bconsx_tmp = avg_sdetr*Bvecxplus(i,j,k) + Bconsy_tmp = avg_sdetr*Bvecyplus(i,j,k) + Bconsz_tmp = avg_sdetr*Bveczplus(i,j,k) call GRHydro_Con2PrimM_pt2(GRHydro_eos_handle, keytemp, GRHydro_eos_rf_prec, GRHydro_perc_ptol, local_gam(1), densplus(i,j,k), & sxplus(i,j,k),syplus(i,j,k),szplus(i,j,k), tauplus(i,j,k), & Bconsx_tmp, Bconsy_tmp, Bconsz_tmp, xye(1), xtemp(1), rhoplus(i,j,k),& velxplus(i,j,k),velyplus(i,j,k),velzplus(i,j,k),epsplus(i,j,k),pressplus(i,j,k),& Bvecxplus(i,j,k), Bvecyplus(i,j,k), Bveczplus(i,j,k),b2plus, w_lorentzplus(i,j,k),& gxxr,gxyr,gxzr,gyyr,gyzr,gzzr, & - uxxr,uxyr,uxzr,uyyr,uyzr,uzzr,avg_detr, & + uxxr,uxyr,uxzr,uyyr,uyzr,uzzr,avg_sdetr, & epsnegative,GRHydro_C2P_failed(i,j,k)) if (epsnegative .ne. 0) then @@ -993,7 +992,7 @@ subroutine Conservative2PrimitiveBoundsAM(CCTK_ARGUMENTS) velxplus(i,j,k),velyplus(i,j,k),velzplus(i,j,k),epsplus(i,j,k),pressplus(i,j,k),& Bvecxplus(i,j,k), Bvecyplus(i,j,k), Bveczplus(i,j,k),b2plus, w_lorentzplus(i,j,k),& gxxr,gxyr,gxzr,gyyr,gyzr,gzzr, & - uxxr,uxyr,uxzr,uyyr,uyzr,uzzr,avg_detr, & + uxxr,uxyr,uxzr,uyyr,uyzr,uzzr,avg_sdetr, & epsnegative,GRHydro_C2P_failed(i,j,k)) end if @@ -1068,7 +1067,7 @@ subroutine Conservative2PrimitivePolytypeAM(CCTK_ARGUMENTS) DECLARE_CCTK_PARAMETERS integer :: i, j, k, itracer, nx, ny, nz - CCTK_REAL :: uxx, uxy, uxz, uyy, uyz, uzz, det,b2 + CCTK_REAL :: uxx, uxy, uxz, uyy, uyz, uzz, sdet,b2 CCTK_INT :: epsnegative CCTK_REAL :: Bconsx_tmp, Bconsy_tmp, Bconsz_tmp @@ -1140,7 +1139,7 @@ subroutine Conservative2PrimitivePolytypeAM(CCTK_ARGUMENTS) !$OMP PARALLEL DO PRIVATE(i,j,k,itracer,& - !$OMP uxx, uxy, uxz, uyy, uyz, uzz, det, epsnegative, & + !$OMP uxx, uxy, uxz, uyy, uyz, uzz, sdet, epsnegative, & !$OMP Bconsx_tmp, Bconsy_tmp, Bconsz_tmp, & !$OMP b2, xrho, xpress, local_K, local_pgam, sc) do k = 1, nz @@ -1151,8 +1150,8 @@ subroutine Conservative2PrimitivePolytypeAM(CCTK_ARGUMENTS) if ((atmosphere_mask(i,j,k) .ne. 0) .or. & GRHydro_enable_internal_excision /= 0 .and. (hydro_excision_mask(i,j,k) .ne. 0)) cycle - det = SPATIAL_DETERMINANT(g11(i,j,k),g12(i,j,k),g13(i,j,k),g22(i,j,k),g23(i,j,k),g33(i,j,k)) - call UpperMetric(uxx,uxy,uxz,uyy,uyz,uzz,det,& + sdet = sdetg(i,j,k) + call UpperMetric(uxx,uxy,uxz,uyy,uyz,uzz,sdet*sdet,& g11(i,j,k),g12(i,j,k),g13(i,j,k),g22(i,j,k),& g23(i,j,k),g33(i,j,k)) @@ -1184,9 +1183,9 @@ subroutine Conservative2PrimitivePolytypeAM(CCTK_ARGUMENTS) xrho,1.0d0,xtemp,xye,xpress,keyerr,anyerr) local_pgam=log(xpress/local_K)/log(xrho) sc = local_K*dens(i,j,k) - Bconsx_tmp = sqrt(det)*Bprim(i,j,k,1) - Bconsy_tmp = sqrt(det)*Bprim(i,j,k,2) - Bconsz_tmp = sqrt(det)*Bprim(i,j,k,3) + Bconsx_tmp = sdet*Bprim(i,j,k,1) + Bconsy_tmp = sdet*Bprim(i,j,k,2) + Bconsz_tmp = sdet*Bprim(i,j,k,3) call GRHydro_Con2PrimM_Polytype_pt(GRHydro_eos_handle, local_pgam, dens(i,j,k), & scon(i,j,k,1),scon(i,j,k,2),scon(i,j,k,3), sc, & @@ -1194,7 +1193,7 @@ subroutine Conservative2PrimitivePolytypeAM(CCTK_ARGUMENTS) vup(i,j,k,1),vup(i,j,k,2),vup(i,j,k,3),eps(i,j,k),press(i,j,k),& Bprim(i,j,k,1), Bprim(i,j,k,2), Bprim(i,j,k,3),b2, w_lorentz(i,j,k),& g11(i,j,k),g12(i,j,k),g13(i,j,k),g22(i,j,k),g23(i,j,k),g33(i,j,k), & - uxx,uxy,uxz,uyy,uyz,uzz,det, & + uxx,uxy,uxz,uyy,uyz,uzz,sdet, & epsnegative,GRHydro_C2P_failed(i,j,k)) end do @@ -1254,8 +1253,8 @@ subroutine Con2PrimBoundsPolytypeAM(CCTK_ARGUMENTS) integer :: i, j, k, nx, ny, nz CCTK_REAL :: uxxl, uxyl, uxzl, uyyl, uyzl, uzzl,& uxxr, uxyr, uxzr, uyyr, uyzr, uzzr - CCTK_REAL :: gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_detl,& - gxxr,gxyr,gxzr,gyyr,gyzr,gzzr,avg_detr + CCTK_REAL :: gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_sdetl,& + gxxr,gxyr,gxzr,gyyr,gyzr,gzzr,avg_sdetr CCTK_REAL :: b2minus, b2plus CCTK_REAL :: Bconsx_tmp, Bconsy_tmp, Bconsz_tmp CCTK_INT :: epsnegative @@ -1328,11 +1327,11 @@ subroutine Con2PrimBoundsPolytypeAM(CCTK_ARGUMENTS) gyzr = 0.5d0 * (g23(i,j,k) + g23(i+xoffset,j+yoffset,k+zoffset)) gzzr = 0.5d0 * (g33(i,j,k) + g33(i+xoffset,j+yoffset,k+zoffset)) - avg_detl = SPATIAL_DETERMINANT(gxxl,gxyl,gxzl,gyyl, gyzl,gzzl) - avg_detr = SPATIAL_DETERMINANT(gxxr,gxyr,gxzr,gyyr, gyzr,gzzr) - call UpperMetric(uxxl,uxyl,uxzl,uyyl,uyzl,uzzl,avg_detl,& + avg_sdetl = sqrt(SPATIAL_DETERMINANT(gxxl,gxyl,gxzl,gyyl, gyzl,gzzl)) + avg_sdetr = sqrt(SPATIAL_DETERMINANT(gxxr,gxyr,gxzr,gyyr, gyzr,gzzr)) + call UpperMetric(uxxl,uxyl,uxzl,uyyl,uyzl,uzzl,avg_sdetl*avg_sdetl,& gxxl,gxyl,gxzl,gyyl,gyzl,gzzl) - call UpperMetric(uxxr,uxyr,uxzr,uyyr,uyzr,uzzr,avg_detr,& + call UpperMetric(uxxr,uxyr,uxzr,uyyr,uyzr,uzzr,avg_sdetr*avg_sdetr,& gxxr,gxyr,gxzr,gyyr,gyzr,gzzr) xrho=10.0d0 @@ -1345,9 +1344,9 @@ subroutine Con2PrimBoundsPolytypeAM(CCTK_ARGUMENTS) local_pgam=log(xpress/local_K)/log(xrho) scminus = local_K*densminus(i,j,k) scplus = local_K*densplus(i,j,k) - Bconsx_tmp = sqrt(avg_detl)*Bvecxminus(i,j,k) - Bconsy_tmp = sqrt(avg_detl)*Bvecyminus(i,j,k) - Bconsz_tmp = sqrt(avg_detl)*Bveczminus(i,j,k) + Bconsx_tmp = avg_sdetl*Bvecxminus(i,j,k) + Bconsy_tmp = avg_sdetl*Bvecyminus(i,j,k) + Bconsz_tmp = avg_sdetl*Bveczminus(i,j,k) call GRHydro_Con2PrimM_Polytype_pt(GRHydro_eos_handle, local_pgam,densminus(i,j,k), & sxminus(i,j,k),syminus(i,j,k),szminus(i,j,k), scminus,& @@ -1355,19 +1354,19 @@ subroutine Con2PrimBoundsPolytypeAM(CCTK_ARGUMENTS) velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k),epsminus(i,j,k),pressminus(i,j,k),& Bvecxminus(i,j,k), Bvecyminus(i,j,k), Bveczminus(i,j,k),b2minus, w_lorentzminus(i,j,k),& gxxl,gxyl,gxzl,gyyl,gyzl,gzzl, & - uxxl,uxyl,uxzl,uyyl,uyzl,uzzl,avg_detl, & + uxxl,uxyl,uxzl,uyyl,uyzl,uzzl,avg_sdetl, & epsnegative,GRHydro_C2P_failed(i,j,k)) - Bconsx_tmp = sqrt(avg_detr)*Bvecxplus(i,j,k) - Bconsy_tmp = sqrt(avg_detr)*Bvecyplus(i,j,k) - Bconsz_tmp = sqrt(avg_detr)*Bveczplus(i,j,k) + Bconsx_tmp = avg_sdetr*Bvecxplus(i,j,k) + Bconsy_tmp = avg_sdetr*Bvecyplus(i,j,k) + Bconsz_tmp = avg_sdetr*Bveczplus(i,j,k) call GRHydro_Con2PrimM_Polytype_pt(GRHydro_eos_handle, local_pgam,densplus(i,j,k), & sxplus(i,j,k),syplus(i,j,k),szplus(i,j,k), scplus,& Bconsx_tmp, Bconsy_tmp, Bconsz_tmp, rhoplus(i,j,k),& velxplus(i,j,k),velyplus(i,j,k),velzplus(i,j,k),epsplus(i,j,k),pressplus(i,j,k),& Bvecxplus(i,j,k), Bvecyplus(i,j,k),Bveczplus(i,j,k),b2plus,w_lorentzplus(i,j,k),& gxxr,gxyr,gxzr,gyyr,gyzr,gzzr, & - uxxr,uxyr,uxzr,uyyr,uyzr,uzzr,avg_detr, & + uxxr,uxyr,uxzr,uyyr,uyzr,uzzr,avg_sdetr, & epsnegative,GRHydro_C2P_failed(i,j,k)) end do end do |