From a9cc0cef2836f9fd2adc258d5ad39cda66a01bd1 Mon Sep 17 00:00:00 2001 From: "Craig M. Hamel" Date: Sat, 6 Jun 2026 15:01:50 -0400 Subject: [PATCH] got residual working with parrays. --- examples/parrays-example/script_v4.jl | 88 +++++-------------------- ext/PartitionedArraysExt.jl | 13 ++++ src/assemblers/Assemblers.jl | 12 ++-- src/assemblers/SparseMatrixAssembler.jl | 6 ++ src/assemblers/Vector.jl | 21 ++++-- test/ext/TestPartitionedArraysExt.jl | 66 ++++--------------- 6 files changed, 68 insertions(+), 138 deletions(-) diff --git a/examples/parrays-example/script_v4.jl b/examples/parrays-example/script_v4.jl index 4cf2cd6b..6fadd6ae 100644 --- a/examples/parrays-example/script_v4.jl +++ b/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 diff --git a/ext/PartitionedArraysExt.jl b/ext/PartitionedArraysExt.jl index d6bae7d2..47872d3f 100644 --- a/ext/PartitionedArraysExt.jl +++ b/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 diff --git a/src/assemblers/Assemblers.jl b/src/assemblers/Assemblers.jl index 4c35bc2b..a52dc6b2 100644 --- a/src/assemblers/Assemblers.jl +++ b/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 """ diff --git a/src/assemblers/SparseMatrixAssembler.jl b/src/assemblers/SparseMatrixAssembler.jl index add3280d..483b8744 100644 --- a/src/assemblers/SparseMatrixAssembler.jl +++ b/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 diff --git a/src/assemblers/Vector.jl b/src/assemblers/Vector.jl index c2a137cc..fb0fc682 100644 --- a/src/assemblers/Vector.jl +++ b/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!( diff --git a/test/ext/TestPartitionedArraysExt.jl b/test/ext/TestPartitionedArraysExt.jl index 867d450b..b33baaba 100644 --- a/test/ext/TestPartitionedArraysExt.jl +++ b/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