← Back to branch summary

~njansson/dolfin/hpc

« back to all changes in this revision

Viewing changes to src/kernel/nls/NewtonSolver.cpp

  • Committer: Johannes Ring
  • Date: 2008-03-05 22:43:06 UTC
  • Revision ID: johannr@simula.no-20080305224306-2npsdyhfdpl2esji
The BIG commit!

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/kernel/nls/NewtonSolver.cpp
1
 
// Copyright (C) 2005-2007 Garth N. Wells.
2
 
// Licensed under the GNU LGPL Version 2.1.
3
 
//
4
 
// Modified by Anders Logg, 2005-2006.
5
 
//
6
 
// First added:  2005-10-23
7
 
// Last changed: 2007-05-15
8
 
 
9
 
#include <dolfin/NewtonSolver.h>
10
 
#include <dolfin/NonlinearProblem.h>
11
 
#include <dolfin/LUSolver.h>
12
 
#include <dolfin/KrylovSolver.h>
13
 
#include <iostream>
14
 
 
15
 
using namespace dolfin;
16
 
 
17
 
//-----------------------------------------------------------------------------
18
 
NewtonSolver::NewtonSolver() : Parametrized()
19
 
{
20
 
  solver = new LUSolver();
21
 
}
22
 
//-----------------------------------------------------------------------------
23
 
/*
24
 
#ifdef HAVE_PETSC_H
25
 
NewtonSolver::NewtonSolver(Matrix::Type matrix_type) : Parametrized()
26
 
{
27
 
  solver = new LUSolver();
28
 
  // FIXME: Need to select appropriate PETSc matrix
29
 
  //A = new Matrix(matrix_type);
30
 
  A = new Matrix;
31
 
}
32
 
#endif
33
 
*/
34
 
//-----------------------------------------------------------------------------
35
 
NewtonSolver::NewtonSolver(KrylovMethod method, Preconditioner pc)
36
 
  : Parametrized()
37
 
{
38
 
  solver = new KrylovSolver(method, pc);
39
 
}
40
 
//-----------------------------------------------------------------------------
41
 
NewtonSolver::~NewtonSolver()
42
 
{
43
 
  delete solver; 
44
 
}
45
 
//-----------------------------------------------------------------------------
46
 
dolfin::uint NewtonSolver::solve(NonlinearProblem& nonlinear_problem, Vector& x)
47
 
{
48
 
  const uint maxit= get("Newton maximum iterations");
49
 
 
50
 
  // FIXME: add option to compute F(u) anf J together or separately
51
 
 
52
 
  begin("Starting Newton solve.");
53
 
 
54
 
  // Compute F(u) and J
55
 
  set("output destination", "silent");
56
 
  nonlinear_problem.form(A, b, x);
57
 
 
58
 
  uint krylov_iterations = 0;
59
 
  newton_iteration = 0;
60
 
  bool newton_converged = false;
61
 
 
62
 
  // Start iterations
63
 
  while( !newton_converged && newton_iteration < maxit )
64
 
  {
65
 
 
66
 
      set("output destination", "silent");
67
 
      // Perform linear solve and update total number of Krylov iterations
68
 
      krylov_iterations += solver->solve(A, dx, b);
69
 
      set("output destination", "terminal");
70
 
 
71
 
      // Compute initial residual
72
 
      if(newton_iteration == 0)
73
 
        newton_converged = converged(b, dx, nonlinear_problem);
74
 
 
75
 
      // Update solution
76
 
      x += dx;
77
 
 
78
 
      ++newton_iteration;
79
 
 
80
 
      set("output destination", "silent");
81
 
      //FIXME: this step is not needed if residual is based on dx and this has converged.
82
 
      // Compute F(u) and J
83
 
      nonlinear_problem.form(A, b, x);
84
 
 
85
 
      set("output destination", "terminal");
86
 
 
87
 
      // Test for convergence 
88
 
      newton_converged = converged(b, dx, nonlinear_problem);
89
 
  }
90
 
 
91
 
  if(newton_converged)
92
 
    message("Newton solver finished in %d iterations and %d linear solver iterations.", 
93
 
        newton_iteration, krylov_iterations);
94
 
  else
95
 
    warning("Newton solver did not converge.");
96
 
 
97
 
  end();
98
 
 
99
 
  return newton_iteration;
100
 
} 
101
 
//-----------------------------------------------------------------------------
102
 
dolfin::uint NewtonSolver::getIteration() const
103
 
{
104
 
  return newton_iteration; 
105
 
}
106
 
//-----------------------------------------------------------------------------
107
 
bool NewtonSolver::converged(const Vector& b, const Vector& dx, 
108
 
      const NonlinearProblem& nonlinear_problem)
109
 
{
110
 
  const std::string convergence_criterion = get("Newton convergence criterion");
111
 
  const real rtol   = get("Newton relative tolerance");
112
 
  const real atol   = get("Newton absolute tolerance");
113
 
  const bool report = get("Newton report");
114
 
 
115
 
  real residual = 1.0;
116
 
 
117
 
  // Compute resdiual
118
 
  if(convergence_criterion == "residual")
119
 
    residual = b.norm();
120
 
  else if (convergence_criterion == "incremental")
121
 
    residual = dx.norm();
122
 
  else
123
 
    error("Unknown Newton convergence criterion");
124
 
 
125
 
  // If this is the first iteration step, set initial residual
126
 
  if(newton_iteration == 0)
127
 
    residual0 = residual;
128
 
 
129
 
  // Relative residual
130
 
  real relative_residual = residual/residual0;
131
 
 
132
 
  // Output iteration number and residual
133
 
  //FIXME: allow precision to be set for dolfin::cout<<
134
 
  std::cout.precision(3);
135
 
  if(report && newton_iteration >0) 
136
 
    std::cout<< "  Iteration = " << newton_iteration  << ", Absolute, relative residual (" 
137
 
    << convergence_criterion  << " criterion) = " << std::scientific << residual 
138
 
    << ", "<< std::scientific << relative_residual << std::endl;
139
 
 
140
 
  // Return true of convergence criterion is met
141
 
  if(relative_residual < rtol || residual < atol)
142
 
    return true;
143
 
 
144
 
  // Otherwise return false
145
 
  return false;
146
 
}
147
 
//-----------------------------------------------------------------------------
148