Suite-CFD: An Array of Fluid Solvers Written in MATLAB and Python

Computational Fluid Dynamics (CFD) models are being rapidly integrated into applications across all sciences and engineering. CFD harnesses the power of computers to solve the equations of fluid dynamics, which otherwise cannot be solved analytically except for very particular cases. Numerical solutions can be interpreted through traditional quantitative techniques as well as visually through qualitative snapshots of the flow data. As pictures are worth a thousand words, in many cases such visualizations are invaluable for understanding the fluid system. Unfortunately, vast mathematical knowledge is required to develop one’s own CFD software and commercial software options are expensive and thereby may be inaccessible to many potential practitioners. To that extent, CFD materials specifically designed for undergraduate education are limited. Here we provide an open-source repository, which contains numerous popular fluid solvers in 2 D (projection, spectral, and Lattice Boltzmann), with full implementations in both MATLAB and Python3. All output data is saved in the . v t k format, which can be visualized (and analyzed) with open-source visualization tools, such as VisIt or ParaView. Beyond the code, we also provide teaching resources, such as tutorials, flow snapshots, measurements, videos, and slides to streamline use of the software.


Introduction
Computational Fluid Dynamics (CFD) models are being applied to problems across all sciences and engineering. From designing more aerodynamic sportswear and vehicles, to understanding animal locomotion, to personalized medicine, to disease transmission, and to predicting hurricanes and atmospheric phenomena, it is difficult to find situations where greater knowledge of the underlying fluid dynamics isn't desired or valuable. Due to its immense importance, many scientists have dedicated their careers to the field, whether they study particular fluid phenomena or develop tools, either experimental or numerical, for other scientists and engineers to use. Possibly because of how vital understanding fluid dynamics is to human flourishing, proving existence, uniqueness, and the possible breakdown of solutions to the system of non-linear partial differential equations that govern fluid dynamics is one of the Millennium Problems, to which the Clay Institute offers a 1 million dollar prize for solving [1]. On the other hand, CFD allows us to solve these equations, the Navier-Stokes equations, which are the equations that detail the conversation of momentum and mass for a fluid. For an incompressible, viscous fluid, they can be written as follows, ρ ∂u ∂t (x, t) + (u(x, t) · ∇) u(x, t) = −∇p(x, t) + µ∆u(x, t) Rather than develop a collection of modules that teach how to implement particular CFD numerical schemes, we offer the scientific community a variety of popular fluid solvers that solve traditional problems in fluid dynamics, with two independent but equal implementations written in MATLAB and Python. The emphasis is not placed on students having to develop or implement numerical schemes, but to allow them the opportunity to get an accelerated start in CFD, and run, tweak, and analyze simulations of traditionally studied problems in fluids courses, as well as explore data visualization and how to present data that is appropriate for fluid physics interpretation. All of this can be done while testing their own hypotheses in familiar programming environments. In the remainder of this manuscript, we will give a high-level overview of the CFD schemes currently implemented in the software (Section 2), guided tutorials on how to run, visualize, and analyze simulations (Section 3), and provide numerous examples to which students may elect to study or modify as well (Section 4).
Furthermore, we provide an in-depth mathematical description of each numerical scheme in Appendix B. The open-source software repository can be accessed at https://github.com/nickabattista/ Holy_Grail. Any simulation data produced from the software can be visualized and analyzed in open-source software, such as VisIt [40], maintained by Lawrence Livermore National Laboratory, or ParaView [41], developed by Kitware, Inc., Los Alamos National Laboratory, and Sandia National Labs. We also provide teaching resources, such as slides, figures, and movies in the Supplemental Materials.

Brief Overview of the Three Fluid Solvers
In this section we will give a brief overview of the three fluid solvers currently implemented in the software -a projection method, a spectral (FFT) method, and the Lattice-Boltzmann method. Our goal is to provide students a high-level overview of the methods and how they compare to one another. Further details regarding their mathematical foundations and implementations are given in Appendix B.

