Merge branch 'master' of https://github.com/JanThorbecke/OpenSource

fdelmodc/Makefile

@@ -34,6 +34,7 @@ SRCC	= $(PRG).c \
 		getWaveletInfo.c  \
 		getModelInfo.c  \
 		applySource.c  \
+		applyRupture.c \
 		getRecTimes.c  \
 		getBeamTimes.c  \
 		writeSnapTimes.c  \

fdelmodc/applyRupture.c

+#include<stdlib.h>
+#include<stdio.h>
+#include<math.h>
+#include<assert.h>
+#include"fdelmodc.h"
+
+#define NINT(x) ((int)((x)>0.0?(x)+0.5:(x)-0.5))
+
+float stresscalc(float Tu, float Ts, float D, float dD, float D0, float V0);
+void slipcalc(float *D, float *dD, float *vx, float dx, int izsrc, int ixsrc1, int ixsrc2, int n1);
+
+void applyRupture(modPar mod, int itime, float *vx, float *txz, int verbose)
+{
+
+    float           *D, *dD, Tout;
+    int             irup;
+    static float    zrup, xrup1, xrup2, dT, T0, Tu, Ts, D0, V0, dx, trup;
+    static int      izrup, ixrup1, ixrup2, itrup, nrup, n1, first=1;
+
+    if (first==1) {
+        if(!getparfloat("zrup", &zrup)) zrup = 0.0;
+        if(!getparfloat("xrup1", &xrup1)) xrup1 = 0.0;
+        if(!getparfloat("xrup2", &xrup2)) xrup2 = 0.0;
+        if(!getparfloat("trup", &trup)) trup = 0.0;
+        if(!getparfloat("V0", &V0)) V0 = 0.0;
+        if(!getparfloat("D0", &D0)) D0 = 0.0;
+        if(!getparfloat("Tu", &Tu)) Tu = 0.0;
+        if(!getparfloat("Ts", &Ts)) Ts = 0.0;
+        if(!getparfloat("dT", &dT)) dT = 0.0;
+
+        dx = mod.dx;
+        n1 = mod.naz;
+
+        ixrup1 = NINT(xrup1/dx);
+        ixrup2 = NINT(xrup2/dx);
+        itrup = NINT(trup/mod.dt);
+        nrup = ixrup2 - ixrup1;
+        first = 0;
+    }
+
+    if (itime>itrup) {
+        D   = (float *)calloc(nrup,sizeof(float));
+        dD  = (float *)calloc(nrup,sizeof(float));
+
+        slipcalc(D, dD, vx, dx, izrup, ixrup1, ixrup2, n1);
+
+        for (irup=ixrup1; irup<ixrup2; irup++){
+            vmess("D[%d]=%.3e",D[irup-ixrup1]);
+            if (fabs(dD[irup-ixrup1])>0.001) {
+                Tout = stresscalc(Tu, Ts, D[irup-ixrup1], dD[irup-ixrup1], D0, V0);
+                txz[irup*n1+izrup-1] += Tout - T0;
+                txz[irup*n1+izrup] += Tout - T0;
+            }
+            else {
+                vx[irup*n1+izrup] = 0.0;
+            }
+        }
+
+        free(D); free(dD);
+    }
+    else if (itime==itrup) {
+        for (irup=ixrup1; irup<ixrup2; irup++){
+            txz[irup*n1+izrup-1] += dT;
+            txz[irup*n1+izrup] += dT;
+        }
+    }
+
+}
+
+float stresscalc(float Tu, float Ts, float D, float dD, float D0, float V0)
+{
+	float   T1, T2, Tout;
+
+    T1 = Tu*(1-D/D0);
+    T2 = Ts*(V0/(V0+dD));
+
+    if (fabs(T1) > fabs(T2)) Tout=T1;
+    else Tout = T2;
+	
+    return Tout;
+}
+
+void slipcalc(float *D, float *dD, float *vx, float dx, int izsrc, int ixsrc1, int ixsrc2, int n1)
+{
+    int     ix;
+
+    for (ix=ixsrc1; ix<ixsrc2; ix++) {
+        dD[ix-ixsrc1] = vx[ix*n1+izsrc+1] - vx[ix*n1+izsrc-1];
+    }
+
+    for (ix=ixsrc1+1; ix<ixsrc2; ix++) {
+        D[ix-ixsrc1] = (dx/2.0)*(dD[ix-ixsrc1-1]+dD[ix-ixsrc1]);
+    }
+    D[0] = (dx/2.0)*dD[0];
+
+}
\ No newline at end of file

