Commits

commits@lists.einsteintoolkit.org

2 participants
6997 discussions

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 132)
by cott＠tapir.caltech.edu 02 Jul '10

02 Jul '10

User: cott Date: 2010/07/02 04:57 PM Modified: /trunk/src/ GRHydro_Con2Prim.F90, GRHydro_Eigenproblem.F90, GRHydro_Eigenproblem_Marquina.F90, GRHydro_EoSChangeGamma.F90, GRHydro_Marquina.F90, GRHydro_Prim2Con.F90 Log: * Changes to incorporate the EOS_Omni interace. GRHydro has now 3 EOS interfaces. Rejoice! To get the new one compiled in, use -DUSE_EOS_OMNI in CPPFLAGS and FPPFLAGS. File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_Con2Prim.F90 Delta lines: +233 -16 =================================================================== --- trunk/src/GRHydro_Con2Prim.F90 2010-07-02 21:51:08 UTC (rev 131) +++ trunk/src/GRHydro_Con2Prim.F90 2010-07-02 21:57:07 UTC (rev 132) @@ -150,6 +150,21 @@ CCTK_REAL :: local_min_tracer +#if USE_EOS_OMNI +! begin EOS Omni vars + integer :: n = 1 + integer :: poly_eoskey = 0 + integer :: keytemp = 0 + integer :: anyerr = 0 + integer :: keyerr(1) = 0 + real*8 :: xpress = 0.0d0 + real*8 :: xeps = 0.0d0 + real*8 :: xtemp = 0.0d0 + real*8 :: xye = 0.0d0 + poly_eoskey = GRHydro_poly_eoskey +! end EOS Omni vars +#endif + call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere") call SpaceMask_GetStateBits(atmosphere, "Hydro_Atmosphere", "in_atmosphere") type2_bits = -1 @@ -163,9 +178,17 @@ else local_min_tracer = 0d0 end if - + +#if USE_EOS_OMNI + ! this is a poly call + call EOS_Omni_press(poly_eoskey,keytemp,n,& + GRHydro_rho_min,xeps,xtemp,xye,pmin,keyerr,anyerr) + call EOS_Omni_EpsFromPress(poly_eoskey,keytemp,n,& + GRHydro_rho_min,xeps,xtemp,xye,pmin,epsmin,keyerr,anyerr) +#else pmin = EOS_Pressure(GRHydro_polytrope_handle, GRHydro_rho_min, 1.0d0) epsmin = EOS_SpecificIntEnergy(GRHydro_polytrope_handle, GRHydro_rho_min, pmin) +#endif !$OMP PARALLEL DO PRIVATE(i,j,itracer,& !$OMP uxx, uxy, uxz, uyy, uyz, uzz, det, psi4pt, epsnegative) @@ -221,8 +244,16 @@ scon(i,j,k,:) = 0.d0 vel(i,j,k,:) = 0.d0 w_lorentz(i,j,k) = 1.d0 +#if USE_EOS_OMNI + call EOS_Omni_press(poly_eoskey,keytemp,n,& + rho(i,j,k),eps(i,j,k),xtemp,xye,press(i,j,k),keyerr,anyerr) + + call EOS_Omni_EpsFromPress(poly_eoskey,keytemp,n,& + rho(i,j,k),eps(i,j,k),xtemp,xye,press(i,j,k),eps(i,j,k),keyerr,anyerr) +#else press(i,j,k) = EOS_Pressure(GRHydro_polytrope_handle, rho(i,j,k), eps(i,j,k)) eps(i,j,k) = EOS_SpecificIntEnergy(GRHydro_polytrope_handle, rho(i,j,k), press(i,j,k)) +#endif ! w_lorentz=1, so the expression for tau reduces to: tau(i,j,k) = sqrt(det) * (rho(i,j,k)+rho(i,j,k)*eps(i,j,k)) - dens(i,j,k) @@ -235,9 +266,9 @@ scon(i,j,k,3),tau(i,j,k),rho(i,j,k),vel(i,j,k,1),vel(i,j,k,2), & vel(i,j,k,3),eps(i,j,k),press(i,j,k),w_lorentz(i,j,k), & uxx,uxy,uxz,uyy,uyz,uzz,det,x(i,j,k),y(i,j,k), & - z(i,j,k),r(i,j,k),epsnegative,GRHydro_rho_min,pmin, epsmin, GRHydro_reflevel, GRHydro_C2P_failed(i,j,k)) + z(i,j,k),r(i,j,k),epsnegative,GRHydro_rho_min,pmin, epsmin, & + GRHydro_reflevel, GRHydro_C2P_failed(i,j,k)) - if (epsnegative) then #if 1 @@ -267,8 +298,16 @@ scon(i,j,k,:) = 0.d0 vel(i,j,k,:) = 0.d0 w_lorentz(i,j,k) = 1.d0 - press(i,j,k) = EOS_Pressure(GRHydro_polytrope_handle, rho(i,j,k), eps(i,j,k)) - eps(i,j,k) = EOS_SpecificIntEnergy(GRHydro_polytrope_handle, rho(i,j,k), press(i,j,k)) +#if USE_EOS_OMNI + call EOS_Omni_press(poly_eoskey,keytemp,n,& + rho(i,j,k),eps(i,j,k),xtemp,xye,press(i,j,k),keyerr,anyerr) + + call EOS_Omni_EpsFromPress(poly_eoskey,keytemp,n,& + rho(i,j,k),eps(i,j,k),xtemp,xye,press(i,j,k),eps(i,j,k),keyerr,anyerr) +#else + press(i,j,k) = EOS_Pressure(GRHydro_polytrope_handle, rho(i,j,k), eps(i,j,k)) + eps(i,j,k) = EOS_SpecificIntEnergy(GRHydro_polytrope_handle, rho(i,j,k), press(i,j,k)) +#endif ! w_lorentz=1, so the expression for tau reduces to: tau(i,j,k) = sqrt(det) * (rho(i,j,k)+rho(i,j,k)*eps(i,j,k)) - dens(i,j,k) @@ -315,12 +354,12 @@ subroutine Con2Prim_pt(handle, dens, sx, sy, sz, tau, rho, velx, vely, & velz, epsilon, press, w_lorentz, uxx, uxy, uxz, uyy, & - uyz, uzz, det, x, y, z, r, epsnegative, GRHydro_rho_min, pmin, epsmin, GRHydro_reflevel, GRHydro_C2P_failed) + uyz, uzz, det, x, y, z, r, epsnegative, GRHydro_rho_min, pmin, epsmin, & + GRHydro_reflevel, GRHydro_C2P_failed) implicit none DECLARE_CCTK_PARAMETERS - CCTK_REAL dens, sx, sy, sz, tau, rho, velx, vely, velz, epsilon, & press, uxx, uxy, uxz, uyy, uyz, uzz, det, w_lorentz, x, & y, z, r, GRHydro_rho_min @@ -336,6 +375,22 @@ #undef _EOS_BASE_INC_ #endif #include "EOS_Base.inc" + + +#if USE_EOS_OMNI +! begin EOS Omni vars + integer :: n = 1 + integer :: eoskey = 0 + integer :: keytemp = 0 + integer :: anyerr = 0 + integer :: keyerr(1) = 0 + real*8 :: xpress = 0.0d0 + real*8 :: xtemp = 0.0d0 + real*8 :: xye = 0.0d0 + eoskey = GRHydro_eoskey +! end EOS Omni vars +#endif + !!$ Undensitize the variables @@ -348,9 +403,15 @@ 2.*usx*usz*uxz + 2.*usy*usz*uyz !!$ Set initial guess for pressure: - +#if USE_EOS_OMNI + call EOS_Omni_press(eoskey,keytemp,n,& + rho,epsilon,xtemp,xye,xpress,keyerr,anyerr) + pold = max(1.d-10,xpress) +#else pold = max(1.d-10,EOS_Pressure(handle, rho, epsilon)) +#endif + !!$ Check that the variables have a chance of being physical if( (utau + pold + udens)**2 - s2 .le. 0.0d0) then @@ -386,8 +447,16 @@ !!$ Calculate the function +#if USE_EOS_OMNI + call EOS_Omni_press(eoskey,keytemp,n,& + rho,epsilon,xtemp,xye,xpress,keyerr,anyerr) + + f = pold - xpress +#else f = pold - EOS_Pressure(handle, rho, epsilon) +#endif + !!$Find the root count = 0 @@ -430,8 +499,18 @@ goto 51 end if +#if USE_EOS_OMNI + call EOS_Omni_DPressByDRho(eoskey,keytemp,n,& + rho,epsilon,xtemp,xye,dpressbydrho,keyerr,anyerr) + + call EOS_Omni_DPressByDEps(eoskey,keytemp,n,& + rho,epsilon,xtemp,xye,dpressbydeps,keyerr,anyerr) +#else dpressbydrho = EOS_DPressByDRho(handle,rho,epsilon) dpressbydeps = EOS_DPressByDEps(handle,rho,epsilon) +#endif + + temp1 = (utau+udens+pnew)**2 - s2 drhobydpress = udens * s2 / (sqrt(temp1)*(udens+utau+pnew)**2) depsbydpress = pnew * s2 / (rho * (udens + utau + pnew) * temp1) @@ -449,7 +528,14 @@ epsilon = (sqrt( (utau + pnew + udens)**2 - s2) - pnew * w_lorentz - & udens)/udens +#if USE_EOS_OMNI + call EOS_Omni_press(eoskey,keytemp,n,& + rho,epsilon,xtemp,xye,xpress,keyerr,anyerr) + + f = pnew - xpress +#else f = pnew - EOS_Pressure(handle, rho, epsilon) +#endif !check whether rho is NaN if (rho .ne. rho) then @@ -466,8 +552,17 @@ do i=1,GRHydro_polish +#if USE_EOS_OMNI + call EOS_Omni_DPressByDRho(eoskey,keytemp,n,& + rho,epsilon,xtemp,xye,dpressbydrho,keyerr,anyerr) + + call EOS_Omni_DPressByDEps(eoskey,keytemp,n,& + rho,epsilon,xtemp,xye,dpressbydeps,keyerr,anyerr) +#else dpressbydrho = EOS_DPressByDRho(handle,rho,epsilon) dpressbydeps = EOS_DPressByDEps(handle,rho,epsilon) +#endif + temp1 = (utau+udens+pnew)**2 - s2 drhobydpress = udens * s2 / (sqrt(temp1)*(udens+utau+pnew)**2) depsbydpress = pnew * s2 / (rho * (udens + utau + pnew) * temp1) @@ -484,8 +579,16 @@ epsilon = (sqrt( (utau + pnew + udens)**2 - s2) - pnew * w_lorentz - & udens)/udens +#if USE_EOS_OMNI + call EOS_Omni_press(eoskey,keytemp,n, & + rho,epsilon,xtemp,xye,xpress,keyerr,anyerr) + + f = pold - xpress +#else f = pold - EOS_Pressure(handle, rho, epsilon) +#endif + enddo !!$ BEGIN: Check for NaN value (2nd check) @@ -575,14 +678,42 @@ CCTK_INT :: type2_bits CCTK_REAL :: local_min_tracer - + #ifdef _EOS_BASE_INC_ #undef _EOS_BASE_INC_ #endif #include "EOS_Base.inc" + +#if USE_EOS_OMNI +! begin EOS omni + CCTK_INT :: eoskey = 0 + CCTK_INT :: poly_eoskey = 0 + CCTK_INT :: keyerr(1) = 0 + CCTK_INT :: anyerr = 0 + CCTK_INT :: keytemp = 0 + CCTK_INT :: n = 1 + CCTK_REAL :: xye = 0.0d0 + CCTK_REAL :: xeps = 0.0d0 + CCTK_REAL :: xtemp = 0.0d0 + eoskey = GRHydro_eoskey + poly_eoskey = GRHydro_poly_eoskey +! end EOS omni +#endif + + +#if USE_EOS_OMNI + ! this is a poly call + call EOS_Omni_press(poly_eoskey,keytemp,n,& + GRHydro_rho_min,1.0d0,xtemp,xye,pmin,keyerr,anyerr) + + call EOS_Omni_EpsFromPress(poly_eoskey,keytemp,n,& + GRHydro_rho_min,epsmin,xtemp,xye,pmin,epsmin,keyerr,anyerr) +#else pmin=EOS_Pressure(GRHydro_polytrope_handle, GRHydro_rho_min, 1.0d0) epsmin = EOS_SpecificIntEnergy(GRHydro_polytrope_handle, GRHydro_rho_min, pmin) +#endif + call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere") call SpaceMask_GetStateBits(atmosphere, "Hydro_Atmosphere", "in_atmosphere") @@ -807,6 +938,7 @@ local_min_tracer = 0d0 end if + !!$ do k = GRHydro_stencil + 1, nz - GRHydro_stencil !!$ do j = GRHydro_stencil + 1, ny - GRHydro_stencil !!$ do i = GRHydro_stencil + 1, nx - GRHydro_stencil @@ -907,6 +1039,22 @@ #endif #include "EOS_Base.inc" + +#if USE_EOS_OMNI +! begin EOS Omni vars + integer :: n = 1 + integer :: poly_eoskey = 0 + integer :: keytemp = 0 + integer :: anyerr = 0 + integer :: keyerr(1) = 0 + real*8 :: xpress = 0.0d0 + real*8 :: xeps = 0.0d0 + real*8 :: xtemp = 0.0d0 + real*8 :: xye = 0.0d0 + poly_eoskey = GRHydro_poly_eoskey +! end EOS Omni vars +#endif + !!$ Undensitize the variables sqrtdet = sqrt(det) @@ -919,10 +1067,21 @@ 2.d0*usx*usz*uxz + 2.d0*usy*usz*uyz !!$ Set initial guess for rho: + rhoold = max(GRHydro_rho_min,rho) - rhoold = max(GRHydro_rho_min,rho) +#if USE_EOS_OMNI +! call EOS_Omni_press(poly_eoskey,keytemp,n,& +! rhoold,xeps,xtemp,xye,xpress,keyerr,anyerr) + + call EOS_Omni_EpsFromPress(poly_eoskey,keytemp,n,& + rhoold,xeps,xtemp,xye,press,xeps,keyerr,anyerr) + + enthalpy = 1.0d0 + xeps + press / rhoold +#else enthalpy = 1.d0 + EOS_SpecificIntEnergy(handle, rhoold, press) + & EOS_Pressure(handle, rhoold, epsilon) / rhoold +#endif + w_lorentz = sqrt(1.d0 + s2 / ( (udens*enthalpy)**2 )) !!$ BEGIN: Check for NaN value (1st check) @@ -934,8 +1093,12 @@ endif !!$ END: Check for NaN value (1st check) +#if USE_EOS_OMNI + call EOS_Omni_press(poly_eoskey,keytemp,n,& + rhoold,epsilon,xtemp,xye,press,keyerr,anyerr) +#else press = EOS_Pressure(handle, rhoold, epsilon) - +#endif !!$ Calculate the function f = rhoold * w_lorentz - udens @@ -968,8 +1131,16 @@ rhonew = GRHydro_rho_min udens = rhonew dens = sqrtdet * rhonew +#if USE_EOS_OMNI + call EOS_Omni_press(poly_eoskey,keytemp,n,& + rhonew,1.0d0,xtemp,xye,press,keyerr,anyerr) + + call EOS_Omni_EpsFromPress(poly_eoskey,keytemp,n,& + rhonew,xeps,xtemp,xye,press,epsilon,keyerr,anyerr) +#else press = EOS_Pressure(handle, rhonew, 1.d0) epsilon = EOS_SpecificIntEnergy(handle, rhonew, press) +#endif sx = 0.d0 sy = 0.d0 sz = 0.d0 @@ -985,7 +1156,13 @@ !!$ Ian has the feeling that this is an error in general. But for the !!$ 2D_Polytrope it gives the right answers. +#if USE_EOS_OMNI + call EOS_Omni_DPressByDrho(poly_eoskey,keytemp,n,& + rhonew,epsilon,xtemp,xye,denthalpy,keyerr,anyerr) + denthalpy = denthalpy / rhonew +#else denthalpy = EOS_DPressByDrho(handle, rhonew, epsilon) / rhonew +#endif df = w_lorentz - rhonew * s2 * denthalpy / & (w_lorentz*(udens**2)*(enthalpy**3)) @@ -1005,12 +1182,22 @@ end if !!$ Recalculate primitive variables and function - - enthalpy = 1.d0 + EOS_SpecificIntEnergy(handle, rhonew, press) + & + +#if USE_EOS_OMNI + call EOS_Omni_press(poly_eoskey,keytemp,n,& + rhonew,epsilon,xtemp,xye,press,keyerr,anyerr) + + call EOS_Omni_EpsFromPress(poly_eoskey,keytemp,n,& + rhonew,xeps,xtemp,xye,press,epsilon,keyerr,anyerr) + + enthalpy = 1.0d0 + epsilon + press / rhonew +#else + enthalpy = 1.d0 + EOS_SpecificIntEnergy(handle, rhonew, press) + & EOS_Pressure(handle, rhonew, epsilon) / rhonew + press = EOS_Pressure(handle, rhonew, epsilon) +#endif + w_lorentz = sqrt(1.d0 + s2 / ( (udens*enthalpy)**2 )) - - press = EOS_Pressure(handle, rhonew, epsilon) tau = sqrtdet * ((rho * enthalpy) * w_lorentz**2 - press) - dens f = rhonew * w_lorentz - udens @@ -1022,10 +1209,21 @@ if (rhonew .le. GRHydro_rho_min*(1.d0+GRHydro_atmo_tolerance) ) then rhonew = GRHydro_rho_min udens = rhonew - ! before 2009/10/07 dens was not reset and was used with its negative value below; this has since been changed, but the change produces significant changes in the results + ! before 2009/10/07 dens was not reset and was used with its negative + ! value below; this has since been changed, but the change produces + ! significant changes in the results dens = sqrtdet * rhonew + +#if USE_EOS_OMNI + call EOS_Omni_press(poly_eoskey,keytemp,n,& + rhonew,1.0d0,xtemp,xye,press,keyerr,anyerr) + + call EOS_Omni_EpsFromPress(poly_eoskey,keytemp,n,& + rhonew,xeps,xtemp,xye,press,epsilon,keyerr,anyerr) +#else press = EOS_Pressure(handle, rhonew, 1.d0) epsilon = EOS_SpecificIntEnergy(handle, rhonew, press) +#endif ! w_lorentz=1, so the expression for utau reduces to: tau = rhonew + rhonew*epsilon - udens sx = 0.d0 @@ -1040,9 +1238,19 @@ 50 rho = rhonew +#if USE_EOS_OMNI + call EOS_Omni_press(poly_eoskey,keytemp,n,& + rhonew,epsilon,xtemp,xye,press,keyerr,anyerr) + + call EOS_Omni_EpsFromPress(poly_eoskey,keytemp,n,& + rhonew,xeps,xtemp,xye,press,epsilon,keyerr,anyerr) + + enthalpy = 1.0d0 + epsilon + press / rhonew +#else enthalpy = 1.d0 + EOS_SpecificIntEnergy(handle, rhonew, press) + & EOS_Pressure(handle, rhonew, epsilon) / rhonew w_lorentz = sqrt(1.d0 + s2 / ( (udens*enthalpy)**2 )) +#endif !!$ BEGIN: Check for NaN value (2nd check) if (w_lorentz .ne. w_lorentz) then @@ -1053,8 +1261,16 @@ endif !!$ END: Check for NaN value (2nd check) +#if USE_EOS_OMNI + call EOS_Omni_press(poly_eoskey,keytemp,n,& + rhonew,epsilon,xtemp,xye,press,keyerr,anyerr) + + call EOS_Omni_EpsFromPress(poly_eoskey,keytemp,n,& + rhonew,xeps,xtemp,xye,press,epsilon,keyerr,anyerr) +#else press = EOS_Pressure(handle, rhonew, epsilon) epsilon = EOS_SpecificIntEnergy(handle, rhonew, press) +#endif tau = sqrtdet * ((rho * enthalpy) * w_lorentz**2 - press) - dens if (epsilon .lt. 0.0d0) then @@ -1116,6 +1332,7 @@ @@*/ + subroutine Con2PrimBoundsPolytype(CCTK_ARGUMENTS) implicit none File [modified]: GRHydro_Eigenproblem.F90 Delta lines: +92 -0 =================================================================== --- trunk/src/GRHydro_Eigenproblem.F90 2010-07-02 21:51:08 UTC (rev 131) +++ trunk/src/GRHydro_Eigenproblem.F90 2010-07-02 21:57:07 UTC (rev 132) @@ -63,16 +63,49 @@ #endif #include "EOS_Base.inc" +#if USE_EOS_OMNI +! begin EOS Omni vars + integer :: n = 1 + integer :: eoskey = 0 + integer :: keytemp = 0 + integer :: anyerr = 0 + integer :: keyerr(1) = 0 + real*8 :: xpress = 0.0d0 + real*8 :: xeps = 0.0d0 + real*8 :: xtemp = 0.0d0 + real*8 :: xye = 0.0d0 + eoskey = GRHydro_eoskey +! end EOS Omni vars +#endif + one = 1.0d0 two = 2.0d0 !!$ Set required fluid quantities + + +#if USE_EOS_OMNI +! call EOS_Omni_cs2(eoskey,keytemp,n,& +! rho,eps,xtemp,xye,cs2,keyerr,anyerr) + call EOS_Omni_press(eoskey,keytemp,n,& + rho,eps,xtemp,xye,press,keyerr,anyerr) + + call EOS_Omni_DPressByDEps(eoskey,keytemp,n,& + rho,eps,xtemp,xye,dpdeps,keyerr,anyerr) + + call EOS_Omni_DPressByDRho(eoskey,keytemp,n,& + rho,eps,xtemp,xye,dpdrho,keyerr,anyerr) + + cs2 = (dpdrho + press * dpdeps / (rho**2))/ & + (1.0d0 + eps + press/rho) +#else press = EOS_Pressure(handle,rho,eps) dpdrho = EOS_DPressByDRho(handle,rho,eps) dpdeps = EOS_DPressByDEps(handle,rho,eps) cs2 = (dpdrho + press * dpdeps / (rho**2))/ & (1.0d0 + eps + press/rho) +#endif vlowx = gxx*velx + gxy*vely + gxz*velz vlowy = gxy*velx + gyy*vely + gyz*velz @@ -100,6 +133,7 @@ lam(4) = lam3 lam(5) = lamp + end subroutine eigenvalues /*@@ @@ -157,16 +191,47 @@ #endif #include "EOS_Base.inc" +#if USE_EOS_OMNI +! begin EOS Omni vars + integer :: n = 1 + integer :: eoskey = 0 + integer :: keytemp = 0 + integer :: anyerr = 0 + integer :: keyerr(1) = 0 + real*8 :: xpress = 0.0d0 + real*8 :: xeps = 0.0d0 + real*8 :: xtemp = 0.0d0 + real*8 :: xye = 0.0d0 + eoskey = GRHydro_eoskey +! end EOS Omni vars +#endif + one = 1.0d0 two = 2.0d0 !!$ Set required fluid quantities +#if USE_EOS_OMNI + call EOS_Omni_press(eoskey,keytemp,n,& + rho,eps,xtemp,xye,press,keyerr,anyerr) + + call EOS_Omni_DPressByDEps(eoskey,keytemp,n,& + rho,eps,xtemp,xye,dpdeps,keyerr,anyerr) + + call EOS_Omni_DPressByDRho(eoskey,keytemp,n,& + rho,eps,xtemp,xye,dpdrho,keyerr,anyerr) + + cs2 = (dpdrho + press * dpdeps / (rho**2))/ & + (1.0d0 + eps + press/rho) +! call EOS_Omni_cs2(eoskey,keytemp,n,& +! rho,eps,xtemp,xye,cs2,keyerr,anyerr) +#else press = EOS_Pressure(handle,rho,eps) dpdrho = EOS_DPressByDRho(handle,rho,eps) dpdeps = EOS_DPressByDEps(handle,rho,eps) cs2 = (dpdrho + press * dpdeps / (rho**2))/ & (1.0d0 + eps + press/rho) +#endif ! if (cs2<0) cs2=0 ! this does not modify the roe crashing problem with shocktube enthalpy = one + eps + press / rho @@ -345,6 +410,7 @@ rflux(4) = abs(lam1) * du(4) + enthalpy * w * (vlowz * vxa) rflux(5) = abs(lam1) * du(5) + enthalpy * w * (velx * vxb + vxa) - vxa + else !!$Form Jacobian matrix in characteristic form from right eigenvectors. @@ -471,6 +537,10 @@ roeflux4 = rflux(4) roeflux5 = rflux(5) + if(roeflux1.ne.roeflux1) then + stop "bad herberd" + endif + end subroutine eigenproblem /*@@ @@ -522,16 +592,38 @@ #endif #include "EOS_Base.inc" +#if USE_EOS_OMNI +! begin EOS Omni vars + integer :: n = 1 + integer :: eoskey = 0 + integer :: keytemp = 0 + integer :: anyerr = 0 + integer :: keyerr(1) = 0 + real*8 :: xpress = 0.0d0 + real*8 :: xeps = 0.0d0 + real*8 :: xtemp = 0.0d0 + real*8 :: xye = 0.0d0 + eoskey = GRHydro_eoskey +! end EOS Omni vars +#endif + one = 1.0d0 two = 2.0d0 !!$ Set required fluid quantities +#if USE_EOS_OMNI + call EOS_Omni_press(eoskey,keytemp,n,& + rho,eps,xtemp,xye,press,keyerr,anyerr) + call EOS_Omni_cs2(eoskey,keytemp,n,& + rho,eps,xtemp,xye,cs2,keyerr,anyerr) +#else press = EOS_Pressure(handle,rho,eps) dpdrho = EOS_DPressByDRho(handle,rho,eps) dpdeps = EOS_DPressByDEps(handle,rho,eps) cs2 = (dpdrho + press * dpdeps / (rho**2))/ & (1.0d0 + eps + press/rho) +#endif ! if (cs2<0) cs2=0 ! this does not modify the roe crashing problem with shocktube enthalpy = one + eps + press / rho File [modified]: GRHydro_Eigenproblem_Marquina.F90 Delta lines: +44 -5 =================================================================== --- trunk/src/GRHydro_Eigenproblem_Marquina.F90 2010-07-02 21:51:08 UTC (rev 131) +++ trunk/src/GRHydro_Eigenproblem_Marquina.F90 2010-07-02 21:57:07 UTC (rev 132) @@ -99,19 +99,46 @@ #endif #include "EOS_Base.inc" +#if USE_EOS_OMNI +! begin EOS Omni vars + integer :: n = 1 + integer :: eoskey = 0 + integer :: keytemp = 0 + integer :: anyerr = 0 + integer :: keyerr(1) = 0 + real*8 :: xpress = 0.0d0 + real*8 :: xeps = 0.0d0 + real*8 :: xtemp = 0.0d0 + real*8 :: xye = 0.0d0 + eoskey = GRHydro_eoskey +! end EOS Omni vars +#endif + one = 1.0d0 two = 2.0d0 !!$ LEFT !!$ Set required fluid quantities - + invrhol = one / rhol +#if USE_EOS_OMNI + call EOS_Omni_press(eoskey,keytemp,n,& + rhol,epsl,xtemp,xye,pressl,keyerr,anyerr) +! call EOS_Omni_cs2(eoskey,keytemp,n,& +! rhol,epsl,xtemp,xye,cs2l,keyerr,anyerr) + call EOS_Omni_DPressByDEps(eoskey,keytemp,n,& + rhol,epsl,xtemp,xye,dpdepsl,keyerr,anyerr) + call EOS_Omni_DPressByDRho(eoskey,keytemp,n,& + rhol,epsl,xtemp,xye,dpdrhol,keyerr,anyerr) + cs2l = (dpdrhol + pressl * dpdepsl * invrhol**2)/ & + (1.0d0 + epsl + pressl*invrhol) +#else pressl = EOS_Pressure(handle,rhol,epsl) dpdrhol = EOS_DPressByDRho(handle,rhol,epsl) dpdepsl = EOS_DPressByDEps(handle,rhol,epsl) - invrhol = one / rhol cs2l = (dpdrhol + pressl * dpdepsl * invrhol**2)/ & (1.0d0 + epsl + pressl*invrhol) +#endif enthalpyl = one + epsl + pressl * invrhol @@ -161,14 +188,25 @@ !!$ RIGHT !!$ Set required fluid quantities - + invrhor = one / rhor +#if USE_EOS_OMNI + call EOS_Omni_press(eoskey,keytemp,n,& + rhor,epsr,xtemp,xye,pressr,keyerr,anyerr) +! call EOS_Omni_cs2(eoskey,keytemp,n,& +! rhor,epsr,xtemp,xye,cs2r,keyerr,anyerr) + call EOS_Omni_DPressByDEps(eoskey,keytemp,n,& + rhor,epsr,xtemp,xye,dpdepsr,keyerr,anyerr) + call EOS_Omni_DPressByDRho(eoskey,keytemp,n,& + rhor,epsr,xtemp,xye,dpdrhor,keyerr,anyerr) + cs2r = (dpdrhor + pressr * dpdepsr * invrhor**2)/ & + (1.0d0 + epsr + pressr*invrhor) +#else pressr = EOS_Pressure(handle,rhor,epsr) dpdrhor = EOS_DPressByDRho(handle,rhor,epsr) dpdepsr = EOS_DPressByDEps(handle,rhor,epsr) - invrhor = one / rhor cs2r = (dpdrhor + pressr * dpdepsr * invrhor**2)/ & (1.0d0 + epsr + pressr*invrhor) - +#endif enthalpyr = one + epsr + pressr * invrhor vlowxr = gxx*velxr + gxy*velyr + gxz*velzr @@ -247,6 +285,7 @@ kappal = dpdepsl / (dpdepsl - rhol * cs2l) + !!$ Right eigenvector # 1 reivec1l(1) = kappal / (enthalpyl * wl) File [modified]: GRHydro_EoSChangeGamma.F90 Delta lines: +90 -3 =================================================================== --- trunk/src/GRHydro_EoSChangeGamma.F90 2010-07-02 21:51:08 UTC (rev 131) +++ trunk/src/GRHydro_EoSChangeGamma.F90 2010-07-02 21:57:07 UTC (rev 132) @@ -38,7 +38,11 @@ subroutine GRHydro_EoSChangeGamma(CCTK_ARGUMENTS) +#if USE_EOS_OMNI + USE EOS_Omni_Module, only: press_gf, inv_rho_gf, poly_k_cgs, rho_gf +#else USE EOS_Polytrope_Scalars +#endif implicit none @@ -53,16 +57,39 @@ #include "EOS_Base.inc" !!$ Set up the fluid constants +#if USE_EOS_OMNI +! begin EOS Omni vars + integer :: n = 1 + integer :: poly_eoskey = 0 + integer :: keytemp = 0 + integer :: anyerr = 0 + integer :: keyerr(1) = 0 + real*8 :: xpress = 0.0d0 + real*8 :: xeps = 0.0d0 + real*8 :: xtemp = 0.0d0 + real*8 :: xye = 0.0d0 + poly_eoskey = GRHydro_poly_eoskey +! end EOS Omni vars + call EOS_Omni_press(poly_eoskey,keytemp,n,& + 1.0d0,1.0d0,xtemp,xye,xpress,keyerr,anyerr) + + call EOS_Omni_EpsFromPress(poly_eoskey,keytemp,n,& + 1.0d0,1.0d0,xtemp,xye,xpress,xeps,keyerr,anyerr) + local_Gamma = 1.0d0 + xpress/xeps + press = press_gf * poly_k_cgs * & + (rho * inv_rho_gf)**local_Gamma + eps = press / (rho * (local_Gamma - 1.d0)) +#else local_Gamma = 1.d0 + EOS_Pressure(GRHydro_polytrope_handle, 1.d0, 1.d0) / & EOS_SpecificIntEnergy(GRHydro_polytrope_handle, 1.d0, 1.d0) - -!!$ Change the pressure and specific internal energy - press = p_geom_factor * eos_k_cgs * & (rho * rho_geom_factor_inv)**local_Gamma eps = press / (rho * (local_Gamma - 1.d0)) +#endif +!!$ Change the pressure and specific internal energy + !!$ Get the conserved variables. Hardwired polytrope EoS!!! !!$ Note that this call also sets pressure and eps @@ -125,11 +152,33 @@ #include "EOS_Base.inc" !!$ Set up the fluid constants +#if USE_EOS_OMNI +! begin EOS Omni vars + integer :: n = 1 + integer :: poly_eoskey = 0 + integer :: keytemp = 0 + integer :: anyerr = 0 + integer :: keyerr(1) = 0 + real*8 :: xpress = 0.0d0 + real*8 :: xeps = 0.0d0 + real*8 :: xtemp = 0.0d0 + real*8 :: xye = 0.0d0 + poly_eoskey = GRHydro_poly_eoskey +! end EOS Omni vars + call EOS_Omni_press(poly_eoskey,keytemp,n,& + 1.0d0,1.0d0,xtemp,xye,xpress,keyerr,anyerr) + + call EOS_Omni_EpsFromPress(poly_eoskey,keytemp,n,& + 1.0d0,1.0d0,xtemp,xye,xpress,xeps,keyerr,anyerr) + local_Gamma = 1.0d0 + xpress/xeps + local_K = xpress +#else local_Gamma = 1.d0 + EOS_Pressure(GRHydro_polytrope_handle, 1.d0, 1.d0) / & EOS_SpecificIntEnergy(GRHydro_polytrope_handle, 1.d0, 1.d0) local_K = EOS_Pressure(GRHydro_polytrope_handle, 1.d0, 1.d0) +#endif if (abs(local_Gamma - 2.d0) < 1.d-10) then @@ -184,7 +233,11 @@ subroutine GRHydro_EoSChangeGammaK_Shibata(CCTK_ARGUMENTS) +#if USE_EOS_OMNI + USE EOS_Omni_Module, only: press_gf, inv_rho_gf, poly_k_cgs, rho_gf +#else USE EOS_Polytrope_Scalars +#endif implicit none @@ -203,7 +256,39 @@ #include "EOS_Base.inc" !!$ Set up the fluid constants +#if USE_EOS_OMNI +! begin EOS Omni vars + integer :: n = 1 + integer :: poly_eoskey = 0 + integer :: keytemp = 0 + integer :: anyerr = 0 + integer :: keyerr(1) = 0 + real*8 :: xpress = 0.0d0 + real*8 :: xeps = 0.0d0 + real*8 :: xtemp = 0.0d0 + real*8 :: xye = 0.0d0 + poly_eoskey = GRHydro_poly_eoskey +! end EOS Omni vars + call EOS_Omni_press(poly_eoskey,keytemp,n,& + 1.0d0,1.0d0,xtemp,xye,xpress,keyerr,anyerr) + + call EOS_Omni_EpsFromPress(poly_eoskey,keytemp,n,& + 1.0d0,1.0d0,xtemp,xye,xpress,xeps,keyerr,anyerr) + local_Gamma = 1.0d0 + xpress/xeps + local_K = xpress + + eos_k_initial_cgs = initial_k * rho_gf**initial_Gamma / press_gf + + press = (local_Gamma - 1.d0) / (initial_Gamma - 1.0d0 ) * press_gf * eos_k_initial_cgs * & + (rho * rho_gf) ** initial_Gamma + + eps = press_gf * eos_k_initial_cgs * & + (rho * inv_rho_gf) ** initial_Gamma / & + (rho * (initial_Gamma - 1.0d0)) + +#else + call CCTK_INFO("Pulling the rug...by changing K and Gamma") local_Gamma = 1.d0 + EOS_Pressure(GRHydro_polytrope_handle, 1.d0, 1.d0) / & @@ -220,6 +305,8 @@ (rho * rho_geom_factor_inv) ** initial_Gamma / & (rho * (initial_Gamma - 1.0d0)) +#endif + do k = 1, cctk_lsh(3) do j = 1, cctk_lsh(2) do i = 1, cctk_lsh(1) File [modified]: GRHydro_Marquina.F90 Delta lines: +1 -7 =================================================================== --- trunk/src/GRHydro_Marquina.F90 2010-07-02 21:51:08 UTC (rev 131) +++ trunk/src/GRHydro_Marquina.F90 2010-07-02 21:57:07 UTC (rev 132) @@ -37,10 +37,6 @@ implicit none -#ifdef _EOS_BASE_INC_ -#undef _EOS_BASE_INC_ -#endif -#include "EOS_Base.inc" DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS @@ -220,7 +216,6 @@ endif !!$ END: Check for NaN value (1st check) -!!$ pressp = EOS_Pressure(GRHydro_eos_handle,primp(1),primp(5)) !!$right state @@ -244,7 +239,6 @@ endif !!$ END: Check for NaN value (2nd check) -!!$ pressm_i = EOS_Pressure(GRHydro_eos_handle,primm_i(1),primm_i(5)) !!$eigenvalues and right eigenvectors @@ -353,7 +347,7 @@ syflux(i,j,k) = f_marquina(3) szflux(i,j,k) = f_marquina(4) tauflux(i,j,k) = f_marquina(5) - + enddo enddo enddo File [modified]: GRHydro_Prim2Con.F90 Delta lines: +46 -3 =================================================================== --- trunk/src/GRHydro_Prim2Con.F90 2010-07-02 21:51:08 UTC (rev 131) +++ trunk/src/GRHydro_Prim2Con.F90 2010-07-02 21:57:07 UTC (rev 132) @@ -132,6 +132,21 @@ #include "EOS_Base.inc" +#if USE_EOS_OMNI +! begin EOS Omni vars + integer :: n = 1 + integer :: eoskey = 0 + integer :: keytemp = 0 + integer :: anyerr = 0 + integer :: keyerr(1) = 0 + real*8 :: xpress = 0.0d0 + real*8 :: xeps = 0.0d0 + real*8 :: xtemp = 0.0d0 + real*8 :: xye = 0.0d0 + eoskey = GRHydro_eoskey +! end EOS Omni vars +#endif + w = 1.d0 / sqrt(1.d0 - (gxx*dvelx*dvelx + gyy*dvely*dvely + gzz & *dvelz*dvelz + 2*gxy*dvelx*dvely + 2*gxz*dvelx *dvelz + 2*gyz& *dvely*dvelz)) @@ -145,7 +160,12 @@ endif !!$ END: Check for NaN value +#if USE_EOS_OMNI + call EOS_Omni_press(eoskey,keytemp,n,& + drho,deps,xtemp,xye,dpress,keyerr,anyerr) +#else dpress = EOS_Pressure(handle, drho, deps) +#endif vlowx = gxx*dvelx + gxy*dvely + gxz*dvelz vlowy = gxy*dvelx + gyy*dvely + gyz*dvelz @@ -339,6 +359,21 @@ #undef _EOS_BASE_INC_ #endif #include "EOS_Base.inc" + +#if USE_EOS_OMNI +! begin EOS Omni vars + integer :: n = 1 + integer :: poly_eoskey = 0 + integer :: keytemp = 0 + integer :: anyerr = 0 + integer :: keyerr(1) = 0 + real*8 :: xpress = 0.0d0 + real*8 :: xeps = 0.0d0 + real*8 :: xtemp = 0.0d0 + real*8 :: xye = 0.0d0 + poly_eoskey = GRHydro_poly_eoskey +! end EOS Omni vars +#endif w_tmp = gxx*dvelx*dvelx + gyy*dvely*dvely + gzz *dvelz*dvelz + & 2*gxy*dvelx*dvely + 2*gxz*dvelx*dvelz + 2*gyz*dvely*dvelz @@ -357,11 +392,19 @@ else w = 1.d0 / sqrt(1.d0 - w_tmp) endif - + +#if USE_EOS_OMNI + call EOS_Omni_press(poly_eoskey,keytemp,n,& + drho,xeps,xtemp,xye,dpress,keyerr,anyerr) + + call EOS_Omni_EpsFromPress(poly_eoskey,keytemp,n,& + drho,xeps,xtemp,xye,dpress,deps,keyerr,anyerr) +#else if (handle .ge. 0) then - dpress = EOS_Pressure(handle, drho, deps) - deps = EOS_SpecificIntEnergy(handle, drho, dpress) + dpress = EOS_Pressure(handle, drho, deps) + deps = EOS_SpecificIntEnergy(handle, drho, dpress) end if +#endif vlowx = gxx*dvelx + gxy*dvely + gxz*dvelz vlowy = gxy*dvelx + gyy*dvely + gyz*dvelz

