MFEM  v4.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
ode.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
2 // the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
3 // reserved. See file COPYRIGHT for details.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability see http://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the GNU Lesser General Public License (as published by the Free
10 // Software Foundation) version 2.1 dated February 1999.
11 
12 #include "operator.hpp"
13 #include "ode.hpp"
14 
15 namespace mfem
16 {
17 
19 {
20  this->f = &f;
22 }
23 
25 {
26  ODESolver::Init(_f);
27  dxdt.SetSize(f->Width(), mem_type);
28 }
29 
30 void ForwardEulerSolver::Step(Vector &x, double &t, double &dt)
31 {
32  f->SetTime(t);
33  f->Mult(x, dxdt);
34  x.Add(dt, dxdt);
35  t += dt;
36 }
37 
38 
40 {
41  ODESolver::Init(_f);
42  int n = f->Width();
43  dxdt.SetSize(n, mem_type);
44  x1.SetSize(n, mem_type);
45 }
46 
47 void RK2Solver::Step(Vector &x, double &t, double &dt)
48 {
49  // 0 |
50  // a | a
51  // ---+--------
52  // | 1-b b b = 1/(2a)
53 
54  const double b = 0.5/a;
55 
56  f->SetTime(t);
57  f->Mult(x, dxdt);
58  add(x, (1. - b)*dt, dxdt, x1);
59  x.Add(a*dt, dxdt);
60 
61  f->SetTime(t + a*dt);
62  f->Mult(x, dxdt);
63  add(x1, b*dt, dxdt, x);
64  t += dt;
65 }
66 
67 
69 {
70  ODESolver::Init(_f);
71  int n = f->Width();
72  y.SetSize(n, mem_type);
73  k.SetSize(n, mem_type);
74 }
75 
76 void RK3SSPSolver::Step(Vector &x, double &t, double &dt)
77 {
78  // x0 = x, t0 = t, k0 = dt*f(t0, x0)
79  f->SetTime(t);
80  f->Mult(x, k);
81 
82  // x1 = x + k0, t1 = t + dt, k1 = dt*f(t1, x1)
83  add(x, dt, k, y);
84  f->SetTime(t + dt);
85  f->Mult(y, k);
86 
87  // x2 = 3/4*x + 1/4*(x1 + k1), t2 = t + 1/2*dt, k2 = dt*f(t2, x2)
88  y.Add(dt, k);
89  add(3./4, x, 1./4, y, y);
90  f->SetTime(t + dt/2);
91  f->Mult(y, k);
92 
93  // x3 = 1/3*x + 2/3*(x2 + k2), t3 = t + dt
94  y.Add(dt, k);
95  add(1./3, x, 2./3, y, x);
96  t += dt;
97 }
98 
99 
101 {
102  ODESolver::Init(_f);
103  int n = f->Width();
104  y.SetSize(n, mem_type);
105  k.SetSize(n, mem_type);
106  z.SetSize(n, mem_type);
107 }
108 
109 void RK4Solver::Step(Vector &x, double &t, double &dt)
110 {
111  // 0 |
112  // 1/2 | 1/2
113  // 1/2 | 0 1/2
114  // 1 | 0 0 1
115  // -----+-------------------
116  // | 1/6 1/3 1/3 1/6
117 
118  f->SetTime(t);
119  f->Mult(x, k); // k1
120  add(x, dt/2, k, y);
121  add(x, dt/6, k, z);
122 
123  f->SetTime(t + dt/2);
124  f->Mult(y, k); // k2
125  add(x, dt/2, k, y);
126  z.Add(dt/3, k);
127 
128  f->Mult(y, k); // k3
129  add(x, dt, k, y);
130  z.Add(dt/3, k);
131 
132  f->SetTime(t + dt);
133  f->Mult(y, k); // k4
134  add(z, dt/6, k, x);
135  t += dt;
136 }
137 
138 ExplicitRKSolver::ExplicitRKSolver(int _s, const double *_a, const double *_b,
139  const double *_c)
140 {
141  s = _s;
142  a = _a;
143  b = _b;
144  c = _c;
145  k = new Vector[s];
146 }
147 
149 {
150  ODESolver::Init(_f);
151  int n = f->Width();
152  y.SetSize(n, mem_type);
153  for (int i = 0; i < s; i++)
154  {
155  k[i].SetSize(n, mem_type);
156  }
157 }
158 
159 void ExplicitRKSolver::Step(Vector &x, double &t, double &dt)
160 {
161  // 0 |
162  // c[0] | a[0]
163  // c[1] | a[1] a[2]
164  // ... | ...
165  // c[s-2] | ... a[s(s-1)/2-1]
166  // --------+---------------------
167  // | b[0] b[1] ... b[s-1]
168 
169  f->SetTime(t);
170  f->Mult(x, k[0]);
171  for (int l = 0, i = 1; i < s; i++)
172  {
173  add(x, a[l++]*dt, k[0], y);
174  for (int j = 1; j < i; j++)
175  {
176  y.Add(a[l++]*dt, k[j]);
177  }
178 
179  f->SetTime(t + c[i-1]*dt);
180  f->Mult(y, k[i]);
181  }
182  for (int i = 0; i < s; i++)
183  {
184  x.Add(b[i]*dt, k[i]);
185  }
186  t += dt;
187 }
188 
190 {
191  delete [] k;
192 }
193 
194 const double RK6Solver::a[] =
195 {
196  .6e-1,
197  .1923996296296296296296296296296296296296e-1,
198  .7669337037037037037037037037037037037037e-1,
199  .35975e-1,
200  0.,
201  .107925,
202  1.318683415233148260919747276431735612861,
203  0.,
204  -5.042058063628562225427761634715637693344,
205  4.220674648395413964508014358283902080483,
206  -41.87259166432751461803757780644346812905,
207  0.,
208  159.4325621631374917700365669070346830453,
209  -122.1192135650100309202516203389242140663,
210  5.531743066200053768252631238332999150076,
211  -54.43015693531650433250642051294142461271,
212  0.,
213  207.0672513650184644273657173866509835987,
214  -158.6108137845899991828742424365058599469,
215  6.991816585950242321992597280791793907096,
216  -.1859723106220323397765171799549294623692e-1,
217  -54.66374178728197680241215648050386959351,
218  0.,
219  207.9528062553893734515824816699834244238,
220  -159.2889574744995071508959805871426654216,
221  7.018743740796944434698170760964252490817,
222  -.1833878590504572306472782005141738268361e-1,
223  -.5119484997882099077875432497245168395840e-3
224 };
225 const double RK6Solver::b[] =
226 {
227  .3438957868357036009278820124728322386520e-1,
228  0.,
229  0.,
230  .2582624555633503404659558098586120858767,
231  .4209371189673537150642551514069801967032,
232  4.405396469669310170148836816197095664891,
233  -176.4831190242986576151740942499002125029,
234  172.3641334014150730294022582711902413315
235 };
236 const double RK6Solver::c[] =
237 {
238  .6e-1,
239  .9593333333333333333333333333333333333333e-1,
240  .1439,
241  .4973,
242  .9725,
243  .9995,
244  1.,
245 };
246 
247 const double RK8Solver::a[] =
248 {
249  .5e-1,
250  -.69931640625e-2,
251  .1135556640625,
252  .399609375e-1,
253  0.,
254  .1198828125,
255  .3613975628004575124052940721184028345129,
256  0.,
257  -1.341524066700492771819987788202715834917,
258  1.370126503900035259414693716084313000404,
259  .490472027972027972027972027972027972028e-1,
260  0.,
261  0.,
262  .2350972042214404739862988335493427143122,
263  .180855592981356728810903963653454488485,
264  .6169289044289044289044289044289044289044e-1,
265  0.,
266  0.,
267  .1123656831464027662262557035130015442303,
268  -.3885046071451366767049048108111244567456e-1,
269  .1979188712522045855379188712522045855379e-1,
270  -1.767630240222326875735597119572145586714,
271  0.,
272  0.,
273  -62.5,
274  -6.061889377376669100821361459659331999758,
275  5.650823198222763138561298030600840174201,
276  65.62169641937623283799566054863063741227,
277  -1.180945066554970799825116282628297957882,
278  0.,
279  0.,
280  -41.50473441114320841606641502701994225874,
281  -4.434438319103725011225169229846100211776,
282  4.260408188586133024812193710744693240761,
283  43.75364022446171584987676829438379303004,
284  .787142548991231068744647504422630755086e-2,
285  -1.281405999441488405459510291182054246266,
286  0.,
287  0.,
288  -45.04713996013986630220754257136007322267,
289  -4.731362069449576477311464265491282810943,
290  4.514967016593807841185851584597240996214,
291  47.44909557172985134869022392235929015114,
292  .1059228297111661135687393955516542875228e-1,
293  -.5746842263844616254432318478286296232021e-2,
294  -1.724470134262485191756709817484481861731,
295  0.,
296  0.,
297  -60.92349008483054016518434619253765246063,
298  -5.95151837622239245520283276706185486829,
299  5.556523730698456235979791650843592496839,
300  63.98301198033305336837536378635995939281,
301  .1464202825041496159275921391759452676003e-1,
302  .6460408772358203603621865144977650714892e-1,
303  -.7930323169008878984024452548693373291447e-1,
304  -3.301622667747079016353994789790983625569,
305  0.,
306  0.,
307  -118.011272359752508566692330395789886851,
308  -10.14142238845611248642783916034510897595,
309  9.139311332232057923544012273556827000619,
310  123.3759428284042683684847180986501894364,
311  4.623244378874580474839807625067630924792,
312  -3.383277738068201923652550971536811240814,
313  4.527592100324618189451265339351129035325,
314  -5.828495485811622963193088019162985703755
315 };
316 const double RK8Solver::b[] =
317 {
318  .4427989419007951074716746668098518862111e-1,
319  0.,
320  0.,
321  0.,
322  0.,
323  .3541049391724448744815552028733568354121,
324  .2479692154956437828667629415370663023884,
325  -15.69420203883808405099207034271191213468,
326  25.08406496555856261343930031237186278518,
327  -31.73836778626027646833156112007297739997,
328  22.93828327398878395231483560344797018313,
329  -.2361324633071542145259900641263517600737
330 };
331 const double RK8Solver::c[] =
332 {
333  .5e-1,
334  .1065625,
335  .15984375,
336  .39,
337  .465,
338  .155,
339  .943,
340  .901802041735856958259707940678372149956,
341  .909,
342  .94,
343  1.,
344 };
345 
346 
348 {
349  ODESolver::Init(_f);
350  k.SetSize(f->Width(), mem_type);
351 }
352 
353 void BackwardEulerSolver::Step(Vector &x, double &t, double &dt)
354 {
355  f->SetTime(t + dt);
356  f->ImplicitSolve(dt, x, k); // solve for k: k = f(x + dt*k, t + dt)
357  x.Add(dt, k);
358  t += dt;
359 }
360 
361 
363 {
364  ODESolver::Init(_f);
365  k.SetSize(f->Width(), mem_type);
366 }
367 
368 void ImplicitMidpointSolver::Step(Vector &x, double &t, double &dt)
369 {
370  f->SetTime(t + dt/2);
371  f->ImplicitSolve(dt/2, x, k);
372  x.Add(dt, k);
373  t += dt;
374 }
375 
376 
378 {
379  if (gamma_opt == 0)
380  {
381  gamma = (3. - sqrt(3.))/6.; // not A-stable, order 3
382  }
383  else if (gamma_opt == 2)
384  {
385  gamma = (2. - sqrt(2.))/2.; // L-stable, order 2
386  }
387  else if (gamma_opt == 3)
388  {
389  gamma = (2. + sqrt(2.))/2.; // L-stable, order 2
390  }
391  else
392  {
393  gamma = (3. + sqrt(3.))/6.; // A-stable, order 3
394  }
395 }
396 
398 {
399  ODESolver::Init(_f);
400  k.SetSize(f->Width(), mem_type);
401  y.SetSize(f->Width(), mem_type);
402 }
403 
404 void SDIRK23Solver::Step(Vector &x, double &t, double &dt)
405 {
406  // with a = gamma:
407  // a | a
408  // 1-a | 1-2a a
409  // ------+-----------
410  // | 1/2 1/2
411  // note: with gamma_opt=3, both solve are outside [t,t+dt] since a>1
412  f->SetTime(t + gamma*dt);
413  f->ImplicitSolve(gamma*dt, x, k);
414  add(x, (1.-2.*gamma)*dt, k, y); // y = x + (1-2*gamma)*dt*k
415  x.Add(dt/2, k);
416 
417  f->SetTime(t + (1.-gamma)*dt);
418  f->ImplicitSolve(gamma*dt, y, k);
419  x.Add(dt/2, k);
420  t += dt;
421 }
422 
423 
425 {
426  ODESolver::Init(_f);
427  k.SetSize(f->Width(), mem_type);
428  y.SetSize(f->Width(), mem_type);
429  z.SetSize(f->Width(), mem_type);
430 }
431 
432 void SDIRK34Solver::Step(Vector &x, double &t, double &dt)
433 {
434  // a | a
435  // 1/2 | 1/2-a a
436  // 1-a | 2a 1-4a a
437  // ------+--------------------
438  // | b 1-2b b
439  // note: two solves are outside [t,t+dt] since c1=a>1, c3=1-a<0
440  const double a = 1./sqrt(3.)*cos(M_PI/18.) + 0.5;
441  const double b = 1./(6.*(2.*a-1.)*(2.*a-1.));
442 
443  f->SetTime(t + a*dt);
444  f->ImplicitSolve(a*dt, x, k);
445  add(x, (0.5-a)*dt, k, y);
446  add(x, (2.*a)*dt, k, z);
447  x.Add(b*dt, k);
448 
449  f->SetTime(t + dt/2);
450  f->ImplicitSolve(a*dt, y, k);
451  z.Add((1.-4.*a)*dt, k);
452  x.Add((1.-2.*b)*dt, k);
453 
454  f->SetTime(t + (1.-a)*dt);
455  f->ImplicitSolve(a*dt, z, k);
456  x.Add(b*dt, k);
457  t += dt;
458 }
459 
460 
462 {
463  ODESolver::Init(_f);
464  k.SetSize(f->Width(), mem_type);
465  y.SetSize(f->Width(), mem_type);
466 }
467 
468 void SDIRK33Solver::Step(Vector &x, double &t, double &dt)
469 {
470  // a | a
471  // c | c-a a
472  // 1 | b 1-a-b a
473  // -----+----------------
474  // | b 1-a-b a
475  const double a = 0.435866521508458999416019;
476  const double b = 1.20849664917601007033648;
477  const double c = 0.717933260754229499708010;
478 
479  f->SetTime(t + a*dt);
480  f->ImplicitSolve(a*dt, x, k);
481  add(x, (c-a)*dt, k, y);
482  x.Add(b*dt, k);
483 
484  f->SetTime(t + c*dt);
485  f->ImplicitSolve(a*dt, y, k);
486  x.Add((1.-a-b)*dt, k);
487 
488  f->SetTime(t + dt);
489  f->ImplicitSolve(a*dt, x, k);
490  x.Add(a*dt, k);
491  t += dt;
492 }
493 
494 
496 {
497  ODESolver::Init(_f);
498  k.SetSize(f->Width(), mem_type);
499  y.SetSize(f->Width(), mem_type);
500  xdot.SetSize(f->Width(), mem_type);
501  xdot = 0.0;
502  first = true;
503 }
504 
506 {
507  rho_inf = (rho_inf > 1.0) ? 1.0 : rho_inf;
508  rho_inf = (rho_inf < 0.0) ? 0.0 : rho_inf;
509 
510  alpha_m = 0.5*(3.0 - rho_inf)/(1.0 + rho_inf);
511  alpha_f = 1.0/(1.0 + rho_inf);
512  gamma = 0.5 + alpha_m - alpha_f;
513 }
514 
516 {
517  out << "Generalized alpha time integrator:" << std::endl;
518  out << "alpha_m = " << alpha_m << std::endl;
519  out << "alpha_f = " << alpha_f << std::endl;
520  out << "gamma = " << gamma << std::endl;
521 
522  if (gamma == 0.5 + alpha_m - alpha_f)
523  {
524  out<<"Second order"<<" and ";
525  }
526  else
527  {
528  out<<"First order"<<" and ";
529  }
530 
531  if ((alpha_m >= alpha_f)&&(alpha_f >= 0.5))
532  {
533  out<<"Stable"<<std::endl;
534  }
535  else
536  {
537  out<<"Unstable"<<std::endl;
538  }
539 }
540 
541 // This routine assumes xdot is initialized.
542 void GeneralizedAlphaSolver::Step(Vector &x, double &t, double &dt)
543 {
544  double dt_fac1 = alpha_f*(1.0 - gamma/alpha_m);
545  double dt_fac2 = alpha_f*gamma/alpha_m;
546  double dt_fac3 = 1.0/alpha_m;
547 
548  // In the first pass xdot is not yet computed. If parameter choices requires
549  // xdot midpoint rule is used instead for the first step only.
550  if (first && (dt_fac1 != 0.0))
551  {
552  dt_fac1 = 0.0;
553  dt_fac2 = 0.5;
554  dt_fac3 = 2.0;
555  first = false;
556  }
557 
558  add(x, dt_fac1*dt, xdot, y);
559  f->SetTime(t + dt_fac2*dt);
560  f->ImplicitSolve(dt_fac2*dt, y, k);
561 
562  add(y, dt_fac2*dt, k, x);
563  k.Add(-1.0, xdot);
564  xdot.Add(dt_fac3, k);
565 
566  t += dt;
567 }
568 
569 
570 void
572 {
573  P_ = &P; F_ = &F;
574 
575  dp_.SetSize(F_->Height());
576  dq_.SetSize(P_->Height());
577 }
578 
579 void
580 SIA1Solver::Step(Vector &q, Vector &p, double &t, double &dt)
581 {
582  F_->SetTime(t);
583  F_->Mult(q,dp_);
584  p.Add(dt,dp_);
585 
586  P_->Mult(p,dq_);
587  q.Add(dt,dq_);
588 
589  t += dt;
590 }
591 
592 void
593 SIA2Solver::Step(Vector &q, Vector &p, double &t, double &dt)
594 {
595  P_->Mult(p,dq_);
596  q.Add(0.5*dt,dq_);
597 
598  F_->SetTime(t+0.5*dt);
599  F_->Mult(q,dp_);
600  p.Add(dt,dp_);
601 
602  P_->Mult(p,dq_);
603  q.Add(0.5*dt,dq_);
604 
605  t += dt;
606 }
607 
609  : order_(order)
610 {
611  a_.SetSize(order);
612  b_.SetSize(order);
613 
614  switch (order_)
615  {
616  case 1:
617  a_[0] = 1.0;
618  b_[0] = 1.0;
619  break;
620  case 2:
621  a_[0] = 0.5;
622  a_[1] = 0.5;
623  b_[0] = 0.0;
624  b_[1] = 1.0;
625  break;
626  case 3:
627  a_[0] = 2.0/3.0;
628  a_[1] = -2.0/3.0;
629  a_[2] = 1.0;
630  b_[0] = 7.0/24.0;
631  b_[1] = 0.75;
632  b_[2] = -1.0/24.0;
633  break;
634  case 4:
635  a_[0] = (2.0+pow(2.0,1.0/3.0)+pow(2.0,-1.0/3.0))/6.0;
636  a_[1] = (1.0-pow(2.0,1.0/3.0)-pow(2.0,-1.0/3.0))/6.0;
637  a_[2] = a_[1];
638  a_[3] = a_[0];
639  b_[0] = 0.0;
640  b_[1] = 1.0/(2.0-pow(2.0,1.0/3.0));
641  b_[2] = 1.0/(1.0-pow(2.0,2.0/3.0));
642  b_[3] = b_[1];
643  break;
644  default:
645  MFEM_ASSERT(false, "Unsupported order in SIAVSolver");
646  };
647 }
648 
649 void
650 SIAVSolver::Step(Vector &q, Vector &p, double &t, double &dt)
651 {
652  for (int i=0; i<order_; i++)
653  {
654  if ( b_[i] != 0.0 )
655  {
656  F_->SetTime(t);
657  if ( F_->isExplicit() )
658  {
659  F_->Mult(q, dp_);
660  }
661  else
662  {
663  F_->ImplicitSolve(b_[i] * dt, q, dp_);
664  }
665  p.Add(b_[i] * dt, dp_);
666  }
667 
668  P_->Mult(p, dq_);
669  q.Add(a_[i] * dt, dq_);
670 
671  t += a_[i] * dt;
672  }
673 }
674 
675 }
virtual void Init(TimeDependentOperator &_f)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:39
void PrintProperties(std::ostream &out=mfem::out)
Definition: ode.cpp:515
virtual void Init(Operator &P, TimeDependentOperator &F)
Definition: ode.cpp:571
void SetRhoInf(double rho_inf)
Definition: ode.cpp:505
void Step(Vector &q, Vector &p, double &t, double &dt)
Definition: ode.cpp:580
virtual void Init(TimeDependentOperator &_f)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:397
virtual void Step(Vector &x, double &t, double &dt)
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in]...
Definition: ode.cpp:109
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:400
Base abstract class for time dependent operators.
Definition: operator.hpp:162
MemoryType mem_type
Definition: ode.hpp:27
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:42
virtual void Init(TimeDependentOperator &f)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:18
virtual void SetTime(const double _t)
Set the current time.
Definition: operator.hpp:192
virtual void Step(Vector &x, double &t, double &dt)
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in]...
Definition: ode.cpp:159
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
virtual void Step(Vector &x, double &t, double &dt)
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in]...
Definition: ode.cpp:468
void add(const Vector &v1, const Vector &v2, Vector &v)
Definition: vector.cpp:238
virtual void Step(Vector &x, double &t, double &dt)
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in]...
Definition: ode.cpp:432
virtual void Init(TimeDependentOperator &_f)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:362
virtual void Init(TimeDependentOperator &_f)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:100
virtual MemoryClass GetMemoryClass() const
Return the MemoryClass preferred by the Operator.
Definition: operator.hpp:56
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:36
virtual ~ExplicitRKSolver()
Definition: ode.cpp:189
virtual void Step(Vector &x, double &t, double &dt)
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in]...
Definition: ode.cpp:47
Vector dq_
Definition: ode.hpp:342
virtual void Init(TimeDependentOperator &_f)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:68
bool isExplicit() const
True if type is EXPLICIT.
Definition: operator.hpp:195
virtual void ImplicitSolve(const double dt, const Vector &x, Vector &k)
Solve the equation: k = f(x + dt k, t), for the unknown k at the current time t.
Definition: operator.hpp:245
virtual void Step(Vector &x, double &t, double &dt)
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in]...
Definition: ode.cpp:368
ExplicitRKSolver(int _s, const double *_a, const double *_b, const double *_c)
Definition: ode.cpp:138
virtual void Step(Vector &x, double &t, double &dt)
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in]...
Definition: ode.cpp:353
virtual void Init(TimeDependentOperator &_f)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:347
virtual void Init(TimeDependentOperator &_f)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:495
virtual void Init(TimeDependentOperator &_f)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:424
virtual void Mult(const Vector &x, Vector &y) const
Perform the action of the operator: y = k = f(x, t), where k solves the algebraic equation F(x...
Definition: operator.hpp:224
virtual void Init(TimeDependentOperator &_f)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:24
TimeDependentOperator * F_
Definition: ode.hpp:338
void SetSize(int nsize)
Change logical size of the array, keep existing entries.
Definition: array.hpp:618
Vector dp_
Definition: ode.hpp:341
Operator * P_
Definition: ode.hpp:339
Vector & Add(const double a, const Vector &Va)
(*this) += a * Va
Definition: vector.cpp:190
void Step(Vector &q, Vector &p, double &t, double &dt)
Definition: ode.cpp:593
void Step(Vector &q, Vector &p, double &t, double &dt)
Definition: ode.cpp:650
MemoryType GetMemoryType(MemoryClass mc)
Return a suitable MemoryType for a given MemoryClass.
Definition: mem_manager.cpp:23
SIAVSolver(int order)
Definition: ode.cpp:608
virtual void Step(Vector &x, double &t, double &dt)
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in]...
Definition: ode.cpp:542
virtual void Step(Vector &x, double &t, double &dt)
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in]...
Definition: ode.cpp:76
Vector data type.
Definition: vector.hpp:48
virtual void Init(TimeDependentOperator &_f)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:461
virtual void Step(Vector &x, double &t, double &dt)
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in]...
Definition: ode.cpp:404
TimeDependentOperator * f
Pointer to the associated TimeDependentOperator.
Definition: ode.hpp:26
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
Definition: globals.hpp:64
Abstract operator.
Definition: operator.hpp:21
virtual void Step(Vector &x, double &t, double &dt)
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in]...
Definition: ode.cpp:30
SDIRK23Solver(int gamma_opt=1)
Definition: ode.cpp:377
virtual void Init(TimeDependentOperator &_f)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:148