[Commits] [svn:einsteintoolkit] EOS_Omni/trunk/ (Rev. 84)
cott at tapir.caltech.edu
cott at tapir.caltech.edu
Thu Jun 13 18:10:44 CDT 2013
User: cott
Date: 2013/06/13 06:10 PM
Modified:
/trunk/
interface.ccl
/trunk/src/
EOS_Omni_MultiVarCalls.F90, EOS_Omni_SingleVarCalls.F90
Log:
* remove OMP pragmas from most routines, because they
slow down pointwise calls
* make a copy of the press routine, call it OMP, add interface
entry; use this routine for GF-wide calls.
* uglyness ensues
File Changes:
Directory: /trunk/src/
======================
File [modified]: EOS_Omni_MultiVarCalls.F90
Delta lines: +0 -30
===================================================================
--- trunk/src/EOS_Omni_MultiVarCalls.F90 2013-06-08 00:23:17 UTC (rev 83)
+++ trunk/src/EOS_Omni_MultiVarCalls.F90 2013-06-13 23:10:44 UTC (rev 84)
@@ -47,8 +47,6 @@
anyerr = 0
keyerr(:) = 0
- !$OMP PARALLEL DO PRIVATE(xrho,xtemp,xye,xenr,xent,xcs2,&
- !$OMP xdedt,xdpderho,xdpdrhoe,xmunu)
do i=1,npoints
xrho = rho(i) * inv_rho_gf
@@ -61,9 +59,7 @@
keytemp,keyerr(i),rf_precision)
if(keyerr(i).ne.0) then
- !$OMP CRITICAL
anyerr = 1
- !$OMP END CRITICAL
endif
if(keytemp.eq.1) then
@@ -84,7 +80,6 @@
munu(i) = xmunu
enddo
- !$OMP END PARALLEL DO
end subroutine EOS_Omni_EOS_short
@@ -123,45 +118,35 @@
case (1)
! polytropic EOS
if(keytemp.eq.1) then
- !$OMP PARALLEL DO
do i=1,npoints
eps(i) = poly_k * &
rho(i)**(poly_gamma) / &
(poly_gamma - 1.0d0) / rho(i)
enddo
- !$OMP END PARALLEL DO
endif
- !$OMP PARALLEL DO
do i=1,npoints
dpdrhoe(i) = poly_k * &
poly_gamma * &
rho(i) ** (poly_gamma - 1.d0)
dpderho(i) = 0.0d0
enddo
- !$OMP END PARALLEL DO
case (2)
! gamma-law EOS
if(keytemp.eq.1) then
- !$OMP PARALLEL DO
do i=1,npoints
eps(i) = gl_k * &
rho(i)**(gl_gamma) / &
(gl_gamma - 1.0d0) / rho(i)
enddo
- !$OMP END PARALLEL DO
endif
- !$OMP PARALLEL DO
do i=1,npoints
dpdrhoe(i) = (gl_gamma-1.0d0) * &
eps(i)
dpderho(i) = (gl_gamma - 1.0d0) * &
rho(i)
enddo
- !$OMP END PARALLEL DO
case (3)
! hybrid EOS
- !$OMP PARALLEL DO PRIVATE(hybrid_local_gamma,hybrid_local_k,&
- !$OMP hybrid_dp_poly,hybrid_dp_th1,hybrid_dp_th2)
do i=1,npoints
if(rho(i).gt.hybrid_rho_nuc) then
hybrid_local_gamma = hybrid_gamma2
@@ -187,9 +172,7 @@
dpdrhoe(i) = hybrid_dp_poly + hybrid_dp_th1 + hybrid_dp_th2
dpderho(i) = (hybrid_gamma_th - 1.0d0) * rho(i)
enddo
- !$OMP END PARALLEL DO
case (4)
- !$OMP PARALLEL DO PRIVATE(xrho,xtemp,xye,xenr,xdpderho,xdpdrhoe)
do i=1,npoints
xrho = rho(i) * inv_rho_gf
xtemp = temp(i)
@@ -200,9 +183,7 @@
keytemp,keyerr(i),rf_precision)
if(keyerr(i).ne.0) then
- !$OMP CRITICAL
anyerr = 1
- !$OMP END CRITICAL
endif
if(keytemp.eq.1) then
@@ -257,45 +238,34 @@
case (1)
! polytropic EOS
if(keytemp.eq.1) then
- !$OMP PARALLEL DO
do i=1,npoints
eps(i) = poly_k * &
rho(i)**(poly_gamma) / &
(poly_gamma - 1.0d0) / rho(i)
enddo
- !$OMP END PARALLEL DO
endif
- !$OMP PARALLEL DO
do i=1,npoints
depsdpress(i) = 1.0d0/(poly_gamma - 1.0d0)/rho(i)
depsdrho(i) = depsdpress(i) * poly_k * &
poly_gamma * &
rho(i) ** (poly_gamma - 1.d0)
enddo
- !$OMP END PARALLEL DO
case (2)
! gamma-law EOS
if(keytemp.eq.1) then
- !$OMP PARALLEL DO
do i=1,npoints
eps(i) = gl_k * &
rho(i)**(gl_gamma) / &
(gl_gamma - 1.0d0) / rho(i)
enddo
- !$OMP END PARALLEL DO
endif
- !$OMP PARALLEL DO
do i=1,npoints
depsdpress(i) = 1.0/( (gl_gamma - 1.0d0) * &
rho(i))
depsdrho(i) = -eps(i)/rho(i)
enddo
- !$OMP END PARALLEL DO
case (3)
! hybrid EOS
- !$OMP PARALLEL DO PRIVATE(hybrid_local_gamma,hybrid_local_k,&
- !$OMP hybrid_dp_poly,hybrid_dp_th1,&
- !$OMP hybrid_dp_th2,xdpdrhoe,xdpderho)
do i=1,npoints
if(rho(i).gt.hybrid_rho_nuc) then
hybrid_local_gamma = hybrid_gamma2
File [modified]: EOS_Omni_SingleVarCalls.F90
Delta lines: +162 -66
===================================================================
--- trunk/src/EOS_Omni_SingleVarCalls.F90 2013-06-08 00:23:17 UTC (rev 83)
+++ trunk/src/EOS_Omni_SingleVarCalls.F90 2013-06-13 23:10:44 UTC (rev 84)
@@ -47,6 +47,156 @@
case (1)
! polytropic EOS
if(keytemp.eq.1) then
+ do i=1,npoints
+ eps(i) = poly_k * &
+ rho(i)**(poly_gamma) / &
+ (poly_gamma - 1.0d0) / rho(i)
+ enddo
+ endif
+ do i=1,npoints
+ press(i) = poly_k * &
+ rho(i)**poly_gamma
+ enddo
+ case (2)
+ ! gamma-law EOS
+ if(keytemp.eq.1) then
+ do i=1,npoints
+ eps(i) = gl_k * &
+ rho(i)**(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 (3)
+ ! hybrid EOS
+ do i=1,npoints
+ if(rho(i).gt.hybrid_rho_nuc) then
+ hybrid_local_gamma = hybrid_gamma2
+ hybrid_local_k = hybrid_k2
+ else
+ hybrid_local_gamma = hybrid_gamma1
+ hybrid_local_k = hybrid_k1
+ endif
+ hybrid_p_poly = hybrid_local_k * &
+ rho(i)**hybrid_local_gamma
+ hybrid_p_th = - hybrid_local_k * (hybrid_gamma_th - 1.d0) / &
+ (hybrid_local_gamma - 1.0d0) * rho(i)**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) * &
+ hybrid_k1 * &
+ 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 (4)
+ ! nuc eos
+ do i=1,npoints
+ xrho = rho(i) * inv_rho_gf
+ xtemp = temp(i)
+ xye = ye(i)
+ xenr = eps(i) * inv_eps_gf
+ call nuc_eos_press_eps(xrho,xtemp,xye,xenr,xprs,&
+ keytemp,keyerr(i),rf_precision)
+
+ if(keyerr(i).ne.0) then
+ anyerr = 1
+ endif
+
+ if(keytemp.eq.1) then
+ eps(i) = xenr * eps_gf
+ else
+ temp(i) = xtemp
+ endif
+
+ press(i) = xprs * press_gf
+ enddo
+
+ case (5)
+ ! cold tabular EOS with gamma law
+ do i=1,npoints
+ if(rho(i).lt.coldeos_rhomin) then
+ press(i) = coldeos_low_kappa * rho(i)**coldeos_low_gamma
+ eps(i) = press(i) / (coldeos_low_gamma - 1.0) / rho(i)
+ cycle
+ else if(rho(i).gt.coldeos_rhomax) then
+ keyerr(i) = 103
+ anyerr = 1
+ else
+ xrho = log10(rho(i))
+ ir = 2 + INT( (xrho - coldeos_logrho(1) - 1.0d-10) * coldeos_dlrhoi )
+ endif
+ ir = max(2, min(ir,coldeos_nrho))
+
+ gamma = (coldeos_gamma(ir) - coldeos_gamma(ir-1)) / &
+ (coldeos_logrho(ir) - coldeos_logrho(ir-1)) * &
+ (xrho - coldeos_logrho(ir-1)) + coldeos_gamma(ir-1)
+
+ eps_cold = (coldeos_eps(ir) - coldeos_eps(ir-1)) / &
+ (coldeos_logrho(ir) - coldeos_logrho(ir-1)) * &
+ (xrho - coldeos_logrho(ir-1)) + coldeos_eps(ir-1)
+
+ if(keytemp.eq.0) then
+ eps_th = max(0.0d0,eps(i) - eps_cold)
+ else if (keytemp.eq.1) then
+ eps_th = 0.0d0
+ eps(i) = eps_cold
+ endif
+
+ press_cold = coldeos_kappa * rho(i)**gamma
+ press_th = coldeos_thfac*(coldeos_gammath - 1.0d0)*rho(i)*eps_th
+ press(i) = press_cold + press_th
+ 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_PressOMP(eoskey,keytemp,rf_precision,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) :: rf_precision
+ 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
+ real*8 :: hybrid_local_gamma, hybrid_local_k, &
+ hybrid_p_poly, hybrid_p_th
+ real*8,parameter :: zero = 0.0d0
+ ! temporary vars for nuc_eos
+ real*8 :: xrho,xye,xtemp,xenr,xent
+ real*8 :: xprs,xmunu,xcs2
+ real*8 :: xdedt,xdpderho,xdpdrhoe
+ ! temporary vars for coldeos + gamma law
+ integer :: ir
+ real*8 :: press_cold, press_th
+ real*8 :: eps_cold, eps_th
+ real*8 :: gamma
+
+ anyerr = 0
+ keyerr(:) = 0
+
+ select case (eoskey)
+ case (1)
+ ! polytropic EOS
+ if(keytemp.eq.1) then
!$OMP PARALLEL DO
do i=1,npoints
eps(i) = poly_k * &
@@ -177,10 +327,9 @@
call CCTK_WARN(0,warnstring)
end select
-end subroutine EOS_Omni_EOS_Press
+end subroutine EOS_Omni_EOS_PressOMP
-
subroutine EOS_Omni_EOS_DPressByDRho(eoskey,keytemp,rf_precision,npoints,&
rho,eps,temp,ye,dpdrhoe,keyerr,anyerr)
@@ -219,43 +368,32 @@
case (1)
! polytropic EOS
if(keytemp.eq.1) then
- !$OMP PARALLEL DO
do i=1,npoints
eps(i) = poly_k * &
rho(i)**(poly_gamma) / &
(poly_gamma - 1.0d0) / rho(i)
enddo
- !$OMP END PARALLEL DO
endif
- !$OMP PARALLEL DO
do i=1,npoints
dpdrhoe(i) = poly_k * &
poly_gamma * &
rho(i) ** (poly_gamma - 1.d0)
enddo
- !$OMP END PARALLEL DO
case (2)
! gamma-law EOS
if(keytemp.eq.1) then
- !$OMP PARALLEL DO
do i=1,npoints
eps(i) = gl_k * &
rho(i)**(gl_gamma) / &
(gl_gamma - 1.0d0) / rho(i)
enddo
- !$OMP END PARALLEL DO
endif
- !$OMP PARALLEL DO
do i=1,npoints
dpdrhoe(i) = (gl_gamma-1.0d0) * &
eps(i)
enddo
- !$OMP END PARALLEL DO
case (3)
! hybrid EOS
- !$OMP PARALLEL DO PRIVATE(hybrid_local_gamma,hybrid_local_k,&
- !$OMP hybrid_dp_poly, hybrid_dp_th1, &
- !$OMP hybrid_dp_th2)
do i=1,npoints
if(rho(i).gt.hybrid_rho_nuc) then
hybrid_local_gamma = hybrid_gamma2
@@ -281,10 +419,8 @@
dpdrhoe(i) = hybrid_dp_poly + max(0.0d0,hybrid_dp_th1 + hybrid_dp_th2)
enddo
- !$OMP END PARALLEL DO
+
case (4)
- !$OMP PARALLEL DO PRIVATE(xrho,xtemp,xye,xenr,xprs,xent,xcs2,xdedt,&
- !$OMP xdpderho,xdpdrhoe,xmunu,h)
do i=1,npoints
xrho = rho(i) * inv_rho_gf
xtemp = temp(i)
@@ -295,9 +431,7 @@
keytemp,keyerr(i),rf_precision)
if(keyerr(i).ne.0) then
- !$OMP CRITICAL
anyerr = 1
- !$OMP END CRITICAL
endif
if(keytemp.eq.1) then
@@ -308,13 +442,11 @@
dpdrhoe(i) = xdpdrhoe * press_gf * inv_rho_gf
enddo
- !$OMP END PARALLEL DO
+
case (5)
! with the cold eos we have to assume P = P(rho), so
! by definition dPdrho is at constant internal energy
! and entropy (the latter, because T=0)
- !$OMP PARALLEL DO PRIVATE(xrho,ir,gamma,eps_cold,eps_th, &
- !$OMP cs2,press_cold,press_th)
do i=1,npoints
if(rho(i).lt.coldeos_rhomin) then
dpdrhoe(i) = coldeos_low_kappa * coldeos_low_gamma * &
@@ -322,9 +454,7 @@
cycle
else if(rho(i).gt.coldeos_rhomax) then
keyerr(i) = 103
- !$OMP CRITICAL
anyerr = 1
- !$OMP END CRITICAL
else
xrho = log10(rho(i))
ir = 2 + INT( (xrho - coldeos_logrho(1) - 1.0d-10) * coldeos_dlrhoi )
@@ -401,49 +531,38 @@
case (1)
! polytropic EOS
if(keytemp.eq.1) then
- !$OMP PARALLEL DO
do i=1,npoints
eps(i) = poly_k * &
rho(i)**(poly_gamma) / &
(poly_gamma - 1.0d0) / rho(i)
enddo
- !$OMP END PARALLEL DO
endif
- !$OMP PARALLEL DO
do i=1,npoints
dpdepsrho(i) = 0.0d0
enddo
- !$OMP END PARALLEL DO
case (2)
! gamma-law EOS
if(keytemp.eq.1) then
- !$OMP PARALLEL DO
do i=1,npoints
eps(i) = gl_k * &
rho(i)**(gl_gamma) / &
(gl_gamma - 1.0d0) / rho(i)
enddo
- !$OMP END PARALLEL DO
endif
- !$OMP PARALLEL DO
do i=1,npoints
dpdepsrho(i) = (gl_gamma - 1.0d0) * &
rho(i)
enddo
- !$OMP END PARALLEL DO
case (3)
! hybrid EOS
- !$OMP PARALLEL DO
do i=1,npoints
dpdepsrho(i) = (hybrid_gamma_th - 1.0d0) * rho(i)
enddo
- !$OMP END PARALLEL DO
+
case (4)
! nuc_eos
- !$OMP PARALLEL DO PRIVATE(xrho,xtemp,xye,xenr,xprs,&
- !$OMP xent,xcs2,xdedt,xdpderho,xdpdrhoe,&
- !$OMP xmunu)
do i=1,npoints
+
xrho = rho(i) * inv_rho_gf
xtemp = temp(i)
xye = ye(i)
@@ -453,9 +572,7 @@
keytemp,keyerr(i),rf_precision)
if(keyerr(i).ne.0) then
- !$OMP CRITICAL
anyerr = 1
- !$OMP END CRITICAL
endif
if(keytemp.eq.1) then
@@ -467,20 +584,17 @@
dpdepsrho(i) = xdpderho * press_gf * inv_eps_gf
enddo
- !$OMP END PARALLEL DO
+
case (5)
! with the cold eos we have to assume P = P(rho), so
! only the gamma law has non-zero dPdeps
- !$OMP PARALLEL DO PRIVATE(xrho,ir,eps_cold,eps_th)
do i=1,npoints
if(rho(i).lt.coldeos_rhomin) then
dpdepsrho(i) = 0.0d0
cycle
else if(rho(i).gt.coldeos_rhomax) then
keyerr(i) = 103
- !$OMP CRITICAL
anyerr = 1
- !$OMP END CRITICAL
else
xrho = log10(rho(i))
ir = 2 + INT( (xrho - coldeos_logrho(1) - 1.0d-10) * coldeos_dlrhoi )
@@ -547,34 +661,27 @@
case (1)
! polytropic EOS
if(keytemp.eq.1) then
- !$OMP PARALLEL DO
do i=1,npoints
eps(i) = poly_k * &
rho(i)**(poly_gamma) / &
(poly_gamma - 1.0d0) / rho(i)
enddo
- !$OMP END PARALLEL DO
endif
- !$OMP PARALLEL DO
do i=1,npoints
xpress = poly_k * &
rho(i)**(poly_gamma)
cs2(i) = poly_gamma * xpress / rho(i) / &
(1 + eps(i) + xpress/rho(i))
enddo
- !$OMP END PARALLEL DO
case (2)
! gamma-law EOS
if(keytemp.eq.1) then
- !$OMP PARALLEL DO
do i=1,npoints
eps(i) = gl_k * &
rho(i)**(gl_gamma) / &
(gl_gamma - 1.0d0) / rho(i)
enddo
- !$OMP END PARALLEL DO
endif
- !$OMP PARALLEL DO PRIVATE(xpress,xdpdrhoe,xdpderho)
do i=1,npoints
xpress = (gl_gamma-1.0d0)*rho(i)*eps(i)
xdpdrhoe = (gl_gamma-1.0d0)*eps(i)
@@ -582,12 +689,8 @@
cs2(i) = (xdpdrhoe + xpress * xdpderho / (rho(i)**2)) / &
(1.0d0 + eps(i) + xpress/rho(i))
enddo
- !$OMP END PARALLEL DO
case(3)
! hybrid EOS
- !$OMP PARALLEL DO PRIVATE(hybrid_local_gamma,hybrid_local_k,&
- !$OMP hybrid_p_poly,hybrid_p_th,&
- !$OMP xpress)
do i=1,npoints
if(rho(i).gt.hybrid_rho_nuc) then
hybrid_local_gamma = hybrid_gamma2
@@ -611,12 +714,11 @@
cs2(i) = (hybrid_local_gamma * hybrid_p_poly + hybrid_gamma_th * hybrid_p_th) / &
rho(i) / (1.0d0 + eps(i) + xpress/rho(i))
enddo
- !$OMP END PARALLEL DO
case(4)
! nuc_eos
- !$OMP PARALLEL DO PRIVATE(xrho,xtemp,xye,xenr,xprs,xent,xcs2,&
- !$OMP xdedt,xdpderho,xdpdrhoe,xmunu)
+
do i=1,npoints
+
xrho = rho(i) * inv_rho_gf
xtemp = temp(i)
xye = ye(i)
@@ -626,9 +728,7 @@
keytemp,keyerr(i),rf_precision)
if(keyerr(i).ne.0) then
- !$OMP CRITICAL
anyerr = 1
- !$OMP END CRITICAL
endif
if(keytemp.eq.1) then
@@ -639,15 +739,12 @@
cs2(i) = xcs2 * cliteinv2 / &
(1.0d0 + eps(i) + xprs * press_gf / rho(i))
+
enddo
- !$OMP END PARALLEL DO
case (5)
! with the cold eos we have to assume P = P(rho), so
! by definition dPdrho is at constant internal energy
! and entropy (the latter, because T=0)
- !$OMP PARALLEL DO PRIVATE(xprs,xrho,ir,gamma,cs2_cold,&
- !$OMP eps_cold,eps_th,press_cold,press_th,&
- !$OMP xdpdrhoe,xdpderho,cs2_th,h,h_cold)
do i=1,npoints
if(rho(i).lt.coldeos_rhomin) then
xprs = coldeos_low_kappa * rho(i)**coldeos_low_gamma
@@ -658,10 +755,8 @@
cycle
else if(rho(i).gt.coldeos_rhomax) then
keyerr(i) = 103
- !$OMP CRITICAL
anyerr = 1
- !$OMP END CRITICAL
- else
+ else
xrho = log10(rho(i))
ir = 2 + INT( (xrho - coldeos_logrho(1) - 1.0d-10) * coldeos_dlrhoi )
endif
@@ -686,6 +781,7 @@
eps_th = 0.0d0
eps(i) = eps_cold
endif
+
press_cold = coldeos_kappa * rho(i)**gamma
press_th = coldeos_thfac*(coldeos_gammath - 1.0d0)*rho(i)*eps_th
@@ -698,7 +794,7 @@
cs2(i) = (cs2_cold * h_cold + cs2_th) / h
enddo
- !$OMP END PARALLEL DO
+
case DEFAULT
write(warnstring,*) "eoskey ",eoskey," not implemented!"
call CCTK_WARN(0,warnstring)
Directory: /trunk/
==================
File [modified]: interface.ccl
Delta lines: +14 -0
===================================================================
--- trunk/interface.ccl 2013-06-08 00:23:17 UTC (rev 83)
+++ trunk/interface.ccl 2013-06-13 23:10:44 UTC (rev 84)
@@ -22,7 +22,21 @@
PROVIDES FUNCTION EOS_Omni_press WITH EOS_Omni_EOS_Press LANGUAGE Fortran
+void FUNCTION EOS_Omni_pressOMP(CCTK_INT IN eoskey, \
+ CCTK_INT IN havetemp, \
+ CCTK_REAL IN rf_precision, \
+ 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_pressOMP WITH EOS_Omni_EOS_PressOMP LANGUAGE Fortran
+
+
void FUNCTION EOS_Omni_DPressByDEps(CCTK_INT IN eoskey, \
CCTK_INT IN havetemp, \
CCTK_REAL IN rf_precision, \
More information about the Commits
mailing list