Skip to content
Snippets Groups Projects
Commit fadc6cab authored by Sander, Oliver's avatar Sander, Oliver
Browse files

Remove the assemblers implemented in adolctest.cc itself 

There is no point in testing assemblers that only exist in the test
itself.  Historically, they may have been useful when developing
the dune-gfe assembler framework because they were simpler (they
compute the Euclidean derivatives, not the Riemannian ones).
Nowadays, they just waste CPU cycles.
parent a7ed4215
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 7 additions and 369 deletions
test/adolctest.cc
+ 7
− 369
View file @ fadc6cab
......@@ -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);
......
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