Commit 07097559 authored by Giovanni La Mura's avatar Giovanni La Mura
Browse files

Divert TEMPTAPE1 swapping to RAM

parent 927f93b7

src/include/tfrfme.h

+105 −6
Original line number Diff line number Diff line
@@ -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:

  /*! \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

src/include/tra_subs.h

+4 −4
Original line number Diff line number Diff line
@@ -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

);

src/libnptm/tfrfme.cpp

+155 −8
Original line number Diff line number Diff line
@@ -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());

  }

src/libnptm/tra_subs.cpp

+20 −17
Original line number Diff line number Diff line
@@ -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);

      }

    }

  }

src/trapping/cfrfme.cpp

+17 −14
Original line number Diff line number Diff line
@@ -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");

}