got residual working with parrays.

examples/parrays-example/script_v4.jl

@@ -6,7 +6,7 @@ using PartitionedArrays
 using StaticArrays
 
 include("../../test/poisson/TestPoissonCommon.jl")
-f(X, _) = 2. * π^2 * cos(π * X[1]) * cos(π * X[2])
+f(X, _) = 2. * π^2 * sin(2π * X[1]) * sin(2π * X[2])
 bc_func(_, _) = 0.
 mesh_file = Base.source_dir() * "/square.g"
 output_file = Base.source_dir() * "/output.e"
@@ -29,97 +29,39 @@ dof     = DofManager(u, dbcs, num_ranks, ranks, mesh_file, mesh)
 asm     = SparseMatrixAssembler(dof)
 p       = create_parameters(mesh, asm, physics, props; dirichlet_bcs = dbcs)
 
-# stuff to still work out below
-mat_pattern_new = FiniteElementContainers.create_matrix_sparsity_pattern(dof)
-vec_pattern_new = FiniteElementContainers.create_vector_sparsity_pattern(dof)
-X = nodal_coordinates(dof)
-
-u_analytic(x) = 0.5 * (x[1] + x[2])
-h = 0.1
-# Ae = (h^2 / 6) * [
-#      4.0  -1.0  -1.0  -2.0
-#     -1.0   4.0  -2.0  -1.0
-#     -1.0  -2.0   4.0  -1.0
-#     -2.0  -1.0  -1.0   4.0
-# ]
-# Ae = SMatrix{4, 4, Float64, 16}([
-#     0.666667  -0.166667  -0.333333  -0.166667
-#     -0.166667   0.666667  -0.166667  -0.333333
-#     -0.333333  -0.166667   0.666667  -0.166667
-#     -0.166667  -0.333333  -0.166667   0.666667
-# ])
-# VVs_new = map(dof.var.fspace) do fspace
-#     VVs = Float64[]
-#     Ae_vec = vec(Ae)
-#     for b in 1:FiniteElementContainers.num_blocks(fspace)
-#         conn = connectivity(fspace, b)
-#         for e in axes(conn, 2)
-#             for i in axes(conn, 1)
-#                 for j in axes(conn, 1)
-#                     push!(VVs, Ae[i, j])
-#                 end
-#             end
-#         end
-#     end
-#     VVs
-# end
-
-Xs = nodal_coordinates(dof)
-Vs_new = map(dof.var.fspace, dof.local_dof_managers, dof.field_partition.parts, partition(Xs)) do fspace, local_dof, field_part, X
-    Vs = Float64[]
-    ue = zeros(size(Ae, 2))
-    ge = similar(ue)
-    field_ltg = local_to_global(field_part)
-    for b in 1:FiniteElementContainers.num_blocks(fspace)
-        block_local_conns = connectivity(fspace, b)
-        block_global_conns = field_ltg[block_local_conns]
-        for e in axes(block_global_conns, 2)
-            fill!(ue, 0.0)
-            for i in axes(block_global_conns, 1)
-                # this is for bcs here
-                # if insorted(block_global_conns[i, e], dof.dirichlet_dofs)
-                if insorted(block_local_conns[i, e], local_dof.dirichlet_dofs)
-                    ue[i] = u_analytic(fspace.coords[:, block_local_conns[i, e]])
-                end
-            end
-
-            mul!(ge, Ae, ue)
-
-            for i in axes(block_global_conns, 1)
-                push!(Vs, -ge[i])
-            end
-        end
-    end
-    Vs
-end
+# TODO
+# first setup a solution vector by doing IterativeSolvers.cg(K, -0)
+# so that we get a solution vector with the right sizes
+# then we can copy that and cache them as a 
+# current solution and increment
 
+# clean up below
 Uu = create_unknowns(asm)
+
 U = create_field(asm)
 update_field_unknowns!(U, dof, Uu)
+
 assemble_stiffness!(asm, stiffness, U, p)
 K = stiffness(asm)
 
-# A_new = psparse(mat_pattern_new, VVs_new) |> fetch
-b_new = pvector(vec_pattern_new, Vs_new) |> fetch
-x_new = IterativeSolvers.cg(K, b_new, verbose = i_am_main(rank))
+assemble_vector!(asm, residual, U, p)
+R = residual(asm)
+
+x_new = IterativeSolvers.cg(K, -R, verbose = i_am_main(rank))
 
 U = create_field(dof)
 # update_field_dirichlet_bcs!(U, dof)
 update_field_unknowns!(U, dof, x_new)
 
-map(partition(U), partition(Xs), dof.var.fspace, mesh, ranks) do U_local, X, fspace, mesh_local, rank
+map(partition(U), dof.var.fspace, mesh, ranks) do U_local, fspace, mesh_local, rank
     u_temp = ScalarFunction(fspace, "u")
-    x_temp = VectorFunction(fspace, "x")
-    mesh_file = mesh_local.mesh_obj.file_name
     output_file_temp = output_file * ".$(num_ranks).$(rank - 1)"
-    pp = PostProcessor(mesh_local, output_file_temp, u_temp, x_temp)
+    pp = PostProcessor(mesh_local, output_file_temp, u_temp)
     write_times(pp, 1, 0.0)
     write_field(pp, 1, ["u"], U_local)
-    write_field(pp, 1, ["x_x", "x_y"], X)
     close(pp)
 end
 
 map_main(ranks) do rank
-    # output_file = splitext(mesh_file)[1] * ".e.$(num_ranks)."
     epu(output_file * ".$(num_ranks).$(rank - 1)")
 end

ext/PartitionedArraysExt.jl

@@ -507,6 +507,12 @@ function FiniteElementContainers.assemble_stiffness!(asm::PSparseMatrixAssembler
 	end
 end
 
+function FiniteElementContainers.assemble_vector!(asm::PSparseMatrixAssembler, func, u, p)
+	map(asm.local_assemblers, partition(u), p.local_parameters) do local_asm, local_u, local_p
+		assemble_vector!(local_asm, func, local_u, local_p)
+	end
+end
+
 function FiniteElementContainers.create_field(asm::PSparseMatrixAssembler)
 	return create_field(asm.vector_pattern.dof)
 end
@@ -515,6 +521,13 @@ function FiniteElementContainers.create_unknowns(asm::PSparseMatrixAssembler)
 	return create_unknowns(asm.vector_pattern.dof)
 end
 
+function FiniteElementContainers.residual(asm::PSparseMatrixAssembler)
+	vals = map(asm.local_assemblers) do local_asm
+		local_asm.residual_unknowns
+	end
+	return pvector(asm.vector_pattern, vals) |> fetch
+end
+
 function FiniteElementContainers.stiffness(asm::PSparseMatrixAssembler)
 	vals = map(asm.local_assemblers) do local_asm
 		local_asm.stiffness_storage

src/assemblers/Assemblers.jl

@@ -97,7 +97,7 @@ function _assemble_element!(
   end_id = start_id + NDOF - 1
   ids = start_id:end_id
   for (i, id) in enumerate(ids)
-    storage[id] += R_el.data[i]
+    storage[id] = R_el.data[i]
   end
   return nothing
 end
@@ -326,10 +326,10 @@ end
 $(TYPEDSIGNATURES)
 assumes assemble_vector! has already been called
 """
-function residual(asm::AbstractAssembler; use_sparse_vector = false)
-  # if use_sparse_vector
-    # return sparsevec(asm.vector_pattern, asm.residual_unknowns)
-  # else
+function residual(asm::AbstractAssembler)
+  if _use_sparse_vector(asm)
+    return sparsevec(asm.vector_pattern, asm.residual_unknowns)
+  else
     if _is_condensed(asm.dof)
       _adjust_vector_entries_for_constraints!(
         asm.residual_storage, asm.constraint_storage
@@ -343,7 +343,7 @@ function residual(asm::AbstractAssembler; use_sparse_vector = false)
       )
       return asm.residual_unknowns
     end
-  # end
+  end
 end
 
 """

src/assemblers/SparseMatrixAssembler.jl

@@ -256,5 +256,11 @@ function update_dofs!(assembler::AbstractAssembler, dirichlet_bcs::DirichletBCs)
     end
     _update_dofs!(assembler.vector_pattern, assembler.dof, ddofs)
   end
+
+  # problaby a better way to order all this logic for different things
+  if _use_sparse_vector(assembler)
+    _, n_entries = _setup_block_sizes(assembler.dof, 1)
+    resize!(assembler.residual_unknowns, n_entries)
+  end
   return nothing
 end

src/assemblers/Vector.jl

@@ -2,10 +2,16 @@
 $(TYPEDSIGNATURES)
 """
 function assemble_vector!(
-  assembler, func::F, Uu, p
+  assembler::AbstractAssembler, func::F, Uu, p
 ) where F <: Function
+  if _use_sparse_vector(assembler)
+    storage = assembler.residual_unknowns
+  else
+    storage = assembler.residual_storage
+  end
   assemble_vector!(
-    assembler.residual_storage, 
+    # assembler.residual_storage, 
+    storage,
     assembler.vector_pattern, assembler.dof,
     func, Uu, p;
     use_inplace_methods = _use_inplace_methods(assembler),
@@ -34,11 +40,12 @@ function assemble_vector!(
   conns = fspace.elem_conns
   # foreach_block(conns, p.physics, p.properties, fspace.ref_fes) do physics, props, ref_fe, b
   foreach_block(fspace, p) do physics, props, ref_fe, b
-    if use_sparse_vector
-      field = block_view(storage, pattern, b)
-    else
-      field = storage
-    end
+    # if use_sparse_vector
+    #   field = block_view(storage, pattern, b)
+    # else
+    #   field = storage
+    # end
+    field = storage
 
     if use_inplace_methods
       _assemble_block!(

test/ext/TestPartitionedArraysExt.jl

@@ -8,7 +8,7 @@
 
     if !Sys.iswindows()
         include("../poisson/TestPoissonCommon.jl")
-        f(X, _) = 2. * π^2 * cos(π * X[1]) * cos(π * X[2])
+        f(X, _) = 2. * π^2 * sin(2π * X[1]) * sin(2π * X[2])
         bc_func(_, _) = 0.
         mesh_file = Base.source_dir() * "/square.g"
         output_file = Base.source_dir() * "/output.e"
@@ -31,74 +31,36 @@
         asm     = SparseMatrixAssembler(dof)
         p       = create_parameters(mesh, asm, physics, props; dirichlet_bcs = dbcs)
 
-        # stuff to still work out below
-        mat_pattern_new = FiniteElementContainers.create_matrix_sparsity_pattern(dof)
-        vec_pattern_new = FiniteElementContainers.create_vector_sparsity_pattern(dof)
-        X = nodal_coordinates(dof)
-
-        u_analytic(x) = 0.5 * (x[1] + x[2])
-
-        # need this for residual for now
-        Ae = SMatrix{4, 4, Float64, 16}([
-            0.666667  -0.166667  -0.333333  -0.166667
-            -0.166667   0.666667  -0.166667  -0.333333
-            -0.333333  -0.166667   0.666667  -0.166667
-            -0.166667  -0.333333  -0.166667   0.666667
-        ])
-
-        Xs = nodal_coordinates(dof)
-        Vs_new = map(dof.var.fspace, dof.local_dof_managers, dof.field_partition.parts, partition(Xs)) do fspace, local_dof, field_part, X
-            Vs = Float64[]
-            ue = zeros(size(Ae, 2))
-            ge = similar(ue)
-            field_ltg = local_to_global(field_part)
-            for b in 1:FiniteElementContainers.num_blocks(fspace)
-                block_local_conns = connectivity(fspace, b)
-                block_global_conns = field_ltg[block_local_conns]
-                for e in axes(block_global_conns, 2)
-                    fill!(ue, 0.0)
-                    for i in axes(block_global_conns, 1)
-                        # this is for bcs here
-                        # if insorted(block_global_conns[i, e], dof.dirichlet_dofs)
-                        if insorted(block_local_conns[i, e], local_dof.dirichlet_dofs)
-                            ue[i] = u_analytic(fspace.coords[:, block_local_conns[i, e]])
-                        end
-                    end
-
-                    mul!(ge, Ae, ue)
-
-                    for i in axes(block_global_conns, 1)
-                        push!(Vs, -ge[i])
-                    end
-                end
-            end
-            Vs
-        end
+        # TODO
+        # first setup a solution vector by doing IterativeSolvers.cg(K, -0)
+        # so that we get a solution vector with the right sizes
+        # then we can copy that and cache them as a 
+        # current solution and increment
 
+        # clean up below
         Uu = create_unknowns(asm)
+
         U = create_field(asm)
         update_field_unknowns!(U, dof, Uu)
 
-        # stiffness matrix now handled by FEC
         assemble_stiffness!(asm, stiffness, U, p)
         K = stiffness(asm)
 
-        b_new = pvector(vec_pattern_new, Vs_new) |> fetch
-        x_new = IterativeSolvers.cg(K, b_new, verbose = i_am_main(rank))
+        assemble_vector!(asm, residual, U, p)
+        R = residual(asm)
+
+        x_new = IterativeSolvers.cg(K, -R, verbose = i_am_main(rank))
 
         U = create_field(dof)
         # update_field_dirichlet_bcs!(U, dof)
         update_field_unknowns!(U, dof, x_new)
 
-        map(partition(U), partition(Xs), dof.var.fspace, mesh, ranks) do U_local, X, fspace, mesh_local, rank
+        map(partition(U), dof.var.fspace, mesh, ranks) do U_local, fspace, mesh_local, rank
             u_temp = ScalarFunction(fspace, "u")
-            x_temp = VectorFunction(fspace, "x")
-            mesh_file = mesh_local.mesh_obj.file_name
             output_file_temp = output_file * ".$(num_ranks).$(rank - 1)"
-            pp = PostProcessor(mesh_local, output_file_temp, u_temp, x_temp)
+            pp = PostProcessor(mesh_local, output_file_temp, u_temp)
             write_times(pp, 1, 0.0)
             write_field(pp, 1, ["u"], U_local)
-            write_field(pp, 1, ["x_x", "x_y"], X)
             close(pp)
         end