← Back to branch summary

~njansson/dolfin/hpc

« back to all changes in this revision

Viewing changes to src/kernel/ode/TimeSlabSolver.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/ode/TimeSlabSolver.cpp
1
 
// Copyright (C) 2005-2006 Anders Logg.
2
 
// Licensed under the GNU LGPL Version 2.1.
3
 
//
4
 
// First added:  2005-01-05
5
 
// Last changed: 2006-04-20
6
 
 
7
 
#include <cmath>
8
 
#include <dolfin/parameters.h>
9
 
#include <dolfin/TimeSlab.h>
10
 
#include <dolfin/TimeSlabSolver.h>
11
 
 
12
 
using namespace dolfin;
13
 
 
14
 
//-----------------------------------------------------------------------------
15
 
TimeSlabSolver::TimeSlabSolver(TimeSlab& timeslab)
16
 
  : ode(timeslab.ode), method(*timeslab.method), tol(0.0), maxiter(0),
17
 
    monitor(dolfin_get("ODE monitor convergence")),
18
 
    num_timeslabs(0), num_global_iterations(0), num_local_iterations(0)
19
 
{
20
 
  // Choose tolerance
21
 
  chooseTolerance();
22
 
 
23
 
  // Get maximum number of iterations
24
 
  maxiter = dolfin_get("ODE maximum iterations");
25
 
}
26
 
//-----------------------------------------------------------------------------
27
 
TimeSlabSolver::~TimeSlabSolver()
28
 
{
29
 
  if ( num_timeslabs > 0 )
30
 
  {
31
 
    const real n = static_cast<real>(num_timeslabs);
32
 
    const real global_average = static_cast<real>(num_global_iterations) / n;
33
 
    const real local_average = static_cast<real>(num_local_iterations) / 
34
 
      static_cast<real>(num_global_iterations);
35
 
    message("Average number of global iterations per step: %.3f",
36
 
                global_average);
37
 
    message("Average number of local iterations per global iteration: %.3f",
38
 
                local_average);
39
 
  }
40
 
 
41
 
  message("Total number of (macro) time steps: %d", num_timeslabs);
42
 
}
43
 
//-----------------------------------------------------------------------------
44
 
bool TimeSlabSolver::solve()
45
 
{
46
 
  for (uint attempt = 0; attempt < maxiter; attempt++)
47
 
  {
48
 
    // Try to solve system
49
 
    if ( solve(attempt) )
50
 
      return true;
51
 
    
52
 
    // Check if we should try again
53
 
    if ( !retry() )
54
 
      return false;
55
 
  }
56
 
 
57
 
  return false;
58
 
}
59
 
//-----------------------------------------------------------------------------
60
 
bool TimeSlabSolver::solve(uint attempt)
61
 
{
62
 
  start();
63
 
 
64
 
  real d0 = 0.0;
65
 
  real d1 = 0.0;
66
 
  for (uint iter = 0; iter < maxiter; iter++)
67
 
  {
68
 
    // Do one iteration
69
 
    real d2 = iteration(tol, iter, d0, d1);
70
 
    if ( monitor )
71
 
      message("--- iter = %d: increment = %.3e", iter, d2);
72
 
    
73
 
    // Check convergenge
74
 
    if ( d2 < tol )
75
 
    {
76
 
      end();
77
 
      num_timeslabs += 1;
78
 
      num_global_iterations += iter + 1;
79
 
      if ( monitor )
80
 
        message("Time slab system of size %d converged in %d iterations.", size(), iter + 1);
81
 
      return true;
82
 
    }
83
 
 
84
 
    // Check divergence
85
 
    // FIXME: implement better check and make this a parameter
86
 
    if ( (iter > 0 && d2 > 1000.0 * d1) || !std::isnormal(d2) )
87
 
    {
88
 
      warning("Time slab system seems to be diverging.");
89
 
      return false;
90
 
    }
91
 
    
92
 
    d0 = d1;
93
 
    d1 = d2;
94
 
  }
95
 
 
96
 
  warning("Time slab system did not converge.");
97
 
  return false;
98
 
}
99
 
//-----------------------------------------------------------------------------
100
 
bool TimeSlabSolver::retry()
101
 
{
102
 
  // By default, we don't know how to make a new attempt
103
 
  return false;
104
 
}
105
 
//-----------------------------------------------------------------------------
106
 
void TimeSlabSolver::start()
107
 
{
108
 
  // Do nothing
109
 
}
110
 
//-----------------------------------------------------------------------------
111
 
void TimeSlabSolver::end()
112
 
{
113
 
  // Do nothing
114
 
}
115
 
//-----------------------------------------------------------------------------
116
 
void TimeSlabSolver::chooseTolerance()
117
 
{
118
 
  const real TOL   = dolfin_get("ODE tolerance");
119
 
  const real alpha = dolfin_get("ODE discrete tolerance factor");
120
 
 
121
 
  tol = dolfin_get("ODE discrete tolerance");
122
 
  if ( !dolfin_get("ODE fixed time step") )
123
 
    tol = std::min(tol, alpha*TOL);
124
 
  cout << "Using discrete tolerance tol = " << tol << "." << endl;
125
 
}
126
 
//-----------------------------------------------------------------------------