In ocean engineering, numerical problems are often large-scale and involve highly complex hydrodynamic phenomena with free surfaces, which pose significant challenges in terms of robustness, computational efficiency, and virtual memory demands due to the extensive mesh requirements. To address the above issues, this work proposes a numerical method based on the lattice Boltzmann flux solver (LBFS) for multiphase flows. This method employs a specialized memory optimization strategy and enables data exchange between different GPUs, efficient parallel computation on multi-GPU platforms is achieved thereby enabling the simulation of flow fields with meshes exceeding hundreds of millions. Meanwhile, by applying the immersed boundary method, the motion of rigid structures can span flow fields stored on different GPUs, realizing efficient interactions between the solid and fluid. The computational efficiency of the single-GPU and multi-GPU parallel computing versions is compared in this work. The validation shows that the multi-GPU numerical method exhibits outstanding parallel efficiency, and the accuracy of algorithm is verified by the simple of water entry of a cylinder. Finally, this method is used to simulate the rolling and advancing motions of a small boat, analyzing both the details of flow field and the time evolution of the boat's motion. In conclusion, this method has demonstrated superior computational efficiency and memory usage efficiency in simulating large-scale numerical problems.