/*@@ @file GRHydro_RegisterVars.c @date Sat Jan 26 01:06:01 2002 @author The GRHydro Developers @desc The routines for converting conservative to primitive variables. @enddesc @@*/ #include "cctk.h" #include "cctk_Parameters.h" #include "cctk_Arguments.h" #include "cctk_Functions.h" #include "SpaceMask.h" module Con2Prim_fortran_interfaces implicit none interface subroutine Con2Prim_pt(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, & GRHydro_reflevel, GRHydro_C2P_failed) implicit none CCTK_INT handle CCTK_REAL dens CCTK_REAL sx, sy, sz CCTK_REAL tau CCTK_REAL rho CCTK_REAL velx, vely, velz CCTK_REAL epsilon, press CCTK_REAL w_lorentz CCTK_REAL uxx, uxy, uxz CCTK_REAL uyy, uyz, uzz CCTK_REAL det CCTK_REAL x, y, z, r logical epsnegative CCTK_REAL GRHydro_rho_min, pmin, epsmin CCTK_INT GRHydro_reflevel CCTK_REAL GRHydro_C2P_failed end subroutine Con2Prim_pt subroutine Con2Prim_ptPolytype(GRHydro_polytrope_handle, & dens, & sx, sy, sz, & tau, & rho, & velx, vely, velz, & eps, press, & w_lorentz, & uxx, uxy, uxz, uyy, uyz, uzz, & det, & x, y, z, r, & GRHydro_rho_min, & GRHydro_reflevel, GRHydro_C2P_failed) implicit none CCTK_INT GRHydro_polytrope_handle CCTK_REAL dens CCTK_REAL sx, sy, sz CCTK_REAL tau CCTK_REAL rho CCTK_REAL velx, vely, velz CCTK_REAL eps, press CCTK_REAL w_lorentz CCTK_REAL uxx, uxy, uxz CCTK_REAL uyy, uyz, uzz CCTK_REAL det CCTK_REAL x, y, z, r CCTK_REAL GRHydro_rho_min CCTK_INT GRHydro_reflevel CCTK_REAL GRHydro_C2P_failed end subroutine Con2Prim_ptPolytype subroutine Con2Prim_ptTracer(cons_tracer, tracer, dens) implicit none CCTK_REAL cons_tracer, tracer, dens end subroutine Con2Prim_ptTracer subroutine Con2Prim_ptBoundsTracer(cons_tracer, tracer, rho, one_over_w_lorentz, det) implicit none CCTK_REAL cons_tracer, tracer, rho, one_over_w_lorentz, det end subroutine Con2Prim_ptBoundsTracer end interface end module Con2Prim_fortran_interfaces /*@@ @routine Conservative2Primitive @date Sat Jan 26 01:08:33 2002 @author Ian Hawke @desc Wrapper routine that converts from conserved to primitive variables at every grid cell centre. @enddesc @calls @calledby @history Trimmed and altered from the GR3D routines, original author Mark Miller. 2007?: Bruno escluded the points in the atmosphere and excision region from the computation. Aug. 2008: Luca added a check on whether a failure at a given point may be disregarded, because that point will then be restricted from a finer level. This should be completely safe only if *regridding happens at times when all levels are evolved.* Feb. 2009: The above procedure proved to be wrong, so Luca implemented another one. When a failure occurs, it is temporarily ignored, except for storing the location of where it occured in a mask. At the end, after a Carpet restriction, the mask is checked and if it still contains failures, the run is aborted with an error message. Only used routines have been updated to use this procedure. @endhistory @@*/ subroutine Conservative2Primitive(CCTK_ARGUMENTS) use Con2Prim_fortran_interfaces implicit none DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS #ifdef _EOS_BASE_INC_ #undef _EOS_BASE_INC_ #endif #include "EOS_Base.inc" integer :: i, j, k, itracer, nx, ny, nz CCTK_REAL :: uxx, uxy, uxz, uyy, uyz, uzz, det, psi4pt, pmin, epsmin logical :: epsnegative character(len=100) warnline CCTK_INT :: type_bits, atmosphere CCTK_INT :: type2_bits CCTK_REAL :: local_min_tracer call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere") call SpaceMask_GetStateBits(atmosphere, "Hydro_Atmosphere", "in_atmosphere") type2_bits = -1 nx = cctk_lsh(1) ny = cctk_lsh(2) nz = cctk_lsh(3) if (use_min_tracer .ne. 0) then local_min_tracer = min_tracer else local_min_tracer = 0d0 end if pmin = EOS_Pressure(GRHydro_polytrope_handle, GRHydro_rho_min, 1.0d0) epsmin = EOS_SpecificIntEnergy(GRHydro_polytrope_handle, GRHydro_rho_min, pmin) !$OMP PARALLEL DO PRIVATE(i,j,itracer,& !$OMP uxx, uxy, uxz, uyy, uyz, uzz, det, psi4pt, epsnegative) do k = 1, nz do j = 1, ny do i = 1, nx !do not compute if in atmosphere or in excised region if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere) .or. & (hydro_excision_mask(i,j,k) .ne. 0)) cycle epsnegative = .false. if (conformal_state .eq. 0) then call SpatialDeterminant(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) call UpperMetric(uxx,uxy,uxz,uyy,uyz,uzz,det,& gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),& gyz(i,j,k),gzz(i,j,k)) else psi4pt = psi(i,j,k)**4 call SpatialDeterminant(psi4pt*gxx(i,j,k),psi4pt*gxy(i,j,k),& psi4pt*gxz(i,j,k),psi4pt*gyy(i,j,k),& psi4pt*gyz(i,j,k),psi4pt*gzz(i,j,k),det) call UpperMetric(uxx,uxy,uxz,uyy,uyz,uzz,det,& psi4pt*gxx(i,j,k),psi4pt*gxy(i,j,k),psi4pt*gxz(i,j,k),& psi4pt*gyy(i,j,k),psi4pt*gyz(i,j,k),psi4pt*gzz(i,j,k)) end if !!$ if (det < 0.e0) then !!$ call CCTK_WARN(0, "nan produced (det negative)") !!$ end if if (evolve_tracer .ne. 0) then do itracer=1,number_of_tracers call Con2Prim_ptTracer(cons_tracer(i,j,k,itracer), tracer(i,j,k,itracer), & dens(i,j,k)) if (use_min_tracer .ne. 0) then if (tracer(i,j,k,itracer) .le. local_min_tracer) then tracer(i,j,k,itracer) = local_min_tracer end if end if enddo endif if ( dens(i,j,k) .le. sqrt(det)*GRHydro_rho_min*(1.d0+GRHydro_atmo_tolerance) ) then dens(i,j,k) = sqrt(det)*GRHydro_rho_min !/(1.d0+GRHydro_atmo_tolerance) rho(i,j,k) = GRHydro_rho_min scon(i,j,k,:) = 0.d0 vel(i,j,k,:) = 0.d0 w_lorentz(i,j,k) = 1.d0 press(i,j,k) = EOS_Pressure(GRHydro_polytrope_handle, rho(i,j,k), eps(i,j,k)) eps(i,j,k) = EOS_SpecificIntEnergy(GRHydro_polytrope_handle, rho(i,j,k), press(i,j,k)) ! 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) cycle end if call Con2Prim_pt(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),vel(i,j,k,1),vel(i,j,k,2), & vel(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), & z(i,j,k),r(i,j,k),epsnegative,GRHydro_rho_min,pmin, epsmin, GRHydro_reflevel, GRHydro_C2P_failed(i,j,k)) if (epsnegative) then #if 1 ! cott 2010/03/30: ! Set point to atmosphere, but continue evolution -- this is better than using ! the poly EOS -- it will lead the code to crash if this happens inside a (neutron) star, ! but will allow the job to continue if it happens in the atmosphere or in a ! zone that contains garbage (i.e. boundary, buffer zones) ! Ultimately, we want this fixed via a new carpet mask presently under development ! GRHydro_C2P_failed(i,j,k) = 1 !$OMP CRITICAL call CCTK_WARN(GRHydro_NaN_verbose+2, 'Specific internal energy just went below 0! ') write(warnline,'(a28,i2)') 'on carpet reflevel: ',GRHydro_reflevel call CCTK_WARN(GRHydro_NaN_verbose+2,warnline) write(warnline,'(a20,3g16.7)') 'xyz location: ',& x(i,j,k),y(i,j,k),z(i,j,k) call CCTK_WARN(GRHydro_NaN_verbose+2,warnline) write(warnline,'(a20,g16.7)') 'radius: ',r(i,j,k) call CCTK_WARN(GRHydro_NaN_verbose+2,warnline) call CCTK_WARN(GRHydro_NaN_verbose+2,"Setting the point to atmosphere") !$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) rho(i,j,k) = GRHydro_rho_min scon(i,j,k,:) = 0.d0 vel(i,j,k,:) = 0.d0 w_lorentz(i,j,k) = 1.d0 press(i,j,k) = EOS_Pressure(GRHydro_polytrope_handle, rho(i,j,k), eps(i,j,k)) eps(i,j,k) = EOS_SpecificIntEnergy(GRHydro_polytrope_handle, rho(i,j,k), press(i,j,k)) ! 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) #else ! cott 2010/03/27: ! Honestly, this should never happen. We need to flag the point where ! this happened as having led to failing con2prim. !$OMP CRITICAL call CCTK_WARN(GRHydro_NaN_verbose+2, 'Specific internal energy just went below 0, trying polytype.') !$OMP END CRITICAL call Con2Prim_ptPolytype(GRHydro_polytrope_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), vel(i,j,k,1), vel(i,j,k,2), & vel(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), & z(i,j,k), r(i,j,k),GRHydro_rho_min, GRHydro_reflevel, GRHydro_C2P_failed(i,j,k)) #endif end if end do end do end do !$OMP END PARALLEL DO return end subroutine Conservative2Primitive /*@@ @routine Con2Prim_pt @date Sat Jan 26 01:09:39 2002 @author The GRHydro Develpoers @desc Given all the appropriate data, recover the primitive variables at a single point. @enddesc @calls @calledby @history @endhistory @@*/ subroutine Con2Prim_pt(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, 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, & y, z, r, GRHydro_rho_min CCTK_REAL s2, c0, c1, c2, c3, c4, f, df, ftol, v2, w, vlowx, vlowy, vlowz CCTK_INT count, i, handle, GRHydro_reflevel, GRHydro_C2P_failed CCTK_REAL udens, usx, usy, usz, utau, pold, pnew, epsold, epsnew, w2, & w2mhalf, temp1, drhobydpress, depsbydpress, dpressbydeps, dpressbydrho, pmin, epsmin character(len=200) warnline logical epsnegative #ifdef _EOS_BASE_INC_ #undef _EOS_BASE_INC_ #endif #include "EOS_Base.inc" !!$ Undensitize the variables udens = dens /sqrt(det) usx = sx /sqrt(det) usy = sy /sqrt(det) usz = sz /sqrt(det) utau = tau /sqrt(det) s2 = usx*usx*uxx + usy*usy*uyy + usz*usz*uzz + 2.*usx*usy*uxy + & 2.*usx*usz*uxz + 2.*usy*usz*uyz !!$ Set initial guess for pressure: pold = max(1.d-10,EOS_Pressure(handle, rho, epsilon)) !!$ Check that the variables have a chance of being physical if( (utau + pold + udens)**2 - s2 .le. 0.0d0) then if (c2p_reset_pressure .ne. 0) then pold = sqrt(s2 + c2p_reset_pressure_to_value) - utau - udens else !$OMP CRITICAL call CCTK_WARN(GRHydro_NaN_verbose, "c2p failed and being told not to reset the pressure") GRHydro_C2P_failed = 1 !$OMP END CRITICAL endif endif !!$ Calculate rho and epsilon !define temporary variables to speed up rho = udens * sqrt( (utau + pold + udens)**2 - s2)/(utau + pold + udens) w_lorentz = (utau + pold + udens) / sqrt( (utau + pold + udens)**2 - s2) epsilon = (sqrt( (utau + pold + udens)**2 - s2) - pold * w_lorentz - & udens)/udens !!$ BEGIN: Check for NaN value (1st check) if (rho .ne. rho) then !$OMP CRITICAL write(warnline,'(a70,7g15.6)') 'NaN produced in sqrt(): (utau, pold, udens, s2, x, y, z)', utau, pold, udens, s2, x, y, z call CCTK_WARN(GRHydro_NaN_verbose, warnline) !$OMP END CRITICAL endif !!$ END: Check for NaN value (1st check) !!$ Calculate the function f = pold - EOS_Pressure(handle, rho, epsilon) !!$Find the root count = 0 pnew = pold do while ( ((abs(pnew - pold)/abs(pnew) .gt. GRHydro_perc_ptol) .and. & (abs(pnew - pold) .gt. GRHydro_del_ptol)) .or. & (count .lt. GRHydro_countmin)) count = count + 1 if (count > GRHydro_countmax) then GRHydro_C2P_failed = 1 !$OMP CRITICAL call CCTK_WARN(1, 'count > GRHydro_countmax! ') write(warnline,'(a28,i2)') 'on carpet reflevel: ',GRHydro_reflevel call CCTK_WARN(1,warnline) write(warnline,'(a20,3g16.7)') 'xyz location: ',x,y,z call CCTK_WARN(1,warnline) write(warnline,'(a20,g16.7)') 'radius: ',r call CCTK_WARN(1,warnline) call CCTK_WARN(1,"Setting the point to atmosphere") !$OMP END CRITICAL ! for safety, let's set the point to atmosphere rho = GRHydro_rho_min udens = rho dens = sqrt(det) * rho pnew = pmin epsilon = epsmin ! w_lorentz=1, so the expression for utau reduces to: utau = rho + rho*epsmin - udens sx = 0.d0 sy = 0.d0 sz = 0.d0 s2 = 0.d0 usx = 0.d0 usy = 0.d0 usz = 0.d0 w_lorentz = 1.d0 goto 51 end if dpressbydrho = EOS_DPressByDRho(handle,rho,epsilon) dpressbydeps = EOS_DPressByDEps(handle,rho,epsilon) temp1 = (utau+udens+pnew)**2 - s2 drhobydpress = udens * s2 / (sqrt(temp1)*(udens+utau+pnew)**2) depsbydpress = pnew * s2 / (rho * (udens + utau + pnew) * temp1) df = 1.0d0 - dpressbydrho*drhobydpress - & dpressbydeps*depsbydpress pold = pnew pnew = max(pold - f/df, pmin) !!$ Recalculate primitive variables and function rho = udens * sqrt( (utau + pnew + udens)**2 - s2)/(utau + pnew + udens) w_lorentz = (utau + pnew + udens) / sqrt( (utau + pnew + udens)**2 - & s2) epsilon = (sqrt( (utau + pnew + udens)**2 - s2) - pnew * w_lorentz - & udens)/udens f = pnew - EOS_Pressure(handle, rho, epsilon) !check whether rho is NaN if (rho .ne. rho) then !$OMP CRITICAL write(warnline,'(a20,g20.7)') 'rho is NaN: ',rho call CCTK_WARN(GRHydro_NaN_verbose,warnline) !$OMP END CRITICAL goto 51 end if enddo !!$ Polish the root do i=1,GRHydro_polish dpressbydrho = EOS_DPressByDRho(handle,rho,epsilon) dpressbydeps = EOS_DPressByDEps(handle,rho,epsilon) temp1 = (utau+udens+pnew)**2 - s2 drhobydpress = udens * s2 / (sqrt(temp1)*(udens+utau+pnew)**2) depsbydpress = pnew * s2 / (rho * (udens + utau + pnew) * temp1) df = 1.0d0 - dpressbydrho*drhobydpress - & dpressbydeps*depsbydpress pold = pnew pnew = pold - f/df !!$ Recalculate primitive variables and function rho = udens * sqrt( (utau + pnew + udens)**2 - s2)/(utau + pnew + udens) w_lorentz = (utau + pnew + udens) / sqrt( (utau + pnew + udens)**2 - & s2) epsilon = (sqrt( (utau + pnew + udens)**2 - s2) - pnew * w_lorentz - & udens)/udens f = pold - EOS_Pressure(handle, rho, epsilon) enddo !!$ BEGIN: Check for NaN value (2nd check) if (rho .ne. rho) then !$OMP CRITICAL write(warnline,'(a70,7g15.6)') 'NaN produced in sqrt(): (utau, pold, udens, s2, x, y, z)', utau, pold, udens, s2, x, y, z call CCTK_WARN(GRHydro_NaN_verbose, warnline) !$OMP END CRITICAL endif !!$ END: Check for NaN value (2nd check) !!$ Calculate primitive variables from root if (rho .le. GRHydro_rho_min*(1.d0+GRHydro_atmo_tolerance) ) then rho = GRHydro_rho_min udens = rho dens = sqrt(det) * rho ! epsilon = (sqrt( (utau + pnew + udens)**2) - pnew - udens)/udens epsilon = epsmin ! w_lorentz=1, so the expression for utau reduces to: utau = rho + rho*epsmin - udens sx = 0.d0 sy = 0.d0 sz = 0.d0 s2 = 0.d0 usx = 0.d0 usy = 0.d0 usz = 0.d0 w_lorentz = 1.d0 end if 51 press = pnew vlowx = usx / ( (rho + rho*epsilon + press) * w_lorentz**2) vlowy = usy / ( (rho + rho*epsilon + press) * w_lorentz**2) vlowz = usz / ( (rho + rho*epsilon + press) * w_lorentz**2) velx = uxx * vlowx + uxy * vlowy + uxz * vlowz vely = uxy * vlowx + uyy * vlowy + uyz * vlowz velz = uxz * vlowx + uyz * vlowy + uzz * vlowz !!$If all else fails, use the polytropic EoS if(epsilon .lt. 0.0d0) then epsnegative = .true. endif end subroutine Con2Prim_pt /*@@ @routine Conservative2PrimitiveBoundaries @date Tue Mar 12 18:04:40 2002 @author The GRHydro Developers @desc This routine is used only if the reconstruction is performed on the conserved variables. It computes the primitive variables on cell boundaries. Since reconstruction on conservative had not proved to be very successful, some of the improvements to the C2P routines (e.g. the check about whether a failure happens in a point that will be restriced anyway) are not implemented here yet. @enddesc @calls @calledby @history @endhistory @@*/ subroutine Conservative2PrimitiveBounds(CCTK_ARGUMENTS) use Con2Prim_fortran_interfaces implicit none DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS integer :: i, j, k, itracer, nx, ny, nz CCTK_REAL :: uxxl, uxyl, uxzl, uyyl, uyzl, uzzl,& uxxr, uxyr, uxzr, uyyr, uyzr, uzzr, pmin, epsmin CCTK_REAL :: gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_detl,& gxxr,gxyr,gxzr,gyyr,gyzr,gzzr,avg_detr,psi4r,psi4l logical :: epsnegative character(len=100) warnline CCTK_INT :: type_bits, atmosphere CCTK_INT :: type2_bits CCTK_REAL :: local_min_tracer #ifdef _EOS_BASE_INC_ #undef _EOS_BASE_INC_ #endif #include "EOS_Base.inc" pmin=EOS_Pressure(GRHydro_polytrope_handle, GRHydro_rho_min, 1.0d0) epsmin = EOS_SpecificIntEnergy(GRHydro_polytrope_handle, GRHydro_rho_min, pmin) call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere") call SpaceMask_GetStateBits(atmosphere, "Hydro_Atmosphere", "in_atmosphere") type2_bits = -1 nx = cctk_lsh(1) ny = cctk_lsh(2) nz = cctk_lsh(3) if (use_min_tracer .ne. 0) then local_min_tracer = min_tracer else local_min_tracer = 0d0 end if do k = GRHydro_stencil, nz - GRHydro_stencil + 1 do j = GRHydro_stencil, ny - GRHydro_stencil + 1 do i = GRHydro_stencil, nx - GRHydro_stencil + 1 !do not compute if in atmosphere or in an excised region if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere) .or. & GRHydro_enable_internal_excision /= 0 .and. (hydro_excision_mask(i,j,k) .ne. 0)) cycle gxxl = 0.5d0 * (gxx(i,j,k) + gxx(i-xoffset,j-yoffset,k-zoffset)) gxyl = 0.5d0 * (gxy(i,j,k) + gxy(i-xoffset,j-yoffset,k-zoffset)) gxzl = 0.5d0 * (gxz(i,j,k) + gxz(i-xoffset,j-yoffset,k-zoffset)) gyyl = 0.5d0 * (gyy(i,j,k) + gyy(i-xoffset,j-yoffset,k-zoffset)) gyzl = 0.5d0 * (gyz(i,j,k) + gyz(i-xoffset,j-yoffset,k-zoffset)) gzzl = 0.5d0 * (gzz(i,j,k) + gzz(i-xoffset,j-yoffset,k-zoffset)) gxxr = 0.5d0 * (gxx(i,j,k) + gxx(i+xoffset,j+yoffset,k+zoffset)) gxyr = 0.5d0 * (gxy(i,j,k) + gxy(i+xoffset,j+yoffset,k+zoffset)) gxzr = 0.5d0 * (gxz(i,j,k) + gxz(i+xoffset,j+yoffset,k+zoffset)) gyyr = 0.5d0 * (gyy(i,j,k) + gyy(i+xoffset,j+yoffset,k+zoffset)) 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)) epsnegative = .false. if (conformal_state .eq. 0) then call SpatialDeterminant(gxxl,gxyl,gxzl,gyyl, gyzl,gzzl,avg_detl) call SpatialDeterminant(gxxr,gxyr,gxzr,gyyr, gyzr,gzzr,avg_detr) call UpperMetric(uxxl,uxyl,uxzl,uyyl,uyzl,uzzl,avg_detl,& gxxl,gxyl,gxzl,gyyl,gyzl,gzzl) call UpperMetric(uxxr,uxyr,uxzr,uyyr,uyzr,uzzr,avg_detr,& gxxr,gxyr,gxzr,gyyr,gyzr,gzzr) else psi4r = (0.5d0*(psi(i,j,k)+psi(i+xoffset,j+yoffset,k+zoffset)))**4 psi4l = (0.5d0*(psi(i,j,k)+psi(i-xoffset,j-yoffset,k-zoffset)))**4 call SpatialDeterminant(psi4l*gxxl,psi4l*gxyl,psi4l*gxzl,& psi4l*gyyl,psi4l*gyzl,psi4l*gzzl,avg_detl) call SpatialDeterminant(psi4r*gxxr,psi4r*gxyr,psi4r*gxzr,& psi4r*gyyr,psi4r*gyzr,psi4r*gzzr,avg_detr) call UpperMetric(uxxl,uxyl,uxzl,uyyl,uyzl,uzzl,avg_detl,& psi4l*gxxl,psi4l*gxyl,psi4l*gxzl,psi4l*gyyl,psi4l*gyzl,& psi4l*gzzl) call UpperMetric(uxxr,uxyr,uxzr,uyyr,uyzr,uzzr,avg_detr,& psi4r*gxxr,psi4r*gxyr,psi4r*gxzr,psi4r*gyyr,psi4r*gyzr,& psi4r*gzzr) end if if (evolve_tracer .ne. 0) then do itracer=1,number_of_tracers call Con2Prim_ptTracer(cons_tracer(i,j,k,itracer), & tracer(i,j,k,itracer), dens(i,j,k)) enddo if (use_min_tracer .ne. 0) then if (tracer(i,j,k,itracer) .le. local_min_tracer) then tracer(i,j,k,itracer) = local_min_tracer end if end if endif call Con2Prim_pt(GRHydro_eos_handle, 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),& 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),& epsnegative,GRHydro_rho_min,pmin, epsmin, GRHydro_reflevel, GRHydro_C2P_failed(i,j,k)) if (epsnegative) then !$OMP CRITICAL call CCTK_WARN(GRHydro_NaN_verbose+2, 'Specific internal energy just went below 0, trying polytype!') !$OMP END CRITICAL call Con2Prim_ptPolytype(GRHydro_polytrope_handle, 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),& 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)) end if if (epsminus(i,j,k) .lt. 0.0d0) then if (GRHydro_reflevel.ge.GRHydro_c2p_warn_from_reflevel) then !$OMP CRITICAL call CCTK_WARN(1,'Con2Prim: stopping the code.') call CCTK_WARN(1, ' specific internal energy just went below 0! ') write(warnline,'(a28,i2)') 'on carpet reflevel: ',GRHydro_reflevel call CCTK_WARN(1,warnline) write(warnline,'(a20,3g16.7)') 'xyz location: ',& x(i,j,k),y(i,j,k),z(i,j,k) call CCTK_WARN(1,warnline) write(warnline,'(a20,g16.7)') 'radius: ',r(i,j,k) call CCTK_WARN(1,warnline) write(warnline,'(a20,3g16.7)') 'velocities: ',& velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k) call CCTK_WARN(1,warnline) call CCTK_WARN(GRHydro_c2p_warnlevel, "Specific internal energy negative") !$OMP END CRITICAL exit endif endif epsnegative = .false. call Con2Prim_pt(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),& velyplus(i,j,k),velzplus(i,j,k),epsplus(i,j,k),& pressplus(i,j,k),w_lorentzplus(i,j,k),& uxxr,uxyr,uxzr,uyyr,uyzr,uzzr,avg_detr,& 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),& epsnegative,GRHydro_rho_min,pmin, epsmin, GRHydro_reflevel,GRHydro_C2P_failed(i,j,k)) if (epsnegative) then !$OMP CRITICAL call CCTK_WARN(GRHydro_NaN_verbose+2, 'Specific internal energy just went below 0, trying polytype!!') !$OMP END CRITICAL call Con2Prim_ptPolytype(GRHydro_polytrope_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),& velyplus(i,j,k),velzplus(i,j,k),epsplus(i,j,k),& pressplus(i,j,k),w_lorentzplus(i,j,k),& uxxr,uxyr,uxzr,uyyr,uyzr,uzzr,avg_detr,& 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)) end if if (epsplus(i,j,k) .lt. 0.0d0) then if (GRHydro_reflevel.ge.GRHydro_c2p_warn_from_reflevel) then !$OMP CRITICAL call CCTK_WARN(1,'Con2Prim: stopping the code.') call CCTK_WARN(1, ' specific internal energy just went below 0! ') write(warnline,'(a28,i2)') 'on carpet reflevel: ',GRHydro_reflevel call CCTK_WARN(1,warnline) write(warnline,'(a20,3g16.7)') 'xyz location: ',& x(i,j,k),y(i,j,k),z(i,j,k) call CCTK_WARN(1,warnline) write(warnline,'(a20,g16.7)') 'radius: ',r(i,j,k) call CCTK_WARN(1,warnline) write(warnline,'(a20,3g16.7)') 'velocities: ',& velxplus(i,j,k),velyplus(i,j,k),velzplus(i,j,k) call CCTK_WARN(1,warnline) call CCTK_WARN(GRHydro_c2p_warnlevel, "Specific internal energy negative") write(warnline,'(a25,4g15.6)') 'coordinates: x,y,z,r:',& x(i,j,k),y(i,j,k),z(i,j,k),r(i,j,k) call CCTK_WARN(1,warnline) !$OMP END CRITICAL endif endif end do end do end do end subroutine Conservative2PrimitiveBounds /*@@ @routine Con2PrimPolytype @date Tue Mar 19 11:43:06 2002 @author Ian Hawke @desc All routines below are identical to those above, just specialised from polytropic type EOS. @enddesc @calls @calledby @history @endhistory @@*/ subroutine Conservative2PrimitivePolytype(CCTK_ARGUMENTS) implicit none DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS integer :: i, j, k, itracer, nx, ny, nz CCTK_REAL :: uxx, uxy, uxz, uyy, uyz, uzz, det, psi4pt CCTK_INT :: type_bits, atmosphere CCTK_INT :: type2_bits CCTK_REAL :: local_min_tracer ! character(len=400) :: warnline call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere") call SpaceMask_GetStateBits(atmosphere, "Hydro_Atmosphere", "in_atmosphere") type2_bits = -1 nx = cctk_lsh(1) ny = cctk_lsh(2) nz = cctk_lsh(3) if (use_min_tracer .ne. 0) then local_min_tracer = min_tracer else local_min_tracer = 0d0 end if !!$ do k = GRHydro_stencil + 1, nz - GRHydro_stencil !!$ do j = GRHydro_stencil + 1, ny - GRHydro_stencil !!$ do i = GRHydro_stencil + 1, nx - GRHydro_stencil !$OMP PARALLEL DO PRIVATE(i,j,itracer,& !$OMP uxx, uxy, uxz, uyy, uyz, uzz, det, psi4pt) do k = 1, nz do j = 1, ny do i = 1, nx !do not compute if in atmosphere or in an excised region if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere) .or. & GRHydro_enable_internal_excision /= 0 .and. (hydro_excision_mask(i,j,k) .ne. 0)) cycle if (conformal_state .eq. 0) then call SpatialDeterminant(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) call UpperMetric(uxx,uxy,uxz,uyy,uyz,uzz,det,& gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),& gyz(i,j,k),gzz(i,j,k)) else psi4pt = psi(i,j,k)**4 call SpatialDeterminant(psi4pt*gxx(i,j,k),psi4pt*gxy(i,j,k),& psi4pt*gxz(i,j,k),psi4pt*gyy(i,j,k),& psi4pt*gyz(i,j,k),psi4pt*gzz(i,j,k),det) call UpperMetric(uxx,uxy,uxz,uyy,uyz,uzz,det,& psi4pt*gxx(i,j,k),psi4pt*gxy(i,j,k),psi4pt*gxz(i,j,k),& psi4pt*gyy(i,j,k),psi4pt*gyz(i,j,k),psi4pt*gzz(i,j,k)) end if if (evolve_tracer .ne. 0) then do itracer=1,number_of_tracers call Con2Prim_ptTracer(cons_tracer(i,j,k,itracer), & tracer(i,j,k,itracer), dens(i,j,k)) enddo if (use_min_tracer .ne. 0) then if (tracer(i,j,k,itracer) .le. local_min_tracer) then tracer(i,j,k,itracer) = local_min_tracer end if end if endif call Con2Prim_ptPolytype(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),vel(i,j,k,1),vel(i,j,k,2), & vel(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), & z(i,j,k),r(i,j,k),GRHydro_rho_min, GRHydro_reflevel, GRHydro_C2P_failed(i,j,k)) end do end do end do !$OMP END PARALLEL DO return end subroutine Conservative2PrimitivePolytype /*@@ @routine Con2Prim_ptPolytype @date Sat Jan 26 01:09:39 2002 @author The GRHydro Developers @desc Given all the appropriate data, recover the primitive variables at a single point. @enddesc @calls @calledby @history @endhistory @@*/ 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) 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, & 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, epsold, epsnew, & enthalpy, denthalpy, sqrtdet, invsqrtdet, invfac, GRHydro_C2P_failed character(len=200) warnline #ifdef _EOS_BASE_INC_ #undef _EOS_BASE_INC_ #endif #include "EOS_Base.inc" !!$ Undensitize the variables sqrtdet = sqrt(det) invsqrtdet = 1.d0/sqrtdet udens = dens*invsqrtdet usx = sx*invsqrtdet usy = sy*invsqrtdet usz = sz*invsqrtdet 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 !!$ Set initial guess for rho: rhoold = max(GRHydro_rho_min,rho) enthalpy = 1.d0 + EOS_SpecificIntEnergy(handle, rhoold, press) + & EOS_Pressure(handle, rhoold, epsilon) / rhoold w_lorentz = sqrt(1.d0 + s2 / ( (udens*enthalpy)**2 )) !!$ BEGIN: Check for NaN value (1st check) if (w_lorentz .ne. w_lorentz) then !$OMP CRITICAL write(warnline,'(a70,8g15.6)') 'NaN produced in sqrt(): (dens, det, s2, udens, enthalpy, x, y, z)', dens, det, s2, udens, enthalpy, x, y, z call CCTK_WARN(GRHydro_NaN_verbose, warnline) !$OMP END CRITICAL endif !!$ END: Check for NaN value (1st check) press = EOS_Pressure(handle, rhoold, epsilon) !!$ Calculate the function f = rhoold * w_lorentz - udens !!$Find the root count = 0 rhonew = rhoold do while ( ((abs(rhonew - rhoold)/abs(rhonew) .gt. GRHydro_perc_ptol) .and. & (abs(rhonew - rhoold) .gt. GRHydro_del_ptol)) .or. & (count .lt. GRHydro_countmin)) count = count + 1 if (count > GRHydro_countmax) then GRHydro_C2P_failed = 1 !$OMP CRITICAL call CCTK_WARN(1, 'count > GRHydro_countmax! ') write(warnline,'(a28,i2)') 'on carpet reflevel: ',GRHydro_reflevel call CCTK_WARN(1,warnline) write(warnline,'(a20,3g16.7)') 'xyz location: ',x,y,z call CCTK_WARN(1,warnline) write(warnline,'(a20,g16.7)') 'radius: ',r call CCTK_WARN(1,warnline) call CCTK_WARN(1,"Setting the point to atmosphere") !$OMP END CRITICAL ! for safety, let's set the point to atmosphere rhonew = GRHydro_rho_min udens = rhonew dens = sqrtdet * rhonew press = EOS_Pressure(handle, rhonew, 1.d0) epsilon = EOS_SpecificIntEnergy(handle, rhonew, press) sx = 0.d0 sy = 0.d0 sz = 0.d0 s2 = 0.d0 usx = 0.d0 usy = 0.d0 usz = 0.d0 w_lorentz = 1.d0 goto 50 end if !!$ Ian has the feeling that this is an error in general. But for the !!$ 2D_Polytrope it gives the right answers. denthalpy = EOS_DPressByDrho(handle, rhonew, epsilon) / rhonew df = w_lorentz - rhonew * s2 * denthalpy / & (w_lorentz*(udens**2)*(enthalpy**3)) rhoold = rhonew rhonew = rhoold - f/df if (rhonew .lt. 0.d0) then #if 0 GRHydro_C2P_failed = 1 !$OMP CRITICAL call CCTK_WARN(GRHydro_NaN_verbose, 'rhonew went below 0') !$OMP END CRITICAL #else rhonew = GRHydro_rho_min/100; #endif end if !!$ Recalculate primitive variables and function enthalpy = 1.d0 + EOS_SpecificIntEnergy(handle, rhonew, press) + & EOS_Pressure(handle, rhonew, epsilon) / rhonew w_lorentz = sqrt(1.d0 + s2 / ( (udens*enthalpy)**2 )) press = EOS_Pressure(handle, rhonew, epsilon) tau = sqrtdet * ((rho * enthalpy) * w_lorentz**2 - press) - dens f = rhonew * w_lorentz - udens enddo !!$ Calculate primitive variables from root if (rhonew .le. GRHydro_rho_min*(1.d0+GRHydro_atmo_tolerance) ) then rhonew = GRHydro_rho_min udens = rhonew ! 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 press = EOS_Pressure(handle, rhonew, 1.d0) epsilon = EOS_SpecificIntEnergy(handle, rhonew, press) ! w_lorentz=1, so the expression for utau reduces to: tau = rhonew + rhonew*epsilon - udens sx = 0.d0 sy = 0.d0 sz = 0.d0 s2 = 0.d0 usx = 0.d0 usy = 0.d0 usz = 0.d0 w_lorentz = 1.d0 end if 50 rho = rhonew enthalpy = 1.d0 + EOS_SpecificIntEnergy(handle, rhonew, press) + & EOS_Pressure(handle, rhonew, epsilon) / rhonew w_lorentz = sqrt(1.d0 + s2 / ( (udens*enthalpy)**2 )) !!$ BEGIN: Check for NaN value (2nd check) if (w_lorentz .ne. w_lorentz) then !$OMP CRITICAL write(warnline,'(a70,6g15.6)') 'NaN produced in sqrt(): (s2, udens, enthalpy, x, y, z)', s2, udens, enthalpy, x, y, z call CCTK_WARN(GRHydro_NaN_verbose, warnline) !$OMP END CRITICAL endif !!$ END: Check for NaN value (2nd check) press = EOS_Pressure(handle, rhonew, epsilon) epsilon = EOS_SpecificIntEnergy(handle, rhonew, press) tau = sqrtdet * ((rho * enthalpy) * w_lorentz**2 - press) - dens if (epsilon .lt. 0.0d0) then GRHydro_C2P_failed = 1 !$OMP CRITICAL call CCTK_WARN(1, 'epsilon < 0! ') write(warnline,'(a28,i2)') 'on carpet reflevel: ',GRHydro_reflevel call CCTK_WARN(1,warnline) write(warnline,'(a20,3g16.7)') 'xyz location: ',x,y,z call CCTK_WARN(1,warnline) write(warnline,'(a20,g16.7)') 'radius: ',r call CCTK_WARN(1,warnline) call CCTK_WARN(1,"Setting the point to atmosphere") !$OMP END CRITICAL rho = GRHydro_rho_min dens = sqrtdet * rho press = EOS_Pressure(handle, rhonew, 1.d0) epsilon = EOS_SpecificIntEnergy(handle, rhonew, press) ! w_lorentz=1, so the expression for tau reduces to: tau = sqrtdet * (rho+rho*epsilon) - dens usx = 0.d0 usy = 0.d0 usz = 0.d0 w_lorentz = 1.d0 end if invfac = 1.d0 / ( (rho + rho*epsilon + press) * w_lorentz**2) vlowx = usx * invfac vlowy = usy * invfac vlowz = usz * invfac velx = uxx * vlowx + uxy * vlowy + uxz * vlowz vely = uxy * vlowx + uyy * vlowy + uyz * vlowz velz = uxz * vlowx + uyz * vlowy + uzz * vlowz return end subroutine Con2Prim_ptPolytype /*@@ @routine Cons2PrimBoundsPolytype @date Tue Mar 12 18:04:40 2002 @author The GRHydro Developers @desc This routine is used only if the reconstruction is performed on the conserved variables. It computes the primitive variables on cell boundaries. Since reconstruction on conservative had not proved to be very successful, some of the improvements to the C2P routines (e.g. the check about whether a failure happens in a point that will be restriced anyway) are not implemented here yet. @enddesc @calls @calledby @history @endhistory @@*/ subroutine Con2PrimBoundsPolytype(CCTK_ARGUMENTS) implicit none DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS 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,psi4r,psi4l CCTK_INT :: type_bits, atmosphere CCTK_INT :: type2_bits call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere") call SpaceMask_GetStateBits(atmosphere, "Hydro_Atmosphere", "in_atmosphere") type2_bits = -1 nx = cctk_lsh(1) ny = cctk_lsh(2) nz = cctk_lsh(3) do k = GRHydro_stencil, nz - GRHydro_stencil + 1 do j = GRHydro_stencil, ny - GRHydro_stencil + 1 do i = GRHydro_stencil, nx - GRHydro_stencil + 1 !do not compute if in atmosphere or in an excised region if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere) .or. & GRHydro_enable_internal_excision /= 0 .and. (hydro_excision_mask(i,j,k) .ne. 0)) cycle gxxl = 0.5d0 * (gxx(i,j,k) + gxx(i-xoffset,j-yoffset,k-zoffset)) gxyl = 0.5d0 * (gxy(i,j,k) + gxy(i-xoffset,j-yoffset,k-zoffset)) gxzl = 0.5d0 * (gxz(i,j,k) + gxz(i-xoffset,j-yoffset,k-zoffset)) gyyl = 0.5d0 * (gyy(i,j,k) + gyy(i-xoffset,j-yoffset,k-zoffset)) gyzl = 0.5d0 * (gyz(i,j,k) + gyz(i-xoffset,j-yoffset,k-zoffset)) gzzl = 0.5d0 * (gzz(i,j,k) + gzz(i-xoffset,j-yoffset,k-zoffset)) gxxr = 0.5d0 * (gxx(i,j,k) + gxx(i+xoffset,j+yoffset,k+zoffset)) gxyr = 0.5d0 * (gxy(i,j,k) + gxy(i+xoffset,j+yoffset,k+zoffset)) gxzr = 0.5d0 * (gxz(i,j,k) + gxz(i+xoffset,j+yoffset,k+zoffset)) gyyr = 0.5d0 * (gyy(i,j,k) + gyy(i+xoffset,j+yoffset,k+zoffset)) 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)) if (conformal_state .eq. 0) then call SpatialDeterminant(gxxl,gxyl,gxzl,gyyl, gyzl,gzzl,avg_detl) call SpatialDeterminant(gxxr,gxyr,gxzr,gyyr, gyzr,gzzr,avg_detr) call UpperMetric(uxxl,uxyl,uxzl,uyyl,uyzl,uzzl,avg_detl,& gxxl,gxyl,gxzl,gyyl,gyzl,gzzl) call UpperMetric(uxxr,uxyr,uxzr,uyyr,uyzr,uzzr,avg_detr,& gxxr,gxyr,gxzr,gyyr,gyzr,gzzr) else psi4r = (0.5d0*(psi(i,j,k)+psi(i+xoffset,j+yoffset,k+zoffset)))**4 psi4l = (0.5d0*(psi(i,j,k)+psi(i-xoffset,j-yoffset,k-zoffset)))**4 call SpatialDeterminant(psi4l*gxxl,psi4l*gxyl,psi4l*gxzl,& psi4l*gyyl,psi4l*gyzl,psi4l*gzzl,avg_detl) call SpatialDeterminant(psi4r*gxxr,psi4r*gxyr,psi4r*gxzr,& psi4r*gyyr,psi4r*gyzr,psi4r*gzzr,avg_detr) call UpperMetric(uxxl,uxyl,uxzl,uyyl,uyzl,uzzl,avg_detl,& psi4l*gxxl,psi4l*gxyl,psi4l*gxzl,psi4l*gyyl,psi4l*gyzl,& psi4l*gzzl) call UpperMetric(uxxr,uxyr,uxzr,uyyr,uyzr,uzzr,avg_detr,& psi4r*gxxr,psi4r*gxyr,psi4r*gxzr,psi4r*gyyr,psi4r*gyzr,& psi4r*gzzr) end if call Con2Prim_ptPolytype(GRHydro_eos_handle, 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),& 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)) 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),& velyplus(i,j,k),velzplus(i,j,k),epsplus(i,j,k),& pressplus(i,j,k),w_lorentzplus(i,j,k),& uxxr,uxyr,uxzr,uyyr,uyzr,uzzr,avg_detr,& 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)) end do end do end do end subroutine Con2PrimBoundsPolytype /*@@ @routine Con2PrimBoundsTracer @date Mon Mar 8 13:41:55 2004 @author Ian Hawke @desc @enddesc @calls @calledby @history @endhistory @@*/ subroutine Con2PrimBoundsTracer(CCTK_ARGUMENTS) implicit none DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS integer :: i, j, k, itracer, 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,psi4r,psi4l nx = cctk_lsh(1) ny = cctk_lsh(2) nz = cctk_lsh(3) do k = GRHydro_stencil, nz - GRHydro_stencil + 1 do j = GRHydro_stencil, ny - GRHydro_stencil + 1 do i = GRHydro_stencil, nx - GRHydro_stencil + 1 gxxl = 0.5d0 * (gxx(i,j,k) + gxx(i-xoffset,j-yoffset,k-zoffset)) gxyl = 0.5d0 * (gxy(i,j,k) + gxy(i-xoffset,j-yoffset,k-zoffset)) gxzl = 0.5d0 * (gxz(i,j,k) + gxz(i-xoffset,j-yoffset,k-zoffset)) gyyl = 0.5d0 * (gyy(i,j,k) + gyy(i-xoffset,j-yoffset,k-zoffset)) gyzl = 0.5d0 * (gyz(i,j,k) + gyz(i-xoffset,j-yoffset,k-zoffset)) gzzl = 0.5d0 * (gzz(i,j,k) + gzz(i-xoffset,j-yoffset,k-zoffset)) gxxr = 0.5d0 * (gxx(i,j,k) + gxx(i+xoffset,j+yoffset,k+zoffset)) gxyr = 0.5d0 * (gxy(i,j,k) + gxy(i+xoffset,j+yoffset,k+zoffset)) gxzr = 0.5d0 * (gxz(i,j,k) + gxz(i+xoffset,j+yoffset,k+zoffset)) gyyr = 0.5d0 * (gyy(i,j,k) + gyy(i+xoffset,j+yoffset,k+zoffset)) 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)) if (conformal_state .eq. 0) then call SpatialDeterminant(gxxl,gxyl,gxzl,gyyl, gyzl,gzzl,avg_detl) call SpatialDeterminant(gxxr,gxyr,gxzr,gyyr, gyzr,gzzr,avg_detr) call UpperMetric(uxxl,uxyl,uxzl,uyyl,uyzl,uzzl,avg_detl,& gxxl,gxyl,gxzl,gyyl,gyzl,gzzl) call UpperMetric(uxxr,uxyr,uxzr,uyyr,uyzr,uzzr,avg_detr,& gxxr,gxyr,gxzr,gyyr,gyzr,gzzr) else psi4r = (0.5d0*(psi(i,j,k)+psi(i+xoffset,j+yoffset,k+zoffset)))**4 psi4l = (0.5d0*(psi(i,j,k)+psi(i-xoffset,j-yoffset,k-zoffset)))**4 call SpatialDeterminant(psi4l*gxxl,psi4l*gxyl,psi4l*gxzl,& psi4l*gyyl,psi4l*gyzl,psi4l*gzzl,avg_detl) call SpatialDeterminant(psi4r*gxxr,psi4r*gxyr,psi4r*gxzr,& psi4r*gyyr,psi4r*gyzr,psi4r*gzzr,avg_detr) call UpperMetric(uxxl,uxyl,uxzl,uyyl,uyzl,uzzl,avg_detl,& psi4l*gxxl,psi4l*gxyl,psi4l*gxzl,psi4l*gyyl,psi4l*gyzl,& psi4l*gzzl) call UpperMetric(uxxr,uxyr,uxzr,uyyr,uyzr,uzzr,avg_detr,& psi4r*gxxr,psi4r*gxyr,psi4r*gxzr,psi4r*gyyr,psi4r*gyzr,& psi4r*gzzr) end if do itracer=1,number_of_tracers call Con2Prim_ptBoundsTracer(cons_tracerplus(i,j,k,itracer), & tracerplus(i,j,k,itracer), & rhoplus(i,j,k), sqrt(1.d0 - & (gxxr * velxplus(i,j,k)**2 + & gyyr * velyplus(i,j,k)**2 + & gzzr * velzplus(i,j,k)**2 + & 2.d0 * (gxyr * velxplus(i,j,k) * velyplus(i,j,k) + & gxzr * velxplus(i,j,k) * velzplus(i,j,k) + & gyzr * velyplus(i,j,k) * velzplus(i,j,k) ) ) ), & avg_detr) call Con2Prim_ptBoundsTracer(cons_tracerminus(i,j,k,itracer), & tracerminus(i,j,k,itracer), & rhominus(i,j,k), sqrt(1.d0 - & (gxxr * velxminus(i,j,k)**2 + & gyyr * velyminus(i,j,k)**2 + & gzzr * velzminus(i,j,k)**2 + & 2.d0 * (gxyr * velxminus(i,j,k) * velyminus(i,j,k) + & gxzr * velxminus(i,j,k) * velzminus(i,j,k) + & gyzr * velyminus(i,j,k) * velzminus(i,j,k) ) ) ), & avg_detl) enddo !!$ tracerplus(i,j,k) = cons_tracerplus(i,j,k) / & !!$ sqrt(avg_detr) / rhoplus(i,j,k) * & !!$ sqrt(1.d0 - & !!$ (gxxr * velxplus(i,j,k)**2 + & !!$ gyyr * velyplus(i,j,k)**2 + & !!$ gzzr * velzplus(i,j,k)**2 + & !!$ 2.d0 * (gxyr * velxplus(i,j,k) * velyplus(i,j,k) + & !!$ gxzr * velxplus(i,j,k) * velzplus(i,j,k) + & !!$ gyzr * velyplus(i,j,k) * velzplus(i,j,k) ) ) ) !!$ tracerminus(i,j,k) = cons_tracerminus(i,j,k) / & !!$ sqrt(avg_detl) / rhominus(i,j,k) * & !!$ sqrt(1.d0 - & !!$ (gxxl * velxminus(i,j,k)**2 + & !!$ gyyl * velyminus(i,j,k)**2 + & !!$ gzzl * velzminus(i,j,k)**2 + & !!$ 2.d0 * (gxyl * velxminus(i,j,k) * velyminus(i,j,k) + & !!$ gxzl * velxminus(i,j,k) * velzminus(i,j,k) + & !!$ gyzl * velyminus(i,j,k) * velzminus(i,j,k) ) ) ) end do end do end do end subroutine Con2PrimBoundsTracer /*@@ @routine Con2Prim_ptTracer @date Mon Mar 8 14:26:29 2004 @author Ian Hawke @desc @enddesc @calls @calledby @history @endhistory @@*/ subroutine Con2Prim_ptTracer(cons_tracer, tracer, dens) implicit none DECLARE_CCTK_PARAMETERS CCTK_REAL cons_tracer, tracer, dens tracer = cons_tracer / dens end subroutine Con2Prim_ptTracer /*@@ @routine Con2Prim_ptTracerBounds @date Mon Mar 8 14:26:29 2004 @author Ian Hawke @desc @enddesc @calls @calledby @history @endhistory @@*/ subroutine Con2Prim_ptBoundsTracer(cons_tracer, tracer, rho, one_over_w_lorentz, det) implicit none DECLARE_CCTK_PARAMETERS CCTK_REAL cons_tracer, tracer, rho, one_over_w_lorentz, det tracer = cons_tracer / sqrt(det) / rho * one_over_w_lorentz end subroutine Con2Prim_ptBoundsTracer /*@@ @routine Conservative2PrimitiveGeneral @date Sat Jan 26 01:08:33 2002 @author Ian Hawke @desc Wrapper routine that converts from conserved to primitive variables at every grid cell centre. Converted to the EOSGeneral format @enddesc @calls @calledby @history @endhistory @@*/ subroutine Conservative2PrimitiveGeneral(CCTK_ARGUMENTS) USE GRHydro_Scalars implicit none DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS CCTK_INT :: i, j, k, itracer, nx, ny, nz, ierr CCTK_REAL :: psi4pt character(len=100) warnline integer :: maxerrloc(3) integer :: ii,jj,kk CCTK_REAL :: maxerr CCTK_INT :: count CCTK_REAL :: s2, df, v2, w, vlowx, vlowy, vlowz CCTK_REAL :: temp1, temp2, temp3 CCTK_REAL :: udens, usx, usy, usz, utau CCTK_REAL :: drhodpress, depsdpress logical :: loop_condition CCTK_REAL, dimension(:,:,:), allocatable :: f logical, dimension(:,:,:), allocatable :: point_loop_condition CCTK_INT :: type_bits, atmosphere CCTK_INT :: type2_bits CCTK_REAL :: local_min_tracer CCTK_INT, dimension(3) :: loc, loc2 allocate(f(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),& point_loop_condition(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),& STAT=ierr) if (ierr .ne. 0) then call CCTK_WARN(0, "Allocation problems with f") end if call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere") call SpaceMask_GetStateBits(atmosphere, "Hydro_Atmosphere", "in_atmosphere") type2_bits = -1 nx = cctk_lsh(1) ny = cctk_lsh(2) nz = cctk_lsh(3) count = 0 press_old = max(eosgeneral_pmin, press) if (use_min_tracer .ne. 0) then local_min_tracer = min_tracer else local_min_tracer = 0d0 end if loop_condition = .true. !!$ Set up rho and epsilon do k = 1, nz do j = 1, ny do i = 1, nx if (evolve_tracer .ne. 0) then do itracer=1,number_of_tracers call Con2Prim_ptTracer(cons_tracer(i,j,k,itracer), & tracer(i,j,k,itracer), dens(i,j,k)) if (use_min_tracer .ne. 0) then if (tracer(i,j,k,itracer) .le. local_min_tracer) then tracer(i,j,k,itracer) = local_min_tracer end if end if enddo endif if ( (dens(i,j,k).le.sqrt(GRHydro_Det(i,j,k)) * & GRHydro_rho_min*(1.0d0+GRHydro_atmo_tolerance)) & .or.(tau(i,j,k) .le. 0d0)) then rho(i,j,k) = GRHydro_rho_min dens(i,j,k) = sqrt(GRHydro_det(i,j,k)) * rho(i,j,k) scon(i,j,k,:) = 0.d0 w_lorentz(i,j,k) = 1.d0 press(i,j,k) = 0.1d0 * eosgeneral_pmin ! eps(i,j,k) = press(i,j,k) / rho(i,j,k) ! Note that this should be improved eps(i,j,k) = 1.0e-12 tau(i,j,k) = sqrt(GRHydro_det(i,j,k)) * rho(i,j,k) * eps(i,j,k) SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bits, atmosphere) atmosphere_mask(i,j,k) = 1 end if udens = dens(i,j,k) / sqrt(GRHydro_det(i,j,k)) usx = scon(i,j,k,1) / sqrt(GRHydro_det(i,j,k)) usy = scon(i,j,k,2) / sqrt(GRHydro_det(i,j,k)) usz = scon(i,j,k,3) / sqrt(GRHydro_det(i,j,k)) utau = tau(i,j,k) / sqrt(GRHydro_det(i,j,k)) s2 = usx*usx*GRHydro_uxx(i,j,k) + & usy*usy*GRHydro_uyy(i,j,k) + & usz*usz*GRHydro_uzz(i,j,k) + & 2.d0*usx*usy*GRHydro_uxy(i,j,k) + & 2.d0*usx*usz*GRHydro_uxz(i,j,k) + & 2.d0*usy*usz*GRHydro_uyz(i,j,k) !!$ Check that the variables have a chance of being physical if( (utau + press_old(i,j,k) + udens)**2 - s2 .le. 0.0d0) then if(GRHydro_reflevel.ge.GRHydro_c2p_warn_from_reflevel) then !$OMP CRITICAL call CCTK_WARN(1,'Con2PrimGeneral: variables unphysical!') write(warnline,'(a28,i2)') 'on carpet reflevel: ',GRHydro_reflevel call CCTK_WARN(1,warnline) write(warnline,'(a25,g15.6)') 'undensitized energy: ',utau call CCTK_WARN(1,warnline) write(warnline,'(a25,g15.6)') 'undensitized density: ',udens call CCTK_WARN(1,warnline) write(warnline,'(a25,g15.6)') 'undensitized s2: ',s2 call CCTK_WARN(1,warnline) write(warnline,'(a25,g15.6)') 'pressure guess: ',press_old(i,j,k) call CCTK_WARN(1,warnline) write(warnline,'(a25,i4)') 'refinement level: ',GRHydro_reflevel call CCTK_WARN(1,warnline) write(warnline,'(a25,4g15.6)') 'coordinates: x,y,z,r:',& x(i,j,k),y(i,j,k),z(i,j,k),r(i,j,k) call CCTK_WARN(1,warnline) call CCTK_WARN(GRHydro_c2p_warnlevel, "Unphysical variables") !$OMP END CRITICAL exit else !$OMP CRITICAL write(warnline,'(a60,i2)') 'Con2Prim: eps negative, but I was told to ignore level: ',GRHydro_reflevel call CCTK_WARN(1,warnline) write(warnline,'(a25,4g15.6)') 'coordinates: x,y,z,r:',& x(i,j,k),y(i,j,k),z(i,j,k),r(i,j,k) call CCTK_WARN(1,warnline) !$OMP END CRITICAL exit endif endif temp1 = (utau + udens + press_old(i,j,k)) temp2 = temp1**2 - s2 temp3 = sqrt(temp2) rho(i,j,k) = udens * temp3 / temp1 w_lorentz(i,j,k) = temp1 / temp3 !!$ Eps not previously set? eps(i,j,k) = & (temp3 - press_old(i,j,k) * w_lorentz(i,j,k) - udens) / & udens end do end do end do ierr = EOS_SetGFs(cctkGH, EOS_Con2PrimCall) where (press .le. eosgeneral_pmin) press = eosgeneral_pmin rho = GRHydro_rho_min dens = sqrt(GRHydro_det) * rho scon(:,:,:,1) = 0.d0 scon(:,:,:,2) = 0.d0 scon(:,:,:,3) = 0.d0 w_lorentz = 1d0 ! eps = press / rho eps = 1.0e-12 tau = sqrt(GRHydro_det) * rho * eps end where f = press_old - press press_new = press_old do while (loop_condition) count = count + 1 if (count > GRHydro_countmax) then if(GRHydro_reflevel.ge.GRHydro_c2p_warn_from_reflevel) then !$OMP CRITICAL call CCTK_WARN(1, 'Con2PrimGeneral: error: did not converge in ') write(warnline,'(a20,i12,a10)') ' ',& GRHydro_countmax,' steps' call CCTK_WARN(1,warnline) write(warnline,'(a28,i2)') 'on carpet reflevel: ',GRHydro_reflevel call CCTK_WARN(1,warnline) maxerr = maxval(abs(press_new - press_old) & / abs(press_new)) maxerrloc = maxloc(abs(press_new - press_old) & / abs(press_new)) ii = maxerrloc(1) jj = maxerrloc(2) kk = maxerrloc(3) write(warnline,'(1P10E15.6)') maxerr call CCTK_WARN(1,warnline) write(warnline,'(1P10E15.6)') press_new(ii,jj,kk), press_old(ii,jj,kk) call CCTK_WARN(1,warnline) write(warnline,'(3i4)') ii,jj,kk call CCTK_WARN(1,warnline) write(warnline,'(1P10E15.6)') x(ii,jj,kk), & y(ii,jj,kk), & z(ii,jj,kk) call CCTK_WARN(1,warnline) call CCTK_WARN(0,"Abort!") !$OMP END CRITICAL exit else !$OMP CRITICAL write(warnline,'(a60,i2)') 'Con2Prim: eps negative, but I was told to ignore level: ',& GRHydro_reflevel call CCTK_WARN(1,warnline) !$OMP END CRITICAL exit endif end if do k = 1, nz do j = 1, ny do i = 1, nx !do not compute if in atmosphere or in an excised region if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere) .or. & GRHydro_enable_internal_excision /= 0 .and. (hydro_excision_mask(i,j,k) .ne. 0)) & then press_new(i,j,k) = press(i,j,k) press_old(i,j,k) = press(i,j,k) ! cycle endif ! All right, is this point already converged?!? loop_condition = abs(press_new(i,j,k)-press_old(i,j,k)) / abs(press_new(i,j,k)) > & GRHydro_perc_ptol loop_condition = loop_condition .and. & abs(press_new(i,j,k)-press_old(i,j,k)) > GRHydro_del_ptol loop_condition = loop_condition .or. (count < GRHydro_countmin) & .or. (count < 2) ! we have to at least go once through c2p. point_loop_condition(i,j,k) = loop_condition if(.not.loop_condition) cycle udens = dens(i,j,k) / sqrt(GRHydro_det(i,j,k)) usx = scon(i,j,k,1) / sqrt(GRHydro_det(i,j,k)) usy = scon(i,j,k,2) / sqrt(GRHydro_det(i,j,k)) usz = scon(i,j,k,3) / sqrt(GRHydro_det(i,j,k)) utau = tau(i,j,k) / sqrt(GRHydro_det(i,j,k)) s2 = usx*usx*GRHydro_uxx(i,j,k) + & usy*usy*GRHydro_uyy(i,j,k) + & usz*usz*GRHydro_uzz(i,j,k) + & 2.d0*usx*usy*GRHydro_uxy(i,j,k) + & 2.d0*usx*usz*GRHydro_uxz(i,j,k) + & 2.d0*usy*usz*GRHydro_uyz(i,j,k) temp1 = (utau + udens + press_new(i,j,k)) temp2 = temp1**2 - s2 temp3 = sqrt(temp2) drhodpress = udens * s2 / (temp3 * temp1**2) depsdpress = press_new(i,j,k) * s2 / (rho(i,j,k) * temp1 * temp2) df = 1.d0 - eos_dpdeps_temp(i,j,k) * depsdpress - & eos_dpdrho_temp(i,j,k) * drhodpress press_old(i,j,k) = press_new(i,j,k) press_new(i,j,k) = press_old(i,j,k) - f(i,j,k) / df if (press_new(i,j,k) .le. eosgeneral_pmin) then press_new(i,j,k) = eosgeneral_pmin rho(i,j,k) = GRHydro_rho_min dens(i,j,k) = sqrt(GRHydro_det(i,j,k)) * rho(i,j,k) scon(i,j,k,:) = 0.d0 w_lorentz(i,j,k) = 1d0 eps(i,j,k) = press(i,j,k) / rho(i,j,k) tau(i,j,k) = sqrt(GRHydro_det(i,j,k)) * rho(i,j,k) * eps(i,j,k) udens = rho(i,j,k) utau = rho(i,j,k) * eps(i,j,k) s2 = 0.d0 end if temp1 = (utau + udens + press_new(i,j,k)) temp2 = temp1**2 - s2 temp3 = sqrt(temp2) rho(i,j,k) = udens * temp3 / temp1 w_lorentz(i,j,k) = temp1 / temp3 eps(i,j,k) = & (temp3 - press_new(i,j,k) * w_lorentz(i,j,k) - udens) / & udens point_loop_condition(i,j,k) = & (abs(press_new(i,j,k) - press_old(i,j,k)) / abs(press_new(i,j,k)) > & GRHydro_perc_ptol) point_loop_condition(i,j,k) = point_loop_condition(i,j,k) .and. & (abs(press_new(i,j,k) - press_old(i,j,k)) > & GRHydro_del_ptol) end do end do end do ierr = EOS_SetGFs(cctkGH, EOS_Con2PrimCall) press = max(press, eosgeneral_pmin) f = press_new - press loop_condition = any(point_loop_condition) loop_condition = loop_condition .or. & (count < GRHydro_countmin) end do deallocate(point_loop_condition) do count = 1, GRHydro_polish do k = 1, nz do j = 1, ny do i = 1, nx !do not compute if in atmosphere or in an excised region if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere) .or. & GRHydro_enable_internal_excision /= 0 .and. (hydro_excision_mask(i,j,k) .ne. 0)) cycle f = press_new(i,j,k) - press(i,j,k) udens = dens(i,j,k) / sqrt(GRHydro_det(i,j,k)) usx = scon(i,j,k,1) / sqrt(GRHydro_det(i,j,k)) usy = scon(i,j,k,2) / sqrt(GRHydro_det(i,j,k)) usz = scon(i,j,k,3) / sqrt(GRHydro_det(i,j,k)) utau = tau(i,j,k) / sqrt(GRHydro_det(i,j,k)) s2 = usx*usx*GRHydro_uxx(i,j,k) + & usy*usy*GRHydro_uyy(i,j,k) + & usz*usz*GRHydro_uzz(i,j,k) + & 2.d0*usx*usy*GRHydro_uxy(i,j,k) + & 2.d0*usx*usz*GRHydro_uxz(i,j,k) + & 2.d0*usy*usz*GRHydro_uyz(i,j,k) temp1 = (utau + udens + press_new(i,j,k)) temp2 = temp1**2 - s2 temp3 = sqrt(temp2) drhodpress = udens * s2 / (temp3 * temp1**2) depsdpress = press_new(i,j,k) * s2 / (rho(i,j,k) * temp1 * temp2) df = 1.d0 - eos_dpdeps_temp(i,j,k) * depsdpress - & eos_dpdrho_temp(i,j,k) * drhodpress press_new(i,j,k) = max(press_new(i,j,k) - f(i,j,k) / df, eosgeneral_pmin) temp1 = (utau + udens + press_new(i,j,k)) temp2 = temp1**2 - s2 temp3 = sqrt(temp2) rho(i,j,k) = udens * temp3 / temp1 w_lorentz(i,j,k) = temp1 / temp3 eps(i,j,k) = & (temp3 - press_new(i,j,k) * w_lorentz(i,j,k) - udens) / & udens end do end do end do ierr = EOS_SetGFs(cctkGH, EOS_Con2PrimCall) press = max(press, eosgeneral_pmin) end do do k = 1, nz do j = 1, ny do i = 1, nx !do not compute if in atmosphere or in an excised region if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere) .or. & GRHydro_enable_internal_excision /= 0 .and. (hydro_excision_mask(i,j,k) .ne. 0)) cycle !!$ Note that in the following we enforce eps > 0 !!$ This makes the later check redundant if ( (rho(i,j,k) .le. GRHydro_rho_min*(1.d0+GRHydro_atmo_tolerance) ) .or. & (eps(i,j,k) .lt. 0d0) ) then rho(i,j,k) = GRHydro_rho_min udens = rho(i,j,k) utau = tau(i,j,k) / sqrt(GRHydro_det(i,j,k)) dens(i,j,k) = sqrt(GRHydro_det(i,j,k)) * rho(i,j,k) eps(i,j,k) = utau / udens scon(i,j,k,:) = 0.d0 w_lorentz(i,j,k) = 1.d0 usx = 0.d0 usy = 0.d0 usz = 0.d0 else udens = dens(i,j,k) / sqrt(GRHydro_det(i,j,k)) usx = scon(i,j,k,1) / sqrt(GRHydro_det(i,j,k)) usy = scon(i,j,k,2) / sqrt(GRHydro_det(i,j,k)) usz = scon(i,j,k,3) / sqrt(GRHydro_det(i,j,k)) utau = tau(i,j,k) / sqrt(GRHydro_det(i,j,k)) end if vlowx = usx / & ( (rho(i,j,k) + rho(i,j,k) * eps(i,j,k) + press(i,j,k)) * & w_lorentz(i,j,k)**2) vlowy = usy / & ( (rho(i,j,k) + rho(i,j,k) * eps(i,j,k) + press(i,j,k)) * & w_lorentz(i,j,k)**2) vlowz = usz / & ( (rho(i,j,k) + rho(i,j,k) * eps(i,j,k) + press(i,j,k)) * & w_lorentz(i,j,k)**2) vel(i,j,k,1) = GRHydro_uxx(i,j,k) * vlowx + & GRHydro_uxy(i,j,k) * vlowy + & GRHydro_uxz(i,j,k) * vlowz vel(i,j,k,2) = GRHydro_uxy(i,j,k) * vlowx + & GRHydro_uyy(i,j,k) * vlowy + & GRHydro_uyz(i,j,k) * vlowz vel(i,j,k,3) = GRHydro_uxz(i,j,k) * vlowx + & GRHydro_uyz(i,j,k) * vlowy + & GRHydro_uzz(i,j,k) * vlowz !!$ If eps negative then complain vociferously if (eps(i,j,k) .le. 0.d0) then if(GRHydro_reflevel.ge.GRHydro_c2p_warn_from_reflevel) then !$OMP CRITICAL call CCTK_WARN(1,'Con2PrimGeneral: stopping the code.') call CCTK_WARN(1, ' specific internal energy just went below 0! ') write(warnline,'(a28,i2)') 'on carpet reflevel: ',GRHydro_reflevel call CCTK_WARN(1,warnline) write(warnline,'(a20,3i6)') 'ijk location: ',& i,j,k call CCTK_WARN(1,warnline) write(warnline,'(a20,1P10E15.6)') 'xyz location: ',& x(i,j,k),y(i,j,k),z(i,j,k) call CCTK_WARN(1,warnline) write(warnline,'(a25,i5)') 'refinement level: ',GRHydro_reflevel call CCTK_WARN(1,warnline) write(warnline,'(a20,1P1E15.6)') 'radius: ',r(i,j,k) call CCTK_WARN(1,warnline) write(warnline,'(a20,1P1E15.6)') 'velocities: ',& vel(i,j,k,1),vel(i,j,k,2),vel(i,j,k,3) call CCTK_WARN(1,warnline) call CCTK_WARN(GRHydro_c2p_warnlevel, \ "Specific internal energy negative") !$OMP END CRITICAL exit else !$OMP CRITICAL write(warnline,'(a60,i2)') 'Con2Prim: eps negative, but \ I was told to ignore level: ',GRHydro_reflevel call CCTK_WARN(1,warnline) write(warnline,'(a25,4g15.6)') 'coordinates: x,y,z,r:',& x(i,j,k),y(i,j,k),z(i,j,k),r(i,j,k) call CCTK_WARN(1,warnline) !$OMP END CRITICAL exit endif endif end do end do end do deallocate(f) end subroutine Conservative2PrimitiveGeneral /*@@ @routine Conservative2PrimitivePolytypeGeneral @date Sat Jan 26 01:08:33 2002 @author Ian Hawke @desc Wrapper routine that converts from conserved to primitive variables at every grid cell centre. Converted to the EOSGeneral format @enddesc @calls @calledby @history @endhistory @@*/ subroutine Con2PrimPolytypeGeneral(CCTK_ARGUMENTS) USE GRHydro_Scalars implicit none DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS CCTK_INT :: i, j, k, itracer, nx, ny, nz CCTK_REAL :: psi4pt character(len=100) warnline CCTK_INT :: count, ierr CCTK_REAL :: s2, f, df, v2, w, vlowx, vlowy, vlowz CCTK_REAL :: temp1, temp2, temp3 CCTK_REAL :: udens, usx, usy, usz CCTK_REAL :: enthalpy, denthalpy logical :: loop_condition CCTK_INT :: type_bits CCTK_INT :: atmosphere CCTK_INT :: type2_bits CCTK_REAL :: local_min_tracer call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere") call SpaceMask_GetStateBits(atmosphere, "Hydro_Atmosphere", "in_atmosphere") type2_bits = -1 nx = cctk_lsh(1) ny = cctk_lsh(2) nz = cctk_lsh(3) if (use_min_tracer .ne. 0) then local_min_tracer = min_tracer else local_min_tracer = 0d0 end if !!$ In what follows press_temp is really rho_temp count = 0 press_new = max(GRHydro_rho_min, rho) loop_condition = .true. ierr = EOS_SetGFs(cctkGH, EOS_Con2PrimCall) !!$ Set up W and epsilon do k = 1, nz do j = 1, ny do i = 1, nx if (evolve_tracer .ne. 0) then do itracer=1,number_of_tracers call Con2Prim_ptTracer(cons_tracer(i,j,k,itracer), & tracer(i,j,k,itracer), dens(i,j,k)) if (use_min_tracer .ne. 0) then if (tracer(i,j,k,itracer) .le. local_min_tracer) then tracer(i,j,k,itracer) = local_min_tracer end if end if enddo end if udens = dens(i,j,k) / sqrt(GRHydro_det(i,j,k)) usx = scon(i,j,k,1) / sqrt(GRHydro_det(i,j,k)) usy = scon(i,j,k,2) / sqrt(GRHydro_det(i,j,k)) usz = scon(i,j,k,3) / sqrt(GRHydro_det(i,j,k)) s2 = usx*usx*GRHydro_uxx(i,j,k) + & usy*usy*GRHydro_uyy(i,j,k) + & usz*usz*GRHydro_uzz(i,j,k) + & 2.d0*usx*usy*GRHydro_uxy(i,j,k) + & 2.d0*usx*usz*GRHydro_uxz(i,j,k) + & 2.d0*usy*usz*GRHydro_uyz(i,j,k) enthalpy = 1.d0 + eps(i,j,k) + press(i,j,k) / press_new(i,j,k) w_lorentz(i,j,k) = sqrt( 1.d0 + s2 / ( (udens*enthalpy)**2 ) ) end do end do end do do while (loop_condition) count = count + 1 if (count > GRHydro_countmax) then if (GRHydro_reflevel.ge.GRHydro_c2p_warn_from_reflevel) then !$OMP CRITICAL call CCTK_WARN(1, 'Con2PrimPolytypeGeneral: error: did not converge in ') write(warnline,'(a20,g12.7,a10)') ' ',& GRHydro_countmax,' steps' call CCTK_WARN(1,warnline) write(warnline,'(a28,i2)') 'on carpet reflevel: ',GRHydro_reflevel call CCTK_WARN(1,warnline) call CCTK_WARN(GRHydro_c2p_warnlevel, "Did not converge") !$OMP END CRITICAL exit else !$OMP CRITICAL write(warnline,'(a60,i2)') 'Con2Prim: eps negative, but I was told to ignore level: ',GRHydro_reflevel call CCTK_WARN(1,warnline) !$OMP END CRITICAL exit endif end if do k = 1, nz do j = 1, ny do i = 1, nx !do not compute if in atmosphere or in an excised region if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere) .or. & GRHydro_enable_internal_excision /= 0 .and. (hydro_excision_mask(i,j,k) .ne. 0)) cycle enthalpy = 1.d0 + eps(i,j,k) + press(i,j,k) / press_new(i,j,k) tau(i,j,k) = sqrt(GRHydro_det(i,j,k)) * & ( press_new(i,j,k) * enthalpy * w_lorentz(i,j,k)**2 - & press(i,j,k) ) - dens(i,j,k) udens = dens(i,j,k) / sqrt(GRHydro_det(i,j,k)) usx = scon(i,j,k,1) / sqrt(GRHydro_det(i,j,k)) usy = scon(i,j,k,2) / sqrt(GRHydro_det(i,j,k)) usz = scon(i,j,k,3) / sqrt(GRHydro_det(i,j,k)) s2 = usx*usx*GRHydro_uxx(i,j,k) + & usy*usy*GRHydro_uyy(i,j,k) + & usz*usz*GRHydro_uzz(i,j,k) + & 2.d0*usx*usy*GRHydro_uxy(i,j,k) + & 2.d0*usx*usz*GRHydro_uxz(i,j,k) + & 2.d0*usy*usz*GRHydro_uyz(i,j,k) f = press_new(i,j,k) * w_lorentz(i,j,k) - udens denthalpy = eos_dpdrho_temp(i,j,k) / press_new(i,j,k) df = w_lorentz(i,j,k) - press_new(i,j,k) * s2 * denthalpy / & ( w_lorentz(i,j,k) * (udens**2) * (enthalpy**3) ) rho(i,j,k) = press_new(i,j,k) press_new(i,j,k) = max(press_new(i,j,k) - f / df, GRHydro_rho_min) enthalpy = 1.d0 + eps(i,j,k) + press(i,j,k) / press_new(i,j,k) w_lorentz(i,j,k) = sqrt( 1.d0 + s2 / ( (udens*enthalpy)**2 ) ) end do end do end do ierr = EOS_SetGFs(cctkGH, EOS_Con2PrimCall) loop_condition = (maxval(abs(rho - press_new) / abs(press_new)) > & GRHydro_perc_ptol) loop_condition = loop_condition .and. & (maxval(abs(rho - press_new)) > GRHydro_del_ptol) loop_condition = loop_condition .or. & (count < GRHydro_countmin) end do do k = 1, nz do j = 1, ny do i = 1, nx !do not compute if in atmosphere or in an excised region if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere) .or. & GRHydro_enable_internal_excision /= 0 .and. (hydro_excision_mask(i,j,k) .ne. 0)) cycle udens = dens(i,j,k) / sqrt(GRHydro_det(i,j,k)) if (rho(i,j,k) .le. GRHydro_rho_min*(1.d0+GRHydro_atmo_tolerance) ) then rho(i,j,k) = GRHydro_rho_min scon(i,j,k,:) = 0.d0 w_lorentz(i,j,k) = 1.d0 usx = 0.d0 usy = 0.d0 usz = 0.d0 s2 = 0.d0 else usx = scon(i,j,k,1) / sqrt(GRHydro_det(i,j,k)) usy = scon(i,j,k,2) / sqrt(GRHydro_det(i,j,k)) usz = scon(i,j,k,3) / sqrt(GRHydro_det(i,j,k)) s2 = usx*usx*GRHydro_uxx(i,j,k) + & usy*usy*GRHydro_uyy(i,j,k) + & usz*usz*GRHydro_uzz(i,j,k) + & 2.d0*usx*usy*GRHydro_uxy(i,j,k) + & 2.d0*usx*usz*GRHydro_uxz(i,j,k) + & 2.d0*usy*usz*GRHydro_uyz(i,j,k) end if enthalpy = 1.d0 + eps(i,j,k) + press(i,j,k) / rho(i,j,k) w_lorentz(i,j,k) = sqrt( 1.d0 + s2 / ( (udens*enthalpy)**2 ) ) tau(i,j,k) = sqrt(GRHydro_det(i,j,k)) * & ( (rho(i,j,k) * enthalpy) * w_lorentz(i,j,k)**2 - & press(i,j,k) ) - dens(i,j,k) vlowx = usx / & ( (rho(i,j,k) + rho(i,j,k) * eps(i,j,k) + press(i,j,k)) * & w_lorentz(i,j,k)**2) vlowy = usy / & ( (rho(i,j,k) + rho(i,j,k) * eps(i,j,k) + press(i,j,k)) * & w_lorentz(i,j,k)**2) vlowz = usz / & ( (rho(i,j,k) + rho(i,j,k) * eps(i,j,k) + press(i,j,k)) * & w_lorentz(i,j,k)**2) vel(i,j,k,1) = GRHydro_uxx(i,j,k) * vlowx + & GRHydro_uxy(i,j,k) * vlowy + & GRHydro_uxz(i,j,k) * vlowz vel(i,j,k,2) = GRHydro_uxy(i,j,k) * vlowx + & GRHydro_uyy(i,j,k) * vlowy + & GRHydro_uyz(i,j,k) * vlowz vel(i,j,k,3) = GRHydro_uxz(i,j,k) * vlowx + & GRHydro_uyz(i,j,k) * vlowy + & GRHydro_uzz(i,j,k) * vlowz !!$ If eps negative then complain vociferously if (eps(i,j,k) .lt. 0.0d0) then call CCTK_WARN(1, "c2p failed: eps went negative") GRHydro_C2P_failed(i,j,k) = 1 endif end do end do end do end subroutine Con2PrimPolytypeGeneral ! subroutines to manage the C2P failure mask subroutine reset_GRHydro_C2P_failed(CCTK_ARGUMENTS) implicit none DECLARE_CCTK_ARGUMENTS GRHydro_C2P_failed = 0 return end subroutine reset_GRHydro_C2P_failed subroutine sync_GRHydro_C2P_failed(CCTK_ARGUMENTS) ! a dummy routine to syncronise GRHydro_C2P_failed implicit none DECLARE_CCTK_ARGUMENTS return end subroutine sync_GRHydro_C2P_failed subroutine check_GRHydro_C2P_failed(CCTK_ARGUMENTS) implicit none DECLARE_CCTK_ARGUMENTS integer :: i, j, k, nx, ny, nz character(len=100) warnline ! call CCTK_INFO("Checking the C2P failure mask.") nx = cctk_lsh(1) ny = cctk_lsh(2) nz = cctk_lsh(3) !$OMP PARALLEL DO PRIVATE(i,j) do k = 1, nz do j = 1, ny do i = 1, nx if (GRHydro_C2P_failed(i,j,k) == 1) then !$OMP CRITICAL call CCTK_WARN(1,'Con2Prim failed; stopping the code.') call CCTK_WARN(1,'Even with mesh refinement, this point is not restricted from a finer level, so this is really an error') write(warnline,'(a28,i2)') 'on carpet reflevel: ',GRHydro_reflevel call CCTK_WARN(1,warnline) write(warnline,'(a20,3g16.7)') 'xyz location: ',& x(i,j,k),y(i,j,k),z(i,j,k) call CCTK_WARN(1,warnline) write(warnline,'(a20,g16.7)') 'radius: ',r(i,j,k) call CCTK_WARN(1,warnline) write(warnline,'(a20,3g16.7)') 'velocities: ',& vel(i,j,k,1),vel(i,j,k,2),vel(i,j,k,3) call CCTK_WARN(0,warnline) !$OMP END CRITICAL end if end do end do end do !$OMP END PARALLEL DO return end subroutine check_GRHydro_C2P_failed