Merge pull request #2553 from nilsleiffischer/pass_all_vars_to_ellipt…

…ic_computers

Pass fluxes to elliptic source computers

src/Elliptic/DiscontinuousGalerkin/InitializeFirstOrderOperator.hpp

@@ -71,7 +71,7 @@ struct InitializeFirstOrderOperator {
                                               vars_tag, PrimalVariables,
                                               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

src/Elliptic/FirstOrderComputeTags.hpp

@@ -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...);
   }
 };

src/Elliptic/FirstOrderOperator.hpp

@@ -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))...,
+        sources_args..., get<PrimalFields>(vars)...,
+        get<::Tags::Flux<PrimalFields, tmpl::size_t<Dim>, Frame::Inertial>>(
+            fluxes)...);
   }
 };
 
@@ -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.
  *
+ * 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;
 }
 // @}

src/Elliptic/Systems/Elasticity/Equations.hpp

@@ -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*/,
+      const tnsr::IJ<DataVector, Dim>& /*stress*/) noexcept {
     for (size_t d = 0; d < Dim; d++) {
       source_for_displacement->get(d) = 0.;
     }

src/Elliptic/Systems/Poisson/Equations.cpp

@@ -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);
@@ -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))

src/Elliptic/Systems/Poisson/Equations.hpp

@@ -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
@@ -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
@@ -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);
   }
@@ -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);
@@ -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.;
   }
 };

tests/Unit/Elliptic/DiscontinuousGalerkin/Test_InitializeFirstOrderOperator.cpp

@@ -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);
   }
 };

tests/Unit/Elliptic/Systems/Elasticity/Test_Equations.cpp

@@ -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>;
@@ -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);
   }