diff options
| author | rhaas <rhaas@c83d129a-5a75-4d5a-9c4d-ed3a5855bf45> | 2013-07-06 18:11:12 +0000 |
|---|---|---|
| committer | rhaas <rhaas@c83d129a-5a75-4d5a-9c4d-ed3a5855bf45> | 2013-07-06 18:11:12 +0000 |
| commit | 669d5b37c19d16380ba6c63a7195cbcea8d08f35 (patch) | |
| tree | 43c35e9b941be21d80ba34711d3bbfe7825b3b6d | |
| parent | fcbaf698714e9594ef35fbba18530a4a1ab164ba (diff) | |
GRHydro: add grid function for sqrt(detg)
* add new 1 tl grid function sdetg that stores the sqrt of the determinent
of the 3-metric.
* replace lots of re-computation of det by use of this grid function
From: Christian Ott <cott@bethe.tapir.caltech.edu>
git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEvolve/GRHydro/trunk@555 c83d129a-5a75-4d5a-9c4d-ed3a5855bf45
| -rw-r--r-- | interface.ccl | 8 | ||||
| -rw-r--r-- | schedule.ccl | 20 | ||||
| -rw-r--r-- | src/GRHydro_BvecfromAvec.F90 | 8 | ||||
| -rw-r--r-- | src/GRHydro_Con2Prim.F90 | 117 | ||||
| -rw-r--r-- | src/GRHydro_Con2PrimAM.F90 | 99 | ||||
| -rw-r--r-- | src/GRHydro_Con2PrimHot.F90 | 46 | ||||
| -rw-r--r-- | src/GRHydro_Con2PrimM.F90 | 79 | ||||
| -rw-r--r-- | src/GRHydro_Con2PrimM_polytype_pt.c | 7 | ||||
| -rw-r--r-- | src/GRHydro_Con2PrimM_pt.c | 9 | ||||
| -rw-r--r-- | src/GRHydro_Con2PrimM_pt_EOSOmni.c | 13 | ||||
| -rw-r--r-- | src/GRHydro_Con2PrimM_pt_EOSOmni.c.hopefullyfixed.c | 9 | ||||
| -rw-r--r-- | src/GRHydro_Con2PrimM_pt_EOSOmniold.c | 9 | ||||
| -rw-r--r-- | src/GRHydro_Con2PrimM_ptee.c | 9 | ||||
| -rw-r--r-- | src/GRHydro_EoSChangeGamma.F90 | 15 | ||||
| -rw-r--r-- | src/GRHydro_FluxSplit.F90 | 26 | ||||
| -rw-r--r-- | src/GRHydro_Macros.h | 2 | ||||
| -rw-r--r-- | src/GRHydro_Prim2Con.F90 | 133 | ||||
| -rw-r--r-- | src/GRHydro_Prim2ConAM.F90 | 114 | ||||
| -rw-r--r-- | src/GRHydro_Prim2ConM.F90 | 138 | ||||
| -rw-r--r-- | src/GRHydro_Source.F90 | 18 | ||||
| -rw-r--r-- | src/GRHydro_SourceAM.F90 | 10 | ||||
| -rw-r--r-- | src/GRHydro_SourceM.F90 | 10 | ||||
| -rw-r--r-- | src/GRHydro_UpdateMask.F90 | 38 | ||||
| -rw-r--r-- | src/GRHydro_UpdateMaskM.F90 | 76 | ||||
| -rw-r--r-- | src/Utils.F90 | 88 |
25 files changed, 517 insertions, 584 deletions
diff --git a/interface.ccl b/interface.ccl index 6fa145d..4659e4b 100644 --- a/interface.ccl +++ b/interface.ccl @@ -423,6 +423,8 @@ real GRHydro_tracers[number_of_tracers] type = GF Timelevels = 3 tags='Prolongat #real w_lorentz type = GF Timelevels = 3 tags='ProlongationParameter="HydroBase::prolongation_type" tensortypealias="Scalar" tensorweight=+1.0 interpolator="matter"' "Lorentz factor" +real sdetg type = GF Timelevels = 1 tags='Prolongation="None" tensortypealias="Scalar" tensorweight=+1.0 interpolator="matter" checkpoint="no"' "Sqrt of the determinant of the 3-metric" + real psidc type = GF Timelevels = 3 tags='ProlongationParameter="HydroBase::prolongation_type" tensortypealias="Scalar" tensorweight=+1.0 tensorparity=-1 jacobian="inverse_jacobian" interpolator="matter"' "Psi parameter for divergence cleaning" real densrhs type = GF Timelevels = 1 tags='Prolongation="None" checkpoint="no"' "Update term for dens" @@ -775,12 +777,6 @@ CCTK_REAL RoeAverage_temps TYPE=GF tags='Prolongation="None" checkpoint="no"' eos_cs2_ave, eos_dpdeps_ave } "Temporaries for the Roe solver" -CCTK_REAL Metric_temps TYPE=GF tags='Prolongation="None" checkpoint="no"' -{ - GRHydro_det - GRHydro_uxx,GRHydro_uxy,GRHydro_uxz,GRHydro_uyy,GRHydro_uyz,GRHydro_uzz -} "Temporaries for the determinant of the 3-metric and the inverse metric" - CCTK_REAL Con2Prim_temps TYPE=GF tags='Prolongation="None" checkpoint="no"' { press_old, press_new diff --git a/schedule.ccl b/schedule.ccl index e01be59..bf4ec6e 100644 --- a/schedule.ccl +++ b/schedule.ccl @@ -88,6 +88,7 @@ STORAGE:evolve_Y_e STORAGE:evolve_temper STORAGE:GRHydro_reflevel STORAGE:InLastMoLPostStep +STORAGE:sdetg STORAGE:densrhs STORAGE:taurhs STORAGE:srhs @@ -1062,6 +1063,25 @@ schedule GRHydro_RefinementLevel IN FluxTerms BEFORE Reconstruct LANG: Fortran } "Calculate current refinement level" +########################### Sqrt(detg) ################# + +# recompute root of metric determinant whenever metric is updated we could alse +# abuse ADMBase_SetADMVars the way GRHydroTransformADMToLocalBasis does. Abuse +# since the group description says to "Set the ADM variables before this group, +# and use them afterwards" it does not say to schedule anything in it. +# If we ever face trouble with the current placement we will ignore the +# description and do so anyway. +schedule GRHydro_SqrtSpatialDeterminant IN HydroBase_Con2Prim BEFORE Con2Prim IF GRHydro::execute_MoL_Step +{ + LANG: Fortran +} "Calculate sdetg" + +schedule GRHydro_SqrtSpatialDeterminant AT CCTK_INITIAL AFTER (HydroBase_Initial,GRHydroTransformADMToLocalBasis) BEFORE HydroBase_Prim2ConInitial +{ + LANG: Fortran +} "Calculate sdetg" + + #debug ###schedule GRHydro_Debug IN HydroBase_PostStep ###{ diff --git a/src/GRHydro_BvecfromAvec.F90 b/src/GRHydro_BvecfromAvec.F90 index ea0a58d..30da906 100644 --- a/src/GRHydro_BvecfromAvec.F90 +++ b/src/GRHydro_BvecfromAvec.F90 @@ -53,7 +53,7 @@ subroutine GRHydro_BvecfromAvec(CCTK_ARGUMENTS) implicit none CCTK_REAL :: Az_y, Ay_z, Ax_z, Az_x, Ay_x, Ax_y - CCTK_REAL :: det, isdet + CCTK_REAL :: sdet, isdet CCTK_REAL :: dx, dy, dz, idx, idy, idz CCTK_INT :: i, j, k, nx, ny, nz, GRHydro_UseGeneralCoordinates @@ -102,7 +102,7 @@ subroutine GRHydro_BvecfromAvec(CCTK_ARGUMENTS) call CCTK_WARN(1,"Bvec from Avec start."); - !$OMP PARALLEL DO PRIVATE( det, isdet, local_spatial_order, & + !$OMP PARALLEL DO PRIVATE( sdet, isdet, local_spatial_order, & !$OMP Az_y, Ay_z, Ax_z, Az_x, Ay_x, Ax_y) do k = GRHydro_stencil, cctk_lsh(3)-GRHydro_stencil+1 do j = GRHydro_stencil, cctk_lsh(2)-GRHydro_stencil+1 @@ -129,8 +129,8 @@ subroutine GRHydro_BvecfromAvec(CCTK_ARGUMENTS) Ax_y = DIFF_Y_4(Avecx) end if - det = SPATIAL_DETERMINANT(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k)) - isdet = 1.d0/sqrt(det) + sdet = sdetg(i,j,k) + isdet = 1.d0/sdet Bvecx(i,j,k) = isdet*( Az_y - Ay_z ) Bvecy(i,j,k) = isdet*( Ax_z - Az_x ) diff --git a/src/GRHydro_Con2Prim.F90 b/src/GRHydro_Con2Prim.F90 index 9d33962..176ba77 100644 --- a/src/GRHydro_Con2Prim.F90 +++ b/src/GRHydro_Con2Prim.F90 @@ -54,7 +54,7 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS) DECLARE_CCTK_PARAMETERS integer :: i, j, k, itracer, nx, ny, nz - CCTK_REAL :: uxx, uxy, uxz, uyy, uyz, uzz, det, pmin, epsmin, dummy1, dummy2 + CCTK_REAL :: uxx, uxy, uxz, uyy, uyz, uzz, pmin, epsmin, dummy1, dummy2 logical :: epsnegative character*256 :: warnline @@ -115,7 +115,7 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS) GRHydro_rho_min,xeps,xtemp,xye,pmin,epsmin,keyerr,anyerr) !$OMP PARALLEL DO PRIVATE(i,j,k,itracer,& - !$OMP uxx, uxy, uxz, uyy, uyz, uzz, det, epsnegative, anyerr, keyerr, keytemp,& + !$OMP uxx, uxy, uxz, uyy, uyz, uzz, epsnegative, anyerr, keyerr, keytemp,& !$OMP warnline, dummy1, dummy2) do k = 1, nz do j = 1, ny @@ -160,15 +160,13 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS) epsnegative = .false. - 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,& + call UpperMetric(uxx,uxy,uxz,uyy,uyz,uzz,sdetg(i,j,k)*sdetg(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)) ! 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), sdetg(i,j,k)*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 @@ -193,7 +191,7 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS) endif ! w_lorentz=1, so the expression for tau reduces to: - tau(i,j,k) = sqrt(det) * (rho(i,j,k)+rho(i,j,k)*eps(i,j,k)) - dens(i,j,k) + tau(i,j,k) = sdetg(i,j,k) * (rho(i,j,k)+rho(i,j,k)*eps(i,j,k)) - dens(i,j,k) cycle endif @@ -223,8 +221,8 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS) endif !if ( dens(i,j,k) .le. sqrt(det)*GRHydro_rho_min*(1.d0+GRHydro_atmo_tolerance) ) then - IF_BELOW_ATMO(dens(i,j,k), sqrt(det)*GRHydro_rho_min, GRHydro_atmo_tolerance, r(i,j,k)) 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) + IF_BELOW_ATMO(dens(i,j,k), sdetg(i,j,k)*GRHydro_rho_min, GRHydro_atmo_tolerance, r(i,j,k)) then + SET_ATMO_MIN(dens(i,j,k), sdetg(i,j,k)*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 @@ -249,7 +247,7 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS) endif ! w_lorentz=1, so the expression for tau reduces to: - tau(i,j,k) = sqrt(det) * (rho(i,j,k)+rho(i,j,k)*eps(i,j,k)) - dens(i,j,k) + tau(i,j,k) = sdetg(i,j,k) * (rho(i,j,k)+rho(i,j,k)*eps(i,j,k)) - dens(i,j,k) cycle @@ -260,7 +258,7 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS) GRHydro_eos_handle, dens(i,j,k),scon(i,j,k,1),scon(i,j,k,2), & scon(i,j,k,3),tau(i,j,k),rho(i,j,k),vup(i,j,k,1),vup(i,j,k,2), & vup(i,j,k,3),eps(i,j,k),press(i,j,k),w_lorentz(i,j,k), & - uxx,uxy,uxz,uyy,uyz,uzz,det,x(i,j,k),y(i,j,k), & + uxx,uxy,uxz,uyy,uyz,uzz,sdetg(i,j,k),x(i,j,k),y(i,j,k), & z(i,j,k),r(i,j,k),epsnegative,GRHydro_rho_min,pmin, epsmin, & GRHydro_reflevel, GRHydro_C2P_failed(i,j,k)) else @@ -292,7 +290,7 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS) !$OMP END CRITICAL ! for safety, let's set the point to atmosphere - dens(i,j,k) = ATMOMIN(sqrt(det)*GRHydro_rho_min, r(i,j,k)) !sqrt(det)*GRHydro_rho_min !/(1.d0+GRHydro_atmo_tolerance) + dens(i,j,k) = ATMOMIN(sdetg(i,j,k)*GRHydro_rho_min, r(i,j,k)) !sqrt(det)*GRHydro_rho_min !/(1.d0+GRHydro_atmo_tolerance) rho(i,j,k) = ATMOMIN(GRHydro_rho_min, r(i,j,k)) !GRHydro_rho_min scon(i,j,k,:) = 0.d0 vup(i,j,k,:) = 0.d0 @@ -305,7 +303,7 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS) rho(i,j,k),eps(i,j,k),xtemp,xye,press(i,j,k),eps(i,j,k),keyerr,anyerr) ! w_lorentz=1, so the expression for tau reduces to: - tau(i,j,k) = sqrt(det) * (rho(i,j,k)+rho(i,j,k)*eps(i,j,k)) - dens(i,j,k) + tau(i,j,k) = sdetg(i,j,k) * (rho(i,j,k)+rho(i,j,k)*eps(i,j,k)) - dens(i,j,k) #else ! cott 2010/03/27: @@ -319,7 +317,7 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS) dens(i,j,k), scon(i,j,k,1), scon(i,j,k,2), scon(i,j,k,3), & tau(i,j,k), rho(i,j,k), vup(i,j,k,1), vup(i,j,k,2), & vup(i,j,k,3), eps(i,j,k), press(i,j,k), w_lorentz(i,j,k), & - uxx, uxy, uxz, uyy, uyz, uzz, det, x(i,j,k), y(i,j,k), & + uxx, uxy, uxz, uyy, uyz, uzz, sdetg(i,j,k), x(i,j,k), y(i,j,k), & z(i,j,k), r(i,j,k),GRHydro_rho_min, GRHydro_reflevel, GRHydro_C2P_failed(i,j,k)) #endif @@ -393,7 +391,7 @@ a single point. subroutine Con2Prim_pt(cctk_iteration,ii,jj,kk,& handle, dens, sx, sy, sz, tau, rho, velx, vely, & velz, epsilon, press, w_lorentz, uxx, uxy, uxz, uyy, & - uyz, uzz, det, x, y, z, r, epsnegative, GRHydro_rho_min, pmin, epsmin, & + uyz, uzz, sdetg, x, y, z, r, epsnegative, GRHydro_rho_min, pmin, epsmin, & GRHydro_reflevel, GRHydro_C2P_failed) implicit none @@ -401,7 +399,7 @@ subroutine Con2Prim_pt(cctk_iteration,ii,jj,kk,& DECLARE_CCTK_PARAMETERS CCTK_REAL dens, sx, sy, sz, tau, rho, velx, vely, velz, epsilon, & - press, uxx, uxy, uxz, uyy, uyz, uzz, det, w_lorentz, x, & + press, uxx, uxy, uxz, uyy, uyz, uzz, sdetg, isdetg, w_lorentz, x, & y, z, r, GRHydro_rho_min CCTK_REAL s2, f, df, vlowx, vlowy, vlowz @@ -424,11 +422,12 @@ subroutine Con2Prim_pt(cctk_iteration,ii,jj,kk,& !!$ Undensitize the variables - udens = dens /sqrt(det) - usx = sx /sqrt(det) - usy = sy /sqrt(det) - usz = sz /sqrt(det) - utau = tau /sqrt(det) + isdetg = 1.0d0/sdetg + udens = dens * isdetg + usx = sx * isdetg + usy = sy * isdetg + usz = sz * isdetg + utau = tau * isdetg s2 = usx*usx*uxx + usy*usy*uyy + usz*usz*uzz + 2.*usx*usy*uxy + & 2.*usx*usz*uxz + 2.*usy*usz*uyz @@ -512,7 +511,7 @@ subroutine Con2Prim_pt(cctk_iteration,ii,jj,kk,& ! for safety, let's set the point to atmosphere SET_ATMO_MIN(rho, GRHydro_rho_min, r) udens = rho - dens = sqrt(det) * rho + dens = sdetg * rho pnew = pmin epsilon = epsmin ! w_lorentz=1, so the expression for utau reduces to: @@ -685,7 +684,7 @@ subroutine Con2Prim_pt(cctk_iteration,ii,jj,kk,& IF_BELOW_ATMO(rho, GRHydro_rho_min, GRHydro_atmo_tolerance, r) then SET_ATMO_MIN(rho, GRHydro_rho_min, r) !GRHydro_rho_min udens = rho - dens = sqrt(det) * rho + dens = isdetg * rho ! epsilon = (sqrt( (utau + pnew + udens)**2) - pnew - udens)/udens epsilon = epsmin ! w_lorentz=1, so the expression for utau reduces to: @@ -874,11 +873,11 @@ subroutine Conservative2PrimitiveBounds(CCTK_ARGUMENTS) epsnegative = .false. - 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_detl = sqrt(SPATIAL_DETERMINANT(gxxl,gxyl,gxzl,gyyl, gyzl,gzzl)) + avg_detr = sqrt(SPATIAL_DETERMINANT(gxxr,gxyr,gxzr,gyyr, gyzr,gzzr)) + call UpperMetric(uxxl,uxyl,uxzl,uyyl,uyzl,uzzl,avg_detl*avg_detl,& 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_detr*avg_detr,& gxxr,gxyr,gxzr,gyyr,gyzr,gzzr) if (evolve_tracer .ne. 0) then @@ -1035,7 +1034,7 @@ subroutine Conservative2PrimitivePolytype(CCTK_ARGUMENTS) DECLARE_CCTK_PARAMETERS integer :: i, j, k, itracer, nx, ny, nz - CCTK_REAL :: uxx, uxy, uxz, uyy, uyz, uzz, det + CCTK_REAL :: uxx, uxy, uxz, uyy, uyz, uzz CCTK_REAL :: local_min_tracer ! character(len=400) :: warnline @@ -1080,7 +1079,7 @@ subroutine Conservative2PrimitivePolytype(CCTK_ARGUMENTS) !!$ do j = GRHydro_stencil + 1, ny - GRHydro_stencil !!$ do i = GRHydro_stencil + 1, nx - GRHydro_stencil !$OMP PARALLEL DO PRIVATE(i,j,k,itracer,& - !$OMP uxx, uxy, uxz, uyy, uyz, uzz, det) + !$OMP uxx, uxy, uxz, uyy, uyz, uzz) do k = 1, nz do j = 1, ny do i = 1, nx @@ -1089,9 +1088,7 @@ subroutine Conservative2PrimitivePolytype(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,& + call UpperMetric(uxx,uxy,uxz,uyy,uyz,uzz,sdetg(i,j,k)*sdetg(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)) @@ -1117,7 +1114,7 @@ subroutine Conservative2PrimitivePolytype(CCTK_ARGUMENTS) scon(i,j,k,1),scon(i,j,k,2), & scon(i,j,k,3),tau(i,j,k),rho(i,j,k),vup(i,j,k,1),vup(i,j,k,2), & vup(i,j,k,3),eps(i,j,k),press(i,j,k),w_lorentz(i,j,k), & - uxx,uxy,uxz,uyy,uyz,uzz,det,x(i,j,k),y(i,j,k), & + uxx,uxy,uxz,uyy,uyz,uzz,sdetg(i,j,k),x(i,j,k),y(i,j,k), & z(i,j,k),r(i,j,k),GRHydro_rho_min, GRHydro_reflevel, GRHydro_C2P_failed(i,j,k)) end do @@ -1149,20 +1146,20 @@ a single point. subroutine Con2Prim_ptPolytype(handle, dens, sx, sy, sz, tau, rho, & velx, vely, velz, epsilon, press, w_lorentz, uxx, uxy, uxz, uyy, & - uyz, uzz, det, x, y, z, r, GRHydro_rho_min, GRHydro_reflevel, GRHydro_C2P_failed) + uyz, uzz, sdetg, x, y, z, r, GRHydro_rho_min, GRHydro_reflevel, GRHydro_C2P_failed) implicit none DECLARE_CCTK_PARAMETERS CCTK_REAL dens, sx, sy, sz, tau, rho, velx, vely, velz, epsilon, & - press, uxx, uxy, uxz, uyy, uyz, uzz, det, w_lorentz, x, & + press, uxx, uxy, uxz, uyy, uyz, uzz, sdetg, isdetg, w_lorentz, x, & y, z, r, GRHydro_rho_min CCTK_REAL s2, f, df, vlowx, vlowy, vlowz CCTK_INT count, handle, GRHydro_reflevel CCTK_REAL udens, usx, usy, usz, rhoold, rhonew, & - enthalpy, denthalpy, sqrtdet, invsqrtdet, invfac, GRHydro_C2P_failed, dummy1, dummy2 + enthalpy, denthalpy, invfac, GRHydro_C2P_failed, dummy1, dummy2 character(len=200) warnline ! begin EOS Omni vars @@ -1175,12 +1172,11 @@ subroutine Con2Prim_ptPolytype(handle, dens, sx, sy, sz, tau, rho, & !!$ Undensitize the variables - sqrtdet = sqrt(det) - invsqrtdet = 1.d0/sqrtdet - udens = dens*invsqrtdet - usx = sx*invsqrtdet - usy = sy*invsqrtdet - usz = sz*invsqrtdet + isdetg = 1.0d0/sdetg + udens = dens*isdetg + usx = sx*isdetg + usy = sy*isdetg + usz = sz*isdetg s2 = usx*usx*uxx + usy*usy*uyy + usz*usz*uzz + 2.d0*usx*usy*uxy + & 2.d0*usx*usz*uxz + 2.d0*usy*usz*uyz @@ -1231,7 +1227,7 @@ subroutine Con2Prim_ptPolytype(handle, dens, sx, sy, sz, tau, rho, & ! for safety, let's set the point to atmosphere SET_ATMO_MIN(rhonew, GRHydro_rho_min, r) udens = rhonew - dens = sqrtdet * rhonew + dens = sdetg * rhonew call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,& rhonew,1.0d0,xtemp,xye,press,keyerr,anyerr) @@ -1286,7 +1282,7 @@ subroutine Con2Prim_ptPolytype(handle, dens, sx, sy, sz, tau, rho, & enthalpy = 1.0d0 + epsilon + press / rhonew w_lorentz = sqrt(1.d0 + s2 / ( (udens*enthalpy)**2 )) - tau = sqrtdet * ((rho * enthalpy) * w_lorentz**2 - press) - dens + tau = sdetg * ((rho * enthalpy) * w_lorentz**2 - press) - dens f = rhonew * w_lorentz - udens @@ -1302,7 +1298,7 @@ subroutine Con2Prim_ptPolytype(handle, dens, sx, sy, sz, tau, rho, & ! before 2009/10/07 dens was not reset and was used with its negative ! value below; this has since been changed, but the change produces ! significant changes in the results - dens = sqrtdet * rhonew + dens = sdetg * rhonew call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,& rhonew,1.0d0,xtemp,xye,press,keyerr,anyerr) @@ -1338,7 +1334,7 @@ subroutine Con2Prim_ptPolytype(handle, dens, sx, sy, sz, tau, rho, & call EOS_Omni_EpsFromPress(handle,keytemp,GRHydro_eos_rf_prec,n,& rhonew,xeps,xtemp,xye,press,epsilon,keyerr,anyerr) - tau = sqrtdet * ((rho * enthalpy) * w_lorentz**2 - press) - dens + tau = sdetg * ((rho * enthalpy) * w_lorentz**2 - press) - dens if (epsilon .lt. 0.0d0) then GRHydro_C2P_failed = 1 @@ -1355,7 +1351,7 @@ subroutine Con2Prim_ptPolytype(handle, dens, sx, sy, sz, tau, rho, & !$OMP END CRITICAL rho = GRHydro_rho_min - dens = sqrtdet * rho + dens = sdetg * rho call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,& rhonew,1.d0,xtemp,xye,press,keyerr,anyerr) @@ -1363,7 +1359,7 @@ subroutine Con2Prim_ptPolytype(handle, dens, sx, sy, sz, tau, rho, & rhonew,xeps,xtemp,xye,press,epsilon,keyerr,anyerr) ! w_lorentz=1, so the expression for tau reduces to: - tau = sqrtdet * (rho+rho*epsilon) - dens + tau = sdetg * (rho+rho*epsilon) - dens usx = 0.d0 usy = 0.d0 usz = 0.d0 @@ -1469,11 +1465,11 @@ subroutine Con2PrimBoundsPolytype(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_detl = sqrt(SPATIAL_DETERMINANT(gxxl,gxyl,gxzl,gyyl, gyzl,gzzl)) + avg_detr = sqrt(SPATIAL_DETERMINANT(gxxr,gxyr,gxzr,gyyr, gyzr,gzzr)) + call UpperMetric(uxxl,uxyl,uxzl,uyyl,uyzl,uzzl,avg_detl*avg_detl,& 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_detr*avg_detr,& gxxr,gxyr,gxzr,gyyr,gyzr,gzzr) call Con2Prim_ptPolytype(GRHydro_eos_handle, densminus(i,j,k),& @@ -1484,7 +1480,8 @@ subroutine Con2PrimBoundsPolytype(CCTK_ARGUMENTS) uxxl,uxyl,uxzl,uyyl,uyzl,uzzl,avg_detl,& x(i,j,k)-0.5d0*CCTK_DELTA_SPACE(1),& y(i,j,k)-0.5d0*CCTK_DELTA_SPACE(2), & - z(i,j,k)-0.5d0*CCTK_DELTA_SPACE(3),r(i,j,k),GRHydro_rho_min, GRHydro_reflevel, GRHydro_C2P_failed(i,j,k)) + z(i,j,k)-0.5d0*CCTK_DELTA_SPACE(3),& + r(i,j,k),GRHydro_rho_min, GRHydro_reflevel, GRHydro_C2P_failed(i,j,k)) call Con2Prim_ptPolytype(GRHydro_eos_handle, densplus(i,j,k),& sxplus(i,j,k),syplus(i,j,k),szplus(i,j,k),& tauplus(i,j,k),rhoplus(i,j,k),velxplus(i,j,k),& @@ -1582,11 +1579,11 @@ subroutine Con2PrimBoundsTracer(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_detl = sqrt(SPATIAL_DETERMINANT(gxxl,gxyl,gxzl,gyyl, gyzl,gzzl)) + avg_detr = sqrt(SPATIAL_DETERMINANT(gxxr,gxyr,gxzr,gyyr, gyzr,gzzr)) + call UpperMetric(uxxl,uxyl,uxzl,uyyl,uyzl,uzzl,avg_detl*avg_detl,& 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_detr*avg_detr,& gxxr,gxyr,gxzr,gyyr,gyzr,gzzr) do itracer=1,number_of_tracers @@ -1681,15 +1678,15 @@ end subroutine Con2Prim_ptTracer @@*/ -subroutine Con2Prim_ptBoundsTracer(cons_tracer, tracer, rho, one_over_w_lorentz, det) +subroutine Con2Prim_ptBoundsTracer(cons_tracer, tracer, rho, one_over_w_lorentz, sdet) implicit none DECLARE_CCTK_PARAMETERS - CCTK_REAL cons_tracer, tracer, rho, one_over_w_lorentz, det + CCTK_REAL cons_tracer, tracer, rho, one_over_w_lorentz, sdet - tracer = cons_tracer / sqrt(det) / rho * one_over_w_lorentz + tracer = cons_tracer / sdet / rho * one_over_w_lorentz end subroutine Con2Prim_ptBoundsTracer 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 diff --git a/src/GRHydro_Con2PrimHot.F90 b/src/GRHydro_Con2PrimHot.F90 index 5725e18..7598e07 100644 --- a/src/GRHydro_Con2PrimHot.F90 +++ b/src/GRHydro_Con2PrimHot.F90 @@ -21,7 +21,7 @@ subroutine Conservative2PrimitiveHot(CCTK_ARGUMENTS) DECLARE_CCTK_FUNCTIONS integer :: i, j, k, itracer, nx, ny, nz, myproc - CCTK_REAL :: uxx, uxy, uxz, uyy, uyz, uzz, det, pmin, epsmin, dummy1, dummy2 + CCTK_REAL :: uxx, uxy, uxz, uyy, uyz, uzz, sdet, pmin, epsmin, dummy1, dummy2 CCTK_REAL :: vlowx, vlowy, vlowz logical :: epsnegative character*512 :: warnline @@ -82,7 +82,7 @@ subroutine Conservative2PrimitiveHot(CCTK_ARGUMENTS) epsmin = xeps(1) !$OMP PARALLEL DO PRIVATE(i,j,k,itracer,& - !$OMP uxx, uxy, uxz, uyy, uyz, uzz, det, epsnegative, anyerr, keyerr, keytemp,& + !$OMP uxx, uxy, uxz, uyy, uyz, uzz, sdet, epsnegative, anyerr, keyerr, keytemp,& !$OMP warnline, dummy1, dummy2,reset_to_atmo) do k = 1, nz do j = 1, ny @@ -94,9 +94,9 @@ subroutine Conservative2PrimitiveHot(CCTK_ARGUMENTS) epsnegative = .false. - 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)) @@ -121,12 +121,12 @@ subroutine Conservative2PrimitiveHot(CCTK_ARGUMENTS) endif reset_to_atmo = 0 - IF_BELOW_ATMO(dens(i,j,k), sqrt(det)*GRHydro_rho_min, GRHydro_atmo_tolerance, r(i,j,k)) then + IF_BELOW_ATMO(dens(i,j,k), sdet*GRHydro_rho_min, GRHydro_atmo_tolerance, r(i,j,k)) then reset_to_atmo = 1 endif if (reset_to_atmo .gt. 0 .or. (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)) + SET_ATMO_MIN(dens(i,j,k), sdet*GRHydro_rho_min, r(i,j,k)) SET_ATMO_MIN(rho(i,j,k), GRHydro_rho_min, r(i,j,k)) scon(i,j,k,:) = 0.d0 vup(i,j,k,:) = 0.d0 @@ -142,7 +142,7 @@ subroutine Conservative2PrimitiveHot(CCTK_ARGUMENTS) press(i,j,k),keyerr,anyerr) keytemp = 0 ! w_lorentz=1, so the expression for tau reduces to: - tau(i,j,k) = sqrt(det) * (rho(i,j,k)+rho(i,j,k)*eps(i,j,k)) - dens(i,j,k) + tau(i,j,k) = sdet * (rho(i,j,k)+rho(i,j,k)*eps(i,j,k)) - dens(i,j,k) cycle @@ -153,7 +153,7 @@ subroutine Conservative2PrimitiveHot(CCTK_ARGUMENTS) scon(i,j,k,2),scon(i,j,k,3),tau(i,j,k),Y_e_con(i,j,k),rho(i,j,k),vup(i,j,k,1),& vup(i,j,k,2),vup(i,j,k,3),eps(i,j,k),temperature(i,j,k),y_e(i,j,k),& press(i,j,k),w_lorentz(i,j,k), & - uxx,uxy,uxz,uyy,uyz,uzz,det,x(i,j,k),y(i,j,k), & + uxx,uxy,uxz,uyy,uyz,uzz,sdet,x(i,j,k),y(i,j,k), & z(i,j,k),r(i,j,k),epsnegative,GRHydro_rho_min,pmin, epsmin, & GRHydro_reflevel, GRHydro_C2P_failed(i,j,k), GRHydro_perc_ptol) @@ -198,7 +198,7 @@ subroutine Conservative2PrimitiveHot(CCTK_ARGUMENTS) scon(i,j,k,2),scon(i,j,k,3),tau(i,j,k),Y_e_con(i,j,k),rho(i,j,k),& vup(i,j,k,1),vup(i,j,k,2), vup(i,j,k,3),eps(i,j,k),& temperature(i,j,k),y_e(i,j,k),press(i,j,k),w_lorentz(i,j,k), & - uxx,uxy,uxz,uyy,uyz,uzz,det,x(i,j,k),y(i,j,k), & + uxx,uxy,uxz,uyy,uyz,uzz,sdet,x(i,j,k),y(i,j,k), & z(i,j,k),r(i,j,k),epsnegative,GRHydro_rho_min,pmin, epsmin, & GRHydro_reflevel, GRHydro_C2P_failed(i,j,k), local_perc_ptol) if( abs(GRHydro_C2P_failed(i,j,k)-2.0d0) .lt. 1.0d-10) then @@ -267,13 +267,13 @@ subroutine Conservative2PrimitiveHot(CCTK_ARGUMENTS) vlowz = gxz(i,j,k)*vel(i,j,k,1) & + gyz(i,j,k)*vel(i,j,k,2) & + gzz(i,j,k)*vel(i,j,k,3) - scon(i,j,k,1) = sqrt(det) * (rho(i,j,k)*(1.0d0+eps(i,j,k)) & + scon(i,j,k,1) = sdet * (rho(i,j,k)*(1.0d0+eps(i,j,k)) & + press(i,j,k))*w_lorentz(i,j,k)**2 * vlowx - scon(i,j,k,2) = sqrt(det) * (rho(i,j,k)*(1.0d0+eps(i,j,k)) & + scon(i,j,k,2) = sdet * (rho(i,j,k)*(1.0d0+eps(i,j,k)) & +press(i,j,k))*w_lorentz(i,j,k)**2 * vlowy - scon(i,j,k,3) = sqrt(det) * (rho(i,j,k)*(1.0d0+eps(i,j,k)) & + scon(i,j,k,3) = sdet * (rho(i,j,k)*(1.0d0+eps(i,j,k)) & + press(i,j,k))*w_lorentz(i,j,k)**2 * vlowz - tau(i,j,k) = sqrt(det) * ((rho(i,j,k)*(1.0d0+eps(i,j,k)) & + tau(i,j,k) = sdet * ((rho(i,j,k)*(1.0d0+eps(i,j,k)) & + press(i,j,k))*w_lorentz(i,j,k)**2 - press(i,j,k)) & - dens(i,j,k) endif @@ -291,7 +291,7 @@ end subroutine Conservative2PrimitiveHot subroutine Con2Prim_pt_hot3(cctk_iteration, myproc, ii,jj,kk,handle, dens, & sx, sy, sz, tau, ye_con, rho, velx, vely, & velz, epsilon, temp, ye, press, w_lorentz, uxx, uxy, uxz, uyy, & - uyz, uzz, det, x, y, z, r, epsnegative, GRHydro_rho_min, pmin, epsmin, & + uyz, uzz, sdet, x, y, z, r, epsnegative, GRHydro_rho_min, pmin, epsmin, & GRHydro_reflevel, GRHydro_C2P_failed, local_perc_ptol) implicit none @@ -299,7 +299,7 @@ subroutine Con2Prim_pt_hot3(cctk_iteration, myproc, ii,jj,kk,handle, dens, & DECLARE_CCTK_PARAMETERS CCTK_REAL dens, sx, sy, sz, tau, ye_con, rho, velx, vely, velz, epsilon, & - press, uxx, uxy, uxz, uyy, uyz, uzz, det, w_lorentz, x, & + press, uxx, uxy, uxz, uyy, uyz, uzz, invsdet, sdet, w_lorentz, x, & y, z, r, GRHydro_rho_min CCTK_REAL temp, ye CCTK_REAL s2, f, df, vlowx, vlowy, vlowz @@ -345,12 +345,12 @@ subroutine Con2Prim_pt_hot3(cctk_iteration, myproc, ii,jj,kk,handle, dens, & endif !!$ Undensitize the variables - - udens = dens /sqrt(det) - usx = sx /sqrt(det) - usy = sy /sqrt(det) - usz = sz /sqrt(det) - utau = tau /sqrt(det) + invsdet = 1.0d0/sdet + udens = dens * invsdet + usx = sx * invsdet + usy = sy * invsdet + usz = sz * invsdet + utau = tau * invsdet s2 = usx*usx*uxx + usy*usy*uyy + usz*usz*uzz + 2.*usx*usy*uxy + & 2.*usx*usz*uxz + 2.*usy*usz*uyz @@ -627,7 +627,7 @@ subroutine Con2Prim_pt_hot3(cctk_iteration, myproc, ii,jj,kk,handle, dens, & IF_BELOW_ATMO(rho, GRHydro_rho_min, GRHydro_atmo_tolerance, r) then SET_ATMO_MIN(rho, GRHydro_rho_min, r) !GRHydro_rho_min udens = rho - dens = sqrt(det) * rho + dens = sdet * rho temp = GRHydro_hot_atmo_temp ye = GRHydro_hot_atmo_Y_e ye_con = dens * ye 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 diff --git a/src/GRHydro_Con2PrimM_polytype_pt.c b/src/GRHydro_Con2PrimM_polytype_pt.c index 108d700..b9e8a89 100644 --- a/src/GRHydro_Con2PrimM_polytype_pt.c +++ b/src/GRHydro_Con2PrimM_polytype_pt.c @@ -196,7 +196,7 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_Polytype_pt) ( CCTK_REAL *gyy, CCTK_REAL *gyz, CCTK_REAL *gzz, CCTK_REAL *uxx, CCTK_REAL *uxy, CCTK_REAL *uxz, CCTK_REAL *uyy, CCTK_REAL *uyz, CCTK_REAL *uzz, - CCTK_REAL *det, + CCTK_REAL *sdet, CCTK_INT *epsnegative, CCTK_REAL *retval) @@ -209,8 +209,7 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_Polytype_pt) ( CCTK_REAL rho0,u,p,w,gammasq,gamma,W_last,W,vsq; CCTK_REAL g_o_WBsq, QdB_o_W; CCTK_REAL rho_g, x_rho[1]; - CCTK_REAL detg = (*det); - CCTK_REAL sqrt_detg = sqrt(detg); + CCTK_REAL sqrt_detg = *sdet; CCTK_REAL inv_sqrt_detg = 1./sqrt_detg; CCTK_REAL t2, rho_gm1; CCTK_INT i_increase,ntries ; @@ -256,7 +255,7 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_Polytype_pt) ( fprintf(stdout," *uyy = %26.16e \n", *uyy ); fprintf(stdout," *uyz = %26.16e \n", *uyz ); fprintf(stdout," *uzz = %26.16e \n", *uzz ); - fprintf(stdout," *det = %26.16e \n", *det ); + fprintf(stdout," *sdet = %26.16e \n", *sdet ); fprintf(stdout," *epsnegative = %10d \n", *epsnegative ); fprintf(stdout," *retval = %26.16e \n", *retval ); fflush(stdout); diff --git a/src/GRHydro_Con2PrimM_pt.c b/src/GRHydro_Con2PrimM_pt.c index bf947d3..06edca8 100644 --- a/src/GRHydro_Con2PrimM_pt.c +++ b/src/GRHydro_Con2PrimM_pt.c @@ -119,7 +119,7 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_ptold) ( CCTK_REAL *gyy, CCTK_REAL *gyz, CCTK_REAL *gzz, CCTK_REAL *uxx, CCTK_REAL *uxy, CCTK_REAL *uxz, CCTK_REAL *uyy, CCTK_REAL *uyz, CCTK_REAL *uzz, - CCTK_REAL *det, + CCTK_REAL *sdet, CCTK_INT *epsnegative, CCTK_REAL *retval); @@ -201,7 +201,7 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_ptold) ( CCTK_REAL *gyy, CCTK_REAL *gyz, CCTK_REAL *gzz, CCTK_REAL *uxx, CCTK_REAL *uxy, CCTK_REAL *uxz, CCTK_REAL *uyy, CCTK_REAL *uyz, CCTK_REAL *uzz, - CCTK_REAL *det, + CCTK_REAL *sdet, CCTK_INT *epsnegative, CCTK_REAL *retval) @@ -213,8 +213,7 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_ptold) ( CCTK_REAL QdotB; CCTK_REAL rho0,u,p,w,gammasq,gamma,gtmp,W_last,W,vsq; CCTK_REAL g_o_WBsq, QdB_o_W; - CCTK_REAL detg = (*det); - CCTK_REAL sqrt_detg = sqrt(detg); + CCTK_REAL sqrt_detg = *sdet; CCTK_REAL inv_sqrt_detg = 1./sqrt_detg; CCTK_INT i_increase; @@ -260,7 +259,7 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_ptold) ( fprintf(stdout," *uyy = %26.16e \n", *uyy ); fprintf(stdout," *uyz = %26.16e \n", *uyz ); fprintf(stdout," *uzz = %26.16e \n", *uzz ); - fprintf(stdout," *det = %26.16e \n", *det ); + fprintf(stdout," *sdet = %26.16e \n", *sdet ); fprintf(stdout," *epsnegative = %10d \n", *epsnegative ); fprintf(stdout," *retval = %26.16e \n", *retval ); fflush(stdout); diff --git a/src/GRHydro_Con2PrimM_pt_EOSOmni.c b/src/GRHydro_Con2PrimM_pt_EOSOmni.c index 96e2bd8..3700b5b 100644 --- a/src/GRHydro_Con2PrimM_pt_EOSOmni.c +++ b/src/GRHydro_Con2PrimM_pt_EOSOmni.c @@ -141,7 +141,7 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt) ( CCTK_REAL *gyy, CCTK_REAL *gyz, CCTK_REAL *gzz, CCTK_REAL *uxx, CCTK_REAL *uxy, CCTK_REAL *uxz, CCTK_REAL *uyy, CCTK_REAL *uyz, CCTK_REAL *uzz, - CCTK_REAL *det, + CCTK_REAL *sdet, CCTK_INT *epsnegative, CCTK_REAL *retval); @@ -225,7 +225,7 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt) ( CCTK_REAL *gyy, CCTK_REAL *gyz, CCTK_REAL *gzz, CCTK_REAL *uxx, CCTK_REAL *uxy, CCTK_REAL *uxz, CCTK_REAL *uyy, CCTK_REAL *uyz, CCTK_REAL *uzz, - CCTK_REAL *det, + CCTK_REAL *sdet, CCTK_INT *epsnegative, CCTK_REAL *retval) @@ -237,8 +237,7 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt) ( CCTK_REAL QdotB,VdotB; CCTK_REAL rho0,u,p,w,gammasq,gamma,gtmp,W_last,W,vsq; CCTK_REAL g_o_WBsq, QdB_o_W; - CCTK_REAL detg = (*det); - CCTK_REAL sqrt_detg = sqrt(detg); + CCTK_REAL sqrt_detg = *sdet; CCTK_REAL inv_sqrt_detg = 1./sqrt_detg; CCTK_INT i,j, i_increase ; CCTK_INT nf, nfudgemax, failwarnmode, failinfomode; @@ -313,7 +312,7 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt) ( fprintf(stdout," *uyy = %26.16e \n", *uyy ); fprintf(stdout," *uyz = %26.16e \n", *uyz ); fprintf(stdout," *uzz = %26.16e \n", *uzz ); - fprintf(stdout," *det = %26.16e \n", *det ); + fprintf(stdout," *sdet = %26.16e \n", *sdet ); fprintf(stdout," *epsnegative = %10d \n", *epsnegative ); fprintf(stdout," *retval = %26.16e \n", *retval ); fflush(stdout); @@ -539,7 +538,7 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt) ( fprintf(stdout," *uyy = %26.16e \n", *uyy ); fprintf(stdout," *uyz = %26.16e \n", *uyz ); fprintf(stdout," *uzz = %26.16e \n", *uzz ); - fprintf(stdout," *det = %26.16e \n", *det ); + fprintf(stdout," *sdet = %26.16e \n", *sdet ); fprintf(stdout," *epsnegative = %10d \n", *epsnegative ); fprintf(stdout," *retval = %26.16e \n", *retval ); fflush(stdout); @@ -723,7 +722,7 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt) ( fprintf(stdout," *uyyo = %26.16e \n", *uyy ); fprintf(stdout," *uyzo = %26.16e \n", *uyz ); fprintf(stdout," *uzzo = %26.16e \n", *uzz ); - fprintf(stdout," *deto = %26.16e \n", *det ); + fprintf(stdout," *sdeto = %26.16e \n", *sdet ); fprintf(stdout," *epsnegativeo = %10d \n", *epsnegative ); diff --git a/src/GRHydro_Con2PrimM_pt_EOSOmni.c.hopefullyfixed.c b/src/GRHydro_Con2PrimM_pt_EOSOmni.c.hopefullyfixed.c index a08e933..cd64fb7 100644 --- a/src/GRHydro_Con2PrimM_pt_EOSOmni.c.hopefullyfixed.c +++ b/src/GRHydro_Con2PrimM_pt_EOSOmni.c.hopefullyfixed.c @@ -139,7 +139,7 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt) ( CCTK_REAL *gyy, CCTK_REAL *gyz, CCTK_REAL *gzz, CCTK_REAL *uxx, CCTK_REAL *uxy, CCTK_REAL *uxz, CCTK_REAL *uyy, CCTK_REAL *uyz, CCTK_REAL *uzz, - CCTK_REAL *det, + CCTK_REAL *sdet, CCTK_INT *epsnegative, CCTK_REAL *retval); @@ -225,7 +225,7 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt) ( CCTK_REAL *gyy, CCTK_REAL *gyz, CCTK_REAL *gzz, CCTK_REAL *uxx, CCTK_REAL *uxy, CCTK_REAL *uxz, CCTK_REAL *uyy, CCTK_REAL *uyz, CCTK_REAL *uzz, - CCTK_REAL *det, + CCTK_REAL *sdet, CCTK_INT *epsnegative, CCTK_REAL *retval) @@ -237,8 +237,7 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt) ( CCTK_REAL QdotB; CCTK_REAL rho0,u,p,w,gammasq,gamma,gtmp,W_last,W,vsq; CCTK_REAL g_o_WBsq, QdB_o_W; - CCTK_REAL detg = (*det); - CCTK_REAL sqrt_detg = sqrt(detg); + CCTK_REAL sqrt_detg = *sdetg; CCTK_REAL inv_sqrt_detg = 1./sqrt_detg; CCTK_INT i,j, i_increase ; @@ -296,7 +295,7 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt) ( fprintf(stdout," *uyy = %26.16e \n", *uyy ); fprintf(stdout," *uyz = %26.16e \n", *uyz ); fprintf(stdout," *uzz = %26.16e \n", *uzz ); - fprintf(stdout," *det = %26.16e \n", *det ); + fprintf(stdout," *sdet = %26.16e \n", *sdet ); fprintf(stdout," *epsnegative = %10d \n", *epsnegative ); fprintf(stdout," *retval = %26.16e \n", *retval ); fflush(stdout); diff --git a/src/GRHydro_Con2PrimM_pt_EOSOmniold.c b/src/GRHydro_Con2PrimM_pt_EOSOmniold.c index 0eb4845..7c89190 100644 --- a/src/GRHydro_Con2PrimM_pt_EOSOmniold.c +++ b/src/GRHydro_Con2PrimM_pt_EOSOmniold.c @@ -139,7 +139,7 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt2) ( CCTK_REAL *gyy, CCTK_REAL *gyz, CCTK_REAL *gzz, CCTK_REAL *uxx, CCTK_REAL *uxy, CCTK_REAL *uxz, CCTK_REAL *uyy, CCTK_REAL *uyz, CCTK_REAL *uzz, - CCTK_REAL *det, + CCTK_REAL *sdet, CCTK_INT *epsnegative, CCTK_REAL *retval); @@ -225,7 +225,7 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt2) ( CCTK_REAL *gyy, CCTK_REAL *gyz, CCTK_REAL *gzz, CCTK_REAL *uxx, CCTK_REAL *uxy, CCTK_REAL *uxz, CCTK_REAL *uyy, CCTK_REAL *uyz, CCTK_REAL *uzz, - CCTK_REAL *det, + CCTK_REAL *sdet, CCTK_INT *epsnegative, CCTK_REAL *retval) @@ -237,8 +237,7 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt2) ( CCTK_REAL QdotB; CCTK_REAL rho0,u,p,w,gammasq,gamma,gtmp,W_last,W,vsq; CCTK_REAL g_o_WBsq, QdB_o_W; - CCTK_REAL detg = (*det); - CCTK_REAL sqrt_detg = sqrt(detg); + CCTK_REAL sqrt_detg = *sdet; CCTK_REAL inv_sqrt_detg = 1./sqrt_detg; CCTK_INT i,j, i_increase ; @@ -297,7 +296,7 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt2) ( fprintf(stdout," *uyy = %26.16e \n", *uyy ); fprintf(stdout," *uyz = %26.16e \n", *uyz ); fprintf(stdout," *uzz = %26.16e \n", *uzz ); - fprintf(stdout," *det = %26.16e \n", *det ); + fprintf(stdout," *sdet = %26.16e \n", *sdet ); fprintf(stdout," *epsnegative = %10d \n", *epsnegative ); fprintf(stdout," *retval = %26.16e \n", *retval ); fflush(stdout); diff --git a/src/GRHydro_Con2PrimM_ptee.c b/src/GRHydro_Con2PrimM_ptee.c index fcf639a..59ce313 100644 --- a/src/GRHydro_Con2PrimM_ptee.c +++ b/src/GRHydro_Con2PrimM_ptee.c @@ -99,7 +99,7 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_ptee) ( CCTK_REAL *gyy, CCTK_REAL *gyz, CCTK_REAL *gzz, CCTK_REAL *uxx, CCTK_REAL *uxy, CCTK_REAL *uxz, CCTK_REAL *uyy, CCTK_REAL *uyz, CCTK_REAL *uzz, - CCTK_REAL *det, + CCTK_REAL *sdet, CCTK_INT *epsnegative, CCTK_REAL *retval) ; @@ -195,7 +195,7 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_ptee) ( CCTK_REAL *gyy, CCTK_REAL *gyz, CCTK_REAL *gzz, CCTK_REAL *uxx, CCTK_REAL *uxy, CCTK_REAL *uxz, CCTK_REAL *uyy, CCTK_REAL *uyz, CCTK_REAL *uzz, - CCTK_REAL *det, + CCTK_REAL *sdet, CCTK_INT *epsnegative, CCTK_REAL *retval) @@ -207,8 +207,7 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_ptee) ( CCTK_REAL QdotB; CCTK_REAL rho0,u,p,w,gammasq,gamma,W,vsq; CCTK_REAL g_o_WBsq, QdB_o_W; - CCTK_REAL detg = (*det); - CCTK_REAL sqrt_detg = sqrt(detg); + CCTK_REAL sqrt_detg = *sdet; CCTK_REAL inv_sqrt_detg = 1./sqrt_detg; CCTK_REAL rho_gm1; CCTK_REAL utsq; @@ -257,7 +256,7 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_ptee) ( fprintf(stdout," *uyy = %26.16e \n", *uyy ); fprintf(stdout," *uyz = %26.16e \n", *uyz ); fprintf(stdout," *uzz = %26.16e \n", *uzz ); - fprintf(stdout," *det = %26.16e \n", *det ); + fprintf(stdout," *sdet = %26.16e \n", *sdet ); fprintf(stdout," *epsnegative = %10d \n", *epsnegative ); fprintf(stdout," *retval = %26.16e \n", *retval ); fflush(stdout); diff --git a/src/GRHydro_EoSChangeGamma.F90 b/src/GRHydro_EoSChangeGamma.F90 index b81a267..8a56430 100644 --- a/src/GRHydro_EoSChangeGamma.F90 +++ b/src/GRHydro_EoSChangeGamma.F90 @@ -51,7 +51,6 @@ subroutine GRHydro_EoSChangeGamma(CCTK_ARGUMENTS) DECLARE_CCTK_PARAMETERS CCTK_INT :: i, j, k - CCTK_REAL :: det CCTK_REAL :: local_gamma @@ -106,11 +105,9 @@ subroutine GRHydro_EoSChangeGamma(CCTK_ARGUMENTS) do j = 1, cctk_lsh(2) do i = 1, cctk_lsh(1) - 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 prim2conpolytype(GRHydro_polytrope_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, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),& + sdetg(i,j,k), dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),& tau(i,j,k),rho(i,j,k),velx(i,j,k),vely(i,j,k),velz(i,j,k),& eps(i,j,k),press(i,j,k),w_lorentz(i,j,k)) @@ -156,7 +153,6 @@ subroutine GRHydro_EoSChangeK(CCTK_ARGUMENTS) DECLARE_CCTK_PARAMETERS CCTK_INT :: i, j, k - CCTK_REAL :: det CCTK_REAL :: local_gamma, local_k @@ -233,11 +229,9 @@ subroutine GRHydro_EoSChangeK(CCTK_ARGUMENTS) do j = 1, cctk_lsh(2) do i = 1, cctk_lsh(1) - 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 prim2conpolytype(GRHydro_polytrope_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, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),& + sdetg(i,j,k), dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),& tau(i,j,k),rho(i,j,k),velx(i,j,k),vely(i,j,k),velz(i,j,k),& eps(i,j,k),press(i,j,k),w_lorentz(i,j,k)) @@ -287,7 +281,6 @@ subroutine GRHydro_EoSChangeGammaK_Shibata(CCTK_ARGUMENTS) DECLARE_CCTK_PARAMETERS CCTK_INT :: i, j, k - CCTK_REAL :: det CCTK_REAL :: local_Gamma, local_k, eos_k_initial @@ -345,11 +338,9 @@ subroutine GRHydro_EoSChangeGammaK_Shibata(CCTK_ARGUMENTS) do j = 1, cctk_lsh(2) do i = 1, cctk_lsh(1) - 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 prim2con(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, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),& + sdetg(i,j,k), dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),& tau(i,j,k),rho(i,j,k),velx(i,j,k),vely(i,j,k),velz(i,j,k),& eps(i,j,k),press(i,j,k),w_lorentz(i,j,k)) diff --git a/src/GRHydro_FluxSplit.F90 b/src/GRHydro_FluxSplit.F90 index 7ff4abc..a05b88f 100644 --- a/src/GRHydro_FluxSplit.F90 +++ b/src/GRHydro_FluxSplit.F90 @@ -40,7 +40,7 @@ subroutine GRHydro_FSAlpha(CCTK_ARGUMENTS) DECLARE_CCTK_FUNCTIONS integer :: nx, ny, nz, i, j, k, ierr, max_handle - CCTK_REAL :: uxx, uxy, uxz, uyy, uyz, uzz, det, beta + CCTK_REAL :: uxx, uxy, uxz, uyy, uyz, uzz, sdet, beta CCTK_REAL, dimension(5) :: lambda CCTK_REAL :: alpha1_local, alpha2_local, alpha3_local, alpha4_local, & alpha5_local @@ -67,9 +67,8 @@ subroutine GRHydro_FSAlpha(CCTK_ARGUMENTS) do j = 1, ny do i = 1, nx - det = SPATIAL_DETERMINANT(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),\ - gyy(i,j,k),gyz(i,j,k),gzz(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,& gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),& gyz(i,j,k),gzz(i,j,k)) @@ -109,9 +108,8 @@ subroutine GRHydro_FSAlpha(CCTK_ARGUMENTS) do j = 1, ny do i = 1, nx - det = SPATIAL_DETERMINANT(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),\ - gyy(i,j,k),gyz(i,j,k),gzz(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,& gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),& gyz(i,j,k),gzz(i,j,k)) @@ -151,9 +149,8 @@ subroutine GRHydro_FSAlpha(CCTK_ARGUMENTS) do j = 1, ny do i = 1, nx - det = SPATIAL_DETERMINANT(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),\ - gyy(i,j,k),gyz(i,j,k),gzz(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,& gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),& gyz(i,j,k),gzz(i,j,k)) @@ -301,8 +298,7 @@ subroutine GRHydro_SplitFlux(CCTK_ARGUMENTS) dummy = betax(:,j,k) do i = 1, cctk_lsh(1) - det = SPATIAL_DETERMINANT(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),\ - gyy(i,j,k),gyz(i,j,k),gzz(i,j,k)) + det(i) = sdetg(i,j,k)*sdetg(i,j,k) call UpperMetric(uxx,uxy,uxz,uyy,uyz,uzz,det(i),& gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),& gyz(i,j,k),gzz(i,j,k)) @@ -335,8 +331,7 @@ subroutine GRHydro_SplitFlux(CCTK_ARGUMENTS) dummy = betay(i,:,k) do j = 1, cctk_lsh(2) - det = SPATIAL_DETERMINANT(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),\ - gyy(i,j,k),gyz(i,j,k),gzz(i,j,k)) + det(j) = sdetg(i,j,k)**2 call UpperMetric(uxx,uxy,uxz,uyy,uyz,uzz,det(j),& gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),& gyz(i,j,k),gzz(i,j,k)) @@ -369,8 +364,7 @@ subroutine GRHydro_SplitFlux(CCTK_ARGUMENTS) dummy = betaz(i,j,:) do k = 1, cctk_lsh(3) - det = SPATIAL_DETERMINANT(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),\ - gyy(i,j,k),gyz(i,j,k),gzz(i,j,k)) + det(k) = sdetg(i,j,k)*sdetg(i,j,k) call UpperMetric(uxx,uxy,uxz,uyy,uyz,uzz,det(k),& gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),& gyz(i,j,k),gzz(i,j,k)) diff --git a/src/GRHydro_Macros.h b/src/GRHydro_Macros.h index 2f1f532..5de138d 100644 --- a/src/GRHydro_Macros.h +++ b/src/GRHydro_Macros.h @@ -1,5 +1,5 @@ #define SPATIAL_DETERMINANT(gxx_,gxy_,gxz_,gyy_,gyz_,gzz_) \ - (-(gxz_)**2*(gyy_) + 2*(gxy_)*(gxz_)*(gyz_) - (gxx_)*(gyz_)**2 - (gxy_)**2*(gzz_) \ + (-(gxz_)**2*(gyy_) + 2.0d0*(gxy_)*(gxz_)*(gyz_) - (gxx_)*(gyz_)**2 - (gxy_)**2*(gzz_) \ + (gxx_)*(gyy_)*(gzz_)) #define DOTP(gxx_,gxy_,gxz_,gyy_,gyz_,gzz_,x1_,y1_,z1_,x2_,y2_,z2_) \ diff --git a/src/GRHydro_Prim2Con.F90 b/src/GRHydro_Prim2Con.F90 index 680ecfc..20a6c9d 100644 --- a/src/GRHydro_Prim2Con.F90 +++ b/src/GRHydro_Prim2Con.F90 @@ -51,8 +51,8 @@ subroutine primitive2conservative(CCTK_ARGUMENTS) integer :: i, j, k integer :: keytemp - CCTK_REAL :: g11l,g12l,g13l,g22l,g23l,g33l,avg_detl,& - g11r,g12r,g13r,g22r,g23r,g33r,avg_detr + CCTK_REAL :: g11l,g12l,g13l,g22l,g23l,g33l,avg_sdetl,& + g11r,g12r,g13r,g22r,g23r,g33r,avg_sdetr CCTK_REAL :: xtemp(1) ! save memory when MP is not used @@ -79,7 +79,7 @@ subroutine primitive2conservative(CCTK_ARGUMENTS) end if if(evolve_temper.ne.1) then - !$OMP PARALLEL DO PRIVATE(i, j, k, avg_detl, avg_detr,& + !$OMP PARALLEL DO PRIVATE(i, j, k, avg_sdetl, avg_sdetr,& !$OMP g11l,g12l,g13l,g22l,g23l,g33l, & !$OMP g11r,g12r,g13r,g22r,g23r,g33r) do k = GRHydro_stencil,cctk_lsh(3)-GRHydro_stencil+1 @@ -99,12 +99,12 @@ subroutine primitive2conservative(CCTK_ARGUMENTS) g23r = 0.5d0 * (g23(i,j,k) + g23(i+xoffset,j+yoffset,k+zoffset)) g33r = 0.5d0 * (g33(i,j,k) + g33(i+xoffset,j+yoffset,k+zoffset)) - avg_detl = SPATIAL_DETERMINANT(g11l,g12l,g13l,g22l, g23l,g33l) - avg_detr = SPATIAL_DETERMINANT(g11r,g12r,g13r,g22r, g23r,g33r) + avg_sdetl = sqrt(SPATIAL_DETERMINANT(g11l,g12l,g13l,g22l, g23l,g33l)) + avg_sdetr = sqrt(SPATIAL_DETERMINANT(g11r,g12r,g13r,g22r, g23r,g33r)) call prim2con(GRHydro_eos_handle, g11l,g12l,g13l,g22l,& g23l,g33l, & - avg_detl,densminus(i,j,k),sxminus(i,j,k),& + avg_sdetl,densminus(i,j,k),sxminus(i,j,k),& syminus(i,j,k),szminus(i,j,k),tauminus(i,j,k),& rhominus(i,j,k), & velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k),& @@ -112,7 +112,7 @@ subroutine primitive2conservative(CCTK_ARGUMENTS) call prim2con(GRHydro_eos_handle, g11r,g12r,g13r,& g22r,g23r,g33r, & - avg_detr, densplus(i,j,k),sxplus(i,j,k),& + avg_sdetr, densplus(i,j,k),sxplus(i,j,k),& syplus(i,j,k),szplus(i,j ,k),tauplus(i,j,k),& rhoplus(i,j,k),velxplus(i,j,k),velyplus(i,j,k),& velzplus(i,j,k),epsplus(i,j,k),pressplus(i,j,k),& @@ -125,7 +125,7 @@ subroutine primitive2conservative(CCTK_ARGUMENTS) else if(reconstruct_temper.ne.0) then keytemp = 1 - !$OMP PARALLEL DO PRIVATE(i, j, k, avg_detl, avg_detr, xtemp,& + !$OMP PARALLEL DO PRIVATE(i, j, k, avg_sdetl, avg_sdetr, xtemp,& !$OMP g11l,g12l,g13l,g22l,g23l,g33l, & !$OMP g11r,g12r,g13r,g22r,g23r,g33r) do k = GRHydro_stencil,cctk_lsh(3)-GRHydro_stencil+1 @@ -145,15 +145,15 @@ subroutine primitive2conservative(CCTK_ARGUMENTS) g23r = 0.5d0 * (g23(i,j,k) + g23(i+xoffset,j+yoffset,k+zoffset)) g33r = 0.5d0 * (g33(i,j,k) + g33(i+xoffset,j+yoffset,k+zoffset)) - avg_detl = SPATIAL_DETERMINANT(g11l,g12l,g13l,g22l, g23l,g33l) - avg_detr = SPATIAL_DETERMINANT(g11r,g12r,g13r,g22r, g23r,g33r) + avg_sdetl = sqrt(SPATIAL_DETERMINANT(g11l,g12l,g13l,g22l, g23l,g33l)) + avg_sdetr = sqrt(SPATIAL_DETERMINANT(g11r,g12r,g13r,g22r, g23r,g33r)) call prim2con_hot(GRHydro_eos_handle, keytemp, GRHydro_reflevel,& cctk_iteration,i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),& r(i,j,k),& g11l,g12l,g13l,g22l,& g23l,g33l, & - avg_detl,densminus(i,j,k),sxminus(i,j,k),& + avg_sdetl,densminus(i,j,k),sxminus(i,j,k),& syminus(i,j,k),szminus(i,j,k),tauminus(i,j,k),& rhominus(i,j,k), & velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k),& @@ -165,7 +165,7 @@ subroutine primitive2conservative(CCTK_ARGUMENTS) r(i,j,k),& g11r,g12r,g13r,& g22r,g23r,g33r,& - avg_detr, densplus(i,j,k),sxplus(i,j,k),& + avg_sdetr, densplus(i,j,k),sxplus(i,j,k),& syplus(i,j,k),szplus(i,j ,k),tauplus(i,j,k),& rhoplus(i,j,k),velxplus(i,j,k),velyplus(i,j,k),& velzplus(i,j,k),epsplus(i,j,k),pressplus(i,j,k),& @@ -178,7 +178,7 @@ subroutine primitive2conservative(CCTK_ARGUMENTS) !$OMP END PARALLEL DO else keytemp = 0 - !$OMP PARALLEL DO PRIVATE(i, j, k, avg_detl, avg_detr, xtemp,& + !$OMP PARALLEL DO PRIVATE(i, j, k, avg_sdetl, avg_sdetr, xtemp,& !$OMP g11l,g12l,g13l,g22l,g23l,g33l, & !$OMP g11r,g12r,g13r,g22r,g23r,g33r) do k = GRHydro_stencil,cctk_lsh(3)-GRHydro_stencil+1 @@ -198,8 +198,8 @@ subroutine primitive2conservative(CCTK_ARGUMENTS) g23r = 0.5d0 * (g23(i,j,k) + g23(i+xoffset,j+yoffset,k+zoffset)) g33r = 0.5d0 * (g33(i,j,k) + g33(i+xoffset,j+yoffset,k+zoffset)) - avg_detl = SPATIAL_DETERMINANT(g11l,g12l,g13l,g22l, g23l,g33l) - avg_detr = SPATIAL_DETERMINANT(g11r,g12r,g13r,g22r, g23r,g33r) + avg_sdetl = sqrt(SPATIAL_DETERMINANT(g11l,g12l,g13l,g22l, g23l,g33l)) + avg_sdetr = sqrt(SPATIAL_DETERMINANT(g11r,g12r,g13r,g22r, g23r,g33r)) xtemp(1) = 0.5d0*(temperature(i,j,k) + & temperature(i-xoffset,j-yoffset,k-zoffset)) @@ -208,7 +208,7 @@ subroutine primitive2conservative(CCTK_ARGUMENTS) r(i,j,k),& g11l,g12l,g13l,g22l,& g23l,g33l, & - avg_detl,densminus(i,j,k),sxminus(i,j,k),& + avg_sdetl,densminus(i,j,k),sxminus(i,j,k),& syminus(i,j,k),szminus(i,j,k),tauminus(i,j,k),& rhominus(i,j,k), & velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k),& @@ -222,7 +222,7 @@ subroutine primitive2conservative(CCTK_ARGUMENTS) r(i,j,k),& g11r,g12r,g13r,& g22r,g23r,g33r,& - avg_detr, densplus(i,j,k),sxplus(i,j,k),& + avg_sdetr, densplus(i,j,k),sxplus(i,j,k),& syplus(i,j,k),szplus(i,j ,k),tauplus(i,j,k),& rhoplus(i,j,k),velxplus(i,j,k),velyplus(i,j,k),& velzplus(i,j,k),epsplus(i,j,k),pressplus(i,j,k),& @@ -254,14 +254,14 @@ end subroutine primitive2conservative @@*/ -subroutine prim2con(handle, gxx, gxy, gxz, gyy, gyz, gzz, det, ddens, & +subroutine prim2con(handle, gxx, gxy, gxz, gyy, gyz, gzz, sdetg, ddens, & dsx, dsy, dsz, dtau , drho, dvelx, dvely, dvelz, deps, dpress, w) implicit none DECLARE_CCTK_PARAMETERS - CCTK_REAL :: gxx, gxy, gxz, gyy, gyz, gzz, det + CCTK_REAL :: gxx, gxy, gxz, gyy, gyz, gzz, sdetg CCTK_REAL :: ddens, dsx, dsy, dsz, dtau, drho, dvelx, dvely, dvelz,& deps, dpress, w, vlowx, vlowy, vlowz CCTK_INT :: handle @@ -284,17 +284,17 @@ subroutine prim2con(handle, gxx, gxy, gxz, gyy, gyz, gzz, det, ddens, & vlowy = gxy*dvelx + gyy*dvely + gyz*dvelz vlowz = gxz*dvelx + gyz*dvely + gzz*dvelz - ddens = sqrt(det) * drho * w - dsx = sqrt(det) * (drho*(1+deps)+dpress)*w*w * vlowx - dsy = sqrt(det) * (drho*(1+deps)+dpress)*w*w * vlowy - dsz = sqrt(det) * (drho*(1+deps)+dpress)*w*w * vlowz - dtau = sqrt(det) * ((drho*(1+deps)+dpress)*w*w - dpress) - ddens + ddens = sdetg * drho * w + dsx = sdetg * (drho*(1+deps)+dpress)*w*w * vlowx + dsy = sdetg * (drho*(1+deps)+dpress)*w*w * vlowy + dsz = sdetg * (drho*(1+deps)+dpress)*w*w * vlowz + dtau = sdetg * ((drho*(1+deps)+dpress)*w*w - dpress) - ddens end subroutine prim2con subroutine prim2con_hot(handle, keytemp, GRHydro_reflevel, cctk_iteration, ii, jj, kk, & - x, y, z, r, gxx, gxy, gxz, gyy, gyz, gzz, det, ddens, & + x, y, z, r, gxx, gxy, gxz, gyy, gyz, gzz, sdetg, ddens, & dsx, dsy, dsz, dtau , drho, dvelx, dvely, dvelz, deps, dpress, w, & temp,ye) @@ -303,12 +303,12 @@ subroutine prim2con_hot(handle, keytemp, GRHydro_reflevel, cctk_iteration, ii, j DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS - CCTK_REAL :: gxx, gxy, gxz, gyy, gyz, gzz, det + CCTK_REAL :: gxx, gxy, gxz, gyy, gyz, gzz, sdetg CCTK_REAL :: ddens, dsx, dsy, dsz, dtau, drho(1), dvelx, dvely, dvelz,& deps(1), dpress(1), w, vlowx, vlowy, vlowz CCTK_REAL :: temp(1),ye(1), x, y, z, r CCTK_INT :: handle, GRHydro_reflevel, cctk_iteration, ii, jj, kk - CCTK_REAL :: sdet, h + CCTK_REAL :: h character(len=512) warnline ! begin EOS Omni vars @@ -439,13 +439,12 @@ subroutine prim2con_hot(handle, keytemp, GRHydro_reflevel, cctk_iteration, ii, j vlowy = gxy*dvelx + gyy*dvely + gyz*dvelz vlowz = gxz*dvelx + gyz*dvely + gzz*dvelz - sdet = sqrt(det) h = drho(1)*(1.0d0+deps(1))+dpress(1) - ddens = sdet * drho(1) * w - dsx = sdet * h*w*w * vlowx - dsy = sdet * h*w*w * vlowy - dsz = sdet * h*w*w * vlowz - dtau = sdet * (h*w*w - dpress(1)) - ddens + ddens = sdetg * drho(1) * w + dsx = sdetg * h*w*w * vlowx + dsy = sdetg * h*w*w * vlowy + dsz = sdetg * h*w*w * vlowz + dtau = sdetg * (h*w*w - dpress(1)) - ddens end subroutine prim2con_hot @@ -483,8 +482,8 @@ subroutine Primitive2ConservativeCells(CCTK_ARGUMENTS) CCTK_INT :: i, j, k CCTK_INT :: keytemp - CCTK_REAL :: det + character(len=512) :: warnline ! save memory when MP is not used CCTK_INT :: GRHydro_UseGeneralCoordinates CCTK_REAL, DIMENSION(:,:,:), POINTER :: g11, g12, g13, g22, g23, g33 @@ -509,41 +508,36 @@ subroutine Primitive2ConservativeCells(CCTK_ARGUMENTS) end if if(evolve_temper.ne.1) then - !$OMP PARALLEL DO PRIVATE(i, j, k, det) + !$OMP PARALLEL DO PRIVATE(i, j, k) do k = 1,cctk_lsh(3) do j = 1,cctk_lsh(2) do i = 1,cctk_lsh(1) - 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 prim2con(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, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),& + sdetg(i,j,k), dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),& tau(i,j,k),rho(i,j,k),velx(i,j,k),vely(i,j,k),velz(i,j,k),& eps(i,j,k),press(i,j,k),w_lorentz(i,j,k)) + end do end do end do !$OMP END PARALLEL DO else keytemp = 0 - !$OMP PARALLEL DO PRIVATE(i, j, k, det) + !$OMP PARALLEL DO PRIVATE(i, j, k) do k = 1,cctk_lsh(3) do j = 1,cctk_lsh(2) do i = 1,cctk_lsh(1) - 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 prim2con_hot(GRHydro_eos_handle,keytemp,& GRHydro_reflevel,cctk_iteration,& i,j,k,& x(i,j,k),y(i,j,k),z(i,j,k),r(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),& - det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),& + sdetg(i,j,k), dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),& tau(i,j,k),rho(i,j,k),velx(i,j,k),vely(i,j,k),velz(i,j,k),& eps(i,j,k),press(i,j,k),w_lorentz(i,j,k), & temperature(i,j,k),y_e(i,j,k)) @@ -599,8 +593,8 @@ subroutine Prim2ConservativePolytype(CCTK_ARGUMENTS) DECLARE_CCTK_PARAMETERS integer :: i, j, k - CCTK_REAL :: g11l,g12l,g13l,g22l,g23l,g33l,avg_detl,& - g11r,g12r,g13r,g22r,g23r,g33r,avg_detr + CCTK_REAL :: g11l,g12l,g13l,g22l,g23l,g33l,avg_sdetl,& + g11r,g12r,g13r,g22r,g23r,g33r,avg_sdetr ! save memory when MP is not used CCTK_INT :: GRHydro_UseGeneralCoordinates @@ -625,8 +619,8 @@ subroutine Prim2ConservativePolytype(CCTK_ARGUMENTS) vup => vel end if - !$OMP PARALLEL DO PRIVATE(i, j, k, g11l,g12l,g13l,g22l,g23l,g33l,avg_detl,& - !$OMP g11r,g12r,g13r,g22r,g23r,g33r,avg_detr) + !$OMP PARALLEL DO PRIVATE(i, j, k, g11l,g12l,g13l,g22l,g23l,g33l,avg_sdetl,& + !$OMP g11r,g12r,g13r,g22r,g23r,g33r,avg_sdetr) do k = GRHydro_stencil,cctk_lsh(3)-GRHydro_stencil+1 do j = GRHydro_stencil,cctk_lsh(2)-GRHydro_stencil+1 do i = GRHydro_stencil,cctk_lsh(1)-GRHydro_stencil+1 @@ -644,19 +638,19 @@ subroutine Prim2ConservativePolytype(CCTK_ARGUMENTS) g23r = 0.5d0 * (g23(i,j,k) + g23(i+xoffset,j+yoffset,k+zoffset)) g33r = 0.5d0 * (g33(i,j,k) + g33(i+xoffset,j+yoffset,k+zoffset)) - avg_detl = SPATIAL_DETERMINANT(g11l,g12l,g13l,g22l, g23l,g33l) - avg_detr = SPATIAL_DETERMINANT(g11r,g12r,g13r,g22r, g23r,g33r) + avg_sdetl = sqrt(SPATIAL_DETERMINANT(g11l,g12l,g13l,g22l, g23l,g33l)) + avg_sdetr = sqrt(SPATIAL_DETERMINANT(g11r,g12r,g13r,g22r, g23r,g33r)) call prim2conpolytype(GRHydro_eos_handle, g11l,g12l,g13l,& g22l,g23l,g33l, & - avg_detl,densminus(i,j,k),sxminus(i,j,k),& + avg_sdetl,densminus(i,j,k),sxminus(i,j,k),& syminus(i,j,k),szminus(i,j,k),tauminus(i,j,k),rhominus(i,j,k), & velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k),& epsminus(i,j,k),pressminus(i,j,k),w_lorentzminus(i, j, k)) call prim2conpolytype(GRHydro_eos_handle, g11r,g12r,g13r,& g22r,g23r,g33r, & - avg_detr, densplus(i,j,k),sxplus(i,j,k),& + avg_sdetr, densplus(i,j,k),sxplus(i,j,k),& syplus(i,j,k),szplus(i,j ,k),tauplus(i,j,k),& rhoplus(i,j,k),velxplus(i,j,k),velyplus(i,j,k),& velzplus(i,j,k),epsplus(i,j,k),pressplus(i,j,k),& @@ -689,14 +683,14 @@ end subroutine Prim2ConservativePolytype @@*/ subroutine prim2conpolytype(handle, gxx, gxy, gxz, gyy, gyz, & - gzz, det, ddens, & + gzz, sdetg, ddens, & dsx, dsy, dsz, dtau , drho, dvelx, dvely, dvelz, deps, dpress, w) implicit none DECLARE_CCTK_PARAMETERS - CCTK_REAL :: gxx, gxy, gxz, gyy, gyz, gzz, det + CCTK_REAL :: gxx, gxy, gxz, gyy, gyz, gzz, sdetg CCTK_REAL :: ddens, dsx, dsy, dsz, dtau, drho, dvelx, dvely, dvelz,& deps, dpress, w_tmp, w, vlowx, vlowy, vlowz, sqrtdet CCTK_INT :: handle @@ -737,12 +731,11 @@ subroutine prim2conpolytype(handle, gxx, gxy, gxz, gyy, gyz, & vlowy = gxy*dvelx + gyy*dvely + gyz*dvelz vlowz = gxz*dvelx + gyz*dvely + gzz*dvelz - sqrtdet = sqrt(det) - ddens = sqrtdet * drho * w - dsx = sqrtdet * (drho*(1+deps)+dpress)*w*w * vlowx - dsy = sqrtdet * (drho*(1+deps)+dpress)*w*w * vlowy - dsz = sqrtdet * (drho*(1+deps)+dpress)*w*w * vlowz - dtau = sqrtdet * ((drho*(1+deps)+dpress)*w*w - dpress) - ddens + ddens = sdetg * drho * w + dsx = sdetg * (drho*(1+deps)+dpress)*w*w * vlowx + dsy = sdetg * (drho*(1+deps)+dpress)*w*w * vlowy + dsz = sdetg * (drho*(1+deps)+dpress)*w*w * vlowz + dtau = sdetg * ((drho*(1+deps)+dpress)*w*w - dpress) - ddens end subroutine prim2conpolytype @@ -778,7 +771,6 @@ subroutine Primitive2ConservativePolyCells(CCTK_ARGUMENTS) DECLARE_CCTK_PARAMETERS CCTK_INT :: i, j, k - CCTK_REAL :: det ! save memory when MP is not used CCTK_INT :: GRHydro_UseGeneralCoordinates @@ -803,17 +795,14 @@ subroutine Primitive2ConservativePolyCells(CCTK_ARGUMENTS) vup => vel end if - !$OMP PARALLEL DO PRIVATE(i, j, k, det) + !$OMP PARALLEL DO PRIVATE(i, j, k) do k = 1,cctk_lsh(3) do j = 1,cctk_lsh(2) do i = 1,cctk_lsh(1) - 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 prim2conpolytype(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, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),& + sdetg(i,j,k), dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),& tau(i,j,k),rho(i,j,k),velx(i,j,k),vely(i,j,k),velz(i,j,k),& eps(i,j,k),press(i,j,k),w_lorentz(i,j,k)) @@ -864,8 +853,8 @@ subroutine Prim2ConservativeTracer(CCTK_ARGUMENTS) DECLARE_CCTK_PARAMETERS integer :: i, j, k - CCTK_REAL :: g11l,g12l,g13l,g22l,g23l,g33l,avg_detl,& - g11r,g12r,g13r,g22r,g23r,g33r,avg_detr + CCTK_REAL :: g11l,g12l,g13l,g22l,g23l,g33l,avg_sdetl,& + g11r,g12r,g13r,g22r,g23r,g33r,avg_sdetr ! save memory when MP is not used CCTK_INT :: GRHydro_UseGeneralCoordinates @@ -904,10 +893,10 @@ subroutine Prim2ConservativeTracer(CCTK_ARGUMENTS) g23r = 0.5d0 * (g23(i,j,k) + g23(i+xoffset,j+yoffset,k+zoffset)) g33r = 0.5d0 * (g33(i,j,k) + g33(i+xoffset,j+yoffset,k+zoffset)) - avg_detl = SPATIAL_DETERMINANT(g11l,g12l,g13l,g22l, g23l,g33l) - avg_detr = SPATIAL_DETERMINANT(g11r,g12r,g13r,g22r, g23r,g33r) + avg_sdetl = sqrt(SPATIAL_DETERMINANT(g11l,g12l,g13l,g22l, g23l,g33l)) + avg_sdetr = sqrt(SPATIAL_DETERMINANT(g11r,g12r,g13r,g22r, g23r,g33r)) cons_tracerplus(i,j,k,:) = tracerplus(i,j,k,:) * & - sqrt(avg_detr) * rhoplus(i,j,k) / & + avg_sdetr * rhoplus(i,j,k) / & sqrt(1.d0 - & (g11r * velxplus(i,j,k)**2 + & g22r * velyplus(i,j,k)**2 + & @@ -916,7 +905,7 @@ subroutine Prim2ConservativeTracer(CCTK_ARGUMENTS) g13r * velxplus(i,j,k) * velzplus(i,j,k) + & g23r * velyplus(i,j,k) * velzplus(i,j,k) ) ) ) cons_tracerminus(i,j,k,:) = tracerminus(i,j,k,:) * & - sqrt(avg_detl) * rhominus(i,j,k) / & + avg_sdetl * rhominus(i,j,k) / & sqrt(1.d0 - & (g11l * velxminus(i,j,k)**2 + & g22l * velyminus(i,j,k)**2 + & diff --git a/src/GRHydro_Prim2ConAM.F90 b/src/GRHydro_Prim2ConAM.F90 index 7bda654..2aee515 100644 --- a/src/GRHydro_Prim2ConAM.F90 +++ b/src/GRHydro_Prim2ConAM.F90 @@ -54,8 +54,8 @@ subroutine primitive2conservativeAM(CCTK_ARGUMENTS) DECLARE_CCTK_PARAMETERS integer :: i, j, k - 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 :: g11l,g12l,g13l,g22l,g23l,g33l,& ! g11r,g12r,g13r,g22r,g23r,g33r CCTK_REAL :: xtemp(1) @@ -76,7 +76,7 @@ subroutine primitive2conservativeAM(CCTK_ARGUMENTS) if(evolve_temper.ne.1) then - !$OMP PARALLEL DO PRIVATE(k,j,i,avg_detl,avg_detr,& + !$OMP PARALLEL DO PRIVATE(k,j,i,avg_sdetl,avg_sdetr,& !$OMP gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,& !$OMP gxxr,gxyr,gxzr,gyyr,gyzr,gzzr) do k = GRHydro_stencil, cctk_lsh(3)-GRHydro_stencil+1 + transport_constraints*(1-zoffset) @@ -96,11 +96,11 @@ subroutine primitive2conservativeAM(CCTK_ARGUMENTS) gyzr = 0.5d0 * (gyz(i,j,k) + gyz(i+xoffset,j+yoffset,k+zoffset)) gzzr = 0.5d0 * (gzz(i,j,k) + gzz(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) + avg_sdetl = sqrt(SPATIAL_DETERMINANT(gxxl,gxyl,gxzl,gyyl, gyzl,gzzl)) + avg_sdetr = sqrt(SPATIAL_DETERMINANT(gxxr,gxyr,gxzr,gyyr, gyzr,gzzr)) call prim2conAM(GRHydro_eos_handle, gxxl,gxyl,gxzl,gyyl,gyzl,gzzl, & - avg_detl,densminus(i,j,k),sxminus(i,j,k),& + avg_sdetl,densminus(i,j,k),sxminus(i,j,k),& syminus(i,j,k),szminus(i,j,k),tauminus(i,j,k),& rhominus(i,j,k), & velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k),& @@ -108,7 +108,7 @@ subroutine primitive2conservativeAM(CCTK_ARGUMENTS) Bvecyminus(i,j,k), Bveczminus(i,j,k), w_lorentzminus(i, j, k)) call prim2conAM(GRHydro_eos_handle, gxxr,gxyr,gxzr,gyyr,gyzr,gzzr, & - avg_detr, densplus(i,j,k),sxplus(i,j,k),& + avg_sdetr, densplus(i,j,k),sxplus(i,j,k),& syplus(i,j,k),szplus(i,j ,k),tauplus(i,j,k),& rhoplus(i,j,k),velxplus(i,j,k),velyplus(i,j,k),& velzplus(i,j,k),epsplus(i,j,k),pressplus(i,j,k),& @@ -122,7 +122,7 @@ subroutine primitive2conservativeAM(CCTK_ARGUMENTS) else - !$OMP PARALLEL DO PRIVATE(i, j, k, avg_detl, avg_detr, xtemp,& + !$OMP PARALLEL DO PRIVATE(i, j, k, avg_sdetl, avg_sdetr, xtemp,& !$OMP gxxl,gxyl,gxzl,gyyl,gyzl,gzzl, & !$OMP gxxr,gxyr,gxzr,gyyr,gyzr,gzzr) do k = GRHydro_stencil, cctk_lsh(3)-GRHydro_stencil+1 + transport_constraints*(1-zoffset) @@ -142,8 +142,8 @@ subroutine primitive2conservativeAM(CCTK_ARGUMENTS) gyzr = 0.5d0 * (gyz(i,j,k) + gyz(i+xoffset,j+yoffset,k+zoffset)) gzzr = 0.5d0 * (gzz(i,j,k) + gzz(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) + avg_sdetl = sqrt(SPATIAL_DETERMINANT(gxxl,gxyl,gxzl,gyyl,gyzl,gzzl)) + avg_sdetr = sqrt(SPATIAL_DETERMINANT(gxxr,gxyr,gxzr,gyyr,gyzr,gzzr)) ! we do not have plus/minus vars for temperature since ! eps is reconstructed. Hence, we do not update the @@ -161,7 +161,7 @@ subroutine primitive2conservativeAM(CCTK_ARGUMENTS) !$OMP END CRITICAL call prim2conAM_hot(GRHydro_eos_handle, GRHydro_reflevel,& i,j,k,x(i,j,k),y(i,j,k),z(i,j,k), gxxl,gxyl,gxzl,gyyl,gyzl,gzzl, & - avg_detl,densminus(i,j,k),sxminus(i,j,k),& + avg_sdetl,densminus(i,j,k),sxminus(i,j,k),& syminus(i,j,k),szminus(i,j,k),tauminus(i,j,k),& rhominus(i,j,k), & velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k),& @@ -177,7 +177,7 @@ subroutine primitive2conservativeAM(CCTK_ARGUMENTS) call prim2conAM_hot(GRHydro_eos_handle, GRHydro_reflevel,& i,j,k,x(i,j,k),y(i,j,k),z(i,j,k), gxxr,gxyr,gxzr,gyyr,gyzr,gzzr, & - avg_detr, densplus(i,j,k),sxplus(i,j,k),& + avg_sdetr, densplus(i,j,k),sxplus(i,j,k),& syplus(i,j,k),szplus(i,j,k),tauplus(i,j,k),& rhoplus(i,j,k),velxplus(i,j,k),velyplus(i,j,k),& velzplus(i,j,k),epsplus(i,j,k),pressplus(i,j,k),& @@ -208,15 +208,17 @@ end subroutine primitive2conservativeAM @@*/ -subroutine prim2conAM_hot(handle, GRHydro_reflevel, ii, jj, kk, x, y, z, gxx, gxy, gxz, gyy, gyz, gzz, det, ddens, & - dsx, dsy, dsz, dtau, drho, dvelx, dvely, dvelz, deps, dpress, dBvcx, dBvcy, dBvcz, w, temp, ye) +subroutine prim2conAM_hot(handle, GRHydro_reflevel, ii, jj, kk, x, y, z, & + gxx, gxy, gxz, gyy, gyz, gzz, sdet, ddens, & + dsx, dsy, dsz, dtau, drho, dvelx, dvely, dvelz, & + deps, dpress, dBvcx, dBvcy, dBvcz, w, temp, ye) implicit none DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS - CCTK_REAL :: gxx, gxy, gxz, gyy, gyz, gzz, det + CCTK_REAL :: gxx, gxy, gxz, gyy, gyz, gzz, sdet CCTK_REAL, dimension(1) :: ddens, dsx, dsy, dsz, dtau, & drho, dvelx, dvely, dvelz,& deps, dpress, dBvcx, dBvcy, dBvcz, vlowx, vlowy, vlowz @@ -333,23 +335,23 @@ subroutine prim2conAM_hot(handle, GRHydro_reflevel, ii, jj, kk, x, y, z, gxx, gx blowz = (gxz*dBvcx + gyz*dBvcy + gzz*dBvcz)/w + & w*Bdotv*vlowz - ddens = sqrt(det) * drho * w - dsx = sqrt(det) * ((drho*(1+deps)+dpress+b2)*w*w * vlowx - & + ddens = sdet * drho * w + dsx = sdet * ((drho*(1+deps)+dpress+b2)*w*w * vlowx - & ab0*blowx) - dsy = sqrt(det) * ((drho*(1+deps)+dpress+b2)*w*w * vlowy - & + dsy = sdet * ((drho*(1+deps)+dpress+b2)*w*w * vlowy - & ab0*blowy) - dsz = sqrt(det) * ((drho*(1+deps)+dpress+b2)*w*w * vlowz - & + dsz = sdet * ((drho*(1+deps)+dpress+b2)*w*w * vlowz - & ab0*blowz) - dtau = sqrt(det) * ((drho*(1+deps)+dpress+b2)*w*w - dpress-b2/2.0-ab0**2) - ddens + dtau = sdet * ((drho*(1+deps)+dpress+b2)*w*w - dpress-b2/2.0-ab0**2) - ddens !!$OMP CRITICAL - !write(NaN_WarnLine,'(a100,3g15.6)') '(dens out, sqrt(det), w:)', ddens,sqrt(det),w + !write(NaN_WarnLine,'(a100,3g15.6)') '(dens out, sdet, w:)', ddens,sdet,w !call CCTK_WARN(1, NaN_WarnLine) !!$OMP END CRITICAL end subroutine prim2conAM_hot -subroutine prim2conAM(handle, gxx, gxy, gxz, gyy, gyz, gzz, det, ddens, & +subroutine prim2conAM(handle, gxx, gxy, gxz, gyy, gyz, gzz, sdet, ddens, & dsx, dsy, dsz, dtau, drho, dvelx, dvely, dvelz, deps, dpress, dBvcx, dBvcy, dBvcz, w) implicit none @@ -357,7 +359,7 @@ subroutine prim2conAM(handle, gxx, gxy, gxz, gyy, gyz, gzz, det, ddens, & DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS - CCTK_REAL :: gxx, gxy, gxz, gyy, gyz, gzz, det + CCTK_REAL :: gxx, gxy, gxz, gyy, gyz, gzz, sdet CCTK_REAL, dimension(1) :: ddens, dsx, dsy, dsz, dtau, & drho, dvelx, dvely, dvelz,& deps, dpress, dBvcx, dBvcy, dBvcz, vlowx, vlowy, vlowz @@ -403,14 +405,14 @@ subroutine prim2conAM(handle, gxx, gxy, gxz, gyy, gyz, gzz, det, ddens, & blowz = (gxz*dBvcx + gyz*dBvcy + gzz*dBvcz)/w + & w*Bdotv*vlowz - ddens = sqrt(det) * drho * w - dsx = sqrt(det) * ((drho*(1+deps)+dpress+b2)*w*w * vlowx - & + ddens = sdet * drho * w + dsx = sdet * ((drho*(1+deps)+dpress+b2)*w*w * vlowx - & ab0*blowx) - dsy = sqrt(det) * ((drho*(1+deps)+dpress+b2)*w*w * vlowy - & + dsy = sdet * ((drho*(1+deps)+dpress+b2)*w*w * vlowy - & ab0*blowy) - dsz = sqrt(det) * ((drho*(1+deps)+dpress+b2)*w*w * vlowz - & + dsz = sdet * ((drho*(1+deps)+dpress+b2)*w*w * vlowz - & ab0*blowz) - dtau = sqrt(det) * ((drho*(1+deps)+dpress+b2)*w*w - dpress-b2/2.0-ab0**2) - ddens + dtau = sdet * ((drho*(1+deps)+dpress+b2)*w*w - dpress-b2/2.0-ab0**2) - ddens end subroutine prim2conAM @@ -439,25 +441,24 @@ subroutine Primitive2ConservativeCellsAM(CCTK_ARGUMENTS) DECLARE_CCTK_PARAMETERS CCTK_INT :: i, j, k - CCTK_REAL :: det + CCTK_REAL :: sdet CCTK_REAL :: maxtau0 maxtau0 = -1.0d60 if(evolve_temper.ne.1) then - !$OMP PARALLEL DO PRIVATE(k,j,i,det), REDUCTION(MAX:maxtau0) + !$OMP PARALLEL DO PRIVATE(k,j,i,sdet), REDUCTION(MAX:maxtau0) do k = GRHydro_stencil,cctk_lsh(3)-GRHydro_stencil+1 do j = GRHydro_stencil,cctk_lsh(2)-GRHydro_stencil+1 do i = GRHydro_stencil,cctk_lsh(1)-GRHydro_stencil+1 - det = SPATIAL_DETERMINANT(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k), \ - gyy(i,j,k),gyz(i,j,k),gzz(i,j,k)) + sdet = sdetg(i,j,k) call prim2conAM(GRHydro_eos_handle,gxx(i,j,k),& gxy(i,j,k),gxz(i,j,k),& gyy(i,j,k),gyz(i,j,k),gzz(i,j,k),& - det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),& + sdet, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),& tau(i,j,k),& rho(i,j,k),velx(i,j,k),vely(i,j,k),velz(i,j,k),& eps(i,j,k),press(i,j,k),Bvecx(i,j,k), & @@ -483,19 +484,18 @@ subroutine Primitive2ConservativeCellsAM(CCTK_ARGUMENTS) !GRHydro_tau_min = GRHydro_tau_min * maxtau0 else - !$OMP PARALLEL DO PRIVATE(k,j,i,det) + !$OMP PARALLEL DO PRIVATE(k,j,i,sdet) do k = GRHydro_stencil,cctk_lsh(3)-GRHydro_stencil+1 do j = GRHydro_stencil,cctk_lsh(2)-GRHydro_stencil+1 do i = GRHydro_stencil,cctk_lsh(1)-GRHydro_stencil+1 - det = SPATIAL_DETERMINANT(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k), \ - gyy(i,j,k),gyz(i,j,k),gzz(i,j,k)) + sdet = sdetg(i,j,k) call prim2conAM_hot(GRHydro_eos_handle,GRHydro_reflevel,& i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),& gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),& gyy(i,j,k),gyz(i,j,k),gzz(i,j,k),& - det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),& + sdet, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),& tau(i,j,k),& rho(i,j,k),velx(i,j,k),vely(i,j,k),velz(i,j,k),& eps(i,j,k),press(i,j,k),Bvecx(i,j,k), & @@ -548,13 +548,13 @@ subroutine Prim2ConservativePolytypeAM(CCTK_ARGUMENTS) DECLARE_CCTK_PARAMETERS integer :: i, j, k - 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 ! constraint transport needs to be able to average fluxes in the directions ! other that flux_direction, which in turn need the primitives on interfaces - !$OMP PARALLEL DO PRIVATE(i, j, k, gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_detl,& - !$OMP gxxr,gxyr,gxzr,gyyr,gyzr,gzzr,avg_detr) + !$OMP PARALLEL DO PRIVATE(i, j, k, gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_sdetl,& + !$OMP gxxr,gxyr,gxzr,gyyr,gyzr,gzzr,avg_sdetr) do k = GRHydro_stencil, cctk_lsh(3)-GRHydro_stencil+1 + transport_constraints*(1-zoffset) do j = GRHydro_stencil, cctk_lsh(2)-GRHydro_stencil+1 + transport_constraints*(1-yoffset) do i = GRHydro_stencil, cctk_lsh(1)-GRHydro_stencil+1 + transport_constraints*(1-xoffset) @@ -572,12 +572,12 @@ subroutine Prim2ConservativePolytypeAM(CCTK_ARGUMENTS) gyzr = 0.5d0 * (gyz(i,j,k) + gyz(i+xoffset,j+yoffset,k+zoffset)) gzzr = 0.5d0 * (gzz(i,j,k) + gzz(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) + avg_sdetl = sqrt(SPATIAL_DETERMINANT(gxxl,gxyl,gxzl,gyyl, gyzl,gzzl)) + avg_sdetr = sqrt(SPATIAL_DETERMINANT(gxxr,gxyr,gxzr,gyyr, gyzr,gzzr)) call prim2conpolytypeAM(GRHydro_eos_handle, gxxl,gxyl,gxzl,& gyyl,gyzl,gzzl, & - avg_detl,densminus(i,j,k),sxminus(i,j,k),& + avg_sdetl,densminus(i,j,k),sxminus(i,j,k),& syminus(i,j,k),szminus(i,j,k),tauminus(i,j,k),& rhominus(i,j,k),& velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k),& @@ -586,7 +586,7 @@ subroutine Prim2ConservativePolytypeAM(CCTK_ARGUMENTS) call prim2conpolytypeAM(GRHydro_eos_handle, gxxr,gxyr,gxzr,& gyyr,gyzr,gzzr, & - avg_detr,densplus(i,j,k),sxplus(i,j,k),& + avg_sdetr,densplus(i,j,k),sxplus(i,j,k),& syplus(i,j,k),szplus(i,j,k),tauplus(i,j,k),& rhoplus(i,j,k),& velxplus(i,j,k),velyplus(i,j,k),velzplus(i,j,k),& @@ -620,14 +620,14 @@ end subroutine Prim2ConservativePolytypeAM @@*/ subroutine prim2conpolytypeAM(handle, gxx, gxy, gxz, gyy, gyz, & - gzz, det, ddens, dsx, dsy, dsz, dtau, & + gzz, sdet, ddens, dsx, dsy, dsz, dtau, & drho, dvelx, dvely, dvelz, deps, dpress, dBvcx, dBvcy, dBvcz, w) implicit none DECLARE_CCTK_PARAMETERS - CCTK_REAL :: gxx, gxy, gxz, gyy, gyz, gzz, det + CCTK_REAL :: gxx, gxy, gxz, gyy, gyz, gzz, sdet CCTK_REAL :: ddens, dsx, dsy, dsz, dtau, & drho, dvelx, dvely, dvelz, deps, dpress, dBvcx, dBvcy, dBvcz, & w_tmp, w, vlowx, vlowy, vlowz @@ -680,14 +680,14 @@ subroutine prim2conpolytypeAM(handle, gxx, gxy, gxz, gyy, gyz, & blowz = (gxz*dBvcx + gyz*dBvcy + gzz*dBvcz)/w + & w*Bdotv*vlowz - ddens = sqrt(det) * drho * w - dsx = sqrt(det) * ((drho*(1+deps)+dpress+b2)*w*w * vlowx - & + ddens = sdet * drho * w + dsx = sdet * ((drho*(1+deps)+dpress+b2)*w*w * vlowx - & ab0*blowx) - dsy = sqrt(det) * ((drho*(1+deps)+dpress+b2)*w*w * vlowy - & + dsy = sdet * ((drho*(1+deps)+dpress+b2)*w*w * vlowy - & ab0*blowy) - dsz = sqrt(det) * ((drho*(1+deps)+dpress+b2)*w*w * vlowz - & + dsz = sdet * ((drho*(1+deps)+dpress+b2)*w*w * vlowz - & ab0*blowz) - dtau = sqrt(det) * ((drho*(1+deps)+dpress+b2)*w*w - dpress-b2/2.0-ab0**2) - ddens + dtau = sdet * ((drho*(1+deps)+dpress+b2)*w*w - dpress-b2/2.0-ab0**2) - ddens end subroutine prim2conpolytypeAM @@ -717,19 +717,19 @@ subroutine Primitive2ConservativePolyCellsAM(CCTK_ARGUMENTS) DECLARE_CCTK_PARAMETERS CCTK_INT :: i, j, k - CCTK_REAL :: det + CCTK_REAL :: sdet - !$OMP PARALLEL DO PRIVATE(k,j,i,det) + !$OMP PARALLEL DO PRIVATE(k,j,i,sdet) do k = GRHydro_stencil,cctk_lsh(3)-GRHydro_stencil+1 do j = GRHydro_stencil,cctk_lsh(2)-GRHydro_stencil+1 do i = GRHydro_stencil,cctk_lsh(1)-GRHydro_stencil+1 - det = SPATIAL_DETERMINANT(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),\ - gyy(i,j,k),gyz(i,j,k),gzz(i,j,k)) + sdet = sdetg(i,j,k) + call prim2conpolytypeAM(GRHydro_eos_handle,gxx(i,j,k),gxy(i,j,k),& gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k),& - det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),& + sdet, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),& tau(i,j,k),& rho(i,j,k),velx(i,j,k),vely(i,j,k),velz(i,j,k),& eps(i,j,k),press(i,j,k),Bvecx(i,j,k), Bvecy(i,j,k), & diff --git a/src/GRHydro_Prim2ConM.F90 b/src/GRHydro_Prim2ConM.F90 index cb459ee..a473163 100644 --- a/src/GRHydro_Prim2ConM.F90 +++ b/src/GRHydro_Prim2ConM.F90 @@ -61,8 +61,8 @@ subroutine primitive2conservativeM(CCTK_ARGUMENTS) DECLARE_CCTK_PARAMETERS integer :: i, j, k - CCTK_REAL :: g11l,g12l,g13l,g22l,g23l,g33l,avg_detl,& - g11r,g12r,g13r,g22r,g23r,g33r,avg_detr + CCTK_REAL :: g11l,g12l,g13l,g22l,g23l,g33l,avg_sdetl,& + g11r,g12r,g13r,g22r,g23r,g33r,avg_sdetr CCTK_REAL :: xtemp(1) character(len=256) NaN_WarnLine ! save memory when MP is not used @@ -96,7 +96,7 @@ subroutine primitive2conservativeM(CCTK_ARGUMENTS) if(evolve_temper.ne.1) then - !$OMP PARALLEL DO PRIVATE(k,j,i,avg_detl,avg_detr,& + !$OMP PARALLEL DO PRIVATE(k,j,i,avg_sdetl,avg_sdetr,& !$OMP g11l,g12l,g13l,g22l,g23l,g33l,& !$OMP g11r,g12r,g13r,g22r,g23r,g33r) do k = GRHydro_stencil, cctk_lsh(3)-GRHydro_stencil+1 + transport_constraints*(1-zoffset) @@ -116,11 +116,11 @@ subroutine primitive2conservativeM(CCTK_ARGUMENTS) g23r = 0.5d0 * (g23(i,j,k) + g23(i+xoffset,j+yoffset,k+zoffset)) g33r = 0.5d0 * (g33(i,j,k) + g33(i+xoffset,j+yoffset,k+zoffset)) - avg_detl = SPATIAL_DETERMINANT(g11l,g12l,g13l,g22l, g23l,g33l) - avg_detr = SPATIAL_DETERMINANT(g11r,g12r,g13r,g22r, g23r,g33r) + avg_sdetl = sqrt(SPATIAL_DETERMINANT(g11l,g12l,g13l,g22l, g23l,g33l)) + avg_sdetr = sqrt(SPATIAL_DETERMINANT(g11r,g12r,g13r,g22r, g23r,g33r)) call prim2conM(GRHydro_eos_handle, g11l,g12l,g13l,g22l,g23l,g33l, & - avg_detl,densminus(i,j,k),sxminus(i,j,k),& + avg_sdetl,densminus(i,j,k),sxminus(i,j,k),& syminus(i,j,k),szminus(i,j,k),tauminus(i,j,k),& Bconsxminus(i,j,k),Bconsyminus(i,j,k),Bconszminus(i,j,k), rhominus(i,j,k), & velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k),& @@ -128,7 +128,7 @@ subroutine primitive2conservativeM(CCTK_ARGUMENTS) Bvecyminus(i,j,k), Bveczminus(i,j,k), w_lorentzminus(i, j, k)) call prim2conM(GRHydro_eos_handle, g11r,g12r,g13r,g22r,g23r,g33r, & - avg_detr, densplus(i,j,k),sxplus(i,j,k),& + avg_sdetr, densplus(i,j,k),sxplus(i,j,k),& syplus(i,j,k),szplus(i,j ,k),tauplus(i,j,k),& Bconsxplus(i,j,k),Bconsyplus(i,j,k),Bconszplus(i,j,k), & rhoplus(i,j,k),velxplus(i,j,k),velyplus(i,j,k),& @@ -137,8 +137,8 @@ subroutine primitive2conservativeM(CCTK_ARGUMENTS) w_lorentzplus(i,j,k)) if(evolve_entropy.ne.0) then - entropyconsminus(i,j,k) = sqrt(avg_detl)*entropyminus(i,j,k)*w_lorentzminus(i,j,k) - entropyconsplus(i,j,k) = sqrt(avg_detr)*entropyplus(i,j,k)*w_lorentzplus(i,j,k) + entropyconsminus(i,j,k) = avg_sdetl*entropyminus(i,j,k)*w_lorentzminus(i,j,k) + entropyconsplus(i,j,k) = avg_sdetr*entropyplus(i,j,k)*w_lorentzplus(i,j,k) end if end do @@ -148,7 +148,7 @@ subroutine primitive2conservativeM(CCTK_ARGUMENTS) else - !$OMP PARALLEL DO PRIVATE(i, j, k, avg_detl, avg_detr, xtemp,& + !$OMP PARALLEL DO PRIVATE(i, j, k, avg_sdetl, avg_sdetr, xtemp,& !$OMP g11l,g12l,g13l,g22l,g23l,g33l, & !$OMP g11r,g12r,g13r,g22r,g23r,g33r) do k = GRHydro_stencil, cctk_lsh(3)-GRHydro_stencil+1 + transport_constraints*(1-zoffset) @@ -168,8 +168,8 @@ subroutine primitive2conservativeM(CCTK_ARGUMENTS) g23r = 0.5d0 * (g23(i,j,k) + g23(i+xoffset,j+yoffset,k+zoffset)) g33r = 0.5d0 * (g33(i,j,k) + g33(i+xoffset,j+yoffset,k+zoffset)) - avg_detl = SPATIAL_DETERMINANT(g11l,g12l,g13l,g22l,g23l,g33l) - avg_detr = SPATIAL_DETERMINANT(g11r,g12r,g13r,g22r,g23r,g33r) + avg_sdetl = sqrt(SPATIAL_DETERMINANT(g11l,g12l,g13l,g22l,g23l,g33l)) + avg_sdetr = sqrt(SPATIAL_DETERMINANT(g11r,g12r,g13r,g22r,g23r,g33r)) ! we do not have plus/minus vars for temperature since ! eps is reconstructed. Hence, we do not update the @@ -187,7 +187,7 @@ subroutine primitive2conservativeM(CCTK_ARGUMENTS) endif call prim2conM_hot(GRHydro_eos_handle, GRHydro_reflevel,& i,j,k,x(i,j,k),y(i,j,k),z(i,j,k), g11l,g12l,g13l,g22l,g23l,g33l, & - avg_detl,densminus(i,j,k),sxminus(i,j,k),& + avg_sdetl,densminus(i,j,k),sxminus(i,j,k),& syminus(i,j,k),szminus(i,j,k),tauminus(i,j,k),& Bconsxminus(i,j,k),Bconsyminus(i,j,k),Bconszminus(i,j,k), rhominus(i,j,k), & velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k),& @@ -203,7 +203,7 @@ subroutine primitive2conservativeM(CCTK_ARGUMENTS) call prim2conM_hot(GRHydro_eos_handle, GRHydro_reflevel,& i,j,k,x(i,j,k),y(i,j,k),z(i,j,k), g11r,g12r,g13r,g22r,g23r,g33r, & - avg_detr, densplus(i,j,k),sxplus(i,j,k),& + avg_sdetr, densplus(i,j,k),sxplus(i,j,k),& syplus(i,j,k),szplus(i,j,k),tauplus(i,j,k),& Bconsxplus(i,j,k),Bconsyplus(i,j,k),Bconszplus(i,j,k),& rhoplus(i,j,k),velxplus(i,j,k),velyplus(i,j,k),& @@ -235,15 +235,17 @@ end subroutine primitive2conservativeM @@*/ -subroutine prim2conM_hot(handle, GRHydro_reflevel, ii, jj, kk, x, y, z, gxx, gxy, gxz, gyy, gyz, gzz, det, ddens, & - dsx, dsy, dsz, dtau , dBconsx, dBconsy, dBconsz, drho, dvelx, dvely, dvelz, deps, dpress, dBvcx, dBvcy, dBvcz, w, temp, ye) +subroutine prim2conM_hot(handle, GRHydro_reflevel, ii, jj, kk, & + x, y, z, gxx, gxy, gxz, gyy, gyz, gzz, sdet, ddens, & + dsx, dsy, dsz, dtau , dBconsx, dBconsy, dBconsz, drho, & + dvelx, dvely, dvelz, deps, dpress, dBvcx, dBvcy, dBvcz, w, temp, ye) implicit none DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS - CCTK_REAL :: gxx, gxy, gxz, gyy, gyz, gzz, det + CCTK_REAL :: gxx, gxy, gxz, gyy, gyz, gzz, sdet CCTK_REAL, dimension(1) :: ddens, dsx, dsy, dsz, dtau, & dBconsx, dBconsy, dBconsz, & drho, dvelx, dvely, dvelz,& @@ -362,18 +364,18 @@ subroutine prim2conM_hot(handle, GRHydro_reflevel, ii, jj, kk, x, y, z, gxx, gxy blowz = (gxz*dBvcx + gyz*dBvcy + gzz*dBvcz)/w + & w*Bdotv*vlowz - ddens = sqrt(det) * drho * w - dsx = sqrt(det) * ((drho*(1+deps)+dpress+b2)*w*w * vlowx - & + ddens = sdet * drho * w + dsx = sdet * ((drho*(1+deps)+dpress+b2)*w*w * vlowx - & ab0*blowx) - dsy = sqrt(det) * ((drho*(1+deps)+dpress+b2)*w*w * vlowy - & + dsy = sdet * ((drho*(1+deps)+dpress+b2)*w*w * vlowy - & ab0*blowy) - dsz = sqrt(det) * ((drho*(1+deps)+dpress+b2)*w*w * vlowz - & + dsz = sdet * ((drho*(1+deps)+dpress+b2)*w*w * vlowz - & ab0*blowz) - dtau = sqrt(det) * ((drho*(1+deps)+dpress+b2)*w*w - dpress-b2/2.0-ab0**2) - ddens + dtau = sdet * ((drho*(1+deps)+dpress+b2)*w*w - dpress-b2/2.0-ab0**2) - ddens - dBconsx = sqrt(det)*dBvcx - dBconsy = sqrt(det)*dBvcy - dBconsz = sqrt(det)*dBvcz + dBconsx = sdet*dBvcx + dBconsy = sdet*dBvcy + dBconsz = sdet*dBvcz !!$OMP CRITICAL !write(NaN_WarnLine,'(a100,3g15.6)') '(dens out, sqrt(det), w:)', ddens,sqrt(det),w @@ -382,15 +384,16 @@ subroutine prim2conM_hot(handle, GRHydro_reflevel, ii, jj, kk, x, y, z, gxx, gxy end subroutine prim2conM_hot -subroutine prim2conM(handle, gxx, gxy, gxz, gyy, gyz, gzz, det, ddens, & - dsx, dsy, dsz, dtau , dBconsx, dBconsy, dBconsz, drho, dvelx, dvely, dvelz, deps, dpress, dBvcx, dBvcy, dBvcz, w) +subroutine prim2conM(handle, gxx, gxy, gxz, gyy, gyz, gzz, sdet, ddens, & + dsx, dsy, dsz, dtau , dBconsx, dBconsy, dBconsz, drho, & + dvelx, dvely, dvelz, deps, dpress, dBvcx, dBvcy, dBvcz, w) implicit none DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS - CCTK_REAL :: gxx, gxy, gxz, gyy, gyz, gzz, det + CCTK_REAL :: gxx, gxy, gxz, gyy, gyz, gzz, sdet CCTK_REAL, dimension(1) :: ddens, dsx, dsy, dsz, dtau, & dBconsx, dBconsy, dBconsz, & drho, dvelx, dvely, dvelz,& @@ -437,18 +440,18 @@ subroutine prim2conM(handle, gxx, gxy, gxz, gyy, gyz, gzz, det, ddens, & blowz = (gxz*dBvcx + gyz*dBvcy + gzz*dBvcz)/w + & w*Bdotv*vlowz - ddens = sqrt(det) * drho * w - dsx = sqrt(det) * ((drho*(1+deps)+dpress+b2)*w*w * vlowx - & + ddens = sdet * drho * w + dsx = sdet * ((drho*(1+deps)+dpress+b2)*w*w * vlowx - & ab0*blowx) - dsy = sqrt(det) * ((drho*(1+deps)+dpress+b2)*w*w * vlowy - & + dsy = sdet * ((drho*(1+deps)+dpress+b2)*w*w * vlowy - & ab0*blowy) - dsz = sqrt(det) * ((drho*(1+deps)+dpress+b2)*w*w * vlowz - & + dsz = sdet * ((drho*(1+deps)+dpress+b2)*w*w * vlowz - & ab0*blowz) - dtau = sqrt(det) * ((drho*(1+deps)+dpress+b2)*w*w - dpress-b2/2.0-ab0**2) - ddens + dtau = sdet * ((drho*(1+deps)+dpress+b2)*w*w - dpress-b2/2.0-ab0**2) - ddens - dBconsx = sqrt(det)*dBvcx - dBconsy = sqrt(det)*dBvcy - dBconsz = sqrt(det)*dBvcz + dBconsx = sdet*dBvcx + dBconsy = sdet*dBvcy + dBconsz = sdet*dBvcz end subroutine prim2conM @@ -484,7 +487,7 @@ subroutine Primitive2ConservativeCellsM(CCTK_ARGUMENTS) DECLARE_CCTK_PARAMETERS CCTK_INT :: i, j, k - CCTK_REAL :: det + CCTK_REAL :: sdet CCTK_REAL :: maxtau0 ! save memory when MP is not used CCTK_INT :: GRHydro_UseGeneralCoordinates @@ -517,25 +520,24 @@ subroutine Primitive2ConservativeCellsM(CCTK_ARGUMENTS) if(evolve_temper.ne.1) then - !$OMP PARALLEL DO PRIVATE(k,j,i,det), REDUCTION(MAX:maxtau0) + !$OMP PARALLEL DO PRIVATE(k,j,i,sdet), REDUCTION(MAX:maxtau0) do k = 1,cctk_lsh(3) do j = 1,cctk_lsh(2) do i = 1,cctk_lsh(1) - 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 = sdetg(i,j,k) 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, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),& + sdet, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),& tau(i,j,k),Bconsx(i,j,k),Bconsy(i,j,k),Bconsz(i,j,k),& rho(i,j,k),velx(i,j,k),vely(i,j,k),velz(i,j,k),& eps(i,j,k),press(i,j,k),Bvecx(i,j,k), & Bvecy(i,j,k), Bvecz(i,j,k), w_lorentz(i,j,k)) if(evolve_entropy.ne.0) then - entropycons(i,j,k) = sqrt(det)*entropy(i,j,k)*w_lorentz(i,j,k) + entropycons(i,j,k) = sdet*entropy(i,j,k)*w_lorentz(i,j,k) end if maxtau0 = max(maxtau0,tau(i,j,k)) @@ -558,19 +560,18 @@ subroutine Primitive2ConservativeCellsM(CCTK_ARGUMENTS) !GRHydro_tau_min = GRHydro_tau_min * maxtau0 else - !$OMP PARALLEL DO PRIVATE(k,j,i,det) + !$OMP PARALLEL DO PRIVATE(k,j,i,sdet) do k = 1,cctk_lsh(3) do j = 1,cctk_lsh(2) do i = 1,cctk_lsh(1) - 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 = sdetg(i,j,k) call prim2conM_hot(GRHydro_eos_handle,GRHydro_reflevel,& i,j,k,x(i,j,k),y(i,j,k),z(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),& - det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),& + sdet, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),& tau(i,j,k),Bconsx(i,j,k),Bconsy(i,j,k),Bconsz(i,j,k),& rho(i,j,k),velx(i,j,k),vely(i,j,k),velz(i,j,k),& eps(i,j,k),press(i,j,k),Bvecx(i,j,k), & @@ -630,8 +631,8 @@ subroutine Prim2ConservativePolytypeM(CCTK_ARGUMENTS) DECLARE_CCTK_PARAMETERS integer :: i, j, k - CCTK_REAL :: g11l,g12l,g13l,g22l,g23l,g33l,avg_detl,& - g11r,g12r,g13r,g22r,g23r,g33r,avg_detr + CCTK_REAL :: g11l,g12l,g13l,g22l,g23l,g33l,avg_sdetl,& + g11r,g12r,g13r,g22r,g23r,g33r,avg_sdetr ! save memory when MP is not used CCTK_INT :: GRHydro_UseGeneralCoordinates CCTK_REAL, DIMENSION(:,:,:), POINTER :: g11, g12, g13, g22, g23, g33 @@ -661,8 +662,8 @@ subroutine Prim2ConservativePolytypeM(CCTK_ARGUMENTS) ! constraint transport needs to be able to average fluxes in the directions ! other that flux_direction, which in turn need the primitives on interfaces - !$OMP PARALLEL DO PRIVATE(i, j, k, g11l,g12l,g13l,g22l,g23l,g33l,avg_detl,& - !$OMP g11r,g12r,g13r,g22r,g23r,g33r,avg_detr) + !$OMP PARALLEL DO PRIVATE(i, j, k, g11l,g12l,g13l,g22l,g23l,g33l,avg_sdetl,& + !$OMP g11r,g12r,g13r,g22r,g23r,g33r,avg_sdetr) do k = GRHydro_stencil, cctk_lsh(3)-GRHydro_stencil+1 + transport_constraints*(1-zoffset) do j = GRHydro_stencil, cctk_lsh(2)-GRHydro_stencil+1 + transport_constraints*(1-yoffset) do i = GRHydro_stencil, cctk_lsh(1)-GRHydro_stencil+1 + transport_constraints*(1-xoffset) @@ -680,12 +681,12 @@ subroutine Prim2ConservativePolytypeM(CCTK_ARGUMENTS) g23r = 0.5d0 * (g23(i,j,k) + g23(i+xoffset,j+yoffset,k+zoffset)) g33r = 0.5d0 * (g33(i,j,k) + g33(i+xoffset,j+yoffset,k+zoffset)) - avg_detl = SPATIAL_DETERMINANT(g11l,g12l,g13l,g22l, g23l,g33l) - avg_detr = SPATIAL_DETERMINANT(g11r,g12r,g13r,g22r, g23r,g33r) + avg_sdetl = sqrt(SPATIAL_DETERMINANT(g11l,g12l,g13l,g22l, g23l,g33l)) + avg_sdetr = sqrt(SPATIAL_DETERMINANT(g11r,g12r,g13r,g22r, g23r,g33r)) call prim2conpolytypeM(GRHydro_eos_handle, g11l,g12l,g13l,& g22l,g23l,g33l, & - avg_detl,densminus(i,j,k),sxminus(i,j,k),& + avg_sdetl,densminus(i,j,k),sxminus(i,j,k),& syminus(i,j,k),szminus(i,j,k),tauminus(i,j,k),& Bconsxminus(i,j,k),Bconsyminus(i,j,k),Bconszminus(i,j,k),& rhominus(i,j,k),& @@ -695,7 +696,7 @@ subroutine Prim2ConservativePolytypeM(CCTK_ARGUMENTS) call prim2conpolytypeM(GRHydro_eos_handle, g11r,g12r,g13r,& g22r,g23r,g33r, & - avg_detr,densplus(i,j,k),sxplus(i,j,k),& + avg_sdetr,densplus(i,j,k),sxplus(i,j,k),& syplus(i,j,k),szplus(i,j,k),tauplus(i,j,k),& Bconsxplus(i,j,k),Bconsyplus(i,j,k),Bconszplus(i,j,k),& rhoplus(i,j,k),& @@ -730,14 +731,14 @@ end subroutine Prim2ConservativePolytypeM @@*/ subroutine prim2conpolytypeM(handle, gxx, gxy, gxz, gyy, gyz, & - gzz, det, ddens, dsx, dsy, dsz, dtau, dBconsx, dBconsy, dBconsz, & + gzz, sdet, ddens, dsx, dsy, dsz, dtau, dBconsx, dBconsy, dBconsz, & drho, dvelx, dvely, dvelz, deps, dpress, dBvcx, dBvcy, dBvcz, w) implicit none DECLARE_CCTK_PARAMETERS - CCTK_REAL :: gxx, gxy, gxz, gyy, gyz, gzz, det + CCTK_REAL :: gxx, gxy, gxz, gyy, gyz, gzz, sdet CCTK_REAL :: ddens, dsx, dsy, dsz, dtau, dBconsx, dBconsy, dBconsz, & drho, dvelx, dvely, dvelz, deps, dpress, dBvcx, dBvcy, dBvcz, & w_tmp, w, vlowx, vlowy, vlowz @@ -790,18 +791,18 @@ subroutine prim2conpolytypeM(handle, gxx, gxy, gxz, gyy, gyz, & blowz = (gxz*dBvcx + gyz*dBvcy + gzz*dBvcz)/w + & w*Bdotv*vlowz - ddens = sqrt(det) * drho * w - dsx = sqrt(det) * ((drho*(1+deps)+dpress+b2)*w*w * vlowx - & + ddens = sdet * drho * w + dsx = sdet * ((drho*(1+deps)+dpress+b2)*w*w * vlowx - & ab0*blowx) - dsy = sqrt(det) * ((drho*(1+deps)+dpress+b2)*w*w * vlowy - & + dsy = sdet * ((drho*(1+deps)+dpress+b2)*w*w * vlowy - & ab0*blowy) - dsz = sqrt(det) * ((drho*(1+deps)+dpress+b2)*w*w * vlowz - & + dsz = sdet * ((drho*(1+deps)+dpress+b2)*w*w * vlowz - & ab0*blowz) - dtau = sqrt(det) * ((drho*(1+deps)+dpress+b2)*w*w - dpress-b2/2.0-ab0**2) - ddens + dtau = sdet * ((drho*(1+deps)+dpress+b2)*w*w - dpress-b2/2.0-ab0**2) - ddens - dBconsx = sqrt(det)*dBvcx - dBconsy = sqrt(det)*dBvcy - dBconsz = sqrt(det)*dBvcz + dBconsx = sdet*dBvcx + dBconsy = sdet*dBvcy + dBconsz = sdet*dBvcz end subroutine prim2conpolytypeM @@ -838,7 +839,7 @@ subroutine Primitive2ConservativePolyCellsM(CCTK_ARGUMENTS) DECLARE_CCTK_PARAMETERS CCTK_INT :: i, j, k - CCTK_REAL :: det + CCTK_REAL :: sdet ! save memory when MP is not used CCTK_INT :: GRHydro_UseGeneralCoordinates CCTK_REAL, DIMENSION(:,:,:), POINTER :: g11, g12, g13, g22, g23, g33 @@ -865,17 +866,16 @@ subroutine Primitive2ConservativePolyCellsM(CCTK_ARGUMENTS) Bprim => Bvec end if - !$OMP PARALLEL DO PRIVATE(k,j,i,det) + !$OMP PARALLEL DO PRIVATE(k,j,i,sdet) do k = 1,cctk_lsh(3) do j = 1,cctk_lsh(2) do i = 1,cctk_lsh(1) - 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 = sdetg(i,j,k) call prim2conpolytypeM(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, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),& + sdet, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),& tau(i,j,k),Bconsx(i,j,k),Bconsy(i,j,k),Bconsz(i,j,k),& rho(i,j,k),velx(i,j,k),vely(i,j,k),velz(i,j,k),& eps(i,j,k),press(i,j,k),Bvecx(i,j,k), Bvecy(i,j,k), & diff --git a/src/GRHydro_Source.F90 b/src/GRHydro_Source.F90 index 553c916..1923c2d 100644 --- a/src/GRHydro_Source.F90 +++ b/src/GRHydro_Source.F90 @@ -71,7 +71,7 @@ subroutine SourceTerms(CCTK_ARGUMENTS) CCTK_INT :: i, j, k, na, nb, nc CCTK_REAL :: one, two, half CCTK_REAL :: t00, t0a, t0b, t0c, taa, tab, tac, tbb, tbc, tcc - CCTK_REAL :: sqrtdet, det, uaa, uab, uac, ubb, ubc, ucc, rhoenthalpyW2 + CCTK_REAL :: uaa, uab, uac, ubb, ubc, ucc, rhoenthalpyW2 CCTK_REAL :: shifta, shiftb, shiftc, velashift, velbshift, velcshift CCTK_REAL :: vlowa, vlowb, vlowc CCTK_REAL :: da_betaa, da_betab, da_betac, db_betaa, db_betab,& @@ -191,7 +191,7 @@ subroutine SourceTerms(CCTK_ARGUMENTS) !$OMP PARALLEL DO PRIVATE(i, j, k, local_spatial_order,& !$OMP localgaa,localgab,localgac,localgbb,localgbc,localgcc,& - !$OMP det,sqrtdet,rhoenthalpyW2,shifta,shiftb,shiftc,& + !$OMP rhoenthalpyW2,shifta,shiftb,shiftc,& !$OMP da_betaa,da_betab,da_betac,db_betaa,db_betab,db_betac,& !$OMP dc_betaa,dc_betab,dc_betac,velashift,velbshift,velcshift,& !$OMP vlowa,vlowb,vlowc,t00,t0a,t0b,t0c,taa,tab,tac,tbb,tbc,tcc,& @@ -221,10 +221,8 @@ subroutine SourceTerms(CCTK_ARGUMENTS) localgbc = g23(i,j,k) localgcc = g33(i,j,k) - det = SPATIAL_DETERMINANT(localgaa, localgab, localgac,\ - localgbb, localgbc, localgcc) - sqrtdet = sqrt(det) - call UpperMetric(uaa, uab, uac, ubb, ubc, ucc, det, localgaa,& + call UpperMetric(uaa, uab, uac, ubb, ubc, ucc, & + sdetg(i,j,k)*sdetg(i,j,k), localgaa,& localgab, localgac, localgbb, localgbc, localgcc) !!$ All the matter variables (except velocity) always appear @@ -442,10 +440,10 @@ subroutine SourceTerms(CCTK_ARGUMENTS) invalp densrhs(i,j,k) = 0.d0 - srhs(i,j,k,1) = alp(i,j,k)*sqrtdet*sa_source - srhs(i,j,k,2) = alp(i,j,k)*sqrtdet*sb_source - srhs(i,j,k,3) = alp(i,j,k)*sqrtdet*sc_source - taurhs(i,j,k) = alp(i,j,k)*sqrtdet*tau_source + srhs(i,j,k,1) = alp(i,j,k)*sdetg(i,j,k)*sa_source + srhs(i,j,k,2) = alp(i,j,k)*sdetg(i,j,k)*sb_source + srhs(i,j,k,3) = alp(i,j,k)*sdetg(i,j,k)*sc_source + taurhs(i,j,k) = alp(i,j,k)*sdetg(i,j,k)*tau_source enddo enddo diff --git a/src/GRHydro_SourceAM.F90 b/src/GRHydro_SourceAM.F90 index 7c4280a..d090fa7 100644 --- a/src/GRHydro_SourceAM.F90 +++ b/src/GRHydro_SourceAM.F90 @@ -76,7 +76,7 @@ subroutine SourceTermsAM(CCTK_ARGUMENTS) CCTK_INT :: i, j, k, nx, ny, nz CCTK_REAL :: one, two, half CCTK_REAL :: t00, t0x, t0y, t0z, txx, txy, txz, tyy, tyz, tzz - CCTK_REAL :: sqrtdet, det, uxx, uxy, uxz, uyy, uyz, uzz + CCTK_REAL :: sqrtdet, uxx, uxy, uxz, uyy, uyz, uzz CCTK_REAL :: shiftx, shifty, shiftz, velxshift, velyshift, velzshift CCTK_REAL :: vlowx, vlowy, vlowz CCTK_REAL :: dx_betax, dx_betay, dx_betaz, dy_betax, dy_betay,& @@ -250,7 +250,7 @@ subroutine SourceTermsAM(CCTK_ARGUMENTS) !$OMP PARALLEL DO PRIVATE(i, j, k, local_spatial_order,& !$OMP localgxx,localgxy,localgxz,localgyy,localgyz,localgzz,& - !$OMP det,sqrtdet,shiftx,shifty,shiftz,& + !$OMP sqrtdet,shiftx,shifty,shiftz,& !$OMP dx_betax,dx_betay,dx_betaz,dy_betax,dy_betay,dy_betaz,& !$OMP dz_betax,dz_betay,dz_betaz,velxshift,velyshift,velzshift,& !$OMP vlowx,vlowy,vlowz,Bvecxlow,Bvecylow,Bveczlow, & @@ -286,10 +286,8 @@ subroutine SourceTermsAM(CCTK_ARGUMENTS) localgyz = g23(i,j,k) localgzz = g33(i,j,k) - det = SPATIAL_DETERMINANT(localgxx, localgxy, localgxz,\ - localgyy, localgyz, localgzz) - sqrtdet = sqrt(det) - call UpperMetric(uxx, uxy, uxz, uyy, uyz, uzz, det, localgxx,& + sqrtdet = sdetg(i,j,k) + call UpperMetric(uxx, uxy, uxz, uyy, uyz, uzz, sqrtdet*sqrtdet, localgxx,& localgxy, localgxz, localgyy, localgyz, localgzz) diff --git a/src/GRHydro_SourceM.F90 b/src/GRHydro_SourceM.F90 index cd734eb..4b1321c 100644 --- a/src/GRHydro_SourceM.F90 +++ b/src/GRHydro_SourceM.F90 @@ -76,7 +76,7 @@ subroutine SourceTermsM(CCTK_ARGUMENTS) CCTK_INT :: i, j, k, nx, ny, nz CCTK_REAL :: one, two, half CCTK_REAL :: t00, t0x, t0y, t0z, txx, txy, txz, tyy, tyz, tzz - CCTK_REAL :: sqrtdet, det, uxx, uxy, uxz, uyy, uyz, uzz + CCTK_REAL :: sqrtdet, uxx, uxy, uxz, uyy, uyz, uzz CCTK_REAL :: shiftx, shifty, shiftz, velxshift, velyshift, velzshift CCTK_REAL :: vlowx, vlowy, vlowz CCTK_REAL :: dx_betax, dx_betay, dx_betaz, dy_betax, dy_betay,& @@ -254,7 +254,7 @@ subroutine SourceTermsM(CCTK_ARGUMENTS) !$OMP PARALLEL DO PRIVATE(i, j, k, local_spatial_order,& !$OMP localgxx,localgxy,localgxz,localgyy,localgyz,localgzz,& - !$OMP det,sqrtdet,shiftx,shifty,shiftz,& + !$OMP sqrtdet,shiftx,shifty,shiftz,& !$OMP dx_betax,dx_betay,dx_betaz,dy_betax,dy_betay,dy_betaz,& !$OMP dz_betax,dz_betay,dz_betaz,velxshift,velyshift,velzshift,& !$OMP vlowx,vlowy,vlowz,Bvecxlow,Bvecylow,Bveczlow, & @@ -290,10 +290,8 @@ subroutine SourceTermsM(CCTK_ARGUMENTS) localgyz = g23(i,j,k) localgzz = g33(i,j,k) - det = SPATIAL_DETERMINANT(localgxx, localgxy, localgxz,\ - localgyy, localgyz, localgzz) - sqrtdet = sqrt(det) - call UpperMetric(uxx, uxy, uxz, uyy, uyz, uzz, det, localgxx,& + sqrtdet = sdetg(i,j,k) + call UpperMetric(uxx, uxy, uxz, uyy, uyz, uzz, sqrtdet*sqrtdet, localgxx,& localgxy, localgxz, localgyy, localgyz, localgzz) diff --git a/src/GRHydro_UpdateMask.F90 b/src/GRHydro_UpdateMask.F90 index 2fd47df..60243de 100644 --- a/src/GRHydro_UpdateMask.F90 +++ b/src/GRHydro_UpdateMask.F90 @@ -250,7 +250,7 @@ subroutine GRHydro_AtmosphereReset(CCTK_ARGUMENTS) DECLARE_CCTK_PARAMETERS CCTK_INT :: i, j, k - CCTK_REAL :: det, dummy1, dummy2 + CCTK_REAL :: sdet, dummy1, dummy2 ! save memory when MP is not used @@ -286,7 +286,7 @@ subroutine GRHydro_AtmosphereReset(CCTK_ARGUMENTS) if (verbose.eq.1) call CCTK_INFO("Entering AtmosphereReset.") -!$OMP PARALLEL DO PRIVATE(det,keytemp,i,j,k,anyerr,keyerr, dummy1, dummy2) +!$OMP PARALLEL DO PRIVATE(sdet,keytemp,i,j,k,anyerr,keyerr, dummy1, dummy2) do k = 1, cctk_lsh(3) do j = 1, cctk_lsh(2) do i = 1, cctk_lsh(1) @@ -297,8 +297,7 @@ subroutine GRHydro_AtmosphereReset(CCTK_ARGUMENTS) velx(i,j,k) = 0.0d0 vely(i,j,k) = 0.0d0 velz(i,j,k) = 0.0d0 - 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 = sdetg(i,j,k) if(evolve_temper.ne.0) then ! ! set the temperature to be relatively low @@ -313,7 +312,7 @@ subroutine GRHydro_AtmosphereReset(CCTK_ARGUMENTS) cctk_iteration,i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),r(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), & - det,dens(i,j,k),scon(i,j,k,1), scon(i,j,k,2), scon(i,j,k,3), & + sdet,dens(i,j,k),scon(i,j,k,1), scon(i,j,k,2), scon(i,j,k,3), & tau(i,j,k), rho(i,j,k), velx(i,j,k), vely(i,j,k), & velz(i,j,k), eps(i,j,k), press(i,j,k), w_lorentz(i,j,k),& temperature(i,j,k),y_e(i,j,k)) @@ -321,7 +320,7 @@ subroutine GRHydro_AtmosphereReset(CCTK_ARGUMENTS) else call prim2conpolytype(GRHydro_polytrope_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, & + 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), rho(i,j,k), velx(i,j,k), vely(i,j,k), & velz(i,j,k), eps(i,j,k), press(i,j,k), w_lorentz(i,j,k)) @@ -361,7 +360,7 @@ subroutine GRHydro_InitialAtmosphereReset(CCTK_ARGUMENTS) DECLARE_CCTK_PARAMETERS CCTK_INT :: i, j, k - CCTK_REAL :: det + CCTK_REAL :: sdet CCTK_REAL :: rho_min, dummy1, dummy2 CCTK_INT :: eos_handle @@ -450,7 +449,7 @@ subroutine GRHydro_InitialAtmosphereReset(CCTK_ARGUMENTS) rho_min = rho_min * initial_atmosphere_factor endif -!$OMP PARALLEL DO PRIVATE(det,keytemp,i,j,k,anyerr,keyerr, dummy1, dummy2) +!$OMP PARALLEL DO PRIVATE(sdet,keytemp,i,j,k,anyerr,keyerr, dummy1, dummy2) do k = 1, cctk_lsh(3) do j = 1, cctk_lsh(2) do i = 1, cctk_lsh(1) @@ -460,8 +459,7 @@ subroutine GRHydro_InitialAtmosphereReset(CCTK_ARGUMENTS) GRHydro_enable_internal_excision /= 0 .and. & hydro_excision_mask(i,j,k) .ne. 0) then SET_ATMO_MIN(rho(i,j,k), dummy2, r(i,j,k)) - 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 = sdetg(i,j,k) velx(i,j,k) = 0.0d0 vely(i,j,k) = 0.0d0 @@ -480,7 +478,7 @@ subroutine GRHydro_InitialAtmosphereReset(CCTK_ARGUMENTS) cctk_iteration,i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),r(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), & - det,dens(i,j,k),scon(i,j,k,1), scon(i,j,k,2), scon(i,j,k,3), & + sdet,dens(i,j,k),scon(i,j,k,1), scon(i,j,k,2), scon(i,j,k,3), & tau(i,j,k), rho(i,j,k), velx(i,j,k), vely(i,j,k), & velz(i,j,k), eps(i,j,k), press(i,j,k), w_lorentz(i,j,k),& temperature(i,j,k),y_e(i,j,k)) @@ -493,7 +491,7 @@ subroutine GRHydro_InitialAtmosphereReset(CCTK_ARGUMENTS) rho(i,j,k),xeps,xtemp,xye,press(i,j,k),eps(i,j,k),keyerr,anyerr) call prim2conpolytype(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, & + 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), rho(i,j,k), velx(i,j,k), vely(i,j,k), & velz(i,j,k), eps(i,j,k), press(i,j,k), w_lorentz(i,j,k)) @@ -506,8 +504,8 @@ subroutine GRHydro_InitialAtmosphereReset(CCTK_ARGUMENTS) vely_p(i,j,k) = 0.0d0 velz_p(i,j,k) = 0.0d0 - det = SPATIAL_DETERMINANT(g11_p(i,j,k), g12_p(i,j,k), g13_p(i,j,k), \ - g22_p(i,j,k), g23_p(i,j,k), g33_p(i,j,k)) + sdet = sqrt(SPATIAL_DETERMINANT(g11_p(i,j,k), g12_p(i,j,k), g13_p(i,j,k), \ + g22_p(i,j,k), g23_p(i,j,k), g33_p(i,j,k))) if(evolve_temper.ne.0) then ! set the temperature to be relatively low @@ -522,7 +520,7 @@ subroutine GRHydro_InitialAtmosphereReset(CCTK_ARGUMENTS) cctk_iteration,i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),r(i,j,k),& g11_p(i,j,k),g12_p(i,j,k),& g13_p(i,j,k),g22_p(i,j,k),g23_p(i,j,k),g33_p(i,j,k), & - det,dens_p(i,j,k),scon_p(i,j,k,1), scon_p(i,j,k,2), scon_p(i,j,k,3), & + sdet,dens_p(i,j,k),scon_p(i,j,k,1), scon_p(i,j,k,2), scon_p(i,j,k,3), & tau_p(i,j,k), rho_p(i,j,k), velx_p(i,j,k), vely_p(i,j,k), & velz_p(i,j,k), eps_p(i,j,k), press_p(i,j,k), w_lorentz_p(i,j,k),& temperature_p(i,j,k),y_e_p(i,j,k)) @@ -535,7 +533,7 @@ subroutine GRHydro_InitialAtmosphereReset(CCTK_ARGUMENTS) rho_p(i,j,k),xeps,xtemp,xye,press_p(i,j,k),eps_p(i,j,k),keyerr,anyerr) call prim2conpolytype(eos_handle, & g11_p(i,j,k), g12_p(i,j,k), g13_p(i,j,k), & - g22_p(i,j,k), g23_p(i,j,k), g33_p(i,j,k), det, & + g22_p(i,j,k), g23_p(i,j,k), g33_p(i,j,k), sdet, & dens_p(i,j,k), scon_p(i,j,k,1), scon_p(i,j,k,2), scon_p(i,j,k,3), & tau_p(i,j,k), rho_p(i,j,k), velx_p(i,j,k), vely_p(i,j,k), & velz_p(i,j,k), eps_p(i,j,k), press_p(i,j,k), w_lorentz_p(i,j,k)) @@ -549,8 +547,8 @@ subroutine GRHydro_InitialAtmosphereReset(CCTK_ARGUMENTS) velx_p_p(i,j,k) = 0.0d0 vely_p_p(i,j,k) = 0.0d0 velz_p_p(i,j,k) = 0.0d0 - det = SPATIAL_DETERMINANT(g11_p_p(i,j,k), g12_p_p(i,j,k), g13_p_p(i,j,k), \ - g22_p_p(i,j,k), g23_p_p(i,j,k), g33_p_p(i,j,k)) + sdet = sqrt(SPATIAL_DETERMINANT(g11_p_p(i,j,k), g12_p_p(i,j,k), g13_p_p(i,j,k), \ + g22_p_p(i,j,k), g23_p_p(i,j,k), g33_p_p(i,j,k))) if(evolve_temper.ne.0) then ! set the temperature to be relatively low @@ -564,7 +562,7 @@ subroutine GRHydro_InitialAtmosphereReset(CCTK_ARGUMENTS) cctk_iteration,i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),r(i,j,k),& g11_p_p(i,j,k),g12_p_p(i,j,k),& g13_p_p(i,j,k),g22_p_p(i,j,k),g23_p_p(i,j,k),g33_p_p(i,j,k), & - det,dens_p_p(i,j,k),scon_p_p(i,j,k,1), scon_p_p(i,j,k,2), scon_p_p(i,j,k,3), & + sdet,dens_p_p(i,j,k),scon_p_p(i,j,k,1), scon_p_p(i,j,k,2), scon_p_p(i,j,k,3), & tau_p_p(i,j,k), rho_p_p(i,j,k), velx_p_p(i,j,k), vely_p_p(i,j,k), & velz_p_p(i,j,k), eps_p_p(i,j,k), press_p_p(i,j,k), w_lorentz_p_p(i,j,k),& temperature_p_p(i,j,k),y_e_p_p(i,j,k)) @@ -577,7 +575,7 @@ subroutine GRHydro_InitialAtmosphereReset(CCTK_ARGUMENTS) rho_p_p(i,j,k),xeps,xtemp,xye,press_p_p(i,j,k),eps_p_p(i,j,k),keyerr,anyerr) call prim2conpolytype(eos_handle, & g11_p_p(i,j,k), g12_p_p(i,j,k), g13_p_p(i,j,k), & - g22_p_p(i,j,k), g23_p_p(i,j,k), g33_p_p(i,j,k), det, & + g22_p_p(i,j,k), g23_p_p(i,j,k), g33_p_p(i,j,k), sdet, & dens_p_p(i,j,k), scon_p_p(i,j,k,1), scon_p_p(i,j,k,2), scon_p_p(i,j,k,3), & tau_p_p(i,j,k), rho_p_p(i,j,k), velx_p_p(i,j,k), vely_p_p(i,j,k), & velz_p_p(i,j,k), eps_p_p(i,j,k), press_p_p(i,j,k), w_lorentz_p_p(i,j,k)) diff --git a/src/GRHydro_UpdateMaskM.F90 b/src/GRHydro_UpdateMaskM.F90 index 84fdec5..2c16ccf 100644 --- a/src/GRHydro_UpdateMaskM.F90 +++ b/src/GRHydro_UpdateMaskM.F90 @@ -59,7 +59,7 @@ subroutine GRHydro_AtmosphereResetM(CCTK_ARGUMENTS) DECLARE_CCTK_PARAMETERS CCTK_INT :: i, j, k - CCTK_REAL :: det + CCTK_REAL :: sdet ! begin EOS Omni vars integer :: n,keytemp,anyerr,keyerr(1) @@ -107,7 +107,7 @@ subroutine GRHydro_AtmosphereResetM(CCTK_ARGUMENTS) if (verbose.eq.1) call CCTK_INFO("Entering AtmosphereReset.") -!$OMP PARALLEL DO PRIVATE(det,keytemp,i,j,k,anyerr,keyerr) +!$OMP PARALLEL DO PRIVATE(sdet,keytemp,i,j,k,anyerr,keyerr) do k = 1, cctk_lsh(3) do j = 1, cctk_lsh(2) do i = 1, cctk_lsh(1) @@ -118,8 +118,7 @@ subroutine GRHydro_AtmosphereResetM(CCTK_ARGUMENTS) vup(i,j,k,1) = 0.0d0 vup(i,j,k,2) = 0.0d0 vup(i,j,k,3) = 0.0d0 - 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 = sdetg(i,j,k) if(evolve_temper.ne.0) then @@ -133,7 +132,7 @@ subroutine GRHydro_AtmosphereResetM(CCTK_ARGUMENTS) call prim2conM_hot(GRHydro_eos_handle, GRHydro_reflevel,& i,j,k,x(i,j,k),y(i,j,k),z(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),& - det, dens(i,j,k),scon(i,j,k,1),scon(i,j,k,2),scon(i,j,k,3),& + 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),& eps(i,j,k),press(i,j,k),Bprim(i,j,k,1), & @@ -145,7 +144,7 @@ subroutine GRHydro_AtmosphereResetM(CCTK_ARGUMENTS) call prim2conpolytypeM(GRHydro_polytrope_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, & + 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), & @@ -199,7 +198,7 @@ subroutine GRHydro_InitialAtmosphereResetM(CCTK_ARGUMENTS) DECLARE_CCTK_PARAMETERS CCTK_INT :: i, j, k - CCTK_REAL :: det + CCTK_REAL :: sdet CCTK_REAL :: rho_min CCTK_INT :: eos_handle @@ -318,7 +317,7 @@ subroutine GRHydro_InitialAtmosphereResetM(CCTK_ARGUMENTS) rho_min = rho_min * initial_atmosphere_factor endif -!$OMP PARALLEL DO PRIVATE(det,keytemp,i,j,k,anyerr,keyerr) +!$OMP PARALLEL DO PRIVATE(sdet,keytemp,i,j,k,anyerr,keyerr) do k = 1, cctk_lsh(3) do j = 1, cctk_lsh(2) do i = 1, cctk_lsh(1) @@ -331,8 +330,7 @@ subroutine GRHydro_InitialAtmosphereResetM(CCTK_ARGUMENTS) vup(i,j,k,2) = 0.0d0 vup(i,j,k,3) = 0.0d0 - 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 = sdetg(i,j,k) if(evolve_temper.ne.0) then ! ! set the temperature to be relatively low @@ -345,7 +343,7 @@ subroutine GRHydro_InitialAtmosphereResetM(CCTK_ARGUMENTS) call prim2conM_hot(GRHydro_eos_handle, GRHydro_reflevel,& i,j,k,x(i,j,k),y(i,j,k),z(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),& - det, dens(i,j,k),scon(i,j,k,1),scon(i,j,k,2),scon(i,j,k,3),& + 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),& eps(i,j,k),press(i,j,k),Bprim(i,j,k,1), & @@ -363,7 +361,7 @@ subroutine GRHydro_InitialAtmosphereResetM(CCTK_ARGUMENTS) call prim2conpolytypeM(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, & + 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), & @@ -378,8 +376,8 @@ subroutine GRHydro_InitialAtmosphereResetM(CCTK_ARGUMENTS) vup_p(i,j,k,2) = 0.0d0 vup_p(i,j,k,3) = 0.0d0 - det = SPATIAL_DETERMINANT(g11_p(i,j,k), g12_p(i,j,k), g13_p(i,j,k), \ - g22_p(i,j,k), g23_p(i,j,k), g33_p(i,j,k)) + sdet = sqrt(SPATIAL_DETERMINANT(g11_p(i,j,k), g12_p(i,j,k), g13_p(i,j,k), \ + g22_p(i,j,k), g23_p(i,j,k), g33_p(i,j,k))) if(evolve_temper.ne.0) then ! ! set the temperature to be relatively low @@ -392,7 +390,7 @@ subroutine GRHydro_InitialAtmosphereResetM(CCTK_ARGUMENTS) call prim2conM_hot(GRHydro_eos_handle, GRHydro_reflevel,& i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),g11_p(i,j,k),& g12_p(i,j,k),g13_p(i,j,k),g22_p(i,j,k),g23_p(i,j,k),g33_p(i,j,k),& - det, dens_p(i,j,k),scon_p(i,j,k,1),scon_p(i,j,k,2),scon_p(i,j,k,3),& + sdet, dens_p(i,j,k),scon_p(i,j,k,1),scon_p(i,j,k,2),scon_p(i,j,k,3),& tau_p(i,j,k),Bcons_p(i,j,k,1),Bcons_p(i,j,k,2),Bcons_p(i,j,k,3),& rho_p(i,j,k),vup_p(i,j,k,1),vup_p(i,j,k,2),vup_p(i,j,k,3),& eps_p(i,j,k),press_p(i,j,k),Bprim_p(i,j,k,1), & @@ -410,7 +408,7 @@ subroutine GRHydro_InitialAtmosphereResetM(CCTK_ARGUMENTS) call prim2conpolytypeM(eos_handle, & g11_p(i,j,k), g12_p(i,j,k), g13_p(i,j,k), & - g22_p(i,j,k), g23_p(i,j,k), g33_p(i,j,k), det, & + g22_p(i,j,k), g23_p(i,j,k), g33_p(i,j,k), sdet, & dens_p(i,j,k), scon_p(i,j,k,1), scon_p(i,j,k,2), scon_p(i,j,k,3), & tau_p(i,j,k), Bcons_p(i,j,k,1),Bcons_p(i,j,k,2),Bcons_p(i,j,k,3),& rho_p(i,j,k), vup_p(i,j,k,1), vup_p(i,j,k,2), & @@ -427,8 +425,8 @@ subroutine GRHydro_InitialAtmosphereResetM(CCTK_ARGUMENTS) vup_p_p(i,j,k,2) = 0.0d0 vup_p_p(i,j,k,3) = 0.0d0 - det = SPATIAL_DETERMINANT(g11_p_p(i,j,k), g12_p_p(i,j,k), g13_p_p(i,j,k), \ - g22_p_p(i,j,k), g23_p_p(i,j,k), g33_p_p(i,j,k)) + sdet = sqrt(SPATIAL_DETERMINANT(g11_p_p(i,j,k), g12_p_p(i,j,k), g13_p_p(i,j,k), \ + g22_p_p(i,j,k), g23_p_p(i,j,k), g33_p_p(i,j,k))) if(evolve_temper.ne.0) then ! ! set the temperature to be relatively low @@ -441,7 +439,7 @@ subroutine GRHydro_InitialAtmosphereResetM(CCTK_ARGUMENTS) call prim2conM_hot(GRHydro_eos_handle, GRHydro_reflevel,& i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),g11_p_p(i,j,k),& g12_p_p(i,j,k),g13_p_p(i,j,k),g22_p_p(i,j,k),g23_p_p(i,j,k),g33_p_p(i,j,k),& - det, dens_p_p(i,j,k),scon_p_p(i,j,k,1),scon_p_p(i,j,k,2),scon_p_p(i,j,k,3),& + sdet, dens_p_p(i,j,k),scon_p_p(i,j,k,1),scon_p_p(i,j,k,2),scon_p_p(i,j,k,3),& tau_p_p(i,j,k),Bcons_p_p(i,j,k,1),Bcons_p_p(i,j,k,2),Bcons_p_p(i,j,k,3),& rho_p_p(i,j,k),vup_p_p(i,j,k,1),vup_p_p(i,j,k,2),vup_p_p(i,j,k,3),& eps_p_p(i,j,k),press_p_p(i,j,k),Bprim_p_p(i,j,k,1), & @@ -459,7 +457,7 @@ subroutine GRHydro_InitialAtmosphereResetM(CCTK_ARGUMENTS) call prim2conpolytypeM(eos_handle, & g11_p_p(i,j,k), g12_p_p(i,j,k), g13_p_p(i,j,k), & - g22_p_p(i,j,k), g23_p_p(i,j,k), g33_p_p(i,j,k), det, & + g22_p_p(i,j,k), g23_p_p(i,j,k), g33_p_p(i,j,k), sdet, & dens_p_p(i,j,k), scon_p_p(i,j,k,1), scon_p_p(i,j,k,2), scon_p_p(i,j,k,3), & tau_p_p(i,j,k), Bcons_p_p(i,j,k,1),Bcons_p_p(i,j,k,2),Bcons_p_p(i,j,k,3),& rho_p_p(i,j,k), vup_p_p(i,j,k,1), vup_p_p(i,j,k,2), & @@ -534,7 +532,7 @@ subroutine GRHydro_AtmosphereResetAM(CCTK_ARGUMENTS) DECLARE_CCTK_PARAMETERS CCTK_INT :: i, j, k - CCTK_REAL :: det + CCTK_REAL :: sdet ! begin EOS Omni vars integer :: n,keytemp,anyerr,keyerr(1) @@ -582,7 +580,7 @@ subroutine GRHydro_AtmosphereResetAM(CCTK_ARGUMENTS) if (verbose.eq.1) call CCTK_INFO("Entering AtmosphereReset.") -!$OMP PARALLEL DO PRIVATE(det,keytemp,i,j,k,anyerr,keyerr) +!$OMP PARALLEL DO PRIVATE(sdet,keytemp,i,j,k,anyerr,keyerr) do k = 1, cctk_lsh(3) do j = 1, cctk_lsh(2) do i = 1, cctk_lsh(1) @@ -593,8 +591,7 @@ subroutine GRHydro_AtmosphereResetAM(CCTK_ARGUMENTS) vup(i,j,k,1) = 0.0d0 vup(i,j,k,2) = 0.0d0 vup(i,j,k,3) = 0.0d0 - 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 = sdetg(i,j,k) if(evolve_temper.ne.0) then @@ -608,7 +605,7 @@ subroutine GRHydro_AtmosphereResetAM(CCTK_ARGUMENTS) call prim2conM_hot(GRHydro_eos_handle, GRHydro_reflevel,& i,j,k,x(i,j,k),y(i,j,k),z(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),& - det, dens(i,j,k),scon(i,j,k,1),scon(i,j,k,2),scon(i,j,k,3),& + 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),& eps(i,j,k),press(i,j,k),Bprim(i,j,k,1), & @@ -620,7 +617,7 @@ subroutine GRHydro_AtmosphereResetAM(CCTK_ARGUMENTS) call prim2conpolytypeM(GRHydro_polytrope_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, & + 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), & @@ -674,7 +671,7 @@ subroutine GRHydro_InitialAtmosphereResetAM(CCTK_ARGUMENTS) DECLARE_CCTK_PARAMETERS CCTK_INT :: i, j, k - CCTK_REAL :: det + CCTK_REAL :: sdet CCTK_REAL :: rho_min CCTK_INT :: eos_handle @@ -793,7 +790,7 @@ subroutine GRHydro_InitialAtmosphereResetAM(CCTK_ARGUMENTS) rho_min = rho_min * initial_atmosphere_factor endif -!$OMP PARALLEL DO PRIVATE(det,keytemp,i,j,k,anyerr,keyerr) +!$OMP PARALLEL DO PRIVATE(sdet,keytemp,i,j,k,anyerr,keyerr) do k = 1, cctk_lsh(3) do j = 1, cctk_lsh(2) do i = 1, cctk_lsh(1) @@ -806,8 +803,7 @@ subroutine GRHydro_InitialAtmosphereResetAM(CCTK_ARGUMENTS) vup(i,j,k,2) = 0.0d0 vup(i,j,k,3) = 0.0d0 - 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 = sdetg(i,j,k) if(evolve_temper.ne.0) then ! ! set the temperature to be relatively low @@ -820,7 +816,7 @@ subroutine GRHydro_InitialAtmosphereResetAM(CCTK_ARGUMENTS) call prim2conM_hot(GRHydro_eos_handle, GRHydro_reflevel,& i,j,k,x(i,j,k),y(i,j,k),z(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),& - det, dens(i,j,k),scon(i,j,k,1),scon(i,j,k,2),scon(i,j,k,3),& + 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),& eps(i,j,k),press(i,j,k),Bprim(i,j,k,1), & @@ -838,7 +834,7 @@ subroutine GRHydro_InitialAtmosphereResetAM(CCTK_ARGUMENTS) call prim2conpolytypeM(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, & + 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), & @@ -853,8 +849,8 @@ subroutine GRHydro_InitialAtmosphereResetAM(CCTK_ARGUMENTS) vup_p(i,j,k,2) = 0.0d0 vup_p(i,j,k,3) = 0.0d0 - det = SPATIAL_DETERMINANT(g11_p(i,j,k), g12_p(i,j,k), g13_p(i,j,k), \ - g22_p(i,j,k), g23_p(i,j,k), g33_p(i,j,k)) + sdet = sqrt(SPATIAL_DETERMINANT(g11_p(i,j,k), g12_p(i,j,k), g13_p(i,j,k), \ + g22_p(i,j,k), g23_p(i,j,k), g33_p(i,j,k))) if(evolve_temper.ne.0) then ! ! set the temperature to be relatively low @@ -867,7 +863,7 @@ subroutine GRHydro_InitialAtmosphereResetAM(CCTK_ARGUMENTS) call prim2conM_hot(GRHydro_eos_handle, GRHydro_reflevel,& i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),g11_p(i,j,k),& g12_p(i,j,k),g13_p(i,j,k),g22_p(i,j,k),g23_p(i,j,k),g33_p(i,j,k),& - det, dens_p(i,j,k),scon_p(i,j,k,1),scon_p(i,j,k,2),scon_p(i,j,k,3),& + sdet, dens_p(i,j,k),scon_p(i,j,k,1),scon_p(i,j,k,2),scon_p(i,j,k,3),& tau_p(i,j,k),Bcons_p(i,j,k,1),Bcons_p(i,j,k,2),Bcons_p(i,j,k,3),& rho_p(i,j,k),vup_p(i,j,k,1),vup_p(i,j,k,2),vup_p(i,j,k,3),& eps_p(i,j,k),press_p(i,j,k),Bprim_p(i,j,k,1), & @@ -885,7 +881,7 @@ subroutine GRHydro_InitialAtmosphereResetAM(CCTK_ARGUMENTS) call prim2conpolytypeM(eos_handle, & g11_p(i,j,k), g12_p(i,j,k), g13_p(i,j,k), & - g22_p(i,j,k), g23_p(i,j,k), g33_p(i,j,k), det, & + g22_p(i,j,k), g23_p(i,j,k), g33_p(i,j,k), sdet, & dens_p(i,j,k), scon_p(i,j,k,1), scon_p(i,j,k,2), scon_p(i,j,k,3), & tau_p(i,j,k), Bcons_p(i,j,k,1),Bcons_p(i,j,k,2),Bcons_p(i,j,k,3),& rho_p(i,j,k), vup_p(i,j,k,1), vup_p(i,j,k,2), & @@ -902,8 +898,8 @@ subroutine GRHydro_InitialAtmosphereResetAM(CCTK_ARGUMENTS) vup_p_p(i,j,k,2) = 0.0d0 vup_p_p(i,j,k,3) = 0.0d0 - det = SPATIAL_DETERMINANT(g11_p_p(i,j,k), g12_p_p(i,j,k), g13_p_p(i,j,k), \ - g22_p_p(i,j,k), g23_p_p(i,j,k), g33_p_p(i,j,k)) + sdet = sqrt(SPATIAL_DETERMINANT(g11_p_p(i,j,k), g12_p_p(i,j,k), g13_p_p(i,j,k), \ + g22_p_p(i,j,k), g23_p_p(i,j,k), g33_p_p(i,j,k))) if(evolve_temper.ne.0) then ! ! set the temperature to be relatively low @@ -916,7 +912,7 @@ subroutine GRHydro_InitialAtmosphereResetAM(CCTK_ARGUMENTS) call prim2conM_hot(GRHydro_eos_handle, GRHydro_reflevel,& i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),g11_p_p(i,j,k),& g12_p_p(i,j,k),g13_p_p(i,j,k),g22_p_p(i,j,k),g23_p_p(i,j,k),g33_p_p(i,j,k),& - det, dens_p_p(i,j,k),scon_p_p(i,j,k,1),scon_p_p(i,j,k,2),scon_p_p(i,j,k,3),& + sdet, dens_p_p(i,j,k),scon_p_p(i,j,k,1),scon_p_p(i,j,k,2),scon_p_p(i,j,k,3),& tau_p_p(i,j,k),Bcons_p_p(i,j,k,1),Bcons_p_p(i,j,k,2),Bcons_p_p(i,j,k,3),& rho_p_p(i,j,k),vup_p_p(i,j,k,1),vup_p_p(i,j,k,2),vup_p_p(i,j,k,3),& eps_p_p(i,j,k),press_p_p(i,j,k),Bprim_p_p(i,j,k,1), & @@ -934,7 +930,7 @@ subroutine GRHydro_InitialAtmosphereResetAM(CCTK_ARGUMENTS) call prim2conpolytypeM(eos_handle, & g11_p_p(i,j,k), g12_p_p(i,j,k), g13_p_p(i,j,k), & - g22_p_p(i,j,k), g23_p_p(i,j,k), g33_p_p(i,j,k), det, & + g22_p_p(i,j,k), g23_p_p(i,j,k), g33_p_p(i,j,k), sdet, & dens_p_p(i,j,k), scon_p_p(i,j,k,1), scon_p_p(i,j,k,2), scon_p_p(i,j,k,3), & tau_p_p(i,j,k), Bcons_p_p(i,j,k,1),Bcons_p_p(i,j,k,2),Bcons_p_p(i,j,k,3),& rho_p_p(i,j,k), vup_p_p(i,j,k,1), vup_p_p(i,j,k,2), & diff --git a/src/Utils.F90 b/src/Utils.F90 index d8b29e4..9f670e2 100644 --- a/src/Utils.F90 +++ b/src/Utils.F90 @@ -71,21 +71,34 @@ subroutine GRHydro_RefinementLevel(CCTK_ARGUMENTS) end subroutine GRHydro_RefinementLevel +subroutine GRHydro_SqrtSpatialDeterminant(CCTK_ARGUMENTS) + + implicit none + DECLARE_CCTK_ARGUMENTS + integer i,j,k + integer nx, ny, nz + + nx = cctk_lsh(1) + ny = cctk_lsh(2) + nz = cctk_lsh(3) + + !$OMP PARALLEL DO + do k=1,nz + do j=1,ny + do i=1,nx + sdetg(i,j,k) = -(gxz(i,j,k)**2)*gyy(i,j,k) + & + 2.0d0*gxy(i,j,k)*gxz(i,j,k)*gyz(i,j,k) - & + gxx(i,j,k)*(gyz(i,j,k)**2) - & + (gxy(i,j,k)**2)*gzz(i,j,k) + & + (gxx(i,j,k)*gyy(i,j,k))*gzz(i,j,k) + sdetg(i,j,k) = sqrt(sdetg(i,j,k)) + enddo + enddo + enddo + !$OMP END PARALLEL DO + +end subroutine GRHydro_SqrtSpatialDeterminant - /*@@ - @routine SpatialDeterminant - @date Sat Jan 26 02:30:23 2002 - @author - @desc - Calculates the determinant of the spatial metric. - PLEASE USE THE MACRO SPATIAL_DETERMINANT from now on!!! - @enddesc - @calls - @calledby - @history - - @endhistory -@@*/ subroutine SpatialDeterminant(gxx,gxy,gxz,gyy,gyz,gzz,det) @@ -145,51 +158,4 @@ subroutine UpperMetric(uxx, uxy, uxz, uyy, uyz, uzz, & end subroutine UpperMetric - /*@@ - @routine SetMetricTemps - @date Wed Mar 23 10:01:08 2005 - @author Ian Hawke - @desc - If using the new EOS store the metric determinant and inverse metric - to avoid repeated computation - @enddesc - @calls - @calledby - @history - - @endhistory - -@@*/ - -subroutine SetMetricTemps(CCTK_ARGUMENTS) - - implicit none - - DECLARE_CCTK_ARGUMENTS - DECLARE_CCTK_PARAMETERS - - CCTK_INT :: i, j, k, nx, ny, nz - - nx = cctk_lsh(1) - ny = cctk_lsh(2) - nz = cctk_lsh(3) - - do k = 1, nz - do j = 1, ny - do i = 1, nx - - GRHydro_det(i,j,k) = SPATIAL_DETERMINANT(gxx(i,j,k),gxy(i,j,k),\ - gxz(i,j,k), gyy(i,j,k),gyz(i,j,k),gzz(i,j,k)) - call UpperMetric(& - GRHydro_uxx(i,j,k),GRHydro_uxy(i,j,k),GRHydro_uxz(i,j,k),& - GRHydro_uyy(i,j,k),GRHydro_uyz(i,j,k),GRHydro_uzz(i,j,k),& - GRHydro_det(i,j,k),& - gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),& - gyy(i,j,k),gyz(i,j,k),gzz(i,j,k)) - - end do - end do - end do - -end subroutine SetMetricTemps |