Favorite Codes and Subroutines
EducationallyDesigned Unstructured 2D (EDU2D) EDU2DAdvDiff (16 files) edu2d_advdiff_basic_package_v0.f90 edu2d_advdiff_solver_nc_v2.f90 < Corrected 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
 
Kahan's sum to minimize the roundoff error Avoiding roundoff errors: kahan_sum_example.f90
 
EducationallyDesigned Unstructured 2D (EDU2D) EDU2D_MMS_TWODNS edu2d_mms_twodns_v2.f90
 
HCH Grid Generation Code Package:
 
EducationallyDesigned Unstructured 3D (EDU3D) EDU3D Euler: Serious bugs reported. Currently under investigation.  
Solver Package: HANIM and GCR:
Serious bugs reported. Currently under investigation. 
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.
112415: Added an irregular grid version (tetgrid_cube_ptb_v5.f90). 100617: Added a subroutine that computes edgebased data used in the edgebased discretization (tetgrid_cube_ptb_v6p1.f90). 
Diffusion equation is solved by 1st/2nd/3rdorder upwind schemes on irregularlyspaced grids.
The idea is to integrate an equivalent hyperbolic system toward
a steady state. This way, we can advance in pseudo time with a large
O(h) time step (not O(h^2)), and compute the solution gradient
with the equal order of accuracy on irregular grids. Compare with a conventional
diffusion scheme to see how accurate the gradient can be and
how fast it converges.
Note: This thirdorder scheme does not require computation/storage of second derivatives. See Ref.[2]. References: [1] JCP2014  Preprint  1st/2nd/3rdorder FV upwind schemes for diffusion [2] JCP2017  Preprint  Thirdorder accuracy without second derivatives. 
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 7point quadrature formula is compared with the exact value 0.53726127666666666668,
up to 15 significant digits (double precision).

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. 
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.

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.

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 yplus value.
Learn how a 3D unstructured grid is generated over a sphere.

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 yplus value.
Learn how a 3D unstructured grid is generated over a hemisphere.

This program can be useful for generating triangular grids for use in CFD or just learning how to manipulate
unstructuredgrid (invert, stretch, desired gridspacing, 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 50yen domain. A tecplot style file is
included to generate a picture such as the one below after the code is run.

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 2stage RungeKutta scheme and an implicit (defect correction) method with
the exact Jacobian for a 1storder scheme, on irregular triangular grids. A grid generation code is included for a
bump problem.
 Nodecentered finitevolume discretization  Fully unstructured grids (triangles, quads, or mixed)  Leastsquares gradients (linear or quadratic LSQ)  Defectcorrection implicit solver (Newton's method for 1storder)  Both explicit and implicit solvers are implemented 
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 RotatedRHLL fluxes, Van Albada limiter, and a 2stage RungeKutta
timestepping 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.
 Nodecentered finitevolume discretization  Fully unstructured grids (triangles, quads, or mixed)  Leastsquares gradients (linear or quadratic LSQ)  Limiter is applied at edges  Explicit timestepping with 2stage RungeKutta scheme 
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/cellcentered finitevolume solver can be implemented, e.g., loop over edges, nodes, elements, and computation of leastsquares 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 EDU2DTemplate code. 
A flow over Rankine's half body is an interesting exact solution to the potential flow
equations: CauchyRiemann 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 xvelocity contours on
an irregular triangular grid .

This is an MPI program for the Jacobi iteration, solving the finitedifference
discretization of the Laplace equation in a square domain. It can be useful
for learning how to write an MPI program.

What is machine zero? Maybe, it is a nonzero 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. (112912) 
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.

Here are 3D Euler numerical fluxes. Download and take a look. Learn
how the standard Roe and a very robust rotatedhybrid 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 RotatedRHLL 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. (043012). 
The code is available at NASA's
Turbulence Modeling Resource (TMR) website.

KarmanTrefftz 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 thirdorder multigrid CauchyRiemann
solver (IJNMF
2004).
Updated (101310). Formula corrected. It works now for cambered airfoils. 
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 
Mixed (tetrahedralprismatic) grid in a unit cube is generated and written in the UGRID format. Learn how a typical viscoustype 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. 
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. 
Hexahedral grid in a unit cube is generated and written as unstructured/finiteelement 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. 
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. 
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 2stage RungeKutta
timestepping. Learn how a secondorder nonoscillatory 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 (122910). 
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 shocktube problems.
Included are LaxFriedrichs, Richtmyer, MacCormack,
StegerWarming, Van Leer, AUSM, ZhaBilgen, Godunov, Osher, Roe, Rusanov,
HLL, HLLL, AUFS flux functions. Bug fixed for Godunov fluxes. 
Want 2D fluxes also? Here you are. Download and take a look. Learn
how the standard Roe and a very robust rotatedhybrid fluxes are implemented.
Included are Roe and RotatedRHLL fluxes (Riemann Solvers).
Bugs fixed for the RotatedRHLL flux. It works very well now. 
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 (122910). 
Exact solution for a flow over a flatplate.
This solution has been used by many people to verify the accuracy
of their NavierStokes 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 (122910). 
Exact solution for a shock wave internal structure to the 1D NavierStokes
equations. This solution has been used by some people to verify the accuracy
of their 1D NavierStokes 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 (122910). 