Viscoelastic fluid flow with immersed boundaries of complex geometry is widely found both in nature and engineering processes. Examples include haemocytes moving in human blood flow, self-propulsion of microscopic organisms in complex liquids, hydraulic fracturing with sand in oil flow, and suspension flow with viscoelastic medium. Computational modelling of such systems is important for understanding complex biological processes and assisting engineering designs. Conventional simulation methods use conformed meshes to resolve the boundaries of complex geometry. Dynamically updating the conformed mesh is computationally expensive and makes parallelization difficult. In comparison, Cartesian grid methods are more promising for large scale parallel simulation. Using Cartesian grid methods to simulate viscoelastic fluid flow with complex boundaries is a relatively unexplored area.In this thesis, a sharp interface Cartesian grid method (SICG) and a smoothed interface immersed boundary method (SIIB) are developed in order to simulate viscoelastic fluids in complex geometries. The SICG method shows a better prediction of the stress on stationary boundaries while the SIIB method shows reduced non-physical oscillations in the computation of drag and lift forces on moving boundaries. Parallel implementations of both solvers are developed. Convergence of the numerical schemes is shown and the implementations are validated with a few benchmark problems with both stationary and moving boundaries.This study also focuses on the simulation of flows past 2D cylindrical or 3D spherical particles. Lateral migration of particles induced by inertial and viscoelastic effects are investigated with different flow types. Equilibrium positions of inertia-induced migration are reported as a function of the particle Reynolds number and the blockage ratio. The migration in the viscoelastic fluid is simulated from zero elastic number to a finite elastic number. The inclusion of both inertial and viscoelastic effects on the lateral migration of a particle is the first of its kind. New findings are reported for the equilibrium positions of a spherical particle in square duct flow, which suggest the need for future experimental and computational work.