From 498577e68f50a5e2c68828906a968fcc08778da3 Mon Sep 17 00:00:00 2001 From: Chris Richardson Date: Thu, 9 Oct 2025 16:03:56 +0100 Subject: [PATCH 1/3] Updates for latest main --- src/cg.h | 2 +- src/cgpoisson_problem.cpp | 60 ++++++++++++++++++-------------------- src/elasticity_problem.cpp | 19 ++++++------ src/poisson_problem.cpp | 16 +++++----- 4 files changed, 46 insertions(+), 51 deletions(-) diff --git a/src/cg.h b/src/cg.h index 78af592..018090e 100644 --- a/src/cg.h +++ b/src/cg.h @@ -20,7 +20,7 @@ void axpy(la::Vector& r, U alpha, const la::Vector& x, const la::Vector& y) { std::transform(x.array().begin(), x.array().end(), y.array().begin(), - r.mutable_array().begin(), + r.array().begin(), [alpha](auto x, auto y) { return alpha * x + y; }); } diff --git a/src/cgpoisson_problem.cpp b/src/cgpoisson_problem.cpp index 3704e6e..bdf172a 100644 --- a/src/cgpoisson_problem.cpp +++ b/src/cgpoisson_problem.cpp @@ -70,7 +70,7 @@ cgpoisson::problem(std::shared_ptr> mesh, int order, common::Timer t2("ZZZ Create boundary conditions"); // Define boundary condition auto u0 = std::make_shared>(V); - u0->x()->set(0); + std::ranges::fill(u0->x()->array(), 0.0); // Find facets with bc applied const int tdim = mesh->topology()->dim(); @@ -147,24 +147,25 @@ cgpoisson::problem(std::shared_ptr> mesh, int order, // Create la::Vector la::Vector b(L->function_spaces()[0]->dofmap()->index_map, L->function_spaces()[0]->dofmap()->index_map_bs()); - b.set(0); + std::ranges::fill(b.array(), 0.0); + common::Timer t5("ZZZ Assemble vector"); const std::vector constants_L = fem::pack_constants(*L); auto coeffs_L = fem::allocate_coefficient_storage(*L); fem::pack_coefficients(*L, coeffs_L); - fem::assemble_vector(b.mutable_array(), *L, constants_L, - fem::make_coefficients_span(coeffs_L)); + fem::assemble_vector(b.array(), *L, std::span(constants_L), + fem::make_coefficients_span(coeffs_L)); // Apply lifting to account for Dirichlet boundary condition // b <- b - A * x_bc - bc->set(un->x()->mutable_array(), std::nullopt, -1.0); - fem::assemble_vector(b.mutable_array(), *M); + bc->set(un->x()->array(), std::nullopt, -1.0); + fem::assemble_vector(b.array(), *M); // Communicate ghost values b.scatter_rev(std::plus()); // Set BC dofs to zero (effectively zeroes columns of A) - bc->set(b.mutable_array(), std::nullopt, 0.0); + bc->set(b.array(), std::nullopt, 0.0); b.scatter_fwd(); // Pack coefficients and constants @@ -185,52 +186,47 @@ cgpoisson::problem(std::shared_ptr> mesh, int order, int bs = V->dofmap()->bs(); common::Scatterer sct(*idx_map, bs); - std::vector local_buffer(sct.local_buffer_size(), 0); - std::vector remote_buffer(sct.remote_buffer_size(), 0); - - common::Scatterer<>::type type; - if (scatterer == "neighbor") - type = common::Scatterer<>::type::neighbor; - if (scatterer == "p2p") - type = common::Scatterer<>::type::p2p; - - std::vector request = sct.create_request_vector(type); + std::vector local_buffer(sct.local_indices().size(), 0); + std::vector remote_buffer(sct.remote_indices().size(), 0); // Create function for computing the action of A on x (y = Ax) auto action = [&](la::Vector& x, la::Vector& y) { // Zero y - y.set(0.0); + std::ranges::fill(y.array(), 0.0); // Update coefficient un (just copy data from x to un) - std::copy(x.array().begin(), x.array().end(), - un->x()->mutable_array().begin()); + std::copy(x.array().begin(), x.array().end(), un->x()->array().begin()); // Compute action of A on x fem::pack_coefficients(*M, coeff); - fem::assemble_vector(y.mutable_array(), *M, std::span(constants), + fem::assemble_vector(y.array(), *M, std::span(constants), fem::make_coefficients_span(coeff)); // Set BC dofs to zero (effectively zeroes rows of A) - bc->set(y.mutable_array(), std::nullopt, 0.0); + bc->set(y.array(), std::nullopt, 0.0); // Accumuate ghost values // y.scatter_rev(std::plus()); const std::int32_t local_size = bs * idx_map->size_local(); const std::int32_t num_ghosts = bs * idx_map->num_ghosts(); - std::span remote_data(y.mutable_array().data() + local_size, - num_ghosts); - std::span local_data(y.mutable_array().data(), local_size); - sct.scatter_rev_begin(remote_data, remote_buffer, local_buffer, - pack_fn, request, type); - sct.scatter_rev_end(local_buffer, local_data, unpack_fn, - std::plus(), request); + std::span remote_data(y.array().data() + local_size, num_ghosts); + std::span local_data(y.array().data(), local_size); + + MPI_Request request = MPI_REQUEST_NULL; + pack_fn(remote_data, sct.remote_indices(), remote_buffer); + sct.scatter_rev_begin(remote_buffer.data(), local_buffer.data(), request); + sct.scatter_end(request); + unpack_fn(local_buffer, sct.local_indices(), local_data, std::plus()); // Update ghost values - sct.scatter_fwd_begin(local_data, local_buffer, remote_buffer, pack_fn, - request, type); - sct.scatter_fwd_end(remote_buffer, remote_data, unpack_fn, request); + request = MPI_REQUEST_NULL; + pack_fn(local_data, sct.local_indices(), local_buffer); + sct.scatter_fwd_begin(local_buffer.data(), remote_buffer.data(), request); + sct.scatter_end(request); + unpack_fn(remote_buffer, sct.remote_indices(), remote_data, + [](auto x, auto y) { return y; }); }; common::Timer tcg; diff --git a/src/elasticity_problem.cpp b/src/elasticity_problem.cpp index 52c409a..6faa04f 100644 --- a/src/elasticity_problem.cpp +++ b/src/elasticity_problem.cpp @@ -44,15 +44,15 @@ MatNullSpace build_near_nullspace(const fem::FunctionSpace& V) std::int32_t length_block = map->size_local() + map->num_ghosts(); for (int k = 0; k < 3; ++k) { - std::span x = basis[k].mutable_array(); + std::span x = basis[k].array(); for (std::int32_t i = 0; i < length_block; ++i) x[bs * i + k] = 1.0; } // Rotations - auto x3 = basis[3].mutable_array(); - auto x4 = basis[4].mutable_array(); - auto x5 = basis[5].mutable_array(); + auto x3 = basis[3].array(); + auto x4 = basis[4].array(); + auto x5 = basis[5].array(); const std::vector x = V.tabulate_dof_coordinates(false); const std::int32_t* dofs = V.dofmap()->map().data_handle(); @@ -220,13 +220,12 @@ elastic::problem(std::shared_ptr> mesh, int order) const std::vector constants_L = fem::pack_constants(*L); auto coeffs_L = fem::allocate_coefficient_storage(*L); fem::pack_coefficients(*L, coeffs_L); - fem::assemble_vector(b.mutable_array(), *L, constants_L, - fem::make_coefficients_span(coeffs_L)); - fem::apply_lifting(b.mutable_array(), {*a}, {constants_L}, - {fem::make_coefficients_span(coeffs_L)}, - {{*bc}}, {}, 1.0); + fem::assemble_vector(b.array(), *L, std::span(constants_L), + fem::make_coefficients_span(coeffs_L)); + fem::apply_lifting(b.array(), {*a}, {constants_L}, + {fem::make_coefficients_span(coeffs_L)}, {{*bc}}, {}, 1.0); b.scatter_rev(std::plus<>()); - bc->set(b.mutable_array(), std::nullopt); + bc->set(b.array(), std::nullopt); t3.stop(); t3.flush(); diff --git a/src/poisson_problem.cpp b/src/poisson_problem.cpp index 39cc21b..c06f688 100644 --- a/src/poisson_problem.cpp +++ b/src/poisson_problem.cpp @@ -51,7 +51,7 @@ poisson::problem(std::shared_ptr> mesh, int order) common::Timer t2("ZZZ Create boundary conditions"); // Define boundary condition auto u0 = std::make_shared>(V); - u0->x()->set(0); + std::ranges::fill(u0->x()->array(), 0.0); // Find facets with bc applied const int tdim = mesh->topology()->dim(); @@ -141,18 +141,18 @@ poisson::problem(std::shared_ptr> mesh, int order) // Create la::Vector la::Vector b(L->function_spaces()[0]->dofmap()->index_map, L->function_spaces()[0]->dofmap()->index_map_bs()); - b.set(0); + std::ranges::fill(b.array(), 0.0); + common::Timer t5("ZZZ Assemble vector"); const std::vector constants_L = fem::pack_constants(*L); auto coeffs_L = fem::allocate_coefficient_storage(*L); fem::pack_coefficients(*L, coeffs_L); - fem::assemble_vector(b.mutable_array(), *L, constants_L, - fem::make_coefficients_span(coeffs_L)); - fem::apply_lifting(b.mutable_array(), {*a}, {constants_L}, - {fem::make_coefficients_span(coeffs_L)}, - {{*bc}}, {}, 1.0); + fem::assemble_vector(b.array(), *L, std::span(constants_L), + fem::make_coefficients_span(coeffs_L)); + fem::apply_lifting(b.array(), {*a}, {constants_L}, + {fem::make_coefficients_span(coeffs_L)}, {{*bc}}, {}, 1.0); b.scatter_rev(std::plus<>()); - bc->set(b.mutable_array(), std::nullopt); + bc->set(b.array(), std::nullopt); t5.stop(); t5.flush(); From ae5fcb01d93f3b669a68f03195aff204f575e586 Mon Sep 17 00:00:00 2001 From: Chris Richardson Date: Thu, 9 Oct 2025 16:06:30 +0100 Subject: [PATCH 2/3] Remove deprecated lines --- src/elasticity_problem.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/elasticity_problem.cpp b/src/elasticity_problem.cpp index 6faa04f..b9059e2 100644 --- a/src/elasticity_problem.cpp +++ b/src/elasticity_problem.cpp @@ -117,7 +117,7 @@ elastic::problem(std::shared_ptr> mesh, int order) // Define boundary condition auto u0 = std::make_shared>(V); - u0->x()->set(0); + std::ranges::fill(u0->x()->array(), 0.0); const int tdim = mesh->topology()->dim(); @@ -215,7 +215,8 @@ elastic::problem(std::shared_ptr> mesh, int order) // Wrap la::Vector with Petsc Vec la::Vector b(L->function_spaces()[0]->dofmap()->index_map, L->function_spaces()[0]->dofmap()->index_map_bs()); - b.set(0); + std::ranges::fill(b.array(), 0.0); + common::Timer t3("ZZZ Assemble vector"); const std::vector constants_L = fem::pack_constants(*L); auto coeffs_L = fem::allocate_coefficient_storage(*L); From 5e2c26ee1e4fd8c2ecf4f78a902ad7833a275340 Mon Sep 17 00:00:00 2001 From: Chris Richardson Date: Thu, 16 Oct 2025 10:19:32 +0100 Subject: [PATCH 3/3] Change auto to std::vector& --- src/elasticity_problem.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/elasticity_problem.cpp b/src/elasticity_problem.cpp index b9059e2..de32001 100644 --- a/src/elasticity_problem.cpp +++ b/src/elasticity_problem.cpp @@ -50,9 +50,9 @@ MatNullSpace build_near_nullspace(const fem::FunctionSpace& V) } // Rotations - auto x3 = basis[3].array(); - auto x4 = basis[4].array(); - auto x5 = basis[5].array(); + std::vector& x3 = basis[3].array(); + std::vector& x4 = basis[4].array(); + std::vector& x5 = basis[5].array(); const std::vector x = V.tabulate_dof_coordinates(false); const std::int32_t* dofs = V.dofmap()->map().data_handle();