Commits

commits@lists.einsteintoolkit.org

2 participants
6997 discussions

[svn:einsteintoolkit] incoming/TOVSolverHot/src/ (Rev. 30)
by cott＠tapir.caltech.edu 04 Dec '10

04 Dec '10

User: cott Date: 2010/12/04 08:05 PM Modified: /TOVSolverHot/src/ tov.c Log: * small modifications to initial guesses for rho -- this helps if one wants to make a hot NS File Changes: Directory: /TOVSolverHot/src/ ============================= File [modified]: tov.c Delta lines: +75 -25 =================================================================== --- TOVSolverHot/src/tov.c 2010-12-05 01:19:49 UTC (rev 29) +++ TOVSolverHot/src/tov.c 2010-12-05 02:05:49 UTC (rev 30) @@ -125,7 +125,7 @@ @history @endhistory @@*/ -void TOV_C_Source_RHS(CCTK_REAL r, CCTK_REAL K, CCTK_REAL Gamma, +void TOV_C_Source_RHS(CCTK_REAL r, CCTK_REAL K, CCTK_REAL Gamma, CCTK_REAL old_data[NUMVARS], CCTK_REAL source_data[NUMVARS]) { CCTK_REAL LOCAL_TINY, PI; @@ -170,7 +170,8 @@ } else { keytemp = 0; } - rho = TOV_Rho_Central[0]; + // guess for rho + rho = pow(press / 100.0, 1.0 / 2.0); EOS_Omni_RhoFromPressEpsTempEnt(eoskey, keytemp, TOVSolverHot_eos_rf_prec, @@ -396,7 +397,7 @@ TOV_C_fill(source_data, 6, 0.0); TOV_C_Source_RHS(TOV_r_1d[i], - TOV_K, TOV_Gamma, + TOV_K, TOV_Gamma, in_data, source_data); RKLOOP k1[rk] = TOV_dr[star] * source_data[rk]; @@ -835,7 +836,10 @@ CCTK_INT keyerr = 0; CCTK_INT anyerr = 0; CCTK_INT n = 1; - CCTK_REAL xeps,xtemp,xye,xent; + CCTK_REAL xeps,xtemp,xye,xent,xprs; + + + if(eoskey==4) { xtemp = nuc_eos_temperature; xent = nuc_eos_entropy; @@ -845,30 +849,76 @@ } else { keytemp = 1; } + if(press_point[star] > 0.0e0) { + // initial guess + rho_point[star] = pow(press_point[star] / 100.0, 1.0 / 2.0); + xprs = press_point[star]; + EOS_Omni_RhoFromPressEpsTempEnt(eoskey, + keytemp, + TOVSolverHot_eos_rf_prec, + n, + &(rho_point[star]), + &(eps_point[star]), + &xtemp, + &xent, + &xye, + &(press_point[star]), + &keyerr, + &anyerr); + + if(anyerr) { + CCTK_VInfo(CCTK_THORNSTRING, "EOS problem in TOV 4: %d", + keyerr); + CCTK_VInfo(CCTK_THORNSTRING, "%15.6E %15.6E %15.6E %15.6E %15.6E", + rho_point[star],eps_point[star],xtemp,xye,xprs); + CCTK_WARN(0,"aborting"); + } + } else { + rho_point[star] = nuc_eos_rho_min; + keytemp = 1; + EOS_Omni_press(eoskey, + keytemp, + TOVSolverHot_eos_rf_prec, + n, + &rho_point[star], + &eps_point[star], + &xtemp, + &xye, + &xprs, + &keyerr, + &anyerr); + if(anyerr) { + CCTK_VInfo(CCTK_THORNSTRING, "EOS problem in TOV 4a: %d", + keyerr); + CCTK_VInfo(CCTK_THORNSTRING, "%15.6E %15.6E %15.6E %15.6E %15.6E", + rho_point[star],eps_point[star],xtemp,xye,press_point[star]); + CCTK_WARN(0,"aborting"); + } + } } else { keytemp = 0; - } - // initial guess - rho_point[star] = 0.001e0; - EOS_Omni_RhoFromPressEpsTempEnt(eoskey, - keytemp, - TOVSolverHot_eos_rf_prec, - n, - &(rho_point[star]), - &(eps_point[star]), + // initial guess + rho_point[star] = pow(press_point[star] / 100.0, 1.0 / 2.0); + EOS_Omni_RhoFromPressEpsTempEnt(eoskey, + keytemp, + TOVSolverHot_eos_rf_prec, + n, + &(rho_point[star]), + &(eps_point[star]), &xtemp, - &xent, - &xye, - &(press_point[star]), - &keyerr, - &anyerr); - - if(anyerr) { - CCTK_VInfo(CCTK_THORNSTRING, "EOS problem in TOV 4: %d", - keyerr); - CCTK_VInfo(CCTK_THORNSTRING, "%15.6E %15.6E %15.6E %15.6E %15.6E", - rho_point[star],eps_point[star],xtemp,xye,press_point[star]); - CCTK_WARN(0,"aborting"); + &xent, + &xye, + &(press_point[star]), + &keyerr, + &anyerr); + + if(anyerr) { + CCTK_VInfo(CCTK_THORNSTRING, "EOS problem in TOV 4: %d", + keyerr); + CCTK_VInfo(CCTK_THORNSTRING, "%15.6E %15.6E %15.6E %15.6E %15.6E", + rho_point[star],eps_point[star],xtemp,xye,press_point[star]); + CCTK_WARN(0,"aborting"); + } } if(rho_point[star] > nuc_eos_rho_min) {

1 0

[svn:einsteintoolkit] EOS_Omni/trunk/src/nuc_eos/ (Rev. 29)
by cott＠tapir.caltech.edu 04 Dec '10

04 Dec '10

User: cott Date: 2010/12/04 07:29 PM Added: /trunk/src/nuc_eos/ findrho.F90 Log: * add missing file File Changes: Directory: /trunk/src/nuc_eos/ ============================== File [added]: findrho.F90 Delta lines: +112 -0 =================================================================== --- trunk/src/nuc_eos/findrho.F90 (rev 0) +++ trunk/src/nuc_eos/findrho.F90 2010-12-05 01:29:37 UTC (rev 29) @@ -0,0 +1,112 @@ +!-*-f90-*- +subroutine findrho_press(lr0,lt,y,lpressin,keyerrr,tol) + + use eosmodule + + implicit none + + !given initial guess of density + real*8 :: lr0 + !working density + real*8 :: lr, lr_lastguess + !given T and Ye + real*8 :: lt,y + !pressure corresponding to lr1,lt,y + real*8 :: lpressin + + !accuracy we need rho to, ~1 part in tol^{-1} + real*8 :: tol + !will be derivatives of pressure wrt rho,T,ye + real*8 :: d1,d2,d3 + !helpers along the way + real*8 :: lpress_of_guess,lpress_of_lastguess,first_lpress,ldr,lr_new + !bounds of density + real*8 :: lrmax,lrmin + !counters & flags + integer :: rl = 0 + integer :: itmax,i,keyerrr + + keyerrr=0 + + itmax=20 ! use at most 20 iterations, then bisection + + lr=lr0 + + lrmax=logrho(nrho) + lrmin=logrho(1) + + !Note: We are using Ewald's Lagrangian interpolator here! + + !first use initial guess to estimate derivatives + call findthis(lr,lt,y,lpress_of_guess,alltables(:,:,:,1),d1,d2,d3) + + first_lpress = lpress_of_guess + + !preconditioning 1: do we already have the right rho? + if (abs(lpressin-lpress_of_guess).lt.tol*abs(lpressin)) then + lr0 = lr + return + endif + + + !otherwise we must newton-raphson to get the real rho + do i=1,itmax + + !d1 is the derivative dlpress/dlogrho, use to guess hopefully better rho +! write(*,*) i,lr,d1,lpressin-lpress_of_guess + ldr= (lpressin-lpress_of_guess)/d1 + if (abs(ldr).gt.5.0d0) then + write(*,*) i,ldr,d1,lr,lr_new + keyerrr = 473 + write(*,*) "dpdrho very small" + return + endif + lr_new = lr+ldr + + !make sure we do not go out of bounds + lr_new = min(lr_new,lrmax) + lr_new = max(lr_new,lrmin) + + !keep old variables incase we want to make our own derivatives + lpress_of_lastguess = lpress_of_guess + lr_lastguess = lr + + !use to get next iteration of NR + lr = lr_new + call findthis(lr,lt,y,lpress_of_guess,alltables(:,:,:,1),d1,d2,d3) + if (abs(lpressin - lpress_of_guess).lt.tol*abs(lpressin)) then + !yes, we got it! + lr0 = lr + return + endif + + ! if we are closer than 10^-3 to the + ! root (lpressin-lpress_of_guess)=0, we are switching to + ! the secant method, since the table is rather coarse and the + ! derivatives may be garbage. + if(abs(lpressin-lpress_of_guess).lt.1.0d-3*abs(lpressin)) then + d1 = (lpress_of_guess-lpress_of_lastguess)/(lr-lr_lastguess) + endif + + enddo + + !we may fail in the NR, after itmax reached. Then we revert to the bisection method. + if(i.ge.itmax) then + keyerrr=667 + call bisection_rho(lr0,lt,y,lpressin,lr,alltables(:,:,:,1),keyerrr,3) + if(keyerrr.eq.667) then + ! total failure + call findthis(lr,lt,y,lpress_of_guess,alltables(:,:,:,1),d1,d2,d3) + write(*,*) "EOS: Did not converge in findrho_press!" + write(*,*) "iteration,logrho0,logtemp,ye,lpressin,lpress_first,rhoreturn,press_of_rhoreturn" + write(*,"(i4,1P10E19.10)") i,lr0,lt,y,lpressin,first_lpress,lr,lpress_of_guess + write(*,*) "Tried calling bisection... didn't help... :-/" + write(*,*) "Bisection error: ",keyerrr + endif + endif + + lr0 = lr + return + + +end subroutine findrho_press

1 0

[svn:einsteintoolkit] incoming/TOVSolverHot/src/ (Rev. 29)
by cott＠tapir.caltech.edu 04 Dec '10

04 Dec '10

User: cott Date: 2010/12/04 07:19 PM Modified: /TOVSolverHot/src/ tov.c Log: * small mapping fix File Changes: Directory: /TOVSolverHot/src/ ============================= File [modified]: tov.c Delta lines: +5 -0 =================================================================== --- TOVSolverHot/src/tov.c 2010-12-05 00:44:50 UTC (rev 28) +++ TOVSolverHot/src/tov.c 2010-12-05 01:19:49 UTC (rev 29) @@ -848,6 +848,8 @@ } else { keytemp = 0; } + // initial guess + rho_point[star] = 0.001e0; EOS_Omni_RhoFromPressEpsTempEnt(eoskey, keytemp, TOVSolverHot_eos_rf_prec, @@ -1058,6 +1060,7 @@ (r_to_star[star_i] * r_to_star[star_i] + 1.0e-30)) / my_psi4; rho[i3D] += rho_point[star_i]; + eps[i3D] += eps_point[star_i]; press[i3D] += press_point[star_i]; if(eoskey==4) { @@ -1077,6 +1080,7 @@ /* Reset atmosphere according to chosen star */ if(rho[i3D] <= TOV_Atmosphere[star_i]) rho[i3D] = TOV_Atmosphere[star_i]; + } if (TOV_Conformal_Flat_Three_Metric) @@ -1141,6 +1145,7 @@ Y_e_con[i3D] = Y_e[i3D] * dens[i3D]; } + } /* if used, recalculate the derivatives of the conformal factor */ if (*conformal_state > 1)

