Automatic Chapman Enskog Analysis of moment-based methods

chapman_enskog/ChapmanEnskog.ipynb

+%% Cell type:code id: tags:
+
+``` python
+import sympy as sp
+
+# Using operator and kets from quantum module to represent differential operators and pdfs
+from  sympy.physics.quantum import Ket as Func
+from  sympy.physics.quantum import Operator
+# Disable Ket notation |f> for functions
+Func.lbracket_latex =''
+Func.rbracket_latex = ''
+
+sp.init_printing()
+```
+
+%% Cell type:markdown id: tags:
+
+# Chapman Enskog analysis
+
+Particle distribution function $f$:
+
+%% Cell type:code id: tags:
+
+``` python
+c = sp.Symbol("c_x")
+dt = sp.Symbol("Delta_t")
+t = sp.symbols("t")
+f = Func("f")
+C = Func("C")
+```
+
+%% Cell type:markdown id: tags:
+
+Differential operators (defined simply as non-commutative symbols here)
+
+%% Cell type:code id: tags:
+
+``` python
+Dx = Operator("\partial_x")
+Dt = Operator("\partial_t")
+```
+
+%% Cell type:code id: tags:
+
+``` python
+taylorOrder = 2
+taylorOperator = sum(dt**k * (Dt + c* Dx)**k / sp.functions.factorial(k)
+                     for k in range(1, taylorOrder+1))
+taylorOperator
+```
+
+%%%% Output: execute_result
+
+    ![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAASYAAAAvBAMAAAC8knqVAAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAZjJ2me8QRFSJqyLN3buI9fnTAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAEq0lEQVRYCc1WTYhbVRT+kiYvefmZGVGKoNCIuivMoIPSuvChVVBoGQQRQWgWunHEyaZQBTHgpi7U6BRhGMWgS4tm0Y3IOKkUBI0QUanVhaMbQRdO7TAiFcd7z7nn/mRehYbIvLt479zvnvPd773cvO8AMuJjD0qYmfsqzmRGiwj5AbM9ibNyv4D1Vla0KB0zrOW1ZnY05S+zlq+zIwmF4/SC4iRDmk4e6Gg1N2RIEpLitpJT6pQamVFVauf+UmIOnnq1mRlNBUAfqOM7O5mRhJMAH6hxJf0PrtQF6ECNKwmTd6VyG8htjS1IFU7elSpajj5Qn+pgnDF5VzqlZRxIEOkPwphj0q7U1TqKV1AejilIlU3Ylcrnj6hxdKt69HxjXFGTdqX6Do9mrW0kVU2bwNPcodU0qQLfpBfTXMknCRmvX+n9B2POX5yWySdBxSqebQUATwSuqKo0V5Jeg7IDxrgX/ZJCqL4ovNH73uIJib+QQN9zfeQXfIBjC0f9dFcyvQZnB4w3A7czHFwt4zMefDjiF5VLPBAFJeiSA8pzHDt42Xell2ym6TVoHjI+AEzP2DxIiWWMzRY65f6YE+s9V8DRNw4QTRph+Cm3CLsBTK/Ba7sY64krEk0aIcZowS0+fwfHaw7SUf5G/OMQq8nCtZZbdZpMr8FLI4y3tKc3XI3VZBk/c4sSzetg+c4Gz/PnFt/Y4lBfRZODg++A3cD0GlwYMh7c/9Fswgv6KiWO8bRblOgnFay1a32en2hgyrTrGhBNDlZn0w3ZIOw1AsbSI8Dv3pGREseo/gSj4zsFPI5Km/Cq0lP/06UYTQbWDhl8/mUDdE2vQR4aMK4rPT96v7cp8Rhfgflk8m1TbfIYULDWV1FRccFoigeDL58YDNR+DJND5nUNjbcGg28Hg891LL0Ge2jA+GETeJTy1cWVeIyzsuruiqG4IdPpPrDUkZn97RimV+Q0qSx5T9Jr8EsMGP8GIt1qyzAlHmOKJvWmSUU8f+u7WEqAe5sUEon57QheJIcsDYVd3UWT6TWMh/qM0RawbxvL9x1ucJ0p8Rhf8AhNqE6k/qa176r0L2BpDrgICmlZNBFMDpl6xrs6WfUa7KE+Y6QPaGdf80y9RYTyGLwR5aec8dPqPHXxZKNRmwGmuqhsgEKiMJoYJocsdJibruahpdcAe2jA+CvwAfI4J1WmxGN8SJbcXb+6levm1NvqAdWvFo+oufVn+d8RTA5Za7hSeWjba7CHBoxr+88mquKKVBlNvBHlm3737LGfJQfWCWbRUIf6daUMFOoMo4lhcsjgI2020Jk02ENDxpd1e5PfzJsUKaGNdH51SCtRF++0TA5yCUXVzbfRYcwLjU8zTg6pPNiNRRdSxB66m3E9HsYmNSjR+fEcrVSa7iME/EZY9N5zh5pc54UMmKt2yKhvJqk346G7GJ++e/nNtAKdb3qV4hCVP2xO0IFZ9CrBVO8qCz58TYz8AKht+5rCTtXnTomp903BA+haGL3et27tJGDb08mL5jzvqYiRzb8fmWdgWkgyIGJEwvzIPAPTcoKPMyAjkHAbcE8A7P2kenHl4eHeywgU1FTnmyFN/wL4tmTZeQ4GZwAAAABJRU5ErkJggg==)
+    $$\frac{\Delta_{t}^{2}}{2} \left(c_{x} \partial_x + \partial_t\right)^{2} + \Delta_{t} \left(c_{x} \partial_x + \partial_t\right)$$
+      2                           2
+    Δₜ ⋅(cₓ⋅\partialₓ + \partialₜ)
+    ─────────────────────────────── + Δₜ⋅(cₓ⋅\partialₓ + \partialₜ)
+                   2
+
+%% Cell type:markdown id: tags:
+
+As right-hand-side we use the abstract collision operator $C$, which corresponds to (4.5)
+
+%% Cell type:code id: tags:
+
+``` python
+eq_4_5 = (taylorOperator * f) - C
+eq_4_5
+```
+
+%%%% Output: execute_result
+
+    ![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAX8AAAA1BAMAAABb3vbqAAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAiUSZq1TvELvdZiIyds1Wk1T5AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAHE0lEQVRoBdVYb4hUVRQ/8/ft/NnZET+E2eYEJfRFJwxTMpwoKip1ooz8EDuV9AdKhwwxFtyt/NIfaDSzLNKBCCGxnfBLSx/mQdCHIHdtpYKKHbI+aDGtmhiZ2Tn33HvfvbM7My9hxvHC3HPuOb/fOfe+9+595w1AQwusuKXBcoUNF8LaK2zGDdN9HyYqDabeHwYy3hx3QdUYefbLoDk5v0nnM7DMYqjol9dp3BKfCQI5AUxeYPwHPmmdh8Vr/nI8zJc8tFnIQNYfqxuol3wlcR5n2A1jaVI+9EXqDmhr3k+ePnnJs9FzCA+mg64fVlcwkYKfNHLTBvPxUwg/sv6Roh9WdzA7/aR5i0EhANoEmy9e9EPqEmao0j5R/DxjbgDgTdCeMjeiE0VITD7ec2dka2yUZQ5AbAIeXULfiSIkTtuyTZvICEAkDxA/0wbb0t2JIsSR76ZWiZezM0yCNsHvrcCtfB0pQuqtMrJvI4v1JMay4Pi4acyY3XegCJluu4uTf/NEciSi5yFS4vGl9B0oQqrldhPp50MocnwZtpVnEiuPu+0ozfydKEKiuWbZlL2vILTYRW7FFG5m0RJlqQgRP7rQHGr90Dxxjw+SYa4ipKyRAHZEyTT8QpWJ4vLR6Wv7QERH7RADkgmfWfaFsC1jGXgQqDgvkxZG1lxFiKpwBdqKqJgcx+tVovvYFJiRri1TamIK63x/572JNIxklYHlYjX8TSkksbBNFkyD1D8B+JJUpzZ3ESIrXERgsyIqJrt0rxNtZ1P8tHTdPpTXIFa+LUNkeRkmJm37MYdXGrcWFioAyN1O8Igi4Rk8UCbLoFmEPEoW0WSFK3Q7omYyUFFCBZkowCnULXROP89A3W97B9XVFWhc2G0BhsR4HRoPIGsmsugFoM5v+60GENRsALjCZd+siGad4FFkIqcgWI58t0ZmzASk/1HEDp+XMZdGXrv+R9bHPRNpyY/gL8/iLeCr/MAo2VMZ6mXTs5EVLpsbIiomOzVFJ/qZHX+yCJRY6j50ltQvAJq9KabIP/iTSwLnf3L4GnktaKgXcOTApxNZslhHqJ6NrHAJAGBH1Ex2KoqXaD47+LLF79ot0rCN+qpY0TDAoopnNLXXcDCeT9XYttiFvgusUq8WEHwW4M1Jslifr2o2YFW4VkSPSWzQT91inQh3ObVNLAZEFtZFv0gZVhcNq6G+jfpjEM4LUwInH1PnAVrUAqoYZX8m5aKlJIDc6QXkZIUr6isrIjM9jqTIRITfwM5nWFRdlrrfrwzr2CTfYyxm0IYfyiFdFoVRixYYCYF6/cSmej2Hw3VFgCdhBPskcUS7tl7fU6//QrqqcLm+siIykxngUTiRwE+wcw+LISGcE3ux7XNxsBmTAmAnF0Ajq2G66KiyDNQARtJq5N2Bf/ENcEpcK28BiFJ3IEwMzMS3x4rITC+ipHAigbcXsMCDsrafRMTFk1SshI1mjzdcTDkw9fX94m03XRSqwMhHiM64/nOH963C13DJIKsFrCfbWFbWV2ZEZsLgqmMu8ySFXqvTw6Ieu5Ede1nsYOH1+AbAZwR/LTYxvaLy34Vru2BkEmAnCFWEUAu4gFsjDb+ibc5NnCMwVripPCm4iXVE+lCJpfuLa2MZcumbxokE3t7EzzHM66tp1K/C37TrGU0NT7FQDra4bqoM0JeD8CgIVWDUJsbrcg8A1bMhiqeavJyqwgWur6yIgpmEkw0UTiTwN7HraSEc4wRhewQLgweLqI9l2NDY0x2c9/EkXrUK1pHvDS/DsS701ALGDxzOQrCArpSLnWpyAbrC5frKiiiYvHbBkhROJPDi5ALApxRbUJ8nYkjdlh13ZEgOiZ40gMMr3mAF+xjOW7QJcHG3XE1DoZJRLQAeOIi7IIevk3EyqyZno4bA9ZUdkZh4diUlSFFEIsInSsKDZwS1wKgQc3VGMefkYFFGYWTplZhZAGm2Gaos+NgemkTAoOKRxAVZjeur2RGrgZKsvGwK4a1i7iF6jpu0gax2hIveYQ/wurA7a647WmSEoWoKKcE1Ljg1VJo2WV/NivjCN4O0DWc1wm9na+QsyXdfnIXRhlhNq9EShP/RI+vzQ1ubKH2VJg7T/L8i8mohQCcEHLrVjGProRk9Tp0zF5Aoa0d75WB7SOMnZWuG+qQMe9NrQrD/+4rN3u1NeF0ye9VAs4RJ8ZAp79K00npEUtHXpm00/bvNQS/ozT5XjLmtNvRQ1hj0hPpq+1kszXiYKU/tDc3Pn7vGp3UkC5/3xsTVLPp9HCrGMYSfyD8oam/IaM7HPJ5QmMTOeTeX1KA35FDFxzz0n+Ip/KDssQUc9zF/SPXc0aNnHXlKqy2UhFf/tEBdFtfWjK+0G3yhLgfoFX9J8dOxN1t81Oe81vjEdRu2xG/CUMYvsqs4524f6f4D09Dm4HjqNRkAAAAASUVORK5CYII=)
+    $$\left(\frac{\Delta_{t}^{2}}{2} \left(c_{x} \partial_x + \partial_t\right)^{2} + \Delta_{t} \left(c_{x} \partial_x + \partial_t\right)\right) {f} - {C}$$
+    ⎛  2                           2                                ⎞
+    ⎜Δₜ ⋅(cₓ⋅\partialₓ + \partialₜ)                                 ⎟
+    ⎜─────────────────────────────── + Δₜ⋅(cₓ⋅\partialₓ + \partialₜ)⎟⋅❘f⟩ - ❘C⟩
+    ⎝               2                                               ⎠
+
+%% Cell type:markdown id: tags:
+
+Following the same steps as in the book, getting rid of second derivative, and discarding $\Delta_t^3$ we get to (4.7)
+
+%% Cell type:code id: tags:
+
+``` python
+eq_4_7 = eq_4_5 - (dt/2)* (Dt + c*Dx) * eq_4_5
+eq_4_7 = eq_4_7.expand().subs(dt**3, 0)
+eq_4_7
+```
+
+%% Cell type:markdown id: tags:
+
+### Chapman Enskog Ansatz
+
+Open Question:
+ why is not everything expanded equally (derivatives start at 1, spatial terminates one earlier...)
+
+%% Cell type:code id: tags:
+
+``` python
+eps = sp.Symbol("epsilon")
+
+ceSymbolsF = [Func("f{}".format(i)) for i in range(3)]
+ceSymbolsDt = [Operator("\partial_t^{{ ({}) }}".format(i)) for i in range(1,3)]
+ceSymbolsDx = [Operator("\partial_x^{{ ({}) }}".format(i)) for i in range(1,2)]
+
+ceF =  sum(eps**k * s for k, s in enumerate(ceSymbolsF, start=0))
+ceDt = sum(eps**k * s for k, s in enumerate(ceSymbolsDt, start=1))
+ceDx = sum(eps**k * s for k, s in enumerate(ceSymbolsDx, start=1))
+
+ceSubstitutions = {
+    Dt : ceDt,
+    Dx : ceDx,
+    f: ceF
+}
+ceSubstitutions
+```
+
+%% Cell type:markdown id: tags:
+
+Inserting the SRT/BGK collision operator
+
+%% Cell type:code id: tags:
+
+``` python
+srtC = -dt / sp.Symbol("tau") * ( ceF - ceSymbolsF[0])
+srtC
+```
+
+%% Cell type:code id: tags:
+
+``` python
+eq_4_7_ce = eq_4_7.subs(ceSubstitutions).subs(C, srtC).expand().collect(eps)
+eq_4_7_ce
+```
+
+%% Cell type:code id: tags:
+
+``` python
+eq_4_9_a = (eq_4_7_ce.coeff(eps) / dt).expand()
+eq_4_9_a
+```
+
+%% Cell type:code id: tags:
+
+``` python
+eq_4_9_b = (eq_4_7_ce.coeff(eps**2) / dt).expand()
+eq_4_9_b
+```
+
+%% Cell type:markdown id: tags:
+
+Computing moments
+
+%% Cell type:code id: tags:
+
+``` python
+import operator
+from functools import reduce
+def prod(factors):
+    return reduce(operator.mul, factors, 1)
+
+def getMomentSymbol(idx, order):
+    momentOrderStr = ['alpha', 'beta', 'gamma', 'delta']
+    momentOrderStr = "_".join(momentOrderStr[:order-1] + ['x'])
+    return sp.Symbol("Pi^(%d)_%s" % (idx, momentOrderStr), commutative=False)
+
+def takeMoments(eqn, order):
+    markerTerms = sp.symbols("c_alpha c_beta c_gamma c_delta")
+    velocityTermSet = set(markerTerms)
+    velocityTermSet.add(c)
+    factor = prod(markerTerms[:order])
+    eqn = (eqn*factor).expand()
+    result = 0
+    for fac in eqn.args:
+        occuringFs = list(fac.atoms(Func))
+        assert len(occuringFs) < 2
+        if len(occuringFs) == 0:
+            result += fac
+            continue
+        fIdx = ceSymbolsF.index(occuringFs[0])
+
+        occuringCs = fac.atoms(sp.Symbol).intersection(velocityTermSet)
+        facWithoutCs = sp.cancel(fac / prod(occuringCs) / occuringFs[0])
+        assert not facWithoutCs.atoms(sp.Symbol).intersection(velocityTermSet), "Non-linear velocity term"
+
+        # TODO handle multiple non-marker terms
+        moment = getMomentSymbol(fIdx, len(occuringCs))
+        #momentOrderStr = ['alpha', 'beta', 'gamma', 'delta']
+        #momentOrderStr = "_".join(momentOrderStr[:len(occuringCs)])
+        #moment = sp.Symbol("Pi^(%d)_%s" % (fIdx, momentOrderStr), commutative=False)
+        result += facWithoutCs * moment
+
+    return result
+```
+
+%% Cell type:code id: tags:
+
+``` python
+rho = Func("rho")
+u_x = Func("u_x")
+momentReplacements = {
+    getMomentSymbol(0, 0): rho,
+    getMomentSymbol(0, 1): rho * u_x,
+
+    getMomentSymbol(1, 0): 0, # According to eq 4.4
+    getMomentSymbol(1, 1): 0,
+    getMomentSymbol(2, 0): 0,
+    getMomentSymbol(2, 1): 0,
+}
+momentReplacements
+```
+
+%% Cell type:markdown id: tags:
+
+### Moments
+
+Moments of $\epsilon$-sorted equations
+
+$O(\epsilon)$
+
+%% Cell type:code id: tags:
+
+``` python
+eq_4_10_a = takeMoments(eq_4_9_a, 0).subs(momentReplacements)
+eq_4_10_a
+```
+
+%% Cell type:code id: tags:
+
+``` python
+eq_4_10_b = takeMoments(eq_4_9_a, 1).subs(momentReplacements)
+eq_4_10_b
+```
+
+%% Cell type:code id: tags:
+
+``` python
+eq_4_10_c = takeMoments(eq_4_9_a, 2).subs(momentReplacements)
+eq_4_10_c
+```
+
+%% Cell type:markdown id: tags:
+
+$O(\epsilon^2)$
+
+%% Cell type:code id: tags:
+
+``` python
+eq_4_12_a = takeMoments(eq_4_9_b, 0).subs(momentReplacements)
+eq_4_12_a
+```
+
+%% Cell type:code id: tags:
+
+``` python
+eq_4_12_b = takeMoments(eq_4_9_b, 1).subs(momentReplacements)
+eq_4_12_b
+```
+
+%% Cell type:markdown id: tags:
+
+Recombination
+
+%% Cell type:code id: tags:
+
+``` python
+(eps * eq_4_10_a + eps**2 * eq_4_12_a).expand()
+```
+
+%% Cell type:code id: tags:
+
+``` python
+(eps * eq_4_10_b + eps**2 * eq_4_12_b).expand()
+```
+
+%% Cell type:code id: tags:
+
+``` python
+stressTensorSubs = -
+```
+
+%% Cell type:code id: tags:
+
+``` python
+from lbmpy.maxwellian_equilibrium import getMomentsOfDiscreteMaxwellianEquilibrium
+from lbmpy.stencils import getStencil
+from lbmpy.moments import momentsOfOrder
+stencil = getStencil("D3Q19")
+dim = len(stencil[0])
+moments = tuple(momentsOfOrder(2, dim=dim))
+print(moments)
+u = Func("u_x"), Func("u_y"), Func("u_z")
+rho = Func("rho")
+pi_ab = getMomentsOfDiscreteMaxwellianEquilibrium(stencil, moments,c_s_sq=sp.Rational(1,3), order=2, u=u, rho=rho)
+pi_ab
+```
+
+%% Cell type:code id: tags:
+
+``` python
+from lbmpy.maxwellian_equilibrium import getMomentsOfContinuousMaxwellianEquilibrium
+pi_ab2 = getMomentsOfContinuousMaxwellianEquilibrium(moments, dim, c_s_sq=sp.Rational(1,3), order=2)
+pi_ab2
+```
+
+%% Cell type:code id: tags:
+
+``` python
+sp.solve(eq_4_10_c, getMomentSymbol(1, 2))[0]
+```
+
+%% Cell type:code id: tags:
+
+``` python
+from lbmpy.chapman_enskog import *
+productRule((ceSymbolsDt[0] * pi_ab[1]).expand())
+```
+
+%% Cell type:code id: tags:
+
+``` python
+dtSubsDict = {
+    ceSymbolsDt[0] * rho: -ceSymbolsDx[0]*(rho*u[0]),
+    productRule((ceSymbolsDt[0] * (rho * u[1])).expand()) : sp.Symbol("bla")
+}
+dtSubsDict
+```
+
+%% Cell type:code id: tags:
+
+``` python
+-ceSymbolsDx[0]*(rho*u[0])
+```
+
+%% Cell type:code id: tags:
+
+``` python
+```

chapman_enskog/__init__.py

+from lbmpy.chapman_enskog.derivative import DiffOperator, Diff, expandUsingLinearity, expandUsingProductRule, \
+    normalizeDiffOrder, chapmanEnskogDerivativeExpansion, chapmanEnskogDerivativeRecombination
+from lbmpy.chapman_enskog.chapman_enskog import getExpandedName, LbMethodEqMoments, insertMoments, takeMoments, \
+    CeMoment, chainSolveAndSubstitute, timeDiffSelector, momentSelector

chapman_enskog/chapman_enskog_old.py

+import sympy as sp
+from sympy.physics.quantum import Ket as Func
+from sympy.physics.quantum import Operator
+from functools import reduce
+import operator
+
+# Disable Ket notation |f> for functions
+Func.lbracket_latex =''
+Func.rbracket_latex = ''
+
+
+def prod(seq):
+    return reduce(operator.mul, seq, 1)
+
+
+def allIn(a, b):
+    """Tests if all elements of a are contained in b"""
+    return all(element in b for element in a)
+
+
+def isDerivative(term):
+    return isinstance(term, Operator) and term.args[0].name.startswith("\\partial")
+
+
+def getDiffOperator(subscript, superscript=None):
+    symbolName = "\partial_" + subscript
+    if superscript:
+        symbolName += "^" + superscript
+    return Operator(symbolName)
+
+
+def normalizeProduct(prod):
+    """Takes a sympy Mul node and returns a list of factors (with Pow nodes resolved)"""
+    def handlePow(pow):
+        if pow.exp.is_integer and pow.exp.is_number and pow.exp > 0:
+            return [pow.base] * pow.exp
+        else:
+            return [pow]
+
+    if prod.func == sp.Symbol or prod.func == Func:
+        return [prod]
+
+    if prod.func == sp.Pow:
+        return handlePow(prod)
+
+    assert prod.func == sp.Mul
+
+    result = []
+    for a in prod.args:
+        if a.func == sp.Pow:
+            result += handlePow(a)
+        else:
+            result.append(a)
+    return result
+
+
+def splitDerivativeProd(expr):
+    """Splits a term a b \partial c d  into three parts: functions before operator (here a b), 
+       the operator itself, and terms the operator acts on (here c d). If the passed expressions does
+       not have above form, None is returned"""
+    productArgs = normalizeProduct(expr)
+    argumentsActedOn = []
+    remainingArguments = []
+    derivOp = None
+    for t in reversed(productArgs):
+        if isDerivative(t) and derivOp is None:
+            derivOp = t
+        elif derivOp is None:
+            argumentsActedOn.append(t)
+        else:
+            remainingArguments.append(t)
+    if derivOp is None:
+        return None
+    return list(reversed(remainingArguments)), derivOp, argumentsActedOn
+
+
+def getName(obj):
+    if type(obj) in [Operator, Func]:
+        return obj.args[0].name
+    else:
+        return obj.name
+
+
+def combineDerivatives(expr):
+    from collections import defaultdict
+
+    def searchAndCombine(termList):
+        if len(termList) <= 1:
+            return termList, False
+
+        newTermList = []
+        handledIdxSet = set()
+
+        changes = False
+        for i in range(len(termList)):
+            if i in handledIdxSet:
+                continue
+            for j in range(i + 1, len(termList)):
+                pre1, post1 = termList[i]
+                pre2, post2 = termList[j]
+                rest1 = [item for item in pre1 if item not in post2]
+                rest2 = [item for item in pre2 if item not in post1]
+                restEqual = (len(rest1) == len(rest2)) and (set(rest1) == set(rest2))
+                if allIn(post1, pre2) and allIn(post2, pre1) and restEqual:
+                    handledIdxSet.add(j)
+                    newTermList.append((rest1, post1 + post2))
+                    changes = True
+                    break
+                else:
+                    newTermList.append(termList[i])
+        lastIdx = len(termList) - 1
+        if lastIdx not in handledIdxSet:
+            newTermList.append(termList[lastIdx])
+
+        return newTermList, changes
+
+    expr = expr.expand()
+    if expr.func != sp.Add:
+        return expr
+
+    operatorDict = defaultdict(list)  # maps the operator to a list of tuples with terms before and after the operator
+    result = 0
+    for term in expr.args:
+        splitted = splitDerivativeProd(term)
+        if splitted:
+            pre, op, post = splitted
+            operatorDict[op].append((pre, post))
+        else:
+            result += term
+
+    for op, termList in operatorDict.items():
+        newTermList, changes = searchAndCombine(termList)
+        while changes:
+            newTermList, changes = searchAndCombine(newTermList)
+
+        for pre, post in newTermList:
+            result += prod(pre) * op * prod(post)
+
+    return result
+
+
+def expandDerivatives(expr):
+    """Fully expands all derivatives by applying product rule"""
+
+    def handleProduct(term):
+        splittedTerm = splitDerivativeProd(term)
+        if splittedTerm is None:
+            return term
+        remainingArguments, derivOp, argumentsActedOn = splittedTerm
+
+        result = 0
+        for i in range(len(argumentsActedOn)):
+            beforeOp = prod([argumentsActedOn[j] for j in range(len(argumentsActedOn)) if j != i])
+            result += beforeOp * derivOp * argumentsActedOn[i]
+        return prod(remainingArguments) * result
+
+    expr = expr.expand()
+    if expr.func != sp.Add:
+        return handleProduct(expr)
+    else:
+        return sum(handleProduct(t) for t in expr.args)
+
+
+def commuteTerms(expr):
+    def commute(term):
+        forward = {obj: sp.Symbol(getName(obj)) for obj in term.atoms(Func)}
+        backward = {sp.Symbol(getName(obj)): obj for obj in term.atoms(Func)}
+        return term.subs(forward).subs(backward)
+
+    def handleProduct(term):
+        splittedTerm = splitDerivativeProd(term)
+        if splittedTerm is None:
+            return term
+        remainingArguments, derivOp, argumentsActedOn = splittedTerm
+
+        return commute(prod(remainingArguments)) * derivOp * commute(prod(argumentsActedOn))
+
+    expr = expr.expand()
+    if expr.func != sp.Add:
+        return handleProduct(expr)
+    else:
+        return sum(handleProduct(t) for t in expr.args)
+
+
+#def splitDerivativeTerm(term):
+#    if term.func == sp.Mul:
+#        derivativeIdx = [i for i, arg in enumerate(term.args) if isDerivative(arg)]
+#        if len(derivativeIdx) > 0:
+#            idx = derivativeIdx[-1]
+#            diffOperator = term.args[idx]
+#            assert not diffOperator.is_commutative
+#            derivedTerms = [t for t in term.args[idx + 1:] if not t.is_commutative]
+#
+#            newDerivedTerms = []
+#            for t in derivedTerms:
+#                if t.func == sp.Pow and t.exp.is_integer and t.exp.is_number and 0 < t.exp < 10:
+#                    newDerivedTerms += [t.base] * t.exp
+#                else:
+#                    newDerivedTerms.append(t)
+#            derivedTerms = newDerivedTerms
+#
+#            notDerivedTerms = list(term.args[:idx]) + [t for t in term.args[idx + 1:] if t.is_commutative]
+#            return diffOperator, notDerivedTerms, derivedTerms
+#    return None, [], []
+#
+#
+#def productRule(expr):
+#    """
+#    Applies product rule to terms with differential operator(s)
+#    :param expr: sympy expression to apply product rule to
+#    :return:
+#    """
+#    def visit(term):
+#        diffOperator, notDerivedTerms, derivedTerms = splitDerivativeTerm(term)
+#        if len(derivedTerms) > 1:
+#            result = 0
+#            for i in range(len(derivedTerms)):
+#                beforeOperator = [t for j, t in enumerate(derivedTerms) if j != i]
+#                result += prod(notDerivedTerms + beforeOperator + [diffOperator, derivedTerms[i]])
+#            return result
+#
+#        return term if not term.args else term.func(*[visit(a) for a in term.args])
+#
+#    return visit(expr)
+#
+#
+#def substituteDerivative(expr, diffOperatorToSearch, termsToSearch, substitution):
+#    def visit(term):
+#        diffOperator, notDerivedTerms, derivedTerms = splitDerivativeTerm(term)
+#        if diffOperator == diffOperatorToSearch:
+#            termMatches = True
+#            for innerTerm in termsToSearch:
+#                if innerTerm not in notDerivedTerms and innerTerm not in derivedTerms:
+#                    termMatches = False
+#                    break
+#            if termMatches:
+#                pass
+#
+#            if len(derivedTerms) > 1:
+#                result = 0
+#                for i in range(len(derivedTerms)):
+#                    beforeOperator = [t for j, t in enumerate(derivedTerms) if j != i]
+#                    result += prod(notDerivedTerms + beforeOperator + [diffOperator, derivedTerms[i]])
+#                return result
+#
+#        return term if not term.args else term.func(*[visit(a) for a in term.args])
+#
+#    return visit(expr)
+#
+import sympy
+import sympy.core
+import sympy as sp
+