OOPS 2D Tutorial: Scalar Wave Equation

OOPS-2D Structure

The most basic OOPS-2D program contains the following components:

• A Domain object, which defines the space for your problem and contains communication information for how the physical domain is split across processors.
• A Grid object, which is a subdivision of the Domain into discrete points to be solved on a particular processor.
• A Solver object, such as RK4, which will numerically solve any PDE fed into it.
• An ODE object, which provides the right-hand side for a system of time-dependent ODEs (such as PDEs discretized with finite-difference methods), boundary conditions, and initial conditions. More complex codes will also implement Parameters and ParamParser classes, but these are not much more difficult to add in. As mentioned before, OOPS-2D works as hard as possible to hide all the gory details behind the scenes.

Setting Up a New Project

The build system for OOPS automatically searches all subdirectories for CMake build lists. Therefore, it's recommended that you set up your projects inside the OOPS directory itself. Begin by building a new directory tree for your project:

mkdir WaveEquation
cd WaveEquation
mkdir include src pars scripts

The include directory is where we will store all our header files (*.h), and src will contain all source files (*.cpp). The scripts and pars directories will not be used for the first section of this tutorial, but they will be used when we add custom parameters to the code.

Go ahead at this time and add a new file, CMakeLists.txt, into OOPS-2D/WaveEquation/. Add the following code to it:

cmake_minimum_required(VERSION 3.0)
project(WaveEquation)

set(WAVE_INCLUDE_FILES
     include/wave.h
   )
set(WAVE_SOURCE_FILES
    src/wave.cpp
    src/main.cpp
   )

set(SOURCE_FILES ${WAVE_INCLUDE_FILES} ${WAVE_SOURCE_FILES})
add_executable(Wave ${SOURCE_FILES})
target_include_directories(Wave PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/include)
target_include_directories(Wave PRIVATE ${CMAKE_SOURCE_DIR}/include)
target_link_libraries(Wave oops2d ${EXTRA_LIBS})

This file contains the information CMake needs to build our new project. The set lines tell CMake exactly which files our project needs to compile. The add_executable command specifies links the source files we've set to Wave, which will be the name of our final executable. Both target_include_directories commands specify additional locations CMake can look for header files. Finally, the target_link_libraries command links our Wave project to the core OOPS-2D library (which is compiled separately) and any additional libraries, such as MPI, that are specified in the primary build list for OOPS-2D.

The Scalar Wave Equation Code

The Scalar Wave Equation

Before we start writing actual code, let's talk about the problem we're solving. Because it's easy, we'll start with a simple wave equation on a two-dimensional sheet with outflow/advection boundaries. In two dimensions, the scalar wave equation takes the form

.

In order for OOPS-2D to solve this equation, we need to reduce the order of the time derivative. So, we introduce a new variable, , which we define as

.

Therefore, we now have a system of two first-order PDEs rather than a single second-order PDE:

,

.

We can now introduce our boundary conditions:

,

where the sign of the right-hand side depends on the boundary itself (positive for left, negative for right). The boundaries also apply to the y-direction, simply substituting y for x, down for left, and up for right. This guarantees that any data that hits the edge of the grid will simply propagate off.

The Wave Class

Now we're in a position to start writing code. In the include directory, create a new file, wave.h, and add the following code.

#ifndef WAVE_H
#define WAVE_H

#include <ode.h>

class Wave : public ODE{
  private:
    // Variable labels
    static const unsigned int U_PHI = 0;
    static const unsigned int U_PI  = 1;

  protected:
    virtual void rhs(std::shared_ptr<FieldMap>& fieldMap);

  public:
    Wave(Domain* d, Solver* s);
    virtual ~Wave();

    virtual void initData();
};
#endif

The ode.h file contains the declaration of the ODE class, which we must inherit from. ODE is an abstract class which provides a useful framework for things like our right-hand side, boundary conditions, and initial data. When OOPS-2D is run with multiple cores, ODE also takes care of transferring data between neighboring grids, and it also takes care of formatting data for output.

Underneath the private label, we define two constant variables, U_PHI and U_PI, which are simply indices. These are not strictly necessary, but they are good practice because they make for more readable code. For larger systems of variables, some users may prefer to declare an enumerator instead.

The protected label contains (for now) one function, rhs. This is a virtual function declared in ODE, which is used by Solver to evolve the system of equations forward in time. The std::shared_ptr<FieldMap> type is a smart pointer (specifically, an std::shared_ptr) to FieldMap, which contains the data for all the various fields our ODE will evolve. (Users unfamiliar with the concept of a "smart pointer" may want to read more about them. They are a useful feature introduced in C++11 which, used properly, offer no noticeable performance decrease over traditional pointers while greatly reducing the risk of introducing memory leaks into the code.)

Inside the public function, we define a constructor (which requires Domain and Solver objects), a destructor (which, as Wave is a derived class, is declared as virtual), and initData, which is an inherited function we can override to set up our problem.

Following this, we need to set up main.cpp inside WaveEquation/src/, where we will define the Wave class's methods. We'll start with the constructor and destructor:

#include <wave.h>
#include <iostream>
#include <cmath>

// Constructor
Wave::Wave(Domain* d, Solver* s) : ODE(2, d, s){
  addField("Evolution", nEqs, true, true);

  reallocateData();
}

// Destructor
Wave::~Wave(){

}

We make sure to include the wave.h class we just created, as well as some utility classes from the Standard Template Library. The destructor is completely empty, as any memory we allocate is handled by ODE or through a smart pointer.

The constructor itself is a little more interesting. Right off the bat, we have to make a call to the ODE constructor, which has the following prototype:

ODE::ODE(unsigned int nEqs, Domain* d, Solver* s)

The first argument is just the number of equations in our ODE system. With the more modern FieldMap interface functionality implemented, this is somewhat obsolete, but it can be a useful piece of information to hold onto for information about the primary evolved field. The second and third arguments provide pointers to the Domain and Solver objects, respectively, that will be used by the ODE object.

The next line makes a call to addField:

Result addField(std::string name, unsigned int nEqs, bool isEvolved, bool isComm)

The addField function creates a new set of variables for the ODE object. The first parameter, name, assigns a unique name to these variables so they can be accessed via FieldMap. The second parameter, nEqs, defines how many equations this particular set of fields has. The third parameter, isEvolved, is a boolean value that tells ODE whether or not this particular field needs to be evolved or not. This can be useful for defining utility fields (such as constraint equations) which shouldn't be evolved by the Solver object and will be handled manually by the user. The final parameter, isComm, tells ODE whether or not this field needs to be communicated between different processors. Again, this adds flexibility for fields whose values may not be independent of a field which is communicated.

In our particular case, we call our field "Evolution", state that it has two variables, and that it will be evolved by the Solver object attached to ODE and communicated across processors.

Next, we need to write rhs. This particular function has a lot of pieces, so we'll take it one step at a time.

void Wave::rhs(std::shared_ptr<FieldMap>& fieldMap){
  const Grid& grid = fieldMap->getGrid();
  unsigned int nx = grid.getSize()[0];
  unsigned int ny = grid.getSize()[1];
  double dx = grid.getSpacing()[0];
  double dy = grid.getSpacing()[1];
  unsigned int nb = domain->getGhostPoints();

  double **dudt = fieldMap->getSolverField("Evolution")->getCurrentRHS();
  double **u = fieldMap->getSolverField("Evolution")->getIntermediateData();

  unsigned int xstart = (domain->hasBoundary(LEFT)) ? nb + 1 : nb;
  unsigned int xend = (domain->hasBoundary(RIGHT)) ? nx - nb - 1 : nx - nb;
  unsigned int ystart = (domain->hasBoundary(DOWN)) ? nb + 1 : nb;
  unsigned int yend = (domain->hasBoundary(UP)) ? ny - nb - 1 : ny - nb;

The first few lines here are just grabbing information we need from the local Grid and Domain objects. The nx and ny variables are the dimensions of the Grid, and dx and dy are the spacing between neighboring points in x and y, respectively. The nb variable is the number of ghost points. Ghost points in OOPS-2D serve two purposes: 1) they allow for data to be communicated between grids, and 2) they allow for fluid-style boundary conditions. The next couple of lines reference the "Evolution" field we created in our constructor. We grab the current right-hand side for our variables (or the right-hand side associated with the current Solver stage) and the "intermediate" data, the temporary data calculated based on the current stage.

The last four lines help us determine where the interior region of the Grid is. Because OOPS-2D, unlike OOPS, can only have one Grid per Domain per processor, boundary information is stored in Domain, not Grid. If we're on a physical boundary, we don't need to include the first physical point in our calculations; it will just be overwritten by the boundary conditions. However, if the edge of the grid is a processor boundary, that first physical point is actually an interior point.

The next lines are

  for(unsigned int j = 0; j < ny; j++){
    for(unsigned int i = 0; i < nx; i++){
      unsigned int xy = grid.getIndex(i, j);
      dudt[U_PHI][xy] = 0.0;
      dudt[U_PI][xy] = 0.0;
    }
  }

This section of code simply clears out the right-hand side routine. While not strictly necessary, particularly for the interior points, it does simplify the boundary calculations because we can write them in a way that doesn't require a special treatment for the corners. The getIndex(unsigned int x, unsigned int y) function is a good way to introduce how OOPS-2D stores grid data. Because of the idiosyncrasies of transferring data via MPI and handling multidimensional arrays in C++, as well as a desire to make the code somewhat more general (should an OOPS-3D or an OOPS-ND get built), all data for a single variable is stored in a single flat array. Every Grid object contains an inlined indexing function to make data access easier. Because all the x-points are stored contiguously, the loop order should always put x as the innermost loop where possible.

  for(unsigned int j = ystart; j < yend; j++){
    for(unsigned int i = xstart; i < xend; i++){
      unsigned int xy = grid.getIndex(i, j);
      dudt[U_PHI][xy] = u[U_PI][xy];
      dudt[U_PI ][xy] = (u[U_PHI][grid.getIndex(i+1,j)] - 2.0*u[U_PHI][xy] +
                         u[U_PHI][grid.getIndex(i-1,j)])/(dx*dx) +
                        (u[U_PHI][grid.getIndex(i,j+1)] - 2.0*u[U_PHI][xy] +
                         u[U_PHI][grid.getIndex(i,j-1)])/(dy*dy);
    }
  }

This next set of loops is where we calculate the right-hand side for the interior points. More experienced readers should recognize the U_PI calculation as a second-order finite-difference Laplacian operator, or the sum of second-order finite-difference second derivative operators for the x and y directions.

The final piece is to apply boundary conditions. This is done as follows:

  // Left side.
  if(domain->hasBoundary(LEFT)){
    for(unsigned int j = 0; j < ny; j++){
      unsigned int xy = grid.getIndex(nb,j);
      unsigned int x1y = grid.getIndex(nb + 1, j);
      unsigned int x2y = grid.getIndex(nb + 2, j);
      dudt[U_PHI][xy] += (-3.0*u[U_PHI][xy] + 4.0*u[U_PHI][x1y] - u[U_PHI][x2y])/(2.0*dx);
      dudt[U_PI ][xy] += (-3.0*u[U_PI ][xy] + 4.0*u[U_PI ][x1y] - u[U_PI ][x2y])/(2.0*dx);
    }
  }
  // Right side.
  if(domain->hasBoundary(RIGHT)){
    for(unsigned int j = 0; j < ny; j++){
      unsigned int xy = grid.getIndex(nx - nb - 1,j);
      unsigned int x1y = grid.getIndex(nx - nb - 2, j);
      unsigned int x2y = grid.getIndex(nx - nb - 3, j);
      dudt[U_PHI][xy] += (-3.0*u[U_PHI][xy] + 4.0*u[U_PHI][x1y] - u[U_PHI][x2y])/(2.0*dx);
      dudt[U_PI ][xy] += (-3.0*u[U_PI ][xy] + 4.0*u[U_PI ][x1y] - u[U_PI ][x2y])/(2.0*dx);
    }
  }
  // Bottom side.
  if(domain->hasBoundary(DOWN)){
    for(unsigned int i = 0; i < nx; i++){
      unsigned int xy = grid.getIndex(i,nb);
      unsigned int x1y = grid.getIndex(i, nb + 1);
      unsigned int x2y = grid.getIndex(i, nb + 2);
      dudt[U_PHI][xy] += (-3.0*u[U_PHI][xy] + 4.0*u[U_PHI][x1y] - u[U_PHI][x2y])/(2.0*dy);
      dudt[U_PI ][xy] += (-3.0*u[U_PI ][xy] + 4.0*u[U_PI ][x1y] - u[U_PI ][x2y])/(2.0*dy);
    }
  }
  // Top side.
  if(domain->hasBoundary(UP)){
    for(unsigned int i = 0; i < nx; i++){
      unsigned int xy = grid.getIndex(i,ny - nb - 1);
      unsigned int x1y = grid.getIndex(i, ny - nb - 2);
      unsigned int x2y = grid.getIndex(i, ny - nb - 3);
      dudt[U_PHI][xy] += (-3.0*u[U_PHI][xy] + 4.0*u[U_PHI][x1y] - u[U_PHI][x2y])/(2.0*dy);
      dudt[U_PI ][xy] += (-3.0*u[U_PI ][xy] + 4.0*u[U_PI ][x1y] - u[U_PI ][x2y])/(2.0*dy);
    }
  }
}

