Repository: sympy/sympy

Issue:
Tensor: abstract tensor calculations (no components): calculating derivative wrt symbol within tensor expression does not work
Hey guys,

First of all: Sympy is a great Python package! Keep up the good work! 

I recently realized that the `tensor` module can be used for performing abstract tensor calculations without relying on any specific component representations. I played around a bit and tried to perform a derivative of a tensor expression by a specific tensor utilizing directional derivatives (sample code below). The directional derivative relies on inserting a small perturbation into the expression and calculate the derivative by a real variable, but this leads to an error. Since I am not a sympy power user, probably I do not understand how to perform the task properly. But maybe it is a bug. Please find my sample code below:
```Python
import sympy
from sympy.tensor.tensor import (TensorIndexType, 
                                 tensor_indices, 
				 tensorhead)
from sympy import Wild

# wild cards
 
w0 = Wild('w0')

# expansion parameter

tau = sympy.Symbol('tau')

Euclidean = TensorIndexType('Euclidean', dummy_fmt='E')
g = Euclidean.metric
a = tensorhead('a', [Euclidean])
b = tensorhead('b', [Euclidean])
i, j, k = tensor_indices('i j k', Euclidean)
db = tensorhead('{{\\Delta}b}', [Euclidean])

print(a(i)*b(-i))
# a(E_0)*b(-E_0)

# Calculate tensor derivative w.r.t b(-i) by using a directional derivative
# https://en.wikipedia.org/wiki/Tensor_derivative_(continuum_mechanics)
# 1. expand b(-i) by b(-i) -> b(-i) + tau*db(-i)

print((a(i)*b(-i)).replace(b(-w0), b(-w0) + tau*db(-w0)))
# a(E_0)*(b(-E_0) + tau*{{\Delta}b}(-E_0))

print(((a(i)*b(-i)).replace(b(-w0), b(-w0) + tau*db(-w0))).expand())
# a(E_0)*b(-E_0) + tau*a(E_0)*{{\Delta}b}(-E_0)

# 2. calculate derivative w.r.t expansion variable tau -> BOOM

((a(i)*b(-i)).replace(b(-w0), b(-w0) + tau*db(-w0))).expand().diff(tau)
# Traceback (most recent call last):
#  File "<stdin>", line 1, in <module>
#  File "/home/joha2/anaconda3/envs/sage/lib/python3.7/site-packages/sympy/core/expr.py", line 3095, in diff
#    return Derivative(self, *symbols, **assumptions)
#  File "/home/joha2/anaconda3/envs/sage/lib/python3.7/site-packages/sympy/core/function.py", line 1198, in __new__
#    it cannot be differentiated.''' % expr))
#ValueError: 
#Since there are no variables in the expression a(E_0)*b(-E_0) +
#tau*a(E_0)*{{\Delta}b}(-E_0), it cannot be differentiated.

# 3. would be: set `tau = 0` and isolate coefficient of db(-i) which 
# gives the tensor derivative of expression by b(-i)
# Should also work for more general expressions 
# and after canonicalization and contraction also for 
# expressions with more complicated monoterm symmetries
```

Also I realized that some of the code examples in the doc strings of the `tensor` module lead to errors. I will file an issue for that later.

Best regards
Johannes

PS: `sympy.__version__` shows 1.4


Discussion:

Patch:
diff --git a/sympy/tensor/tensor.py b/sympy/tensor/tensor.py
index 3d3e7fce0e..e07b600745 100644
--- a/sympy/tensor/tensor.py
+++ b/sympy/tensor/tensor.py
@@ -37,7 +37,7 @@
 from collections import defaultdict
 import operator
 import itertools
-from sympy import Rational, prod, Integer
+from sympy import Rational, prod, Integer, default_sort_key
 from sympy.combinatorics import Permutation
 from sympy.combinatorics.tensor_can import get_symmetric_group_sgs, \
     bsgs_direct_product, canonicalize, riemann_bsgs
@@ -2252,6 +2252,14 @@ def _check_add_Sum(self, expr, index_symbols):
             expr = Sum(expr, *sum_indices)
         return expr
 
+    def _expand_partial_derivative(self):
+        # simply delegate the _expand_partial_derivative() to
+        # its arguments to expand a possibly found PartialDerivative
+        return self.func(*[
+                    a._expand_partial_derivative()
+                    if isinstance(a, TensExpr) else a
+                    for a in self.args])
+
 
 class TensAdd(TensExpr, AssocOp):
     """
@@ -2298,6 +2306,7 @@ class TensAdd(TensExpr, AssocOp):
     def __new__(cls, *args, **kw_args):
         args = [_sympify(x) for x in args if x]
         args = TensAdd._tensAdd_flatten(args)
+        args.sort(key=default_sort_key)
         if not args:
             return S.Zero
         if len(args) == 1:
@@ -2545,7 +2554,6 @@ def _print(self):
         args = self.args
         for x in args:
             a.append(str(x))
-        a.sort()
         s = ' + '.join(a)
         s = s.replace('+ -', '- ')
         return s
@@ -2590,6 +2598,18 @@ def __iter__(self):
     def _eval_rewrite_as_Indexed(self, *args):
         return Add.fromiter(args)
 
+    def _eval_partial_derivative(self, s):
+        # Evaluation like Add
+        list_addends = []
+        for a in self.args:
+            if isinstance(a, TensExpr):
+                list_addends.append(a._eval_partial_derivative(s))
+            # do not call diff if s is no symbol
+            elif s._diff_wrt:
+                list_addends.append(a._eval_derivative(s))
+
+        return self.func(*list_addends)
+
 
 class Tensor(TensExpr):
     """
