ssantisi
/
simulacion-permeabilidad
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Debugging fourt

Browse Source
milestone_5
chortas 3 years ago
parent c64a61e071
commit 7ef3c2d2e5
4 changed files with 80 additions and 54 deletions
Show Stats Download Patch File Download Diff File

2
fftma_module/gen/include/chunk_array.h
Unescape Escape View File

@ -7,7 +7,7 @@
#include <string.h> #include <string.h>
#include <stdbool.h> #include <stdbool.h>
#define MAX_CHUNK_SIZE 512 #define MAX_CHUNK_SIZE 1500
typedef struct chunk_array { typedef struct chunk_array {
size_t init_pos; size_t init_pos;

8
fftma_module/gen/lib_src/fftma2.c
Unescape Escape View File

@ -51,10 +51,12 @@ void FFTMA2(struct vario_mod variogram, struct grid_mod grid, int n[3], struct r
NXYZ = grid.NX * grid.NY * grid.NZ; NXYZ = grid.NX * grid.NY * grid.NZ;
nxyz = NXYZ + 1; nxyz = NXYZ + 1;
printf("ntot es %d\n", ntot);
/*array initialization*/ /*array initialization*/
covar = chunk_array_create("covar.txt", ntot, 512); covar = chunk_array_create("covar.txt", ntot, 1500);
ireal = chunk_array_create("ireal.txt", ntot, 512); ireal = chunk_array_create("ireal.txt", ntot, 1500);
realization = chunk_array_create("realization.txt", ntot, 512); realization = chunk_array_create("realization.txt", ntot, 1500);
workr = (double*)malloc(nmax * sizeof(double)); workr = (double*)malloc(nmax * sizeof(double));
testmemory(workr); testmemory(workr);

119
fftma_module/gen/lib_src/fourt.c
Unescape Escape View File

@ -101,14 +101,9 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int
ntot *= nn[idim]; ntot *= nn[idim];
} }
//printf("104\n");
chunk_array_read(datar); chunk_array_read(datar);
//printf("107\n");
chunk_array_read(datai); chunk_array_read(datai);
//printf("110\n");
/*main loop for each dimension*/ /*main loop for each dimension*/
np1 = 1; np1 = 1;
for (idim = 1; idim <= ndim; idim++) { for (idim = 1; idim <= ndim; idim++) {
@ -194,12 +189,12 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int
ntot /= 2; ntot /= 2;
i = 1; i = 1;
for (j = 1; j <= ntot; j++) { for (j = 1; j <= ntot; j++) {
//printf("196\n");
chunk_array_get(datar, i, &valueri); chunk_array_get(datar, i, &valueri);
printf("datar[i] = %f\n", valueri);
chunk_array_get(datar, i+1, &valueri1); chunk_array_get(datar, i+1, &valueri1);
printf("datar[i1] = %f\n", valueri1);
chunk_array_save(datar, j, valueri); chunk_array_save(datar, j, valueri);
chunk_array_save(datai, j, valueri1); chunk_array_save(datai, j, valueri1);
//printf("201\n");
i += 2; i += 2;
} }
@ -217,22 +212,23 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int
for (i3 = i1; i3 <= ntot; i3 += np2) { for (i3 = i1; i3 <= ntot; i3 += np2) {
j3 = j + i3 - i2; j3 = j + i3 - i2;
//printf("219\n");
chunk_array_get(datar, i3, &valueri3); chunk_array_get(datar, i3, &valueri3);
chunk_array_get(datai, i3, &valueii3); chunk_array_get(datai, i3, &valueii3);
chunk_array_get(datar, j3, &valuerj3); chunk_array_get(datar, j3, &valuerj3);
chunk_array_get(datai, j3, &valueij3); chunk_array_get(datai, j3, &valueij3);
//printf("224\n");
printf("datar[i3] = %f\n", valueri3);
printf("datai[i3] = %f\n", valueii3);
printf("datar[j3] = %f\n", valuerj3);
printf("datai[j3] = %f\n", valueij3);
tempr = valueri3; tempr = valueri3;
tempi = valueii3; tempi = valueii3;
//printf("229\n");
chunk_array_save(datar, i3, valuerj3); chunk_array_save(datar, i3, valuerj3);
chunk_array_save(datai, i3, valueij3); chunk_array_save(datai, i3, valueij3);
chunk_array_save(datar, j3, tempr); chunk_array_save(datar, j3, tempr);
chunk_array_save(datai, j3, tempi); chunk_array_save(datai, j3, tempi);
//printf("234\n");
} }
} }
@ -257,10 +253,12 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int
for (i3 = i1; i3 <= ntot; i3 += np2) { for (i3 = i1; i3 <= ntot; i3 += np2) {
j = i3; j = i3;
for (i = 1; i <= n; i++) { for (i = 1; i <= n; i++) {
//printf("259\n");
chunk_array_get(datar, j, &valuerj); chunk_array_get(datar, j, &valuerj);
chunk_array_get(datai, j, &valueij); chunk_array_get(datai, j, &valueij);
//printf("262\n");
printf("datar[j] = %f\n", valuerj);
printf("datai[j] = %f\n", valueij);
if (icase != 3) { if (icase != 3) {
workr[i] = valuerj; workr[i] = valuerj;
worki[i] = valueij; worki[i] = valueij;
@ -284,10 +282,8 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int
i2max = i3 + np2 - np1; i2max = i3 + np2 - np1;
i = 1; i = 1;
for (i2 = i3; i2 <= i2max; i2 += np1) { for (i2 = i3; i2 <= i2max; i2 += np1) {
//printf("286\n");
chunk_array_save(datar, i2, workr[i]); chunk_array_save(datar, i2, workr[i]);
chunk_array_save(datai, i2, worki[i]); chunk_array_save(datai, i2, worki[i]);
//printf("289\n");
i++; i++;
} }
} }
@ -308,17 +304,20 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int
L310: L310:
j = i1; j = i1;
for (i = imin; i <= ntot; i += istep) { for (i = imin; i <= ntot; i += istep) {
//printf("310\n");
chunk_array_get(datar, i, &tempr); chunk_array_get(datar, i, &tempr);
chunk_array_get(datai, i, &tempi); chunk_array_get(datai, i, &tempi);
chunk_array_get(datar, j, &valuerj); chunk_array_get(datar, j, &valuerj);
chunk_array_get(datai, j, &valueij); chunk_array_get(datai, j, &valueij);
printf("datar[i] = %f\n", tempr);
printf("datai[i] = %f\n", tempi);
printf("datar[j] = %f\n", valuerj);
printf("datai[j] = %f\n", valueij);
chunk_array_save(datar, i, valuerj - tempr); chunk_array_save(datar, i, valuerj - tempr);
chunk_array_save(datai, i, valueij - tempi); chunk_array_save(datai, i, valueij - tempi);
chunk_array_save(datar, j, valuerj + tempr); chunk_array_save(datar, j, valuerj + tempr);
chunk_array_save(datai, j, valueij + tempi); chunk_array_save(datai, j, valueij + tempi);
//printf("320\n");
j += istep; j += istep;
} }
imin = 2 * imin - i1; imin = 2 * imin - i1;
@ -337,31 +336,37 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int
j = imin - istep / 2; j = imin - istep / 2;
for (i = imin; i <= ntot; i += istep) { for (i = imin; i <= ntot; i += istep) {
if (ifrwd != 0) { if (ifrwd != 0) {
//printf("339\n");
chunk_array_get(datai, i, &tempr); chunk_array_get(datai, i, &tempr);
chunk_array_get(datar, i, &tempi); chunk_array_get(datar, i, &tempi);
//printf("342\n");
printf("datai[i] = %f\n", tempr);
printf("datar[i] = %f\n", tempi);
tempi = -tempi; tempi = -tempi;
} else { } else {
//printf("345\n");
chunk_array_get(datai, i, &tempr); chunk_array_get(datai, i, &tempr);
printf("datai[i] = %f\n", tempr);
tempr = -tempr; tempr = -tempr;
chunk_array_get(datar, i, &tempi); chunk_array_get(datar, i, &tempi);
//printf("349\n"); printf("datar[i] = %f\n", tempi);
} }
//printf("351\n");
chunk_array_get(datar, j, &valuerj); chunk_array_get(datar, j, &valuerj);
chunk_array_get(datai, j, &valueij); chunk_array_get(datai, j, &valueij);
printf("datar[j] = %f\n", valuerj);
printf("datai[j] = %f\n", valueij);
chunk_array_save(datar, i, valuerj - tempr); chunk_array_save(datar, i, valuerj - tempr);
chunk_array_save(datai, i, valueij - tempi); chunk_array_save(datai, i, valueij - tempi);
chunk_array_get(datar, j, &valuerj); chunk_array_get(datar, j, &valuerj);
chunk_array_get(datai, j, &valueij); chunk_array_get(datai, j, &valueij);
printf("datar[j] = %f\n", valuerj);
printf("datai[j] = %f\n", valueij);
chunk_array_save(datar, j, valuerj + tempr); chunk_array_save(datar, j, valuerj + tempr);
chunk_array_save(datai, j, valueij + tempi); chunk_array_save(datai, j, valueij + tempi);
//printf("363\n");
j += istep; j += istep;
} }
@ -398,19 +403,22 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int
L510: L510:
j = imin - istep / 2; j = imin - istep / 2;
for (i = imin; i <= ntot; i += istep) { for (i = imin; i <= ntot; i += istep) {
//printf("400\n");
chunk_array_get(datar, i, &valueri); chunk_array_get(datar, i, &valueri);
chunk_array_get(datar, j, &valuerj); chunk_array_get(datar, j, &valuerj);
chunk_array_get(datai, i, &valueii); chunk_array_get(datai, i, &valueii);
chunk_array_get(datai, j, &valueij); chunk_array_get(datai, j, &valueij);
printf("datar[i] = %f\n", valueri);
printf("datar[j] = %f\n", valuerj);
printf("datai[i] = %f\n", valueii);
printf("datai[j] = %f\n", valueij);
tempr = valueri * wr - valueii * wi; tempr = valueri * wr - valueii * wi;
tempi = valueri * wi + valueii * wr; tempi = valueri * wi + valueii * wr;
chunk_array_save(datar, i, valuerj - tempr); chunk_array_save(datar, i, valuerj - tempr);
chunk_array_save(datai, i, valueij - tempi); chunk_array_save(datai, i, valueij - tempi);
chunk_array_save(datar, j, valuerj + tempr); chunk_array_save(datar, j, valuerj + tempr);
chunk_array_save(datai, j, valueij + tempi); chunk_array_save(datai, j, valueij + tempi);
//printf("412\n");
j += istep; j += istep;
} }
imin = 2 * imin - i1; imin = 2 * imin - i1;
@ -463,20 +471,24 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int
j3max = j2 + np2 - ifp2; j3max = j2 + np2 - ifp2;
for (j3 = j2; j3 <= j3max; j3 += ifp2) { for (j3 = j2; j3 <= j3max; j3 += ifp2) {
j = jmin + ifp2 - ifp1; j = jmin + ifp2 - ifp1;
//printf("465\n");
chunk_array_get(datar, j, &sr); chunk_array_get(datar, j, &sr);
chunk_array_get(datai, j, &si); chunk_array_get(datai, j, &si);
//printf("468\n");
printf("datar[j] = %f\n", sr);
printf("datai[j] = %f\n", si);
oldsr = 0.; oldsr = 0.;
oldsi = 0.; oldsi = 0.;
j -= ifp1; j -= ifp1;
L620: L620:
stmpr = sr; stmpr = sr;
stmpi = si; stmpi = si;
//printf("475\n");
chunk_array_get(datar, j, &valuerj); chunk_array_get(datar, j, &valuerj);
chunk_array_get(datai, j, &valueij); chunk_array_get(datai, j, &valueij);
//printf("478\n");
printf("datar[j] = %f\n", valuerj);
printf("datai[j] = %f\n", valueij);
sr = twowr * sr - oldsr + valuerj; sr = twowr * sr - oldsr + valuerj;
si = twowr * si - oldsi + valueij; si = twowr * si - oldsi + valueij;
@ -498,10 +510,8 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int
for (j2 = i3; j2 <= j2max; j2 += ifp1) { for (j2 = i3; j2 <= j2max; j2 += ifp1) {
j3max = j2 + np2 - ifp2; j3max = j2 + np2 - ifp2;
for (j3 = j2; j3 <= j3max; j3 += ifp2) { for (j3 = j2; j3 <= j3max; j3 += ifp2) {
//printf("500\n");
chunk_array_save(datar, j3, workr[i]); chunk_array_save(datar, j3, workr[i]);
chunk_array_save(datai, j3, worki[i]); chunk_array_save(datai, j3, worki[i]);
//printf("503\n");
i++; i++;
} }
} }
@ -545,12 +555,16 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int
L710: L710:
j = jmin; j = jmin;
for (i = imin; i <= ntot; i += np2) { for (i = imin; i <= ntot; i += np2) {
//printf("547\n");
chunk_array_get(datar, i, &valueri); chunk_array_get(datar, i, &valueri);
chunk_array_get(datai, i, &valueii); chunk_array_get(datai, i, &valueii);
chunk_array_get(datar, j, &valuerj); chunk_array_get(datar, j, &valuerj);
chunk_array_get(datai, j, &valueij); chunk_array_get(datai, j, &valueij);
printf("datar[i] = %f\n", valueri);
printf("datai[i] = %f\n", valueii);
printf("datar[j] = %f\n", valuerj);
printf("datai[j] = %f\n", valueij);
sumr = (valueri + valuerj) / 2.; sumr = (valueri + valuerj) / 2.;
sumi = (valueii + valueij) / 2.; sumi = (valueii + valueij) / 2.;
difr = (valueri - valuerj) / 2.; difr = (valueri - valuerj) / 2.;
@ -578,10 +592,9 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int
if (ifrwd == 0) if (ifrwd == 0)
goto L740; goto L740;
for (i = imin; i <= ntot; i += np2) { for (i = imin; i <= ntot; i += np2) {
//printf("580\n");
chunk_array_get(datai, i, &valueii); chunk_array_get(datai, i, &valueii);
printf("datai[i] = %f\n", valueii);
chunk_array_save(datai, i, -valueii); chunk_array_save(datai, i, -valueii);
//printf("583\n");
} }
L740: L740:
np2 *= 2; np2 *= 2;
@ -593,61 +606,67 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int
i = imin; i = imin;
goto L755; goto L755;
L750: ; L750: ;
//printf("595\n");
chunk_array_get(datar, i, &valueri); chunk_array_get(datar, i, &valueri);
chunk_array_get(datai, i, &valueii); chunk_array_get(datai, i, &valueii);
printf("datar[i] = %f\n", valueri);
printf("datai[i] = %f\n", valueii);
chunk_array_save(datar, j, valueri); chunk_array_save(datar, j, valueri);
chunk_array_save(datai, j, -valueii); chunk_array_save(datai, j, -valueii);
//printf("600\n");
L755: L755:
i++; i++;
j--; j--;
if (i < imax) if (i < imax)
goto L750; goto L750;
//printf("606\n");
chunk_array_get(datar, imin, &valuerimin); chunk_array_get(datar, imin, &valuerimin);
chunk_array_get(datai, imin, &valueiimin); chunk_array_get(datai, imin, &valueiimin);
printf("datar[imin] = %f\n", valuerimin);
printf("datai[imin] = %f\n", valueiimin);
chunk_array_save(datar, j, valuerimin - valueiimin); chunk_array_save(datar, j, valuerimin - valueiimin);
chunk_array_save(datai, j, 0.); chunk_array_save(datai, j, 0.);
//printf("612\n");
if (i >= j) { if (i >= j) {
goto L780; goto L780;
} else { } else {
goto L770; goto L770;
} }
L765: L765:
//printf("619\n");
chunk_array_get(datar, i, &valueri); chunk_array_get(datar, i, &valueri);
chunk_array_get(datai, i, &valueii); chunk_array_get(datai, i, &valueii);
printf("datar[i] = %f\n", valueri);
printf("datai[i] = %f\n", valueii);
chunk_array_save(datar, j, valueri); chunk_array_save(datar, j, valueri);
chunk_array_save(datai, j, valueii); chunk_array_save(datai, j, valueii);
//printf("625\n");
L770: L770:
i--; i--;
j--; j--;
if (i > imin) if (i > imin)
goto L765; goto L765;
//printf("632\n");
chunk_array_get(datar, imin, &valuerimin); chunk_array_get(datar, imin, &valuerimin);
chunk_array_get(datai, imin, &valueiimin); chunk_array_get(datai, imin, &valueiimin);
printf("datar[imin] = %f\n", valuerimin);
printf("datai[imin] = %f\n", valueiimin);
chunk_array_save(datar, j, valuerimin + valueiimin); chunk_array_save(datar, j, valuerimin + valueiimin);
chunk_array_save(datai, j, 0.); chunk_array_save(datai, j, 0.);
//printf("638\n");
imax = imin; imax = imin;
goto L745; goto L745;
L780: ; L780: ;
//printf("643\n");
chunk_array_get(datai, 1, &valuei1); chunk_array_get(datai, 1, &valuei1);
chunk_array_get(datar, 1, &valuer1); chunk_array_get(datar, 1, &valuer1);
printf("datai[1] = %f\n", valuei1);
printf("datar[1] = %f\n", valuer1);
chunk_array_save(datar, 1, valuei1 + valuer1); chunk_array_save(datar, 1, valuei1 + valuer1);
chunk_array_save(datai, 1, 0.); chunk_array_save(datai, 1, 0.);
//printf("649\n");
goto L900; goto L900;
/*complete a real transform for the 2nd, 3rd, ... dimension by conjugate symmetries*/ /*complete a real transform for the 2nd, 3rd, ... dimension by conjugate symmetries*/
@ -665,25 +684,27 @@ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int
if (idim > 2) { if (idim > 2) {
j = jmax + np0; j = jmax + np0;
for (i = imin; i <= imax; i++) { for (i = imin; i <= imax; i++) {
//printf("667\n");
chunk_array_get(datar, j, &valuerj); chunk_array_get(datar, j, &valuerj);
chunk_array_get(datai, j, &valueij); chunk_array_get(datai, j, &valueij);
printf("datar[j] = %f\n", valuerj);
printf("datai[j] = %f\n", valueij);
chunk_array_save(datar, i, valuerj); chunk_array_save(datar, i, valuerj);
chunk_array_save(datai, i, -valueij); chunk_array_save(datai, i, -valueij);
//printf("673\n");
j--; j--;
} }
} }
j = jmax; j = jmax;
for (i = imin; i <= imax; i += np0) { for (i = imin; i <= imax; i += np0) {
//printf("679\n");
chunk_array_get(datar, j, &valuerj); chunk_array_get(datar, j, &valuerj);
chunk_array_get(datai, j, &valueij); chunk_array_get(datai, j, &valueij);
printf("datar[j] = %f\n", valuerj);
printf("datai[j] = %f\n", valueij);
chunk_array_save(datar, i, valuerj); chunk_array_save(datar, i, valuerj);
chunk_array_save(datai, i, -valueij); chunk_array_save(datai, i, -valueij);
//printf("685\n");
j -= np0; j -= np0;
} }

5
fftma_module/gen/test.py
Unescape Escape View File

@ -29,6 +29,9 @@ variance=3.5682389
typ=3 typ=3
k=gen(nx, ny, nz, dx, dy, dz, seed, variograms, mean, variance, typ, 8) k=gen(nx, ny, nz, dx, dy, dz, seed, variograms, mean, variance, typ, 8)
np.save(f"out_{N}.npy",k)
#np.save(f"out_{N}.npy",k)
#k2 = np.load("out_2_16.npy")
del k del k
gc.collect() gc.collect()
Write Preview
Loading…
Cancel
Save