Numerical Methods and Diffpack Programming

## Contents |

The target audience of this book is students and researchers in computational sciences who need to develop computer codes for solving partial differential equations. The exposition is focused on numerics and software related to mathematical models in solid and fluid mechanics. The book teaches finite element methods, and basic finite difference methods from a computational point of view. The main emphasis regards development of flexible computer programs, using the numerical library Diffpack. The application of Diffpack is explained in detail for problems including model equations in applied mathematics, heat transfer, elasticity, and viscous fluid flow. Diffpack is a modern software development environment based on C++ and object-oriented programming.

During the last decades there has been a tremendous advancement of computer hardware, numerical algorithms, and scientific software. Engineers and scientists are now equipped with tools that make it possible to explore real-world applications of high complexity by means of computer simulation. Experimentation based on numerical simulation has become fundamental in engineering and many of the traditional sciences. A common feature of mathematical models in physics, geology, astrophysics, mechanics, geophysics, as well as in most engineering disciplines, is the appearance of systems of partial differential equations (PDEs). This text aims at equipping the reader with tools and skills for formulating solution methods for PDEs and producing associated running code.

Successful problem solving by means of mathematical models in science and engineering often demands a synthesis of knowledge from several fields. Besides the physical application itself, one must master the tools of mathematical modeling, numerical methods, as well as software design and implementation. In addition, physical experiments or field measurements might play an important role in the derivation and the validation of models. This book is written in the spirit of computational sciences as interdisciplinary activities. Although it would be attractive to integrate subjects like mathematics, physics, numerics, and software in book form, few readers would have the necessary broad background to approach such a text. We have therefore chosen to focus the present book on numerics and software, with some optional material on the physical background for models from fluid and solid mechanics.

The main goal of the text is to educate the reader in developing simulation programs for a range of different applications, using a common set of generic algorithms and software tools. This means that we mainly address readers who are or want to become professional programmers of numerical applications. As the resulting codes for solving PDEs tend to be very large and complicated, the implementational work is indeed nontrivial and time consuming. This fact calls for a careful choice of programming techniques.

During the 90s the software industry has experienced a change in programming technologies towards modern techniques such as object-oriented programming. In a number of contexts this has proved to increase the human efficiency of the software development and maintenance process considerably. The interest in these new techniques has grown significantly also in the numerical community. The software tools and programming style in this book reflect this modern trend. One of our main goals with the present text is in fact to explore the advantages of programming with objects in numerical contexts. The resulting programming is mainly on a high abstraction level, close to the numerical formulation of the PDE problem. We can do this because we build our PDE solvers on the Diffpack software.

Diffpack is a set of libraries containing building blocks in numerical methods for PDEs, for example, arrays, linear systems, linear and nonlinear solvers, grids, scalar and vector fields over grids, finite elements, and visualization support. Diffpack utilizes object-oriented programming techniques to a large extent and is coded in the C++ programming language. This means that we must write the PDE solvers in C++. If you do not already know C++, this text will motivate you to pick up the perhaps most popular programming language of the 90s. Most of the knowledge as a C++ programmer can be reused as a C, Java, or even Fortran 90/95 programmer as well, so we speak about a fortunate investment. You do not need to study a textbook on C++ before continuing with the present text, because experience shows that one can get started as a Diffpack programmer without any experience in C++ and learn the language gently as one proceeds with numerical algorithms, the software guide, and more advanced example codes. This book comes with a large number of complete Diffpack solvers for a range of PDE problems, and one can often adapt an existing solver to one's particular problem at hand.

In the past, Diffpack was distributed with a pure software guide as the only documentation. It soon became apparent that successful utilization of numerical software like Diffpack requires (i) a good understanding of the particular formulation of numerical methods that form the theoretical foundations of the package and (ii) a guide to the software tools, exemplified first in detail on simple problems and then extended, with small modifications, to more advanced engineering applications. The software guide can be significantly improved by providing precise references to the suitable generic description of various numerical algorithms. Although this generic view on methods is to a large extent available in the literature, the material is scattered around in textbooks, journal papers, and conference proceedings, mostly written for specialists. It was therefore advantageous to write down the most important numerical topics that must be mastered before one can develop PDE solvers in a programming environment like Diffpack, and tailor the exposition to a programmer.