1 0

[svn:einsteintoolkit] EOS_Omni/trunk/src/nuc_eos/ (Rev. 28)
by cott＠tapir.caltech.edu 04 Dec '10

04 Dec '10

User: cott Date: 2010/12/04 07:19 PM Modified: /trunk/src/nuc_eos/ nuc_eos.F90 Log: * small fix File Changes: Directory: /trunk/src/nuc_eos/ ============================== File [modified]: nuc_eos.F90 Delta lines: +1 -1 =================================================================== --- trunk/src/nuc_eos/nuc_eos.F90 2010-12-05 00:46:52 UTC (rev 27) +++ trunk/src/nuc_eos/nuc_eos.F90 2010-12-05 01:19:33 UTC (rev 28) @@ -238,7 +238,7 @@ keyerr = 0 -if(keytemp.eq.0) then + if(keytemp.eq.0) then !need to find temperature based on xeps call findtemp(lr,lt,y,leps,keyerrt,rfeps) if(keyerrt.ne.0) then

1 0

[svn:einsteintoolkit] EOS_Omni/trunk/ (Rev. 27)
by cott＠tapir.caltech.edu 04 Dec '10

04 Dec '10

User: cott Date: 2010/12/04 06:46 PM Modified: /trunk/ interface.ccl /trunk/src/ EOS_Omni_SingleVarCalls.F90 /trunk/src/nuc_eos/ bisection.F90, make.code.defn, nuc_eos.F90 Log: * updates needed to make TOVSolverHot work File Changes: Directory: /trunk/src/nuc_eos/ ============================== File [modified]: bisection.F90 Delta lines: +83 -2 =================================================================== --- trunk/src/nuc_eos/bisection.F90 2010-11-26 19:45:37 UTC (rev 26) +++ trunk/src/nuc_eos/bisection.F90 2010-12-05 00:46:52 UTC (rev 27) @@ -7,7 +7,7 @@ real*8 lr,lt0,y,eps0,lt integer keyerrt - integer keybisect + integer keybisect !this doesn't do anything... ! 1 -> operate on energy ! 2 -> operate on entropy @@ -19,7 +19,7 @@ integer bcount,i,itmax,maxbcount - real*8 bivar(*) + real*8 bivar bcount = 0 maxbcount = 80 @@ -80,3 +80,84 @@ end subroutine bisection + +subroutine bisection_rho(lr0,lt,y,lpressin,lr,bivar,keyerrr,keybisect) + + use eosmodule + + implicit none + + real*8 lr,lr0,y,lpressin,lt + integer keyerrr + + integer keybisect !this doesn't do anything + ! 3 -> operate on pressure + + !temporary vars + real*8 lr1,lr2,lrmin,lrmax + real*8 f1,f2,fmid,dlr,lrmid + real*8 d1,d2,d3,tol + real*8 f1a,f2a + + integer bcount,i,itmax,maxbcount + + real*8 bivar + + bcount = 0 + maxbcount = 80 + + tol=1.d-9 ! need to find energy to less than 1 in 10^-9 + itmax=50 + + lrmax=logrho(nrho) + lrmin=logrho(1) + + lr = lr0 + lr1 = dlog10(min(10.0d0**lrmax,1.10d0*(10.0d0**lr0))) + lr2 = dlog10(max(10.0d0**lrmin,0.90d0*(10.0d0**lr0))) + + call findthis(lr1,lt,y,f1a,bivar,d1,d2,d3) + call findthis(lr2,lt,y,f2a,bivar,d1,d2,d3) + + f1=f1a-lpressin + f2=f2a-lpressin + + keyerrr=0 + + do while(f1*f2.ge.0.0d0) + bcount=bcount+1 + lr1=dlog10(min(10.0d0**lrmax,1.2d0*(10.0d0**lr1))) + lr2=dlog10(max(10.0d0**lrmin,0.8d0*(10.0d0**lr2))) + call findthis(lr1,lt,y,f1a,bivar,d1,d2,d3) + call findthis(lr2,lt,y,f2a,bivar,d1,d2,d3) + f1=f1a-lpressin + f2=f2a-lpressin + if(bcount.ge.maxbcount) then + keyerrr=667 + return + endif + enddo + + if(f1.lt.0.0d0) then + lr=lr1 + dlr=dlog10((10.0D0**lr2)-(10.0d0**lr1)) + else + lr=lr2 + dlr=dlog10((10.0D0**lr1)-(10.0d0**lr2)) + endif + + do i=1,itmax + dlr=dlog10((10.0d0**dlr)*0.5d0) + lrmid=dlog10(10.0d0**lr+10.0d0**dlr) + call findthis(lrmid,lt,y,f2a,bivar,d1,d2,d3) + fmid=f2a-lpressin + if(fmid.le.0.0d0) lr=lrmid + if(abs(1.0d0-f2a/lpressin).lt.tol) then + lr=lrmid + return + endif + enddo + + + +end subroutine bisection_rho File [modified]: make.code.defn Delta lines: +1 -1 =================================================================== --- trunk/src/nuc_eos/make.code.defn 2010-11-26 19:45:37 UTC (rev 26) +++ trunk/src/nuc_eos/make.code.defn 2010-12-05 00:46:52 UTC (rev 27) @@ -1,5 +1,5 @@ SRCS = eosmodule.F90 nuc_eos.F90 bisection.F90 \ - findtemp.F90 linterp_many.F90 readtable.F90 \ + findtemp.F90 findrho.F90 linterp_many.F90 readtable.F90 \ linterp.f SUBDIRS = File [modified]: nuc_eos.F90 Delta lines: +67 -36 =================================================================== --- trunk/src/nuc_eos/nuc_eos.F90 2010-11-26 19:45:37 UTC (rev 26) +++ trunk/src/nuc_eos/nuc_eos.F90 2010-12-05 00:46:52 UTC (rev 27) @@ -24,7 +24,7 @@ use eosmodule implicit none - real*8, intent(in) :: xrho,xye + real*8, intent(inout) :: xrho,xye real*8, intent(inout) :: xtemp,xenr,xent real*8, intent(out) :: xprs,xcs2,xdedt real*8, intent(out) :: xdpderho,xdpdrhoe,xxa,xxh,xxn,xxp @@ -35,10 +35,11 @@ ! local variables - real*8 :: lr,lt,y,xx,xeps,leps,xs + real*8 :: lr,lt,y,xx,xeps,leps,xs,xpressure real*8 :: d1,d2,d3 real*8 :: ff(nvars) integer :: keyerrt = 0 + integer :: keyerrr = 0 if(xrho.gt.eos_rhomax) then stop "nuc_eos: rho > rhomax" @@ -49,11 +50,13 @@ endif if(xye.gt.eos_yemax) then - stop "nuc_eos: ye > yemax" + keyerr = 101 + return endif if(xye.lt.eos_yemin) then - stop "nuc_eos: ye < yemin" + keyerr = 102 + return endif if(keytemp.eq.1) then @@ -87,19 +90,33 @@ xs = xent call findtemp_entropy(lr,lt,y,xs,keyerrt,rfeps) xtemp = 10.0d0**lt + + elseif(keytemp.eq.3) then + !need to find rho based on xprs + xpressure = log10(xprs) + call findrho_press(lr,lt,y,xpressure,keyerrr,rfeps) + if(keyerrr.ne.0) then + write(*,*) "Problem in findrho_press:", keyerr + keyerr = keyerrr + return + endif + xrho = 10.0d0**lr + endif ! have temperature, proceed: call findall(lr,lt,y,ff) - xprs = 10.0d0**ff(1) + !unless we want xprs to be constant (keytemp==3), reset xprs + if(.not.keytemp.eq.3) then + xprs = 10.0d0**ff(1) + endif - xeps = 10.0d0**ff(2) - energy_shift - if(keytemp.eq.1.or.keytemp.eq.2) then - xenr = xeps + if(.not.keytemp.eq.0) then + xenr = 10.0d0**ff(2) - energy_shift endif - if(keytemp.eq.1.or.keytemp.eq.0) then + if(.not.keytemp.eq.2) then xent = ff(3) endif @@ -167,7 +184,7 @@ use eosmodule implicit none - real*8, intent(in) :: xrho,xye + real*8, intent(inout) :: xrho,xye real*8, intent(inout) :: xtemp,xenr,xent real*8, intent(in) :: rfeps real*8, intent(out) :: xprs,xmunu,xcs2,xdedt @@ -176,9 +193,10 @@ integer, intent(out) :: keyerr ! local variables - real*8 :: lr,lt,y,xx,xeps,leps,xs + real*8 :: lr,lt,y,xx,xeps,leps,xs,xpressure real*8 :: d1,d2,d3,ff(8) integer :: keyerrt = 0 + integer :: keyerrr = 0 if(xrho.gt.eos_rhomax) then stop "nuc_eos: rho > rhomax" @@ -191,13 +209,13 @@ endif if(xye.gt.eos_yemax) then - stop "nuc_eos: ye > yemax" + keyerr = 101 + return endif if(xye.lt.eos_yemin) then keyerr = 102 - write(6,*) "ye: ", xye - stop "nuc_eos ye < yemin" + return endif if(keytemp.eq.1) then @@ -220,7 +238,7 @@ keyerr = 0 - if(keytemp.eq.0) then +if(keytemp.eq.0) then !need to find temperature based on xeps call findtemp(lr,lt,y,leps,keyerrt,rfeps) if(keyerrt.ne.0) then @@ -239,18 +257,33 @@ endif xtemp = 10.0d0**lt + elseif(keytemp.eq.3) then + !need to find rho based on xprs + xpressure = log10(xprs) + call findrho_press(lr,lt,y,xpressure,keyerrr,rfeps) + if (keyerrr.ne.0) then + keyerr = keyerrr + write(*,*) "Problem in findrho_press:", keyerr + return + endif + xrho = 10.0d0**lr endif - ! have temperature, proceed: + ! have rho,temp,ye; proceed: call findall_short(lr,lt,y,ff) - xprs = 10.0d0**ff(1) - xeps = 10.0d0**ff(2) - energy_shift - if((keytemp.eq.1).or.(keytemp.eq.2)) then - xenr = xeps + !unless we want xprs to be constant (keytemp==3), reset xprs + if(.not.keytemp.eq.3) then + xprs = 10.0d0**ff(1) endif - if(keytemp.eq.1.or.keytemp.eq.0) then + !unless we want xenr to be constant (keytemp==0), reset xenr + if(.not.keytemp.eq.0) then + xenr = 10.0d0**ff(2) - energy_shift + endif + + !unless we want xent to be constant (keytemp==2), reset xent + if(.not.keytemp.eq.2) then xent = ff(3) endif @@ -296,13 +329,13 @@ endif if(xye.gt.eos_yemax) then - stop "nuc_eos: ye > yemax" + keyerr = 101 + return endif if(xye.lt.eos_yemin) then keyerr = 102 - write(6,*) "ye: ", xye - stop "nuc_eos ye < yemin" + return endif if(keytemp.eq.1) then @@ -333,17 +366,16 @@ endif xtemp = 10.0d0**lt - elseif(keytemp.eq.2) then - stop "eos_nuc_press does not support keytemp.eq.2" + elseif(keytemp.gt.1) then + stop "eos_nuc_press does not support keytemp other than 0 and 1" endif ! have temperature, proceed: call findall_press_eps(lr,lt,y,ff) xprs = 10.0d0**ff(1) - xeps = 10.0d0**ff(2) - energy_shift - if((keytemp.eq.1).or.(keytemp.eq.2)) then - xenr = xeps + if(keytemp.eq.1) then + xenr = 10.0d0**ff(2) - energy_shift endif end subroutine nuc_eos_press_eps @@ -378,13 +410,13 @@ endif if(xye.gt.eos_yemax) then - stop "nuc_eos: ye > yemax" + keyerr = 101 + return endif if(xye.lt.eos_yemin) then keyerr = 102 - write(6,*) "ye: ", xye - stop "nuc_eos ye < yemin" + return endif if(keytemp.eq.1) then @@ -415,8 +447,8 @@ endif xtemp = 10.0d0**lt - elseif(keytemp.eq.2) then - stop "eos_nuc_press does not support keytemp.eq.2" + elseif(keytemp.gt.1) then + stop "eos_nuc_press does not support keytemp > 1" endif ! have temperature, proceed: @@ -424,8 +456,7 @@ xdpdrhoe = 10.0d0**ff(1) xdpderho = 10.0d0**ff(2) - - if((keytemp.eq.1).or.(keytemp.eq.2)) then + if(keytemp.eq.1) then call findthis(lr,lt,y,xeps,alltables(:,:,:,2),d1,d2,d3) xeps = 10.0d0**xeps - energy_shift xenr = xeps Directory: /trunk/src/ ====================== File [modified]: EOS_Omni_SingleVarCalls.F90 Delta lines: +53 -17 =================================================================== --- trunk/src/EOS_Omni_SingleVarCalls.F90 2010-11-26 19:45:37 UTC (rev 26) +++ trunk/src/EOS_Omni_SingleVarCalls.F90 2010-12-05 00:46:52 UTC (rev 27) @@ -509,38 +509,43 @@ end subroutine EOS_Omni_EOS_eps_from_press -subroutine EOS_Omni_EOS_RestMassDensityFromEpsPress(eoskey,keytemp,rf_precision,& - npoints,rho,eps,temp,ye,press,keyerr,anyerr) +subroutine EOS_Omni_EOS_RhoFromPressEpsTempEnt(eoskey,& + ikeytemp,rf_precision,& + npoints,rho,eps,temp,entropy,ye,press,keyerr,anyerr) use EOS_Omni_Module implicit none DECLARE_CCTK_PARAMETERS - CCTK_INT, intent(in) :: eoskey,keytemp,npoints + CCTK_INT, intent(in) :: eoskey,ikeytemp,npoints CCTK_INT, intent(out) :: keyerr(npoints) CCTK_INT, intent(out) :: anyerr CCTK_REAL, intent(in) :: rf_precision - CCTK_REAL, intent(in) :: ye(npoints),press(npoints),eps(npoints) - CCTK_REAL, intent(inout) :: rho(npoints),temp(npoints) + CCTK_REAL, intent(in) :: ye(npoints),press(npoints) + CCTK_REAL, intent(inout) :: rho(npoints),temp(npoints),eps(npoints) + CCTK_REAL, intent(inout) :: entropy(npoints) ! local vars integer :: i character(256) :: warnstring + integer :: keytemp + ! temporary vars for nuc_eos + real*8 :: xrho,xye,xtemp,xenr,xent + real*8 :: xprs,xmunu,xcs2 + real*8 :: xdedt,xdpderho,xdpdrhoe + + keytemp = ikeytemp - if(keytemp.eq.1) then - anyerr = 1 - keyerr(:) = -1 - else - anyerr = 0 - keyerr(:) = 0 - endif + anyerr = 0 + keyerr(:) = 0 select case (eoskey) case (1) ! polytropic EOS do i=1,npoints - rho(i) = press(i) / ((poly_gamma - 1.0d0)*eps(i)) + rho(i) = (press(i) / poly_k)**(1.0/poly_gamma) + eps(i) = press(i) / (poly_gamma - 1.0d0) / rho(i) enddo case (2) ! gamma-law EOS @@ -549,15 +554,46 @@ enddo case (3) ! hybrid EOS - write(warnstring,*) "EOS_Omni_RestMassDensityFromEpsPress not supported for hybrid EOS" + write(warnstring,*) "EOS_Omni_RestMassDensityFromPressEpsTemp not supported for hybrid EOS" call CCTK_WARN(0,warnstring) case (4) ! nuc EOS - write(warnstring,*) "EOS_Omni_RestMassDensityFromEpsPress call not supported for nuc_eos yet" - call CCTK_WARN(0,warnstring) + if(ikeytemp.eq.2) then + call CCTK_WARN(0,"This function does not work yet when coming in with entropy") + else if(ikeytemp.eq.1) then + keytemp = 3 + else + call CCTK_WARN(0,"This function does not work yet when coming in with this keytemp") + endif + + do i=1,npoints + ! initial guess + xprs = press(i) * inv_press_gf + xrho = rho(i) * inv_rho_gf + xtemp = temp(i) + xent = entropy(i) + xye = ye(i) + xenr = eps(i) * inv_eps_gf + call nuc_eos_short(xrho,xtemp,xye,xenr,xprs,& + xent,xcs2,xdedt,xdpderho,xdpdrhoe,xmunu,& + keytemp,keyerr(i),rf_precision) + + if(keyerr(i).ne.0) then + anyerr = 1 + endif + + if(keytemp.eq.3) then + eps(i) = xenr * eps_gf + rho(i) = xrho * rho_gf + entropy(i) = xent + else + call CCTK_WARN(0, "This keytemp is not implemented yet") + endif + + enddo case DEFAULT write(warnstring,*) "eoskey ",eoskey," not implemented!" call CCTK_WARN(0,warnstring) end select -end subroutine EOS_Omni_EOS_RestMassDensityFromEpsPress +end subroutine EOS_Omni_EOS_RhoFromPressEpsTempEnt Directory: /trunk/ ================== File [modified]: interface.ccl Delta lines: +4 -3 =================================================================== --- trunk/interface.ccl 2010-11-26 19:45:37 UTC (rev 26) +++ trunk/interface.ccl 2010-12-05 00:46:52 UTC (rev 27) @@ -100,19 +100,20 @@ PROVIDES FUNCTION EOS_Omni_EpsFromPress WITH EOS_Omni_EOS_eps_from_press LANGUAGE Fortran -void FUNCTION EOS_Omni_RestMassDensityFromEpsPress(CCTK_INT IN eoskey, \ +void FUNCTION EOS_Omni_RhoFromPressEpsTempEnt(CCTK_INT IN eoskey, \ CCTK_INT IN havetemp, \ CCTK_REAL IN rf_precision, \ CCTK_INT IN npoints, \ CCTK_REAL OUT ARRAY rho, \ - CCTK_REAL IN ARRAY eps, \ + CCTK_REAL INOUT ARRAY eps, \ CCTK_REAL INOUT ARRAY temp, \ + CCTK_REAL INOUT ARRAY ent, \ CCTK_REAL IN ARRAY ye, \ CCTK_REAL IN ARRAY press, \ CCTK_INT OUT ARRAY keyerr, \ CCTK_INT OUT anyerr) -PROVIDES FUNCTION EOS_Omni_RestMassDensityFromEpsPress WITH EOS_Omni_EOS_RestMassDensityFromEpsPress LANGUAGE Fortran +PROVIDES FUNCTION EOS_Omni_RhoFromPressEpsTempEnt WITH EOS_Omni_EOS_RhoFromPressEpsTempEnt LANGUAGE Fortran ################################################################################ # short and long and specific calls

1 0

[svn:einsteintoolkit] incoming/TOVSolverHot/src/ (Rev. 28)
by cott＠tapir.caltech.edu 04 Dec '10

04 Dec '10

User: cott Date: 2010/12/04 06:44 PM Modified: /TOVSolverHot/src/ tov.c Log: * fix stupid bug introduced by myself in mapping 1D->3D File Changes: Directory: /TOVSolverHot/src/ ============================= File [modified]: tov.c Delta lines: +12 -7 =================================================================== --- TOVSolverHot/src/tov.c 2010-12-05 00:09:30 UTC (rev 27) +++ TOVSolverHot/src/tov.c 2010-12-05 00:44:50 UTC (rev 28) @@ -869,10 +869,15 @@ CCTK_WARN(0,"aborting"); } - temperature_point[star] = xtemp; - entropy_point[star] = xent; - ye_point[star] = xye; - + if(rho_point[star] > nuc_eos_rho_min) { + temperature_point[star] = xtemp; + entropy_point[star] = xent; + ye_point[star] = xye; + } else { + temperature_point[star] = nuc_eos_temperature; + entropy_point[star] = nuc_eos_entropy; + ye_point[star] = nuc_eos_ye; + } } else if(!use_EOS_Omni) { @@ -1056,9 +1061,9 @@ eps[i3D] += eps_point[star_i]; press[i3D] += press_point[star_i]; if(eoskey==4) { - temperature[i3D] += temperature_point[star]; - entropy[i3D] += entropy_point[star]; - Y_e[i3D] += ye_point[star]; + temperature[i3D] += temperature_point[star_i]; + entropy[i3D] += entropy_point[star_i]; + Y_e[i3D] += ye_point[star_i]; } /* we still have to know if we are inside one star - and which */

1 0

[svn:einsteintoolkit] incoming/TOVSolverHot/ (Rev. 27)
by cott＠tapir.caltech.edu 04 Dec '10

04 Dec '10

User: cott Date: 2010/12/04 06:09 PM Modified: /TOVSolverHot/ param.ccl /TOVSolverHot/src/ tov.c Log: * this TOV solver can not make microphysical NSs File Changes: Directory: /TOVSolverHot/ ========================= File [modified]: param.ccl Delta lines: +15 -0 =================================================================== --- TOVSolverHot/param.ccl 2010-12-04 22:52:47 UTC (rev 26) +++ TOVSolverHot/param.ccl 2010-12-05 00:09:30 UTC (rev 27) @@ -155,11 +155,23 @@ 2 :: "entropy" } 1 +CCTK_REAL nuc_eos_rho_min "Minimum density to use with nuc_eos [Cactus]" +{ + 0.0:* :: "Better be positive or zero" +} 1.0e-13 + CCTK_REAL nuc_eos_temperature "Temperature to use when using nuc_eos [MeV]" { 0.1:* :: "Below 0.1 MeV tables get bad" } 1.0e0 + +CCTK_REAL nuc_eos_ye "Y_e to use when using nuc_eos" +{ + 0.1:* :: "Below 0.1 things get bad" +} 0.1e0 + + CCTK_REAL nuc_eos_entropy "Entropy to use when using nuc_eos [k_B/baryon]" { 0.3:* :: "Below 0.3 k_B/baryon things get dicey" @@ -179,6 +191,9 @@ "tov" :: "TOV star initial hydrobase variables" } +USES KEYWORD Y_e_evolution_method +USES KEYWORD temperature_evolution_method + shares:admbase EXTENDS KEYWORD initial_data Directory: /TOVSolverHot/src/ ============================= File [modified]: tov.c Delta lines: +164 -29 =================================================================== --- TOVSolverHot/src/tov.c 2010-12-04 22:52:47 UTC (rev 26) +++ TOVSolverHot/src/tov.c 2010-12-05 00:09:30 UTC (rev 27) @@ -159,10 +159,18 @@ CCTK_REAL xye; CCTK_REAL xent; if(eoskey==4) { - CCTK_WARN(0,"fix me"); + xtemp = nuc_eos_temperature; + xent = nuc_eos_entropy; + xye = nuc_eos_ye; + if(nuc_eos_keytemp==2) { + CCTK_WARN(0,"nuc_eos_keytemp=2 not yet handled"); + } else { + keytemp = 1; + } } else { keytemp = 0; } + rho = TOV_Rho_Central[0]; EOS_Omni_RhoFromPressEpsTempEnt(eoskey, keytemp, TOVSolverHot_eos_rf_prec, @@ -176,6 +184,15 @@ &keyerr, &anyerr); + + if(anyerr) { + CCTK_VInfo(CCTK_THORNSTRING, "EOS problem in TOV 1: %d %d", + keyerr,anyerr); + CCTK_VInfo(CCTK_THORNSTRING, "%15.6E %15.6E %15.6E %15.6E %15.6E", + rho,eps,xtemp,xye,press); + CCTK_WARN(0,"aborting"); + } + } mu = rho * (1.0 + eps); @@ -228,6 +245,7 @@ k1[NUMVARS], k2[NUMVARS], k3[NUMVARS], k4[NUMVARS]; CCTK_REAL Surface_Mass, factor, local_rho; CCTK_REAL rho_central; + CCTK_INT have_temp_ye; LOCAL_TINY = 1.0e-20; @@ -244,6 +262,15 @@ assert(TOV_mbary_1d!=0); assert(TOV_rprop_1d!=0); + if(use_EOS_Omni && eoskey == 4) { + if(CCTK_EQUALS(Y_e_evolution_method,"none") || + CCTK_EQUALS(temperature_evolution_method,"none")) { + CCTK_WARN(0,"temperature and/or Y_e evolution not active"); + } else { + have_temp_ye = 1; + } + } + /* do it for all stars */ for (star=0; star < TOV_Num_TOVs; star++) { @@ -305,7 +332,14 @@ CCTK_INT n = 1; CCTK_REAL xeps,xtemp,xye,xent; if(eoskey==4) { - CCTK_WARN(0,"fix me"); + xtemp = nuc_eos_temperature; + xent = nuc_eos_entropy; + xye = nuc_eos_ye; + if(nuc_eos_keytemp==2) { + CCTK_WARN(0,"nuc_eos_keytemp=2 not yet handled"); + } else { + keytemp = 1; + } } else { keytemp = 0; } @@ -320,6 +354,14 @@ &(TOV_press_1d[star_i]), &keyerr, &anyerr); + if(anyerr) { + CCTK_VInfo(CCTK_THORNSTRING, "EOS problem in TOV 2: %d", + keyerr); + CCTK_VInfo(CCTK_THORNSTRING, "%15.6E %15.6E %15.6E %15.6E %15.6E", + rho_central,eps,xtemp,xye,TOV_press_1d[star_i]); + CCTK_WARN(0,"aborting"); + } + } @@ -398,10 +440,18 @@ CCTK_INT n = 1; CCTK_REAL xeps,xtemp,xye,xent; if(eoskey==4) { - CCTK_WARN(0,"fix me"); + xye = nuc_eos_ye; + xtemp = nuc_eos_temperature; + xent = nuc_eos_entropy; + if(nuc_eos_keytemp==2) { + CCTK_WARN(0,"nuc_eos_keytemp=2 not yet handled"); + } else { + keytemp = 1; + } } else { keytemp = 0; } + local_rho = TOV_Rho_Central[0]; EOS_Omni_RhoFromPressEpsTempEnt(eoskey, keytemp, TOVSolverHot_eos_rf_prec, @@ -415,16 +465,36 @@ &keyerr, &anyerr); + if(anyerr) { + CCTK_VInfo(CCTK_THORNSTRING, "EOS problem in TOV 3: %d", + keyerr); + CCTK_VInfo(CCTK_THORNSTRING, "%15.6E %15.6E %15.6E %15.6E %15.6E", + local_rho,xeps,xtemp,xye,TOV_press_1d[i+1]); + CCTK_WARN(0,"aborting"); + } + // CCTK_VInfo(CCTK_THORNSTRING, "%6d %15.6E %15.6E", + // i,local_rho,TOV_rbar_1d[i]); + } /* scan for the surface */ - if ( (local_rho <= 0.0) || - (TOV_press_1d[i+1] <= 0.0) ) - { - TOV_Surface[star] = TOV_r_1d[i]; - TOV_R_Surface[star] = TOV_rbar_1d[i]; - TOV_RProp_Surface[star] = TOV_rprop_1d[i]; - TOV_Surface_Index = i; + if(!(use_EOS_Omni && eoskey==4)) { + if ( (local_rho <= 0.0) || + (TOV_press_1d[i+1] <= 0.0) ) + { + TOV_Surface[star] = TOV_r_1d[i]; + TOV_R_Surface[star] = TOV_rbar_1d[i]; + TOV_RProp_Surface[star] = TOV_rprop_1d[i]; + TOV_Surface_Index = i; + } + } else { + if ( (local_rho <= nuc_eos_rho_min) ) + { + TOV_Surface[star] = TOV_r_1d[i]; + TOV_R_Surface[star] = TOV_rbar_1d[i]; + TOV_RProp_Surface[star] = TOV_rprop_1d[i]; + TOV_Surface_Index = i; + } } } if (TOV_Surface[star] < 0.0) @@ -476,16 +546,31 @@ pow(CONSTANT_c_cgi,4.0), TOV_Gamma); } else { - CCTK_VInfo("", "TOV radius mass bary_mass cent.rho rho(cgs)"); - for (i=0; i<TOV_Num_TOVs; i++) - CCTK_VInfo(""," %d %8g %8g %8g %8.3g %8g", - i+1, TOV_R_Surface[i], - TOV_m_1d[(i+1)*TOV_Num_Radial-1], - TOV_mbary_1d[(i+1)*TOV_Num_Radial-1], - TOV_Rho_Central[i], - TOV_Rho_Central[i]/pow(CONSTANT_G_cgi,3.0)/ - pow(CONSTANT_Msolar_cgi,2.0)* + if(eoskey != 4 ) { + CCTK_VInfo("", "TOV radius mass bary_mass cent.rho rho(cgs)"); + for (i=0; i<TOV_Num_TOVs; i++) + CCTK_VInfo(""," %d %8g %8g %8g %8.3g %8g", + i+1, TOV_R_Surface[i], + TOV_m_1d[(i+1)*TOV_Num_Radial-1], + TOV_mbary_1d[(i+1)*TOV_Num_Radial-1], + TOV_Rho_Central[i], + TOV_Rho_Central[i]/pow(CONSTANT_G_cgi,3.0)/ + pow(CONSTANT_Msolar_cgi,2.0)* pow(CONSTANT_c_cgi,6.0)); + } else { + CCTK_VInfo(CCTK_THORNSTRING, "Made a hot polytrope: T = %5.2f MeV", + nuc_eos_temperature); + CCTK_VInfo(CCTK_THORNSTRING, "TOV radius mass bary_mass cent.rho rho(cgs)"); + for (i=0; i<TOV_Num_TOVs; i++) + CCTK_VInfo(CCTK_THORNSTRING," %d %8g %8g %8g %8.3g %8g", + i+1, TOV_R_Surface[i], + TOV_m_1d[(i+1)*TOV_Num_Radial-1], + TOV_mbary_1d[(i+1)*TOV_Num_Radial-1], + TOV_Rho_Central[i], + TOV_Rho_Central[i]/pow(CONSTANT_G_cgi,3.0)/ + pow(CONSTANT_Msolar_cgi,2.0)* + pow(CONSTANT_c_cgi,6.0)); + } } } @@ -593,7 +678,8 @@ DECLARE_CCTK_PARAMETERS CCTK_REAL *press_point, *rho_point, *eps_point, - *mu_point, *phi_point, *r_point; + *mu_point, *phi_point, *r_point, *temperature_point, + *ye_point,*entropy_point; CCTK_INT LSH_MAX_I; @@ -625,14 +711,18 @@ assert(TOV_m_1d!=0); /* allocate local arrays */ - r_to_star = (CCTK_REAL *) calloc (TOV_Num_TOVs, sizeof(CCTK_REAL)); - press_point = (CCTK_REAL *) calloc (TOV_Num_TOVs, sizeof(CCTK_REAL)); - rho_point = (CCTK_REAL *) calloc (TOV_Num_TOVs, sizeof(CCTK_REAL)); - eps_point = (CCTK_REAL *) calloc (TOV_Num_TOVs, sizeof(CCTK_REAL)); - mu_point = (CCTK_REAL *) calloc (TOV_Num_TOVs, sizeof(CCTK_REAL)); - phi_point = (CCTK_REAL *) calloc (TOV_Num_TOVs, sizeof(CCTK_REAL)); - r_point = (CCTK_REAL *) calloc (TOV_Num_TOVs, sizeof(CCTK_REAL)); + r_to_star = (CCTK_REAL *) calloc (TOV_Num_TOVs, sizeof(CCTK_REAL)); + press_point = (CCTK_REAL *) calloc (TOV_Num_TOVs, sizeof(CCTK_REAL)); + rho_point = (CCTK_REAL *) calloc (TOV_Num_TOVs, sizeof(CCTK_REAL)); + eps_point = (CCTK_REAL *) calloc (TOV_Num_TOVs, sizeof(CCTK_REAL)); + temperature_point = (CCTK_REAL *) calloc (TOV_Num_TOVs, sizeof(CCTK_REAL)); + ye_point = (CCTK_REAL *) calloc (TOV_Num_TOVs, sizeof(CCTK_REAL)); + entropy_point = (CCTK_REAL *) calloc (TOV_Num_TOVs, sizeof(CCTK_REAL)); + mu_point = (CCTK_REAL *) calloc (TOV_Num_TOVs, sizeof(CCTK_REAL)); + phi_point = (CCTK_REAL *) calloc (TOV_Num_TOVs, sizeof(CCTK_REAL)); + r_point = (CCTK_REAL *) calloc (TOV_Num_TOVs, sizeof(CCTK_REAL)); + /* clear initial data */ if (TOV_Clear_Initial_Data > 0 && !(TOV_Use_Old_Initial_Data)) { @@ -692,6 +782,12 @@ TOV_C_fill(velx, LSH_MAX_I+1, 0.0); TOV_C_fill(vely, LSH_MAX_I+1, 0.0); TOV_C_fill(velz, LSH_MAX_I+1, 0.0); + if(eoskey == 4) { + TOV_C_fill(temperature, LSH_MAX_I+1, 0.0); + TOV_C_fill(entropy, LSH_MAX_I+1, 0.0); + TOV_C_fill(Y_e, LSH_MAX_I+1, 0.0); + } + } /* use the fast interpolation? only useful for testing this */ if (TOV_Fast_Interpolation == 0) @@ -741,7 +837,14 @@ CCTK_INT n = 1; CCTK_REAL xeps,xtemp,xye,xent; if(eoskey==4) { - CCTK_WARN(0,"fix me"); + xtemp = nuc_eos_temperature; + xent = nuc_eos_entropy; + xye = nuc_eos_ye; + if(nuc_eos_keytemp==2) { + CCTK_WARN(0,"nuc_eos_keytemp=2 not yet handled"); + } else { + keytemp = 1; + } } else { keytemp = 0; } @@ -758,6 +861,18 @@ &keyerr, &anyerr); + if(anyerr) { + CCTK_VInfo(CCTK_THORNSTRING, "EOS problem in TOV 4: %d", + keyerr); + CCTK_VInfo(CCTK_THORNSTRING, "%15.6E %15.6E %15.6E %15.6E %15.6E", + rho_point[star],eps_point[star],xtemp,xye,press_point[star]); + CCTK_WARN(0,"aborting"); + } + + temperature_point[star] = xtemp; + entropy_point[star] = xent; + ye_point[star] = xye; + } else if(!use_EOS_Omni) { @@ -899,7 +1014,12 @@ rho[i3D] = TOV_Atmosphere[star]; eps[i3D] = eps_point[star]; press[i3D] = press_point[star]; - + + if(eoskey==4) { + temperature[i3D] = temperature_point[star]; + entropy[i3D] = entropy_point[star]; + Y_e[i3D] = ye_point[star]; + } } else if (CCTK_EQUALS(TOV_Combine_Method, "average")) { @@ -935,6 +1055,12 @@ rho[i3D] += rho_point[star_i]; eps[i3D] += eps_point[star_i]; press[i3D] += press_point[star_i]; + if(eoskey==4) { + temperature[i3D] += temperature_point[star]; + entropy[i3D] += entropy_point[star]; + Y_e[i3D] += ye_point[star]; + } + /* we still have to know if we are inside one star - and which */ if ( (rho_point[star_i] > max_rho) && (rho_point[star_i] > TOV_Atmosphere[star_i])) @@ -1005,6 +1131,11 @@ (1.0 + eps[i3D]) + press[i3D]*(w_lorentz[i3D]*w_lorentz[i3D]-1.0) ) - dens[i3D]; + + if(eoskey==4) { + Y_e_con[i3D] = Y_e[i3D] * dens[i3D]; + } + } /* if used, recalculate the derivatives of the conformal factor */ if (*conformal_state > 1) @@ -1095,6 +1226,10 @@ free(mu_point); free(phi_point); free(r_point); + free(temperature_point); + free(entropy_point); + free(ye_point); + } void TOV_Prepare_Fake_Evolution(CCTK_ARGUMENTS)

