From cb0d5f4cd7bff0f1e7bd9b4159931089036467c2 Mon Sep 17 00:00:00 2001 From: Jooheon Yoo Date: Mon, 4 Dec 2023 06:49:17 -0800 Subject: [PATCH] Add an option to reflect both outgoing and incoming (temp) --- .../BoundaryConditions/Reflective.cpp | 165 ++++++++++++------ .../BoundaryConditions/Reflective.hpp | 24 ++- .../BoundaryConditions/Reflective.py | 106 +++++++++-- .../BoundaryConditions/Test_Reflective.cpp | 23 ++- .../BoundaryConditions.hpp | 131 ++++++++++++-- 5 files changed, 358 insertions(+), 91 deletions(-) diff --git a/src/Evolution/Systems/GrMhd/ValenciaDivClean/BoundaryConditions/Reflective.cpp b/src/Evolution/Systems/GrMhd/ValenciaDivClean/BoundaryConditions/Reflective.cpp index 63089cec8cf65..fb3ec7d4c335a 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,76 @@ 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) { - 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]; + + // 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) { + 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]; + } } - } 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] + + 2.0 * get(normal_dot_interior_magnetic_field)[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]; + } + + // 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]; + } } - } 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]; + 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 +269,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 +434,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/Reflective.py b/tests/Unit/Evolution/Systems/GrMhd/ValenciaDivClean/BoundaryConditions/Reflective.py index f5221343da725..1b4f07e287765 100644 --- a/tests/Unit/Evolution/Systems/GrMhd/ValenciaDivClean/BoundaryConditions/Reflective.py +++ b/tests/Unit/Evolution/Systems/GrMhd/ValenciaDivClean/BoundaryConditions/Reflective.py @@ -10,32 +10,59 @@ 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 ) - return np.where( - dot_interior_spatial_velocity > 0.0, - interior_spatial_velocity - - 2 * outward_directed_normal_vector * dot_interior_spatial_velocity, - 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 ) - return np.where( - dot_interior_magnetic_field > 0.0, - interior_magnetic_field - - 2 * outward_directed_normal_vector * dot_interior_magnetic_field, - 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): @@ -44,6 +71,7 @@ def _spatial_metric(inv_spatial_metric): def error( face_mesh_velocity, + reflect_both, outward_directed_normal_covector, outward_directed_normal_vector, interior_rest_mass_density, @@ -62,6 +90,7 @@ def error( def tilde_d( face_mesh_velocity, + reflect_both, outward_directed_normal_covector, outward_directed_normal_vector, interior_rest_mass_density, @@ -79,11 +108,13 @@ def tilde_d( 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)) @@ -104,6 +135,7 @@ def tilde_d( def tilde_ye( face_mesh_velocity, + reflect_both, outward_directed_normal_covector, outward_directed_normal_vector, interior_rest_mass_density, @@ -121,11 +153,13 @@ def tilde_ye( 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)) @@ -146,6 +180,7 @@ def tilde_ye( def tilde_tau( face_mesh_velocity, + reflect_both, outward_directed_normal_covector, outward_directed_normal_vector, interior_rest_mass_density, @@ -163,11 +198,13 @@ def tilde_tau( 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)) @@ -188,6 +225,7 @@ def tilde_tau( def tilde_s( face_mesh_velocity, + reflect_both, outward_directed_normal_covector, outward_directed_normal_vector, interior_rest_mass_density, @@ -205,11 +243,13 @@ def tilde_s( 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)) @@ -230,6 +270,7 @@ def tilde_s( def tilde_b( face_mesh_velocity, + reflect_both, outward_directed_normal_covector, outward_directed_normal_vector, interior_rest_mass_density, @@ -247,11 +288,13 @@ def tilde_b( 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)) @@ -272,6 +315,7 @@ def tilde_b( def tilde_phi( face_mesh_velocity, + reflect_both, outward_directed_normal_covector, outward_directed_normal_vector, interior_rest_mass_density, @@ -289,11 +333,13 @@ def tilde_phi( 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)) @@ -314,6 +360,7 @@ def tilde_phi( def _return_cons_vars( face_mesh_velocity, + reflect_both, outward_directed_normal_covector, outward_directed_normal_vector, interior_rest_mass_density, @@ -330,6 +377,7 @@ def _return_cons_vars( return { "tilde_d": tilde_d( face_mesh_velocity, + reflect_both, outward_directed_normal_covector, outward_directed_normal_vector, interior_rest_mass_density, @@ -345,6 +393,7 @@ def _return_cons_vars( ), "tilde_ye": tilde_ye( face_mesh_velocity, + reflect_both, outward_directed_normal_covector, outward_directed_normal_vector, interior_rest_mass_density, @@ -360,6 +409,7 @@ def _return_cons_vars( ), "tilde_tau": tilde_tau( face_mesh_velocity, + reflect_both, outward_directed_normal_covector, outward_directed_normal_vector, interior_rest_mass_density, @@ -375,6 +425,7 @@ def _return_cons_vars( ), "tilde_s": tilde_s( face_mesh_velocity, + reflect_both, outward_directed_normal_covector, outward_directed_normal_vector, interior_rest_mass_density, @@ -390,6 +441,7 @@ def _return_cons_vars( ), "tilde_b": tilde_b( face_mesh_velocity, + reflect_both, outward_directed_normal_covector, outward_directed_normal_vector, interior_rest_mass_density, @@ -405,6 +457,7 @@ def _return_cons_vars( ), "tilde_phi": tilde_phi( face_mesh_velocity, + reflect_both, outward_directed_normal_covector, outward_directed_normal_vector, interior_rest_mass_density, @@ -423,6 +476,7 @@ def _return_cons_vars( def flux_tilde_d( face_mesh_velocity, + reflect_both, outward_directed_normal_covector, outward_directed_normal_vector, interior_rest_mass_density, @@ -440,24 +494,27 @@ def flux_tilde_d( 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, + reflect_both, outward_directed_normal_covector, outward_directed_normal_vector, interior_rest_mass_density, interior_electron_fraction, interior_specific_internal_energy, interior_spatial_velocity, - exterior_magnetic_field, + interior_magnetic_field, interior_lorentz_factor, interior_pressure, shift, @@ -486,6 +543,7 @@ def flux_tilde_d( def flux_tilde_ye( face_mesh_velocity, + reflect_both, outward_directed_normal_covector, outward_directed_normal_vector, interior_rest_mass_density, @@ -503,17 +561,20 @@ def flux_tilde_ye( 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, + reflect_both, outward_directed_normal_covector, outward_directed_normal_vector, interior_rest_mass_density, @@ -549,6 +610,7 @@ def flux_tilde_ye( def flux_tilde_tau( face_mesh_velocity, + reflect_both, outward_directed_normal_covector, outward_directed_normal_vector, interior_rest_mass_density, @@ -566,17 +628,20 @@ def flux_tilde_tau( 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, + reflect_both, outward_directed_normal_covector, outward_directed_normal_vector, interior_rest_mass_density, @@ -612,6 +677,7 @@ def flux_tilde_tau( def flux_tilde_s( face_mesh_velocity, + reflect_both, outward_directed_normal_covector, outward_directed_normal_vector, interior_rest_mass_density, @@ -629,17 +695,20 @@ def flux_tilde_s( 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, + reflect_both, outward_directed_normal_covector, outward_directed_normal_vector, interior_rest_mass_density, @@ -675,6 +744,7 @@ def flux_tilde_s( def flux_tilde_b( face_mesh_velocity, + reflect_both, outward_directed_normal_covector, outward_directed_normal_vector, interior_rest_mass_density, @@ -692,17 +762,20 @@ def flux_tilde_b( 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, + reflect_both, outward_directed_normal_covector, outward_directed_normal_vector, interior_rest_mass_density, @@ -738,6 +811,7 @@ def flux_tilde_b( def flux_tilde_phi( face_mesh_velocity, + reflect_both, outward_directed_normal_covector, outward_directed_normal_vector, interior_rest_mass_density, @@ -755,17 +829,20 @@ def flux_tilde_phi( 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, + reflect_both, outward_directed_normal_covector, outward_directed_normal_vector, interior_rest_mass_density, @@ -801,6 +878,7 @@ def flux_tilde_phi( def lapse( face_mesh_velocity, + reflect_both, outward_directed_normal_covector, outward_directed_normal_vector, interior_rest_mass_density, @@ -819,6 +897,7 @@ def lapse( def shift( face_mesh_velocity, + reflect_both, outward_directed_normal_covector, outward_directed_normal_vector, interior_rest_mass_density, @@ -837,6 +916,7 @@ def shift( def inv_spatial_metric( face_mesh_velocity, + reflect_both, outward_directed_normal_covector, outward_directed_normal_vector, interior_rest_mass_density, 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..3bd430d98ae2c 100644 --- a/tests/Unit/Evolution/Systems/GrMhd/ValenciaDivClean/BoundaryConditions/Test_Reflective.cpp +++ b/tests/Unit/Evolution/Systems/GrMhd/ValenciaDivClean/BoundaryConditions/Test_Reflective.cpp @@ -24,10 +24,8 @@ #include "Utilities/TaggedTuple.hpp" namespace helpers = ::TestHelpers::evolution::dg; - -SPECTRE_TEST_CASE("Unit.GrMhd.BoundaryConditions.Reflective", - "[Unit][GrMhd]") { - pypp::SetupLocalPythonEnvironment local_python_env{""}; +namespace { +void test_stuffs(bool reflect_both) { MAKE_GENERATOR(gen); const auto face_mesh_index = Index<2>{2}; @@ -44,6 +42,9 @@ SPECTRE_TEST_CASE("Unit.GrMhd.BoundaryConditions.Reflective", gr::Tags::SqrtDetSpatialMetric>>( spatial_metric, sqrt_det_spatial_metric); + // for factory string + std::string reflect_both_str = (reflect_both) ? "true" : "false"; + helpers::test_boundary_condition_with_python< grmhd::ValenciaDivClean::BoundaryConditions::Reflective, grmhd::ValenciaDivClean::BoundaryConditions::BoundaryCondition, @@ -94,6 +95,18 @@ 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, + "Reflective:\n" + " ReflectBoth: " + + reflect_both_str + "\n", + face_mesh_index, box_with_gridless_tags, reflect_both, 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(true); +} diff --git a/tests/Unit/Helpers/Evolution/DiscontinuousGalerkin/BoundaryConditions.hpp b/tests/Unit/Helpers/Evolution/DiscontinuousGalerkin/BoundaryConditions.hpp index be4338b54ccc0..6f8f0a4ed60b8 100644 --- a/tests/Unit/Helpers/Evolution/DiscontinuousGalerkin/BoundaryConditions.hpp +++ b/tests/Unit/Helpers/Evolution/DiscontinuousGalerkin/BoundaryConditions.hpp @@ -180,6 +180,7 @@ void test_boundary_condition_with_python_impl( const BoundaryCondition& boundary_condition, const Index& face_points, const db::DataBox& box_of_volume_data, + std::optional member_bool, const tuples::TaggedTuple...>& ranges, const bool use_moving_mesh, tmpl::list /*meta*/, tmpl::list /*meta*/, @@ -414,7 +415,7 @@ void test_boundary_condition_with_python_impl( Variables exterior_face_fields{ number_of_points_on_face}; auto apply_bc = [&boundary_condition, &box_of_volume_data, epsilon, - &exterior_face_fields, &face_mesh_velocity, + member_bool, &exterior_face_fields, &face_mesh_velocity, &interior_normal_covector, &python_boundary_condition_functions, &python_module]( const auto&... interior_face_and_volume_args) { @@ -439,11 +440,19 @@ void test_boundary_condition_with_python_impl( const std::string& python_error_msg_function = get_python_error_message_function( python_boundary_condition_functions); - const auto python_error_message = - call_for_error_message( - python_module, python_error_msg_function, face_mesh_velocity, - interior_normal_covector, interior_face_and_volume_args..., - db::get(box_of_volume_data)...); + std::optional python_error_message{}; + if (member_bool.has_value()) { + python_error_message = call_for_error_message( + python_module, python_error_msg_function, face_mesh_velocity, + member_bool, interior_normal_covector, + interior_face_and_volume_args..., + db::get(box_of_volume_data)...); + } else { + python_error_message = call_for_error_message( + python_module, python_error_msg_function, face_mesh_velocity, + interior_normal_covector, interior_face_and_volume_args..., + db::get(box_of_volume_data)...); + } CAPTURE(python_error_msg_function); CAPTURE(python_error_message.value_or("")); CAPTURE(error_msg.value_or("")); @@ -470,12 +479,22 @@ void test_boundary_condition_with_python_impl( CAPTURE(pretty_type::short_name()); typename BoundaryCorrectionTag::type python_result{}; try { - python_result = pypp::call( - python_module, python_tag_function, face_mesh_velocity, - interior_normal_covector, interior_face_and_volume_args..., - db::get( - box_of_volume_data)...); + if (member_bool.has_value()) { + python_result = pypp::call( + python_module, python_tag_function, face_mesh_velocity, + member_bool, interior_normal_covector, + interior_face_and_volume_args..., + db::get( + box_of_volume_data)...); + } else { + python_result = pypp::call( + python_module, python_tag_function, face_mesh_velocity, + interior_normal_covector, interior_face_and_volume_args..., + db::get( + box_of_volume_data)...); + } } catch (const std::exception& e) { INFO("Failed python call with '" << e.what() << "'"); // Use REQUIRE(false) to print all the CAPTURE variables @@ -623,8 +642,88 @@ void test_boundary_condition_with_python( System, ConversionClassList, BoundaryCorrection>( generator, python_module, python_boundary_condition_functions, dynamic_cast(*boundary_condition), - face_points, box_of_volume_data, ranges, use_moving_mesh, - variables_tags{}, package_data_input_tags{}, + face_points, box_of_volume_data, std::nullopt, ranges, + use_moving_mesh, variables_tags{}, package_data_input_tags{}, + boundary_condition_dg_gridless_tags{}, + ExtraTagsForPythonFromDataBox{}, epsilon); + // Now serialize and deserialize and test again + INFO("Test boundary condition after serialization."); + const auto deserialized_bc = + serialize_and_deserialize(boundary_condition); + detail::test_boundary_condition_with_python_impl< + System, ConversionClassList, BoundaryCorrection>( + generator, python_module, python_boundary_condition_functions, + dynamic_cast(*deserialized_bc), + face_points, box_of_volume_data, std::nullopt, ranges, + use_moving_mesh, variables_tags{}, package_data_input_tags{}, + boundary_condition_dg_gridless_tags{}, + ExtraTagsForPythonFromDataBox{}, epsilon); + } + }); +} +// Variant that takes extra boolean argument. (For, Reflective BC) +template , + typename ExtraTagsForPythonFromDataBox = tmpl::list<>, + typename Metavariables = NoSuchType, size_t FaceDim, + typename DbTagsList, typename... RangeTags, + typename... PythonFunctionNameTags> +void test_boundary_condition_with_python( + const gsl::not_null generator, + const std::string& python_module, + const tuples::TaggedTuple& + python_boundary_condition_functions, + const std::string& factory_string, const Index& face_points, + const db::DataBox& box_of_volume_data, const bool member_bool, + const tuples::TaggedTuple...>& ranges, + const double epsilon = 1.0e-12) { + PUPable_reg(BoundaryCondition); + static_assert(std::is_final_v>, + "All boundary condition classes must be marked `final`."); + static_assert(tt::is_a_v); + using variables_tags = typename System::variables_tag::tags_list; + using flux_variables = typename System::flux_variables; + using fluxes_tags = + db::wrap_tags_in<::Tags::Flux, flux_variables, tmpl::size_t, + Frame::Inertial>; + const std::unique_ptr + boundary_condition = [&factory_string]() { + if constexpr (std::is_same_v) { + return TestHelpers::test_factory_creation( + factory_string); + } else { + return TestHelpers::test_creation< + std::unique_ptr, Metavariables>( + factory_string); + } + }(); + + REQUIRE_FALSE( + domain::BoundaryConditions::is_periodic(boundary_condition->get_clone())); + + tmpl::for_each( + [&boundary_condition, &box_of_volume_data, epsilon, member_bool, + &face_points, &generator, &python_boundary_condition_functions, + &python_module, &ranges](auto boundary_correction_v) { + using BoundaryCorrection = + tmpl::type_from; + using package_data_input_tags = tmpl::append< + variables_tags, fluxes_tags, + typename BoundaryCorrection::dg_package_data_temporary_tags, + typename ::evolution::dg::Actions::detail::get_primitive_vars< + System::has_primitive_and_conservative_vars>:: + template f>; + using boundary_condition_dg_gridless_tags = + typename BoundaryCondition::dg_gridless_tags; + for (const auto use_moving_mesh : {false, true}) { + detail::test_boundary_condition_with_python_impl< + System, ConversionClassList, BoundaryCorrection>( + generator, python_module, python_boundary_condition_functions, + dynamic_cast(*boundary_condition), + face_points, box_of_volume_data, member_bool, ranges, + use_moving_mesh, variables_tags{}, package_data_input_tags{}, boundary_condition_dg_gridless_tags{}, ExtraTagsForPythonFromDataBox{}, epsilon); // Now serialize and deserialize and test again @@ -635,8 +734,8 @@ void test_boundary_condition_with_python( System, ConversionClassList, BoundaryCorrection>( generator, python_module, python_boundary_condition_functions, dynamic_cast(*deserialized_bc), - face_points, box_of_volume_data, ranges, use_moving_mesh, - variables_tags{}, package_data_input_tags{}, + face_points, box_of_volume_data, member_bool, ranges, + use_moving_mesh, variables_tags{}, package_data_input_tags{}, boundary_condition_dg_gridless_tags{}, ExtraTagsForPythonFromDataBox{}, epsilon); }