Introduction to Modern Scientific Programming and Numerical Methods
It was essentially these skills that I tried to expose students to in my ASTE-499 Applied Scientific Computing course at USC. That course was taught mainly off slides, but in the future, students taking this course – and similar courses at different universities – will have a textbook at their disposal. We recently got a contract from CRC Press (the publisher of my Plasma Simulations by Example text) for Introduction to Modern Scientific Programming and Numerical Methods, due in July 2021! This book is co-authored with Prof. Joseph Wang at USC, and Dr. Robert Martin at the Air Force Research Laboratory (AFRL). I put this page together to solicit feedback. Below you will find the proposed table of contents. Some material will likely end up getting re-ordered but it summarizes the topics we plan to include. Please leave a comment below (or just email me) with any suggestions. Also let us know if you would like to review draft copies.
Table of Contents
1. Introduction to Scientific Computing
This chapter introduces basic concepts from scientific computing by demonstrating how to simulate a ball dropped from some height and falling under the action of gravity (and drag) and bouncing on surface impact (example differential equation).
1.1 Problem Description
We start by describing the problem. Here we introduce limits, Taylor series, differentiation, and numerical round off errors.
1.2 Pseudo code
Next, a pseudo code is developed. Basic programming concepts, such as variables, data types, for loops, if statements, functions, and input/output are introduced. We also discuss binary representation and finite-precision computer arithmetic.
1.3 Programming Languages
1.4 Integration Methods
Numerical results are compared to theoretical model. We also discuss visualization options, including use of Python/Matlab plotting and the use of Tecplot and Paraview. We discuss convergence, consistency, and alternate integration methods including forward and backward Euler, Leapfrog, Crank-Nicolson, multistep, Runge Kutta. These examples are developed primarily in Matlab and Python.
1.5 Additional Physics
Additional physics is included to increase complexity. We mainly focus on adding aerodynamic drag and surface reflection. Root finding is introduced.
1.6 Modularity and Data Reduction
The simulation is expanded to include multiple bouncing balls. This section allows us to introduce arrays, functions, and function arguments. We also discuss pass-by-reference and pass-by-value. Additional visualization options such as scatter plots and animations are discussed.
2. Numerical Analysis and Linear Algebra
In this chapter, we cover many elementary topics from Numerical Analysis, including equation solving and linear algebra.
2.1 Root Finding
We start by discussing root finding (solving equations of one variable). We cover Bisection Method, Fixed point iteration, Newton’s method. Error analysis and rate of convergence is also discussed.
2.2 Differentiation and Integration
We discuss Taylor Series in more detail and discuss “big O” notation. We then cover numerical integration (quadrature) including Trapezoidal and Simpson rule, and Gaussian Quadrature. Double (or triple) integrals are also covered.
2.3 Matrix Algebra
We switch to Linear Algebra, and introduce the concept of a matrix and a vector. We learn about matrix multiplication, and see how to implement it in Matlab and also languages where it is not natively supported (such as C++). Einstein (index) notation is also discussed. We also learn about different matrix types, such as banded, sparse, singular, symmetric, diagonally dominant, and positive definite. We demonstrate how to use rotation and translation matrix for coordinate transformation.
2.4 Solving Systems of Equations
We introduce matrix inverse for solving systems of equations. We also cover topics such as determinants and cofactors and see how to compute inverse of small matrices.
2.5 Direct Methods
We cover Gaussian elimination, a direct solution method for larger systems for which computing a matrix inverse is not practical. We discuss the Thomas algorithm for tridiagonal systems. These examples are worked out in Matlab, Python and C++.
2.6 Iterative Methods
We discuss iterative methods for solving matrix problems, including Jacobi, and Gauss-Seidel, and Successive Over Relaxation. We learn about residue norm and convergence.
2.7 Newton-Raphson Method
Non-linear systems of equations are introduced and we discuss the popular Newton-Raphson method
2.8 Conjugate Gradient
The preconditioned conjugate gradient method is discussed and convergence is compared to Jacobi
The multigrid method is covered next. We cover high and low frequency errors.
2.10 LU and Cholesky Decomposition
LU and Cholesky matrix decomposition is discussed and illustrated with Matlab and C++ example
Eigenvalues and eigenvectors are covered. Some popular algorithms such as power method and QR decomposition are illustrated.
2.12 Numerical Libraries
We review popular numerical libraries such as NumPy/SciPy and BLAS/LAPACK and see how to use them in our codes
3. Interpolation and Signal Processing
This chapter focus on methods for fitting curves to large amounts of data and reducing data trends
3.1 Polynomial Fits
We cover Lagrange polynomials, Hermite interpolation, Bezier splines, Richardson extrapolation
3.2 Least Squares Fit
The Least Squares method is covered next
3.3 Fourier Series
Fourier Series and the FFT method is covered. We demonstrate the method using built-in functions from Matlab and Python and with a custom C++ code.
Methods for filtering data are covered.
3.5 Probability Distributions
Probability distribution functions are covered. We learn about cumulative distribution function and stochastic (random) sampling. Earth Mover’s method is discussed for comparing distribution functions.
4. Object Oriented Programming
This chapter we introduce advanced C++ syntax and learn about object oriented programming
4.1 Data Encapsulation
Here we discuss benefits of object oriented programming, including data encapsulation, access control, storage, and operator overloading without getting into the actual implementation details
4.2 Structures and Classes
Structs and classes are introduced here. They form the foundation of object oriented programming. We learn about access control, constructors, and destructors.
4.3 Pointers, References, and Dynamic Memory Allocation
We then start with a crash course in C++. We start by covering pointers, references, and dynamic memory allocation. These topics, while not difficult, are are a source of confusion to many new students of C++.
4.4 Operator Overloading
We learn about implementing operators for manipulating our custom classes. These operators allow us to write the code in a more concise, mathematical way.
We see how templates allow us to generalize the functionality of a class to operate on different data types. We implement a “vec3” data object for storing 3D vectors.
4.6 Storage Containers
Next we discuss various storage containers and see how they are implemented inthe C++11 standard library. These include vector, linked list, map,and set. We also learn about octrees.
4.7 Lambda Functions and Algorithms
We discuss anonymous lambda functions and see how they can be used with standard library algorithms to perform manipulations such as data sorting.
4.8 Putting it Together
We close the chapter by utilizing the new skills to write a multiple-bouncy ball simulation from Chap. 1 as well as a matrix solver. We compare the performance to Python/MATLAB version.
5. Interactive Web Applications
We start by describing the basic syntax of HTML files. We learn about basic element types such as <head>, <body>, <p>, <h?>, <a>, <ul>, <li>, <input>, and <div>.
We next learn how to change the visualization representation using styles. We learn about the rule selectors in CSS files.
We then introduce the <canvas> element that can be used for 2D and 3D rendering. We mainly focus on the 2D “context”. We learn how to draw shaded circles, lines, and rectangles. Discussion of webGL (for parallel 3D rendering) is left for Chapter 10.
We use the above material to develop several interactive simulations covering methods discussed previously. We start with a simulation that adds a new bouncing ball at mouse click location. We then write code for filtering and interpolating .csv data that is drag-and-dropped from the file system.
6. Documentation and Code Testing
This chapter covers additional programming topics which are not commonly discussed but knowledge of which is critical for developing large-scale codes that can be utilized and expanded by other researchers.
6.1 Testing and Visualization
Here we start by introducing code validation and verification, and sensitivity and uncertainty analysis. We also discuss various methods for visualizing data, including isosurfaces, scatter plots, and virtual probes
Techniques for code debugging are discussed. We also introduce the use of automated tools such as valgrind for finding memory leaks.
6.3 Coding Standards
We discuss the need for coding standards for code maintainability. Several popular standards are introduced.
6.4 Test Suites
GTest and JUnit is discussed here. We learn about benefits of automated unit testing
6.5 Version Control
We introduce version control using CVS,SVN, Mercurial, but mainly focusing on Git (current favorite).
6.6 Documentation Systems
In-line documentation using Doxygen is discussed. We include a brief description of LaTeX, focusing on mathematical formulas. We see how to include LaTeX syntax in Doxygen.
7. Partial Differential Equations
This chapter introduces basic solution approaches for selected partial differential equations. We illustrate the basic methods in Matlab and Python, but the more advanced methods are illustrated mainly using C++.
We learn about elliptic, parabolic, and hyperbollic partial differential equations. We discuss their derivation, and learn about initial and boundary conditions
7.2 Finite Difference Method
The finite different method is discussed. We learn about discretization, mesh nodes and cells, and matrix representation. First and second order differencing is discussed.
7.3 Heat Equation
7.4 Finite Volume Method
The finite volume method is covered. We learn about the Divergence and Stokes Theorem, and also see how it naturally retrieves coefficients in cylindrical coordinates. We see that for uniform Cartesian meshes, the method reduces to the Finite Difference. A FVM example is developed in C++.
7.5 Finite Element Method
The FEM is discussed briefly. We use the method to develop a simple FEM solver of the heat equation in C++.
7.6 Wave Equation
We then discuss the parabolic wave equation and cover several popular solution methods.
7.7 Beam Equation
The beam equation (example of a higher-order PDE) is covered and several popular methods are illustrated.
7.8 Advection-Diffusion Equation and the SIMPLE Method
Another model equation from aerodynamics is discussed. We discuss some approaches for the Advective Term starting with the Upwind method. We then cover the SIMPLE method used for incompressible fluids.
7.9 Vlasov Equation
Vlasov Equation covering the evolution of a velocity distribution function is introduced. We cover a simple first order integration scheme.
8. Particle Methods
This chapter introduces numerical methods based on Lagrangian Mechanics and stochastic (particle) methods.
8.1 Random Numbers and Monte Carlo
We learn about using random numbers to sample values from Maxwellian velocity distribution function
8.2 Gas Kinetics
We modify the “bouncy-balls” example from Ch. 4 to simulate free molecular flow gas dynamics. We learn about using a computational mesh to sample velocity moments for computing density, stream velocity, and temperature.
8.3 Direct Simulation Monte Carlo
The Direct Simulation Monte Carlo (DSMC) method is introduced next to add intermolecular collisions
8.4 Particle in Cell
We cover Poisson’s equation (which has the same form as steady heat equation covered in Ch. 7) and electrostatics. We see how to use the PIC method to simulate ionized gas (plasma) dynamics.
8.5 Molecular Dynamics
We discuss the molecular dynamics (MD) method that is used in computational chemistry to model molecular structures
8.6 Smoothed Particle Hydrodynamics
In this section we describe the SPH method that blends Eulerian and Lagrangian description of fluid motion
9. Parallel Processing
This chapter introduces the reader to code optimization and high performance computing. These examples are worked out in C++ and CUDA
9.1 Profiling and Memory Access
We learn how to use timers and profilers to find parts of the code that run slowly. We also demonstrate how re-arranging memory access can lead to drastic performance gain, and learn about CPU cache.
We discuss vector operations and see how to split them up using threads. We cover race condition, atomics, and locks (and how to avoid needing them). We see how some library functions are not thread safe and lead to serialization and slow performance (i.e. Random). We also demonstrate that there is finite penalty for thread creation that may make multithreading not be efficient for small computations. Parallel efficiency is compared.
9.3 Message Passing Interface
We learn about MPI and discuss several popular matrix solution methods such as red-black ordering and domain decomposition. We also discuss approaches for particle simulations. We briefly describe setting up a cluster using Amazon Web Services and provide a short list of basic Linux commands.
Graphics card processing using CUDA is discussed next. We cover memory transfer, pinned memory, and streams.
10. Optimization and Machine Learning
This chapter introduces optimization methods for large problems and machine learning
10.1 Optimization Basics
We discuss the need for optimization – often we do not have the full parameter space for the simulation and need to perform parameter space search. We cover the brute force and shooting methods. We see this approach is inefficient for a large parameter space.
10.2 Gradient Descent
The gradient descent method is described next. This follows up on the Conjugate Gradient discussion for matrix solving.
10.3 Genetic Algorithms
Genetic algorithms are described next. We develop simulation to find the shortest path between points.
10.4 Neural Networks
Deep neural networks, which form the foundation of machine learning, are covered. We cover basic activation functions such as ReLU and logistic curve. Cost function is also covered.
We next see how back-propagation is used to update activation function weights
11. Embedded Systems
This chapter discusses programming microcontrollers and FPGAs, as they allow us to develop data acquisition sensors for code validation.
11.1 Introduction to Arduino
We learn about the basics of the popular Arduino microcontroller platform, which is programmed using a simplified form of C++.
11.2 Sensors and Electronics
Here we cover basic electronic components including resistors, capacitors, transistors, diodes, and LEDs. We also discuss breadboards and soldering.
11.3 Basic DAQ system
Putting it all together, we learn how to create a simple data acquisition system that collects pressure and humidity data from multiple sensors and stores the data onto an SD card
11.4 Edge Computing
We learn how to use TensorFlow Lite to add machine learning to microcontrollers to enable “edge computing”
11.5 FPGAs and Verilog
Here we cover what FPGAs are and introduce the syntax of Verilog. We learn how to write our own hardware instructions for performing rapid math operations. We write code for rapidly “mux”-ing data passed by an Arduino.
Appendix: Programming Syntax
Various code examples are posted here:
Initial Progress on Micro Cathodic Arc Thruster Plume Simulations
Implementation of VTK-based 3D visualization capability in a solver GUI
This looks great. As a former CS guy in aero, this is pretty much everything I would hope to see. A few random suggestions: validation cases for testing consistency, versioning (i.e. bake in to your results what build produced the result), C++11 style pointers, build systems, file formats like json & xml for configuration (please do not invent new text file formats for your solver, it makes me cry…). I would say also that some discussion of UI options wouldn’t be amiss. PyQT, etc. Maybe interprocess comm too. There’s a lot of merit to doing solvers in one language and your problem setup and viz in another (C++ and python being a common example). I also wouldn’t mind seeing some EM-specific examples, and maybe some inside baseball stuff like using ghost cells, specific treatment of BCs, etc to help improve results. Though this may be in your other book!
Anyhow these are all relatively minor points. Looking forward to this one!
Thanks Philip for many great suggestions! I plan to have a description of various common file formats in the text, somewhere. However, utilizing either is not trivial in C++, at least not without using external libraries. At least Java has an XML parser (as used by my 2D code Starfish), but I ended up writing my own JSON-like parser for CTSP.
Here are my wild thoughts – no need to respond. Good to hear you are doing well.
Looks like you didn’t leave anything out. How many thousands of pages are you planning on? I think you need two books. One for programming principles for engineers/scientists that includes all the computer stuff and basic numerical analysis techniques that are the most commonly used in engineering. The other is numerical analysis book with all the other numerical techniques you mention using a single programming language. The second book is similar to a few already on the market.
I’m not sure where the chapter on embedded systems fits.
You may want to spend some time on “make” and tool chains. Coming from matlab linking will be a foreign concept.
I’m not sure about the logic of putting optimization in the same chapter as machine learning. I think you are talking about regression machine learning which fits more with the section on interpolation. Also, I’m not sure how interpolation fits with signal processing.
You may want to add an example of using a large off-the-shelf code like Sparta / Paraview?
Is the audience for the book people who need to understand what is going on inside of a code like Sparta or people who can build and augment large complex codes for a particular domain?
Thanks Ken. So we actually went back and forth on whether to do a single book or a multi-volume, and we settled on a single volume as otherwise, there is not much difference from having multiple books on related subjects. The total number of pages is going to be only about 450 (so about 40 pages per chapter). I realize this may sound too little, but the point of this text is to introduce the reader to these various topics (“I didn’t even know what I didn’t know”) and refer them to other texts for the actual details. I do plan on having a brief section on data visualization with Paraview. For your two other questions, I guess I don’t see the topics to be that incompatible. Interpolation (or perhaps what I meant is smoothing/fitting) involves looking for trends in noisy data which is essentially signal processing, no? And with optimization and machine learning, the idea is that we are looking for some parameter – or a set of relationships – that will get us a desired answer given some set of inputs, with us letting the computer find it on its own. But perhaps there is a deeper philosophical distinction between the two I am not aware of.
When will the book be out for sale? Or it is and I missed reading that somewhere
Not until late next year. It is due to the publisher in July 2021, so likely targeting November/December release.
Hi, Chemical Engineer self taught programmer here. My opinion is similar to the one given by Ken Daniel. My biggest concern is that you try to cover too many topics, that could leave you in some cases with a single page to “cover” a topic.
But everyone will clear out once the printed version is ready, I hope to see how you manage to solve this problem.
Thank you Francisco. Yes, this is definitely a concern. This is why I am not trying the book to be fully inclusive of all methods/etc, but to serve as the proper introductory material, with references to more detailed texts. For example, there are 500+ page books covering CFD, and here we will try to just cover basics of Eulerian methods in about 40 pages.
I feel like Newton-Raphson is a bit out of order… If I were to structure chapter 2, I might have something like:
a) Solving Linear Systems of equations
1) Direct methods
2) Stationary methods
3) Steepest descent and conjugate gradient methods
4) Krylov methods
5) Geometric MultiGrid
6) Algebraic Multigrid
b) Solving Nonlinear Systems of Equations
3) Nonlinear conjugate gradient (I know you have this in optimization, but could fit here as well).
4) (others as relevant)
Also, some other things to consider might be automatic differentiation (aka algorithmic differentiation), complex-step, symbolic calculations, scientific file formats like HDF5 & NetCDF, Verification & Validation (at least what they are, and how they’re different).
Perhaps also some discussion on how to present results / data, e.g. when to use log scales, when to use a diverging colormap vs when to use a “single-sided” colormap (e.g. white-to-red), when to use a contour plot vs a surface plot, etc.
Would be good to present Amdahl’s Law in the parallel section.
Thanks Greg, these are all good suggestions. But I think that Krylov methods are much more complicated to implement than N-R, no? For instance, I have implemented PCG many times, but every time it was just following a recipe in a numerical book. I don’t actually know how to derive the solver (yet?), but N-R is fairly straightforward.
I think for the most part, the majority of users of scientific computing will be calling a library that already exists for the more complicated methods. I think what’s more important for these more involved solvers is to give an introduction to *how* they work, *when* to use them, and *why* they fail. To your earlier point, it seems like this book is a bit of an introductory book – with references to the vast 500+ page treatises on the subject. Thus I think it’s more important to give the reader a full overview of linear solvers THEN give a full overview of nonlinear solvers, rather than jump between linear -> nonlinear -> linear -> nonlinear (optimization). Especially since Newton methods are effectively just wrappers around linear solvers.
Might also be a good idea to talk preconditioning…
Sounds good. I will send you the draft copy of this chapter once ready, if interested, to get additional feedback.
I’d be happy to review a draft of the chapter. Feel free to direct email me (I presume you have my email) if you’d like to.
Please mention semi-discretisation schemes for PDEs(Heat PDE, Wave PDE) along with the state Laplacian matrix generation for Cartesian, semi-Cartesian and polar coordinates evaluated with examples involving various initial and boundary conditions.