Combustion remains a critical technology for electricity generation, heating, transportation and other industrial processes. Turbulent combustion lies at the heart of many of these processes. The accurate, robust and efficient computational modeling of turbulent combustion is necessary to design clean, efficient and safe combustion devices and processes. For practical combustion problems the direct numerical simulation (DNS) of the governing equations is computationally intractable. The Reynolds averaged Navier-Stokes simulation (RANS) and large eddy simulation (LES) techniques have emerged as powerful tools to simulate turbulent reacting flows. RANS and LES methodologies require closure of the unclosed terms arising from the averaging or filtering the governing equations. Even with adequate closure, RANS and LES remain computationally infeasible for simulating many combustion processes in engineering applications, further simplifications regarding flame thickness, flow and chemical reaction time scales are required. The high Damk\"{o}hler (Da) number flames can be modeled using a reduced chemistry model. A flamelet derived reduced chemistry model like Flamelet Generated Manifolds (FGM) accounts for finite rate chemistry while it greatly simplifies the simulation of turbulent combustion as it decouples the turbulent transport and flame structure. The interaction between the turbulence and the flame front in non-premixed combustion is described by the probability density function (PDF) of the composition variables. In this work, a framework for turbulent combustion modeling is presented for both RANS and, LES with FGM reduced chemistry model. This framework consists of implementation of presumed and transported PDF models and is developed within the open source CFD software OpenFOAM. The simulation of the well-known piloted methane-air jet flames (Sandia flames) is conducted in RANS context with both presumed and transported PDF models. An "A priori" analysis is conducted based on the RANS/TPDF simulation data. The analysis quantifies the extent of errors in PPDF models, specifically errors in choice of presumed PDF, statistical independence and the number moments and cross moments considered. A new PPDF model based on the Gaussian copula approach for correlation of the composition variables is developed and analyzed. The implementation of RANS/TPDF solver incorporates robust algorithms for particle tracking, position and number control. The LES/TFDF simulation of Sandia flame D is conducted to showcase the capability of the developed framework.