C-Poisson-with-PETSc

Repository of FEM Poisson MMS code for 2D and 3D

Q2 Poisson (2D/3D) — MATLAB & PETSc/MPI

This repo contains compact, reproducible Q2 (quadratic) finite-element solvers for Poisson with zero Dirichlet BCs and manufactured sine solutions.

• MATLAB: clean reference, matrix-free(ish) application of the stiffness using sum-factorization; assembles only the RHS.
• C++/PETSc/MPI: parallel assembled FEM + CG; robust true (L^2) error via Gauss quadrature (numerically stable at fine meshes). 3D code also prints region timings.

---

Repo layout

matlab/
  FEM2D_sumfact2.m          # 2D reference on [0,π]×[0,π], m=n=3 by default
  FEM3D_sumfact2.m          # 3D reference on [0,1]^3, mx=2,my=3,mz=4 by default

cpp/
  poisson_q2_petsc_cg.cpp         # 2D PETSc/MPI (CG; robust L2)
  poisson_q2_3d_petsc_cg.cpp      # 3D PETSc/MPI (CG; robust L2 + timings)

If your MATLAB filenames differ, just adjust the paths; the descriptions still apply.

---

Requirements

MATLAB

• MATLAB R2018a+ (no toolboxes required).

C++/PETSc

• MPI toolchain (mpicxx, mpirun).

• PETSc with C/C++: tested with PETSc 3.15.5 (Debian/Ubuntu packages). Also works on newer PETSc (3.18–3.20+) with the options shown below.

• Optional PETSc features (if you want them):

  • Hypre BoomerAMG: PETSc built --with-hypre=1.
  • MUMPS: PETSc built --with-mumps=1.

On Debian/Ubuntu PETSc 3.15, use pkg-config --cflags --libs petsc. If you built PETSc yourself, set PETSC_DIR/PETSC_ARCH.

---

Build (C++)

From cpp/:

# 2D
mpicxx -std=c++11 poisson_q2_petsc_cg.cpp \
  $(pkg-config --cflags --libs petsc) -o poisson_q2_2d

# 3D
mpicxx -std=c++11 poisson_q2_3d_petsc_cg.cpp \
  $(pkg-config --cflags --libs petsc) -o poisson_q2_3d

If pkg-config can’t find PETSc:

export PETSC_DIR=/path/to/petsc
export PETSC_ARCH=arch-linux-c-opt   # your PETSc build arch name
mpicxx -std=c++11 -I$PETSC_DIR/include -I$PETSC_DIR/$PETSC_ARCH/include \
  poisson_q2_3d_petsc_cg.cpp \
  -L$PETSC_DIR/$PETSC_ARCH/lib -lpetsc -o poisson_q2_3d

---

Running tips

• Pin ranks & avoid threaded BLAS inside MPI ranks:

  export OMP_NUM_THREADS=1
  # Open MPI (example: 15 hwthreads on an 8c/16t CPU)
  mpirun -n 15 --use-hwthread-cpus --bind-to hwthread --map-by hwthread <program> ...

• Solver tolerance rule of thumb for Q2 Poisson: -ksp_rtol ≈ h^2. (E.g., for Nx=400 in 2D, h≈π/400 ⇒ h^2≈6e-5; so 1e-6…1e-5 is plenty.)

---

2D (C++) — poisson_q2_petsc_cg.cpp

Problem: (\Omega=[0,\pi]^2), (u=\sin(mx,x)\sin(ny,y)), default (m=n=3). Unknowns: ((2N_x+1)\times(2N_y+1)) nodes (Q2).

Match MATLAB (no preconditioner, true residual)

mpirun -n 15 --use-hwthread-cpus --bind-to hwthread --map-by hwthread \
  ./poisson_q2_2d -Nx 400 -Ny 400 \
  -ksp_type cg -pc_type none -ksp_norm_type unpreconditioned -ksp_rtol 1e-6 \
  -ksp_converged_reason -ksp_monitor_true_residual

Faster with Hypre BoomerAMG (recommended for larger, general cases)

