diff --git a/Project.toml b/Project.toml index eaba750..fcf9155 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "FinEtoolsDDMethods" uuid = "093fdad3-4fce-4574-aa49-c2a24de4f00d" authors = ["Petr Krysl "] -version = "0.3.2" +version = "0.4.0" [deps] CoNCMOR = "3a989f29-cff0-4f1c-87be-5dcdd41bd240" diff --git a/README.md b/README.md index ed4d7b1..abb7758 100644 --- a/README.md +++ b/README.md @@ -4,9 +4,23 @@ Domain decomposition methods: solvers and algorithms used with FinEtools applications. +## Capabilities and limitations + +Only linear problems are addressed at the moment. `FinEtools` discretization only. + +- Schur-based decomposition for conjugate gradient iteration. +- Overlapping Schwarz two-level preconditioning of conjugate gradient iteration. + +Both of the above can now be expressed as MPI-parallel algorithms. + +## Examples + Please refer to the `README` in the `examples` folder for instructions on how to get the examples to run. +## References + Paper [submitted](https://papers.ssrn.com/sol3/papers.cfm?abstract_id=4902156) for the coherent node cluster conjugate gradient preconditioner with applications -to almost incompressible elasticity and general shells. \ No newline at end of file +to almost incompressible elasticity and general shells. + diff --git a/examples/README.md b/examples/README.md index c1a6ed1..5be9670 100644 --- a/examples/README.md +++ b/examples/README.md @@ -33,7 +33,8 @@ On Windows 11, the following would work: mpiexec -n 3 julia --project=. .\conc\heat\Poisson2D_overlapped_mpi_examples.jl ``` -## How to run a CoNC example + +## How to run a CoNC (sequential) example Two classes of problems solved here: - Three dimensional elasticity (composite of matrix with embedded fibers). @@ -105,4 +106,15 @@ Run ``` $ julia conc/shells/barrel.jl --help ``` -to see the available options. \ No newline at end of file +to see the available options. + +## How to run a CoNC (MPI-parallel) example + +At the moment only the heat conduction examples have been cast in this form. Try +``` +mpiexec -n 5 julia --project=. .\conc\heat\Poisson2D_overlapped_mpi_examples.jl +``` + +## To do + +- MPI parallelization: When making the stiffness matrices of the partitions, for the shell it will be necessary to associate the geometry (normals!) with the entire geometry. diff --git a/examples/conc/heat/Poisson2D_overlapped_mpi_examples.jl b/examples/conc/heat/Poisson2D_overlapped_mpi_examples.jl index fbd7d5e..5d3bcc3 100644 --- a/examples/conc/heat/Poisson2D_overlapped_mpi_examples.jl +++ b/examples/conc/heat/Poisson2D_overlapped_mpi_examples.jl @@ -69,31 +69,34 @@ function _execute(N, mesher, volrule, nelperpart, nbf1max, overlap, itmax, relre List = vcat(l1, l2, l3, l4, ) setebc!(Temp, List, true, 1, tempf(geom.values[List, :])[:]) numberdofs!(Temp) - fr = dofrange(Temp, DOF_KIND_FREE) - dr = dofrange(Temp, DOF_KIND_DATA) - rank == 0 && @info("Number of free degrees of freedom: $(nfreedofs(Temp)) ($(round(time() - t1, digits=3)) [s])") - t1 = time() + material = MatHeatDiff(thermal_conductivity) - femm = FEMMHeatDiff(IntegDomain(fes, volrule, 1.0), material) - K = conductivity(femm, geom, Temp) - rank == 0 && @info("Conductivity ($(round(time() - t1, digits=3)) [s])") - K_ff = K[fr, fr] - t1 = time() - fi = ForceIntensity(Float64[Q]) - F1 = distribloads(femm, geom, Temp, fi, 3) - rank == 0 && @info("Internal heat generation ($(round(time() - t1, digits=3)) [s])") - t1 = time() - T_d = gathersysvec(Temp, DOF_KIND_DATA) - F_f = (F1 .- K[:, dr] * T_d)[fr] - rank == 0 && @info("Right hand side ($(round(time() - t1, digits=3)) [s])") + + # fr = dofrange(Temp, DOF_KIND_FREE) + # dr = dofrange(Temp, DOF_KIND_DATA) + # rank == 0 && @info("Number of free degrees of freedom: $(nfreedofs(Temp)) ($(round(time() - t1, digits=3)) [s])") + # t1 = time() + + # femm = FEMMHeatDiff(IntegDomain(fes, volrule, 1.0), material) + # K = conductivity(femm, geom, Temp) + # rank == 0 && @info("Conductivity ($(round(time() - t1, digits=3)) [s])") + # K_ff = K[fr, fr] + # t1 = time() + # fi = ForceIntensity(Float64[Q]) + # F1 = distribloads(femm, geom, Temp, fi, 3) + # rank == 0 && @info("Internal heat generation ($(round(time() - t1, digits=3)) [s])") + # t1 = time() + # T_d = gathersysvec(Temp, DOF_KIND_DATA) + # F_f = (F1 .- K[:, dr] * T_d)[fr] + # rank == 0 && @info("Right hand side ($(round(time() - t1, digits=3)) [s])") t1 = time() cpartitioning, ncpartitions = FinEtoolsDDMethods.cluster_partitioning(fens, fes, fes.label, nelperpart) rank == 0 && @info("Number of clusters (coarse grid partitions): $(ncpartitions)") mor = CoNCData(fens, cpartitioning) Phi = transfmatrix(mor, LegendreBasis, nbf1max, Temp) - Phi_f = Phi[fr, :] - rank == 0 && @info("Size of the reduced problem: $(size(Phi_f, 2))") + Phi = Phi[freedofs(Temp), :] + rank == 0 && @info("Size of the reduced problem: $(size(Phi, 2))") rank == 0 && @info "Clustering ($(round(time() - t1, digits=3)) [s])" function make_matrix(fes) @@ -101,25 +104,31 @@ function _execute(N, mesher, volrule, nelperpart, nbf1max, overlap, itmax, relre return conductivity(femm2, geom, Temp) end + function make_interior_load(fes) + femm2 = FEMMHeatDiff(IntegDomain(fes, volrule, 1.0), material) + fi = ForceIntensity(Float64[Q]) + return distribloads(femm2, geom, Temp, fi, 3) + end + t1 = time() partition = nothing if rank > 0 cpi = CoNCPartitioningInfo(fens, fes, nfpartitions, overlap, Temp) - partition = CoNCPartitionData(cpi, rank, fes, Phi, make_matrix) + partition = CoNCPartitionData(cpi, rank, fes, Phi, make_matrix, make_interior_load) end + MPI.Barrier(comm) rank == 0 && (@info "Create partitions time: $(time() - t1)") t1 = time() - F_f_1 = zeros(size(K_ff, 2)) + F_f = zeros(nfreedofs(Temp)) if rank > 0 - F_f_1 .= partition.data_rhs + F_f .= partition.rhs end - MPI.Reduce!(F_f_1, MPI.SUM, comm; root=0) # Reduce the data-determined rhs - rank == 0 && (@show norm(F_f_1 + (K[:, dr] * T_d)[fr])) + MPI.Reduce!(F_f, MPI.SUM, comm; root=0) # Reduce the data-determined rhs rank == 0 && (@info "Compute RHS: $(time() - t1)") t1 = time() - Kr_ff = spzeros(size(Phi_f, 2), size(Phi_f, 2)) + Kr_ff = spzeros(size(Phi, 2), size(Phi, 2)) if rank > 0 Kr_ff = partition.reduced_K end @@ -140,7 +149,7 @@ function _execute(N, mesher, volrule, nelperpart, nbf1max, overlap, itmax, relre (q, p) -> partition_multiply!(q, partition, p), F_f, zeros(size(F_f)), - (q, p) -> precondition_global_solve!(q, Krfactor, Phi_f, p), + (q, p) -> precondition_global_solve!(q, Krfactor, Phi, p), (q, p) -> precondition_local_solve!(q, partition, p); itmax=itmax, atol=relrestol * norm_F_f, rtol=0, peeksolution=peeksolution) diff --git a/examples/conc/lindef/fibers_overlapped_examples.jl b/examples/conc/lindef/fibers_overlapped_examples.jl index 488965f..49903ab 100644 --- a/examples/conc/lindef/fibers_overlapped_examples.jl +++ b/examples/conc/lindef/fibers_overlapped_examples.jl @@ -6,6 +6,7 @@ using FinEtools.MeshTetrahedronModule: tetv using FinEtoolsDeforLinear using FinEtoolsDDMethods using FinEtoolsDDMethods.CGModule: pcg_seq +using FinEtoolsDDMethods.CompatibilityModule using SymRCM using Metis using Test @@ -343,7 +344,7 @@ function _execute(label, kind, Em, num, Ef, nuf, nelperpart, nbf1max, nfpartitio # VTK.vtkexportmesh("fibers-tet-red-sol.vtk", fens, fes; vectors=[("u", deepcopy(u.values),)]) # fpartitioning, nfpartitions = fine_grid_partitioning(fens, nfpartitions) - nodelists = fine_grid_node_lists(fens, fes, nfpartitions, overlap) + nodelists = CompatibilityModule.fine_grid_node_lists(fens, fes, nfpartitions, overlap) @assert length(nodelists) == nfpartitions partitions = [] diff --git a/examples/conc/shells/LE5_Z_cantilever_overlapped_examples.jl b/examples/conc/shells/LE5_Z_cantilever_overlapped_examples.jl index 2bf0fa6..547a334 100644 --- a/examples/conc/shells/LE5_Z_cantilever_overlapped_examples.jl +++ b/examples/conc/shells/LE5_Z_cantilever_overlapped_examples.jl @@ -26,7 +26,8 @@ using FinEtoolsFlexStructures.FEMMShellT3FFModule using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_field! using FinEtoolsDDMethods using FinEtoolsDDMethods.CGModule: pcg_seq -using FinEtoolsDDMethods.PartitionCoNCDDModule: patch_coordinates +using FinEtoolsDDMethods.CoNCUtilitiesModule: patch_coordinates +using FinEtoolsDDMethods.CompatibilityModule using SymRCM using Metis using Test @@ -179,7 +180,7 @@ function _execute(ncoarse, aspect, nelperpart, nbf1max, nfpartitions, overlap, r # VTK.vtkexportmesh("fibers-tet-red-sol.vtk", fens, fes; vectors=[("u", deepcopy(u.values),)]) # fpartitioning, nfpartitions = fine_grid_partitioning(fens, nfpartitions) - nodelists = fine_grid_node_lists(fens, fes, nfpartitions, overlap) + nodelists = CompatibilityModule.fine_grid_node_lists(fens, fes, nfpartitions, overlap) @assert length(nodelists) == nfpartitions partitions = [] diff --git a/examples/conc/shells/barrel_ilu_examples.jl b/examples/conc/shells/barrel_ilu_examples.jl index 75c15a8..f57448a 100644 --- a/examples/conc/shells/barrel_ilu_examples.jl +++ b/examples/conc/shells/barrel_ilu_examples.jl @@ -17,7 +17,7 @@ using FinEtoolsFlexStructures.FEMMShellT3FFModule using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_field! using FinEtoolsDDMethods using FinEtoolsDDMethods.CGModule: pcg_seq -using FinEtoolsDDMethods.PartitionCoNCDDModule: patch_coordinates +using FinEtoolsDDMethods.CoNCUtilitiesModule: patch_coordinates using SymRCM using Metis using Test diff --git a/examples/conc/shells/barrel_overlapped_examples.jl b/examples/conc/shells/barrel_overlapped_examples.jl index 35f7391..d86c7ff 100644 --- a/examples/conc/shells/barrel_overlapped_examples.jl +++ b/examples/conc/shells/barrel_overlapped_examples.jl @@ -17,7 +17,8 @@ using FinEtoolsFlexStructures.FEMMShellT3FFModule using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_field! using FinEtoolsDDMethods using FinEtoolsDDMethods.CGModule: pcg_seq -using FinEtoolsDDMethods.PartitionCoNCDDModule: patch_coordinates +using FinEtoolsDDMethods.CoNCUtilitiesModule: patch_coordinates +using FinEtoolsDDMethods.CompatibilityModule using SymRCM using Metis using Test @@ -194,7 +195,7 @@ function _execute(ncoarse, nelperpart, nbf1max, nfpartitions, overlap, ref, itma # VTK.vtkexportmesh("fibers-tet-red-sol.vtk", fens, fes; vectors=[("u", deepcopy(u.values),)]) - nodelists = fine_grid_node_lists(fens, fes, nfpartitions, overlap) + nodelists = CompatibilityModule.fine_grid_node_lists(fens, fes, nfpartitions, overlap) @assert length(nodelists) == nfpartitions partitions = [] @@ -224,34 +225,6 @@ function _execute(ncoarse, nelperpart, nbf1max, nfpartitions, overlap, ref, itma q end - # peeksolution(iter, x) = begin - # println("Iteration: $(iter)") - # if iter % 10 == 0 - # if iter > 10 - # f = "barrel_overlapped" * - # "-rf=$(ref)" * - # "-ne=$(nelperpart)" * - # "-n1=$(nbf1max)" * - # "-nf=$(nfpartitions)" * - # "-ov=$(overlap)" * - # "-cg-sol-iter$(iter-10)" - # xp = DataDrop.retrieve_matrix(f * ".h5", "x") - # scattersysvec!(dchi, x - xp) - # VTK.vtkexportmesh("barrel-xincr-iter=$(iter).vtk", fens, fes; - # vectors=[("u", deepcopy(dchi.values[:, 1:3]),)]) - # end - # f = "barrel_overlapped" * - # "-rf=$(ref)" * - # "-ne=$(nelperpart)" * - # "-n1=$(nbf1max)" * - # "-nf=$(nfpartitions)" * - # "-ov=$(overlap)" * - # "-cg-sol-iter$(iter)" - # DataDrop.empty_hdf5_file(f * ".h5") - # DataDrop.store_matrix(f * ".h5", "x", x) - # end - # end - peeksolution(iter, x, resnorm) = begin println("Iteration: $(iter)") println("Residual Norm: $(resnorm)") diff --git a/examples/conc/shells/cos_2t_press_hyperboloid_free_overlapped_examples.jl b/examples/conc/shells/cos_2t_press_hyperboloid_free_overlapped_examples.jl index f7f6fbb..0d95929 100644 --- a/examples/conc/shells/cos_2t_press_hyperboloid_free_overlapped_examples.jl +++ b/examples/conc/shells/cos_2t_press_hyperboloid_free_overlapped_examples.jl @@ -8,7 +8,8 @@ using FinEtoolsFlexStructures.FEMMShellT3FFModule using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_field! using FinEtoolsDDMethods using FinEtoolsDDMethods.CGModule: pcg_seq -using FinEtoolsDDMethods.PartitionCoNCDDModule: patch_coordinates +using FinEtoolsDDMethods.CoNCUtilitiesModule: patch_coordinates +using FinEtoolsDDMethods.CompatibilityModule using SymRCM using Metis using Test @@ -149,7 +150,7 @@ function _execute(ncoarse, aspect, nelperpart, nbf1max, nfpartitions, overlap, r # VTK.vtkexportmesh("fibers-tet-red-sol.vtk", fens, fes; vectors=[("u", deepcopy(u.values),)]) # fpartitioning, nfpartitions = fine_grid_partitioning(fens, nfpartitions) - nodelists = fine_grid_node_lists(fens, fes, nfpartitions, overlap) + nodelists = CompatibilityModule.fine_grid_node_lists(fens, fes, nfpartitions, overlap) @assert length(nodelists) == nfpartitions partitions = [] diff --git a/src/CompatibilityModule.jl b/src/CompatibilityModule.jl new file mode 100644 index 0000000..383cf1c --- /dev/null +++ b/src/CompatibilityModule.jl @@ -0,0 +1,74 @@ +""" + CompatibilityModule + +Module for compatibility with the old version that went into the review. +""" +module CompatibilityModule + +__precompile__(true) + +using FinEtools +using CoNCMOR +using SparseArrays +using Krylov +using Metis +using LinearOperators +using SparseArrays +using LinearAlgebra +import Base: size, eltype +import LinearAlgebra: mul!, eigen +using Statistics: mean +using ..FENodeToPartitionMapModule: FENodeToPartitionMap +using ShellStructureTopo +using DataDrop + +function _make_overlapping_partition(fes, n2e, overlap, element_1st_partitioning, i) + enl = findall(x -> x == i, element_1st_partitioning) + touched = fill(false, count(fes)) + touched[enl] .= true + for ov in 1:overlap + sfes = subset(fes, enl) + bsfes = meshboundary(sfes) + addenl = Int[] + for e in eachindex(bsfes) + for n in bsfes.conn[e] + for ne in n2e.map[n] + if !touched[ne] + touched[ne] = true + push!(addenl, ne) + end + end + end + end + enl = cat(enl, addenl; dims=1) + end + # vtkexportmesh("i=$i" * "-enl" * "-final" * ".vtk", fens, subset(fes, enl)) + return connectednodes(subset(fes, enl)) +end + +""" + fine_grid_node_lists(fens, fes, npartitions, overlap) + +Make node lists for all partitions of the fine grid. + +The fine grid is partitioned into `npartitions` using the Metis library. The +partitions are extended by the given overlap. Then for each extended partition +of elements, the list of nodes is collected. + +List of these node-lists is returned. +""" +function fine_grid_node_lists(fens, fes, npartitions, overlap) + femm = FEMMBase(IntegDomain(fes, PointRule())) + C = dualconnectionmatrix(femm, fens, nodesperelem(boundaryfe(fes))) + g = Metis.graph(C; check_hermitian=true) + element_1st_partitioning = Metis.partition(g, npartitions; alg=:KWAY) + npartitions = maximum(element_1st_partitioning) + n2e = FENodeToFEMap(fes, count(fens)) + nodelists = [] + for i in 1:npartitions + push!(nodelists, _make_overlapping_partition(fes, n2e, overlap, element_1st_partitioning, i)) + end + nodelists +end + +end # module CompatibilityModule diff --git a/src/FinEtoolsDDMethods.jl b/src/FinEtoolsDDMethods.jl index 2ba9d9d..8d323be 100644 --- a/src/FinEtoolsDDMethods.jl +++ b/src/FinEtoolsDDMethods.jl @@ -26,6 +26,8 @@ using .PartitionCoNCDDMPIModule: partition_multiply!, precondition_global_solve! export CoNCPartitioningInfo, CoNCPartitionData export partition_multiply!, precondition_global_solve!, precondition_local_solve! +include("CompatibilityModule.jl") + include("cg.jl") end # module FinEtoolsDDMethods diff --git a/src/PartitionCoNCDDMPIModule.jl b/src/PartitionCoNCDDMPIModule.jl index e12f203..101d6b8 100644 --- a/src/PartitionCoNCDDMPIModule.jl +++ b/src/PartitionCoNCDDMPIModule.jl @@ -145,7 +145,7 @@ struct CoNCPartitionData{T, IT, FACTOR} nonoverlapping_K::SparseMatrixCSC{T, IT} reduced_K::SparseMatrixCSC{T, IT} overlapping_K_factor::FACTOR - data_rhs::Vector{T} + rhs::Vector{T} ndof::Vector{IT} ntempq::Vector{T} ntempp::Vector{T} @@ -154,7 +154,10 @@ struct CoNCPartitionData{T, IT, FACTOR} otempp::Vector{T} end -function CoNCPartitionData(cpi::CPI, i, fes, Phi, make_matrix) where {CPI<:CoNCPartitioningInfo} +function CoNCPartitionData(cpi::CPI, i, fes, Phi, + make_matrix, + make_interior_load = nothing + ) where {CPI<:CoNCPartitioningInfo} element_lists = cpi.element_lists dof_lists = cpi.dof_lists fr = dofrange(cpi.u, DOF_KIND_FREE) @@ -168,7 +171,11 @@ function CoNCPartitionData(cpi::CPI, i, fes, Phi, make_matrix) where {CPI<:CoNCP Kn_fd = Kn[fr, dr] # Compute the right hand side contribution u_d = gathersysvec(cpi.u, DOF_KIND_DATA) - data_rhs = - Kn_fd * u_d + rhs = - Kn_fd * u_d + if make_interior_load !== nothing + rhs .+= make_interior_load(subset(fes, el))[fr] + end + # Compute contribution to the reduced matrix Kr_ff = Phi' * Kn_ff * Phi # Compute the matrix for the remaining (overlapping - nonoverlapping) elements el = setdiff(element_lists[i].all_connected, element_lists[i].nonoverlapping) @@ -185,7 +192,7 @@ function CoNCPartitionData(cpi::CPI, i, fes, Phi, make_matrix) where {CPI<:CoNCP otempp = zeros(eltype(cpi.u.values), length(odof)) ntempq = zeros(eltype(cpi.u.values), length(ndof)) ntempp = zeros(eltype(cpi.u.values), length(ndof)) - return CoNCPartitionData(Kn_ff, Kr_ff, lu(Ko), data_rhs, ndof, ntempq, ntempp, odof, otempq, otempp) + return CoNCPartitionData(Kn_ff, Kr_ff, lu(Ko), rhs, ndof, ntempq, ntempp, odof, otempq, otempp) end function partition_multiply!(q, partition, p)