diff --git a/src/include/tfrfme.h b/src/include/tfrfme.h index 19174277cdf07c1d118cfca9ac6dd2f34d9e5ae0..daa6e359fe6f79946ce1bd999ed990987b9e78a3 100644 --- a/src/include/tfrfme.h +++ b/src/include/tfrfme.h @@ -6,16 +6,108 @@ #ifndef INCLUDE_TFRFME_H_ #define INCLUDE_TFRFME_H_ +/*! \brief Class to represent the first group of trapping swap data. + */ +class Swap1 { +protected: + //! NLMMT = 2 * LM * (LM + 2) + int nlmmt; + + //! QUESTION: definition? + std::complex<double> *wk; + + /*! \brief Load a Swap1 instance from a HDF5 binary file. + * + * \param file_name: `string` Name of the file to be loaded. + * \return instance: `Swap1 *` Pointer to a new Swap1 instance. + */ + static Swap1 *from_hdf5(std::string file_name); + + /*! \brief Load a Swap1 instance from a legacy binary file. + * + * \param file_name: `string` Name of the file to be loaded. + * \return instance: `Swap1 *` Pointer to a new Swap1 instance. + */ + static Swap1 *from_legacy(std::string file_name); + + /*! \brief Save a Swap1 instance to a HDF5 binary file. + * + * \param file_name: `string` Name of the file to be written. + */ + void write_hdf5(std::string file_name); + + /*! \brief Save a Swap1 instance to a legacy binary file. + * + * \param file_name: `string` Name of the file to be written. + */ + void write_legacy(std::string file_name); + +public: + /*! \brief Swap1 instance constructor. + * + * \param lm: `int` Maximum field expansion order. + */ + Swap1(int lm); + + /*! \brief Swap1 instance destroyer. + */ + ~Swap1() { delete[] wk; } + + /*! \brief Load a Swap1 instance from binary file. + * + * \param file_name: `string` Name of the file. + * \param mode: `string` Format of the file (can be either "HDF5" + * or "LGEACY". Default is "LEGACY"). + * \return instance: `Swap1 *` Pointer to a newly created Swap1 instance. + */ + static Swap1* from_binary(std::string file_name, std::string mode="LEGACY"); + + /*! \brief Get an element from the WK vector. + * + * \param index: `int` Index of the desired element. + * \return value: `complex<double>` The value of the requested element. + */ + std::complex<double> get_element(int index) { return wk[index]; } + + /*! \brief Set an element in the WK vector. + * + * \param index: `int` Index of the desired element. + * \param value: `complex<double>` The value to be stored in the vector. + */ + void set_element(int index, std::complex<double> value) { wk[index] = value; } + + /*! \brief Write a Swap1 instance to binary file. + * + * \param file_name: `string` Name of the file. + * \param mode: `string` Format of the file (can be either "HDF5" + * or "LGEACY". Default is "LEGACY"). + */ + void write_binary(std::string file_name, std::string mode="LEGACY"); + + /*! \brief Test whether two instances of Swap1 are equal. + * + * \param other: `Swap1 &` Reference to the instance to be compared + * with. + * \return result: `bool` True, if the two instances are equal, + * false otherwise. + */ + bool operator ==(Swap1 &other); +}; + +/*! \brief Class to represent the second group of trapping swap data. + */ +class Swap2 { +}; + /*! \brief Class to represent the trapping configuration. */ class TFRFME { - private: +protected: //! NLMMT = 2 * LM * (LM + 2) int nlmmt; //! NRVC = NXV * NYV * NZV int nrvc; - protected: //! Field expansion mode identifier. int lmode; //! Maximim field expansion order; @@ -71,18 +163,25 @@ class TFRFME { /*! \brief Save a configuration instance to a HDF5 binary file. * - * \param file_name: `string` Name of the file to be loaded. + * \param file_name: `string` Name of the file to be written. */ void write_hdf5(std::string file_name); /*! \brief Save a configuration instance to a legacy binary file. * - * \param file_name: `string` Name of the file to be loaded. + * \param file_name: `string` Name of the file to be written. */ void write_legacy(std::string file_name); - public: +public: /*! \brief Trapping configuration instance constructor. + * + * \param _lmode: `int` Order expansion mode flag. + * \param _lm: `int` Maximum field expansion order. + * \param _nkv: `int` Number of wave vector coordinates. QUESTION: correct? + * \param _nxv: `int` Number of computed X coordinates. + * \param _nyv: `int` Number of computed Y coordinates. + * \param _nzv: `int` Number of computed Z coordinates. */ TFRFME( int _lmode, int _lm, int _nkv, int _nxv, int _nyv, int _nzv @@ -97,7 +196,7 @@ class TFRFME { * \param file_name: `string` Name of the file. * \param mode: `string` Format of the file (can be either "HDF5" * or "LGEACY". Default is "LEGACY"). - * \return instance: `TFRFME *` Pointer to anewly created configuration + * \return instance: `TFRFME *` Pointer to a newly created configuration * instance. */ static TFRFME* from_binary(std::string file_name, std::string mode="LEGACY"); diff --git a/src/include/tra_subs.h b/src/include/tra_subs.h index 82594bd0744878f8109bf16003686c2f992a7341..9aece525aa42b28f3ac933317c2b913c95269b6f 100644 --- a/src/include/tra_subs.h +++ b/src/include/tra_subs.h @@ -117,13 +117,13 @@ void ffrt( * \param le: `int` QUESTION: definition? * \param lmode: `int` QUESTION: definition? * \param pmf: `double` QUESTION: definition? - * \param tt1: `fstream &` Handle to first temporary binary file. + * \param tt1: `Swap1 *` Pointer to first swap object. * \param tt2: `fstream &` Handle to second temporary binary file. */ void frfmer( int nkv, double vkm, double *vkv, double vknmx, double apfafa, double tra, double spd, double rir, double ftcn, int le, int lmode, double pmf, - std::fstream &tt1, std::fstream &tt2 + Swap1 *tt1, std::fstream &tt2 ); /*! C++ porting of PWMALP @@ -167,7 +167,7 @@ void sampoa( /*! C++ porting of WAMFF * - * \param wk: Vector of complex. QUESTION: definition? + * \param tt1: `Swap1 *` Pointer to the swap object containing the WK vector. * \param x: `double` * \param y: `double` * \param z: `double` @@ -181,6 +181,6 @@ void sampoa( * \param pmf: `double` QUESTION: definition? */ void wamff( - std::complex<double> *wk, double x, double y, double z, int lm, double apfafa, + Swap1 *tt1, double x, double y, double z, int lm, double apfafa, double tra, double spd, double rir, double ftcn, int lmode, double pmf ); diff --git a/src/libnptm/tfrfme.cpp b/src/libnptm/tfrfme.cpp index 2b649f86e58968167b5fb3deb442fb952c180e65..ac45463bb9ddfcca18ae16f71dcb4d84abf776a3 100644 --- a/src/libnptm/tfrfme.cpp +++ b/src/libnptm/tfrfme.cpp @@ -27,6 +27,153 @@ using namespace std; +Swap1::Swap1(int lm) { + nlmmt = 2 * lm * (lm + 2); + wk = new complex<double>[nlmmt](); +} + +Swap1* Swap1::from_binary(string file_name, string mode) { + Swap1 *instance = NULL; + if (mode.compare("LEGACY") == 0) { + instance = from_legacy(file_name); + } else if (mode.compare("HDF5") == 0) { + instance = from_hdf5(file_name); + } else { + string message = "Unknown format mode: \"" + mode + "\""; + throw UnrecognizedFormatException(message); + } + return instance; +} + +Swap1* Swap1::from_hdf5(string file_name) { + Swap1 *instance = NULL; + unsigned int flags = H5F_ACC_RDONLY; + HDFFile *hdf_file = new HDFFile(file_name, flags); + herr_t status = hdf_file->get_status(); + double *elements; + string str_type; + int _nlmmt, lm, num_elements, index; + complex<double> value; + if (status == 0) { + status = hdf_file->read("NLMMT", "INT32", &_nlmmt); + lm = (int)((-2.0 + sqrt(4.0 + 2.0 * _nlmmt)) / 2.0); + num_elements = 2 * _nlmmt; + instance = new Swap1(lm); + elements = new double[num_elements](); + str_type = "FLOAT64_(" + to_string(num_elements) + ")"; + status = hdf_file->read("WK", str_type, elements); + for (int wi = 0; wi < _nlmmt; wi++) { + index = 2 * wi; + value = complex<double>(elements[index], elements[index + 1]); + instance->set_element(wi, value); + } // wi loop + delete[] elements; + status = hdf_file->close(); + delete hdf_file; + } + return instance; +} + +Swap1* Swap1::from_legacy(string file_name) { + fstream input; + Swap1 *instance = NULL; + int _nlmmt, lm; + double rval, ival; + input.open(file_name.c_str(), ios::in | ios::binary); + if (input.is_open()) { + input.read(reinterpret_cast<char *>(&_nlmmt), sizeof(int)); + lm = (int)((-2.0 + sqrt(4.0 + 2.0 * _nlmmt)) / 2.0); + instance = new Swap1(lm); + for (int j = 0; j < _nlmmt; j++) { + input.read(reinterpret_cast<char *>(&rval), sizeof(double)); + input.read(reinterpret_cast<char *>(&ival), sizeof(double)); + instance->set_element(j, complex<double>(rval, ival)); + } + input.close(); + } else { + printf("ERROR: could not open input file \"%s\"\n", file_name.c_str()); + } + return instance; +} + +void Swap1::write_binary(string file_name, string mode) { + if (mode.compare("LEGACY") == 0) { + write_legacy(file_name); + } else if (mode.compare("HDF5") == 0) { + write_hdf5(file_name); + } else { + string message = "Unknown format mode: \"" + mode + "\""; + throw UnrecognizedFormatException(message); + } +} + +void Swap1::write_hdf5(string file_name) { + List<string> rec_name_list(1); + List<string> rec_type_list(1); + List<void *> rec_ptr_list(1); + herr_t status; + string str_type; + int num_elements = 2 * nlmmt; + rec_name_list.set(0, "NLMMT"); + rec_type_list.set(0, "INT32_(1)"); + rec_ptr_list.set(0, &nlmmt); + rec_name_list.append("WK"); + str_type = "FLOAT64_(" + to_string(num_elements) + ")"; + rec_type_list.append(str_type); + double *ptr_elements = new double[num_elements](); + for (int wi = 0; wi < nlmmt; wi++) { + ptr_elements[2 * wi] = wk[wi].real(); + ptr_elements[2 * wi + 1] = wk[wi].imag(); + } + rec_ptr_list.append(ptr_elements); + + string *rec_names = rec_name_list.to_array(); + string *rec_types = rec_type_list.to_array(); + void **rec_pointers = rec_ptr_list.to_array(); + const int rec_num = rec_name_list.length(); + FileSchema schema(rec_num, rec_types, rec_names); + HDFFile *hdf_file = HDFFile::from_schema(schema, file_name, H5F_ACC_TRUNC); + for (int ri = 0; ri < rec_num; ri++) + status = hdf_file->write(rec_names[ri], rec_types[ri], rec_pointers[ri]); + status = hdf_file->close(); + + delete[] ptr_elements; + delete[] rec_names; + delete[] rec_types; + delete[] rec_pointers; + delete hdf_file; +} + +void Swap1::write_legacy(string file_name) { + fstream output; + double rval, ival; + output.open(file_name.c_str(), ios::out | ios::binary); + if (output.is_open()) { + output.write(reinterpret_cast<char *>(&nlmmt), sizeof(int)); + for (int j = 0; j < nlmmt; j++) { + rval = wk[j].real(); + ival = wk[j].imag(); + output.write(reinterpret_cast<char *>(&rval), sizeof(double)); + output.write(reinterpret_cast<char *>(&ival), sizeof(double)); + } + output.close(); + } else { // Should never occur + printf("ERROR: could not open output file \"%s\"\n", file_name.c_str()); + } +} + +bool Swap1::operator ==(Swap1 &other) { + if (nlmmt != other.nlmmt) { + return false; + } + for (int i = 0; i < nlmmt; i++) { + if (wk[i] != other.wk[i]) { + return false; + } + } + return true; +} + TFRFME::TFRFME(int _lmode, int _lm, int _nkv, int _nxv, int _nyv, int _nzv) { lmode = _lmode; lm = _lm; @@ -67,10 +214,10 @@ TFRFME* TFRFME::from_binary(string file_name, string mode) { instance = from_legacy(file_name); } else if (mode.compare("HDF5") == 0) { instance = from_hdf5(file_name); - } /* else { + } else { string message = "Unknown format mode: \"" + mode + "\""; throw UnrecognizedFormatException(message); - } */ + } return instance; } @@ -202,6 +349,7 @@ TFRFME* TFRFME::from_legacy(string file_name) { instance->set_matrix_element(wi, wj, value); } // wj loop } // wi loop + input.close(); } else { printf("ERROR: could not open input file \"%s\"\n", file_name.c_str()); } @@ -230,8 +378,7 @@ double TFRFME::get_param(string param_name) { else if (param_name.compare("nzv") == 0) value = 1.0 * nzv; else { string message = "Unrecognized parameter name \"" + param_name + "\""; - UnrecognizedParameterException ex(message); - throw ex; + throw UnrecognizedParameterException(message); } return value; } @@ -251,8 +398,7 @@ void TFRFME::set_param(string param_name, double value) { else if (param_name.compare("exril") == 0) exril = value; else { string message = "Unrecognized parameter name \"" + param_name + "\""; - UnrecognizedParameterException ex(message); - throw ex; + throw UnrecognizedParameterException(message); } } @@ -261,10 +407,10 @@ void TFRFME::write_binary(string file_name, string mode) { write_legacy(file_name); } else if (mode.compare("HDF5") == 0) { write_hdf5(file_name); - } /* else { + } else { string message = "Unknown format mode: \"" + mode + "\""; throw UnrecognizedFormatException(message); - } */ + } } void TFRFME::write_hdf5(string file_name) { @@ -392,6 +538,7 @@ void TFRFME::write_legacy(string file_name) { output.write(reinterpret_cast<char *>(&ival), sizeof(double)); } // wj loop } // wi loop + output.close(); } else { // Should never occur printf("ERROR: could not open output file \"%s\"\n", file_name.c_str()); } diff --git a/src/libnptm/tra_subs.cpp b/src/libnptm/tra_subs.cpp index aef2de470a5c3fe0b68aa95f05c0dbdac9d63615..e7a5287315f0c4b651a180e0f483f6ee783f841e 100644 --- a/src/libnptm/tra_subs.cpp +++ b/src/libnptm/tra_subs.cpp @@ -14,6 +14,10 @@ #include "../include/sph_subs.h" #endif +#ifndef INCLUDE_TFRFME_H_ +#include "../include/tfrfme.h" +#endif + #ifndef INCLUDE_TRA_SUBS_H_ #include "../include/tra_subs.h" #endif @@ -243,11 +247,10 @@ void ffrt( void frfmer( int nkv, double vkm, double *vkv, double vknmx, double apfafa, double tra, double spd, double rir, double ftcn, int le, int lmode, double pmf, - std::fstream &tt1, std::fstream &tt2 + Swap1 *tt1, std::fstream &tt2 ) { const int nlemt = le * (le + 2) * 2; const complex<double> cc0(0.0, 0.0); - complex<double> *wk = new complex<double>[nlemt](); double sq = vkm * vkm; for (int jy90 = 0; jy90 < nkv; jy90++) { double vky = vkv[jy90]; @@ -259,27 +262,27 @@ void frfmer( double vkn = sqrt(sqn); if (vkn <= vknmx) { double vkz = sqrt(sq - sqn); - wamff(wk, vkx, vky, vkz, le, apfafa, tra, spd, rir, ftcn, lmode, pmf); - for (int j = 0; j < nlemt; j++) { + wamff(tt1, vkx, vky, vkz, le, apfafa, tra, spd, rir, ftcn, lmode, pmf); + /*for (int j = 0; j < nlemt; j++) { double vreal = wk[j].real(); double vimag = wk[j].imag(); tt1.write(reinterpret_cast<char *>(&vreal), sizeof(double)); tt1.write(reinterpret_cast<char *>(&vimag), sizeof(double)); - } + }*/ tt2.write(reinterpret_cast<char *>(&vkz), sizeof(double)); } else { // label 50 for (int j = 0; j < nlemt; j++) { - double vreal = 0.0; - double vimag = 0.0; - tt1.write(reinterpret_cast<char *>(&vreal), sizeof(double)); - tt1.write(reinterpret_cast<char *>(&vimag), sizeof(double)); + //double vreal = 0.0; + //double vimag = 0.0; + //tt1.write(reinterpret_cast<char *>(&vreal), sizeof(double)); + //tt1.write(reinterpret_cast<char *>(&vimag), sizeof(double)); + tt1->set_element(j, cc0); } double vkz = 0.0; tt2.write(reinterpret_cast<char *>(&vkz), sizeof(double)); } } // jx80 loop } // jy90 loop - delete[] wk; } void pwmalp(std::complex<double> **w, double *up, double *un, std::complex<double> *ylm, int lw) { @@ -364,7 +367,7 @@ void sampoa( } void wamff( - std::complex<double> *wk, double x, double y, double z, int lm, double apfafa, + Swap1 *tt1, double x, double y, double z, int lm, double apfafa, double tra, double spd, double rir, double ftcn, int lmode, double pmf ) { const int nlmm = lm * (lm + 2); @@ -459,13 +462,13 @@ void wamff( } else if (lmode == 2) { // label 50 ftc1 = 2.0 * cthl / (cthl * rir + cth); cfam *= (sth * pmf * ftc1); - for (int i52 = 0; i52 < nlmmt; i52++) wk[i52] = w[i52][0] * cfam; + for (int i52 = 0; i52 < nlmmt; i52++) tt1->set_element(i52, w[i52][0] * cfam); // returns skip62 = true; } else if (lmode == 3) { // label 53 ftc2 = 2.0 * cthl / (cthl + cth * rir); cfam *= (sth * pmf * ftc2); - for (int i55 = 0; i55 < nlmmt; i55++) wk[i55] = w[i55][1] * cfam; + for (int i55 = 0; i55 < nlmmt; i55++) tt1->set_element(i55, w[i55][1] * cfam); // returns skip62 = true; } @@ -480,7 +483,7 @@ void wamff( if (lmode == 0) { double f1 = cph * fam; double f2 = -sph * fam; - for (int i58 = 0; i58 < nlmmt; i58++) wk[i58] = w[i58][0] * f1 + w[i58][1] * f2; + for (int i58 = 0; i58 < nlmmt; i58++) tt1->set_element(i58, w[i58][0] * f1 + w[i58][1] * f2); // returns skip62 = true; } else if (lmode == 1) { // label 60 @@ -491,19 +494,19 @@ void wamff( skip62 = false; } else if (lmode == 2) { // label 65 fam *= (pmf * sth); - for (int i67 = 0; i67 < nlmmt; i67++) wk[i67] = w[i67][0] * fam; + for (int i67 = 0; i67 < nlmmt; i67++) tt1->set_element(i67, w[i67][0] * fam); // returns skip62 = true; } else if (lmode == 3) { // label 68 fam *= (pmf * sth); - for (int i70 = 0; i70 < nlmmt; i70++) wk[i70] = w[i70][1] * fam; + for (int i70 = 0; i70 < nlmmt; i70++) tt1->set_element(i70, w[i70][1] * fam); // returns skip62 = true; } } if (!skip62) { if (lmode == 0 || lmode == 1) { // label 62 - for (int i63 = 0; i63 < nlmmt; i63++) wk[i63] = w[i63][0] * cf1 + w[i63][1] * cf2; + for (int i63 = 0; i63 < nlmmt; i63++) tt1->set_element(i63, w[i63][0] * cf1 + w[i63][1] * cf2); } } } diff --git a/src/trapping/cfrfme.cpp b/src/trapping/cfrfme.cpp index 29d21e8f18711313c521dc75d980870018b0939e..7e2915171054ea250a6e7cb3ccbedb75798a0671 100644 --- a/src/trapping/cfrfme.cpp +++ b/src/trapping/cfrfme.cpp @@ -42,12 +42,12 @@ using namespace std; */ void frfme(string data_file, string output_path) { string tfrfme_name = output_path + "/c_TFRFME.hd5"; - //fstream tfrfme; TFRFME *tfrfme = NULL; + Swap1 *tt1 = NULL; char namef[7]; char more; double *vkv = NULL, **vkzm = NULL; - complex<double> *wk = NULL, **w = NULL; + complex<double> **w = NULL; const complex<double> cc0(0.0, 0.0); const complex<double> uim(0.0, 1.0); int line_count = 0, last_read_line = 0; @@ -272,14 +272,16 @@ void frfme(string data_file, string output_path) { tfrfme->set_param("frsh", frsh); tfrfme->set_param("exril", exril); fstream temptape1, temptape2; - string temp_name1 = output_path + "/c_TEMPTAPE1"; + string temp_name1 = output_path + "/c_TEMPTAPE1.hd5"; string temp_name2 = output_path + "/c_TEMPTAPE2"; - temptape1.open(temp_name1.c_str(), ios::out | ios::binary); + //temptape1.open(temp_name1.c_str(), ios::out | ios::binary); + tt1 = new Swap1(lm); temptape2.open(temp_name2.c_str(), ios::out | ios::binary); for (int jx = 0; jx < nkv; jx++) temptape2.write(reinterpret_cast<char *>(&(vkv[jx])), sizeof(double)); - frfmer(nkv, vkm, vkv, vknmx, apfafa, tra, spd, rir, ftcn, lm, lmode, pmf, temptape1, temptape2); - temptape1.close(); + frfmer(nkv, vkm, vkv, vknmx, apfafa, tra, spd, rir, ftcn, lm, lmode, pmf, tt1, temptape2); + //temptape1.close(); + tt1->write_binary(temp_name1, "HDF5"); temptape2.write(reinterpret_cast<char *>(&apfafa), sizeof(double)); temptape2.write(reinterpret_cast<char *>(&pmf), sizeof(double)); temptape2.write(reinterpret_cast<char *>(&spd), sizeof(double)); @@ -316,21 +318,21 @@ void frfme(string data_file, string output_path) { } w = new complex<double>*[nkv]; for (int wi = 0; wi < nkv; wi++) w[wi] = new complex<double>[nkv](); - if (wk != NULL) delete[] wk; - wk = new complex<double>[nlmmt](); - temptape1.open(temp_name1.c_str(), ios::in | ios::binary); + //if (wk != NULL) delete[] wk; + //wk = new complex<double>[nlmmt](); + //temptape1.open(temp_name1.c_str(), ios::in | ios::binary); for (int jy50 = 0; jy50 < nkv; jy50++) { for (int jx50 = 0; jx50 < nkv; jx50++) { - for (int i = 0; i < nlmmt; i++) { + /*for (int i = 0; i < nlmmt; i++) { double vreal, vimag; temptape1.read(reinterpret_cast<char *>(&vreal), sizeof(double)); temptape1.read(reinterpret_cast<char *>(&vimag), sizeof(double)); wk[i] = complex<double>(vreal, vimag); - } - w[jx50][jy50] = wk[j80]; + }*/ + w[jx50][jy50] = tt1->get_element(j80); } // jx50 } // jy50 loop - temptape1.close(); + //temptape1.close(); int ixyz = 0; for (int wj = 0; wj < nrvc; wj++) tfrfme->set_matrix_element(j80, wj, cc0); for (int iz75 = 0; iz75 < nzv; iz75++) { @@ -413,6 +415,7 @@ void frfme(string data_file, string output_path) { for (int wi = nkv - 1; wi > -1; wi--) delete[] w[wi]; delete[] w; } - if (wk != NULL) delete[] wk; + //if (wk != NULL) delete[] wk; + if (tt1 != NULL) delete tt1; printf("Done.\n"); }