Our decision to include brief material on the background for and derivation of a model helps the reader with physical knowledge and interest to more clearly see the link between our software-oriented numerical descriptions and comprehensive specialized text on the physics of a problem. Much of the literature on applications, and also on numerical analysis, works with formulations of equations that are not directly suited for numerical implementation. Our chapters on mechanical applications therefore emphasizes a combined physical, mathematical, and numerical framework that aids a flexible software implementation.

The present text has been used in a computationally oriented PDE course at the University of Oslo. Experience from the course shows that whether the aim is to teach numerical methods or software issues, both subjects can benefit greatly from an integrated approach where theory, algorithms, programming, and experimentation are combined. The result is that people use less time to grasp the theory and much less time to produce running code than what we have experienced in the past.

1 | Getting Started | ||

1.1 | The First Diffpack Encounter | ||

1.1.1 | What is Diffpack? | ||

1.1.2 | A Trivial C++ Program | ||

1.1.3 | A Trivial Diffpack Program | ||

1.2 | Overview of
Application Examples | ||

1.1.1 | Very Simple Introductory Program Examples | ||

1.1.2 | Finite Difference Simulators | ||

1.1.3 | Finite Element Simulators | ||

1.1.3 | More Advanced Applications | ||

1.3 | Steady 1D Heat Conduction | ||

1.3.1 | The Physical and Mathematical Model | ||

1.3.2 | A Finite Difference Method | ||

1.3.3 | Implementation in Diffpack | ||

1.3.4 | Dissection of the Program | ||

1.3.5 | Tridiagonal Matrices | ||

1.3.6 | Variable Coefficients | ||

1.3.7 | A Nonlinear Heat Conduction Problem | ||

1.4 | Simulation of Waves | ||

1.4.1 | Modeling Vibrations of a String | ||

1.4.2 | A Finite Difference Method | ||

1.4.3 | Implementation | ||

1.4.4 | Visualizing the Results | ||

1.4.5 | A 2D Wave Equation with Variable Wave Velocity | ||

1.4.6 | A Model for Water Waves | ||

1.5 | Projects | ||

1.5.1 | A Uni-Directional Wave Equation | ||

1.5.2 | Centered Differences for a Boundary-Layer Problem | ||

1.5.3 | Upwind Differences for a Boundary-Layer Problem | ||

1.6 | About Programming with Objects | ||

1.6.1 | Motivation for the Object Concept | ||

1.6.2 | Example: Implementation of a Vector Class in C++ | ||

1.6.3 | Arrays in Diffpack | ||

1.6.4 | Example: Design of an ODE Solver Environment | ||

1.6.5 | Abstractions for Grids and Fields | ||

1.7 | Coding the PDE Simulator as a Class | ||

1.7.1 | Steady 1D Heat Conduction Revisited | ||

1.7.2 | Nonlinear 1D Heat Conduction Revisited | ||

1.7.3 | Empirical Investigation of the Numerical Method | ||

1.7.4 | Simulation of 1D Waves Revisited | ||

1.7.5 | Simulation of 2D Waves Revisited | ||

1.7.6 | Transient Heat Conduction | ||

1.8 | Projects | ||

1.8.1 | Transient Flow Between Moving Plates | ||

1.8.2 | Transient Channel Flow | ||

1.8.3 | Coupled Heat and Fluid Flow | ||

1.8.4 | Difference Schemes for Transport Equations | ||

1.8.5 | 3D Sound Waves |

2 | Introduction to Finite Element Discretization | ||

2.1 | Weighted Residual Methods | ||

2.1.1 | Basic Principles | ||

2.1.2 | Example: A 1D Poisson Equation | ||

2.1.3 | Treatment of Boundary Conditions | ||

2.2 | Time Dependent Problems | ||

2.2.1 | A Wave Equation | ||

2.2.2 | A Heat Equation | ||

