From ee48514cf06abab48e4bc3ab1345fb29a4d5b159 Mon Sep 17 00:00:00 2001
From: Nils Leif Fischer <nils.fischer@aei.mpg.de>
Date: Mon, 1 Feb 2021 12:19:01 +0100
Subject: [PATCH] Loosen restrictions on divergence return buffer tags

Allow the divergence function to write the result into any Variables
that holds compatible tags, not only ::Tags::div<FluxTags>.
---
 .../LinearOperators/Divergence.hpp            | 16 +++++++-------
 .../LinearOperators/Divergence.tpp            | 21 ++++++++++---------
 2 files changed, 18 insertions(+), 19 deletions(-)

diff --git a/src/NumericalAlgorithms/LinearOperators/Divergence.hpp b/src/NumericalAlgorithms/LinearOperators/Divergence.hpp
index 43610a23cfeb..be9cd8e3230b 100644
--- a/src/NumericalAlgorithms/LinearOperators/Divergence.hpp
+++ b/src/NumericalAlgorithms/LinearOperators/Divergence.hpp
@@ -63,11 +63,11 @@ auto divergence(
         inverse_jacobian) noexcept
     -> Variables<db::wrap_tags_in<Tags::div, FluxTags>>;
 
-template <typename FluxTags, size_t Dim, typename DerivativeFrame>
+template <typename... DivTags, typename... FluxTags, size_t Dim,
+          typename DerivativeFrame>
 void divergence(
-    gsl::not_null<Variables<db::wrap_tags_in<Tags::div, FluxTags>>*>
-        divergence_of_F,
-    const Variables<FluxTags>& F, const Mesh<Dim>& mesh,
+    gsl::not_null<Variables<tmpl::list<DivTags...>>*> divergence_of_F,
+    const Variables<tmpl::list<FluxTags...>>& F, const Mesh<Dim>& mesh,
     const InverseJacobian<DataVector, Dim, Frame::Logical, DerivativeFrame>&
         inverse_jacobian) noexcept;
 // @}
@@ -127,11 +127,9 @@ struct DivVariablesCompute : db::add_tag_prefix<div, Tag>, db::ComputeTag {
  public:
   using base = db::add_tag_prefix<div, Tag>;
   using return_type = typename base::type;
-  static constexpr void (*function)(const gsl::not_null<return_type*>,
-                                    const typename Tag::type&, const Mesh<dim>&,
-                                    const typename InverseJacobianTag::type&) =
-      divergence<typename Tag::type::tags_list, dim,
-                 typename tmpl::back<inv_jac_indices>::Frame>;
+  static constexpr void (*function)(
+      const gsl::not_null<return_type*>, const typename Tag::type&,
+      const Mesh<dim>&, const typename InverseJacobianTag::type&) = divergence;
   using argument_tags =
       tmpl::list<Tag, domain::Tags::Mesh<dim>, InverseJacobianTag>;
 };
diff --git a/src/NumericalAlgorithms/LinearOperators/Divergence.tpp b/src/NumericalAlgorithms/LinearOperators/Divergence.tpp
index 684cd2c1c5e1..b66e451fd74d 100644
--- a/src/NumericalAlgorithms/LinearOperators/Divergence.tpp
+++ b/src/NumericalAlgorithms/LinearOperators/Divergence.tpp
@@ -22,11 +22,11 @@ Variables<db::wrap_tags_in<Tags::div, FluxTags>> divergence(
   return divergence_of_F;
 }
 
-template <typename FluxTags, size_t Dim, typename DerivativeFrame>
+template <typename... DivTags, typename... FluxTags, size_t Dim,
+          typename DerivativeFrame>
 void divergence(
-    const gsl::not_null<Variables<db::wrap_tags_in<Tags::div, FluxTags>>*>
-        divergence_of_F,
-    const Variables<FluxTags>& F, const Mesh<Dim>& mesh,
+    const gsl::not_null<Variables<tmpl::list<DivTags...>>*> divergence_of_F,
+    const Variables<tmpl::list<FluxTags...>>& F, const Mesh<Dim>& mesh,
     const InverseJacobian<DataVector, Dim, Frame::Logical, DerivativeFrame>&
         inverse_jacobian) noexcept {
   if (UNLIKELY(divergence_of_F->number_of_grid_points() !=
@@ -35,13 +35,13 @@ void divergence(
   }
 
   const auto logical_partial_derivatives_of_F =
-      logical_partial_derivatives<FluxTags>(F, mesh);
+      logical_partial_derivatives<tmpl::list<FluxTags...>>(F, mesh);
 
-  tmpl::for_each<FluxTags>([
+  const auto apply_div = [
     &divergence_of_F, &inverse_jacobian, &logical_partial_derivatives_of_F
-  ](auto tag) noexcept {
-    using FluxTag = tmpl::type_from<decltype(tag)>;
-    using DivFluxTag = Tags::div<FluxTag>;
+  ](auto flux_tag_v, auto div_tag_v) noexcept {
+    using FluxTag = std::decay_t<decltype(flux_tag_v)>;
+    using DivFluxTag = std::decay_t<decltype(div_tag_v)>;
 
     using first_index =
         tmpl::front<typename FluxTag::type::index_list>;
@@ -66,5 +66,6 @@ void divergence(
         }
       }
     }
-  });
+  };
+  EXPAND_PACK_LEFT_TO_RIGHT(apply_div(FluxTags{}, DivTags{}));
 }