Because this Domain/Grid pair may not be on a processor with a physical boundary, we have to manually check for physical boundaries. If we are on a boundary, we calculate indices for the boundary point itself (xy), the first interior point (x1y), and the second interior point (x2y). We then apply a second-order finite-difference first derivative operator:

.

Because the sign of the operator flips for the backward/left-going derivative, this cancels out the sign change required by the boundary condition on the right boundary, and the calculation is identical for every boundary. Also, note that the boundaries are added, rather than just copied. This ensures that our corner points, which may have both x and y data propagating off, are handled correctly.

After that mouthful, we have one more (easy) function to worry about: the initData function, which applies our initial conditions.

void Wave::initData(){
  pair2<double> bounds = domain->getBounds();
  double x0 = 0.5*(bounds[0][1] + bounds[0][0]);
  double y0 = 0.5*(bounds[1][1] + bounds[1][0]);
  double amp = 1.0;
  double sigma = 0.125;
  // Just assume a Gaussian for the time being.
  unsigned int nx = domain->getGrid()->getSize()[0];
  unsigned int ny = domain->getGrid()->getSize()[1];

  auto evol = (*fieldData)["Evolution"];
  auto points = domain->getGrid()->getPoints();
  const double *x = points[0];
  const double *y = points[1];
  double **u = evol->getData();
  for(unsigned int j = 0; j < ny; j++){
    for(unsigned int i = 0; i < nx; i++){
      unsigned int xy = domain->getGrid()->getIndex(i,j);
      u[U_PHI][xy] = amp*std::exp(-((x[i] - x0)*(x[i] - x0) + (y[j] - y0)*(y[j] - y0))/(sigma*sigma));
      u[U_PI ][xy] = 0.0;
    }
  }
}

Much like the right-hand side function, the first few lines are largely just collecting some data we'll be using throughout. We set up some parameters for a Gaussian (we'll generalize this with an actual Parameters object later), which we want placed at the center of the physical domain. Notice that the "Evolution" field is accessed slightly differently, this time by the [] operator instead of getSolverField. The FieldMap class actually stores two std::map objects: one for smart pointers to SolverData objects, like "Evolution", and one for smart pointers to ODEData objects. Because SolverData is a derived class from ODEData, all fields, including "Evolution", are stored in the ODEData map and can be accessed via the [] operator like a traditional array, but they will be returned as ODEData objects. Accessing their SolverData methods (such as getCurrentRHS and getIntermediateData) requires a potentially dangerous and aesthetically displeasing dynamic_cast or returning the pointer via getSolverField instead. However, we only need access to the getData method here, so the [] operator will work fine. The rest of the code is more straightforward: we loop over all the data, setting U_PHI to a Gaussian distribution and U_PI to zero.

The main Function

The last piece of code is very straightforward. Create a file in OOPS-2D/WaveEquation/src titled main.cpp:

#include <mpicommunicator.h>
#include <wave.h>
#include <rk4.h>
#include <iostream>
#include <cstdio>
#include <cmath>

int main(int argc, char *argv[]){
  MPICommunicator *comm = MPICommunicator::getInstance();
  Result result = comm->init();
  if(result != SUCCESS){
    std::cout << "There was an error initializing MPICommunicator.\n";
    return 0;
  }

  // Set up the PDE
  Domain domain = Domain();
  pair2<double> bounds;
  bounds[0][0] = -1.0;
  bounds[0][1] = 1.0;
  bounds[1][0] = -1.0;
  bounds[1][1] = 1.0;
  domain.setBounds(bounds);
  domain.setGhostPoints(1);

  unsigned int shp[2] = {101, 101};

  result = domain.buildMesh(shp);

  if(result != SUCCESS){
    std::cout << "There was an error building the mesh.\n";
    std::cout << "  Error Code: " << result << "\n";
  }

  // Set up the ODE system.
  RK4 rk4 = RK4();
  Wave ode = Wave(&domain, &rk4);
  ode.initData();
  
  double ti = 0.0;
  double tf = 2.0;
  auto dx = domain.getGrid()->getSpacing();
  double dt = domain.getCFL()*fmin(dx[0], dx[1]);
  unsigned int M = (tf - ti)/dt;

  ode.dumpField("Evolution", "phi00000.csv", 0, 0);
  for(unsigned int i = 0; i < M; i++){
    double t = (i + 1)*dt;
    ode.evolveStep(dt);

    char buffer[12];
    sprintf(buffer, "phi%05d.csv", i+1);
    ode.dumpField("Evolution", buffer, t, 0);
  }
  
  result = comm->cleanup();
  if(result != SUCCESS){
    std::cout << "There was an error cleaning up MPICommunicator.\n";
  }
  return 0;
}

