diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 7e9e488..f0b4a7b 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -62,6 +62,8 @@ add_library(tasmet_src material/solid.cpp material/adiabatictemp.cpp + export/export.cpp + # sys/enginesystem.cpp sys/globalconf.cpp sys/tasmet_variable.cpp @@ -73,4 +75,4 @@ add_library(tasmet_src add_subdirectory(gui) add_subdirectory(protobuf) -target_link_libraries(tasmet_src Qt5::Widgets superlu dl pthread) +target_link_libraries(tasmet_src Qt5::Widgets hdf5 hdf5_hl dl superlu pthread) diff --git a/src/duct/duct.cpp b/src/duct/duct.cpp index e02ceee..453e481 100644 --- a/src/duct/duct.cpp +++ b/src/duct/duct.cpp @@ -10,6 +10,7 @@ #include "tasmet_assert.h" #include "tasmet_evalscript.h" #include "perfectgas.h" +#include "../export/export.h" Duct::Duct(const TaSystem& sys, const us id, @@ -47,7 +48,7 @@ Duct::Duct(const TaSystem& sys,const Duct& other): _Tsprescribed(other._Tsprescribed) { // Do something with the equations here - TRACE(15,"Duct::~Duct"); + TRACE(15,"Duct::Duct"); } @@ -109,6 +110,8 @@ void Duct::residualJac(SegPositionMapper& residual, pp = getvart(constants::p,gp+1); cont = ((rhop%up)-(rho%u))/dx; + // cont += sys.DDTtd()*rho; + *residual.at(i) = cont; jacrows.at(i)[{_id,gp*nvars_per_gp+constants::u}] = @@ -119,7 +122,7 @@ void Duct::residualJac(SegPositionMapper& residual, diagmat(up)/dx; jacrows.at(i)[{_id,(gp+1)*nvars_per_gp+constants::u}] = diagmat(rhop)/dx; - TRACE(15,"SFSG"); + i++; mom = (rhop%up%up - rho%u%u + pp - p)/dx; @@ -136,7 +139,7 @@ void Duct::residualJac(SegPositionMapper& residual, jacrows.at(i)[{_id,(gp+1)*nvars_per_gp+constants::p}] = eye(Ns,Ns)/dx; *residual.at(i) = mom; - TRACE(15,"SFSG"); + i++; switch (_ductpb.htmodel()) { @@ -155,7 +158,7 @@ void Duct::residualJac(SegPositionMapper& residual, eye(Ns,Ns)/p0; jacrows.at(i)[{_id,gp*nvars_per_gp+constants::rho}] = -(gamma0/rho0)*diagmat(pow(rho/rho0,gamma0-1)); - TRACE(15,"SFSG"); + } break; default: @@ -163,10 +166,10 @@ void Duct::residualJac(SegPositionMapper& residual, } i++; - TRACE(15,"SFSG"); + // Equation of state st = gas.rho(T,p) - rho; - TRACE(15,"SFSG"); + jacrows.at(i)[{_id,gp*nvars_per_gp+constants::rho}] = -eye(Ns,Ns); jacrows.at(i)[{_id,gp*nvars_per_gp+constants::p}] = @@ -175,7 +178,7 @@ void Duct::residualJac(SegPositionMapper& residual, diagmat(gas.drhodT(T,p)); *residual.at(i) = st; - TRACE(15,"SFSG"); + i++; } @@ -330,5 +333,78 @@ void Duct::show(us verbosity_level) const { } +void Duct::exportHDF5(const hid_t group_id) const { + TRACE(15,"Duct::exportHDF5"); + + PosData x; + x.name = "Position"; + x.unit = "m"; + x.symbol = "x"; + x.x = this->x; + x.exportHDF5(group_id); + + PosData S; + S.name = "Cross-sectional area"; + S.unit = "m^2"; + S.symbol = "S"; + S.x = this->S; + S.exportHDF5(group_id); + + PosData phi; + phi.name = "Porosity"; + phi.unit = "-"; + phi.symbol = "phi"; + phi.x = this->phi; + phi.exportHDF5(group_id); + + PosData rh; + rh.name = "Hydraulic radius"; + rh.unit = "m"; + rh.symbol = "rh"; + rh.x = this->rh; + rh.exportHDF5(group_id); + + TXData u; + u.name = "Velocity"; + u.unit = "m/s"; + u.symbol = "u"; + u.x = dmat(sys.Ns(),ngp()); + for(us gp=0;gp + +#if TASMET_FLOAT == 64 +#define make_dataset H5LTmake_dataset_double +#elif TASMET_FLOAT == 32 +#define make_dataset H5LTmake_dataset_float +#else +#error "Invalid TASMET_FLOAT" +#endif + +namespace { + /** + * Set the attributes of a dataset + * + * @param data the data from which to take the name, unit, symbol + * @param dsetname the dataset name to add the attributes to + * @param hfgroup the group + */ + void setAttributes(const ExportData& data, + const string& dsetname, + const hid_t hfgroup) { + + H5LTset_attribute_string(hfgroup, + dsetname.c_str(), + "name", + data.name.c_str()); + + H5LTset_attribute_string(hfgroup, + dsetname.c_str(), + "unit", + data.unit.c_str()); + + H5LTset_attribute_string(hfgroup, + dsetname.c_str(), + "symbol", + data.symbol.c_str()); + + + } +} + +void PointData::exportHDF5(hid_t hfgroup) { + + // We use the lite API for setting the data + + hsize_t dims[] = {1}; // First dimension has lenght 1 + hsize_t rank = 1; + herr_t status; + + string dsetname = (string("/") + symbol); // name + + status = make_dataset(hfgroup, // Group + dsetname.c_str(), + rank, + dims, + &x); + + + setAttributes(*this,dsetname,hfgroup); + + +} +void TimeData::exportHDF5(hid_t hfgroup) { + + // We use the lite API for setting the data + + hsize_t dims[] = {x.size()}; // First dimension has lenght 1 + hsize_t rank = 1; + + string dsetname = (string("/") + symbol); // name + herr_t status; + status = make_dataset(hfgroup, // Group + dsetname.c_str(), // Dataset name + rank, + dims, + x.memptr()); + + + setAttributes(*this,dsetname,hfgroup); + +} +void PosData::exportHDF5(hid_t hfgroup) { + + TRACE(15,"PosData::exportHDF5"); + + // We use the lite API for setting the data + + hsize_t dims[] = {x.size()}; // First dimension has lenght 1 + hsize_t rank = 1; + + + string dsetname = (string("/") + symbol); // name + herr_t status; + status = make_dataset(hfgroup, // Group + dsetname.c_str(), + rank, + dims, + x.memptr()); + + + setAttributes(*this,dsetname,hfgroup); + + H5Eprint(status,0); + +} +void TXData::exportHDF5(hid_t hfgroup) { + + // We use the lite API for setting the data + + // HDF5 Stores data in row-major ordering. Armadillo uses + // column-major ordering. Therefore, we store the transpose of the data. + dmat x = this->x.t(); + + hsize_t dims[] = {x.n_cols,x.n_rows}; // First dimension has lenght 1 + hsize_t rank = 2; + + d* xptr = x.memptr(); + + std::vector col_ptrs(x.n_cols); + + for(us col=0;col +#endif + +DECLARE_ENUM(FileType,FileTypeTXT) + +struct ExportData +{ + string name; /**< Full name of the quantity */ + string unit; /**< Unit of the quantity */ + string symbol; /**< Used symbol / abbreviation */ + + /** + * Export the data to a text file + * + * @param FileType + * @param filename + * @param group + * + * + */ + // void export(const string& filename)=0; + + /** + * Export the result to an existing HDF5 + * + * @param hfgroup + */ + #ifdef TASMET_HDF5 + virtual void exportHDF5(hid_t hfgroup)=0; + #endif + + + // void show() const = 0; + + virtual ~ExportData(){} +}; + +/** + * A physical quantity result, not a function of time or position + * + */ +struct PointData: public ExportData { + d x; /**< The value */ + + #ifdef TASMET_HDF5 + void exportHDF5(hid_t hfgroup); + #endif +}; + +/** + * A physical quantity at a single position, as a function of time. + * + */ +struct TimeData: public ExportData +{ + vd t; /**< The time instances of the + record */ + vd x; /**< The quantity at the time + instance */ + + #ifdef TASMET_HDF5 + void exportHDF5(hid_t hfgroup); + #endif +}; + +/** + * A physical quantity at a single time instance, as a function of + * position. + * + */ +struct PosData: public ExportData +{ + vd X; /**< The positions of the + record */ + vd x; /**< The quantity as function of the + position */ + + #ifdef TASMET_HDF5 + void exportHDF5(hid_t hfgroup); + #endif +}; + + +/** + * A physical quantity at a single position, as a function of time and + * position (X) + * + */ +struct TXData : public ExportData +{ + + vd t; /**< The time instances of the + record */ + vd X; /**< The positions of the record */ + + dmat x; /**< The quantity at the time and position instance. Note: + the first axis denotes time, the second axis denotes + position! */ + + #ifdef TASMET_HDF5 + void exportHDF5(hid_t hfgroup); + #endif +}; + +#endif // EXPORT_H +////////////////////////////////////////////////////////////////////// diff --git a/src/solver/system.h b/src/solver/system.h index fe7ecbe..b43d590 100644 --- a/src/solver/system.h +++ b/src/solver/system.h @@ -83,7 +83,7 @@ public: virtual vd getSolution() const=0; virtual void updateSolution(const vd& new_guess)=0; - + virtual ~GradientNonlinearSystem(){} }; diff --git a/src/sys/globalconf.cpp b/src/sys/globalconf.cpp index cc9d5f0..45547d3 100644 --- a/src/sys/globalconf.cpp +++ b/src/sys/globalconf.cpp @@ -74,7 +74,7 @@ void GlobalConf::setomg(d omg){ _DDTfd(2*i-1,2*i )=-double(i)*_omg; _DDTfd(2*i ,2*i-1)=double(i)*_omg; } - + _DDTtd = _iDFT * _DDTfd * _fDFT; } ////////////////////////////////////////////////////////////////////// diff --git a/src/sys/globalconf.h b/src/sys/globalconf.h index aed3a7f..b25a12d 100644 --- a/src/sys/globalconf.h +++ b/src/sys/globalconf.h @@ -22,7 +22,7 @@ class GlobalConf{ d _omg; // The "base" frequency in rad/s us _Nf; // Number of frequencies to solve for - dmat _fDFT,_iDFT,_DDTfd; + dmat _fDFT,_iDFT,_DDTfd,_DDTtd; // d Wfo_=0; // First order 'upwind' factor. If // Wfo=-1, interpolation is done from // the left side. If Wfo=0, @@ -35,6 +35,7 @@ public: const dmat& iDFT = _iDFT; //inverse discrete Fourier transform matrix const dmat& fDFT = _fDFT; //forward discrete Fourier transform matrix const dmat& DDTfd = _DDTfd; //Derivative in frequency domain + const dmat& DDTtd = _DDTtd; //Derivative in time domain us Nf() const {return _Nf;} us Ns() const {return 2*_Nf+1;} diff --git a/src/sys/segment.h b/src/sys/segment.h index 8ffd1df..85921f7 100644 --- a/src/sys/segment.h +++ b/src/sys/segment.h @@ -12,13 +12,12 @@ #pragma once #ifndef SEGMENT_H #define SEGMENT_H - +#include #include "tasmet_types.h" #include "tasmet_exception.h" #include "jacobian.h" // #include "phaseconstraint.h" - class GlobalConf; class TaSystem; @@ -83,6 +82,8 @@ public: // Reset amplitude data in higher harmonics // virtual void resetHarmonics() = 0; + virtual void exportHDF5(const hid_t group_id) const {} + }; diff --git a/src/sys/tasystem.cpp b/src/sys/tasystem.cpp index 6685b0a..7aa5aa1 100644 --- a/src/sys/tasystem.cpp +++ b/src/sys/tasystem.cpp @@ -5,6 +5,8 @@ // Description: // ////////////////////////////////////////////////////////////////////// +#include + #include "segment.h" #include "tasystem.h" #include "triplets.h" @@ -16,6 +18,7 @@ #include "duct.h" #include "ductbc.h" + TaSystem::TaSystem(const pb::System& sys): GlobalConf(sys.nf(),sys.freq()) { @@ -413,4 +416,42 @@ void TaSystem::cleanup() { } +void TaSystem::exportHDF5(const string& filename) const { + + TRACE(15,"TaSystem::exportHDF5"); + hid_t file_id; + herr_t status; + + file_id = H5Fcreate (filename.c_str(), + H5F_ACC_TRUNC, + H5P_DEFAULT, + H5P_DEFAULT); + + + + for(const auto& seg_ : _segs) { + // Create a group for each segment + hid_t grp_id = H5Gcreate(file_id, + (string("/") + std::to_string(seg_.first)).c_str(), + H5P_DEFAULT, + H5P_DEFAULT, + H5P_DEFAULT); + + seg_.second->exportHDF5(grp_id); + + H5Gclose(grp_id); + + // ^^ uncommenting above results in segfault. Why? + // seg_.second->exportHDF5(file_id); + + } + + status = H5Fclose(file_id); + H5Eprint(status,0); + + +} + + + ////////////////////////////////////////////////////////////////////// diff --git a/src/sys/tasystem.h b/src/sys/tasystem.h index d091905..2ba4474 100644 --- a/src/sys/tasystem.h +++ b/src/sys/tasystem.h @@ -112,6 +112,17 @@ public: const Segment& getSegment(const us id) const {return *_segs.at(id);} + /** + * Store the current results to a HDF5 file. For each segment, a + * group is created with the name corresponding to the segment + * ID. Then, the group contents is filled based on what the + * segment wants to store as result data. *Warning*: overwrites + * existing files without confirmation! + * + * @param filename The filename of the output file + */ + void exportHDF5(const string& filename) const; + protected: virtual int getArbitrateMassEq(const vus& neqs) const; private: diff --git a/src/tasmet_solvemodel.cpp b/src/tasmet_solvemodel.cpp index 2d75ddf..b6d373f 100644 --- a/src/tasmet_solvemodel.cpp +++ b/src/tasmet_solvemodel.cpp @@ -58,6 +58,15 @@ int main(int argc, char *argv[]) { progress_callback cb = simple_progress_callback(1e-6,1e-6); solver.start(&cb); + + + vd sol = solver.getSolution(); + + // Put solution in our copy of system + system->updateSolution(sol); + + system->exportHDF5(filename+".h5"); + }