thermodynamiclatticeoutput.cc

Go to the documentation of this file.

/*
 *
 *    Copyright (c) 2014-2022
 *      SMASH Team
 *
 *    GNU General Public License (GPLv3 or later)
 *
 */
 
#include "smash/thermodynamiclatticeoutput.h"
 
#include <fstream>
#include <memory>
 
#include <boost/filesystem.hpp>
 
#include "smash/clock.h"
#include "smash/config.h"
#include "smash/density.h"
#include "smash/energymomentumtensor.h"
#include "smash/experimentparameters.h"
#include "smash/thermodynamicoutput.h"
#include "smash/vtkoutput.h"
 
namespace smash {
 
/* initialization of the static member version */
const double_t ThermodynamicLatticeOutput::version = 1.0;
 
ThermodynamicLatticeOutput::ThermodynamicLatticeOutput(
    const bf::path &path, const std::string &name,
    const OutputParameters &out_par, const bool enable_ascii,
    const bool enable_binary)
    : OutputInterface(name),
      out_par_(out_par),
      base_path_(std::move(path)),
      enable_ascii_(enable_ascii),
      enable_binary_(enable_binary) {
  if (enable_ascii_ || enable_binary_) {
    enable_output_ = true;
  } else {
    enable_output_ = false;
  }
  if (enable_ascii_) {
    if (out_par_.td_rho_eckart) {
      output_ascii_files_[ThermodynamicQuantity::EckartDensity] =
          std::make_shared<std::ofstream>(nullptr);
    }
    if (out_par_.td_tmn_landau) {
      output_ascii_files_[ThermodynamicQuantity::TmnLandau] =
          std::make_shared<std::ofstream>(nullptr);
    }
    if (out_par_.td_tmn) {
      output_ascii_files_[ThermodynamicQuantity::Tmn] =
          std::make_shared<std::ofstream>(nullptr);
    }
    if (out_par_.td_v_landau) {
      output_ascii_files_[ThermodynamicQuantity::LandauVelocity] =
          std::make_shared<std::ofstream>(nullptr);
    }
    if (out_par_.td_jQBS) {
      output_ascii_files_[ThermodynamicQuantity::j_QBS] =
          std::make_shared<std::ofstream>(nullptr);
    }
  }
  if (enable_binary_) {
    if (out_par_.td_rho_eckart) {
      output_binary_files_[ThermodynamicQuantity::EckartDensity] =
          std::make_shared<std::ofstream>(nullptr);
    }
    if (out_par_.td_tmn_landau) {
      output_binary_files_[ThermodynamicQuantity::TmnLandau] =
          std::make_shared<std::ofstream>(nullptr);
    }
    if (out_par_.td_tmn) {
      output_binary_files_[ThermodynamicQuantity::Tmn] =
          std::make_shared<std::ofstream>(nullptr);
    }
    if (out_par_.td_v_landau) {
      output_binary_files_[ThermodynamicQuantity::LandauVelocity] =
          std::make_shared<std::ofstream>(nullptr);
    }
    if (out_par_.td_jQBS) {
      output_binary_files_[ThermodynamicQuantity::j_QBS] =
          std::make_shared<std::ofstream>(nullptr);
    }
  }
}
 
ThermodynamicLatticeOutput::~ThermodynamicLatticeOutput() {}
 
void ThermodynamicLatticeOutput::at_eventstart(
    const int event_number, const ThermodynamicQuantity tq,
    const DensityType dens_type, RectangularLattice<DensityOnLattice> lattice) {
  if (!enable_output_) {
    return;
  }
  assert((tq == ThermodynamicQuantity::EckartDensity) ||
         (tq == ThermodynamicQuantity::j_QBS));
  // at the next refactoring of the code,
  // this piece should go in the constructor
  const auto dim = lattice.n_cells();
  const auto cs = lattice.cell_sizes();
  const auto orig = lattice.origin();
  for (int l = 0; l < 3; l++) {
    nodes_[l] = dim[l];
    sizes_[l] = cs[l];
    origin_[l] = orig[l];
  }
  std::shared_ptr<std::ofstream> fp(nullptr);
  std::string varname;
  std::string filename;
  varname = make_varname(tq, dens_type);
  if (tq == ThermodynamicQuantity::EckartDensity) {
    if (enable_ascii_) {
      filename = make_filename(varname, event_number, 'a');
      try {
        output_ascii_files_[ThermodynamicQuantity::EckartDensity]->open(
            filename, std::ios::out);
      } catch (std::ofstream::failure &e) {
        logg[LogArea::Main::id].fatal()
            << "Error in opening " << filename << std::endl;
        throw std::runtime_error(
            "Not possible to write thermodynamic "
            "lattice output to file.");
      }
      fp = output_ascii_files_[ThermodynamicQuantity::EckartDensity];
      write_therm_lattice_ascii_header(fp, tq);
    }
    if (enable_binary_) {
      filename = make_filename(varname, event_number, 'b');
      try {
        output_binary_files_[ThermodynamicQuantity::EckartDensity]->open(
            filename, std::ios::out | std::ios::binary);
      } catch (std::ofstream::failure &e) {
        logg[LogArea::Main::id].fatal()
            << "Error in opening " << filename << std::endl;
        throw std::runtime_error(
            "Not possible to write thermodynamic "
            "lattice output to file.");
      }
      fp = output_binary_files_[ThermodynamicQuantity::EckartDensity];
      write_therm_lattice_binary_header(fp, tq);
    }
  } else {
    if (enable_ascii_) {
      filename = make_filename(varname, event_number, 'a');
      try {
        output_ascii_files_[ThermodynamicQuantity::j_QBS]->open(filename,
                                                                std::ios::out);
      } catch (std::ofstream::failure &e) {
        logg[LogArea::Main::id].fatal()
            << "Error in opening " << filename << std::endl;
        throw std::runtime_error(
            "Not possible to write thermodynamic "
            "lattice output to file.");
      }
      fp = output_ascii_files_[ThermodynamicQuantity::j_QBS];
      write_therm_lattice_ascii_header(fp, tq);
    }
    if (enable_binary_) {
      filename = make_filename(varname, event_number, 'b');
      try {
        output_binary_files_[ThermodynamicQuantity::j_QBS]->open(
            filename, std::ios::out | std::ios::binary);
      } catch (std::ofstream::failure &e) {
        logg[LogArea::Main::id].fatal()
            << "Error in opening " << filename << std::endl;
        throw std::runtime_error(
            "Not possible to write thermodynamic "
            "lattice output to file.");
      }
      fp = output_binary_files_[ThermodynamicQuantity::j_QBS];
      write_therm_lattice_binary_header(fp, tq);
    }
  }
}
 
void ThermodynamicLatticeOutput::at_eventstart(
    const int event_number, const ThermodynamicQuantity tq,
    const DensityType dens_type,
    RectangularLattice<EnergyMomentumTensor> lattice) {
  if (!enable_output_) {
    return;
  }
  const auto dim = lattice.n_cells();
  const auto cs = lattice.cell_sizes();
  const auto orig = lattice.origin();
  for (int l = 0; l < 3; l++) {
    nodes_[l] = dim[l];
    sizes_[l] = cs[l];
    origin_[l] = orig[l];
  }
  std::shared_ptr<std::ofstream> fp(nullptr);
  std::string varname;
  std::string filename;
  varname = make_varname(tq, dens_type);
  if (enable_ascii_) {
    filename = make_filename(varname, event_number, 'a');
    if (tq == ThermodynamicQuantity::TmnLandau) {
      try {
        output_ascii_files_[ThermodynamicQuantity::TmnLandau]->open(
            filename, std::ios::out);
      } catch (std::ofstream::failure &e) {
        logg[LogArea::Main::id].fatal()
            << "Error in opening " << filename << std::endl;
        throw std::runtime_error(
            "Not possible to write thermodynamic "
            "lattice output to file.");
      }
      fp = output_ascii_files_[ThermodynamicQuantity::TmnLandau];
    } else if (tq == ThermodynamicQuantity::Tmn) {
      try {
        output_ascii_files_[ThermodynamicQuantity::Tmn]->open(filename,
                                                              std::ios::out);
      } catch (std::ofstream::failure &e) {
        logg[LogArea::Main::id].fatal()
            << "Error in opening " << filename << std::endl;
        throw std::runtime_error(
            "Not possible to write thermodynamic "
            "lattice output to file.");
      }
      fp = output_ascii_files_[ThermodynamicQuantity::Tmn];
    } else if (tq == ThermodynamicQuantity::LandauVelocity) {
      try {
        output_ascii_files_[ThermodynamicQuantity::LandauVelocity]->open(
            filename, std::ios::out);
      } catch (std::ofstream::failure &e) {
        logg[LogArea::Main::id].fatal()
            << "Error in opening " << filename << std::endl;
        throw std::runtime_error(
            "Not possible to write thermodynamic "
            "lattice output to file.");
      }
      fp = output_ascii_files_[ThermodynamicQuantity::LandauVelocity];
    } else {
      try {
        output_ascii_files_[ThermodynamicQuantity::j_QBS]->open(filename,
                                                                std::ios::out);
      } catch (std::ofstream::failure &e) {
        logg[LogArea::Main::id].fatal()
            << "Error in opening " << filename << std::endl;
        throw std::runtime_error(
            "Not possible to write thermodynamic "
            "lattice output to file.");
      }
    }
    write_therm_lattice_ascii_header(fp, tq);
  }
  if (enable_binary_) {
    filename = make_filename(varname, event_number, 'b');
    if (tq == ThermodynamicQuantity::TmnLandau) {
      try {
        output_binary_files_[ThermodynamicQuantity::TmnLandau]->open(
            filename, std::ios::out | std::ios::binary);
      } catch (std::ofstream::failure &e) {
        logg[LogArea::Main::id].fatal()
            << "Error in opening " << filename << std::endl;
        throw std::runtime_error(
            "Not possible to write thermodynamic "
            "lattice output to file.");
      }
      fp = output_binary_files_[ThermodynamicQuantity::TmnLandau];
    } else if (tq == ThermodynamicQuantity::Tmn) {
      try {
        output_binary_files_[ThermodynamicQuantity::Tmn]->open(
            filename, std::ios::out | std::ios::binary);
      } catch (std::ofstream::failure &e) {
        logg[LogArea::Main::id].fatal()
            << "Error in opening " << filename << std::endl;
        throw std::runtime_error(
            "Not possible to write thermodynamic "
            "lattice output to file.");
      }
      fp = output_binary_files_[ThermodynamicQuantity::Tmn];
    } else if (tq == ThermodynamicQuantity::LandauVelocity) {
      try {
        output_binary_files_[ThermodynamicQuantity::LandauVelocity]->open(
            filename, std::ios::out | std::ios::binary);
      } catch (std::ofstream::failure &e) {
        logg[LogArea::Main::id].fatal()
            << "Error in opening " << filename << std::endl;
        throw std::runtime_error(
            "Not possible to write thermodynamic "
            "lattice output to file.");
      }
      fp = output_binary_files_[ThermodynamicQuantity::LandauVelocity];
    } else {
      try {
        output_binary_files_[ThermodynamicQuantity::j_QBS]->open(filename,
                                                                 std::ios::out);
      } catch (std::ofstream::failure &e) {
        logg[LogArea::Main::id].fatal()
            << "Error in opening " << filename << std::endl;
        throw std::runtime_error(
            "Not possible to write thermodynamic "
            "lattice output to file.");
      }
    }
    write_therm_lattice_binary_header(fp, tq);
  }
}
 
void ThermodynamicLatticeOutput::at_eventend(const ThermodynamicQuantity tq) {
  if (!enable_output_) {
    return;
  }
  if (tq == ThermodynamicQuantity::EckartDensity) {
    if (enable_ascii_) {
      output_ascii_files_[ThermodynamicQuantity::EckartDensity]->close();
    }
    if (enable_binary_) {
      output_binary_files_[ThermodynamicQuantity::EckartDensity]->close();
    }
    return;
  }
  if (tq == ThermodynamicQuantity::Tmn) {
    if (enable_ascii_) {
      output_ascii_files_[ThermodynamicQuantity::Tmn]->close();
    }
    if (enable_binary_) {
      output_binary_files_[ThermodynamicQuantity::Tmn]->close();
    }
    return;
  }
  if (tq == ThermodynamicQuantity::TmnLandau) {
    if (enable_ascii_) {
      output_ascii_files_[ThermodynamicQuantity::TmnLandau]->close();
    }
    if (enable_binary_) {
      output_binary_files_[ThermodynamicQuantity::TmnLandau]->close();
    }
    return;
  }
  if (tq == ThermodynamicQuantity::LandauVelocity) {
    if (enable_ascii_) {
      output_ascii_files_[ThermodynamicQuantity::LandauVelocity]->close();
    }
    if (enable_binary_) {
      output_binary_files_[ThermodynamicQuantity::LandauVelocity]->close();
    }
    return;
  }
  if (tq == ThermodynamicQuantity::j_QBS) {
    if (enable_ascii_) {
      output_ascii_files_[ThermodynamicQuantity::j_QBS]->close();
    }
    if (enable_binary_) {
      output_binary_files_[ThermodynamicQuantity::j_QBS]->close();
    }
    return;
  }
}
 
void ThermodynamicLatticeOutput::thermodynamics_lattice_output(
    RectangularLattice<DensityOnLattice> &lattice, double ctime) {
  double result;
  const auto dim = lattice.n_cells();
  std::shared_ptr<std::ofstream> fp(nullptr);
  if (enable_ascii_) {
    fp = output_ascii_files_[ThermodynamicQuantity::EckartDensity];
    *fp << std::setprecision(14);
    *fp << std::scientific;
    *fp << ctime << std::endl;
  }
  if (enable_binary_) {
    fp = output_binary_files_[ThermodynamicQuantity::EckartDensity];
    assert(sizeof(ctime) == sizeof(double));
    fp->write(reinterpret_cast<char *>(&ctime), sizeof(ctime));
  }
  lattice.iterate_sublattice(
      {0, 0, 0}, dim, [&](DensityOnLattice &node, int ix, int, int) {
        if (enable_ascii_) {
          *fp << node.rho() << " ";
          if (ix == dim[0] - 1) {
            *fp << "\n";
          }
        }
        if (enable_binary_) {
          result = node.rho();
          fp->write(reinterpret_cast<char *>(&result), sizeof(double));
        }
      });
}
 
void ThermodynamicLatticeOutput::thermodynamics_lattice_output(
    RectangularLattice<DensityOnLattice> &lattice, double ctime,
    const std::vector<Particles> &ensembles,
    const DensityParameters &dens_param) {
  if (!enable_output_) {
    return;
  }
  double result;
  const auto dim = lattice.n_cells();
  std::shared_ptr<std::ofstream> fp(nullptr);
  FourVector jQ = FourVector(), jB = FourVector(), jS = FourVector();
  constexpr bool compute_gradient = false;
  if (enable_ascii_) {
    fp = output_ascii_files_[ThermodynamicQuantity::j_QBS];
    *fp << std::setprecision(14);
    *fp << std::scientific;
    *fp << ctime << std::endl;
  }
  if (enable_binary_) {
    fp = output_binary_files_[ThermodynamicQuantity::j_QBS];
    assert(sizeof(ctime) == sizeof(double));
    fp->write(reinterpret_cast<char *>(&ctime), sizeof(ctime));
  }
  lattice.iterate_sublattice(
      {0, 0, 0}, dim, [&](DensityOnLattice &, int ix, int iy, int iz) {
        const ThreeVector position = lattice.cell_center(ix, iy, iz);
        jQ.reset();
        jB.reset();
        jS.reset();
        for (const Particles &particles : ensembles) {
          jQ += std::get<1>(current_eckart(
              position, particles, dens_param, DensityType::Charge,
              compute_gradient, out_par_.td_smearing));
          jB += std::get<1>(current_eckart(
              position, particles, dens_param, DensityType::Baryon,
              compute_gradient, out_par_.td_smearing));
          jS += std::get<1>(current_eckart(
              position, particles, dens_param, DensityType::Strangeness,
              compute_gradient, out_par_.td_smearing));
        }
        if (enable_ascii_) {
          *fp << jQ[0];
          for (int l = 1; l < 4; l++) {
            *fp << " " << jQ[l];
          }
          for (int l = 0; l < 4; l++) {
            *fp << " " << jB[l];
          }
          for (int l = 0; l < 4; l++) {
            *fp << " " << jS[l];
          }
          *fp << "\n";
        }
        if (enable_binary_) {
          for (int l = 0; l < 4; l++) {
            result = jQ[l];
            fp->write(reinterpret_cast<char *>(&result), sizeof(double));
          }
          for (int l = 0; l < 4; l++) {
            result = jB[l];
            fp->write(reinterpret_cast<char *>(&result), sizeof(double));
          }
          for (int l = 0; l < 4; l++) {
            result = jS[l];
            fp->write(reinterpret_cast<char *>(&result), sizeof(double));
          }
        }
      });
}
 
void ThermodynamicLatticeOutput::thermodynamics_lattice_output(
    const ThermodynamicQuantity tq,
    RectangularLattice<EnergyMomentumTensor> &lattice, double ctime) {
  if (!enable_output_) {
    return;
  }
  double result;
  const auto dim = lattice.n_cells();
  std::shared_ptr<std::ofstream> fp(nullptr);
  if (enable_ascii_) {
    switch (tq) {
      case ThermodynamicQuantity::Tmn:
        fp = output_ascii_files_[ThermodynamicQuantity::Tmn];
        break;
      case ThermodynamicQuantity::TmnLandau:
        fp = output_ascii_files_[ThermodynamicQuantity::TmnLandau];
        break;
      case ThermodynamicQuantity::LandauVelocity:
        fp = output_ascii_files_[ThermodynamicQuantity::LandauVelocity];
        break;
      default:
        return;
    }
    *fp << std::setprecision(14);
    *fp << std::scientific;
    *fp << ctime << std::endl;
  }
  if (enable_binary_) {
    switch (tq) {
      case ThermodynamicQuantity::Tmn:
        fp = output_binary_files_[ThermodynamicQuantity::Tmn];
        break;
      case ThermodynamicQuantity::TmnLandau:
        fp = output_binary_files_[ThermodynamicQuantity::TmnLandau];
        break;
      case ThermodynamicQuantity::LandauVelocity:
        fp = output_binary_files_[ThermodynamicQuantity::LandauVelocity];
        break;
      default:
        return;
    }
    assert(sizeof(ctime) == sizeof(double));
    fp->write(reinterpret_cast<char *>(&ctime), sizeof(double));
  }
  switch (tq) {
    case ThermodynamicQuantity::Tmn:
      for (int i = 0; i < 4; i++) {
        for (int j = i; j < 4; j++) {
          lattice.iterate_sublattice(
              {0, 0, 0}, dim,
              [&](EnergyMomentumTensor &node, int ix, int, int) {
                if (enable_ascii_) {
                  *fp << node[EnergyMomentumTensor::tmn_index(i, j)] << " ";
                  if (ix == dim[0] - 1) {
                    *fp << "\n";
                  }
                }
                if (enable_binary_) {
                  result = node[EnergyMomentumTensor::tmn_index(i, j)];
                  fp->write(reinterpret_cast<char *>(&result), sizeof(double));
                }
              });
        }
      }
      break;
    case ThermodynamicQuantity::TmnLandau:
      for (int i = 0; i < 4; i++) {
        for (int j = i; j < 4; j++) {
          lattice.iterate_sublattice(
              {0, 0, 0}, dim,
              [&](EnergyMomentumTensor &node, int ix, int, int) {
                if (enable_ascii_) {
                  const FourVector u = node.landau_frame_4velocity();
                  const EnergyMomentumTensor Tmn_L = node.boosted(u);
                  *fp << Tmn_L[EnergyMomentumTensor::tmn_index(i, j)] << " ";
                  if (ix == dim[0] - 1) {
                    *fp << "\n";
                  }
                }
                if (enable_binary_) {
                  const FourVector u = node.landau_frame_4velocity();
                  const EnergyMomentumTensor Tmn_L = node.boosted(u);
                  result = Tmn_L[EnergyMomentumTensor::tmn_index(i, j)];
                  fp->write(reinterpret_cast<char *>(&result), sizeof(double));
                }
              });
        }
      }
      break;
    case ThermodynamicQuantity::LandauVelocity:
      lattice.iterate_sublattice(
          {0, 0, 0}, dim, [&](EnergyMomentumTensor &node, int, int, int) {
            if (enable_ascii_) {
              const FourVector u = node.landau_frame_4velocity();
              const ThreeVector v = -u.velocity();
              *fp << v.x1() << " " << v.x2() << " " << v.x3() << "\n";
            }
            if (enable_binary_) {
              const FourVector u = node.landau_frame_4velocity();
              ThreeVector v = -u.velocity();
              fp->write(reinterpret_cast<char *>(&v), 3 * sizeof(double));
            }
          });
      break;
    default:
      return;
  }
}
 
int ThermodynamicLatticeOutput::to_int(const ThermodynamicQuantity &tq) {
  switch (tq) {
    case ThermodynamicQuantity::EckartDensity:
      return 0;
    case ThermodynamicQuantity::Tmn:
      return 1;
    case ThermodynamicQuantity::TmnLandau:
      return 2;
    case ThermodynamicQuantity::LandauVelocity:
      return 3;
    case ThermodynamicQuantity::j_QBS:
      return 4;
    default:
      throw std::runtime_error(
          "Error when converting a thermodynamic quantity "
          "to an int, unknown quantity.");
  }
}
 
std::string ThermodynamicLatticeOutput::make_filename(const std::string &descr,
                                                      const int event_number,
                                                      const char type) {
  char suffix[13];
  assert((type == 'a') || (type == 'b'));
  if (type == 'a') {
    snprintf(suffix, sizeof(suffix), "_%07i.dat", event_number);
  } else {
    snprintf(suffix, sizeof(suffix), "_%07i.bin", event_number);
  }
  return base_path_.string() + std::string("/") + descr + std::string(suffix);
}
 
std::string ThermodynamicLatticeOutput::make_varname(
    const ThermodynamicQuantity tq, const DensityType dens_type) {
  return std::string(to_string(dens_type)) + std::string("_") +
         std::string(to_string(tq));
}
 
void ThermodynamicLatticeOutput::write_therm_lattice_ascii_header(
    std::shared_ptr<std::ofstream> fp, const ThermodynamicQuantity &tq) {
  *fp << std::setprecision(2);
  *fp << std::fixed;
  *fp << "#Thermodynamic Lattice Output version: "
      << ThermodynamicLatticeOutput::version << std::endl;
  *fp << std::setprecision(6);
  *fp << "#Quantity:"
      << " " << std::string(to_string(tq)) << std::endl;
  *fp << "#Grid dimensions: " << nodes_[0] << " " << nodes_[1] << " "
      << nodes_[2] << std::endl;
  *fp << "#Grid spacing: " << sizes_[0] << " " << sizes_[1] << " " << sizes_[2]
      << std::endl;
  *fp << "#Grid origin: " << origin_[0] << " " << origin_[1] << " "
      << origin_[2] << std::endl;
}
 
void ThermodynamicLatticeOutput::write_therm_lattice_binary_header(
    std::shared_ptr<std::ofstream> fp, const ThermodynamicQuantity &tq) {
  auto variable_id = to_int(tq);
  fp->write(
      reinterpret_cast<const char *>(&ThermodynamicLatticeOutput::version),
      sizeof(double));
  fp->write(reinterpret_cast<char *>(&variable_id), sizeof(int));
  fp->write(reinterpret_cast<char *>(&nodes_), sizeof(nodes_));
  fp->write(reinterpret_cast<char *>(&sizes_), sizeof(sizes_));
  fp->write(reinterpret_cast<char *>(&origin_), sizeof(origin_));
}
}  // namespace smash