From 0299c743a5faf89dc65ef50a414a4ca8b1f05025 Mon Sep 17 00:00:00 2001 From: Jean-Paul Pelteret Date: Wed, 25 Jul 2018 18:55:13 +0200 Subject: [PATCH] AD Helpers: Introduce base class for AD helper classes --- .../automatic_and_symbolic_differentiation.h | 2 + .../deal.II/differentiation/ad/ad_helpers.h | 723 +++++++++++++++++ source/differentiation/ad/CMakeLists.txt | 3 + source/differentiation/ad/ad_helpers.cc | 743 ++++++++++++++++++ source/differentiation/ad/ad_helpers.inst1.in | 55 ++ source/differentiation/ad/ad_helpers.inst2.in | 67 ++ 6 files changed, 1593 insertions(+) create mode 100644 include/deal.II/differentiation/ad/ad_helpers.h create mode 100644 source/differentiation/ad/ad_helpers.cc create mode 100644 source/differentiation/ad/ad_helpers.inst1.in create mode 100644 source/differentiation/ad/ad_helpers.inst2.in diff --git a/doc/doxygen/headers/automatic_and_symbolic_differentiation.h b/doc/doxygen/headers/automatic_and_symbolic_differentiation.h index 15ea93c26e..ae5178e2c8 100644 --- a/doc/doxygen/headers/automatic_and_symbolic_differentiation.h +++ b/doc/doxygen/headers/automatic_and_symbolic_differentiation.h @@ -411,6 +411,8 @@ * - ad_drivers.h: Provides classes that act as drivers to the interface of internally supported * automatic differentiation libraries. These are used internally as an intermediary to the * helper classes that we provide. + * - ad_helpers.h: Provides a set of classes to help perform automatic differentiation in a + * number of different contexts. These are detailed in @ref auto_diff_1_3. * - ad_number_types.h: Introduces an enumeration (called a type code) for the * auto-differentiable number combinations that will be supported by the driver classes. * The rationale behind the use of this somewhat restrictive mechanism is discussed below. diff --git a/include/deal.II/differentiation/ad/ad_helpers.h b/include/deal.II/differentiation/ad/ad_helpers.h new file mode 100644 index 0000000000..a74c699b94 --- /dev/null +++ b/include/deal.II/differentiation/ad/ad_helpers.h @@ -0,0 +1,723 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2017 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +#ifndef dealii_differentiation_ad_ad_helpers_h +#define dealii_differentiation_ad_ad_helpers_h + +#include + +#if defined(DEAL_II_WITH_ADOLC) || defined(DEAL_II_TRILINOS_WITH_SACADO) + +# include + +# include +# include +# include +# include +# include +# include + +# include +# include +# include +# include +# include + +DEAL_II_NAMESPACE_OPEN + +namespace Differentiation +{ + namespace AD + { + /** + * A base helper class that facilitates the evaluation of the derivative(s) + * of a number of user-defined dependent variables $\mathbf{f}(\mathbf{X})$ + * with respect to a set of independent variables $\mathbf{X}$, that is + * $\dfrac{d^{i} \mathbf{f}(\mathbf{X})}{d \mathbf{X}^{i}}$. + * + * In addition to the dimension @p dim, this class is templated on the + * floating point type @p scalar_type of the number that we'd like to + * differentiate, as well as an enumeration indicating the @p ADNumberTypeCode . + * The @p ADNumberTypeCode dictates which auto-differentiation library is + * to be used, and what the nature of the underlying auto-differentiable + * number is. Refer to the @ref auto_symb_diff module for more details in + * this regard. + * + * For all of the classes derived from this base class, there are two + * possible ways that the code in which they are used can be structured. + * The one approach is effectively a subset of the other, and which might + * be necessary to use depends on the nature of the chosen + * auto-differentiable number. + * + * When "tapeless" numbers are employed, the most simple code structure + * would be of the following form: + * + * @code + * // Initialize AD helper + * ADHelperType ad_helper (...); + * + * // Register independent variables + * ad_helper.register_independent_variable(...); + * + * // Extract the sensitive equivalent of the independent variables. + * // They are the auto-differentiable counterparts to the values + * // used as arguments to the register_independent_variable() function. + * // The operations conducted with these AD numbers will be tracked. + * const auto ad_independent_variables + * = ad_helper.get_sensitive_variables(...); + * + * // Use the sensitive variables to compute the dependent variables. + * const auto ad_dependent_variables = func(sensitive_variables); + * + * // Register the dependent variables with the helper class + * ad_helper.register_dependent_variables(ad_dependent_variables); + * + * // Compute derivatives of the dependent variables + * const auto derivatives = ad_helper.compute_gradients(); + * @endcode + * + * Note that since the specialized classes interpret the independent + * variables in different ways, above represents only an outline of the + * steps taken to compute derivatives. More specific examples are outlined + * in the individual classes that specialize this base class. + * + * When "taped" numbers are to be used, the above code should be wrapped by + * a few more lines of code to manage the taping procedure: + * + * @code + * // Initialize AD helper + * ADHelperType ad_helper (...); + * + * // An optional call to set the amount of memory to be allocated to + * // storing taped data + * ad_helper.set_tape_buffer_sizes(); + * + * // 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); + * if (is_recording == true) + * { + * // This is the "recording" phase of the operations. + * // In this block one places the majority of the operations described + * // in the previous code block. The set of operations that are + * // conducted here therefore includes the following steps: + * // - Register independent variables + * // - Extract the sensitive equivalent of the independent variables + * // - Use the sensitive variables to compute the dependent variables + * // - Register the dependent variables with the helper class + * + * // Indicate that 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 revaluate 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_no); + * + * // Set the new values of the independent variables where the recorded + * // dependent functions are to be evaluated (and differentiated + * // around). + * ad_helper.set_independent_variable(...); + * } + * + * // Compute derivatives of the dependent variables + * const auto derivatives = ad_helper.compute_gradients(); + * @endcode + * + * The second approach outlined here is more general than the first, and + * will work equally well for both taped and tapeless auto-differentiable + * numbers. + * + * @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. + * + * @todo Make this class thread safe for Sacado number and ADOL-C tapeless + * numbers (if supported). + * + * @author Jean-Paul Pelteret, 2016, 2017, 2018 + */ + template + class ADHelperBase + { + public: + /** + * Type definition for the floating point number type that is used in, + * and result from, all computations. + */ + using scalar_type = + typename AD::NumberTraits::scalar_type; + + /** + * Type definition for the auto-differentiation number type that is used + * in all computations. + */ + using ad_type = + typename AD::NumberTraits::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{f}(\mathbf{X})$, this will be the number of inputs + * $\mathbf{X}$, i.e. the dimension of the range 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{f}(\mathbf{X})$, this will + * be the number of outputs $\mathbf{f}$, i.e. the dimension of the + * domain or image space. + */ + ADHelperBase(const unsigned int n_independent_variables, + const unsigned int n_dependent_variables); + + /** + * Destructor + */ + virtual ~ADHelperBase(); + + //@} + + /** + * @name Interrogation of internal information + */ + //@{ + + /** + * Returns the number of independent variables that this object expects to + * work with. This is the dimension of the range space. + */ + std::size_t + n_independent_variables() const; + + /** + * Returns the number of dependent variables that this object expects to + * operate on. This is the dimension of the domain or image space. + */ + std::size_t + n_dependent_variables() const; + + /** + * Prints the status of all queriable data. Exactly what is printed and + * its format depends on the @p ad_type, as is determined by the + * @p ADNumberTypeCode template parameter. + * + * @param[in] stream The output stream to which the values are to be + * written. + */ + void + print(std::ostream &stream) const; + + /** + * Prints the values currently assigned to the independent variables. + * + * @param[in] stream The output stream to which the values are to be + * written. + */ + void + print_values(std::ostream &stream) const; + + /** + * Prints the statistics regarding the usage of the tapes. + * + * @param[in] tape_index The index of the tape to get the statistics of. + * @param[out] stream The output stream to which the values are to be + * written. + * + * @note This function only produces meaningful output when @p ad_type + * is a taped auto-differentiable number. + */ + void + print_tape_stats(const types::tape_index tape_index, + std::ostream & stream) const; + + //@} + + /** + * @name Operations specific to tapeless mode + */ + //@{ + + /** + * Pre-specify the number of @p independent_variables to be used in + * tapeless mode. + * + * Although this function is called internally in the ADHelperBase + * constructor, there may be occasions when ADOL-C tapeless numbers + * (adtl::adoubles) are created before an instance of this class + * is created. This function therefore allows one to declare at the + * earliest possible instance how many directional derivatives will be + * considered in tapeless mode. + * + * @warning With @p ensure_persistent_setting set to true when + * the @p ad_type is an ADOL-C tapeless number, calling this function + * leaves the set number of directional derivatives in a persistent state. + * It will therefore not be possible to further modify the number of + * directional derivatives to be tracked by adtl::adoubles's + * during course of the program's execution. + */ + static void + configure_tapeless_mode(const unsigned int n_independent_variables, + const bool ensure_persistent_setting = true); + + //@} + + /** + * @name Operations specific to taped mode: Recording tapes + */ + //@{ + + /** + * Reset the state of the helper class. + * + * When an instance of an ADHelperBase is stored as a class member object + * with the intention to reuse its instance, it may be necessary to reset + * the state of the object before use. This is because, internally, there + * is error checking performed to ensure that the correct + * auto-differentiable data is being tracked and used only when + * appropriate. This function clears all member data and, therefore, + * allows the state of all internal flags to be safely reset to their + * initial state. + * + * In the rare case that the number of independent or dependent variables + * has changed, this can also reconfigured by passing in the appropriate + * arguments to the function. + * + * @note This also resets the active tape number to an invalid number, and + * deactivates the recording mode for taped variables. + */ + virtual void + reset(const unsigned int n_independent_variables = + dealii::numbers::invalid_unsigned_int, + const unsigned int n_dependent_variables = + dealii::numbers::invalid_unsigned_int, + const bool clear_registered_tapes = true); + + /** + * Returns whether or not this class is tracking calculations performed + * with its marked independent variables. + */ + bool + is_recording() const; + + /** + * Returns the tape number which is currently activated for recording or + * reading. + */ + types::tape_index + active_tape() const; + + /** + * Returns whether or not a tape number has already been used + * or registered. + */ + bool + is_registered_tape(const types::tape_index tape_index) const; + + /** + * Set the buffer sizes for the next active tape. + * + * This function must be called before start_recording_operations() for + * it to have any influence on the memory allocated to the next recorded + * tape. + * + * @note This function only has an effect when using ADOL-C numbers. As + * stated by the ADOL-C manual, it may be desirable to create a file + * ".adolcrc" in the program run directory and set the buffer size + * therein. + * Alternatively, this function can be used to override the settings for + * any given tape, or can be used in the event that no ".adolcrc" file + * exists. + * The default value for each buffer is set at 64 megabytes, a + * heuristically chosen value thought to be appropriate for use within + * the context of finite element analysis when considering coupled + * problems with multiple vector-valued fields discretised by higher + * order shape functions, as well as complex constitutive laws. + * + * @param[in] obufsize ADOL-C operations buffer size + * @param[in] lbufsize ADOL-C locations buffer size + * @param[in] vbufsize ADOL-C value buffer size + * @param[in] tbufsize ADOL-C Taylor buffer size + */ + void + set_tape_buffer_sizes(const types::tape_buffer_sizes obufsize = 67108864, + const types::tape_buffer_sizes lbufsize = 67108864, + const types::tape_buffer_sizes vbufsize = 67108864, + const types::tape_buffer_sizes tbufsize = 67108864); + + /** + * Enable recording mode for a given tape. The use of this function is + * mandatory if the auto-differentiable number is a taped type. However, + * for the purpose of developing generic code, it can also be safely + * called for tapeless auto-differentiable numbers. + * + * The operations that take place between this function call and that + * of stop_recording_operations() are recorded to the tape and can + * be replayed and reevaluated as necessary. + * + * The typical set of operations to be performed during this "recording" + * phase (between the calls to start_recording_operations() and + * stop_recording_operations() ) are: + * - Definition of some independent variables via + * register_independent_variable() or + * register_independent_variables(). These define the branch of + * operations tracked by the tape. If the @p keep flag is set to + * true then these represent precisely the point about which + * the function derivatives are to be computed. If the @p keep flag is + * set to false then these only represent dummy values, and + * the point at which the function derivatives are to be computed must + * be set by calling set_independent_variables() again. + * - Extraction of a set of independent variables of auto-differentiable + * type using get_sensitive_variables(). These are then tracked during + * later computations. + * - Defining the dependent variables via register_dependent_variable() + * or register_dependent_variables(). These are the functions that + * will be differentiated with respect to the independent variables. + * + * @param[in] tape_index The index of the tape to be written + * @param[in] overwrite_tape Express whether tapes are allowed to be + * overwritten. If true then any existing tape with a given + * @p tape_index will be destroyed and a new tape traced over it. + * @param[in] keep_values Determines whether the numerical values of all + * independent variables are recorded in the tape buffer. If true, then + * the tape can be immediately used to perform computations after + * recording is complete. + * + * @note During the recording phase, no value(), gradient(), hessian(), + * or jacobian() operations can be performed. + * + * @note The chosen tape index must be greater than + * numbers::invalid_tape_index and less than numbers::max_tape_index. + */ + bool + start_recording_operations(const types::tape_index tape_index, + const bool overwrite_tape = false, + const bool keep_values = true); + + /** + * Disable recording mode for a given tape. The use of this function is + * mandatory if the auto-differentiable number is a taped type. However, + * for the purpose of developing generic code, it can also be safely + * called for tapeless auto-differentiable numbers. + * + * @note After this function call, the tape is considered ready for use and + * operations such as value(), gradient() or hessian() can be executed. + * + * @note For taped AD numbers, this operation is only valid in recording mode. + */ + void + stop_recording_operations(const bool write_tapes_to_file = false); + + /** + * Select a pre-recorded tape to read from. + * + * @param[in] tape_index The index of the tape to be read from. + * + * @note The chosen tape index must be greater than + * numbers::invalid_tape_index and less than numbers::max_tape_index. + */ + void + activate_recorded_tape(const types::tape_index tape_index); + + //@} + + protected: + /** + * @name Taping + */ + //@{ + + /** + * Index of the tape that is currently in use. It is this tape that will + * be recorded to or read from when performing computations using "taped" + * auto-differentiable numbers. + */ + types::tape_index active_tape_index; + + /** + * A collection of tapes that have been recorded to on this process. + * + * It is important to keep track of this so that one doesn't accidentally + * record over a tape (unless specifically instructed to) and that one + * doesn't try to use a tape that doesn't exist. + */ + static std::set registered_tapes; + + /** + * Mark whether we're going to inform taped data structures to retain + * the coefficients ("Taylors" in ADOL-C nomenclature) stored on the + * tape so that they can be evaluated again at a later stage. + */ + bool keep_values; + + /** + * Mark whether we're currently recording a tape. Dependent on the state + * of this flag, only a restricted set of operations are allowable. + */ + bool is_recording_flag; + + /** + * A flag indicating that we should preferentially use the user-defined + * taped buffer sizes as opposed to either the default values selected + * by the AD library (or, in the case of ADOL-C, defined in an + * ".adolcrc" file). + */ + bool use_stored_taped_buffer_sizes; + + /** + * ADOL-C operations buffer size. + */ + types::tape_buffer_sizes obufsize; + + /** + * ADOL-C locations buffer size. + */ + types::tape_buffer_sizes lbufsize; + + /** + * ADOL-C value buffer size. + */ + types::tape_buffer_sizes vbufsize; + + /** + * ADOL-C Taylor buffer size. + */ + types::tape_buffer_sizes tbufsize; + + /** + * Select a tape to record to or read from. + * + * This function activates a tape, but depending on whether @p read_mode + * is set, the tape is either taken as previously written to (and put + * into read-only mode), or cleared for (re-)taping. + * + * @param[in] tape_index The index of the tape to be written to/read + * from. + * @param[in] read_mode A flag that marks whether or not we expect to + * read data from a preexisting tape. + * + * @note The chosen tape index must be greater than + * numbers::invalid_tape_index and less than numbers::max_tape_index. + */ + void + activate_tape(const types::tape_index tape_index, const bool read_mode); + + //@} + + /** + * @name Independent variables + */ + //@{ + + /** + * A set of independent variables $\mathbf{X}$ that differentiation will + * be performed with respect to. + * + * The gradients and Hessians of dependent variables will be computed + * at these finite values. + */ + mutable std::vector independent_variable_values; + + /** + * A set of sensitive independent variables $\mathbf{X}$ that + * differentiation will be performed with respect to. + * + * The gradients and Hessians of dependent variables will be computed + * using these configured AD numbers. Note that only reverse-mode AD + * requires that the sensitive independent variables be stored. + */ + mutable std::vector independent_variables; + + /** + * A list of registered independent variables that have been manipulated + * for a given set of operations. + */ + std::vector registered_independent_variable_values; + + /** + * A list of registered independent variables that have been extracted and + * their sensitivities marked. + */ + mutable std::vector registered_marked_independent_variables; + + /** + * Reset the boolean vector @p registered_independent_variable_values that + * indicates which independent variables we've been manipulating for + * the current set of operations. + */ + void + reset_registered_independent_variables(); + + /** + * Set the actual value of the independent variable $X_{i}$. + * + * @param[in] index The index in the vector of independent variables. + * @param[in] value The value to set the index'd independent variable to. + */ + void + set_sensitivity_value(const unsigned int index, const scalar_type &value); + + /** + * Initialize an independent variable $X_{i}$ such that subsequent + * operations performed with it are tracked. + * + * @note Care must be taken to mark each independent variable only once. + * + * @note The order in which the independent variables are marked defines the + * order of all future internal operations. They must be manipulated in + * the same order as that in which they are first marked. If not then, + * for example, ADOL-C won't throw an error, but rather it might complain + * nonsensically during later computations or produce garbage results. + * + * @param[in] index The index in the vector of independent variables. + * @param[out] out An auto-differentiable number that is ready for use in + * computations. The operations that are performed with it are recorded on + * the tape and will be considered in the computation of dependent + * variable values. + */ + void + mark_independent_variable(const unsigned int index, ad_type &out) const; + + /** + * Finalize the state of the independent variables before use. + * + * This step and the storage of the independent variables is done + * separately because some derived classes may offer the capability + * to add independent variables in a staggered manner. This function + * is to be triggered when these values are considered finalized + * and we can safely initialize the sensitive equivalents of those + * values. + */ + void + finalize_sensitive_independent_variables() const; + + /** + * Initialize an independent variable $X_{i}$. + * + * @param[out] out An auto-differentiable number that is ready for use in + * standard computations. The operations that are performed with it are + * not recorded on the tape, and so should only be used when not in + * recording mode. + * @param[in] index The index in the vector of independent variables. + */ + void + initialize_non_sensitive_independent_variable(const unsigned int index, + ad_type &out) const; + + /** + * The number of independent variables that have been manipulated within a + * set of operations. + */ + unsigned int + n_registered_independent_variables() const; + //@} + + /** + * @name Dependent variables + */ + //@{ + + /** + * The set of dependent variables $\mathbf{f}(\mathbf{X})$ of which the + * derivatives with respect to $\mathbf{X}$ will be computed. + * + * The gradients and Hessians of these dependent variables will be + * computed at the values $\mathbf{X}$ set with the + * set_sensitivity_value() function. + * + * @note These are stored as an @p ad_type so that we can use them to + * compute function values and directional derivatives in the case that + * tapeless numbers are used + */ + std::vector dependent_variables; + + /** + * A list of registered dependent variables. + */ + std::vector registered_marked_dependent_variables; + + /** + * Reset the boolean vector @p registered_marked_dependent_variables that + * indicates which independent variables have been manipulated by the + * current set of operations. All entries in the vector are set to the + * value of the @p flag. + */ + void + reset_registered_dependent_variables(const bool flag = false); + + /** + * The number of dependent variables that have been registered. + */ + unsigned int + n_registered_dependent_variables() const; + + /** + * Register the definition of the index'th dependent variable + * $f(\mathbf{X})$. + * + * @param[in] index The index of the entry in the global list of dependent + * variables that this function belongs to. + * @param[in] func The recorded function that defines a dependent + * variable. + * + * @note Each dependent variable must only be registered once. + */ + void + register_dependent_variable(const unsigned int index, + const ad_type & func); + //@} + + private: + /** + * A counter keeping track of the number of helpers in existence. + * + * This is only important information for when we use taped number types. + * As the tapes are stored in some global register, they exist independent + * of these helpers. However, it is assumed that when all helpers go + * out of scope then the tapes can be written over. + */ + static unsigned int n_helpers; + + }; // class ADHelperBase + + + } // namespace AD +} // namespace Differentiation + + +DEAL_II_NAMESPACE_CLOSE + +#endif // defined(DEAL_II_WITH_ADOLC) || defined(DEAL_II_TRILINOS_WITH_SACADO) + +#endif // dealii__adolc_helpers_h diff --git a/source/differentiation/ad/CMakeLists.txt b/source/differentiation/ad/CMakeLists.txt index 10781942ed..81e83bb633 100644 --- a/source/differentiation/ad/CMakeLists.txt +++ b/source/differentiation/ad/CMakeLists.txt @@ -16,11 +16,14 @@ INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_BINARY_DIR}) SET(_src + ad_helpers.cc adolc_number_types.cc sacado_number_types.cc ) SET(_inst + ad_helpers.inst1.in + ad_helpers.inst2.in adolc_number_types.inst.in sacado_number_types.inst1.in sacado_number_types.inst2.in diff --git a/source/differentiation/ad/ad_helpers.cc b/source/differentiation/ad/ad_helpers.cc new file mode 100644 index 0000000000..db9ed5325f --- /dev/null +++ b/source/differentiation/ad/ad_helpers.cc @@ -0,0 +1,743 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2017 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +#include + +#if defined(DEAL_II_WITH_ADOLC) || defined(DEAL_II_TRILINOS_WITH_SACADO) + +# include +# include + +# include + + +DEAL_II_NAMESPACE_OPEN + + +namespace Differentiation +{ + namespace AD + { + /* -------------------------- ADHelperBase -------------------------- */ + + + + template + std::set + ADHelperBase::registered_tapes; + + + + template + unsigned int ADHelperBase::n_helpers = 0; + + + + template + ADHelperBase::ADHelperBase( + const unsigned int n_independent_variables, + const unsigned int n_dependent_variables) + : active_tape_index(numbers::invalid_tape_index) + , keep_values(true) + , is_recording_flag(false) + , use_stored_taped_buffer_sizes(false) + , obufsize(0u) + , lbufsize(0u) + , vbufsize(0u) + , tbufsize(0u) + , independent_variable_values( + n_independent_variables, + dealii::internal::NumberType::value(0.0)) + , registered_independent_variable_values(n_independent_variables, false) + , registered_marked_independent_variables(n_independent_variables, false) + , registered_marked_dependent_variables(n_dependent_variables, false) + { + ++n_helpers; + + // Tapeless mode must be configured before any active live + // variables are created. + if (AD::is_tapeless_ad_number::value) + { + configure_tapeless_mode(n_independent_variables, + false /*ensure_persistent_setting*/); + } + + // For safety, we ensure that the entries in this vector *really* are + // initialized correctly by sending in the constructed zero-valued + // initializer. + dependent_variables.resize(n_dependent_variables, + dealii::internal::NumberType::value( + 0.0)); + } + + + + template + ADHelperBase::~ADHelperBase() + { + --n_helpers; + + // Clear any static data when there are no more helpers in scope. + // This means that we've effectively marked any of these tapes to + // be safely overwritten. + if (n_helpers == 0) + registered_tapes.clear(); + } + + + + template + void + ADHelperBase::reset_registered_independent_variables() + { + for (typename std::vector::iterator it = + registered_independent_variable_values.begin(); + it != registered_independent_variable_values.end(); + ++it) + *it = false; + } + + + + template + void + ADHelperBase:: + reset_registered_dependent_variables(const bool flag) + { + for (typename std::vector::iterator it = + registered_marked_dependent_variables.begin(); + it != registered_marked_dependent_variables.end(); + ++it) + *it = flag; + } + + + + template + void + ADHelperBase::set_sensitivity_value( + const unsigned int index, + const scalar_type &value) + { + if (ADNumberTraits::is_tapeless == true) + { + // A dummy call in case the user does not encapsulate a set of + // calls to tapeless ADHelpers with the initial and final calls to + // [start,stop]_recording_operations. + if (this->is_recording() == false) + start_recording_operations(1 /*tape index*/); + + Assert(this->is_recording() == true, + ExcMessage( + "Cannot change the value of an independent variable " + "of the tapeless variety while this class is not set" + "in recording operations.")); + } + if (ADNumberTraits::is_taped == true) + { + Assert(this->active_tape() != numbers::invalid_tape_index, + ExcMessage("Invalid tape index")); + } + Assert( + index < n_independent_variables(), + ExcMessage( + "Trying to set the value of a non-existent independent variable.")); + + independent_variable_values[index] = value; + registered_independent_variable_values[index] = true; + } + + + + template + void + ADHelperBase::mark_independent_variable( + const unsigned int index, + ad_type & out) const + { + Assert(index < n_independent_variables(), ExcInternalError()); + Assert(registered_independent_variable_values[index] == true, + ExcInternalError()); + + if (index > 0) + { + Assert( + registered_marked_independent_variables[index - 1] == true, + ExcMessage( + "Need to extract sensitivities in the order they're created.")); + } + + if (ADNumberTraits::is_taped == true) + { + Assert(active_tape() != numbers::invalid_tape_index, + ExcMessage("Invalid tape index")); + Assert(is_recording() == true, + ExcMessage( + "The marking of independent variables is only valid " + "during recording.")); + } + + internal::Marking::independent_variable( + independent_variable_values[index], + index, + this->n_independent_variables(), + out); + registered_marked_independent_variables[index] = true; + } + + + + template + void + ADHelperBase::finalize_sensitive_independent_variables() const + { + // Double check that we've actually registered all DoFs + Assert(n_registered_independent_variables() == n_independent_variables(), + ExcMessage("Not all values of sensitivities have been recorded!")); + + // This should happen only once + if (this->independent_variables.size() == 0) + { + this->independent_variables.resize( + this->n_independent_variables(), + dealii::internal::NumberType::value(0.0)); + + // Indicate the sensitivity that each entry represents + for (unsigned int i = 0; i < this->n_independent_variables(); ++i) + this->mark_independent_variable(i, this->independent_variables[i]); + } + } + + + + template + void + ADHelperBase:: + initialize_non_sensitive_independent_variable(const unsigned int index, + ad_type &out) const + { + if (ADNumberTraits::is_taped == true) + { + Assert(active_tape() != numbers::invalid_tape_index, + ExcMessage("Invalid tape index")); + } + Assert(is_recording() == false, + ExcMessage( + "The initialization of non-sensitive independent variables is " + "only valid outside of recording operations.")); + + Assert(index < n_independent_variables(), ExcInternalError()); + Assert(registered_independent_variable_values[index] == true, + ExcInternalError()); + + out = independent_variable_values[index]; + } + + + + template + unsigned int + ADHelperBase::n_registered_independent_variables() const + { + return std::count(registered_independent_variable_values.begin(), + registered_independent_variable_values.end(), + true); + } + + + + template + std::size_t + ADHelperBase::n_independent_variables() const + { + return independent_variable_values.size(); + } + + + + template + unsigned int + ADHelperBase::n_registered_dependent_variables() const + { + return std::count(registered_marked_dependent_variables.begin(), + registered_marked_dependent_variables.end(), + true); + } + + + + template + std::size_t + ADHelperBase::n_dependent_variables() const + { + return dependent_variables.size(); + } + + + + template + bool + ADHelperBase::is_recording() const + { + return is_recording_flag; + } + + + + template + types::tape_index + ADHelperBase::active_tape() const + { + return active_tape_index; + } + + + + template + bool + ADHelperBase::is_registered_tape( + const types::tape_index tape_index) const + { + return registered_tapes.find(tape_index) != registered_tapes.end(); + } + + + + template + void + ADHelperBase::print( + std::ostream &stream) const + { + // Store stream flags + const std::ios_base::fmtflags stream_flags(stream.flags()); + // Set stream to print booleans as "true"/"false" + stream.setf(std::ios_base::boolalpha); + + stream << "Active tape index: " << active_tape_index << "\n"; + stream << "Registered tapes: "; + auto it_registered_tape = registered_tapes.begin(); + for (unsigned int i = 0; i < registered_tapes.size(); + ++i, ++it_registered_tape) + stream << *it_registered_tape + << (i < (registered_tapes.size() - 1) ? "," : ""); + stream << "\n"; + stream << "Recording? " << is_recording() << "\n"; + stream << "Keep values? " << keep_values << "\n"; + stream << "Use stored tape buffer sizes? " + << use_stored_taped_buffer_sizes << "\n"; + stream << std::flush; + + stream << "Registered independent variables: " + << "\n"; + for (unsigned int i = 0; i < n_independent_variables(); i++) + stream << registered_independent_variable_values[i] + << (i < (n_independent_variables() - 1) ? "," : ""); + stream << std::endl; + + stream << "Independent variable values: " + << "\n"; + print_values(stream); + + stream << "Registered marked independent variables: " + << "\n"; + for (unsigned int i = 0; i < n_independent_variables(); i++) + stream << registered_marked_independent_variables[i] + << (i < (n_independent_variables() - 1) ? "," : "") + << std::flush; + stream << std::endl; + + stream << "Dependent variable values: " + << "\n"; + for (unsigned int i = 0; i < n_dependent_variables(); i++) + stream << dependent_variables[i] + << (i < (n_dependent_variables() - 1) ? "," : ""); + stream << std::endl; + + stream << "Registered dependent variables: " + << "\n"; + for (unsigned int i = 0; i < n_dependent_variables(); i++) + stream << registered_marked_dependent_variables[i] + << (i < (n_dependent_variables() - 1) ? "," : ""); + stream << std::endl; + + // Restore stream flags + stream.flags(stream_flags); + } + + + + template + void + ADHelperBase::print_values( + std::ostream &stream) const + { + for (unsigned int i = 0; i < n_independent_variables(); i++) + stream << independent_variable_values[i] + << (i < (n_independent_variables() - 1) ? "," : "") + << std::flush; + stream << std::endl; + } + + + + template + void + ADHelperBase::print_tape_stats( + const types::tape_index tape_index, + std::ostream & stream) const + { + if (ADNumberTraits::is_tapeless == true) + return; + + Assert(is_registered_tape(tape_index), + ExcMessage("Tape number not registered")); + + TapedDrivers::print_tape_stats(tape_index, stream); + } + + + + template + void + ADHelperBase::reset( + const unsigned int n_independent_variables, + const unsigned int n_dependent_variables, + const bool clear_registered_tapes) + { + const unsigned int new_n_independent_variables = + (n_independent_variables != dealii::numbers::invalid_unsigned_int ? + n_independent_variables : + this->n_independent_variables()); + const unsigned int new_n_dependent_variables = + (n_dependent_variables != dealii::numbers::invalid_unsigned_int ? + n_dependent_variables : + this->n_dependent_variables()); + + // Here we clear our vectors of AD data with a reallocation of memory + // (i.e. we *nuke* them entirely). Why we do this differs for each + // AD type: + // - ADOL-C taped mode must have their tapes fully cleared of marked data + // before the tapes can be overwritten. + // - ADOL-C tapeless mode must be configured for before any active live + // variables are created. + // - Reverse-mode Sacado numbers to must be destroyed to reset their + // accumulations. + // - Forward-mode Sacado numbers have no specific requirements, but it + // doesn't really hurt to perform this operation anyway. + { + std::vector().swap(independent_variables); + std::vector().swap(dependent_variables); + } + + // Tapeless mode must be configured before any active live + // variables are created. + if (AD::is_tapeless_ad_number::value) + { + configure_tapeless_mode(n_independent_variables, + false /*ensure_persistent_setting*/); + } + + active_tape_index = numbers::invalid_tape_index; + is_recording_flag = false; + if (clear_registered_tapes) + registered_tapes.clear(); + + independent_variable_values = std::vector( + new_n_independent_variables, + dealii::internal::NumberType::value(0.0)); + registered_independent_variable_values = + std::vector(new_n_independent_variables, false); + registered_marked_independent_variables = + std::vector(new_n_independent_variables, false); + dependent_variables = + std::vector(new_n_dependent_variables, + dealii::internal::NumberType::value(0.0)); + registered_marked_dependent_variables = + std::vector(new_n_dependent_variables, false); + } + + + + template + void + ADHelperBase::configure_tapeless_mode( + const unsigned int n_independent_variables, + const bool ensure_persistent_setting) + { + if (ADNumberTraits::is_taped == true) + return; + + // Try to safely initialize the global environment + TapelessDrivers::initialize_global_environment( + n_independent_variables); + + if (ensure_persistent_setting == true) + if (is_adolc_tapeless_number::value) + { + // In order to ensure that the settings remain for the entire + // duration of the simulation, we create a global live variable + // that doesn't go out of scope. + static ad_type num = 0.0; + } + } + + + + template + void + ADHelperBase::activate_recorded_tape( + const types::tape_index tape_index) + { + activate_tape(tape_index, true /*read_mode*/); + } + + + + template + void + ADHelperBase::activate_tape( + const types::tape_index tape_index, + const bool read_mode) + { + if (ADNumberTraits::is_taped == true) + { + Assert(tape_index != numbers::invalid_tape_index, + ExcMessage("Invalid tape index")); + Assert(tape_index < numbers::max_tape_index, + ExcMessage("Tape index exceeds maximum allowable value")); + Assert(is_registered_tape(tape_index), + ExcMessage("Tape number not registered")); + active_tape_index = tape_index; + reset_registered_independent_variables(); + + // A tape may have been defined by a different ADHelper, so in this + // case we ignore the fact that any dependent variables within the + // current data structure have not been marked as dependents + if (read_mode == true) + { + reset_registered_dependent_variables(true); + Assert(n_registered_dependent_variables() == + n_dependent_variables(), + ExcMessage("Not all dependent variables have been set!")); + } + } + else + { + Assert(ADNumberTraits::is_tapeless == true, + ExcInternalError()); + active_tape_index = tape_index; // This is, in effect, a dummy value + } + } + + + + template + void + ADHelperBase::set_tape_buffer_sizes( + const types::tape_buffer_sizes in_obufsize, + const types::tape_buffer_sizes in_lbufsize, + const types::tape_buffer_sizes in_vbufsize, + const types::tape_buffer_sizes in_tbufsize) + { + // When valid for the chosen AD number type, these values will be used the + // next time start_recording_operations() is called. + obufsize = in_obufsize; + lbufsize = in_lbufsize; + vbufsize = in_vbufsize; + tbufsize = in_tbufsize; + use_stored_taped_buffer_sizes = true; + } + + + + template + bool + ADHelperBase::start_recording_operations( + const types::tape_index tape_index, + const bool overwrite_tape, + const bool keep) + { + // Define this here for clarity when this flag is used later. + const bool read_mode = false; + + if (ADNumberTraits::is_taped == true) + { + if (overwrite_tape != true) + { + Assert(is_recording() == false, + ExcMessage("Already recording...")); + } + if (is_registered_tape(tape_index) == false || overwrite_tape == true) + { + // Register the tape... + registered_tapes.insert(tape_index); + // ... and setup the data structures for this class in the + // appropriate manner + activate_tape(tape_index, read_mode); + + // Start taping + if (use_stored_taped_buffer_sizes) + { + Assert(obufsize > 0u, + ExcMessage("Buffer value not initialized.")); + Assert(lbufsize > 0u, + ExcMessage("Buffer value not initialized.")); + Assert(vbufsize > 0u, + ExcMessage("Buffer value not initialized.")); + Assert(tbufsize > 0u, + ExcMessage("Buffer value not initialized.")); + TapedDrivers::start_taping( + active_tape(), + keep, + obufsize, + lbufsize, + vbufsize, + tbufsize); + + // Reset this, as we don't assume that the same tape buffer + // sizes are used for every tape + use_stored_taped_buffer_sizes = false; + obufsize = 0u; + lbufsize = 0u; + vbufsize = 0u; + tbufsize = 0u; + } + else + { + TapedDrivers::start_taping( + active_tape(), keep); + } + + // Clear the flags that state which independent and + // dependent variables have been registered + reset_registered_independent_variables(); + reset_registered_dependent_variables(); + + // Set some other flags to their indicated / required values + keep_values = keep; + is_recording_flag = true; + } + } + else + { + Assert(ADNumberTraits::is_tapeless == true, + ExcInternalError()); + is_recording_flag = true; + + // Dummy call to ensure that the intuitively correct + // value for the active tape (whether "valid" or not) + // is always returned to the user. + activate_tape(tape_index, read_mode); + } + + return is_recording_flag; + } + + + + template + void + ADHelperBase::stop_recording_operations( + const bool write_tapes_to_file) + { + Assert(is_recording() == true, ExcMessage("Not currently recording...")); + + // Double check that we've actually registered all DoFs + Assert(n_registered_independent_variables() == n_independent_variables(), + ExcMessage("Not all values of sensitivities have been recorded!")); + + if (ADNumberTraits::is_taped == true) + { + // Stop taping + TapedDrivers::stop_taping(active_tape(), + write_tapes_to_file); + + // Now that we've turned tracing off, we've definitely + // stopped all tape recording. + is_recording_flag = false; + + // If the keep_values flag is set, then we expect the user to use this + // tape immediately after recording it. There is therefore no need to + // invalidate it. However, there is now also no way to double-check + // that the newly recorded tape is indeed the active tape. + if (keep_values == false) + active_tape_index = numbers::invalid_tape_index; + } + else + { + Assert(ADNumberTraits::is_tapeless == true, + ExcInternalError()); + // Double check that we've actually registered dependent variables + Assert(n_registered_dependent_variables() == n_dependent_variables(), + ExcMessage("Not all dependent variables have been set!")); + + // By changing this flag, we ensure that the we can no longer + // legally alter the values of the dependent variables using + // set_independent_variable(). This is important because the value of + // the tapeless independent variables are set and finalized when + // mark_independent_variable() is called. So we cannot allow this to + // be done when not in the "recording" phase + is_recording_flag = false; + active_tape_index = numbers::invalid_tape_index; + } + } + + + + template + void + ADHelperBase::register_dependent_variable( + const unsigned int index, + const ad_type & func) + { + Assert(index < n_dependent_variables(), ExcMessage("Index out of range")); + Assert(registered_marked_dependent_variables[index] == false, + ExcMessage( + "This dependent variable has already been registered.")); + + if (ADNumberTraits::is_taped == true) + { + Assert(active_tape() != numbers::invalid_tape_index, + ExcMessage("Invalid tape index")); + Assert(is_recording() == true, + ExcMessage( + "Must be recording when registering dependent variables.")); + } + + // Register the given dependent variable + internal::Marking::dependent_variable(dependent_variables[index], + func); + registered_marked_dependent_variables[index] = true; + } + } // namespace AD +} // namespace Differentiation + + +/* --- Explicit instantiations --- */ +# ifdef DEAL_II_WITH_ADOLC +# include "ad_helpers.inst1" +# endif +# ifdef DEAL_II_TRILINOS_WITH_SACADO +# include "ad_helpers.inst2" +# endif + + +DEAL_II_NAMESPACE_CLOSE + +#endif // defined(DEAL_II_WITH_ADOLC) || defined(DEAL_II_TRILINOS_WITH_SACADO) diff --git a/source/differentiation/ad/ad_helpers.inst1.in b/source/differentiation/ad/ad_helpers.inst1.in new file mode 100644 index 0000000000..b91a64ed0b --- /dev/null +++ b/source/differentiation/ad/ad_helpers.inst1.in @@ -0,0 +1,55 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2017 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +// TODO: Include complex types + +for (number : REAL_SCALARS) +{ + namespace Differentiation + \{ + namespace AD + \{ + // -------------------------- ADHelperBase ---------------------- + + template + class ADHelperBase; + + template + class ADHelperBase; + + \} + \} +} + +// Instantiations for ADHelpers for which the underlying number type is fixed +for () +{ + namespace Differentiation + \{ + namespace AD + \{ + + // -------------------------- ADHelperBase ---------------------- + + template + class ADHelperBase::ad_type>; + + template + class ADHelperBase::ad_type>; + + + \} + \} +} diff --git a/source/differentiation/ad/ad_helpers.inst2.in b/source/differentiation/ad/ad_helpers.inst2.in new file mode 100644 index 0000000000..793f417da1 --- /dev/null +++ b/source/differentiation/ad/ad_helpers.inst2.in @@ -0,0 +1,67 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2017 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +// TODO: Include complex types + +for (number : REAL_SCALARS) +{ + namespace Differentiation + \{ + namespace AD + \{ + // -------------------------- ADHelperBase ---------------------- + + template + class ADHelperBase; + + template + class ADHelperBase; + + template + class ADHelperBase; + + template + class ADHelperBase; + + \} + \} +} + +// Instantiations for ADHelpers for which the underlying number type is fixed +for () +{ + namespace Differentiation + \{ + namespace AD + \{ + + // -------------------------- ADHelperBase ---------------------- + + template + class ADHelperBase::ad_type>; + + template + class ADHelperBase::ad_type>; + + template + class ADHelperBase::ad_type>; + + template + class ADHelperBase::ad_type>; + + + \} + \} +} -- 2.39.5