From d4e21af8debdc3dba9f7876b789f42592aa27c12 Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Sat, 5 Dec 2015 10:58:01 +0200 Subject: [PATCH] petsc interface --- src/directsolver.jl | 42 +++++++++-------- src/elements.jl | 4 +- src/integrate.jl | 2 +- src/interfaces.jl | 20 +++----- src/mortar.jl | 6 ++- src/petsc.jl | 98 +++++++++++++++++++++++++++++++++++++++ test/test_directsolver.jl | 7 +-- test/test_mortar.jl | 12 ++++- test/test_solver.jl | 48 ++++++++++++++++++- 9 files changed, 194 insertions(+), 45 deletions(-) create mode 100644 src/petsc.jl diff --git a/src/directsolver.jl b/src/directsolver.jl index a25f1d6..3912118 100644 --- a/src/directsolver.jl +++ b/src/directsolver.jl @@ -1,6 +1,12 @@ # This file is a part of JuliaFEM. # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md +#= Solution norms for piston model + +piston_19611_P2.inp iter 1 2.048090408266966 + +=# + ## Direct solver using JuliaFEM @@ -32,7 +38,7 @@ function DirectSolver(name="DirectSolver") 1.0e-6, # convergence tolerance false, # dump matrices true, # reduce stiffness matrix - :LDLt # method: LDLt or LU ? + :CHOLMOD # method: CHOLMOD, UMFPACK, PETSc_GMRES ) end @@ -64,10 +70,9 @@ Solve problem Cu = g """ -function solve(K, f, C, g, ::Type{Val{:LDLt}}) +function solve(K, f, C, g, ::Type{Val{:CHOLMOD}}) t0 = time() - # make sure K is symmetric # K = Symmetric(K) s = maximum(abs(1/2*(K + K') - K)) @@ -84,55 +89,52 @@ function solve(K, f, C, g, ::Type{Val{:LDLt}}) all_dofs = unique(rowvals(K)) interior_dofs = setdiff(all_dofs, boundary_dofs) - info("all dofs = $(length(all_dofs))") - info("interior dofs = $(length(interior_dofs))") - info("boundary dofs = $(length(boundary_dofs))") - info("preparation in ", time()-t0, " seconds") + info("CHOLMOD: all dofs = $(length(all_dofs))") + info("CHOLMOD: interior dofs = $(length(interior_dofs))") + info("CHOLMOD: boundary dofs = $(length(boundary_dofs))") # solve displacement on known boundary - t0 = time() LUF = lufact(C[boundary_dofs, boundary_dofs]) u = zeros(dim) u[boundary_dofs] = LUF \ full(g[boundary_dofs]) - info("displacement on boundary solved.") + info("CHOLMOD: displacement on boundary solved.") normub = norm(u[boundary_dofs]) - info("norm[u_boundary_dofs] = ", normub) if isapprox(normub, 0.0) - info("homogeneous dirichlet boundary") + info("CHOLMOD: homogeneous dirichlet boundary") end - info("solve boundary = ", time()-t0) # factorize interior domain using cholmod - t0 = time() + t = time() CF = cholfact(K[interior_dofs, interior_dofs]) Kib = K[interior_dofs, boundary_dofs] Kbb = K[boundary_dofs, boundary_dofs] fi = f[interior_dofs] - info("factorizations done in ", time()-t0, " seconds") + info("CHOLMOD: LDLt factorization done in ", time()-t, " seconds") # solve interior domain + lagrange multipliers - t0 = time() u[interior_dofs] = CF \ (fi - Kib*u[boundary_dofs]) la = zeros(dim) la[boundary_dofs] = LUF \ full(Kib'*u[interior_dofs] - Kbb*u[boundary_dofs]) - info("solved interior in ", time()-t0, " seconds. norm = ", norm(u)) + info("CHOLMOD: solved in ", time()-t0, " seconds. norm = ", norm(u)) return u, la end -function solve(K, f, C, g, ::Type{Val{:LU}}) +function solve(K, f, C, g, ::Type{Val{:UMFPACK}}) + t0 = time() dim = size(K, 1) A = nothing try A = [K C'; C spzeros(dim, dim)] catch - info("size(K) = ", size(K)) - info("size(C) = ", size(C)) - error("Failed to construct problem. dim = $dim") + info("UMFPACK: size(K) = ", size(K)) + info("UMFPACK: size(C) = ", size(C)) + error("UMFPACK: Failed to construct problem. dim = $dim") end b = [f; g] nz = sort(unique(rowvals(A))) u = zeros(length(b)) u[nz] = lufact(A[nz,nz]) \ full(b[nz]) + info("UMFPACK: solved in ", time()-t0, " seconds. norm = ", norm(u[1:dim])) return u[1:dim], u[dim+1:end] end diff --git a/src/elements.jl b/src/elements.jl index fbdfdf4..2cff190 100644 --- a/src/elements.jl +++ b/src/elements.jl @@ -13,8 +13,8 @@ function convert{E}(::Type{Element{E}}, connectivity::Vector{Int}) return Element{E}(connectivity, Dict()) end -function get_integration_points{E}(element::Element{E}) - return get_integration_points(E) +function get_integration_points{E}(element::Element{E}, args...) + return get_integration_points(E, args...) end function update_gauss_fields!(element::Element, data::Vector{IntegrationPoint}, time::Real) diff --git a/src/integrate.jl b/src/integrate.jl index 16aed6a..102e00b 100644 --- a/src/integrate.jl +++ b/src/integrate.jl @@ -36,7 +36,7 @@ function get_integration_points(::Type{Seg3}, ::Type{Val{4}}) ] end -function get_integration_points(::Type{Seg3}, ::Type{Val{5}}) +function get_integration_points(::Union{Type{Seg2}, Type{Seg3}}, ::Type{Val{5}}) [ IntegrationPoint([-1/3*sqrt(5 + 2*sqrt(10/7))], (322-13*sqrt(70))/900), IntegrationPoint([-1/3*sqrt(5 - 2*sqrt(10/7))], (322+13*sqrt(70))/900), diff --git a/src/interfaces.jl b/src/interfaces.jl index b65082e..2491bb0 100644 --- a/src/interfaces.jl +++ b/src/interfaces.jl @@ -8,7 +8,7 @@ using JuliaFEM.API: Model Function for creating solver and all the necessary components for the calculation """ -function get_solver(model::Model, case_name::ASCIIString, time::Float64) +function get_solver(model::Model, case_name::ASCIIString) case = model.load_cases[case_name] # Create core elements @@ -41,15 +41,9 @@ function create_solver(model, case, core_elements, dirichlet_arr) end # Creating the solver and pushing problems and - # boundary conditions -# if case.solver == :LinearSolver -# solver = JuliaFEM.Core.(case.solver)(field_problem, -# dirichlet_arr...) -# else - solver = JuliaFEM.Core.(case.solver)() - push!(solver, field_problem) - push!(solver, dirichlet_arr...) -# end + solver = JuliaFEM.Core.(case.solver)() + push!(solver, field_problem) + push!(solver, dirichlet_arr...) return solver end @@ -76,7 +70,7 @@ function create_dirichlet_bcs(model, case, core_elements) end push!(dirichlet_arr, problem) end - dirichlet_arr + return dirichlet_arr end """ @@ -100,7 +94,7 @@ function create_core_elements(model) core_elements[el_id] = core_element model.elements[el_id].results = core_element end - core_elements + return core_elements end """ @@ -126,7 +120,7 @@ end """ function solve!(model::Model, case_name::ASCIIString, time::Float64) # Create solver - solver = get_solver(model, case_name, time) + solver = get_solver(model, case_name) # Solve problem at given time solver(time) diff --git a/src/mortar.jl b/src/mortar.jl index 30a5a8d..7820087 100644 --- a/src/mortar.jl +++ b/src/mortar.jl @@ -159,7 +159,7 @@ function assemble!(assembly::Assembly, problem::BoundaryProblem{MortarProblem}, continue # no contribution end master_dofs = get_gdofs(master_element, field_dim) - for ip in get_integration_points(slave_element) + for ip in get_integration_points(slave_element, Val{5}) w = ip.weight*det(slave_element, ip, time)*l # integration point on slave side segment @@ -171,7 +171,9 @@ function assemble!(assembly::Assembly, problem::BoundaryProblem{MortarProblem}, N1 = slave_element(xi_gauss, time) N2 = master_element(xi_projected, time) S = w*N1'*N1 - M = w*N1'*N2 + M = w*(N1'*N2)' + # FIXME: why this needs now to be transpose? + # assembly / repeat for i=1:field_dim sd = slave_dofs[i:field_dim:end] md = master_dofs[i:field_dim:end] diff --git a/src/petsc.jl b/src/petsc.jl new file mode 100644 index 0000000..ba0cbde --- /dev/null +++ b/src/petsc.jl @@ -0,0 +1,98 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md + +# PETSc interface for solver + +using PETSc + +import JuliaFEM.Core: solve + +""" +Parameters +---------- +preconditioner : "jacobi" +ksp_type: "bcgs", "gmres"? +""" +function solve(K, f, C, g, ::Type{Val{:PETSc_GMRES}}; preconditioner=nothing) + t0 = time() + dim = size(K, 1) + + # make sure C is square + boundary_dofs = unique(rowvals(C)) + boundary_dofs2 = unique(rowvals(C')) + @assert length(boundary_dofs) == length(boundary_dofs2) + @assert setdiff(Set(boundary_dofs), Set(boundary_dofs2)) == Set() + all_dofs = unique(rowvals(K)) + interior_dofs = setdiff(all_dofs, boundary_dofs) + info("PETSc: all dofs = $(length(all_dofs))") + info("PETSc: interior dofs = $(length(interior_dofs))") + info("PETSc: boundary dofs = $(length(boundary_dofs))") + # solve displacement on known boundary + LUF = lufact(C[boundary_dofs, boundary_dofs]) + u = zeros(dim) + u[boundary_dofs] = LUF \ full(g[boundary_dofs]) + info("PETSc: displacement on boundary solved.") + normub = norm(u[boundary_dofs]) + if isapprox(normub, 0.0) + info("PETSc: homogeneous dirichlet boundary") + end + + # interior domain and lagrange multipliers + + t = time() + # this is completely unnecessary step and will be removed in future. + # --> + info("PETSc: creating matrices in PETSc format.") + ninterior_dofs = length(interior_dofs) + + # nz, see https://github.com/JuliaParallel/PETSc.jl/issues/52 + d = Dict{Int64, Int64}() + for i in rowvals(K) + haskey(d, i) ? (d[i] += 1) : (d[i] = 1) + end + nz = maximum(values(d)) + + A = PETSc.Mat(Float64, ninterior_dofs, ninterior_dofs; nz=nz) + info("PETSc: $ninterior_dofs interior dofs, assembling to PETSc Mat") + for (i, j, v) in zip(findnz(K[interior_dofs, interior_dofs])...) + A[i, j] = v + end + + fi = f[interior_dofs] + + b = PETSc.Vec(Float64, ninterior_dofs, PETSc.C.VECMPI) + for (i, j, v) in zip(findnz(sparse(f[interior_dofs]))...) + b[i] = v + end + + info("PETSc: initialization of matrices in ", time()-t, " seconds") + # <-- + + kspg = PETSc.KSP(A, ksp_monitor="") + + # apply preconditioner if defined + if !isa(preconditioner, Void) + info("PETSc: preconditioner: $preconditioner") + pc = PETSc.PC(Float64, comm=PETSc.comm(kspg), pc_type=preconditioner) + PETSc.chk(PETSc.C.PCSetOperators(pc.p, A.p, A.p)) + kspg = PETSc.KSP(pc, ksp_monitor="") + end + + info("PETSc: performing ksp GMRES solve") + x = kspg \ b + info("PETSc: finished ksp solve") + info("PETSc: ksp info:\n",petscview(kspg)) + for (i, d) in enumerate(interior_dofs) + u[d] = x[i] + end + + la = zeros(dim) + Kib = K[interior_dofs, boundary_dofs] + Kbb = K[boundary_dofs, boundary_dofs] + la[boundary_dofs] = LUF \ full(Kib'*u[interior_dofs] - Kbb*u[boundary_dofs]) + + info("PETSc: solved in ", time()-t0, " seconds. norm = ", norm(u)) + return u, la +end + +info("PETSc interface loaded.") diff --git a/test/test_directsolver.jl b/test/test_directsolver.jl index f4585af..51797b9 100644 --- a/test/test_directsolver.jl +++ b/test/test_directsolver.jl @@ -98,10 +98,10 @@ function test_direct_cholesky_with_non_homogeneous_dirichlet_conditions() push!(solver, boundary) # launch solver -# solver.method = :LU -# solver.dump_matrices = true + solver.method = :UMFPACK + solver.dump_matrices = true solver.max_iterations = 1 - solver(0.0) + iters, status = solver(0.0) # FIXME: solver gives no convergence warning when all dofs are fixed. n1disp = e1("displacement", [-1.0, -1.0], 0.0) n2disp = e1("displacement", [ 1.0, -1.0], 0.0) @@ -113,6 +113,7 @@ function test_direct_cholesky_with_non_homogeneous_dirichlet_conditions() @test isapprox(n3disp, [-0.1, 0.1]) @test isapprox(n2disp, [ 0.2, -0.2]) @test isapprox(n4disp, [ 0.2, -0.2]) + @test status == true end #test_direct_cholesky_with_non_homogeneous_dirichlet_conditions() diff --git a/test/test_mortar.jl b/test/test_mortar.jl index c4f7b00..4bb88ad 100644 --- a/test/test_mortar.jl +++ b/test/test_mortar.jl @@ -23,6 +23,13 @@ function get_test_2d_model() master2 = Seg2([8, 9]) master2["geometry"] = Vector[N[8], N[9]] +#= + master1 = Seg2([9, 8]) + master1["geometry"] = Vector[N[9], N[8]] + master2 = Seg2([8, 7]) + master2["geometry"] = Vector[N[8], N[7]] +=# + slave1 = Seg2([10, 11]) slave1["geometry"] = Vector[N[10], N[11]] # should be n = [0 -1]' and t = [1 0]' @@ -202,8 +209,9 @@ function test_2d_mortar_multiple_bodies_multiple_dirichlet_bc() push!(solver, boundary2) push!(solver, boundary3) + solver.name = "test_2d_mortar_multiple_bodies_multiple_dirichlet_bcs" solver.dump_matrices = true - solver.method = :LU + solver.method = :UMFPACK # launch solver solver(0.0) @@ -324,7 +332,7 @@ function test_2d_mortar_three_bodies_shared_nodes() push!(solver, bc5) # launch solver - solver.method = :LU + solver.method = :UMFPACK call(solver, 0.0) disp = e2("displacement", [1.0, 1.0], 0.0) diff --git a/test/test_solver.jl b/test/test_solver.jl index 2070b89..4c5ef8c 100644 --- a/test/test_solver.jl +++ b/test/test_solver.jl @@ -8,6 +8,7 @@ using JuliaFEM.Test using JuliaFEM.Core: Seg2, Quad4 using JuliaFEM.Core: DirichletProblem, HeatProblem using JuliaFEM.Core: LinearSolver +using JuliaFEM.Core: solve function test_linearsolver() el1 = Quad4([1, 2, 3, 4]) @@ -35,7 +36,9 @@ function test_linearsolver() push!(boundary_problem, el3) # Create a solver for a set of problems - solver = LinearSolver(field_problem, boundary_problem) + solver = LinearSolver() + push!(solver, field_problem) + push!(solver, boundary_problem) # Solve problem at time t=1.0 and update fields solver(1.0) @@ -48,7 +51,48 @@ function test_linearsolver() info("Temperature at point X = $X is T = $T") @test isapprox(T, 100.0) end +#test_basic() -test_basic() +function test_solvers() + K = [ + 440.0 150.0 -260.0 -30.0 40.0 30.0 -220.0 -150.0 + 150.0 440.0 30.0 40.0 -30.0 -260.0 -150.0 -220.0 + -260.0 30.0 440.0 -150.0 -220.0 150.0 40.0 -30.0 + -30.0 40.0 -150.0 440.0 150.0 -220.0 30.0 -260.0 + 40.0 -30.0 -220.0 150.0 440.0 -150.0 -260.0 30.0 + 30.0 -260.0 150.0 -220.0 -150.0 440.0 -30.0 40.0 + -220.0 -150.0 40.0 30.0 -260.0 -30.0 440.0 150.0 + -150.0 -220.0 -30.0 -260.0 30.0 40.0 150.0 440.0] + + C = 1/3*[ + 1.0 0.0 0.0 0.0 0.5 0.0 0.0 0.0 + 0.0 1.0 0.0 0.0 0.0 0.5 0.0 0.0 + 0.0 0.0 1.0 0.0 0.0 0.0 0.5 0.0 + 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.5 + 0.5 0.0 0.0 0.0 1.0 0.0 0.0 0.0 + 0.0 0.5 0.0 0.0 0.0 1.0 0.0 0.0 + 0.0 0.0 0.5 0.0 0.0 0.0 1.0 0.0 + 0.0 0.0 0.0 0.5 0.0 0.0 0.0 1.0] + + f = zeros(8) + + g = 1/100 * [-5.0, 5.0, 10.0, -10.0, -5.0, 5.0, 10.0, -10.0] + + K = sparse(K) + C = sparse(C) + f = sparse(f) + g = sparse(g) + + expected = [-0.1, 0.1, 0.2, -0.2, -0.1, 0.1, 0.2, -0.2] + u1, la1 = solve(K, f, C, g, Val{:UMFPACK}) + @test isapprox(u1, expected) + u2, la2 = solve(K, f, C, g, Val{:CHOLMOD}) + @test isapprox(u2, expected) + include(Pkg.dir("JuliaFEM"*"/src/petsc.jl")) + u3, la3 = solve(K, f, C, g, Val{:PETSc_GMRES}) + @test isapprox(u3, expected) +end + +# test_solvers() end