diff options
| author | rhaas <rhaas@c83d129a-5a75-4d5a-9c4d-ed3a5855bf45> | 2013-05-11 06:29:34 +0000 |
|---|---|---|
| committer | rhaas <rhaas@c83d129a-5a75-4d5a-9c4d-ed3a5855bf45> | 2013-05-11 06:29:34 +0000 |
| commit | 007c18aa5e6aed27c7edbf71b1f03267c0ffef9d (patch) | |
| tree | 6d9c4331bcf20e8fc4f7fc2ab4fbb6802b2293a0 | |
| parent | d73ce95f867e1a041cc51b82a954b6578e24610c (diff) | |
GRHydro: Fixing Con2PrimM to call the proper point-wise routines
This re-introduces routines that work for hybrid/hot EOS and
corresponding changes in pointwise routines for hot EOS error checking
and temperature treatment by adding old EOSOmni pointwise routine.
From: Philipp Moesta <pmoesta@tapir.caltech.edu>
git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEvolve/GRHydro/trunk@511 c83d129a-5a75-4d5a-9c4d-ed3a5855bf45
| -rw-r--r-- | interface.ccl | 40 | ||||
| -rw-r--r-- | src/GRHydro_Con2PrimAM.F90 | 20 | ||||
| -rw-r--r-- | src/GRHydro_Con2PrimM.F90 | 26 | ||||
| -rw-r--r-- | src/GRHydro_Con2PrimM_pt.c | 75 | ||||
| -rw-r--r-- | src/GRHydro_Con2PrimM_pt_EOSOmni.c | 448 | ||||
| -rw-r--r-- | src/GRHydro_Con2PrimM_pt_EOSOmniold.c | 1056 | ||||
| -rw-r--r-- | src/GRHydro_Con2PrimM_ptee.c | 66 | ||||
| -rw-r--r-- | src/GRHydro_InterfacesM.h | 94 | ||||
| -rw-r--r-- | src/make.code.defn | 2 |
9 files changed, 1649 insertions, 178 deletions
diff --git a/interface.ccl b/interface.ccl index db6615d..e074c2e 100644 --- a/interface.ccl +++ b/interface.ccl @@ -67,7 +67,10 @@ void FUNCTION Con2PrimGen(CCTK_INT INOUT handle, CCTK_REAL INOUT dens, \ CCTK_INT INOUT GRHydro_reflevel, \ CCTK_REAL OUT retval) -void FUNCTION Con2PrimGenM(CCTK_INT INOUT handle, CCTK_INT INOUT keytemp, CCTK_REAL INOUT prec, CCTK_REAL INOUT gamma_eos, CCTK_REAL INOUT dens, \ +void FUNCTION Con2PrimGenM2(CCTK_INT IN handle, CCTK_INT IN rl, CCTK_INT IN ii, CCTK_INT IN jj, \ + CCTK_INT IN kk, CCTK_REAL IN x, CCTK_REAL IN y, CCTK_REAL IN z, \ + CCTK_INT IN keytemp, CCTK_REAL IN eos_prec, CCTK_REAL IN prec, \ + CCTK_REAL INOUT gamma_eos, CCTK_REAL INOUT dens, \ CCTK_REAL INOUT sx, CCTK_REAL INOUT sy, CCTK_REAL INOUT sz, \ CCTK_REAL INOUT tau, CCTK_REAL INOUT Bconsx, CCTK_REAL INOUT Bconsy, CCTK_REAL INOUT Bconsz, \ CCTK_REAL INOUT y_e, CCTK_REAL INOUT temp, CCTK_REAL INOUT rho, \ @@ -105,22 +108,22 @@ void FUNCTION Con2PrimGenMee(CCTK_INT INOUT handle, CCTK_INT INOUT keytemp, \ CCTK_INT OUT epsnegative, \ CCTK_REAL OUT retval) -#void FUNCTION Con2PrimGenM(CCTK_INT INOUT handle, CCTK_INT INOUT keytemp, CCTK_REAL INOUT prec,CCTK_REAL INOUT gamma_eos, CCTK_REAL INOUT dens, \ -# CCTK_REAL INOUT sx, CCTK_REAL INOUT sy, CCTK_REAL INOUT sz, \ -# CCTK_REAL INOUT tau, CCTK_REAL INOUT Bconsx, CCTK_REAL INOUT Bconsy, CCTK_REAL INOUT Bconsz, \ -# CCTK_REAL INOUT y_e, CCTK_REAL INOUT temp, CCTK_REAL INOUT rho, \ -# CCTK_REAL INOUT velx, CCTK_REAL INOUT vely, CCTK_REAL INOUT velz, \ -# CCTK_REAL INOUT epsilon, CCTK_REAL INOUT pressure, \ -# CCTK_REAL INOUT Bvecx, CCTK_REAL INOUT Bvecy, CCTK_REAL INOUT Bvecz, \ -# CCTK_REAL INOUT Bvecsq, \ -# CCTK_REAL INOUT w_lorentz, \ -# CCTK_REAL INOUT gxx, CCTK_REAL INOUT gxy, CCTK_REAL INOUT gxz, \ -# CCTK_REAL INOUT gyy, CCTK_REAL INOUT gyz, CCTK_REAL INOUT gzz, \ -# CCTK_REAL INOUT uxx, CCTK_REAL INOUT uxy, CCTK_REAL INOUT uxz, \ -# CCTK_REAL INOUT uyy, CCTK_REAL INOUT uyz, CCTK_REAL INOUT uzz, \ -# CCTK_REAL INOUT det, \ -# CCTK_INT OUT epsnegative, \ -# CCTK_REAL OUT retval) +void FUNCTION Con2PrimGenM(CCTK_INT INOUT handle, CCTK_INT INOUT keytemp, CCTK_REAL INOUT prec,CCTK_REAL INOUT gamma_eos, CCTK_REAL INOUT dens, \ + CCTK_REAL INOUT sx, CCTK_REAL INOUT sy, CCTK_REAL INOUT sz, \ + CCTK_REAL INOUT tau, CCTK_REAL INOUT Bconsx, CCTK_REAL INOUT Bconsy, CCTK_REAL INOUT Bconsz, \ + CCTK_REAL INOUT y_e, CCTK_REAL INOUT temp, CCTK_REAL INOUT rho, \ + CCTK_REAL INOUT velx, CCTK_REAL INOUT vely, CCTK_REAL INOUT velz, \ + CCTK_REAL INOUT epsilon, CCTK_REAL INOUT pressure, \ + CCTK_REAL INOUT Bvecx, CCTK_REAL INOUT Bvecy, CCTK_REAL INOUT Bvecz, \ + CCTK_REAL INOUT Bvecsq, \ + CCTK_REAL INOUT w_lorentz, \ + CCTK_REAL INOUT gxx, CCTK_REAL INOUT gxy, CCTK_REAL INOUT gxz, \ + CCTK_REAL INOUT gyy, CCTK_REAL INOUT gyz, CCTK_REAL INOUT gzz, \ + CCTK_REAL INOUT uxx, CCTK_REAL INOUT uxy, CCTK_REAL INOUT uxz, \ + CCTK_REAL INOUT uyy, CCTK_REAL INOUT uyz, CCTK_REAL INOUT uzz, \ + CCTK_REAL INOUT det, \ + CCTK_INT OUT epsnegative, \ + CCTK_REAL OUT retval) void FUNCTION Con2PrimPoly(CCTK_INT IN handle, CCTK_REAL OUT dens, \ CCTK_REAL OUT sx, CCTK_REAL OUT sy, \ @@ -227,7 +230,8 @@ PROVIDES FUNCTION SpatialDet WITH SpatialDeterminant LANGUAGE Fortran PROVIDES FUNCTION UpperMet WITH UpperMetric LANGUAGE Fortran #PROVIDES FUNCTION Con2Prim WITH Con2Prim_pt LANGUAGE Fortran PROVIDES FUNCTION Con2PrimPoly WITH Con2Prim_ptPolytype LANGUAGE Fortran -PROVIDES FUNCTION Con2PrimGenM WITH GRHydro_Con2PrimM_pt LANGUAGE Fortran +PROVIDES FUNCTION Con2PrimGenM WITH GRHydro_Con2PrimM_ptold LANGUAGE Fortran +PROVIDES FUNCTION Con2PrimGenM2 WITH GRHydro_Con2PrimM_pt LANGUAGE Fortran PROVIDES FUNCTION Con2PrimGenMee WITH GRHydro_Con2PrimM_ptee LANGUAGE Fortran PROVIDES FUNCTION Con2PrimGen WITH Con2Prim_pt LANGUAGE Fortran PROVIDES FUNCTION Con2PrimPolyM WITH GRHydro_Con2PrimM_Polytype_pt LANGUAGE Fortran diff --git a/src/GRHydro_Con2PrimAM.F90 b/src/GRHydro_Con2PrimAM.F90 index 5bdf02a..6b8d894 100644 --- a/src/GRHydro_Con2PrimAM.F90 +++ b/src/GRHydro_Con2PrimAM.F90 @@ -187,7 +187,6 @@ subroutine Conservative2PrimitiveAM(CCTK_ARGUMENTS) det = SPATIAL_DETERMINANT(g11(i,j,k),g12(i,j,k),g13(i,j,k),g22(i,j,k),g23(i,j,k),g33(i,j,k)) sdet = sqrt(det) - 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)) @@ -345,9 +344,12 @@ subroutine Conservative2PrimitiveAM(CCTK_ARGUMENTS) Bconsz_tmp = sdet*Bvecz_tmp keytemp = 0 - - call GRHydro_Con2PrimM_pt(GRHydro_eos_handle, keytemp, & - GRHydro_eos_rf_prec, local_gam(1), dens(i,j,k), & + !Watch out for the values returned to b2. Here b2 is the Bprim^2 + !while inside the point-wise con2prim routines it is the square + !of the comoving B-field, b^{\mu} b_{\mu}. It is overwritten + !in this routine, but we may need to find a better notation + !avoid future confusions. + call GRHydro_Con2PrimM_pt2(GRHydro_eos_handle, keytemp, GRHydro_eos_rf_prec, GRHydro_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), & Bconsx_tmp, Bconsy_tmp, Bconsz_tmp,xye(1), & xtemp(1),rho_tmp,velx_tmp,vely_tmp,velz_tmp,& @@ -488,8 +490,7 @@ subroutine Conservative2PrimitiveAM(CCTK_ARGUMENTS) keytemp = 0 - call GRHydro_Con2PrimM_pt(GRHydro_eos_handle, keytemp, & - GRHydro_eos_rf_prec, local_gam(1), dens(i,j,k), & + call GRHydro_Con2PrimM_pt(GRHydro_eos_handle, GRHydro_reflevel, i ,j ,k, x(i,j,k), y(i,j,k), z(i,j,k), keytemp, GRHydro_eos_rf_prec, GRHydro_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), & Bconsx_tmp, Bconsy_tmp, Bconsz_tmp, Y_e(i,j,k), & temperature(i,j,k),rho_tmp,velx_tmp,vely_tmp,velz_tmp,& @@ -518,8 +519,7 @@ subroutine Conservative2PrimitiveAM(CCTK_ARGUMENTS) Bconsx_tmp = sdet*Bvecx_tmp Bconsy_tmp = sdet*Bvecy_tmp Bconsz_tmp = sdet*Bvecz_tmp - call GRHydro_Con2PrimM_pt(GRHydro_eos_handle, keytemp, & - local_perc_ptol, local_gam(1), dens(i,j,k), & + call GRHydro_Con2PrimM_pt(GRHydro_eos_handle, GRHydro_reflevel, i ,j ,k, x(i,j,k), y(i,j,k), z(i,j,k), keytemp, GRHydro_eos_rf_prec, GRHydro_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), & Bconsx_tmp, Bconsy_tmp, Bconsz_tmp, Y_e(i,j,k), & temperature(i,j,k),rho_tmp,velx_tmp,vely_tmp,velz_tmp,& @@ -873,7 +873,7 @@ subroutine Conservative2PrimitiveBoundsAM(CCTK_ARGUMENTS) Bconsx_tmp = sqrt(avg_detl)*Bvecxminus(i,j,k) Bconsy_tmp = sqrt(avg_detl)*Bvecyminus(i,j,k) Bconsz_tmp = sqrt(avg_detl)*Bveczminus(i,j,k) - call GRHydro_Con2PrimM_pt(GRHydro_eos_handle, keytemp, GRHydro_eos_rf_prec, local_gam(1),densminus(i,j,k), & + call GRHydro_Con2PrimM_pt2(GRHydro_eos_handle, keytemp, GRHydro_eos_rf_prec, GRHydro_perc_ptol, local_gam(1), densminus(i,j,k), & sxminus(i,j,k),syminus(i,j,k),szminus(i,j,k), tauminus(i,j,k), & Bconsx_tmp, Bconsy_tmp, Bconsz_tmp, xye(1), xtemp(1), rhominus(i,j,k),& velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k),epsminus(i,j,k),pressminus(i,j,k),& @@ -934,7 +934,7 @@ subroutine Conservative2PrimitiveBoundsAM(CCTK_ARGUMENTS) Bconsx_tmp = sqrt(avg_detr)*Bvecxplus(i,j,k) Bconsy_tmp = sqrt(avg_detr)*Bvecyplus(i,j,k) Bconsz_tmp = sqrt(avg_detr)*Bveczplus(i,j,k) - call GRHydro_Con2PrimM_pt(GRHydro_eos_handle, keytemp, GRHydro_eos_rf_prec, local_gam(1),densplus(i,j,k), & + call GRHydro_Con2PrimM_pt2(GRHydro_eos_handle, keytemp, GRHydro_eos_rf_prec, GRHydro_perc_ptol, local_gam(1), densplus(i,j,k), & sxplus(i,j,k),syplus(i,j,k),szplus(i,j,k), tauplus(i,j,k), & Bconsx_tmp, Bconsy_tmp, Bconsz_tmp, xye(1), xtemp(1), rhoplus(i,j,k),& velxplus(i,j,k),velyplus(i,j,k),velzplus(i,j,k),epsplus(i,j,k),pressplus(i,j,k),& diff --git a/src/GRHydro_Con2PrimM.F90 b/src/GRHydro_Con2PrimM.F90 index 0960b0d..720e826 100644 --- a/src/GRHydro_Con2PrimM.F90 +++ b/src/GRHydro_Con2PrimM.F90 @@ -339,8 +339,19 @@ subroutine Conservative2PrimitiveM(CCTK_ARGUMENTS) !of the comoving B-field, b^{\mu} b_{\mu}. It is overwritten !in this routine, but we may need to find a better notation !avoid future confusions. - call GRHydro_Con2PrimM_pt(GRHydro_eos_handle, keytemp, & - GRHydro_eos_rf_prec, local_gam(1), dens(i,j,k), & + if(GRHydro_eos_handle .eq. 1 .or. GRHydro_eos_handle .eq. 2) then + + call GRHydro_Con2PrimM_ptold(GRHydro_eos_handle, 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_tmp, & + velx_tmp,vely_tmp,velz_tmp,eps_tmp,press_tmp, & + Bvecx_tmp,Bvecy_tmp,Bvecz_tmp,b2, w_lorentz_tmp,& + g11(i,j,k),g12(i,j,k),g13(i,j,k),g22(i,j,k),g23(i,j,k),g33(i,j,k), & + uxx,uxy,uxz,uyy,uyz,uzz,det, & + epsnegative,GRHydro_C2P_failed(i,j,k)) + + else + call GRHydro_Con2PrimM_pt2(GRHydro_eos_handle, keytemp, GRHydro_eos_rf_prec, GRHydro_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),xye(1), & xtemp(1),rho_tmp,velx_tmp,vely_tmp,velz_tmp,& @@ -349,6 +360,7 @@ subroutine Conservative2PrimitiveM(CCTK_ARGUMENTS) g22(i,j,k),g23(i,j,k),g33(i,j,k), & uxx,uxy,uxz,uyy,uyz,uzz,det, & epsnegative,GRHydro_C2P_failed(i,j,k)) + endif if(evolve_entropy.ne.0) then if(GRHydro_C2P_failed(i,j,k).ne.0) then @@ -492,8 +504,7 @@ subroutine Conservative2PrimitiveM(CCTK_ARGUMENTS) keytemp = 0 - call GRHydro_Con2PrimM_pt(GRHydro_eos_handle, keytemp, & - GRHydro_eos_rf_prec, local_gam(1), dens(i,j,k), & + call GRHydro_Con2PrimM_pt(GRHydro_eos_handle, GRHydro_reflevel, i ,j ,k, x(i,j,k), y(i,j,k), z(i,j,k), keytemp, GRHydro_eos_rf_prec, GRHydro_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_tmp,velx_tmp,vely_tmp,velz_tmp,& @@ -519,8 +530,7 @@ subroutine Conservative2PrimitiveM(CCTK_ARGUMENTS) Bvecx_tmp = Bprim(i,j,k,1) Bvecy_tmp = Bprim(i,j,k,2) Bvecz_tmp = Bprim(i,j,k,3) - call GRHydro_Con2PrimM_pt(GRHydro_eos_handle, keytemp, & - local_perc_ptol, local_gam(1), dens(i,j,k), & + call GRHydro_Con2PrimM_pt(GRHydro_eos_handle, GRHydro_reflevel, i ,j ,k, x(i,j,k), y(i,j,k), z(i,j,k), keytemp, GRHydro_eos_rf_prec, GRHydro_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_tmp,velx_tmp,vely_tmp,velz_tmp,& @@ -864,7 +874,7 @@ 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, keytemp, GRHydro_eos_rf_prec, local_gam(1),densminus(i,j,k), & + call GRHydro_Con2PrimM_pt2(GRHydro_eos_handle, keytemp, GRHydro_eos_rf_prec, GRHydro_perc_ptol, local_gam(1), densminus(i,j,k), & sxminus(i,j,k),syminus(i,j,k),szminus(i,j,k), tauminus(i,j,k), & 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),& @@ -922,7 +932,7 @@ subroutine Conservative2PrimitiveBoundsM(CCTK_ARGUMENTS) epsnegative = 0 - call GRHydro_Con2PrimM_pt(GRHydro_eos_handle, keytemp, GRHydro_eos_rf_prec, local_gam(1),densplus(i,j,k), & + call GRHydro_Con2PrimM_pt2(GRHydro_eos_handle, keytemp, GRHydro_eos_rf_prec, GRHydro_perc_ptol, local_gam(1), densplus(i,j,k), & sxplus(i,j,k),syplus(i,j,k),szplus(i,j,k), tauplus(i,j,k), & 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),& diff --git a/src/GRHydro_Con2PrimM_pt.c b/src/GRHydro_Con2PrimM_pt.c index b980887..bf947d3 100644 --- a/src/GRHydro_Con2PrimM_pt.c +++ b/src/GRHydro_Con2PrimM_pt.c @@ -61,7 +61,7 @@ #define NEWT_TOL2 (1.0e-15) /* TOL of new 1D^*_{v^2} gnr2 method */ #define MIN_NEWT_TOL2 (1.0e-10) /* TOL of new 1D^*_{v^2} gnr2 method */ -#define W_TOO_BIG (1.e20) /* \gamma^2 (\rho_0 + u + p) is assumed +#define W_TOO_BIG (1.e20) /* \gamma^2 (\rho_0 + u + p) is assumed to always be smaller than this. This is used to detect solver failures */ @@ -104,13 +104,12 @@ static CCTK_REAL pressure_rho0_w(CCTK_REAL rho0, CCTK_REAL w,CCTK_REAL gammaeos) return((gammaeos-1.)*(w - rho0)/gammaeos) ; } -void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt) ( - CCTK_INT *handle, CCTK_INT *keytemp, CCTK_REAL *prec, - CCTK_REAL *gamma_eos, + +void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_ptold) ( + CCTK_INT *handle, CCTK_REAL *gamma_eos, CCTK_REAL *dens_in, CCTK_REAL *sx_in, CCTK_REAL *sy_in, CCTK_REAL *sz_in, CCTK_REAL *tau_in, CCTK_REAL *Bconsx_in, CCTK_REAL *Bconsy_in, CCTK_REAL *Bconsz_in, - CCTK_REAL *y_e_in, CCTK_REAL *temp_in, CCTK_REAL *rho, CCTK_REAL *velx, CCTK_REAL *vely, CCTK_REAL *velz, CCTK_REAL *epsilon, CCTK_REAL *pressure, @@ -186,15 +185,12 @@ 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, +void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_ptold) ( + CCTK_INT *handle, CCTK_REAL *gamma_eos, CCTK_REAL *dens_in, CCTK_REAL *sx_in, CCTK_REAL *sy_in, CCTK_REAL *sz_in, CCTK_REAL *tau_in, CCTK_REAL *Bconsx_in, CCTK_REAL *Bconsy_in, CCTK_REAL *Bconsz_in, - CCTK_REAL *y_e_in, CCTK_REAL* temp_in, CCTK_REAL *rho, CCTK_REAL *velx, CCTK_REAL *vely, CCTK_REAL *velz, CCTK_REAL *epsilon, CCTK_REAL *pressure, @@ -229,43 +225,42 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt) ( gammaeos = *gamma_eos; /* Assume ok initially: */ - // BCM: But let the driver function take care of its initialization - //*retval = 0.; + *retval = 0.; *epsnegative = 0; #if(DEBUG_CON2PRIMM) fprintf(stdout," *dens = %26.16e \n", *dens_in ); - fprintf(stdout," *sx = %26.16e \n", *sx_in ); - fprintf(stdout," *sy = %26.16e \n", *sy_in ); - fprintf(stdout," *sz = %26.16e \n", *sz_in ); - fprintf(stdout," *tau = %26.16e \n", *tau_in ); - fprintf(stdout," *Bconsx = %26.16e \n", *Bconsx_in ); - fprintf(stdout," *Bconsy = %26.16e \n", *Bconsy_in ); + fprintf(stdout," *sx = %26.16e \n", *sx_in ); + fprintf(stdout," *sy = %26.16e \n", *sy_in ); + fprintf(stdout," *sz = %26.16e \n", *sz_in ); + fprintf(stdout," *tau = %26.16e \n", *tau_in ); + fprintf(stdout," *Bconsx = %26.16e \n", *Bconsx_in ); + fprintf(stdout," *Bconsy = %26.16e \n", *Bconsy_in ); fprintf(stdout," *Bconsz = %26.16e \n", *Bconsz_in ); - fprintf(stdout," *rho = %26.16e \n", *rho ); - fprintf(stdout," *velx = %26.16e \n", *velx ); - fprintf(stdout," *vely = %26.16e \n", *vely ); - fprintf(stdout," *velz = %26.16e \n", *velz ); + fprintf(stdout," *rho = %26.16e \n", *rho ); + fprintf(stdout," *velx = %26.16e \n", *velx ); + fprintf(stdout," *vely = %26.16e \n", *vely ); + fprintf(stdout," *velz = %26.16e \n", *velz ); fprintf(stdout," *epsilon = %26.16e \n", *epsilon ); fprintf(stdout," *pressure = %26.16e \n", *pressure ); - fprintf(stdout," *Bx = %26.16e \n", *Bx ); - fprintf(stdout," *By = %26.16e \n", *By ); - fprintf(stdout," *Bz = %26.16e \n", *Bz ); - fprintf(stdout," *bsq = %26.16e \n", *bsq ); + fprintf(stdout," *Bx = %26.16e \n", *Bx ); + fprintf(stdout," *By = %26.16e \n", *By ); + fprintf(stdout," *Bz = %26.16e \n", *Bz ); + fprintf(stdout," *bsq = %26.16e \n", *bsq ); fprintf(stdout," *w_lorentz = %26.16e \n", *w_lorentz ); - fprintf(stdout," *gxx = %26.16e \n", *gxx ); - fprintf(stdout," *gxy = %26.16e \n", *gxy ); - fprintf(stdout," *gxz = %26.16e \n", *gxz ); - fprintf(stdout," *gyy = %26.16e \n", *gyy ); - fprintf(stdout," *gyz = %26.16e \n", *gyz ); - fprintf(stdout," *gzz = %26.16e \n", *gzz ); - fprintf(stdout," *uxx = %26.16e \n", *uxx ); - fprintf(stdout," *uxy = %26.16e \n", *uxy ); - fprintf(stdout," *uxz = %26.16e \n", *uxz ); - fprintf(stdout," *uyy = %26.16e \n", *uyy ); - fprintf(stdout," *uyz = %26.16e \n", *uyz ); - fprintf(stdout," *uzz = %26.16e \n", *uzz ); - fprintf(stdout," *det = %26.16e \n", *det ); + fprintf(stdout," *gxx = %26.16e \n", *gxx ); + fprintf(stdout," *gxy = %26.16e \n", *gxy ); + fprintf(stdout," *gxz = %26.16e \n", *gxz ); + fprintf(stdout," *gyy = %26.16e \n", *gyy ); + fprintf(stdout," *gyz = %26.16e \n", *gyz ); + fprintf(stdout," *gzz = %26.16e \n", *gzz ); + fprintf(stdout," *uxx = %26.16e \n", *uxx ); + fprintf(stdout," *uxy = %26.16e \n", *uxy ); + fprintf(stdout," *uxz = %26.16e \n", *uxz ); + fprintf(stdout," *uyy = %26.16e \n", *uyy ); + fprintf(stdout," *uyz = %26.16e \n", *uyz ); + fprintf(stdout," *uzz = %26.16e \n", *uzz ); + fprintf(stdout," *det = %26.16e \n", *det ); fprintf(stdout," *epsnegative = %10d \n", *epsnegative ); fprintf(stdout," *retval = %26.16e \n", *retval ); fflush(stdout); @@ -472,7 +467,6 @@ static CCTK_REAL x1_of_x0(CCTK_REAL x0, struct LocGlob *lgp ) } -#if(DEBUG_CON2PRIMM) /******************************************************************** validate_x(): @@ -497,7 +491,6 @@ static void validate_x(CCTK_REAL x[2], CCTK_REAL x0[2] ) return; } -#endif /************************************************************ diff --git a/src/GRHydro_Con2PrimM_pt_EOSOmni.c b/src/GRHydro_Con2PrimM_pt_EOSOmni.c index 1d7e15c..96e2bd8 100644 --- a/src/GRHydro_Con2PrimM_pt_EOSOmni.c +++ b/src/GRHydro_Con2PrimM_pt_EOSOmni.c @@ -61,7 +61,7 @@ #define NEWT_TOL2 (1.0e-15) /* TOL of new 1D^*_{v^2} gnr2 method */ #define MIN_NEWT_TOL2 (1.0e-10) /* TOL of new 1D^*_{v^2} gnr2 method */ -#define W_TOO_BIG (1.e20) /* \gamma^2 (\rho_0 + u + p) is assumed +#define W_TOO_BIG (1.e20) /* \gamma^2 (\rho_0 + u + p) is assumed to always be smaller than this. This is used to detect solver failures */ @@ -73,13 +73,15 @@ /* Local Globals */ struct LocGlob { - CCTK_REAL Bsq,QdotBsq,Qtsq,Qdotn,D,half_Bsq ; + CCTK_REAL Bsq,bsq,QdotBsq,VdotBsq,Qtsq,Qdotn,D,half_Bsq,tau ; } ; struct eosomnivars { CCTK_INT eoshandle,eoskeytemp; - CCTK_REAL eosprec,eos_y_e[1],eos_temp[1]; + CCTK_REAL prec,eosprec,eos_y_e[1],eos_temp[1]; CCTK_INT eoskeyerr[1],eosanyerr[1]; + CCTK_INT eosrl,eosit,eosi,eosj,eosk; + CCTK_REAL eosx,eosy,eosz; } ; // Declarations: @@ -122,9 +124,9 @@ static CCTK_REAL pressure_rho0_w(CCTK_REAL rho0, CCTK_REAL w,CCTK_REAL gammaeos) { return((gammaeos-1.)*(w - rho0)/gammaeos) ; } - +//CCTK_INT *handle, CCTK_INT GRHydro_reflevel, CCTK_INT *keytemp, CCTK_REAL *prec, void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt) ( - CCTK_INT *handle, CCTK_INT *keytemp, CCTK_REAL *prec, + CCTK_INT *handle, CCTK_INT *rl, CCTK_INT *ii, CCTK_INT *jj, CCTK_INT *kk, CCTK_REAL *x, CCTK_REAL *y, CCTK_REAL *z, CCTK_INT *keytemp, CCTK_REAL *eos_prec, CCTK_REAL *prec, CCTK_REAL *gamma_eos, CCTK_REAL *dens_in, CCTK_REAL *sx_in, CCTK_REAL *sy_in, CCTK_REAL *sz_in, @@ -172,7 +174,7 @@ INPUT: (using GRHydro variable defintions) tau = \alpha^2 \sqrt(\gamma) T^{00} - D g[x,y,z][x,y,x] = spatial metric corresponding to \gamma u[x,y,z][x,y,z] = inverse of the spatial metric, g[x,y,z][x,y,x] - det = gamma + det = sqrt(\gamma) B[x,y,z] = Bvec[0,1,2] bsq = b^\mu b_\mu @@ -205,10 +207,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_INT *handle, CCTK_INT *rl, CCTK_INT *ii, CCTK_INT *jj, CCTK_INT *kk, CCTK_REAL *x, CCTK_REAL *y, CCTK_REAL *z, CCTK_INT *keytemp, CCTK_REAL *eos_prec, CCTK_REAL *prec, CCTK_REAL *gamma_eos, CCTK_REAL *dens_in, CCTK_REAL *sx_in, CCTK_REAL *sy_in, CCTK_REAL *sz_in, @@ -234,14 +234,14 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt) ( CCTK_REAL sx, sy, sz; CCTK_REAL usx, usy, usz; CCTK_REAL tau, dens, gammaeos; - CCTK_REAL QdotB; + CCTK_REAL QdotB,VdotB; CCTK_REAL rho0,u,p,w,gammasq,gamma,gtmp,W_last,W,vsq; CCTK_REAL g_o_WBsq, QdB_o_W; CCTK_REAL detg = (*det); CCTK_REAL sqrt_detg = sqrt(detg); CCTK_REAL inv_sqrt_detg = 1./sqrt_detg; CCTK_INT i,j, i_increase ; - + CCTK_INT nf, nfudgemax, failwarnmode, failinfomode; DECLARE_CCTK_PARAMETERS; struct LocGlob lg; @@ -250,53 +250,70 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt) ( eosvars.eoshandle = *handle; //printf("handle = %i\n",*handle); eosvars.eoskeytemp = *keytemp; - eosvars.eosprec = *prec; + eosvars.eosprec = *eos_prec; + eosvars.prec = *prec; eosvars.eos_y_e[0] = *y_e_in; eosvars.eos_temp[0] = *temp_in; eosvars.eoskeyerr[0] = 0; eosvars.eosanyerr[0] = 0; - + eosvars.eosrl = *rl; + eosvars.eosit = 1; + eosvars.eosi = *ii; + eosvars.eosj = *jj; + eosvars.eosk = *kk; + eosvars.eosx = *x; + eosvars.eosy = *y; + eosvars.eosz = *z; + gammaeos = *gamma_eos; + nf =0; + nfudgemax = 30; + failinfomode = 1; + failwarnmode = 0; + /* Assume ok initially: */ *retval = 0.; *epsnegative = 0; #if(DEBUG_CON2PRIMM) + fprintf(stdout," *x1 = %26.16e \n", *x ); + fprintf(stdout," *y1 = %26.16e \n", *y ); + fprintf(stdout," *z1 = %26.16e \n", *z ); fprintf(stdout," *dens = %26.16e \n", *dens_in ); - fprintf(stdout," *sx = %26.16e \n", *sx_in ); - fprintf(stdout," *sy = %26.16e \n", *sy_in ); - fprintf(stdout," *sz = %26.16e \n", *sz_in ); - fprintf(stdout," *tau = %26.16e \n", *tau_in ); - fprintf(stdout," *Bconsx = %26.16e \n", *Bconsx_in ); - fprintf(stdout," *Bconsy = %26.16e \n", *Bconsy_in ); + fprintf(stdout," *sx = %26.16e \n", *sx_in ); + fprintf(stdout," *sy = %26.16e \n", *sy_in ); + fprintf(stdout," *sz = %26.16e \n", *sz_in ); + fprintf(stdout," *tau = %26.16e \n", *tau_in ); + fprintf(stdout," *Bconsx = %26.16e \n", *Bconsx_in ); + fprintf(stdout," *Bconsy = %26.16e \n", *Bconsy_in ); fprintf(stdout," *Bconsz = %26.16e \n", *Bconsz_in ); - fprintf(stdout," *rho = %26.16e \n", *rho ); - fprintf(stdout," *velx = %26.16e \n", *velx ); - fprintf(stdout," *vely = %26.16e \n", *vely ); - fprintf(stdout," *velz = %26.16e \n", *velz ); + fprintf(stdout," *rho = %26.16e \n", *rho ); + fprintf(stdout," *velx = %26.16e \n", *velx ); + 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 ); - fprintf(stdout," *Bz = %26.16e \n", *Bz ); - fprintf(stdout," *bsq = %26.16e \n", *bsq ); + fprintf(stdout," *Bx = %26.16e \n", *Bx ); + fprintf(stdout," *By = %26.16e \n", *By ); + fprintf(stdout," *Bz = %26.16e \n", *Bz ); + fprintf(stdout," *bsq = %26.16e \n", *bsq ); fprintf(stdout," *w_lorentz = %26.16e \n", *w_lorentz ); - fprintf(stdout," *gxx = %26.16e \n", *gxx ); - fprintf(stdout," *gxy = %26.16e \n", *gxy ); - fprintf(stdout," *gxz = %26.16e \n", *gxz ); - fprintf(stdout," *gyy = %26.16e \n", *gyy ); - fprintf(stdout," *gyz = %26.16e \n", *gyz ); - fprintf(stdout," *gzz = %26.16e \n", *gzz ); - fprintf(stdout," *uxx = %26.16e \n", *uxx ); - fprintf(stdout," *uxy = %26.16e \n", *uxy ); - fprintf(stdout," *uxz = %26.16e \n", *uxz ); - fprintf(stdout," *uyy = %26.16e \n", *uyy ); - fprintf(stdout," *uyz = %26.16e \n", *uyz ); - fprintf(stdout," *uzz = %26.16e \n", *uzz ); - fprintf(stdout," *det = %26.16e \n", *det ); + fprintf(stdout," *gxx = %26.16e \n", *gxx ); + fprintf(stdout," *gxy = %26.16e \n", *gxy ); + fprintf(stdout," *gxz = %26.16e \n", *gxz ); + fprintf(stdout," *gyy = %26.16e \n", *gyy ); + fprintf(stdout," *gyz = %26.16e \n", *gyz ); + fprintf(stdout," *gzz = %26.16e \n", *gzz ); + fprintf(stdout," *uxx = %26.16e \n", *uxx ); + fprintf(stdout," *uxy = %26.16e \n", *uxy ); + fprintf(stdout," *uxz = %26.16e \n", *uxz ); + fprintf(stdout," *uyy = %26.16e \n", *uyy ); + fprintf(stdout," *uyz = %26.16e \n", *uyz ); + fprintf(stdout," *uzz = %26.16e \n", *uzz ); + fprintf(stdout," *det = %26.16e \n", *det ); fprintf(stdout," *epsnegative = %10d \n", *epsnegative ); fprintf(stdout," *retval = %26.16e \n", *retval ); fflush(stdout); @@ -329,6 +346,13 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt) ( (*gyz) * (*By) * (*Bz) ); QdotB = (sx * (*Bx) + sy * (*By) + sz * (*Bz)) ; + + VdotB = ((*gxx) * (*velx) + (*gxy) * (*vely) + (*gxz) * (*velz))* (*Bx) + + ((*gxy) * (*velx) + (*gyy) * (*vely) + (*gyz) * (*velz))* (*By) + + ((*gxz) * (*velx) + (*gyz) * (*vely) + (*gzz) * (*velz))* (*Bz); + + lg.VdotBsq = VdotB*VdotB; + lg.QdotBsq = QdotB*QdotB ; lg.Qdotn = -(tau + dens) ; @@ -338,6 +362,8 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt) ( lg.D = dens; lg.half_Bsq = 0.5*lg.Bsq; + + lg.tau = tau; /* calculate W from last timestep and use for guess */ vsq = @@ -354,7 +380,8 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt) ( } if(vsq < 0. || vsq > 1. ) { *retval = 2.; - fprintf(stdout," *retval = %26.16e \n", *retval ); + fprintf(stdout," *retval (inside loop) = %26.16e \n", *retval ); + fflush(stdout); return; } @@ -368,18 +395,86 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt) ( CCTK_REAL uold = u; CCTK_REAL dum1,dum2; +// eosvars.eosanyerr[0]=5; p = pressure_rho0_eps_eosomni(rho0,*epsilon,&dum1,&dum2,&eosvars) ; // EOSOMNI - CCTK_REAL p2 = pressure_rho0_u(rho0,u,gammaeos) ; // EOS + // p = pressure_rho0_u(rho0,u,gammaeos) ; // EOS + + w = rho0 + u + p; - if (fabs(p-p2) > 1.0e-8){ - printf("p = %20.16g, p2 = %20.16g, delta pressure = %20.16g \n",p, p2, p-p2); - } + // we need this for the hot EOS failsafe tau reset later - w = rho0 + u + p ; + lg.bsq = lg.Bsq/(w*w) + lg.VdotBsq; + + W_last = w*gammasq; + +// fprintf(stdout," anyerr[0] (before loop) = %d \n", eosvars.eosanyerr[0]); +// fprintf(stdout,"vsq: rho0=%f, epsilon=%f, temp=%f, ye=%f\n", rho0,epsilon,eosvars->eos_temp[0],eosvars->eos_y_e[0]); +// fprintf(stdout,"rho0=%f, epsilon=%f, temp=%f, ye=%f\n", rho0,*epsilon,eosvars.eos_temp[0],eosvars.eos_y_e[0]); +// fprintf(stdout,"keyerr=%d, anyerr=%d \n", eosvars.eoskeyerr[0], eosvars.eosanyerr[0]); +// #pragma omp critical + if (eosvars.eosanyerr[0] != 0) { + fprintf(stdout,"inside 1: anyerr=%d\n",eosvars.eosanyerr[0]); + if (eosvars.eosrl >= GRHydro_c2p_warn_from_reflevel) { + fprintf(stdout,"inside 2: anyerr=%d\n",eosvars.eosanyerr[0]); + if(eosvars.eoskeyerr[0] == 667 && GRHydro_eos_hot_eps_fix != 0) { + fprintf(stdout,"inside 3: anyerr=%d, keyerr=%d\n",eosvars.eosanyerr[0],eosvars.eoskeyerr[0]); + // Handling of the case when no new temperature can be + // found for a given epsilon. The amount of times + // we get to this place should be monitored, as this + // 'failsafe' mode creates artificial heat. + + nf = 0; + while(eosvars.eosanyerr[0] != 0 && nf<= nfudgemax) { + eosvars.eosanyerr[0] = 0; + tau = tau + rho0*abs((*epsilon))*0.05*w*w; + *epsilon = (*epsilon) + abs((*epsilon))*0.05; + u = (*epsilon) * rho0; + EOS_Omni_press(eosvars.eoshandle,eosvars.eoskeytemp,eosvars.eosprec,1, + &rho0,epsilon,eosvars.eos_temp, + eosvars.eos_y_e,&p,eosvars.eoskeyerr,eosvars.eosanyerr); + nf = nf + 1; + } + fprintf(stdout, "Iterations of heat injection: %d \n",nf); + fprintf(stdout, "it=%d, rl=%d, i=%d, j=%d, k=%d, x=%f,y=%f,z=%f\n",eosvars.eosit,eosvars.eosrl,eosvars.eosi,eosvars.eosj,eosvars.eosk,eosvars.eosx,eosvars.eosy,eosvars.eosz); + if(nf >= nfudgemax) { + if(GRHydro_c2p_reset_eps_tau_hot_eos != 1) { + CCTK_WARN(failinfomode,"EOS error in c2p 0: injected heat too many times"); + fprintf(stdout, "it=%d, rl=%d, i=%d, j=%d, k=%d, x=%f,y=%f,z=%f\n",eosvars.eosit,eosvars.eosrl,eosvars.eosi,eosvars.eosj,eosvars.eosk,eosvars.eosx,eosvars.eosy,eosvars.eosz); + //CCTK_WARN(failinfomode,warnline); + fprintf(stdout,"rho=%f, epsilon=%f, temp=%f, ye=%f\n", rho0,*epsilon,eosvars.eos_temp[0],eosvars.eos_y_e[0]); + //call CCTK_WARN(failinfomode,warnline) + fprintf(stdout,"code=%d\n",eosvars.eoskeyerr[0]); + //call CCTK_WARN(failinfomode,warnline) + CCTK_WARN(failwarnmode,"Aborting!!!"); + } + else { + CCTK_WARN(failinfomode,"EOS error in c2p 0: LAST RESORT -- reset eps and tau based on old temp"); + fprintf(stdout, "it=%d, rl=%d, i=%d, j=%d, k=%d, x=%f,y=%f,z=%f\n",eosvars.eosit,eosvars.eosrl,eosvars.eosi,eosvars.eosj,eosvars.eosk,eosvars.eosx,eosvars.eosy,eosvars.eosz); + fprintf(stdout,"rho=%f, epsilon=%f, temp=%f, ye=%f\n", rho0,*epsilon,eosvars.eos_temp[0],eosvars.eos_y_e[0]); + eosvars.eoskeytemp =1; + EOS_Omni_press(eosvars.eoshandle,eosvars.eoskeytemp,eosvars.eosprec,1, + &rho0,epsilon,eosvars.eos_temp,eosvars.eos_y_e,&p,eosvars.eoskeyerr,eosvars.eosanyerr); + tau = (rho0 + rho0*(*epsilon) + lg.bsq + p)*w*w - (p + lg.bsq/2.0) - w*w*lg.VdotBsq - rho0*w; + u = (*epsilon) * rho0; + //tau = ( (rho + rho*epsilon) + xpress ) * w_lorentz**2 - xpress - rho*w_lorentz; + fprintf(stdout,"rho=%f, epsilon=%f, temp=%f, ye=%f\n", rho0,*epsilon,eosvars.eos_temp[0],eosvars.eos_y_e[0]); + eosvars.eoskeytemp =0; + } + } + } + else { + CCTK_WARN(failinfomode,"EOS error in c2p 0"); + fprintf(stdout, "it=%d, rl=%d, i=%d, j=%d, k=%d, x=%f,y=%f,z=%f\n",eosvars.eosit,eosvars.eosrl,eosvars.eosi,eosvars.eosj,eosvars.eosk,eosvars.eosx,eosvars.eosy,eosvars.eosz); + fprintf(stdout,"rho=%f, dens=%f, epsilon=%f, temp=%f, ye=%f\n", rho0,dens,*epsilon,eosvars.eos_temp[0],eosvars.eos_y_e[0]); + fprintf(stdout,"code=%d\n",eosvars.eoskeyerr[0]); + CCTK_WARN(failwarnmode,"Aborting!!!"); + } + } + } - W_last = w*gammasq ; + lg.tau = tau; - //fprintf(stdout," p = %26.16e \n", p ); + //fprintf(stdout," p = %26.16e \n", p ); // Make sure that W is large enough so that v^2 < 1 : i_increase = 0; @@ -410,16 +505,51 @@ 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 ); - fprintf(stdout," x_3d[%d] = %26.16e \n", 0, x_3d[0] ); - fprintf(stdout," x_3d[%d] = %26.16e \n", 1, x_3d[1] ); - fprintf(stdout," x_3d[%d] = %26.16e \n", 2, x_3d[2] ); + fprintf(stdout," *retval = %26.16e \n", *retval ); + fprintf(stdout," *dens = %26.16e \n", *dens_in ); + fprintf(stdout," *sx = %26.16e \n", *sx_in ); + fprintf(stdout," *sy = %26.16e \n", *sy_in ); + fprintf(stdout," *sz = %26.16e \n", *sz_in ); + fprintf(stdout," *tau = %26.16e \n", *tau_in ); + fprintf(stdout," *Bconsx = %26.16e \n", *Bconsx_in ); + fprintf(stdout," *Bconsy = %26.16e \n", *Bconsy_in ); + fprintf(stdout," *Bconsz = %26.16e \n", *Bconsz_in ); + fprintf(stdout," *rho = %26.16e \n", *rho ); + fprintf(stdout," *velx = %26.16e \n", *velx ); + 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 ); + fprintf(stdout," *Bz = %26.16e \n", *Bz ); + fprintf(stdout," *bsq = %26.16e \n", *bsq ); + fprintf(stdout," *w_lorentz = %26.16e \n", *w_lorentz ); + fprintf(stdout," *gxx = %26.16e \n", *gxx ); + fprintf(stdout," *gxy = %26.16e \n", *gxy ); + fprintf(stdout," *gxz = %26.16e \n", *gxz ); + fprintf(stdout," *gyy = %26.16e \n", *gyy ); + fprintf(stdout," *gyz = %26.16e \n", *gyz ); + fprintf(stdout," *gzz = %26.16e \n", *gzz ); + fprintf(stdout," *uxx = %26.16e \n", *uxx ); + fprintf(stdout," *uxy = %26.16e \n", *uxy ); + fprintf(stdout," *uxz = %26.16e \n", *uxz ); + fprintf(stdout," *uyy = %26.16e \n", *uyy ); + fprintf(stdout," *uyz = %26.16e \n", *uyz ); + fprintf(stdout," *uzz = %26.16e \n", *uzz ); + fprintf(stdout," *det = %26.16e \n", *det ); + fprintf(stdout," *epsnegative = %10d \n", *epsnegative ); + fprintf(stdout," *retval = %26.16e \n", *retval ); + fflush(stdout); return; } else{ if(W <= 0. || W > W_TOO_BIG) { *retval = 3.; fprintf(stdout," *retval = %26.16e \n", *retval ); + fflush(stdout); return; } } @@ -428,6 +558,7 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt) ( if( vsq >= 1. ) { *retval = 4.; fprintf(stdout," *retval = %26.16e \n", *retval ); + fflush(stdout); return; } @@ -445,9 +576,80 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt) ( } else { u=x_3d[2]; *epsilon = u/rho0; + + // Update tau in case we have changed it during the NR iterations + + tau = lg.tau; + +// eosvars.eosanyerr[0]=5; + p = pressure_rho0_eps_eosomni(rho0,*epsilon,&dum1,&dum2,&eosvars) ; // EOSOMNI - // printf("%g %g %g %g\n",rho0,u,*epsilon,p); - } + +// #pragma omp critical + if (eosvars.eosanyerr[0] != 0) { + fprintf(stdout,"inside 4: anyerr=%d\n",eosvars.eosanyerr[0]); + if (eosvars.eosrl >= GRHydro_c2p_warn_from_reflevel) { + if(eosvars.eoskeyerr[0] == 667 && GRHydro_eos_hot_eps_fix != 0) { + // Handling of the case when no new temperature can be + // found for a given epsilon. The amount of times + // we get to this place should be monitored, as this + // 'failsafe' mode creates artificial heat. + nf = 0; + while(eosvars.eosanyerr[0] != 0 && nf<= nfudgemax) { + eosvars.eosanyerr[0] = 0; + tau = tau + rho0*abs((*epsilon))*0.05*w*w; + *epsilon = (*epsilon) + abs((*epsilon))*0.05; + u = (*epsilon) * rho0; + EOS_Omni_press(eosvars.eoshandle,eosvars.eoskeytemp,eosvars.eosprec,1, + &rho0,epsilon,eosvars.eos_temp, + eosvars.eos_y_e,&p,eosvars.eoskeyerr,eosvars.eosanyerr); + nf = nf + 1; + } + fprintf(stdout, "Iterations of heat injection: %d \n",nf); + fprintf(stdout, "it=%d, rl=%d, i=%d, j=%d, k=%d, x=%f,y=%f,z=%f\n",eosvars.eosit,eosvars.eosrl,eosvars.eosi,eosvars.eosj,eosvars.eosk,eosvars.eosx,eosvars.eosy,eosvars.eosz); + if(nf >= nfudgemax) { + if(GRHydro_c2p_reset_eps_tau_hot_eos != 1) { + CCTK_WARN(failinfomode,"EOS error in c2p 0: injected heat too many times"); + fprintf(stdout, "it=%d, rl=%d, i=%d, j=%d, k=%d, x=%f,y=%f,z=%f\n",eosvars.eosit,eosvars.eosrl,eosvars.eosi,eosvars.eosj,eosvars.eosk,eosvars.eosx,eosvars.eosy,eosvars.eosz); + //CCTK_WARN(failinfomode,warnline); + fprintf(stdout,"rho=%f, epsilon=%f, temp=%f, ye=%f\n", rho0,*epsilon,eosvars.eos_temp[0],eosvars.eos_y_e[0]); + //call CCTK_WARN(failinfomode,warnline) + fprintf(stdout,"code=%d\n",eosvars.eoskeyerr[0]); + //call CCTK_WARN(failinfomode,warnline) + CCTK_WARN(failwarnmode,"Aborting!!!"); + } + else { + CCTK_WARN(failinfomode,"EOS error in c2p 0: LAST RESORT -- reset eps and tau based on old temp"); + fprintf(stdout, "it=%d, rl=%d, i=%d, j=%d, k=%d, x=%f,y=%f,z=%f\n",eosvars.eosit,eosvars.eosrl,eosvars.eosi,eosvars.eosj,eosvars.eosk,eosvars.eosx,eosvars.eosy,eosvars.eosz); + fprintf(stdout,"rho=%f, epsilon=%f, temp=%f, ye=%f\n", rho0,*epsilon,eosvars.eos_temp[0],eosvars.eos_y_e[0]); + eosvars.eoskeytemp =1; + EOS_Omni_press(eosvars.eoshandle,eosvars.eoskeytemp,eosvars.eosprec,1, + &rho0,epsilon,eosvars.eos_temp,eosvars.eos_y_e,&p,eosvars.eoskeyerr,eosvars.eosanyerr); + // tau = (rho0 + rho0*(*epsilon) + lg.bsq + p)*w*w - (p + lg.bsq/2.0) - rho0*w; + tau = (rho0 + rho0*(*epsilon) + lg.bsq + p)*w*w - (p + lg.bsq/2.0) - w*w*lg.VdotBsq - rho0*w; + u = (*epsilon) * rho0; + //tau = ( (rho + rho*epsilon) + xpress ) * w_lorentz**2 - xpress - rho*w_lorentz; + //fprintf(stdout,"rho=%f, epsilon=%f, temp=%f, ye=%f\n", rho0,*epsilon,eosvars.eos_temp[0],eosvars.eos_y_e[0]); + eosvars.eoskeytemp =0; + } + } + } + else { + CCTK_WARN(failinfomode,"EOS error in c2p 0"); + fprintf(stdout, "it=%d, rl=%d, i=%d, j=%d, k=%d, x=%f,y=%f,z=%f\n",eosvars.eosit,eosvars.eosrl,eosvars.eosi,eosvars.eosj,eosvars.eosk,eosvars.eosx,eosvars.eosy,eosvars.eosz); + fprintf(stdout,"rho=%f, dens=%f, epsilon=%f, temp=%f, ye=%f\n", rho0,dens,*epsilon,eosvars.eos_temp[0],eosvars.eos_y_e[0]); + fprintf(stdout,"code=%d\n",eosvars.eoskeyerr[0]); + CCTK_WARN(0,"Aborting!!!"); + } + } + } + + // printf("%g %g %g %g\n",rho0,u,*epsilon,p); + } + + // update tau in case we have modified it in C2P + + *tau_in = tau*sqrt_detg; // 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: @@ -455,11 +657,12 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt) ( 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 ); + 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 ); + fflush(stdout); return; } @@ -482,8 +685,48 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt) ( *velz = 0.0; *w_lorentz = 1.0; } - #if(DEBUG_CON2PRIMM) + fprintf(stdout," *retvalo = %26.16e \n", *retval ); + fprintf(stdout," *x2 = %26.16e \n", *x ); + fprintf(stdout," *y2 = %26.16e \n", *y ); + fprintf(stdout," *z2 = %26.16e \n", *z ); + fprintf(stdout," *denso = %26.16e \n", *dens_in ); + fprintf(stdout," *sxo = %26.16e \n", *sx_in ); + fprintf(stdout," *syo = %26.16e \n", *sy_in ); + fprintf(stdout," *szo = %26.16e \n", *sz_in ); + fprintf(stdout," *tauo = %26.16e \n", *tau_in ); + fprintf(stdout," *Bconsxo = %26.16e \n", *Bconsx_in ); + fprintf(stdout," *Bconsyo = %26.16e \n", *Bconsy_in ); + fprintf(stdout," *Bconszo = %26.16e \n", *Bconsz_in ); + fprintf(stdout," *rhoo = %26.16e \n", *rho ); + fprintf(stdout," *velxo = %26.16e \n", *velx ); + fprintf(stdout," *velyo = %26.16e \n", *vely ); + fprintf(stdout," *velzo = %26.16e \n", *velz ); + fprintf(stdout," *epsilono = %26.16e \n", *epsilon ); + fprintf(stdout," *temp_ino = %26.16e \n", *temp_in ); + fprintf(stdout," *y_e_ino = %26.16e \n", *y_e_in ); + fprintf(stdout," *pressureo = %26.16e \n", *pressure ); + fprintf(stdout," *Bxo = %26.16e \n", *Bx ); + fprintf(stdout," *Byo = %26.16e \n", *By ); + fprintf(stdout," *Bzo = %26.16e \n", *Bz ); + fprintf(stdout," *bsqo = %26.16e \n", *bsq ); + fprintf(stdout," *w_lorentzo = %26.16e \n", *w_lorentz ); + fprintf(stdout," *gxxo = %26.16e \n", *gxx ); + fprintf(stdout," *gxyo = %26.16e \n", *gxy ); + fprintf(stdout," *gxzo = %26.16e \n", *gxz ); + fprintf(stdout," *gyyo = %26.16e \n", *gyy ); + fprintf(stdout," *gyzo = %26.16e \n", *gyz ); + fprintf(stdout," *gzzo = %26.16e \n", *gzz ); + fprintf(stdout," *uxxo = %26.16e \n", *uxx ); + fprintf(stdout," *uxyo = %26.16e \n", *uxy ); + fprintf(stdout," *uxzo = %26.16e \n", *uxz ); + fprintf(stdout," *uyyo = %26.16e \n", *uyy ); + fprintf(stdout," *uyzo = %26.16e \n", *uyz ); + fprintf(stdout," *uzzo = %26.16e \n", *uzz ); + fprintf(stdout," *deto = %26.16e \n", *det ); + fprintf(stdout," *epsnegativeo = %10d \n", *epsnegative ); + + fprintf(stdout,"rho = %26.16e \n",*rho ); fprintf(stdout,"epsilon = %26.16e \n",*epsilon ); fprintf(stdout,"pressure = %26.16e \n",*pressure ); @@ -759,13 +1002,13 @@ static CCTK_INT threed_newton_raphson_omni( CCTK_REAL x[], struct eosomnivars *e /* before stopping */ /*****************************************************************************/ - if( (fabs(errx) <= NEWT_TOL) && (doing_extra == 0) && (EXTRA_NEWT_ITER > 0) ) { + if( (fabs(errx) <= eosvars->prec) && (doing_extra == 0) && (EXTRA_NEWT_ITER > 0) ) { doing_extra = 1; } if( doing_extra == 1 ) i_extra++ ; - if( ((fabs(errx) <= NEWT_TOL)&&(doing_extra == 0)) + if( ((fabs(errx) <= eosvars->prec)&&(doing_extra == 0)) || (i_extra > EXTRA_NEWT_ITER) || (n_iter >= (MAX_NEWT_ITER-1)) ) { keep_iterating = 0; } @@ -781,10 +1024,10 @@ static CCTK_INT threed_newton_raphson_omni( CCTK_REAL x[], struct eosomnivars *e } - if( fabs(errx) <= NEWT_TOL ){ + if( fabs(errx) <= eosvars->prec ){ return(0); } - else if( (fabs(errx) <= MIN_NEWT_TOL) && (fabs(errx) > NEWT_TOL) ){ + else if( (fabs(errx) <= eosvars->prec) && (fabs(errx) > eosvars->prec) ){ return(0); } else { @@ -904,23 +1147,94 @@ static void func_vsq_eosomni(CCTK_REAL x[], CCTK_REAL dx[], CCTK_REAL resid[], struct eosomnivars *eosvars, struct LocGlob *lgp) { - CCTK_REAL W, vsq, u, p_tmp,epsilon,Wsq, dPdvsq, dPdW; - CCTK_REAL res0, LorInv, rho0, Winv, QB2Winv2, drho0_dv, res1, j11, detJ,detJinv; + CCTK_REAL p_tmp; + CCTK_REAL W, vsq, u, epsilon,Wsq, dPdvsq, dPdW; + CCTK_REAL res0, Lor, LorInv, rho0, Winv, QB2Winv2, drho0_dv, res1, j11, detJ,detJinv; CCTK_REAL QB2Winv3, j10,B2plusW,detJ_gcf,t36,B2plusW_sq, j00,j01,j12,j20,j21,j22; CCTK_REAL dpress_dv,dpress_du,dpdrho,dpdeps,c00,c01,c02,c10,c11,c12,c20,c21,c22; + CCTK_INT nf, nfudgemax, failwarnmode, failinfomode; + DECLARE_CCTK_PARAMETERS; + + nf =0; + nfudgemax =30; + failinfomode =2; + failwarnmode =1; W = x[0]; vsq = x[1]; u = x[2]; Wsq = W*W; - + Lor = 1.0/sqrt(1.0-vsq); LorInv = sqrt(1.0-vsq); rho0 = lgp->D * LorInv; epsilon = u/rho0; p_tmp = pressure_rho0_eps_eosomni(rho0,epsilon,&dpdrho,&dpdeps,eosvars); + +// fprintf(stdout,"vsq: keyerr=%d, eosrl = %d, anyerr = %d \n",eosvars->eoskeyerr[0],eosvars->eosrl, eosvars->eosanyerr[0]); +// fprintf(stdout,"vsq: rho0=%f, epsilon=%f, temp=%f, ye=%f\n", rho0,epsilon,eosvars->eos_temp[0],eosvars->eos_y_e[0]); + #pragma omp critical + if (eosvars->eosanyerr[0] != 0) { + fprintf(stdout,"inside vsq_omni: anyerr=%d rl=%d\n",eosvars->eosanyerr[0], eosvars->eosrl); + if (eosvars->eosrl >= GRHydro_c2p_warn_from_reflevel) { + if(eosvars->eoskeyerr[0] == 667 && GRHydro_eos_hot_eps_fix != 0) { + // Handling of the case when no new temperature can be + // found for a given epsilon. The amount of times + // we get to this place should be monitored, as this + // 'failsafe' mode creates artificial heat. + nf = 0; + while(eosvars->eosanyerr[0] != 0 && nf<= nfudgemax) { + eosvars->eosanyerr[0] = 0; + lgp->tau = lgp->tau + rho0*abs(epsilon)*0.05*Lor*Lor; + epsilon = epsilon + abs(epsilon)*0.05; + u = epsilon*rho0; + CCTK_WARN(failwarnmode,"Injecting heat inside vsq_omni without resetting tau"); + fprintf(stdout,"Injecting heat inside vsq_omni without resetting tau: anyerr=%d\n",eosvars->eosanyerr[0]); + EOS_Omni_press(eosvars->eoshandle,eosvars->eoskeytemp,eosvars->eosprec,1, + &rho0,&epsilon,eosvars->eos_temp, + eosvars->eos_y_e,&p_tmp,eosvars->eoskeyerr,eosvars->eosanyerr); + nf = nf + 1; + } + fprintf(stdout, "Iterations of heat injection: %d \n",nf); + fprintf(stdout, "it=%d, rl=%d, i=%d, j=%d, k=%d, x=%f,y=%f,z=%f\n",eosvars->eosit,eosvars->eosrl,eosvars->eosi,eosvars->eosj,eosvars->eosk,eosvars->eosx,eosvars->eosy,eosvars->eosz); + if(nf >= nfudgemax) { + if(GRHydro_c2p_reset_eps_tau_hot_eos != 1) { + CCTK_WARN(failinfomode,"EOS error in c2p 0: injected heat too many times"); + fprintf(stdout, "it=%d, rl=%d, i=%d, j=%d, k=%d, x=%f,y=%f,z=%f\n",eosvars->eosit,eosvars->eosrl,eosvars->eosi,eosvars->eosj,eosvars->eosk,eosvars->eosx,eosvars->eosy,eosvars->eosz); + //CCTK_WARN(failinfomode,warnline); + fprintf(stdout,"rho=%f, epsilon=%f, temp=%f, ye=%f\n", rho0,epsilon,eosvars->eos_temp[0],eosvars->eos_y_e[0]); + //call CCTK_WARN(failinfomode,warnline) + fprintf(stdout,"code=%d\n",eosvars->eoskeyerr[0]); + //call CCTK_WARN(failinfomode,warnline) + CCTK_WARN(failwarnmode,"Aborting!!!"); + } + else { + CCTK_WARN(failinfomode,"EOS error in c2p 0: LAST RESORT -- reset eps and tau based on old temp"); + fprintf(stdout, "it=%d, rl=%d, i=%d, j=%d, k=%d, x=%f,y=%f,z=%f\n",eosvars->eosit,eosvars->eosrl,eosvars->eosi,eosvars->eosj,eosvars->eosk,eosvars->eosx,eosvars->eosy,eosvars->eosz); + fprintf(stdout,"rho=%f, epsilon=%f, temp=%f, ye=%f\n", rho0,epsilon,eosvars->eos_temp[0],eosvars->eos_y_e[0]); + eosvars->eoskeytemp =1; + EOS_Omni_press(eosvars->eoshandle,eosvars->eoskeytemp,eosvars->eosprec,1, + &rho0,&epsilon,eosvars->eos_temp,eosvars->eos_y_e,&p_tmp,eosvars->eoskeyerr,eosvars->eosanyerr); + //tau = (rho0 + rho0*(*epsilon) + lg.bsq + p)*w*w - (p + lg.bsq/2.0) - w*w*lg.VdotBsq - rho0*w; + lgp->tau = (rho0 + rho0*epsilon + lgp->bsq + p_tmp)*Lor*Lor - (p_tmp + lgp->bsq/2.0) - Lor*Lor*lgp->VdotBsq - rho0*Lor; + u = epsilon*rho0; + eosvars->eoskeytemp = 0; + fprintf(stdout,"rho=%f, epsilon=%f, temp=%f, ye=%f\n", rho0,epsilon,eosvars->eos_temp[0],eosvars->eos_y_e[0]); + } + } + } + else { + CCTK_WARN(failinfomode,"EOS error in c2p 0"); + fprintf(stdout, "it=%d, rl=%d, i=%d, j=%d, k=%d, x=%f,y=%f,z=%f\n",eosvars->eosit,eosvars->eosrl,eosvars->eosi,eosvars->eosj,eosvars->eosk,eosvars->eosx,eosvars->eosy,eosvars->eosz); + fprintf(stdout,"rho=%f, epsilon=%f, temp=%f, ye=%f\n", rho0,epsilon,eosvars->eos_temp[0],eosvars->eos_y_e[0]); + fprintf(stdout,"code=%d\n",eosvars->eoskeyerr[0]); + CCTK_WARN(failwarnmode,"Aborting!!!"); + } + } + } + lgp->Qdotn = -(lgp->tau + lgp->D); B2plusW = lgp->Bsq+W; B2plusW_sq = B2plusW*B2plusW; Winv = 1/W; diff --git a/src/GRHydro_Con2PrimM_pt_EOSOmniold.c b/src/GRHydro_Con2PrimM_pt_EOSOmniold.c new file mode 100644 index 0000000..0eb4845 --- /dev/null +++ b/src/GRHydro_Con2PrimM_pt_EOSOmniold.c @@ -0,0 +1,1056 @@ +/*********************************************************************************** + Copyright 2006 Scott C. Noble, Charles F. Gammie, Jonathan C. McKinney, + and Luca Del Zanna. + + PVS_GRMHD + + This file was derived from PVS_GRMHD. The authors of PVS_GRMHD include + Scott C. Noble, Charles F. Gammie, Jonathan C. McKinney, and Luca Del Zanna. + PVS_GRMHD is available under the GPL from: + http://rainman.astro.uiuc.edu/codelib/ + + You are morally obligated to cite the following paper in his/her + scientific literature that results from use of this file: + + [1] Noble, S. C., Gammie, C. F., McKinney, J. C., \& Del Zanna, L. \ 2006, + Astrophysical Journal, 641, 626. + + PVS_GRMHD is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + PVS_GRMHD is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with PVS_GRMHD; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + + If the user has any questions, please direct them to Scott C. Noble at + scn@astro.rit.edu . + +***********************************************************************************/ + + + +#include <stdio.h> +#include <stdlib.h> +#include <stdarg.h> +#include <string.h> +#include <math.h> +#include <float.h> +#include <complex.h> + +#include "cctk.h" +#include "cctk_Parameters.h" + +/* Set this to be 1 if you want debug output */ +#define DEBUG_CON2PRIMM (0) + + +/* Adiabatic index used for the state equation */ + +#define MAX_NEWT_ITER (30) /* Max. # of Newton-Raphson iterations for find_root_2D(); */ +#define NEWT_TOL (1.0e-10) /* Min. of tolerance allowed for Newton-Raphson iterations */ +#define MIN_NEWT_TOL (1.0e-10) /* Max. of tolerance allowed for Newton-Raphson iterations */ +#define EXTRA_NEWT_ITER (2) + +#define NEWT_TOL2 (1.0e-15) /* TOL of new 1D^*_{v^2} gnr2 method */ +#define MIN_NEWT_TOL2 (1.0e-10) /* TOL of new 1D^*_{v^2} gnr2 method */ + +#define W_TOO_BIG (1.e20) /* \gamma^2 (\rho_0 + u + p) is assumed + to always be smaller than this. This + is used to detect solver failures */ + +#define FAIL_VAL (1.e30) /* Generic value to which we set variables when a problem arises */ + +/************************************************** + The following functions assume a Gamma-law EOS: +***************************************************/ + +/* Local Globals */ +struct LocGlob { + CCTK_REAL Bsq,QdotBsq,Qtsq,Qdotn,D,half_Bsq ; +} ; + +struct eosomnivars { + CCTK_INT eoshandle,eoskeytemp; + CCTK_REAL prec,eosprec,eos_y_e[1],eos_temp[1]; + CCTK_INT eoskeyerr[1],eosanyerr[1]; +} ; + +// Declarations: +static CCTK_REAL vsq_calc(CCTK_REAL W, struct LocGlob *lgp); + +static CCTK_INT twod_newton_raphson( CCTK_REAL x[], CCTK_REAL gammaeos, struct LocGlob *lgp, + void (*funcd) (CCTK_REAL [], CCTK_REAL [], + CCTK_REAL [], CCTK_REAL [][2], + CCTK_REAL *, CCTK_REAL *, CCTK_REAL, struct LocGlob *) ); + +static CCTK_INT threed_newton_raphson_omni( CCTK_REAL x[], struct eosomnivars *eosvars, struct LocGlob *lgp, + void (*funcd) (CCTK_REAL [], CCTK_REAL [], + CCTK_REAL [], CCTK_REAL [][3], + CCTK_REAL *, CCTK_REAL *, + struct eosomnivars *eosvars, struct LocGlob *) ); + +static void func_vsq( CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], CCTK_REAL [][2], + CCTK_REAL *f, CCTK_REAL *df, CCTK_REAL gammaeos, struct LocGlob *lgp); + +static void func_vsq_eosomni( CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], CCTK_REAL [][3], + CCTK_REAL *f, CCTK_REAL *df, struct eosomnivars *eosvars, struct LocGlob *lgp); + +static CCTK_REAL x1_of_x0(CCTK_REAL x0, struct LocGlob *lgp ) ; + + +// EOS STUFF: +static CCTK_REAL eos_info(CCTK_REAL W, CCTK_REAL vsq, CCTK_REAL *dpdw, CCTK_REAL *dpdvsq, CCTK_REAL gammaeos, struct LocGlob *lgp); +static CCTK_REAL eos_info_eosomni(CCTK_REAL W, CCTK_REAL vsq, CCTK_REAL eps0, struct LocGlob *lgp); + +/* pressure as a function of rho0 and u */ +//static CCTK_REAL pressure_rho0_u(CCTK_REAL rho0, CCTK_REAL u, CCTK_REAL gammaeos) +//{ +// return((gammaeos - 1.)*u) ; +//} + +static CCTK_REAL pressure_rho0_eps_eosomni(CCTK_REAL rho0, CCTK_REAL eps, CCTK_REAL* dpdrho, CCTK_REAL* dpdeps, struct eosomnivars *eosvars); + +/* Pressure as a function of rho0 and w = rho0 + u + p */ +static CCTK_REAL pressure_rho0_w(CCTK_REAL rho0, CCTK_REAL w,CCTK_REAL gammaeos) +{ + return((gammaeos-1.)*(w - rho0)/gammaeos) ; +} + +void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt2) ( + CCTK_INT *handle, CCTK_INT *keytemp, CCTK_REAL *eos_prec, CCTK_REAL *prec, + CCTK_REAL *gamma_eos, + CCTK_REAL *dens_in, + CCTK_REAL *sx_in, CCTK_REAL *sy_in, CCTK_REAL *sz_in, + CCTK_REAL *tau_in, CCTK_REAL *Bconsx_in, CCTK_REAL *Bconsy_in, CCTK_REAL *Bconsz_in, + CCTK_REAL *y_e_in, CCTK_REAL *temp_in, + CCTK_REAL *rho, + CCTK_REAL *velx, CCTK_REAL *vely, CCTK_REAL *velz, + CCTK_REAL *epsilon, CCTK_REAL *pressure, + CCTK_REAL *Bx, CCTK_REAL *By, CCTK_REAL *Bz, CCTK_REAL *bsq, + CCTK_REAL *w_lorentz, + CCTK_REAL *gxx, CCTK_REAL *gxy, CCTK_REAL *gxz, + CCTK_REAL *gyy, CCTK_REAL *gyz, CCTK_REAL *gzz, + CCTK_REAL *uxx, CCTK_REAL *uxy, CCTK_REAL *uxz, + CCTK_REAL *uyy, CCTK_REAL *uyz, CCTK_REAL *uzz, + CCTK_REAL *det, + CCTK_INT *epsnegative, + CCTK_REAL *retval); + +/**********************************************************************/ +/********************************************************************************** + + Con2PrimM_pt(): + ----------------------------- + + -- Attempts an inversion from GRMHD conserved variables to primitive variables assuming a guess. + + -- Uses the 2D method of Noble et al. (2006): + -- Solves for two independent variables (W,v^2) via a 2D + Newton-Raphson method + -- Can be used (in principle) with a general equation of state. + + -- Minimizes two residual functions using a homemade Newton-Raphson routine. + -- It is homemade so that it can catch exceptions and handle them correctly, plus it is + optimized for this particular problem. + + -- Note that the notation used herein is that of Noble et al. (2006) except for the argument + list. + + +INPUT: (using GRHydro variable defintions) + + s[x,y,z] = scons[0,1,2] = \alpha \sqrt(\gamma) T^0_i + dens, tau = as defined in GRHydro and are assumed to be densitized (i.e. with sqrt(\gamma)) + dens = D = \sqrt(\gamma) W \rho + tau = \alpha^2 \sqrt(\gamma) T^{00} - D + g[x,y,z][x,y,x] = spatial metric corresponding to \gamma + u[x,y,z][x,y,z] = inverse of the spatial metric, g[x,y,z][x,y,x] + det = sqrt(\gamma) + B[x,y,z] = Bvec[0,1,2] + bsq = b^\mu b_\mu + + epsnegative = (integer) + = 0 if rho and epsilon are positive + != 0 otherwise + + + -- (users should set B[x,y,z] = 0 for hydrodynamic runs) + + +OUTPUT: (using GRHydro variable defintions) + rho, eps = as defined in GRHydro, primitive variables + vel[x,y,z] = as defined in GRHydro, primitive variables + + +RETURN VALUE: of retval = (i*100 + j) where + i = 0 -> Newton-Raphson solver either was not called (yet or not used) + or returned successfully; + 1 -> Newton-Raphson solver did not converge to a solution with the + given tolerances; + 2 -> Newton-Raphson procedure encountered a numerical divergence + (occurrence of "nan" or "+/-inf" ; + + j = 0 -> success + 1 -> failure: some sort of failure in Newton-Raphson; + 2 -> failure: unphysical vsq = v^2 value at initial guess; + 3 -> failure: W<0 or W>W_TOO_BIG + 4 -> failure: v^2 > 1 + ( used to be 5 -> failure: rho,uu <= 0 but now sets epsnegative to non-zero ) + +**********************************************************************************/ + + +void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt2) ( + CCTK_INT *handle, CCTK_INT *keytemp, CCTK_REAL *eos_prec, CCTK_REAL *prec, + CCTK_REAL *gamma_eos, + CCTK_REAL *dens_in, + CCTK_REAL *sx_in, CCTK_REAL *sy_in, CCTK_REAL *sz_in, + CCTK_REAL *tau_in, + CCTK_REAL *Bconsx_in, CCTK_REAL *Bconsy_in, CCTK_REAL *Bconsz_in, + CCTK_REAL *y_e_in, CCTK_REAL* temp_in, + CCTK_REAL *rho, + CCTK_REAL *velx, CCTK_REAL *vely, CCTK_REAL *velz, + CCTK_REAL *epsilon, CCTK_REAL *pressure, + CCTK_REAL *Bx, CCTK_REAL *By, CCTK_REAL *Bz, + CCTK_REAL *bsq, + CCTK_REAL *w_lorentz, + CCTK_REAL *gxx, CCTK_REAL *gxy, CCTK_REAL *gxz, + CCTK_REAL *gyy, CCTK_REAL *gyz, CCTK_REAL *gzz, + CCTK_REAL *uxx, CCTK_REAL *uxy, CCTK_REAL *uxz, + CCTK_REAL *uyy, CCTK_REAL *uyz, CCTK_REAL *uzz, + CCTK_REAL *det, + CCTK_INT *epsnegative, + CCTK_REAL *retval) + +{ + CCTK_REAL x_3d[3]; + CCTK_REAL sx, sy, sz; + CCTK_REAL usx, usy, usz; + CCTK_REAL tau, dens, gammaeos; + CCTK_REAL QdotB; + CCTK_REAL rho0,u,p,w,gammasq,gamma,gtmp,W_last,W,vsq; + CCTK_REAL g_o_WBsq, QdB_o_W; + CCTK_REAL detg = (*det); + CCTK_REAL sqrt_detg = sqrt(detg); + CCTK_REAL 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); + eosvars.eoskeytemp = *keytemp; + eosvars.eosprec = *eos_prec; + eosvars.prec = *prec; + eosvars.eos_y_e[0] = *y_e_in; + eosvars.eos_temp[0] = *temp_in; + eosvars.eoskeyerr[0] = 0; + eosvars.eosanyerr[0] = 0; + + gammaeos = *gamma_eos; + + /* Assume ok initially: */ + *retval = 0.; + *epsnegative = 0; + +#if(DEBUG_CON2PRIMM) + fprintf(stdout," *dens = %26.16e \n", *dens_in ); + fprintf(stdout," *sx = %26.16e \n", *sx_in ); + fprintf(stdout," *sy = %26.16e \n", *sy_in ); + fprintf(stdout," *sz = %26.16e \n", *sz_in ); + fprintf(stdout," *tau = %26.16e \n", *tau_in ); + fprintf(stdout," *Bconsx = %26.16e \n", *Bconsx_in ); + fprintf(stdout," *Bconsy = %26.16e \n", *Bconsy_in ); + fprintf(stdout," *Bconsz = %26.16e \n", *Bconsz_in ); + fprintf(stdout," *rho = %26.16e \n", *rho ); + fprintf(stdout," *velx = %26.16e \n", *velx ); + 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 ); + fprintf(stdout," *Bz = %26.16e \n", *Bz ); + fprintf(stdout," *bsq = %26.16e \n", *bsq ); + fprintf(stdout," *w_lorentz = %26.16e \n", *w_lorentz ); + fprintf(stdout," *gxx = %26.16e \n", *gxx ); + fprintf(stdout," *gxy = %26.16e \n", *gxy ); + fprintf(stdout," *gxz = %26.16e \n", *gxz ); + fprintf(stdout," *gyy = %26.16e \n", *gyy ); + fprintf(stdout," *gyz = %26.16e \n", *gyz ); + fprintf(stdout," *gzz = %26.16e \n", *gzz ); + fprintf(stdout," *uxx = %26.16e \n", *uxx ); + fprintf(stdout," *uxy = %26.16e \n", *uxy ); + fprintf(stdout," *uxz = %26.16e \n", *uxz ); + fprintf(stdout," *uyy = %26.16e \n", *uyy ); + fprintf(stdout," *uyz = %26.16e \n", *uyz ); + fprintf(stdout," *uzz = %26.16e \n", *uzz ); + fprintf(stdout," *det = %26.16e \n", *det ); + fprintf(stdout," *epsnegative = %10d \n", *epsnegative ); + fprintf(stdout," *retval = %26.16e \n", *retval ); + fflush(stdout); +#endif + + /* First undensitize all conserved variables : */ + sx = ( *sx_in) * inv_sqrt_detg; + sy = ( *sy_in) * inv_sqrt_detg; + sz = ( *sz_in) * inv_sqrt_detg; + tau = ( *tau_in) * inv_sqrt_detg; + dens = (*dens_in) * inv_sqrt_detg; + + usx = (*uxx)*sx + (*uxy)*sy + (*uxz)*sz; + usy = (*uxy)*sx + (*uyy)*sy + (*uyz)*sz; + usz = (*uxz)*sx + (*uyz)*sy + (*uzz)*sz; + + *Bx = (*Bconsx_in) * inv_sqrt_detg; + *By = (*Bconsy_in) * inv_sqrt_detg; + *Bz = (*Bconsz_in) * inv_sqrt_detg; + + // Calculate various scalars (Q.B, Q^2, etc) from the conserved variables: + + lg.Bsq = + (*gxx) * (*Bx) * (*Bx) + + (*gyy) * (*By) * (*By) + + (*gzz) * (*Bz) * (*Bz) + + 2*( + (*gxy) * (*Bx) * (*By) + + (*gxz) * (*Bx) * (*Bz) + + (*gyz) * (*By) * (*Bz) ); + + QdotB = (sx * (*Bx) + sy * (*By) + sz * (*Bz)) ; + lg.QdotBsq = QdotB*QdotB ; + + lg.Qdotn = -(tau + dens) ; + + lg.Qtsq = (usx * sx + usy * sy + usz * sz) ; + + lg.D = dens; + + lg.half_Bsq = 0.5*lg.Bsq; + + /* calculate W from last timestep and use for guess */ + vsq = + (*gxx) * (*velx) * (*velx) + + (*gyy) * (*vely) * (*vely) + + (*gzz) * (*velz) * (*velz) + + 2*( + (*gxy) * (*velx) * (*vely) + + (*gxz) * (*velx) * (*velz) + + (*gyz) * (*vely) * (*velz) ); + + if( (vsq < 0.) && (fabs(vsq) < 1.0e-13) ) { + vsq = fabs(vsq); + } + if(vsq < 0. || vsq > 1. ) { + *retval = 2.; + fprintf(stdout," *retval = %26.16e \n", *retval ); + return; + } + + gammasq = 1. / (1. - vsq); + gamma = sqrt(gammasq); + + // Always calculate rho from D and gamma so that using D in EOS remains consistent + // 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 + w = rho0 + u + p ; + + 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; + while( (( W_last*W_last*W_last * ( W_last + 2.*lg.Bsq ) + - lg.QdotBsq*(2.*W_last + lg.Bsq) ) <= W_last*W_last*(lg.Qtsq-lg.Bsq*lg.Bsq)) + && (i_increase < 10) ) { + W_last *= 10.; + i_increase++; + } + + // Calculate W and vsq: + x_3d[0] = fabs( W_last ); + x_3d[1] = x1_of_x0( W_last, &lg ) ; + + //Use 2d NR for polytropes! + if (*handle==1 || *handle==2) { + *retval = 1.0*twod_newton_raphson( x_3d, gammaeos, &lg, func_vsq ) ; + + } else { + //USE 3d NR for non-polytropes! + x_3d[2] = u; + *retval = 1.0*threed_newton_raphson_omni( x_3d, &eosvars, &lg, func_vsq_eosomni ) ; + } + + W = x_3d[0]; + vsq = x_3d[1]; + + /* 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; + } + } + + // Calculate v^2: + if( vsq >= 1. ) { + *retval = 4.; + fprintf(stdout," *retval = %26.16e \n", *retval ); + return; + } + + // Recover the primitive variables from the scalars and conserved variables: + gtmp = sqrt(1. - vsq); + gamma = 1./gtmp ; + rho0 = lg.D * gtmp; + + w = W * (1. - vsq) ; + + if (*handle==1 || *handle==2) { + p = pressure_rho0_w(rho0,w,gammaeos) ; // EOS + u = w - (rho0 + p) ; + *epsilon = u / rho0; + } else { + 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); + } + + // 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; + 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; + } + + *rho = rho0; + *w_lorentz = gamma; + *pressure = p ; + + g_o_WBsq = 1./(W+lg.Bsq); + QdB_o_W = QdotB / W; + *bsq = lg.Bsq * (1.-vsq) + QdB_o_W*QdB_o_W; + + *velx = g_o_WBsq * ( usx + QdB_o_W*(*Bx) ) ; + *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 ); + fprintf(stdout,"epsilon = %26.16e \n",*epsilon ); + fprintf(stdout,"pressure = %26.16e \n",*pressure ); + fprintf(stdout,"w_lorentz = %26.16e \n",*w_lorentz); + fprintf(stdout,"bsq = %26.16e \n",*bsq ); + fprintf(stdout,"velx = %26.16e \n",*velx ); + fprintf(stdout,"vely = %26.16e \n",*vely ); + fprintf(stdout,"velz = %26.16e \n",*velz ); + fprintf(stdout,"gammaeos = %26.16e \n",gammaeos ); + fflush(stdout); +#endif + + /* done! */ + return; + +} + + +/**********************************************************************/ +/**************************************************************************** + vsq_calc(): + + -- evaluate v^2 (spatial, normalized velocity) from + W = \gamma^2 w + +****************************************************************************/ +static CCTK_REAL vsq_calc(CCTK_REAL W, struct LocGlob *lgp) +{ + CCTK_REAL Wsq,Xsq,Bsq_W; + + Wsq = W*W ; + Bsq_W = (lgp->Bsq + W); + Xsq = Bsq_W * Bsq_W; + + return( ( Wsq * lgp->Qtsq + lgp->QdotBsq * (Bsq_W + W)) / (Wsq*Xsq) ); +} + + +/******************************************************************** + + x1_of_x0(): + + -- calculates v^2 from W with some physical bounds checking; + -- asumes x0 is already physical + -- makes v^2 physical if not; + +*********************************************************************/ + +static CCTK_REAL x1_of_x0(CCTK_REAL x0, struct LocGlob *lgp ) +{ + CCTK_REAL x1,vsq; + CCTK_REAL dv = 1.e-15; + + vsq = fabs(vsq_calc(x0,lgp)) ; // guaranteed to be positive + + return( ( vsq > 1. ) ? (1.0 - dv) : vsq ); + +} + +/******************************************************************** + + validate_x(): + + -- makes sure that x[0,1] have physical values, based upon + their definitions: + +*********************************************************************/ + +static void validate_x(CCTK_REAL x[2], CCTK_REAL x0[2] ) +{ + + const CCTK_REAL dv = 1.e-15; + + /* Always take the absolute value of x[0] and check to see if it's too big: */ + x[0] = fabs(x[0]); + x[0] = (x[0] > W_TOO_BIG) ? x0[0] : x[0]; + + x[1] = (x[1] < 0.) ? 0. : x[1]; /* if it's too small */ + x[1] = (x[1] > 1.) ? (1. - dv) : x[1]; /* if it's too big */ + + return; + +} + +/************************************************************ + + twod_newton_raphson(): + + -- performs Newton-Rapshon method on an 2d system for polytropes. + + -- inspired in part by Num. Rec.'s routine newt(); + +*****************************************************************/ +static CCTK_INT twod_newton_raphson( CCTK_REAL x[], CCTK_REAL gammaeos, struct LocGlob *lgp, + void (*funcd) (CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], + CCTK_REAL [][2], CCTK_REAL *, + CCTK_REAL *, CCTK_REAL, struct LocGlob *) ) +{ + CCTK_REAL f, df, dx[2], x_old[2]; + CCTK_REAL resid[2], jac[2][2]; + CCTK_REAL errx, x_orig[2]; + CCTK_INT n_iter, id, jd, i_extra, doing_extra; + CCTK_REAL dW,dvsq,vsq_old,vsq,W,W_old; + const CCTK_REAL dv = (1.-1.e-15); + + CCTK_INT keep_iterating; + + + // Initialize various parameters and variables: + errx = 1. ; + df = f = 1.; + i_extra = doing_extra = 0; + x_old[0] = x_orig[0] = x[0] ; + x_old[1] = x_orig[1] = x[1] ; + + vsq_old = vsq = W = W_old = 0.; + n_iter = 0; + + /* Start the Newton-Raphson iterations : */ + keep_iterating = 1; + while( keep_iterating ) { + + (*funcd) (x, dx, resid, jac, &f, &df, gammaeos, lgp); /* returns with new dx, f, df */ + + + /* Save old values before calculating the new: */ + errx = 0.; + x_old[0] = x[0] ; + x_old[1] = x[1] ; + + /* Make the newton step: */ + x[0] += dx[0] ; + x[1] += dx[1] ; + + /****************************************/ + /* Calculate the convergence criterion */ + /****************************************/ + errx = (x[0]==0.) ? fabs(dx[0]) : fabs(dx[0]/x[0]); + + + /****************************************/ + /* Make sure that the new x[] is physical : */ + /****************************************/ + if( x[0] < 0. ) { x[0] = fabs(x[0]); } + else { + if(x[0] > W_TOO_BIG) { x[0] = x_old[0] ; } + } + + if( x[1] < 0. ) { x[1] = 0.; } + else { + if( x[1] > 1. ) { x[1] = dv; } + } + + /*****************************************************************************/ + /* If we've reached the tolerance level, then just do a few extra iterations */ + /* before stopping */ + /*****************************************************************************/ + + if( (fabs(errx) <= NEWT_TOL) && (doing_extra == 0) && (EXTRA_NEWT_ITER > 0) ) { + doing_extra = 1; + } + + if( doing_extra == 1 ) i_extra++ ; + + if( ((fabs(errx) <= NEWT_TOL)&&(doing_extra == 0)) + || (i_extra > EXTRA_NEWT_ITER) || (n_iter >= (MAX_NEWT_ITER-1)) ) { + keep_iterating = 0; + } + + n_iter++; + + } // END of while(keep_iterating) + + + /* Check for bad untrapped divergences : */ + if( (!finite(f)) || (!finite(df)) ) { + return(2); + } + + + if( fabs(errx) <= NEWT_TOL ){ + return(0); + } + else if( (fabs(errx) <= MIN_NEWT_TOL) && (fabs(errx) > NEWT_TOL) ){ + return(0); + } + else { + return(1); + } + + return(0); + +} + + +/************************************************************ + + threed_newton_raphson_omni(): + + -- performs Newton-Rapshon method on an 2d system for polytropes. + + -- inspired in part by Num. Rec.'s routine newt(); + +*****************************************************************/ +static CCTK_INT threed_newton_raphson_omni( CCTK_REAL x[], struct eosomnivars *eosvars, struct LocGlob *lgp, + void (*funcd) (CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], + CCTK_REAL [][3], CCTK_REAL *, + CCTK_REAL *, struct eosomnivars *, struct LocGlob *) ) +{ + CCTK_REAL f, df, dx[3], x_old[3]; + CCTK_REAL resid[3], jac[3][3]; + CCTK_REAL errx, x_orig[3]; + CCTK_INT n_iter, id, jd, i_extra, doing_extra; + CCTK_REAL dW,dvsq,du,vsq_old,vsq,W,W_old,u,u_old; + const CCTK_REAL dv = (1.-1.e-15); + + CCTK_INT keep_iterating; + + + // Initialize various parameters and variables: + errx = 1. ; + df = f = 1.; + i_extra = doing_extra = 0; + x_old[0] = x_orig[0] = x[0] ; + x_old[1] = x_orig[1] = x[1] ; + x_old[2] = x_orig[2] = x[2] ; + + vsq_old = vsq = W = W_old = u = u_old = 0.; + n_iter = 0; + + /* Start the Newton-Raphson iterations : */ + keep_iterating = 1; + while( keep_iterating ) { + + (*funcd) (x, dx, resid, jac, &f, &df, eosvars, lgp); /* returns with new dx, f, df */ + + /* Save old values before calculating the new: */ + errx = 0.; + x_old[0] = x[0] ; + x_old[1] = x[1] ; + x_old[2] = x[2] ; + + /* Make the newton step: */ + 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 */ + /****************************************/ + errx = (x[0]==0.) ? fabs(dx[0]) : fabs(dx[0]/x[0]); + + + /****************************************/ + /* Make sure that the new x[] is physical : */ + /****************************************/ + if( x[0] < 0. ) { x[0] = fabs(x[0]); } + else { + if(x[0] > W_TOO_BIG) { x[0] = x_old[0] ; } + } + + if( x[1] < 0. ) { x[1] = 0.; } + else { + if( x[1] > 1. ) { x[1] = dv; } + } + + if( x[2] < 0. ) { x[2] = 0.; } + + /*****************************************************************************/ + /* If we've reached the tolerance level, then just do a few extra iterations */ + /* before stopping */ + /*****************************************************************************/ + + if( (fabs(errx) <= NEWT_TOL) && (doing_extra == 0) && (EXTRA_NEWT_ITER > 0) ) { + doing_extra = 1; + } + + if( doing_extra == 1 ) i_extra++ ; + + if( ((fabs(errx) <= NEWT_TOL)&&(doing_extra == 0)) + || (i_extra > EXTRA_NEWT_ITER) || (n_iter >= (MAX_NEWT_ITER-1)) ) { + keep_iterating = 0; + } + + n_iter++; + + } // END of while(keep_iterating) + + + /* Check for bad untrapped divergences : */ + if( (!finite(f)) || (!finite(df)) ) { + return(2); + } + + + if( fabs(errx) <= eosvars->prec ){ + return(0); + } + else if( (fabs(errx) <= eosvars->prec) && (fabs(errx) > eosvars->prec) ){ + return(0); + } + //if( fabs(errx) <= NEWT_TOL ){ + // return(0); + //} + //else if( (fabs(errx) <= MIN_NEWT_TOL) && (fabs(errx) > NEWT_TOL) ){ + // return(0); + //} + else { + return(1); + } + + return(0); + +} + +/**********************************************************************/ +/********************************************************************************* + func_vsq(): + + -- calculates the residuals, and Newton step for general_newton_raphson(); + -- for this method, x=W,vsq here; + + Arguments: + x = current value of independent var's (on input & output); + dx = Newton-Raphson step (on output); + resid = residuals based on x (on output); + jac = Jacobian matrix based on x (on output); + f = resid.resid/2 (on output) + df = -2*f; (on output) + n = dimension of x[]; + *********************************************************************************/ + +static void func_vsq(CCTK_REAL x[], CCTK_REAL dx[], CCTK_REAL resid[], + CCTK_REAL jac[][2], CCTK_REAL *f, CCTK_REAL *df, CCTK_REAL gammaeos, struct LocGlob *lgp) +{ + + + CCTK_REAL W, vsq, Wsq, p_tmp, dPdvsq, dPdW; + CCTK_REAL res0, QB2Winv2,res1,j11,detJinv, Winv, QB2Winv3, j10, B2plusW, detJ, t36, mj01, j00; + + + W = x[0]; + vsq = x[1]; + + Wsq = W*W; + + p_tmp = eos_info(W, vsq, &dPdW, &dPdvsq, gammaeos, lgp); + + // These expressions were calculated using Mathematica, but made into efficient + // code using Maple. Since we know the analytic form of the equations, we can + // explicitly calculate the Newton-Raphson step: + + //j11 = dP/dv^2-B^2/2 + j11 = -lgp->half_Bsq+dPdvsq; + + B2plusW = lgp->Bsq+W; + + //mj01 is B2plusW squared = - (partial Eq. 4 / partial v^2) + mj01 = B2plusW*B2plusW; + + Winv = 1/W; + + QB2Winv2 = lgp->QdotBsq*Winv*Winv; + + //Eq. 4 - Residual 0: Qtsq - v^2(B^2+W)^2 -QdotBsq(B^2+2W)/W^2 + res0 = lgp->Qtsq-vsq*mj01+QB2Winv2*(lgp->Bsq+W+W); + + //Eq. 5 - Residual 1: -Qdotn - B^2/2(1+v^2)+1/2 QdotBsq/W^2 - W+p + res1 = -lgp->Qdotn-lgp->half_Bsq*(1.0+vsq)+0.5*QB2Winv2-W+p_tmp; + + QB2Winv3 = QB2Winv2*Winv; + + //j10 is -QB2Winv3 - 1 + dp/dW - (partial Eq. 5 / partial W) + j10 = -1.0+dPdW-QB2Winv3; + + //This is detJ: j10*mj01/B2W + -2 j11 * -j00/2 + detJ = B2plusW*(j10*B2plusW+(lgp->Bsq-2.0*dPdvsq)*(QB2Winv2+vsq*W)*Winv); + + detJinv = 1/detJ; + + // - (Jinv00 * res0 + Jinv01 * res 1)/detJ + dx[0] = -(j11*res0+mj01*res1)*detJinv; + + //j00 is -2v^2(B^2+W)-2QB2 (B2+W)/W^3 - (partial Eq. 4 / partial W) + j00 = -2*(vsq+QB2Winv3)*B2plusW; + + // (-Jinv10 * res0 -Jinv11 * res1) / DetJ + dx[1] = (j10*res0-j00*res1)*detJinv; + // detJ = B2plusW*detJ_gcf; + jac[0][0] = j00; + jac[0][1] = -mj01; + jac[1][0] = j10; + jac[1][1] = j11; + resid[0] = res0; + resid[1] = res1; + + *df = -resid[0]*resid[0] - resid[1]*resid[1]; + + *f = -0.5 * ( *df ); + +} + +/**********************************************************************/ +/********************************************************************************* + func_vsq_eosomni(): + + -- calculates the residuals, and Newton step for general_newton_raphson(); + -- for this method, x=W,vsq,u here; + + Arguments: + x = current value of independent var's (on input & output); + dx = Newton-Raphson step (on output); + resid = residuals based on x (on output); + jac = Jacobian matrix based on x (on output); + f = resid.resid/2 (on output) + df = -2*f; (on output) + n = dimension of x[]; + *********************************************************************************/ + +static void func_vsq_eosomni(CCTK_REAL x[], CCTK_REAL dx[], CCTK_REAL resid[], + CCTK_REAL jac[][3], CCTK_REAL *f, CCTK_REAL *df, + struct eosomnivars *eosvars, struct LocGlob *lgp) +{ + + CCTK_REAL W, vsq, u, p_tmp,epsilon,Wsq, dPdvsq, dPdW; + CCTK_REAL res0, LorInv, rho0, Winv, QB2Winv2, drho0_dv, res1, j11, detJ,detJinv; + CCTK_REAL QB2Winv3, j10,B2plusW,detJ_gcf,t36,B2plusW_sq, j00,j01,j12,j20,j21,j22; + CCTK_REAL dpress_dv,dpress_du,dpdrho,dpdeps,c00,c01,c02,c10,c11,c12,c20,c21,c22; + + + W = x[0]; + vsq = x[1]; + u = x[2]; + + Wsq = W*W; + + LorInv = sqrt(1.0-vsq); + rho0 = lgp->D * LorInv; + epsilon = u/rho0; + p_tmp = pressure_rho0_eps_eosomni(rho0,epsilon,&dpdrho,&dpdeps,eosvars); + + B2plusW = lgp->Bsq+W; + B2plusW_sq = B2plusW*B2plusW; + Winv = 1/W; + QB2Winv2 = lgp->QdotBsq*Winv*Winv; + QB2Winv3=QB2Winv2*Winv; + + //Eq. 4: Qtsq-v^2(B^2+W)^2-QdotBsq(B^2+2W)/W^2=0 <-No u or p dependence + + resid[0] = lgp->Qtsq - vsq*B2plusW_sq+QB2Winv2*(lgp->Bsq+W+W); + + j00 = -2*(vsq+QB2Winv3)*B2plusW; + j01 = -1.0*B2plusW_sq; + + + //Eq. 5: -Qdotn - B^2(1+vsq)/2 + QdotBsq/2/W^2 - vsq W - rho_0(vsq) - u = 0 + //rho0 = D * sqrt(1-vsq) + + resid[1] = -lgp->Qdotn - lgp->half_Bsq*(1.0+vsq) + 0.5*QB2Winv2 - vsq*W - rho0 - u; + + drho0_dv = -0.5*lgp->D / LorInv; + + j10 = -vsq-QB2Winv3; + j11 = -lgp->half_Bsq - W - drho0_dv; + + //Eq. 6: u+ p - W(1-vsq) + rho0 = 0 => p-W = -W vsq - rho0 - u + + resid[2] = u + p_tmp - W*(1.0-vsq) + rho0; + + //dp/dv = (dp/drho)u * drho/dv + //dp/drho_u = dp/drho_eps - eps/rho0 dpdeps + dpress_dv = drho0_dv*dpdrho+ + u/2.0/lgp->D*pow(1.0-vsq,-1.5)*dpdeps; + + //dp/du = 1/rho dp deps, since rho0 is function of v only + dpress_du=dpdeps/rho0; + + jac[0][0] = j00; + jac[0][1] = j01; + jac[0][2] = 0.0; + jac[1][0] = j10; + jac[1][1] = j11; + jac[1][2] = -1.0; + jac[2][0] = vsq-1.0; + jac[2][1] = dpress_dv + W + drho0_dv; + jac[2][2] = dpress_du + 1.0; + + c00 = j11*jac[2][2]+jac[2][1]; + c01 = -1*j01*jac[2][2]; + c02 = j01*jac[1][2]; + + c10 = -jac[2][0]-j10*jac[2][2]; + c11 = j00*jac[2][2]; + c12 = j00; + + c20 = j10*jac[2][1]-j11*jac[2][0]; + c21 = j01*jac[2][0]-j00*jac[2][1]; + c22 = j00*j11-j01*j10; + + detJ=j00*c00+j01*c10; + detJinv = 1/detJ; + + dx[0] = -(c00*resid[0]+c01*resid[1]+c02*resid[2])*detJinv; + dx[1] = -(c10*resid[0]+c11*resid[1]+c12*resid[2])*detJinv; + dx[2] = -(c20*resid[0]+c21*resid[1]+c22*resid[2])*detJinv; + + *df = -resid[0]*resid[0] - resid[1]*resid[1] - resid[2]*resid[2]; + + *f = -0.5 * ( *df ); + +} + + +/********************************************************************** + ********************************************************************** + + The following routines specify the equation of state. All routines + above here should be indpendent of EOS. If the user wishes + to use another equation of state, the below functions must be replaced + by equivalent routines based upon the new EOS. + + ********************************************************************** +**********************************************************************/ + +/**********************************************************************/ +/********************************************************************** + eos_info(): + + -- returns with all the EOS-related values needed; + **********************************************************************/ +static CCTK_REAL eos_info(CCTK_REAL W, CCTK_REAL vsq, CCTK_REAL *dpdw, CCTK_REAL *dpdvsq, CCTK_REAL gammaeos, struct LocGlob *lgp) +{ + register CCTK_REAL ftmp,gtmp; + + ftmp = 1. - vsq; + gtmp = sqrt(ftmp); + + CCTK_REAL gam_m1_o_gam = ((gammaeos-1.)/gammaeos); + + *dpdw = gam_m1_o_gam * ftmp ; + *dpdvsq = gam_m1_o_gam * ( 0.5 * lgp->D/gtmp - W ) ; + + return( gam_m1_o_gam * ( W * ftmp - lgp->D * gtmp ) ); // p + +} + +static CCTK_REAL pressure_rho0_eps_eosomni(CCTK_REAL rho,CCTK_REAL epsilon, CCTK_REAL* dpdrho, CCTK_REAL* dpdeps, struct eosomnivars *eosvars) +{ + + CCTK_REAL rhopt[1],epspt[1],press[1]; + rhopt[0]=rho; + epspt[0]=epsilon; + + EOS_Omni_press(eosvars->eoshandle,eosvars->eoskeytemp,eosvars->eosprec,1, + &rho,&epsilon,eosvars->eos_temp, + eosvars->eos_y_e,press,eosvars->eoskeyerr,eosvars->eosanyerr); + + EOS_Omni_DPressByDRho(eosvars->eoshandle,eosvars->eoskeytemp,eosvars->eosprec,1, + &rho,&epsilon,eosvars->eos_temp, + eosvars->eos_y_e,dpdrho,eosvars->eoskeyerr,eosvars->eosanyerr); + + EOS_Omni_DPressByDEps(eosvars->eoshandle,eosvars->eoskeytemp,eosvars->eosprec,1, + &rho,&epsilon,eosvars->eos_temp, + eosvars->eos_y_e,dpdeps,eosvars->eoskeyerr,eosvars->eosanyerr); + + return press[0]; + +} + + +/****************************************************************************** + END + ******************************************************************************/ + + +#undef DEBUG_CON2PRIMM diff --git a/src/GRHydro_Con2PrimM_ptee.c b/src/GRHydro_Con2PrimM_ptee.c index de58d4c..fcf639a 100644 --- a/src/GRHydro_Con2PrimM_ptee.c +++ b/src/GRHydro_Con2PrimM_ptee.c @@ -61,9 +61,9 @@ #define NEWT_TOL2 (1.0e-15) /* TOL of new 1D^*_{v^2} gnr2 method */ #define MIN_NEWT_TOL2 (1.0e-10) /* TOL of new 1D^*_{v^2} gnr2 method */ -#define W_TOO_BIG (1.e20) /* \gamma^2 (\rho_0 + u + p) is assumed +#define W_TOO_BIG (1.e20) /* \gamma^2 (\rho_0 + u + p) is assumed to always be smaller than this. This - is used to detect solver failures */ + is used to detect solver failures */ #define FAIL_VAL (1.e30) /* Generic value to which we set variables when a problem arises */ @@ -166,11 +166,11 @@ RETURN VALUE: of retval = (i*100 + j) where given tolerances; 2 -> Newton-Raphson procedure encountered a numerical divergence (occurrence of "nan" or "+/-inf" ; - + j = 0 -> success 1 -> failure: some sort of failure in Newton-Raphson; 2 -> failure: unphysical vsq = v^2 value at initial guess; - 3 -> failure: W<0 or W>W_TOO_BIG + 3 -> failure: W<0 or W>W_TOO_BIG 4 -> failure: v^2 > 1 ( used to be 5 -> failure: rho,uu <= 0 but now sets epsnegative to non-zero ) @@ -226,38 +226,38 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_ptee) ( #if(DEBUG_CON2PRIMM) fprintf(stdout," *dens = %26.16e \n", *dens_in ); - fprintf(stdout," *sx = %26.16e \n", *sx_in ); - fprintf(stdout," *sy = %26.16e \n", *sy_in ); - fprintf(stdout," *sz = %26.16e \n", *sz_in ); - fprintf(stdout," *tau = %26.16e \n", *tau_in ); - fprintf(stdout," *Bconsx = %26.16e \n", *Bconsx_in ); - fprintf(stdout," *Bconsy = %26.16e \n", *Bconsy_in ); + fprintf(stdout," *sx = %26.16e \n", *sx_in ); + fprintf(stdout," *sy = %26.16e \n", *sy_in ); + fprintf(stdout," *sz = %26.16e \n", *sz_in ); + fprintf(stdout," *tau = %26.16e \n", *tau_in ); + fprintf(stdout," *Bconsx = %26.16e \n", *Bconsx_in ); + fprintf(stdout," *Bconsy = %26.16e \n", *Bconsy_in ); fprintf(stdout," *Bconsz = %26.16e \n", *Bconsz_in ); fprintf(stdout," *entropycons = %26.16e \n", *entropycons_in ); - fprintf(stdout," *rho = %26.16e \n", *rho ); - fprintf(stdout," *velx = %26.16e \n", *velx ); - fprintf(stdout," *vely = %26.16e \n", *vely ); - fprintf(stdout," *velz = %26.16e \n", *velz ); + fprintf(stdout," *rho = %26.16e \n", *rho ); + fprintf(stdout," *velx = %26.16e \n", *velx ); + fprintf(stdout," *vely = %26.16e \n", *vely ); + fprintf(stdout," *velz = %26.16e \n", *velz ); fprintf(stdout," *epsilon = %26.16e \n", *epsilon ); fprintf(stdout," *pressure = %26.16e \n", *pressure ); - fprintf(stdout," *Bx = %26.16e \n", *Bx ); - fprintf(stdout," *By = %26.16e \n", *By ); - fprintf(stdout," *Bz = %26.16e \n", *Bz ); - fprintf(stdout," *bsq = %26.16e \n", *bsq ); + fprintf(stdout," *Bx = %26.16e \n", *Bx ); + fprintf(stdout," *By = %26.16e \n", *By ); + fprintf(stdout," *Bz = %26.16e \n", *Bz ); + fprintf(stdout," *bsq = %26.16e \n", *bsq ); fprintf(stdout," *w_lorentz = %26.16e \n", *w_lorentz ); - fprintf(stdout," *gxx = %26.16e \n", *gxx ); - fprintf(stdout," *gxy = %26.16e \n", *gxy ); - fprintf(stdout," *gxz = %26.16e \n", *gxz ); - fprintf(stdout," *gyy = %26.16e \n", *gyy ); - fprintf(stdout," *gyz = %26.16e \n", *gyz ); - fprintf(stdout," *gzz = %26.16e \n", *gzz ); - fprintf(stdout," *uxx = %26.16e \n", *uxx ); - fprintf(stdout," *uxy = %26.16e \n", *uxy ); - fprintf(stdout," *uxz = %26.16e \n", *uxz ); - fprintf(stdout," *uyy = %26.16e \n", *uyy ); - fprintf(stdout," *uyz = %26.16e \n", *uyz ); - fprintf(stdout," *uzz = %26.16e \n", *uzz ); - fprintf(stdout," *det = %26.16e \n", *det ); + fprintf(stdout," *gxx = %26.16e \n", *gxx ); + fprintf(stdout," *gxy = %26.16e \n", *gxy ); + fprintf(stdout," *gxz = %26.16e \n", *gxz ); + fprintf(stdout," *gyy = %26.16e \n", *gyy ); + fprintf(stdout," *gyz = %26.16e \n", *gyz ); + fprintf(stdout," *gzz = %26.16e \n", *gzz ); + fprintf(stdout," *uxx = %26.16e \n", *uxx ); + fprintf(stdout," *uxy = %26.16e \n", *uxy ); + fprintf(stdout," *uxz = %26.16e \n", *uxz ); + fprintf(stdout," *uyy = %26.16e \n", *uyy ); + fprintf(stdout," *uyz = %26.16e \n", *uyz ); + fprintf(stdout," *uzz = %26.16e \n", *uzz ); + fprintf(stdout," *det = %26.16e \n", *det ); fprintf(stdout," *epsnegative = %10d \n", *epsnegative ); fprintf(stdout," *retval = %26.16e \n", *retval ); fflush(stdout); @@ -323,7 +323,7 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_ptee) ( gammasq = 1. / (1. - vsq); gamma = sqrt(gammasq); - + // Always calculate rho from D and gamma so that using D in EOS remains consistent // i.e. you don't get positive values for dP/d(vsq) . rho0 = lg.D / gamma ; @@ -339,7 +339,7 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_ptee) ( // *retval = 1.0*twod_newton_raphson( x_2d, gammaeos, &lg, func_vsq ) ; *retval = general_newton_raphson( x_1d, gammaeos, &lg, func_rho ) ; rho0 = x_1d[0]; - + /* Problem with solver, so return denoting error before doing anything further */ if( ((*retval) != 0.) || (rho0 == FAIL_VAL) ) { *retval = *retval*100.+1.; diff --git a/src/GRHydro_InterfacesM.h b/src/GRHydro_InterfacesM.h index 758667f..6d90adf 100644 --- a/src/GRHydro_InterfacesM.h +++ b/src/GRHydro_InterfacesM.h @@ -2,8 +2,96 @@ module Con2PrimM_fortran_interfaces implicit none interface + subroutine GRHydro_Con2PrimM_ptold( handle, & + local_gam, dens, & + sx, sy, sz, & + tau, & + Bconsx, Bconsy, Bconsz, & + rho, & + velx, vely, velz,& + epsilon, pressure,& + Bx, By, Bz, & + bsq,& + w_lorentz, & + gxx, gxy, gxz, & + gyy, gyz, gzz, & + uxx, uxy, uxz,& + uyy, uyz, uzz,& + det,& + epsnegative, & + retval) + + implicit none + CCTK_INT handle + CCTK_REAL local_gam + CCTK_REAL dens + CCTK_REAL sx, sy, sz + CCTK_REAL tau + CCTK_REAL Bconsx, Bconsy, Bconsz + CCTK_REAL rho + CCTK_REAL velx, vely, velz + CCTK_REAL epsilon, pressure + CCTK_REAL Bx, By, Bz + CCTK_REAL bsq + CCTK_REAL w_lorentz + CCTK_REAL gxx, gxy, gxz + CCTK_REAL gyy, gyz, gzz + CCTK_REAL uxx, uxy, uxz + CCTK_REAL uyy, uyz, uzz + CCTK_REAL det + CCTK_INT epsnegative + CCTK_REAL retval + end subroutine GRHydro_Con2PrimM_ptold + + subroutine GRHydro_Con2PrimM_pt2( handle, keytemp, eos_prec, prec, & + local_gam, dens, & + sx, sy, sz, & + tau, & + Bconsx, Bconsy, Bconsz, & + xye, xtemp, & + rho, & + velx, vely, velz,& + epsilon, pressure,& + Bx, By, Bz, & + bsq,& + w_lorentz, & + gxx, gxy, gxz, & + gyy, gyz, gzz, & + uxx, uxy, uxz,& + uyy, uyz, uzz,& + det,& + epsnegative, & + retval) + + implicit none + CCTK_INT handle + CCTK_INT keytemp + CCTK_REAL eos_prec + 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 + CCTK_REAL Bx, By, Bz + CCTK_REAL bsq + CCTK_REAL w_lorentz + CCTK_REAL gxx, gxy, gxz + CCTK_REAL gyy, gyz, gzz + CCTK_REAL uxx, uxy, uxz + CCTK_REAL uyy, uyz, uzz + CCTK_REAL det + CCTK_INT epsnegative + CCTK_REAL retval + end subroutine GRHydro_Con2PrimM_pt2 + - subroutine GRHydro_Con2PrimM_pt( handle, keytemp, prec, & + subroutine GRHydro_Con2PrimM_pt( handle, GRHydro_reflevel, i, j, k, x, y, z, keytemp, eos_prec, prec, & local_gam, dens, & sx, sy, sz, & tau, & @@ -25,7 +113,11 @@ module Con2PrimM_fortran_interfaces implicit none CCTK_INT handle + CCTK_INT GRHydro_reflevel + CCTK_INT i, j, k + CCTK_REAL x, y, z CCTK_INT keytemp + CCTK_REAL eos_prec CCTK_REAL prec CCTK_REAL local_gam CCTK_REAL dens diff --git a/src/make.code.defn b/src/make.code.defn index 7bdd9f6..e7f495c 100644 --- a/src/make.code.defn +++ b/src/make.code.defn @@ -45,6 +45,8 @@ SRCS = Utils.F90 \ GRHydro_EoSChangeGamma.F90 \ GRHydro_Tmunu.F90 \ GRHydro_Con2PrimM.F90 \ + GRHydro_Con2PrimM_pt_EOSOmniold.c \ + GRHydro_Con2PrimM_pt_EOSOmni.c \ GRHydro_Con2PrimM_pt.c \ GRHydro_Con2PrimM_ptee.c \ GRHydro_Con2PrimM_polytype_pt.c \ |