]> https://gitweb.dealii.org/ - dealii.git/commitdiff
AD Helpers: Add helper for residual linearisation (cell-level)
authorJean-Paul Pelteret <jppelteret@gmail.com>
Wed, 25 Jul 2018 21:13:12 +0000 (23:13 +0200)
committerJean-Paul Pelteret <jppelteret@gmail.com>
Fri, 12 Oct 2018 08:39:50 +0000 (10:39 +0200)
doc/news/changes/major/20180930Jean-PaulPelteret [new file with mode: 0644]
include/deal.II/differentiation/ad/ad_helpers.h
source/differentiation/ad/ad_helpers.cc
source/differentiation/ad/ad_helpers.inst1.in
source/differentiation/ad/ad_helpers.inst2.in

diff --git a/doc/news/changes/major/20180930Jean-PaulPelteret b/doc/news/changes/major/20180930Jean-PaulPelteret
new file mode 100644 (file)
index 0000000..7b6ec47
--- /dev/null
@@ -0,0 +1,5 @@
+New: A new class Differentiation::AD::ADHelperResidualLinearisation has been
+added to help compute the linearization of a residual vector defined on the 
+level of a cell (a finite element residual), or for local nonlinear equations. 
+<br>
+(Jean-Paul Pelteret, 2018/09/30)
index c131f31fdc359d57a3e39991f95cfe4e0f32e8e7..75bd125da65575c839f445bae9e7acd6320918ed 100644 (file)
@@ -1277,6 +1277,300 @@ namespace Differentiation
     }; // class ADHelperEnergyFunctional
 
 
