Tsunami Simulator

Tsunami Simulator#

The following program simulates a tsunami water wave.
Adapted from: “Modern Fortran: Building Efficient Parallel Applications” by Milan Curcic (Manning)

The tsunami water wave height and speed will be calculated using the advection equation.

From wikipedia:

Note

In the field of physics, engineering, and earth sciences, advection is the transport of a substance or quantity by bulk motion of a fluid.

The Advection Equation#

\[ \Large \frac{\delta u}{\delta t} + c\frac{\delta u}{\delta x} = 0 \]

Program to Predict the Linear Advection of an Object, with Periodic Boundary Conditions#

program tsunami

    implicit none
    integer :: i, n
    
    integer, parameter :: grid_size = 100
    integer, parameter :: num_time_steps = 100

    real, parameter :: dt = 1 ! time step [s]
    real, parameter :: dx = 1 ! grid spacing [m]
    real, parameter :: c  = 1 ! phase speed [m/s]

    real :: h(grid_size), dh(grid_size)
    
    integer, parameter :: icenter = 25
    real, parameter    :: decay = 0.02
    
    logical :: file_exists
    
    open(9, file = 'tsunami_simulator_data.txt')
    
    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'
    if ( c <= 0 )  stop 'background flow speed c must be > 0'
    
    do concurrent (i = 1:grid_size)
        h(i) = exp(-decay * (i - icenter)**2)
    end do
    
    !print *, 0, h
    write (9, *) 0, h
    close(9)
    
    time_loop: do n=1,num_time_steps    
        
        dh(1) = h(1) - h(grid_size)
        
        do concurrent (i = 2:grid_size)
            dh(i) = h(i) - h(i-1)       
        end do
        
        do concurrent (i = 1:grid_size)
            h(i)=h(i)-c*dh(i) / dx * dt
        end do
        
        !print *, n, h
        inquire(file = 'tsunami_simulator_data.txt', exist = file_exists)
        if (file_exists) then
            open(9, file = 'tsunami_simulator_data.txt', status = 'old', position = 'append', action = 'write')
        else
            open(9, file = 'tsunami_simulator_data.txt', status = "new", action = 'write')
        end if
        
        write (9, *) n, h

    end do time_loop

    close (9)
    
end program tsunami

Program Explanation#

The Fortran program tsunami simulates the propagation of a wave in a one-dimensional grid using a basic numerical scheme akin to the finite difference method. The purpose of the program is to approximate how a disturbance, initialized in the middle of the domain, evolves over time given the specified parameters.

Key Elements of the Program#

Variable Declarations and Constants:
- grid_size: The number of grid points in the simulation domain.
- num_time_steps: The number of time steps for which the simulation will run.
- dt: Time step in seconds.
- dx: Spatial grid spacing in meters.
- c: Phase speed of the wave in meters per second.
- h: Height of the wave at each grid point.
- dh: Temporal change in height at each grid point between consecutive steps.
Initialization:
- An initial wave profile is set using a Gaussian-like exponential decay formula centered at icenter, influencing how sharply the wave amplitude decreases away from the center.
File Handling:
- The program opens a file named ‘tsunami_simulator_data.txt’ for output and checks for its existence before each write operation to handle potential errors or file state changes during execution.
Numerical Scheme:
- The simulation uses a simple difference scheme to approximate the temporal evolution of the wave:
  - Boundary Condition: A cyclic boundary condition is assumed (dh(1) = h(1) - h(grid_size)), which allows the wave to wrap around the domain, simulating an infinite or periodic space.
  - Interior Points: The change in height (dh) at each internal grid point is computed based on the difference with its immediate left neighbor.
  - Update Step: The wave height at each grid point is updated using the computed differences, scaled by the wave speed, grid spacing, and time step.
Looping Structure:
- The main loop (time_loop) iterates over the number of time steps, updating the wave’s profile at each step and writing the results to the output file.
- An inquire statement checks if the output file still exists before attempting to append new data, which helps in managing output data integrity.
Output:
- For each time step, the wave heights across the grid are written to the file. This output can be used for further analysis or visualization of the wave’s propagation.

Program Execution Details#

Initialization of h:
- The wave height at each grid point is initialized to create a localized disturbance. The exponential decay formula makes the initial wave amplitude highest at the center and decreases moving away from the center.
Time Evolution:
- At each time step, the program computes the spatial gradient of the wave heights, then adjusts the heights based on this gradient, effectively simulating the wave propagation.
File Operations:
- Before writing to the file, the program ensures it’s in the correct state (open and ready for appending), handling the possibility that the file might not be available or already open from previous operations.

