!********************************************************************************
!* Educationally-Designed Unstructured 2D (EDU2D) Code
!*
!*
!* --- This is EDU2D-Euler-RK2 ---
!*
!*
!* EDU2D-Euler-RK2: An Euler code with
!*
!* - Node-centered finite-volume discretization
!* - 2-stage Runge-Kutta explicit time-stepping scheme (RK2)
!*
!*
!*
!* specially set up for a shock-diffraction problem
!*
!* Wall
!* --------------------
!* Post-shock (inflow) | |
!* |->Shock | o: Corner node
!* | Mach=5.09 |
!* .......o |Outflow
!* Wall | |
!* | |
!* | |
!* --------------------
!* Outflow
!*
!* - Node-centered finite-volume method for unstructured grids (quad/tri/mixed)
!* - Roe flux with an entropy fix and Rotated-RHLL flux
!* - Gradient reconstruction by unweighted least-squares method
!* - Van Albada slope limiter to the primitive variable gradients
!* - 2-Stage Runge-Kutta global time-stepping towards the final time
!* - All quantities are nondimensionalized; velocity and pressure are
!* nondimensionalized based on the free stream speed of sound
!* (see "I do like CFD, VOL.1", whose PDF can be downloaded at cfdbooks.com).
!*
!* written by Dr. Katate Masatsuka (info[at]cfdbooks.com),
!*
!* the author of useful CFD books, "I do like CFD" (http://www.cfdbooks.com).
!*
!*
!* This is Version 0 (July 2015)
!*
!* ------------------------------------------------------------------------------
!* Files: There are 3 files.
!*
!* ------------------------------------------
!* - Main driver program file : This reads grid and BC files, and call dummy NC/CC programs.
!*
!* edu2d_euler_rk2_main.f90 (This file), which contains a main driver program
!* -- edu2d_euler_rk2 : Main driver code, which calls an Euler solver
!*
!* ------------------------------------------
!* - Basic EDU2D package file : Arranged for a 2D Euler code
!*
!* edu2d_basic_package_euler_rk2.f90, which contains the following modules.
!* -- edu2d_constants : Numerical values defined
!* -- edu2d_grid_data_type : Grid data types defined
!* -- edu2d_main_data : Main grid data and parameters declared
!* -- edu2d_grid_data : Read/construct/check grid data
!*
!* ------------------------------------------
!* - Euler solver file : This computes a solution to the shock diffraction problem.
!*
!* edu2d_euler_rk2_main.f90, which contains a 2D Euler solver with RK2.
!* -- edu2d_euler_rk2_solver : Node-centered Explicit Euler solver with RK2
!*
!* ------------------------------------------------------------------------------
!* Notes:
!*
!* The purpose of this code is to give a beginner an opportunity to learn how to
!* write an unstructured CFD code. Hence, the focus here is on the simplicity.
!* The code is not optimized for efficiency.
!*
!* This code is set up specifically for a shock diffraction problem.
!* It can be modified easily to solve other problems:
!*
!* 1. Define your own free stream Mach number, M_inf, at the beginning of the main.
!* 2. Modify the subroutine, initial_solution_shock_diffraction(), which is in this
!* file, to set up an appropriate initial condition for your problem.
!* 3. Delete the special treatment at the corner in euler_solver_main.f90
!* (Special treatment is done in two places in that file.)
!*
!* If the code is not simple enough to understand, please send questions to Hiro
!* at sunmasen(at)hotmail.com. I'll greatly appreciate it and revise the code.
!*
!* If the code helps you understand how to write your own code that is more
!* efficient and has more features, it'll have served its purpose.
!*
!* ------------------------------------------------------------------------------
!* Examples of additional features you might want to add.
!*
!* 1. Local time-stepping (faster convergence for steady problems)
!* 2. Implicit time-stepping (remove CFL restriction)
!* 3. More boundary conditions (periodic, symmetry, suction, etc.)
!* 4. Other reconstruction (Van Leer's kappa scheme)
!* 5. Other limiters (Venkat/Barth limiter,etc.)
!* 6. Other flux functions (HLL, LDFSS, AUSM, etc.)
!* 7. Local-preconditioning (low-Mach number accuracy and stiffness removal)
!* 8. CFL ramping (increase CFL gradually for a stable start-up)
!* 9. More output (convergence history, etc.)
!* 10. Parallelization (large-scale problems)
!* 11. Grid adaptation (h-refinement, steady or unsteady)
!* 12. Adjoint capability (aerodynamic design, adaptation, etc.)
!* 13. Moving grid (sliding mesh, overset grid, etc.)
!* 14. Multigrid (grid-independent convergence)
!* ------------------------------------------------------------------------------
!*
!* Katate Masatsuka, http://www.cfdbooks.com
!********************************************************************************
!********************************************************************************
!* Main program: Node-centered finite-volume Euler code
!*
!* This code computes an unsteady solution of the Euler equations.
!* It is set up to solve a shock diffraction problem.
!* So, it is not really (but almost) a general purpose code.
!*
!* Input -------------------------------------------------------
!*
!* project.grid = grid file containing boundary information
!* project.bcmap = file that specifies boundary conditions
!*
!* (Note: See the subroutine "read_grid", which is in this file,
!* for the format of these files.)
!*
!* Parameters to be specified inside the main program:
!*
!* M_inf = Upstream Mach number
!* gamma = Ratio of specific heats (1.4 for air)
!* CFL = CFL number (global time step)
!* t_final = Final time to stop the calculation
!* time_step_max = Max iterations (just a big enough number)
!* inviscid_flux = Inviscid flux selection (1 = Roe, 2 = Rotated-RHLL)
!* limiter_type = Limiter selection (1=Van Albada limiter, 0=No limiter)
!*
!*
!* Output ------------------------------------------------------
!*
!* "project_tecplot.dat" = Tecplot file containing the grid and the solution.
!*
!*
!* NOTE: The variables are nondimensionalized values (compressible scale),
!* rho=rho/rho_inf, u=u/a_inf, v=v/a_inf, rho=p/(rho_inf*a_inf^2)
!*
!* NOTE: Many variables are passed to subroutines via USE statement.
!* Each module contains specific data, and they are accessed by USE.
!*
!*
!* Katate Masatsuka, http://www.cfdbooks.com
!********************************************************************************
program edu2d_euler_rk2
use edu2d_constants , only : p2, one
use edu2d_my_allocation , only : my_alloc_p2_ptr, my_alloc_p2_matrix_ptr
use edu2d_my_main_data , only : nq, M_inf, gamma, time_step_max, t_final, &
CFL, inviscid_flux, limiter_type, &
gradient_type, gradient_weight, gradient_weight_p, &
node, nnodes
use edu2d_grid_data , only : read_grid, construct_grid_data, check_grid_data
use edu2d_euler_rk2_solver, only : euler_solver_main, compute_lsq_coeff_nc, check_lsq_coeff_nc
implicit none
!Inout data files
character(80) :: datafile_grid_in !Grid file
character(80) :: datafile_bcmap_in !Boundary condition file
!Output data file
character(80) :: datafile_tec !Tecplot file for viewing the result.
integer :: i
! Set file names
datafile_grid_in = 'project.grid'
datafile_bcmap_in = 'project.bcmap'
datafile_tec = 'project_tecplot.dat'
!--------------------------------------------------------------------------------
! Input Parameters
M_inf = 0.0_p2 ! Freestream Mach number to be set in the subroutine
! -> "initial_solution_shock_diffraction"
! (Specify M_inf here for other problems.)
gamma = 1.4_p2 ! Ratio of specific heats
CFL = 0.95_p2 ! CFL number
t_final = 0.18_p2 ! Final time to stop the calculation.
time_step_max = 5000 ! Max time steps (just a big enough number)
inviscid_flux = "rhll" ! = Rotated-RHLL , "roe" = Roe flux
limiter_type = "vanalbada" ! = Van Albada limiter, "none" = No limiter
nq = 4 ! The number of equtaions/variables in the target equtaion.
gradient_type = "linear" ! or "quadratic2 for a quadratic LSQ.
gradient_weight = "none" ! or "inverse_distance"
gradient_weight_p = one ! or any other real value
!--------------------------------------------------------------------------------
! Solve the Euler equations and write the output datafile.
!
! (1) Read grid files
call read_grid(datafile_grid_in, datafile_bcmap_in)
!Allocate arrays
do i = 1, nnodes
allocate( node(i)%u( nq ) )
allocate( node(i)%du( nq ) )
allocate( node(i)%w( nq ) )
allocate( node(i)%gradw(nq,2) ) !<- 2: x and y components.
allocate( node(i)%res( nq ) )
end do
! (2) Construct grid data
call construct_grid_data
! (3) Check the grid data (It is always good to check them before use!)
call check_grid_data
! (4) Prepare LSQ gradients
call compute_lsq_coeff_nc
call check_lsq_coeff_nc
! (5) Set initial solution for a shock diffraction problem
! (Re-write or replace it by your own subroutine for other problems.)
call initial_solution_shock_diffraction
! (6) Compute the solution (March in time to the final time)
call euler_solver_main
! (7) Write out the tecplot data file (Solutions at nodes)
call write_tecplot_file(datafile_tec)
!--------------------------------------------------------------------------------
write(*,*) "Successfully completed. Stop."
stop
contains
!********************************************************************************
!* Initial solution for the shock diffraction problem:
!*
!* NOTE: So, this is NOT a general purpose subroutine.
!* For other problems, specify M_inf in the main program, and
!* modify this subroutine to set up an appropriate initial solution.
!
! Shock Diffraction Problem:
!
! Wall
! --------------------
! Post-shock (inflow) | |
! (rho,u,v,p)_inf |->Shock (M_shock) | o: Corner node
! M_inf | |
! .......o Pre-shock |Outflow
! Wall | (rho0,u0,v0,p0) |
! | |
! | |
! --------------------
! Outflow
!
!********************************************************************************
subroutine initial_solution_shock_diffraction
use edu2d_constants , only : p2, zero, one, two
use edu2d_my_main_data , only : nnodes, node, gamma, M_inf, rho_inf, u_inf, v_inf, p_inf
use edu2d_euler_rk2_solver, only : w2u
implicit none
!Local variables
integer :: i
real(p2) :: M_shock, u_shock, rho0, u0, v0, p0
do i = 1, nnodes
! Pre-shock state: uniform state; no disturbance has reahced yet.
rho0 = one
u0 = zero
v0 = zero
p0 = one/gamma
! Incoming shock speed
M_shock = 5.09_p2
u_shock = M_shock * sqrt(gamma*p0/rho0)
! Post-shock state: These values will be used in the inflow boundary condition.
rho_inf = rho0 * (gamma + one)*M_shock**2/( (gamma - one)*M_shock**2 + two )
p_inf = p0 * ( two*gamma*M_shock**2 - (gamma - one) )/(gamma + one)
u_inf = (one - rho0/rho_inf)*u_shock
M_inf = u_inf / sqrt(gamma*p_inf/rho_inf)
v_inf = zero
! Set the initial solution: set the pre-shock state inside the domain.
node(i)%w = (/ rho0, u0, v0, p0 /)
node(i)%u = w2u(node(i)%w)
end do
end subroutine initial_solution_shock_diffraction
!********************************************************************************
!* Write a tecplot file: grid and solution
!********************************************************************************
subroutine write_tecplot_file(datafile_tec)
use edu2d_my_main_data, only : nnodes, node, elm, nelms, gamma
implicit none
character(80), intent(in) :: datafile_tec
integer :: i, k, os
!--------------------------------------------------------------------------------
open(unit=1, file=datafile_tec, status="unknown", iostat=os)
write(1,*) 'title = "grid"'
write(1,'(a80)') 'variables = "x","y","rho","u","v","p","Mach"'
write(1,*) 'zone n=',nnodes,' e =', nelms,' et=quadrilateral, f=fepoint'
!--------------------------------------------------------------------------------
! Nodal quantities: x, y, rho, u, v, p, Mach number
do i = 1, nnodes
write(1,*) node(i)%x, node(i)%y, (node(i)%w(k),k=1,4), &
sqrt( (node(i)%w(2)**2+node(i)%w(3)**2)/(gamma*node(i)%w(4)/node(i)%w(1)) )
end do
!--------------------------------------------------------------------------------
! Both quad and tria elements in quad format:
do i = 1, nelms
if (elm(i)%nvtx == 3) then
write(1,*) elm(i)%vtx(1), elm(i)%vtx(2), elm(i)%vtx(3), elm(i)%vtx(3)
elseif (elm(i)%nvtx == 4) then
write(1,*) elm(i)%vtx(1), elm(i)%vtx(2), elm(i)%vtx(3), elm(i)%vtx(4)
else
!Impossible
write(*,*) " Error in elm%vtx data... Stop..: elm(i)%nvtx=",elm(i)%nvtx
stop
endif
end do
!--------------------------------------------------------------------------------
close(1)
end subroutine write_tecplot_file
!********************************************************************************
end program edu2d_euler_rk2