Problem solving for scientific computing: data modelling instead of algorithms

Problem solving for scientic computing:

data modelling instead of algorithms?

Victor Ganzha, Dmytro Chibisov

1

and Evgenii Vorozhtsov

2 1 Institute of Informatics, Technical University of Munich, Munich, Germany

e-mail: ganzha@informatik.tu-muenchen.de

e-mail: chibisov@informatik.tu-muenchen.de

2 Institute of Theoretical and Applied Mechanics, SB RAS, 630090 Novosibirsk,

Russia

e-mail: vorozh@itam.nsc.ru

Received: December 20, 2001

Summary.

We present an approach to rapid solving of modelling problems for uid dynamics. We show how the problem solving pro-cess in uid dynamics can be expressed through the incremental suc-cessive data modelling on the basis of graphical description tech-niques like class diagrams in UML. As will be shown, in contrast to algorithmic approach the data modelling in terms of objects and relations among them enable one to describe the physical and math-ematical properties and relations of the uid ow. At the same time, it provides the necessary facilities to obtain the numerical code au-tomatically from an abstract problem description. Moreover the tool prototype called GROOME to numerical code generation on the ba-sis of graphical object oriented description techniques is presented. In order to demonstrate these facilities we consider a uid dynamics problem involving shock waves in two-dimensional inviscid gas ow.

Key words:

Euler equations, nite volume method, object-oriented modelling, problem solving environment

and boundary-value problems for PDE to ecient, documented and executable code, which is typically generated in C or Fortran. But the development of new numerical methods is not simply an automatable process of transforming one method into another it in-volves complex synthesis and analysis tasks in order to understand the multilevel relationships between dierent objects in the problem domain.

This paper presents an approach to develop the numerical methods based on graphical data modelling techniques. One of such techniques called "Entity-Relationship-Modell" was proposed by Peter Chen in 1976 1]. Another example is the well known Unied Modelling Lan-guage (UML).

We will show how these techniques can eciently be used to solve modelling problemsfor uid dynamics. As will be shown in the present paper, the traditional methods used by mathematical modelling sep-arate mathematical foundations of a particular numerical method from the resulting computer code. Such approach leads to redundant work in order to translate a mathematical model into software. Our approach in contrast to them is an integration of both aspects of scientic computing in the form of a single description technique to describe at the same time mathematical models, data structures and algorithms. This technique is a combination of graphical object ori-ented diagrams like class diagrams in UML and symbolic expressions, for example, formulas. As a result, the computer code for a particular problem to solve can be obtained automatically.

Moreover, we propose some techniques to get over the well known limitations of the graphical data modelling techniques. For example, we will show how the notion of periodicity and symmetry can be modelled through UML like class diagrams.

Modelling Environment) supporting these features has been present-ed in 7]. In 7] we have shown how the package GROOME functions at a relatively simple example of the one-dimensional nite element method.

In the present paper, we present an extension of the general GROOME methodology for computationally much more complex two-dimensio-nal uid dynamics problems. In this extension, the problem solving is considered as an

incremental

successive process from a problem specication in the form of conservation laws that are widely used in uid dynamics to the nal stage { numerical computer code.

Our approach is motivated by published theories of the cognitive challenges of software development 11] and, of course, by the nature of mathematical models of uid dynamics tasks.

1.1. From Problem Specication to Executable Code

In this section we describe some basic ideas of Problem Solving En-vironments generally. The termProblem Solving Environment (PSE) has dierent meanings to dierent people. According to Gallopoulos, Houstis and Rice (4]), a problem solving environment can generally be characterized as

A PSE is a computer system that provides all the computa-tional facilities necessary to solve a target class of problems. These features include advanced solution methods, automatic and semiautomatic selection of solution methods, and ways to easily incorporate novel solution methods. Moreover, PSEs use the language of the target class of problems, so user can run them without specialized knowledge of the underlying computer hardware or software. Overall, they create a framework that is all things to all people: they solve simple or complex problems, support rapid prototyping or detailed analysis,and can be used in introductory education or at the frontiers of science.

In the next section we will discuss some psychological theories of problem solving and implications for the tools supporting problem solving.

theory.

For example, by ow modelling the problem situation model consists of knowledge about the physical eects of the ow and the equations that describe it. The solutions model consists of knowledge about the discretization methods, data structures, algorithms, syntax and semantics of a programming language needed to implement the solu-tion of the discretized equasolu-tions 5].

Implications for PSE In order to help the developer to map his mental problem situation model onto the solution model in the form of computer program we propose the integration of both traditional mathematical notations to describe laws and behaviour of gas ow and graphical description techniques like UML to describe software structure. Let us explain it with the aid of an example.

