Rıdvan Örsvuran

Einstein toplam gösterimi tekrar eden indislerin toplandığını varsayarak toplam işaretlerini yazma zorunluluğunu kaldırıyor.

Örneğin sismolojide moment tensör yoğunluğu elastik tensör (C), kayma vektörü (u) ve fay normali (v) cinsinden aşağıdaki gibi yazılabiliyor.

$$ M_{kl} = C_{ijkl} u_i v_j $$

Elastik tensör izotropik ortam için aşağıdaki gibi yazılabiliyor.

$$ C_{ijkl} = \lambda \delta_{ij} \delta_{kl}+ \mu(\delta_{ik}\delta_{jl} + \delta_{il}\delta_{jk}) $$

Bunları kodlamak isteğinizde ise yazılmayan toplam işaretleri size for döngüleri olarak geri dönüyor. İlk yazdığım versiyonda elastik tensörü fonksiyon olarak yazarak onun elemanlarını tek tek hesaplamaktan kurtulmuşum. Sıra moment tensöre geldiğinde ise 4 tane iç içe döngü yazmak zorunda kalmışım.

def generate_moment_tensor(u, v):
   """Generates and returns moment tensor from given u and v vectors
   for isotropic medium.

   Args:
       u (list): slip vector
       v (list): fault plane normal
   """
   mu = 1
   lmbda = 1

   M = np.zeros((3, 3))

   d = np.eye(3)  # kronecker delta

   C = lambda i, j, k, l: lmbda*d[i, j]*d[k, l]  \
       + mu*(d[i, k]*d[j, l] + d[i, l]*d[j, k])

   # M_kl = u_i v_j C_ijkl
   for k in range(3):
       for l in range(3):
           total = 0
           for i in range(3):
               for j in range(3):
                   total += u[i]*v[j]*C(i, j, k, l)
           M[k, l] = total

    return M

Sonradan numpy in einsum adlı bir fonksiyonu olduğunu öğrendim. Bu fonksiyonu kullanarak yazınca en azından görünüşte for döngülerini yazmak zorunda kalmıyorsunuz. Yalnız C için birkaç einsum kullanmak zorunda kaldım. Yazdığım şekilde çalışıyor ama bunu kronecker deltanın simetrisine borçluyum sanırım. Böyle birkaç toplam şeklinde yazılan değerler için belki fonksiyon yazmak daha doğru. C okunabilirliği de düştü diye düşünüyorum ama moment tensörünki fena durmuyor.

import numpy as np

def generate_moment_tensor_compact(u, v):
    """Generates and returns moment tensor from given u and v vectors
    for isotropic medium.

    Args:
        u (list): slip vector
        v (list): fault plane normal
    """
    mu = 1
    lmbda = 1
    d = eye(3) #kronecker delta
    C = lmbda*np.einsum('ij,kl', d, d) + mu*(np.einsum('ik,jl', d, d)+np.einsum('il,jk', d, d))

    return np.einsum('ijkl, i, j', C, u, v)