Difference between revisions of "Guide to programming with madagascar"

From Madagascar
Jump to navigation Jump to search
(update + python)
 
(2 intermediate revisions by one other user not shown)
Line 33: Line 33:
  
 
[[Image:wavec.png|frame|center|Wave propagation snapshot.]]
 
[[Image:wavec.png|frame|center|Wave propagation snapshot.]]
 +
 +
==C program==
  
 
<syntaxhighlight lang="c">
 
<syntaxhighlight lang="c">
Line 95: Line 97:
  
 
     /* MAIN LOOP */
 
     /* MAIN LOOP */
     if(verb) fprintf(stderr,"\n");
+
     if(verb) fprintf(stderr,"
 +
");
 
     for (it=0; it<nt; it++) {
 
     for (it=0; it<nt; it++) {
if(verb) fprintf(stderr,"\b\b\b\b\b%d",it);
+
if(verb) fprintf(stderr,"\b\b\b\b\b %d",it);
 
 
 
/* 4th order laplacian */
 
/* 4th order laplacian */
Line 147: Line 150:
 
</syntaxhighlight>
 
</syntaxhighlight>
  
#Declare input, output and auxiliary file tags: <tt>Fw</tt> for input wavelet,  <tt>Fv</tt> for velocity, <tt>Fr</tt> for reflectivity, and <tt>Fo</tt> for output wavefield.  
+
 
<syntaxhighlight lang="c">
+
#Declare input, output and auxiliary file tags: <tt>Fw</tt> for input wavelet,  <tt>Fv</tt> for velocity, <tt>Fr</tt> for reflectivity, and <tt>Fo</tt> for output wavefield. <syntaxhighlight lang="c">
 
     sf_file Fw=NULL,Fv=NULL,Fr=NULL,Fo=NULL; /* I/O files */
 
     sf_file Fw=NULL,Fv=NULL,Fr=NULL,Fo=NULL; /* I/O files */
</syntaxhighlight>  
+
</syntaxhighlight>  
#Declare RSF cube axes: <tt>at</tt> time axis, <tt>ax</tt> space axis, <tt>az</tt> depth axis.  
+
#Declare RSF cube axes: <tt>at</tt> time axis, <tt>ax</tt> space axis, <tt>az</tt> depth axis. \tiny <syntaxhighlight lang="c">
<syntaxhighlight lang="c">
 
 
     sf_axis at,az,ax;    /* cube axes */
 
     sf_axis at,az,ax;    /* cube axes */
 
</syntaxhighlight>   
 
</syntaxhighlight>   
#Declare multi-dimensional arrays for input, output and computations. <syntaxhighlight lang="c">
+
#Declare multi-dimensional arrays for input, output and computations. <syntaxhighlight lang="c">
 +
 
 
     float  *ww,**vv,**rr;    /* I/O arrays*/
 
     float  *ww,**vv,**rr;    /* I/O arrays*/
 
</syntaxhighlight>   
 
</syntaxhighlight>   
#Open files for input/output.  
+
#Open files for input/output. <syntaxhighlight lang="c">
<syntaxhighlight lang="c">
 
 
     Fw = sf_input ("in" );
 
     Fw = sf_input ("in" );
 
     Fo = sf_output("out");
 
     Fo = sf_output("out");
 
     Fv = sf_input ("vel");
 
     Fv = sf_input ("vel");
 
     Fr = sf_input ("ref");
 
     Fr = sf_input ("ref");
</syntaxhighlight>
+
</syntaxhighlight>  
#Read axes from input files; write axes to output file. \tiny <syntaxhighlight lang="c">
+
#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);
 
     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);
 
     az = sf_iaxa(Fv,1); nz = sf_n(az); dz = sf_d(az);
Line 173: Line 175:
 
     sf_oaxa(Fo,ax,2);  
 
     sf_oaxa(Fo,ax,2);  
 
     sf_oaxa(Fo,at,3);
 
     sf_oaxa(Fo,at,3);
</syntaxhighlight>
+
</syntaxhighlight>  
#Allocate arrays and read wavelet, velocity and reflectivity. \tiny <syntaxhighlight lang="c">
+
#Allocate arrays and read wavelet, velocity and reflectivity. <syntaxhighlight lang="c">
 
     ww=sf_floatalloc(nt);    sf_floatread(ww  ,nt  ,Fw);
 
     ww=sf_floatalloc(nt);    sf_floatread(ww  ,nt  ,Fw);
 
     vv=sf_floatalloc2(nz,nx); sf_floatread(vv[0],nz*nx,Fv);
 
     vv=sf_floatalloc2(nz,nx); sf_floatread(vv[0],nz*nx,Fv);
 
     rr=sf_floatalloc2(nz,nx); sf_floatread(rr[0],nz*nx,Fr);
 
     rr=sf_floatalloc2(nz,nx); sf_floatread(rr[0],nz*nx,Fr);
</syntaxhighlight>  
+
</syntaxhighlight>  
#Allocate temporary arrays.  
+
#Allocate temporary arrays. <syntaxhighlight lang="c">
<syntaxhighlight lang="c">
 
 
     um=sf_floatalloc2(nz,nx);
 
     um=sf_floatalloc2(nz,nx);
 
     uo=sf_floatalloc2(nz,nx);
 
     uo=sf_floatalloc2(nz,nx);
Line 186: Line 187:
 
     ud=sf_floatalloc2(nz,nx);
 
     ud=sf_floatalloc2(nz,nx);
 
</syntaxhighlight>   
 
</syntaxhighlight>   
#Loop over time.  
+
#Loop over time. <syntaxhighlight lang="c">
<syntaxhighlight lang="c">
 
 
     for (it=0; it<nt; it++) {
 
     for (it=0; it<nt; it++) {
</syntaxhighlight>  
+
</syntaxhighlight>
#Compute Laplacian: <math>\Delta U</math>.  
+
#Compute Laplacian: <math>\Delta U</math>. <syntaxhighlight lang="c">
<syntaxhighlight lang="c">
 
 
for (ix=2; ix<nx-2; ix++) {
 
for (ix=2; ix<nx-2; ix++) {
 
    for (iz=2; iz<nz-2; iz++) {
 
    for (iz=2; iz<nz-2; iz++) {
Line 217: Line 216:
 
}
 
}
 
</syntaxhighlight>   
 
</syntaxhighlight>   
#Time step:  <math>U_{i+1} = \left[ \Delta U -f(t) \right] v^2 \Delta t^2 + 2 U_{i} - U_{i-1}</math>  
+
#Time step:  <math>U_{i+1} = \left[ \Delta U -f(t) \right] v^2 \Delta t^2 + 2 U_{i} - U_{i-1}</math> <syntaxhighlight lang="c">
<syntaxhighlight lang="c">
 
 
for (ix=0; ix<nx; ix++) {
 
for (ix=0; ix<nx; ix++) {
 
    for (iz=0; iz<nz; iz++) {
 
    for (iz=0; iz<nz; iz++) {
Line 230: Line 228:
 
    }
 
    }
 
}
 
}
</syntaxhighlight>
+
</syntaxhighlight>  
 +
 
  
 +
\newpage
 
==C++ program==
 
==C++ program==
  
Line 240: Line 240:
 
#include <rsf.hh>
 
#include <rsf.hh>
 
#include <cub.hh>
 
#include <cub.hh>
#include <vai.hh>
+
 
 +
#include "vai.hh"
 +
 
 
using namespace std;
 
using namespace std;
  
Line 256: Line 258:
 
     CUB Fv("vel","i"); Fv.headin(); //Fv.report();
 
     CUB Fv("vel","i"); Fv.headin(); //Fv.report();
 
     CUB Fr("ref","i"); Fr.headin(); //Fr.report();
 
     CUB Fr("ref","i"); Fr.headin(); //Fr.report();
     CUB Fo("out","o"); Fo.setup(3,Fv.esize());  
+
     CUB Fo("out","o"); Fo.setup(3);  
  
 
     // Read/Write axes
 
     // Read/Write axes
Line 292: Line 294:
  
 
// 4th order laplacian
 
// 4th order laplacian
for (int iz=2; iz<nz-2; iz++) {
+
for (int ix=2; ix<nx-2; ix++) {
    for (int ix=2; ix<nx-2; ix++) {
+
    for (int iz=2; iz<nz-2; iz++) {
 
ud[k(iz,ix)] =  
 
ud[k(iz,ix)] =  
 
    c0* uo[ k(iz  ,ix  )] * (idx+idz) +
 
    c0* uo[ k(iz  ,ix  )] * (idx+idz) +
Line 321: Line 323:
 
     exit(0);
 
     exit(0);
 
}
 
}
 +
 
</syntaxhighlight>
 
</syntaxhighlight>
  
*Declare input, output and auxiliary file cubes  (of type <tt>CUB</tt>).
+
#Declare input, output and auxiliary file cubes  (of type <tt>CUB</tt>). <syntaxhighlight lang="cpp">
<syntaxhighlight lang="cpp">
 
 
     CUB Fw("in", "i"); Fw.headin(); //Fw.report();
 
     CUB Fw("in", "i"); Fw.headin(); //Fw.report();
 
     CUB Fv("vel","i"); Fv.headin(); //Fv.report();
 
     CUB Fv("vel","i"); Fv.headin(); //Fv.report();
 
     CUB Fr("ref","i"); Fr.headin(); //Fr.report();
 
     CUB Fr("ref","i"); Fr.headin(); //Fr.report();
     CUB Fo("out","o"); Fo.setup(3,Fv.esize());  
+
     CUB Fo("out","o"); Fo.setup(3);  
</syntaxhighlight>  
+
</syntaxhighlight>
*Declare, read and write RSF cube axes: <tt>at</tt> time axis, <tt>ax</tt> space axis, <tt>az</tt> depth axis.   
+
#Declare, read and write RSF cube axes: <tt>at</tt> time axis, <tt>ax</tt> space axis, <tt>az</tt> depth axis.  <syntaxhighlight lang="cpp">
<syntaxhighlight lang="cpp">
 
 
     sf_axis at = Fw.getax(0); int nt = sf_n(at); float dt = sf_d(at);
 
     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 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);
 
     sf_axis ax = Fv.getax(1); int nx = sf_n(ax); float dx = sf_d(ax);
 
+
</syntaxhighlight>  
    Fo.putax(0,az);
+
#Declare multi-dimensional  <tt>valarrays</tt> for input, output and read data. <syntaxhighlight lang="cpp">
    Fo.putax(1,ax);
 
    Fo.putax(2,at);
 
    Fo.headou();
 
</syntaxhighlight>
 
*Declare multi-dimensional  <tt>valarrays</tt> for input, output and read data.
 
<syntaxhighlight lang="cpp">
 
 
     valarray<float> ww( nt    ); ww=0; Fw >> ww;
 
     valarray<float> ww( nt    ); ww=0; Fw >> ww;
 
     valarray<float> vv( nz*nx ); vv=0; Fv >> vv;
 
     valarray<float> vv( nz*nx ); vv=0; Fv >> vv;
 
     valarray<float> rr( nz*nx ); rr=0; Fr >> rr;
 
     valarray<float> rr( nz*nx ); rr=0; Fr >> rr;
</syntaxhighlight>  
+
</syntaxhighlight>  
*Declare multi-dimensional  <tt>valarrays</tt> for temporary storage.
+
#Declare multi-dimensional  <tt>valarrays</tt> for temporary storage. <syntaxhighlight lang="cpp">
<syntaxhighlight lang="cpp">
 
 
     valarray<float> um(nz*nx); um=0;
 
     valarray<float> um(nz*nx); um=0;
 
     valarray<float> uo(nz*nx); uo=0;
 
     valarray<float> uo(nz*nx); uo=0;
Line 354: Line 348:
 
     valarray<float> ud(nz*nx); ud=0;
 
     valarray<float> ud(nz*nx); ud=0;
 
</syntaxhighlight>   
 
</syntaxhighlight>   
*Initialize multidimensional <tt>valarray</tt>  index counter (of type <tt>VAI</tt>).
+
#Initialize multidimensional <tt>valarray</tt>  index counter (of type <tt>VAI</tt>). <syntaxhighlight lang="cpp">
<syntaxhighlight lang="cpp">
 
 
     VAI k(nz,nx);
 
     VAI k(nz,nx);
 
</syntaxhighlight>   
 
</syntaxhighlight>   
*Loop over time.
+
#Loop over time. <syntaxhighlight lang="cpp">
<syntaxhighlight lang="cpp">
 
 
     for (int it=0; it<nt; it++) {
 
     for (int it=0; it<nt; it++) {
</syntaxhighlight>  
+
</syntaxhighlight>  
*Compute Laplacian: <math>\Delta U</math>.  
+
#Compute Laplacian: <math>\Delta U</math>. <syntaxhighlight lang="cpp">
<syntaxhighlight lang="cpp">
+
for (int ix=2; ix<nx-2; ix++) {
for (int iz=2; iz<nz-2; iz++) {
+
    for (int iz=2; iz<nz-2; iz++) {
    for (int ix=2; ix<nx-2; ix++) {
 
 
ud[k(iz,ix)] =  
 
ud[k(iz,ix)] =  
 
    c0* uo[ k(iz  ,ix  )] * (idx+idz) +
 
    c0* uo[ k(iz  ,ix  )] * (idx+idz) +
Line 374: Line 365:
 
    }
 
    }
 
}
 
}
</syntaxhighlight>
+
</syntaxhighlight>  
*Inject source wavelet: <math>\left[ \Delta U - f(t) \right]</math>  
+
#Inject source wavelet: <math>\left[ \Delta U - f(t) \right]</math> <syntaxhighlight lang="cpp">
<syntaxhighlight lang="cpp">
 
 
ud -= ww[it] * rr;
 
ud -= ww[it] * rr;
</syntaxhighlight>
+
</syntaxhighlight>  
*Scale by velocity: <math>\left[ \Delta U - f(t) \right] v^2</math>  
+
#Scale by velocity: <math>\left[ \Delta U - f(t) \right] v^2</math> <syntaxhighlight lang="cpp">
<syntaxhighlight lang="cpp">
 
 
ud *= vv*vv;
 
ud *= vv*vv;
</syntaxhighlight>  
+
</syntaxhighlight>
*Time step:  <math>U_{i+1} = \left[ \Delta U -f(t) \right] v^2 \Delta t^2 + 2 U_{i} - U_{i-1}</math>  
+
#Time step:  <math>U_{i+1} = \left[ \Delta U -f(t) \right] v^2 \Delta t^2 + 2 U_{i} - U_{i-1}</math> <syntaxhighlight lang="cpp">
<syntaxhighlight lang="cpp">
 
 
up=(float)2 * uo - um + ud * dt2;
 
up=(float)2 * uo - um + ud * dt2;
 
um =  uo;
 
um =  uo;
 
uo =  up;
 
uo =  up;
</syntaxhighlight>
+
</syntaxhighlight>
  
 
==Fortran 90 program==
 
==Fortran 90 program==
Line 406: Line 394:
 
   type(axa)  :: at,az,ax    ! cube axes
 
   type(axa)  :: at,az,ax    ! cube axes
 
   integer    :: it,iz,ix    ! index variables
 
   integer    :: it,iz,ix    ! index variables
 +
  integer    :: nt,nz,nx
 +
  real      :: dt,dz,dx
 
   real      :: idx,idz,dt2
 
   real      :: idx,idz,dt2
  
Line 421: Line 411:
  
 
   ! Read/Write axes
 
   ! Read/Write axes
   call iaxa(Fw,at,1); call iaxa(Fv,az,1); call iaxa(Fv,ax,2)
+
   call iaxa(Fw,at,1); nt = at%n; dt = at%d
   call oaxa(Fo,az,1); call oaxa(Fo,ax,2); call oaxa(Fo,at,3)
+
  call iaxa(Fv,az,1); nz = az%n; dz = az%d
 +
  call iaxa(Fv,ax,2); nx = ax%n; dx = ax%d
 +
 
 +
   call oaxa(Fo,az,1)
 +
  call oaxa(Fo,ax,2)
 +
  call oaxa(Fo,at,3)
  
   dt2 =    at%d*at%d
+
   dt2 =    dt*dt
   idz = 1/(az%d*az%d)
+
   idz = 1/(dz*dz)
   idx = 1/(ax%d*ax%d)  
+
   idx = 1/(dx*dx)  
  
 
   ! read wavelet, velocity & reflectivity
 
   ! read wavelet, velocity & reflectivity
   allocate(ww(at%n    )); ww=0.; call rsf_read(Fw,ww)
+
   allocate(ww(nt));   call rsf_read(Fw,ww)
   allocate(vv(az%n,ax%n)); vv=0.; call rsf_read(Fv,vv)
+
   allocate(vv(nz,nx)); call rsf_read(Fv,vv)
   allocate(rr(az%n,ax%n)); rr=0.; call rsf_read(Fr,rr)
+
   allocate(rr(nz,nx)); call rsf_read(Fr,rr)
  
 
   ! allocate temporary arrays
 
   ! allocate temporary arrays
   allocate(um(az%n,ax%n)); um=0.
+
   allocate(um(nz,nx)); um=0.
   allocate(uo(az%n,ax%n)); uo=0.
+
   allocate(uo(nz,nx)); uo=0.
   allocate(up(az%n,ax%n)); up=0.
+
   allocate(up(nz,nx)); up=0.
   allocate(ud(az%n,ax%n)); ud=0.
+
   allocate(ud(nz,nx)); ud=0.
  
 
   ! MAIN LOOP
 
   ! MAIN LOOP
   do it=1,at%n
+
   do it=1,nt
 
     if(verb) write (0,*) it
 
     if(verb) write (0,*) it
  
     ! 4th order laplacian
+
     ud(3:nz-2,3:nx-2) = &
    do iz=2,az%n-2
+
          c0* uo(3:nz-2,3:nx-2) * (idx + idz)           + &
        do ix=2,ax%n-2
+
          c1*(uo(3:nz-2,2:nx-3) + uo(3:nz-2,4:nx-1))*idx + &
          ud(iz,ix) = &
+
          c2*(uo(3:nz-2,1:nx-4) + uo(3:nz-2,5:nx ))*idx + &
                c0* uo(iz, ix  ) * (idx + idz)       + &
+
          c1*(uo(2:nz-3,3:nx-2) + uo(4:nz-1,3:nx-2))*idz + &
                c1*(uo(iz  ,ix-1) + uo(iz  ,ix+1))*idx + &
+
          c2*(uo(1:nz-4,3:nx-2) + uo(5:nz  ,3:nx-2))*idz
                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
 
     ! inject wavelet
Line 473: Line 463:
 
end program AFDMf90
 
end program AFDMf90
 
</syntaxhighlight>
 
</syntaxhighlight>
 
+
 
*Declare input, output and auxiliary file tags.
+
#Declare input, output and auxiliary file tags. <syntaxhighlight lang="fortran">
<syntaxhighlight lang="fortran">
 
 
   type(file) :: Fw,Fv,Fr,Fo  ! I/O files
 
   type(file) :: Fw,Fv,Fr,Fo  ! I/O files
 
</syntaxhighlight>   
 
</syntaxhighlight>   
*Declare RSF cube axes: <tt>at</tt> time axis, <tt>ax</tt> space axis, <tt>az</tt> depth axis.
+
#Declare RSF cube axes: <tt>at</tt> time axis, <tt>ax</tt> space axis, <tt>az</tt> depth axis. <syntaxhighlight lang="fortran">
<syntaxhighlight lang="fortran">
 
 
   type(axa)  :: at,az,ax    ! cube axes
 
   type(axa)  :: at,az,ax    ! cube axes
</syntaxhighlight>  
+
</syntaxhighlight>  
*Declare multi-dimensional arrays for input, output and computations.
+
#Declare multi-dimensional arrays for input, output and computations. <syntaxhighlight lang="fortran">
<syntaxhighlight lang="fortran">
 
 
   real, allocatable :: vv(:,:),rr(:,:),ww(:)          ! I/O arrays
 
   real, allocatable :: vv(:,:),rr(:,:),ww(:)          ! I/O arrays
 
   real, allocatable :: um(:,:),uo(:,:),up(:,:),ud(:,:) ! tmp arrays
 
   real, allocatable :: um(:,:),uo(:,:),up(:,:),ud(:,:) ! tmp arrays
</syntaxhighlight>  
+
</syntaxhighlight>
*Open files for input/output.
+
#Open files for input/output. <syntaxhighlight lang="fortran">
<syntaxhighlight lang="fortran">
 
 
   Fw=rsf_input ("in")
 
   Fw=rsf_input ("in")
 
   Fv=rsf_input ("vel")
 
   Fv=rsf_input ("vel")
 
   Fr=rsf_input ("ref")
 
   Fr=rsf_input ("ref")
 
   Fo=rsf_output("out")
 
   Fo=rsf_output("out")
 +
</syntaxhighlight> 
 +
#Read axes from input files; write axes to output file. <syntaxhighlight lang="fortran">
 +
  call iaxa(Fw,at,1); nt = at%n; dt = at%d
 +
  call iaxa(Fv,az,1); nz = az%n; dz = az%d
 +
  call iaxa(Fv,ax,2); nx = ax%n; dx = ax%d
 +
 
 +
  call oaxa(Fo,az,1)
 +
  call oaxa(Fo,ax,2)
 +
  call oaxa(Fo,at,3)
 
</syntaxhighlight>   
 
</syntaxhighlight>   
*Read axes from input files; write axes to output file.
+
#Allocate arrays and read wavelet, velocity and reflectivity. <syntaxhighlight lang="fortran">
<syntaxhighlight lang="fortran">
+
   allocate(ww(nt));    call rsf_read(Fw,ww)
   call iaxa(Fw,at,1); call iaxa(Fv,az,1); call iaxa(Fv,ax,2)
+
  allocate(vv(nz,nx)); call rsf_read(Fv,vv)
   call oaxa(Fo,az,1); call oaxa(Fo,ax,2); call oaxa(Fo,at,3)
+
   allocate(rr(nz,nx)); call rsf_read(Fr,rr)
 
</syntaxhighlight>   
 
</syntaxhighlight>   
*Allocate arrays and read wavelet, velocity and reflectivity.
+
#Allocate temporary arrays. <syntaxhighlight lang="fortran">
<syntaxhighlight lang="fortran">
+
   allocate(um(nz,nx)); um=0.
   allocate(ww(at%n    )); ww=0.; call rsf_read(Fw,ww)
+
  allocate(uo(nz,nx)); uo=0.
   allocate(vv(az%n,ax%n)); vv=0.; call rsf_read(Fv,vv)
+
   allocate(up(nz,nx)); up=0.
   allocate(rr(az%n,ax%n)); rr=0.; call rsf_read(Fr,rr)
+
   allocate(ud(nz,nx)); ud=0.
 
</syntaxhighlight>   
 
</syntaxhighlight>   
*Allocate temporary arrays. 
+
#Loop over time. <syntaxhighlight lang="fortran">
<syntaxhighlight lang="fortran">
+
   do it=1,nt
  allocate(um(az%n,ax%n)); um=0.
+
</syntaxhighlight>  
  allocate(uo(az%n,ax%n)); uo=0.
+
#Compute Laplacian: <math>\Delta U</math>. <syntaxhighlight lang="fortran">
  allocate(up(az%n,ax%n)); up=0.
+
     ud(3:nz-2,3:nx-2) = &
  allocate(ud(az%n,ax%n)); ud=0.
+
          c0* uo(3:nz-2,3:nx-2) * (idx + idz)           + &
</syntaxhighlight> 
+
          c1*(uo(3:nz-2,2:nx-3) + uo(3:nz-2,4:nx-1))*idx + &
*Loop over time.
+
          c2*(uo(3:nz-2,1:nx-4) + uo(3:nz-2,5:nx ))*idx + &
<syntaxhighlight lang="fortran">
+
          c1*(uo(2:nz-3,3:nx-2) + uo(4:nz-1,3:nx-2))*idz + &
   do it=1,at%n
+
          c2*(uo(1:nz-4,3:nx-2) + uo(5:nz  ,3:nx-2))*idz
</syntaxhighlight>  
+
</syntaxhighlight>
*Compute Laplacian: <math>\Delta U</math>.
+
#Inject source wavelet: <math>\left[ \Delta U - f(t) \right]</math> <syntaxhighlight lang="fortran">
<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: <math>\left[ \Delta U - f(t) \right]</math>
 
<syntaxhighlight lang="fortran">
 
 
     ud = ud - ww(it) * rr
 
     ud = ud - ww(it) * rr
</syntaxhighlight>
+
</syntaxhighlight>
*Scale by velocity: <math>\left[ \Delta U - f(t) \right] v^2</math>
+
#Scale by velocity: <math>\left[ \Delta U - f(t) \right] v^2</math> <syntaxhighlight lang="fortran">
<syntaxhighlight lang="fortran">
 
 
     ud= ud *vv*vv
 
     ud= ud *vv*vv
</syntaxhighlight>
+
</syntaxhighlight>
*Time step:  <math>U_{i+1} = \left[ \Delta U -f(t) \right] v^2 \Delta t^2 + 2 U_{i} - U_{i-1}</math>  
+
#Time step:  <math>U_{i+1} = \left[ \Delta U -f(t) \right] v^2 \Delta t^2 + 2 U_{i} - U_{i-1}</math> <syntaxhighlight lang="fortran">
<syntaxhighlight lang="fortran">
 
 
     up = 2*uo - um + ud * dt2
 
     up = 2*uo - um + ud * dt2
 
     um =  uo
 
     um =  uo
 
     uo =  up
 
     uo =  up
 +
</syntaxhighlight> 
 +
 +
==Python program==
 +
 +
<syntaxhighlight lang="python">
 +
#!/usr/bin/env python
 +
 +
import sys
 +
import numpy
 +
import m8r
 +
 +
c0=-30./12.
 +
c1=+16./12.
 +
c2=- 1./12.
 +
 +
par = m8r.Par()
 +
verb = par.bool("verb",False) # verbosity
 +
 +
# setup I/O files
 +
Fw=m8r.Input()
 +
Fv=m8r.Input ("vel")
 +
Fr=m8r.Input ("ref")
 +
Fo=m8r.Output()
 +
 +
# Read/Write axes
 +
at = Fw.axis(1); nt = at['n']; dt = at['d']
 +
az = Fv.axis(1); nz = az['n']; dz = az['d']
 +
ax = Fv.axis(2); nx = ax['n']; dx = ax['d']
 +
 +
Fo.putaxis(az,1)
 +
Fo.putaxis(ax,2)
 +
Fo.putaxis(at,3)
 +
 +
dt2 =    dt*dt
 +
idz = 1/(dz*dz)
 +
idx = 1/(dx*dx)
 +
 +
# read wavelet, velocity & reflectivity
 +
ww = numpy.zeros(nt,'f');      Fw.read(ww)
 +
vv = numpy.zeros([nz,nx],'f'); Fv.read(vv)
 +
rr = numpy.zeros([nz,nx],'f'); Fr.read(rr)
 +
 +
# allocate temporary arrays
 +
um = numpy.zeros([nz,nx],'f')
 +
uo = numpy.zeros([nz,nx],'f')
 +
up = numpy.zeros([nz,nx],'f')
 +
ud = numpy.zeros([nz,nx],'f')
 +
 +
# MAIN LOOP
 +
for it in range(nt):
 +
    if verb:
 +
        sys.stderr.write("\b\b\b\b\b %d" % it)
 +
 +
    ud[2:-2,2:-2] = \
 +
    c0* uo[2:-2,2:-2] * (idx + idz)        + \
 +
    c1*(uo[2:-2,1:-3] + uo[2:-2,3:-1])*idx + \
 +
    c2*(uo[2:-2, :-4] + uo[2:-2,4:  ])*idx + \
 +
    c1*(uo[1:-3,2:-2] + uo[3:-1,2:-2])*idz + \
 +
    c2*(uo[ :-4,2:-2] + uo[4:  ,2:-2])*idz
 +
 +
    # 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
 +
 +
if verb:
 +
    sys.stderr.write("\n")
 +
    Fo.write(uo)
 +
 +
sys.exit(0)
 
</syntaxhighlight>
 
</syntaxhighlight>
 +
 +
[[Image:wavepython.png|frame|center|Wave propagation snapshot.]]
 +
 +
#Open files for input/output. <syntaxhighlight lang="python">
 +
Fw=m8r.Input()
 +
Fv=m8r.Input ("vel")
 +
Fr=m8r.Input ("ref")
 +
Fo=m8r.Output()
 +
</syntaxhighlight> 
 +
#Read axes from input files; write axes to output file. <syntaxhighlight lang="python">
 +
at = Fw.axis(1); nt = at['n']; dt = at['d']
 +
az = Fv.axis(1); nz = az['n']; dz = az['d']
 +
ax = Fv.axis(2); nx = ax['n']; dx = ax['d']
 +
 +
Fo.putaxis(az,1)
 +
Fo.putaxis(ax,2)
 +
Fo.putaxis(at,3)
 +
</syntaxhighlight> 
 +
#Allocate arrays and read wavelet, velocity and reflectivity. <syntaxhighlight lang="python">
 +
ww = numpy.zeros(nt,'f');      Fw.read(ww)
 +
vv = numpy.zeros([nz,nx],'f'); Fv.read(vv)
 +
rr = numpy.zeros([nz,nx],'f'); Fr.read(rr)
 +
</syntaxhighlight> 
 +
#Allocate temporary arrays. <syntaxhighlight lang="python">
 +
um = numpy.zeros([nz,nx],'f')
 +
uo = numpy.zeros([nz,nx],'f')
 +
up = numpy.zeros([nz,nx],'f')
 +
ud = numpy.zeros([nz,nx],'f')
 +
</syntaxhighlight> 
 +
#Loop over time. <syntaxhighlight lang="python">
 +
for it in range(nt):
 +
</syntaxhighlight> 
 +
#Compute Laplacian: <math>\Delta U</math>. <syntaxhighlight lang="python">
 +
    ud[2:-2,2:-2] = \
 +
    c0* uo[2:-2,2:-2] * (idx + idz)        + \
 +
    c1*(uo[2:-2,1:-3] + uo[2:-2,3:-1])*idx + \
 +
    c2*(uo[2:-2, :-4] + uo[2:-2,4:  ])*idx + \
 +
    c1*(uo[1:-3,2:-2] + uo[3:-1,2:-2])*idz + \
 +
    c2*(uo[ :-4,2:-2] + uo[4:  ,2:-2])*idz
 +
</syntaxhighlight> 
 +
#Inject source wavelet: <math>\left[ \Delta U - f(t) \right]</math> <syntaxhighlight lang="python">
 +
    ud = ud - ww[it] * rr
 +
</syntaxhighlight> 
 +
#Scale by velocity: <math>\left[ \Delta U - f(t) \right] v^2</math> <syntaxhighlight lang="python">
 +
    ud= ud *vv*vv
 +
</syntaxhighlight> 
 +
#Time step:  <math>U_{i+1} = \left[ \Delta U -f(t) \right] v^2 \Delta t^2 + 2 U_{i} - U_{i-1}</math>  <syntaxhighlight lang="python">
 +
    up = 2*uo - um + ud * dt2
 +
    um =  uo
 +
    uo =  up
 +
</syntaxhighlight> 
  
 
==References==
 
==References==
 
<references/>
 
<references/>

Latest revision as of 18:29, 5 May 2021

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.

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

can be written as

is the Laplacian symbol, is the source wavelet, is the velocity, and is a scalar wavefield. A discrete time-step involves the following computations:

where , and represent the propagating wavefield at various time steps.

C program

Wave propagation snapshot.

C program

/* time-domain acoustic FD modeling */
#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; /* verbose flag */

    /* 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 (ix=0; ix<nx; ix++) {
	for (iz=0; iz<nz; iz++) {
	    um[ix][iz]=0;
	    uo[ix][iz]=0;
	    up[ix][iz]=0;
	    ud[ix][iz]=0;
	}
    }

    /* MAIN LOOP */
    if(verb) fprintf(stderr,"
");
    for (it=0; it<nt; it++) {
	if(verb) fprintf(stderr,"\b\b\b\b\b %d",it);
	
	/* 4th order laplacian */
	for (ix=2; ix<nx-2; ix++) {
	    for (iz=2; iz<nz-2; iz++) {
		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 (ix=0; ix<nx; ix++) {
	    for (iz=0; iz<nz; iz++) {
		ud[ix][iz] -= ww[it] * rr[ix][iz];
	    }
	}

	/* scale by velocity */
	for (ix=0; ix<nx; ix++) {
	    for (iz=0; iz<nz; iz++) {
		ud[ix][iz] *= vv[ix][iz]*vv[ix][iz];
	    }
	}
	
	/* time step */
	for (ix=0; ix<nx; ix++) {
	    for (iz=0; iz<nz; iz++) {
		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");

    exit (0);
}


  1. Declare input, output and auxiliary file tags: Fw for input wavelet, Fv for velocity, Fr for reflectivity, and Fo for output wavefield.
        sf_file Fw=NULL,Fv=NULL,Fr=NULL,Fo=NULL; /* I/O files */
    
  2. Declare RSF cube axes: at time axis, ax space axis, az depth axis. \tiny
        sf_axis at,az,ax;    /* cube axes */
    
  3. Declare multi-dimensional arrays for input, output and computations.
        float  *ww,**vv,**rr;     /* I/O arrays*/
    
  4. Open files for input/output.
        Fw = sf_input ("in" );
        Fo = sf_output("out");
        Fv = sf_input ("vel");
        Fr = sf_input ("ref");
    
  5. Read axes from input files; write axes to output file.
        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);
    
  6. Allocate arrays and read wavelet, velocity and 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);
    
  7. Allocate temporary arrays.
        um=sf_floatalloc2(nz,nx);
        uo=sf_floatalloc2(nz,nx);
        up=sf_floatalloc2(nz,nx);
        ud=sf_floatalloc2(nz,nx);
    
  8. Loop over time.
        for (it=0; it<nt; it++) {
    
  9. Compute Laplacian: .
    	for (ix=2; ix<nx-2; ix++) {
    	    for (iz=2; iz<nz-2; iz++) {
    		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;	  
    	    }
    	}
    
  10. Inject source wavelet:
    	for (ix=0; ix<nx; ix++) {
    	    for (iz=0; iz<nz; iz++) {
    		ud[ix][iz] -= ww[it] * rr[ix][iz];
    	    }
    	}
    
  11. Scale by velocity:
    	for (ix=0; ix<nx; ix++) {
    	    for (iz=0; iz<nz; iz++) {
    		ud[ix][iz] *= vv[ix][iz]*vv[ix][iz];
    	    }
    	}
    
  12. Time step:
    	for (ix=0; ix<nx; ix++) {
    	    for (iz=0; iz<nz; iz++) {
    		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];
    	    }
    	}
    


\newpage

C++ program

// time-domain acoustic FD modeling
#include <valarray>
#include <iostream>
#include <rsf.hh>
#include <cub.hh>

#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); 

    // 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 ix=2; ix<nx-2; ix++) {
	    for (int iz=2; iz<nz-2; iz++) {
		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);
}
  1. Declare input, output and auxiliary file cubes (of type CUB).
        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);
    
  2. Declare, read and write RSF cube axes: at time axis, ax space axis, az depth axis.
        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);
    
  3. Declare multi-dimensional valarrays for input, output and read data.
        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;
    
  4. Declare multi-dimensional valarrays for temporary storage.
        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;
    
  5. Initialize multidimensional valarray index counter (of type VAI).
        VAI k(nz,nx);
    
  6. Loop over time.
        for (int it=0; it<nt; it++) {
    
  7. Compute Laplacian: .
    	for (int ix=2; ix<nx-2; ix++) {
    	    for (int iz=2; iz<nz-2; iz++) {
    		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;
    	    }
    	}
    
  8. Inject source wavelet:
    	ud -= ww[it] * rr;
    
  9. Scale by velocity:
    	ud *= vv*vv;
    
  10. Time step:
    	up=(float)2 * uo - um + ud * dt2;
    	um =   uo;
    	uo =   up;
    

Fortran 90 program

! 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
  integer    :: nt,nz,nx
  real       :: dt,dz,dx
  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); nt = at%n; dt = at%d
  call iaxa(Fv,az,1); nz = az%n; dz = az%d
  call iaxa(Fv,ax,2); nx = ax%n; dx = ax%d
  
  call oaxa(Fo,az,1)
  call oaxa(Fo,ax,2)
  call oaxa(Fo,at,3)

  dt2 =    dt*dt
  idz = 1/(dz*dz)
  idx = 1/(dx*dx) 

  ! read wavelet, velocity & reflectivity
  allocate(ww(nt));    call rsf_read(Fw,ww)
  allocate(vv(nz,nx)); call rsf_read(Fv,vv)
  allocate(rr(nz,nx)); call rsf_read(Fr,rr)

  ! allocate temporary arrays
  allocate(um(nz,nx)); um=0.
  allocate(uo(nz,nx)); uo=0.
  allocate(up(nz,nx)); up=0.
  allocate(ud(nz,nx)); ud=0.

  ! MAIN LOOP
  do it=1,nt
     if(verb) write (0,*) it

     ud(3:nz-2,3:nx-2) = &
          c0* uo(3:nz-2,3:nx-2) * (idx + idz)            + &
          c1*(uo(3:nz-2,2:nx-3) + uo(3:nz-2,4:nx-1))*idx + &
          c2*(uo(3:nz-2,1:nx-4) + uo(3:nz-2,5:nx  ))*idx + &
          c1*(uo(2:nz-3,3:nx-2) + uo(4:nz-1,3:nx-2))*idz + &
          c2*(uo(1:nz-4,3:nx-2) + uo(5:nz  ,3:nx-2))*idz

     ! 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
  1. Declare input, output and auxiliary file tags.
      type(file) :: Fw,Fv,Fr,Fo  ! I/O files
    
  2. Declare RSF cube axes: at time axis, ax space axis, az depth axis.
      type(axa)  :: at,az,ax     ! cube axes
    
  3. Declare multi-dimensional arrays for input, output and computations.
      real, allocatable :: vv(:,:),rr(:,:),ww(:)           ! I/O arrays
      real, allocatable :: um(:,:),uo(:,:),up(:,:),ud(:,:) ! tmp arrays
    
  4. Open files for input/output.
      Fw=rsf_input ("in")
      Fv=rsf_input ("vel")
      Fr=rsf_input ("ref")
      Fo=rsf_output("out")
    
  5. Read axes from input files; write axes to output file.
      call iaxa(Fw,at,1); nt = at%n; dt = at%d
      call iaxa(Fv,az,1); nz = az%n; dz = az%d
      call iaxa(Fv,ax,2); nx = ax%n; dx = ax%d
      
      call oaxa(Fo,az,1)
      call oaxa(Fo,ax,2)
      call oaxa(Fo,at,3)
    
  6. Allocate arrays and read wavelet, velocity and reflectivity.
      allocate(ww(nt));    call rsf_read(Fw,ww)
      allocate(vv(nz,nx)); call rsf_read(Fv,vv)
      allocate(rr(nz,nx)); call rsf_read(Fr,rr)
    
  7. Allocate temporary arrays.
      allocate(um(nz,nx)); um=0.
      allocate(uo(nz,nx)); uo=0.
      allocate(up(nz,nx)); up=0.
      allocate(ud(nz,nx)); ud=0.
    
  8. Loop over time.
      do it=1,nt
    
  9. Compute Laplacian: .
         ud(3:nz-2,3:nx-2) = &
              c0* uo(3:nz-2,3:nx-2) * (idx + idz)            + &
              c1*(uo(3:nz-2,2:nx-3) + uo(3:nz-2,4:nx-1))*idx + &
              c2*(uo(3:nz-2,1:nx-4) + uo(3:nz-2,5:nx  ))*idx + &
              c1*(uo(2:nz-3,3:nx-2) + uo(4:nz-1,3:nx-2))*idz + &
              c2*(uo(1:nz-4,3:nx-2) + uo(5:nz  ,3:nx-2))*idz
    
  10. Inject source wavelet:
         ud = ud - ww(it) * rr
    
  11. Scale by velocity:
         ud= ud *vv*vv
    
  12. Time step:
         up = 2*uo - um + ud * dt2
         um =   uo
         uo =   up
    

Python program

#!/usr/bin/env python

import sys
import numpy
import m8r

c0=-30./12.
c1=+16./12.
c2=- 1./12.

par = m8r.Par()
verb = par.bool("verb",False) # verbosity

# setup I/O files
Fw=m8r.Input()
Fv=m8r.Input ("vel")
Fr=m8r.Input ("ref")
Fo=m8r.Output()

# Read/Write axes
at = Fw.axis(1); nt = at['n']; dt = at['d']
az = Fv.axis(1); nz = az['n']; dz = az['d']
ax = Fv.axis(2); nx = ax['n']; dx = ax['d']

Fo.putaxis(az,1)
Fo.putaxis(ax,2)
Fo.putaxis(at,3)

dt2 =    dt*dt
idz = 1/(dz*dz)
idx = 1/(dx*dx) 

# read wavelet, velocity & reflectivity
ww = numpy.zeros(nt,'f');      Fw.read(ww)
vv = numpy.zeros([nz,nx],'f'); Fv.read(vv)
rr = numpy.zeros([nz,nx],'f'); Fr.read(rr)

# allocate temporary arrays
um = numpy.zeros([nz,nx],'f')
uo = numpy.zeros([nz,nx],'f')
up = numpy.zeros([nz,nx],'f')
ud = numpy.zeros([nz,nx],'f')

# MAIN LOOP
for it in range(nt):
    if verb:
        sys.stderr.write("\b\b\b\b\b %d" % it)

    ud[2:-2,2:-2] = \
    c0* uo[2:-2,2:-2] * (idx + idz)        + \
    c1*(uo[2:-2,1:-3] + uo[2:-2,3:-1])*idx + \
    c2*(uo[2:-2, :-4] + uo[2:-2,4:  ])*idx + \
    c1*(uo[1:-3,2:-2] + uo[3:-1,2:-2])*idz + \
    c2*(uo[ :-4,2:-2] + uo[4:  ,2:-2])*idz

    # 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

if verb:
    sys.stderr.write("\n")
    Fo.write(uo)

sys.exit(0)
Wave propagation snapshot.
  1. Open files for input/output.
    Fw=m8r.Input()
    Fv=m8r.Input ("vel")
    Fr=m8r.Input ("ref")
    Fo=m8r.Output()
    
  2. Read axes from input files; write axes to output file.
    at = Fw.axis(1); nt = at['n']; dt = at['d']
    az = Fv.axis(1); nz = az['n']; dz = az['d']
    ax = Fv.axis(2); nx = ax['n']; dx = ax['d']
    
    Fo.putaxis(az,1)
    Fo.putaxis(ax,2)
    Fo.putaxis(at,3)
    
  3. Allocate arrays and read wavelet, velocity and reflectivity.
    ww = numpy.zeros(nt,'f');      Fw.read(ww)
    vv = numpy.zeros([nz,nx],'f'); Fv.read(vv)
    rr = numpy.zeros([nz,nx],'f'); Fr.read(rr)
    
  4. Allocate temporary arrays.
    um = numpy.zeros([nz,nx],'f')
    uo = numpy.zeros([nz,nx],'f')
    up = numpy.zeros([nz,nx],'f')
    ud = numpy.zeros([nz,nx],'f')
    
  5. Loop over time.
    for it in range(nt):
    
  6. Compute Laplacian: .
        ud[2:-2,2:-2] = \
        c0* uo[2:-2,2:-2] * (idx + idz)        + \
        c1*(uo[2:-2,1:-3] + uo[2:-2,3:-1])*idx + \
        c2*(uo[2:-2, :-4] + uo[2:-2,4:  ])*idx + \
        c1*(uo[1:-3,2:-2] + uo[3:-1,2:-2])*idz + \
        c2*(uo[ :-4,2:-2] + uo[4:  ,2:-2])*idz
    
  7. Inject source wavelet:
        ud = ud - ww[it] * rr
    
  8. Scale by velocity:
        ud= ud *vv*vv
    
  9. Time step:
        up = 2*uo - um + ud * dt2
        um =   uo
        uo =   up
    

References

  1. "Hello world" of seismic imaging.