From 5e58dbdcfabace40665e20ff97818cad4c58d594 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Fri, 30 Sep 2016 16:56:35 +0200 Subject: [PATCH] Augment documentation of FEEvaluation --- include/deal.II/matrix_free/fe_evaluation.h | 373 +++++++++++++++++--- include/deal.II/matrix_free/matrix_free.h | 2 + 2 files changed, 322 insertions(+), 53 deletions(-) diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h index d98bc41696..f7f5302966 100644 --- a/include/deal.II/matrix_free/fe_evaluation.h +++ b/include/deal.II/matrix_free/fe_evaluation.h @@ -1376,9 +1376,9 @@ protected: /** * The class that provides all functions necessary to evaluate functions at * quadrature points and cell integrations. In functionality, this class is - * similar to FEValues, however, it includes a lot of specialized - * functions that make it much faster (between 5 and 500, depending on the - * polynomial order). + * similar to FEValues, however, it includes a lot of specialized functions + * that make it much faster (between 5 and 500, depending on the polynomial + * degree). * *

Usage and initialization

* @@ -1387,18 +1387,121 @@ protected: * The first and foremost way of usage is to initialize this class from a * MatrixFree object that caches everything related to the degrees of freedom * and the mapping information. This way, it is possible to use vectorization - * for applying a vector operation for several cells at once. This setting is - * explained in the step-37 and step-48 tutorial programs. For vector-valued - * problems, the deal.II test suite includes a few additional examples as - * well, e.g. the Stokes operator found at - * https://github.com/dealii/dealii/blob/master/tests/matrix_free/matrix_vector_stokes_noflux.cc + * for applying a differential operator for several cells at once. + * + * The capabilities of FEEvaluation span a large spectrum of integration tasks + * for weak forms. In general, there are two classes of tasks that get + * done. One is the @p evaluate path that interpolates from a solution vector + * to quadrature points: + * + * @code + * FEEvaluation phi(matrix_free); + * for (unsigned int cell_index = cell_range.first; + * cell_index < cell_range.second; ++cell_index) + * { + * phi.reinit(cell_index); + * phi.read_dof_values(vector); + * phi.evaluate(true, false); // interpolate values, but not gradients + * for (unsigned int q_index=0; q_index val = phi.get_value(q_index); + * // do something with val + * } + * } + * @endcode * - * For most operator evaluation tasks, this path provides the most efficient - * solution by combining pre-computed data for the mapping (Jacobian - * transformations for the geometry description) with on-the-fly evaluation of - * basis functions. In other words, the framework provides a trade-off between - * memory usage and initialization of objects that is suitable for matrix-free - * operator evaluation. + * Likewise, a gradient of the finite element solution represented by @p + * vector can be interpolated to the quadrature points by @p + * phi.get_gradient(q). The combination of read_dof_values(), evaluate() and + * get_value() is similar to what FEValues::get_function_values or + * FEValues::get_function_gradients does, but it is in general much faster + * because it makes use of the tensor product, see the description of the + * evaluation routines below, and can do this operation for several cells at + * once through vectorization. + * + * The second class of tasks done by FEEvaluation are integration tasks for + * right hand sides. In finite element computations, these typically consist + * of multiplying a quantity on quadrature points (a function value, or a + * field interpolated by the finite element space itself) by a set of test + * functions and integrating over the cell through summation of the values in + * each quadrature point, multiplied by the quadrature weight and the Jacobian + * determinant of the transformation. If a generic Function object is given + * and we want to compute $v_i = \int_\Omega \varphi_i f dx$, this is done by + * the following cell-wise integration: + * + * @code + * FEEvaluation phi(matrix_free); + * Function &function = ...; + * for (unsigned int cell_index = cell_range.first; + * cell_index < cell_range.second; ++cell_index) + * { + * phi.reinit(cell_index); + * for (unsigned int q_index=0; q_index > p_vect = + * phi.quadrature_point(q_index); + * // Need to evaluate function for each component in VectorizedArray + * VectorizedArray f_value; + * for (unsigned int v=0; v::n_array_elements; ++v) + * { + * Point p; + * for (unsigned int d=0; d(1) in 2D) is computed, which gives the vector entries to be + * summed into the @p dst vector. Note that the above code needs to explicitly + * loop over the components in the vectorized array for evaluating the + * function, which is necessary for interfacing with a generic Function object + * with double arguments. Simple functions can also be implemented in + * VectorizedArray form directly as VectorizedArray provides the basic math + * operations. + * + * For evaluating a bilinear form, the evaluation on a source vector is + * combined with the integration involving test functions that get written + * into a result vector. This setting is the context of matrix-free operator + * evaluation and explained in the step-37 and step-48 tutorial programs. + * + * Note that the two vector accesses through FEEvaluation::read_dof_values and + * FEEvaluation::distribute_local_to_global resolve constraints on the fly, + * based on the ConstraintMatrix object specified at the MatrixFree::reinit() + * call. In case the values in the degrees of freedom are of interest (usually + * only the values in quadrature points are necessary), these can be accessed + * through FEEvaluation::get_dof_value(i), where i is the index of the basis + * function. Note that the numbering of the degrees of freedom for continuous + * elements in FEEvaluation is different from the ordering in FE_Q (or + * FEValues) because FEEvaluation needs to access them in lexicographic order, + * which is the ordering used in FE_DGQ, for instance. Re-indexing would be + * too expensive because the access inside evaluate() and integrate() is on + * the critical path in the tensorial evaluation parts. An alternative to + * filling the DoF values by read_dof_values() before an evaluate() call is to + * manually assign a value by a set_dof_value() call. Likewise, if the local + * result of integration should be further processed rather than scattered + * into a vector by distribute_local_to_global(), one can access it by + * get_dof_value() after an integrate() call. An example for using the values + * of an integral in a different context is fast assembly of matrices as shown + * in the next subsection. + * + * For most operator evaluation tasks that repeatedly go through the mesh, the + * realization by MatrixFree that combines pre-computed data for the mapping + * (Jacobian transformations for the geometry description) with on-the-fly + * evaluation of basis functions is the most efficient way of doing things. In + * other words, the framework selects a trade-off between memory usage and + * initialization of objects that is suitable for replacement of matrix-vector + * products or explicit time integration in a matrix-free way. * *

Usage without pre-initialized MatrixFree object

* @@ -1412,11 +1515,11 @@ protected: * MatrixFree class may not be desired. In such a case, the cost to initialize * the necessary geometry data on the fly is comparably low and thus avoiding * a global object MatrixFree can be useful. When used in this way, reinit - * methods reminiscent from FEValues with a cell iterator are to be used. - * However, note that this model results in working on a single cell at a - * time, with geometry data duplicated in all components of the vectorized - * array. Thus, vectorization is only useful when it can apply the same - * operation on different data, e.g. when performing matrix assembly. + * methods reminiscent from FEValues with a cell iterator are used. However, + * note that this model results in working on a single cell at a time, with + * geometry data duplicated in all components of the vectorized array. Thus, + * vectorization is only useful when it can apply the same operation on + * different data, e.g. when performing matrix assembly. * * As an example, consider the following code to assemble the contributions to * the Laplace matrix: @@ -1438,11 +1541,16 @@ protected: * * // Set n_items unit vectors * for (unsigned int j=0; j(); + * fe_eval.set_dof_value(VectorizedArray(), j); * for (unsigned int v=0; v one_value = VectorizedArray(); + * one_value[v] = 1.; + * fe_eval.set_dof_value(one_value, i+v); + * } * - * // Apply operator on unit vector + * // Apply operator on unit vector to generate the next few matrix + * // columns * fe_eval.evaluate(true, true); * for (unsigned int q=0; q(fe_degree+1)2*dim+1 rather than + * (fe_degree+1)3*dim . * - * Due to vectorization, we can actually generate matrix data for several unit + * Due to vectorization, we can generate matrix columns for several unit * vectors at a time (e.g. 4). The variable @p n_items make sure that we do * the last iteration where the number of cell dofs is not divisible by the * vectorization length correctly. Also note that we need to get the internal * dof numbering applied by fe_eval because FEEvaluation internally uses a - * lexicographic numbering of degrees of freedom. This is necessary to - * efficiently work with tensor products where all degrees of freedom along a - * dimension must be laid out in a regular way. + * lexicographic numbering of degrees of freedom as explained above. * *

Description of evaluation routines

* - * This class contains specialized evaluation routines for several elements - * based on tensor-product quadrature formulas and tensor-product-like shape - * functions, including standard FE_Q or FE_DGQ elements and quadrature points - * symmetric around 0.5 (like Gauss quadrature), FE_DGP elements based on - * truncated tensor products as well as the faster case of Gauss-Lobatto - * elements with Gauss-Lobatto quadrature which give diagonal mass matrices - * and quicker evaluation internally. The main benefit of this class is the - * evaluation of all shape functions in all quadrature or integration over all - * shape functions in dim (fe_degree+1)dim+1 - * operations instead of the slower - * (fe_degree+1)2*dim complexity in the evaluation routines - * of FEValues. + * This class contains specialized evaluation routines for elements based on + * tensor-product quadrature formulas and tensor-product-like shape functions, + * including standard FE_Q or FE_DGQ elements and quadrature points symmetric + * around 0.5 (like Gauss quadrature), FE_DGP elements based on truncated + * tensor products as well as the faster case of Gauss-Lobatto elements with + * Gauss-Lobatto quadrature which give diagonal mass matrices and quicker + * evaluation internally. The main benefit of this class is the evaluation of + * all shape functions in all quadrature or integration over all shape + * functions in dim (fe_degree+1)dim+1 operations + * instead of the slower (fe_degree+1)2*dim + * complexity in the evaluation routines of FEValues. This is done by an + * algorithm called sum factorization which factors out constant factors + * during the evaluation along a coordinate direction. This algorithm is the + * basis of many spectral element algorithms. * * Note that many of the operations available through this class are inherited * from the base class FEEvaluationBase, in particular reading from and @@ -1500,27 +1610,184 @@ protected: * implements access to values, gradients and Hessians of the finite element * function on quadrature points. * - * This class assumes that shape functions of the FiniteElement under - * consideration do not depend on the actual shape of the cells in - * real space. Currently, other finite elements cannot be treated with the + * This class assumes that the shape functions of the FiniteElement under + * consideration do not depend on the geometry of the cells in real + * space. Currently, other finite elements cannot be treated with the * matrix-free concept. * - * This class has five template arguments: + *

Handling multi-component systems

* - * @param dim Dimension in which this class is to be used + * FEEvaluation also allows for treating vector-valued problems through a + * template parameter on the number of components: + * + * @code + * FEEvaluation phi(matrix_free); + * @endcode + * + * If used this way, the components can be gathered from several components of + * an @p std::vector through the call + * + * @code + * phi.read_dof_values(src, 0); + * @endcode * - * @param fe_degree Degree of the tensor product finite element with + * where the 0 means that the vectors starting from the zeroth vector in the + * @p std::vector should be used, src[0], src[1], ..., + * src[n_components-1]. + * + * An alternative way for reading multi-component systems is possible if the + * DoFHandler underlying the MatrixFree data is based on an FESystem of @p + * n_components entries. In that case, a single vector is provided for the + * read_dof_values() and distribute_local_to_global() calls. + * + * An important property of FEEvaluation in multi-component systems is the + * layout of multiple components in the get_value(), get_gradient(), or + * get_dof_value() calls. In this case, instead of a scalar return field + * VectorizedArray@ a tensor is returned, + * + * @code + * get_value -> Tensor<1,n_components,VectorizedArray > + * get_gradient -> Tensor<1,n_components,Tensor<1,dim,VectorizedArray > > + * @endcode + * + * In a similar vein, the submit_value() and submit_gradient() calls take + * tensors of values. Note that there exist specializations for @p + * n_components=1 and @p n_components=dim, which are provided through the base + * class FEEvaluationAccess. In the scalar case, these provide the scalar + * return types described above. In the vector-valued case, the gradient is + * converted from Tensor@<1,dim,Tensor@<1,dim,VectorizedArray@ + * @> @> to Tensor@<2,dim,VectorizedArray@ + * @>. Furthermore, additional operations such as the diveregence or + * curl are available. + * + * In case different shape functions are combined, for example mixed finite + * element formulations in Stokes flow, two FEEvaluation objects are created, + * one for the velocity and one for the pressure. Those are then combined on + * quadrature points: + * + * @code + * FEEvaluation velocity (data, 0); + * FEEvaluation pressure (data, 1); + * + * for (unsigned int cell=cell_range.first; cell > sym_grad_u = + * velocity.get_symmetric_gradient (q); + * VectorizedArray pres = pressure.get_value(q); + * VectorizedArray div = -trace(sym_grad_u); + * pressure.submit_value (div, q); + * + * // subtract p * I + * for (unsigned int d=0; dHandling several integration tasks and data storage + * + * The design of FEEvaluation and MatrixFree separates the geometry from the + * basis functions. Therefore, several DoFHandler objects (or the same + * DoFHandler equipped with different constraint objects) can share the same + * geometry information like in the Stokes example above. All geometry is + * cached once in MatrixFree, so FEEvaluation does not need to do expensive + * initialization calls and rather sets a few pointers. This realization is + * based on the idea that the geometry information is needed only once also + * when several fields are evaluated, in a departure from FEValues which sets + * up the internal mapping data for each field. If for example a + * multi-component PDE involves the shape values on one component and the + * shape gradient on the other, no efficiency is lost if both are based on the + * same MatrixFree object where the update flags specify that both @p + * update_values , @p update_gradients , and @p update_jxw_values are + * given. The selection of desired quantities of shape values is through the + * flags in the evaluate() or integrate calls and the access at quadrature + * points: + * + * @code + * phi1.evaluate(true, false); + * phi2.evaluate(false, true); + * for (unsigned int q_index=0; q_index val1 = phi1.get_value(q); + * Tensor<1,dim,VectorizedArray > grad2 = phi2.get_gradient(q); + * Point > point = phi1.get_quadrature_point(q); + * // ... some complicated formula combining those three... + * } + * @endcode + * + * In the loop over quadrature points, one can ask any of the two FEEvaluation + * objects — it does not really matter which one because they only keep + * pointers to the quadrature point data — to provide the quadrature + * point location. + * + * This observation also translates to the case when different differential + * operators are implemented in a program, for example the action of a mass + * matrix for one phase of the algorithm and the action of a stiffness matrix + * in another one. Only a single MatrixFree object is necessary, maintaining + * full efficiency by using different local functions with the respective + * implementation in separate FEEvaluation objects. In other words, a user + * does not need to bother about being conservative when providing + * update_flags to the initialization of MatrixFree for efficiency reasons - + * no overhead incurs inside FEEvaluation, except for at most one or two more + * @p if statements inside the FEEvaluation::reinit() call. Rather, the + * largest set of flags necessary among all calls is perfectly fine from an + * efficiency point of view. + * + * For the combination of different fields, including different solution + * vectors that come from different time steps, it is mandatory that all + * FEEvaluation objects share the same MatrixFree object. This is because the + * way cells are looped by MatrixFree::cell_loop() can be different for + * different DoFHandler or ConstraintMatrix arguments. More precisely, even + * though the layout is going to be the same in serial, there is no guarantee + * about the ordering for different DoFHandler/Constraints in the MPI + * case. The reason is that the algorithm detects cells that need data + * exchange with MPI and those can change for different elements — FE_Q + * with hanging node constraints connects to more neighbors than a FE_DGQ + * element, for instance, and cells which need data exchange are put in + * different positions inside the cell loop. Of course, if the exact same + * DoFHandler, ConstraintMatrix, and options (such as the setting for thread + * parallelism) are set, then the order is going to be the same because the + * algorithm is deterministic. + * + * @tparam dim Dimension in which this class is to be used + * + * @tparam fe_degree Degree of the tensor product finite element with * fe_degree+1 degrees of freedom per coordinate direction * - * @param n_q_points_1d Number of points in the quadrature formula in 1D, + * @tparam n_q_points_1d Number of points in the quadrature formula in 1D, * defaults to fe_degree+1 * - * @param n_components Number of vector components when solving a system of + * @tparam n_components Number of vector components when solving a system of * PDEs. If the same operation is applied to several components of a PDE (e.g. * a vector Laplace equation), they can be applied simultaneously with one * call (and often more efficiently). Defaults to 1. * - * @param Number Number format, usually @p double or @p float. Defaults to @p + * @tparam Number Number format, usually @p double or @p float. Defaults to @p * double * * @author Katharina Kormann and Martin Kronbichler, 2010, 2011 diff --git a/include/deal.II/matrix_free/matrix_free.h b/include/deal.II/matrix_free/matrix_free.h index 04f8ca4dba..34c56a5c12 100644 --- a/include/deal.II/matrix_free/matrix_free.h +++ b/include/deal.II/matrix_free/matrix_free.h @@ -97,6 +97,8 @@ DEAL_II_NAMESPACE_OPEN * operations for several cells with one CPU instruction and is one of the * main features of this framework. * + * For details on usage of this class, see the description of FEEvaluation. + * * @author Katharina Kormann, Martin Kronbichler, 2010, 2011 */ -- 2.39.5