r/ScientificComputing • u/TheMakpu • Apr 29 '26

Heat2D: a C++ heat equation solver in 2D

https://github.com/paurierap/heat-2D

Hi all!

I recently built a small C++ project to solve the generalized heat equation in 2D rectangular domains, mainly as a way to better understand PDE solver design and experiment with separating between spatial discretization and time integration.

The architecture of the project is modular:

- Spatial discretization (currently finite differences with Eigen sparse matrices)

- Time integration (Explicit/Implicit Euler, Crank–Nicolson)

The idea is to make it easy to plug different methods in each part (I’m planning to add finite elements soon). There are also a few example simulations (moving heat source, colliding pulses, etc.) with GIFs in the README.

I’d really appreciate some feedback, especially on:

- The overall design/architecture

- Performance considerations

- Whether this approach makes sense vs more template-heavy designs

28 Upvotes

100% Upvoted

3 comments sorted by

5

u/glvz Apr 29 '26

You could look into GridAp Written in Julia to get a sense of the architecture. This looks nice. Have you profiled the code? What's the bottleneck

1

u/TheMakpu Apr 29 '26

Funny you should mention that, some years ago I worked on GridapMakie during Google Summer of Code, but never really got into the design of it…

In terms of profiling, and speaking about the Crank-Nicolson integrator only, the bottleneck appears when solving the linear system at every timestep. I make sure to cache the factorization of the LHS matrix and use an “appropriate” solver for that (whether it is SPD and conditioning number is low enough for an iterative solver), but the performance depends heavily on both space resolution and the timestep. There is a benchmark executable that provides time measurements for the factorization and every timestep too!

3

u/Organic-Scratch109 Apr 29 '26

There many different approaches to speed up the linear solver for this problem:

  • Use an iterative solver (cg if the matrix is SPD). This is not necessary for two d problems, but very necessary in 3d due to the curse of dimensionality.
  • Good preconditioner: Geometric multigrid methods are easy to code in rectangular domains.
  • Good stopping criteria (you don't need 1e-8 error from your linear solver if the discretization error is 1e-3).
  • Alternatively, look into exponential integrators with Krylov subspaces.