Merge pull request #19 from PetrKryslUCSD/development

Development

Project.toml

@@ -1,7 +1,7 @@
 name = "FinEtoolsDDMethods"
 uuid = "093fdad3-4fce-4574-aa49-c2a24de4f00d"
 authors = ["Petr Krysl <[email protected]>"]
-version = "0.3.2"
+version = "0.4.0"
 
 [deps]
 CoNCMOR = "3a989f29-cff0-4f1c-87be-5dcdd41bd240"

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.
+to almost incompressible elasticity and general shells.
+

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.
+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.

examples/conc/heat/Poisson2D_overlapped_mpi_examples.jl

@@ -69,57 +69,66 @@ 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)
         femm2 = FEMMHeatDiff(IntegDomain(fes, volrule, 1.0), material)
         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)

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 = []

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 = []

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

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)")

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 = []

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

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