17#if defined(DEAL_II_WITH_ADOLC) || defined(DEAL_II_TRILINOS_WITH_SACADO)
22# include <type_traits>
36 template <enum AD::NumberTypes ADNumberTypeCode,
typename ScalarType>
54 "Floating point/arithmetic numbers have no derivatives."));
58 "The AD number type does not support the calculation of any derivatives."));
78 template <enum AD::NumberTypes ADNumberTypeCode,
typename ScalarType>
90 template <enum AD::NumberTypes ADNumberTypeCode,
typename ScalarType>
102 template <enum AD::NumberTypes ADNumberTypeCode,
typename ScalarType>
105 const unsigned int index,
118 "Cannot change the value of an independent variable "
119 "of the tapeless variety while this class is not set "
120 "in recording operations."));
131 "Trying to set the value of a non-existent independent variable."));
139 template <enum AD::NumberTypes ADNumberTypeCode,
typename ScalarType>
142 const unsigned int index,
154 "Need to extract sensitivities in the order they're created."));
163 "The marking of independent variables is only valid "
164 "during recording."));
177 template <enum AD::NumberTypes ADNumberTypeCode,
typename ScalarType>
184 ExcMessage(
"Not all values of sensitivities have been recorded!"));
201 template <enum AD::NumberTypes ADNumberTypeCode,
typename ScalarType>
214 "The initialization of non-sensitive independent variables is "
215 "only valid outside of recording operations."));
226 template <enum AD::NumberTypes ADNumberTypeCode,
typename ScalarType>
238 template <enum AD::NumberTypes ADNumberTypeCode,
typename ScalarType>
247 template <enum AD::NumberTypes ADNumberTypeCode,
typename ScalarType>
259 template <enum AD::NumberTypes ADNumberTypeCode,
typename ScalarType>
268 template <enum AD::NumberTypes ADNumberTypeCode,
typename ScalarType>
280 template <enum AD::NumberTypes ADNumberTypeCode,
typename ScalarType>
293 template <enum AD::NumberTypes ADNumberTypeCode,
typename ScalarType>
306 template <enum AD::NumberTypes ADNumberTypeCode,
typename ScalarType>
311 const std::ios_base::fmtflags stream_flags(stream.flags());
313 stream.setf(std::ios_base::boolalpha);
317 stream << std::flush;
322 stream <<
"Registered independent variables: " <<
'\n';
328 stream <<
"Independent variable values: " <<
'\n';
331 stream <<
"Registered marked independent variables: " <<
'\n';
338 stream <<
"Dependent variable values: " <<
'\n';
344 stream <<
"Registered dependent variables: " <<
'\n';
351 stream.flags(stream_flags);
356 template <enum AD::NumberTypes ADNumberTypeCode,
typename ScalarType>
359 std::ostream &stream)
const
370 template <enum AD::NumberTypes ADNumberTypeCode,
typename ScalarType>
374 std::ostream &stream)
const
382 this->
taped_driver.print_tape_stats(tape_index, stream);
387 template <enum AD::NumberTypes ADNumberTypeCode,
typename ScalarType>
392 const bool clear_registered_tapes)
394 const unsigned int new_n_independent_variables =
397 this->n_independent_variables());
398 const unsigned int new_n_dependent_variables =
401 this->n_dependent_variables());
430 new_n_independent_variables,
433 std::vector<bool>(new_n_independent_variables,
false);
435 std::vector<bool>(new_n_independent_variables,
false);
437 std::vector<ad_type>(new_n_dependent_variables,
440 std::vector<bool>(new_n_dependent_variables,
false);
445 template <enum AD::NumberTypes ADNumberTypeCode,
typename ScalarType>
449 const bool ensure_persistent_setting)
458 if (ensure_persistent_setting ==
true)
471 template <enum AD::NumberTypes ADNumberTypeCode,
typename ScalarType>
481 template <enum AD::NumberTypes ADNumberTypeCode,
typename ScalarType>
494 template <enum AD::NumberTypes ADNumberTypeCode,
typename ScalarType>
507 template <enum AD::NumberTypes ADNumberTypeCode,
typename ScalarType>
519 template <enum AD::NumberTypes ADNumberTypeCode,
typename ScalarType>
523 const bool read_mode)
530 ExcMessage(
"Tape index exceeds maximum allowable value"));
537 if (read_mode ==
true)
544 ExcMessage(
"Not all dependent variables have been set!"));
556 template <enum AD::NumberTypes ADNumberTypeCode,
typename ScalarType>
575 template <enum AD::NumberTypes ADNumberTypeCode,
typename ScalarType>
579 const bool overwrite_tape,
580 const bool keep_independent_values)
583 const bool read_mode =
false;
587 if (overwrite_tape !=
true)
602 keep_independent_values);
613 "Tape recording is unexpectedly still enabled."));
622 Assert(ADNumberTraits<ad_type>::is_tapeless ==
true,
627 tapeless_driver.allow_dependent_variable_marking();
632 activate_tape(tape_index, read_mode);
635 return is_recording();
640 template <enum AD::NumberTypes ADNumberTypeCode,
typename ScalarType>
643 const bool write_tapes_to_file)
649 ExcMessage(
"Not all values of sensitivities have been recorded!"));
662 ExcMessage(
"Not all dependent variables have been set!"));
676 template <enum AD::NumberTypes ADNumberTypeCode,
typename ScalarType>
679 const unsigned int index,
685 "This dependent variable has already been registered."));
693 "Must be recording when registering dependent variables."));
699 registered_marked_dependent_variables[index] =
true;
708 template <enum AD::NumberTypes ADNumberTypeCode,
typename ScalarType>
718 template <enum AD::NumberTypes ADNumberTypeCode,
typename ScalarType>
721 const std::vector<scalar_type> &dof_values)
727 Assert(dof_values.size() == this->n_independent_variables(),
729 "Vector size does not match number of independent variables"));
730 for (
unsigned int i = 0; i < this->n_independent_variables(); ++i)
733 ExcMessage(
"Independent variable value already registered."));
735 set_dof_values(dof_values);
740 template <enum AD::NumberTypes ADNumberTypeCode,
typename ScalarType>
757 this->n_independent_variables(),
759 this->n_independent_variables()));
766 template <enum AD::NumberTypes ADNumberTypeCode,
typename ScalarType>
769 const std::vector<scalar_type> &values)
777 Assert(values.size() == this->n_independent_variables(),
779 "Vector size does not match number of independent variables"));
780 for (
unsigned int i = 0; i < this->n_independent_variables(); ++i)
791 template <enum AD::NumberTypes ADNumberTypeCode,
typename ScalarType>
799 template <enum AD::NumberTypes ADNumberTypeCode,
typename ScalarType>
811 template <enum AD::NumberTypes ADNumberTypeCode,
typename ScalarType>
816 this->
taped_driver.keep_independent_values() ==
false) ||
821 this->n_independent_variables(),
823 "Not all values of sensitivities have been registered or subsequently set!"));
827 ExcMessage(
"Not all dependent variables have been registered."));
832 "The EnergyFunctional class expects there to be only one dependent variable."));
841 "Cannot compute value while tape is being recorded."));
843 this->n_independent_variables(),
845 this->n_independent_variables()));
855 this->n_independent_variables(),
857 this->n_independent_variables()));
865 template <enum AD::NumberTypes ADNumberTypeCode,
typename ScalarType>
871 this->
taped_driver.keep_independent_values() ==
false) ||
876 this->n_independent_variables(),
878 "Not all values of sensitivities have been registered or subsequently set!"));
880 Assert(this->n_registered_dependent_variables() ==
881 this->n_dependent_variables(),
882 ExcMessage(
"Not all dependent variables have been registered."));
885 this->n_dependent_variables() == 1,
887 "The EnergyFunctional class expects there to be only one dependent variable."));
892 if (gradient.size() != this->n_independent_variables())
893 gradient.reinit(this->n_independent_variables(),
898 Assert(this->active_tape_index() !=
901 Assert(this->is_recording() ==
false,
903 "Cannot compute gradient while tape is being recorded."));
905 this->n_independent_variables(),
907 this->n_independent_variables()));
917 Assert(this->independent_variables.size() ==
918 this->n_independent_variables(),
920 this->n_independent_variables()));
922 this->tapeless_driver.gradient(this->independent_variables,
923 this->dependent_variables,
930 template <enum AD::NumberTypes ADNumberTypeCode,
typename ScalarType>
937 "Cannot computed function Hessian: AD number type does "
938 "not support the calculation of second order derivatives."));
945 this->n_independent_variables(),
947 "Not all values of sensitivities have been registered or subsequently set!"));
951 ExcMessage(
"Not all dependent variables have been registered."));
956 "The EnergyFunctional class expects there to be only one dependent variable."));
961 if (hessian.m() != this->n_independent_variables() ||
962 hessian.n() != this->n_independent_variables())
963 hessian.reinit({this->n_independent_variables(),
964 this->n_independent_variables()},
969 Assert(this->active_tape_index() !=
972 Assert(this->is_recording() ==
false,
974 "Cannot compute hessian while tape is being recorded."));
975 Assert(this->independent_variable_values.size() ==
976 this->n_independent_variables(),
978 this->n_independent_variables()));
988 Assert(this->independent_variables.size() ==
989 this->n_independent_variables(),
991 this->n_independent_variables()));
993 this->tapeless_driver.hessian(this->independent_variables,
994 this->dependent_variables,
1004 template <enum AD::NumberTypes ADNumberTypeCode,
typename ScalarType>
1014 template <enum AD::NumberTypes ADNumberTypeCode,
typename ScalarType>
1019 Assert(residual.size() == this->n_dependent_variables(),
1021 "Vector size does not match number of dependent variables"));
1029 template <enum AD::NumberTypes ADNumberTypeCode,
typename ScalarType>
1035 this->
taped_driver.keep_independent_values() ==
false) ||
1040 this->n_independent_variables(),
1042 "Not all values of sensitivities have been registered or subsequently set!"));
1046 ExcMessage(
"Not all dependent variables have been registered."));
1051 if (values.size() != this->n_dependent_variables())
1062 "Cannot compute values while tape is being recorded."));
1064 this->n_independent_variables(),
1066 this->n_independent_variables()));
1083 template <enum AD::NumberTypes ADNumberTypeCode,
typename ScalarType>
1089 this->
taped_driver.keep_independent_values() ==
false) ||
1094 this->n_independent_variables(),
1096 "Not all values of sensitivities have been registered or subsequently set!"));
1100 ExcMessage(
"Not all dependent variables have been registered."));
1105 if (jacobian.
m() != this->n_dependent_variables() ||
1106 jacobian.
n() != this->n_independent_variables())
1108 this->n_independent_variables()},
1118 "Cannot compute hessian while tape is being recorded."));
1120 this->n_independent_variables(),
1122 this->n_independent_variables()));
1134 this->n_independent_variables(),
1136 this->n_independent_variables()));
1152 typename ScalarType>
1165 typename ScalarType>
1170 const bool clear_registered_tapes)
1174 clear_registered_tapes);
1176 const unsigned int new_n_independent_variables =
1179 this->n_independent_variables());
1181 std::vector<bool>(new_n_independent_variables,
false);
1188 typename ScalarType>
1202 typename ScalarType>
1216 typename ScalarType>
1225 Assert(values.size() == this->n_independent_variables(),
1227 "Vector size does not match number of independent variables"));
1228 for (
unsigned int i = 0; i < this->n_independent_variables(); ++i)
1231 ExcMessage(
"Independent variable value already registered."));
1240 typename ScalarType>
1243 ScalarType>::ad_type> &
1259 this->n_independent_variables(),
1261 this->n_independent_variables()));
1270 typename ScalarType>
1274 const bool symmetric_component,
1282 "Trying to set the symmetry flag of a non-existent independent variable."));
1292 typename ScalarType>
1303 Assert(values.size() == this->n_independent_variables(),
1305 "Vector size does not match number of independent variables"));
1306 for (
unsigned int i = 0; i < this->n_independent_variables(); ++i)
1319 typename ScalarType>
1331 typename ScalarType>
1345 typename ScalarType>
1350 this->
taped_driver.keep_independent_values() ==
false) ||
1355 this->n_independent_variables(),
1357 "Not all values of sensitivities have been registered or subsequently set!"));
1361 ExcMessage(
"Not all dependent variables have been registered."));
1366 "The ScalarFunction class expects there to be only one dependent variable."));
1375 "Cannot compute values while tape is being recorded."));
1377 this->n_independent_variables(),
1379 this->n_independent_variables()));
1395 typename ScalarType>
1401 this->
taped_driver.keep_independent_values() ==
false) ||
1406 this->n_independent_variables(),
1408 "Not all values of sensitivities have been registered or subsequently set!"));
1412 ExcMessage(
"Not all dependent variables have been registered."));
1417 "The ScalarFunction class expects there to be only one dependent variable."));
1422 if (gradient.size() != this->n_independent_variables())
1423 gradient.reinit(this->n_independent_variables(),
1433 "Cannot compute gradient while tape is being recorded."));
1435 this->n_independent_variables(),
1437 this->n_independent_variables()));
1448 this->n_independent_variables(),
1450 this->n_independent_variables()));
1458 for (
unsigned int i = 0; i < this->n_independent_variables(); ++i)
1469 typename ScalarType>
1476 "Cannot computed function Hessian: AD number type does "
1477 "not support the calculation of second order derivatives."));
1480 this->
taped_driver.keep_independent_values() ==
false))
1484 this->n_independent_variables(),
1486 "Not all values of sensitivities have been registered or subsequently set!"));
1490 ExcMessage(
"Not all dependent variables have been registered."));
1495 "The ScalarFunction class expects there to be only one dependent variable."));
1500 if (hessian.m() != this->n_independent_variables() ||
1501 hessian.n() != this->n_independent_variables())
1502 hessian.reinit({this->n_independent_variables(),
1503 this->n_independent_variables()},
1513 "Cannot compute Hessian while tape is being recorded."));
1515 this->n_independent_variables(),
1517 this->n_independent_variables()));
1528 this->n_independent_variables(),
1530 this->n_independent_variables()));
1538 for (
unsigned int i = 0; i < this->n_independent_variables(); ++i)
1539 for (
unsigned int j = 0; j < i + 1; ++j)
1544 hessian[i][j] *= 0.25;
1546 hessian[j][i] *= 0.25;
1553 hessian[i][j] *= 0.5;
1555 hessian[j][i] *= 0.5;
1564 typename ScalarType>
1583 const std::vector<unsigned int> row_index_set(
1585 const std::vector<unsigned int> col_index_set(
1592 hessian[row_index_set[0]][col_index_set[0]]);
1601 typename ScalarType>
1621 const std::vector<unsigned int> row_index_set(
1623 const std::vector<unsigned int> col_index_set(
1626 for (
unsigned int r = 0; r < row_index_set.size(); ++r)
1627 for (
unsigned int c = 0; c < col_index_set.size(); ++c)
1630 out, r, c, hessian[row_index_set[r]][col_index_set[c]]);
1644 typename ScalarType>
1657 typename ScalarType>
1662 Assert(funcs.size() == this->n_dependent_variables(),
1664 "Vector size does not match number of dependent variables"));
1674 typename ScalarType>
1680 this->
taped_driver.keep_independent_values() ==
false) ||
1685 this->n_independent_variables(),
1687 "Not all values of sensitivities have been registered or subsequently set!"));
1691 ExcMessage(
"Not all dependent variables have been registered."));
1696 if (values.size() != this->n_dependent_variables())
1707 "Cannot compute values while tape is being recorded."));
1709 this->n_independent_variables(),
1711 this->n_independent_variables()));
1730 typename ScalarType>
1736 this->
taped_driver.keep_independent_values() ==
false) ||
1741 this->n_independent_variables(),
1743 "Not all values of sensitivities have been registered or subsequently set!"));
1747 ExcMessage(
"Not all dependent variables have been registered."));
1752 if (jacobian.
m() != this->n_dependent_variables() ||
1753 jacobian.
n() != this->n_independent_variables())
1755 this->n_independent_variables()},
1765 "Cannot compute Jacobian while tape is being recorded."));
1767 this->n_independent_variables(),
1769 this->n_independent_variables()));
1781 this->n_independent_variables(),
1783 this->n_independent_variables()));
1790 for (
unsigned int j = 0; j < this->n_independent_variables(); ++j)
1797 jacobian[i][j] *= 0.5;
1805 typename ScalarType>
1825 const std::vector<unsigned int> row_index_set(
1827 const std::vector<unsigned int> col_index_set(
1834 jacobian[row_index_set[0]][col_index_set[0]]);
1843 typename ScalarType>
1863 const std::vector<unsigned int> row_index_set(
1865 const std::vector<unsigned int> col_index_set(
1868 for (
unsigned int r = 0; r < row_index_set.size(); ++r)
1869 for (
unsigned int c = 0; c < col_index_set.size(); ++c)
1872 out, r, c, jacobian[row_index_set[r]][col_index_set[c]]);
1884# include "ad_helpers.inst"
1886# ifdef DEAL_II_WITH_ADOLC
1887# include "ad_helpers.inst1"
1889# ifdef DEAL_II_TRILINOS_WITH_SACADO
1890# include "ad_helpers.inst2"
const std::vector< ad_type > & get_sensitive_dof_values() const
void set_dof_values(const std::vector< scalar_type > &dof_values)
void register_dof_values(const std::vector< scalar_type > &dof_values)
typename HelperBase< ADNumberTypeCode, ScalarType >::ad_type ad_type
CellLevelBase(const unsigned int n_independent_variables, const unsigned int n_dependent_variables)
void register_energy_functional(const ad_type &energy)
EnergyFunctional(const unsigned int n_independent_variables)
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
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
typename AD::NumberTraits< ScalarType, ADNumberTypeCode >::ad_type ad_type
void stop_recording_operations(const bool write_tapes_to_file=false)
std::vector< bool > registered_marked_independent_variables
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
typename AD::NumberTraits< ScalarType, ADNumberTypeCode >::scalar_type scalar_type
void set_sensitivity_value(const unsigned int index, const scalar_type &value)
bool active_tape_requires_retaping() const
std::vector< ad_type > independent_variables
HelperBase(const unsigned int n_independent_variables, const unsigned int n_dependent_variables)
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
void activate_recorded_tape(const typename Types< ad_type >::tape_index tape_index)
bool is_symmetric_independent_variable(const unsigned int index) const
PointLevelFunctionsBase(const unsigned int n_independent_variables, const unsigned int n_dependent_variables)
void set_independent_variables(const std::vector< scalar_type > &values)
std::vector< bool > symmetric_independent_variables
unsigned int n_symmetric_independent_variables() const
void register_independent_variables(const std::vector< scalar_type > &values)
typename HelperBase< ADNumberTypeCode, ScalarType >::scalar_type scalar_type
void set_sensitivity_value(const unsigned int index, const bool symmetric_component, const scalar_type &value)
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
ResidualLinearization(const unsigned int n_independent_variables, const unsigned int n_dependent_variables)
virtual void compute_residual(Vector< scalar_type > &residual) const override
virtual void compute_linearization(FullMatrix< scalar_type > &linearization) const override
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)
void register_dependent_variable(const ad_type &func)
typename HelperBase< ADNumberTypeCode, ScalarType >::scalar_type scalar_type
ScalarFunction(const unsigned int n_independent_variables)
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 >::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)
VectorFunction(const unsigned int n_independent_variables, const unsigned int n_dependent_variables)
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
#define Assert(cond, exc)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
const bool IsBlockVector< VectorType >::value
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)
static const unsigned int invalid_unsigned_int
static const Types< ADNumberType >::tape_index max_tape_index
static const Types< ADNumberType >::tape_index invalid_tape_index
static void initialize_global_environment(const unsigned int n_independent_variables)
unsigned int tape_buffer_sizes
static constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE const T & value(const T &t)