Skip to content

Commit

Permalink
Merge pull request #5562 from nilsdeppe/wcns5z_ghmhd
Browse files Browse the repository at this point in the history
Add Wcns5z for ghmhd, add Enable options to LimitLorentz & FixConservatives
  • Loading branch information
nilsdeppe authored Oct 13, 2023
2 parents 7a2ee66 + 30cfebd commit 71112b4
Show file tree
Hide file tree
Showing 24 changed files with 852 additions and 64 deletions.
Original file line number Diff line number Diff line change
Expand Up @@ -273,7 +273,7 @@ struct get_thermodynamic_dim;
template <typename InitialData>
struct get_thermodynamic_dim<InitialData, true> {
// Controls the thermodynamic dim used for numeric initial data.
static constexpr size_t value = 1;
static constexpr size_t value = 2;
};

template <typename InitialData>
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,7 @@ spectre_target_sources(
PositivityPreservingAdaptiveOrder.cpp
Reconstructor.cpp
RegisterDerivedWithCharm.cpp
Wcns5z.cpp
)

spectre_target_headers(
Expand All @@ -30,4 +31,5 @@ spectre_target_headers(
Reconstructor.hpp
RegisterDerivedWithCharm.hpp
Tag.hpp
Wcns5z.hpp
)
Original file line number Diff line number Diff line change
Expand Up @@ -6,3 +6,4 @@
#include "Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/MonotonisedCentral.hpp"
#include "Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/PositivityPreservingAdaptiveOrder.hpp"
#include "Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/Reconstructor.hpp"
#include "Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/Wcns5z.hpp"
Original file line number Diff line number Diff line change
Expand Up @@ -66,9 +66,11 @@ PositivityPreservingAdaptiveOrderPrim::PositivityPreservingAdaptiveOrderPrim(
PARSE_ERROR(context, "None is not an allowed low-order reconstructor.");
}
if (alpha_7.has_value()) {
PARSE_ERROR(context, "Alpha7 hasn't been tested.");
six_to_the_alpha_7_ = pow(6.0, alpha_7.value());
}
if (alpha_9.has_value()) {
PARSE_ERROR(context, "Alpha9 hasn't been tested.");
eight_to_the_alpha_9_ = pow(8.0, alpha_9.value());
}
set_function_pointers();
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,7 @@ namespace grmhd::GhValenciaDivClean::fd {
/// \cond
class MonotonisedCentralPrim;
class PositivityPreservingAdaptiveOrderPrim;
class Wcns5zPrim;
/// \endcond

/*!
Expand All @@ -38,8 +39,9 @@ class Reconstructor : public PUP::able {
WRAPPED_PUPable_abstract(Reconstructor); // NOLINT
/// \endcond

using creatable_classes = tmpl::list<MonotonisedCentralPrim,
PositivityPreservingAdaptiveOrderPrim>;
using creatable_classes =
tmpl::list<MonotonisedCentralPrim, PositivityPreservingAdaptiveOrderPrim,
Wcns5zPrim>;

virtual std::unique_ptr<Reconstructor> get_clone() const = 0;

Expand Down

Large diffs are not rendered by default.

Original file line number Diff line number Diff line change
@@ -0,0 +1,204 @@
// Distributed under the MIT License.
// See LICENSE.txt for details.

#pragma once

#include <array>
#include <cstddef>
#include <limits>
#include <memory>
#include <utility>

#include "DataStructures/FixedHashMap.hpp"
#include "DataStructures/Variables.hpp"
#include "DataStructures/VariablesTag.hpp"
#include "Domain/Structure/MaxNumberOfNeighbors.hpp"
#include "Domain/Tags.hpp"
#include "Evolution/DgSubcell/Tags/GhostDataForReconstruction.hpp"
#include "Evolution/DgSubcell/Tags/Mesh.hpp"
#include "Evolution/Systems/GeneralizedHarmonic/Tags.hpp"
#include "Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/Reconstructor.hpp"
#include "Evolution/Systems/GrMhd/GhValenciaDivClean/System.hpp"
#include "Evolution/Systems/GrMhd/ValenciaDivClean/Tags.hpp"
#include "NumericalAlgorithms/FiniteDifference/FallbackReconstructorType.hpp"
#include "Options/Auto.hpp"
#include "Options/Context.hpp"
#include "Options/String.hpp"
#include "PointwiseFunctions/Hydro/Tags.hpp"
#include "Utilities/TMPL.hpp"

/// \cond
class DataVector;
template <size_t Dim>
class Direction;
template <size_t Dim>
class Element;
template <size_t Dim>
class ElementId;
namespace EquationsOfState {
template <bool IsRelativistic, size_t ThermodynamicDim>
class EquationOfState;
} // namespace EquationsOfState
template <size_t Dim>
class Mesh;
namespace gsl {
template <typename T>
class not_null;
} // namespace gsl
namespace PUP {
class er;
} // namespace PUP
template <typename TagsList>
class Variables;
namespace evolution::dg::subcell {
class GhostData;
} // namespace evolution::dg::subcell
/// \endcond

namespace grmhd::GhValenciaDivClean::fd {
/*!
* \brief Fifth order weighted nonlinear compact scheme reconstruction using the
* Z oscillation indicator. See ::fd::reconstruction::wcns5z() for details.
*
*/
class Wcns5zPrim : public Reconstructor {
private:
using prims_to_reconstruct_tags =
tmpl::list<hydro::Tags::RestMassDensity<DataVector>,
hydro::Tags::ElectronFraction<DataVector>,
hydro::Tags::Temperature<DataVector>,
hydro::Tags::LorentzFactorTimesSpatialVelocity<DataVector, 3>,
hydro::Tags::MagneticField<DataVector, 3>,
hydro::Tags::DivergenceCleaningField<DataVector>>;

using FallbackReconstructorType =
::fd::reconstruction::FallbackReconstructorType;

public:
static constexpr size_t dim = 3;

struct NonlinearWeightExponent {
using type = size_t;
static constexpr Options::String help = {
"The exponent q to which the oscillation indicator term is raised"};
};
struct Epsilon {
using type = double;
static constexpr Options::String help = {
"The parameter added to the oscillation indicators to avoid division "
"by zero"};
};
struct FallbackReconstructor {
using type = FallbackReconstructorType;
static constexpr Options::String help = {
"A reconstruction scheme to fallback to adaptively. Finite difference "
"will switch to this reconstruction scheme if there are more extrema "
"in a FD stencil than a specified number. See also the option "
"'MaxNumberOfExtrema' below. Adaptive fallback is disabled if 'None'."};
};
struct MaxNumberOfExtrema {
using type = size_t;
static constexpr Options::String help = {
"The maximum allowed number of extrema in FD stencil for using Wcns5z "
"reconstruction before switching to a low-order reconstruction. If "
"FallbackReconstructor=None, this option is ignored"};
};

using options = tmpl::list<NonlinearWeightExponent, Epsilon,
FallbackReconstructor, MaxNumberOfExtrema>;

static constexpr Options::String help{
"WCNS 5Z reconstruction scheme using primitive variables."};

Wcns5zPrim() = default;
Wcns5zPrim(Wcns5zPrim&&) = default;
Wcns5zPrim& operator=(Wcns5zPrim&&) = default;
Wcns5zPrim(const Wcns5zPrim&) = default;
Wcns5zPrim& operator=(const Wcns5zPrim&) = default;
~Wcns5zPrim() override = default;

Wcns5zPrim(size_t nonlinear_weight_exponent, double epsilon,
FallbackReconstructorType fallback_reconstructor,
size_t max_number_of_extrema);

explicit Wcns5zPrim(CkMigrateMessage* msg);

WRAPPED_PUPable_decl_base_template(Reconstructor, Wcns5zPrim);

auto get_clone() const -> std::unique_ptr<Reconstructor> override;

static constexpr bool use_adaptive_order = false;

void pup(PUP::er& p) override;

size_t ghost_zone_size() const override { return 3; }

using reconstruction_argument_tags =
tmpl::list<::Tags::Variables<hydro::grmhd_tags<DataVector>>,
typename System::variables_tag,
hydro::Tags::EquationOfStateBase, domain::Tags::Element<dim>,
evolution::dg::subcell::Tags::GhostDataForReconstruction<dim>,
evolution::dg::subcell::Tags::Mesh<dim>>;

template <size_t ThermodynamicDim, typename TagsList>
void reconstruct(
gsl::not_null<std::array<Variables<TagsList>, dim>*> vars_on_lower_face,
gsl::not_null<std::array<Variables<TagsList>, dim>*> vars_on_upper_face,
const Variables<hydro::grmhd_tags<DataVector>>& volume_prims,
const Variables<typename System::variables_tag::type::tags_list>&
volume_spacetime_and_cons_vars,
const EquationsOfState::EquationOfState<true, ThermodynamicDim>& eos,
const Element<dim>& element,
const FixedHashMap<
maximum_number_of_neighbors(dim),
std::pair<Direction<dim>, ElementId<dim>>,
evolution::dg::subcell::GhostData,
boost::hash<std::pair<Direction<dim>, ElementId<dim>>>>& ghost_data,
const Mesh<dim>& subcell_mesh) const;

template <size_t ThermodynamicDim, typename TagsList>
void reconstruct_fd_neighbor(
gsl::not_null<Variables<TagsList>*> vars_on_face,
const Variables<hydro::grmhd_tags<DataVector>>& subcell_volume_prims,
const Variables<
grmhd::GhValenciaDivClean::Tags::spacetime_reconstruction_tags>&
subcell_volume_spacetime_metric,
const EquationsOfState::EquationOfState<true, ThermodynamicDim>& eos,
const Element<dim>& element,
const FixedHashMap<
maximum_number_of_neighbors(dim),
std::pair<Direction<dim>, ElementId<dim>>,
evolution::dg::subcell::GhostData,
boost::hash<std::pair<Direction<dim>, ElementId<dim>>>>& ghost_data,
const Mesh<dim>& subcell_mesh,
const Direction<dim>& direction_to_reconstruct) const;

private:
// NOLINTNEXTLINE(readability-redundant-declaration)
friend bool operator==(const Wcns5zPrim& lhs, const Wcns5zPrim& rhs);
friend bool operator!=(const Wcns5zPrim& lhs, const Wcns5zPrim& rhs);

size_t nonlinear_weight_exponent_ = 0;
double epsilon_ = std::numeric_limits<double>::signaling_NaN();
FallbackReconstructorType fallback_reconstructor_ =
FallbackReconstructorType::None;
size_t max_number_of_extrema_ = 0;

void (*reconstruct_)(gsl::not_null<std::array<gsl::span<double>, dim>*>,
gsl::not_null<std::array<gsl::span<double>, dim>*>,
const gsl::span<const double>&,
const DirectionMap<dim, gsl::span<const double>>&,
const Index<dim>&, size_t, double, size_t) = nullptr;
void (*reconstruct_lower_neighbor_)(gsl::not_null<DataVector*>,
const DataVector&, const DataVector&,
const Index<dim>&, const Index<dim>&,
const Direction<dim>&, const double&,
const size_t&) = nullptr;
void (*reconstruct_upper_neighbor_)(gsl::not_null<DataVector*>,
const DataVector&, const DataVector&,
const Index<dim>&, const Index<dim>&,
const Direction<dim>&, const double&,
const size_t&) = nullptr;
};

} // namespace grmhd::GhValenciaDivClean::fd
Original file line number Diff line number Diff line change
Expand Up @@ -104,7 +104,7 @@ FixConservatives::FixConservatives(
const double safety_factor_for_magnetic_field,
const double safety_factor_for_momentum_density,
const double safety_factor_for_momentum_density_cutoff_d,
const double safety_factor_for_momentum_density_slope,
const double safety_factor_for_momentum_density_slope, const bool enable,
const Options::Context& context)
: minimum_rest_mass_density_times_lorentz_factor_(
minimum_rest_mass_density_times_lorentz_factor),
Expand All @@ -119,7 +119,8 @@ FixConservatives::FixConservatives(
safety_factor_for_momentum_density_cutoff_d_(
safety_factor_for_momentum_density_cutoff_d),
safety_factor_for_momentum_density_slope_(
safety_factor_for_momentum_density_slope) {
safety_factor_for_momentum_density_slope),
enable_(enable) {
if (minimum_rest_mass_density_times_lorentz_factor_ >
rest_mass_density_times_lorentz_factor_cutoff_) {
PARSE_ERROR(context,
Expand Down Expand Up @@ -161,6 +162,7 @@ void FixConservatives::pup(PUP::er& p) {
p | one_minus_safety_factor_for_momentum_density_;
p | safety_factor_for_momentum_density_cutoff_d_;
p | safety_factor_for_momentum_density_slope_;
p | enable_;
}

// WARNING!
Expand All @@ -182,6 +184,9 @@ bool FixConservatives::operator()(
const tnsr::II<DataVector, 3, Frame::Inertial>& inv_spatial_metric,
const Scalar<DataVector>& sqrt_det_spatial_metric) const {
bool needed_fixing = false;
if (not enable_) {
return needed_fixing;
}
const size_t size = get<0>(tilde_b).size();
Variables<tmpl::list<::Tags::TempScalar<0>, ::Tags::TempScalar<1>,
::Tags::TempScalar<2>, ::Tags::TempScalar<3>,
Expand Down Expand Up @@ -395,7 +400,8 @@ bool operator==(const FixConservatives& lhs, const FixConservatives& rhs) {
lhs.safety_factor_for_momentum_density_cutoff_d_ ==
rhs.safety_factor_for_momentum_density_cutoff_d_ and
lhs.safety_factor_for_momentum_density_slope_ ==
rhs.safety_factor_for_momentum_density_slope_;
rhs.safety_factor_for_momentum_density_slope_ and
lhs.enable_ == rhs.enable_;
}

bool operator!=(const FixConservatives& lhs, const FixConservatives& rhs) {
Expand Down
30 changes: 19 additions & 11 deletions src/Evolution/Systems/GrMhd/ValenciaDivClean/FixConservatives.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -154,21 +154,28 @@ class FixConservatives {
"Slope of safety factor for momentum density bound below "
"SafetyFactorForSCutoffD, express as a function of log10(rho*W)."};
};
/// Whether or not the limiting is enabled
struct Enable {
using type = bool;
static constexpr Options::String help = {
"If true then the limiting is applied."};
};

using options = tmpl::list<MinimumValueOfD, CutoffD, MinimumValueOfYe,
CutoffYe, SafetyFactorForB, SafetyFactorForS,
SafetyFactorForSCutoffD, SafetyFactorForSSlope>;
using options =
tmpl::list<MinimumValueOfD, CutoffD, MinimumValueOfYe, CutoffYe,
SafetyFactorForB, SafetyFactorForS, SafetyFactorForSCutoffD,
SafetyFactorForSSlope, Enable>;
static constexpr Options::String help = {
"Variable fixing used in Foucart's thesis.\n"};

FixConservatives(const double minimum_rest_mass_density_times_lorentz_factor,
const double rest_mass_density_times_lorentz_factor_cutoff,
const double minimum_electron_fraction,
const double electron_fraction_cutoff,
const double safety_factor_for_magnetic_field,
const double safety_factor_for_momentum_density,
const double safety_factor_for_momentum_density_cutoff_d,
const double safety_factor_for_momentum_density_slope,
FixConservatives(double minimum_rest_mass_density_times_lorentz_factor,
double rest_mass_density_times_lorentz_factor_cutoff,
double minimum_electron_fraction,
double electron_fraction_cutoff,
double safety_factor_for_magnetic_field,
double safety_factor_for_momentum_density,
double safety_factor_for_momentum_density_cutoff_d,
double safety_factor_for_momentum_density_slope, bool enable,
const Options::Context& context = {});

FixConservatives() = default;
Expand Down Expand Up @@ -222,6 +229,7 @@ class FixConservatives {
std::numeric_limits<double>::signaling_NaN()};
double safety_factor_for_momentum_density_slope_{
std::numeric_limits<double>::signaling_NaN()};
bool enable_{true};
};

bool operator!=(const FixConservatives& lhs, const FixConservatives& rhs);
Expand Down
14 changes: 11 additions & 3 deletions src/Evolution/VariableFixing/LimitLorentzFactor.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -16,20 +16,27 @@

namespace VariableFixing {
LimitLorentzFactor::LimitLorentzFactor(const double max_density_cutoff,
const double lorentz_factor_cap)
const double lorentz_factor_cap,
const bool enable)
: max_density_cuttoff_(max_density_cutoff),
lorentz_factor_cap_(lorentz_factor_cap) {}
lorentz_factor_cap_(lorentz_factor_cap),
enable_(enable) {}

void LimitLorentzFactor::pup(PUP::er& p) {
p | max_density_cuttoff_;
p | lorentz_factor_cap_;
p | enable_;
}

void LimitLorentzFactor::operator()(
const gsl::not_null<Scalar<DataVector>*> lorentz_factor,
const gsl::not_null<tnsr::I<DataVector, 3, Frame::Inertial>*>
spatial_velocity,
const Scalar<DataVector>& rest_mass_density) const {
if (not enable_) {
return;
}

constexpr size_t dim = 3;
const size_t number_of_grid_points = get(rest_mass_density).size();
for (size_t s = 0; s < number_of_grid_points; ++s) {
Expand All @@ -49,7 +56,8 @@ void LimitLorentzFactor::operator()(

bool operator==(const LimitLorentzFactor& lhs, const LimitLorentzFactor& rhs) {
return lhs.max_density_cuttoff_ == rhs.max_density_cuttoff_ and
lhs.lorentz_factor_cap_ == rhs.lorentz_factor_cap_;
lhs.lorentz_factor_cap_ == rhs.lorentz_factor_cap_ and
lhs.enable_ == rhs.enable_;
}

bool operator!=(const LimitLorentzFactor& lhs, const LimitLorentzFactor& rhs) {
Expand Down
Loading

0 comments on commit 71112b4

Please sign in to comment.