Commit 81ace1e1 authored by Gräser, Carsten's avatar Gräser, Carsten
Browse files

Ensure and track monotonicity

This is a collection of many interacting patches that help to ensure
monotonicity of the TNNMG iteration and add more diagnostic output to
track the latter:

* Implement `ScalarSumFunctional::operator()` in a way that is
  consistent with the global energy evaluation. Unfortunately
  this means we cannot use the evaluation of the individual
  constrained functionals.
* Implement `derivative(ScalarSumFunctional)`. The functional
  is in general not differentiable, hence we return the minimal
  norm entry of the subdifferential or, equivalently, the projection
  of zero into the subdifferential. Finding a root of the
  subdifferential is equivalent to find a root of this single-valued
  map. Furthermore a subdifferential-aware bisection is equivalent
  to a bisection for the single-valued map.
* To ensure monotonicity of the coarse grid correction add a
  postprocessing step to the line search: After an approximate
  root of the derivative has been found, check if this decreases energy.
  If not, then half the step size until it does. Since we have a descent
  direction this must terminate. This step is necessary because the
  functional is not convex.
* When comparing energies take care to always evaluate it exactly the
  same way, to avoid non-monotonicity due to inconsistent evaluation
  algorithms.
* Track energy decrease of nonlinear smoother and coarse grid
  correction. While the latter is guaranteed by the above changes,
  it may happen that the smoother is not monotone - either because
  of local solves with uncontrolled inexactness (called exact in our terminalogy, arg...)
  or because of round-off-errors.
* Add more diagnostic output to the table:
  * Number of subdifferential evaluations during line search.
    Hint on interpretation: This will be at least 2+#bisectionSteps.
  * Number of functional evaluations during line search.
    Hint on interpretation: This will be 2+#reductionSteps,
    where #reductionSteps is the number of halving steps done
    to ensure monotonicity (see above).
  * Energy after TNNMG step with 16 digits.
  * Energy after nonlinear pre-smoothing with 16 digits.
  * If the nonlinear pre-smoother is not monotone:
    Relative monotonicity violation of the smoother with an additional
    `X` marker for easier grepping.
parent d1a5767b
Pipeline #8773 passed with stage
in 4 minutes and 17 seconds
......@@ -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