@@ -2997,6 +3017,50 @@ def _eval_rewrite_as_Indexed(self, tens, indices):
         expr = Indexed(tens.args[0], *index_symbols)
         return self._check_add_Sum(expr, index_symbols)
 
+    def _eval_partial_derivative(self, s):  # type: (Tensor) -> Expr
+
+        if not isinstance(s, Tensor):
+            return S.Zero
+        else:
+
+            # @a_i/@a_k = delta_i^k
+            # @a_i/@a^k = g_ij delta^j_k
+            # @a^i/@a^k = delta^i_k
+            # @a^i/@a_k = g^ij delta_j^k
+            # TODO: if there is no metric present, the derivative should be zero?
+
+            if self.head != s.head:
+                return S.Zero
+
+            # if heads are the same, provide delta and/or metric products
+            # for every free index pair in the appropriate tensor
+            # assumed that the free indices are in proper order
+            # A contravariante index in the derivative becomes covariant
+            # after performing the derivative and vice versa
+
+            kronecker_delta_list = [1]
+
+            # not guarantee a correct index order
+
+            for (count, (iself, iother)) in enumerate(zip(self.get_free_indices(), s.get_free_indices())):
+                if iself.tensor_index_type != iother.tensor_index_type:
+                    raise ValueError("index types not compatible")
+                else:
+                    tensor_index_type = iself.tensor_index_type
+                    tensor_metric = tensor_index_type.metric
+                    dummy = TensorIndex("d_" + str(count), tensor_index_type,
+                                        is_up=iself.is_up)
+                    if iself.is_up == iother.is_up:
+                        kroneckerdelta = tensor_index_type.delta(iself, -iother)
+                    else:
+                        kroneckerdelta = (
+                            TensMul(tensor_metric(iself, dummy),
+                                    tensor_index_type.delta(-dummy, -iother))
+                        )
+                    kronecker_delta_list.append(kroneckerdelta)
+            return TensMul.fromiter(kronecker_delta_list).doit()
+            # doit necessary to rename dummy indices accordingly
+
 
 class TensMul(TensExpr, AssocOp):
     """
@@ -3770,6 +3834,24 @@ def _eval_rewrite_as_Indexed(self, *args):
         expr = Mul.fromiter(args)
         return self._check_add_Sum(expr, index_symbols)
 
+    def _eval_partial_derivative(self, s):
+        # Evaluation like Mul
+        terms = []
+        for i, arg in enumerate(self.args):
+            # checking whether some tensor instance is differentiated
+            # or some other thing is necessary, but ugly
+            if isinstance(arg, TensExpr):
+                d = arg._eval_partial_derivative(s)
+            else:
+                # do not call diff is s is no symbol
+                if s._diff_wrt:
+                    d = arg._eval_derivative(s)
+                else:
+                    d = S.Zero
+            if d:
+                terms.append(TensMul.fromiter(self.args[:i] + (d,) + self.args[i + 1:]))
+        return TensAdd.fromiter(terms)
+
 
 class TensorElement(TensExpr):
     """
@@ -3870,6 +3952,7 @@ def canon_bp(p):
         return p.canon_bp()
     return p
 
