applyDirichletBCs now automatically called in endAssembly

benchmarks/Diffusion3DBenchmark.cpp

@@ -113,10 +113,9 @@ int main(int argc, char* argv[])
     auto alg_system = makeAlgebraicSystem(comm, my_partition, probdef_ctwrpr, dirichletdef_ctwrpr, algpar_ctwrpr);
     alg_system.beginAssembly();
     alg_system.assembleProblem(diffusion_kernel3d, {domain_id});
+    alg_system.setDirichletBCValues(dirichlet_bc_kernel, boundary_ids, T_inds);
     alg_system.endAssembly();
     alg_system.describe(comm);
-    alg_system.setDirichletBCValues(dirichlet_bc_kernel, boundary_ids, T_inds);
-    alg_system.applyDirichletBCs();
 
     constexpr auto solver_opts  = IterSolverOpts{.verbosity = {.summary = true, .timing = true}};
     constexpr auto precond_opts = ChebyshevOpts{.degree = 3};
@@ -131,7 +130,7 @@ int main(int argc, char* argv[])
 
     L3STER_PROFILE_REGION_BEGIN("Compute solution error");
     const auto fval_getter = solution_manager.makeFieldValueGetter(field_inds);
-    const auto error = computeNormL2(comm, error_kernel, *my_partition, {domain_id}, fval_getter);
+    const auto error       = computeNormL2(comm, error_kernel, *my_partition, {domain_id}, fval_getter);
     L3STER_PROFILE_REGION_END("Compute solution error");
 
     if (comm.getRank() == 0)

benchmarks/UtilBenchmarks.cpp

@@ -1,6 +1,6 @@
 #include "Common.hpp"
 
-static auto kernelImpl(size_t x, size_t rand) -> size_t
+static auto kernelImpl(size_t in, size_t rand) -> size_t
 {
     thread_local const auto hash = std::hash< size_t >{};
     return hash(hash(hash(in) % rand)) * 953u + rand + 1u;

include/l3ster/glob_asm/AlgebraicSystem.hpp

@@ -60,7 +60,6 @@ class AlgebraicSystem
                    const SolutionManager::FieldValueGetter< n_fields >&          field_val_getter = {},
                    val_t                                                         time             = 0.) const;
 
-    inline void applyDirichletBCs();
     template < ResidualKernel_c Kernel, std::integral dofind_t = size_t, size_t n_fields = 0 >
     void setDirichletBCValues(const Kernel&                                                 kernel,
                               const util::ArrayOwner< d_id_t >&                             domain_ids,
@@ -74,6 +73,7 @@ class AlgebraicSystem
 
 private:
     inline void setToZero();
+    inline void applyDirichletBCs();
 
     enum struct State
     {
@@ -247,7 +247,7 @@ AlgebraicSystem< max_dofs_per_node, CP, n_rhs, orders... >::AlgebraicSystem(
         L3STER_PROFILE_REGION_END("Dirichlet BCs");
     }
     if (m_dirichlet_bcs.has_value())
-        m_dirichlet_values = util::makeTeuchosRCP< tpetra_multivector_t >(m_sparsity_graph->getRowMap(), 1u);
+        m_dirichlet_values = util::makeTeuchosRCP< tpetra_multivector_t >(m_sparsity_graph->getRowMap(), n_rhs);
 }
 
 template < size_t max_dofs_per_node, CondensationPolicy CP, size_t n_rhs, el_o_t... orders >
@@ -288,6 +288,8 @@ void AlgebraicSystem< max_dofs_per_node, CP, n_rhs, orders... >::endAssembly()
 #ifdef L3STER_PROFILE_EXECUTION
     m_rhs->getMap()->getComm()->barrier();
 #endif
+    if (m_dirichlet_bcs.has_value())
+        applyDirichletBCs();
 }
 
 template < size_t max_dofs_per_node, CondensationPolicy CP, size_t n_rhs, el_o_t... orders >
@@ -331,10 +333,6 @@ template < size_t max_dofs_per_node, CondensationPolicy CP, size_t n_rhs, el_o_t
 void AlgebraicSystem< max_dofs_per_node, CP, n_rhs, orders... >::applyDirichletBCs()
 {
     L3STER_PROFILE_FUNCTION;
-    util::throwingAssert(m_dirichlet_bcs.has_value(),
-                         "`applyDirichletBCs()` was called, but no Dirichlet BCs were defined");
-    assertState(State::Closed, "`applyDirichletBCs()` was called before `endAssembly()`");
-
     m_matrix->beginModify();
     m_rhs->beginModify();
     m_dirichlet_bcs->apply(*m_dirichlet_values, *m_matrix, *m_rhs);

tests/Diffusion2DTest.cpp

@@ -16,8 +16,7 @@ constexpr auto node_dist = std::array{0., 1., 2., 3., 4., 5., 6.};
 auto makeMesh(const MpiComm& comm, auto probdef_ctwrpr)
 {
     constexpr auto mesh_order = 2;
-    return generateAndDistributeMesh< mesh_order >(
-        comm, [&] { return makeSquareMesh(node_dist); }, {}, probdef_ctwrpr);
+    return generateAndDistributeMesh< mesh_order >(comm, [&] { return makeSquareMesh(node_dist); }, {}, probdef_ctwrpr);
 }
 
 template < CondensationPolicy CP >
@@ -76,25 +75,22 @@ void test()
         alg_sys.assembleProblem(neumann_bc_kernel, adiabatic_bound_ids);
     };
 
+    constexpr auto dirichlet_bound_ids = std::array{left_boundary, right_boundary};
+    alg_sys.setDirichletBCValues(dirichlet_bc_kernel, dirichlet_bound_ids, std::array{0});
+
     // Check constraints on assembly state
     alg_sys.beginAssembly();
     assembleDomainProblem();
     assembleBoundaryProblem();
-    CHECK_THROWS(alg_sys.applyDirichletBCs());
     alg_sys.endAssembly();
     alg_sys.describe(comm);
     CHECK_THROWS(assembleDomainProblem());
     CHECK_THROWS(assembleBoundaryProblem());
     CHECK_THROWS(alg_sys.endAssembly());
 
-    constexpr auto dirichlet_bound_ids = std::array{left_boundary, right_boundary};
-    alg_sys.setDirichletBCValues(dirichlet_bc_kernel, dirichlet_bound_ids, std::array{0});
-    alg_sys.applyDirichletBCs();
-
     {
         auto dummy_problem = makeAlgebraicSystem(comm, mesh, probdef_ctwrpr, {}, algparams_ctwrpr);
         dummy_problem.endAssembly();
-        CHECK_THROWS(dummy_problem.applyDirichletBCs());
     }
 
     constexpr auto dof_inds = util::makeIotaArray< size_t, 3 >();