I do like CFD.
Useful books on Computational Fluid Dynamics
Free CFD Codes
Facebook share

More sample programs are coming up: Delaunay triangulation code, a panel code for an airfoil, 2D unstructured Navier-Stokes code, etc.
[Note: As the name, Katate Masatsuka, implies, I write only when I find time.]

"I do like CFD, VOL.1, Second Edition" is now available in both printed and PDF versions. Download PDF (FREE) at cfd-boook page.

CFD seminar at National Institute of Aerospace is now broadcast online!
Check out the website for seminar videos and files - NIA CFD Seminar


Favorite Codes and Subroutines      


Educationally-Designed Unstructured 2D (EDU2D)

EDU2D_Quadrature_Triangle

  edu2d_quadrature_triangle_v1.f90

This program illustrates the use of a Gaussian quadrature, which integrates a function over a triangle. An example function, x^5 + y^5 -32*x^2*y^3 + x^3*y^2, is integrated over a triangle defined by the vertices: (0.1,0.3), (1.3,-0.2), (-0.9,1.7), and the numerical value obtained by a 7-point quadrature formula is compared with the exact value 0.53726127666666666668, up to 15 significant digits (double precision).



Educationally-Designed Unstructured 2D (EDU2D)

EDU2D-Airfoil-Spline (3 files)

  edu2d_airfoil_spline_v1.f90
  om6_wing_section_sharp.dat
  edu2d_airfoil_plot.m

Compile edu2d_airfoil_spline_v1.f90 and run it. This f90 program reads OM6 wing section data contained in the file 'om6_wing_section_sharp.dat', constructs a cubic spline (piecewise cubic polynomial interpolation), and writes out data files for plotting the airfoil by Matlab or Tecplot. I remember I used a cubic spline to represent a Joukowsky airfoil, so that boundary nodes can be adaptively moved along the surface to better represent numerical solutions [ IJNMF2002 ]


The OM6 section data are from [ AIAA J. Vol. 54, No. 9, September 2016 ], and the trailing edge is closed as described in the AIAA jounal paper.



Educationally-Designed Unstructured 2D (EDU2D)

EDU2D-Twod2Threed (5 files)

  edu2d_twod2threed_v2.f90
  edu2d_twod_rankine_half_body_v1.f90
  tria_rankine_input
  quad_rankine_input
  readme_v2.txt


This program reads a 2D mixed grid, and generates a 3D grid by extending the 2D grid to the third dimension. A hexahedral grid will be generated for a pure quadrilateral grid, and a prismatic grid will be generated for a pure triangular grid as below. Grids for a Rankine's half body are used as examples.



Educationally-Designed Unstructured 2D (EDU2D)

EDU2D-Tecplot2Grid (5 files)

  edu2d_tecplot2grid_v1.f90
  sample_input
  tria_rankine_grid_tecplot.dat
  tria_rankine_grid_tecplot.data
  readme_v1.txt


This program reads a Tecplot data file, generates a boundary information, and writes .grid and .bcmap file for edu2d solvers. Tecplot file does not contain boundary information, and so the boundary information is genearted inside the program. Not an elegant code, but it may be useful to learn how to identify boundary nodes for a given grid, which is not really trivial. A grid for a Rankine's half body is used as an example.



3D Grid Generation for Sphere:  

sphere_grid_v05.f90


A 3D grid is generated for a sphere and written in the UGRID format. A tetrahedarl, prismatic, or mixed grid can be generated with a specified y-plus value. Learn how a 3D unstructured grid is generated over a sphere.




3D Grid Generation for Hemisphere:  

hemisphere_grid_v09.f90


A 3D grid is generated for a hemisphere (with two configurations) and written in the UGRID format. A tetrahedarl, prismatic, or mixed grid can be generated with a specified y-plus value. Learn how a 3D unstructured grid is generated over a hemisphere.




3D Tetrahedral Grid Generation:  

tetgrid_cube_v3.f90 (regular grid)
tetgrid_cube_ptb_v5.f90 (irregular grid with nodal perturbation)

Tetrahedral grid in a unit cube is generated and written in the UGRID format. Learn how the hexahedron can be divided into six tetrahedra. Mdify the code to divide the hexahedron into five tetrahedra.

