Skip to content
Snippets Groups Projects
Commit 32b93767 authored by Patrick Jaap's avatar Patrick Jaap
Browse files

Add assembleGradient() method to FE assemblers

This adds a method to compute only the gradient.

On the pro side it saves
a lot of time if only the gradient is needed (for first order methods).

On the con side it introduces some duplicated code (but only for hyberdual
numbers).

A test case is added.
parent 950f4a03
No related branches found
No related tags found
1 merge request!72Add assembleGradient() method to FE assemblers
...@@ -54,6 +54,11 @@ public: ...@@ -54,6 +54,11 @@ public:
localStiffness_(localStiffness) localStiffness_(localStiffness)
{} {}
/** \brief Assemble the functional gradient */
void assembleGradient(const VectorType& configuration,
VectorType& gradient) const;
/** \brief Assemble the tangent stiffness matrix and the functional gradient together /** \brief Assemble the tangent stiffness matrix and the functional gradient together
* *
* This may be more efficient than computing them separately * This may be more efficient than computing them separately
...@@ -71,6 +76,47 @@ public: ...@@ -71,6 +76,47 @@ public:
}; // end class }; // end class
template <class Basis, class VectorType>
void FEAssembler<Basis,VectorType>::
assembleGradient(const VectorType& configuration,
VectorType& gradient) const
{
// create backends for multi index access
auto configurationBackend = Dune::Functions::istlVectorBackend(configuration);
auto gradientBackend = Dune::Functions::istlVectorBackend(gradient);
gradientBackend.resize(basis_);
gradient = 0;
auto localView = basis_.localView();
for (const auto& element : elements(basis_.gridView(), partitionType))
{
localView.bind(element);
const auto size = localView.size();
// Extract local configuration
std::vector<field_type> localConfiguration(size);
std::vector<field_type> localGradient(size);
for (size_t i=0; i<size; i++)
localConfiguration[i] = configurationBackend[ localView.index(i) ];
// setup local matrix and gradient
localStiffness_->assembleGradient(localView, localConfiguration, localGradient);
for (size_t i=0; i<size; i++)
{
auto row = localView.index(i);
gradientBackend[row] += localGradient[i];
}
}
}
template <class Basis, class VectorType> template <class Basis, class VectorType>
template <class MatrixType> template <class MatrixType>
void FEAssembler<Basis,VectorType>:: void FEAssembler<Basis,VectorType>::
......
...@@ -41,6 +41,14 @@ public: ...@@ -41,6 +41,14 @@ public:
virtual double energy (const LocalView& localView, virtual double energy (const LocalView& localView,
const std::vector<double>& localConfiguration) const; const std::vector<double>& localConfiguration) const;
/** \brief Assemble the local gradient at the current position
*
* This uses the automatic differentiation toolbox ADOL_C.
*/
virtual void assembleGradient(const LocalView& localView,
const std::vector<double>& localConfiguration,
std::vector<double>& localGradient);
/** \brief Assemble the local stiffness matrix at the current position /** \brief Assemble the local stiffness matrix at the current position
* *
* This uses the automatic differentiation toolbox ADOL_C. * This uses the automatic differentiation toolbox ADOL_C.
...@@ -88,6 +96,31 @@ energy(const LocalView& localView, ...@@ -88,6 +96,31 @@ energy(const LocalView& localView,
} }
template<class LocalView>
void LocalADOLCStiffness<LocalView>::
assembleGradient(const LocalView& localView,
const std::vector<double>& localConfiguration,
std::vector<double>& localGradient
)
{
int rank = Dune::MPIHelper::getCommunication().rank();
// Tape energy computation. We may not have to do this every time, but it's comparatively cheap.
energy(localView, localConfiguration);
/////////////////////////////////////////////////////////////////
// Compute the energy gradient
/////////////////////////////////////////////////////////////////
// Compute the actual gradient
size_t nDoubles = localConfiguration.size();
// Compute gradient
localGradient.resize(nDoubles);
gradient(rank,nDoubles,localConfiguration.data(),localGradient.data());
}
template<class LocalView> template<class LocalView>
void LocalADOLCStiffness<LocalView>:: void LocalADOLCStiffness<LocalView>::
assembleGradientAndHessian(const LocalView& localView, assembleGradientAndHessian(const LocalView& localView,
...@@ -97,24 +130,19 @@ assembleGradientAndHessian(const LocalView& localView, ...@@ -97,24 +130,19 @@ assembleGradientAndHessian(const LocalView& localView,
) )
{ {
int rank = Dune::MPIHelper::getCommunication().rank(); int rank = Dune::MPIHelper::getCommunication().rank();
// Tape energy computation. We may not have to do this every time, but it's comparatively cheap.
energy(localView, localConfiguration);
///////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////
// Compute the energy gradient // Compute the energy gradient
///////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////
// Compute the actual gradient // this tapes the energy computation
size_t nDoubles = localConfiguration.size(); assembleGradient(localView, localConfiguration, localGradient);
// Compute gradient
localGradient.resize(nDoubles);
gradient(rank,nDoubles,localConfiguration.data(),localGradient.data());
///////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////
// Compute Hessian // Compute Hessian
///////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////
size_t nDoubles = localConfiguration.size();
localHessian.setSize(nDoubles,nDoubles); localHessian.setSize(nDoubles,nDoubles);
// allocate raw doubles: ADOL-C requires "**double" arguments // allocate raw doubles: ADOL-C requires "**double" arguments
......
...@@ -16,6 +16,13 @@ class LocalFEStiffness ...@@ -16,6 +16,13 @@ class LocalFEStiffness
public: public:
/** \brief Assemble the local gradient at the current position
*/
virtual void assembleGradient(const LocalView& localView,
const std::vector<RT>& localConfiguration,
std::vector<RT>& localGradient
) = 0;
/** \brief Assemble the local gradient and tangent matrix at the current position /** \brief Assemble the local gradient and tangent matrix at the current position
*/ */
virtual void assembleGradientAndHessian(const LocalView& localView, virtual void assembleGradientAndHessian(const LocalView& localView,
......
...@@ -27,6 +27,15 @@ public: ...@@ -27,6 +27,15 @@ public:
virtual double energy (const LocalView& localView, virtual double energy (const LocalView& localView,
const std::vector<double>& localConfiguration) const; const std::vector<double>& localConfiguration) const;
/** \brief Assemble the local gradient at the current position
*
* This uses the automatic differentiation using hyperdual numbers
*/
virtual void assembleGradient(const LocalView& localView,
const std::vector<double>& localConfiguration,
std::vector<double>& localGradient);
/** \brief Assemble the local stiffness matrix at the current position /** \brief Assemble the local stiffness matrix at the current position
* *
* This uses the automatic differentiation using hyperdual numbers * This uses the automatic differentiation using hyperdual numbers
...@@ -62,6 +71,37 @@ energy(const LocalView& localView, ...@@ -62,6 +71,37 @@ energy(const LocalView& localView,
} }
template<class LocalView>
void LocalHyperDualStiffness<LocalView>::
assembleGradient(const LocalView& localView,
const std::vector<double>& localConfiguration,
std::vector<double>& localGradient
)
{
size_t nDoubles = localConfiguration.size();
localGradient.resize(nDoubles);
// Create hyperdual vector localHyperDualConfiguration = double vector localConfiguration
std::vector<hyperdual> localHyperDualConfiguration(nDoubles);
for (size_t i=0; i<nDoubles; i++)
localHyperDualConfiguration[i] = hyperdual(localConfiguration[i]);
// Compute gradient and hessian (symmetry is used) using hyperdual function evaluation
// localHyperDualConfiguration puts Ones in the eps1 and eps2 component to evaluate the different partial derivatives
for (size_t i=0; i<nDoubles; i++)
{
localHyperDualConfiguration[i].seteps1(1.0);
hyperdual temp = localEnergy_->energy(localView, localHyperDualConfiguration);
localGradient[i] = temp.eps1();
localHyperDualConfiguration[i].seteps1(0.0);
}
}
template<class LocalView> template<class LocalView>
void LocalHyperDualStiffness<LocalView>:: void LocalHyperDualStiffness<LocalView>::
assembleGradientAndHessian(const LocalView& localView, assembleGradientAndHessian(const LocalView& localView,
......
...@@ -85,10 +85,10 @@ TestSuite hyperdualEqualsADOL_C() ...@@ -85,10 +85,10 @@ TestSuite hyperdualEqualsADOL_C()
Timer timer; Timer timer;
hyperdualAssembler.assembleGradientAndHessian(x,hyperdualGradient,hyperdualHessian); hyperdualAssembler.assembleGradientAndHessian(x,hyperdualGradient,hyperdualHessian);
std::cout << dim <<"D: hyperdual assembler took " << timer.elapsed() << " sec." << std::endl; std::cout << dim <<"D: hyperdual hessian and gradient assembler took " << timer.elapsed() << " sec." << std::endl;
timer.reset(); timer.reset();
adolcAssembler.assembleGradientAndHessian(x,adolcGradient,adolcHessian); adolcAssembler.assembleGradientAndHessian(x,adolcGradient,adolcHessian);
std::cout << dim << "D: ADOL-C assembler took " << timer.elapsed() << " sec." << std::endl; std::cout << dim << "D: ADOL-C hessian and gradient assembler took " << timer.elapsed() << " sec." << std::endl;
hyperdualHessian -= adolcHessian; hyperdualHessian -= adolcHessian;
...@@ -97,6 +97,18 @@ TestSuite hyperdualEqualsADOL_C() ...@@ -97,6 +97,18 @@ TestSuite hyperdualEqualsADOL_C()
// check the relative error // check the relative error
t.check(hyperdualHessian.frobenius_norm()/adolcHessian.frobenius_norm() < 1e-12 && hyperdualGradient.two_norm()/adolcGradient.two_norm() < 1e-12); t.check(hyperdualHessian.frobenius_norm()/adolcHessian.frobenius_norm() < 1e-12 && hyperdualGradient.two_norm()/adolcGradient.two_norm() < 1e-12);
timer.reset();
hyperdualAssembler.assembleGradient(x,hyperdualGradient);
std::cout << dim <<"D: hyperdual gradient assembler took " << timer.elapsed() << " sec." << std::endl;
timer.reset();
adolcAssembler.assembleGradient(x,adolcGradient);
std::cout << dim << "D: ADOL-C gradient assembler took " << timer.elapsed() << " sec." << std::endl;
hyperdualGradient -= adolcGradient;
// check the relative error
t.check(hyperdualGradient.two_norm()/adolcGradient.two_norm() < 1e-12);
return t; return t;
} }
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment