Title  Object Oriented Numerics
Participant  Gerhard Zumbusch
Key words  Object Oriented Design, Problem Solving Environments, Software Engineering
Description  The goal for the development of our finite difference and sparse grid codes is to be able to test and to verify different types of discretizations on sparse grids and to tackle different types of partial differential equations. Hence a flexible and modular design is a must. Techniques such as abstract data types and object oriented programming provide such a flexibility. However, they often lead to slow inefficient code due to an over-use of design features.

For example, overloading the arithmetic operators `+' and `*' for vectors in C++ usually is less efficient than proving directly a saxpy operation for expressions of   av + w  type. Performing tree traversal for the operations above is even slower. Hence, splitting a large code into many small functions and loops may lead to inefficiencies, which cannot be resolved by an optimizing compiler. numerical vector

We have based our code on several fundamental abstractions, which are well separated both in functionality and in implementation. This guarantees that we do not loose efficiency in a substantial way due to this separation. We have identified the following building blocks.

They are ordered from low level, computationally expensive and efficient, to high level routines, where efficiency is achieved through call of some efficient subroutines. Similar abstractions can be found in other object oriented software packages for partial differetial equations such as Diffpack.
• grid
• vector
• field
• operator
• solver

grid

geometric description of an adaptively refined sparse grid and provides i/o, refinement and addressing, or e.g. a full grid for debugging purposes

vector

a large container for real numbers, including BLAS level 1

data: Fortran compatible linear storage

class Vec {
protected:
real *r;
int n;
...
};

methods: Blas level 1 (with compiler inlining)

void Vec:: add(const Vec &a, const Vec &b)
{
assert((a.n==n)&&(b.n==n));
for (int i=0; i<n; i++)
r[i] = a.r[i] + b.r[i];
}

including bit (flag) vectors and short vectors for coordinates etc.

field

a (solution) function on a  grid. It provides a mapping between a grid and a vector and interprets the data as (collocated) scalar or vector field
including special purpose flag fields ect.

operator

the finite difference operators and operates on grids
and non-linear operators

solver

different iterative (Krylov) solvers, uses a (differential) operator: iterative solvers for linear problems,

linear solver
data: operator, parameters
methods: solve

void ConjGrad:: solve(const Vec &b, Vec &x)
{
assert(b.dim()==x.dim());
int n = x.dim();
Vec r(n), p(n), Ap(n);
iter = 0;
real beta, alpha, rtr, rtrold;
op->apply(x, Ap);
p.sub(b, Ap);
r.copy(p);
rtr = r.prod(r);

while (rtr>tol) {
op->apply(p, Ap);
alpha = rtr / p.prod(Ap);
rtrold = rtr;
rtr = r.prod(r);
beta = rtr / rtrold;
iter++;
printIter("ConjGrad it ", " res= ", rtr);
if ((iter>=maxiter)||(error>=HUGE)) break;
}
error = rtr;
}

for eigen-value problems,

and for non-linear problems

main

putting the components together leads to a simple main program, which may look like this:

calling sequence main

uchar dim = 2;
SVec x0(dim), x1(dim);
x0.set(-1);
x1.set(1);

SpGrid grid(dim, x0, x1);         //create sparse grid
grid.refineAll();                           //refine regularly

uint n = grid.nodes();
Vec x(n), f(n);
x.set(0);
f.set(1);
SpField field(x, grid);
SpLaplace lap(field, 1);     //create Poisson problem

Bicgstab j;
j.attach(lap);
j.setTol(1e-8);
j.setMaxIter(10000);
j.setVerbose(1);
j.solve(f, x);                                // solve equation system

cout<<j<<endl;
ofstream of("lap.wrl");
field.print(of, Grid::vrml); // dump solution

For the implementation of such finite difference methods, we have discussed some software concepts. This included an object oriented layout of the basic ingredients of such a code for a powerfull, highly flexible and modular and still efficient code development.

Bibliography
• E. Arge, A. M. Bruaset, H. P. Langtangen, Object-oriented Numerics in Numerical Methods and Software Tools in Industrial Mathematics, eds.: M. Dæhlen, A. Tveito, Birkhäuser, Basel, 1997.
• A. M. Bruaset, H. P. Langtangen, G. Zumbusch, Domain Decomposition and Multilevel Methods in Diffpack. Proceedings of Domain Decomposition Methods 9, Wiley, 1998, pp. 655-662
• T. Schiekofer, G. Zumbusch, Software Concepts of a Sparse Grid Finite Difference Code. , Proceedings of the 14th GAMM-Seminar Kiel on Concepts of Numerical Software, editor: W. Hackbusch, G. Wittum, Notes on Numerical Fluid Mechanics, Vieweg, 1998
• G. Zumbusch, A Sparse Grid PDE Solver. , Advances in Software Tools for Scientific Computing (Proceedings SciTools '98), editor: H. P. Langtangen, A. M. Bruaset, E. Quak, Springer, 2000, chapter 4 in Lecture Notes in Computational Science and Engineering 10, pp. 133-178
• Related projects
• Adaptive Parallel Multigrid Methods with Space-Filling Curves and Hash Storage
• Finite Difference Schemes on Sparse Grids for Time Dependent Problems