Guide to programming with madagascar

From Madagascar
Jump to: navigation, search
This page was created from the LaTeX source in book/rsf/rsf/demo.tex using latex2wiki

This guide demonstrates a simple time-domain finite-differences modeling code in RSF.

Contents

Introduction

This section presents time-domain finite-difference modeling [1] written with the RSF library. The program is demonstrated with the C, C++ and Fortran 90 interfaces. The acoustic wave-equation


\Delta U - \frac{1}{v^2} \frac{\partial^2 U}{\partial t^2} = f(t)

can be written as


\left[ \Delta U - f(t) \right] v^2 =
\frac{\partial^2 U}{\partial t^2} \;.

\Delta is the Laplacian symbol, f(t) is the source wavelet, v is the velocity, and U is a scalar wavefield. A discrete time-step involves the following computations:


U_{i+1} = \left[ \Delta U -f(t) \right] v^2 \Delta t^2 + 2 U_{i} - U_{i-1} \;,

where U_{i-1}, U_{i} and U_{i+1} represent the propagating wavefield at various time steps.

C program

<syntaxhighlight lang="c"> /* time-domain acoustic FD modeling */

  1. include <rsf.h>

int main(int argc, char* argv[]) {

   /* Laplacian coefficients */
   float c0=-30./12.,c1=+16./12.,c2=- 1./12.;
   bool verb;           /* verbose flag */
   sf_file Fw=NULL,Fv=NULL,Fr=NULL,Fo=NULL; /* I/O files */
   sf_axis at,az,ax;    /* cube axes */
   int it,iz,ix;        /* index variables */
   int nt,nz,nx;
   float dt,dz,dx,idx,idz,dt2;
   float  *ww,**vv,**rr;     /* I/O arrays*/
   float **um,**uo,**up,**ud;/* tmp arrays */
   sf_init(argc,argv);
   if(! sf_getbool("verb",&verb)) verb=0;
   /* setup I/O files */
   Fw = sf_input ("in" );
   Fo = sf_output("out");
   Fv = sf_input ("vel");
   Fr = sf_input ("ref");
   /* Read/Write axes */
   at = sf_iaxa(Fw,1); nt = sf_n(at); dt = sf_d(at);
   az = sf_iaxa(Fv,1); nz = sf_n(az); dz = sf_d(az);
   ax = sf_iaxa(Fv,2); nx = sf_n(ax); dx = sf_d(ax);
   sf_oaxa(Fo,az,1); 
   sf_oaxa(Fo,ax,2); 
   sf_oaxa(Fo,at,3);
   dt2 =    dt*dt;
   idz = 1/(dz*dz);
   idx = 1/(dx*dx);
   /* read wavelet, velocity & reflectivity */
   ww=sf_floatalloc(nt);     sf_floatread(ww   ,nt   ,Fw);
   vv=sf_floatalloc2(nz,nx); sf_floatread(vv[0],nz*nx,Fv);
   rr=sf_floatalloc2(nz,nx); sf_floatread(rr[0],nz*nx,Fr);
   /* allocate temporary arrays */
   um=sf_floatalloc2(nz,nx);
   uo=sf_floatalloc2(nz,nx);
   up=sf_floatalloc2(nz,nx);
   ud=sf_floatalloc2(nz,nx);
   
   for (iz=0; iz<nz; iz++) {

for (ix=0; ix<nx; ix++) { um[ix][iz]=0; uo[ix][iz]=0; up[ix][iz]=0; ud[ix][iz]=0; }

   }
   
   /* MAIN LOOP */
   if(verb) fprintf(stderr,"\n");
   for (it=0; it<nt; it++) {

if(verb) fprintf(stderr,"\b\b\b\b\b%d",it);

/* 4th order laplacian */ for (iz=2; iz<nz-2; iz++) { for (ix=2; ix<nx-2; ix++) { ud[ix][iz] = c0* uo[ix ][iz ] * (idx+idz) + c1*(uo[ix-1][iz ] + uo[ix+1][iz ])*idx + c2*(uo[ix-2][iz ] + uo[ix+2][iz ])*idx + c1*(uo[ix ][iz-1] + uo[ix ][iz+1])*idz + c2*(uo[ix ][iz-2] + uo[ix ][iz+2])*idz; } }

/* inject wavelet */ for (iz=0; iz<nz; iz++) { for (ix=0; ix<nx; ix++) { ud[ix][iz] -= ww[it] * rr[ix][iz]; } }

/* scale by velocity */ for (iz=0; iz<nz; iz++) { for (ix=0; ix<nx; ix++) { ud[ix][iz] *= vv[ix][iz]*vv[ix][iz]; } }

/* time step */ for (iz=0; iz<nz; iz++) { for (ix=0; ix<nx; ix++) { up[ix][iz] = 2*uo[ix][iz] - um[ix][iz] + ud[ix][iz] * dt2;

um[ix][iz] = uo[ix][iz]; uo[ix][iz] = up[ix][iz]; } }

/* write wavefield to output */ sf_floatwrite(uo[0],nz*nx,Fo);

   }
   if(verb) fprintf(stderr,"\n");    
   sf_close()
   exit (0);

} </syntaxhighlight>

  • Declare input, output and auxiliary file tags: Fw for input wavelet, Fv for velocity, Fr for reflectivity, and Fo for output wavefield.

