Commit b9056d1e authored by niklas.thiel's avatar niklas.thiel

parallised wavefield separation

parent d6833638
......@@ -816,6 +816,17 @@ int main(int argc, char **argv){
MPI_Barrier(MPI_COMM_WORLD);
/* MPI split for processors used for pup parallelisation */
int myid_pup, pupid=0;
MPI_Comm MPI_COMM_PUP;
if (MYID<=7) pupid = 1;
else pupid = 0;
MPI_Comm_split(MPI_COMM_WORLD, pupid, MYID, &MPI_COMM_PUP);
MPI_Comm_rank(MPI_COMM_PUP, &myid_pup);
/* end of MPI split for processors with ntr>0 */
SHOTINC=1;
RECINC=1;
......@@ -1384,7 +1395,7 @@ int main(int argc, char **argv){
/* wavefield separation for forward modelling of STF*/
if (WAVESEP==1) {
if (MYID == 0) fprintf(FP,"\n start wavefield separation...");
if (MYID == 0 || MYID==1) pup(fulldata_p, fulldata_vy, ntr_glob,recpos,ishot,ns,iter);
if (pupid) pup(fulldata_p, fulldata_vy, ntr_glob,recpos,ishot,ns,iter,MPI_COMM_PUP,myid_pup);
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 ");
......@@ -1939,7 +1950,7 @@ int main(int argc, char **argv){
/* wavefield separation*/
if (WAVESEP==1) {
if (MYID == 0) fprintf(FP,"\n start wavefield separation...");
if (MYID == 0 || MYID==1) pup(fulldata_p, fulldata_vy, ntr_glob,recpos,ishot,ns,iter);
if (pupid) pup(fulldata_p, fulldata_vy, ntr_glob,recpos,ishot,ns,iter,MPI_COMM_PUP,myid_pup);
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 ");
......@@ -3542,7 +3553,7 @@ 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 || MYID==1) pup(fulldata_p, fulldata_vy, ntr_glob,recpos,ishot,ns,iter);
if (pupid) pup(fulldata_p, fulldata_vy, ntr_glob,recpos,ishot,ns,iter,MPI_COMM_PUP,myid_pup);
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 ");
......
......@@ -19,7 +19,7 @@ EXEC= ../bin
# LINUX with OpenMPI / IntelMPI and INTEL Compiler
# Use icc whenever possible, this will be much faster than gcc
CC=mpicc
CC=mpiicc
LFLAGS=-lm -lcseife -lstfinv -laff -lfourierxx -lfftw3 -lstdc++
CFLAGS=-O3
SFLAGS=-L./../contrib/libcseife -L./../contrib/bin
......@@ -183,8 +183,7 @@ IFOS2D= \
create_trkill_table.c \
filter_frequencies.c \
pup.c \
fft2_filt.c \
fft2d.c
fft.c
# -------------
# Targes
......
......@@ -180,8 +180,10 @@ double calc_misfit(float **sectiondata, float **section, int ntr, int ns, int LN
abs_sectiondata+=intseis_sectiondata[i][j]*intseis_sectiondata[i][j];
abs_section+=intseis_section[i][j]*intseis_section[i][j];
}
abs_sectiondata=sqrt(abs_sectiondata);
abs_section=sqrt(abs_section);
if (abs_sectiondata==0) abs_sectiondata=1;
else abs_sectiondata=sqrt(abs_sectiondata);
if (abs_section==0) abs_section==1;
else abs_section=sqrt(abs_section);
}
/* calculate L2 residuals */
......
......@@ -207,8 +207,10 @@ double calc_res(float **sectiondata, float **section, float **sectiondiff, float
abs_section+=intseis_section[i][j]*intseis_section[i][j];
sectiondata_mult_section+=intseis_sectiondata[i][j]*intseis_section[i][j]; /* calculation of dot product for measured (section) and synthetic (sectiondata) data*/
}
abs_sectiondata=sqrt(abs_sectiondata);
abs_section=sqrt(abs_section);
if (abs_sectiondata==0) abs_sectiondata=1;
else abs_sectiondata=sqrt(abs_sectiondata);
if (abs_section==0) abs_section==1;
else abs_section=sqrt(abs_section);
}
/* calculate residual seismograms and norm */
......
......@@ -25,29 +25,6 @@
#define REQUEST_COUNT 4
#define WORKFLOW_MAX_VAR 13
/***************/
/* FFT2D stuff */
/***************/
typedef struct {float re; float im;} COMPLEX;
typedef struct {double re; double im;} DCOMPLEX;
#ifndef ERROR
#define ERROR -1
#define NO_ERROR 0
#endif
#define FFT_FORWARD 0
#define FFT_INVERSE 1
#define TWOPI 6.2831853071795865 /* 2.0 * PI */
#define HALFPI 1.5707963267948966 /* PI / 2.0 */
#define PI8 0.392699081698724 /* PI / 8.0 */
#define RT2 1.4142135623731 /* sqrt(2.0) */
#define IRT2 0.707106781186548 /* 1.0/sqrt(2.0) */
extern int forward_fft2f(COMPLEX *array, int rows, int cols); /* Perform forward 2D transform on a COMPLEX array. */
extern int inverse_fft2f(COMPLEX *array, int rows, int cols); /* Perform inverse 2D transform on a COMPLEX array. */
/****************************/
/* declaration of functions */
/****************************/
......@@ -156,9 +133,9 @@ void FFT_filt(float ** data, float freqshift, int ntr, int ns,int method);
/*short FFT(short int dir,long m,float *x,float *y);*/
void FFT(int isign, unsigned long nlog, float *re, float *im); /* NR version*/
// void FFT(int isign, unsigned long nlog, float *re, float *im); /* NR version*/
void fft2(float **array, float **arrayim, int NYG, int NXG, int dir); /* 2d fft filtering*/
void fft(float *array, float *arrayim, int npad, int dir); /* 1d fft*/
float *filter_frequencies(int *nfrq);
......@@ -236,7 +213,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_vy, int ntr_glob, int** recpos, int ishot, int ns, int iter );
void pup(float** data_p, float** data_vy, int ntr_glob, int** recpos, int ishot, int ns, int iter, MPI_Comm newcomm_nodentr, int myid_pup );
float *rd_sour(int *nts,FILE* fp_source);
......@@ -500,7 +477,6 @@ void declare_error(char err_text[]);
void warning(char warn_text[]);
double maximum(float **a, int nx, int ny);
float *vector(int nl, int nh);
COMPLEX *cplxvector(int nstart, int nend);
int *ivector(int nl, int nh);
double *dvector(int nl, int nh);
float **fmatrix(int nrl, int nrh, int ncl, int nch);
......@@ -510,7 +486,6 @@ float **matrix(int nrl, int nrh, int ncl, int nch);
int **imatrix(int nrl, int nrh, int ncl, int nch);
float ***f3tensor(int nrl, int nrh, int ncl, int nch,int ndl, int ndh);
void free_vector(float *v, int nl, int nh);
void free_cplxvector(COMPLEX *v, int nstart, int nend);
void free_dvector(double *v, int nl, int nh);
void free_ivector(int *v, int nl, int nh);
void free_matrix(float **m, int nrl, int nrh, int ncl, int nch);
......
/*-----------------------------------------------------------------------------------------
* Copyright (C) 2016 For the list of authors, see file AUTHORS.
*
* This file is part of IFOS.
*
* IFOS is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, version 2.0 of the License only.
*
* IFOS is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with IFOS. See file COPYING and/or <http://www.gnu.org/licenses/gpl-2.0.html>.
-----------------------------------------------------------------------------------------*/
/*------------------------------------------------------------------------
* 1D-FFT
* ----------------------------------------------------------------------*/
#include "fd.h"
#include <fftw3.h>
void fft(float *array, float *arrayim, int npad, int dir) {
/* declaration of local variables */
int i,j;
/* declaration of variables for FFTW3*/
fftw_complex *in,*out;
fftw_plan p;
in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * npad);
out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * npad);
for (j=0;j<npad;j++){
in[j][0]=array[j+1];
in[j][1]=arrayim[j+1];
}
if (dir==1) { p = fftw_plan_dft_1d(npad, in, out, FFTW_FORWARD, FFTW_ESTIMATE); }
else { p = fftw_plan_dft_1d(npad, in, out, FFTW_BACKWARD, FFTW_ESTIMATE); }
fftw_execute(p);
/* write output into data matrix */
for(j=0;j<npad;j++){
array[j+1] = out[j][0];
arrayim[j+1] = out[j][1];
}
fftw_free(in);
fftw_free(out);
fftw_destroy_plan(p);
}
......@@ -27,207 +27,460 @@
#include "fd.h"
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,MYID;
extern FILE *FP;
void pup(float **data_p, float **data_vy, int ntr_glob, int **recpos, int ishot, int ns, int iter, MPI_Comm newcomm_nodentr,int myid_pup) {
extern float DT,DH;
extern float VEL, DENS, ANGTAPI, ANGTAPO;
extern char SEIS_FILE[STRING_SIZE];
extern int SEIS_FORMAT, VERBOSE,MYID;
extern FILE *FP;
int i, j, k, nx_pot, ny_pot;
double t0,t1,t2,t3,t4,t5,t6,t7,t8;
int no_myid,ntr_start,ntr_end,ns_start,ns_end;
float *col,*col_im,*row,*row_im;
float nyq_k, nyq_f, dk, df, angi, ango, ka, dtr;
float *w, *kx;
float **filtwk, **filtwk_full;
float **data_pim,**data_im1,**data_vyim,**data_pup,**data_pupim,**data_pt,**data_t1,**data_vyt,**data_vyt1;
char pf[STRING_SIZE], pf1[STRING_SIZE], pf2[STRING_SIZE];
FILE *file;
MPI_Status status;
/* MPI split for processors used for pup parallelisation */
int myid_pup1, myid_pup2, pupid1=0, pupid2=0;
MPI_Comm MPI_COMM_PUP1,MPI_COMM_PUP2;
/* number of processors for every fft */
no_myid=4;
if (MYID<=3) pupid1 = 1;
else pupid2 = 1;
MPI_Comm_split(newcomm_nodentr, pupid1, MYID, &MPI_COMM_PUP1);
MPI_Comm_rank(MPI_COMM_PUP1, &myid_pup1);
MPI_Comm_split(newcomm_nodentr, pupid2, MYID, &MPI_COMM_PUP2);
MPI_Comm_rank(MPI_COMM_PUP2, &myid_pup2);
t1=MPI_Wtime();
/* 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)))));
fprintf(FP,"nx_pot: %i; ny_pot: %i",nx_pot,ny_pot);
data_pim=matrix(1,ny_pot,1,nx_pot);
data_pt=matrix(1,ny_pot,1,nx_pot);
data_im1=matrix(1,ny_pot,1,nx_pot);
data_t1=matrix(1,ny_pot,1,nx_pot);
col = vector(1,ny_pot);
col_im = vector(1,ny_pot);
row = vector(1,nx_pot);
row_im = vector(1,nx_pot);
data_pup=matrix(1,ny_pot,1,nx_pot);
data_pupim=matrix(1,ny_pot,1,nx_pot);
data_vyt=matrix(1,ny_pot,1,nx_pot);
data_vyim=matrix(1,ny_pot,1,nx_pot);
int i, j, k, nx_pot, ny_pot, dtr;
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;
/* 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){
/* 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_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_vyt=matrix(1,ny_pot,1,nx_pot);
w=vector(1,ny_pot/2);
kx=vector(1,nx_pot/2);
/* calculate variable for wave field separation */
nyq_k = 1/(2*dtr); /* Nyquist of data in first dimension */
nyq_f = 1/(2*DT); /* Nyquist of data in time dimension */
dk = 1/((nx_pot/2)*dtr); /* Wavenumber increment */
df = 1/((ny_pot/2)*DT); /* Frequency increment */
for (j=1; j<=ny_pot/2; j++) { w[j]=df*(j-1)*2*PI; } /* vector of angular frequencies*/
for (j=1; j<=nx_pot/2; j++) { kx[j]=dk*(j-1)*2*PI; } /* vector of angular wavenumbers*/
angi = ANGTAPI * PI/180.0; /* angle at which taper begin */
ango = ANGTAPO * PI/180.0; /* angle at which taper ends */
/* calculate Filter vor wavefield separation for first quadrant (after Kluever, 2008) */
if (VERBOSE==1) { fprintf(FP,"Calculate Filter for wavefield separation... \n "); }
for (k=1; k<=ny_pot/2; k++) {
ka = fabs(w[k])/VEL;
for (i=1; i<=nx_pot/2; i++) {
if (kx[i] < ka*sin(angi)) {
filtwk[k][i] = (DENS*fabs(w[k]))/sqrt((ka*ka)-(kx[i]*kx[i]));
} else {
if (kx[i] < ka*sin(ango) && kx[i] > ka*sin(angi)) {
filtwk[k][i] = 0.5*(1.0+cos(PI/((ka*sin(ango))-(ka*sin(angi))) * (kx[i]-(ka*sin(angi))))) * (DENS*fabs(w[k]))/sqrt((ka*ka)-(kx[i]*kx[i]));
} else {
filtwk[k][i] = 0.0;
if (MYID==0) {
/* 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);
w=vector(1,ny_pot/2);
kx=vector(1,nx_pot/2);
data_vyt1=matrix(1,ny_pot,1,nx_pot);
/* calculate variable for wave field separation */
nyq_k = 1/(2*dtr); /* Nyquist of data in first dimension */
nyq_f = 1/(2*DT); /* Nyquist of data in time dimension */
dk = 1/((nx_pot/2)*dtr); /* Wavenumber increment */
df = 1/((ny_pot/2)*DT); /* Frequency increment */
for (j=1; j<=ny_pot/2; j++) {
w[j]=df*(j-1)*2*PI; /* vector of angular frequencies*/
}
for (j=1; j<=nx_pot/2; j++) {
kx[j]=dk*(j-1)*2*PI; /* vector of angular wavenumbers*/
}
angi = ANGTAPI * PI/180.0; /* angle at which taper begin */
ango = ANGTAPO * PI/180.0; /* angle at which taper ends */
/* calculate Filter vor wavefield separation for first quadrant (after Kluever, 2008) */
if (VERBOSE==1) {
fprintf(FP,"Calculate Filter for wavefield separation... \n ");
}
for (k=1; k<=ny_pot/2; k++) {
ka = fabs(w[k])/VEL;
for (i=1; i<=nx_pot/2; i++) {
if (kx[i] < ka*sin(angi)) {
filtwk[k][i] = (DENS*fabs(w[k]))/sqrt((ka*ka)-(kx[i]*kx[i]));
} else {
if (kx[i] < ka*sin(ango) && kx[i] > ka*sin(angi)) {
filtwk[k][i] = 0.5*(1.0+cos(PI/((ka*sin(ango))-(ka*sin(angi))) * (kx[i]-(ka*sin(angi))))) * (DENS*fabs(w[k]))/sqrt((ka*ka)-(kx[i]*kx[i]));
} else {
filtwk[k][i] = 0.0;
}
}
}
}
/* write filter to file*/
if (VERBOSE==1 && iter==1 && ishot==1) {
sprintf(pf1,"%s_spectrum.filter_q.bin",SEIS_FILE);
file = fopen(pf1,"wb");
for (i=1; i<=nx_pot/2; i++)
for (j=1; j<=ny_pot/2; j++) {
fwrite(&filtwk[j][i],sizeof(float),1,file);
}
fclose(file);
}
/* expand filter to all 4 quadrants */
for (k=1; k<=ny_pot; k++) {
for (i=1; i<=nx_pot; i++) {
if (k<=ny_pot/2 && i<=nx_pot/2) {
filtwk_full[k][i]=filtwk[(ny_pot/2)-(k-1)][(nx_pot/2)-(i-1)];
} else if (k<=ny_pot/2 && i>nx_pot/2) {
filtwk_full[k][i]=filtwk[(ny_pot/2)-(k-1)][i-(nx_pot/2)];
} else if (k>ny_pot/2 && i<=nx_pot/2) {
filtwk_full[k][i]=filtwk[k-(ny_pot/2)][(nx_pot/2)-(i-1)];
} else if (k>ny_pot/2 && i>nx_pot/2) {
filtwk_full[k][i]=filtwk[k-(ny_pot/2)][i-(nx_pot/2)];
}
}
}
// for (k=1; k<=ny_pot; k++) {
// for (i=1; i<=nx_pot; i++) {
// if (k<=ny_pot/2 && i<=nx_pot/2) {
// filtwk_full[k][i]=filtwk[k][i];
// } else if (k<=ny_pot/2 && i>nx_pot/2) {
// filtwk_full[k][i]=filtwk[k][(nx_pot)-(i-1)];
// } else if (k>ny_pot/2 && i<=nx_pot/2) {
// filtwk_full[k][i]=filtwk[(ny_pot)-(k-1)][i];
// } else if (k>ny_pot/2 && i>nx_pot/2) {
// filtwk_full[k][i]=filtwk[(ny_pot)-(k-1)][(nx_pot)-(i-1)];
// }
//
// }
// }
t0=MPI_Wtime();
fprintf(FP,"\ntime for calculating filter: %f\n",t0-t1);
/* write filter to file*/
if (VERBOSE==1 && iter==1 && ishot==1) {
sprintf(pf1,"%s_spectrum.filter.bin",SEIS_FILE);
file = fopen(pf1,"wb");
for (i=1; i<=nx_pot; i++)
for (j=1; j<=ny_pot; j++) {
fwrite(&filtwk_full[j][i],sizeof(float),1,file);
}
fclose(file);
}
if (VERBOSE==1) {
fprintf(FP,"2D FFT of p-data and vz-data... \n");
}
}
// MPI_Barrier(newcomm_nodentr);
t2=MPI_Wtime();
if (pupid1) { /* copy data to temp array */
ntr_start=nx_pot/(no_myid)*(myid_pup1)+1;
ntr_end=nx_pot/(no_myid)*(myid_pup1+1);
ns_start=ny_pot/(no_myid)*(myid_pup1)+1;
ns_end=ny_pot/(no_myid)*(myid_pup1+1);
/* transformation into fk domain */
/* columns */
for (k=ntr_start; k<=ntr_end; k++) {
if (k<=ntr_glob) {
for (i=1; i<=ns; i++) {
col[i]=data_p[k][i];
};
fft(col, col_im, ny_pot,1);
for (i=1; i<=ny_pot; i++) {
data_t1[i][k]=col[i];
data_im1[i][k]=col_im[i];
col[i]=0;
col_im[i]=0;
}
}
}
MPI_Allreduce(MPI_IN_PLACE, &data_t1[1][1], nx_pot*ny_pot, MPI_FLOAT, MPI_SUM, MPI_COMM_PUP1);
MPI_Allreduce(MPI_IN_PLACE, &data_im1[1][1], nx_pot*ny_pot, MPI_FLOAT, MPI_SUM, MPI_COMM_PUP1);
/* rows */
for (i=ns_start; i<=ns_end; i++) {
for (k=1; k<=nx_pot; k++) {
row[k]=data_t1[i][k];
row_im[k]=data_im1[i][k];
}
fft(row, row_im, nx_pot,1);
for (k=1; k<=nx_pot; k++) {
// data_pt[i][k]=row[k];
// data_pim[i][k]=row_im[k];
if (i<=ny_pot/2 && k<=nx_pot/2) {
data_pt[i+ny_pot/2][k+nx_pot/2]=row[k];
data_pim[i+ny_pot/2][k+nx_pot/2]=row_im[k];
} else if (i<=ny_pot/2 && k>nx_pot/2) {
data_pt[i+ny_pot/2][k-nx_pot/2]=row[k];
data_pim[i+ny_pot/2][k-nx_pot/2]=row_im[k];
} else if (i>ny_pot/2 && k<=nx_pot/2) {
data_pt[i-ny_pot/2][k+nx_pot/2]=row[k];
data_pim[i-ny_pot/2][k+nx_pot/2]=row_im[k];
} else if (i>ny_pot/2 && k>nx_pot/2) {
data_pt[i-ny_pot/2][k-nx_pot/2]=row[k];
data_pim[i-ny_pot/2][k-nx_pot/2]=row_im[k];
}
}
}
}
MPI_Allreduce(MPI_IN_PLACE, &data_pt[1][1], nx_pot*ny_pot, MPI_FLOAT, MPI_SUM, newcomm_nodentr);
MPI_Allreduce(MPI_IN_PLACE, &data_pim[1][1], nx_pot*ny_pot, MPI_FLOAT, MPI_SUM, newcomm_nodentr);
if (pupid2) { /* copy data to temp array */
ntr_start=nx_pot/(no_myid)*(myid_pup2)+1;
ntr_end=nx_pot/(no_myid)*(myid_pup2+1);
ns_start=ny_pot/(no_myid)*(myid_pup2)+1;
ns_end=ny_pot/(no_myid)*(myid_pup2+1);
/* transformation into fk domain */
/* columns */
for (k=ntr_start; k<=ntr_end; k++) {
if (k<=ntr_glob) {
for (i=1; i<=ns; i++) {
col[i]=data_vy[k][i];
};
fft(col, col_im, ny_pot,1);
for (i=1; i<=ny_pot; i++) {
data_t1[i][k]=col[i];
data_im1[i][k]=col_im[i];
col[i]=0;
col_im[i]=0;
}
}
}
MPI_Allreduce(MPI_IN_PLACE, &data_t1[1][1], nx_pot*ny_pot, MPI_FLOAT, MPI_SUM, MPI_COMM_PUP2);
MPI_Allreduce(MPI_IN_PLACE, &data_im1[1][1], nx_pot*ny_pot, MPI_FLOAT, MPI_SUM, MPI_COMM_PUP2);
/* rows */
for (i=ns_start; i<=ns_end; i++) {
for (k=1; k<=nx_pot; k++) {
row[k]=data_t1[i][k];
row_im[k]=data_im1[i][k];
}
fft(row, row_im, nx_pot,1);
for (k=1; k<=nx_pot; k++) {
if (i<=ny_pot/2 && k<=nx_pot/2) {
data_vyt[i+ny_pot/2][k+nx_pot/2]=row[k];
data_vyim[i+ny_pot/2][k+nx_pot/2]=row_im[k];
} else if (i<=ny_pot/2 && k>nx_pot/2) {
data_vyt[i+ny_pot/2][k-nx_pot/2]=row[k];
data_vyim[i+ny_pot/2][k-nx_pot/2]=row_im[k];
} else if (i>ny_pot/2 && k<=nx_pot/2) {
data_vyt[i-ny_pot/2][k+nx_pot/2]=row[k];
data_vyim[i-ny_pot/2][k+nx_pot/2]=row_im[k];
} else if (i>ny_pot/2 && k>nx_pot/2) {
data_vyt[i-ny_pot/2][k-nx_pot/2]=row[k];
data_vyim[i-ny_pot/2][k-nx_pot/2]=row_im[k];
}
}
}
}
MPI_Allreduce(MPI_IN_PLACE, &data_vyt[1][1], nx_pot*ny_pot, MPI_FLOAT, MPI_SUM, newcomm_nodentr);
MPI_Allreduce(MPI_IN_PLACE, &data_vyim[1][1], nx_pot*ny_pot, MPI_FLOAT, MPI_SUM, newcomm_nodentr);
if (MYID==0) {
t4=MPI_Wtime();
fprintf(FP,"\ntime for parallised ffts: %f\n",t4-t2);
/* write spectra to file */
if (VERBOSE==1 && ishot==1 && iter==1) {
/* p spectrum */
sprintf(pf2,"%s_spectrum.p.shot%i.it%i.%i.%i.bin",SEIS_FILE,ishot, iter, ny_pot, nx_pot);
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);
/* vz spectrum */
sprintf(pf2,"%s_spectrum.vz.shot%i.it%i.%i.%i.bin",SEIS_FILE,ishot, iter, ny_pot, nx_pot);
file = fopen(pf2,"wb");
for (i=1; i<=nx_pot; i++)
for (j=1; j<=ny_pot; j++) { fwrite(&data_vyt[j][i],sizeof(float),1,file);}
fclose(file);
}
t5=MPI_Wtime();
if (VERBOSE==1) { fprintf(FP,"Perform wavefield separation... \n"); }
/* perform wavefield separation */
for (k=1; k<=nx_pot; k++) {
for (i=ny_pot/2-1000; i<=ny_pot/2+1000; i++) {
data_pup[i][k]=0.5*(data_pt[i][k]-(filtwk_full[i][k]*data_vyt[i][k]));
data_pupim[i][k]=0.5*(data_pim[i][k]-(filtwk_full[i][k]*data_vyim[i][k]));
// if (k<=nx_pot/2 && i<=ny_pot/2) {
// data_pup[i+ny_pot/2][k+nx_pot/2]=0.5*(data_pt[i][k]-(filtwk_full[i][k]*data_vyt[i][k]));
// data_pupim[i+ny_pot/2][k+nx_pot/2]=0.5*(data_pim[i][k]-(filtwk_full[i][k]*data_vyim[i][k]));
// } else if (k>nx_pot/2 && i<=ny_pot/2) {
// data_pup[i+ny_pot/2][k-nx_pot/2]=0.5*(data_pt[i][k]-(filtwk_full[i][k]*data_vyt[i][k]));
// data_pupim[i+ny_pot/2][k-nx_pot/2]=0.5*(data_pim[i][k]-(filtwk_full[i][k]*data_vyim[i][k]));
// } else if (k<=nx_pot/2 && i>ny_pot/2) {
// data_pup[i-ny_pot/2][k+nx_pot/2]=0.5*(data_pt[i][k]-(filtwk_full[i][k]*data_vyt[i][k]));
// data_pupim[i-ny_pot/2][k+nx_pot/2]=0.5*(data_pim[i][k]-(filtwk_full[i][k]*data_vyim[i][k]));
// } else if (k>nx_pot/2 && i>ny_pot/2) {
// data_pup[i-ny_pot/2][k-nx_pot/2]=0.5*(data_pt[i][k]-(filtwk_full[i][k]*data_vyt[i][k]));
// data_pupim[i-ny_pot/2][k-nx_pot/2]=0.5*(data_pim[i][k]-(filtwk_full[i][k]*data_vyim[i][k]));
// }
data_t1[i][k]=0;
data_im1[i][k]=0;
}
}
}
MPI_Bcast(&data_pup[1][1],nx_pot*ny_pot,MPI_FLOAT,0,newcomm_nodentr);
MPI_Bcast(&data_pupim[1][1],nx_pot*ny_pot,MPI_FLOAT,0,newcomm_nodentr);
t6=MPI_Wtime();
fprintf(FP,"\ntime for wavefield separation: %f\n",t6-t5);
/* write pup spectrum to file */
if (VERBOSE==1 && ishot==1 && iter==1 && MYID==0) {
sprintf(pf2,"%s_spectrum.pup.shot%i.it%i.%i.%i.bin",SEIS_FILE,ishot, iter, ny_pot, nx_pot);
file = fopen(pf2,"wb");
for (i=1; i<=nx_pot; i++)
for (j=1; j<=ny_pot; j++) { fwrite(&data_pup[j][i],sizeof(float),1,file); }
fclose(file);
}
// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
t7=MPI_Wtime();
if (VERBOSE==1) { fprintf(FP,"Reverse 2D FFT...\n "); }
ntr_start=nx_pot/(no_myid*2)*(myid_pup)+1;
ntr_end=nx_pot/(no_myid*2)*(myid_pup+1);
ns_start=ny_pot/(no_myid*2)*(myid_pup)+1;
ns_end=ny_pot/(no_myid*2)*(myid_pup+1);
/* transformation into fk domain */
/* columns */
for (k=ntr_start; k<=ntr_end; k++) {
for (i=1; i<=ny_pot; i++) {
// col[i]=data_pt[i][k];
// col_im[i]=data_pim[i][k];
if (i<=ny_pot/2 && k<=nx_pot/2) {
col[i]=data_pt[i+ny_pot/2][k+nx_pot/2];
col_im[i]=data_pim[i+ny_pot/2][k+nx_pot/2];
} else if (i<=ny_pot/2 && k>nx_pot/2) {
col[i]=data_pt[i+ny_pot/2][k-nx_pot/2];
col_im[i]=data_pim[i+ny_pot/2][k-nx_pot/2];
} else if (i>ny_pot/2 && k<=nx_pot/2) {
col[i]=data_pt[i-ny_pot/2][k+nx_pot/2];
col_im[i]=data_pim[i-ny_pot/2][k+nx_pot/2];
} else if (i>ny_pot/2 && k>nx_pot/2) {
col[i]=data_pt[i-ny_pot/2][k-nx_pot/2];
col_im[i]=data_pim[i-ny_pot/2][k-nx_pot/2];
}
}
}
/* expand filter to all 4 quadrants */
for (k=1; k<=ny_pot; k++) {
for (i=1; i<=nx_pot; i++) {
if (k<=ny_pot/2 && i<=nx_pot/2) {
filtwk_full[k][i]=filtwk[(ny_pot/2)-(k-1)][(nx_pot/2)-(i-1)];
} else if (k<=ny_pot/2 && i>nx_pot/2) {
filtwk_full[k][i]=filtwk[(ny_pot/2)-(k-1)][i-(nx_pot/2)];
} else if (k>ny_pot/2 && i<=nx_pot/2) {
filtwk_full[k][i]=filtwk[k-(ny_pot/2)][(nx_pot/2)-(i-1)];
} else if (k>ny_pot/2 && i>nx_pot/2) {
filtwk_full[k][i]=filtwk[k-(ny_pot/2)][i-(nx_pot/2)];
}
}
}
/* write filter to file*/
if (VERBOSE==1 && iter==1 && ishot==1) {
sprintf(pf1,"%s_spectrum.filter.bin",SEIS_FILE);
file = fopen(pf1,"wb");
};
for (i=1; i<=nx_pot; i++)
for (j=1; j<=ny_pot; j++) { fwrite(&filtwk_full[j][i],sizeof(float),1,file); }
fft(col, col_im, ny_pot,-1);
fclose(file);
}
if (VERBOSE==1) { fprintf(FP,"2D FFT of p-data and vz-data... \n"); }
}