mirror of
https://github.com/JuliaFEM/JuliaFEM.jl.git
synced 2026-09-20 01:59:59 +00:00
added performance study between FAD and analytical stiffness matrices
This commit is contained in:
Jukka Aho
@@ -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 <jukka.aho@kapsi.fi>\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
|
||||
}
|
||||
Reference in New Issue
Block a user