Commit ed739722 authored by Sander, Oliver's avatar Sander, Oliver
Browse files

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

Make TNNMG faster

See merge request !47
parents 9f5be6e5 28e01732
Pipeline #8486 passed with stage
in 7 minutes and 13 seconds
......@@ -13,6 +13,8 @@ install(FILES
damageassembler.hh
elementquantityassembler.hh
globalquadrature.hh
gradienttostrainvaluemap.hh
gradientvaluejetevaluator.hh
hfielddamageassembler.hh
integraldensity.hh
integralfunctional.hh
......@@ -22,4 +24,5 @@ install(FILES
damagedelasticenergydensity_nondecomposed.hh
ductilefracturejet.hh
smallstrainbrittlefracturejetevaluator.hh
unblockeddiagonalmatrix.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/fracture-phasefields)
......@@ -21,10 +21,19 @@ public:
DuctileElasticEnergyConstrainedLinearization(const F& f, const BitVector& ignore)
: f_(f),
ignore_(ignore),
ignore_(&ignore),
truncationTolerance_(1e-10)
{}
/** \brief Specify what degrees of freedom to ignore
*/
void setIgnore(const BitVector& ignore)
{
ignore_ = &ignore;
}
/** \brief Evaluate the linearization at a given configuration
*/
void bind(const Vector& x)
{
const static int nDamage=1;
......@@ -114,7 +123,7 @@ public:
////////////////////////////////////////////////////
// The ignore_ field contains the Dirichlet dofs
truncationFlags_ = ignore_;
truncationFlags_ = *ignore_;
// determine which components to truncate for damage
for (std::size_t i=0; i<x.size(); i++)
......@@ -259,7 +268,7 @@ private:
const F& f_;
const BitVector& ignore_;
const BitVector* ignore_;
double truncationTolerance_;
......
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_FRACTUREPHASEFIELDS_GRADIENTTOSTRAINVALUEMAP_HH
#define DUNE_FRACTUREPHASEFIELDS_GRADIENTTOSTRAINVALUEMAP_HH
#include <dune/fracture-phasefields/unblockeddiagonalmatrix.hh>
namespace Dune::FracturePhasefields
{
/** \brief Take a matrix with a gradient and a value und turn it into a matrix with strains and a value
*
* This is a policy class for the `MappedMatrix` class from the `dune-fufem` module.
* It takes a matrix with entries of type FieldMatrix<dim+1,1> and replaces each entry
* with an object of type FieldMatrix<dim*(dim+1)/2+1,dim+1>.
*
* Call the entries of the input matrix are (v_0, ..., v_(dim-1), a)^T. The output matrix
* is a block-diagonal matrix with block sizes (dim*(dim+1)+1) x dim and 1.
*
* The i-th column the upper left matrix for is the symmetric strain corresponding to the vector v,
* when considered as pointing into the i-th canonical direction. The symmetric strain
* is a symmetric matrix of size dim x dim. We store only the lower triangular part,
* which gives a vector of size dim*(dim+1)/2.
*
* The lower-right matrix block is simply the value 'a' from the input matrix.
*
* \param GradientValueMatrix The input matrix whose entries are to be transformed
* \param dim The number of space dimensions
* \param transposed If this is true, then the returned matrix is not what is
* described above, but its transposed.
*/
template <typename GradientValueMatrix, int dim, bool transposed>
class GradientToStrainValueMap
{
static_assert(GradientValueMatrix::block_type::rows == ((transposed) ? 1 : dim+1),
"GradientValueMatrix entries have the wrong number of rows!");
static_assert(GradientValueMatrix::block_type::cols == ((transposed) ? dim+1 : 1),
"GradientValueMatrix entries have the wrong number of columns!");
using field_type = typename GradientValueMatrix::field_type;
/** \brief Return view onto symmetric matrix that has operator(i,j)
*
* Expected memory layout is v00, v10, v11, v20, v21, v22, v30, ...
*/
static std::size_t linearIndex(unsigned int i, unsigned int j)
{
return (i >= j) ? (i*(i+1)/2+j) : (j*(j+1)/2+i);
}
public:
static constexpr int rowFactor = 1; // number of rows will be this times larger
static constexpr int colFactor = 1; // number of columns will be this times larger
using Matrix = GradientValueMatrix; // type of matrix to be transformed
// Type of entries for transformed matrix
// A strain is a vector of four entries, and there are three shape functions
using StrainMatrix = std::conditional_t<transposed,
FieldMatrix<field_type, dim, dim*(dim+1)/2>,
FieldMatrix<field_type, dim*(dim+1)/2, dim> >;
using block_type = UnblockedDiagonalMatrix<StrainMatrix,FieldMatrix<field_type,1,1> >;
/** \brief Transform (virtualRow,virtualCol)-th copy of (row,col)-th entry A of original matrix
*/
block_type apply(const typename Matrix::block_type& A, int row, int col, int virtualRow, int virtualCol) const
{
StrainMatrix strain;
// Compute strain in i-th direction
for (unsigned int i=0; i<dim; i++)
{
for (unsigned int j=0; j<dim; j++)
if (transposed)
strain[i][linearIndex(i,j)] = ((i==j) ? 1.0 : 0.5) * A[0][j];
else
strain[linearIndex(i,j)][i] = ((i==j) ? 1.0 : 0.5) * A[j][0];
}
FieldMatrix<field_type,1,1> value = (transposed) ? A[0][A.cols-1] : A[A.rows-1][0];
return block_type(strain, value);
}
};
} // namespace Dune::FracturePhasefields
#endif // DUNE_FRACTUREPHASEFIELDS_GRADIENTTOSTRAINVALUEMAP_HH
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_FRACTUREPHASEFIELDS_GRADIENTVALUEJETEVALUATOR_HH
#define DUNE_FRACTUREPHASEFIELDS_GRADIENTVALUEJETEVALUATOR_HH
#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
namespace Dune::FracturePhasefields
{
/** \brief Evaluate gradient and value of a scalar basis at a given point
*
* \tparam Basis A scalar-valued finite element basis
*/
template <class Basis>
class GradientValueJetEvaluator
{
static constexpr std::size_t dim = Basis::GridView::dimension;
using Element = typename Basis::GridView::template Codim<0>::Entity;
public:
/** \brief Data structure holding the gradient and the value of a scalar function
*/
using Jet = FieldVector<double, dim+1>;
/** \brief Default constructor
*/
GradientValueJetEvaluator()
{}
/** \brief Evaluate the jet at a given point for all basis functions
*
* \todo Properly return 'jet' (as method return value)
*/
void operator()(const typename Basis::LocalView& localView,
std::vector<Jet>& jet,
const typename Element::Geometry::LocalCoordinate& position) const
{
const auto& tree = localView.tree();
const auto& element = localView.element();
const auto& localFiniteElement = tree.finiteElement();
using FiniteElement = decltype(localFiniteElement);
using Range = typename std::decay_t<FiniteElement>::Traits::LocalBasisType::Traits::RangeType;
using Jacobian = typename std::decay_t<FiniteElement>::Traits::LocalBasisType::Traits::JacobianType;
const auto& geometry = element.geometry();
const auto invGeometryJacobian = geometry.jacobianInverseTransposed(position);
std::vector<Range> values(localFiniteElement.localBasis().size());
localFiniteElement.localBasis().evaluateFunction(position, values);
auto referenceJacobians = std::vector<Jacobian>(localFiniteElement.size());
localFiniteElement.localBasis().evaluateJacobian(position, referenceJacobians);
for (std::size_t i=0; i<localFiniteElement.size(); ++i)
{
Jacobian jacobian = referenceJacobians[i] * transpose(invGeometryJacobian);
std::size_t index = tree.localIndex(i);
// Evaluate deformation gradient
for (std::size_t j=0; j<dim; j++)
jet[index][j] = jacobian[0][j];
jet[index][dim] = values[i];
}
}
};
} // namespace Dune::FracturePhasefields
#endif
......@@ -8,11 +8,19 @@
#include <dune/common/exceptions.hh>
#include <dune/common/typeutilities.hh>
#include <dune/solvers/common/interval.hh>
#include <dune/solvers/common/resize.hh>
namespace Dune {
namespace FracturePhasefields {
namespace Impl {
// Helper to test whether T is of type UnblockedDiagonal
template<typename T>
using isUnblockedDiagonalMatrix_t = decltype( &T::isUnblockedDiagonalMatrix );
}
/** \brief Restriction of an integral functional to a low-dimensional axis-aligned subspace
*
* In other words, if the functional is block-separable we restrict to a single block.
......@@ -113,8 +121,6 @@ public:
BT_ij.umtv(v, jet_j);
auto g_j = densityGradient(jet_j, block);
// Multiply with the entire gradient, even though we know that some entries are zero.
// I think exploiting the zeros here would complicate the code too much.
BT_ij.usmv(weights[j], g_j, gradient);
}
return gradient;
......@@ -137,11 +143,16 @@ public:
Hj *= weights[j];
// compute BT_ij * Hj * (BT_ij)^T
for (auto&& [BT_ij_k, k] : sparseRange(BT_ij))
for (auto&& [BT_ij_l, l] : sparseRange(BT_ij))
for (auto&& [Hj_r, r] : sparseRange(Hj))
for (auto&& [Hj_rs, s] : sparseRange(Hj_r))
result[k][l] += BT_ij_k[r] * Hj_rs * BT_ij_l[s];
if constexpr (Std::is_detected_v<Impl::isUnblockedDiagonalMatrix_t, decltype(BT_ij)>)
BT_ij.leftRightMultiply(Hj, result);
else
{
for (auto&& [BT_ij_k, k] : sparseRange(BT_ij))
for (auto&& [BT_ij_l, l] : sparseRange(BT_ij))
for (auto&& [Hj_r, r] : sparseRange(Hj))
for (auto&& [Hj_rs, s] : sparseRange(Hj_r))
result[k][l] += BT_ij_k[r] * Hj_rs * BT_ij_l[s];
}
}
}
......@@ -166,22 +177,28 @@ public:
auto jet_j = f_->jet()[j];
BT_ij.umtv(v, jet_j);
auto rowRange = densityHessian.getBlockRange(row);
auto colRange = densityHessian.getBlockRange(col);
auto Hj = densityHessian(jet_j, row, col);
Hj *= weights[j];
// compute BT_ij * Hj * (BT_ij)^T
for (size_t k : rowRange)
if constexpr (Std::is_detected_v<Impl::isUnblockedDiagonalMatrix_t, decltype(BT_ij)>)
BT_ij.leftRightMultiply(Hj, result, row, col);
else
{
for (size_t l : colRange)
{
auto&& BT_ij_k = BT_ij[k];
auto&& BT_ij_l = BT_ij[l];
auto rowRange = densityHessian.getBlockRange(row);
auto colRange = densityHessian.getBlockRange(col);
for (auto&& [Hj_r, r] : sparseRange(Hj))
for (auto&& [Hj_rs, s] : sparseRange(Hj_r))
result[k][l] += BT_ij_k[r] * Hj_rs * BT_ij_l[s];
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))
result[k][l] += BT_ij_k[r] * Hj_rs * BT_ij_l[s];
}
}
}
}
......@@ -310,8 +327,13 @@ public:
const GlobalVector& direction)
: f_(f),
origin_(&origin),
direction_(&direction)
{}
direction_(&direction),
evaluationMatrixTimesDirection_(f_->evaluationMatrix().N())
{
// Premultiply evaluationMatrix and direction vector.
// This is used later in the subDifferential method.
f_->evaluationMatrix().mv(*direction_, evaluationMatrixTimesDirection_);
}
Range operator()(const Scalar& v) const
{
......@@ -351,6 +373,8 @@ public:
typename Functional::Density::Domain BiXv;
double directionalDerivative = 0;
for (std::size_t i=0; i<f_->evaluationMatrix().N(); ++i)
{
BiXv = 0.0;
......@@ -360,12 +384,10 @@ public:
auto Dgamma = densityDerivative(BiXv);
Dgamma *= f_->quadratureWeights()[i];
Functional::addRowTransposedXVector(f_->evaluationMatrix()[i], Dgamma, gradient);
directionalDerivative += Dgamma * evaluationMatrixTimesDirection_[i];
}
Solvers::Interval<double> dJ = gradient * (*direction_);
return dJ;
return directionalDerivative;
}
/** \brief (Temporary?) alias for the subDifferential method
......@@ -382,6 +404,11 @@ protected:
const Functional* f_;
const GlobalVector* origin_;
const GlobalVector* direction_;
// The evaluation matrix multiplied by the direction vector
// This product is used in the subDifferential method,
// and it pays to precompute it.
BlockVector<typename Functional::Density::Domain> evaluationMatrixTimesDirection_;
};
......
......@@ -5,6 +5,8 @@
#include <dune/common/rangeutilities.hh>
#include <dune/fracture-phasefields/integralfunctional.hh>
namespace Dune::FracturePhasefields {
/**
......@@ -30,9 +32,8 @@ namespace Dune::FracturePhasefields {
*
* The constructor assembles the pattern of the Hesse matrix
*/
IntegralFunctionalConstrainedLinearization(const F& f, const BitVector& ignore)
: f_(f),
ignore_(ignore)
IntegralFunctionalConstrainedLinearization(const F& f)
: f_(f)
{
// Set the pattern for the Hesse matrix
auto nBlocks = f_.evaluationMatrix().M();
......@@ -50,6 +51,23 @@ namespace Dune::FracturePhasefields {
hessian_ = 0;
}
/** \brief Construct the linearization of a given IntegralFunctional
*
* The constructor assembles the pattern of the Hesse matrix
*/
IntegralFunctionalConstrainedLinearization(const F& f, const BitVector& ignore)
: IntegralFunctionalConstrainedLinearization(f)
{
ignore_ = &ignore;
}
/** \brief Specify what degrees of freedom to ignore
*/
void setIgnore(const BitVector& ignore)
{
ignore_ = &ignore;
}
/** \brief Evaluate the linearization at a given configuration
*/
void bind(const Vector& x)
......@@ -73,7 +91,6 @@ namespace Dune::FracturePhasefields {
addRowTransposedXMatrixXRow(f_.evaluationMatrix()[i],
DDgamma,
f_.evaluationMatrix()[i],
hessian_);
}
......@@ -83,7 +100,7 @@ namespace Dune::FracturePhasefields {
// The ignore_ field contains the Dirichlet dofs
// TODO: Should we truncate anything else?
truncationFlags_ = ignore_;
truncationFlags_ = *ignore_;
// Do the actual truncation
for (std::size_t i=0; i<x.size(); ++i)
......@@ -139,38 +156,52 @@ namespace Dune::FracturePhasefields {
entry.umv(x[i], y);
}
static void addRowTransposedXMatrixXRow(const typename EvaluationMatrix::row_type& rowA,
static void addRowTransposedXMatrixXRow(const typename EvaluationMatrix::row_type& row,
const FunctionalHessian& X,
const typename EvaluationMatrix::row_type& rowB,
Matrix& Y)
{
static const int M = EvaluationMatrix::block_type::rows;
for (auto [rowA_i, i] : sparseRange(rowA))
for (auto [rowA_i, i] : sparseRange(row))
{
for (auto [rowB_j, j] : sparseRange(rowB))
// Precompute a matrix-matrix product
auto precomp = X * rowA_i;
// Iterate only over the upper right triangle of Y.
for (auto [rowB_j, j] : sparseRange(row))
{
// Y[i][j] += rowA_i' * X * rowB_j;
typename Matrix::block_type& Y_ij = Y[i][j];
for (int k=0; k<M; ++k)
// Compute matrix entry in upper right triangle part of Y
if (i<=j)
{
const auto& rowA_i_k = rowA_i[k];
for(int l=0; l<M; ++l)
if constexpr (Std::is_detected_v<Impl::isUnblockedDiagonalMatrix_t, decltype(rowB_j)>)
rowB_j.leftMultiplyTransposed(precomp, Y_ij);
else
{
const auto& rowB_j_l = rowB_j[l];
for (auto [rowA_i_k_ii, ii] : sparseRange(rowA_i_k))
for (auto [rowB_j_l_jj, jj] : sparseRange(rowB_j_l))
Y_ij[ii][jj] += rowA_i_k_ii * X[k][l] * rowB_j_l_jj;
for (int ii=0; ii<Y_ij.rows; ++ii)
for (int jj=0; jj<Y_ij.cols; ++jj)
for (int l=0; l<M; ++l)
Y_ij[ii][jj] += precomp[l][ii] * rowB_j[l][jj];
}
}
else
// Do not compute matrix entry in left triangle part of Y.
// Copy it from the upper right part instead.
{
typename Matrix::block_type& Y_ji = Y[j][i];
for (int ii=0; ii<Y_ji.rows; ++ii)
for (int jj=0; jj<Y_ji.cols; ++jj)
Y_ij[ii][jj] = Y_ji[jj][ii];
}
}
}
}
const F& f_;
const BitVector& ignore_;
const BitVector* ignore_;
Vector negativeGradient_;
Matrix hessian_;
......
......@@ -7,3 +7,5 @@ dune_add_test(SOURCES damagedductileelasticenergydensitytest.cc)
dune_add_test(SOURCES integralfunctionaltest.cc)
dune_add_test(SOURCES largestraindamagedelasticenergydensitytest.cc)
dune_add_test(SOURCES unblockeddiagonalmatrixtest.cc)
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#include "config.h"
#include <iostream>
#include <dune/common/exceptions.hh>
#include <dune/common/fmatrix.hh>
#include <dune/matrix-vector/transpose.hh>
#include <dune/fracture-phasefields/unblockeddiagonalmatrix.hh>
using namespace Dune;
template <class TestMatrix, class ReferenceMatrix>
bool compareBehavior(TestMatrix testMatrix,
ReferenceMatrix referenceMatrix,
const std::vector<size_t>& rowBlocks,
const std::vector<size_t>& colBlocks)
{
bool passed = true;
// Test whether the matrix dimensions match
if (testMatrix.rows != referenceMatrix.rows)
DUNE_THROW(Exception, "testMatrix does not have the correct number of rows!");
if (testMatrix.cols != referenceMatrix.cols)
DUNE_THROW(Exception, "testMatrix does not have the correct number of columns!");
// Test whether method umv computes the right thing
{
FieldVector<double, testMatrix.cols> v;
for (size_t i=0; i<v.size(); ++i)
{
v = 0;
v[i] = 1.0;
FieldVector<double, testMatrix.rows> testResult(0);
FieldVector<double, referenceMatrix.rows> referenceResult(0);
testMatrix.umv(v,testResult);
referenceMatrix.umv(v,referenceResult);
if ( (testResult - referenceResult).infinity_norm() > 1e-6 )
{
std::cout << "Method umv did not compute the correct result!" << std::endl;
std::cout << "Result: " << testResult << std::endl;
std::cout << "Expected: " << referenceResult << std::endl;
passed = false;
}
}
}
// Test whether method usmv computes the right thing
{
FieldVector<double, testMatrix.cols> v;
for (size_t i=0; i<v.size(); ++i)
{
v = 0;
v[i] = 1.0;
FieldVector<double, testMatrix.rows> testResult(0);
FieldVector<double, referenceMatrix.rows> referenceResult(0);
testMatrix.usmv(0.5, v,testResult);
referenceMatrix.usmv(0.5, v,referenceResult);
if ( (testResult - referenceResult).infinity_norm() > 1e-6 )
{
std::cout << "Method usmv did not compute the correct result!" << std::endl;
std::cout << "Result: " << testResult << std::endl;
std::cout << "Expected: " << referenceResult << std::endl;
passed = false;
}
}
}
// Test whether method umtv computes the right thing
{
FieldVector<double, testMatrix.rows> v;
for (size_t i=0; i<v.size(); ++i)
{
v = 0;
v[i] = 1.0;
FieldVector<double, testMatrix.cols> testResult(0);
FieldVector<double, referenceMatrix.cols> referenceResult(0);
testMatrix.umtv(v,testResult);
referenceMatrix.umtv(v,referenceResult);
if ( (testResult - referenceResult).infinity_norm() > 1e-6 )
{
std::cout << "Method umtv did not compute the correct result!" << std::endl;