Commits March 2014

commits@lists.einsteintoolkit.org

8 participants
172 discussions

[ianhinder/Kranc] 7a62fd: Handle namespaces correctly for DGFE
by Erik Schnetter 13 Mar '14

13 Mar '14

Branch: refs/heads/master Home: https://github.com/ianhinder/Kranc Commit: 7a62fd9774e25ff383d1cc493cce79adfb387b11 https://github.com/ianhinder/Kranc/commit/7a62fd9774e25ff383d1cc493cce79adf… Author: Erik Schnetter <schnetter(a)gmail.com> Date: 2014-03-13 (Thu, 13 Mar 2014) Changed paths: M Tools/CodeGen/CodeGenCalculation.m M Tools/CodeGen/DGFE.m Log Message: ----------- Handle namespaces correctly for DGFE #include statements must occur outside namespaces. Template instantiations require switching forth and back between namespaces.

1 0

[McLachlan]McLachlan branch, dgfe2, updated. fd6af7e85e9b1d2307923e4b06fcfd6d4a12b1a5
by Carpet git repositories 13 Mar '14

13 Mar '14

This is an automated email from the git hooks/post-receive script. It was generated because a ref change was pushed to the repository containing the project "McLachlan". The branch, dgfe2 has been updated via fd6af7e85e9b1d2307923e4b06fcfd6d4a12b1a5 (commit) via 299573595547c5dbe886e38438038266ee710406 (commit) via 0607bf6e17ee84c1da9da751062a007e7da8d43d (commit) from 75a09aafd77e08d5dc809fbfdb627dfc82e120a7 (commit) Those revisions listed above that are new to this repository have not appeared on any other notification email; so we list those revisions in full, below. - Log ----------------------------------------------------------------- commit fd6af7e85e9b1d2307923e4b06fcfd6d4a12b1a5 Author: Erik Schnetter <schnetter(a)gmail.com> Date: Thu Mar 13 10:53:46 2014 -0400 Add six-patch DG examples commit 299573595547c5dbe886e38438038266ee710406 Author: Erik Schnetter <schnetter(a)gmail.com> Date: Thu Mar 13 10:52:35 2014 -0400 Regenerate code commit 0607bf6e17ee84c1da9da751062a007e7da8d43d Author: Erik Schnetter <schnetter(a)gmail.com> Date: Thu Mar 13 10:41:24 2014 -0400 Re-calculate radius r after changing coordinates ----------------------------------------------------------------------- Summary of changes: ML_Test_DG4/src/ML_Coordinates.cc | 7 ++ ML_Test_FD4/src/ML_Coordinates.cc | 7 ++ ML_WaveToy_DG4/src/ML_Coordinates.cc | 7 ++ ML_WaveToy_DG8/src/ML_Coordinates.cc | 7 ++ ML_WaveToy_FD4/src/ML_Coordinates.cc | 7 ++ ML_WaveToy_FD8/src/ML_Coordinates.cc | 7 ++ m/Test.merge/src/ML_Coordinates.cc | 7 ++ m/WaveToy.merge/src/ML_Coordinates.cc | 7 ++ ...o_dg4_sinex.par => wavetoy_fo_dg4_gaussian.par} | 22 +++--- ....par => wavetoy_fo_dg4_sixpatches_gaussian.par} | 73 ++++++++++++-------- ...nex.par => wavetoy_fo_dg4_sixpatches_sinex.par} | 70 +++++++++++-------- 11 files changed, 156 insertions(+), 65 deletions(-) copy par/{wavetoy_fo_dg4_sinex.par => wavetoy_fo_dg4_gaussian.par} (86%) copy par/{wavetoy_fo_dg4_distorted_sinex.par => wavetoy_fo_dg4_sixpatches_gaussian.par} (69%) copy par/{wavetoy_fo_dg4_distorted_sinex.par => wavetoy_fo_dg4_sixpatches_sinex.par} (70%) hooks/post-receive -- McLachlan

1 0

[svn:einsteintoolkit] EOS_Omni/trunk/doc/ (Rev. 106)
by rhaas＠tapir.caltech.edu 12 Mar '14

12 Mar '14

