A domain decomposition strategy and a parallel gradient‐type iterative solution scheme have been developed and implemented for the computation of complex three‐dimensional viscous flow problems involving heat transfer and surface tension effects. Special attention has been paid to the kernels for the computationally intensive matrix–vector products and dot products, to memory management, and to overlapping communication and computation. Details of these implementation issues are described together with associated performance and scalability studies. Representative Rayleigh–Bénard and microgravity Marangoni flow calculations and performance results on the Cray T3D and T3E are presented. Performance studies have been recently carried out and sustained rates above 50 gigaflops and 100 gigaflops have been achieved on the 512‐node T3E‐600 and 1024‐node T3E‐900 configurations respectively. The work is currently being extended to tightly‐coupled parallel ‘Beowulf‐type’ PC clusters and some preliminary performance results on this platform are presented. Copyright © 1999 John Wiley & Sons, Ltd.