I don't think that to understand convolution you need to dig into the nested code of huge libraries, since the code becomes quickly really hard to understand and convoluted (ba dum tsss!). Joking apart, in PyTorch Conv2d is a layer that applies another low level function, conv2d, written in c++.
Luckily enough, the guys from PyTorch wrote the general idea of how convolution is implemented in the documentation:
From this paragraph, we already have some important information like the input and output dimensions. The number of channels should be easy to understand, if we have an RGB image, for example, the channels are 3, one for each color, so they are just different matrices representing different features.
The next important element is the reference to the cross correlation, the function applied to our input images through the kernel k. Why cross-correlation? Because it is almost identical to a convolution, as you can see comparing their formulas:
The only difference consists in the way the indexing of the width is implemented, which causes the operation to start from the bottom right of the input matrix for the convolution and from the top left of the input image from the cross-correlation (the circles in the squares in the previous pic). Since, in most programming languages, matrix indexing starts from top left, cross-correlation is the most common choice to implement.
But how do these formulas work in practice? Here's another picture taken from Chapter 9 of Deep Learning (Goodfellow, Bengio, Courville), which I strongly suggest you to read.
Basically, from the input matrix, a submatrix is extracted, with the same dimension of the kernel, then the sub-matrix and the kernel are multiplied elementwise and all the resulting product summed together to produce a single output element that will form a 'pixel' of the resulting feature map (the output matrix).
Here's another example with fake numbers that I made. I hope the double notation for filter/kernel doesn't generate confusion, I actually found that sometimes it is inconsistent (in the chapter I linked they don't even use the filter at all). In practice they actually mean the same thing, I usually call kernel the actual matrix that is multiplied to the input and with filter I refer to the sliding window on the input image (that has, of course, must have the same dimension of the kernel).
Lastly, when you apply padding, the filter can actually move also outside the 'edges' of the input matrix, in which case all the elements outside would be considered to be zero. The computation is exactly the same, but since there are more splitting steps, the output matrix will have bigger dimension.
Please note that with multiple input channels you can perform either 2d convolution or 3d convolution, the difference relies in the filter dimension, in case of 2d convolution it would be a square whereas in 3d convolution it would be a cube. This means that for an RGB image a 2d convolution would treat each color layer independently, mixing the information from each channel with further computations like pooling (averaging the resulting feature maps of each color or select the max value among the feature maps for each pixel, etc..) while 3d convolution would mix together the color layers already during convolution thanks to the 3d kernel which will sum together elements from different layers.
