Robustness of FEEvaluation's pre-compiled evaluation routines

[kronbichler]

In exadg/exadg#241, we discussed the following problem:

• The ExaDG code wants more polynomial degrees inside the FEEvaluation evaluators to be pre-compiled with high efficiency (polynomial degree and number of quadrature points are compile-time constants), so it compiles deal.II's file https://github.com/dealii/dealii/blob/master/include/deal.II/matrix_free/evaluation_template_factory.templates.h with a different setting of

  dealii/include/deal.II/matrix_free/evaluation_template_factory_internal.h

  Lines 23 to 25 in dfde63a

  	#ifndef FE_EVAL_FACTORY_DEGREE_MAX
  	# define FE_EVAL_FACTORY_DEGREE_MAX 6
  	#endif

• The duplicate symbols within libdeal_II.so and libexadg.so are resolved differently between different code components. This resulted in run time errors, because the check

  dealii/include/deal.II/matrix_free/fe_evaluation.h

  Lines 8643 to 8645 in cfebfb5

  	if (this->data->data.front().fe_degree > 0 &&
  	fast_evaluation_supported(this->data->data.front().fe_degree,
  	this->data->data.front().n_q_points_1d) &&

  which goes to

  dealii/include/deal.II/matrix_free/evaluation_template_factory.templates.h

  Lines 32 to 40 in dfde63a

  	struct FastEvaluationSupported
  	{
  	template <int fe_degree, int n_q_points_1d>
  	static bool
  	run()
  	{
  	return fe_degree != -1;
  	}
  	};

  went to code in libexadg.so, whereas the code in the same function just a few lines below,

  dealii/include/deal.II/matrix_free/fe_evaluation.h

  Lines 8673 to 8682 in cfebfb5

  	internal::FEFaceEvaluationGatherFactory<
  	dim,
  	typename VectorType::value_type,
  	VectorizedArrayType>::evaluate(n_components,
  	evaluation_flag,
  	internal::get_beginning<
  	typename VectorType::value_type>(
  	input_vector),
  	shared_vector_data,
  	*this);

  called into the code of libdeal_II.so. This results in contradictions between the assumptions made at these two places, creating run-time errors.
• Interestingly, this worked well up to a few months ago. In the ExaDG investigation, nobody made a bisection for the cause, but I am pretty sure that the cause is Matrixfree evaluation with FE_RaviartThomasNodal #13591, where the code in the two said functions was split into different compilation units, whereas it was in the same compilation unit before.
• I know too little about how to ensure priorities in terms of the library selection when two shared libraries, libdeal_II.so and libexadg.so, provide the same symbols (as we really need to hook into a function name we know in the deal.II headers), except for LD_PRELOAD. However, we still recommend this in our documentation

  dealii/include/deal.II/matrix_free/fe_evaluation.h

  Lines 1703 to 1715 in dfde63a

  	* <h4>Pre-compiling code for more polynomial degrees</h4>
  	*
  	* It is also possible to pre-compile the code in FEEvaluation for a different
  	* maximal polynomial degree. This is controlled by the class
  	* internal::FEEvaluationFactory and the implementation in
  	* `include/deal.II/matrix_free/evaluation_template_factory.templates.h`. By
  	* setting the macro `FE_EVAL_FACTORY_DEGREE_MAX` to the desired integer and
  	* instantiating the classes FEEvaluationFactory and FEFaceEvaluationFactory
  	* (the latter for FEFaceEvaluation) creates paths to templated functions for
  	* a possibly larger set of degrees. You can check if fast
  	* evaluation/integration for a given degree/n_quadrature_points pair by calling
  	* FEEvaluation::fast_evaluation_supported() or
  	* FEFaceEvaluation::fast_evaluation_supported().

This issue is here to discuss solutions, also given that my knowledge on that end is rather limited. I think that it is good to have some polynomial degrees compiled in deal.II already as we make use of them in various places, and probably we get more. At the same time, we would like to enable the user to compile more degrees or combinations because that can be critical for an application (it is for ExaDG, for instance). So how to proceed:

• Unless we find someone with knowledge how to ensure the right selection, we might need to expose the maximal degree that gets compiled at the configure stage of deal.II. We do this already through the respective define flag

  dealii/.github/workflows/windows.yml

  Line 35 in dfde63a

  cmake -DCMAKE_BUILD_TYPE=Debug -DCMAKE_INSTALL_PREFIX=c:/project -DDEAL_II_WITH_ZLIB=off -DDEAL_II_CXX_FLAGS="-WX /D FE_EVAL_FACTORY_DEGREE_MAX=2" -T host=x64 -A x64 ..

  but pull through a setting as a define flag to the compile line seems a bit of a hack. Do we want to enable this via a configure option?
• I think we might still allow a user code to compile something separate in the evaluation template factory, but we should try hard to not expose any inconsistency. Referring to the code in the beginning of this post, I believe we should make sure that the decision whether some setting is supported in the pre-compiled degree happen in the same function call as where we rely on this information. This won't change the risks a user takes when mixing the same symbols we provide from libdeal_II.so, but at least we do so consistently.

[kronbichler]

I will try to address this issue in the next week.