mpirun -n 15 --use-hwthread-cpus --bind-to hwthread --map-by hwthread \
  ./poisson_q2_2d -Nx 400 -Ny 400 \
  -ksp_type cg -ksp_rtol 1e-8 \
  -pc_type hypre -pc_hypre_type boomeramg \
  -pc_hypre_boomeramg_coarsen_type HMIS \
  -pc_hypre_boomeramg_interp_type ext+i \
  -pc_hypre_boomeramg_strong_threshold 0.5 \
  -pc_hypre_boomeramg_max_iter 1 \
  -pc_hypre_boomeramg_relax_type_all symmetric-SOR/Jacobi \
  -ksp_cg_single_reduction -ksp_converged_reason

Direct solve (ground truth, smaller N)

./poisson_q2_2d -Nx 120 -Ny 120 \
  -ksp_type preonly -pc_type lu -pc_factor_mat_solver_type mumps

---

3D (C++) — poisson_q2_3d_petsc_cg.cpp

Problem: (\Omega=[0,1]^3), (u=\sin(m_x\pi x)\sin(m_y\pi y)\sin(m_z\pi z)), default (m_x=2, m_y=3, m_z=4). Unknowns: ((2N_x+1)(2N_y+1)(2N_z+1)) nodes (Q2).

Unpreconditioned (to mirror MATLAB), with region timings

mpirun -n 15 --use-hwthread-cpus --bind-to hwthread --map-by hwthread \
  ./poisson_q2_3d -Nx 20 -Ny 20 -Nz 20 \
  -ksp_type cg -pc_type none -ksp_norm_type unpreconditioned -ksp_rtol 1e-6 \
  -ksp_converged_reason -ksp_monitor_true_residual

The program prints:

TIMINGS (s): Assembly ...  Solve ...  L2 ...  Total ...

With Hypre BoomerAMG (fast & robust)

mpirun -n 15 --use-hwthread-cpus --bind-to hwthread --map-by hwthread \
  ./poisson_q2_3d -Nx 40 -Ny 40 -Nz 40 \
  -ksp_type cg -ksp_rtol 1e-8 \
  -pc_type hypre -pc_hypre_type boomeramg \
  -pc_hypre_boomeramg_coarsen_type HMIS \
  -pc_hypre_boomeramg_interp_type ext+i \
  -pc_hypre_boomeramg_strong_threshold 0.5 \
  -pc_hypre_boomeramg_max_iter 1 \
  -pc_hypre_boomeramg_relax_type_all symmetric-SOR/Jacobi \
  -ksp_cg_single_reduction -ksp_converged_reason

---

(L^2) error (C++): what we compute & why it’s stable

• We compute the true (|u_h-u|_{L^2}) by direct Gauss quadrature:

  • 2D: 3×3 points/element; 3D: 3×3×3 — exact for Q2 on affine quads/hexes.

• To be MPI-robust and numerically stable at fine meshes:

  • Each rank performs one-shot VecScatter to gather all needed nodal values (no VecGetValues in loops).
  • Accumulation of ((u_h-u)^2) uses long double to avoid cancellation/NaNs.

The reported (L^2) matches MATLAB to within roundoff.

---

Preconditioner notes & options

• None (-pc_type none): often converges in a handful of steps for this sine MMS; good for reproducing MATLAB behavior.
• Hypre BoomerAMG: excellent for Poisson; use the symmetric relaxation shown above to keep CG safe (SPD).
• GAMG: available without external deps; for Q2 Poisson, consider providing a near-nullspace in code (MatSetNearNullSpace) for best results.
• Direct (MUMPS): convenient for verification on smaller meshes.

---

Version specifics (important)

• Primary target: PETSc 3.15.5 (Ubuntu/Debian packages). The following options are supported in this version and newer:

  • -ksp_norm_type unpreconditioned, -ksp_monitor_true_residual, -ksp_cg_single_reduction
  • Hypre BoomerAMG knobs used above (-pc_hypre_boomeramg_*)

• If you use a much newer PETSc (≥3.20), the same options work; behavior differences should be minor.

---

Profiling

• PETSc profiler:

  ... -log_view

• The 3D binary additionally prints concise wall-times for Assembly / Solve / L2 / Total.

---

License

MIT

---

If you spot a mismatch or want a preset like -match_matlab / -use_amg baked into the binaries, open an issue or PR. Happy solving! 🧮⚡