Projection Method
The projection method was first introduced by Chorin in 1967 [2] and independently by Temam in 1968 [4], to solve the viscous, incompressible Navier-Stokes equations. Projection methods are finite difference based numerical schemes [42], in which the fluid equations are solved in an Eulerian framework, i.e., the computational grid is static and fluid dynamics variables, like velocity, pressure, etc., are measured at particular locations on the grid over time during a simulation. Projection methods have been extensively used throughout the fluid dynamics community for decades, in numerous applications across many fields, while being continually improved for efficiency and accuracy [5,[43][44][45].
In a nutshell, this method decouples the velocity and pressure fields, using operator splitting and a Helmholtz-Hodge decomposition. This makes it possible to explicitly solve Equations (1) and (2) in a few steps [2] while also increasing computational efficiency. Please see Appendix B.1 for further mathematical details.
Furthermore, projection methods allow one to define explicit boundary conditions (BCs) for velocity on edges of the computational domain, whether they are Dirichlet, Neumann, or Robin boundary conditions [5,46]. In our software, we allow users to impose tangential velocity boundary conditions along the edges of the rectangular domain, while the normal direction boundary conditions for velocity are explicitly set to zero. The choice to set normal BCs to zero was made in order to simplify the code by ensuring that volume conservation would be automatically satisfied, i.e., situations will not arise in which a specified amount of fluid is being pumped into the domain while a different amount is leaving the domain per unit time.
Having defined the fluid velocity as u = (u, v), one can set the boundary conditions for u on the top and bottom of the domain (v = 0 along these edges) or v on the left and right sides (u = 0 along these edges). See Figure 1 for an illustration. Later in Section 4, we will showcase two different traditional problems in fluid dynamics with a projection method:

Cavity Flow 2. Circular Flow in a Square Domain
The examples are different due to the changes in boundary conditions that the user can impose, see Figure 1. In Figure 1 we illustrate the boundary conditions for each example: cavity flow (a) and circular flow (b). This figure also illustrates that the domain could be constructed to be either rectangular or square and that the boundary conditions do not have to be uniform among all sides, nor across the domain. We will also compare different fluid scales in both examples, using the Reynolds Number, Re. The Reynolds Number is given by where ρ and µ are the fluid's density and dynamic viscosity, respectively, while L and V are a characteristic length-and velocity-scale of the system. In essence, Re is a ratio of inertial forces to viscous forces in a system. If Re >> 1, inertia dominates, if Re << 1, viscous forces dominate, and if Re = 1, the inertia and viscous forces are balanced. An example of high Re flow would be the flow around a marble quickly falling through air. An example of low Re flow would be the fluid flow around a marble slowly falling through honey.

Spectral Method (FFT)
Spectral methods are a different class of differential equation solvers. They are well-known for achieving highly accurate solutions, i.e., spectral accuracy [47][48][49]. Spectral accuracy occurs when error decreases exponentially with a small change in grid resolution, e.g., there is a linear relationship between the logarithm of error and grid resolution, i.e., where N is the number of grid points. These methods approximate solutions using orthogonal expansions, such as Fourier Series or other orthogonal basis of functions, like Chebyshev Polynomials or Legendre Polynomials. In contrast to finite difference based schemes, like the projection method, where solutions are approximated at specific locations in space, in spectral methods a particular orthogonal series expansion is chosen and the earnest is on finding the coefficients of the expansion.
Here we chose to use the basis of Discrete Fourier Transform (DFT). The DFT can be used to express functions that are not periodic, unlike Fourier Series expansions, which are used to represent periodic functions. Moreover, we used the Fast Fourier Transform (FFT), which yields the same numerical values as a DFT, but is considerable more computationally efficient, i.e., fast.
Thus, we implemented a FFT based spectral scheme to solve the viscous, incompressible Navier-Stokes equations in 2D. We also chose to solve these equations in their vorticity formulation (see Appendix B.2 for mathematical details). The vorticity formulation will lead to a few numerical benefits.
In a nutshell, this will recast our problem into solving for the vorticity, ω ω ω = ωk, and the streamfunction, ψ. Velocity can be defined in terms of the curl of the streamfunction, i.e., a vector potential, see below Hence this allows us to find the velocity field, u = (u, v), once we find ψ, e.g., Moreover, the incompressibility condition (Equation (2)) will be automatically satisfied since u is defined to be the curl of a scalar, and the divergence of the curl of a vector is always identically zero. Furthermore, working in the vorticity formulation requires us to only solve a Poisson equation for the streamfunction in terms of the vorticity, ω ω ω, i.e., ∆ψ = −ω. (5) Notice that Equation (5) is a linear equation. Moreover, upon taking the FFT of the Equation (5), the computation is transformed into frequency (Fourier) space, which offers two advantages-increased speed and accuracy.
In practice, we will solve for the streamfunction at the next time-step based on the previously computed vorticity and then update the vorticity using information from the newly updated streamfunction. We chose to use a Crank-Nicholson time-stepping scheme [50] to update the vorticity to the next time-step. The Crank-Nicholson scheme is well-known for being unconditionally-stable for diffusion problems and as well as for being 2 nd order accurate in time and space, as an implicit method [51]. These accuracy and stability properties make the Crank-Nicholson an ideal candidate for numerically solving such equations.
Although using a FFT (spectral method) preserves high accuracy of the spatial information at a particular time-step, numerical error still unfortunately creeps in due to the time-stepping nature of solving an evolution equation for the fluid's momentum (see Equation (A15) in Appendix B.2). That is, evolving the system forward in time is where the brunt of the error takes place, rather than the spatial discretization.
In contrary to the projection method, this FFT-based approach allows us to explore periodic boundary conditions. For example, if fluid flow is moving vertically-upwards through the top of the computational domain, it will re-enter moving vertically-upwards through the bottom of the domain. This has advantages and disadvantages. Some advantages are that explicit boundary conditions do not need to be satisfied nor enforced, which aids in computational speed-up. However, to that extent, one disadvantage is that this makes it more difficult to enforce particular boundary conditions when desired. FFT-based fluid solvers have been used in the fluid-structure interaction community, in particular within immersed boundary methods [32,36,52], when one wants to study the interactions of an object and the fluid to which it's immersed. As the focus is on the object and fluid interactions, it is typically modeled away from the domain boundaries in the middle of the computational domain to avoid boundary artifacts. For a mathematically detailed description of this FFT-based spectral method see Appendix B.2.
For this particular FFT-based fluid solver for the vorticity-formulation of the viscous, incompressible Navier-Stokes equations, we will highlight the following examples:
Evolution of Vorticity from an Initial Velocity Field Contrary to the examples shown for the projection method in Section 2.1, these examples are not differentiated by different boundary conditions, but rather different initial vorticity configurations in the computational domain. Figure 2 gives the initial vorticity configurations for the two cases considered: side-by-side vorticity region interactions (left) and an initial vorticity field that defined by a velocity vector field (right). In Figure 2a, counterclockwise (CCW) and clockwise (CW) correspond to regions of uniform vorticity, where vorticity initialized as a positive or negative constant for CCW and CW, respectively. In all other regions of the computational domain, the vorticity is either initialized to zero (Figure 2a). Although not shown, one could also initialize simulations on a rectangular grid. In Figure 2b, the initial vorticity is defined by computing the vorticity from a simulation's velocity vector field. The velocity vector field was taken from a time-point in a Pulsing_Heart simulation [53] in the open-source IB2d software [32,36]. In this example, we illustrate that if a snapshot of the velocity field is known, it is possible to evolve the fluid dynamics forward using this spectral (FFT) method, even though it is based on the vorticity formulation of the Navier-Stokes equations.

Lattice Boltzmann Method
The Lattice Boltzmann method (LBM) does not explicitly (or implicitly) solve the incompressible, Navier-Stokes equations, rather it uses discrete Boltzmann equations to model the flow of a fluid [11]. In a nutshell, the LBM tracks fictitious particles of fluid flow, thinking of the problem more as a transport equation, e.g., where f (x, t) is the particle distribution function, i.e., a probability density, u is the fluid particle's velocity, and Ω is what is called the collision operator. However, rather than these particles moving in a Lagrangian framework, the Lattice Boltzmann method simplifies this assumption and restricts the particles to the nodes on a lattice. In traditional CFD methods, one typically discretizes the fluid equations in an Eulerian framework, as in a projection or spectral method. When you are free of such restriction, as in the LBM, it allows one to more readily deal with complex boundaries, model porous structures and multi-phase flow, as well as incorporate microscopic interactions, while also allowing for massive parallelization of the algorithm, leading to increased computational efficiency [12,[54][55][56].
The LBM involves collision and streaming steps, where fictional fluid particles consecutively propagate (and collide) over a discretized lattice, and the fluid density is evolved forward. There are multiple ways to handle boundary conditions in the LBM. One such way is to use bounce-back boundary conditions, which are executed by masking points in the domain via boolean expressions, thus defining regions where the fluid can and cannot flow [55]. In short, the propagating (streaming) directions are simply reversed when they hit a boundary node, see Figure A3 in Appendix B.3. The bounce-back conditions can be used to enforce the necessary and physical no-slip boundary conditions for fluid flow as well as defining solid objects in the interior of the grid (for an example see Figure 33). For a detailed mathematical description of the LBM see Appendix B.3.
To highlight the LBM fluid solver, we showcase the following examples: 1.
Flow past one or more cylinders 2.
Flow past porous cylinder

How to Run the Simulations, Visualize, and Analyze
In this section we will give instructions on how to run, visualize, and analyze the simulations produced in this software. We will explain how to change parameters below as well. In particular, we will provide both overviews of what is possible to visualize and analyze as well as a guided walk-through of the steps required, designated by Guide in their title. We broke this section into the following subsections:
Visualizing the Data in VisIt 3.

Running a Simulation
Here we will briefly describe how to run each fluid solver's corresponding simulations. To illustrate the software structure and how to run the simulations, we will use the MATLAB (https: //www.mathworks.com/products/matlab.html) [29] version. Note that the Python implementation is consistent in its structure as well as naming conventions.
To run the simulation, one would need to do the following:

1.
Either clone the Holy_Grail repository or download the Holy_Grail zip file at https://github.com/nickabattista/Holy_Grail/ to your local machine. Note you can download or clone this repository to any directory on your local machine.

2.
Open MATLAB and go to the appropriate sub-directory for the particular fluid solver example you wish to run, e.g., 3.
The main script name for each fluid solver is named accordingly, e.g.,

•
Projection_Method.m -main script for the projection method • FFT_NS_Solver.m -main script for the spectral (FFT) solver • lets_do_LBM.m -main script for the Lattice Boltzmann solver Each of these scripts contains all the functions and modules necessary to run each example. The file print_vtk_files.m produces the .vtk data files during each simulation. The user should not have to change this file, unless they do not wish to save some of the data for storage reasons. If this is the case, this can be accomplished by commenting out the lines corresponding to whichever data is not to be saved.

4.
Depending on the subfolder that you are in, to run an example, you would type its main script name into the MATLAB command window and click enter.

5.
The simulations are not instantaneous; they may take a few minutes. Upon starting the simulation, a folder named vtk_data is produced, which will be used for storing the simulation data files in .vtk format. As the simulation runs, it will dump more data into this folder. Note that depending on the simulation you choose and the grid resolution, the simulations may take from a few minutes to ∼30 min. Moreover, if you wish to run multiple simulations, note that the default folder name is set to vtk_data, and so running additional simulations will overwrite any previously saved data in the previous vtk_data folder. Please manually change the name of the folder after each simulation to avoid this.
Digging deeper into each fluid solver script: 1. If one desires to change parameters for a particular simulation, e.g., grid resolution, viscosity (either µ (Projection) or ν (FFT) or τ (LBM), etc.), options are available to the user without digging deep into any of the sub-functions or individual modules beyond the main function itself, see 3. Note that each example will run as described in Section 4 if all the parameters are selected to be those in the parameter

Visualizing the Data in VisIt
In this section we will briefly describe how produce visualizations, such as those Section 4. We will briefly detail the steps to visualize the data using the open-source software VisIt (https: //visit.llnl.gov) [40]. Note that each simulation produces data in the .vtk format and so one could alternatively use the open-source software ParaView [41] for visualization purposes as well. Later, in Section 3.3 we will guide a user through running an example (the spectral (FFT) method's bubble3 example) followed by a step-by-step guide to visualizing its corresponding data in Section 3.4. We used VisIt version 2.13.3.
To visualize the data, one would need to do the following:
Open the desired Eulerian data (magnitude of velocity, vorticity, velocity vector field, etc.) from the vtk_data folder.

3.
A table of the possible data stored and their corresponding name when saved is found below in Table 1. To add streamlines of the velocity field, first make sure that your Active Source is set to u (see Figure 4a), then click Add→Pseudocolor→operators→IntegralCurve→u (see Figure 4b).
Note you have the option to put numerous seed points for streamlines to stem in the domain.
Here we have chosen to use a line of points, sampled it at 4 evenly spaced locations within that line, and specified particular integration tolerances, and a maximum number of steps (see Figure 4c). Click Draw. You can also adjust the streamline aesthetics to how you desire, e.g., color, thickness, etc..

6.
To Visualize the Eulerian scalar data (e.g., Vorticity, Magnitude of Velocity, etc.): Go to the vtk_data data folder that the simulation produced (c) Click on the grouping of the desired Eulerian scalar data, for example, Omega (for Vorticity), click OK (d) In VisIt, click on Add then Pseudocolor→Omega. (e) Then click Draw (f) You can elect to change the colormap and/or colormap scaling by double clicking on Omega in the VisIt database listing window. The following table (Table 2) details the specific colormap used for each scalar variable in the manuscript,  Figure 5a) must be on the velocity vectors, u. Then click Add→Pseudocolor→ operators→LCS→u (see Figure 5b). We used the attributes in Figure 5c for computing the FTLE. In particular, we limited the maximum advection time to 0.05 and maximum number of steps to 10. Once those are changed, click Draw. (h) To add contours of the desired scalar variable, first make sure that your Active Source (see Figure 5a) is set to the correct variable, then click Add→Contour→<desired variable name>. Note that you have the option to scale the contours separate from the colormap of the variable itself. Moreover, for FTLE your active source must be u and you can modify how FTLE are computed, e.g., Figure 5c; however, we used consistent values for both FTLE computations.

Guide: Running the Spectral (FFT) Method's 'bubble3' Example
In this guided walk-through, we wish to illustrate how to run a simulation and change its parameters while also observing what information is being printed to the screen. We will use the Spectral (FFT) Method's bubble3 example and run two different examples where each uses a different fluid kinematic viscosity, ν.
To run the simulation(s), do the following: 1.
First we must make sure that we are inside of MATLAB the corresponding directory for the MATLAB version of the spectral-based solver.
To do this click through: Holy_Grail→FFT_NS_Solver→MATLAB, see Figure 6. Note that I have downloaded (or cloned) the software into my Documents folder.

2.
To run this simulation, type FFT_NS_Solver into MATLAB's command window and click enter. See Figure 7 for what should print to the command window shortly thereafter. Note that the bubble3 example is the default built-in example upon downloading the software (see Figure 3b).
Upon starting the simulation, information regarding the simulation solver and simulation is printed to the screen. Moreover, the current_time within the simulation gets printed to the screen at the time-points in which the data is stored. The data is stored in the vtk_data folder, which is made upon starting the simulation.

Figure 7.
Showing the output on the screen when this example begins running. Note that the vtk_data folder is produced, in which stores the simulation data at predetermined time-points.

3.
Once the simulation has finished running, we will change the vtk_folder's name to Simulation_1_Data. This folder contains all of the data that was stored during the simulation, i.e., the vorticity, fluid velocity field, magnitude of velocity, etc.. It is imperative to change this folder name after each simulation if more simulations are to be run, or else the new data will printed over the previous simulation's data.

4.
Once the folder's name has been changed, we will open the FFT_NS_Solver.m script to change one (or more) of the parameters. Here we will only change the fluid's kinematic viscosity, ν; however, there are other options in which you could change, such as the computational domain size and/or grid resolution, see Figure 3b. As an example change nu from 1.0 × 10 −3 → 1.0 × 10 −2 , see Figure 8. Once you have changed ν, you can close the script.
Note that if you chose to change the grid resolution and domain size, make sure that the resolution in x and y are equal, i.e., dx = L x /N x = L y /N y = dy. Moreover if you increase the resolution (increase N x , N y ) the simulations will take longer to run and the resulting data will require more storage.

5.
Finally, run this new simulation case by typing FFT_NS_Solver into the command window and pressing enter. Once that simulation has finished running, change its folder name to Simulation_2_Data. Now that we have a couple of simulations run, and have their corresponding data saved, we will next focus on visualizing each simulation's dynamics using the open-source software VisIt [40].

Guide: Visualizing the Spectral (FFT) Method's 'bubble3' Data
Here we will guide the user through visualizing the bubble3 simulation data from Section 3.3. In particular we will visualize the magnitude of velocity data. Note that other Eulerian scalar data, like vorticity, the xor y-component of velocity, etc., may be visualized analogously, as we are following the protocol from the previous section on data visualization, Section 3.2's Visualizing the Eulerian scalar data. To visualize the fluid's velocity vector field, please follow Section 3.2's steps for To Visualize Velocity Vectors.
To visualize the uMag data, do the following: 1.
First we must open the desired magnitude of velocity data. Recall that the fluid's magnitude of velocity data was saved under the name uMag. Click open in the VisIt toolbar and in the dialog box that pops up, locate the desired data directory, here it is Simulation_Data_1 from earlier, and select the uMag group. Figure 9 is provided as a visual aid for these steps.

2.
Although we have opened the data, nothing will appear visualized, yet. However, now uMag is the Active Source in the window, see Figure 10. We must now tell VisIt how to visualize the data. We will visualize magnitude of velocity as a background colormap, as it is a scalar quantity. To do this click 'Add' and then select Pseudocolor from the menu that appears, followed by uMag, i.e., Add→Pseudocolor→uMag.

3.
To finally visualize the data, click Draw. Window 1 will then show the fluid's magnitude of velocity data at time-point 0, i.e., the initial magnitude of velocity, see Figure 11.

4.
Furthermore, we note that certain choices of colormaps are better for both conveying the data, but also for people who are colorblind. Generally, the default colormap, which was used in Figure 11, is a bad choice, as it is impossible for people who are red-green colorblind to interpret the field [57]. Thus, we will change the colormap to another, e.g., the RdYlBu (red-yellow-blue) colormap, by double-clicking on the uMag.*.vtk database:Pseudocolor uMag in the database list and appropriately selecting the desired colormap. Other properties may also be changed in this box, such as the colormap's scaling.
Also, you can visualize the magnitude of velocity at the different time-points that were saved during the simulation by moving the time-slider. Figure 12 shows the data at time-point 30. Moreover, the colormap and its properties may be modified by double-clicking on uMag.*.vtk database:Pseudocolor uMag in the database list.

5.
Note that the time in the simulation does not generally correspond to the time-point stored. For example, in this case the simulation modeled 30 s of time and 60 total data time-points were stored. Thus, data was saved every 0.5 s of simulation time. Furthermore, you could repeat Steps 1-4 for data in Simulation_2_Data from the second simulation you ran. Instead of comparing the 30th time-point's data, we will compare the 60th, i.e., the last time-point. Repeating this entire process for the other simulation's data yields the comparison as shown in Figure 13, below. However, we cannot qualitatively compare the magnitude of velocity between each of these simulations as they are in Figure 13. Each simulation's colormap has a different scale. In order to compare, we must make both scales equivalent.

6.
To change the scale, double-click on each of the databases. This will bring up a dialog window (as described earlier) where you can change the scaling, e.g., minimum and maximum, as well as the colormap itself. Change the minimum and maximum values to 0.0 and 0.2, respectively. Figure 14 illustrates where to change the minimum and maximum for the colormap's scale. Once those values are changed, click Apply. Note that you must do this for both of the databases individually.

7.
After changing the scale, the comparison is significantly different than when first visualized, see Figure 15. Furthermore, we also provide visualizations of other time-points and/or scalar data in Appendix C.

Guide: Analyzing the Spectral (FFT) Method's 'bubble3' Data
To analyze the data in VisIt, we will provide a walk-through of the steps that necessary to measure a fluid quantity (particular variable) across a line of interest within the computational domain using the Spectral (FFT) Method's bubble3 example from Section 3.3. We will assume that you have visualized the data for the magnitude of velocity that corresponds to both of the simulations. In this tutorial we will compare the magnitude of velocity for both simulations across a vertical line through the center of the computational domain at the last time-point stored, i.e., the 60th time-point. See Figure 16 for an idea of where the measurement will take place.
To analyze the uMag data, do the following: 1.
Make sure that Active source is the uMag data corresponding to the data from Simulation_1_Data in VisIt, see Figure 16.
Once the Active Source has been appropriately selected, in the command bar, in the VisIt menu go to Controls→Query (see Figure 16).

Figure 16.
To analyze data in VisIt, the active source and active time slider must be set to specific data that is desired to be analyzed, e.g., here it is the magnitude of velocity, uMag, for the Simulation_1_Data.

2.
Next select Lineout under Queries and then under Variables select Scalars and then uMag, i.e., Variables→Scalars→uMag, see Figure 17.

3.
Next we will define the line segment in which we intend to measure the data across. To do this the user gets to choose the starting point and ending point of the line segment, see Figure 17b. We will measure across a vertical line through the center of the domain, hence we choose (x, y, z) = (0.5, 0.0, 0.0) and (x, y, z) = (0.5, 1.0, 0.0) as the starting and ending point, respectively. We also will measure the value of uMag at 250 sampled points. Thus, change Sample Points to 250, see Figure 17b and select Use sampling. VisIt will automatically interpolate the data so that it is measured at number of equally-spaced sample points that is designated by the user.

5.
To get back to the colormap window, select Active Window 1, rather than 2.

6.
You can repeat this procedure and overlay multiple uMag measured data on the same plot, each corresponding to a different simulation, e.g., repeat this process for Simulation_2_Data. Be sure to choose the correct data for the Active source.

7.
If you repeat this entire process for the Simulation_2_Data and move the time-slider such that the data being shown in Active window 2 corresponds to the last time-point stored (the 60th), you should see what is give in Figure 19. Also, note that we have changed each line's thickness individually.

Built-in Examples (for Each Fluid Solver)
In this section we will highlight a subset of the built-in examples that can be run immediately upon download (or git clone) of the software. These examples were selected as they are either traditional problems in fluid dynamics or problems that could naturally lead to fruitful discussion and analysis of the underlying dynamics. The examples we will discuss are:

1.
Lid-Driven Cavity Flow via the Projection Method 2.
Circular Flow in a Square Domain via the Projection Method 3.
Side-by-side Vorticity Region Interactions via the Spectral (FFT) Method 4.
Evolution of Vorticity from an Initial Velocity Field via the Spectral (FFT) Method 5.
Flow Past One or More Cylinders via the Lattice Boltzmann Method 6.
Flow Past a Porous Cylinders via the Lattice Boltzmann Method For each of these examples, we will provide the necessary simulation details, including a brief background about problem as well as the numerical parameters that were used in each simulation. Furthermore, we will also make observations regarding each simulation's dynamics, investigate some of the data, and suggest future ideas for how students could either modify or explore each example further.
These examples are available within both the MATLAB and Python implementations in the software. Note that Section 3 described how to run, visualize, and analyze the simulations, and thus the visualizations and analysis contained here can be recreated by students and other practitioners.

Cavity Flow (via Projection Method)
In this example we will use the Projection Method in the software to run a cavity flow problem. Note that this example is contained in the Projection Method folder and may be run by selecting the 'cavity_top' option (see line 58 in Figure 3a).
For this lid-driven cavity flow problem, the non-zero horizontal velocity on the top wall, u = u T , is set to u T = 4.0 m/s, while all other tangential (and normal) velocities along the boundaries are set to zero (see Figure 1a). Note that the top wall velocity is not immediately set at u T , but rather the flow ramps up from 0 to the preset u T along this wall, to avoid instantaneous acceleration artifacts. For the simulations shown below with Re = 4000, use the computational parameters listed in Table 3. The Reynolds Number was set by defining the characteristic length scale to be L = L x , the horizontal length of the domain, and the characteristic velocity to be V = u T , e.g., Re = ρL x u T µ . To change the Reynolds Number, the dynamic viscosity was varied. For Re = 4000, 400, 40 and 4, we used µ = 1, 10, 100, and 1000 kg/(m · s) , respectively. Note that for the cases with Re = 4 and 40, we decreased the time-step to dt = 5 × 10 −5 s to ensure numerical stability of solutions. If we had not, errors would have been magnified every time-step of the simulation and/or the numerical solutions may have blown up. In CFD it is common that one may need to vary the time-step depending on the Re or other parameters of the system. Situations in which very small time-steps need to be taken to ensure numerical stability of solutions are known as stiff equations, for more information please see [42]. Table 3. Numerical parameters for the Cavity simulation for Re = 4000.

Parameter
Variable Units Value Total Simulation Time T s 6 Fluid Density ρ kg/m 3 1000 Upon having run the cavity model for Re = 4000, we visualized its corresponding simulation data using VisIt [40], as illustrated by Figure 20, which gives colormaps of vorticity, magnitude of velocity, velocity in the horizontal direction throughout the domain, and pressure, as well as depictions of the velocity field, both as a scaled and non-scaled quantity. For the scaled velocity field snapshot, the scale factor was equal to the largest magnitude of velocity across the entire computational domain. In the non-scaled plot, streamlines are also given. Streamlines show the path that a passive particle would take in the flow at a particular moment in time. The data shown here was taken at t = 6.0 s, when the simulation ended. A fully-formed vortex developed by that time near the top of the cavity, while a smaller vortex appears to be forming below (see Vorticity pane in Figure 20). Snapshots illustrating the formation of such vortices from this case of Re = 4000 are presented in Figure 21, which gives a colormap of vorticity and the velocity vector field. As fluid is moving along the top edge of the domain from left-to-right, the fluid nearest to the top begins moving in the same direction. Eventually fluid down in the cavity begins moving in the same direction as well, until it reaches the right boundary, where it must move downwards. It is easier for the fluid to move downward because of the faster flows moving towards the right above it, due to the boundary condition. As the fluid moves downward along the wall, fluid towards the middle of the cavity begins moving downwards as well. In tandem, with the fluid moving left-to-right along the top and top-to-bottom along the right, these fluid patterns initiate the formation of a clockwise-spinning vortex.
Now that there is a vortex spinning clockwise, it causes the fluid below the vortex (moving right-to-left) to begin moving right-to-left until it reaches the wall on the left and an analogous situation occurs to the above, except an oppositely-spinning (counter-clockwise) vortex begins to form. Recall that these vortices can both be traced back to the fluid moving along the top of the domain left-to-right. Within the region that the second vortex forms, much of the energy (velocity) has dissipated away, resulting in a smaller, less strong vortex.
We can explicitly measure the horizontal velocity descending down the cavity. This data is presented in Figure 22, which gives the horizontal velocity vs. depth in the cavity during four different snapshots in the simulation. We also measure the horizontal velocity across three different lines descending into the cavity. Near the surface of the cavity (at depth = 0 m), the velocity is equal to the boundary condition. As one descends into the cavity, there are spikes in both the positive and negative horizontal velocities. These spikes correspond to locations near the edges of the vortices, where there is faster moving fluid.  Finally, one can change the Re by varying the dynamic viscosity, µ to observe how vortex formation changes within the cavity, as well as, how quickly energy gets dissipated in higher viscosity settings. Recall that if µ is changed to 40 or 4 the time-step, dt, was adjusted for as described above. Figure 23 illustrates how vortex formation changes within the cavity for different Re at three different snapshots during the simulations. Moreover, it gives the horizontal velocity vs. cavity depth, as measured across the middle of the cavity. Note that the Re = 4 and Re = 40 data virtually overlaps in the figure, and velocity (and system's energy) quickly dissipates to zero going further into the cavity. Students may elect to try the following: 1.
Recreate the above results for cavity flow with the geometric and simulation parameters given in Table 3 for one or more Re.

2.
Change the domain width (cavity width) to observe differences as the cavity gets wider for a given Re.

3.
Vary the lid-velocity to observe differences (e.g., varying Re for different velocities as opposed to differing viscosities, µ).

Circular Flow in a Square Domain (via Projection Method)
In this example we will use the Projection Method in the software to run an example involving circular flow in a square domain. Note that this example is in the Projection Method folder and may be run by selecting the 'whirlwind' option (see line 58 in Figure 3a).
In this case, all the tangential boundary conditions along the computational domain were set to be U = u top = −u bot = v le f t = −v right = 1.0 m/s, see Figure 1b. As mentioned in Section 4.1 above, the tangential wall velocities are not immediately set to U, but rather the flow ramps up along each edge from 0 to the preset u top , u bot , v le f t , or v right along this wall, to avoid instantaneous acceleration artifacts. For the simulation shown below with Re = 4000, use the computational parameters listed in Table 4. Note that as Re = ρLU µ , with L = L x = L y and U = 1.0 m/s, one can vary the dynamic viscosity, µ, to change the Re. If you decide to make Re smaller (e.g., increasing µ), you may also need to decrease the time-step, dt, significantly to ensure numerical stability, as mentioned in Section 4.1.
The requirement that very small time-steps are necessary to ensure numerical stability denotes what are called stiff equations. For more information on stiff equations, please see [42]. Having run the circular flow simulation for Re = 4000, we visualized its corresponding simulation data using VisIt [40], see Figure 24, which gives colormaps of vorticity, magnitude of velocity, velocity in the horizontal and vertical directions, and the finite-time Lyapunov exponent (FTLE). The FTLE can be used to find the rate of separation in the trajectories of two infinitesimally close parcels of fluid. Maxima in the FTLE have been used to determine approximations to the true Lagrangian Coherent Structures (LCSs). LCSs are used to determine distinct flow structures in the fluid [58][59][60][61] can be used to divide the fluid's complex dynamics into distinct regions to better understand transport properties of flow [62][63][64]. True LCSs would require knowledge of the infinite-time Lyapunov exponent (ITLE); however, the FTLE can serve as a proxy to the ITLE. In this paper, we computed the forward-time FTLE field, whose maximal ridges give approximate LCSs corresponding to regions of repelling fluid trajectories and whose low values give rise to regions of attraction [61]. Thus, the FTLE can be used to help find regions of high fluid mixing [58,59,61,65,66]. Figure 24 also provides contours for each quantity. It also includes a depiction of the fluid velocity field, as scaled by the maximum in the entire domain's magnitude of velocity, which also includes its associated streamlines. Recall that streamlines are curves that depict instantaneous tangent lines along direction of the flow velocity. They show the direction that a neutrally-buoyant particle would travel at a particular snapshot in time. These are different than contours (also known as level-sets or isolines), which are curves where a function has a constant value. This snapshot was taken at t = 24.0 s, when the simulation ended. As the velocity boundary conditions continually increased, the interior of the domain began to move clockwise, in the same direction as the flow at the boundaries (see Figure 1b). The colormaps and contours in each velocity-related panel reveal that flow speeds decrease towards the interior of the domain. Moreover, as all the boundary conditions are uniform, the vertical and horizontal velocity plots are identical under a 90 degree rotation. Also, the vorticity plot illustrates that vorticity is still low by the ending of the simulation in the center of the domain, even though flow across the whole domain is rotating clockwise. The FTLE plot suggests that the majority of fluid mixing occurs near the edge of the domain, rather than the interior, as higher values of FTLE suggest nearby fluid trajectories move away from each other at higher rates. This is due to larger spatial gradients in the velocity in these regions, where flow is moving at different speeds in different directions.
We also present data comparing simulations for Re = 400, 1000, and 4000 (for µ = 2.5, 1.0, and 0.25 kg/(m · s), respectively). Figure 25 illustrates colormaps of the magnitude of velocity (with corresponding contours) at various time-points during the simulation. Every colormap uses the same scaling in the figure. Higher velocities appear to creep into the interior of the domain quicker in the lower Re cases. By t = 24 s it appears that the Re = 400 and 1000 cases look almost identical, the Re = 4000 case tells a different story -the interior of the domain is still moving considerable slower than the other two cases. This might seem counter-intuitive at first, as one might generally think that higher Re tends to lead to more fluid motion, especially when we are lowering viscosity to increase Re. However, this is due to differences in time-scales for the diffusion of momentum through the fluid. The viscous diffusion time can be thought of as the time it takes for a fluid parcel to diffuse a particular distance on average. Ift is the diffusion time, ν = µ/ρ is the kinematic viscosity, andl is the mean-squared distance, thent if all of the simulation geometries are identical and only µ (and hence ν) is varied. This concept is based off viewing diffusion as a random walk process [67]. Hence for the simulations presented here the time-scales vary betweeñ t Re=400 ∼ 1 2.5/1000 = 400 s andt Re=4000 ∼ 1 0.25/1000 = 4000 s. Therefore in the Re = 4000 case, the diffusion time-scale is 10x longer than that of the Re = 400 case, and thus the dynamics evolve much slower in the Re = 4000 case! Moreover, this phenomena is shown in Figure 26, which illustrates the fluid velocity field along with its associated streamlines at the same snapshots in time. This is more apparent if we compare flow profiles within the domain. Figure 27 compares the horizontal flow profile for Re = 400, 1000, and 4000 across the vertical mid-line of the domain. As mentioned above the dynamics in the Re = 4000 case evolves about 4x and 10x slower than the Re = 1000 and Re = 400 cases, respectively. By 24.0 s, the flow profiles in the Re = 400 and Re = 1000 cases look almost identical, while the flow profile in the Re = 4000 resembles that of the the flow profile in the Re = 1000 at ∼6.4 s-approximately a factor of 4 different in time. Furthermore Figure 28 shows the flow profiles for all three cases of Re when the diffusion of momentum has reached the same depth. This occurs at different simulation times due do the variations in viscous diffusion time-scales.  Students may elect to try the following: 1.
Recreate the above results with the parameters listed in Table 4 2.
Change the computational geometry, e.g., rectangular instead of square 3.
Vary the input velocity along each side of the domain to make them non-uniform 4.
Vary the Re by changing dynamic viscosity for suggestion (2) or (3)

Side-by-side Voritices (via Spectral Method)
In this example we will use a FFT-based fluid solver in the software to run a simulation of multiple regions of vorticity interacting with one another. Note that the first example given here can be run by going into in the FFT_NS_Solver script selecting the 'qtrs' option (see line 47 in Figure 3b). The vorticity is initialized as in Figure 2a. We will also highlight the 'half' option in this section as well.
Recall that for this solver, an initial vorticity configuration is needed. For the simulations of four interacting vorticity regions, the CW and CCW vorticity were initialized as ω = −1 and 1 rad/s, respectively. In other cases with only two interacting voricity regions, the strengths of each regions were varied between −1, −0.5, 0.5 and 1 rad/s, as appropriate. For the simulations shown below with either 2 or 4 interacting vorticity regions, use the computational parameters listed in Table 5.  Figure 29 provides the simulation data for the case with 4 interacting vorticity regions, as described above. It presents colormaps for vorticity, magnitude of velocity, horizontal velocity, vertical velocity, and the finite-time Lyapunov exponent (FTLE), which is used to determine approximations to the true Lagrangian Coherent Structures (LCS) to which can be used to distinguish fluid mixing regions [58,59,61,65,66]. Their corresponding contours are also given. Recall that contours are curves to which the value of a quantity is constant. Figure 29 also gives a snapshot of the velocity vector field. This data is from the last snapshot of the simulation at t = 5.0 s.
The top two and bottom two vorticity regions are initialized the same, e.g., clockwise (CW) and counterclockwise (CCW), respectively, by this point in the simulation they are beginning to lean in the direction of rotation. Moreover, it can be seen that fluid is moving quickly horizontally through the top, middle, and bottom of the domain, as illustrated by the magnitude of velocity plot. When this is compared to the horizontal velocity plot, one observes that the fluid along the top and bottom of the domain both are moving left-to-right, while fluid moving horizontally through the middle of the domain is moving right-to-left. These directions align with the direction of the spinning vortices in each of these horizontally-aligned regions -top, middle, and bottom, which can be seen in the velocity vector field depiction.
The high values of FTLE illustrate regions of higher fluid mixing, as higher FTLE suggest faster rates of separation of close fluid blobs. In Figure 29   One could elect to change the size or strength of a subset of these vorticity regions to simulate asymmetric vorticity interactions in this configuration. However, rather than study 4 interacting vorticity regions, we compared 2 interacting regions at a time, and varied the strength (magnitude of vorticity), size, and initial vorticity values (spinning-direction) between them, see Figure 30. Figure 30 provides comparison data of snapshots at t = 6.0 s in each simulation for different quantities -either vorticity with velocity vectors, magnitude of velocity with its corresponding contours, and FTLE with its corresponding contours. We use the language that a "strong" vorticity region has a |ω| = 1 rad/s, while a "weaker" region has a vorticity of |ω| = 0.5 rad/s. Each initial size of the vorticity regions had a radius of 0.15 m or 0.3 m.
In general, one can notice the following:

1.
When same-spinning and size vorticity regions interact, more fluid mixing occurs directly between them (1st row) 2.
When oppositely-spinning and same size vortex regions interact, there is enhanced vertical fluid flow between the vorticity regions, but less mixing between them (2nd row) 3.
When oppositely-spinning vorticity regions of different strengths interact, the stronger region (whose magnitude of vorticity is greater), influences fluid flow to move in the same direction further away from it. There is still some enhanced vertical flow through the region between the regions. (3rd row) 4.
If there are asymmetric region sizes and strengths, where the larger region is weaker (smaller magnitude of vorticity), the region with the larger vorticity magnitude may be able to keep the other region from influencing much of the flow around it, while imposing its flow direction on the weaker region. (4th row) 5.
Two regions of vorticity with the same sign, but differing strengths, leads to enhanced fluid mixing between the vortices, while flow around the smaller region are significantly less than the other. (5th row)

6.
If there are asymmetric region sizes but uniform strengths, enhanced vertical flow will be seen between the regions, although the larger vorticity region exerts more influence on flow patterns closer to the smaller region. (6th row) 7.
If there are asymmetric region sizes and strengths, where the smaller region is weaker (smaller magnitude of vorticity), enhanced vertical flow will be seen between the regions, although the stronger vorticity region exerts more influence closer to the smaller region. (7th row) Figure 30. Vorticity and velocity field (left column) and magnitude of velocity with its associated contours (right column) for cases of two vorticity regions with different sizes, strengths, and vorticity initialization.
Students may elect to try the following: 1.
Recreate the above results for interacting side-by-side vorticity regions with the geometric and simulation parameters given in Table 5 or for different ν or computational domains.

2.
Add more vorticity regions or make the case with four asymmetric in terms of placement and/or size.

3.
Change the distance between the vorticity regions to observe how magnitude of velocity or FTLE contours change.

Evolution of Vorticity from an Initial Velocity Field
In this example using a spectral (FFT) solver, we begin with a snapshot of the velocity vector field from an independent CFD simulation, see Figure 31a,b. Note that this snapshot was taken from an open-source fluid-structure interaction example of a pulsing cartoon heart [53]. We used a stored time-point's velocity field data, to which we extracted the horizontal and vertical velocity components, and then computed its associated vorticity at that particular time-step, see Figure 31c. We note that there are numerical errors in computing the vorticity, ω, as a centered finite difference scheme [42] was used to compute each first-order partial derivative in calculating vorticity, i.e., ω = ∂v ∂x − ∂u ∂y . Moreover, since this fluid solver uses periodic boundary conditions, we computed the partial derivatives at each boundary using data from the opposite side of the computational domain. Note this example may be run by selecting the 'jets' option in the FFT_NS_Solver script (see line 47 in Figure 3b). After having computed the vorticity from a velocity field, the simulation could begin. This simulation used the computational parameters listed in Table 6. We further note that the velocity field had an original grid resolution of [N x , N y ] = [512, 512], so we down-sampled it to a [256, 256] grid for the example shown here. Furthermore, as the original velocity field was computed from an independent fluid-structure simulation, we comment that there is no complex, moving boundary within the computational domain; however, at the initial point you can see the remnants of the coinciding heart structure, see Figure 31.  Figure 32 provides snapshots over the simulation to illustrate how the vorticity and magnitude of velocity evolve. Both Figure 32a,b illustrate the effect of periodic boundary conditions along each edge of the domain, e.g., vortices moving through the top or right side of the domain continue through the bottom or left side of the domain, respectively. Furthermore, in the snapshots provided, the overall vorticity and magnitude of velocity appears to dissipate within the domain . This is because there are no source terms (or explicit boundary conditions) providing additional flow (energy) into the system to induce more flow. Although the original velocity vector field that defined the initial vorticity of the system came from a simulation of pulsing hearts, without the heart moving and thus pushing on the fluid, the fluid will eventually stop moving and reach zero velocity due to the fluid's viscous nature. Therefore initializing the vorticity out of a time-point from another simulation's velocity field, is akin to studying the evolution of the fluid's motion from a singular impulse in the flow. Students may elect to try the following: 1.
Recreate the above results with the parameters listed in Table 6 or for different ν.

2.
Take a velocity field from one of the Projection method examples, e.g., from Section 4.1 or 4.3, and use it to initialize the vorticity for a simulation using this spectral (FFT) method.

3.
Find how long it takes the simulation to reach a 1/10 or 1/100 of the original maximum magnitude of velocity.

4.
Discover how that time varies with changes in viscosity, ν.

Flow Past One or More Cylinders (via Lattice Boltzmann)
Flow past objects is one of the most heavily studied problems in aerodynamics and thus is a standard example used in many fluid dynamics courses [68][69][70]. While many texts only consider the problem for inviscid flows and use potential flow methods to solve for exact solutions [69,70], here we include effects of viscosity to illustrate transitions to vortical flow, and upon doing so, a subset of the possible flow separation cases. Furthermore, flow past cylinders is also still an active area of contemporary research [71][72][73][74][75][76][77].
We performed simulations for both a single cylinder as well as multiple cylinders. The computational geometry for a single cylinder case is given in Figure 33. Figure 33a gives the boundary conditions considered for the problem for flow past a rigid cylinder contained within a rigid channel. The walls of the channel are rigid and so bounce-back boundary conditions are used. The ends of the channel use periodic boundary conditions, e.g., what goes out the right side, comes in through the left. Note that this is the first example discussed in which uses a mix of explicit boundary conditions and periodic boundary conditions. Moreover, the cylinder itself is modeled as a rigid structure and so bounce-back conditions are used on it. Figure 33b provides visual details how the inflow is modeled in the simulation. Every iteration, the horizontal velocity at the inflow is increased by an amount of ∆U x to drive fluid flow towards the right. All simulation parameters (fluid, grid, and geometry) for both cases involving either a single cylinder or multiple rigid cylinders are given in Table 7. Note that the cases with multiple cylinders simulate flow around cylinders with larger radii, use a larger value for ∆U x , and a different final number of steps and grid resolutions. Furthermore, due to these differences we are not comparing the case of one single cylinder to the case with multiple cylinders. The simulation data obtained at step n = 49, 600 is given in Figure 34. It presents colormaps (and corresponding contours) for vorticity, magnitude of velocity, horizontal velocity, vertical velocity, and the finite time Lyapunov exponent (FTLE), to which knowledge of the approximate Lagrangian Coherent Structures (LCS) can be extracted [58,59,61,65,66]. We attribute regions with high FTLE to regions of high fluid mixing, as higher FTLE suggests that nearby fluid parcels separate at a much faster rate than those in lower FTLE regions. From the vorticity panel, significant flow separation is seen as vortices are being shed off the cylinder and flow pushes past it left-to-right. Oppositely spinning vortices are shed off cylinder one after another. Moreover, we observe that there is significantly more horizontal fluid motion than vertical within the channel. Since flow is being pushed left-to-right from the inflow condition, it would require more energy for fluid to move vertically, than to simply. . . go with the flow. Figure 35 provides snapshots of vorticity (left) and FTLE (right) during the simulation to highlight how smooth flow developed into vortical flow patterns once the horizontal velocity increased past a threshold. Such threshold appears to have occurred around simulation step ∼ 47000.
In particular, flow instabilities began to manifest by step 46000, as seen by both the vorticity and FTLE snapshots. The majority of fluid mixing occurred in large contours behind the cylinder up to this point. Once flow separation occurs, some mixing occurs in the wake in small patches, but no geometrically long regions of higher fluid mixing are present, as seen by the contours. Note that in the last two panels, for steps 49200 and 49600, large FTLE-valued regions do not correspond to locations of shed vortices, rather higher FTLE valued regions appear between oppositely spinning vortices.  Finally we performed simulations with multiple cylinders in a channel. Each cylinder in these simulations had a radii 66% larger than the single cylinder shown above. We performed three separate simulations with different configurations of the cylinders, each based off a triangle formation. Figure 36 illustrates how flow evolved for each geometric configuration, in particular for (a) vorticity and (b) FTLE.
Vortices are observed shedding off each cylinder, possibly at an enhanced rate due to flow interactions with the other cylinders (see suggested activity below). In the symmetric configurations, e.g., the triangle-formation and vertical-line formation, flow symmetry is preserved, while in the case with only one trailing cylinder, flow asymmetries arise. FTLE plots at the last step shown (5000) illustrate different geometrical configurations not only can lead to different patterns of fluid mixing, but also enhanced mixing, even among cases with only two cylinders. The vertically-aligned cylinder case appears to elicit more downstream mixing in the wake than the case of asymmetrically placed cylinders; however, this may also be an artifact of overall geometry of the system, e.g., the size of the cylinders and the proximity of cylinders to the channel walls. Students may elect to try the following: 1.
Recreate the above results with the parameters listed in Table 7 2.
Using the radius for the larger geometry case, run a simulation with only one cylinder present in the middle and compare vorticity and FTLE evolution plots.

3.
Vary τ (relaxation parameter relate to viscosity) to test the system's sensitivity to τ 4.
Change the number, placement, or size of each cylinder in the domain.

5.
Change the shape of the object in the channel, e.g., try an ellipse or airfoil-like shape

Flow Past a Porous Cylinder (via Lattice Boltzmann)
While flow around solid objects is a popular problem in aerodynamics, flow around porous objects has only recently begun to be investigated within in the past few decades, using either high fidelity numerical simulation or experiments [78][79][80][81][82]. Here, we present an example of flow past a porous cylinder using the LBM. The cylinder geometry, computational setup, and inflow/boundary conditions are identical to the single cylinder case of Section 4.5 (see Table 7). The only difference being that the cylinder geometry itself is porous.
The porosity, or void fraction (VF), of the cylinder was chosen to be VF ∈ {0%, 10%, 25%, 50%, 62.5%}. To model the void fraction, each grid cell inside the cylinder was given a random number between [0, 1] and then depending on the specific VF, if that grid cell's randomly assigned number was above the VF in decimal form, it was assigned as a solid boundary. For example, for VF = 25%, each grid cell with a randomly assigned number greater than 0.25 was declared a solid boundary point. Figure 37 depicts the porous cylinders considered for each VF case studied below.
Although, the cylinder is no longer uniformly solid, in this example we will not focus our efforts to study flow through the cylinder itself. For example, in Figure 37's VF = 10% case, there are holes in the interior of the cylinder domain that are not connected to the fluid region outside of the cylinder. Those void regions will not influence the flow outside of the cylinder, as they are blocked from outside fluid penetrating them. However, if you wanted to study the fluid dynamics through a porous structure, you could design a specific void network structure within the cylinder or in another desired geometry. What we will emphasize here is that the porous regions connected to the outside of the cylinder cause asymmetries to develop in the overall flow pattern at an accelerated rate as compared to the non-porous (solid) case. Such asymmetries then lead to quicker vortex shedding.  Figure 38 shows the evolution of vortex shedding in cases of differing porosities (void fractions) of {0%, 10%, 25%, 50%, 62.5%}. Every case shows that more porosity leads to faster development of vortex shedding. Furthermore, even the case with VF = 10% has pronounced vortices being shed before the VF = 0% case even shows any signs of significant flow asymmetry. The increased porosity accelerates flow asymmetries to manifest, leading to vortex shedding occurring earlier on. This example highlights how small perturbations in fluid dynamics problems (here small differences in geometric structure) can lead to the system having significantly different time-scales for flow structures to develop or for a bifurcation in the resulting dynamics to emerge. Students may elect to try the following: 1.
Recreate the above results with the parameters listed in Table 7 2.
Create the accompanying FTLE plots and discuss how fluid mixing changes due to the porous cylinder 3.
Vary τ (relaxation parameter relate to viscosity) to test the system's sensitivity to τ 4.
Increase the number of porous cylinder in the domain and vary their placement or size 5.
Change the shape of the porous object in the channel, e.g., try an ellipse or airfoil-like shape 6.
Create a specific porous network structure through the cylinder and test how flows directly through a porous cylinder may affect vortex shedding

Discussion
In this work we provide software that was developed for the specific purpose to make CFD accessible to undergraduate (and graduate) students in order to provide them an opportunity to perform traditional and contemporary CFD simulations as a valuable learning experience. To that end all the fluid solvers contained within were written in both MATLAB and Python, two languages that are commonly familiar to most science and engineering students [24][25][26][27][28]31]. Students also have the opportunity to modify existing examples or create their own examples within the software to test hypothesis or for stimulate further scientific curiosity. To that extent, we provided an overview of how the code is structured (Section 3) and offered interesting variations to each of the examples that we showcased in this paper (Section 4).
Not only can students run CFD simulations, they also have the chance the practice the art of scientific visualization, and use open-source software that is used by many researchers (VisIt [40] or Paraview [41]), to visualize the corresponding data. These software packages are both open-source and have easy-to-use and learn graphical user interfaces (GUIs) that streamline the visualization process. They can also be used to perform data analysis as well, which further reduces the computational learning curve for students to be able probe the data produced. We provided step-by-step guides on how to use VisIt for these tasks in Section 3.
By performing their own CFD experiments, fluid dynamics students have the opportunity to vary parameters (e.g., fluid scale (Re), boundary conditions, geometry, etc.) of a given system and observe the resulting dynamics. More specifically, they are able to decrypt how the dynamics may vary across parameter regimes, both qualitatively and quantitatively. This grants them the ability to challenge their own developing intuition of the resulting flow physics to further accelerate their learning. They also gain more computer experience and witness first-hand some of the advantages that complex computer simulations offer.
If students wish to go a step further and study the underlying mathematical structure of the code, we provide insightful comments everywhere possible within each fluid solver script. To complement this, we provided a detailed mathematical overviews of each fluid solver used here in Appendix B. While it is not our mission here to teach students how implement their own numerical fluid solvers, we believe this work contributes to providing a foundation for capturing the importance (and utility) of different fluid solver schemes. For each solver that was introduced for particular applications in Sections 2.1-2.3, we give a high-level overview of the method. If students wish to continue learning about numerical methods for solving the fluid equations or numerical methods for partial differential equations in greater depth, we suggest Barba and G. Forsyth's CFD Python: the 12 steps to Navier-Stokes equations [33], L.A. Barba and O. Mesnard's Aero Python: classical aerodynamics of potential flow using Python [34], both of which use Jupyter Notebooks [37], or Pawar and San's [38] modules for developing fluid solvers in the Julia programming language [39].
In this day and age, as the importance of computer literacy increases rapidly [83][84][85][86], integrating more hands-on computer-based activities into course structures will be of critical importance. However, students also seeing successful execution of code is paramount [85]. Here we provide gateway software for students to dive into the world of CFD, in familiar programming environments. While proprietary software offers immense advantages to running simulations and analyzing data, students may see a disconnect between the programming knowledge they've acquired and such polished software. We hope to provide them a unique opportunity to see fluid solvers written in familiar, non-intimidating environments that they may be able to tweak and modify comfortably. Our hope is to inspire further ownership of their learning and to stimulate more interest in (lucrative, transferable) computer programming skills. Acknowledgments: The author would like to thank Laura Miller for introducing him to the joys of computational fluid dynamics. He would also like to thank Christina Battista, Robert Booth, Christina Hamlet, Alexander Hoover, Matthew Mizuhara, Arvind Santhanakrishnan, Emily Slesinger, and Lindsay Waldrop for comments and discussion. He also would like to sincerely thank the Reviewers, especially Reviewer 2, whose above and beyond effort in providing incredibly thorough, insightful, and constructive feedback led to the manuscript becoming significantly stronger. divergence-free, and hence the necessary incompressible condition will not be satisfied (Equation (2)). In discretization terms, this step takes the form of where u * is an intermediate (auxiliary) velocity field, which is not divergence-free. The second step is known as the projection step, where the pressure gets reintroduced to give a final velocity field, u n+1 that satisfies the incompressibility condition (Equation (2)). This discretized step takes the following form, where p n+1 is the pressure field at the next time-step. Because it requires an updated pressure term, we must first find such a pressure. To do this we recall Helmholtz-Hodge Decomposition [87,88], which says any vector field, that is twice continuously differentiable on a bounded domain, say v, can be decomposed into a solenoidal part (divergence-free) and an irrotational part (curl-free), i.e., where v sol is the solenoidal part and v irr is the irrotational part. We note that an irrotational vector field can be written as the gradient of a scalar, e.g., v irr = ∇φ, where φ is some scalar function (sometimes φ is referred to as a potential).
Note that if we take the divergence of (A3), we obtain, It is then possible to find the divergence-free part of the vector field v by solving the above Poisson problem in (A4). This motivates the form of the second step for a projection method given in (A2). However, to find the pressure, we take a divergence of (A2) and note that we require that u n+1 be divergence-free, i.e., ∇ · u n+1 = 0.
Taking the divergence of (A2) and requiring the condition in (A5), we obtain the following Poisson problem for the pressure, p n+1 , in terms of the intermediate velocity field u * , Note that this equation can be solved explicitly. Hence once p n+1 is found, we can then solve for u n+1 using (A2), e.g.,

. Spectral Methods via Fast Fourier Transform (FFT)
For the spectral (FFT) method fluid solver we choose to work in the vorticity formulation of the viscous, incompressible Navier-Stokes equations. Here we will begin by deriving such equations from the previously written viscous, incompressible Navier-Stokes equations, i.e., Equations (1) and (2), ρ ∂u ∂t + (u · ∇)u = −∇p + µ∆u and ∇ · u = 0.
We will use a lot of identities from vector calculus. First, recall the following identity for any vector F, Substituting the above identity into the conservation of momentum equation (Equation (1) and dividing by ρ gives us, where ν = µ ρ is the kinematic viscosity. Note that u · u is a scalar quantity, and thus we define it to be U 2 = u · u. Moreover, we also get to define a new quantity called the vorticity, It is tempting to think of ω ω ω as describing the global rotation of the fluid, but this is misleading. Although many flows can be characterized by local regions of intense rotation, such as smoke rings, whirlpools, tornadoes, or even the red spot on Jupiter, some flows have no global rotation, but do have vorticity. Vorticity describes the local spinning of a fluid near a fixed point in space, as seen by an observer in the Eulerian framework.
Substituting the definitions of vorticity and U 2 into Equation (A8), we obtain Now taking the curl of the above equation we get an equation for the evolution of the vorticity, yields since ∇ × (∆u) = ∆(∇ × u) = ∆ω ω ω. Furthermore we note that the pressure terms drop out as the resulting force from pressure only acts perpendicular to the surface of a fluid blob and not parallel to it, i.e., ∇ × (∇Φ) = 0 for any scalar field Φ. The viscous, incompressible Navier-Stokes equations can thus be written as follows, Shortly we will introduce a vector potential to strive towards introducing a streamfunction, ψ, into our formulation, but first we will use a vector calculus identity to re-write Equation (A11) into a more traditional looking advection-diffusion equation. Using the following identity from vector calculus, we can mathematically massage Equation (A11) into the following form, ∂ω ω ω ∂t + (u · ∇)ω ω ω = (ω ω ω · ∇)u + ν∆ω ω ω.
Note that the evolution equation for vorticity, Equation (A14), now looks like an advection-diffusion equation, but with an additional extra term, (ω ω ω · ∇)u. For 2D flows, recall that u = (u, v, 0) and hence ω ω ω = (0, 0, ω). Moreover, in 2D flows, all partial derivatives with respect to z are zero; hence the (ω ω ω · ∇)u term becomes zero, e.g., Thus, in 2D, we have the following form of the momentum equation in terms of vorticity Note that the form of Equation (A15) suggests that if ν ≡ 0 and if ω ω ω = 0 everywhere at any moment in time, then ω ω ω = 0 for all future time. Moreover, since ω = ∇ × u = 0, we have irrotational flow. Thus, we would be studying an incompressible potential flow problem [70]. Next we introduce the streamfunction, ψ, as part of the vector potential for u, Note that if the streamfunction, ψ = ψ(x, y), is known, it is possible to extract the components of the 2D fluid velocity field, u = (u, v), from it e.g., Furthermore, taking the curl of (A16), we are able to get a Poisson problem for ψ in terms of ω ω ω, where we have used the following vector calculus identity, and the fact that ∂ψ ∂z = 0. At this point, the idea is that if we are able to solve for the streamfunction, ψ, from the vorticity, ω, we can then get the fluid velocity ,u, and it will automatically satisfy the incompressibility condition by definition of the vector potential, i.e., Equation (A16). In essence this is the algorithm; however, within this algorithm, we will work as much as possible in the Fourier frequency space, granted to us by taking the Fast Fourier Transform (FFT). Before diving into the 4 main steps of this scheme, we will introduce some notation involving Discrete Fourier Transforms (DFT) and hence FFT. Note that the FFT yields the same results as the DFT but does so in a more computationally efficient, i.e., fast, manner.

•
Taking the Discrete Fourier Transform (DFT) of a set of complex numbers produces another set of complex numbers. The notation F {z} denotes taking the Discrete Fourier Transform of a set of N-values, {z n } N−1 n=0 , e.g., for k = 0, 1, . . . , N − 1, we definê z n e −2πi kn N .
• Variables with a hat, such asẑ k , denote variables that have been transformed into frequency space via the Discrete Fourier Transform. Moreover, the indices k are known as the DFT's wave-numbers.

•
We can also take the Inverse Discrete Fourier Transform (IDFT) to return our quantities of interest from frequency space to real space. We denote the IDFT as for all n = 0, 1, 2, . . . , N − 1. • Next we will define the discretized quantities in the algorithm. A quantity such as f n ij denotes that quantity's value at the nth time-step at spatial location (x i , y j ) within the rectangular computational grid. Hence we have a set of numbers { f n ij } (or matrix that changes as n → n + 1), and can transform it in analogous manner using a DFT. • We define K X and K Y to be matrices of the DFT's wave-numbers, where K X , K Y ∈ R N x ×N y . They are defined as: Note that each row of the matrix K X is a vector that we denote k X with components k X = (0, 1, 2, . . . , N y − 1) T . Similarly each column of the matrix K Y is a vector that we denote k Y = (0, 1, 2, . . . , N x − 1).

•
In order to benefit from the FFT we will assume our spatial grid has a resolution of N x × N y , where both N x and N y are powers of 2 [89].
We will now dive into the the 4 main steps in this spectral (FFT) method's algorithm. They are as follows:

1.
Update the streamfunction to the current time-step, n: From the previous time-step's vorticity, ω n , we can solve the Poisson problem (Equation (A18)) for the streamfunction at the current time-step, ψ n , i.e., where k X i and k Y j are the Fourier wave-numbers, e.g., the ith and jth components of the vectors k X and k Y , respectively.

2.
Obtain the components of velocity and the gradient of vorticity: With the newly updatedψ n as well asω n from the previous time-step, we are able to compute the components of the velocity field at the current time-step, u n = (u n , v n ). To do this, we take derivatives of the streamfunction and vorticity in frequency space. This then gives us the components of velocity (see Equation (A17)) and the gradient of vorticity, in frequency space respectively. We can then use the IDFT to transform these quantities back into in real space, e.g., where the operation A • B between two matrices of equal size is called the Hadamard product. The Hadamard product is element-wise multiplication, e.g., if A and B are matrices of the same size, for all components i, j, (A • B) ij = A ij B ij .
Once you have the velocity field (u n , v n ) and partial derivatives of vorticity ω n x = ∂ω n ∂x and ω y = ∂ω n ∂y at the current time-step, it is now possible to compute the advection term from Equation (A15), i.e., u · ∇ω ω ω. We define F n adv ij to be the above advection term, and hence get that F n adv ij = u n ij · ω n x ij + v n ij · ω n y ij .
Furthermore, we can apply the DFT to Equation (A25) to transform it into frequency space, e.g., and these streaming velocities, {e i }, are called the microscopic velocities. The directions illustrated in Figure A1 is commonly called the D2Q9 Lattice Boltzmann Model.

2.
The second step involves finding what is referred to as the equilibrium distribution. This step is a part of the collision step, where you want to relax the particle density distributions towards a local equilibrium. The local equilibrium is denoted f eq i (x, t). First we must compute macroscopic the properties (density and velocity) from the intermediate particle distributions f * i using (A29) and (A30).
Before we mention how to handle boundary conditions we will briefly discuss some of the advantages of the LBM. One of the biggest advantages of LBM is its implementation lends itself toward massive GPU or CPU parallelization. Due to parallelization it can be an incredibly fast way of solving fluid problems that are coupled with equations that model heat transfer or chemical processes [91]. Moreover, the algorithm also prides itself for the ability to compute flows through complex geometries and porous structures rather easily and efficiently [55]. From the structure of the streaming step, one can easily prescribe boundary conditions and regions in the grid where fluid is not allowed to flow easily. For our considerations here we only will introduce what are referred to as bounce-back boundary conditions [55].
The bounce-back boundary conditions are used to enforce no-slip conditions; however, as we will show, they are not only used on the edges of the domain, but can be implemented on the interior to create complex geometries. In a nutshell the incoming streaming directions of the distribution functions are simply reversed when they hit a boundary node. This idea is depicted in Figure A3. In practice, one can simply mask these boundary points on the domain using boolean logic. Figure A3. Illustration of bounce-back boundary conditions. During the pre-streaming step there are microscopic velocities set on the boundary and then they are reversed during the streaming step

Appendix C. Extra Visualizations of Data from Section 3: Spectral (FFT) Method's bubble3 Example
These visualizations are to complement those already presented in the guided tutorials from Section 3. We will also give a bit of background for the simulation as well. This example used the spectral (FFT-based) fluid solver in the software and models multiple regions of vorticity 'overlapping' at the beginning. Note that first example given here can be run by going into in the FFT_NS_Solver script selecting the 'bubble3' option. The vorticity is initialized as in Figure A4. Recall that counterclockwise (CCW) and clockwise (CW) correspond to regions of uniform vorticity, where vorticity initialized as a positive or negative constant for CCW and CW, respectively. In this example, three circular regions of vorticity are placed, partially on top of one another. The largest region begins with a uniform value of +0.4, followed by the smaller regions with −0.5 and +0.5, respectively. The remainder of the computational domain is initialized with a random value of vorticity between [−1, 1]. We note that initializing a random values of vorticity will generally not satisfy the incompressibility condition (Equation (2)); however, here we do it only to illustrate that the solver is able to handle initial random noise. This simulation used the computational parameters found in Table A1.  Figure A5 provides the simulation data at the beginning of the simulation (a) and at the last time-step (b). It presents colormaps (and corresponding contours) for vorticity, magnitude of velocity, horizontal velocity, vertical velocity, and the finite-time Lyapunov exponent (FTLE). It also gives a snapshot of the velocity vector field. The last snapshot of the simulation was from t = 30.0 s. From Figure A5, it is evident that the random vorticity values at the start of the simulation eventually interact and dissipate in the flow, as observed in the vorticity panel. Initially there appears to be a lot of noise in the background vorticity (it was initialized to be random between [-1,1]) and by the end it appears virtually averaged out. Moreover, due to the background vorticity noise at the beginning, there are a lot of tiny patches of oppositely moving fluid. This leads to an initial background of high FTLE values, which suggests there is significant fluid mixing occurring. Similarly to vorticity, by the end, the background has significantly less mixing (e.g., smaller FTLE values) overall in areas away from the interacting vortical structures, which started off as overlapping vorticity regions. The area of high mixing in the FTLE panel at the last time-step is in a region where there is a lot of oppositely moving fluid, both horizontally as well as vertically.  Figure A6 provides snapshots over the simulation to illustrate how the magnitude of velocity and vorticity evolved over time. Figure A6a shows that the overlapping regions of vorticity induce the highest flows overall, as quantified by magnitude of velocity. That is, although the background was initialized to random values between [−1,1], which may include values that are higher than initial vorticity of the overlapping regions (+0.4, −0.5, and +0.5 for largest to smallest, respectively), it did not significantly contribute to bulk flow within the domain. Moreover, as suggested earlier, the initial random vorticity configuration quickly dissipates itself out, see Figure A6b. Students may elect to try the following: 1.
Recreate the above results with the parameters listed in Table A1 or for different ν or computational domains.

2.
Change the initial vorticity in each 'overlapping' region.

3.
Change the placement of where each vorticity region is. 4.
Modify the example to have more/fewer overlapping regions of vorticity.

5.
Try initializing a completely random background vorticity without any other vorticity structures.