Skip to content

Commit

Permalink
Allow for general magnetic field configs in MagnetizedTovStar
Browse files Browse the repository at this point in the history
  • Loading branch information
@nilsdeppe
nilsdeppe committed Oct 7, 2023
1 parent 3756857 commit dc6c932
Show file tree
Hide file tree
Showing 5 changed files with 178 additions and 129 deletions.
120 changes: 56 additions & 64 deletions src/PointwiseFunctions/AnalyticData/GrMhd/MagnetizedTovStar.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,27 +16,48 @@
#include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
#include "PointwiseFunctions/Hydro/EquationsOfState/Factory.hpp"
#include "PointwiseFunctions/Hydro/Tags.hpp"
#include "Utilities/Algorithm.hpp"
#include "Utilities/CloneUniquePtrs.hpp"
#include "Utilities/ConstantExpressions.hpp"
#include "Utilities/ContainerHelpers.hpp"
#include "Utilities/ErrorHandling/Assert.hpp"
#include "Utilities/GenerateInstantiations.hpp"
#include "Utilities/MakeWithValue.hpp"

namespace grmhd::AnalyticData {
MagnetizedTovStar::MagnetizedTovStar() = default;
MagnetizedTovStar::MagnetizedTovStar(MagnetizedTovStar&& /*rhs*/) = default;
MagnetizedTovStar& MagnetizedTovStar::operator=(MagnetizedTovStar&& /*rhs*/) =
default;
MagnetizedTovStar::~MagnetizedTovStar() = default;

MagnetizedTovStar::MagnetizedTovStar(
const double central_rest_mass_density,
std::unique_ptr<MagnetizedTovStar::equation_of_state_type>
equation_of_state,
const RelativisticEuler::Solutions::TovCoordinates coordinate_system,
const size_t pressure_exponent, const double cutoff_pressure_fraction,
const double vector_potential_amplitude)
std::vector<std::unique_ptr<
grmhd::AnalyticData::InitialMagneticFields::InitialMagneticField>>
magnetic_fields)
: tov_star(central_rest_mass_density, std::move(equation_of_state),
coordinate_system),
pressure_exponent_(pressure_exponent),
cutoff_pressure_(cutoff_pressure_fraction *
get(this->equation_of_state().pressure_from_density(
Scalar<double>{central_rest_mass_density}))),
vector_potential_amplitude_(vector_potential_amplitude) {}
magnetic_fields_(std::move(magnetic_fields)) {}

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_)) {}

MagnetizedTovStar& MagnetizedTovStar::operator=(const MagnetizedTovStar& rhs) {
if (this == &rhs) {
return *this;
}
static_cast<RelativisticEuler::Solutions::TovStar&>(*this) =
static_cast<const RelativisticEuler::Solutions::TovStar&>(rhs);
magnetic_fields_ = clone_unique_ptrs(rhs.magnetic_fields_);
return *this;
}

