test-gradient-sample-nonsmooth.cc

/* Checks if the descent direction is computed correctly
   using the analytic solution */

#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

#include <cassert>

#include <dune/common/shared_ptr.hh>

#include <dune/tectonic/ellipticenergy.hh>

#include "test-gradient-method-nicefunction.hh"
#include "test-gradient-method-helper.hh"

int main() {
  int const dim = 2;
  typedef Dune::EllipticEnergy<dim> Functional;
  typedef Functional::SmallMatrix SmallMatrix;
  typedef Functional::SmallVector SmallVector;

  SmallMatrix const A = { { 3, 1.5 }, { 1.5, 4 } };
  SmallVector const b = { 1, 2 };

  auto const f = Dune::make_shared<Dune::SampleFunction<2> const>();
  auto const phi = Dune::make_shared<Functional::NonlinearityType const>(f);
  Functional const J(A, b, phi);

  SmallVector start;
  /*
    for x = b/|b|:

    j_(x)
    = Ax - b + 1/|x| x
    = 1/sqrt(5) (6, 9.5) - (1, 2) + 1/sqrt(5) (1, 2)
    = (7/sqrt(5) - 1, 11.5/sqrt(5) - 2)

    j+(x)
    = Ax - b + 2/|x| x
    = 1/sqrt(5) (6, 9.5) - (1, 2) + 2/sqrt(5) (1, 2)
    = (8/sqrt(5) - 1, 13.5/sqrt(5) - 2)
  */
  {
    start = b;
    start /= (b.two_norm() + 1e-12); // Make sure the norm is below 1;
    assert(start.two_norm() < 1);

    SmallVector numerical_descent;
    J.descentDirection(start, numerical_descent);
    SmallVector analytic_descent = { -(7 / std::sqrt(5) - 1),
                                     -(11.5 / std::sqrt(5) - 2) };
    assert(two_distance<dim>(analytic_descent, numerical_descent) < 1e-10);

    functionTester(J, start, 6);
  }
  {
    start = b;
    start /= (b.two_norm() - 1e-12); // Make sure the norm is above 1;
    assert(start.two_norm() > 1);

    SmallVector numerical_descent;
    J.descentDirection(start, numerical_descent);
    SmallVector analytic_descent = { -(8 / std::sqrt(5) - 1),
                                     -(13.5 / std::sqrt(5) - 2) };
    assert(two_distance<dim>(analytic_descent, numerical_descent) < 1e-10);

    functionTester(J, start, 6);
  }
}