← Back to branch summary

~njansson/dolfin/hpc

« back to all changes in this revision

Viewing changes to src/kernel/ode/MultiAdaptiveJacobian.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/MultiAdaptiveJacobian.cpp
1
 
// Copyright (C) 2005-2006 Anders Logg.
2
 
// Licensed under the GNU LGPL Version 2.1.
3
 
//
4
 
// First added:  2005-01-27
5
 
// Last changed: 2006-08-08
6
 
 
7
 
#include <dolfin/dolfin_math.h>
8
 
#include <dolfin/uBlasVector.h>
9
 
#include <dolfin/ODE.h>
10
 
#include <dolfin/Method.h>
11
 
#include <dolfin/MultiAdaptiveTimeSlab.h>
12
 
#include <dolfin/MultiAdaptiveNewtonSolver.h>
13
 
#include <dolfin/MultiAdaptiveJacobian.h>
14
 
 
15
 
using namespace dolfin;
16
 
 
17
 
//-----------------------------------------------------------------------------
18
 
MultiAdaptiveJacobian::MultiAdaptiveJacobian(MultiAdaptiveNewtonSolver& newton,
19
 
                                             MultiAdaptiveTimeSlab& timeslab)
20
 
  : TimeSlabJacobian(timeslab), newton(newton), ts(timeslab),
21
 
    Jvalues(0), Jindices(0), Jlookup(0)
22
 
{
23
 
  // Allocate Jacobian row indices
24
 
  Jindices = new uint[ode.size()];
25
 
  
26
 
  // Compute start of each row
27
 
  uint sum = 0;
28
 
  for (uint i = 0; i < ode.size(); i++)
29
 
  {
30
 
    Jindices[i] = sum;
31
 
    sum += ode.dependencies[i].size();
32
 
  }
33
 
 
34
 
  // Allocate Jacobian values
35
 
  Jvalues = new real[sum];
36
 
  for (uint pos = 0; pos < sum; pos++)
37
 
    Jvalues[pos] = 0.0;
38
 
 
39
 
  message("Generated Jacobian data structure for %d dependencies.", sum);
40
 
 
41
 
  // Compute maximum number of dependencies
42
 
  uint maxsize = 0;
43
 
  for (uint i = 0; i < ode.size(); i++)
44
 
  {
45
 
    const uint size = ode.dependencies[i].size();
46
 
    if ( size > maxsize )
47
 
      maxsize = size;
48
 
  }
49
 
 
50
 
  // Allocate lookup table for dependencies to components with small time steps
51
 
  Jlookup = new real[std::max(static_cast<unsigned int>(1), maxsize - 1)];
52
 
}
53
 
//-----------------------------------------------------------------------------
54
 
MultiAdaptiveJacobian::~MultiAdaptiveJacobian()
55
 
{
56
 
  if ( Jvalues ) delete [] Jvalues;
57
 
  if ( Jindices ) delete [] Jindices;
58
 
  if ( Jlookup ) delete [] Jlookup;
59
 
}
60
 
//-----------------------------------------------------------------------------
61
 
dolfin::uint MultiAdaptiveJacobian::size(uint dim) const
62
 
{
63
 
  return ts.nj;
64
 
}
65
 
//-----------------------------------------------------------------------------
66
 
void MultiAdaptiveJacobian::mult(const uBlasVector& x, uBlasVector& y) const
67
 
{
68
 
  // We iterate over all degrees of freedom j in the time slab and compute
69
 
  // y_j = (Ax)_j for each degree of freedom of the system.
70
 
 
71
 
  // Start with y = x, accounting for the derivative dF_j/dx_j = 1
72
 
  y = x;
73
 
 
74
 
  // Choose method
75
 
  if ( method.type() == Method::cG )
76
 
    cGmult(x, y);
77
 
  else
78
 
    dGmult(x, y);
79
 
}
80
 
//-----------------------------------------------------------------------------
81
 
void MultiAdaptiveJacobian::init()
82
 
