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>; } // 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 +#include +#include +#include + +#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 +Reflective::get_clone() const { + return std::make_unique(*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 Reflective::dg_ghost( + const gsl::not_null*> tilde_d, + const gsl::not_null*> tilde_ye, + const gsl::not_null*> tilde_tau, + const gsl::not_null*> tilde_s, + const gsl::not_null*> tilde_b, + const gsl::not_null*> tilde_phi, + + const gsl::not_null*> tilde_d_flux, + const gsl::not_null*> tilde_ye_flux, + const gsl::not_null*> + tilde_tau_flux, + const gsl::not_null*> tilde_s_flux, + const gsl::not_null*> tilde_b_flux, + const gsl::not_null*> + tilde_phi_flux, + + const gsl::not_null*> lapse, + const gsl::not_null*> shift, + const gsl::not_null*> + inv_spatial_metric, + + const std::optional>& + /*face_mesh_velocity*/, + const tnsr::i& + outward_directed_normal_covector, + const tnsr::I& + outward_directed_normal_vector, + + const Scalar& interior_rest_mass_density, + const Scalar& interior_electron_fraction, + const Scalar& interior_specific_internal_energy, + const tnsr::I& interior_spatial_velocity, + const tnsr::I& interior_magnetic_field, + const Scalar& interior_lorentz_factor, + const Scalar& interior_pressure, + + const tnsr::I& interior_shift, + const Scalar& interior_lapse, + 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 + // * 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, hydro::Tags::SpatialVelocity, + ::Tags::TempScalar<1>, hydro::Tags::MagneticField, + ::Tags::TempScalar<2>, gr::Tags::SpatialMetric, + gr::Tags::SqrtDetSpatialMetric>> + temp_buffer{number_of_grid_points}; + auto& normal_dot_interior_spatial_velocity = + get<::Tags::TempScalar<0>>(temp_buffer); + auto& exterior_spatial_velocity = + get>(temp_buffer); + auto& normal_dot_interior_magnetic_field = + get<::Tags::TempScalar<1>>(temp_buffer); + auto& exterior_magnetic_field = + get>(temp_buffer); + auto& exterior_divergence_cleaning_field = + get<::Tags::TempScalar<2>>(temp_buffer); + auto& interior_spatial_metric = + get>(temp_buffer); + auto& interior_sqrt_det_spatial_metric = + get>(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*> rest_mass_density, + const gsl::not_null*> electron_fraction, + const gsl::not_null*> temperature, + const gsl::not_null*> + lorentz_factor_times_spatial_velocity, + const gsl::not_null*> + magnetic_field, + const gsl::not_null*> divergence_cleaning_field, + + const gsl::not_null>>*> + cell_centered_ghost_fluxes, + + const Direction<3>& direction, + + // fd_interior_temporary_tags + const Mesh<3>& subcell_mesh, + const tnsr::I& interior_shift, + const Scalar& interior_lapse, + const tnsr::ii& interior_spatial_metric, + + // fd_interior_primitive_variables_tags + const Scalar& interior_rest_mass_density, + const Scalar& interior_electron_fraction, + const Scalar& interior_temperature, + const Scalar& interior_pressure, + const Scalar& interior_specific_internal_energy, + const Scalar& interior_lorentz_factor, + const tnsr::I& interior_spatial_velocity, + const tnsr::I& interior_magnetic_field, + + // fd_gridless_tags + const fd::Reconstructor& reconstructor) const { + Variables> + temp_vars{get(*rest_mass_density).size()}; + fd_ghost_impl(rest_mass_density, electron_fraction, temperature, + make_not_null(&get(temp_vars)), + make_not_null(&get(temp_vars)), + lorentz_factor_times_spatial_velocity, + make_not_null(&get(temp_vars)), + make_not_null(&get(temp_vars)), magnetic_field, + divergence_cleaning_field, + make_not_null(&get(temp_vars)), + make_not_null(&get(temp_vars)), + make_not_null(&get(temp_vars)), + make_not_null(&get(temp_vars)), + make_not_null(&get(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(temp_vars)), + make_not_null(&get(temp_vars)), + make_not_null(&get(temp_vars)), + make_not_null(&get>(temp_vars)), + make_not_null(&get>(temp_vars)), + make_not_null(&get(temp_vars)), + + // Note: Only the spatial velocity changes. + *rest_mass_density, *electron_fraction, + get(temp_vars), get(temp_vars), + get(temp_vars), get(temp_vars), + *magnetic_field, + + get(temp_vars), get(temp_vars), + *divergence_cleaning_field); + + ComputeFluxes::apply( + make_not_null( + &get>(cell_centered_ghost_fluxes->value())), + make_not_null( + &get>(cell_centered_ghost_fluxes->value())), + make_not_null( + &get>(cell_centered_ghost_fluxes->value())), + make_not_null( + &get>>(cell_centered_ghost_fluxes->value())), + make_not_null( + &get>>(cell_centered_ghost_fluxes->value())), + make_not_null( + &get>(cell_centered_ghost_fluxes->value())), + + get(temp_vars), get(temp_vars), + get(temp_vars), get>(temp_vars), + get>(temp_vars), get(temp_vars), + + get(temp_vars), get(temp_vars), + get(temp_vars), get(temp_vars), + get(temp_vars), get(temp_vars), + get(temp_vars), get(temp_vars), + *magnetic_field); + } +} + +void Reflective::fd_ghost_impl( + const gsl::not_null*> rest_mass_density, + const gsl::not_null*> electron_fraction, + const gsl::not_null*> temperature, + const gsl::not_null*> pressure, + const gsl::not_null*> specific_internal_energy, + const gsl::not_null*> + lorentz_factor_times_spatial_velocity, + const gsl::not_null*> + spatial_velocity, + const gsl::not_null*> lorentz_factor, + const gsl::not_null*> + magnetic_field, + const gsl::not_null*> divergence_cleaning_field, + const gsl::not_null*> + spatial_metric, + const gsl::not_null*> + inv_spatial_metric, + const gsl::not_null*> sqrt_det_spatial_metric, + const gsl::not_null*> lapse, + const gsl::not_null*> shift, + + const Direction<3>& direction, + + // fd_interior_temporary_tags + const Mesh<3>& subcell_mesh, + + // fd_interior_primitive_variables_tags + const Scalar& interior_rest_mass_density, + const Scalar& interior_electron_fraction, + const Scalar& interior_temperature, + const Scalar& interior_pressure, + const Scalar& interior_specific_internal_energy, + const Scalar& interior_lorentz_factor, + const tnsr::I& interior_spatial_velocity, + const tnsr::I& interior_magnetic_field, + const tnsr::ii& interior_spatial_metric, + const Scalar& interior_lapse, + const tnsr::I& 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; + + // Create a single large DV to reduce the number of Variables allocations + using fluxes_tags = + tmpl::list; + const size_t buffer_size_for_fluxes = + need_tags_for_fluxes + ? Variables::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 outermost_prim_vars{ + buffer_for_vars.data(), num_face_pts * buffer_size_per_grid_pts}; + Variables 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(outermost_prim_vars) = + get_boundary_val(interior_rest_mass_density); + get(outermost_prim_vars) = + get_boundary_val(interior_electron_fraction); + get(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(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) = + -1.0 * 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); + } + } + } + // 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(outermost_prim_vars).get(i) = + get(get_boundary_val(interior_lorentz_factor)) * + min(normal_spatial_velocity_at_boundary, + -normal_spatial_velocity_at_boundary); + 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); + 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)) * + get_boundary_val(interior_spatial_velocity).get(i); + get(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(ghost_prim_vars); + *electron_fraction = get(ghost_prim_vars); + *temperature = get(ghost_prim_vars); + *lorentz_factor_times_spatial_velocity = + get(ghost_prim_vars); + *magnetic_field = get(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 outermost_fluxes_vars{ + std::next(ghost_prim_vars.data(), + static_cast(ghost_prim_vars.size())), + num_face_pts * buffer_size_for_fluxes}; + Variables ghost_fluxes_vars{ + std::next(outermost_fluxes_vars.data(), + static_cast(outermost_fluxes_vars.size())), + num_face_pts * buffer_size_for_fluxes * ghost_zone_size}; + + get(outermost_fluxes_vars) = get_boundary_val(interior_pressure); + get(outermost_fluxes_vars) = + get_boundary_val(interior_specific_internal_energy); + get(outermost_fluxes_vars) = + get_boundary_val(interior_spatial_metric); + get(outermost_fluxes_vars) = get_boundary_val(interior_lapse); + get(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(ghost_fluxes_vars); + *specific_internal_energy = get(ghost_fluxes_vars); + *spatial_metric = get(ghost_fluxes_vars); + *lapse = get(ghost_fluxes_vars); + *shift = get(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( + 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 +#include +#include +#include + +#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 +class Direction; +namespace gsl { +template +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; + using ElectronFraction = hydro::Tags::ElectronFraction; + using Temperature = hydro::Tags::Temperature; + using Pressure = hydro::Tags::Pressure; + using LorentzFactorTimesSpatialVelocity = + hydro::Tags::LorentzFactorTimesSpatialVelocity; + using MagneticField = hydro::Tags::MagneticField; + using DivergenceCleaningField = + hydro::Tags::DivergenceCleaningField; + using SpecificInternalEnergy = + hydro::Tags::SpecificInternalEnergy; + using SpatialVelocity = hydro::Tags::SpatialVelocity; + using LorentzFactor = hydro::Tags::LorentzFactor; + using SqrtDetSpatialMetric = gr::Tags::SqrtDetSpatialMetric; + using SpatialMetric = gr::Tags::SpatialMetric; + using InvSpatialMetric = gr::Tags::InverseSpatialMetric; + using Lapse = gr::Tags::Lapse; + using Shift = gr::Tags::Shift; + + using prim_tags_for_reconstruction = + tmpl::list; + + template + using Flux = ::Tags::Flux, 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; + 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::ElectronFraction, + hydro::Tags::SpecificInternalEnergy, + hydro::Tags::SpatialVelocity, + hydro::Tags::MagneticField, + hydro::Tags::LorentzFactor, + hydro::Tags::Pressure>; + using dg_interior_temporary_tags = tmpl::list; + using dg_gridless_tags = tmpl::list<>; + + std::optional dg_ghost( + gsl::not_null*> tilde_d, + gsl::not_null*> tilde_ye, + gsl::not_null*> tilde_tau, + gsl::not_null*> tilde_s, + gsl::not_null*> tilde_b, + gsl::not_null*> tilde_phi, + + gsl::not_null*> tilde_d_flux, + gsl::not_null*> tilde_ye_flux, + gsl::not_null*> tilde_tau_flux, + gsl::not_null*> tilde_s_flux, + gsl::not_null*> tilde_b_flux, + gsl::not_null*> tilde_phi_flux, + + gsl::not_null*> lapse, + gsl::not_null*> shift, + gsl::not_null*> + inv_spatial_metric, + + const std::optional>& + /*face_mesh_velocity*/, + const tnsr::i& + outward_directed_normal_covector, + const tnsr::I& + outward_directed_normal_vector, + + const Scalar& interior_rest_mass_density, + const Scalar& interior_electron_fraction, + const Scalar& interior_specific_internal_energy, + const tnsr::I& interior_spatial_velocity, + const tnsr::I& interior_magnetic_field, + const Scalar& interior_lorentz_factor, + const Scalar& interior_pressure, + + const tnsr::I& interior_shift, + const Scalar& interior_lapse, + const tnsr::II& + interior_inv_spatial_metric) const; + + using fd_interior_evolved_variables_tags = tmpl::list<>; + using fd_interior_temporary_tags = + tmpl::list, Shift, Lapse, + SpatialMetric>; + using fd_interior_primitive_variables_tags = + tmpl::list, + hydro::Tags::SpecificInternalEnergy, + hydro::Tags::LorentzFactor, + hydro::Tags::SpatialVelocity, MagneticField>; + + using fd_gridless_tags = tmpl::list; + + void fd_ghost( + gsl::not_null*> rest_mass_density, + gsl::not_null*> electron_fraction, + gsl::not_null*> temperature, + gsl::not_null*> + lorentz_factor_times_spatial_velocity, + gsl::not_null*> magnetic_field, + gsl::not_null*> divergence_cleaning_field, + + gsl::not_null>>*> + cell_centered_ghost_fluxes, + + const Direction<3>& direction, + + // interior temporary tags + const Mesh<3>& subcell_mesh, + const tnsr::I& interior_shift, + const Scalar& interior_lapse, + const tnsr::ii& interior_spatial_metric, + + // interior prim vars tags + const Scalar& interior_rest_mass_density, + const Scalar& interior_electron_fraction, + const Scalar& interior_temperature, + const Scalar& interior_pressure, + const Scalar& interior_specific_internal_energy, + const Scalar& interior_lorentz_factor, + const tnsr::I& interior_spatial_velocity, + const tnsr::I& 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*> rest_mass_density, + gsl::not_null*> electron_fraction, + gsl::not_null*> temperature, + gsl::not_null*> pressure, + gsl::not_null*> specific_internal_energy, + gsl::not_null*> + lorentz_factor_times_spatial_velocity, + gsl::not_null*> spatial_velocity, + gsl::not_null*> lorentz_factor, + gsl::not_null*> magnetic_field, + gsl::not_null*> divergence_cleaning_field, + gsl::not_null*> spatial_metric, + gsl::not_null*> + inv_spatial_metric, + gsl::not_null*> sqrt_det_spatial_metric, + gsl::not_null*> lapse, + gsl::not_null*> shift, + + const Direction<3>& direction, + + // fd_interior_temporary_tags + const Mesh<3>& subcell_mesh, + + // fd_interior_primitive_variables_tags + const Scalar& interior_rest_mass_density, + const Scalar& interior_electron_fraction, + const Scalar& interior_temperature, + const Scalar& interior_pressure, + const Scalar& interior_specific_internal_energy, + const Scalar& interior_lorentz_factor, + const tnsr::I& interior_spatial_velocity, + const tnsr::I& interior_magnetic_field, + const tnsr::ii& interior_spatial_metric, + const Scalar& interior_lapse, + const tnsr::I& 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 + +#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, + gr::Tags::SqrtDetSpatialMetric, + 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, + tmpl::list<>, tmpl::list>( + 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>, + helpers::Tags::PythonFunctionName< + grmhd::ValenciaDivClean::Tags::TildeB>, + helpers::Tags::PythonFunctionName< + grmhd::ValenciaDivClean::Tags::TildePhi>, + + helpers::Tags::PythonFunctionName< + ::Tags::Flux, Frame::Inertial>>, + helpers::Tags::PythonFunctionName< + ::Tags::Flux, Frame::Inertial>>, + helpers::Tags::PythonFunctionName< + ::Tags::Flux, Frame::Inertial>>, + helpers::Tags::PythonFunctionName<::Tags::Flux< + grmhd::ValenciaDivClean::Tags::TildeS, + tmpl::size_t<3>, Frame::Inertial>>, + helpers::Tags::PythonFunctionName<::Tags::Flux< + grmhd::ValenciaDivClean::Tags::TildeB, + tmpl::size_t<3>, Frame::Inertial>>, + helpers::Tags::PythonFunctionName< + ::Tags::Flux, Frame::Inertial>>, + + helpers::Tags::PythonFunctionName>, + helpers::Tags::PythonFunctionName>, + helpers::Tags::PythonFunctionName< + gr::Tags::InverseSpatialMetric>>{ + "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