← Back to branch summary

~njansson/dolfin/hpc

« back to all changes in this revision

Viewing changes to bench/ode/reaction/main.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:
bench/ode/reaction/main.cpp
 
1
// Copyright (C) 2005-2006 Anders Logg.
 
2
// Licensed under the GNU LGPL Version 2.1.
 
3
//
 
4
// First added:  2005-10-14
 
5
// Last changed: 2006-08-22
 
6
 
 
7
#include <dolfin.h>
 
8
 
 
9
using namespace dolfin;
 
10
 
 
11
/// Test problem taken from "Multirate time stepping for parabolic PDEs"
 
12
/// by Savcenco, Hundsdorfer and Verwer:
 
13
///
 
14
///    u' - epsilon u'' = gamma u^2 (1 - u)  in (0, L) x (0, T).
 
15
///
 
16
/// The solution is a reaction front sweeping across the domain.
 
17
 
 
18
class Reaction : public ODE
 
19
{
 
20
public:
 
21
 
 
22
  /// Constructor
 
23
  Reaction(unsigned int N, real T, real L, real epsilon, real gamma)
 
24
    : ODE(N, T), L(L), epsilon(epsilon), gamma(gamma)
 
25
  {
 
26
    // Compute parameters
 
27
    h = L / static_cast<real>(N - 1);
 
28
    lambda = 0.5*sqrt(2.0*gamma/epsilon);
 
29
    v = 0.5*sqrt(2.0*gamma*epsilon);
 
30
    
 
31
    // Set sparse dependency pattern
 
32
    for (unsigned int i = 0; i < N; i++)
 
33
    {
 
34
      dependencies.setsize(i, 3);
 
35
      if ( i > 0 ) dependencies.set(i, i - 1);
 
36
      dependencies.set(i, i);
 
37
      if ( i < (N - 1) ) dependencies.set(i, i + 1);
 
38
    }
 
39
  }
 
40
 
 
41
  /// Initial condition
 
42
  void u0(uBlasVector& u)
 
43
  {
 
44
    for (unsigned int i = 0; i < N; i++)
 
45
    {
 
46
      const real x = static_cast<real>(i)*h;
 
47
      u(i) = 1.0 / (1.0 + exp(lambda*(x - 1.0)));
 
48
    }
 
49
  }
 
50
 
 
51
  /// Right-hand side, mono-adaptive version
 
52
  void f(const uBlasVector& u, real t, uBlasVector& y)
 
53
  {
 
54
    printf("MONO!\n");
 
55
    for (unsigned int i = 0; i < N; i++)
 
56
    {
 
57
      const real ui = u(i);
 
58
 
 
59
      real sum = 0.0;
 
60
      if ( i == 0 )
 
61
        sum = u(i + 1) - ui;
 
62
      else if ( i == (N - 1) )
 
63
        sum = u(i - 1) - ui;
 
64
      else
 
65
        sum = u(i + 1) - 2.0*ui + u(i - 1);
 
66
 
 
67
      y(i) = epsilon * sum / (h*h) + gamma * ui*ui * (1.0 - ui);
 
68
    }
 
69
  }
 
70
 
 
71
  /// Right-hand side, multi-adaptive version
 
72
  real f(const uBlasVector& u, real t, unsigned int i)
 
73
  {
 
74
    printf("MULTI!\n");
 
75
    const real ui = u(i);
 
76
    
 
77
    real sum = 0.0;
 
78
    if ( i == 0 )
 
79
      sum = u(i + 1) - ui;
 
80
    else if ( i == (N - 1) )
 
81
      sum = u(i - 1) - ui;
 
82
    else
 
83
      sum = u(i + 1) - 2.0*ui + u(i - 1);
 
84
    
 
85
    return epsilon * sum / (h*h) + gamma * ui*ui * (1.0 - ui);
 
86
  }
 
87
 
 
88
public:
 
89
 
 
90
  real L;       // Length of domain
 
91
  real epsilon; // Diffusivity
 
92
  real gamma;   // Reaction rate
 
93
  real h;       // Mesh size
 
94
  real lambda;  // Parameter for initial data
 
95
  real v;       // Speed of reaction front
 
96
 
 
97
};
 
98
 
 
99
int main(int argc, char* argv[])
 
100
{
 
101
  // Parse command line arguments
 
102
  if ( argc != 7 )
 
103
  {
 
104
    message("Usage: dolfin-ode-reaction method solver tol kmax N L params");
 
105
    message("");
 
106
    message("method - 'cg' or 'mcg'");
 
107
    message("solver - 'fixed-point' or 'newton'");
 
108
    message("tol    - tolerance");
 
109
    message("N      - number of components");
 
110
    message("L      - length of domain");
 
111
    message("params - name of parameter file");
 
112
    return 1;
 
113
  }
 
114
  const char* method = argv[1];
 
115
  const char* solver = argv[2];
 
116
  const real tol = static_cast<real>(atof(argv[3]));
 
117
  const unsigned int N = static_cast<unsigned int>(atoi(argv[4]));
 
118
  const real L = static_cast<unsigned int>(atof(argv[5]));
 
119
  const char* params = argv[6];
 
120
  
 
121
  // Load solver parameters from file
 
122
  File file(params);
 
123
  file >> ParameterSystem::parameters;
 
124
 
 
125
  // Set remaining solver parameters from command-line arguments
 
126
  set("ODE method", method);
 
127
  set("ODE nonlinear solver", solver);
 
128
  set("ODE tolerance", tol);
 
129
 
 
130
  // Set fixed parameters for test problem
 
131
  const real T = 1.0;
 
132
  const real epsilon = 0.01;
 
133
  const real gamma = 1000.0;
 
134
 
 
135
  // Solve system of ODEs
 
136
  Reaction ode(N, T, L, epsilon, gamma);
 
137
  ode.solve();
 
138
 
 
139
  return 0;
 
140
}