Finally, we get to the end! The complete code for the main application is shown in the following listing.
Listing 4.10 Complete code of the main program of the tsunami simulator
program tsunami
use iso_fortran_env, only: int32, real32 ❶
use mod_diff, only: diff => diff_centered ❷
use mod_initial, only: set_gaussian ❸
implicit none
integer(int32) :: n
integer(int32), parameter :: grid_size = 100 ❹
integer(int32), parameter :: num_time_steps = 5000 ❹
real(real32), parameter :: dt = 0.02, dx = 1, g = 9.8 ❺
real(real32), parameter :: hmean = 10
real(real32) :: h(grid_size), u(grid_size) ❻
integer(int32), parameter :: icenter = 25 ❼
real(real32), parameter :: decay = 0.02 ❼
if (grid_size <= 0) stop 'grid_size must be > 0' ❽
if (dt <= 0) stop 'time step dt must be > 0' ❽
if (dx <= 0) stop 'grid spacing dx must be > 0' ❽
call set_gaussian(h, icenter, decay) ❾
u = 0 ❾
print *, 0, h
time_loop: do n = 1, num_time_steps ❿
u = u - (u * diff(u) + g * diff(h)) / dx * dt ⓫
h = h - diff(u * (hmean + h)) / dx * dt ⓫
print *, n, h
end do time_loop
end program tsunami❶ Imports type kind parameters for numeric data
❷ Imports the diff_centered function from the mod_diff module and renames it as diff
❸ Imports the set_gaussian procedure from the mod_initial module
❹ Sets grid size and number of time steps
❺ Sets time step, grid spacing, and gravitational acceleration
❻ Declares arrays to use in the simulation
❼ Parameters to use to set the initial condition for water height
❽ Checks input parameter values
❾ Initializes water height and velocity arrays
❿ Loops for num_time_steps time steps
⓫ Computes the water velocity in the next step and updates
⓫ Computes the water height in the next step and updates
There are quite a few changes in this program relative to our previous version.
First, we’re now using a new function (diff_centered) for the finite difference calculation, which allows us to account for changes coming from either direction. On import, we rename this function to just diff so that we don’t have to change the existing terms from the previous version of the simulator. The other imports from modules are also new and should be familiar by now, if you’ve worked through this chapter in order.
Second, our main simulation loop now has two equations instead of one. We’re now solving for both water velocity (u) and water height (h). These two equations are coupled: water velocity appears on the right side of the equation for water height, and vice versa. This coupling is exactly what leads to a more realistic fluid flow. However, more realistic physics requires higher resolution in time, and this is reflected in our time step dt now being smaller than before. For an analogy, think about watching an extremely slow-paced video in which not much is going on. You may not even notice if the video is playing at a low frame rate. However, if you’re watching a fast-paced video of a car chase or a basketball game, you’ll easily notice any lag or drop in frame rate. The same goes for fluid dynamics–a more realistic simulation (more physics terms) and a finer computational grid (level of detail) require a shorter time step, which translates to a higher frame rate.
Congratulations, we made it! If you now compile the simulator, run it, and plot the results, you’ll be able to see the flow similar to that in figure 4.4. To account for the shorter time step, I increased the number of time steps (now 5,000) to simulate for a long enough period of time. I encourage you to play with a few of the simulation parameters, such as the grid size, number of time steps, and the initial shape and perturbation of the water.