Error when USE_FE_Q_iso_Q1 enabled

[vovannikov]

When trying to run this:

mpirun -n 1 ./applications/sintering/sintering-2D --circle 2 ../pf-applications/applications/sintering/analysis_examples/default.prm

I get the following error:

Create subdivided hyperrectangle [37.500000x22.500000] with 2 refinements and 16x10 subdivisions

System statistics:
  - n cell:                    2560 (min: 2560, max: 2560)
  - n cell w hanging nodes:    0 (min: 0, max: 0)
  - n cell wo hanging nodes:   2560 (min: 2560, max: 2560)
  - n levels:                  3
  - n dofs:                    10449

Execute refinement/coarsening:
System statistics:
  - n cell:                    5908 (min: 5908, max: 5908)
  - n cell w hanging nodes:    658 (min: 658, max: 658)
  - n cell wo hanging nodes:   5250 (min: 5250, max: 5250)
  - n levels:                  4
  - n dofs:                    24481
  - n components:              4

Execute refinement/coarsening:
System statistics:
  - n cell:                    5944 (min: 5944, max: 5944)
  - n cell w hanging nodes:    532 (min: 532, max: 532)
  - n cell wo hanging nodes:   5412 (min: 5412, max: 5412)
  - n levels:                  4
  - n dofs:                    24481
  - n components:              4

Execute refinement/coarsening:
System statistics:
  - n cell:                    5944 (min: 5944, max: 5944)
  - n cell w hanging nodes:    532 (min: 532, max: 532)
  - n cell wo hanging nodes:   5412 (min: 5412, max: 5412)
  - n levels:                  4
  - n dofs:                    24481
  - n components:              4

Execute grain tracker:
Grain tracker CPU time = 0.004787 sec
Number of order parameters: 2
Number of grains: 2
op_index_current = 0 | op_index_old = 0 | segments = 1 | grain_index = 0
    segment: center = 7.55077 0 | radius = 9.15176
op_index_current = 1 | op_index_old = 1 | segments = 1 | grain_index = 1
    segment: center = -7.55077 0 | radius = 9.15176
Outputing data at t = 0 (./solution.0.vtu)
Outputing data at t = 0 (./contour_solution.0.vtu)

it      res_abs      res_rel   ch_rel_abs   ch_res_rel   ac_rel_abs   ac_res_rel
 0 3.456635e+00 ------------ 2.251159e-14 ------------ 1.069416e+00 ------------

Create sparsity pattern (cahn_hilliard_op) with:
 - NNZ: 883428

terminate called after throwing an instance of 'dealii::TrilinosWrappers::SparseMatrix::ExcTrilinosError'
  what():
--------------------------------------------------------
An error occurred in line <1681> of file </home/vladimir/development/dealii/source/lac/trilinos_sparse_matrix.cc> in function
    void dealii::TrilinosWrappers::SparseMatrix::add(dealii::TrilinosWrappers::SparseMatrix::size_type, dealii::TrilinosWrappers::SparseMatrix::size_type, const size_type*, const TrilinosScalar*, bool, bool)
The violated condition was:
    ierr == 0
Additional information:
    An error with error number 2 occurred while calling a Trilinos
    function.

    For historical reasons, many Trilinos functions express errors by
    returning specific integer values to indicate certain errors.
    Unfortunately, different Trilinos functions often use the same integer
    values for different kinds of errors, and in most cases it is also not
    documented what each error code actually means. As a consequence, it
    is often difficult to say what a particular error (in this case, the
    error with integer code '2') represents and how one should fix a code
    to avoid it. The best one can often do is to look up the call stack to
    see which deal.II function generated the error, and which Trilinos
    function the error code had originated from; then look up the Trilinos
    source code of that function (for example on github) to see what code
    path set that error code. Short of going through all of that, the only
    other option is to guess the cause of the error from the context in
    which the error appeared.

[peterrum]

Have you updated deal.II? We need dealii/dealii#14197. I have stated to recompile it, since I was on an other branch.

[peterrum]

It's running on my side.

[vovannikov]

Ah, I haven't updated, missed that. Thank you for the hint.

[peterrum]

@vovannikov Any news?

[vovannikov]

Everything works after the update, cool!

But I have a few additional questions.

1. To make FE_Q_iso_Q1 work, I need to have Approximation settings to be defined like this, right?

subsection Approximation
    set FEDegree = 2
    set NPoints1D = 4
    set Subdivisions = 2
end

If I do not alter the first 2 options, then the corresponding asserts will get triggered.

2. If in the section above I set set Subdivisions = 1, then the asserts are still in play but FE_Q is still used, am I right?

3. Maybe we should output somewhere the approximation type and integration rule which is currently used? As it is now, given the above implicit restictions coming from FE_DEGREE and N_Q_POINTS_1D these things are not explicitly clear.

[peterrum]

diff --git a/applications/sintering/sintering.cc b/applications/sintering/sintering.cc
index afbb06e..32df6d2 100644
--- a/applications/sintering/sintering.cc
+++ b/applications/sintering/sintering.cc
@@ -23,7 +23,7 @@ static_assert(false, "No dimension has been given!");
 static_assert(false, "No grains number has been given!");
 #endif
 
-//#define USE_FE_Q_iso_Q1
+#define USE_FE_Q_iso_Q1
 
 #ifdef USE_FE_Q_iso_Q1
 #  define FE_DEGREE 2

diff --git a/applications/sintering/analysis_examples/49particles.json b/applications/sintering/analysis_examples/49particles.json
index ff5bad5..270c806 100644
--- a/applications/sintering/analysis_examples/49particles.json
+++ b/applications/sintering/analysis_examples/49particles.json
@@ -4,13 +4,13 @@
     "Approximation": {
         "FEDegree": "1",
         "NPoints1D": "2",
-        "NSubdivisions": "1"
+        "NSubdivisions": "2"
     },
 }

[peterrum]

For testing purposes, one could disable these lines:

#ifdef USE_FE_Q_iso_Q1
#  define FE_DEGREE 2
#  define N_Q_POINTS_1D FE_DEGREE * 2
#else
#  define FE_DEGREE 1
#  define N_Q_POINTS_1D FE_DEGREE + 1
#endif

In this case, FEEvaluation would take the non-templated path. Which is a bit slower but good enough.

[peterrum]

Maybe we should output somewhere the approximation type and integration rule which is currently used? As it is now, given the above implicit restictions coming from FE_DEGREE and N_Q_POINTS_1D these things are not explicitly clear.

Good idea.