Skip to content

Commit

Permalink
minor improvements
Browse files Browse the repository at this point in the history
  • Loading branch information
@marcoalopez
marcoalopez committed Jul 4, 2024
1 parent 48a7239 commit 64db3af
Show file tree
Hide file tree
Showing 2 changed files with 47 additions and 41 deletions.
22 changes: 11 additions & 11 deletions src/christoffel.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -185,7 +185,7 @@ def _christoffel_matrix(wavevectors: np.ndarray, Cijkl: np.ndarray) -> np.ndarra
Notes
-----
The Christoffel matrix is calculated using the formula
M = k @ Cijkl @ k, where M is the Christoffel matrix, k is a
Mil = qj @ Cijkl @ qk, where M is the Christoffel matrix, q is a
wave vector, and Cijkl is the elastic tensor (stiffness matrix).
"""

Expand All @@ -197,15 +197,15 @@ def _christoffel_matrix(wavevectors: np.ndarray, Cijkl: np.ndarray) -> np.ndarra
n = wavevectors.shape[0]

# initialize array (pre-allocate)
Mij = np.zeros((n, 3, 3))
Mil = np.zeros((n, 3, 3))

for i in range(n):
Mij[i, :, :] = np.dot(wavevectors[i, :], np.dot(wavevectors[i, :], Cijkl))
Mil[i, :, :] = np.dot(wavevectors[i, :], np.dot(wavevectors[i, :], Cijkl))

return Mij
return Mil


def _calc_eigen(Mij: np.ndarray):
def _calc_eigen(Mil: np.ndarray):
"""Return the eigenvalues and eigenvectors of the Christoffel
matrix sorted from low to high. The eigenvalues are related to
primary (P) and secondary (S-fast, S-slow) wave speeds. The
Expand All @@ -228,15 +228,15 @@ def _calc_eigen(Mij: np.ndarray):
"""

# get the number of Christoffel matrices to proccess
n = Mij.shape[0]
n = Mil.shape[0]

# preallocate arrays
eigen_values = np.zeros((n, 3))
eigen_vectors = np.zeros((n, 3, 3))

# TO REIMPLEMENT ASSUMING A SHAPE (n, 3, 3).
for i in range(n):
eigen_values[i, :], eigen_vectors[i, :, :] = np.linalg.eigh(Mij[i, :, :])
eigen_values[i, :], eigen_vectors[i, :, :] = np.linalg.eigh(Mil[i, :, :])

return eigen_values, eigen_vectors

Expand Down Expand Up @@ -328,16 +328,16 @@ def _christoffel_gradient_matrix(
n = wavevectors.shape[0]

# initialize array (pre-allocate)
delta_Mij = np.zeros((n, 3, 3, 3))
M_ilk = np.zeros((n, 3, 3, 3))

# calculate the gradient matrices from the Christoffel matrices
for i in range(n):
delta_Mij_temp = np.dot(
Milk_temp = np.dot(
wavevectors[i, :], Cijkl + np.transpose(Cijkl, (0, 2, 1, 3))
)
delta_Mij[i, :, :, :] = np.transpose(delta_Mij_temp, (1, 0, 2))
M_ilk[i, :, :, :] = np.transpose(Milk_temp, (1, 0, 2))

return delta_Mij
return M_ilk


def _eigenvalue_derivatives(
Expand Down
66 changes: 36 additions & 30 deletions src/test_christoffel_stepbystep.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,9 +31,9 @@
"\n",
"## Input data\n",
"\n",
"The necessary input data to estimate the seismic properties in any direction of a given material using the Christoffel equation are the density and the stiffness tensor of the material.\n",
"The data required to estimate the seismic properties in any direction of a given material using the Christoffel equation are the density and the stiffness tensor of the material.\n",
"\n",
"The stiffness tensor $C$ is a fundamental property of a material that generalises Hooke's law in three dimensions, relating strains to stresses in the elastic regime. This is usually given (or abbreviated) in a 6x6 matrix ($C_{ij}$) using Voigt's notation with values given in $GPa$. Density is usually given in $g/cm^3$.\n",
"The stiffness tensor $C$ is a fundamental property of a material that generalises Hooke's law in three dimensions, relating strains to stresses in the elastic regime. This is usually given (or abbreviated) in a 6x6 matrix ($C_{ij}$) using Voigt notation with values given in $GPa$. The density is usually given in $g/cm^3$.\n",
"\n",
"> _Other material properties unrelated to seismic waves that can be derived from these data are: Young modulus, shear modulus, Poisson's ratio, linear compresibility and Debye temperature._"
]
Expand Down Expand Up @@ -61,11 +61,11 @@
"source": [
"## Step 0: create an array with wave vectors\n",
"\n",
"For most materials, the seismic properties are anisotropic, i.e. they depend on the direction. Accordingly, the Christoffel equation must be solved for each propagation direction of interest, and we must first define such a direction, denoted as wave vectors or $q$. Wave vectors can be defined in spherical coordinates (polar and azimuthal angles) or in Cartesian coordinates ($x, y, z$) normalised to the unit sphere (i.e. vector of size 1). Cartesian coordinates are used for calculations. In this case we will use the following procedure:\n",
"For most materials, the seismic properties are anisotropic, i.e. they depend on the direction. Accordingly, the Christoffel equation must be solved for each propagation direction of interest, and we must first define these directions, denoted as wave vectors or $\\vec{q}$. Wave vectors can be defined in spherical coordinates (polar and azimuthal angles) or in Cartesian coordinates ($x, y, z$) normalised to the unit sphere (i.e. vectors of size 1 or unit vectors). Cartesian coordinates are typically used for calculations. In this case we will use the following procedure:\n",
"\n",
"- Create equispaced wavevectors in spherical coordinates\n",
"- Convert from spherical to 3D Cartesian unit vectors\n",
"- Store the wavevectors in a single array with three columns ($x, y, z$) for calculations"
"- Store the wave vectors in a single array with three columns ($x, y, z$) for calculations"
]
},
{
Expand Down Expand Up @@ -171,10 +171,10 @@
"q_n \\, C_{inmj} \\, q_m = \\rho \\frac{\\omega^2}{k^2} s_j\n",
"$$\n",
"\n",
"TODO\n",
"The Christoffel matrix for a given direction is determined as follows:\n",
"\n",
"$$\n",
"M_{ij} = \\sum_{nm} q_n \\, C_{inmj} \\, q_m\n",
"M_{il} = \\sum_{nm} q_j \\, C_{ijkl} \\, q_k\n",
"$$"
]
},
Expand All @@ -195,9 +195,9 @@
}
],
"source": [
"Mij = ch._christoffel_matrix(wavevectors, Cijkl)\n",
"Mil = ch._christoffel_matrix(wavevectors, Cijkl)\n",
"\n",
"Mij.shape"
"Mil.shape"
]
},
{
Expand All @@ -207,7 +207,7 @@
"## Step 3: Normalize the Christoffel matrices\n",
"\n",
"$$\n",
"M_{ij(norm)} = \\frac{M_{ij}}{\\rho}\n",
"M_{il(norm)} = \\frac{M_{il}}{\\rho}\n",
"$$"
]
},
Expand All @@ -218,7 +218,7 @@
"outputs": [],
"source": [
"scaling_factor = 1 / density\n",
"norm_Mij = Mij * scaling_factor"
"norm_Mil = Mil * scaling_factor"
]
},
{
Expand All @@ -244,7 +244,7 @@
],
"source": [
"# check\n",
"Mij[0], norm_Mij[0]"
"Mil[0], norm_Mil[0]"
]
},
{
Expand Down Expand Up @@ -273,7 +273,7 @@
}
],
"source": [
"eigenvalues, eigenvectors = ch._calc_eigen(norm_Mij)\n",
"eigenvalues, eigenvectors = ch._calc_eigen(norm_Mil)\n",
"\n",
"eigenvalues.shape, eigenvectors.shape"
]
Expand All @@ -284,21 +284,21 @@
"source": [
"## Step 5: Calculate phase velocities ($v_p$)\n",
"\n",
"The eigenvalues $\\lambda_1$, $\\lambda_2$ and $\\lambda_3$ of the Christoffel tensor $M_{ij}$ are always real, positive and related to the phase wave velocities $V_p$, $V_{s1}$, $V_{s2}$ propagating in the direction $\\vec{n}$ by the formulae\n",
"**Phase velocities** represents the speed at which individual wave (like the crest or trough) propagate through a medium. It's the velocity of a single frequency component within a wave packet. Pphase velocity depends on the material (rock type) and the wave's frequency due to dispersion in Earth materials. It’s defined as the ratio of the wave’s frequency ($\\omega$) to its wavelength ($\\kappa$).\n",
"\n",
"The eigenvalues $\\lambda_1$, $\\lambda_2$ and $\\lambda_3$ of the Christoffel tensor ($M_{il}$) are always real, positive and relate to the phase wave velocities $V_p$, $V_{s1}$, $V_{s2}$ propagating in the direction $\\vec{n}$ by the formulae\n",
"\n",
"$$\n",
"V_p = \\sqrt{ \\frac{\\lambda_1}{\\rho} }, \\quad V_{s1} = \\sqrt{ \\frac{\\lambda_2}{\\rho} }, \\quad V_{s2} = \\sqrt{ \\frac{\\lambda_3}{\\rho} }\n",
"$$\n",
"\n",
"where $\\rho$ denotes the material density. Note that if the Christoffel tensor is already normalized to the density this is just\n",
"where $\\rho$ is the material density. Note that if the Christoffel tensor is already normalized to the density so that\n",
"\n",
"$$\n",
"v_p = \\sqrt{\\lambda_n}, \\quad \\lambda_n = v^2_p\n",
"$$\n",
"\n",
"where $v_p$ denotes phase velocities and $\\lambda_n$ the various eigenvalues. Due to this, the output will be an array of shape $(n, 3)$, where $n$ represents the number of wavevectors (orientations) considered.\n",
"\n",
"> TODO -> define phase velocities"
"where $v_p$ denotes the phase velocities and $\\lambda_n$ the different eigenvalues. The output will be an array of shape $(n, 3)$, where $n$ represents the number of wave vectors (orientations) considered."
]
},
{
Expand Down Expand Up @@ -351,20 +351,25 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"## Towards group velocities ($v_g$)\n",
"## Determining the group velocities ($v_g$)\n",
"\n",
"> TODO -> Explain group velocities\n",
"**Group velocity** is the speed at which a wave packet (the \"entire energy pulse\") travels through a medium. It represents the velocity of the entire wave group, comprising multiple frequencies travelling together. In dispersive media, phase and group velocities can differ. In seismology, group velocities are often more relevant than phase velocities, as they represent the speed of energy propagation and are used for predicting the arrival times of seismic phases. Group velocity is defined as the derivative of the wave’s angular frequency with respect to its wave number and is related to phase velocity as follows:\n",
"\n",
"$$\n",
"v_g = \\vec{\\varDelta} v_p = \\vec{\\varDelta} \\sqrt{\\lambda} = \\frac{\\vec{\\varDelta} \\lambda}{2 \\sqrt{\\lambda}} = \\frac{\\vec{\\varDelta} \\lambda}{2 v_p}\n",
"v_g = \\frac{\\partial \\omega}{\\partial \\kappa} = \\vec{\\varDelta} v_p\n",
"$$\n",
"\n",
"In order to estimate the group velocities, it is first necessary to estimate the derivatives (gradients) of the Christoffel matrix and the eigenvalues for each direction.\n",
"\n",
"## Step 6: Calculate the derivatives of the Christoffel matrices ($∇M_{ij}$)\n",
"## Step 6: Calculate the derivatives of the Christoffel matrices ($∇M_{ilk}$)\n",
"\n",
"The gradient of the Christoffel matrix for a given direction is determined as follows (e.g. Jaeken and Cottenier, 2016):\n",
"\n",
"$$\n",
"\\frac{\\partial M_{ij}}{\\partial q_k} = \\sum_m (C_{ikmj} + C_{imkj}) q_m\n",
"$$"
"∇M_{ilk} = \\frac{\\partial M_{il}}{\\partial q_k} = \\sum_m (C_{ikmj} + C_{imkj}) q_m\n",
"$$\n",
"\n",
"For each direction, the result is reshaped to a 3D matrix with shape (3, 3, 3) resuting in an array of shape (n, 3, 3, 3)."
]
},
{
Expand All @@ -384,21 +389,21 @@
}
],
"source": [
"dMijk = ch._christoffel_gradient_matrix(wavevectors, Cijkl)\n",
"dMilk = ch._christoffel_gradient_matrix(wavevectors, Cijkl)\n",
"\n",
"dMijk.shape"
"dMilk.shape"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Step 7: Calculate the derivative of the eigenvalues\n",
"## Step 7: Calculate the derivative of the eigenvalues $\\vec{\\varDelta} \\lambda$\n",
"\n",
"The first derivative of an eigenvalue $\\lambda$ of a matrix $M$ is given by:\n",
"\n",
"$$\n",
"\\frac{\\partial \\lambda_i}{\\partial q_k} = \\hat{s}_i \\cdot \\frac{\\partial M_{ij}}{\\partial q_k} \\cdot \\hat{s}_i\n",
"\\vec{\\varDelta} \\lambda = \\frac{\\partial \\lambda_i}{\\partial q_k} = \\hat{s}_i \\cdot \\frac{\\partial M_{ij}}{\\partial q_k} \\cdot \\hat{s}_i\n",
"$$\n",
"\n",
"where $\\hat{s}_i$ is the normalized eigenvector corresponding to eigenvalue $\\lambda_i$."
Expand All @@ -421,7 +426,7 @@
}
],
"source": [
"eigenvalue_gradients = ch._eigenvalue_derivatives(eigenvectors, dMijk)\n",
"eigenvalue_gradients = ch._eigenvalue_derivatives(eigenvectors, dMilk)\n",
"\n",
"eigenvalue_gradients.shape"
]
Expand Down Expand Up @@ -518,6 +523,7 @@
}
],
"source": [
"# display the first ten directions\n",
"group_velocity_magnitudes[:10]"
]
},
Expand Down Expand Up @@ -995,7 +1001,7 @@
" # dv/dq = dv^2/dq / (2v)\n",
" group_vel[pol][cart] = grad_eig[pol][cart] / (2 * phase_vel[pol])\n",
" _group_abs[pol] = np.linalg.norm(group_vel[pol])\n",
" _group_dir[pol] = group_vel[pol] / _group_abs[pol]\n",
" _group_dir[pol] = group_vel[pol] / _group_abs[pol] # why?\n",
"\n",
" x = _group_dir[pol][0]\n",
" z = _group_dir[pol][2]\n",
Expand Down Expand Up @@ -1024,7 +1030,7 @@
"metadata": {},
"outputs": [],
"source": [
"res1, res2, res3, res4, res5, thetas, phis = set_group_velocity(phase_velocities[5], eigenvectors[5], dMijk[5], wavevectors[5])"
"res1, res2, res3, res4, res5, thetas, phis = set_group_velocity(phase_velocities[5], eigenvectors[5], dMilk[5], wavevectors[5])"
]
},
{
Expand Down

0 comments on commit 64db3af

Please sign in to comment.