Adding initial arldatastream with ARLIndexHeader, ARLHeader and ARLRecord.

radixio/CMakeLists.txt

@@ -2,6 +2,7 @@ TRIBITS_SUBPACKAGE(io)
 
 
 SET(SOURCE
+   arldatastream.cc
    csvfile.cc
    cfgfile.cc
    eafstream.cc
@@ -15,6 +16,7 @@ SET(SOURCE
 )
 
 SET(HEADERS
+   arldatastream.hh
    csvfile.hh
    cfgfile.hh
    eafstream.hh

radixio/arldatastream.cc

+#include "radixio/arldatastream.hh"
+
+#include "radixbug/bug.hh"
+#include "radixcore/stringfunctions.hh"
+
+#include <fstream>
+namespace radix
+{
+class ARLDataStream::PImpl
+{
+ public:
+  std::string file;
+  std::shared_ptr<std::fstream> stream;
+};
+
+void ARLDataStream::expand(const std::string& val, ARLIndexHeader& index)
+{
+  if (val.size() < 50)
+  {
+    throw std::runtime_error(
+        "Incorrect size for string expandsion to ARLIndexHeader.");
+  }
+  index.year  = std::atoi(val.substr(0, 2).c_str());
+  index.month = std::atoi(val.substr(2, 2).c_str());
+  index.day   = std::atoi(val.substr(4, 2).c_str());
+  index.hour  = std::atoi(val.substr(6, 2).c_str());
+  index.ic    = std::atoi(val.substr(8, 2).c_str());
+  index.il    = std::atoi(val.substr(10, 2).c_str());
+  index.cgrid = val.substr(12, 2).c_str();
+  index.kvar  = val.substr(14, 4);
+  index.nexp  = std::atoi(val.substr(18, 4).c_str());
+  index.prec  = float(std::atof(val.substr(22, 14).c_str()));
+  index.var1  = float(std::atof(val.substr(36, 14).c_str()));
+}
+void ARLDataStream::expand(const std::string& val, const ARLIndexHeader& index,
+                           ARLHeader& header)
+{
+  if (val.size() < 108)
+  {
+    throw std::runtime_error(
+        "Incorrect size for string expansion to ARLHeader.");
+  }
+  header.model_id = val.substr(0, 4);
+  header.icx      = std::atoi(val.substr(4, 3).c_str());
+  header.mn       = std::atoi(val.substr(7, 2).c_str());
+  header.pole_lat = float(std::atof(val.substr(9, 7).c_str()));
+  header.pole_lon = float(std::atof(val.substr(16, 7).c_str()));
+  header.ref_lat  = float(std::atof(val.substr(23, 7).c_str()));
+  header.ref_lon  = float(std::atof(val.substr(30, 7).c_str()));
+  header.size     = float(std::atof(val.substr(37, 7).c_str()));
+  header.orient   = float(std::atof(val.substr(44, 7).c_str()));
+  header.tang_lat = float(std::atof(val.substr(51, 7).c_str()));
+  header.sync_xp  = float(std::atof(val.substr(58, 7).c_str()));
+  header.sync_yp  = float(std::atof(val.substr(65, 7).c_str()));
+  header.sync_lat = float(std::atof(val.substr(72, 7).c_str()));
+  header.sync_lon = float(std::atof(val.substr(79, 7).c_str()));
+  header.dummy    = float(std::atof(val.substr(86, 7).c_str()));
+  header.nx       = std::atoi(val.substr(93, 3).c_str());
+  header.ny       = std::atoi(val.substr(96, 3).c_str());
+  header.nz       = std::atoi(val.substr(99, 3).c_str());
+  header.z_flag   = std::atoi(val.substr(102, 2).c_str());
+  header.lenh     = std::atoi(val.substr(104, 4).c_str());
+
+  int knx = ordinal(index.cgrid[0]);
+  int kny = ordinal(index.cgrid[1]);
+  // Check for the grid domain extending beyond 3 digits
+  if (knx >= 64 || kny >= 64)
+  {
+    header.nx = (knx - 64) * 1000 + header.nx;
+    header.ny = (kny - 64) * 1000 + header.ny;
+  }
+}
+
+ARLDataStream::ARLDataStream(const std::string& file)
+    : p(new PImpl(), [](PImpl* impl) { delete impl; })
+{
+  p->file = file;
+}
+
+/**
+ * @brief The ARLRecord::PImpl class
+ */
+class ARLRecord::PImpl
+{
+ public:
+};
+ARLRecord::ARLRecord()
+    : p(new ARLRecord::PImpl(), [](ARLRecord::PImpl* impl) { delete impl; })
+{
+}
+
+ARLRecord::ARLRecord(const ARLRecord& orig)
+    : p(new ARLRecord::PImpl(), [](ARLRecord::PImpl* impl) { delete impl; })
+{
+}
+
+}  // namespace radix

radixio/arldatastream.hh

+#ifndef RADIX_RADIXIO_ARLDATASTREAM_HH_
+#define RADIX_RADIXIO_ARLDATASTREAM_HH_
+
+#include <algorithm>
+#include <cmath>
+#include <iomanip>
+#include <memory>
+#include <sstream>
+#include <string>
+#include <vector>
+
+#include "radixcore/visibility.hh"
+
+namespace radix
+{
+//
+// Forward declare ARLRecord and ARLHeader
+struct ARLIndexHeader;
+struct ARLHeader;
+class ARLRecord;
+
+class RADIX_PUBLIC ARLDataStream
+{
+  // forward declare private implementation
+  class PImpl;
+  // unique pointer to private implmentation
+  std::unique_ptr<PImpl, void (*)(PImpl*)> p;
+
+ private:
+  /**
+   * @brief expand Expands string of 50 characters to an ARLIndexHeader
+   * @param val String for expansion
+   * @param index ARLIndexHeader
+   */
+  static void expand(const std::string& val, ARLIndexHeader& index);
+
+  /**
+   * @brief expand Expands string of 108 characters to ARLHeader
+   * @param val String for expansion
+   * @param index ALRIndexHeader required for determining grid domain extending
+   * beyond 3 digits
+   * @param header ARLHeader
+   */
+  static void expand(const std::string& val, const ARLIndexHeader& index,
+                     ARLHeader& header);
+
+ public:
+  ARLDataStream() = delete;
+  ARLDataStream(const std::string& file);
+
+  /**
+   * @brief read_index_header Reads a single index header from the stream
+   * @param iheader ARLIndexHeader&
+   * @return
+   */
+  bool read_index_header(ARLIndexHeader& iheader);
+
+  bool read_record_header(ARLHeader& rheader);
+  /**
+   * @brief read_record Reads a single record from the stream
+   * @param ARLRecord& record
+   *
+   * @return bool on if inventory was populated
+   */
+  bool read_record(ARLRecord& record);
+
+  // TODO
+  // bool write_record(const ARLRecord& record);
+
+};  // class ARLDataStream
+
+struct RADIX_PUBLIC ARLHeader
+{
+ public:
+  bool prime;
+  bool latlon;
+  bool gbldat;
+  bool global;
+  int icx;
+  int mn;
+  float pole_lat;
+  float pole_lon;
+  float ref_lon;
+  float ref_lat;
+  float size;
+  float orient;
+  float tang_lat;
+  float sync_xp;
+  float sync_yp;
+  float sync_lat;
+  float sync_lon;
+  float dummy;
+  int nx;
+  int ny;
+  int nz;
+  int z_flag;
+  int lenh;
+  std::string model_id;
+};  // ARLHeader
+
+struct RADIX_PUBLIC ARLIndexHeader
+{
+ public:
+  int year;
+  int month;
+  int day;
+  int hour;
+  int ic;
+  int il;
+  int nexp;
+  float prec;
+  float var1;
+  std::string cgrid;
+  std::string kvar;
+};
+
+class RADIX_PUBLIC ARLRecord
+{
+  // forward declare private implementation
+  class PImpl;
+  // unique pointer to private implmentation
+  std::unique_ptr<PImpl, void (*)(PImpl*)> p;
+
+ public:
+  ARLRecord();
+  ARLRecord(const ARLRecord& orig);
+
+};  // ARLRecord
+}  // namespace radix
+
+#endif /** RADIX_RADIXIO_ARLDATASTREAM_HH_ */

radixio/arldatastream.i.hh

+
+#ifndef RADIX_RADIXIO_F71STREAM_I_HH_
+#define RADIX_RADIXIO_F71STREAM_I_HH_
+
+#include "radixio/f71stream.hh"
+
+#include "radixbug/bug.hh"
+#include "radixcore/stringfunctions.hh"
+#include "radixio/eafstream.hh"
+
+namespace radix
+{
+template <typename f71case_type>
+bool F71Stream::read_case(f71case_type& inventory)
+{
+  //
+  // Load the inventory with processed data
+  //  ITXX - total number of nuclides or ‘0’to signal an end-of-file
+  //  ILITE - number of nuclides in light element library - 0 if not using
+  //  IACT - number of nuclides in actinide library - 0 if not using
+  //  IFP - number of nuclides in fission product library - 0 if not using
+  //  NRFLAG - flag that determines position numbers for data set to be saved
+  //  MSUB1 - time step number of the data to be saved
+  //  NSTEPO - position of requested data in concentration file
+  //  KASEPO - position of the requested case in the concentration file
+  //  JOBPOS - position in the file of the job containing the requested data
+  //  NOCS - subcase number in the case containing the requested data
+  //  NOBLND - 0/1/N – job termination/no blending/N streams blended
+  //  NDSET - nuclear data library unit number
+  //  NTYPE - the type of library on NDSET, 0/>0 card-image / binary library
+  //  NGRP - the number of neutron energy groups
+  //  NELEM - maximum atomic number of nuclides in library
+  //  NVERT - −N/0/N – request integral option and writing N/not wanted/ reading
+  //  N NG -the number of photon energy groups in the master photon library MMN
+  //  -  number of irradiation time-step intervals in a subcase MOUT - controls
+  //  writing concentrations on unit NXTR
+  //          If MOUT=0 and MSUB=m, concentration for time-step m is written.
+  //          If MOUT= −N and MSUB>= 0, concentrations for time-steps MSUB
+  //          through N, inclusive, are written.
+  //  INDEX - 0/1 – power read in 58* array/ flux read in 59* array
+  //  MSTAR - time step at which concentration print cutoffs are applied
+  //  NUNIT - 1/2/3/4/5/6 – time (in 60*) units s/min/h/d/y/other
+  //  KBLEND - time-step number indicating when blending is to occur
+  //  NENLE - number of groups in light element photon library
+  //  NENAC - number of groups in actinide photon library
+  //  NENFP - number of groups in fission product photon library
+  //  L1 - max0(1, NENLE)
+  //  L2 - max0(1, NENAC)
+  //  L3 - max0(1, NANFP)
+  //  L4 - max0(1, IACT)
+  //  TMO - time at start of first time-step
+  //  FRACPW - ratio of subcase irradiation time to total case irradiation time
+  //  TCONST - If NUNIT=6, number of seconds in the special time unit.
+  //  TUNIT - time unit printed in output (4-character variable)
+  //  TWRITE - value of time-step accumulative time of concentrations
+  //           written in file
+  //  PWRITE - value of time-step power of concentrations written in file
+  //  FWRITE - value of time-step flux of concentrations written in the file
+
+  //  The second record type description follows:
+
+  //  NUCL(ITXX) - ITXX values of nuclide ID’s as ZAS
+  //  (Z=atomic number, A=mass number, S=metastable state, S=0 for ground state)
+  //  X(ITXX) - ITXX values of gram-atom concentrations, at step requested MSUB.
+  //  TITLE(20) - 80 character title for this subcase if NTI=1 or 3
+  //  BASIS(10) - 40 character title for basis of calculation
+  //              - may be blank if NTI=2 or 3
+  //  EGAMLE(NENLE+1)- Gamma group boundaries for light elements
+  //  EGAMAC(NENAC+1)- Gamma group boundaries for actinides
+  //  EGAMFP(NENFP +1)- Gamma group boundaries for fission products
+  //  SPECLE(L1) - L1 values of gamma source spectrum for
+  //               light elements in photons/s
+  //  SPECAC(L2) - L2 values of gamma source spectrum for actinides in photons/s
+  //  SPECFP(L3) - L3 values of gamma source spectrum for
+  //               fission products in photons/s
+  //  ENER(NG+1) - NG +1
+  //  values of gamma energy group boundaries on master photon library
+  //  DSAV(NG) -  NG values of master photon spectrum in photons/s
+  //  ESAV(NG) - NG values of master photon spectrum in MeV/s
+  //  ITSAV(20) - 80-character photon spectrum title
+  //  SPNNUC(L4) - L4 values of spontaneous fission neutron source in neutrons/s
+  //  ALPNUC(L4) - L4 values of actinide nuclide (alpha,n) source in neutrons/s
+  //  SPNEUT(NGRP) - NGRP values of total neutron spectrum in neutrons/s
+  //  SPECAN(NGRP) - NGRP values of (alpha,n) spectrum in neutrons/s
+  //  SPECSP(NGRP) - NGRP values of spontaneous fission neutron
+  //  spectrum in neutrons/s
+  //  ENEUTS(NGRP+1) - Neutron energy-group boundaries
+
+  // check if the file is open
+  if (!mStream.is_open())
+  {
+    mStream.open(mFile.c_str(), std::ios::binary | std::ios::in);
+  }
+  if (!mStream.is_open())
+  {
+    throw std::runtime_error("Failed to open file: '" + mFile + "'");
+  }
+  if (!mStream.good()) return false;
+
+  // header
+  mStream.readHeader();
+  // read record content
+  int num_nuclides = mStream.readInt();  // ITXX
+  if (num_nuclides == 0) return false;   // EOF
+
+  int num_light_nuclides   = mStream.readInt();        // ILITE
+  int num_actinides        = mStream.readInt();        // IACT   =
+  int num_fission_products = mStream.readInt();        // IFP    =
+  mStream.readInt();                                   // NRFLAG =
+  mStream.readInt();                                   // MSUB1  =
+  mStream.readInt();                                   // NSTEPO =
+  mStream.readInt();                                   // KASEPO =
+  mStream.readInt();                                   // JOBPOS =
+  mStream.readInt();                                   // NOCS   =
+  mStream.readInt();                                   // NOBLND =
+  mStream.readInt();                                   // NDSET  =
+  mStream.readInt();                                   // NTYPE  =
+  int ngrp = mStream.readInt();                        // NGRP   =
+  mStream.readInt();                                   // NELEM  =
+  mStream.readInt();                                   // NVERT  =
+  int ng = mStream.readInt();                          // NG     =
+  mStream.readInt();                                   // MMN    =
+  mStream.readInt();                                   // MOUT   =
+  mStream.readInt();                                   // INDEX  =
+  mStream.readInt();                                   // MSTAR  =
+  int nunit = mStream.readInt();                       // NUNIT  =
+  mStream.readInt();                                   // KBLEND =
+  int num_light_element_groups   = mStream.readInt();  // NENLE  =
+  int num_actinide_groups        = mStream.readInt();  // NENAC  =
+  int num_fission_product_groups = mStream.readInt();  // NENFP  =
+  int l1                         = mStream.readInt();  // L1     =
+  int l2                         = mStream.readInt();  // L2     =
+  int l3                         = mStream.readInt();  // L3     =
+  int l4                         = mStream.readInt();  // L4     =
+  mStream.readFloat();                                 // TMO    =
+  mStream.readFloat();                                 // FRACPW =
+  mStream.readFloat();                                 // TCONST =
+  mStream.readString(4);                               // TUNIT  =
+  float twrite = mStream.readFloat();                  // TWRITE =
+  float pwrite = mStream.readFloat();                  // PWRITE =
+  float fwrite = mStream.readFloat();                  // FWRITE =
+
+  // read end of record
+  mStream.readFooter();
+  radix_line("Total nuclides " << num_nuclides << "=[" << num_light_nuclides
+                               << ", " << num_actinides << ", "
+                               << num_fission_products << "]");
+  radix_line("Time: " << twrite << " " << nunit << " power(" << pwrite
+                      << ") flux(" << fwrite << ")");
+
+  // read record 2: the nuclides and their amounts
+  // read header
+  mStream.readHeader();
+  // break iza_list into light/actinide/fission product nuclides lists
+  std::vector<int> light_nuclides((size_t(num_light_nuclides)));
+  std::vector<int> actinide_nuclides((size_t(num_actinides)));
+  std::vector<int> fission_product_nuclides((size_t(num_fission_products)));
+  mStream >> light_nuclides >> actinide_nuclides >> fission_product_nuclides;
+
+  // break abundances by light/actinide/fission product
+  std::vector<float> light_abs((size_t(num_light_nuclides)));
+  std::vector<float> actinide_abs((size_t(num_actinides)));
+  std::vector<float> fission_product_abs((size_t(num_fission_products)));
+  mStream >> light_abs >> actinide_abs >> fission_product_abs;
+
+  std::string title = trim_string(mStream.readString(80));
+  std::string basis = trim_string(mStream.readString(40));
+  radix_line("Title(" << title << ") Basis(" << basis << ")");
+
+  std::vector<float> gamma_le_bounds((size_t(num_light_element_groups)));
+  mStream >> gamma_le_bounds;
+
+  std::vector<float> gamma_ae_bounds((size_t(num_actinide_groups)));
+  mStream >> gamma_ae_bounds;
+
+  std::vector<float> gamma_fpe_bounds((size_t(num_fission_product_groups)));
+  mStream >> gamma_fpe_bounds;
+
+  std::vector<float> spectrum_le((size_t(l1)));
+  mStream >> spectrum_le;
+
+  std::vector<float> spectrum_ac((size_t(l2)));
+  mStream >> spectrum_ac;
+
+  std::vector<float> spectrum_fp((size_t(l3)));
+  mStream >> spectrum_fp;
+
+  std::vector<float> energy_bounds((size_t(ng + 1)));
+  mStream >> energy_bounds;
+
+  std::vector<float> dsav((size_t(ng)));
+  mStream >> dsav;
+
+  std::vector<float> esav((size_t(ng)));
+  mStream >> esav;
+
+  // photon spectrum title
+  std::string itsav = trim_string(mStream.readString(80));
+
+  // values of total neutron spectrum in neutron/s
+  std::vector<float> spnnuc((size_t(l4)));
+  mStream >> spnnuc;
+
+  std::vector<float> alpnuc((size_t(l4)));
+  mStream >> alpnuc;
+
+  // values of total neutron spectrum in neutrons/s
+  std::vector<float> spneut((size_t(ngrp)));
+  mStream >> spneut;
+
+  // values of (alpha,n) spectrum in neutrons/s
+  std::vector<float> specan((size_t(ngrp)));
+  mStream >> specan;
+
+  // values of spontaneous fission neutron spectrum in neutrons/s
+  std::vector<float> specsp((size_t(ngrp)));
+  mStream >> specsp;
+
+  std::vector<float> eneuts((size_t(ngrp + 1)));
+  mStream >> eneuts;
+
+  // read end of record
+  mStream.readFooter();
+
+  mStream.peek();  // peek into the next bit. This ensures EOF is seen in good()
+                   // call
+  inventory.setLightElementNuclides(light_nuclides);
+  inventory.setLightElementNuclideAbundances(light_abs);
+  inventory.setActinideNuclides(actinide_nuclides);
+  inventory.setActinideNuclideAbundances(actinide_abs);
+  inventory.setFissionProductNuclides(fission_product_nuclides);
+  inventory.setFissionProductNuclideAbundances(fission_product_abs);
+
+  switch (nunit)
+  {
+    case 1:  // seconds
+      inventory.setTime(twrite);
+      break;
+    case 2:  // minutes
+      inventory.setTime(twrite * 60.);
+      break;
+    case 3:  // hours
+      inventory.setTime(twrite * 3600.);
+      break;
+    case 4:  // days
+      inventory.setTime(twrite * 86400);
+      break;
+    case 5:  // years
+      inventory.setTime(twrite * 3.154e+7);
+      break;
+  }
+  inventory.setPower(pwrite);
+  inventory.setFlux(fwrite);
+  inventory.setTitle(title);
+  inventory.setBasis(basis);
+  inventory.setLightElementGammaEnergyBounds(gamma_le_bounds);
+  inventory.setActinideGammaEnergyBounds(gamma_ae_bounds);
+  inventory.setFissionProductGammaEnergyBounds(gamma_fpe_bounds);
+  inventory.setLightElementSpectrum(spectrum_le);
+  inventory.setActinideSpectrum(spectrum_ac);
+  inventory.setFissionProductSpectrum(spectrum_fp);
+  inventory.setGammaEnergyBounds(energy_bounds);
+  inventory.setPhotonPerSSpectrum(dsav);
+  inventory.setMevPerSSpectrum(esav);
+  inventory.setPhotonSpectrumTitle(itsav);
+  inventory.setSpontaneousFissionNeutronSource(spnnuc);
+  inventory.setActinideAlphaNSource(alpnuc);
+  inventory.setTotalNeutronSpectrum(spneut);
+  inventory.setAlphaNSpectrum(specan);
+  inventory.setSpontaneousFissionSpectrum(specsp);
+  inventory.setNeutronEnergyBounds(eneuts);
+  return true;
+}
+}  // namespace radix
+
+#endif /** RADIX_RADIXIO_F71STREAM_I_HH_ */