Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Add perturbation option to MagnetizedTovStar #6016

Open
wants to merge 1 commit into
base: develop
Choose a base branch
from
+78 −27
Open

Add perturbation option to MagnetizedTovStar #6016

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
36 changes: 33 additions & 3 deletions src/PointwiseFunctions/AnalyticData/GrMhd/MagnetizedTovStar.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -38,16 +38,19 @@ MagnetizedTovStar::MagnetizedTovStar(
const RelativisticEuler::Solutions::TovCoordinates coordinate_system,
std::vector<std::unique_ptr<
grmhd::AnalyticData::InitialMagneticFields::InitialMagneticField>>
magnetic_fields)
magnetic_fields,
const double perturbation)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe a more specific name, like radial_velocity_perturbation?

: tov_star(central_rest_mass_density, std::move(equation_of_state),
coordinate_system),
magnetic_fields_(std::move(magnetic_fields)) {}
magnetic_fields_(std::move(magnetic_fields)),
perturbation_(perturbation) {}

MagnetizedTovStar::MagnetizedTovStar(const MagnetizedTovStar& rhs)
: evolution::initial_data::InitialData{rhs},
RelativisticEuler::Solutions::TovStar(
static_cast<const RelativisticEuler::Solutions::TovStar&>(rhs)),
magnetic_fields_(clone_unique_ptrs(rhs.magnetic_fields_)) {}
magnetic_fields_(clone_unique_ptrs(rhs.magnetic_fields_)),
perturbation_(rhs.perturbation_) {}

MagnetizedTovStar& MagnetizedTovStar::operator=(const MagnetizedTovStar& rhs) {
if (this == &rhs) {
Expand All @@ -56,6 +59,7 @@ MagnetizedTovStar& MagnetizedTovStar::operator=(const MagnetizedTovStar& rhs) {
static_cast<RelativisticEuler::Solutions::TovStar&>(*this) =
static_cast<const RelativisticEuler::Solutions::TovStar&>(rhs);
magnetic_fields_ = clone_unique_ptrs(rhs.magnetic_fields_);
perturbation_ = rhs.perturbation_;
return *this;
}

Expand All @@ -69,6 +73,7 @@ MagnetizedTovStar::MagnetizedTovStar(CkMigrateMessage* msg) : tov_star(msg) {}
void MagnetizedTovStar::pup(PUP::er& p) {
tov_star::pup(p);
p | magnetic_fields_;
p | perturbation_;
}

namespace magnetized_tov_detail {
Expand All @@ -94,6 +99,31 @@ void MagnetizedTovVariables<DataType, Region>::operator()(
sqrt_det_spatial_metric, deriv_pressure);
}
}

template <typename DataType, StarRegion Region>
void MagnetizedTovVariables<DataType, Region>::operator()(
const gsl::not_null<tnsr::I<DataType, 3>*> spatial_velocity,
[[maybe_unused]] const gsl::not_null<Cache*> cache,
hydro::Tags::SpatialVelocity<DataType, 3> /*meta*/) const {
if constexpr (Region == StarRegion::Center or
Region == StarRegion::Exterior) {
get<0>(*spatial_velocity) = 0.;
get<1>(*spatial_velocity) = 0.;
get<2>(*spatial_velocity) = 0.;
} else {
const auto& areal_radius =
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hm, I'm a bit worried about the normalization here. I think probably the number that gets passed in should be the norm of the spatial velocity. |v| =sqrt(v^iv^j\gamma_{ij}). This way users would know if the perturbation they're passing in is reasonable or not. In this case this is pretty straightforward, for Schwarzschild coordinares you just have to multiply the spatial velocity by sqrt(\gamma^{rr}) = sqrt(1-2M(r)/r). I think for isotropic coordinates what you have is already correct, because \gamma_rr = \psi^4= (areal_radius/coord_radius)^2, so just dividing by the areal radius is all you have to do to get the spatial velocity to be normalized correctly. You should be able to check this by starting a run, and the max value of the spatial velocity in the domain should be the value you input to the input file.

radial_solution.coordinate_system() ==
RelativisticEuler::Solutions::TovCoordinates::Isotropic
? get(cache->get_var(
*this,
RelativisticEuler::Solutions::tov_detail::Tags::ArealRadius<
DataType>{}))
: radius;
get<0>(*spatial_velocity) = perturbation * coords.get(0) / areal_radius;
get<1>(*spatial_velocity) = perturbation * coords.get(1) / areal_radius;
get<2>(*spatial_velocity) = perturbation * coords.get(2) / areal_radius;
}
}
Comment on lines +103 to +126
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Here's where the perturbation is actually applied.

} // namespace magnetized_tov_detail

PUP::able::PUP_ID MagnetizedTovStar::my_PUP_ID = 0;
Expand Down
27 changes: 22 additions & 5 deletions src/PointwiseFunctions/AnalyticData/GrMhd/MagnetizedTovStar.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -47,21 +47,28 @@ struct MagnetizedTovVariables
const std::vector<std::unique_ptr<
grmhd::AnalyticData::InitialMagneticFields::InitialMagneticField>>&
magnetic_fields;
const double perturbation;

MagnetizedTovVariables(
const tnsr::I<DataType, 3>& local_x, const DataType& local_radius,
const RelativisticEuler::Solutions::TovSolution& local_radial_solution,
const EquationsOfState::EquationOfState<true, 1>& local_eos,
const std::vector<std::unique_ptr<
grmhd::AnalyticData::InitialMagneticFields::InitialMagneticField>>&
mag_fields)
mag_fields,
const double perturb)
: Base(local_x, local_radius, local_radial_solution, local_eos),
magnetic_fields(mag_fields) {}
magnetic_fields(mag_fields),
perturbation(perturb) {}

void operator()(
gsl::not_null<tnsr::I<DataType, 3>*> magnetic_field,
gsl::not_null<Cache*> cache,
hydro::Tags::MagneticField<DataType, 3> /*meta*/) const override;
void operator()(
gsl::not_null<tnsr::I<DataType, 3>*> spatial_velocity,
gsl::not_null<Cache*> cache,
hydro::Tags::SpatialVelocity<DataType, 3> /*meta*/) const override;
};

} // namespace magnetized_tov_detail
Expand Down Expand Up @@ -128,7 +135,15 @@ class MagnetizedTovStar : public virtual evolution::initial_data::InitialData,
"Magnetic fields to superpose on the TOV solution."};
};

using options = tmpl::push_back<tov_star::options, MagneticFields>;
struct Perturbation {
using type = double;
static constexpr Options::String help = {
"Size of radially symmetric spatial velocity perturbation for driving "
"oscillations."};
};

using options =
tmpl::push_back<tov_star::options, MagneticFields, Perturbation>;

static constexpr Options::String help = {"Magnetized TOV star."};

Expand All @@ -151,7 +166,8 @@ class MagnetizedTovStar : public virtual evolution::initial_data::InitialData,
RelativisticEuler::Solutions::TovCoordinates coordinate_system,
std::vector<std::unique_ptr<
grmhd::AnalyticData::InitialMagneticFields::InitialMagneticField>>
magnetic_fields);
magnetic_fields,
double perturbation);

auto get_clone() const
-> std::unique_ptr<evolution::initial_data::InitialData> override;
Expand All @@ -170,7 +186,7 @@ class MagnetizedTovStar : public virtual evolution::initial_data::InitialData,
tuples::TaggedTuple<Tags...> variables(const tnsr::I<DataType, 3>& x,
tmpl::list<Tags...> /*meta*/) const {
return variables_impl<magnetized_tov_detail::MagnetizedTovVariables>(
x, tmpl::list<Tags...>{}, magnetic_fields_);
x, tmpl::list<Tags...>{}, magnetic_fields_, perturbation_);
}

// NOLINTNEXTLINE(google-runtime-references)
Expand All @@ -184,6 +200,7 @@ class MagnetizedTovStar : public virtual evolution::initial_data::InitialData,
std::vector<std::unique_ptr<
grmhd::AnalyticData::InitialMagneticFields::InitialMagneticField>>
magnetic_fields_{};
double perturbation_ = std::numeric_limits<double>::signaling_NaN();
};

