← Back to branch summary

~njansson/dolfin/hpc

« back to all changes in this revision

Viewing changes to src/kernel/ode/ComplexODE.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/ComplexODE.cpp
1
 
// Copyright (C) 2005-2006 Anders Logg.
2
 
// Licensed under the GNU LGPL Version 2.1.
3
 
//
4
 
// First added:  2005-02-02
5
 
// Last changed: 2006-08-21
6
 
 
7
 
#include <dolfin/Array.h>
8
 
#include <dolfin/ComplexODE.h>
9
 
 
10
 
using namespace dolfin;
11
 
 
12
 
//-----------------------------------------------------------------------------
13
 
ComplexODE::ComplexODE(uint n, real T) : ODE(2*n, T), n(n), j(0.0, 1.0),
14
 
                                         zvalues(0), fvalues(0), yvalues(0)
15
 
{
16
 
  message("Creating complex ODE of size %d (%d complex components).", N, n);
17
 
  
18
 
  // Initialize complex solution vector and right-hand side
19
 
  zvalues = new complex[n];
20
 
  fvalues = new complex[n];
21
 
  for (uint i = 0; i < n; i++)
22
 
  {
23
 
    zvalues[i] = 0.0;
24
 
    fvalues[i] = 0.0;
25
 
  }
26
 
}
27
 
//-----------------------------------------------------------------------------
28
 
ComplexODE::~ComplexODE()
29
 
{
30
 
  if ( zvalues ) delete [] zvalues;
31
 
  if ( fvalues ) delete [] fvalues;
32
 
  if ( yvalues ) delete [] yvalues;
33
 
}
34
 
//-----------------------------------------------------------------------------
35
 
complex ComplexODE::f(const complex z[], real t, uint i)
36
 
{
37
 
  error("Right-hand side for complex ODE not supplied by user.");
38
 
  
39
 
  complex zvalue;
40
 
  return zvalue;
41
 
}
42
 
//-----------------------------------------------------------------------------
43
 
void ComplexODE::f(const complex z[], real t, complex y[])
44
 
{
45
 
  // If a user of the mono-adaptive solver does not supply this function,
46
 
  // then call f() for each component.
47
 
 
48
 
  for (uint i = 0; i < n; i++)
49
 
    y[i] = this->f(z, t, i);
50
 
}
51
 
//-----------------------------------------------------------------------------
52
 
void ComplexODE::M(const complex x[], complex y[], const complex z[], real t)
53
 
{
54
 
  // Assume M is the identity if not supplied by user: y = x
55
 
  for (uint i = 0; i < n; i++)
56
 
    y[i] = x[i];
57
 
}
58
 
//-----------------------------------------------------------------------------
59
 
void ComplexODE::J(const complex x[], complex y[], const complex z[], real t)
60
 
{
61
 
  // If a user does not supply J, then compute it by the approximation
62
 
  //
63
 
  //     Jx = ( f(z + hx) - f(z - hx) ) / 2h
64
 
 
65
 
  error("Not implemented yet...");
66
 
}
67
 
//-----------------------------------------------------------------------------
68
 
real ComplexODE::k(uint i)
69
 
{
70
 
  return default_timestep;
71
 
}
72
 
//-----------------------------------------------------------------------------
73
 
bool ComplexODE::update(const complex z[], real t, bool end)
74
 
{
75
 
  return true;
76
 
}
77
 
//-----------------------------------------------------------------------------
78
 
void ComplexODE::u0(uBlasVector& u)
79
 
{
80
 
  // Translate initial value from complex to real
81
 
  z0(zvalues);
82
 
  for (uint i = 0; i < N; i++)
83
 
  {
84
 
    complex z = zvalues[i / 2];
85
 
    u(i) = ( i % 2 == 0 ? z.real() : z.imag() );
86
 
  }
87
 
}
88
 
//-----------------------------------------------------------------------------
89
 
real ComplexODE::f(const uBlasVector& u, real t, uint i)
90
 
{
91
 
  // Translate right-hand side from complex to real, assuming that if
92
 
  // u_i depends on u_j, then u_i depends on both the real and
93
 
  // imaginary parts of the corresponding z_i
94
 
 
95
 
  // Update zvalues for correct components
96
 
  const Array<uint>& deps = dependencies[i];
97
 
  for (uint pos = 0; pos < deps.size(); pos++)
98
 
  {
99
 
    // Element of u that needs to be updated
100
 
    const uint ju = deps[pos];
101
 
 
102
 
    // Corresponding element of z
103
 
    const uint jz = ju / 2;
104
 
 
105
 
    // Update value
106
 
    const complex zvalue(u(2*jz), u(2*jz + 1));
107
 
    zvalues[jz] = zvalue;
108
 
  }
109
 
 
110
 
  // Call user-supplied function f(z, t, i)
111
 
  const complex fvalue = f(zvalues, t, i / 2);
112
 
  
113
 
  // Return value
114
 
  return ( i % 2 == 0 ? fvalue.real() : fvalue.imag() );
115
 
}
116
 