3. Equations to Solve

Consider the Euler equation in the following conservation form:

(1) @w @t +@f(w) @x + @f(w) @y = 0

wherex and y are Cartesian coordinates and w= 0 B B @  u v E 1 C C A  f(w) = 0 B B @ u u 2+ p uv uH 1 C C A  g(w) = 0 B B @ v vu v 2+ p vH 1 C C A :

q

q

_q

q

q

jk;1 j;1k jk j+ 1k jk+ 1     1 2 3 4

Fig.1. The centered scheme

Here puvE and H denote the pressure, density, Cartesian

ve-locity components, total energy and total enthalpy. For a perfect gas

E = p (;1)+ 12( u 2+ v 2)  H=E+ p  

where is the ratio of specic heats.

The above Euler equation can be written in integral form for a given volumeV as follows: @ @t ZZ V wdxdy=; I ; (fdy;gdx):

What it means is that the change of w in the volume V equals the

ux ofw through its contour.

3.1. Discretization

Let us take an arbitrary cell of curvilinear grid (Fig. 1). Let the values ofw at the cell center be denoted byw

jk, and let ;

jk and V

jk be the

cell contour and the control volume bounded by contour ;

jk. The

Euler equations (1) can be written in integral form for the regionV jk with boundary ; jk as (2) @ @t ZZ V jk wdxdy+ I ; jk (fdy;gdx) = 0:

The above notation is used by mathematicianssince hundreds of years but in order to construct computer programs it is not precise enough.

Fig. 2. Class Diagram shows two relationships betweenV and Vjk: 1)V

con-sists of 80 times 40 Vjk 2) EachVjk inherits the physical properties, such as

conservation of mass, momentum and energy, fromV

For instance, the relationship "is part of" between the control volume

V

jk and volume

V is not implicit in the formula (2). Furthermore,

the information about relationships between the coordinates of V jk

needed for the evaluation of integrals is missed.

In order to deal with such data incompleteness we propose to use the data modelling techniques on the graphical basis like class diagrams in UML.

4. Graphical Object Oriented Problem

Domain Description

As we have seen above, the traditional notation says us nothing about the relationships between mathematical objects. Consider, for exam-ple, the control volume V

jk and volume

relation-Fig. 3. Class Diagram representing the continuity of a ow as relationship

be-tweenVjk and;jk

ships among them. From topological point of view, V

consists of

a nite set fV jk

g ordered in space, for example i = 1:::80j =

1:::40. From physical viewpoint each nite volume V

jk has the

same physical properties as V, for example, it conserves mass,

mo-mentum, and energy. What it means in terms of object oriented mod-elling is that V

jk

inherits

its properties from V.

Such facts can be expressed with the aid of a UML class diagram as shown in Fig. 2.

The arrow means that the volume "V(j,k)" inherits the conservation properties of "V". The second relation means that the volume "V" consists of 80 times 40 nite volumes "V(j,k)".

Let us note the conservation properties in terms of object oriented relations. Equation (2) means that the change of w in the volume V

jk is equal to the ux of

wthrough the contour ;

jk. This relation

is shown in Figure 3. It shows how the attributewof the objectV jk

depends on the attributesf(w) and g(w) of the object; jk.

In the result of the discretization of (2) the user can obtain the semi-discrete equation 9,8,6]: (3) @A jk w @t +Qw= 0 where A

jk is the cell area, and the operator

Q represents an

ap-proximation to the boundary integral dened by the second term of (2).

Fig. 4. Static class diagram to represent the relations between a cell V jk, its

boundary;jkand its sides

Fig.5. Parameterized static class diagram shows the periodical grid structure

For example, the ux balance for the x momentum component is

represented in (3) as (4) @ @t (A jk u) + 4 X (Q k u k+ y k p k) = 0 

Fig.6. Flux betweenV(jk) andV(j+ 1k) as relation between the attributes

of the appropriate classes

where the ux velocity

(5) Q k = y k u k ;x k v k

and the sum in (4) is over four sides of the cell, see Fig. 1. The values

x k and

y

k in (4) and (5) are the increments of

xand yalong side k of the cell, with appropriate signs. Each quantity in (4) and (5)

such as u 2 or (

u)

2 is evaluated as the average of the values in the

cells on the two sides of the face, (u) 2= 12( u) jk+ ( u) j+1k] :

We have mentioned above that we consider the problem solving such as above discretization as incremental modifying of the relations be-tween objects through the better understanding of them. For in-stance, to compute the uxes Q we need the values of w on each

cell side. Therefore, we introduce a new object "Side" (Fig. 4). The relation between "V(jk)" and "Side" expresses that each "Side"

be-longs to 1 (in case of boundary sides) or 2 (in case of inferior sides) nite volumes "V(jk)" and in other direction each "V(jk)" consists

of 4 objects "Side".

Traditional graphical data modelling techniques, for example, Entity-Relationship-Diagrams or Static Class Diagrams in UML (except for Object Constrain Language, an additional textual notation included in UML) do not provide any possibility to express such case depen-dent relations. Another example that can not be expressed is the simple symmetry: a right hand side of the cellV(jk) is a left hand

side of the cell V(jk+ 1).

To come over this limitationwe introduce the

parameterized

or

Fig.7. Flux through theV(jk) contour represented as relation betweenV(jk), V(j;1k),V(j+ 1k),V(jk;1),V(jk+ 1)

means that the data structure can be described depending on one or several parameters, for example, spatial coordinates. The periodical or parameterized phenomena are often to nd in the nature (waves, crystal structures, etc.). For example, consider the object oriented description of the grid in Fig. 5.

The ux relation between two neighbour control volumes can now be expressed as relation between the attributes of classes as shown in Fig. 6 for V(jk) and V(j+ 1k). However, the ux relation can

be expressed as relation between V(jk), V(j + 1k), V(j ;1k), V(jk+1),V(jk;1) only. It means that objects "Horizontal Side"

and ";

jk can be eliminated (Fig. 7). The time stepping scheme can

be considered as relation between the ow-objects at dierent time points.

Let us construct it as an explicit three-stage scheme of the Runge-Kutta type: (6) w (0)= w n w (1)= w (0) ; 1 (P h w (0) ;D w (0)) w (2)= w (0) ; 2 (P h w (1) ;D w (0)) w (3)= w (0) ; (P h w (2) ;D w (0)) 

Fig.8. Artical dissipation as relationships between objects where w n = w(xyt n) t n = t n;1+ n  n = 12:::, t 0 = 0  n =  P h w = (1=A jk) Qw jk  and  1 

2 are the nondimensional weight

parameters.

The above Jameson scheme reduces to a central dierence scheme on a Cartesian grid, therefore, it is second order accurate in space. As is known, the second-order schemes produce spurious oscillations of the numerical solution at the shock wave fronts and in their vicinity. In this connection it was proposed in 9] to introduce the articial dissipation terms in the nite volume scheme in order to damp the spurious oscillations. These terms are added to equation (3) as fol-lows: dw dt + (Qw =A jk) ;D w= 0

whereDis a dissipative operator the structure of which is similar for

each of the four dependent variables:

(7) D w jk = D x w jk+ D y w jk  where (8) D x w jk= (1 =A jk)(

d

j+1=2k ;

d

j;1=2k) 

