Specialized binary functions for Tensor<1, 1, Number> and Number

[peterrum]

@mschreter and I would like to propose following binary function specializations to be able operate on dealii::Tensor<1, 1, dealii::VectorizedArray<Number, N>> and dealii::VectorizedArray<Number, N> easily:

namespace dealii
{
  template <typename Number, std::size_t N>
  dealii::Tensor<1, 1, dealii::VectorizedArray<Number, N>>
  operator+=(const dealii::Tensor<1, 1, dealii::VectorizedArray<Number, N>> &vec,
             const dealii::VectorizedArray<Number, N> &                      scalar)
  {
    auto temp = vec;
    temp[0] += scalar;

    return temp;
  }

  template <typename Number, std::size_t N>
  dealii::Tensor<1, 1, dealii::VectorizedArray<Number, N>>
  operator-=(const dealii::Tensor<1, 1, dealii::VectorizedArray<Number, N>> &vec,
             const dealii::VectorizedArray<Number, N> &                      scalar)
  {
    auto temp = vec;
    temp[0] -= scalar;

    return temp;
  }

  template <typename Number, std::size_t N>
  dealii::VectorizedArray<Number, N>
  scalar_product(const dealii::VectorizedArray<Number, N> &                      scalar,
                 const dealii::Tensor<1, 1, dealii::VectorizedArray<Number, N>> &vec)
  {
    return vec[0] * scalar;
  }

  template <typename Number, std::size_t N>
  dealii::VectorizedArray<Number, N>
  scalar_product(const dealii::Tensor<1, 1, dealii::VectorizedArray<Number, N>> &vec,
                 const dealii::VectorizedArray<Number, N> &                      scalar)
  {
    return vec[0] * scalar;
  }
} // namespace dealii

We use such specializations in the context of dim-independent MatrixFree programming when we also want to enable algorithms in 1D without needing to write the specializations (see: https://github.com/MeltPoolDG/MeltPoolDG/pull/98). We are not sure if this should be part of deal.II and if yes what the appropriate place is: vectoriazation.h or tensor.h or somewhere else!

Note: I guess the above discussion is also true for dealii::Tensor<1, 1, Number> and Number with Number being double/float.

[bangerth]

I don't particularly like this kind of change. My objection is that these things represent different concepts. A Tensor<1,1> is a vector that happens to live in a 1d world. A double on the other hand is a scalar. Nothing good can come off if you add a scalar to a vector -- you can probably define what is supposed to happen in that case, but it is conflating concepts and we have generally been quite careful in deal.II to concisely define concepts in the sense that, for example, shape_grad() in 1d returns a Tensor<1,1> and not a double. I think you will also realize that the existence of these kinds of overloads makes it more difficult to write truly dimension-independent codes.

So, in summary, I'm not in favor.

[peterrum]

@bangerth I agree with your comments. But, I am afraid the current interface of MatrixFree and FEEvaluation makes dimension- and component-independent programing unnecessarily complex due to many specializations of the FEEvaluation class (see https://github.com/dealii/dealii/blob/master/include/deal.II/matrix_free/fe_evaluation.h). My guess is that the interface was chosen in such a way that it is nicely usable and to be efficient (which it is).

Let me quickly show the some specializations of the functions FEEvaluation::get_value():

FEEvaluation<dim, fe_degree, n_q_points_1d, n_components>::get_value() -> Tensor<1, n_components>
FEEvaluation<dim, fe_degree, n_q_points_1d, 1>::get_value()            -> double

and FEEvaluation::get_gradient():

FEEvaluation<dim, fe_degree, n_q_points_1d, n_components>::get_gradient() -> Tensor<1, n_components_, Tensor<1, dim>>;
FEEvaluation<dim, fe_degree, n_q_points_1d, 1>::get_gradient()            -> Tensor<1, dim>;
FEEvaluation<dim, fe_degree, n_q_points_1d, dim>::get_gradient()          -> Tensor<2, dim>;

Now lets assume that want to compute scalar_product(FEEvaluation<dim, fe_degree, n_q_points_1d, dim>::get_value(), FEEvaluation<dim, fe_degree, n_q_points_1d, 1>::get_gradient()) (as in the context of scalar transport problems). This will work for 2D and 3D, since the return type of both functions is Tensor<1, dim>. But for 1D, we get the case that one is returning a double and the other function a Tensor<1, 1>.

There are other instances where this might happen.

Together with @mschreter, we used to have a free function that converts a double to a tensor:

    template <int dim, typename number = double>
    static Tensor<1, dim, VectorizedArray<number>>
    convert_to_vector(const VectorizedArray<number> &in)
    {
      Tensor<1, dim, VectorizedArray<number>> vec;

      for (unsigned int v = 0; v < VectorizedArray<number>::size(); ++v)
        vec[0][v] = in[v];

      return vec;
    }

    template <int dim, typename number = double>
    static Tensor<1, dim, VectorizedArray<number>>
    convert_to_vector(const Tensor<1, dim, VectorizedArray<number>> &in)
    {
      return in;
    }

which had to be called explicitly (and made dimension-independent matrix-free code hard to read).

This is why we introduced the specializations (proposed in the description of this PR) in our user code; since we like it there, we would like to make part of the library.

[bangerth]

Now lets assume that want to compute |scalar_product(FEEvaluation<dim, fe_degree, n_q_points_1d, dim>::get_value(), FEEvaluation<dim, fe_degree, n_q_points_1d, 1>::get_gradient())| (as in the context of scalar transport problems). This will work for 2D and 3D, since the return type of both functions is |Tensor<1, dim>|. But for 1D, we get the case that one is returning a |double| and the other function a |Tensor<1, 1>|.

There lies your problem. Let the 1d function return |Tensor<1,1>|. I continue to think that softening semantic notions for the purpose of convenience is the wrong approach. I get that the code you have to write for dimension-independent programming is cumbersome, but I think you are addressing the problem at the wrong end. If the MatrixFree framework is giving you the wrong/inconvenient types, then change things at that end rather than in the downstream places. If that happens to require incompatible changes, then so be it -- better an incompatibility now than inconsistent semantics forever.

[peterrum]

Let's close this PR. 10 years ago, decisions were made to make MatrixFree as fast as possible for many configurations 1D/2D/3D and scalar/vectorial. The style is inconsistent with the rest of the library, which is not tuned for performance. This inconsistency cannot be changed anymore. All we wanted is to improve data conversion at a reasonable place, which is not performance critical, to be able to simply switch between dimensions (1D to 3D and back). Currently, it is NOT possible to write nice dimension-independent code with MatrixFree, which ends up having to write specializations like this ones above in user code. It would have been nice to move such specializations from our user code to a place where everyone profits from.