From 540ae7e5aaa0fa8a234ce17ae15d88ad5a07fd2f Mon Sep 17 00:00:00 2001 From: Jean-Paul Pelteret Date: Wed, 25 Jul 2018 23:12:29 +0200 Subject: [PATCH] AD Helpers: Add helper for variational formulations (cell-level) --- .../changes/major/20180827Jean-PaulPelteret | 6 + .../deal.II/differentiation/ad/ad_helpers.h | 293 ++++++++++++++++++ source/differentiation/ad/ad_helpers.cc | 202 ++++++++++++ source/differentiation/ad/ad_helpers.inst1.in | 16 + source/differentiation/ad/ad_helpers.inst2.in | 28 ++ 5 files changed, 545 insertions(+) create mode 100644 doc/news/changes/major/20180827Jean-PaulPelteret diff --git a/doc/news/changes/major/20180827Jean-PaulPelteret b/doc/news/changes/major/20180827Jean-PaulPelteret new file mode 100644 index 0000000000..d0476e4b3b --- /dev/null +++ b/doc/news/changes/major/20180827Jean-PaulPelteret @@ -0,0 +1,6 @@ +New: A new class Differentiation::AD::ADHelperEnergyFunctional has been +added to help implement (incremental) variational formulations using automatic +differentiation. In particular, this class is designed to compute the finite +element residuals and their linearizations. +
+(Jean-Paul Pelteret, 2018/08/27) diff --git a/include/deal.II/differentiation/ad/ad_helpers.h b/include/deal.II/differentiation/ad/ad_helpers.h index b6fa81fd42..71f09f9069 100644 --- a/include/deal.II/differentiation/ad/ad_helpers.h +++ b/include/deal.II/differentiation/ad/ad_helpers.h @@ -1031,6 +1031,299 @@ namespace Differentiation }; // class ADHelperCellLevelBase + + /** + * A helper class that facilitates the implementation of a generic + * (incremental) variational formulation from which the computation of the + * residual vector, as well as its linearization, is automated. This class + * would typically be used to derive the residual vector and tangent matrix + * (defined on the level of a cell), or a linearized system of + * equations, starting from a scalar energy functional. + * + * An example of its usage in the case of a residual and tangent + * computations might be as follows (in this case we'll compute the + * linearization of a finite-strain hyperelastic solid from a stored/strain + * energy density function): + * + * @code + * // Existing data structures: + * Vector solution (...); // Or another vector type + * std::vector local_dof_indices (...); + * const FEValuesExtractors::Vector u_fe (...); + * 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(); + * + * // Create some aliases for the AD helper. + * // In the example, the AD_typecode used for the template argument can + * // be refer to either a taped or tapeless type. + * 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); + * + * // 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_independent_variables); + * cell_matrix.reinit(n_independent_variables,n_independent_variables); + * + * // An optional call to set the amount of memory to be allocated to + * // storing taped data. + * // If using a taped AD number then we would likely want to increase + * // the buffer size from the default values as the expression for each + * // residual component will likely be lengthy, and therefore memory + * // intensive. + * ad_helper.set_tape_buffer_sizes(...); + * + * // If using a taped AD number, then at this point we would initiate + * // taping of the expression for the energy for this FE type and + * // material combination: + * + * // Select a tape number to record to + * const types::tape_index tape_index = ...; + * + * // Indicate that we are about to start tracing the operations for + * // function evaluation on the tape. If this tape has already been + * // used (i.e. the operations are already recorded) then we + * // (optionally) load the tape and reuse this data. + * const bool is_recording + * = ad_helper.start_recording_operations(tape_index); + * + * // The steps that follow in the recording phase are required for + * // tapeless methods as well. + * if (is_recording == true) + * { + * // This is the "recording" phase of the operations. + * // 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. + * std::vector> Grad_u( + * n_q_points, Tensor<2, dim, ADNumberType>()); + * fe_values[u_fe].get_function_gradients_from_local_dof_values( + * dof_values_ad, Grad_u); + * + * // This variable stores the cell total energy. + * // IMPORTANT: Note that it 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. + * ADNumberType energy_ad = ADNumberType(0.0); + * + * // Compute the cell total energy = (internal + external) energies + * for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + * { + * // Calculate the deformation gradient at this quadrature point + * const Tensor<2, dim, ADNumberType> F = + * unit_symmetric_tensor() + Grad_u[q_point]; + * Assert(numbers::value_is_greater_than(determinant(F), 0.0), + * ExcMessage("Negative jacobian detected!")); + * + * // Add contribution of the internal energy: + * // Integrate the stored energy density function with the current + * // solution. + * energy_ad += get_Psi(F) * fe_values.JxW(q_point); + * } + * + * // Add contribution from external energy: + * // Loop over faces and accumulate external energy into cell + * // total energy + * // energy_ad += ... + * + * // Register the definition of the total cell energy + * ad_helper.register_energy_functional(energy_ad); + * + * // Indicate that we have completed tracing the operations onto + * // the tape. + * ad_helper.stop_recording_operations(false); // write_tapes_to_file + * } + * else + * { + * // This is the "tape reuse" phase of the operations. + * // Here we will leverage the already traced operations that reside + * // on a tape, and simply re-evaluate the tape at a different point + * // to get the function values and their derivatives. + * + * // Load the existing tape to be reused + * ad_helper.activate_recorded_tape(tape_index); + * + * // Set the new values of the independent variables where the + * // recorded dependent functions are to be evaluated (and + * // differentiated around). + * ad_helper.set_dof_values(solution, local_dof_indices); + * } + * + * // 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 + * + * @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 ADHelperEnergyFunctional + : 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 + * $\Psi(\mathbf{X})$, this will be the number of inputs + * $\mathbf{X}$, i.e. the dimension of the domain space. + * + * @note There is only one dependent variable associated with the total + * energy attributed to the local finite element. That is to say, this + * class assumes that the (local) right hand side and matrix contribution + * is computed from the first and second derivatives of a scalar + * function $\Psi(\mathbf{X})$. + */ + ADHelperEnergyFunctional(const unsigned int n_independent_variables); + + /** + * Destructor + */ + virtual ~ADHelperEnergyFunctional() = default; + + //@} + + /** + * @name Dependent variables + */ + //@{ + + /** + * Register the definition of the total cell energy + * $\Psi(\mathbf{X})$. + * + * @param[in] energy A recorded function that defines the total cell + * energy. This represents the single dependent variable from which both + * the residual and its linearization are to be computed. + * + * @note For this class that expects only a single scalar dependent + * variable, this function must only be called once per tape. + * + * @note For taped AD numbers, this operation is only valid in recording mode. + */ + void + register_energy_functional(const ad_type &energy); + + /** + * Evaluation of the total scalar energy functional for a chosen set of + * degree-of-freedom values, i.e. + * @f[ + * \Psi(\mathbf{X}) \vert_{\mathbf{X}} + * @f] + * + * The values at the evaluation point $\mathbf{X}$ are by calling + * ADHelperCellLevelBase::set_dof_values(). + * + * @return The value of the energy functional at the evaluation point + * corresponding to a chosen set of local degree-of freedom values. + */ + scalar_type + compute_energy() const; + + /** + * Evaluation of the residual for a chosen set of degree-of-freedom + * values. Underlying this is the computation of the gradient (first + * derivative) of the scalar function $\Psi$ with respect to all + * independent variables, i.e. + * @f[ + * \mathbf{r}(\mathbf{X}) = + * \frac{\partial\Psi(\mathbf{X})}{\partial\mathbf{X}} + * \Big\vert_{\mathbf{X}} + * @f] + * + * The values at the evaluation point $\mathbf{X}$ are 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_independent_variables. + */ + 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 + * Hessian (second derivative) of the scalar function $\Psi$ with respect + * to all independent variables, i.e. + * @f[ + * \frac{\partial\mathbf{r}(\mathbf{X})}{\partial\mathbf{X}} + * = + * \frac{\partial^{2}\Psi(\mathbf{X})}{\partial\mathbf{X} + * \otimes \partial\mathbf{X}} \Big\vert_{\mathbf{X}} + * @f] + * + * The values at the evaluation point $\mathbf{X}$ are 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_independent_variables$\times$n_independent_variables. + */ + virtual void + compute_linearization( + FullMatrix &linearization) const override; + + //@} + + }; // class ADHelperEnergyFunctional + + } // namespace AD } // namespace Differentiation diff --git a/source/differentiation/ad/ad_helpers.cc b/source/differentiation/ad/ad_helpers.cc index fed7fe4d93..60e4488adb 100644 --- a/source/differentiation/ad/ad_helpers.cc +++ b/source/differentiation/ad/ad_helpers.cc @@ -830,6 +830,208 @@ namespace Differentiation } + + /* ------------------ ADHelperEnergyFunctional ------------------ */ + + + + template + ADHelperEnergyFunctional:: + ADHelperEnergyFunctional(const unsigned int n_independent_variables) + : ADHelperCellLevelBase( + n_independent_variables, + 1) + {} + + + + template + void + ADHelperEnergyFunctional:: + register_energy_functional(const ad_type &energy) + { + Assert(this->n_dependent_variables() == 1, ExcInternalError()); + ADHelperBase::register_dependent_variable( + 0, energy); + } + + + + template + typename ADHelperEnergyFunctional::scalar_type + ADHelperEnergyFunctional::compute_energy() + const + { + if (this->keep_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.")); + + Assert( + this->n_dependent_variables() == 1, + ExcMessage( + "The ADHelperEnergyFunctional class expects there to be only one dependent variable.")); + + if (ADNumberTraits::is_taped == true) + { + Assert(this->active_tape() != numbers::invalid_tape_index, + ExcMessage("Invalid tape index")); + Assert(this->is_recording() == false, + ExcMessage( + "Cannot compute value 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())); + + return TapedDrivers::value( + this->active_tape(), this->independent_variable_values); + } + else + { + Assert(ADNumberTraits::is_tapeless == true, + ExcInternalError()); + Assert(this->independent_variables.size() == + this->n_independent_variables(), + ExcInternalError()); + + return TapelessDrivers::value( + this->dependent_variables); + } + } + + + + template + void + ADHelperEnergyFunctional::compute_residual( + Vector &gradient) const + { + if (this->keep_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.")); + + Assert( + this->n_dependent_variables() == 1, + ExcMessage( + "The ADHelperEnergyFunctional class expects there to be only one dependent variable.")); + + // We can neglect correctly initializing the entries as + // we'll be overwriting them immediately. + if (gradient.size() != this->n_independent_variables()) + gradient.reinit(this->n_independent_variables(), + true /*omit_zeroing_entries*/); + + if (ADNumberTraits::is_taped == true) + { + Assert(this->active_tape() != numbers::invalid_tape_index, + ExcMessage("Invalid tape index")); + Assert(this->is_recording() == false, + ExcMessage( + "Cannot compute gradient 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::gradient( + this->active_tape(), this->independent_variable_values, gradient); + } + else + { + Assert(ADNumberTraits::is_tapeless == true, + ExcInternalError()); + Assert(this->independent_variables.size() == + this->n_independent_variables(), + ExcInternalError()); + + TapelessDrivers::gradient( + this->independent_variables, this->dependent_variables, gradient); + } + } + + + + template + void + ADHelperEnergyFunctional:: + compute_linearization(FullMatrix &hessian) const + { + Assert(AD::ADNumberTraits::n_supported_derivative_levels >= 2, + ExcMessage( + "Cannot computed function Hessian: AD number type does" + "not support the calculation of second order derivatives.")); + + if (this->keep_values == false) + { + 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.")); + + Assert( + this->n_dependent_variables() == 1, + ExcMessage( + "The ADHelperEnergyFunctional class expects there to be only one dependent variable.")); + + // We can neglect correctly initializing the entries as + // we'll be overwriting them immediately. + if (hessian.m() != this->n_independent_variables() && + hessian.n() != this->n_independent_variables()) + hessian.reinit({this->n_independent_variables(), + this->n_independent_variables()}, + true /*omit_default_initialization*/); + + if (ADNumberTraits::is_taped == true) + { + Assert(this->active_tape() != 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::hessian( + this->active_tape(), this->independent_variable_values, hessian); + } + else + { + Assert(ADNumberTraits::is_tapeless == true, + ExcInternalError()); + Assert(this->independent_variables.size() == + this->n_independent_variables(), + ExcInternalError()); + TapelessDrivers::hessian( + this->independent_variables, this->dependent_variables, hessian); + } + } + + } // namespace AD } // namespace Differentiation diff --git a/source/differentiation/ad/ad_helpers.inst1.in b/source/differentiation/ad/ad_helpers.inst1.in index 9ede62db6e..418cdf3b78 100644 --- a/source/differentiation/ad/ad_helpers.inst1.in +++ b/source/differentiation/ad/ad_helpers.inst1.in @@ -37,6 +37,14 @@ for (number : REAL_SCALARS) template class ADHelperCellLevelBase; + // -------------------------- ADHelperEnergyFunctional ---------------------- + + template + class ADHelperEnergyFunctional; + + template + class ADHelperEnergyFunctional; + \} \} } @@ -65,6 +73,14 @@ for () template class ADHelperCellLevelBase::ad_type>; + // -------------------------- ADHelperEnergyFunctional ---------------------- + + template + class ADHelperEnergyFunctional::ad_type>; + + template + class ADHelperEnergyFunctional::ad_type>; + \} \} diff --git a/source/differentiation/ad/ad_helpers.inst2.in b/source/differentiation/ad/ad_helpers.inst2.in index 0b436715cd..be1bc98baa 100644 --- a/source/differentiation/ad/ad_helpers.inst2.in +++ b/source/differentiation/ad/ad_helpers.inst2.in @@ -49,6 +49,20 @@ for (number : REAL_SCALARS) template class ADHelperCellLevelBase; + // -------------------------- ADHelperEnergyFunctional ---------------------- + + template + class ADHelperEnergyFunctional; + + template + class ADHelperEnergyFunctional; + + template + class ADHelperEnergyFunctional; + + template + class ADHelperEnergyFunctional; + \} \} } @@ -89,6 +103,20 @@ for () template class ADHelperCellLevelBase::ad_type>; + // -------------------------- ADHelperEnergyFunctional ---------------------- + + template + class ADHelperEnergyFunctional::ad_type>; + + template + class ADHelperEnergyFunctional::ad_type>; + + template + class ADHelperEnergyFunctional::ad_type>; + + template + class ADHelperEnergyFunctional::ad_type>; + \} \} -- 2.39.5