Practical Applications#

This program is a simplistic representation useful for educational purposes in teaching numerical methods for solving partial differential equations, specifically hyperbolic equations like the wave equation. It can also serve as a basic model in physical sciences and engineering for understanding wave dynamics in a controlled setting.

Program Compilation and Execution#

The above program is compiled and run using Fortran Package Manager (fpm):

import os
root_dir = os.getcwd()

code_dir = root_dir + "/" + "Fortran_Code/Tsunami_Simulator"

os.chdir(code_dir)

build_status = os.system("fpm build 2>/dev/null")

exec_status = os.system("fpm run 2>/dev/null")

The above Fortran code writes the calculations to a file called tsunami_simulator_data.txt

import pandas as pd
output_filename = 'tsunami_simulator_data.txt'
data_file = code_dir + "/" + output_filename
table = pd.read_fwf(data_file, header=None)
table

	0	1	2	3	4	5	6	7	8	9	...	91	92	93	94	95	96	97	98	99	100
0	0	9.929509e-06	0.000025	0.000063	0.000148	0.000335	0.000732	0.001534	0.003089	0.005976	...	1.459711e-38	1.021038e-39	6.862018e-41	4.430906e-42	2.746545e-43	1.681558e-44	1.401298e-45	0.000000e+00	0.000000e+00	0.000000e+00
1	1	0.000000e+00	0.000010	0.000025	0.000063	0.000148	0.000335	0.000732	0.001534	0.003089	...	2.005009e-37	1.459711e-38	1.021038e-39	6.862018e-41	4.430906e-42	2.746545e-43	1.681558e-44	1.401298e-45	0.000000e+00	0.000000e+00
2	2	0.000000e+00	0.000000	0.000010	0.000025	0.000063	0.000148	0.000335	0.000732	0.001534	...	2.646043e-36	2.005009e-37	1.459711e-38	1.021038e-39	6.862018e-41	4.430906e-42	2.746545e-43	1.681558e-44	1.401298e-45	0.000000e+00
3	3	0.000000e+00	0.000000	0.000000	0.000010	0.000025	0.000063	0.000148	0.000335	0.000732	...	3.355098e-35	2.646043e-36	2.005009e-37	1.459711e-38	1.021038e-39	6.862018e-41	4.430906e-42	2.746545e-43	1.681558e-44	1.401298e-45
4	4	1.401298e-45	0.000000	0.000000	0.000000	0.000010	0.000025	0.000063	0.000148	0.000335	...	4.087346e-34	3.355098e-35	2.646043e-36	2.005009e-37	1.459711e-38	1.021038e-39	6.862018e-41	4.430906e-42	2.746545e-43	1.681558e-44
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
96	96	3.354626e-04	0.000732	0.001534	0.003089	0.005976	0.011109	0.019841	0.034047	0.056135	...	2.746545e-43	1.681558e-44	1.401298e-45	0.000000e+00	0.000000e+00	0.000000e+00	9.929508e-06	2.541935e-05	6.252155e-05	1.477484e-04
97	97	1.477484e-04	0.000335	0.000732	0.001534	0.003089	0.005976	0.011109	0.019841	0.034047	...	4.430906e-42	2.746545e-43	1.681558e-44	1.401298e-45	0.000000e+00	0.000000e+00	0.000000e+00	9.929508e-06	2.541935e-05	6.252155e-05
98	98	6.252155e-05	0.000148	0.000335	0.000732	0.001534	0.003089	0.005976	0.011109	0.019841	...	6.862018e-41	4.430906e-42	2.746545e-43	1.681558e-44	1.401298e-45	0.000000e+00	0.000000e+00	0.000000e+00	9.929508e-06	2.541935e-05
99	99	2.541935e-05	0.000063	0.000148	0.000335	0.000732	0.001534	0.003089	0.005976	0.011109	...	1.021038e-39	6.862018e-41	4.430906e-42	2.746545e-43	1.681558e-44	1.401298e-45	0.000000e+00	0.000000e+00	0.000000e+00	9.929508e-06
100	0	9.929508e-06	0.000025	0.000063	0.000148	0.000335	0.000732	0.001534	0.003089	0.005976	...	1.459711e-38	1.021038e-39	6.862018e-41	4.430906e-42	2.746545e-43	1.681558e-44	1.401298e-45	0.000000e+00	0.000000e+00	0.000000e+00

101 rows × 101 columns