## 2 Methodology
### Calculation of RTM SOGs with Multi-GPUs  
In the RTM calculation, the kernel of the wavefield propagation is a finite-difference calculation. Among them, the best GPU calculation strategy for spatial difference has been achieved in previous studies with a single GPU [25]; it can achieve nearly 100% occupancy, defined as the ratio of active warps on an SM to the maximum number of active warps supported by the SM. However, with an increasing demand on computation, the memory of a single GPU cannot meet the requirements for larger storage space. A GPU enables two solutions: (1) it provides a unified memory space in which a memory pool is built using larger CPU memory, thereby enabling data in the memory pool to be read by the CPU and GPU simultaneouslysolving the demand for large memory, but reading the memory slowly; (2) it introduces a multi-GPU calculation, which achieves the common calculation by dividing the computational domain into different sub-cells, placing them in different GPUs within the same network for calculation, and exchanging data in the overlapped domain when necessary. We used two GPUs as an example and focused on expounding the second solution. In the calculation, the wavefield space to be calculated was divided into two mutually overlapped subspaces, in the z direction, and they were calculated on two GPUs, respectively. For the single GPUs in the multi-GPU calculation, the finite-difference calculation was still based on the conventional methods, and the only distinction was that after the finite difference of each time step was completed, data exchange in the overlapping area was needed. The difference calculation space was divided in two GPUs as shown in Figure 6, in which Figure 6a shows the entire computational space, and Figure 6b shows the computational space partitioned into two equal subspaces. However, based on the order in the finite-difference calculation, an overlapped domain was needed for data exchange and updating, as the green part shows in Figure 6b. Figure 6c is a slice along the z direction, in which the red is the computational domain, and the green is the overlapped domain for data exchange. Figure 5: The process for RTM-based calculation to output SOGs. In terms of memory management methods, for multi-GPU computations, the computational process can be illustrated as shown in Figure 7. The entire imaging space was divided into two independent parts in the z direction, which were calculated with GPU0 and GPU1, respectively. Some of the overlapped domain was used to exchange boundary data at different times of propagation, and the radius was determined by the finite-difference order. For each GPU computation, we focused on the management and calculation of data in the overlapped domain, which involved data transmission between two GPUs. To realize the overlap of calculation and data transmission, it was necessary to adopt different methods when dividing the computational domain. For example, for GPU1 computations, the computational domain partitioned to it was further divided into three parts, of which, d_wf11 and d_wf21 were calculated in different job streams, respectively, and after d_wf20 was computed, it was directly exchanged to the corresponding domain of GPU1 through point-to-point transmission, and GPU1 was calculated using the same process as for GPU0. Here, the magnitude of d_wf20 and d_wf11 was by the order of spatial finite difference and was much smaller than that of d_wf10 and d_wf20. If the computation and transmission of a small domain on a single GPU are placed in the same stream, when the computation is not completed in the larger domain, the computation and transmission are completed in the small area, so the transmission of data between different GPUs is completely concealed.