diff --git a/meta/Velocity gradient.ipynb b/meta/Velocity gradient.ipynb new file mode 100644 index 0000000000000000000000000000000000000000..deaa7e596a699113671997b70b4d5e3e3517b722 --- /dev/null +++ b/meta/Velocity gradient.ipynb @@ -0,0 +1,197 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 31, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "-Axx**2/2 - Axy*Ayx - Axz*Azx - Ayy**2/2 - Ayz*Azy - Azz**2/2\n", + "-Axx*(Axx**2/3 + Axy*Ayx + Axz*Azx) - Axy*(Ayx*Ayy + Ayz*Azx) - Axz*(Ayx*Azy + Azx*Azz) - Ayy*(Ayy**2/3 + Ayz*Azy) - Ayz*Azy*Azz - Azz**3/3\n", + "Axx**2 + Axy*(Axy/2 + Ayx) + Axz*(Axz/2 + Azx) + Ayx**2/2 + Ayy**2 + Ayz*(Ayz/2 + Azy) + Azx**2/2 + Azy**2/2 + Azz**2\n" + ] + } + ], + "source": [ + "import sympy as sp\n", + "\n", + "A = []\n", + "for deriv in ['x', 'y', 'z']:\n", + " A.append([])\n", + " for field in ['x', 'y', 'z']:\n", + " A[-1].append(sp.Symbol('A' + deriv + field))\n", + "\n", + "A = sp.Matrix(A)\n", + "\n", + "A2 = A**2\n", + "A3 = A**3\n", + "Q = -sp.horner(sp.simplify(sum(A2[i, i] for i in range(3))/2))\n", + "R = -sp.horner(sp.simplify(sum(A3[i, i] for i in range(3))/3))\n", + "print(Q)\n", + "print(R)\n", + "\n", + "S = (A + A.T)/2\n", + "S2 = S**2\n", + "trS2 = sp.horner(sp.simplify(sum(S2[i, i] for i in range(3))))\n", + "print(trS2)" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "-Axx**2/2 - Axy*Ayx - Axz*Azx - Ayy**2/2 - Ayz*Azy - Azz**2/2\n", + "3*DIV + 3*MUL + NEG + 3*POW + 5*SUB\n", + "-Axx**3/3 - Axx*Axy*Ayx - Axx*Axz*Azx - Axy*Ayx*Ayy - Axy*Ayz*Azx - Axz*Ayx*Azy - Axz*Azx*Azz - Ayy**3/3 - Ayy*Ayz*Azy - Ayz*Azy*Azz - Azz**3/3\n", + "3*DIV + 16*MUL + NEG + 3*POW + 10*SUB\n" + ] + } + ], + "source": [ + "def alt_measure(expr):\n", + " POW = sp.Symbol('POW')\n", + " count = sp.count_ops(expr, visual=True).subs(POW, 10)\n", + " count = count.replace(sp.Symbol, type(sp.S.One))\n", + " return count\n", + "\n", + "Qalt = sp.simplify(Q, measure = alt_measure)\n", + "print(Qalt)\n", + "print(sp.count_ops(Qalt, visual=True))\n", + "Ralt = sp.simplify(R, measure = alt_measure)\n", + "print(Ralt)\n", + "print(sp.count_ops(Ralt, visual=True))" + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "-Axx**2/2 - Axy*Ayx - Axz*Azx - Ayy**2/2 - Ayz*Azy - Azz**2/2\n", + "0\n" + ] + } + ], + "source": [ + "Qalt = - sum(A[i, (i+1)%3] * A[(i+1)%3, i] for i in range(3)) - sum(A[i, i]**2 for i in range(3))/2\n", + "print Qalt\n", + "print(sp.simplify(Qalt - Q))" + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "-Axx*(Axx**2/3 + Axy*Ayx + Axz*Azx) - Axy*Ayz*Azx - Axz*Ayx*Azy - Ayy*(Axy*Ayx + Ayy**2/3 + Ayz*Azy) - Azz*(Axz*Azx + Ayz*Azy + Azz**2/3)\n", + "6*ADD + 3*DIV + 13*MUL + NEG + 3*POW + 4*SUB\n", + "0\n" + ] + } + ], + "source": [ + "Ralt = - (sum(A[i, i]*(A[i, i]**2/3 + sum(A[i, (i+j)%3]*A[(i+j)%3, i]\n", + " for j in range(1, 3)))\n", + " for i in range(3)) +\n", + " A[0, 1]*A[1, 2]*A[2, 0] + A[0, 2]*A[1, 0]*A[2, 1])\n", + "print Ralt\n", + "print(sp.count_ops(Ralt, visual=True))\n", + "print(sp.simplify(Ralt - R))" + ] + }, + { + "cell_type": "code", + "execution_count": 34, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "0\n", + "0\n", + "0\n" + ] + } + ], + "source": [ + "AxxAxx = A[0, 0]**2\n", + "AyyAyy = A[1, 1]**2\n", + "AzzAzz = A[2, 2]**2\n", + "AxyAyx = A[0, 1]*A[1, 0]\n", + "AyzAzy = A[1, 2]*A[2, 1]\n", + "AzxAxz = A[2, 0]*A[0, 2]\n", + "\n", + "Qalt = - (AxxAxx + AyyAyy + AzzAzz)/2 - AxyAyx -AyzAzy - AzxAxz\n", + "print(sp.simplify(Qalt - Q))\n", + "Ralt = - (A[0, 0]*(AxxAxx/3 + AxyAyx + AzxAxz) +\n", + " A[1, 1]*(AyyAyy/3 + AxyAyx + AyzAzy) +\n", + " A[2, 2]*(AzzAzz/3 + AzxAxz + AyzAzy) +\n", + " A[0, 1]*A[1, 2]*A[2, 0] +\n", + " A[0, 2]*A[1, 0]*A[2, 1])\n", + "print sp.simplify(Ralt - R)\n", + "trS2alt = (AxxAxx + AyyAyy + AzzAzz +\n", + " ((A[0, 1] + A[1, 0])**2 +\n", + " (A[1, 2] + A[2, 1])**2 +\n", + " (A[2, 0] + A[0, 2])**2)/2)\n", + "print sp.simplify(trS2alt - trS2)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 2", + "language": "python", + "name": "python2" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 2 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython2", + "version": "2.7.9" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +}