Skip to content
Permalink
Browse files

Merge pull request #1461 from carmaza/ne_cleanup

Minor cleanup to evolve NewtonianEuler
  • Loading branch information...
kidder committed Apr 9, 2019
2 parents de4d5e4 + 0257a40 commit 589397ee1fa5fcc96e42a71bb6965c09818befc8
@@ -9,18 +9,23 @@
#include "Utilities/GenerateInstantiations.hpp"
#include "Utilities/Gsl.hpp"

// IWYU pragma: no_include <array>

// IWYU pragma: no_forward_declare Tensor

/// \cond
namespace NewtonianEuler {

template <size_t Dim>
void ConservativeFromPrimitive<Dim>::apply(
const gsl::not_null<Scalar<DataVector>*> mass_density_cons,
const gsl::not_null<tnsr::I<DataVector, Dim>*> momentum_density,
const gsl::not_null<Scalar<DataVector>*> energy_density,
const Scalar<DataVector>& mass_density,
const tnsr::I<DataVector, Dim>& velocity,
const Scalar<DataVector>& specific_internal_energy) noexcept {
get(*mass_density_cons) = get(mass_density);

for (size_t i = 0; i < Dim; ++i) {
momentum_density->get(i) = get(mass_density) * velocity.get(i);
}
@@ -11,6 +11,7 @@

// IWYU pragma: no_forward_declare NewtonianEuler::Tags::EnergyDensity
// IWYU pragma: no_forward_declare NewtonianEuler::Tags::MassDensity
// IWYU pragma: no_forward_declare NewtonianEuler::Tags::MassDensityCons
// IWYU pragma: no_forward_declare NewtonianEuler::Tags::MomentumDensity
// IWYU pragma: no_forward_declare NewtonianEuler::Tags::SpecificInternalEnergy
// IWYU pragma: no_forward_declare NewtonianEuler::Tags::Velocity
@@ -39,17 +40,20 @@ namespace NewtonianEuler {
* where \f$S^i\f$ is the momentum density, \f$e\f$ is the energy density,
* \f$\rho\f$ is the mass density, \f$v^i\f$ is the velocity, \f$v^2\f$ is its
* magnitude squared, and \f$\epsilon\f$ is the specific internal energy.
* In addition, this method returns the mass density as a conservative.
*/
template <size_t Dim>
struct ConservativeFromPrimitive {
using return_tags = tmpl::list<Tags::MomentumDensity<DataVector, Dim>,
using return_tags = tmpl::list<Tags::MassDensityCons<DataVector>,
Tags::MomentumDensity<DataVector, Dim>,
Tags::EnergyDensity<DataVector>>;

using argument_tags =
tmpl::list<Tags::MassDensity<DataVector>, Tags::Velocity<DataVector, Dim>,
Tags::SpecificInternalEnergy<DataVector>>;

static void apply(
gsl::not_null<Scalar<DataVector>*> mass_density_cons,
gsl::not_null<tnsr::I<DataVector, Dim>*> momentum_density,
gsl::not_null<Scalar<DataVector>*> energy_density,
const Scalar<DataVector>& mass_density,
@@ -17,7 +17,7 @@ namespace NewtonianEuler {

template <size_t Dim>
void ComputeFluxes<Dim>::apply(
const gsl::not_null<tnsr::I<DataVector, Dim>*> mass_density_flux,
const gsl::not_null<tnsr::I<DataVector, Dim>*> mass_density_cons_flux,
const gsl::not_null<tnsr::IJ<DataVector, Dim>*> momentum_density_flux,
const gsl::not_null<tnsr::I<DataVector, Dim>*> energy_density_flux,
const tnsr::I<DataVector, Dim>& momentum_density,
@@ -27,7 +27,7 @@ void ComputeFluxes<Dim>::apply(
const DataVector enthalpy_density = get(energy_density) + get(pressure);

for (size_t i = 0; i < Dim; ++i) {
mass_density_flux->get(i) = momentum_density.get(i);
mass_density_cons_flux->get(i) = momentum_density.get(i);
for (size_t j = 0; j < Dim; ++j) {
momentum_density_flux->get(i, j) =
momentum_density.get(i) * velocity.get(j);
@@ -10,6 +10,12 @@
#include "Evolution/Systems/NewtonianEuler/TagsDeclarations.hpp" // IWYU pragma: keep
#include "Utilities/TMPL.hpp"

// IWYU pragma: no_forward_declare NewtonianEuler::Tags::EnergyDensity
// IWYU pragma: no_forward_declare NewtonianEuler::Tags::MomentumDensity
// IWYU pragma: no_forward_declare NewtonianEuler::Tags::MassDensityCons
// IWYU pragma: no_forward_declare NewtonianEuler::Tags::Pressure
// IWYU pragma: no_forward_declare NewtonianEuler::Tags::Velocity
// IWYU pragma: no_forward_declare Tags::Flux
// IWYU pragma: no_forward_declare Tensor

/// \cond
@@ -48,8 +54,8 @@ namespace NewtonianEuler {
template <size_t Dim>
struct ComputeFluxes {
using return_tags =
tmpl::list<::Tags::Flux<Tags::MassDensity<DataVector>, tmpl::size_t<Dim>,
Frame::Inertial>,
tmpl::list<::Tags::Flux<Tags::MassDensityCons<DataVector>,
tmpl::size_t<Dim>, Frame::Inertial>,
::Tags::Flux<Tags::MomentumDensity<DataVector, Dim>,
tmpl::size_t<Dim>, Frame::Inertial>,
::Tags::Flux<Tags::EnergyDensity<DataVector>,
@@ -61,7 +67,7 @@ struct ComputeFluxes {
Tags::Velocity<DataVector, Dim>, Tags::Pressure<DataVector>>;

static void apply(
gsl::not_null<tnsr::I<DataVector, Dim>*> mass_density_flux,
gsl::not_null<tnsr::I<DataVector, Dim>*> mass_density_cons_flux,
gsl::not_null<tnsr::IJ<DataVector, Dim>*> momentum_density_flux,
gsl::not_null<tnsr::I<DataVector, Dim>*> energy_density_flux,
const tnsr::I<DataVector, Dim>& momentum_density,
@@ -6,42 +6,66 @@
#include "DataStructures/DataVector.hpp" // IWYU pragma: keep
#include "DataStructures/Tensor/EagerMath/DotProduct.hpp"
#include "DataStructures/Tensor/Tensor.hpp"
#include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp"
#include "PointwiseFunctions/Hydro/Tags.hpp" // IWYU pragma: keep
#include "Utilities/GenerateInstantiations.hpp"
#include "Utilities/Gsl.hpp"
#include "Utilities/Overloader.hpp"

// IWYU pragma: no_include <array>

// IWYU pragma: no_forward_declare EquationsOfState::EquationOfState
// IWYU pragma: no_forward_declare Tensor

/// \cond
namespace NewtonianEuler {

template <size_t Dim>
void PrimitiveFromConservative<Dim>::apply(
template <size_t Dim, size_t ThermodynamicDim>
void PrimitiveFromConservative<Dim, ThermodynamicDim>::apply(
const gsl::not_null<Scalar<DataVector>*> mass_density,
const gsl::not_null<tnsr::I<DataVector, Dim>*> velocity,
const gsl::not_null<Scalar<DataVector>*> specific_internal_energy,
const Scalar<DataVector>& mass_density,
const gsl::not_null<Scalar<DataVector>*> pressure,
const Scalar<DataVector>& mass_density_cons,
const tnsr::I<DataVector, Dim>& momentum_density,
const Scalar<DataVector>& energy_density) noexcept {
const DataVector one_over_mass_density = 1.0 / get(mass_density);
const Scalar<DataVector>& energy_density,
const EquationsOfState::EquationOfState<false, ThermodynamicDim>&
equation_of_state) noexcept {
get(*mass_density) = get(mass_density_cons);
const DataVector one_over_mass_density = 1.0 / get(mass_density_cons);

for (size_t i = 0; i < Dim; ++i) {
velocity->get(i) = momentum_density.get(i) * one_over_mass_density;
}

get(*specific_internal_energy) = one_over_mass_density * get(energy_density) -
0.5 * get(dot_product(*velocity, *velocity));

*pressure = make_overloader(
[&mass_density_cons](const EquationsOfState::EquationOfState<false, 1>&
the_equation_of_state) noexcept {
return the_equation_of_state.pressure_from_density(mass_density_cons);
},
[&mass_density_cons, &specific_internal_energy ](
const EquationsOfState::EquationOfState<false, 2>&
the_equation_of_state) noexcept {
return the_equation_of_state.pressure_from_density_and_energy(
mass_density_cons, *specific_internal_energy);
})(equation_of_state);
}

} // namespace NewtonianEuler

#define DIM(data) BOOST_PP_TUPLE_ELEM(0, data)
#define THERMO_DIM(data) BOOST_PP_TUPLE_ELEM(1, data)

#define INSTANTIATE(_, data) \
template struct NewtonianEuler::PrimitiveFromConservative<DIM(data)>;
#define INSTANTIATE(_, data) \
template struct NewtonianEuler::PrimitiveFromConservative<DIM(data), \
THERMO_DIM(data)>;

GENERATE_INSTANTIATIONS(INSTANTIATE, (1, 2, 3))
GENERATE_INSTANTIATIONS(INSTANTIATE, (1, 2, 3), (1, 2))

#undef DIM
#undef THERMO_DIM
#undef INSTANTIATE
/// \endcond
@@ -7,13 +7,19 @@

#include "DataStructures/Tensor/TypeAliases.hpp"
#include "Evolution/Systems/NewtonianEuler/TagsDeclarations.hpp" // IWYU pragma: keep
#include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp" // IWYU pragma: keep
#include "PointwiseFunctions/Hydro/Tags.hpp"
#include "Utilities/TMPL.hpp"

// IWYU pragma: no_forward_declare EquationsOfState::EquationOfState
// IWYU pragma: no_forward_declare NewtonianEuler::Tags::EnergyDensity
// IWYU pragma: no_forward_declare NewtonianEuler::Tags::MassDensity
// IWYU pragma: no_forward_declare NewtonianEuler::Tags::MassDensityCons
// IWYU pragma: no_forward_declare NewtonianEuler::Tags::MomentumDensity
// IWYU pragma: no_forward_declare NewtonianEuler::Tags::Pressure
// IWYU pragma: no_forward_declare NewtonianEuler::Tags::SpecificInternalEnergy
// IWYU pragma: no_forward_declare NewtonianEuler::Tags::Velocity
// IWYU pragma: no_forward_declare hydro::Tags::EquationOfStateBase

/// \cond
namespace gsl {
@@ -40,20 +46,30 @@ namespace NewtonianEuler {
* internal energy, \f$e\f$ is the energy density, \f$\rho\f$
* is the mass density, \f$S^i\f$ is the momentum density, and
* \f$S^2\f$ is the momentum density squared.
*
* This routine also returns the mass density as a primitive, and the pressure
* from a generic equation of state \f$p = p(\rho, \epsilon)\f$.
*/
template <size_t Dim>
template <size_t Dim, size_t ThermodynamicDim>
struct PrimitiveFromConservative {
using return_tags = tmpl::list<Tags::Velocity<DataVector, Dim>,
Tags::SpecificInternalEnergy<DataVector>>;

using argument_tags = tmpl::list<Tags::MassDensity<DataVector>,
Tags::MomentumDensity<DataVector, Dim>,
Tags::EnergyDensity<DataVector>>;

static void apply(gsl::not_null<tnsr::I<DataVector, Dim>*> velocity,
gsl::not_null<Scalar<DataVector>*> specific_internal_energy,
const Scalar<DataVector>& mass_density,
const tnsr::I<DataVector, Dim>& momentum_density,
const Scalar<DataVector>& energy_density) noexcept;
using return_tags =
tmpl::list<Tags::MassDensity<DataVector>, Tags::Velocity<DataVector, Dim>,
Tags::SpecificInternalEnergy<DataVector>,
Tags::Pressure<DataVector>>;

using argument_tags = tmpl::list<
Tags::MassDensityCons<DataVector>, Tags::MomentumDensity<DataVector, Dim>,
Tags::EnergyDensity<DataVector>, hydro::Tags::EquationOfStateBase>;

static void apply(
gsl::not_null<Scalar<DataVector>*> mass_density,
gsl::not_null<tnsr::I<DataVector, Dim>*> velocity,
gsl::not_null<Scalar<DataVector>*> specific_internal_energy,
gsl::not_null<Scalar<DataVector>*> pressure,
const Scalar<DataVector>& mass_density_cons,
const tnsr::I<DataVector, Dim>& momentum_density,
const Scalar<DataVector>& energy_density,
const EquationsOfState::EquationOfState<false, ThermodynamicDim>&
equation_of_state) noexcept;
};
} // namespace NewtonianEuler
@@ -48,5 +48,17 @@ struct SoundSpeedSquaredCompute : SoundSpeedSquared<DataType>, db::ComputeTag {
tmpl::list<MassDensity<DataType>, SpecificInternalEnergy<DataType>,
hydro::Tags::EquationOfStateBase>;
};

/// Compute item for the sound speed \f$c_s\f$.
///
/// Can be retrieved using `NewtonianEuler::Tags::SoundSpeed`
template <typename DataType>
struct SoundSpeedCompute : SoundSpeed<DataType>, db::ComputeTag {
static Scalar<DataType> function(
const Scalar<DataType>& sound_speed_squared) noexcept {
return Scalar<DataType>{sqrt(get(sound_speed_squared))};
}
using argument_tags = tmpl::list<SoundSpeedSquared<DataType>>;
};
} // namespace Tags
} // namespace NewtonianEuler
@@ -28,6 +28,13 @@ struct MassDensity : db::SimpleTag {
static std::string name() noexcept { return "MassDensity"; }
};

/// The mass density of the fluid (as a conservative variable).
template <typename DataType>
struct MassDensityCons : db::SimpleTag {
using type = Scalar<DataType>;
static std::string name() noexcept { return "MassDensityCons"; }
};

/// The momentum density of the fluid.
template <typename DataType, size_t Dim, typename Fr>
struct MomentumDensity : db::SimpleTag {
@@ -14,6 +14,8 @@ template <size_t Dim>
struct CharacteristicSpeeds;
template <typename DataType>
struct MassDensity;
template <typename DataType>
struct MassDensityCons;
template <typename DataType, size_t Dim, typename Fr = Frame::Inertial>
struct MomentumDensity;
template <typename DataType>
@@ -27,6 +29,8 @@ struct Pressure;
template <typename DataType>
struct SoundSpeed;
template <typename DataType>
struct SoundSpeedCompute;
template <typename DataType>
struct SoundSpeedSquared;
template <typename DataType>
struct SoundSpeedSquaredCompute;
@@ -110,6 +110,16 @@ IsentropicVortex<Dim>::variables(
variables(tmpl::list<Tags::MassDensity<DataType>>{}, vars)));
}

template <size_t Dim>
template <typename DataType>
tuples::TaggedTuple<Tags::Pressure<DataType>> IsentropicVortex<Dim>::variables(
tmpl::list<Tags::Pressure<DataType>> /*meta*/,
const IntermediateVariables<DataType>& vars) const noexcept {
return equation_of_state_.pressure_from_density(
get<Tags::MassDensity<DataType>>(
variables(tmpl::list<Tags::MassDensity<DataType>>{}, vars)));
}

template <size_t Dim>
bool operator==(const IsentropicVortex<Dim>& lhs,
const IsentropicVortex<Dim>& rhs) noexcept {
@@ -151,7 +161,8 @@ GENERATE_INSTANTIATIONS(INSTANTIATE_CLASS, (2, 3))
const IntermediateVariables<DTYPE(data)>&) const noexcept;

GENERATE_INSTANTIATIONS(INSTANTIATE_SCALARS, (2, 3), (double, DataVector),
(Tags::MassDensity, Tags::SpecificInternalEnergy))
(Tags::MassDensity, Tags::SpecificInternalEnergy,
Tags::Pressure))

#define INSTANTIATE_VELOCITY(_, data) \
template tuples::TaggedTuple<TAG(data) < DTYPE(data), DIM(data), \
@@ -183,6 +183,11 @@ class IsentropicVortex {
auto variables(tmpl::list<Tags::SpecificInternalEnergy<DataType>> /*meta*/,
const IntermediateVariables<DataType>& vars) const noexcept
-> tuples::TaggedTuple<Tags::SpecificInternalEnergy<DataType>>;

template <typename DataType>
auto variables(tmpl::list<Tags::Pressure<DataType>> /*meta*/,
const IntermediateVariables<DataType>& vars) const noexcept
-> tuples::TaggedTuple<Tags::Pressure<DataType>>;
// @}

// Intermediate variables needed to compute the primitives
@@ -18,6 +18,10 @@ def characteristic_speeds(velocity, sound_speed, normal):


# Functions for testing ConservativeFromPrimitive.cpp
def mass_density_cons(mass_density, velocity, specific_internal_energy):
return mass_density


def momentum_density(mass_density, velocity, specific_internal_energy):
return mass_density * velocity

@@ -31,7 +35,8 @@ def energy_density(mass_density, velocity, specific_internal_energy):


# Functions for testing Fluxes.cpp
def mass_density_flux(momentum_density, energy_density, velocity, pressure):
def mass_density_cons_flux(momentum_density, energy_density,
velocity, pressure):
return momentum_density


@@ -49,13 +54,30 @@ def energy_density_flux(momentum_density, energy_density, velocity, pressure):


# Functions for testing PrimitiveFromConservative.cpp
def velocity(mass_density, momentum_density, energy_density):
return (momentum_density / mass_density)
def mass_density(mass_density_cons, momentum_density, energy_density):
return mass_density_cons


def velocity(mass_density_cons, momentum_density, energy_density):
return (momentum_density / mass_density_cons)


def specific_internal_energy(mass_density_cons, momentum_density,
energy_density):
veloc = velocity(mass_density_cons, momentum_density, energy_density)
return (energy_density / mass_density_cons - 0.5 * np.dot(veloc, veloc))


# Hard-coded for PolytropicFluid (representative ThermoDim = 1 case)
def pressure_1d(mass_density_cons, momentum_density, energy_density):
return 1.4 * np.power(mass_density_cons, 5.0 / 3.0)


def specific_internal_energy(mass_density, momentum_density, energy_density):
veloc = velocity(mass_density, momentum_density, energy_density)
return (energy_density / mass_density - 0.5 * np.dot(veloc, veloc))
# Hard-coded for IdealFluid (representative ThermoDim = 2 case)
def pressure_2d(mass_density_cons, momentum_density, energy_density):
return ((2.0 / 3.0) * mass_density_cons *
specific_internal_energy(mass_density_cons, momentum_density,
energy_density))


# End functions for testing PrimitiveFromConservative.cpp
@@ -16,7 +16,7 @@ template <size_t Dim>
void test_conservative_from_primitive(const DataVector& used_for_size) {
pypp::check_with_random_values<3>(
&NewtonianEuler::ConservativeFromPrimitive<Dim>::apply, "TestFunctions",
{"momentum_density", "energy_density"},
{"mass_density_cons", "momentum_density", "energy_density"},
{{{-1.0, 1.0}, {-2.0, 2.0}, {-3.0, 3.0}}}, used_for_size);
}

Oops, something went wrong.

0 comments on commit 589397e

Please sign in to comment.
You can’t perform that action at this time.