# include <deal.II/differentiation/ad/sacado_number_types.h>
# include <deal.II/differentiation/ad/sacado_product_types.h>
+# include <deal.II/lac/full_matrix.h>
+# include <deal.II/lac/vector.h>
+
# include <algorithm>
# include <iostream>
# include <iterator>
}; // class ADHelperBase
+
+ /**
+ * A general helper class that facilitates the evaluation of a vector of
+ * functions, as well as its first derivatives (their Jacobian).
+ * This class would typically be used to compute the linearization of a
+ * set of local nonlinear equations, but can also be used as the basis of
+ * the linearization of the residual vector defined on the level of a finite
+ * element (for example, in order to compute the Jacobian matrix necessary
+ * in Newton-type solvers for nonlinear problems).
+ *
+ * @note When using the cell-level taped AD methods in 3d and/or with higher
+ * order elements, it is incredibly easy to exceed the tape buffer size.
+ * The reason for this is two-fold:
+ * 1. there are are many independent variables (the local
+ * degrees-of-freedom) to take the derivatives with respect to, and
+ * 2. the expressions for the dependent variables (each being a component
+ * of the residual vector) in terms of all of the independent variables
+ * are lengthy, especially when non-trivial constitutive laws are
+ * considered.
+ * These buffer variables dictate the amount of memory allocated to a tape
+ * before it is written to file (at a significant performance loss).
+ * Therefore for ADOL-C taped AD numbers, it may be desirable to
+ * create a file ".adolcrc" in the program run directory and set the buffer
+ * size therein (as is suggested by the ADOL-C manual). For example, the
+ * following settings increase the default buffer size by 128 times:
+ * @code
+ * "OBUFSIZE" "67108864"
+ * "LBUFSIZE" "67108864"
+ * "VBUFSIZE" "67108864"
+ * "TBUFSIZE" "67108864"
+ * @endcode
+ * Note that the quotation marks are mandatory.
+ * An alternative approach that allows for run-time decision making is to
+ * use the ADHelperBase::set_tape_buffer_sizes() function before starting
+ * taping (as done via the ADHelperBase::start_recording_operations()
+ * function).
+ *
+ * @warning ADOL-C does not support the standard threading models used by
+ * deal.II, so this class should @b not be embedded within a multithreaded
+ * function when using ADOL-C number types. It is, however, suitable for use
+ * in both serial and MPI routines.
+ *
+ * @author Jean-Paul Pelteret, 2016, 2017, 2018
+ */
+ template <enum AD::NumberTypes ADNumberTypeCode,
+ typename ScalarType = double>
+ class ADHelperCellLevelBase
+ : public ADHelperBase<ADNumberTypeCode, ScalarType>
+ {
+ public:
+ /**
+ * Type definition for the floating point number type that is used in,
+ * and results from, all computations.
+ */
+ using scalar_type =
+ typename ADHelperBase<ADNumberTypeCode, ScalarType>::scalar_type;
+
+ /**
+ * Type definition for the auto-differentiation number type that is used
+ * in all computations.
+ */
+ using ad_type =
+ typename ADHelperBase<ADNumberTypeCode, ScalarType>::ad_type;
+
+ /**
+ * @name Constructor / destructor
+ */
+ //@{
+
+ /**
+ * The constructor for the class.
+ *
+ * @param[in] n_independent_variables The number of independent variables
+ * that will be used in the definition of the functions that it is
+ * desired to compute the sensitivities of. In the computation of
+ * $\mathbf{f}(\mathbf{X})$, this will be the number of inputs
+ * $\mathbf{X}$, i.e. the dimension of the domain space.
+ * @param[in] n_dependent_variables The number of scalar functions to be
+ * defined that will have a sensitivity to the given independent
+ * variables. In the computation of $\mathbf{f}(\mathbf{X})$, this will
+ * be the number of outputs $\mathbf{f}$, i.e. the dimension of the
+ * image space.
+ */
+ ADHelperCellLevelBase(const unsigned int n_independent_variables,
+ const unsigned int n_dependent_variables);
+
+ /**
+ * Destructor
+ */
+ virtual ~ADHelperCellLevelBase() = default;
+
+ //@}
+
+ /**
+ * @name Independent variables
+ */
+ //@{
+
+ /**
+ * Register the complete set of independent variables $\mathbf{X}$ that
+ * represent the local degree-of-freedom values.
+ *
+ * @param[in] dof_values A vector field associated with local
+ * degree-of-freedom values on the current finite element. These define
+ * the values of all independent variables. When considering taped AD
+ * numbers with branching functions, to avoid potential issues with branch
+ * switching it may be a good idea to choose these values close or equal
+ * to those that will be later evaluated and linearized around.
+ *
+ * @note The input value type must correspond to this class's @p scalar_type.
+ * Depending on the selected @p ADNumberTypeCode, this may or may not
+ * correspond with the @p ScalarType prescribed as a template argument.
+ *
+ * @note For taped AD numbers, this operation is only valid in recording mode.
+ */
+ void
+ register_dof_values(const std::vector<scalar_type> &dof_values);
+
+ /**
+ * Register the complete set of independent variables $\mathbf{X}$ that
+ * represent the local degree-of-freedom values.
+ *
+ * @param[in] values A global field from which the values of all
+ * independent variables will be extracted. This typically will be the
+ * solution vector around which point a residual vector is to be
+ * computed and around which linearization is to occur.
+ * When considering taped AD numbers with branching functions, to avoid
+ * potential issues with branch switching it may be a good idea to choose
+ * these values close or equal to those that will be later evaluated and
+ * linearized around.
+ * @param[in] local_dof_indices A vector of degree-of-freedom indices from
+ * which to extract the local degree-of-freedom values. This would
+ * typically obtained by calling <code>cell->get_dof_indices()</code>.
+ *
+ * @note For taped AD numbers, this operation is only valid in recording mode.
+ */
+ template <typename VectorType>
+ void
+ register_dof_values(
+ const VectorType & values,
+ const std::vector<dealii::types::global_dof_index> &local_dof_indices);
+
+ /**
+ * Returns the complete set of degree-of-freedom values as represented by
+ * auto-differentiable numbers. These are the independent
+ * variables $\mathbf{X}$ about which the solution is linearized.
+ *
+ * It is indicated to the AD library that operations performed with these
+ * numbers are to be tracked, so they are considered "sensitive"
+ * variables. This is, therefore, the set of variables with which one
+ * would then perform computations, and based on which one can then
+ * extract both the value of the function and its derivatives with the
+ * member functions below. The values of the components of the returned
+ * object are initialized to the values set with
+ * register_independent_variable().
+ *
+ * @return An array of auto-differentiable type numbers representing the
+ * local degree-of-freedom values.
+ *
+ * @note For taped AD numbers, this operation is only valid in recording mode.
+ */
+ const std::vector<ad_type> &
+ get_sensitive_dof_values();
+
+ //@}
+
+ /**
+ * @name Post-processing
+ */
+ //@{
+
+ /*
+ * Returns the complete set of degree-of-freedom values of
+ * auto-differentiable number type. These store the same scalar values as
+ * the independent variables $\mathbf{X}$ about which the solution is
+ * linearized.
+ *
+ * Operations performed with these numbers are not tracked by the AD,
+ * libraries so they are considered "non-sensitive" variables.
+ * The values of the components of the returned object are initialized to
+ * the values set with register_dof_values().
+ *
+ * @return An array of auto-differentiable type numbers representing the
+ * local degree-of-freedom values.
+ *
+ * @note This function is not typically used within the context of automatic
+ * differentation computations, but can make performing substutitions in
+ * other post-processing computations more straight forward.
+ *
+ * @note For taped AD numbers, this operation is only valid outside recording mode.
+ */
+ std::vector<ad_type>
+ get_non_sensitive_dof_values() const;
+
+ //@}
+
+ /**
+ * @name Operations specific to taped mode: Reusing tapes
+ */
+ //@{
+
+ /**
+ * Set the values for the independent variables $\mathbf{X}$, i.e. the
+ * linearization point.
+ *
+ * @param[in] dof_values A vector field associated with local
+ * degree-of-freedom values on the current finite element. These define
+ * the values of all independent variables.
+ *
+ * @note The input value type must correspond to this class's @p scalar_type.
+ * Depending on the selected @p ADNumberTypeCode, this may or may not
+ * correspond with the @p ScalarType prescribed as a template argument.
+ *
+ * @note If the keep flag has been set when
+ * ADHelperBase::start_recording_operations() is called then the tape is
+ * immediately usable after creation, and the values of the independent
+ * variables set by register_dof_values() are those at which the function
+ * is to be evaluated. In this case, a separate call to this function is
+ * not strictly necessary.
+ */
+ void
+ set_dof_values(const std::vector<scalar_type> &dof_values);
+
+ /**
+ * Set the values for the independent variables $\mathbf{X}$, i.e. the
+ * linearization point.
+ *
+ * @param[in] values A vector field from which the values of all
+ * independent variables is to be extracted.
+ * @param[in] local_dof_indices A vector of degree-of-freedom indices from
+ * which to extract the local degree-of-freedom values. This would
+ * typically obtained by calling <code>cell->get_dof_indices()</code>.
+ *
+ * @note If the keep flag has been set when
+ * ADHelperBase::start_recording_operations() is called then the tape is
+ * immediately usable after creation, and the values of the independent
+ * variables set by register_dof_values() are those at which the function
+ * is to be evaluated. In this case, a separate call to this function is
+ * not strictly necessary.
+ */
+ template <typename VectorType>
+ void
+ set_dof_values(
+ const VectorType & values,
+ const std::vector<dealii::types::global_dof_index> &local_dof_indices);
+
+ //@}
+
+ /**
+ * @name Dependent variables
+ */
+ //@{
+
+ /**
+ * Computes the value of the residual vector field
+ * $\mathbf{r}(\mathbf{X})$.
+ *
+ * @param[out] residual A Vector object with the value for each component
+ * of the vector field evaluated at the point defined by the independent
+ * variable values.
+ *
+ * @note The size of the @p residual vector is determined by the derived
+ * classes, as it depends on the order of the dependent variable(s)
+ * derivative(s) that it represents. Code examples that show how to use
+ * this interface will be provided in the documentation of the derived
+ * classes.
+ */
+ virtual void
+ compute_residual(Vector<scalar_type> &residual) const = 0;
+
+ /**
+ * Computes the gradient (first derivative) of the residual vector field
+ * with respect to all independent variables, i.e.
+ * @f[
+ * \frac{\partial\mathbf{r}(\mathbf{X})}{\partial\mathbf{X}}
+ * @f]
+ *
+ * @param[out] linearization A FullMatrix with the gradient of each
+ * component of the vector field evaluated at the point defined by the
+ * independent variable values.
+ *
+ * @note The dimensions of the @p linearization matrix is determined by
+ * the derived classes, as it depends on the order of the dependent
+ * variable(s) derivative(s) that it represents. Code examples that show
+ * how to use this interface will be provided in the documentation of
+ * the derived classes.
+ */
+ virtual void
+ compute_linearization(FullMatrix<scalar_type> &linearization) const = 0;
+
+ //@}
+
+ }; // class ADHelperCellLevelBase
+
+
+ } // namespace AD
+} // namespace Differentiation
+
+
+/* ----------------- inline and template functions ----------------- */
+
+
+# ifndef DOXYGEN
+
+namespace Differentiation
+{
+ namespace AD
+ {
+ /* ----------------- ADHelperCellLevelBase ----------------- */
+
+
+
+ template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
+ template <typename VectorType>
+ void
+ ADHelperCellLevelBase<ADNumberTypeCode, ScalarType>::register_dof_values(
+ const VectorType & values,
+ const std::vector<dealii::types::global_dof_index> &local_dof_indices)
+ {
+ // This is actually the same thing the set_dof_values() function,
+ // in the sense that we simply populate our array of independent values
+ // with a meaningful number. However, in this case we need to double check
+ // that we're not registering these variables twice
+ Assert(
+ local_dof_indices.size() == this->n_independent_variables(),
+ ExcMessage(
+ "Degree-of-freedom index vector size does not match number of independent variables"));
+# ifdef DEBUG
+ for (unsigned int i = 0; i < this->n_independent_variables(); ++i)
+ {
+ Assert(this->registered_independent_variable_values[i] == false,
+ ExcMessage("Independent variables already registered."));
+ }
+# endif
+ set_dof_values(values, local_dof_indices);
+ }
+
+
+
+ template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
+ template <typename VectorType>
+ void
+ ADHelperCellLevelBase<ADNumberTypeCode, ScalarType>::set_dof_values(
+ const VectorType & values,
+ const std::vector<dealii::types::global_dof_index> &local_dof_indices)
+ {
+ Assert(local_dof_indices.size() == this->n_independent_variables(),
+ ExcMessage(
+ "Vector size does not match number of independent variables"));
+ for (unsigned int i = 0; i < this->n_independent_variables(); ++i)
+ ADHelperBase<ADNumberTypeCode, ScalarType>::set_sensitivity_value(
+ i, values[local_dof_indices[i]]);
+ }
+
+
} // namespace AD
} // namespace Differentiation
+# endif // DOXYGEN
+
+
DEAL_II_NAMESPACE_CLOSE
#endif // defined(DEAL_II_WITH_ADOLC) || defined(DEAL_II_TRILINOS_WITH_SACADO)
func);
registered_marked_dependent_variables[index] = true;
}
+
+
+
+ /* -------------------- ADHelperCellLevelBase -------------------- */
+
+
+
+ template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
+ ADHelperCellLevelBase<ADNumberTypeCode, ScalarType>::ADHelperCellLevelBase(
+ const unsigned int n_independent_variables,
+ const unsigned int n_dependent_variables)
+ : ADHelperBase<ADNumberTypeCode, ScalarType>(n_independent_variables,
+ n_dependent_variables)
+ {}
+
+
+
+ template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
+ void
+ ADHelperCellLevelBase<ADNumberTypeCode, ScalarType>::register_dof_values(
+ const std::vector<scalar_type> &dof_values)
+ {
+ // This is actually the same thing the set_independent_variable function,
+ // in the sense that we simply populate our array of independent values
+ // with a meaningful number. However, in this case we need to double check
+ // that we're not registering these variables twice
+ Assert(dof_values.size() == this->n_independent_variables(),
+ ExcMessage(
+ "Vector size does not match number of independent variables"));
+ for (unsigned int i = 0; i < this->n_independent_variables(); ++i)
+ {
+ Assert(this->registered_independent_variable_values[i] == false,
+ ExcMessage("Independent variable value already registered."));
+ }
+ set_dof_values(dof_values);
+ }
+
+
+
+ template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
+ const std::vector<
+ typename ADHelperCellLevelBase<ADNumberTypeCode, ScalarType>::ad_type> &
+ ADHelperCellLevelBase<ADNumberTypeCode,
+ ScalarType>::get_sensitive_dof_values()
+ {
+ if (ADNumberTraits<ad_type>::is_taped == true)
+ {
+ Assert(this->active_tape() != numbers::invalid_tape_index,
+ ExcMessage("Invalid tape index"));
+ }
+
+ // If necessary, initialize the internally stored vector of
+ // AD numbers that represents the independent variables
+ this->finalize_sensitive_independent_variables();
+ Assert(this->independent_variables.size() ==
+ this->n_independent_variables(),
+ ExcInternalError());
+
+ return this->independent_variables;
+ }
+
+
+
+ template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
+ std::vector<
+ typename ADHelperCellLevelBase<ADNumberTypeCode, ScalarType>::ad_type>
+ ADHelperCellLevelBase<ADNumberTypeCode,
+ ScalarType>::get_non_sensitive_dof_values() const
+ {
+ if (ADNumberTraits<ad_type>::is_taped == true)
+ {
+ Assert(this->active_tape() != numbers::invalid_tape_index,
+ ExcMessage("Invalid tape index"));
+ }
+
+ std::vector<ad_type> out(this->n_independent_variables(),
+ dealii::internal::NumberType<ad_type>::value(
+ 0.0));
+ for (unsigned int i = 0; i < this->n_independent_variables(); ++i)
+ this->initialize_non_sensitive_independent_variable(i, out[i]);
+
+ return out;
+ }
+
+
+
+ template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
+ void
+ ADHelperCellLevelBase<ADNumberTypeCode, ScalarType>::set_dof_values(
+ const std::vector<scalar_type> &values)
+ {
+ if (ADNumberTraits<ad_type>::is_taped == true)
+ {
+ Assert(this->active_tape() != numbers::invalid_tape_index,
+ ExcMessage("Invalid tape index"));
+ }
+ Assert(values.size() == this->n_independent_variables(),
+ ExcMessage(
+ "Vector size does not match number of independent variables"));
+ for (unsigned int i = 0; i < this->n_independent_variables(); ++i)
+ ADHelperBase<ADNumberTypeCode, ScalarType>::set_sensitivity_value(
+ i, values[i]);
+ }
+
+
} // namespace AD
} // namespace Differentiation
projects / dealii.git / commitdiff