komatits
commented
7 years ago Reply from Jeroen Tromp and from Daniel Peter @danielpeter : (the problem is thus specific to the 2D case, no problem in 3D, thus closing the 3D issue).
On 8/28/17 8:37 AM, Jeroen Tromp wrote:
Hello all:
In 2D, the "deviator" of a second-order tensor T_{ij} should be defined as
T{ij}-1/2 tr(T{ij}) delta_{ij}
and, in general, in n dimensions
T{ij}-1/n tr(T{ij}) delta_{ij}
For SH I start with defining a wave polarized in the, say, y direction:
s(x,z) \hat{y}
The 2D gradient is \hat{x}\partial_x + \hat{z}\partial_z, so that the displacement gradient is
\partial_x s \hat{x}\hat{y}+\partial_z s \hat{z}\hat{y}
and the strain is
1/2 \partial_x s (\hat{x}\hat{y}+\hat{y}\hat{x})+1/2 \partial_z s (\hat{z}\hat{y}+\hat{y}\hat{z})
The trace of this strain tensor is zero.
The stress is then
\mu \partial_x s (\hat{x}\hat{y}+\hat{y}\hat{x})+\mu \partial_z s (\hat{z}\hat{y}+\hat{y}\hat{z})
and the divergence of this stress is
[\partial_x(\mu \partial_x s)+\partial_z(\mu\partial_z s)] \hat{y}
which is in the y direction. Thus the SH equation of motion reduces to
\rho \partial_t^2 s = \partial_x(\mu \partial_x s)+\partial_z(\mu\partial_z s)
It's now easy to see that the \mu kernel is
K_\mu = int \mu [(\partial_x s) (\partial_x s^\dagger) + (\partial_z s) (\partial_z s^\dagger)] dt
i.e., no TWO.
P.S. The P-SV case should be done in a similar fashion, but starting with the two-component displacement
s_x(x,z) \hat{x}+ s_z(x,z) \hat{z}
Best regards,
Jeroen
On 8/28/17 12:38 AM, Daniel B. Peter wrote:
hi Jeannot, hi all,
good question…
short answer: yes there might be a bug in the code. let me test it with the example and see if the bug is in the kernel line you point out or the strain computation.
long answer: SPECFEM2D uses a plane strain convention, i.e. the off-plane strain is zero. thus, the question for the shear sensitivity kernel K_mu in the P-SV case (and SH case) is then what is the mean strain and the traceless, deviatoric strain for this problem?
the K_mu kernel (as defined in Tromp et al. 2005) uses a formulation with deviatoric strains: D = 1/2 ( grad s + (grad s)T) - 1/3 (div s)
here the mean strain is defined as 1/3 (div s) for the 3D case.
for SPECFEM2D, I was trying to find some good derivations for P-SV and SH cases, but only found in Peter Shearer’s book (Introduction to Seismology) some useful hints how you would set up the problem. in principle, that derivation follows the 3D case but uses a vector which assumes that the plane is the x-z plane and that displacement is either (u,0,w) or (0,v,0) for P-SV and SH waves, respectively. for the P-SV case, it follows that strains e_yy, e_xy and e_zy are all zero, which is the plane strain convention. thus, the strain tensor could be described by a 2D tensor.
and that’s the confusing part, in 2D the mean strain would be 1/2 [e_xx + e_zz] while in the 3D case it would be 1/3 [e_xx + e_yy(=0) + e_zz]. that’s ugly. so, for a 2D problem, the P-SV displacement vector (u,0,w) and no y-components would shrink to a 2x2 tensor and one could also take 1/2 (div s) = 1/2 (d/dx u + d/dz w) as mean strain… (then the factor 1/3 in the code would need to be corrected to be 1/2).
now, given Hooke’s law for isotropic materials: sigma_ij = lambda * delta_ij e_kk + 2 mu e_ij, you notice that although the strain e_yy vanishes, the stress component sigma_yy (transverse stress) is non-zero (and equal to lambda(e_xx+e_zz)). this transverse stress is often neglected and the whole problem treated as 2D, which then also often requires a change to “effective” material properties.
however, i guess the mean strain needs a coherent convention with the mean stress. pressure is related to mean stress by p = - 1/3 tr(sigma) in 3D. for a plane strain convention and the mean stress, Dimitri pointed out (in compute_pressure.f90):
! to compute pressure in 2D in an elastic solid in the plane strain convention i.e. in the P-SV case, ! one still uses pressure = - trace(sigma) / 3 but taking into account the fact ! that the off-plane strain epsilon_zz is zero by definition of the plane strain convention ! but thus the off-plane stress sigma_zz is not equal to zero, ! one has instead: sigma_zz = lambda (epsilon_xx + epsilon_yy), thus ! sigma_ij = lambda delta_ij trace(epsilon) + 2 mu epsilon_ij ! = lambda (epsilon_xx + epsilon_yy) + 2 mu epsilon_ij ! sigma_xx = lambda (epsilon_xx + epsilon_yy) + 2 mu epsilon_xx ! sigma_yy = lambda (epsilon_xx + epsilon_yy) + 2 mu epsilon_yy ! sigma_zz = lambda (epsilon_xx + epsilon_yy) ! pressure = - trace(sigma) / 3 = - (lambda + 2*mu/3) (epsilon_xx + epsilon_yy)
and so the code uses a factor 1/3.
now going back to the kernel computations, what i see is that the 2D code does implement the kernel K_mu using deviatoric strains.
- for the P-SV case:
line 229: ! shear modulus kernel mu_kl(i,j,ispec) = mu_kl(i,j,ispec) - TWO mul mu_k(iglob) * deltat
is based on mu_k(..), which is calculated starting on line 161: ! P-SV waves dsxx = dux_dxl dsxz = HALF * (duz_dxl + dux_dzl) dszz = duz_dzl
b_dsxx = b_dux_dxl b_dsxz = HALF * (b_duz_dxl + b_dux_dzl) b_dszz = b_duz_dzl kappa_k(iglob) = (dsxx + dszz) * (b_dsxx + b_dszz) mu_k(iglob) = dsxx * b_dsxx + dszz * b_dszz + & 2._CUSTOM_REAL * dsxz * b_dsxz - 1._CUSTOM_REAL/3._CUSTOM_REAL * kappa_k(iglob)which looks like the implementation of the product of the two deviatoric strains for forward and adjoint (') fields: D' : D = [1/2 ( grad s' + (grad s')T) - 1/3 (div s’)I ] : [1/2 (grad s + (grad s)T) - 1/3 (div s)I ] however, the code line might be a bit off. anyway, this would be the formula taken for a 3D case where the mean strain is 1/3 (div s) and removed to obtain the traceless deviatoric strain.
regarding the kernel K_mu, as defined in e.g. Tromp et al. 2005: K_mu = - integral 2 mu D' : D
is implemented by: mu_kl(i,j,ispec) = mu_kl(i,j,ispec) - TWO mul mu_k(iglob) * deltat
which seems ok on line 229, if one uses a parameterization for (rho,mu,bulk) as in the 3D case. some recent changes in the 2D code have now introduced a new bulk parameter: kappa = lambda + mu instead of kappa = lambda + 2/3 mu as in the 3D case. this could affect the kernel expression as well, making this factor TWO superfluous. so, it needs some more checking...
- for the SH-case, starting on line 174: ! SH (membrane) waves mu_k(iglob) = dux_dxl b_dux_dxl + dux_dzl b_dux_dzl
means for an SH-displacement vector (0,v,0) in the x/z-plane the derivatives: mu_k = d/dx v' d/dx v + d/dz v’ d/dz v
here strain e would be only given by e_xy = e_yx = 1/2 (d/dx v) and e_yz = e_zy = 1/2 (d/dz v), with zero trace. thus, it looks like mu_k is again the product D’ : D. however, again the code line looks a bit off by missing the 1/2 factors.
well, well, i’ll try to check with the example you tested…
best wishes, daniel
See https://github.com/geodynamics/specfem2d/issues/816 .