|
|
SG 2 is a library for the numerical solution of partial differential equations on regular and sparse grids. For the sparse grids, the combination technique is used. SG 2 is intended to solve discretized systems of PDEs of arbitrary dimensionalities. The application area of SG 2 are in particular the numerical simulations in the financial mathematics and population dynamics, as well as computations in the medical imaging methods etc. For the combination technique, SG 2 supports parallel computations. The project page can be found here
SG 2 is a template, class and function library written in C++. It implements the discretization methods, the solvers for the discretized non-linear and linear systems as well as some tools for the output of the results. SG 2 provides its own classes for the representation of the numerical solutions and intermediate results in the memory. These are templates depending on the basic number type T that can be double or float (or in some cases other types). The class hierarchy of the discretization methods and the solvers supports the simple choice of subordinated coupled algorithms (like the linear solver in the Newton method or the preconditioner in the CG-iteration). Most of these classes are capable of reading of parameters from user-provided ini-files that allows the quick fine tuning of the methods.
Application of SG 2 requires some basic knowledge in the field of numerical methods for the partial differential equations. Nevertheless, the class hierarchy of the package allows you to use the methods not going into the details of the implementation. Furthermore, you can implement your own methods (for ex. a new discretization or a preconditioner for the CG-iteration) and use them together with the existing ones.
Examples of the use of the package are included into the distribution (cf. subdirectory Examples). They implement the numerical solution of PDEs of different types. You may use these examples as templates for your applications.
Up to now, SG 2 can be compiled and installed under Unix-type systems including MacOS X. You need a C++-compiler like g++. The installation procedure is standard:
When writing your application that uses SG 2, you #include the SG-2-headers for data structures, methods and tools that you use in your modules. These headers are described below for every particular class. The names of the headers usually coinside with the names of the classes they provide: MultiVector.h contains MultiVector<T> etc. If a method uses some data structures then its header #includes the header files for these data structures. Thus, if you #include MGPrecond.h (for the multigrid method), you have no need to #include MultiVector.h explicitly.
You should also specify SG-2-libraries for the linker. These libraries are:
You may link all the libraries or exclude those that you do not need.
SG 2 uses own classes for the representation of regular and sparse grids, grid functions, sparse linear grid operators and non-linear grid operators. These classes can be divided into 5 categories:
Below in this section, we describe a part of the class hierarchy, implementing the representation of data, in detail. For the detailed information of each class, see the Doxygen documentation of the package.
One of the base classes of SG 2 in its class hierarchy is Array<T>. This class implements a (one-dimensional) array of an arbitrary length. This array is automatically dynamically allocated in the memory. The array may be resized or have the zero length. T may be a an arbitrary type or class. Elements of an Array<T>-object can be accessed using operator [] as for the regular C-arrays. Furthermore, these objects can be assigned to each other using operator = that means setting the new size of the destination object and copying the elements. (There is also a variant of operator = that sets the same value of type T to all the entries of the object.)
Multidimensional arrays of arbitrary dimensionality and arbitrary size in each dimension are stored in objects of MultiArray<T>. This class is used to represent generic grid functions, with no assumptions on T. Every object of MultiArray<T> keeps an dynamically allocated array of all the elements as well as an Array<int> of the sizes. The dimensionality of such an object as well as its sizes may be changed any time.
To access elements of a MultiArray<T>, there is the auxiliary class MultiIndex derived from Array<int>. An object of MultiIndex stores a set of indices (i.e. a multi-index). MultiArray<T> implements operator [] with one MultiIndex as argument. This returns a reference (T &) to the element with these indices. The multi-index should keep exactly the same number of indices as the dimensionality of the MultiArray<T>. Furthermore, there is operator () with a MultiIndex argument. This returns a pointer (T *) to a sub-array indexed by the multi-index. Then the multi-index may be shorter as for the operator []. The elements of a MultiArray<T> may also be accessed by one raw index (of type long) in the lexicographic ordering: Use overloaded operator [] for this.
Besides the access to the elements, MultiIndex implements some facilities to organize loops over all the elements of one MultiArray<T>. The typical problem of implementing such the loops in the cases of arbitrary dimensionality is that you can not merely use nested cycles. One possible solution is to use recursion. The another one is to use functions MultiArray<T>::NextPos and MultiArray<T>::PrevPos like in this example (gfunc is an object of MultiArray<double>):
MultiIndex I (gfunc.Dimension ());
I.SetZero ();
do
cout << "gfunc [" << I << "] =" << gfunc [I] << endl;
while (gfunc.NextPos (I));
In this example, NextPos moves I to the next element in gfunc. As soon as the last element has been reached, NextPos returns false. Otherwise it returns true. There is also a possibility to loop only selected indices in a multi-index when the other indices are fixed.
Classes Array<T> and MultiArray<T> support index range checking. Note that this verification can be switch on or off by the preprocessor. On default, this checking is off. To turn it on, pass --enable-range to configure. If you run configure with the --enable-debug option, the index range checking is turned on.
Checking of the range of the indices is done at every call of operator [] in Array<T> and MultiArray<T>, as well as in some other functions and classes. This helps a lot when debugging your application. But it makes computations significantly slower. Therefore we suggest to compile the release version of your application without the range checking.
If the index range checking is on, SG 2 throws exceptions on every violation of the indexation. You can catch and process these exceptions in your application. You can also make your debugger to stop at throwing exceptions and watch the program stack to find out which function causes these problems. We consider the SG 2 exception classes in Section 2.5.
SG 2 stores grid functions in MultiArrays. But MultiArray<T> is a generic class allowing any type T. For this, there are no algebraic operations defined for MultiArray<T>: Not all types support arithmetics. MultiVector<T> is a class template derived from MultiArray<T> and equipped with the arithmetical operations introducing the structure of a Euclidean vector spase. These operations include addition, subtraction, multiplication by a scalar in different combinations occuring in usual implementations of numerical algorithms. Throughout, where the arithmetics of grid functions is assumed, one should use MultiVector<T>.
Note that in MultiVector<T>, T may not be an arbitrary type. It should be an arithmetical type like double or float.
Besides these arithmetical operations, MultiVector<T> implements the multi-linear interpolation.
In the implicit time stepping schemes, one deals with linear grid operators represented by large sparse matrices. For the efficiency of implementaion of the numerical algorithms, these matrices must be stored in a special format making use of the sparsity. This format should provide the easy matrix-grid-function multiplication and some arithmetical operations between matrices (like the addition or the multiplication by a scalar).
SG 2 uses the Harwell-Boeing-like format for the sparse matrices. These matrices are represented by class SparseMatrix<T>. The arithmetics is defined between SparseMatrix<T> and MultiVector<T> for the same T. Note that as well as for MultiVector<T>, T should be an arithmetical type.
An object of SparseMatrix<T> is a MultiArray of lines of the matrix, i.e. MultiArray<SparseMatrixLine<T>>. For a SparseMatrix<T> A, expression A [I] returns a reference to the matrix line (the object of SparseMatrixLine<T>) which correspond to the grid position with multiindex I. Objects of SparseMatrixLine<T> store potentially nonzero entries of the matrix lines in a compact form. Users of the SG2-library should think of the objects of class SparseMatrixLine<T> as of black boxes. It is not possible to create objects of SparseMatrixLine<T> anywhare except of the friend class SparseMatrix<T>, but you may declare and use references and pointers to these objects. Public member functions of SparseMatrixLine<T> allow you to get the stored entries: There are NEntries () entries, and for each i, 0 <= i < NEntries (), Data (i) returns the reference to the entry and RefOffset (i) returns the offset (the raw index) of the connected element of the grid function. The latter element of the grid function is multiplied by the value of Data (i) in the matrix-grid-function multiplication. Further public member functions allow to set the sparse matrix line to zero, reset its size etc. Class SparseMatrix<T> also provides some functions for working with the matrix lines.
But the user usually does not need to deal directly with the matrix lines and the class SparseMatrixLine<T>. Class SparseMatrix<T> provides functions and operators for the most important matrix-grid-function and matrix-matrix operations. For assembling sparse matrices in discretization routines, there are functions Connection, SetConnection and AddToConnection that provide access to the entries avoiding operations with the SparseMatrix<T>-objects. In this notation, for a matrix A, a ’connection’ is the matrix entry AIJ, where I and J are multiindices of the grid positions.
When creating a object of SparseMatrix<T>, one should specify a number of the entries in the matrix lines. This number is used to allocate memory in the SparseMatrixLine<T>-objects. If you specify more then you really set, some memory is unused. In you specify less entries then your matrix needs, functions like SetConnection and AddToConnection reallocate the memory automatically so that your program should work any way. But the reallocation of the memory every type you add a new entry can slow the program. Thus, a good estimate of the initial number of the entries is an issue of efficiency.
MultiVector<T> represents a (possibly vector-valued) grid function on one grid. Some discretizations of systems of PDEs deal with several grid functions representing different unknowns on grids of several types (vertex-centered, cell-centered etc). A typical example is the discretization of the Navier-Stokes equations on the staggered grids. MultiVector<T> and SparseMatrix<T> do not support implementations of such discretizations, and SG 2 implements special classes for this. Class GridVec<T> is an Array<MultiVector<T>> and can therefore keep a bundle of different grid functions. Each of these grid functions may be defined on it own grid. An object of GridVec<T> can be considered as a block-vector. The algebraic operations implemented for class MultiVector<T> are generalized for GridVec<T>. Clearly, the GridVec<T>-objects in arguments of these operations should be of the same structure.
Linear operators on the space of grid functions represented by GridVec<T>-objects are objects of class SysMatrix<T>. SysMatrix<T> is a square array of SparseMatrix<T>-objects. The matrices represented by SysMatrix<T> can be considered as block-matrices whose blocks are large sparse matrices. Note that some of these blocks may be absend. For this, the corresponding SparseMatrix<T>-object must be left uninitialized, or the SparseMatrix<T>::Reset () method should be called for it. This object represents a zero matrix.
To access the MultiVector<T>-objects in a GridVec<T> object, use GridVec<T>::operator []. Its single argument is the index of the MultiVector<T>. The blocks in a SysMatrix<T>-object are indexed by two integers. To access these blocks, use SysMatrix::operator ().
GridVec<T> and SysMatrix<T> are more general classes as MultiVector<T> and SparseMatrix<T>, and they are prefered base classes for linear and non-linear functions etc. Objects of GridVec<T> and SysMatrix<T> may also store only one MultiVector<T> and a single SparseMatrix<T> and this is often used to make standard implementations of the solvers to work in the most general cases.
Class SysMatrix<T> has one further important feature: It is designed to store hierarchies of sparse matrices. This feature is used in multigrid methods. After declaring several SysMatrix<T>-objects associated with a sequence of refined grids you may attach one to the other using the SysMatrix<T>::AttachCoarse method. This method makes a SysMatrix<T>-object associated with a finer grid to reference a SysMatrix<T>-object associated with the refined grid. Then you get a list of the SysMatrix<T>-objects whose head corresponds to the finest grid.
Unlike linear operators, there is no way to represent a general non-linear operator by merely storing its coefficients in the memory. But non-linear solvers for general systems of non-linear equations do not need such a representation. Instead, the solver should be able to compute the defect of the system and its Jacobian. This is the way that SG 2 uses to work with such systems.
Class NonLinSys<T> provides an interface for this representation. This is an abstract base class which is used in all SG 2 non-linear solvers. The user should derived his own class from this one and implement at least virtual functions NonLinSys<T>::PrintParam, NonLinSys<T>::Defect and NonLinSys<T>::Jacobian.
A numerical solver for a PDE consists of several parts: a discretization, a time-stepping scheme, non-linear and the linear solvers. For some of them, one uses standard algorithms that are sufficiently robust and efficient for a large class of problems (like the implicit Euler method). The other parts are numerical methods that should be carefully tuned, and there are a lot of cases in which a special choice of some methods is required (often, this concerns the linear solvers). Furthermore, every class of PDEs usually requires its own discretization. Although all these discretizations belong to a small number of classes, the differences in details are essential. Therefore, a flexible library for numerical simulations should provide a possibility to incorporate used-defined routines to the entire solver.
The object-oriented structure of SG 2 allows to include the user-defined methods in the solver in the same way as the standard ones. Every type of numerical methods (like the time-stepping schemes, discretizations for stationary and time-dependent PDEs, non-linear and linear solvers) has its own base class. An implementation of a numerical method should be a class derived from one of these base classes.
An example of the interconnections of objects of such the classes for a solver for a non-linear time-dependent problem is presented in Fig. 2.1. The arrows are labeled with the names of the base classes. For example, the implicit Euler method requires a non-linear solver. Here, the Newton method is used. This method requires a linear solver for the Jacobian. Here, the BiCGStab method is used, preconditioned with the geometric multigrid method with the ILU-smoothers. If the Jacobian were too ill-conditioned, we could easily replace the Newton method by the fixed-point iteration that is derived from the same class. Alternatively, we could try other smoothers, and even implement our own smoothers. We could also implement some algebraic multigrid method by deriving a new class from SysPrecond<T>.
The methods are plugged into each other at declaring the objects. To use the CG-iteration preconditioned with the Gustafsson (ILU-1-) preconditioner, you need to declare 2 objects: one for the ILUβ-preconditioner and one for the BiCGStab iteration:
LexILUPrecond<double> lin_precond (SG_DISP_NO, -1);
CGIter<double> lin_solver (SG_DISP_RED, 1e-6, 100, 1e-15, &lin_precond);
Object lin_precond (the ILU-1-preconditioner) is referenced by lin_sover (the CG-iteration) and used from it implicitly. Analogously, the object lin_solver can be referenced in a non-linear solver (or applied by the user to solve sparse linear systems). The object lin_precond may be referenced in futher methods where a linear preconditioner is required. But the user may declare a further object of the class LexILUPrecond<double> (for example, with a different β).
The first argument in both the constructors in the example (and for many other numerical methods in SG 2) above sets a style of output. Numerical method may print some information like convergence rates that can be important for tuning and debugging. To manage the styles of the output, SG 2 has class SGDisplay which is a virtual base class for the numerical methods. There are several standard values for this style:
Some types of numerical methods (like the preconditioners) are not expected to print anything (but may do it with SG_DISP_DEBUG).
The main base classes for the numerical methods are:
 | (2.1) |
