Commit 1b1d9ed0 authored by niklas.thiel's avatar niklas.thiel

small parallelisation

parent 958d3001
......@@ -1427,11 +1427,16 @@ int main(int argc, char **argv){
catseis(sectionp, fulldata_p, recswitch, ntr_glob, MPI_COMM_WORLD);
if (MYID==0) saveseis_glob(FP,fulldata_vx,fulldata_vy,fulldata_vz,fulldata_p,fulldata_curl,fulldata_div,recpos,recpos_loc,ntr_glob,srcpos,ishot,ns,iter,1);
/* wavefield separation*/
/* wavefield separation for forward modelling of STF*/
if (WAVESEP==1) {
if (MYID == 0) fprintf(FP,"\n start wavefield separation...");
if (MYID == 0) pup(fulldata_p, fulldata_vy, FP, ntr_glob,recpos,recpos_loc,srcpos,ishot,ns,iter,1);
if ((VERBOSE==1)&&(MYID == 0)) fprintf(FP,"\n wavefield separation finished");
if (MYID == 0 || MYID==1) pup(fulldata_p, fulldata_vy, ntr_glob,recpos,ishot,ns,iter);
if ((VERBOSE==1)&&(MYID == 0)) fprintf(FP,"\n wavefield separation finished"); /* save pup*/
if (VERBOSE==1 && MYID==0) {
fprintf(FP,"Write PUP data to file...\n ");
sprintf(pf,"%s_pup.su.stf.shot%i",SEIS_FILE,ishot);
outseis_glob(FP,fopen(pf,"w"), 0, fulldata_p,recpos,recpos_loc,ntr_glob,srcpos,1,ns,SEIS_FORMAT,ishot,0);
}
MPI_Barrier(MPI_COMM_WORLD);
MPI_Bcast( &fulldata_p[1][1], ntr_glob*ns, MPI_FLOAT, 0, MPI_COMM_WORLD);
h=1;
......@@ -1971,10 +1976,16 @@ int main(int argc, char **argv){
}
catseis(sectionp, fulldata_p, recswitch, ntr_glob, MPI_COMM_WORLD);
/* wavefield separation */
/* wavefield separation*/
if (WAVESEP==1) {
if (MYID == 0) fprintf(FP,"\n start wavefield separation...");
if (MYID == 0) pup(fulldata_p, fulldata_vy, FP, ntr_glob,recpos,recpos_loc,srcpos,ishot,ns,iter,1);
if (MYID == 0 || MYID==1) pup(fulldata_p, fulldata_vy, ntr_glob,recpos,ishot,ns,iter);
if ((VERBOSE==1)&&(MYID == 0)) fprintf(FP,"\n wavefield separation finished"); /* save pup*/
if (VERBOSE==1 && MYID==0) {
fprintf(FP,"Write PUP data to file...\n ");
sprintf(pf,"%s_pup.su.shot%i",SEIS_FILE,ishot);
outseis_glob(FP,fopen(pf,"w"), 0, fulldata_p,recpos,recpos_loc,ntr_glob,srcpos,1,ns,SEIS_FORMAT,ishot,0);
}
if ((VERBOSE==1) && (MYID == 0)) fprintf(FP,"wavefield separation finished.");
MPI_Barrier(MPI_COMM_WORLD);
MPI_Bcast( &fulldata_p[1][1], ntr_glob*ns, MPI_FLOAT, 0, MPI_COMM_WORLD);
......@@ -3519,7 +3530,13 @@ int main(int argc, char **argv){
catseis(sectionvy, fulldata_vy, recswitch, ntr_glob, MPI_COMM_WORLD);
catseis(sectionp, fulldata_p, recswitch, ntr_glob, MPI_COMM_WORLD);
if (MYID == 0) fprintf(FP,"\n start wavefield separation...");
if (MYID == 0) pup(fulldata_p, fulldata_vy, FP, ntr_glob,recpos,recpos_loc,srcpos,ishot,ns,iter,2);
if (MYID == 0 || MYID==1) pup(fulldata_p, fulldata_vy, ntr_glob,recpos,ishot,ns,iter);
if ((VERBOSE==1)&&(MYID == 0)) fprintf(FP,"\n wavefield separation finished"); /* save pup*/
if (VERBOSE==1 && MYID==0) {
fprintf(FP,"Write PUP data to file...\n ");
sprintf(pf,"%s_pup.su.step.shot%i",SEIS_FILE,ishot);
outseis_glob(FP,fopen(pf,"w"), 0, fulldata_p,recpos,recpos_loc,ntr_glob,srcpos,1,ns,SEIS_FORMAT,ishot,0);
}
if ((VERBOSE==1) && (MYID == 0)) fprintf(FP,"\n wavefield separation finished");
MPI_Barrier(MPI_COMM_WORLD);
MPI_Bcast( &fulldata_p[1][1], ntr_glob*ns, MPI_FLOAT, 0, MPI_COMM_WORLD);
......
......@@ -228,7 +228,7 @@ void psource(int nt, float ** sxx, float ** syy, float ** sp,
void psource_rsg(int nt, float ** sxx, float ** syy,
float ** srcpos_loc, float ** signals, int nsrc);
void pup(float **data_p, float **data_vz, FILE *fp, int ntr_glob, int **recpos, int **recpos_loc, float ** srcpos, int ishot, int ns, int iter, int sw );
void pup(float** data_p, float** data_vy, int ntr_glob, int** recpos, int ishot, int ns, int iter );
float *rd_sour(int *nts,FILE* fp_source);
......
......@@ -24,39 +24,29 @@
void fft2(float **array, float **arrayim, int NYG, int NXG, int dir) {
extern FILE *FP;
COMPLEX *C;
int j, i, k, nx, ny, nxy;
float **RE, **IM;
char pf1[STRING_SIZE];
FILE *file;
extern char SEIS_FILE[STRING_SIZE];
nx = NXG;
ny = NYG;
int j, i, k;
C = cplxvector(1, ny*nx);
C = cplxvector(1, NYG*NXG);
if (dir==1) {
k = 1;
for (i=1; i<=nx; i++)
for (j=1; j<=ny; j++) {
for (i=1; i<=NXG; i++)
for (j=1; j<=NYG; j++) {
C[k].re = array[j][i];
C[k].im = 0.0;
k++;
}
forward_fft2f(C, ny, nx);
forward_fft2f(C, NXG, NYG);
k = 1;
for (i=1; i<=nx; i++)
for (j=1; j<=ny; j++) {
for (i=1; i<=NXG; i++)
for (j=1; j<=NYG; j++) {
array[j][i] = C[k].re;
arrayim[j][i] = C[k].im;
k++;
......@@ -66,24 +56,24 @@ void fft2(float **array, float **arrayim, int NYG, int NXG, int dir) {
k = 1;
for (i=1; i<=nx; i++)
for (j=1; j<=ny; j++) {
for (i=1; i<=NXG; i++)
for (j=1; j<=NYG; j++) {
C[k].re = array[j][i];
C[k].im = arrayim[j][i];
k++;
}
inverse_fft2f(C, ny, nx);
inverse_fft2f(C, NYG, NXG);
k = 1;
for (i=1; i<=nx; i++)
for (j=1; j<=ny; j++) {
for (i=1; i<=NXG; i++)
for (j=1; j<=NYG; j++) {
array[j][i] = C[k].re;
arrayim[j][i] = 0.0;
k++;
}
}
free_cplxvector(C, 1, ny*nx);
free_cplxvector(C, 1, NYG*NXG);
}
......@@ -27,38 +27,41 @@
#include "fd.h"
void pup(float **data_p, float **data_vy, FILE *fp, int ntr_glob, int **recpos, int **recpos_loc, float **srcpos, int ishot, int ns, int iter, int sw) {
void pup(float **data_p, float **data_vy, int ntr_glob, int **recpos, int ishot, int ns, int iter) {
extern float DT,DH;
extern float VEL, DENS, ANGTAPI, ANGTAPO;
extern char SEIS_FILE[STRING_SIZE];
extern int SEIS_FORMAT, VERBOSE;
extern int SEIS_FORMAT, VERBOSE,MYID;
extern FILE *FP;
int i, j, k, nx_pot, ny_pot;
float nyq_k, nyq_f, dk, df, angi, ango, ka, dtr;
int i, j, k, nx_pot, ny_pot, dtr;
double t1,t2,t3,t5,t6;
float nyq_k, nyq_f, dk, df, angi, ango, ka;
float *w, *kx;
float **filtwk, **filtwk_full;
float **data_pim,**data_vyim,**data_pup,**data_pupim,**data_pt,**data_vyt;
char pf[STRING_SIZE], pf1[STRING_SIZE], pf2[STRING_SIZE];
FILE *file;
MPI_Status status;
/* calculate receiver sampling */
dtr=(recpos[1][ntr_glob]-recpos[1][1])/(ntr_glob-1)*DH;
if (VERBOSE==1) { fprintf(FP,"\n Receiver sampling was calculatet to dtr: %f m \n ",dtr); }
/* double the dimensions to avoid artefacts and recompute the dimensions to become values of 2^m */
nx_pot = (int)(pow(2,(ceil(log(ntr_glob*2)/log(2)))));
ny_pot = (int)(pow(2,(ceil(log(ns*2)/log(2)))));
data_pim=matrix(1,ny_pot,1,nx_pot);
data_pt=matrix(1,ny_pot,1,nx_pot);
if (MYID==0){
t1=MPI_Wtime();
/* calculate receiver sampling */
dtr=(recpos[1][ntr_glob]-recpos[1][1])/(ntr_glob-1)*DH;
if (VERBOSE==1) { fprintf(FP,"\n Receiver sampling was calculatet to dtr: %i m \n ",dtr); }
/* temporary array */
filtwk=matrix(1,ny_pot/2,1,nx_pot/2);
filtwk_full=matrix(1,ny_pot,1,nx_pot);
data_pim=matrix(1,ny_pot,1,nx_pot);
data_vyim=matrix(1,ny_pot,1,nx_pot);
data_pup=matrix(1,ny_pot,1,nx_pot);
data_pupim=matrix(1,ny_pot,1,nx_pot);
data_pt=matrix(1,ny_pot,1,nx_pot);
data_vyt=matrix(1,ny_pot,1,nx_pot);
w=vector(1,ny_pot/2);
......@@ -113,7 +116,7 @@ void pup(float **data_p, float **data_vy, FILE *fp, int ntr_glob, int **recpos,
}
/* write filter to file*/
if (iter==1 && ishot==1) {
if (VERBOSE==1 && iter==1 && ishot==1) {
sprintf(pf1,"%s_spectrum.filter.bin",SEIS_FILE);
file = fopen(pf1,"wb");
......@@ -122,33 +125,46 @@ void pup(float **data_p, float **data_vy, FILE *fp, int ntr_glob, int **recpos,
fclose(file);
}
if (VERBOSE==1) { fprintf(FP,"2D FFT of p-data and vz-data... \n"); }
t2=MPI_Wtime();
fprintf(FP,"\ntime for calculating filter: %f\n",t2-t1);
}
/* copy data to temp array */
if (MYID==1){ /* copy data to temp array */
for (k=1; k<=ntr_glob; k++) {
for (i=1; i<=ns; i++) {
data_pt[i][k]=data_p[k][i];
data_vyt[i][k]=data_vy[k][i];
}
}
if (VERBOSE==1) { fprintf(FP,"2D FFT of p-data and vz-data... \n"); }
/* transformation into fk domain */
fft2(data_pt, data_pim, ny_pot, nx_pot,1);
MPI_Send(&data_pt[1][1],ny_pot*nx_pot , MPI_FLOAT, 0, 1, MPI_COMM_WORLD);
}
if (MYID==0){
/* copy data to temp array */
for (k=1; k<=ntr_glob; k++) {
for (i=1; i<=ns; i++) {
data_vyt[i][k]=data_vy[k][i];
}
}
fft2(data_vyt, data_vyim, ny_pot, nx_pot,1);
MPI_Recv(&data_pt[1][1], ny_pot*nx_pot, MPI_FLOAT, 1, 1, MPI_COMM_WORLD, &status);
/* write spectrum to file */
if (ishot==1 && iter==1) {
sprintf(pf2,"%s_spectrum.p.real.shot%i.it%i.bin",SEIS_FILE,ishot, iter);
file = fopen(pf2,"wb");
for (i=1; i<=nx_pot; i++)
for (j=1; j<=ny_pot; j++) { fwrite(&data_pt[j][i],sizeof(float),1,file); }
fclose(file);
sprintf(pf2,"%s_spectrum.p.real.shot%i.it%i.bin",SEIS_FILE,ishot, iter);
file = fopen(pf2,"wb");
for (i=1; i<=nx_pot; i++)
for (j=1; j<=ny_pot; j++) { fwrite(&data_pt[j][i],sizeof(float),1,file); }
fclose(file);
}
fft2(data_vyt, data_vyim, ny_pot, nx_pot,1);
t3=MPI_Wtime();
fprintf(FP,"\ntime for one 2d FFT: %f s\n",t3-t2);
/* write spectrum to file */
if (ishot==1 && iter==1) {
sprintf(pf2,"%s_spectrum.vz.real.shot%i.it%i.bin",SEIS_FILE,ishot, iter);
......@@ -169,7 +185,8 @@ void pup(float **data_p, float **data_vy, FILE *fp, int ntr_glob, int **recpos,
data_pupim[i][k]=0.5*(data_pim[i][k]-(filtwk_full[i][k]*data_vyim[i][k]));
}
}
t5=MPI_Wtime();
fprintf(FP,"\ntime for wavefield separation: %f s\n",t5-t3);
/* write spectrum to file */
if (ishot==1 && iter==1) {
sprintf(pf2,"%s_spectrum.pup.real.shot%i.it%i.bin",SEIS_FILE,ishot, iter);
......@@ -197,16 +214,19 @@ void pup(float **data_p, float **data_vy, FILE *fp, int ntr_glob, int **recpos,
data_p[k][1]=0; /*first sample need to be set to 0 to avoid artefacts after integrating from velocity to displacement*/
}
/* save pup*/
if (VERBOSE==1) {
fprintf(FP,"Write PUP data to file...\n ");
if (sw == 1) sprintf(pf,"%s_pup.su.shot%i",SEIS_FILE,ishot);
if (sw == 2) sprintf(pf,"%s_pup.su.step.shot%i",SEIS_FILE,ishot);
// /* save pup*/
// if (VERBOSE==1) {
// fprintf(FP,"Write PUP data to file...\n ");
//
// if (sw == 1) sprintf(pf,"%s_pup.su.shot%i",SEIS_FILE,ishot);
// if (sw == 2) sprintf(pf,"%s_pup.su.step.shot%i",SEIS_FILE,ishot);
//
// outseis_glob(fp,fopen(pf,"w"), 0, data_p,recpos,recpos_loc,ntr_glob,srcpos,1,ns,SEIS_FORMAT,ishot,0);
// }
t6=MPI_Wtime();
fprintf(FP,"\ntime for reverse fft and writing: %f s\n",t6-t5);
fprintf(FP,"\ntime for pup.c: %f s\n",t6-t1);
outseis_glob(fp,fopen(pf,"w"), 0, data_p,recpos,recpos_loc,ntr_glob,srcpos,1,ns,SEIS_FORMAT,ishot,0);
}
free_matrix(filtwk, 1,ny_pot/2,1,nx_pot/2);
free_matrix(filtwk_full, 1,ny_pot,1,nx_pot);
free_matrix(data_pim, 1,ny_pot,1,nx_pot);
......@@ -217,6 +237,7 @@ void pup(float **data_p, float **data_vy, FILE *fp, int ntr_glob, int **recpos,
free_matrix(data_vyt, 1,ny_pot,1,nx_pot);
free_vector(w,1,ny_pot/2);
free_vector(kx,1,ny_pot/2);
}
}
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment