Overriding operators in the tsunami simulator

You now have a good idea of how defining an arithmetic operator for use with derived types works. Let’s apply this knowledge to implement the full arithmetic for the Field type in the tsunami simulator. We implemented this derived type back in chapter 8, where we used it to model the physical properties of the water flow, such as water surface height and velocity. Using a derived type instead of plain Fortran arrays allowed us to abstract away much of the boilerplate code into the constructor function and type-bound methods. Although I didn’t explain how we made the arithmetic operators to work with derived types, we used this concept to code our equations of motion and continuity, just like we did with regular arrays, as shown in the following listing.

Listing 10.7 Excerpt from the Tsunami simulator main time loop

time_loop: do n = 1, 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()                       ❹
  ...
end do time_loop

❶ Iterates for num_time_steps time steps

❷ Calculates the new value of velocity in the x axis and synchronizes

❸ Calculates the new value of velocity in the y axis and synchronizes

❹ Calculates the new value of water height and synchronizes

This is the core of the tsunami simulator. First, two equations evaluate the water velocity components u and v in the x and y axis, respectively, based on the slope of the water surface. The higher diffx(h) / dx and diffy(h) / dy are, the steeper the water surface in the x and y directions, respectively, becomes. A steeper water surface leads to an increase in water velocity. The third equation evaluates the water surface height, depending on how much water is coming into or out of the grid cell.

After updating each Field variable, we make a call to sync_edges. This method updates (or, as we also call it, synchronizes) the edge values with those from neighboring parallel tiles. sync_edges is what allows the features of the flow to propagate through from one parallel CPU to another. There are three important points to raise here:

• In serial mode (single core), sync_edges does nothing.
• sync_edges is required to make the parallel program work in a mechanical sense but is not otherwise meaningful for the physics of the simulator.
• Parallel synchronization of a field is required only when we update the values of that field, not otherwise.

Considering these points, can we find a way to abstract away the explicit synchronization via sync_edges in parallel mode? If we’re designing tsunami as a library to allow users to design custom fluid dynamics solvers, we shouldn’t expect them to also write low-level synchronization calls. Natively parallel code should look the same as its serial counterpart.

Having worked through this chapter up to this point, you have all you need to know to implement the arithmetic operators for the Field derived type. We’ll go over the details of this implementation in this section. This won’t involve any new Fortran features–we’ll just apply the existing knowledge from implementing the countdown app to the tsunami simulator. However, we’ll also use this opportunity to take it one step further and build the synchronization method (sync_edges) into the assignment operator (=) itself! While this involves no new syntax elements, it’s a powerful little hack that will allow you to build intricate functionality into mundane statements such as assignments or arithmetic.