Skip to content

Commit

Permalink
Generalize Poisson system to complex numbers
Browse files Browse the repository at this point in the history
  • Loading branch information
@nilsvu
nilsvu committed Aug 28, 2024
1 parent 028e94c commit 173f836
Show file tree
Hide file tree
Showing 4 changed files with 132 additions and 115 deletions.
128 changes: 66 additions & 62 deletions src/Elliptic/Systems/Poisson/Equations.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@

#include <cstddef>

#include "DataStructures/ComplexDataVector.hpp"
#include "DataStructures/DataBox/PrefixHelpers.hpp"
#include "DataStructures/DataVector.hpp"
#include "DataStructures/Tensor/EagerMath/DotProduct.hpp"
Expand All @@ -15,111 +16,114 @@

namespace Poisson {

template <size_t Dim>
template <typename DataType, size_t Dim>
void flat_cartesian_fluxes(
const gsl::not_null<tnsr::I<DataVector, Dim>*> flux_for_field,
const tnsr::i<DataVector, Dim>& field_gradient) {
const gsl::not_null<tnsr::I<DataType, Dim>*> flux_for_field,
const tnsr::i<DataType, Dim>& field_gradient) {
for (size_t d = 0; d < Dim; d++) {
flux_for_field->get(d) = field_gradient.get(d);
}
}

template <size_t Dim>
void curved_fluxes(
const gsl::not_null<tnsr::I<DataVector, Dim>*> flux_for_field,
const tnsr::II<DataVector, Dim>& inv_spatial_metric,
const tnsr::i<DataVector, Dim>& field_gradient) {
template <typename DataType, size_t Dim>
void curved_fluxes(const gsl::not_null<tnsr::I<DataType, Dim>*> flux_for_field,
const tnsr::II<DataVector, Dim>& inv_spatial_metric,
const tnsr::i<DataType, Dim>& field_gradient) {
raise_or_lower_index(flux_for_field, field_gradient, inv_spatial_metric);
}

template <size_t Dim>
void fluxes_on_face(gsl::not_null<tnsr::I<DataVector, Dim>*> flux_for_field,
template <typename DataType, size_t Dim>
void fluxes_on_face(gsl::not_null<tnsr::I<DataType, Dim>*> flux_for_field,
const tnsr::I<DataVector, Dim>& face_normal_vector,
const Scalar<DataVector>& field) {
*flux_for_field = face_normal_vector;
const Scalar<DataType>& field) {
std::copy(face_normal_vector.begin(), face_normal_vector.end(),
flux_for_field->begin());
for (size_t d = 0; d < Dim; d++) {
flux_for_field->get(d) *= get(field);
}
}

template <size_t Dim>
void add_curved_sources(
const gsl::not_null<Scalar<DataVector>*> source_for_field,
const tnsr::i<DataVector, Dim>& christoffel_contracted,
const tnsr::I<DataVector, Dim>& flux_for_field) {
template <typename DataType, size_t Dim>
void add_curved_sources(const gsl::not_null<Scalar<DataType>*> source_for_field,
const tnsr::i<DataVector, Dim>& christoffel_contracted,
const tnsr::I<DataType, Dim>& flux_for_field) {
get(*source_for_field) -=
get(dot_product(christoffel_contracted, flux_for_field));
}

template <size_t Dim>
void Fluxes<Dim, Geometry::FlatCartesian>::apply(
const gsl::not_null<tnsr::I<DataVector, Dim>*> flux_for_field,
const Scalar<DataVector>& /*field*/,
const tnsr::i<DataVector, Dim>& field_gradient) {
template <size_t Dim, typename DataType>
void Fluxes<Dim, Geometry::FlatCartesian, DataType>::apply(
const gsl::not_null<tnsr::I<DataType, Dim>*> flux_for_field,
const Scalar<DataType>& /*field*/,
const tnsr::i<DataType, Dim>& field_gradient) {
flat_cartesian_fluxes(flux_for_field, field_gradient);
}

template <size_t Dim>
void Fluxes<Dim, Geometry::FlatCartesian>::apply(
const gsl::not_null<tnsr::I<DataVector, Dim>*> flux_for_field,
template <size_t Dim, typename DataType>
void Fluxes<Dim, Geometry::FlatCartesian, DataType>::apply(
const gsl::not_null<tnsr::I<DataType, Dim>*> flux_for_field,
const tnsr::i<DataVector, Dim>& /*face_normal*/,
const tnsr::I<DataVector, Dim>& face_normal_vector,
const Scalar<DataVector>& field) {
const Scalar<DataType>& field) {
fluxes_on_face(flux_for_field, face_normal_vector, field);
}

template <size_t Dim>
void Fluxes<Dim, Geometry::Curved>::apply(
const gsl::not_null<tnsr::I<DataVector, Dim>*> flux_for_field,
template <size_t Dim, typename DataType>
void Fluxes<Dim, Geometry::Curved, DataType>::apply(
const gsl::not_null<tnsr::I<DataType, Dim>*> flux_for_field,
const tnsr::II<DataVector, Dim>& inv_spatial_metric,
const Scalar<DataVector>& /*field*/,
const tnsr::i<DataVector, Dim>& field_gradient) {
const Scalar<DataType>& /*field*/,
const tnsr::i<DataType, Dim>& field_gradient) {
curved_fluxes(flux_for_field, inv_spatial_metric, field_gradient);
}

template <size_t Dim>
void Fluxes<Dim, Geometry::Curved>::apply(
const gsl::not_null<tnsr::I<DataVector, Dim>*> flux_for_field,
template <size_t Dim, typename DataType>
void Fluxes<Dim, Geometry::Curved, DataType>::apply(
const gsl::not_null<tnsr::I<DataType, Dim>*> flux_for_field,
const tnsr::II<DataVector, Dim>& /*inv_spatial_metric*/,
const tnsr::i<DataVector, Dim>& /*face_normal*/,
const tnsr::I<DataVector, Dim>& face_normal_vector,
const Scalar<DataVector>& field) {
const Scalar<DataType>& field) {
fluxes_on_face(flux_for_field, face_normal_vector, field);
}

template <size_t Dim>
void Sources<Dim, Geometry::Curved>::apply(
const gsl::not_null<Scalar<DataVector>*> equation_for_field,
template <size_t Dim, typename DataType>
void Sources<Dim, Geometry::Curved, DataType>::apply(
const gsl::not_null<Scalar<DataType>*> equation_for_field,
const tnsr::i<DataVector, Dim>& christoffel_contracted,
const Scalar<DataVector>& /*field*/,
const tnsr::I<DataVector, Dim>& field_flux) {
const Scalar<DataType>& /*field*/,
const tnsr::I<DataType, Dim>& field_flux) {
add_curved_sources(equation_for_field, christoffel_contracted, field_flux);
}

} // namespace Poisson

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

#define INSTANTIATE(_, data) \
template void Poisson::flat_cartesian_fluxes<DIM(data)>( \
const gsl::not_null<tnsr::I<DataVector, DIM(data)>*>, \
const tnsr::i<DataVector, DIM(data)>&); \
template void Poisson::curved_fluxes<DIM(data)>( \
const gsl::not_null<tnsr::I<DataVector, DIM(data)>*>, \
const tnsr::II<DataVector, DIM(data)>&, \
const tnsr::i<DataVector, DIM(data)>&); \
template void Poisson::fluxes_on_face<DIM(data)>( \
const gsl::not_null<tnsr::I<DataVector, DIM(data)>*>, \
const tnsr::I<DataVector, DIM(data)>&, const Scalar<DataVector>&); \
template void Poisson::add_curved_sources<DIM(data)>( \
const gsl::not_null<Scalar<DataVector>*>, \
const tnsr::i<DataVector, DIM(data)>&, \
const tnsr::I<DataVector, DIM(data)>&); \
template class Poisson::Fluxes<DIM(data), Poisson::Geometry::FlatCartesian>; \
template class Poisson::Fluxes<DIM(data), Poisson::Geometry::Curved>; \
template class Poisson::Sources<DIM(data), Poisson::Geometry::Curved>;

GENERATE_INSTANTIATIONS(INSTANTIATE, (1, 2, 3))
#define DTYPE(data) BOOST_PP_TUPLE_ELEM(0, data)
#define DIM(data) BOOST_PP_TUPLE_ELEM(1, data)

#define INSTANTIATE(_, data) \
template void Poisson::flat_cartesian_fluxes( \
const gsl::not_null<tnsr::I<DTYPE(data), DIM(data)>*>, \
const tnsr::i<DTYPE(data), DIM(data)>&); \
template void Poisson::curved_fluxes( \
const gsl::not_null<tnsr::I<DTYPE(data), DIM(data)>*>, \
const tnsr::II<DataVector, DIM(data)>&, \
const tnsr::i<DTYPE(data), DIM(data)>&); \
template void Poisson::fluxes_on_face( \
const gsl::not_null<tnsr::I<DTYPE(data), DIM(data)>*>, \
const tnsr::I<DataVector, DIM(data)>&, const Scalar<DTYPE(data)>&); \
template void Poisson::add_curved_sources( \
const gsl::not_null<Scalar<DTYPE(data)>*>, \
const tnsr::i<DataVector, DIM(data)>&, \
const tnsr::I<DTYPE(data), DIM(data)>&); \
template class Poisson::Fluxes<DIM(data), Poisson::Geometry::FlatCartesian, \
DTYPE(data)>; \
template class Poisson::Fluxes<DIM(data), Poisson::Geometry::Curved, \
DTYPE(data)>; \
template class Poisson::Sources<DIM(data), Poisson::Geometry::Curved, \
DTYPE(data)>;

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

#undef INSTANTIATE
#undef DIM
68 changes: 35 additions & 33 deletions src/Elliptic/Systems/Poisson/Equations.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,9 +18,11 @@ namespace PUP {
class er;
} // namespace PUP
namespace Poisson {
template <size_t Dim, Geometry BackgroundGeometry>
template <size_t Dim, Geometry BackgroundGeometry,
typename DataType = DataVector>
struct Fluxes;
template <size_t Dim, Geometry BackgroundGeometry>
template <size_t Dim, Geometry BackgroundGeometry,
typename DataType = DataVector>
struct Sources;
} // namespace Poisson
/// \endcond
Expand All @@ -31,30 +33,30 @@ namespace Poisson {
* \brief Compute the fluxes \f$F^i=\partial_i u(x)\f$ for the Poisson
* equation on a flat spatial metric in Cartesian coordinates.
*/
template <size_t Dim>
template <typename DataType, size_t Dim>
void flat_cartesian_fluxes(
gsl::not_null<tnsr::I<DataVector, Dim>*> flux_for_field,
const tnsr::i<DataVector, Dim>& field_gradient);
gsl::not_null<tnsr::I<DataType, Dim>*> flux_for_field,
const tnsr::i<DataType, Dim>& field_gradient);

/*!
* \brief Compute the fluxes \f$F^i=\gamma^{ij}\partial_j u(x)\f$
* for the curved-space Poisson equation on a spatial metric \f$\gamma_{ij}\f$.
*/
template <size_t Dim>
void curved_fluxes(gsl::not_null<tnsr::I<DataVector, Dim>*> flux_for_field,
template <typename DataType, size_t Dim>
void curved_fluxes(gsl::not_null<tnsr::I<DataType, Dim>*> flux_for_field,
const tnsr::II<DataVector, Dim>& inv_spatial_metric,
const tnsr::i<DataVector, Dim>& field_gradient);
const tnsr::i<DataType, Dim>& field_gradient);

/*!
* \brief Compute the fluxes $F^i=\gamma^{ij} n_j u$ where $n_j$ is the
* `face_normal`.
*
* The `face_normal_vector` is $\gamma^{ij} n_j$.
*/
template <size_t Dim>
void fluxes_on_face(gsl::not_null<tnsr::I<DataVector, Dim>*> flux_for_field,
template <typename DataType, size_t Dim>
void fluxes_on_face(gsl::not_null<tnsr::I<DataType, Dim>*> flux_for_field,
const tnsr::I<DataVector, Dim>& face_normal_vector,
const Scalar<DataVector>& field);
const Scalar<DataType>& field);

/*!
* \brief Add the sources \f$S=-\Gamma^i_{ij}v^j\f$
Expand All @@ -63,31 +65,31 @@ void fluxes_on_face(gsl::not_null<tnsr::I<DataVector, Dim>*> flux_for_field,
* These sources arise from the non-principal part of the Laplacian on a
* non-Euclidean background.
*/
template <size_t Dim>
void add_curved_sources(gsl::not_null<Scalar<DataVector>*> source_for_field,
template <typename DataType, size_t Dim>
void add_curved_sources(gsl::not_null<Scalar<DataType>*> source_for_field,
const tnsr::i<DataVector, Dim>& christoffel_contracted,
const tnsr::I<DataVector, Dim>& flux_for_field);
const tnsr::I<DataType, Dim>& flux_for_field);

/*!
* \brief Compute the fluxes \f$F^i\f$ for the Poisson equation on a flat
* metric in Cartesian coordinates.
*
* \see Poisson::FirstOrderSystem
*/
template <size_t Dim>
struct Fluxes<Dim, Geometry::FlatCartesian> {
template <size_t Dim, typename DataType>
struct Fluxes<Dim, Geometry::FlatCartesian, DataType> {
using argument_tags = tmpl::list<>;
using volume_tags = tmpl::list<>;
using const_global_cache_tags = tmpl::list<>;
static constexpr bool is_trivial = true;
static constexpr bool is_discontinuous = false;
static void apply(gsl::not_null<tnsr::I<DataVector, Dim>*> flux_for_field,
const Scalar<DataVector>& field,
const tnsr::i<DataVector, Dim>& field_gradient);
static void apply(gsl::not_null<tnsr::I<DataVector, Dim>*> flux_for_field,
static void apply(gsl::not_null<tnsr::I<DataType, Dim>*> flux_for_field,
const Scalar<DataType>& field,
const tnsr::i<DataType, Dim>& field_gradient);
static void apply(gsl::not_null<tnsr::I<DataType, Dim>*> flux_for_field,
const tnsr::i<DataVector, Dim>& face_normal,
const tnsr::I<DataVector, Dim>& face_normal_vector,
const Scalar<DataVector>& field);
const Scalar<DataType>& field);
};

/*!
Expand All @@ -96,23 +98,23 @@ struct Fluxes<Dim, Geometry::FlatCartesian> {
*
* \see Poisson::FirstOrderSystem
*/
template <size_t Dim>
struct Fluxes<Dim, Geometry::Curved> {
template <size_t Dim, typename DataType>
struct Fluxes<Dim, Geometry::Curved, DataType> {
using argument_tags =
tmpl::list<gr::Tags::InverseSpatialMetric<DataVector, Dim>>;
using volume_tags = tmpl::list<>;
using const_global_cache_tags = tmpl::list<>;
static constexpr bool is_trivial = true;
static constexpr bool is_discontinuous = false;
static void apply(gsl::not_null<tnsr::I<DataVector, Dim>*> flux_for_field,
static void apply(gsl::not_null<tnsr::I<DataType, Dim>*> flux_for_field,
const tnsr::II<DataVector, Dim>& inv_spatial_metric,
const Scalar<DataVector>& field,
const tnsr::i<DataVector, Dim>& field_gradient);
static void apply(gsl::not_null<tnsr::I<DataVector, Dim>*> flux_for_field,
const Scalar<DataType>& field,
const tnsr::i<DataType, Dim>& field_gradient);
static void apply(gsl::not_null<tnsr::I<DataType, Dim>*> flux_for_field,
const tnsr::II<DataVector, Dim>& inv_spatial_metric,
const tnsr::i<DataVector, Dim>& face_normal,
const tnsr::I<DataVector, Dim>& face_normal_vector,
const Scalar<DataVector>& field);
const Scalar<DataType>& field);
};

/*!
Expand All @@ -121,14 +123,14 @@ struct Fluxes<Dim, Geometry::Curved> {
*
* \see Poisson::FirstOrderSystem
*/
template <size_t Dim>
struct Sources<Dim, Geometry::Curved> {
template <size_t Dim, typename DataType>
struct Sources<Dim, Geometry::Curved, DataType> {
using argument_tags = tmpl::list<
gr::Tags::SpatialChristoffelSecondKindContracted<DataVector, Dim>>;
static void apply(gsl::not_null<Scalar<DataVector>*> equation_for_field,
static void apply(gsl::not_null<Scalar<DataType>*> equation_for_field,
const tnsr::i<DataVector, Dim>& christoffel_contracted,
const Scalar<DataVector>& field,
const tnsr::I<DataVector, Dim>& field_flux);
const Scalar<DataType>& field,
const tnsr::I<DataType, Dim>& field_flux);
};

} // namespace Poisson
26 changes: 18 additions & 8 deletions src/Elliptic/Systems/Poisson/FirstOrderSystem.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -53,20 +53,30 @@ namespace Poisson {
* The fluxes and sources simplify significantly when the background metric is
* flat and we employ Cartesian coordinates so \f$\gamma_{ij} = \delta_{ij}\f$
* and \f$\Gamma^i_{jk} = 0\f$. Set the template parameter `BackgroundGeometry`
* to `Poisson::Geometry::FlatCartesian` to specialise the system for this case.
* to `Poisson::Geometry::FlatCartesian` to specialize the system for this case.
* Set it to `Poisson::Geometry::Curved` for the general case.
*
* ## Complex Poisson equation
*
* This system can also be used to solve the complex Poisson equation where
* $u(x)$ and $f(x)$ are complex-valued, by setting the `DataType` template
* parameter to `ComplexDataVector`. Note that the real and imaginary sectors of
* the equations decouple, so they are essentially two independent Poisson
* equations. This is useful for testing the elliptic solver with complex-valued
* fields, and also as building blocks for other Poisson-like systems of
* equations that have additional complex-valued source terms.
*/
template <size_t Dim, Geometry BackgroundGeometry>
template <size_t Dim, Geometry BackgroundGeometry,
typename DataType = DataVector>
struct FirstOrderSystem
: tt::ConformsTo<elliptic::protocols::FirstOrderSystem> {
static constexpr size_t volume_dim = Dim;

using primal_fields = tmpl::list<Tags::Field<DataVector>>;
using primal_fields = tmpl::list<Tags::Field<DataType>>;
// We just use the standard `Flux` prefix because the fluxes don't have
// symmetries and we don't need to give them a particular meaning.
using primal_fluxes =
tmpl::list<::Tags::Flux<Tags::Field<DataVector>, tmpl::size_t<Dim>,
Frame::Inertial>>;
using primal_fluxes = tmpl::list<
::Tags::Flux<Tags::Field<DataType>, tmpl::size_t<Dim>, Frame::Inertial>>;

using background_fields = tmpl::conditional_t<
BackgroundGeometry == Geometry::FlatCartesian, tmpl::list<>,
Expand All @@ -77,10 +87,10 @@ struct FirstOrderSystem
tmpl::conditional_t<BackgroundGeometry == Geometry::FlatCartesian, void,
gr::Tags::InverseSpatialMetric<DataVector, Dim>>;

using fluxes_computer = Fluxes<Dim, BackgroundGeometry>;
using fluxes_computer = Fluxes<Dim, BackgroundGeometry, DataType>;
using sources_computer =
tmpl::conditional_t<BackgroundGeometry == Geometry::FlatCartesian, void,
Sources<Dim, BackgroundGeometry>>;
Sources<Dim, BackgroundGeometry, DataType>>;

using boundary_conditions_base =
elliptic::BoundaryConditions::BoundaryCondition<Dim>;
Expand Down
Loading

0 comments on commit 173f836

Please sign in to comment.