<syntaxhighlight lang="c">

   sf_file Fw,Fv,Fr,Fo; /* I/O files */

</syntaxhighlight>

  • Declare RSF cube axes: at time axis, ax space axis, az depth axis.

<syntaxhighlight lang="c">

   sf_axis at,az,ax;    /* cube axes */

</syntaxhighlight>

  • Declare multi-dimensional arrays for input, output and computations.

<syntaxhighlight lang="c">

   float  *ww,**vv,**rr;     /* I/O arrays*/

</syntaxhighlight>

  • Open files for input/output.

<syntaxhighlight lang="c">

   Fw = sf_input ("in" );
   Fo = sf_output("out");
   Fv = sf_input ("vel");
   Fr = sf_input ("ref");

</syntaxhighlight>

  • Read axes from input files; write axes to output file.

<syntaxhighlight lang="c">

   at = sf_iaxa(Fw,1); nt = sf_n(at); dt = sf_d(at);
   az = sf_iaxa(Fv,1); nz = sf_n(az); dz = sf_d(az);
   ax = sf_iaxa(Fv,2); nx = sf_n(ax); dx = sf_d(ax);
   sf_oaxa(Fo,az,1); 
   sf_oaxa(Fo,ax,2); 
   sf_oaxa(Fo,at,3);

</syntaxhighlight>

  • Allocate arrays and read wavelet, velocity and reflectivity.

<syntaxhighlight lang="c">

   ww=sf_floatalloc(nt);     sf_floatread(ww   ,nt   ,Fw);
   vv=sf_floatalloc2(nz,nx); sf_floatread(vv[0],nz*nx,Fv);
   rr=sf_floatalloc2(nz,nx); sf_floatread(rr[0],nz*nx,Fr);

</syntaxhighlight>

  • Allocate temporary arrays.

<syntaxhighlight lang="c">

   um=sf_floatalloc2(nz,nx);
   uo=sf_floatalloc2(nz,nx);
   up=sf_floatalloc2(nz,nx);
   ud=sf_floatalloc2(nz,nx);

</syntaxhighlight>

  • Loop over time.

