From 6c4b80b739477290930b020eb529982c442fe9b4 Mon Sep 17 00:00:00 2001 From: Jean-Paul Pelteret Date: Wed, 25 Jul 2018 23:13:12 +0200 Subject: [PATCH] AD Helpers: Add helper for residual linearisation (cell-level) --- .../changes/major/20180930Jean-PaulPelteret | 5 + .../deal.II/differentiation/ad/ad_helpers.h | 294 ++++++++++++++++++ source/differentiation/ad/ad_helpers.cc | 144 +++++++++ source/differentiation/ad/ad_helpers.inst1.in | 16 + source/differentiation/ad/ad_helpers.inst2.in | 28 ++ 5 files changed, 487 insertions(+) create mode 100644 doc/news/changes/major/20180930Jean-PaulPelteret diff --git a/doc/news/changes/major/20180930Jean-PaulPelteret b/doc/news/changes/major/20180930Jean-PaulPelteret new file mode 100644 index 0000000000..7b6ec4731c --- /dev/null +++ b/doc/news/changes/major/20180930Jean-PaulPelteret @@ -0,0 +1,5 @@ +New: A new class Differentiation::AD::ADHelperResidualLinearisation has been +added to help compute the linearization of a residual vector defined on the +level of a cell (a finite element residual), or for local nonlinear equations. +
+(Jean-Paul Pelteret, 2018/09/30) diff --git a/include/deal.II/differentiation/ad/ad_helpers.h b/include/deal.II/differentiation/ad/ad_helpers.h index c131f31fdc..75bd125da6 100644 --- a/include/deal.II/differentiation/ad/ad_helpers.h +++ b/include/deal.II/differentiation/ad/ad_helpers.h @@ -1277,6 +1277,300 @@ namespace Differentiation }; // class ADHelperEnergyFunctional + + /** + * A helper class that facilitates the evaluation and automated + * linearization of a vector of functions that represents a residual vector + * (as computed from some corresponding local degree of freedom values). + * This class would typically be used to compute the linearization of a + * residual vector defined on the level of a cell, or for local + * nonlinear equations. + * + * An example of its usage in the case of a residual linearization + * might be as follows (in this case we'll compute the + * linearization of a finite-strain magneto-elastic solid from the residual, + * as constructed from the Piola-Kirchoff stress and magnetic induction + * assuming the magnetic scalar potential formulation): + * + * @code + * // Existing data structures: + * Vector solution (...); // Or another vector type + * std::vector local_dof_indices (...); + * const FEValuesExtractors::Vector u_fe (...); + * const FEValuesExtractors::Scalar msp_fe (...); + * const unsigned int u_block (...); + * const unsigned int msp_block (...); + * FEValues fe_values (...); + * const unsigned int n_q_points (...); + * FullMatrix cell_matrix (...); + * Vector cell_rhs (...); + * + * // Assembly loop: + * for (auto cell & : ...) + * { + * cell->get_dof_indices(local_dof_indices); + * const unsigned int n_independent_variables + * = local_dof_indices.size(); + * const unsigned int n_dependent_variables + * = local_dof_indices.size(); + * + * // Create some aliases for the AD helper. + * // In this example, we strictly assume that we're using tapeless + * // AD types, and so the AD_typecode used in the template argument + * // must refer to one of these types. See the example for the + * // ADHelperEnergyFunctional for details on how to extend + * // support to taped AD numbers. + * using ADHelper = AD::ADHelperEnergyFunctional<...>; + * using ADNumberType = typename ADHelper::ad_type; + * + * // Create and initialize an instance of the helper class. + * ADHelper ad_helper(n_independent_variables,n_dependent_variables); + * + * // Initialize the local data structures for assembly. + * // This is also taken care of by the ADHelper, so this step could + * // be skipped. + * cell_rhs.reinit(n_dependent_variables); + * cell_matrix.reinit(n_independent_variables,n_dependent_variables); + * + * // This next code block corresponds to the "recording" phase, where + * // the operations performed with the AD numbers are tracked and + * // differentiation is performed. + * { + * // First, we set the values for all DoFs. + * ad_helper.register_dof_values(solution, local_dof_indices); + * + * // Then we get the complete set of degree of freedom values as + * // represented by auto-differentiable numbers. The operations + * // performed with these variables are tracked by the AD library + * // from this point until stop_recording_operations() is called. + * const std::vector dof_values_ad + * = ad_helper.get_sensitive_dof_values(); + * + * // Then we do some problem specific tasks, the first being to + * // compute all values, gradients, etc. based on sensitive AD DoF + * // values. Here we are fetching the displacement gradients at each + * // quadrature point, as well as the gradients of the magnetic + * // scalar potential field. + * std::vector> Grad_u( + * n_q_points, Tensor<2, dim, ADNumberType>()); + * std::vector> Grad_msp( + * n_q_points, Tensor<1, dim, ADNumberType>()); + * fe_values[u_fe].get_function_gradients_from_local_dof_values( + * dof_values_ad, Grad_u); + * fe_values[msp_fe].get_function_gradients_from_local_dof_values( + * dof_values_ad, Grad_msp); + * + * // This variable stores the cell residual vector contributions. + * // IMPORTANT: Note that each entry is hand-initialized with a value + * // of zero. This is a highly recommended practise, as some AD + * // numbers appear not to safely initialize their internal data + * // structures. + * std::vector residual_ad ( + * n_dependent_variables, ADNumberType(0.0)); + * + * // Compute the cell total residual + * // = (internal + external) contributions + * for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + * { + * // Calculate the deformation gradient and magnetic field at this + * // quadrature point + * const Tensor<2, dim, ADNumberType> F = + * unit_symmetric_tensor() + Grad_u[q_point]; + * const Tensor<1, dim, ADNumberType> H = -Grad_msp[q_point]; + * Assert(numbers::value_is_greater_than(determinant(F), 0.0), + * ExcMessage("Negative determinant of the deformation " + * "gradient detected!")); + * + * // Extract some configuration dependent variables from our + * // nonlinear constitutive law for the current quadrature point. + * // In this way they only have to be computed once per quadrature + * // point. + * const SymmetricTensor<2,dim,ad_type> S = get_S(F,H); + * const Tensor<1,dim,ad_type> B = get_B(F,H); + * + * // Define some position-dependent aliases, to make the assembly + * // process easier to follow. + * const double JxW = fe_values.JxW(q_point); + * + * // Add contribution of the internal forces: + * // Note that we assemble the residual contribution directly + * // as it is this vector that is to be automatically linearized. + * for (unsigned int I = 0; I < n_dofs_per_cell; ++I) + * { + * // Determine the component and block associated with + * // the I'th local degree of freedom. + * const unsigned int block_I = + * fe.system_to_base_index(I).first.first; + * + * if (block_I == u_block) // u-terms + * { + * // Variation of the Green-Lagrange strain tensor + * // associated with the I'th vector-valued basis function. + * const SymmetricTensor<2,dim,double> dE_I + * = symmetrize(transpose(F) + * * fe_values[u_fe].gradient(I,q_point)); + * + * residual_ad[I] += (dE_I*S) * JxW; + * } + * else if (block_I == msp_block) + * { + * // Variation of the magnetic field vector associated with + * // the I'th scalar-valued basis function + * const Tensor<1,dim,double> dH_I + * = -fe_values[msp_fe].gradient(I, q_point); + * + * residual_ad[I] -= (dH_I*B) * JxW; + * } + * } + * } + * + * // Add contribution from external sources. If these contributions + * // are also solution dependent then they will also be consistently + * // linearized. + * // Loop over faces and accumulate external contributions into the + * // cell total residual. + * for (unsigned int face : ...) + * if (cell->face(face)->at_boundary()) + * residual_ad[I] += ... + * + * // Register the definition of the cell residual + * ad_helper.register_residual_vector(residual_ad); + * } + * + * // Compute the residual values and their Jacobian at the + * // evaluation point + * ad_helper.compute_residual(cell_rhs); + * cell_rhs *= -1.0; // RHS = - residual + * ad_helper.compute_linearization(cell_matrix); + * } + * @endcode + * + * In most use cases, and in particular in the code example shown above, + * both the number of independent and dependent variables equals the + * number of dofs_per_cell for the used finite element. + * + * @warning ADOL-C does not support the standard threading models used by + * deal.II, so this class should @b not be embedded within a multithreaded + * function when using ADOL-C number types. It is, however, suitable for use + * in both serial and MPI routines. + * + * @author Jean-Paul Pelteret, 2016, 2017, 2018 + */ + template + class ADHelperResidualLinearisation + : public ADHelperCellLevelBase + { + public: + /** + * Type definition for the floating point number type that is used in, + * and results from, all computations. + */ + using scalar_type = + typename ADHelperBase::scalar_type; + + /** + * Type definition for the auto-differentiation number type that is used + * in all computations. + */ + using ad_type = + typename ADHelperBase::ad_type; + + /** + * @name Constructor / destructor + */ + //@{ + + /** + * The constructor for the class. + * + * @param[in] n_independent_variables The number of independent variables + * that will be used in the definition of the functions that it is + * desired to compute the sensitivities of. In the computation of + * $\mathbf{r}(\mathbf{X})$, this will be the number of inputs + * $\mathbf{X}$, i.e. the dimension of the domain space. + * @param[in] n_dependent_variables The number of scalar functions to be + * defined that will have a sensitivity to the given independent + * variables. In the computation of $\mathbf{r}(\mathbf{X})$, this will + * be the number of outputs $\mathbf{r}$, i.e. the dimension of the + * image space. + */ + ADHelperResidualLinearisation(const unsigned int n_independent_variables, + const unsigned int n_dependent_variables); + + /** + * Destructor + */ + virtual ~ADHelperResidualLinearisation() = default; + + //@} + + /** + * @name Dependent variables + */ + //@{ + + /** + * Register the definition of the cell residual vector + * $\mathbf{r}(\mathbf{X})$. + * + * @param[in] residual A vector of recorded functions that defines the + * residual. The components of this vector represents the dependent + * variables. + * + * @note For this class that expects only vector fields of dependent + * variables, this function must only be called once per tape. + * + * @note For taped AD numbers, this operation is only valid in recording mode. + */ + void + register_residual_vector(const std::vector &residual); + + /** + * Evaluation of the residual for a chosen set of degree of freedom + * values. This corresponds to the computation of the residual vector, + * i.e. + * @f[ + * \mathbf{r}(\mathbf{X}) \vert_{\mathbf{X}} + * @f] + * + * The values at the evaluation point $\mathbf{X}$ are obtained by calling + * ADHelperCellLevelBase::set_dof_values(). + * + * @param[out] residual A Vector object, for which the value for each + * entry represents the residual value for the corresponding local + * degree of freedom. The output @p residual vector has a length + * corresponding to @p n_dependent_variables. + */ + virtual void + compute_residual(Vector &residual) const override; + + /** + * Computes the linearization of the residual vector around a chosen set + * of degree of freedom values. Underlying this is the computation of the + * gradient (first derivative) of the residual vector $\mathbf{r}$ with + * respect to all independent variables, i.e. + * @f[ + * \frac{\partial\mathbf{r}(\mathbf{X})}{\partial\mathbf{X}} + * @f] + * + * The values at the evaluation point $\mathbf{X}$ are obtained by calling + * ADHelperCellLevelBase::set_dof_values(). + * + * @param[out] linearization A FullMatrix representing the linearization + * of the residual vector. The output @p linearization matrix has + * dimensions corresponding to + * n_dependent_variables$\times$n_independent_variables. + */ + virtual void + compute_linearization( + FullMatrix &linearization) const override; + + //@} + + }; // class ADHelperResidualLinearisation + + } // namespace AD } // namespace Differentiation diff --git a/source/differentiation/ad/ad_helpers.cc b/source/differentiation/ad/ad_helpers.cc index ea732e626a..844b7653f7 100644 --- a/source/differentiation/ad/ad_helpers.cc +++ b/source/differentiation/ad/ad_helpers.cc @@ -977,6 +977,150 @@ namespace Differentiation } + /* ------------------- ADHelperResidualLinearisation ------------------- */ + + + + template + ADHelperResidualLinearisation:: + ADHelperResidualLinearisation(const unsigned int n_independent_variables, + const unsigned int n_dependent_variables) + : ADHelperCellLevelBase( + n_independent_variables, + n_dependent_variables) + {} + + + + template + void + ADHelperResidualLinearisation:: + register_residual_vector(const std::vector &residual) + { + Assert(residual.size() == this->n_dependent_variables(), + ExcMessage( + "Vector size does not match number of dependent variables")); + for (unsigned int i = 0; i < this->n_dependent_variables(); ++i) + ADHelperBase::register_dependent_variable( + i, residual[i]); + } + + + + template + void + ADHelperResidualLinearisation:: + compute_residual(Vector &values) const + { + if ((ADNumberTraits::is_taped == true && + this->taped_driver.keep_independent_values() == false) || + ADNumberTraits::is_tapeless == true) + { + Assert( + this->n_registered_independent_variables() == + this->n_independent_variables(), + ExcMessage( + "Not all values of sensitivities have been registered or subsequently set!")); + } + Assert(this->n_registered_dependent_variables() == + this->n_dependent_variables(), + ExcMessage("Not all dependent variables have been registered.")); + + // We can neglect correctly initializing the entries as + // we'll be overwriting them immediately. + if (values.size() != this->n_dependent_variables()) + values.reinit(this->n_dependent_variables(), + true /*omit_zeroing_entries*/); + + if (ADNumberTraits::is_taped == true) + { + Assert(this->active_tape_index() != + Numbers::invalid_tape_index, + ExcMessage("Invalid tape index")); + Assert(this->is_recording() == false, + ExcMessage( + "Cannot compute values while tape is being recorded.")); + Assert(this->independent_variable_values.size() == + this->n_independent_variables(), + ExcDimensionMismatch(this->independent_variable_values.size(), + this->n_independent_variables())); + + TapedDrivers::values( + this->active_tape_index(), + this->n_dependent_variables(), + this->independent_variable_values, + values); + } + else + { + Assert(ADNumberTraits::is_tapeless == true, + ExcInternalError()); + TapelessDrivers::values( + this->dependent_variables, values); + } + } + + + + template + void + ADHelperResidualLinearisation:: + compute_linearization(FullMatrix &jacobian) const + { + if ((ADNumberTraits::is_taped == true && + this->taped_driver.keep_independent_values() == false) || + ADNumberTraits::is_tapeless == true) + { + Assert( + this->n_registered_independent_variables() == + this->n_independent_variables(), + ExcMessage( + "Not all values of sensitivities have been registered or subsequently set!")); + } + Assert(this->n_registered_dependent_variables() == + this->n_dependent_variables(), + ExcMessage("Not all dependent variables have been registered.")); + + // We can neglect correctly initializing the entries as + // we'll be overwriting them immediately. + if (jacobian.m() != this->n_dependent_variables() || + jacobian.n() != this->n_independent_variables()) + jacobian.reinit({this->n_dependent_variables(), + this->n_independent_variables()}, + true /*omit_default_initialization*/); + + if (ADNumberTraits::is_taped == true) + { + Assert(this->active_tape_index() != + Numbers::invalid_tape_index, + ExcMessage("Invalid tape index")); + Assert(this->is_recording() == false, + ExcMessage( + "Cannot compute hessian while tape is being recorded.")); + Assert(this->independent_variable_values.size() == + this->n_independent_variables(), + ExcDimensionMismatch(this->independent_variable_values.size(), + this->n_independent_variables())); + + TapedDrivers::jacobian( + this->active_tape_index(), + this->n_dependent_variables(), + this->independent_variable_values, + jacobian); + } + else + { + Assert(ADNumberTraits::is_tapeless == true, + ExcInternalError()); + Assert(this->independent_variables.size() == + this->n_independent_variables(), + ExcInternalError()); + TapelessDrivers::jacobian( + this->independent_variables, this->dependent_variables, jacobian); + } + } + + } // namespace AD } // namespace Differentiation diff --git a/source/differentiation/ad/ad_helpers.inst1.in b/source/differentiation/ad/ad_helpers.inst1.in index 418cdf3b78..79bc953f1c 100644 --- a/source/differentiation/ad/ad_helpers.inst1.in +++ b/source/differentiation/ad/ad_helpers.inst1.in @@ -45,6 +45,14 @@ for (number : REAL_SCALARS) template class ADHelperEnergyFunctional; + // -------------------------- ADHelperResidualLinearisation ---------------------- + + template + class ADHelperResidualLinearisation; + + template + class ADHelperResidualLinearisation; + \} \} } @@ -81,6 +89,14 @@ for () template class ADHelperEnergyFunctional::ad_type>; + // -------------------------- ADHelperResidualLinearisation ---------------------- + + template + class ADHelperResidualLinearisation::ad_type>; + + template + class ADHelperResidualLinearisation::ad_type>; + \} \} diff --git a/source/differentiation/ad/ad_helpers.inst2.in b/source/differentiation/ad/ad_helpers.inst2.in index be1bc98baa..1e64944e57 100644 --- a/source/differentiation/ad/ad_helpers.inst2.in +++ b/source/differentiation/ad/ad_helpers.inst2.in @@ -63,6 +63,20 @@ for (number : REAL_SCALARS) template class ADHelperEnergyFunctional; + // -------------------------- ADHelperResidualLinearisation ---------------------- + + template + class ADHelperResidualLinearisation; + + template + class ADHelperResidualLinearisation; + + template + class ADHelperResidualLinearisation; + + template + class ADHelperResidualLinearisation; + \} \} } @@ -117,6 +131,20 @@ for () template class ADHelperEnergyFunctional::ad_type>; + // -------------------------- ADHelperResidualLinearisation ---------------------- + + template + class ADHelperResidualLinearisation::ad_type>; + + template + class ADHelperResidualLinearisation::ad_type>; + + template + class ADHelperResidualLinearisation::ad_type>; + + template + class ADHelperResidualLinearisation::ad_type>; + \} \} -- 2.39.5