+
 def tensor_mul(*a):
     """
     product of tensors
diff --git a/sympy/tensor/toperators.py b/sympy/tensor/toperators.py
index abe30585ad..40702bf1fc 100644
--- a/sympy/tensor/toperators.py
+++ b/sympy/tensor/toperators.py
@@ -1,5 +1,7 @@
+from sympy import Symbol, Number, sympify
 from sympy import MutableDenseNDimArray, S
-from sympy.tensor.tensor import (TensExpr, TensMul, TensorIndex)
+from sympy.tensor.tensor import (Tensor, TensExpr, TensAdd, TensMul,
+                                 TensorIndex)
 
 
 class PartialDerivative(TensExpr):
@@ -41,8 +43,6 @@ def __new__(cls, expr, *variables):
             variables = expr.variables + variables
             expr = expr.expr
 
-        # TODO: check that all variables have rank 1.
-
         args, indices, free, dum = cls._contract_indices_for_derivative(
             expr, variables)
 
@@ -64,16 +64,23 @@ def nocoeff(self):
     @classmethod
     def _contract_indices_for_derivative(cls, expr, variables):
         variables_opposite_valence = []
+
         for i in variables:
-            i_free_indices = i.get_free_indices()
-            variables_opposite_valence.append(i.xreplace({k: -k for k in i_free_indices}))
+            if isinstance(i, Tensor):
+                i_free_indices = i.get_free_indices()
+                variables_opposite_valence.append(
+                        i.xreplace({k: -k for k in i_free_indices}))
+            elif isinstance(i, Symbol):
+                variables_opposite_valence.append(i)
 
         args, indices, free, dum = TensMul._tensMul_contract_indices(
             [expr] + variables_opposite_valence, replace_indices=True)
 
         for i in range(1, len(args)):
-            i_indices = args[i].get_free_indices()
-            args[i] = args[i].xreplace({k: -k for k in i_indices})
+            args_i = args[i]
+            if isinstance(args_i, Tensor):
+                i_indices = args[i].get_free_indices()
+                args[i] = args[i].xreplace({k: -k for k in i_indices})
 
         return args, indices, free, dum
 
@@ -84,8 +91,63 @@ def doit(self):
         obj._indices = indices
         obj._free = free
         obj._dum = dum
+
         return obj
 
+    def _expand_partial_derivative(self):
+        args, indices, free, dum = self._contract_indices_for_derivative(self.expr, self.variables)
+
+        obj = self.func(*args)
+        obj._indices = indices
+        obj._free = free
+        obj._dum = dum
+
+        result = obj
+
+        if isinstance(obj.expr, TensAdd):
+            # take care of sums of multi PDs
+            result = obj.expr.func(*[
+                    self.func(a, *obj.variables)._expand_partial_derivative()
+                    for a in result.expr.args])
+        elif isinstance(obj.expr, TensMul):
+            # take care of products of multi PDs
+            if len(obj.variables) == 1:
+                # derivative with respect to single variable
+                terms = []
+                mulargs = list(obj.expr.args)
+                for ind in range(len(mulargs)):
+                    if not isinstance(sympify(mulargs[ind]), Number):
+                        # a number coefficient is not considered for
+                        # expansion of PartialDerivative
+                        d = self.func(mulargs[ind], *obj.variables)._expand_partial_derivative()
+                        terms.append(TensMul(*(mulargs[:ind]
+                                               + [d]
+                                               + mulargs[(ind + 1):])))
+                result = TensAdd.fromiter(terms)
+            else:
+                # derivative with respect to multiple variables
+                # decompose:
+                # partial(expr, (u, v))
+                # = partial(partial(expr, u).doit(), v).doit()
+                result = obj.expr  # init with expr
+                for v in obj.variables:
+                    result = self.func(result, v)._expand_partial_derivative()
+                    # then throw PD on it
+
+        return result
+
+    def _perform_derivative(self):
+        result = self.expr
+        for v in self.variables:
+            if isinstance(result, TensExpr):
+                result = result._eval_partial_derivative(v)
+            else:
+                if v._diff_wrt:
+                    result = result._eval_derivative(v)
+                else:
+                    result = S.Zero
+        return result
+
     def get_indices(self):
         return self._indices
 


Test patch:
diff --git a/sympy/printing/pretty/tests/test_pretty.py b/sympy/printing/pretty/tests/test_pretty.py
index 5b04588b24..802f160f9f 100644
--- a/sympy/printing/pretty/tests/test_pretty.py
+++ b/sympy/printing/pretty/tests/test_pretty.py
@@ -6543,14 +6543,14 @@ def test_pretty_print_tensor_expr():
     expr = A(i) + 3*B(i)
     ascii_str = \
 """\
- i      i\n\
-A  + 3*B \n\
+   i    i\n\
+3*B  + A \n\
          \
 """
     ucode_str = \
 u("""\
- i      i\n\
-A  + 3⋅B \n\
+   i    i\n\
+3⋅B  + A \n\
          \
 """)
     assert pretty(expr) == ascii_str
@@ -6611,16 +6611,16 @@ def test_pretty_print_tensor_partial_deriv():
     expr = A(i)*PartialDerivative(B(k)*C(-i) + 3*H(k, -i), A(j))
     ascii_str = \
 """\
- L_0  d / k           k   \\\n\
-A   *---|B *C    + 3*H    |\n\
+ L_0  d /   k       k     \\\n\
+A   *---|3*H     + B *C   |\n\
        j\\    L_0       L_0/\n\
      dA                    \n\
                            \
 """
     ucode_str = \
 u("""\
- L₀  ∂ ⎛ k          k  ⎞\n\
-A  ⋅───⎜B ⋅C   + 3⋅H   ⎟\n\
+ L₀  ∂ ⎛   k      k    ⎞\n\
+A  ⋅───⎜3⋅H    + B ⋅C  ⎟\n\
       j⎝    L₀       L₀⎠\n\
     ∂A                  \n\
                         \
diff --git a/sympy/tensor/tests/test_tensor.py b/sympy/tensor/tests/test_tensor.py
index c5b9bf9159..ab1e0377d8 100644
--- a/sympy/tensor/tests/test_tensor.py
+++ b/sympy/tensor/tests/test_tensor.py
@@ -516,6 +516,7 @@ def test_TensorHead():
     assert A.symmetry == TensorSymmetry.no_symmetry(2)
     assert A.comm == 0
 
+
 def test_add1():
     assert TensAdd().args == ()
     assert TensAdd().doit() == 0
@@ -536,7 +537,7 @@ def test_add1():
     t2b = t2 + t1
     assert str(t2b) == 'A(a, L_0)*A(b, -L_0) + A(b, -L_0)*B(L_0, a) + A(b, L_0)*B(a, -L_0)'
     t2b = t2b.canon_bp()
-    assert str(t2b) == '2*A(b, L_0)*B(a, -L_0) + A(a, L_0)*A(b, -L_0)'
+    assert str(t2b) == 'A(a, L_0)*A(b, -L_0) + 2*A(b, L_0)*B(a, -L_0)'
     p, q, r = tensor_heads('p,q,r', [Lorentz])
     t = q(d0)*2
     assert str(t) == '2*q(d0)'
@@ -550,16 +551,16 @@ def test_add1():
     t3 = t1*t2
     assert str(t3) == 'p(L_0)*(2*q(-L_0) + p(-L_0))'
     t3 = t3.expand()
-    assert str(t3) == '2*p(L_0)*q(-L_0) + p(L_0)*p(-L_0)'
+    assert str(t3) == 'p(L_0)*p(-L_0) + 2*p(L_0)*q(-L_0)'
     t3 = t2*t1
     t3 = t3.expand()
-    assert str(t3) == '2*q(-L_0)*p(L_0) + p(-L_0)*p(L_0)'
+    assert str(t3) == 'p(-L_0)*p(L_0) + 2*q(-L_0)*p(L_0)'
     t3 = t3.canon_bp()
-    assert str(t3) == '2*p(L_0)*q(-L_0) + p(L_0)*p(-L_0)'
+    assert str(t3) == 'p(L_0)*p(-L_0) + 2*p(L_0)*q(-L_0)'
     t1 = p(d0) + 2*q(d0)
     t3 = t1*t2
     t3 = t3.canon_bp()
-    assert str(t3) == '4*p(L_0)*q(-L_0) + 4*q(L_0)*q(-L_0) + p(L_0)*p(-L_0)'
+    assert str(t3) == 'p(L_0)*p(-L_0) + 4*p(L_0)*q(-L_0) + 4*q(L_0)*q(-L_0)'
     t1 = p(d0) - 2*q(d0)
     assert str(t1) == '-2*q(d0) + p(d0)'
     t2 = p(-d0) + 2*q(-d0)
@@ -652,16 +653,16 @@ def test_add3():
     B = TensorHead('B', [Lorentz])
 
     expr1 = A(i0)*A(-i0) - (E**2 - px**2 - py**2 - pz**2)
-    assert expr1.args == (px**2, py**2, pz**2, -E**2, A(i0)*A(-i0))
+    assert expr1.args == (-E**2, px**2, py**2, pz**2, A(i0)*A(-i0))
 
     expr2 = E**2 - px**2 - py**2 - pz**2 - A(i0)*A(-i0)
     assert expr2.args == (E**2, -px**2, -py**2, -pz**2, -A(i0)*A(-i0))
 
     expr3 = A(i0)*A(-i0) - E**2 + px**2 + py**2 + pz**2
-    assert expr3.args == (px**2, py**2, pz**2, -E**2, A(i0)*A(-i0))
+    assert expr3.args == (-E**2, px**2, py**2, pz**2, A(i0)*A(-i0))
 
     expr4 = B(i1)*B(-i1) + 2*E**2 - 2*px**2 - 2*py**2 - 2*pz**2 - A(i0)*A(-i0)
-    assert expr4.args == (-2*px**2, -2*py**2, -2*pz**2, 2*E**2, -A(i0)*A(-i0), B(i1)*B(-i1))
+    assert expr4.args == (2*E**2, -2*px**2, -2*py**2, -2*pz**2, B(i1)*B(-i1), -A(i0)*A(-i0))
 
 
 def test_mul():
@@ -1732,9 +1733,9 @@ def test_index_iteration():
     assert list(e5._iterate_dummy_indices) == [(L0, (0, 1, 0)), (L1, (0, 1, 1)), (-L0, (1, 1, 0)), (-L1, (1, 1, 1))]
     assert list(e5._iterate_indices) == [(L0, (0, 1, 0)), (L1, (0, 1, 1)), (-L0, (1, 1, 0)), (-L1, (1, 1, 1))]
 
-    assert list(e6._iterate_free_indices) == [(i0, (0, 1, 0)), (i2, (0, 1, 1)), (i0, (1, 0, 1, 0)), (i2, (1, 1, 1, 0))]
-    assert list(e6._iterate_dummy_indices) == [(L0, (1, 0, 1, 1)), (-L0, (1, 1, 1, 1))]
-    assert list(e6._iterate_indices) == [(i0, (0, 1, 0)), (i2, (0, 1, 1)), (i0, (1, 0, 1, 0)), (L0, (1, 0, 1, 1)), (i2, (1, 1, 1, 0)), (-L0, (1, 1, 1, 1))]
+    assert list(e6._iterate_free_indices) == [(i0, (0, 0, 1, 0)), (i2, (0, 1, 1, 0)), (i0, (1, 1, 0)), (i2, (1, 1, 1))]
+    assert list(e6._iterate_dummy_indices) == [(L0, (0, 0, 1, 1)), (-L0, (0, 1, 1, 1))]
+    assert list(e6._iterate_indices) == [(i0, (0, 0, 1, 0)), (L0, (0, 0, 1, 1)), (i2, (0, 1, 1, 0)), (-L0, (0, 1, 1, 1)), (i0, (1, 1, 0)), (i2, (1, 1, 1))]
 
     assert e1.get_indices() == [i0, i2]
     assert e1.get_free_indices() == [i0, i2]
@@ -1772,7 +1773,7 @@ def test_tensor_expand():
     assert expr != A(-i)*A(i)*A(j) + A(-i)*A(i)*B(j)*C(k)*C(-k)
     assert expr.expand() == A(-i)*A(i)*A(j) + A(-i)*A(i)*B(j)*C(k)*C(-k)
     assert str(expr) == "A(-L_0)*(A(L_0)*A(j) + A(L_0)*B(j)*C(L_1)*C(-L_1))"
-    assert str(expr.canon_bp()) == 'A(L_0)*A(-L_0)*B(j)*C(L_1)*C(-L_1) + A(j)*A(L_0)*A(-L_0)'
+    assert str(expr.canon_bp()) == 'A(j)*A(L_0)*A(-L_0) + A(L_0)*A(-L_0)*B(j)*C(L_1)*C(-L_1)'
 
     expr = A(-i)*(2*A(i)*A(j) + A(i)*B(j))
     assert expr.expand() == 2*A(-i)*A(i)*A(j) + A(-i)*A(i)*B(j)
diff --git a/sympy/tensor/tests/test_tensor_operators.py b/sympy/tensor/tests/test_tensor_operators.py
index 813fab94be..2fa11017c4 100644
--- a/sympy/tensor/tests/test_tensor_operators.py
+++ b/sympy/tensor/tests/test_tensor_operators.py
@@ -1,14 +1,18 @@
 from sympy.testing.pytest import raises
 
 from sympy.tensor.toperators import PartialDerivative
-from sympy.tensor.tensor import TensorIndexType, tensor_indices, TensorHead, \
-    tensor_heads
+from sympy.tensor.tensor import (TensorIndexType,
+                                 tensor_indices,
+                                 TensorHead, tensor_heads)
 from sympy import symbols, diag
-from sympy import Array
+from sympy import Array, Rational
+
+from sympy import sympify
+from random import randint
 
 
 L = TensorIndexType("L")
-i, j, k = tensor_indices("i j k", L)
+i, j, k, m, m1, m2, m3, m4 = tensor_indices("i j k m m1 m2 m3 m4", L)
 i0 = tensor_indices("i0", L)
 L_0, L_1 = tensor_indices("L_0 L_1", L)
 
@@ -117,3 +121,328 @@ def test_replace_arrays_partial_derivative():
     expr = A(i)*PartialDerivative(C(-j), D(i))
     assert expr.get_indices() == [L_0, -j, -L_0]
     assert expr.get_free_indices() == [-j]
