How to compute flux Jacobians:
edu3d_exact_flux_jacobian_v0.f90
module_ddt5.f90
[To compile: gfortran edu3d_exact_flux_jacobian_v0.f90]
This code shows two ways to compute flux Jacobians:
automatic differentiation and handcoded subroutines for
the Roe flux (inviscid) and the alphadamping flux (viscous) in
three dimensions. This serves also as verification of the handcoded
Jacobians (which are complicated). See how the Roe flux is exactly differentiated.
It may be useful also for learning how the viscous flux is computed.
Uploaded. (040719)


Machine zero:
machinezero_v3.f90
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.
Comments corrected. (092318)
Generalized. Now it finds `machine zero' for a given value. (112912)


EducationallyDesigned Unstructured 1D (EDU1D)
First/Second/ThirdOrder Diffusion Schemes in 1D:
edu1d_oned_first_order_diffusion_v1.f90
edu1d_oned_second_order_diffusion_v1.f90
edu1d_oned_third_order_diffusion_v1.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.


Solver Package: GCR solver module:
gcr_module_v3.f90
This module is a GCR (Generalized Conjugate Residual) linear solver, which
can be used to construct a JacobianFree NewtonKrylov solver (JFNK).
The code illustrates how a JFNK solver can be implemented.


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
This is an unstructured advectiondiffusion solver.
 Nodecentered finitevolume discretization
 Fully unstructured grids (triangles, quads, or mixed)
 Leastsquares gradients (linear or quadratic LSQ)
 Diffusion is discretized as a hyperbolic system for accurate gradient prediction
 Implicit timestepping with BDF2
 JacobianFree NewtonKrylov with GCR
 Defectcorrection 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).


Kahan's sum to minimize the roundoff error
Avoiding roundoff errors:
kahan_sum_example.f90
This code shows how the roundoff error can be minimized
in computing a sum of small values, by using Kahan's roundoffcorrected summation.


EducationallyDesigned Unstructured 2D (EDU2D)
EDU2D_MMS_TWODNS
edu2d_mms_twodns_v2.f90
This code shows how to compute the source terms in the method of
manufactured solutions (MMS) for the 2D NavierStokes equations. Any function can be made
an exact solution to the 2D NavierStokes equations with suitable source terms. That is,
any function v(x,y) is an exact solution to the following equation:
MMS has been widely used for verifying the order of accuracy of a CFD code.
Grab this code, and learn how it can be implemented. I always use MMS to verify the
implementation. It is a very useful tool for debugging.


HCH Grid Generation Code Package:
This package contains two f90 codes.
(1)Structured/unstructured grid generation code for
a hemispherecylinderhemisphere body. It generates
hexprism structured, prismatic, tetrahedral, mixed prismtetra, and
mixed prismhex grids.
The shape of the body can be defined by a shape function.
Four different functions are implemented in the code: a suboff model,
a suboff with a hemisphere front, a sine function, and an ellipsoid.
(2)A coarsening program. It generates a series of coarse grids
from the grid generated in (1).
12172017: Several bugs fixed.


EducationallyDesigned Unstructured 3D (EDU3D)
EDU3D Euler:
Serious bugs reported. Currently under investigation.

EducationallyDesigned Unstructured 3D (EDU3D)
3D Tetrahedral Grid Generation:
edu3d_tetgrid_cube_v3.f90 (regular grid)
edu3d_tetgrid_cube_ptb_v6p1.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.
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).


EducationallyDesigned 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 7point quadrature formula is compared with the exact value 0.53726127666666666668,
up to 15 significant digits (double precision).


EducationallyDesigned Unstructured 2D (EDU2D)
EDU2DAirfoilSpline (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.


EducationallyDesigned Unstructured 2D (EDU2D)
EDU2DTwod2Threed (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.


EducationallyDesigned Unstructured 2D (EDU2D)
EDU2DTecplot2Grid (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 yplus 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 yplus value.
Learn how a 3D unstructured grid is generated over a hemisphere.


EducationallyDesigned Unstructured 2D (EDU2D)
Special triangular grid generator: EDU2D50yenTriaGrids
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
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.


EducationallyDesigned Unstructured 2D (EDU2D)
EDU2DEulerSteady (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 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


EducationallyDesigned Unstructured 2D (EDU2D)
EDU2DEulerUnsteady (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 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


EducationallyDesigned Unstructured 2D (EDU2D)
EDU2DTemplate: 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/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.


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


MPI program for Jacobi Iteration Ver.1:
jacobi_mpi_v1.f90
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.


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


3D Grid Generation for Hemishpere Cylinder:

Flow around a KarmanTrefftz airfoil (generate grids and exact solutions):
vkt_airfoil_v3.f90
,
vkt_airfoil_v1_display.m
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.


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_v4.f90
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. (08222018: A bug fixed
about the boundary grid information.)


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


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


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


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 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 (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 (122910).


Blasius solutions for a flow over a flatplate (compute the exact solution) :
blasius_v1.f90
,
blasius_plot_v1.m
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).


Viscous shock structure NSsolution (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 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).

