Skip to content
Snippets Groups Projects
Commit bb0e72ac authored by Oliver Sander's avatar Oliver Sander Committed by sander
Browse files

Generalize to allow FD approximation using boost multiprecision number types 

Currently, the only use of FD is to double check the derivatives computed
either by hand-written code or (more frequently) by automatic differentiation.
However, only multiprecision number types allow to really be sure that the
rounding/cut off errors of FD do not hide subtle errors.

Unfortunately, this patch doesn't go all the way.  It now won't run with 'double'
anymore.  Also, user code needs to be modified in small but inelegant ways
(e.g., no 'std::' prefixes for math functions).  Since the FD tests don't
have to be run every day that's okay for now.

[[Imported from SVN: r9884]]
parent ced0bc6a
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
dune/gfe/localgeodesicfefdstiffness.hh
+ 69
− 48
View file @ bb0e72ac
......@@ -8,16 +8,15 @@
/** \brief Assembles energy gradient and Hessian with ADOL-C (automatic differentiation)
*/
template<class GridView, class LocalFiniteElement, class TargetSpace>
template<class GridView, class LocalFiniteElement, class TargetSpace, class field_type=double>
class LocalGeodesicFEFDStiffness
: public LocalGeodesicFEStiffness<GridView,LocalFiniteElement,TargetSpace>
{
// grid types
typedef typename GridView::Grid::ctype DT;
typedef typename TargetSpace::ctype RT;
typedef typename GridView::template Codim<0>::Entity Entity;
typedef typename TargetSpace::template rebind<double>::other ATargetSpace;
typedef typename TargetSpace::template rebind<field_type>::other ATargetSpace;
// some other sizes
enum {gridDim=GridView::dimension};
......@@ -35,7 +34,7 @@ public:
{}
/** \brief Compute the energy at the current configuration */
virtual RT energy (const Entity& element,
virtual double energy (const Entity& element,
const LocalFiniteElement& localFiniteElement,
const std::vector<TargetSpace>& localSolution) const
{
......@@ -73,8 +72,8 @@ public:
};
template <class GridView, class LocalFiniteElement, class TargetSpace>
void LocalGeodesicFEFDStiffness<GridView, LocalFiniteElement, TargetSpace>::
template <class GridView, class LocalFiniteElement, class TargetSpace, class field_type>
void LocalGeodesicFEFDStiffness<GridView, LocalFiniteElement, TargetSpace, field_type>::
assembleGradient(const Entity& element,
const LocalFiniteElement& localFiniteElement,
const std::vector<TargetSpace>& localSolution,
......@@ -85,35 +84,45 @@ assembleGradient(const Entity& element,
// Compute gradient by finite-difference approximation
// ///////////////////////////////////////////////////////////
double eps = 1e-6;
field_type eps = 1e-6;
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
}
localGradient.resize(localSolution.size());
std::vector<TargetSpace> forwardSolution = localSolution;
std::vector<TargetSpace> backwardSolution = localSolution;
std::vector<ATargetSpace> forwardSolution = localASolution;
std::vector<ATargetSpace> backwardSolution = localASolution;
for (size_t i=0; i<localSolution.size(); i++) {
// basis vectors of the tangent space of the i-th entry of localSolution
const Dune::FieldMatrix<RT,blocksize,embeddedBlocksize> B = localSolution[i].orthonormalFrame();
const Dune::FieldMatrix<field_type,blocksize,embeddedBlocksize> B = localSolution[i].orthonormalFrame();
for (int j=0; j<blocksize; j++) {
typename TargetSpace::EmbeddedTangentVector forwardCorrection = B[j];
typename ATargetSpace::EmbeddedTangentVector forwardCorrection = B[j];
forwardCorrection *= eps;
typename TargetSpace::EmbeddedTangentVector backwardCorrection = B[j];
typename ATargetSpace::EmbeddedTangentVector backwardCorrection = B[j];
backwardCorrection *= -eps;
forwardSolution[i] = TargetSpace::exp(localSolution[i], forwardCorrection);
backwardSolution[i] = TargetSpace::exp(localSolution[i], backwardCorrection);
forwardSolution[i] = ATargetSpace::exp(localASolution[i], forwardCorrection);
backwardSolution[i] = ATargetSpace::exp(localASolution[i], backwardCorrection);
localGradient[i][j] = (energy(element,localFiniteElement,forwardSolution) - energy(element,localFiniteElement, backwardSolution)) / (2*eps);
field_type foo = (localEnergy_->energy(element,localFiniteElement,forwardSolution) - localEnergy_->energy(element,localFiniteElement, backwardSolution)) / (2*eps);
localGradient[i][j] = foo.template convert_to<double>();
}
forwardSolution[i] = localSolution[i];
backwardSolution[i] = localSolution[i];
forwardSolution[i] = localASolution[i];
backwardSolution[i] = localASolution[i];
}
......@@ -123,8 +132,8 @@ assembleGradient(const Entity& element,
// ///////////////////////////////////////////////////////////
// Compute gradient by finite-difference approximation
// ///////////////////////////////////////////////////////////
template <class GridType, class LocalFiniteElement, class TargetSpace>
void LocalGeodesicFEFDStiffness<GridType, LocalFiniteElement, TargetSpace>::
template <class GridType, class LocalFiniteElement, class TargetSpace, class field_type>
void LocalGeodesicFEFDStiffness<GridType, LocalFiniteElement, TargetSpace, field_type>::
assembleGradientAndHessian(const Entity& element,
const LocalFiniteElement& localFiniteElement,
const std::vector<TargetSpace>& localSolution,
......@@ -138,7 +147,16 @@ assembleGradientAndHessian(const Entity& element,
this->A_ = 0;
const double eps = 5e-4;
const field_type eps = 1e-10;
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];
}
std::vector<Dune::FieldMatrix<double,blocksize,embeddedBlocksize> > B(localSolution.size());
for (size_t i=0; i<B.size(); i++)
......@@ -146,28 +164,28 @@ assembleGradientAndHessian(const Entity& element,
// Precompute negative energy at the current configuration
// (negative because that is how we need it as part of the 2nd-order fd formula)
RT centerValue = -energy(element, localFiniteElement, localSolution);
field_type centerValue = -localEnergy_->energy(element, localFiniteElement, localASolution);
// Precompute energy infinitesimal corrections in the directions of the local basis vectors
std::vector<Dune::array<RT,blocksize> > forwardEnergy(nDofs);
std::vector<Dune::array<RT,blocksize> > backwardEnergy(nDofs);
std::vector<Dune::array<field_type,blocksize> > forwardEnergy(nDofs);
std::vector<Dune::array<field_type,blocksize> > backwardEnergy(nDofs);
//#pragma omp parallel for schedule (dynamic)
for (size_t i=0; i<localSolution.size(); i++) {
for (size_t i2=0; i2<blocksize; i2++) {
typename TargetSpace::EmbeddedTangentVector epsXi = B[i][i2];
typename ATargetSpace::EmbeddedTangentVector epsXi = B[i][i2];
epsXi *= eps;
typename TargetSpace::EmbeddedTangentVector minusEpsXi = epsXi;
typename ATargetSpace::EmbeddedTangentVector minusEpsXi = epsXi;
minusEpsXi *= -1;
std::vector<TargetSpace> forwardSolution = localSolution;
std::vector<TargetSpace> backwardSolution = localSolution;
std::vector<ATargetSpace> forwardSolution = localASolution;
std::vector<ATargetSpace> backwardSolution = localASolution;
forwardSolution[i] = TargetSpace::exp(localSolution[i],epsXi);
backwardSolution[i] = TargetSpace::exp(localSolution[i],minusEpsXi);
forwardSolution[i] = ATargetSpace::exp(localASolution[i],epsXi);
backwardSolution[i] = ATargetSpace::exp(localASolution[i],minusEpsXi);
forwardEnergy[i][i2] = energy(element, localFiniteElement, forwardSolution);
backwardEnergy[i][i2] = energy(element, localFiniteElement, backwardSolution);
forwardEnergy[i][i2] = localEnergy_->energy(element, localFiniteElement, forwardSolution);
backwardEnergy[i][i2] = localEnergy_->energy(element, localFiniteElement, backwardSolution);
}
......@@ -181,8 +199,10 @@ assembleGradientAndHessian(const Entity& element,
for (size_t i=0; i<localSolution.size(); i++)
for (int j=0; j<blocksize; j++)
localGradient[i][j] = (forwardEnergy[i][j] - backwardEnergy[i][j]) / (2*eps);
{
field_type foo = (forwardEnergy[i][j] - backwardEnergy[i][j]) / (2*eps);
localGradient[i][j] = foo.template convert_to<double>();
}
///////////////////////////////////////////////////////////////////////////
// Compute Riemannian Hesse matrix by finite-difference approximation.
......@@ -195,33 +215,34 @@ assembleGradientAndHessian(const Entity& element,
for (size_t j=0; j<=i; j++) {
for (size_t j2=0; j2<((i==j) ? i2+1 : blocksize); j2++) {
std::vector<TargetSpace> forwardSolutionXiEta = localSolution;
std::vector<TargetSpace> backwardSolutionXiEta = localSolution;
std::vector<ATargetSpace> forwardSolutionXiEta = localASolution;
std::vector<ATargetSpace> backwardSolutionXiEta = localASolution;
typename TargetSpace::EmbeddedTangentVector epsXi = B[i][i2]; epsXi *= eps;
typename TargetSpace::EmbeddedTangentVector epsEta = B[j][j2]; epsEta *= eps;
typename ATargetSpace::EmbeddedTangentVector epsXi = B[i][i2]; epsXi *= eps;
typename ATargetSpace::EmbeddedTangentVector epsEta = B[j][j2]; epsEta *= eps;
typename TargetSpace::EmbeddedTangentVector minusEpsXi = epsXi; minusEpsXi *= -1;
typename TargetSpace::EmbeddedTangentVector minusEpsEta = epsEta; minusEpsEta *= -1;
typename ATargetSpace::EmbeddedTangentVector minusEpsXi = epsXi; minusEpsXi *= -1;
typename ATargetSpace::EmbeddedTangentVector minusEpsEta = epsEta; minusEpsEta *= -1;
if (i==j)
forwardSolutionXiEta[i] = TargetSpace::exp(localSolution[i],epsXi+epsEta);
forwardSolutionXiEta[i] = ATargetSpace::exp(localASolution[i],epsXi+epsEta);
else {
forwardSolutionXiEta[i] = TargetSpace::exp(localSolution[i],epsXi);
forwardSolutionXiEta[j] = TargetSpace::exp(localSolution[j],epsEta);
forwardSolutionXiEta[i] = ATargetSpace::exp(localASolution[i],epsXi);
forwardSolutionXiEta[j] = ATargetSpace::exp(localASolution[j],epsEta);
}
if (i==j)
backwardSolutionXiEta[i] = TargetSpace::exp(localSolution[i],minusEpsXi+minusEpsEta);
backwardSolutionXiEta[i] = ATargetSpace::exp(localASolution[i],minusEpsXi+minusEpsEta);
else {
backwardSolutionXiEta[i] = TargetSpace::exp(localSolution[i],minusEpsXi);
backwardSolutionXiEta[j] = TargetSpace::exp(localSolution[j],minusEpsEta);
backwardSolutionXiEta[i] = ATargetSpace::exp(localASolution[i],minusEpsXi);
backwardSolutionXiEta[j] = ATargetSpace::exp(localASolution[j],minusEpsEta);
}
RT forwardValue = energy(element, localFiniteElement, forwardSolutionXiEta) - forwardEnergy[i][i2] - forwardEnergy[j][j2];
RT backwardValue = energy(element, localFiniteElement, backwardSolutionXiEta) - backwardEnergy[i][i2] - backwardEnergy[j][j2];
field_type forwardValue = localEnergy_->energy(element, localFiniteElement, forwardSolutionXiEta) - forwardEnergy[i][i2] - forwardEnergy[j][j2];
field_type backwardValue = localEnergy_->energy(element, localFiniteElement, backwardSolutionXiEta) - backwardEnergy[i][i2] - backwardEnergy[j][j2];
this->A_[i][j][i2][j2] = this->A_[j][i][j2][i2] = 0.5 * (forwardValue - 2*centerValue + backwardValue) / (eps*eps);
field_type foo = 0.5 * (forwardValue - 2*centerValue + backwardValue) / (eps*eps);
this->A_[i][j][i2][j2] = this->A_[j][i][j2][i2] = foo.template convert_to<double>();
}
}
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment