sudrefsvmaster.c

su 的源代码库

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

#include "refMaster.h"

segy trR, trZ;     /* SEGY trace */
/*********************** self documentation **********************/
char *sdoc[] = {
"                                                                        ",
" SUDREFSVMASTER - Distributed REFlectivity modeling MASTER code for     ",
"                  the generation of the vertical and radial component   ",
"                  (SV) of the displacement field due to a point         ",
"                  force. In this code a master-slave algorithm is       ",
"                  implemented. A master processor sends several         ",
"                  frequency partitions to different processors. After   ",
"                  modeling its frequency range the slave processor      ",
"                  sends the result back to the master. The frequency    ",
"                  partitions that are still to be modeled, if any, are  ",
"                  sent by the master to the available slaves. After     ",
"                  composition all the frequency partitions the master   ",
"                  outputs the modeled time seismograms to stdout.       ",
"                                                                        ",
"                  The code SUDREFSVSLAVE implements the slave functions ",
"                  and should be installed with SUDREFSVMASTER. The      ",
"                  communication and synchronization routines are        ",
"                  provided by the PVM library. The version used by      ",
"                  SUDREFSVMASTER is the version 3.3.9 or higher. The user ",
"                  MUST install PVM prior to the installation of this    ",
"                  code.                                                 ",
"                                                                        ",
" Many of the input parameters will be better understood if the reader   ",
" is familiar with the following two references:                         ",
"                                                                        ",
" Muller, G: The Reflectivity Method: a Tutorial, Journal of Geophysics  ",
"            (1985) 58: 153-174                                          ",
"                                                                        ",
" Mallick S. and Fraser N.: Practical Aspects of Reflectivity Modeling   ",
"            Geophysics (1987) 52 no 10: 1355-1364                       ",
"                                                                        ",
"                                                                        ",
" SUDREFSVMAIN requires the following input parameters (Values between   ",
" brackets represent the default value of the parameter).                ",
"                                                                        ",
" model=[model]  Elastic mode ASCII file. This file should contain one   ",
"                row per layer of the model in the following format:     ",
"                                                                        ",
"                Thick(km)  Rho(g/cm^3)  Vp(km/s)  Qp  Vs(km/s)  Qs      ",
"                                                                        ",
"                Where:                                                  ",
"                Thick: Layer thickness                                  ",
"                Rho:   Layer density                                    ",
"                Vp:    P-wave velocity                                  ",
"                Qp:    P-wave quality factor                            ",
"                Vs:    S-wave velocity                                  ",
"                Qs:    S-wave quality factor                            ",
" recfile=[NULL] File with receiver offsets. The format is one offset    ",
"                line. If none is specified (default), the acquisition   ",
"                is end on                                               ",
" directwave=1   Direct wave flag. If zero direct wave is not generated  ",
" r1=[0]         Minimum source-receiver offset (km)                     ",
" nr=[48]        Number of receivers                                     ",
" dr=[0.025]     Receiver spacing (km)                                   ",
" zs=[0.001]     Source depth (km)                                       ",
" u1=[0.0002]    Initial horizontal slowness (s/km)                      ",
" u2=[1]         Final horizontal slowness (s/km)                        ",
" nu=[1000]      Number of slownesses                                    ",
" f1=[2]         Initial frequency for modeling (Hz)                     ",
" f2=[50]        Final frequency for modeling (Hz)                       ",
" dt=[0.004]     Time sampling interval (0.004)                          ",
" tmax=[8]       Maximum modeling time (s)                               ",
" F1=[0]         EAST-WEST point force coordinate                        ",
" F2=[0]         NORTH-SOUTH point force coordinate                      ",
" F3=[1]         VERTICAL point force coordinate                         ",
" Hanning=[1]    Final trace is convolved with a hanning window          ",
" wu=[5]         Percentual of slowness interval that is windowed        ",
" ww=[5]         Percentual of frequency interval that is windowed       ",
" fr=[1]         Reference frequency used in absorption model            ",
" tau=[50]       Wrap-around attenuation parameter                       ",
" nProc=[1]      Number of processors used in the computation            ",
" nfreqproc=[0]  Number of frequencies allocated to each processor       ",
"                The default value means that the frequency              ",
"                components will be evenly split among the slave         ",
"                processors                                              ",
" verbose=[0]    =1 print useful information                             ",
"                                                                        ",
" SUDREFSVMASTER output the synthetic data in SU format to standard      ",
" output. First the vertical then the horizontal component.              ",
"                                                                        ",
" Reminder: YOU SHOULD INSTALL PVM (version 3.3.9 or higher) BEFORE THE  ",
"           INSTALLATION OF THESE CODES                                  ",
"                                                                        ",
" Notes on execution:                                                    ",
"                                                                        ",
" You should start PVM before the execution of SUDREFSVMASTER. Each slave",
" will generate a report file names PSU[XXX], where [XXX] is the process ",
" id number. Those files are placed in the directory specified by the    ",
" environment variable PVMLOG.                                           ",
"                                                                        ",
" Limitations of this version:                                           ",
"                                                                        ",
" 1) No SH displacement field is generated                               ",
"                                                                        ",
" 2) Point source (explosion like) is not implemented                    ",
"                                                                        ",
"                                 Author: Wences Gouveia - 95-11-29      ",
"                                 Center For Wave Phenomena              ",
"                                 Colorado School of Mines               ",
NULL};   
/************************ end self doc ***********************************/ 

void
main (int argc, char **argv)
{
   /* declaration of variables */
   FILE *fp;                     /* file pointer */
   char *auxChar;                /* auxiliar character */
   char *modelFile = " ";        /* elastic model file */
                                 /* THICK - RHO - VP - QP - VS - QS */
   int i, k, iProc, iR;          /* counters */
   int initF, lastF;             /* initial and final frequencies */
   int apl_pid;                  /* PVM process id control */
   int nSamplesOrig;             /* time series length */
   int die;                      /* flag used to kill processes */
   int pid;                      /* process id */
   int nProc;                    /* number of processes */
   int processControl;           /* monitoring PVM start */
   int *processes;               /* array with process ids */
   int FReceived;                /* number of frequencies processed */
   int nFreqProc;                /* number of frequencies per process */
   int nFreqPart;                /* number of frequency partitions */
   int **statusFreq;             /* monitors processed frequencies */
   int FInfo[2];                 /* frequency delimiters */
   int **procInfo;               /* frequency limits for each processor */ 
   float wallcpu;                /* wall clock time */
   float dt;                     /* time sampling interval */
   float f;                      /* current frequency */
   float fR;                     /* reference frequency */
   float tMax;                   /* maximum recording time */
   float *thick, *alpha, *beta,
   *rho, *qP, *qS;               /* elastic constants and thickness */
   complex **freqPart;           /* frequency arrays sent by the slaves */
   complex **uRF, **uZF;         /* final frequency components */
   INFO info[1];                 /* basic information for slaves */
   
   /* Logging information */
   CleanLog();

   /* getting input */
   initargs(argc, argv);
   requestdoc(0);
   
   if (!getparstring("model", &modelFile)) modelFile = "model";
   if (!getparstring("recfile", &auxChar)) auxChar = " ";
   sprintf(info->recFile, "%s", auxChar);
   if (!getparint("directwave", &info->directWave)) info->directWave = 1;
   if (!getparfloat("r1", &info->r1)) info->r1 = 0;
   if (!getparint("nr", &info->nR)) info->nR = 148;
   if (!getparfloat("dr", &info->dR)) info->dR = .025;
   if (!getparfloat("zs", &info->zs)) info->zs = 0.001;
   if (info->zs <= 0) info->zs = 0.001;
   if (!getparfloat("u1", &info->u1)) info->u1 = 0.0002;
   if (!getparfloat("u2", &info->u2)) info->u2 = 1.;
   if (!getparint("nu", &info->nU)) info->nU = 1000;
   if (!getparfloat("f1", &info->f1)) info->f1 = 2;
   if (!getparfloat("f2", &info->f2)) info->f2 = 50;
   if (!getparfloat("dt", &dt)) dt = 0.004;
   if (!getparfloat("tmax", &tMax)) tMax = 8;
   if (!getparfloat("F1", &info->F1)) info->F1 = 0;
   if (!getparfloat("F2", &info->F2)) info->F2 = 0;
   if (!getparfloat("F3", &info->F3)) info->F3 = 1;
   if (!getparint("hanning", &info->hanningFlag)) info->hanningFlag = 0;
   if (!getparfloat("wu", &info->percU)) info->percU = 5; info->percU /= 100;
   if (!getparfloat("ww", &info->percW)) info->percW = 5; info->percW /= 100;
   if (!getparfloat("fr", &fR)) fR = 1; info->wR = 2 * PI * fR;
   if (!getparfloat("tau", &info->tau)) info->tau = 50;
   if (!getparint("nproc", &nProc)) nProc = 1;
   if (!getparint("nfreqproc", &nFreqProc) || nProc == 1) nFreqProc = 0;
   if (!getparint("verbose", &info->verbose)) info->verbose = 0;

   /* how many layers */
   fp = fopen(modelFile,"r");
   if (fp == NULL)
      err("No model file!\n");

   info->nL = 0;
   while (fscanf(fp, "%f %f %f %f %f %f\n", 
		 &f, &f, &f, &f, &f, &f) != EOF)
      info->nL++;
   info->nL--;
   fclose(fp);

   if (info->verbose)
      fprintf(stderr,"Number of layers in model %s : %d\n", 
	      modelFile, info->nL + 1); 
   
   /* if specific geometry, count number of receivers */
   fp = fopen(info->recFile, "r");
   if (fp != NULL)
   {
      info->nR = 0;
      while (fscanf(fp, "%f\n", &f) != EOF)
	 info->nR++;
   fclose(fp);
   }

   /* memory allocation */
   alpha = alloc1float(info->nL + 1);
   beta = alloc1float(info->nL + 1);
   rho = alloc1float(info->nL + 1);
   qP = alloc1float(info->nL + 1);
   qS = alloc1float(info->nL + 1);
   thick = alloc1float(info->nL + 1);
   processes = alloc1int(nProc);
   procInfo = alloc2int(2, nProc);

   /* reading the file */
   fp = fopen(modelFile,"r");
   if (info->verbose)
      fprintf(stderr,"Thickness     rho     vP     qP    vS     qS\n");
   for (i = 0; i < info->nL + 1; i++)
   {
      fscanf(fp, "%f %f %f %f %f %f\n", &thick[i], &rho[i], &alpha[i], 
	     &qP[i], &beta[i], &qS[i]);
      if (info->verbose)
	 fprintf(stderr,"   %7.4f      %4.3f   %3.2f  %5.1f  %3.2f  %5.1f\n",
		 thick[i], rho[i], alpha[i], qP[i], beta[i], qS[i]);
   }
   fclose(fp);

   /* computing frequency interval */
   info->nSamples = NINT(tMax / dt) + 1;
   nSamplesOrig = info->nSamples;
   info->nSamples = npfar(info->nSamples);

   /* slowness increment */
   info->dU = (info->u2 - info->u1) / (float) info->nU;

   /* computing more frequency related quatities */
   tMax = dt * (info->nSamples - 1);
   info->dF = 1. / (tMax);   
   f = info->dF;
   while (f < info->f1) f += info->dF;
   info->f1 = f;
   while (f < info->f2) f += info->dF;
   info->f2 = f; 
   initF = NINT(info->f1 / info->dF);
   lastF = NINT(info->f2 / info->dF);

   info->nF = NINT(info->f2 / info->dF) - NINT(info->f1 / info->dF) + 1;

   if (info->nF%2 == 0)
   {
      info->f2 += info->dF;
      info->nF++;