diff options
author | rhaas <rhaas@c83d129a-5a75-4d5a-9c4d-ed3a5855bf45> | 2012-08-27 19:19:17 +0000 |
---|---|---|
committer | rhaas <rhaas@c83d129a-5a75-4d5a-9c4d-ed3a5855bf45> | 2012-08-27 19:19:17 +0000 |
commit | 9df45a235ce39bfd46d83369f2c4ae859088abf9 (patch) | |
tree | 902146e0f8ccf3f96bb2a8cc01cda3b61224cdce /src | |
parent | aa3f5cca41b9b67e16cde545cf14a93633f3d6cf (diff) |
GRHydro: Changes to make GRHydro MHD work with nuclear/hot equation of state:
1) Switch to EOSOmni pointwise C2P routine and modify where necessary.
2) Modify Con2PrimM.F90 to allow for the evolution of temperature and adjust the wrapper routine.
3) Create EigenProblemM_hot pointwise routine and call that from HLLEM.F90 when temperature is evolved.
Additionally adjust HLLEM where necessary.
4) Adjust InterfacesM.h to incorporate the newly created functions.
5) Fix a loop problem (not taking into account constraint transport) in PPM reconstruction of Y_e
6) Introduce Prim2ConM_hot and call this pointwise routine from Prim2ConM.F90 when temperature is evolved.
Additionally also make this routine available to initial data routine in GRHydro_InitData
7) Adjust loops in GRHydro_PoloidalMagFieldM.F90 to not set boundary points it cannot set but instead call boundary group afterwards! Pay attention as this will not work with boundary conditions set to "none" in MHD case anymore but is the correct thing to do.
8) Allow StarMapper to extend HydroBase::initial_hydro = "starmapper".
9) Smaller fixes.
From: Philipp Moesta <pmoesta@tapir.caltech.edu>
git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEvolve/GRHydro/trunk@410 c83d129a-5a75-4d5a-9c4d-ed3a5855bf45
Diffstat (limited to 'src')
-rw-r--r-- | src/GRHydro_Analysis.F90 | 4 | ||||
-rw-r--r-- | src/GRHydro_CalcUpdate.F90 | 18 | ||||
-rw-r--r-- | src/GRHydro_Con2PrimM.F90 | 141 | ||||
-rw-r--r-- | src/GRHydro_Con2PrimM_pt_EOSOmni.c | 43 | ||||
-rw-r--r-- | src/GRHydro_EigenproblemM.F90 | 85 | ||||
-rw-r--r-- | src/GRHydro_HLLEM.F90 | 173 | ||||
-rw-r--r-- | src/GRHydro_InterfacesM.h | 7 | ||||
-rw-r--r-- | src/GRHydro_PPMMReconstruct_drv.F90 | 12 | ||||
-rw-r--r-- | src/GRHydro_Prim2ConM.F90 | 271 | ||||
-rw-r--r-- | src/GRHydro_UpdateMaskM.F90 | 205 | ||||
-rw-r--r-- | src/make.code.defn | 6 |
11 files changed, 790 insertions, 175 deletions
diff --git a/src/GRHydro_Analysis.F90 b/src/GRHydro_Analysis.F90 index 918f91a..3bed392 100644 --- a/src/GRHydro_Analysis.F90 +++ b/src/GRHydro_Analysis.F90 @@ -76,9 +76,9 @@ subroutine GRHydro_CalcDivB(CCTK_ARGUMENTS) Bcons(i,j+1,k,2)) Bcons_r3 = 0.5d0 * (Bcons(i,j,k,3) + & Bcons(i,j,k+1,3)) - !write(*,*) "Bcons_r1 = ", Bcons_r1 + divB(i,j,k) = divB(i,j,k) + (Bcons_l1 - Bcons_r1) * idx + (Bcons_l2 - Bcons_r2) * idy + (Bcons_l3 - Bcons_r3) * idz - !divB(i,j,k) = 0.0d0 + endif endif endif diff --git a/src/GRHydro_CalcUpdate.F90 b/src/GRHydro_CalcUpdate.F90 index f0eca19..32ac21d 100644 --- a/src/GRHydro_CalcUpdate.F90 +++ b/src/GRHydro_CalcUpdate.F90 @@ -110,24 +110,6 @@ subroutine UpdateCalculation(CCTK_ARGUMENTS) (alp_l * psidcflux(i-xoffset,j-yoffset,k-zoffset) - & alp_r * psidcflux(i,j,k)) * idx endif - if(track_divB.ne.0) then - if(transport_constraints.ne.0) then - ! edge based divergence (see WhiskyMHD & Bruno's thesis, Eq. 7.27) - ! FIXME: this uses the Bcons before the current MoL substep update, should be in MoL_PseudoEvolution instead - divB(i,j,k) = divB(i,j,k) + & - 0.25d0*(Bcons(i+xoffset,j+yoffset ,k+zoffset,flux_direction)-Bcons(i ,j ,k ,flux_direction)+ & - Bcons(i+1 ,j+1-zoffset,k+zoffset,flux_direction)-Bcons(i+1-xoffset,j+xoffset ,k ,flux_direction)+ & - Bcons(i+xoffset,j+1-xoffset,k+1 ,flux_direction)-Bcons(i ,j+zoffset ,k+1-zoffset,flux_direction)+ & - Bcons(i+1 ,j+1 ,k+1 ,flux_direction)-Bcons(i+1-xoffset,j+1-yoffset,k+1-zoffset,flux_direction))*idx - else - Bcons_l = 0.5d0 * (Bcons(i,j,k,flux_direction) + & - Bcons(i-xoffset,j-yoffset,k-zoffset,flux_direction)) - Bcons_r = 0.5d0 * (Bcons(i,j,k,flux_direction) + & - Bcons(i+xoffset,j+yoffset,k+zoffset,flux_direction)) - divB(i,j,k) = divB(i,j,k) + (Bcons_l - Bcons_r ) * idx - endif - endif - endif if (evolve_tracer .ne. 0) then 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, & diff --git a/src/GRHydro_Con2PrimM_pt_EOSOmni.c b/src/GRHydro_Con2PrimM_pt_EOSOmni.c index f0b6a39..4a52414 100644 --- a/src/GRHydro_Con2PrimM_pt_EOSOmni.c +++ b/src/GRHydro_Con2PrimM_pt_EOSOmni.c @@ -45,6 +45,7 @@ #include <complex.h> #include "cctk.h" +#include "cctk_Parameters.h" /* Set this to be 1 if you want debug output */ #define DEBUG_CON2PRIMM (0) @@ -123,7 +124,7 @@ static CCTK_REAL pressure_rho0_w(CCTK_REAL rho0, CCTK_REAL w,CCTK_REAL gammaeos) } void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt) ( - CCTK_INT *handle, CCTK_INT *keytemp, CCTK_REAL *prec, + CCTK_INT *handle, CCTK_INT *keytemp, CCTK_REAL *prec, CCTK_REAL *gamma_eos, CCTK_REAL *dens_in, CCTK_REAL *sx_in, CCTK_REAL *sy_in, CCTK_REAL *sz_in, @@ -204,6 +205,8 @@ RETURN VALUE: of retval = (i*100 + j) where ( used to be 5 -> failure: rho,uu <= 0 but now sets epsnegative to non-zero ) **********************************************************************************/ + + void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt) ( CCTK_INT *handle, CCTK_INT *keytemp, CCTK_REAL *prec, CCTK_REAL *gamma_eos, @@ -239,11 +242,13 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt) ( CCTK_REAL inv_sqrt_detg = 1./sqrt_detg; CCTK_INT i,j, i_increase ; + DECLARE_CCTK_PARAMETERS; + struct LocGlob lg; struct eosomnivars eosvars; eosvars.eoshandle = *handle; - printf("handle = %i\n",*handle); + //printf("handle = %i\n",*handle); eosvars.eoskeytemp = *keytemp; eosvars.eosprec = *prec; eosvars.eos_y_e[0] = *y_e_in; @@ -271,6 +276,8 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt) ( fprintf(stdout," *vely = %26.16e \n", *vely ); fprintf(stdout," *velz = %26.16e \n", *velz ); fprintf(stdout," *epsilon = %26.16e \n", *epsilon ); + fprintf(stdout," *temp_in = %26.16e \n", *temp_in ); + fprintf(stdout," *y_e_in = %26.16e \n", *y_e_in ); fprintf(stdout," *pressure = %26.16e \n", *pressure ); fprintf(stdout," *Bx = %26.16e \n", *Bx ); fprintf(stdout," *By = %26.16e \n", *By ); @@ -347,6 +354,7 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt) ( } if(vsq < 0. || vsq > 1. ) { *retval = 2.; + fprintf(stdout," *retval = %26.16e \n", *retval ); return; } @@ -357,7 +365,8 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt) ( // i.e. you don't get positive values for dP/d(vsq) . rho0 = lg.D / gamma ; u = (*epsilon) * rho0; - + CCTK_REAL uold = u; + CCTK_REAL dum1,dum2; p = pressure_rho0_eps_eosomni(rho0,*epsilon,&dum1,&dum2,&eosvars) ; // EOSOMNI // p = pressure_rho0_u(rho0,u,gammaeos) ; // EOS @@ -365,6 +374,7 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt) ( W_last = w*gammasq ; + //fprintf(stdout," p = %26.16e \n", p ); // Make sure that W is large enough so that v^2 < 1 : i_increase = 0; @@ -395,11 +405,13 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt) ( /* Problem with solver, so return denoting error before doing anything further */ if( ((*retval) != 0.) || (W == FAIL_VAL) ) { *retval = *retval*100.+1.; + fprintf(stdout," *retval = %26.16e \n", *retval ); return; } else{ if(W <= 0. || W > W_TOO_BIG) { *retval = 3.; + fprintf(stdout," *retval = %26.16e \n", *retval ); return; } } @@ -407,6 +419,7 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt) ( // Calculate v^2: if( vsq >= 1. ) { *retval = 4.; + fprintf(stdout," *retval = %26.16e \n", *retval ); return; } @@ -425,13 +438,20 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt) ( u=x_3d[2]; *epsilon = u/rho0; p = pressure_rho0_eps_eosomni(rho0,*epsilon,&dum1,&dum2,&eosvars) ; // EOSOMNI - printf("%g %g %g %g\n",rho0,u,*epsilon,p); + // printf("%g %g %g %g\n",rho0,u,*epsilon,p); } // User may want to handle this case differently, e.g. do NOT return upon // a negative rho/u, calculate v^i so that rho/u can be floored by other routine: - if( (rho0 <= 0.) || (u <= 0.) ) { - *epsnegative = 1; + // HOT: if( (rho0 <= 0.) || (u <= 0.) ) { + if(rho0 <= 0.) { +// *epsnegative = 1; + fprintf(stdout," *epsnegative = %10d \n", *epsnegative ); + fprintf(stdout," rho0 = %26.16e \n", rho0 ); + fprintf(stdout," u = %26.16e \n", u ); + fprintf(stdout," W = %26.16e \n", W ); + fprintf(stdout," vsq = %26.16e \n", vsq ); + fprintf(stdout," uold = %26.16e \n", uold ); return; } @@ -447,6 +467,13 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt) ( *vely = g_o_WBsq * ( usy + QdB_o_W*(*By) ) ; *velz = g_o_WBsq * ( usz + QdB_o_W*(*Bz) ) ; + if (*rho <= rho_abs_min*(1.0+GRHydro_atmo_tolerance) ) { + *rho = rho_abs_min; + *velx = 0.0; + *vely = 0.0; + *velz = 0.0; + *w_lorentz = 1.0; + } #if(DEBUG_CON2PRIMM) fprintf(stdout,"rho = %26.16e \n",*rho ); @@ -696,7 +723,7 @@ static CCTK_INT threed_newton_raphson_omni( CCTK_REAL x[], struct eosomnivars *e x[1] += dx[1] ; x[2] += dx[2] ; - printf("Updating vars: %g %g %g %g %g %g\n",x[0],dx[0],x[1],dx[1],x[2],dx[2]); + //printf("Updating vars: %g %g %g %g %g %g\n",x[0],dx[0],x[1],dx[1],x[2],dx[2]); /****************************************/ /* Calculate the convergence criterion */ @@ -717,7 +744,7 @@ static CCTK_INT threed_newton_raphson_omni( CCTK_REAL x[], struct eosomnivars *e if( x[1] > 1. ) { x[1] = dv; } } - if( x[2] < 0. ) { x[2] = 0.; } + // HOT: if( x[2] < 0. ) { x[2] = 0.; } /*****************************************************************************/ /* If we've reached the tolerance level, then just do a few extra iterations */ diff --git a/src/GRHydro_EigenproblemM.F90 b/src/GRHydro_EigenproblemM.F90 index 31a6c51..ec4d44a 100644 --- a/src/GRHydro_EigenproblemM.F90 +++ b/src/GRHydro_EigenproblemM.F90 @@ -138,6 +138,91 @@ subroutine eigenvaluesM(handle,rho,velx,vely,velz,eps,w_lorentz,& end subroutine eigenvaluesM +subroutine eigenvaluesM_hot(handle,ii,jj,kk,rho,velx,vely,velz,eps,w_lorentz,& + Bvcx,Bvcy,Bvcz,temp,ye,lam,gxx,gxy,gxz,gyy,gyz,gzz,u,alp,beta) + implicit none + + DECLARE_CCTK_PARAMETERS + + CCTK_REAL rho,velx,vely,velz,eps,w_lorentz + CCTK_REAL Bvcx,Bvcy,Bvcz + CCTK_REAL lam(5) + CCTK_REAL gxx,gxy,gxz,gyy,gyz,gzz + CCTK_REAL alp,beta,u + CCTK_REAL temp,ye + + CCTK_REAL cs2,one,two,U2 + CCTK_REAL vlowx,vlowy,vlowz,v2,w + CCTK_REAL lam1,lam2,lam3,lamm,lamp,lamm_nobeta,lamp_nobeta + CCTK_INT handle, ii,jj,kk + CCTK_REAL dpdrho,dpdeps,press + + CCTK_REAL Bvcxlow,Bvcylow,Bvczlow,Bvc2,rhohstar,va2 + CCTK_REAL Bdotv,b2 + +! begin EOS Omni vars + integer :: n,keytemp,anyerr,keyerr(1) + real*8 :: xpress,xeps,xtemp,xye + 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 + + one = 1.0d0 + two = 2.0d0 + +!!$ Set required fluid quantities + + call EOS_Omni_cs2(handle,keytemp,GRHydro_eos_rf_prec,n,& + rho,eps,temp,ye,cs2,keyerr,anyerr) + + vlowx = gxx*velx + gxy*vely + gxz*velz + vlowy = gxy*velx + gyy*vely + gyz*velz + vlowz = gxz*velx + gyz*vely + gzz*velz + v2 = vlowx*velx + vlowy*vely + vlowz*velz + +!!$ Lower the B-field, and square of the magnitude + Bvcxlow = gxx*Bvcx + gxy*Bvcy + gxz*Bvcz + Bvcylow = gxy*Bvcx + gyy*Bvcy + gyz*Bvcz + Bvczlow = gxz*Bvcx + gyz*Bvcy + gzz*Bvcz + Bvc2 = Bvcxlow*Bvcx + Bvcylow*Bvcy + Bvczlow*Bvcz + + Bdotv = Bvcxlow*velx + Bvcylow*vely + Bvczlow*velz + w = w_lorentz + + b2=Bvc2/w**2+Bdotv**2 + + +!!$ rhohstar is the term that appears in Tmunu as well = rho*enth + b^2 + rhohstar = rho*(1.0+eps)+press+b2 + +!!$ Alfven velocity squared + va2 = b2/(rhohstar) + +!!$ The following combination always comes up in the wavespeed calculation: +!!$ U2 = v_a^2 + c_s^2(1-v_a^2) +!!$ In the unmagnetized case, it reduces to cs2 + U2 = va2+cs2*(1.d0-va2) + +!!$ Calculate eigenvalues + + lam1 = velx - beta/alp + lam2 = velx - beta/alp + lam3 = velx - beta/alp + lamp_nobeta = (velx*(one-U2) + sqrt(U2*(one-v2)*& + (u*(one-v2*U2) - velx**2*(one-U2))))/(one-v2*U2) + lamm_nobeta = (velx*(one-U2) - sqrt(U2*(one-v2)*& + (u*(one-v2*U2) - velx**2*(one-U2))))/(one-v2*U2) + + lamp = lamp_nobeta - beta/alp + lamm = lamm_nobeta - beta/alp + + lam(1) = lamm + lam(2) = lam1 + lam(3) = lam2 + lam(4) = lam3 + lam(5) = lamp + +end subroutine eigenvaluesM_hot !!$ WE need to implement eigenproblem and eigenproblem_leftright!!!! diff --git a/src/GRHydro_HLLEM.F90 b/src/GRHydro_HLLEM.F90 index 3c694ca..87df4d9 100644 --- a/src/GRHydro_HLLEM.F90 +++ b/src/GRHydro_HLLEM.F90 @@ -55,9 +55,10 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS) CCTK_REAL :: pressstarp,pressstarm,velxlowp,velxlowm,velylowp,velylowm,velzlowp,velzlowm CCTK_REAL :: psidcp, psidcm, psidcf, psidcdiff, psidcfp, psidcfm - CCTK_REAL :: charmax_dc, charmin_dc, charpm_dc - + CCTK_REAL :: charmax_dc, charmin_dc, charpm_dc + CCTK_INT :: type_bits, trivial, not_trivial + CCTK_REAL :: xtemp if (flux_direction == 1) then call SpaceMask_GetTypeBits(type_bits, "Hydro_RiemannProblemX") @@ -227,9 +228,9 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS) vxtp,vytp,vztp,pressstarp,bxlowp,bylowp,bzlowp,ab0p,wp, & avg_det,avg_alp,avg_beta) if(clean_divergence.ne.0) then - f1(6)=f1(6) + 1.0*sdet*uxxh*psidcp - cons_p(6)*avg_betax/avg_alp - f1(7)=f1(7) + 1.0*sdet*uxyh*psidcp - cons_p(6)*avg_betay/avg_alp - f1(8)=f1(8) + 1.0*sdet*uxzh*psidcp - cons_p(6)*avg_betaz/avg_alp + f1(6)=f1(6) + 1.0d0*sdet*uxxh*psidcp - cons_p(6)*avg_betax/avg_alp + f1(7)=f1(7) + 1.0d0*sdet*uxyh*psidcp - cons_p(6)*avg_betay/avg_alp + f1(8)=f1(8) + 1.0d0*sdet*uxzh*psidcp - cons_p(6)*avg_betaz/avg_alp psidcf = cons_p(6)/sdet-psidcp*avg_betax/avg_alp endif @@ -240,9 +241,9 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS) vytp,vztp,vxtp,pressstarp,bylowp,bzlowp,bxlowp,ab0p,wp, & avg_det,avg_alp,avg_beta) if(clean_divergence.ne.0) then - f1(6)=f1(6) + 1.0*sdet*uxyh*psidcp - cons_p(7)*avg_betax/avg_alp - f1(7)=f1(7) + 1.0*sdet*uyyh*psidcp - cons_p(7)*avg_betay/avg_alp - f1(8)=f1(8) + 1.0*sdet*uyzh*psidcp - cons_p(7)*avg_betaz/avg_alp + f1(6)=f1(6) + 1.0d0*sdet*uxyh*psidcp - cons_p(7)*avg_betax/avg_alp + f1(7)=f1(7) + 1.0d0*sdet*uyyh*psidcp - cons_p(7)*avg_betay/avg_alp + f1(8)=f1(8) + 1.0d0*sdet*uyzh*psidcp - cons_p(7)*avg_betaz/avg_alp psidcf = cons_p(7)/sdet-psidcp*avg_betay/avg_alp endif @@ -253,9 +254,9 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS) vztp,vxtp,vytp,pressstarp,bzlowp,bxlowp,bylowp,ab0p,wp, & avg_det,avg_alp,avg_beta) if(clean_divergence.ne.0) then - f1(6)=f1(6) + 1.0*sdet*uxzh*psidcp - cons_p(8)*avg_betax/avg_alp - f1(7)=f1(7) + 1.0*sdet*uyzh*psidcp - cons_p(8)*avg_betay/avg_alp - f1(8)=f1(8) + 1.0*sdet*uzzh*psidcp - cons_p(8)*avg_betaz/avg_alp + f1(6)=f1(6) + 1.0d0*sdet*uxzh*psidcp - cons_p(8)*avg_betax/avg_alp + f1(7)=f1(7) + 1.0d0*sdet*uyzh*psidcp - cons_p(8)*avg_betay/avg_alp + f1(8)=f1(8) + 1.0d0*sdet*uzzh*psidcp - cons_p(8)*avg_betaz/avg_alp psidcf = cons_p(8)/sdet-psidcp*avg_betaz/avg_alp endif @@ -293,16 +294,34 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS) !!$ Eigenvalues and fluxes either side of the cell interface if (flux_direction == 1) then - call eigenvaluesM(GRHydro_eos_handle,& - prim_m(1),prim_m(2),prim_m(3),prim_m(4),prim_m(5),prim_m(7), & - prim_m(8),prim_m(9),prim_m(10),& - lamminus,gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,& - usendh,avg_alp,avg_beta) - call eigenvaluesM(GRHydro_eos_handle, & - prim_p(1),prim_p(2),prim_p(3),prim_p(4),prim_p(5),prim_p(7), & - prim_p(8),prim_p(9),prim_p(10),& - lamplus,gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,& - usendh,avg_alp,avg_beta) + if(evolve_temper.ne.1) then + call eigenvaluesM(GRHydro_eos_handle,& + prim_m(1),prim_m(2),prim_m(3),prim_m(4),prim_m(5),prim_m(7), & + prim_m(8),prim_m(9),prim_m(10),& + lamminus,gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,& + usendh,avg_alp,avg_beta) + call eigenvaluesM(GRHydro_eos_handle, & + prim_p(1),prim_p(2),prim_p(3),prim_p(4),prim_p(5),prim_p(7), & + prim_p(8),prim_p(9),prim_p(10),& + lamplus,gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,& + usendh,avg_alp,avg_beta) + else + xtemp = temperature(i+xoffset,j+yoffset,k+zoffset) + call eigenvaluesM_hot(GRHydro_eos_handle,i,j,k,& + prim_m(1),prim_m(2),prim_m(3),prim_m(4),prim_m(5),prim_m(7), & + prim_m(8),prim_m(9),prim_m(10),& + xtemp,y_e_minus(i+xoffset,j+yoffset,k+zoffset),& + lamminus,gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,& + usendh,avg_alp,avg_beta) + xtemp = temperature(i,j,k) + call eigenvaluesM_hot(GRHydro_eos_handle,i,j,k,& + prim_p(1),prim_p(2),prim_p(3),prim_p(4),prim_p(5),prim_p(7), & + prim_p(8),prim_p(9),prim_p(10),& + xtemp,y_e_plus(i,j,k),& + lamplus,gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,& + usendh,avg_alp,avg_beta) + endif + call num_x_fluxM(cons_p(1),cons_p(2),cons_p(3),cons_p(4),cons_p(5),& cons_p(6),cons_p(7),cons_p(8),& fplus(1),fplus(2),fplus(3),fplus(4),fplus(5),fplus(6),fplus(7),fplus(8),& @@ -316,27 +335,45 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS) avg_det,avg_alp,avg_beta) if(clean_divergence.ne.0) then - fminus(6)=fminus(6) + 1.0*sdet*uxxh*psidcm - cons_m(6)*avg_betax/avg_alp - fminus(7)=fminus(7) + 1.0*sdet*uxyh*psidcm - cons_m(6)*avg_betay/avg_alp - fminus(8)=fminus(8) + 1.0*sdet*uxzh*psidcm - cons_m(6)*avg_betaz/avg_alp - fplus(6)=fplus(6) + 1.0*sdet*uxxh*psidcp - cons_p(6)*avg_betax/avg_alp - fplus(7)=fplus(7) + 1.0*sdet*uxyh*psidcp - cons_p(6)*avg_betay/avg_alp - fplus(8)=fplus(8) + 1.0*sdet*uxzh*psidcp - cons_p(6)*avg_betaz/avg_alp + fminus(6)=fminus(6) + 1.0d0*sdet*uxxh*psidcm - cons_m(6)*avg_betax/avg_alp + fminus(7)=fminus(7) + 1.0d0*sdet*uxyh*psidcm - cons_m(6)*avg_betay/avg_alp + fminus(8)=fminus(8) + 1.0d0*sdet*uxzh*psidcm - cons_m(6)*avg_betaz/avg_alp + fplus(6)=fplus(6) + 1.0d0*sdet*uxxh*psidcp - cons_p(6)*avg_betax/avg_alp + fplus(7)=fplus(7) + 1.0d0*sdet*uxyh*psidcp - cons_p(6)*avg_betay/avg_alp + fplus(8)=fplus(8) + 1.0d0*sdet*uxzh*psidcp - cons_p(6)*avg_betaz/avg_alp psidcfp = cons_p(6)/sdet-avg_betax*psidcp/avg_alp psidcfm = cons_m(6)/sdet-avg_betax*psidcm/avg_alp endif else if (flux_direction == 2) then - call eigenvaluesM(GRHydro_eos_handle,& - prim_m(1),prim_m(3),prim_m(4),prim_m(2),prim_m(5),prim_m(7), & - prim_m(9),prim_m(10),prim_m(8),& - lamminus,gyyh,gyzh,gxyh,gzzh,gxzh,gxxh,& - usendh,avg_alp,avg_beta) - call eigenvaluesM(GRHydro_eos_handle, & - prim_p(1),prim_p(3),prim_p(4),prim_p(2),prim_p(5),prim_p(7), & - prim_p(9),prim_p(10),prim_p(8),& - lamplus,gyyh,gyzh,gxyh,gzzh,gxzh,gxxh,& - usendh,avg_alp,avg_beta) + if(evolve_temper.ne.1) then + call eigenvaluesM(GRHydro_eos_handle,& + prim_m(1),prim_m(3),prim_m(4),prim_m(2),prim_m(5),prim_m(7), & + prim_m(9),prim_m(10),prim_m(8),& + lamminus,gyyh,gyzh,gxyh,gzzh,gxzh,gxxh,& + usendh,avg_alp,avg_beta) + call eigenvaluesM(GRHydro_eos_handle, & + prim_p(1),prim_p(3),prim_p(4),prim_p(2),prim_p(5),prim_p(7), & + prim_p(9),prim_p(10),prim_p(8),& + lamplus,gyyh,gyzh,gxyh,gzzh,gxzh,gxxh,& + usendh,avg_alp,avg_beta) + else + xtemp = temperature(i+xoffset,j+yoffset,k+zoffset) + call eigenvaluesM_hot(GRHydro_eos_handle,i,j,k,& + prim_m(1),prim_m(3),prim_m(4),prim_m(2),prim_m(5),prim_m(7), & + prim_m(9),prim_m(10),prim_m(8),& + xtemp,y_e_minus(i+xoffset,j+yoffset,k+zoffset),& + lamminus,gyyh,gyzh,gxyh,gzzh,gxzh,gxxh,& + usendh,avg_alp,avg_beta) + xtemp = temperature(i,j,k) + call eigenvaluesM_hot(GRHydro_eos_handle,i,j,k,& + prim_p(1),prim_p(3),prim_p(4),prim_p(2),prim_p(5),prim_p(7), & + prim_p(9),prim_p(10),prim_p(8),& + xtemp,y_e_plus(i,j,k),& + lamplus,gyyh,gyzh,gxyh,gzzh,gxzh,gxxh,& + usendh,avg_alp,avg_beta) + endif + call num_x_fluxM(cons_p(1),cons_p(3),cons_p(4),cons_p(2),cons_p(5),& cons_p(7),cons_p(8),cons_p(6),& fplus(1),fplus(3),fplus(4),fplus(2),fplus(5),fplus(7),fplus(8),fplus(6),& @@ -350,27 +387,45 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS) avg_det,avg_alp,avg_beta) if(clean_divergence.ne.0) then - fminus(6)=fminus(6) + 1.0*sdet*uxyh*psidcm - cons_m(7)*avg_betax/avg_alp - fminus(7)=fminus(7) + 1.0*sdet*uyyh*psidcm - cons_m(7)*avg_betay/avg_alp - fminus(8)=fminus(8) + 1.0*sdet*uyzh*psidcm - cons_m(7)*avg_betaz/avg_alp - fplus(6)=fplus(6) + 1.0*sdet*uxyh*psidcp - cons_p(7)*avg_betax/avg_alp - fplus(7)=fplus(7) + 1.0*sdet*uyyh*psidcp - cons_p(7)*avg_betay/avg_alp - fplus(8)=fplus(8) + 1.0*sdet*uyzh*psidcp - cons_p(7)*avg_betaz/avg_alp + fminus(6)=fminus(6) + 1.0d0*sdet*uxyh*psidcm - cons_m(7)*avg_betax/avg_alp + fminus(7)=fminus(7) + 1.0d0*sdet*uyyh*psidcm - cons_m(7)*avg_betay/avg_alp + fminus(8)=fminus(8) + 1.0d0*sdet*uyzh*psidcm - cons_m(7)*avg_betaz/avg_alp + fplus(6)=fplus(6) + 1.0d0*sdet*uxyh*psidcp - cons_p(7)*avg_betax/avg_alp + fplus(7)=fplus(7) + 1.0d0*sdet*uyyh*psidcp - cons_p(7)*avg_betay/avg_alp + fplus(8)=fplus(8) + 1.0d0*sdet*uyzh*psidcp - cons_p(7)*avg_betaz/avg_alp psidcfp = cons_p(7)/sdet-avg_betay*psidcp/avg_alp psidcfm = cons_m(7)/sdet-avg_betay*psidcm/avg_alp endif else if (flux_direction == 3) then - call eigenvaluesM(GRHydro_eos_handle,& - prim_m(1),prim_m(4),prim_m(2),prim_m(3),prim_m(5),prim_m(7), & - prim_m(10),prim_m(8),prim_m(9),& - lamminus,gzzh,gxzh,gyzh,gxxh,gxyh,gyyh,& - usendh,avg_alp,avg_beta) - call eigenvaluesM(GRHydro_eos_handle,& - prim_p(1),prim_p(4),prim_p(2),prim_p(3),prim_p(5),prim_p(7), & - prim_p(10),prim_p(8),prim_p(9),& - lamplus,gzzh,gxzh,gyzh,gxxh,gxyh,gyyh,& - usendh,avg_alp,avg_beta) + if(evolve_temper.ne.1) then + call eigenvaluesM(GRHydro_eos_handle,& + prim_m(1),prim_m(4),prim_m(2),prim_m(3),prim_m(5),prim_m(7), & + prim_m(10),prim_m(8),prim_m(9),& + lamminus,gzzh,gxzh,gyzh,gxxh,gxyh,gyyh,& + usendh,avg_alp,avg_beta) + call eigenvaluesM(GRHydro_eos_handle, & + prim_p(1),prim_p(4),prim_p(2),prim_p(3),prim_p(5),prim_p(7), & + prim_p(10),prim_p(8),prim_p(9),& + lamplus,gzzh,gxzh,gyzh,gxxh,gxyh,gyyh,& + usendh,avg_alp,avg_beta) + else + xtemp = temperature(i+xoffset,j+yoffset,k+zoffset) + call eigenvaluesM_hot(GRHydro_eos_handle,i,j,k,& + prim_m(1),prim_m(4),prim_m(2),prim_m(3),prim_m(5),prim_m(7), & + prim_m(10),prim_m(8),prim_m(9),& + xtemp,y_e_minus(i+xoffset,j+yoffset,k+zoffset),& + lamminus,gzzh,gxzh,gyzh,gxxh,gxyh,gyyh,& + usendh,avg_alp,avg_beta) + xtemp = temperature(i,j,k) + call eigenvaluesM_hot(GRHydro_eos_handle,i,j,k,& + prim_p(1),prim_p(4),prim_p(2),prim_p(3),prim_p(5),prim_p(7), & + prim_p(10),prim_p(8),prim_p(9),& + xtemp,y_e_plus(i,j,k),& + lamplus,gzzh,gxzh,gyzh,gxxh,gxyh,gyyh,& + usendh,avg_alp,avg_beta) + endif + call num_x_fluxM(cons_p(1),cons_p(4),cons_p(2),cons_p(3),cons_p(5),& cons_p(8),cons_p(6),cons_p(7),& fplus(1),fplus(4),fplus(2),fplus(3),fplus(5),fplus(8),fplus(6),fplus(7), & @@ -384,12 +439,12 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS) avg_det,avg_alp,avg_beta) if(clean_divergence.ne.0) then - fminus(6)=fminus(6) + 1.0*sdet*uxzh*psidcm - cons_m(8)*avg_betax/avg_alp - fminus(7)=fminus(7) + 1.0*sdet*uyzh*psidcm - cons_m(8)*avg_betay/avg_alp - fminus(8)=fminus(8) + 1.0*sdet*uzzh*psidcm - cons_m(8)*avg_betaz/avg_alp - fplus(6)=fplus(6) + 1.0*sdet*uxzh*psidcp - cons_p(8)*avg_betax/avg_alp - fplus(7)=fplus(7) + 1.0*sdet*uyzh*psidcp - cons_p(8)*avg_betay/avg_alp - fplus(8)=fplus(8) + 1.0*sdet*uzzh*psidcp - cons_p(8)*avg_betaz/avg_alp + fminus(6)=fminus(6) + 1.0d0*sdet*uxzh*psidcm - cons_m(8)*avg_betax/avg_alp + fminus(7)=fminus(7) + 1.0d0*sdet*uyzh*psidcm - cons_m(8)*avg_betay/avg_alp + fminus(8)=fminus(8) + 1.0d0*sdet*uzzh*psidcm - cons_m(8)*avg_betaz/avg_alp + fplus(6)=fplus(6) + 1.0d0*sdet*uxzh*psidcp - cons_p(8)*avg_betax/avg_alp + fplus(7)=fplus(7) + 1.0d0*sdet*uyzh*psidcp - cons_p(8)*avg_betay/avg_alp + fplus(8)=fplus(8) + 1.0d0*sdet*uzzh*psidcp - cons_p(8)*avg_betaz/avg_alp psidcfp = cons_p(8)/sdet-avg_betaz*psidcp/avg_alp psidcfm = cons_m(8)/sdet-avg_betaz*psidcm/avg_alp endif diff --git a/src/GRHydro_InterfacesM.h b/src/GRHydro_InterfacesM.h index 8cd85fa..5463cc7 100644 --- a/src/GRHydro_InterfacesM.h +++ b/src/GRHydro_InterfacesM.h @@ -3,11 +3,12 @@ module Con2PrimM_fortran_interfaces interface - subroutine GRHydro_Con2PrimM_pt( handle, & + subroutine GRHydro_Con2PrimM_pt( handle, keytemp, prec, & local_gam, dens, & sx, sy, sz, & tau, & Bconsx, Bconsy, Bconsz, & + xye, xtemp, & rho, & velx, vely, velz,& epsilon, pressure,& @@ -24,11 +25,15 @@ module Con2PrimM_fortran_interfaces implicit none CCTK_INT handle + CCTK_INT keytemp + CCTK_REAL prec CCTK_REAL local_gam CCTK_REAL dens CCTK_REAL sx, sy, sz CCTK_REAL tau CCTK_REAL Bconsx, Bconsy, Bconsz + CCTK_REAL xye + CCTK_REAL xtemp CCTK_REAL rho CCTK_REAL velx, vely, velz CCTK_REAL epsilon, pressure diff --git a/src/GRHydro_PPMMReconstruct_drv.F90 b/src/GRHydro_PPMMReconstruct_drv.F90 index 40b1c93..93538b3 100644 --- a/src/GRHydro_PPMMReconstruct_drv.F90 +++ b/src/GRHydro_PPMMReconstruct_drv.F90 @@ -243,8 +243,8 @@ subroutine GRHydro_PPMMReconstruct_drv(CCTK_ARGUMENTS) if(evolve_Y_e.ne.0) then !$OMP PARALLEL DO PRIVATE(j, k) - do k = GRHydro_stencil, nz - GRHydro_stencil + 1 - do j = GRHydro_stencil, ny - GRHydro_stencil + 1 + do k = GRHydro_stencil, nz - GRHydro_stencil + 1 + transport_constraints + do j = GRHydro_stencil, ny - GRHydro_stencil + 1 + transport_constraints call SimplePPM_ye_1d(nx,CCTK_DELTA_SPACE(1),rho(:,j,k), & velx(:,j,k),vely(:,j,k),velz(:,j,k), & Y_e(:,j,k),Y_e_minus(:,j,k),Y_e_plus(:,j,k), & @@ -315,8 +315,8 @@ subroutine GRHydro_PPMMReconstruct_drv(CCTK_ARGUMENTS) end if if(evolve_Y_e.ne.0) then !$OMP PARALLEL DO PRIVATE(j, k) - do k = GRHydro_stencil, nz - GRHydro_stencil + 1 - do j = GRHydro_stencil, nx - GRHydro_stencil + 1 + do k = GRHydro_stencil, nz - GRHydro_stencil + 1 + transport_constraints + do j = GRHydro_stencil, nx - GRHydro_stencil + 1 + transport_constraints call SimplePPM_ye_1d(ny,CCTK_DELTA_SPACE(2),rho(j,:,k), & velx(j,:,k),vely(j,:,k),velz(j,:,k), & Y_e(j,:,k),Y_e_minus(j,:,k),Y_e_plus(j,:,k), & @@ -388,8 +388,8 @@ subroutine GRHydro_PPMMReconstruct_drv(CCTK_ARGUMENTS) end if if(evolve_Y_e.ne.0) then !$OMP PARALLEL DO PRIVATE(j, k) - do k = GRHydro_stencil, ny - GRHydro_stencil + 1 - do j = GRHydro_stencil, nx - GRHydro_stencil + 1 + do k = GRHydro_stencil, ny - GRHydro_stencil + 1 + transport_constraints + do j = GRHydro_stencil, nx - GRHydro_stencil + 1 + transport_constraints call SimplePPM_ye_1d(nz,CCTK_DELTA_SPACE(3),rho(j,k,:), & velx(j,k,:),vely(j,k,:),velz(j,k,:), & Y_e(j,k,:),Y_e_minus(j,k,:),Y_e_plus(j,k,:), & diff --git a/src/GRHydro_Prim2ConM.F90 b/src/GRHydro_Prim2ConM.F90 index 59acf0d..a6aef7d 100644 --- a/src/GRHydro_Prim2ConM.F90 +++ b/src/GRHydro_Prim2ConM.F90 @@ -56,9 +56,29 @@ subroutine primitive2conservativeM(CCTK_ARGUMENTS) integer :: i, j, k CCTK_REAL :: gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_detl,& gxxr,gxyr,gxzr,gyyr,gyzr,gzzr,avg_detr + !CCTK_REAL :: g11l,g12l,g13l,g22l,g23l,g33l,& + ! g11r,g12r,g13r,g22r,g23r,g33r + CCTK_REAL :: xtemp(1) + character(len=256) NaN_WarnLine + !TARGET gxx, gxy, gxz, gyy, gyz, gzz + + !CCTK_REAL, DIMENSION(:,:,:), POINTER :: g11, g12, g13, g22, g23, g33 + !g11 => gxx + !g12 => gxy + !g13 => gxz + !g22 => gyy + !g23 => gyz + !g33 => gzz + ! constraint transport needs to be able to average fluxes in the directions ! other that flux_direction, which in turn need the primitives on interfaces + + if(evolve_temper.ne.1) then + + !$OMP PARALLEL DO PRIVATE(k,j,i,avg_detl,avg_detr,& + !$OMP gxxl,gxyl,gxzl,gyyl,gyzl,gzzl, & + !$OMP gxxr,gxyr,gxzr,gyyr,gyzr,gzzr) do k = GRHydro_stencil, cctk_lsh(3)-GRHydro_stencil+1 + transport_constraints*(1-zoffset) do j = GRHydro_stencil, cctk_lsh(2)-GRHydro_stencil+1 + transport_constraints*(1-yoffset) do i = GRHydro_stencil, cctk_lsh(1)-GRHydro_stencil+1 + transport_constraints*(1-xoffset) @@ -99,6 +119,79 @@ subroutine primitive2conservativeM(CCTK_ARGUMENTS) end do end do end do + !$OMP END PARALLEL DO + + else + + !$OMP PARALLEL DO PRIVATE(i, j, k, avg_detl, avg_detr, xtemp,& + !$OMP gxxl,gxyl,gxzl,gyyl,gyzl,gzzl, & + !$OMP gxxr,gxyr,gxzr,gyyr,gyzr,gzzr) + do k = GRHydro_stencil, cctk_lsh(3)-GRHydro_stencil+1 + transport_constraints*(1-zoffset) + do j = GRHydro_stencil, cctk_lsh(2)-GRHydro_stencil+1 + transport_constraints*(1-yoffset) + do i = GRHydro_stencil, cctk_lsh(1)-GRHydro_stencil+1 + transport_constraints*(1-xoffset) + + 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)) + + avg_detl = SPATIAL_DETERMINANT(gxxl,gxyl,gxzl,gyyl,gyzl,gzzl) + avg_detr = SPATIAL_DETERMINANT(gxxr,gxyr,gxzr,gyyr,gyzr,gzzr) + + ! we do not have plus/minus vars for temperature since + ! eps is reconstructed. Hence, we do not update the + ! variable 'temperature' in prim2con at the interfaces + ! We will instead use an average temperature as an initial + ! guess. + xtemp(1) = 0.5d0*(temperature(i,j,k) + & + temperature(i-xoffset,j-yoffset,k-zoffset)) + + !$OMP CRITICAL + if (y_e_minus(i,j,k) .le. 0.0d0 .or. y_e_plus(i,j,k) .le. 0.0d0) then + write(NaN_WarnLine,'(a100,7g15.6)') '(y_e_minus,y_e_plus,x,y,z,rho)', y_e(i,j,k), y_e_minus(i,j,k), y_e_plus(i,j,k), x(i,j,k),y(i,j,k),z(i,j,k),rho(i,j,k) + call CCTK_WARN(1, NaN_WarnLine) + endif + !$OMP END CRITICAL + call prim2conM_hot(GRHydro_eos_handle, GRHydro_reflevel,& + i,j,k,x(i,j,k),y(i,j,k),z(i,j,k), gxxl,gxyl,gxzl,gyyl,gyzl,gzzl, & + avg_detl,densminus(i,j,k),sxminus(i,j,k),& + syminus(i,j,k),szminus(i,j,k),tauminus(i,j,k),& + Bconsxminus(i,j,k),Bconsyminus(i,j,k),Bconszminus(i,j,k), rhominus(i,j,k), & + velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k),& + epsminus(i,j,k),pressminus(i,j,k),Bvecxminus(i,j,k), & + Bvecyminus(i,j,k), Bveczminus(i,j,k), w_lorentzminus(i, j, k), xtemp, y_e_minus(i,j,k)) + ! we do not have plus/minus vars for temperature since + ! eps is reconstructed. Hence, we do not update the + ! variable 'temperature' in prim2con at the interfaces + ! We will instead use an average temperature as an initial + ! guess. + xtemp(1) = 0.5d0*(temperature(i,j,k) + & + temperature(i+xoffset,j+yoffset,k+zoffset)) + + call prim2conM_hot(GRHydro_eos_handle, GRHydro_reflevel,& + i,j,k,x(i,j,k),y(i,j,k),z(i,j,k), gxxr,gxyr,gxzr,gyyr,gyzr,gzzr, & + avg_detr, densplus(i,j,k),sxplus(i,j,k),& + syplus(i,j,k),szplus(i,j,k),tauplus(i,j,k),& + Bconsxplus(i,j,k),Bconsyplus(i,j,k),Bconszplus(i,j,k),& + rhoplus(i,j,k),velxplus(i,j,k),velyplus(i,j,k),& + velzplus(i,j,k),epsplus(i,j,k),pressplus(i,j,k),& + Bvecxplus(i,j,k), Bvecyplus(i,j,k), Bveczplus(i,j,k), & + w_lorentzplus(i,j,k), xtemp, y_e_plus(i,j,k)) + + end do + end do + end do + !$OMP END PARALLEL DO + endif + end subroutine primitive2conservativeM @@ -117,6 +210,152 @@ end subroutine primitive2conservativeM @@*/ +subroutine prim2conM_hot(handle, GRHydro_reflevel, ii, jj, kk, x, y, z, gxx, gxy, gxz, gyy, gyz, gzz, det, ddens, & + dsx, dsy, dsz, dtau , dBconsx, dBconsy, dBconsz, drho, dvelx, dvely, dvelz, deps, dpress, dBvcx, dBvcy, dBvcz, w, temp, ye) + + implicit none + + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_FUNCTIONS + + CCTK_REAL :: gxx, gxy, gxz, gyy, gyz, gzz, det + CCTK_REAL, dimension(1) :: ddens, dsx, dsy, dsz, dtau, & + dBconsx, dBconsy, dBconsz, & + drho, dvelx, dvely, dvelz,& + deps, dpress, dBvcx, dBvcy, dBvcz, vlowx, vlowy, vlowz + CCTK_REAL :: temp(1),ye(1), yein(1), x, y, z + CCTK_INT :: handle, GRHydro_reflevel, ii, jj, kk + CCTK_REAL :: w + CCTK_REAL, dimension(1) :: Bdotv,ab0,b2,blowx,blowy,blowz + character(len=256) NaN_WarnLine + character(len=512) warnline + +! begin EOS Omni vars + integer :: n, keytemp, anyerr, keyerr(1) + ! real*8 :: xpress(1),xeps(1),xtemp(1),xye(1) + real*8 :: temp0(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 + + temp0 = temp + w = 1.d0 / sqrt(1.d0 - DOT2(dvelx(1),dvely(1),dvelz(1))) + +!!$ BEGIN: Check for NaN value + if (w .ne. w) then + !$OMP CRITICAL + write(NaN_WarnLine,'(a100,6g15.6)') 'NaN produced in sqrt(): (gxx,gxy,gxz,gyy,gyz,gzz)', gxx, gxy, gxz, gyy, gyz, gzz + call CCTK_WARN(1, NaN_WarnLine) + write(NaN_WarnLine,'(a100,3g15.6)') 'NaN produced in sqrt(): (dvelx,dvely,dvelz)', dvelx, dvely, dvelz + call CCTK_WARN(1, NaN_WarnLine) + write(NaN_WarnLine,'(a100,3g15.6)') 'NaN produced in sqrt(): (x,y,z)', x, y, z + call CCTK_WARN(1, NaN_WarnLine) + !write(NaN_WarnLine,'(a100,g15.6)') 'NaN produced in sqrt(): w', w + !call CCTK_WARN(GRHydro_NaN_verbose, NaN_WarnLine) + !$OMP END CRITICAL + endif +!!$ END: Check for NaN value + + ye = max(min(ye,GRHydro_Y_e_max),GRHydro_Y_e_min) + + yein = ye + + call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,& + drho,deps,temp,ye,dpress,keyerr,anyerr) + ! error handling + if(anyerr.ne.0) then + if(GRHydro_reflevel.lt.GRHydro_c2p_warn_from_reflevel) then + ! in this case (coarse grid error that is hopefully restricted + ! away), we use the average temperature between cells and call + ! the EOS with keytemp=1 + keytemp=1 + call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,& + drho,deps,temp,ye,dpress,keyerr,anyerr) + keytemp=0 + if(anyerr.ne.0) then + !$OMP CRITICAL + call CCTK_WARN(1,"EOS error in prim2con_hot: lev 2") + write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z + call CCTK_WARN(1,warnline) + write(warnline,"(1P10E15.6)") drho,deps,temp,ye,yein + 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) + !$OMP END CRITICAL + endif + else + ! This is a way of recovering even on finer refinement levels: + ! Use the average temperature at the interface instead of the + ! reconstructed specific internal energy. + !$OMP CRITICAL + call CCTK_WARN(1,"EOS error in prim2con_hot: NOW using averaged temp!") + write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z + call CCTK_WARN(1,warnline) + write(warnline,"(1P10E15.6)") drho,deps,temp0,ye + 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) + !$OMP END CRITICAL + keytemp=1 + temp = temp0 + call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,& + drho,deps,temp,ye,dpress,keyerr,anyerr) + keytemp=0 + if(anyerr.ne.0) then + !$OMP CRITICAL + call CCTK_WARN(1,"EOS error in prim2con_hot") + write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z + call CCTK_WARN(1,warnline) + write(warnline,"(1P10E15.6)") drho,deps,temp,ye + 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 + endif + endif + + vlowx = gxx*dvelx + gxy*dvely + gxz*dvelz + vlowy = gxy*dvelx + gyy*dvely + gyz*dvelz + vlowz = gxz*dvelx + gyz*dvely + gzz*dvelz + + Bdotv=DOT(dvelx,dvely,dvelz,dBvcx,dBvcy,dBvcz) + ab0=w*Bdotv + b2 = DOT2(dBvcx,dBvcy,dBvcz)/w**2+Bdotv**2 + blowx = (gxx*dBvcx + gxy*dBvcy + gxz*dBvcz)/w + & + w*Bdotv*vlowx + blowy = (gxy*dBvcx + gyy*dBvcy + gyz*dBvcz)/w + & + w*Bdotv*vlowy + blowz = (gxz*dBvcx + gyz*dBvcy + gzz*dBvcz)/w + & + w*Bdotv*vlowz + + ddens = sqrt(det) * drho * w + dsx = sqrt(det) * ((drho*(1+deps)+dpress+b2)*w*w * vlowx - & + ab0*blowx) + dsy = sqrt(det) * ((drho*(1+deps)+dpress+b2)*w*w * vlowy - & + ab0*blowy) + dsz = sqrt(det) * ((drho*(1+deps)+dpress+b2)*w*w * vlowz - & + ab0*blowz) + dtau = sqrt(det) * ((drho*(1+deps)+dpress+b2)*w*w - dpress-b2/2.0-ab0**2) - ddens + + dBconsx = sqrt(det)*dBvcx + dBconsy = sqrt(det)*dBvcy + dBconsz = sqrt(det)*dBvcz + + !!$OMP CRITICAL + !write(NaN_WarnLine,'(a100,3g15.6)') '(dens out, sqrt(det), w:)', ddens,sqrt(det),w + !call CCTK_WARN(1, NaN_WarnLine) + !!$OMP END CRITICAL +end subroutine prim2conM_hot + + subroutine prim2conM(handle, gxx, gxy, gxz, gyy, gyz, gzz, det, ddens, & dsx, dsy, dsz, dtau , dBconsx, dBconsy, dBconsz, drho, dvelx, dvely, dvelz, deps, dpress, dBvcx, dBvcy, dBvcz, w) @@ -202,7 +441,6 @@ end subroutine prim2conM @@*/ - subroutine Primitive2ConservativeCellsM(CCTK_ARGUMENTS) implicit none @@ -210,9 +448,13 @@ subroutine Primitive2ConservativeCellsM(CCTK_ARGUMENTS) DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS + CCTK_REAL :: xtemp(1) CCTK_INT :: i, j, k CCTK_REAL :: det + if(evolve_temper.ne.1) then + + !$OMP PARALLEL DO PRIVATE(k,j,i,det) do k = GRHydro_stencil,cctk_lsh(3)-GRHydro_stencil+1 do j = GRHydro_stencil,cctk_lsh(2)-GRHydro_stencil+1 do i = GRHydro_stencil,cctk_lsh(1)-GRHydro_stencil+1 @@ -232,6 +474,31 @@ subroutine Primitive2ConservativeCellsM(CCTK_ARGUMENTS) end do end do end do + !$OMP END PARALLEL DO + else + !$OMP PARALLEL DO PRIVATE(k,j,i,det) + do k = GRHydro_stencil,cctk_lsh(3)-GRHydro_stencil+1 + do j = GRHydro_stencil,cctk_lsh(2)-GRHydro_stencil+1 + do i = GRHydro_stencil,cctk_lsh(1)-GRHydro_stencil+1 + + det = SPATIAL_DETERMINANT(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k), \ + gyy(i,j,k),gyz(i,j,k),gzz(i,j,k)) + + call prim2conM_hot(GRHydro_eos_handle,GRHydro_reflevel,& + i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),& + gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),& + gyy(i,j,k),gyz(i,j,k),gzz(i,j,k),& + det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),& + tau(i,j,k),Bconsx(i,j,k),Bconsy(i,j,k),Bconsz(i,j,k),& + rho(i,j,k),velx(i,j,k),vely(i,j,k),velz(i,j,k),& + eps(i,j,k),press(i,j,k),Bvecx(i,j,k), & + Bvecy(i,j,k), Bvecz(i,j,k), w_lorentz(i,j,k), temperature(i,j,k), Y_e(i,j,k)) + + end do + end do + end do + !$OMP END PARALLEL DO + endif if(evolve_Y_e.ne.0) then !$OMP PARALLEL DO PRIVATE(i, j, k) @@ -451,6 +718,7 @@ subroutine Primitive2ConservativePolyCellsM(CCTK_ARGUMENTS) CCTK_INT :: i, j, k CCTK_REAL :: det + !$OMP PARALLEL DO PRIVATE(k,j,i,det) do k = GRHydro_stencil,cctk_lsh(3)-GRHydro_stencil+1 do j = GRHydro_stencil,cctk_lsh(2)-GRHydro_stencil+1 do i = GRHydro_stencil,cctk_lsh(1)-GRHydro_stencil+1 @@ -469,6 +737,7 @@ subroutine Primitive2ConservativePolyCellsM(CCTK_ARGUMENTS) end do end do end do + !$OMP END PARALLEL DO if(evolve_Y_e.ne.0) then !$OMP PARALLEL DO PRIVATE(i, j, k) diff --git a/src/GRHydro_UpdateMaskM.F90 b/src/GRHydro_UpdateMaskM.F90 index d6a6bb8..8c3feef 100644 --- a/src/GRHydro_UpdateMaskM.F90 +++ b/src/GRHydro_UpdateMaskM.F90 @@ -44,11 +44,18 @@ subroutine GRHydro_AtmosphereResetM(CCTK_ARGUMENTS) CCTK_INT :: i, j, k CCTK_REAL :: det, psi4pt +! begin EOS Omni vars + integer :: n,keytemp,anyerr,keyerr(1) + real*8 :: xpress,xeps,xtemp,xye + n=1;keytemp=0;anyerr=0;keyerr(1)=0 + xpress=0.0d0;xye=0.0d0;xeps=0.0d0;xtemp=0.0d0 +! end EOS Omni vars + do k = 1, cctk_lsh(3) do j = 1, cctk_lsh(2) do i = 1, cctk_lsh(1) - - if ( atmosphere_mask(i, j, k) .eq. 1) then + + if (atmosphere_mask(i, j, k) .ne. 0) then rho(i,j,k) = GRHydro_rho_min vel(i,j,k,1) = 0.0d0 @@ -56,20 +63,43 @@ subroutine GRHydro_AtmosphereResetM(CCTK_ARGUMENTS) vel(i,j,k,3) = 0.0d0 det = SPATIAL_DETERMINANT(gxx(i,j,k), gxy(i,j,k), gxz(i,j,k), \ gyy(i,j,k), gyz(i,j,k), gzz(i,j,k)) - call prim2conpolytypeM(GRHydro_polytrope_handle, & - gxx(i,j,k), gxy(i,j,k), gxz(i,j,k), & - gyy(i,j,k), gyz(i,j,k), gzz(i,j,k), det, & - dens(i,j,k), 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), 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),w_lorentz(i,j,k)) - if (wk_atmosphere .eq. 0) then - atmosphere_mask(i, j, k) = 0 - atmosphere_mask_real(i, j, k) = 0 - end if - end if + if(evolve_temper.ne.0) then + +! ! set the temperature to be relatively low + temperature(i,j,k) = grhydro_hot_atmo_temp + y_e(i,j,k) = grhydro_hot_atmo_Y_e + 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) + call prim2conM_hot(GRHydro_eos_handle, GRHydro_reflevel,& + i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),gxx(i,j,k),& + gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k),& + det, dens(i,j,k),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),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), w_lorentz(i,j,k),& + temperature(i,j,k),y_e(i,j,k)) + y_e_con(i,j,k) = dens(i,j,k) * y_e(i,j,k) + + else + + call prim2conpolytypeM(GRHydro_polytrope_handle, & + gxx(i,j,k), gxy(i,j,k), gxz(i,j,k), & + gyy(i,j,k), gyz(i,j,k), gzz(i,j,k), det, & + dens(i,j,k), 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), 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),w_lorentz(i,j,k)) + if (wk_atmosphere .eq. 0) then + atmosphere_mask(i, j, k) = 0 + end if + + end if + endif end do end do @@ -88,6 +118,7 @@ subroutine GRHydro_InitialAtmosphereResetM(CCTK_ARGUMENTS) CCTK_INT :: i, j, k CCTK_REAL :: det, psi4pt + CCTK_REAL :: rho_min CCTK_INT :: eos_handle @@ -100,23 +131,50 @@ subroutine GRHydro_InitialAtmosphereResetM(CCTK_ARGUMENTS) eos_handle = GRHydro_polytrope_handle + rho_min = GRHydro_rho_min + if (initial_atmosphere_factor .gt. 0) then + rho_min = rho_min * initial_atmosphere_factor + endif + do k = 1, cctk_lsh(3) do j = 1, cctk_lsh(2) do i = 1, cctk_lsh(1) - if (rho(i,j,k) .le. GRHydro_rho_min) then - rho(i,j,k) = GRHydro_rho_min + if (rho(i,j,k) .le. rho_min) then + rho(i,j,k) = rho_min vel(i,j,k,1) = 0.0d0 vel(i,j,k,2) = 0.0d0 vel(i,j,k,3) = 0.0d0 + det = SPATIAL_DETERMINANT(gxx(i,j,k), gxy(i,j,k), gxz(i,j,k), \ + gyy(i,j,k), gyz(i,j,k), gzz(i,j,k)) + + if(evolve_temper.ne.0) then +! ! set the temperature to be relatively low + temperature(i,j,k) = grhydro_hot_atmo_temp + y_e(i,j,k) = grhydro_hot_atmo_Y_e + 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) + call prim2conM_hot(GRHydro_eos_handle, GRHydro_reflevel,& + i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),gxx(i,j,k),& + gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k),& + det, dens(i,j,k),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),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), w_lorentz(i,j,k),& + temperature(i,j,k),y_e(i,j,k)) + y_e_con(i,j,k) = dens(i,j,k) * y_e(i,j,k) + + else + call EOS_Omni_press(eos_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(eos_handle,keytemp,GRHydro_eos_rf_prec,n,& rho(i,j,k),xeps,xtemp,xye,press(i,j,k),eps(i,j,k),keyerr,anyerr) - det = SPATIAL_DETERMINANT(gxx(i,j,k), gxy(i,j,k), gxz(i,j,k), \ - gyy(i,j,k), gyz(i,j,k), gzz(i,j,k)) call prim2conpolytypeM(eos_handle, & gxx(i,j,k), gxy(i,j,k), gxz(i,j,k), & gyy(i,j,k), gyz(i,j,k), gzz(i,j,k), det, & @@ -125,6 +183,7 @@ subroutine GRHydro_InitialAtmosphereResetM(CCTK_ARGUMENTS) 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),w_lorentz(i,j,k)) + end if end if if (timelevels .gt. 1) then if (rho_p(i,j,k) .le. GRHydro_rho_min) then @@ -132,47 +191,93 @@ subroutine GRHydro_InitialAtmosphereResetM(CCTK_ARGUMENTS) vel_p(i,j,k,1) = 0.0d0 vel_p(i,j,k,2) = 0.0d0 vel_p(i,j,k,3) = 0.0d0 + + det = SPATIAL_DETERMINANT(gxx_p(i,j,k), gxy_p(i,j,k), gxz_p(i,j,k), \ + gyy_p(i,j,k), gyz_p(i,j,k), gzz_p(i,j,k)) + if(evolve_temper.ne.0) then +! ! set the temperature to be relatively low + temperature_p(i,j,k) = grhydro_hot_atmo_temp + y_e_p(i,j,k) = grhydro_hot_atmo_Y_e + keytemp = 1 + call EOS_Omni_press(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,& + rho_p(i,j,k),eps_p(i,j,k),temperature_p(i,j,k),y_e_p(i,j,k),& + press_p(i,j,k),keyerr,anyerr) + call prim2conM_hot(GRHydro_eos_handle, GRHydro_reflevel,& + i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),gxx_p(i,j,k),& + gxy_p(i,j,k),gxz_p(i,j,k),gyy_p(i,j,k),gyz_p(i,j,k),gzz_p(i,j,k),& + det, dens_p(i,j,k),scon_p(i,j,k,1),scon_p(i,j,k,2),scon_p(i,j,k,3),& + tau_p(i,j,k),Bcons_p(i,j,k,1),Bcons_p(i,j,k,2),Bcons_p(i,j,k,3),& + rho_p(i,j,k),vel_p(i,j,k,1),vel_p(i,j,k,2),vel_p(i,j,k,3),& + eps_p(i,j,k),press_p(i,j,k),Bvec_p(i,j,k,1), & + Bvec_p(i,j,k,2), Bvec_p(i,j,k,3), w_lorentz_p(i,j,k),& + temperature_p(i,j,k),y_e_p(i,j,k)) + y_e_con_p(i,j,k) = dens_p(i,j,k) * y_e_p(i,j,k) - call EOS_Omni_press(eos_handle,keytemp,GRHydro_eos_rf_prec,n,& - rho_p(i,j,k),eps_p(i,j,k),xtemp,xye,press_p(i,j,k),keyerr,anyerr) - call EOS_Omni_EpsFromPress(eos_handle,keytemp,GRHydro_eos_rf_prec,n,& - rho_p(i,j,k),xeps,xtemp,xye,press_p(i,j,k),eps_p(i,j,k),keyerr,anyerr) + else - det = SPATIAL_DETERMINANT(gxx_p(i,j,k), gxy_p(i,j,k), gxz_p(i,j,k), \ - gyy_p(i,j,k), gyz_p(i,j,k), gzz_p(i,j,k)) - call prim2conpolytypeM(eos_handle, & - gxx_p(i,j,k), gxy_p(i,j,k), gxz_p(i,j,k), & - gyy_p(i,j,k), gyz_p(i,j,k), gzz_p(i,j,k), det, & - dens_p(i,j,k), scon_p(i,j,k,1), scon_p(i,j,k,2), scon_p(i,j,k,3), & - tau_p(i,j,k), Bcons_p(i,j,k,1),Bcons_p(i,j,k,2),Bcons_p(i,j,k,3),& - rho_p(i,j,k), vel_p(i,j,k,1), vel_p(i,j,k,2), & - vel_p(i,j,k,3), eps_p(i,j,k), press_p(i,j,k), & - Bvec_p(i,j,k,1),Bvec_p(i,j,k,2),Bvec_p(i,j,k,3),w_lorentz_p(i,j,k)) - endif + call EOS_Omni_press(eos_handle,keytemp,GRHydro_eos_rf_prec,n,& + rho_p(i,j,k),eps_p(i,j,k),xtemp,xye,press_p(i,j,k),keyerr,anyerr) + call EOS_Omni_EpsFromPress(eos_handle,keytemp,GRHydro_eos_rf_prec,n,& + rho_p(i,j,k),xeps,xtemp,xye,press_p(i,j,k),eps_p(i,j,k),keyerr,anyerr) + call prim2conpolytypeM(eos_handle, & + gxx_p(i,j,k), gxy_p(i,j,k), gxz_p(i,j,k), & + gyy_p(i,j,k), gyz_p(i,j,k), gzz_p(i,j,k), det, & + dens_p(i,j,k), scon_p(i,j,k,1), scon_p(i,j,k,2), scon_p(i,j,k,3), & + tau_p(i,j,k), Bcons_p(i,j,k,1),Bcons_p(i,j,k,2),Bcons_p(i,j,k,3),& + rho_p(i,j,k), vel_p(i,j,k,1), vel_p(i,j,k,2), & + vel_p(i,j,k,3), eps_p(i,j,k), press_p(i,j,k), & + Bvec_p(i,j,k,1),Bvec_p(i,j,k,2),Bvec_p(i,j,k,3),w_lorentz_p(i,j,k)) + + end if + end if end if + if (timelevels .gt. 2) then if (rho_p_p(i,j,k) .le. GRHydro_rho_min) then rho_p_p(i,j,k) = GRHydro_rho_min vel_p_p(i,j,k,1) = 0.0d0 vel_p_p(i,j,k,2) = 0.0d0 vel_p_p(i,j,k,3) = 0.0d0 - call EOS_Omni_press(eos_handle,keytemp,GRHydro_eos_rf_prec,n,& - rho_p_p(i,j,k),eps_p_p(i,j,k),xtemp,xye,press_p_p(i,j,k),keyerr,anyerr) - call EOS_Omni_EpsFromPress(eos_handle,keytemp,GRHydro_eos_rf_prec,n,& - rho_p_p(i,j,k),xeps,xtemp,xye,press_p_p(i,j,k),eps_p_p(i,j,k),keyerr,anyerr) - det = SPATIAL_DETERMINANT(gxx_p_p(i,j,k), gxy_p_p(i,j,k), gxz_p_p(i,j,k), \ - gyy_p_p(i,j,k), gyz_p_p(i,j,k), gzz_p_p(i,j,k)) - call prim2conpolytypeM(eos_handle, & - gxx_p_p(i,j,k), gxy_p_p(i,j,k), gxz_p_p(i,j,k), & - gyy_p_p(i,j,k), gyz_p_p(i,j,k), gzz_p_p(i,j,k), det, & - dens_p_p(i,j,k), scon_p_p(i,j,k,1), scon_p_p(i,j,k,2), scon_p_p(i,j,k,3), & - tau_p_p(i,j,k), Bcons_p_p(i,j,k,1),Bcons_p_p(i,j,k,2),Bcons_p_p(i,j,k,3),& - rho_p_p(i,j,k), vel_p_p(i,j,k,1), vel_p_p(i,j,k,2), & - vel_p_p(i,j,k,3), eps_p_p(i,j,k), press_p_p(i,j,k), & - Bvec_p_p(i,j,k,1),Bvec_p_p(i,j,k,2),Bvec_p_p(i,j,k,3),w_lorentz_p_p(i,j,k)) - endif - endif + gyy_p_p(i,j,k), gyz_p_p(i,j,k), gzz_p_p(i,j,k)) + + if(evolve_temper.ne.0) then +! ! set the temperature to be relatively low + temperature_p_p(i,j,k) = grhydro_hot_atmo_temp + y_e_p_p(i,j,k) = grhydro_hot_atmo_Y_e + keytemp = 1 + call EOS_Omni_press(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,& + rho_p_p(i,j,k),eps_p_p(i,j,k),temperature_p_p(i,j,k),y_e_p_p(i,j,k),& + press_p_p(i,j,k),keyerr,anyerr) + call prim2conM_hot(GRHydro_eos_handle, GRHydro_reflevel,& + i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),gxx_p_p(i,j,k),& + gxy_p_p(i,j,k),gxz_p_p(i,j,k),gyy_p_p(i,j,k),gyz_p_p(i,j,k),gzz_p_p(i,j,k),& + det, dens_p_p(i,j,k),scon_p_p(i,j,k,1),scon_p_p(i,j,k,2),scon_p_p(i,j,k,3),& + tau_p_p(i,j,k),Bcons_p_p(i,j,k,1),Bcons_p_p(i,j,k,2),Bcons_p_p(i,j,k,3),& + rho_p_p(i,j,k),vel_p_p(i,j,k,1),vel_p_p(i,j,k,2),vel_p_p(i,j,k,3),& + eps_p_p(i,j,k),press_p_p(i,j,k),Bvec_p_p(i,j,k,1), & + Bvec_p_p(i,j,k,2), Bvec_p_p(i,j,k,3), w_lorentz_p_p(i,j,k),& + temperature_p_p(i,j,k),y_e_p_p(i,j,k)) + y_e_con_p_p(i,j,k) = dens_p_p(i,j,k) * y_e_p_p(i,j,k) + + else + + call EOS_Omni_press(eos_handle,keytemp,GRHydro_eos_rf_prec,n,& + rho_p_p(i,j,k),eps_p_p(i,j,k),xtemp,xye,press_p_p(i,j,k),keyerr,anyerr) + call EOS_Omni_EpsFromPress(eos_handle,keytemp,GRHydro_eos_rf_prec,n,& + rho_p_p(i,j,k),xeps,xtemp,xye,press_p_p(i,j,k),eps_p_p(i,j,k),keyerr,anyerr) + call prim2conpolytypeM(eos_handle, & + gxx_p_p(i,j,k), gxy_p_p(i,j,k), gxz_p_p(i,j,k), & + gyy_p_p(i,j,k), gyz_p_p(i,j,k), gzz_p_p(i,j,k), det, & + dens_p_p(i,j,k), scon_p_p(i,j,k,1), scon_p_p(i,j,k,2), scon_p_p(i,j,k,3), & + tau_p_p(i,j,k), Bcons_p_p(i,j,k,1),Bcons_p_p(i,j,k,2),Bcons_p_p(i,j,k,3),& + rho_p_p(i,j,k), vel_p_p(i,j,k,1), vel_p_p(i,j,k,2), & + vel_p_p(i,j,k,3), eps_p_p(i,j,k), press_p_p(i,j,k), & + Bvec_p_p(i,j,k,1),Bvec_p_p(i,j,k,2),Bvec_p_p(i,j,k,3),w_lorentz_p_p(i,j,k)) + + end if + end if + end if end do end do diff --git a/src/make.code.defn b/src/make.code.defn index 25b0aa1..b1042c7 100644 --- a/src/make.code.defn +++ b/src/make.code.defn @@ -9,7 +9,7 @@ SRCS = Utils.F90 \ GRHydro_CalcBcom.F90 \ GRHydro_CalcUpdate.F90 \ GRHydro_Con2Prim.F90 \ - GRHydro_DivergenceClean.F90 \ + GRHydro_DivergenceClean.F90 \ GRHydro_Eigenproblem.F90 \ GRHydro_Eigenproblem_Marquina.F90 \ GRHydro_ENOReconstruct.F90 \ @@ -45,12 +45,12 @@ SRCS = Utils.F90 \ GRHydro_EoSChangeGamma.F90 \ GRHydro_Tmunu.F90 \ GRHydro_Con2PrimM.F90 \ - GRHydro_Con2PrimM_pt.c \ + GRHydro_Con2PrimM_pt_EOSOmni.c \ GRHydro_Con2PrimM_polytype_pt.c \ GRHydro_EigenproblemM.F90 \ GRHydro_FluxM.F90 \ GRHydro_HLLEM.F90 \ - GRHydro_PPMM.F90 \ + GRHydro_PPMM.F90 \ GRHydro_Prim2ConM.F90 \ GRHydro_RiemannSolveM.F90 \ GRHydro_SourceM.F90 \ |