}; // 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<double> solution (...); // Or another vector type
+ * std::vector<types::global_dof_index> local_dof_indices (...);
+ * const FEValuesExtractors::Vector u_fe (...);
+ * const FEValuesExtractors::Scalar msp_fe (...);
+ * const unsigned int u_block (...);
+ * const unsigned int msp_block (...);
+ * FEValues<dim> fe_values (...);
+ * const unsigned int n_q_points (...);
+ * FullMatrix<double> cell_matrix (...);
+ * Vector<double> 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<ADNumberType> 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<Tensor<2, dim, ADNumberType>> Grad_u(
+ * n_q_points, Tensor<2, dim, ADNumberType>());
+ * std::vector<Tensor<1, dim, ADNumberType>> 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<ADNumberType> 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<dim>() + 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 <code>dofs_per_cell</code> 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 <enum AD::NumberTypes ADNumberTypeCode,
+ typename ScalarType = double>
+ class ADHelperResidualLinearisation
+ : public ADHelperCellLevelBase<ADNumberTypeCode, ScalarType>
+ {
+ public:
+ /**
+ * Type definition for the floating point number type that is used in,
+ * and results from, all computations.
+ */
+ using scalar_type =
+ typename ADHelperBase<ADNumberTypeCode, ScalarType>::scalar_type;
+
+ /**
+ * Type definition for the auto-differentiation number type that is used
+ * in all computations.
+ */
+ using ad_type =
+ typename ADHelperBase<ADNumberTypeCode, ScalarType>::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<ad_type> &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<scalar_type> &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
+ * <code>n_dependent_variables</code>$\times$<code>n_independent_variables</code>.
+ */
+ virtual void
+ compute_linearization(
+ FullMatrix<scalar_type> &linearization) const override;
+
+ //@}
+
+ }; // class ADHelperResidualLinearisation
+
+
} // namespace AD
} // namespace Differentiation
}
+ /* ------------------- ADHelperResidualLinearisation ------------------- */
+
+
+
+ template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
+ ADHelperResidualLinearisation<ADNumberTypeCode, ScalarType>::
+ ADHelperResidualLinearisation(const unsigned int n_independent_variables,
+ const unsigned int n_dependent_variables)
+ : ADHelperCellLevelBase<ADNumberTypeCode, ScalarType>(
+ n_independent_variables,
+ n_dependent_variables)
+ {}
+
+
+
+ template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
+ void
+ ADHelperResidualLinearisation<ADNumberTypeCode, ScalarType>::
+ register_residual_vector(const std::vector<ad_type> &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<ADNumberTypeCode, ScalarType>::register_dependent_variable(
+ i, residual[i]);
+ }
+
+
+
+ template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
+ void
+ ADHelperResidualLinearisation<ADNumberTypeCode, ScalarType>::
+ compute_residual(Vector<scalar_type> &values) const
+ {
+ if ((ADNumberTraits<ad_type>::is_taped == true &&
+ this->taped_driver.keep_independent_values() == false) ||
+ ADNumberTraits<ad_type>::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<ad_type>::is_taped == true)
+ {
+ Assert(this->active_tape_index() !=
+ Numbers<ad_type>::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<ad_type, scalar_type>::values(
+ this->active_tape_index(),
+ this->n_dependent_variables(),
+ this->independent_variable_values,
+ values);
+ }
+ else
+ {
+ Assert(ADNumberTraits<ad_type>::is_tapeless == true,
+ ExcInternalError());
+ TapelessDrivers<ad_type, scalar_type>::values(
+ this->dependent_variables, values);
+ }
+ }
+
+
+
+ template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
+ void
+ ADHelperResidualLinearisation<ADNumberTypeCode, ScalarType>::
+ compute_linearization(FullMatrix<scalar_type> &jacobian) const
+ {
+ if ((ADNumberTraits<ad_type>::is_taped == true &&
+ this->taped_driver.keep_independent_values() == false) ||
+ ADNumberTraits<ad_type>::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<ad_type>::is_taped == true)
+ {
+ Assert(this->active_tape_index() !=
+ Numbers<ad_type>::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<ad_type, scalar_type>::jacobian(
+ this->active_tape_index(),
+ this->n_dependent_variables(),
+ this->independent_variable_values,
+ jacobian);
+ }
+ else
+ {
+ Assert(ADNumberTraits<ad_type>::is_tapeless == true,
+ ExcInternalError());
+ Assert(this->independent_variables.size() ==
+ this->n_independent_variables(),
+ ExcInternalError());
+ TapelessDrivers<ad_type, scalar_type>::jacobian(
+ this->independent_variables, this->dependent_variables, jacobian);
+ }
+ }
+
+
} // namespace AD
} // namespace Differentiation
projects / dealii.git / commitdiff