Method solve

cpp/dolfinx/nls/NewtonSolver.cpp:149–276

-----------------------------------------------------------------------------

}
//-----------------------------------------------------------------------------
std::pair<int, bool> nls::petsc::NewtonSolver::solve(Vec x)
{
  // Reset iteration counts
  _iteration = 0;
  _krylov_iterations = 0;
  _residual = -1;
  _residual0 = 0;

  if (!_fnF)
  {
    throw std::runtime_error("Function for computing residual vector has not "
                             "been provided to the NewtonSolver.");
  }

  if (!_fnJ)
  {
    throw std::runtime_error("Function for computing Jacobian has not "
                             "been provided to the NewtonSolver.");
  }

  if (_system)
    _system(x);
  assert(_b);
  _fnF(x, _b);

  // Check convergence
  bool newton_converged = false;
  if (convergence_criterion == "residual")
    std::tie(_residual, newton_converged) = this->_converged(*this, _b);
  else if (convergence_criterion == "incremental")
  {
    // We need to do at least one Newton step with the ||dx||-stopping
    // criterion
    newton_converged = false;
  }
  else
  {
    throw std::runtime_error("Unknown convergence criterion: "
                             + convergence_criterion);
  }

  // FIXME: check that this is efficient if A and/or P are unchanged
  // Set operators
  if (_matP)
    _solver.set_operators(_matJ, _matP);
  else
    _solver.set_operators(_matJ, _matJ);

  if (!_dx)
    MatCreateVecs(_matJ, &_dx, nullptr);

  // Start iterations
  while (!newton_converged and _iteration < max_it)
  {
    // Compute Jacobian
    assert(_matJ);
    _fnJ(x, _matJ);