compute_state.cc

#include <cassert>

#include "LambertW.h"

#include <gsl/gsl_sf_lambert.h>

#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
#include <dune/fufem/interval.hh>

#include <dune/tectonic/mydirectionalconvexfunction.hh>
#include <dune/tnnmg/problem-classes/bisection.hh>

#include "compute_state.hh"

// The nonlinearity exp(-x)
class DecayingExponential {
public:
  typedef Dune::FieldVector<double, 1> VectorType;
  typedef Dune::FieldMatrix<double, 1, 1> MatrixType;

  double operator()(VectorType const &u) const { return std::exp(-u[0]); }

  void directionalSubDiff(VectorType const &u, VectorType const &v,
                          Interval<double> &D) const {
    D[0] = D[1] = v[0] * (-std::exp(-u[0]));
  }

  void directionalDomain(VectorType const &, VectorType const &,
                         Interval<double> &dom) const {
    dom[0] = -std::numeric_limits<double>::max();
    dom[1] = std::numeric_limits<double>::max();
  }
};

double compute_state_update_bisection(double h, double unorm, double L,
                                      double old_state) {
  MyDirectionalConvexFunction<DecayingExponential> const J(
      1.0 / h, (old_state - unorm / L) / h, DecayingExponential(), 0, 1);
  int bisectionsteps = 0;
  Bisection const bisection(0.0, 1.0, 1e-12, true, 0);    // TODO
  return bisection.minimize(J, 0.0, 0.0, bisectionsteps); // TODO
}

double compute_state_update_lambert(double h, double unorm, double L,
                                    double old_state) {
  double const rhs = unorm / L - old_state;
  return LambertW(0, h * std::exp(rhs)) - rhs;
}

double compute_state_update_lambert_gsl(double h, double unorm, double L,
                                        double old_state) {
  double const rhs = unorm / L - old_state;
  return gsl_sf_lambert_W0(h * std::exp(rhs)) - rhs;
}

double compute_state_update(double h, double unorm, double L,
                            double old_state) {
  double ret1 = compute_state_update_bisection(h, unorm, L, old_state);
  assert(std::abs(1.0 / h * ret1 - (old_state - unorm / L) / h -
                  std::exp(-ret1)) < 1e-10);

  double ret2 = compute_state_update_lambert(h, unorm, L, old_state);
  assert(std::abs(1.0 / h * ret2 - (old_state - unorm / L) / h -
                  std::exp(-ret2)) < 1e-10);

  double ret3 = compute_state_update_lambert_gsl(h, unorm, L, old_state);
  assert(std::abs(1.0 / h * ret3 - (old_state - unorm / L) / h -
                  std::exp(-ret3)) < 1e-10);

  assert(std::abs(ret1 - ret2) < 1e-14);
  assert(std::abs(ret1 - ret3) < 1e-14);

  return ret1;
}