newton_solver.c

一个用来实现偏微分方程中网格的计算库

                  << current_residual << std::endl;

      // Make sure our linear tolerance is low enough
      current_linear_tolerance = std::min (current_linear_tolerance,
        current_residual * linear_tolerance_multiplier);

      // But don't let it be too small
      if (current_linear_tolerance < minimum_linear_tolerance)
        {
          current_linear_tolerance = minimum_linear_tolerance;
        }

      // At this point newton_iterate is the current guess, and
      // linear_solution is now about to become the NEGATIVE of the next
      // Newton step.

      // Our best initial guess for the linear_solution is zero!
      linear_solution.zero();

      if (!quiet)
        std::cout << "Linear solve starting, tolerance " 
                  << current_linear_tolerance << std::endl;

      // Solve the linear system.  Two cases:
      const std::pair<unsigned int, Real> rval =
        (_system.have_matrix("Preconditioner")) ?
      // 1.) User-supplied preconditioner
        linear_solver->solve (matrix, _system.get_matrix("Preconditioner"),
                              linear_solution, rhs, current_linear_tolerance,
                              max_linear_iterations) :
      // 2.) Use system matrix for the preconditioner
        linear_solver->solve (matrix, linear_solution, rhs,
                              current_linear_tolerance, 
                              max_linear_iterations);
      // We may need to localize a parallel solution
      _system.update ();
      // The linear solver may not have fit our constraints exactly
#ifdef ENABLE_AMR
      _system.get_dof_map().enforce_constraints_exactly(_system, &linear_solution);
#endif

      const unsigned int linear_steps = rval.first;
      libmesh_assert(linear_steps <= max_linear_iterations);
      _inner_iterations += linear_steps;

      const bool linear_solve_finished = 
        !(linear_steps == max_linear_iterations);

      if (!quiet)
        std::cout << "Linear solve finished, step " << linear_steps
                  << ", residual " << rval.second
                  << std::endl;

      // Compute the l2 norm of the nonlinear update
      Real norm_delta = linear_solution.l2_norm();

      if (!quiet)
        std::cout << "Trying full Newton step" << std::endl;
      // Take a full Newton step
      newton_iterate.add (-1., linear_solution);
      newton_iterate.close();

      // Check residual with full Newton step
      _system.assembly(true, false);

      rhs.close();
      current_residual = rhs.l2_norm();
      if (!quiet)
        std::cout << "  Current Residual: "
                  << current_residual << std::endl;

// A potential method for avoiding oversolving?
/*
      Real predicted_absolute_error =
        current_residual * norm_delta / last_residual;

      Real predicted_relative_error =
        predicted_absolute_error / max_solution_norm;

      std::cout << "Predicted absolute error = " <<
        predicted_absolute_error << std::endl;

      std::cout << "Predicted relative error = " <<
        predicted_relative_error << std::endl;
*/

      // don't fiddle around if we've already converged
      if (test_convergence(current_residual, norm_delta,
                           linear_solve_finished))
        {
          if (!quiet)
            print_convergence(_outer_iterations, current_residual,
                              norm_delta, linear_solve_finished);
          _outer_iterations++;
          break; // out of _outer_iterations for loop
        }

      // otherwise, backtrack if necessary
      Real steplength =
        this->line_search(std::sqrt(TOLERANCE),
                          last_residual, current_residual,
                          newton_iterate, linear_solution);
      norm_delta *= steplength;

      // Check to see if backtracking failed,
      // and break out of the nonlinear loop if so...
      if (_solve_result == DiffSolver::DIVERGED_BACKTRACKING_FAILURE)
        {
          _outer_iterations++;
	  break; // out of _outer_iterations for loop
        }

      // Compute the l2 norm of the whole solution
      norm_total = newton_iterate.l2_norm();

      max_solution_norm = std::max(max_solution_norm, norm_total);

      // Print out information for the 
      // nonlinear iterations.
      if (!quiet)
        std::cout << "  Nonlinear step: |du|/|u| = "
                  << norm_delta / norm_total
                  << ", |du| = " << norm_delta
                  << std::endl;

      // Terminate the solution iteration if the difference between
      // this iteration and the last is sufficiently small.
      if (!quiet)
        print_convergence(_outer_iterations, current_residual,
                          norm_delta / steplength,
                          linear_solve_finished);
      if (test_convergence(current_residual, norm_delta / steplength,
                           linear_solve_finished))
        {
          _outer_iterations++;
	  break; // out of _outer_iterations for loop
        }

      if (_outer_iterations >= max_nonlinear_iterations - 1)
        {
          std::cout << "  Nonlinear solver FAILED TO CONVERGE by step "
                    << _outer_iterations << " with norm "
                    << norm_total << std::endl;
          if (continue_after_max_iterations)
	    {
	      _solve_result = DiffSolver::DIVERGED_MAX_NONLINEAR_ITERATIONS;
	      std::cout << "  Continuing anyway..." << std::endl;
	    }
          else
	    {
	      libmesh_error();
	    }
          continue;
        }
    } // end nonlinear loop

  // The linear solver may not have fit our constraints exactly
