Skip to content
Draft
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
33 changes: 33 additions & 0 deletions eospac-wrapper/eospac_wrapper.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -335,6 +335,14 @@ void eosSafeTableCmnts(EOS_INTEGER *table, EOS_CHAR *comments, Verbosity eospacW
eosCheckError(errorCode, "eos_GetTableCmnts", eospacWarn);
}

void eosSafeTableMetaData(EOS_INTEGER *table, EOS_INTEGER infoItem, EOS_CHAR *infoString,
Verbosity eospacWarn) {
EOS_INTEGER errorCode = EOS_OK;
EOS_INTEGER infoTypes[1] = {infoItem};
eos_GetTableMetaData(table, infoTypes, infoString, &errorCode);
eosCheckError(errorCode, "eos_GetTableMetaData", eospacWarn);
}

void eosCheckError(EOS_INTEGER errorCode, const std::string &name, Verbosity eospacWarn) {
EOS_CHAR errorMessage[EOS_MaxErrMsgLen];
if (errorCode != EOS_OK && eospacWarn != Verbosity::Quiet) {
Expand All @@ -343,6 +351,31 @@ void eosCheckError(EOS_INTEGER errorCode, const std::string &name, Verbosity eos
}
}

int eosCheckTableExistence(EOS_INTEGER tableType_, int matid_, Verbosity eospacWarn) {
// Note this opens and closes the table
int exists = 1;
int NT = 1;
EOS_INTEGER matid[1] = {matid_};
EOS_INTEGER tableType[1] = {tableType_};
EOS_INTEGER tableHandle[1];
EOS_INTEGER errorCode = EOS_OK;

eos_CreateTables(&NT, tableType, matid, tableHandle, &errorCode);
eosCheckError(errorCode, "eos_CreateTables", eospacWarn);
eos_LoadTables(&NT, &tableHandle[0], &errorCode);
if (errorCode != EOS_OK) {
if (eosErrorCodesEqual(EOS_NO_DATA_TABLE, errorCode)) {
exists = 0;
} else {
// something else happended
eosCheckError(errorCode, "eos_LoadTables", eospacWarn);
exists = -1;
}
}
eosSafeDestroy(NT, tableHandle, eospacWarn);
return exists;
}

std::string eosErrorString(EOS_INTEGER errorCode) {
// Algorithmicallly generated by parsing the EOSPAC docs
// I'm sorry. It's gross. ~JMM
Expand Down
5 changes: 3 additions & 2 deletions eospac-wrapper/eospac_wrapper.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -122,13 +122,14 @@ void eosSafeTableInfo(EOS_INTEGER *table, EOS_INTEGER numInfoItems,
EOS_INTEGER infoItems[], EOS_REAL infoVals[], Verbosity eospacWarn);

void eosSafeTableCmnts(EOS_INTEGER *table, EOS_CHAR *comments, Verbosity eospacWarn);

void eosSafeTableMetaData(EOS_INTEGER *table, EOS_INTEGER infoItem, EOS_CHAR *infoString,
Verbosity eospacWarn);
void eosCheckError(EOS_INTEGER errorCode, const std::string &name, Verbosity eospacWarn);
std::string eosErrorString(EOS_INTEGER errorCode);
void eosSafeDestroy(int ntables, EOS_INTEGER tableHandles[], Verbosity eospacWarn);
std::string getName(std::string comment);
bool eosErrorCodesEqual(EOS_INTEGER lhs, EOS_INTEGER rhs);

int eosCheckTableExistence(EOS_INTEGER tableType_, int matid_, Verbosity eospacWarn);
} // namespace EospacWrapper

#endif // _EOSPAC_WRAPPER_EOSPAC_WRAPPER_HPP_
21 changes: 20 additions & 1 deletion sesame2spiner/generate_files.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -53,7 +53,7 @@ herr_t saveMaterial(hid_t loc, const SesameMetadata &metadata, const Bounds &lRh
std::string sMatid = std::to_string(matid);

herr_t status = 0;
hid_t matGroup, lTGroup, leGroup, coldGroup;
hid_t matGroup, lTGroup, leGroup, coldGroup, mfGroup;

matGroup = H5Gcreate(loc, sMatid.c_str(), H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
status += H5Lcreate_soft(sMatid.c_str(), loc, name.c_str(), H5P_DEFAULT, H5P_DEFAULT);
Expand Down Expand Up @@ -116,6 +116,25 @@ herr_t saveMaterial(hid_t loc, const SesameMetadata &metadata, const Bounds &lRh
// status += transitionMask.saveHDF(coldGroup, SP5::Fields::transitionMask);
}

{
DataBox mf, mask;
Bounds nphBounds;
std::string phase_names;

if (eosMassFraction(matid, lRhoBounds, lTBounds, nphBounds, mf, mask, phase_names,
eospacWarn)) {
mfGroup = H5Gcreate(matGroup, SP5::Depends::massFrac, H5P_DEFAULT, H5P_DEFAULT,
H5P_DEFAULT);

status += mf.saveHDF(mfGroup, SP5::Fields::massFrac);
const int numphases = mf.dim(1);
status += H5LTset_attribute_int(mfGroup, ".", "numphases", &numphases, 1);
status +=
H5LTset_attribute_string(mfGroup, ".", "phase names", phase_names.c_str());
status += H5Gclose(mfGroup);
}
}

if (addSubtables) {
int i = 0;
std::vector<TableSplit> splits = {TableSplit::ElectronOnly, TableSplit::IonCold};
Expand Down
74 changes: 74 additions & 0 deletions sesame2spiner/io_eospac.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -328,6 +328,80 @@ void eosColdCurveMask(int matid, const Bounds &lRhoBounds, const int numSie,
eosSafeDestroy(NT, tableHandle, eospacWarn);
}

bool eosMassFraction(int matid, const Bounds &lRhoBounds, const Bounds &lTBounds,
Bounds &nphBounds, DataBox &Ms, DataBox &mask,
std::string &phase_names, Verbosity eospacWarn) {
using namespace EospacWrapper;

EOS_INTEGER errorCode = EOS_OK;
constexpr int NT = 2;
int ntables = NT;
EOS_INTEGER tableHandle[NT];
EOS_INTEGER tableType[NT] = {EOS_M_DT, EOS_Comment};
std::vector<EOS_INTEGER> matid_v(NT, matid);
std::vector<std::string> table_names = {"EOS_M_DT", "Eos_Material_Phases"};

const int exists = eosCheckTableExistence(EOS_M_DT, matid, eospacWarn);
if (exists > 0) {
eosSafeLoad(NT, matid, tableType, tableHandle, {"EOS_M_DT", "Eos_Material_Phases"},
eospacWarn);
} else {
phase_names = std::string("");
return false;
}
// indep vars
// Reuses the EOS log temp and log rho bounds
std::vector<EOS_REAL> rhos, Ts, phs;
makeInterpPoints(rhos, lRhoBounds);
makeInterpPoints(Ts, lTBounds);
Comment on lines +354 to +356
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am re-using the EOS rho-T grid for this but in testing we are seeing that maybe that isn't great, e.g., instead of clustering around STP, we probably want to cluster the points around phase transitions. I could add another set of parameters to the inputs that describe the grid decomp for just mass fractions?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think that makes a lot of sense. Yeah there's no reason the two grids have to be the same.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If you're going to use different grids, I would just use the sesame grid. (piecewise linear). eospac just does bilinear interpolation for mass fractions so we should not need higher density than sesame.


EOS_REAL infoVals[1];
EOS_INTEGER infoItems[1] = {EOS_NUM_PHASES};
eosSafeTableInfo(&tableHandle[0], 1, infoItems, infoVals, eospacWarn);
const int nph = static_cast<int>(infoVals[0]);

EOS_CHAR infoString[EOS_META_DATA_STRLEN];
infoString[0] = '\0';
eosSafeTableMetaData(&tableHandle[1], EOS_Material_Phases, infoString, eospacWarn);
phase_names = std::string(infoString);

DataBox dMdt, dMdr;
Ms.resize(rhos.size(), Ts.size(), nph);
Ms.setRange(1, lTBounds.grid);
Ms.setRange(2, lRhoBounds.grid);

dMdr.copyMetadata(Ms);
dMdt.copyMetadata(Ms);
mask.copyMetadata(Ms);

// Interpolatable vars
EOS_INTEGER nXYPairs = rhos.size() * Ts.size();
std::vector<EOS_REAL> M_pack(nXYPairs * nph), DMDR_T(nXYPairs), DMDT_R(nXYPairs),
rho_flat(nXYPairs), T_flat(nXYPairs);

// prepare flat data structures
for (std::size_t j = 0, iflat = 0; j < rhos.size(); ++j) {
for (std::size_t i = 0; i < Ts.size(); ++i, iflat++) {
rho_flat[iflat] = densityToSesame(rhos[j]);
T_flat[iflat] = temperatureToSesame(Ts[i]);
}
}

const bool no_errors = eosSafeInterpolate(&tableHandle[0], nXYPairs, rho_flat.data(),
T_flat.data(), M_pack.data(), DMDR_T.data(),
DMDT_R.data(), "EOS_M_DT", eospacWarn);

for (size_t j = 0, iflat = 0; j < rhos.size(); j++) {
for (size_t i = 0; i < Ts.size(); i++, iflat++) {
for (size_t k = 0; k < nph; k++) {
Ms(j, i, k) = M_pack[k * nXYPairs + iflat];
}
}
}
eosSafeDestroy(NT, tableHandle, eospacWarn);
return true;
}

void makeInterpPoints(std::vector<EOS_REAL> &v, const Bounds &b) {
v.resize(b.grid.nPoints());
for (size_t i = 0; i < v.size(); i++) {
Expand Down
3 changes: 3 additions & 0 deletions sesame2spiner/io_eospac.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -76,6 +76,9 @@ void eosColdCurves(int matid, const Bounds &lRhoBounds, DataBox &Ps, DataBox &si
void eosColdCurveMask(int matid, const Bounds &lRhoBounds, const int numSie,
const DataBox &sieColdCurve, DataBox &mask,
Verbosity eospacWarn = Verbosity::Quiet);
bool eosMassFraction(int matid, const Bounds &lRhoBounds, const Bounds &lTBounds,
Bounds &nphBounds, DataBox &mf, DataBox &mask,
std::string &phase_names, Verbosity eospacWarn = Verbosity::Quiet);

void makeInterpPoints(std::vector<EOS_REAL> &v, const Bounds &b);

Expand Down
2 changes: 2 additions & 0 deletions singularity-eos/base/sp5/singularity_eos_sp5.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,7 @@ namespace Depends {
constexpr char logRhoLogSie[] = "dependsLogRhoLogSie";
constexpr char logRhoLogT[] = "dependsLogRhoLogT";
constexpr char coldCurve[] = "coldCurve";
constexpr char massFrac[] = "massFrac";
} // namespace Depends

namespace SubTable {
Expand Down Expand Up @@ -68,6 +69,7 @@ constexpr char dEdRho[] = "dEdRho";
constexpr char dEdT[] = "dEdT";
constexpr char mask[] = "mask";
constexpr char transitionMask[] = "transition mask";
constexpr char massFrac[] = "mass fractions";
} // namespace Fields

} // namespace SP5
Expand Down
Loading