Here are my arguments for why functional programming can, and should be utilized for computational science. The benefits are vast, and the cons are quickly going away. In my mind there's only one con:
Con: lack of language support in C / C++ / Fortran
At least in C++, this con is vanishing - as C++14/17 has added powerful facilities to support functional programming. You might need to write some library / support code yourself, but the language will be your friend. As an example, here is a (warning: plug) library that does immutable multi-dimensional arrays in C++: https://github.com/jzrake/ndarray-v2.
Also, here is a link to a good book on functional programming in C++, although it's not focused on science applications.
Here is my summary of what I believe are the pro's:
Pros:
- Correctness
- Understandability
- Performance
In terms of correctness, functional programs are manifestly well-posed: they force you to properly define the minimal state of your physics variables, and the function that advances that state forward in time:
int main()
{
auto state = initial_condition();
while (should_continue(state))
{
state = advance(state);
side_effects(state);
}
return 0;
}
Solving a partial differential equation (or ODE) is perfect for functional programming; you're just applying a pure function (advance) to the current solution to generate the next one.
In my experience, physics simulation software is by and large, burdened by poor state management. Usually, each stage of the algorithm operates on some piece of a shared (effectively global) state. This makes it difficult, or even impossible, to ensure the correct order of operations, leaving the software vulnerable to bugs that can manifest as seg-faults, or worse, error terms that do not crash your code but silently compromise the integrity of its science output. Attempting to manage shared state in a physics simulation also inhibits multi-threading - which is a problem for the future, as supercomputers are moving toward higher core counts, and scaling with MPI often tops out at ~100k tasks. In contrast, functional programming makes shared-memory parallelism trivial, because of immutability.
Performance is also improved in functional programming due to the lazy evaluation of algorithms (in C++, this means generating many types at compile time - often one for each application of a function). But it reduces the overhead of memory accesses and allocations, as well as eliminating virtual dispatch - allowing the compiler to optimize an entire algorithm by seeing at once all the function objects comprising it. In practice, you'll experiment with different arrangements of the evaluation points (where the algorithm result is cached to a memory buffer) to optimize use of CPU vs. memory allocations. This is rather easy due to the high locality (see the example below) of the algorithm stages compared to what you'll typically see in a module or class-based code.
Functional programs are easier to understand insofar as they trivialize the physics state. That is not to say their syntax is readily understandable by all of your colleagues! Authors should be careful to use well-named functions, and researchers in general should get accustomed to seeing algorithms expressed functionally rather than procedurally. I'll admit that the absence of control structures can be off-putting to some, but I don't think that should stop us from going into the future able to do better quality science on computers.
Below is a sample advance function, adapted from a finite-volume code using the ndarray-v2 package. Note the to_shared operators - these are the evaluation points I was alluding to earlier.
auto advance(const solution_state_t& state)
{
auto dt = determine_time_step_size(state);
auto du = state.u
| divide(state.vertices | volume_from_vertices)
| nd::map(recover_primitive)
| extrapolate_boundary_on_axis(0)
| nd::to_shared()
| compute_intercell_flux(0)
| nd::to_shared()
| nd::difference_on_axis(0)
| nd::multiply(-dt * mara::make_area(1.0));
return solution_state_t {
state.time + dt,
state.iteration + 1,
state.vertices,
state.u + du | nd::to_shared() };
}