+
+
+def test_expand_partial_derivative_sum_rule():
+    tau = symbols("tau")
+
+    # check sum rule for D(tensor, symbol)
+    expr1aa = PartialDerivative(A(i), tau)
+
+    assert expr1aa._expand_partial_derivative() == PartialDerivative(A(i), tau)
+
+    expr1ab = PartialDerivative(A(i) + B(i), tau)
+
+    assert (expr1ab._expand_partial_derivative() ==
+            PartialDerivative(A(i), tau) +
+            PartialDerivative(B(i), tau))
+
+    expr1ac = PartialDerivative(A(i) + B(i) + C(i), tau)
+
+    assert (expr1ac._expand_partial_derivative() ==
+            PartialDerivative(A(i), tau) +
+            PartialDerivative(B(i), tau) +
+            PartialDerivative(C(i), tau))
+
+    # check sum rule for D(tensor, D(j))
+    expr1ba = PartialDerivative(A(i), D(j))
+
+    assert expr1ba._expand_partial_derivative() ==\
+        PartialDerivative(A(i), D(j))
+    expr1bb = PartialDerivative(A(i) + B(i), D(j))
+
+    assert (expr1bb._expand_partial_derivative() ==
+            PartialDerivative(A(i), D(j)) +
+            PartialDerivative(B(i), D(j)))
+
+    expr1bc = PartialDerivative(A(i) + B(i) + C(i), D(j))
+    assert expr1bc._expand_partial_derivative() ==\
+        PartialDerivative(A(i), D(j))\
+        + PartialDerivative(B(i), D(j))\
+        + PartialDerivative(C(i), D(j))
+
+    # check sum rule for D(tensor, H(j, k))
+    expr1ca = PartialDerivative(A(i), H(j, k))
+    assert expr1ca._expand_partial_derivative() ==\
+        PartialDerivative(A(i), H(j, k))
+    expr1cb = PartialDerivative(A(i) + B(i), H(j, k))
+    assert (expr1cb._expand_partial_derivative() ==
+            PartialDerivative(A(i), H(j, k))
+            + PartialDerivative(B(i), H(j, k)))
+    expr1cc = PartialDerivative(A(i) + B(i) + C(i), H(j, k))
+    assert (expr1cc._expand_partial_derivative() ==
+            PartialDerivative(A(i), H(j, k))
+            + PartialDerivative(B(i), H(j, k))
+            + PartialDerivative(C(i), H(j, k)))
+
+    # check sum rule for D(D(tensor, D(j)), H(k, m))
+    expr1da = PartialDerivative(A(i), (D(j), H(k, m)))
+    assert expr1da._expand_partial_derivative() ==\
+        PartialDerivative(A(i), (D(j), H(k, m)))
+    expr1db = PartialDerivative(A(i) + B(i), (D(j), H(k, m)))
+    assert expr1db._expand_partial_derivative() ==\
+        PartialDerivative(A(i), (D(j), H(k, m)))\
+        + PartialDerivative(B(i), (D(j), H(k, m)))
+    expr1dc = PartialDerivative(A(i) + B(i) + C(i), (D(j), H(k, m)))
+    assert expr1dc._expand_partial_derivative() ==\
+        PartialDerivative(A(i), (D(j), H(k, m)))\
+        + PartialDerivative(B(i), (D(j), H(k, m)))\
+        + PartialDerivative(C(i), (D(j), H(k, m)))
+
+
+def test_expand_partial_derivative_constant_factor_rule():
+    pos_random_int1 = sympify(randint(0, 1000))
+    pos_random_int2 = sympify(randint(0, 1000))
+    neg_random_int = sympify(randint(-1000, -1))
+
+    c1 = Rational(pos_random_int1, pos_random_int2)
+    c2 = Rational(neg_random_int, pos_random_int2)
+    c3 = Rational(pos_random_int1, neg_random_int)
+
+    expr2a = PartialDerivative(pos_random_int1*A(i), D(j))
+    assert expr2a._expand_partial_derivative() ==\
+        pos_random_int1*PartialDerivative(A(i), D(j))
+
+    expr2b = PartialDerivative(neg_random_int*A(i), D(j))
+    assert expr2b._expand_partial_derivative() ==\
+        neg_random_int*PartialDerivative(A(i), D(j))
+
+    expr2ca = PartialDerivative(c1*A(i), D(j))
+    assert expr2ca._expand_partial_derivative() ==\
+        c1*PartialDerivative(A(i), D(j))
+
+    expr2cb = PartialDerivative(c2*A(i), D(j))
+    assert expr2cb._expand_partial_derivative() ==\
+        c2*PartialDerivative(A(i), D(j))
+
+    expr2cc = PartialDerivative(c3*A(i), D(j))
+    assert expr2cc._expand_partial_derivative() ==\
+        c3*PartialDerivative(A(i), D(j))
+
+
+def test_expand_partial_derivative_full_linearity():
+    pos_random_int1 = sympify(randint(0, 1000))
+    pos_random_int2 = sympify(randint(0, 1000))
+    neg_random_int = sympify(randint(-1000, -1))
+
+    c1 = Rational(pos_random_int1, pos_random_int2)
+    c2 = Rational(neg_random_int, pos_random_int2)
+    c3 = Rational(pos_random_int1, neg_random_int)
+
+    # check full linearity
+
+    expr3a = PartialDerivative(pos_random_int1*A(i) + pos_random_int2*B(i), D(j))
+    assert expr3a._expand_partial_derivative() ==\
+        pos_random_int1*PartialDerivative(A(i), D(j))\
+        + pos_random_int2*PartialDerivative(B(i), D(j))
+
+    expr3b = PartialDerivative(pos_random_int1*A(i) + neg_random_int*B(i), D(j))
+    assert expr3b._expand_partial_derivative() ==\
+        pos_random_int1*PartialDerivative(A(i), D(j))\
+        + neg_random_int*PartialDerivative(B(i), D(j))
+
+    expr3c = PartialDerivative(neg_random_int*A(i) + pos_random_int2*B(i), D(j))
+    assert expr3c._expand_partial_derivative() ==\
+        neg_random_int*PartialDerivative(A(i), D(j))\
+        + pos_random_int2*PartialDerivative(B(i), D(j))
+
+    expr3d = PartialDerivative(c1*A(i) + c2*B(i), D(j))
+    assert expr3d._expand_partial_derivative() ==\
+        c1*PartialDerivative(A(i), D(j))\
+        + c2*PartialDerivative(B(i), D(j))
+
+    expr3e = PartialDerivative(c2*A(i) + c1*B(i), D(j))
+    assert expr3e._expand_partial_derivative() ==\
+        c2*PartialDerivative(A(i), D(j))\
+        + c1*PartialDerivative(B(i), D(j))
+
+    expr3f = PartialDerivative(c2*A(i) + c3*B(i), D(j))
+    assert expr3f._expand_partial_derivative() ==\
+        c2*PartialDerivative(A(i), D(j))\
+        + c3*PartialDerivative(B(i), D(j))
+
+    expr3g = PartialDerivative(c3*A(i) + c2*B(i), D(j))
+    assert expr3g._expand_partial_derivative() ==\
+        c3*PartialDerivative(A(i), D(j))\
+        + c2*PartialDerivative(B(i), D(j))
+
+    expr3h = PartialDerivative(c3*A(i) + c1*B(i), D(j))
+    assert expr3h._expand_partial_derivative() ==\
+        c3*PartialDerivative(A(i), D(j))\
+        + c1*PartialDerivative(B(i), D(j))
+
+    expr3i = PartialDerivative(c1*A(i) + c3*B(i), D(j))
+    assert expr3i._expand_partial_derivative() ==\
+        c1*PartialDerivative(A(i), D(j))\
+        + c3*PartialDerivative(B(i), D(j))
+
+
+def test_expand_partial_derivative_product_rule():
+    # check product rule
+    expr4a = PartialDerivative(A(i)*B(j), D(k))
+
+    assert expr4a._expand_partial_derivative() == \
+        PartialDerivative(A(i), D(k))*B(j)\
+        + A(i)*PartialDerivative(B(j), D(k))
+
+    expr4b = PartialDerivative(A(i)*B(j)*C(k), D(m))
+    assert expr4b._expand_partial_derivative() ==\
+        PartialDerivative(A(i), D(m))*B(j)*C(k)\
+        + A(i)*PartialDerivative(B(j), D(m))*C(k)\
+        + A(i)*B(j)*PartialDerivative(C(k), D(m))
+
+    expr4c = PartialDerivative(A(i)*B(j), C(k), D(m))
+    assert expr4c._expand_partial_derivative() ==\
+        PartialDerivative(A(i), C(k), D(m))*B(j) \
+        + PartialDerivative(A(i), C(k))*PartialDerivative(B(j), D(m))\
+        + PartialDerivative(A(i), D(m))*PartialDerivative(B(j), C(k))\
+        + A(i)*PartialDerivative(B(j), C(k), D(m))
+
+
+def test_eval_partial_derivative_expr_by_symbol():
+
+    tau, alpha = symbols("tau alpha")
+
+    expr1 = PartialDerivative(tau**alpha, tau)
+    assert expr1._perform_derivative() == alpha * 1 / tau * tau ** alpha
+
+    expr2 = PartialDerivative(2*tau + 3*tau**4, tau)
+    assert expr2._perform_derivative() == 2 + 12 * tau ** 3
+
+    expr3 = PartialDerivative(2*tau + 3*tau**4, alpha)
+    assert expr3._perform_derivative() == 0
+
+
+def test_eval_partial_derivative_single_tensors_by_scalar():
+
+    tau, mu = symbols("tau mu")
+
+    expr = PartialDerivative(tau**mu, tau)
+    assert expr._perform_derivative() == mu*tau**mu/tau
+
+    expr1a = PartialDerivative(A(i), tau)
+    assert expr1a._perform_derivative() == 0
+
+    expr1b = PartialDerivative(A(-i), tau)
+    assert expr1b._perform_derivative() == 0
+
+    expr2a = PartialDerivative(H(i, j), tau)
+    assert expr2a._perform_derivative() == 0
+
+    expr2b = PartialDerivative(H(i, -j), tau)
+    assert expr2b._perform_derivative() == 0
+
+    expr2c = PartialDerivative(H(-i, j), tau)
+    assert expr2c._perform_derivative() == 0
+
+    expr2d = PartialDerivative(H(-i, -j), tau)
+    assert expr2d._perform_derivative() == 0
+
+
+def test_eval_partial_derivative_single_1st_rank_tensors_by_tensor():
+
+    expr1 = PartialDerivative(A(i), A(j))
+    assert expr1._perform_derivative() - L.delta(i, -j) == 0
+
+    expr2 = PartialDerivative(A(i), A(-j))
+    assert expr2._perform_derivative() - L.metric(i, L_0) * L.delta(-L_0, j) == 0
+
+    expr3 = PartialDerivative(A(-i), A(-j))
+    assert expr3._perform_derivative() - L.delta(-i, j) == 0
+
+    expr4 = PartialDerivative(A(-i), A(j))
+    assert expr4._perform_derivative() - L.metric(-i, -L_0) * L.delta(L_0, -j) == 0
+
+    expr5 = PartialDerivative(A(i), B(j))
+    expr6 = PartialDerivative(A(i), C(j))
+    expr7 = PartialDerivative(A(i), D(j))
+    expr8 = PartialDerivative(A(i), H(j, k))
+    assert expr5._perform_derivative() == 0
+    assert expr6._perform_derivative() == 0
+    assert expr7._perform_derivative() == 0
+    assert expr8._perform_derivative() == 0
+
+    expr9 = PartialDerivative(A(i), A(i))
+    assert expr9._perform_derivative() - L.delta(L_0, -L_0) == 0
+
+    expr10 = PartialDerivative(A(-i), A(-i))
+    assert expr10._perform_derivative() - L.delta(-L_0, L_0) == 0
+
+
+def test_eval_partial_derivative_single_2nd_rank_tensors_by_tensor():
+
+    expr1 = PartialDerivative(H(i, j), H(m, m1))
+    assert expr1._perform_derivative() - L.delta(i, -m) * L.delta(j, -m1) == 0
+
+    expr2 = PartialDerivative(H(i, j), H(-m, m1))
+    assert expr2._perform_derivative() - L.metric(i, L_0) * L.delta(-L_0, m) * L.delta(j, -m1) == 0
+
+    expr3 = PartialDerivative(H(i, j), H(m, -m1))
+    assert expr3._perform_derivative() - L.delta(i, -m) * L.metric(j, L_0) * L.delta(-L_0, m1) == 0
+
+    expr4 = PartialDerivative(H(i, j), H(-m, -m1))
+    assert expr4._perform_derivative() - L.metric(i, L_0) * L.delta(-L_0, m) * L.metric(j, L_1) * L.delta(-L_1, m1) == 0
+
+def test_eval_partial_derivative_divergence_type():
+    expr1a = PartialDerivative(A(i), A(i))
+    expr1b = PartialDerivative(A(i), A(k))
+    expr1c = PartialDerivative(L.delta(-i, k) * A(i), A(k))
+
+    assert (expr1a._perform_derivative()
+            - (L.delta(-i, k) * expr1b._perform_derivative())).contract_delta(L.delta) == 0
+
+    assert (expr1a._perform_derivative()
+            - expr1c._perform_derivative()).contract_delta(L.delta) == 0
+
+    expr2a = PartialDerivative(H(i, j), H(i, j))
+    expr2b = PartialDerivative(H(i, j), H(k, m))
+    expr2c = PartialDerivative(L.delta(-i, k) * L.delta(-j, m) * H(i, j), H(k, m))
+
+    assert (expr2a._perform_derivative()
+            - (L.delta(-i, k) * L.delta(-j, m) * expr2b._perform_derivative())).contract_delta(L.delta) == 0
+
+    assert (expr2a._perform_derivative()
+            - expr2c._perform_derivative()).contract_delta(L.delta) == 0
+
+
+
+def test_eval_partial_derivative_expr1():
+
+    tau, alpha = symbols("tau alpha")
+
+    # this is only some special expression
+    # tested: vector derivative
+    # tested: scalar derivative
+    # tested: tensor derivative
+    base_expr1 = A(i)*H(-i, j) + A(i)*A(-i)*A(j) + tau**alpha*A(j)
+
+    tensor_derivative = PartialDerivative(base_expr1, H(k, m))._perform_derivative()
+    vector_derivative = PartialDerivative(base_expr1, A(k))._perform_derivative()
+    scalar_derivative = PartialDerivative(base_expr1, tau)._perform_derivative()
+
+    assert (tensor_derivative - A(L_0)*L.metric(-L_0, -L_1)*L.delta(L_1, -k)*L.delta(j, -m)) == 0
+
+    assert (vector_derivative - (tau**alpha*L.delta(j, -k) +
+        L.delta(L_0, -k)*A(-L_0)*A(j) +
+        A(L_0)*L.metric(-L_0, -L_1)*L.delta(L_1, -k)*A(j) +
+        A(L_0)*A(-L_0)*L.delta(j, -k) +
+        L.delta(L_0, -k)*H(-L_0, j))).expand() == 0
+
+    assert (vector_derivative.contract_metric(L.metric).contract_delta(L.delta) -
+        (tau**alpha*L.delta(j, -k) + A(L_0)*A(-L_0)*L.delta(j, -k) + H(-k, j) + 2*A(j)*A(-k))).expand() == 0
+
+    assert scalar_derivative - alpha*1/tau*tau**alpha*A(j) == 0
+
+
+def test_eval_partial_derivative_mixed_scalar_tensor_expr2():
+
+    tau, alpha = symbols("tau alpha")
+
+    base_expr2 = A(i)*A(-i) + tau**2
+
+    vector_expression = PartialDerivative(base_expr2, A(k))._perform_derivative()
+    assert  (vector_expression -
+        (L.delta(L_0, -k)*A(-L_0) + A(L_0)*L.metric(-L_0, -L_1)*L.delta(L_1, -k))).expand() == 0
+
+    scalar_expression = PartialDerivative(base_expr2, tau)._perform_derivative()
+    assert scalar_expression == 2*tau