taper_grad_shot.c 13.7 KB
Newer Older
Tilman Steinweg's avatar
Tilman Steinweg committed
1
/*-----------------------------------------------------------------------------------------
Florian Wittkamp's avatar
Florian Wittkamp committed
2
 * Copyright (C) 2016  For the list of authors, see file AUTHORS.
Tilman Steinweg's avatar
Tilman Steinweg committed
3
4
 *
 * This file is part of DENISE.
5
 *
Tilman Steinweg's avatar
Tilman Steinweg committed
6
7
8
 * DENISE 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.
9
 *
Tilman Steinweg's avatar
Tilman Steinweg committed
10
11
12
13
 * DENISE 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.
14
 *
Tilman Steinweg's avatar
Tilman Steinweg committed
15
16
 * You should have received a copy of the GNU General Public License
 * along with DENISE. See file COPYING and/or <http://www.gnu.org/licenses/gpl-2.0.html>.
17
 -----------------------------------------------------------------------------------------*/
Tilman Steinweg's avatar
Tilman Steinweg committed
18
19
20
21
22
23

/*
 * taper gradient with a gaussian frame to damp inversion artefacts near the sources and receivers
 sws == 1 vertical taper (for tomography geometry)
 sws == 2 horizontal taper (for reflection geometry)
 sws == 3 local radial taper at the source and receiver positions
24
 sws == 4
Tilman Steinweg's avatar
Tilman Steinweg committed
25
26
27
28
29
30
 */

#include "fd.h"