a linear solver finds u that satisfyes this equation for a given f, or an approximation to this u with a given precision. A typical example of a linear solver is the linear iteration that computes a sequence u(0), …, u(k), … of appxoximations of u using the recurrence
Here, W is a preconditioner matrix. Note that the computation of the product of f - Au(k-1) by W-1 is not a part of the solver. It is implemented in a separate class called preconditioner. For this, the solver may be used with any preconditioner.
Although these base classes are different in details, there are some common rules for all them. All these classes have virtual functions PreProcess and PostProcess and a function that implements the computation of the solution (like Run in SysPrecond<T>). The user should first call PreProcess. Additionally to some technical preparations (for ex. allocation of memory) this function may compute some auxiliary data. In the implementation of the ILUβ-preconditioner, PreProcess computes the ILU-decomposition. After PreProcess the user may call the function of the computation of the solution any numer of times. As soon as the user is done with the preconditioner, PostProcess should be called. PostProcess releases the allocated memory etc. Furthemore, all these classes have the virtual function PrintParam that prints current parameters of the object (like the specified precision).
Note that all the base classes for the numerical methods are templates depending on T. This T should denote an arithmetic type like double or float.
SG 2 is capable of reading settings stored in text format in initialization files (from now on called INI-files). The settings are stored in form of named parameters. The text in the INI-files consists of a series of assignements of values to these parameters. A typical statement in an INI-file is
<parameter name> = <value>;
for example,
N_op = 10;
Note that the semicolon is a part of the statement. Only constants are allowed in the right-hand side, no arithmetical expressions and no references to the other parameters. The value should correspond to the type of the parameter (integer, boolean etc.).
Parameters can be collected in directories. If parameters prec and N_op are in directory LinSolver then their values should be specified by the following sequence of the statements:
LinSolver {prec = 1e-12; N_op = 100;}
The directories may have subdirectories. Formally, all the parameters and directories are in the root directory. The structure of the parameters and the directories is similar to that of file systems. In this analogy, parameters are files.
Any text from character # up to the end of the current line is skipped (i.e. is treated as a comment). In the other text, the white space characters like spaces, tabs and the new-like charachters are skipped, too.
The structure of the parameters and directories is fixed by the application that uses the SG 2 library. In an INI-file, there must be no other parameters than those whose names are declared in the application. All the parameters should be in their directories. To parameters may have equal names only if they are in the different directories. A parameter may not be set two times. But the order of the assignments inside a directory is not fixed and may be arbitrary.
The SG 2 application may declare default values for the parameters. These parameters may be omitted in the INI-file. Then they are considered as assigned to the default value. All the other parameters should be present in the INI-file.
In the application, all the parameters are objects of classes derived from SGParam. The derived classes correspond to different types of the parameters. In particular, there are
The directories of parameters are parameters, too: They are objects of class SGParamDir derived from SGParam. All these objects should be created (for example, declared as static objects) in the applications before they are read from the INI-file. The name of the parameter is specifyed as the first argument of the constructors of the derived classes. The second argument is a pointer to the previously declared parameter in its directory. (For the first object of a parameter in some directory, this pointer should be NULL. The object of the directory itself should be created after all its member parameters.) By this pointers, the parameter objects of every directory, in particular of the root directory, are arranged in a single-linked list. The last parameter of every list (i.e. the head of this list) is then specifyed as the third argument of the constructor of its directory object. This results in a tree structure with the root in the last declared parameter object of the root directory.
For example, assume your application requires real parameter D for a diffusion coefficient and a directory LinSolver with three parameters: DispType for the type of the output of the linear solver, LinMaxIter for maximum number of iterations of the linear solver and LinReduction for the reduction factor of the norm of the residual. Then you declare the following objects:
SGRealParam D_param ("D", NULL);
SGDisplayParam lin_display ("DispType", NULL);
SGIntParam lin_max_iter ("LinMaxIter", &lin_display);
SGRealParam lin_reduction ("LinReduction", &lin_max_iter);
SGParamDir lin_solver_dir ("LinSolver", D_param, lin_reduction);
Here, objects lin_display, lin_max_iter and lin_reduction constitute the list of the subdirectory LinSolver (object lin_solver_dir). The list of the root directory consists of objects D_param and lin_solver_dir. The head of the root directory is lin_solver_dir. The specification of these parameters in the INI-file may look like
D = 0.1;
LinSolver {
LinMaxIter = 100; LinReduction = 1e-8; DispType = full;
}
To read the parameters from the INI-file, the application should call static function SGParam::Read. The first argument of this function is a C-string with the file name. The second is the pointer to the head of the root directory of the parameters. In the example above, we should call
SGParam::Read ("my_ini_file", &lin_solver_dir);
If the function fails to read the parameters (if no such file found or for illegal syntax, parameter type mismatch etc.), this function may print a message in the standard output and throw exceptions as described in Section 2.5. One may also print the read values of the parameters in the INI-file-format compatible form using function SGParam::Print.
Every class derived from SGParam provides its own means to access the value read from the INI-file. Cf. the Doxygen documentation for the details. Classes like SGIntParam and SGRealParam implement member function value that returns the value. Other classes have the second base class that corresponds to the type of the value. For example, SGIntArrayParam is derived also from MultiIndex so that its objects may be directly used where one needs a MultiIndex or an Array<int>.
There is also the second way to use the parameters. For many standard classes of the numerical procedures (cf. Section 2.2), the parameters of the method can be automatically initialized from the SG-parameters. This is done by the constructor of the class of the numerical procedures. This constructor has only one argument: the list of the parameters. (This list can correspond to the list of some subdirectory.) This parameters should be created as described above and have names predefined by this class. They should belong to one list. If some of them is not found by the class constructor in the given list, some default value is set. (Cf. the Doxygen documentation of the classes of the standard numerical procedures like the Newton method or the BiCGStab-method).
SG 2 has its own interface for the arithmetics of small dense matrices and small vectors. A small matrix is represented by an object of class small_matrix<T>. It can represent both small matrices and small vectors. For the small vectors, there is also class small_vector<T> derived from small_matrix<T>. Class small_matrix<T> supplies all arithmetic operations with matrices and vectors including computation of the Euclidean and the maximum norm, the Euclidean scalar product, the inverse matrix, as well as the solution of linear systems. There are also some further functions, for ex. the QR-method for the eigenvalue problems.
Small matrices and vectors arise when a discretization deals with vector-valued grid functions. These functions are then represented by objects MultiVector<T> (or MultiArray<T>) with an additional dimension and T being an arithmetical type. The index in the last dimension counts the components of the vector-valued degrees of freedom. For ex., if your problem is 2d and you need a grid function defined on a 20 × 30-grid with values in R2 then you should declare a three-dimensional MultiVector<double>. Its first two sizes are 20 and 30 whereas the last one is 2. If gf is such an object then for a two-dimensional MultiIndex I, expression gf (I) returns a pointer to the R2-vector at the grid point referenced by I. This pointer can be then used to initialized an object of small_vector<double> so that this vector can be involved in the arithmetic operations. (Note that it would be wrong to declare MultiVector< small_vector<double> > to describe such a grid function).
For such problems, the stiffness and mass matrices consists of entries that are small matrices. The entries of these small matrices are stored by the objects of SparseMatrixLine<T> as raw arrays. Initializing objects of small_matrix<T> with pointers to these arrays, you may involve these small matrices in the matrix-vector arithmetics.
Objects of class small_matrix<T> do not store the entries of the matrices. Instead, they store the sizes and a pointer to an array of the entries. The latter pointer is passed to the class constructor at the initialization of the object. Class SmallMatrix<T> derived from small_matrix<T> allocates the array itself. This class can be used whenever a temporary small matrix or a temporary small vector is required. Nevertheless note that the allocation and the disposition of the arrays of the entries is a time-consuming operation, and this class should not be used frequently.
The template parameter T of the template small_matrix and the derived template classes should denote an arithmetical type like double, float, a type for complex numbers etc.
On unexpected situations like misspecified parameters or lack of memory, SG 2 throws C++-exceptions. All these exceptions are objects of classes derived from sg_fatal_error. This base class contains at least a pointer to a static string with an error message. But if the library is configured with the option --enable-fullfaterr (cf. Section 1.1), it contains also a static string with the file name of an SG-module and the line number where the exception was thrown. This option is also activated by --enable-debug.
Recognition of the unexpected situations requires time and slows down the computations. A particular example of this is checking the range of indices when accessing elements of an Array<T> or a MultiArray<T>. To prevent the loss of efficiency, SG 2 allow to switch the verification off. The verification is on only if macro CHECK_ARRAY_RANGE is defined. This can be done by the configure options --enable-range or --enable-debug. Note that if this macro is off, an out-of-range index induces other unexpected situations in which other exceptions may be thrown. For this, we recommend to keep the range checking on at the debugging phase.
In the SG 2 modules, the exceptions are usually thrown using macro SG_FATAL_ERROR. This is done to specify the additional arguments for the constructor of sg_fatal_error depending on the specified configuration options.
When an SG 2 exception is cought, the user may call function process for its object. This function prints the error message and halts the program. For a parallel program, this function properly halts the parallel computation. But the user may process this exception in any other way.