11-24-15: Added an irregular grid version (tetgrid_cube_ptb_v5.f90).

Educationally-Designed Unstructured 2D (EDU2D)

Special triangular grid generator: EDU2D-50yen-Tria-Grids

  edu2d_50yen_tria_grids_v0.f90
  edu2d_50yen_tria_grids_tecplot_v0.lay

This program can be useful for generating triangular grids for use in CFD or just learning how to manipulate unstructured-grid (invert, stretch, desired grid-spacing, perturb, diagonal swapping). It generates a simple triangulation of a sector, and then use it to generate 6 different triangular grids in a disk or a 50-yen domain. A tecplot style file is included to generate a picture such as the one below after the code is run.




Educationally-Designed Unstructured 2D (EDU2D)

EDU2D-AdvDiff (16 files)

  edu2d_advdiff_basic_package_v0.f90
  edu2d_advdiff_solver_nc_v1.f90
  edu2d_advdiff_main_v1.f90
  edu2d_advdiff_lsq_grad_nc_v0.f90
  edu2d_ddt3_v0.f90
  generate_grids_for_edu2d_steady_v0.f90
  generate_grids_for_edu2d_unsteady_v0.f90
  Makefile

  readme_v0.txt
  steady_screen_out.txt
  unsteady_screen_out.txt

  steady_plot_solution.lay
  steady_plot_boundary_solution.lay
  unsteady_initial_solution.lay
  unsteady_final_solution.lay
  unsteady_initial_final_solutions.lay


This is an unstructured advection-diffusion solver.

- Node-centered finite-volume discretization
- Fully unstructured grids (triangles, quads, or mixed)
- Least-squares gradients (linear or quadratic LSQ)
- Diffusion is discretized as a hyperbolic system for accurate gradient prediction
- Implicit time-stepping with BDF2
- Jacobian-Free Newton-Krylov with GCR
- Defect-correction implicit solver as a variable preconditioner
- Designed to solve both steady and unsteady problems.


Unseady advection diffusion problem.


Highly accurate normal gradient prediction on a boundary (steady).


Educationally-Designed Unstructured 2D (EDU2D)

EDU2D-Euler-Steady (9 files)

  edu2d_euler_steady_basic_package_v0.f90
  edu2d_euler_steady_solver_v0.f90
  edu2d_euler_steady_main_v0.f90
  edu2d_euler_jacobian_v1.f90
  edu2d_euler_linear_solve_v0.f90
  twod_bump_irregular_grid_v0.f90
  Makefile
  readme_v0.txt
  bump_screen_out.txt


Want to learn how to write an implicit unstructured CFD code? Grab this code, look inside to see how it is written, get good understanding, and then write your own. This code computes a steady flow over a bump with the Roe flux by two solution methods: an explicit 2-stage Runge-Kutta scheme and an implicit (defect correction) method with the exact Jacobian for a 1st-order scheme, on irregular triangular grids. A grid generation code is included for a bump problem.

