Merge pull request #3204 from nilsdeppe/subcell_grmhd_tci_on_dg_grid

Subcell: Add TciOnDgGrid to ValenciaDivClean

src/Evolution/Systems/GrMhd/ValenciaDivClean/CMakeLists.txt

@@ -46,6 +46,7 @@ target_link_libraries(
   PUBLIC
   Boost::boost
   DataStructures
+  DgSubcell
   ErrorHandling
   GeneralRelativity
   Hydro

src/Evolution/Systems/GrMhd/ValenciaDivClean/Subcell/CMakeLists.txt

@@ -4,6 +4,7 @@
 spectre_target_sources(
   ${LIBRARY}
   PRIVATE
+  TciOnDgGrid.cpp
   TciOptions.cpp
   )
 
@@ -12,5 +13,6 @@ spectre_target_headers(
   INCLUDE_DIRECTORY ${CMAKE_SOURCE_DIR}/src
   HEADERS
   Subcell.hpp
+  TciOnDgGrid.hpp
   TciOptions.hpp
   )

src/Evolution/Systems/GrMhd/ValenciaDivClean/Subcell/TciOnDgGrid.cpp

@@ -0,0 +1,170 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#include "Evolution/Systems/GrMhd/ValenciaDivClean/Subcell/TciOnDgGrid.hpp"
+
+#include <cstddef>
+
+#include "DataStructures/DataVector.hpp"
+#include "DataStructures/Tags/TempTensor.hpp"
+#include "DataStructures/Tensor/EagerMath/DotProduct.hpp"
+#include "DataStructures/Tensor/Tensor.hpp"
+#include "DataStructures/Variables.hpp"
+#include "Evolution/DgSubcell/PerssonTci.hpp"
+#include "Evolution/Systems/GrMhd/ValenciaDivClean/KastaunEtAl.hpp"
+#include "Evolution/Systems/GrMhd/ValenciaDivClean/NewmanHamlin.hpp"
+#include "Evolution/Systems/GrMhd/ValenciaDivClean/PalenzuelaEtAl.hpp"
+#include "Evolution/Systems/GrMhd/ValenciaDivClean/PrimitiveFromConservative.hpp"
+#include "Evolution/Systems/GrMhd/ValenciaDivClean/Subcell/TciOptions.hpp"
+#include "NumericalAlgorithms/Spectral/Mesh.hpp"
+#include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp"
+#include "Utilities/GenerateInstantiations.hpp"
+#include "Utilities/Gsl.hpp"
+
+namespace grmhd::ValenciaDivClean::subcell {
+template <typename RecoveryScheme>
+template <size_t ThermodynamicDim>
+bool TciOnDgGrid<RecoveryScheme>::apply(
+    const gsl::not_null<Variables<hydro::grmhd_tags<DataVector>>*> dg_prim_vars,
+    const Scalar<DataVector>& tilde_d, const Scalar<DataVector>& tilde_tau,
+    const tnsr::i<DataVector, 3, Frame::Inertial>& tilde_s,
+    const tnsr::I<DataVector, 3, Frame::Inertial>& tilde_b,
+    const Scalar<DataVector>& tilde_phi,
+    const tnsr::ii<DataVector, 3, Frame::Inertial>& spatial_metric,
+    const tnsr::II<DataVector, 3, Frame::Inertial>& inv_spatial_metric,
+    const Scalar<DataVector>& sqrt_det_spatial_metric,
+    const EquationsOfState::EquationOfState<true, ThermodynamicDim>& eos,
+    const Mesh<3>& dg_mesh, const TciOptions& tci_options,
+    const double persson_exponent) noexcept {
+  constexpr double persson_tci_epsilon = 1.0e-18;
+  const size_t number_of_points = dg_mesh.number_of_grid_points();
+  Variables<hydro::grmhd_tags<DataVector>> temp_prims(number_of_points);
+
+  // require: tilde_d/sqrt{gamma} >= 0.0 (or some positive user-specified value)
+  if (min(get(tilde_d) / get(sqrt_det_spatial_metric)) <
+      tci_options.minimum_rest_mass_density_times_lorentz_factor) {
+    return true;
+  }
+
+  // Check if we are in atmosphere (but not negative tildeD), and if so, then we
+  // continue using DG on this element.
+  if (max(get(tilde_d) / get(sqrt_det_spatial_metric) /
+          get(get<hydro::Tags::LorentzFactor<DataVector>>(*dg_prim_vars))) <
+          tci_options.atmosphere_density and
+      max(get(get<hydro::Tags::RestMassDensity<DataVector>>(*dg_prim_vars))) <
+          tci_options.atmosphere_density) {
+    // In atmosphere, we only need to recover the primitive variables for the
+    // magnetic field and divergence cleaning field
+    for (size_t i = 0; i < 3; ++i) {
+      get<hydro::Tags::MagneticField<DataVector, 3>>(*dg_prim_vars).get(i) =
+          tilde_b.get(i) / get(sqrt_det_spatial_metric);
+    }
+    get(get<hydro::Tags::DivergenceCleaningField<DataVector>>(*dg_prim_vars)) =
+        get(tilde_phi) / get(sqrt_det_spatial_metric);
+
+    return false;
+  }
+
+  {
+    // require: tilde{B}^2 <= 2sqrt{gamma}(1-epsilon_B)\tilde{tau}
+    Scalar<DataVector>& tilde_b_squared =
+        get<hydro::Tags::RestMassDensity<DataVector>>(temp_prims);
+    dot_product(make_not_null(&tilde_b_squared), tilde_b, tilde_b,
+                spatial_metric);
+    for (size_t i = 0; i < number_of_points; ++i) {
+      if (get(tilde_b_squared)[i] >
+          (1.0 - tci_options.safety_factor_for_magnetic_field) * 2.0 *
+              get(tilde_tau)[i] * get(sqrt_det_spatial_metric)[i]) {
+        return true;
+      }
+    }
+  }
+
+  // Try to recover the primitive variables.
+  // We assign them to a temporary so that if recovery fails at any of the
+  // points we can use the valid primitives at the current time to provide a
+  // high-order initial guess for the recovery on the subcells.
+  //
+  // Copy over the pressure since it's used as an initial guess in some
+  // recovery schemes.
+  get<hydro::Tags::Pressure<DataVector>>(temp_prims) =
+      get<hydro::Tags::Pressure<DataVector>>(*dg_prim_vars);
+  if (not grmhd::ValenciaDivClean::
+          PrimitiveFromConservative<tmpl::list<RecoveryScheme>, false>::apply(
+              make_not_null(
+                  &get<hydro::Tags::RestMassDensity<DataVector>>(temp_prims)),
+              make_not_null(
+                  &get<hydro::Tags::SpecificInternalEnergy<DataVector>>(
+                      temp_prims)),
+              make_not_null(&get<hydro::Tags::SpatialVelocity<DataVector, 3>>(
+                  temp_prims)),
+              make_not_null(
+                  &get<hydro::Tags::MagneticField<DataVector, 3>>(temp_prims)),
+              make_not_null(
+                  &get<hydro::Tags::DivergenceCleaningField<DataVector>>(
+                      temp_prims)),
+              make_not_null(
+                  &get<hydro::Tags::LorentzFactor<DataVector>>(temp_prims)),
+              make_not_null(
+                  &get<hydro::Tags::Pressure<DataVector>>(temp_prims)),
+              make_not_null(
+                  &get<hydro::Tags::SpecificEnthalpy<DataVector>>(temp_prims)),
+              tilde_d, tilde_tau, tilde_s, tilde_b, tilde_phi, spatial_metric,
+              inv_spatial_metric, sqrt_det_spatial_metric, eos)) {
+    return true;
+  }
+
+  // Check if we are in atmosphere after recovery. Unlikely we'd hit this and
+  // not the check before the recovery, but just in case.
+  if (max(get(get<hydro::Tags::RestMassDensity<DataVector>>(temp_prims))) <
+      tci_options.atmosphere_density) {
+    *dg_prim_vars = std::move(temp_prims);
+    return false;
+  }
+
+  // Check that tilde_d and tilde_tau satisfy the Persson TCI
+  if (evolution::dg::subcell::persson_tci(tilde_d, dg_mesh, persson_exponent,
+                                          persson_tci_epsilon) or
+      evolution::dg::subcell::persson_tci(tilde_tau, dg_mesh, persson_exponent,
+                                          persson_tci_epsilon)) {
+    return true;
+  }
+
+  *dg_prim_vars = std::move(temp_prims);
+  return false;
+}
+
+#define RECOVERY(data) BOOST_PP_TUPLE_ELEM(0, data)
+#define INSTANTIATION(r, data) template class TciOnDgGrid<RECOVERY(data)>;
+GENERATE_INSTANTIATIONS(
+    INSTANTIATION,
+    (grmhd::ValenciaDivClean::PrimitiveRecoverySchemes::KastaunEtAl,
+     grmhd::ValenciaDivClean::PrimitiveRecoverySchemes::NewmanHamlin,
+     grmhd::ValenciaDivClean::PrimitiveRecoverySchemes::PalenzuelaEtAl))
+#undef INSTANTIATION
+
+#define THERMO_DIM(data) BOOST_PP_TUPLE_ELEM(1, data)
+#define INSTANTIATION(r, data)                                                \
+  template bool TciOnDgGrid<RECOVERY(data)>::apply<THERMO_DIM(data)>(         \
+      const gsl::not_null<Variables<hydro::grmhd_tags<DataVector>>*>          \
+          dg_prim_vars,                                                       \
+      const Scalar<DataVector>& tilde_d, const Scalar<DataVector>& tilde_tau, \
+      const tnsr::i<DataVector, 3, Frame::Inertial>& tilde_s,                 \
+      const tnsr::I<DataVector, 3, Frame::Inertial>& tilde_b,                 \
+      const Scalar<DataVector>& tilde_phi,                                    \
+      const tnsr::ii<DataVector, 3, Frame::Inertial>& spatial_metric,         \
+      const tnsr::II<DataVector, 3, Frame::Inertial>& inv_spatial_metric,     \
+      const Scalar<DataVector>& sqrt_det_spatial_metric,                      \
+      const EquationsOfState::EquationOfState<true, THERMO_DIM(data)>& eos,   \
+      const Mesh<3>& dg_mesh, const TciOptions& tci_options,                  \
+      const double persson_exponent) noexcept;
+GENERATE_INSTANTIATIONS(
+    INSTANTIATION,
+    (grmhd::ValenciaDivClean::PrimitiveRecoverySchemes::KastaunEtAl,
+     grmhd::ValenciaDivClean::PrimitiveRecoverySchemes::NewmanHamlin,
+     grmhd::ValenciaDivClean::PrimitiveRecoverySchemes::PalenzuelaEtAl),
+    (1, 2))
+#undef INSTANTIATION
+#undef THERMO_DIM
+#undef RECOVERY
+}  // namespace grmhd::ValenciaDivClean::subcell

src/Evolution/Systems/GrMhd/ValenciaDivClean/Subcell/TciOnDgGrid.hpp

@@ -0,0 +1,96 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#pragma once
+
+#include <cstddef>
+
+#include "DataStructures/Tensor/TypeAliases.hpp"
+#include "DataStructures/VariablesTag.hpp"
+#include "Domain/Tags.hpp"
+#include "Evolution/Systems/GrMhd/ValenciaDivClean/Subcell/TciOptions.hpp"
+#include "Evolution/Systems/GrMhd/ValenciaDivClean/Tags.hpp"
+#include "PointwiseFunctions/GeneralRelativity/TagsDeclarations.hpp"
+#include "PointwiseFunctions/Hydro/Tags.hpp"
+#include "Utilities/TMPL.hpp"
+
+/// \cond
+class DataVector;
+namespace EquationsOfState {
+template <bool IsRelativistic, size_t ThermodynamicDim>
+class EquationOfState;
+}  // namespace EquationsOfState
+template <size_t Dim>
+class Mesh;
+namespace gsl {
+template <typename T>
+class not_null;
+}  // namespace gsl
+template <typename TagsList>
+class Variables;
+/// \endcond
+
+namespace grmhd::ValenciaDivClean::subcell {
+/*!
+ * \brief The troubled-cell indicator run on the DG grid to check if the
+ * solution is admissible.
+ *
+ * We denote variables at the candidate solution's time level by a superscript
+ * \f$n+1\f$ and at the time level where the solution is known to be admissible
+ * by a superscript \f$n\f$.
+ *
+ * The following checks are done in the order they are listed:
+ *
+ * - if \f$\min(\tilde{D}^{n+1}/\sqrt{\gamma^{n}})\f$ is less than
+ *   `tci_options.minimum_rest_mass_density_times_lorentz_factor` then we have a
+ *   negative (or extremely small) density and the cell is troubled. Note that
+ *   if this `tci_option` is approximately equal to or larger than the
+ *   `atmosphere_density`, the atmosphere will be flagged as troubled.
+ * - if \f$\max(\tilde{D}^{n+1}/(\sqrt{\gamma^n}W^n))\f$ and \f$\max(\rho^n)\f$
+ *   are less than `tci_options.atmosphere_density` then the entire DG element
+ *   is in atmosphere and it is _not_ troubled.
+ * - if
+ *   \f$(\tilde{B}^{n+1})^2>2\sqrt{\gamma^n}(1-\epsilon_B)\tilde{\tau}^{n+1}\f$
+ *   at any grid point, then the cell is troubled
+ * - attempt a primitive recovery using the `RecoveryScheme` from the template
+ *   parameter. The cell is marked as troubled if the primitive recovery fails
+ *   at any grid point.
+ * - if \f$\max(\rho^{n+1})\f$ is below `tci_options.atmosphere_density` then
+ *   the cell is in atmosphere and not marked as troubled. Note that the
+ *   magnetic field is still freely evolved.
+ * - apply the Persson TCI to \f$\tilde{D}^{n+1}\f$ and \f$\tilde{\tau}^{n+1}\f$
+ *
+ * If the cell is not flagged as troubled then the primitives are computed at
+ * time level `n+1`.
+ */
+template <typename RecoveryScheme>
+class TciOnDgGrid {
+ public:
+  using return_tags =
+      tmpl::list<::Tags::Variables<hydro::grmhd_tags<DataVector>>>;
+  using argument_tags =
+      tmpl::list<grmhd::ValenciaDivClean::Tags::TildeD,
+                 grmhd::ValenciaDivClean::Tags::TildeTau,
+                 grmhd::ValenciaDivClean::Tags::TildeS<>,
+                 grmhd::ValenciaDivClean::Tags::TildeB<>,
+                 grmhd::ValenciaDivClean::Tags::TildePhi,
+                 gr::Tags::SpatialMetric<3>, gr::Tags::InverseSpatialMetric<3>,
+                 gr::Tags::SqrtDetSpatialMetric<>,
+                 hydro::Tags::EquationOfStateBase, domain::Tags::Mesh<3>,
+                 Tags::TciOptions>;
+
+  template <size_t ThermodynamicDim>
+  static bool apply(
+      gsl::not_null<Variables<hydro::grmhd_tags<DataVector>>*> dg_prim_vars,
+      const Scalar<DataVector>& tilde_d, const Scalar<DataVector>& tilde_tau,
+      const tnsr::i<DataVector, 3, Frame::Inertial>& tilde_s,
+      const tnsr::I<DataVector, 3, Frame::Inertial>& tilde_b,
+      const Scalar<DataVector>& tilde_phi,
+      const tnsr::ii<DataVector, 3, Frame::Inertial>& spatial_metric,
+      const tnsr::II<DataVector, 3, Frame::Inertial>& inv_spatial_metric,
+      const Scalar<DataVector>& sqrt_det_spatial_metric,
+      const EquationsOfState::EquationOfState<true, ThermodynamicDim>& eos,
+      const Mesh<3>& dg_mesh, const TciOptions& tci_options,
+      double persson_exponent) noexcept;
+};
+}  // namespace grmhd::ValenciaDivClean::subcell

tests/Unit/Evolution/Systems/GrMhd/ValenciaDivClean/CMakeLists.txt

@@ -7,6 +7,7 @@ set(LIBRARY_SOURCES
   BoundaryConditions/Test_DirichletAnalytic.cpp
   BoundaryConditions/Test_Periodic.cpp
   BoundaryCorrections/Test_Rusanov.cpp
+  Subcell/Test_TciOnDgGrid.cpp
   Subcell/Test_TciOptions.cpp
   Test_Characteristics.cpp
   Test_ConservativeFromPrimitive.cpp