void taper_grad_shot(float ** waveconv,float ** taper_coeff, float **srcpos, int nshots, int **recpos, int ntr, int ishot, int sws)
{
31
32
33
34
    
    /* extern variables */
    
    extern float DH;
Tilman Steinweg's avatar
Tilman Steinweg committed
35
	extern int FREE_SURF, NX, NY, NXG, NYG;
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
    extern int NPROCX, NPROCY, MYID, POS[3];
    extern FILE *FP;
    extern char TAPER_FILE_NAME[STRING_SIZE], TAPER_FILE_NAME_U[STRING_SIZE], TAPER_FILE_NAME_RHO[STRING_SIZE];
    extern int VERBOSE;
    /* local variables */
    int i, j, h, ifw, ii, jj, n, xb, yb, xe, ye, taperlength,taperlength2, VTON, SRTON;
    int ijc, iy, ix, iii, jjj, xx, yy, srctaper_gridpt, i1, j1;
    
    /*extern int GRADT1, GRADT2, GRADT3, GRADT4;*/
    float amp, a, *window, grad_tap, **waveconvtmp;
    char modfile[STRING_SIZE], taper_file[STRING_SIZE];
    
    extern float SRTRADIUS;
    extern int SRTSHAPE, FILTSIZE;
    float **m, **edgemat, **mm, **msum, minm, maxm, x, y, rad, **taper_coeff_glob;
    float maxrad;
    FILE *fp_taper;
    
    
    
    
    if(sws==1){
        
Tilman Steinweg's avatar
Tilman Steinweg committed
59
60
        /* =================================== */
        /* taper source and receiver positions */
61
62
        /* =================================== */
        
Tilman Steinweg's avatar
Tilman Steinweg committed
63
64
65
        /* Convert from meters to gridpoints -> minimum 5x5 gridpoints */
        srctaper_gridpt = (int)(ceil(2.0*SRTRADIUS/DH));
        if (srctaper_gridpt<5)  srctaper_gridpt = 5;
66
        
Tilman Steinweg's avatar
Tilman Steinweg committed
67
68
69
70
71
        m               = matrix(1,srctaper_gridpt,1,srctaper_gridpt);
        edgemat         = matrix(1,4,1,1);
        mm              = matrix(1,NYG,1,NXG);
        msum            = matrix(1,NYG,1,NXG);
        taper_coeff_glob= matrix(1,NYG,1,NXG);
72
73
        waveconvtmp     = matrix(0,NY+1,0,NX+1);
        
Tilman Steinweg's avatar
Tilman Steinweg committed
74
        for (iy=1;iy<=NYG;iy++)
75
76
            for (ix=1;ix<=NXG;ix++)  msum[iy][ix] = 1.0;
        
Tilman Steinweg's avatar
Tilman Steinweg committed
77
        MPI_Barrier(MPI_COMM_WORLD);
78
        
Tilman Steinweg's avatar
Tilman Steinweg committed
79
80
81
        /*****************************/
        /* Taper at source positions */
        /*****************************/
82
83
        
        a = 1.0;
Tilman Steinweg's avatar
Tilman Steinweg committed
84
85
        maxrad = sqrt(2.0*SRTRADIUS*SRTRADIUS);
        for (j=1;j<=srctaper_gridpt;j++) {
86
87
88
89
90
91
92
93
94
95
96
97
98
            for (i=1;i<=srctaper_gridpt;i++) {
                x = ((float)i-((float)srctaper_gridpt)/2.0-0.5)*DH;
                y = ((float)j-((float)srctaper_gridpt)/2.0-0.5)*DH;
                rad = sqrt(x*x+y*y);
                
                switch (SRTSHAPE) {
                    case 1:
                        m[j][i] = erf(a*rad/maxrad);
                        break;
                    case 2:
                        if (rad>0)      m[j][i] = log(rad);
                        else            m[j][i] = 0.0;
                        break;
Tilman Steinweg's avatar
Tilman Steinweg committed
99
                }
100
            }
Tilman Steinweg's avatar
Tilman Steinweg committed
101
        }
102
        
Tilman Steinweg's avatar
Tilman Steinweg committed
103
104
105
        /* generate local taper matrix */
        minm = minimum_m(m,srctaper_gridpt,srctaper_gridpt);
        for (j=1;j<=srctaper_gridpt;j++)
106
107
            for (i=1;i<=srctaper_gridpt;i++)  m[j][i] -= minm;
        
Tilman Steinweg's avatar
Tilman Steinweg committed
108
109
110
111
112
113
114
115
        /* normalize taper matrix to max of values at the centre of all 4 taper area edges,     */
        /* not the global maximum, which is located at the corners                              */
        edgemat[1][1] = m[1][srctaper_gridpt/2];
        edgemat[2][1] = m[srctaper_gridpt/2][1];
        edgemat[3][1] = m[srctaper_gridpt/2][srctaper_gridpt];
        edgemat[4][1] = m[srctaper_gridpt][srctaper_gridpt/2];
        maxm = maximum_m(edgemat,1,4);
        for (j=1;j<=srctaper_gridpt;j++)
116
117
118
119
            for (i=1;i<=srctaper_gridpt;i++) {
                m[j][i] /= maxm;
                if (m[j][i]>1.0)  m[j][i] = 1.0;
            }
Tilman Steinweg's avatar
Tilman Steinweg committed
120
121
        /* get central position within the taper */
        ijc = (int)(ceil((float)srctaper_gridpt/2));
122
        
Tilman Steinweg's avatar
Tilman Steinweg committed
123
124
125
        /*********************/
        /* loop over sources */
        /*for (n=1;n<=nshots;n++) {*/
126
        n=ishot;
Tilman Steinweg's avatar
Tilman Steinweg committed
127
        for (iy=1;iy<=NYG;iy++)
128
129
            for (ix=1;ix<=NXG;ix++)  mm[iy][ix] = 1.0;
        
Tilman Steinweg's avatar
Tilman Steinweg committed
130
131
132
        i = iround(srcpos[1][n]/DH);
        j = iround(srcpos[2][n]/DH);
        for (iy=1;iy<=srctaper_gridpt;iy++) {
133
134
135
136
137
138
            for (ix=1;ix<=srctaper_gridpt;ix++) {
                xx = i + ix - ijc;
                yy = j + iy - ijc;
                if ((xx<1) || (xx>NXG) || (yy<1) || (yy>NYG))  continue;
                mm[yy][xx] = m[iy][ix];
            }
Tilman Steinweg's avatar
Tilman Steinweg committed
139
        }
140
141
142
143
        
        /*                      for (iy=1;iy<=NYG;iy++)
         for (ix=1;ix<=NXG;ix++)  msum[iy][ix] += mm[iy][ix];
         */
Tilman Steinweg's avatar
Tilman Steinweg committed
144
        for (iy=1;iy<=NYG;iy++)
145
146
147
148
149
150
            for (ix=1;ix<=NXG;ix++)
                if (msum[iy][ix] > mm[iy][ix])
                    msum[iy][ix] = mm[iy][ix];
        
        /* }*/
        
Tilman Steinweg's avatar
Tilman Steinweg committed
151
152
153
        /***********************/
        /* loop over receivers */
        /*for (n=1;n<=ntr;n++) {
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
         for (iy=1;iy<=NYG;iy++)
         for (ix=1;ix<=NXG;ix++)  mm[iy][ix] = 1.0;
         i = recpos[1][n];
         j = recpos[2][n];
         for (iy=1;iy<=srctaper_gridpt;iy++) {
         for (ix=1;ix<=srctaper_gridpt;ix++) {
         xx = i + ix - ijc;
         yy = j + iy - ijc;
         if ((xx<1) || (xx>NXG) || (yy<1) || (yy>NYG))  continue;
         mm[yy][xx] = m[iy][ix];
         }
         }*/
        
        /*                      for (iy=1;iy<=NYG;iy++)    Die kommenden zwei Zeilen wurden von Daniel auskommentiert.
         for (ix=1;ix<=NXG;ix++)  msum[iy][ix] += mm[iy][ix];
         */
        /* for (iy=1;iy<=NYG;iy++)
         for (ix=1;ix<=NXG;ix++)
         if (msum[iy][ix] > mm[iy][ix])
         msum[iy][ix] = mm[iy][ix];
         
         }*/
        
        
Tilman Steinweg's avatar
Tilman Steinweg committed
178
179
        minm = minimum_m(msum,NXG,NYG);
        for (iy=1;iy<=NYG;iy++)
180
181
            for (ix=1;ix<=NXG;ix++)  msum[iy][ix] -= minm;
        
Tilman Steinweg's avatar
Tilman Steinweg committed
182
183
        maxm = maximum_m(msum,NXG,NYG);
        for (iy=1;iy<=NYG;iy++)
184
185
186
187
188
189
190
191
            for (ix=1;ix<=NXG;ix++) {
                taper_coeff_glob[iy][ix] = msum[iy][ix]/maxm;
                if ((POS[1]==((ix-1)/NX)) && (POS[2]==((iy-1)/NY))){
                    ii = ix-POS[1]*NX;
                    jj = iy-POS[2]*NY;
                    /*Diese Zeile wurde von Daniel auskommentiert. taper_coeff[jj][ii] = taper_coeff_glob[iy][ix] * ((float)(iy*DH));*/
                    /*taper_coeff[jj][ii] = ((float)(iy*DH)) * ((float)(iy*DH)) * ((float)(iy*DH));*/
                    taper_coeff[jj][ii]=msum[iy][ix]/maxm;
Tilman Steinweg's avatar
Tilman Steinweg committed
192
                }
193
194
195
196
            }
        
        
        /* apply taper on local gradient */
Tilman Steinweg's avatar
Tilman Steinweg committed
197
        for (j=1;j<=NY;j++){
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
            for (i=1;i<=NX;i++){
                waveconv[j][i]*=taper_coeff[j][i];
                waveconvtmp[j][i] = waveconv[j][i];
            }
        }
        
        
        /* apply filter at shot and receiver points */
        /*for (n=1;n<=nshots;n++)
         {*/
        n=ishot;
        i1 = iround(srcpos[1][n]/DH);
        j1 = iround(srcpos[2][n]/DH);
        /*fprintf(FP,"\n Shot %d (printed by PE %d):\n",n,MYID);
         fprintf(FP,"\n i1: %d (printed by PE %d):\n",i1,MYID);
         fprintf(FP,"\n j1: %d (printed by PE %d):\n",j1,MYID);*/
        
        for (i=i1-FILTSIZE;i<=i1+FILTSIZE;i++){
            for (j=j1-FILTSIZE;j<=j1+FILTSIZE;j++){
                if ((POS[1]==((i-1)/NX)) && (POS[2]==((j-1)/NY))){
                    ii = i-POS[1]*NX;
                    jj = j-POS[2]*NY;
                    /*printf("\n ii: %d (printed by PE %d):\n",ii,MYID);
                     printf("\n jj: %d (printed by PE %d):\n",jj,MYID);*/
                    /*waveconvtmp[jj][ii] = 1*waveconv[jj][ii]
                     + 8*(waveconv[jj-1][ii] + waveconv[jj+1][ii] + waveconv[jj][ii-1] + waveconv[jj][ii+1])
                     + 4*(waveconv[jj-1][ii-1] + waveconv[jj-1][ii+1] + waveconv[jj+1][ii+1] + waveconv[jj+1][ii-1]);
                     waveconvtmp[jj][ii] = waveconvtmp[jj][ii]/49;*/
                    if (jj>0){
                        waveconvtmp[jj][ii] = 0.0;
                        taper_coeff[jj][ii] = 0.0;
                    }
                }
            }
        }
        
        
        /*}*/
        
        
        /*for (n=1;n<=ntr;n++)
         {
         i1 = recpos[1][n];
         j1 = recpos[2][n];
         
         for (i=i1-filtsize;i<=i1+filtsize;i++){
         for (j=j1-filtsize;j<=j1+filtsize;j++){
         
Tilman Steinweg's avatar
Tilman Steinweg committed
246
		       if ((POS[1]==((i-1)/NX)) && (POS[2]==((j-1)/NY))){
247
248
249
250
251
         ii = i-POS[1]*NX;
         jj = j-POS[2]*NY;   */
        /* Die kommenden 4 Zeilen wurden von Daniel auskommentiert. waveconvtmp[jj][ii] = 1*waveconv[jj][ii]
         + 8*(waveconv[jj-1][ii] + waveconv[jj+1][ii] + waveconv[jj][ii-1] + waveconv[jj][ii+1])
         + 4*(waveconv[jj-1][ii-1] + waveconv[jj-1][ii+1] + waveconv[jj+1][ii+1] + waveconv[jj+1][ii-1]);
Tilman Steinweg's avatar
Tilman Steinweg committed
252
			      waveconvtmp[jj][ii] = waveconvtmp[jj][ii]/49;*/
253
254
255
        
        /*	waveconvtmp[jj][ii] = 0.0;
         
Tilman Steinweg's avatar
Tilman Steinweg committed
256
		       }
257
258
259
260
261
         }
         }
         }*/
        
        
Tilman Steinweg's avatar
Tilman Steinweg committed
262
263
        /* apply taper on local gradient */
        for (j=1;j<=NY;j++){
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
            for (i=1;i<=NX;i++){
                waveconv[j][i] = waveconvtmp[j][i];
            }
        }
        
        
        free_matrix(m,1,srctaper_gridpt,1,srctaper_gridpt);
        free_matrix(edgemat,1,4,1,1);
        free_matrix(mm,1,NYG,1,NXG);
        free_matrix(msum,1,NYG,1,NXG);
        free_matrix(taper_coeff_glob,1,NYG,1,NXG);
        free_matrix(waveconvtmp,0,NX+1,0,NY+1);
        
        
        
        
        
        
//        MPI_Barrier(MPI_COMM_WORLD);
//        sprintf(modfile,"taper_coeff_%d.bin",ishot);
//        writemod(modfile,taper_coeff,3);
//        MPI_Barrier(MPI_COMM_WORLD);
//        if (MYID==0) mergemod(modfile,3);
//        MPI_Barrier(MPI_COMM_WORLD);
//        sprintf(modfile,"taper_coeff_%d.bin.%i.%i",ishot,POS[1],POS[2]);
//        remove(modfile);
//        MPI_Barrier(MPI_COMM_WORLD);
//        if(VERBOSE) fprintf(FP,"\n removing successful\n");
    }	/* end of sws==1 */
    
    
    
    
    
    if((sws>=2)&&(sws<=4)){
	if (MYID==0&&VERBOSE)
        {
            fprintf(FP,"\n **Message from taper_grad_shot (printed by PE %d):\n",MYID);
            fprintf(FP," Coefficients for gradient tapers are now read in.\n");
Tilman Steinweg's avatar
Tilman Steinweg committed
303
304
305
        }
        
        if(sws==2){
306
307
308
309
310
            sprintf(taper_file,"%s.shot%d",TAPER_FILE_NAME,ishot);
		if (MYID==0&&VERBOSE)
                fprintf(FP," Coefficients for vp or lambda gradient is read from %s.\n",taper_file);
            fp_taper=fopen(taper_file,"r");}
        
Tilman Steinweg's avatar
Tilman Steinweg committed
311
        if(sws==3){
312
313
314
315
316
            sprintf(taper_file,"%s.shot%i",TAPER_FILE_NAME_U,ishot);
		if (MYID==0&&VERBOSE)
                fprintf(FP," Coefficients for vs or mu gradient is read from %s.\n",taper_file);
            fp_taper=fopen(taper_file,"r");}
        
Tilman Steinweg's avatar
Tilman Steinweg committed
317
        if(sws==4){
318
319
320
321
322
323
            sprintf(taper_file,"%s.shot%i",TAPER_FILE_NAME_RHO,ishot);
		if (MYID==0&&VERBOSE)
                fprintf(FP," Coefficients for rho gradient is read from %s.\n",taper_file);
            fp_taper=fopen(taper_file,"r");}
        
        
Tilman Steinweg's avatar
Tilman Steinweg committed
324
325
326
        /* loop over global grid */
        for (i=1;i<=NXG;i++){
            for (j=1;j<=NYG;j++){
327
                
Tilman Steinweg's avatar
Tilman Steinweg committed
328
                fread(&grad_tap, sizeof(float), 1, fp_taper);
329
                
Tilman Steinweg's avatar
Tilman Steinweg committed
330
                if ((POS[1]==((i-1)/NX)) && (POS[2]==((j-1)/NY))){
331
332
333
334
335
                    ii=i-POS[1]*NX;
                    jj=j-POS[2]*NY;
                    
                    taper_coeff[jj][ii]=grad_tap;
                    
Tilman Steinweg's avatar
Tilman Steinweg committed
336
337
338
                }
            }
        }
339
340
341
342
343
344
345
        
        for (j=1;j<=NY;j++){
            for (i=1;i<=NX;i++){
                waveconv[j][i]*=taper_coeff[j][i];
            }
        }
        
Tilman Steinweg's avatar
Tilman Steinweg committed
346
347
        fclose(fp_taper);
        
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
//        /* output of taper files for testing */
//        if(sws==2){
//            sprintf(modfile,"taper_coeff_file.bin.shot%d",ishot);
//            writemod(modfile,taper_coeff,3);
//            MPI_Barrier(MPI_COMM_WORLD);
//            if (MYID==0) mergemod(modfile,3);
//            MPI_Barrier(MPI_COMM_WORLD);
//            sprintf(modfile,"taper_coeff_file.bin.shot%d.%i.%i",ishot,POS[1],POS[2]);
//            remove(modfile);}
//        
//        if(sws==3){
//            sprintf(modfile,"taper_coeff_file_u.bin.shot%d",ishot);
//            writemod(modfile,taper_coeff,3);
//            MPI_Barrier(MPI_COMM_WORLD);
//            if (MYID==0) mergemod(modfile,3);
//            MPI_Barrier(MPI_COMM_WORLD);
//            sprintf(modfile,"taper_coeff_file_u.bin.shot%d.%i.%i",ishot,POS[1],POS[2]);
//            remove(modfile);}
//        
//        if(sws==4){
//            sprintf(modfile,"taper_coeff_file_rho.bin.shot%d",ishot);
//            writemod(modfile,taper_coeff,3);
//            MPI_Barrier(MPI_COMM_WORLD);
//            if (MYID==0) mergemod(modfile,3);
//            MPI_Barrier(MPI_COMM_WORLD); 
//            sprintf(modfile,"taper_coeff_file_rho.bin.shot%d.%i.%i",ishot,POS[1],POS[2]);
//            remove(modfile);}   
        
    }  /* end of sws==2 || sws==3 || sws==4 */
    
    
    
Tilman Steinweg's avatar
Tilman Steinweg committed
380
381
382
383
}