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

Commit to debug

Browse Source
milestone_5
chortas 3 years ago
parent 8ab0002e44
commit 0d40eae129
10 changed files with 206 additions and 104 deletions
Show Stats Download Patch File Download Diff File

4
fftma_module/gen/include/geostat.h
Unescape Escape View File

@ -327,7 +327,7 @@ void coordinates(int maille, int i[3], struct grid_mod grid);
/*variogram: structure defined above */ /*variogram: structure defined above */
/*grid: structure defined above */ /*grid: structure defined above */
/*n: number of gridblocks along X,Y and Z*/ /*n: number of gridblocks along X,Y and Z*/
void covariance(double* covar, struct vario_mod variogram, struct grid_mod grid, int n[3], int cores); void covariance(chunk_array_t* covar, struct vario_mod variogram, struct grid_mod grid, int n[3], int cores);
/*computation of the covariance matrix for the well data*/ /*computation of the covariance matrix for the well data*/
/*well coordinates are given as a number of cells */ /*well coordinates are given as a number of cells */
@ -403,7 +403,7 @@ double exponential(double h);
/*workr: utility real part vector for storage */ /*workr: utility real part vector for storage */
/*worki: utility imaginary part vector for storage */ /*worki: utility imaginary part vector for storage */
/*The transformed data are returned in datar and datai*/ /*The transformed data are returned in datar and datai*/
void fourt(double* datar, double* datai, int nn[3], int ndim, int ifrwd, int icplx, double* workr, double* worki, int cores); void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int ifrwd, int icplx, double* workr, double* worki, int cores);
/*calculates F(x) = (1/a)*exp(-x*x/2)*/ /*calculates F(x) = (1/a)*exp(-x*x/2)*/
double funtrun1(double x); double funtrun1(double x);

6
fftma_module/gen/include/toolsFFTMA.h
Unescape Escape View File