2.3 | Finite Elements in One Space Dimension | ||

2.3.1 | Piecewise Polynomials | ||

2.3.2 | Handling of Essential Boundary Conditions | ||

2.3.3 | Direct Computation of the Linear System | ||

2.3.4 | Element-By-Element Formulation | ||

2.3.5 | Extending the Concepts to Quadratic Elements | ||

2.3.6 | Summary of the Element-By-Element Algorithm | ||

2.4 | Example: A 1D Wave Equation | ||

2.4.1 | The Finite Element Equations | ||

2.4.2 | Interpretation of the Discrete Equations | ||

2.4.3 | Accuracy and Stability | ||

2.5 | Naive Implementation | ||

2.6 | Projects | ||

2.6.1 | Steady Heat Conduction with Cooling Law | ||

2.6.2 | Stationary Pipe Flow | ||

2.6.3 | Transient Pipe Flow | ||

2.6.4 | Retardation of a Well-Bore | ||

2.7 | Higher-Dimensional Finite Elements | ||

2.7.1 | The Bilinear Element and Generalizations | ||

2.7.2 | The Linear Triangle | ||

2.7.3 | Example: A 2D Wave Equation | ||

2.7.4 | Other Two-Dimensional Element Types | ||

2.7.5 | Three-Dimensional Elements | ||

2.8 | Calculation of derivatives | ||

2.8.1 | Global Least-Squares Smoothing | ||

2.8.2 | Flux Computations in Heterogeneous Media | ||

2.9 | Convection-Diffusion Equations | ||

2.9.1 | A One-Dimensional Model Problem | ||

2.9.2 | Multi-Dimensional Equations | ||

2.9.3 | Time-Dependent Problems | ||

2.10 | Analysis of the Finite Element Method | ||

2.10.1 | Weak Formulations | ||

2.10.2 | Variational Problems | ||

2.10.3 | Results for Continuous Problems | ||

2.10.4 | Results for Discrete Problems | ||

2.10.5 | A Priori Error Estimates | ||

2.10.6 | Numerical Experiments | ||

2.10.7 | Adaptive Finite Element Methods |

3 | Programming of Finite Element Solvers | ||

3.1 | A Simple Program for the Poisson Equation | ||

3.1.1 | Discretization | ||

3.1.2 | Basic Parts of a Simulator Class | ||

3.2 | Increasing the Flexibility | ||

3.2.1 | A Generalized Model Problem | ||

3.2.2 | Using the Menu System | ||

3.2.3 | Creating the Grid Object | ||

3.3 | Some Visualization Tools | ||

3.3.1 | Storing Fields for Later Visualization | ||

3.3.2 | Filtering Simres Data | ||

3.3.3 | Visualizing Diffpack Data in Plotmtv | ||

3.3.4 | Visualizing Diffpack Data in Gnuplot | ||

3.3.5 | Visualizing Diffpack Data in Matlab | ||

3.3.6 | Visualizing Diffpack Data in Vtk | ||

3.3.7 | Visualizing Diffpack Data in IRIS Explorer | ||

3.3.8 | Plotting Fields Along Lines | ||

3.4 | Some Useful Diffpack Features | ||

3.4.1 | The Menu System | ||

3.4.2 | Multiple Loops | ||

3.4.3 | Computing Numerical Errors | ||

3.4.4 | Functors | ||

3.4.5 | Computing Derivatives of Finite Element Fields | ||

3.4.6 | Specializing Code in Subclass Solvers | ||

3.5 | Introducing More Flexibility | ||

3.5.1 | Setting Boundary Condition Information in the Grid | ||

3.5.2 | Line and Surface Integrals | ||

3.5.3 | Simple Mesh Generation Tools | ||

3.5.4 | Grid Generation by Super Elements | ||

3.5.5 | Debugging | ||

3.5.6 | Automatic Report Generation | ||

3.5.7 | Specializing Code in Subclass Solvers | ||

3.5.8 | Overriding Menu Answers in the Program | ||

3.5.9 | Estimating Convergence Rates | ||

