← Back to branch summary

~njansson/dolfin/hpc

« back to all changes in this revision

Viewing changes to sandbox/la/hypre/main.cpp

  • 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:
sandbox/la/hypre/main.cpp
1
 
// Copyright (C) 2008 Anders Logg.
2
 
// Licensed under the GNU LGPL Version 2.1.
3
 
//
4
 
// Simple test program for solving a singular system
5
 
// with GMRES + AMG (PETSc + Hypre BoomerAMG)
6
 
//
7
 
// This system has a zero eigenvalue with corresponding
8
 
// eigenvector v = (1, 1, 1, ...). This gives a compatibility
9
 
// condition (corresponding to the integral of the right-hand
10
 
// side being zero for the Neumann problem).
11
 
//
12
 
// To solve the linear system, we must therefore either make the
13
 
// system nonsingular by adding a constraint, like zero average for x,
14
 
// or modify the right-hand side to satisfy the compatibility
15
 
// condition.
16
 
 
17
 
#include <dolfin.h>
18
 
#include "Poisson.h"
19
 
  
20
 
using namespace dolfin;
21
 
 
22
 
int main(int argc, char* argv[])
23
 
{
24
 
  // Initialize PETSc with command-line arguments
25
 
  SubSystemsManager::initPETSc(argc, argv);
26
 
 
27
 
  // Monitor convergence
28
 
  dolfin_set("Krylov monitor convergence", true);
29
 
 
30
 
  // Source term
31
 
  class Source : public Function
32
 
  {
33
 
  public:
34
 
    
35
 
    Source(Mesh& mesh) : Function(mesh) {}
36
 
 
37
 
    real eval(const real* x) const
38
 
    {
39
 
      real dx = x[0] - 0.5;
40
 
      real dy = x[1] - 0.5;
41
 
      return 500.0*exp(-(dx*dx + dy*dy)/0.02);
42
 
    }
43
 
 
44
 
  };
45
 
 
46
 
  // Create mesh
47
 
  UnitSquare mesh(2, 2);
48
 
 
49
 
  // Create functions
50
 
  Source f(mesh);
51
 
 
52
 
  // Variational forms
53
 
  PoissonBilinearForm a;
54
 
  PoissonLinearForm L(f);
55
 
 
56
 
  // Assemble linear system
57
 
  Matrix A;
58
 
  Vector x, b;
59
 
  assemble(A, a, mesh);
60
 
  assemble(b, L, mesh);
61
 
 
62
 
  // Solve linear system using ordinary linear solver
63
 
  LinearSolver s0(gmres, amg);
64
 
  s0.solve(A, x, b);
65
 
  cout << "Residual: " << residual(A, x, b) << endl;
66
 
 
67
 
  // Solve linear system using special singular solver
68
 
  SingularSolver s1(gmres, amg);
69
 
  s1.solve(A, x, b);
70
 
  cout << "Residual: " << residual(A, x, b) << endl;
71
 
 
72
 
  // Solve modified linear system using ordinary linear solver
73
 
  LinearSolver s2(gmres, amg);
74
 
  normalize(b, average);
75
 
  s2.solve(A, x, b);
76
 
  cout << "Residual: " << residual(A, x, b) << endl;
77
 
 
78
 
  return 0;
79
 
}