#include <deal.II/differentiation/ad/ad_helpers.h>

Inheritance diagram for Differentiation::AD::ResidualLinearization< ADNumberTypeCode, ScalarType >:

Public Types
using	scalar_type = typename HelperBase< ADNumberTypeCode, ScalarType >::scalar_type

using	ad_type = typename HelperBase< ADNumberTypeCode, ScalarType >::ad_type

Public Types inherited from Differentiation::AD::CellLevelBase< ADNumberTypeCode, ScalarType >
using	scalar_type = typename HelperBase< ADNumberTypeCode, ScalarType >::scalar_type

using	ad_type = typename HelperBase< ADNumberTypeCode, ScalarType >::ad_type

Public Types inherited from Differentiation::AD::HelperBase< ADNumberTypeCode, ScalarType >
using	scalar_type = typename AD::NumberTraits< ScalarType, ADNumberTypeCode >::scalar_type

using	ad_type = typename AD::NumberTraits< ScalarType, ADNumberTypeCode >::ad_type

Public Member Functions
Constructor / destructor
	ResidualLinearization (const unsigned int n_independent_variables, const unsigned int n_dependent_variables)

virtual	~ResidualLinearization ()=default

Dependent variables
void	register_residual_vector (const std::vector< ad_type > &residual)

virtual void	compute_residual (Vector< scalar_type > &residual) const override

virtual void	compute_linearization (FullMatrix< scalar_type > &linearization) const override

Public Member Functions inherited from Differentiation::AD::CellLevelBase< ADNumberTypeCode, ScalarType >
	CellLevelBase (const unsigned int n_independent_variables, const unsigned int n_dependent_variables)

virtual	~CellLevelBase ()=default

void	register_dof_values (const std::vector< scalar_type > &dof_values)

template<typename VectorType >
void	register_dof_values (const VectorType &values, const std::vector<::types::global_dof_index > &local_dof_indices)

const std::vector< ad_type > &	get_sensitive_dof_values () const

void	set_dof_values (const std::vector< scalar_type > &dof_values)

template<typename VectorType >
void	set_dof_values (const VectorType &values, const std::vector<::types::global_dof_index > &local_dof_indices)

Public Member Functions inherited from Differentiation::AD::HelperBase< ADNumberTypeCode, ScalarType >
	HelperBase (const unsigned int n_independent_variables, const unsigned int n_dependent_variables)

virtual	~HelperBase ()=default

std::size_t	n_independent_variables () const

std::size_t	n_dependent_variables () const

void	print (std::ostream &stream) const

void	print_values (std::ostream &stream) const

void	print_tape_stats (const typename Types< ad_type >::tape_index tape_index, std::ostream &stream) const

virtual void	reset (const unsigned int n_independent_variables=::numbers::invalid_unsigned_int, const unsigned int n_dependent_variables=::numbers::invalid_unsigned_int, const bool clear_registered_tapes=true)

bool	is_recording () const

Types< ad_type >::tape_index	active_tape_index () const

bool	is_registered_tape (const typename Types< ad_type >::tape_index tape_index) const

void	set_tape_buffer_sizes (const typename Types< ad_type >::tape_buffer_sizes obufsize=64 1024 1024, const typename Types< ad_type >::tape_buffer_sizes lbufsize=64 1024 1024, const typename Types< ad_type >::tape_buffer_sizes vbufsize=64 1024 1024, const typename Types< ad_type >::tape_buffer_sizes tbufsize=64 1024 1024)

bool	start_recording_operations (const typename Types< ad_type >::tape_index tape_index, const bool overwrite_tape=false, const bool keep_independent_values=true)

void	stop_recording_operations (const bool write_tapes_to_file=false)

void	activate_recorded_tape (const typename Types< ad_type >::tape_index tape_index)

bool	recorded_tape_requires_retaping (const typename Types< ad_type >::tape_index tape_index) const

bool	active_tape_requires_retaping () const

void	clear_active_tape ()

Additional Inherited Members
Static Public Member Functions inherited from Differentiation::AD::HelperBase< ADNumberTypeCode, ScalarType >
static void	configure_tapeless_mode (const unsigned int n_independent_variables, const bool ensure_persistent_setting=true)

Protected Member Functions inherited from Differentiation::AD::HelperBase< ADNumberTypeCode, ScalarType >
void	activate_tape (const typename Types< ad_type >::tape_index tape_index, const bool read_mode)

void	reset_registered_independent_variables ()

void	set_sensitivity_value (const unsigned int index, const scalar_type &value)

void	mark_independent_variable (const unsigned int index, ad_type &out) const

void	finalize_sensitive_independent_variables () const

void	initialize_non_sensitive_independent_variable (const unsigned int index, ad_type &out) const

unsigned int	n_registered_independent_variables () const

void	reset_registered_dependent_variables (const bool flag=false)

unsigned int	n_registered_dependent_variables () const

void	register_dependent_variable (const unsigned int index, const ad_type &func)

Protected Attributes inherited from Differentiation::AD::HelperBase< ADNumberTypeCode, ScalarType >
TapedDrivers< ad_type, scalar_type >	taped_driver

TapelessDrivers< ad_type, scalar_type >	tapeless_driver

std::vector< scalar_type >	independent_variable_values

std::vector< ad_type >	independent_variables