1 0

[svn:einsteintoolkit] GRHydro/trunk/ (Rev. 131)
by cott＠tapir.caltech.edu 02 Jul '10

02 Jul '10

User: cott Date: 2010/07/02 04:51 PM Modified: /trunk/ interface.ccl Log: * hooks for EOS omni functions. File Changes: Directory: /trunk/ ================== File [modified]: interface.ccl Delta lines: +78 -0 =================================================================== --- trunk/interface.ccl 2010-07-02 21:50:18 UTC (rev 130) +++ trunk/interface.ccl 2010-07-02 21:51:08 UTC (rev 131) @@ -157,7 +157,84 @@ CCTK_STRING IN cactus_group_name) USES FUNCTION GZPatchSystem_register_cxform + ####################################################### +### Functions provided by the Omni EOS interface ### +####################################################### + +void FUNCTION EOS_Omni_press(CCTK_INT IN eoskey, \ + CCTK_INT IN havetemp, \ + CCTK_INT IN npoints, \ + CCTK_REAL IN ARRAY rho, \ + CCTK_REAL INOUT ARRAY eps, \ + CCTK_REAL INOUT ARRAY temp, \ + CCTK_REAL IN ARRAY ye, \ + CCTK_REAL OUT ARRAY press, \ + CCTK_INT OUT ARRAY keyerr, \ + CCTK_INT OUT anyerr) + +USES FUNCTION EOS_Omni_press + + +void FUNCTION EOS_Omni_DPressByDEps(CCTK_INT IN eoskey, \ + CCTK_INT IN havetemp, \ + CCTK_INT IN npoints, \ + CCTK_REAL IN ARRAY rho, \ + CCTK_REAL INOUT ARRAY eps, \ + CCTK_REAL INOUT ARRAY temp, \ + CCTK_REAL IN ARRAY ye, \ + CCTK_REAL OUT ARRAY dpdepsrho, \ + CCTK_INT OUT ARRAY keyerr, \ + CCTK_INT OUT anyerr) + +USES FUNCTION EOS_Omni_DPressByDEps + + +void FUNCTION EOS_Omni_DPressByDRho(CCTK_INT IN eoskey, \ + CCTK_INT IN havetemp, \ + CCTK_INT IN npoints, \ + CCTK_REAL IN ARRAY rho, \ + CCTK_REAL INOUT ARRAY eps, \ + CCTK_REAL INOUT ARRAY temp, \ + CCTK_REAL IN ARRAY ye, \ + CCTK_REAL OUT ARRAY dpdrhoeps, \ + CCTK_INT OUT ARRAY keyerr, \ + CCTK_INT OUT anyerr) + +USES FUNCTION EOS_Omni_DPressByDRho + + +void FUNCTION EOS_Omni_cs2(CCTK_INT IN eoskey, \ + CCTK_INT IN havetemp, \ + CCTK_INT IN npoints, \ + CCTK_REAL IN ARRAY rho, \ + CCTK_REAL INOUT ARRAY eps, \ + CCTK_REAL INOUT ARRAY temp, \ + CCTK_REAL IN ARRAY ye, \ + CCTK_REAL OUT ARRAY cs2, \ + CCTK_INT OUT ARRAY keyerr, \ + CCTK_INT OUT anyerr) + +USES FUNCTION EOS_Omni_cs2 + +void FUNCTION EOS_Omni_EpsFromPress(CCTK_INT IN eoskey, \ + CCTK_INT IN havetemp, \ + CCTK_INT IN npoints, \ + CCTK_REAL IN ARRAY rho, \ + CCTK_REAL INOUT ARRAY eps, \ + CCTK_REAL INOUT ARRAY temp, \ + CCTK_REAL IN ARRAY ye, \ + CCTK_REAL IN ARRAY press, \ + CCTK_REAL OUT ARRAY xeps, \ + CCTK_INT OUT ARRAY keyerr, \ + CCTK_INT OUT anyerr) + +USES FUNCTION EOS_Omni_EpsFromPress + + + + +####################################################### ### Functions provided by the general EOS interface ### ####################################################### @@ -378,6 +455,7 @@ GRHydro_mppm_xwind } "debug variable for flux eigenvalues in mppm" + CCTK_REAL particles TYPE=ARRAY DIM=1 TIMELEVELS=3 SIZE=number_of_particles GHOSTSIZE=0 DISTRIB=DEFAULT { particle_x

