Add executable to export EOS for RotNS

docs/GroupDefs.hpp

@@ -442,6 +442,11 @@
  * <td> \ref ParallelInfoExecutablePage "ParallelInfo" </td>
  * <td> Executable for checking number of nodes, cores, etc.</td>
  * </tr>
+ * <tr>
+ * <td> ExportEquationOfStateForRotNS </td>
+ * <td> Exports a 1d equation of state in a table format that the RotNS fortran
+ * code can read in.</td>
+ * </tr>
  * </table>
  */
 

src/Executables/CMakeLists.txt

@@ -5,6 +5,7 @@ add_subdirectory(Benchmark)
 add_subdirectory(CombineH5)
 add_subdirectory(ConvertComposeTable)
 add_subdirectory(DebugPreprocessor)
+add_subdirectory(ExportEquationOfStateForRotNS)
 add_subdirectory(Examples)
 add_subdirectory(ExportCoordinates)
 add_subdirectory(FindHorizons)

src/Executables/ExportEquationOfStateForRotNS/CMakeLists.txt

@@ -0,0 +1,21 @@
+# Distributed under the MIT License.
+# See LICENSE.txt for details.
+
+set(EXECUTABLE ExportEquationOfStateForRotNS)
+
+add_spectre_executable(
+  ${EXECUTABLE}
+  EXCLUDE_FROM_ALL
+  ExportEquationOfStateForRotNS.cpp
+  )
+
+target_link_libraries(
+  ${EXECUTABLE}
+  PRIVATE
+  Boost::boost
+  Boost::program_options
+  DataStructures
+  Hydro
+  Parallel
+  Utilities
+  )

src/Executables/ExportEquationOfStateForRotNS/ExportEquationOfStateForRotNS.cpp

