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

full parallelisation of pup.c

parent b9056d1e
......@@ -818,10 +818,11 @@ int main(int argc, char **argv){
/* MPI split for processors used for pup parallelisation */
int myid_pup, pupid=0;
int myid_pup, pupid=0, no_myid;
MPI_Comm MPI_COMM_PUP;
if (MYID<=7) pupid = 1;
no_myid=9; /* must be power of 2 +1 */
if (MYID<=no_myid-1) pupid = 1;
else pupid = 0;
MPI_Comm_split(MPI_COMM_WORLD, pupid, MYID, &MPI_COMM_PUP);
MPI_Comm_rank(MPI_COMM_PUP, &myid_pup);
......@@ -1395,7 +1396,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 (pupid) pup(fulldata_p, fulldata_vy, ntr_glob,recpos,ishot,ns,iter,MPI_COMM_PUP,myid_pup);
if (pupid) pup(fulldata_p, fulldata_vy, ntr_glob,recpos,ishot,ns,iter,MPI_COMM_PUP,myid_pup,no_myid);
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 ");
......@@ -1950,7 +1951,7 @@ int main(int argc, char **argv){
/* wavefield separation*/
if (WAVESEP==1) {
if (MYID == 0) fprintf(FP,"\n start wavefield separation...");
if (pupid) pup(fulldata_p, fulldata_vy, ntr_glob,recpos,ishot,ns,iter,MPI_COMM_PUP,myid_pup);
if (pupid) pup(fulldata_p, fulldata_vy, ntr_glob,recpos,ishot,ns,iter,MPI_COMM_PUP,myid_pup,no_myid);
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 ");
......@@ -3553,7 +3554,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 (pupid) pup(fulldata_p, fulldata_vy, ntr_glob,recpos,ishot,ns,iter,MPI_COMM_PUP,myid_pup);
if (pupid) pup(fulldata_p, fulldata_vy, ntr_glob,recpos,ishot,ns,iter,MPI_COMM_PUP,myid_pup,no_myid);
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 ");
......
......@@ -213,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, MPI_Comm newcomm_nodentr, int myid_pup );
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, int no_myid );
float *rd_sour(int *nts,FILE* fp_source);
......
......@@ -35,20 +35,40 @@ void fft(float *array, float *arrayim, int npad, int dir) {
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);
if (dir==1) {
/* write data into complex vector */
for (j=0;j<npad;j++){
in[j][0]=array[j+1];
in[j][1]=arrayim[j+1];
}
/* write output into data matrix */
for(j=0;j<npad;j++){
array[j+1] = out[j][0];
arrayim[j+1] = out[j][1];
p = fftw_plan_dft_1d(npad, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
fftw_execute(p);
/* scramble data (reorder seismogram to put frequency 0 in the middle) */
for (i=0,j=npad/2;j<npad;j++,i++){
array[j+1]=out[i][0];
array[i+1]=out[j][0];
arrayim[j+1]=out[i][1];
arrayim[i+1]=out[j][1];
}
} else {
/* unscramble data */
for (i=0,j=npad/2;j<npad;j++,i++){
in[i][0]=array[j+1];
in[j][0]=array[i+1];
in[i][1]=arrayim[j+1];
in[j][1]=arrayim[i+1];
}
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);
......
This diff is collapsed.
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