smesh.cc

二维射线追踪地震层析成像

/*
 * smesh.cc - slowness mesh implementation
 *
 * Jun Korenaga, MIT/WHOI
 * December 1998
 */

#include <iostream>
#include <fstream>
#include <cstdlib>
#include "smesh.h"
#include "interface.h"
#include "d_nr.h"

SlownessMesh2d::SlownessMesh2d(const char* fn)
    : eps(1e-6)
{
    ifstream s(fn);
    if (!s){
	cerr << "SlownessMesh2d::cannot open " << fn << "\n";
	exit(1);
    }

    s >> nx >> nz >> p_water >> p_air;
    if (p_water <= 0.0)
	cerr << "SlownessMesh2d::invalid water velocity\n";
    if (p_air <= 0.0)
	cerr << "SlownessMesh2d::invalid air velocity\n";
    p_water = 1.0/p_water;
    p_air = 1.0/p_air;
    nnodes = nx*nz;
    ncells = (nx-1)*(nz-1);

    pgrid.resize(nx,nz);
    vgrid.resize(nx,nz);
    ser_index.resize(nx,nz);
    node_index.resize(nnodes);
    cell_index.resize(ncells);
    index2cell.resize(nx-1,nz-1);
    xpos.resize(nx);
    topo.resize(nx);
    zpos.resize(nz);

    for (int i=1; i<=nx; i++) s >> xpos(i);
    for (int i=1; i<=nx; i++) s >> topo(i);
    for (int k=1; k<=nz; k++) s >> zpos(k);
    
    int N=1;
    vgrid = 0.0;
    for (int i=1; i<=nx; i++){
	for (int k=1; k<=nz; k++){
	    double vel;
	    s >> vel;
	    vgrid(i,k) = vel;
	    pgrid(i,k) = 1.0/vel;
	    ser_index(i,k) = N;
	    node_index(N).set(i,k);
	    N++;
	}
    }
    int icell=1;
    for (int i=1; i<nx; i++){
	for (int k=1; k<nz; k++){
	    index2cell(i,k) = icell;
	    cell_index(icell++) = ser_index(i,k);
	    
	}
    }

    // sanity check
    for (int i=2; i<=nx; i++){
	if (xpos(i)<=xpos(i-1)){
	    cerr << "SlownessMesh2d: illegal ordering of x nodes at i=" << i << '\n';
	    exit(1);
	}
    }
    for (int i=2; i<=nz; i++){
	if (zpos(i)<=zpos(i-1)){
	    cerr << "SlownessMesh2d: illegal ordering of z nodes at i=" << i << '\n';
	    exit(1);
	}
    }
    for (int i=1; i<=nx; i++){
	for (int k=1; k<=nz; k++){
	    if (vgrid(i,k)<=0){
		cerr << "SlownessMesh2d: non-positive velocity is found at (i,k)="
		     << i << "," << k << '\n';
		exit(1);
	    }
	}
    }

    dx_vec.resize(nx-1); rdx_vec.resize(nx-1);
    dz_vec.resize(nz-1); rdz_vec.resize(nz-1);
    b_vec.resize(nx-1);
    for (int i=1; i<nx; i++){
	dx_vec(i) = xpos(i+1)-xpos(i);
	rdx_vec(i) = 1.0/dx_vec(i);
	b_vec(i) = topo(i+1)-topo(i);
    }
    for (int k=1; k<nz; k++){
	dz_vec(k) = zpos(k+1)-zpos(k);
	rdz_vec(k) = 1.0/dz_vec(k);
    }

    Sm_H1.resize(4,4); Sm_H2.resize(4,4); Sm_V.resize(4,4);
    T_common.resize(4,4);
    commonGradientKernel();
    commonNormKernel();
}

void SlownessMesh2d::set(const Array1d<double>& u)
{
    if (u.size() != nnodes) error("SlownessMesh2d::set - size mismatch");

    int N=1;
    for (int i=1; i<=nx; i++){
	for (int k=1; k<=nz; k++){
	    pgrid(i,k) = u(N);
	    vgrid(i,k) = 1.0/pgrid(i,k);
	    N++;
	}
    }
}

