The wind flow through an urban built environment has a significant impact on the safety and comfort of pedestrians and inhabitants. Simulation of urban wind flow presents a formidable challenge to computational engineering field due to the complex geometry of the built environment and the multiscale nature of the flow. Numerical analysis via computational fluid dynamics (CFD) based on the Navier-Stokes equations is now commonplace in the industry, usually run on high performance computer clusters of CPU nodes. However, resolution of all the turbulent scales of motion in the entire domain is likely to remain beyond the reach of such hardware for the foreseable future; thus, so called Direct Numerical Simulation is not practical for industrial applications. Moreover, the region of interest usually represents a small percentage of the total volume of the domain. The objective of this thesis is to achieve a time-dependent simulation that resolves large to medium turbulence scales in the region of interest at a significantly reduced computational cost compared to turbulent scale resolving Navier-Stokes methods. To do so, we couple a lattice Boltzmann (LB) solver running on graphics processing units (GPUs) with a Navier-Stokes (NS) solver running on CPUs. The LB solver incorporates the mean turbulent flow information from the Navier-Stokes model into its resolved velocity via a synthetic eddy method (SEM) implemented at the inlet of the LB domain. The resulting coupled Navier-Stokes lattice Boltzmann (NSLB) solver combines the accuracy and computing speed of the GPU implementation of the LB model with the stability, low memory consumption and mesh flexibility of the NS solver. Moreover, the NSLB model exploits the widespread availability of CPU and GPU hardware on desktop, and workstation computers. Validation results of the two-way coupled NSLB solver for laminar flow demonstrate that the coupled solver is able to reproduce the results of the full domain single solver (either LB or NS) independently of the position of the interface between the solvers. For the application to urban wind flow, we embedded our lattice Boltzmann large eddy simulation (LB-LES) solver within a pre-calculated Reynolds averaged Navier-Stokes (RANS) simulation of flow around a single building at $Re_H = 47893$. The LB-LES model accurately predicts the mean and fluctuating velocity around the building, increasing the flow data available and its accuracy with respect to the underlying RANS simulation. The RANS LB-LES coupling allows fully resolved LES results to be achieved in practical time scales with a single desktop based GPU, which shows potential to run industry applications on consumer devices. The original contributions of this work include the development of the coupling framework, adaptation of the SEM to LB and refinement of the embedded boundary conditions. The resulting solver fully realises the objective of a coupled Navier-Stokes method with a GPU accelerated lattice Boltzmann method in order to accurately simulate high Reynolds number flow at affordable computational cost.