else{FP=fopen(LOG_FILE,"w");}fprintf(FP," This is the log-file generated by PE %d \n\n",MYID);/* domain decomposition */initproc();NT=iround(TIME/DT);/* number of timesteps *//*ns=iround(NT/NDT);*//* number of samples per trace */ns=NT;/* in a FWI one has to keep all samples of the forward modeled data at the receiver positions to calculate the adjoint sources and to do the backpropagation; look at function saveseis_glob.c to see that every NDT sample for the forward modeled wavefield is written to su files*/lsnap=iround(TSNAP1/DT);/* first snapshot at this timestep */lsamp=NDT;/* output of parameters to log-file or stdout */if(MYID==0)write_par(FP);/* NXG, NYG denote size of the entire (global) grid */NXG=NX;NYG=NY;/* In the following, NX and NY denote size of the local grid ! */NX=IENDX;NY=IENDY;if(SEISMO){recpos=receiver(FP,&ntr);recswitch=ivector(1,ntr);recpos_loc=splitrec(recpos,&ntr_loc,ntr,recswitch);ntr_glob=ntr;ntr=ntr_loc;}/* memory allocation for abort criterion*/L2_hist=vector(1,1000);if(INV_STF)fulldata=matrix(1,ntr_glob,1,NT);/* estimate memory requirement of the variables in megabytes*/switch(SEISMO){case1:/* particle velocities only */nseismograms=2;break;case2:/* pressure only */nseismograms=1;break;case3:/* curl and div only */nseismograms=2;break;case4:/* everything */nseismograms=5;break;case5:/* everything except curl and div */nseismograms=3;break;}/* use only every DTINV time sample for the inversion *//*DTINV=15;*/DTINV_help=ivector(1,NT);NTDTINV=ceil((float)NT/(float)DTINV);/* round towards next higher integer value *//* save every IDXI and IDYI spatial point during the forward modelling */IDXI=1;IDYI=1;/*allocate memory for dynamic, static and buffer arrays */fac1=(NX+FDORDER)*(NY+FDORDER);fac2=sizeof(float)*pow(2.0,-20.0);nd=FDORDER/2+1;// decide how much space for exchange is neededswitch(WAVETYPE){case1:fdo3=2*nd;break;case2:fdo3=1*nd;break;case3:fdo3=3*nd;break;default:fdo3=2*nd;break;}if(L){memdyn=(5.0+3.0*(float)L)*fac1*fac2;memmodel=(12.0+3.0*(float)L)*fac1*fac2;}else{memdyn=5.0*fac1*fac2;memmodel=6.0*fac1*fac2;}memseismograms=nseismograms*ntr*ns*fac2;memfwt=5.0*((NX/IDXI)+FDORDER)*((NY/IDYI)+FDORDER)*NTDTINV*fac2;memfwt1=20.0*NX*NY*fac2;memfwtdata=6.0*ntr*ns*fac2;membuffer=2.0*fdo3*(NY+NX)*fac2;buffsize=2.0*2.0*fdo3*(NX+NY)*sizeof(MPI_FLOAT);memtotal=memdyn+memmodel+memseismograms+memfwt+memfwt1+memfwtdata+membuffer+(buffsize*pow(2.0,-20.0));if(MYID==0&&WAVETYPE==1){fprintf(FP,"\n **Message from main (printed by PE %d):\n",MYID);fprintf(FP," Size of local grids: NX=%d \t NY=%d\n",NX,NY);fprintf(FP," Each process is now trying to allocate memory for:\n");fprintf(FP," Dynamic variables: \t\t %6.2f MB\n",memdyn);fprintf(FP," Static variables: \t\t %6.2f MB\n",memmodel);fprintf(FP," Seismograms: \t\t\t %6.2f MB\n",memseismograms);fprintf(FP," Buffer arrays for grid exchange:%6.2f MB\n",membuffer);fprintf(FP," Network Buffer for MPI_Bsend: \t %6.2f MB\n",buffsize*pow(2.0,-20.0));fprintf(FP," ------------------------------------------------ \n");fprintf(FP," Total memory required: \t %6.2f MB.\n\n",memtotal);}/* allocate buffer for buffering messages */buff_addr=malloc(buffsize);

fprintf(FP,"\n\n\n ------------------------------------------------------------------\n");}countstep=0;if(GRAD_METHOD==1){FWI_run=1;steplength_search=0;gradient_optimization=1;}/*-----------------------------------------------------*//* While loop for Wolfe step length search *//*-----------------------------------------------------*/while(FWI_run||steplength_search||gradient_optimization){/*-----------------------------------------------------*//* Calculate Misfit and gradient *//*-----------------------------------------------------*/if(FWI_run){/* For the calculation of the material parameters between gridpoints they have to be averaged. For this, values lying at 0 and NX+1, for example, are required on the local grid. These are now copied from the neighbouring grids */

}else{if(!ACOUSTIC){matcopy_elastic(prho,ppi,pu);}else{matcopy_acoustic(prho,ppi);}}MPI_Barrier(MPI_COMM_WORLD);/* MPI split for processors with ntr>0 */intmyid_ntr,group_id=0,groupsize;MPI_CommMPI_COMM_NTR;if(ntr)group_id=1;elsegroup_id=0;MPI_Comm_split(MPI_COMM_WORLD,group_id,MYID,&MPI_COMM_NTR);MPI_Comm_rank(MPI_COMM_NTR,&myid_ntr);/* end of MPI split for processors with ntr>0 */

/*------------------------------------------------------------------------------*//*----------- Start of inversion of source time function -----------------------*//*------------------------------------------------------------------------------*/

/*------------------------------------------------------------------------------*//*----------- Start of inversion of source time function -----------------------*//*------------------------------------------------------------------------------*/

}printf("\n ====================================================================================================== \n");printf("\n Time Domain Filter is used for the inversion: lowpass filter, corner frequency of %.2f Hz, order %d\n",FC,ORDER);printf("\n ====================================================================================================== \n");if(iter==1){printf("\n ====================================================== \n");printf("\n MYID = %d: STF inversion at first iteration \n",MYID);}else{printf("\n ================================================================================================ \n");printf("\n MYID = %d: STF inversion because of frequency step at the end of the last iteration \n",MYID);}