<syntaxhighlight lang="c">

   for (it=0; it<nt; it++) {

</syntaxhighlight>

  • Compute Laplacian: \Delta U.

<syntaxhighlight lang="c"> for (iz=2; iz<nz-2; iz++) { for (ix=2; ix<nx-2; ix++) { ud[ix][iz] = c0* uo[ix ][iz ] * (idx+idz) + c1*(uo[ix-1][iz ] + uo[ix+1][iz ])*idx + c2*(uo[ix-2][iz ] + uo[ix+2][iz ])*idx + c1*(uo[ix ][iz-1] + uo[ix ][iz+1])*idz + c2*(uo[ix ][iz-2] + uo[ix ][iz+2])*idz; } } </syntaxhighlight>

  • Inject source wavelet: \left[ \Delta U - f(t) \right]

<syntaxhighlight lang="c"> for (iz=0; iz<nz; iz++) { for (ix=0; ix<nx; ix++) { ud[ix][iz] -= ww[it] * rr[ix][iz]; } } </syntaxhighlight>

  • Scale by velocity: \left[ \Delta U - f(t) \right] v^2

<syntaxhighlight lang="c"> for (iz=0; iz<nz; iz++) { for (ix=0; ix<nx; ix++) { ud[ix][iz] *= vv[ix][iz]*vv[ix][iz]; } } </syntaxhighlight>

  • Time step: U_{i+1} = \left[ \Delta U -f(t) \right] v^2 \Delta t^2 + 2 U_{i} - U_{i-1}

<syntaxhighlight lang="c"> for (iz=0; iz<nz; iz++) { for (ix=0; ix<nx; ix++) { up[ix][iz] = 2*uo[ix][iz] - um[ix][iz] + ud[ix][iz] * dt2;

um[ix][iz] = uo[ix][iz]; uo[ix][iz] = up[ix][iz]; } } </syntaxhighlight>

C++ program

<syntaxhighlight lang="cpp"> // time-domain acoustic FD modeling

  1. include <valarray>
  2. include <iostream>
  3. include <rsf.hh>
  4. include <cub.hh>
  5. include <vai.hh>

using namespace std;

int main(int argc, char* argv[]) {

   // Laplacian coefficients
   float c0=-30./12.,c1=+16./12.,c2=- 1./12.;
   sf_init(argc,argv);// init RSF
   bool verb;         // vebose flag
   if(! sf_getbool("verb",&verb)) verb=0;
   // setup I/O files
   CUB Fw("in", "i"); Fw.headin(); //Fw.report();
   CUB Fv("vel","i"); Fv.headin(); //Fv.report();
   CUB Fr("ref","i"); Fr.headin(); //Fr.report();
   CUB Fo("out","o"); Fo.setup(3,Fv.esize()); 
   // Read/Write axes
   sf_axis at = Fw.getax(0); int nt = sf_n(at); float dt = sf_d(at);
   sf_axis az = Fv.getax(0); int nz = sf_n(az); float dz = sf_d(az);
   sf_axis ax = Fv.getax(1); int nx = sf_n(ax); float dx = sf_d(ax);
   Fo.putax(0,az); 
   Fo.putax(1,ax); 
   Fo.putax(2,at);
   Fo.headou();
   float dt2 =    dt*dt;
   float idz = 1/(dz*dz);
   float idx = 1/(dx*dx);
   // read wavelet, velocity and reflectivity
   valarray<float> ww( nt    ); ww=0; Fw >> ww;
   valarray<float> vv( nz*nx ); vv=0; Fv >> vv;
   valarray<float> rr( nz*nx ); rr=0; Fr >> rr;
  
   // allocate temporary arrays
   valarray<float> um(nz*nx); um=0;
   valarray<float> uo(nz*nx); uo=0;
   valarray<float> up(nz*nx); up=0;
   valarray<float> ud(nz*nx); ud=0;
   // init ValArray Index counter
   VAI k(nz,nx);
   // MAIN LOOP
   if(verb) cerr << endl;
   for (int it=0; it<nt; it++) {

if(verb) cerr << "\b\b\b\b\b" << it;

// 4th order laplacian for (int iz=2; iz<nz-2; iz++) { for (int ix=2; ix<nx-2; ix++) { ud[k(iz,ix)] = c0* uo[ k(iz ,ix )] * (idx+idz) + c1*(uo[ k(iz ,ix-1)]+uo[ k(iz ,ix+1)]) * idx + c1*(uo[ k(iz-1,ix )]+uo[ k(iz+1,ix )]) * idz + c2*(uo[ k(iz ,ix-2)]+uo[ k(iz ,ix+2)]) * idx + c2*(uo[ k(iz-2,ix )]+uo[ k(iz+2,ix )]) * idz; } }

// inject wavelet ud -= ww[it] * rr;

// scale by velocity ud *= vv*vv;

// time step up=(float)2 * uo - um + ud * dt2; um = uo; uo = up;

// write wavefield to output output Fo << uo;

   }
   if(verb) cerr << endl;
   exit(0);

} </syntaxhighlight>

  • Declare input, output and auxiliary file cubes (of type CUB).

<syntaxhighlight lang="cpp">

   CUB Fw("in", "i"); Fw.headin(); //Fw.report();
   CUB Fv("vel","i"); Fv.headin(); //Fv.report();
   CUB Fr("ref","i"); Fr.headin(); //Fr.report();
   CUB Fo("out","o"); Fo.setup(3,Fv.esize()); 

</syntaxhighlight>

  • Declare, read and write RSF cube axes: at time axis, ax space axis, az depth axis.

<syntaxhighlight lang="cpp">

   sf_axis at = Fw.getax(0); int nt = sf_n(at); float dt = sf_d(at);
   sf_axis az = Fv.getax(0); int nz = sf_n(az); float dz = sf_d(az);
   sf_axis ax = Fv.getax(1); int nx = sf_n(ax); float dx = sf_d(ax);
   Fo.putax(0,az); 
   Fo.putax(1,ax); 
   Fo.putax(2,at);
   Fo.headou();

</syntaxhighlight>

  • Declare multi-dimensional valarrays for input, output and read data.

<syntaxhighlight lang="cpp">

   valarray<float> ww( nt    ); ww=0; Fw >> ww;
   valarray<float> vv( nz*nx ); vv=0; Fv >> vv;
   valarray<float> rr( nz*nx ); rr=0; Fr >> rr;

</syntaxhighlight>

  • Declare multi-dimensional valarrays for temporary storage.

<syntaxhighlight lang="cpp">

   valarray<float> um(nz*nx); um=0;
   valarray<float> uo(nz*nx); uo=0;
   valarray<float> up(nz*nx); up=0;
   valarray<float> ud(nz*nx); ud=0;

</syntaxhighlight>

  • Initialize multidimensional valarray index counter (of type VAI).

<syntaxhighlight lang="cpp">

   VAI k(nz,nx);

</syntaxhighlight>

  • Loop over time.

<syntaxhighlight lang="cpp">

   for (int it=0; it<nt; it++) {

</syntaxhighlight>

  • Compute Laplacian: \Delta U.

<syntaxhighlight lang="cpp"> for (int iz=2; iz<nz-2; iz++) { for (int ix=2; ix<nx-2; ix++) { ud[k(iz,ix)] = c0* uo[ k(iz ,ix )] * (idx+idz) + c1*(uo[ k(iz ,ix-1)]+uo[ k(iz ,ix+1)]) * idx + c1*(uo[ k(iz-1,ix )]+uo[ k(iz+1,ix )]) * idz + c2*(uo[ k(iz ,ix-2)]+uo[ k(iz ,ix+2)]) * idx + c2*(uo[ k(iz-2,ix )]+uo[ k(iz+2,ix )]) * idz; } } </syntaxhighlight>

  • Inject source wavelet: \left[ \Delta U - f(t) \right]

<syntaxhighlight lang="cpp"> ud -= ww[it] * rr; </syntaxhighlight>

  • Scale by velocity: \left[ \Delta U - f(t) \right] v^2

<syntaxhighlight lang="cpp"> ud *= vv*vv; </syntaxhighlight>

  • Time step: U_{i+1} = \left[ \Delta U -f(t) \right] v^2 \Delta t^2 + 2 U_{i} - U_{i-1}

<syntaxhighlight lang="cpp"> up=(float)2 * uo - um + ud * dt2; um = uo; uo = up; </syntaxhighlight>

Fortran 90 program

<syntaxhighlight lang="fortran"> ! time-domain acoustic FD modeling program AFDMf90

 use rsf
 implicit none
 ! Laplacian coefficients
 real :: c0=-30./12.,c1=+16./12.,c2=- 1./12.
 logical    :: verb         ! verbose flag
 type(file) :: Fw,Fv,Fr,Fo  ! I/O files
 type(axa)  :: at,az,ax     ! cube axes
 integer    :: it,iz,ix     ! index variables
 real       :: idx,idz,dt2
 real, allocatable :: vv(:,:),rr(:,:),ww(:)           ! I/O arrays
 real, allocatable :: um(:,:),uo(:,:),up(:,:),ud(:,:) ! tmp arrays
 call sf_init() ! init RSF
 call from_par("verb",verb,.false.)
 ! setup I/O files
 Fw=rsf_input ("in")
 Fv=rsf_input ("vel")
 Fr=rsf_input ("ref")
 Fo=rsf_output("out")
 ! Read/Write axes
 call iaxa(Fw,at,1); call iaxa(Fv,az,1); call iaxa(Fv,ax,2)
 call oaxa(Fo,az,1); call oaxa(Fo,ax,2); call oaxa(Fo,at,3)
 dt2 =    at%d*at%d
 idz = 1/(az%d*az%d)
 idx = 1/(ax%d*ax%d) 
 ! read wavelet, velocity & reflectivity
 allocate(ww(at%n     )); ww=0.; call rsf_read(Fw,ww)
 allocate(vv(az%n,ax%n)); vv=0.; call rsf_read(Fv,vv)
 allocate(rr(az%n,ax%n)); rr=0.; call rsf_read(Fr,rr)
 ! allocate temporary arrays
 allocate(um(az%n,ax%n)); um=0.
 allocate(uo(az%n,ax%n)); uo=0.
 allocate(up(az%n,ax%n)); up=0.
 allocate(ud(az%n,ax%n)); ud=0.
 ! MAIN LOOP
 do it=1,at%n
    if(verb) write (0,*) it
    ! 4th order laplacian
    do iz=2,az%n-2
       do ix=2,ax%n-2
          ud(iz,ix) = &
               c0* uo(iz,  ix  ) * (idx + idz)        + &
               c1*(uo(iz  ,ix-1) + uo(iz  ,ix+1))*idx + &
               c2*(uo(iz  ,ix-2) + uo(iz  ,ix+2))*idx + &
               c1*(uo(iz-1,ix  ) + uo(iz+1,ix  ))*idz + &
               c2*(uo(iz-2,ix  ) + uo(iz+2,ix  ))*idz
       end do
    end do
    ! inject wavelet
    ud = ud - ww(it) * rr
    ! scale by velocity
    ud= ud *vv*vv
    ! time step
    up = 2*uo - um + ud * dt2
    um =   uo
    uo =   up
    ! write wavefield to output
    call rsf_write(Fo,uo)
 end do
 call exit(0)

end program AFDMf90 </syntaxhighlight>

  • Declare input, output and auxiliary file tags.

<syntaxhighlight lang="fortran">

 type(file) :: Fw,Fv,Fr,Fo  ! I/O files

</syntaxhighlight>

  • Declare RSF cube axes: at time axis, ax space axis, az depth axis.

<syntaxhighlight lang="fortran">

 type(axa)  :: at,az,ax     ! cube axes

</syntaxhighlight>

  • Declare multi-dimensional arrays for input, output and computations.

<syntaxhighlight lang="fortran">

 real, allocatable :: vv(:,:),rr(:,:),ww(:)           ! I/O arrays
 real, allocatable :: um(:,:),uo(:,:),up(:,:),ud(:,:) ! tmp arrays

</syntaxhighlight>

  • Open files for input/output.

<syntaxhighlight lang="fortran">

 Fw=rsf_input ("in")
 Fv=rsf_input ("vel")
 Fr=rsf_input ("ref")
 Fo=rsf_output("out")

</syntaxhighlight>

  • Read axes from input files; write axes to output file.

<syntaxhighlight lang="fortran">

 call iaxa(Fw,at,1); call iaxa(Fv,az,1); call iaxa(Fv,ax,2)
 call oaxa(Fo,az,1); call oaxa(Fo,ax,2); call oaxa(Fo,at,3)

</syntaxhighlight>

  • Allocate arrays and read wavelet, velocity and reflectivity.

<syntaxhighlight lang="fortran">

 allocate(ww(at%n     )); ww=0.; call rsf_read(Fw,ww)
 allocate(vv(az%n,ax%n)); vv=0.; call rsf_read(Fv,vv)
 allocate(rr(az%n,ax%n)); rr=0.; call rsf_read(Fr,rr)

</syntaxhighlight>

  • Allocate temporary arrays.

<syntaxhighlight lang="fortran">

 allocate(um(az%n,ax%n)); um=0.
 allocate(uo(az%n,ax%n)); uo=0.
 allocate(up(az%n,ax%n)); up=0.
 allocate(ud(az%n,ax%n)); ud=0.

</syntaxhighlight>

  • Loop over time.

<syntaxhighlight lang="fortran">

 do it=1,at%n

</syntaxhighlight>

  • Compute Laplacian: \Delta U.

<syntaxhighlight lang="fortran">

    do iz=2,az%n-2
       do ix=2,ax%n-2
          ud(iz,ix) = &
               c0* uo(iz,  ix  ) * (idx + idz)        + &
               c1*(uo(iz  ,ix-1) + uo(iz  ,ix+1))*idx + &
               c2*(uo(iz  ,ix-2) + uo(iz  ,ix+2))*idx + &
               c1*(uo(iz-1,ix  ) + uo(iz+1,ix  ))*idz + &
               c2*(uo(iz-2,ix  ) + uo(iz+2,ix  ))*idz
       end do
    end do

</syntaxhighlight>

  • Inject source wavelet: \left[ \Delta U - f(t) \right]

<syntaxhighlight lang="fortran">

    ud = ud - ww(it) * rr

</syntaxhighlight>

  • Scale by velocity: \left[ \Delta U - f(t) \right] v^2

<syntaxhighlight lang="fortran">

    ud= ud *vv*vv

</syntaxhighlight>

  • Time step: U_{i+1} = \left[ \Delta U -f(t) \right] v^2 \Delta t^2 + 2 U_{i} - U_{i-1}

<syntaxhighlight lang="fortran">

    up = 2*uo - um + ud * dt2
    um =   uo
    uo =   up

</syntaxhighlight>

References

  1. "Hello world" of seismic imaging.