Skip to content
Snippets Groups Projects
Commit d547efbf authored by Jan Thorbecke's avatar Jan Thorbecke
Browse files

fixed for general Make

parent 6acb6d35
No related branches found
No related tags found
No related merge requests found
#include<stdlib.h>
#include<stdio.h>
#include<math.h>
#include<assert.h>
#include"fdelmodc.h"
#define MAX(x,y) ((x) > (y) ? (x) : (y))
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 storeSourceOnSurface(modPar mod, srcPar src, bndPar bnd, int ixsrc, int izsrc, float *vx, float *vz, float *tzz, float *txx, float *txz, int verbose);
int reStoreSourceOnSurface(modPar mod, srcPar src, bndPar bnd, int ixsrc, int izsrc, float *vx, float *vz, float *tzz, float *txx, float *txz, int verbose);
int boundariesP(modPar mod, bndPar bnd, float *vx, float *vz, float *tzz, float *txx, float *txz, float *rox, float *roz, float *l2m, float *lam, float *mul, int itime, int verbose);
int boundariesV(modPar mod, bndPar bnd, float *vx, float *vz, float *tzz, float *txx, float *txz, float *rox, float *roz, float *l2m, float *lam, float *mul, int itime, int verbose);
int elastic4dc(modPar mod, srcPar src, wavPar wav, bndPar bnd, int itime, int ixsrc, int izsrc, float **src_nwav, float *vx,
float *vz, float *tzz, float *txx, float *txz, float *rox, float *roz, float
*l2m, float *lam, float *mul, int verbose)
{
/*********************************************************************
COMPUTATIONAL OVERVIEW OF THE 4th ORDER STAGGERED GRID:
The captial symbols T (=Txx,Tzz) Txz,Vx,Vz represent the actual grid
The indices ix,iz are related to the T grid, so the capital T
symbols represent the actual modelled grid.
one cel (iz,ix)
|
V extra column of vx,txz
|
------- V
| txz vz| txz vz txz vz txz vz txz vz txz vz txz
| |
| vx t | vx t vx t vx t vx t vx t vx
-------
txz vz txz vz txz vz txz vz txz vz txz vz txz
vx t vx T---Vx--T---Vx--T---Vx--T vx t vx
| | | | | | |
txz vz txz Vz--Txz-Vz--Txz-Vz Txz-Vz txz vz txz
| | | | | | |
vx t vx T---Vx--T---Vx--T---Vx--T vx t vx
| | | | | | |
txz vz txz Vz Txz-Vz Txz-Vz Txz-Vz txz vz txz
| | | | | | |
vx t vx T---Vx--T---Vx--T---Vx--T vx t vx
| | | | | | |
txz vz txz Vz Txz-Vz Txz-Vz Txz-Vz txz vz txz
| | | | | | |
vx t vx T---Vx--T---Vx--T---Vx--T vx t vx
txz vz txz vz txz vz txz vz txz vz txz vz txz
vx t vx t vx t vx t vx t vx t vx
txz vz txz vz txz vz txz vz txz vz txz vz txz <--|
|
extra row of txz/vz |
AUTHOR:
Jan Thorbecke (janth@xs4all.nl)
The Netherlands
***********************************************************************/
float c1, c2;
float dvx, dvz;
int ix, iz;
int n1;
c1 = 9.0/8.0;
c2 = -1.0/24.0;
n1 = mod.naz;
/* calculate vx for all grid points except on the virtual boundary*/
#pragma omp for private (ix, iz) nowait schedule(guided,1)
for (ix=mod.ioXx; ix<mod.ieXx; ix++) {
#pragma ivdep
for (iz=mod.ioXz; iz<mod.ieXz; iz++) {
vx[ix*n1+iz] -= rox[ix*n1+iz]*(
c1*(txx[ix*n1+iz] - txx[(ix-1)*n1+iz] +
txz[ix*n1+iz+1] - txz[ix*n1+iz]) +
c2*(txx[(ix+1)*n1+iz] - txx[(ix-2)*n1+iz] +
txz[ix*n1+iz+2] - txz[ix*n1+iz-1]) );
}
}
/* calculate vz for all grid points except on the virtual boundary */
#pragma omp for private (ix, iz) schedule(guided,1)
for (ix=mod.ioZx; ix<mod.ieZx; ix++) {
#pragma ivdep
for (iz=mod.ioZz; iz<mod.ieZz; iz++) {
vz[ix*n1+iz] -= roz[ix*n1+iz]*(
c1*(tzz[ix*n1+iz] - tzz[ix*n1+iz-1] +
txz[(ix+1)*n1+iz] - txz[ix*n1+iz]) +
c2*(tzz[ix*n1+iz+1] - tzz[ix*n1+iz-2] +
txz[(ix+2)*n1+iz] - txz[(ix-1)*n1+iz]) );
}
}
/* Add force source */
if (src.type > 5) {
applySource(mod, src, wav, bnd, itime, ixsrc, izsrc, vx, vz, tzz, txx, txz, rox, roz, l2m, src_nwav, verbose);
}
/* boundary condition clears velocities on boundaries */
boundariesP(mod, bnd, vx, vz, tzz, txx, txz, rox, roz, l2m, lam, mul, itime, verbose);
/* calculate Txx/tzz for all grid points except on the virtual boundary */
#pragma omp for private (ix, iz, dvx, dvz) nowait schedule(guided,1)
for (ix=mod.ioPx; ix<mod.iePx; ix++) {
#pragma ivdep
for (iz=mod.ioPz; iz<mod.iePz; iz++) {
dvx = c1*(vx[(ix+1)*n1+iz] - vx[ix*n1+iz]) +
c2*(vx[(ix+2)*n1+iz] - vx[(ix-1)*n1+iz]);
dvz = c1*(vz[ix*n1+iz+1] - vz[ix*n1+iz]) +
c2*(vz[ix*n1+iz+2] - vz[ix*n1+iz-1]);
txx[ix*n1+iz] -= l2m[ix*n1+iz]*dvx + l2m[ix*n1+iz]*dvz;
tzz[ix*n1+iz] -= l2m[ix*n1+iz]*dvz + l2m[ix*n1+iz]*dvx;
}
}
/* calculate Txz for all grid points except on the virtual boundary */
/*
#pragma omp for private (ix, iz) schedule(guided,1)
for (ix=mod.ioTx; ix<mod.ieTx; ix++) {
#pragma ivdep
for (iz=mod.ioTz; iz<mod.ieTz; iz++) {
txz[ix*n1+iz] -= mul[ix*n1+iz]*(
c1*(vx[ix*n1+iz] - vx[ix*n1+iz-1] +
vz[ix*n1+iz] - vz[(ix-1)*n1+iz]) +
c2*(vx[ix*n1+iz+1] - vx[ix*n1+iz-2] +
vz[(ix+1)*n1+iz] - vz[(ix-2)*n1+iz]) );
}
}
*/
/* Add stress source */
if (src.type < 6) {
applySource(mod, src, wav, bnd, itime, ixsrc, izsrc, vx, vz, tzz, txx, txz, rox, roz, l2m, src_nwav, verbose);
}
/* check if there are sources placed on the boundaries */
storeSourceOnSurface(mod, src, bnd, ixsrc, izsrc, vx, vz, tzz, txx, txz, verbose);
/* Free surface: calculate free surface conditions for stresses */
boundariesV(mod, bnd, vx, vz, tzz, txx, txz, rox, roz, l2m, lam, mul, itime, verbose);
/* restore source positions on the edge */
reStoreSourceOnSurface(mod, src, bnd, ixsrc, izsrc, vx, vz, tzz, txx, txz, verbose);
return 0;
}
......@@ -7,7 +7,7 @@ LIBS += -L$L -lgenfft -lm $(LIBSM)
#ALL: fmute marchenko marchenko2
ALL: fmute marchenko marchenko_tshift
ALL: fmute marchenko
SRCJ = fmute.c \
getFileInfo.c \
......@@ -37,23 +37,6 @@ SRCH = marchenko.c \
OBJJ = $(SRCJ:%.c=%.o)
SRCT = marchenko_tshift.c \
getFileInfo.c \
readData.c \
readShotData.c \
readTinvData.c \
applyMute_tshift.c \
writeData.c \
writeDataIter.c \
wallclock_time.c \
name_ext.c \
verbosepkg.c \
atopkge.c \
docpkge.c \
getpars.c
OBJT = $(SRCT:%.c=%.o)
fmute: $(OBJJ)
$(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o fmute $(OBJJ) $(LIBS)
......@@ -62,15 +45,12 @@ OBJH = $(SRCH:%.c=%.o)
marchenko: $(OBJH)
$(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o marchenko $(OBJH) $(LIBS)
marchenko_tshift: $(OBJT)
$(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o marchenko_tshift $(OBJT) $(LIBS)
OBJH2 = $(SRCH2:%.c=%.o)
marchenko2: $(OBJH2)
$(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o marchenko2 $(OBJH2) $(LIBS)
install: fmute marchenko marchenko_tshift
install: fmute marchenko
cp fmute $B
cp marchenko $B
......
......@@ -966,3 +966,4 @@ int main(int argc, char **argv)
return 0;
}
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment