$ valgrind trackorigins=yes ./program.bin
==50190== Conditional jump or move depends on uninitialised value(s)
==50190== at 0x4CB87BC: sqrt (in /usr/lib/libm2.32.so)
==50190== by 0x10A3A8: absComp(double (*) [2], double*) (arithmeticFunctions.cpp:254)
==50190== by 0x10A536: max_absComp(double (*) [2]) (arithmeticFunctions.cpp:292)
==50190== by 0x10BEFD: potentialk(double (*) [2], double (*) [2], double (*) [2], double (*) [2], double (*) [2], double*, double*, double*, double, int) (spectralFunctions.cpp:491)
==50190== by 0x10D93E: main (testBed.cpp:304)
==50190== Uninitialised value was created by a heap allocation
==50190== at 0x483D118: memalign (vg_replace_malloc.c:907)
==50190== by 0x10D8DD: main (testBed.cpp:300) 
1 2 3 4 5 6 7 8 9 10 11 12

//potentialk()
// Calculate maximum of absolute value of previous phi
phik_max_old = max_absComp(phik); //485
// Multiply by ninvksqu to get the updated phi
Arr3DArr2DMult(RHS, ninvksqu, phik);
// Calculate maximum of absolute value of updated phi
phik_max = max_absComp(phik); //491
//main()
potSourcek = (fftw_complex*) fftw_malloc(nx*nyk* sizeof(fftw_complex)); //testBed.cpp:300
 
notice:
 it complains in line 491, but not on 485, however both are the same call:
max_absComp(phik)
 the uninitialised value was pointed out to `potSourcek', not `phik'
so, phik was fine at line 485, you called `Arr3DArr2DMult(RHS, ninvksqu, phik)' and now it's uninitialised
checking out `RHS', it's filled like
RHS[j + nyk*i][k] = potSourcek[j + nyk*i][k]  gradNgradPhi_x[j + nyk*i][k]  gradNgradPhi_y[j + nyk*i][k]; //476
using the invalid values from `potSourcek'
so there you are, initialise `potSourcek'
I was hoping that valgrind would detect faulty line 476, however
apart:
 `Comp' makes me think on «compare/comparison», not «computing», I'll suggest to drop the Comp suffix on absComp(), max_absComp()
 «Arr3DArr2DMult() [This function multiplies a 3D Fourier array by a 2D real valued array.]» no, it multiplies a "complex" array by a "real" array. They are both 2d (xy plane) but are stored in a single dimension array, so the function may work for 1d signals, just pass a size argument

int potentialk(double invnk[16][ncomp], double dndxk[16][ncomp], double dndyk[16][ncomp] ...)
first dimension is meaningless to the compiler, it may be used for documentation but that's not your case, as your arrays are not size 16 (or limited to size 16)
 test first with toy examples, things that you may see at a glance that are wrong, a 512x257 matrix is too big
 a little explanation of what the program does would be nice (¿poisson heat equation? ¿what conditions?)
 you are not charged by vowels, don't compress your variable names (you are not the only one that sees your code, not all are immersed in your study subject)
 read and fix the warnings that your compiler throws (
W{all,extra,pedantic} to compile with warnings)
regards