← Back to branch summary

~njansson/dolfin/hpc

« back to all changes in this revision

Viewing changes to src/kernel/ode/MultiAdaptivity.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/MultiAdaptivity.cpp
1
 
// Copyright (C) 2005-2006 Anders Logg.
2
 
// Licensed under the GNU LGPL Version 2.1.
3
 
//
4
 
// First added:  2005-01-29
5
 
// Last changed: 2006-04-20
6
 
 
7
 
#include <cmath>
8
 
#include <dolfin/parameters.h>
9
 
#include <dolfin/Array.h>
10
 
#include <dolfin/ODE.h>
11
 
#include <dolfin/Method.h>
12
 
#include <dolfin/MultiAdaptiveTimeSlab.h>
13
 
#include <dolfin/MultiAdaptivity.h>
14
 
 
15
 
using namespace dolfin;
16
 
 
17
 
//-----------------------------------------------------------------------------
18
 
MultiAdaptivity::MultiAdaptivity(const ODE& ode, const Method& method)
19
 
  : Adaptivity(ode, method), 
20
 
    timesteps(0), residuals(0), ktmp(0), f(0), rmax(0), emax(0)
21
 
{
22
 
  // Initialize time steps and residuals
23
 
  timesteps = new real[ode.size()];
24
 
  residuals = new real[ode.size()];
25
 
  ktmp = new real[ode.size()];
26
 
 
27
 
  // Initialize local array for quadrature
28
 
  f = new real[method.qsize()];
29
 
  for (unsigned int i = 0; i < method.qsize(); i++)
30
 
    f[i] = 0.0;
31
 
 
32
 
  // Set initial time steps
33
 
  real k0 = dolfin_get("ODE initial time step");
34
 
  if ( kfixed )
35
 
  {
36
 
    for (uint i = 0; i < ode.size(); i++)
37
 
    {
38
 
      timesteps[i] = ode.timestep(0.0, i, k0);
39
 
    }
40
 
  }
41
 
  else
42
 
  {
43
 
    real k = k0;
44
 
    if ( k > _kmax )
45
 
    {
46
 
      k = _kmax;
47
 
      warning("Initial time step larger than maximum time step, using k = %.3e.", k);
48
 
    }
49
 
    for (uint i = 0; i < ode.size(); i++)
50
 
      timesteps[i] = k;
51
 
  }
52
 
 
53
 
  // Initialize arrays
54
 
  for (uint i = 0; i < ode.size(); i++)
55
 
  {
56
 
    residuals[i] = 0.0;
57
 
    ktmp[i] = 0.0;
58
 
  }
59
 
}
60
 
//-----------------------------------------------------------------------------
61
 
MultiAdaptivity::~MultiAdaptivity()
62
 
{
63
 
  if ( timesteps ) delete [] timesteps;
64
 
  if ( residuals ) delete [] residuals;
65
 
  if ( ktmp ) delete [] ktmp;
66
 
}
67
 
//-----------------------------------------------------------------------------
68
 
real MultiAdaptivity::timestep(uint i) const
69
 
{
70
 
  return timesteps[i];
71
 
}
72
 
//-----------------------------------------------------------------------------
73
 
real MultiAdaptivity::residual(uint i) const
74
 
{
75
 
  return residuals[i];
76
 
}
77
 
//-----------------------------------------------------------------------------
78
 
void MultiAdaptivity::update(MultiAdaptiveTimeSlab& ts, real t, bool first)
79
 
{
80
 
  _accept = false;
81
 
 
82
 
  // Check if time step is fixed
83
 
  if ( kfixed )
84
 
  {
85
 
    for (uint i = 0; i < ts.N; i++)
86
 
      timesteps[i] = ode.timestep(t, i, timesteps[i]);
87
 
   
88
 
    _accept = true;
89
 
    return;
90
 
  }
91
 
 
92
 
  // Compute maximum residuals for all components in time slab
93
 
  computeResiduals(ts);
94
 
  
95
 
  // Accept if error is small enough
96
 
  if ( emax <= tol )
97
 
    _accept = true;
98
 
 
99
 
  // Update time steps for all components
100
 
  for (uint i = 0; i < ode.size(); i++)
101
 
  {
102
 
    // Previous time step
103
 
    const real k0 = timesteps[i];
104
 
    
105
 
    // Include dynamic safety factor
106
 
    real used_tol = safety*tol;
107
 
    
108
 
    // Compute new time step
109
 
    real k = method.timestep(residuals[i], used_tol, k0, _kmax);
110
 
 
111
 
    // Apply time step regulation
112
 
    k = Controller::updateHarmonic(k, timesteps[i], _kmax);
113
 
    
114
 
    // Make sure to decrease the time step if not accepted
115
 
    if ( !_accept )
116
 
    {
117
 
      k = std::min(k, 0.9*k0);
118
 
    }
119
 
    
120
 
    // Save time step for component
121
 
    timesteps[i] = k;
122
 
  }
123
 
 
124
 
  // Propagate time steps according to dependencies
125
 
  propagateDependencies();
126
 
}
127
 
