From 88be3018d8a26ad832e96975ea0fcec000ed0a71 Mon Sep 17 00:00:00 2001 From: JBrackenhoff <J.A.Brackenhoff@tudelft.nl> Date: Tue, 30 Jun 2020 15:49:24 +0200 Subject: [PATCH] rupture --- fdelmodc/Makefile | 1 - fdelmodc/applyRupture.c | 96 ------------------------ fdelmodc/demo/fdelmodc_moment_tensor.scr | 2 +- fdelmodc/elastic4.c | 7 -- 4 files changed, 1 insertion(+), 105 deletions(-) delete mode 100644 fdelmodc/applyRupture.c diff --git a/fdelmodc/Makefile b/fdelmodc/Makefile index 37f8880..ed9746c 100644 --- a/fdelmodc/Makefile +++ b/fdelmodc/Makefile @@ -35,7 +35,6 @@ SRCC = $(PRG).c \ getWaveletInfo.c \ getModelInfo.c \ applySource.c \ - applyRupture.c \ getRecTimes.c \ getBeamTimes.c \ writeSnapTimes.c \ diff --git a/fdelmodc/applyRupture.c b/fdelmodc/applyRupture.c deleted file mode 100644 index 12212df..0000000 --- a/fdelmodc/applyRupture.c +++ /dev/null @@ -1,96 +0,0 @@ -#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 diff --git a/fdelmodc/demo/fdelmodc_moment_tensor.scr b/fdelmodc/demo/fdelmodc_moment_tensor.scr index 6fa9912..a0df514 100755 --- a/fdelmodc/demo/fdelmodc_moment_tensor.scr +++ b/fdelmodc/demo/fdelmodc_moment_tensor.scr @@ -22,7 +22,7 @@ export OMP_NUM_THREADS=1 ../fdelmodc \ file_cp=$filecp file_den=$filero \ ischeme=5 \ - src_type=11 Mxx=0.0 Mzz=0.0 Mxz=-1.0 tmod=2.0 \ + src_type=11 Mxx=1.0 Mzz=0.2 Mxz=-1.0 tmod=2.0 \ file_src=wavelet.su verbose=2 \ file_rcv=rec4S.su \ file_snap=snap4S0.su \ diff --git a/fdelmodc/elastic4.c b/fdelmodc/elastic4.c index 532f1cf..033ace2 100644 --- a/fdelmodc/elastic4.c +++ b/fdelmodc/elastic4.c @@ -17,8 +17,6 @@ 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) @@ -146,11 +144,6 @@ 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); -- GitLab