std::vector< bool >	registered_independent_variable_values

std::vector< bool >	registered_marked_independent_variables

std::vector< ad_type >	dependent_variables

std::vector< bool >	registered_marked_dependent_variables

Detailed Description

template<enum AD::NumberTypes ADNumberTypeCode, typename ScalarType = double>
class Differentiation::AD::ResidualLinearization< ADNumberTypeCode, ScalarType >

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):

// 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
  // EnergyFunctional for details on how to extend
  // support to taped AD numbers.
  using ADHelper = AD::ResidualLinearization<...>;
  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);
}

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 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

Definition at line 1539 of file ad_helpers.h.

Member Typedef Documentation

◆ scalar_type

template<enum AD::NumberTypes ADNumberTypeCode, typename ScalarType = double>

using Differentiation::AD::ResidualLinearization< ADNumberTypeCode, ScalarType >::scalar_type = typename HelperBase<ADNumberTypeCode, ScalarType>::scalar_type

Type definition for the floating point number type that is used in, and results from, all computations.

Definition at line 1548 of file ad_helpers.h.

◆ ad_type

template<enum AD::NumberTypes ADNumberTypeCode, typename ScalarType = double>

using Differentiation::AD::ResidualLinearization< ADNumberTypeCode, ScalarType >::ad_type = typename HelperBase<ADNumberTypeCode, ScalarType>::ad_type

Type definition for the auto-differentiation number type that is used in all computations.

Definition at line 1555 of file ad_helpers.h.

Constructor & Destructor Documentation

◆ ResidualLinearization()

template<enum AD::NumberTypes ADNumberTypeCode, typename ScalarType >

Differentiation::AD::ResidualLinearization< ADNumberTypeCode, ScalarType >::ResidualLinearization	(	const unsigned int	n_independent_variables,
		const unsigned int	n_dependent_variables
	)

The constructor for the class.

Parameters

[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.
[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.

Definition at line 1015 of file ad_helpers.cc.

◆ ~ResidualLinearization()

template<enum AD::NumberTypes ADNumberTypeCode, typename ScalarType = double>

virtual Differentiation::AD::ResidualLinearization< ADNumberTypeCode, ScalarType >::~ResidualLinearization ( )

virtualdefault

Destructor

Member Function Documentation

◆ register_residual_vector()

template<enum AD::NumberTypes ADNumberTypeCode, typename ScalarType >

void Differentiation::AD::ResidualLinearization< ADNumberTypeCode, ScalarType >::register_residual_vector ( const std::vector< ad_type > & residual )

Register the definition of the cell residual vector \(\mathbf{r}(\mathbf{X})\).

Parameters

[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.; For taped AD numbers, this operation is only valid in recording mode.

Definition at line 1027 of file ad_helpers.cc.

◆ compute_residual()

template<enum AD::NumberTypes ADNumberTypeCode, typename ScalarType >

void Differentiation::AD::ResidualLinearization< ADNumberTypeCode, ScalarType >::compute_residual ( Vector< scalar_type > & residual ) const

overridevirtual

Evaluation of the residual for a chosen set of degree of freedom values. This corresponds to the computation of the residual vector, i.e.

\[ \mathbf{r}(\mathbf{X}) \vert_{\mathbf{X}} \]

The values at the evaluation point \(\mathbf{X}\) are obtained by calling CellLevelBase::set_dof_values().

Parameters

[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 residual vector has a length corresponding to n_dependent_variables.

Implements Differentiation::AD::CellLevelBase< ADNumberTypeCode, ScalarType >.

Definition at line 1041 of file ad_helpers.cc.

◆ compute_linearization()

template<enum AD::NumberTypes ADNumberTypeCode, typename ScalarType >

void Differentiation::AD::ResidualLinearization< ADNumberTypeCode, ScalarType >::compute_linearization ( FullMatrix< scalar_type > & linearization ) const

overridevirtual

Compute 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.

\[ \frac{\partial\mathbf{r}(\mathbf{X})}{\partial\mathbf{X}} \]

The values at the evaluation point \(\mathbf{X}\) are obtained by calling CellLevelBase::set_dof_values().

Parameters

[out] linearization A FullMatrix representing the linearization of the residual vector. The output linearization matrix has dimensions corresponding to n_dependent_variables \(\times\)n_independent_variables.

Implements Differentiation::AD::CellLevelBase< ADNumberTypeCode, ScalarType >.

Definition at line 1095 of file ad_helpers.cc.

The documentation for this class was generated from the following files:

deal.II/differentiation/ad/ad_helpers.h
source/differentiation/ad/ad_helpers.cc

Public Types

Public Member Functions

Additional Inherited Members

Detailed Description

template<enum AD::NumberTypes ADNumberTypeCode, typename ScalarType = double> class Differentiation::AD::ResidualLinearization< ADNumberTypeCode, ScalarType >

Member Typedef Documentation

◆ scalar_type

◆ ad_type

Constructor & Destructor Documentation

◆ ResidualLinearization()

◆ ~ResidualLinearization()

Member Function Documentation

◆ register_residual_vector()

◆ compute_residual()

◆ compute_linearization()

template<enum AD::NumberTypes ADNumberTypeCode, typename ScalarType = double>
class Differentiation::AD::ResidualLinearization< ADNumberTypeCode, ScalarType >