+
+    /**
+     * A helper class that facilitates the evaluation and automated
+     * linearization of a vector of functions that represents a residual vector
+     * (as computed from some corresponding local degree of freedom values).
+     * This class would typically be used to compute the linearization of a
+     * residual vector defined on the level of a cell, or for local
+     * nonlinear equations.
+     *
+     * An example of its usage in the case of a residual linearization
+     * might be as follows (in this case we'll compute the
+     * linearization of a finite-strain magneto-elastic solid from the residual,
+     * as constructed from the Piola-Kirchoff stress and magnetic induction
+     * assuming the magnetic scalar potential formulation):
+     *
+     * @code
+     *   // Existing data structures:
+     *   Vector<double> solution (...); // Or another vector type
+     *   std::vector<types::global_dof_index> local_dof_indices (...);
+     *   const FEValuesExtractors::Vector u_fe (...);
+     *   const FEValuesExtractors::Scalar msp_fe (...);
+     *   const unsigned int u_block (...);
+     *   const unsigned int msp_block (...);
+     *   FEValues<dim> fe_values (...);
+     *   const unsigned int n_q_points (...);
+     *   FullMatrix<double> cell_matrix (...);
+     *   Vector<double> cell_rhs (...);
+     *
+     *   // Assembly loop:
+     *   for (auto cell & : ...)
+     *   {
+     *     cell->get_dof_indices(local_dof_indices);
+     *     const unsigned int n_independent_variables
+     *       = local_dof_indices.size();
+     *     const unsigned int n_dependent_variables
+     *       = local_dof_indices.size();
+     *
+     *     // Create some aliases for the AD helper.
+     *     // In this example, we strictly assume that we're using tapeless
+     *     // AD types, and so the AD_typecode used in the template argument
+     *     // must refer to one of these types. See the example for the
+     *     // ADHelperEnergyFunctional for details on how to extend
+     *     // support to taped AD numbers.
+     *     using ADHelper = AD::ADHelperEnergyFunctional<...>;
+     *     using ADNumberType = typename ADHelper::ad_type;
+     *
+     *     // Create and initialize an instance of the helper class.
+     *     ADHelper ad_helper(n_independent_variables,n_dependent_variables);
+     *
+     *     // Initialize the local data structures for assembly.
+     *     // This is also taken care of by the ADHelper, so this step could
+     *     // be skipped.
+     *     cell_rhs.reinit(n_dependent_variables);
+     *     cell_matrix.reinit(n_independent_variables,n_dependent_variables);
+     *
+     *     // This next code block corresponds to the "recording" phase, where
+     *     // the operations performed with the AD numbers are tracked and
+     *     // differentiation is performed.
+     *     {
+     *       // First, we set the values for all DoFs.
+     *       ad_helper.register_dof_values(solution, local_dof_indices);
+     *
+     *       // Then we get the complete set of degree of freedom values as
+     *       // represented by auto-differentiable numbers. The operations
+     *       // performed with these variables are tracked by the AD library
+     *       // from this point until stop_recording_operations() is called.
+     *       const std::vector<ADNumberType> dof_values_ad
+     *         = ad_helper.get_sensitive_dof_values();
+     *
+     *       // Then we do some problem specific tasks, the first being to
+     *       // compute all values, gradients, etc. based on sensitive AD DoF
+     *       // values. Here we are fetching the displacement gradients at each
+     *       // quadrature point, as well as the gradients of the magnetic
+     *       // scalar potential field.
+     *       std::vector<Tensor<2, dim, ADNumberType>> Grad_u(
+     *         n_q_points, Tensor<2, dim, ADNumberType>());
+     *       std::vector<Tensor<1, dim, ADNumberType>> Grad_msp(
+     *         n_q_points, Tensor<1, dim, ADNumberType>());
+     *       fe_values[u_fe].get_function_gradients_from_local_dof_values(
+     *         dof_values_ad, Grad_u);
+     *       fe_values[msp_fe].get_function_gradients_from_local_dof_values(
+     *         dof_values_ad, Grad_msp);
+     *
+     *       // This variable stores the cell residual vector contributions.
+     *       // IMPORTANT: Note that each entry is hand-initialized with a value
+     *       // of zero. This is a highly recommended practise, as some AD
+     *       // numbers appear not to safely initialize their internal data
+     *       // structures.
+     *       std::vector<ADNumberType> residual_ad (
+     *         n_dependent_variables, ADNumberType(0.0));
+     *
+     *       // Compute the cell total residual
+     *       //   = (internal + external) contributions
+     *       for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
+     *       {
+     *         // Calculate the deformation gradient and magnetic field at this
+     *         // quadrature point
+     *         const Tensor<2, dim, ADNumberType> F =
+     *           unit_symmetric_tensor<dim>() + Grad_u[q_point];
+     *         const Tensor<1, dim, ADNumberType> H = -Grad_msp[q_point];
+     *         Assert(numbers::value_is_greater_than(determinant(F), 0.0),
+     *                ExcMessage("Negative determinant of the deformation "
+     *                           "gradient detected!"));
+     *
+     *         // Extract some configuration dependent variables from our
+     *         // nonlinear constitutive law for the current quadrature point.
+     *         // In this way they only have to be computed once per quadrature
+     *         // point.
+     *         const SymmetricTensor<2,dim,ad_type> S = get_S(F,H);
+     *         const Tensor<1,dim,ad_type>          B = get_B(F,H);
+     *
+     *         // Define some position-dependent aliases, to make the assembly
+     *         // process easier to follow.
+     *         const double JxW = fe_values.JxW(q_point);
+     *
+     *         // Add contribution of the internal forces:
+     *         // Note that we assemble the residual contribution directly
+     *         // as it is this vector that is to be automatically linearized.
+     *         for (unsigned int I = 0; I < n_dofs_per_cell; ++I)
+     *         {
+     *           // Determine the component and block associated with
+     *           // the I'th local degree of freedom.
+     *           const unsigned int block_I =
+     *             fe.system_to_base_index(I).first.first;
+     *
+     *           if (block_I == u_block) // u-terms
+     *           {
+     *             // Variation of the Green-Lagrange strain tensor
+     *             // associated with the I'th vector-valued basis function.
+     *             const SymmetricTensor<2,dim,double> dE_I
+     *               = symmetrize(transpose(F)
+     *               * fe_values[u_fe].gradient(I,q_point));
+     *
+     *             residual_ad[I] += (dE_I*S) * JxW;
+     *           }
+     *           else if (block_I == msp_block)
+     *           {
+     *             // Variation of the magnetic field vector associated with
+     *             // the I'th scalar-valued basis function
+     *             const Tensor<1,dim,double> dH_I
+     *               = -fe_values[msp_fe].gradient(I, q_point);
+     *
+     *             residual_ad[I] -= (dH_I*B) * JxW;
+     *           }
+     *         }
+     *       }
+     *
+     *       // Add contribution from external sources. If these contributions
+     *       // are also solution dependent then they will also be consistently
+     *       // linearized.
+     *       // Loop over faces and accumulate external contributions into the
+     *       // cell total residual.
+     *       for (unsigned int face : ...)
+     *         if (cell->face(face)->at_boundary())
+     *           residual_ad[I] += ...
+     *
+     *       // Register the definition of the cell residual
+     *       ad_helper.register_residual_vector(residual_ad);
+     *     }
+     *
+     *     // Compute the residual values and their Jacobian at the
+     *     // evaluation point
+     *     ad_helper.compute_residual(cell_rhs);
+     *     cell_rhs *= -1.0; // RHS = - residual
+     *     ad_helper.compute_linearization(cell_matrix);
+     *   }
+     * @endcode
+     *
+     * In most use cases, and in particular in the code example shown above,
+     * both the number of independent and dependent variables equals the
+     * number of <code>dofs_per_cell</code> for the used finite element.
+     *
+     * @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 ADHelperResidualLinearisation
+      : public ADHelperCellLevelBase<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{r}(\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{r}(\mathbf{X})$, this will
+       * be the number of outputs $\mathbf{r}$, i.e. the dimension of the
+       * image space.
+       */
+      ADHelperResidualLinearisation(const unsigned int n_independent_variables,
+                                    const unsigned int n_dependent_variables);
+
+      /**
+       * Destructor
+       */
+      virtual ~ADHelperResidualLinearisation() = default;
+
+      //@}
+
+      /**
+       * @name Dependent variables
+       */
+      //@{
+
+      /**
+       * Register the definition of the cell residual vector
+       * $\mathbf{r}(\mathbf{X})$.
+       *
+       * @param[in] residual A vector of recorded functions that defines the
+       * residual. The components of this vector represents the dependent
+       * variables.
+       *
+       * @note For this class that expects only vector fields of dependent
+       * variables, this function must only be called once per tape.
+       *
+       * @note For taped AD numbers, this operation is only valid in recording mode.
+       */
+      void
+      register_residual_vector(const std::vector<ad_type> &residual);
+
+      /**
+       * Evaluation of the residual for a chosen set of degree of freedom
+       * values. This corresponds to the computation of the residual vector,
+       * i.e.
+       * @f[
+       *   \mathbf{r}(\mathbf{X}) \vert_{\mathbf{X}}
+       * @f]
+       *
+       * The values at the evaluation point $\mathbf{X}$ are obtained by calling
+       * ADHelperCellLevelBase::set_dof_values().
+       *
+       * @param[out] residual A Vector object, for which the value for each
+       * entry represents the residual value for the corresponding local
+       * degree of freedom. The output @p residual vector has a length
+       * corresponding to @p n_dependent_variables.
+       */
+      virtual void
+      compute_residual(Vector<scalar_type> &residual) const override;
+
+      /**
+       * Computes the linearization of the residual vector around a chosen set
+       * of degree of freedom values. Underlying this is the computation of the
+       * gradient (first derivative) of the residual vector $\mathbf{r}$ with
+       * respect to all independent variables, i.e.
+       * @f[
+       *   \frac{\partial\mathbf{r}(\mathbf{X})}{\partial\mathbf{X}}
+       * @f]
+       *
+       * The values at the evaluation point $\mathbf{X}$ are obtained by calling
+       * ADHelperCellLevelBase::set_dof_values().
+       *
+       * @param[out] linearization A FullMatrix representing the linearization
+       * of the residual vector. The output @p linearization matrix has
+       * dimensions corresponding to
+       * <code>n_dependent_variables</code>$\times$<code>n_independent_variables</code>.
+       */
+      virtual void
+      compute_linearization(
+        FullMatrix<scalar_type> &linearization) const override;
+
+      //@}
+
+    }; // class ADHelperResidualLinearisation
+
+
   } // namespace AD
 } // namespace Differentiation
 