1 0

[svn:einsteintoolkit] GRHydro/trunk/ (Rev. 130)
by cott＠tapir.caltech.edu 02 Jul '10

02 Jul '10

User: cott Date: 2010/07/02 04:50 PM Modified: /trunk/ param.ccl Log: * add two parameters to get EOS Omni working. File Changes: Directory: /trunk/ ================== File [modified]: param.ccl Delta lines: +17 -0 =================================================================== --- trunk/param.ccl 2010-06-30 06:49:09 UTC (rev 129) +++ trunk/param.ccl 2010-07-02 21:50:18 UTC (rev 130) @@ -242,7 +242,24 @@ .* :: "Can be anything" } "Ideal_Fluid" +####################### EOS Omni Parameters ############################## +int GRHydro_poly_eoskey "eoskey for poly EOS" +{ + 1:* :: "Must be larger than 0" +} 1 + + +int GRHydro_eoskey "eoskey determining which EOS to use" +{ + 1 :: "Polytropic EOS" + 2 :: "Gamma-Law EOS" + 3 :: "Hybrid EOS" + 4:* :: "something else" +} 2 + +####################### Other Parameters ############################## + #Parameters controlling conservative <-> primitive change. int GRHydro_c2p_warnlevel "Warnlevel for Con2Prim warnings" STEERABLE=ALWAYS

