Integrators for DG face integrals with tests. #5870

kkormann · 2018-02-06T20:10:18Z

Integrator part of #5667.

New functionality EvaluatorTensorProduct::apply_face providing face integrals. Introduces structs FEFaceEvaluationImpl for operations within the face and FEFaceNormalEvaluationImpl for interpolation from the cell to the face.
Changed template parameters in EvaluatorTensorProduct from fe_degree and n_q_points_1d to n_rows and n_columns because the kernels are a linear algebra abstraction and do not necessarily need to be associated to quadrature and elements.
Prepare EvaluatorTensorProduct for dimension-independent case (dim>3).
Add DEAL_II_RESTRICT to main worker kernels where appropriate.
Ensure that EvaluatorTensorProduct::apply also works when in and out point to same data.
Tests for the various in-place cases of EvaluatorTensorProduct.
Make number type of FEEvaluationImpl a true abstract type Number, rather than constructing VectorizedArray<Number>. This removes the hardcoding of VectorizedArray in this class.
Extend FEEvaluationImplTransformToCollocation for n_q_points_1d > fe_degree+1.

Since all changes are in the internal namespace, breaking the old interface is uncritical.

kronbichler · 2018-02-06T20:10:53Z

/run-tests

kronbichler · 2018-02-07T08:52:18Z

The evaluator part of face integrals FEFaceEvaluationImpl and FEFaceNormalEvaluationImpl will be covered by tests that come once the whole chain of features listed in #5667 is supported, i.e., once the evaluators are introduced in a class FEFaceEvaluation and we actually have high-level support for them with mapping data and indices.

Rombur · 2018-02-07T14:28:57Z

I'd like to review this patch but it will probably take me a couple of days to get to it.

bangerth · 2018-02-07T14:40:25Z

I'm glad @Rombur wants to do it, because 3800 new lines would have blown a substantial hole into one of my next days :-)

kronbichler · 2018-02-07T14:58:51Z

@Rombur great!

@kkormann or I might open another PR with other unrelated functionality in the meantime, though.

tjhei

I only skimmed it and I hope somebody else does a thorough review, but I have a few comments...

tjhei · 2018-02-09T17:42:37Z

source/multigrid/mg_transfer_matrix_free.cc

+              }
+            else
+              {
+                evaluator.template values<2,false,false>(t0, t2);


can you please leave a comment to explain this reverse ordering? (also above)

It is because the EvaluatorTensorProduct kernels have gained the possibility to do in-place operations, which necessitated a change in the hardcoded strides/offsets in the generalized matrix-matrix multiplication to make sure it does not overwrite things. In particular, we need to start with the outermost dimension in the no-transpose path with more columns than rows (second false template argument). By the way, this reminds me that we could now do all this in a single array t0 without temporary arrays t1,t2.

We have replaced the whole function by the new FEEvaluationImplBasisChange::do_forward/do_backward call.

tjhei · 2018-02-09T17:49:26Z

include/deal.II/matrix_free/evaluation_kernels.h

+                                       Utilities::fixed_power<dim>(shape_info.fe_degree+1) : shape_info.dofs_per_component_on_cell;
+    Number *values_dofs = (type == MatrixFreeFunctions::truncated_tensor) ?
+                          scratch_data+2*(std::max(shape_info.dofs_per_component_on_cell, shape_info.n_q_points)) :
+                          const_cast<Number *> (values_dofs_actual);


I don't like the const cast, wouldn't it be cleaner to move this variable below the if and make it const and operate on a different (non-const) variable inside the if? (this assumes this is the only place you modify it)

note: I hope I understood the code correctly.

We need to check this carefully. I agree that it is not pretty right now.

tjhei · 2018-02-09T17:52:34Z

include/deal.II/matrix_free/evaluation_kernels.h

+   * evaluation, spectral collocation or simply collocation, meaning the same
+   * location for shape functions and evaluation space (quadrature points).
+   *
+   * @author Katharina Kormann, 2012


you have been sitting on this since 2012? ;-)

No, it is just github that gets confused by moving this class FEEvaluationImplCollocation above FEEvaluationImplTransformToCollocation - the latter now uses the former after having changed the basis, as promised by the name.

And just as a side note, the changes in that file are actually mostly 2017.

tjhei · 2018-02-09T17:57:02Z

include/deal.II/matrix_free/evaluation_kernels.h

+                                     EvalGeneric::dofs_per_cell : EvalGeneric::n_q_points;
+      static_assert(temp_size > 0, "temp_size should not be zero");
+
+      Number  temp_data[temp_size < 100 ? temp_size : 1];