{
83
 
  // Compute Jacobian at the beginning of the slab
84
 
  real t = ts.starttime();
85
 
  //message("Recomputing Jacobian matrix at t = %f.", t);
86
 
  
87
 
  // Compute Jacobian
88
 
  for (uint i = 0; i < ode.size(); i++)
89
 
  {
90
 
    const Array<uint>& deps = ode.dependencies[i];
91
 
    for (uint pos = 0; pos < deps.size(); pos++)
92
 
      Jvalues[Jindices[i] + pos] = ode.dfdu(ts.u0, t, i, deps[pos]);
93
 
  }
94
 
 
95
 
  /*
96
 
  // Compute Jacobian at the end of the slab
97
 
  real t = ts.endtime();
98
 
  //message("Recomputing Jacobian matrix at t = %f.", t);
99
 
  
100
 
  // Compute Jacobian
101
 
  for (uint i = 0; i < ode.size(); i++)
102
 
  {
103
 
    const Array<uint>& deps = ode.dependencies[i];
104
 
    for (uint pos = 0; pos < deps.size(); pos++)
105
 
      Jvalues[Jindices[i] + pos] = ode.dfdu(ts.u, t, i, deps[pos]);
106
 
  }
107
 
  */
108
 
}
109
 
//-----------------------------------------------------------------------------
110
 
void MultiAdaptiveJacobian::cGmult(const uBlasVector& x, uBlasVector& y) const
111
 
{
112
 
  // Reset current sub slab
113
 
  int s0 = -1;
114
 
 
115
 
  // Reset elast
116
 
  for (uint i = 0; i < ts.N; i++)
117
 
    ts.elast[i] = -1;
118
 
 
119
 
  // Iterate over all elements
120
 
  for (uint e0 = 0; e0 < ts.ne; e0++)
121
 
  {
122
 
    // Cover all elements in current sub slab
123
 
    s0 = ts.coverNext(s0, e0);
124
 
    
125
 
    // Get element data
126
 
    const uint i0 = ts.ei[e0];
127
 
    const real a0 = ts.sa[s0];
128
 
    const real b0 = ts.sb[s0];
129
 
    const real k0 = b0 - a0;
130
 
    const uint j0 = e0 * method.nsize();
131
 
    
132
 
    // Add dependency on predecessor for all dofs of element
133
 
    const int ep = ts.ee[e0];
134
 
    if ( ep != -1 )
135
 
    {
136
 
      const real xp = x(ep*method.nsize() + method.nsize() - 1);
137
 
      for (uint n = 0; n < method.nsize(); n++)
138
 
        y(j0 + n) -= xp;
139
 
    }
140
 
 
141
 
    // Reset Jpos
142
 
    uint Jpos = 0;
143
 
 
144
 
    // Iterate over dependencies for the current component
145
 
    const Array<uint>& deps = ode.dependencies[i0];
146
 
    for (uint pos = 0; pos < deps.size(); pos++)
147
 
    {
148
 
      // Get derivative
149
 
      const real dfdu = Jvalues[Jindices[i0] + pos];
150
 
 
151
 
      // Skip elements which have not been covered
152
 
      const uint i1 = deps[pos];
153
 
      const int e1 = ts.elast[i1];      
154
 
      if ( e1 == -1 )
155
 
      {
156
 
        // Save dependency for later
157
 
        Jlookup[Jpos++] = dfdu;
158
 
        continue;
159
 
      }
160
 
 
161
 
      // Skip elements with smaller time steps
162
 
      const uint s1 = ts.es[e1];
163
 
      const real b1 = ts.sb[s1];
164
 
      if ( b1 < (a0 + DOLFIN_EPS) )
165
 
      {
166
 
        // Save dependency for later
167
 
        Jlookup[Jpos++] = dfdu;
168
 
 
169
 
        // Add dependency to initial value for all dofs
170
 
        const real tmp = k0 * dfdu * x(e1 * method.nsize() + method.nsize() - 1);
171
 
        for (uint n = 0; n < method.nsize(); n++)
172
 
          y(j0 + n) -= tmp * method.nweight(n, 0);
173
 
 
174
 
        continue;
175
 
      }
176
 
      
177
 
      // Get first dof for other element
178
 
      const uint j1 = e1 * method.nsize();
179
 
      
180
 
      // Use fast evaluation for elements in the same sub slab
181
 
      if ( s0 == static_cast<int>(s1) )
182
 
      {
183
 
        // Add dependency to dof of initial value if any
184
 
        const int ep = ts.ee[e1];
185
 
        const real tmp0 = k0 * dfdu;
186
 
        if ( ep != -1 )
187
 
        {
188
 
          const real tmp1 = tmp0 * x(ep * method.nsize() + method.nsize() - 1);
189
 
          for (uint n = 0; n < method.nsize(); n++)
190
 
            y(j0 + n) -= tmp1 * method.nweight(n, 0);
191
 
        }
192
 
        
193
 
        // Add dependencies to internal dofs
194
 
        for (uint n = 0; n < method.nsize(); n++)
195
 
        {
196
 
          real sum = 0.0;
197
 
          for (uint m = 0; m < method.nsize(); m++)
198
 
            sum += method.nweight(n, m + 1) * x(j1 + m);
199
 
          y(j0 + n) -= tmp0 * sum;
200
 
        }
201
 
      }
202
 
      else
203
 
      {
204
 
        const real a1 = ts.sa[s1];
205
 
        const real k1 = b1 - a1;
206
 
        
207
 
        // Iterate over dofs of element
208
 
        const real tmp0 = k0 * dfdu;
209
 
        for (uint n = 0; n < method.nsize(); n++)
210
 
        {
211
 
          // Iterate over quadrature points
212
 
          real sum = 0.0;
213
 
          for (uint m = 0; m < method.qsize(); m++)
214
 
          {
215
 
            const real tau = (a0 + k0*method.qpoint(m) - a1) / k1;
216
 
            const real tmp1 = method.nweight(n, m);
217
 
            dolfin_assert(tau >= -DOLFIN_EPS);
218
 
            dolfin_assert(tau <= 1.0 + DOLFIN_EPS);
219
 
            
220
 
            // Add dependency to dof of initial value if any
221
 
            const int ep = ts.ee[e1];
222
 
            if ( ep != -1 )
223
 
            {
224
 
              const real x0 = x(ep * method.nsize() + method.nsize() - 1);
225
 
              sum += tmp1 * method.eval(0, tau) * x0;
226
 
            }
227
 
            
228
 
            // Iterate over dofs of other element and add dependencies
229
 
            for (uint l = 0; l < method.nsize(); l++)
230
 
              sum += tmp1 * method.eval(l + 1, tau) * x(j1 + l);
231
 
          }
232
 
            
233
 
          // Add dependencies
234
 
          y(j0 + n) -= tmp0 * sum;
235
 
        }
236
 
      }      
237
 
    }
238
 
    
239
 
    // At this point, we have added dependencies to the predecessor,
240
 
    // to dofs in the same sub slab and to components with larger time
241
 
    // steps. It remains to add dependencies to components with
242
 
    // smaller time steps. We need to do this by iterating over
243
 
    // quadrature points, since this is the order in which the
244
 
    // dependencies are stored.
245
 
 
246
 
    // Get first dependency to components with smaller time steps for element
247
 
    uint d = ts.ed[e0];
248
 
    
249
 
    // Compute number of such dependencies for each nodal point
250
 
    const uint end = ( e0 < (ts.ne - 1) ? ts.ed[e0 + 1] : ts.nd );
251
 
    const uint ndep = (end - d) / method.nsize();
252
 
    dolfin_assert(ndep * method.nsize() == (end - d));
253
 
 
254
 
    // Iterate over quadrature points of current element
255
 
    for (uint m = 1; m < method.qsize(); m++)
256
 
    {
257
 
      // Compute quadrature point
258
 
      const real t = a0 + k0*method.qpoint(m);
259
 
 
260
 
      // Iterate over dependencies
261
 
      for (uint dep = 0; dep < ndep; dep++)
262
 
      {
263
 
        // Get element data
264
 
        const uint e1 = ts.de[d++];
265
 
        const uint j1 = e1 * method.nsize();
266
 
        const uint s1 = ts.es[e1];
267
 
        //const uint i1 = ts.ei[e1];
268
 
        const real a1 = ts.sa[s1];
269
 
        const real b1 = ts.sb[s1];
270
 
        const real k1 = b1 - a1;
271
 
        const real tau = (t - a1) / k1;
272
 
        
273
 
        // We don't know how to index Jvalues here and want to avoid
274
 
        // searching, but we were clever enough to pick out the value
275
 
        // before when we had the chance... :-)
276
 
        const real dfdu = Jlookup[dep];
277
 
        //message("Looks like df_%d/du_%d = %f", i0, i1, dfdu);      
278
 
        
279
 
        // Iterate over quadrature points of other element
280
 
        const real tmp0 = k0 * dfdu;
281
 
        for (uint l = 0; l < method.qsize(); l++)
282
 
        {
283
 
          real tmp1 = tmp0 * method.eval(l, tau);
284
 
          if ( l == 0 )
285
 
          {
286
 
            const int ep = ts.ee[e1];
287
 
            if ( ep != -1 )
288
 
            {
289
 
              // Iterate over dofs of current element
290
 
              tmp1 *= x(ep*method.nsize() + method.nsize() - 1);
291
 
              for (uint n = 0; n < method.nsize(); n++)
292
 
                y(j0 + n) -= tmp1 * method.nweight(n, m);               
293
 
            }
294
 
          }
295
 
          else
296
 
          {
297
 
            // Iterate over dofs of current element
298
 
            tmp1 *= x(j1 + l - 1);
299
 
            for (uint n = 0; n < method.nsize(); n++)
300
 
              y(j0 + n) -= tmp1 * method.nweight(n, m);         
301
 
          }
302
 
        }
303
 
      }
304
 
    }
305
 
  }
306
 
}
307
 
//-----------------------------------------------------------------------------
308
 
void MultiAdaptiveJacobian::dGmult(const uBlasVector& x, uBlasVector& y) const
309
 
{
310
 
  // Reset current sub slab
311
 
  int s0 = -1;
312
 
 
313
 
  // Reset elast
314
 
  for (uint i = 0; i < ts.N; i++)
315
 
    ts.elast[i] = -1;
316
 
 
317
 
  // Iterate over all elements
318
 
  for (uint e0 = 0; e0 < ts.ne; e0++)
319
 
  {
320
 
    // Cover all elements in current sub slab
321
 
    s0 = ts.coverNext(s0, e0);
322
 
    
323
 
    // Get element data
324
 
    const uint i0 = ts.ei[e0];
325
 
    const real a0 = ts.sa[s0];
326
 
    const real b0 = ts.sb[s0];
327
 
    const real k0 = b0 - a0;
328
 
    const uint j0 = e0 * method.nsize();
329
 
    
330
 
    // Add dependency on predecessor for all dofs of element
331
 
    const int ep = ts.ee[e0];
332
 
    if ( ep != -1 )
333
 
    {
334
 
      const real xp = x(ep*method.nsize() + method.nsize() - 1);
335
 
      for (uint n = 0; n < method.nsize(); n++)
336
 
        y(j0 + n) -= xp;
337
 
    }
338
 
 
339
 
    // Reset Jpos
340
 
    uint Jpos = 0;
341
 
 
342
 
    // Iterate over dependencies for the current component
343
 
    const Array<uint>& deps = ode.dependencies[i0];
344
 
    for (uint pos = 0; pos < deps.size(); pos++)
345
 
    {
346
 
      // Get derivative
347
 
      const real dfdu = Jvalues[Jindices[i0] + pos];
348
 
 
349
 
      // Skip elements which have not been covered
350
 
      const uint i1 = deps[pos];
351
 
      const int e1 = ts.elast[i1];      
352
 
      if ( e1 == -1 )
353
 
      {
354
 
        // Save dependency for later
355
 
        Jlookup[Jpos++] = dfdu;
356
 
        continue;
357
 
      }
358
 
 
359
 
      // Skip elements with smaller time steps
360
 
      const uint s1 = ts.es[e1];
361
 
      const real b1 = ts.sb[s1];
362
 
      if ( b1 < (a0 + DOLFIN_EPS) )
363
 
      {
364
 
        // Save dependency for later
365
 
        Jlookup[Jpos++] = dfdu;
366
 
        continue;
367
 
      }
368
 
      
369
 
      // Get first dof for other element
370
 
      const uint j1 = e1 * method.nsize();
371
 
      
372
 
      // Use fast evaluation for elements in the same sub slab
373
 
      if ( s0 == static_cast<int>(s1) )
374
 
      {
375
 
        const real tmp = k0 * dfdu;
376
 
        for (uint n = 0; n < method.nsize(); n++)
377
 
        {
378
 
          real sum = 0.0;
379
 
          for (uint m = 0; m < method.qsize(); m++)
380
 
            sum += method.nweight(n, m) * x(j1 + m);
381
 
          y(j0 + n) -= tmp * sum;
382
 
        }
383
 
      }
384
 
      else
385
 
      {
386
 
        const real a1 = ts.sa[s1];
387
 
        const real k1 = b1 - a1;
388
 
        
389
 
        // Iterate over dofs of element
390
 
        const real tmp0 = k0 * dfdu;
391
 
        for (uint n = 0; n < method.nsize(); n++)
392
 
        {
393
 
          // Iterate over quadrature points
394
 
          real sum = 0.0;
395
 
          for (uint m = 0; m < method.qsize(); m++)
396
 
          {
397
 
            const real tau = (a0 + k0*method.qpoint(m) - a1) / k1;
398
 
            const real tmp1 = method.nweight(n, m);
399
 
            
400
 
            // Iterate over dofs of other element and add dependencies
401
 
            for (uint l = 0; l < method.nsize(); l++)
402
 
              sum += tmp1 * method.eval(l, tau) * x(j1 + l);
403
 
          }
404
 
          
405
 
          // Add dependencies
406
 
          y(j0 + n) -= tmp0 * sum;
407
 
        }
408
 
      }
409
 
    }
410
 
    
411
 
    // At this point, we have added dependencies to the predecessor,
412
 
    // to dofs in the same sub slab and to components with larger time
413
 
    // steps. It remains to add dependencies to components with
414
 
    // smaller time steps. We need to do this by iterating over
415
 
    // quadrature points, since this is the order in which the
416
 
    // dependencies are stored.
417
 
 
418
 
    // Get first dependency to components with smaller time steps for element
419
 
    uint d = ts.ed[e0];
420
 
    
421
 
    // Compute number of such dependencies for each nodal point
422
 
    const uint end = ( e0 < (ts.ne - 1) ? ts.ed[e0 + 1] : ts.nd );
423
 
    const uint ndep = (end - d) / method.nsize();
424
 
    dolfin_assert(ndep * method.nsize() == (end - d));
425
 
 
426
 
    // Iterate over quadrature points of current element
427
 
    for (uint m = 0; m < method.qsize(); m++)
428
 
    {
429
 
      // Compute quadrature point
430
 
      const real t = a0 + k0*method.qpoint(m);
431
 
 
432
 
      // Iterate over dependencies
433
 
      for (uint dep = 0; dep < ndep; dep++)
434
 
      {
435
 
        // Get element data
436
 
        const uint e1 = ts.de[d++];
437
 
        const uint j1 = e1 * method.nsize();
438
 
        const uint s1 = ts.es[e1];
439
 
        //const uint i1 = ts.ei[e1];
440
 
        const real a1 = ts.sa[s1];
441
 
        const real b1 = ts.sb[s1];
442
 
        const real k1 = b1 - a1;
443
 
        const real tau = (t - a1) / k1;
444
 
        
445
 
        // We don't know how to index Jvalues here and want to avoid
446
 
        // searching, but we were clever enough to pick out the value
447
 
        // before when we had the chance... :-)
448
 
        const real dfdu = Jlookup[dep];
449
 
        //message("Looks like df_%d/du_%d = %f", i0, i1, dfdu);      
450
 
        
451
 
        // Iterate over quadrature points of other element
452
 
        const real tmp0 = k0 * dfdu;
453
 
        for (uint l = 0; l < method.qsize(); l++)
454
 
        {
455
 
          real tmp1 = tmp0 * method.eval(l, tau) * x(j1 + l);
456
 
          for (uint n = 0; n < method.nsize(); n++)
457
 
            y(j0 + n) -= tmp1 * method.nweight(n, m);
458
 
        }
459
 
      }
460
 
    }
461
 
  }
462
 
}
463
 
//-----------------------------------------------------------------------------