Commits March 2014

commits@lists.einsteintoolkit.org

8 participants
172 discussions

[svn:einsteintoolkit] EOS_Omni/trunk/src/nuc_eos_cxx/ (Rev. 98)
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.hh Log: EOS_Omni: add comments and TODOs From: Roland Haas <rhaas(a)tapir.caltech.edu> File Changes: Directory: /trunk/src/nuc_eos_cxx/ ================================== File [modified]: nuc_eos.hh Delta lines: +5 -1 =================================================================== --- trunk/src/nuc_eos_cxx/nuc_eos.hh 2014-03-13 03:01:42 UTC (rev 97) +++ trunk/src/nuc_eos_cxx/nuc_eos.hh 2014-03-13 03:01:46 UTC (rev 98) @@ -3,6 +3,8 @@ #include "cctk.h" +// TODO: remove hard coded constants +// TODO: introduce defines for table index of variables #define HAVEGR 1 #define MAX(x, y) (((x) > (y)) ? (x) : (y)) #define MIN(x, y) (((x) < (y)) ? (x) : (y)) @@ -70,7 +72,9 @@ // 16 Abar // 17 Zbar // 18 Gamma - + enum eos_var {i_logpress=0, i_logenergy, i_entropy, i_munu, i_cs2, i_dedt, + i_dpdrhoe, i_dpderho, i_muhat, i_mu_e, i_mu_p, i_mu_n, i_Xa, + i_Xh, i_Xn, i_Xp, i_Abar, i_Zbar, i_Gamma}; } #endif // NUC_EOS_HH

1 0

[svn:einsteintoolkit] EOS_Omni/trunk/src/nuc_eos_cxx/ (Rev. 97)
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.hh, readtable.cc Log: EOS_Omni: declare global variables restrict should help if we ever directly access them in an inner loop From: Roland Haas <rhaas(a)tapir.caltech.edu> File Changes: Directory: /trunk/src/nuc_eos_cxx/ ================================== File [modified]: nuc_eos.hh Delta lines: +7 -5 =================================================================== --- trunk/src/nuc_eos_cxx/nuc_eos.hh 2014-03-13 03:01:39 UTC (rev 96) +++ trunk/src/nuc_eos_cxx/nuc_eos.hh 2014-03-13 03:01:42 UTC (rev 97) @@ -1,6 +1,8 @@ #ifndef NUC_EOS_HH #define NUC_EOS_HH +#include "cctk.h" + #define HAVEGR 1 #define MAX(x, y) (((x) > (y)) ? (x) : (y)) #define MIN(x, y) (((x) < (y)) ? (x) : (y)) @@ -20,16 +22,16 @@ extern int ntemp; extern int nye; - extern double *alltables; - extern double *epstable; - extern double *logrho; - extern double *logtemp; + 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 *yes; + extern double * restrict yes; extern double energy_shift; extern double dtemp, dtempi; extern double drho, drhoi; File [modified]: readtable.cc Delta lines: +5 -5 =================================================================== --- trunk/src/nuc_eos_cxx/readtable.cc 2014-03-13 03:01:39 UTC (rev 96) +++ trunk/src/nuc_eos_cxx/readtable.cc 2014-03-13 03:01:42 UTC (rev 97) @@ -36,16 +36,16 @@ int ntemp; int nye; - double *alltables; - double *epstable; - double *logrho; - double *logtemp; + double * restrict alltables; + double * restrict epstable; + double * restrict logrho; + double * restrict logtemp; double temp0, temp1; double dlintemp, dlintempi; double drholintempi; double dlintempyei; double drholintempyei; - double *yes; + double * restrict yes; double energy_shift; double dtemp, dtempi; double drho, drhoi;

1 0

[svn:einsteintoolkit] EOS_Omni/trunk/src/nuc_eos_cxx/ (Rev. 96)
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/ readtable.cc Log: EOS_Omni: use Cactus error reporting functions in readtable From: Roland Haas <rhaas(a)tapir.caltech.edu> File Changes: Directory: /trunk/src/nuc_eos_cxx/ ================================== File [modified]: readtable.cc Delta lines: +22 -22 =================================================================== --- trunk/src/nuc_eos_cxx/readtable.cc 2014-03-13 03:01:36 UTC (rev 95) +++ trunk/src/nuc_eos_cxx/readtable.cc 2014-03-13 03:01:39 UTC (rev 96) @@ -9,11 +9,11 @@ #define HDF5_ERROR(fn_call) \ do { \ int _error_code = fn_call; \ - if (_error_code < 0) { \ - fprintf(stderr, \ + if (_error_code < 0) { \ + CCTK_VError(__LINE__, __FILE__, CCTK_THORNSTRING, \ "HDF5 call '%s' returned error code %d", \ #fn_call, _error_code); \ - abort(); } \ + } \ } while (0) static int file_is_readable(const char* filename); @@ -67,6 +67,7 @@ namespace nuc_eos { +// TODO: replace with version in ET EOS_Omni. NOTE: table arrangement changed. // Cactus calls this function. It reads in the table and calls a fortran // function to setup values for the fortran eos module @@ -75,17 +76,16 @@ { using namespace nuc_eos; - fprintf(stdout,"*******************************\n"); - fprintf(stdout,"Reading nuc_eos table file:\n"); - fprintf(stdout,"%s\n",nuceos_table_name); - fprintf(stdout,"*******************************\n"); + CCTK_VInfo(CCTK_THORNSTRING,"*******************************"); + CCTK_VInfo(CCTK_THORNSTRING,"Reading nuc_eos table file:"); + CCTK_VInfo(CCTK_THORNSTRING,"%s",nuceos_table_name); + CCTK_VInfo(CCTK_THORNSTRING,"*******************************"); hid_t file; if (!file_is_readable(nuceos_table_name)) { - fprintf(stderr, + CCTK_VError(__LINE__, __FILE__, CCTK_THORNSTRING, "Could not read nuceos_table_name %s \n", nuceos_table_name); - abort(); } HDF5_ERROR(file = H5Fopen(nuceos_table_name, H5F_ACC_RDONLY, H5P_DEFAULT)); @@ -103,7 +103,7 @@ do { \ hsize_t offset[2] = {OFF,0}; \ H5Sselect_hyperslab(mem3, H5S_SELECT_SET, offset, NULL, var3, NULL); \ - READ_EOS_HDF5(NAME,alltables_temp,H5T_NATIVE_DOUBLE,mem3); \ + READ_EOS_HDF5(NAME,alltables_temp,H5T_NATIVE_DOUBLE,mem3); \ } while (0) // Read size of tables @@ -115,20 +115,20 @@ // Allocate memory for tables double* alltables_temp; if (!(alltables_temp = (double*)malloc(nrho * ntemp * nye * NTABLES * sizeof(double)))) { - fprintf(stderr, "Cannot allocate memory for EOS table\n"); - abort(); + CCTK_VError(__LINE__, __FILE__, CCTK_THORNSTRING, + "Cannot allocate memory for EOS table"); } if (!(logrho = (double*)malloc(nrho * sizeof(double)))) { - fprintf(stderr, "Cannot allocate memory for EOS table\n"); - abort(); + CCTK_VError(__LINE__, __FILE__, CCTK_THORNSTRING, + "Cannot allocate memory for EOS table"); } if (!(logtemp = (double*)malloc(ntemp * sizeof(double)))) { - fprintf(stderr, "Cannot allocate memory for EOS table\n"); - abort(); + CCTK_VError(__LINE__, __FILE__, CCTK_THORNSTRING, + "Cannot allocate memory for EOS table"); } if (!(yes = (double*)malloc(nye * sizeof(double)))) { - fprintf(stderr, "Cannot allocate memory for EOS table\n"); - abort(); + CCTK_VError(__LINE__, __FILE__, CCTK_THORNSTRING, + "Cannot allocate memory for EOS table"); } // Prepare HDF5 to read hyperslabs into alltables_temp @@ -174,8 +174,8 @@ // the table kind is the fastest changing index if (!(alltables = (double*)malloc(nrho * ntemp * nye * NTABLES * sizeof(double)))) { - fprintf(stderr, "Cannot allocate memory for EOS table\n"); - abort(); + CCTK_VError(__LINE__, __FILE__, CCTK_THORNSTRING, + "Cannot allocate memory for EOS table"); } for(int iv = 0;iv<NTABLES;iv++) for(int k = 0; k<nye;k++) @@ -205,8 +205,8 @@ // that allows us to extrapolate to negative eps if (!(epstable = (double*)malloc(nrho * ntemp * nye * sizeof(double)))) { - fprintf(stderr, "Cannot allocate memory for eps table\n"); - abort(); + CCTK_VError(__LINE__, __FILE__, CCTK_THORNSTRING, + "Cannot allocate memory for eps table\n"); } // convert units

1 0

