From 32b93767a8b094e8468381248cb5c6f8a7c3cfe2 Mon Sep 17 00:00:00 2001 From: Patrick Jaap <patrick.jaap@tu-dresden.de> Date: Wed, 10 Aug 2022 16:24:06 +0200 Subject: [PATCH] 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. --- dune/elasticity/assemblers/feassembler.hh | 46 +++++++++++++++++++ .../assemblers/localadolcstiffness.hh | 44 ++++++++++++++---- .../elasticity/assemblers/localfestiffness.hh | 7 +++ .../assemblers/localhyperdualstiffness.hh | 40 ++++++++++++++++ test/localhyperdualstiffnesstest.cc | 16 ++++++- 5 files changed, 143 insertions(+), 10 deletions(-) diff --git a/dune/elasticity/assemblers/feassembler.hh b/dune/elasticity/assemblers/feassembler.hh index 19a71f1..8b922c6 100644 --- a/dune/elasticity/assemblers/feassembler.hh +++ b/dune/elasticity/assemblers/feassembler.hh @@ -54,6 +54,11 @@ public: 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 * * This may be more efficient than computing them separately @@ -71,6 +76,47 @@ public: }; // 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 MatrixType> void FEAssembler<Basis,VectorType>:: diff --git a/dune/elasticity/assemblers/localadolcstiffness.hh b/dune/elasticity/assemblers/localadolcstiffness.hh index 3dd6512..ec6d51e 100644 --- a/dune/elasticity/assemblers/localadolcstiffness.hh +++ b/dune/elasticity/assemblers/localadolcstiffness.hh @@ -41,6 +41,14 @@ public: virtual double energy (const LocalView& localView, 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 * * This uses the automatic differentiation toolbox ADOL_C. @@ -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> void LocalADOLCStiffness<LocalView>:: assembleGradientAndHessian(const LocalView& localView, @@ -97,24 +130,19 @@ assembleGradientAndHessian(const LocalView& localView, ) { 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()); + // this tapes the energy computation + assembleGradient(localView, localConfiguration, localGradient); ///////////////////////////////////////////////////////////////// // Compute Hessian ///////////////////////////////////////////////////////////////// + size_t nDoubles = localConfiguration.size(); localHessian.setSize(nDoubles,nDoubles); // allocate raw doubles: ADOL-C requires "**double" arguments diff --git a/dune/elasticity/assemblers/localfestiffness.hh b/dune/elasticity/assemblers/localfestiffness.hh index bcb5ef8..d622727 100644 --- a/dune/elasticity/assemblers/localfestiffness.hh +++ b/dune/elasticity/assemblers/localfestiffness.hh @@ -16,6 +16,13 @@ class LocalFEStiffness 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 */ virtual void assembleGradientAndHessian(const LocalView& localView, diff --git a/dune/elasticity/assemblers/localhyperdualstiffness.hh b/dune/elasticity/assemblers/localhyperdualstiffness.hh index 0589359..398e18c 100644 --- a/dune/elasticity/assemblers/localhyperdualstiffness.hh +++ b/dune/elasticity/assemblers/localhyperdualstiffness.hh @@ -27,6 +27,15 @@ public: virtual double energy (const LocalView& localView, 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 * * This uses the automatic differentiation using hyperdual numbers @@ -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> void LocalHyperDualStiffness<LocalView>:: assembleGradientAndHessian(const LocalView& localView, diff --git a/test/localhyperdualstiffnesstest.cc b/test/localhyperdualstiffnesstest.cc index f0c8887..f60dcef 100644 --- a/test/localhyperdualstiffnesstest.cc +++ b/test/localhyperdualstiffnesstest.cc @@ -85,10 +85,10 @@ TestSuite hyperdualEqualsADOL_C() Timer timer; 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(); 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; @@ -97,6 +97,18 @@ TestSuite hyperdualEqualsADOL_C() // check the relative error 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; } -- GitLab