std::unique_ptr<evolution::initial_data::InitialData>
MagnetizedTovStar::get_clone() const {
Expand All @@ -47,9 +68,7 @@ MagnetizedTovStar::MagnetizedTovStar(CkMigrateMessage* msg) : tov_star(msg) {}

void MagnetizedTovStar::pup(PUP::er& p) {
tov_star::pup(p);
p | pressure_exponent_;
p | cutoff_pressure_;
p | vector_potential_amplitude_;
p | magnetic_fields_;
}

namespace magnetized_tov_detail {
Expand All @@ -58,68 +77,41 @@ void MagnetizedTovVariables<DataType, Region>::operator()(
const gsl::not_null<tnsr::I<DataType, 3>*> magnetic_field,
const gsl::not_null<Cache*> cache,
hydro::Tags::MagneticField<DataType, 3> /*meta*/) const {
const size_t num_pts = get_size(get<0>(coords));
const auto& pressure_profile =
get(cache->get_var(*this, hydro::Tags::Pressure<DataType>{}));
const auto& dr_pressure_profile = get(cache->get_var(
*this,
RelativisticEuler::Solutions::tov_detail::Tags::DrPressure<DataType>{}));
const auto& sqrt_det_spatial_metric =
get(cache->get_var(*this, gr::Tags::SqrtDetSpatialMetric<DataType>{}));
for (size_t i = 0; i < num_pts; ++i) {
const double pressure = get_element(pressure_profile, i);
if (LIKELY(get_element(radius, i) > 1.0e-16)) {
if (pressure < cutoff_pressure) {
get_element(get<0>(*magnetic_field), i) = 0.0;
get_element(get<1>(*magnetic_field), i) = 0.0;
get_element(get<2>(*magnetic_field), i) = 0.0;
continue;
}

const double x = get_element(get<0>(coords), i);
const double y = get_element(get<1>(coords), i);
const double z = get_element(get<2>(coords), i);
const double radius_i = get_element(radius, i);
const double dr_pressure = get_element(dr_pressure_profile, i);
const double pressure_term =
pow(pressure - cutoff_pressure, pressure_exponent);
const double deriv_pressure_term =
pressure_exponent *
pow(pressure - cutoff_pressure,
static_cast<int>(pressure_exponent) - 1) *
dr_pressure;

get_element(get<0>(*magnetic_field), i) =
-x * z / radius_i * deriv_pressure_term;

get_element(get<1>(*magnetic_field), i) =
-y * z / radius_i * deriv_pressure_term;

get_element(get<2>(*magnetic_field), i) =
2.0 * pressure_term +
(square(x) + square(y)) / radius_i * deriv_pressure_term;
} else {
get_element(get<0>(*magnetic_field), i) = 0.0;
get_element(get<1>(*magnetic_field), i) = 0.0;
get_element(get<2>(*magnetic_field), i) =
2.0 * pow(pressure - cutoff_pressure, pressure_exponent);
}
}
const Scalar<DataType>& sqrt_det_spatial_metric =
cache->get_var(*this, gr::Tags::SqrtDetSpatialMetric<DataType>{});
const Scalar<DataType>& pressure =
cache->get_var(*this, hydro::Tags::Pressure<DataType>{});
const auto& deriv_pressure =
cache->get_var(*this, ::Tags::deriv<hydro::Tags::Pressure<DataType>,
tmpl::size_t<3>, Frame::Inertial>{});
set_number_of_grid_points(magnetic_field, get_size(get<0>(coords)));
for (size_t i = 0; i < 3; ++i) {
magnetic_field->get(i) *=
vector_potential_amplitude / sqrt_det_spatial_metric;
magnetic_field->get(i) = 0.0;
}

for (const auto& magnetic_field_compute : magnetic_fields) {
magnetic_field_compute->variables(magnetic_field, coords, pressure,
sqrt_det_spatial_metric, deriv_pressure);
}
}
} // namespace magnetized_tov_detail

PUP::able::PUP_ID MagnetizedTovStar::my_PUP_ID = 0;

bool operator==(const MagnetizedTovStar& lhs, const MagnetizedTovStar& rhs) {
return static_cast<const typename MagnetizedTovStar::tov_star&>(lhs) ==
static_cast<const typename MagnetizedTovStar::tov_star&>(rhs) and
lhs.pressure_exponent_ == rhs.pressure_exponent_ and
lhs.cutoff_pressure_ == rhs.cutoff_pressure_ and
lhs.vector_potential_amplitude_ == rhs.vector_potential_amplitude_;
bool equal =
static_cast<const typename MagnetizedTovStar::tov_star&>(lhs) ==
static_cast<const typename MagnetizedTovStar::tov_star&>(rhs) and
lhs.magnetic_fields_.size() == rhs.magnetic_fields_.size();
if (not equal) {
return false;
}
for (size_t i = 0; i < lhs.magnetic_fields_.size(); ++i) {
if (not lhs.magnetic_fields_[i]->is_equal(*rhs.magnetic_fields_[i])) {
return false;
}
}
return true;
}

bool operator!=(const MagnetizedTovStar& lhs, const MagnetizedTovStar& rhs) {
Expand Down
72 changes: 27 additions & 45 deletions src/PointwiseFunctions/AnalyticData/GrMhd/MagnetizedTovStar.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,8 @@
#include "Options/String.hpp"
#include "PointwiseFunctions/AnalyticData/AnalyticData.hpp"
#include "PointwiseFunctions/AnalyticData/GrMhd/AnalyticData.hpp"
#include "PointwiseFunctions/AnalyticData/GrMhd/InitialMagneticFields/Factory.hpp"
#include "PointwiseFunctions/AnalyticData/GrMhd/InitialMagneticFields/InitialMagneticField.hpp"
#include "PointwiseFunctions/AnalyticSolutions/RelativisticEuler/TovStar.hpp"
#include "PointwiseFunctions/Hydro/EquationsOfState/Factory.hpp"
#include "PointwiseFunctions/Hydro/TagsDeclarations.hpp"
Expand Down Expand Up @@ -42,20 +44,19 @@ struct MagnetizedTovVariables
using Base::radial_solution;
using Base::radius;

size_t pressure_exponent;
double cutoff_pressure;
double vector_potential_amplitude;
const std::vector<std::unique_ptr<
grmhd::AnalyticData::InitialMagneticFields::InitialMagneticField>>&
magnetic_fields;

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,
size_t local_pressure_exponent, double local_cutoff_pressure,
double local_vector_potential_amplitude)
const std::vector<std::unique_ptr<
grmhd::AnalyticData::InitialMagneticFields::InitialMagneticField>>&
mag_fields)
: Base(local_x, local_radius, local_radial_solution, local_eos),
pressure_exponent(local_pressure_exponent),
cutoff_pressure(local_cutoff_pressure),
vector_potential_amplitude(local_vector_potential_amplitude) {}
magnetic_fields(mag_fields) {}

