Commit bb2ab533 authored by lisa_julia.nebel_at_tu-dresden.de's avatar lisa_julia.nebel_at_tu-dresden.de
Browse files

Parallelize the gradient calculation for hyperdual numbers with openmp but continue to use ADOLC

parent ff213c8a
......@@ -33,7 +33,9 @@ public:
LocalADOLCStiffness(std::shared_ptr<Dune::Elasticity::LocalEnergy<LocalView, adouble>> energy, bool useHessian2 = true)
: localEnergy_(energy)
, useHessian2_(useHessian2)
{}
{
std::cout << "Using ADOL-C for differentiation." << std::endl;
}
virtual ~LocalADOLCStiffness() {}
......
#pragma once
#include <omp.h>
#include <dune/common/fmatrix.hh>
#include <dune/istl/matrix.hh>
......@@ -21,7 +22,9 @@ public:
// accepts hyperdual only
LocalHyperDualStiffness(std::shared_ptr<Dune::Elasticity::LocalEnergy<LocalView, hyperdual>> energy)
: localEnergy_(energy)
{}
{
std::cout << "Using Hyperdual differentiation with " << omp_get_num_procs() << " threads." << std::endl;
}
/** \brief Compute the energy at the current configuration */
virtual double energy (const LocalView& localView,
......@@ -51,7 +54,7 @@ energy(const LocalView& localView,
std::vector<hyperdual> localHyperDualConfiguration(localConfiguration.size());
hyperdual energy;
for (size_t i=0; i<localConfiguration.size(); i++)
for (int i=0; i<localConfiguration.size(); i++)
localHyperDualConfiguration[i] = hyperdual(localConfiguration[i]);
energy = localEnergy_->energy(localView,localHyperDualConfiguration);
......@@ -77,12 +80,14 @@ assembleGradientAndHessian(const LocalView& localView,
localHessian.setSize(nDoubles,nDoubles);
// Create hyperdual vector localHyperDualConfiguration = double vector localConfiguration
std::vector<hyperdual> localHyperDualConfiguration(nDoubles);
for (size_t i=0; i<nDoubles; i++)
#pragma omp parallel for
for (int i=0; i<nDoubles; i++)
localHyperDualConfiguration[i] = hyperdual(localConfiguration[i]);
// Compute gradient and hessian (symmetry is used) using hyperdual function evaluation
// localHyperDualConfiguration puts Ones in the eps1 and eps2 component to evaluate the different partial derivatives
for (size_t i=0; i<nDoubles; i++)
#pragma omp parallel for
for (int i=0; i<nDoubles; i++)
{
localHyperDualConfiguration[i].seteps1(1.0);
localHyperDualConfiguration[i].seteps2(1.0);
......@@ -93,7 +98,8 @@ assembleGradientAndHessian(const LocalView& localView,
localHyperDualConfiguration[i].seteps2(0.0);
for (size_t j=i+1; j<nDoubles; j++)
#pragma omp parallel for
for (int j=i+1; j<nDoubles; j++)
{
localHyperDualConfiguration[j].seteps2(1.0);
localHessian[i][j] = localEnergy_->energy(localView, localHyperDualConfiguration).eps1eps2();
......
......@@ -616,6 +616,9 @@ namespace std
template<>
struct numeric_limits<hyperdual> : std::numeric_limits<double>
{};
bool isnormal(const hyperdual a) {
return std::isnormal(a.real());
}
}
......
#define ADOLC 1
#include <config.h>
#include <functional>
......@@ -13,6 +14,7 @@
#include <dune/common/parametertree.hh>
#include <dune/common/parametertreeparser.hh>
#include <dune/elasticity/common/hyperdual.hh>
#include <dune/grid/uggrid.hh>
#include <dune/grid/utility/structuredgridfactory.hh>
......@@ -33,6 +35,7 @@
#include <dune/elasticity/common/trustregionsolver.hh>
#include <dune/elasticity/assemblers/localadolcstiffness.hh>
#include <dune/elasticity/assemblers/localhyperdualstiffness.hh>
#include <dune/elasticity/assemblers/feassembler.hh>
#include <dune/elasticity/materials/localintegralenergy.hh>
#include <dune/elasticity/materials/stvenantkirchhoffdensity.hh>
......@@ -261,35 +264,43 @@ int main (int argc, char *argv[]) try
materialParameters.report();
}
// Assembler using ADOL-C
// Assembler using ADOL-C or Hyperdual
#if ADOLC
typedef adouble ValueType;
#else
typedef hyperdual ValueType;
#endif
if (mpiHelper.rank()==0)
std::cout << "Selected energy is: " << parameterSet.get<std::string>("energy") << std::endl;
std::cout << "Selected energy is: " << parameterSet.get<std::string>("energy") << std::endl;
std::shared_ptr<Elasticity::LocalDensity<dim,adouble>> elasticDensity;
std::shared_ptr<Elasticity::LocalDensity<dim,ValueType>> elasticDensity;
if (parameterSet.get<std::string>("energy") == "stvenantkirchhoff")
elasticDensity = std::make_shared<Elasticity::StVenantKirchhoffDensity<dim,adouble>>(materialParameters);
elasticDensity = std::make_shared<Elasticity::StVenantKirchhoffDensity<dim,ValueType>>(materialParameters);
if (parameterSet.get<std::string>("energy") == "neohooke")
elasticDensity = std::make_shared<Elasticity::NeoHookeDensity<dim,adouble>>(materialParameters);
elasticDensity = std::make_shared<Elasticity::NeoHookeDensity<dim,ValueType>>(materialParameters);
if (parameterSet.get<std::string>("energy") == "hencky")
elasticDensity = std::make_shared<Elasticity::HenckyDensity<dim,adouble>>(materialParameters);
elasticDensity = std::make_shared<Elasticity::HenckyDensity<dim,ValueType>>(materialParameters);
if (parameterSet.get<std::string>("energy") == "exphencky")
elasticDensity = std::make_shared<Elasticity::ExpHenckyDensity<dim,adouble>>(materialParameters);
elasticDensity = std::make_shared<Elasticity::ExpHenckyDensity<dim,ValueType>>(materialParameters);
if (parameterSet.get<std::string>("energy") == "mooneyrivlin")
elasticDensity = std::make_shared<Elasticity::MooneyRivlinDensity<dim,adouble>>(materialParameters);
elasticDensity = std::make_shared<Elasticity::MooneyRivlinDensity<dim,ValueType>>(materialParameters);
if(!elasticDensity)
DUNE_THROW(Exception, "Error: Selected energy not available!");
auto elasticEnergy = std::make_shared<Elasticity::LocalIntegralEnergy<LocalView, adouble>>(elasticDensity);
auto neumannEnergy = std::make_shared<Elasticity::NeumannEnergy<LocalView,adouble>>(neumannBoundary,neumannFunction);
auto elasticEnergy = std::make_shared<Elasticity::LocalIntegralEnergy<LocalView, ValueType>>(elasticDensity);
auto totalEnergy = std::make_shared<Elasticity::SumEnergy<LocalView,adouble>>(elasticEnergy, neumannEnergy);
auto neumannEnergy = std::make_shared<Elasticity::NeumannEnergy<LocalView,ValueType>>(neumannBoundary,neumannFunction);
auto totalEnergy = std::make_shared<Elasticity::SumEnergy<LocalView,ValueType>>(elasticEnergy, neumannEnergy);
#if ADOLC
auto localADOLCStiffness = std::make_shared<Elasticity::LocalADOLCStiffness<LocalView>>(totalEnergy);
Elasticity::FEAssembler<PowerBasis,SolutionType> assembler(basis, localADOLCStiffness);
#else
auto localHyperDualStiffness = std::make_shared<Elasticity::LocalHyperDualStiffness<LocalView>>(totalEnergy);
Elasticity::FEAssembler<PowerBasis,SolutionType> assembler(basis, localHyperDualStiffness);
#endif
// /////////////////////////////////////////////////
// Create a Riemannian trust-region solver
......
Markdown is supported
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