diff --git a/test/adolctest.cc b/test/adolctest.cc index 0dd4ec71436d8578f4d1b070c6735886b53fea92..22944f5efa2496a5b228ab76e95051d599eaba02 100644 --- a/test/adolctest.cc +++ b/test/adolctest.cc @@ -54,329 +54,6 @@ const int dim = 2; // Image space of the geodesic fe functions using TargetSpace = GFE::ProductManifold<RealTuple<double,3>,Rotation<double,3> >; -/** \brief Assembles energy gradient and Hessian with ADOL-C - */ -template<class Basis> -class LocalADOLCStiffness -{ - // grid types - typedef typename Basis::GridView GridView; - typedef typename GridView::ctype DT; - typedef typename TargetSpace::ctype RT; - typedef typename GridView::template Codim<0>::Entity Entity; - - typedef typename TargetSpace::template rebind<adouble>::other ATargetSpace; - - // some other sizes - constexpr static int gridDim = GridView::dimension; - -public: - - //! Dimension of the embedding space - constexpr static int embeddedBlocksize = TargetSpace::EmbeddedTangentVector::dimension; - - LocalADOLCStiffness(const GFE::LocalEnergy<Basis, ATargetSpace>* energy) - : localEnergy_(energy) - {} - - /** \brief Compute the energy at the current configuration */ - virtual RT energy (const typename Basis::LocalView& localView, - const std::vector<TargetSpace>& localSolution) const; - - /** \brief Assemble the local stiffness matrix at the current position - - This uses the automatic differentiation toolbox ADOL_C. - */ - virtual void assembleGradientAndHessian(const typename Basis::LocalView& localView, - const std::vector<TargetSpace>& localSolution, - std::vector<Dune::FieldVector<double, embeddedBlocksize> >& localGradient, - Dune::Matrix<Dune::FieldMatrix<RT,embeddedBlocksize,embeddedBlocksize> >& localHessian, - bool vectorMode); - - const GFE::LocalEnergy<Basis, ATargetSpace>* localEnergy_; - -}; - - -template <class Basis> -typename LocalADOLCStiffness<Basis>::RT -LocalADOLCStiffness<Basis>:: -energy(const typename Basis::LocalView& localView, - const std::vector<TargetSpace>& localSolution) const -{ - double pureEnergy; - - std::vector<ATargetSpace> localASolution(localSolution.size()); - - trace_on(1); - - adouble energy = 0; - - // The following loop is not quite intuitive: we copy the localSolution into an - // array of FieldVector<double>, go from there to FieldVector<adouble> and - // only then to ATargetSpace. - // Rationale: The constructor/assignment-from-vector of TargetSpace frequently - // contains a projection onto the manifold from the surrounding Euclidean space. - // ADOL-C needs a function on the whole Euclidean space, hence that projection - // is part of the function and needs to be taped. - - // The following variable cannot be declared inside of the loop, or ADOL-C will report wrong results - // (Presumably because several independent variables use the same memory location.) - std::vector<typename ATargetSpace::CoordinateType> aRaw(localSolution.size()); - for (size_t i=0; i<localSolution.size(); i++) { - typename TargetSpace::CoordinateType raw = localSolution[i].globalCoordinates(); - for (size_t j=0; j<raw.size(); j++) - aRaw[i][j] <<= raw[j]; - localASolution[i] = aRaw[i]; // may contain a projection onto M -- needs to be done in adouble - } - - energy = localEnergy_->energy(localView,localASolution); - - energy >>= pureEnergy; - - trace_off(); - return pureEnergy; -} - - - -// /////////////////////////////////////////////////////////// -// Compute gradient and Hessian together -// To compute the Hessian we need to compute the gradient anyway, so we may -// as well return it. This saves assembly time. -// /////////////////////////////////////////////////////////// -template <class Basis> -void LocalADOLCStiffness<Basis>:: -assembleGradientAndHessian(const typename Basis::LocalView& localView, - const std::vector<TargetSpace>& localSolution, - std::vector<Dune::FieldVector<double,embeddedBlocksize> >& localGradient, - Dune::Matrix<Dune::FieldMatrix<RT,embeddedBlocksize,embeddedBlocksize> >& localHessian, - bool vectorMode) -{ - // Tape energy computation. We may not have to do this every time, but it's comparatively cheap. - energy(localView, localSolution); - - ///////////////////////////////////////////////////////////////// - // Compute the gradient. - ///////////////////////////////////////////////////////////////// - - // Copy data from Dune data structures to plain-C ones - size_t nDofs = localSolution.size(); - size_t nDoubles = nDofs*embeddedBlocksize; - std::vector<double> xp(nDoubles); - int idx=0; - for (size_t i=0; i<nDofs; i++) - for (size_t j=0; j<embeddedBlocksize; j++) - xp[idx++] = localSolution[i].globalCoordinates()[j]; - - // Compute gradient - std::vector<double> g(nDoubles); - gradient(1,nDoubles,xp.data(),g.data()); // gradient evaluation - - // Copy into Dune type - std::vector<typename TargetSpace::EmbeddedTangentVector> localEmbeddedGradient(localSolution.size()); - - idx=0; - for (size_t i=0; i<nDofs; i++) - for (size_t j=0; j<embeddedBlocksize; j++) - localGradient[i][j] = g[idx++]; - - ///////////////////////////////////////////////////////////////// - // Compute Hessian - ///////////////////////////////////////////////////////////////// - - localHessian.setSize(nDofs,nDofs); - - double* rawHessian[nDoubles]; - for(size_t i=0; i<nDoubles; i++) - rawHessian[i] = (double*)malloc((i+1)*sizeof(double)); - - if (vectorMode) - hessian2(1,nDoubles,xp.data(),rawHessian); - else - hessian(1,nDoubles,xp.data(),rawHessian); - - // Copy Hessian into Dune data type - for(size_t i=0; i<nDoubles; i++) - for (size_t j=0; j<nDoubles; j++) - { - double value = (i>=j) ? rawHessian[i][j] : rawHessian[j][i]; - localHessian[j/embeddedBlocksize][i/embeddedBlocksize][j%embeddedBlocksize][i%embeddedBlocksize] = value; - } - - for(size_t i=0; i<nDoubles; i++) - free(rawHessian[i]); - -} - -/** \brief Assembles energy gradient and Hessian with finite differences - */ -template<class Basis, class field_type=double> -class LocalFDStiffness -{ - // grid types - typedef typename Basis::GridView GridView; - typedef typename GridView::Grid::ctype DT; - typedef typename GridView::template Codim<0>::Entity Entity; - - typedef typename TargetSpace::template rebind<field_type>::other ATargetSpace; - - -public: - - //! Dimension of a tangent space - constexpr static int blocksize = TargetSpace::TangentVector::dimension; - - //! Dimension of the embedding space - constexpr static int embeddedBlocksize = TargetSpace::EmbeddedTangentVector::dimension; - - LocalFDStiffness(const GFE::LocalEnergy<Basis, ATargetSpace>* energy) - : localEnergy_(energy) - {} - - virtual void assembleGradientAndHessian(const typename Basis::LocalView& localView, - const std::vector<TargetSpace>& localSolution, - std::vector<Dune::FieldVector<double,embeddedBlocksize> >& localGradient, - Dune::Matrix<Dune::FieldMatrix<double,embeddedBlocksize,embeddedBlocksize> >& localHessian); - - const GFE::LocalEnergy<Basis, ATargetSpace>* localEnergy_; -}; - -// /////////////////////////////////////////////////////////// -// Compute gradient by finite-difference approximation -// /////////////////////////////////////////////////////////// -template <class Basis, class field_type> -void LocalFDStiffness<Basis, field_type>:: -assembleGradientAndHessian(const typename Basis::LocalView& localView, - const std::vector<TargetSpace>& localSolution, - std::vector<Dune::FieldVector<double, embeddedBlocksize> >& localGradient, - Dune::Matrix<Dune::FieldMatrix<double,embeddedBlocksize,embeddedBlocksize> >& localHessian) -{ - // Number of degrees of freedom for this element - size_t nDofs = localSolution.size(); - - // Clear assemble data - localHessian.setSize(nDofs, nDofs); - localHessian = 0; - -#ifdef MULTIPRECISION - const field_type eps = 1e-10; -#else - const field_type eps = 1e-4; -#endif - - std::vector<ATargetSpace> localASolution(localSolution.size()); - std::vector<typename ATargetSpace::CoordinateType> aRaw(localSolution.size()); - for (size_t i=0; i<localSolution.size(); i++) { - typename TargetSpace::CoordinateType raw = localSolution[i].globalCoordinates(); - for (size_t j=0; j<raw.size(); j++) - aRaw[i][j] = raw[j]; - localASolution[i] = aRaw[i]; // may contain a projection onto M -- needs to be done in adouble - } - - std::vector<Dune::FieldMatrix<field_type,embeddedBlocksize,embeddedBlocksize> > B(localSolution.size()); - for (size_t i=0; i<B.size(); i++) - { - B[i] = 0; - for (int j=0; j<embeddedBlocksize; j++) - B[i][j][j] = 1.0; - } - - // Precompute negative energy at the current configuration - // (negative because that is how we need it as part of the 2nd-order fd formula) - field_type centerValue = -localEnergy_->energy(localView, localASolution); - - // Precompute energy infinitesimal corrections in the directions of the local basis vectors - std::vector<std::array<field_type,embeddedBlocksize> > forwardEnergy(nDofs); - std::vector<std::array<field_type,embeddedBlocksize> > backwardEnergy(nDofs); - - for (size_t i=0; i<localSolution.size(); i++) { - for (size_t i2=0; i2<embeddedBlocksize; i2++) { - typename ATargetSpace::EmbeddedTangentVector epsXi = B[i][i2]; - epsXi *= eps; - typename ATargetSpace::EmbeddedTangentVector minusEpsXi = epsXi; - minusEpsXi *= -1; - - std::vector<ATargetSpace> forwardSolution = localASolution; - std::vector<ATargetSpace> backwardSolution = localASolution; - - forwardSolution[i] = ATargetSpace(localASolution[i].globalCoordinates() + epsXi); - backwardSolution[i] = ATargetSpace(localASolution[i].globalCoordinates() + minusEpsXi); - - forwardEnergy[i][i2] = localEnergy_->energy(localView, forwardSolution); - backwardEnergy[i][i2] = localEnergy_->energy(localView, backwardSolution); - - } - - } - - ////////////////////////////////////////////////////////////// - // Compute gradient by finite-difference approximation - ////////////////////////////////////////////////////////////// - - localGradient.resize(localSolution.size()); - - for (size_t i=0; i<localSolution.size(); i++) - for (int j=0; j<embeddedBlocksize; j++) - { - field_type foo = (forwardEnergy[i][j] - backwardEnergy[i][j]) / (2*eps); -#ifdef MULTIPRECISION - localGradient[i][j] = foo.template convert_to<double>(); -#else - localGradient[i][j] = foo; -#endif - } - - /////////////////////////////////////////////////////////////////////////// - // Compute Riemannian Hesse matrix by finite-difference approximation. - // We loop over the lower left triangular half of the matrix. - // The other half follows from symmetry. - /////////////////////////////////////////////////////////////////////////// - //#pragma omp parallel for schedule (dynamic) - for (size_t i=0; i<localSolution.size(); i++) { - for (size_t i2=0; i2<embeddedBlocksize; i2++) { - for (size_t j=0; j<=i; j++) { - for (size_t j2=0; j2<((i==j) ? i2+1 : embeddedBlocksize); j2++) { - - std::vector<ATargetSpace> forwardSolutionXiEta = localASolution; - std::vector<ATargetSpace> backwardSolutionXiEta = localASolution; - - typename ATargetSpace::EmbeddedTangentVector epsXi = B[i][i2]; epsXi *= eps; - typename ATargetSpace::EmbeddedTangentVector epsEta = B[j][j2]; epsEta *= eps; - - typename ATargetSpace::EmbeddedTangentVector minusEpsXi = epsXi; minusEpsXi *= -1; - typename ATargetSpace::EmbeddedTangentVector minusEpsEta = epsEta; minusEpsEta *= -1; - - if (i==j) - forwardSolutionXiEta[i] = ATargetSpace(localASolution[i].globalCoordinates() + epsXi+epsEta); - else { - forwardSolutionXiEta[i] = ATargetSpace(localASolution[i].globalCoordinates() + epsXi); - forwardSolutionXiEta[j] = ATargetSpace(localASolution[j].globalCoordinates() + epsEta); - } - - if (i==j) - backwardSolutionXiEta[i] = ATargetSpace(localASolution[i].globalCoordinates() + minusEpsXi+minusEpsEta); - else { - backwardSolutionXiEta[i] = ATargetSpace(localASolution[i].globalCoordinates() + minusEpsXi); - backwardSolutionXiEta[j] = ATargetSpace(localASolution[j].globalCoordinates() + minusEpsEta); - } - - field_type forwardValue = localEnergy_->energy(localView, forwardSolutionXiEta) - forwardEnergy[i][i2] - forwardEnergy[j][j2]; - field_type backwardValue = localEnergy_->energy(localView, backwardSolutionXiEta) - backwardEnergy[i][i2] - backwardEnergy[j][j2]; - - field_type foo = 0.5 * (forwardValue - 2*centerValue + backwardValue) / (eps*eps); -#ifdef MULTIPRECISION - localHessian[i][j][i2][j2] = localHessian[j][i][j2][i2] = foo.template convert_to<double>(); -#else - localHessian[i][j][i2][j2] = localHessian[j][i][j2][i2] = foo; -#endif - } - } - } - } -} - - // Compare two matrices template <int N> void compareMatrices(const Matrix<FieldMatrix<double,N,N> >& matrixA, std::string nameA, @@ -417,7 +94,6 @@ int main (int argc, char *argv[]) try MPIHelper::instance(argc, argv); typedef std::vector<TargetSpace> SolutionType; - constexpr static int embeddedBlocksize = TargetSpace::EmbeddedTangentVector::dimension; constexpr static int blocksize = TargetSpace::TangentVector::dimension; // /////////////////////////////////////// @@ -498,29 +174,21 @@ int main (int argc, char *argv[]) try materialParameters["b2"] = "1"; materialParameters["b3"] = "1"; - /////////////////////////////////////////////////////////////////////// - // Assemblers for the Euclidean derivatives in an embedding space - /////////////////////////////////////////////////////////////////////// - - // Assembler using ADOL-C - CosseratEnergyLocalStiffness<FEBasis, - 3,adouble> cosseratEnergyADOLCLocalStiffness(materialParameters, nullptr, nullptr, nullptr); - - LocalADOLCStiffness<FEBasis> localADOLCStiffness(&cosseratEnergyADOLCLocalStiffness); - - CosseratEnergyLocalStiffness<FEBasis, - 3,FDType> cosseratEnergyFDLocalStiffness(materialParameters, nullptr, nullptr, nullptr); - - LocalFDStiffness<FEBasis,FDType> localFDStiffness(&cosseratEnergyFDLocalStiffness); /////////////////////////////////////////////////////////////////////// - // Assemblers for the Riemannian derivatives without embedding space + // Assemblers for the Riemannian derivatives /////////////////////////////////////////////////////////////////////// // Assembler using ADOL-C + CosseratEnergyLocalStiffness<FEBasis, 3,adouble> + cosseratEnergyADOLCLocalStiffness(materialParameters, nullptr, nullptr, nullptr); + LocalGeodesicFEADOLCStiffness<FEBasis, TargetSpace> localGFEADOLCStiffness(&cosseratEnergyADOLCLocalStiffness); + // Assembler using finite differences + CosseratEnergyLocalStiffness<FEBasis, 3,FDType> + cosseratEnergyFDLocalStiffness(materialParameters, nullptr, nullptr, nullptr); LocalGeodesicFEFDStiffness<FEBasis, TargetSpace, FDType> localGFEFDStiffness(&cosseratEnergyFDLocalStiffness); @@ -541,36 +209,6 @@ int main (int argc, char *argv[]) try for (int i=0; i<numOfBaseFct; i++) localSolution[i] = x[localView.index(i)]; - std::vector<Dune::FieldVector<double,embeddedBlocksize> > localADGradient(numOfBaseFct); - std::vector<Dune::FieldVector<double,embeddedBlocksize> > localADVMGradient(numOfBaseFct); // VM: vector-mode - std::vector<Dune::FieldVector<double,embeddedBlocksize> > localFDGradient(numOfBaseFct); - - Matrix<FieldMatrix<double,embeddedBlocksize,embeddedBlocksize> > localADHessian; - Matrix<FieldMatrix<double,embeddedBlocksize,embeddedBlocksize> > localADVMHessian; // VM: vector-mode - Matrix<FieldMatrix<double,embeddedBlocksize,embeddedBlocksize> > localFDHessian; - - // Assemble Euclidean derivatives - localADOLCStiffness.assembleGradientAndHessian(localView, - localSolution, - localADGradient, - localADHessian, - false); // 'true' means 'vector mode' - - localADOLCStiffness.assembleGradientAndHessian(localView, - localSolution, - localADGradient, - localADVMHessian, - true); // 'true' means 'vector mode' - - localFDStiffness.assembleGradientAndHessian(localView, - localSolution, - localFDGradient, - localFDHessian); - - // compare - compareMatrices(localADHessian, "AD", localFDHessian, "FD"); - compareMatrices(localADHessian, "AD scalar", localADVMHessian, "AD vector"); - // Assemble Riemannian derivatives std::vector<double> localRiemannianADGradient(numOfBaseFct*blocksize); std::vector<double> localRiemannianFDGradient(numOfBaseFct*blocksize);