← Back to branch summary

~ubuntu-branches/ubuntu/maverick/dolfin/maverick

« back to all changes in this revision

Viewing changes to dolfin/ode/MultiAdaptiveNewtonSolver.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Johannes Ring
  • Date: 2009-10-12 14:13:18 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20091012141318-hkbxl0sq555vqv7d
Tags: 0.9.4-1
* New upstream release. This version cleans up the design of the
  function class by adding a new abstraction for user-defined
  functions called Expression. A number of minor bugfixes and
  improvements have also been made.
* debian/watch: Update for new flat directory structure.
* Update debian/copyright.
* debian/rules: Use explicit paths to PETSc 3.0.0 and SLEPc 3.0.0.

Show diffs side-by-side

added added

removed removed

Lines of Context:
dolfin/ode/MultiAdaptiveNewtonSolver.cpp
1
 
// Copyright (C) 2005-2008 Anders Logg.
 
1
// Copyright (C) 2005-2009 Anders Logg.
2
2
// Licensed under the GNU LGPL Version 2.1.
3
3
//
4
4
// First added:  2005-01-27
5
 
// Last changed: 2008-04-22
 
5
// Last changed: 2009-09-08
6
6
 
7
7
#include <dolfin/common/constants.h>
8
8
#include <dolfin/log/dolfin_log.h>
9
9
#include <dolfin/math/dolfin_math.h>
10
 
#include <dolfin/parameter/parameters.h>
11
10
#include <dolfin/la/uBLASSparseMatrix.h>
12
11
#include "Alloc.h"
13
12
#include "ODE.h"
24
23
(MultiAdaptiveTimeSlab& timeslab)
25
24
  : TimeSlabSolver(timeslab), ts(timeslab), A(0),
26
25
    mpc(timeslab, method), solver(mpc), f(0), u(0), num_elements(0), num_elements_mono(0),
27
 
    updated_jacobian(ode.get("ODE updated jacobian"))
 
26
    updated_jacobian(ode.parameters["updated_jacobian"])
28
27
{
29
28
  // Initialize local arrays
30
29
  f = new real[method.qsize()];
31
30
  u = new real[method.nsize()];
32
 
  
33
 
  // Don't report number of GMRES iteration if not asked to
34
 
  solver.set("Krylov report", monitor);
35
 
  solver.set("Krylov absolute tolerance", 0.01);
36
 
  solver.set("Krylov relative tolerance", 0.01 * to_double(tol));
37
 
  
 
31
 
 
32
  // Set parameters for Krylov solver
 
33
  solver.parameters["report"] = monitor;
 
34
  solver.parameters["absolute_tolerance"] = 0.01;
 
35
  solver.parameters["relative_tolerance"] = 0.01 * to_double(tol);
 
36
 
38
37
  // Initialize Jacobian
39
38
  if ( updated_jacobian )
40
39
    A = new UpdatedMultiAdaptiveJacobian(*this, timeslab);
48
47
  if ( num_elements > 0 )
49
48
  {
50
49
    const real alpha = num_elements_mono / static_cast<real>(num_elements);
51
 
    message("Multi-adaptive efficiency index: %.3f", to_double(alpha));
 
50
    info("Multi-adaptive efficiency index: %.3f", to_double(alpha));
52
51
  }
53
 
  
 
52
 
54
53
  // Delete local arrays
55
54
  if ( f ) delete [] f;
56
55
  if ( u ) delete [] u;
77
76
  // Initialize right-hand side
78
77
  b.resize(nj);
79
78
  b.zero();
80
 
  
 
79
 
81
80
  // Recompute Jacobian on each time slab
82
81
  A->init();
83
82
 
84
83
  //debug();
85
 
  //A->disp(true, 10);
 
84
  //info(A);
86
85
}
87
86
//-----------------------------------------------------------------------------
88
 
real MultiAdaptiveNewtonSolver::iteration(real tol, uint iter,
89
 
                                          real d0, real d1)
 
87
real MultiAdaptiveNewtonSolver::iteration(const real& tol, uint iter, const real& d0, const real& d1)
90
88
{
91
89
  // Evaluate b = -F(x) at current x
92
90
  Feval(b);
93
 
 
 
91
 
94
92
  // FIXME: Scaling needed for PETSc Krylov solver, but maybe not for uBLAS?
95
93
 
96
94
  // Save norm of old solution
97
95
  xnorm = 0.0;
98
96
  for (uint j = 0; j < ts.nj; j++)
99
 
    xnorm = max(xnorm, abs(ts.jx[j]));
100
 
  
 
97
    xnorm = real_max(xnorm, real_abs(ts.jx[j]));
 
98
 
101
99
  // Solve linear system
102
 
  const double r = b.norm(linf) + to_double( real_epsilon() );
 
100
  const double r = b.norm("linf") + to_double( real_epsilon() );
103
101
  b /= r;
104
102
  num_local_iterations += solver.solve(*A, dx, b);
105
103
  dx *= r;
112
110
  real max_increment = 0.0;
113
111
  for (uint j = 0; j < ts.nj; j++)
114
112
  {
115
 
    const real increment = fabs(dx[j]);
116
 
    if ( increment > max_increment )
 
113
    const real increment = real_abs(dx[j]);
 
114
    if (increment > max_increment)
117
115
      max_increment = increment;
118
116
  }
119
117
 
141
139
  for (uint e = 0; e < ts.ne; e++)
142
140
  {
143
141
    // Cover all elements in current sub slab
144
 
    s = ts.coverNext(s, e);
 
142
    s = ts.cover_next(s, e);
145
143
 
146
144
    // Get element data
147
145
    const uint i = ts.ei[e];
155
153
 
156
154
    // Evaluate right-hand side at quadrature points of element
157
155
    if ( method.type() == Method::cG )
158
 
      ts.cGfeval(f, s, e, i, a, b, k);
 
156
      ts.cg_feval(f, s, e, i, a, b, k);
159
157
    else
160
 
      ts.dGfeval(f, s, e, i, a, b, k);  
161
 
    //cout << "f = "; Alloc::disp(f, method.qsize());
 
158
      ts.dg_feval(f, s, e, i, a, b, k);
162
159
 
163
160
    // Update values on element using fixed-point iteration
164
161
    method.update(x0, f, k, u);
165
 
    
 
162
 
166
163
    // Subtract current values
167
164
    for (uint n = 0; n < method.nsize(); n++)
168
165
      F[j + n] = to_double(u[j] - ts.jx[j + n]);
183
180
  for (uint j = 0; j < n; j++)
184
181
  {
185
182
    const real xj = ts.jx[j];
186
 
    real dx = max(DOLFIN_SQRT_EPS, DOLFIN_SQRT_EPS * abs(xj));
187
 
                  
 
183
    real dx = real_max(DOLFIN_SQRT_EPS, DOLFIN_SQRT_EPS * real_abs(xj));
 
184
 
188
185
    ts.jx[j] -= 0.5*dx;
189
186
    Feval(F1);
190
187
 
195
192
 
196
193
    for (uint i = 0; i < n; i++)
197
194
    {
198
 
      real dFdx = (F1[i] - F2[i]) / dx;
199
 
      if ( abs(dFdx) > real_epsilon() )
200
 
        _B(i, j) = to_double(dFdx);
 
195
      real df_dx = (F1[i] - F2[i]) / dx;
 
196
      if (real_abs(df_dx) > real_epsilon())
 
197
        _B(i, j) = to_double(df_dx);
201
198
    }
202
199
  }
203
200
 
204
 
  B.disp();
 
201
  info(B);
205
202
}
206
203
//-----------------------------------------------------------------------------