# Difference between revisions of "Guide to programming with madagascar"

Jump to navigation Jump to 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.

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

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

can be written as

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

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

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

where ${\displaystyle U_{i-1}}$, ${\displaystyle U_{i}}$ and ${\displaystyle U_{i+1}}$ 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: ${\displaystyle \Delta U}$.
	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: ${\displaystyle \left[\Delta U-f(t)\right]}$
	for (ix=0; ix<nx; ix++) {
for (iz=0; iz<nz; iz++) {
ud[ix][iz] -= ww[it] * rr[ix][iz];
}
}

11. Scale by velocity: ${\displaystyle \left[\Delta U-f(t)\right]v^{2}}$
	for (ix=0; ix<nx; ix++) {
for (iz=0; iz<nz; iz++) {
ud[ix][iz] *= vv[ix][iz]*vv[ix][iz];
}
}

12. Time step: ${\displaystyle U_{i+1}=\left[\Delta U-f(t)\right]v^{2}\Delta t^{2}+2U_{i}-U_{i-1}}$
	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: ${\displaystyle \Delta U}$.
	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: ${\displaystyle \left[\Delta U-f(t)\right]}$
	ud -= ww[it] * rr;

9. Scale by velocity: ${\displaystyle \left[\Delta U-f(t)\right]v^{2}}$
	ud *= vv*vv;

10. Time step: ${\displaystyle U_{i+1}=\left[\Delta U-f(t)\right]v^{2}\Delta t^{2}+2U_{i}-U_{i-1}}$
	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: ${\displaystyle \Delta U}$.
     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: ${\displaystyle \left[\Delta U-f(t)\right]}$
     ud = ud - ww(it) * rr

11. Scale by velocity: ${\displaystyle \left[\Delta U-f(t)\right]v^{2}}$
     ud= ud *vv*vv

12. Time step: ${\displaystyle U_{i+1}=\left[\Delta U-f(t)\right]v^{2}\Delta t^{2}+2U_{i}-U_{i-1}}$
     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: ${\displaystyle \Delta U}$.
    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: ${\displaystyle \left[\Delta U-f(t)\right]}$
    ud = ud - ww[it] * rr

8. Scale by velocity: ${\displaystyle \left[\Delta U-f(t)\right]v^{2}}$
    ud= ud *vv*vv

9. Time step: ${\displaystyle U_{i+1}=\left[\Delta U-f(t)\right]v^{2}\Delta t^{2}+2U_{i}-U_{i-1}}$
    up = 2*uo - um + ud * dt2
um =   uo
uo =   up


## References

1. "Hello world" of seismic imaging.