- Node-centered finite-volume discretization
- Fully unstructured grids (triangles, quads, or mixed)
- Least-squares gradients (linear or quadratic LSQ)
- Defect-correction implicit solver (Newton's method for 1st-order)
- Both explicit and implicit solvers are implemented






Educationally-Designed Unstructured 2D (EDU2D)

EDU2D-Euler-Unsteady (8 files)

  edu2d_euler_unsteady_basic_package_v0.f90
  edu2d_euler_unsteady_solver_v0.f90
  edu2d_euler_unsteady_main_v0.f90
  twod_rectangular_grid_v0.f90
  Makefile
  readme_v0.txt
  project_quad_screen_out.txt
  project_tria_screen_out.txt


Want to learn how to write an unstructured CFD code? Grab this code, look inside to see how it is written, get good understanding, and then write your own. This code has Roe and Rotated-RHLL fluxes, Van Albada limiter, and a 2-stage Runge-Kutta time-stepping for solving a shock diffraction problem. It works for quadrilateral grids, triangular grids, and mixed grids also. It is set up to solve a shock diffraction problem. You can easily modify it for solving other problems. A grid generation code is included.

- Node-centered finite-volume discretization
- Fully unstructured grids (triangles, quads, or mixed)
- Least-squares gradients (linear or quadratic LSQ)
- Limiter is applied at edges
- Explicit time-stepping with 2-stage Runge-Kutta scheme






Educationally-Designed Unstructured 2D (EDU2D)

EDU2D-Template: 2D Example Unstructured CFD Code (8 files)

  edu2d_basic_package_v0.f90
  edu2d_template_main_v0.f90
  edu2d_template_nc_v0.f90
  edu2d_template_cc_v0.f90
  generate_grids_for_edu2d_v0.f90
  Makefile
  readme_v0.txt
  any_name_screen_out.txt


This code reads an unstructured grid file, generate various grid data, go through some dummy CFD solvers, and then writes out Tecplot data files for viewing the solution and the grid. Solvers are dummies and so do not solve anything, but you'll see how a node/cell-centered finite-volume solver can be implemented, e.g., loop over edges, nodes, elements, and computation of least-squares gradients, etc. You can generate your own unstructured CFD code by replacing the dummy solver by your own solver. A grid generation file is included, which can be used to generate grid files you need to run the EDU2D-Template code.

Flow over Rankine's half body (generate grids and exact solutions):

  twod_rankine_half_body_grid_v0.f90
A flow over Rankine's half body is an interesting exact solution to the potential flow equations: Cauchy-Riemann systems, Lapace equations, the potential equation, Euler equations. You can download, compile, and run this code to generate a (quadrilateral or triangular) grid, with the exact solution computed at each grid point (stored in Tecplot files). At some point, I'd like to use this solution for verification. Shown below is the exact x-velocity contours on an irregular triangular grid .





MPI program for Jacobi Iteration Ver.1:   jacobi_mpi_v1.f90
This is an MPI program for the Jacobi iteration, solving the finite-difference discretization of the Laplace equation in a square domain. It can be useful for learning how to write an MPI program.

Machine zero:   machinezero_v2.f90
What is machine zero? Maybe, it is a non-zero number which cannot be recognized by a machine. So, it depends on machines. Download, compile, and run it to find out the `zero' in your machine.

Generalized. Now it finds `machine zero' for a given value. (11-29-12)

Automatic Differentiation Ver.1:   ad_driver_v1.f90   module_ddt3.f90
Automatic differentiation is nice. It is sometimes used to compute flux Jacobians for implicit formulaitons in CFD. To those who are curious about it, here is a set of files with which you can experience and learn how it works. Simply compile the file "ad_driver.f90" and run it (the other file will be automatically included in the program) to get a feeling of automatic differentiation.

Numerical Fluxes (3D Euler) Ver.3: threed_euler_fluxes_v3.f90
Here are 3D Euler numerical fluxes. Download and take a look. Learn how the standard Roe and a very robust rotated-hybrid fluxes are implemented for the 3D Euler equations, and also how the Roe flux can be implemented without tangent vectors. Included are Roe with/without tangent vectors, and Rotated-RHLL fluxes. I'll be extremely happy if you kindly report bugs. Thank you, arigatou!

Bugs fixed for the RHLL flux, more comments added, Roe without tangent vectors added. (04-30-12).

3D Grid Generation for Hemishpere Cylinder:
The code is available at NASA's Turbulence Modeling Resource (TMR) website.

Flow around a Karman-Trefftz airfoil (generate grids and exact solutions):
  vkt_airfoil_v3.f90 ,   vkt_airfoil_v1_display.m
Karman-Trefftz airfoil is an intriguing airfoil for which a complete set of exact solutions can be computed. This solution has been used by many people to verify the accuracy of their inviscid code (incompressible limit). If you have not, use this code now to generate a (quadrilateral or triangular) grid, run your code on it, compute the error, and verify the accuracy of your code. Me? I have used it for my third-order multigrid Cauchy-Riemann solver (IJNMF 2004).

Updated (10-13-10). Formula corrected. It works now for cambered airfoils.

1D Hyperbolic Diffusion Scheme:   oned_upwind_diffusion.f90
Believe or not, the diffusion equation is solved by an upwind scheme. The idea is to integrate an equivalent hyperbolic system toward a steady state. This way, we can advance in time with a large O(h) time step (not O(h^2)), and compute the solution gradient with the equal order of accuracy. Compare with a common scheme (Galerkin) for 512 nodes to see how fast the upwind scheme can be.

Reference: JCP2007 | Preprint

3D Mixed Grid Generation:   mixgrid_cube_v3.f90
Mixed (tetrahedral-prismatic) grid in a unit cube is generated and written in the UGRID format. Learn how a typical viscous-type grid can be generated. You may want to modify the code to apply stretching to the prismatic layer for a smooth transition to the isotropic tetrahedral region.

3D Prismatic Grid Generation:   przgrid_cube_v3.f90
Prismatic grid in a unit cube is generated and written in the UGRID format. It's just a cubic domain, but may be useful in learning unstructured grid data.

3D Hex Grid Generation (Unstructured Format):   hexgrid_cube_v3.f90
Hexahedral grid in a unit cube is generated and written as unstructured/finite-element data in a UGRID 3D unstructured grid file. It may be useful for those who want to learn a typical data structure of unstructured grids.

Genetic Algorithm, Ver.1:   ga_v1.f90
This is a simple genetic algorithm program for finding minimizers of a function. I wrote this in 1999 when I was interested in aerodynamic optimization problems. I think that genetic algorithm is a very interesting optimization algorithm.

1D Euler Code Ver.1:   oned_euler_v1.f90 ,   oned_euler_plot_v1.m
Here is a 1D Euler code (1D shock tube code) for solving Sod's shock tube problem, using Roe's Approximate Riemann solver, minmod limiter, and 2-stage Runge-Kutta time-stepping. Learn how a second-order non-oscillatory Euler code is written, or just run it to see how it is capable of computing discontinuous solutions. Incorporate various flux subroutines given below to explore other methods.

Minor bugs fixed (12-29-10).

Numerical Fluxes (1D Euler) Ver.5:   oned_euler_fluxes_v5.f90
All numerical fluxes are functions of the left and right states (and possibly dt and dx). Learn how those famous fluxes can be implemented, or just use them to see how they work for various shock-tube problems. Included are Lax-Friedrichs, Richtmyer, MacCormack, Steger-Warming, Van Leer, AUSM, Zha-Bilgen, Godunov, Osher, Roe, Rusanov, HLL, HLLL, AUFS flux functions.

Bug fixed for Godunov fluxes.

Numerical Fluxes (2D Euler) Ver.2:   twod_euler_fluxes_v2.f90
Want 2D fluxes also? Here you are. Download and take a look. Learn how the standard Roe and a very robust rotated-hybrid fluxes are implemented. Included are Roe and Rotated-RHLL fluxes (Riemann Solvers).

Bugs fixed for the Rotated-RHLL flux. It works very well now.

Ringleb's Flow (generate grids and exact solutions):   ringleb_v1.f90
Ringleb's flow is a famous exact solution of the compressible Euler equations with a smooth transition from subsonic to supersonic without any shock waves. This solution has been used actually by many people to verify the accuracy of their Euler code. If you have not, use this code now to generate a (quadrilateral or triangular) grid, run your Euler code on it, compute the error, and verify the accuracy of your Euler code.

Minor bugs fixed (12-29-10).

Blasius solutions for a flow over a flat-plate (compute the exact solution) :
  blasius_v1.f90 ,   blasius_plot_v1.m
Exact solution for a flow over a flat-plate. This solution has been used by many people to verify the accuracy of their Navier-Stokes code. If you have not, use this code now to compute the exact solution at any point on your grid, compute the error, and verify the accuracy of your code.

Minor bugs fixed (12-29-10).

Viscous shock structure NS-solution (compute the exact solution) :
  ns_shock_structure_v1.f90 ,   ns_shock_structure_plot_v1.m
Exact solution for a shock wave internal structure to the 1D Navier-Stokes equations. This solution has been used by some people to verify the accuracy of their 1D Navier-Stokes code. If you have not, use this code now to generate a data file for the exact solution, use it as an initial solution for your code, converge to a steady state, and verify the accuracy of your code.

Minor bugs fixed (12-29-10).