Skip to content

Commit

Permalink
Add an option to reflect both outgoing and incoming (temp)
Browse files Browse the repository at this point in the history
  • Loading branch information
jyoo1042 committed Dec 12, 2023
1 parent 289ee7e commit 6c9ba04
Show file tree
Hide file tree
Showing 3 changed files with 111 additions and 53 deletions.
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand All @@ -24,20 +25,25 @@
#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<domain::BoundaryConditions::BoundaryCondition>
Reflective::get_clone() const {
return std::make_unique<Reflective>(*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;
Expand Down Expand Up @@ -81,8 +87,8 @@ std::optional<std::string> Reflective::dg_ghost(

const tnsr::I<DataVector, 3, Frame::Inertial>& interior_shift,
const Scalar<DataVector>& interior_lapse,
const tnsr::II<DataVector, 3, Frame::Inertial>&
interior_inv_spatial_metric) {
const tnsr::II<DataVector, 3, Frame::Inertial>& interior_inv_spatial_metric)
const {
const size_t number_of_grid_points = get(interior_rest_mass_density).size();

// temp buffer to store
Expand Down Expand Up @@ -139,35 +145,56 @@ std::optional<std::string> Reflective::dg_ghost(
outward_directed_normal_covector, interior_spatial_velocity);
dot_product(make_not_null(&normal_dot_interior_magnetic_field),
outward_directed_normal_covector, interior_magnetic_field);
for (size_t i = 0; i < number_of_grid_points; ++i) {
if (get(normal_dot_interior_spatial_velocity)[i] > 0.0) {

// Reflect both outward and inward normal component of
// spatial velocity and magnetic field
if (reflect_both_) {
for (size_t i = 0; i < number_of_grid_points; ++i) {
for (size_t spatial_index = 0; spatial_index < 3; ++spatial_index) {
exterior_spatial_velocity.get(spatial_index)[i] =
interior_spatial_velocity.get(spatial_index)[i] -
2.0 * get(normal_dot_interior_spatial_velocity)[i] *
outward_directed_normal_vector.get(spatial_index)[i];
}
} else {
for (size_t spatial_index = 0; spatial_index < 3; ++spatial_index) {
exterior_spatial_velocity.get(spatial_index)[i] =
interior_spatial_velocity.get(spatial_index)[i];
}
}
if (get(normal_dot_interior_magnetic_field)[i] > 0.0) {
for (size_t spatial_index = 0; spatial_index < 3; ++spatial_index) {

exterior_magnetic_field.get(spatial_index)[i] =
interior_magnetic_field.get(spatial_index)[i] -
2.0 * get(normal_dot_interior_magnetic_field)[i] *
outward_directed_normal_vector.get(spatial_index)[i];
}
} else {
for (size_t spatial_index = 0; spatial_index < 3; ++spatial_index) {
exterior_magnetic_field.get(spatial_index)[i] =
interior_magnetic_field.get(spatial_index)[i];
}
}
// Reflect only the outward normal component of
// spatial velocity and magnetic field
else {
for (size_t i = 0; i < number_of_grid_points; ++i) {
if (get(normal_dot_interior_spatial_velocity)[i] > 0.0) {
for (size_t spatial_index = 0; spatial_index < 3; ++spatial_index) {
exterior_spatial_velocity.get(spatial_index)[i] =
interior_spatial_velocity.get(spatial_index)[i] -
2.0 * get(normal_dot_interior_spatial_velocity)[i] *
outward_directed_normal_vector.get(spatial_index)[i];
}
} else {
for (size_t spatial_index = 0; spatial_index < 3; ++spatial_index) {
exterior_spatial_velocity.get(spatial_index)[i] =
interior_spatial_velocity.get(spatial_index)[i];
}
}
if (get(normal_dot_interior_magnetic_field)[i] > 0.0) {
for (size_t spatial_index = 0; spatial_index < 3; ++spatial_index) {
exterior_magnetic_field.get(spatial_index)[i] =
interior_magnetic_field.get(spatial_index)[i] -
2.0 * get(normal_dot_interior_magnetic_field)[i] *
outward_directed_normal_vector.get(spatial_index)[i];
}
} else {
for (size_t spatial_index = 0; spatial_index < 3; ++spatial_index) {
exterior_magnetic_field.get(spatial_index)[i] =
interior_magnetic_field.get(spatial_index)[i];
}
}
}
}

get(exterior_divergence_cleaning_field) = 0.0;

ConservativeFromPrimitive::apply(
Expand Down Expand Up @@ -222,7 +249,7 @@ void Reflective::fd_ghost(
const tnsr::I<DataVector, 3, Frame::Inertial>& interior_magnetic_field,

// fd_gridless_tags
const fd::Reconstructor& reconstructor) {
const fd::Reconstructor& reconstructor) const {
Variables<tmpl::push_back<typename System::variables_tag::tags_list,
SpatialVelocity, LorentzFactor, Pressure,
SpecificInternalEnergy, SqrtDetSpatialMetric,
Expand Down Expand Up @@ -340,7 +367,7 @@ void Reflective::fd_ghost_impl(
const Scalar<DataVector>& interior_lapse,
const tnsr::I<DataVector, 3, Frame::Inertial>& interior_shift,

const size_t ghost_zone_size, const bool need_tags_for_fluxes) {
const 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()};
Expand Down Expand Up @@ -387,38 +414,56 @@ void Reflective::fd_ghost_impl(
get_boundary_val(interior_temperature);

{
// Invert the outgoing components of spatial velocity and compute Wv^i
// and the outgoing components of magnetic field.
//
// Note : Here we require the grid to be Cartesian, therefore we will need
// to change the implementation below once subcell supports curved mesh.
const auto normal_spatial_velocity_at_boundary =
get_boundary_val(interior_spatial_velocity).get(dim_direction);
for (size_t i = 0; i < 3; ++i) {
if (i == dim_direction) {
if (direction.sign() > 0.0) {
// Reflect both the outgoing and ingoing components of spatial
// velocity and the magnetic field.
if (reflect_both_) {
for (size_t i = 0; i < 3; ++i) {
if (i == dim_direction) {
get<LorentzFactorTimesSpatialVelocity>(outermost_prim_vars).get(i) =
-1.0 * get(get_boundary_val(interior_lorentz_factor)) *
get_boundary_val(interior_spatial_velocity).get(i);
get<MagneticField>(outermost_prim_vars).get(i) =
-1.0 * get_boundary_val(interior_magnetic_field).get(i);
} else {
get<LorentzFactorTimesSpatialVelocity>(outermost_prim_vars).get(i) =
get(get_boundary_val(interior_lorentz_factor)) *
min(normal_spatial_velocity_at_boundary,
normal_spatial_velocity_at_boundary * -1.0);
get_boundary_val(interior_spatial_velocity).get(i);
get<MagneticField>(outermost_prim_vars).get(i) =
min(get_boundary_val(interior_magnetic_field).get(i),
-1.0 * get_boundary_val(interior_magnetic_field).get(i));
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<LorentzFactorTimesSpatialVelocity>(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<MagneticField>(outermost_prim_vars).get(i) =
min(get_boundary_val(interior_magnetic_field).get(i),
-1.0 * get_boundary_val(interior_magnetic_field).get(i));
} else {
get<LorentzFactorTimesSpatialVelocity>(outermost_prim_vars).get(i) =
get(get_boundary_val(interior_lorentz_factor)) *
max(normal_spatial_velocity_at_boundary,
normal_spatial_velocity_at_boundary * -1.0);
get<MagneticField>(outermost_prim_vars).get(i) =
max(get_boundary_val(interior_magnetic_field).get(i),
-1.0 * get_boundary_val(interior_magnetic_field).get(i));
}
} else {
get<LorentzFactorTimesSpatialVelocity>(outermost_prim_vars).get(i) =
get(get_boundary_val(interior_lorentz_factor)) *
max(normal_spatial_velocity_at_boundary,
normal_spatial_velocity_at_boundary * -1.0);
get_boundary_val(interior_spatial_velocity).get(i);
get<MagneticField>(outermost_prim_vars).get(i) =
max(get_boundary_val(interior_magnetic_field).get(i),
-1.0 * get_boundary_val(interior_magnetic_field).get(i));
get_boundary_val(interior_magnetic_field).get(i);
}
} else {
get<LorentzFactorTimesSpatialVelocity>(outermost_prim_vars).get(i) =
get(get_boundary_val(interior_lorentz_factor)) *
get_boundary_val(interior_spatial_velocity).get(i);
get<MagneticField>(outermost_prim_vars).get(i) =
get_boundary_val(interior_magnetic_field).get(i);
}
}
}
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -61,6 +61,8 @@ namespace grmhd::ValenciaDivClean::BoundaryConditions {
*/
class Reflective final : public BoundaryCondition {
private:
bool reflect_both_{false};

using RestMassDensity = hydro::Tags::RestMassDensity<DataVector>;
using ElectronFraction = hydro::Tags::ElectronFraction<DataVector>;
using Temperature = hydro::Tags::Temperature<DataVector>;
Expand Down Expand Up @@ -89,7 +91,13 @@ class Reflective final : public BoundaryCondition {
using Flux = ::Tags::Flux<T, tmpl::size_t<3>, 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<ReflectBoth>;
static constexpr Options::String help{
"Reflective boundary conditions, inverting the sign "
"of outgoing normal component of spatial velocity "
Expand All @@ -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(
Expand All @@ -127,7 +137,7 @@ class Reflective final : public BoundaryCondition {
using dg_interior_temporary_tags = tmpl::list<Shift, Lapse, InvSpatialMetric>;
using dg_gridless_tags = tmpl::list<>;

static std::optional<std::string> dg_ghost(
std::optional<std::string> dg_ghost(
gsl::not_null<Scalar<DataVector>*> tilde_d,
gsl::not_null<Scalar<DataVector>*> tilde_ye,
gsl::not_null<Scalar<DataVector>*> tilde_tau,
Expand Down Expand Up @@ -165,7 +175,7 @@ class Reflective final : public BoundaryCondition {
const tnsr::I<DataVector, 3, Frame::Inertial>& interior_shift,
const Scalar<DataVector>& interior_lapse,
const tnsr::II<DataVector, 3, Frame::Inertial>&
interior_inv_spatial_metric);
interior_inv_spatial_metric) const;

using fd_interior_evolved_variables_tags = tmpl::list<>;
using fd_interior_temporary_tags =
Expand All @@ -180,7 +190,7 @@ class Reflective final : public BoundaryCondition {

using fd_gridless_tags = tmpl::list<fd::Tags::Reconstructor>;

static void fd_ghost(
void fd_ghost(
gsl::not_null<Scalar<DataVector>*> rest_mass_density,
gsl::not_null<Scalar<DataVector>*> electron_fraction,
gsl::not_null<Scalar<DataVector>*> temperature,
Expand Down Expand Up @@ -212,10 +222,10 @@ class Reflective final : public BoundaryCondition {
const tnsr::I<DataVector, 3, Frame::Inertial>& 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<Scalar<DataVector>*> rest_mass_density,
gsl::not_null<Scalar<DataVector>*> electron_fraction,
gsl::not_null<Scalar<DataVector>*> temperature,
Expand Down Expand Up @@ -252,6 +262,6 @@ class Reflective final : public BoundaryCondition {
const Scalar<DataVector>& interior_lapse,
const tnsr::I<DataVector, 3, Frame::Inertial>& 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
Original file line number Diff line number Diff line change
Expand Up @@ -94,6 +94,9 @@ SPECTRE_TEST_CASE("Unit.GrMhd.BoundaryConditions.Reflective",
"tilde_phi", "flux_tilde_d", "flux_tilde_ye", "flux_tilde_tau",
"flux_tilde_s", "flux_tilde_b", "flux_tilde_phi", "lapse", "shift",
"inv_spatial_metric"},
"Reflective:\n", face_mesh_index, box_with_gridless_tags,
tuples::TaggedTuple<>{}, 1.0e-10);
"Reflective:\n"
" ReflectBoth: false\n",
face_mesh_index, box_with_gridless_tags, tuples::TaggedTuple<>{},
1.0e-10);
// FIXME: only testing ReflectBoth==false right now.
}

0 comments on commit 6c9ba04

Please sign in to comment.