diff --git a/docs/References.bib b/docs/References.bib
index 644cc10e6fd6..02b048d421ff 100644
--- a/docs/References.bib
+++ b/docs/References.bib
@@ -2198,6 +2198,21 @@ @article{Shibata2003
url = "https://link.aps.org/doi/10.1103/PhysRevD.68.104020"
}
+@article{Shiokawa2011ih,
+ author = "Shiokawa, Hotaka and Dolence, Joshua C. and Gammie, Charles F.
+ and Noble, Scott C.",
+ title = "{Global GRMHD Simulations of Black Hole Accretion Flows:
+ a Convergence Study}",
+ eprint = "1111.0396",
+ archivePrefix = "arXiv",
+ primaryClass = "astro-ph.HE",
+ doi = "10.1088/0004-637X/744/2/187",
+ journal = "Astrophys. J.",
+ volume = "744",
+ pages = "187",
+ year = "2012"
+}
+
@article{Shu1988439,
title = "Efficient implementation of essentially non-oscillatory
shock-capturing schemes",
diff --git a/src/Evolution/Systems/GrMhd/ValenciaDivClean/BoundaryConditions/CMakeLists.txt b/src/Evolution/Systems/GrMhd/ValenciaDivClean/BoundaryConditions/CMakeLists.txt
index 3ae01b6ef84d..1f5e5506c620 100644
--- a/src/Evolution/Systems/GrMhd/ValenciaDivClean/BoundaryConditions/CMakeLists.txt
+++ b/src/Evolution/Systems/GrMhd/ValenciaDivClean/BoundaryConditions/CMakeLists.txt
@@ -8,6 +8,7 @@ spectre_target_sources(
DemandOutgoingCharSpeeds.cpp
DirichletAnalytic.cpp
HydroFreeOutflow.cpp
+ Reflective.cpp
)
spectre_target_headers(
@@ -19,4 +20,5 @@ spectre_target_headers(
DirichletAnalytic.hpp
Factory.hpp
HydroFreeOutflow.hpp
+ Reflective.hpp
)
diff --git a/src/Evolution/Systems/GrMhd/ValenciaDivClean/BoundaryConditions/Factory.hpp b/src/Evolution/Systems/GrMhd/ValenciaDivClean/BoundaryConditions/Factory.hpp
index 5b7db2701acc..48561e5f7ad7 100644
--- a/src/Evolution/Systems/GrMhd/ValenciaDivClean/BoundaryConditions/Factory.hpp
+++ b/src/Evolution/Systems/GrMhd/ValenciaDivClean/BoundaryConditions/Factory.hpp
@@ -8,10 +8,12 @@
#include "Evolution/Systems/GrMhd/ValenciaDivClean/BoundaryConditions/DemandOutgoingCharSpeeds.hpp"
#include "Evolution/Systems/GrMhd/ValenciaDivClean/BoundaryConditions/DirichletAnalytic.hpp"
#include "Evolution/Systems/GrMhd/ValenciaDivClean/BoundaryConditions/HydroFreeOutflow.hpp"
+#include "Evolution/Systems/GrMhd/ValenciaDivClean/BoundaryConditions/Reflective.hpp"
namespace grmhd::ValenciaDivClean::BoundaryConditions {
/// Typelist of standard BoundaryConditions
using standard_boundary_conditions =
tmpl::list<DemandOutgoingCharSpeeds, DirichletAnalytic, HydroFreeOutflow,
+ Reflective,
domain::BoundaryConditions::Periodic<BoundaryCondition>>;
} // namespace grmhd::ValenciaDivClean::BoundaryConditions
diff --git a/src/Evolution/Systems/GrMhd/ValenciaDivClean/BoundaryConditions/Reflective.cpp b/src/Evolution/Systems/GrMhd/ValenciaDivClean/BoundaryConditions/Reflective.cpp
new file mode 100644
index 000000000000..c9cdfc546a6d
--- /dev/null
+++ b/src/Evolution/Systems/GrMhd/ValenciaDivClean/BoundaryConditions/Reflective.cpp
@@ -0,0 +1,553 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#include "Evolution/Systems/GrMhd/ValenciaDivClean/BoundaryConditions/Reflective.hpp"
+
+#include <cstddef>
+#include <memory>
+#include <optional>
+#include <pup.h>
+
+#include "DataStructures/DataVector.hpp"
+#include "DataStructures/Index.hpp"
+#include "DataStructures/SliceVariables.hpp"
+#include "DataStructures/Tags/TempTensor.hpp"
+#include "DataStructures/Tensor/EagerMath/DeterminantAndInverse.hpp"
+#include "DataStructures/Tensor/EagerMath/DotProduct.hpp"
+#include "DataStructures/Tensor/EagerMath/RaiseOrLowerIndex.hpp"
+#include "DataStructures/Tensor/Expressions/Evaluate.hpp"
+#include "DataStructures/Tensor/Expressions/TensorExpression.hpp"
+#include "DataStructures/Tensor/Tensor.hpp"
+#include "DataStructures/Variables.hpp"
+#include "Domain/Structure/Direction.hpp"
+#include "Evolution/DgSubcell/SliceTensor.hpp"
+#include "Evolution/Systems/GrMhd/ValenciaDivClean/ConservativeFromPrimitive.hpp"
+#include "Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/Reconstructor.hpp"
+#include "Evolution/Systems/GrMhd/ValenciaDivClean/Fluxes.hpp"
+#include "NumericalAlgorithms/Spectral/Mesh.hpp"
+#include "PointwiseFunctions/Hydro/LorentzFactor.hpp"
+#include "PointwiseFunctions/Hydro/Tags.hpp"
+#include "Utilities/Gsl.hpp"
+
+namespace grmhd::ValenciaDivClean::BoundaryConditions {
+
+Reflective::Reflective(bool reflect_both) : reflect_both_(reflect_both) {}
+
+Reflective::Reflective(CkMigrateMessage* const msg) : BoundaryCondition(msg) {}
+
+std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>
+Reflective::get_clone() const {
+ return std::make_unique<Reflective>(*this);
+}
+
+void Reflective::pup(PUP::er& p) {
+ BoundaryCondition::pup(p);
+ p | reflect_both_;
+}
+
+// NOLINTNEXTLINE
+PUP::able::PUP_ID Reflective::my_PUP_ID = 0;
+
+std::optional<std::string> Reflective::dg_ghost(
+ const gsl::not_null<Scalar<DataVector>*> tilde_d,
+ const gsl::not_null<Scalar<DataVector>*> tilde_ye,
+ const gsl::not_null<Scalar<DataVector>*> tilde_tau,
+ const gsl::not_null<tnsr::i<DataVector, 3, Frame::Inertial>*> tilde_s,
+ const gsl::not_null<tnsr::I<DataVector, 3, Frame::Inertial>*> tilde_b,
+ const gsl::not_null<Scalar<DataVector>*> tilde_phi,
+
+ const gsl::not_null<tnsr::I<DataVector, 3, Frame::Inertial>*> tilde_d_flux,
+ const gsl::not_null<tnsr::I<DataVector, 3, Frame::Inertial>*> tilde_ye_flux,
+ const gsl::not_null<tnsr::I<DataVector, 3, Frame::Inertial>*>
+ tilde_tau_flux,
+ const gsl::not_null<tnsr::Ij<DataVector, 3, Frame::Inertial>*> tilde_s_flux,
+ const gsl::not_null<tnsr::IJ<DataVector, 3, Frame::Inertial>*> tilde_b_flux,
+ const gsl::not_null<tnsr::I<DataVector, 3, Frame::Inertial>*>
+ tilde_phi_flux,
+
+ const gsl::not_null<Scalar<DataVector>*> lapse,
+ const gsl::not_null<tnsr::I<DataVector, 3, Frame::Inertial>*> shift,
+ const gsl::not_null<tnsr::II<DataVector, 3, Frame::Inertial>*>
+ inv_spatial_metric,
+
+ const std::optional<tnsr::I<DataVector, 3, Frame::Inertial>>&
+ /*face_mesh_velocity*/,
+ const tnsr::i<DataVector, 3, Frame::Inertial>&
+ outward_directed_normal_covector,
+ const tnsr::I<DataVector, 3, Frame::Inertial>&
+ outward_directed_normal_vector,
+
+ const Scalar<DataVector>& interior_rest_mass_density,
+ const Scalar<DataVector>& interior_electron_fraction,
+ const Scalar<DataVector>& interior_specific_internal_energy,
+ const tnsr::I<DataVector, 3, Frame::Inertial>& interior_spatial_velocity,
+ const tnsr::I<DataVector, 3, Frame::Inertial>& interior_magnetic_field,
+ const Scalar<DataVector>& interior_lorentz_factor,
+ const Scalar<DataVector>& interior_pressure,
+
+ const tnsr::I<DataVector, 3, Frame::Inertial>& interior_shift,
+ const Scalar<DataVector>& interior_lapse,
+ const tnsr::II<DataVector, 3, Frame::Inertial>& interior_inv_spatial_metric)
+ const {
+ const size_t number_of_grid_points = get(interior_rest_mass_density).size();
+
+ // temp buffer to store
+ // * n_i v^i
+ // * v_{ghost}^i
+ // * n_i B^i
+ // * B_{ghost}^i
+ // * divergence cleaning field on ghost zone
+ // * spatial metric
+ // * sqrt determinant of spatial metric
+ Variables<tmpl::list<
+ ::Tags::TempScalar<0>, hydro::Tags::SpatialVelocity<DataVector, 3>,
+ ::Tags::TempScalar<1>, hydro::Tags::MagneticField<DataVector, 3>,
+ ::Tags::TempScalar<2>, gr::Tags::SpatialMetric<DataVector, 3>,
+ gr::Tags::SqrtDetSpatialMetric<DataVector>>>
+ temp_buffer{number_of_grid_points};
+ auto& normal_dot_interior_spatial_velocity =
+ get<::Tags::TempScalar<0>>(temp_buffer);
+ auto& exterior_spatial_velocity =
+ get<hydro::Tags::SpatialVelocity<DataVector, 3>>(temp_buffer);
+ auto& normal_dot_interior_magnetic_field =
+ get<::Tags::TempScalar<1>>(temp_buffer);
+ auto& exterior_magnetic_field =
+ get<hydro::Tags::MagneticField<DataVector, 3>>(temp_buffer);
+ auto& exterior_divergence_cleaning_field =
+ get<::Tags::TempScalar<2>>(temp_buffer);
+ auto& interior_spatial_metric =
+ get<gr::Tags::SpatialMetric<DataVector, 3>>(temp_buffer);
+ auto& interior_sqrt_det_spatial_metric =
+ get<gr::Tags::SqrtDetSpatialMetric<DataVector>>(temp_buffer);
+
+ get(*lapse) = get(interior_lapse);
+ for (size_t i = 0; i < 3; ++i) {
+ (*shift).get(i) = interior_shift.get(i);
+ for (size_t j = 0; j < 3; ++j) {
+ (*inv_spatial_metric).get(i, j) = interior_inv_spatial_metric.get(i, j);
+ }
+ }
+
+ // spatial metric and sqrt determinant of spatial metric can be retrived from
+ // Databox but only as gridless_tags with whole volume data (unlike all the
+ // other arguments which are face tensors). Rather than doing expensive tensor
+ // slicing operations on those, we just compute those two quantities from
+ // inverse spatial metric as below.
+ determinant_and_inverse(make_not_null(&interior_sqrt_det_spatial_metric),
+ make_not_null(&interior_spatial_metric),
+ interior_inv_spatial_metric);
+ get(interior_sqrt_det_spatial_metric) =
+ 1.0 / sqrt(get(interior_sqrt_det_spatial_metric));
+
+ // copy-paste interior spatial velocity to exterior spatial velocity, but
+ // kill ingoing normal component to zero
+ dot_product(make_not_null(&normal_dot_interior_spatial_velocity),
+ outward_directed_normal_covector, interior_spatial_velocity);
+ dot_product(make_not_null(&normal_dot_interior_magnetic_field),
+ outward_directed_normal_covector, interior_magnetic_field);
+
+ // reflect both outward and inward normal components of
+ // spatial velocity and magnetic field.
+ if (reflect_both_) {
+ for (size_t i = 0; i < number_of_grid_points; ++i) {
+ if (get(normal_dot_interior_spatial_velocity)[i] > 0.0) {
+ for (size_t spatial_index = 0; spatial_index < 3; ++spatial_index) {
+ exterior_spatial_velocity.get(spatial_index)[i] =
+ interior_spatial_velocity.get(spatial_index)[i] -
+ 2.0 * get(normal_dot_interior_spatial_velocity)[i] *
+ outward_directed_normal_vector.get(spatial_index)[i];
+ }
+ } else {
+ for (size_t spatial_index = 0; spatial_index < 3; ++spatial_index) {
+ exterior_spatial_velocity.get(spatial_index)[i] =
+ interior_spatial_velocity.get(spatial_index)[i] +
+ 2.0 * get(normal_dot_interior_spatial_velocity)[i] *
+ outward_directed_normal_vector.get(spatial_index)[i];
+ }
+ }
+ if (get(normal_dot_interior_magnetic_field)[i] > 0.0) {
+ for (size_t spatial_index = 0; spatial_index < 3; ++spatial_index) {
+ exterior_magnetic_field.get(spatial_index)[i] =
+ interior_magnetic_field.get(spatial_index)[i] -
+ 2.0 * get(normal_dot_interior_magnetic_field)[i] *
+ outward_directed_normal_vector.get(spatial_index)[i];
+ }
+ } else {
+ for (size_t spatial_index = 0; spatial_index < 3; ++spatial_index) {
+ exterior_magnetic_field.get(spatial_index)[i] =
+ interior_magnetic_field.get(spatial_index)[i] +
+ 2.0 * get(normal_dot_interior_magnetic_field)[i] *
+ outward_directed_normal_vector.get(spatial_index)[i];
+ }
+ }
+ }
+ }
+
+ // reflect only the outward normal component of
+ // spatial velocity and magnetic field.
+ else {
+ for (size_t i = 0; i < number_of_grid_points; ++i) {
+ if (get(normal_dot_interior_spatial_velocity)[i] > 0.0) {
+ for (size_t spatial_index = 0; spatial_index < 3; ++spatial_index) {
+ exterior_spatial_velocity.get(spatial_index)[i] =
+ interior_spatial_velocity.get(spatial_index)[i] -
+ 2.0 * get(normal_dot_interior_spatial_velocity)[i] *
+ outward_directed_normal_vector.get(spatial_index)[i];
+ }
+ } else {
+ for (size_t spatial_index = 0; spatial_index < 3; ++spatial_index) {
+ exterior_spatial_velocity.get(spatial_index)[i] =
+ interior_spatial_velocity.get(spatial_index)[i];
+ }
+ }
+ if (get(normal_dot_interior_magnetic_field)[i] > 0.0) {
+ for (size_t spatial_index = 0; spatial_index < 3; ++spatial_index) {
+ exterior_magnetic_field.get(spatial_index)[i] =
+ interior_magnetic_field.get(spatial_index)[i] -
+ 2.0 * get(normal_dot_interior_magnetic_field)[i] *
+ outward_directed_normal_vector.get(spatial_index)[i];
+ }
+ } else {
+ for (size_t spatial_index = 0; spatial_index < 3; ++spatial_index) {
+ exterior_magnetic_field.get(spatial_index)[i] =
+ interior_magnetic_field.get(spatial_index)[i];
+ }
+ }
+ }
+ }
+ get(exterior_divergence_cleaning_field) = 0.0;
+
+ ConservativeFromPrimitive::apply(
+ tilde_d, tilde_ye, tilde_tau, tilde_s, tilde_b, tilde_phi,
+ interior_rest_mass_density, interior_electron_fraction,
+ interior_specific_internal_energy, interior_pressure,
+ exterior_spatial_velocity, interior_lorentz_factor,
+ exterior_magnetic_field, interior_sqrt_det_spatial_metric,
+ interior_spatial_metric, exterior_divergence_cleaning_field);
+
+ ComputeFluxes::apply(tilde_d_flux, tilde_ye_flux, tilde_tau_flux,
+ tilde_s_flux, tilde_b_flux, tilde_phi_flux, *tilde_d,
+ *tilde_ye, *tilde_tau, *tilde_s, *tilde_b, *tilde_phi,
+ *lapse, *shift, interior_sqrt_det_spatial_metric,
+ interior_spatial_metric, *inv_spatial_metric,
+ interior_pressure, exterior_spatial_velocity,
+ interior_lorentz_factor, exterior_magnetic_field);
+
+ return {};
+}
+
+void Reflective::fd_ghost(
+ const gsl::not_null<Scalar<DataVector>*> rest_mass_density,
+ const gsl::not_null<Scalar<DataVector>*> electron_fraction,
+ const gsl::not_null<Scalar<DataVector>*> temperature,
+ const gsl::not_null<tnsr::I<DataVector, 3, Frame::Inertial>*>
+ lorentz_factor_times_spatial_velocity,
+ const gsl::not_null<tnsr::I<DataVector, 3, Frame::Inertial>*>
+ magnetic_field,
+ const gsl::not_null<Scalar<DataVector>*> divergence_cleaning_field,
+
+ const gsl::not_null<std::optional<Variables<db::wrap_tags_in<
+ Flux, typename grmhd::ValenciaDivClean::System::flux_variables>>>*>
+ cell_centered_ghost_fluxes,
+
+ const Direction<3>& direction,
+
+ // fd_interior_temporary_tags
+ const Mesh<3>& subcell_mesh,
+ const tnsr::I<DataVector, 3, Frame::Inertial>& interior_shift,
+ const Scalar<DataVector>& interior_lapse,
+ const tnsr::ii<DataVector, 3, Frame::Inertial>& interior_spatial_metric,
+
+ // fd_interior_primitive_variables_tags
+ const Scalar<DataVector>& interior_rest_mass_density,
+ const Scalar<DataVector>& interior_electron_fraction,
+ const Scalar<DataVector>& interior_temperature,
+ const Scalar<DataVector>& interior_pressure,
+ const Scalar<DataVector>& interior_specific_internal_energy,
+ const Scalar<DataVector>& interior_lorentz_factor,
+ const tnsr::I<DataVector, 3, Frame::Inertial>& interior_spatial_velocity,
+ const tnsr::I<DataVector, 3, Frame::Inertial>& interior_magnetic_field,
+
+ // fd_gridless_tags
+ const fd::Reconstructor& reconstructor) const {
+ Variables<tmpl::push_back<typename System::variables_tag::tags_list,
+ SpatialVelocity, LorentzFactor, Pressure,
+ SpecificInternalEnergy, SqrtDetSpatialMetric,
+ SpatialMetric, InvSpatialMetric, Lapse, Shift>>
+ temp_vars{get(*rest_mass_density).size()};
+ fd_ghost_impl(rest_mass_density, electron_fraction, temperature,
+ make_not_null(&get<Pressure>(temp_vars)),
+ make_not_null(&get<SpecificInternalEnergy>(temp_vars)),
+ lorentz_factor_times_spatial_velocity,
+ make_not_null(&get<SpatialVelocity>(temp_vars)),
+ make_not_null(&get<LorentzFactor>(temp_vars)), magnetic_field,
+ divergence_cleaning_field,
+ make_not_null(&get<SpatialMetric>(temp_vars)),
+ make_not_null(&get<InvSpatialMetric>(temp_vars)),
+ make_not_null(&get<SqrtDetSpatialMetric>(temp_vars)),
+ make_not_null(&get<Lapse>(temp_vars)),
+ make_not_null(&get<Shift>(temp_vars)),
+
+ direction,
+
+ // fd_interior_temporary_tags
+ subcell_mesh,
+
+ // fd_interior_primitive_variables_tags
+ interior_rest_mass_density, interior_electron_fraction,
+ interior_temperature, interior_pressure,
+ interior_specific_internal_energy, interior_lorentz_factor,
+ interior_spatial_velocity, interior_magnetic_field,
+ interior_spatial_metric, interior_lapse, interior_shift,
+
+ reconstructor.ghost_zone_size(),
+
+ cell_centered_ghost_fluxes->has_value());
+ if (cell_centered_ghost_fluxes->has_value()) {
+ ConservativeFromPrimitive::apply(
+ make_not_null(&get<Tags::TildeD>(temp_vars)),
+ make_not_null(&get<Tags::TildeYe>(temp_vars)),
+ make_not_null(&get<Tags::TildeTau>(temp_vars)),
+ make_not_null(&get<Tags::TildeS<>>(temp_vars)),
+ make_not_null(&get<Tags::TildeB<>>(temp_vars)),
+ make_not_null(&get<Tags::TildePhi>(temp_vars)),
+
+ // Note: Only the spatial velocity changes.
+ *rest_mass_density, *electron_fraction,
+ get<SpecificInternalEnergy>(temp_vars), get<Pressure>(temp_vars),
+ get<SpatialVelocity>(temp_vars), get<LorentzFactor>(temp_vars),
+ *magnetic_field,
+
+ get<SqrtDetSpatialMetric>(temp_vars), get<SpatialMetric>(temp_vars),
+ *divergence_cleaning_field);
+
+ ComputeFluxes::apply(
+ make_not_null(
+ &get<Flux<Tags::TildeD>>(cell_centered_ghost_fluxes->value())),
+ make_not_null(
+ &get<Flux<Tags::TildeYe>>(cell_centered_ghost_fluxes->value())),
+ make_not_null(
+ &get<Flux<Tags::TildeTau>>(cell_centered_ghost_fluxes->value())),
+ make_not_null(
+ &get<Flux<Tags::TildeS<>>>(cell_centered_ghost_fluxes->value())),
+ make_not_null(
+ &get<Flux<Tags::TildeB<>>>(cell_centered_ghost_fluxes->value())),
+ make_not_null(
+ &get<Flux<Tags::TildePhi>>(cell_centered_ghost_fluxes->value())),
+
+ get<Tags::TildeD>(temp_vars), get<Tags::TildeYe>(temp_vars),
+ get<Tags::TildeTau>(temp_vars), get<Tags::TildeS<>>(temp_vars),
+ get<Tags::TildeB<>>(temp_vars), get<Tags::TildePhi>(temp_vars),
+
+ get<Lapse>(temp_vars), get<Shift>(temp_vars),
+ get<SqrtDetSpatialMetric>(temp_vars), get<SpatialMetric>(temp_vars),
+ get<InvSpatialMetric>(temp_vars), get<Pressure>(temp_vars),
+ get<SpatialVelocity>(temp_vars), get<LorentzFactor>(temp_vars),
+ *magnetic_field);
+ }
+}
+
+void Reflective::fd_ghost_impl(
+ const gsl::not_null<Scalar<DataVector>*> rest_mass_density,
+ const gsl::not_null<Scalar<DataVector>*> electron_fraction,
+ const gsl::not_null<Scalar<DataVector>*> temperature,
+ const gsl::not_null<Scalar<DataVector>*> pressure,
+ const gsl::not_null<Scalar<DataVector>*> specific_internal_energy,
+ const gsl::not_null<tnsr::I<DataVector, 3, Frame::Inertial>*>
+ lorentz_factor_times_spatial_velocity,
+ const gsl::not_null<tnsr::I<DataVector, 3, Frame::Inertial>*>
+ spatial_velocity,
+ const gsl::not_null<Scalar<DataVector>*> lorentz_factor,
+ const gsl::not_null<tnsr::I<DataVector, 3, Frame::Inertial>*>
+ magnetic_field,
+ const gsl::not_null<Scalar<DataVector>*> divergence_cleaning_field,
+ const gsl::not_null<tnsr::ii<DataVector, 3, Frame::Inertial>*>
+ spatial_metric,
+ const gsl::not_null<tnsr::II<DataVector, 3, Frame::Inertial>*>
+ inv_spatial_metric,
+ const gsl::not_null<Scalar<DataVector>*> sqrt_det_spatial_metric,
+ const gsl::not_null<Scalar<DataVector>*> lapse,
+ const gsl::not_null<tnsr::I<DataVector, 3, Frame::Inertial>*> shift,
+
+ const Direction<3>& direction,
+
+ // fd_interior_temporary_tags
+ const Mesh<3>& subcell_mesh,
+
+ // fd_interior_primitive_variables_tags
+ const Scalar<DataVector>& interior_rest_mass_density,
+ const Scalar<DataVector>& interior_electron_fraction,
+ const Scalar<DataVector>& interior_temperature,
+ const Scalar<DataVector>& interior_pressure,
+ const Scalar<DataVector>& interior_specific_internal_energy,
+ const Scalar<DataVector>& interior_lorentz_factor,
+ const tnsr::I<DataVector, 3, Frame::Inertial>& interior_spatial_velocity,
+ const tnsr::I<DataVector, 3, Frame::Inertial>& interior_magnetic_field,
+ const tnsr::ii<DataVector, 3, Frame::Inertial>& interior_spatial_metric,
+ const Scalar<DataVector>& interior_lapse,
+ const tnsr::I<DataVector, 3, Frame::Inertial>& interior_shift,
+
+ const size_t ghost_zone_size, const bool need_tags_for_fluxes) const {
+ const size_t dim_direction{direction.dimension()};
+
+ const auto subcell_extents{subcell_mesh.extents()};
+ const size_t num_face_pts{
+ subcell_extents.slice_away(dim_direction).product()};
+
+ using prim_tags_without_cleaning_field =
+ tmpl::list<RestMassDensity, ElectronFraction, Temperature,
+ LorentzFactorTimesSpatialVelocity, MagneticField>;
+
+ // Create a single large DV to reduce the number of Variables allocations
+ using fluxes_tags =
+ tmpl::list<Pressure, SpecificInternalEnergy, SpatialMetric, Lapse, Shift>;
+ const size_t buffer_size_for_fluxes =
+ need_tags_for_fluxes
+ ? Variables<fluxes_tags>::number_of_independent_components
+ : 0;
+ const size_t buffer_size_per_grid_pts = Variables<
+ prim_tags_without_cleaning_field>::number_of_independent_components;
+ DataVector buffer_for_vars{
+ num_face_pts * ((1 + ghost_zone_size) *
+ (buffer_size_per_grid_pts + buffer_size_for_fluxes)),
+ 0.0};
+
+ // Assign two Variables object to the buffer
+ Variables<prim_tags_without_cleaning_field> outermost_prim_vars{
+ buffer_for_vars.data(), num_face_pts * buffer_size_per_grid_pts};
+ Variables<prim_tags_without_cleaning_field> ghost_prim_vars{
+ outermost_prim_vars.data() + outermost_prim_vars.size(),
+ num_face_pts * buffer_size_per_grid_pts * ghost_zone_size};
+
+ // Slice values on the outermost grid points and store them to
+ // `outermost_prim_vars`
+ auto get_boundary_val = [&direction, &subcell_extents](auto volume_tensor) {
+ return evolution::dg::subcell::slice_tensor_for_subcell(
+ volume_tensor, subcell_extents, 1, direction, {});
+ };
+
+ get<RestMassDensity>(outermost_prim_vars) =
+ get_boundary_val(interior_rest_mass_density);
+ get<ElectronFraction>(outermost_prim_vars) =
+ get_boundary_val(interior_electron_fraction);
+ get<Temperature>(outermost_prim_vars) =
+ get_boundary_val(interior_temperature);
+
+ {
+ const auto normal_spatial_velocity_at_boundary =
+ get_boundary_val(interior_spatial_velocity).get(dim_direction);
+ // reflect both the outgoing and ingoing normal components of
+ // spatial velocity and the magnetic field.
+ if (reflect_both_) {
+ for (size_t i = 0; i < 3; ++i) {
+ if (i == dim_direction) {
+ get<LorentzFactorTimesSpatialVelocity>(outermost_prim_vars).get(i) =
+ -get(get_boundary_val(interior_lorentz_factor)) *
+ get_boundary_val(interior_spatial_velocity).get(i);
+ get<MagneticField>(outermost_prim_vars).get(i) =
+ -1.0 * get_boundary_val(interior_magnetic_field).get(i);
+ } else {
+ get<LorentzFactorTimesSpatialVelocity>(outermost_prim_vars).get(i) =
+ get(get_boundary_val(interior_lorentz_factor)) *
+ get_boundary_val(interior_spatial_velocity).get(i);
+ get<MagneticField>(outermost_prim_vars).get(i) =
+ get_boundary_val(interior_magnetic_field).get(i);
+ }
+ }
+ }
+ // reflect only the outgoing normal component of spatial
+ // velocity and the magnetic field.
+ else {
+ for (size_t i = 0; i < 3; ++i) {
+ if (i == dim_direction) {
+ if (direction.sign() > 0.0) {
+ get<LorentzFactorTimesSpatialVelocity>(outermost_prim_vars).get(i) =
+ get(get_boundary_val(interior_lorentz_factor)) *
+ min(normal_spatial_velocity_at_boundary,
+ -normal_spatial_velocity_at_boundary);
+ get<MagneticField>(outermost_prim_vars).get(i) =
+ min(get_boundary_val(interior_magnetic_field).get(i),
+ -1.0 * get_boundary_val(interior_magnetic_field).get(i));
+ } else {
+ get<LorentzFactorTimesSpatialVelocity>(outermost_prim_vars).get(i) =
+ get(get_boundary_val(interior_lorentz_factor)) *
+ max(normal_spatial_velocity_at_boundary,
+ -normal_spatial_velocity_at_boundary);
+ get<MagneticField>(outermost_prim_vars).get(i) =
+ max(get_boundary_val(interior_magnetic_field).get(i),
+ -1.0 * get_boundary_val(interior_magnetic_field).get(i));
+ }
+ } else {
+ get<LorentzFactorTimesSpatialVelocity>(outermost_prim_vars).get(i) =
+ get(get_boundary_val(interior_lorentz_factor)) *
+ get_boundary_val(interior_spatial_velocity).get(i);
+ get<MagneticField>(outermost_prim_vars).get(i) =
+ get_boundary_val(interior_magnetic_field).get(i);
+ }
+ }
+ }
+ }
+
+ // Now copy `outermost_prim_vars` into each slices of `ghost_prim_vars`.
+ Index<3> ghost_data_extents = subcell_extents;
+ ghost_data_extents[dim_direction] = ghost_zone_size;
+
+ for (size_t i_ghost = 0; i_ghost < ghost_zone_size; ++i_ghost) {
+ add_slice_to_data(make_not_null(&ghost_prim_vars), outermost_prim_vars,
+ ghost_data_extents, dim_direction, i_ghost);
+ }
+
+ *rest_mass_density = get<RestMassDensity>(ghost_prim_vars);
+ *electron_fraction = get<ElectronFraction>(ghost_prim_vars);
+ *temperature = get<Temperature>(ghost_prim_vars);
+ *lorentz_factor_times_spatial_velocity =
+ get<LorentzFactorTimesSpatialVelocity>(ghost_prim_vars);
+ *magnetic_field = get<MagneticField>(ghost_prim_vars);
+
+ // divergence cleaning scalar field is just set to zero in the ghost zone.
+ get(*divergence_cleaning_field) = 0.0;
+
+ if (need_tags_for_fluxes) {
+ Variables<fluxes_tags> outermost_fluxes_vars{
+ std::next(ghost_prim_vars.data(),
+ static_cast<std::ptrdiff_t>(ghost_prim_vars.size())),
+ num_face_pts * buffer_size_for_fluxes};
+ Variables<fluxes_tags> ghost_fluxes_vars{
+ std::next(outermost_fluxes_vars.data(),
+ static_cast<std::ptrdiff_t>(outermost_fluxes_vars.size())),
+ num_face_pts * buffer_size_for_fluxes * ghost_zone_size};
+
+ get<Pressure>(outermost_fluxes_vars) = get_boundary_val(interior_pressure);
+ get<SpecificInternalEnergy>(outermost_fluxes_vars) =
+ get_boundary_val(interior_specific_internal_energy);
+ get<SpatialMetric>(outermost_fluxes_vars) =
+ get_boundary_val(interior_spatial_metric);
+ get<Lapse>(outermost_fluxes_vars) = get_boundary_val(interior_lapse);
+ get<Shift>(outermost_fluxes_vars) = get_boundary_val(interior_shift);
+
+ for (size_t i_ghost = 0; i_ghost < ghost_zone_size; ++i_ghost) {
+ add_slice_to_data(make_not_null(&ghost_fluxes_vars),
+ outermost_fluxes_vars, ghost_data_extents,
+ dim_direction, i_ghost);
+ }
+ // Need pressure for high-order finite difference
+ *pressure = get<Pressure>(ghost_fluxes_vars);
+ *specific_internal_energy = get<SpecificInternalEnergy>(ghost_fluxes_vars);
+ *spatial_metric = get<SpatialMetric>(ghost_fluxes_vars);
+ *lapse = get<Lapse>(ghost_fluxes_vars);
+ *shift = get<Shift>(ghost_fluxes_vars);
+
+ determinant_and_inverse(sqrt_det_spatial_metric, inv_spatial_metric,
+ *spatial_metric);
+ get(*sqrt_det_spatial_metric) = sqrt(get(*sqrt_det_spatial_metric));
+ tenex::evaluate(
+ lorentz_factor,
+ sqrt(1.0 + (*spatial_metric)(ti::i, ti::j) *
+ (*lorentz_factor_times_spatial_velocity)(ti::I) *
+ (*lorentz_factor_times_spatial_velocity)(ti::J)));
+ tenex::evaluate<ti::I>(
+ spatial_velocity,
+ (*lorentz_factor_times_spatial_velocity)(ti::I) / (*lorentz_factor)());
+ }
+}
+} // namespace grmhd::ValenciaDivClean::BoundaryConditions
diff --git a/src/Evolution/Systems/GrMhd/ValenciaDivClean/BoundaryConditions/Reflective.hpp b/src/Evolution/Systems/GrMhd/ValenciaDivClean/BoundaryConditions/Reflective.hpp
new file mode 100644
index 000000000000..341eeb5298ed
--- /dev/null
+++ b/src/Evolution/Systems/GrMhd/ValenciaDivClean/BoundaryConditions/Reflective.hpp
@@ -0,0 +1,271 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#pragma once
+
+#include <memory>
+#include <optional>
+#include <pup.h>
+#include <string>
+
+#include "DataStructures/DataBox/PrefixHelpers.hpp"
+#include "DataStructures/DataBox/Prefixes.hpp"
+#include "DataStructures/Tensor/TypeAliases.hpp"
+#include "Domain/BoundaryConditions/BoundaryCondition.hpp"
+#include "Evolution/BoundaryConditions/Type.hpp"
+#include "Evolution/DgSubcell/Tags/Mesh.hpp"
+#include "Evolution/Systems/GrMhd/ValenciaDivClean/BoundaryConditions/BoundaryCondition.hpp"
+#include "Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/Tag.hpp"
+#include "Evolution/Systems/GrMhd/ValenciaDivClean/System.hpp"
+#include "Options/String.hpp"
+#include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
+#include "PointwiseFunctions/Hydro/Tags.hpp"
+#include "PointwiseFunctions/Hydro/Temperature.hpp"
+#include "Utilities/TMPL.hpp"
+
+/// \cond
+class DataVector;
+template <size_t VolumeDim>
+class Direction;
+namespace gsl {
+template <typename T>
+class not_null;
+} // namespace gsl
+/// \endcond
+
+namespace grmhd::ValenciaDivClean::BoundaryConditions {
+/*!
+ * \brief Apply "soft" reflective boundary condition as described in
+ * \cite Shiokawa2011ih.
+ *
+ * All primitive variables at the boundary are copied into ghost zone except :
+ *
+ * - If \f$n_iv^i \leq 0\f$ where \f$v^i\f$ is spatial velocity and \f$n_i\f$
+ * is outward directed normal covector, copy the values of \f$v^i\f$ at the
+ * boundary to ghost zone. If \f$n_iv^i>0\f$, spatial velocity in the ghost zone
+ * is modified such that the sign of normal component is inverted at the
+ * interface i.e. \f$v_\text{ghost}^i = v^i - 2*(n_jv^j)n^i\f$.
+ *
+ * - If \f$n_iB^i \leq 0\f$ where \f$B^i\f$ is magnetic field and \f$n_i\f$
+ * is outward directed normal covector, copy the values of \f$B^i\f$ at the
+ * boundary to ghost zone. If \f$n_iB^i>0\f$, magnetic field in the ghost zone
+ * is modified such that the sign of normal component is inverted at the
+ * interface i.e. \f$B_\text{ghost}^i = B^i - 2*(n_jB^j)n^i\f$.
+ *
+ * - If reflect_both is true, then spatial velocity and magnetic field are
+ * are inverted regardless of whether the normal component is pointing outward
+ * or inward.
+ *
+ * - However, regardless of whether the normal component of the spatial
+ * velocity $n_iv^i$ is pointing outward or inward, the lorentz factor \f$W\f$
+ * is copied into ghost zone without any modification.
+ *
+ * - Divergence cleaning scalar field \f$\Phi\f$ is set to zero in ghost zone.
+ *
+ */
+class Reflective final : public BoundaryCondition {
+ private:
+ bool reflect_both_{false};
+
+ using RestMassDensity = hydro::Tags::RestMassDensity<DataVector>;
+ using ElectronFraction = hydro::Tags::ElectronFraction<DataVector>;
+ using Temperature = hydro::Tags::Temperature<DataVector>;
+ using Pressure = hydro::Tags::Pressure<DataVector>;
+ using LorentzFactorTimesSpatialVelocity =
+ hydro::Tags::LorentzFactorTimesSpatialVelocity<DataVector, 3>;
+ using MagneticField = hydro::Tags::MagneticField<DataVector, 3>;
+ using DivergenceCleaningField =
+ hydro::Tags::DivergenceCleaningField<DataVector>;
+ using SpecificInternalEnergy =
+ hydro::Tags::SpecificInternalEnergy<DataVector>;
+ using SpatialVelocity = hydro::Tags::SpatialVelocity<DataVector, 3>;
+ using LorentzFactor = hydro::Tags::LorentzFactor<DataVector>;
+ using SqrtDetSpatialMetric = gr::Tags::SqrtDetSpatialMetric<DataVector>;
+ using SpatialMetric = gr::Tags::SpatialMetric<DataVector, 3>;
+ using InvSpatialMetric = gr::Tags::InverseSpatialMetric<DataVector, 3>;
+ using Lapse = gr::Tags::Lapse<DataVector>;
+ using Shift = gr::Tags::Shift<DataVector, 3>;
+
+ using prim_tags_for_reconstruction =
+ tmpl::list<RestMassDensity, ElectronFraction, Temperature,
+ LorentzFactorTimesSpatialVelocity, MagneticField,
+ DivergenceCleaningField>;
+
+ template <typename T>
+ using Flux = ::Tags::Flux<T, tmpl::size_t<3>, Frame::Inertial>;
+
+ public:
+ struct ReflectBoth {
+ using type = bool;
+ static constexpr Options::String help = {
+ "Reflect both outgoing and incoming normal component of "
+ "spatial velocity and magnetic field."};
+ };
+ using options = tmpl::list<ReflectBoth>;
+ static constexpr Options::String help{
+ "Reflective boundary conditions, inverting the sign "
+ "of outgoing normal component of spatial velocity "
+ "and magnetic field."};
+
+ Reflective() = default;
+ Reflective(Reflective&&) = default;
+ Reflective& operator=(Reflective&&) = default;
+ Reflective(const Reflective&) = default;
+ Reflective& operator=(const Reflective&) = default;
+ ~Reflective() override = default;
+
+ explicit Reflective(bool reflect_both);
+
+ explicit Reflective(CkMigrateMessage* msg);
+
+ WRAPPED_PUPable_decl_base_template(
+ domain::BoundaryConditions::BoundaryCondition, Reflective);
+
+ auto get_clone() const -> std::unique_ptr<
+ domain::BoundaryConditions::BoundaryCondition> override;
+
+ static constexpr evolution::BoundaryConditions::Type bc_type =
+ evolution::BoundaryConditions::Type::Ghost;
+
+ void pup(PUP::er& p) override;
+
+ using dg_interior_evolved_variables_tags = tmpl::list<>;
+ using dg_interior_primitive_variables_tags =
+ tmpl::list<hydro::Tags::RestMassDensity<DataVector>,
+ hydro::Tags::ElectronFraction<DataVector>,
+ hydro::Tags::SpecificInternalEnergy<DataVector>,
+ hydro::Tags::SpatialVelocity<DataVector, 3>,
+ hydro::Tags::MagneticField<DataVector, 3>,
+ hydro::Tags::LorentzFactor<DataVector>,
+ hydro::Tags::Pressure<DataVector>>;
+ using dg_interior_temporary_tags = tmpl::list<Shift, Lapse, InvSpatialMetric>;
+ using dg_gridless_tags = tmpl::list<>;
+
+ std::optional<std::string> dg_ghost(
+ gsl::not_null<Scalar<DataVector>*> tilde_d,
+ gsl::not_null<Scalar<DataVector>*> tilde_ye,
+ gsl::not_null<Scalar<DataVector>*> tilde_tau,
+ gsl::not_null<tnsr::i<DataVector, 3, Frame::Inertial>*> tilde_s,
+ gsl::not_null<tnsr::I<DataVector, 3, Frame::Inertial>*> tilde_b,
+ gsl::not_null<Scalar<DataVector>*> tilde_phi,
+
+ gsl::not_null<tnsr::I<DataVector, 3, Frame::Inertial>*> tilde_d_flux,
+ gsl::not_null<tnsr::I<DataVector, 3, Frame::Inertial>*> tilde_ye_flux,
+ gsl::not_null<tnsr::I<DataVector, 3, Frame::Inertial>*> tilde_tau_flux,
+ gsl::not_null<tnsr::Ij<DataVector, 3, Frame::Inertial>*> tilde_s_flux,
+ gsl::not_null<tnsr::IJ<DataVector, 3, Frame::Inertial>*> tilde_b_flux,
+ gsl::not_null<tnsr::I<DataVector, 3, Frame::Inertial>*> tilde_phi_flux,
+
+ gsl::not_null<Scalar<DataVector>*> lapse,
+ gsl::not_null<tnsr::I<DataVector, 3, Frame::Inertial>*> shift,
+ gsl::not_null<tnsr::II<DataVector, 3, Frame::Inertial>*>
+ inv_spatial_metric,
+
+ const std::optional<tnsr::I<DataVector, 3, Frame::Inertial>>&
+ /*face_mesh_velocity*/,
+ const tnsr::i<DataVector, 3, Frame::Inertial>&
+ outward_directed_normal_covector,
+ const tnsr::I<DataVector, 3, Frame::Inertial>&
+ outward_directed_normal_vector,
+
+ const Scalar<DataVector>& interior_rest_mass_density,
+ const Scalar<DataVector>& interior_electron_fraction,
+ const Scalar<DataVector>& interior_specific_internal_energy,
+ const tnsr::I<DataVector, 3, Frame::Inertial>& interior_spatial_velocity,
+ const tnsr::I<DataVector, 3, Frame::Inertial>& interior_magnetic_field,
+ const Scalar<DataVector>& interior_lorentz_factor,
+ const Scalar<DataVector>& interior_pressure,
+
+ const tnsr::I<DataVector, 3, Frame::Inertial>& interior_shift,
+ const Scalar<DataVector>& interior_lapse,
+ const tnsr::II<DataVector, 3, Frame::Inertial>&
+ interior_inv_spatial_metric) const;
+
+ using fd_interior_evolved_variables_tags = tmpl::list<>;
+ using fd_interior_temporary_tags =
+ tmpl::list<evolution::dg::subcell::Tags::Mesh<3>, Shift, Lapse,
+ SpatialMetric>;
+ using fd_interior_primitive_variables_tags =
+ tmpl::list<RestMassDensity, ElectronFraction, Temperature,
+ hydro::Tags::Pressure<DataVector>,
+ hydro::Tags::SpecificInternalEnergy<DataVector>,
+ hydro::Tags::LorentzFactor<DataVector>,
+ hydro::Tags::SpatialVelocity<DataVector, 3>, MagneticField>;
+
+ using fd_gridless_tags = tmpl::list<fd::Tags::Reconstructor>;
+
+ void fd_ghost(
+ gsl::not_null<Scalar<DataVector>*> rest_mass_density,
+ gsl::not_null<Scalar<DataVector>*> electron_fraction,
+ gsl::not_null<Scalar<DataVector>*> temperature,
+ gsl::not_null<tnsr::I<DataVector, 3, Frame::Inertial>*>
+ lorentz_factor_times_spatial_velocity,
+ gsl::not_null<tnsr::I<DataVector, 3, Frame::Inertial>*> magnetic_field,
+ gsl::not_null<Scalar<DataVector>*> divergence_cleaning_field,
+
+ gsl::not_null<std::optional<Variables<db::wrap_tags_in<
+ Flux, typename grmhd::ValenciaDivClean::System::flux_variables>>>*>
+ cell_centered_ghost_fluxes,
+
+ const Direction<3>& direction,
+
+ // interior temporary tags
+ const Mesh<3>& subcell_mesh,
+ const tnsr::I<DataVector, 3, Frame::Inertial>& interior_shift,
+ const Scalar<DataVector>& interior_lapse,
+ const tnsr::ii<DataVector, 3, Frame::Inertial>& interior_spatial_metric,
+
+ // interior prim vars tags
+ const Scalar<DataVector>& interior_rest_mass_density,
+ const Scalar<DataVector>& interior_electron_fraction,
+ const Scalar<DataVector>& interior_temperature,
+ const Scalar<DataVector>& interior_pressure,
+ const Scalar<DataVector>& interior_specific_internal_energy,
+ const Scalar<DataVector>& interior_lorentz_factor,
+ const tnsr::I<DataVector, 3, Frame::Inertial>& interior_spatial_velocity,
+ const tnsr::I<DataVector, 3, Frame::Inertial>& interior_magnetic_field,
+
+ // gridless tags
+ const fd::Reconstructor& reconstructor) const;
+
+ // have an impl to make sharing code with GH+GRMHD easy
+ void fd_ghost_impl(
+ gsl::not_null<Scalar<DataVector>*> rest_mass_density,
+ gsl::not_null<Scalar<DataVector>*> electron_fraction,
+ gsl::not_null<Scalar<DataVector>*> temperature,
+ gsl::not_null<Scalar<DataVector>*> pressure,
+ gsl::not_null<Scalar<DataVector>*> specific_internal_energy,
+ gsl::not_null<tnsr::I<DataVector, 3, Frame::Inertial>*>
+ lorentz_factor_times_spatial_velocity,
+ gsl::not_null<tnsr::I<DataVector, 3, Frame::Inertial>*> spatial_velocity,
+ gsl::not_null<Scalar<DataVector>*> lorentz_factor,
+ gsl::not_null<tnsr::I<DataVector, 3, Frame::Inertial>*> magnetic_field,
+ gsl::not_null<Scalar<DataVector>*> divergence_cleaning_field,
+ gsl::not_null<tnsr::ii<DataVector, 3, Frame::Inertial>*> spatial_metric,
+ gsl::not_null<tnsr::II<DataVector, 3, Frame::Inertial>*>
+ inv_spatial_metric,
+ gsl::not_null<Scalar<DataVector>*> sqrt_det_spatial_metric,
+ gsl::not_null<Scalar<DataVector>*> lapse,
+ gsl::not_null<tnsr::I<DataVector, 3, Frame::Inertial>*> shift,
+
+ const Direction<3>& direction,
+
+ // fd_interior_temporary_tags
+ const Mesh<3>& subcell_mesh,
+
+ // fd_interior_primitive_variables_tags
+ const Scalar<DataVector>& interior_rest_mass_density,
+ const Scalar<DataVector>& interior_electron_fraction,
+ const Scalar<DataVector>& interior_temperature,
+ const Scalar<DataVector>& interior_pressure,
+ const Scalar<DataVector>& interior_specific_internal_energy,
+ const Scalar<DataVector>& interior_lorentz_factor,
+ const tnsr::I<DataVector, 3, Frame::Inertial>& interior_spatial_velocity,
+ const tnsr::I<DataVector, 3, Frame::Inertial>& interior_magnetic_field,
+ const tnsr::ii<DataVector, 3, Frame::Inertial>& interior_spatial_metric,
+ const Scalar<DataVector>& interior_lapse,
+ const tnsr::I<DataVector, 3, Frame::Inertial>& interior_shift,
+
+ size_t ghost_zone_size, bool need_tags_for_fluxes) const;
+ };
+} // namespace grmhd::ValenciaDivClean::BoundaryConditions
diff --git a/tests/Unit/Evolution/Systems/GrMhd/ValenciaDivClean/BoundaryConditions/Reflective.py b/tests/Unit/Evolution/Systems/GrMhd/ValenciaDivClean/BoundaryConditions/Reflective.py
new file mode 100644
index 000000000000..494da7f2580d
--- /dev/null
+++ b/tests/Unit/Evolution/Systems/GrMhd/ValenciaDivClean/BoundaryConditions/Reflective.py
@@ -0,0 +1,934 @@
+# Distributed under the MIT License.
+# See LICENSE.txt for details.
+
+import Evolution.Systems.GrMhd.ValenciaDivClean.Fluxes as fluxes
+import Evolution.Systems.GrMhd.ValenciaDivClean.TestFunctions as cons
+import numpy as np
+
+
+def _exterior_spatial_velocity(
+ outward_directed_normal_covector,
+ outward_directed_normal_vector,
+ interior_spatial_velocity,
+ reflect_both,
+):
+ dot_interior_spatial_velocity = np.dot(
+ outward_directed_normal_covector, interior_spatial_velocity
+ )
+
+ if reflect_both:
+ return np.where(
+ dot_interior_spatial_velocity > 0.0,
+ interior_spatial_velocity
+ - 2
+ * outward_directed_normal_vector
+ * dot_interior_spatial_velocity,
+ interior_spatial_velocity
+ + 2
+ * outward_directed_normal_vector
+ * dot_interior_spatial_velocity,
+ )
+
+ else:
+ return np.where(
+ dot_interior_spatial_velocity > 0.0,
+ interior_spatial_velocity
+ - 2
+ * outward_directed_normal_vector
+ * dot_interior_spatial_velocity,
+ interior_spatial_velocity,
+ )
+
+
+def _exterior_magnetic_field(
+ outward_directed_normal_covector,
+ outward_directed_normal_vector,
+ interior_magnetic_field,
+ reflect_both,
+):
+ dot_interior_magnetic_field = np.dot(
+ outward_directed_normal_covector, interior_magnetic_field
+ )
+ if reflect_both:
+ return np.where(
+ dot_interior_magnetic_field > 0.0,
+ interior_magnetic_field
+ - 2 * outward_directed_normal_vector * dot_interior_magnetic_field,
+ interior_magnetic_field
+ + 2 * outward_directed_normal_vector * dot_interior_magnetic_field,
+ )
+ else:
+ return np.where(
+ dot_interior_magnetic_field > 0.0,
+ interior_magnetic_field
+ - 2 * outward_directed_normal_vector * dot_interior_magnetic_field,
+ interior_magnetic_field,
+ )
+
+
+def _spatial_metric(inv_spatial_metric):
+ return np.linalg.inv(inv_spatial_metric)
+
+
+def error(
+ face_mesh_velocity,
+ outward_directed_normal_covector,
+ outward_directed_normal_vector,
+ interior_rest_mass_density,
+ interior_electron_fraction,
+ interior_specific_internal_energy,
+ interior_spatial_velocity,
+ interior_magnetic_field,
+ interior_lorentz_factor,
+ interior_pressure,
+ shift,
+ lapse,
+ inv_spatial_metric,
+ reflect_both,
+):
+ return None
+
+
+def tilde_d(
+ face_mesh_velocity,
+ outward_directed_normal_covector,
+ outward_directed_normal_vector,
+ interior_rest_mass_density,
+ interior_electron_fraction,
+ interior_specific_internal_energy,
+ interior_spatial_velocity,
+ interior_magnetic_field,
+ interior_lorentz_factor,
+ interior_pressure,
+ shift,
+ lapse,
+ inv_spatial_metric,
+ reflect_both,
+):
+ exterior_spatial_velocity = _exterior_spatial_velocity(
+ outward_directed_normal_covector,
+ outward_directed_normal_vector,
+ interior_spatial_velocity,
+ reflect_both,
+ )
+ exterior_magnetic_field = _exterior_magnetic_field(
+ outward_directed_normal_covector,
+ outward_directed_normal_vector,
+ interior_magnetic_field,
+ reflect_both,
+ )
+ spatial_metric = _spatial_metric(inv_spatial_metric)
+ sqrt_det_spatial_metric = np.sqrt(np.linalg.det(spatial_metric))
+
+ return cons.tilde_d(
+ interior_rest_mass_density,
+ interior_electron_fraction,
+ interior_specific_internal_energy,
+ interior_pressure,
+ exterior_spatial_velocity,
+ interior_lorentz_factor,
+ exterior_magnetic_field,
+ sqrt_det_spatial_metric,
+ spatial_metric,
+ 0.0,
+ )
+
+
+def tilde_ye(
+ face_mesh_velocity,
+ outward_directed_normal_covector,
+ outward_directed_normal_vector,
+ interior_rest_mass_density,
+ interior_electron_fraction,
+ interior_specific_internal_energy,
+ interior_spatial_velocity,
+ interior_magnetic_field,
+ interior_lorentz_factor,
+ interior_pressure,
+ shift,
+ lapse,
+ inv_spatial_metric,
+ reflect_both,
+):
+ exterior_spatial_velocity = _exterior_spatial_velocity(
+ outward_directed_normal_covector,
+ outward_directed_normal_vector,
+ interior_spatial_velocity,
+ reflect_both,
+ )
+ exterior_magnetic_field = _exterior_magnetic_field(
+ outward_directed_normal_covector,
+ outward_directed_normal_vector,
+ interior_magnetic_field,
+ reflect_both,
+ )
+ spatial_metric = _spatial_metric(inv_spatial_metric)
+ sqrt_det_spatial_metric = np.sqrt(np.linalg.det(spatial_metric))
+
+ return cons.tilde_ye(
+ interior_rest_mass_density,
+ interior_electron_fraction,
+ interior_specific_internal_energy,
+ interior_pressure,
+ exterior_spatial_velocity,
+ interior_lorentz_factor,
+ exterior_magnetic_field,
+ sqrt_det_spatial_metric,
+ spatial_metric,
+ 0.0,
+ )
+
+
+def tilde_tau(
+ face_mesh_velocity,
+ outward_directed_normal_covector,
+ outward_directed_normal_vector,
+ interior_rest_mass_density,
+ interior_electron_fraction,
+ interior_specific_internal_energy,
+ interior_spatial_velocity,
+ interior_magnetic_field,
+ interior_lorentz_factor,
+ interior_pressure,
+ shift,
+ lapse,
+ inv_spatial_metric,
+ reflect_both,
+):
+ exterior_spatial_velocity = _exterior_spatial_velocity(
+ outward_directed_normal_covector,
+ outward_directed_normal_vector,
+ interior_spatial_velocity,
+ reflect_both,
+ )
+ exterior_magnetic_field = _exterior_magnetic_field(
+ outward_directed_normal_covector,
+ outward_directed_normal_vector,
+ interior_magnetic_field,
+ reflect_both,
+ )
+ spatial_metric = _spatial_metric(inv_spatial_metric)
+ sqrt_det_spatial_metric = np.sqrt(np.linalg.det(spatial_metric))
+
+ return cons.tilde_tau(
+ interior_rest_mass_density,
+ interior_electron_fraction,
+ interior_specific_internal_energy,
+ interior_pressure,
+ exterior_spatial_velocity,
+ interior_lorentz_factor,
+ exterior_magnetic_field,
+ sqrt_det_spatial_metric,
+ spatial_metric,
+ 0.0,
+ )
+
+
+def tilde_s(
+ face_mesh_velocity,
+ outward_directed_normal_covector,
+ outward_directed_normal_vector,
+ interior_rest_mass_density,
+ interior_electron_fraction,
+ interior_specific_internal_energy,
+ interior_spatial_velocity,
+ interior_magnetic_field,
+ interior_lorentz_factor,
+ interior_pressure,
+ shift,
+ lapse,
+ inv_spatial_metric,
+ reflect_both,
+):
+ exterior_spatial_velocity = _exterior_spatial_velocity(
+ outward_directed_normal_covector,
+ outward_directed_normal_vector,
+ interior_spatial_velocity,
+ reflect_both,
+ )
+ exterior_magnetic_field = _exterior_magnetic_field(
+ outward_directed_normal_covector,
+ outward_directed_normal_vector,
+ interior_magnetic_field,
+ reflect_both,
+ )
+ spatial_metric = _spatial_metric(inv_spatial_metric)
+ sqrt_det_spatial_metric = np.sqrt(np.linalg.det(spatial_metric))
+
+ return cons.tilde_s(
+ interior_rest_mass_density,
+ interior_electron_fraction,
+ interior_specific_internal_energy,
+ interior_pressure,
+ exterior_spatial_velocity,
+ interior_lorentz_factor,
+ exterior_magnetic_field,
+ sqrt_det_spatial_metric,
+ spatial_metric,
+ 0.0,
+ )
+
+
+def tilde_b(
+ face_mesh_velocity,
+ outward_directed_normal_covector,
+ outward_directed_normal_vector,
+ interior_rest_mass_density,
+ interior_electron_fraction,
+ interior_specific_internal_energy,
+ interior_spatial_velocity,
+ interior_magnetic_field,
+ interior_lorentz_factor,
+ interior_pressure,
+ shift,
+ lapse,
+ inv_spatial_metric,
+ reflect_both,
+):
+ exterior_spatial_velocity = _exterior_spatial_velocity(
+ outward_directed_normal_covector,
+ outward_directed_normal_vector,
+ interior_spatial_velocity,
+ reflect_both,
+ )
+ exterior_magnetic_field = _exterior_magnetic_field(
+ outward_directed_normal_covector,
+ outward_directed_normal_vector,
+ interior_magnetic_field,
+ reflect_both,
+ )
+ spatial_metric = _spatial_metric(inv_spatial_metric)
+ sqrt_det_spatial_metric = np.sqrt(np.linalg.det(spatial_metric))
+
+ return cons.tilde_b(
+ interior_rest_mass_density,
+ interior_electron_fraction,
+ interior_specific_internal_energy,
+ interior_pressure,
+ exterior_spatial_velocity,
+ interior_lorentz_factor,
+ exterior_magnetic_field,
+ sqrt_det_spatial_metric,
+ spatial_metric,
+ 0.0,
+ )
+
+
+def tilde_phi(
+ face_mesh_velocity,
+ outward_directed_normal_covector,
+ outward_directed_normal_vector,
+ interior_rest_mass_density,
+ interior_electron_fraction,
+ interior_specific_internal_energy,
+ interior_spatial_velocity,
+ interior_magnetic_field,
+ interior_lorentz_factor,
+ interior_pressure,
+ shift,
+ lapse,
+ inv_spatial_metric,
+ reflect_both,
+):
+ exterior_spatial_velocity = _exterior_spatial_velocity(
+ outward_directed_normal_covector,
+ outward_directed_normal_vector,
+ interior_spatial_velocity,
+ reflect_both,
+ )
+ exterior_magnetic_field = _exterior_magnetic_field(
+ outward_directed_normal_covector,
+ outward_directed_normal_vector,
+ interior_magnetic_field,
+ reflect_both,
+ )
+ spatial_metric = _spatial_metric(inv_spatial_metric)
+ sqrt_det_spatial_metric = np.sqrt(np.linalg.det(spatial_metric))
+
+ return cons.tilde_phi(
+ interior_rest_mass_density,
+ interior_electron_fraction,
+ interior_specific_internal_energy,
+ interior_pressure,
+ exterior_spatial_velocity,
+ interior_lorentz_factor,
+ exterior_magnetic_field,
+ sqrt_det_spatial_metric,
+ spatial_metric,
+ 0.0,
+ )
+
+
+def _return_cons_vars(
+ face_mesh_velocity,
+ outward_directed_normal_covector,
+ outward_directed_normal_vector,
+ interior_rest_mass_density,
+ interior_electron_fraction,
+ interior_specific_internal_energy,
+ interior_spatial_velocity,
+ interior_magnetic_field,
+ interior_lorentz_factor,
+ interior_pressure,
+ shift,
+ lapse,
+ inv_spatial_metric,
+ reflect_both,
+):
+ return {
+ "tilde_d": tilde_d(
+ face_mesh_velocity,
+ outward_directed_normal_covector,
+ outward_directed_normal_vector,
+ interior_rest_mass_density,
+ interior_electron_fraction,
+ interior_specific_internal_energy,
+ interior_spatial_velocity,
+ interior_magnetic_field,
+ interior_lorentz_factor,
+ interior_pressure,
+ shift,
+ lapse,
+ inv_spatial_metric,
+ reflect_both,
+ ),
+ "tilde_ye": tilde_ye(
+ face_mesh_velocity,
+ outward_directed_normal_covector,
+ outward_directed_normal_vector,
+ interior_rest_mass_density,
+ interior_electron_fraction,
+ interior_specific_internal_energy,
+ interior_spatial_velocity,
+ interior_magnetic_field,
+ interior_lorentz_factor,
+ interior_pressure,
+ shift,
+ lapse,
+ inv_spatial_metric,
+ reflect_both,
+ ),
+ "tilde_tau": tilde_tau(
+ face_mesh_velocity,
+ outward_directed_normal_covector,
+ outward_directed_normal_vector,
+ interior_rest_mass_density,
+ interior_electron_fraction,
+ interior_specific_internal_energy,
+ interior_spatial_velocity,
+ interior_magnetic_field,
+ interior_lorentz_factor,
+ interior_pressure,
+ shift,
+ lapse,
+ inv_spatial_metric,
+ reflect_both,
+ ),
+ "tilde_s": tilde_s(
+ face_mesh_velocity,
+ outward_directed_normal_covector,
+ outward_directed_normal_vector,
+ interior_rest_mass_density,
+ interior_electron_fraction,
+ interior_specific_internal_energy,
+ interior_spatial_velocity,
+ interior_magnetic_field,
+ interior_lorentz_factor,
+ interior_pressure,
+ shift,
+ lapse,
+ inv_spatial_metric,
+ reflect_both,
+ ),
+ "tilde_b": tilde_b(
+ face_mesh_velocity,
+ outward_directed_normal_covector,
+ outward_directed_normal_vector,
+ interior_rest_mass_density,
+ interior_electron_fraction,
+ interior_specific_internal_energy,
+ interior_spatial_velocity,
+ interior_magnetic_field,
+ interior_lorentz_factor,
+ interior_pressure,
+ shift,
+ lapse,
+ inv_spatial_metric,
+ reflect_both,
+ ),
+ "tilde_phi": tilde_phi(
+ face_mesh_velocity,
+ outward_directed_normal_covector,
+ outward_directed_normal_vector,
+ interior_rest_mass_density,
+ interior_electron_fraction,
+ interior_specific_internal_energy,
+ interior_spatial_velocity,
+ interior_magnetic_field,
+ interior_lorentz_factor,
+ interior_pressure,
+ shift,
+ lapse,
+ inv_spatial_metric,
+ reflect_both,
+ ),
+ }
+
+
+def flux_tilde_d(
+ face_mesh_velocity,
+ outward_directed_normal_covector,
+ outward_directed_normal_vector,
+ interior_rest_mass_density,
+ interior_electron_fraction,
+ interior_specific_internal_energy,
+ interior_spatial_velocity,
+ interior_magnetic_field,
+ interior_lorentz_factor,
+ interior_pressure,
+ shift,
+ lapse,
+ inv_spatial_metric,
+ reflect_both,
+):
+ exterior_spatial_velocity = _exterior_spatial_velocity(
+ outward_directed_normal_covector,
+ outward_directed_normal_vector,
+ interior_spatial_velocity,
+ reflect_both,
+ )
+ exterior_magnetic_field = _exterior_magnetic_field(
+ outward_directed_normal_covector,
+ outward_directed_normal_vector,
+ interior_magnetic_field,
+ reflect_both,
+ )
+ spatial_metric = _spatial_metric(inv_spatial_metric)
+ sqrt_det_spatial_metric = np.sqrt(np.linalg.det(spatial_metric))
+
+ cons_vars = _return_cons_vars(
+ face_mesh_velocity,
+ outward_directed_normal_covector,
+ outward_directed_normal_vector,
+ interior_rest_mass_density,
+ interior_electron_fraction,
+ interior_specific_internal_energy,
+ interior_spatial_velocity,
+ interior_magnetic_field,
+ interior_lorentz_factor,
+ interior_pressure,
+ shift,
+ lapse,
+ inv_spatial_metric,
+ reflect_both,
+ )
+
+ return fluxes.tilde_d_flux(
+ cons_vars["tilde_d"],
+ cons_vars["tilde_ye"],
+ cons_vars["tilde_tau"],
+ cons_vars["tilde_s"],
+ cons_vars["tilde_b"],
+ cons_vars["tilde_phi"],
+ lapse,
+ shift,
+ sqrt_det_spatial_metric,
+ spatial_metric,
+ inv_spatial_metric,
+ interior_pressure,
+ exterior_spatial_velocity,
+ interior_lorentz_factor,
+ exterior_magnetic_field,
+ )
+
+
+def flux_tilde_ye(
+ face_mesh_velocity,
+ outward_directed_normal_covector,
+ outward_directed_normal_vector,
+ interior_rest_mass_density,
+ interior_electron_fraction,
+ interior_specific_internal_energy,
+ interior_spatial_velocity,
+ interior_magnetic_field,
+ interior_lorentz_factor,
+ interior_pressure,
+ shift,
+ lapse,
+ inv_spatial_metric,
+ reflect_both,
+):
+ exterior_spatial_velocity = _exterior_spatial_velocity(
+ outward_directed_normal_covector,
+ outward_directed_normal_vector,
+ interior_spatial_velocity,
+ reflect_both,
+ )
+ exterior_magnetic_field = _exterior_magnetic_field(
+ outward_directed_normal_covector,
+ outward_directed_normal_vector,
+ interior_magnetic_field,
+ reflect_both,
+ )
+ spatial_metric = _spatial_metric(inv_spatial_metric)
+ sqrt_det_spatial_metric = np.sqrt(np.linalg.det(spatial_metric))
+
+ cons_vars = _return_cons_vars(
+ face_mesh_velocity,
+ outward_directed_normal_covector,
+ outward_directed_normal_vector,
+ interior_rest_mass_density,
+ interior_electron_fraction,
+ interior_specific_internal_energy,
+ interior_spatial_velocity,
+ interior_magnetic_field,
+ interior_lorentz_factor,
+ interior_pressure,
+ shift,
+ lapse,
+ inv_spatial_metric,
+ reflect_both,
+ )
+
+ return fluxes.tilde_ye_flux(
+ cons_vars["tilde_d"],
+ cons_vars["tilde_ye"],
+ cons_vars["tilde_tau"],
+ cons_vars["tilde_s"],
+ cons_vars["tilde_b"],
+ cons_vars["tilde_phi"],
+ lapse,
+ shift,
+ sqrt_det_spatial_metric,
+ spatial_metric,
+ inv_spatial_metric,
+ interior_pressure,
+ exterior_spatial_velocity,
+ interior_lorentz_factor,
+ exterior_magnetic_field,
+ )
+
+
+def flux_tilde_tau(
+ face_mesh_velocity,
+ outward_directed_normal_covector,
+ outward_directed_normal_vector,
+ interior_rest_mass_density,
+ interior_electron_fraction,
+ interior_specific_internal_energy,
+ interior_spatial_velocity,
+ interior_magnetic_field,
+ interior_lorentz_factor,
+ interior_pressure,
+ shift,
+ lapse,
+ inv_spatial_metric,
+ reflect_both,
+):
+ exterior_spatial_velocity = _exterior_spatial_velocity(
+ outward_directed_normal_covector,
+ outward_directed_normal_vector,
+ interior_spatial_velocity,
+ reflect_both,
+ )
+ exterior_magnetic_field = _exterior_magnetic_field(
+ outward_directed_normal_covector,
+ outward_directed_normal_vector,
+ interior_magnetic_field,
+ reflect_both,
+ )
+ spatial_metric = _spatial_metric(inv_spatial_metric)
+ sqrt_det_spatial_metric = np.sqrt(np.linalg.det(spatial_metric))
+
+ cons_vars = _return_cons_vars(
+ face_mesh_velocity,
+ outward_directed_normal_covector,
+ outward_directed_normal_vector,
+ interior_rest_mass_density,
+ interior_electron_fraction,
+ interior_specific_internal_energy,
+ interior_spatial_velocity,
+ interior_magnetic_field,
+ interior_lorentz_factor,
+ interior_pressure,
+ shift,
+ lapse,
+ inv_spatial_metric,
+ reflect_both,
+ )
+
+ return fluxes.tilde_tau_flux(
+ cons_vars["tilde_d"],
+ cons_vars["tilde_ye"],
+ cons_vars["tilde_tau"],
+ cons_vars["tilde_s"],
+ cons_vars["tilde_b"],
+ cons_vars["tilde_phi"],
+ lapse,
+ shift,
+ sqrt_det_spatial_metric,
+ spatial_metric,
+ inv_spatial_metric,
+ interior_pressure,
+ exterior_spatial_velocity,
+ interior_lorentz_factor,
+ exterior_magnetic_field,
+ )
+
+
+def flux_tilde_s(
+ face_mesh_velocity,
+ outward_directed_normal_covector,
+ outward_directed_normal_vector,
+ interior_rest_mass_density,
+ interior_electron_fraction,
+ interior_specific_internal_energy,
+ interior_spatial_velocity,
+ interior_magnetic_field,
+ interior_lorentz_factor,
+ interior_pressure,
+ shift,
+ lapse,
+ inv_spatial_metric,
+ reflect_both,
+):
+ exterior_spatial_velocity = _exterior_spatial_velocity(
+ outward_directed_normal_covector,
+ outward_directed_normal_vector,
+ interior_spatial_velocity,
+ reflect_both,
+ )
+ exterior_magnetic_field = _exterior_magnetic_field(
+ outward_directed_normal_covector,
+ outward_directed_normal_vector,
+ interior_magnetic_field,
+ reflect_both,
+ )
+ spatial_metric = _spatial_metric(inv_spatial_metric)
+ sqrt_det_spatial_metric = np.sqrt(np.linalg.det(spatial_metric))
+
+ cons_vars = _return_cons_vars(
+ face_mesh_velocity,
+ outward_directed_normal_covector,
+ outward_directed_normal_vector,
+ interior_rest_mass_density,
+ interior_electron_fraction,
+ interior_specific_internal_energy,
+ interior_spatial_velocity,
+ interior_magnetic_field,
+ interior_lorentz_factor,
+ interior_pressure,
+ shift,
+ lapse,
+ inv_spatial_metric,
+ reflect_both,
+ )
+
+ return fluxes.tilde_s_flux(
+ cons_vars["tilde_d"],
+ cons_vars["tilde_ye"],
+ cons_vars["tilde_tau"],
+ cons_vars["tilde_s"],
+ cons_vars["tilde_b"],
+ cons_vars["tilde_phi"],
+ lapse,
+ shift,
+ sqrt_det_spatial_metric,
+ spatial_metric,
+ inv_spatial_metric,
+ interior_pressure,
+ exterior_spatial_velocity,
+ interior_lorentz_factor,
+ exterior_magnetic_field,
+ )
+
+
+def flux_tilde_b(
+ face_mesh_velocity,
+ outward_directed_normal_covector,
+ outward_directed_normal_vector,
+ interior_rest_mass_density,
+ interior_electron_fraction,
+ interior_specific_internal_energy,
+ interior_spatial_velocity,
+ interior_magnetic_field,
+ interior_lorentz_factor,
+ interior_pressure,
+ shift,
+ lapse,
+ inv_spatial_metric,
+ reflect_both,
+):
+ exterior_spatial_velocity = _exterior_spatial_velocity(
+ outward_directed_normal_covector,
+ outward_directed_normal_vector,
+ interior_spatial_velocity,
+ reflect_both,
+ )
+ exterior_magnetic_field = _exterior_magnetic_field(
+ outward_directed_normal_covector,
+ outward_directed_normal_vector,
+ interior_magnetic_field,
+ reflect_both,
+ )
+ spatial_metric = _spatial_metric(inv_spatial_metric)
+ sqrt_det_spatial_metric = np.sqrt(np.linalg.det(spatial_metric))
+
+ cons_vars = _return_cons_vars(
+ face_mesh_velocity,
+ outward_directed_normal_covector,
+ outward_directed_normal_vector,
+ interior_rest_mass_density,
+ interior_electron_fraction,
+ interior_specific_internal_energy,
+ interior_spatial_velocity,
+ interior_magnetic_field,
+ interior_lorentz_factor,
+ interior_pressure,
+ shift,
+ lapse,
+ inv_spatial_metric,
+ reflect_both,
+ )
+
+ return fluxes.tilde_b_flux(
+ cons_vars["tilde_d"],
+ cons_vars["tilde_ye"],
+ cons_vars["tilde_tau"],
+ cons_vars["tilde_s"],
+ cons_vars["tilde_b"],
+ cons_vars["tilde_phi"],
+ lapse,
+ shift,
+ sqrt_det_spatial_metric,
+ spatial_metric,
+ inv_spatial_metric,
+ interior_pressure,
+ exterior_spatial_velocity,
+ interior_lorentz_factor,
+ exterior_magnetic_field,
+ )
+
+
+def flux_tilde_phi(
+ face_mesh_velocity,
+ outward_directed_normal_covector,
+ outward_directed_normal_vector,
+ interior_rest_mass_density,
+ interior_electron_fraction,
+ interior_specific_internal_energy,
+ interior_spatial_velocity,
+ interior_magnetic_field,
+ interior_lorentz_factor,
+ interior_pressure,
+ shift,
+ lapse,
+ inv_spatial_metric,
+ reflect_both,
+):
+ exterior_spatial_velocity = _exterior_spatial_velocity(
+ outward_directed_normal_covector,
+ outward_directed_normal_vector,
+ interior_spatial_velocity,
+ reflect_both,
+ )
+ exterior_magnetic_field = _exterior_magnetic_field(
+ outward_directed_normal_covector,
+ outward_directed_normal_vector,
+ interior_magnetic_field,
+ reflect_both,
+ )
+ spatial_metric = _spatial_metric(inv_spatial_metric)
+ sqrt_det_spatial_metric = np.sqrt(np.linalg.det(spatial_metric))
+
+ cons_vars = _return_cons_vars(
+ face_mesh_velocity,
+ outward_directed_normal_covector,
+ outward_directed_normal_vector,
+ interior_rest_mass_density,
+ interior_electron_fraction,
+ interior_specific_internal_energy,
+ interior_spatial_velocity,
+ interior_magnetic_field,
+ interior_lorentz_factor,
+ interior_pressure,
+ shift,
+ lapse,
+ inv_spatial_metric,
+ reflect_both,
+ )
+
+ return fluxes.tilde_phi_flux(
+ cons_vars["tilde_d"],
+ cons_vars["tilde_ye"],
+ cons_vars["tilde_tau"],
+ cons_vars["tilde_s"],
+ cons_vars["tilde_b"],
+ cons_vars["tilde_phi"],
+ lapse,
+ shift,
+ sqrt_det_spatial_metric,
+ spatial_metric,
+ inv_spatial_metric,
+ interior_pressure,
+ exterior_spatial_velocity,
+ interior_lorentz_factor,
+ exterior_magnetic_field,
+ )
+
+
+def lapse(
+ face_mesh_velocity,
+ outward_directed_normal_covector,
+ outward_directed_normal_vector,
+ interior_rest_mass_density,
+ interior_electron_fraction,
+ interior_specific_internal_energy,
+ interior_spatial_velocity,
+ interior_magnetic_field,
+ interior_lorentz_factor,
+ interior_pressure,
+ shift,
+ lapse,
+ inv_spatial_metric,
+ reflect_both,
+):
+ return lapse
+
+
+def shift(
+ face_mesh_velocity,
+ outward_directed_normal_covector,
+ outward_directed_normal_vector,
+ interior_rest_mass_density,
+ interior_electron_fraction,
+ interior_specific_internal_energy,
+ interior_spatial_velocity,
+ interior_magnetic_field,
+ interior_lorentz_factor,
+ interior_pressure,
+ shift,
+ lapse,
+ inv_spatial_metric,
+ reflect_both,
+):
+ return shift
+
+
+def inv_spatial_metric(
+ face_mesh_velocity,
+ outward_directed_normal_covector,
+ outward_directed_normal_vector,
+ interior_rest_mass_density,
+ interior_electron_fraction,
+ interior_specific_internal_energy,
+ interior_spatial_velocity,
+ interior_magnetic_field,
+ interior_lorentz_factor,
+ interior_pressure,
+ shift,
+ lapse,
+ inv_spatial_metric,
+ reflect_both,
+):
+ return inv_spatial_metric
diff --git a/tests/Unit/Evolution/Systems/GrMhd/ValenciaDivClean/BoundaryConditions/Test_Reflective.cpp b/tests/Unit/Evolution/Systems/GrMhd/ValenciaDivClean/BoundaryConditions/Test_Reflective.cpp
new file mode 100644
index 000000000000..44efaa0ac264
--- /dev/null
+++ b/tests/Unit/Evolution/Systems/GrMhd/ValenciaDivClean/BoundaryConditions/Test_Reflective.cpp
@@ -0,0 +1,118 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#include "Framework/TestingFramework.hpp"
+
+#include <array>
+
+#include "DataStructures/DataBox/DataBox.hpp"
+#include "DataStructures/DataVector.hpp"
+#include "DataStructures/Index.hpp"
+#include "DataStructures/Tensor/EagerMath/Determinant.hpp"
+#include "Evolution/Systems/GrMhd/ValenciaDivClean/BoundaryConditions/BoundaryCondition.hpp"
+#include "Evolution/Systems/GrMhd/ValenciaDivClean/BoundaryConditions/Reflective.hpp"
+#include "Evolution/Systems/GrMhd/ValenciaDivClean/BoundaryCorrections/Rusanov.hpp"
+#include "Evolution/Systems/GrMhd/ValenciaDivClean/System.hpp"
+#include "Evolution/Systems/GrMhd/ValenciaDivClean/Tags.hpp"
+#include "Framework/SetupLocalPythonEnvironment.hpp"
+#include "Framework/TestHelpers.hpp"
+#include "Helpers/DataStructures/MakeWithRandomValues.hpp"
+#include "Helpers/Evolution/DiscontinuousGalerkin/BoundaryConditions.hpp"
+#include "Helpers/PointwiseFunctions/GeneralRelativity/TestHelpers.hpp"
+#include "Utilities/Gsl.hpp"
+#include "Utilities/TMPL.hpp"
+#include "Utilities/TaggedTuple.hpp"
+
+namespace helpers = ::TestHelpers::evolution::dg;
+namespace {
+void test_stuffs(const bool reflect_both) {
+ MAKE_GENERATOR(gen);
+
+ const auto face_mesh_index = Index<2>{2};
+ const DataVector used_for_size{face_mesh_index.product()};
+ std::uniform_real_distribution<> dist(0.5, 1.0);
+
+ const auto spatial_metric =
+ ::TestHelpers::gr::random_spatial_metric<3, DataVector, Frame::Inertial>(
+ make_not_null(&gen), used_for_size);
+ const auto sqrt_det_spatial_metric = determinant(spatial_metric);
+
+ struct ReflectBoth : db::SimpleTag {
+ using type = bool;
+ };
+
+ const auto box_with_gridless_tags =
+ db::create<db::AddSimpleTags<gr::Tags::SpatialMetric<DataVector, 3>,
+ gr::Tags::SqrtDetSpatialMetric<DataVector>,
+ ReflectBoth>>(
+ spatial_metric, sqrt_det_spatial_metric, reflect_both);
+
+ // for factory string
+ const std::string reflect_both_str = (reflect_both) ? "true" : "false";
+
+ helpers::test_boundary_condition_with_python<
+ grmhd::ValenciaDivClean::BoundaryConditions::Reflective,
+ grmhd::ValenciaDivClean::BoundaryConditions::BoundaryCondition,
+ grmhd::ValenciaDivClean::System,
+ tmpl::list<grmhd::ValenciaDivClean::BoundaryCorrections::Rusanov>,
+ tmpl::list<>, tmpl::list<ReflectBoth>>(
+ make_not_null(&gen),
+ "Evolution.Systems.GrMhd.ValenciaDivClean.BoundaryConditions."
+ "Reflective",
+ tuples::TaggedTuple<
+ helpers::Tags::PythonFunctionForErrorMessage<>,
+ helpers::Tags::PythonFunctionName<
+ grmhd::ValenciaDivClean::Tags::TildeD>,
+ helpers::Tags::PythonFunctionName<
+ grmhd::ValenciaDivClean::Tags::TildeYe>,
+ helpers::Tags::PythonFunctionName<
+ grmhd::ValenciaDivClean::Tags::TildeTau>,
+ helpers::Tags::PythonFunctionName<
+ grmhd::ValenciaDivClean::Tags::TildeS<Frame::Inertial>>,
+ helpers::Tags::PythonFunctionName<
+ grmhd::ValenciaDivClean::Tags::TildeB<Frame::Inertial>>,
+ helpers::Tags::PythonFunctionName<
+ grmhd::ValenciaDivClean::Tags::TildePhi>,
+
+ helpers::Tags::PythonFunctionName<
+ ::Tags::Flux<grmhd::ValenciaDivClean::Tags::TildeD,
+ tmpl::size_t<3>, Frame::Inertial>>,
+ helpers::Tags::PythonFunctionName<
+ ::Tags::Flux<grmhd::ValenciaDivClean::Tags::TildeYe,
+ tmpl::size_t<3>, Frame::Inertial>>,
+ helpers::Tags::PythonFunctionName<
+ ::Tags::Flux<grmhd::ValenciaDivClean::Tags::TildeTau,
+ tmpl::size_t<3>, Frame::Inertial>>,
+ helpers::Tags::PythonFunctionName<::Tags::Flux<
+ grmhd::ValenciaDivClean::Tags::TildeS<Frame::Inertial>,
+ tmpl::size_t<3>, Frame::Inertial>>,
+ helpers::Tags::PythonFunctionName<::Tags::Flux<
+ grmhd::ValenciaDivClean::Tags::TildeB<Frame::Inertial>,
+ tmpl::size_t<3>, Frame::Inertial>>,
+ helpers::Tags::PythonFunctionName<
+ ::Tags::Flux<grmhd::ValenciaDivClean::Tags::TildePhi,
+ tmpl::size_t<3>, Frame::Inertial>>,
+
+ helpers::Tags::PythonFunctionName<gr::Tags::Lapse<DataVector>>,
+ helpers::Tags::PythonFunctionName<gr::Tags::Shift<DataVector, 3>>,
+ helpers::Tags::PythonFunctionName<
+ gr::Tags::InverseSpatialMetric<DataVector, 3>>>{
+ "error", "tilde_d", "tilde_ye", "tilde_tau", "tilde_s", "tilde_b",
+ "tilde_phi", "flux_tilde_d", "flux_tilde_ye", "flux_tilde_tau",
+ "flux_tilde_s", "flux_tilde_b", "flux_tilde_phi", "lapse", "shift",
+ "inv_spatial_metric"},
+ "Reflective:\n"
+ " ReflectBoth: " +
+ reflect_both_str + "\n",
+ face_mesh_index, box_with_gridless_tags, tuples::TaggedTuple<>{},
+ 1.0e-10);
+}
+} // namespace
+
+SPECTRE_TEST_CASE("Unit.GrMhd.BoundaryConditions.Reflective", "[Unit][GrMhd]") {
+ pypp::SetupLocalPythonEnvironment local_python_env{""};
+ // Test for two cases of Reflective.reflect_both
+
+ test_stuffs(true);
+ test_stuffs(false);
+}
diff --git a/tests/Unit/Evolution/Systems/GrMhd/ValenciaDivClean/CMakeLists.txt b/tests/Unit/Evolution/Systems/GrMhd/ValenciaDivClean/CMakeLists.txt
index f9dfb390f259..3f1bd89f0d5d 100644
--- a/tests/Unit/Evolution/Systems/GrMhd/ValenciaDivClean/CMakeLists.txt
+++ b/tests/Unit/Evolution/Systems/GrMhd/ValenciaDivClean/CMakeLists.txt
@@ -9,6 +9,7 @@ set(LIBRARY_SOURCES
BoundaryConditions/Test_DemandOutgoingCharSpeeds.cpp
BoundaryConditions/Test_HydroFreeOutflow.cpp
BoundaryConditions/Test_Periodic.cpp
+ BoundaryConditions/Test_Reflective.cpp
BoundaryCorrections/Test_Hll.cpp
BoundaryCorrections/Test_Rusanov.cpp
FiniteDifference/Test_BoundaryConditionGhostData.cpp