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.

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

Pass fluxes to elliptic source computers #2553

Merged
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
Original file line number Diff line number Diff line change
Expand Up @@ -71,7 +71,7 @@ struct InitializeFirstOrderOperator {
vars_tag, PrimalVariables,
Copy link
Member

Choose a reason for hiding this comment

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

The commit message is incorrect. The systems don't need the fluxes to compute sources, but that having access to the fluxes is an optimization.

AuxiliaryVariables>;
using sources_compute_tag = elliptic::Tags::FirstOrderSourcesCompute<
SourcesComputer, vars_tag, PrimalVariables, AuxiliaryVariables>;
Dim, SourcesComputer, vars_tag, PrimalVariables, AuxiliaryVariables>;

using exterior_tags = tmpl::list<
// On exterior (ghost) boundary faces we compute the fluxes from the
Expand Down
14 changes: 8 additions & 6 deletions src/Elliptic/FirstOrderComputeTags.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -46,23 +46,25 @@ struct FirstOrderFluxesCompute
}
};

template <typename SourcesComputer, typename VariablesTag,
template <size_t Dim, typename SourcesComputer, typename VariablesTag,
typename PrimalVariables, typename AuxiliaryVariables>
struct FirstOrderSourcesCompute
: db::add_tag_prefix<::Tags::Source, VariablesTag>,
db::ComputeTag {
using base = db::add_tag_prefix<::Tags::Source, VariablesTag>;
using argument_tags =
tmpl::push_front<typename SourcesComputer::argument_tags, VariablesTag>;
tmpl::push_front<typename SourcesComputer::argument_tags, VariablesTag,
db::add_tag_prefix<::Tags::Flux, VariablesTag,
tmpl::size_t<Dim>, Frame::Inertial>>;
using volume_tags = get_volume_tags<SourcesComputer>;
using return_type = typename base::type;
template <typename... SourcesArgs>
template <typename Vars, typename Fluxes, typename... SourcesArgs>
static void function(const gsl::not_null<return_type*> sources,
const typename VariablesTag::type& vars,
const Vars& vars, const Fluxes& fluxes,
const SourcesArgs&... sources_args) noexcept {
*sources = return_type{vars.number_of_grid_points()};
elliptic::first_order_sources<PrimalVariables, AuxiliaryVariables,
SourcesComputer>(sources, vars,
elliptic::first_order_sources<Dim, PrimalVariables, AuxiliaryVariables,
SourcesComputer>(sources, vars, fluxes,
sources_args...);
}
};
Expand Down
59 changes: 40 additions & 19 deletions src/Elliptic/FirstOrderOperator.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -43,34 +43,44 @@ struct FirstOrderFluxesImpl<Dim, tmpl::list<PrimalFields...>,
}
};

template <typename PrimalFields, typename AuxiliaryFields,
template <size_t Dim, typename PrimalFields, typename AuxiliaryFields,
typename SourcesComputer>
struct FirstOrderSourcesImpl;