1 0

[svn:einsteintoolkit] incoming/EOS_Omni/ (Rev. 3)
by cott＠tapir.caltech.edu 02 Jul '10

02 Jul '10

User: cott Date: 2010/07/02 04:41 PM Modified: /EOS_Omni/ param.ccl /EOS_Omni/src/ EOS_Omni_Module.F90, EOS_Omni_SingleVarCalls.F90, EOS_Omni_Startup.F90 Log: * add hybrid EOS File Changes: Directory: /EOS_Omni/ ===================== File [modified]: param.ccl Delta lines: +33 -2 =================================================================== --- EOS_Omni/param.ccl 2010-07-02 21:41:36 UTC (rev 2) +++ EOS_Omni/param.ccl 2010-07-02 21:41:51 UTC (rev 3) @@ -4,18 +4,23 @@ restricted: # poly EOS -CCTK_REAL poly_gamma "Adiabatic Index for poly EOS" +REAL poly_gamma "Adiabatic Index for poly EOS" { : :: "" } 2.0 +REAL poly_gamma_ini "Initial Adiabatic Index for poly EOS" +{ + : :: "" +} 2.0 + REAL poly_k "Polytropic constant in c=G=Msun=1" { : :: "" } 100.0 # gamma-law EOS -CCTK_REAL gl_gamma "Adiabatic Index for gamma-law EOS" +REAL gl_gamma "Adiabatic Index for gamma-law EOS" { : :: "" } 2.0 @@ -24,3 +29,29 @@ { : :: "" } 100.0 + +# hybrid EOS +REAL hybrid_gamma1 "subnuclear adiabatic Index for hybrid EOS" +{ + : :: "" +} 1.325 + +REAL hybrid_gamma2 "subnuclear adiabatic Index for hybrid EOS" +{ + : :: "" +} 2.5 + +REAL hybrid_gamma_th "Thermal gamma for hybrid EOS" +{ + : :: "" +} 1.5 + +REAL hybrid_k1 "Polytropic constant in c=G=Msun=1 for hybrid EOS" +{ + : :: "" +} 0.4640517 + +REAL hybrid_rho_nuc "Density at which to switch betwen gammas; c=G=Msun=1" +{ + : :: "" +} 3.238607e-4 Directory: /EOS_Omni/src/ ========================= File [modified]: EOS_Omni_Module.F90 Delta lines: +11 -4 =================================================================== --- EOS_Omni/src/EOS_Omni_Module.F90 2010-07-02 21:41:36 UTC (rev 2) +++ EOS_Omni/src/EOS_Omni_Module.F90 2010-07-02 21:41:51 UTC (rev 3) @@ -2,15 +2,19 @@ implicit none - real*8,parameter :: rho_gf = 1.61930347d-18 - real*8,parameter :: press_gf = 1.80171810d-39 + real*8,parameter :: rho_gf = 1.61620075314614d-18 +! real*8,parameter :: rho_gf = 1.61930347d-18 + real*8,parameter :: press_gf = 1.7982953469278d-39 +! real*8,parameter :: press_gf = 1.80171810d-39 real*8,parameter :: eps_gf = 1.11265006d-21 real*8,parameter :: time_gf = 2.03001708d+05 real*8,parameter :: mass_gf = 5.02765209d-34 real*8,parameter :: length_gf = 6.77140812d-06 - real*8,parameter :: inv_rho_gf = 6.17549470205236d17 - real*8,parameter :: inv_press_gf = 5.55025783445257d38 + real*8,parameter :: inv_rho_gf = 6.18735016707159d17 +! real*8,parameter :: inv_rho_gf = 6.17549470205236d17 + real*8,parameter :: inv_press_gf = 5.56082181777535d38 +! real*8,parameter :: inv_press_gf = 5.55025783445257d38 real*8,parameter :: inv_eps_gf = 8.98755175549085d20 real*8,parameter :: inv_time_gf = 4.92606692747629d-6 real*8,parameter :: inv_mass_gf = 1.98899999860571d33 @@ -19,4 +23,7 @@ real*8 :: poly_k_cgs = 0.0d0 real*8 :: gl_k_cgs = 0.0d0 + real*8 :: hybrid_k1_cgs = 0.0d0 + real*8 :: hybrid_k2_cgs = 0.0d0 + end module EOS_Omni_Module File [modified]: EOS_Omni_SingleVarCalls.F90 Delta lines: +106 -13 =================================================================== --- EOS_Omni/src/EOS_Omni_SingleVarCalls.F90 2010-07-02 21:41:36 UTC (rev 2) +++ EOS_Omni/src/EOS_Omni_SingleVarCalls.F90 2010-07-02 21:41:51 UTC (rev 3) @@ -24,8 +24,11 @@ CCTK_REAL, intent(out) :: press(npoints) ! local vars - integer :: i - character(256) :: warnstring + integer :: i + character(256) :: warnstring + real*8 :: hybrid_local_gamma, hybrid_local_k_cgs, & + hybrid_p_poly, hybrid_p_th + real*8,parameter :: zero = 0.0d0 anyerr = 0 keyerr(:) = 0 @@ -57,6 +60,29 @@ press(i) = (gl_gamma - 1.0d0) * rho(i) * eps(i) enddo + case (3) + ! hybrid EOS + do i=1,npoints + if(rho(i).gt.hybrid_rho_nuc) then + hybrid_local_gamma = hybrid_gamma2 + hybrid_local_k_cgs = hybrid_k2_cgs + else + hybrid_local_gamma = hybrid_gamma1 + hybrid_local_k_cgs = hybrid_k1_cgs + endif + hybrid_p_poly = press_gf * hybrid_local_k_cgs * & + (rho(i) * inv_rho_gf)**hybrid_local_gamma + hybrid_p_th = - press_gf * hybrid_local_k_cgs * (hybrid_gamma_th - 1.d0) / & + (hybrid_local_gamma - 1.0d0) * (rho(i) * inv_rho_gf)**hybrid_local_gamma + & + (hybrid_gamma_th - 1.0d0) * rho(i) * eps(i) - & + (hybrid_gamma_th - 1.d0) * (hybrid_local_gamma - hybrid_gamma1) / & + (hybrid_gamma1 - 1.d0) / (hybrid_gamma2 - 1.d0) * & + press_gf * hybrid_k1_cgs * inv_rho_gf**hybrid_gamma1 * & + hybrid_rho_nuc**(hybrid_gamma1 - 1.d0) * rho(i) + hybrid_p_th = max(zero, hybrid_p_th) + press(i) = hybrid_p_poly + hybrid_p_th + enddo + case DEFAULT write(warnstring,*) "eoskey ",eoskey," not implemented!" call CCTK_WARN(0,warnstring) @@ -81,8 +107,11 @@ CCTK_REAL, intent(out) :: dpdrhoe(npoints) ! local vars - integer :: i - character(256) :: warnstring + integer :: i + character(256) :: warnstring + real*8 :: hybrid_local_gamma, hybrid_local_k_cgs, & + hybrid_dp_poly, hybrid_dp_th1, hybrid_dp_th2 + real*8,parameter :: zero = 0.0d0 anyerr = 0 keyerr(:) = 0 @@ -115,7 +144,34 @@ dpdrhoe(i) = (gl_gamma-1.0d0) * & eps(i) enddo + case (3) + ! hybrid EOS + do i=1,npoints + if(rho(i).gt.hybrid_rho_nuc) then + hybrid_local_gamma = hybrid_gamma2 + hybrid_local_k_cgs = hybrid_k2_cgs + else + hybrid_local_gamma = hybrid_gamma1 + hybrid_local_k_cgs = hybrid_k1_cgs + endif + hybrid_dp_poly = hybrid_local_gamma * press_gf * & + hybrid_local_k_cgs * rho(i)**(hybrid_local_gamma - 1.0d0) * & + inv_rho_gf**hybrid_local_gamma + hybrid_dp_th1 = - hybrid_local_gamma * press_gf * hybrid_local_k_cgs * & + (hybrid_gamma_th - 1.d0) / (hybrid_local_gamma - 1.d0) * & + rho(i)**(hybrid_local_gamma - 1.d0) * inv_rho_gf**hybrid_local_gamma + + hybrid_dp_th2 = (hybrid_gamma_th - 1.d0) * eps(i) & + - (hybrid_gamma_th - 1.d0) * (hybrid_local_gamma - hybrid_gamma1) / & + (hybrid_gamma1 - 1.d0) / (hybrid_gamma2 - 1.d0) * & + press_gf * hybrid_k1_cgs * inv_rho_gf**hybrid_gamma1 * & + hybrid_rho_nuc**(hybrid_gamma1 - 1.d0) + + dpdrhoe(i) = hybrid_dp_poly + hybrid_dp_th1 + hybrid_dp_th2 + + enddo + case DEFAULT write(warnstring,*) "eoskey ",eoskey," not implemented!" call CCTK_WARN(0,warnstring) @@ -138,9 +194,8 @@ CCTK_REAL, intent(out) :: dpdepsrho(npoints) ! local vars - integer :: i - character(256) :: warnstring - + integer :: i + character(256) :: warnstring anyerr = 0 keyerr(:) = 0 @@ -170,6 +225,11 @@ dpdepsrho(i) = (gl_gamma - 1.0d0) * & rho(i) enddo + case (3) + ! hybrid EOS + do i=1,npoints + dpdepsrho(i) = (hybrid_gamma_th - 1.0d0) * rho(i) + enddo case DEFAULT write(warnstring,*) "eoskey ",eoskey," not implemented!" @@ -193,9 +253,12 @@ CCTK_REAL, intent(out) :: cs2(npoints) ! local vars - integer :: i - character(256) :: warnstring - real*8 :: xpress,xdpdrhoe,xdpderho + integer :: i + character(256) :: warnstring + real*8 :: xpress,xdpdrhoe,xdpderho + real*8 :: hybrid_local_gamma, hybrid_local_k_cgs, & + hybrid_p_poly, hybrid_p_th + real*8,parameter :: zero = 0.0d0 anyerr = 0 keyerr(:) = 0 @@ -232,7 +295,31 @@ cs2(i) = (xdpdrhoe + xpress * xdpderho / (rho(i)**2)) / & (1.0d0 + eps(i) + xpress/rho(i)) enddo - + case(3) + ! hybrid EOS + do i=1,npoints + if(rho(i).gt.hybrid_rho_nuc) then + hybrid_local_gamma = hybrid_gamma2 + hybrid_local_k_cgs = hybrid_k2_cgs + else + hybrid_local_gamma = hybrid_gamma1 + hybrid_local_k_cgs = hybrid_k1_cgs + endif + ! first calculate the pressure + hybrid_p_poly = press_gf * hybrid_local_k_cgs * & + (rho(i) * inv_rho_gf)**hybrid_local_gamma + hybrid_p_th = - press_gf * hybrid_local_k_cgs * (hybrid_gamma_th - 1.d0) / & + (hybrid_local_gamma - 1.0d0) * (rho(i) * inv_rho_gf)**hybrid_local_gamma + & + (hybrid_gamma_th - 1.0d0) * rho(i) * eps(i) - & + (hybrid_gamma_th - 1.d0) * (hybrid_local_gamma - hybrid_gamma1) / & + (hybrid_gamma1 - 1.d0) / (hybrid_gamma2 - 1.d0) * & + press_gf * hybrid_k1_cgs * inv_rho_gf**hybrid_gamma1 * & + hybrid_rho_nuc**(hybrid_gamma1 - 1.d0) * rho(i) + hybrid_p_th = max(zero, hybrid_p_th) + xpress = hybrid_p_poly + hybrid_p_th + cs2(i) = (hybrid_local_gamma * hybrid_p_poly + hybrid_gamma_th * hybrid_p_th) / & + rho(i) / (1.0d0 + eps(i) + xpress/rho(i)) + enddo case DEFAULT write(warnstring,*) "eoskey ",eoskey," not implemented!" call CCTK_WARN(0,warnstring) @@ -272,14 +359,16 @@ ! polytropic EOS do i=1,npoints xeps(i) = press(i) / (poly_gamma - 1.0d0) / rho(i) - eps(i) = xeps(i) enddo case (2) ! gamma-law EOS do i=1,npoints xeps(i) = press(i) / (gl_gamma - 1.0d0) / rho(i) - eps(i) = xeps(i) enddo + case (3) + ! hybrid EOS + write(warnstring,*) "EOS_Omni_EpsFromPress call not supported for hybrid EOS" + call CCTK_WARN(0,warnstring) case DEFAULT write(warnstring,*) "eoskey ",eoskey," not implemented!" call CCTK_WARN(0,warnstring) @@ -325,6 +414,10 @@ do i=1,npoints rho(i) = press(i) / ((gl_gamma - 1.0d0)*eps(i)) enddo + case (3) + ! hybrid EOS + write(warnstring,*) "EOS_Omni_RestMassDensityFromEpsPress not supported for hybrid EOS" + call CCTK_WARN(0,warnstring) case DEFAULT write(warnstring,*) "eoskey ",eoskey," not implemented!" call CCTK_WARN(0,warnstring) File [modified]: EOS_Omni_Startup.F90 Delta lines: +7 -2 =================================================================== --- EOS_Omni/src/EOS_Omni_Startup.F90 2010-07-02 21:41:36 UTC (rev 2) +++ EOS_Omni/src/EOS_Omni_Startup.F90 2010-07-02 21:41:51 UTC (rev 3) @@ -11,8 +11,13 @@ DECLARE_CCTK_PARAMETERS DECLARE_CCTK_ARGUMENTS - poly_k_cgs = poly_k * rho_gf**poly_gamma / press_gf + poly_k_cgs = poly_k * rho_gf**poly_gamma_ini / press_gf - gl_k_cgs = gl_k * rho_gf**gl_gamma / press_gf + gl_k_cgs = gl_k * rho_gf**poly_gamma_ini / press_gf + hybrid_k1_cgs = hybrid_k1 * rho_gf**poly_gamma_ini / press_gf + + hybrid_k2_cgs = hybrid_k1 * (hybrid_rho_nuc * inv_rho_gf)**(hybrid_gamma1-hybrid_gamma2) + + end subroutine EOS_Omni_Startup

