1
// Copyright (C) 2005-2006 Anders Logg.
2
// Licensed under the GNU LGPL Version 2.1.
4
// First added: 2005-01-29
5
// Last changed: 2006-04-20
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>
15
using namespace dolfin;
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)
22
// Initialize time steps and residuals
23
timesteps = new real[ode.size()];
24
residuals = new real[ode.size()];
25
ktmp = new real[ode.size()];
27
// Initialize local array for quadrature
28
f = new real[method.qsize()];
29
for (unsigned int i = 0; i < method.qsize(); i++)
32
// Set initial time steps
33
real k0 = dolfin_get("ODE initial time step");
36
for (uint i = 0; i < ode.size(); i++)
38
timesteps[i] = ode.timestep(0.0, i, k0);
47
warning("Initial time step larger than maximum time step, using k = %.3e.", k);
49
for (uint i = 0; i < ode.size(); i++)
54
for (uint i = 0; i < ode.size(); i++)
60
//-----------------------------------------------------------------------------
61
MultiAdaptivity::~MultiAdaptivity()
63
if ( timesteps ) delete [] timesteps;
64
if ( residuals ) delete [] residuals;
65
if ( ktmp ) delete [] ktmp;
67
//-----------------------------------------------------------------------------
68
real MultiAdaptivity::timestep(uint i) const
72
//-----------------------------------------------------------------------------
73
real MultiAdaptivity::residual(uint i) const
77
//-----------------------------------------------------------------------------
78
void MultiAdaptivity::update(MultiAdaptiveTimeSlab& ts, real t, bool first)
82
// Check if time step is fixed
85
for (uint i = 0; i < ts.N; i++)
86
timesteps[i] = ode.timestep(t, i, timesteps[i]);
92
// Compute maximum residuals for all components in time slab
95
// Accept if error is small enough
99
// Update time steps for all components
100
for (uint i = 0; i < ode.size(); i++)
102
// Previous time step
103
const real k0 = timesteps[i];
105
// Include dynamic safety factor
106
real used_tol = safety*tol;
108
// Compute new time step
109
real k = method.timestep(residuals[i], used_tol, k0, _kmax);
111
// Apply time step regulation
112
k = Controller::updateHarmonic(k, timesteps[i], _kmax);
114
// Make sure to decrease the time step if not accepted
117
k = std::min(k, 0.9*k0);
120
// Save time step for component
124
// Propagate time steps according to dependencies
125
propagateDependencies();
127
//-----------------------------------------------------------------------------
128
void MultiAdaptivity::computeResiduals(MultiAdaptiveTimeSlab& ts)
134
for (uint i = 0; i < ts.N; i++)
138
for (uint i = 0; i < ts.N; i++)
141
// Reset maximum local residual and error
145
// Iterate over all sub slabs
148
for (uint s = 0; s < ts.ns; s++)
150
// Cover all elements in current sub slab
151
e1 = ts.coverSlab(s, e0);
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;
158
// Iterate over all elements in current sub slab
159
for (uint e = e0; e < e1; e++)
162
const uint i = ts.ei[e];
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) );
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);
172
ts.dGfeval(f, s, e, i, a, b, k);
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));
178
// Update maximum residual and error
179
rmax = std::max(rmax, r);
180
emax = std::max(emax, method.error(k, r));
186
// Step to next sub slab
190
//-----------------------------------------------------------------------------
191
void MultiAdaptivity::propagateDependencies()
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.
201
// Don't propagate dependencies if not sparse
202
if ( !ode.dependencies.sparse() )
206
for (uint i = 0; i < ode.size(); i++)
207
ktmp[i] = timesteps[i];
209
// Iterate over components
210
for (uint i = 0; i < ode.size(); i++)
212
// Get time step for current component
213
const real k = ktmp[i];
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);
221
//-----------------------------------------------------------------------------