← Back to branch summary

~njansson/dolfin/hpc

« back to all changes in this revision

Viewing changes to demo/ode/stiff/TestProblem5.h

  • Committer: Garth N. Wells
  • Date: 2008-03-29 09:34:25 UTC
  • Revision ID: gnw20@cam.ac.uk-20080329093425-3ea2vhjvccq56zvi
Add some basic build & install instructions.

Show diffs side-by-side

added added

removed removed

Lines of Context:
demo/ode/stiff/TestProblem5.h
 
1
// Copyright (C) 2004-2006 Anders Logg.
 
2
// Licensed under the GNU LGPL Version 2.1.
 
3
//
 
4
// First added:  2004
 
5
// Last changed: 2006-08-21
 
6
 
 
7
#include <cmath>
 
8
#include <dolfin.h>
 
9
 
 
10
using namespace dolfin;
 
11
 
 
12
class TestProblem5 : public ODE
 
13
{
 
14
public:
 
15
  
 
16
  TestProblem5() : ODE(6, 180.0)
 
17
  {
 
18
    message("The Chemical Akzo-Nobel problem.");
 
19
 
 
20
    k1  = 18.7;
 
21
    k2  = 0.58;
 
22
    k3  = 0.09;
 
23
    k4  = 0.42;
 
24
    K   = 34.4;
 
25
    klA = 3.3;
 
26
    Ks  = 115.83;
 
27
    p   = 0.9;
 
28
    H   = 737.0;
 
29
  }
 
30
 
 
31
  void u0(uBlasVector& u)
 
32
  {
 
33
    u(0) = 0.444;
 
34
    u(1) = 0.00123;
 
35
    u(2) = 0.0;
 
36
    u(3) = 0.007;
 
37
    u(4) = 0.0;
 
38
    u(5) = 0.36;
 
39
  }
 
40
 
 
41
  void f(const uBlasVector& u, real t, uBlasVector& y)
 
42
  {
 
43
    y(0) = -2.0*r1(u) + r2(u) - r3(u) - r4(u);
 
44
    y(1) = -0.5*r1(u) - r4(u) - 0.5*r5(u) + F(u);
 
45
    y(2) = r1(u) - r2(u) + r3(u);
 
46
    y(3) = -r2(u) + r3(u) - 2.0*r4(u);
 
47
    y(4) = r2(u) - r3(u) + r5(u);
 
48
    y(5) = Ks*u(0)*u(3) - u(5);
 
49
  }
 
50
 
 
51
private:
 
52
 
 
53
  real r1(const uBlasVector& u)
 
54
  {
 
55
    return k1*pow(u(0), 4.0)*sqrt(u(1));
 
56
  }
 
57
  
 
58
  real r2(const uBlasVector& u)
 
59
  {
 
60
    return k2*u(2)*u(3);
 
61
  }
 
62
 
 
63
  real r3(const uBlasVector& u)
 
64
  {
 
65
    return (k2/K)*u(0)*u(4);
 
66
  }
 
67
 
 
68
  real r4(const uBlasVector& u)
 
69
  {
 
70
    return k3*u(0)*pow(u(3), 2.0);
 
71
  }
 
72
 
 
73
  real r5(const uBlasVector& u)
 
74
  {
 
75
    return k4*pow(u(5), 2.0)*sqrt(u(1));
 
76
  }
 
77
 
 
78
  real F(const uBlasVector& u)
 
79
  {
 
80
    return klA * (p/H - u(1));
 
81
  }
 
82
 
 
83
  real k1, k2, k3, k4, K, klA, Ks, p, H;
 
84
 
 
85
};