diff options
Diffstat (limited to 'src/GRHydro_Con2PrimM.F90')
| -rw-r--r-- | src/GRHydro_Con2PrimM.F90 | 79 |
1 files changed, 39 insertions, 40 deletions
diff --git a/src/GRHydro_Con2PrimM.F90 b/src/GRHydro_Con2PrimM.F90 index 6517d41..a871b7e 100644 --- a/src/GRHydro_Con2PrimM.F90 +++ b/src/GRHydro_Con2PrimM.F90 @@ -60,7 +60,7 @@ subroutine Conservative2PrimitiveM(CCTK_ARGUMENTS) DECLARE_CCTK_FUNCTIONS integer :: i, j, k, itracer, nx, ny, nz - CCTK_REAL :: uxx, uxy, uxz, uyy, uyz, uzz, det, sdet, pmin(1), epsmin(1) + CCTK_REAL :: uxx, uxy, uxz, uyy, uyz, uzz, sdet, pmin(1), epsmin(1) CCTK_REAL :: oob, b2, d2, s2, bscon, bxhat, byhat, bzhat, bhatscon CCTK_REAL :: Wm, Wm0, Wm_resid, Wmold CCTK_REAL :: s2m, s2m0, s2m_resid, s2mold, s2max, taum, dummy1, dummy2 @@ -161,7 +161,7 @@ subroutine Conservative2PrimitiveM(CCTK_ARGUMENTS) local_gam = local_gam+1.d0 !$OMP PARALLEL DO PRIVATE(i,j,k,itracer,& - !$OMP uxx, uxy, uxz, uyy, uyz, uzz, det, epsnegative, & + !$OMP uxx, uxy, uxz, uyy, uyz, uzz, epsnegative, & !$OMP b2,xrho,xeps,xpress,xtemp,local_K,local_pgam,sc,keyerr,anyerr,keytemp, & !$OMP local_perc_ptol,posdef,g11c,g12c,g13c,g22c,g23c,g33c,tmp1, & !$OMP sdet,d2,s2,oob,bscon,bxhat,byhat,bzhat, & @@ -181,10 +181,9 @@ subroutine Conservative2PrimitiveM(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)) @@ -206,7 +205,7 @@ subroutine Conservative2PrimitiveM(CCTK_ARGUMENTS) ! In excised region, set to atmosphere! if (GRHydro_enable_internal_excision /= 0 .and. (hydro_excision_mask(i,j,k) .gt. 0)) then - SET_ATMO_MIN(dens(i,j,k), sqrt(det)*GRHydro_rho_min, r(i,j,k)) !sqrt(det)*GRHydro_rho_min !/(1.d0+GRHydro_atmo_tolerance) + SET_ATMO_MIN(dens(i,j,k), sdet*GRHydro_rho_min, r(i,j,k)) !sqrt(det)*GRHydro_rho_min !/(1.d0+GRHydro_atmo_tolerance) SET_ATMO_MIN(rho(i,j,k), GRHydro_rho_min, r(i,j,k)) !GRHydro_rho_min scon(i,j,k,:) = 0.d0 vup(i,j,k,:) = 0.d0 @@ -386,7 +385,7 @@ subroutine Conservative2PrimitiveM(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 @@ -397,7 +396,7 @@ subroutine Conservative2PrimitiveM(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 @@ -443,13 +442,13 @@ subroutine Conservative2PrimitiveM(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)) if(GRHydro_C2P_failed(i,j,k).ne.0) then GRHydro_C2P_failed(i,j,k) = 0 call prim2conM(GRHydro_eos_handle,g11(i,j,k),g12(i,j,k), & - g13(i,j,k),g22(i,j,k),g23(i,j,k),g33(i,j,k),det, & + g13(i,j,k),g22(i,j,k),g23(i,j,k),g33(i,j,k),sdet, & dens(i,j,k),scon(i,j,k,1),scon(i,j,k,2),scon(i,j,k,3), & tau(i,j,k),Bcons(i,j,k,1),Bcons(i,j,k,2),Bcons(i,j,k,3), & rho(i,j,k),vup(i,j,k,1),vup(i,j,k,2),vup(i,j,k,3), & @@ -501,7 +500,7 @@ subroutine Conservative2PrimitiveM(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, & @@ -510,7 +509,7 @@ subroutine Conservative2PrimitiveM(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 @@ -550,7 +549,7 @@ subroutine Conservative2PrimitiveM(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. @@ -576,7 +575,7 @@ subroutine Conservative2PrimitiveM(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 @@ -621,7 +620,7 @@ subroutine Conservative2PrimitiveM(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 @@ -639,7 +638,7 @@ subroutine Conservative2PrimitiveM(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 @@ -677,7 +676,7 @@ subroutine Conservative2PrimitiveM(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 @@ -793,8 +792,8 @@ subroutine Conservative2PrimitiveBoundsM(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 character(len=100) warnline @@ -887,11 +886,11 @@ subroutine Conservative2PrimitiveBoundsM(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! @@ -919,7 +918,7 @@ subroutine Conservative2PrimitiveBoundsM(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)) if (epsnegative .ne. 0) then @@ -943,7 +942,7 @@ subroutine Conservative2PrimitiveBoundsM(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 @@ -977,7 +976,7 @@ subroutine Conservative2PrimitiveBoundsM(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)) if (epsnegative .ne. 0) then @@ -1001,7 +1000,7 @@ subroutine Conservative2PrimitiveBoundsM(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 @@ -1076,7 +1075,7 @@ subroutine Conservative2PrimitivePolytypeM(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 @@ -1144,7 +1143,7 @@ subroutine Conservative2PrimitivePolytypeM(CCTK_ARGUMENTS) !!$ do i = GRHydro_stencil + 1, nx - GRHydro_stencil !$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 b2, xrho, xpress, local_K, local_pgam, sc) do k = 1, nz do j = 1, ny @@ -1154,8 +1153,8 @@ subroutine Conservative2PrimitivePolytypeM(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)) @@ -1194,7 +1193,7 @@ subroutine Conservative2PrimitivePolytypeM(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 Con2PrimBoundsPolytypeM(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_INT :: epsnegative @@ -1327,11 +1326,11 @@ subroutine Con2PrimBoundsPolytypeM(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 @@ -1351,7 +1350,7 @@ subroutine Con2PrimBoundsPolytypeM(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)) call GRHydro_Con2PrimM_Polytype_pt(GRHydro_eos_handle, local_pgam,densplus(i,j,k), & @@ -1360,7 +1359,7 @@ subroutine Con2PrimBoundsPolytypeM(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 do end do |