From 0787ae4a453ffc9191fa33aa2c5c178ccfe2daac Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Sun, 5 Oct 2014 19:10:26 +0000
Subject: [PATCH] Also test the Riemannian Hesse matrix
[[Imported from SVN: r9912]]
---
test/adolctest.cc | 122 +++++++++++++++++++++++++++++++++++-----------
1 file changed, 93 insertions(+), 29 deletions(-)
diff --git a/test/adolctest.cc b/test/adolctest.cc
index 0d710f5a..129a9931 100644
--- a/test/adolctest.cc
+++ b/test/adolctest.cc
@@ -6,8 +6,13 @@
#include <boost/multiprecision/mpfr.hpp>
-//typedef double FDType;
+//#define MULTIPRECISION
+
+#ifdef MULTIPRECISION
typedef boost::multiprecision::mpfr_float_50 FDType;
+#else
+typedef double FDType;
+#endif
// Includes for the ADOL-C automatic differentiation library
// Need to come before (almost) all others.
@@ -31,6 +36,8 @@ typedef boost::multiprecision::mpfr_float_50 FDType;
#include <dune/gfe/localgeodesicfestiffness.hh>
#include <dune/gfe/localgeodesicfefunction.hh>
#include <dune/gfe/cosseratenergystiffness.hh>
+#include <dune/gfe/localgeodesicfeadolcstiffness.hh>
+#include <dune/gfe/localgeodesicfefdstiffness.hh>
// grid dimension
const int dim = 2;
@@ -45,7 +52,7 @@ using namespace Dune;
/** \brief Assembles energy gradient and Hessian with ADOL-C
*/
template<class GridView, class LocalFiniteElement>
-class LocalGeodesicFEADOLCStiffness
+class LocalADOLCStiffness
{
// grid types
typedef typename GridView::Grid::ctype DT;
@@ -62,7 +69,7 @@ public:
//! Dimension of the embedding space
enum { embeddedBlocksize = TargetSpace::EmbeddedTangentVector::dimension };
- LocalGeodesicFEADOLCStiffness(const LocalGeodesicFEStiffness<GridView, LocalFiniteElement, ATargetSpace>* energy)
+ LocalADOLCStiffness(const LocalGeodesicFEStiffness<GridView, LocalFiniteElement, ATargetSpace>* energy)
: localEnergy_(energy)
{}
@@ -88,8 +95,8 @@ public:
template <class GridView, class LocalFiniteElement>
-typename LocalGeodesicFEADOLCStiffness<GridView, LocalFiniteElement>::RT
-LocalGeodesicFEADOLCStiffness<GridView, LocalFiniteElement>::
+typename LocalADOLCStiffness<GridView, LocalFiniteElement>::RT
+LocalADOLCStiffness<GridView, LocalFiniteElement>::
energy(const Entity& element,
const LocalFiniteElement& localFiniteElement,
const std::vector<TargetSpace>& localSolution) const
@@ -136,7 +143,7 @@ energy(const Entity& element,
// as well return it. This saves assembly time.
// ///////////////////////////////////////////////////////////
template <class GridType, class LocalFiniteElement>
-void LocalGeodesicFEADOLCStiffness<GridType, LocalFiniteElement>::
+void LocalADOLCStiffness<GridType, LocalFiniteElement>::
assembleGradientAndHessian(const Entity& element,
const LocalFiniteElement& localFiniteElement,
const std::vector<TargetSpace>& localSolution,
@@ -203,7 +210,7 @@ assembleGradientAndHessian(const Entity& element,
/** \brief Assembles energy gradient and Hessian with finite differences
*/
template<class GridView, class LocalFiniteElement, class field_type=double>
-class LocalGeodesicFEFDStiffness
+class LocalFDStiffness
{
// grid types
typedef typename GridView::Grid::ctype DT;
@@ -220,7 +227,7 @@ public:
//! Dimension of the embedding space
enum { embeddedBlocksize = TargetSpace::EmbeddedTangentVector::dimension };
- LocalGeodesicFEFDStiffness(const LocalGeodesicFEStiffness<GridView, LocalFiniteElement, ATargetSpace>* energy)
+ LocalFDStiffness(const LocalGeodesicFEStiffness<GridView, LocalFiniteElement, ATargetSpace>* energy)
: localEnergy_(energy)
{}
@@ -237,7 +244,7 @@ public:
// Compute gradient by finite-difference approximation
// ///////////////////////////////////////////////////////////
template <class GridType, class LocalFiniteElement, class field_type>
-void LocalGeodesicFEFDStiffness<GridType, LocalFiniteElement, field_type>::
+void LocalFDStiffness<GridType, LocalFiniteElement, field_type>::
assembleGradientAndHessian(const Entity& element,
const LocalFiniteElement& localFiniteElement,
const std::vector<TargetSpace>& localSolution,
@@ -251,7 +258,11 @@ assembleGradientAndHessian(const Entity& element,
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());
@@ -308,7 +319,11 @@ assembleGradientAndHessian(const Entity& element,
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
}
///////////////////////////////////////////////////////////////////////////
@@ -349,8 +364,11 @@ assembleGradientAndHessian(const Entity& element,
field_type backwardValue = localEnergy_->energy(element, localFiniteElement, 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
}
}
}
@@ -359,15 +377,16 @@ assembleGradientAndHessian(const Entity& element,
// Compare two matrices
-void compareMatrices(const Matrix<FieldMatrix<double,7,7> >& matrixA, std::string nameA,
- const Matrix<FieldMatrix<double,7,7> >& matrixB, std::string nameB)
+template <int N>
+void compareMatrices(const Matrix<FieldMatrix<double,N,N> >& matrixA, std::string nameA,
+ const Matrix<FieldMatrix<double,N,N> >& matrixB, std::string nameB)
{
double maxAbsDifference = -1;
double maxRelDifference = -1;
- for(int i=0; i<matrixA.N(); i++) {
+ for(size_t i=0; i<matrixA.N(); i++) {
- for (int j=0; j<matrixA.M(); j++ ) {
+ for (size_t j=0; j<matrixA.M(); j++ ) {
for (int ii=0; ii<matrixA[i][j].N(); ii++)
for (int jj=0; jj<matrixA[i][j].M(); jj++)
@@ -396,6 +415,7 @@ int main (int argc, char *argv[]) try
{
typedef std::vector<TargetSpace> SolutionType;
enum { embeddedBlocksize = TargetSpace::EmbeddedTangentVector::dimension };
+ enum { blocksize = TargetSpace::TangentVector::dimension };
// ///////////////////////////////////////
// Create the grid
@@ -444,28 +464,46 @@ int main (int argc, char *argv[]) try
// ////////////////////////////////////////////////////////////
ParameterTree materialParameters;
- materialParameters["thickness"] = "2.5e-5";
- materialParameters["mu"] = "5.6452e+09";
- materialParameters["lambda"] = "2.1796e+09";
- materialParameters["mu_c"] = "5.6452e+09";
+ materialParameters["thickness"] = "1";
+ materialParameters["mu"] = "1";
+ materialParameters["lambda"] = "1";
+ materialParameters["mu_c"] = "1";
materialParameters["L_c"] = "1";
materialParameters["q"] = "2";
materialParameters["kappa"] = "1";
+ ///////////////////////////////////////////////////////////////////////
+ // Assemblers for the Euclidean derivatives in an embedding space
+ ///////////////////////////////////////////////////////////////////////
+
// Assembler using ADOL-C
CosseratEnergyLocalStiffness<GridView,
FEBasis::LocalFiniteElement,
3,adouble> cosseratEnergyADOLCLocalStiffness(materialParameters, nullptr, nullptr);
- LocalGeodesicFEADOLCStiffness<GridView,
- FEBasis::LocalFiniteElement> localGFEADOLCStiffness(&cosseratEnergyADOLCLocalStiffness);
+ LocalADOLCStiffness<GridView,
+ FEBasis::LocalFiniteElement> localADOLCStiffness(&cosseratEnergyADOLCLocalStiffness);
CosseratEnergyLocalStiffness<GridView,
FEBasis::LocalFiniteElement,
3,FDType> cosseratEnergyFDLocalStiffness(materialParameters, nullptr, nullptr);
+ LocalFDStiffness<GridView,
+ FEBasis::LocalFiniteElement,FDType> localFDStiffness(&cosseratEnergyFDLocalStiffness);
+
+ ///////////////////////////////////////////////////////////////////////
+ // Assemblers for the Riemannian derivatives without embedding space
+ ///////////////////////////////////////////////////////////////////////
+
+ // Assembler using ADOL-C
+ LocalGeodesicFEADOLCStiffness<GridView,
+ FEBasis::LocalFiniteElement,
+ TargetSpace> localGFEADOLCStiffness(&cosseratEnergyADOLCLocalStiffness);
+
LocalGeodesicFEFDStiffness<GridView,
- FEBasis::LocalFiniteElement,FDType> localGFEFDStiffness(&cosseratEnergyFDLocalStiffness);
+ FEBasis::LocalFiniteElement,
+ TargetSpace,
+ FDType> localGFEFDStiffness(&cosseratEnergyFDLocalStiffness);
// Compute and compare matrices
auto it = gridView.template begin<0>();
@@ -477,7 +515,7 @@ int main (int argc, char *argv[]) try
const int numOfBaseFct = feBasis.getLocalFiniteElement(*it).localBasis().size();
- // Extract local solution
+ // Extract local configuration
std::vector<TargetSpace> localSolution(numOfBaseFct);
for (int i=0; i<numOfBaseFct; i++)
@@ -488,28 +526,54 @@ int main (int argc, char *argv[]) try
std::vector<Dune::FieldVector<double,embeddedBlocksize> > localFDGradient(numOfBaseFct);
Matrix<FieldMatrix<double,embeddedBlocksize,embeddedBlocksize> > localADHessian;
- Matrix<FieldMatrix<double,7,7> > localADVMHessian; // VM: vector-mode
- Matrix<FieldMatrix<double,7,7> > localFDHessian;
-
- // setup local matrix and gradient
- localGFEADOLCStiffness.assembleGradientAndHessian(*it,
+ Matrix<FieldMatrix<double,embeddedBlocksize,embeddedBlocksize> > localADVMHessian; // VM: vector-mode
+ Matrix<FieldMatrix<double,embeddedBlocksize,embeddedBlocksize> > localFDHessian;
+ if (false) {
+ // Assemble Euclidean derivatives
+ localADOLCStiffness.assembleGradientAndHessian(*it,
feBasis.getLocalFiniteElement(*it),
localSolution,
localADGradient,
localADHessian,
false); // 'true' means 'vector mode'
- localGFEADOLCStiffness.assembleGradientAndHessian(*it,
+ localADOLCStiffness.assembleGradientAndHessian(*it,
feBasis.getLocalFiniteElement(*it),
localSolution,
localADGradient,
localADVMHessian,
true); // 'true' means 'vector mode'
- localGFEFDStiffness.assembleGradientAndHessian(*it, feBasis.getLocalFiniteElement(*it), localSolution, localFDGradient, localFDHessian);
+ localFDStiffness.assembleGradientAndHessian(*it,
+ feBasis.getLocalFiniteElement(*it),
+ localSolution,
+ localFDGradient,
+ localFDHessian);
+ // compare
compareMatrices(localADHessian, "AD", localFDHessian, "FD");
compareMatrices(localADHessian, "AD scalar", localADVMHessian, "AD vector");
+ }
+ // Assemble Riemannian derivatives
+ std::vector<Dune::FieldVector<double,blocksize> > localRiemannianADGradient(numOfBaseFct);
+ std::vector<Dune::FieldVector<double,blocksize> > localRiemannianFDGradient(numOfBaseFct);
+
+ Matrix<FieldMatrix<double,blocksize,blocksize> > localRiemannianADHessian;
+ Matrix<FieldMatrix<double,blocksize,blocksize> > localRiemannianFDHessian;
+
+ localGFEADOLCStiffness.assembleGradientAndHessian(*it,
+ feBasis.getLocalFiniteElement(*it),
+ localSolution,
+ localRiemannianADGradient);
+
+ localGFEFDStiffness.assembleGradientAndHessian(*it,
+ feBasis.getLocalFiniteElement(*it),
+ localSolution,
+ localRiemannianFDGradient);
+
+ // compare
+ compareMatrices(localGFEADOLCStiffness.A_, "Riemannian AD", localGFEFDStiffness.A_, "Riemannian FD");
+
}
// //////////////////////////////
--
GitLab