intermediate commit
parent
de5f4c8fe6
commit
028d0a8d9f
@ -1,110 +0,0 @@
|
|||||||
#include <stdlib.h>
|
|
||||||
#include <string.h>
|
|
||||||
#include <math.h>
|
|
||||||
#include "geostat.h"
|
|
||||||
|
|
||||||
|
|
||||||
/*FAST FOURIER TRANSFORM MOVING AVERAGE METHOD */
|
|
||||||
/*Turns a Gaussian white noise vector into a */
|
|
||||||
/*spatially correlated vector */
|
|
||||||
/*input: */
|
|
||||||
/*variogram: structure defining the variogram */
|
|
||||||
/* model */
|
|
||||||
/*grid: structure defining the grid */
|
|
||||||
/*n: vector with the number of cells along the */
|
|
||||||
/* X, Y and Z axes for the underlying grid */
|
|
||||||
/* i = [0 1 2] */
|
|
||||||
/* --> 0 0 0 : n will be computed and */
|
|
||||||
/* updated as output */
|
|
||||||
/* --> nx ny nz: these dimensions are used */
|
|
||||||
/*realin: structure defining a realization - */
|
|
||||||
/* must be a Gaussian white noise */
|
|
||||||
/*output: */
|
|
||||||
/*realout: structure defining a realization - */
|
|
||||||
|
|
||||||
void FFTMA2(struct vario_mod variogram,struct grid_mod grid,int n[3],struct realization_mod *realin,struct realization_mod *realout)
|
|
||||||
{
|
|
||||||
int NTOT,i,j,k,NMAX,NDIM,ntot,nmax,NXYZ,nxyz,maille0,maille1;
|
|
||||||
int solver;
|
|
||||||
double temp;
|
|
||||||
double *ireal,*covar,*workr,*worki,*realization;
|
|
||||||
string variable;
|
|
||||||
|
|
||||||
/*covariance axis normalization*/
|
|
||||||
axes(variogram.ap,variogram.scf,variogram.Nvario);
|
|
||||||
|
|
||||||
/*pseudo-grid definition*/
|
|
||||||
cgrid(variogram,grid,n);
|
|
||||||
|
|
||||||
/*constant definition*/
|
|
||||||
NTOT = n[0]*n[1]*n[2];
|
|
||||||
ntot = NTOT+1;
|
|
||||||
NMAX = n[0];
|
|
||||||
NDIM = 3;
|
|
||||||
for (i=1;i<3;i++) {
|
|
||||||
if (n[i] > NMAX)
|
|
||||||
NMAX = n[i];
|
|
||||||
if (n[i] == 1)
|
|
||||||
NDIM--;
|
|
||||||
}
|
|
||||||
nmax = NMAX+1;
|
|
||||||
NXYZ = grid.NX*grid.NY*grid.NZ;
|
|
||||||
nxyz = NXYZ+1;
|
|
||||||
|
|
||||||
/*array initialization*/
|
|
||||||
covar = (double *) malloc(ntot * sizeof(double));
|
|
||||||
/* testmemory(covar); */
|
|
||||||
allouememoire(covar,'covar');
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
ireal = (double *) malloc(ntot * sizeof(double));
|
|
||||||
testmemory(ireal);
|
|
||||||
|
|
||||||
realization = (double *) malloc(ntot * sizeof(double));
|
|
||||||
testmemory(realization);
|
|
||||||
|
|
||||||
workr = (double *) malloc(nmax * sizeof(double));
|
|
||||||
testmemory(workr);
|
|
||||||
|
|
||||||
worki = (double *) malloc(nmax * sizeof(double));
|
|
||||||
testmemory(worki);
|
|
||||||
|
|
||||||
/*covariance function creation*/
|
|
||||||
covariance(covar,variogram,grid,n);
|
|
||||||
|
|
||||||
/*power spectrum*/
|
|
||||||
fourt(covar,ireal,n,NDIM,1,0,workr,worki);
|
|
||||||
|
|
||||||
/*organization of the input Gaussian white noise*/
|
|
||||||
solver=0;
|
|
||||||
prebuild_gwn(grid,n,realin,realization,solver);
|
|
||||||
|
|
||||||
/*forward fourier transform of the GWN*/
|
|
||||||
fourt(realization,ireal,n,NDIM,1,0,workr,worki);
|
|
||||||
|
|
||||||
/* build realization in spectral domain */
|
|
||||||
build_real(n,NTOT,covar,realization,ireal);
|
|
||||||
|
|
||||||
free(covar);
|
|
||||||
|
|
||||||
/*backward fourier transform*/
|
|
||||||
fourt(realization,ireal,n,NDIM,0,1,workr,worki);
|
|
||||||
|
|
||||||
|
|
||||||
free(ireal);
|
|
||||||
free(workr);
|
|
||||||
free(worki);
|
|
||||||
|
|
||||||
/*output realization*/
|
|
||||||
clean_real(realin,n,grid,realization,realout);
|
|
||||||
|
|
||||||
free(realization);
|
|
||||||
|
|
||||||
return;
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
@ -1,162 +0,0 @@
|
|||||||
#include <stdlib.h>
|
|
||||||
#include <math.h>
|
|
||||||
#include "geostat.h"
|
|
||||||
#include "toolsludo.h"
|
|
||||||
|
|
||||||
|
|
||||||
/*FAST FOURIER TRANSFORM Pressure Simulation */
|
|
||||||
/*Turns a Gaussian white noise vector into a */
|
|
||||||
/*spatially correlated vector */
|
|
||||||
/*input: */
|
|
||||||
/*variogram: structure defining the variogram */
|
|
||||||
/* model */
|
|
||||||
/*grid: structure defining the grid */
|
|
||||||
/*n: vector with the number of cells along the */
|
|
||||||
/* X, Y and Z axes for the underlying grid */
|
|
||||||
/* i = [0 1 2] */
|
|
||||||
/* --> 0 0 0 : n will be computed and */
|
|
||||||
/* updated as output */
|
|
||||||
/* --> nx ny nz: these dimensions are used */
|
|
||||||
/*realin: structure defining a realization - */
|
|
||||||
/* must be a Gaussian white noise */
|
|
||||||
/*gradient: macroscopic gradient pression vector */
|
|
||||||
/*output: */
|
|
||||||
/*realout: structure defining a realization - */
|
|
||||||
/*realout2: structure defining a pressure field */
|
|
||||||
/*realout3: structure defining a xvelocity field */
|
|
||||||
/*realout4: structure defining a yvelocity field */
|
|
||||||
/*realout5: structure defining a zvelocity field */
|
|
||||||
|
|
||||||
|
|
||||||
void FFTPSim(struct vario_mod variogram,struct statistic_mod stat,struct grid_mod grid,int n[3],struct realization_mod *realin,struct pressure_mod gradient,struct realization_mod *realout,struct realization_mod *realout2,struct realization_mod *realout3,struct realization_mod *realout4,struct realization_mod *realout5)
|
|
||||||
|
|
||||||
{
|
|
||||||
int NTOT,i,j,k,NMAX,NDIM,ntot,nmax,NXYZ,nxyz,maille0,maille1;
|
|
||||||
double *workr,*worki,temp,temp2,coeff;
|
|
||||||
double *covar,*realization,*pressure;
|
|
||||||
double *icovar,*ireal,*ipressure;
|
|
||||||
double *xvelocity,*ixvelocity,*yvelocity,*iyvelocity,*zvelocity,*izvelocity;
|
|
||||||
double ki,kj,kk;
|
|
||||||
|
|
||||||
/*covariance axis normalization*/
|
|
||||||
axes(variogram.ap,variogram.scf,variogram.Nvario);
|
|
||||||
|
|
||||||
/*pseudo-grid definition*/
|
|
||||||
cgrid(variogram,grid,n);
|
|
||||||
|
|
||||||
/*constant definition*/
|
|
||||||
NTOT = n[0]*n[1]*n[2];
|
|
||||||
ntot = NTOT+1;
|
|
||||||
NMAX = n[0];
|
|
||||||
NDIM = 3;
|
|
||||||
for (i=1;i<3;i++) {
|
|
||||||
if (n[i] > NMAX)
|
|
||||||
NMAX = n[i];
|
|
||||||
if (n[i] == 1)
|
|
||||||
NDIM--;
|
|
||||||
}
|
|
||||||
nmax = NMAX+1;
|
|
||||||
NXYZ = grid.NX*grid.NY*grid.NZ;
|
|
||||||
nxyz = NXYZ+1;
|
|
||||||
|
|
||||||
/*array initialization*/
|
|
||||||
covar = (double *) malloc(ntot * sizeof(double));
|
|
||||||
testmemory(covar);
|
|
||||||
|
|
||||||
icovar = (double *) malloc(ntot * sizeof(double));
|
|
||||||
testmemory(icovar);
|
|
||||||
|
|
||||||
realization = (double *) malloc(ntot * sizeof(double));
|
|
||||||
testmemory(realization);
|
|
||||||
|
|
||||||
ireal = (double *) malloc(ntot * sizeof(double));
|
|
||||||
testmemory(ireal);
|
|
||||||
|
|
||||||
pressure = (double *) malloc(ntot * sizeof(double));
|
|
||||||
testmemory(pressure);
|
|
||||||
|
|
||||||
ipressure = (double *) malloc(ntot * sizeof(double));
|
|
||||||
testmemory(ipressure);
|
|
||||||
|
|
||||||
xvelocity = (double *) malloc(ntot * sizeof(double));
|
|
||||||
testmemory(xvelocity);
|
|
||||||
|
|
||||||
ixvelocity = (double *) malloc(ntot * sizeof(double));
|
|
||||||
testmemory(ixvelocity);
|
|
||||||
|
|
||||||
yvelocity = (double *) malloc(ntot * sizeof(double));
|
|
||||||
testmemory(yvelocity);
|
|
||||||
|
|
||||||
iyvelocity = (double *) malloc(ntot * sizeof(double));
|
|
||||||
testmemory(iyvelocity);
|
|
||||||
|
|
||||||
zvelocity = (double *) malloc(ntot * sizeof(double));
|
|
||||||
testmemory(zvelocity);
|
|
||||||
|
|
||||||
izvelocity = (double *) malloc(ntot * sizeof(double));
|
|
||||||
testmemory(izvelocity);
|
|
||||||
|
|
||||||
workr = (double *) malloc(nmax * sizeof(double));
|
|
||||||
testmemory(workr);
|
|
||||||
|
|
||||||
worki = (double *) malloc(nmax * sizeof(double));
|
|
||||||
testmemory(worki);
|
|
||||||
|
|
||||||
/*covariance function creation*/
|
|
||||||
covariance(covar,variogram,grid,n);
|
|
||||||
|
|
||||||
/*power spectrum*/
|
|
||||||
fourt(covar,icovar,n,NDIM,1,0,workr,worki);
|
|
||||||
|
|
||||||
/*organization of the input Gaussian white noise*/
|
|
||||||
|
|
||||||
prebuild_gwn(grid,n,realin,realization);
|
|
||||||
|
|
||||||
/*forward fourier transform of the GWN*/
|
|
||||||
fourt(realization,ireal,n,NDIM,1,0,workr,worki);
|
|
||||||
|
|
||||||
/* build realization in spectral domain */
|
|
||||||
build_real(n,NTOT,covar,realization,ireal);
|
|
||||||
|
|
||||||
free(covar);
|
|
||||||
|
|
||||||
/* pressure calculation in the spectral domain*/
|
|
||||||
|
|
||||||
build_pressure(n,grid,gradient,realization,ireal,pressure,ipressure);
|
|
||||||
build_velocity(n,grid,stat,gradient,realization,ireal,xvelocity,ixvelocity,1);
|
|
||||||
build_velocity(n,grid,stat,gradient,realization,ireal,yvelocity,iyvelocity,2);
|
|
||||||
build_velocity(n,grid,stat,gradient,realization,ireal,zvelocity,izvelocity,3);
|
|
||||||
|
|
||||||
/*backward fourier transform*/
|
|
||||||
fourt(realization,ireal,n,NDIM,0,1,workr,worki);
|
|
||||||
fourt(pressure,ipressure,n,NDIM,0,1,workr,worki);
|
|
||||||
fourt(xvelocity,ixvelocity,n,NDIM,0,1,workr,worki);
|
|
||||||
fourt(yvelocity,iyvelocity,n,NDIM,0,1,workr,worki);
|
|
||||||
fourt(zvelocity,izvelocity,n,NDIM,0,1,workr,worki);
|
|
||||||
|
|
||||||
free(ireal);
|
|
||||||
free(ipressure);
|
|
||||||
free(ixvelocity);
|
|
||||||
free(iyvelocity);
|
|
||||||
free(izvelocity);
|
|
||||||
free(workr);
|
|
||||||
free(worki);
|
|
||||||
|
|
||||||
/*output realization*/
|
|
||||||
/*is the output realization already allocated?*/
|
|
||||||
/*if not, memory allocation*/
|
|
||||||
|
|
||||||
clean_real(realin,n,grid,realization,realout);
|
|
||||||
clean_real(realin,n,grid,pressure,realout2);
|
|
||||||
clean_real(realin,n,grid,xvelocity,realout3);
|
|
||||||
clean_real(realin,n,grid,yvelocity,realout4);
|
|
||||||
clean_real(realin,n,grid,zvelocity,realout5);
|
|
||||||
|
|
||||||
free(realization);
|
|
||||||
free(pressure);
|
|
||||||
free(xvelocity);
|
|
||||||
free(yvelocity);
|
|
||||||
free(zvelocity);
|
|
||||||
|
|
||||||
return;
|
|
||||||
}
|
|
@ -1,157 +0,0 @@
|
|||||||
#include <stdio.h>
|
|
||||||
#include <stddef.h>
|
|
||||||
#include <string.h>
|
|
||||||
#include <stdarg.h>
|
|
||||||
#include <stdlib.h>
|
|
||||||
#include <math.h>
|
|
||||||
#include "geostat.h"
|
|
||||||
#include "pressure.h"
|
|
||||||
|
|
||||||
|
|
||||||
/*FAST FOURIER TRANSFORM Pressure Simulation */
|
|
||||||
/*Turns a Gaussian white noise vector into a */
|
|
||||||
/*spatially correlated vector */
|
|
||||||
/*input: */
|
|
||||||
/*variogram: structure defining the variogram */
|
|
||||||
/* model */
|
|
||||||
/*grid: structure defining the grid */
|
|
||||||
/*n: vector with the number of cells along the */
|
|
||||||
/* X, Y and Z axes for the underlying grid */
|
|
||||||
/* i = [0 1 2] */
|
|
||||||
/* --> 0 0 0 : n will be computed and */
|
|
||||||
/* updated as output */
|
|
||||||
/* --> nx ny nz: these dimensions are used */
|
|
||||||
/*realin: structure defining a realization - */
|
|
||||||
/* must be a Gaussian white noise */
|
|
||||||
/*gradient: macroscopic gradient pression vector */
|
|
||||||
/*output: */
|
|
||||||
/*realout: structure defining a realization - */
|
|
||||||
/*realout2: structure defining a pressure field */
|
|
||||||
/*realout3: structure defining a xvelocity field */
|
|
||||||
/*realout4: structure defining a yvelocity field */
|
|
||||||
/*realout5: structure defining a zvelocity field */
|
|
||||||
|
|
||||||
|
|
||||||
void FFTPressure(int n[3],struct grid_mod grid,struct realization_mod *realin,struct statistic_mod stat,struct pressure_mod gradient,struct realization_mod *realout,struct realization_mod *realout2,struct realization_mod *realout3,struct realization_mod *realout4,int solver)
|
|
||||||
|
|
||||||
{
|
|
||||||
int NTOT,i,j,k,NMAX,NDIM,ntot,nmax,NXYZ,nxyz,maille0,maille1;
|
|
||||||
double *workr,*worki,temp,temp2,coeff;
|
|
||||||
double *realization,*pressure;
|
|
||||||
double *ireal,*ipressure;
|
|
||||||
double *xvelocity,*ixvelocity,*yvelocity,*iyvelocity,*zvelocity,*izvelocity;
|
|
||||||
double ki,kj,kk;
|
|
||||||
FILE *fp;
|
|
||||||
/* string nomfichier; */
|
|
||||||
|
|
||||||
/*constant definition*/
|
|
||||||
NTOT = n[0]*n[1]*n[2];
|
|
||||||
ntot = NTOT+1;
|
|
||||||
NMAX = n[0];
|
|
||||||
NDIM = 3;
|
|
||||||
for (i=1;i<3;i++) {
|
|
||||||
if (n[i] > NMAX)
|
|
||||||
NMAX = n[i];
|
|
||||||
if (n[i] == 1)
|
|
||||||
NDIM--;
|
|
||||||
}
|
|
||||||
nmax = NMAX+1;
|
|
||||||
NXYZ = grid.NX*grid.NY*grid.NZ;
|
|
||||||
nxyz = NXYZ+1;
|
|
||||||
|
|
||||||
realization = (double *) malloc(ntot * sizeof(double));
|
|
||||||
testmemory(realization);
|
|
||||||
|
|
||||||
ireal = (double *) malloc(ntot * sizeof(double));
|
|
||||||
testmemory(ireal);
|
|
||||||
|
|
||||||
pressure = (double *) malloc(ntot * sizeof(double));
|
|
||||||
testmemory(pressure);
|
|
||||||
|
|
||||||
ipressure = (double *) malloc(ntot * sizeof(double));
|
|
||||||
testmemory(ipressure);
|
|
||||||
|
|
||||||
xvelocity = (double *) malloc(ntot * sizeof(double));
|
|
||||||
testmemory(xvelocity);
|
|
||||||
|
|
||||||
ixvelocity = (double *) malloc(ntot * sizeof(double));
|
|
||||||
testmemory(ixvelocity);
|
|
||||||
|
|
||||||
yvelocity = (double *) malloc(ntot * sizeof(double));
|
|
||||||
testmemory(yvelocity);
|
|
||||||
|
|
||||||
iyvelocity = (double *) malloc(ntot * sizeof(double));
|
|
||||||
testmemory(iyvelocity);
|
|
||||||
|
|
||||||
zvelocity = (double *) malloc(ntot * sizeof(double));
|
|
||||||
testmemory(zvelocity);
|
|
||||||
|
|
||||||
izvelocity = (double *) malloc(ntot * sizeof(double));
|
|
||||||
testmemory(izvelocity);
|
|
||||||
|
|
||||||
workr = (double *) malloc(nmax * sizeof(double));
|
|
||||||
testmemory(workr);
|
|
||||||
|
|
||||||
worki = (double *) malloc(nmax * sizeof(double));
|
|
||||||
testmemory(worki);
|
|
||||||
|
|
||||||
/*organization of the realization*/
|
|
||||||
prebuild_gwn(grid,n,realin,realization,solver);
|
|
||||||
|
|
||||||
/*forward fourier transform of the GWN*/
|
|
||||||
fourt(realization,ireal,n,NDIM,1,0,workr,worki);
|
|
||||||
|
|
||||||
/* pressure calculation in the spectral domain*/
|
|
||||||
|
|
||||||
build_pressure(n,grid,gradient,realization,ireal,pressure,ipressure);
|
|
||||||
|
|
||||||
build_velocity(n,grid,stat,gradient,realization,ireal,xvelocity,ixvelocity,1);
|
|
||||||
build_velocity(n,grid,stat,gradient,realization,ireal,yvelocity,iyvelocity,2);
|
|
||||||
build_velocity(n,grid,stat,gradient,realization,ireal,zvelocity,izvelocity,3);
|
|
||||||
|
|
||||||
/*backward fourier transform*/
|
|
||||||
fourt(realization,ireal,n,NDIM,0,1,workr,worki);
|
|
||||||
fourt(pressure,ipressure,n,NDIM,0,1,workr,worki);
|
|
||||||
fourt(xvelocity,ixvelocity,n,NDIM,0,1,workr,worki);
|
|
||||||
fourt(yvelocity,iyvelocity,n,NDIM,0,1,workr,worki);
|
|
||||||
fourt(zvelocity,izvelocity,n,NDIM,0,1,workr,worki);
|
|
||||||
|
|
||||||
free(ireal);
|
|
||||||
free(ipressure);
|
|
||||||
free(ixvelocity);
|
|
||||||
free(iyvelocity);
|
|
||||||
free(izvelocity);
|
|
||||||
free(workr);
|
|
||||||
free(worki);
|
|
||||||
|
|
||||||
for (i=1;i<=NTOT;i++)
|
|
||||||
realization[i]=realization[i]/(double) NTOT;
|
|
||||||
|
|
||||||
fp = fopen("realization.test", "w");
|
|
||||||
for (i=1;i<=NTOT;i++)
|
|
||||||
fprintf(fp,"%f\n",realization[i]);
|
|
||||||
fclose(fp);
|
|
||||||
|
|
||||||
/* nomfichier="pression.test"; */
|
|
||||||
fp = fopen("pression.test", "w");
|
|
||||||
for (i=1;i<=NTOT;i++)
|
|
||||||
fprintf(fp,"%f\n",pressure[i]);
|
|
||||||
fclose(fp);
|
|
||||||
|
|
||||||
/*output realization*/
|
|
||||||
/*is the output realization already allocated?*/
|
|
||||||
/*if not, memory allocation*/
|
|
||||||
|
|
||||||
clean_real2(realin,n,grid,solver,pressure,realout);
|
|
||||||
clean_real2(realin,n,grid,solver,xvelocity,realout2);
|
|
||||||
clean_real2(realin,n,grid,solver,yvelocity,realout3);
|
|
||||||
clean_real2(realin,n,grid,solver,zvelocity,realout4);
|
|
||||||
|
|
||||||
/* free(realization); */
|
|
||||||
/* free(pressure); */
|
|
||||||
/* free(xvelocity); */
|
|
||||||
/* free(yvelocity); */
|
|
||||||
/* free(zvelocity); */
|
|
||||||
|
|
||||||
return;
|
|
||||||
}
|
|
@ -1,68 +0,0 @@
|
|||||||
#include <stdio.h>
|
|
||||||
#include <stddef.h>
|
|
||||||
#include <string.h>
|
|
||||||
#include <stdarg.h>
|
|
||||||
#include <stdlib.h>
|
|
||||||
#include <math.h>
|
|
||||||
#include "geostat.h"
|
|
||||||
|
|
||||||
|
|
||||||
/*FAST FOURIER TRANSFORM Test */
|
|
||||||
|
|
||||||
void FFTtest(int n[3],struct grid_mod grid,struct realization_mod *realin,struct realization_mod *realout,int solver)
|
|
||||||
|
|
||||||
{
|
|
||||||
int NTOT,i,j,k,NMAX,NDIM,ntot,nmax,NXYZ,nxyz,maille0,maille1;
|
|
||||||
double *workr,*worki,temp,temp2,coeff;
|
|
||||||
double *realization;
|
|
||||||
double *ireal;
|
|
||||||
double ki,kj,kk;
|
|
||||||
FILE *fp;
|
|
||||||
|
|
||||||
/*constant definition*/
|
|
||||||
NTOT = n[0]*n[1]*n[2];
|
|
||||||
ntot = NTOT+1;
|
|
||||||
NMAX = n[0];
|
|
||||||
NDIM = 3;
|
|
||||||
for (i=1;i<3;i++) {
|
|
||||||
if (n[i] > NMAX)
|
|
||||||
NMAX = n[i];
|
|
||||||
if (n[i] == 1)
|
|
||||||
NDIM--;
|
|
||||||
}
|
|
||||||
nmax = NMAX+1;
|
|
||||||
NXYZ = grid.NX*grid.NY*grid.NZ;
|
|
||||||
nxyz = NXYZ+1;
|
|
||||||
|
|
||||||
realization = (double *) malloc(ntot * sizeof(double));
|
|
||||||
testmemory(realization);
|
|
||||||
|
|
||||||
ireal = (double *) malloc(ntot * sizeof(double));
|
|
||||||
testmemory(ireal);
|
|
||||||
|
|
||||||
workr = (double *) malloc(nmax * sizeof(double));
|
|
||||||
testmemory(workr);
|
|
||||||
|
|
||||||
worki = (double *) malloc(nmax * sizeof(double));
|
|
||||||
testmemory(worki);
|
|
||||||
|
|
||||||
/*organization of the realization*/
|
|
||||||
prebuild_gwn(grid,n,realin,realization,solver);
|
|
||||||
|
|
||||||
fp = fopen("perm.test", "w");
|
|
||||||
for (i=1;i<=NTOT;i++)
|
|
||||||
fprintf(fp,"%f\n",realization[i]);
|
|
||||||
fclose(fp);
|
|
||||||
|
|
||||||
/*forward fourier transform of the GWN*/
|
|
||||||
fourt(realization,ireal,n,NDIM,1,0,workr,worki);
|
|
||||||
|
|
||||||
/*backward fourier transform*/
|
|
||||||
fourt(realization,ireal,n,NDIM,0,1,workr,worki);
|
|
||||||
|
|
||||||
fp = fopen("perm2.test", "w");
|
|
||||||
for (i=1;i<=NTOT;i++)
|
|
||||||
fprintf(fp,"%f\n",realization[i]);
|
|
||||||
fclose(fp);
|
|
||||||
return;
|
|
||||||
}
|
|
@ -1,13 +0,0 @@
|
|||||||
CC = cc
|
|
||||||
CCFLAG =
|
|
||||||
|
|
||||||
LIB = libCS106X_${ARCH}.a
|
|
||||||
%.o: %.c
|
|
||||||
$(CC) $(CCFLAG) -c $<
|
|
||||||
|
|
||||||
|
|
||||||
$(LIB) : genlib.o random.o simpio.o strlib.o symtab.o scanadt.o stack.o
|
|
||||||
ar r $@ $?
|
|
||||||
|
|
||||||
clean :
|
|
||||||
rm *.o
|
|
@ -1,30 +0,0 @@
|
|||||||
|
|
||||||
########################################################################
|
|
||||||
#
|
|
||||||
# Makefile for library
|
|
||||||
#
|
|
||||||
########################################################################
|
|
||||||
|
|
||||||
|
|
||||||
INTERFACE = ../include
|
|
||||||
INCLUDE = -I${INTERFACE}
|
|
||||||
LIBS = -lm -L../lib
|
|
||||||
CC = cc
|
|
||||||
CCFLAGS =
|
|
||||||
|
|
||||||
LIB = libFFTMA_${ARCH}.a
|
|
||||||
|
|
||||||
NOBJS= gammf.o fftma.o addstat.o axes.o cgrid.o covariance.o fourt.o length.o maxfactor.o test_fact.o cov_value.o generate.o gasdev.o ran2.o stable.o gaussian.o power.o cubic.o spherical.o nugget.o exponential.o cardsin.o nor2log.o
|
|
||||||
|
|
||||||
|
|
||||||
.c.o:
|
|
||||||
${CC} $(INCLUDE) $(CCFLAGS) -c $<
|
|
||||||
|
|
||||||
# LIBRARY
|
|
||||||
$(LIB) : $(NOBJS)
|
|
||||||
ar cr $(LIB) $(NOBJS)
|
|
||||||
ranlib $(LIB)
|
|
||||||
|
|
||||||
clean :
|
|
||||||
rm *.o
|
|
||||||
|
|
@ -1,30 +0,0 @@
|
|||||||
|
|
||||||
########################################################################
|
|
||||||
#
|
|
||||||
# Makefile for library
|
|
||||||
#
|
|
||||||
########################################################################
|
|
||||||
|
|
||||||
|
|
||||||
INTERFACE = ../include
|
|
||||||
|
|
||||||
INCLUDE = -I${INTERFACE}
|
|
||||||
LIBS = -lm -L../lib
|
|
||||||
CC = cc
|
|
||||||
CCFLAGS =
|
|
||||||
|
|
||||||
LIB = libFFTMA2_${ARCH}.a
|
|
||||||
|
|
||||||
NOBJS= kgeneration.o kgeneration2.o fftma2.o prebuild_gwn.o build_real.o addstat2.o clean_real.o
|
|
||||||
|
|
||||||
.c.o:
|
|
||||||
${CC} $(INCLUDE) $(CCFLAGS) -c $<
|
|
||||||
|
|
||||||
# LIBRARY
|
|
||||||
$(LIB) : $(NOBJS)
|
|
||||||
ar cr $(LIB) $(NOBJS)
|
|
||||||
ranlib $(LIB)
|
|
||||||
|
|
||||||
clean :
|
|
||||||
rm *.o
|
|
||||||
|
|
@ -1,34 +0,0 @@
|
|||||||
|
|
||||||
########################################################################
|
|
||||||
#
|
|
||||||
# Makefile for library
|
|
||||||
#
|
|
||||||
########################################################################
|
|
||||||
|
|
||||||
|
|
||||||
#INTERFACE = ${BOSSE}/CODES/LIBS_C/Interface
|
|
||||||
#INTERFACE2= ${BOSSE}/CODES/Geostat/FFTPSim/Version_binaire2/src/Interface
|
|
||||||
|
|
||||||
INTERFACE = ../include
|
|
||||||
#INTERFACE2= ../Interface
|
|
||||||
|
|
||||||
INCLUDE = -I${INTERFACE}
|
|
||||||
LIBS = -lm -L../lib
|
|
||||||
CCFLAGS =
|
|
||||||
CC = cc
|
|
||||||
|
|
||||||
LIB = libFFTPSIM_${ARCH}.a
|
|
||||||
|
|
||||||
NOBJS= pgeneration.o pgeneration2.o FFTPressure.o FFTtest.o build_pressure.o build_velocity.o total_pressure.o total_velocity.o clean_real2.o waveVectorCompute3D.o mat_vec.o derivReal.o
|
|
||||||
|
|
||||||
.c.o:
|
|
||||||
${CC} $(INCLUDE) $(CCFLAGS) -c $<
|
|
||||||
|
|
||||||
# LIBRARY
|
|
||||||
$(LIB) : $(NOBJS)
|
|
||||||
ar cr $(LIB) $(NOBJS)
|
|
||||||
ranlib $(LIB)
|
|
||||||
|
|
||||||
clean :
|
|
||||||
rm *.o
|
|
||||||
|
|
@ -1,31 +0,0 @@
|
|||||||
#
|
|
||||||
########################################################################
|
|
||||||
#
|
|
||||||
# Makefile for library
|
|
||||||
#
|
|
||||||
########################################################################
|
|
||||||
#
|
|
||||||
|
|
||||||
#INTERFACE = ${BOSSE}/CODES/LIBS_C/Interface
|
|
||||||
INTERFACE = ../include
|
|
||||||
|
|
||||||
INCLUDE = -I${INTERFACE}
|
|
||||||
CC = cc
|
|
||||||
CCFLAGS =
|
|
||||||
|
|
||||||
|
|
||||||
LIB = libIO_${ARCH}.a
|
|
||||||
|
|
||||||
NOBJS= inputdata.o inputfiledata.o debuginput.o readdata.o readfile_bin.o writefile.o writefile_bin.o testmemory.o testopenfile.o readdata2.o readdata3.o
|
|
||||||
|
|
||||||
.c.o:
|
|
||||||
${CC} $(INCLUDE) $(CCFLAGS) -c $<
|
|
||||||
|
|
||||||
# LIBRARY
|
|
||||||
$(LIB) : $(NOBJS)
|
|
||||||
ar cr $(LIB) $(NOBJS)
|
|
||||||
ranlib $(LIB)
|
|
||||||
|
|
||||||
clean :
|
|
||||||
rm *.o
|
|
||||||
|
|
@ -1,21 +0,0 @@
|
|||||||
# CC = gcc
|
|
||||||
# CCFLAG =
|
|
||||||
|
|
||||||
ifeq (${ARCH},sun)
|
|
||||||
CC = cc
|
|
||||||
CCFLAG = -g
|
|
||||||
endif
|
|
||||||
ifeq (${ARCH},linux)
|
|
||||||
CC = gcc
|
|
||||||
CCFLAG = -g
|
|
||||||
endif
|
|
||||||
LIB = libCS106X_debug_${ARCH}.a
|
|
||||||
%.o: %.c
|
|
||||||
$(CC) $(CCFLAG) -c $<
|
|
||||||
|
|
||||||
|
|
||||||
$(LIB) : genlib.o random.o simpio.o strlib.o symtab.o scanadt.o stack.o
|
|
||||||
ar r $@ $?
|
|
||||||
|
|
||||||
clean :
|
|
||||||
rm *.o
|
|
@ -1,106 +0,0 @@
|
|||||||
#include <stdio.h>
|
|
||||||
#include <math.h>
|
|
||||||
#include "geostat.h"
|
|
||||||
#include "genlib.h"
|
|
||||||
|
|
||||||
|
|
||||||
/*addstat */
|
|
||||||
/*adds mean and variance effects to the N(0,1) realization*/
|
|
||||||
/*input: */
|
|
||||||
/*realin: structure defining a realization - */
|
|
||||||
/* must have zeio mean and unit variance */
|
|
||||||
/*stat: structure defining the mean and variance */
|
|
||||||
/*output: */
|
|
||||||
/*realout: structure defining a realization - */
|
|
||||||
void addstat(struct realization_mod *realin,struct statistic_mod stat ,struct realization_mod *realout)
|
|
||||||
|
|
||||||
{
|
|
||||||
|
|
||||||
int i,nblockm,nblockv;
|
|
||||||
double r,moy,var;
|
|
||||||
|
|
||||||
/*Is the output realization allocated ?*/
|
|
||||||
/*is its length equal to realin.n?*/
|
|
||||||
if ((*realout).vector == NULL || (*realout).n != (*realin).n) {
|
|
||||||
(*realout).vector = (double *) malloc((*realin).n * sizeof(double));
|
|
||||||
if ((*realout).vector == NULL)
|
|
||||||
Error("No memory available");
|
|
||||||
}
|
|
||||||
(*realout).n = (*realin).n;
|
|
||||||
|
|
||||||
|
|
||||||
/*test over the input realization*/
|
|
||||||
switch ((*realin).code) {
|
|
||||||
case 0:
|
|
||||||
case 1:
|
|
||||||
(*realout).code = 2;
|
|
||||||
break;
|
|
||||||
case 6:
|
|
||||||
(*realout).code = 7;
|
|
||||||
break;
|
|
||||||
default:
|
|
||||||
(*realout).code = (*realin).code;
|
|
||||||
break;
|
|
||||||
}
|
|
||||||
|
|
||||||
for (i = 0; i < (*realin).n; i++) {
|
|
||||||
|
|
||||||
|
|
||||||
/*mean*/
|
|
||||||
switch (stat.nblock_mean) {
|
|
||||||
case 1:
|
|
||||||
/*stationary case*/
|
|
||||||
nblockm = 1;
|
|
||||||
break;
|
|
||||||
default:
|
|
||||||
/*number of the cell - nonstationary case*/
|
|
||||||
nblockm = i+1;
|
|
||||||
break;
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
/*variance*/
|
|
||||||
switch (stat.nblock_var) {
|
|
||||||
case 1:
|
|
||||||
/*stationary case*/
|
|
||||||
nblockv = 1;
|
|
||||||
break;
|
|
||||||
default:
|
|
||||||
/*number of the cell - nonstationary case*/
|
|
||||||
nblockv = i+1;
|
|
||||||
break;
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
switch (stat.type) {
|
|
||||||
case 0:
|
|
||||||
/*normal case*/
|
|
||||||
moy = stat.mean[nblockm-1];
|
|
||||||
var = stat.variance[nblockv-1];
|
|
||||||
break;
|
|
||||||
case 1:
|
|
||||||
/*lognormal (natural) case*/
|
|
||||||
r = (double)sqrt(stat.variance[nblockv-1])/stat.mean[nblockm-1];
|
|
||||||
r *= r;
|
|
||||||
moy = (double)log(stat.mean[nblockm-1]/sqrt(1.0+r));
|
|
||||||
var = (double)log(1.0+r);
|
|
||||||
break;
|
|
||||||
case 2:
|
|
||||||
/*lognormal (log10) case*/
|
|
||||||
r = (double)sqrt(stat.variance[nblockv-1])/stat.mean[nblockm-1];
|
|
||||||
r *= r;
|
|
||||||
moy = (double)log10(stat.mean[nblockm-1]/sqrt(1.0+r));
|
|
||||||
var = (double)log10(1.0+r);
|
|
||||||
break;
|
|
||||||
default:
|
|
||||||
Error("Type not defined in addstat");
|
|
||||||
break;
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
(*realout).vector[i] = (double)sqrt(var)*(*realin).vector[i]+moy;
|
|
||||||
}
|
|
||||||
|
|
||||||
return;
|
|
||||||
}
|
|
@ -1,119 +0,0 @@
|
|||||||
#include <stdio.h>
|
|
||||||
#include <math.h>
|
|
||||||
#include "genlib.h"
|
|
||||||
#include "geostat.h"
|
|
||||||
|
|
||||||
/*addstat */
|
|
||||||
/*adds mean and variance effects to the N(0,1) realization*/
|
|
||||||
/*input: */
|
|
||||||
/*realin: structure defining a realization - */
|
|
||||||
/* must have zeio mean and unit variance */
|
|
||||||
/*stat: structure defining the mean and variance */
|
|
||||||
/*output: */
|
|
||||||
/*realout: structure defining a realization - */
|
|
||||||
|
|
||||||
void addstat2(struct realization_mod *realin,struct statistic_mod stat ,struct realization_mod *realout,struct realization_mod *realout2)
|
|
||||||
|
|
||||||
{
|
|
||||||
|
|
||||||
int i,nblockm,nblockv;
|
|
||||||
double r,moy,var;
|
|
||||||
|
|
||||||
/*Is the output realization allocated ?*/
|
|
||||||
/*is its length equal to realin.n?*/
|
|
||||||
if ((*realout).vector == NULL || (*realout).n != (*realin).n) {
|
|
||||||
(*realout).vector = (double *) malloc((*realin).n * sizeof(double));
|
|
||||||
if ((*realout).vector == NULL)
|
|
||||||
Error("No memory available");
|
|
||||||
}
|
|
||||||
(*realout).n = (*realin).n;
|
|
||||||
|
|
||||||
if ((*realout2).vector == NULL || (*realout2).n != (*realin).n) {
|
|
||||||
(*realout2).vector = (double *) malloc((*realin).n * sizeof(double));
|
|
||||||
if ((*realout2).vector == NULL)
|
|
||||||
Error("No memory available");
|
|
||||||
}
|
|
||||||
(*realout2).n = (*realin).n;
|
|
||||||
|
|
||||||
|
|
||||||
/*test over the input realization*/
|
|
||||||
switch ((*realin).code) {
|
|
||||||
case 0:
|
|
||||||
case 1:
|
|
||||||
(*realout).code = 2;
|
|
||||||
(*realout2).code = 2;
|
|
||||||
break;
|
|
||||||
case 6:
|
|
||||||
(*realout).code = 7;
|
|
||||||
(*realout2).code = 7;
|
|
||||||
break;
|
|
||||||
default:
|
|
||||||
(*realout).code = (*realin).code;
|
|
||||||
(*realout2).code = (*realin).code;
|
|
||||||
break;
|
|
||||||
}
|
|
||||||
|
|
||||||
for (i = 0; i < (*realin).n; i++) {
|
|
||||||
|
|
||||||
|
|
||||||
/*mean*/
|
|
||||||
switch (stat.nblock_mean) {
|
|
||||||
case 1:
|
|
||||||
/*stationary case*/
|
|
||||||
nblockm = 1;
|
|
||||||
break;
|
|
||||||
default:
|
|
||||||
/*number of the cell - nonstationary case*/
|
|
||||||
nblockm = i+1;
|
|
||||||
break;
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
/*variance*/
|
|
||||||
switch (stat.nblock_var) {
|
|
||||||
case 1:
|
|
||||||
/*stationary case*/
|
|
||||||
nblockv = 1;
|
|
||||||
break;
|
|
||||||
default:
|
|
||||||
/*number of the cell - nonstationary case*/
|
|
||||||
nblockv = i+1;
|
|
||||||
break;
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
/* switch (stat.type) { */
|
|
||||||
/* case 0: */
|
|
||||||
|
|
||||||
/*normal case*/
|
|
||||||
moy = stat.mean[nblockm-1];
|
|
||||||
var = stat.variance[nblockv-1];
|
|
||||||
|
|
||||||
/* break; */
|
|
||||||
/* case 1: */
|
|
||||||
/*lognormal (natural) case*/
|
|
||||||
/* r = (double)sqrt(stat.variance[nblockv-1])/stat.mean[nblockm-1]; */
|
|
||||||
/* r *= r; */
|
|
||||||
/* moy = (double)log(stat.mean[nblockm-1]/sqrt(1.0+r)); */
|
|
||||||
/* var = (double)log(1.0+r); */
|
|
||||||
/* break; */
|
|
||||||
/* case 2: */
|
|
||||||
/*lognormal (log10) case*/
|
|
||||||
/* r = (double)sqrt(stat.variance[nblockv-1])/stat.mean[nblockm-1]; */
|
|
||||||
/* r *= r; */
|
|
||||||
/* moy = (double)log10(stat.mean[nblockm-1]/sqrt(1.0+r)); */
|
|
||||||
/* var = (double)log10(1.0+r); */
|
|
||||||
/* break; */
|
|
||||||
/* default: */
|
|
||||||
/* Error("Type not defined in addstat"); */
|
|
||||||
/* break; */
|
|
||||||
/* } */
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
(*realout).vector[i] = (double)sqrt(var)*(*realin).vector[i];
|
|
||||||
(*realout2).vector[i] = (double)sqrt(var)*(*realin).vector[i]+moy;
|
|
||||||
}
|
|
||||||
|
|
||||||
return;
|
|
||||||
}
|
|
@ -1,16 +0,0 @@
|
|||||||
#include <stdio.h>
|
|
||||||
#include <stddef.h>
|
|
||||||
#include <string.h>
|
|
||||||
#include <stdarg.h>
|
|
||||||
#include <stdlib.h>
|
|
||||||
#include <math.h>
|
|
||||||
|
|
||||||
void allouememoire(double *realint, string variable)
|
|
||||||
{
|
|
||||||
|
|
||||||
if (realint == NULL) {
|
|
||||||
printf("Testmemory.c: No memory available for %s\n",variable);
|
|
||||||
exit;
|
|
||||||
}
|
|
||||||
return;
|
|
||||||
}
|
|
@ -1,50 +0,0 @@
|
|||||||
#include <stdio.h>
|
|
||||||
#include <stddef.h>
|
|
||||||
#include <string.h>
|
|
||||||
#include <stdarg.h>
|
|
||||||
#include <stdlib.h>
|
|
||||||
#include <math.h>
|
|
||||||
#include "geostat.h"
|
|
||||||
|
|
||||||
void clean_real2(struct realization_mod *realin,int n[3],struct grid_mod grid,int solver,double *vectorresult,struct realization_mod *realout)
|
|
||||||
{
|
|
||||||
int i,j,k,maille0,maille1;
|
|
||||||
double NTOT;
|
|
||||||
|
|
||||||
NTOT=n[0]*n[1]*n[2];
|
|
||||||
/*is the output realization already allocated?*/
|
|
||||||
/*if not, memory allocation*/
|
|
||||||
|
|
||||||
if ((*realout).vector == NULL || (*realout).n != (*realin).n)
|
|
||||||
{
|
|
||||||
(*realout).vector = (double *) malloc((*realin).n * sizeof(double));
|
|
||||||
if ((*realout).vector == NULL)
|
|
||||||
{
|
|
||||||
printf("Clean_real2.c: No memory available\n");
|
|
||||||
exit;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
(*realout).n = (*realin).n;
|
|
||||||
(*realout).code = 1;
|
|
||||||
if (solver==1)
|
|
||||||
{
|
|
||||||
for ( k = 1; k <= NTOT;k++)
|
|
||||||
{
|
|
||||||
(*realout).vector[k-1] = vectorresult[k]/(double) NTOT;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
else
|
|
||||||
{
|
|
||||||
for ( k = 1; k <= grid.NZ; k++) {
|
|
||||||
for (j = 1; j <= grid.NY; j++) {
|
|
||||||
for (i = 1; i <= grid.NX; i++) {
|
|
||||||
maille1 = i+(j-1+(k-1)*n[1])*n[0];
|
|
||||||
maille0 = i-1+(j-1+(k-1)*grid.NY)*grid.NX;
|
|
||||||
(*realout).vector[maille0] = vectorresult[maille1]/(double) NTOT;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
return;
|
|
||||||
}
|
|
@ -1,216 +0,0 @@
|
|||||||
#include <math.h>
|
|
||||||
#include <stdlib.h>
|
|
||||||
#include <stdio.h>
|
|
||||||
#include "condor.h"
|
|
||||||
#include "geostat.h"
|
|
||||||
|
|
||||||
/* Private functions */
|
|
||||||
void normAxes(double *vec, double *normed);
|
|
||||||
|
|
||||||
|
|
||||||
void derivReal(struct realization_mod *Z, struct realization_mod *dZ, double *dir, struct grid_mod *grid, struct vario_mod vario) {
|
|
||||||
|
|
||||||
int n[3],i,j,k,maille,extMaille;
|
|
||||||
int NTOT,ntot,NMAX,nmax,NDIM;
|
|
||||||
// int NXYZ,nxyz;
|
|
||||||
double nDir[3];
|
|
||||||
double *table,*workr,*worki,*realization,*waveVect;
|
|
||||||
|
|
||||||
|
|
||||||
int Ngrid;
|
|
||||||
double tmp;
|
|
||||||
double *ExtendedWaveVect;
|
|
||||||
|
|
||||||
/* Test the input real*/
|
|
||||||
/* if ((*Z).code != 1) { */
|
|
||||||
/* printf("Realization should be Standard Normal\n"); */
|
|
||||||
/* return; */
|
|
||||||
/* } */
|
|
||||||
|
|
||||||
printf("grid.Nx = %d\n",(*grid).NX);
|
|
||||||
printf("grid.Ny = %d\n",(*grid).NY);
|
|
||||||
printf("grid.Nz = %d\n",(*grid).NZ);
|
|
||||||
|
|
||||||
printf("vario.Nvario = %d\n",vario.Nvario);
|
|
||||||
for(i=0; i<vario.Nvario;i++) {
|
|
||||||
printf("Component %d:\n",i);
|
|
||||||
printf("Type: %d\n",vario.vario[i]);
|
|
||||||
printf("Exponent: %f\n",vario.alpha[i]);
|
|
||||||
printf("Anisotropy axes:");
|
|
||||||
for(j=0; j<6; j++) {
|
|
||||||
printf(" %f",vario.ap[6*i+j]);
|
|
||||||
}
|
|
||||||
printf("\nCorrelation length:\n");
|
|
||||||
for(j=0;j<3;j++) {
|
|
||||||
printf("axe%d :%f\n",j+1,vario.scf[3*i+j]);
|
|
||||||
}
|
|
||||||
printf("Variance: %f\n",vario.var[i]);
|
|
||||||
}
|
|
||||||
|
|
||||||
/* covariance axis normalisation */
|
|
||||||
axes(vario.ap,vario.scf,vario.Nvario);
|
|
||||||
|
|
||||||
n[0]=0;
|
|
||||||
n[1]=0;
|
|
||||||
n[2]=0;
|
|
||||||
|
|
||||||
/* pseudo-grid definition */
|
|
||||||
cgrid(vario,*grid,n);
|
|
||||||
|
|
||||||
printf("n[0] = %d\n", n[0]);
|
|
||||||
printf("n[1] = %d\n", n[1]);
|
|
||||||
printf("n[2] = %d\n", n[2]);
|
|
||||||
|
|
||||||
/* constants definition */
|
|
||||||
NTOT = n[0]*n[1]*n[2];
|
|
||||||
ntot = NTOT + 1;
|
|
||||||
NDIM = 3;
|
|
||||||
NMAX = n[0];
|
|
||||||
|
|
||||||
Ngrid = (*grid).NX*(*grid).NY*(*grid).NZ;
|
|
||||||
|
|
||||||
for (i=1;i<3;i++) {
|
|
||||||
if (n[i] > NMAX) NMAX = n[i];
|
|
||||||
if (n[i] == 1) NDIM = NDIM-1;
|
|
||||||
}
|
|
||||||
nmax = NMAX+1;
|
|
||||||
|
|
||||||
printf("NTOT = %d, ntot = %d, NMAX = %d\n",NTOT,ntot,NMAX);
|
|
||||||
|
|
||||||
/* wave vector computation */
|
|
||||||
normAxes(dir,nDir);
|
|
||||||
printf("Derivation direction (normed) %f %f %f\n",nDir[0],nDir[1],nDir[2]);
|
|
||||||
waveVect = (double *) malloc(Ngrid*sizeof(double));
|
|
||||||
if (waveVect == NULL) {
|
|
||||||
printf("waveVectCompute.cpp : No memory available for waveVect\n");
|
|
||||||
return;
|
|
||||||
}
|
|
||||||
waveVectorCompute3D((*grid).NX,(*grid).NY,(*grid).NZ,nDir,waveVect);
|
|
||||||
|
|
||||||
/* memory allocation */
|
|
||||||
table = (double *) malloc(ntot * sizeof(double));
|
|
||||||
if (table == NULL) {
|
|
||||||
printf("derivReal.cpp: No memory availble for table\n");
|
|
||||||
return;
|
|
||||||
}
|
|
||||||
|
|
||||||
realization = (double *) malloc(ntot * sizeof(double));
|
|
||||||
if (realization == NULL) {
|
|
||||||
printf("derivReal.cpp: No memory availble for realization\n");
|
|
||||||
return;
|
|
||||||
}
|
|
||||||
|
|
||||||
ExtendedWaveVect = (double *) malloc(ntot * sizeof(double));
|
|
||||||
if (ExtendedWaveVect == NULL) {
|
|
||||||
printf("derivReal.cpp: No memory availble for realization\n");
|
|
||||||
return;
|
|
||||||
}
|
|
||||||
|
|
||||||
workr = (double *) malloc(nmax * sizeof(double));
|
|
||||||
if (workr == NULL) {
|
|
||||||
printf("derivReal.cpp: No memory available for workr\n");
|
|
||||||
return;
|
|
||||||
}
|
|
||||||
|
|
||||||
worki = (double *) malloc(nmax * sizeof(double));
|
|
||||||
if (worki == NULL) {
|
|
||||||
printf("derivReal.cpp: No memory available for worki\n");
|
|
||||||
return;
|
|
||||||
}
|
|
||||||
|
|
||||||
extMaille =0;
|
|
||||||
/* organization of the extended realization */
|
|
||||||
for (k=1;k<=n[2];k++) {
|
|
||||||
for (j=1;j<=n[1];j++) {
|
|
||||||
for (i=1;i<=n[0];i++) {
|
|
||||||
extMaille = i + ((j-1) + (((k-1) * n[1]) * n[0]));
|
|
||||||
if (i <= (*grid).NX && j <= (*grid).NY && k <= (*grid).NZ) {
|
|
||||||
maille = i-1 + ((j-1) + ((k-1) * (*grid).NY) * (*grid).NX);
|
|
||||||
realization[extMaille] = (*Z).vector[maille];
|
|
||||||
ExtendedWaveVect[extMaille] = waveVect[maille];
|
|
||||||
} else {
|
|
||||||
realization[extMaille] = 0.0;
|
|
||||||
ExtendedWaveVect[extMaille] = 0.0;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
/* forward fourier transform of the realization */
|
|
||||||
fourt(realization,table,n,NDIM,1,0,workr,worki);
|
|
||||||
FILE *wave;
|
|
||||||
wave = fopen("/home/irsrvshare1/R03/UPS_FLEX/waveVector.eas","w");
|
|
||||||
|
|
||||||
for (i=1;i<ntot;i++) {
|
|
||||||
tmp = realization[i];
|
|
||||||
realization[i] = - ExtendedWaveVect[i]*table[i];
|
|
||||||
table[i] = - ExtendedWaveVect[i]*tmp;
|
|
||||||
fprintf(wave,"%f\n", ExtendedWaveVect[i]);
|
|
||||||
}
|
|
||||||
fclose(wave);
|
|
||||||
free(waveVect);
|
|
||||||
|
|
||||||
/* backward fourier tranform */
|
|
||||||
fourt(realization,table,n,NDIM,0,1,workr,worki);
|
|
||||||
|
|
||||||
/* free the memory */
|
|
||||||
free(table);
|
|
||||||
free(workr);
|
|
||||||
free(worki);
|
|
||||||
|
|
||||||
printf("after second fourt \n");
|
|
||||||
|
|
||||||
(*dZ).n = (*Z).n;
|
|
||||||
(*dZ).code = (*Z).code;
|
|
||||||
(*dZ).vector = (double *) malloc(Ngrid*sizeof(double));
|
|
||||||
|
|
||||||
FILE *real;
|
|
||||||
FILE *derReal;
|
|
||||||
|
|
||||||
real = fopen("/home/irsrvshare1/R03/UPS_FLEX/realization.eas","w");
|
|
||||||
derReal = fopen("/home/irsrvshare1/R03/UPS_FLEX/derivative.eas","w");
|
|
||||||
|
|
||||||
printf("before saving file \n");
|
|
||||||
|
|
||||||
for(k=1;k<=(*grid).NZ;k++) {
|
|
||||||
for(j=1;j<=(*grid).NY;j++) {
|
|
||||||
for(i=1;i<=(*grid).NX;i++) {
|
|
||||||
extMaille = i+(j-1+(k-1)*n[1])*n[0];
|
|
||||||
maille = i-1+(j-1+(k-1)*(*grid).NY)*(*grid).NX;
|
|
||||||
(*dZ).vector[maille] = realization[extMaille]/(double)NTOT;
|
|
||||||
// printf("Z[%d] = %f, dZ[%d] = %f\n",maille,(*Z).vector[maille],maille,(*dZ).vector[maille]);
|
|
||||||
fprintf(real,"%f\n",(*Z).vector[maille]);
|
|
||||||
fprintf(derReal,"%f\n",(*dZ).vector[maille]);
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
printf("after saving file \n");
|
|
||||||
|
|
||||||
fclose(real);
|
|
||||||
fclose(derReal);
|
|
||||||
free(realization);
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
void normAxes(double *vec, double *normed) {
|
|
||||||
|
|
||||||
int i;
|
|
||||||
double r;
|
|
||||||
|
|
||||||
r = sqrt(vec[0]*vec[0]+vec[1]*vec[1]+vec[2]*vec[2]);
|
|
||||||
|
|
||||||
if (normed == NULL) {
|
|
||||||
normed = (double *)malloc(3*sizeof(double));
|
|
||||||
if (normed == NULL) {
|
|
||||||
printf("derivReal.c: No memory available for normed\n");
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
for (i=0;i<3;i++) {
|
|
||||||
normed[i] = vec[i]/r;
|
|
||||||
}
|
|
||||||
|
|
||||||
}
|
|
||||||
|
|
@ -1,185 +0,0 @@
|
|||||||
|
|
||||||
#include <stdlib.h>
|
|
||||||
#include <math.h>
|
|
||||||
#include "geostat.h"
|
|
||||||
|
|
||||||
|
|
||||||
/*FAST FOURIER TRANSFORM MOVING AVERAGE METHOD */
|
|
||||||
/*Turns a Gaussian white noise vector into a */
|
|
||||||
/*spatially correlated vector */
|
|
||||||
/*input: */
|
|
||||||
/*variogram: structure defining the variogram */
|
|
||||||
/* model */
|
|
||||||
/*grid: structure defining the grid */
|
|
||||||
/*n: vector with the number of cells along the */
|
|
||||||
/* X, Y and Z axes for the underlying grid */
|
|
||||||
/* i = [0 1 2] */
|
|
||||||
/* --> 0 0 0 : n will be computed and */
|
|
||||||
/* updated as output */
|
|
||||||
/* --> nx ny nz: these dimensions are used */
|
|
||||||
/*realin: structure defining a realization - */
|
|
||||||
/* must be a Gaussian white noise */
|
|
||||||
/*output: */
|
|
||||||
/*realout: structure defining a realization - */
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
void FFTMA(struct vario_mod variogram,struct grid_mod grid,int n[3],struct realization_mod *realin,struct realization_mod *realout)
|
|
||||||
|
|
||||||
{
|
|
||||||
int NTOT,i,j,k,NMAX,NDIM,ntot,nmax,NXYZ,nxyz,maille0,maille1;
|
|
||||||
double *table,*covar,*workr,*worki,*realization,temp;
|
|
||||||
|
|
||||||
|
|
||||||
/*test over the input realization*/
|
|
||||||
/*if ((*realin).code != 0) {
|
|
||||||
printf("Input realizations in FFTMA must be Gaussian white noises");
|
|
||||||
exit;
|
|
||||||
}*/
|
|
||||||
|
|
||||||
/*covariance axis normalization*/
|
|
||||||
axes(variogram.ap,variogram.scf,variogram.Nvario);
|
|
||||||
|
|
||||||
|
|
||||||
/*pseudo-grid definition*/
|
|
||||||
cgrid(variogram,grid,n);
|
|
||||||
|
|
||||||
|
|
||||||
/*constant definition*/
|
|
||||||
NTOT = n[0]*n[1]*n[2];
|
|
||||||
ntot = NTOT+1;
|
|
||||||
NMAX = n[0];
|
|
||||||
NDIM = 3;
|
|
||||||
for (i=1;i<3;i++) {
|
|
||||||
if (n[i] > NMAX)
|
|
||||||
NMAX = n[i];
|
|
||||||
if (n[i] == 1)
|
|
||||||
NDIM--;
|
|
||||||
}
|
|
||||||
nmax = NMAX+1;
|
|
||||||
NXYZ = grid.NX*grid.NY*grid.NZ;
|
|
||||||
nxyz = NXYZ+1;
|
|
||||||
|
|
||||||
|
|
||||||
/*array initialization*/
|
|
||||||
covar = (double *) malloc(ntot * sizeof(double));
|
|
||||||
if (covar == NULL) {
|
|
||||||
printf("FFTMA.c: No memory available for covar");
|
|
||||||
exit;
|
|
||||||
}
|
|
||||||
|
|
||||||
table = (double *) malloc(ntot * sizeof(double));
|
|
||||||
if (table == NULL) {
|
|
||||||
printf("FFTMA.c: No memory available for table");
|
|
||||||
exit;
|
|
||||||
}
|
|
||||||
|
|
||||||
realization = (double *) malloc(ntot * sizeof(double));
|
|
||||||
if (realization == NULL) {
|
|
||||||
printf("FFTMA.c: No memory available for realization");
|
|
||||||
exit;
|
|
||||||
}
|
|
||||||
|
|
||||||
workr = (double *) malloc(nmax * sizeof(double));
|
|
||||||
if (workr == NULL) {
|
|
||||||
printf("FFTMA.c: No memory available for workr");
|
|
||||||
exit;
|
|
||||||
}
|
|
||||||
|
|
||||||
worki = (double *) malloc(nmax * sizeof(double));
|
|
||||||
if (worki == NULL) {
|
|
||||||
printf("FFTMA.c: No memory available for worki");
|
|
||||||
exit;
|
|
||||||
}
|
|
||||||
|
|
||||||
/*covariance function creation*/
|
|
||||||
covariance(covar,variogram,grid,n);
|
|
||||||
|
|
||||||
|
|
||||||
/*power spectrum*/
|
|
||||||
fourt(covar,table,n,NDIM,1,0,workr,worki);
|
|
||||||
|
|
||||||
|
|
||||||
/*organization of the input Gaussian white noise*/
|
|
||||||
for ( k = 1; k <= n[2]; k++) {
|
|
||||||
for (j = 1; j <= n[1]; j++) {
|
|
||||||
for (i = 1; i <= n[0]; i++) {
|
|
||||||
maille1 = i+(j-1+(k-1)*n[1])*n[0];
|
|
||||||
if (i <= grid.NX && j <= grid.NY && k <= grid.NZ) {
|
|
||||||
maille0 = i-1+(j-1+(k-1)*grid.NY)*grid.NX;
|
|
||||||
realization[maille1] = (*realin).vector[maille0];
|
|
||||||
} else {
|
|
||||||
realization[maille1] = 0.;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
/*forward fourier transform of the GWN*/
|
|
||||||
fourt(realization,table,n,NDIM,1,0,workr,worki);
|
|
||||||
|
|
||||||
|
|
||||||
/*decomposition and multiplication in the spectral domain*/
|
|
||||||
for ( k = 1; k <= n[2]; k++) {
|
|
||||||
for (j = 1; j <= n[1]; j++) {
|
|
||||||
for (i = 1; i <= n[0]; i++) {
|
|
||||||
maille1 = i+(j-1+(k-1)*n[1])*n[0];
|
|
||||||
temp = covar[maille1];
|
|
||||||
if (temp > 0.) {
|
|
||||||
temp = sqrt(temp)/(double) NTOT;
|
|
||||||
} else if (temp < 0.) {
|
|
||||||
temp = sqrt(-temp)/(double) NTOT;
|
|
||||||
}
|
|
||||||
realization[maille1] *= temp;
|
|
||||||
table[maille1] *= temp;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
free(covar);
|
|
||||||
|
|
||||||
|
|
||||||
/*backward fourier transform*/
|
|
||||||
fourt(realization,table,n,NDIM,0,1,workr,worki);
|
|
||||||
|
|
||||||
|
|
||||||
free(table);
|
|
||||||
free(workr);
|
|
||||||
free(worki);
|
|
||||||
|
|
||||||
|
|
||||||
/*output realization*/
|
|
||||||
|
|
||||||
/*is the output realization already allocated?*/
|
|
||||||
/*if not, memory allocation*/
|
|
||||||
if ((*realout).vector == NULL || (*realout).n != (*realin).n) {
|
|
||||||
(*realout).vector = (double *) malloc((*realin).n * sizeof(double));
|
|
||||||
if ((*realout).vector == NULL) {
|
|
||||||
printf("FFTMA.c: No memory available");
|
|
||||||
exit;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
(*realout).n = (*realin).n;
|
|
||||||
(*realout).code = 1;
|
|
||||||
for ( k = 1; k <= grid.NZ; k++) {
|
|
||||||
for (j = 1; j <= grid.NY; j++) {
|
|
||||||
for (i = 1; i <= grid.NX; i++) {
|
|
||||||
maille1 = i+(j-1+(k-1)*n[1])*n[0];
|
|
||||||
maille0 = i-1+(j-1+(k-1)*grid.NY)*grid.NX;
|
|
||||||
(*realout).vector[maille0] = realization[maille1];
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
free(realization);
|
|
||||||
|
|
||||||
|
|
||||||
return;
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
@ -1,36 +0,0 @@
|
|||||||
#
|
|
||||||
########################################################################
|
|
||||||
#
|
|
||||||
# Makefile for library
|
|
||||||
#
|
|
||||||
########################################################################
|
|
||||||
#
|
|
||||||
|
|
||||||
INCLUDE = -I./
|
|
||||||
LIBS = -lm -L./
|
|
||||||
# CCFLAGS = -fast
|
|
||||||
ifeq (${ARCH},sun)
|
|
||||||
CC = cc
|
|
||||||
CCFLAGS = -g
|
|
||||||
endif
|
|
||||||
ifeq (${ARCH},linux)
|
|
||||||
CC = gcc
|
|
||||||
CCFLAGS = -g
|
|
||||||
endif
|
|
||||||
|
|
||||||
LIB = libFFTMA_debug_${ARCH}.a
|
|
||||||
|
|
||||||
NOBJS= gammf.o FFTMA.o addstat.o axes.o cgrid.o covariance.o fourt.o length.o maxfactor.o test_fact.o cov_value.o generate.o gasdev.o ran2.o stable.o gaussian.o power.o cubic.o spherical.o nugget.o exponential.o cardsin.o nor2log.o
|
|
||||||
|
|
||||||
|
|
||||||
.c.o:
|
|
||||||
${CC} $(INCLUDE) $(CCFLAGS) -c $<
|
|
||||||
|
|
||||||
# LIBRARY
|
|
||||||
$(LIB) : $(NOBJS)
|
|
||||||
ar cr $(LIB) $(NOBJS)
|
|
||||||
ranlib $(LIB)
|
|
||||||
|
|
||||||
clean :
|
|
||||||
rm *.o
|
|
||||||
|
|
Loading…
Reference in New Issue