1 0

[svn:einsteintoolkit] incoming/TOVSolverHot/ (Rev. 26)
by cott＠tapir.caltech.edu 04 Dec '10

04 Dec '10

User: cott Date: 2010/12/04 04:52 PM Modified: /TOVSolverHot/ interface.ccl, param.ccl /TOVSolverHot/src/ external.inc, tov.c Log: * this now works for polytropic NSs using EOS omni. File Changes: Directory: /TOVSolverHot/ ========================= File [modified]: interface.ccl Delta lines: +14 -12 =================================================================== --- TOVSolverHot/interface.ccl 2010-12-04 21:53:41 UTC (rev 25) +++ TOVSolverHot/interface.ccl 2010-12-04 22:52:47 UTC (rev 26) @@ -91,17 +91,19 @@ USES FUNCTION EOS_Omni_press -void FUNCTION EOS_Omni_RhoFromPress(CCTK_INT IN eoskey, \ - CCTK_INT IN havetemp, \ - CCTK_REAL IN rf_precision, \ - CCTK_INT IN npoints, \ - CCTK_REAL INOUT ARRAY rho, \ - CCTK_REAL INOUT ARRAY eps, \ - CCTK_REAL IN ARRAY temp, \ - CCTK_REAL IN ARRAY ye, \ - CCTK_REAL IN ARRAY press, \ - CCTK_REAL OUT ARRAY xeps, \ - CCTK_INT OUT ARRAY keyerr, \ + +void FUNCTION EOS_Omni_RhoFromPressEpsTempEnt( \ + CCTK_INT IN eoskey, \ + CCTK_INT IN havetemp, \ + CCTK_REAL IN rf_precision, \ + CCTK_INT IN npoints, \ + CCTK_REAL OUT ARRAY rho, \ + CCTK_REAL INOUT ARRAY eps, \ + CCTK_REAL INOUT ARRAY temp, \ + CCTK_REAL INOUT ARRAY entropy, \ + CCTK_REAL IN ARRAY ye, \ + CCTK_REAL IN ARRAY press, \ + CCTK_INT OUT ARRAY keyerr, \ CCTK_INT OUT anyerr) -USES FUNCTION EOS_Omni_RhoFromPress +USES FUNCTION EOS_Omni_RhoFromPressEpsTempEnt File [modified]: param.ccl Delta lines: +41 -0 =================================================================== --- TOVSolverHot/param.ccl 2010-12-04 21:53:41 UTC (rev 25) +++ TOVSolverHot/param.ccl 2010-12-04 22:52:47 UTC (rev 26) @@ -131,6 +131,47 @@ ".*" :: "Any filename, not used if empty" } "" +######################################################## +# EOS Omni Parameters +######################################################## +BOOLEAN use_EOS_Omni "should we use omni?" +{ +} no + +CCTK_REAL TOVSolverHot_eos_rf_prec "precision for root finder" +{ + 1.0e-12:* :: "anything above crazy" +} 1.0e-10 + +CCTK_INT eoskey "Which EOS are we using" +{ + 1 :: "Polytropic" + 4 :: "nuc_eos" +} 1 + +CCTK_INT nuc_eos_keytemp "When using nuc_eos: fix temperature or entropy?" +{ + 1 :: "temperature" + 2 :: "entropy" +} 1 + +CCTK_REAL nuc_eos_temperature "Temperature to use when using nuc_eos [MeV]" +{ + 0.1:* :: "Below 0.1 MeV tables get bad" +} 1.0e0 + +CCTK_REAL nuc_eos_entropy "Entropy to use when using nuc_eos [k_B/baryon]" +{ + 0.3:* :: "Below 0.3 k_B/baryon things get dicey" +} 1.0e0 + + + + +######################################################## +# Other Parameters +######################################################## + shares:HydroBase EXTENDS KEYWORD initial_hydro "" Directory: /TOVSolverHot/src/ ============================= File [modified]: external.inc Delta lines: +21 -0 =================================================================== --- TOVSolverHot/src/external.inc 2010-12-04 21:53:41 UTC (rev 25) +++ TOVSolverHot/src/external.inc 2010-12-04 22:52:47 UTC (rev 26) @@ -65,6 +65,9 @@ /* TOV_debug_input_points(size,x,y,z); */ /* loop over all points we got */ + if(use_EOS_Omni) + CCTK_WARN(0,"TOV_Set_Rho_ADM does not work with EOS_Omni"); + #pragma omp parallel for for (CCTK_INT i=0; i<size; i++) { @@ -139,6 +142,9 @@ /* CCTK_INFO("Request for Momentum-Sources from TOVSolver"); */ /* loop over all points we got */ + if(use_EOS_Omni) + CCTK_WARN(0,"TOV_Set_Momentum_Source does not work with EOS_Omni"); + for (i=0; i<size; i++) { /* calculate the distance to the star */ @@ -209,6 +215,9 @@ CCTK_INFO("Using initial guess from TOVSolver"); /* TOV_debug_input_points(size,x,y,z); */ + if(use_EOS_Omni) + CCTK_WARN(0,"TOV_Set_Initial_Guess_for_u does not work with EOS_Omni"); + /* loop over all points we got */ for (i=0; i<size; i++) { @@ -251,9 +260,14 @@ CCTK_REAL my_psi4, CCTK_REAL x, CCTK_REAL y, CCTK_REAL z, CCTK_REAL rho_orig) { + DECLARE_CCTK_PARAMETERS; CCTK_REAL lb, ub, f; lb=0.0; ub=2*rho_orig; + + if(use_EOS_Omni) + CCTK_WARN(0,"bisection does not work with EOS_Omni"); + do { *rhonew = (ub-lb)/2.; @@ -289,8 +303,12 @@ CCTK_REAL my_psi4, CCTK_REAL x, CCTK_REAL y, CCTK_REAL z, CCTK_REAL rho_orig) { + DECLARE_CCTK_PARAMETERS int count=0; CCTK_REAL d_w_lorentz_2, f, df; + if(use_EOS_Omni) + CCTK_WARN(0,"newton_raphson does not work with EOS_Omni"); + do { count++; @@ -371,6 +389,9 @@ static int debug = 1; FILE *debugfile = NULL; + if(use_EOS_Omni) + CCTK_WARN(0,"Rescaling Sources does not work with EOS_Omni"); + CCTK_INFO("Rescaling Sources"); /* get GRHydro variables */ File [modified]: tov.c Delta lines: +148 -27 =================================================================== --- TOVSolverHot/src/tov.c 2010-12-04 21:53:41 UTC (rev 25) +++ TOVSolverHot/src/tov.c 2010-12-04 22:52:47 UTC (rev 26) @@ -11,6 +11,7 @@ #include <cctk.h> #include <cctk_Arguments.h> +#include <cctk_Functions.h> #include <cctk_Parameters.h> #include "constants.h" @@ -131,6 +132,8 @@ CCTK_REAL press, rho, eps, mu, m, phi, mbary, rprop; CCTK_REAL log_rbar_over_r, r_minus_two_m; + DECLARE_CCTK_PARAMETERS; + LOCAL_TINY = 1.0e-35; PI=4.0*atan(1.0); @@ -143,8 +146,37 @@ mbary = old_data[4]; rprop = old_data[5]; - rho = pow(press / K, 1.0 / Gamma); - eps = press / (Gamma - 1.0) / rho; + if(!use_EOS_Omni) { + rho = pow(press / K, 1.0 / Gamma); + eps = press / (Gamma - 1.0) / rho; + } else { + // EOS_Omni code + CCTK_INT keytemp; + CCTK_INT keyerr = 0; + CCTK_INT anyerr = 0; + CCTK_INT n = 1; + CCTK_REAL xtemp; + CCTK_REAL xye; + CCTK_REAL xent; + if(eoskey==4) { + CCTK_WARN(0,"fix me"); + } else { + keytemp = 0; + } + EOS_Omni_RhoFromPressEpsTempEnt(eoskey, + keytemp, + TOVSolverHot_eos_rf_prec, + n, + &rho, + &eps, + &xtemp, + &xent, + &xye, + &press, + &keyerr, + &anyerr); + + } mu = rho * (1.0 + eps); r_minus_two_m = r - 2.0 * m; @@ -262,8 +294,35 @@ TOV_C_fill(&(TOV_rprop_1d[star_i]), TOV_Num_Radial, 0.0); /* set start values */ - TOV_press_1d[star_i] = TOV_K * - pow(rho_central, TOV_Gamma); + if(!use_EOS_Omni) { + TOV_press_1d[star_i] = TOV_K * + pow(rho_central, TOV_Gamma); + } else { + // EOS_Omni code + CCTK_INT keytemp; + CCTK_INT keyerr = 0; + CCTK_INT anyerr = 0; + CCTK_INT n = 1; + CCTK_REAL xeps,xtemp,xye,xent; + if(eoskey==4) { + CCTK_WARN(0,"fix me"); + } else { + keytemp = 0; + } + EOS_Omni_press(eoskey, + keytemp, + TOVSolverHot_eos_rf_prec, + n, + &rho_central, + &xeps, + &xtemp, + &xye, + &(TOV_press_1d[star_i]), + &keyerr, + &anyerr); + } + + /* TOV_r_1d [star_i] = LOCAL_TINY; TOV_rbar_1d [star_i] = LOCAL_TINY;*/ @@ -329,8 +388,35 @@ if (TOV_press_1d[i+1] < 0.0) TOV_press_1d[i+1] = 0.0; - local_rho = pow(TOV_press_1d[i+1] / TOV_K, 1.0 / TOV_Gamma); + if(!use_EOS_Omni) { + local_rho = pow(TOV_press_1d[i+1] / TOV_K, 1.0 / TOV_Gamma); + } else { + // EOS_Omni code + CCTK_INT keytemp; + CCTK_INT keyerr = 0; + CCTK_INT anyerr = 0; + CCTK_INT n = 1; + CCTK_REAL xeps,xtemp,xye,xent; + if(eoskey==4) { + CCTK_WARN(0,"fix me"); + } else { + keytemp = 0; + } + EOS_Omni_RhoFromPressEpsTempEnt(eoskey, + keytemp, + TOVSolverHot_eos_rf_prec, + n, + &local_rho, + &xeps, + &xtemp, + &xent, + &xye, + &(TOV_press_1d[i+1]), + &keyerr, + &anyerr); + } + /* scan for the surface */ if ( (local_rho <= 0.0) || (TOV_press_1d[i+1] <= 0.0) ) @@ -373,14 +459,13 @@ CCTK_INFO("Integrated TOV equation"); /* do some info */ CCTK_VInfo(CCTK_THORNSTRING, "Information about the TOVs used:"); - CCTK_VInfo("", "TOV radius mass bary_mass mass(g) cent.rho rho(cgi) K K(cgi) Gamma"); - for (i=0; i<TOV_Num_TOVs; i++) - if (fabs(TOV_Gamma - 2.0) < LOCAL_TINY) - CCTK_VInfo(""," %d %8g %8g %8g %8.3g %8g %8.3g %8g %8.3g %8g", + if(!use_EOS_Omni) { + CCTK_VInfo("", "TOV radius mass bary_mass cent.rho rho(cgs) K K(cgs) Gamma"); + for (i=0; i<TOV_Num_TOVs; i++) + CCTK_VInfo(""," %d %8g %8g %8g %8g %8.3g %8g %8.3g %8g", i+1, TOV_R_Surface[i], TOV_m_1d[(i+1)*TOV_Num_Radial-1], TOV_mbary_1d[(i+1)*TOV_Num_Radial-1], - TOV_m_1d[(i+1)*TOV_Num_Radial-1]*CONSTANT_Msolar_cgi, TOV_Rho_Central[i], TOV_Rho_Central[i]/pow(CONSTANT_G_cgi,3.0)/ pow(CONSTANT_Msolar_cgi,2.0)* @@ -390,16 +475,18 @@ pow(CONSTANT_Msolar_cgi,2.0)/ pow(CONSTANT_c_cgi,4.0), TOV_Gamma); - else - CCTK_VInfo(""," %d %8g %8g %8.3g %8g %8.3g %8g %8g", + } else { + CCTK_VInfo("", "TOV radius mass bary_mass cent.rho rho(cgs)"); + for (i=0; i<TOV_Num_TOVs; i++) + CCTK_VInfo(""," %d %8g %8g %8g %8.3g %8g", i+1, TOV_R_Surface[i], TOV_m_1d[(i+1)*TOV_Num_Radial-1], - TOV_m_1d[(i+1)*TOV_Num_Radial-1]*CONSTANT_Msolar_cgi, + TOV_mbary_1d[(i+1)*TOV_Num_Radial-1], TOV_Rho_Central[i], TOV_Rho_Central[i]/pow(CONSTANT_G_cgi,3.0)/ pow(CONSTANT_Msolar_cgi,2.0)* - pow(CONSTANT_c_cgi,6.0), - TOV_K, TOV_Gamma); + pow(CONSTANT_c_cgi,6.0)); + } } @@ -643,20 +730,54 @@ /* is some perturbation wanted? */ if (Perturb[star] == 0) - rho_point[star] = pow(press_point[star]/TOV_K, - 1.0/TOV_Gamma); + if(!use_EOS_Omni) { + rho_point[star] = pow(press_point[star]/TOV_K, + 1.0/TOV_Gamma); + } else { + // EOS_Omni code + CCTK_INT keytemp; + CCTK_INT keyerr = 0; + CCTK_INT anyerr = 0; + CCTK_INT n = 1; + CCTK_REAL xeps,xtemp,xye,xent; + if(eoskey==4) { + CCTK_WARN(0,"fix me"); + } else { + keytemp = 0; + } + EOS_Omni_RhoFromPressEpsTempEnt(eoskey, + keytemp, + TOVSolverHot_eos_rf_prec, + n, + &(rho_point[star]), + &(eps_point[star]), + &xtemp, + &xent, + &xye, + &(press_point[star]), + &keyerr, + &anyerr); + + } else - rho_point[star] = pow(press_point[star]/TOV_K, - 1.0/TOV_Gamma) * - (1.0 + - Pert_Amplitude[star] * - cos(PI/2.0 * r[i3D] / TOV_R_Surface[star])); + if(!use_EOS_Omni) { + rho_point[star] = pow(press_point[star]/TOV_K, + 1.0/TOV_Gamma) * + (1.0 + + Pert_Amplitude[star] * + cos(PI/2.0 * r[i3D] / TOV_R_Surface[star])); + } else { + CCTK_WARN(0, "Perturbation handling not implemented for EOS_Omni"); + } - if (rho_point[star] > local_tiny) - eps_point[star] = press_point[star] / (TOV_Gamma - 1.0) - / rho_point[star]; - else - eps_point[star] = 0.0; + if(!use_EOS_Omni) { + if (rho_point[star] > local_tiny) + eps_point[star] = press_point[star] / (TOV_Gamma - 1.0) + / rho_point[star]; + else + eps_point[star] = 0.0; + } // eps already handled by EOS_Omni call above + mu_point[star] = rho_point[star] * (1.0 + eps_point[star]); } /* find out from which star we want to have the data */

