Tsunami Simulator with Modules

Tsunami Simulator with Modules#

Adapted from: “Modern Fortran: Building Efficient Parallel Applications” by Milan Curcic (Manning)

Program to Calculate the Temperature due to a Cold Front Using a Loop#

In this section the Tsunami Simulator is rewritten using Fortran Modules.

The file mod_inital.f08 contains:#

module mod_initial
    
    private
    public :: set_gaussian

contains
    
    pure subroutine set_gaussian(x, icenter, decay)
        real, intent(in out) :: x(:)
        integer, intent(in)  :: icenter
        real, intent(in)     :: decay
        
        integer :: i
    
        do concurrent(i = 1:size(x))
            x(i) = exp(-decay * (i - icenter)**2)
        end do
    end subroutine set_gaussian
    
end module mod_initial

The file mod_diff.f08 contains:#

! %module magic necessary so that jupyter fortran kernel knows that 
! the code below is to be compiled as a module and not executable code.
!%module: mod_diff 
module mod_diff

    use iso_fortran_env, only: int32, real32
    implicit none
    
    private
    public :: diff_upwind, diff_centered
    
contains

    pure function diff_upwind(x) result(dx)
        real(real32), intent(in) :: x(:)
        real(real32)             :: dx(size(x))
        integer(int32)           :: im
        
        im       = size(x)
        dx(1)    = x(1) - x(im)
        dx(2:im) = x(2:im) - x(1:im - 1)
    end function diff_upwind
    
    pure function diff_centered(x) result(dx)
        real(real32), intent(in) :: x(:)
        real(real32)             :: dx(size(x))
        integer(int32)           :: im
        
        im         = size(x)
        dx(1)      = x(2) - x(im)
        dx(im)     = x(1) - x(im - 1)
        dx(2:im-1) = x(3:im) - x(1:im-2)
    end function diff_centered
    
end module mod_diff