void SlownessMesh2d::get(Array1d<double>& u) const
{
    if (u.size() != nnodes) error("SlownessMesh2d::set - size mismatch");

    int N=1;
    for (int i=1; i<=nx; i++){
	for (int k=1; k<=nz; k++){
	    u(N) = pgrid(i,k);
	    N++;
	}
    }
}

void SlownessMesh2d::vget(Array1d<double>& v) const
{
    if (v.size() != nnodes) error("SlownessMesh2d::set - size mismatch");

    int N=1;
    for (int i=1; i<=nx; i++){
	for (int k=1; k<=nz; k++){
	    v(N) = vgrid(i,k);
	    N++;
	}
    }
}

double SlownessMesh2d::xmin() const { return xpos.front(); }

double SlownessMesh2d::xmax() const { return xpos.back(); }

double SlownessMesh2d::zmin() const
{
    double tmin=1e30;
    for (int i=1; i<=topo.size(); i++)
	if (topo(i)<tmin) tmin = topo(i);

    return tmin+zpos.front();
}

double SlownessMesh2d::zmax() const
{
    double tmax=-1e30;
    for (int i=1; i<=topo.size(); i++)
	if (topo(i)>tmax) tmax = topo(i);
    return tmax+zpos.back();
}

const Index2d&
SlownessMesh2d::nodeIndex(int i) const { return node_index(i); }

int SlownessMesh2d::nodeIndex(int i, int k) const { return ser_index(i,k); }

Point2d SlownessMesh2d::nodePos(int i) const
{
    int node_i = node_index(i).i();
    int node_k = node_index(i).k();
    Point2d p(xpos(node_i),zpos(node_k)+topo(node_i));

    return p;
}

double SlownessMesh2d::calc_ttime(int nsrc, int nrcv) const
{
    if (nsrc<1 || nsrc>nnodes || nrcv<1 || nrcv>nnodes)
	error("SlownessMesh2d::calc_ttime - invalid input");
    Index2d src=node_index(nsrc);
    Index2d rcv=node_index(nrcv);
    int isrc=src.i(), ksrc=src.k();
    int ircv=rcv.i(), krcv=rcv.k();
    int di = ircv-isrc;
    int dk = krcv-ksrc;
    int dii = di > 0 ? 1 : -1;
    int dkk = dk > 0 ? 1 : -1;
    
    double ttime=0.0;
    if (di==0){ // ray path on a vertical grid
	int k1=ksrc, k2=ksrc+dkk;
	while (dk!=0){
	    double dz=abs(zpos(k1)-zpos(k2));
	    double pave = 0.5*(pgrid(isrc,k1)+pgrid(isrc,k2));
	    ttime += dz*pave;
	    dk -= dkk; k1 += dkk; k2 += dkk;
	}
    }else if (abs(di)==1){ // ray path within a single sheared column
	double dx = abs(xpos(isrc)-xpos(ircv));
	double pave1 = pgrid(isrc,ksrc);
	double dtopo = abs(topo(isrc)-topo(ircv));
	if (dk==0){
	    double dist = sqrt(dtopo*dtopo+dx*dx);
	    double pave2 = pgrid(ircv,ksrc);
	    double pave = 0.5*(pave1+pave2);
	    ttime += dist*pave;
	}else{
	    double ddx = dx/abs(dk);
	    double ddx2 = ddx*ddx;
	    double dratio = 1.0/abs(dk);

	    double ratio = dratio;
	    int k1=ksrc, k2=ksrc+dkk;
	    while (dk!=0){
		double dz=abs(zpos(k1)-zpos(k2))+dtopo;
		double dist = sqrt(dz*dz+ddx2);
		double pave2 = (1.0-ratio)*pgrid(isrc,k2)+ratio*pgrid(ircv,k2);
		double pave = 0.5*(pave1+pave2);
		ttime += dist*pave;
		dk -= dkk; k1 += dkk; k2 += dkk;
		ratio += dratio;
		pave1 = pave2; 
	    }
	}
    }else{
	double xs=xpos(isrc), xr=xpos(ircv);
	double dx = xr-xs;
	double zs=topo(isrc)+zpos(ksrc), zr=topo(ircv)+zpos(krcv);
	double dz = zr-zs;
	double x1=xs, z1=zs;
	double pave1 = pgrid(isrc,ksrc);
	int i2 = isrc+dii;
	double pave2;
	while (di!=0){
	    double x2 = xpos(i2);
	    double z2 = zs+dz*(x2-xs)/dx;
	    if (topo(i2)+zpos(1) > z2){ // above the model domain
		pave2 = p_water;
	    }else if (topo(i2)+zpos(nz) < z2){ // below the model domain
		pave2 = pgrid(i2,nz);
	    }else{
		// search for enclosing nodes
		int k1=-1, k2=-1;
		double zsrc = topo(i2)+zpos(ksrc);
		if (abs(zsrc-z2)<eps){
		    k1 = k2 = ksrc;
		}else if (zsrc<z2){ // search downward
		    if (ksrc<nz){
			for (int k=ksrc+1; k<=nz; k++){
			    double ztest = topo(i2)+zpos(k);
			    if (abs(ztest-z2)<eps){
				k1 = k2 = k;
				break;
			    }else if (ztest>z2){
				k2 = k;
				k1 = k-1;
				break;
			    }
			}
		    }
		}else{		// search upward
		    if (ksrc>1){
			for (int k=ksrc-1; k>-1; k--){
			    double ztest = topo(i2)+zpos(k);
			    if (abs(ztest-z2)<eps){
				k1 = k2 = k;
				break;
			    }else if (ztest<z2){
				k1 = k;
				k2 = k+1;
				break;
			    }
			}
		    }
		}
		if (k1==k2){
		    pave2 = pgrid(i2,k1);
		}else{
		    double ztest1=topo(i2)+zpos(k1);
		    double ztest2=topo(i2)+zpos(k2);
		    double dztest=ztest2-ztest1;
		    double ratio = (z2-ztest1)/dztest;
		    pave2 = ratio*pgrid(i2,k1)+(1.0-ratio)*pgrid(i2,k2);
		}
	    }
	    double ddx=x2-x1;
	    double ddz=z2-z1;
	    double dist = sqrt(ddx*ddx+ddz*ddz);
	    double pave = 0.5*(pave1+pave2);
	    ttime += dist*pave;

	    di-=dii; i2+=dii;
	    x1=x2; z1=z2;
	    pave1=pave2;
	}
    }

    if (ttime<0){
	cerr << "SlownessMesh2d::calc_ttime - negative traveltime encountered for ("
	     << nsrc << ", " << nrcv << ")\n";
	exit(1);
    }
    return ttime;
}