3.5.10 | Axisymmetric Formulations and Cartesian 2D Code | ||

3.5.11 | Summary | ||

3.6 | Step-by-Step Development of a Diffpack Solver | ||

3.6.1 | Physical and Mathematical Problem | ||

3.6.2 | Editing and Writing Source Code | ||

3.6.3 | A Simplified Test Case | ||

3.6.4 | Creating the Grid | ||

3.6.5 | Running Some Initial 2D Simulations | ||

3.6.6 | Running Real Simulations | ||

3.7 | Adaptive Grids | ||

3.7.1 | Grid Classes with Local Mesh Refinements | ||

3.7.2 | How to Extend an Existing Simulator | ||

3.7.3 | Organization of Refinement Criteria | ||

3.7.4 | Grid Refinements as a Preprocessor | ||

3.7.5 | Example: Corner-Flow Singularity | ||

3.7.6 | User-Defined Refinement Criteria | ||

3.7.7 | Transient Problems | ||

3.8 | Projects | ||

3.8.1 | Flow in an Open Inclined Channel | ||

3.8.2 | Stress Concentration due to Geometric Imperfections | ||

3.8.3 | A Poisson Problem with Pure Neumann Conditions | ||

3.8.4 | Lifting Airfoil | ||

3.9 | A Convection-Diffusion Solver | ||

3.10 | A Heat Equation Solver | ||

3.10.1 | Discretization | ||

3.10.2 | Implementation | ||

3.11 | A More Flexible Heat Equation Solver | ||

3.11.1 | About the Model Problem and the Simulator | ||

3.11.2 | Variable Time Step size | ||

3.11.3 | Applying a Transient Solver to a Stationary PDE | ||

3.11.4 | Thermal Conditions During Welding | ||

3.12 | Visualization of Time-Dependent Fields | ||

3.12.1 | Filtering Time-Dependent Simres Data | ||

3.12.2 | Storing Fields at Selected Time Points | ||

3.12.3 | Time Series at Selected Spatial Points | ||

3.12.4 | Using ImageMagick Tools | ||

3.12.5 | Animation Using Plotmtv | ||

3.12.6 | Animation Using Vtk | ||

3.12.7 | Animation Using Matlab | ||

3.12.8 | Real-Time Visualization | ||

3.12.9 | Handling Simulation and Visualization from a Script | ||

3.12.10 | Heat Transfer Exercises | ||

3.13 | A Transient Heat Transfer Application | ||

3.13.1 | The Mathematical and Physical Model | ||

3.13.2 | Implementation | ||

3.13.3 | Testing and Debugging the Initial State | ||

3.13.4 | Creating the Grid | ||

3.13.5 | Running Time-Dependent Simulations | ||

3.13.6 | A Scripting Interface for Automating Simulations | ||

3.14 | Projects | ||

3.14.1 | Transient Heat Transfer in a Two-Material Structure | ||

3.14.2 | Transient Flow with Non-Circular Cross Section | ||

3.14.3 | Transient Groundwater Flow | ||

3.15 | Efficient Solution of the Wave Equation | ||

3.15.1 | Discretization | ||

3.15.2 | Implementation | ||

3.15.3 | Extensions of the Model Problem | ||

3.15.4 | Flexible Representation of Variable Coefficients |

4 | Nonlinear Problems | ||

4.1 | Discretization and Solution of Nonlinear PDEs | ||

4.1.1 | Finite Difference Discretization | ||

4.1.2 | Finite Element Discretization | ||

4.1.3 | The Group Finite Element Method | ||

4.1.4 | Successive Substitutions | ||

4.1.5 | Newton-Raphson's Method | ||

4.1.6 | A Transient Nonlinear Heat Conduction Problem | ||

4.1.7 | Iteration Methods at the PDE Level | ||

4.1.8 | Continuation Methods | ||

4.2 | Software Tools for Nonlinear Finite Element Problems | ||

4.2.1 | A Solver for a Nonlinear Heat Equation | ||

4.2.2 | Extending the Solver | ||

4.3 | Projects | ||