1 0

[svn:einsteintoolkit] incoming/EOS_Omni/doc/ (Rev. 2)
by cott＠tapir.caltech.edu 02 Jul '10

02 Jul '10

User: cott Date: 2010/07/02 04:41 PM Modified: /EOS_Omni/doc/ documentation.tex Log: * add prelim documentation File Changes: Directory: /EOS_Omni/doc/ ========================= File [modified]: documentation.tex Delta lines: +69 -14 =================================================================== --- EOS_Omni/doc/documentation.tex 2010-06-29 14:59:02 UTC (rev 1) +++ EOS_Omni/doc/documentation.tex 2010-07-02 21:41:36 UTC (rev 2) @@ -99,7 +99,9 @@ % Add an abstract for this thorn's documentation \begin{abstract} - +This Thorn provides a unified EOS interface and implements +multiple EOS. These currently are the polytropic EOS, the +gamma-law EOS, and the hybrid EOS. \end{abstract} % The following sections are suggestive only. @@ -107,36 +109,89 @@ \section{Introduction} -\section{Physical System} +Equations of State (EOS) are crucial for hydrodynamics and hydro +codes (as well as other codes needing/providing microphysics) +are closely coupled to EOS and call EOS routines many times for +each grid point during the caculation of a time update. -\section{Numerical Implementation} +{\tt EOS\_Omni} is presently coded for cold and hot EOS, +including those based on microphysical models. It does currently +assume nuclear statistical equilibrium (NSE) with density, +specific internal energy (or temperature), and electron fraction +being the independent variables. {\tt EOS\_Omni} may be called +on vectors or on single points. \section{Using This Thorn} -\subsection{Obtaining This Thorn} \subsection{Basic Usage} -\subsection{Special Behaviour} +{\tt EOS\_Omni} works via the aliased-function interface and +EOS functions to be used must be declared in {\tt interface.ccl}. +Here is an example call/{interface.ccl} entry: +\begin{verbatim} +void FUNCTION EOS_Omni_press(CCTK_INT IN eoskey, \ + CCTK_INT IN havetemp, \ + CCTK_INT IN npoints, \ + CCTK_REAL IN ARRAY rho, \ + CCTK_REAL INOUT ARRAY eps, \ + CCTK_REAL INOUT ARRAY temp, \ + CCTK_REAL IN ARRAY ye, \ + CCTK_REAL OUT ARRAY press, \ + CCTK_INT OUT ARRAY keyerr, \ + CCTK_INT OUT anyerr) +\end{verbatim} -\subsection{Interaction With Other Thorns} +Here, +\begin{itemize} +\item {\tt eoskey} is the type of EOS to be used in this call. + \begin{itemize} + \item {\tt eoskey = 1}: Polytropic EOS + \item {\tt eoskey = 2}: Gamma-Law EOS + \item {\tt eoskey = 3}: Hybrid EOS (2 Polys, 1 Gamma-Law), used + for core-collapse simulations. + \end{itemize} +\item {\tt havetemp} tells the EOS if the call is $X=X(\rho,T,Y_e)$ + ({\tt havetemp = 1}) or $X=X(\rho,\epsilon,Y_e)$ ({\tt havetemp = 0}), +\item {\tt npoints} tells the EOS how many data points are passed + in, +\item {\tt rho,eps,temp,ye,press} are obvious, +\item {\tt keyerr} is an array (with n entries for n data points) + containing error codes (relevant only for tabular EOS), +\item {\tt anyerr} is an integer that is .gt. 0 in the case + any error occured. +\end{itemize} -\subsection{Examples} +\subsection{Parameter Settings} -\subsection{Support and Feedback} +Many hydro codes require a fallback EOS in case something goes +wrong. This is also true for the Einstein Toolkit GR hydro code {\tt +EinsteinEvolve/GRHydro}. If you want to use {\tt EOS\_Omni} with {\tt +GRHydro}, you must set parameters for the EOS of your choice and, {\bf +in addition}, the following parameters must be set to sensible values: -\section{History} +\begin{verbatim} +eos_omni::poly_gamma +eos_omni::poly_gamma_ini +eos_omni::poly_k +\end{verbatim} -\subsection{Thorn Source Code} +The only non-obvious parameter here is {\tt poly\_gamma\_ini}. +In most simulations it should be set equal to {\tt poly\_gamma}. +In simulations that are run with a different adiabatic index than +what was used to set up the initial data, {\tt poly\_gamma} should +be the evolution value, and {\tt poly\_gamma\_ini} should be the initial +data value. -\subsection{Thorn Documentation} +Check out {\tt param.ccl} for parameters for other EOS. -\subsection{Acknowledgements} +\subsection{Special Behavior} +If special behavior occurs, please report it to {\tt cott(a)tapir.caltech.edu}. -\begin{thebibliography}{9} +%\begin{thebibliography}{9} -\end{thebibliography} +%\end{thebibliography} % Do not delete next line % END CACTUS THORNGUIDE

1 0

[svn:einsteintoolkit] manifest/trunk/ (Rev. 32)
by schnetter＠cct.lsu.edu 30 Jun '10

30 Jun '10

