← Back to branch summary

~ubuntu-branches/ubuntu/maverick/dolfin/maverick

« back to all changes in this revision

Viewing changes to demo/ode/reaction/cpp/main.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Johannes Ring
  • Date: 2009-10-12 14:13:18 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20091012141318-hkbxl0sq555vqv7d
Tags: 0.9.4-1
* New upstream release. This version cleans up the design of the
  function class by adding a new abstraction for user-defined
  functions called Expression. A number of minor bugfixes and
  improvements have also been made.
* debian/watch: Update for new flat directory structure.
* Update debian/copyright.
* debian/rules: Use explicit paths to PETSc 3.0.0 and SLEPc 3.0.0.

Show diffs side-by-side

added added

removed removed

Lines of Context:
demo/ode/reaction/cpp/main.cpp
 
1
// Copyright (C) 2005-2009 Anders Logg.
 
2
// Licensed under the GNU LGPL Version 2.1.
 
3
//
 
4
// First added:  2005-10-14
 
5
// Last changed: 2009-09-08
 
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
/// This is a simplified (only multi-adaptive) version of the ODE
 
19
/// benchmark found under bench/ode/reaction
 
20
 
 
21
class Reaction : public ODE
 
22
{
 
23
public:
 
24
 
 
25
  /// Constructor
 
26
  Reaction(unsigned int N, real T, real L, real epsilon, real gamma)
 
27
    : ODE(N, T), L(L), epsilon(epsilon), gamma(gamma)
 
28
  {
 
29
    // Compute parameters
 
30
    h = L / static_cast<real>(N - 1);
 
31
    lambda = 0.5*sqrt(2.0*gamma/epsilon);
 
32
    v = 0.5*sqrt(2.0*gamma*epsilon);
 
33
 
 
34
    // Set sparse dependency pattern
 
35
    for (unsigned int i = 0; i < N; i++)
 
36
    {
 
37
      dependencies.setsize(i, 3);
 
38
      if (i > 0) dependencies.set(i, i - 1);
 
39
      dependencies.set(i, i);
 
40
      if (i < (N - 1)) dependencies.set(i, i + 1);
 
41
    }
 
42
  }
 
43
 
 
44
  /// Initial condition
 
45
  void u0(real* u)
 
46
  {
 
47
    for (unsigned int i = 0; i < N; i++)
 
48
    {
 
49
      const real x = static_cast<real>(i)*h;
 
50
      u[i] = 1.0 / (1.0 + real_exp(lambda*(x - 1.0)));
 
51
    }
 
52
  }
 
53
 
 
54
  /// Right-hand side, multi-adaptive version
 
55
  real f(const real* u, real t, unsigned int 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
    return epsilon * sum / (h*h) + gamma * ui*ui * (1.0 - ui);
 
68
  }
 
69
 
 
70
public:
 
71
 
 
72
  real L;       // Length of domain
 
73
  real epsilon; // Diffusivity
 
74
  real gamma;   // Reaction rate
 
75
  real h;       // Mesh size
 
76
  real lambda;  // Parameter for initial data
 
77
  real v;       // Speed of reaction front
 
78
 
 
79
};
 
80
 
 
81
int main()
 
82
{
 
83
  info("Reaction ODE demo needs to be fixed.");
 
84
 
 
85
  // FIXME: Does not work with GMP enabled
 
86
 
 
87
  // Set some parameters
 
88
  const real T = 0.01;
 
89
  const real epsilon = 0.01;
 
90
  const real gamma = 1000.0;
 
91
  const real L = 1.0;
 
92
  const unsigned int N = 5000;
 
93
 
 
94
  // Create ODE
 
95
  Reaction ode(N, T, L, epsilon, gamma);
 
96
  ode.parameters["method"] = "mcg";
 
97
  ode.parameters["order"] = 1;
 
98
  ode.parameters["nonlinear_solver"] = "fixed-point";
 
99
  ode.parameters["tolerance"] = 1e-3;
 
100
  ode.parameters["partitioning_threshold"] = 0.5;
 
101
  ode.parameters["initial_time_step"] = 1e-5;
 
102
  ode.parameters["maximum_time_step"] = 1e-3;
 
103
  ode.parameters["adaptive_samples"] = true;
 
104
 
 
105
  // Solve ODE
 
106
  ode.solve();
 
107
 
 
108
  return 0;
 
109
}