bool operator!=(const MagnetizedTovStar& lhs, const MagnetizedTovStar& rhs);
Expand Down
2 changes: 1 addition & 1 deletion src/PointwiseFunctions/AnalyticSolutions/RelativisticEuler/TovStar.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -215,7 +215,7 @@ struct TovVariables {
void operator()(gsl::not_null<Scalar<DataType>*> lorentz_factor,
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this also needs to be modified, right now I'm pretty sure it's being used to compute prim2con here , because initialization is just using the primitive variables tag list here, which points to this. I would love to be wrong about this though, because we don't really need the Lorentz factor if we have the spatial velocity and the spatial metric.

gsl::not_null<Cache*> cache,
hydro::Tags::LorentzFactor<DataType> /*meta*/) const;
void operator()(gsl::not_null<tnsr::I<DataType, 3>*> spatial_velocity,
virtual void operator()(gsl::not_null<tnsr::I<DataType, 3>*> spatial_velocity,
gsl::not_null<Cache*> cache,
hydro::Tags::SpatialVelocity<DataType, 3> /*meta*/) const;
virtual void operator()(
Expand Down
5 changes: 3 additions & 2 deletions .../Evolution/Systems/GrMhd/GhValenciaDivClean/BoundaryConditions/Test_DirichletAnalytic.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -362,7 +362,8 @@ SPECTRE_TEST_CASE(
" CutoffPressure: 0.04\n"
" VectorPotentialAmplitude: 2500\n"
" Center: [0.0, 0.0, 0.0]\n"
" MaxDistanceFromCenter: 100.0\n")
" MaxDistanceFromCenter: 100.0\n"
" Perturbation: 0.0\n")
->get_clone();

const gh::Solutions::WrappedGr<grmhd::AnalyticData::MagnetizedTovStar>
Expand All @@ -376,7 +377,7 @@ SPECTRE_TEST_CASE(
InitialMagneticField>>(
std::make_unique<
grmhd::AnalyticData::InitialMagneticFields::Poloidal>(
2, 0.04, 2500.0, std::array{0.0, 0.0, 0.0}, 100.0))};
2, 0.04, 2500.0, std::array{0.0, 0.0, 0.0}, 100.0)), 0.0};
const auto serialized_and_deserialized_condition =
serialize_and_deserialize(
*dynamic_cast<grmhd::GhValenciaDivClean::BoundaryConditions::
Expand Down
5 changes: 3 additions & 2 deletions ...olution/Systems/GrMhd/GhValenciaDivClean/BoundaryConditions/Test_DirichletFreeOutflow.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -562,7 +562,8 @@ SPECTRE_TEST_CASE(
" CutoffPressure: 0.04\n"
" VectorPotentialAmplitude: 2500\n"
" Center: [0.0, 0.0, 0.0]\n"
" MaxDistanceFromCenter: 100.0\n")
" MaxDistanceFromCenter: 100.0\n"
" Perturbation: 0.0\n")
->get_clone();
const gh::Solutions::WrappedGr<grmhd::AnalyticData::MagnetizedTovStar>
analytic_solution_or_data{
Expand All @@ -575,7 +576,7 @@ SPECTRE_TEST_CASE(
InitialMagneticField>>(
std::make_unique<
grmhd::AnalyticData::InitialMagneticFields::Poloidal>(
2, 0.04, 2500.0, std::array{0.0, 0.0, 0.0}, 100.0))};
2, 0.04, 2500.0, std::array{0.0, 0.0, 0.0}, 100.0)), 0.0};
const auto serialized_and_deserialized_condition =
serialize_and_deserialize(
*dynamic_cast<grmhd::GhValenciaDivClean::BoundaryConditions::
Expand Down
30 changes: 16 additions & 14 deletions tests/Unit/PointwiseFunctions/AnalyticData/GrMhd/Test_MagnetizedTovStar.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -44,26 +44,27 @@ void test_equality() {
1.28e-3,
std::make_unique<EquationsOfState::PolytropicFluid<true>>(100.0, 2.0),
TovCoordinates::Schwarzschild,
{}};
{},
0.0};
const auto mag_tov = serialize_and_deserialize(mag_tov_original);
CHECK(
mag_tov ==
MagnetizedTovStar(
1.28e-3,
std::make_unique<EquationsOfState::PolytropicFluid<true>>(100.0, 2.0),
TovCoordinates::Schwarzschild, {}));
TovCoordinates::Schwarzschild, {}, 0.0));
CHECK(
mag_tov !=
MagnetizedTovStar(
2.28e-3,
std::make_unique<EquationsOfState::PolytropicFluid<true>>(100.0, 2.0),
TovCoordinates::Schwarzschild, {}));
TovCoordinates::Schwarzschild, {}, 0.0));
CHECK(
mag_tov !=
MagnetizedTovStar(
1.28e-3,
std::make_unique<EquationsOfState::PolytropicFluid<true>>(100.0, 2.0),
TovCoordinates::Isotropic, {}));
TovCoordinates::Isotropic, {}, 0.0));
CHECK(
mag_tov !=
MagnetizedTovStar(
Expand All @@ -75,7 +76,7 @@ void test_equality() {
InitialMagneticField>>(
std::make_unique<
grmhd::AnalyticData::InitialMagneticFields::Poloidal>(
2, 0.04, 2500.0, std::array{0.0, 0.0, 0.0}, 100.0))));
2, 0.04, 2500.0, std::array{0.0, 0.0, 0.0}, 100.0)), 0.0));
CHECK(
MagnetizedTovStar(
1.28e-3,
Expand All @@ -86,7 +87,7 @@ void test_equality() {
InitialMagneticField>>(
std::make_unique<
grmhd::AnalyticData::InitialMagneticFields::Poloidal>(
2, 0.04, 2500.0, std::array{0.0, 0.0, 0.0}, 100.0))) ==
2, 0.04, 2500.0, std::array{0.0, 0.0, 0.0}, 100.0)), 0.0) ==
MagnetizedTovStar(
1.28e-3,
std::make_unique<EquationsOfState::PolytropicFluid<true>>(100.0, 2.0),
Expand All @@ -96,7 +97,7 @@ void test_equality() {
InitialMagneticField>>(
std::make_unique<
grmhd::AnalyticData::InitialMagneticFields::Poloidal>(
2, 0.04, 2500.0, std::array{0.0, 0.0, 0.0}, 100.0))));
2, 0.04, 2500.0, std::array{0.0, 0.0, 0.0}, 100.0)), 0.0));
CHECK(
MagnetizedTovStar(
1.28e-3,
Expand All @@ -107,7 +108,7 @@ void test_equality() {
InitialMagneticField>>(
std::make_unique<
grmhd::AnalyticData::InitialMagneticFields::Toroidal>(
2, 0.04, 2500.0, std::array{0.0, 0.0, 0.0}, 100.0))) !=
2, 0.04, 2500.0, std::array{0.0, 0.0, 0.0}, 100.0)), 0.0) !=
MagnetizedTovStar(
1.28e-3,
std::make_unique<EquationsOfState::PolytropicFluid<true>>(100.0, 2.0),
Expand All @@ -117,7 +118,7 @@ void test_equality() {
InitialMagneticField>>(
std::make_unique<
grmhd::AnalyticData::InitialMagneticFields::Poloidal>(
2, 0.04, 2500.0, std::array{0.0, 0.0, 0.0}, 100.0))));
2, 0.04, 2500.0, std::array{0.0, 0.0, 0.0}, 100.0)), 0.0));
CHECK(
MagnetizedTovStar(
1.28e-3,
Expand All @@ -128,7 +129,7 @@ void test_equality() {
InitialMagneticField>>(
std::make_unique<
grmhd::AnalyticData::InitialMagneticFields::Poloidal>(
2, 0.04, 2500.0, std::array{0.0, 0.0, 0.0}, 100.0))) !=
2, 0.04, 2500.0, std::array{0.0, 0.0, 0.0}, 100.0)), 0.0) !=
MagnetizedTovStar(
1.28e-3,
std::make_unique<EquationsOfState::PolytropicFluid<true>>(100.0, 2.0),
Expand All @@ -141,7 +142,7 @@ void test_equality() {
2, 0.04, 2500.0, std::array{0.0, 0.0, 0.0}, 100.0),
std::make_unique<
grmhd::AnalyticData::InitialMagneticFields::Toroidal>(
2, 0.04, 2500.0, std::array{0.0, 0.0, 0.0}, 100.0))));
2, 0.04, 2500.0, std::array{0.0, 0.0, 0.0}, 100.0)), 0.0));
// Check order of magnetic fields doesn't matter.
const MagnetizedTovStar toroidal_plus_poloidal(
1.28e-3,
Expand All @@ -154,7 +155,7 @@ void test_equality() {
2, 1.0e-6, 2500.0, std::array{0.0, 0.0, 0.0}, 100.0),
std::make_unique<
grmhd::AnalyticData::InitialMagneticFields::Poloidal>(
2, 1.0e-6, 2500.0, std::array{0.0, 0.0, 0.0}, 100.0)));
2, 1.0e-6, 2500.0, std::array{0.0, 0.0, 0.0}, 100.0)), 0.0);
const MagnetizedTovStar poloidal_plus_toroidal(
1.28e-3,
std::make_unique<EquationsOfState::PolytropicFluid<true>>(100.0, 2.0),
Expand All @@ -166,7 +167,7 @@ void test_equality() {
2, 1.0e-6, 2500.0, std::array{0.0, 0.0, 0.0}, 100.0),
std::make_unique<
grmhd::AnalyticData::InitialMagneticFields::Toroidal>(
2, 1.0e-6, 2500.0, std::array{0.0, 0.0, 0.0}, 100.0)));
2, 1.0e-6, 2500.0, std::array{0.0, 0.0, 0.0}, 100.0)), 0.0);
const DataVector coord{2.0};
const tnsr::I<DataVector, 3, Frame::Inertial> coords{{coord, coord, coord}};
CHECK_ITERABLE_APPROX(
Expand Down Expand Up @@ -204,7 +205,8 @@ void test_magnetized_tov_star(const TovCoordinates coord_system) {
" VectorPotentialAmplitude: 2500\n"
" CutoffPressure: 6.5536e-06\n" // 0.04 * 100 * (1.28e-3)**2
" Center: [0.0, 0.0, 0.0]\n"
" MaxDistanceFromCenter: 100.0\n")
" MaxDistanceFromCenter: 100.0\n"
" Perturbation: 0.0\n")
->get_clone();
const auto deserialized_option_solution =
serialize_and_deserialize(option_solution);
Expand Down
Loading