New feature: LVFZ in FINITE=1 loading conditions.

[bidini]

1D faults with a low-velocity fault zone are now available for FINITE=1 loading conditions (See documentation).

[martijnende]

Nice job Ben. I ran the Python test suite to ensure that your changes do not affect other parts of the code (all tests passed), and I modified a few comments. I suggest we merge this pull request into release/2.0.1 instead of directly into master.

[jpampuero]

I would not merge this version until we know the kernel is correct. Ben, did you derive your new kernel by following the same derivation steps as section 5 of Cochard and Rice (1997 http://esag.harvard.edu/rice/182_Cochard_Rice_JMPS_97.pdf)? I actually don't know how to do that, because the derivation starts with the expression of the kernel in space domain (eq 38) but we only know the LVFZ kernel in wavenumber domain (its Fourier transform). If you have achieved such a derivation please share it for our review. I used a different method in a matlab script "fz_crack.m" you have already used. That script can be adapted to pre-compute the LVFZ kernel and store it in a file, in the same way as src/TabKernelFiniteFlt.m does for the homogeneous kernel. If you decide to use that method, we can discuss the details by email.

[bidini]

I haven't written a new kernel. My implementation simply takes the homogeneous-medium kernel of FINITE=1 that was already there and multiplies it by:

(1-D) cotanh(Hk + atanh(1-D)),

where D and H are parameters of the LVFZ and k is the wavenumber. If H = D = 0, my implementation does nothing.

[jpampuero]

Without a mathematical derivation, there is no guarantee your kernel is correct. So let's use my alternative method. In matlab: pad=32; k = StaticKernel(pad * N,3,H/(pad * L), 1-D); k = 0.5 * mu/(pad * L) * real(ifft(k)); k = [k(1:N) ; k(pad * N-N+1:pad * N)]; k = real(fft(k)); Note that k(x) is a real and even function, its Fourier transform is real and even, but I still use real(...) above to avoid roundoff errors. Export k(1:N) into a file. Then modify subroutine init_kernel_2D to read the file as: k%kernel(1) = 0d0 do i=1,k%nnfft/2-1 read(57, * ,end=100) k%kernel(2 * i+1) k%kernel(2 * i+2) = k%kernel(2 * i+1) enddo read(57, * ,end=100) k%kernel(2) ! Nyquist

Every time you change mu, D, H, L or N you will have to generate a new kernel file. You can avoid that by porting to fortran my matlab code above (that would be more user-friendly).

[martijnende]

Every time you change mu, D, H, L or N you will have to generate a new kernel file. You can avoid that by porting to fortran my matlab code above (that would be more user-friendly).

Generating the kernel file could be done either in the wrapper or in the FORTRAN code. I would prefer an implementation in FORTRAN, to avoid duplicate code (since we have two wrappers: MATLAB and Python), and prevent users from (accidentally) messing with it.

How can we test the implementation of the kernel? Are there analytical solutions or asymptotic approximations that we can incorporate in the test suite?

[jpampuero]

Ben and I have used the kernel I proposed in matlab before for a different purpose (to compute slip induced by imposing a constant stress drop over a finite fault length). The results seemed reasonable. For instance, they approached the expected homogeneous limit when H was very small or very large. We could use that as a test. There is no analytical test, but we could also validate it by comparison to a finite/spectral element simulation.

I recommend to first implement it by writing a kernel file via matlab. Once tested, we can move it inside the fortran code, computing the kernel directly in module fault_stress without writing it on a file.

[bidini]

My understanding is that the kernel from FINITE=1 is indeed the Fourier transform of Eq. (40) in Cochard & Rice (1997). The file kernel_I.tab provides a numerical solution to the term in brackets. This is a kernel for a finite fault in a homogeneous medium. The LVFZ kernel in Ampuero et al. (2002) can be broken into the kernel of a periodic fault in a homogeneous medium, K(k), times a term related to the LVFZ properties:

Kd(k) = K(k)(1-D)cotanh(H*k + atanh(1-D) )

My idea was that we could apply such composition to the FINITE=1 implementation. Hence, the LVFZ kernel would be the Fourier transform of Eq. (40) in Cochard & Rice (1997) (already implemented by someone else in FINITE=1) times the LVFZ term. @jpampuero : Is this the step you think is not rigorous enough?

Moving to your suggestion, would it be wrong if instead of using finite=3 in StaticKernel.m we use finite=1(i.e., to use the CR97 static kernel for a finite fault instead of the trim trick)? Both implementations seem to be quite similar in that case. I can do some numerical comparison of the kernels in Fourier space if it helps.

[bidini]

A comparison between both implementations, borrowing extensively from Pablo's MATLAB scripts. You will need a kernel_I file large enough to run it.

clear all, close all Number of elements in the fault N = 1024;

LVFZ parameters H = 5e-4; D = 0.9;

My implementation n = abs(fftshift([-N:N-1]')); load kernel_I K1 = pi*n .*I(n+(n==0))*(1-D).*coth(pi*n*H + atanh(1-D));

Pablo's suggestion (trim trick) pad = 32; % factor used to increase the periodic fault length M = pad*N; % increase fault length k = 2*pi*fftshift([-M/2:M/2-1]'); K2 = abs(k)*(1-D).*coth(H/pad*abs(k) + atanh(1-D));

K2 = 1/(pad) * real(ifft(K2)); K2 = [K2(1:N) ; K2(M-N+1:M)]; % trim the fault length K2 = real(fft(K2));

Fig. 1: I played with different LVFZ parameters and I don't see a difference between the two kernels.

[jpampuero]

The kernel in your figure looks very linear. This suggests that the non-linearity of the coth term is not significant in your example, or at the scale of your plot. Can you zoom in close to zero wavenumber to do a more critical test? If that works it would be great news.

Your implementation is based on a guess, not on a rigorous mathematical derivation. I read CR97 again and tried to introduce the LVFZ, but I did not find anything like your proposed kernel. If the test shows that your intuitive guess is correct, we can try harder to derive it.

[martijnende]

@bidini What is the status of this pull request? Is it ready for merging now?

[jpampuero]

Some time ago @bidini showed me a figure that confirms his kernel implementation is numerically accurate.