fdelmodc/elastic4.c

@@ -17,6 +17,8 @@ int boundariesP(modPar mod, bndPar bnd, float *vx, float *vz, float *tzz, float
 
 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);
 
+void applyRupture(modPar mod, int itime, float *vx, float *txz, int verbose);
+
 int elastic4(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)
@@ -103,7 +105,7 @@ float *vz, float *tzz, float *txx, float *txz, float *rox, float *roz, float
 	}
 
 	/* Add force source */
-	if (src.type > 5) {
+	if (src.type > 5 && src.type!=20) {
 		 applySource(mod, src, wav, bnd, itime, ixsrc, izsrc, vx, vz, tzz, txx, txz, rox, roz, l2m, src_nwav, verbose);
 	}
 
@@ -144,6 +146,11 @@ float *vz, float *tzz, float *txx, float *txz, float *rox, float *roz, float
 	if (src.type < 6) {
 		 applySource(mod, src, wav, bnd, itime, ixsrc, izsrc, vx, vz, tzz, txx, txz, rox, roz, l2m, src_nwav, verbose);
 	}
+
+	/* Add rupture */
+	if (src.type==20) {
+		applyRupture(mod, itime, vx, txz, verbose);
+	}
     
 	/* check if there are sources placed on the boundaries */
     storeSourceOnSurface(mod, src, bnd, ixsrc, izsrc, vx, vz, tzz, txx, txz, verbose);

marchenko3D/zfpmar.h

+#ifndef ZFPMAR_H
+#define ZFPMAR_H
+#define MARBYTES        48
+#define TOPBYTES        80
+
+/* Size of the different types
+long    -   8 bytes
+int     -   4 bytes
+char    -   1 byte
+float   -   4 bytes
+*/
+
+/* TYPEDEFS */
+typedef struct {	/* zfpmar - headers for the compression and decompression of data for marchenko */
+
+    long nx; /* Number of samples in x-direction */
+
+    long ny; /* Number of samples in y-direction */
+
+    int sx; /* Source coordinate in x-direction */
+    
+    int sy; /* Source coordinate in y-direction */
+    
+    int sz; /* Source coordinate in z-direction */
+    
+    int gx; /* Receiver coordinate in x-direction */
+    
+    long gy; /* Receiver coordinate in y-direction */
+
+    long compsize; /* Size of the compressed data */
+
+} zfpmar;
+
+typedef struct {	/* zfptop - headers for the compression and decompression of data for marchenko */
+    
+    long ndim; /* Number of dimension is between 1 and 4 */
+
+    long nz; /* Number of samples in z-direction */
+
+    long ns; /* Number of shots */
+
+    long nt; /* Number of time samples */
+    
+    float dx; /* Sampling distance in x-direction */
+    
+    float dy; /* Sampling distance in x-direction */
+    
+    float dz; /* Sampling distance in x-direction */
+
+    float scale; /* Scaling of the coordinates and the sampling */
+
+    double tolerance; /* Set the tolerance of the zfp (de)compression */
+
+    float fmin; /* Minimum frequency of the signal */
+
+    float fmax; /* Maximum frequency of the singal */
+
+    float fx; /* First location of receiver in x-direction */
+
+    float fy; /* First location of receiver in y-direction */
+
+    float fz; /* First location of receiver in y-direction */
+
+} zfptop;
+#endif
\ No newline at end of file