template <typename... PrimalFields, typename... AuxiliaryFields,
template <size_t Dim, typename... PrimalFields, typename... AuxiliaryFields,
typename SourcesComputer>
struct FirstOrderSourcesImpl<tmpl::list<PrimalFields...>,
struct FirstOrderSourcesImpl<Dim, tmpl::list<PrimalFields...>,
tmpl::list<AuxiliaryFields...>, SourcesComputer> {
template <typename VarsTags, typename... SourcesArgs>
static constexpr void apply(
const gsl::not_null<
Variables<db::wrap_tags_in<::Tags::Source, VarsTags>>*>
sources,
const Variables<VarsTags>& vars,
const Variables<db::wrap_tags_in<
::Tags::Flux, VarsTags, tmpl::size_t<Dim>, Frame::Inertial>>& fluxes,
const SourcesArgs&... sources_args) noexcept {
// Compute sources for primal fields
SourcesComputer::apply(
make_not_null(&get<::Tags::Source<PrimalFields>>(*sources))...,
sources_args..., get<PrimalFields>(vars)...);
// Compute sources for auxiliary fields. They are just the auxiliary field
// values.
// Set auxiliary field sources to the auxiliary field values to begin with.
// This is the standard choice, since the auxiliary equations define the
// auxiliary variables. Source computers can adjust or add to these sources.
tmpl::for_each<tmpl::list<AuxiliaryFields...>>(
[&sources, &vars ](const auto auxiliary_field_tag_v) noexcept {
[&sources, &vars](const auto auxiliary_field_tag_v) noexcept {
using auxiliary_field_tag =
tmpl::type_from<decltype(auxiliary_field_tag_v)>;
get<::Tags::Source<auxiliary_field_tag>>(*sources) =
get<auxiliary_field_tag>(vars);
});
// Call into the sources computer to set primal field sources and possibly
// adjust auxiliary field sources. The sources depend on the primal and the
// auxiliary variables. However, we pass the volume fluxes instead of the
// auxiliary variables to the source computer as an optimization so they
// don't have to be re-computed.
SourcesComputer::apply(
make_not_null(&get<::Tags::Source<PrimalFields>>(*sources))...,
make_not_null(&get<::Tags::Source<AuxiliaryFields>>(*sources))...,
Copy link
Member

Choose a reason for hiding this comment

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

You keep saying "instead of the auxiliary fields" but here you pass the aux fields. Please correct this

Copy link
Member Author

Choose a reason for hiding this comment

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

No this passes the sources for the auxiliary fields so terms can be added to them, such as Christoffel symbol terms (see this comment #2553 (comment)). Does that make sense?

Copy link
Member

Choose a reason for hiding this comment

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

Okay, that makes sense.

I still don't understand why you say "instead of" because the code didn't previously pass the aux vars. Why wouldn't you need the aux vars (I realize they weren't passed before)? Is that because technically they are just (roughly) the derivatives of the primal fields? If so, I'm a bit confused that stability and well-posedness aren't adversely affected by never using the aux variables, though explaining that is beyond the scope of a GitHub comment, I think ;)

Copy link
Member Author

Choose a reason for hiding this comment

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

Previously I didn't pass
(i) the auxiliary sources, because they are "usually" just the auxiliary vars, i.e. when no Christoffel-terms need to be added (e.g. consider the Poisson auxiliary equation -grad(u) + v = 0 where the sources are just v. Remember the auxiliary equation essentially defines the aux vars, which is why the aux vars are always its sources). The aux sources are set to the aux vars just a few lines above.
(ii) the auxiliary variables, because none of the systems that are currently implemented needs to compute non-trivial sources: Both the Poisson and the (linear) elasticity equations are only sourced by a fixed function f(x).

In the comments I say "instead of" because the sources are a function of the variables (primal and auxiliary), but instead of the aux vars I pass the primal fluxes as an optimization, as we discussed before.

Does this clear things up?

Copy link
Member

Choose a reason for hiding this comment

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

Yes, but I think this lengthy discussion suggests that you should either remove the instead of or expand on it in both the docs and the git message

Copy link
Member Author

Choose a reason for hiding this comment

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

Adjusted the commit message and the code comment in the fixup commit, please take a look.

sources_args..., get<PrimalFields>(vars)...,
get<::Tags::Flux<PrimalFields, tmpl::size_t<Dim>, Frame::Inertial>>(
fluxes)...);
}
};