#ifdef notdef

double SlownessMesh2d::calc_ttime3(int nsrc, int nrcv) const
{
    const double tol_vdiff=1e-4;
    
    if (nsrc<1 || nsrc>nnodes || nrcv<1 || nrcv>nnodes)
	error("SlownessMesh2d::calc_ttime - invalid input");
    Index2d src=node_index(nsrc);
    Index2d rcv=node_index(nrcv);
    int isrc=src.i(), ksrc=src.k();
    int ircv=rcv.i(), krcv=rcv.k();
    int di = ircv-isrc;
    int dk = krcv-ksrc;
    int dii = di > 0 ? 1 : -1;
    int dkk = dk > 0 ? 1 : -1;
    
    double ttime=0.0;
    if (di==0){ // ray path on a vertical grid
	int k1=ksrc, k2=ksrc+dkk;
	while (dk!=0){
	    double dz=abs(zpos(k1)-zpos(k2));
	    double v1=vgrid(isrc,k1), v2=vgird(isrc,k2);
	    ttime += almost_exact_ttime(v1,v2,dz);
	    dk -= dkk; k1 += dkk; k2 += dkk;
	}
    }else if (abs(di)==1){ // ray path within a single sheared column
	double dx = abs(xpos(isrc)-xpos(ircv));
	double v1 = vgrid(isrc,ksrc);
	double dtopo = abs(topo(isrc)-topo(ircv));
	if (dk==0){
	    double dist = sqrt(dtopo*dtopo+dx*dx);
	    double v2 = vgrid(ircv,ksrc);
	    ttime += almost_exact_ttime(v1,v2,dist);
	}else{
	    double ddx = dx/abs(dk);
	    double ddx2 = ddx*ddx;
	    double dratio = 1.0/abs(dk);

	    double ratio = dratio;
	    int k1=ksrc, k2=ksrc+dkk;
	    while (dk!=0){
		double dz=abs(zpos(k1)-zpos(k2))+dtopo;
		double dist = sqrt(dz*dz+ddx2);
		double v2 = (1.0-ratio)*vgrid(isrc,k2)+ratio*vgrid(ircv,k2);
		ttime += almost_exact_ttime(v1,v2,dist);
		dk -= dkk; k1 += dkk; k2 += dkk;
		ratio += dratio;
	    }
	}
    }else{
	double xs=xpos(isrc), xr=xpos(ircv);
	double dx = xr-xs;
	double zs=topo(isrc)+zpos(ksrc), zr=topo(ircv)+zpos(krcv);
	double dz = zr-zs;
	double x1=xs, z1=zs;
	double v1 = vgrid(isrc,ksrc);
	double i2 = isrc+dii;
	double pave2;
	while (di!=0){
	    double x2 = xpos(i2);
	    double z2 = zs+dz*(x2-xs)/dx;
	    if (topo(i2)+zpos(1) > z2){ // above the model domain
		pave2 = p_water;
	    }else if (topo(i2)+zpos(nz) < z2){ // below the model domain
		pave2 = pgrid(i2,nz);
	    }else{
		// search for enclosing nodes
		int k1=-1, k2=-1;
		double zsrc = topo(i2)+zpos(ksrc);
		if (abs(zsrc-z2)<eps){
		    k1 = k2 = ksrc;
		}else if (zsrc<z2){ // search downward
		    if (ksrc<nz){
			for (int k=ksrc+1; k<=nz; k++){
			    double ztest = topo(i2)+zpos(k);
			    if (abs(ztest-z2)<eps){
				k1 = k2 = k;
				break;
			    }else if (ztest>z2){
				k2 = k;
				k1 = k-1;
				break;
			    }
			}
		    }
		}else{		// search upward
		    if (ksrc>1){
			for (int k=ksrc-1; k>-1; k--){
			    double ztest = topo(i2)+zpos(k);
			    if (abs(ztest-z2)<eps){
				k1 = k2 = k;
				break;
			    }else if (ztest<z2){
				k1 = k;
				k2 = k+1;
				break;
			    }
			}
		    }
		}
		if (k1==k2){
		    pave2 = pgrid(i2,k1);
		}else{
		    double ztest1=topo(i2)+zpos(k1);
		    double ztest2=topo(i2)+zpos(k2);
		    double dztest=ztest2-ztest1;
		    double ratio = (z2-ztest1)/dztest;
		    pave2 = ratio*pgrid(i2,k1)+(1.0-ratio)*pgrid(i2,k2);
		}
	    }
	    double ddx=x2-x1;
	    double ddz=z2-z1;
	    double dist = sqrt(ddx*ddx+ddz*ddz);
	    double pave = 0.5*(pave1+pave2);
	    ttime += dist*pave;

	    di-=dii; i2+=dii;
	    x1=x2; z1=z2; pave1=pave2;
	}
    }

    if (ttime<0){
	cerr << "SlownessMesh2d::calc_ttime - negative traveltime encountered for ("
	     << nsrc << ", " << nrcv << ")\n";
	exit(1);
    }
    return ttime;
}

double SlownessMesh2d::almost_exact_ttime(double v1, double v2,
					  double dpath) const
{
    const double tol_vdiff=1e-4;
    
    double dv = v2-v1;
    if (abs(dv)<tol_vdiff){
	return dpath*log(v2/v1)/dv;
    }else{
	return dpath/v1;
    }
}

#endif

void SlownessMesh2d::upperleft(const Point2d& p, Index2d& guess) const
{
    // note: this code guarantees that the final guess index is bounded by
    //       valid range (1...nx-1)(1...nz-1).
    int i_guess=guess.i(), k_guess=guess.k();

    if (xpos(i_guess)<=p.x()){
	if (i_guess < nx) i_guess++;
	while (xpos(i_guess)<=p.x() && i_guess<nx) i_guess++;
	i_guess--;
    }else{
	while (xpos(i_guess)>p.x() && i_guess>1) i_guess--;
    }
    double r = (p.x()-xpos(i_guess))*rdx_vec(i_guess);
    double dz = r*b_vec(i_guess)+topo(i_guess);

    if (zpos(k_guess)+dz<=p.y()){
	if (k_guess < nz) k_guess++;