-
Joeri Brackenhoff authoredJoeri Brackenhoff authored
applySource.c 10.94 KiB
#include<stdlib.h>
#include<stdio.h>
#include<math.h>
#include<assert.h>
#include"fdelmodc.h"
void vmess(char *fmt, ...);
#define c1 (9.0/8.0)
#define c2 (-1.0/24.0)
/*********************************************************************
*
* Add's the source amplitude(s) to the grid.
*
* For the acoustic schemes, the source-type must not be txx tzz or txz.
*
* AUTHOR:
* Jan Thorbecke (janth@xs4all.nl)
* The Netherlands
*
**********************************************************************/
int applySource(modPar mod, srcPar src, wavPar wav, bndPar bnd, int itime, int ixsrc, int izsrc, float *vx, float *vz, float *tzz, float *txx, float *txz, float *rox, float *roz, float *l2m, float **src_nwav, int verbose)
{
int is0, ibndz, ibndx;
int isrc, ix, iz, n1;
int id1, id2;
float src_ampl, time, scl, dt, sdx;
float Mxx, Mzz, Mxz, rake;
static int first=1;
if (src.type==6) {
ibndz = mod.ioXz;
ibndx = mod.ioXx;
}
else if (src.type==7) {
ibndz = mod.ioZz;
ibndx = mod.ioZx;
}
else if (src.type==2) {
ibndz = mod.ioTz;
ibndx = mod.ioTx;
if (bnd.lef==4 || bnd.lef==2) ibndx += bnd.ntap;
if (bnd.top==4 || bnd.top==2) ibndz += bnd.ntap;
}
else {
ibndz = mod.ioPz;
ibndx = mod.ioPx;
if (bnd.lef==4 || bnd.lef==2) ibndx += bnd.ntap;
if (bnd.top==4 || bnd.top==2) ibndz += bnd.ntap;
}
n1 = mod.naz;
dt = mod.dt;
sdx = 1.0/mod.dx;
/* special txz source activated? */
if ((bnd.top==1) && (src.type==2)) {
iz = izsrc + ibndz;
if (iz==ibndz) {
if (src.orient != 1) {
if (first) {
vmess("Only monopole Txz source allowed at surface. Reset to monopole");
first = 0;
}
src.orient=1;
}
} }
/*
* for plane wave sources the sources are placed
* around the central shot position
* the first source position has an offset in x of is0
*
* itime = 0 corresponds with time=0
* itime = 1 corresponds with time=dt
* src[0] (the first sample) corresponds with time = 0
*/
is0 = -1*floor((src.n-1)/2);
#pragma omp for private (isrc, src_ampl, ix, iz, time, id1, id2, scl)
for (isrc=0; isrc<src.n; isrc++) {
src_ampl=0.0;
/* calculate the source position */
if (src.random || src.multiwav) {
ix = src.x[isrc] + ibndx;
iz = src.z[isrc] + ibndz;
}
else { /* plane wave and point sources */
ix = ixsrc + ibndx + is0 + isrc;
iz = izsrc + ibndz;
}
time = itime*dt - src.tbeg[isrc];
id1 = floor(time/dt);
id2 = id1+1;
/* delay not reached or no samples left in source wavelet? */
if ( (time < 0.0) || ( (itime*dt) >= src.tend[isrc]) ) continue;
// fprintf(stderr,"isrc=%d ix=%d iz=%d src.x=%d src.z=%d\n", isrc, ix, iz, src.x[isrc], src.z[isrc]);
if (!src.multiwav) { /* only one wavelet for all sources */
src_ampl = src_nwav[0][id1]*(id2-time/dt) + src_nwav[0][id2]*(time/dt-id1);
}
else { /* multi-wavelet sources */
src_ampl = src_nwav[isrc][id1]*(id2-time/dt) + src_nwav[isrc][id2]*(time/dt-id1);
}
if (src_ampl==0.0) continue;
if ( ((ix-ibndx)<0) || ((ix-ibndx)>mod.nx) ) continue; /* source outside grid */
if (verbose>=4 && itime==0) {
vmess("Source %d positioned at grid ix=%d iz=%d",isrc, ix, iz);
}
/* cosine squared windowing to reduce edge effects on shot arrays */
if ( (src.n>1) && src.window) {
scl = 1.0;
if (isrc < src.window) {
scl = cos(0.5*M_PI*(src.window - isrc)/src.window);
}
else if (isrc > src.n-src.window+1) {
scl = cos(0.5*M_PI*(src.window - (src.n-isrc+1))/src.window);
}
src_ampl *= scl*scl;
}
/* source scaling factor to compensate for discretisation */
/* old amplitude setting does not obey reciprocity */
// src_ampl *= rox[ix*n1+iz]*l2m[ix*n1+iz]/(dt);
/* in older version added factor 2.0 to be compliant with defined Green's functions in Marchenko algorithm */
/* this is now set to 1.0 */
src_ampl *= (1.0/mod.dx)*l2m[ix*n1+iz];
if (verbose>5) { vmess("Source %d at grid [ix=%d,iz=%d] at itime %d has value %e",isrc, ix,iz, itime, src_ampl);
}
/* Force source */
if (src.type == 6) {
vx[ix*n1+iz] += src_ampl*rox[ix*n1+iz]/(l2m[ix*n1+iz]);
/* stable implementation from "Numerical Techniques for Conservation Laws with Source Terms" by Justin Hudson */
//vx[ix*n1+iz] = 0.5*(vx[(ix+1)*n1+iz]+vx[(ix-1)*n1+iz])+src_ampl*rox[ix*n1+iz]/(l2m[ix*n1+iz]);
}
else if (src.type == 7) {
vz[ix*n1+iz] += src_ampl*roz[ix*n1+iz]/(l2m[ix*n1+iz]);
/* stable implementation from "Numerical Techniques for Conservation Laws with Source Terms" by Justin Hudson */
/* stable implementation changes amplitude and more work is needed */
//vz[ix*n1+iz] = 0.5*(vz[ix*n1+iz-1]+vz[ix*n1+iz+1])+src_ampl*roz[ix*n1+iz]/(l2m[ix*n1+iz]);
//vz[ix*n1+iz] = 0.25*(vz[ix*n1+iz-2]+vz[ix*n1+iz-1]+vz[ix*n1+iz]+vz[ix*n1+iz+1])+src_ampl*roz[ix*n1+iz]/(l2m[ix*n1+iz]);
} /* src.type */
/* Stress source */
if (mod.ischeme <= 2) { /* Acoustic scheme */
/* Compressional source */
if (src.type == 1) {
if (src.orient != 1) src_ampl=src_ampl/mod.dx;
if (src.orient==1) { /* monopole */
tzz[ix*n1+iz] += src_ampl;
}
else if (src.orient==2) { /* dipole +/- */
tzz[ix*n1+iz] += src_ampl;
tzz[ix*n1+iz+1] -= src_ampl;
}
else if (src.orient==3) { /* dipole - + */
tzz[ix*n1+iz] += src_ampl;
tzz[(ix-1)*n1+iz] -= src_ampl;
}
else if (src.orient==4) { /* dipole +/0/- */
if (iz > ibndz)
tzz[ix*n1+iz-1]+= 0.5*src_ampl;
if (iz < mod.nz+ibndz-1)
tzz[ix*n1+iz+1] -= 0.5*src_ampl;
}
else if (src.orient==5) { /* dipole + - */
tzz[ix*n1+iz] += src_ampl;
tzz[(ix+1)*n1+iz] -= src_ampl;
}
}
}
else { /* Elastic scheme */
/* Compressional source */
if (src.type == 1) {
if (src.orient==1) { /* monopole */
txx[ix*n1+iz] += src_ampl;
tzz[ix*n1+iz] += src_ampl;
}
else if (src.orient==2) { /* dipole +/- */
txx[ix*n1+iz] += src_ampl;
tzz[ix*n1+iz] += src_ampl;
txx[ix*n1+iz+1] -= src_ampl;
tzz[ix*n1+iz+1] -= src_ampl;
}
else if (src.orient==3) { /* dipole - + */
txx[ix*n1+iz] += src_ampl;
tzz[ix*n1+iz] += src_ampl;
txx[(ix-1)*n1+iz] -= src_ampl;
tzz[(ix-1)*n1+iz] -= src_ampl;
}
else if (src.orient==4) { /* dipole +/0/- */
if (iz > ibndz) { txx[ix*n1+iz-1]+= 0.5*src_ampl;
tzz[ix*n1+iz-1]+= 0.5*src_ampl;
}
if (iz < mod.nz+ibndz-1) {
txx[ix*n1+iz+1] -= 0.5*src_ampl;
tzz[ix*n1+iz+1] -= 0.5*src_ampl;
}
}
else if (src.orient==5) { /* dipole + - */
txx[ix*n1+iz] += src_ampl;
tzz[ix*n1+iz] += src_ampl;
txx[(ix+1)*n1+iz] -= src_ampl;
tzz[(ix+1)*n1+iz] -= src_ampl;
}
}
else if (src.type == 2) {
/* Txz source */
if ((iz == ibndz) && bnd.top==1) {
txz[(ix-1)*n1+iz-1] += src_ampl;
txz[ix*n1+iz-1] += src_ampl;
}
else {
txz[ix*n1+iz] += src_ampl;
}
/* possible dipole orientations for a txz source */
if (src.orient == 2) { /* dipole +/- */
txz[ix*n1+iz+1] -= src_ampl;
}
else if (src.orient == 3) { /* dipole - + */
txz[(ix-1)*n1+iz] -= src_ampl;
}
else if (src.orient == 4) { /* dipole +/O/- */
/* correction: subtrace previous value to prevent z-1 values. */
txz[ix*n1+iz] -= 2.0*src_ampl;
txz[ix*n1+iz+1] += src_ampl;
}
else if (src.orient == 5) { /* dipole + - */
txz[(ix+1)*n1+iz] -= src_ampl;
}
}
/* Tzz source */
else if(src.type == 3) {
tzz[ix*n1+iz] += src_ampl;
}
/* Txx source */
else if(src.type == 4) {
txx[ix*n1+iz] += src_ampl;
}
/***********************************************************************
* pure potential shear S source (experimental)
* Curl S-pot = CURL(F) = dF_x/dz - dF_z/dx
***********************************************************************/
else if(src.type == 5) {
src_ampl = src_ampl*rox[ix*n1+iz]/(l2m[ix*n1+iz]);
if (src.orient == 3) src_ampl = -src_ampl;
/* first order derivatives */
vx[ix*n1+iz] += src_ampl*sdx;
vx[ix*n1+iz-1] -= src_ampl*sdx;
vz[ix*n1+iz] -= src_ampl*sdx;
vz[(ix-1)*n1+iz] += src_ampl*sdx;
/* second order derivatives */
/*
vx[ix*n1+iz] += c1*src_ampl*sdx;
vx[ix*n1+iz-1] -= c1*src_ampl*sdx;
vx[ix*n1+iz+1] += c2*src_ampl*sdx;
vx[ix*n1+iz-2] -= c2*src_ampl*sdx;
vz[ix*n1+iz] -= c1*src_ampl*sdx; vz[(ix-1)*n1+iz] += c1*src_ampl*sdx;
vz[(ix+1)*n1+iz] -= c2*src_ampl*sdx;
vz[(ix-2)*n1+iz] += c2*src_ampl*sdx;
*/
/* determine second position of dipole */
if (src.orient == 2) { /* dipole +/- vertical */
iz += 1;
vx[ix*n1+iz] -= src_ampl*sdx;
vx[ix*n1+iz-1] += src_ampl*sdx;
vz[ix*n1+iz] += src_ampl*sdx;
vz[(ix-1)*n1+iz] -= src_ampl*sdx;
}
else if (src.orient == 3) { /* dipole - + horizontal */
ix += 1;
vx[ix*n1+iz] -= src_ampl*sdx;
vx[ix*n1+iz-1] += src_ampl*sdx;
vz[ix*n1+iz] += src_ampl*sdx;
vz[(ix-1)*n1+iz] -= src_ampl*sdx;
}
}
/***********************************************************************
* pure potential pressure P source (experimental)
* Divergence P-pot = DIV(F) = dF_x/dx + dF_z/dz
***********************************************************************/
else if(src.type == 8) {
src_ampl = src_ampl*rox[ix*n1+iz]/(l2m[ix*n1+iz]);
if (src.orient == 3) src_ampl = -src_ampl;
vx[(ix+1)*n1+iz] += src_ampl*sdx;
vx[ix*n1+iz] -= src_ampl*sdx;
vz[ix*n1+iz+1] += src_ampl*sdx;
vz[ix*n1+iz] -= src_ampl*sdx;
/* determine second position of dipole */
if (src.orient == 2) { /* dipole +/- */
iz += 1;
vx[(ix+1)*n1+iz] -= src_ampl*sdx;
vx[ix*n1+iz] += src_ampl*sdx;
vz[ix*n1+iz+1] -= src_ampl*sdx;
vz[ix*n1+iz] += src_ampl*sdx;
}
else if (src.orient == 3) { /* dipole - + */
ix += 1;
vx[(ix+1)*n1+iz] -= src_ampl*sdx;
vx[ix*n1+iz] += src_ampl*sdx;
vz[ix*n1+iz+1] -= src_ampl*sdx;
vz[ix*n1+iz] += src_ampl*sdx;
}
}
else if(src.type == 9) {
rake = 0.5*M_PI;
Mxx = -1.0*(sin(src.dip)*cos(rake)*sin(2.0*src.strike)+sin(src.dip*2.0)*sin(rake)*sin(src.strike)*sin(src.strike));
Mxz = -1.0*(cos(src.dip)*cos(rake)*cos(src.strike)+cos(src.dip*2.0)*sin(rake)*sin(src.strike));
Mzz = sin(src.dip*2.0)*sin(rake);
txx[ix*n1+iz] -= Mxx*src_ampl;
tzz[ix*n1+iz] -= Mzz*src_ampl;
txz[ix*n1+iz] -= Mxz*src_ampl;
} /* src.type */
} /* ischeme */
} /* loop over isrc */
return 0;
}