Was trying to implement a poisson 2d solver using Conjugate Gradient Method, so from 10x10 grid the matrix becomes 100x100 (since we have 100 nodes to find the values at), 100x100 grid goes to 10000x10000 matrix, but matrix is sparse with 5 elements on each row (not pentadiagonal!).
Hence, receiving segmentation fault, perhaps due to bad memory allocation, but still, quite a big matrix to store and use, because zeros are still considered as doubles, which must potentially slow down the solver.
Thought of implementing data structure for sparse matrix to overcome the memory issues. Is there any other ideas on the memory management in solvers?
Watching these lectures And reading the book of the video lecturer "Numerical Methods for Partial Differential Equations: Finite Difference and Finite Volume Methods" by Sandip Mazumder.
Below is the draft C code to fill up the coefficient matrix for all nodes into dense matrix (not sparse).
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
int main(){
double L = 1;
int N = 4, M = 3;
int kmax = M*N;
double sum = 0;
int fi = 1;
double dx = L/(N-1);
double dy = L/(M-1);
double dx2 = dx*dx;
double dy2 = dy*dy;
double *S;
S = (double *) malloc(N*M * sizeof(double));
double *Q;
Q = (double *) malloc(N*M * sizeof(double));
int** Acoord = (int**)malloc(kmax * sizeof(int*));
for (int i = 0; i < kmax; i++)
Acoord[i] = (int*)malloc(5 * sizeof(int));
double** Aval = (double**)malloc(kmax * sizeof(double*));
for (int i = 0; i < kmax; i++)
Aval[i] = (double*)malloc(5 * sizeof(double));
// fill S(ource) rhs f function
for (int i = 0; i< N; i++){
for (int j = 0; j < M ; j++){
int k = j*N+i;
S[k] = -exp(-pow(i*dx-L/2,2)-pow(j*dy-L/2,2));
sum = sum + S[k];
}
}
sum = sum/N/M;
for (int i = 0; i < N; i++){
for (int j = 0; j < M; j++){
int k = j*N+i;
S[k] = S[k] - sum;
}
}
for (int i = 1; i < N-1; i++){
for (int j = 1; j < M-1; j++){
int k = j*N+i;
Aval[k][2] = -(2/dx2+2/dy2);
Aval[k][1] = 1/dx2;
Aval[k][3] = 1/dx2;
Aval[k][0] = 1/dy2;
Aval[k][4] = 1/dy2;
Acoord[k][0] = k - N;
Acoord[k][1] = k - 1;
Acoord[k][2] = k;
Acoord[k][3] = k + 1;
Acoord[k][4] = k + N;
Q[k] = S[k];
}
}
// fill boundary part
int i = N-1; // RIGHT
double JR = 0; // FLUX
for (int j = 1; j < M-1; j++){
int k = j*N+i;
Aval[k][0] = 1/dy2;
Aval[k][1] = 1/dx2;
Aval[k][2] = -1/dx2-2/dy2;
Aval[k][3] = 0;
Aval[k][4] = 1/dy2;
Acoord[k][0] = k - N;
Acoord[k][1] = k - 1;
Acoord[k][2] = k;
Acoord[k][3] = -1;
Acoord[k][4] = k + N;
Q[k] = S[k];
}
i = 0; // LEFT
double JL = 0; // FLUX
for (int j = 1; j < M-1; j++){
int k = j*N+i;
Aval[k][0] = 1/dy2;
Aval[k][1] = 0;
Aval[k][2] = -1/dx2-2/dy2;
Aval[k][3] = 1/dx2;
Aval[k][4] = 1/dy2;
Acoord[k][0] = k - N;
Acoord[k][1] = - 1;
Acoord[k][2] = k;
Acoord[k][3] = k + 1;
Acoord[k][4] = k + N;
Q[k] = S[k];
}
int j = 0; // BOT
double JB = 0; // FLUX
for (int i = 1; i < N-1; i++){
int k = j*N+i;
Aval[k][0] = 0;
Aval[k][1] = 1/dx2;
Aval[k][2] = -2/dx2-1/dy2;
Aval[k][3] = 1/dx2;
Aval[k][4] = 1/dy2;
Acoord[k][0] = -1;
Acoord[k][1] = k - 1;
Acoord[k][2] = k;
Acoord[k][3] = k + 1;
Acoord[k][4] = k + N;
Q[k] = S[k];
// printf("k = %d kmax = %d\n",k,kmax);
}
j = M-1; // TOP
double JT = 0; // FLUX
for (int i = 1; i < N-1; i++){
int k = j*N+i;
Aval[k][0] = 1/dy2;
Aval[k][1] = 1/dx2;
Aval[k][2] = -2/dx2-1/dy2;
Aval[k][3] = 1/dx2;
Aval[k][4] = 0;
Acoord[k][0] = k - N;
Acoord[k][1] = k - 1;
Acoord[k][2] = k;
Acoord[k][3] = k + 1;
Acoord[k][4] = -1;
Q[k] = S[k];
}
i = 0; // LEFT
j = 0; // BOT
int k = j*N+i;
Aval[k][0] = 0;
Aval[k][1] = 0;
Aval[k][2] = -1/dx2 - 1/dy2;
Aval[k][3] = 1/dx2;
Aval[k][4] = 1/dy2;
Acoord[k][0] = -1;
Acoord[k][1] = -1;
Acoord[k][2] = k;
Acoord[k][3] = k + 1;
Acoord[k][4] = k + N;
Q[k] = S[k];
i = 0; // LEFT
j = M-1; // TOP
k = j*N+i;
Aval[k][0] = 1/dy2;
Aval[k][1] = 0;
Aval[k][2] = -1/dx2 - 1/dy2;
Aval[k][3] = 1/dx2;
Aval[k][4] = 0;
Acoord[k][0] = k - N;
Acoord[k][1] = -1;
Acoord[k][2] = k;
Acoord[k][3] = k + 1;
Acoord[k][4] = -1;
Q[k] = S[k];
i = N-1; // RIGHT
j = M-1; // TOP
k = j*N+i;
Aval[k][0] = 1/dy2;
Aval[k][1] = 1/dx2;
Aval[k][2] = -1/dx2 - 1/dy2;
Aval[k][3] = 0;
Aval[k][4] = 0;
Acoord[k][0] = k - N;
Acoord[k][1] = k - 1;
Acoord[k][2] = k;
Acoord[k][3] = -1;
Acoord[k][4] = -1;
Q[k] = S[k];
i = N-1; // RIGHT
j = 0; // BOT
k = (j)*M+i;
Aval[k][0] = 0;
Aval[k][1] = 1/dx2;
Aval[k][2] = -1/dx2 - 1/dy2;
Aval[k][3] = 0;
Aval[k][4] = 1/dy2;
Acoord[k][0] = -1;
Acoord[k][1] = k - 1;
Acoord[k][2] = k;
Acoord[k][3] = -1;
Acoord[k][4] = k + N;
Q[k] = S[k];
// CJ solver
// CJ solver ends
FILE *f1 = fopen("testCJ_coord.txt", "w+");
if (f1 == NULL)
{
printf("Error opening file!\n");
exit(1);
}
for (int i = 0; i < kmax; i++){
fprintf(f1,"%d\t",i);
for (int j = 0; j < 5; j++){
fprintf(f1,"%d\t",Acoord[i][j]);
}
fprintf(f1,"\n");
}
fclose(f1);
FILE *f2 = fopen("testCJ_val.txt", "w+");
if (f2 == NULL)
{
printf("Error opening file!\n");
exit(1);
}
for (int i = 0; i < kmax; i++){
fprintf(f2,"%d\t",i);
for (int j = 0; j < 5; j++){
fprintf(f2,"%.4f\t",Aval[i][j]);
}
fprintf(f2,"\n");
}
fclose(f2);
FILE *f3 = fopen("testCJ_Q.txt", "w+");
if (f3 == NULL)
{
printf("Error opening file!\n");
exit(1);
}
for (int i = 0; i < kmax; i++){
fprintf(f3,"%d\t",i);
fprintf(f3,"%12.5f",Q[i]);
fprintf(f3,"\n");
}
fclose(f3);
free(S);
free(Q);
for (int i = 0; i < kmax; i++)
free(Aval[i]);
free(Aval);
for (int i = 0; i < kmax; i++)
free(Acoord[i]);
free(Acoord);
return 0;
}