Inside our main function, the first thing we have to do is initialize MPI. OOPS-2D makes this very easy with the use of MPICommunicator. MPICommunicator uses a singleton design pattern to ensure that there is only ever one instance of the object, which must be accessed via the getInstance static method. We call the init function from the communicator and quit if there were any issues.

After that, we set up the Domain and its associated Grid object. This includes setting the bounds to the region [-1,1]x[-1,1], and the number of ghost points. We set this to 1, which guarantees that there is enough data at the interprocessor boundary to get a full three-point stencil for our Laplacian operator.

We next set how many points our Grid should have and call buildMesh(unsigned int shp[2]). Adding a Grid in OOPS is a fairly simple operation; the Domain just checks for overlaps and boundary issues, then stores it in a sorted set. In OOPS-2D, however, this is a more delicate operation. Each processor calculates its position on the Domain, then uses the supplied parameters to determine the boundaries and number of points for its Grid object.

After this is done, we can (finally) set up our Wave ODE. We create an RK4 solver object, pass it and our Domain into the constructor for Wave, and apply the initial data.

The next piece of the code is the evolution loop. We figure out our time step, the total number of steps, and then run evolveStep inside a loop over all time. The dumpField method allows us to select a particular field ("Evolution", in this case) and save a particular variable from that field. The output data is in a .csv format, with one file per frame and the columns being (time, x-coordinate, y-coordinate, value). The prototype for this function is

dumpField(std::string field, char *name, double time, unsigned int var),

where field is the name of the field stored in the ODE object's FieldMap object, name is the filename, time is the current time, and var is the index of the variable inside field that should be output.

Finally, we have to shut down MPI, which is done with a simple call to the MPICommunicator object's cleanup function.

Wrapping Up

Return to the OOPS-2D/build and compile again. The build system should automatically find WaveEquation and try to compile the code. If you're lucky and put everything in correctly, you should have a directory titled WaveEquation containing an executable named Wave. Run it (mpirun -np <number of cores> Wave) and make sure that it generates a whole mess of data files. These .csv data files can be plotted using a tool like Paraview or Python's Matplotlib module. Make sure that everything looks like what you expect.

At this point, you're probably saying, "Wait a minute! That took forever! You said this would be easy!"

If all you wanted was a simple wave equation, you're probably right. This was a lot of work. You could probably have solved the same problem in MATLAB or Python in a matter of minutes. But let's talk about what you didn't have to do. You didn't have to write your own ODE integrator; you just used the built-in RK4 routine. You didn't have to worry about discretizing the Domain yourself; while not particularly difficult, it can be a common source of errors. You didn't have to worry about allocating or cleaning up any memory yourself. OOPS-2D did it all for you. Lastly, and perhaps most importantly, you didn't really have to do anything with MPI at all. You made a quick initialization call at the beginning and a cleanup call at the end, and that was it. OOPS did all the heavy lifting and split up the Grid objects appropriately, handled all the data transfer, and even collated the individual .csv files produced by every processor during dumpField into a single file.

Now, let's talk about what you can do with OOPS-2D and where it really shines. Later on, you'll learn how to add parameters to your code, which is a royal pain to do yourself but turned into a nearly trivial task by OOPS-2D. If you had to deal with a more complicated system of equations that involved separate unevolved constraint variables or utility quantities, you would just call addField in your ODE derivative's constructor for the extra fields. You could even

Extensions and Exercises

To make sure that you understand how OOPS-2D works, here are some exercises worth trying:

1. Inside the rhs definition in wave.cpp, we looped over the entire right-hand side to clear it. Try to modify this loop so it only applies to the boundaries.
2. The ODE class has an empty function, virtual void applyBoundaries(), which automatically gets called inside evolveStep. Try overriding this function in Wave and using it to apply fluid-style boundaries (i.e., boundaries using ghost points) instead of writing an advection condition in rhs. (A hint for OOPS users: there is no longer an intermediate flag passed into this function; all boundaries should be applied to the intermediate data.)
3. Try rewriting the code to use a higher-order derivative operator. This website may be useful.