#ifdef ENABLE_AMR
  _system.get_dof_map().enforce_constraints_exactly(_system);
#endif

  // We may need to localize a parallel solution
  _system.update ();

  STOP_LOG("solve()", "NewtonSolver");

  // Make sure we are returning something sensible as the
  // _solve_result.
  libmesh_assert (_solve_result != DiffSolver::INVALID_SOLVE_RESULT);
  
  return _solve_result;
}



bool NewtonSolver::test_convergence(Real current_residual,
                                    Real step_norm,
                                    bool linear_solve_finished)
{
  // We haven't converged unless we pass a convergence test
  bool has_converged = false;

  // Is our absolute residual low enough?
  if (current_residual < absolute_residual_tolerance)
    {
      _solve_result |= CONVERGED_ABSOLUTE_RESIDUAL;
      has_converged = true;
    }
  
  // Is our relative residual low enough?
  if ((current_residual / max_residual_norm) <
      relative_residual_tolerance)
    {
      _solve_result |= CONVERGED_RELATIVE_RESIDUAL;
      has_converged = true;
    }
  
  // For incomplete linear solves, it's not safe to test step sizes
  if (!linear_solve_finished)
    {
      return has_converged;
    }
  
  // Is our absolute Newton step size small enough?
  if (step_norm < absolute_step_tolerance)
    {
      _solve_result |= CONVERGED_ABSOLUTE_STEP;
      has_converged = true;
    }
  
  // Is our relative Newton step size small enough?
  if (step_norm / max_solution_norm <
      relative_step_tolerance)
    {
      _solve_result |= CONVERGED_RELATIVE_STEP;
      has_converged = true;
    }
  
  return has_converged;
}


void NewtonSolver::print_convergence(unsigned int step_num,
                                     Real current_residual,
                                     Real step_norm,
                                     bool linear_solve_finished)
{
  // Is our absolute residual low enough?
  if (current_residual < absolute_residual_tolerance)
    {
      std::cout << "  Nonlinear solver converged, step " << step_num
                << ", residual " << current_residual
                << std::endl;
    }
  else if (absolute_residual_tolerance)
    {
      std::cout << "  Nonlinear solver current_residual "
                << current_residual << " > "
                << (absolute_residual_tolerance) << std::endl;
    }

  // Is our relative residual low enough?
  if ((current_residual / max_residual_norm) <
      relative_residual_tolerance)
    {
      std::cout << "  Nonlinear solver converged, step " << step_num
                << ", residual reduction "
                << current_residual / max_residual_norm
                << " < " << relative_residual_tolerance
                << std::endl;
    }
  else if (relative_residual_tolerance)
    {
      if (!quiet)
        std::cout << "  Nonlinear solver relative residual "
                  << (current_residual / max_residual_norm)
                  << " > " << relative_residual_tolerance
                  << std::endl;
    }

  // For incomplete linear solves, it's not safe to test step sizes
  if (!linear_solve_finished)
    return;

  // Is our absolute Newton step size small enough?
  if (step_norm < absolute_step_tolerance)
    {
      std::cout << "  Nonlinear solver converged, step " << step_num
                << ", absolute step size "
                << step_norm
                << " < " << absolute_step_tolerance
                << std::endl;
    }
  else if (absolute_step_tolerance)
    {
      std::cout << "  Nonlinear solver absolute step size "
                << step_norm
                << " > " << absolute_step_tolerance
                << std::endl;
    }

  // Is our relative Newton step size small enough?
  if (step_norm / max_solution_norm <
      relative_step_tolerance)
    {
      std::cout << "  Nonlinear solver converged, step " << step_num
                << ", relative step size "
                << (step_norm / max_solution_norm)
                << " < " << relative_step_tolerance
                << std::endl;
    }
  else if (relative_step_tolerance)
    {
      std::cout << "  Nonlinear solver relative step size "
                << (step_norm / max_solution_norm)
                << " > " << relative_step_tolerance
                << std::endl;
    }
}