Merge branch 'master' of git://github.com/JuliaFEM/JuliaFEM.jl

This commit is contained in:
Olli Väinölä
2015-12-03 08:38:38 +02:00
6 changed files with 514 additions and 273 deletions
@@ -59,11 +59,33 @@
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": true
"collapsed": false
},
"outputs": [],
"source": [
"using JuliaFEM"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"addprocs(7)'\n",
"@everywhere using JuliaFEM"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"using JuliaFEM\n",
"using JuliaFEM.Core: Element, Tri3, Tet4, Tri6, Tet10\n",
"using JuliaFEM.Core: ElasticityProblem, DirichletProblem\n",
"using JuliaFEM.Core: DirectSolver"
@@ -71,7 +93,7 @@
},
{
"cell_type": "code",
"execution_count": 2,
"execution_count": 4,
"metadata": {
"collapsed": false
},
@@ -82,26 +104,25 @@
"text": [
"INFO: Registered handlers: Any[\"ELEMENT\",\"NODE\",\"NSET\",\"ELSET\"]\n",
"INFO: Parsing elements\n",
"INFO: 216652 elements found\n",
"INFO: Creating ELSET PISTON\n",
"INFO: Parsing elements\n",
"INFO: 24783 elements found\n"
"INFO: 62454 elements found\n",
"INFO: Creating ELSET PISTON\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
" "
" "
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO: Parsing elements\n",
"INFO: 3494 elements found\n",
"INFO: Creating element set BC1\n",
"INFO: Creating element set BC2\n",
"INFO: Creating element set BC3\n",
"INFO: model loaded.\n"
]
},
@@ -109,7 +130,7 @@
"name": "stdout",
"output_type": "stream",
"text": [
"5.760825 seconds (44.46 M allocations: 1.301 GB, 30.79% gc time)\n"
" 2.269207 seconds (14.52 M allocations: 433.635 MB, 13.38% gc time)\n"
]
}
],
@@ -126,7 +147,7 @@
" # Quadratic models\n",
" #model = open(JuliaFEM.Core.parse_abaqus, \"../geometry/piston/piston_19611_P2.inp\")\n",
" #model = open(JuliaFEM.Core.parse_abaqus, \"../geometry/piston/piston_55950_P2.inp\")\n",
" #model = open(JuliaFEM.Core.parse_abaqus, \"../geometry/piston/piston_107168_P2.inp\")\n",
" model = open(JuliaFEM.Core.parse_abaqus, \"/Temp/piston_107168_P2.inp\")\n",
" #model = open(JuliaFEM.Core.parse_abaqus, \"../geometry/piston/piston_345757_P2.inp\")\n",
"\n",
" info(\"model loaded.\")\n",
@@ -135,7 +156,7 @@
},
{
"cell_type": "code",
"execution_count": 3,
"execution_count": 5,
"metadata": {
"collapsed": false
},
@@ -144,7 +165,7 @@
"name": "stderr",
"output_type": "stream",
"text": [
"INFO: 241435 elements created.\n",
"INFO: 65948 elements created.\n",
"INFO: all ready for solver\n"
]
}
@@ -202,7 +223,7 @@
},
{
"cell_type": "code",
"execution_count": 4,
"execution_count": 6,
"metadata": {
"collapsed": true
},
@@ -214,7 +235,7 @@
},
{
"cell_type": "code",
"execution_count": 5,
"execution_count": 7,
"metadata": {
"collapsed": false
},
@@ -237,23 +258,385 @@
"INFO: Assemble: 80 % done\n",
"INFO: Assemble: 90 % done\n",
"INFO: Assemble: 100 % done\n",
"INFO: # of interface dofs: 19848\n",
"INFO: # of interface dofs: 8568\n",
"INFO: Assembling field problems...\n",
"INFO: Assembling body 1...\n",
"INFO: Assemble: 10 % done\n"
"INFO: Assemble: 10 % done\n",
"INFO: Assemble: 10 % done\n",
"INFO: Assemble: 10 % done\n",
"INFO: Assemble: 10 % done\n",
"INFO: Assemble: 10 % done\n",
"INFO: Assemble: 10 % done\n",
"INFO: Assemble: 10 % done\n",
"INFO: Assemble: 20 % done\n",
"INFO: Assemble: 20 % done\n",
"INFO: Assemble: 20 % done\n",
"INFO: Assemble: 20 % done\n",
"INFO: Assemble: 20 % done\n",
"INFO: Assemble: 20 % done\n",
"INFO: Assemble: 20 % done\n",
"INFO: Assemble: 30 % done\n",
"INFO: Assemble: 30 % done\n",
"INFO: Assemble: 30 % done\n",
"INFO: Assemble: 30 % done\n",
"INFO: Assemble: 30 % done\n",
"INFO: Assemble: 30 % done\n",
"INFO: Assemble: 30 % done\n",
"INFO: Assemble: 40 % done\n",
"INFO: Assemble: 40 % done\n",
"INFO: Assemble: 40 % done\n",
"INFO: Assemble: 40 % done\n",
"INFO: Assemble: 40 % done\n",
"INFO: Assemble: 40 % done\n",
"INFO: Assemble: 40 % done\n",
"INFO: Assemble: 50 % done\n",
"INFO: Assemble: 50 % done\n",
"INFO: Assemble: 50 % done\n",
"INFO: Assemble: 50 % done\n",
"INFO: Assemble: 50 % done\n",
"INFO: Assemble: 50 % done\n",
"INFO: Assemble: 50 % done\n",
"INFO: Assemble: 60 % done\n",
"INFO: Assemble: 60 % done\n",
"INFO: Assemble: 60 % done\n",
"INFO: Assemble: 60 % done\n",
"INFO: Assemble: 60 % done\n",
"INFO: Assemble: 60 % done\n",
"INFO: Assemble: 60 % done\n",
"INFO: Assemble: 70 % done\n",
"INFO: Assemble: 70 % done\n",
"INFO: Assemble: 70 % done\n",
"INFO: Assemble: 70 % done\n",
"INFO: Assemble: 70 % done\n",
"INFO: Assemble: 70 % done\n",
"INFO: Assemble: 70 % done\n",
"INFO: Assemble: 80 % done\n",
"INFO: Assemble: 80 % done\n",
"INFO: Assemble: 80 % done\n",
"INFO: Assemble: 80 % done\n",
"INFO: Assemble: 80 % done\n",
"INFO: Assemble: 80 % done\n",
"INFO: Assemble: 90 % done\n",
"INFO: Assemble: 80 % done\n",
"INFO: Assemble: 90 % done\n",
"INFO: Assemble: 90 % done\n",
"INFO: Assemble: 90 % done\n",
"INFO: Assemble: 90 % done\n",
"INFO: Assemble: 90 % done\n",
"INFO: Assemble: 90 % done\n",
"INFO: # of dofs in problem 1: 321504\n",
"INFO: Solving system\n",
"INFO: Solved, calculating interior dofs...\n",
"INFO: Problem solved. solution norm: 11.391177449683427\n",
"INFO: timing info for non-linear iteration:\n",
"INFO: boundary assembly : 3.686000108718872\n",
"INFO: field assembly : 162.32500004768372\n",
"INFO: reduce stiffness matrix : 5.069999933242798\n",
"INFO: create sparse matrices : 0.1399998664855957\n",
"INFO: dump matrices to disk : 0.0\n",
"INFO: solve problem : 29.73900008201599\n",
"INFO: back substitute : 0.0\n",
"INFO: update element data : 1.621999979019165\n",
"INFO: non-linear iteration : 204.1119999885559\n",
"INFO: Starting iteration 2\n",
"INFO: Assembling boundary problems...\n",
"INFO: Assemble: 10 % done\n",
"INFO: Assemble: 20 % done\n",
"INFO: Assemble: 30 % done\n",
"INFO: Assemble: 40 % done\n",
"INFO: Assemble: 50 % done\n",
"INFO: Assemble: 60 % done\n",
"INFO: Assemble: 70 % done\n",
"INFO: Assemble: 80 % done\n",
"INFO: Assemble: 90 % done\n",
"INFO: Assemble: 100 % done\n",
"INFO: # of interface dofs: 8568\n",
"INFO: Assembling field problems...\n",
"INFO: Assembling body 1...\n",
"INFO: Assemble: 10 % done\n",
"INFO: Assemble: 10 % done\n",
"INFO: Assemble: 10 % done\n",
"INFO: Assemble: 10 % done\n",
"INFO: Assemble: 10 % done\n",
"INFO: Assemble: 10 % done\n",
"INFO: Assemble: 10 % done\n",
"INFO: Assemble: 20 % done\n",
"INFO: Assemble: 20 % done\n",
"INFO: Assemble: 20 % done\n",
"INFO: Assemble: 20 % done\n",
"INFO: Assemble: 20 % done\n",
"INFO: Assemble: 20 % done\n",
"INFO: Assemble: 20 % done\n",
"INFO: Assemble: 30 % done\n",
"INFO: Assemble: 30 % done\n",
"INFO: Assemble: 30 % done\n",
"INFO: Assemble: 30 % done\n",
"INFO: Assemble: 30 % done\n",
"INFO: Assemble: 30 % done\n",
"INFO: Assemble: 30 % done\n",
"INFO: Assemble: 40 % done\n",
"INFO: Assemble: 40 % done\n",
"INFO: Assemble: 40 % done\n",
"INFO: Assemble: 40 % done\n",
"INFO: Assemble: 40 % done\n",
"INFO: Assemble: 40 % done\n",
"INFO: Assemble: 40 % done\n",
"INFO: Assemble: 50 % done\n",
"INFO: Assemble: 50 % done\n",
"INFO: Assemble: 50 % done\n",
"INFO: Assemble: 50 % done\n",
"INFO: Assemble: 50 % done\n",
"INFO: Assemble: 50 % done\n",
"INFO: Assemble: 50 % done\n",
"INFO: Assemble: 60 % done\n",
"INFO: Assemble: 60 % done\n",
"INFO: Assemble: 60 % done\n",
"INFO: Assemble: 60 % done\n",
"INFO: Assemble: 60 % done\n",
"INFO: Assemble: 60 % done\n",
"INFO: Assemble: 60 % done\n",
"INFO: Assemble: 70 % done\n",
"INFO: Assemble: 70 % done\n",
"INFO: Assemble: 70 % done\n",
"INFO: Assemble: 70 % done\n",
"INFO: Assemble: 70 % done\n",
"INFO: Assemble: 70 % done\n",
"INFO: Assemble: 70 % done\n",
"INFO: Assemble: 80 % done\n",
"INFO: Assemble: 80 % done\n",
"INFO: Assemble: 80 % done\n",
"INFO: Assemble: 90 % done\n",
"INFO: Assemble: 80 % done\n",
"INFO: Assemble: 80 % done\n",
"INFO: Assemble: 80 % done\n",
"INFO: Assemble: 80 % done\n",
"INFO: Assemble: 90 % done\n",
"INFO: Assemble: 90 % done\n",
"INFO: Assemble: 90 % done\n",
"INFO: Assemble: 90 % done\n",
"INFO: Assemble: 90 % done\n",
"INFO: Assemble: 90 % done\n",
"INFO: # of dofs in problem 1: 321504\n",
"INFO: Solving system\n",
"INFO: Solved, calculating interior dofs...\n",
"INFO: Problem solved. solution norm: 0.11391526987712108\n",
"INFO: timing info for non-linear iteration:\n",
"INFO: boundary assembly : 1.9660000801086426\n",
"INFO: field assembly : 160.2829999923706\n",
"INFO: reduce stiffness matrix : 5.057999849319458\n",
"INFO: create sparse matrices : 0.1420001983642578\n",
"INFO: dump matrices to disk : 0.0\n",
"INFO: solve problem : 60.674999952316284\n",
"INFO: back substitute : 0.0\n",
"INFO: update element data : 1.6159999370574951\n",
"INFO: non-linear iteration : 231.24799990653992\n",
"INFO: Starting iteration 3\n",
"INFO: Assembling boundary problems...\n",
"INFO: Assemble: 10 % done\n",
"INFO: Assemble: 20 % done\n",
"INFO: Assemble: 30 % done\n",
"INFO: Assemble: 40 % done\n",
"INFO: Assemble: 50 % done\n",
"INFO: Assemble: 60 % done\n",
"INFO: Assemble: 70 % done\n",
"INFO: Assemble: 80 % done\n",
"INFO: Assemble: 90 % done\n",
"INFO: Assemble: 100 % done\n",
"INFO: # of interface dofs: 8568\n",
"INFO: Assembling field problems...\n",
"INFO: Assembling body 1...\n",
"INFO: Assemble: 10 % done\n",
"INFO: Assemble: 10 % done\n",
"INFO: Assemble: 10 % done\n",
"INFO: Assemble: 10 % done\n",
"INFO: Assemble: 10 % done\n",
"INFO: Assemble: 10 % done\n",
"INFO: Assemble: 10 % done\n",
"INFO: Assemble: 20 % done\n",
"INFO: Assemble: 20 % done\n",
"INFO: Assemble: 20 % done\n",
"INFO: Assemble: 20 % done\n",
"INFO: Assemble: 20 % done\n",
"INFO: Assemble: 20 % done\n",
"INFO: Assemble: 20 % done\n",
"INFO: Assemble: 30 % done\n",
"INFO: Assemble: 30 % done\n",
"INFO: Assemble: 30 % done\n",
"INFO: Assemble: 30 % done\n",
"INFO: Assemble: 30 % done\n",
"INFO: Assemble: 30 % done\n",
"INFO: Assemble: 30 % done\n",
"INFO: Assemble: 40 % done\n",
"INFO: Assemble: 40 % done\n",
"INFO: Assemble: 40 % done\n",
"INFO: Assemble: 40 % done\n",
"INFO: Assemble: 40 % done\n",
"INFO: Assemble: 40 % done\n",
"INFO: Assemble: 40 % done\n",
"INFO: Assemble: 50 % done\n",
"INFO: Assemble: 50 % done\n",
"INFO: Assemble: 50 % done\n",
"INFO: Assemble: 50 % done\n",
"INFO: Assemble: 50 % done\n",
"INFO: Assemble: 50 % done\n",
"INFO: Assemble: 50 % done\n",
"INFO: Assemble: 60 % done\n",
"INFO: Assemble: 60 % done\n",
"INFO: Assemble: 60 % done\n",
"INFO: Assemble: 60 % done\n",
"INFO: Assemble: 60 % done\n",
"INFO: Assemble: 60 % done\n",
"INFO: Assemble: 60 % done\n",
"INFO: Assemble: 70 % done\n",
"INFO: Assemble: 70 % done\n",
"INFO: Assemble: 70 % done\n",
"INFO: Assemble: 70 % done\n",
"INFO: Assemble: 70 % done\n",
"INFO: Assemble: 70 % done\n",
"INFO: Assemble: 70 % done\n",
"INFO: Assemble: 80 % done\n",
"INFO: Assemble: 80 % done\n",
"INFO: Assemble: 80 % done\n",
"INFO: Assemble: 80 % done\n",
"INFO: Assemble: 90 % done\n",
"INFO: Assemble: 80 % done\n",
"INFO: Assemble: 80 % done\n",
"INFO: Assemble: 80 % done\n",
"INFO: Assemble: 90 % done\n",
"INFO: Assemble: 90 % done\n",
"INFO: Assemble: 90 % done\n",
"INFO: Assemble: 90 % done\n",
"INFO: Assemble: 90 % done\n",
"INFO: Assemble: 90 % done\n",
"INFO: # of dofs in problem 1: 321504\n",
"INFO: Solving system\n",
"INFO: Solved, calculating interior dofs...\n",
"INFO: Problem solved. solution norm: 0.0004840672223822576\n",
"INFO: timing info for non-linear iteration:\n",
"INFO: boundary assembly : 2.003999948501587\n",
"INFO: field assembly : 168.06199979782104\n",
"INFO: reduce stiffness matrix : 5.039999961853027\n",
"INFO: create sparse matrices : 0.14100003242492676\n",
"INFO: dump matrices to disk : 0.0\n",
"INFO: solve problem : 60.27800011634827\n",
"INFO: back substitute : 0.0\n",
"INFO: update element data : 1.61899995803833\n",
"INFO: non-linear iteration : 238.61800003051758\n",
"INFO: Starting iteration 4\n",
"INFO: Assembling boundary problems...\n",
"INFO: Assemble: 10 % done\n",
"INFO: Assemble: 20 % done\n",
"INFO: Assemble: 30 % done\n",
"INFO: Assemble: 40 % done\n",
"INFO: Assemble: 50 % done\n",
"INFO: Assemble: 60 % done\n",
"INFO: Assemble: 70 % done\n",
"INFO: Assemble: 80 % done\n",
"INFO: Assemble: 90 % done\n",
"INFO: Assemble: 100 % done\n",
"INFO: # of interface dofs: 8568\n",
"INFO: Assembling field problems...\n",
"INFO: Assembling body 1...\n",
"INFO: Assemble: 10 % done\n",
"INFO: Assemble: 10 % done\n",
"INFO: Assemble: 10 % done\n",
"INFO: Assemble: 10 % done\n",
"INFO: Assemble: 10 % done\n",
"INFO: Assemble: 10 % done\n",
"INFO: Assemble: 10 % done\n",
"INFO: Assemble: 20 % done\n",
"INFO: Assemble: 20 % done\n",
"INFO: Assemble: 20 % done\n",
"INFO: Assemble: 20 % done\n",
"INFO: Assemble: 20 % done\n",
"INFO: Assemble: 20 % done\n",
"INFO: Assemble: 20 % done\n",
"INFO: Assemble: 30 % done\n",
"INFO: Assemble: 30 % done\n",
"INFO: Assemble: 30 % done\n",
"INFO: Assemble: 30 % done\n",
"INFO: Assemble: 30 % done\n",
"INFO: Assemble: 30 % done\n",
"INFO: Assemble: 30 % done\n",
"INFO: Assemble: 40 % done\n",
"INFO: Assemble: 40 % done\n",
"INFO: Assemble: 40 % done\n",
"INFO: Assemble: 40 % done\n",
"INFO: Assemble: 40 % done\n",
"INFO: Assemble: 40 % done\n",
"INFO: Assemble: 40 % done\n",
"INFO: Assemble: 50 % done\n",
"INFO: Assemble: 50 % done\n",
"INFO: Assemble: 50 % done\n",
"INFO: Assemble: 50 % done\n",
"INFO: Assemble: 50 % done\n",
"INFO: Assemble: 50 % done\n",
"INFO: Assemble: 50 % done\n",
"INFO: Assemble: 60 % done\n",
"INFO: Assemble: 60 % done\n",
"INFO: Assemble: 60 % done\n",
"INFO: Assemble: 60 % done\n",
"INFO: Assemble: 60 % done\n",
"INFO: Assemble: 60 % done\n",
"INFO: Assemble: 60 % done\n",
"INFO: Assemble: 70 % done\n",
"INFO: Assemble: 70 % done\n",
"INFO: Assemble: 70 % done\n",
"INFO: Assemble: 70 % done\n",
"INFO: Assemble: 70 % done\n",
"INFO: Assemble: 70 % done\n",
"INFO: Assemble: 70 % done\n",
"INFO: Assemble: 80 % done\n",
"INFO: Assemble: 80 % done\n",
"INFO: Assemble: 90 % done\n",
"INFO: Assemble: 80 % done\n",
"INFO: Assemble: 80 % done\n",
"INFO: Assemble: 80 % done\n",
"INFO: Assemble: 80 % done\n",
"INFO: Assemble: 80 % done\n",
"INFO: Assemble: 90 % done\n",
"INFO: Assemble: 90 % done\n",
"INFO: Assemble: 90 % done\n",
"INFO: Assemble: 90 % done\n",
"INFO: Assemble: 90 % done\n",
"INFO: Assemble: 90 % done\n",
"INFO: # of dofs in problem 1: 321504\n",
"INFO: Solving system\n",
"INFO: Solved, calculating interior dofs...\n",
"INFO: Problem solved. solution norm: 5.318293313386342e-10\n"
]
},
{
"ename": "LoadError",
"evalue": "LoadError: InterruptException:\nwhile loading In[5], in expression starting on line 1",
"output_type": "error",
"traceback": [
"LoadError: InterruptException:\nwhile loading In[5], in expression starting on line 1",
""
"data": {
"text/plain": [
"(4,true)"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO: timing info for non-linear iteration:\n",
"INFO: boundary assembly : 2.0319998264312744\n",
"INFO: field assembly : 166.31699991226196\n",
"INFO: reduce stiffness matrix : 5.141999959945679\n",
"INFO: create sparse matrices : 0.1419999599456787\n",
"INFO: dump matrices to disk : 0.0\n",
"INFO: solve problem : 60.019999980926514\n",
"INFO: back substitute : 0.0\n",
"INFO: update element data : 1.621999979019165\n",
"INFO: non-linear iteration : 236.77899980545044\n",
"INFO: solver finished in 912.1029999256134 seconds.\n"
]
}
],
"source": [
"solver.reduce_stiffness_matrix = false\n",
"#solver.dump_matrices = true\n",
"call(solver, 0.0)"
]
},
+22 -7
View File
@@ -13,6 +13,12 @@ type CAssembly
fi :: SparseMatrixCSC
end
function optimize!(assembly::Assembly)
optimize!(assembly.mass_matrix)
optimize!(assembly.stiffness_matrix)
optimize!(assembly.force_vector)
end
function assemble!(assembly::Assembly, problem::AllProblems, time::Float64, empty_assembly::Bool=true)
if empty_assembly
empty!(assembly)
@@ -22,19 +28,28 @@ function assemble!(assembly::Assembly, problem::AllProblems, time::Float64, empt
end
end
function assemble(problem::AllProblems, time::Float64)
function assemble(problem::AllProblems, elrange::UnitRange{Int64}, time::Real)
elements = get_elements(problem)[elrange]
assembly = Assembly()
ne = length(get_elements(problem))
ne = length(elrange)
p = ne > 10 ? round(Int, ne/10) : ne
for (i, element) in enumerate(get_elements(problem))
for (i, element) in enumerate(elements)
mod(i, p) == 0 && info("Assemble: ", round(Int, i/ne*100), " % done")
assemble!(assembly, problem, element, time)
end
optimize!(assembly)
return assembly
end
""" Calculate reduced stiffness matrix. """
function reduce(assembly::Assembly, boundary_dofs_::Vector{Int})
function assemble(problem::AllProblems, time::Real)
ne = length(get_elements(problem))
assemble(problem, 1:ne, time)
end
""" Calculate reduced stiffness matrix.
mindofs: if dofs < mindofs, do not reduce
"""
function reduce(assembly::Assembly, boundary_dofs_::Vector{Int}, mindofs=100000)
all_dofs = unique(assembly.stiffness_matrix.I)
boundary_dofs = intersect(all_dofs, boundary_dofs_)
interior_dofs = setdiff(all_dofs, boundary_dofs_)
@@ -49,7 +64,7 @@ function reduce(assembly::Assembly, boundary_dofs_::Vector{Int})
empty!(assembly.force_vector)
gc()
if dim < 100000
if dim < mindofs
# no need to do any reduction of matrix size at all, just \ it.
return CAssembly([], all_dofs, Matrix{Float64}(), K, f, spzeros(0, 0), spzeros(0,1))
end
@@ -108,7 +123,7 @@ function reduce(assembly::Assembly, boundary_dofs_::Vector{Int})
end
info("Reduction: ", round(k/chunks*100, 0), " % done")
end
Kc = spzeros(dim, dim)
Kc[boundary_dofs, boundary_dofs] = Kbb - Kd
=#
+2
View File
@@ -31,6 +31,8 @@ using ForwardDiff
autodiffcache = ForwardDiffCache()
# export derivative, jacobian, hessian
using JLD
""" Simple linspace extension to arrays.
Examples
+76 -242
View File
@@ -3,6 +3,10 @@
## Direct solver
using JuliaFEM
@everywhere using JuliaFEM
@everywhere assemble = JuliaFEM.Core.assemble
type DirectSolver <: Solver
field_problems :: Vector{Problem}
boundary_problems :: Vector{BoundaryProblem}
@@ -10,6 +14,22 @@ type DirectSolver <: Solver
nonlinear_problem :: Bool
max_iterations :: Int64
tol :: Float64
dump_matrices :: Bool
reduce_stiffness_matrix :: Bool
end
""" Default initializer. """
function DirectSolver()
DirectSolver(
[], # field problems
[], # boundary problems
false, # parallel run?
true, # nonlinear problem?
10, # max nonlinear iterations
1.0e-6, # convergence tolerance
false, # dump matrices
true # reduce stiffness matrix
)
end
function push!(solver::DirectSolver, problem::Problem)
@@ -20,11 +40,6 @@ function push!(solver::DirectSolver, problem::BoundaryProblem)
push!(solver.boundary_problems, problem)
end
""" Default initializer. """
function DirectSolver()
DirectSolver([], [], false, true, 10, 1.0e-6)
end
function tic(timing, what::ASCIIString)
timing[what * " start"] = time()
end
@@ -37,205 +52,12 @@ function time_elapsed(timing, what::ASCIIString)
return timing[what * " finish"] - timing[what * " start"]
end
""" Call solver to solve a set of problems. """
function call(solver::DirectSolver, ::Type{Val{:noreduce}}, time::Number=0.0)
#@assert length(solver.field_problems) == 1
info("# of field problems: $(length(solver.field_problems))")
info("# of boundary problems: $(length(solver.boundary_problems))")
@assert solver.nonlinear_problem == true
timing = Dict{ASCIIString, Float64}()
tic(timing, "solver")
tic(timing, "initialization")
# check that all problems are "same kind"
field_name = get_unknown_field_name(solver.field_problems[1])
field_dim = get_unknown_field_dimension(solver.field_problems[1])
for field_problem in solver.field_problems
get_unknown_field_name(field_problem) == field_name || error("several different fields not supported yet")
get_unknown_field_dimension(field_problem) == field_dim || error("several different field dimensions not supported yet")
end
# create initial fields for this increment
# i.e., copy last known values as initial guess
# for this increment
for field_problem in solver.field_problems
for element in get_elements(field_problem)
gdofs = get_gdofs(element, field_dim)
if haskey(element, field_name)
if !isapprox(last(element[field_name]).time, time)
last_data = copy(last(element[field_name]).data)
push!(element[field_name], time => last_data)
end
else
data = Vector{Float64}[zeros(field_dim) for i in 1:length(element)]
element[field_name] = (time => data)
end
end
end
for boundary_problem in solver.boundary_problems
for element in get_elements(boundary_problem)
gdofs = get_gdofs(element, field_dim)
data = Vector{Float64}[zeros(field_dim) for i in 1:length(element)]
if haskey(element, "reaction force")
if !isapprox(last(element["reaction force"]).time, time)
push!(element["reaction force"], time => data)
end
else
element["reaction force"] = (time => data)
end
end
end
toc(timing, "initialization")
dim = 0
for iter=1:solver.max_iterations
info("Starting iteration $iter")
tic(timing, "non-linear iteration")
mapper = solver.parallel ? pmap : map
info("Assembling problems.")
# assemble boundary problems
tic(timing, "boundary assembly")
boundary_assembly = sum(mapper((p)->assemble(p, time), solver.boundary_problems))
boundary_dofs = unique(boundary_assembly.stiffness_matrix.I)
# boundary_dofs = collect(range(1, 12))
info("# of interface dofs: $(length(boundary_dofs))")
#info("dofs = $boundary_dofs")
toc(timing, "boundary assembly")
#static_condensation = false
# assemble field problems
dim = 0
assemblies = []
for (i, problem) in enumerate(solver.field_problems)
tic(timing, "field assembly")
field_assembly = assemble(problem, time)
#info("full assembly body $i")
#info(round(full(field_assembly.stiffness_matrix), 3))
toc(timing, "field assembly")
field_dofs = unique(field_assembly.stiffness_matrix.I)
#info("# of dofs in problem $i: $(length(field_dofs))")
dim = maximum([dim, maximum(field_dofs)])
#tic(timing, "condensate")
#cfield_assembly = condensate(field_assembly, boundary_dofs)
#toc(timing, "condensate")
#push!(assemblies, cfield_assembly)
push!(assemblies, field_assembly)
end
#info("dim = $dim")
#info("assembly done")
tic(timing, "create sparse matrices")
# create sparse matrices and saddle point problem
#K = sparse(field_assembly.stiffness_matrix)
#dim = size(K, 1)
#r = sparse(field_assembly.force_vector, dim, 1)
K = spzeros(dim, dim)
r = spzeros(dim, 1)
for (i, assembly) in enumerate(assemblies)
#info("body $i")
#info(round(full(assembly.Kc), 3))
#resize!(assembly.stiffness_matrix, dim, dim)
#resize!(assembly.force_vector, dim, 1)
K += sparse(assembly.stiffness_matrix, dim, dim)
r += sparse(assembly.force_vector, dim, 1)
end
#info(round(full(K), 3))
C = sparse(boundary_assembly.stiffness_matrix, dim, dim)
g = sparse(boundary_assembly.force_vector, dim, 1)
A = [K C'; C spzeros(dim, dim)]
b = [r; g]
toc(timing, "create sparse matrices")
#info("problem size = ", size(A))
info("Solving system")
tic(timing, "solution of system")
# solve increment for linearized problem
nz = unique(rowvals(A)) # take only non-zero rows
sol = zeros(length(b))
sol[nz] = full(A[nz,nz]) \ full(b[nz])
#info("solution vector before reconstruction")
#info(full(sol)')
la = sol[dim+1:end]
#for assembly in assemblies
# reconstruct!(assembly, sol)
#end
# la = vec(full(la))
# sol = vec(full(sol))
#info("la = ", la')
#info("sol = ", sol')
info("solved. solution norm: $(norm(sol[1:dim]))")
toc(timing, "solution of system")
info("Updating element data")
tic(timing, "update element data")
# update elements in field problems
for field_problem in solver.field_problems
for element in get_elements(field_problem)
gdofs = get_gdofs(element, field_dim)
local_sol = sol[gdofs] # incremental data for element
local_sol = reshape(local_sol, field_dim, length(element))
local_sol = Vector{Float64}[local_sol[:,i] for i=1:length(element)]
last(element[field_name]).data += local_sol # <-- added
end
end
# update elements in boundary problems
for boundary_problem in solver.boundary_problems
for element in get_elements(boundary_problem)
gdofs = get_gdofs(element, field_dim)
local_sol = la[gdofs]
local_sol = reshape(local_sol, field_dim, length(element))
local_sol = Vector{Float64}[local_sol[:,i] for i=1:length(element)]
last(element["reaction force"]).data = local_sol # <-- replaced
end
end
toc(timing, "update element data")
toc(timing, "non-linear iteration")
if true
info("timing info for non-linear iteration:")
info("boundary assembly : ", time_elapsed(timing, "boundary assembly"))
info("field assembly : ", time_elapsed(timing, "field assembly"))
# info("condensate : ", time_elapsed(timing, "condensate"))
info("create sparse matrices : ", time_elapsed(timing, "create sparse matrices"))
info("solution of system : ", time_elapsed(timing, "solution of system"))
info("update element data : ", time_elapsed(timing, "update element data"))
info("non-linear iteration : ", time_elapsed(timing, "non-linear iteration"))
end
if norm(sol[1:dim]) < solver.tol
toc(timing, "solver")
info("solver finished in ", time_elapsed(timing, "solver"), " seconds.")
return (iter, true)
end
end
info("Warning: did not coverge in $(solver.max_iterations) iterations!")
return (solver.max_iterations, false)
end
""" Call solver to solve a set of problems. """
function call(solver::DirectSolver, time::Number=0.0)
info("# of field problems: $(length(solver.field_problems))")
info("# of boundary problems: $(length(solver.boundary_problems))")
@assert solver.nonlinear_problem == true
timing = Dict{ASCIIString, Float64}()
tic(timing, "solver")
tic(timing, "initialization")
@@ -296,6 +118,10 @@ function call(solver::DirectSolver, time::Number=0.0)
boundary_assembly = sum(mapper((p)->assemble(p, time), solver.boundary_problems))
boundary_dofs = unique(boundary_assembly.stiffness_matrix.I)
info("# of interface dofs: $(length(boundary_dofs))")
C = sparse(boundary_assembly.stiffness_matrix)
g = sparse(boundary_assembly.force_vector)
boundary_assembly = nothing
gc()
toc(timing, "boundary assembly")
info("Assembling field problems...")
@@ -303,20 +129,34 @@ function call(solver::DirectSolver, time::Number=0.0)
assemblies = []
for (i, problem) in enumerate(solver.field_problems)
info("Assembling body $i...")
tic(timing, "field assembly")
field_assembly = assemble(problem, time)
nchunks = length(workers())
ne = length(get_elements(problem))
kk = round(Int, collect(linspace(0, ne, nchunks+1)))
slices = [kk[j]+1:kk[j+1] for j=1:length(kk)-1]
field_assembly = sum(pmap((s) -> assemble(problem, s, time), slices))
#field_assembly = assemble(problem, time)
toc(timing, "field assembly")
field_dofs = unique(field_assembly.stiffness_matrix.I)
info("# of dofs in problem $i: $(length(field_dofs))")
info("Eliminating interior dofs for body $i...")
dim = maximum([dim, maximum(field_dofs)])
tic(timing, "condensate")
cfield_assembly = reduce(field_assembly, boundary_dofs)
toc(timing, "condensate")
tic(timing, "reduce stiffness matrix")
cfield_assembly = nothing
if solver.reduce_stiffness_matrix && (nnz(boundary) != 0)
info("Eliminating interior dofs for body $i...")
cfield_assembly = reduce(field_assembly, boundary_dofs)
else
cfield_assembly = reduce(field_assembly, boundary_dofs, Inf)
end
toc(timing, "reduce stiffness matrix")
push!(assemblies, cfield_assembly)
end
tic(timing, "create sparse matrices")
K = spzeros(dim, dim)
f = spzeros(dim, 1)
@@ -327,45 +167,37 @@ function call(solver::DirectSolver, time::Number=0.0)
f += assembly.fc
end
C = sparse(boundary_assembly.stiffness_matrix, dim, dim)
g = sparse(boundary_assembly.force_vector, dim, 1)
A = [K C'; C spzeros(dim, dim)]
b = [f; g]
resize!(C, dim, dim)
resize!(g, dim, 1)
toc(timing, "create sparse matrices")
info("Solving interface system")
tic(timing, "solution of system")
nz = sort(unique(rowvals(A))) # take only non-zero rows
sol = zeros(b)
sol[nz] = A[nz,nz] \ full(b[nz])
toc(timing, "solution of system")
#=
try
catch
dump(round(full(A[nz,nz]), 3))
dump(round(full(b[nz]'), 3))
for (i, assembly) in enumerate(assemblies)
info("assembly $i dump")
dump(round(full(assembly.Kc), 3))
dump(round(full(assembly.fc), 3)')
end
info("matrix K")
dump(round(full(K), 3))
info("interface matrix")
dump(round(full(C), 3))
info("final assembly to solve:")
dump(round(full(A), 3))
dump(round(full(b'), 3))
info("nonzero dofs: $nz")
info("nonzero dofs removed:")
dump(round(full(A[nz,nz]), 3))
dump(round(full(b[nz]'), 3))
detsys = det(A[nz,nz])
info("determinant of system: $detsys")
error("Solving system failed.")
tic(timing, "dump matrices to disk")
if solver.dump_matrices
save("host_$(myid())_iteration_$(iter)_matrices.jld",
"stiffness matrix", K, "force vector", f,
"constraint matrix lhs", C, "constraint matrix rhs", g)
end
=#
toc(timing, "dump matrices to disk")
all_dofs = sort(unique(rowvals(K)))
field_dofs = setdiff(all_dofs, boundary_dofs)
info("Solving system")
tic(timing, "solution of system")
sol = zeros(2*dim)
if nnz(g) != 0
A = [K C'; C spzeros(dim, dim)]
b = [f; g]
K = 0
C = 0
gc()
nz = sort(unique(rowvals(A))) # take only non-zero rows
sol[nz] = A[nz,nz] \ full(b[nz])
else
K = 1/2*(K + K')
sol[field_dofs] = cholfact(K[field_dofs, field_dofs]) \ f[field_dofs]
end
toc(timing, "solution of system")
info("Solved, calculating interior dofs...")
tic(timing, "back substitute")
@@ -410,9 +242,11 @@ function call(solver::DirectSolver, time::Number=0.0)
info("timing info for non-linear iteration:")
info("boundary assembly : ", time_elapsed(timing, "boundary assembly"))
info("field assembly : ", time_elapsed(timing, "field assembly"))
info("reduce stiffness matrix : ", time_elapsed(timing, "condensate"))
info("reduce stiffness matrix : ", time_elapsed(timing, "reduce stiffness matrix"))
info("create sparse matrices : ", time_elapsed(timing, "create sparse matrices"))
info("solution of system : ", time_elapsed(timing, "solution of system"))
info("dump matrices to disk : ", time_elapsed(timing, "dump matrices to disk"))
info("solve problem : ", time_elapsed(timing, "solution of system"))
info("back substitute : ", time_elapsed(timing, "back substitute"))
info("update element data : ", time_elapsed(timing, "update element data"))
info("non-linear iteration : ", time_elapsed(timing, "non-linear iteration"))
end
+1 -1
View File
@@ -42,7 +42,7 @@ References
https://en.wikipedia.org/wiki/Heat_equation
"""
function assemble!{E<:CG}(assembly::Assembly, problem::Problem{HeatProblem}, element::Element{E}, time::Number)
function assemble!(assembly::Assembly, problem::Problem{HeatProblem}, element::Element, time::Number)
gdofs = get_gdofs(element, problem.dim)
for ip in get_integration_points(element)
+7
View File
@@ -10,6 +10,8 @@ type SparseMatrixIJV
V :: Vector{Float64}
end
typealias SparseMatrixCOO SparseMatrixIJV
function SparseMatrixIJV()
SparseMatrixIJV([], [], [])
end
@@ -89,3 +91,8 @@ function add!(A::SparseMatrixIJV, dofs::Vector{Int}, data::Array{Float64})
append!(A.V, vec(data))
end
function optimize!(A::SparseMatrixIJV)
I, J, V = findnz(sparse(A))
A = SparseMatrixCOO(I, J, V)
gc()
end