15#ifndef dealii_differentiation_ad_ad_helpers_h
16#define dealii_differentiation_ad_ad_helpers_h
20#if defined(DEAL_II_WITH_ADOLC) || defined(DEAL_II_TRILINOS_WITH_SACADO)
242 print(std::ostream &stream)
const;
265 std::ostream &stream)
const;
452 const bool keep_independent_values =
true);
923 template <
typename VectorType>
926 const VectorType &values,
927 const std::vector<::types::global_dof_index> &local_dof_indices);
949 const std::vector<ad_type> &
998 template <
typename VectorType>
1001 const VectorType &values,
1002 const std::vector<::types::global_dof_index> &local_dof_indices);
1647 template <
int dim,
typename ExtractorType>
1662 static const unsigned int n_components = 1;
1667 static const unsigned int rank = 0;
1672 template <
typename NumberType>
1677 "The number of components doesn't match that of the corresponding tensor type.");
1680 "The rank doesn't match that of the corresponding tensor type.");
1689 template <
typename NumberType>
1695 template <
typename NumberType>
1701 static inline unsigned int
1728 template <
typename IndexType =
unsigned int,
int rank_in>
1752 static const unsigned int n_components = dim;
1757 static const unsigned int rank = 1;
1762 template <
typename NumberType>
1767 "The number of components doesn't match that of the corresponding tensor type.");
1770 "The rank doesn't match that of the corresponding tensor type.");
1775 template <
typename NumberType>
1781 template <
typename NumberType>
1787 static inline unsigned int
1812 template <
int rank_in>
1836 template <
typename IndexType =
unsigned int,
int rank_in>
1843 "Cannot extract more table indices than the input table has!");
1862 static const unsigned int n_components =
1868 static const unsigned int rank = 1;
1873 template <
typename NumberType>
1879 template <
typename NumberType>
1885 template <
typename NumberType>
1891 static inline unsigned int
1915 template <
int rank_in>
1939 template <
typename IndexType =
unsigned int,
int rank_in>
1946 "Cannot extract more table indices than the input table has!");
1965 static const unsigned int n_components =
1976 template <
typename NumberType>
1982 template <
typename NumberType>
1988 template <
typename NumberType>
1994 static inline unsigned int
2018 template <
int rank_in>
2044 template <
typename IndexType =
unsigned int,
int rank_in>
2051 "Cannot extract more table indices than the input table has!");
2070 static const unsigned int n_components =
2081 template <
typename NumberType>
2087 template <
typename NumberType>
2093 template <
typename NumberType>
2099 static inline unsigned int
2125 template <
int rank_in>
2151 template <
typename IndexType =
unsigned int,
int rank_in>
2158 "Cannot extract more table indices than the input table has!");
2171 template <
int dim,
typename NumberType,
typename ExtractorType>
2193 template <
typename ExtractorType_Row,
typename ExtractorType_Col>
2207 template <
int rank,
int dim,
typename NumberType>
2231 template <
int ,
int dim,
typename NumberType>
2256 template <
int ,
int dim,
typename NumberType>
2280 template <
int ,
int dim,
typename NumberType>
2322 template <
int dim,
typename NumberType,
typename ExtractorType_Field>
2354 std::vector<IndexType>
2363 std::vector<IndexType> indices(n_components);
2364 std::iota(indices.begin(), indices.end(),
comp_first);
2377 template <
int dim,
typename IndexType =
unsigned int>
2378 std::vector<IndexType>
2392 std::vector<IndexType> indices(n_components);
2393 std::iota(indices.begin(), indices.end(),
comp_first);
2401 std::vector<IndexType> indices =
2406 for (
unsigned int i = 0; i < indices.size(); ++i)
2433 template <
typename TensorType,
typename NumberType>
2449 template <
int dim,
typename NumberType>
2466 template <
typename NumberType>
2483 template <
int dim,
typename NumberType>
2512 typename NumberType,
2513 template <
int,
int,
typename>
2531 typename NumberType,
2532 template <
int,
int,
typename>
2549 template <
typename NumberType>
2550 inline const NumberType &
2565 typename NumberType,
2566 template <
int,
int,
typename>
2584 typename NumberType,
2585 template <
int,
int,
typename>
2589 const unsigned int index)
2602 template <
typename NumberType>
2640 :
public HelperBase<ADNumberTypeCode, ScalarType>
2785 template <
typename ValueType,
typename ExtractorType>
2810 const std::vector<ad_type> &
2838 template <
typename ExtractorType>
2896 template <
typename ValueType,
typename ExtractorType>
2921 const bool symmetric_component,
3252 template <
typename ExtractorType_Row>
3253 static typename internal::
3254 ScalarFieldGradient<dim, scalar_type, ExtractorType_Row>::type
3296 template <
typename ExtractorType_Row,
typename ExtractorType_Col>
3592 template <
typename ValueType,
typename ExtractorType>
3638 template <
typename ExtractorType_Row>
3639 static typename internal::
3640 VectorFieldValue<dim, scalar_type, ExtractorType_Row>::type
3691 template <
typename ExtractorType_Row,
typename ExtractorType_Col>
3762 template <enum AD::NumberTypes ADNumberTypeCode,
typename ScalarType>
3763 template <
typename VectorType>
3766 const VectorType &values,
3767 const std::vector<::types::global_dof_index> &local_dof_indices)
3774 local_dof_indices.size() ==
this->n_independent_variables(),
3776 "Degree of freedom index vector size does not match number of independent variables"));
3777 for (
unsigned int i = 0; i < this->n_independent_variables(); ++i)
3779 Assert(this->registered_independent_variable_values[i] ==
false,
3780 ExcMessage(
"Independent variables already registered."));
3787 template <enum AD::NumberTypes ADNumberTypeCode,
typename ScalarType>
3788 template <
typename VectorType>
3791 const VectorType &values,
3792 const std::vector<::types::global_dof_index> &local_dof_indices)
3794 Assert(local_dof_indices.size() ==
this->n_independent_variables(),
3796 "Vector size does not match number of independent variables"));
3797 for (
unsigned int i = 0; i < this->n_independent_variables(); ++i)
3799 i, values[local_dof_indices[i]]);
3811 template <
typename ValueType,
typename ExtractorType>
3822 const std::vector<unsigned int> index_set(
3823 internal::extract_field_component_indices<dim>(extractor));
3824 for (
const unsigned int index : index_set)
3827 this->registered_independent_variable_values[index] ==
false,
3829 "Overlapping indices for independent variables. "
3830 "One or more indices associated with the field that "
3831 "is being registered as an independent variable have "
3832 "already been associated with another field. This suggests "
3833 "that the component offsets used to construct their counterpart "
3834 "extractors are incompatible with one another. Make sure that "
3835 "the first component for each extractor properly takes into "
3836 "account the dimensionality of the preceding fields."));
3839 set_independent_variable(value, extractor);
3847 template <
typename ValueType,
typename ExtractorType>
3853 const std::vector<unsigned int> index_set(
3854 internal::extract_field_component_indices<dim>(extractor));
3855 for (
unsigned int i = 0; i < index_set.size(); ++i)
3857 set_sensitivity_value(
3859 internal::Extractor<dim, ExtractorType>::symmetric_component(i),
3869 template <
typename ExtractorType>
3870 typename internal::Extractor<dim, ExtractorType>::template tensor_type<
3875 if (ADNumberTraits<ad_type>::is_taped ==
true)
3877 Assert(this->active_tape_index() !=
3884 this->finalize_sensitive_independent_variables();
3885 Assert(this->independent_variables.size() ==
3886 this->n_independent_variables(),
3888 this->n_independent_variables()));
3890 const std::vector<unsigned int> index_set(
3891 internal::extract_field_component_indices<dim>(extractor));
3892 typename internal::Extractor<dim,
3896 for (
unsigned int i = 0; i < index_set.size(); ++i)
3898 const unsigned int index = index_set[i];
3900 Assert(this->registered_independent_variable_values[index] ==
true,
3903 this->independent_variables[
index];
3918 template <
typename ExtractorType_Row>
3919 typename internal::ScalarFieldGradient<
3930 ScalarFieldGradient<dim, scalar_type, ExtractorType_Row>::type
out;
3934 const std::vector<unsigned int> row_index_set(
3935 internal::extract_field_component_indices<dim>(
extractor_row));
3937 ExcMessage(
"Not all tensor components have been extracted!"));
3938 for (
unsigned int r = 0; r < row_index_set.size(); ++r)
3949 template <
typename ExtractorType_Row,
typename ExtractorType_Col>
3950 typename internal::ScalarFieldHessian<
3966 using HessianType =
typename InternalHessian::type;
3982 const std::vector<unsigned int> row_index_set(
3983 internal::extract_field_component_indices<dim>(
3986 internal::extract_field_component_indices<dim>(
3989 for (
unsigned int index = 0;
3990 index < HessianType::n_independent_components;
3994 HessianType::unrolled_to_component_indices(index);
3995 const unsigned int r =
3996 InternalExtractorRow::local_component(
ti_out, 0);
3997 const unsigned int c =
3998 InternalExtractorCol::local_component(
ti_out,
4017 template <
typename ValueType,
typename ExtractorType>
4023 const std::vector<unsigned int> index_set(
4024 internal::extract_field_component_indices<dim>(extractor));
4025 for (
unsigned int i = 0; i < index_set.size(); ++i)
4027 Assert(this->registered_marked_dependent_variables[index_set[i]] ==
4029 ExcMessage(
"Overlapping indices for dependent variables."));
4040 template <
typename ExtractorType_Row>
4051 typename internal::VectorFieldValue<dim, scalar_type, ExtractorType_Row>
::
4056 const std::vector<unsigned int> row_index_set(
4057 internal::extract_field_component_indices<dim>(
extractor_row));
4059 ExcMessage(
"Not all tensor components have been extracted!"));
4060 for (
unsigned int r = 0; r < row_index_set.size(); ++r)
4071 template <
typename ExtractorType_Row,
typename ExtractorType_Col>
4104 const std::vector<unsigned int> row_index_set(
4105 internal::extract_field_component_indices<dim>(
4108 internal::extract_field_component_indices<dim>(
4111 for (
unsigned int index = 0;
4117 const unsigned int r =
4118 InternalExtractorRow::local_component(
ti_out, 0);
4119 const unsigned int c =
4120 InternalExtractorCol::local_component(
ti_out,
void register_dof_values(const VectorType &values, const std::vector<::types::global_dof_index > &local_dof_indices)
virtual ~CellLevelBase()=default
const std::vector< ad_type > & get_sensitive_dof_values() const
typename HelperBase< ADNumberTypeCode, ScalarType >::ad_type ad_type
void set_dof_values(const std::vector< scalar_type > &dof_values)
virtual void compute_residual(Vector< scalar_type > &residual) const =0
void register_dof_values(const std::vector< scalar_type > &dof_values)
virtual void compute_linearization(FullMatrix< scalar_type > &linearization) const =0
void set_dof_values(const VectorType &values, const std::vector<::types::global_dof_index > &local_dof_indices)
typename HelperBase< ADNumberTypeCode, ScalarType >::scalar_type scalar_type
void register_energy_functional(const ad_type &energy)
typename HelperBase< ADNumberTypeCode, ScalarType >::ad_type ad_type
void compute_residual(Vector< scalar_type > &residual) const override
virtual void compute_linearization(FullMatrix< scalar_type > &linearization) const override
virtual ~EnergyFunctional()=default
scalar_type compute_energy() const
typename HelperBase< ADNumberTypeCode, ScalarType >::scalar_type scalar_type
unsigned int n_registered_dependent_variables() const
TapedDrivers< ad_type, scalar_type > taped_driver
void stop_recording_operations(const bool write_tapes_to_file=false)
std::vector< bool > registered_marked_independent_variables
typename AD::NumberTraits< ScalarType, ADNumberTypeCode >::scalar_type scalar_type
std::vector< bool > registered_independent_variable_values
Types< ad_type >::tape_index active_tape_index() const
void print_tape_stats(const typename Types< ad_type >::tape_index tape_index, std::ostream &stream) const
std::size_t n_independent_variables() const
bool start_recording_operations(const typename Types< ad_type >::tape_index tape_index, const bool overwrite_tape=false, const bool keep_independent_values=true)
std::vector< scalar_type > independent_variable_values
void reset_registered_independent_variables()
void mark_independent_variable(const unsigned int index, ad_type &out) const
bool is_registered_tape(const typename Types< ad_type >::tape_index tape_index) const
std::vector< bool > registered_marked_dependent_variables
void reset_registered_dependent_variables(const bool flag=false)
void activate_tape(const typename Types< ad_type >::tape_index tape_index, const bool read_mode)
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)
unsigned int n_registered_independent_variables() const
void print(std::ostream &stream) const
void set_sensitivity_value(const unsigned int index, const scalar_type &value)
bool active_tape_requires_retaping() const
std::vector< ad_type > independent_variables
virtual ~HelperBase()=default
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)
TapelessDrivers< ad_type, scalar_type > tapeless_driver
bool recorded_tape_requires_retaping(const typename Types< ad_type >::tape_index tape_index) const
std::vector< ad_type > dependent_variables
bool is_recording() const
std::size_t n_dependent_variables() const
void print_values(std::ostream &stream) const
void register_dependent_variable(const unsigned int index, const ad_type &func)
void finalize_sensitive_independent_variables() const
static void configure_tapeless_mode(const unsigned int n_independent_variables, const bool ensure_persistent_setting=true)
void initialize_non_sensitive_independent_variable(const unsigned int index, ad_type &out) const
typename AD::NumberTraits< ScalarType, ADNumberTypeCode >::ad_type ad_type
void activate_recorded_tape(const typename Types< ad_type >::tape_index tape_index)
void set_independent_variable(const ValueType &value, const ExtractorType &extractor)
bool is_symmetric_independent_variable(const unsigned int index) const
virtual ~PointLevelFunctionsBase()=default
typename HelperBase< ADNumberTypeCode, ScalarType >::ad_type ad_type
static constexpr unsigned int dimension
void set_independent_variables(const std::vector< scalar_type > &values)
typename HelperBase< ADNumberTypeCode, ScalarType >::scalar_type scalar_type
std::vector< bool > symmetric_independent_variables
unsigned int n_symmetric_independent_variables() const
void register_independent_variables(const std::vector< scalar_type > &values)
void set_sensitivity_value(const unsigned int index, const bool symmetric_component, const scalar_type &value)
void register_independent_variable(const ValueType &value, const ExtractorType &extractor)
const std::vector< ad_type > & get_sensitive_variables() 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) override
internal::Extractor< dim, ExtractorType >::template tensor_type< typename HelperBase< ADNumberTypeCode, ScalarType >::ad_type > get_sensitive_variables(const ExtractorType &extractor) const
virtual ~ResidualLinearization()=default
typename HelperBase< ADNumberTypeCode, ScalarType >::scalar_type scalar_type
virtual void compute_residual(Vector< scalar_type > &residual) const override
virtual void compute_linearization(FullMatrix< scalar_type > &linearization) const override
typename HelperBase< ADNumberTypeCode, ScalarType >::ad_type ad_type
void register_residual_vector(const std::vector< ad_type > &residual)
static internal::ScalarFieldHessian< dim, scalar_type, ExtractorType_Row, ExtractorType_Col >::type extract_hessian_component(const FullMatrix< scalar_type > &hessian, const ExtractorType_Row &extractor_row, const ExtractorType_Col &extractor_col)
typename HelperBase< ADNumberTypeCode, ScalarType >::ad_type ad_type
void register_dependent_variable(const ad_type &func)
static internal::ScalarFieldGradient< dim, scalar_type, ExtractorType_Row >::type extract_gradient_component(const Vector< scalar_type > &gradient, const ExtractorType_Row &extractor_row)
void compute_gradient(Vector< scalar_type > &gradient) const
scalar_type compute_value() const
void compute_hessian(FullMatrix< scalar_type > &hessian) const
typename HelperBase< ADNumberTypeCode, ScalarType >::scalar_type scalar_type
virtual ~ScalarFunction()=default
virtual ~VectorFunction()=default
typename HelperBase< ADNumberTypeCode, ScalarType >::ad_type ad_type
void compute_values(Vector< scalar_type > &values) const
void register_dependent_variables(const std::vector< ad_type > &funcs)
void compute_jacobian(FullMatrix< scalar_type > &jacobian) const
typename HelperBase< ADNumberTypeCode, ScalarType >::scalar_type scalar_type
static internal::VectorFieldJacobian< dim, scalar_type, ExtractorType_Row, ExtractorType_Col >::type extract_jacobian_component(const FullMatrix< scalar_type > &jacobian, const ExtractorType_Row &extractor_row, const ExtractorType_Col &extractor_col)
static internal::VectorFieldValue< dim, scalar_type, ExtractorType_Row >::type extract_value_component(const Vector< scalar_type > &values, const ExtractorType_Row &extractor_row)
void register_dependent_variable(const ValueType &funcs, const ExtractorType &extractor)
static DEAL_II_HOST constexpr unsigned int component_to_unrolled_index(const TableIndices< rank_ > &indices)
static constexpr unsigned int rank
static constexpr unsigned int component_to_unrolled_index(const TableIndices< rank_ > &indices)
static constexpr TableIndices< rank_ > unrolled_to_component_indices(const unsigned int i)
static constexpr unsigned int n_independent_components
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
#define Assert(cond, exc)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
ScalarFieldHessian< dim, NumberType, ExtractorType_Field, ExtractorType_Derivative > VectorFieldJacobian
NumberType get_tensor_entry(const TensorType< rank, dim, NumberType > &t, const unsigned int unrolled_index)
void set_tensor_entry(TensorType &t, const unsigned int unrolled_index, const NumberType &value)
std::vector< IndexType > extract_field_component_indices(const ExtractorType &extractor, const bool ignore_symmetries=true)
ScalarFieldGradient< dim, NumberType, ExtractorType_Field > VectorFieldValue
void set_dof_values(const DoFCellAccessor< dim, spacedim, lda > &cell, const Vector< number > &local_values, OutputVector &values, const bool perform_check)
static const unsigned int invalid_unsigned_int
const InputIterator OutputIterator out
static const Types< ADNumberType >::tape_index invalid_tape_index
typename Extractor< dim, ExtractorType >::template tensor_type< NumberType > type
typename HessianType< ExtractorType_Row, ExtractorType_Col >::template type< rank, dim, NumberType > type