A general formulation of the spherical harmonics (PN) methods was developed recently to expand the method to high orders of PN. The set of N(N + 1)/2 three-dimensional second-order elliptic PDEs formulation and their Marshak boundary conditions for arbitrary geometries are implemented in the openfoam finite volume based cfd software. The results are verified for four cases, including a 1D slab, a 2D square enclosure, a 3D cylindrical enclosure, and an axisymmetric flame. All cases have strongly varying radiative properties, and the results are compared with exact solutions and solutions from the photon Monte Carlo method (PMC).