← Back to branch summary

~njansson/dolfin/hpc

« back to all changes in this revision

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

  • Committer: Anders Logg
  • Date: 2008-05-15 16:02:44 UTC
  • Revision ID: logg@simula.no-20080515160244-e2a4wkfxfzmts16s
Add new class SingularSolver:

This class provides a linear solver for singular linear systems
Ax = b where A has a one-dimensional null-space (kernel). This
may happen for example when solving Poisson's equation with
pure Neumann boundary conditions.

The solver attempts to create an extended non-singular system
by adding the constraint [1, 1, 1, ...]^T x = 0.

If an optional mass matrix M is supplied, the solver attempts
to create an extended non-singular system by adding the
constraint m^T x = 0 where m is the lumped mass matrix. This
corresponds to setting the average (integral) of the finite
element function with coefficients x to zero.

The solver makes not attempt to check that the null-space is
indeed one-dimensional. It is also assumed that the system
Ax = b retains its sparsity pattern between calls to solve().

Show diffs side-by-side

added added

removed removed

Lines of Context:
sandbox/la/hypre/main.cpp
1
 
// Copyright (C) 2006-2007 Anders Logg.
 
1
// Copyright (C) 2008 Anders Logg.
2
2
// Licensed under the GNU LGPL Version 2.1.
3
3
//
4
4
// Simple test program for solving a singular system
5
5
// with GMRES + AMG (PETSc + Hypre BoomerAMG)
 
6
//
 
7
// Singular solver converges much better
6
8
 
7
9
#include <dolfin.h>
8
10
#include "Poisson.h"
34
36
  };
35
37
 
36
38
  // Create mesh
37
 
  UnitSquare mesh(32, 32);
 
39
  UnitSquare mesh(2, 2);
38
40
 
39
41
  // Create functions
40
42
  Source f(mesh);
49
51
  assemble(A, a, mesh);
50
52
  assemble(b, L, mesh);
51
53
 
52
 
  // Solve linear sysem
53
 
  solve(A, x, b, gmres, amg);
 
54
  // Solve linear system using ordinary linear solver
 
55
  LinearSolver s0(gmres, amg);
 
56
  s0.solve(A, x, b);
 
57
 
 
58
  // Solve linear system using special singular solver
 
59
  SingularSolver s1(gmres, amg);
 
60
  s1.solve(A, x, b);
54
61
 
55
62
  return 0;
56
63
}