Commit 3c44697c authored by Gräser, Carsten's avatar Gräser, Carsten
Browse files

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

Add damping to local Newton solver

See merge request !52
parents 749a5bbe 0f1afe6b
Pipeline #8791 passed with stage
in 4 minutes and 10 seconds
......@@ -51,8 +51,40 @@ class BrittleFractureLocalSolver
// solve for upper left block
FieldVector<double,dim> correction;
H00.solve(correction, g0);
// Extend vector to full size
FieldVector<double,dim+1> direction;
for(std::size_t i=0; i<dim; ++i)
direction[i] = correction[i];
direction[dim] = 0;
// Constrained functional along correction
auto f = [&](auto t) {
auto z = x + t*direction;
return quadraticFunctional(z) + nonlinearFunctional(z);
};
// Reduce step size if needed to ensure monotonicity
double t = 1;
auto f0 = f(0);
auto ft = f(t);
while (std::max(ft-f0, 0.0) > 1e-14 * std::fabs(f0))
{
t *= .5;
ft = f(t);
}
// if (t<1.0)
// {
// std::cout << "++++++++++++++++++++++++++++++++++++++++++++++++++++" << std::endl;
// std::cout << " |direction| = " << direction.infinity_norm() << std::endl;
// std::cout << " |f(0)-f(1)|/|f(0)| = " << std::fabs((f(1)-f0)/f0) << std::endl;
// std::cout << " f(1) = " << f(1) << std::endl;
// std::cout << " f(" << t << ") = " << f(t) << std::endl;
// }
// Apply damped correction
for(std::size_t i=0; i<dim; ++i)
x[i] += correction[i];
x[i] += t*correction[i];
}
template <class SumFunctional, class BitVector>
......
......@@ -900,10 +900,6 @@ int main (int argc, char *argv[]) try
return dSumFunctional(z);
};
auto contains = [](const Solvers::Interval<double>& J, double t) {
return (J[0] <= t) and (t <= J[1]);
};
postGSEnergy = g(0);
......@@ -923,7 +919,7 @@ int main (int argc, char *argv[]) try
// If the domain contains 1, try to find a better initial interval guess.
// Otherwise stay with the domain as initial guess.
if (contains(I, 1))
if ((I[0] <= 1.0) and (1.0 <= I[1]))
{
// Start from I[1]=1 and successively double the step length
// until we have found an upper bound with dg(I[1])>0.
......@@ -960,7 +956,7 @@ int main (int argc, char *argv[]) try
}
postMGEnergy = g(xi);
while(postMGEnergy>postGSEnergy)
while (std::max(postMGEnergy-postGSEnergy, 0.0) > 1e-14 * std::fabs(postGSEnergy))
{
xi *= .5;
postMGEnergy = g(xi);
......
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