Step 1 — From the Mathematical Problem to the Solver
Before adding OpenMP, MPI, or CUDA, we need a numerical problem that is well defined and a sequential solver whose behavior we understand. The course handout develops this part mainly from the mathematical side. Here we will follow the same formulation while connecting each step to the objects and operations that eventually appear in the program.
The complete path is
\[ \text{Poisson equation} \longrightarrow \text{fictitious domain} \longrightarrow \text{finite-difference operator} \longrightarrow Aw=B \longrightarrow \text{preconditioned conjugate gradient}. \]Each transformation has a practical purpose. By the end of this section, the original PDE will have become a small collection of array operations, stencil evaluations, vector updates, and reductions. Those operations are also the units that we will parallelize later.
1. Define the physical domain
The original problem is a two-dimensional Poisson equation
\[ -\Delta u = f(x,y), \qquad (x,y)\in D, \]with homogeneous Dirichlet boundary conditions
\[ u(x,y)=0,\qquad (x,y)\in\partial D. \]The assignment fixes
\[ f(x,y)=1 \]and provides several possible geometries for \(D\). The reference implementation uses the trapezoidal variant with vertices
\[ A=(0,0),\qquad B=(3,0),\qquad C=(2,3),\qquad D=(0,3). \]This is variant 4 in the assignment.
In code, the geometry can be represented by a function that answers one simple question: does the point \((x,y)\) belong to the physical domain?
bool region(double x, double y) {
return (x > 0 && y > 0 && y < 3 && (y + 3 * x < 9));
}
double f_func(double x, double y) {
return 1.0;
}
The inequality
\[ y+3x<9 \]describes the slanted edge connecting \((3,0)\) and \((2,3)\). Together with the remaining bounds, it gives exactly the trapezoidal region used by the solver.
This small region() function becomes surprisingly important. The numerical code does not need special cases for triangles, trapezoids, ellipses, or other geometries throughout the solver. It only needs a way to evaluate whether a point lies inside \(D\). A different assignment variant can therefore reuse most of the numerical solver after replacing the geometry function.
2. Embed the domain into a rectangular computational region
The trapezoid is inconvenient for a regular finite-difference grid. The course therefore introduces the fictitious-domain method.
We place the physical region \(D\) inside a rectangular computational domain
\[ \Pi=(0,3)\times(0,3), \]and define the part of this rectangle outside the original geometry as the fictitious domain
\[ \widehat D=\Pi\setminus D. \]The computation can now use a regular rectangular grid across the whole of \(\Pi\). The irregular geometry is encoded through coefficients instead of through the shape of the grid. This is the main computational advantage of the fictitious-domain construction.
The modified problem introduces a piecewise coefficient
\[ k(x,y)= \begin{cases} 1, & (x,y)\in D,\\ \dfrac{1}{\varepsilon}, & (x,y)\in\widehat D, \end{cases} \]and a modified right-hand side
\[ F(x,y)= \begin{cases} f(x,y), & (x,y)\in D,\\ 0, & (x,y)\in\widehat D. \end{cases} \]The problem solved over the rectangle becomes
\[ -\frac{\partial}{\partial x} \left( k(x,y)\frac{\partial v}{\partial x} \right) - \frac{\partial}{\partial y} \left( k(x,y)\frac{\partial v}{\partial y} \right) = F(x,y), \]with
\[ v=0 \]on the outer boundary of \(\Pi\).
The parameter \(\varepsilon\) determines the approximation introduced by the fictitious region. The handout states that the solution of the rectangular problem approximates the original solution in \(D\) with an error bounded proportionally to \(\varepsilon\). For this assignment, the prescribed choice is
\[ \varepsilon=h^2, \qquad h=\max(h_x,h_y). \]The reference implementation follows this choice directly:
hx = (x_max - x_min) / M;
hy = (y_max - y_min) / N;
eps = std::max(hx, hy) * std::max(hx, hy);
It then stores the coefficient and right-hand side on the grid as two arrays:
f[i][j] = inside_region(x, y) ? source_func(x, y) : 0.0;
k[i][j] = inside_region(x, y) ? 1.0 : 1.0 / eps;
At this point the geometry has been reduced to numerical data. Every grid point carries a source value \(F\) and a coefficient \(k\), and the remaining solver can operate on a rectangular array.
3. Put a finite-difference grid over the rectangle
We divide the computational rectangle into an \(M\times N\) grid,
\[ x_i=A_1+i h_x, \qquad y_j=A_2+j h_y, \]where
\[ h_x=\frac{B_1-A_1}{M}, \qquad h_y=\frac{B_2-A_2}{N}. \]The unknown continuous function \(v(x,y)\) is represented by grid values
\[ w_{ij}\approx v(x_i,y_j). \]Only the interior points
\[ i=1,\ldots,M-1, \qquad j=1,\ldots,N-1 \]are unknown. The values on the outer boundary remain zero because of the Dirichlet condition. The course handout then expresses the discrete problem in operator form,
\[ Aw=B. \]This notation is convenient mathematically, although the implementation does not need to materialize \(A\) as one enormous matrix. We only need an operation that evaluates \(Aw\) for a supplied grid function \(w\).
For the variable-coefficient operator, the handout uses coefficients on grid-cell edges. In the \(x\) direction,
\[ a_{ij} = \frac{1}{h_y} \int_{y_{j-1/2}}^{y_{j+1/2}} k(x_{i-1/2},t)\,dt, \]and in the \(y\) direction,
\[ b_{ij} = \frac{1}{h_x} \int_{x_{i-1/2}}^{x_{i+1/2}} k(t,y_{j-1/2})\,dt. \]The resulting finite-difference equation has the form
\[ -\frac{ a_{i+1,j}(w_{i+1,j}-w_{ij}) - a_{ij}(w_{ij}-w_{i-1,j}) }{h_x^2} - \frac{ b_{i,j+1}(w_{i,j+1}-w_{ij}) - b_{ij}(w_{ij}-w_{i,j-1}) }{h_y^2} = F_{ij}. \]This is a local stencil calculation: the value at \((i,j)\) depends on the center point and its four immediate neighbors.
That locality will matter throughout the rest of the project. OpenMP can distribute different grid points among threads. MPI will eventually place neighboring values on different processes and require halo exchange. CUDA can assign these stencil evaluations to GPU threads.
4. Handle cells close to the physical boundary
The boundary between \(D\) and \(\widehat D\) can cut through a grid cell. The mathematical formulation pays special attention to these cells.
For the coefficients \(a_{ij}\) and \(b_{ij}\), the handout asks us to integrate the piecewise-constant coefficient \(k\) analytically along the corresponding grid edges. If an edge lies entirely inside \(D\), its effective coefficient is \(1\). If it lies entirely inside the fictitious domain, its coefficient is \(1/\varepsilon\). If the boundary crosses the edge, the coefficient depends on the length of the part lying in each region.
The same idea appears in the right-hand side. Formally,
\[ F_{ij} = \frac{1}{h_xh_y} \iint_{\Pi_{ij}}F(x,y)\,dx\,dy. \]For a cell completely inside \(D\), the handout permits using the value of \(f\) at the grid point. For a cell intersected by the physical boundary, the approximation uses the area
\[ S_{ij} = \operatorname{mes}(\Pi_{ij}\cap D). \]There is an important implementation detail here.
The reference implementation uses a simpler discretization. It stores \(k\) directly at grid nodes and evaluates edge coefficients with arithmetic averages:
double kx_plus = 0.5 * (k[i][j] + k[i+1][j]);
double kx_minus = 0.5 * (k[i][j] + k[i-1][j]);
double ky_plus = 0.5 * (k[i][j] + k[i][j+1]);
double ky_minus = 0.5 * (k[i][j] + k[i][j-1]);
The source term is also sampled according to whether the grid point itself belongs to region().
So there are two levels worth keeping separate when following this project:
| Part | Course formulation | Current SM25 implementation |
|---|---|---|
| \(k\) near the boundary | analytic edge integration | nodal \(k\) with arithmetic edge averages |
| \(F_{ij}\) near the boundary | cell integral / intersection area | value at the grid point |
| operator | variable-coefficient finite difference | matrix-free variable-coefficient stencil |
The reference code therefore captures the same fictitious-domain structure and parallel workload, while using a simpler treatment of cut cells. If the assignment is being implemented literally from the handout, the analytical edge coefficients and cell-area treatment should be implemented as specified there.
5. Implement \(A\) as an operator
Once the coefficients have been initialized, the core numerical operation becomes
\[ v\longmapsto Av. \]The reference implementation places this operation in apply_A():
out[i][j] = -1.0 * (
(kx_plus * (v[i+1][j] - v[i][j])
- kx_minus * (v[i][j] - v[i-1][j])) / hx2
+
(ky_plus * (v[i][j+1] - v[i][j])
- ky_minus * (v[i][j] - v[i][j-1])) / hy2
);
where hx2 = h_x^2 and hy2 = h_y^2 are the squared spacings precomputed by
the solver.
This matrix-free organization is useful for the entire project. We store the vectors needed by the iterative method while evaluating the action of \(A\) only when required.
The solver keeps arrays for
\[ u,\quad f,\quad r,\quad z,\quad p,\quad Ap, \]along with \(k\) and the diagonal preconditioner.
These objects have direct numerical meanings:
\[ u \]is the current approximation to the solution,
\[ r=B-Au \]is the residual,
\[ z=M^{-1}r \]is the preconditioned residual,
\[ p \]is the current conjugate search direction, and
\[ Ap \]is the operator applied to that direction.
Once the implementation reaches this form, the mathematical problem has effectively become a collection of structured array kernels.
6. Add diagonal preconditioning
The handout solves the linear system with the conjugate-gradient method and introduces diagonal preconditioning to improve convergence.
The diagonal of the finite-difference operator at one grid point can be written from the coefficients surrounding that point. The reference implementation computes
\[ d_{ij} = \frac{k_{x+}+k_{x-}}{h_x^2} + \frac{k_{y+}+k_{y-}}{h_y^2} \]and stores
\[ M^{-1}_{ij}=\frac{1}{d_{ij}}. \]In code:
double diag =
(kx_plus + kx_minus) / hx2
+ (ky_plus + ky_minus) / hy2;
M_inv[i][j] = 1.0 / diag;
Applying the preconditioner then becomes a pointwise operation,
z[i][j] = M_inv[i][j] * r[i][j];
This is computationally simple and particularly convenient for later parallelization because every point can perform the multiplication independently.
7. Turn conjugate gradient into a sequence of kernels
With \(A\) and the preconditioner available, we can assemble the actual iterative solver.
The initial approximation in the reference implementation is
\[ u^{(0)}=0, \]because all solution arrays are initialized to zero. We first compute
\[ r^{(0)}=B-Au^{(0)}, \]then apply the preconditioner,
\[ z^{(0)}=M^{-1}r^{(0)}, \]and choose
\[ p^{(0)}=z^{(0)}. \]The implementation then repeatedly performs the following PCG update.
First compute
\[ Ap^{(k)}. \]Then evaluate
\[ \alpha_k = \frac{(r^{(k)},z^{(k)})} {(p^{(k)},Ap^{(k)})}. \]Update the solution,
\[ u^{(k+1)} = u^{(k)}+\alpha_k p^{(k)}, \]and the residual,
\[ r^{(k+1)} = r^{(k)}-\alpha_k Ap^{(k)}. \]Apply the diagonal preconditioner again,
\[ z^{(k+1)} = M^{-1}r^{(k+1)}, \]compute
\[ \beta_k = \frac{(r^{(k+1)},z^{(k+1)})} {(r^{(k)},z^{(k)})}, \]and obtain the next search direction,
\[ p^{(k+1)} = z^{(k+1)}+\beta_k p^{(k)}. \]This loop maps almost one-to-one to the source:
apply_A(p, A_p);
double p_Ap = compute_p_Ap();
double alpha = rz_prev / p_Ap;
update_u(alpha);
update_r(alpha);
r_norm = compute_l2_norm();
apply_preconditioner();
double rz = compute_rz();
double beta = rz / rz_prev;
update_p(beta);
rz_prev = rz;
This decomposition is one of the most useful things to understand before moving to parallel programming. The solver is built from three recurring computational patterns:
Stencil operations appear in apply_A.
Pointwise vector updates appear in update_u, update_r, update_p, and apply_preconditioner.
Global reductions appear in inner products such as
\[ (p,Ap) \]and
\[ (r,z), \]as well as in the residual norm.
These three patterns will behave very differently once the program moves from shared memory to distributed memory and then to GPUs.
8. Decide when the iteration has converged
The course notes propose monitoring convergence through the change between successive approximations,
\[ \left\| w^{(k+1)}-w^{(k)} \right\|_E<\delta, \]and also discuss monitoring the behavior of the associated functional to detect numerical deterioration.
The reference implementation uses a different practical stopping test. It computes the Euclidean norm of the residual,
\[ \|r\|_2 = \sqrt{ \sum_{i,j}r_{ij}^2 }, \]and continues while
\[ \|r\|_2>\text{tolerance}. \]The default tolerance in solve() is \(10^{-4}\), with a maximum of 100,000 iterations. The residual is periodically written to debug.log, making it possible to inspect convergence during larger runs.
Again, this is a point where the reference implementation and the handout use different stopping criteria. A course submission should make the chosen criterion explicit and keep it consistent when comparing sequential and parallel versions.
9. From the Mathematics to the Program
At this point, the mathematical formulation has given us all of the pieces needed to build the first solver.
The physical domain \(D\) defines the geometry of the problem. The fictitious-domain construction embeds that geometry into a rectangular computational region \(\Pi\), allowing us to work on a regular grid. The finite-difference scheme then replaces the differential equation with a discrete operator
\[ Aw = B, \]and the preconditioned conjugate-gradient method gives us an iterative way to solve that system without explicitly constructing the full matrix \(A\).
For implementation, these mathematical objects translate into a small set of data structures and computational operations:
| Mathematics | Program |
|---|---|
| computational grid | \(M \times N\) arrays |
| physical domain \(D\) | region(x, y) |
| right-hand side \(F\) | f |
| fictitious-domain coefficient \(k\) | k |
| approximate solution \(w\) | u |
| discrete operator \(A\) | apply_A() |
| residual \(B-Au\) | r |
| preconditioned residual | z |
| conjugate direction | p |
| \(Ap\) | A_p |
| diagonal preconditioner | M_inv |
The computation can also be viewed in terms of the three recurring patterns from Section 7 — stencil operations, pointwise updates, and reductions.
Together, these operations form the computational core of the solver:
\[ \text{initialize the grid} \rightarrow \text{build } F \text{ and } k \rightarrow \text{apply } A \rightarrow \text{update the PCG vectors} \rightarrow \text{check convergence}. \]There is still an important step between understanding this algorithm and running it on multiple threads or processes. We first need to assemble these pieces into a complete sequential program and verify that it produces a stable numerical solution.
That sequential implementation will serve as the reference for every later version of the assignment. It gives us a known execution order, a numerical result to compare against, and a baseline runtime for measuring speedup.
In the next step, we will build that reference solver from the ground up: organize its state, initialize the grid and fictitious domain, implement the finite-difference operator, assemble the PCG iteration, save the solution, and check that the program behaves correctly on progressively finer grids.