User: eschnett Date: 2010/06/30 09:55 PM Modified: /trunk/ einsteintoolkit.th Log: Add CactusNumerical/InterpToArray File Changes: Directory: /trunk/ ================== File [modified]: einsteintoolkit.th Delta lines: +1 -0 =================================================================== --- trunk/einsteintoolkit.th 2010-06-18 15:50:44 UTC (rev 31) +++ trunk/einsteintoolkit.th 2010-07-01 02:55:19 UTC (rev 32) @@ -82,6 +82,7 @@ CactusNumerical/MoL CactusNumerical/Cartoon2D CactusNumerical/Dissipation +CactusNumerical/InterpToArray CactusNumerical/Noise CactusNumerical/Norms CactusNumerical/Periodic

1 0

[svn:einsteintoolkit] IDConstraintViolate/trunk/ (Rev. 14)
by schnetter＠cct.lsu.edu 30 Jun '10

30 Jun '10

User: eschnett Date: 2010/06/30 08:07 PM Modified: /trunk/ interface.ccl, param.ccl /trunk/par/ Gaussian.par Log: Update the example parameter file in thorn IDConstraintViolate to using McLachlan. (Note: the test case should also be updated.) File Changes: Directory: /trunk/par/ ====================== File [modified]: Gaussian.par Delta lines: +239 -23 =================================================================== --- trunk/par/Gaussian.par 2010-06-17 19:59:27 UTC (rev 13) +++ trunk/par/Gaussian.par 2010-07-01 01:07:10 UTC (rev 14) @@ -1,35 +1,251 @@ +Cactus::cctk_run_title = "Constraint violating Gaussian initial data" -!DESC "Constraint violating Gaussian initial data" +Cactus::cctk_full_warnings = yes +Cactus::highlight_warning_messages = no -ActiveThorns = "SpaceMask Boundary ADMBase ADM ADMMacros ADMAnalysis ADMConstraints CoordGauge StaticConformal CartGrid3D IOBasic IOUtil IOASCII IDConstraintViolate PsiKadelia Pugh PughReduce PughSlab Time" - -#------------------------------------------------------------------------------ +Cactus::terminate = "time" +Cactus::cctk_final_time = 10.0 -driver::global_nsize = 40 -cactus::cctk_itlast = 1000 -IOBasic::outInfo_vars = "admbase::alp admbase::gxx admconstraints::ham" -IOBasic::outInfo_every = 5 -IOASCII::out1D_vars = "admbase::alp admbase::gxx admanalysis::grr psikadelia::psi0re admconstraints::ham admconstraints::momentum" -IOASCII::out1D_every = 5 +ActiveThorns = "IOUtil" -#------------------------------------------------------------------------------ +IO::out_dir = $parfile -IO::out_dir = "Gaussian" -grid::type = "box" -time::dtfac = 0.25 - -admbase::metric_type = "physical" -admbase::initial_data = "constraint violating gaussian" -admbase::evolution_method = "ADM" -admbase::lapse_evolution_method = "1+log" -admbase::initial_lapse = "one" +ActiveThorns = "Fortran" -adm::bound = "radiative" -adm::method = "stagleap" +ActiveThorns = "GSL" -#------------------------------------------------------------------------------ +ActiveThorns = "GenericFD" +ActiveThorns = "LoopControl" + + + +ActiveThorns = "InitBase" + + + +ActiveThorns = "Carpet CarpetLib CarpetReduce" + +Carpet::domain_from_coordbase = yes + +driver::ghost_size = 3 + +Carpet::init_fill_timelevels = yes + +Carpet::poison_new_timelevels = yes +CarpetLib::poison_new_memory = yes + + + +ActiveThorns = "NaNChecker" + +NaNChecker::check_every = 512 +NaNChecker::action_if_found = "just warn" +NaNChecker::check_vars = " + ML_BSSN::ML_log_confac + ML_BSSN::ML_metric + ML_BSSN::ML_trace_curv + ML_BSSN::ML_curv + ML_BSSN::ML_Gamma + ML_BSSN::ML_lapse + ML_BSSN::ML_shift + ML_BSSN::ML_dtlapse + ML_BSSN::ML_dtshift + ADMBase::metric + ADMBase::curv + ADMBase::lapse + ADMBase::shift + ADMBase::dtlapse + ADMBase::dtshift +" + + + +ActiveThorns = "Boundary CartGrid3D CoordBase SymBase" + +CoordBase::domainsize = "minmax" + +CoordBase::xmin = -1.0 +CoordBase::ymin = -1.0 +CoordBase::zmin = -1.0 +CoordBase::xmax = +1.0 +CoordBase::ymax = +1.0 +CoordBase::zmax = +1.0 +CoordBase::dx = 0.05 +CoordBase::dy = 0.05 +CoordBase::dz = 0.05 + +CoordBase::boundary_size_x_lower = 3 +CoordBase::boundary_size_y_lower = 3 +CoordBase::boundary_size_z_lower = 3 +CoordBase::boundary_size_x_upper = 3 +CoordBase::boundary_size_y_upper = 3 +CoordBase::boundary_size_z_upper = 3 + +GenericFD::boundary_width = 3 + +CartGrid3D::type = "coordbase" + + + +ActiveThorns = "SphericalSurface" + + + +ActiveThorns = "MoL Time" + +MoL::ODE_Method = "RK4" +MoL::MoL_Intermediate_Steps = 4 +MoL::MoL_Num_Scratch_Levels = 1 + +Time::dtfac = 0.25 + + + +ActiveThorns = "ADMBase ADMCoupling ADMMacros CoordGauge SpaceMask StaticConformal TmunuBase" + + + +ActiveThorns = "IDConstraintViolate" + +ADMBase::initial_data = "constraint violating gaussian" +ADMBase::initial_lapse = "one" +ADMBase::initial_shift = "zero" +ADMBase::initial_dtlapse = "zero" +ADMBase::initial_dtshift = "zero" + + + +ActiveThorns = "ML_BSSN ML_BSSN_Helper NewRad" + +ADMBase::evolution_method = "ML_BSSN" +ADMBase::lapse_evolution_method = "ML_BSSN" +ADMBase::shift_evolution_method = "ML_BSSN" +ADMBase::dtlapse_evolution_method = "ML_BSSN" +ADMBase::dtshift_evolution_method = "ML_BSSN" + +ML_BSSN::timelevels = 3 + +ML_BSSN::harmonicF = 2.0 # 1+log +ML_BSSN::harmonicN = 1 # 1+log +ML_BSSN::ShiftGammaCoeff = 0.75 +ML_BSSN::LapseACoeff = 1.0 +ML_BSSN::ShiftBCoeff = 1.0 +ML_BSSN::AlphaDriver = 1.0 +ML_BSSN::BetaDriver = 1.0 +ML_BSSN::LapseAdvectionCoeff = 0.0 +ML_BSSN::ShiftAdvectionCoeff = 0.0 + +ML_BSSN::MinimumLapse = 1.0e-8 + +ML_BSSN::dt_lapse_shift_method = "noLapseShiftAdvection" + +ML_BSSN::my_initial_boundary_condition = "extrapolate-gammas" +ML_BSSN::my_rhs_boundary_condition = "NewRad" +Boundary::radpower = 2 + +ML_BSSN::ML_log_confac_bound = "none" +ML_BSSN::ML_metric_bound = "none" +ML_BSSN::ML_Gamma_bound = "none" +ML_BSSN::ML_trace_curv_bound = "none" +ML_BSSN::ML_curv_bound = "none" +ML_BSSN::ML_lapse_bound = "none" +ML_BSSN::ML_dtlapse_bound = "none" +ML_BSSN::ML_shift_bound = "none" +ML_BSSN::ML_dtshift_bound = "none" + + + +ActiveThorns = "Dissipation" + +Dissipation::order = 5 +Dissipation::vars = " + ML_BSSN::ML_log_confac + ML_BSSN::ML_metric + ML_BSSN::ML_trace_curv + ML_BSSN::ML_curv + ML_BSSN::ML_Gamma + ML_BSSN::ML_lapse + ML_BSSN::ML_shift + ML_BSSN::ML_dtlapse + ML_BSSN::ML_dtshift +" + + + +ActiveThorns = "ML_ADMConstraints" + + + +ActiveThorns = "CarpetIOBasic" + +IOBasic::outInfo_every = 1 +IOBasic::outInfo_reductions = "norm2" +IOBasic::outInfo_vars = " + Carpet::physical_time_per_hour + ML_ADMConstraints::H +" + + + +ActiveThorns = "CarpetIOScalar" + +IOScalar::one_file_per_group = yes + +IOScalar::outScalar_every = 1 +IOScalar::outScalar_vars = " + CarpetReduce::weight + ADMBase::metric + ADMBase::curv + ADMBase::lapse + ADMBase::shift + ADMBase::dtlapse + ADMBase::dtshift + ML_ADMConstraints::ML_Ham + ML_ADMConstraints::ML_mom +" + + + +ActiveThorns = "CarpetIOASCII" + +IOASCII::one_file_per_group = yes + +IOASCII::output_symmetry_points = no +IOASCII::out3D_ghosts = no + +IOASCII::out0D_every = 1 +IOASCII::out0D_vars = " + Carpet::timing +" + +IOASCII::out1D_every = 10 +IOASCII::out1D_vars = " + CarpetReduce::weight + ADMBase::metric + ADMBase::curv + ADMBase::lapse + ADMBase::shift + ADMBase::dtlapse + ADMBase::dtshift + ML_ADMConstraints::ML_Ham + ML_ADMConstraints::ML_mom +" + + + +ActiveThorns = "Formaline" + + + +ActiveThorns = "TimerReport" + +TimerReport::out_every = 10 +TimerReport::out_filename = "TimerReport" +TimerReport::output_all_timers_together = yes +TimerReport::output_all_timers_readable = yes +TimerReport::n_top_timers = 20 Directory: /trunk/ ================== File [modified]: interface.ccl Delta lines: +1 -1 =================================================================== --- trunk/interface.ccl 2010-06-17 19:59:27 UTC (rev 13) +++ trunk/interface.ccl 2010-07-01 01:07:10 UTC (rev 14) @@ -2,4 +2,4 @@ # $Header$ implements: idconstraintviolate -inherits: admbase, grid \ No newline at end of file +inherits: admbase, grid File [modified]: param.ccl Delta lines: +5 -5 =================================================================== --- trunk/param.ccl 2010-06-17 19:59:27 UTC (rev 13) +++ trunk/param.ccl 2010-07-01 01:07:10 UTC (rev 14) @@ -7,19 +7,19 @@ extends KEYWORD initial_data { - "constraint violating gaussian" :: "Constraint violating gaussian" + "constraint violating gaussian" :: "Constraint violating Gaussian" } private: -REAL radius "The radial position of the gaussian waves" +REAL radius "The radial position of the Gaussian wave" { - 0:* :: "Positive" + 0:* :: "Positive" } 0.0 -REAL sigma "The sigma for the gaussian wave" +REAL sigma "The sigma for the Gaussian wave" { - 0:* :: "Positive" + (0:* :: "Positive" } 0.1 REAL amplitude "The amplitude of the waves"

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 129)
by cott＠tapir.caltech.edu 30 Jun '10

30 Jun '10

User: cott Date: 2010/06/30 01:49 AM Modified: /trunk/src/ GRHydro_PPM.F90 Log: * remove unused EOS code from GRHydro_PPM.F90 File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_PPM.F90 Delta lines: +0 -18 =================================================================== --- trunk/src/GRHydro_PPM.F90 2010-06-26 13:47:35 UTC (rev 128) +++ trunk/src/GRHydro_PPM.F90 2010-06-30 06:49:09 UTC (rev 129) @@ -64,11 +64,6 @@ DECLARE_CCTK_PARAMETERS -#ifdef _EOS_BASE_INC_ -#undef _EOS_BASE_INC_ -#endif -#include "EOS_Base.inc" - CCTK_INT :: handle,poly,nx CCTK_REAL :: dx CCTK_REAL, dimension(nx) :: rho,velx,vely,velz,eps @@ -109,19 +104,6 @@ trivial_rp = .true. -!!$ Loop to set the pressure and eps -!!$ poly=0 -!!$ if (poly .eq. 0) then -!!$ do i = 1, nx -!!$ press(i) = EOS_Pressure(handle, rho(i), eps(i)) -!!$ end do -!!$ else -!!$ do i = 1, nx -!!$ press(i) = EOS_Pressure(handle, rho(i), eps(i)) -!!$ eps(i) = EOS_SpecificIntEnergy(handle, rho(i), press(i)) -!!$ end do -!!$ end if - !!$ Average slopes delta_m(a). See (1.7) of Colella and Woodward, p.178 !!$ This is the expression for an even grid.

1 0

[svn:einsteintoolkit] incoming/EOS_Omni/ (Rev. 1)
by cott＠tapir.caltech.edu 29 Jun '10

29 Jun '10

