← Back to branch summary

~njansson/dolfin/hpc

« back to all changes in this revision

Viewing changes to dolfin/ode/ReducedModel.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:
dolfin/ode/ReducedModel.cpp
 
1
// Copyright (C) 2004-2005 Anders Logg.
 
2
// Licensed under the GNU LGPL Version 2.1.
 
3
//
 
4
// First added:  2004-04-04
 
5
// Last changed: 2005-12-19
 
6
 
 
7
/*
 
8
 
 
9
// FIXME: BROKEN
 
10
 
 
11
#include <cmath>
 
12
#include <dolfin/parameter/parameters.h>
 
13
#include "ReducedModel.h"
 
14
 
 
15
using namespace dolfin;
 
16
 
 
17
//-----------------------------------------------------------------------------
 
18
ReducedModel::ReducedModel(ODE& ode) 
 
19
  : ODE(ode.size(), ode.endtime()), ode(ode), g(ode.size()), reduced(false)
 
20
{
 
21
  warning("Automatic modeling is EXPERIMENTAL.");
 
22
 
 
23
  tau     = get("ODE average length");
 
24
  samples = get("ODE average samples");
 
25
  tol     = get("ODE average tolerance");
 
26
 
 
27
  // Copy the sparsity
 
28
  //sparsity = ode.sparsity;
 
29
 
 
30
  // Adjust the maximum allowed time step to the initial time step
 
31
  real kmax = get("ODE initial time step");
 
32
  set("ODE maximum time step", kmax);
 
33
}
 
34
//-----------------------------------------------------------------------------
 
35
ReducedModel::~ReducedModel()
 
36
{
 
37
  // Do nothing
 
38
}
 
39
//-----------------------------------------------------------------------------
 
40
real ReducedModel::f(const Vector& u, real t, unsigned int i)
 
41
{
 
42
  if ( g[i].active() )
 
43
    return ode.f(u, t, i) + g[i]();
 
44
  else
 
45
    return 0.0;
 
46
 
 
47
  return 0.0;
 
48
}
 
49
//-----------------------------------------------------------------------------
 
50
real ReducedModel::u0(unsigned int i)
 
51
{
 
52
  return ode.u0(i);
 
53
}
 
54
//-----------------------------------------------------------------------------
 
55
void ReducedModel::update(RHS& f, Function& u, real t)
 
56
{
 
57
  // Update for the given model if used
 
58
  ode.update(f, u, t);
 
59
 
 
60
  // Create the reduced model only one time
 
61
  if ( reduced )
 
62
    return;
 
63
 
 
64
  // Check if we have reached the beyond twice the average length
 
65
  if ( t < (2.0*tau) )
 
66
    return;
 
67
 
 
68
  // Write a message
 
69
  message("Creating reduced model at t = %.1e.", 2*tau);
 
70
 
 
71
  // Compute averages of u and f
 
72
  computeAverages(f, u, fbar, ubar);
 
73
 
 
74
  // Compute reduced model
 
75
  for (unsigned int i = 0; i < ode.size(); i++)
 
76
    g[i].computeModel(ubar, fbar, i, tau, ode);
 
77
}
 
78
//-----------------------------------------------------------------------------
 
79
void ReducedModel::update(Solution& u, Adaptivity& adaptivity, real t)
 
80
{
 
81
  // Update for the given model if used
 
82
  ode.update(u, adaptivity, t);
 
83
 
 
84
  // Check if we have reached the beyond twice the average length
 
85
  if ( t < (2.0*tau) )
 
86
    return;
 
87
 
 
88
  // Create the reduced model only one time
 
89
  if ( reduced )
 
90
    return;  
 
91
 
 
92
  // FIXME: Choose better maximum time step
 
93
  // Adjust (increase) maximum allowed time step
 
94
  adaptivity.adjustMaximumTimeStep(0.1);
 
95
 
 
96
  // Adjust end values for inactive components
 
97
  for (unsigned int i = 0; i < ode.size(); i++)
 
98
    if ( !g[i].active() )
 
99
      u.setlast(i, ubar(i));
 
100
 
 
101
  // Clear averages
 
102
  ubar.clear();
 
103
  fbar.clear();
 
104
 
 
105
  // Remember that we have created the model
 
106
  reduced = true;
 
107
}
 
108
//-----------------------------------------------------------------------------
 
109
void ReducedModel::save(Sample& sample)
 
110
{
 
111
  ode.save(sample);
 
112
}
 
113
//-----------------------------------------------------------------------------
 
114
void ReducedModel::computeAverages(RHS& f, Function& u,
 
115
                                   Vector& fbar, Vector& ubar)
 
116
{
 
117
error("This function needs to be updated to the new format.");
 
118
 
 
119
  // Sample length
 
120
  real k = tau / static_cast<real>(samples);
 
121
 
 
122
  // Weight for each sample
 
123
  real w = 1.0 / static_cast<real>(samples);
 
124
 
 
125
  // Initialize averages
 
126
  ubar.init(ode.size()); // Average on first half (and then both...)
 
127
  fbar.init(ode.size()); // Average on both intervals
 
128
  ubar = 0.0;
 
129
  fbar = 0.0;
 
130
 
 
131
  // Minimum and maximum values
 
132
  Vector umin(ode.size());
 
133
  Vector umax(ode.size());
 
134
  umin = 0.0;
 
135
  umax = 0.0;
 
136
 
 
137
  // Additional average on second interval for u
 
138
  Vector ubar2(ode.size());
 
139
  ubar2 = 0.0;
 
140
 
 
141
  Progress p("Computing averages", 2*samples);
 
142
 
 
143
  // Compute average values on first half of the averaging interval
 
144
  for (unsigned int j = 0; j < samples; j++)
 
145
  {
 
146
    for (unsigned int i = 0; i < ode.size(); i++)
 
147
    {
 
148
      // Sample time
 
149
      real t = 0.5*k + static_cast<real>(j)*k;
 
150
      
 
151
      // Compute u and f
 
152
      real uu = u(i, t);
 
153
      real ff = f(i, t);
 
154
 
 
155
      // Compute average values of u and f
 
156
      ubar(i) += w * uu;
 
157
      fbar(i) += 0.5 * w * ff;
 
158
      
 
159
      // Compute minimum and maximum values of u and f
 
160
      if ( j == 0 )
 
161
      {
 
162
        umin(i) = uu;
 
163
        umax(i) = uu;
 
164
      }
 
165
      else
 
166
      {
 
167
        umin(i) = std::min(umin(i), uu);
 
168
        umax(i) = std::max(umax(i), uu);
 
169
      }
 
170
    }
 
171
    
 
172
    ++p;
 
173
  }
 
174
    
 
175
  // Compute average values on second half of the averaging interval
 
176
  for (unsigned int j = 0; j < samples; j++)
 
177
  {
 
178
    for (unsigned int i = 0; i < ode.size(); i++)
 
179
    {
 
180
      // Sample time
 
181
      real t = tau + 0.5*k + static_cast<real>(j)*k;
 
182
      
 
183
      // Compute u and f
 
184
      real uu = u(i, t);
 
185
      real ff = f(i, t);
 
186
      
 
187
      // Compute average values of u and f
 
188
      ubar2(i) += w * uu;
 
189
      fbar(i)  += 0.5 * w * ff;
 
190
      
 
191
      // Compute minimum and maximum values of u and f
 
192
      umin(i) = std::min(umin(i), uu);
 
193
      umax(i) = std::max(umax(i), uu);
 
194
    }
 
195
 
 
196
    ++p;
 
197
  }
 
198
 
 
199
  // Compute relative changes in u
 
200
  for (unsigned int i = 0; i < ode.size(); i++)
 
201
  {
 
202
    real du = fabs(ubar(i) - ubar2(i)) / (umax(i) - umin(i) + DOLFIN_EPS);
 
203
 
 
204
    cout << "Component " << i << ":" << endl;
 
205
    cout << "  ubar1 = " << ubar(i) << endl;
 
206
    cout << "  ubar2 = " << ubar2(i) << endl;
 
207
    cout << "  fbar  = " << fbar(i) << endl;
 
208
    cout << "  du/u  = " << du << endl;
 
209
 
 
210
    // Inactivate components whose average is constant
 
211
    if ( du < tol )
 
212
    {
 
213
      cout << "Inactivating component " << i << endl;
 
214
      g[i].inactivate();
 
215
    }
 
216
 
 
217
    // Save average
 
218
    ubar(i) = 0.5*(ubar(i) + ubar2(i));
 
219
  }
 
220
 
 
221
}
 
222
//-----------------------------------------------------------------------------
 
223
// ReducedModel::Model
 
224
//-----------------------------------------------------------------------------
 
225
ReducedModel::Model::Model() : g(0), _active(true)
 
226
{
 
227
  // Do nothing
 
228
}
 
229
//-----------------------------------------------------------------------------
 
230
ReducedModel::Model::~Model()
 
231
{
 
232
  // Do nothing
 
233
}
 
234
//-----------------------------------------------------------------------------
 
235
real ReducedModel::Model::operator() () const
 
236
{
 
237
  return g;
 
238
}
 
239
//-----------------------------------------------------------------------------
 
240
bool ReducedModel::Model::active() const
 
241
{
 
242
  return _active;
 
243
}
 
244
//-----------------------------------------------------------------------------
 
245
void ReducedModel::Model::inactivate()
 
246
{
 
247
  _active = false;
 
248
}
 
249
//-----------------------------------------------------------------------------
 
250
void ReducedModel::Model::computeModel(Vector& ubar, Vector& fbar,
 
251
                                       unsigned int i, real tau, ODE& ode)
 
252
{
 
253
  // FIXME: BROKEN
 
254
  //g = fbar(i) - ode.f(ubar, tau, i);
 
255
 
 
256
  cout << "Modeling term for component " << i << ": " << g << endl;
 
257
}
 
258
//-----------------------------------------------------------------------------
 
259
*/