Home

Rudolf Beck Bodo Erdmann Rainer Roitzsch KASKADE User's Guide

1. Figure 24 Final triangulations for the first and the last time step in example 2 There are further applications in metallurgy soil mechanics decision and control theory etc For more details we refer to Crank Cr88 Elliot and Ockendon EO82 Hoffmann and Sprekels HS90 Kornhuber Kor95 and the literature cited therein For discretization in time we use a backward Euler scheme with a uniform step size The spatial problems at the different time levels are similar to those solved in the obstacle problem 1 In figure 24 you see the final triangulations for the first and the last time step 3 Nonlinear Problem using the files pmedia cmd pmedia geo pmedia mat In this example we consider the degenerate parabolic initial boundary value problem 00 0t A0 0 0920 m gt 1 in Q x 0 7 0 0 6 in Q 0 g on Ep x 0 7 00 0n 0 on Ty x 0 7 4 PROBLEM CLASSES 3T of type 6 in appendix A describing the adiabatic flow of a homogeneous gas through a porous medium Mus37 This example in two spac
2. EE T N D u l Se ES 7 Figure 25 Final mesh and isolines of solution after 1 20 and 40 time steps in example 3 5 GEOMETRY AND MATERIAL PROPERTIES 39 5 Geometry and Material Properties For a complete definition of a problem we have to specify the type of the partial differential equation see preceding chapter the space dimension of the domain a triangulation of the domain and the coefficients e g k and f in appendix A 1 which usually describe material properties 5 1 Definition of the Geometry The space dimension is selected by the command spaceDim Here stands for 1 2 or 3 The default value is 2 The triangulation of the domain 2 is given in a file defining a set of intervals 1D triangles 2D or tetrahedra 3D We use the extension geo for this file The filename is specified by the command file We offer two formats to code the geometrical data By the command read0ldFormat 1 you select the format we used in former versions up to 2 x of KASKADE see Roitzsch et al Roi93 We recommend as default a simpler format readZIBFormat 1 which we describe here e 1D r
3. Figure 14 Mesh and solution after seven refinement steps in example 12 3 STARTING THE PROGRAM 27 14 Static problem using the files 1ayer 3d cmd layer 3d geo layer 3d mat Here we solve the analogous equation to Problem 13 in three space dimensions lO0 uz Uyy uz 1 in 0 25 0 75 x 0 25 0 75 x 0 25 0 75 Use Uyy cuu 1 elsewhere in the unit cube u 10 0 on the boundary By the command innerBoundary 1 we additionally integrate the term k ko x y z on the layer between the two different materials which must be marked by type 1 and 2 in the mat file k 1 and kz 100 denote the material coefficients Static problem using the files flow 2d cmd flow 2d geo flow 2d mat The same equation and region as in Problem 5 but with Neumann boundary conditions on each boundary edge Only on two points we prescribe Dirichlet values Ugg tty 1 in 0 1 x 0 1 Ou On 0 0 on the boundary u 10 in 0 0 u 50 in 1 0 Figure 15 shows the initial coarse mesh and the mesh after some adaptive refi nement steps Static problem using the files flow 2d a cmd flow 2d a geo flow 2d a mat Here we solve an elliptic equation with a mass term and only Neumann condition on the boundary Urr dap ee 1 on a circle of radius v2 Ou On 1 0 on the boundary Figure 16 shows the mesh and the isolines of the approximate solution after so
4. Konrad Zuse Zentrum f r Informationstechnik Berlin Heilbronner Str 10 D 10711 Berlin Wilmersdorf Rudolf Beck Bodo Erdmann Rainer Roitzsch KASKADE User s Guide Version 3 x Technical Report TR 95 11 December 1995 KASKADE User s Guide Version 3 x Rudolf Beck Bodo Erdmann Rainer Roitzsch Abstract KASKADE 3 x was developed for the solution of partial differential equati ons in one two or three space dimensions Its object oriented implementation concept is based on the programming language C Adaptive finite element techniques are employed to provide solution procedures of optimal computatio nal complexity This implies a posteriori error estimation local mesh refinement and multilevel preconditioning The program was designed both as a platform for further developments of ad aptive multilevel codes and as a tool to tackle practical problems Up to now we have implemented scalar problem types like stationary or transient heat con duction The latter one is solved with the Rothe method enabling adaptivity both in space and time Some nonlinear phenomena like obstacle problems or two phase Stefan problems are incorporated as well Extensions to vector valued functions and complex arithmetic are provided This report helps to work with KASKADE Especially we study a set of examples explain how to define a user s problem and introduce a graphical user interface
5. cube with faces marked to approximate curvilinear faces of a cylinder during the refinement process points maxIndex 10 10 1 11 2 111 111 1 1 1 LoI 8 zi b 1 1 1 de FE NoN AOU END elements 10 2 3 1035 3 56 1034 4 57 3 45 END NOON DD ARA boundary faces 10 10 10 2 3 5 10 10 4 4 3 5 END C or c marks a curvilinear face 2000 PRPRPrPrPRPRPRP PrP RE ann qm NON OKRWAAWOANDND ONWDANKRNNADDHDW 2 oongnononononononnununmn As an example look at figure 32 52 5 GEOMETRY AND MATERIAL PROPERTIES 53 Figure 32 Approximation of a cylindrical geometry during 3 uniform refinement steps In the future we will realize a more general concept for curvilinear boundaries 5 2 Ouput of Geometrical Data We offer a set of formats to store the geometrical data of the actual refinement level into a file The user can select one of the following formats e writeZIBFormat This format recommended by the ZIB is described in the preceding section e writeGRAPEFormat This format can be used by the graphical visualization software GRAPE e writeAVSFormat This format is evaluable by the graphical visualization software AVS The output is initiated by the command writesolution 1 In additon to the geometrical data we write the values of the approximate solution in the nodes into the file A file generated in the ZIB format can be reread by the file command In this ca
6. 3 STARTING THE PROGRAM Figure 20 Mesh and solution on a cut after some refinement steps in example 19 height h T 200 T 300 1 r 50 Figure 21 Axialsymmetric problem in example 20 32 3 STARTING THE PROGRAM 33 0 Neumann boundary dT dn T 50 300 1 r z 1 0 r 00 T 200 r 1 0 r 0 5 Figure 22 Axialsymmetric problem in example 20 u r Or or Oz OT OT rAT rr Tu T apl Te in Z zZ By substitutions r x and z gt y this results in the usual 2d problem with a variable elliptic coefficient k x and the right hand side multiplied by z ee o OT 0 T OT CS x ox Oy Oy At the inner boundary we have the Neumann boundary condition OT On 0 Analogously we can treat the Laplace operator in spherical coordinates 140 7 uuo 1 1 E Par dr e r2cos20 09 e r cosd 295839 3 STARTING THE PROGRAM 34 Multiplying with r and cos0 yields 0 1 0 O o 35 095825 r codA E r2cos0 jr ar Bg cos DG Again by a simple multiplication of the transformed partial differential equation with a special factor depending on the used coordinate system we arrive at the standard problem formulation 1 see appendix A Note that
7. indices of points of boundary edge type classification indices of points of boundary edge type classification END this section is optional Points lt index of point on boundary type classification index of point on boundary type classification END In the section after the keyword Boundary you have to specify the boundary condition for each boundary edge All the points on an edge inherit the boundary property from the edge e g if a point belongs to a Neumann edge it becomes a Neumann point itself If it lies both on a Neumann edge and on a Dirichlet edge it becomes a Dirichlet point The keyword Points after the keyword Boundary opens the optional section for specifying boundary conditions on points This special treatment of points may be necessary in some cases e g if we have only Neumann edges and want to specify one point of Dirichlet type During the refinement process of our algorithm a new point inherits the proper ties of the edge on which it lies New edges or triangles inherit the properties of their fathers As an example we present a unit square divided into two triangles see figure 27 5 GEOMETRY AND MATERIAL PROPERTIES point with index 4 x 0 0 y 1 0 point with index 3 x 1 0 y 1 0 Dirichlet boundary with points 3 and 4 classification 1 1 1 element with points 1 3 and 4 and classification 1 element wi
8. Parameter Name infoErrorEstimator infoLinSystem infoMG infoPrecond infoRefinement infoSolution levelOdirect linSolver localExtend localOnTop localSmooth material maxUrthoGMRes maxRefSteps minRefRatio minRefSteps nPostSmooth nPreSmooth omega Default set in file kaskade init O false 25 O false O false 1 cg 1 true 1 true 1 true defaultMaterial 84 Description info about error estimation print information on every 25th iteration step info about multi grid preconditioning info about preconditioning mesh info levels ranging from 0 to 2 gi ving statistics about refinement steps info after each solution step limit for the direct sparse matrix solver on level 0 i e the solver will be used for a matrix dimension up to 3000 determines the type of iterative solver extend smoothing pattern 0 only new nodes 1 new nodes neighbours local smoothing of top level matrix local multigrid name identifying the material type defaultMaterial defines constant coefh cients on the elements maximal number of vectors to be orthogo nalized in GMRES maximal number of adaptive refinement steps in solve elements marked for refinement at least 5 minimal number of adaptive refinement steps in solve number of post smoothing steps in multigrid number of pre smoothing steps in multigrid relaxation parame
9. ug on Pp o E qN on Uy 1 On o kc Bou qc on Po On where k and a denote material parameters thermal conductivity and transfer coef ficient Q may be a one two or three dimensional region boundary conditions of Dirichlet Neumann or Cauchy type are applied on the corresponding boundary sections I p Py and Pc The source f and the coefficients k q and a may depend on space coordinates An extension to anisotropic materials with tensors k kj is quite straightforward The weak formulation associated with 1 is Find Hp Q satisfying a v G v vv Hi Q 2 where kVuVv quv dxdydz auvds To fodzdydz qcvds J oves o TN There is the possibility in KASKADE to compute one additional boundary integral on the right hand side of the equation 2 e g due to a special source in the region It has the form G v farvas 3 Pr Here I denotes an interior layer e g the face between two different materials and qr is a function defined on this layer A PROBLEM CLASSES 81 Transient Heat Conduction This class comprises linear parabolic problems of the form e De YkVu qu f in Q 4 The boundary conditions are like in 1 but may be time dependent like the source term f Additionally we have to define initial values for t 0 u z 0 uo x in Q uo x has to be supplied by the user or can be calculated via the solution of a static heat conduction problem Nonlinear
10. VILLAIN NER DIENEN TEST LIU ANSA X 1 3 STARTING THE PROGRAM Figure 12 Initial mesh and mesh after some refinement steps in example 11 quasi 3d plot after some adaptive refinement steps and in Figure 14 you can see the mesh and solution after seven refinement steps 13 Static problem using the files layer 2d cmd layer 2d geo layer 2d mat Here we want to solve the same equation as in Problem 12 But additionally we consider a term of type 3 in the weak formulation of the problem see appendix A For example problems of this kind appear in electrodynamics in 0 25 0 75 x 0 25 0 75 1 1 elsewhere in the unit square 100 U zx Uyy on the boundary 10 0 On the layer between the two materials which must be marked by type 1 and 2 in the mat file we integrate the term k ka x y where the material 100 l and ka integration must be initiated by the command innerBoundary coefficients are k This 3 STARTING THE PROGRAM 26 Figure 13 Initial mesh and solution plot after some refinement steps in example 12 Min 10 Max 10 05
11. targets in kaskade make to obtain separate versions for different space dimensions and one that comprises all of them gt make gt make gt make gt make f kaskade make k1 1D code k1 f kaskade make k2 2D code k2 f kaskade make k3 3D code k3 f kaskade make k6 code included for all space dimensions Or shorter if you have made a copy makefile of kaskade make make make gt make make k1 1D code k1 k2 gt 2D code k2 k3 gt 3D code k3 k6 code included for all space dimensions Additionally the desired space dimension has to be set in the file dimension h k6 is the default target 2 INSTALLATION 10 The default compiler is GNU s g version 2 7 2 If you prefer GNU s g version 2 6 3 you have to change the type of complex in the file general h Addi tionally we use the C Set for AIX 6000 version 2 1 on IBM RS 6000 work stations and DEC C V5 0 3 for DEC OSF 1 2 3 More about the Makefile The makefile helping to install KASKADE is called kaskade make It includes the configuration of system resources necessary to compile and link the program Here you can select your compilers C Fortran loaders and the pathnames for the libraries you need on your machine e g X11 Fortran We mentioned examples for some platforms Maybe they are not up to date In the makefile you can select your compile and loader options e g FORFLAGS or OPTFLAG Finally all th
12. 3d step mat with the suitable material and boundary condition e obstacle Here we present a special stationary obstacle problem In this class of problems we have to solve variational inequalities with constraints In chapter 3 nonlinear problem 1 we give a more detailed description of this 2D problem e stefan The stefan problem describes the transition between two phases of a material e g a melting iceberg 12 FREQUENTLY ASKED QUESTIONS 79 Problem This time dependent problem is formulated in form of variational inequali ties See chapter 3 nonlinear problem 2 for a detailed description of this 2D problem Numerics For this class of nonlinear problems we provide a special error estimator and two suitable multigrid methods as linear solver e porousMedia This example shows the time dependent diffusion in porous media In chapter 3 nonlinear problem 3 you find more informations on this 2D problem 12 Frequently Asked Questions Here we present a list of frequently asked questions e Which file is responsible of what e How to add a convective term e How can we handle periodic boundary conditions Shortly you can read here the answers too A PROBLEM CLASSES 80 A Problem Classes Here we give a description of the installed problem classes Static Heat Conduction and Related Problems This problem class involves elliptic partial differential equations of the type VkVutqu f in Q u
13. 5 steps too But less than 500 nodes This example shows the efficiency of adaptive methodes e boundaryLayer The solution of this elliptic 2D problem has a boundary layer along the sides x and y 1 Problem equation Au 100u 0 domain unit square cosh 10x cosh 10y boundary condition Dirichlet defined by 2 cosh 10 Numerics In this example it is not reasonable to vary the material data but the user can test the different numerical parameters like in the default example e corner The solution of this Laplace equation on a reentrant corner 2D problem has a singularity in the origin Problem equation Au 0 domain part of 2D hexagon with reentrant corner boundary condition On parts of the boundary we have homogeneous Neu mann conditions which are defined in the file corner mat On the rest of the boundary there are non constant Dirichlet conditions u r y r with r yr Fy They are defined in the source file dirichletA cc In near future we will present a function in the ZIBGui in order to modify interactively parts of the code and to recompile the executable automatically Numerics In this example the user can test the different numerical parameters like in the default example 11 GRAPHICAL USER INTERFACE TT e staticPeak The solution of this stationary Poisson problem has a peak at a point Problem equation Au f the source f has a peak of the form
14. A dad 11 Graphical User Interface 11 1 How to install the ZIBGui environment 1 152 Examples A E E A 12 Frequently Asked Questions A Problem Classes B Command Listing References 67 67 67 70 71 71 72 79 80 83 87 LIST OF TABLES 5 List of Tables PON ant parameters siTe RT EE e e Rv AS Be ER 13 problem peace La edi RO RULES E EUER EMEN 37 maternal types sasia da a dea eee Do RE Doo ue dee dde A Id dE 54 Dirichlet boundary conditions wav ar ove Re DR eI Ere Res 55 rara e cas eu doa erea de ek ie T an 58 parameters for linear solvers soy uoc oa Pe ee AL IR 59 preconditioners for linear problems llle 60 preconditioners for nonlinear problems 60 Error OSLO S gt urere ie elite od n e edu e d eue I Hg Roof uf xD eng 61 refinement strategies ser CU ERE RIP EI E EIEI TE XR o 61 info and print commands ux mcs qw qub Svid erm 68 plot commeands ae x ee eS Ge a UR E E E dde Ee den 69 commands for time informations a nn 70 LIST OF FIGURES 6 List of Figures 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 main iteration loop ira AEG POE RAUS REQUE RA Nub 12 example peak ld zu Lo af a e o BUR Barren 16 example pesada qan are aun deca qo C e Sod uf d th dd 17 example peak 3d ses ese aa aa ra a BOR eS 19 examples unit 2d A ouo a A E UE SX 20 example SI 20 dB Lo ADE Y O O RA AA d 20
15. FOR PROGRAMMERS 69 Parameter Name graphic plotBoundary plotElements plotIsoLines plotSolution plot3d plotKeep plotLevels plotPointNodes plotSize plotTimeStep plotTriangleNodes postScript scaleX scaleY scaleZ Default set in file kaskade init false true false true true false 0 false 10 0 false 0 4 1 0 false 0 false Description switch graphic support on off plot no boundary when solution is plotted plot no elements when the solution is plotted plot isolines of solution plot solution data after each space step plot 2D results as a surface plot quasi 3d plot keep all plots in separate windows desired number of iso lines plot number of nodes at points size of plot window plot the solution after each time step plot number of nodes at triangle generate postscript pictures of mesh and solution scaling x coordinate for quasi 3d plots scaling y coordinate for quasi 3d plots scaling z coordinate for quasi 3d plots Table 12 Plot commands 10 TECHNICAL DETAILS FOR PROGRAMMERS 70 Parameter Name Default Description set in file kaskade init accTime false print all accumulated cpu times time false print cpu time for all essential modules timeAssembler false print cpu time for matrix assembling timeErrorEstimator false print cpu time for error estimation timeLinSystem false print cpu time for system solution timeR
16. Gauss Seidel preconditioner Table 8 Preconditioners for nonlinear problems 7 ERROR ESTIMATION AND MESH REFINEMENT 61 The local information gives the error on each element on the actual mesh We use it to determine a set of elements to be marked for refinement An estimator is selected by the parameter errorEstimator x where x represents one of the keywords in the left column of table 9 Keyword Estimator DLY estimator due to DLY89 see also BEK92 modDLY modification of DLY see BEK93 in case of constant coef ficients and constant right hand side identical with the stan dard DLY estimator 3D estimator due to BEK93 estimator for nonlinear problems due to Kor95 estimator for nonlinear problems including one multigrid cycle for hierarchical defect correction due to Kor95 Table 9 Keywords for Error Estimators We offer some refinement strategies The parameter refStrategy invokes one of these where represents one of the keywords in the left column of table 10 Keyword Strategy maxValue refinement with respect to the maximum of local errors extrapolation refinement with respect to extrapolated errors uniform uniform refinement Table 10 Keywords for Refinement Strategies The refinement strategy computes a threshold and marks those elements where the estimated local error is beyond this threshold Use the command minRefRatio 8 HOW TO DEFINE A NEW PROBLEM 62
17. Problems Here the problem classes have an additional function describing the nonlinearity It defines the nonlinear characteristics of the problems e g obstacle functions or a nonlinear heat capacitance Obstacle Problems Two obstacle functions Pup 1 Plow x can be added to a differential equation of type 1 u r Pup x 5 u x gt Plow x Such constraints may occur in various types of applications The corresponding va riational inequality yields a nonlinear system which is solved by a single level Gauss Seidel or a monotone multigrid method Kor95 Stefan Problems and Related Nonlinear Equations These problems may be regarded as degenerate parabolic initial value problems of the form a E00 VkVu f in Q MEN Mole in Q 6 H x 0 Holz in Q A PROBLEM CLASSES 82 combined with boundary conditions like in 1 H is the heat content or a generalized enthalpy In KASKADE H is the discontinuous derivative of a continuous piecewise quadratic function of u Thus H must be specified as a piecewise linear function of u a discontinuity in A defines a change of phase with a certain amount of latent heat Other nonlinear problems may be formulated in a similar way e g the porous media equation see Kor95 Of course we may have additional obstacle functions of the form 5 B COMMAND LISTING 83 B Command Listing On the following pages we present an alphabetical list of the command parameters Para
18. Qaagaganaiaaaaa a a PRR PRP RP RP RP RB SR GS EP END We present as example for a little bit more complicated 3D geometry a cube with a hole inside The domain is Q x y z 0 lt z y z lt 1 0 4 ayy 2 0 25 lt x y z lt 0 75 and a coarse tetrahedral mesh is given in holes 3d geo see also figure 30 showing a cut through z 0 5 The faces outside are marked as Dirichlet type the inner ones as Neumann type In layer 3d geo you find the geometry similar to the last example but the hole is filled with another material The triangulation of the interior is generated by an additional point in 0 5 0 5 0 5 and 48 tetrahedra Figure 31 shows the configuration On the boundary we defined Dirichlet conditions The generation of a sphere in 3D is analogous to the 2D case We also use the character S to mark a boundary face for curvilinear adaption during refinement The ball must be centered in the origin For example see the files circle3d cmd circle3d geo circle3d mat We can approximate a cylindrical region too On the initial grid we have to mark by a C those boundary faces which shall change during refinement to the curvili near faces of a cylinder For example we start from the following cube compare cylindrical 3d geo 51 5 GEOMETRY AND MATERIAL PROPERTIES Figure 30 Cut through a cube with a hole Figure 31 Cut through a cube consisting of two materials 5 GEOMETRY AND MATERIAL PROPERTIES
19. STARTING THE PROGRAM 30 19 COSAS DAVA VA Figure 18 Surface of the initial triangulation of a skull 18 position of the antennas and frequency of the electric field have to be chosen properly to achieve a local heating of the tumor compare HYP89 This requires the efficient and robust solution of the stationary or transient BHT bio heat transfer equation Pen48 EU IET O with coefficients k and q depending of the tissue Qe describing the power of the electrical field generated by radiowaves in the tissue We refer to Seeba See90 for a more detailed explanation and the numerical data of the following example We consider the heating of the tumor by the electric field resulting from three implanted antennas A suitable description of the computational domain Q is obtained by computer tomography Figure 18 shows the surface of and the initial triangulation To from Seba See90 Note that To is chosen such that different tissues as bone brain and tumor belong to different tetrahedra The adaptive refinement concentrates on the neighborhood of the antennas and the tumor in the interior of the skull where we have high gradients of the tem perature This is illustrated by the figure 19 showing the final triangulation 7 and the contour lines of the temperature on some clipping plane The different tissues are indicated by black bone white brain and grey tumor inside the brain Static problem
20. coefficient k the mass term q and the source f Real UserStaticMaterial E int type Vector lt Real gt x int i int j Vector lt Real gt amp xx x switch spaceDim case 1 return Eid xx 1 case 2 return E2d xx 1 xx 2 case 3 return E3d xx 1 xx 2 xx 3 8 HOW TO DEFINE A NEW PROBLEM 65 default abort return 0 Real UserStaticMaterial S int type Vector lt Real gt x Real time Vector lt Real gt amp xx x switch spaceDim case 1 return Sid xx 1 case 2 return S2d xx 1 xx 2 case 3 return S3d xx 1 xx 2 xx 3 default abortO return O Real UserStaticMaterial M int type Vector lt Real gt x int i int j 1 Vector lt Real gt amp xx x switch spaceDim t case 1 return Mid xx 1 case 2 return M2d xx 1 xx 21 case 3 return M3d xx 1 xx 2 xx 3 default abort return 0 The coefficient k is defined in the function Real UserStaticMaterial E which calls one of the functions E1d E2d or E3d depending on the space dimension The coefficient q is computed in the function M and the source term f in the function S Available are functions for variable Neumann Cauchy or Dirichlet boundary condi tions To specify a Cauchy condition we define two functions the coefficient a of u is provided by the function Real UserStaticMaterial Cauchy and the right hand side by the function Real UserStaticMaterial Neumann In order to disti
21. ex e It changes parameters and generates a new parameter file e It starts the KASKADE executable The interface is based on the tool Tel Tk due to J Ousterhout Ouster94 for develo ping graphical interfaces It is independent of the numerical code KASKADE Each ex file is recognized as an example and the parameters in it can be interpreted and changed However the use of the parameters is more restrictive than for the parser in the KASKADE program Only those parameters are understood correctly which con cur fully with the notation of the table in appendix B Abbreviations are not allowed capital letters are significant You find a detailed description in German in the manuals NPR96 and NPRW96 including examples for customizing the user interface for your changes in the numerical application e g you have a new parameter 11 1 How to install the ZIBGui environment First install the public domain software Tcl Tk version Tcl 7 4 and Tk4 0 on your system Then execute the following steps in order to translate the ZIBGui sources which you find in subdirectory zgui e Goto subdirectory zgui e Set in the Makefile the path to the directories of your Tcl Tk libraries and include files Additionally adapt the name of your operating system for instance OS SUNOS This name will be used to store object files and libraries in the subdirectories obj OSname and 1ib OSname e Type make to compile and link the files The name
22. f c e p z a in the 1D case and analogously in higher dimensions It is defined in the source file materialsA cc as a function of the coordinates domain unit interval in 1D unit square in 2D cube in 3D boundary condition homogeneous Dirichlet Numerics In this example it is not reasonable to vary the material data or boundary conditions but the user can test the different numerical parameters like in the default example It is very interesting to study the differences between adaptive and uniform refinement Variation of parameters change the region by the geo file button e g select peak 2d geo peak 3d geo change the numerical methods similarly to the tests in the default example e transientPeak The solution of this parabolic problem has a peak at a point of the domain Problem equation u Au f the source f has a time dependent peak of the form foc e p a t in the 1D case and analogously in higher dimensions It is defined in the source file materialstr cc as function of the coordinates and the time domain unit interval in 1D unit square in 2D unit cube in 3D boundary condition homogeneous Dirichlet conditions Numerics In this example it is not reasonable to vary the material data or boundary conditions but the user can test the different numerical parameters like in the default example Variation of parameters Change the region by the geo file button e
23. g select transpeak 2d geo unit 3d geo We recommend to request less accuracy factor globalPrecision in higher space dimensions 11 GRAPHICAL USER INTERFACE 78 Change the numerical methods like in the default example e transient Diffusion The solution of this parabolic problem shows a diffusion process starting from a sharp initial profile Problem equation 1 0 for t lt 10 0 POSS l 0 0 fort 10 0 The coefficients are described in the file u 1d step mat u 2d step mat resp u 3d step mat depending on the space dimension domain unit interval in 1D unit square in 2D cube in 3D boundary condition u 0 on the left side u 1 opposite to it and homoge neous Neumann condition on the other parts of the boundary initial value is given as function in the file dirichlettr cc and is available under the name Step Numerics This example demonstrates the parabolic solver which is adaptive in space and time In each time step the computation of the spatial solution starts from the initial grid In this example it is not reasonable to vary the material data or boundary conditions but the user can test the different numerical parameters like in the default example Variation of parameters change the region by the geo file button e g select u 1d step geo u 2d step geo or u 3d step geo If the user changes the geometry in this example he has also to select the mat file u 1d step mat u 2d step mat or u
24. of the generated executable is zgui e Set in the file zgui k6 tclthe path to the directory with the KASKADE ap plication k6 In the standard case the path is the parent directory of the zgui 11 GRAPHICAL USER INTERFACE 72 hd ZIB GUI Kaskade 3 0 Example Settings Information Be ee k6 in home newton bzferdmaA 0 For Heip point on object press shit and the fem mouse button Il II ld Figure 34 Starting menu of the graphical user interface directory Further you can modify the path to the library with the tcl files of zgui By default it is zgui library e If you have not already installed the KASKADE executible k6 please do it now e Simply start the graphical interface by typing gt zgui 11 2 Examples When ZIBGui is started it opens a window similar to that shown in figure 34 Please note If the following pictures of some ZIBGui windows are a little bit blurred try to print these pages on a 600dpi printer The menus of the graphical user interface ZIBGui reflect the structure of the nume rical software Before starting the KASKADE code by clicking on the Start button in the application control board you may change your mathematical problem type material and boundary conditions the components of the multilevel solver linear solver preconditioner convergence criteria error estimator and select some graphical and protocol parameters in the submenus of the Settings menu in the ro
25. provide appropriate functions in the source code The material class is specified by the parameter material e Dirichlet boundary conditions The command DirichletBCs selects the con stant or variable Dirichlet boundary conditions 3 STARTING THE PROGRAM Parameter Name dirichletBCs file globalPrecision material linSolver pause problem SpaceDim Default set in file kaskade init constDirichlet unit 2d geo defaultMaterial cg 1 true staticHeat 13 Description identifies the type of Dirichlet boundary condition name of geometry file the extension geo need not be specified The material file is expected to have the same name but the extension mat desired relative precision maximum dis cretization error with respect to global energy if absPrecision 0 otherwise de sired absolute precision name identifying the material type defaultMaterial defines constant coeffi cients on the elements determines the type of iterative solver stop when function Pause is called in the code and wait for input lt CarriageReturn gt continue until next Pause is encountered c continue and disable the function Pause p gt generate a picture in postscript for mat of the actual mesh and the approxi mate solution q gt quit program specifies the problem type to be allocated and solved space dimension Table 1 Some important parameters 3 STAR
26. strategy uses the local error information to mark elements for refinement it can be chosen by the command refStrategy 4 We repeat the steps 2 and 3 until a convergence or break condition stops this cycle Normally the process ends if the desired precision defined by command globalPrecision or a maximal number of refinement steps defined by com mand maxRefSteps is reached The command minRefSteps prescribes a minimal number of refinement steps By default we look at the estimated error relative to the energy norm of the approxi mate solution Only if the parameter absPrecision is set to true we compute until the absolute error is below the threshold globalPrecision This break condition is reasonable if we want to approximate a solution with a range near zero Figure 1 illustrates the main iteration loop of this adaptive process The self adaptive grid refinement minimizes the number of nodes to achieve the re quested precision Additionally the refinement history allows an efficient multilevel preconditioning of the linear systems For linear parabolic partial differential equations we use the adaptive approach of F Bornemann see Bor90 which is able to handle complicated geometries and in consistent initial data Here the time is discretized first allowing an adaptive stepsize 3 STARTING THE PROGRAM 15 control The errors of the arising spatial elliptic subproblems are controlled indepen dently Thus advanced adap
27. term recurrence relax relaxation routine Jacobi SSOR etc for technical reasons the type is determined by the choice of a preconditioner applied to the Richardson iteration gmRes GMRES of Saad and Schultz SS86 stdBiCG bi conjugate gradients standard version biCGStab BiCGStab of van der Vorst vdV92 stabilized version cgs the CGS algorithm of Sonneveld Son89 1sqCG conjugate gradients for the normal equations CGNR nonLinRelax nonlinear relaxation routine to be used in conjunction with a nonlinear preconditioner Table 5 Keywords for iterative solvers In table 6 we list the commands which are relevant for the iterative solution procedures 6 2 Preconditioners for Linear Problems The possible choices are summarized in table 7 6 3 Preconditioners for Nonlinear Problems These types are summarized in table 8 7 Error Estimation and Mesh Refinement The error estimator provides estimates for the global and the local discretization errors The global error is used to stop the algorithm when the requested precision is reached 7 ERROR ESTIMATION AND MESH REFINEMENT 59 Parameter Name directSolverLimit extPrecFactor infoLinSystem levelOdirect linSolver maxUrthoGMRes omega preConditioner timeLinSystem Default set in file kaskade init MGsgs O false Description limit for the direct sparse matrix solver on levels gt 0 see also levelOdirect ex
28. to request a minimal number of elements to be refined The default is 0 05 which means that at least 5 96 of all elements will be refined If no estimator or no refinement strategy is selected all elements will be refined 8 How to Define a New Problem In the preceding chapters we already learned how to solve a problem class of type 1 or 4 see appendix A with piecewise constant coefficients and piecewise constant boundary conditions But a lot of applications also treat material properties varying in space i e k k x e or f f x We summarize the steps for defining a new problem e Define the triangulation of the domain in a geo file e Define boundary conditions and material properties the coefficients of the equa tion in a mat file if piecewise constant or in the file userStatic cc resp userTrans cc if variable e Specify your problem in a cmd file 8 1 Everything is variable You can handle equations with variable coefficients and boundary conditions very ea sily by defining the coefficient functions in the file userStatic cc for static problems 1 or in userTransient cc for transient problems 4 In these files we prepared sample functions which can be overwritten by the user They are named EC for the elliptic term k MO for the mass term q or P for the coefficient of the time derivate in 4 and SO for the source term f Each of these functions calls a subfunction depending on the space dimensio
29. using the files cylindrical 3d cmd cylindrical 3d geo cylindrical 3d mat Poisson equation with constant coefficients in a cylindrical region Z Z r 9 z with radius r 0 y2 0 and ze 1 1 Uss Uyy tty 4 0 z y e7 in the cylinder Z 3 STARTING THE PROGRAM 31 SE SELLER DATE Figure 19 Clipping of the mesh on refinement level 7 and temperature on this plane in example 18 20 u y e on the boundary We start on a coarse triangulation of the cube 1 1 x 1 1 x 1 1 During refinement we approximate the shape of Z by placing new point coordinates on the surface of the cylinder Figure 20 shows the first three steps of uniform refinement and figure 20 the mesh of the cylinder and isolines of the solution on a vertical cut through it Static problem using the files cylindrical 2d cmd cylindrical 2d geo cylindrical 2d mat Poisson equation in a cylindrical domain Z formulated in cylindrical coordinates r and z The cylinder Z Z r z has radius r 1 and height z between and 1 AT e with 19 OT 18T OF ra oU c amc oe UB The boundary conditions are shown in figure 21 On the top and on the cur vilinear faces we have a temperature of 50 C On the bottom there is an area r lt 0 5 with temperature of 200 C and 50 300 1 r C elsewhere This problem has axial symmetry and can be treated as a 2d problem as des cribed in figure 22
30. D objects writeCpuTime false write cpu time information to info file writeGrapeFormat false set format for geometrical output it cor responds to GRAPE visualization tools only for 3D objects writelterations false write iteration numbers to info file writeZIBFormat false set format for geometrical output it cor responds to the standard ZIB format writeSolution false alow output on a file the format of output has to be defined by one of the commands writeZIBFormat writeGRA PEFormat writeAVSFormat writeTimeStep false write the solution after each time step REFERENCES 87 References Bor90 F A Bornemann An Adaptive Multilevel Approach to Parabolic Equations I II IMPACT Comput Sci Engrg 2 1990 and IMPACT Comput Sci Engrg 3 1991 BEK92 F Bornemann B Erdmann R Kornhuber Adaptive Multilevel Methods in Three Space Dimensions Int J Numer Meths in Eng 36 1993 BEK93 F Bornemann B Erdmann R Kornhuber A Posteriori Error Estimates for Elliptic Problems in Two and Three Dimensions Preprint SC 93 29 Konrad Zuse Zentrum f r Informationstechnik Berlin 1993 BER95 R Beck B Erdmann and R Roitzsch KASKADE 3 0 An Object Oriented Adaptive Finite Element Codes Technical Report TR 95 4 Konrad Zuse Zentrum f r Informationstechnik Berlin 1995 BPX90 J H Bramble J E Pasciak and J Xhu Parallel Multilevel Preconditioners Math Comp Vol 55 1990 Cr88 J Cran
31. E T ir E ss Baling u Sea E Je El 13 ae i 4 2j Figure 23 Final mesh and isolines of solution after 8 refinement steps in example 1 3 3 3 Nonlinear Problems 1 Nonlinear Problem using the files obstacle cmd obstacle geo obstacle mat Here we treat a stationary 2D problem of type 5 in appendix A with cons traints u x lt dist x O0 see figure 23 This obstacle problem is modelling a form of elasto plastic torsion of a cylindrical bar with cross section Q The region along the diagonals of Q is inactive i e not on the obstacle Figure 23 shows the final mesh and the contour lines of the solution after 8 adaptive refinement steps Nonlinear Problem using the files stefan cmd stefan geo stefan mat The two dimensional Stefan problem 6 in appendix A describes the transition between two phases of a material e g a melting iceberg Two phase Stefan problems of this form are arising in the simulation of the heat conduction in a substance undergoing a change of pase at the nominal phase change temperature 01 In this case the function U is a generalized temperature resulting from a standard Kirchhoff transformation U Ko for 0 lt 6 and U k10 for O gt 0 of the physical temperature 0 3 STARTING THE PROGRAM 36
32. TING THE PROGRAM 14 e Space dimension It is set by parameter spaceDim e A coarse triangulation for the region 2 in a geo file This file is specified by the parameter file 2 On the actual triangulation of the region 2 the weak formulation of the partial differential equation is discretized by linear finite elements The resulting linear system is usually solved by a direct sparse matrix Cholesky me thod This direct sparse solver is only efficient as long as the dimension of the matrix is not too large We use the command level0direct to select an iterative solver instead of Cholesky s method i e the direct solver will be used on the first grid level 0 for a matrix dimension up to The switch between direct and iterative sol ver on higher refinement levels is defined by the command directSolverLimit where stands for the maximal matrix dimension for the direct solver on levels gt 0 The iterative solver may include a preconditioner to reduce the number of iterations for the new approximate solution The iterative solver must be specified with the command linSolver the preconditioner with the command precond 3 The global discretization error of the approximate solution is estimated If it is below the required precision the process stops Otherwise we use local error indicators to refine the mesh where the errors are beyond a threshold The estimator is selected by the command errorEstimator A refinement
33. We are extending this guide continuously The latest version is available by network CONTENTS 3 Contents 1 Introduction 8 2 Installation 8 2 1 Obtaining the Source Code o e 8 2 2 Compiling and Linking the Code v a ERE eR 9 2 3 More about the Makefile 10 3 Starting the Program 11 3 1 Parameters and Parameter Files 11 32 Algorithmic Kernel oy 2 4 a anor Auris EEE ER OS AA DS 12 2 07 ERA IOS uat se be rite a e atur ds or rate ee 15 3 3 1 Static Problems eh 15 3 8 2 Transient Problems 2 2 2 om on 34 3 3 3 Nonlinear Problems 35 4 Problem Classes 37 5 Geometry and Material Properties 39 51 Definition of the Geometry oaoa ar aa y mud enn 39 5 2 Ouput of Geometrical Data rs 53 5 3 Material Coefficients Boundary Conditions and Initial Values 54 6 System Solution and Preconditioning 57 6 1 The Iterative Solvers 2 oo oo s 57 6 2 Preconditioners for Linear Problems 58 6 3 Preconditioners for Nonlinear Problems 58 7 Error Estimation and Mesh Refinement 58 8 How to Define a New Problem 62 8 1 Everything is variable anes dat due du ide res 62 8 2 Something is variable something is constant 67 CONTENTS 9 Obtaining Run Time Information 10 Technical Details for Programmers 10 1 Vector and Matrix Classes o aaa e note RES 10 2 CAOS A RL de d ue IA O
34. cation indices of points of boundary triangle type classification 5 GEOMETRY AND MATERIAL PROPERTIES 49 END optional section Points lt index of point on boundary lt type gt lt classification gt lt index of point on boundary lt type gt lt classification gt END In the section after the keyword Boundary you have to specify the boundary condition for each triangle All the points on an triangle inherit the boundary property of the triangle If a point belongs to a Neumann and to a Dirichlet triangle it becomes a Dirichlet point Edges are treated analogously The keyword Points after the keyword Boundary opens the section for specify ing boundary conditions on points This special treatment of points may be necessary in some cases e g if we have only Neumann triangles and want to specify one point of Dirichlet type This data section can be omitted During the refinement process of our algorithm a new point inherits the pro perties of the edge on which it lies New edges triangles or tetrahedra inherit the properties of their fathers As an example we present a unit cube divided into six tetrahedra unit cube points maxIndex 10 1 001 2 011 3 111 4 101 5 000 6 010 7 110 8 100 END elements 1236 1 1356 1 3567 1 1345 1 5 GEOMETRY AND MATERIAL PROPERTIES 50 w A e q Cc N N END boundary faces N w oOowrPAPRPRRP OWN RP 5 015000 0 0 IN NOONAN Y 0 0 O
35. cess by placing new boundary points onto the circle In figure 11 the meshes after 4 9 and 12 refinement steps are presented Static problem using the files unit 3d cmd unit 3d geo unit 3d mat Poisson equation with constant coefficients on the unit cube Ure Uyy uu 1 in 0 1 x 0 1 x 0 1 u 0 on the boundary 3 STARTING THE PROGRAM 24 Figure 11 Initial mesh and meshes after 4 9 and 12 refinement steps in example 9 11 12 Static problem using the files 1 3d cmd 1 3d geo 1 3d mat The same equation as in Problem 10 but on a region with reentrant corner as described in 1 3d geo see figure 12 The singularity of the solution in the re entrant corner is well resolved by the error estimator Static problem using the files jump cmd jump geo jump mat The same equation and region as in Problem 5 but with a jump in the coeffi cients 100 uzs uy 1 in 0 25 0 75 x 0 25 0 75 Use uy 1 elsewhere in the unit square u 10 0 on the boundary In this example it is recommended to work with the option plot3d 1 for gra phical illustration Figure 13 shows the initial coarse mesh and the solution as 25 N WSN NISSAN LODOS POSE SIS EROS BSA IN AS CH ETT PL PEA K BOSON TARTA LF DESSA NERO ANGE ARSON NIK ORO LLL LG 2 gt A EN ESS A SSISSISSSSISIISSISA ISI 77 AZ T ELCH DEE NND FIT LIAN J LED LENA DS LIRA RDA ADA RR
36. d the name transHeat to linear scalar parabolic problem because we cannot solve general elliptic or parabolic equations e Use the geo file button to select the geo file which describes the triangulation 1D 2D or 3D of the region e The mat file button offers the possibility to select a mat file describing pie cewise constant material coefficients or constant boundary conditions If you have variable coefficients and boundary conditions the specification of a mat file is not necessary e The global precision parameter specifies the error tolerance for the solution pro cess The adaptive process will be continued until the estimated error lies under this threshold The value for maximal level provides a possibility to stop the multilevel solver after a certain number of steps independently of the achieved precision e Now you can choose the setup for the multilevel solver First you may select the iterative method to compute the solution of linear systems The system may be 11 GRAPHICAL USER INTERFACE 74 preconditioned by one of the offered methods Depending on the problem type you have some suitable error estimators which can provide an estimation of the global error This is important to control the multilevel process In addition the estimator computes local errors which can be used together with a refinement strategy to generate an adaptive grid e The coupling of requested precision in the iterative solver of the lin
37. e dimensions shows the time dependent diffusion in porous media In our example we select m 2 and the initial condition 6 9 21 23 0 4 r 1 0 55in 149 using r x 22 and q atan x2 x1 We solve on Q 0 1 x 0 1 and T 0 0 05 prescribing homogeneous Neumann conditions at x 0 22 0 and homogeneous Dirichlet conditions at x 1 2 1 Again we use only a fixed stepsize for time discretization 4 Problem Classes The program handles several problem types Anyone of these is selected by the com mand problem x x x where represents one of the keywords in the left column of table 2 The problem type is used internally to describe the structure of the differential equa tion and to select the suitable solver routines Keyword problems staticHeat linear elliptic problem e g static heat conduction transientHeat linear parabolic problem e g transient heat conduction obstacle obstacle problem stefan Stefan problem porousMedia porous media equation Table 2 Keywords for problem types See appendix A for a desription of basic problem classes 38 4 PROBLEM CLASSES eA A b E SCA E H z E i
38. e file dependencies are listed Thus we are sure that changes in one of the source files initiate a new compilation user is recommended to update this list after he added new dependencies e g include files That can be done by make f kaskade make dependencies generating a new makefile makefile with all discovered dependencies 3 STARTING THE PROGRAM 11 3 Starting the Program 3 1 Parameters and Parameter Files Please see section 2 how to obtain the source files and how to compile and link the code The name of the program depends on the chosen program version 1D ki 2D K2 3D K3 all three space dimensions K6 The program is started by specifying a parameter file which defines a problem and sets relevant parameters e g gt k6 cmd unit 1d cmd Here unit 1d cmd is the command file for the selected problem 1D static heat con duction Parameters may also be set on the command line e g gt k6 cmd unit 1d cmd graphics 1 globalPrecision 1 0e 2 Commands have the following syntax 1 command example pause 2 command number example spaceDim 3 3 command keyword example problem staticHeatConduction The case 1 is only possible for commands of boolean type and is equivalent to command 1 where 1 means true The value is set to false by command 0 Some general remarks on commands and command files e All commands are read at startup Thus the program execution cannot be altered during run t
39. ear system and the estimated discretization error is tuned by the convergence test strategy Both the residual and the cascadic iteration method are used for elliptic problems In general the former one is recommended see BER95 e For parabolic equations you find entries in the problem definition window for the start time and final time e The graphics setting menu provides different graphical presentations of the solution process e g plots of the solution and the grid The graphic output is X11 based e The amount of information of the program during runtime is tuned by the protocol settings For further evaluations of the protocol information we in troduce the Special Statistic button The tool which uses this additional output is not yet part of the KASKADE package For 3 D examples there is no interactive graphic support integrated in this version of the program Some examples e default This is the example which is automatically selected when the program starts Problem Static heat conduction in two space dimensions constant coefficients k 1 and constant right hand side s 1 and homogeneous boundary conditions Problem settings Type staticHeat Material constant Dirichlet constant geo file selection file unit 2d geo unit square mat file selection file unit 2d mat Numerics Conjugate gradient method with a multigrid preconditioner sym Gauss Seidel smoother error est
40. efinement false print cpu time for mesh refinement timeTransport false print cpu time for data transport between different grids timeUpdate false print cpu time for update operation in node management Table 13 Time information such control is expensive in time Therfore we provided a switch to turn on or off the array bound check by setting the flag CheckBoundsFlag in the files vector h and matrix hto True or False default A change of the flag is effective after recompiling the code 10 2 Batchjobs When the program fulfills a break condition e g the requested precision is reached it stops with the prompt REGULAR PROGRAM TERMINATION CR It terminates not before the user types the RETURN key Sometimes this pause is not very practical e g if you want to run a sequence of batchjobs controlled by a shell script In such applications you can initiate an automatic termination of each job by setting the parameter batchJob 1 in the corresponding cmd file 11 GRAPHICAL USER INTERFACE 71 11 Graphical User Interface Often it is difficult to learn the handling of a program and to understand the meaning of the parameters To make these first steps easier we support the user by a graphi cal user interface ZIBGui The ZIBGui is a manager only for the most important commands and parameters e It reads the usual cmd files and displays the parameters For practical reasons we changed the extension of those files to
41. egion comments Points maxIndex mazIndex lt index of point gt lt a coordinate of point gt lt index of point gt lt a coordinate of point end Elements indices of points of element classification material id 7 indices of points of element classification material id 7 END 5 GEOMETRY AND MATERIAL PROPERTIES 40 element with points 3 and 4 element with points 2 and 4 ne and classification 2 and classification 1 e e e point with index 2 x 0 0 point with index 4 x 0 5 point with index 3 x 1 0 Figure 26 Partition of unit interval Boundary index of point on boundary boundary conditon type classification index of point on boundary boundary conditon type classification END The classification parameter has to be greater than 0 and less than 128 The same is valid in the 2d and 3d case As an example we present a unit interval divided into two intervals see figure 26 1D unit interval Points maxIndex 4 2 0 0 4 0 5 3 1 0 end Elements 24 1 34 2 end Boundary 2 D 3D1 end e 2D region 5 GEOMETRY AND MATERIAL PROPERTIES A comments Points maxIndex mazIndex lt index of point gt lt ax y coordinates of point gt lt index of point gt lt xzx y coordinates of point end Elements indices of points of element classification indices of points of element classification END Boundary
42. er compareSolution 1 we can compute STARTING THE PROGRAM 16 Figure 2 Final mesh and solution in example 1 the maximum error of the approximate solution in the nodes In the chapter 8 we explain how you can compare the approximate solution of your problem with a given function e g the true solution 2 Static problem using the files peak 2d cmd peak 2d geo peak 2d mat This problem is the 2D analogue to problem 1 usr uy f in 0 1 x 0 1 u 0 on the boundary Here the source function f exp 100 0 x x 0 5 y 0 5 has a peak in the center of the unit square The inital grid consists of only four triangles In figure 3 you see the grid after 6 refinement steps and the solution in a quasi 3d plot 3 Static problem using the files peak 3d cmd peak 3d geo peak 3d mat 3 STARTING THE PROGRAM 17 Figure 3 Mesh and quasi 3d plot of the solution after 6 refinement steps in example 2 3 STARTING THE PROGRAM 18 This problem is the 3D analogue to problem 1 uas Uyy tios y in 0 1 x 0 1 x 0 1 u 0 on the boundary Here the source function f exp 100 0 x x 0 5 y 0 5 2 0 5 has a peak in the center of the unit cube The algorithm starts with a partition of the cube into 6 tetrahedra In
43. example slit 2d 45 sd a a ar ees 21 example Side rta m esee e a ek ted 22 example SIUS AC rane IA root boe d Kr ne amp ee vf ge 22 example slit solution dar a re E Po TE QUEUE ETUR 23 example slit 2d a wx Ee ER Gite Pe E US EX US EM EE 24 example corner 3d amp 2 dcs oak dn SS Bae e ar Aaa ech 25 example jump 2d S dd eo eme eth a x ceptus aus 26 example jump 2da da la vd ER A 26 example flowe2d 4 2 4s RAEE EE NAAA 28 examples OW ue u ap na do Oe v O A do WIR Yo Mee IS 28 example flow 3d uuu era X omnee Vie ER e UE en X A en 29 example skull coarse grid a aa 30 example skull refined grid 34 ww ef Er es ras 31 example cylindrical mesh y cerros ee 32 example axisymmetrical problem 1 2 222 2 nennen 32 example axisymmetrical problem 2 2 2 2 2 nn nenn 33 example obstacle problem 0 000 epee eee 35 example Stefan problem ars quie eta 2 3 4 2m AA Ee 36 example porous media La Dis rare bees 38 partition of the unit interval a de a doas ogg cd urgeat de geld 40 partition of the unit squares 23024 dais Dar S oa 42 circlular geometry sa I cr hod Sd 43 2D region milk holes a ps qtu eC RR he 48 3D region with holes rar Axa mx as mec LE ta ura 51 LIST OF FIGURES 3l 32 33 34 35 7 Tver 1 OQ os A AS AA M EORUM ERA M EIU 51 cylindrical geometry tt ue cet dora mg pe he age 53 userdefined proble ss 49 0 ire S NER a a 63 graphical user interface 1 er Re A all 72 gra
44. figure 4 we show a cut through the finite element mesh provided after some steps of the adaptive solver This picture and some of the following were produced by GRAPE For 3D geometries there is no online graphics integrated in KASKADE We recommend to use one of the popular visualization programs GRAPE EXPLORER or AVS and to use our interfaces to write data on a file in an appropriate format 4 Static problem using the files unit 1d cmd unit 1d geo unit 1d mat Poisson equation with constant coefficients on the unit interval 0 0001u 5 5 in 0 1 u 0 in the points 0 and 1 The coefficient of uz is arbitrary and chosen so mysteriously only to illustrate the use of keywords e g Factors in the mat file Compare the chapter 5 and the file unit 1d mat 5 Static problem using the files unit 2d cmd unit 2d geo unit 2d mat Poisson equation with constant coefficients on the unit square Urt Uy 1 in 0 1 x 0 1 u 10 0 on the boundary 6 Static problem using the files unit 2d a cmd unit 2d a geo unit 2d a mat IGRAPE Graphical Programming Environment Copyright Sonderforschungsbereich 256 Uni versity of Bonn Germany Copyright The Numerical Algorithms Group Ltd Oxford UK 3Copyright Advanced Visual Systems Inc 3 STARTING THE PROGRAM 19 Figure 4 Cut through the adaptive mesh in example 3 Here we solve the same elliptic equations as in example 5 but with additional Neumann and Cauc
45. for timestep which fulfills timestep postTimeStep 0 and time step gt 0 select one of the preconditioners print linear system print multi level matrices print stiffness matrix A in a format reada ble by MatLab print triangulation lists print command parameters specifies the problem type to be allocated and solved B COMMAND LISTING 86 Parameter Name Default Description set in file kaskade init read0ldFormat O false read the geometry input in old ZIB format readZIBFormat 1 true read the geometry input in new ZIB format refStrategy extrapolation select the refinement strategy refTetrahedron shortestEdge method for refinement of a tetrahedron scaleX scaling x coordinate for quasi 3d plots scaleY scaling y coordinate for quasi 3d plots scaleZ scaling z coordinate for quasi 3d plots spaceDim space dimension time false print cpu time for all essential modules timeAssembler false print cpu time for matrix assembling timeErrorEstimator false print cpu time for error estimation timeLinSystem false print cpu time for system solution timeRefinement false print cpu time for mesh refinement O O O O O O N e FP n false print cpu time for data transport between different grids timeTransport timeUpdate false print cpu time for update operation in node management writeAVSFormat false set format for geometrical output it cor responds to AVS visualization tools only for 3
46. ggestion to improve the code or this guide Please send your contribution by e mail to erdmann zib berlin de or roitzsch zib berlin de Please note that there is another way to exchange experiences about KASKADE the mailing list kaskade 10zib berlin de You can subsribe by sending a mail to the listserver Majordomo zib berlin de with the message subscribe kaskade 1 in the body By this list you are always informed on essential news about the code and its appli cations and you can use this channel for discussing questions of general interest 2 Installation 2 1 Obtaining the Source Code The KASKADE source code is part of the CodeLib a collection of numerical codes developed at the Konrad Zuse Zentrum It is available in the electronical library elib 2 INSTALLATION 9 by anonymous ftp gt ftp elib Internet 130 73 108 11 gt username anonymous gt password e mail address gt cd pub kaskade 3 x In this directory you find a compressed tar file 3 x tar Z and a README with some hints for installation x stands for the number of the newest version By the commands gt binary get README gt get 3 x tar Z you fetch copies into your local directory The tar file 3 x tar Z contains the source code and examples Use the following commands to extract these files gt uncompress 3 x tar Z gt tar xf 3 x tar 2 2 Compiling and Linking the Code The make file is kaskade make There are four
47. hlet and one with Neumann condition compare figure 8 Additionally we approximate a circular boundary during the refinement process du n 0 0 on the boundary y 0 1 gt 0 Neumann u 10 elsewhere on the boundary Dirichlet In figure 9 you see the mesh and the isolines of the solution after adaptive refinement 9 Static problem using the files s1it 2d a cmd slit 2d a geo slit 2d a mat Consider Laplace s equation A u 0 0 3 STARTING THE PROGRAM 0 1 4142 D Dirichlet N Neumann D D 1 1 1 1 1 1 D D 1 4142 0 N Poit S D 1 1 Cep E D 0 1 4142 Figure 8 Initial triangulation in example 8 Min 1 Max 1 352 22 Figure 9 Mesh and solution after adaptive refinement steps in example 8 3 STARTING THE PROGRAM 23 Min 0 Max 1 091 Figure 10 Contour lines of the solution in example 9 on a circle of radius 1 centered at the origin with a slit along the positive x axis Homogeneous Dirichlet boundary conditions are imposed on the top and homogeneous Neumann boundary conditions at the bottom of the slit The solution is singular and has the leading term u r sin JA This function is used as Dirichlet boundary condition on the remaining part of the boundary yielding it as exact solution The associated contour plot is given in figure 10 Starting with a simple polygonial region we approximate the circular boundary during the refinement pro
48. hy boundary conditions 10 0 on the boundary z 0 or y 0 Dirichlet uy 0 0 on the boundary y 1 Neumann u Ustu 20 on the boundary x 1 Cauchy Figure 5 illustrates the grid and solution after some refinement steps Static problem using the files slit 2d 45 cmd slit 2d 45 geo slit 2d 45 mat The problem is defined on the polygonial region in figure 6 by the same equation as in example 5 and the boundary conditions Ou On u 10 elsewhere on the boundary Dirichlet 0 0 on the boundary y 0 2 gt 0 Neumann 3 STARTING THE PROGRAM Figure 5 Mesh and solution after some refinement steps in example 6 0 1 4142 D D 1 1 eun 1 1 D Dirichlet N Neumann D D N 1 4142 0 1 4142 0 D D 1 1 1 D D 0 1 4142 Figure 6 Initial triangulation in example 7 3 STARTING THE PROGRAM 21 Min 1 Max 1 31 Figure 7 Mesh and solution after 13 refinement steps in example 7 The adaptive refinement generates a fine mesh in the neighbourhood of the origin where the solution has a singularity In figure 7 we show the mesh and the isolines of the solution after 13 refinement steps 8 Static problem using the files slit 2d cmd slit 2d geo slit 2d mat The same problem as in 7 but on a different domain Here the angle of the slit has been set to 0 degree That means we have two edges from 0 0 to 1 4142 0 one with Diric
49. imation Deuflhard Leinen Yserentant 1989 see DLY89 and adaptive refinement Numerical parameters 11 GRAPHICAL USER INTERFACE 75 Graphic Parameters Fa ZIB GUI Kaskade 3 0 Example Settings Information Size Plot Window 4 sl Start Step Pause Resume Finish a Show Indices K6 in mome newton bzferdma O oK oo Point TF Graphics Fa Numerical Parameters FI Linear System Triangle Preconditioner Convergence Test CG SSOR A Residual Bicostab v Jacobi a E Problem Definition yes ym Type Material Dirichlet GMRes BPX Extend Ite Prec x x StaticHeat Constant Constant y NonLinHelax HB 1 0 ES TransientHeat v PeakSource w Hexagon x TrcNonLinMLGS Frosmoottistega T se PorousMada Y EA lapa KC NORALGS PostSmoothSteps 1 wv Obstacle MultiPeakSource v Step S NONE MGNodeProgFactor 3 0 x Stefan w StefanSource x TransPeak x Nome 7 Local multigrid y Stefan I Local smooth on top level x Porous Media Extend smooth pattern Geo file unit 2d geo Mesh Refinement Mat file Error Estimator Refinement Strategy Min Ref Ratio Dimension 2 DLY gt MaxValue 0 05 global precision 1e 3 s FastDLY ExtraPolation minimal level O wv AKI x Uniform maximal level 50 wv None AAA AE Figure 35 Menus for setting the problem parameters linear solver cg preconditioner MLsgs error estimator DLY
50. ime with a few exceptions e g for output routines e The symbol starts a comment up to the rest of the line in any input file e The main command file kaskade init is always read as the first It sets default values for most commands e Parameter values are overwritten by subsequent commands 3 STARTING THE PROGRAM 12 Read Coarse Initial Grid Y Assemble and Solve Linear System 7 y Estimate Discretization Error Postprocessing Yes BreakCondition No Y Refine Grid Figure 1 Main iteration loop of the adaptive solution procedure e If the characters lt CR gt appear on the output screen at run time the function Pause has been activated by the command pause and the program expects input from the user see the description of pause in table 1 Some important parameters are shown in Table 1 Any problem like the one specified by unit 1d cmd expects a geometry definition file e g unit 1d geo and a material file e g unit 1d mat 3 2 Algorithmic Kernel The finite element methods in KASKADE use adaptive multilevel techniques to achieve of optimal complexity see DLY89 Here we give a short outline of the algorithm in the elliptic case which is selected by the command problem staticHeat 1 The user has to define the following input e Coefficients of the equation in a mat file for piecewise constant values or he has
51. in is x y 0 lt v y lt 10 1 y l4 lt r lt 6A2 lt y lt A4 Ha yd lt r lt 6A6 lt y lt 8 5 GEOMETRY AND MATERIAL PROPERTIES 46 and one coarse triangulation is given by the following description see also figure 28 holes 2d geo Points maxIndex 24 1 0 0 0 0 4 0 0 0 6 0 0 0 10 0 0 0 oo oo 5 WN Rp RO HAHAHAHAHA O MOAN OOF CO th Oi ey Oe opes ax Ooo o 00 00 00 gt O O OO o 10 0 8 0 21 0 0 10 0 22 4 0 10 0 23 6 0 10 0 24 10 0 10 0 END N o Elements 151 NO oP J NNN oo 000 O 0 YX Y O O 0 HAHAHA 8 12 11 1 9 10 13 1 10 14 13 10 11 14 11 15 14 11 12 15 12 16 15 13 14 17 14 17 18 15 19 16 16 20 19 17 21 18 18 21 22 18 19 22 19 22 23 19 20 23 20 23 24 END Boundary i214 23d1 34d1 48d1 8 12d 1 12 16 16 20 20 24 24 23 23 22 aaaaaaAaQA RE 14 15d 3 15 19d 3 19 18 a 3 5 GEOMETRY AND MATERIAL PROPERTIES BPRPRPPP RPP RPP ERP PR EP EE AT 5 GEOMETRY AND MATERIAL PROPERTIES 48 EINE Figure 29 Triangulation of a quadratic geometry with two holes 18 14d 3 END e 3D region comments Points maxIndex mazIndex index of point lt zx y 2 coordinates of point lt index of point lt zx y 2 coordinates of point end Elements indices of points of element classification indices of points of element classification END Boundary indices of points of boundary triangle type classifi
52. k Free and Moving Boundary Problems Oxford University Press Oxford 1988 DLY89 P Deuflhard P Leinen and H Yserentant Concepts of an Adaptive Hier archical Finite Element Code IMPACT Comput Sci Engrg 1 1989 EO82 C M Elliot J R Ockendon Weak and Variational Methods for Moving Boundary Problems Reseach Notes in Mathematics 53 Pitman London 1982 GRA GRAPE Graphical Programming Environment Sonderforschungsbereich 256 University of Bonn Germany HS90 K H Hoffmann and J Sprekels editors Free Boundary Value Problems Birkh user Basel 1990 HYP89 P Wust J Nadobny R Felix P Deuflhard W John A Louis Numerical approaches to treatment planning in deep RF Hyperthermia Strahlenther Onkol 165 p 751 757 1989 Kor95 R Kornhuber Monotone Multigrid Methods for Nonlinear Variational Pro blems Habilitationsschrift Freie Universit t Berlin Berlin 1995 Mus37 M Muskat The Flow of Homogeneous Fluids Through Porous Media McGraw Hill New York 1937 REFERENCES 88 NPR96 U Nowak U P hle R Roitzsch Eine graphische Oberfl che f r nume rische Programme Internal report Konrad Zuse Zentrum f r Informationstech nik Berlin 1996 Available as postscript file zgui german ps Z in subdirectory pub kaskade Manuals 3 x by anonymous ftp to machine elib zib berlin de NPRW 96 U Nowak U P hle R Roitzsch R Werk ZGUI Manual Internal report Konrad Zuse Zentrum f
53. me refinement steps In chapter 5 we explain how we handle arcs of a circle as boundaries 3 STARTING THE PROGRAM Figure 15 Initial mesh and mesh after some refinement steps in example 15 Figure 16 Mesh and contour map of the solution in example 16 28 3 STARTING THE PROGRAM Figure 17 Mesh and contour map of the solution on a cutting plane in example 17 17 Static problem using the files flow 3d cmd flow 3d geo flow 3d mat 29 18 As in example 16 we solve an elliptic equation with a mass term and only Neumann condition on the boundary U Uyy Uz u 1 in a ball of radius v3 Ou On 1 0 on the boundary Figure 17 shows the mesh and the isolines of the approximate solution on a cutting plane after some refinement steps In the chapter 5 about geometrical input we explain how we handle curvilinear boundaries Static problem a calculation from Hyperthermia In order to demonstrate the applicability of our code to real life problems we consider a problem arising in hyperthermia Hyperthermia is a cancer therapy based on the observation that local heating may slow down the growth of a tumor especially if it is applied in combination with other methods like chemotherapy or radiotherapy The heating of the tissue is obtained by radiowaves generated by an antenna array The antennas are either fixed on the skin or implanted in the tissue itself Of course the 3
54. meter Name Default Description set in file kaskade init accTime false print all accumulated cpu times absPrecision false see parameter globalPrecision batchJob false the program terminates automatically when fulfilling a break condition cmd input command file compareSolution false compare approximate solution with the true solution in the grid nodes cycle false new boundary points are placed on a circle 2d or ball 3d around the origin cylinder false new boundary points are placed on a circle around the origin parallel to the plane z 0 directSolverLimit limit for the direct sparse matrix solver on levels gt 0 see also levelOdirect dirichletBCs constDirichlet identifies the type of Dirichlet boundary condition errorEstimator DLY select an error estimator extPrecFactor 1 0 external precision factor manipulates re quested precision for convergence test in iterative solvers unit 2d geo name of geometry file the extension geo need not be specified The material file is expected to have the same name but the extension mat globalPrecision desired relative precision maximum dis cretization error with respect to global energy if absPrecision 0 otherwise de sired absolute precision graphic 0 false switch graphic support on off info O false all information levels are activated or de activated infoAb 0 false general info about linear system B COMMAND LISTING
55. n You only have to define the functions for the space dimension of your problem Furthermore it may be that you want to com pare the approximate solution with a given function e g for a simple test problem the known true solution Such a function can be defined in the member function trueSollnPoint and will be used if the parameter compareSolution is set to 1 in your cmd file In userStatic cc we realized the following examples e 1 D 8 HOW TO DEFINE A NEW PROBLEM 63 0 1 Cauchy boundary 1 1 Dirichlet E a d 3 E E 3 3 E E 2 e s z Ss 5 E A Z 0 0 Dirichlet boundary 1 0 Figure 33 Geometry and boundary condition in example user_ static 2d VeVu zu 1 4x 4z x 1 in Q u emo on Pp 25 z 2x 1l on Py e 2 D VayVut atyju r y rcty z z y y in Q ry on Tp u 1 in 1 1 o zy xy on Py On Ou 2 ay trxu 2x y on Po On I is defined in the file user static 2d geo see figure 33 and the following description unit square Points maxIndex 4 1 0 0 0 0 8 HOW TO DEFINE A NEW PROBLEM 64 e UN Orere ooo BRO ooo end Elements 123 1 134 1 END Boundary 12 D1 Dirichlet 23 1 i Neumann 34 2 Cauchy 4 1 1 Dirichlet End vos Points 3D1 Dirichlet END e 3 D VeVu ze y z u zxyzlr y 2 yz in Q w us on Pp du da DU on Py On As an example we give the definition of dummy functions for the elliptic
56. nd 5 GEOMETRY AND MATERIAL PROPERTIES 56 Factors e le 4 end The first integer in each row is the classification of the region or boundary It cor responds to the classification in the geometry file set for each element of the region or each face of the boundary The classification parameter must be an integer greater than 0 In our example we define values for boundary sections of classification 1 and 2 On sections with classification 1 or 2 and boundary type Dirichlet D we have the constant values 10 and 12 On sections with boundary type Neumann N we have homogeneous or inhomogeneous Neumann conditions On elements marked with 1 the elliptic coefficient is k 1 and there is no source f 0 or mass term q 0 On elements marked with 2 we have an elliptic coefficient k 2 a constant source f 1 and no mass term q 0 Material terms which are set to zero may be omitted If the keyword Factors occurs then the values following after one of the keywords e elleptic s source m mass co convection c Cauchy n Neumann or d Dirichlet are used for scaling the corresponding values of the material or boundary defined in the Materials section before We introduced this scaling to simplify the handling of physical units The next example compare the preceding example 5 shows how to handle Cauchy boundary conditions We have the unit square as geometry described in the file unit 2d A geo
57. nguish this Neumann part of the Cauchy boundary from other Neumann boun dary faces we have to use a common class identifier to the boundary face in the geo 8 HOW TO DEFINE A NEW PROBLEM 66 file Though we make no use of it in our examples we provided a function Real UserStaticMaterial Inner which is only evaluated if the command pa rameter innerBoundary is set to 1 This is necessary if you want to compute an additional term of the form 3 appendix A Please compare the examples 13 and 14 in chapter 3 The point 1 1 is declared to be of type Dirichlet in the user static 2d geo file The values for this point and the other Dirichlet points are computed by the function Real UserDirichletBCs setBC For invoking your user problem of type 1 in appendix A you just have to set the parameters e problem staticHeat e material userStatic e dirichletBCs userStatic We provided cmd files with these commands user static 1d cmd user static 2d cmd and user static 3d cmd The geo files are named user static 1d geo user static 2d geo or user static 3d geo Of course you can choose other names In user Transient cc you find the following 2 D problem a V Vutu 2 0 x 1 0 y y 1 0 exp t in Q u zx x 1 0 y y 1 0 exp t on Pp u r y 0 x x 1 0 y y 1 0 in Q A problem of type 4 in appendix A is invoked by the commands e problem transientHeat e material userTransien
58. o PP ww anaana PrPrRPrPR EB nn nm nn an nm n END The character s marks the edge for being approximating as circular arc during further refinement steps by placing the coordinates of a new point on this edge onto a circle centered in the origin x 0 y 0 through the boundary points of the edge Omitting the specification of the boundary shape or the character P provides the usual calculation of the coordinates as midpoint between the two boundary points of an edge Though our implementation of curvilinear boundaries is rather special e g cir cular arcs are centered in the origin it does not seem very difficult to generalize the concept As an example for a more complicated geometry we present the geo file of example 9 among the static problems slit 2d a geo Points maxIndex 10 1 1 000000e 00 1 000000e 00 2 0 000000e 00 1 414213562e 00 3 1 000000e 00 1 000000e 00 4 1 414213562e 00 0 000000e 00 5 0 000000e 00 0 000000e 00 5 GEOMETRY AND MATERIAL PROPERTIES 45 1 414213562e 00 0 000000e 00 1 000000e 00 1 000000e 00 0 000000e 00 1 414213562e 00 1 000000e 00 1 000000e 00 10 1 414213562e 00 0 000000e 00 END o ON OD Elements 125 1 235 1 aa ana an w PN OOF OD Riba 0 c co ee END Boundary N a ARANHA nun nun un un un RH XOwoQgmn0anaNn Hr r o Qaaqaaadqapaa EN CO o END Points 6 D2 END We present one more 2D example showing how to define regions with holes The doma
59. ot window 11 GRAPHICAL USER INTERFACE 73 There are default values which configure a simple static heat conduction problem on the unit square in 2 space dimensions Meanwhile we added a Postscript button to the application control board which can be used to generate postscript pictures of the actual mesh and approximate solu tion In the root window you find the Example menu It offers a set of prepared examples Each example is a collection of parameter settings stored in a file extension ex To create your own example file you can use the saveAs or save field in this menu By save you can overwrite only the example you selected before Maybe you have overwritten an example but want to reestablish the example settings recommended by the authors then you can select Default in the Example menu It may be that for a special problem type the selection of some other parameters materials boundary conditions or numerical methods makes no sense In this case these widgets should be disabled Changes in one of the settings windows have to be confirmed with the Apply button before being valid for the next run The Default button resets all options to the default values you get when you start the program By the Reset button the state before the last Apply click is reestablished A path through the program e Select the type of your problem In the user interface we use the name staticHeat synonymous to linear scalar elliptic problem an
60. phical user interface 2 cs da x e ale EI de deus ES T5 1 INTRODUCTION 8 1 Introduction The finite element code KASKADE version 3 x solves some types of linear elliptic or parabolic equations and some nonlinear problems A more detailed description of the related problem classes is given in appendix A Our main goal is to present an easy way to customize this software package for solving a wide range of partial differential equations Of course there are a lot of applications of similar problem types which cannot be handled directly with the included algorithms But we hope that our additional advices for extending the code by numerical methods or problem classes will help to use the code as a base for the development of new prototypes After reading the first two chapters the reader should be able to install the code and to run a set of provided examples In the following chapters he can learn how to implement his own problem This user s guide does neither provide a treatment of mathematical concepts nor does it describe the object oriented implementation of these Such issues are covered in the technical report BER95 We are extending this tutorial continuously The latest version is available by network as described in the next chapter If there are any users extending KASKADE we would be pleased to get notice of their work If possible we will publish their solution in the appendix The authors are very thankful for any su
61. r Informationstechnik Berlin 1996 Ouster94 J Ousterhout Tel and the Tk Tollkit Addison Wesley Publishing Com pany Professional Computing Series 1994 Pen48 H H Pennes Analysis of tissue and arterial blood temperatures in the resting human forearm J Appl Physiol 1 p 93 122 1948 Rodr87 J F Rodrigues Obstacle Problems in Mathematical Physics Mathematical Studies 134 North Holland Amsterdam 1987 Roi93 B Erdmann J Lang and R Roitzsch KASKADE Manual Version 2 0 FEM for 2 and 3 Space Dimensions Technical Report TR 93 5 Konrad Zuse Zentrum f r Informationstechnik Berlin 1993 See90 M Seeba 3D Computersimulation der interstitiellen Mikrowellen Hyperthermie von Hirntumoren Bericht Nr 1 90 Inst f Radiologie und Patho physiologie Deutsches Krebsforschungszentrum Heidelberg 1990 Son89 P Sonneveld CGS a fast Lanczos type solver for nonsymmetric linear sy stems SIAM J Sci Stat Comp 10 1 1989 SS86 Y Saad and M H Schultz GMRES a generalized minimal residual algorithm for solving nonsymmetric linear systems SIAM J Sci Stat Comp 7 3 July 1986 vdV92 Henk A van der Vorst BI CGSTAB a fast and smoothly converging variant of BI CG for the solution of nonsymmetric linear systems SIAM J Sci Stat Comp 13 3 March 1992
62. refinement strategy extrapolation Variation of parameters change the region by the geo file button e g select unit 1d geo unit 3d geo hexagon geo 1 2d geo 1 3d geo slit 2d geo or corner geo The space dimension of the problem is automatically provided by analyzing the the geo file Then the program selects the appropriate numerical algo rithms Of course the material data must be reasonable for the geometry e g if the region includes a Neumann boundary the values of it must be defined in the corresponding mat file if values are piecewise constant or in a member function of the selected material if values are not piecewise constant We recommend to use always a mat file with the same root name as the geo file 11 GRAPHICAL USER INTERFACE 76 change the coefficients in the equation e g select Material PeakSource peak in the source term or Material MultiPeakSource change the numerical parameters in the Numerics window e g select different linear solvers and combine them with some preconditioners or test different refinement strategies For example a computation without error estimation error estimator None or with uniform refinement refinement strategy uniform is quite interesting if you have selec ted the problem with a peak source Material peakSource In the uniform case you need 5 refinement steps and about 2000 nodes to achieve an accuracy of 296 In the adaptive case you need
63. se the values in the nodes are ignored 5 GEOMETRY AND MATERIAL PROPERTIES 54 5 3 Material Coefficients Boundary Conditions and Initial Values We select the coefficients usually describing material properties of the equation the Neumann and the Cauchy boundary conditions by the command material where represents one of the keywords in the left column of table 3 Keyword materials DefaultMaterial all coefficients are constant PeakSource coefficients in the differential operator are constant the peak shaped source function depends on spatial coordina tes e g used in example peakld cmd MultiPeakSource similar to PeakSource but a different source function TransPeakSource in transient problems variable source function all the other coefficients are constant StepSource in transient problems all coefficients are piecewise constant e g used in example u2step cmd StefanSource variable source function all coefficients are piecewise con stant e g used in example stefan cmd UserStaticMaterial all coefficients and boundary conditions are variable and gi ven by a function e g used in example user static 1d cmd UserTransientMaterial all coefficients and boundary conditions are variable and gi ven by a function e g used in example user trans 1d cmd Table 3 Keywords for material types The boundary condition is chosen by the command DirichletBCs where represents one of the keywords in
64. ss 9 Obtaining Run Time Information In the tables 11 12 and 13 we list commands which are relevant for getting runtime information 10 Technical Details for Programmers 10 1 Vector and Matrix Classes We realized vector and matrix structures as special classes allowing index control during runtime This is a very helpful feature during the development of new code because the violation of array bounds is one of the most frequent errors Of course 10 TECHNICAL DETAILS FOR PROGRAMMERS 68 Parameter Name info infoAb infoErrorEstimator infoLinSystem infoMG infoPrecond infoRefinement infoSolution printAb printAL printMatLabFormat printMesh printParameters Default set in file kaskade init O false O false O false 25 false false true false false false false false Description all information levels are activated or de activated general info about linear system info about error estimation print information on every 25th iteration step info about multi grid preconditioning info about preconditioning mesh info levels ranging from 0 to 2 gi ving statistics about refinement steps info after each solution step print linear system print multi level matrices print stiffness matrix in a format reada ble by MatLab print triangulation lists print command parameters Table 11 Info and print commands 10 TECHNICAL DETAILS
65. t e dirichletBCs userTransient You find sample cmd files setting these parameters user trans 1d cmd user trans 2d cmd and user trans 3d cmd The geo files are named user trans 1d geo user trans 2d geo or user trans 3d geo 9 OBTAINING RUN TIME INFORMATION 67 8 2 Something is variable something is constant In the last section we made the assumption that all the functions defining problem coefficients are variable This allows an easy handling of the program similar to the case of constant functions Of course you can mix problems of the form constant Dirichlet boundary and variable material coefficients including Neumann and Cauchy boundary or vice versa You have to select either dirichletBCS ConstDirichlet and material userStaticMaterial or dirichletBCS userDirichlet and material DefaultMaterial It is a bit more complicated if some coefficients are piecewise constant and others are variable For example we derived a material type UserVarSource for problems with piecewise constant coefficients k q or c and a variable right hand side This material can be selected by command material userVarSource The source function f has to be specified in the class UserVarSource within file materialA cc The other coefficients the Neumann or Cauchy boundary are ex pected to be constant with values defined in the mat file Other constellations of coefficients can be handled analogously by implementing a new material cla
66. ter for preconditioning with Jacobi or SSOR smoothing B COMMAND LISTING Parameter Name plotBoundary plotElements plotSolution plot3d plotIsoLines plotKeep plotLevels plotPointNodes plotSize plotTimeStep plotTriangleNodes postScript postTimeStep preConditioner printAb printAL printMatLabFormat printMesh printParameters problem Default set in file kaskade init 1 true true false true false true false false 4 false false MGsgs 0 false 0 false false O false O false staticHeat 85 Description stop when function Pause is called in the code and wait for input lt CarriageReturn gt continue until next Pause is encountered c continue and disable the function Pause p gt generate a picture in postscript for mat of the actual mesh and the approxi mate solution q gt quit program plot no boundary when solution is plotted plot no elements when the solution is plotted plot solution data after each space step plot 2D results as a surface plot quasi 3d plot plot isolines of solution keep all plots in separate windows desired number of iso lines plot number of nodes at points size of plot window plot the solution after each time step plot number of nodes at triangle generate postscript pictures of mesh and solution generate postscript pictures of mesh and solution
67. ternal precision factor manipulates re quested precision for convergence test in iterative solvers print information on every 25th iteration step limit for the direct sparse matrix solver on level 0 i e the solver will be used for a matrix dimension up to 3000 determines the type of iterative solver maximal number of vectors to be orthogo nalized in GMRES relaxation parameter for preconditioning with Jacobi or SSOR smoothing select one of the preconditioners print cpu time for system solution Table 6 Parameters for linear solvers and preconditioners 7 ERROR ESTIMATION AND MESH REFINEMENT 60 Keyword Preconditioner Jacobi Jacobi type preconditioner SGS symmetric Gauss Seidel preconditioner ILU preconditioning by incomplete LU factorization MLJacobi multilevel preconditioner with Jacobi type smoothing The default relaxation parameter w is 2 3 MLSGS multilevel preconditioner with symmetric Gauss Seidel smoothing Forward Gauss seidel is applied for pre smoothing the backward operation for post smoothing AMLJacobi Additive multilevel preconditioner with Jacobi type smoo thing also called BPX preconditioner BPX90 AMLSGS Additive multilevel preconditioner with Gauss Seidel smoothing Table 7 Preconditioners for linear problems Keyword Preconditioner nonlinSGGS single grid Gauss Seidel preconditioner nonlinMLGS multi level Gauss Seidel preconditioner trcNonlinMLGS truncated multi level
68. th points 1 2 and 3 and classification 1 Dirichlet boundary with points 1 and 4 classification Dirichlet boundary with points 2 and 3 classification Dirichlet boundary with points 1 and 2 classification 1 e point with index 1 x 0 0 y 0 0 point with index 2 x 1 0 y 0 0 Figure 27 Partition of unit square 2D unit square Points maxIndex 4 1 0 0 0 0 e WON Orr ooo ereo ooo end Elements 123 1 134 1 END Boundary 12 D1 e WA e Ae 0 UU J hip 5 GEOMETRY AND MATERIAL PROPERTIES 43 Figure 28 Approximation of a circular geometry by 3 uniform refinement steps END In the 2D case we implemented a simple mechanism to approximate a circular geometry The refinement algorithm computes new boundary points on the circle centered in the origin After some steps of refinement the polygonal region is very similar to a circle As an example look at figure 28 This special handling of boundary edges is initiated by specifying an additional character for each curvilinear boundary edge see for example the parameter file circle cmd with the geometry description in circle geo or circle 2d geo circle geo Points maxIndex 10 1 1 0 0 0 2 0 5 0 866025404 3 0 5 0 866025404 4 1 0 0 0 5 0 5 0 866025404 6 0 5 0 866025404 5 GEOMETRY AND MATERIAL PROPERTIES 44 7 0 0 0 0 end Elements 2 7 AHHH 1 2 3 4 5 6 FO A 0 NNNNWN END Boundary N Q OUO AUNE H
69. the left column of table 4 In general the coefficients of the equation the Dirichlet boundary values and the initial values in transient problems are provided by member functions of the classes Material or DirichletBCS which are implemented in materials cc materialsA cc dirichlet cc dirichletA cc In the case of piecewise constant material coefficients or piecewise constant Dirichlet boundaries it is not necessary to modify the source code The user only has to set the 5 GEOMETRY AND MATERIAL PROPERTIES 55 Keyword dirichletBCs ConstDirichletBCs the solution is piecewise constant on the Dirichlet boundary UserDirichlet the Dirichlet boundary values are given by a function UserTransient the Dirichlet boundary and initial values are given by a function Table 4 Keywords for Dirichlet boundary types parameters material defaultMaterial and dirichletBCs constDirichlet Then the description of these values is read from a mat file in which the material coefficients and the boundary conditions are specified like in the following example The classifi cations of the region or boundary have to correspond to those of the geometry file If there is no mat file specified by the command matFile x x x we suppose as default name the root name of the geo file but with the extension mat This is a simple example for a mat file Boundary 1 D 10 0 1 N 0 0 2 D 12 0 2N 1 0 end Materials isotropic lelsOend 2e2slend e
70. this conversion has also to be done for the Cauchy boundary conditions 3 3 2 Transient Problems The following examples demonstrate the parabolic solver which is adaptive in space and time In each time step the computation of the spatial solution starts from the initial grid 1 Transient Problem using the files u id step cmd u 1d step geo u 1d step mat The solution of this parabolic problem shows a transient diffusion process star ting from a sharp initial profile 1000u Au 100u 0 0 oe 1 0 for x lt 0 5 ze for x gt 0 5 The coefficients are described in the file u 1d step mat the domain unit in terval in u 1d step geo The boundary conditions are u 0 in x 0 and u 1 in x 1 The initial values are provided by a function in the file dirichlettr cc 2 Transient Problem using the files u 2d step cmd u 2d step geo u 2d step mat Analogous to example 1 but on the unit square The boundary condition are u 0 on the left side u 1 opposite to it and homogeneous Neumann condition on the other parts of the boundary 3 Transient Problem using the files u 3d step cmd u 3d step geo u 3d step mat Analogous to example 2 but on the unit cube 3 STARTING THE PROGRAM 35 Min 0 Max 0 4709 TE Tir Ir yl m 8 T
71. tive techniques for ordinary differential equations and elliptic boundary value problems are combined R Kornhuber contributes the algorithms for problems formulated as variational ine qualities see appendix A and Kor95 3 3 Examples The KASKADE package includes some example problems each of them is defined by setting parameters in a command file extension cmd Just type gt k6 cmd x xkx to compute one of them with the executable k6 The stands for the name of a command file e g gt k6 cmd unit 1d cmd In the next chapter you learn how to formulate these command files 3 3 1 Static Problems 1 Static problem using the files peak 1d cmd peak 1d geo peak 1d mat Static heat conduction of type 1 in appendix A on the one dimensional unit interval Sx uf in 0 1 u 0 in the points 0 and 1 This Poisson problem has the constant coefficient k 1 in the Laplacian and a source function f exp 100 0x x 0 5 with a central peak and zero values on the boundary points The solver stops when the relative global error of the approximated solution is less than 1 0e 5 relative to the energy norm The conjugate gradient method that solves the linear systems is preconditioned by a multigrid algorithm The adaptive refinement is tuned by an error estimator Figure 2 shows the final grid In this and in the two following examples we know the true solution of the problems and by setting the paramet
72. unit square Points maxIndex 4 1 0 0 0 0 e NM Orere ooo BRO ooo end Elements 123 1 134 1 END 6 SYSTEM SOLUTION AND PRECONDITIONING 57 Boundary 12 D1 Dirichlet 23 C3 Cauchy 34 N2 Neumann 41 Di Dirichlet END The corresponding values of the boundary conditions are set in the unit 2d A mat file Boundary 1 D 10 0 Dirichlet 2 NO0 0 Neumann 3 C1 0 2 0 Cauchy end Materials isotropic 1 el s1 end end The Cauchy boundary condition is given by two values compare 1 in appendix A the first is a the second is the right hand side qc Here we have u 10 0 on the boundary x 0 or y 0 Dirichlet uy 0 0 on the boundary y 1 Neumann Ustu 2 0 on the boundary z 1 Cauchy The values specified in the mat file are only relevant for the program in case of constant materials or constant boundary conditions In one of the following sections we explain where to define nonconstant coefficients or boundary values in the source code 6 System Solution and Preconditioning 6 1 The Iterative Solvers We have implemented several iterative solvers Anyone of these is invoked by the command 7 ERROR ESTIMATION AND MESH REFINEMENT 58 linSolver where represents one of the keywords in the left column of table 5 Keyword Solver cg conjugate gradient solver cgODir conjugate gradient solver with 3 term recurrence cr conjugate residual solver cr Dir conjugate residual solver with 3

Rudolf Beck Bodo Erdmann Rainer Roitzsch KASKADE User's Guide

Contents

Download Pdf Manuals

Related Search

Related Contents