Commit 749a5bbe authored by Gräser, Carsten's avatar Gräser, Carsten
Browse files

Merge branch 'feature/ensure-monotonocity' into 'master'

Ensure and track monotonicity

See merge request !51
parents b7e0f05b 81ace1e1
Pipeline #8774 passed with stage
in 4 minutes and 13 seconds
......@@ -384,13 +384,31 @@ public:
// Premultiply evaluationMatrix and direction vector.
// This is used later in the subDifferential method.
f_->evaluationMatrix().mv(*direction_, evaluationMatrixTimesDirection_);
f_->evaluationMatrix().mv(*origin_, evaluationMatrixTimesOrigin_);
}
Range operator()(const Scalar& v) const
/** \brief Compute the value of the functional a given point
*
* \param z Where to evaluate the subdifferential
* (The coordinate is a single number because we are in a one-dimensional space).
*/
Range operator()(const Scalar& z) const
{
DUNE_THROW(NotImplemented, "DirectionalRestriction::operator()");
Range result(0.0);
return result;
// The global position of where the directional derivative
// is to be evaluated is v = origin_ + z*direction. However,
// in the following we only need B*(origin + z*direction) = Bv
auto BXv = evaluationMatrixTimesDirection_;
BXv *= z;
BXv += evaluationMatrixTimesOrigin_;
auto&& density = f_->density();
auto&& weights = f_->quadratureWeights();
Range value = 0;
for (std::size_t i=0; i<f_->evaluationMatrix().N(); ++i)
value += density(BXv[i])*weights[i];
return value;;
}
/** \brief The set of points where the functional takes finite values
......@@ -415,26 +433,21 @@ public:
*/
Solvers::Interval<double> subDifferential(Scalar z) const
{
// The global position of where the directional derivative is to be evaluated
auto evaluationPoint = *origin_;
evaluationPoint.axpy(z,*direction_);
GlobalVector gradient;
Solvers::resizeInitializeZero(gradient,(*origin_));
// The global position of where the directional derivative
// is to be evaluated is v = origin_ + z*direction. However,
// in the following we only need B*(origin + z*direction) = Bv
auto BXv = evaluationMatrixTimesDirection_;
BXv *= z;
BXv += evaluationMatrixTimesOrigin_;
typename Functional::Density::Domain BiXv;
auto densityDerivative = derivative(f_->density());
auto&& weights = f_->quadratureWeights();
double directionalDerivative = 0;
for (std::size_t i=0; i<f_->evaluationMatrix().N(); ++i)
{
BiXv = 0.0;
Functional::addRowXVector(f_->evaluationMatrix()[i], evaluationPoint, BiXv);
auto densityDerivative = derivative(f_->density());
auto Dgamma = densityDerivative(BiXv);
Dgamma *= f_->quadratureWeights()[i];
auto Dgamma = densityDerivative(BXv[i]);
Dgamma *= weights[i];
directionalDerivative += Dgamma * evaluationMatrixTimesDirection_[i];
}
......@@ -451,15 +464,26 @@ public:
s = subDifferential(z);
}
const GlobalVector& origin() const
{
return *origin_;
}
const GlobalVector& direction() const
{
return *direction_;
}
protected:
const Functional* f_;
const GlobalVector* origin_;
const GlobalVector* direction_;
// The evaluation matrix multiplied by the direction vector
// The evaluation matrix multiplied by the origin and direction vector
// This product is used in the subDifferential method,
// and it pays to precompute it.
BlockVector<typename Functional::Density::Domain> evaluationMatrixTimesDirection_;
BlockVector<typename Functional::Density::Domain> evaluationMatrixTimesOrigin_;
};
......
......@@ -77,19 +77,38 @@
using namespace Dune;
// In principle this could be implemented as a simple lambda,
// but this would lead to an incomplete type warning.
template <class SSF>
class ScalarSumFunctionalDerivative
{
public:
ScalarSumFunctionalDerivative(const SSF& ssf) : ssf_(ssf) {}
double operator()(const double& z) const {
return ssf_.subDifferential(z).projectIn(0);
}
const SSF& ssf_;
};
/** \brief Sum of two functionals defined on the set of real numbers
*
* \tparam DamageFunctional Expected to model BoxConstrainedQuadraticFunctional
* \tparam DamagedElasticityFunctional Expected to be smooth but nonquadratic
*/
template <class DamageFunctional, class DamagedElasticityFunctional>
template <class DamageFunctional, class DamagedElasticityFunctional, class GlobalFunctional>
class ScalarSumFunctional
{
public:
ScalarSumFunctional(const DamageFunctional& damageFunctional,
const DamagedElasticityFunctional& damagedElasticityFunctional)
const DamagedElasticityFunctional& damagedElasticityFunctional, const GlobalFunctional& J)
: damageFunctional_(damageFunctional),
damagedElasticityFunctional_(damagedElasticityFunctional)
damagedElasticityFunctional_(damagedElasticityFunctional),
J_(J)
{}
Solvers::Interval<double> domain() const
......@@ -125,8 +144,27 @@ public:
return dJ;
}
double operator()(const double& z) const
{
const auto& origin = damagedElasticityFunctional_.origin();
const auto& direction = damagedElasticityFunctional_.direction();
temp = direction;
temp *= z;
temp += origin;
return J_(temp);
}
friend auto derivative(const ScalarSumFunctional& sumFunctional)
{
// Using a lambda would lead to an incomplete type warning
return ScalarSumFunctionalDerivative(sumFunctional);
}
const DamageFunctional& damageFunctional_;
const DamagedElasticityFunctional& damagedElasticityFunctional_;
const GlobalFunctional& J_;
mutable typename DamagedElasticityFunctional::GlobalVector temp;
};
......@@ -822,10 +860,17 @@ int main (int argc, char *argv[]) try
using Functional = TNNMG::SumFunctional<QuadraticPartFunctional,DissipationFunctional>;
auto J = Functional(quadraticPartFunctional,dissipationFunctional);
auto bisectionSolver = [](auto& xi, const auto& restriction, const auto& ignoreI) {
std::size_t lineSearchFunctionalEvaluation = 0;
std::size_t lineSearchDerivativeEvaluation = 0;
double postGSEnergy = 0;
double postMGEnergy = 0;
double oldPostMGEnergy = std::numeric_limits<double>::max();
auto bisectionSolver = [&](auto& xi, const auto& restriction, const auto& ignoreI) {
ScalarSumFunctional sumFunctional(std::get<0>(restriction.functions()), std::get<1>(restriction.functions()), J);
ScalarSumFunctional<decltype(std::get<0>(restriction.functions())),
decltype(std::get<1>(restriction.functions()))> sumFunctional(std::get<0>(restriction.functions()), std::get<1>(restriction.functions()));
#if 0
Bisection bisection;
......@@ -837,73 +882,42 @@ int main (int argc, char *argv[]) try
// It should probably by one of the following:
// xi = bisection.minimize(sumFunctional, 0, 0, bisectionSteps);
// xi = bisection.minimize(sumFunctional, 1, 0, bisectionSteps);
#elif 0
auto f = [&](auto z) {
return sumFunctional.subDifferential(z).projectIn(0);
};
auto I = sumFunctional.domain();
I[0] = 0;
#else
auto fI = I;
fI[0] = f(I[0]);
auto dSumFunctional = derivative(sumFunctional);
auto dom = sumFunctional.domain();
if (f(0) < 0)
{
fI[1] = f(I[1]);
lineSearchFunctionalEvaluation = 0;
lineSearchDerivativeEvaluation = 0;
// Absolute error tolerance 0
// Relative error tolerance 0.1
// If possible, start from interval [1,3].
NewBisection bisection(0, 0.2);
if (I[1] >= 1)
{
auto f1 = f(1);
if (f1 <= 0)
{
I[0] = 1;
fI[0] = f1;
}
}
if (I[1] >= 3)
{
auto f3 = f(3);
if (f3 >= 0)
{
I[1] = 3;
fI[1] = f3;
}
}
xi = bisection.solve(f, I, fI, 0);
}
else
xi = 0;
#else
auto g = [&](auto z) {
++lineSearchFunctionalEvaluation;
return sumFunctional(z);
};
auto f = [&](auto z) {
return sumFunctional.subDifferential(z).projectIn(0);
// auto result = sumFunctional.subDifferential(z).projectIn(0);
// std::cout << "f(" << z << ") = " << result << std::endl;
// return result;
auto dg = [&](auto z) {
++lineSearchDerivativeEvaluation;
return dSumFunctional(z);
};
auto contains = [](const Solvers::Interval<double>& J, double t) {
return (J[0] <= t) and (t <= J[1]);
};
auto dom = sumFunctional.domain();
postGSEnergy = g(0);
auto I = dom;
auto fI = I;
auto dgI = I;
I[0] = 0;
fI[0] = f(I[0]);
dgI[0] = dg(I[0]);
// Ensure that we have a descent direction. Otherwise use step length = 0
if (fI[0] >= 0)
if (dgI[0] >= 0)
{
xi = 0;
postMGEnergy = postGSEnergy;
return;
}
......@@ -912,38 +926,46 @@ int main (int argc, char *argv[]) try
if (contains(I, 1))
{
// Start from I[1]=1 and successively double the step length
// until we have found an upper bound with f(I[1])>0.
// until we have found an upper bound with dg(I[1])>0.
// While doing so, replace the lower bound I[0] by I[1]
// whenever f(I[1])<0. That way I[0] will always be a lower
// bound with f(I[0])<0.
// In the special case f(I[1])=0, we're done and can instantly
// whenever dg(I[1])<0. That way I[0] will always be a lower
// bound with dg(I[0])<0.
// In the special case dg(I[1])=0, we're done and can instantly
// return the found root I[1].
I[1] = 1;
fI[1] = f(I[1]);
while(fI[1] < 0)
dgI[1] = dg(I[1]);
while(dgI[1] < 0)
{
I[0] = I[1];
fI[0] = fI[1];
dgI[0] = dgI[1];
I[1] = dom.projectIn(I[1]*2);
fI[1] = f(I[1]);
}
if (fI[1] == 0)
{
xi = I[1];
return;
dgI[1] = dg(I[1]);
}
}
else
fI[1] = f(I[1]);
// Start bisection with:
// * absolute error tolerance 0
// * relative error tolerance 0.2
// * residual tolerance 0
// * initial interval guess I
// * x_ref=0 as refenence value to determine the relative error
NewBisection bisection(0, 0.2);
xi = bisection.solve(f, I, fI, 0);
dgI[1] = dg(I[1]);
if (dgI[1] == 0)
xi = I[1];
else
{
// Start bisection with:
// * absolute error tolerance 0
// * relative error tolerance 0.2
// * residual tolerance 0
// * initial interval guess I
// * x_ref=0 as refenence value to determine the relative error
NewBisection bisection(0, 0.2);
xi = bisection.solve(dg, I, dgI, 0);
}
postMGEnergy = g(xi);
while(postMGEnergy>postGSEnergy)
{
xi *= .5;
postMGEnergy = g(xi);
}
#endif
};
......@@ -993,21 +1015,44 @@ int main (int argc, char *argv[]) try
solver.addCriterion(
[&](){
return Dune::formatString(" % 12d", step.linearization().truncated().count());
return Dune::formatString("% 10d", step.linearization().truncated().count());
},
" truncated ");
" truncated");
solver.addCriterion(
[&](){
return Dune::formatString(" % 12.5e", step.lastDampingFactor());
},
" step size ");
" step size");
solver.addCriterion(
[&](){
return Dune::formatString(" % 3d % 3d", lineSearchDerivativeEvaluation, lineSearchFunctionalEvaluation);
},
" bisection");
solver.addCriterion(
[&](){
return Dune::formatString(" % 22.15e", postMGEnergy);
},
" energy");
solver.addCriterion(
[&](){
return Dune::formatString(" % 22.15e", postGSEnergy);
},
" post GS energy");
solver.addCriterion(
[&](){
return Dune::formatString(" % 12.5e", J(x));
double relEnergyViolation = std::max(postGSEnergy-oldPostMGEnergy, 0.0)/std::fabs(oldPostMGEnergy);
oldPostMGEnergy = postMGEnergy;
if (relEnergyViolation==0)
return std::string(" ");
else
return Dune::formatString(" % 12.5e X", relEnergyViolation);
},
" energy ");
" rel.GS.mono.err");
solver.preprocess();
Timer timer;
......
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