From b262aa049d6ac2052c7183cdc3e747435f960920 Mon Sep 17 00:00:00 2001 From: "Craig M. Hamel" Date: Mon, 24 Aug 2026 22:23:51 -0400 Subject: [PATCH] Another in place scatter method. --- examples/cooks_membrane/script.jl | 72 +++++++++++++---------- examples/mpi-example/example.jl | 48 --------------- examples/mpi-example/square.g | Bin 8392 -> 0 bytes examples/poisson/script.jl | 2 +- src/Formulations.jl | 93 ++++++++++++++++++++++++++++++ test/poisson/TestPoissonCommon.jl | 3 +- 6 files changed, 138 insertions(+), 80 deletions(-) delete mode 100644 examples/mpi-example/example.jl delete mode 100644 examples/mpi-example/square.g diff --git a/examples/cooks_membrane/script.jl b/examples/cooks_membrane/script.jl index dc70f67a..f8fcb672 100644 --- a/examples/cooks_membrane/script.jl +++ b/examples/cooks_membrane/script.jl @@ -165,8 +165,10 @@ FEC.update_dofs!( (p_u.periodic_bcs, p_p.periodic_bcs) ) -pp_u = PostProcessor(mesh_u, "u.exo", u) -pp_p = PostProcessor(mesh_p, "p.exo", p) +# pp_u = PostProcessor(mesh_u, "u.exo", u) +# pp_p = PostProcessor(mesh_p, "p.exo", p) + + # solver = NewtonSolver(DirectLinearSolver(asm)) # integrator = QuasiStaticIntegrator(solver) @@ -186,39 +188,49 @@ Uu = create_unknowns(asm) # temp = K_up .- K_pu' # display(K_up .- K_pu') -for n in 1:10 - FiniteElementContainers.update_time!(params) - FiniteElementContainers.update_bc_values!(params, asm) - r0 = -1e6 - for iter in 1:10 - assemble_vector!(asm, residual, Uu, params) - R = residual(asm) +function test_func(asm, u, p) + FiniteElementContainers.update_time!(p) + FiniteElementContainers.update_bc_values!(p, asm) + assemble_vector!(asm, residual, u, p) + return nothing +end + +@time test_func(asm, Uu, params) +@time test_func(asm, Uu, params) + +# for n in 1:10 +# FiniteElementContainers.update_time!(params) +# FiniteElementContainers.update_bc_values!(params, asm) +# r0 = -1e6 +# for iter in 1:10 +# assemble_vector!(asm, residual, Uu, params) +# R = residual(asm) - rnorm = norm(R) - if iter == 1 - r0 = rnorm - end - - if rnorm / r0 < 1e-8 - break - end +# rnorm = norm(R) +# if iter == 1 +# r0 = rnorm +# end + +# if rnorm / r0 < 1e-8 +# break +# end - assemble_stiffness!(asm, stiffness, Uu, params) - K = stiffness(asm) +# assemble_stiffness!(asm, stiffness, Uu, params) +# K = stiffness(asm) - ΔU = K \ R - Uu .-= ΔU +# ΔU = K \ R +# Uu .-= ΔU - println("iter = $iter, |R| = $(rnorm / r0), |ΔU| = $(norm(ΔU))") +# println("iter = $iter, |R| = $(rnorm / r0), |ΔU| = $(norm(ΔU))") - end +# end - write_times(pp_u, n + 1, params[1].times.time_current) - write_times(pp_p, n + 1, params[2].times.time_current) - write_field(pp_u, n + 1, ("displ_x", "displ_y"), params[1].field) - write_field(pp_p, n + 1, ("pressure",), params[2].field) -end -close(pp_u) -close(pp_p) +# write_times(pp_u, n + 1, params[1].times.time_current) +# write_times(pp_p, n + 1, params[2].times.time_current) +# write_field(pp_u, n + 1, ("displ_x", "displ_y"), params[1].field) +# write_field(pp_p, n + 1, ("pressure",), params[2].field) +# end +# close(pp_u) +# close(pp_p) # Δu = K \ R diff --git a/examples/mpi-example/example.jl b/examples/mpi-example/example.jl deleted file mode 100644 index 3a6cfc1b..00000000 --- a/examples/mpi-example/example.jl +++ /dev/null @@ -1,48 +0,0 @@ -import FiniteElementContainers: - communication_graph, - cpu, - decompose_mesh, - distribute, - global_colorings, - DistributedDevice, - getdata, - num_dofs_per_rank, - ParVector, - rank_devices, - scatter_ghosts! -using Exodus -using KernelAbstractions -using MPI - -backend = CPU() -num_dofs = 1 -num_ranks = 4 -ranks = 1:4 |> collect -mesh_file = Base.source_dir() * "/square.g" - -decompose_mesh(mesh_file, num_ranks) -global_dofs_to_colors = global_colorings(mesh_file, num_dofs, num_ranks) - -comm = MPI.COMM_WORLD -# comm = nothing - -ranks = distribute(ranks, comm) -# n_dofs_per_rank = distribute(num_dofs_per_rank(global_dofs_to_colors), comm) - -# shards = shard_indices(global_dofs_to_colors, ) -comm_graphs = map(ranks) do rank - mesh_file = Base.source_dir() * "/square.g" - mesh_file = mesh_file * ".$num_ranks" * ".$(lpad(rank - 1, Exodus.exodus_pad(num_ranks |> Int32), '0'))" - comm_graph = communication_graph(mesh_file, global_dofs_to_colors, num_dofs, rank) -end - -par_vec = ParVector(comm_graphs) -parts = map(par_vec.parts, ranks) do part, rank - part .= rank - part -end -par_vec = ParVector(parts) -par_vec_2 = scatter_ghosts!(par_vec) -MPI.Barrier(comm) -# par_vec_2 = par_vec_2.parts -print("part on $(MPI.Comm_rank(comm) + 1) is $par_vec_2\n") diff --git a/examples/mpi-example/square.g b/examples/mpi-example/square.g deleted file mode 100644 index f41996e896f4eddf675e0d719918e7fa845726b7..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 8392 zcmeI0%X1Xh8O0?)LL*B&jBUUe7#_C44~>`^fh9jiAi&0NCxJN6W;AL^4KvfjOb-bE z1zDz2u9Q%vQdwk~e<1%sT!p={H{RLhR4S>;!iy}j$Rg+0{oRq&Fd2wdNnGJp9qzfO zZ-4jP*C%>r&Ys`Ymy}^uqryEH_CC_iF8bb4{KAq--;iT=&R{QQvqeaOYCx8gX*q>Tt$c~))&-gUp( z3c|XYeFH-O>AnxXRK2x}zGuD~mLrYeHNv1Cc{jr5oYx9I&^pE*nAebynGdR3&*Z_H zXFT6kD=!5#)>7;5$@D_BD7o`~R(wRkPZl+sa@1zYm>=&_2BILUGN20ka<0;z4WeVY z#jxh*D$R1Rkei7X!+H?5a_5%9xpr&pN;TKI)-E^wv4!JP+40F-VLX=~&reKc^M$F2 z@u}iu98GJYdZ#cPpNAWZp$){?ZhX293d~a_`v?7#@n=`fs!dr8qnuH^(gbM;$h^1}cM$oAl*&B=hOgH*;IhSqG-e9G=&kR+cJGN#T$TOe+ON2OoS5*! zTDZdV-z1F#e^xy@7qlX8zFdhs%VdfFweM)Z`=^t44*pCvF*_)IojfW1^)}_>PVMx+ z+@h?MT>Yo{xRhM|RUs}VSASvgT>W|c{&Fd)Z((K?Z|U{@cjHc^&RKeWyT$Fc4sLzjI^21F-FsH_D zy_1h! zzVA~1*7V<{o_}cVblu4lE^j`lH)aK9!*+eH+)VzLyfGoaRZ*6{Rc<%l?Eme?uP!D3 zH|JZ3Zv_6o59j-j?>02J`%XfGyKfx6o#1@y@D1R;d-!c48hn@V zh-vUgg~v?86FznedyL~_uW4|6JZ>5sANx#$h!~Y8o7GubBqN+Xd6$_;}qk__XkbfCk6Mn*tgfZ*K``c*5J; z!aK(C_O5AgyuD`{9Bh#VM4P#I<^mhgwdm*F;uzPbTNyu@ z%b$vW1@6)#KR$J0T)z%?d91_#%{_E}iS9Cpex5E&OchBJcQ5}T4=+sR|Jxqc{!D+@ zwdj|+N!ToG5e9@oVXH7CWQ1YiAz?%q6}Ac6g&o39;bCEy@QCoJ@R+b$*dy!}9vAir zPYC;k1HwTeD;yFI3rB>b!jr;N!qdVr;kfXO@T~Bh@Vqc4pXME zF9}D4Ugw%zz}Y@5To%&j+m|*kEC}R`Md6ST2v>xwLRF{i|C?;JUZ{6Zcg*N}6_J>(#A5jlz6v|AVvxO2S2IK0muvA^so`^et0XRI%4 zfp6x)y?#ev%_jx+oi!ph=Y;W@8|UMca7JLw-xn$Z^TRXVnFDiSPRxxtGT*TC+x6FK LsXS=k-#Yvg7G2h$ diff --git a/examples/poisson/script.jl b/examples/poisson/script.jl index 1ed8f848..07aba849 100644 --- a/examples/poisson/script.jl +++ b/examples/poisson/script.jl @@ -35,7 +35,7 @@ physics = Poisson(f) props = create_properties(physics) u = ScalarFunction(V, "u") # asm = SparseMatrixAssembler(u; use_condensed = true) -asm = SparseMatrixAssembler(u) +asm = SparseMatrixAssembler(u; use_inplace_methods = true) # dbcs = nothing diff --git a/src/Formulations.jl b/src/Formulations.jl index e47a0832..b84e510a 100644 --- a/src/Formulations.jl +++ b/src/Formulations.jl @@ -191,6 +191,99 @@ function scatter_with_gradients_and_gradients!( return nothing end +""" +General stiffness like matrix scatter method +""" +function scatter_with_gradients_and_gradients!( + storage::AbstractVector, + form::AbstractElementFormulation{ND, NF}, + e, + conns, + ∇N_X, + K::AbstractMatrix{T}, +) where {ND, NF, T <: Number} + + @assert size(K, 1) == NF^2 + @assert size(K, 2) == NF^2 + @assert ND == size(∇N_X, 2) + + N = size(∇N_X, 1) + NDPE = N * NF + + start_id = (e - 1) * NDPE^2 + 1 + ids = start_id:(start_id + NDPE^2 - 1) + + inc = 1 + + for n2 in 1:N + for d2 in 1:NF + for n1 in 1:N + for d1 in 1:NF + contrib = zero(T) + + for J in 1:ND + a = (J - 1) * ND + d1 + + for L in 1:ND + b = (L - 1) * ND + d2 + contrib += ∇N_X[n1, J] * K[a, b] * ∇N_X[n2, L] + end + end + + storage[ids[inc]] += contrib + inc += 1 + end + end + end + end + + return nothing +end + +function scatter_gradient_gradient!( + storage::AbstractField, + form::AbstractElementFormulation{ND, NF}, + e, + conns, + ∇N_X, + A::AbstractMatrix{T}, + v_el::AbstractVector{T}, +) where {ND, NF, T <: Number} + + @assert ND == size(∇N_X, 2) + @assert size(A) == (ND * NF, ND * NF) + + N = size(∇N_X, 1) + NDPE = N * NF + + @assert length(v_el) == NDPE + + for n1 in 1:N + for d1 in 1:NF + global_id = NF * (conns[n1] - 1) + d1 + contrib = zero(T) + + for n2 in 1:N + for d2 in 1:NF + local_id = NF * (n2 - 1) + d2 + for j1 in 1:ND + a = (j1 - 1) * NF + d1 + + for j2 in 1:ND + b = (j2 - 1) * NF + d2 + contrib += ∇N_X[n1, j1] * A[a, b] * ∇N_X[n2, j2] * v_el[local_id] + end + end + end + end + + fec_atomic_add!(storage, global_id, contrib) + end + end + + return nothing +end + # implement for those that have it """ Scalar equation specialization diff --git a/test/poisson/TestPoissonCommon.jl b/test/poisson/TestPoissonCommon.jl index 69b07d55..b48df482 100644 --- a/test/poisson/TestPoissonCommon.jl +++ b/test/poisson/TestPoissonCommon.jl @@ -1,6 +1,8 @@ # since this type of physics has no state # just return the old state which will be empty +using ReferenceFiniteElements + struct Poisson{F <: Function} <: AbstractPhysics{1, 0, 0} func::F end @@ -91,7 +93,6 @@ end interps = map_interpolants(interps, x_el) (; X_q, N, ∇N_X, JxW) = interps ∇u_q = interpolate_field_gradients(physics, interps, u_el) - # R_q = ∇u_q * ∇N_X' - N' * physics.func(X_q, 0.0) form = GeneralFormulation{size(X_q, 1), num_fields(physics)}() scatter_with_gradients!(storage, form, e, conn, ∇N_X, JxW * ∇u_q) scatter_with_values!(storage, form, e, conn, N, -JxW * physics.func(X_q, 0.0))