User: cott Date: 2010/06/29 09:59 AM Added: /EOS_Omni/ README, interface.ccl, param.ccl, schedule.ccl /EOS_Omni/doc/ documentation.tex /EOS_Omni/par/ /EOS_Omni/src/ EOS_Omni_Module.F90, EOS_Omni_SingleVarCalls.F90, EOS_Omni_Startup.F90, make.code.defn /EOS_Omni/test/ Log: * EOS_Omni thorn: Initial implementation providing a polytrope and a gamma-law EOS for testing. File Changes: Directory: /EOS_Omni/ ===================== File [added]: README Delta lines: +24 -0 =================================================================== --- EOS_Omni/README (rev 0) +++ EOS_Omni/README 2010-06-29 14:59:02 UTC (rev 1) @@ -0,0 +1,24 @@ +Cactus Code Thorn EOS_Omni +Author(s) : Christian D. Ott <cott(a)tapir.caltech.edu> +Maintainer(s): Christian D. Ott <cott(a)tapir.caltech.edu> +Licence : LGPL +-------------------------------------------------------------------------- + +This Thorn provides a one-fits-all EOS handler. +All EOS are actually implemented right inside this single +thorn without crazy interfaces and lots of function calls. + +eoskey controls the type of EOS being used: + +eoskey: + +1 --- polytropic EOS +2 --- gamma-law EOS +3 --- hybrid EOS +4 --- finite-T microphysical NSE EOS + + +keyerr settings: + +-1 -- called EOS routine with keytemp=1, but + routine does not allow such a call File [added]: interface.ccl Delta lines: +89 -0 =================================================================== --- EOS_Omni/interface.ccl (rev 0) +++ EOS_Omni/interface.ccl 2010-06-29 14:59:02 UTC (rev 1) @@ -0,0 +1,89 @@ +# Interface definition for thorn EOS_Omni + +implements: EOS_Omni + +inherits: + +void FUNCTION EOS_Omni_press(CCTK_INT IN eoskey, \ + CCTK_INT IN havetemp, \ + CCTK_INT IN npoints, \ + CCTK_REAL IN ARRAY rho, \ + CCTK_REAL INOUT ARRAY eps, \ + CCTK_REAL INOUT ARRAY temp, \ + CCTK_REAL IN ARRAY ye, \ + CCTK_REAL OUT ARRAY press, \ + CCTK_INT OUT ARRAY keyerr, \ + CCTK_INT OUT anyerr) + +PROVIDES FUNCTION EOS_Omni_press WITH EOS_Omni_EOS_Press LANGUAGE Fortran + + +void FUNCTION EOS_Omni_DPressByDEps(CCTK_INT IN eoskey, \ + CCTK_INT IN havetemp, \ + CCTK_INT IN npoints, \ + CCTK_REAL IN ARRAY rho, \ + CCTK_REAL INOUT ARRAY eps, \ + CCTK_REAL INOUT ARRAY temp, \ + CCTK_REAL IN ARRAY ye, \ + CCTK_REAL OUT ARRAY DPressByDEps, \ + CCTK_INT OUT ARRAY keyerr, \ + CCTK_INT OUT anyerr) + +PROVIDES FUNCTION EOS_Omni_DPressByDEps WITH EOS_Omni_EOS_DPressByDEps LANGUAGE Fortran + + +void FUNCTION EOS_Omni_DPressByDRho(CCTK_INT IN eoskey, \ + CCTK_INT IN havetemp, \ + CCTK_INT IN npoints, \ + CCTK_REAL IN ARRAY rho, \ + CCTK_REAL INOUT ARRAY eps, \ + CCTK_REAL INOUT ARRAY temp, \ + CCTK_REAL IN ARRAY ye, \ + CCTK_REAL OUT ARRAY DPressByDRho, \ + CCTK_INT OUT ARRAY keyerr, \ + CCTK_INT OUT anyerr) + +PROVIDES FUNCTION EOS_Omni_DPressByDRho WITH EOS_Omni_EOS_DPressByDRho LANGUAGE Fortran + + +void FUNCTION EOS_Omni_cs2(CCTK_INT IN eoskey, \ + CCTK_INT IN havetemp, \ + CCTK_INT IN npoints, \ + CCTK_REAL IN ARRAY rho, \ + CCTK_REAL INOUT ARRAY eps, \ + CCTK_REAL INOUT ARRAY temp, \ + CCTK_REAL IN ARRAY ye, \ + CCTK_REAL OUT ARRAY cs2, \ + CCTK_INT OUT ARRAY keyerr, \ + CCTK_INT OUT anyerr) + +PROVIDES FUNCTION EOS_Omni_cs2 WITH EOS_Omni_EOS_cs2 LANGUAGE Fortran + + +void FUNCTION EOS_Omni_EpsFromPress(CCTK_INT IN eoskey, \ + CCTK_INT IN havetemp, \ + CCTK_INT IN npoints, \ + CCTK_REAL IN ARRAY rho, \ + CCTK_REAL INOUT ARRAY eps, \ + CCTK_REAL INOUT ARRAY temp, \ + CCTK_REAL IN ARRAY ye, \ + CCTK_REAL IN ARRAY press, \ + CCTK_REAL OUT ARRAY xeps, \ + CCTK_INT OUT ARRAY keyerr, \ + CCTK_INT OUT anyerr) + +PROVIDES FUNCTION EOS_Omni_EpsFromPress WITH EOS_Omni_EOS_eps_from_press LANGUAGE Fortran + + +void FUNCTION EOS_Omni_RestMassDensityFromEpsPress(CCTK_INT IN eoskey, \ + CCTK_INT IN havetemp, \ + CCTK_INT IN npoints, \ + CCTK_REAL OUT ARRAY rho, \ + CCTK_REAL IN ARRAY eps, \ + CCTK_REAL INOUT ARRAY temp, \ + 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 File [added]: param.ccl Delta lines: +26 -0 =================================================================== --- EOS_Omni/param.ccl (rev 0) +++ EOS_Omni/param.ccl 2010-06-29 14:59:02 UTC (rev 1) @@ -0,0 +1,26 @@ +# Parameter definitions for thorn EOS_Omni + + +restricted: + +# poly EOS +CCTK_REAL poly_gamma "Adiabatic Index for poly EOS" +{ + : :: "" +} 2.0 + +REAL poly_k "Polytropic constant in c=G=Msun=1" +{ + : :: "" +} 100.0 + +# gamma-law EOS +CCTK_REAL gl_gamma "Adiabatic Index for gamma-law EOS" +{ + : :: "" +} 2.0 + +REAL gl_k "Polytropic constant in c=G=Msun=1 for gamma-law EOS" +{ + : :: "" +} 100.0 File [added]: schedule.ccl Delta lines: +8 -0 =================================================================== --- EOS_Omni/schedule.ccl (rev 0) +++ EOS_Omni/schedule.ccl 2010-06-29 14:59:02 UTC (rev 1) @@ -0,0 +1,8 @@ +# Schedule definitions for thorn EOS_Omni + + +schedule EOS_Omni_Startup AT WRAGH +{ + LANG: Fortran + OPTIONS: global +} "Set up conversion factors and other fun stuff" Directory: /EOS_Omni/doc/ ========================= File [added]: documentation.tex Delta lines: +144 -0 =================================================================== --- EOS_Omni/doc/documentation.tex (rev 0) +++ EOS_Omni/doc/documentation.tex 2010-06-29 14:59:02 UTC (rev 1) @@ -0,0 +1,144 @@ +% *======================================================================* +% Cactus Thorn template for ThornGuide documentation +% Author: Ian Kelley +% Date: Sun Jun 02, 2002 +% $Header$ +% +% Thorn documentation in the latex file doc/documentation.tex +% will be included in ThornGuides built with the Cactus make system. +% The scripts employed by the make system automatically include +% pages about variables, parameters and scheduling parsed from the +% relevant thorn CCL files. +% +% This template contains guidelines which help to assure that your +% documentation will be correctly added to ThornGuides. More +% information is available in the Cactus UsersGuide. +% +% Guidelines: +% - Do not change anything before the line +% % START CACTUS THORNGUIDE", +% except for filling in the title, author, date, etc. fields. +% - Each of these fields should only be on ONE line. +% - Author names should be separated with a \\ or a comma. +% - You can define your own macros, but they must appear after +% the START CACTUS THORNGUIDE line, and must not redefine standard +% latex commands. +% - To avoid name clashes with other thorns, 'labels', 'citations', +% 'references', and 'image' names should conform to the following +% convention: +% ARRANGEMENT_THORN_LABEL +% For example, an image wave.eps in the arrangement CactusWave and +% thorn WaveToyC should be renamed to CactusWave_WaveToyC_wave.eps +% - Graphics should only be included using the graphicx package. +% More specifically, with the "\includegraphics" command. Do +% not specify any graphic file extensions in your .tex file. This +% will allow us to create a PDF version of the ThornGuide +% via pdflatex. +% - References should be included with the latex "\bibitem" command. +% - Use \begin{abstract}...\end{abstract} instead of \abstract{...} +% - Do not use \appendix, instead include any appendices you need as +% standard sections. +% - For the benefit of our Perl scripts, and for future extensions, +% please use simple latex. +% +% *======================================================================* +% +% Example of including a graphic image: +% \begin{figure}[ht] +% \begin{center} +% \includegraphics[width=6cm]{MyArrangement_MyThorn_MyFigure} +% \end{center} +% \caption{Illustration of this and that} +% \label{MyArrangement_MyThorn_MyLabel} +% \end{figure} +% +% Example of using a label: +% \label{MyArrangement_MyThorn_MyLabel} +% +% Example of a citation: +% \cite{MyArrangement_MyThorn_Author99} +% +% Example of including a reference +% \bibitem{MyArrangement_MyThorn_Author99} +% {J. Author, {\em The Title of the Book, Journal, or periodical}, 1 (1999), +% 1--16. {\tt http://www.nowhere.com/}} +% +% *======================================================================* + +% If you are using CVS use this line to give version information +% $Header$ + +\documentclass{article} + +% Use the Cactus ThornGuide style file +% (Automatically used from Cactus distribution, if you have a +% thorn without the Cactus Flesh download this from the Cactus +% homepage at www.cactuscode.org) +\usepackage{../../../../doc/latex/cactus} + +\begin{document} + +% The author of the documentation +\author{Christian D. Ott \textless cott(a)tapir.caltech.edu\textgreater} + +% The title of the document (not necessarily the name of the Thorn) +\title{EOS\_Omni} + +% the date your document was last changed, if your document is in CVS, +% please use: +% \date{$ $Date: 2004-01-07 12:12:39 -0800 (Wed, 07 Jan 2004) $ $} +\date{June 27 2010} + +\maketitle + +% Do not delete next line +% START CACTUS THORNGUIDE + +% Add all definitions used in this documentation here +% \def\mydef etc + +% Add an abstract for this thorn's documentation +\begin{abstract} + +\end{abstract} + +% The following sections are suggestive only. +% Remove them or add your own. + +\section{Introduction} + +\section{Physical System} + +\section{Numerical Implementation} + +\section{Using This Thorn} + +\subsection{Obtaining This Thorn} + +\subsection{Basic Usage} + +\subsection{Special Behaviour} + +\subsection{Interaction With Other Thorns} + +\subsection{Examples} + +\subsection{Support and Feedback} + +\section{History} + +\subsection{Thorn Source Code} + +\subsection{Thorn Documentation} + +\subsection{Acknowledgements} + + +\begin{thebibliography}{9} + +\end{thebibliography} + +% Do not delete next line +% END CACTUS THORNGUIDE + +\end{document} Directory: /EOS_Omni/src/ ========================= File [added]: EOS_Omni_Module.F90 Delta lines: +22 -0 =================================================================== --- EOS_Omni/src/EOS_Omni_Module.F90 (rev 0) +++ EOS_Omni/src/EOS_Omni_Module.F90 2010-06-29 14:59:02 UTC (rev 1) @@ -0,0 +1,22 @@ +module EOS_Omni_Module + + implicit none + + real*8,parameter :: rho_gf = 1.61930347d-18 + real*8,parameter :: press_gf = 1.80171810d-39 + real*8,parameter :: eps_gf = 1.11265006d-21 + real*8,parameter :: time_gf = 2.03001708d+05 + real*8,parameter :: mass_gf = 5.02765209d-34 + real*8,parameter :: length_gf = 6.77140812d-06 + + real*8,parameter :: inv_rho_gf = 6.17549470205236d17 + real*8,parameter :: inv_press_gf = 5.55025783445257d38 + real*8,parameter :: inv_eps_gf = 8.98755175549085d20 + real*8,parameter :: inv_time_gf = 4.92606692747629d-6 + real*8,parameter :: inv_mass_gf = 1.98899999860571d33 + real*8,parameter :: inv_length_gf = 147679.77092481d0 + + real*8 :: poly_k_cgs = 0.0d0 + real*8 :: gl_k_cgs = 0.0d0 + +end module EOS_Omni_Module File [added]: EOS_Omni_SingleVarCalls.F90 Delta lines: +333 -0 =================================================================== --- EOS_Omni/src/EOS_Omni_SingleVarCalls.F90 (rev 0) +++ EOS_Omni/src/EOS_Omni_SingleVarCalls.F90 2010-06-29 14:59:02 UTC (rev 1) @@ -0,0 +1,333 @@ +#include "cctk.h" +#include "cctk_Parameters.h" +#include "cctk_Arguments.h" +#include "cctk_Functions.h" + +! eoskey: +! 1 --- polytropic EOS +! 2 --- gamma-law EOS +! 3 --- hybrid EOS +! 4 --- finite-T microphysical NSE EOS + +subroutine EOS_Omni_EOS_Press(eoskey,keytemp,npoints,& + rho,eps,temp,ye,press,keyerr,anyerr) + + use EOS_Omni_Module + implicit none + DECLARE_CCTK_PARAMETERS + + CCTK_INT, intent(in) :: eoskey,keytemp,npoints + CCTK_INT, intent(out) :: keyerr(npoints) + CCTK_INT, intent(out) :: anyerr + CCTK_REAL, intent(in) :: rho(npoints),ye(npoints) + CCTK_REAL, intent(inout) :: eps(npoints), temp(npoints) + CCTK_REAL, intent(out) :: press(npoints) + + ! local vars + integer :: i + character(256) :: warnstring + + anyerr = 0 + keyerr(:) = 0 + + select case (eoskey) + case (1) + ! polytropic EOS + if(keytemp.eq.1) then + do i=1,npoints + eps(i) = press_gf * poly_k_cgs * & + (rho(i)*inv_rho_gf)**(poly_gamma) / & + (poly_gamma - 1.0d0) / rho(i) + enddo + endif + do i=1,npoints + press(i) = press_gf * poly_k_cgs * & + (rho(i)*inv_rho_gf)**poly_gamma + enddo + case (2) + ! gamma-law EOS + if(keytemp.eq.1) then + do i=1,npoints + eps(i) = press_gf * gl_k_cgs * & + (rho(i)*inv_rho_gf)**(gl_gamma) / & + (gl_gamma - 1.0d0) / rho(i) + enddo + endif + do i=1,npoints + press(i) = (gl_gamma - 1.0d0) * rho(i) * eps(i) + enddo + + case DEFAULT + write(warnstring,*) "eoskey ",eoskey," not implemented!" + call CCTK_WARN(0,warnstring) + end select + +end subroutine EOS_Omni_EOS_Press + + + +subroutine EOS_Omni_EOS_DPressByDRho(eoskey,keytemp,npoints,& + rho,eps,temp,ye,dpdrhoe,keyerr,anyerr) + + use EOS_Omni_Module + implicit none + DECLARE_CCTK_PARAMETERS + + CCTK_INT, intent(in) :: eoskey,keytemp,npoints + CCTK_INT, intent(out) :: keyerr(npoints) + CCTK_INT, intent(out) :: anyerr + CCTK_REAL, intent(in) :: rho(npoints),ye(npoints) + CCTK_REAL, intent(inout) :: eps(npoints), temp(npoints) + CCTK_REAL, intent(out) :: dpdrhoe(npoints) + + ! local vars + integer :: i + character(256) :: warnstring + + anyerr = 0 + keyerr(:) = 0 + + select case (eoskey) + case (1) + ! polytropic EOS + if(keytemp.eq.1) then + do i=1,npoints + eps(i) = press_gf * poly_k_cgs * & + (rho(i)*inv_rho_gf)**(poly_gamma) / & + (poly_gamma - 1.0d0) / rho(i) + enddo + endif + do i=1,npoints + dpdrhoe(i) = press_gf * poly_k_cgs * & + poly_gamma * inv_rho_gf * & + (rho(i)*inv_rho_gf) ** (poly_gamma - 1.d0) + enddo + case (2) + ! gamma-law EOS + if(keytemp.eq.1) then + do i=1,npoints + eps(i) = press_gf * gl_k_cgs * & + (rho(i)*inv_rho_gf)**(gl_gamma) / & + (gl_gamma - 1.0d0) / rho(i) + enddo + endif + do i=1,npoints + dpdrhoe(i) = (gl_gamma-1.0d0) * & + eps(i) + enddo + + case DEFAULT + write(warnstring,*) "eoskey ",eoskey," not implemented!" + call CCTK_WARN(0,warnstring) + end select + +end subroutine EOS_Omni_EOS_DPressByDRho + +subroutine EOS_Omni_EOS_DPressByDEps(eoskey,keytemp,npoints,& + rho,eps,temp,ye,dpdepsrho,keyerr,anyerr) + + use EOS_Omni_Module + implicit none + DECLARE_CCTK_PARAMETERS + + CCTK_INT, intent(in) :: eoskey,keytemp,npoints + CCTK_INT, intent(out) :: keyerr(npoints) + CCTK_INT, intent(out) :: anyerr + CCTK_REAL, intent(in) :: rho(npoints),ye(npoints) + CCTK_REAL, intent(inout) :: eps(npoints), temp(npoints) + CCTK_REAL, intent(out) :: dpdepsrho(npoints) + + ! local vars + integer :: i + character(256) :: warnstring + + anyerr = 0 + keyerr(:) = 0 + + select case (eoskey) + case (1) + ! polytropic EOS + if(keytemp.eq.1) then + do i=1,npoints + eps(i) = press_gf * poly_k_cgs * & + (rho(i)*inv_rho_gf)**(poly_gamma) / & + (poly_gamma - 1.0d0) / rho(i) + enddo + endif + do i=1,npoints + dpdepsrho(i) = 0.0d0 + enddo + case (2) + ! gamma-law EOS + if(keytemp.eq.1) then + do i=1,npoints + eps(i) = press_gf * gl_k_cgs * & + (rho(i)*inv_rho_gf)**(gl_gamma) / & + (gl_gamma - 1.0d0) / rho(i) + enddo + endif + do i=1,npoints + dpdepsrho(i) = (gl_gamma - 1.0d0) * & + rho(i) + enddo + + case DEFAULT + write(warnstring,*) "eoskey ",eoskey," not implemented!" + call CCTK_WARN(0,warnstring) + end select + +end subroutine EOS_Omni_EOS_DPressByDEps + +subroutine EOS_Omni_EOS_cs2(eoskey,keytemp,npoints,& + rho,eps,temp,ye,cs2,keyerr,anyerr) + + use EOS_Omni_Module + implicit none + DECLARE_CCTK_PARAMETERS + + CCTK_INT, intent(in) :: eoskey,keytemp,npoints + CCTK_INT, intent(out) :: keyerr(npoints) + CCTK_INT, intent(out) :: anyerr + CCTK_REAL, intent(in) :: rho(npoints),ye(npoints) + CCTK_REAL, intent(inout) :: eps(npoints), temp(npoints) + CCTK_REAL, intent(out) :: cs2(npoints) + + ! local vars + integer :: i + character(256) :: warnstring + real*8 :: xpress,xdpdrhoe,xdpderho + + anyerr = 0 + keyerr(:) = 0 + + select case (eoskey) + case (1) + ! polytropic EOS + if(keytemp.eq.1) then + do i=1,npoints + eps(i) = press_gf * poly_k_cgs * & + (rho(i)*inv_rho_gf)**(poly_gamma) / & + (poly_gamma - 1.0d0) / rho(i) + enddo + endif + do i=1,npoints + xpress = press_gf*poly_k_cgs * & + (rho(i)*inv_rho_gf)**(poly_gamma) + cs2(i) = poly_gamma * xpress / rho(i) / & + (1 + eps(i) + xpress/rho(i)) + enddo + case (2) + ! gamma-law EOS + if(keytemp.eq.1) then + do i=1,npoints + eps(i) = press_gf * gl_k_cgs * & + (rho(i)*inv_rho_gf)**(gl_gamma) / & + (gl_gamma - 1.0d0) / rho(i) + enddo + endif + do i=1,npoints + xpress = (gl_gamma-1.0d0)*rho(i)*eps(i) + xdpdrhoe = (gl_gamma-1.0d0)*eps(i) + xdpderho = (gl_gamma-1.0d0)*rho(i) + cs2(i) = (xdpdrhoe + xpress * xdpderho / (rho(i)**2)) / & + (1.0d0 + eps(i) + xpress/rho(i)) + enddo + + case DEFAULT + write(warnstring,*) "eoskey ",eoskey," not implemented!" + call CCTK_WARN(0,warnstring) + end select + +end subroutine EOS_Omni_EOS_cs2 + + +subroutine EOS_Omni_EOS_eps_from_press(eoskey,keytemp,npoints,& + rho,eps,temp,ye,press,xeps,keyerr,anyerr) + + use EOS_Omni_Module + implicit none + DECLARE_CCTK_PARAMETERS + + CCTK_INT, intent(in) :: eoskey,keytemp,npoints + CCTK_INT, intent(out) :: keyerr(npoints) + CCTK_INT, intent(out) :: anyerr + CCTK_REAL, intent(in) :: rho(npoints),ye(npoints),press(npoints) + CCTK_REAL, intent(inout) :: eps(npoints), temp(npoints) + CCTK_REAL, intent(out) :: xeps(npoints) + + ! local vars + integer :: i + character(256) :: warnstring + + if(keytemp.eq.1) then + anyerr = 1 + keyerr(:) = -1 + else + anyerr = 0 + keyerr(:) = 0 + endif + + select case (eoskey) + case (1) + ! polytropic EOS + do i=1,npoints + xeps(i) = press(i) / (poly_gamma - 1.0d0) / rho(i) + eps(i) = xeps(i) + enddo + case (2) + ! gamma-law EOS + do i=1,npoints + xeps(i) = press(i) / (gl_gamma - 1.0d0) / rho(i) + eps(i) = xeps(i) + enddo + case DEFAULT + write(warnstring,*) "eoskey ",eoskey," not implemented!" + call CCTK_WARN(0,warnstring) + end select + +end subroutine EOS_Omni_EOS_eps_from_press + + +subroutine EOS_Omni_EOS_RestMassDensityFromEpsPress(eoskey,keytemp,npoints,& + rho,eps,temp,ye,press,keyerr,anyerr) + + use EOS_Omni_Module + implicit none + DECLARE_CCTK_PARAMETERS + + CCTK_INT, intent(in) :: eoskey,keytemp,npoints + CCTK_INT, intent(out) :: keyerr(npoints) + CCTK_INT, intent(out) :: anyerr + CCTK_REAL, intent(in) :: ye(npoints),press(npoints),eps(npoints) + CCTK_REAL, intent(inout) :: rho(npoints),temp(npoints) + + + ! local vars + integer :: i + character(256) :: warnstring + + if(keytemp.eq.1) then + anyerr = 1 + keyerr(:) = -1 + else + anyerr = 0 + keyerr(:) = 0 + endif + + select case (eoskey) + case (1) + ! polytropic EOS + do i=1,npoints + rho(i) = press(i) / ((poly_gamma - 1.0d0)*eps(i)) + enddo + case (2) + ! gamma-law EOS + do i=1,npoints + rho(i) = press(i) / ((gl_gamma - 1.0d0)*eps(i)) + enddo + case DEFAULT + write(warnstring,*) "eoskey ",eoskey," not implemented!" + call CCTK_WARN(0,warnstring) + end select + +end subroutine EOS_Omni_EOS_RestMassDensityFromEpsPress File [added]: EOS_Omni_Startup.F90 Delta lines: +18 -0 =================================================================== --- EOS_Omni/src/EOS_Omni_Startup.F90 (rev 0) +++ EOS_Omni/src/EOS_Omni_Startup.F90 2010-06-29 14:59:02 UTC (rev 1) @@ -0,0 +1,18 @@ +#include "cctk.h" +#include "cctk_Parameters.h" +#include "cctk_Arguments.h" +#include "cctk_Functions.h" + +subroutine EOS_Omni_Startup(CCTK_ARGUMENTS) + + use EOS_Omni_Module + implicit none + + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_ARGUMENTS + + poly_k_cgs = poly_k * rho_gf**poly_gamma / press_gf + + gl_k_cgs = gl_k * rho_gf**gl_gamma / press_gf + +end subroutine EOS_Omni_Startup File [added]: make.code.defn Delta lines: +7 -0 =================================================================== --- EOS_Omni/src/make.code.defn (rev 0) +++ EOS_Omni/src/make.code.defn 2010-06-29 14:59:02 UTC (rev 1) @@ -0,0 +1,7 @@ +# Main make.code.defn file for thorn EOS_Omni + +# Source files in this directory +SRCS = EOS_Omni_Module.F90 EOS_Omni_Startup.F90 EOS_Omni_SingleVarCalls.F90 + +# Subdirectories containing source files +SUBDIRS =

1 0

[svn:einsteintoolkit] www/info/support/ (Rev. 308)
by gallen＠cct.lsu.edu 28 Jun '10

28 Jun '10

User: gallen Date: 2010/06/28 08:55 AM Modified: /info/support/ index.php Log: Tweak to phone details File Changes: Directory: /info/support/ ========================= File [modified]: index.php Delta lines: +2 -2 =================================================================== --- info/support/index.php 2010-06-27 04:34:14 UTC (rev 307) +++ info/support/index.php 2010-06-28 13:55:19 UTC (rev 308) @@ -38,8 +38,8 @@ call is 10am (central US time) every Monday. Instructions how to dial in are given in the following: <h4>Dial-in Instructions</h4> -Each participant of the teleconference will call 225-578-4942 -(8-4942 from LSU-campus or 866-573-0359 toll free). +Each participant of the teleconference will call (+1) 225-578-4942 +((+1) 866-573-0359 toll free from inside the USA). They will hear: <pre>

1 0

Jump to page:

2026

2025

2024

2023

2022

2021

2020

2019

2018

2017

2016

2015

2014

2013

2012

2011

2010

Commits