4.3.1 | Operator Splitting for a Reaction-Diffusion Model | ||

4.3.2 | Compressible Potential Flow |

5 | Solid Mechanics Applications | ||

5.1 | Linear Thermo-Elasticity | ||

5.1.1 | The Physical and Mathematical Problem | ||

5.1.2 | A Finite Element Method | ||

5.1.3 | Engineering Finite Element Notation | ||

5.1.4 | Implementation | ||

5.1.5 | Examples | ||

5.1.6 | Elastic Vibrations | ||

5.2 | Elasto-Viscoplasticity | ||

5.2.1 | Basic Physical Features of Elasto-Viscoplasticity | ||

5.2.2 | A Three-Dimensional Elasto-Viscoplastic Model | ||

5.2.3 | Simplification; a Forward Scheme in Time | ||

5.2.4 | Numerical Handling of Yield Criteria | ||

5.2.5 | Implementation | ||

5.2.6 | Examples |

6 | Fluid Mechanics Applications | ||

6.1 | Convection-Diffusion Equations | ||

6.1.1 | The Physical and Mathematical Model | ||

6.1.2 | A Finite Element Method | ||

6.1.3 | Incorporation of Nonlinearities | ||

6.1.4 | Software Tools | ||

6.1.5 | Melting and Solidification | ||

6.2 | Shallow Water Equations | ||

6.2.1 | The Physical and Mathematical Model | ||

6.2.2 | Finite Difference Methods on Staggered Grids | ||

6.2.3 | Implementation | ||

6.2.4 | Nonlinear and Dispersive Terms | ||

6.2.5 | Finite Element Methods | ||

6.3 | An Implicit Finite Element Navier-Stokes Solver | ||

6.3.1 | The Physical and Mathematical Model | ||

6.3.2 | A Finite Element Method | ||

6.3.3 | Solution of the Nonlinear Systems | ||

6.3.4 | Implementation | ||

6.4 | A Classical Finite Difference Navier-Stokes Solver | ||

6.4.1 | Operator Splitting | ||

6.4.2 | Finite Differences on 3D Staggered Grids | ||

6.4.3 | A Multigrid Solver for the Pressure Equation | ||

6.4.4 | Implementation | ||

6.5 | A Fast Finite Element Navier-Stokes Solver | ||

6.5.1 | Operator Splitting and Finite Element Discretization | ||

6.5.2 | An Optimized Implementation | ||

6.6 | Projects | ||

6.6.1 | Analysis of Discrete Shallow Water Waves | ||

6.6.2 | Approximating the Navier-Stokes Equations by a Laplace Equation |

7 | Coupled Problems | ||

7.1 | Fluid-Structure Interaction; Squeeze-Film Damping | ||

7.1.1 | The Physical and Mathematical Model | ||

7.1.2 | Numerical Methods | ||

7.1.3 | Implementation | ||

7.2 | Fluid Flow and Heat Conduction in Pipes | ||

7.2.1 | The Physical and Mathematical Model | ||

7.2.2 | Numerical Methods | ||

7.2.3 | Implementation | ||

7.3 | Projects | ||

7.3.1 | Transient Spherical-Symmetric Thermo-Elasticity | ||

7.3.2 | Transient 2D/3D Thermo-Elasticity | ||

7.3.3 | Convective-Diffusive Transport in Viscous Flow | ||

7.3.4 | Chemically Reacting Fluid |

A | Mathematical Topics | ||

A.1 | Scaling and Dimensionless Variables | ||

A.2 | Indicial Notation | ||

A.3 | Compact Notation for Difference Equations | ||

A.4 | Stability and Accuracy of Difference Approximations | ||

A.4.1 | Typical Solutions of Simple Prototype PDEs | ||

A.4.2 | Physical Significance of Parameters in the Solution | ||

A.4.3 | Analytical Dispersion Relations | ||

A.4.4 | Solution of Discrete Equations | ||

A.4.5 | Numerical Dispersion Relations | ||

A.4.6 | Convergence | ||

A.4.7 | Stability | ||

A.4.8 | Accuracy | ||

