next up previous contents [pdf] index

Next: Data types Up: An example: Finite-Difference modeling Previous: C program

Explanation of the code

  2    #include <rsf.h>
  4    int main(int argc, char* argv[])
Line 2 is a preprocessor directive to include the rsf.h header file which contains the RSF library functions. Line 4 has parameters in the main function. This is to enable the program to take command line arguments. char* argv[] defines the pointer to the array of type char and int argc is the length of that array.
  7        float c0=-30./12.,c1=+16./12.,c2=- 1./12.;
As was mentioned earlier, the Laplacian is being evaluated with an accuracy of up to the fourth order. These coefficients arise as a result of using five terms in the discrete form of the Laplacian.
  9        bool verb;           /* verbose flag */
 10        sf_file Fw=NULL,Fv=NULL,Fr=NULL,Fo=NULL; /* I/O files */
 11        sf_axis at,az,ax;    /* cube axes */
 12        int it,iz,ix;        /* index variables */
 13        int nt,nz,nx;
 14        float dt,dz,dx,idx,idz,dt2;
Line 9 defines a variable verb of type bool. This variable will be used in the program to check for verbosity flag. Lines 10-11 define the variables of the abstract data type provided by the RSF API. These will be used to store the input and output files. Lines 12-14 are the variables of integer and float type defined to be used as running variables (it, iz, ix) for the main loop, length of the axes (nt, nz, nx), the sampling of the axes (dt, dx, dz) and the squares and inverse squares of the samples (dt2, idz, idx).
 16        float  *ww,**vv,**rr;     /* I/O arrays*/
 17        float **um,**uo,**up,**ud;/* tmp arrays */
Lines 16-17 define pointers to the arrays, which will be used for input (*ww , **vv , **rr) and for temporary storage (**um,**uo,**up,**ud).
 19	        sf_init(argc,argv);
 20         if(! sf_getbool("verb",&verb)) verb=0;
Line 19 initializes the symbol tables used to store the argument from the command line. Line 20 tests the verbosity flag specified in the command line arguments. If the verbosity flag in the command line is set to n, the variable verb (of type bool) is set to zero. This would allow the verbose output to be printed only if the user set the verbosity flag to y in the command line.
 22        /* setup I/O files */
 23   	   Fw = sf_input ("in" );
 24   	   Fo = sf_output("out");
 25   	   Fv = sf_input ("vel");
 26   	   Fr = sf_input ("ref");
In these lines we use the [sec:sf_input]sf_input (see p. [*]) and [sec:sf_output]sf_output (see p. [*]) functions of the RSF API. These functions take a string as argument and return a variable of type sf_file, we had already defined this type of variables earlier in the program.
 28         /* Read/Write axes */
 29         at = sf_iaxa(Fw,1); nt = sf_n(at); dt = sf_d(at);
 30         az = sf_iaxa(Fv,1); nz = sf_n(az); dz = sf_d(az);
 31         ax = sf_iaxa(Fv,2); nx = sf_n(ax); dx = sf_d(ax);