void operator()(
gsl::not_null<tnsr::I<DataType, 3>*> magnetic_field,
Expand Down Expand Up @@ -194,32 +195,14 @@ class MagnetizedTovStar : public virtual evolution::initial_data::InitialData,
using tov_star = RelativisticEuler::Solutions::TovStar;

public:
struct PressureExponent {
using type = size_t;
struct MagneticFields {
using type = std::vector<std::unique_ptr<
grmhd::AnalyticData::InitialMagneticFields::InitialMagneticField>>;
static constexpr Options::String help = {
"The exponent n_s controlling the smoothness of the field"};
"Magnetic fields to superpose on the TOV solution."};
};

struct VectorPotentialAmplitude {
using type = double;
static constexpr Options::String help = {
"The amplitude A_b of the phi-component of the vector potential. This "
"controls the magnetic field strength."};
static type lower_bound() { return 0.0; }
};

struct CutoffPressureFraction {
using type = double;
static constexpr Options::String help = {
"The fraction of the central pressure below which there is no magnetic "
"field. p_cut = Fraction * p_max."};
static type lower_bound() { return 0.0; }
static type upper_bound() { return 1.0; }
};

using options =
tmpl::push_back<tov_star::options, PressureExponent,
CutoffPressureFraction, VectorPotentialAmplitude>;
using options = tmpl::push_back<tov_star::options, MagneticFields>;

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

Expand All @@ -228,20 +211,21 @@ class MagnetizedTovStar : public virtual evolution::initial_data::InitialData,
template <typename DataType>
using tags = typename tov_star::template tags<DataType>;

MagnetizedTovStar() = default;
MagnetizedTovStar(const MagnetizedTovStar& /*rhs*/) = default;
MagnetizedTovStar& operator=(const MagnetizedTovStar& /*rhs*/) = default;
MagnetizedTovStar(MagnetizedTovStar&& /*rhs*/) = default;
MagnetizedTovStar& operator=(MagnetizedTovStar&& /*rhs*/) = default;
~MagnetizedTovStar() override = default;
MagnetizedTovStar();
MagnetizedTovStar(const MagnetizedTovStar& rhs);
MagnetizedTovStar& operator=(const MagnetizedTovStar& rhs);
MagnetizedTovStar(MagnetizedTovStar&& /*rhs*/);
MagnetizedTovStar& operator=(MagnetizedTovStar&& /*rhs*/);
~MagnetizedTovStar() override;

MagnetizedTovStar(
double central_rest_mass_density,
std::unique_ptr<EquationsOfState::EquationOfState<true, 1>>
equation_of_state,
RelativisticEuler::Solutions::TovCoordinates coordinate_system,
size_t pressure_exponent, double cutoff_pressure_fraction,
double vector_potential_amplitude);
std::vector<std::unique_ptr<
grmhd::AnalyticData::InitialMagneticFields::InitialMagneticField>>
magnetic_fields);

