Using Eigen to Solve the Pressure Poisson Equation
Published:
When developing Computational Fluid Dynamics (CFD) solvers, one frequently encounters the pressure Poisson equation. This equation typically results in a large sparse system of linear equations (Ax=b) that needs to be solved efficiently. In this post, I’ll explain how to use the Eigen C++ library to solve such systems, particularly when your matrix is stored in the Compressed Row Storage (CRS) format.
Why BiCGSTAB?
The pressure Poisson equation typically results in a matrix that is:
- Sparse (most elements are zero)
- Potentially asymmetric (due to boundary conditions)
- Large (especially for 3D problems)
While the basic Poisson operator is symmetric, the incorporation of boundary conditions often introduces asymmetry. For this reason, we use the BiCGSTAB (Biconjugate Gradient Stabilized) solver, which can handle asymmetric matrices efficiently.
Code Structure
Let’s break down the solution process into several key components:
1. Required Headers and Setup
#include <Eigen/Sparse>
#include <iostream>
#include <vector>
#include <fstream>
#include <string>
using namespace Eigen;
This section includes the necessary Eigen headers for sparse matrix operations, along with standard C++ headers for file I/O.
2. File Reading Utility
template<typename T>
std::vector<T> readVectorFromFile(const std::string& filename) {
std::vector<T> data;
std::ifstream file(filename);
T value;
if (!file.is_open()) {
std::cerr << "Failed to open file: " << filename << std::endl;
exit(1);
}
while (file >> value) {
data.push_back(value);
}
file.close();
return data;
}
This templated function reads numerical data from files, handling both integer (for indices) and double (for values) inputs. It’s used to read the CRS format components and the right-hand side vector.
3. Matrix Assembly
// Read matrix dimensions
int rows, cols;
std::ifstream dimFile("dimensions.txt");
dimFile >> rows >> cols;
// Read CRS format data
std::vector<double> values = readVectorFromFile<double>("values.txt");
std::vector<int> column_indices = readVectorFromFile<int>("column_indices.txt");
std::vector<int> row_pointers = readVectorFromFile<int>("row_pointers.txt");
std::vector<double> b = readVectorFromFile<double>("rhs.txt");
This section reads the matrix data stored in CRS format. The data is split across multiple files for better organization:
- Matrix dimensions
- Non-zero values
- Column indices for each value
- Row pointers indicating where each row starts
- Right-hand side vector
4. Converting to Eigen’s Format
typedef Triplet<double> T;
std::vector<T> tripletList;
tripletList.reserve(values.size());
for(int i = 0; i < rows; i++) {
for(int j = row_pointers[i]; j < row_pointers[i+1]; j++) {
tripletList.push_back(T(i, column_indices[j], values[j]));
}
}
SparseMatrix<double> A(rows, cols);
A.setFromTriplets(tripletList.begin(), tripletList.end());
A.makeCompressed();
Eigen uses a triplet format for matrix construction. This section converts our CRS data into Eigen’s format and creates the sparse matrix.
5. Solver Setup and Solution
BiCGSTAB<SparseMatrix<double>> solver;
solver.setTolerance(1e-6);
solver.setMaxIterations(1000);
solver.compute(A);
VectorXd x = solver.solve(b_eigen);
Here we:
- Create the BiCGSTAB solver
- Set the convergence tolerance
- Set maximum iterations
- Compute the matrix factorization
- Solve the system
6. Solution Output and Verification
std::cout << "Number of iterations: " << solver.iterations() << std::endl;
std::cout << "Estimated error: " << solver.error() << std::endl;
std::ofstream solFile("solution.txt");
solFile.precision(12);
for(int i = 0; i < x.size(); i++) {
solFile << x[i] << "\n";
}
double residual = (A * x - b_eigen).norm();
std::cout << "Residual: " << residual << std::endl;
Finally, we:
- Print solver statistics
- Save the solution to a file
- Compute and print the residual norm for verification
Performance Considerations
- Matrix Construction
- Using
reserve()
for the triplet list prevents reallocations makeCompressed()
ensures efficient storage
- Using
- Solver Parameters
- Tolerance affects accuracy vs. speed
- Maximum iterations prevents infinite loops
- Consider adjusting these based on your specific problem
Conclusion
This approach provides an efficient way to solve the pressure Poisson equation using Eigen’s sparse solvers. The code is designed to work with external data in CRS format, making it easy to integrate with existing CFD codes. The BiCGSTAB solver handles the potential asymmetry from boundary conditions, while maintaining good performance for large sparse systems.
Remember to compile with proper flags:
g++ -I /usr/include/eigen3 poisson_solver.cpp -o poisson_solver
Leave a Comment