A Basic Example
For this project, we'll create a simple program using OOPS that solves the wave equation on a domain with outflow boundaries. Being by building a new project directory for the program:
cd <path to OOPS base directory>
mkdir WaveEquation
cd WaveEquation
mkdir include
mkdir src
The include folder will contain any header files specific to the WaveEquation project, and src will hold the main file and any additional source files required. Go ahead and navigate to OOPS/WaveEquation and add a new file, CMakeLists.txt:
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 oops ${EXTRA_LIBS})CMake is a build system that generates custom Makefiles, which are used to compile your code. Although its syntax is somewhat baroque, it is very easy to write one configuration file that works on nearly all systems rather than manually writing Makefiles for each machine. Very simply, we define a project with the project property, then set some convenience variables WAVE_INCLUDE_FILES and WAVE_SOURCE_FILES that state what the header and source files are for this project. We then state that we want to add a new target in the form of an executable, Wave, which should be compiled from the header and source files we defined earlier. The target_include_directories command simply says where Wave should look for the source files we've included, and then target_link_libraries will link Wave to the OOPS library and any additional libraries (like SDF) selected during configuration.
Before we start writing our code, remember that the wave equation in one dimension takes the form
OOPS can't solve this in its current form because the time integrator expects it to be first order in space. Therefore, we make the definition , which suggests the ODE system
,
.
OOPS can actually solve this as is because we can just define a second-order derivative operator, but it's also easy to reduce it one more time to a set of three first-order equations with the definition , yielding the three equations
,
,
,
where the third equation is derived by recognizing that partial derivatives must commute. Despite the two systems being mathematically identical and producing valid solutions to the wave equation, they exhibit very different numerical properties, with the first-order system being more dispersive and the second-order system being more diffusive. (Additionally, notice that in the first-order formulation, neither nor
depend on
any longer, reducing the solution to the latter to a simple integration problem.)
To solve these equations, we also need to define adequate boundary conditions. For an outflow condition, we want the wave to move off the grid smoothly without affecting the interior solution at all. This corresponds to the condition
at the left (+) and right (-) boundaries, with similar conditions for and
. The astute reader may say, "Wait a second! We just had a second-order system that only needed two boundaries, and now we have six! Where are the extra boundaries coming from? Isn't the system overdetermined?" Technically speaking, this is true. We only have two unique boundary conditions; the remaining four, to be valid boundaries, must satisfy a set of constraint equations. In the case of the scalar wave equation, these constraints are just the additional definitions we used to reduce the order of the original PDE. Remember that
is the solution to the wave equation -- we can impose boundaries on the wave equation, then use
and
to derive their boundaries (which happen to be exactly the same in our case). For other systems, such as BSSN formulation of the Einstein equations, these constraints are considerably more complicated.
For the system to be well-defined, the one last piece comes in the form of initial conditions. Like when we defined our boundaries, our system is overdetermined: we only have two initial conditions in the original wave equation, but now we have three first-order equations in time. Therefore, two of our boundaries cannot be independent. The usual initial conditions for the wave equation are for and
, so we are free to choose initial conditions for
and
. However, the condition for
must satisfy the constraint
, so if we define
,
,
We must also define
.
We can now talk about writing down our PDE in OOPS. In the Wave/include directory, create a new file, wave.h, that reads as follows:
#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;
static const unsigned int U_CHI = 2;
protected:
virtual void rhs(std::shared_ptr<FieldMap>& fieldMap);
public:
Wave(Domain& d, Solver& s);
virtual ~Wave();
virtual void initData();
};
#endifThe ode.h file defines the ODE class, which is an abstract class representing a general system of ODEs (or, rather, a discretized system of PDEs). Our own PDE, the wave function, will inherit from this class. So, we define a class, Wave, which inherits directly from ODE, and declare its contents. The variable labels under the private section are for convenience -- rather than accessing our variables as a numerical index, we will access them using their names.
Under the protected section, we declare the function rhs(). Because this is an abstract function coming from ODE, we declare it as virtual. This is the meat of our ODE object, and it takes a smart pointer, fieldMap, which will contain the data for all our fields (U_PHI, U_PI, and U_CHI). (If you're not familiar with the concept of a "smart pointer", it's worth reading up on; they're basically C++ pointers that handle allocation and deallocation themselves to prevent memory leaks and access errors.) The rhs() function defines the right-hand side of our equations, so it's where we will discretize our spatial derivatives and apply our boundary conditions.
Finally, underneath the public section, we set up a constructor, a destructor (which should be virtual because ODE is abstract), and define initData, which is where we apply the initial data. Note that the constructor takes a Domain object, which defines the physical domain for the problem, and a Solver object, which specifies the numerical integrator that we should use, like a 4th-order Runge-Kutta integrator.
We now need to implement these methods. Inside Wave/src, create another new file, wave.cpp. The implementation is quite a bit longer, so we'll take it one piece at a time. First, we need the following includes:
#include <wave.h>
#include <operators.h>
#include <solverdata.h>
#include <iostream>
#include <cmath>
#include <memory>The wave.h file is what we just created above, operators.h contains some finite-difference operators, and solverdata.h will allow us to access both our solution data and the righthand side.
Next, we'll write the constructor and destructor:
Wave::Wave(Domain &d, Solver &s) : ODE(3, 1){
domain = &d;
solver = &s;
addField("Evolution",nEqs, true, true);
reallocateData();
}
Wave::~Wave(){}The destructor is straightforward; because we handle everything through smart pointers, there's nothing to free at the end, and we don't need to do anything.
The constructor merits a little more discussion. We first call ODE's constructor explicitly. The ODE class takes two const unsigned int arguments, the first of which defines the number of variables (n = 3 here), and the second an id (id = 1 here). This second argument is a relic from an older version of OOPS and a much more perverse parameter system. Consequently, it doesn't really matter what its value is, and it's likely to disappear in a future version of OOPS. Next, we take the Domain and Solver objects passed into the constructor and save them to Wave so that we can reference them later.
The last two functions deserve a little bit more discussion. All data in OOPS is stored in a FieldMap object. Therefore, all the equations we're solving need to be added to the FieldMap. The function addField() takes four arguments: an std::string representing a name (which must be unique), an unsigned int stating how many equations are in our system, and two bool parameters, isEvolved and isComm. The first parameter, isEvolved, indicates whether or not the field should be evolved, i.e., whether or not it's going to be integrated directly. The second parameter, isComm, has to do with communication across multiple Grid objects, which will be discussed in a different tutorial. For now, just set it to true and be done with it. Therefore, the function call
addField("Evolution", nEqs, true, true)indicates that we want to add a new field labeled "Evolution" that has 3 variables (nEqs is set by the ODE constructor), and we intend to evolve these variables forward in time and communicate them across Grid boundaries.
Now we can write the rhs() function. This is a bit more complicated, so we'll take it one piece at a time.
void Wave::rhs(std::shared_ptr<FieldMap>& fieldMap){
unsigned int nb = domain->getGhostPoints();
// Check that the grid is actually big enough
if(grid.getSize() < 3){
printf("Grid is too small. Need at least 3 points.\n");
return;
}
double **dudt = fieldMap->getSolverField("Evolution")->getCurrentRHS();
double **u = fieldMap->getSolverField("Evolution")->getIntermediateData();
// Define some variables we'll need.
double stencil[3] = {0.0, 0.0, 0.0};
double dx = grid.getSpacing():
int shp = grid.getSize();
}The first section has to deal with the computational domain. variable nb represents the number of "ghost points" in the system. Ghost points are additional points added to either end of a Grid, and they're used for a variety of things, including communication and special boundary conditions. The grid object grabs a reference to the Grid that the FieldMap is defined on. Note that while a Domain defines the physical domain of the PDE, it does nothing else; the Grid object defines the actual coordinate points that we can use to perform calculations. If this Grid doesn't have enough points, we can't do any calculations, so we have to quit.
Next, we grab the data from the FieldMap so that we can perform our calculations. Notice that both the solution (u) and the righthand side (dudt) are stored in 2d arrays: the first index corresponds to the variable (U_PHI, U_PI, or U_CHI), and the second to the index along the Grid. Users familiar with older versions of OOPS may notice that we are grabbing "intermediate data" for the solution rather than checking the stage of the Solver. OOPS 1.2 takes care of that for you by copying the first stage into the intermediate data, so users should always grab the intermediate data to get the current solution.
The last piece defines some convenience variables. The stencil array is used to help calculate our spatial derivatives, dx defines the grid spacing between points, and shp says how many points are stored in the Grid total.
Next, we can worry about actually calculating the righthand side for our PDE:
// Interior points
for(int i = nb + 1; i < shp - nb - 1; i++){
dudt[U_PHI][i] = u[U_PI][i];
for(int j = 0; j < 3; j++){
stencil3[j] = u[U_CHI][i - 1 + j];
}
dudt[U_PI][i] = operators::dx_2(stencil3, dx);
for(int j = 0; j < 3; j++){
stencil3[j] = u[U_PI][i - 1 + j];
}
dudt[U_CHI][i] = operators::dx_2(stencil3, dx);
}We define a loop over all the interior points, which excludes the ghost points plus the boundary points. The proceeding mess is just our equations that we wrote above earlier. The operators::dx_2() function indicates that we want to calculate a second-order-accurate finite-difference approximation to the first derivative with a spacing dx, which takes the mathematical form
,
where the subscript i indicates the spatial index on the grid. Therefore, a second-order finite-difference to the _i_th point on the grid requires the points on either side. We store all three of these points in stencil, which then gets passed into operators::dx_2(). For a simple derivative like this, it is just as simple to write it out explicitly, but this becomes much more difficult for higher-order derivative operators, so it is much more convenient to use the functions defined in operators.
The last part of rhs involves boundary conditions:
// Boundary points
for(unsigned int m = 0; m < nEqs; m++){
// Left boundary
for(unsigned int i = 0; i < 3; i++){
stencil3[i] = u[m][nb + i];
}
dudt[m][nb] = operators::dx_2off(stencil3, dx);
// Right boundary
for(unsigned int i = 0; i < 3; i++){
stencil3[i] = u[m][shp - nb - 1 - i];
}
dudt[m][shp - nb - 1] = operators::dx_2off(stencil3, dx);
}
}For our particular situation, the ghost points that may or may not exist don't carry physical information, so we can't calculate a derivative using points on either side of the boundary. Therefore, we use an asymmetric finite-difference operator, operators::dx_2off, which takes the _i_th point and the two points to its right to calculate the first derivative. This gives us the left boundary. The operator is antisymmetric, so if we decide to use the two points to the left instead, we get a left-going first derivative, but with an extra negative. This happens to be exactly what we need for the right boundary, so we push the last three physical points in reverse order ([i], [i-1], and [i-2]) into our stencil and calculate the derivative.
Now that we have the righthand side written, the last step in our Wave class is the initial conditions. We will use a steep Gaussian defined as follows:
,
which has the spatial derivative
,
and set our initial wave velocity to zero. This is written in our code as
void Wave::initData(){
double x0 = 0.5*(domain->getBounds()[1] + domain->getBounds()[0]);
double A = 1.0;
double wsq = 64.0;
for (auto it = fieldData.begin(); it != fieldData.end(); ++it){
auto evol = (**it)["Evolution"];
const double *x = evol->getGrid().getPoints();
unsigned int nx = evol->getGrid().getSize();
double **u = evol->getData();
for (unsigned int i = 0; i < nx; i++){
double f = A*std::exp(-wsq*(x[i] - x0)*(x[i] - x0));
u[U_PHI][i] = f;
u[U_PI ][i] = 0.0;
u[U_CHI][i] = -2.0*wsq*(x[i] - x0)*f;
}
}
}The first three lines just define our constants; these are things that should be generalized via a Parameters object, but for now we'll hardcode them to fixed values. We set the center of the Gaussian, x0, to be in exactly the middle of the computational domain, then define our amplitude A as 1 and the width w of the Gaussian to 8.
The exterior loop requires a little bit of explanation. Note that OOPS supports multiple Grid objects on a single Domain. There is a unique data set for each Grid, so we have to loop over all of them. The line that follows, which has the rather baroque construction (**it)["Evolution"], is just part of grabbing the right data. We can then grab our Grid coordinates, x, the size of the Grid, nx, and the actual points, u. The loop that follows is where we apply the initial data, as discussed above.
Before we can compile our code, we need to define a main function that sets up some options and runs the main loop. In Wave/src, create a new file, main.cpp, which we define as follows:
#include <wave.h>
#include <domain.h>
#include <grid.h>
#include <rk4.h>
#include <polynomialinterpolator.h>
int main(int argc, char* argv[]){
Domain domain = Domain();
domain.setGhostPoints(0);
domain.setCFL(0.25);
int N0 = 101;
double bounds[2] = {0.0, 1.0};
domain.addGrid(bounds, N0);
RK4 rk4 = RK4();
PolynomialInterpolator interpolator = PolynomialInterpolator(4);
Wave ode = Wave(domain, rk4);
ode.setInterpolator(&interpolator);
ode.initData();
double ti = 0.0;
double tf = 2.5;
double dt = domain.getCFL()*(domain.getGrids().begin())->getSpacing();
unsigned int M = (tf - ti)/dt;
ode.dumpCSV(std::string("Evolution"), std::string("Phi00000.csv"), 0, 0);
for(unsigned int i = 0; i < M; i++){
double t = (i + 1)*dt;
ode.evolveStep(dt);
char name[32];
sprintf(name, "Phi%05d.csv", i+1);
ode.dumpCSV(std::string("Evolution"), std::string(name), t, 0);
}
return 0;
}The includes are all for OOPS-related things: wave.h, of course, is our own Wave function, and domain.h, grid.h, rk4.h, and polynomialinterpolator.h give us access to the Domain, Grid, RK4, and PolynomialInterpolator objects.
Inside our main() function, we first set up our Domain object. Because we're only using a single Grid and have no need for special boundary conditions, we set the number of ghost points to 0. We also set the Courant-Friedrichs-Lewy (CFL) condition, which defines the desired ratio between our time step and our spatial interval
. This is incredibly important for stability with our schemes, and unfamiliar readers should consult a standard text on numerical methods for differential equations. Finally, we define that we want a
Grid of 101 points (or 100 cells), and we set our domain to cover the interval [0, 1]. Finally, we construct a Grid object with these bounds and size.
The next section involves setting up the ODE object itself. We first need a numerical integrator, which is part of the Solver class. For this example, we'll use the RK4 integrator, which is an implementation of the classic fourth-order Runge-Kutta method. We next define an Interpolator object, in this case PolynomialInterpolator, which tells the ODE how to interpolate between Grid objects. It serves no purpose for our example, but it's a necessary part of the ODE object setup. Finally, we define an ODE of type Wave using our Domain and Solver objects, attach the Interpolator, then apply our initial conditions.
The last segment has to do with evolution and data output. The dumpCSV() function takes four parameters:
dumpCSV(std::string field, std::string filename, double time, unsigned int var)The field argument states which field we're interested pulling data from, and should be called by the unique field name that we setup in the Wave class's constructor. The filename argument says what the output file should be named, and the time argument says what the output time should be. The last argument, var, is the variable index. In our case, 0 corresponds to U_PHI, so we'll output U_PHI.
The loop simply increments the frame index in our filename, calls ode.evolveStep(dt) to evaluate the next step in our evolution, and outputs the new frame of data. It will loop until the time exceeds tf.
For a Unix or Unix-like system, compilation is very easy thanks to CMake. Navigate to the OOPS base directory and do the following:
mkdir build
ccmake ../(This may require you to install the curses interface for CMake.) Simply follow the instructions in the terminal to configure and generate the Makefiles, then run make to compile. A directory named Wave should appear, which will contain your Wave binary.
Make sure that your version of Visual Studio is configured to use CMake. This is the default option as long as you enabled the "Desktop development with C/C++" or "Linux Development with C++" options during installation. After that, compiling OOPS is relatively simple.
- Visual Studio should automatically detect the
CMakeLists.txtfile and generate a debug CMake build configuration labeledx86-Debug. The default options should be sufficient for testing purposes, but you can use the CMake Settings wizard to modify and add new configurations in CMakeSettings.json. - Compile the code by clicking on "Build > Build All" or hitting F7.
For either system, the CMakeLists.txt file in the OOPS base directory should automatically detect the Wave project and add it to the build procedure. If all goes well, there should be no build errors, and running the code is as simple as running the Wave executable. A very large number of Phi*.csv files should start appearing in the directory, which is the output specified during the main loop.
The data can be viewing with the plotData.py Python script located in the OOPS/scripts directory. Copy it to your output directory, then run the command python3 plotData.py *.csv. It will generate a new folder, plots_data, which will generate a simple plot for each data file, and should look like a Gaussian that splits into two waves (one propagating to the left, the other to the right) and moves off the screen. Additional oscillations, numerical garbage, or incomprehensible plots are an indication that mistakes were made.
Alternatively, you can punish yourself with something like gnuplot. Instructions are left as an exercise for the reader.