@ -51,7 +51,7 @@ void FFTMA2(struct vario_mod variogram, struct grid_mod grid, int n[3], struct r
/* must be a Gaussian white noise */ /* must be a Gaussian white noise */
/*realization: structure defining a realization*/ /*realization: structure defining a realization*/
void prebuild_gwn(struct grid_mod grid, int n[3], struct realization_mod* realin, double* realization, int solver, int cores, long* seed); void prebuild_gwn(struct grid_mod grid, int n[3], struct realization_mod* realin, chunk_array_t* realization, int solver, int cores, long* seed);
/* build_real */ /* build_real */
/* build a realization in the spectral domain */ /* build a realization in the spectral domain */
@ -65,8 +65,8 @@ void prebuild_gwn(struct grid_mod grid, int n[3], struct realization_mod* realin
/*realization: vector defining the real part */ /*realization: vector defining the real part */
/*ireal: vector defining the i-part */ /*ireal: vector defining the i-part */
void build_real(int n[3], int NTOT, double* covar, double* realization, double* ireal, int cores); void build_real(int n[3], int NTOT, chunk_array_t* covar, chunk_array_t* realization, double* ireal, int cores);
void clean_real(struct realization_mod* realin, int n[3], struct grid_mod grid, double* vectorresult, struct realization_mod* realout); void clean_real(struct realization_mod* realin, int n[3], struct grid_mod grid, chunk_array_t* vectorresult, struct realization_mod* realout, int cores);
#endif // define _TOOLSFFTMA_H #endif // define _TOOLSFFTMA_H

2
fftma_module/gen/lib_src/Py_kgeneration.c
Unescape Escape View File

@ -27,6 +27,8 @@ void Py_kgeneration(long seed, struct grid_mod grid, struct statistic_mod stat,
n[1] = 0; n[1] = 0;
n[2] = 0; n[2] = 0;
printf("aaa");
generate(&seed, N, Z, cores); generate(&seed, N, Z, cores);
/*FFTMA*/ /*FFTMA*/

11
fftma_module/gen/lib_src/build_real.c
Unescape Escape View File

@ -1,4 +1,5 @@
#include "geostat.h" #include "geostat.h"
#include "chunk_array.h"
#include <math.h> #include <math.h>
#include <stdarg.h> #include <stdarg.h>
#include <stddef.h> #include <stddef.h>
@ -19,22 +20,26 @@
/*realization: vector defining the real part */ /*realization: vector defining the real part */
/*ireal: vector defining the i-part */ /*ireal: vector defining the i-part */
void build_real(int n[3], int NTOT, double* covar, double* realization, double* ireal, int cores) { void build_real(int n[3], int NTOT, chunk_array_t* covar, chunk_array_t* realization, double* ireal, int cores) {
int i, j, k, maille1; int i, j, k, maille1;
double temp; double temp;
chunk_array_read(realization);
/*decomposition and multiplication in the spectral domain*/ /*decomposition and multiplication in the spectral domain*/
for (k = 1; k <= n[2]; k++) { for (k = 1; k <= n[2]; k++) {
for (j = 1; j <= n[1]; j++) { for (j = 1; j <= n[1]; j++) {
for (i = 1; i <= n[0]; i++) { for (i = 1; i <= n[0]; i++) {
maille1 = i + (j - 1 + (k - 1) * n[1]) * n[0]; maille1 = i + (j - 1 + (k - 1) * n[1]) * n[0];
temp = covar[maille1]; chunk_aray_get(covar, maille1, &temp);
if (temp > 0.) { if (temp > 0.) {
temp = sqrt(temp) / (double)NTOT; temp = sqrt(temp) / (double)NTOT;
} else if (temp < 0.) { } else if (temp < 0.) {
temp = sqrt(-temp) / (double)NTOT; temp = sqrt(-temp) / (double)NTOT;
} }
realization[maille1] *= temp; double valuemaille1;
chunk_array_get(realization, maille1, &valuemaille1);
chunk_array_save(realization, maille1, temp * valuemaille1);
ireal[maille1] *= temp; ireal[maille1] *= temp;
} }
} }

9
fftma_module/gen/lib_src/clean_real.c
Unescape Escape View File

@ -1,4 +1,5 @@
#include "geostat.h" #include "geostat.h"
#include "chunk_array.h"
#include <math.h> #include <math.h>
#include <stdarg.h> #include <stdarg.h>
#include <stddef.h> #include <stddef.h>
@ -7,7 +8,7 @@
#include <string.h> #include <string.h>
#include <time.h> #include <time.h>
void clean_real(struct realization_mod* realin, int n[3], struct grid_mod grid, double* vectorresult, struct realization_mod* realout, int cores) { void clean_real(struct realization_mod* realin, int n[3], struct grid_mod grid, chunk_array_t* vectorresult, struct realization_mod* realout, int cores) {
int i, j, k, maille0, maille1; int i, j, k, maille0, maille1;
double NTOT; double NTOT;
@ -15,6 +16,8 @@ void clean_real(struct realization_mod* realin, int n[3], struct grid_mod grid,
/*is the output realization already allocated?*/ /*is the output realization already allocated?*/
/*if not, memory allocation*/ /*if not, memory allocation*/
chunk_array_read(vectorresult);
if (realout->vector == NULL || realout->n != realin->n) { if (realout->vector == NULL || realout->n != realin->n) {
realout->vector = (double*)malloc(realin->n * sizeof(double)); realout->vector = (double*)malloc(realin->n * sizeof(double));
if (realout->vector == NULL) { if (realout->vector == NULL) {
@ -31,7 +34,9 @@ void clean_real(struct realization_mod* realin, int n[3], struct grid_mod grid,
maille0 = i - 1 + (j - 1 + (k - 1) * grid.NY) * grid.NX; maille0 = i - 1 + (j - 1 + (k - 1) * grid.NY) * grid.NX;
/* Modif du 18 juin 2003 */ /* Modif du 18 juin 2003 */
/*realout->vector[maille0] = vectorresult[maille1]/(double) NTOT;*/ /*realout->vector[maille0] = vectorresult[maille1]/(double) NTOT;*/
realout->vector[maille0] = vectorresult[maille1]; double valuemaille1;
chunk_array_get(vectorresult, maille1, &valuemaille1);
realout->vector[maille0] = valuemaille1;
} }
} }
} }

29
fftma_module/gen/lib_src/covariance.c
Unescape Escape View File

@ -1,15 +1,18 @@
#include "geostat.h" #include "geostat.h"
#include "chunk_array.h"
#include <time.h> #include <time.h>
/*builds the sampled covariance function*/ /*builds the sampled covariance function*/
/*dimensions are even*/ /*dimensions are even*/
void covariance(double* covar, struct vario_mod variogram, struct grid_mod mesh, int n[3], int cores) { void covariance(chunk_array_t* covar, struct vario_mod variogram, struct grid_mod mesh, int n[3], int cores) {
int i, j, k, maille, n2[3], symmetric; int i, j, k, maille, n2[3], symmetric;
double di, dj, dk; double di, dj, dk;
for (i = 0; i < 3; i++) for (i = 0; i < 3; i++)
n2[i] = n[i] / 2; n2[i] = n[i] / 2;
chunk_array_read(covar);
for (i = 0; i <= n2[0]; i++) { for (i = 0; i <= n2[0]; i++) {
for (j = 0; j <= n2[1]; j++) { for (j = 0; j <= n2[1]; j++) {
for (k = 0; k <= n2[2]; k++) { for (k = 0; k <= n2[2]; k++) {
@ -19,12 +22,14 @@ void covariance(double* covar, struct vario_mod variogram, struct grid_mod mesh,
di = (double)i * mesh.DX; di = (double)i * mesh.DX;
dj = (double)j * mesh.DY; dj = (double)j * mesh.DY;
dk = (double)k * mesh.DZ; dk = (double)k * mesh.DZ;
covar[maille] = (double)cov_value(variogram, di, dj, dk, cores); chunk_array_save(covar, maille, (double)cov_value(variogram, di, dj, dk, cores));
if (k > 0 && k < n2[2] && j > 0 && j < n2[1] && i > 0 && i < n2[0]) { if (k > 0 && k < n2[2] && j > 0 && j < n2[1] && i > 0 && i < n2[0]) {
/*area 2*/ /*area 2*/
symmetric = 1 + n[0] - i + n[0] * (n[1] - j + n[1] * (n[2] - k)); symmetric = 1 + n[0] - i + n[0] * (n[1] - j + n[1] * (n[2] - k));
covar[symmetric] = covar[maille]; double value;
chunk_array_get(covar, maille, &value);
chunk_array_save(covar, symmetric, value);
} }
if (i > 0 && i < n2[0]) { if (i > 0 && i < n2[0]) {
@ -33,13 +38,15 @@ void covariance(double* covar, struct vario_mod variogram, struct grid_mod mesh,
dj = (double)j * mesh.DY; dj = (double)j * mesh.DY;
dk = (double)k * mesh.DZ; dk = (double)k * mesh.DZ;
maille = 1 + (n[0] - i) + n[0] * (j + n[1] * k); maille = 1 + (n[0] - i) + n[0] * (j + n[1] * k);
covar[maille] = (double)cov_value(variogram, di, dj, dk, cores); chunk_array_save(covar, maille, (double)cov_value(variogram, di, dj, dk, cores));
} }
if (k > 0 && k < n2[2] && j > 0 && j < n2[1]) { if (k > 0 && k < n2[2] && j > 0 && j < n2[1]) {
/*area 8*/ /*area 8*/
symmetric = 1 + i + n[0] * (n[1] - j + n[1] * (n[2] - k)); symmetric = 1 + i + n[0] * (n[1] - j + n[1] * (n[2] - k));
covar[symmetric] = covar[maille]; double value;
chunk_array_get(covar, maille, &value);
chunk_array_save(covar, symmetric, value);
} }
if (i > 0 && i < n2[0] && j > 0 && j < n2[1]) { if (i > 0 && i < n2[0] && j > 0 && j < n2[1]) {
@ -48,13 +55,15 @@ void covariance(double* covar, struct vario_mod variogram, struct grid_mod mesh,
dj = -(double)j * mesh.DY; dj = -(double)j * mesh.DY;
dk = (double)k * mesh.DZ; dk = (double)k * mesh.DZ;
maille = 1 + (n[0] - i) + n[0] * (n[1] - j + n[1] * k); maille = 1 + (n[0] - i) + n[0] * (n[1] - j + n[1] * k);
covar[maille] = (double)cov_value(variogram, di, dj, dk, cores); chunk_array_save(covar, maille, (double)cov_value(variogram, di, dj, dk, cores));
} }
if (k > 0 && k < n2[2]) { if (k > 0 && k < n2[2]) {
/*area 6*/ /*area 6*/
symmetric = 1 + i + n[0] * (j + n[1] * (n[2] - k)); symmetric = 1 + i + n[0] * (j + n[1] * (n[2] - k));
covar[symmetric] = covar[maille]; double value;
chunk_array_get(covar, maille, &value);
chunk_array_save(covar, symmetric, value);
} }
if (j > 0 && j < n2[1]) { if (j > 0 && j < n2[1]) {
@ -63,13 +72,15 @@ void covariance(double* covar, struct vario_mod variogram, struct grid_mod mesh,
dj = -(double)j * mesh.DY; dj = -(double)j * mesh.DY;
dk = (double)k * mesh.DZ; dk = (double)k * mesh.DZ;
maille = 1 + i + n[0] * (n[1] - j + n[1] * k); maille = 1 + i + n[0] * (n[1] - j + n[1] * k);
covar[maille] = (double)cov_value(variogram, di, dj, dk, cores); chunk_array_save(covar, maille, (double)cov_value(variogram, di, dj, dk, cores));
} }
if (k > 0 && k < n2[2] && i > 0 && i < n2[0]) { if (k > 0 && k < n2[2] && i > 0 && i < n2[0]) {
/*area 7*/ /*area 7*/
symmetric = 1 + n[0] - i + n[0] * (j + n[1] * (n[2] - k)); symmetric = 1 + n[0] - i + n[0] * (j + n[1] * (n[2] - k));
covar[symmetric] = covar[maille]; double value;
chunk_array_get(covar, maille, &value);
chunk_array_save(covar, symmetric, value);
} }
} }
} }

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

@ -1,4 +1,5 @@
#include "geostat.h" #include "geostat.h"
#include "chunk_array.h"
#include <math.h> #include <math.h>
#include <stdlib.h> #include <stdlib.h>
#include <string.h> #include <string.h>
@ -26,7 +27,10 @@ void FFTMA2(struct vario_mod variogram, struct grid_mod grid, int n[3], struct r
int NTOT, i, j, k, NMAX, NDIM, ntot, nmax, NXYZ, nxyz; int NTOT, i, j, k, NMAX, NDIM, ntot, nmax, NXYZ, nxyz;
int solver; int solver;
double temp; double temp;
double *ireal, *covar, *workr, *worki, *realization; chunk_array_t *covar, *ireal, *realization;
double *workr, *worki;
printf("Ayuda");
/*covariance axis normalization*/ /*covariance axis normalization*/
axes(variogram.ap, variogram.scf, variogram.Nvario); axes(variogram.ap, variogram.scf, variogram.Nvario);
@ -50,14 +54,9 @@ void FFTMA2(struct vario_mod variogram, struct grid_mod grid, int n[3], struct r
nxyz = NXYZ + 1; nxyz = NXYZ + 1;
/*array initialization*/ /*array initialization*/
covar = (double*)malloc(ntot * sizeof(double)); covar = chunk_array_create("covar.txt", ntot, 512);
testmemory(covar); ireal = chunk_array_create("ireal.txt", ntot, 512);
realization = chunk_array_create("realization.txt", ntot, 512);
ireal = (double*)malloc(ntot * sizeof(double));
testmemory(ireal);
realization = (double*)malloc(ntot * sizeof(double));
testmemory(realization);
workr = (double*)malloc(nmax * sizeof(double)); workr = (double*)malloc(nmax * sizeof(double));
testmemory(workr); testmemory(workr);
@ -83,12 +82,12 @@ void FFTMA2(struct vario_mod variogram, struct grid_mod grid, int n[3], struct r
/* build realization in spectral domain */ /* build realization in spectral domain */
build_real(n, NTOT, covar, realization, ireal, cores); build_real(n, NTOT, covar, realization, ireal, cores);
free(covar); chunk_array_free(covar);
/*backward fourier transform*/ /*backward fourier transform*/
fourt(realization, ireal, n, NDIM, 0, 1, workr, worki, cores); fourt(realization, ireal, n, NDIM, 0, 1, workr, worki, cores);
free(ireal); chunk_array_free(ireal);
free(workr); free(workr);
free(worki); free(worki);

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

@ -1,6 +1,7 @@
#include <math.h> #include <math.h>
#include <stdio.h> #include <stdio.h>
#include <time.h> #include <time.h>
#include "chunk_array.h"
/*fast fourier transform */ /*fast fourier transform */
/* THE COOLEY-TUKEY FAST FOURIER TRANSFORM */ /* THE COOLEY-TUKEY FAST FOURIER TRANSFORM */
@ -90,7 +91,7 @@
/* PROGRAM MODIFIED FROM A SUBROUTINE OF BRENNER */ /* PROGRAM MODIFIED FROM A SUBROUTINE OF BRENNER */
/* 10-06-2000, MLR */ /* 10-06-2000, MLR */
void fourt(double* datar, double* datai, int nn[3], int ndim, int ifrwd, int icplx, double* workr, double* worki, int cores) { void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int ifrwd, int icplx, double* workr, double* worki, int cores) {
int ifact[21], ntot, idim, np1, n, np2, m, ntwo, iff, idiv, iquot, irem, inon2, non2p, np0, nprev, icase, ifmin, i, j, jmax, np2hf, i2, i1max, i3, j3, i1, ifp1, ifp2, i2max, i1rng, istep, imin, imax, mmax, mmin, mstep, j1, j2max, j2, jmin, j3max, nhalf; int ifact[21], ntot, idim, np1, n, np2, m, ntwo, iff, idiv, iquot, irem, inon2, non2p, np0, nprev, icase, ifmin, i, j, jmax, np2hf, i2, i1max, i3, j3, i1, ifp1, ifp2, i2max, i1rng, istep, imin, imax, mmax, mmin, mstep, j1, j2max, j2, jmin, j3max, nhalf;
double theta, wstpr, wstpi, wminr, wmini, wr, wi, wtemp, thetm, wmstr, wmsti, twowr, sr, si, oldsr, oldsi, stmpr, stmpi, tempr, tempi, difi, difr, sumr, sumi, TWOPI = 6.283185307179586476925286766559; double theta, wstpr, wstpi, wminr, wmini, wr, wi, wtemp, thetm, wmstr, wmsti, twowr, sr, si, oldsr, oldsi, stmpr, stmpi, tempr, tempi, difi, difr, sumr, sumi, TWOPI = 6.283185307179586476925286766559;
@ -99,6 +100,9 @@ void fourt(double* datar, double* datai, int nn[3], int ndim, int ifrwd, int icp
ntot *= nn[idim]; ntot *= nn[idim];
} }
chunk_array_read(datar);
chunk_array_read(datai);
/*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++) {
@ -184,8 +188,11 @@ void fourt(double* datar, double* datai, int nn[3], int ndim, int ifrwd, int icp
ntot /= 2; ntot /= 2;
i = 1; i = 1;
for (j = 1; j <= ntot; j++) { for (j = 1; j <= ntot; j++) {
datar[j] = datar[i]; double valuei, valuei1;
datai[j] = datar[i + 1]; chunk_array_get(datar, i, &valuei);
chunk_array_get(datar, i+1, &valuei1);
chunk_array_save(datar, j, valuei);
chunk_array_save(datai, j, valuei1);
i += 2; i += 2;
} }
@ -202,12 +209,20 @@ void fourt(double* datar, double* datai, int nn[3], int ndim, int ifrwd, int icp
for (i1 = i2; i1 <= i1max; i1++) { for (i1 = i2; i1 <= i1max; i1++) {
for (i3 = i1; i3 <= ntot; i3 += np2) { for (i3 = i1; i3 <= ntot; i3 += np2) {
j3 = j + i3 - i2; j3 = j + i3 - i2;
tempr = datar[i3]; double valueri3, valueii3, valuerj3, valueij3;
tempi = datai[i3];
datar[i3] = datar[j3]; chunk_array_get(datar, i3, &valueri3);
datai[i3] = datai[j3]; chunk_array_get(datai, i3, &valueii3);
datar[j3] = tempr; chunk_array_get(datar, j3, &valuerj3);
datai[j3] = tempi; chunk_array_get(datai, j3, &valueij3);
tempr = valueri3;
tempi = valueii3;
chunk_array_save(datar, i3, valuerj3);
chunk_array_save(datai, i3, valueij3);
chunk_array_save(datar, j3, tempr);
chunk_array_save(datai, j3, tempi);
} }
} }
@ -232,11 +247,14 @@ void fourt(double* datar, double* datai, int nn[3], int ndim, int ifrwd, int icp
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++) {
double valuerj, valueij;
chunk_array_get(datar, j, &valuerj);
chunk_array_get(datai, j, &valueij);
if (icase != 3) { if (icase != 3) {
workr[i] = datar[j]; workr[i] = valuerj;
worki[i] = datai[j]; worki[i] = valueij;
} else { } else {
workr[i] = datar[j]; workr[i] = valuerj;
worki[i] = 0.; worki[i] = 0.;
} }
ifp2 = np2; ifp2 = np2;
@ -255,8 +273,8 @@ void fourt(double* datar, double* datai, int nn[3], int ndim, int ifrwd, int icp
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) {
datar[i2] = workr[i]; chunk_save(datar, i2, workr[i]);
datai[i2] = worki[i]; chunk_save(datai, i2, worki[i]);
i++; i++;
} }
} }
@ -277,12 +295,16 @@ void fourt(double* datar, double* datai, int nn[3], int ndim, int ifrwd, int icp
L310: L310:
j = i1; j = i1;
for (i = imin; i <= ntot; i += istep) { for (i = imin; i <= ntot; i += istep) {
tempr = datar[i]; chunk_array_get(datar, i, &tempr);
tempi = datai[i]; chunk_array_get(datai, i, &tempi);
datar[i] = datar[j] - tempr;
datai[i] = datai[j] - tempi; double valuerj, valueij;
datar[j] = datar[j] + tempr; chunk_array_get(datar, j, &valuerj);
datai[j] = datai[j] + tempi; chunk_array_get(datai, j, &valueij);
chunk_array_save(datar, i, valuerj - tempr);
chunk_array_save(datai, i, valueij - tempi);
chunk_array_save(datar, j, valuerj + tempr);
chunk_array_save(datai, j, valueij + tempi);
j += istep; j += istep;
} }
imin = 2 * imin - i1; imin = 2 * imin - i1;
@ -301,16 +323,26 @@ void fourt(double* datar, double* datai, int nn[3], int ndim, int ifrwd, int icp
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) {
tempr = datai[i]; chunk_array_get(datai, i, &tempr);
tempi = -datar[i]; chunk_array_get(datar, i, &tempi);
tempi = -tempi;
} else { } else {
tempr = -datai[i]; chunk_array_get(datai, i, &tempr);
tempi = datar[i]; chunk_array_get(datar, i, &tempi);
} }
datar[i] = datar[j] - tempr; double valuerj, valueij;
datai[i] = datai[j] - tempi; chunk_array_get(datar, j, &valuerj);
datar[j] += tempr; chunk_array_get(datai, j, &valueij);
datai[j] += tempi;
chunk_array_save(datar, i, valuerj - tempr);
chunk_array_save(datai, i, valueij - tempi);
chunk_array_get(datar, j, &valuerj);
chunk_array_get(datai, j, &valueij);
chunk_array_save(datar, j, valuerj + tempr);
chunk_array_save(datai, j, valueij + tempi);
j += istep; j += istep;
} }
@ -347,12 +379,18 @@ void fourt(double* datar, double* datai, int nn[3], int ndim, int ifrwd, int icp
L510: L510:
j = imin - istep / 2; j = imin - istep / 2;
for (i = imin; i <= ntot; i += istep) { for (i = imin; i <= ntot; i += istep) {
tempr = datar[i] * wr - datai[i] * wi; double valueri, valueii, valuerj, valueij;
tempi = datar[i] * wi + datai[i] * wr; chunk_array_get(datar, i, &valueri);
datar[i] = datar[j] - tempr; chunk_array_get(datar, j, &valuerj);
datai[i] = datai[j] - tempi; chunk_array_get(datai, i, &valueii);
datar[j] += tempr; chunk_array_get(datai, j, &valueij);
datai[j] += tempi;
tempr = valueri * wr - valueii * wi;
tempi = valueri * wi + valueii * wr;
chunk_data_save(datar, i, valuerj - tempr);
chunk_data_save(datai, i, valueij - tempi);
chunk_data_save(datar, j, valuerj + tempr);
chunk_data_save(datai, j, valueij + tempi);
j += istep; j += istep;
} }
imin = 2 * imin - i1; imin = 2 * imin - i1;
@ -405,23 +443,27 @@ void fourt(double* datar, double* datai, int nn[3], int ndim, int ifrwd, int icp
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;
sr = datar[j]; chunk_array_get(datar, j, &sr);
si = datai[j]; chunk_array_get(datai, j, &si);
oldsr = 0.; oldsr = 0.;
oldsi = 0.; oldsi = 0.;
j -= ifp1; j -= ifp1;
L620: L620:
stmpr = sr; stmpr = sr;
stmpi = si; stmpi = si;
sr = twowr * sr - oldsr + datar[j]; double valuerj, valueij;
si = twowr * si - oldsi + datai[j]; chunk_array_get(datar, j, &valuerj);
chunk_array_get(datai, j, &valueij);
sr = twowr * sr - oldsr + valuerj;
si = twowr * si - oldsi + valueij;
oldsr = stmpr; oldsr = stmpr;
oldsi = stmpi; oldsi = stmpi;
j -= ifp1; j -= ifp1;
if (j > jmin) if (j > jmin)
goto L620; goto L620;
workr[i] = wr * sr - wi * si - oldsr + datar[j]; workr[i] = wr * sr - wi * si - oldsr + valuerj;
worki[i] = wi * sr + wr * si - oldsi + datai[j]; worki[i] = wi * sr + wr * si - oldsi + valueij;
jmin += ifp2; jmin += ifp2;
i++; i++;
} }
@ -433,8 +475,8 @@ void fourt(double* datar, double* datai, int nn[3], int ndim, int ifrwd, int icp
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) {
datar[j3] = workr[i]; chunk_array_save(datar, j3, workr[i]);
datai[j3] = worki[i]; chunk_array_save(datai, j3, worki[i]);
i++; i++;
} }
} }
@ -478,16 +520,22 @@ void fourt(double* datar, double* datai, int nn[3], int ndim, int ifrwd, int icp
L710: L710:
j = jmin; j = jmin;
for (i = imin; i <= ntot; i += np2) { for (i = imin; i <= ntot; i += np2) {
sumr = (datar[i] + datar[j]) / 2.; double valueri, valueii, valuerj, valueij;
sumi = (datai[i] + datai[j]) / 2.; chunk_array_get(datar, i, &valueri);
difr = (datar[i] - datar[j]) / 2.; chunk_array_get(datai, i, &valueii);
difi = (datai[i] - datai[j]) / 2.; chunk_array_get(datar, j, &valuerj);
chunk_array_get(datai, j, &valueij);
sumr = (valueri + valuerj) / 2.;
sumi = (valueii + valueij) / 2.;
difr = (valueri - valuerj) / 2.;
difi = (valueii - valueij) / 2.;
tempr = wr * sumi + wi * difr; tempr = wr * sumi + wi * difr;
tempi = wi * sumi - wr * difr; tempi = wi * sumi - wr * difr;
datar[i] = sumr + tempr; chunk_data_save(datar, i, sumr + tempr);
datai[i] = difi + tempi; chunk_data_save(datai, i, difi + tempi);
datar[j] = sumr - tempr; chunk_data_save(datar, j, sumr - tempr);
datai[j] = tempi - difi; chunk_data_save(datai, j, tempi - difi);
j += np2; j += np2;
} }
imin++; imin++;
@ -504,7 +552,9 @@ void fourt(double* datar, double* datai, int nn[3], int ndim, int ifrwd, int icp
if (ifrwd == 0) if (ifrwd == 0)
goto L740; goto L740;
for (i = imin; i <= ntot; i += np2) { for (i = imin; i <= ntot; i += np2) {
datai[i] = -datai[i]; double valueii;
chunk_data_get(datai, i, &valueii);
chunk_data_save(datai, i, -valueii);
} }
L740: L740:
np2 *= 2; np2 *= 2;
@ -515,36 +565,56 @@ void fourt(double* datar, double* datai, int nn[3], int ndim, int ifrwd, int icp
imin = imax - nhalf; imin = imax - nhalf;
i = imin; i = imin;
goto L755; goto L755;
L750: L750: ;
datar[j] = datar[i]; double valueri, valueii;
datai[j] = -datai[i]; chunk_data_get(datar, i, &valueri);
chunk_data_get(datai, i, &valueii);
chunk_data_save(datar, j, valueri);
chunk_data_save(datai, j, -valueii);
L755: L755:
i++; i++;
j--; j--;
if (i < imax) if (i < imax)
goto L750; goto L750;
datar[j] = datar[imin] - datai[imin]; double valuerimin, valueiimin;
datai[j] = 0.; chunk_data_get(datar, imin, &valuerimin);
chunk_data_get(datai, imin, &valueiimin);
chunk_data_save(datar, j, valuerimin - valueiimin);
chunk_data_save(datai, j, 0.);
if (i >= j) { if (i >= j) {
goto L780; goto L780;
} else { } else {
goto L770; goto L770;
} }
L765: L765:
datar[j] = datar[i]; chunk_data_get(datar, i, &valueri);
datai[j] = datai[i]; chunk_data_get(datai, i, &valueii);
chunk_data_save(datar, j, valueri);
chunk_data_save(datai, j, valueii);
L770: L770:
i--; i--;
j--; j--;
if (i > imin) if (i > imin)
goto L765; goto L765;
datar[j] = datar[imin] + datai[imin];
datai[j] = 0.; chunk_data_get(datar, imin, &valuerimin);
chunk_data_get(datai, imin, &valueiimin);
chunk_data_save(datar, j, valuerimin + valueiimin);
chunk_data_save(datai, j, 0.);
imax = imin; imax = imin;
goto L745; goto L745;
L780: L780: ;
datar[1] += datai[1]; double valuei1, valuer1;
datai[1] = 0.; chunk_data_get(datai, 1, &valuei1);
chunk_data_get(datar, 1, &valuer1);
chunk_data_save(datar, 1, valuei1 + valuer1);
chunk_data_save(datai, 1, 0.);
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*/
@ -562,15 +632,24 @@ void fourt(double* datar, double* datai, int nn[3], int ndim, int ifrwd, int icp
if (idim > 2) { if (idim > 2) {
j = jmax + np0; j = jmax + np0;
for (i = imin; i <= imax; i++) { for (i = imin; i <= imax; i++) {
datar[i] = datar[j]; double valuerj, valueij;
datai[i] = -datai[j]; chunk_data_get(datar, j, &valuerj);
chunk_data_get(datai, j, &valueij);
chunk_data_save(datar, i, valuerj);
chunk_data_save(datai, i, -valueij);
j--; j--;
} }
} }
j = jmax; j = jmax;
for (i = imin; i <= imax; i += np0) { for (i = imin; i <= imax; i += np0) {
datar[i] = datar[j]; double valuerj, valueij;
datai[i] = -datai[j]; chunk_data_get(datar, j, &valuerj);
chunk_data_get(datai, j, &valueij);
chunk_data_save(datar, i, valuerj);
chunk_data_save(datai, i, -valueij);
j -= np0; j -= np0;
} }
} }

4
fftma_module/gen/lib_src/geostat.h
Unescape Escape View File

@ -328,7 +328,7 @@ void coordinates(int maille, int i[3], struct grid_mod grid);
/*variogram: structure defined above */ /*variogram: structure defined above */
/*grid: structure defined above */ /*grid: structure defined above */
/*n: number of gridblocks along X,Y and Z*/ /*n: number of gridblocks along X,Y and Z*/
void covariance(double* covar, struct vario_mod variogram, struct grid_mod grid, int n[3], int cores); void covariance(chunk_array_t* covar, struct vario_mod variogram, struct grid_mod grid, int n[3], int cores);
/*computation of the covariance matrix for the well data*/ /*computation of the covariance matrix for the well data*/
/*well coordinates are given as a number of cells */ /*well coordinates are given as a number of cells */
@ -404,7 +404,7 @@ double exponential(double h);
/*workr: utility real part vector for storage */ /*workr: utility real part vector for storage */
/*worki: utility imaginary part vector for storage */ /*worki: utility imaginary part vector for storage */
/*The transformed data are returned in datar and datai*/ /*The transformed data are returned in datar and datai*/
void fourt(double* datar, double* datai, int nn[3], int ndim, int ifrwd, int icplx, double* workr, double* worki, int cores); void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int ifrwd, int icplx, double* workr, double* worki, int cores);
/*calculates F(x) = (1/a)*exp(-x*x/2)*/ /*calculates F(x) = (1/a)*exp(-x*x/2)*/
double funtrun1(double x); double funtrun1(double x);

11
fftma_module/gen/lib_src/prebuild_gwn.c
Unescape Escape View File

@ -1,4 +1,5 @@
#include "geostat.h" #include "geostat.h"
#include "chunk_array.h"
#include <math.h> #include <math.h>
#include <stdarg.h> #include <stdarg.h>
#include <stddef.h> #include <stddef.h>
@ -20,7 +21,7 @@
/* must be a Gaussian white noise */ /* must be a Gaussian white noise */
/*realization: structure defining a realization*/ /*realization: structure defining a realization*/
void prebuild_gwn(struct grid_mod grid, int n[3], struct realization_mod* realin, double* realization, int solver, int cores, long* seed) { void prebuild_gwn(struct grid_mod grid, int n[3], struct realization_mod* realin, chunk_array_t* realization, int solver, int cores, long* seed) {
int i, j, k, maille0, maille1; int i, j, k, maille0, maille1;
int ntot; int ntot;
@ -32,14 +33,14 @@ void prebuild_gwn(struct grid_mod grid, int n[3], struct realization_mod* realin
*seed = -(*seed); *seed = -(*seed);
ntot = n[0] * n[1] * n[2]; ntot = n[0] * n[1] * n[2];
realization[0] = 0.; chunk_array_save(realization, 0, 0.);
/*printf("Antes de llamar a chunkarray read\n"); /*printf("Antes de llamar a chunkarray read\n");
chunk_array_read((*realin).vector_2); chunk_array_read((*realin).vector_2);
printf("Despues de llamar a chunkarray read\n");*/ printf("Despues de llamar a chunkarray read\n");*/
if (solver == 1) { if (solver == 1) {
for (i = 0; i < ntot; i++) { for (i = 0; i < ntot; i++) {
double value = gasdev(seed, &idum2, &iy, iv, cores); double value = gasdev(seed, &idum2, &iy, iv, cores);
realization[i+1] = value; chunk_array_save(realization, i+1, value);
//chunk_array_get((*realin).vector_2, i, &realization[i + 1]); //chunk_array_get((*realin).vector_2, i, &realization[i + 1]);
} }
} else { } else {
@ -51,10 +52,10 @@ void prebuild_gwn(struct grid_mod grid, int n[3], struct realization_mod* realin
//maille0 = i - 1 + (j - 1 + (k - 1) * grid.NY) * grid.NX; //maille0 = i - 1 + (j - 1 + (k - 1) * grid.NY) * grid.NX;
//printf("Maille0 es %d", maille0); //printf("Maille0 es %d", maille0);
double value = gasdev(seed, &idum2, &iy, iv, cores); double value = gasdev(seed, &idum2, &iy, iv, cores);
realization[maille1] = value; chunk_array_save(realization, maille1, value);
//chunk_array_get((*realin).vector_2, maille0, &realization[maille1]); //chunk_array_get((*realin).vector_2, maille0, &realization[maille1]);
} else { } else {
realization[maille1] = 0.; chunk_array_save(realization, maille1, 0.);
} }
} }
} }

Write Preview
Loading…
Cancel
Save