Memory leak Valgrind doesn't report but code using so much memory

I have posted recently about a segmentation fault issue that I have been getting in my code when running with Valgrind, I have got rid of the segmentation fault problem however, now the code is running, but is "eating" so much memory which indicates a memory leak somewhere that valgrind doesn't necessarily see? I am not sure really. For example it would run to the 40th iteration before the kernel kills my process.

previous post: Segmentation Fault with Valgrind

If I run the code with the corrections to the code: (initializing my arrays and changes size of arrays to be large enough)

char save [16] = {0};
snprintf(save, 16,  "%d", saveNum);
const char *type = " .txt";
            
char Nfilename[16] = {0};
strcat(Nfilename, "N");
strcat(Nfilename, save);

char KEfilename[16] = {0};
strcat(KEfilename, "KE");
strcat(KEfilename, save);

char KIfilename[16] = {0};
strcat(KIfilename, "KI");
strcat(KIfilename, save);

char Pfilename[16] = {0};
strcat(Pfilename, "P");
strcat(Pfilename, save);

print2DArrf(Nfilename, N); 
print2DArrf(KEfilename, KE);
print2DArrf(KIfilename, KI);
print2DArrf(Pfilename, P);

This has fixed the segmentation fault issue I was having at the exact iteration every single time. However, if I run the code for few iterations my computer would run out of RAM, I tried running the code with the activity memory monitor on Mac (since Valgrind takes days to run on debugging mode and give no memory leak issues) and the result is the code reaches to 16GB within few minutes!! I made sure to free everything, so something else must be happening.

I am wondering what's going on, and why is it taking so much memory in the code and how to fix it? I thought looking at each one of the major functions (shown below) in my code and run each one separately to pinpoint the issue but that seems tedious and might not really help. Is it actually a memory leak?

These are the major functions that I was looking to run individually to tackle this "memory" issue:

// calculate some values 
        convolution(Nk, KEk, Pek);
        Fourier3D(Pek, kb, Pek);

        
        convolution(Nk, KIk, Pik);
        Fourier3D(Pik, kb, Pik);
        
        Freqk(Nk, KIk,KEk , kb, eps0, mi, me, ri, rn, nn, Oci, Oce, e, nuink, nuiek, nuiik, nuenk, nueek, nueik, isigPk, NNk);
        
        dk(Nk, kx, dndxk);
        dk(Nk, ky, dndyk);

        calcPSk(dndxk, dndyk, Pik, Pek, nuink, nuiek, nuenk, nueik, isigPk, Oci, Oce, u, B, ksqu, potSourcek);

Example of the main code:

#include<math.h>
#include<stdio.h>
#include "Functions.h"
#include "fftw3.h"
#include <cstring>

int main(){         
    
    // Setting parameters 
    double frequency = 200.0;    
    double dt = 0.1;        
    double tend = 5000.; 
    double err_max = 1e-7;          

    int P_iter_max = 500;      

    int iter_max = tend / dt;
    int saveNumber = 1;

//Initialize Parameters 

double *N;
    N = (double*) fftw_malloc(nx*ny*sizeof(double));

fftw_complex *Nk;
    Nk = (fftw_complex*) fftw_malloc(nx*nyk*sizeof(fftw_complex)); 

double *KIk;
    KIk = (double*) fftw_malloc(nx*ny*sizeof(double));

double *KE;
    KE = (double*) fftw_malloc(nx*ny*sizeof(double));
    
fftw_complex *KEk;
    KEk = (fftw_complex*) fftw_malloc(nx*nyk* sizeof(fftw_complex)); 
    
    
double *P;
    P = (double*) fftw_malloc(nx*ny*sizeof(double));
    
    fftw_complex *Pk;
    Pk = (fftw_complex*) fftw_malloc(nx*nyk* sizeof(fftw_complex)); 

//Save I.Cs 

    int P_iter = Potk(invnk, dndxk, dndyk, Pk, PSk, kx, ky, ninvksq, err_max, P_iter_max);



    for (int iter = 0; iter < iter_max; iter++){

         // calculate some values 
        convolution(Nk, KEk, Pek);
        Fourier3D(Pek, kb, Pek);

        
        convolution(Nk, KIk, Pik);
        Fourier3D(Pik, kb, Pik);
        
        Freqk(Nk, KIk,KEk , kb, eps0, mi, me, ri, rn, nn, Oci, Oce, e, nuink, nuiek, nuiik, nuenk, nueek, nueik, isigPk, NNk);
        
        dk(Nk, kx, dndxk);
        dk(Nk, ky, dndyk);

        calcPSk(dndxk, dndyk, Pik, Pek, nuink, nuiek, nuenk, nueik, isigPk, Oci, Oce, u, B, ksqu, potSourcek);
    



         //Check for convergence

         // Save I.C loop

        P_iter = Pk(NNk, dndxk, dndyk, phik, potSourcek, kx, ky, ninvksqu, err_max, P_iter_max);
    } // seems to be the correct location of for loop end (see comment)
    double time[iter_max + 1];
    time[0] = 0;
    
    c2rfft(Nk, N);
    c2rfft(KEk, KE);    
    c2rfft(KIk, KI);
    c2rfft(Pk, P);
    
    char Ninitial[] = "N_initial.txt";
    char KEinitial[] = "KE_initial.txt";
    char KIinitial[] = "KI_initial.txt";
    char Pinitial[] = "P_initial.txt";
    

    print2DArrf(Ninitial, N);
    print2DArrf(KEinitial, KE);
    print2DArrf(KIinitial, KI);
    print2DArrf(Pinitial, P);
} // close of main() (??)

I Would really appreciate some insight on how to tackle this problem. I am new to this and I am finding some difficulty understanding the issue.

Solution

Make sure you free() all malloc'ed memory from heap

double *N;
    N = (double*) fftw_malloc(nx*ny*sizeof(double));

fftw_complex *Nk;
    Nk = (fftw_complex*) fftw_malloc(nx*nyk*sizeof(fftw_complex)); 

double *KIk;
    KIk = (double*) fftw_malloc(nx*ny*sizeof(double));

double *KE;
    KE = (double*) fftw_malloc(nx*ny*sizeof(double));
    
fftw_complex *KEk;
    KEk = (fftw_complex*) fftw_malloc(nx*nyk* sizeof(fftw_complex)); 
    
    
double *P;
    P = (double*) fftw_malloc(nx*ny*sizeof(double));
    
    fftw_complex *Pk;
    Pk = (fftw_complex*) fftw_malloc(nx*nyk* sizeof(fftw_complex)); 


    .
    .
    .
    free(N);
    free(Nk);
    free(KIk);
    free(KE);
    free(KEk);
    free(P);
    free(Pk);