{ "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": [ { "data": { "text/plain": [ "2015-06-15" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "using Dates\n", "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: 11.336512664 seconds (2020893880 bytes allocated, 26.15% 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 in 6 iterations\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 in \", i, \" iterations\")\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: 8.447644445 seconds (1554800080 bytes allocated, 29.45% gc time)\n" ] } ], "source": [ "@time test_algo2()" ] } ], "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 }