From 89eaee659fe189b752b1ca0afaf207ac3bb4bd4c Mon Sep 17 00:00:00 2001
From: daniel-koehn <dkoehn@geophysik.uni-kiel.de>
Date: Thu, 16 Feb 2017 02:32:14 +0100
Subject: [PATCH] Added elastic TTI PSV RTM module
---
include/fd.h | 7 +
src/Makefile | 2 +
src/TTI/RTM_TTI.c | 570 +++++++++++++++++++++++++++++++++++++++++
src/TTI/grad_obj_TTI.c | 327 +++++++++++++++++++++++
src/TTI/physics_TTI.c | 10 +-
5 files changed, 911 insertions(+), 5 deletions(-)
create mode 100644 src/TTI/RTM_TTI.c
create mode 100644 src/TTI/grad_obj_TTI.c
diff --git a/include/fd.h b/include/fd.h
index 0a932fb..5c39e49 100644
--- a/include/fd.h
+++ b/include/fd.h
@@ -306,6 +306,11 @@ void checkfd_ssg_TTI(FILE *fp, float ** prho, float ** c11, float ** c13, float
void FD_TTI();
+float grad_obj_TTI(struct wavePSV *wavePSV, struct wavePSV_PML *wavePSV_PML, struct matTTI *matTTI, struct fwiPSV *fwiPSV, struct mpiPSV *mpiPSV,
+ struct seisPSV *seisPSV, struct seisPSVfwi *seisPSVfwi, struct acq *acq, float *hc, int iter, int nsrc, int ns, int ntr, int ntr_glob, int nsrc_glob,
+ int nsrc_loc, int ntr_loc, int nstage, float **We, float **Ws, float **Wr, float ** taper_coeff, int hin, int *DTINV_help,
+ MPI_Request * req_send, MPI_Request * req_rec);
+
void model_elastic_TTI(float ** rho, float ** c11, float ** c13, float ** c33, float ** c44, float ** theta);
void physics_TTI();
@@ -314,6 +319,8 @@ void readmod_elastic_TTI(float ** rho, float ** c11, float ** c13, float **
void rot_el_tensor_TTI(struct matTTI *matTTI);
+void RTM_TTI();
+
void update_s_elastic_PML_TTI(int nx1, int nx2, int ny1, int ny2,
float ** vx, float ** vy, float ** ux, float ** uy, float ** uxy, float ** uyx, float ** sxx, float ** syy,
float ** sxy, float ** d11, float ** d13, float ** d15, float ** d15h, float ** d33, float ** d35, float ** d35h, float ** d55h,
diff --git a/src/Makefile b/src/Makefile
index 7392a88..34ebbc0 100644
--- a/src/Makefile
+++ b/src/Makefile
@@ -133,10 +133,12 @@ DENISE += alloc_matVTI.c \
DENISE += alloc_matTTI.c \
checkfd_ssg_TTI.c \
FD_TTI.c \
+ grad_obj_TTI.c \
$(MODEL_TTI) \
physics_TTI.c \
readmod_elastic_TTI.c \
rot_el_tensor_TTI.c \
+ RTM_TTI.c \
update_s_elastic_PML_TTI.c \
TTI.c
diff --git a/src/TTI/RTM_TTI.c b/src/TTI/RTM_TTI.c
new file mode 100644
index 0000000..b2a5aa9
--- /dev/null
+++ b/src/TTI/RTM_TTI.c
@@ -0,0 +1,570 @@
+/*
+ * Elastic Reverse Time Migration RTM (2D PSV TTI problem)
+ *
+ * Daniel Koehn
+ * Kiel, 16/02/2017
+ */
+
+#include "fd.h"
+
+void RTM_TTI(){
+
+/* global variables */
+/* ---------------- */
+
+/* forward modelling */
+extern int MYID, FDORDER, NX, NY, NT, L, READMOD, QUELLART, RUN_MULTIPLE_SHOTS, TIME_FILT, READREC;
+extern int LOG, SEISMO, N_STREAMER, FW, NXG, NYG, IENDX, IENDY, NTDTINV, IDXI, IDYI, NXNYI, INV_STF, DTINV;
+extern float FC_SPIKE_1, FC_SPIKE_2, FC, FC_START, TIME, DT;
+extern char LOG_FILE[STRING_SIZE], MFILE[STRING_SIZE];
+extern FILE *FP;
+
+/* gravity modelling/inversion */
+extern int GRAVITY, NZGRAV, NGRAVB, GRAV_TYPE, BACK_DENSITY;
+extern char GRAV_DATA_OUT[STRING_SIZE], GRAV_DATA_IN[STRING_SIZE], GRAV_STAT_POS[STRING_SIZE], DFILE[STRING_SIZE];
+extern float LAM_GRAV, GAMMA_GRAV, LAM_GRAV_GRAD, L2_GRAV_IT1;
+
+/* full waveform inversion */
+extern int GRAD_METHOD, NLBFGS, ITERMAX, IDX, IDY, INVMAT1, EPRECOND;
+extern int GRAD_FORM, POS[3], QUELLTYPB, MIN_ITER, MODEL_FILTER, RTM_SHOT;
+extern float FC_END, PRO, C_vp, C_vs, C_rho;
+extern char MISFIT_LOG_FILE[STRING_SIZE], JACOBIAN[STRING_SIZE];
+extern char *FILEINP1;
+
+/* local variables */
+int ns, nseismograms=0, nt, nd, fdo3, j, i, iter, h, hin, iter_true, SHOTINC, s=0;
+int buffsize, ntr=0, ntr_loc=0, ntr_glob=0, nsrc=0, nsrc_loc=0, nsrc_glob=0, ishot, nshots=0, itestshot;
+
+float sum, eps_scale, opteps_vp, opteps_vs, opteps_rho, Vp_avg, Vs_avg, rho_avg, Vs_sum, Vp_sum, rho_sum;
+char *buff_addr, ext[10], *fileinp, jac[225], source_signal_file[STRING_SIZE];
+
+double time1, time2, time7, time8, time_av_v_update=0.0, time_av_s_update=0.0, time_av_v_exchange=0.0, time_av_s_exchange=0.0, time_av_timestep=0.0;
+
+float L2sum, *L2t;
+
+float ** taper_coeff, * epst1, *hc=NULL;
+int * DTINV_help;
+
+MPI_Request *req_send, *req_rec;
+MPI_Status *send_statuses, *rec_statuses;
+
+/* Variables for step length calculation */
+int step1, step3=0;
+float eps_true, tmp;
+
+/* Variables for the L-BFGS method */
+float * rho_LBFGS, * alpha_LBFGS, * beta_LBFGS;
+float * y_LBFGS, * s_LBFGS, * q_LBFGS, * r_LBFGS;
+int NLBFGS_class, LBFGS_pointer, NLBFGS_vec;
+
+/* Variables for energy weighted gradient */
+float ** Ws, **Wr, **We;
+
+/* parameters for FWI-workflow */
+int stagemax=0, nstage;
+
+/* help variable for MIN_ITER */
+int min_iter_help=0;
+
+/* parameters for gravity inversion */
+char jac_grav[STRING_SIZE];
+
+FILE *FP_stage, *FP_GRAV;
+
+if (MYID == 0){
+ time1=MPI_Wtime();
+ clock();
+}
+
+/* open log-file (each PE is using different file) */
+/* fp=stdout; */
+sprintf(ext,".%i",MYID);
+strcat(LOG_FILE,ext);
+
+if ((MYID==0) && (LOG==1)) FP=stdout;
+else FP=fopen(LOG_FILE,"w");
+fprintf(FP," This is the log-file generated by PE %d \n\n",MYID);
+
+/* ----------------------- */
+/* define FD grid geometry */
+/* ----------------------- */
+
+/* domain decomposition */
+initproc();
+
+NT=iround(TIME/DT); /* number of timesteps */
+
+/* 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;
+
+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;
+
+NXNYI=(NX/IDXI)*(NY/IDYI);
+SHOTINC=1;
+
+/* use only every DTINV time sample for the inversion */
+DTINV_help=ivector(1,NT);
+
+/* Check if RTM workflow-file is defined (stdin) */
+FP_stage=fopen(FILEINP1,"r");
+if(FP_stage==NULL) {
+ if (MYID == 0){
+ printf("\n==================================================================\n");
+ printf(" Cannot open Denise workflow input file %s \n",FILEINP1);
+ printf("\n==================================================================\n\n");
+ err(" --- ");
+ }
+}
+
+fclose(FP_stage);
+
+/* define data structures for PSV problem */
+struct wavePSV;
+struct wavePSV_PML;
+struct matTTI;
+struct fwiPSV;
+struct mpiPSV;
+struct seisPSV;
+struct seisPSVfwi;
+struct acq;
+
+nd = FDORDER/2 + 1;
+fdo3 = 2*nd;
+buffsize=2.0*2.0*fdo3*(NX +NY)*sizeof(MPI_FLOAT);
+
+/* allocate buffer for buffering messages */
+buff_addr=malloc(buffsize);
+if (!buff_addr) err("allocation failure for buffer for MPI_Bsend !");
+MPI_Buffer_attach(buff_addr,buffsize);
+
+/* allocation for request and status arrays */
+req_send=(MPI_Request *)malloc(REQUEST_COUNT*sizeof(MPI_Request));
+req_rec=(MPI_Request *)malloc(REQUEST_COUNT*sizeof(MPI_Request));
+send_statuses=(MPI_Status *)malloc(REQUEST_COUNT*sizeof(MPI_Status));
+rec_statuses=(MPI_Status *)malloc(REQUEST_COUNT*sizeof(MPI_Status));
+
+/* --------- add different modules here ------------------------ */
+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*/
+
+if (SEISMO&&(READREC!=2)){
+
+ acq.recpos=receiver(FP, &ntr, ishot);
+ acq.recswitch = ivector(1,ntr);
+ acq.recpos_loc = splitrec(acq.recpos,&ntr_loc, ntr, acq.recswitch);
+ ntr_glob=ntr;
+ ntr=ntr_loc;
+
+ if(N_STREAMER>0){
+ free_imatrix(acq.recpos,1,3,1,ntr_glob);
+ if(ntr>0) free_imatrix(acq.recpos_loc,1,3,1,ntr);
+ free_ivector(acq.recswitch,1,ntr_glob);
+ }
+
+}
+
+if((N_STREAMER==0)&&(READREC!=2)){
+
+ /* Memory for seismic data */
+ alloc_seisPSV(ntr,ns,&seisPSV);
+
+ /* Memory for FWI seismic data */
+ alloc_seisPSVfwi(ntr,ntr_glob,ns,&seisPSVfwi);
+
+ /* Memory for full data seismograms */
+ alloc_seisPSVfull(&seisPSV,ntr_glob);
+
+}
+
+/* estimate memory requirement of the variables in megabytes*/
+
+switch (SEISMO){
+case 1 : /* particle velocities only */
+ nseismograms=2;
+ break;
+case 2 : /* pressure only */
+ nseismograms=1;
+ break;
+case 3 : /* curl and div only */
+ nseismograms=2;
+ break;
+case 4 : /* everything */
+ nseismograms=5;
+ break;
+}
+
+/* calculate memory requirements for PSV forward problem */
+mem_fwiPSV(nseismograms,ntr,ns,fdo3,nd,buffsize,ntr_glob);
+
+/* Define gradient formulation */
+/* GRAD_FORM = 1 - stress-displacement gradients */
+/* GRAD_FORM = 2 - stress-velocity gradients for decomposed impedance matrix */
+GRAD_FORM = 1;
+
+/* allocate memory for PSV forward problem */
+alloc_PSV(&wavePSV,&wavePSV_PML);
+
+/* calculate damping coefficients for CPMLs (PSV problem)*/
+if(FW>0){PML_pro(wavePSV_PML.d_x, wavePSV_PML.K_x, wavePSV_PML.alpha_prime_x, wavePSV_PML.a_x, wavePSV_PML.b_x, wavePSV_PML.d_x_half, wavePSV_PML.K_x_half, wavePSV_PML.alpha_prime_x_half, wavePSV_PML.a_x_half,
+ wavePSV_PML.b_x_half, wavePSV_PML.d_y, wavePSV_PML.K_y, wavePSV_PML.alpha_prime_y, wavePSV_PML.a_y, wavePSV_PML.b_y, wavePSV_PML.d_y_half, wavePSV_PML.K_y_half, wavePSV_PML.alpha_prime_y_half,
+ wavePSV_PML.a_y_half, wavePSV_PML.b_y_half);
+}
+
+/* allocate memory for TTI material parameters */
+alloc_matTTI(&matTTI);
+
+/* allocate memory for PSV FWI parameters */
+alloc_fwiPSV(&fwiPSV);
+
+/* allocate memory for PSV MPI variables */
+alloc_mpiPSV(&mpiPSV);
+
+taper_coeff= matrix(1,NY,1,NX);
+
+/* memory for source position definition */
+acq.srcpos1=fmatrix(1,8,1,1);
+
+/* memory of L2 norm */
+L2t = vector(1,4);
+epst1 = vector(1,3);
+
+fprintf(FP," ... memory allocation for PE %d was successfull.\n\n", MYID);
+
+/* Holberg coefficients for FD operators*/
+hc = holbergcoeff();
+
+MPI_Barrier(MPI_COMM_WORLD);
+
+/* Reading source positions from SOURCE_FILE */
+acq.srcpos=sources(&nsrc);
+nsrc_glob=nsrc;
+
+
+/* create model grids */
+if (READMOD) readmod_elastic_TTI(matTTI.prho,matTTI.c11,matTTI.c13,matTTI.c33,matTTI.c44,matTTI.theta);
+else model_elastic_TTI(matTTI.prho,matTTI.c11,matTTI.c13,matTTI.c33,matTTI.c44,matTTI.theta);
+
+/* rotate elastic tensor components */
+rot_el_tensor_TTI(&matTTI);
+
+/* check if the FD run will be stable and free of numerical dispersion */
+/* checkfd_ssg_VTI(FP,matVTI.prho,matVTI.c11,matVTI.c13,matVTI.c33,matVTI.c44,hc); */
+
+SHOTINC=1;
+
+/* For RTM read only first line from FWI workflow file and set iter = 1 */
+stagemax = 1;
+iter_true = 1;
+iter = 1;
+
+/* Begin of FWI-workflow */
+for(nstage=1;nstage<=stagemax;nstage++){
+
+/* read workflow input file *.inp */
+FP_stage=fopen(FILEINP1,"r");
+read_par_inv(FP_stage,nstage,stagemax);
+/*fclose(FP_stage);*/
+
+if((EPRECOND==1)||(EPRECOND==3)){
+ Ws = matrix(1,NY,1,NX); /* total energy of the source wavefield */
+ Wr = matrix(1,NY,1,NX); /* total energy of the receiver wavefield */
+ We = matrix(1,NY,1,NX); /* total energy of source and receiver wavefield */
+}
+
+FC=FC_END;
+
+if (MYID==0)
+ {
+ time2=MPI_Wtime();
+ fprintf(FP,"\n\n\n ------------------------------------------------------------------\n");
+ fprintf(FP,"\n\n\n Elastic Reverse Time Migration RTM \n");
+ fprintf(FP,"\n\n\n ------------------------------------------------------------------\n");
+ fprintf(FP,"\n\n");
+ if(RTM_SHOT==0){fprintf(FP,"Output of RTM image after stacking of shots\n");}
+ if(RTM_SHOT==1){fprintf(FP,"Output of RTM images for each shot\n");}
+ }
+
+/* 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. */
+ matcopy_elastic_VTI(matTTI.prho,matTTI.c11,matTTI.d15);
+ matcopy_elastic_VTI(matTTI.prho,matTTI.c11,matTTI.d35);
+ matcopy_elastic_VTI(matTTI.prho,matTTI.c11,matTTI.d55);
+
+MPI_Barrier(MPI_COMM_WORLD);
+
+/* harmonic average of d15, d35, d55 */
+av_harm(matTTI.d15,matTTI.d15h);
+av_harm(matTTI.d35,matTTI.d35h);
+av_harm(matTTI.d55,matTTI.d55h);
+
+/* arithmetic average of density */
+av_rho(matTTI.prho,matTTI.prip,matTTI.prjp);
+
+
+/* Preparing memory variables for update_s (viscoelastic) */
+/*if (L) prepare_update_s_visc_PSV(matPSV.etajm,matPSV.etaip,matPSV.peta,matPSV.fipjp,matPSV.pu,matPSV.puipjp,matPSV.ppi,matPSV.prho,matPSV.ptaus,matPSV.ptaup,matPSV.ptausipjp,matPSV.f,matPSV.g,
+ matPSV.bip,matPSV.bjm,matPSV.cip,matPSV.cjm,matPSV.dip,matPSV.d,matPSV.e);*/
+
+C_vp = 1.0;
+C_vs = 1.0;
+C_rho = 1.0;
+
+/* ---------------------------------------------------------------------------------------------------- */
+/* --------- Calculate RTM P- and S- wave image using the adjoint state method ------------------------ */
+/* ---------------------------------------------------------------------------------------------------- */
+
+L2sum = grad_obj_TTI(&wavePSV, &wavePSV_PML, &matTTI, &fwiPSV, &mpiPSV, &seisPSV, &seisPSVfwi, &acq, hc, iter, nsrc, ns, ntr, ntr_glob,
+nsrc_glob, nsrc_loc, ntr_loc, nstage, We, Ws, Wr, taper_coeff, hin, DTINV_help, req_send, req_rec);
+
+/* Interpolate missing spatial gradient values in case IDXI > 1 || IDXY > 1 */
+/* ------------------------------------------------------------------------ */
+
+if((IDXI>1)||(IDYI>1)){
+
+ interpol(IDXI,IDYI,fwiPSV.waveconv,1);
+ interpol(IDXI,IDYI,fwiPSV.waveconv_u,1);
+ interpol(IDXI,IDYI,fwiPSV.waveconv_rho,1);
+
+}
+
+/* Preconditioning of gradients after shot summation */
+precond_PSV(&fwiPSV,&acq,nsrc,ntr_glob,taper_coeff,FP_GRAV);
+
+/* Output of RTM results */
+if(RTM_SHOT==0){RTM_PSV_out(&fwiPSV);}
+
+} /* End of RTM-workflow loop */
+
+/* deallocate memory for PSV forward problem */
+dealloc_PSV(&wavePSV,&wavePSV_PML);
+
+/* deallocation of memory */
+free_matrix(fwiPSV.Vp0,-nd+1,NY+nd,-nd+1,NX+nd);
+free_matrix(fwiPSV.Vs0,-nd+1,NY+nd,-nd+1,NX+nd);
+free_matrix(fwiPSV.Rho0,-nd+1,NY+nd,-nd+1,NX+nd);
+
+free_matrix(matTTI.prho,-nd+1,NY+nd,-nd+1,NX+nd);
+free_matrix(matTTI.prip,-nd+1,NY+nd,-nd+1,NX+nd);
+free_matrix(matTTI.prjp,-nd+1,NY+nd,-nd+1,NX+nd);
+
+free_matrix(matTTI.c11,-nd+1,NY+nd,-nd+1,NX+nd);
+free_matrix(matTTI.c13,-nd+1,NY+nd,-nd+1,NX+nd);
+free_matrix(matTTI.c33,-nd+1,NY+nd,-nd+1,NX+nd);
+free_matrix(matTTI.c44,-nd+1,NY+nd,-nd+1,NX+nd);
+
+free_matrix(matTTI.d11,-nd+1,NY+nd,-nd+1,NX+nd);
+free_matrix(matTTI.d13,-nd+1,NY+nd,-nd+1,NX+nd);
+free_matrix(matTTI.d15,-nd+1,NY+nd,-nd+1,NX+nd);
+free_matrix(matTTI.d33,-nd+1,NY+nd,-nd+1,NX+nd);
+free_matrix(matTTI.d35,-nd+1,NY+nd,-nd+1,NX+nd);
+free_matrix(matTTI.d55,-nd+1,NY+nd,-nd+1,NX+nd);
+
+free_matrix(matTTI.d15h,-nd+1,NY+nd,-nd+1,NX+nd);
+free_matrix(matTTI.d35h,-nd+1,NY+nd,-nd+1,NX+nd);
+free_matrix(matTTI.d55h,-nd+1,NY+nd,-nd+1,NX+nd);
+
+free_matrix(matTTI.theta,-nd+1,NY+nd,-nd+1,NX+nd);
+
+free_matrix(fwiPSV.prho_old,-nd+1,NY+nd,-nd+1,NX+nd);
+free_matrix(fwiPSV.ppi_old,-nd+1,NY+nd,-nd+1,NX+nd);
+free_matrix(fwiPSV.pu_old,-nd+1,NY+nd,-nd+1,NX+nd);
+free_matrix(fwiPSV.waveconv,-nd+1,NY+nd,-nd+1,NX+nd);
+free_matrix(fwiPSV.waveconv_lam,-nd+1,NY+nd,-nd+1,NX+nd);
+free_matrix(fwiPSV.waveconv_shot,-nd+1,NY+nd,-nd+1,NX+nd);
+
+free_matrix(mpiPSV.bufferlef_to_rig,1,NY,1,fdo3);
+free_matrix(mpiPSV.bufferrig_to_lef,1,NY,1,fdo3);
+free_matrix(mpiPSV.buffertop_to_bot,1,NX,1,fdo3);
+free_matrix(mpiPSV.bufferbot_to_top,1,NX,1,fdo3);
+
+free_vector(hc,0,6);
+
+free_matrix(fwiPSV.gradg,-nd+1,NY+nd,-nd+1,NX+nd);
+free_matrix(fwiPSV.gradp,-nd+1,NY+nd,-nd+1,NX+nd);
+free_matrix(fwiPSV.gradg_rho,-nd+1,NY+nd,-nd+1,NX+nd);
+free_matrix(fwiPSV.gradp_rho,-nd+1,NY+nd,-nd+1,NX+nd);
+free_matrix(fwiPSV.waveconv_rho,-nd+1,NY+nd,-nd+1,NX+nd);
+free_matrix(fwiPSV.waveconv_rho_s,-nd+1,NY+nd,-nd+1,NX+nd);
+free_matrix(fwiPSV.waveconv_rho_shot,-nd+1,NY+nd,-nd+1,NX+nd);
+free_matrix(fwiPSV.gradg_u,-nd+1,NY+nd,-nd+1,NX+nd);
+free_matrix(fwiPSV.gradp_u,-nd+1,NY+nd,-nd+1,NX+nd);
+free_matrix(fwiPSV.waveconv_u,-nd+1,NY+nd,-nd+1,NX+nd);
+free_matrix(fwiPSV.waveconv_mu,-nd+1,NY+nd,-nd+1,NX+nd);
+free_matrix(fwiPSV.waveconv_u_shot,-nd+1,NY+nd,-nd+1,NX+nd);
+
+free_vector(fwiPSV.forward_prop_x,1,NY*NX*NT);
+free_vector(fwiPSV.forward_prop_y,1,NY*NX*NT);
+free_vector(fwiPSV.forward_prop_rho_x,1,NY*NX*NT);
+free_vector(fwiPSV.forward_prop_rho_y,1,NY*NX*NT);
+free_vector(fwiPSV.forward_prop_u,1,NY*NX*NT);
+
+if (nsrc_loc>0){
+ free_matrix(acq.signals,1,nsrc_loc,1,NT);
+ free_matrix(acq.srcpos_loc,1,8,1,nsrc_loc);
+ free_matrix(acq.srcpos_loc_back,1,6,1,nsrc_loc);
+}
+
+ /* free memory for global source positions */
+ free_matrix(acq.srcpos,1,8,1,nsrc);
+
+ /* free memory for source position definition */
+ free_matrix(acq.srcpos1,1,8,1,1);
+
+ /* free memory for abort criterion */
+
+ free_vector(L2t,1,4);
+ free_vector(epst1,1,3);
+
+ if((N_STREAMER==0)||(READREC!=2)){
+
+ if (SEISMO) free_imatrix(acq.recpos,1,3,1,ntr_glob);
+
+ if ((ntr>0) && (SEISMO)){
+
+ free_imatrix(acq.recpos_loc,1,3,1,ntr);
+ acq.recpos_loc = NULL;
+
+ switch (SEISMO){
+ case 1 : /* particle velocities only */
+ free_matrix(seisPSV.sectionvx,1,ntr,1,ns);
+ free_matrix(seisPSV.sectionvy,1,ntr,1,ns);
+ seisPSV.sectionvx=NULL;
+ seisPSV.sectionvy=NULL;
+ break;
+ case 2 : /* pressure only */
+ free_matrix(seisPSV.sectionp,1,ntr,1,ns);
+ break;
+ case 3 : /* curl and div only */
+ free_matrix(seisPSV.sectioncurl,1,ntr,1,ns);
+ free_matrix(seisPSV.sectiondiv,1,ntr,1,ns);
+ break;
+ case 4 : /* everything */
+ free_matrix(seisPSV.sectionvx,1,ntr,1,ns);
+ free_matrix(seisPSV.sectionvy,1,ntr,1,ns);
+ free_matrix(seisPSV.sectionp,1,ntr,1,ns);
+ free_matrix(seisPSV.sectioncurl,1,ntr,1,ns);
+ free_matrix(seisPSV.sectiondiv,1,ntr,1,ns);
+ break;
+
+ }
+
+ }
+
+ free_matrix(seisPSVfwi.sectionread,1,ntr_glob,1,ns);
+ free_ivector(acq.recswitch,1,ntr);
+
+ if((QUELLTYPB==1)||(QUELLTYPB==3)||(QUELLTYPB==5)||(QUELLTYPB==7)){
+ free_matrix(seisPSVfwi.sectionvxdata,1,ntr,1,ns);
+ free_matrix(seisPSVfwi.sectionvxdiff,1,ntr,1,ns);
+ free_matrix(seisPSVfwi.sectionvxdiffold,1,ntr,1,ns);
+ }
+
+ if((QUELLTYPB==1)||(QUELLTYPB==2)||(QUELLTYPB==6)||(QUELLTYPB==7)){
+ free_matrix(seisPSVfwi.sectionvydata,1,ntr,1,ns);
+ free_matrix(seisPSVfwi.sectionvydiff,1,ntr,1,ns);
+ free_matrix(seisPSVfwi.sectionvydiffold,1,ntr,1,ns);
+ }
+
+ if(QUELLTYPB>=4){
+ free_matrix(seisPSVfwi.sectionpdata,1,ntr,1,ns);
+ free_matrix(seisPSVfwi.sectionpdiff,1,ntr,1,ns);
+ free_matrix(seisPSVfwi.sectionpdiffold,1,ntr,1,ns);
+ }
+
+ if(SEISMO){
+ free_matrix(seisPSV.fulldata,1,ntr_glob,1,NT);
+ }
+
+ if(SEISMO==1){
+ free_matrix(seisPSV.fulldata_vx,1,ntr_glob,1,NT);
+ free_matrix(seisPSV.fulldata_vy,1,ntr_glob,1,NT);
+ }
+
+ if(SEISMO==2){
+ free_matrix(seisPSV.fulldata_p,1,ntr_glob,1,NT);
+ }
+
+ if(SEISMO==3){
+ free_matrix(seisPSV.fulldata_curl,1,ntr_glob,1,NT);
+ free_matrix(seisPSV.fulldata_div,1,ntr_glob,1,NT);
+ }
+
+ if(SEISMO==4){
+ free_matrix(seisPSV.fulldata_vx,1,ntr_glob,1,NT);
+ free_matrix(seisPSV.fulldata_vy,1,ntr_glob,1,NT);
+ free_matrix(seisPSV.fulldata_p,1,ntr_glob,1,NT);
+ free_matrix(seisPSV.fulldata_curl,1,ntr_glob,1,NT);
+ free_matrix(seisPSV.fulldata_div,1,ntr_glob,1,NT);
+ }
+
+ }
+
+ free_ivector(DTINV_help,1,NT);
+
+ /* free memory for viscoelastic modeling variables */
+ /*if (L) {
+ free_matrix(matPSV.ptaus,-nd+1,NY+nd,-nd+1,NX+nd);
+ free_matrix(matPSV.ptausipjp,-nd+1,NY+nd,-nd+1,NX+nd);
+ free_matrix(matPSV.ptaup,-nd+1,NY+nd,-nd+1,NX+nd);
+ free_vector(matPSV.peta,1,L);
+ free_vector(matPSV.etaip,1,L);
+ free_vector(matPSV.etajm,1,L);
+ free_vector(matPSV.bip,1,L);
+ free_vector(matPSV.bjm,1,L);
+ free_vector(matPSV.cip,1,L);
+ free_vector(matPSV.cjm,1,L);
+ free_matrix(matPSV.f,-nd+1,NY+nd,-nd+1,NX+nd);
+ free_matrix(matPSV.g,-nd+1,NY+nd,-nd+1,NX+nd);
+ free_matrix(matPSV.fipjp,-nd+1,NY+nd,-nd+1,NX+nd);
+ free_f3tensor(matPSV.dip,-nd+1,NY+nd,-nd+1,NX+nd,1,L);
+ free_f3tensor(matPSV.d,-nd+1,NY+nd,-nd+1,NX+nd,1,L);
+ free_f3tensor(matPSV.e,-nd+1,NY+nd,-nd+1,NX+nd,1,L);
+}*/
+
+/* de-allocate buffer for messages */
+MPI_Buffer_detach(buff_addr,&buffsize);
+
+MPI_Barrier(MPI_COMM_WORLD);
+
+if (MYID==0){
+ fprintf(FP,"\n **Info from main (written by PE %d): \n",MYID);
+ fprintf(FP," CPU time of program per PE: %li seconds.\n",clock()/CLOCKS_PER_SEC);
+ time8=MPI_Wtime();
+ fprintf(FP," Total real time of program: %4.2f seconds.\n",time8-time1);
+ time_av_v_update=time_av_v_update/(double)NT;
+ time_av_s_update=time_av_s_update/(double)NT;
+ time_av_v_exchange=time_av_v_exchange/(double)NT;
+ time_av_s_exchange=time_av_s_exchange/(double)NT;
+ time_av_timestep=time_av_timestep/(double)NT;
+ /* fprintf(FP," Average times for \n");
+ fprintf(FP," velocity update: \t %5.3f seconds \n",time_av_v_update);
+ fprintf(FP," stress update: \t %5.3f seconds \n",time_av_s_update);
+ fprintf(FP," velocity exchange: \t %5.3f seconds \n",time_av_v_exchange);
+ fprintf(FP," stress exchange: \t %5.3f seconds \n",time_av_s_exchange);
+ fprintf(FP," timestep: \t %5.3f seconds \n",time_av_timestep);*/
+
+}
+
+fclose(FP);
+
+
+}
+
+
+
diff --git a/src/TTI/grad_obj_TTI.c b/src/TTI/grad_obj_TTI.c
new file mode 100644
index 0000000..af09deb
--- /dev/null
+++ b/src/TTI/grad_obj_TTI.c
@@ -0,0 +1,327 @@
+/*------------------------------------------------------------------------
+ * Calculate gradient and objective function for 2D PSV TTI problem
+ *
+ *
+ * D. Koehn
+ * Kiel, 16.02.2017
+ * ----------------------------------------------------------------------*/
+
+#include "fd.h"
+
+
+float grad_obj_TTI(struct wavePSV *wavePSV, struct wavePSV_PML *wavePSV_PML, struct matTTI *matTTI, struct fwiPSV *fwiPSV, struct mpiPSV *mpiPSV,
+ struct seisPSV *seisPSV, struct seisPSVfwi *seisPSVfwi, struct acq *acq, float *hc, int iter, int nsrc, int ns, int ntr, int ntr_glob, int nsrc_glob,
+ int nsrc_loc, int ntr_loc, int nstage, float **We, float **Ws, float **Wr, float ** taper_coeff, int hin, int *DTINV_help,
+ MPI_Request * req_send, MPI_Request * req_rec){
+
+ /* global variables */
+ extern int MYID, TIME_FILT, IDX, IDY, NX, NY, NT, RUN_MULTIPLE_SHOTS, INV_STF, QUELLART;
+ extern int TESTSHOT_START, TESTSHOT_END, TESTSHOT_INCR, SEISMO, EPRECOND, LNORM, READREC;
+ extern int N_STREAMER, SWS_TAPER_CIRCULAR_PER_SHOT, QUELLTYPB, QUELLTYP, LOG;
+ extern int ORDER_SPIKE, ORDER, SHOTINC, RTM_SHOT;
+ extern float EPSILON, FC, FC_START, FC_SPIKE_1, FC_SPIKE_2;
+ extern float C_vp, C_vs, C_rho;
+ extern char MFILE[STRING_SIZE];
+
+ /* local variables */
+ int i, j, nshots, ishot, nt, lsnap, itestshot, swstestshot;
+ float L2sum, L2_all_shots, energy_all_shots, energy_tmp, L2_tmp;
+ char source_signal_file[STRING_SIZE];
+
+ FILE *FP;
+
+ if ((MYID==0) && (LOG==1)) FP=stdout;
+
+ /* initialization of L2 calculation */
+ (*seisPSVfwi).L2=0.0;
+ (*seisPSVfwi).energy=0.0;
+ L2_all_shots=0.0;
+ energy_all_shots=0.0;
+
+ EPSILON=0.0; /* test step length */
+
+ /* set gradient and preconditioning matrices 0 before next iteration*/
+ init_grad((*fwiPSV).waveconv);
+ init_grad((*fwiPSV).waveconv_rho);
+ init_grad((*fwiPSV).waveconv_u);
+
+ itestshot=TESTSHOT_START;
+ swstestshot=0;
+ SHOTINC=1;
+
+ if (RUN_MULTIPLE_SHOTS) nshots=nsrc; else nshots=1;
+
+ for (ishot=1;ishot<=nshots;ishot+=SHOTINC){
+ /*for (ishot=1;ishot<=1;ishot+=1){*/
+
+ /*initialize gradient matrices for each shot with zeros*/
+ init_grad((*fwiPSV).waveconv_shot);
+ init_grad((*fwiPSV).waveconv_u_shot);
+ init_grad((*fwiPSV).waveconv_rho_shot);
+
+ if((EPRECOND==1)||(EPRECOND==3)){
+ init_grad(Ws);
+ init_grad(Wr);
+ init_grad(We);
+ }
+
+ if((N_STREAMER>0)||(READREC==2)){
+
+ if (SEISMO){
+ (*acq).recpos=receiver(FP, &ntr, ishot);
+ (*acq).recswitch = ivector(1,ntr);
+ (*acq).recpos_loc = splitrec((*acq).recpos,&ntr_loc, ntr, (*acq).recswitch);
+ ntr_glob=ntr;
+ ntr=ntr_loc;
+ }
+
+ /* Memory for seismic data */
+ alloc_seisPSV(ntr,ns,seisPSV);
+
+ /* Memory for full data seismograms */
+ alloc_seisPSVfull(seisPSV,ntr_glob);
+
+ /* Memory for FWI seismic data */
+ alloc_seisPSVfwi(ntr,ntr_glob,ns,seisPSVfwi);
+
+ }
+
+ for (nt=1;nt<=8;nt++) (*acq).srcpos1[nt][1]=(*acq).srcpos[nt][ishot];
+
+ /* set QUELLTYP for each shot */
+ QUELLTYP = (*acq).srcpos[8][ishot];
+
+ if (RUN_MULTIPLE_SHOTS){
+
+ /* find this single source positions on subdomains */
+ if (nsrc_loc>0) free_matrix((*acq).srcpos_loc,1,8,1,1);
+ (*acq).srcpos_loc=splitsrc((*acq).srcpos1,&nsrc_loc, 1);
+ }
+
+
+ else{
+ /* Distribute multiple source positions on subdomains */
+ (*acq).srcpos_loc = splitsrc((*acq).srcpos,&nsrc_loc, nsrc);
+ }
+
+ MPI_Barrier(MPI_COMM_WORLD);
+
+ /*==================================================================================
+ Estimate source time function by Wiener deconvolution
+ ==================================================================================*/
+
+ /*if((INV_STF)&&(iter==1)){
+ stf_psv(wavePSV,wavePSV_PML,matPSV,fwiPSV,mpiPSV,seisPSV,seisPSVfwi,acq,hc,ishot,nshots,nsrc_loc,nsrc,ns,ntr,ntr_glob,iter,Ws,Wr,hin,DTINV_help,req_send,req_rec);
+ }*/
+
+ /*==================================================================================
+ Starting simulation (forward model)
+ ==================================================================================*/
+
+ /* calculate wavelet for each source point */
+ (*acq).signals=NULL;
+ (*acq).signals=wavelet((*acq).srcpos_loc,nsrc_loc,ishot);
+
+ if (nsrc_loc){if(QUELLART==6){
+
+ /* time domain filtering of the source signal */
+ apply_tdfilt((*acq).signals,nsrc_loc,ns,ORDER_SPIKE,FC_SPIKE_2,FC_SPIKE_1);
+
+ }
+ }
+
+ /* time domain filtering*/
+ if ((TIME_FILT)&&(INV_STF==0)){
+
+ /* time domain filtering of the source signal */
+ apply_tdfilt((*acq).signals,nsrc_loc,ns,ORDER,FC,FC_START);
+
+ }
+
+ /*printf("MYID=%d, nsrc_loc = %d \n",MYID,nsrc_loc);*/
+
+ /*char source_signal_file[STRING_SIZE];
+ sprintf(source_signal_file,"source_signal.%d.su.shot%d.it%d",MYID,ishot,iter);
+ fprintf(stdout,"\n PE %d outputs source time function in SU format to %s \n ", MYID, source_signal_file);
+ output_source_signal(fopen(source_signal_file,"w"),signals,NT,3);*/
+
+ /* output source signal e.g. for cross-correlation of comparison with analytical solutions */
+ if(RUN_MULTIPLE_SHOTS){
+
+ if(nsrc_loc>0){
+ sprintf(source_signal_file,"%s_source_signal.%d.su.shot%d", MFILE, MYID,ishot);
+ fprintf(stdout,"\n PE %d outputs source time function in SU format to %s \n ", MYID, source_signal_file);
+ output_source_signal(fopen(source_signal_file,"w"),(*acq).signals,NT,1);
+ }
+
+ MPI_Barrier(MPI_COMM_WORLD);
+ }
+
+ /* solve forward problem */
+ TTI(wavePSV,wavePSV_PML,matTTI,fwiPSV,mpiPSV,seisPSV,seisPSVfwi,acq,hc,ishot,nshots,nsrc_loc,ns,ntr,Ws,Wr,hin,DTINV_help,0,req_send,req_rec);
+
+
+ /* ===============================================
+ Calculate objective function and data residuals
+ =============================================== */
+
+ if((ishot==itestshot)&&(ishot<=TESTSHOT_END)){swstestshot=1;}
+
+ if (ntr > 0){
+ calc_res_PSV(seisPSV,seisPSVfwi,(*acq).recswitch,(*acq).recpos,(*acq).recpos_loc,ntr_glob,ntr,nsrc_glob,(*acq).srcpos,ishot,ns,iter,swstestshot);
+ }
+
+ if((ishot==itestshot)&&(ishot<=TESTSHOT_END)){
+ swstestshot=0;
+ itestshot+=TESTSHOT_INCR;
+ }
+
+ /* output of time reversed residual seismograms */
+ if ((SEISMO)&&(iter==1)&&(ishot==1)){
+ outseis_PSVres(seisPSV,seisPSVfwi,(*acq).recswitch,(*acq).recpos,(*acq).recpos_loc,ntr_glob,(*acq).srcpos,ishot,ns,nstage,FP);
+ }
+
+ /*================================================================================
+ Starting simulation (backward model)
+ ==================================================================================*/
+
+ /* Distribute multiple source positions on subdomains */
+ /* define source positions at the receivers */
+ (*acq).srcpos_loc_back = matrix(1,6,1,ntr);
+ for (i=1;i<=ntr;i++){
+ (*acq).srcpos_loc_back[1][i] = ((*acq).recpos_loc[1][i]);
+ (*acq).srcpos_loc_back[2][i] = ((*acq).recpos_loc[2][i]);
+ }
+
+ /* backpropagate time-reversed data */
+ TTI(wavePSV,wavePSV_PML,matTTI,fwiPSV,mpiPSV,seisPSV,seisPSVfwi,acq,hc,ishot,nshots,ntr,ns,ntr,Ws,Wr,hin,DTINV_help,1,req_send,req_rec);
+
+ /* assemble PSV VTI gradients for each shot */
+ /*ass_gradVTI(fwiPSV,matVTI,iter);*/
+
+ if((EPRECOND==1)||(EPRECOND==3)){
+ /* calculate energy weights */
+ eprecond1(We,Ws,Wr);
+
+ /* scale gradient with energy weights*/
+ for(i=1;i<=NX;i=i+IDX){
+ for(j=1;j<=NY;j=j+IDY){
+
+ (*fwiPSV).waveconv_shot[j][i] = (*fwiPSV).waveconv_shot[j][i]/(We[j][i]*C_vp*C_vp);
+ (*fwiPSV).waveconv_u_shot[j][i] = (*fwiPSV).waveconv_u_shot[j][i]/(We[j][i]*C_vs*C_vs);
+ if(C_vs==0.0){(*fwiPSV).waveconv_u_shot[j][i] = 0.0;}
+ (*fwiPSV).waveconv_rho_shot[j][i] = (*fwiPSV).waveconv_rho_shot[j][i]/(We[j][i]*C_rho*C_rho);
+
+ }
+ }
+ }
+
+ if (SWS_TAPER_CIRCULAR_PER_SHOT){ /* applying a circular taper at the source position to the gradient of each shot */
+
+ /* applying the preconditioning */
+ taper_grad_shot((*fwiPSV).waveconv_shot,taper_coeff,(*acq).srcpos,nsrc,(*acq).recpos,ntr_glob,ishot);
+ taper_grad_shot((*fwiPSV).waveconv_rho_shot,taper_coeff,(*acq).srcpos,nsrc,(*acq).recpos,ntr_glob,ishot);
+ taper_grad_shot((*fwiPSV).waveconv_u_shot,taper_coeff,(*acq).srcpos,nsrc,(*acq).recpos,ntr_glob,ishot);
+
+ } /* end of SWS_TAPER_CIRCULAR_PER_SHOT == 1 */
+
+ for(i=1;i<=NX;i=i+IDX){
+ for(j=1;j<=NY;j=j+IDY){
+ (*fwiPSV).waveconv[j][i] += (*fwiPSV).waveconv_shot[j][i];
+ (*fwiPSV).waveconv_rho[j][i] += (*fwiPSV).waveconv_rho_shot[j][i];
+ (*fwiPSV).waveconv_u[j][i] += (*fwiPSV).waveconv_u_shot[j][i];
+ }
+ }
+
+ if(RTM_SHOT==1){RTM_PSV_out_shot(fwiPSV,ishot);}
+
+ if((N_STREAMER>0)||(READREC==2)){
+
+ if (SEISMO) free_imatrix((*acq).recpos,1,3,1,ntr_glob);
+
+ if ((ntr>0) && (SEISMO)){
+
+ free_imatrix((*acq).recpos_loc,1,3,1,ntr);
+ (*acq).recpos_loc = NULL;
+
+ switch (SEISMO){
+ case 1 : /* particle velocities only */
+ free_matrix((*seisPSV).sectionvx,1,ntr,1,ns);
+ free_matrix((*seisPSV).sectionvy,1,ntr,1,ns);
+ (*seisPSV).sectionvx=NULL;
+ (*seisPSV).sectionvy=NULL;
+ break;
+ case 2 : /* pressure only */
+ free_matrix((*seisPSV).sectionp,1,ntr,1,ns);
+ break;
+ case 3 : /* curl and div only */
+ free_matrix((*seisPSV).sectioncurl,1,ntr,1,ns);
+ free_matrix((*seisPSV).sectiondiv,1,ntr,1,ns);
+ break;
+ case 4 : /* everything */
+ free_matrix((*seisPSV).sectionvx,1,ntr,1,ns);
+ free_matrix((*seisPSV).sectionvy,1,ntr,1,ns);
+ free_matrix((*seisPSV).sectionp,1,ntr,1,ns);
+ free_matrix((*seisPSV).sectioncurl,1,ntr,1,ns);
+ free_matrix((*seisPSV).sectiondiv,1,ntr,1,ns);
+ break;
+
+ }
+
+ }
+
+ free_matrix((*seisPSVfwi).sectionread,1,ntr_glob,1,ns);
+ free_ivector((*acq).recswitch,1,ntr);
+
+ if((QUELLTYPB==1)||(QUELLTYPB==3)||(QUELLTYPB==5)||(QUELLTYPB==7)){
+ free_matrix((*seisPSVfwi).sectionvxdata,1,ntr,1,ns);
+ free_matrix((*seisPSVfwi).sectionvxdiff,1,ntr,1,ns);
+ free_matrix((*seisPSVfwi).sectionvxdiffold,1,ntr,1,ns);
+ }
+
+ if((QUELLTYPB==1)||(QUELLTYPB==2)||(QUELLTYPB==6)||(QUELLTYPB==7)){
+ free_matrix((*seisPSVfwi).sectionvydata,1,ntr,1,ns);
+ free_matrix((*seisPSVfwi).sectionvydiff,1,ntr,1,ns);
+ free_matrix((*seisPSVfwi).sectionvydiffold,1,ntr,1,ns);
+ }
+
+ if(QUELLTYPB>=4){
+ free_matrix((*seisPSVfwi).sectionpdata,1,ntr,1,ns);
+ free_matrix((*seisPSVfwi).sectionpdiff,1,ntr,1,ns);
+ free_matrix((*seisPSVfwi).sectionpdiffold,1,ntr,1,ns);
+ }
+
+ ntr=0;
+ ntr_glob=0;
+
+ }
+
+ nsrc_loc=0;
+
+ } /* end of loop over shots (forward and backpropagation) */
+
+ /* calculate L2 norm of all CPUs*/
+ L2sum = 0.0;
+ L2_tmp = (*seisPSVfwi).L2;
+ MPI_Allreduce(&L2_tmp,&L2sum,1,MPI_FLOAT,MPI_SUM,MPI_COMM_WORLD);
+
+ /* calculate L2 norm of all CPUs*/
+ energy_all_shots = 0.0;
+ energy_tmp = (*seisPSVfwi).energy;
+ MPI_Allreduce(&energy_tmp,&energy_all_shots,1,MPI_FLOAT,MPI_SUM,MPI_COMM_WORLD);
+
+ /*if(MYID==0){
+ printf("L2sum: %e\n", L2sum);
+ printf("energy_sum: %e\n\n", energy_all_shots);
+
+ }*/
+
+ if(LNORM==2){
+ L2sum = L2sum/energy_all_shots;
+ }
+ else{L2sum=L2sum;}
+
+return L2sum;
+
+}
+
diff --git a/src/TTI/physics_TTI.c b/src/TTI/physics_TTI.c
index 18dbf32..e1ca43f 100644
--- a/src/TTI/physics_TTI.c
+++ b/src/TTI/physics_TTI.c
@@ -12,7 +12,7 @@ void physics_TTI(){
/* global variables */
extern int MODE;
- /* 2D TTI Forward Problem */
+ /* 2D PSV TTI Forward Problem */
if(MODE==0){
FD_TTI();
}
@@ -22,10 +22,10 @@ void physics_TTI(){
FWI_PSV();
}*/
- /* 2D PSV Reverse Time Migration */
- /*if(MODE==2){
- RTM_PSV();
- }*/
+ /* 2D PSV TTI Reverse Time Migration */
+ if(MODE==2){
+ RTM_TTI();
+ }
}