diff --git a/src/Evolution/Systems/GrMhd/ValenciaDivClean/BoundaryConditions/Reflective.cpp b/src/Evolution/Systems/GrMhd/ValenciaDivClean/BoundaryConditions/Reflective.cpp index 63089cec8cf65..5a38377b4102d 100644 --- a/src/Evolution/Systems/GrMhd/ValenciaDivClean/BoundaryConditions/Reflective.cpp +++ b/src/Evolution/Systems/GrMhd/ValenciaDivClean/BoundaryConditions/Reflective.cpp @@ -14,6 +14,7 @@ #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" @@ -24,12 +25,14 @@ #include "Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/Reconstructor.hpp" #include "Evolution/Systems/GrMhd/ValenciaDivClean/Fluxes.hpp" #include "NumericalAlgorithms/Spectral/Mesh.hpp" -#include "PointwiseFunctions/GeneralRelativity/IndexManipulation.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 @@ -37,7 +40,10 @@ Reflective::get_clone() const { return std::make_unique(*this); } -void Reflective::pup(PUP::er& p) { BoundaryCondition::pup(p); } +void Reflective::pup(PUP::er& p) { + BoundaryCondition::pup(p); + p | reflect_both_; +} // NOLINTNEXTLINE PUP::able::PUP_ID Reflective::my_PUP_ID = 0; @@ -81,8 +87,8 @@ std::optional Reflective::dg_ghost( const tnsr::I& interior_shift, const Scalar& interior_lapse, - const tnsr::II& - interior_inv_spatial_metric) { + const tnsr::II& interior_inv_spatial_metric) + const { const size_t number_of_grid_points = get(interior_rest_mass_density).size(); // temp buffer to store @@ -139,35 +145,56 @@ std::optional Reflective::dg_ghost( outward_directed_normal_covector, interior_spatial_velocity); dot_product(make_not_null(&normal_dot_interior_magnetic_field), outward_directed_normal_covector, interior_magnetic_field); - for (size_t i = 0; i < number_of_grid_points; ++i) { - if (get(normal_dot_interior_spatial_velocity)[i] > 0.0) { + + // Reflect both outward and inward normal component of + // spatial velocity and magnetic field + if (reflect_both_) { + for (size_t i = 0; i < number_of_grid_points; ++i) { 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]; + } + } + // 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( @@ -222,7 +249,7 @@ void Reflective::fd_ghost( const tnsr::I& interior_magnetic_field, // fd_gridless_tags - const fd::Reconstructor& reconstructor) { + const fd::Reconstructor& reconstructor) const { Variables& interior_lapse, const tnsr::I& interior_shift, - const size_t ghost_zone_size, const bool need_tags_for_fluxes) { + 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()}; @@ -387,38 +414,56 @@ void Reflective::fd_ghost_impl( get_boundary_val(interior_temperature); { - // Invert the outgoing components of spatial velocity and compute Wv^i - // and the outgoing components of magnetic field. - // - // Note : Here we require the grid to be Cartesian, therefore we will need - // to change the implementation below once subcell supports curved mesh. const auto normal_spatial_velocity_at_boundary = get_boundary_val(interior_spatial_velocity).get(dim_direction); - for (size_t i = 0; i < 3; ++i) { - if (i == dim_direction) { - if (direction.sign() > 0.0) { + // Reflect both the outgoing and ingoing components of spatial + // velocity and the magnetic field. + if (reflect_both_) { + for (size_t i = 0; i < 3; ++i) { + if (i == dim_direction) { + get(outermost_prim_vars).get(i) = + -1.0 * get(get_boundary_val(interior_lorentz_factor)) * + get_boundary_val(interior_spatial_velocity).get(i); + get(outermost_prim_vars).get(i) = + -1.0 * get_boundary_val(interior_magnetic_field).get(i); + } else { get(outermost_prim_vars).get(i) = get(get_boundary_val(interior_lorentz_factor)) * - min(normal_spatial_velocity_at_boundary, - normal_spatial_velocity_at_boundary * -1.0); + get_boundary_val(interior_spatial_velocity).get(i); get(outermost_prim_vars).get(i) = - min(get_boundary_val(interior_magnetic_field).get(i), - -1.0 * get_boundary_val(interior_magnetic_field).get(i)); + get_boundary_val(interior_magnetic_field).get(i); + } + } + } + // Reflect just the outgoing components 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(outermost_prim_vars).get(i) = + get(get_boundary_val(interior_lorentz_factor)) * + min(normal_spatial_velocity_at_boundary, + normal_spatial_velocity_at_boundary * -1.0); + get(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(outermost_prim_vars).get(i) = + get(get_boundary_val(interior_lorentz_factor)) * + max(normal_spatial_velocity_at_boundary, + normal_spatial_velocity_at_boundary * -1.0); + get(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(outermost_prim_vars).get(i) = get(get_boundary_val(interior_lorentz_factor)) * - max(normal_spatial_velocity_at_boundary, - normal_spatial_velocity_at_boundary * -1.0); + get_boundary_val(interior_spatial_velocity).get(i); get(outermost_prim_vars).get(i) = - max(get_boundary_val(interior_magnetic_field).get(i), - -1.0 * get_boundary_val(interior_magnetic_field).get(i)); + get_boundary_val(interior_magnetic_field).get(i); } - } else { - get(outermost_prim_vars).get(i) = - get(get_boundary_val(interior_lorentz_factor)) * - get_boundary_val(interior_spatial_velocity).get(i); - get(outermost_prim_vars).get(i) = - get_boundary_val(interior_magnetic_field).get(i); } } } diff --git a/src/Evolution/Systems/GrMhd/ValenciaDivClean/BoundaryConditions/Reflective.hpp b/src/Evolution/Systems/GrMhd/ValenciaDivClean/BoundaryConditions/Reflective.hpp index 6031ea76bd015..ccc913c2bea3a 100644 --- a/src/Evolution/Systems/GrMhd/ValenciaDivClean/BoundaryConditions/Reflective.hpp +++ b/src/Evolution/Systems/GrMhd/ValenciaDivClean/BoundaryConditions/Reflective.hpp @@ -61,6 +61,8 @@ namespace grmhd::ValenciaDivClean::BoundaryConditions { */ class Reflective final : public BoundaryCondition { private: + bool reflect_both_{false}; + using RestMassDensity = hydro::Tags::RestMassDensity; using ElectronFraction = hydro::Tags::ElectronFraction; using Temperature = hydro::Tags::Temperature; @@ -89,7 +91,13 @@ class Reflective final : public BoundaryCondition { using Flux = ::Tags::Flux, Frame::Inertial>; public: - using options = tmpl::list<>; + 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; static constexpr Options::String help{ "Reflective boundary conditions, inverting the sign " "of outgoing normal component of spatial velocity " @@ -102,6 +110,8 @@ class Reflective final : public BoundaryCondition { Reflective& operator=(const Reflective&) = default; ~Reflective() override = default; + explicit Reflective(bool reflect_both); + explicit Reflective(CkMigrateMessage* msg); WRAPPED_PUPable_decl_base_template( @@ -127,7 +137,7 @@ class Reflective final : public BoundaryCondition { using dg_interior_temporary_tags = tmpl::list; using dg_gridless_tags = tmpl::list<>; - static std::optional dg_ghost( + std::optional dg_ghost( gsl::not_null*> tilde_d, gsl::not_null*> tilde_ye, gsl::not_null*> tilde_tau, @@ -165,7 +175,7 @@ class Reflective final : public BoundaryCondition { const tnsr::I& interior_shift, const Scalar& interior_lapse, const tnsr::II& - interior_inv_spatial_metric); + interior_inv_spatial_metric) const; using fd_interior_evolved_variables_tags = tmpl::list<>; using fd_interior_temporary_tags = @@ -180,7 +190,7 @@ class Reflective final : public BoundaryCondition { using fd_gridless_tags = tmpl::list; - static void fd_ghost( + void fd_ghost( gsl::not_null*> rest_mass_density, gsl::not_null*> electron_fraction, gsl::not_null*> temperature, @@ -212,10 +222,10 @@ class Reflective final : public BoundaryCondition { const tnsr::I& interior_magnetic_field, // gridless tags - const fd::Reconstructor& reconstructor); + const fd::Reconstructor& reconstructor) const; // have an impl to make sharing code with GH+GRMHD easy - static void fd_ghost_impl( + void fd_ghost_impl( gsl::not_null*> rest_mass_density, gsl::not_null*> electron_fraction, gsl::not_null*> temperature, @@ -252,6 +262,6 @@ class Reflective final : public BoundaryCondition { const Scalar& interior_lapse, const tnsr::I& interior_shift, - size_t ghost_zone_size, bool need_tags_for_fluxes); + 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/Test_Reflective.cpp b/tests/Unit/Evolution/Systems/GrMhd/ValenciaDivClean/BoundaryConditions/Test_Reflective.cpp index d54e579f29f22..f5f2320dce910 100644 --- a/tests/Unit/Evolution/Systems/GrMhd/ValenciaDivClean/BoundaryConditions/Test_Reflective.cpp +++ b/tests/Unit/Evolution/Systems/GrMhd/ValenciaDivClean/BoundaryConditions/Test_Reflective.cpp @@ -94,6 +94,9 @@ SPECTRE_TEST_CASE("Unit.GrMhd.BoundaryConditions.Reflective", "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", face_mesh_index, box_with_gridless_tags, - tuples::TaggedTuple<>{}, 1.0e-10); + "Reflective:\n" + " ReflectBoth: false\n", + face_mesh_index, box_with_gridless_tags, tuples::TaggedTuple<>{}, + 1.0e-10); + // FIXME: only testing ReflectBoth==false right now. }