Expand Down Expand Up @@ -118,30 +128,41 @@ auto first_order_fluxes(const Variables<VarsTags>& vars,
* \brief Compute the sources \f$S(u)\f$ for the first-order formulation of
* elliptic systems.
*
Copy link
Member

Choose a reason for hiding this comment

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

add a discussion as to why the sources terms are getting the fluxes

* This function takes the `fluxes` as an argument in addition to the variables
* as an optimization. The fluxes will generally be computed before the sources
* anyway, so we pass them to the source computers to avoid having to re-compute
* them for source-terms that have the same form as the fluxes.
*
* \see `elliptic::first_order_operator`
*/
template <typename PrimalFields, typename AuxiliaryFields,
template <size_t Dim, typename PrimalFields, typename AuxiliaryFields,
typename SourcesComputer, typename VarsTags, typename... SourcesArgs>
void first_order_sources(
gsl::not_null<Variables<db::wrap_tags_in<::Tags::Source, VarsTags>>*>
sources,
const Variables<VarsTags>& vars,
const Variables<db::wrap_tags_in<::Tags::Flux, VarsTags, tmpl::size_t<Dim>,
Frame::Inertial>>& fluxes,
const SourcesArgs&... sources_args) noexcept {
first_order_operator_detail::FirstOrderSourcesImpl<
PrimalFields, AuxiliaryFields, SourcesComputer>::apply(std::move(sources),
vars,
sources_args...);
Dim, PrimalFields, AuxiliaryFields,
SourcesComputer>::apply(std::move(sources), vars, fluxes,
sources_args...);
}

template <typename PrimalFields, typename AuxiliaryFields,
template <size_t Dim, typename PrimalFields, typename AuxiliaryFields,
typename SourcesComputer, typename VarsTags, typename... SourcesArgs>
auto first_order_sources(const Variables<VarsTags>& vars,
const SourcesArgs&... sources_args) noexcept {
auto first_order_sources(
const Variables<VarsTags>& vars,
const Variables<db::wrap_tags_in<::Tags::Flux, VarsTags, tmpl::size_t<Dim>,
Frame::Inertial>>& fluxes,
const SourcesArgs&... sources_args) noexcept {
Variables<db::wrap_tags_in<::Tags::Source, VarsTags>> sources{
vars.number_of_grid_points()};
first_order_operator_detail::FirstOrderSourcesImpl<
PrimalFields, AuxiliaryFields,
SourcesComputer>::apply(make_not_null(&sources), vars, sources_args...);
Dim, PrimalFields, AuxiliaryFields,
SourcesComputer>::apply(make_not_null(&sources), vars, fluxes,
sources_args...);
return sources;
}
// @}
Expand Down
4 changes: 3 additions & 1 deletion src/Elliptic/Systems/Elasticity/Equations.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -94,7 +94,9 @@ struct Sources {
using argument_tags = tmpl::list<>;
static void apply(
const gsl::not_null<tnsr::I<DataVector, Dim>*> source_for_displacement,
const tnsr::I<DataVector, Dim>& /*displacement*/) noexcept {
const gsl::not_null<tnsr::ii<DataVector, Dim>*> /*source_for_strain*/,
const tnsr::I<DataVector, Dim>& /*displacement*/,
Copy link
Member

Choose a reason for hiding this comment

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

aren't the new arguments fluxes? Maybe make that clear in the name?

Copy link
Member Author

Choose a reason for hiding this comment

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

  • The new source_for_strain argument is the source for the auxiliary equation (think the equation grad(u) = v that defines the auxiliary variable for a Poisson equation. The auxiliary variable for the elasticity system is the symmetric "strain", i.e. the symmetrized gradient of the displacement vector). The source computers gain this argument to get the opportunity to add non-principal contributions to the gradient (Christoffel-symbol terms) to the auxiliary equation.
  • The new stress argument is indeed the flux related to the displacement vector field. It's the strain contracted with the constitutive relation of the elastic material, and that's called "stress" in elasticity terms. I think in the context of the elasticity equations it makes sense to call it by its name.

Copy link
Member

Choose a reason for hiding this comment

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

Okay, I can't speak to what would be most intuitive for people used to the elasticity literature, so I'll have to take your word for it :)

const tnsr::IJ<DataVector, Dim>& /*stress*/) noexcept {
for (size_t d = 0; d < Dim; d++) {
source_for_displacement->get(d) = 0.;
}
Expand Down
39 changes: 18 additions & 21 deletions src/Elliptic/Systems/Poisson/Equations.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -16,31 +16,29 @@ namespace Poisson {

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

template <size_t Dim>
void non_euclidean_fluxes(
const gsl::not_null<tnsr::I<DataVector, Dim, Frame::Inertial>*>
flux_for_field,
const tnsr::II<DataVector, Dim, Frame::Inertial>& inv_spatial_metric,
const gsl::not_null<tnsr::I<DataVector, Dim>*> flux_for_field,
const tnsr::II<DataVector, Dim>& inv_spatial_metric,
const Scalar<DataVector>& det_spatial_metric,
const tnsr::i<DataVector, Dim, Frame::Inertial>& field_gradient) noexcept {
const tnsr::i<DataVector, Dim>& field_gradient) noexcept {
raise_or_lower_index(flux_for_field, field_gradient, inv_spatial_metric);
for (size_t i = 0; i < Dim; i++) {
flux_for_field->get(i) *= sqrt(get(det_spatial_metric));
}
}

template <size_t Dim>
void auxiliary_fluxes(gsl::not_null<tnsr::Ij<DataVector, Dim, Frame::Inertial>*>
flux_for_gradient,
const Scalar<DataVector>& field) noexcept {
void auxiliary_fluxes(
gsl::not_null<tnsr::Ij<DataVector, Dim>*> flux_for_gradient,
const Scalar<DataVector>& field) noexcept {
std::fill(flux_for_gradient->begin(), flux_for_gradient->end(), 0.);
for (size_t d = 0; d < Dim; d++) {
flux_for_gradient->get(d, d) = get(field);
Expand All @@ -51,17 +49,16 @@ void auxiliary_fluxes(gsl::not_null<tnsr::Ij<DataVector, Dim, Frame::Inertial>*>

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

#define INSTANTIATE(_, data) \
template void Poisson::euclidean_fluxes<DIM(data)>( \
const gsl::not_null<tnsr::I<DataVector, DIM(data), Frame::Inertial>*>, \
const tnsr::i<DataVector, DIM(data), Frame::Inertial>&) noexcept; \
template void Poisson::non_euclidean_fluxes<DIM(data)>( \
const gsl::not_null<tnsr::I<DataVector, DIM(data), Frame::Inertial>*>, \
const tnsr::II<DataVector, DIM(data), Frame::Inertial>&, \
const Scalar<DataVector>&, \
const tnsr::i<DataVector, DIM(data), Frame::Inertial>&) noexcept; \
template void Poisson::auxiliary_fluxes<DIM(data)>( \
gsl::not_null<tnsr::Ij<DataVector, DIM(data), Frame::Inertial>*>, \
#define INSTANTIATE(_, data) \
template void Poisson::euclidean_fluxes<DIM(data)>( \
const gsl::not_null<tnsr::I<DataVector, DIM(data)>*>, \
const tnsr::i<DataVector, DIM(data)>&) noexcept; \
template void Poisson::non_euclidean_fluxes<DIM(data)>( \
const gsl::not_null<tnsr::I<DataVector, DIM(data)>*>, \
const tnsr::II<DataVector, DIM(data)>&, const Scalar<DataVector>&, \
const tnsr::i<DataVector, DIM(data)>&) noexcept; \
template void Poisson::auxiliary_fluxes<DIM(data)>( \
gsl::not_null<tnsr::Ij<DataVector, DIM(data)>*>, \
const Scalar<DataVector>&) noexcept;

GENERATE_INSTANTIATIONS(INSTANTIATE, (1, 2, 3))
Expand Down
48 changes: 23 additions & 25 deletions src/Elliptic/Systems/Poisson/Equations.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -25,20 +25,19 @@ namespace Poisson {
* equation on a flat spatial metric in Cartesian coordinates.
*/
template <size_t Dim>
void euclidean_fluxes(
gsl::not_null<tnsr::I<DataVector, Dim, Frame::Inertial>*> flux_for_field,
const tnsr::i<DataVector, Dim, Frame::Inertial>& field_gradient) noexcept;
void euclidean_fluxes(gsl::not_null<tnsr::I<DataVector, Dim>*> flux_for_field,
const tnsr::i<DataVector, Dim>& field_gradient) noexcept;

/*!
* \brief Compute the fluxes \f$F^i=\sqrt{\gamma}\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 non_euclidean_fluxes(
gsl::not_null<tnsr::I<DataVector, Dim, Frame::Inertial>*> flux_for_field,
const tnsr::II<DataVector, Dim, Frame::Inertial>& inv_spatial_metric,
gsl::not_null<tnsr::I<DataVector, Dim>*> flux_for_field,
const tnsr::II<DataVector, Dim>& inv_spatial_metric,
const Scalar<DataVector>& det_spatial_metric,
const tnsr::i<DataVector, Dim, Frame::Inertial>& field_gradient) noexcept;
const tnsr::i<DataVector, Dim>& field_gradient) noexcept;

/*!
* \brief Compute the fluxes \f$F^i_j=\delta^i_j u(x)\f$ for the auxiliary
Expand All @@ -47,9 +46,9 @@ void non_euclidean_fluxes(
* \see Poisson::FirstOrderSystem
*/
template <size_t Dim>
void auxiliary_fluxes(gsl::not_null<tnsr::Ij<DataVector, Dim, Frame::Inertial>*>
flux_for_gradient,
const Scalar<DataVector>& field) noexcept;
void auxiliary_fluxes(
gsl::not_null<tnsr::Ij<DataVector, Dim>*> flux_for_gradient,
const Scalar<DataVector>& field) noexcept;

/*!
* \brief Compute the fluxes \f$F^i_A\f$ for the Poisson equation on a flat
Expand All @@ -61,15 +60,12 @@ template <size_t Dim>
struct EuclideanFluxes {
using argument_tags = tmpl::list<>;
static void apply(
const gsl::not_null<tnsr::I<DataVector, Dim, Frame::Inertial>*>
flux_for_field,
const tnsr::i<DataVector, Dim, Frame::Inertial>&
field_gradient) noexcept {
const gsl::not_null<tnsr::I<DataVector, Dim>*> flux_for_field,
const tnsr::i<DataVector, Dim>& field_gradient) noexcept {
euclidean_fluxes(flux_for_field, field_gradient);
}
static void apply(
const gsl::not_null<tnsr::Ij<DataVector, Dim, Frame::Inertial>*>
flux_for_gradient,
const gsl::not_null<tnsr::Ij<DataVector, Dim>*> flux_for_gradient,
const Scalar<DataVector>& field) noexcept {
auxiliary_fluxes(flux_for_gradient, field);
}
Expand All @@ -89,19 +85,16 @@ struct NonEuclideanFluxes {
gr::Tags::InverseSpatialMetric<Dim, Frame::Inertial, DataVector>,
gr::Tags::DetSpatialMetric<DataVector>>;
static void apply(
const gsl::not_null<tnsr::I<DataVector, Dim, Frame::Inertial>*>
flux_for_field,
const tnsr::II<DataVector, Dim, Frame::Inertial>& inv_spatial_metric,
const gsl::not_null<tnsr::I<DataVector, Dim>*> flux_for_field,
const tnsr::II<DataVector, Dim>& inv_spatial_metric,
const Scalar<DataVector>& det_spatial_metric,
const tnsr::i<DataVector, Dim, Frame::Inertial>&
field_gradient) noexcept {
const tnsr::i<DataVector, Dim>& field_gradient) noexcept {
non_euclidean_fluxes(flux_for_field, inv_spatial_metric, det_spatial_metric,
field_gradient);
}
static void apply(
const gsl::not_null<tnsr::Ij<DataVector, Dim, Frame::Inertial>*>
flux_for_gradient,
const tnsr::II<DataVector, Dim, Frame::Inertial>& /*inv_spatial_metric*/,
const gsl::not_null<tnsr::Ij<DataVector, Dim>*> flux_for_gradient,
const tnsr::II<DataVector, Dim>& /*inv_spatial_metric*/,
const Scalar<DataVector>& /*det_spatial_metric*/,
const Scalar<DataVector>& field) noexcept {
auxiliary_fluxes(flux_for_gradient, field);
Expand All @@ -117,8 +110,13 @@ struct NonEuclideanFluxes {
*/
struct Sources {
using argument_tags = tmpl::list<>;
static void apply(const gsl::not_null<Scalar<DataVector>*> source_for_field,
const Scalar<DataVector>& /*field*/) noexcept {
template <size_t Dim>
static void apply(
const gsl::not_null<Scalar<DataVector>*> source_for_field,
const gsl::not_null<
tnsr::i<DataVector, Dim>*> /*source_for_field_gradient*/,
const Scalar<DataVector>& /*field*/,
const tnsr::I<DataVector, Dim>& /*field_flux*/) noexcept {
get(*source_for_field) = 0.;
}
};
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -86,8 +86,12 @@ struct Fluxes {

struct Sources {
using argument_tags = tmpl::list<>;
static void apply(const gsl::not_null<Scalar<DataVector>*> source_for_field,
const Scalar<DataVector>& field) {
template <size_t Dim>
static void apply(
const gsl::not_null<Scalar<DataVector>*> source_for_field,
const gsl::not_null<tnsr::i<DataVector, Dim>*> /*source_for_aux_field*/,
const Scalar<DataVector>& field,
const tnsr::I<DataVector, Dim>& /*field_flux*/) {
get(*source_for_field) = get(field);
}
};
Expand Down
21 changes: 14 additions & 7 deletions tests/Unit/Elliptic/Systems/Elasticity/Test_Equations.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -63,6 +63,7 @@ void test_computers(const DataVector& used_for_size) {
using auxiliary_flux_tag =
::Tags::Flux<auxiliary_field_tag, tmpl::size_t<Dim>, Frame::Inertial>;
using field_source_tag = ::Tags::Source<field_tag>;
using auxiliary_source_tag = ::Tags::Source<auxiliary_field_tag>;
using constitutive_relation_tag = Elasticity::Tags::ConstitutiveRelation<
Elasticity::ConstitutiveRelations::IsotropicHomogeneous<Dim>>;
using coordinates_tag = domain::Tags::Coordinates<Dim, Frame::Inertial>;
Expand Down Expand Up @@ -103,17 +104,23 @@ void test_computers(const DataVector& used_for_size) {
}
{
INFO("Sources" << Dim << "D");
auto box = db::create<db::AddSimpleTags<field_tag, field_source_tag>>(
tnsr::I<DataVector, Dim>{num_points,
std::numeric_limits<double>::signaling_NaN()},
tnsr::I<DataVector, Dim>{num_points,
std::numeric_limits<double>::signaling_NaN()});
auto box =
db::create<db::AddSimpleTags<field_source_tag, auxiliary_source_tag>>(
tnsr::I<DataVector, Dim>{
num_points, std::numeric_limits<double>::signaling_NaN()},
tnsr::ii<DataVector, Dim>{
num_points, std::numeric_limits<double>::signaling_NaN()});

const Elasticity::Sources<Dim> sources_computer{};
using argument_tags = typename Elasticity::Sources<Dim>::argument_tags;

db::mutate_apply<tmpl::list<field_source_tag>, argument_tags>(
sources_computer, make_not_null(&box), get<field_tag>(box));
db::mutate_apply<tmpl::list<field_source_tag, auxiliary_source_tag>,
argument_tags>(
sources_computer, make_not_null(&box),
tnsr::I<DataVector, Dim>{num_points,
std::numeric_limits<double>::signaling_NaN()},
tnsr::IJ<DataVector, Dim>{
num_points, std::numeric_limits<double>::signaling_NaN()});
auto expected_field_source = tnsr::I<DataVector, Dim>{num_points, 0.};
CHECK(get<field_source_tag>(box) == expected_field_source);
}
Expand Down