index ea732e626a194e8a99b0df1b1e16a56c1bc86382..844b7653f744aefea31ca5d2455d1b96bb35a4ca 100644 (file)
@@ -977,6 +977,150 @@ namespace Differentiation
     }
 
 
+    /* ------------------- ADHelperResidualLinearisation ------------------- */
+
+
+
+    template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
+    ADHelperResidualLinearisation<ADNumberTypeCode, ScalarType>::
+      ADHelperResidualLinearisation(const unsigned int n_independent_variables,
+                                    const unsigned int n_dependent_variables)
+      : ADHelperCellLevelBase<ADNumberTypeCode, ScalarType>(
+          n_independent_variables,
+          n_dependent_variables)
+    {}
+
+
+
+    template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
+    void
+    ADHelperResidualLinearisation<ADNumberTypeCode, ScalarType>::
+      register_residual_vector(const std::vector<ad_type> &residual)
+    {
+      Assert(residual.size() == this->n_dependent_variables(),
+             ExcMessage(
+               "Vector size does not match number of dependent variables"));
+      for (unsigned int i = 0; i < this->n_dependent_variables(); ++i)
+        ADHelperBase<ADNumberTypeCode, ScalarType>::register_dependent_variable(
+          i, residual[i]);
+    }
+
+
+
+    template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
+    void
+    ADHelperResidualLinearisation<ADNumberTypeCode, ScalarType>::
+      compute_residual(Vector<scalar_type> &values) const
+    {
+      if ((ADNumberTraits<ad_type>::is_taped == true &&
+           this->taped_driver.keep_independent_values() == false) ||
+          ADNumberTraits<ad_type>::is_tapeless == true)
+        {
+          Assert(
+            this->n_registered_independent_variables() ==
+              this->n_independent_variables(),
+            ExcMessage(
+              "Not all values of sensitivities have been registered or subsequently set!"));
+        }
+      Assert(this->n_registered_dependent_variables() ==
+               this->n_dependent_variables(),
+             ExcMessage("Not all dependent variables have been registered."));
+
+      // We can neglect correctly initializing the entries as
+      // we'll be overwriting them immediately.
+      if (values.size() != this->n_dependent_variables())
+        values.reinit(this->n_dependent_variables(),
+                      true /*omit_zeroing_entries*/);
+
+      if (ADNumberTraits<ad_type>::is_taped == true)
+        {
+          Assert(this->active_tape_index() !=
+                   Numbers<ad_type>::invalid_tape_index,
+                 ExcMessage("Invalid tape index"));
+          Assert(this->is_recording() == false,
+                 ExcMessage(
+                   "Cannot compute values while tape is being recorded."));
+          Assert(this->independent_variable_values.size() ==
+                   this->n_independent_variables(),
+                 ExcDimensionMismatch(this->independent_variable_values.size(),
+                                      this->n_independent_variables()));
+
+          TapedDrivers<ad_type, scalar_type>::values(
+            this->active_tape_index(),
+            this->n_dependent_variables(),
+            this->independent_variable_values,
+            values);
+        }
+      else
+        {
+          Assert(ADNumberTraits<ad_type>::is_tapeless == true,
+                 ExcInternalError());
+          TapelessDrivers<ad_type, scalar_type>::values(
+            this->dependent_variables, values);
+        }
+    }
+
+
+
+    template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
+    void
+    ADHelperResidualLinearisation<ADNumberTypeCode, ScalarType>::
+      compute_linearization(FullMatrix<scalar_type> &jacobian) const
+    {
+      if ((ADNumberTraits<ad_type>::is_taped == true &&
+           this->taped_driver.keep_independent_values() == false) ||
+          ADNumberTraits<ad_type>::is_tapeless == true)
+        {
+          Assert(
+            this->n_registered_independent_variables() ==
+              this->n_independent_variables(),
+            ExcMessage(
+              "Not all values of sensitivities have been registered or subsequently set!"));
+        }
+      Assert(this->n_registered_dependent_variables() ==
+               this->n_dependent_variables(),
+             ExcMessage("Not all dependent variables have been registered."));
+
+      // We can neglect correctly initializing the entries as
+      // we'll be overwriting them immediately.
+      if (jacobian.m() != this->n_dependent_variables() ||
+          jacobian.n() != this->n_independent_variables())
+        jacobian.reinit({this->n_dependent_variables(),
+                         this->n_independent_variables()},
+                        true /*omit_default_initialization*/);
+
+      if (ADNumberTraits<ad_type>::is_taped == true)
+        {
+          Assert(this->active_tape_index() !=
+                   Numbers<ad_type>::invalid_tape_index,
+                 ExcMessage("Invalid tape index"));
+          Assert(this->is_recording() == false,
+                 ExcMessage(
+                   "Cannot compute hessian while tape is being recorded."));
+          Assert(this->independent_variable_values.size() ==
+                   this->n_independent_variables(),
+                 ExcDimensionMismatch(this->independent_variable_values.size(),
+                                      this->n_independent_variables()));
+
+          TapedDrivers<ad_type, scalar_type>::jacobian(
+            this->active_tape_index(),
+            this->n_dependent_variables(),
+            this->independent_variable_values,
+            jacobian);
+        }
+      else
+        {
+          Assert(ADNumberTraits<ad_type>::is_tapeless == true,
+                 ExcInternalError());
+          Assert(this->independent_variables.size() ==
+                   this->n_independent_variables(),
+                 ExcInternalError());
+          TapelessDrivers<ad_type, scalar_type>::jacobian(
+            this->independent_variables, this->dependent_variables, jacobian);
+        }
+    }
+
+
   } // namespace AD
 } // namespace Differentiation
 
index 418cdf3b78c8d50547392f60b1b908d308775ff1..79bc953f1cf27377c8816c65e8e7f963cb85dad4 100644 (file)
@@ -45,6 +45,14 @@ for (number : REAL_SCALARS)
     template
     class ADHelperEnergyFunctional<NumberTypes::adolc_tapeless,number>;
 
+    // -------------------------- ADHelperResidualLinearisation ----------------------
+
+    template
+    class ADHelperResidualLinearisation<NumberTypes::adolc_taped,number>;
+
+    template
+    class ADHelperResidualLinearisation<NumberTypes::adolc_tapeless,number>;
+
     \}
     \}
 }
@@ -81,6 +89,14 @@ for ()
     template
     class ADHelperEnergyFunctional<NumberTypes::adolc_tapeless,typename NumberTraits<double,NumberTypes::adolc_tapeless>::ad_type>;
 
+    // -------------------------- ADHelperResidualLinearisation ----------------------
+
+    template
+    class ADHelperResidualLinearisation<NumberTypes::adolc_taped,typename NumberTraits<double,NumberTypes::adolc_taped>::ad_type>;
+
+    template
+    class ADHelperResidualLinearisation<NumberTypes::adolc_tapeless,typename NumberTraits<double,NumberTypes::adolc_tapeless>::ad_type>;
+
 
     \}
     \}
