Skip to content

Commit

Permalink
Merge branch 'eosx' into eosx
Browse files Browse the repository at this point in the history
  • Loading branch information
@acwen11
acwen11 authored Jul 25, 2024
2 parents 1bf3157 + b5720ef commit 08319ed
Show file tree
Hide file tree
Showing 2 changed files with 183 additions and 70 deletions.
2 changes: 1 addition & 1 deletion AsterX/par/Balsara1_shocktube_xdir.par
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,7 @@ ActiveThorns = "
"

$nlevels = 1
$ncells = 400
$ncells = 100

Cactus::cctk_show_schedule = yes

Expand Down
251 changes: 182 additions & 69 deletions EOSX/src/eos_3p_tabulated3d.hxx
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,16 @@
#ifndef EOS_3P_TABULATED3D_HXX
#define EOS_3P_TABULATED3D_HXX

#define NTABLES 19
#define LENGTHGF 6.77269222552442e-06
#define TIMEGF 2.03040204956746e05
#define RHOGF 1.61887093132742e-18
#define PRESSGF 1.80123683248503e-39
#define EPSGF 1.11265005605362e-21
#define INVRHOGF 6.17714470405638e17
#define INVEPSGF 8.98755178736818e20
#define INVPRESSGF 5.55174079257738e38

#include <cctk.h>
#include <cmath>
#include <hdf5.h>
Expand All @@ -13,16 +23,20 @@ using namespace std;

