frechetslave.c

su 的源代码库

/* Copyright (c) Colorado School of Mines, 2006.*/
/* All rights reserved.                       */

/*                                                              */
/*  Function gradient()                                         */
/*                                                              */
/*  Compute the FRECHET derivatives for the elastic             */
/*  modeling. Distributed implementation.                       */
/*                                                              */
/*  Input parameters:                                           */
/*  alpha..................p-wave velocities of the model       */
/*                         global variable                      */
/*  beta...................s-wave velocities of the model       */
/*                         global variable                      */
/*  rho....................densities of the elastic model       */
/*                         global variable                      */
/*  thick..................thicknesses of the model             */
/*                         global variable                      */
/*  res....................synthetic model generated for model  */
/*                         global                               */
/*  CD.....................data covariance matrix               */
/*                         global                               */
/*                                                              */
/*  Output parameters:                                          */
/*  grad...................the gradient, a vector which         */
/*                         dimension is the number of layers    */
/*                                                              */
/*                                              Wences Gouveia  */
/*                                              September 1995  */
#include "frechetSlave.h"
void main()
{
   /* declaration of variables */
   FILE *fp;                       /* file pointer for receiver array */
   char *recFile = " ";            /* file with receiver coordinates */
   int i, indexF, iF, iR, iU, iDer, iL;
                                   /* counters */
   int masterId;                   /* master id */
   int die;                        /* if die=1 process stops */
   int directWave;                 /* direct wave flag */
   int nU;                         /* number of slownesses */
   int wL;                         /* taper length */
   int wInfo[2];                   /* frequency limits */
   int verbose;                    /* dialogue flag */ 
   float *grad;                    /* partial gradient */
   float cte;                      /* a constant */
   float tau;                      /* magnitude of wrap-around attenuation */
   float percU;                    /* percentual of slowness */
                                   /* domain that are windowed */
   float percW;                    /* percentual of frequency */
                                   /* domain that are windowed */
   float r1, dR;                   /* defines receiver array */
   float u1, u2;                   /* slowness window */
   float u, w, f;                  /* slowness & frequency */
   float phi;                      /* azimuth angle */
   float F1, F2, F3;               /* source components */
   float epslon1, epslon2;         /* auxiliary quantities */
   float dU;                       /* slowness interval */
   float wRef;                     /* reference frequency */     
   complex irr[2][2];              /* = (I - R-R+) */
   complex irrI[2][2];             /* = (I - R-R+)^-1 */                     
   complex muC;                    /* auxiliar variable */
   complex dUC;                    /* slowness complex interval */
   complex **resCD;                /* current residual dotted */
                                   /* into covariance */ 
   complex dUCEp1, dUCEp2;         /* dUC * epslon1 and dUC * epslon2 */     
   complex wCCte;                  /* auxiliar variable */ 
   complex am;                     /* vertical P-wave slownesses */
   complex amInv;                  /* 1. / am */
   complex amI;                    /* amI = am * I */
   complex bm;                     /* vertical S-wave slownesses */
   complex bmInv;                  /* 1. / bm */
   complex bmI;                    /* bmI = bm * I */
   complex As1, As2;               /* amplitudes of plane wave components (P)*/
   complex Cs1, Cs2;               /* amplitudes of plane wave components (S)*/
                                   /* downgoing waves */
   complex Bs1, Bs2;               /* amplitudes of plane wave components (P)*/
   complex Ds1, Ds2;               /* amplitudes of plane wave components (S)*/
                                   /* upgoing waves */
   complex g[2];                   /* phase-shift vector */ 
   complex displ;                  /* Frechet derivative of displacements */
   INFO info[1];                   /* basic information for slaves */
   
   SeisSlave logInfo;              /* report structure */

   /* logging information */
   InitLog(&logInfo, PROCESS_FRECHET);
    
   /* receiving INFO data structure */
   masterId = RecvINFO(info, 1, -1, GENERAL_INFORMATION);

   recFile = info->recFile;
   directWave = info->directWave;
   nL = info->nL;
   nF = info->nF;
   nR = info->nR;
   nU = info->nU;
   hanningFlag = info->hanningFlag;
   nSamples = info->nSamples;
   numberPar = info->numberPar;
   limRange = info->limRange;
   lim[0] = info->lim[0];
   lim[1] = info->lim[1];
   VERTICAL = info->vertical;
   RADIAL = info->radial;
   r1 = info->r1;
   dR = info->dR;
   zs = info->zs;
   u1 = info->u1;
   u2 = info->u2;
   dU = info->dU;
   f1 = info->f1;
   f2 = info->f2;
   dF = info->dF;
   F1 = info->F1;
   F2 = info->F2;
   F3 = info->F3;
   percU = info->percU;
   percW = info->percW;
   wR = info->wR;
   tau = info->tau;
   vpFrechet = info->vpFrechet;
   vsFrechet = info->vsFrechet;
   rhoFrechet = info->rhoFrechet;

   if (verbose)
   {
      fprintf(logInfo.fp_log, "Starting modeling for %s\n", info->id);
      fflush(logInfo.fp_log);
   }

   /* DD 
   fprintf(logInfo.fp_log, "numberPar %d\n", numberPar);
   fprintf(logInfo.fp_log, "lim0 %d lim1 %d\n", lim[0], lim[1]);
   fprintf(logInfo.fp_log, "limRange %d\n", limRange);*/
   
   /* memory allocation */
   thick = alloc1float(nL + 1);
   alpha = alloc1float(nL + 1);
   beta = alloc1float(nL + 1);
   rho = alloc1float(nL + 1);
   qP = alloc1float(nL + 1);
   qS = alloc1float(nL + 1);
   recArray = alloc1float(nR);
   grad = alloc1float(numberPar * limRange);
   v1 = alloc2complex(2, numberPar * limRange + 1);
   v2 = alloc2complex(2, numberPar * limRange + 1);
   DmB = alloc3complex(4, numberPar * (limRange + 2), nL);
   derFactor = alloc2complex(2, nL + 1);
   aux11 = alloc2complex(nR, numberPar * limRange);
   aux12 = alloc2complex(nR, numberPar * limRange);
   aux21 = alloc2complex(nR, numberPar * limRange);
   aux22 = alloc2complex(nR, numberPar * limRange);
   aux11Old = alloc2complex(nR, numberPar * limRange);
   aux12Old = alloc2complex(nR, numberPar * limRange);
   aux21Old = alloc2complex(nR, numberPar * limRange);
   aux22Old = alloc2complex(nR, numberPar * limRange);
        
   PSlowness = alloc2complex(2, nL + 1);
   SSlowness = alloc2complex(2, nL + 1);
   S2Velocity = alloc2complex(2, nL + 1);                  
   
   /* DD
   fprintf(logInfo.fp_log,  
	   "receiving model info \n");
   fflush(logInfo.fp_log);*/

   /* receiving elastic model */
   RecvFloat(thick, nL + 1, -1, THICKNESS);
   RecvFloat(rho, nL + 1, -1, DENSITY);
   RecvFloat(alpha, nL + 1, -1, ALPHAS);
   RecvFloat(qP, nL + 1, -1, QALPHA);
   RecvFloat(beta, nL + 1, -1, BETAS);
   RecvFloat(qS, nL + 1, -1, QBETA);

   /* DD 
   fprintf(logInfo.fp_log,  
	   "received model info \n");
   fflush(logInfo.fp_log);*/

   /* source is at 1st layer */
   thick[0] -= zs;

   /* computing the window length for the slowness domain */
   epslon1 = (u2 - u1) * percU;
   wL = NINT(epslon1 / dU);   wL = 2 * wL + 1;
   u2 += epslon1;
   nU = NINT((u2 - u1) / dU);   /* new nU to preserve last slowness */
                                /* w/o being windowed */     
  /* building window for slowness integration */
   taper = alloc1float(nU);
   for (i = (wL - 1) / 2, iU = 0; iU < nU; iU++)
   {
      taper[iU] = 1;
      if (iU >= nU - (wL - 1) / 2)
      {
	 i++;
	 taper[iU] =
	    .42 - .5 * cos(2 * PI * (float) i / ((float) (wL - 1))) +
	       .08 * cos(4 * PI * (float) i / ((float) (wL - 1)));
      }
   }
         
   /* building windowing operator */
   /* DD
      fprintf(logInfo.fp_log, "going to filter\n");
      fflush(logInfo.fp_log);*/
   filter(percW);

   /* I will assume that the receivers are in line (at z = 0) so phi = 0 */
   phi = 0;
   epslon1 = F3;
   epslon2 = F1 * cos(phi) + F2 * sin(phi);

   /* auxiliar constant */
   cte = 1. / (4 * PI * rho[0]);

   /* normalization for the complex slowness */
   if (f1 > 7.5)
      wRef = f1 * 2 * PI;
   else
      wRef = 7.5 * 2 * PI;
   
   /* specific geometry */
   fp = fopen(recFile, "r");
   if (fp == NULL)
   {
      /* standard end-on */
      for (iR = 0; iR < nR; iR++)
      {
	 recArray[iR] = r1 + iR * dR;
      }
   }
   else
   {
      for (iR = 0; iR < nR; iR++)
      {
	 fscanf(fp, "%f\n", &recArray[iR]);
      }
   }
   fclose(fp);
   
   /* loop over frequencies */
   do 
   {
      /* reseting gradient */
      for (iDer = 0; iDer < numberPar * limRange; iDer++)
	 grad[iDer] = 0;
   
      /* DD 
      fprintf(logInfo.fp_log,  
	      "receiving frequency limits and correlation matrix\n");
      fflush(logInfo.fp_log);*/

      /* receiving frequency limits and covariance matrix */
      masterId = RecvInt(wInfo, 2, -1, FREQUENCY_LIMITS);
      resCD = alloc2complex(wInfo[1] - wInfo[0] + 1, nR);
      masterId = RecvCplx(resCD[0], nR * (wInfo[1] - wInfo[0] + 1), -1,
			  COVARIANCE_PARTITION);
      if (verbose)
      {
	 fprintf(logInfo.fp_log,
		 "Frequency limits received: [%d, %d] : [%f, %f]\n",
		 wInfo[0], wInfo[1], wInfo[0] * dF, wInfo[1] * dF);
	 fflush(logInfo.fp_log);
      }

      /* first frequency */
      f1 = wInfo[0] * dF;

      for (indexF = NINT(f1 / dF), f = f1, iF = 0;
	   iF < (wInfo[1] - wInfo[0] + 1); iF++, f += dF, indexF++)
      {
	 if (verbose)
	 {
	    fprintf(logInfo.fp_log, 
		    "Slave modeling Frechet frequency (Hz) : %f: \n", f);
	    fflush(logInfo.fp_log);
	 }
	 
	 /* reseting */   
	 for (i = 0; i < numberPar * limRange; i++)
	 {
	    for (iR = 0; iR < nR; iR++)
	    {  
	       aux11[i][iR] = zeroC;	         aux12[i][iR] = zeroC;
	       aux21[i][iR] = zeroC;	         aux22[i][iR] = zeroC;
	       aux11Old[i][iR] = zeroC;	         aux12Old[i][iR] = zeroC;
	       aux21Old[i][iR] = zeroC;	         aux22Old[i][iR] = zeroC;
	    }
	 }
      
	 w = 2 * PI * f;
	 wC.r = w; wC.i = -tau;
      
	 /* module and phase of complex frequency */
	 wCR = sqrt(wC.r * wC.r + wC.i * wC.i);
	 wCP = atan2(wC.i, wC.r);
	 
	 /* complex slowness step */
	 dUC.r = w * dU / wCR;
	 dUC.i = info->tau * info->dU / wCR;
	 
	 /* wCR / wR */
	 wCRwR = wCR / info->wR;
      
	 /* auxiliary variable */
	 wCCte.r = wC.r * cte;
	 wCCte.i = wC.i * cte;

	 /* compute frequency-dependent horizontal slownesses (squared) */
	 /* and also the s-wave VELOCITIES (squared) for all layers */
	 horSlownessFrechet();
      
	 for (u = u1, iU = 0; iU < nU; iU++, u += dU, uC.r += dUC.r, 
	      uC.i += dUC.i)
	 {
	    uC.r = u;
	    uC.i = u * tau / wRef;
	    
	    uC2.r = 2 * uC.r;
	    uC2.i = 2 * uC.i;
	    
	    aux = uC.r * uC.r - uC.i * uC.i;
	    uuC.i = 2 * uC.r * uC.i;
	    uuC.r = aux;