can you use boost::container::small_vector here instead?

Fixed (with using the provided scratch_data array rather than a separate one).

Rombur

I am not done reviewing this PR yet. This only the first part.

Rombur · 2018-02-11T19:57:01Z

include/deal.II/matrix_free/evaluation_kernels.h

            if (evaluate_hessians == true)
-              eval.template hessians<0,true,false> (values_dofs[c], hessians_quad[c][0]);
+              eval.template hessians<0,true,false> (values_dofs, hessians_quad);
+            values_dofs += dofs_per_comp;


Can you add a comment that we are moving in the 1D array. This was confusing at first.

Rombur · 2018-02-11T20:24:45Z

include/deal.II/matrix_free/evaluation_kernels.h


    switch (dim)
      {
      case 1:
        for (unsigned int c=0; c<n_components; c++)
          {
            if (integrate_values == true)
-              eval.template values<0,false,false> (values_quad[c], values_dofs[c]);
+              {
+                if (overwrite_values == true)


Is overwrite_values the appropriate name? Is it just the opposite of the add template parameters in the tensor product kernels? If so, I think that it would makes sense to rename the variable because it is unclear what overwrite_values==false means.

Also you may want to make overwrite_values a template parameter. I feel that it would simplify the code quite a bit here.

Changed the name to add_into_values_array which better reflects the purpose. We decided not to have a function template within the already templated class for now because it is shorter to write at the call.

Rombur · 2018-02-11T20:52:07Z

include/deal.II/matrix_free/evaluation_kernels.h

+              ::do_forward(transformation_matrix,
+                           values_in + (q-1)*Utilities::fixed_int_power<n_points_1,next_dim>::value,
+                           values_out + (q-1)*Utilities::fixed_int_power<n_points_2,next_dim>::value);
+          if (n_points_2 == n_points_1 && dim == 2)


So if I understand this correctly, the end of the recursion is for dim==1 except if n_points_2==n_points_1 && dim==2 where we stop early.

Rombur · 2018-02-11T21:23:51Z

include/deal.II/matrix_free/evaluation_kernels.h

-    const unsigned int d3 = d1+d2;
-    const unsigned int d4 = dim>2?4:d3;
-    const unsigned int d5 = dim>2?5:d4;
+    Assert(n_q_points_1d > fe_degree,


Shouldn't this be a static assert?

The compiler does not like the static assert because the function where we select the different paths sees the struct also in the invalid case - even though the code then does not enter the function.

Rombur · 2018-02-11T21:25:47Z

include/deal.II/matrix_free/evaluation_kernels.h

-    // gradients_quad[2] only for dim==3.
-    const unsigned int d1 = dim>1?1:0;
-    const unsigned int d2 = dim>2?2:0;
+    Assert(n_q_points_1d > fe_degree,


static assert?

Rombur · 2018-02-11T21:32:41Z

include/deal.II/matrix_free/evaluation_kernels.h

+                            const Number *input,
+                            Number *output,
+                            const bool               do_gradients,
+                            const unsigned int       face_no)


The indentation is pretty strange here

Rombur · 2018-02-11T21:39:11Z

include/deal.II/matrix_free/evaluation_selector.h

+          ::evaluate(shape_info, values_dofs_actual, values_quad,
+                     gradients_quad, hessians_quad, scratch_data,
+                     evaluate_values, evaluate_gradients, evaluate_hessians);
+        else if (degree < n_q_points_1d)


Is this the same condition as before or is it more general?

Before we only did it for degree+1==n_q_points_1d but now we do it for degree<n_q_points_1d. It works because a lower polynomial degree can be exactly embedded into a space of a higher degree.

jwitte08 · 2018-02-16T16:43:41Z

include/deal.II/matrix_free/tensor_product_kernels.h

+      Assert(shape_values.size() == 0 ||
+             shape_values.size() == n_rows*n_columns ||
+             shape_values.size() % n_rows == 0,
+             ExcDimensionMismatch(shape_values.size(), n_rows*n_columns));


Assert can be reduced to shape_values.size() % n_rows == 0?

Can you elaborate on shape_values.size() % n_rows == 0?

jwitte08 · 2018-02-16T16:44:30Z

include/deal.II/matrix_free/tensor_product_kernels.h

-    static const unsigned int dofs_per_cell = Utilities::fixed_int_power<fe_degree+1,dim>::value;
-    static const unsigned int n_q_points = Utilities::fixed_int_power<n_q_points_1d,dim>::value;
+    static const unsigned int dofs_per_cell = Utilities::fixed_int_power<n_rows,dim>::value;
+    static const unsigned int n_q_points = Utilities::fixed_int_power<n_columns,dim>::value;


Do we need dofs_per_cell and n_q_points only in the context of FEEvaluation?

Right now it's only used there, but again it is a general concept. It would probably make sense to give the variables more descriptive names now that we don't use fe_degree any more. How about n_rows_of_product and n_columns_of_product to indicate that it is the size accumulated over all dimensions, n_rows^dim?

I like that! It generically describes the size of the resulting matrix with respect to the underlying Kronecker product.

Do you want to change the name of the variables in this PR? If not, I'll merge.

It is fixed now. The PR should be ready for merging.

jwitte08 · 2018-02-16T16:46:53Z

include/deal.II/matrix_free/tensor_product_kernels.h

+      Assert(shape_values.size() == 0 ||
+             shape_values.size() == n_rows*n_columns ||
+             shape_values.size() % n_rows == 0,
+             ExcDimensionMismatch(shape_values.size(), n_rows*n_columns));


Can you elaborate on shape_values.size() % n_rows == 0?

jwitte08 · 2018-02-16T17:04:46Z

include/deal.II/matrix_free/tensor_product_kernels.h

+     * not correct.
+     *
+     * @tparam direction Direction that is evaluated
+     * @tparam contract_over_rows If true, the tensor contraction goes sums


typo ~~goes~~

jwitte08 · 2018-02-16T17:15:06Z

include/deal.II/matrix_free/tensor_product_kernels.h

+     * the one currently applied. In other words, apply_face() must be called
+     * before calling any interpolation within the face.
+     *
+     * @tparam direction Direction of the normal vector (0=x, 1=y, etc)


is it better to call it face_direction?

Fixed. (This was a typo in the comment, the actual parameter was already named this way.)

jwitte08 · 2018-02-16T17:32:04Z

include/deal.II/matrix_free/tensor_product_kernels.h

+     * @tparam contract_onto_face If true, the input vector is of size n_rows^dim
+     *                            and interpolation into n_rows^(dim-1) points
+     *                            is performed. This is a typical scenario in
+     *                            FEEvaluation::evaluate() calls. If false, the


Do you mean FEFaceEvaluation::evaluate() here?

jwitte08 · 2018-02-16T17:36:30Z

include/deal.II/matrix_free/tensor_product_kernels.h

-  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
-  struct EvaluatorTensorProduct<evaluate_general,dim,fe_degree,n_q_points_1d,Number>
+  template <int dim, int n_rows, int n_columns, typename Number, typename Number2>
+  struct EvaluatorTensorProduct<evaluate_general,dim,n_rows,n_columns,Number,Number2>


Great! I like the generalization to n_rows and n_columns.

Rombur

Looks good. Next time please split the patch in more PR/commits it's really hard to review such a large PR

Rombur · 2018-02-18T22:03:53Z

include/deal.II/matrix_free/tensor_product_kernels.h

+   * @tparam variant Variant of evaluation used for creating template
+   *                 specializations
+   * @tparam dim Dimension of the function
+   * @tparam n_rows Number rows in the transformation matrix, which corresponds


Number of rows

Rombur · 2018-02-18T22:04:13Z

include/deal.II/matrix_free/tensor_product_kernels.h

+   * tensor product form of the basis functions.
+   *
+   * @tparam dim Space dimension in which this class is applied
+   * @tparam n_rows Number rows in the transformation matrix, which corresponds


Number of rows

Rombur · 2018-02-18T22:09:22Z

include/deal.II/matrix_free/tensor_product_kernels.h

+     *                            together
+     * @tparam add If true, the result is added to the output vector, else
+     *             the computed values overwrite the content in the output
+     * @tparam max_derivatve Sets the number of derivatives that should be


typo max_derivative

Rombur · 2018-02-18T22:45:31Z

include/deal.II/matrix_free/tensor_product_kernels.h

@@ -456,48 +799,70 @@ namespace internal
   * tensor-product based elements for "symmetric" finite elements, i.e., when
   * the shape functions are symmetric about 0.5 and the quadrature points
   * are, too.
+   *
+   * @tparam dim Space dimension in which this class is applied
+   * @tparam n_rows Number rows in the transformation matrix, which corresponds


Number of rows

Rombur · 2018-02-18T22:49:41Z

include/deal.II/matrix_free/tensor_product_kernels.h

@@ -1032,99 +1340,200 @@ namespace internal
   * experiments in the book say that the method is not efficient for N<20, it
   * is more efficient in the context where the loop bounds are compile-time
   * constants (templates).
+   *
+   * @tparam dim Space dimension in which this class is applied
+   * @tparam n_rows Number rows in the transformation matrix, which corresponds


Rombur · 2018-02-18T22:52:40Z

include/deal.II/matrix_free/tensor_product_kernels.h

+      shape_gradients (nullptr),
+      shape_hessians (nullptr)
+    {
+      AssertDimension(shape_values.size(), n_rows*((n_columns+1)/2));


Does this assert work if shape_values is empty? You check that it is not empty below.

In this case, it makes no sense if shape_values is empty since it is the only argument.

I see, this function is internal right? If it is then it's fine otherwise we should check that shape_values is not empty, so that user get a nicer error message.

Yes, the function is internal (at least for now). We might think about making some of those generally usable at some point in the future but certainly not for the 9.0 release.

Rombur · 2018-02-18T23:05:24Z

include/deal.II/matrix_free/tensor_product_kernels.h

+   * difference is in the way the input and output quantities are symmetrized.
+   *
+   * @tparam dim Space dimension in which this class is applied
+   * @tparam n_rows Number rows in the transformation matrix, which corresponds


Rombur · 2018-02-18T23:06:37Z

include/deal.II/matrix_free/tensor_product_kernels.h

+     * Empty constructor. Does nothing. Be careful when using 'values' and
+     * related methods because they need to be filled with the other
+     * constructor passing in at least an array for the values.
+     */


Does this mean that the object constructed is never usable?

You can only use the (static) function apply.

kronbichler · 2018-02-20T09:37:14Z

Next time please split the patch in more PR/commits it's really hard to review such a large PR

This PR was already split out from the overall matrix-free DG support as laid out in #5667 and it would have been difficult to split this up further without considerably re-working on the code. I think we have a few smaller patches ahead though.

On the plus side, I would argue that pulling in code from a development branch that we have worked on for more than three years has the advantage that the master branch has skipped a lot of initially bad decisions that we made under the way and thus reduced the noise so far.

Rombur · 2018-02-20T13:20:37Z

This PR was already split out from the overall matrix-free DG support as laid out in #5667 and it would have been difficult to split this up further without considerably re-working on the code. I think we have a few smaller patches ahead though.

I understand but for example replacing VectorizedArray by Number could have been done in its own commit. Then it is easier to see what is really important but I agree that it's not always possible to split a big PR.

davydden added the Matrix-free label Feb 6, 2018

tjhei reviewed Feb 9, 2018

View reviewed changes

Rombur reviewed Feb 11, 2018

View reviewed changes

kkormann force-pushed the issue_5667_integrators branch from 4c45d17 to 296a914 Compare February 12, 2018 22:21

jwitte08 reviewed Feb 16, 2018

View reviewed changes

drwells mentioned this pull request Feb 17, 2018

Deprecate fixed_int_power #5919

Merged

kkormann force-pushed the issue_5667_integrators branch 2 times, most recently from 779cf96 to 8dc3978 Compare February 18, 2018 20:35

Rombur reviewed Feb 18, 2018

View reviewed changes

kkormann force-pushed the issue_5667_integrators branch from 8dc3978 to 34609c3 Compare February 19, 2018 20:35

davydden mentioned this pull request Feb 20, 2018

Group matrix-free cell blocks based on predicate (ie material ID) #5197

Closed

Integrators for DG face integrals with tests.

13fc76d

kkormann force-pushed the issue_5667_integrators branch from 34609c3 to 13fc76d Compare February 21, 2018 07:13

Rombur merged commit 8cb62c3 into dealii:master Feb 21, 2018

Integrators for DG face integrals with tests. #5870

Integrators for DG face integrals with tests. #5870

Conversation

kkormann commented Feb 6, 2018

kronbichler commented Feb 6, 2018

kronbichler commented Feb 7, 2018

Rombur commented Feb 7, 2018

bangerth commented Feb 7, 2018

kronbichler commented Feb 7, 2018

tjhei left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Rombur left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

kkormann Feb 12, 2018 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

kkormann Feb 18, 2018 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Rombur left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

kronbichler commented Feb 20, 2018

Rombur commented Feb 20, 2018

kkormann Feb 12, 2018 •

edited

kkormann Feb 18, 2018 •

edited