mirror of
https://github.com/JuliaFEM/JuliaFEM.jl.git
synced 2026-08-06 04:21:33 +00:00
Jukka Aho
383 lines
8.7 KiB
Plaintext
383 lines
8.7 KiB
Plaintext
{
|
||
"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 <jukka.aho@kapsi.fi>\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
|
||
}
|