@@ -0,0 +1,190 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#include <boost/program_options.hpp>
+#include <cmath>
+#include <cstddef>
+#include <fstream>
+#include <iomanip>
+#include <iterator>
+#include <string>
+#include <utility>
+#include <vector>
+
+#include "DataStructures/DataVector.hpp"
+#include "Options/Options.hpp"
+#include "Options/ParseOptions.hpp"
+#include "Parallel/Printf.hpp"
+#include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp"
+#include "PointwiseFunctions/Hydro/Tags.hpp"
+#include "PointwiseFunctions/Hydro/Units.hpp"
+#include "Utilities/ErrorHandling/Error.hpp"
+#include "Utilities/FileSystem.hpp"
+#include "Utilities/TMPL.hpp"
+
+// Charm looks for this function but since we build without a main function or
+// main module we just have it be empty
+extern "C" void CkRegisterMainModule(void) {}
+
+namespace {
+void dump_barotropic_eos(
+    const EquationsOfState::EquationOfState<true, 1>& eos,
+    const size_t number_of_log10_number_density_points_for_dump,
+    const std::string& output_file_name,
+    const double lower_bound_rest_mass_density_cgs,
+    const double upper_bound_rest_mass_density_cgs) {
+  using std::log10;
+  using std::pow;
+  // Baryon mass, used to go from number density to rest mass
+  // density. I.e. `rho_cgs = n_cgs * baryon_mass`, where `n_gcs` is the number
+  // density in CGS units. This is the baryon mass that RotNS uses. This
+  // might be different from the baryon mass that the EoS uses.
+  const double baryon_mass_of_rotns_cgs = 1.659e-24;
+  const double log10_lower_bound_number_density_cgs =
+      log10(lower_bound_rest_mass_density_cgs / baryon_mass_of_rotns_cgs);
+  const double log10_upper_bound_number_density_cgs =
+      log10(upper_bound_rest_mass_density_cgs / baryon_mass_of_rotns_cgs);
+  const double delta_log_number_density_cgs =
+      (log10_upper_bound_number_density_cgs -
+       log10_lower_bound_number_density_cgs) /
+      static_cast<double>(number_of_log10_number_density_points_for_dump - 1);
+
+  if (file_system::check_if_file_exists(output_file_name)) {
+    ERROR_NO_TRACE("File " << output_file_name
+                           << " already exists. Refusing to overwrite.");
+  }
+  std::ofstream outfile(output_file_name.c_str());
+
+  for (size_t log10_number_density_index = 0;
+       log10_number_density_index <
+       number_of_log10_number_density_points_for_dump;
+       ++log10_number_density_index) {
+    using std::pow;
+    const double number_density_cgs =
+        pow(10.0, log10_lower_bound_number_density_cgs +
+                      static_cast<double>(log10_number_density_index) *
+                          delta_log_number_density_cgs);
+    const double rest_mass_density_cgs =
+        number_density_cgs * baryon_mass_of_rotns_cgs;
+
+    // Note: we will want to add the baryon mass to our EOS interface.
+    const double baryon_mass_of_eos_cgs = baryon_mass_of_rotns_cgs;
+    const Scalar<double> rest_mass_density_geometric{rest_mass_density_cgs /
+                                                     baryon_mass_of_eos_cgs};
+    const Scalar<double> pressure_geometric =
+        eos.pressure_from_density(rest_mass_density_geometric);
+    const Scalar<double> specific_internal_energy_geometric =
+        eos.specific_internal_energy_from_density(rest_mass_density_geometric);
+    const Scalar<double> total_energy_density_geometric{
+        get(rest_mass_density_geometric) *
+        (1.0 + get(specific_internal_energy_geometric))};
+
+    // Note: the energy density is divided by c^2, so the rest-mass part is rho
+    // c^2
+    const double total_energy_density_cgs =
+        get(total_energy_density_geometric) *
+        hydro::units::cgs::rest_mass_density_unit;
+
+    // should be dyne cm^(-3)
+    const double pressure_cgs =
+        get(pressure_geometric) * hydro::units::cgs::pressure_unit;
+
+    outfile << std::scientific << std::setw(24) << std::setprecision(14)
+            << log10(number_density_cgs) << std::setw(24)
+            << std::setprecision(14) << log10(total_energy_density_cgs)
+            << std::setw(24) << std::setprecision(14) << log10(pressure_cgs)
+            << std::endl;
+  }
+  outfile.close();
+}
+
+namespace OptionTags {
+struct NumberOfPoints {
+  using type = size_t;
+  static constexpr Options::String help = {
+      "Number of points at which to dump the EoS"};
+};
+
+struct OutputFileName {
+  using type = std::string;
+  static constexpr Options::String help = {
+      "Name of the output file to dump the EoS to, including file extension."};
+};
+
+struct LowerBoundRestMassDensityCgs {
+  using type = double;
+  static constexpr Options::String help = {
+      "Lower bound of rest mass density in CGS units."};
+};
+
+struct UpperBoundRestMassDensityCgs {
+  using type = double;
+  static constexpr Options::String help = {
+      "Upper bound of rest mass density in CGS units."};
+};
+}  // namespace OptionTags
+}  // namespace
+
+int main(int argc, char** argv) {
+  namespace bpo = boost::program_options;
+  bpo::positional_options_description pos_desc;
+
+  const std::string help_string =
+      "Dump a relativistic barotropic equation of state to disk.\n"
+      "All options controlling input and output are read from the input file.\n"
+      "This executable can be extended to support 2d and 3d EoS, where "
+      "temperature and electron fraction dependence is available.";
+
+  bpo::options_description desc(help_string);
+  desc.add_options()("help,h,", "show this help message")(
+      "input-file", bpo::value<std::string>()->required(), "Input file name")(
+      "check-options", "Check input file options");
+
+  bpo::variables_map vars;
+
+  bpo::store(bpo::command_line_parser(argc, argv)
+                 .positional(pos_desc)
+                 .options(desc)
+                 .run(),
+             vars);
+
+  if (vars.count("help") != 0u or vars.count("input-file") == 0u) {
+    Parallel::printf("%s\n", desc);
+    return 1;
+  }
+
+  using option_list =
+      tmpl::list<hydro::OptionTags::EquationOfState<true, 1>,
+                 OptionTags::NumberOfPoints, OptionTags::OutputFileName,
+                 OptionTags::LowerBoundRestMassDensityCgs,
+                 OptionTags::UpperBoundRestMassDensityCgs>;
+
+  Options::Parser<option_list> option_parser(help_string);
+  option_parser.parse_file(vars["input-file"].as<std::string>());
+
+  if (vars.count("check-options") != 0) {
+    // Force all the options to be created.
+    option_parser.template apply<option_list>([](auto... args) {
+      (void)std::initializer_list<char>{((void)args, '0')...};
+    });
+    Parallel::printf("\n%s parsed successfully!\n",
+                     vars["input-file"].as<std::string>());
+
+    return 0;
+  }
+
+  const auto options =
+      option_parser.template apply<option_list>([](auto... args) {
+        return tuples::tagged_tuple_from_typelist<option_list>(
+            std::move(args)...);
+      });
+
+  dump_barotropic_eos(
+      *get<hydro::OptionTags::EquationOfState<true, 1>>(options),
+      get<OptionTags::NumberOfPoints>(options),
+      get<OptionTags::OutputFileName>(options),
+      get<OptionTags::LowerBoundRestMassDensityCgs>(options),
+      get<OptionTags::UpperBoundRestMassDensityCgs>(options));
+
+  return 0;
+}

tests/InputFiles/ExportEquationOfStateForRotNS/ExportEosForRotNS.yaml

@@ -0,0 +1,56 @@
+# Distributed under the MIT License.
+# See LICENSE.txt for details.
+
+# Executable: ExportEquationOfStateForRotNS
+# Check: parse;execute
+# ExpectedOutput:
+#   EOS.GN
+
+EquationOfState:
+  Enthalpy(Enthalpy(Enthalpy(Spectral))):
+    ReferenceDensity: 0.001822038675716006
+    MinimumDensity: 0.001822038675716006
+    MaximumDensity: 0.0034
+    TrigScaling: 0.0
+    PolynomialCoefficients:
+      [1.2730598378768818, 1.2730598378768818, 0.6365299189384409,
+      0.21217663964614697, 0.05304415991153674, 0.01060883198230735,
+      0.0017681386637178915, 0.0002525912376739845, 3.157390470924806e-05,
+      3.5082116343608955e-06]
+    SinCoefficients: [0.0]
+    CosCoefficients: [0.0]
+    TransitionDeltaEpsilon: 0.0
+    Enthalpy:
+      ReferenceDensity: 0.0011334511674839399
+      MinimumDensity: 0.0011334511674839399
+      MaximumDensity: 0.001822038675716006
+      TrigScaling: 0.0
+      PolynomialCoefficients: [1.272585148960061, 0.001]
+      SinCoefficients: [0.0]
+      CosCoefficients: [0.0]
+      TransitionDeltaEpsilon: 0.0
+      Enthalpy:
+        ReferenceDensity: 0.00022669023349678794
+        MinimumDensity: 0.0004533804669935759
+        MaximumDensity: 0.0011334511674839399
+        TrigScaling: 1.26426988871305
+        PolynomialCoefficients:
+          [1.0,0.08063293075870805,4.406887319408924e-26,
+          8.177895241388924e-22,0.013558242085066733,0.004117320982626606,
+          9.757362504479485e-26,1.5646573325075753e-30,0.00016253964205058317]
+        SinCoefficients:
+          [0.0003763514388305583,0.017968749910748837,0.008140052979970034,
+          -0.003067418379116628,-0.0008236601907322793]
+        CosCoefficients:
+          [-0.01080996024705052,-0.003421193490191067,0.012325774692378716,
+          0.004367136076912163,-0.00020374276952538073]
+        TransitionDeltaEpsilon: 0.0
+        Spectral:
+          ReferenceDensity: 4.533804669935759e-05
+          ReferencePressure: 9.970647727158039e-08
+          Coefficients: [1.2, 0.0, 1.34440187653529, -0.46098357752567365]
+          UpperDensity: 0.0004533804669935759
+OutputFileName: ./EOS.GN
+NumberOfPoints: 40000
+LowerBoundRestMassDensityCgs: 1.0e3
+UpperBoundRestMassDensityCgs: 1.0e16