Derived type implementation of the tsunami solver

Finally, we’re getting close to the home stretch. The following listing provides the (almost) complete code of the derived type implementation of the tsunami simulator. I’ve omitted the declaration section for brevity.

Listing 8.13 Derived type implementation of the tsunami solver

u = Field('u', [im, jm])                     ❶
v = Field('v', [im, jm])                     ❶
h = Field('h', [im, jm])                     ❶
hm = Field('hm', [im, jm])                   ❶
 
call h % init_gaussian(decay, ic, jc)        ❷
call h % sync_edges()                        ❷
 
hm = 10.                                     ❸
 
call h % write(0)                            ❹
 
time_loop: do n = 1, num_time_steps
 
  if (this_image() == 1) &
    print *, 'Computing time step', n, '/', num_time_steps
 
  u = u - (u * diffx(u) / dx &               ❺
         + v * diffy(u) / dy &               ❺
         + g * diffx(h) / dx) * dt           ❺
  call u % sync_edges()                      ❺
 
  v = v - (u * diffx(v) / dx &               ❻
         + v * diffy(v) / dy &               ❻
         + g * diffy(h) / dy) * dt           ❻
  call v % sync_edges()                      ❻
 
  h = h - (diffx(u * (hm + h)) / dx &        ❼
         + diffy(v * (hm + h)) / dy) * dt    ❼
  call h % sync_edges()                      ❼
 
  call h % write(n)
 
end do time_loop

❶ Initializes Fields

❷ Sets initial water height perturbation and sync

❸ Sets a constant mean water depth

❹ Writes the initial height to the file

❺ Solves for u and syncs with neighbors

❻ Solves for v and syncs with neighbors

❼ Solves for h and syncs with neighbors

How did we achieve with derived types the same form of the code for evaluating u, v, and h inside the time loop as we did with plain arrays? Before we used familiar arithmetic operators +, -, *, and / and applied them on whole arrays at once. Now, since u, v, and h are derived type instances and not arrays, something else must be going on here. The answer is in user-defined operators for derived types, which we haven’t covered yet and will explore in detail in chapter 10.

Why diffx(u) and not u % diffx()?

You may be wondering why we made all external procedures such as init_gaussian, sync_edges, and write type-bound methods of Field, except the finite difference functions diffx and diffy. The answer is simply style preference! Since diffx(u) more closely resembles the mathematical form than u % diffx(), which we certainly could’ve done, I chose to leave diffx and diffy as regular functions that take an instance of Field as an input argument.

Although I didn’t go into detail with the specific implementation of methods such as Field % sync_edges and Field % write, feel free to explore the code in src/ch08/ mod_field.f 90. You can find the main program in src/ch08/tsunami.f 90.

Running the tsunami simulator in two dimensions now produces a circular ripple. (See figure 8.2.) This will be our end result for this chapter. The effect of the initial blob in the center of the domain is the same as if we dropped a pebble in a pond. The perturbation creates a circular ripple that radiates away from the center. If you let the ripple go long enough, it will propagate through the edge and appear on the other side because of the periodic (circular) boundary conditions that are built into the Field % sync_edges method. Like in the one-dimensional case, the pond is 10 meters deep, and the perturbation is about 20 meters wide. After 3 seconds, the wave has propagated about 30 meters.

Did you know?

The theoretical phase speed of a shallow water wave–that is, the speed at which its crest moves–is equal to the square root of the gravitational acceleration (about 9.8 m/s2 in most places on Earth) times water depth. For our wave, this gives us a phase speed of 9.9 m/s, consistent with what we saw in figure 8.2. This is just one example of an emerging pattern from basic laws of physics implemented in code. Nowhere in the code did we specify how fast the wave should move, but expressing the physical laws in code on a grid-point level brought up a natural phenomenon on a larger scale. This is exactly how sophisticated dynamic models can predict a hurricane’s track days ahead without there even being a concept of a hurricane in the equation set.

That’s it for now! In the next chapter, we’ll explore generic procedures, which will allow us to use the same procedure with different input data types. We’ll also dig into redefining built-in arithmetic operators (+, -, *, /, **) for arbitrary derived types. For example, in the tsunami simulator, this will allow us to treat Field instances just like numeric arrays, and do arithmetic operations directly on them.

Run it yourself!

If you’ve cloned the application’s Git repository on GitHub, you can compile and run it like this:

make ch08