diff --git a/notebooks/2015-06-15-performance-studies.ipynb b/notebooks/2015-06-15-performance-studies.ipynb new file mode 100644 index 0000000..36b033a --- /dev/null +++ b/notebooks/2015-06-15-performance-studies.ipynb @@ -0,0 +1,389 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Compare analytical and autodiffed stiffness matrix\n", + "\n", + "Here we compare how much autodiffed solution is slower than analytical.\n", + "\n", + "Author(s): Jukka Aho \n", + "\n", + "Last updated:" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "ename": "LoadError", + "evalue": "today not defined\nwhile loading In[1], in expression starting on line 1", + "output_type": "error", + "traceback": [ + "today not defined\nwhile loading In[1], in expression starting on line 1", + "" + ] + } + ], + "source": [ + "today()" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "160" + ] + }, + "execution_count": 2, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "using ForwardDiff\n", + "ENV[\"COLUMNS\"] = 160" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "(36.0,24.0)" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Partial derivatives of bilinear Lagrange polynomials\n", + "dNdξ(ξ) = [[-(1-ξ[2])/4.0 -(1-ξ[1])/4.0],\n", + " [ (1-ξ[2])/4.0 -(1+ξ[1])/4.0],\n", + " [ (1+ξ[2])/4.0 (1+ξ[1])/4.0],\n", + " [-(1+ξ[2])/4.0 (1-ξ[1])/4.0]] \n", + "\n", + "a = 1/sqrt(3)\n", + "ipoints = [[-a -a], [a -a], [a a], [-a a]]\n", + "iweights = [1 1 1 1]\n", + "\n", + "E = 90\n", + "ν = 0.25\n", + "μ = E/(2*(1+ν))\n", + "λ = E*ν/((1+ν)*(1-2*ν))\n", + "λ = 2*λ*μ/(λ + 2*μ)\n", + "μ, λ" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Version using automatic differentiation" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "calc_local_matrices! (generic function with 1 method)" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "function calc_local_matrices!(X, u, R, Kt; dim=2)\n", + " I = eye(dim)\n", + " \n", + " function calc_Wint!(u, Wint)\n", + " for m = 1:length(iweights)\n", + " w = iweights[m]\n", + " ξ = ipoints[m, :]\n", + " Jᵀ = X*dNdξ(ξ)\n", + " ∇N = inv(Jᵀ)*dNdξ(ξ)'\n", + " ∇u = u*∇N'\n", + " F = I + ∇u # Deformation gradient\n", + " E = 1/2*(∇u' + ∇u + ∇u'*∇u) # Green-Lagrange strain tensor\n", + " S = λ*trace(E)*I + 2*μ*E # PK2 stress tensor\n", + " P = F*S # PK1 stress tensor\n", + " Wint[:,:] += w*P*∇N*det(Jᵀ)\n", + " end\n", + " end\n", + "\n", + " # herlper for tangent stiffness matrix\n", + " function R!(u, R)\n", + " R[:] = 0\n", + " calc_Wint!(reshape(u, 2, 4), reshape(R, 2, 4))\n", + " #calc_Wext!(reshape(u, 2, 4), reshape(R, 2, 4))\n", + " end\n", + " Jacobian = ForwardDiff.forwarddiff_jacobian(R!, Float64, fadtype=:dual, n=8, m=8)\n", + "\n", + " Kt[:,:] = Jacobian(reshape(u, 8))\n", + " R!(reshape(u, 8), reshape(R, 8))\n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Converged in 6 iterations\n", + "[0.0 -0.3991450609547433 -0.07228582695592461 0.0\n", + " 0.0 -2.1779892317073504 -2.222244754401764 0.0]\n" + ] + } + ], + "source": [ + "# validation\n", + "X = [0 0; 10 0; 10 1; 0 1]'\n", + "u = zeros(2,4)\n", + "R = zeros(2,4)\n", + "Kt = zeros(8,8)\n", + "\n", + "free_dofs = [3, 4, 5, 6]\n", + "for i in 1:10\n", + " calc_local_matrices!(X, u, R, Kt)\n", + " R[2,3] += 2\n", + " du = Kt[free_dofs, free_dofs] \\ -reshape(R, 8)[free_dofs]\n", + " u[free_dofs] += du\n", + " if norm(du) < 1.0e-9\n", + " println(\"Converged in \", i, \" iterations\")\n", + " break\n", + " end\n", + "end\n", + "println(u)" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "function test_algo1(N=10000)\n", + " for i=1:N\n", + " calc_local_matrices!(X, u, R, Kt)\n", + " end\n", + "end\n", + "test_algo1()" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "elapsed time: 12.592200168 seconds (2020893880 bytes allocated, 24.21% gc time)\n" + ] + } + ], + "source": [ + "@time test_algo1()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Analytical tangent stiffness" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "calc_local_matrices2! (generic function with 1 method)" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "function calc_local_matrices2!(X, u, R, Kt; dim=2)\n", + " I = eye(dim)\n", + " R[:,:] = 0.0\n", + " Kt[:,:] = 0.0\n", + " N = 4 # number of shape functions\n", + "\n", + " dF = zeros(2, 2)\n", + "\n", + " for m = 1:length(iweights)\n", + " w = iweights[m]\n", + " ξ = ipoints[m, :]\n", + " Jᵀ = X*dNdξ(ξ)\n", + " detJ = det(Jᵀ)\n", + " ∇N = inv(Jᵀ)*dNdξ(ξ)'\n", + " ∇u = u*∇N'\n", + " F = I + ∇u # Deformation gradient\n", + " E = 1/2*(∇u' + ∇u + ∇u'*∇u) # Green-Lagrange strain tensor\n", + " S = λ*trace(E)*I + 2*μ*E # PK2 stress tensor\n", + " P = F*S # PK1 stress tensor\n", + " R[:,:] += w*P*∇N*detJ\n", + "\n", + " for p = 1:N\n", + " for i = 1:dim\n", + " dF[:,:] = 0.0\n", + " dF[i,:] = ∇N[:,p]\n", + " dE = 1/2*(F'*dF + dF'*F)\n", + " dS = λ*trace(dE)*I + 2*μ*dE\n", + " dP = dF*S + F*dS\n", + " for q = 1:N\n", + " for j = 1:dim\n", + " Kt[dim*(p-1)+i,dim*(q-1)+j] += w*(dP[j,:]*∇N[:,q])[1]*detJ\n", + " end\n", + " end\n", + " end\n", + " end\n", + "\n", + " end\n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Converged\n", + "[0.0 -0.39914506095474317 -0.07228582695592449 0.0\n", + " 0.0 -2.1779892317073504 -2.222244754401764 0.0]\n" + ] + } + ], + "source": [ + "# validation\n", + "X = [0 0; 10 0; 10 1; 0 1]'\n", + "u = zeros(2,4)\n", + "R = zeros(2,4)\n", + "Kt = zeros(8,8)\n", + "\n", + "free_dofs = [3, 4, 5, 6]\n", + "for i in 1:10\n", + " calc_local_matrices2!(X, u, R, Kt)\n", + " R[2,3] += 2\n", + " du = Kt[free_dofs, free_dofs] \\ -reshape(R, 8)[free_dofs]\n", + " u[free_dofs] += du\n", + " if norm(du) < 1.0e-9\n", + " println(\"Converged\")\n", + " break\n", + " end\n", + "end\n", + "println(u)" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "function test_algo2(N=10000)\n", + " for i=1:N\n", + " calc_local_matrices2!(X, u, R, Kt)\n", + " end\n", + "end\n", + "test_algo2()" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "elapsed time: 9.96484533 seconds (1554800080 bytes allocated, 27.40% gc time)\n" + ] + } + ], + "source": [ + "@time test_algo2()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Julia 0.3.8", + "language": "julia", + "name": "julia-0.3" + }, + "language_info": { + "name": "julia", + "version": "0.3.8" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +}