← Back to branch summary

~njansson/dolfin/hpc

« back to all changes in this revision

Viewing changes to src/kernel/ode/MultiAdaptiveJacobian.cpp

  • Committer: Anders Logg
  • Date: 2007-01-10 09:04:44 UTC
  • mfrom: (1689.1.221 trunk)
  • Revision ID: logg@simula.no-20070110090444-ecyux3n1qnei4i8h
Remove old head

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/kernel/ode/MultiAdaptiveJacobian.cpp
1
 
// Copyright (C) 2005 Anders Logg.
 
1
// Copyright (C) 2005-2006 Anders Logg.
2
2
// Licensed under the GNU GPL Version 2.
3
3
//
4
4
// First added:  2005-01-27
5
 
// Last changed: 2005-11-10
 
5
// Last changed: 2006-08-08
6
6
 
7
7
#include <dolfin/dolfin_math.h>
 
8
#include <dolfin/uBlasVector.h>
8
9
#include <dolfin/ODE.h>
9
 
#include <dolfin/Vector.h>
10
10
#include <dolfin/Method.h>
11
11
#include <dolfin/MultiAdaptiveTimeSlab.h>
12
12
#include <dolfin/MultiAdaptiveNewtonSolver.h>
58
58
  if ( Jlookup ) delete [] Jlookup;
59
59
}
60
60
//-----------------------------------------------------------------------------
61
 
void MultiAdaptiveJacobian::update()
 
61
dolfin::uint MultiAdaptiveJacobian::size(const 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()
62
82
{
63
83
  // Compute Jacobian at the beginning of the slab
64
84
  real t = ts.starttime();
87
107
  */
88
108
}
89
109
//-----------------------------------------------------------------------------
90
 
void MultiAdaptiveJacobian::mult(const Vector& x, Vector& y) const
91
 
{
92
 
  // We iterate over all degrees of freedom j in the time slab and compute
93
 
  // y_j = (Ax)_j for each degree of freedom of the system. Note that this
94
 
  // implementation will probably not work with parallel vectors since we
95
 
  // use VecGetArray to access the local arrays of the vectors
96
 
 
97
 
  // Start with y = x, accounting for the derivative dF_j/dx_j = 1
98
 
  y = x;
99
 
 
100
 
  // Get data arrays (assumes uniprocessor case)
101
 
  const real* xx = x.array();
102
 
  real* yy = y.array();
103
 
  
104
 
  // Choose method
105
 
  if ( method.type() == Method::cG )
106
 
  cGmult(xx, yy);
107
 
  else
108
 
  dGmult(xx, yy);
109
 
 
110
 
  // Restore data arrays
111
 
  x.restore(xx);
112
 
  y.restore(yy);
113
 
}
114
 
//-----------------------------------------------------------------------------
115
 
void MultiAdaptiveJacobian::cGmult(const real x[], real y[]) const
 
110
void MultiAdaptiveJacobian::cGmult(const uBlasVector& x, uBlasVector& y) const
116
111
{
117
112
  // Reset current sub slab
118
113
  int s0 = -1;
138
133
    const int ep = ts.ee[e0];
139
134
    if ( ep != -1 )
140
135
    {
141
 
      const real xp = x[ep*method.nsize() + method.nsize() - 1];
 
136
      const real xp = x(ep*method.nsize() + method.nsize() - 1);
142
137
      for (uint n = 0; n < method.nsize(); n++)
143
 
        y[j0 + n] -= xp;
 
138
        y(j0 + n) -= xp;
144
139
    }
145
140
 
146
141
    // Reset Jpos
172
167
        Jlookup[Jpos++] = dfdu;
173
168
 
174
169
        // Add dependency to initial value for all dofs
175
 
        const real tmp = k0 * dfdu * x[e1 * method.nsize() + method.nsize() - 1];
 
170
        const real tmp = k0 * dfdu * x(e1 * method.nsize() + method.nsize() - 1);
176
171
        for (uint n = 0; n < method.nsize(); n++)
177
 
          y[j0 + n] -= tmp * method.nweight(n, 0);
 
172
          y(j0 + n) -= tmp * method.nweight(n, 0);
178
173
 
179
174
        continue;
180
175
      }
190
185
        const real tmp0 = k0 * dfdu;
191
186
        if ( ep != -1 )
192
187
        {
193
 
          const real tmp1 = tmp0 * x[ep * method.nsize() + method.nsize() - 1];
 
188
          const real tmp1 = tmp0 * x(ep * method.nsize() + method.nsize() - 1);
194
189
          for (uint n = 0; n < method.nsize(); n++)
195
 
            y[j0 + n] -= tmp1 * method.nweight(n, 0);
 
190
            y(j0 + n) -= tmp1 * method.nweight(n, 0);
196
191
        }
197
192
        
198
193
        // Add dependencies to internal dofs
200
195
        {
201
196
          real sum = 0.0;
202
197
          for (uint m = 0; m < method.nsize(); m++)
203
 
            sum += method.nweight(n, m + 1) * x[j1 + m];
204
 
          y[j0 + n] -= tmp0 * sum;
 
198
            sum += method.nweight(n, m + 1) * x(j1 + m);
 
199
          y(j0 + n) -= tmp0 * sum;
205
200
        }
206
201
      }
