diff options
Diffstat (limited to 'src/GRHydro_Con2PrimHot.F90')
-rw-r--r-- | src/GRHydro_Con2PrimHot.F90 | 658 |
1 files changed, 658 insertions, 0 deletions
diff --git a/src/GRHydro_Con2PrimHot.F90 b/src/GRHydro_Con2PrimHot.F90 new file mode 100644 index 0000000..58cfa29 --- /dev/null +++ b/src/GRHydro_Con2PrimHot.F90 @@ -0,0 +1,658 @@ +#include "cctk.h" +#include "cctk_Parameters.h" +#include "cctk_Arguments.h" +#include "cctk_Functions.h" +#include "SpaceMask.h" +#include "GRHydro_Macros.h" + + +subroutine Conservative2PrimitiveHot(CCTK_ARGUMENTS) + + implicit none + + ! save memory when MP is not used + ! TARGET as to be before DECLARE_CCTK_ARGUMENTS for gcc 4.1 + TARGET gaa, gab, gac, gbb, gbc, gcc + TARGET gxx, gxy, gxz, gyy, gyz, gzz + TARGET lvel, vel + + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + 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 :: vlowx, vlowy, vlowz + logical :: epsnegative + character*512 :: warnline + + CCTK_REAL :: local_min_tracer + CCTK_REAL :: local_perc_ptol + + ! save memory when MP is not used + CCTK_INT :: GRHydro_UseGeneralCoordinates + CCTK_REAL, DIMENSION(:,:,:), POINTER :: g11, g12, g13, g22, g23, g33 + CCTK_REAL, DIMENSION(:,:,:,:), POINTER :: vup + +! begin EOS Omni vars + integer :: n,keytemp,anyerr,keyerr(1) + real*8 :: xpress(1),xeps(1),xtemp(1),xye(1),xrho(1) + n = 1;keytemp = 0;anyerr = 0;keyerr(1) = 0 + xpress = 0.0d0;xeps = 0.0d0;xtemp = 0.0d0;xye = 0.0d0 +! end EOS Omni vars + + myproc = CCTK_MyProc(cctkGH) + + ! save memory when MP is not used + if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then + g11 => gaa + g12 => gab + g13 => gac + g22 => gbb + g23 => gbc + g33 => gcc + vup => lvel + else + g11 => gxx + g12 => gxy + g13 => gxz + g22 => gyy + g23 => gyz + g33 => gzz + vup => vel + end if + + 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 = 0.0d0 + end if + + ! this is a poly call + xrho(1) = GRHydro_rho_min + call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& + xrho,xeps,xtemp,xye,xpress,keyerr,anyerr) + call EOS_Omni_EpsFromPress(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& + xrho,xeps,xtemp,xye,xpress,xeps,keyerr,anyerr) + pmin = xpress(1) + epsmin = xeps(1) + + !$OMP PARALLEL DO PRIVATE(i,j,k,itracer,& + !$OMP uxx, uxy, uxz, uyy, uyz, uzz, det, epsnegative, anyerr, keyerr, keytemp,& + !$OMP warnline, dummy1, dummy2) + do k = 1, nz + do j = 1, ny + do i = 1, nx + + + !do not compute if in atmosphere or in excised region + if ((atmosphere_mask(i,j,k) .ne. 0) .or. & + (hydro_excision_mask(i,j,k) .ne. 0)) cycle + + 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,& + g11(i,j,k),g12(i,j,k),g13(i,j,k),g22(i,j,k),& + g23(i,j,k),g33(i,j,k)) + + 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(evolve_Y_e.ne.0) then + Y_e(i,j,k) = max(min(Y_e_con(i,j,k) / dens(i,j,k),GRHydro_Y_e_max),& + GRHydro_Y_e_min) + endif + + 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)) + 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 + w_lorentz(i,j,k) = 1.d0 + + ! set hot atmosphere values + temperature(i,j,k) = grhydro_hot_atmo_temp + y_e(i,j,k) = grhydro_hot_atmo_Y_e + y_e_con(i,j,k) = y_e(i,j,k) * dens(i,j,k) + keytemp = 1 + call EOS_Omni_press(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,& + rho(i,j,k),eps(i,j,k),temperature(i,j,k),y_e(i,j,k),& + 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) + + cycle + + end if + + call Con2Prim_pt_hot3(cctk_iteration,myproc,i,j,k,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),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), & + 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) + + if(temperature(i,j,k).gt.GRHydro_max_temp) then + !$OMP CRITICAL + call CCTK_WARN(1,"C2P: Temperature too high") + write(warnline,"(i8)") cctk_iteration + call CCTK_WARN(1,warnline) + write(warnline,"(4i5,1P10E15.6)") GRHydro_Reflevel,i,j,k,x(i,j,k),y(i,j,k),z(i,j,k) + call CCTK_WARN(1,warnline) + write(warnline,"(1P10E15.6)") dens(i,j,k),scon(i,j,k,1:3),& + tau(i,j,k),w_lorentz(i,j,k) + call CCTK_WARN(1,warnline) + write(warnline,"(1P10E15.6)") rho(i,j,k),eps(i,j,k),& + temperature(i,j,k),Y_e(i,j,k) + call CCTK_WARN(1,warnline) + write(warnline,"(A7,i8)") "code: ",keyerr(1) + call CCTK_WARN(1,warnline) + write(warnline,"(A10,i5)") "reflevel: ", GRHydro_reflevel + call CCTK_WARN(1,warnline) + call CCTK_WARN(0,"Aborting!!!") + !$OMP END CRITICAL + endif + + + if( abs(GRHydro_C2P_failed(i,j,k)-2.0d0) .lt. 1.0d-10) then + ! this means c2p did not converge. + ! In this case, we attempt to call c2p with a reduced + ! accuracy requirement; if it fails again, we abort + GRHydro_C2P_failed(i,j,k) = 0 + if(temperature(i,j,k) .gt. GRHydro_hot_atmo_temp) then + local_perc_ptol = GRHydro_perc_ptol*100.0d0 + else + ! If we are in the extrapolation regime for the EOS, + ! we accept a larger c2p error, since we will be resetting + ! the temperature anyway. The error scale here is chosen + ! somewhat arbitrarily to be 0.01% in pressnew/pressold. + local_perc_ptol = 1.0d-4 + endif + call Con2Prim_pt_hot3(cctk_iteration,myproc,i,j,k,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),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), & + 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 + !$OMP CRITICAL + if (GRHydro_reflevel.ge.GRHydro_c2p_warn_from_reflevel) then + call CCTK_WARN(1,"Convergence problem in c2p") + write(warnline,"(A10,i5)") "reflevel: ",GRHydro_reflevel + call CCTK_WARN(1,warnline) + write(warnline,"(A10,i10)") "iteration: ",cctk_iteration + call CCTK_WARN(1,warnline) + write(warnline,"(A10,F5.2)") "error: ", GRHydro_C2P_failed(i,j,k) + call CCTK_WARN(1,warnline) + write(warnline,"(3i5,1P10E15.6)") i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),r(i,j,k) + call CCTK_WARN(1,warnline) + write(warnline,"(1P10E15.6)") rho(i,j,k),dens(i,j,k),eps(i,j,k),& + temperature(i,j,k),y_e(i,j,k) + call CCTK_WARN(1,warnline) + call CCTK_WARN(0,"Aborting!!!") + endif + !$OMP END CRITICAL + endif + endif + + if( abs(GRHydro_C2P_failed(i,j,k)-3.0d0) .lt. 1.0d-10 .or. & + temperature(i,j,k).lt.GRHydro_hot_atmo_temp) then + ! dropped out off the EOS table in temperature or below + ! the temperature of the atmosphere. + ! Now reset this point to minimum temperature. + GRHydro_C2P_failed(i,j,k) = 0 + !$OMP CRITICAL + write(warnline,"(A10,i7,4i3,1P10E15.6)") "reset T:",& + cctk_iteration,GRHydro_Reflevel,& + i,j,k,rho(i,j,k),& + eps(i,j,k),temperature(i,j,k),y_e(i,j,k) + call CCTK_WARN(1,warnline) + !$OMP END CRITICAL + temperature(i,j,k) = GRHydro_hot_atmo_temp + keytemp = 1 + call EOS_Omni_press(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,& + rho(i,j,k),eps(i,j,k),temperature(i,j,k),y_e(i,j,k),& + press(i,j,k),keyerr,anyerr) + keytemp = 0 + if(anyerr.ne.0) then + !$OMP CRITICAL + call CCTK_WARN(0,"EOS Problem in C2P hot!!!") + !$OMP END CRITICAL + endif + ! prim2con + w_lorentz(i,j,k) = & + 1.d0 / sqrt(1.d0 - (gxx(i,j,k)*vel(i,j,k,1)**2 & + + gyy(i,j,k)*vel(i,j,k,2)**2 & + + gzz(i,j,k) *vel(i,j,k,3)**2 & + + 2.0d0*gxy(i,j,k)*vel(i,j,k,1)*vel(i,j,k,2) & + + 2.0d0*gxz(i,j,k)*vel(i,j,k,1)*vel(i,j,k,3) & + + 2.0d0*gyz(i,j,k)*vel(i,j,k,2)*vel(i,j,k,3))) + vlowx = gxx(i,j,k)*vel(i,j,k,1) & + + gxy(i,j,k)*vel(i,j,k,2) & + + gxz(i,j,k)*vel(i,j,k,3) + vlowy = gxy(i,j,k)*vel(i,j,k,1) & + + gyy(i,j,k)*vel(i,j,k,2) & + + gyz(i,j,k)*vel(i,j,k,3) + 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)) & + + 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)) & + +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)) & + + 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)) & + + press(i,j,k))*w_lorentz(i,j,k)**2 - press(i,j,k)) & + - dens(i,j,k) + endif + + end do + end do + end do + !$OMP END PARALLEL DO + + return + +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, & + GRHydro_reflevel, GRHydro_C2P_failed, local_perc_ptol) + + implicit none + + 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, & + y, z, r, GRHydro_rho_min + CCTK_REAL temp, ye + CCTK_REAL s2, f, df, vlowx, vlowy, vlowz + CCTK_INT cctk_iteration, ii,jj,kk,count, i, handle, GRHydro_reflevel + CCTK_REAL GRHydro_C2P_failed + CCTK_REAL udens, usx, usy, usz, utau, pold, pnew, & + temp1, drhobydpress, depsbydpress, dpressbydeps, dpressbydrho, pmin, epsmin + CCTK_REAL pminl,plow,tmp, dummy1, dummy2 + CCTK_REAL local_perc_ptol + CCTK_REAL dpf + + integer :: myproc + + character(len=256) warnline + logical epsnegative, mustbisect + + integer :: failwarnmode + integer :: failinfomode + +! begin EOS Omni vars + integer :: n,keytemp,anyerr,keyerr(1) + real*8 :: xpress,temp0 + integer :: nfudgemax,nf + n=1;keytemp=0;anyerr=0;keyerr(1)=0 + temp0 = 0.0d0;xpress = 0.0d0 + nf=0;nfudgemax=30 +! end EOS Omni vars + + mustbisect = .false. + + failinfomode = 1 + failwarnmode = 0 + +! set pmin to something sensible: + pminl = 1.0d-28 + + if(con2prim_oct_hack.ne.0.and.& + (x .lt. -1.0d-12 .or.& + y .lt. -1.0d-12 .or.& + z .lt. -1.0d-12)) then + failwarnmode = 2 + failinfomode = 2 + endif + +!!$ 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: + temp0 = temp + call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,& + rho,epsilon,temp,ye,xpress,keyerr,anyerr) + ! This is the lowest admissible pressure, otherwise we are off the physical regime + ! Sometimes, we may end up with bad initial guesses. + ! In that case we need to start off with something + ! reasonable at least + plow = max(pminl, sqrt(s2) - utau - udens) + ! Start out with some reasonable initial guess + pold = max(plow+1.d-10,xpress) + ! error handling + if(anyerr.ne.0) then + if(keyerr(1).eq.668) then + continue +! GRHydro_C2P_failed = 3.0d0 +! write(warnline,"(A3,3i3,1P10E15.6)") "*1", ii,jj,kk,rho,epsilon,temp,ye,xpress +! call CCTK_WARN(1,warnline) + else + !$OMP CRITICAL + call CCTK_WARN(failinfomode,"EOS error in c2p 0") + write(warnline,"(4i5,1P10E15.6)") GRHydro_Reflevel,ii,jj,kk,x,y,z + call CCTK_WARN(failinfomode,warnline) + write(warnline,"(1P10E15.6)") rho,dens,epsilon,temp,temp0,ye + call CCTK_WARN(failinfomode,warnline) + write(warnline,"(A7,i8)") "code: ",keyerr(1) + call CCTK_WARN(failinfomode,warnline) + write(warnline,"(A10,i5)") "reflevel: ", GRHydro_reflevel + call CCTK_WARN(failinfomode,warnline) + call CCTK_WARN(failwarnmode,"Aborting!!!") + !$OMP END CRITICAL + endif + endif + +!!$ 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.0d0 + !$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 + +!!$ Calculate the function + + call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,& + rho,epsilon,temp,ye,xpress,keyerr,anyerr) + ! error handling + if(anyerr.ne.0) then + if(keyerr(1).eq.668) then + continue +! GRHydro_C2P_failed = 3.0d0 +! write(warnline,"(A3,3i3,1P10E15.6)") "*2", ii,jj,kk,rho,epsilon,temp,ye,xpress +! call CCTK_WARN(1,warnline) + else + !$OMP CRITICAL + call CCTK_WARN(failinfomode,"EOS error in c2p 1") + write(warnline,"(4i5,1P10E15.6)") GRHydro_Reflevel,ii,jj,kk,x,y,z + call CCTK_WARN(failinfomode,warnline) + write(warnline,"(1P10E15.6)") rho,dens,epsilon,temp,temp0,ye + call CCTK_WARN(failinfomode,warnline) + write(warnline,"(A7,i8)") "code: ",keyerr(1) + call CCTK_WARN(failinfomode,warnline) + write(warnline,"(A10,i5)") "reflevel: ", GRHydro_reflevel + call CCTK_WARN(failinfomode,warnline) + call CCTK_WARN(failwarnmode,"Aborting!!!") + !$OMP END CRITICAL + endif + endif + + f = pold - xpress + + if (f .ne. f) then + ! Ok, this yielded nonsense, let's enforce bisection! + mustbisect = .true. + endif + +!!$Find the root + + count = 0 + pnew = pold + do while ( ((abs(pnew - pold)/abs(pnew) .gt. local_perc_ptol) .and. & + (abs(pnew - pold) .gt. GRHydro_del_ptol)) .or. & + (count .lt. GRHydro_countmin)) + count = count + 1 + +#if 0 + if(myproc.eq.3.and.& + cctk_iteration.eq.11929.and.ii.eq.64.and.jj.eq.19.and.kk.eq.29) then + !$OMP CRITICAL + write(warnline,"(i5,1P10E18.9)") count, (abs(pnew - pold)/abs(pnew)), epsilon, temp, press + call CCTK_WARN(1,warnline) + !$OMP END CRITICAL + endif +#endif + + if (count > GRHydro_countmax) then + ! non-convergence is now handled outside of this + ! routine + GRHydro_C2P_failed = 2.0d0 + temp = temp0 + return + end if + + call EOS_Omni_dpderho_dpdrhoe(handle,keytemp,GRHydro_eos_rf_prec,& + n,rho,epsilon,temp,ye,dpressbydeps,dpressbydrho,keyerr,anyerr) + + 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 + + ! Try to obtain new pressure via Newton-Raphson. + + pnew = pold - f/df + + ! Check if Newton-Raphson resulted in something reasonable! + if (c2p_resort_to_bisection.ne.0) then + + tmp = (utau + pnew + udens)**2 - s2 + plow = max(pminl, sqrt(s2) - utau - udens) + + if (pnew .lt. plow .or. tmp .le. 0.0d0 .or. mustbisect) then + + ! Ok, Newton-Raphson ended up finding something unphysical. + ! Let's try to find our root via bisection (which converges slower but is more robust) + + pnew = (plow + pold) / 2 + tmp = (utau + pnew + udens)**2 - s2 + + mustbisect = .false. + end if + + else + + ! This is the standard algorithm without resorting to bisection. + + pnew = max(pminl, pnew) + tmp = (utau + pnew + udens)**2 - s2 + + endif + + ! Check if we are still in the physical range now. + ! If not, we set the C2P failed mask, and set pnew = pold (which will exit the loop). + + if ((tmp .le. 0.0d0) .and. GRHydro_C2P_failed .eq. 0) then + GRHydro_C2P_failed = 1.0d0 + pnew = pold + endif + +!!$ 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 + + call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,& + rho,epsilon,temp,ye,xpress,keyerr,anyerr) + ! error handling + if(anyerr.ne.0) then + if(keyerr(1).eq.668) then + continue +! GRHydro_C2P_failed = 3.0d0 +! write(warnline,"(A3,3i3,1P10E15.6)") "*3", ii,jj,kk,rho,epsilon,temp,ye,xpress +! call CCTK_WARN(1,warnline) + else + !$OMP CRITICAL + call CCTK_WARN(failinfomode,"EOS error in c2p 2") + write(warnline,"(4i5,1P10E15.6)") GRHydro_Reflevel,ii,jj,kk,x,y,z + call CCTK_WARN(failinfomode,warnline) + write(warnline,"(1P10E15.6)") rho,dens,epsilon,temp,temp0,ye + call CCTK_WARN(failinfomode,warnline) + write(warnline,"(A7,i8)") "code: ",keyerr(1) + call CCTK_WARN(failinfomode,warnline) + write(warnline,"(A10,i5)") "reflevel: ", GRHydro_reflevel + call CCTK_WARN(failinfomode,warnline) + call CCTK_WARN(failwarnmode,"Aborting!!!") + !$OMP END CRITICAL + endif + endif + + f = pnew - xpress + + if (f .ne. f) then + ! Ok, this yielded nonsense, let's enforce bisection! + mustbisect = .true. + endif + + enddo + + +!!$ Polish the root + + do i=1,GRHydro_polish + + call EOS_Omni_DPressByDRho(handle,keytemp,GRHydro_eos_rf_prec,n,& + rho,epsilon,temp,ye,dpressbydrho,keyerr,anyerr) + + call EOS_Omni_DPressByDEps(handle,keytemp,GRHydro_eos_rf_prec,n,& + rho,epsilon,temp,ye,dpressbydeps,keyerr,anyerr) + + 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 + + call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n, & + rho,epsilon,temp,ye,xpress,keyerr,anyerr) + + ! error handling + if(anyerr.ne.0) then + if(keyerr(1).eq.668) then + continue +! GRHydro_C2P_failed = 3.0d0 +! write(warnline,"(A3,3i3,1P10E15.6)") "*4", ii,jj,kk,rho,epsilon,temp,ye,xpress +! call CCTK_WARN(1,warnline) + else + !$OMP CRITICAL + call CCTK_WARN(failinfomode,"EOS error in c2p 3") + write(warnline,"(4i5,1P10E15.6)") GRHydro_Reflevel,ii,jj,kk,x,y,z + call CCTK_WARN(failinfomode,warnline) + write(warnline,"(1P10E15.6)") rho,dens,epsilon,temp,temp0,ye + call CCTK_WARN(failinfomode,warnline) + write(warnline,"(A7,i8)") "code: ",keyerr(1) + call CCTK_WARN(failinfomode,warnline) + write(warnline,"(A10,i5)") "reflevel: ", GRHydro_reflevel + call CCTK_WARN(failinfomode,warnline) + call CCTK_WARN(failwarnmode,"Aborting!!!") + !$OMP END CRITICAL + endif + endif + + f = pold - xpress + + enddo + +!!$ Calculate primitive variables from root + + !if (rho .le. GRHydro_rho_min*(1.d0+GRHydro_atmo_tolerance) ) then + 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 + temp = GRHydro_hot_atmo_temp + ye = GRHydro_hot_atmo_Y_e + ye_con = dens * ye + keytemp=1 + call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n, & + rho,epsilon,temp,ye,xpress,keyerr,anyerr) + keytemp=0 + ! w_lorentz=1, so the expression for utau reduces to: + utau = rho + rho*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 + + 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 + + ! indicate to wrapper routine that we dropped out of the + ! EOS table + if(keyerr(1).eq.668) then + GRHydro_C2P_failed = 3.0d0 + endif + +end subroutine Con2Prim_pt_hot3 + |