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
m (Demo moved to Guide to programming with madagascar) |
(update + python) |
||
(11 intermediate revisions by 2 users not shown) | |||
Line 1: | Line 1: | ||
− | <center><font size="-1">''This page was created from the LaTeX source in [ | + | <center><font size="-1">''This page was created from the LaTeX source in [https://github.com/ahay/src/blob/master/book/rsf/rsf/demo.tex book/rsf/rsf/demo.tex] using [[latex2wiki]]''</font></center> |
This guide demonstrates a simple time-domain | This guide demonstrates a simple time-domain | ||
Line 32: | Line 32: | ||
==C program== | ==C program== | ||
+ | [[Image:wavec.png|frame|center|Wave propagation snapshot.]] | ||
− | <c> | + | ==C program== |
+ | |||
+ | <syntaxhighlight lang="c"> | ||
/* time-domain acoustic FD modeling */ | /* time-domain acoustic FD modeling */ | ||
#include <rsf.h> | #include <rsf.h> | ||
Line 42: | Line 45: | ||
bool verb; /* verbose flag */ | bool verb; /* verbose flag */ | ||
− | sf_file Fw,Fv,Fr,Fo; /* I/O files */ | + | sf_file Fw=NULL,Fv=NULL,Fr=NULL,Fo=NULL; /* I/O files */ |
sf_axis at,az,ax; /* cube axes */ | sf_axis at,az,ax; /* cube axes */ | ||
int it,iz,ix; /* index variables */ | int it,iz,ix; /* index variables */ | ||
Line 52: | Line 55: | ||
sf_init(argc,argv); | sf_init(argc,argv); | ||
− | if(! sf_getbool("verb",&verb)) verb=0; | + | if(! sf_getbool("verb",&verb)) verb=0; /* verbose flag */ |
/* setup I/O files */ | /* setup I/O files */ | ||
Line 83: | Line 86: | ||
up=sf_floatalloc2(nz,nx); | up=sf_floatalloc2(nz,nx); | ||
ud=sf_floatalloc2(nz,nx); | ud=sf_floatalloc2(nz,nx); | ||
− | + | ||
− | for ( | + | for (ix=0; ix<nx; ix++) { |
− | for ( | + | for (iz=0; iz<nz; iz++) { |
um[ix][iz]=0; | um[ix][iz]=0; | ||
uo[ix][iz]=0; | uo[ix][iz]=0; | ||
Line 92: | Line 95: | ||
} | } | ||
} | } | ||
− | + | ||
/* MAIN LOOP */ | /* MAIN LOOP */ | ||
− | if(verb) fprintf(stderr," | + | 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 */ | ||
− | for ( | + | for (ix=2; ix<nx-2; ix++) { |
− | for ( | + | for (iz=2; iz<nz-2; iz++) { |
ud[ix][iz] = | ud[ix][iz] = | ||
c0* uo[ix ][iz ] * (idx+idz) + | c0* uo[ix ][iz ] * (idx+idz) + | ||
Line 111: | Line 115: | ||
/* inject wavelet */ | /* inject wavelet */ | ||
− | for ( | + | for (ix=0; ix<nx; ix++) { |
− | for ( | + | for (iz=0; iz<nz; iz++) { |
ud[ix][iz] -= ww[it] * rr[ix][iz]; | ud[ix][iz] -= ww[it] * rr[ix][iz]; | ||
} | } | ||
Line 118: | Line 122: | ||
/* scale by velocity */ | /* scale by velocity */ | ||
− | for ( | + | for (ix=0; ix<nx; ix++) { |
− | for ( | + | for (iz=0; iz<nz; iz++) { |
ud[ix][iz] *= vv[ix][iz]*vv[ix][iz]; | ud[ix][iz] *= vv[ix][iz]*vv[ix][iz]; | ||
} | } | ||
Line 125: | Line 129: | ||
/* time step */ | /* time step */ | ||
− | for ( | + | for (ix=0; ix<nx; ix++) { |
− | for ( | + | for (iz=0; iz<nz; iz++) { |
up[ix][iz] = | up[ix][iz] = | ||
2*uo[ix][iz] | 2*uo[ix][iz] | ||
Line 140: | Line 144: | ||
sf_floatwrite(uo[0],nz*nx,Fo); | sf_floatwrite(uo[0],nz*nx,Fo); | ||
} | } | ||
− | if(verb) fprintf(stderr,"\n"); | + | if(verb) fprintf(stderr,"\n"); |
exit (0); | exit (0); | ||
} | } | ||
− | </ | + | </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"> | |
− | + | sf_file Fw=NULL,Fv=NULL,Fr=NULL,Fo=NULL; /* I/O files */ | |
− | <c> | + | </syntaxhighlight> |
− | sf_file Fw,Fv,Fr,Fo; /* I/O files */ | + | #Declare RSF cube axes: <tt>at</tt> time axis, <tt>ax</tt> space axis, <tt>az</tt> depth axis. \tiny <syntaxhighlight lang="c"> |
− | </ | ||
− | |||
− | <c> | ||
sf_axis at,az,ax; /* cube axes */ | sf_axis at,az,ax; /* cube axes */ | ||
− | </ | + | </syntaxhighlight> |
− | + | #Declare multi-dimensional arrays for input, output and computations. <syntaxhighlight lang="c"> | |
− | <c> | + | |
float *ww,**vv,**rr; /* I/O arrays*/ | float *ww,**vv,**rr; /* I/O arrays*/ | ||
− | </ | + | </syntaxhighlight> |
− | + | #Open files for input/output. <syntaxhighlight lang="c"> | |
− | <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> |
− | + | #Read axes from input files; write axes to output file. <syntaxhighlight lang="c"> | |
− | <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 176: | Line 175: | ||
sf_oaxa(Fo,ax,2); | sf_oaxa(Fo,ax,2); | ||
sf_oaxa(Fo,at,3); | sf_oaxa(Fo,at,3); | ||
− | </ | + | </syntaxhighlight> |
− | + | #Allocate arrays and read wavelet, velocity and reflectivity. <syntaxhighlight lang="c"> | |
− | <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> |
− | + | #Allocate temporary arrays. <syntaxhighlight lang="c"> | |
− | <c> | ||
um=sf_floatalloc2(nz,nx); | um=sf_floatalloc2(nz,nx); | ||
uo=sf_floatalloc2(nz,nx); | uo=sf_floatalloc2(nz,nx); | ||
up=sf_floatalloc2(nz,nx); | up=sf_floatalloc2(nz,nx); | ||
ud=sf_floatalloc2(nz,nx); | ud=sf_floatalloc2(nz,nx); | ||
− | </ | + | </syntaxhighlight> |
− | + | #Loop over time. <syntaxhighlight lang="c"> | |
− | <c> | ||
for (it=0; it<nt; it++) { | for (it=0; it<nt; it++) { | ||
− | </ | + | </syntaxhighlight> |
− | + | #Compute Laplacian: <math>\Delta U</math>. <syntaxhighlight lang="c"> | |
− | <c> | + | for (ix=2; ix<nx-2; ix++) { |
− | for ( | + | for (iz=2; iz<nz-2; iz++) { |
− | for ( | ||
ud[ix][iz] = | ud[ix][iz] = | ||
c0* uo[ix ][iz ] * (idx+idz) + | c0* uo[ix ][iz ] * (idx+idz) + | ||
Line 206: | Line 201: | ||
} | } | ||
} | } | ||
− | </ | + | </syntaxhighlight> |
− | + | #Inject source wavelet: <math>\left[ \Delta U - f(t) \right]</math> <syntaxhighlight lang="c"> | |
− | <c> | + | for (ix=0; ix<nx; ix++) { |
− | for ( | + | for (iz=0; iz<nz; iz++) { |
− | for ( | ||
ud[ix][iz] -= ww[it] * rr[ix][iz]; | ud[ix][iz] -= ww[it] * rr[ix][iz]; | ||
} | } | ||
} | } | ||
− | </ | + | </syntaxhighlight> |
− | + | #Scale by velocity: <math>\left[ \Delta U - f(t) \right] v^2</math> <syntaxhighlight lang="c"> | |
− | <c> | + | for (ix=0; ix<nx; ix++) { |
− | for ( | + | for (iz=0; iz<nz; iz++) { |
− | for ( | ||
ud[ix][iz] *= vv[ix][iz]*vv[ix][iz]; | ud[ix][iz] *= vv[ix][iz]*vv[ix][iz]; | ||
} | } | ||
} | } | ||
− | </ | + | </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="c"> | |
− | <c> | + | for (ix=0; ix<nx; ix++) { |
− | for ( | + | for (iz=0; iz<nz; iz++) { |
− | for ( | ||
up[ix][iz] = | up[ix][iz] = | ||
2*uo[ix][iz] | 2*uo[ix][iz] | ||
Line 236: | Line 228: | ||
} | } | ||
} | } | ||
− | </ | + | </syntaxhighlight> |
+ | |||
+ | \newpage | ||
==C++ program== | ==C++ program== | ||
− | + | <syntaxhighlight lang="cpp"> | |
− | <cpp> | ||
// time-domain acoustic FD modeling | // time-domain acoustic FD modeling | ||
#include <valarray> | #include <valarray> | ||
Line 247: | Line 240: | ||
#include <rsf.hh> | #include <rsf.hh> | ||
#include <cub.hh> | #include <cub.hh> | ||
− | #include | + | |
+ | #include "vai.hh" | ||
+ | |||
using namespace std; | using namespace std; | ||
Line 263: | 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 | + | CUB Fo("out","o"); Fo.setup(3); |
// Read/Write axes | // Read/Write axes | ||
Line 296: | Line 291: | ||
if(verb) cerr << endl; | if(verb) cerr << endl; | ||
for (int it=0; it<nt; it++) { | for (int it=0; it<nt; it++) { | ||
− | if(verb) cerr << " | + | if(verb) cerr << "\b\b\b\b\b" << it; |
// 4th order laplacian | // 4th order laplacian | ||
− | for (int | + | for (int ix=2; ix<nx-2; ix++) { |
− | for (int | + | 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 329: | Line 324: | ||
} | } | ||
− | </ | + | </syntaxhighlight> |
− | + | #Declare input, output and auxiliary file cubes (of type <tt>CUB</tt>). <syntaxhighlight lang="cpp"> | |
− | <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 | + | CUB Fo("out","o"); Fo.setup(3); |
− | </ | + | </syntaxhighlight> |
− | + | #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"> | |
− | <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> | |
− | + | #Declare multi-dimensional <tt>valarrays</tt> for input, output and read data. <syntaxhighlight lang="cpp"> | |
− | |||
− | |||
− | |||
− | </ | ||
− | |||
− | <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> |
− | + | #Declare multi-dimensional <tt>valarrays</tt> for temporary storage. <syntaxhighlight lang="cpp"> | |
− | <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; | ||
valarray<float> up(nz*nx); up=0; | valarray<float> up(nz*nx); up=0; | ||
valarray<float> ud(nz*nx); ud=0; | valarray<float> ud(nz*nx); ud=0; | ||
− | </ | + | </syntaxhighlight> |
− | + | #Initialize multidimensional <tt>valarray</tt> index counter (of type <tt>VAI</tt>). <syntaxhighlight lang="cpp"> | |
− | <cpp> | ||
VAI k(nz,nx); | VAI k(nz,nx); | ||
− | </ | + | </syntaxhighlight> |
− | + | #Loop over time. <syntaxhighlight lang="cpp"> | |
− | <cpp> | ||
for (int it=0; it<nt; it++) { | for (int it=0; it<nt; it++) { | ||
− | </ | + | </syntaxhighlight> |
− | + | #Compute Laplacian: <math>\Delta U</math>. <syntaxhighlight lang="cpp"> | |
− | <cpp> | + | for (int ix=2; ix<nx-2; ix++) { |
− | for (int | + | for (int iz=2; iz<nz-2; iz++) { |
− | for (int | ||
ud[k(iz,ix)] = | ud[k(iz,ix)] = | ||
c0* uo[ k(iz ,ix )] * (idx+idz) + | c0* uo[ k(iz ,ix )] * (idx+idz) + | ||
Line 382: | Line 365: | ||
} | } | ||
} | } | ||
− | </ | + | </syntaxhighlight> |
− | + | #Inject source wavelet: <math>\left[ \Delta U - f(t) \right]</math> <syntaxhighlight lang="cpp"> | |
− | <cpp> | ||
ud -= ww[it] * rr; | ud -= ww[it] * rr; | ||
− | </ | + | </syntaxhighlight> |
− | + | #Scale by velocity: <math>\left[ \Delta U - f(t) \right] v^2</math> <syntaxhighlight lang="cpp"> | |
− | <cpp> | ||
ud *= vv*vv; | 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="cpp"> | |
− | <cpp> | ||
up=(float)2 * uo - um + ud * dt2; | up=(float)2 * uo - um + ud * dt2; | ||
um = uo; | um = uo; | ||
uo = up; | uo = up; | ||
− | </ | + | </syntaxhighlight> |
==Fortran 90 program== | ==Fortran 90 program== | ||
− | + | <syntaxhighlight lang="fortran"> | |
− | <fortran> | ||
! time-domain acoustic FD modeling | ! time-domain acoustic FD modeling | ||
program AFDMf90 | program AFDMf90 | ||
Line 415: | 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 430: | 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 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 = | + | dt2 = dt*dt |
− | idz = 1/( | + | idz = 1/(dz*dz) |
− | idx = 1/( | + | idx = 1/(dx*dx) |
! read wavelet, velocity & reflectivity | ! read wavelet, velocity & reflectivity | ||
− | allocate(ww( | + | allocate(ww(nt)); call rsf_read(Fw,ww) |
− | allocate(vv( | + | allocate(vv(nz,nx)); call rsf_read(Fv,vv) |
− | allocate(rr( | + | allocate(rr(nz,nx)); call rsf_read(Fr,rr) |
! allocate temporary arrays | ! allocate temporary arrays | ||
− | allocate(um( | + | allocate(um(nz,nx)); um=0. |
− | allocate(uo( | + | allocate(uo(nz,nx)); uo=0. |
− | allocate(up( | + | allocate(up(nz,nx)); up=0. |
− | allocate(ud( | + | allocate(ud(nz,nx)); ud=0. |
! MAIN LOOP | ! MAIN LOOP | ||
− | do it=1, | + | do it=1,nt |
if(verb) write (0,*) it | 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 | ! inject wavelet | ||
Line 481: | Line 462: | ||
call exit(0) | call exit(0) | ||
end program AFDMf90 | end program AFDMf90 | ||
− | </ | + | </syntaxhighlight> |
− | + | ||
− | + | #Declare input, output and auxiliary file tags. <syntaxhighlight lang="fortran"> | |
− | <fortran> | ||
type(file) :: Fw,Fv,Fr,Fo ! I/O files | type(file) :: Fw,Fv,Fr,Fo ! I/O files | ||
− | </ | + | </syntaxhighlight> |
− | + | #Declare RSF cube axes: <tt>at</tt> time axis, <tt>ax</tt> space axis, <tt>az</tt> depth axis. <syntaxhighlight lang="fortran"> | |
− | <fortran> | ||
type(axa) :: at,az,ax ! cube axes | type(axa) :: at,az,ax ! cube axes | ||
− | </ | + | </syntaxhighlight> |
− | + | #Declare multi-dimensional arrays for input, output and computations. <syntaxhighlight lang="fortran"> | |
− | <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> |
− | + | #Open files for input/output. <syntaxhighlight lang="fortran"> | |
− | <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"> | |
− | <fortran> | + | call iaxa(Fw,at,1); nt = at%n; dt = at%d |
− | call iaxa(Fw,at,1); call iaxa(Fv,az,1); call iaxa(Fv,ax,2) | + | call iaxa(Fv,az,1); nz = az%n; dz = az%d |
− | call oaxa(Fo,az,1) | + | call iaxa(Fv,ax,2); nx = ax%n; dx = ax%d |
− | </ | + | |
− | + | call oaxa(Fo,az,1) | |
− | <fortran> | + | call oaxa(Fo,ax,2) |
− | allocate(ww( | + | call oaxa(Fo,at,3) |
− | allocate(vv( | + | </syntaxhighlight> |
− | allocate(rr( | + | #Allocate arrays and read wavelet, velocity and reflectivity. <syntaxhighlight lang="fortran"> |
− | </ | + | allocate(ww(nt)); call rsf_read(Fw,ww) |
− | + | allocate(vv(nz,nx)); call rsf_read(Fv,vv) | |
− | <fortran> | + | allocate(rr(nz,nx)); call rsf_read(Fr,rr) |
− | allocate(um( | + | </syntaxhighlight> |
− | allocate(uo( | + | #Allocate temporary arrays. <syntaxhighlight lang="fortran"> |
− | allocate(up( | + | allocate(um(nz,nx)); um=0. |
− | allocate(ud( | + | allocate(uo(nz,nx)); uo=0. |
− | </ | + | allocate(up(nz,nx)); up=0. |
− | + | allocate(ud(nz,nx)); ud=0. | |
− | <fortran> | + | </syntaxhighlight> |
− | do it=1, | + | #Loop over time. <syntaxhighlight lang="fortran"> |
− | </ | + | do it=1,nt |
− | + | </syntaxhighlight> | |
− | <fortran> | + | #Compute Laplacian: <math>\Delta U</math>. <syntaxhighlight lang="fortran"> |
− | + | 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 | |
− | + | </syntaxhighlight> | |
− | + | #Inject source wavelet: <math>\left[ \Delta U - f(t) \right]</math> <syntaxhighlight lang="fortran"> | |
− | |||
− | |||
− | </ | ||
− | |||
− | <fortran> | ||
ud = ud - ww(it) * rr | ud = ud - ww(it) * rr | ||
− | </ | + | </syntaxhighlight> |
− | + | #Scale by velocity: <math>\left[ \Delta U - f(t) \right] v^2</math> <syntaxhighlight lang="fortran"> | |
− | <fortran> | ||
ud= ud *vv*vv | 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="fortran"> | |
− | <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> | ||
+ | |||
+ | [[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 19:29, 5 May 2021
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
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
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);
}
- 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 */
- Declare RSF cube axes: at time axis, ax space axis, az depth axis. \tiny
sf_axis at,az,ax; /* cube axes */
- Declare multi-dimensional arrays for input, output and computations.
float *ww,**vv,**rr; /* I/O arrays*/
- Open files for input/output.
Fw = sf_input ("in" ); Fo = sf_output("out"); Fv = sf_input ("vel"); Fr = sf_input ("ref");
- 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);
- 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);
- Allocate temporary arrays.
um=sf_floatalloc2(nz,nx); uo=sf_floatalloc2(nz,nx); up=sf_floatalloc2(nz,nx); ud=sf_floatalloc2(nz,nx);
- Loop over time.
for (it=0; it<nt; it++) {
- 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; } }
- Inject source 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]; } }
\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);
}
- 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);
- 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);
- 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;
- 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;
- Initialize multidimensional valarray index counter (of type VAI).
VAI k(nz,nx);
- Loop over time.
for (int it=0; it<nt; it++) {
- 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; } }
- Inject source wavelet:
ud -= ww[it] * rr;
- Scale by velocity:
ud *= vv*vv;
- 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
- Declare input, output and auxiliary file tags.
type(file) :: Fw,Fv,Fr,Fo ! I/O files
- Declare RSF cube axes: at time axis, ax space axis, az depth axis.
type(axa) :: at,az,ax ! cube axes
- 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
- Open files for input/output.
Fw=rsf_input ("in") Fv=rsf_input ("vel") Fr=rsf_input ("ref") Fo=rsf_output("out")
- 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)
- 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)
- 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.
- Loop over time.
do it=1,nt
- 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
- Inject source wavelet:
ud = ud - ww(it) * rr
- Scale by velocity:
ud= ud *vv*vv
- Time step: Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle U_{i+1} = \left[ \Delta U -f(t) \right] v^2 \Delta t^2 + 2 U_{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)
- Open files for input/output.
Fw=m8r.Input() Fv=m8r.Input ("vel") Fr=m8r.Input ("ref") Fo=m8r.Output()
- 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)
- 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)
- 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')
- Loop over time.
for it in range(nt):
- 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
- Inject source 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
References
- ↑ "Hello world" of seismic imaging.