Main Program#

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
    
    logical :: file_exists
    
    open(9, file = 'tsunami_simulator_with_modules_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'
    
    call set_gaussian(h, icenter, decay)
    u = 0
    
    !print *, 0, h
    write (9, *) 0, h
    close(9)
    
    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
        
        inquire(file = 'tsunami_simulator_with_modules_data.txt', exist = file_exists)
        if (file_exists) then
            open(9, file = 'tsunami_simulator_with_modules_data.txt', status = 'old', position = 'append', action = 'write')
        else
            open(9, file = 'tsunami_simulator_with_modules_data.txt', status = "new", action = 'write')
        end if
        
        write (9, *) n, h

    end do time_loop
     
    close(9)
    
end program tsunami

Program Explanation#

The following explanation is given by “Code Copilot” within ChatGPT 4.

The Fortran module mod_diff.f90 defines a module named mod_diff that provides numerical differentiation functionalities tailored for use in simulations or numerical methods. Here’s an overview of its content:

Module: mod_diff#

  • Dependencies: It uses the standard Fortran module iso_fortran_env for access to predefined kinds like int32 for integers and real32 for single-precision floating-point numbers.

  • Functions:

    • diff_upwind(x): This function implements the upwind differencing scheme. It is used to approximate derivatives in problems where information propagates in one direction (like fluid flow in one direction). The output dx is a vector where each element is calculated as the difference between consecutive elements of the input vector x, except the first element, which is the difference between the first and last elements, simulating a wrap-around or periodic boundary condition.

    • diff_centered(x): This function implements the centered differencing scheme, which is more accurate but also more sensitive to numerical errors. The function computes differences such that each element of the result dx is the difference between elements of x that are one step apart, with special handling for the first and last elements to manage boundary effects.

These functions are declared as pure, which in Fortran ensures that they do not have side-effects, an essential property for functions that are intended to be used in parallel computations or where determinism is critical.

Both functions are designed to operate on arrays of single-precision floating-point numbers, indicating their intended use in applications where performance and memory usage are balanced against the precision of calculations.

This module can be used in larger Fortran programs or simulations where numerical differentiation is required, such as in fluid dynamics, weather modeling, or other computational physics simulations.

Let’s proceed to analyze the other files for a complete understanding of the entire program.

The Fortran module mod_initial.f90 defines a module named mod_initial that includes functionality for initializing arrays based on specific patterns. Here’s an overview of its content:

Module: mod_initial#

  • Subroutine:

    • set_gaussian(x, icenter, decay): This subroutine initializes an array x with a Gaussian (or normal distribution) curve centered at the index icenter with a specified decay rate. This pattern is commonly used in simulations for setting initial conditions where phenomena like wave packets or other localized features are described.

Details:#

  • Input/Output: The subroutine operates directly on the array x passed to it, modifying it in place (as indicated by intent(in out)). This approach is typical in Fortran where arrays are often modified directly rather than returned.

  • Parameters:

    • icenter: An integer specifying the center position of the Gaussian distribution within the array.

    • decay: A real number that controls the width of the bell curve; higher values make the curve narrower.

  • Loop:

    • The Gaussian values are computed using a concurrent loop (do concurrent), which allows for parallel computation where supported. Each element x(i) is set to exp(-decay * (i - icenter)**2), which is the formula for a Gaussian function.

This module is particularly useful in simulations where initial conditions need to be defined in a spatially varying manner, such as setting up initial wave conditions in a medium or initial concentration distributions in simulations of diffusion processes.

We will next examine the main program file tsunami_simulator_with_modules.f90 to understand how these modules are utilized in a larger context.

The Fortran program tsunami_simulator_with_modules.f90 implements a simple model for simulating tsunami waves using numerical methods. Here’s a detailed explanation of the program:

Program Structure:#

  • Modules Used:

    • mod_diff: Provides the diff_centered function for numerical differentiation.

    • mod_initial: Contains the set_gaussian subroutine for initializing wave heights.

  • Main Variables:

    • grid_size: Defines the size of the grid (number of spatial points).

    • num_time_steps: Specifies the number of time steps for the simulation.

    • dt: Time step increment.

    • dx: Spatial step increment.

    • g: Gravitational constant (approximately 9.8 m/s²).

    • hmean: Mean height of the ocean layer.

    • h: Array representing the height of the wave at each grid point.

    • u: Array representing the velocity of the wave at each grid point.

    • icenter and decay: Parameters for the initial Gaussian distribution of the wave.

Operations:#

  1. Initial Conditions: The water height array h is initialized using a Gaussian distribution centered at icenter with a decay specified by decay. The velocity array u is initialized to zero.

  2. Output File: Data is initially written to a file tsunami_simulator_with_modules_data.txt, and the wave height at each time step is appended to this file.

  3. Simulation Loop (time_loop):

    • The velocities (u) are updated based on the gradient of the velocities and heights. The updates are scaled by dt/dx to maintain stability and accuracy.

    • The heights (h) are updated based on the divergence of the product of velocity and the sum of mean height and wave height. This mimics the conservation of mass in a fluid.

  4. File Handling: Checks if the output file exists and writes data accordingly. If the file does not exist, it is created.

Numerical Methods:#

  • The simulation uses centered differences (diff_centered from mod_diff) to compute spatial derivatives. Centered differences provide a good balance between accuracy and numerical stability.

  • The time integration method appears to be a straightforward forward Euler method, which is simple but conditionally stable.

This program is an example of how modular Fortran programming can be used to build clean and maintainable scientific simulations. The use of separate modules for different functionalities (like initial condition setting and differentiation) makes the code easier to manage and modify.

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_with_Modules"

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

import pandas as pd
output_filename = 'tsunami_simulator_with_modules_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 0.000010 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.000012 0.000029 0.000071 0.000166 0.000372 0.000801 0.001661 0.003311 0.006346 ... 1.171803e-37 8.800796e-39 6.354580e-40 4.410867e-41 2.942727e-42 1.891753e-43 1.261169e-44 0.000000e+00 3.892368e-07 9.964396e-07
2 2 0.000015 0.000037 0.000089 0.000204 0.000447 0.000944 0.001921 0.003764 0.007097 ... 9.342945e-37 7.469310e-38 5.746058e-39 4.253277e-40 3.029047e-41 2.075323e-42 1.525808e-08 3.906043e-08 1.202751e-06 3.060115e-06
3 3 0.000021 0.000050 0.000118 0.000264 0.000566 0.001169 0.002327 0.004463 0.008249 ... 6.785309e-36 5.755973e-37 4.703170e-38 3.701251e-39 5.981168e-10 1.531169e-09 7.646778e-08 1.950150e-07 2.517209e-06 6.346270e-06
4 4 0.000030 0.000070 0.000161 0.000351 0.000736 0.001487 0.002896 0.005436 0.009836 ... 4.514822e-35 4.051931e-36 2.344618e-11 6.002181e-11 4.146878e-09 1.058688e-08 2.303804e-07 5.845143e-07 4.464170e-06 1.112178e-05
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
4996 96 0.178242 0.145647 0.101621 0.100381 0.122926 0.120063 0.123888 0.116525 -0.085472 ... 3.023120e+00 5.041199e+00 4.517421e+00 2.471580e+00 9.274478e-01 2.414372e-01 -1.293256e-02 1.105624e-01 1.427006e-01 1.615230e-01
4997 97 0.166096 0.137091 -0.096517 0.101448 0.132967 0.122498 0.113014 0.109682 -0.084616 ... 1.779129e+00 4.351381e+00 5.147950e+00 3.494045e+00 1.534592e+00 4.983789e-01 9.174558e-02 -7.311907e-02 1.347228e-01 1.539230e-01
4998 98 0.150033 0.127765 -0.095495 0.104681 0.140813 0.123673 0.101828 0.102024 -0.084990 ... 7.355422e-01 3.197397e+00 5.167335e+00 4.501298e+00 2.378713e+00 9.002837e-01 2.534103e-01 -1.492602e-02 1.189377e-01 1.425529e-01
4999 99 0.130618 0.118058 -0.098317 0.109524 0.145394 0.123321 -0.091511 -0.094122 -0.086058 ... 2.072203e-02 1.919918e+00 4.503167e+00 5.186921e+00 3.417265e+00 1.503764e+00 5.057256e-01 7.775510e-02 -8.935884e-02 1.256339e-01
5000 0 0.108276 0.108241 0.104191 0.115224 0.146075 0.121337 -0.083074 -0.086573 -0.087238 ... 3.938838e-01 8.310862e-01 3.336946e+00 5.263032e+00 4.469062e+00 2.348540e+00 8.970855e-01 2.271172e-01 -3.692843e-02 1.003862e-01

5001 rows × 101 columns