diff --git a/notebooks/2015-06-25-elasticity-solver-example.ipynb b/notebooks/2015-06-25-elasticity-solver-example.ipynb index 076f47c..04ef7d5 100644 --- a/notebooks/2015-06-25-elasticity-solver-example.ipynb +++ b/notebooks/2015-06-25-elasticity-solver-example.ipynb @@ -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)" ] }, diff --git a/src/assembly.jl b/src/assembly.jl index b1dfb92..f5ca3a5 100644 --- a/src/assembly.jl +++ b/src/assembly.jl @@ -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 =# diff --git a/src/core.jl b/src/core.jl index 0cf93ce..a5fb9ca 100644 --- a/src/core.jl +++ b/src/core.jl @@ -31,6 +31,8 @@ using ForwardDiff autodiffcache = ForwardDiffCache() # export derivative, jacobian, hessian +using JLD + """ Simple linspace extension to arrays. Examples diff --git a/src/directsolver.jl b/src/directsolver.jl index 4148dfb..93edffc 100644 --- a/src/directsolver.jl +++ b/src/directsolver.jl @@ -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 diff --git a/src/heat.jl b/src/heat.jl index b72cc73..3f1537b 100644 --- a/src/heat.jl +++ b/src/heat.jl @@ -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) diff --git a/src/sparse.jl b/src/sparse.jl index 95a0390..54b2487 100644 --- a/src/sparse.jl +++ b/src/sparse.jl @@ -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