I am solving a 3D heat transfer equation with variable boundaries (insulated, convective, radiative or free) using a F.D.M. technique. My geometry of choice is a cube. The purpose of my work is to get more familiar with such problems -- out of interest -- and their parallelizations using both distributed (MPI) and shared memory (OpenMP/Pthreads).
In order to make this question as neat as possible, I will avoid tediously long expressions -- as is usually the case in 3D -- and go straight to the point. My stencil is a 7-point kernel. I am solving on my cube by traversing the discretization of the inner points using an SOR approach (i.e. excluding boundary surfaces, vertices and corners). Hence I find myself confused about how one usually solves such equations via Finite Differences.
To set-up a (mass) matrix, instead of the naive point-by-point stencil approach is greatly perplexing. I say this because per each row of the mass matrix I end up having only 7 non-zero elements.
In other words, if my cube is discretized equally in x, y and z directions while having 100 points in each respectively, then my overall mass matrix dimensions would become: $(100x100x100)^2 = 1,000,000-by-1,000,000$ with only 7 nnz elements per row. This is too sparse and inefficient from the perspective of memory and certain algorithm applications.
This makes me wonder how inexperienced I currently am when trying to formulate a solution for such problems. That being said, how does one usually formulate such a huge sparse matrix to tackle such problems ? Surely my naive approach isn't the best nor the only solution out there (duh). Any explanation and/or examples would be very welcome.