Alfonsus Mario
114081187
Abstract:
This paper exploits object-oriented implementation techniques to facilitate the development of computer codes for solving systems of coupled partial differential equations. We show how to build a simulator for equation systems by merging independent solvers for each equation that enters the system. The main goal is to obtain a rapid, robust, and reliable software development process with extensive reuse of implemented code. Coupled heat and fluid flow in pipes is used as example for illustrating the implementation techniques. We also present some results for the particular case of temperature-dependent generalized Newtonian fluid flow between two nonconcentric cylinders. The general applicability of the approach is discussed.Categories and Subject Descriptors: D.1.5 [Programming Techniques]: Object-Oriented Programming; G.1.8 [Numerical Analysis]: Partial Differential Equations--Finite-element methods
INTRODUCTION
Programming and debugging simulation codes that involve numerical solution of partial differential equations (PDEs) often take an undesirable large amount of time. To reduce the efforts spent on software issues, one can advantageously apply modern software development techniques to make the implementational aspects cleaner, simpler, and more effective. For example, object-oriented programming (OOP) offers features that may significantly improve the design, implementation, and maintenance of large simulation codes. This has been recognized for a long time in the computer science community, but the severe efficiency requirements of numerical codes have prevented widespread application of OOP among computational scientists until recent years. The C++ language, better optimization modules in compilers, and knowledge about inefficient C++ constructs have contributed to increased interest in exploring C++ and OOP techniques for scie! ntific computing. This is reflected in the many conferences an! d litera ture contributions on the subject [Arge et al. 1997a; 1997b; Norton 1996]. Most of the contributions deal with OOP techniques applied to the various steps in certain numerical methods for partial differential equations (PDEs), with the main emphasis on finite-element methods. The PDE solver is then a short program written at a high abstraction level.
The present paper also treats OOP techniques for PDEs, but in contrast to most literature contributions, we focus on a new technique where a solver for a PDE is an independent standalone object that can be combined with other independent solver objects in order to easily assemble software for systems of PDEs in a flexible manner. The use of OOP techniques, as presented in this paper, provides a foundation for building advanced engineering and scientific applications in a modular approach, with extensive reuse of existing code. The first steps in this direction were taken by Bruaset et al. [1997], but in addition to de! velop their ideas further, we also suggest how to handle various types of constitutive relations. The exposition in the present paper is mainly restricted to an operator-splitting technique for solving system of PDEs, where each equation is solved in sequence. The basic ideas of the paper can be extended to handle fully implicit methods as well [Odegard et al. 2000].
The physical problem of consideration is pipe flow problems for generalized Newtonian fluids exhibiting thermal effects. The properties of these fluids are often sensitive to the temperature, and the effect of viscous dissipation may be an important heat generation mechanism. Many industrial processes, such as different molding techniques, fuel technology, lubrication applications, polymer melt flows, etc., involve flow of viscous non-Newtonian fluids through pipes having various cross-section shapes. For example, mathematical models for these flow phenomena can be valuable for engineers in predictin! g temperature distribution, regions of low shear rate and/or p! oor flow -through, or to determine the relation between the pressure gradient and the volume flux, and how these quantities depend on some key characteristics of the fluid's constitutive properties. The numerous applications are reflected in the many studies on the topic reported in the literature [Ahmadi 1985; Gupta and Massoudi 1993; Heyes and March 1994; Johnson et al. 1991a; 1991b; Massoudi and Christie 1993; 1995; Papachristodoulou and Trass 1987; Ramkissoon and Easwaran 1993; Skelland 1969; Szeri and Rajagopal 1985; Winther 1977].
The particular problem is chosen so as to obtain a mathematical model that is complicated enough to demonstrate several key aspects of our programming technique in general, yet simple enough to keep the attention to the basic implementation issues. Although the technique is beneficial even in the present problem, the impact of the suggested approach will be more substantial in PDE models with greater complexity [Samuelsson 1996], as we outl! ine in the final part of the paper.
The need for a flexible and extensible design of computer codes is evident in many ways. For example, considering the present physical problem, a generalized Newtonian fluid model may adequately fit the shear stress and shear rate measurements for many non-Newtonian fluids, but may be inappropriate to predict other flow phenomena [Massoudi and Christie 1995; Schowalter 1978]. To obtain useful results it may therefore be necessary to modify the mathematical model, for example to higher-order [Massoudi and Christie 1993; 1995; Ramkissoon and Easwaran 1993] non-Newtonian fluids, resulting in the need for either to extend or to replace parts of a simulator for the overall problem of interest. Moreover, a mathematical model may also include additional PDEs representing chemical reactions, turbulence, multiphase flow, etc. With an appropriate software design for solving systems of PDEs, one can plug and play with different modeling a! pproaches in a simple manner, and hence the time spent on impl! ementati onal work can be minimized. Furthermore, the numerical solution of the coupled nonlinear partial differential equations that arises may also easily face convergence problems, depending strongly on some problem-dependent parameters. To obtain reliable conclusions regarding the simulation results, it is therefore important to have access to several numerical formulations and solution strategies. This calls for flexible software with a wide range of methods available, and where new solution approaches and constitutive relations can easily be incorporated in a simple manner. Hence, flexible software is a key issue to make progress in the present and various other fields. We shall show that OOP is a convenient tool for achieving the necessary flexibility.
The pipeflow solver has been implemented using utilities available from the Diffpack software library.(1) Some simplified examples from this implementation appear in the text, and basic knowledge of OOP and the C++ la! nguage will be beneficial for the reader.
The paper is organized as follows. In the next section, we present the mathematical problem and numerical formulation. Thereafter, in Section 3, we give a detailed description of our implementational ideas, and in Section 4 we use our simulator to obtain some results for a particular case. In Section 5 we conclude with a discussion of how the present implementation ideas can be easily extended to address more complex problems.
(2.) PHYSICAL AND MATHEMATICAL PROBLEM
2.1 Partial Differential Equations
We consider stationary single-phase flow of an incompressible generalized Newtonian fluid in a straight pipe with cross-section [Omega]. The basic equations, involving mass, energy, and momentum balance, coupled through constitutive laws, read the following:
(1) [v.sub.k, k] = 0,
(2) ?? ([v.sub.i, t] + [v.sub.k] [v.sub.i, k]) = [-p.sub.,i] + [[Tau].sub.ij, j],
(3) ?? ! [C.sub.p]([T.sub.,t] + [v.sub.j][T.sub.j]) = [Kappa][T.sub.,jj! ] + [Mu] [[Gamma].sup.2],
(4) [[Tau].sub.ij] = [Mu]([v.sub.i,j] + [v.sub.j,i]),
(5) [Mu] = [[Mu].sub.T](T)[[Mu].sub.w]([Gamma]),
(6) [Gamma] = [square root of 1/2 ([v.sub.i,j] + [v.sub.j,i])([v.sub.i,j] + [v.sub.j, i])]
Here [v.sub.i] is the velocity field; T is the temperature; p is the pressure; ?? is the fluid density; [C.sub.p] is the heat capacity at constant pressure; [Kappa] is the (constant) heat conduction coefficient; [[Tau].sub.ij] is the deviatoric stress tensor; and [Gamma] is the shear rate measure. In these expressions, comma denotes derivation, and sum over repeated indices is implied. The function [[Mu].sub.w] ([Gamma]) depends on the constitutive properties of the fluid, and the function [[Mu].sub.T](T) represents the temperature dependency of the fluid viscosity and couples the momentum and energy equations.
For flow in straight pipes it is natural to restrict the velocity to be unidirectional, i.e., {[v.sub.i]! } = (0, 0, w). Together with the assumption of stationary flow, this leads to dramatic simplifications of the momentum and energy equations. Appropriate assumptions on the temperature boundary conditions lead to the associated assumption T = T(x, y). We then achieve the simplified system
(7) [Delta] [multiplied by] [[Mu] [Delta] w] = -[Beta],
(8) [[Delta].sup.2] T = - [[Kappa].sup.-1] [Mu] [[Gamma].sup.2],
(9) [Mu] = [[Mu].sub.T](T) [[Mu].sub.w]([Gamma])
(10) [Gamma] = [square root of [([w.sub.,x]).sup.2] + [([w.sub.,x]).sup.2].
Here, [Delta] = ([differential]/[differential]x, [differential]/[differential]y), and the parameter [Beta] = -[p.sub.,z] is the constant pressure drop along the pipe.
There are numerous generalized Newtonian fluid models in the literature [Schowalter 1978; Skelland 1969]. A quite general and widespread model for [[Mu].sub.w] is the Sisko model:
(11) [Mu] - [[Mu].sub.[infini! ty]] + [[Mu].sub.0] [[Gamma].sup.n-1],
where [[Mu].! sub.[inf inity]] is the viscosity at extreme shear rates; [[Mu].sub.0] is a reference viscosity; and n is the "power-law exponent." When [[Mu].sub.[infinity]] = 0, this model reduces to the standard power-law model, and the choice [[Mu].sub.0] = 0 leads to Newtonian flow with viscosity [[Mu].sub.[infinity]]. Common values of n are in the interval (0.15, 0.6), but there also exist materials for which n is greater than unity (dilational fluids). The Herschel-Bulkley model is more general:
(12) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]
where [Tau] = 2 [Mu] [Gamma] effective stress for plastic flow (similar to von Mises stress in solid mechanics), and [[Tau].sub.0] = is a critical value of [Tau] for the transition between a rigid-body movement of the fluid ([[Mu].sub.[infinity]] [right arrow] [infinity]) and a modified power-law behavior. Notice that [[Tau].sub.0] = 0 recovers the standard power-law model, whereas n = 1 corresponds to a Bingham model. The ! Cross model takes the form
(13) [Mu] = 2[[Mu].sub.0]/1 + [([[Mu].sub.0][Gamma]/[[Tau].sub.0]).sup.1-n],
where [[Tau].sub.0] now is the shear stress level at which the flow undergoes a transition from the Newtonian nature to the power-law region.
The function [[Mu].sub.T] can also be chosen in various ways. The Reynolds model for [[Mu].sub.T] reads
(14) [[Mu].sub.T](T) = exp(-[Alpha](T - [T.sub.0])),
where [[Tau].sub.0] is a reference temperature. An other alternative is the Vogel model
(15) [[Mu].sub.T](T) = exp([Alpha]/[Beta] + (T - [T.sub.0]).
In the case where [[Mu].sub.T] is independent of T, the temperature does not enter the momentum equation. This allows us to first solve for w (x, y) from the momentum equation and thereafter solve a Poisson equation, with known right-hand side, for T. Under certain circumstances, one can then find analytical expressions for w and T. Such solutions are funda! mental for the verification of computer implementations.
!
2. 2 Boundary Conditions
The physical problem was taken to be pipe flow in the above derivation, but related flow cases lead to the same PDEs, though with some different boundary conditions. We, therefore, introduce some different boundaries to increase the flexibility of the problem formulation. Let [differential][[Omega].sub.1] be the part of the boundary corresponding to a fixed wall, where w = 0; let [differential][[Omega].sub.2] be a wall moving with velocity [W.sub.1] in z-direction; let [differential][[Omega].sub.3] be a plane, where [differential]w/[differential]n = 0 or [differential]T/[differential]n = 0; let [differential][[Omega].sub.4] be a wall with fixed temperature T = 0; let [differential][[Omega].sub.5] be another wall with fixed temperature T = [T.sub.1]; and finally let [differential][[Omega].sub.6] be a wall where a cooling condition [differential]T/[differential]n = - [h.sub.T](T - [T.sub.s]) applies. Here, [h.sub.T] is a coefficient that reflec! ts the heat transfer through the channel wall to the pipe's surroundings, which is assumed to have a constant temperature [T.sub.s]. Some of the boundary parts [differential][[Omega].sub.i] may be overlapping, and some parts may be empty. The basic requirement is that we need exactly one condition on w and one condition on T at every point on the boundary.
2.3 Numerical Solution Method
We solve the Eqs. (7)-(8) by a standard Galerkin finite-element method. Let
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII],
where [N.sub.j](x, y) are finite-element basis functions. The discrete equations then take the form
(16) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII],
(17) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII],
for i = 1, ... , m. This is a set of 2m coupled nonlinear algebraic equations. The expressions for [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] and [MATHEMATICAL EXPRESSI! ON NOT REPRODUCIBLE IN ASCII] are given below.
(19)! [MATHEM ATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII],
(20) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].
We will consider two basic solution strategies: (1) a Gauss-Seidel or Jacobi iteration scheme between the two equations, and (2) a fully implicit method. The latter strategy consists in forming a system of 2m nonlinear algebraic equations directly from the [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] and [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] equations and then solving this system by a standard method like Newton-Raphson or Picard iterations.
In the Gauss-Seidel approach we first solve the momentum equations [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], with respect to [w.sub.j], using the [T.sub.j] values from the previous iteration. Thereafter we solve the energy equations [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], with respect to [T.sub.j], using the recently computed [w.sub.j] values. If we introduce! a superscript (q) for the approximation of a quantity in the qth iteration, the Gauss-Seidel algorithm can be written more precisely as
for q = 1, 2, ... until convergence
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]
Each problem in the Gauss-Seidel approach is nonlinear and can be solved by basically two strategies; either a fully nonlinear solution process based on Newton-Raphson or Picard iteration, or just one iteration of a Picard scheme. The latter strategy will most probably lead to more iterations in the outer Gauss-Seidel (q) loop, but can in some cases be essential for obtaining convergence. The Gauss-Seidel approach with only one Picard iteration pass in each q-loop step is the simplest method to implement, but the fundamental question is whether the convergence rates are satisfactory enough. The Jacobi method is similar to the Gauss-Seidel scheme except that we use the "old" values [w.sup.(q-1)] when solving the [MATHEMAT! ICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] equation in iterati! on q.
2.2.1 Initiating the Iterative Process. To obtain a rapid convergence, or in many cases convergence at all, one has to initialize the iteration scheme with an initial guess for the solution that in some sense is close to the real solution. A good initial guess may be difficult to come up with in some cases, and therefore we also explore some other strategies for establishing convergence in difficult cases.
Newton iterations generally give the best convergence performance if the initial guess is good, but the Picard method is often more robust, i.e., less dependent on the initial guess. In order to obtain the best from both these strategies, we allow for iterating with the Picard method until a given criteria is fulfilled, and then switch over to the Newton iterations using the Picard solution as the initial guess. It is also allowed to go back to the Picard scheme again if that should be desirable for some reason.
Another alternative is to try ! the continuation method. In this method, one or more parameters in the problem are first chosen such that the problem becomes easy to solve. The values of these parameters are then gradually changed toward the desired values. In the present case, the obvious candidates for such a strategy are the parameters [Alpha] and n. With n = 1 and [Alpha] = 0, we know that the problem is easy to solve, whereas convergence problems are expected as n and [Alpha] move away from these values. A possible formulation of a continuation method in terms of n in our problem will then be using [Lambda] = (1 - n)/(1 - [n.sub.0]) as continuation parameter, with [n.sub.0] being the target value of n for the computations. Now [Lambda] = 0 is easily solved, and [Lambda] = 1 is the target problem. By defining a set of proper values [[Lambda].sub.0] = 0 [is less than] [[Lambda].sub.1] [is less than] ... [is less than] [[Lambda].sub.p] = 1 of [Lambda], and using the solution obtained with [[Lambda].sub.! i-1] as initial guess for the nonlinear solvers when solving t! he probl em corresponding to [[Lambda].sub.i], we might hope to establish convergence for extreme n values. With the same definitions as above, the corresponding formulation for a will then be [Lambda] = [Alpha]/[[Alpha].sub.0].
2.2.2 Summary of Solution Strategies. Let FI denote the fully implicit methods where we solve 2m equations simultaneously, either by Newton-Raphson's method (FI-NR) or by Picard iteration (FI-P). The number of iterations to reach a stopping criterion is denoted by Q. The family of Gauss-Seidel methods consists of one-pass Picard iterations for each equation (P1), and full Newton-Raphson or Picard iteration for each equation (NR or P). In the latter case we can choose Newton-Raphson or Picard iteration freely for each equation. Hence, we introduce the notation NRm-Pe to represent the outer Gauss-Seidel iteration with a full inner Newton-Raphson method for the momentum equation (lowercase m) and a full inner Picard iteration scheme for the energy equ! ation (lowercase e). The prefix J denotes the Jacobi-type sequential solution method. The meaning of notations like J-Pm-Pe should then be clear. Switching between the Picard and Newton schemes is denoted by an extra hyphen. For example, (Pm-NRe)-NR denotes Gauss-Seidel outer iterations starting with the Newton scheme in the energy equation and the Picard scheme in the momentum equation and then finish with Newton iterations in both equations after some switch criterion is fulfilled.
3. IMPLEMENTATION
The goal of this section is to provide a set of safe and efficient software development steps toward a simulator that is capable of solving the present flow problem using all the numerical methods listed in the previous section. Choosing the specific solution strategy and the type of constitutive law at run-time is an important requirement of the desired flexibility.
3.1 PDE Solvers in Diffpack
Diffpack is a programming environment f! or scientific computing, with special emphasis on numerical so! lution o f PDEs. The package makes use of object-oriented design and programming techniques and is implemented in C++. In the Diffpack libraries one can find useful building blocks for PDE solvers, realized as C++ classes. Arrays, linear systems and solvers, preconditioners, finite-element grids and associated scalar and vector fields, interfaces to visualization and grid generation, finite-element assembly algorithms, etc. are examples on class modules that are used when building PDE solvers. The PDE solver itself is also a C++ class.
Figure 1 shows the main class relations for a typical finite-element PDE solver in Diffpack. Among the internal variables (boxes being pointed to by dashed arrows) in this class we find a pointer to a finite-element grid (GridFE). The GridFE class holds information about the nodal coordinates, the element topology, the element types, as well as boundary nodes. A separate class FieldFE represents a scalar finite-element field, and offers vari! ous functionality for interpolating the field at a given spatial point. The field class has typically a pointer to a finite-element grid object (GridFE) and an array of nodal values. We mention that our term "pointer" actually refers to a kind of smart C++ pointer with built-in reference counting and garbage collection. Another class, DegFreeFE, keeps track of the mapping between local (elementwise) and global degrees of freedom in a linear system. Representation and solution of linear systems are taken care of by the LinEqAdm class. With this class, at run time the user can set a basic solver, preconditioner, and matrix storage scheme appropriate for the size and type of problem at hand. Our PDE solver class also has a NonLinEqSolver component for solving systems of nonlinear equations. This nonlinear solver requires the PDE solver to be derived from a base class NonLinEqSolverUDC, which we will explain later. Another base class of PDESolver is class FEM, containing tools ! and interfaces for finite-element methods.
[ILLUSTR! ATION OM ITTED]
Figure 2 lists the major components of a typical definition of the PDE solver class in the case where we solve a Poisson equation [Delta] [multiplied by] [Lambda][Delta]u = 0 with some Neumann and Dirichlet conditions, and [Lambda] is a prescribed coefficient. In addition to the pointers described above and some internal variables, this class has some functions, referred to as methods, for carrying out the steps in the solution process. The integrands method evaluates the integrands in the weak formulation of the boundary-value problem. The lambda method evaluates the [Lambda] function at a spatial (integration) point. The define method defines a set of menu items for all input data to the simulator. The scan method loads the answers from the menu and initializes the internal data structures. The solve method calls a general finite-element assembly routine in FEM, which jumps back to the solver class in the integrands routine, and thereafter the linear syst! em solver in LinEqAdm to compute the field u. A method fillEssBC sets the essential (Dirichlet) boundary conditions. Natural boundary conditions appearing in the weak formulation are integrated in a integrands4side method which is similar to the integrands method, except that it is only evaluated at the selected boundaries of the elements. For nonstationary problems, an additional method solveProblem is introduced for running the time loop.
Fig. 2. Sketch of typical PDE solvers: (a) linear Poisson equation
[Delta] [multiplied by] [Lambda][Delta]u = 0, (b) nonlinear Poisson
equation [Delta] [multiplied by] ([Lambda](u)[Delta]u) = 0. The
asterisk indicates a pointer (but in our implementation this is
actually a smart pointer with reference counting).
(a)
Solver1
variables:
GridFE* grid
FieldFE* u
LinEqAdm* lineq
methods:
integrands
lambda
define
s! can
solve
fillEssBC
(b)
! Solver2< br />
variables:
GridFE* grid
FieldFE* u
LinEqAdm* lineq
NonLinEqSolver_prm pm
NonLinEqSolver* nlsolver
methods:
integrands
lambda
define
scan
solve
fillEssBC
makeAndSolveLinearSystem
When dealing with nonlinear PDEs, we need additional internal variables pointing to a nonlinear solver class (NonLinEqSolver) and its parameters (NonLinEqSolver_prm) (Figure 2(b)). Most nonlinear solution algorithms perform some kind of iteration, where a linear problem is formulated and solved in each pass of the iteration loop. Class NonLinEqSolver is designed such that the formulation and solution of the linear problem is left to the PDE solver. In Diffpack this is accomplished by viewing all PDE solvers, or any other class needing to solve systems of nonlinear algebraic equations, as special instances of NonLinEqSolverUDC (NonLinEqSolver's User Defined Code). The NonLinEqSolve! rUDC class defines an interface for all nonlinear problems. The most important function in this interface is makeAndSolveLinearSystem, which requires the user's class, here PDESolver, to formulate and solve the linear subsystem in a nonlinear iteration. The NonLinEqSolver class has a pointer to NonLinEqSolverUDC object (which in our example here happens to be PDESolver). The nonlinear solution algorithm can just ask the NonLinEqSolverUDC pointer to call makeAndSolveLinearSystem to solve the linear subsystem of the current iteration, without knowing anything about the user's class where this function is defined. In other words, the nonlinear solution algorithm can be coupled to any class defining a nonlinear problem as long as the class implements a function with the name makeAndSolveLinearSystem and declares NonLinEqSolverUDC as one of its base classes. This is the heart of object-oriented programming and is a convenient tool for separating general library functionality fro! m application code. Without object-orientation, one must rely ! on sendi ng a function pointer for makeAndSolveLinearSystem to the nonlinear solution algorithm. The disadvantage with such an approach is that it becomes tricky to transfer application-dependent data structures needed in makeAndSolveLinearSystem.
Let us look at the NonLinEqSolver class to see another example on object-oriented programming. The NonLinEqSolver is actually not a specific nonlinear solver, it is just a definition of a common set of functions and data needed in nonlinear solvers. Specific solution algorithms are implemented as subclasses of NonLinEqSolver. The use of OOP techniques makes it possible to work exclusively with the general NonLinEqSolver interface in our PDE solver without knowing explicitly the specific nonlinear solution algorithm that is currently in use. For example, NewtonRaphson is a subclass of NonLinEqSolver and implements a Newton-Raphson method. The subclass inherits data and functions from its base class; however, one function solve is ! "empty" in NonLinEqSolver, and this function must be reimplemented in its appropriate form in the NewtonRaphson class. In other words, to implement a new nonlinear solver, we just derive a subclass and put the solution algorithm inside the solve function. Common operations such as convergence criteria and management of data structures are left to the base class. Another subclass Picard implements simple Picard iteration (also called successive substitution). To make the program capable to distinguish between the nonlinear solvers at run time, the solve function must be virtual.
For the compiler, instances of the above subclasses are also considered as NonLinEqSolver objects. This means that NewtonRaphson and Picard objects can be viewed as NonLinEqSolver objects from our simulator code. The advantage is that we hide details concerning different nonlinear solution algorithms, and we can write our PDESolver code in terms of some common abstract nonlinear solver. Si! nce the solve function is virtual, C++ will translate a statem! ent nlso lver->solve into a call to the specific solve function implemented in the specific nonlinear solution algorithm chosen at run time, say Picard. That is, writing nlsolver->solve makes C++ call the Picard's class solve function. In other words, C++ keeps track of some details at run time and lets the programmer write code at a higher abstraction level.
The use of subclasses and virtual functions is crucial when realizing our suggested software framework. Building the code in this way makes, in the authors' opinion, the implementation cleaner and more flexible. Diffpack uses the same principles for handling linear solvers (class LinEqAdm), matrix formats, finite elements, and so on. Programming in terms of abstract interfaces to numerical methods and data structures hides implementational details and simplifies the software development and maintenance process significantly.
We should mention that while the C++ language provides the concepts of subclasses a! nd virtual functions, these concepts are respectively known in object-oriented terminology as inheritance and dynamic binding. The realization of these concepts varies between the different object-oriented programming languages. For example, our proposed design is easily implemented in Java, but not in Fortran 90.
Some Remarks on the Computational Efficiency. In terms of computational efficiency, a sensitive part of the finite-element solution process is to compute and assemble the linear equation system. This process easily becomes more inefficient because of the suggested flexibility. This is mainly due to computations and use of variables that are not always necessary or could be strongly optimized for a specific application and a particular element type (in the general case we are forced to work with an unspecified number of space dimensions, arbitrary elements, etc.).
However, steps for speeding up the assembly process can be taken. Experience wit! h optimization of this type of C++ code shows that performance! improve ments can be obtained by tailoring the finite-element assembly process at the particular PDE and element in question. This reduces the flexibility and specializes the code, but can easily be accomplished by introducing small subclasses of each solver, where the virtual assembly function is hardcoded for the given PDE and element choice. See Langtangen [1999] for detailed descriptions.
Another computationally intensive task is the process of solving linear systems of equations. Since this task mainly consists of matrix-vector arithmetic that can be carried out using optimized libraries, the efficiency will not be influenced by our use of object-oriented implementation techniques. In fact, the use of such techniques allows for a flexible and simple use of strategies for improving the efficiency of the solution process, such as parallel processing, adaptivity, multigrid, and domain decomposition methods [Langtangen 1999].
3.2 Assembling Independent Solver ! Components to Solve a PDE System
In principle, we can use the general set-up of a nonlinear PDE solver for most stationary PDE problems. For example, when solving the coupled w-T problem from Section 2, we can let the solver class have two fields, w and T, and in integrands and other functions we can have an indicator telling whether to deal with the momentum equation or the energy equation (or a simultaneous simulation of both). However, our experience shows that the software development process can be made faster and more reliable by first creating a solver for the momentum equation, like the one sketched above, and debug this solver for a simple choice of the variable coefficient (viscosity [Mu]). Similarly, we develop a solver class for the energy equation, where a known function is used for the right-hand side (this energy equation solver class will then typically have a method rhs that replaces the variable coefficient method lambda).
Ideally, we! should then be able to stick the well-tested, completely unco! upled mo mentum and energy classes together, without touching the source code, and thereby obtain a solver for the compound system. The reason why we want to keep the original code untouched is three-fold:
(1) Any modification increases the danger of introducing errors.
(2) When testing and debugging the solver for the w-T system, we want the possibility to easily pull the software apart and retest the momentum and energy equation solvers independently.
(3) When more sophisticated, highly-efficient, adaptive or parallel Poisson solvers become available as Diffpack classes, we would like to incorporate these into our w-T system without any need to touch the inner delicate numerics in these solvers (e.g., integrands functions).
Assume that we have a momentum equation solver available for [Delta] [multiplied by] [Lambda][Delta]w, as sketched above, but only tested with a simple choice of [Lambda], for which we can find analytical solutions. Le! t Momentum1 be the name of this solver class. Similarly, we have developed and tested a solver class Energy1 for the energy equation [[Delta].sup.2]T = f(T). In this class we have made a simple choice of the right-hand-side function f(T) for testing purposes only.
The development of a solver for the complete equation system then follows three fundamental design principles.
(1) The variable coefficients, including the right-hand sides, in the PDEs (e.g., the lambda method) are represented by virtual functions. Subclasses of a PDE solver override these functions and implement the physically relevant versions, when the coefficients are coupled to other unknown fields in the PDE system.
(2) The equations become coupled through variable coefficients. These coefficients are often built of a common set of relations (constitutive laws). The common relations are collected in class hierarchies and accessed from the variable coefficients through a base! class interface (pointer). Hence, constitutive laws, etc. can! be easi ly changed without affecting the code in the PDE solvers.
(3) A manager class, called Manager, acts as the solver class for the whole PDE system. The class contains two-way pointers to base classes for solving the momenta and energy equations, shown in Figure 3 as FlowSolver and EnergySolver, respectively, as well as a pointer to the common relation hierarchy. The manager is also responsible for creating a grid, and perhaps allocate a common linear system and solver object, and distributing these to the solver classes for the various PDEs in the system. The code in the manager class is usually very small, as the goal is to simply glue the various component solvers and the common relations.
[ILLUSTRATION OMITTED]
The third point requires that the initializing (scan) methods in Momentum1 and Energy1 can either make the grid and allocate the linear/nonlinear system and solvers themselves, or simply set their pointers to point to structures made ! outside the classes (i.e., Manager in this case). These pointers are provided as parameters to the PDE solver's scan method, with a null pointer indicating local allocation of a data structure. Inside a PDE solver we cannot see whether we are working on external or local data structures; in a sense, all data structures appear to be local.
Figure 3 shows a sketch of the relation between the classes. The subclasses Momentum2 and Energy2, which enables the coupling of the standalone solvers Momentum1 and Energy1 to the system, are very short; their contents are limited to overrided versions of the virtual functions in the solver classes that evaluate the variable coefficients in the PDE at a point. Here, Momentum2 will have a lambda method that corresponds to the real one in the w-T system, i.e., [Lambda] = [[Mu].sub.T](T) [[Mu].sub.w]([Gamma]). The information about the constitutive relations [[Mu].sub.T] and [[Mu].sub.w]([Gamma]) are managed by the CommonRel class! ; see Figure 4(b). The lambda method therefore needs to access! CommonR el to compute this variable coefficient. Recall that the Manager class has a pointer to the CommonRel class (see also Figure 4(a)), such that an implementation of a function lambda (x) in Momentum2, where x is a spatial point, may in C++ read return manager->re1->mu (x). Looking at class CommonRel in Figure 4(b), we see that one of the variables is indicated as tables.... This is meant to represent a table of, for example, T, w, [Delta]w, [Gamma], [[Mu].sub.w]([Gamma]), [[Mu].sub.T], etc. at a spatial point such that the functions mu_w and mu_T can be called repeatedly without the need to reinterpolate the T and w fields. At each (integration) point where one needs to evaluate the constitutive relations, one first calls the tabulate method to make internal evaluations in CommonRel, and thereafter the various relations mu_w and mu_T can be called very efficiently. In the case where various constitutive relations have common subexpressions, the tabulate method is even more imp! ortant for the computational efficiency than it is in the present application.
Fig. 4. (a): Sketch of the Manager class for administering the
solution of the w-T system of PDEs, using standalone momentum and
energy equation solvers. (b) Sketch of the base class CommonRel, with
a pool of common relations for the system of PDEs, including the
constitutive laws for [[Mu].sub.w] and [[Mu].sub.T].
(a)
Manager
variables:
CommonRel* rel
Momentum2* weq
Energy2* Teq
GridFE* grid
methods:
define
scan
solve
{ while not converged
weq>solve
Teq> solve
}
(b)
CommonRel
variables:
FieldFE* w
FieldFE* T
tables ...
methods:
muw
muT
dmuw_dwj
dmuT_dTj
srate
dsrate_dwj
tabulate
Different! constitutive laws are implemented as subclasses of CommonRel,! e.g., a power-law model can be represented by a class PowerLaw as sketched in Figure 3. Observe that the handling of constitutive relations in terms of a CommonRel hierarchy is a general idea that is applicable to much more complicated systems of PDEs. An additional advantage is that it also separates the work of, for example, mechanical engineers and chemists. The latter can develop the CommonRel hierarchies without any knowledge of the PDE solver or the numerical solution process, and numerical experts can work on, for example, the Energy1 class and fast adaptive multigrid methods without addressing all the details of constitutive relations in the dissipation term.
Newton-like methods require differentiation of the quantities in the constitutive relations with respect to the nodal unknowns [w.sub.j] and [T.sub.j]. For example, in the present case we need
[differential][[Mu].sub.w]/[differential][w.sub.j], [differential][[Mu].sub.T]/[differential][T.sub.j], [! differential][[Mu].sub.w]/[differential][T.sub.j], and [differential][[Mu].sub.T]/[differential][w.sub.j].
In a general case with a nonlinear function f(u, [u.sub.,1], [u.sub.,2], [u.sub.,3]), where
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII],
one needs
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII],
which is conveniently implemented in a function calculating the vector [differential]f/[differential][u.sub.j] at an evaluation point inside the element. Here, j = 1, ..., [n.sub.e], where [n.sub.e] is the number of degrees of freedom of u in an element. Class CommonRel must hence offer the functions dmuw_dwj for [differential][[Mu].sub.w]/[differential][w.sub.j], dmuT_dTj for [differential][[Mu].sub.T]/[differential][T.sub.j] and so on for analytical computation of the Jacobian in Newton-like methods.
The Manager class in Figure 4(a) is a simple class with pointers to the common relations and the so! lver for each PDE in the system. Its define and scan functions! make ca lls to the define and scan functions in each solver class, in addition to handling input parameters for the overall numerical solution strategy (Gauss-Seidel vs. Jacobi iteration) and the grid. We have made a brief sketch of the solve function, which is simply a loop over calls to the solve function in each solver class.
3.2.1 Making the Solver More Flexible. There are obvious arguments against our suggested approach. One such argument is that class Manager is a specialized class that cannot be reused, i.e., there is no general interface between class Manager and the momentum and energy solvers. Since the manager class is quite short, a tailored glue between the PDE solvers is not necessarily a disadvantage. However, in other applications the manager class might be comprehensive, and one may want to choose between different versions of, for example, momentum equation solvers. In such cases, one can introduce a common interface FlowSolver to momentum equation solve! rs and a similar interface EnergySolver to energy equation solvers. Our special PDE solvers, Momentum2 and Energy2, must then also inherit the interfaces from FlowSolver and EnergySolver, respectively. The resulting class diagram is shown in Figure 5. A typical interface in FlowSolver and EnergySolver would be a scan function for initialization and a solve function for solving the PDEs. In an alternative design the Momentum1 solver could be derived from FlowSolver, thus avoiding the multiple inheritance in Figure. A disadvantage of this approach is the implied coupling of a standalone Poisson-equation solver Momentum1 and our framework for solving systems of PDEs. For example, with our suggested approach, one could easily replace Momentum1 and Energy1 by some class Poisson that solves a general Poisson equation, and derive Momentum2 and Energy2 from Poisson. Anyway, the discussion in this paragraph should reveal that it is quite easy to alter the design in various ways with! no or minor impact on the well-tested numerical parts of the ! code.
[ILLUSTRATION OMITTED]
3.2.2 Fully Implicit Solution of the Equation System. Another argument against our approach is that the design is tailored to sequential solution methods and does not allow fully implicit methods. The fully implicit approach has a linear equation system in each nonlinear iteration with size 2m x 2m, and the elemental matrix and vector now involve full coupling between the two equations. That means that the element matrix of the momentum equation solver contains terms directly related to the energy equation and vice versa. From a pedagogical point of view, we have found it cleaner to develop a fully implicit solver separately from the momentum and energy PDE solver classes, but reusing the CommonRel class. In this case the manager class also has a pointer to a FullyImplicit object, and a switch set at run time determines the solution strategy.
We have extended the suggested design for sequential solution methods to allow f! or fully implicit schemes as well. The class design is pretty much the same as for a sequential solution method (cf. Figures 3 or 5); however, the manager class must be responsible for the linear system assembly and the integrands functions must generate rectangular contributions to the PDE system's element matrix. This is readily implemented using Diffpack's support for mixed finite elements, but the detailed explanation of the design and its implications are somewhat comprehensive and left to another manuscript [Odegard et al. 2000].
Our group's experience from implementation of numerous implicit solvers for various applications has shown that first developing a prototype solver based on a sequential solution approach and the software design herein proves to be very advantageous and productivity enhancing when coding and debugging a separate fully implicit solver.
3.3 Further Remarks
One should observe that neither Momentum1 nor Energy1 kn! ows of Manager or CommonRel and can easily be tested completel! y indepe ndently or applied in other problems where the same equations are encountered. The CommonRel class is also independent of Manager and the PDE solvers. This means that debugging of the constitutive laws also can be performed separately. Perhaps most importantly, the solvers Momentum1 and Energy1 are not over complicated with all the details of the constitutive laws. In the sketch of the class design above we have mainly focused on the numerical functionality and how this is distributed among the classes. However, a large portion of the code will often be devoted to input from a user interface and output to visualization tools. This type of code is conveniently located in the Momentum1 and Energy1 classes and completely reusable when the system of PDEs is to solved. In practice, most of the code for solving this system of PDEs is found in Momentum1, Energy1, and CommonRel. The extra glue provided by Momentum2, Energy2, and Manager is just a small amount of code, which mainly d! eals with setting up pointers and calling functionality in other classes. This type of code is therefore quite easy to debug. The main efforts of the software development for a system of PDEs is then reduced to developing and testing independent solvers and a set of constitutive relations.
We suggested to use a simple expression for [Lambda] when debugging the Momentum1 class. Before one carries out the coupling to the Manager class in a subclass of Momentum1, one can perform one or more smaller extensions of Momentum1. For example, if [Lambda] is constant in Momentum1, one can have a Momentum2 class where [Lambda] is a known spatial function and test the solver for a linear Poisson equation with variable coefficients. A further extension in terms of a subclass Momentum3 can have [Lambda] as a simple function of the primary unknown, with the purpose of testing nonlinear solver capabilities in Momentum1. Finally, a subclass Momentum4 can have communication with th! e Manager class and evaluate [Lambda] according to the physics! in the problem. For more complicated systems of PDEs such a stepwise extension of the primary solver class has proven to be very effective [Bruaset et al. 1997]. If strange errors arise later in the coupled system, it is efficient and straightforward to test the various uncoupled classes to regain reliability in each software component in the system.
Our software design has so far employed two classes for the momentum and the energy equation, although both these equations are actually on the generic form [Delta] [multiplied by] [Lambda] [Delta] u = f, if we introduce two variable coefficients [Lambda] and f. We could then easily create one class, say with name Poisson, where two methods lambda and rhs could be defined. The classes Momentum1 and Energy1 in Figure 3 would then be replaced with the Poisson class. For the thermal pipeflow problem, a common base class for the momentum and the energy equation reduces the amount of programming. However, the flexibility is also ! reduced; for efficiency reasons it can be handy to optimize the operators [Delta] [multiplied by] [Lambda] [Delta] and [[Delta].sup.2] differently, and two separate classes give us more freedom for such purposes. More complicated systems of PDEs will usually have momentum and energy equations with different structure, so in those cases one is forced to work with separate classes for each equation.
The present design relies much on the use of pointers. Admitting that "the spaghetti must go somewhere," we state that these pointers actually represent the intricate coupling of the PDEs and the constitutive relations. However, once the pointers are correctly set by the Manager class, the local code in the various classes become simple and intuitive. In practical implementation in the Diffpack system, one uses "smart pointers" (so-called Handles in Diffpack terminology), which avoid most of the programming problems normally associated with plain C or C++ pointers.
Diffpack allows the user to assign multiple values to i! tems on the input menu. The program will then enter a loop and realize all combinations of the multiple values as input data. This is very convenient for executing parameter analysis studies. The feature can be combined with automatic report generation such that it is easy to browse through the large amount of results that may be generated by multiple input data. The reports in HTML format allows plots to be inserted as well as access to visualization tools. The additional code necessary for the automatic report generation can be distributed by the various classes in Figure 3. This will be especially convenient for a system of PDEs with many physical parameters. The Momentum1, Energy1 and CommonRel can dump their data to some report class. The manager can then just call up the report-writing functionality in the solver and the common relation classes to accomplish a (big) complete report for the system of PDEs. We see that the idea of working with independent units not only simplifi! es the numerics, but also all the nonnumerical tasks that are an intrinsic part of modern PDE software.
We emphasize that no significant efficiency penalty is associated with the use of OOP techniques in our framework. Our use of OOP and advanced C++ constructions will mainly concern the program management, and not the CPU-sensitive arithmetic. An unavoidable overhead is of course associated with the extensive use of virtual functions. Since these functions cannot be "inlined" (inlining refers to copying the body of the routine into the invoking code at compile time), one has to look up the function address at each call. In our case, however, this overhead is small because of the large number of arithmetic operations relative to the number of function calls.
However, a given simulator may easily be optimized toward some target problem at hand. The point is that we first rapidly develop a prototype solver. Then, after having verified its feasibility and! performed a profiling of the program to see where most comput! ational time is spent, we optimize this solver in a subclass. Old and new code then coexist, and the debugging process of the optimized solver is hence significantly facilitated, since one can easily compare intermediate results. Moreover, the major work in the assembly process is to calculate the integrals of the weak formulation in question. In these integrals, the spatially varying part of the integrands may be written as finite-element fields and represented in terms of basis functions. Hence, the finite-element solution process can be speeded up significantly by calculating these integrals only once. The updating of the linear system to be solved then only consists of matrix-vector arithmetic that can be performed using existing, highly optimized software libraries.
The use of OOP techniques as a tool is crucial for a simple and flexible handling of the above and other efficiency improvements. With this tool at hand one can easily create "intelligent" code that merge! s coexisting modules to obtain the most efficient approach to solve a given target problem.
4. SIMULATION RESULTS
In this section we apply the simulator to investigate the volume flux Q and the convergence behavior as a function of n and [Alpha]. We only discuss the results obtained using the power-law viscosity with a temperature dependency according to (6) and the parameters as shown in Table I.
Table I. The Fixed Material Constants
[Beta] 200 Kg[(ms).sup.-2]
[[Mu].sub.0] 50 [Kgm.sup.-1][s.sup.-(2-n)]
[T.sub.0] 50 C
[Kappa] 1 J[(msK).sup.-1]
4.1 Problem Description
The physical problem we consider is flow between two nonconcentric cylinders. The outer cylinder is centered about the origin and has radius [R.sub.o] = 1. The inner cylinder is centered about (0.3, 0), and it has radius [R.sub.i] = 0.5.
The following boundary con! ditions are prescribed. The velocity vanishes at all boundarie! s ([diff erential][[Omega].sub.1]). For the temperature, the inner cylinder is marked as [differential][[Omega].sub.5] with [T.sub.1] = 100; the upper half of the outer cylinder is marked as [differential][[Omega].sub.4], i.e., T = 0; and the lower half of the outer cylinder is assumed to be insulated, i.e., marked as [differential][[Omega].sub.3].
Due to the [Gamma] variable in the equations, it is important to provide a reasonable initial guess for the velocity. In the present case we assume that at each point the velocity is that of a Poiseuille flow between two flat planes distanced apart with the normal distance from the inner cylinder to the outer cylinder. The temperature is initially zero. In cases where the continuation method is used, we have fixed the continuation parameters as [Lambda] = {0.1, 0.5, 0.75, 1}, but in the case where convergence is not reached with a specific [[Lambda].sub.i], another attempt is made using [[Lambda].sub.i] = [[Lambda].sub.i] + [[La! mbda].sub.i-1]/2.
In the present case we have chosen the convergence criterion as
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]
with [Epsilon] = 1.0 x [10.sup.-4]. In the cases where we switch between the Picard and Newton iterations, this is performed when [Epsilon] = 1.0 x [10.sup.-2]. For the implicit solvers FI-NR and FI-P, we run as many iterations that are necessary to obtain convergence. The same is true for the P1 solver. For the other explicit solvers, we perform as many global iterations as necessary, but the number of internal iterations in both the momentum equation and the energy equation is limited to 30 at each global iteration.
4.2 Discussion of the Results
Figures 6 and 7 compares the execution times and the volume flux as a function of n for some different sequential solvers. As expected, the Gauss-Seidel iterations performs better than the Jacobi iterations for challenging values for [Alpha] and! n, i.e., when temperature effects increases and the interacti! on betwe en the equations becomes more significant. As seen from Figure 7 (right) the use of Gauss-Seidel iterations leads to shorter execution times and makes it feasible to obtain convergence for larger values of n for [Alpha] = 0.01. For the sequential solver, the use of Newton iterations is preferable in the momentum equation, whereas the choice of scheme in the energy equation is not important.
[GRAPHS OMITTED]
It is interesting to note, that for the fully implicit schemes, the Picard iterations surprisingly seem to be more appropriate than the Newton iterations, especially for extreme values of n and for all nonzero values of [Alpha]. The equation system arising with the FI-P scheme is better conditioned than the one obtained with the FI-NR scheme. The Conjugate Gradient method could not be used in the FI-NR scheme for n [is less than] 1, but instead a BiCGStab iteration technique had to be used. This leads to an increase in the execution times (see Figure! 9, left). As the temperature effects increases, this problem becomes more expressed, and Newton iteration fails. In the case of [Alpha] [is greater than or equal to] 0.01 using the FI-NR scheme, the equation system becomes singular when approaching the solution, and this scheme hence does not converge even for n = 1. The FI-P scheme, however, which is algebraically equivalent to the J-PI scheme, converges in this case for 0.2 [is less than or equal to] n [is less than or equal to] 1.7. Switching between the Picard and Newton iterations has no significant effect in the present case.
[GRAPH OMITTED]
Figures 10 and 11 show the results obtained in case where n = 0.2 and [Alpha] = 0 and 1. As seen from these figures and the above Volume flux results, the effect of the temperature is significant when [Alpha] = 0.01.
[ILLUSTRATIONS OMITTED]
5. INCREASING THE PROBLEM COMPLEXITY
The particular model problem chosen for illustrat! ion of our proposed implementation technique is relatively sim! ple from a mathematical point of view. In most real-world engineering applications, however, one encounters more challenging problems. For example, assume that we will solve the full 3D case of the equation system (6). We then face a problem that is significantly more complex, including the need to solve incompressible Navier-Stokes-like equations with the associated problems [Chorin 1967; Reddy 1982; Taylor and Hood 1973; Zienkiewicz and Wu 1991].
Figure 12 shows a suggested class hierarchy for dealing with this case. The Manager class is unaltered; besides that it now contains pointers to the more general abstractions FlowSolver and EnergySolver for solving the momenta and energy equations, respectively. The EnergySolver class is now derived from a class ConvDiff that is similar to the Energy2 class, but deals with a full convection-diffusion equation instead of a Poisson equation. The actual flow solver PredCorr, derived from the FlowSolver, however, is somewhat more c! omplicated. Since solving the momenta equations now consists of solving either four scalar equations or one vector and one scalar equation, the PredCorr class cooperates with the Manager class as some kind of a submanager that is responsible for this part of the solution process.
[ILLUSTRATION OMITTED]
We remark that while we in the present case have assumed that the flow equations are solved using a predictor-corrector strategy in a similar manner as in Chorin [1967], the class FlowSolver may equally well represent any other solution strategy or any other efficiently implemented special case similar to the one presented in this paper. The use of OOP lets the user test and replace software components in a flexible manner without affecting any other parts of the code. This gives unique possibilities to plug and play with different solution approaches in a given challenging problem or to rapidly modify the simulator to solve special-case problems in an e! fficient manner.
In general, OOP allows us, in a si! mple fas hion, to build our software in terms of independent modules that can be easily combined to create sophisticated simulators for a wide variety of real-world physical problems. Interesting future applications that may be coupled through the CommonRel class will be turbulence or multiphase simulation. Since such modeling will affect the equation parameters, this should be accomplished in a subclass, symbolized as SomeRel in Figure 12, in the CommonRel-hierarchy. This class may include pointers to data structures that in a similar manner as the FlowSolver class cooperates with the CommonRel class as another submanager for turbulence and other modeling.
6. CONCLUDING REMARKS
This paper has proposed a strategy for the development of computer code for solving systems of nonlinear partial differential equations using object-oriented programming techniques. While the implementation technique, designed for rapid building of simulators for systems of partial diffe! rential equations by merging together independent solvers for alone-standing equations, has shown useful also in the present coupled heat-fluid flow problem, the benefit of a well-structured, modular implementation will become more expressed as the problem complexity grows.
Although the present paper has focused more at human efficiency than on pure computational efficiency, this issue becomes important when the problem size increases. Relative to a highly optimized simulator, and looking at computational efficiency exclusively, some degree of performance loss is necessary to achieve the necessary flexibility. As discussed, however, various steps can be taken to maintain a high degree of computational efficiency in spite of the flexibility that is provided. Also, one can easily trade flexibility and generality for computational efficiency by a special design toward a given application.
[GRAPH OMITTED]
(1) World Wide Web home page: http://www! .diffpack.com/
REFERENCES
AHMADI, G. 198! 5. A gen eralized continuum theory for multiphase suspension flows. Int. J. Engng. Sci. 23, 1.
ARGE, E., BRUASET, A. M., AND LANGTANGEN, H. P., EDS. 1997a. Modern Software Tools for Scientific Computing. Birkhauser Boston Inc., Cambridge, MA.
ARGE, E., BRAUSET, A. M., AND LANGTANGEN, H. P. 1997b. Object-oriented numerics. In Mathematical Models and Software Tools in Industrial Mathematics, M. Daehlen and A. Tveito, Eds. Birkhauser, Basel.
BRAUSET, A. M., HOLM, E. J., AND LANGTANGEN, H. P. 1997. Increasing the reliability and efficiency of numerical software development. In Modern Software Tools for Scientific Computing, E. Arge, A. M. Bruaset, and H. P. Langtangen, Eds. Birkhauser Boston Inc., Cambridge, MA.
CHORIN, A. 1967. A numerical method for solving incompressible viscous flow problems. J. Comp. Phys. 2, 12-26.
GUPTA, G. AND MASSOUDI, M. 1993. Flow of a generalized 2nd-grade fluid between heated plates. Acta Mech. 99, 21-23! .
HEYES, D. M. AND MARCH, N. H. 1994. Mechanical properties of liquids-Newtonian and beyond. Phys. Chem. Liq. 28, 1, 1-27.
JOHNSON, G., MASSOUDI, M., AND RAJAGOPAL, K. R. 1991. Flow of a fluid infused with solid particles through a pipe. Int. J. Eng. Sci. 29, 649.
JOHNSON, G., MASSOUDI, M., AND RAJAGOPAL, K. R. 1991. Flow of a fluid-solid mixture between flat plates. Chem. Eng. Sci. 46, 1713.
LANGTANGEN, H. P. 1999. Computational Partial Differential Equations: Numerical Methods and Diffpack Programming. Springer Lecture Notes in Computational Science and Engineering, vol. 2. Springer-Verlag, Vienna, Austria.
MASSOUDI, M. AND CHRISTIE, I. 1993. Heat transfer and flow of a third grade fluid in a pipe. Math. Model. Sci. Comput. 2, 1275.
MASSOUDI, M. AND CHRISTIE, I. 1995. Effects of variable viscosity and viscous dissipation on the flow of a 3rd grade fluid in a pipe. Int. J. Non-Lin. Mech. 30, 5, 687-699.
NORTON, C.D. 1998. Object-oriented programming paradigm! s in sci entific computing. Ph.D. Dissertation. Computer Science Department, Rensselaer Polytechnic Institute, Troy, NY.
ODEGARD, A., LANGTANGEN, H. P., AND TVEITO, A. 2000. Designing implicit solvers for PDE systems from independent component solvers. In Computational Partial Differential Equations: Advanced Methods and Applications, H. P. Langtangen and A. Tveito, Eds. Springer-Verlag, New York, NY.
PAPACHRISTODOULOU, G. AND TRASS, O. 1987. Coal slurry fuel technology. Can. J. Chem. Engng. 65, 177.
RAMKISSOON, H. AND EASWARAN, C. 1993. Non-newtonian flow between concentric cylinders. Zeit. Ange. Math. Mech. 73, 329-332.
REDDY, J. N. 1982. On penalty function methods in the finite element analysis of flow problems. Int. J. Num. Meth. Flu. 18, 853-870.
SAMUELSSON, K. 1996. Adaptive algorithms for finite element methods approximating flow problems. Ph.D. Dissertation.
SCHOWALTER, R. 1978. Mechanics of Non-Newtonian Flui! ds. Pergamon Press, Inc., Tarrytown, NY.
SKELLAND, A. H. P. 1969. Non-Newtonian Flow and Heat Transfer. John Wiley and Sons, Inc., New York, NY.
SZERI, A. AND RAJAGOPAL, K. 1985. Flow of a non-Newtonian fluid between heated parallel plates. Int. J. Non-Lin. Mech. 20.
TAYLOR, C. AND HOOD, P. 1973. A numerical solution of the Navier-Stokes equations using the finite element method. J. Comp. Phys. 1, 73-100.
WINTHER, H. 1977. Viscous dissipation in shear flows of molten polymers. Adv. Heat Trans. 13, 205.
ZIENKIEWICZ, O. C. AND WU, J. 1991. Incompressibility without tears: How to avoid restrictions of mixed formulation. Int. J. Num. Meth. Flu. 32, 1189-1203.
Received: January 1999; revised: May 2000; accepted: July 2000
This work has received support from The Research Council of Norway (NFR) under Grant 110673/420 (Numerical Computations in Applied Mathematics).
Authors' address: Dept. of M! athematics, University of Oslo, P.O. Box 1053, Blindern, N
B&B Italia via Marcus Design
warymeyers via Alkemie 
via Brume