(11) jk = j+1k 2 jk+ j;1k ( j+1k + 2 jk + j;1k ) Then (12) " (2) j+1=2k = (2)max( j j+1k jj jk j) " (4) j+1=2k = max(0 ( (4) ;" (2) j+1=2k)) :

Typical values of the constants (2)and (4)are 9] (2)= 1 =4 (4)= 1=256. The values

d

jk 1=2in (7) are computed similarly to (8){(10)

for example,A jk +1=2= (1 =2)(A jk+ A jk +1).

As can be seen in (7){(12), the dissipative term can be represented as relationship between the appropriate quantities ofV

jk, V j+1k, V j+2k, V j;1k, V jk +1, V jk +2, V

jk ;2. Let us describe this relationship using

the diagram in Fig. 8.

The Runge-Kutta scheme can be considered as a relationship between the same object V

jk at two dierent time points (Fig. 9).

Now we can substitute the vector consisting ofuvpandE instead

of wand appropriatef(w) as shown in Fig. 10.

As a more specic example we consider the shock waves in a gas and demonstrate how the boundary and initial conditions can be modelled using the above technique.

5. Shock Waves in Gas

As a ow problem we have chosen a simple problem of inviscid ow developed by an oblique stationary shock wave reecting from a rigid surface (Fig. 11). This test problem is often used for the validation of new numerical methods in computational uid dynamics. The advan-tage of this test is that it is possible to obtain the exact solution for it

Fig.9. Runge-Kutta relationship between the same object at two dierent time

points

Fig. 10. Class box

with attributes

by using the theory of stationary oblique shocks. This solution repre-sents a piecewise constant function. We have used the value'==6

for the angle between the incident shock wave front and the x axis

con-q

q

q

k=;1 k= 0 k= 1

Image cells for symmetry techniques

stants of the exact solution in subregions 1, 2, and 3 indicated in Fig. 11 are as follows:

Subregion 1 Subregion 2 Subregion 3

u 1= 1 :0 u 2 = 0 :890755053 u 3= 0 :806645743 v 1 = 0 :0 v 2= ;0:189217798 v 3 = 0 :0 p 1= 0 :084932903 p 2 = 0 :194177850 p 3= 0 :390838939  1= 1 :0  2 = 1 :776135164  3 = 2 :898621574 M 1 = 2 :90 M 2= 2 :327642861 M 3 = 1 :856588584 HereM 1 M 2 M

3are the values of the Mach number M = p u 2+ v 2 =c

in subregions 1, 2, and 3, respectively c is the sound velocity, c = p

p=. The angle between the reected shock wave front and the

positive direction of the x-axis is equal to' 3 = 0

:418279545.

Initial conditions.Initially, the entire ow eld is set equal to the free stream supersonic inow values given above for subregion 1, that is the initial gas ow is parallel with thex-axis.

Analytical boundary conditions. The spatial region has the size 4 along thex-axis and the size 1 along the y-axis. The boundary conditions

are given as follows:

(a) supersonic inow atx= 0, 0y1, which allows the values of uvpto be xed at free stream conditions given above as subregion

1

(c) supersonic outow atx= 4, 0y1

(d) a rigid at surface at y = 0, 0  x  4, which is properly

represented by the condition v = 0 with the additional condition @p=@y= 0 aty= 0 from the normaly-momentum equation.

Numerical boundary conditions.The supersonic outow values w Jk, k= 1:::K, are obtained by zeroth-order extrapolation, i.e.,w

Jk = w

J;1k,

k= 1:::K.

As it follows from formulas (7){(12), for the computation of the ar-ticial dissipation terms in the Jameson's method, the numerical so-lution values are needed in image cells adhering to the spatial region boundaries (see Fig. 12). We have specied the values of the compo-nents of the solution vectorw in image cells in accordance with the

symmetry technique presented in 3]:

Rigid wall Outow boundary x= 4 p j0= p j1 p j;1= p j2 w Jk = w J;1k  j0=  j1  j;1=  j2 w J+1k = w J;2k u j0= u j1 u j;1= u j2 v j0= ;v j1 v j;1= ;v j2

The numerical solution values in two rows of cells adhering to the inow boundaries x = 0 and y = 1 were specied similarly to the

case of the outow boundary x= 4.

In order to specify the boundary and initialconditions using graphical data modelling let us introduce the classes "V(1,k,0)", "V(80,k,0)", "V(i,1,0)" and "V(i,40,0)" corresponding to each of four sides of the considered region at t= 0. These classes are derived from the class

"V(i,j,t)" as shown in Fig. 13.

The inow values on the left region side correspond to object "V(1,k,0)" and can be specied as shown in Fig. 14.

Fig.14. The specication

of inow parameters

6. GROOME - Tool Supported Graphical

Object Oriented Data Modelling

We have seen above how the computational methods can be rep-resented in the declarative way as objects and relationships among them. In order to create and edit such diagrams we have developed a tool prototype called GROOME 7]. GROOME supports the user by creating and editing a hierarchy of such graphical documents that enable a computational problem to be considered at dierent abstrac-tions levels. Furthermore, GROOME enables one to link graphical documents together and generates numerical code automatically.

The above presented graphical description of objects in conjunction with symbolic given relations among their attributes allows the nu-merical code for the Jameson scheme to be generated automatically. Consider, for example, the fragments of the generated class "Vjkt" for inferior control volumes. Note that this Maple code uses our package that provides object oriented programming in CAS Maple (2]). The class Vjkt below needs as parameter for instantiation8 neighbour volumes in the space (as explained above we need 4 neighbour vol-umes to calculate the uxes (Vkt1t- Vkt4t) plus 4 volumes (Vkt5t -Vkt8t) to calculate the dissipative term) and 1 neighbour volume in

the time (Vkt9t): Vjkt:=proc(j,k,t,Vjkt1,Vjkt2,Vjkt3,Vjkt4,Vjkt5,Vjkt6, Vjkt7,Vjkt8,Vkt9t) local this global alpha,dt: this:=NEW_INSTANCE() #************************************************* # ro0,v0,u0,p0,E0 as computational results of the* # previous time step (object Vkt9t) * #************************************************* obj||this||ro0:=obj||(obj||this||Vjkt9)||ro3() obj||this||v0:=obj||(obj||this||Vjkt9)||v3() ...

# ********************************************** # computaton of ro(i) at the appropriate * # stage of Runge-Kutta method (i=1..3) * # ********************************************** obj||this||ro1:=proc() use oop in RETUNR(obj||this||ro0 - alpha*dt*(obj||this||Qjk0()/obj||this||Ajk()-obj||this||Djk())) end use: end proc:

end use: end proc:

# ***************************************************** # computaton of u(i),v(i),E(i).p(i) at the * # appropriate stage of Runge-Kutta method (i=1..3) * # ***************************************************** obj||this||v1:=proc() use oop in obj||this||v:=(obj||this||ro0)*(obj||this||v0)-alpha*dt* (obj||this||Qjk0/obj||this||Ajk()-obj||this||Djk() end use: end proc: ... obj||this||u1:=proc() use oop in ... end use: end proc: ... # ************************************************* # computaton of ro fluxes roQjkt(i),uQjkt(i), * # vQjkt(i), EQjkt(i) at the appropriate stage of * # Runge-Kutta method (i=1..3) * # *************************************************

use oop in RETURN ((obj||this||ro0*obj||this||v0 + obj||this||(obj||this||Vjkt1)||ro0)/2 -(obj||this||ro0*obj||this||v1 + obj||this||(obj||this||Vjkt2)||ro0)/2+ (obj||this||ro0*obj||this||v1 + obj||this||(obj||this||Vjkt3)||ro0)/2 -(obj||this||ro0*obj||this||v1 + obj||this||(obj||this||Vjkt4)||ro0)/2) end use: end proc: ... # ********************************************** # computaton of cell area Ajkt * # **********************************************

...

RETURN(this) end proc:

The classes for boundary volumes will be generated in the same way according to the above described symmetry techniques.

0.1 0.15 0.2 0.25 0.3 0.35 0.4 P 0 1 2 3 4 x 0 0.2 0.4 0.6 0.8 1 0 1 2 3 4 (a) (b)

Fig. 15. Shock reection problem: (a) pressure prole in the section y= 0:45,

(|) the exact solution, () the numerical solution (b) predicted Mach number

contours.

We show in Fig. 15 some computational results obtained by the above described Maple OOP code. These results were obtained on a rela-tively crude mesh of 4111 grid lines (4010 cells). It may be seen

1093{1108.

4. Gallopoulos, E., Houstis, E. and Rice, J.R. (1992)q: Future Research Direc-tions in Problem Solving Environments for Computational Science, Report of a Workshop on Research Directions in Integrating Numerical Analysis, Symbolic Computing, Computational Geometry, and Artical Intelligence for Computational Science. (Technical Report 1259, Center for Supercomputing Research and Development, University of Illinois at Urbana-Champaign). 5. Ganzha, V.G. and Vorozhtsov, E.V. (2000): Diagram design and analysis of

aerodynamics problems with Mathematica, Selcuk J. Appl. Math.1, 21{46.

6. Ganzha, V.G. and Vorozhtsov, E.V. (2001): Stability investigation of Runge{ Kutta schemes with articial dissipator on curvilinear grids for the Euler equations, Math. Comp. Simulat.58, 1{36.

7. Ganzha, V.G., Chibisov, D. and Vorozhtsov, E.V. (2001): GROOME { tool supported graphical object oriented modelling for computer algebra and scien-tic computing, in: Computer Algebra in Scienscien-tic Computing/ CASC 2001, V.G. Ganzha, E.W. Mayr and E.V. Vorozhtsov, eds., Springer-Verlag, Berlin, 213{232.

8. Jameson, A. and Schmidt, W. (1985): Some recent developments in numerical methods for transonic ows, Comp. Meth. Appl. Mech. Eng.51, 467{493.

9. Jameson, A., Schmidt, W. and Turkel, E. (1981):Numerical Solution of the Euler Equations by Finite Volume Methods Using Runge{Kutta Time Step-ping Schemes, AIAA Paper 81-1259.

10. Kintsch, W. and Greeno, J.G. (1995): Understanding and solving word arith-metic problems, Psychological Review.92, 109{129.

11. Robbins, J. E. (1999):Cognitive Support Features for Software Development Tools, Dissertation. University of California, Irvine.

12. Steinberg, S., Akers, R.L., Kant, E., Randall, C. J. and Young, R.L. (1997): SciNapse: A problem-solving environment for partial dierential equations, IEEE Comp. Sci. Eng.4, 32{42.