mirror of
https://github.com/JuliaFEM/JuliaFEM.jl.git
synced 2026-09-19 09:54:55 +00:00
petsc interface
This commit is contained in:
Jukka Aho
+22
-20
@@ -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
|
||||
|
||||
|
||||
+2
-2
@@ -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)
|
||||
|
||||
+1
-1
@@ -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),
|
||||
|
||||
+7
-13
@@ -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)
|
||||
|
||||
+4
-2
@@ -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]
|
||||
|
||||
@@ -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.")
|
||||
@@ -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()
|
||||
|
||||
|
||||
+10
-2
@@ -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)
|
||||
|
||||
+46
-2
@@ -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
|
||||
|
||||
Reference in New Issue
Block a user