auto get_clone() const
-> std::unique_ptr<evolution::initial_data::InitialData> override;
Expand All @@ -260,8 +244,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...>{}, pressure_exponent_, cutoff_pressure_,
vector_potential_amplitude_);
x, tmpl::list<Tags...>{}, magnetic_fields_);
}

// NOLINTNEXTLINE(google-runtime-references)
Expand All @@ -272,10 +255,9 @@ class MagnetizedTovStar : public virtual evolution::initial_data::InitialData,
const MagnetizedTovStar& rhs);

protected:
size_t pressure_exponent_ = std::numeric_limits<size_t>::max();
double cutoff_pressure_ = std::numeric_limits<double>::signaling_NaN();
double vector_potential_amplitude_ =
std::numeric_limits<double>::signaling_NaN();
std::vector<std::unique_ptr<
grmhd::AnalyticData::InitialMagneticFields::InitialMagneticField>>
magnetic_fields_{};
};

bool operator!=(const MagnetizedTovStar& lhs, const MagnetizedTovStar& rhs);
Expand Down
10 changes: 7 additions & 3 deletions .../Evolution/Systems/GrMhd/GhValenciaDivClean/BoundaryConditions/Test_DirichletAnalytic.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,6 +36,7 @@
#include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
#include "PointwiseFunctions/Hydro/Tags.hpp"
#include "Utilities/Gsl.hpp"
#include "Utilities/MakeVector.hpp"
#include "Utilities/Serialization/RegisterDerivedClassesWithCharm.hpp"
#include "Utilities/TMPL.hpp"
#include "Utilities/TaggedTuple.hpp"
Expand Down Expand Up @@ -336,9 +337,12 @@ SPECTRE_TEST_CASE(
std::make_unique<EquationsOfState::PolytropicFluid<true>>(100.0,
2.0),
RelativisticEuler::Solutions::TovCoordinates::Schwarzschild,
2,
0.04,
2500.0};
make_vector<
std::unique_ptr<grmhd::AnalyticData::InitialMagneticFields::
InitialMagneticField>>(
std::make_unique<
grmhd::AnalyticData::InitialMagneticFields::Poloidal>(
2, 0.04, 2500.0, std::array{0.0, 0.0, 0.0}, 100.0))};
const auto serialized_and_deserialized_condition =
serialize_and_deserialize(
*dynamic_cast<grmhd::GhValenciaDivClean::BoundaryConditions::
Expand Down
10 changes: 7 additions & 3 deletions ...olution/Systems/GrMhd/GhValenciaDivClean/BoundaryConditions/Test_DirichletFreeOutflow.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -44,6 +44,7 @@
#include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
#include "PointwiseFunctions/Hydro/Tags.hpp"
#include "Utilities/Gsl.hpp"
#include "Utilities/MakeVector.hpp"
#include "Utilities/Serialization/RegisterDerivedClassesWithCharm.hpp"
#include "Utilities/TMPL.hpp"
#include "Utilities/TaggedTuple.hpp"
Expand Down Expand Up @@ -532,9 +533,12 @@ SPECTRE_TEST_CASE(
std::make_unique<EquationsOfState::PolytropicFluid<true>>(100.0,
2.0),
RelativisticEuler::Solutions::TovCoordinates::Schwarzschild,
2,
0.04,
2500.0};
make_vector<
std::unique_ptr<grmhd::AnalyticData::InitialMagneticFields::
InitialMagneticField>>(
std::make_unique<
grmhd::AnalyticData::InitialMagneticFields::Poloidal>(
2, 0.04, 2500.0, std::array{0.0, 0.0, 0.0}, 100.0))};
const auto serialized_and_deserialized_condition =
serialize_and_deserialize(
*dynamic_cast<grmhd::GhValenciaDivClean::BoundaryConditions::
Expand Down
Loading

0 comments on commit dc6c932

Please sign in to comment.