//-----------------------------------------------------------------------------
128
 
void MultiAdaptivity::computeResiduals(MultiAdaptiveTimeSlab& ts)
129
 
{
130
 
  // Reset dof
131
 
  uint j = 0;
132
 
 
133
 
  // Reset elast
134
 
  for (uint i = 0; i < ts.N; i++)
135
 
    ts.elast[i] = -1;
136
 
 
137
 
  // Reset residuals
138
 
  for (uint i = 0; i < ts.N; i++)
139
 
    residuals[i] = 0.0;
140
 
  
141
 
  // Reset maximum local residual and error
142
 
  rmax = 0.0;
143
 
  emax = 0.0;
144
 
 
145
 
  // Iterate over all sub slabs
146
 
  uint e0 = 0;
147
 
  uint e1 = 0;
148
 
  for (uint s = 0; s < ts.ns; s++)
149
 
  {
150
 
    // Cover all elements in current sub slab
151
 
    e1 = ts.coverSlab(s, e0);
152
 
    
153
 
    // Get data for sub slab
154
 
    const real a = ts.sa[s];
155
 
    const real b = ts.sb[s];
156
 
    const real k = b - a;
157
 
 
158
 
    // Iterate over all elements in current sub slab
159
 
    for (uint e = e0; e < e1; e++)
160
 
    {
161
 
      // Get element data
162
 
      const uint i = ts.ei[e];
163
 
 
164
 
      // Get initial value for element
165
 
      const int ep = ts.ee[e];
166
 
      const real x0 = ( ep != -1 ? ts.jx[ep*method.nsize() + method.nsize() - 1] : ts.u0(i) );
167
 
      
168
 
      // Evaluate right-hand side at quadrature points of element
169
 
      if ( method.type() == Method::cG )
170
 
        ts.cGfeval(f, s, e, i, a, b, k);
171
 
      else
172
 
        ts.dGfeval(f, s, e, i, a, b, k);
173
 
 
174
 
      // Update maximum residual for component
175
 
      const real r = method.residual(x0, ts.jx + j, f[method.nsize()], k);
176
 
      residuals[i] = std::max(residuals[i], fabs(r));
177
 
 
178
 
      // Update maximum residual and error
179
 
      rmax = std::max(rmax, r);
180
 
      emax = std::max(emax, method.error(k, r));
181
 
 
182
 
      // Update dof
183
 
      j += method.nsize();
184
 
    }
185
 
 
186
 
    // Step to next sub slab
187
 
    e0 = e1;
188
 
  }
189
 
}
190
 
//-----------------------------------------------------------------------------
191
 
void MultiAdaptivity::propagateDependencies()
192
 
{
193
 
  // This is a poor man's dual weighting function. For each component,
194
 
  // we look at all other components that the component in question
195
 
  // depends on and tell these other components to reduce their time
196
 
  // steps to the same level. A similar effect would be accomplished
197
 
  // by weighting with the proper weight from the dual solution, but
198
 
  // until we (re-)implement the solution of the dual, this seems to
199
 
  // be a good solution.
200
 
 
201
 
  // Don't propagate dependencies if not sparse
202
 
  if ( !ode.dependencies.sparse() )
203
 
    return;
204
 
 
205
 
  // Copy time steps
206
 
  for (uint i = 0; i < ode.size(); i++)
207
 
    ktmp[i] = timesteps[i];
208
 
 
209
 
  // Iterate over components
210
 
  for (uint i = 0; i < ode.size(); i++)
211
 
  {
212
 
    // Get time step for current component
213
 
    const real k = ktmp[i];
214
 
 
215
 
    // Propagate time step to dependencies
216
 
    const Array<uint>& deps = ode.dependencies[i];
217
 
    for (uint pos = 0; pos < deps.size(); pos++)
218
 
      timesteps[deps[pos]] = std::min(timesteps[deps[pos]], k);
219
 
  }
220
 
}
221
 
//-----------------------------------------------------------------------------