[svn:einsteintoolkit] EOS_Omni/trunk/src/nuc_eos_cxx/ (Rev. 95)
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 Log: EOS_Omni: remove __restrict__ in favor of restrict since Cactus handles restrict for us. From: Roland Haas <rhaas(a)tapir.caltech.edu> File Changes: Directory: /trunk/src/nuc_eos_cxx/ ================================== File [modified]: helpers.hh Delta lines: +19 -18 =================================================================== --- trunk/src/nuc_eos_cxx/helpers.hh 2014-03-13 03:01:30 UTC (rev 94) +++ trunk/src/nuc_eos_cxx/helpers.hh 2014-03-13 03:01:36 UTC (rev 95) @@ -28,6 +28,7 @@ return 101; } if(xye < eos_yemin) { + // this is probably not pure and should be removed fprintf(stderr,"xye: %15.6E eos_yemin: %15.6E\n",xye,eos_yemin); return 102; } @@ -72,10 +73,10 @@ void get_interp_spots(const double x, const double y, const double z, - double* __restrict__ delx, - double* __restrict__ dely, - double* __restrict__ delz, - int* __restrict__ idx) + double* restrict delx, + double* restrict dely, + double* restrict delz, + int* restrict idx) { using namespace nuc_eos; @@ -108,10 +109,10 @@ void get_interp_spots_linT_low(const double x, const double y, const double z, - double* __restrict__ delx, - double* __restrict__ dely, - double* __restrict__ delz, - int* __restrict__ idx) + double* restrict delx, + double* restrict dely, + double* restrict delz, + int* restrict idx) { using namespace nuc_eos; @@ -143,10 +144,10 @@ void get_interp_spots_linT_low_eps(const double x, const double y, const double z, - double* __restrict__ delx, - double* __restrict__ dely, - double* __restrict__ delz, - int* __restrict__ idx) + double* restrict delx, + double* restrict dely, + double* restrict delz, + int* restrict idx) { using namespace nuc_eos; @@ -176,11 +177,11 @@ static inline __attribute__((always_inline)) -void nuc_eos_C_linterp_one(const int* __restrict__ idx, +void nuc_eos_C_linterp_one(const int* restrict idx, const double delx, const double dely, const double delz, - double* __restrict__ f, + double* restrict f, const int iv) { using namespace nuc_eos; @@ -221,11 +222,11 @@ } static inline __attribute__((always_inline)) -void nuc_eos_C_linterp_one_linT_low(const int* __restrict__ idx, +void nuc_eos_C_linterp_one_linT_low(const int* restrict idx, const double delx, const double dely, const double delz, - double* __restrict__ f, + double* restrict f, const int iv) { using namespace nuc_eos; @@ -267,11 +268,11 @@ static inline __attribute__((always_inline)) -void nuc_eos_C_linterp_one_linT_low_eps(const int* __restrict__ idx, +void nuc_eos_C_linterp_one_linT_low_eps(const int* restrict idx, const double delx, const double dely, const double delz, - double* __restrict__ f) + double* restrict f) { using namespace nuc_eos; File [modified]: nuc_eos.hh Delta lines: +0 -1 =================================================================== --- trunk/src/nuc_eos_cxx/nuc_eos.hh 2014-03-13 03:01:30 UTC (rev 94) +++ trunk/src/nuc_eos_cxx/nuc_eos.hh 2014-03-13 03:01:36 UTC (rev 95) @@ -13,7 +13,6 @@ #define INVRHOGF 6.17714470405638e17 #define INVEPSGF 8.98755178736818e20 #define INVPRESSGF 5.55174079257738e38 -#define restrict __restrict__ namespace nuc_eos {

1 0

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

12 Mar '14

User: rhaas Date: 2014/03/12 10:01 PM Removed: /trunk/src/ make.code.deps /trunk/src/nuc_eos_cxx/ make.code.deps Log: EOS_Omni: remove unnecessary make.code.deps Cactus tracks dependencies on include files on its own using calls to the C compiler From: Roland Haas <rhaas(a)tapir.caltech.edu> File Changes: Directory: /trunk/src/nuc_eos_cxx/ ================================== File [removed]: make.code.deps Delta lines: +0 -18 =================================================================== --- trunk/src/nuc_eos_cxx/make.code.deps 2014-03-13 03:01:27 UTC (rev 93) +++ trunk/src/nuc_eos_cxx/make.code.deps 2014-03-13 03:01:30 UTC (rev 94) @@ -1,18 +0,0 @@ -nuc_eos_dpdrhoe_dpderho.cc: nuc_eos.hh -nuc_eos.cc: nuc_eos.hh -nuc_eos_press.cc: nuc_eos.hh -nuc_eos_press_cs2.cc: nuc_eos.hh -nuc_eos_short.cc: nuc_eos.hh -nuc_eos_full.cc: nuc_eos.hh -readtable.cc: nuc_eos.hh -nuc_eos_dpdrhoe_dpderho.cc: helpers.hh -nuc_eos.cc: helpers.hh -nuc_eos_press.cc: helpers.hh -nuc_eos_press_cs2.cc: helpers.hh -nuc_eos_short.cc: helpers.hh -nuc_eos_full.cc: helpers.hh -readtable.cc: helpers.hh - - - - Directory: /trunk/src/ ====================== File [removed]: make.code.deps Delta lines: +0 -1 =================================================================== --- trunk/src/make.code.deps 2014-03-13 03:01:27 UTC (rev 93) +++ trunk/src/make.code.deps 2014-03-13 03:01:30 UTC (rev 94) @@ -1 +0,0 @@ -EOS_Omni_Startup.F90: EOS_Omni_Module.o

1 0

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

12 Mar '14

User: rhaas Date: 2014/03/12 10:01 PM Modified: /trunk/src/ EOS_Omni_BarotropicReadTable.F90, EOS_Omni_ColdEOSReadTable.F90, EOS_Omni_MultiVarCalls.F90, EOS_Omni_SingleVarCalls.F90, EOS_Omni_SingleVarCalls_harm.F90 Log: EOS_Omni: use CCTK_ERROR, add STOP after CCTK_ERROR From: Roland Haas <rhaas(a)tapir.caltech.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: EOS_Omni_BarotropicReadTable.F90 Delta lines: +2 -1 =================================================================== --- trunk/src/EOS_Omni_BarotropicReadTable.F90 2014-03-13 03:01:16 UTC (rev 92) +++ trunk/src/EOS_Omni_BarotropicReadTable.F90 2014-03-13 03:01:27 UTC (rev 93) @@ -30,7 +30,8 @@ if(.not.tablethere) then write(warnline,"(A10,A,A15)") "EOS file ", trim(adjustl(eosfilename)), " does not exist!" - call CCTK_WARN(0,warnline) + call CCTK_ERROR(warnline) + STOP endif ! read header File [modified]: EOS_Omni_ColdEOSReadTable.F90 Delta lines: +8 -4 =================================================================== --- trunk/src/EOS_Omni_ColdEOSReadTable.F90 2014-03-13 03:01:16 UTC (rev 92) +++ trunk/src/EOS_Omni_ColdEOSReadTable.F90 2014-03-13 03:01:27 UTC (rev 93) @@ -31,7 +31,8 @@ if(.not.tablethere) then write(warnline,"(A10,A,A15)") "EOS file ", trim(adjustl(eosfilename)), " does not exist!" - call CCTK_WARN(0,warnline) + call CCTK_ERROR(warnline) + STOP endif ! read header @@ -39,12 +40,14 @@ read(667,"(A)") buffer1 ! check if okay if(.not. (trim(adjustl(buffer1)).eq."EoSType = Tabulated")) then - call CCTK_WARN(0,"Wrong EOS table format -- check header") + call CCTK_ERROR("Wrong EOS table format -- check header") + STOP endif ! now read number of rho entries read(667,"(A7,i7,A6,i6,A5,i6)") buffer1,lnrho,buffer2,lnye,buffer3,lntemp if(lnye.ne.1.or.lntemp.ne.1) then - call CCTK_WARN(0,"Wrong EOS table format -- cannot handle Nye!=1 or NT!=`") + call CCTK_ERROR("Wrong EOS table format -- cannot handle Nye!=1 or NT!=`") + STOP endif coldeos_nrho = lnrho @@ -80,7 +83,8 @@ ! make sure density is in log spacing read(667,"(A30)") buffer1 if(.not.(trim(adjustl(buffer1)).eq."RhoSpacing = Log")) then - call CCTK_WARN(0,"Density spacing not log? Check table format!") + call CCTK_ERROR("Density spacing not log? Check table format!") + STOP endif ! now read everything File [modified]: EOS_Omni_MultiVarCalls.F90 Delta lines: +20 -15 =================================================================== --- trunk/src/EOS_Omni_MultiVarCalls.F90 2014-03-13 03:01:16 UTC (rev 92) +++ trunk/src/EOS_Omni_MultiVarCalls.F90 2014-03-13 03:01:27 UTC (rev 93) @@ -33,19 +33,17 @@ CCTK_REAL, intent(out) :: munu(npoints) ! local vars - integer :: i character(256) :: warnstring - ! temporary vars for nuc_eos - real*8 :: xrho,xye,xtemp,xenr,xent - real*8 :: xprs,xmunu,xcs2 - real*8 :: xdedt,xdpderho,xdpdrhoe if(eoskey.ne.4) then + ! CRITICAL is not really required but might prevent multiple threads from + ! outputting error messages at the same time and cluttering the log file + !$OMP CRITICAL write(warnstring,"(A8,i5)") "eoskey: ", eoskey - !$OMP CRITICAL call CCTK_WARN(1,warnstring) - call CCTK_WARN(0,"EOS_Omni_EOS_short currently does not work for this eoskey") + call CCTK_ERROR("EOS_Omni_EOS_short currently does not work for this eoskey") !$OMP END CRITICAL + STOP endif anyerr = 0 @@ -64,8 +62,9 @@ keyerr,anyerr) else !$OMP CRITICAL - call CCTK_WARN(0,"This keytemp is not supported") + call CCTK_ERROR("This keytemp is not supported") !$OMP END CRITICAL + STOP endif end subroutine EOS_Omni_EOS_short @@ -106,11 +105,12 @@ character(256) :: warnstring if(eoskey.ne.4) then + !$OMP CRITICAL write(warnstring,"(A8,i5)") "eoskey: ", eoskey - !$OMP CRITICAL call CCTK_WARN(1,warnstring) - call CCTK_WARN(0,"EOS_Omni_EOS_full currently does not work for this eoskey") + call CCTK_ERROR("EOS_Omni_EOS_full currently does not work for this eoskey") !$OMP END CRITICAL + STOP endif anyerr = 0 @@ -127,8 +127,9 @@ keyerr,anyerr) else !$OMP CRITICAL - call CCTK_WARN(0,"This keytemp is not supported") + call CCTK_ERROR("This keytemp is not supported") !$OMP END CRITICAL + STOP endif end subroutine EOS_Omni_EOS_full @@ -224,7 +225,8 @@ enddo case (4) if(keytemp.ne.0) then - call CCTK_WARN(0,"Keytemp other than 0 not supported for dpdrhoe, dpderho") + call CCTK_ERROR("Keytemp other than 0 not supported for dpdrhoe, dpderho") + STOP else call nuc_eos_m_kt0_dpdrhoe_dpderho(npoints,& rho,temp,ye,eps,dpdrhoe,dpderho,rf_precision,& @@ -232,7 +234,8 @@ endif case DEFAULT write(warnstring,*) "eoskey ",eoskey," not implemented!" - call CCTK_WARN(0,warnstring) + call CCTK_ERROR(warnstring) + STOP end select end subroutine EOS_Omni_EOS_dpderho_dpdrhoe @@ -330,10 +333,12 @@ enddo case (4) write(warnstring,*) "depsdrho and depsdpress not implemented yet for hot nuclear EOS" - call CCTK_WARN(0,warnstring) + call CCTK_ERROR(warnstring) + STOP case DEFAULT write(warnstring,*) "eoskey ",eoskey," not implemented!" - call CCTK_WARN(0,warnstring) + call CCTK_ERROR(warnstring) + STOP end select end subroutine EOS_Omni_EOS_DEpsByDRho_DEpsByDPress File [modified]: EOS_Omni_SingleVarCalls.F90 Delta lines: +48 -24 =================================================================== --- trunk/src/EOS_Omni_SingleVarCalls.F90 2014-03-13 03:01:16 UTC (rev 92) +++ trunk/src/EOS_Omni_SingleVarCalls.F90 2014-03-13 03:01:27 UTC (rev 93) @@ -102,7 +102,8 @@ call nuc_eos_m_kt0_press(npoints,rho,temp,ye,eps,press,& rf_precision,keyerr,anyerr) else - call CCTK_WARN(0,"This keytemp is not suppported!") + call CCTK_ERROR("This keytemp is not suppported!") + STOP endif case (5) ! cold tabular EOS with gamma law @@ -190,7 +191,8 @@ case DEFAULT write(warnstring,*) "eoskey ",eoskey," not implemented!" - call CCTK_WARN(0,warnstring) + call CCTK_ERROR(warnstring) + STOP end select end subroutine EOS_Omni_EOS_Press @@ -424,7 +426,8 @@ case DEFAULT write(warnstring,*) "eoskey ",eoskey," not implemented!" - call CCTK_WARN(0,warnstring) + call CCTK_ERROR(warnstring) + STOP end select end subroutine EOS_Omni_EOS_PressOMP @@ -530,7 +533,8 @@ eps,xprs,xent,xcs2,xdedt,xdpderho,dpdrhoe,xmunu,rf_precision,& keyerr,anyerr) else - call CCTK_WARN(0,"This keytemp is not supported!") + call CCTK_ERROR("This keytemp is not supported!") + STOP endif case (5) @@ -583,7 +587,8 @@ case DEFAULT write(warnstring,*) "eoskey ",eoskey," not implemented!" - call CCTK_WARN(0,warnstring) + call CCTK_ERROR(warnstring) + STOP end select end subroutine EOS_Omni_EOS_DPressByDRho @@ -658,7 +663,8 @@ eps,xprs,xent,xcs2,xdedt,dpdepsrho,xdpdrhoe,xmunu,rf_precision,& keyerr,anyerr) else - call CCTK_WARN(0,"This keytemp is not supported!") + call CCTK_ERROR("This keytemp is not supported!") + STOP endif case (5) ! with the cold eos we have to assume P = P(rho), so @@ -692,7 +698,8 @@ enddo case DEFAULT write(warnstring,*) "eoskey ",eoskey," not implemented!" - call CCTK_WARN(0,warnstring) + call CCTK_ERROR(warnstring) + STOP end select end subroutine EOS_Omni_EOS_DPressByDEps @@ -798,7 +805,8 @@ call nuc_eos_m_kt0_press_cs2(npoints,rho,temp,ye,& eps,xprs,cs2,rf_precision,keyerr,anyerr) else - call CCTK_WARN(0,"This keytemp is not supported!") + call CCTK_ERROR("This keytemp is not supported!") + STOP endif case (5) ! with the cold eos we have to assume P = P(rho), so @@ -856,7 +864,8 @@ case DEFAULT write(warnstring,*) "eoskey ",eoskey," not implemented!" - call CCTK_WARN(0,warnstring) + call CCTK_ERROR(warnstring) + STOP end select end subroutine EOS_Omni_EOS_cs2 @@ -907,14 +916,17 @@ case (3) ! hybrid EOS write(warnstring,*) "EOS_Omni_EpsFromPress call not supported for hybrid EOS" - call CCTK_WARN(0,warnstring) + call CCTK_ERROR(warnstring) + STOP case (4) ! nuc EOS write(warnstring,*) "EOS_Omni_EpsFromPress call not supported for nuc_eos yet" - call CCTK_WARN(0,warnstring) + call CCTK_ERROR(warnstring) + STOP case DEFAULT write(warnstring,*) "eoskey ",eoskey," not implemented!" - call CCTK_WARN(0,warnstring) + call CCTK_ERROR(warnstring) + STOP end select end subroutine EOS_Omni_EOS_eps_from_press @@ -980,18 +992,22 @@ case (3) ! hybrid EOS write(warnstring,*) "EOS_Omni_RestMassDensityFromPressEpsTemp not supported for hybrid EOS" - call CCTK_WARN(0,warnstring) + call CCTK_ERROR(warnstring) + STOP case (4) ! nuc EOS if(ikeytemp.eq.2) then - call CCTK_WARN(0,"This function does not work yet when coming in with entropy") + call CCTK_ERROR("This function does not work yet when coming in with entropy") + STOP else if(ikeytemp.eq.1) then keytemp = 1 else - call CCTK_WARN(0,"This function does not work yet when coming in with this keytemp") + call CCTK_ERROR("This function does not work yet when coming in with this keytemp") + STOP endif - call CCTK_WARN(0,"This routine does not work with the new nuc_eos_cxx [yet]") + call CCTK_ERROR("This routine does not work with the new nuc_eos_cxx [yet]") + STOP keytemp = 1 do i=1,npoints @@ -1018,8 +1034,9 @@ !$OMP CRITICAL write(warnstring,"(A25,1P10E15.6)") "Issue with table, dpdrho.lt.0",& rho_guess,xtemp,xye - call CCTK_WARN(0,warnstring) + call CCTK_ERROR(warnstring) !$OMP END CRITICAL + STOP endif if (abs(1.0d0-press_guess/xprs).lt.rf_precision) then @@ -1050,7 +1067,8 @@ ! cold tabulated EOS ! deal only with case in which thermal pressure is zero if(keytemp.ne.1) then - call CCTK_WARN(0,"finding rho(press) for tabulated cold EOS only possible if keytemp=1") + call CCTK_ERROR("finding rho(press) for tabulated cold EOS only possible if keytemp=1") + STOP endif do i=1,npoints @@ -1110,8 +1128,9 @@ !$OMP CRITICAL write(warnstring,"(A25,1P10E15.6)") "Issue with table, dpdrho.lt.0",& rho_guess,rho_guess2 - call CCTK_WARN(0,warnstring) + call CCTK_ERROR(warnstring) !$OMP END CRITICAL + STOP endif if (abs(1.0d0-press_guess/xprs).lt.rf_precision) then @@ -1145,7 +1164,8 @@ ! barotropic tabulated EOS ! deal only with case in which thermal pressure is zero if(keytemp.ne.1) then - call CCTK_WARN(0,"finding rho(press) for tabulated barotropic EOS only possible if keytemp=1") + call CCTK_ERROR("finding rho(press) for tabulated barotropic EOS only possible if keytemp=1") + STOP endif do i=1,npoints @@ -1206,8 +1226,9 @@ !$OMP CRITICAL write(warnstring,"(A25,1P10E15.6)") "Issue with table, dpdrho.lt.0",& rho_guess,rho_guess2 - call CCTK_WARN(0,warnstring) + call CCTK_ERROR(warnstring) !$OMP END CRITICAL + STOP endif if (abs(1.0d0-press_guess/press(i)).lt.rf_precision) then @@ -1232,7 +1253,8 @@ if (rho_guess.lt.1.0d-15.or.counter.gt.100000) then ! !$OMP CRITICAL ! write(warnstring,"(A25,1P10E15.6)") "rho(p) issue", rho_guess,press(i),press_guess -! call CCTK_WARN(0,warnstring) +! call CCTK_ERROR(warnstring) +! STOP ! !$OMP END CRITICAL keyerr(i) = 473 anyerr = 1 @@ -1244,7 +1266,8 @@ case DEFAULT write(warnstring,*) "eoskey ",eoskey," not implemented!" - call CCTK_WARN(0,warnstring) + call CCTK_ERROR(warnstring) + STOP end select end subroutine EOS_Omni_EOS_RhoFromPressEpsTempEnt @@ -1332,7 +1355,8 @@ case DEFAULT write(warnstring,*) "eoskey ",eoskey," not implemented for EOS_Omni_EOS_PressEpsTempYe_from_Rho!" - call CCTK_WARN(0,warnstring) + call CCTK_ERROR(warnstring) + STOP end select end subroutine EOS_Omni_EOS_PressEpsTempYe_from_Rho File [modified]: EOS_Omni_SingleVarCalls_harm.F90 Delta lines: +24 -12 =================================================================== --- trunk/src/EOS_Omni_SingleVarCalls_harm.F90 2014-03-13 03:01:16 UTC (rev 92) +++ trunk/src/EOS_Omni_SingleVarCalls_harm.F90 2014-03-13 03:01:27 UTC (rev 93) @@ -40,11 +40,13 @@ select case (eoskey) case (1) ! polytropic EOS - call CCTK_WARN(0,"Polytropic EOS not implemented for press_f_hro_v2_rhoW") + call CCTK_ERROR("Polytropic EOS not implemented for press_f_hro_v2_rhoW") + STOP case (2) ! gamma-law EOS if(keytemp.eq.1) then - call CCTK_WARN(0,"keytemp.eq.1 not implemented for press_f_hro_v2_rhoW") + call CCTK_ERROR("keytemp.eq.1 not implemented for press_f_hro_v2_rhoW") + STOP endif do i=1,npoints gtmp1 = 1.0d0 - v2(i) @@ -55,11 +57,13 @@ enddo case (3) ! hybrid EOS - call CCTK_WARN(0,"Hybrid EOS not implemented for press_f_hro_v2_rhoW") + call CCTK_ERROR("Hybrid EOS not implemented for press_f_hro_v2_rhoW") + STOP case DEFAULT write(warnstring,*) "eoskey ",eoskey," not implemented!" - call CCTK_WARN(0,warnstring) + call CCTK_ERROR(warnstring) + STOP end select end subroutine EOS_Omni_EOS_Press_f_hrho_v2_rhoW @@ -95,22 +99,26 @@ select case (eoskey) case (1) ! polytropic EOS - call CCTK_WARN(0,"Polytropic EOS not implemented for press_f_hro_v2_rhoW") + call CCTK_ERROR("Polytropic EOS not implemented for press_f_hro_v2_rhoW") + STOP case (2) ! gamma-law EOS if(keytemp.eq.1) then - call CCTK_WARN(0,"keytemp.eq.1 not implemented for press_f_hro_v2_rhoW") + call CCTK_ERROR("keytemp.eq.1 not implemented for press_f_hro_v2_rhoW") + STOP endif do i=1,npoints dpdhrho(i) = (gl_gamma - 1.0d0) * (1.0d0 - v2(i)) / (gl_gamma) enddo case (3) ! hybrid EOS - call CCTK_WARN(0,"Hybrid EOS not implemented for press_f_hro_v2_rhoW") + call CCTK_ERROR("Hybrid EOS not implemented for press_f_hro_v2_rhoW") + STOP case DEFAULT write(warnstring,*) "eoskey ",eoskey," not implemented!" - call CCTK_WARN(0,warnstring) + call CCTK_ERROR(warnstring) + STOP end select end subroutine EOS_Omni_EOS_dpdhrho_f_hrho_v2_rhoW @@ -146,11 +154,13 @@ select case (eoskey) case (1) ! polytropic EOS - call CCTK_WARN(0,"Polytropic EOS not implemented for press_f_hro_v2_rhoW") + call CCTK_ERROR("Polytropic EOS not implemented for press_f_hro_v2_rhoW") + STOP case (2) ! gamma-law EOS if(keytemp.eq.1) then - call CCTK_WARN(0,"keytemp.eq.1 not implemented for press_f_hro_v2_rhoW") + call CCTK_ERROR("keytemp.eq.1 not implemented for press_f_hro_v2_rhoW") + STOP endif do i=1,npoints dpdv2(i) = (gl_gamma - 1.0d0) * (0.5d0 * rhoW(i) & @@ -158,11 +168,13 @@ enddo case (3) ! hybrid EOS - call CCTK_WARN(0,"Hybrid EOS not implemented for press_f_hro_v2_rhoW") + call CCTK_ERROR("Hybrid EOS not implemented for press_f_hro_v2_rhoW") + STOP case DEFAULT write(warnstring,*) "eoskey ",eoskey," not implemented!" - call CCTK_WARN(0,warnstring) + call CCTK_ERROR(warnstring) + STOP end select end subroutine EOS_Omni_EOS_dpdv2_f_hrho_v2_rhoW

1 0

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

12 Mar '14

User: rhaas Date: 2014/03/12 10:01 PM Modified: /trunk/src/ EOS_Omni_MultiVarCalls.F90, EOS_Omni_SingleVarCalls.F90 Log: * improved checking of keytemp From: Christian Ott <cott(a)tapir.caltech.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: EOS_Omni_MultiVarCalls.F90 Delta lines: +2 -2 =================================================================== --- trunk/src/EOS_Omni_MultiVarCalls.F90 2014-03-13 03:01:12 UTC (rev 91) +++ trunk/src/EOS_Omni_MultiVarCalls.F90 2014-03-13 03:01:16 UTC (rev 92) @@ -223,8 +223,8 @@ dpderho(i) = (hybrid_gamma_th - 1.0d0) * rho(i) enddo case (4) - if(keytemp.eq.1) then - call CCTK_WARN(0,"keytemp=1 not supported for dpdrhoe, dpderho") + if(keytemp.ne.0) then + call CCTK_WARN(0,"Keytemp other than 0 not supported for dpdrhoe, dpderho") else call nuc_eos_m_kt0_dpdrhoe_dpderho(npoints,& rho,temp,ye,eps,dpdrhoe,dpderho,rf_precision,& File [modified]: EOS_Omni_SingleVarCalls.F90 Delta lines: +12 -4 =================================================================== --- trunk/src/EOS_Omni_SingleVarCalls.F90 2014-03-13 03:01:12 UTC (rev 91) +++ trunk/src/EOS_Omni_SingleVarCalls.F90 2014-03-13 03:01:16 UTC (rev 92) @@ -98,9 +98,11 @@ if(keytemp.eq.1) then call nuc_eos_m_kt1_press_eps(npoints,rho,temp,ye,& eps,press,keyerr,anyerr) - else + else if(keytemp.eq.0) then call nuc_eos_m_kt0_press(npoints,rho,temp,ye,eps,press,& rf_precision,keyerr,anyerr) + else + call CCTK_WARN(0,"This keytemp is not suppported!") endif case (5) ! cold tabular EOS with gamma law @@ -523,10 +525,12 @@ call nuc_eos_m_kt1_short(npoints,rho,temp,ye,& eps,xprs,xent,xcs2,xdedt,xdpderho,dpdrhoe,xmunu,& keyerr,anyerr) - else + else if(keytemp.eq.0) then call nuc_eos_m_kt0_short(npoints,rho,temp,ye,& eps,xprs,xent,xcs2,xdedt,xdpderho,dpdrhoe,xmunu,rf_precision,& keyerr,anyerr) + else + call CCTK_WARN(0,"This keytemp is not supported!") endif case (5) @@ -649,10 +653,12 @@ call nuc_eos_m_kt1_short(npoints,rho,temp,ye,& eps,xprs,xent,xcs2,xdedt,dpdepsrho,xdpdrhoe,xmunu,& keyerr,anyerr) - else + else if(keytemp.eq.0) then call nuc_eos_m_kt0_short(npoints,rho,temp,ye,& eps,xprs,xent,xcs2,xdedt,dpdepsrho,xdpdrhoe,xmunu,rf_precision,& keyerr,anyerr) + else + call CCTK_WARN(0,"This keytemp is not supported!") endif case (5) ! with the cold eos we have to assume P = P(rho), so @@ -788,9 +794,11 @@ if(keytemp.eq.1) then call nuc_eos_m_kt1_press_eps_cs2(npoints,rho,temp,ye,& eps,xprs,cs2,keyerr,anyerr) - else + else if(keytemp.eq.0) then call nuc_eos_m_kt0_press_cs2(npoints,rho,temp,ye,& eps,xprs,cs2,rf_precision,keyerr,anyerr) + else + call CCTK_WARN(0,"This keytemp is not supported!") endif case (5) ! with the cold eos we have to assume P = P(rho), so

1 0

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

12 Mar '14

User: rhaas Date: 2014/03/12 10:01 PM Added: /trunk/src/ make.code.deps /trunk/src/nuc_eos_cxx/ helpers.hh, make.code.defn, make.code.deps, nuc_eos.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, readtable_cactus_wrapper.cc Removed: /trunk/src/nuc_eos/ bisection.F90, dumpASCIItable.F90, eosmodule.F90, findrho.F90, findtemp.F90, linterp.F, linterp_many.F90, make.code.defn, make.code.deps, nuc_eos.F90, readtable.c Modified: /trunk/ interface.ccl, schedule.ccl /trunk/src/ EOS_Omni_Module.F90, EOS_Omni_MultiVarCalls.F90, EOS_Omni_SingleVarCalls.F90, EOS_Omni_Startup.F90, make.code.defn Log: * new nuc_eos backend (in c++) From: Christian Ott <cott(a)tapir.caltech.edu> File Changes: Directory: /trunk/src/nuc_eos/ ============================== File [removed]: bisection.F90 Delta lines: <Unavailable> File [removed]: dumpASCIItable.F90 Delta lines: <Unavailable> File [removed]: eosmodule.F90 Delta lines: <Unavailable> File [removed]: findrho.F90 Delta lines: <Unavailable> File [removed]: findtemp.F90 Delta lines: <Unavailable> File [removed]: linterp.F Delta lines: <Unavailable> File [removed]: linterp_many.F90 Delta lines: <Unavailable> File [removed]: make.code.defn Delta lines: <Unavailable> File [removed]: make.code.deps Delta lines: <Unavailable> File [removed]: nuc_eos.F90 Delta lines: <Unavailable> File [removed]: readtable.c Delta lines: <Unavailable> Directory: /trunk/src/nuc_eos_cxx/ ================================== File [added]: helpers.hh Delta lines: <Unavailable> File [added]: make.code.defn Delta lines: <Unavailable> File [added]: make.code.deps Delta lines: <Unavailable> File [added]: nuc_eos.hh Delta lines: <Unavailable> File [added]: nuc_eos_dpdrhoe_dpderho.cc Delta lines: <Unavailable> File [added]: nuc_eos_full.cc Delta lines: <Unavailable> File [added]: nuc_eos_press.cc Delta lines: <Unavailable> File [added]: nuc_eos_press_cs2.cc Delta lines: <Unavailable> File [added]: nuc_eos_short.cc Delta lines: <Unavailable> File [added]: readtable.cc Delta lines: <Unavailable> File [added]: readtable_cactus_wrapper.cc Delta lines: <Unavailable> Directory: /trunk/src/ ====================== File [modified]: EOS_Omni_Module.F90 Delta lines: <Unavailable> File [modified]: EOS_Omni_MultiVarCalls.F90 Delta lines: <Unavailable> File [modified]: EOS_Omni_SingleVarCalls.F90 Delta lines: <Unavailable> File [modified]: EOS_Omni_Startup.F90 Delta lines: <Unavailable> File [modified]: make.code.defn Delta lines: <Unavailable> File [added]: make.code.deps Delta lines: <Unavailable> Directory: /trunk/ ================== File [modified]: interface.ccl Delta lines: <Unavailable> File [modified]: schedule.ccl Delta lines: <Unavailable>

1 0

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

12 Mar '14

User: rhaas Date: 2014/03/12 10:01 PM Modified: /trunk/ param.ccl /trunk/src/ EOS_Omni_Module.F90, EOS_Omni_MultiVarCalls.F90, EOS_Omni_SingleVarCalls.F90, EOS_Omni_Startup.F90 Log: EOS_Omni: re-add poly_gamma_initial From: Roland Haas <rhaas(a)tapir.caltech.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: EOS_Omni_Module.F90 Delta lines: +9 -0 =================================================================== --- trunk/src/EOS_Omni_Module.F90 2014-03-13 03:01:00 UTC (rev 89) +++ trunk/src/EOS_Omni_Module.F90 2014-03-13 03:01:05 UTC (rev 90) @@ -26,7 +26,13 @@ ! These values are set by EOS_Omni_Startup real*8 :: hybrid_k2 = 0.0d0 + 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 + ! stuff for the cold, tabulated EOS with a gamma law ! set by the reader routine integer :: coldeos_nrho = 0 @@ -63,4 +69,7 @@ real*8, allocatable :: barotropiceos_temp(:) real*8, allocatable :: barotropiceos_ye(:) + ! actual poly_gamma_ini used when using different EOS for initial data and run + real*8 :: poly_gamma_ini + end module EOS_Omni_Module File [modified]: EOS_Omni_MultiVarCalls.F90 Delta lines: +32 -32 =================================================================== --- trunk/src/EOS_Omni_MultiVarCalls.F90 2014-03-13 03:01:00 UTC (rev 89) +++ trunk/src/EOS_Omni_MultiVarCalls.F90 2014-03-13 03:01:05 UTC (rev 90) @@ -104,7 +104,7 @@ integer :: i character(256) :: warnstring real*8 :: hybrid_local_gamma - real*8 :: hybrid_local_k + 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 @@ -119,23 +119,23 @@ ! polytropic EOS if(keytemp.eq.1) then do i=1,npoints - eps(i) = poly_k * & - rho(i)**(poly_gamma) / & + 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) = poly_k * & - poly_gamma * & - rho(i) ** (poly_gamma - 1.d0) + dpdrhoe(i) = press_gf * poly_k_cgs * & + poly_gamma * inv_rho_gf * & + (rho(i)*inv_rho_gf) ** (poly_gamma - 1.d0) dpderho(i) = 0.0d0 enddo case (2) ! gamma-law EOS if(keytemp.eq.1) then do i=1,npoints - eps(i) = gl_k * & - rho(i)**(gl_gamma) / & + eps(i) = press_gf * gl_k_cgs * & + (rho(i)*inv_rho_gf)**(gl_gamma) / & (gl_gamma - 1.0d0) / rho(i) enddo endif @@ -150,23 +150,23 @@ do i=1,npoints if(rho(i).gt.hybrid_rho_nuc) then hybrid_local_gamma = hybrid_gamma2 - hybrid_local_k = hybrid_k2 + hybrid_local_k_cgs = hybrid_k2_cgs else hybrid_local_gamma = hybrid_gamma1 - hybrid_local_k = hybrid_k1 + hybrid_local_k_cgs = hybrid_k1_cgs endif - hybrid_dp_poly = hybrid_local_gamma * & - hybrid_local_k * rho(i)**(hybrid_local_gamma - 1.0d0) - + 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 * hybrid_local_k * & + 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) + 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) * & - hybrid_k1 * & + 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 @@ -223,7 +223,7 @@ integer :: i character(256) :: warnstring real*8 :: hybrid_local_gamma - real*8 :: hybrid_local_k + 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 @@ -239,23 +239,23 @@ ! polytropic EOS if(keytemp.eq.1) then do i=1,npoints - eps(i) = poly_k * & - rho(i)**(poly_gamma) / & + 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 depsdpress(i) = 1.0d0/(poly_gamma - 1.0d0)/rho(i) - depsdrho(i) = depsdpress(i) * poly_k * & - poly_gamma * & - rho(i) ** (poly_gamma - 1.d0) + depsdrho(i) = depsdpress(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) = gl_k * & - rho(i)**(gl_gamma) / & + eps(i) = press_gf * gl_k_cgs * & + (rho(i)*inv_rho_gf)**(gl_gamma) / & (gl_gamma - 1.0d0) / rho(i) enddo endif @@ -269,23 +269,23 @@ do i=1,npoints if(rho(i).gt.hybrid_rho_nuc) then hybrid_local_gamma = hybrid_gamma2 - hybrid_local_k = hybrid_k2 + hybrid_local_k_cgs = hybrid_k2_cgs else hybrid_local_gamma = hybrid_gamma1 - hybrid_local_k = hybrid_k1 + hybrid_local_k_cgs = hybrid_k1_cgs endif - hybrid_dp_poly = hybrid_local_gamma * & - hybrid_local_k * rho(i)**(hybrid_local_gamma - 1.0d0) - + 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 * hybrid_local_k * & + 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) + 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) * & - hybrid_k1 * & + press_gf * hybrid_k1_cgs * inv_rho_gf**hybrid_gamma1 * & hybrid_rho_nuc**(hybrid_gamma1 - 1.d0) xdpdrhoe = hybrid_dp_poly + hybrid_dp_th1 + hybrid_dp_th2 File [modified]: EOS_Omni_SingleVarCalls.F90 Delta lines: +63 -63 =================================================================== --- trunk/src/EOS_Omni_SingleVarCalls.F90 2014-03-13 03:01:00 UTC (rev 89) +++ trunk/src/EOS_Omni_SingleVarCalls.F90 2014-03-13 03:01:05 UTC (rev 90) @@ -27,7 +27,7 @@ ! local vars integer :: i character(256) :: warnstring - real*8 :: hybrid_local_gamma, hybrid_local_k, & + 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 @@ -48,21 +48,21 @@ ! polytropic EOS if(keytemp.eq.1) then do i=1,npoints - eps(i) = poly_k * & - rho(i)**(poly_gamma) / & + 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) = poly_k * & - rho(i)**poly_gamma + 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) = gl_k * & - rho(i)**(gl_gamma) / & + eps(i) = press_gf * gl_k_cgs * & + (rho(i)*inv_rho_gf)**(gl_gamma) / & (gl_gamma - 1.0d0) / rho(i) enddo endif @@ -75,19 +75,19 @@ do i=1,npoints if(rho(i).gt.hybrid_rho_nuc) then hybrid_local_gamma = hybrid_gamma2 - hybrid_local_k = hybrid_k2 + hybrid_local_k_cgs = hybrid_k2_cgs else hybrid_local_gamma = hybrid_gamma1 - hybrid_local_k = hybrid_k1 + hybrid_local_k_cgs = hybrid_k1_cgs 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_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) * & - hybrid_k1 * & + 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 @@ -225,7 +225,7 @@ ! local vars integer :: i character(256) :: warnstring - real*8 :: hybrid_local_gamma, hybrid_local_k, & + 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 @@ -247,16 +247,16 @@ if(keytemp.eq.1) then !$OMP PARALLEL DO do i=1,npoints - eps(i) = poly_k * & - rho(i)**(poly_gamma) / & + eps(i) = press_gf * poly_k_cgs * & + (rho(i)*inv_rho_gf)**(poly_gamma) / & (poly_gamma - 1.0d0) / rho(i) enddo !$OMP END PARALLEL DO endif !$OMP PARALLEL DO do i=1,npoints - press(i) = poly_k * & - rho(i)**poly_gamma + press(i) = press_gf * poly_k_cgs * & + (rho(i)*inv_rho_gf)**poly_gamma enddo !$OMP END PARALLEL DO case (2) @@ -264,8 +264,8 @@ if(keytemp.eq.1) then !$OMP PARALLEL DO do i=1,npoints - eps(i) = gl_k * & - rho(i)**(gl_gamma) / & + eps(i) = press_gf * gl_k_cgs * & + (rho(i)*inv_rho_gf)**(gl_gamma) / & (gl_gamma - 1.0d0) / rho(i) enddo !$OMP END PARALLEL DO @@ -278,24 +278,24 @@ case (3) ! hybrid EOS - !$OMP PARALLEL DO PRIVATE(hybrid_local_gamma,hybrid_local_k,& + !$OMP PARALLEL DO PRIVATE(hybrid_local_gamma,hybrid_local_k_cgs,& !$OMP hybrid_p_poly, hybrid_p_th) do i=1,npoints if(rho(i).gt.hybrid_rho_nuc) then hybrid_local_gamma = hybrid_gamma2 - hybrid_local_k = hybrid_k2 + hybrid_local_k_cgs = hybrid_k2_cgs else hybrid_local_gamma = hybrid_gamma1 - hybrid_local_k = hybrid_k1 + hybrid_local_k_cgs = hybrid_k1_cgs 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_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) * & - hybrid_k1 * & + 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 @@ -448,7 +448,7 @@ ! local vars integer :: i character(256) :: warnstring - real*8 :: hybrid_local_gamma, hybrid_local_k, & + real*8 :: hybrid_local_gamma, hybrid_local_k_cgs, & hybrid_dp_poly, hybrid_dp_th1, hybrid_dp_th2 real*8,parameter :: zero = 0.0d0 ! temporary vars for nuc_eos @@ -469,22 +469,22 @@ ! polytropic EOS if(keytemp.eq.1) then do i=1,npoints - eps(i) = poly_k * & - rho(i)**(poly_gamma) / & + 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) = poly_k * & - poly_gamma * & - rho(i) ** (poly_gamma - 1.d0) + 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) = gl_k * & - rho(i)**(gl_gamma) / & + eps(i) = press_gf * gl_k_cgs * & + (rho(i)*inv_rho_gf)**(gl_gamma) / & (gl_gamma - 1.0d0) / rho(i) enddo endif @@ -497,23 +497,23 @@ do i=1,npoints if(rho(i).gt.hybrid_rho_nuc) then hybrid_local_gamma = hybrid_gamma2 - hybrid_local_k = hybrid_k2 + hybrid_local_k_cgs = hybrid_k2_cgs else hybrid_local_gamma = hybrid_gamma1 - hybrid_local_k = hybrid_k1 + hybrid_local_k_cgs = hybrid_k1_cgs endif - hybrid_dp_poly = hybrid_local_gamma * & - hybrid_local_k * rho(i)**(hybrid_local_gamma - 1.0d0) - + 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 * hybrid_local_k * & + 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) + 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) * & - hybrid_k1 * & + press_gf * hybrid_k1_cgs * inv_rho_gf**hybrid_gamma1 * & hybrid_rho_nuc**(hybrid_gamma1 - 1.d0) dpdrhoe(i) = hybrid_dp_poly + max(0.0d0,hybrid_dp_th1 + hybrid_dp_th2) @@ -632,8 +632,8 @@ ! polytropic EOS if(keytemp.eq.1) then do i=1,npoints - eps(i) = poly_k * & - rho(i)**(poly_gamma) / & + eps(i) = press_gf * poly_k_cgs * & + (rho(i)*inv_rho_gf)**(poly_gamma) / & (poly_gamma - 1.0d0) / rho(i) enddo endif @@ -644,8 +644,8 @@ ! gamma-law EOS if(keytemp.eq.1) then do i=1,npoints - eps(i) = gl_k * & - rho(i)**(gl_gamma) / & + eps(i) = press_gf * gl_k_cgs * & + (rho(i)*inv_rho_gf)**(gl_gamma) / & (gl_gamma - 1.0d0) / rho(i) enddo endif @@ -741,7 +741,7 @@ integer :: i character(256) :: warnstring real*8 :: xpress,xdpdrhoe,xdpderho - real*8 :: hybrid_local_gamma, hybrid_local_k, & + 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 @@ -762,14 +762,14 @@ ! polytropic EOS if(keytemp.eq.1) then do i=1,npoints - eps(i) = poly_k * & - rho(i)**(poly_gamma) / & + 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 = poly_k * & - rho(i)**(poly_gamma) + 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 @@ -777,8 +777,8 @@ ! gamma-law EOS if(keytemp.eq.1) then do i=1,npoints - eps(i) = gl_k * & - rho(i)**(gl_gamma) / & + eps(i) = press_gf * gl_k_cgs * & + (rho(i)*inv_rho_gf)**(gl_gamma) / & (gl_gamma - 1.0d0) / rho(i) enddo endif @@ -794,20 +794,20 @@ do i=1,npoints if(rho(i).gt.hybrid_rho_nuc) then hybrid_local_gamma = hybrid_gamma2 - hybrid_local_k = hybrid_k2 + hybrid_local_k_cgs = hybrid_k2_cgs else hybrid_local_gamma = hybrid_gamma1 - hybrid_local_k = hybrid_k1 + hybrid_local_k_cgs = hybrid_k1_cgs endif ! first calculate the pressure - 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_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) * & - hybrid_k1 * & + 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 File [modified]: EOS_Omni_Startup.F90 Delta lines: +14 -2 =================================================================== --- trunk/src/EOS_Omni_Startup.F90 2014-03-13 03:01:00 UTC (rev 89) +++ trunk/src/EOS_Omni_Startup.F90 2014-03-13 03:01:05 UTC (rev 90) @@ -11,7 +11,19 @@ DECLARE_CCTK_PARAMETERS DECLARE_CCTK_ARGUMENTS - hybrid_k2 = hybrid_k1 * & - hybrid_rho_nuc**(hybrid_gamma1-hybrid_gamma2) + if(poly_gamma_initial .gt. 0d0) then + poly_gamma_ini = poly_gamma_initial + else + poly_gamma_ini = poly_gamma + end if + poly_k_cgs = poly_k * rho_gf**poly_gamma_ini / 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_cgs * & + (hybrid_rho_nuc * inv_rho_gf)**(hybrid_gamma1-hybrid_gamma2) + end subroutine EOS_Omni_Startup Directory: /trunk/ ================== File [modified]: param.ccl Delta lines: +7 -0 =================================================================== --- trunk/param.ccl 2014-03-13 03:01:00 UTC (rev 89) +++ trunk/param.ccl 2014-03-13 03:01:05 UTC (rev 90) @@ -9,6 +9,13 @@ : :: "" } 2.0 + +REAL poly_gamma_initial "Initial Adiabatic Index for poly EOS" STEERABLE=RECOVER +{ + -1 :: "use poly_gamma, ie no change in gamma during ID" + (0:* :: "assume that ID used this adiabiatic index, keep poly_k constant in cgs units" +} -1 + REAL poly_k "Polytropic constant in c=G=Msun=1" STEERABLE=RECOVER { : :: ""

1 0

[svn:einsteintoolkit] EOS_Omni/trunk/src/nuc_eos/ (Rev. 89)
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/ eosmodule.F90, findtemp.F90, linterp.F, linterp_many.F90, nuc_eos.F90 Log: * Major bug fix in EOS_Omni/nuc_eos. In the code (before the fix), offending specific internal energies where limited and the EOS returned seemingly fine, but thermodynamically inconsistent data. This has been fixed. This fix also makes isolated microphysical NS simulations not work anymore. Why this is the case, we are still investigating From: Christian Ott <cott(a)tapir.caltech.edu> File Changes: Directory: /trunk/src/nuc_eos/ ============================== File [modified]: eosmodule.F90 Delta lines: +4 -0 =================================================================== --- trunk/src/nuc_eos/eosmodule.F90 2014-03-13 03:00:46 UTC (rev 88) +++ trunk/src/nuc_eos/eosmodule.F90 2014-03-13 03:01:00 UTC (rev 89) @@ -25,6 +25,7 @@ ! basics integer, parameter :: nvars = 19 real*8,pointer :: alltables(:,:,:,:) + real*8,allocatable :: epstable(:,:,:) ! index variable mapping: ! 1 -> logpress ! 2 -> logenergy @@ -87,6 +88,8 @@ logtemp => logtemp_ ye => ye_ + allocate(epstable(nrho,ntemp,nye)) + epstable(1:nrho,1:ntemp,1:nye) = 10.0d0**alltables(1:nrho,1:ntemp,1:nye,2) energy_shift = energy_shift_ if(do_energy_shift.ne.1) then @@ -103,5 +106,6 @@ eos_tempmin = 10.0d0**logtemp(1) eos_tempmax = 10.0d0**logtemp(ntemp) + end subroutine setup_eosmodule File [modified]: findtemp.F90 Delta lines: +126 -65 =================================================================== --- trunk/src/nuc_eos/findtemp.F90 2014-03-13 03:00:46 UTC (rev 88) +++ trunk/src/nuc_eos/findtemp.F90 2014-03-13 03:01:00 UTC (rev 89) @@ -1,3 +1,107 @@ +subroutine findtemp_deb(lr,lt0,y,epsin,keyerrt,rfeps) + + use eosmodule + + implicit none + + real*8 lr,lt0,y,epsin + real*8 eps,lt,ldt + real*8 tol + real*8 d1,d2,d3 + real*8 eps0,eps1,lt1 + + real*8 ltn,ltmax,ltmin + real*8 tinput,rfeps + + integer :: rl = 0 + integer itmax,i,keyerrt + integer ii,jj,kk + + keyerrt=0 + + tol=rfeps ! precision to which we need to find temp + itmax=20 ! use at most 20 iterations, then bomb + + lt=lt0 + lt1=lt + + eps0=epsin + eps1=eps0 + + ltmax=logtemp(ntemp) + ltmin=logtemp(1) + + ! Note: We are using Ewald's Lagrangian interpolator here! + + !preconditioning 1: do we already have the right temperature? + call findthis(lr,lt,y,eps,alltables(:,:,:,2),d1,d2,d3) + if (abs(eps-eps0).lt.tol*abs(eps0)) then + return + endif + lt1=lt + eps1=eps + + do i=1,itmax + !d2 is the derivative deps/dlogtemp; + ldt = -(eps - eps0)/d2 + ltn = lt+ldt + ltn = min(ltn,ltmax) + ltn = max(ltn,ltmin) + lt1=lt + lt=ltn + eps1=eps + call findthis(lr,lt,y,eps,alltables(:,:,:,2),d1,d2,d3) + if (abs(eps - eps0).lt.tol*abs(eps0)) then + exit + endif + !setup new d2 + + ! if we are closer than 10^-2 to the + ! root (eps-eps0)=0, we are switching to + ! the secant method, since the table is rather coarse and the + ! derivatives may be garbage. + if(abs(eps-eps0).lt.1.0d-3*abs(eps0)) then + d2 = (eps-eps1)/(lt-lt1) + endif + enddo + + + if(i.ge.itmax) then + keyerrt=667 + call bisection(lr,lt0,y,eps0,lt,alltables(:,:,:,2),keyerrt,1) + if(keyerrt.eq.667) then + if(lt.ge.log10(t_max_hack)) then + ! handling for too high temperatures + lt = log10(t_max_hack) + keyerrt=0 + goto 12 + else if(abs(lt-log10(t_max_hack))/log10(t_max_hack).lt.0.025d0) then + lt0 = min(lt,log10(t_max_hack)) + keyerrt=0 + goto 12 + else +#if 0 + ! total failure + write(*,*) "EOS: Did not converge in findtemp!" + write(*,*) "rl,logrho,logtemp0,ye,lt,eps,eps0,abs(eps-eps0)/eps0" + write(*,"(i4,i4,1P10E19.10)") i,rl,lr,lt0,y,lt,eps,eps0,abs(eps-eps0)/eps0 + write(*,*) "Tried calling bisection... didn't help... :-/" + write(*,*) "Bisection error: ",keyerrt +#endif + endif + endif + + lt0=min(lt,log10(t_max_hack)) + return + endif + +12 continue + + lt0=min(lt,log10(t_max_hack)) + + +end subroutine findtemp_deb + subroutine findtemp(lr,lt0,y,epsin,keyerrt,rfeps) use eosmodule @@ -206,77 +310,34 @@ implicit none - real*8 lr,t0,y,epsin - real*8 eps,t,dt - real*8 tol - real*8 dlepsdt - real*8 eps0,eps1,t1 + real*8 lr,t0,y,epsin,rfeps,dlepsdt + integer :: keyerrt - real*8 tn,tmin,tmax - real*8 tinput,rfeps - integer :: rl = 0 - integer itmax,i,keyerrt - integer ii,jj,kk - - keyerrt=0 + real*8 :: eps2, eps1 + real*8 :: t2, t1 + real*8 :: m - tol=rfeps ! need to find energy to less than 1 in 10^-10 - itmax=100 ! use at most 20 iterations, then bomb + keyerrt = 0 + t2 = 10.0d0**logtemp(2) + t1 = 10.0d0**logtemp(1) - t=t0 - t1=t + call intp3d_linearTlow(lr,t2,y,eps2,1,epstable,nrho,& + ntemp,nye,logrho,logtemp,ye,dlepsdt) - eps0=epsin - eps1=eps0 - - tmax=10.0d0**logtemp(2) - tmin=-20.0d0 + call intp3d_linearTlow(lr,t1,y,eps1,1,epstable,nrho,& + ntemp,nye,logrho,logtemp,ye,dlepsdt) - !preconditioning 1: do we already have the right temperature? - call intp3d_linearTlow(lr,t,y,eps,1,alltables(:,:,:,2),nrho,& - ntemp,nye,logrho,logtemp,ye,dlepsdt) + m = (t2-t1)/(eps2-eps1) + t0 = (epsin-eps1) * m + t1 - if (abs(eps-eps0).lt.tol*abs(eps0)) then - return - endif - t1=t - eps1=eps - - do i=1,itmax - dt = -(eps - eps0)/dlepsdt - tn = t+dt - if(tn >= tmax) then - tn = tmax*0.9d0 - endif - tn = min(tn,tmax) - tn = max(tn,tmin) - t1=t - t=tn - eps1=eps +#if 0 + !$OMP CRITICAL + write(6,*) t0,epsin,eps2,eps1 + call intp3d_linearTlow(lr,t0,y,eps2,1,epstable,nrho,& + ntemp,nye,logrho,logtemp,ye,dlepsdt) + write(6,"(1P10E15.6)") abs(epsin-eps2)/epsin + !$OMP END CRITICAL +#endif - call intp3d_linearTlow(lr,t,y,eps,1,alltables(:,:,:,2),nrho,& - ntemp,nye,logrho,logtemp,ye,dlepsdt) - - if (abs(eps - eps0).lt.tol*abs(eps0)) then - exit - endif - - ! if we are closer than 10^-2 to the - ! root (eps-eps0)=0, we are switching to - ! the secant method, since the table is rather coarse and the - ! derivatives may be garbage. - if(abs(eps-eps0).lt.1.0d-3*abs(eps0)) then - dlepsdt = (eps-eps1)/(t-t1) - endif - enddo - - t0 = t - - if(i.ge.itmax) then - keyerrt=667 - return - endif - - end subroutine findtemp_low File [modified]: linterp.F Delta lines: +3 -3 =================================================================== --- trunk/src/nuc_eos/linterp.F 2014-03-13 03:00:46 UTC (rev 88) +++ trunk/src/nuc_eos/linterp.F 2014-03-13 03:01:00 UTC (rev 89) @@ -61,9 +61,9 @@ dy = (yt(ny) - yt(1)) / FLOAT(ny-1) dz = (zt(nz) - zt(1)) / FLOAT(nz-1) c - dxi = 1. / dx - dyi = 1. / dy - dzi = 1. / dz + dxi = 1.0d0 / dx + dyi = 1.0d0 / dy + dzi = 1.0d0 / dz c dxyi = dxi * dyi dxzi = dxi * dzi File [modified]: linterp_many.F90 Delta lines: +9 -9 =================================================================== --- trunk/src/nuc_eos/linterp_many.F90 2014-03-13 03:00:46 UTC (rev 88) +++ trunk/src/nuc_eos/linterp_many.F90 2014-03-13 03:01:00 UTC (rev 89) @@ -54,9 +54,9 @@ dy = (yt(ny) - yt(1)) / FLOAT(ny-1) dz = (zt(nz) - zt(1)) / FLOAT(nz-1) ! - dxi = 1. / dx - dyi = 1. / dy - dzi = 1. / dz + dxi = 1.0d0 / dx + dyi = 1.0d0 / dy + dzi = 1.0d0 / dz ! dxyi = dxi * dyi dxzi = dxi * dzi @@ -181,9 +181,9 @@ dy = 10.0d0**yt(2) - 10.0d0**yt(1) dz = (zt(nz) - zt(1)) / FLOAT(nz-1) ! - dxi = 1. / dx - dyi = 1. / dy - dzi = 1. / dz + dxi = 1.0d0 / dx + dyi = 1.0d0 / dy + dzi = 1.0d0 / dz ! dxyi = dxi * dyi dxzi = dxi * dzi @@ -304,9 +304,9 @@ dy = 10.0d0**yt(2) - 10.0d0**yt(1) dz = (zt(nz) - zt(1)) / FLOAT(nz-1) ! - dxi = 1. / dx - dyi = 1. / dy - dzi = 1. / dz + dxi = 1.0d0 / dx + dyi = 1.0d0 / dy + dzi = 1.0d0 / dz ! dxyi = dxi * dyi dxzi = dxi * dzi File [modified]: nuc_eos.F90 Delta lines: +48 -15 =================================================================== --- trunk/src/nuc_eos/nuc_eos.F90 2014-03-13 03:00:46 UTC (rev 88) +++ trunk/src/nuc_eos/nuc_eos.F90 2014-03-13 03:01:00 UTC (rev 89) @@ -83,13 +83,17 @@ lt = log10(max(xtemp,eos_tempmin)) y = xye xeps = xenr + energy_shift - leps = log10(max(xeps,1.0d0)) keyerr = 0 if(keytemp.eq.0) then !need to find temperature based on xeps - call findtemp(lr,lt,y,leps,keyerrt,rfeps) + if(xeps .gt. 0.0d0) then + leps = log10(xeps) + call findtemp(lr,lt,y,leps,keyerrt,rfeps) + else + keyerrt = 667 + endif if(keyerrt.ne.0) then call CCTK_WARN (0, "Did not find temperature") endif @@ -272,18 +276,23 @@ if(keytemp.eq.0) then !need to find temperature based on xeps xeps = xenr + energy_shift - leps = log10(max(xeps,1.0d0)) - call findtemp(lr,lt,y,leps,keyerrt,rfeps) - if(keyerrt.eq.667 .and. lt .lt. 10.0d0) then + if(xeps .gt. 0.0d0) then + leps = log10(xeps) + call findtemp(lr,lt,y,leps,keyerrt,rfeps) + else + keyerrt = 667 + endif + if(keyerrt.eq.667) then ! turns out, we can still converge, but ! too a bogus, possibly negative temperature xt = xtemp - call findtemp_low(lr,xt,y,leps,keyerrt,rfeps) + call findtemp_low(lr,xt,y,xeps,keyerrt,rfeps) if(keyerrt.eq.0) then keyerrt = 668 xtemp = xt else keyerrt = 667 + keyerr = keyerrt return endif keyerr = keyerrt @@ -365,6 +374,7 @@ real*8 :: lr,lt,y,xt,xx,xeps,leps,xs real*8 :: d1,d2,d3,ff(8) integer :: keyerrt + character(len=256) :: warnline keyerrt = 0 @@ -411,20 +421,38 @@ if(keytemp.eq.0) then !need to find temperature based on xeps xeps = xenr + energy_shift - leps = log10(max(xeps,1.0d0)) - call findtemp(lr,lt,y,leps,keyerrt,rfeps) - if(keyerrt.eq.667 .and. lt .lt. 10.0d0) then + if(xeps .gt. 0.0d0) then + leps = log10(xeps) + call findtemp(lr,lt,y,leps,keyerrt,rfeps) + else + keyerrt = 667 + endif +! if(xenr .lt. -9.00725347549139d20) then +! if(xenr.lt.-2.28659899d20) then +! write(warnline,"(A4,i5,1P10E15.6)") "UGGx",keyerrt,lr,lt,xt,xtemp,y,xeps,xenr +! call CCTK_WARN(1,warnline) +! endif + if(keyerrt.eq.667) then ! turns out, we can still converge, but ! too a bogus, possibly negative temperature xt = xtemp - call findtemp_low(lr,xt,y,leps,keyerrt,rfeps) + call findtemp_low(lr,xt,y,xeps,keyerrt,rfeps) if(keyerrt.eq.0) then keyerrt = 668 + if(xeps.lt.0.0d0) then + keyerrt = 668 ! impossible to do anything + endif xtemp = xt else keyerrt = 667 + keyerr = keyerrt return endif +! if(xenr .lt. -9.00725347549139d20) then +! if(xenr.lt.-2.28659899d20) then +! write(warnline,"(A4,i5,1P10E15.6)") "UGG",keyerrt,lr,lt,xt,xtemp,y,xeps,xenr +! call CCTK_WARN(0,warnline) +! endif keyerr = keyerrt else xtemp = 10.0d0**lt @@ -502,20 +530,25 @@ if(keytemp.eq.0) then !need to find temperature based on xeps xeps = xenr + energy_shift - leps = log10(max(xeps,1.0d0)) - call findtemp(lr,lt,y,leps,keyerrt,rfeps) - if(keyerrt.eq.667 .and. lt .lt. 10.0d0) then + if(xeps .gt. 0.0d0) then + leps = log10(xeps) + call findtemp(lr,lt,y,leps,keyerrt,rfeps) + else + keyerrt = 667 + endif + if(keyerrt.eq.667) then ! turns out, we can still converge, but ! too a bogus, possibly negative temperature xt = xtemp - call findtemp_low(lr,xt,y,leps,keyerrt,rfeps) + call findtemp_low(lr,xt,y,xeps,keyerrt,rfeps) if(keyerrt.eq.0) then keyerrt = 668 xtemp = xt else keyerrt = 667 + keyerr = keyerrt return - endif + endif keyerr = keyerrt else xtemp = 10.0d0**lt

1 0

← Newer
1
...
8
9
10
11
12
13
14
...
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