Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

vectorize anisotropy functions #577

Merged
merged 1 commit into from
Mar 24, 2024
Merged

vectorize anisotropy functions #577

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
85 changes: 55 additions & 30 deletions burnman/classes/anisotropy.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -167,7 +167,7 @@ def isentropic_universal_elastic_anisotropy(self):
def isentropic_isotropic_poisson_ratio(self):
"""
:returns: The isotropic Poisson ratio (mu) [unitless].
A metric of the laterial response to loading.
A metric of the lateral response to loading.
:rtype: float
"""
return (
Expand All @@ -187,14 +187,11 @@ def christoffel_tensor(self, propagation_direction):
:rtype: float
"""
propagation_direction = unit_normalize(propagation_direction)
Tik = np.tensordot(
np.tensordot(
self.full_isentropic_stiffness_tensor,
propagation_direction,
axes=([1], [0]),
),
Tik = np.einsum(
"ijkl, ...j, ...l",
self.full_isentropic_stiffness_tensor,
propagation_direction,
propagation_direction,
axes=([2], [0]),
)
return Tik

Expand All @@ -205,8 +202,12 @@ def isentropic_linear_compressibility(self, direction):
:rtype: float
"""
direction = unit_normalize(direction)
Sijkk = np.einsum("ijkk", self.full_isentropic_compliance_tensor)
beta = Sijkk.dot(direction).dot(direction)
beta = np.einsum(
"ijkk, ...i, ...j",
self.full_isentropic_compliance_tensor,
direction,
direction,
)
return beta

def isentropic_youngs_modulus(self, direction):
Expand All @@ -216,8 +217,14 @@ def isentropic_youngs_modulus(self, direction):
:rtype: float
"""
direction = unit_normalize(direction)
Sijkl = self.full_isentropic_compliance_tensor
S = Sijkl.dot(direction).dot(direction).dot(direction).dot(direction)
S = np.einsum(
"ijkl, ...i, ...j, ...k, ...l",
self.full_isentropic_compliance_tensor,
direction,
direction,
direction,
direction,
)
return 1.0 / S

def isentropic_shear_modulus(self, plane_normal, shear_direction):
Expand All @@ -232,13 +239,16 @@ def isentropic_shear_modulus(self, plane_normal, shear_direction):
assert (
np.abs(plane_normal.dot(shear_direction)) < np.finfo(float).eps
), "plane_normal and shear_direction must be orthogonal"
Sijkl = self.full_isentropic_compliance_tensor
G = (
Sijkl.dot(shear_direction)
.dot(plane_normal)
.dot(shear_direction)
.dot(plane_normal)

G = np.einsum(
"ijkl, ...i, ...j, ...k, ...l",
self.full_isentropic_compliance_tensor,
shear_direction,
plane_normal,
shear_direction,
plane_normal,
)

return 0.25 / G

def isentropic_poissons_ratio(self, axial_direction, lateral_direction):
Expand All @@ -250,14 +260,28 @@ def isentropic_poissons_ratio(self, axial_direction, lateral_direction):

axial_direction = unit_normalize(axial_direction)
lateral_direction = unit_normalize(lateral_direction)
assert (
np.abs(axial_direction.dot(lateral_direction)) < np.finfo(float).eps
), "axial_direction and lateral_direction must be orthogonal"

Sijkl = self.full_isentropic_compliance_tensor
x = axial_direction
y = lateral_direction
nu = -(Sijkl.dot(y).dot(y).dot(x).dot(x) / Sijkl.dot(x).dot(x).dot(x).dot(x))
if len(axial_direction.shape) == 1:
assert (
np.abs(np.einsum("i, i", axial_direction, lateral_direction)) < 1.0e-12
), "axial_direction and lateral_direction must be orthogonal"
else:
dot = np.einsum("...i, ...i->...", axial_direction, lateral_direction)
assert np.all(
np.abs(dot) < 1.0e-12
), "axial_direction and lateral_direction must be orthogonal"

nu = -(
self.isentropic_youngs_modulus(axial_direction)
* np.einsum(
"ijkl, ...i, ...j, ...k, ...l",
self.full_isentropic_compliance_tensor,
lateral_direction,
lateral_direction,
axial_direction,
axial_direction,
)
)
return nu

def wave_velocities(self, propagation_direction):
Expand All @@ -267,15 +291,16 @@ def wave_velocities(self, propagation_direction):
:rtype: list, containing the wave speeds and directions
of particle motion relative to the stiffness tensor
"""
propagation_direction = unit_normalize(propagation_direction)

Tik = self.christoffel_tensor(propagation_direction)

eigenvalues, eigenvectors = np.linalg.eig(Tik)

idx = eigenvalues.argsort()[::-1]
eigenvalues = np.real(eigenvalues[idx])
eigenvectors = eigenvectors[:, idx]
# sort eigs by decreasing eigenvalue
idx = np.flip(eigenvalues.argsort(axis=-1), axis=-1)
vidx = np.expand_dims(idx, -1)
eigenvalues = np.take_along_axis(eigenvalues, idx, axis=-1)
eigenvectors = np.take_along_axis(eigenvectors, vidx, axis=-2)
eigenvalues = np.real(eigenvalues)
velocities = np.sqrt(eigenvalues / self.density)

return velocities, eigenvectors
Expand Down
7 changes: 4 additions & 3 deletions tests/test_anisotropy.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -41,13 +41,14 @@ def test_isotropic(self):
d1 = [1.0, 0.0, 0.0]
d2 = [0.0, 1.0, 0.0]

beta_100 = m.isentropic_linear_compressibility(direction=d1)
E_100 = m.isentropic_youngs_modulus(direction=d1)
ds = np.array([d1, d2])
beta_100 = m.isentropic_linear_compressibility(direction=ds)[0]
E_100 = m.isentropic_youngs_modulus(direction=ds)[0]
G_100_010 = m.isentropic_shear_modulus(plane_normal=d1, shear_direction=d2)
nu_100_010 = m.isentropic_poissons_ratio(
axial_direction=d1, lateral_direction=d2
)
wave_speeds, wave_directions = m.wave_velocities(propagation_direction=d1)
wave_speeds, _ = m.wave_velocities(propagation_direction=d1)
Vp, Vs1, Vs2 = wave_speeds

array_2 = [
Expand Down
Loading