1 0

[svn:einsteintoolkit] incoming/TOVSolverHot/ (Rev. 25)
by cott＠tapir.caltech.edu 04 Dec '10

04 Dec '10

User: cott Date: 2010/12/04 03:53 PM Modified: /TOVSolverHot/ param.ccl /TOVSolverHot/src/ external.inc, tov.c Log: * there may be only one K and Gamma; Make parameters scalar. File Changes: Directory: /TOVSolverHot/ ========================= File [modified]: param.ccl Delta lines: +2 -2 =================================================================== --- TOVSolverHot/param.ccl 2010-12-04 21:45:57 UTC (rev 24) +++ TOVSolverHot/param.ccl 2010-12-04 21:53:41 UTC (rev 25) @@ -24,12 +24,12 @@ 0.0:* :: "Must be positive" } 0.0 -CCTK_REAL TOV_Gamma[10] "The polytropic constant in P = K rho^Gamma" +CCTK_REAL TOV_Gamma "The polytropic constant in P = K rho^Gamma" { 1.0: :: "The physical range at high Lorentz factors is [1,2], but otherwise higher values of gamma can also be used" } 2.0 -CCTK_REAL TOV_K[10] "The polytropic constant in P = K rho^Gamma" +CCTK_REAL TOV_K "The polytropic constant in P = K rho^Gamma" { (0.0:* :: "Greater than 0" } 100.0 Directory: /TOVSolverHot/src/ ============================= File [modified]: external.inc Delta lines: +36 -36 =================================================================== --- TOVSolverHot/src/external.inc 2010-12-04 21:45:57 UTC (rev 24) +++ TOVSolverHot/src/external.inc 2010-12-04 21:53:41 UTC (rev 25) @@ -87,8 +87,8 @@ &r_to_star, TOV_Surface[star], &press, &phi, &r); press = (press > 0.0) ? press : 0.0; - rho = pow(press/TOV_K[star], 1.0/TOV_Gamma[star]); - eps = (rho > 0.0) ? press/(TOV_Gamma[star] - 1.0) / rho : 0.0; + rho = pow(press/TOV_K, 1.0/TOV_Gamma); + eps = (rho > 0.0) ? press/(TOV_Gamma - 1.0) / rho : 0.0; gxx = r / (r_to_star + 1.0e-30) * r / (r_to_star + 1.0e-30); /* velocity components as in gr-qc/9811015 */ w_lorentz_2 = 1. / (1. - gxx*TOV_Velocity_x[star]*TOV_Velocity_x[star] @@ -97,7 +97,7 @@ if (rho > 0.0) /* source[i]=gxx*gxx*rho*(1.0+eps); */ source[i]=gxx*gxx* (w_lorentz_2*rho + - (w_lorentz_2*TOV_Gamma[star]/(TOV_Gamma[star]-1.)-1.)* + (w_lorentz_2*TOV_Gamma/(TOV_Gamma-1.)-1.)* press); else source[i]=0.0; @@ -156,14 +156,14 @@ &(TOV_rbar_1d[star_i]), &(TOV_r_1d[star_i]), &r_to_star, TOV_Surface[star], &press, &phi, &r); - rho = pow(press/TOV_K[star], 1.0/TOV_Gamma[star]); + rho = pow(press/TOV_K, 1.0/TOV_Gamma); gxx = r / (r_to_star + 1.0e-30) * r / (r_to_star + 1.0e-30); psip = pow(gxx, TOV_Momentum_Psi_Power/4.); /* velocity components as in gr-qc/9811015 */ w_lorentz_2 = 1./(1. - gxx*TOV_Velocity_x[star]*TOV_Velocity_x[star] - gxx*TOV_Velocity_y[star]*TOV_Velocity_y[star] - gxx*TOV_Velocity_z[star]*TOV_Velocity_z[star]); - rho_eps = w_lorentz_2 * (rho + TOV_Gamma[star]/(TOV_Gamma[star]-1.) * press); + rho_eps = w_lorentz_2 * (rho + TOV_Gamma/(TOV_Gamma-1.) * press); switch(dir) { case 0: source[i]=psip*rho_eps*gxx*TOV_Velocity_x[star]; break; @@ -228,9 +228,9 @@ &(TOV_rbar_1d[star_i]), &(TOV_r_1d[star_i]), &r_to_star, TOV_Surface[star], &press, &phi, &r); - rho = pow(press/TOV_K[star], 1.0/TOV_Gamma[star]); + rho = pow(press/TOV_K, 1.0/TOV_Gamma); if (rho>0.0) - eps = press/(TOV_Gamma[star] - 1.0) / rho; + eps = press/(TOV_Gamma - 1.0) / rho; else eps = 0.0; gxx = r / (r_to_star + 1.0e-30) * r / (r_to_star + 1.0e-30); @@ -460,11 +460,11 @@ CCTK_REAL f, df; count++; vx = mom_source[0][i3D]/my_psi4/ - (source[i3D]/my_psi4/my_psi4+TOV_K[0] * pow(rhonew, TOV_Gamma[0])); + (source[i3D]/my_psi4/my_psi4+TOV_K * pow(rhonew, TOV_Gamma[0])); vy = mom_source[1][i3D]/my_psi4/ - (source[i3D]/my_psi4/my_psi4+TOV_K[0] * pow(rhonew, TOV_Gamma[0])); + (source[i3D]/my_psi4/my_psi4+TOV_K * pow(rhonew, TOV_Gamma[0])); vz = mom_source[2][i3D]/my_psi4/ - (source[i3D]/my_psi4/my_psi4+TOV_K[0] * pow(rhonew, TOV_Gamma[0])); + (source[i3D]/my_psi4/my_psi4+TOV_K * pow(rhonew, TOV_Gamma[0])); v_2 = my_psi4*(vx*vx+vy*vy+vz*vz); if (count > 100) CCTK_VWarn(count<110, __LINE__, __FILE__, CCTK_THORNSTRING, \ @@ -475,36 +475,36 @@ /* velocity components as in gr-qc/9811015 */ w_lorentz_2 = 1./(1.- v_2); /* - f = TOV_K[0] / (TOV_Gamma[0]-1.0) * pow(rhonew, TOV_Gamma[0]) + + f = TOV_K / (TOV_Gamma-1.0) * pow(rhonew, TOV_Gamma) + rhonew - source[i3D]; - df = TOV_K[0] * TOV_Gamma[0] / (TOV_Gamma[0]-1.0) * - pow(rhonew, TOV_Gamma[0]-1.0) + 1.0; + df = TOV_K * TOV_Gamma / (TOV_Gamma-1.0) * + pow(rhonew, TOV_Gamma-1.0) + 1.0; */ f = my_psi4*my_psi4*( w_lorentz_2 * rhonew + - (w_lorentz_2*TOV_Gamma[0]/(TOV_Gamma[0]-1.)-1.) * - TOV_K[0] * pow(rhonew, TOV_Gamma[0])) - source[i3D]; + (w_lorentz_2*TOV_Gamma/(TOV_Gamma-1.)-1.) * + TOV_K * pow(rhonew, TOV_Gamma)) - source[i3D]; /* df= w_lorentz_2 + - (w_lorentz_2*TOV_Gamma[0]/(TOV_Gamma[0]-1.)-1.) * - TOV_K[0] * TOV_Gamma[0] * pow(rhonew, TOV_Gamma[0]-1.); + (w_lorentz_2*TOV_Gamma/(TOV_Gamma-1.)-1.) * + TOV_K * TOV_Gamma * pow(rhonew, TOV_Gamma-1.); */ /* d_w_lorentz_2/drhonew */ CCTK_REAL d_w_lorentz_2 = -2 * w_lorentz_2*w_lorentz_2 * - v_2 * TOV_K[0]*TOV_Gamma[0]*pow(rhonew, TOV_Gamma[0]-1) / - (source[i3D]/my_psi4/my_psi4 + TOV_K[0]*pow(rhonew, TOV_Gamma[0])); + v_2 * TOV_K*TOV_Gamma*pow(rhonew, TOV_Gamma-1) / + (source[i3D]/my_psi4/my_psi4 + TOV_K*pow(rhonew, TOV_Gamma)); df = my_psi4*my_psi4*( d_w_lorentz_2*rhonew + w_lorentz_2 + - d_w_lorentz_2*TOV_Gamma[0]/(TOV_Gamma[0]-1.)* - TOV_K[0] * pow(rhonew, TOV_Gamma[0]) + - (w_lorentz_2*TOV_Gamma[0]/(TOV_Gamma[0]-1.)-1.) * - TOV_K[0] * TOV_Gamma[0] * pow(rhonew, TOV_Gamma[0]-1.)); + d_w_lorentz_2*TOV_Gamma/(TOV_Gamma-1.)* + TOV_K * pow(rhonew, TOV_Gamma) + + (w_lorentz_2*TOV_Gamma/(TOV_Gamma-1.)-1.) * + TOV_K * TOV_Gamma * pow(rhonew, TOV_Gamma-1.)); /* df= 1. + my_psi4*v_2/(rhonew*rhonew) + - ((TOV_Gamma[0]+(2.-TOV_Gamma[0])*my_psi4*v_2/(rhonew*rhonew)) - / (TOV_Gamma[0]-1.) - 1.) - * TOV_K[0] * TOV_Gamma[0] * pow(rhonew, TOV_Gamma[0]-1.); + ((TOV_Gamma+(2.-TOV_Gamma)*my_psi4*v_2/(rhonew*rhonew)) + / (TOV_Gamma-1.) - 1.) + * TOV_K * TOV_Gamma * pow(rhonew, TOV_Gamma-1.); */ rhoold = rhonew; @@ -516,22 +516,22 @@ if (count==1) printf("Fast rescale! %g %g %g %g\n", my_psi4, my_psi4, mom_source[0][i3D]/my_psi4/ - (source[i3D]/my_psi4/my_psi4+TOV_K[0]* - pow(rhonew, TOV_Gamma[0])), + (source[i3D]/my_psi4/my_psi4+TOV_K* + pow(rhonew, TOV_Gamma)), mom_source[0][i3D]/my_psi4/ - (source[i3D]/my_psi4/my_psi4+TOV_K[0]* - pow(rhonew, TOV_Gamma[0]))); + (source[i3D]/my_psi4/my_psi4+TOV_K* + pow(rhonew, TOV_Gamma))); #else newton_raphson(&vx,&vy,&vz,&rhoold,&rhonew,&w_lorentz_2,&v_2, - TOV_Gamma[0], TOV_K[0], + TOV_Gamma, TOV_K, source[i3D], mom_source[0][i3D],mom_source[1][i3D],mom_source[2][i3D], my_psi4, x[i3D], y[i3D], z[i3D], rho[i3D]); #endif rho[i3D] = rhonew; - press[i3D] = TOV_K[0] * pow(rhonew, TOV_Gamma[0]); - eps[i3D] = TOV_K[0] * pow(rhonew, TOV_Gamma[0] - 1.0) / - (TOV_Gamma[0] - 1.0); + press[i3D] = TOV_K * pow(rhonew, TOV_Gamma); + eps[i3D] = TOV_K * pow(rhonew, TOV_Gamma - 1.0) / + (TOV_Gamma - 1.0); sqrt_det = pow(my_psi4, 1.5); w_lorentz[i3D] = sqrt(w_lorentz_2); dens[i3D] = sqrt_det * w_lorentz[i3D] * rho[i3D]; @@ -558,11 +558,11 @@ "%.15g %.15g %.15g %.15g\n", /*7*/ gxx[i3D]*gxx[i3D]*( /* ham - source*/ w_lorentz_2*rho[i3D]+ - (w_lorentz_2*TOV_Gamma[0]/(TOV_Gamma[0]-1.)-1.)*press[i3D]), + (w_lorentz_2*TOV_Gamma/(TOV_Gamma-1.)-1.)*press[i3D]), gxx[i3D]*gxx[i3D]*( /* ham - source */ w_lorentz_2*(rho[i3D]*(1.+eps[i3D])+press[i3D])-press[i3D]), (w_lorentz_2 * rho[i3D] + /* momx - source */ - w_lorentz_2*TOV_Gamma[0]/(TOV_Gamma[0]-1.) * press[i3D])* + w_lorentz_2*TOV_Gamma/(TOV_Gamma-1.) * press[i3D])* gxx[i3D]*vel0[i3D], w_lorentz_2 * (rho[i3D]* /* momx - source */ (1.+eps[i3D])+ press[i3D])* File [modified]: tov.c Delta lines: +17 -17 =================================================================== --- TOVSolverHot/src/tov.c 2010-12-04 21:45:57 UTC (rev 24) +++ TOVSolverHot/src/tov.c 2010-12-04 21:53:41 UTC (rev 25) @@ -262,8 +262,8 @@ TOV_C_fill(&(TOV_rprop_1d[star_i]), TOV_Num_Radial, 0.0); /* set start values */ - TOV_press_1d[star_i] = TOV_K[star] * - pow(rho_central, TOV_Gamma[star]); + TOV_press_1d[star_i] = TOV_K * + pow(rho_central, TOV_Gamma); /* TOV_r_1d [star_i] = LOCAL_TINY; TOV_rbar_1d [star_i] = LOCAL_TINY;*/ @@ -295,25 +295,25 @@ TOV_C_fill(source_data, 6, 0.0); TOV_C_Source_RHS(TOV_r_1d[i], - TOV_K[star], TOV_Gamma[star], + TOV_K, TOV_Gamma, in_data, source_data); RKLOOP k1[rk] = TOV_dr[star] * source_data[rk]; RKLOOP in_data[rk] = old_data[rk] + 0.5 * k1[rk]; TOV_C_Source_RHS(TOV_r_1d[i]+ 0.5 * TOV_dr[star], - TOV_K[star], TOV_Gamma[star], + TOV_K, TOV_Gamma, in_data, source_data); RKLOOP k2[rk] = TOV_dr[star] * source_data[rk]; RKLOOP in_data[rk] = old_data[rk] + 0.5 * k2[rk]; TOV_C_Source_RHS(TOV_r_1d[i]+ 0.5 * TOV_dr[star], - TOV_K[star], TOV_Gamma[star], + TOV_K, TOV_Gamma, in_data, source_data); RKLOOP k3[rk] = TOV_dr[star] * source_data[rk]; RKLOOP in_data[rk] = old_data[rk] + k3[rk]; TOV_C_Source_RHS(TOV_r_1d[i]+ TOV_dr[star], - TOV_K[star], TOV_Gamma[star], + TOV_K, TOV_Gamma, in_data, source_data); RKLOOP k4[rk] = TOV_dr[star] * source_data[rk]; RKLOOP new_data[rk] = old_data[rk] + (k1[rk] + k4[rk] + 2.0 * (k2[rk] + k3[rk])) /6.0; @@ -329,7 +329,7 @@ if (TOV_press_1d[i+1] < 0.0) TOV_press_1d[i+1] = 0.0; - local_rho = pow(TOV_press_1d[i+1] / TOV_K[star], 1.0 / TOV_Gamma[star]); + local_rho = pow(TOV_press_1d[i+1] / TOV_K, 1.0 / TOV_Gamma); /* scan for the surface */ if ( (local_rho <= 0.0) || @@ -375,7 +375,7 @@ CCTK_VInfo(CCTK_THORNSTRING, "Information about the TOVs used:"); CCTK_VInfo("", "TOV radius mass bary_mass mass(g) cent.rho rho(cgi) K K(cgi) Gamma"); for (i=0; i<TOV_Num_TOVs; i++) - if (fabs(TOV_Gamma[i] - 2.0) < LOCAL_TINY) + if (fabs(TOV_Gamma - 2.0) < LOCAL_TINY) CCTK_VInfo(""," %d %8g %8g %8g %8.3g %8g %8.3g %8g %8.3g %8g", i+1, TOV_R_Surface[i], TOV_m_1d[(i+1)*TOV_Num_Radial-1], @@ -385,11 +385,11 @@ TOV_Rho_Central[i]/pow(CONSTANT_G_cgi,3.0)/ pow(CONSTANT_Msolar_cgi,2.0)* pow(CONSTANT_c_cgi,6.0), - TOV_K[i], - TOV_K[i]*pow(CONSTANT_G_cgi,3.0)* + TOV_K, + TOV_K*pow(CONSTANT_G_cgi,3.0)* pow(CONSTANT_Msolar_cgi,2.0)/ pow(CONSTANT_c_cgi,4.0), - TOV_Gamma[i]); + TOV_Gamma); else CCTK_VInfo(""," %d %8g %8g %8.3g %8g %8.3g %8g %8g", i+1, TOV_R_Surface[i], @@ -399,7 +399,7 @@ TOV_Rho_Central[i]/pow(CONSTANT_G_cgi,3.0)/ pow(CONSTANT_Msolar_cgi,2.0)* pow(CONSTANT_c_cgi,6.0), - TOV_K[i], TOV_Gamma[i]); + TOV_K, TOV_Gamma); } @@ -643,17 +643,17 @@ /* is some perturbation wanted? */ if (Perturb[star] == 0) - rho_point[star] = pow(press_point[star]/TOV_K[star], - 1.0/TOV_Gamma[star]); + rho_point[star] = pow(press_point[star]/TOV_K, + 1.0/TOV_Gamma); else - rho_point[star] = pow(press_point[star]/TOV_K[star], - 1.0/TOV_Gamma[star]) * + rho_point[star] = pow(press_point[star]/TOV_K, + 1.0/TOV_Gamma) * (1.0 + Pert_Amplitude[star] * cos(PI/2.0 * r[i3D] / TOV_R_Surface[star])); if (rho_point[star] > local_tiny) - eps_point[star] = press_point[star] / (TOV_Gamma[star] - 1.0) + eps_point[star] = press_point[star] / (TOV_Gamma - 1.0) / rho_point[star]; else eps_point[star] = 0.0;

1 0

[svn:einsteintoolkit] incoming/TOVSolverHot/ (Rev. 24)
by cott＠tapir.caltech.edu 04 Dec '10

04 Dec '10

User: cott Date: 2010/12/04 03:45 PM Added: /TOVSolverHot/ COPYING, GPLv2, README, configuration.ccl, interface.ccl, param.ccl, schedule.ccl /TOVSolverHot/doc/ Baumgarte-notes-on-TOV.ps.gz, documentation.tex /TOVSolverHot/src/ external.inc, make.code.defn, make.code.deps, tov.c, tov.h, utils.inc /TOVSolverHot/test/ Log: * add a hot TOV solver File Changes: Directory: /TOVSolverHot/ ========================= File [added]: COPYING Delta lines: <Unavailable> File [added]: GPLv2 Delta lines: <Unavailable> File [added]: README Delta lines: <Unavailable> File [added]: configuration.ccl Delta lines: <Unavailable> File [added]: interface.ccl Delta lines: <Unavailable> File [added]: param.ccl Delta lines: <Unavailable> File [added]: schedule.ccl Delta lines: <Unavailable> Directory: /TOVSolverHot/doc/ ============================= File [added]: Baumgarte-notes-on-TOV.ps.gz Delta lines: <Unavailable> File [added]: documentation.tex Delta lines: <Unavailable> Directory: /TOVSolverHot/src/ ============================= File [added]: external.inc Delta lines: <Unavailable> File [added]: make.code.defn Delta lines: <Unavailable> File [added]: make.code.deps Delta lines: <Unavailable> File [added]: tov.c Delta lines: <Unavailable> File [added]: tov.h Delta lines: <Unavailable> File [added]: utils.inc Delta lines: <Unavailable>

1 0

Jump to page:

2026

2025

2024

2023

2022

2021

2020

2019

2018

2017

2016

2015

2014

2013

2012

2011

2010

Commits