Merge pull request #6134 from ctrandrsn/3dEosSpecID

Add support for 3d EOSes for Spec initial data
sxs-collaboration · Jun 29, 2024 · 09e5715 · 09e5715
2 parents e9e6cc7 + 392867c
commit 09e5715
Show file tree

Hide file tree

Showing 2 changed files with 22 additions and 3 deletions.
diff --git a/src/Evolution/Systems/GrMhd/GhValenciaDivClean/AllSolutions.hpp b/src/Evolution/Systems/GrMhd/GhValenciaDivClean/AllSolutions.hpp
@@ -10,7 +10,8 @@
 #ifdef HAS_SPEC_EXPORTER
 #include "PointwiseFunctions/AnalyticData/GrMhd/SpecInitialData.hpp"
 using SpecInitialDataList = tmpl::list<grmhd::AnalyticData::SpecInitialData<1>,
-                                       grmhd::AnalyticData::SpecInitialData<2>>;
+                                       grmhd::AnalyticData::SpecInitialData<2>,
+                                       grmhd::AnalyticData::SpecInitialData<3>>;
 #else
 using SpecInitialDataList = NoSuchType;
 #endif

diff --git a/src/PointwiseFunctions/AnalyticData/GrMhd/SpecInitialData.cpp b/src/PointwiseFunctions/AnalyticData/GrMhd/SpecInitialData.cpp
@@ -106,11 +106,15 @@ void SpecInitialData<ThermodynamicDim>::VariablesComputer<DataType>::operator()(
       get<hydro::Tags::RestMassDensity<DataType>>(interpolated_data);
   const auto& temperature =
       cache->get_var(*this, hydro::Tags::Temperature<DataType>{});
+  const auto& electron_fraction =
+      cache->get_var(*this, hydro::Tags::ElectronFraction<DataType>{});
   const size_t num_points = get_size(get(rest_mass_density));
   for (size_t i = 0; i < num_points; ++i) {
     const double local_rest_mass_density =
         get_element(get(rest_mass_density), i);
     const double local_temperature = get_element(get(temperature), i);
+    const double local_electron_fraction =
+        get_element(get(electron_fraction), i);
     if constexpr (ThermodynamicDim == 1) {
       get_element(get(*specific_internal_energy), i) =
           get(eos.specific_internal_energy_from_density(
@@ -121,7 +125,12 @@ void SpecInitialData<ThermodynamicDim>::VariablesComputer<DataType>::operator()(
               Scalar<double>(local_rest_mass_density),
               Scalar<double>(local_temperature)));
     } else {
-      ERROR("Only 1d & 2d EOS is currently supported for SpEC ID");
+      static_assert(ThermodynamicDim == 3);
+      get_element(get(*specific_internal_energy), i) =
+          get(eos.specific_internal_energy_from_density_and_temperature(
+              Scalar<double>(local_rest_mass_density),
+              Scalar<double>(local_temperature),
+              Scalar<double>(local_electron_fraction)));
     }
   }
 }
@@ -136,11 +145,15 @@ void SpecInitialData<ThermodynamicDim>::VariablesComputer<DataType>::operator()(
       get<hydro::Tags::RestMassDensity<DataType>>(interpolated_data);
   const auto& temperature =
       cache->get_var(*this, hydro::Tags::Temperature<DataType>{});
+  const auto& electron_fraction =
+      cache->get_var(*this, hydro::Tags::ElectronFraction<DataType>{});
   const size_t num_points = get_size(get(rest_mass_density));
   for (size_t i = 0; i < num_points; ++i) {
     const double local_rest_mass_density =
         get_element(get(rest_mass_density), i);
     const double local_temperature = get_element(get(temperature), i);
+    const double local_electron_fraction =
+        get_element(get(electron_fraction), i);
     if constexpr (ThermodynamicDim == 1) {
       get_element(get(*pressure), i) = get(
           eos.pressure_from_density(Scalar<double>(local_rest_mass_density)));
@@ -152,7 +165,12 @@ void SpecInitialData<ThermodynamicDim>::VariablesComputer<DataType>::operator()(
                   Scalar<double>(local_rest_mass_density),
                   Scalar<double>(local_temperature))))));
     } else {
-      ERROR("Only 1d & 2d EOS is currently supported for SpEC ID");
+      static_assert(ThermodynamicDim == 3);
+      get_element(get(*pressure), i) =
+          get(eos.pressure_from_density_and_temperature(
+              Scalar<double>(local_rest_mass_density),
+              Scalar<double>(local_temperature),
+              Scalar<double>(local_electron_fraction)));
     }
   }
 }