207
202
      else
226
221
            const int ep = ts.ee[e1];
227
222
            if ( ep != -1 )
228
223
            {
229
 
              const real x0 = x[ep * method.nsize() + method.nsize() - 1];
 
224
              const real x0 = x(ep * method.nsize() + method.nsize() - 1);
230
225
              sum += tmp1 * method.eval(0, tau) * x0;
231
226
            }
232
227
            
233
228
            // Iterate over dofs of other element and add dependencies
234
229
            for (uint l = 0; l < method.nsize(); l++)
235
 
              sum += tmp1 * method.eval(l + 1, tau) * x[j1 + l];
 
230
              sum += tmp1 * method.eval(l + 1, tau) * x(j1 + l);
236
231
          }
237
232
            
238
233
          // Add dependencies
239
 
          y[j0 + n] -= tmp0 * sum;
 
234
          y(j0 + n) -= tmp0 * sum;
240
235
        }
241
236
      }      
242
237
    }
292
287
            if ( ep != -1 )
293
288
            {
294
289
              // Iterate over dofs of current element
295
 
              tmp1 *= x[ep*method.nsize() + method.nsize() - 1];
 
290
              tmp1 *= x(ep*method.nsize() + method.nsize() - 1);
296
291
              for (uint n = 0; n < method.nsize(); n++)
297
 
                y[j0 + n] -= tmp1 * method.nweight(n, m);               
 
292
                y(j0 + n) -= tmp1 * method.nweight(n, m);               
298
293
            }
299
294
          }
300
295
          else
301
296
          {
302
297
            // Iterate over dofs of current element
303
 
            tmp1 *= x[j1 + l - 1];
 
298
            tmp1 *= x(j1 + l - 1);
304
299
            for (uint n = 0; n < method.nsize(); n++)
305
 
              y[j0 + n] -= tmp1 * method.nweight(n, m);         
 
300
              y(j0 + n) -= tmp1 * method.nweight(n, m);         
306
301
          }
307
302
        }
308
303
      }
310
305
  }
311
306
}
312
307
//-----------------------------------------------------------------------------
313
 
void MultiAdaptiveJacobian::dGmult(const real x[], real y[]) const
 
308
void MultiAdaptiveJacobian::dGmult(const uBlasVector& x, uBlasVector& y) const
314
309
{
315
310
  // Reset current sub slab
316
311
  int s0 = -1;
336
331
    const int ep = ts.ee[e0];
337
332
    if ( ep != -1 )
338
333
    {
339
 
      const real xp = x[ep*method.nsize() + method.nsize() - 1];
 
334
      const real xp = x(ep*method.nsize() + method.nsize() - 1);
340
335
      for (uint n = 0; n < method.nsize(); n++)
341
 
        y[j0 + n] -= xp;
 
336
        y(j0 + n) -= xp;
342
337
    }
343
338
 
344
339
    // Reset Jpos
382
377
        {
383
378
          real sum = 0.0;
384
379
          for (uint m = 0; m < method.qsize(); m++)
385
 
            sum += method.nweight(n, m) * x[j1 + m];
386
 
          y[j0 + n] -= tmp * sum;
 
380
            sum += method.nweight(n, m) * x(j1 + m);
 
381
          y(j0 + n) -= tmp * sum;
387
382
        }
388
383
      }
389
384
      else
404
399
            
405
400
            // Iterate over dofs of other element and add dependencies
406
401
            for (uint l = 0; l < method.nsize(); l++)
407
 
              sum += tmp1 * method.eval(l, tau) * x[j1 + l];
 
402
              sum += tmp1 * method.eval(l, tau) * x(j1 + l);
408
403
          }
409
404
          
410
405
          // Add dependencies
411
 
          y[j0 + n] -= tmp0 * sum;
 
406
          y(j0 + n) -= tmp0 * sum;
412
407
        }
413
408
      }
414
409
    }
457
452
        const real tmp0 = k0 * dfdu;
458
453
        for (uint l = 0; l < method.qsize(); l++)
459
454
        {
460
 
          real tmp1 = tmp0 * method.eval(l, tau) * x[j1 + l];
 
455
          real tmp1 = tmp0 * method.eval(l, tau) * x(j1 + l);
461
456
          for (uint n = 0; n < method.nsize(); n++)
462
 
            y[j0 + n] -= tmp1 * method.nweight(n, m);
 
457
            y(j0 + n) -= tmp1 * method.nweight(n, m);
463
458
        }
464
459
      }
465
460
    }