Here we input axes (at, az, ax) using [sec:sf_iaxa]sf_iaxa (p. [*] of the RSF API. sf_iaxa accepts a variables of type sf_file (RSF API) and an integer. The first argument in sf_iaxa is the input file and the second is the axis which we want to input. In the second column we use [sec:sf_n]sf_n (p. [*])from RSF API to get the lengths of the respective axes. In the third column we use [sec:sf_d]sf_d (p. [*]) of the RSF API to get the sampling interval of the respective axes.
 33         sf_oaxa(Fo,az,1); 
 34         sf_oaxa(Fo,ax,2); 
 35         sf_oaxa(Fo,at,3);
Here we output axes (at, az, ax) using [sec:sf_oaxa]sf_oaxa (p. [*]) of the RSF API. sf_oaxa accepts variables of type sf_file (RSF API), sf_axis (RSF API) and an integer. First argument is the output file, second argument is the name of the axis which we want to output and the third is the number of the axis in the output file (n1 is the fastest axis).
 37       dt2 =    dt*dt;
 38       idz = 1/(dz*dz);
 39       idx = 1/(dx*dx);
These lines define the square of the time sampling interval(dt2) and the inverse squares of the sampling interval of the spatial axes.
 41       /* read wavelet, velocity & reflectivity */
 42       ww = sf_floatalloc(nt);     sf_floatread(ww   ,nt   ,Fw);
 43       vv = sf_floatalloc2(nz,nx); sf_floatread(vv[0],nz*nx,Fv);
 44       rr = sf_floatalloc2(nz,nx); sf_floatread(rr[0],nz*nx,Fr);
In the first column we allocate the memory required to hold the input wavelet, velocity and reflectivity. This is done using [sec:sf_floatalloc]sf_floatalloc (p. [*]) and [sec:sf_floatalloc2]sf_floatalloc2 (p. [*]) of the RSF API. sf_floatalloc takes integers as arguments and from these integers it calculates an allocates a block of memory of appropriate size. sf_floatalloc2 is the same as sf_floatalloc except for the fact that the former allocates an array of two dimensions, size of the memory block assigned in this case is the product of the two integers given as arguments (e.g. nz*nx in this case). Then [sec:sf_floatread]sf_floatread (p. [*]) of the RSF API is used to read the data from the files into the allocated memory blocks (arrays). The sf_floatread takes the arrays, integers and files as arguments and returns arrays filled with the data from the files.
 46        /* allocate temporary arrays */
 47        um = sf_floatalloc2(nz,nx);
 48        uo = sf_floatalloc2(nz,nx);
 49        up = sf_floatalloc2(nz,nx);
 50        ud = sf_floatalloc2(nz,nx);
Just like the memory blocks were allocated for input files to be read in to, we now allocate memory for the temporary arrays which will be used just for the calculation, using [sec:sf_floatalloc2]sf_floatalloc2 (p. [*]).
 52        for (iz=0; iz<nz; iz++) {
 53            for (ix=0; ix<nx; ix++) {
 54                um[ix][iz]=0;
 55                uo[ix][iz]=0;
 56                up[ix][iz]=0;
 57                ud[ix][iz]=0;
 58            }
 59        }
Lines 52-59 initialize the temporary arrays by assigning 0 to each element of every array.
 61        /* MAIN LOOP */
 62        if(verb) fprintf(stderr,"\n");
 63        for (it=0; it<nt; it++) {
 64            if(verb) fprintf(stderr,"\b\b\b\b\b%d",it);
Now the main loop starts. The if condition in line 61 prints the message specified in the \texttt{fprintf} argument. The \texttt{stderr} is a stream in C which is used to direct the output to the screen. In this case the input is just an escape sequence \texttt{\\n}, which will bring the cursor to the next line if the user opted \texttt{y} or \texttt{1} to verbose flag in the command line (\texttt{verb=y} of \texttt{verb=1}).

Then the loop over time starts. Right after the for statement (within the body of the loop) there          is another if condition like the first one but this time it prints the the current value of \texttt{it}. This has escape sequence \texttt{\b} occurring several times. This is when the loop starts the value of \texttt{it} which is \texttt{0}, is printed on the screen, when the loop returns to the start the new value of it is \texttt{1}, so \texttt{\\b} (backspace) removes the previous value \texttt{0}, which is already on the screen, and puts \texttt{1} instead. 

\item [\bf 66-76:]
 66             /* 4th order laplacian */
 67             for (iz=2; iz<nz-2; iz++) {
 68                 for (ix=2; ix<nx-2; ix++) {
 69                     ud[ix][iz] = 
 70                       c0* uo[ix  ][iz  ] * (idx+idz) + 
 71                       c1*(uo[ix-1][iz  ] + uo[ix+1][iz  ])*idx +
 72                       c2*(uo[ix-2][iz  ] + uo[ix+2][iz  ])*idx +
 73                       c1*(uo[ix  ][iz-1] + uo[ix  ][iz+1])*idz +
 74                       c2*(uo[ix  ][iz-2] + uo[ix  ][iz+2])*idz;	  
 75                 }
 76             }
This is the calculation for the fourth order laplacian. By the term ``4th order'' we mean the order of the approximation not the order of the PDE itself which of course is a second order PDE. A second order partial derivative discretized to second order approximation is written as:

$\displaystyle \frac{\partial^2U}{\partial x^2} = \frac{U_{i+1} + 2U_i + U_{i-1}}{\Delta x^2}$    

This is the central difference formula for the second order partial derivative with pivot at $ i-$ th value of $ U$ .Similarly for the $ z$ direction we have:

$\displaystyle \frac{\partial^2U}{\partial z^2} = \frac{U_{i+1} + 2U_i + U_{i-1}}{\Delta z^2}$    

By adding these two we get the central difference formula accurate to the second order for the Laplacian. But we are using a central difference accurate up to the fourth order so for that we have:

$\displaystyle \frac{\partial^2U}{\partial x^2} = \frac{1}{\Delta x^2} \left[-\f...
...2}U_{i+1} -\frac{30}{12}U_i + \frac{16}{12}U_{i-i} - \frac{1}{12}U_{i-2}\right]$    

By writing down a similar equation for $ z$ and adding the two we get the fourth order approximation of the Laplacian or as we refer to it here ``4th order laplacian''. Now returning back to the code, the first line is the start of the loop in the $ z$ direction. Within the body of the z loop there is another loop which runs through all the values of x for one value of z. The second line is start of the for-loop for the x direction. Then in the body of the loop for x direction we use the $ 2\times2$ arrays which we defined earlier. This is just the equation of the Laplacian accurate up to the fourth order, as discussed above, with the common coefficients factored out. Note that the loops for x and z start two units after 0 and end two units before nx and nz. This is because to evaluate the Laplacian at a particular point $ (x,z)$ the farthest values which we are using are two units behind and two units ahead of the current point $ (x,z)$ if we include the points iz=0,1 ; iz=nz-1, nz and ix=0,1;ix=nx-1,nx we will run out of bounds. To fill these we will need a boundary condition which we will get from the next loop for inserting the wavelet.
 78             /* inject wavelet */
 79             for (iz=0; iz<nz; iz++) {
 80                 for (ix=0; ix<nx; ix++) {
 81                     ud[ix][iz] -= ww[it] * rr[ix][iz];
 82                 }
 83             }
These lines insert the wavelet, which means evaluating the expression $ \Delta U - f(t)$ . $ \Delta U$ was already calculated in the previous loop and is stored as the array ud. ww is the array of the wavelet but before subtracting it form the Laplacian (ud) we multiply the wavelet amplitude at current time with the reflectivity at every point in space $ (x,z)$ . This amounts to an initial condition:

$\displaystyle f(x,z,0) = g(x,z) = ww(0)rr(x,z),$    

and thus serves the purpose of filling the values at ix,iz=0 and ix-2,ix-1=0 and iz-2,iz-1=0. But the source wavelet is not an ideal impulse so it has amplitudes at future times so for each time the wavelet will be multiplied by the reflectivity at every point $ (x,z)$ . Why multiply the wavelet with reflectivity? Well, this model assumes a hypothetical situation that the source was set off at each and every point in space $ (x,z)$ under consideration and scaled by the reflectivity at that point $ (x,z)$ . What this means is that the source was set off at all the points where there is a change in the acoustic impedance (because reflectivity is the ratio of the difference and sum of the acoustic impedances across an interface).
 85             /* scale by velocity */
 86             for (iz=0; iz<nz; iz++) {
 87                 for (ix=0; ix<nx; ix++) {
 88                     ud[ix][iz] *= vv[ix][iz]*vv[ix][iz];
 89                 }
 90             }
Here we just multiply $ \Delta U - f(t)$ by the velocity, that is, we evaluate $ (\Delta U - f(t))v^2$
 92             /* time step */
 93             for (iz=0; iz<nz; iz++) {
 94                 for (ix=0; ix<nx; ix++) {
 95                     up[ix][iz] = 2*uo[ix][iz] 
 96                                  - um[ix][iz] 
 97                                  + ud[ix][iz] * dt2; 
 99                     um[ix][iz] = uo[ix][iz];
100                     uo[ix][iz] = up[ix][iz];
101                   }
102               }
Here we calculate the time step, that is,

$\displaystyle U_{i+1} = [\Delta U - f(t)]v^2\Delta t^2 + 2U_i - U_{i-1}.$    

The first for-loop is for the $ z$ direction and within the body of this loop is another for-loop for the $ x$ direction. up is the array which holds the amplitude of the wave at the current time in the time loop. uo is the array which contains the amplitude at a time one unit before the current time and the array um holds the amplitude two units before. ud is the array we calculated earlier in the program, now it gets multiplied by $ \Delta t^2$ (dt2) and included in the final equation. This completes the calculation for one value of it. Now the arrays need to be updated to represent the next time step. This is done in the last two: The first one says $ U_{i-1} \rightarrow U_i $ and the second one says $ U_i \rightarrow U_{i+1}$ , that is, the array um is updated by uo and then the array uo itself gets updated by up.
104               /* write wavefield to output */
105               sf_floatwrite(uo[0],nz*nx,Fo);
106         }
After the calculations for one time step are complete we write the array uo (remember that uo was made equal to up, which is the current time step, in the previous line). To write the array in the output file we use [sec:sf_floatwrite]sf_floatwrite (p. [*]) exactly the same way we used sf_floatread to read in from the input files, only difference is that the array given as the argument is written into the file given in the last argument. The bracket close is for the time loop, after this the time loop will start all over again for the next time value.
107         if(verb) fprintf(stderr,"\n");    
108         sf_close()
109         exit(0);
110     }
The first line puts the cursor in the new line on the screen after the time loop has run through all the time values. The second line uses [sec:sf_close]sf_close (p. [*]) from RSF API to remove the temporary files. The third line uses the exit() function in C language to close the streams and return the control to the host environment. The 0 in the argument indicates normal termination of the program. The last bracket closes the main function.

next up previous contents [pdf] index

Next: Data types Up: An example: Finite-Difference modeling Previous: C program