diff options
Diffstat (limited to 'src/GRHydro_Con2PrimM.F90')
| -rw-r--r-- | src/GRHydro_Con2PrimM.F90 | 141 |
1 files changed, 114 insertions, 27 deletions
diff --git a/src/GRHydro_Con2PrimM.F90 b/src/GRHydro_Con2PrimM.F90 index 2712397..1b0e831 100644 --- a/src/GRHydro_Con2PrimM.F90 +++ b/src/GRHydro_Con2PrimM.F90 @@ -57,11 +57,11 @@ subroutine Conservative2PrimitiveM(CCTK_ARGUMENTS) character(len=100) warnline CCTK_REAL :: local_min_tracer, local_gam(1), local_pgam,local_K, sc + CCTK_REAL :: local_perc_ptol ! begin EOS Omni vars integer :: n,keytemp,anyerr,keyerr(1) real*8 :: xpress(1),xtemp(1),xye(1),xeps(1),xrho(1),one(1)=1.0d0 - n=1;keytemp=0;anyerr=0;keyerr(1)=0 xpress(1)=0.0d0;xtemp(1)=0.0d0;xye(1)=0.0d0;xeps(1)=0.0d0 ! end EOS Omni vars @@ -76,12 +76,16 @@ subroutine Conservative2PrimitiveM(CCTK_ARGUMENTS) local_min_tracer = 0d0 end if +! call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& +! GRHydro_rho_min,xeps,xtemp,xye,pmin,keyerr,anyerr) +! call EOS_Omni_EpsFromPress(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& +! GRHydro_rho_min,xeps,xtemp,xye,pmin,epsmin,keyerr,anyerr) ! 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,pmin,keyerr,anyerr) call EOS_Omni_EpsFromPress(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& - xrho,xeps,xtemp,xye,xpress,epsmin,keyerr,anyerr) + xrho,xeps,xtemp,xye,pmin,epsmin,keyerr,anyerr) call EOS_Omni_press(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,& one,one,xtemp,xye,local_gam,keyerr,anyerr) @@ -89,7 +93,7 @@ subroutine Conservative2PrimitiveM(CCTK_ARGUMENTS) !$OMP PARALLEL DO PRIVATE(i,j,k,itracer,& !$OMP uxx, uxy, uxz, uyy, uyz, uzz, det, epsnegative, & - !$OMP b2,xrho,xeps,xpress,xtemp,local_K,local_pgam,sc,keyerr,anyerr ) + !$OMP b2,xrho,xeps,xpress,xtemp,local_K,local_pgam,sc,keyerr,anyerr,keytemp,local_perc_ptol) do k = 1, nz do j = 1, ny do i = 1, nx @@ -122,11 +126,19 @@ subroutine Conservative2PrimitiveM(CCTK_ARGUMENTS) endif if(evolve_Y_e.ne.0) then - Y_e(i,j,k) = Y_e_con(i,j,k) / dens(i,j,k) + 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 ( dens(i,j,k) .le. sqrt(det)*GRHydro_rho_min*(1.d0+GRHydro_atmo_tolerance) ) then + !call CCTK_WARN(1,"Con2Prim: Resetting to atmosphere") + !write(warnline,"(3i5,1P10E15.6)") i,j,k,x(i,j,k),y(i,j,k),z(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) + b2=gxx(i,j,k)*Bvec(i,j,k,1)**2+gyy(i,j,k)*Bvec(i,j,k,2)**2+gzz(i,j,k)*Bvec(i,j,k,3)**2+ & 2.0*(gxy(i,j,k)*Bvec(i,j,k,1)*Bvec(i,j,k,2)+gxz(i,j,k)*Bvec(i,j,k,1)*Bvec(i,j,k,3)+ & gyz(i,j,k)*Bvec(i,j,k,2)*Bvec(i,j,k,3)) @@ -136,18 +148,34 @@ subroutine Conservative2PrimitiveM(CCTK_ARGUMENTS) scon(i,j,k,:) = 0.d0 vel(i,j,k,:) = 0.d0 w_lorentz(i,j,k) = 1.d0 - xtemp = 0.0d0 - - call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& - rho(i,j,k),eps(i,j,k),xtemp,xye,press(i,j,k),keyerr,anyerr) - call EOS_Omni_EpsFromPress(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& + if(evolve_temper.ne.0) then + ! 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) + else + keytemp = 0 + call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& + rho(i,j,k),eps(i,j,k),xtemp,xye,press(i,j,k),keyerr,anyerr) + call EOS_Omni_EpsFromPress(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& rho(i,j,k),eps(i,j,k),xtemp,xye,press(i,j,k),eps(i,j,k),keyerr,anyerr) + endif + + !call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& + ! rho(i,j,k),eps(i,j,k),xtemp,xye,press(i,j,k),keyerr,anyerr) + ! + !call EOS_Omni_EpsFromPress(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& + ! 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 does need to take into account the existing B-field -!!$ with w_lorentz=1, we find tau = sqrtdet*(rho (1+eps+b^2/2)) - dens [Press drops out] + !!$ tau does need to take into account the existing B-field + !!$ with w_lorentz=1, we find tau = sqrtdet*(rho (1+eps+b^2/2)) - dens [Press drops out] tau(i,j,k) = sqrt(det) * (rho(i,j,k)*(1.0+eps(i,j,k)+b2/2.0)) - dens(i,j,k) cycle @@ -155,17 +183,65 @@ subroutine Conservative2PrimitiveM(CCTK_ARGUMENTS) end if if(evolve_temper.eq.0) then - call GRHydro_Con2PrimM_pt(GRHydro_eos_handle, local_gam(1), dens(i,j,k), & + call GRHydro_Con2PrimM_pt(GRHydro_eos_handle, keytemp, GRHydro_eos_rf_prec, local_gam(1), 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),& + Bcons(i,j,k,1),Bcons(i,j,k,2),Bcons(i,j,k,3), xye(1), xtemp(1), 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),& Bvec(i,j,k,1), Bvec(i,j,k,2), Bvec(i,j,k,3),b2, w_lorentz(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), & uxx,uxy,uxz,uyy,uyz,uzz,det, & epsnegative,GRHydro_C2P_failed(i,j,k)) else - stop "Please implement finite T con2prim routine in MHD part!" - endif + keytemp = 1 + !call Con2Prim_pt_hot(cctk_iteration,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),vel(i,j,k,1),& + ! vel(i,j,k,2),vel(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) + ! Bvec(i,j,k,1) = 0.0d0 + ! Bvec(i,j,k,2) = 0.0d0 + ! Bvec(i,j,k,3) = 0.0d0 + call GRHydro_Con2PrimM_pt(GRHydro_eos_handle, keytemp, GRHydro_eos_rf_prec, local_gam(1), 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), Y_e(i,j,k), temperature(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),& + Bvec(i,j,k,1), Bvec(i,j,k,2), Bvec(i,j,k,3), b2, w_lorentz(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), & + uxx,uxy,uxz,uyy,uyz,uzz,det, & + epsnegative,GRHydro_C2P_failed(i,j,k)) + if(GRHydro_C2P_failed(i,j,k).eq.2) 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 + local_perc_ptol = GRHydro_eos_rf_prec*100.0d0 + call GRHydro_Con2PrimM_pt(GRHydro_eos_handle, keytemp, local_perc_ptol, local_gam(1), 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), Y_e(i,j,k), temperature(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),& + Bvec(i,j,k,1), Bvec(i,j,k,2), Bvec(i,j,k,3), b2, w_lorentz(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), & + uxx,uxy,uxz,uyy,uyz,uzz,det, & + epsnegative,GRHydro_C2P_failed(i,j,k)) + if(GRHydro_C2P_failed(i,j,k).eq.2) 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,"(3i5,1P10E15.6)") i,j,k,x(i,j,k),y(i,j,k),z(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 + endif if (epsnegative .ne. 0) then @@ -244,7 +320,8 @@ subroutine Conservative2PrimitiveM(CCTK_ARGUMENTS) end if - if ( rho(i,j,k) .le. GRHydro_rho_min*(1.d0+GRHydro_atmo_tolerance) .or. GRHydro_C2P_failed(i,j,k) .ge. 1) then + if ( rho(i,j,k) .le. GRHydro_rho_min*(1.d0+GRHydro_atmo_tolerance)) then +! if ( rho(i,j,k) .le. GRHydro_rho_min*(1.d0+GRHydro_atmo_tolerance) .or. GRHydro_C2P_failed(i,j,k) .ge. 1) then b2=gxx(i,j,k)*Bvec(i,j,k,1)**2+gyy(i,j,k)*Bvec(i,j,k,2)**2+gzz(i,j,k)*Bvec(i,j,k,3)**2+ & 2.0*(gxy(i,j,k)*Bvec(i,j,k,1)*Bvec(i,j,k,2)+gxz(i,j,k)*Bvec(i,j,k,1)*Bvec(i,j,k,3)+ & @@ -255,13 +332,23 @@ subroutine Conservative2PrimitiveM(CCTK_ARGUMENTS) scon(i,j,k,:) = 0.d0 vel(i,j,k,:) = 0.d0 w_lorentz(i,j,k) = 1.d0 - xtemp = 0.0d0 - - call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& - rho(i,j,k),eps(i,j,k),xtemp,xye,press(i,j,k),keyerr,anyerr) - call EOS_Omni_EpsFromPress(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& + if(evolve_temper.ne.0) then + ! 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) + else + keytemp = 0 + call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& + rho(i,j,k),eps(i,j,k),xtemp,xye,press(i,j,k),keyerr,anyerr) + call EOS_Omni_EpsFromPress(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& rho(i,j,k),eps(i,j,k),xtemp,xye,press(i,j,k),eps(i,j,k),keyerr,anyerr) + endif ! w_lorentz=1, so the expression for tau reduces to: @@ -401,10 +488,9 @@ subroutine Conservative2PrimitiveBoundsM(CCTK_ARGUMENTS) Y_e(i,j,k) = Y_e_con(i,j,k) / dens(i,j,k) endif - - call GRHydro_Con2PrimM_pt(GRHydro_eos_handle, local_gam(1),densminus(i,j,k), & + call GRHydro_Con2PrimM_pt(GRHydro_eos_handle, keytemp, GRHydro_eos_rf_prec, local_gam(1),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),& + Bconsxminus(i,j,k), Bconsyminus(i,j,k), Bconszminus(i,j,k), 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, & @@ -459,9 +545,10 @@ subroutine Conservative2PrimitiveBoundsM(CCTK_ARGUMENTS) endif epsnegative = 0 - call GRHydro_Con2PrimM_pt(GRHydro_eos_handle, local_gam(1), 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),& + + call GRHydro_Con2PrimM_pt(GRHydro_eos_handle, keytemp, GRHydro_eos_rf_prec, local_gam(1),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), 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, & |