User: rhaas Date: 2014/03/12 10:02 PM Modified: /trunk/doc/ documentation.tex Log: EOS_Omni: update documentation From: Roland Haas <rhaas(a)tapir.caltech.edu> File Changes: Directory: /trunk/doc/ ====================== File [modified]: documentation.tex Delta lines: +16 -3 =================================================================== --- trunk/doc/documentation.tex 2014-03-13 03:02:09 UTC (rev 105) +++ trunk/doc/documentation.tex 2014-03-13 03:02:11 UTC (rev 106) @@ -113,9 +113,10 @@ interface and implements multiple analytic EOS, and also provides table reader and interpolation routines for finite-temperature microphysical EOS available from - \url{http://www.stellarcollapse.org}. Currently, the implemented - analytic EOS are the polytropic EOS, the gamma-law EOS, and a hybrid EOS - consisting of a 2-piece piecewise-polytrope with an a thermal, + \url{http://www.stellarcollapse.org}. In addition cold and barotropic + tabulated EOS are provided for use in initial data thorns. Currently, the + implemented analytic EOS are the polytropic EOS, the gamma-law EOS, and a + hybrid EOS consisting of a 2-piece piecewise-polytrope with a thermal, gamma-law component. \end{abstract} @@ -173,6 +174,7 @@ for stellar core collapse simulations. \item \texttt{eoskey = 4}: Finite-temperature microphysical EOS \item \texttt{eoskey = 5}: Cold tabulated EOS with Gamma-Law + \item \texttt{eoskey = 6}: Tabulated barotropic EOS \end{itemize} \item \texttt{havetemp} determines whether the EOS is to be called as a function of $(\rho,\epsilon,Y_e)$ (\texttt{havetemp = 0}), or as a @@ -211,9 +213,20 @@ \begin{verbatim} eos_omni::poly_gamma +eos_omni::poly_gamma_initial eos_omni::poly_k \end{verbatim} +The only non-obvious parameter here is \texttt{poly\_gamma\_initial}. +In most simulations it should not be set. +In simulations that are run with a different adiabatic index than +what was used to set up the initial data, \texttt{poly\_gamma} should +be the evolution value, and \texttt{poly\_gamma\_initial} should be the +initial data value. \texttt{EOS\_Omni} then rescales \texttt{poly\_k} such +that the \emph{cgs} value of \texttt{poly\_k} is the same for initial data and +evolution. Since the units of \texttt{poly\_k} +($[K] = \frac{\rm N}{\rm m^2} \left(\frac{\rm m^3}{\rm kg}\right)^{\Gamma}$) +depend on \texttt{poly\_gamma} this is not fully trivial. Check \texttt{param.ccl} for parameters for the other EOS\@.

1 0

[svn:einsteintoolkit] EOS_Omni/trunk/src/nuc_eos_cxx/ (Rev. 105)
by rhaas＠tapir.caltech.edu 12 Mar '14

12 Mar '14

User: rhaas Date: 2014/03/12 10:02 PM Modified: /trunk/src/nuc_eos_cxx/ nuc_eos_press_cs2.cc Log: EOS_Omni nuc_eos_cxx: ensure cs2 is 0 or positive From: Christian Ott <cott(a)tapir.caltech.edu> File Changes: Directory: /trunk/src/nuc_eos_cxx/ ================================== File [modified]: nuc_eos_press_cs2.cc Delta lines: +8 -2 =================================================================== --- trunk/src/nuc_eos_cxx/nuc_eos_press_cs2.cc 2014-03-13 03:02:05 UTC (rev 104) +++ trunk/src/nuc_eos_cxx/nuc_eos_press_cs2.cc 2014-03-13 03:02:09 UTC (rev 105) @@ -63,7 +63,10 @@ prs[i] = exp(prs[i]); eps[i] = exp(eps[i]) - energy_shift; #if HAVEGR - cs2[i] = rho[i] * cs2[i] / (rho[i] + rho[i] * eps[i] + prs[i]); + cs2[i] = MAX(rho[i] * cs2[i] / (rho[i] + rho[i] * eps[i] + prs[i]), + 0.0e0); +#else + cs2[i] = MAX(cs2[i], 0.0e0); #endif } @@ -175,7 +178,10 @@ for(int i=0;i<n;i++) { prs[i] = exp(prs[i]); #if HAVEGR - cs2[i] = rho[i] * cs2[i] / (rho[i] + rho[i]*eps[i] + prs[i]); + cs2[i] = MAX(rho[i] * cs2[i] / (rho[i] + rho[i] * eps[i] + prs[i]), + 0.0e0); +#else + cs2[i] = MAX(cs2[i], 0.0e0); #endif }

1 0

[svn:einsteintoolkit] EOS_Omni/trunk/ (Rev. 104)
by rhaas＠tapir.caltech.edu 12 Mar '14

12 Mar '14

User: rhaas Date: 2014/03/12 10:02 PM Added: /trunk/src/nuc_eos_cxx/ dumpASCIItable.cc Modified: /trunk/ schedule.ccl /trunk/src/nuc_eos_cxx/ make.code.defn /trunk/test/ test.ccl Log: EOS_Omni: add code for test case to C++ version From: Roland Haas <rhaas(a)tapir.caltech.edu> File Changes: Directory: /trunk/src/nuc_eos_cxx/ ================================== File [added]: dumpASCIItable.cc Delta lines: +110 -0 =================================================================== --- trunk/src/nuc_eos_cxx/dumpASCIItable.cc (rev 0) +++ trunk/src/nuc_eos_cxx/dumpASCIItable.cc 2014-03-13 03:02:05 UTC (rev 104) @@ -0,0 +1,110 @@ +#include "cctk.h" +#include "cctk_Parameters.h" +#include "cctk_Arguments.h" +#include "cctk_Functions.h" + +#include <string> +#include <cstdio> +#include <cassert> + +#include "nuc_eos.hh" + +#define DIM(v) (sizeof(v)/sizeof((v)[0])) + +extern "C" +void EOS_OMNI_dumptable(CCTK_ARGUMENTS) +{ + using namespace nuc_eos; + using namespace nuc_eos_private; + + DECLARE_CCTK_PARAMETERS; + DECLARE_CCTK_ARGUMENTS; + + int irho,itemp,iye,n; + int strlength1,strlength2; + std::string fullpath; + + if (CCTK_MyProc(cctkGH) != 0) + return; + + fullpath = std::string(out_dir)+"/"+std::string(dump_nuceos_table_name); + + CCTK_INFO("*******************************"); + CCTK_INFO("Dumping nuc_eos table file in ASCII:"); + CCTK_INFO(fullpath.c_str()); + CCTK_INFO("*******************************"); + + FILE *fh = fopen(fullpath.c_str(), "w"); + assert(fh); + + fprintf(fh, "# % 20s\n% 4d\n", "nrho:",nrho); + fprintf(fh, "# % 20s\n% 4d\n", "ntemp:",ntemp); + fprintf(fh, "# % 20s\n% 4d\n", "nye:",nye); + + fprintf(fh, "# % 20s\n%18.9E\n", "energy shift:",energy_shift); + + fprintf(fh, "# % 20s\n%18.9E %18.9E\n", "rho min and max:",eos_rhomin/RHOGF,eos_rhomax/RHOGF); + fprintf(fh, "# % 20s\n%18.9E %18.9E\n", "ye min and max:",eos_yemin,eos_yemax); + fprintf(fh, "# % 20s\n%18.9E %18.9E\n", "temp min and max:",eos_tempmin,eos_tempmax); + + fprintf(fh, "# % 20s\n\n", "log rho points:"); + for(irho=0;irho<nrho;++irho) { + fprintf(fh, "%18.9E\n", (logrho[irho] - log(RHOGF))/log(10.)); + } + fprintf(fh, "# % 20s\n", "log temp points:"); + for(itemp=0;itemp<ntemp;++itemp) { + fprintf(fh, "%18.9E\n", logtemp[itemp]/log(10.)); + } + fprintf(fh, "# % 20s\n", "ye points:"); + for(iye=0;iye<nye;++iye) { + fprintf(fh, "%18.9E\n", yes[iye]); + } + + fprintf(fh, "# % 20s\n% 4d\n", "nvars:",NTABLES); + fprintf(fh, "# % 20s\n", "table mappings:"); + fprintf(fh, "# % 20s\n", " 1 -> logpress"); + fprintf(fh, "# % 20s\n", " 2 -> logenergy"); + fprintf(fh, "# % 20s\n", " 3 -> entropy"); + fprintf(fh, "# % 20s\n", " 4 -> munu"); + fprintf(fh, "# % 20s\n", " 5 -> cs2"); + fprintf(fh, "# % 20s\n", " 6 -> dedT"); + fprintf(fh, "# % 20s\n", " 7 -> dpdrhoe"); + fprintf(fh, "# % 20s\n", " 8 -> dpderho"); + fprintf(fh, "# % 20s\n", " 9 -> muhat"); + fprintf(fh, "# % 20s\n", "10 -> mu_e"); + fprintf(fh, "# % 20s\n", "11 -> mu_p"); + fprintf(fh, "# % 20s\n", "12 -> mu_n"); + fprintf(fh, "# % 20s\n", "13 -> xa"); + fprintf(fh, "# % 20s\n", "14 -> xh"); + fprintf(fh, "# % 20s\n", "15 -> xn"); + fprintf(fh, "# % 20s\n", "16 -> xp"); + fprintf(fh, "# % 20s\n", "17 -> abar"); + fprintf(fh, "# % 20s\n", "18 -> zbar"); + fprintf(fh, "# % 20s\n", "19 -> gamma"); + + const double ctable[][2] = { // un-convert units + { log(10.0), log(PRESSGF) }, // pressure + { log(10.0), log(EPSGF) }, // eps + { 1. }, + { 1. }, + { LENGTHGF*LENGTHGF/TIMEGF/TIMEGF }, // cs2 + { EPSGF }, // dedT + { PRESSGF/RHOGF }, // dpdrhoe + { PRESSGF/EPSGF },// dpderho + }; + for(irho=0;irho<nrho;++irho) { + for(itemp=0;itemp<ntemp;++itemp) { + for(iye=0;iye<nye;++iye) { + for(n=0;n<NTABLES;++n) { + // transpose nvar to front as we go along + const int idx = n + NTABLES*(irho + nrho*(itemp + ntemp*iye)); + double val = alltables[idx]; + if(n < DIM(ctable)) val = (val - ctable[n][1]) / ctable[n][0]; + fprintf(fh, "% 4d % 4d % 4d % 4d %18.9E\n", irho+1,itemp+1,iye+1,n+1,val); + } + } + } + } + + fclose(fh); +} File [modified]: make.code.defn Delta lines: +1 -1 =================================================================== --- trunk/src/nuc_eos_cxx/make.code.defn 2014-03-13 03:02:01 UTC (rev 103) +++ trunk/src/nuc_eos_cxx/make.code.defn 2014-03-13 03:02:05 UTC (rev 104) @@ -1,4 +1,4 @@ SRCS = readtable.cc nuc_eos_short.cc nuc_eos_press_cs2.cc \ nuc_eos_press.cc nuc_eos_full.cc nuc_eos_dpdrhoe_dpderho.cc \ - readtable_cactus_wrapper.cc + readtable_cactus_wrapper.cc dumpASCIItable.cc SUBDIRS = Directory: /trunk/test/ ======================= File [modified]: test.ccl Delta lines: +2 -1 =================================================================== --- trunk/test/test.ccl 2014-03-13 03:02:01 UTC (rev 103) +++ trunk/test/test.ccl 2014-03-13 03:02:05 UTC (rev 104) @@ -1 +1,2 @@ -RELTOL 1e-14 +# we only output 9 digits +RELTOL 1e-8 Directory: /trunk/ ================== File [modified]: schedule.ccl Delta lines: +8 -0 =================================================================== --- trunk/schedule.ccl 2014-03-13 03:02:01 UTC (rev 103) +++ trunk/schedule.ccl 2014-03-13 03:02:05 UTC (rev 104) @@ -21,6 +21,14 @@ OPTIONS: global } "setup energy_shift in EOS_Omni_Module" + if (dump_nuceos_table) + { + SCHEDULE EOS_OMNI_dumptable after nuc_eos_readtable_cactus_wrapper AT CCTK_BASEGRID + { + LANG: C + OPTIONS: global + } "Dump EOS HDF5 table in ASCII" + } } if (coldeos_read_table)

1 0

[svn:einsteintoolkit] EOS_Omni/trunk/ (Rev. 103)
by rhaas＠tapir.caltech.edu 12 Mar '14

12 Mar '14

User: rhaas Date: 2014/03/12 10:02 PM Modified: /trunk/ interface.ccl /trunk/src/ EOS_Omni_MultiVarCalls.F90 Log: make new EOS call that returns press and cs2 [and eps if needed] From: Christian D. Ott <cott(a)tapir.caltech.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: EOS_Omni_MultiVarCalls.F90 Delta lines: +177 -8 =================================================================== --- trunk/src/EOS_Omni_MultiVarCalls.F90 2014-03-13 03:01:59 UTC (rev 102) +++ trunk/src/EOS_Omni_MultiVarCalls.F90 2014-03-13 03:02:01 UTC (rev 103) @@ -101,7 +101,6 @@ CCTK_REAL, intent(out) :: muhat(npoints) ! local vars - integer :: i character(256) :: warnstring if(eoskey.ne.4) then @@ -157,10 +156,6 @@ real*8 :: hybrid_local_gamma real*8 :: hybrid_local_k_cgs real*8 :: hybrid_dp_poly,hybrid_dp_th1,hybrid_dp_th2 - ! temporary vars for nuc_eos - real*8 :: xrho,xye,xtemp,xenr,xent - real*8 :: xprs,xmunu,xcs2 - real*8 :: xdedt,xdpderho,xdpdrhoe anyerr = 0 keyerr(:) = 0 @@ -264,9 +259,7 @@ real*8 :: hybrid_local_k_cgs real*8 :: hybrid_dp_poly,hybrid_dp_th1,hybrid_dp_th2 ! temporary vars for nuc_eos - real*8 :: xrho,xye,xtemp,xenr,xent - real*8 :: xprs,xmunu,xcs2 - real*8 :: xdedt,xdpderho,xdpdrhoe + real*8 :: xdpderho,xdpdrhoe anyerr = 0 @@ -342,3 +335,179 @@ end select end subroutine EOS_Omni_EOS_DEpsByDRho_DEpsByDPress + + +subroutine EOS_Omni_EOS_Press_cs2(eoskey,keytemp,rf_precision,npoints,& + rho,eps,temp,ye,press,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) :: rf_precision + CCTK_REAL, intent(in) :: rho(npoints),ye(npoints) + CCTK_REAL, intent(inout) :: eps(npoints), temp(npoints) + CCTK_REAL, intent(out) :: press(npoints) + CCTK_REAL, intent(out) :: cs2(npoints) + + ! local vars + 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 + ! temporary vars for nuc_eos + real*8 :: xrho + real*8 :: 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, cs2_cold, cs2_th, h, h_cold + + 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 + cs2(i) = poly_gamma * press(i) / rho(i) / & + (1 + eps(i) + press(i)/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 + press(i) = (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 + press(i) * xdpderho / (rho(i)**2)) / & + (1.0d0 + eps(i) + press(i)/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 + 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 + + cs2(i) = (hybrid_local_gamma * hybrid_p_poly + hybrid_gamma_th * hybrid_p_th) / & + rho(i) / (1.0d0 + eps(i) + press(i)/rho(i)) + enddo + + case (4) + ! nuc eos + if(keytemp.eq.1) then + call nuc_eos_m_kt1_press_eps_cs2(npoints,rho,temp,ye,& + eps,press,cs2,keyerr,anyerr) + else if(keytemp.eq.0) then + call nuc_eos_m_kt0_press_cs2(npoints,rho,temp,ye,eps,press,& + cs2,rf_precision,keyerr,anyerr) + else + call CCTK_ERROR("This keytemp is not suppported!") + STOP + endif + 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) + 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) + cs2(i) = coldeos_low_kappa * coldeos_low_gamma * & + rho(i)**(coldeos_low_gamma-1.0d0) / & + (1.0d0 + eps(i) + press(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) + + ! this is the cold speed of sound squared + cs2_cold = (coldeos_cs2(ir) - coldeos_cs2(ir-1)) / & + (coldeos_cs2(ir) - coldeos_cs2(ir-1)) * & + (xrho - coldeos_logrho(ir-1)) + coldeos_cs2(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 + + xdpdrhoe = coldeos_thfac*(coldeos_gammath - 1.0d0)*eps_th + xdpderho = coldeos_thfac*(coldeos_gammath - 1.0d0)*rho(i) + cs2_th = (xdpdrhoe + press_th * xdpderho / (rho(i)**2)) + + h = 1.0d0 + eps(i) + (press_cold+press_th) / rho(i) + h_cold = 1.0d0 + eps_cold + press_cold / rho(i) + + cs2(i) = (cs2_cold * h_cold + cs2_th) / h + + enddo + + case (6) + ! barotropic tabular EOS with gamma law + write(warnstring,*) "eoskey ",eoskey," for barotropic EOS not implemented!" + call CCTK_ERROR(warnstring) + STOP + case DEFAULT + write(warnstring,*) "eoskey ",eoskey," not implemented!" + call CCTK_ERROR(warnstring) + STOP + end select + +end subroutine EOS_Omni_EOS_Press_cs2 + Directory: /trunk/ ================== File [modified]: interface.ccl Delta lines: +15 -0 =================================================================== --- trunk/interface.ccl 2014-03-13 03:01:59 UTC (rev 102) +++ trunk/interface.ccl 2014-03-13 03:02:01 UTC (rev 103) @@ -22,6 +22,21 @@ PROVIDES FUNCTION EOS_Omni_press WITH EOS_Omni_EOS_Press LANGUAGE Fortran +void FUNCTION EOS_Omni_press_cs2(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_REAL OUT ARRAY cs2, \ + CCTK_INT OUT ARRAY keyerr, \ + CCTK_INT OUT anyerr) + +PROVIDES FUNCTION EOS_Omni_press_cs2 WITH EOS_Omni_EOS_Press_cs2 LANGUAGE Fortran + void FUNCTION EOS_Omni_pressOMP(CCTK_INT IN eoskey, \ CCTK_INT IN havetemp, \ CCTK_REAL IN rf_precision, \

1 0

[svn:einsteintoolkit] EOS_Omni/trunk/src/nuc_eos_cxx/ (Rev. 102)
by rhaas＠tapir.caltech.edu 12 Mar '14

12 Mar '14

User: rhaas Date: 2014/03/12 10:01 PM Modified: /trunk/src/nuc_eos_cxx/ helpers.hh, nuc_eos.hh, readtable.cc Log: EOS_Omni: move actual table data into namespace of its own From: Roland Haas <rhaas(a)tapir.caltech.edu> File Changes: Directory: /trunk/src/nuc_eos_cxx/ ================================== File [modified]: helpers.hh Delta lines: +9 -0 =================================================================== --- trunk/src/nuc_eos_cxx/helpers.hh 2014-03-13 03:01:55 UTC (rev 101) +++ trunk/src/nuc_eos_cxx/helpers.hh 2014-03-13 03:01:59 UTC (rev 102) @@ -79,6 +79,7 @@ int* restrict idx) { using namespace nuc_eos; + using namespace nuc_eos_private; int ix = 1 + (int)( (x - logrho[0] - 1.0e-10) * drhoi ); int iy = 1 + (int)( (y - logtemp[0] - 1.0e-10) * dtempi ); @@ -115,6 +116,7 @@ int* restrict idx) { using namespace nuc_eos; + using namespace nuc_eos_private; int ix = 1 + (int)( (x - logrho[0] - 1.0e-10) * drhoi ); int iy = 1; @@ -150,6 +152,7 @@ int* restrict idx) { using namespace nuc_eos; + using namespace nuc_eos_private; int ix = 1 + (int)( (x - logrho[0] - 1.0e-10) * drhoi ); int iy = 1; @@ -185,6 +188,7 @@ const int iv) { using namespace nuc_eos; + using namespace nuc_eos_private; // helper variables double fh[8], a[8]; @@ -230,6 +234,7 @@ const int iv) { using namespace nuc_eos; + using namespace nuc_eos_private; // helper variables double fh[8], a[8]; @@ -276,6 +281,7 @@ { using namespace nuc_eos; + using namespace nuc_eos_private; // helper variables double fh[8], a[8]; @@ -351,6 +357,7 @@ // iv is the index of the variable we do the bisection on using namespace nuc_eos; + using namespace nuc_eos_private; int bcount = 0; int maxbcount = 80; @@ -464,6 +471,7 @@ int *keyerrt) { using namespace nuc_eos; + using namespace nuc_eos_private; // local variables const int itmax = 200; // use at most 10 iterations, then go to bisection @@ -617,6 +625,7 @@ int *keyerrt) { using namespace nuc_eos; + using namespace nuc_eos_private; // local variables const int itmax = 200; // use at most 10 iterations, then go to bisection File [modified]: nuc_eos.hh Delta lines: +27 -20 =================================================================== --- trunk/src/nuc_eos_cxx/nuc_eos.hh 2014-03-13 03:01:55 UTC (rev 101) +++ trunk/src/nuc_eos_cxx/nuc_eos.hh 2014-03-13 03:01:59 UTC (rev 102) @@ -20,28 +20,8 @@ namespace nuc_eos { - extern int nrho; - extern int ntemp; - extern int nye; - - extern double * restrict alltables; - extern double * restrict epstable; - extern double * restrict logrho; - extern double * restrict logtemp; extern double temp0, temp1; - extern double dlintemp,dlintempi; - extern double drholintempi; - extern double dlintempyei; - extern double drholintempyei; - extern double * restrict yes; extern double energy_shift; - extern double dtemp, dtempi; - extern double drho, drhoi; - extern double dye, dyei; - extern double drhotempi; - extern double drhoyei; - extern double dtempyei; - extern double drhotempyei; // min and max values @@ -77,4 +57,31 @@ i_Xh, i_Xn, i_Xp, i_Abar, i_Zbar, i_Gamma}; } +namespace nuc_eos_private { + +// table data + + extern int nrho; + extern int ntemp; + extern int nye; + + extern double * restrict alltables; + extern double * restrict epstable; + extern double * restrict logrho; + extern double * restrict logtemp; + extern double dlintemp,dlintempi; + extern double drholintempi; + extern double dlintempyei; + extern double drholintempyei; + extern double * restrict yes; + extern double dtemp, dtempi; + extern double drho, drhoi; + extern double dye, dyei; + extern double drhotempi; + extern double drhoyei; + extern double dtempyei; + extern double drhotempyei; + +} + #endif // NUC_EOS_HH File [modified]: readtable.cc Delta lines: +13 -10 =================================================================== --- trunk/src/nuc_eos_cxx/readtable.cc 2014-03-13 03:01:55 UTC (rev 101) +++ trunk/src/nuc_eos_cxx/readtable.cc 2014-03-13 03:01:59 UTC (rev 102) @@ -32,6 +32,18 @@ // define the variables namespace nuc_eos { + double temp0, temp1; + double energy_shift; + + double eos_rhomax, eos_rhomin; + double eos_tempmin, eos_tempmax; + double eos_yemin, eos_yemax; + + double c2p_tempmin; + double c2p_tempmax; + +} +namespace nuc_eos_private { int nrho; int ntemp; int nye; @@ -40,13 +52,11 @@ double * restrict epstable; double * restrict logrho; double * restrict logtemp; - double temp0, temp1; double dlintemp, dlintempi; double drholintempi; double dlintempyei; double drholintempyei; double * restrict yes; - double energy_shift; double dtemp, dtempi; double drho, drhoi; double dye, dyei; @@ -54,14 +64,6 @@ double drhoyei; double dtempyei; double drhotempyei; - - double eos_rhomax, eos_rhomin; - double eos_tempmin, eos_tempmax; - double eos_yemin, eos_yemax; - - double c2p_tempmin; - double c2p_tempmax; - } @@ -75,6 +77,7 @@ void nuc_eos_C_ReadTable(char* nuceos_table_name) { using namespace nuc_eos; + using namespace nuc_eos_private; CCTK_VInfo(CCTK_THORNSTRING,"*******************************"); CCTK_VInfo(CCTK_THORNSTRING,"Reading nuc_eos table file:");

1 0

[svn:einsteintoolkit] EOS_Omni/trunk/src/nuc_eos_cxx/ (Rev. 101)
by rhaas＠tapir.caltech.edu 12 Mar '14

12 Mar '14

User: rhaas Date: 2014/03/12 10:01 PM Modified: /trunk/src/nuc_eos_cxx/ helpers.hh, nuc_eos_dpdrhoe_dpderho.cc, nuc_eos_full.cc, nuc_eos_press.cc, nuc_eos_press_cs2.cc, nuc_eos_short.cc, readtable.cc Log: EOS_Omni: low-level optimizations * use log-rules to transform base 10 logs to natural ones faster and hopefully more accurate * optimize away a division when computing cs2 * optimize bracketing test to (a-a1)*(a-a2)<0 * replace divisions by multiplication by inverse * use CCTK_BUILTIN_EXPECT to indicate likely outcome this should help branch prediction since the compiler can put the the less likely branch into the unfavorable location. * remove superfluous if statements whenever we leave a itmax loop regularly we iterated for too long * explicitly replace x/y by x*1./y From: Roland Haas <rhaas(a)tapir.caltech.edu> File Changes: Directory: /trunk/src/nuc_eos_cxx/ ================================== File [modified]: helpers.hh Delta lines: +49 -52 =================================================================== --- trunk/src/nuc_eos_cxx/helpers.hh 2014-03-13 03:01:52 UTC (rev 100) +++ trunk/src/nuc_eos_cxx/helpers.hh 2014-03-13 03:01:55 UTC (rev 101) @@ -18,24 +18,24 @@ // 105 -- rho too high // 106 -- rho too low - if(xrho > eos_rhomax) { + if(CCTK_BUILTIN_EXPECT(xrho > eos_rhomax, false)) { return 105; } - if(xrho < eos_rhomin) { + if(CCTK_BUILTIN_EXPECT(xrho < eos_rhomin, false)) { return 106; } - if(xye > eos_yemax) { + if(CCTK_BUILTIN_EXPECT(xye > eos_yemax, false)) { return 101; } - if(xye < eos_yemin) { + if(CCTK_BUILTIN_EXPECT(xye < eos_yemin, false)) { // this is probably not pure and should be removed fprintf(stderr,"xye: %15.6E eos_yemin: %15.6E\n",xye,eos_yemin); return 102; } - if(xtemp > eos_tempmax) { + if(CCTK_BUILTIN_EXPECT(xtemp > eos_tempmax, false)) { return 103; } - if(xtemp < eos_tempmin) { + if(CCTK_BUILTIN_EXPECT(xtemp < eos_tempmin, false)) { return 104; } return 0; @@ -53,16 +53,16 @@ // 105 -- rho too high // 106 -- rho too low - if(xrho > eos_rhomax) { + if(CCTK_BUILTIN_EXPECT(xrho > eos_rhomax, false)) { return 105; } - if(xrho < eos_rhomin) { + if(CCTK_BUILTIN_EXPECT(xrho < eos_rhomin, false)) { return 106; } - if(xye > eos_yemax) { + if(CCTK_BUILTIN_EXPECT(xye > eos_yemax, false)) { return 101; } - if(xye < eos_yemin) { + if(CCTK_BUILTIN_EXPECT(xye < eos_yemin, false)) { return 102; } return 0; @@ -331,10 +331,12 @@ // then interpolate in y // assume rectangular grid - double t1 = (fs[1]-fs[0])/(xs[1]-xs[0]) * (x - xs[0]) + fs[0]; - double t2 = (fs[3]-fs[2])/(xs[1]-xs[0]) * (x - xs[0]) + fs[2]; + double dxi = 1./(xs[1]-xs[0]); + double dyi = 1./(ys[1]-ys[0]); // x*1./y uses faster instructions than x/y + double t1 = (fs[1]-fs[0])*dxi * (x - xs[0]) + fs[0]; + double t2 = (fs[3]-fs[2])*dxi * (x - xs[0]) + fs[2]; - return (t2 - t1)/(ys[1]-ys[0]) * (y-ys[0]) + t1; + return (t2 - t1)*dyi * (y-ys[0]) + t1; } static inline __attribute__((always_inline)) @@ -359,6 +361,8 @@ const double dltp = log(1.2); const double dltm = log(0.8); + const double leps0_prec = leps0*prec; + // temporary local vars double lt, lt1, lt2; double ltmin = logtemp[0]; @@ -404,7 +408,7 @@ #endif bcount++; - if(bcount >= maxbcount) { + if(CCTK_BUILTIN_EXPECT(bcount >= maxbcount, false)) { *keyerrt = 667; return; } @@ -432,22 +436,20 @@ nuc_eos_C_linterp_one(idx,delx,dely,delz,&f2a,iv); fmid=f2a-leps0; - if(fmid <= 0.0) lt=ltmid; + if(CCTK_BUILTIN_EXPECT(fmid <= 0.0, false)) lt=ltmid; #if DEBUG fprintf(stderr,"bisection step 2 it %d, fmid: %15.6E f2a: %15.6E lt: %15.6E ltmid: %15.6E dlt: %15.6E\n", it,fmid,f2a,lt,ltmid,dlt); #endif - if(fabs(1.0-f2a/leps0) <= prec) { + if(CCTK_BUILTIN_EXPECT(fabs(leps0-f2a) <= leps0_prec, false)) { *ltout = ltmid; return; } } // for it = 0 - if(it >= itmax-1) { - *keyerrt = 667; - return; - } + *keyerrt = 667; + return; } // bisection @@ -465,7 +467,7 @@ // local variables const int itmax = 200; // use at most 10 iterations, then go to bisection - double dlepsdlt; // derivative dlogeps/dlogT + double dlepsdlti; // 1 / derivative dlogeps/dlogT double ldt; double leps,leps0; // temp vars for eps double ltn, lt; // temp vars for temperature @@ -489,6 +491,7 @@ fprintf(stderr,"it: %d t: %15.6E leps: %15.6E eps0: %15.6E del: %15.6E\n", it,lt, leps,leps0,fabs(leps-leps0)/(fabs(leps0))); #endif + // TODO: profile this to see which outcome is more likely if(fabs(leps-leps0) < prec*fabs(leps0)) { *ltout = lt0; return; @@ -548,8 +551,7 @@ // Check if we are already bracketing the input internal // energy. If so, interpolate for new T. - if( (leps0 >= epst1 && leps0 <= epst2) - || (leps0 <= epst1 && leps0 >=epst2)) { + if(CCTK_BUILTIN_EXPECT((leps0 - epst1) * (leps0 - epst2) <= 0., false)) { *ltout = (logtemp[itemp]-logtemp[itemp-1]) / (epst2 - epst1) * (leps0 - epst1) + logtemp[itemp-1]; @@ -559,8 +561,8 @@ // well, then do a Newton-Raphson step // first, guess the derivative - dlepsdlt = (epst2-epst1)/(logtemp[itemp]-logtemp[itemp-1]); - ldt = -(leps - leps0) / dlepsdlt * fac; + dlepsdlti = (logtemp[itemp]-logtemp[itemp-1])/(epst2-epst1); + ldt = -(leps - leps0) * dlepsdlti * fac; ltn = MIN(MAX(lt + ldt,ltmin),ltmax); lt = ltn; @@ -578,7 +580,7 @@ if(oerr < err) fac *= 0.9; oerr = err; - if(err < prec*fabs(leps0)) { + if(CCTK_BUILTIN_EXPECT(err < prec*fabs(leps0), false)) { *ltout = lt; return; } @@ -587,20 +589,18 @@ } // while(it < itmax) - if(it >= itmax - 1) { - // try bisection - // bisection(lr, lt0, ye, leps0, ltout, 1, prec, keyerrt); + // try bisection + // bisection(lr, lt0, ye, leps0, ltout, 1, prec, keyerrt); #if DEBUG - fprintf(stderr, "Failed to converge. This is bad. Trying bisection!\n"); + fprintf(stderr, "Failed to converge. This is bad. Trying bisection!\n"); #endif - bisection(lr,lt0,ye,leps0,prec,ltout,1,keyerrt); + bisection(lr,lt0,ye,leps0,prec,ltout,1,keyerrt); #if DEBUG - if(*keyerrt==667) { - fprintf(stderr,"This is worse. Bisection failed!\n"); - abort(); - } + if(*keyerrt==667) { + fprintf(stderr,"This is worse. Bisection failed!\n"); + abort(); + } #endif - } return; @@ -620,7 +620,7 @@ // local variables const int itmax = 200; // use at most 10 iterations, then go to bisection - double dentdlt; // derivative dentropy/dlogT + double dentdlti; // 1 / derivative dentropy/dlogT double ldt; double ent,ent0; // temp vars for eps double ltn, lt; // temp vars for temperature @@ -644,7 +644,7 @@ fprintf(stderr,"it: %d t: %15.6E leps: %15.6E eps0: %15.6E del: %15.6E\n", it,lt,ent,ent0,fabs(ent-ent0)/(fabs(ent0))); #endif - if(fabs(ent-ent0) < prec*fabs(ent0)) { + if(CCTK_BUILTIN_EXPECT(fabs(ent-ent0) < prec*fabs(ent0), false)) { *ltout = lt0; return; } @@ -703,8 +703,7 @@ // Check if we are already bracketing the input internal // energy. If so, interpolate for new T. - if( (ent0 >= ent1 && ent0 <= ent2) - || (ent0 <= ent1 && ent0 >= ent2)) { + if(CCTK_BUILTIN_EXPECT((ent0 - ent1) * (ent0 - ent2) <= 0., false)) { *ltout = (logtemp[itemp]-logtemp[itemp-1]) / (ent2 - ent1) * (ent0 - ent1) + logtemp[itemp-1]; @@ -714,8 +713,8 @@ // well, then do a Newton-Raphson step // first, guess the derivative - dentdlt = (ent2-ent1)/(logtemp[itemp]-logtemp[itemp-1]); - ldt = -(ent - ent0) / dentdlt * fac; + dentdlti = (logtemp[itemp]-logtemp[itemp-1])/(ent2-ent1); + ldt = -(ent - ent0) * dentdlti * fac; ltn = MIN(MAX(lt + ldt,ltmin),ltmax); lt = ltn; @@ -733,7 +732,7 @@ if(oerr < err) fac *= 0.9; oerr = err; - if(err < prec*fabs(ent0)) { + if(CCTK_BUILTIN_EXPECT(err < prec*fabs(ent0), false)) { *ltout = lt; return; } @@ -742,19 +741,17 @@ } // while(it < itmax) - if(it >= itmax - 1) { - // try bisection + // try bisection #if DEBUG - fprintf(stderr, "Failed to converge. This is bad. Trying bisection!\n"); + fprintf(stderr, "Failed to converge. This is bad. Trying bisection!\n"); #endif - bisection(lr,lt0,ye,ent0,prec,ltout,2,keyerrt); + bisection(lr,lt0,ye,ent0,prec,ltout,2,keyerrt); #if DEBUG - if(*keyerrt==667) { - fprintf(stderr,"This is worse. Bisection failed!\n"); - abort(); - } + if(*keyerrt==667) { + fprintf(stderr,"This is worse. Bisection failed!\n"); + abort(); + } #endif - } return; File [modified]: nuc_eos_dpdrhoe_dpderho.cc Delta lines: +3 -3 =================================================================== --- trunk/src/nuc_eos_cxx/nuc_eos_dpdrhoe_dpderho.cc 2014-03-13 03:01:52 UTC (rev 100) +++ trunk/src/nuc_eos_cxx/nuc_eos_dpdrhoe_dpderho.cc 2014-03-13 03:01:55 UTC (rev 101) @@ -30,7 +30,7 @@ // Note that this code now requires that the // temperature guess be within the table bounds keyerr[i] = checkbounds_kt0_noTcheck(rho[i], ye[i]); - if(keyerr[i] != 0) { + if(CCTK_BUILTIN_EXPECT(keyerr[i] != 0, false)) { *anyerr = 1; } } @@ -45,7 +45,7 @@ const double lt = log(MIN(MAX(temp[i],eos_tempmin),eos_tempmax)); double ltout; const double epstot = eps[i]+energy_shift; - if(epstot>0.0e0) { + if(CCTK_BUILTIN_EXPECT(epstot>0.0e0, true)) { // this is the standard scenario; eps is larger than zero // and we can operate with logarithmic tables const double lxeps = log(epstot); @@ -61,7 +61,7 @@ keyerr[i] = 667; } // epstot > 0.0 - if(keyerr[i]) { + if(CCTK_BUILTIN_EXPECT(keyerr[i] != 0, false)) { // now try negative temperature treatment double eps0, eps1; int idx[8]; File [modified]: nuc_eos_full.cc Delta lines: +6 -6 =================================================================== --- trunk/src/nuc_eos_cxx/nuc_eos_full.cc 2014-03-13 03:01:52 UTC (rev 100) +++ trunk/src/nuc_eos_cxx/nuc_eos_full.cc 2014-03-13 03:01:55 UTC (rev 101) @@ -39,7 +39,7 @@ for(int i=0;i<n;i++) { // check if we are fine keyerr[i] = checkbounds(rho[i], temp[i], ye[i]); - if(keyerr[i] != 0) { + if(CCTK_BUILTIN_EXPECT(keyerr[i] != 0, false)) { *anyerr = 1; } } @@ -148,7 +148,7 @@ prs[i] = exp(prs[i]); eps[i] = exp(eps[i]) - energy_shift; #if HAVEGR - cs2[i] = cs2[i] / (1.0 + eps[i] + prs[i]/rho[i]); + cs2[i] = rho[i] * cs2[i] / (rho[i] + rho[i] * eps[i] + prs[i]); #endif } @@ -194,7 +194,7 @@ // Note that this code now requires that the // temperature guess be within the table bounds keyerr[i] = checkbounds_kt0_noTcheck(rho[i], ye[i]); - if(keyerr[i] != 0) { + if(CCTK_BUILTIN_EXPECT(keyerr[i] != 0, false)) { *anyerr = 1; } } @@ -209,7 +209,7 @@ const double lt = log(MIN(MAX(temp[i],eos_tempmin),eos_tempmax)); double ltout; const double epstot = eps[i]+energy_shift; - if(epstot>0.0e0) { + if(CCTK_BUILTIN_EXPECT(epstot>0.0e0, true)) { // this is the standard scenario; eps is larger than zero // and we can operate with logarithmic tables const double lxeps = log(epstot); @@ -225,7 +225,7 @@ } // epstot > 0.0 - if(keyerr[i]) { + if(CCTK_BUILTIN_EXPECT(keyerr[i] != 0, false)) { *anyerr=1; } else { temp[i] = exp(ltout); @@ -318,7 +318,7 @@ for(int i=0;i<n;i++) { prs[i] = exp(prs[i]); #if HAVEGR - cs2[i] = cs2[i] / (1.0 + eps[i] + prs[i]/rho[i]); + cs2[i] = rho[i]*cs2[i] / (rho[i] + rho[i]*eps[i] + prs[i]); #endif } File [modified]: nuc_eos_press.cc Delta lines: +4 -4 =================================================================== --- trunk/src/nuc_eos_cxx/nuc_eos_press.cc 2014-03-13 03:01:52 UTC (rev 100) +++ trunk/src/nuc_eos_cxx/nuc_eos_press.cc 2014-03-13 03:01:55 UTC (rev 101) @@ -25,7 +25,7 @@ for(int i=0;i<n;i++) { // check if we are fine keyerr[i] = checkbounds(rho[i], temp[i], ye[i]); - if(keyerr[i] != 0) { + if(CCTK_BUILTIN_EXPECT(keyerr[i] != 0, false)) { *anyerr = 1; } } @@ -87,7 +87,7 @@ // Note that this code now requires that the // temperature guess be within the table bounds keyerr[i] = checkbounds_kt0_noTcheck(rho[i], ye[i]); - if(keyerr[i] != 0) { + if(CCTK_BUILTIN_EXPECT(keyerr[i] != 0, false)) { *anyerr = 1; } } @@ -103,7 +103,7 @@ const double lt = log(MIN(MAX(temp[i],eos_tempmin),eos_tempmax)); double ltout; const double epstot = eps[i]+energy_shift; - if(epstot>0.0e0) { + if(CCTK_BUILTIN_EXPECT(epstot>0.0e0, true)) { // this is the standard scenario; eps is larger than zero // and we can operate with logarithmic tables const double lxeps = log(epstot); @@ -119,7 +119,7 @@ keyerr[i] = 667; } // epstot > 0.0 - if(keyerr[i]) { + if(CCTK_BUILTIN_EXPECT(keyerr[i] != 0, false)) { // now try negative temperature treatment double eps0, eps1; int idx[8]; File [modified]: nuc_eos_press_cs2.cc Delta lines: +6 -6 =================================================================== --- trunk/src/nuc_eos_cxx/nuc_eos_press_cs2.cc 2014-03-13 03:01:52 UTC (rev 100) +++ trunk/src/nuc_eos_cxx/nuc_eos_press_cs2.cc 2014-03-13 03:01:55 UTC (rev 101) @@ -26,7 +26,7 @@ for(int i=0;i<n;i++) { // check if we are fine keyerr[i] = checkbounds(rho[i], temp[i], ye[i]); - if(keyerr[i] != 0) { + if(CCTK_BUILTIN_EXPECT(keyerr[i] != 0, false)) { *anyerr = 1; } } @@ -63,7 +63,7 @@ prs[i] = exp(prs[i]); eps[i] = exp(eps[i]) - energy_shift; #if HAVEGR - cs2[i] = cs2[i] / (1.0 + eps[i] + prs[i]/rho[i]); + cs2[i] = rho[i] * cs2[i] / (rho[i] + rho[i] * eps[i] + prs[i]); #endif } @@ -96,7 +96,7 @@ // Note that this code now requires that the // temperature guess be within the table bounds keyerr[i] = checkbounds_kt0_noTcheck(rho[i], ye[i]); - if(keyerr[i] != 0) { + if(CCTK_BUILTIN_EXPECT(keyerr[i] != 0, false)) { *anyerr = 1; } } @@ -112,7 +112,7 @@ const double lt = log(MIN(MAX(temp[i],eos_tempmin),eos_tempmax)); double ltout; const double epstot = eps[i]+energy_shift; - if(epstot>0.0e0) { + if(CCTK_BUILTIN_EXPECT(epstot>0.0e0, true)) { // this is the standard scenario; eps is larger than zero // and we can operate with logarithmic tables const double lxeps = log(epstot); @@ -124,7 +124,7 @@ keyerr[i] = 667; } // epstot > 0.0 - if(keyerr[i]) { + if(CCTK_BUILTIN_EXPECT(keyerr[i] != 0, false)) { // now try negative temperature treatment double eps0, eps1; int idx[8]; @@ -175,7 +175,7 @@ for(int i=0;i<n;i++) { prs[i] = exp(prs[i]); #if HAVEGR - cs2[i] = cs2[i] / (1.0 + eps[i] + prs[i]/rho[i]); + cs2[i] = rho[i] * cs2[i] / (rho[i] + rho[i]*eps[i] + prs[i]); #endif } File [modified]: nuc_eos_short.cc Delta lines: +9 -9 =================================================================== --- trunk/src/nuc_eos_cxx/nuc_eos_short.cc 2014-03-13 03:01:52 UTC (rev 100) +++ trunk/src/nuc_eos_cxx/nuc_eos_short.cc 2014-03-13 03:01:55 UTC (rev 101) @@ -31,7 +31,7 @@ for(int i=0;i<n;i++) { // check if we are fine keyerr[i] = checkbounds(rho[i], temp[i], ye[i]); - if(keyerr[i] != 0) { + if(CCTK_BUILTIN_EXPECT(keyerr[i] != 0, false)) { *anyerr = 1; } } @@ -94,7 +94,7 @@ prs[i] = exp(prs[i]); eps[i] = exp(eps[i]) - energy_shift; #if HAVEGR - cs2[i] = cs2[i] / (1.0 + eps[i] + prs[i]/rho[i]); + cs2[i] = rho[i] * cs2[i] / (rho[i] + rho[i]*eps[i] + prs[i]); #endif } @@ -131,7 +131,7 @@ // Note that this code now requires that the // temperature guess be within the table bounds keyerr[i] = checkbounds_kt0_noTcheck(rho[i], ye[i]); - if(keyerr[i] != 0) { + if(CCTK_BUILTIN_EXPECT(keyerr[i] != 0, false)) { *anyerr = 1; } } @@ -146,7 +146,7 @@ const double lt = log(MIN(MAX(temp[i],eos_tempmin),eos_tempmax)); double ltout; const double epstot = eps[i]+energy_shift; - if(epstot>0.0e0) { + if(CCTK_BUILTIN_EXPECT(epstot>0.0e0, true)) { // this is the standard scenario; eps is larger than zero // and we can operate with logarithmic tables const double lxeps = log(epstot); @@ -163,7 +163,7 @@ } // epstot > 0.0 - if(keyerr[i]) { + if(CCTK_BUILTIN_EXPECT(keyerr[i] != 0, false)) { *anyerr = 1; } else { temp[i] = exp(ltout); @@ -211,7 +211,7 @@ for(int i=0;i<n;i++) { prs[i] = exp(prs[i]); #if HAVEGR - cs2[i] = cs2[i] / (1.0 + eps[i] + prs[i]/rho[i]); + cs2[i] = rho[i] * cs2[i] / (rho[i] + rho[i] * eps[i] + prs[i]); #endif } @@ -247,7 +247,7 @@ // Note that this code now requires that the // temperature guess be within the table bounds keyerr[i] = checkbounds_kt0_noTcheck(rho[i], ye[i]); - if(keyerr[i] != 0) { + if(CCTK_BUILTIN_EXPECT(keyerr[i] != 0, false)) { *anyerr = 1; } } @@ -264,7 +264,7 @@ nuc_eos_findtemp_entropy(lr,lt,ye[i],ent[i],*prec, (double *restrict)(&ltout),&keyerr[i]); - if(keyerr[i]) { + if(CCTK_BUILTIN_EXPECT(keyerr[i] != 0, false)) { *anyerr = 1; } else { temp[i] = exp(ltout); @@ -313,7 +313,7 @@ prs[i] = exp(prs[i]); eps[i] = exp(eps[i]) - energy_shift; #if HAVEGR - cs2[i] = cs2[i] / (1.0 + eps[i] + prs[i]/rho[i]); + cs2[i] = rho[i] * cs2[i] / (rho[i] + rho[i] * eps[i] + prs[i]); #endif } File [modified]: readtable.cc Delta lines: +8 -4 =================================================================== --- trunk/src/nuc_eos_cxx/readtable.cc 2014-03-13 03:01:52 UTC (rev 100) +++ trunk/src/nuc_eos_cxx/readtable.cc 2014-03-13 03:01:55 UTC (rev 101) @@ -194,11 +194,15 @@ // pressure energy_shift = energy_shift * EPSGF; for(int i=0;i<nrho;i++) { - logrho[i] = log(pow(10.0,logrho[i]) * RHOGF); + // rewrite: + //logrho[i] = log(pow(10.0,logrho[i]) * RHOGF); + // by using log(a^b*c) = b*log(a)+log(c) + logrho[i] = logrho[i] * log(10.) + log(RHOGF); } for(int i=0;i<ntemp;i++) { - logtemp[i] = log(pow(10.0,logtemp[i])); + //logtemp[i] = log(pow(10.0,logtemp[i])); + logtemp[i] = logtemp[i]*log(10.0); } // allocate epstable; a linear-scale eps table @@ -214,12 +218,12 @@ { // pressure int idx = 0 + NTABLES*i; - alltables[idx] = log(pow(10.0,alltables[idx])*PRESSGF); + alltables[idx] = alltables[idx] * log(10.0) + log(PRESSGF); } { // eps int idx = 1 + NTABLES*i; - alltables[idx] = log(pow(10.0,alltables[idx])*EPSGF); + alltables[idx] = alltables[idx] * log(10.0) + log(EPSGF); epstable[i] = exp(alltables[idx]); }

1 0

[svn:einsteintoolkit] EOS_Omni/trunk/src/nuc_eos_cxx/ (Rev. 100)
by rhaas＠tapir.caltech.edu 12 Mar '14

12 Mar '14

User: rhaas Date: 2014/03/12 10:01 PM Modified: /trunk/src/nuc_eos_cxx/ nuc_eos_dpdrhoe_dpderho.cc, nuc_eos_full.cc, nuc_eos_press.cc, nuc_eos_press_cs2.cc, nuc_eos_short.cc Log: EOS_Omni: always set anyerr when an error is detected consistently abort the EOS call when checkbounds fails From: Roland Haas <rhaas(a)tapir.caltech.edu> File Changes: Directory: /trunk/src/nuc_eos_cxx/ ================================== File [modified]: nuc_eos_dpdrhoe_dpderho.cc Delta lines: +2 -2 =================================================================== --- trunk/src/nuc_eos_cxx/nuc_eos_dpdrhoe_dpderho.cc 2014-03-13 03:01:48 UTC (rev 99) +++ trunk/src/nuc_eos_cxx/nuc_eos_dpdrhoe_dpderho.cc 2014-03-13 03:01:52 UTC (rev 100) @@ -57,11 +57,11 @@ nuc_eos_findtemp(lr,lt,ye[i],lxeps,*prec, (double *restrict)(&ltout),&keyerr[i]); } else { + // will be overwritten further down, only marks error keyerr[i] = 667; - *anyerr = 1; } // epstot > 0.0 - if(keyerr[i] != 0) { + if(keyerr[i]) { // now try negative temperature treatment double eps0, eps1; int idx[8]; File [modified]: nuc_eos_full.cc Delta lines: +8 -3 =================================================================== --- trunk/src/nuc_eos_cxx/nuc_eos_full.cc 2014-03-13 03:01:48 UTC (rev 99) +++ trunk/src/nuc_eos_cxx/nuc_eos_full.cc 2014-03-13 03:01:52 UTC (rev 100) @@ -41,10 +41,14 @@ keyerr[i] = checkbounds(rho[i], temp[i], ye[i]); if(keyerr[i] != 0) { *anyerr = 1; - return; } } + // Abort if there is any error in checkbounds. + // This should never happen and the program should abort with + // a fatal error anyway. No point in doing any further EOS calculations. + if(*anyerr) return; + for(int i=0;i<n;i++) { int idx[8]; @@ -218,11 +222,12 @@ (double *restrict)(&ltout),&keyerr[i]); } else { keyerr[i] = 667; - *anyerr=1; } // epstot > 0.0 - if(!keyerr[i]) { + if(keyerr[i]) { + *anyerr=1; + } else { temp[i] = exp(ltout); int idx[8]; double delx,dely,delz; File [modified]: nuc_eos_press.cc Delta lines: +6 -1 =================================================================== --- trunk/src/nuc_eos_cxx/nuc_eos_press.cc 2014-03-13 03:01:48 UTC (rev 99) +++ trunk/src/nuc_eos_cxx/nuc_eos_press.cc 2014-03-13 03:01:52 UTC (rev 100) @@ -27,10 +27,14 @@ keyerr[i] = checkbounds(rho[i], temp[i], ye[i]); if(keyerr[i] != 0) { *anyerr = 1; - return; } } + // Abort if there is any error in checkbounds. + // This should never happen and the program should abort with + // a fatal error anyway. No point in doing any further EOS calculations. + if(*anyerr) return; + for(int i=0;i<n;i++) { int idx[8]; double delx,dely,delz; @@ -111,6 +115,7 @@ nuc_eos_findtemp(lr,lt,ye[i],lxeps,*prec, (double *restrict)(&ltout),&keyerr[i]); } else { + // will be overwritten further down, only marks error keyerr[i] = 667; } // epstot > 0.0 File [modified]: nuc_eos_press_cs2.cc Delta lines: +7 -3 =================================================================== --- trunk/src/nuc_eos_cxx/nuc_eos_press_cs2.cc 2014-03-13 03:01:48 UTC (rev 99) +++ trunk/src/nuc_eos_cxx/nuc_eos_press_cs2.cc 2014-03-13 03:01:52 UTC (rev 100) @@ -28,10 +28,14 @@ keyerr[i] = checkbounds(rho[i], temp[i], ye[i]); if(keyerr[i] != 0) { *anyerr = 1; - return; } } + // Abort if there is any error in checkbounds. + // This should never happen and the program should abort with + // a fatal error anyway. No point in doing any further EOS calculations. + if(*anyerr) return; + for(int i=0;i<n;i++) { int idx[8]; double delx,dely,delz; @@ -116,11 +120,11 @@ nuc_eos_findtemp(lr,lt,ye[i],lxeps,*prec, (double *restrict)(&ltout),&keyerr[i]); } else { + // will be overwritten further down, only marks error keyerr[i] = 667; - *anyerr = 1; } // epstot > 0.0 - if(keyerr[i] != 0) { + if(keyerr[i]) { // now try negative temperature treatment double eps0, eps1; int idx[8]; File [modified]: nuc_eos_short.cc Delta lines: +12 -4 =================================================================== --- trunk/src/nuc_eos_cxx/nuc_eos_short.cc 2014-03-13 03:01:48 UTC (rev 99) +++ trunk/src/nuc_eos_cxx/nuc_eos_short.cc 2014-03-13 03:01:52 UTC (rev 100) @@ -33,10 +33,14 @@ keyerr[i] = checkbounds(rho[i], temp[i], ye[i]); if(keyerr[i] != 0) { *anyerr = 1; - return; } } + // Abort if there is any error in checkbounds. + // This should never happen and the program should abort with + // a fatal error anyway. No point in doing any further EOS calculations. + if(*anyerr) return; + for(int i=0;i<n;i++) { int idx[8]; @@ -154,12 +158,14 @@ nuc_eos_findtemp(lr,lt,ye[i],lxeps,*prec, (double *restrict)(&ltout),&keyerr[i]); } else { + // will be overwritten further down, only marks error keyerr[i] = 667; - *anyerr=1; } // epstot > 0.0 - if(!keyerr[i]) { + if(keyerr[i]) { + *anyerr = 1; + } else { temp[i] = exp(ltout); int idx[8]; double delx,dely,delz; @@ -258,7 +264,9 @@ nuc_eos_findtemp_entropy(lr,lt,ye[i],ent[i],*prec, (double *restrict)(&ltout),&keyerr[i]); - if(!keyerr[i]) { + if(keyerr[i]) { + *anyerr = 1; + } else { temp[i] = exp(ltout); int idx[8]; double delx,dely,delz;

1 0

[svn:einsteintoolkit] EOS_Omni/trunk/src/nuc_eos_cxx/ (Rev. 99)
by rhaas＠tapir.caltech.edu 12 Mar '14

12 Mar '14

User: rhaas Date: 2014/03/12 10:01 PM Modified: /trunk/src/nuc_eos_cxx/ nuc_eos_dpdrhoe_dpderho.cc, nuc_eos_full.cc, nuc_eos_press.cc, nuc_eos_press_cs2.cc, nuc_eos_short.cc Log: EOS_Omni: use explicit temporary for loop upper bound not sure if restrict is propagated to called routines or if the compiler must assume that all pointed to variables change when a function is called. From: Roland Haas <rhaas(a)tapir.caltech.edu> File Changes: Directory: /trunk/src/nuc_eos_cxx/ ================================== File [modified]: nuc_eos_dpdrhoe_dpderho.cc Delta lines: +4 -3 =================================================================== --- trunk/src/nuc_eos_cxx/nuc_eos_dpdrhoe_dpderho.cc 2014-03-13 03:01:46 UTC (rev 98) +++ trunk/src/nuc_eos_cxx/nuc_eos_dpdrhoe_dpderho.cc 2014-03-13 03:01:48 UTC (rev 99) @@ -7,7 +7,7 @@ namespace nuc_eos { extern "C" -void CCTK_FNAME(nuc_eos_m_kt0_dpdrhoe_dpderho)(const int *restrict n, +void CCTK_FNAME(nuc_eos_m_kt0_dpdrhoe_dpderho)(const int *restrict n_in, const double *restrict rho, double *restrict temp, const double *restrict ye, @@ -18,12 +18,13 @@ int *restrict keyerr, int *restrict anyerr) { + const int n = *n_in; using namespace nuc_eos; *anyerr = 0; - for(int i=0;i<*n;i++) { + for(int i=0;i<n;i++) { // check if we are fine // Note that this code now requires that the @@ -39,7 +40,7 @@ // a fatal error anyway. No point in doing any further EOS calculations. if(*anyerr) return; - for(int i=0;i<*n;i++) { + for(int i=0;i<n;i++) { const double lr = log(rho[i]); const double lt = log(MIN(MAX(temp[i],eos_tempmin),eos_tempmax)); double ltout; File [modified]: nuc_eos_full.cc Delta lines: +12 -8 =================================================================== --- trunk/src/nuc_eos_cxx/nuc_eos_full.cc 2014-03-13 03:01:46 UTC (rev 98) +++ trunk/src/nuc_eos_cxx/nuc_eos_full.cc 2014-03-13 03:01:48 UTC (rev 99) @@ -6,7 +6,7 @@ namespace nuc_eos { extern "C" -void CCTK_FNAME(nuc_eos_m_kt1_full)(const int *restrict n, +void CCTK_FNAME(nuc_eos_m_kt1_full)(const int *restrict n_in, const double *restrict rho, const double *restrict temp, const double *restrict ye, @@ -33,8 +33,10 @@ using namespace nuc_eos; + const int n = *n_in; + *anyerr = 0; - for(int i=0;i<*n;i++) { + for(int i=0;i<n;i++) { // check if we are fine keyerr[i] = checkbounds(rho[i], temp[i], ye[i]); if(keyerr[i] != 0) { @@ -43,7 +45,7 @@ } } - for(int i=0;i<*n;i++) { + for(int i=0;i<n;i++) { int idx[8]; double delx,dely,delz; @@ -138,7 +140,7 @@ } } - for(int i=0;i<*n;i++) { + for(int i=0;i<n;i++) { prs[i] = exp(prs[i]); eps[i] = exp(eps[i]) - energy_shift; #if HAVEGR @@ -150,7 +152,7 @@ } extern "C" -void CCTK_FNAME(nuc_eos_m_kt0_full)(const int *restrict n, +void CCTK_FNAME(nuc_eos_m_kt0_full)(const int *restrict n_in, const double *restrict rho, double *restrict temp, const double *restrict ye, @@ -178,9 +180,11 @@ using namespace nuc_eos; + const int n = *n_in; + *anyerr = 0; - for(int i=0;i<*n;i++) { + for(int i=0;i<n;i++) { // check if we are fine // Note that this code now requires that the @@ -196,7 +200,7 @@ // a fatal error anyway. No point in doing any further EOS calculations. if(*anyerr) return; - for(int i=0;i<*n;i++) { + for(int i=0;i<n;i++) { const double lr = log(rho[i]); const double lt = log(MIN(MAX(temp[i],eos_tempmin),eos_tempmax)); double ltout; @@ -306,7 +310,7 @@ } } - for(int i=0;i<*n;i++) { + for(int i=0;i<n;i++) { prs[i] = exp(prs[i]); #if HAVEGR cs2[i] = cs2[i] / (1.0 + eps[i] + prs[i]/rho[i]); File [modified]: nuc_eos_press.cc Delta lines: +12 -8 =================================================================== --- trunk/src/nuc_eos_cxx/nuc_eos_press.cc 2014-03-13 03:01:46 UTC (rev 98) +++ trunk/src/nuc_eos_cxx/nuc_eos_press.cc 2014-03-13 03:01:48 UTC (rev 99) @@ -7,7 +7,7 @@ namespace nuc_eos { extern "C" -void CCTK_FNAME(nuc_eos_m_kt1_press_eps)(const int *restrict n, +void CCTK_FNAME(nuc_eos_m_kt1_press_eps)(const int *restrict n_in, const double *restrict rho, const double *restrict temp, const double *restrict ye, @@ -19,8 +19,10 @@ using namespace nuc_eos; + const int n = *n_in; + *anyerr = 0; - for(int i=0;i<*n;i++) { + for(int i=0;i<n;i++) { // check if we are fine keyerr[i] = checkbounds(rho[i], temp[i], ye[i]); if(keyerr[i] != 0) { @@ -29,7 +31,7 @@ } } - for(int i=0;i<*n;i++) { + for(int i=0;i<n;i++) { int idx[8]; double delx,dely,delz; const double xrho = log(rho[i]); @@ -49,7 +51,7 @@ } // now get rid of ln: - for(int i=0;i<*n;i++) { + for(int i=0;i<n;i++) { prs[i] = exp(prs[i]); eps[i] = exp(eps[i]) - energy_shift; } @@ -59,7 +61,7 @@ } extern "C" -void CCTK_FNAME(nuc_eos_m_kt0_press)(const int *restrict n, +void CCTK_FNAME(nuc_eos_m_kt0_press)(const int *restrict n_in, const double *restrict rho, double *restrict temp, const double *restrict ye, @@ -72,9 +74,11 @@ using namespace nuc_eos; + const int n = *n_in; + *anyerr = 0; - for(int i=0;i<*n;i++) { + for(int i=0;i<n;i++) { // check if we are fine // Note that this code now requires that the // temperature guess be within the table bounds @@ -90,7 +94,7 @@ if(*anyerr) return; // first must find the temperature - for(int i=0;i<*n;i++) { + for(int i=0;i<n;i++) { const double lr = log(rho[i]); const double lt = log(MIN(MAX(temp[i],eos_tempmin),eos_tempmax)); double ltout; @@ -148,7 +152,7 @@ } // loop i<n // now get rid of ln: - for(int i=0;i<*n;i++) { + for(int i=0;i<n;i++) { prs[i] = exp(prs[i]); } File [modified]: nuc_eos_press_cs2.cc Delta lines: +13 -9 =================================================================== --- trunk/src/nuc_eos_cxx/nuc_eos_press_cs2.cc 2014-03-13 03:01:46 UTC (rev 98) +++ trunk/src/nuc_eos_cxx/nuc_eos_press_cs2.cc 2014-03-13 03:01:48 UTC (rev 99) @@ -7,7 +7,7 @@ namespace nuc_eos { extern "C" -void CCTK_FNAME(nuc_eos_m_kt1_press_eps_cs2)(const int *restrict n, +void CCTK_FNAME(nuc_eos_m_kt1_press_eps_cs2)(const int *restrict n_in, const double *restrict rho, const double *restrict temp, const double *restrict ye, @@ -20,8 +20,10 @@ using namespace nuc_eos; + const int n = *n_in; + *anyerr = 0; - for(int i=0;i<*n;i++) { + for(int i=0;i<n;i++) { // check if we are fine keyerr[i] = checkbounds(rho[i], temp[i], ye[i]); if(keyerr[i] != 0) { @@ -30,7 +32,7 @@ } } - for(int i=0;i<*n;i++) { + for(int i=0;i<n;i++) { int idx[8]; double delx,dely,delz; const double xrho = log(rho[i]); @@ -53,7 +55,7 @@ } // now get rid of ln: - for(int i=0;i<*n;i++) { + for(int i=0;i<n;i++) { prs[i] = exp(prs[i]); eps[i] = exp(eps[i]) - energy_shift; #if HAVEGR @@ -66,7 +68,7 @@ } extern "C" -void CCTK_FNAME(nuc_eos_m_kt0_press_cs2)(const int *restrict n, +void CCTK_FNAME(nuc_eos_m_kt0_press_cs2)(const int *restrict n_in, const double *restrict rho, double *restrict temp, const double *restrict ye, @@ -80,14 +82,16 @@ using namespace nuc_eos; + const int n = *n_in; + *anyerr = 0; - for(int i=0;i<*n;i++) { + for(int i=0;i<n;i++) { // check if we are fine // Note that this code now requires that the // temperature guess be within the table bounds - keyerr[i] = checkbounds_kt0_noTcheck(rho[i], ye[i]);; + keyerr[i] = checkbounds_kt0_noTcheck(rho[i], ye[i]); if(keyerr[i] != 0) { *anyerr = 1; } @@ -99,7 +103,7 @@ if(*anyerr) return; // first must find the temperature - for(int i=0;i<*n;i++) { + for(int i=0;i<n;i++) { const double lr = log(rho[i]); const double lt = log(MIN(MAX(temp[i],eos_tempmin),eos_tempmax)); double ltout; @@ -164,7 +168,7 @@ } // now get rid of ln: - for(int i=0;i<*n;i++) { + for(int i=0;i<n;i++) { prs[i] = exp(prs[i]); #if HAVEGR cs2[i] = cs2[i] / (1.0 + eps[i] + prs[i]/rho[i]); File [modified]: nuc_eos_short.cc Delta lines: +18 -12 =================================================================== --- trunk/src/nuc_eos_cxx/nuc_eos_short.cc 2014-03-13 03:01:46 UTC (rev 98) +++ trunk/src/nuc_eos_cxx/nuc_eos_short.cc 2014-03-13 03:01:48 UTC (rev 99) @@ -7,7 +7,7 @@ namespace nuc_eos { extern "C" -void CCTK_FNAME(nuc_eos_m_kt1_short)(const int *restrict n, +void CCTK_FNAME(nuc_eos_m_kt1_short)(const int *restrict n_in, const double *restrict rho, const double *restrict temp, const double *restrict ye, @@ -25,8 +25,10 @@ using namespace nuc_eos; + const int n = *n_in; + *anyerr = 0; - for(int i=0;i<*n;i++) { + for(int i=0;i<n;i++) { // check if we are fine keyerr[i] = checkbounds(rho[i], temp[i], ye[i]); if(keyerr[i] != 0) { @@ -35,7 +37,7 @@ } } - for(int i=0;i<*n;i++) { + for(int i=0;i<n;i++) { int idx[8]; double delx,dely,delz; @@ -84,7 +86,7 @@ } } - for(int i=0;i<*n;i++) { + for(int i=0;i<n;i++) { prs[i] = exp(prs[i]); eps[i] = exp(eps[i]) - energy_shift; #if HAVEGR @@ -96,7 +98,7 @@ } extern "C" - void CCTK_FNAME(nuc_eos_m_kt0_short)(const int *restrict n, + void CCTK_FNAME(nuc_eos_m_kt0_short)(const int *restrict n_in, const double *restrict rho, double *restrict temp, const double *restrict ye, @@ -115,9 +117,11 @@ using namespace nuc_eos; + const int n = *n_in; + *anyerr = 0; - for(int i=0;i<*n;i++) { + for(int i=0;i<n;i++) { // check if we are fine // Note that this code now requires that the @@ -133,7 +137,7 @@ // a fatal error anyway. No point in doing any further EOS calculations. if(*anyerr) return; - for(int i=0;i<*n;i++) { + for(int i=0;i<n;i++) { const double lr = log(rho[i]); const double lt = log(MIN(MAX(temp[i],eos_tempmin),eos_tempmax)); double ltout; @@ -198,7 +202,7 @@ } } - for(int i=0;i<*n;i++) { + for(int i=0;i<n;i++) { prs[i] = exp(prs[i]); #if HAVEGR cs2[i] = cs2[i] / (1.0 + eps[i] + prs[i]/rho[i]); @@ -209,7 +213,7 @@ } extern "C" - void CCTK_FNAME(nuc_eos_m_kt2_short)(const int *restrict n, + void CCTK_FNAME(nuc_eos_m_kt2_short)(const int *restrict n_in, const double *restrict rho, double *restrict temp, const double *restrict ye, @@ -228,9 +232,11 @@ using namespace nuc_eos; + const int n = *n_in; + *anyerr = 0; - for(int i=0;i<*n;i++) { + for(int i=0;i<n;i++) { // check if we are fine // Note that this code now requires that the // temperature guess be within the table bounds @@ -245,7 +251,7 @@ // a fatal error anyway. No point in doing any further EOS calculations. if(*anyerr) return; - for(int i=0;i<*n;i++) { + for(int i=0;i<n;i++) { const double lr = log(rho[i]); const double lt = log(MIN(MAX(temp[i],eos_tempmin),eos_tempmax)); double ltout; @@ -295,7 +301,7 @@ } } - for(int i=0;i<*n;i++) { + for(int i=0;i<n;i++) { prs[i] = exp(prs[i]); eps[i] = exp(eps[i]) - energy_shift; #if HAVEGR

1 0

← Newer
1
...
7
8
9
10
11
12
13
...
18
Older →

Jump to page:

2026

2025

2024

2023

2022

2021

2020

2019

2018

2017

2016

2015

2014

2013

2012

2011

2010

Commits March 2014