index be1bc98baaf88c4a787229ee168b3b8e8a81c61f..1e64944e57246939dc8e566cddaa9b144a45a2ea 100644 (file)
@@ -63,6 +63,20 @@ for (number : REAL_SCALARS)
     template
     class ADHelperEnergyFunctional<NumberTypes::sacado_rad_dfad,number>;
 
+    // -------------------------- ADHelperResidualLinearisation ----------------------
+
+    template
+    class ADHelperResidualLinearisation<NumberTypes::sacado_dfad_dfad,number>;
+
+    template
+    class ADHelperResidualLinearisation<NumberTypes::sacado_dfad,number>;
+
+    template
+    class ADHelperResidualLinearisation<NumberTypes::sacado_rad,number>;
+
+    template
+    class ADHelperResidualLinearisation<NumberTypes::sacado_rad_dfad,number>;
+
     \}
     \}
 }
@@ -117,6 +131,20 @@ for ()
     template
     class ADHelperEnergyFunctional<NumberTypes::sacado_rad_dfad,typename NumberTraits<double,NumberTypes::sacado_rad_dfad>::ad_type>;
 
+    // -------------------------- ADHelperResidualLinearisation ----------------------
+
+    template
+    class ADHelperResidualLinearisation<NumberTypes::sacado_dfad_dfad,typename NumberTraits<double,NumberTypes::sacado_dfad_dfad>::ad_type>;
+
+    template
+    class ADHelperResidualLinearisation<NumberTypes::sacado_dfad,typename NumberTraits<double,NumberTypes::sacado_dfad>::ad_type>;
+
+    template
+    class ADHelperResidualLinearisation<NumberTypes::sacado_rad,typename NumberTraits<double,NumberTypes::sacado_rad>::ad_type>;
+
+    template
+    class ADHelperResidualLinearisation<NumberTypes::sacado_rad_dfad,typename NumberTraits<double,NumberTypes::sacado_rad_dfad>::ad_type>;
+
 
     \}
     \}

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.