Skip to content
Snippets Groups Projects
Commit f07d6fd8 authored by oliver.sander_at_tu-dresden.de's avatar oliver.sander_at_tu-dresden.de
Browse files

Merge branch 'feature/adolc-alternative-modes' into 'master'

LocalADOLCStiffness: use "hessian2" instead of "hess_vec"

See merge request agnumpde/dune-elasticity!39
parents 27528420 18229443
Branches
No related tags found
No related merge requests found
......@@ -27,8 +27,12 @@ class LocalADOLCStiffness
public:
// accepts ADOL_C's "adouble" only
LocalADOLCStiffness(std::shared_ptr<Dune::Elasticity::LocalEnergy<LocalView, adouble>> energy)
// useHessian2: there are different basic drivers for computing the hessian by ADOL-C, namely "hessian" and "hessian2".
// They yield the same results, but "hessian2" seems to be faster. The user can switch to "hessian" with
// useHessian2=false. See ADOL-C spec for more info about the drivers
LocalADOLCStiffness(std::shared_ptr<Dune::Elasticity::LocalEnergy<LocalView, adouble>> energy, bool useHessian2 = true)
: localEnergy_(energy)
, useHessian2_(useHessian2)
{}
/** \brief Compute the energy at the current configuration */
......@@ -45,6 +49,7 @@ public:
Dune::Matrix<double>& localHessian);
std::shared_ptr<Dune::Elasticity::LocalEnergy<LocalView, adouble>> localEnergy_;
const bool useHessian2_;
};
......@@ -110,29 +115,41 @@ assembleGradientAndHessian(const LocalView& localView,
localHessian.setSize(nDoubles,nDoubles);
std::vector<double> v(nDoubles);
std::vector<double> w(nDoubles);
std::fill(v.begin(), v.end(), 0.0);
for (size_t i=0; i<nDoubles; i++)
// allocate raw doubles: ADOL-C requires "**double" arguments
double* rawHessian[nDoubles];
for(size_t i=0; i<nDoubles; i++)
{
// Evaluate Hessian multiplied with the i-th canonical basis vector (i.e.: the i-th column of the Hessian)
v[i] = 1.;
// we need only the lower half of the hessian, thus allocate i+1 doubles only
rawHessian[i] = (double*)malloc((i+1)*sizeof(double));
}
// ADOL-C error code
int rc;
int rc= 3;
MINDEC(rc, hess_vec(rank, nDoubles, const_cast<double*>(localConfiguration.data()), v.data(), w.data())); // ADOL-C won't accept pointers with const values
if (rc < 0)
DUNE_THROW(Dune::Exception, "ADOL-C has returned with error code " << rc << "!");
// ADOL-C won't accept pointers with const values, therefore the const_cast
if ( useHessian2_ )
rc = hessian2(rank,nDoubles,const_cast<double*>(localConfiguration.data()),rawHessian);
else
rc = hessian(rank,nDoubles,const_cast<double*>(localConfiguration.data()),rawHessian);
for (size_t j=0; j<nDoubles; j++)
localHessian[i][j] = w[j];
if (rc < 0)
DUNE_THROW(Dune::Exception, "ADOL-C has returned with error code " << rc << "!");
// Make v the null vector again
v[i] = 0.;
// copy the raw values to the localHessian
for (size_t i=0; i<nDoubles; i++)
{
for (size_t j=0; j<i; j++)
{
// we can only access the lower half of rawHessian!
localHessian[i][j] = rawHessian[i][j];
localHessian[j][i] = rawHessian[i][j];
}
localHessian[i][i] = rawHessian[i][i];
}
// TODO: implement ADOLC_VECTOR_MODE variant
// don't forget the clean everything up
for(size_t i=0; i<nDoubles; i++)
free(rawHessian[i]);
}
} // namespace Dune::Elasticity
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment