Commit 923098cb authored by Sander, Oliver's avatar Sander, Oliver
Browse files

Merge branch 'spectral-speedup-patches' into 'master'

Patches to speed up TNNMG for materials with a spectral split.

See merge request !48
parents 30bc76c2 fa151a36
Pipeline #8685 passed with stage
in 4 minutes and 4 seconds
......@@ -25,8 +25,7 @@ class BrittleFractureLocalSolver
auto negativeGradient = derivative(quadraticFunctional)(x);
negativeGradient *= -1;
auto hessian = quadraticFunctional.quadraticPart();
negativeGradient -= nonlinearFunctional.gradient(x, 0);
nonlinearFunctional.addHessian(x, hessian, 0, 0);
nonlinearFunctional.addNegativeGradientHessian(x, negativeGradient, hessian, 0, 0);
// extract upper left block
FieldMatrix<double, dim, dim> H00(0);
......
......@@ -130,9 +130,6 @@ public:
DUNE_THROW(RangeError, "Not implemented for dim = " << dim);
}
const double degradationValue = f_->degradation(jet);
FieldVector<RangeField, size> stresses(0.0);
// psi,{eps eps}
if (row==0 && col==0)
{
......@@ -143,12 +140,13 @@ public:
for (std::size_t i=0; i<size-1; ++i)
result[i][i] += 2.0*f_->mu_*II[i];
result *= degradationValue;
result *= f_->degradation(jet);
}
// psi,{eps d}
if (row==0 && col==1)
{
FieldVector<RangeField, size> stresses(0.0);
const double trEps = f_->traceEps(jet);
for (auto&& diagonal : diagonals())
......@@ -163,6 +161,7 @@ public:
if (row==1 && col==0)
{
FieldVector<RangeField, size> stresses(0.0);
const double trEps = f_->traceEps(jet);
for (auto&& diagonal : diagonals())
......@@ -338,6 +337,159 @@ public:
return energy;
}
/** \brief Compute gradient and Hesse matrix together
*
* \return std::tuple of FieldVector for the gradient and FieldMatrix for the Hesse matrix
*/
auto gradientHessian(const FieldVector<RangeField, size>& state) const
{
std::tuple<typename DamagedElasticEnergyDensityGradient::Range, typename DamagedElasticEnergyDensityHessian::Range> result;
auto& gradient = std::get<0>(result);
auto& hessian = std::get<1>(result);
gradient = 0.0;
hessian = 0.0;
const double degradationValue = degradation(state);
const double trEps = traceEps(state);
FieldVector<double,size-1> II;
switch (dim)
{
case 1:
II = {1.0};
break;
case 2:
II = {1.0, 2.0, 1.0};
break;
case 3:
II = {1.0, 2.0, 1.0, 2.0, 2.0, 1.0};
break;
default:
DUNE_THROW(RangeError, "Not implemented for dim = " << dim);
}
// Compute gradient
for (auto&& diagonal : diagonals())
gradient[diagonal] = degradationValue * lambda_ * trEps;
for (std::size_t i=0; i<size-1; i++)
gradient[i] += degradationValue * 2.0 * mu_ * state[i] * II[i];
gradient[size-1] = -2.0 * (1.0 - state[size-1]) * psi_0(state);
// linear part of critical energy density for model with elastic phase
if (crackSurfaceDensity_.type() == CrackSurfaceDensity::AT1) // psi_c = g_c / (4*c_w*l)
gradient[size-1] += crackSurfaceDensity_.g_c() / (4.0*crackSurfaceDensity_.c_w()*crackSurfaceDensity_.l());
// Compute Hesse matrix
// Displacement-displacement coupling
for (auto&& d1 : diagonals())
for (auto&& d2: diagonals())
hessian[d1][d2] += lambda_;
for (std::size_t i=0; i<size-1; ++i)
hessian[i][i] += 2.0 * mu_ * II[i];
hessian *= degradationValue;
// Displacement-damage coupling
FieldVector<RangeField, size> stresses(0.0);
for (auto&& diagonal : diagonals())
stresses[diagonal] = lambda_*trEps;
for (std::size_t i=0; i<size-1; i++)
stresses[i] += 2.0 * mu_*state[i];
for (std::size_t i=0; i<size-1; ++i)
hessian[i][size-1] = hessian[size-1][i] = -2.0*(1 - state[size-1])*stresses[i] * II[i];
// Damage-damage coupling
hessian[size-1][size-1] = 2.0 * psi_0(state);
return result;
}
/** \brief Compute a block of gradient and Hesse matrix together
*
* \return std::tuple of FieldVector for the gradient and FieldMatrix for the Hesse matrix
*/
auto gradientHessian(const FieldVector<RangeField, size>& state, size_t row, size_t col) const
{
if (row!=col) {
std::cerr << "row != col not implemented!" << std::endl;
std::terminate();
}
// Gradient part
std::tuple<typename DamagedElasticEnergyDensityGradient::Range, typename DamagedElasticEnergyDensityHessian::Range> result;
auto& gradient = std::get<0>(result);
auto& hessian = std::get<1>(result);
gradient = 0.0;
hessian = 0.0;
const double degradationValue = degradation(state);
const double trEps = traceEps(state);
// Displacement part
if (row==0)
{
// Gradient part
for (auto&& diagonal : diagonals())
gradient[diagonal] = degradationValue * lambda_ * trEps;
FieldVector<double,size-1> II;
switch (dim)
{
case 1:
II = {1.0};
break;
case 2:
II = {1.0, 2.0, 1.0};
break;
case 3:
II = {1.0, 2.0, 1.0, 2.0, 2.0, 1.0};
break;
default:
DUNE_THROW(RangeError, "Not implemented for dim = " << dim);
}
for (std::size_t i=0; i<size-1; i++)
gradient[i] += degradationValue * 2.0 * mu_ * state[i] * II[i];
// Hesse matrix part
for (auto&& d1 : diagonals())
for (auto&& d2: diagonals())
hessian[d1][d2] += lambda_;
for (std::size_t i=0; i<size-1; ++i)
hessian[i][i] += 2.0 * mu_ * II[i];
hessian *= degradationValue;
}
// Damage part
if (row==1)
{
// Gradient part
gradient[size-1] = -2.0 * (1.0 - state[size-1]) * psi_0(state);
// linear part of critical energy density for model with elastic phase
if (crackSurfaceDensity_.type() == CrackSurfaceDensity::AT1) // psi_c = g_c / (4*c_w*l)
gradient[size-1] += crackSurfaceDensity_.g_c() / (4.0*crackSurfaceDensity_.c_w()*crackSurfaceDensity_.l());
// Hesse matrix part
hessian[size-1][size-1] = 2.0 * psi_0(state);
}
return result;
}
private:
RangeField traceEps(const FieldVector<RangeField, size>& stateVector) const
{
......
......@@ -204,11 +204,24 @@ public:
return (*pimpl_)(jet);
}
auto gradientHessian(const Jet& jet) const
{
return pimpl_->gradientHessian(jet);
}
auto gradientHessian(const Jet& jet, size_t row, size_t col) const
{
return pimpl_->gradientHessian(jet, row, col);
}
private:
struct Concept
{
virtual ~Concept() {}
virtual RangeField operator() (const Jet& jet) const = 0;
virtual std::tuple<typename IntegralDensityGradient<Jet>::Range, typename IntegralDensityHessian<Jet>::Range> gradientHessian(const Jet& jet) const = 0;
virtual std::tuple<typename IntegralDensityGradient<Jet>::Range, typename IntegralDensityHessian<Jet>::Range> gradientHessian(const Jet& jet, size_t row, size_t col) const = 0;
virtual IntegralDensityGradient<Jet> derivative_() const = 0;
virtual Concept* clone() const = 0;
};
......@@ -218,6 +231,8 @@ private:
{
Model( T const& value ) : object( value ) {}
RangeField operator() (const Jet& jet) const override {return object(jet); }
std::tuple<typename IntegralDensityGradient<Jet>::Range, typename IntegralDensityHessian<Jet>::Range> gradientHessian(const Jet& jet) const override {return object.gradientHessian(jet);}
std::tuple<typename IntegralDensityGradient<Jet>::Range, typename IntegralDensityHessian<Jet>::Range> gradientHessian(const Jet& jet, size_t row, size_t col) const override {return object.gradientHessian(jet,row,col);}
IntegralDensityGradient<Jet> derivative_() const override {return IntegralDensityGradient<Jet>(derivative(object)); }
Concept* clone() const override{ return new Model( object ); }
T object;
......
......@@ -204,6 +204,57 @@ public:
}
}
/** \brief Add only a sub-block of the Hessian
*
* Block decomposition of the Hessian depends on the integral density
*
* \param row Block row of the gradient vector and Hesse matrix block to be added
* \param col Block column of the Hesse matrix block to be added
*/
template<class Hessian>
void addNegativeGradientHessian(const Vector& v, Vector& negativeGradient, Hessian& hessian,
size_t row, size_t col) const
{
const auto& BT = f_->originalFunctional().transposedEvaluationMatrix();
const auto& weights = quadratureWeights();
const auto& densityGradient = derivative(density());
const auto& densityHessian = derivative(derivative(density()));
for (auto&& [BT_ij, j] : sparseRange(BT[index()]))
{
auto jet_j = f_->jet()[j];
BT_ij.umtv(v, jet_j);
// Get gradient and Hesse matrix of the density
auto [g_j, Hj] = density().gradientHessian(jet_j, row, col);
BT_ij.usmv(-1.0 * weights[j], g_j, negativeGradient);
Hj *= weights[j];
// compute BT_ij * Hj * (BT_ij)^T
if constexpr (Std::is_detected_v<Impl::isUnblockedDiagonalMatrix_t, decltype(BT_ij)>)
BT_ij.leftRightMultiply(Hj, hessian, row, col);
else
{
auto rowRange = densityHessian.getBlockRange(row);
auto colRange = densityHessian.getBlockRange(col);
for (size_t k : rowRange)
{
for (size_t l : colRange)
{
auto&& BT_ij_k = BT_ij[k];
auto&& BT_ij_l = BT_ij[l];
for (auto&& [Hj_r, r] : sparseRange(Hj))
for (auto&& [Hj_rs, s] : sparseRange(Hj_r))
hessian[k][l] += BT_ij_k[r] * Hj_rs * BT_ij_l[s];
}
}
}
}
}
protected:
const GlobalFunctional* f_;
......
......@@ -72,23 +72,25 @@ namespace Dune::FracturePhasefields {
*/
void bind(const Vector& x)
{
// Compute first derivate of the functional
negativeGradient_ = f_.gradient(x);
negativeGradient_ *= -1;
typename LocalEnergy::Domain BiXv;
// Compute second derivative of the functional
// Compute first and second derivate of the functional
Solvers::resizeInitializeZero(negativeGradient_, x);
hessian_ = 0;
typename LocalEnergy::Domain BiXv;
for (std::size_t i=0; i<f_.evaluationMatrix().N(); ++i)
{
BiXv = 0.0;
BiXv = 0.0;
addRowXVector(f_.evaluationMatrix()[i], x, BiXv);
auto densityHessian = derivative(derivative(f_.density()));
auto DDgamma = densityHessian(BiXv);
auto [Dgamma, DDgamma] = f_.density().gradientHessian(BiXv);
Dgamma *= -1.0 * f_.quadratureWeights()[i];
DDgamma *= f_.quadratureWeights()[i];
addRowTransposedXVector(f_.evaluationMatrix()[i],
Dgamma,
negativeGradient_);
addRowTransposedXMatrixXRow(f_.evaluationMatrix()[i],
DDgamma,
hessian_);
......@@ -156,6 +158,14 @@ namespace Dune::FracturePhasefields {
entry.umv(x[i], y);
}
static void addRowTransposedXVector(const typename EvaluationMatrix::row_type& evaluationMatrixRow,
const typename LocalEnergy::Domain& x, // Actually: the range of the gradient
Vector& y)
{
for (auto [entry, i] : sparseRange(evaluationMatrixRow))
entry.umtv(x, y[i]);
}
static void addRowTransposedXMatrixXRow(const typename EvaluationMatrix::row_type& row,
const FunctionalHessian& X,
Matrix& Y)
......
......@@ -76,6 +76,17 @@ public:
addHessian(x, hessian, 1, 1);
}
void addNegativeGradientHessian(const FieldVector<double,size>& x,
FieldVector<double,size>& negativeGradient,
FieldMatrix<double,size,size>& hessian,
std::size_t row, std::size_t col) const
{
auto grad = gradient(x, row);
negativeGradient -= grad;
addHessian(x, hessian, row, col);
}
FieldVector<double,size> center_;
double p_;
};
......
......@@ -99,6 +99,17 @@ public:
return jet.two_norm2();
}
auto gradientHessian(const Jet& jet) const
{
std::tuple<FieldVector<double,dim>, ScaledIdentityMatrix<double,dim> > result;
std::get<0>(result) = derivative(*this)(jet);
std::get<1>(result) = derivative(derivative(*this))(jet);
return result;
}
friend auto derivative(const DirichletEnergyDensity& f) {
return DirichletEnergyDensityGradient(f);
}
......
Supports Markdown
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