//-----------------------------------------------------------------------------
117
 
void ComplexODE::f(const uBlasVector& u, real t, uBlasVector& y)
118
 
{
119
 
  // Update zvalues for all components
120
 
  for (uint i = 0; i < n; i++)
121
 
  {
122
 
    const complex zvalue(u(2*i), u(2*i + 1));
123
 
    zvalues[i] = zvalue;
124
 
  }
125
 
 
126
 
  // Call user-supplied function f(z, t, y)
127
 
  f(zvalues, t, fvalues);
128
 
  
129
 
  // Translate values into f
130
 
  for (uint i = 0; i < n; i++)
131
 
  {
132
 
    const complex fvalue = fvalues[i];
133
 
    y(2*i) = fvalue.real();
134
 
    y(2*i + 1) = fvalue.imag();
135
 
  }
136
 
}
137
 
//-----------------------------------------------------------------------------
138
 
void ComplexODE::M(const uBlasVector& x, uBlasVector& y,
139
 
                   const uBlasVector& u, real t)
140
 
{
141
 
  // Update zvalues and fvalues for all components
142
 
  for (uint i = 0; i < n; i++)
143
 
  {
144
 
    const complex zvalue(u(2*i), u(2*i + 1));
145
 
    const complex xvalue(x(2*i), x(2*i + 1));
146
 
    zvalues[i] = zvalue;
147
 
    fvalues[i] = xvalue; // Use fvalues for x
148
 
  }
149
 
 
150
 
  // Use additional array for y, initialize the first time
151
 
  if ( !yvalues )
152
 
  {
153
 
    yvalues = new complex[n];
154
 
    for (uint i = 0; i < n; i++)
155
 
      yvalues[i] = 0.0;
156
 
  }
157
 
  
158
 
  // Call user-supplied function M(x, y, z, t)
159
 
  M(fvalues, yvalues, zvalues, t);
160
 
 
161
 
  // Copy values to y
162
 
  for (uint i = 0; i < n; i++)
163
 
  {
164
 
    const complex yvalue = yvalues[i];
165
 
    y(2*i) = yvalue.real();
166
 
    y(2*i + 1) = yvalue.imag();
167
 
  }
168
 
}
169
 
//-----------------------------------------------------------------------------
170
 
void ComplexODE::J(const uBlasVector& x, uBlasVector& y,
171
 
                   const uBlasVector& u, real t)
172
 
{
173
 
  // Update zvalues and fvalues for all components
174
 
  for (uint i = 0; i < n; i++)
175
 
  {
176
 
    const complex zvalue(u(2*i), u(2*i + 1));
177
 
    const complex xvalue(x(2*i), x(2*i + 1));
178
 
    zvalues[i] = zvalue;
179
 
    fvalues[i] = xvalue; // Use fvalues for x
180
 
  }
181
 
 
182
 
  // Use additional array for y, initialize the first time
183
 
  if ( !yvalues )
184
 
  {
185
 
    yvalues = new complex[n];
186
 
    for (uint i = 0; i < n; i++)
187
 
      yvalues[i] = 0.0;
188
 
  }
189
 
  
190
 
  // Call user-supplied function J(x, y, z, t)
191
 
  J(fvalues, yvalues, zvalues, t);
192
 
 
193
 
  // Copy values to y
194
 
  for (uint i = 0; i < n; i++)
195
 
  {
196
 
    const complex yvalue = yvalues[i];
197
 
    y(2*i) = yvalue.real();
198
 
    y(2*i + 1) = yvalue.imag();
199
 
  }
200
 
}
201
 
//-----------------------------------------------------------------------------
202
 
real ComplexODE::timestep(uint i)
203
 
{
204
 
  // Translate time step
205
 
  return k(i / 2);
206
 
}
207
 
//-----------------------------------------------------------------------------
208
 
bool ComplexODE::update(const uBlasVector& u, real t, bool end)
209
 
{
210
 
  // Update zvalues for all components
211
 
  for (uint i = 0; i < n; i++)
212
 
  {
213
 
    const complex zvalue(u(2*i), u(2*i + 1));
214
 
    zvalues[i] = zvalue;
215
 
  }
216
 
 
217
 
  // Call user-supplied function update(z, t)
218
 
  return update(zvalues, t, end);
219
 
}
220
 
//-----------------------------------------------------------------------------