A.4.9 | Truncation Error | ||

A.4.10 | Traditional von Neumann Stability Analysis | ||

A.4.11 | Example: Analysis of the Heat Equation | ||

A.5 | Exploring the Nature of Some PDEs | ||

A.5.1 | A Hyperbolic Equation | ||

A.5.2 | An Elliptic Equation | ||

A.5.3 | A Parabolic Equation | ||

A.5.4 | The Laplace Equation Solved by a Wave Simulator | ||

A.5.5 | Well Posed Problems |

B | Diffpack Topics | ||

B.1 | Brief Overview of Important Diffpack Classes | ||

B.2 | Diffpack-Related Operating System Interaction | ||

B.2.1 | Unix | ||

B.2.2 | Windows | ||

B.3 | Combining Diffpack with Other Types of Software | ||

B.3.1 | Calling Other Software Packages from Diffpack | ||

B.3.2 | Calling Diffpack from Other Types of Software | ||

B.4 | Basic Diffpack Features | ||

B.4.1 | Diffpack man pages | ||

B.4.2 | Standard Command-Line Options | ||

B.4.3 | Generalized Input and Output | ||

B.4.4 | Automatic Verification of a Code | ||

B.5 | Visualization Support | ||

B.5.1 | Curves | ||

B.5.2 | Scalar and Vector Fields | ||

B.6 | Details on Finite Element Programming | ||

B.6.1 | Basic Functions for Finite Element Assembly | ||

B.6.2 | Using Functors for the Integrands | ||

B.6.3 | Integrating Quantities over the Grid or the Boundary | ||

B.6.4 | Class Relations in the Finite Element Engine | ||

B.7 | Optimizing Diffpack Codes | ||

B.7.1 | Avoiding Repeated Matrix Factorizations | ||

B.7.2 | Optimizing the Assembly Process | ||

B.7.3 | Optimizing Array Indexing |

C | Iterative Methods for Sparse Linear Systems | ||

C.1 | Classical Iterative Methods | ||

C.1.1 | A General Framework | ||

C.1.2 | Jacobi, Gauss-Seidel, SOR, and SSOR Iteration | ||

C.2 | Conjugate Gradient-Like Iterative Methods | ||

C.2.1 | Galerkin and Least-Squares Methods | ||

C.2.2 | Summary of the Algorithms | ||

C.2.3 | A Framework Based on the Error | ||

C.3 | Preconditioning | ||

C.3.1 | Motivation and Basic Principles | ||

C.3.2 | Classical Iterative Methods as Preconditioners | ||

C.3.3 | Incomplete Factorization Preconditioners | ||

C.4 | Multigrid and Domain Decomposition Methods | ||

C.4.1 | Domain Decomposition | ||

C.4.2 | Multigrid Methods |

D | Software Tools for Solving Linear Systems | ||

D.1 | Storing and Initializing Linear Systems | ||

D.1.1 | Vector and Matrix Formats | ||

D.1.2 | Detailed Matrix Examples | ||

D.1.3 | Representation of Linear Systems | ||

D.2 | Programming with Linear Solvers | ||

D.2.1 | Gaussian Elimination | ||

D.2.2 | A Simple Demo Program | ||

D.2.3 | A 3D Poisson Equation Solver | ||

D.3 | Classical Iterative Methods | ||

D.4 | Conjugate Gradient-like Methods | ||

D.4.1 | Symmetric Systems | ||

D.4.2 | Nonsymmetric Systems | ||

D.5 | Preconditioning Strategies | ||

D.6 | Convergence History and Stopping Criteria | ||

D.7 | Example: Implicit Methods for Transient Diffusion | ||

D.8 | High-Level Stencil Programming of Finite Difference Schemes | ||

D.8.1 | Finite Difference Stencils | ||

D.8.2 | Basic Structure of a Stencil-Based Simulator | ||

D.8.3 | Defining the Stencils |

Copyright © inuTech GmbH. All rights reserved. Questions and comments regarding this web site should be directed to diffpack@inutech.de

contact | privacy policy | site notice

contact | privacy policy | site notice