namespace EOSX {

using namespace amrex;

class eos_3p_tabulated3d : public eos_3p {

public:
CCTK_REAL gamma; // FIXME: get rid of this
range rgeps;

CCTK_INT ntemp, nrho, nye;
// AMREX_GPU_MANAGED CCTK_REAL *logrho, *logtemp, *ye;

// amrex::FArrayBox logtemp, logrho, ye;
// amrex::Array4<CCTK_REAL> logpress, ...;

CCTK_REAL *logrho, *logtemp, *yes; // FIXME: AMREX_GPU_MANAGED?
CCTK_REAL *alltables;
CCTK_REAL *epstable;
CCTK_REAL energy_shift;

CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline void
init(const range &rgeps_, const range &rgrho_, const range &rgye_) {
Expand All @@ -34,30 +48,16 @@ public:
set_range_temp(rgeps_);
}

// Routine reading an HDF5 simple dataset consisting of one integer element
// FIXME: make this more general: any type (template), any number of elements
// (but still 1D-shaped)
CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline CCTK_INT
get_hdf5_simple_int(const hid_t &file_id, const string &dset_name) {
// Routine reading an HDF5 integer dataset
CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline
void get_hdf5_int_dset(const hid_t &file_id,
const string &dset_name,
const int npoints,
int* var) {

const auto dset_id = H5Dopen(file_id, dset_name.c_str(), H5P_DEFAULT);
assert(dset_id >= 0);

const auto dspace_id = H5Dget_space(dset_id);
assert(dspace_id >= 0);

const auto dspacetype_id = H5Sget_simple_extent_type(dspace_id);
assert(dspacetype_id == H5S_SIMPLE);

const auto ndims = H5Sget_simple_extent_ndims(dspace_id);
assert(ndims == 1);

hsize_t size;
const auto ndims_again =
H5Sget_simple_extent_dims(dspace_id, &size, nullptr);
CHECK_ERROR(H5Sclose(dspace_id));
assert(ndims_again == 1);
assert(size == 1);

const auto dtype_id = H5Dget_type(dset_id);
assert(dtype_id >= 0);

Expand All @@ -74,9 +74,7 @@ public:
dxpl_id = H5P_DEFAULT;
#endif

CCTK_INT buffer;
CHECK_ERROR(
H5Dread(dset_id, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, dxpl_id, &buffer));
CHECK_ERROR(H5Dread(dset_id, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, dxpl_id, var));
CHECK_ERROR(H5Dclose(dset_id));

#ifdef H5_HAVE_PARALLEL
Expand Down Expand Up @@ -112,39 +110,18 @@ public:
#endif // H5_HAVE_PARALLEL

CHECK_ERROR(H5Pclose(dxpl_id));
return buffer;
}

// Routine reading a 1D HDF5 simple dataset storing real numbers
CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline void
get_hdf5_simple_1Darray(const hid_t &file_id, const string &dset_name,
amrex::FArrayBox &var) {
// Routine reading an HDF5 real number dataset
CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline
void get_hdf5_real_dset(const hid_t &file_id,
const string &dset_name,
const int npoints,
double* var) {

const auto dset_id = H5Dopen(file_id, dset_name.c_str(), H5P_DEFAULT);
assert(dset_id >= 0);

const auto dspace_id = H5Dget_space(dset_id);
assert(dspace_id >= 0);

const auto dspacetype_id = H5Sget_simple_extent_type(dspace_id);
assert(dspacetype_id == H5S_SIMPLE);

const auto ndims = H5Sget_simple_extent_ndims(dspace_id);
assert(ndims == 1);

hsize_t size;
const auto ndims_again =
H5Sget_simple_extent_dims(dspace_id, &size, nullptr);
CHECK_ERROR(H5Sclose(dspace_id));
assert(ndims_again == 1);
assert(size > 0);

const auto dtype_id = H5Dget_type(dset_id);
assert(dtype_id >= 0);

const auto dtypeclass = H5Tget_class(dtype_id);
assert(dtypeclass == H5T_FLOAT);
CHECK_ERROR(H5Tclose(dtype_id));

auto dxpl_id = H5Pcreate(H5P_DATASET_XFER);
assert(dxpl_id >= 0);

Expand All @@ -154,10 +131,7 @@ public:
dxpl_id = H5P_DEFAULT;
#endif

CCTK_REAL buffer[size];
// CCTK_REAL *buffer = new CCTK_REAL[size];
CHECK_ERROR(
H5Dread(dset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, dxpl_id, buffer));
CHECK_ERROR(H5Dread(dset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, dxpl_id, var));
CHECK_ERROR(H5Dclose(dset_id));

#ifdef H5_HAVE_PARALLEL
Expand Down Expand Up @@ -194,13 +168,9 @@ public:

CHECK_ERROR(H5Pclose(dxpl_id));

// amrex::Box box1d(amrex::IntVect{0, 0, 0}, amrex::IntVect{size - 1, 0,
// 0}); var = amrex::FArrayBox(box1d, 1, amrex::The_Managed_Arena());
// TODO: fill var with buffer

return;
}


// Routine reading the EOS table and filling the corresponding object
CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline void
read_eos_table(const string &filename) {
Expand All @@ -226,15 +196,77 @@ public:

assert(file_id >= 0);

// TODO: finish reading the EOS table and filling the EOS object
ntemp = get_hdf5_simple_int(file_id, "pointstemp");
nrho = get_hdf5_simple_int(file_id, "pointsrho");
nye = get_hdf5_simple_int(file_id, "pointsye");
// Get number of points
get_hdf5_int_dset(file_id, "pointstemp", 1, &ntemp);
get_hdf5_int_dset(file_id, "pointsrho", 1, &nrho);
get_hdf5_int_dset(file_id, "pointsye", 1, &nye);

const int npoints = ntemp * nrho * nye;

CCTK_VINFO("EOS table dimensions: ntemp = %d, nrho = %d, nye = %d", ntemp,
nrho, nye);

// Allocate memory for tables
double* alltables_temp;
if (!(alltables_temp = (double*)The_Managed_Arena()->alloc(npoints * NTABLES * sizeof(double)))) {
CCTK_VError(__LINE__, __FILE__, CCTK_THORNSTRING,
"Cannot allocate memory for EOS table");
}
if (!(logrho = (double*)The_Managed_Arena()->alloc(nrho * sizeof(double)))) {
CCTK_VError(__LINE__, __FILE__, CCTK_THORNSTRING,
"Cannot allocate memory for EOS table");
}
if (!(logtemp = (double*)The_Managed_Arena()->alloc(ntemp * sizeof(double)))) {
CCTK_VError(__LINE__, __FILE__, CCTK_THORNSTRING,
"Cannot allocate memory for EOS table");
}
if (!(yes = (double*)The_Managed_Arena()->alloc(nye * sizeof(double)))) {
CCTK_VError(__LINE__, __FILE__, CCTK_THORNSTRING,
"Cannot allocate memory for EOS table");
}

// Prepare HDF5 to read hyperslabs into alltables_temp
hsize_t table_dims[2] = {NTABLES, (hsize_t)npoints};
hsize_t var3[2] = { 1, (hsize_t)npoints};
hid_t mem3 = H5Screate_simple(2, table_dims, NULL);

// hydro (and munu)
get_hdf5_real_dset(file_id, "logpress", npoints, &alltables_temp[0 * npoints]);
get_hdf5_real_dset(file_id, "logenergy", npoints, &alltables_temp[1 * npoints]);
get_hdf5_real_dset(file_id, "entropy", npoints, &alltables_temp[2 * npoints]);
get_hdf5_real_dset(file_id, "munu", npoints, &alltables_temp[3 * npoints]);
get_hdf5_real_dset(file_id, "cs2", npoints, &alltables_temp[4 * npoints]);
get_hdf5_real_dset(file_id, "dedt", npoints, &alltables_temp[5 * npoints]);
get_hdf5_real_dset(file_id, "dpdrhoe", npoints, &alltables_temp[6 * npoints]);
get_hdf5_real_dset(file_id, "dpderho", npoints, &alltables_temp[7 * npoints]);

// chemical potentials
get_hdf5_real_dset(file_id, "muhat", npoints, &alltables_temp[8 * npoints]);
get_hdf5_real_dset(file_id, "mu_e", npoints, &alltables_temp[9 * npoints]);
get_hdf5_real_dset(file_id, "mu_p", npoints, &alltables_temp[10 * npoints]);
get_hdf5_real_dset(file_id, "mu_n", npoints, &alltables_temp[11 * npoints]);

// compositions
get_hdf5_real_dset(file_id, "Xa", npoints, &alltables_temp[12 * npoints]);
get_hdf5_real_dset(file_id, "Xh", npoints, &alltables_temp[13 * npoints]);
get_hdf5_real_dset(file_id, "Xn", npoints, &alltables_temp[14 * npoints]);
get_hdf5_real_dset(file_id, "Xp", npoints, &alltables_temp[15 * npoints]);

// average nucleus
get_hdf5_real_dset(file_id, "Abar", npoints, &alltables_temp[16 * npoints]);
get_hdf5_real_dset(file_id, "Zbar", npoints, &alltables_temp[17 * npoints]);

// Gamma
get_hdf5_real_dset(file_id, "gamma", npoints, &alltables_temp[18 * npoints]);

// Read additional tables and variables
get_hdf5_real_dset(file_id, "logrho", nrho, logrho);
get_hdf5_real_dset(file_id, "logtemp", ntemp, logtemp);
get_hdf5_real_dset(file_id, "ye", nye, yes);
get_hdf5_real_dset(file_id, "energy_shift", 1, &energy_shift);

CHECK_ERROR(H5Pclose(fapl_id));
CHECK_ERROR(H5Sclose(mem3));

#ifdef H5_HAVE_PARALLEL
CHECK_ERROR(H5Fclose(file_id));
Expand All @@ -244,6 +276,87 @@ public:
}
#endif

// Fill actual table
if (!(alltables = (double*)The_Managed_Arena()->alloc(npoints * NTABLES * sizeof(double)))) {
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++)
for(int j = 0; j<ntemp; j++)
for(int i = 0; i<nrho; i++) {
int indold = i + nrho*(j + ntemp*(k + nye*iv));
int indnew = iv + NTABLES*(i + nrho*(j + ntemp*k));

// Maybe swap temp axis?
// int indnew = iv + NTABLES*(j + ntemp*(i + nrho*k));
alltables[indnew] = alltables_temp[indold];
}

// free memory of temporary array
free(alltables_temp);

// allocate epstable; a linear-scale eps table
// that allows us to extrapolate to negative eps
if (!(epstable = (double*)The_Managed_Arena()->alloc(npoints * sizeof(double)))) {
CCTK_VError(__LINE__, __FILE__, CCTK_THORNSTRING,
"Cannot allocate memory for EOS table");
}

// convert units, convert logs to natural log
// The latter is great, because exp() is way faster than pow()
// pressure
energy_shift = energy_shift * EPSGF;
for(int i=0;i<nrho;i++) {
// 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] = logtemp[i]*log(10.0);
}

// convert units
for(int i=0;i<npoints;i++) {

{ // pressure
int idx = 0 + NTABLES*i;
alltables[idx] = alltables[idx] * log(10.0) + log(PRESSGF);
}

{ // eps
int idx = 1 + NTABLES*i;
alltables[idx] = alltables[idx] * log(10.0) + log(EPSGF);
epstable[i] = exp(alltables[idx]);
}

{ // cs2
int idx = 4 + NTABLES*i;
alltables[idx] *= LENGTHGF*LENGTHGF/TIMEGF/TIMEGF;
}

{ // dedT
int idx = 5 + NTABLES*i;
alltables[idx] *= EPSGF;
}

{ // dpdrhoe
int idx = 6 + NTABLES*i;
alltables[idx] *= PRESSGF/RHOGF;
}

{ // dpderho
int idx = 7 + NTABLES*i;
alltables[idx] *= PRESSGF/EPSGF;
}

}

// set up steps, mins, maxes here?

return;
}

Expand Down

0 comments on commit 08319ed

Please sign in to comment.