Research Group of Prof. Dr. M. Griebel
Institute for Numerical Simulation
maximize
 
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); 
    x.add(x,  alpha, p); 
    r.add(r, -alpha, Ap); 
    rtrold = rtr; 
    rtr = r.prod(r); 
    beta = rtr / rtrold; 
    p.add(r, beta, p); 
    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
  • Parallel Adaptive Sparse Grids
    		•
    In cooperation with 
  • SFB 256 "Nonlinear Partial Differential Equations" 
    		•