From 88477c75a4b87053f08220388a23f71b05a6585f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jakub=20Ga=C5=82ecki?= Date: Wed, 3 Jan 2024 17:12:31 +0100 Subject: [PATCH] Unified Domain and Boundary kernel interface --- include/l3ster/common/KernelInterface.hpp | 2 +- include/l3ster/glob_asm/AlgebraicSystem.hpp | 181 ++++---------------- 2 files changed, 34 insertions(+), 149 deletions(-) diff --git a/include/l3ster/common/KernelInterface.hpp b/include/l3ster/common/KernelInterface.hpp index 31413e33..238f3fcb 100644 --- a/include/l3ster/common/KernelInterface.hpp +++ b/include/l3ster/common/KernelInterface.hpp @@ -222,6 +222,6 @@ template < typename T > concept EquationKernel_c = detail::is_domain_equation_kernel< T > or detail::is_boundary_equation_kernel< T >; template < typename T > -concept ResidualnKernel_c = detail::is_domain_residual_kernel< T > or detail::is_boundary_residual_kernel< T >; +concept ResidualKernel_c = detail::is_domain_residual_kernel< T > or detail::is_boundary_residual_kernel< T >; } // namespace lstr #endif // L3STER_COMMON_INTERFACE_HPP diff --git a/include/l3ster/glob_asm/AlgebraicSystem.hpp b/include/l3ster/glob_asm/AlgebraicSystem.hpp index 6554a4a5..4a7bbfbf 100644 --- a/include/l3ster/glob_asm/AlgebraicSystem.hpp +++ b/include/l3ster/glob_asm/AlgebraicSystem.hpp @@ -40,59 +40,33 @@ class AlgebraicSystem inline void beginAssembly(); inline void endAssembly(); - template < typename Kernel, - KernelParams params, + template < EquationKernel_c Kernel, ArrayOf_c< size_t > auto field_inds = util::makeIotaArray< size_t, max_dofs_per_node >(), size_t n_fields = 0, AssemblyOptions asm_opts = AssemblyOptions{} > - void assembleProblem(const DomainEquationKernel< Kernel, params >& kernel, + void assembleProblem(const Kernel& kernel, const util::ArrayOwner< d_id_t >& domain_ids, const SolutionManager::FieldValueGetter< n_fields >& fval_getter = {}, util::ConstexprValue< field_inds > field_inds_ctwrpr = {}, util::ConstexprValue< asm_opts > assembly_options = {}, val_t time = 0.) - requires(params.n_rhs == n_rhs); - template < typename Kernel, - KernelParams params, - ArrayOf_c< size_t > auto field_inds = util::makeIotaArray< size_t, max_dofs_per_node >(), - size_t n_fields = 0, - AssemblyOptions asm_opts = AssemblyOptions{} > - void assembleProblem(const BoundaryEquationKernel< Kernel, params >& kernel, - const util::ArrayOwner< d_id_t >& boundary_ids, - const SolutionManager::FieldValueGetter< n_fields >& fval_getter = {}, - util::ConstexprValue< field_inds > field_inds_ctwrpr = {}, - util::ConstexprValue< asm_opts > assembly_options = {}, - val_t time = 0.) - requires(params.n_rhs == n_rhs); - - template < typename Kernel, KernelParams params, std::integral dofind_t = size_t, size_t n_fields = 0 > - void setValues(const Teuchos::RCP< tpetra_femultivector_t >& solution, - const ResidualDomainKernel< Kernel, params >& kernel, - const util::ArrayOwner< d_id_t >& domain_ids, - const std::array< dofind_t, params.n_equations >& dof_inds, - const SolutionManager::FieldValueGetter< n_fields >& field_val_getter = {}, - val_t time = 0.) const; - template < typename Kernel, KernelParams params, std::integral dofind_t = size_t, size_t n_fields = 0 > - void setValues(const Teuchos::RCP< tpetra_femultivector_t >& solution, - const ResidualBoundaryKernel< Kernel, params >& kernel, - const util::ArrayOwner< d_id_t >& boundary_ids, - const std::array< dofind_t, params.n_equations >& dof_inds, - const SolutionManager::FieldValueGetter< n_fields >& field_val_getter = {}, - val_t time = 0.) const; + requires(Kernel::parameters.n_rhs == n_rhs); + + template < ResidualKernel_c Kernel, std::integral dofind_t = size_t, size_t n_fields = 0 > + void setValues(const Teuchos::RCP< tpetra_femultivector_t >& solution, + const Kernel& kernel, + const util::ArrayOwner< d_id_t >& domain_ids, + const std::array< dofind_t, Kernel::parameters.n_equations >& dof_inds, + const SolutionManager::FieldValueGetter< n_fields >& field_val_getter = {}, + val_t time = 0.) const; inline void applyDirichletBCs(); - template < typename Kernel, KernelParams params, std::integral dofind_t = size_t, size_t n_fields = 0 > - void setDirichletBCValues(const ResidualDomainKernel< Kernel, params >& kernel, - const util::ArrayOwner< d_id_t >& domain_ids, - const std::array< dofind_t, params.n_equations >& dof_inds, - const SolutionManager::FieldValueGetter< n_fields >& field_val_getter = {}, - val_t time = 0.); - template < typename Kernel, KernelParams params, std::integral dofind_t = size_t, size_t n_fields = 0 > - void setDirichletBCValues(const ResidualBoundaryKernel< Kernel, params >& kernel, - const util::ArrayOwner< d_id_t >& boundary_ids, - const std::array< dofind_t, params.n_equations >& dof_inds, - const SolutionManager::FieldValueGetter< n_fields >& field_val_getter = {}, - val_t time = 0.); + 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, + const std::array< dofind_t, Kernel::parameters.n_equations >& dof_inds, + const SolutionManager::FieldValueGetter< n_fields >& field_val_getter = {}, + val_t time = 0.); void solve(solvers::Solver_c auto& solver, const Teuchos::RCP< tpetra_femultivector_t >& solution) const; @@ -162,14 +136,14 @@ void AlgebraicSystem< max_dofs_per_node, CP, n_rhs, orders... >::updateSolution( } template < size_t max_dofs_per_node, CondensationPolicy CP, size_t n_rhs, el_o_t... orders > -template < typename Kernel, KernelParams params, std::integral dofind_t, size_t n_fields > +template < ResidualKernel_c Kernel, std::integral dofind_t, size_t n_fields > void AlgebraicSystem< max_dofs_per_node, CP, n_rhs, orders... >::setValues( - const Teuchos::RCP< tpetra_femultivector_t >& solution, - const ResidualDomainKernel< Kernel, params >& kernel, - const util::ArrayOwner< d_id_t >& domain_ids, - const std::array< dofind_t, params.n_equations >& dof_inds, - const SolutionManager::FieldValueGetter< n_fields >& field_val_getter, - val_t time) const + const Teuchos::RCP< tpetra_femultivector_t >& solution, + const Kernel& kernel, + const util::ArrayOwner< d_id_t >& domain_ids, + const std::array< dofind_t, Kernel::parameters.n_equations >& dof_inds, + const SolutionManager::FieldValueGetter< n_fields >& field_val_getter, + val_t time) const { util::throwingAssert(util::isValidIndexRange(dof_inds, max_dofs_per_node), "The DOF indices are out of bounds for the problem"); @@ -190,41 +164,13 @@ void AlgebraicSystem< max_dofs_per_node, CP, n_rhs, orders... >::setValues( } template < size_t max_dofs_per_node, CondensationPolicy CP, size_t n_rhs, el_o_t... orders > -template < typename Kernel, KernelParams params, std::integral dofind_t, size_t n_fields > -void AlgebraicSystem< max_dofs_per_node, CP, n_rhs, orders... >::setValues( - const Teuchos::RCP< tpetra_femultivector_t >& solution, - const ResidualBoundaryKernel< Kernel, params >& kernel, - const util::ArrayOwner< d_id_t >& boundary_ids, - const std::array< dofind_t, params.n_equations >& dof_inds, - const SolutionManager::FieldValueGetter< n_fields >& field_val_getter, - val_t time) const -{ - util::throwingAssert(util::isValidIndexRange(dof_inds, max_dofs_per_node), - "The DOF indices are out of bounds for the problem"); - util::throwingAssert(n_rhs == solution->getNumVectors(), "The number of columns and number of RHS must match"); - - const auto vals_view = solution->getLocalViewHost(Tpetra::Access::ReadWrite); - computeValuesAtNodes(kernel, - *m_mesh, - boundary_ids, - m_node_dof_map, - dof_inds, - field_val_getter, - util::asSpans< n_rhs >(vals_view), - time); - solution->switchActiveMultiVector(); - solution->doOwnedToOwnedPlusShared(Tpetra::CombineMode::REPLACE); - solution->switchActiveMultiVector(); -} - -template < size_t max_dofs_per_node, CondensationPolicy CP, size_t n_rhs, el_o_t... orders > -template < typename Kernel, KernelParams params, std::integral dofind_t, size_t n_fields > +template < ResidualKernel_c Kernel, std::integral dofind_t, size_t n_fields > void AlgebraicSystem< max_dofs_per_node, CP, n_rhs, orders... >::setDirichletBCValues( - const ResidualDomainKernel< Kernel, params >& kernel, - const util::ArrayOwner< d_id_t >& domain_ids, - const std::array< dofind_t, params.n_equations >& dof_inds, - const SolutionManager::FieldValueGetter< n_fields >& field_val_getter, - val_t time) + const Kernel& kernel, + const util::ArrayOwner< d_id_t >& domain_ids, + const std::array< dofind_t, Kernel::parameters.n_equations >& dof_inds, + const SolutionManager::FieldValueGetter< n_fields >& field_val_getter, + val_t time) { util::throwingAssert(util::isValidIndexRange(dof_inds, max_dofs_per_node), "The DOF indices are out of bounds for the problem"); @@ -240,29 +186,6 @@ void AlgebraicSystem< max_dofs_per_node, CP, n_rhs, orders... >::setDirichletBCV time); } -template < size_t max_dofs_per_node, CondensationPolicy CP, size_t n_rhs, el_o_t... orders > -template < typename Kernel, KernelParams params, std::integral dofind_t, size_t n_fields > -void AlgebraicSystem< max_dofs_per_node, CP, n_rhs, orders... >::setDirichletBCValues( - const ResidualBoundaryKernel< Kernel, params >& kernel, - const util::ArrayOwner< d_id_t >& boundary_ids, - const std::array< dofind_t, params.n_equations >& dof_inds, - const SolutionManager::FieldValueGetter< n_fields >& field_val_getter, - val_t time) -{ - util::throwingAssert(util::isValidIndexRange(dof_inds, max_dofs_per_node), - "The DOF indices are out of bounds for the problem"); - - const auto vals_view = m_dirichlet_values->getLocalViewHost(Tpetra::Access::OverwriteAll); - computeValuesAtNodes(kernel, - *m_mesh, - boundary_ids, - m_node_dof_map, - dof_inds, - field_val_getter, - util::asSpans< n_rhs >(vals_view), - time); -} - template < size_t max_dofs_per_node, CondensationPolicy CP, size_t n_rhs, el_o_t... orders > auto AlgebraicSystem< max_dofs_per_node, CP, n_rhs, orders... >::getMatrix() const -> Teuchos::RCP< const tpetra_crsmatrix_t > @@ -375,19 +298,15 @@ void AlgebraicSystem< max_dofs_per_node, CP, n_rhs, orders... >::setToZero() } template < size_t max_dofs_per_node, CondensationPolicy CP, size_t n_rhs, el_o_t... orders > -template < typename Kernel, - KernelParams params, - ArrayOf_c< size_t > auto field_inds, - size_t n_fields, - AssemblyOptions asm_opts > +template < EquationKernel_c Kernel, ArrayOf_c< size_t > auto field_inds, size_t n_fields, AssemblyOptions asm_opts > void AlgebraicSystem< max_dofs_per_node, CP, n_rhs, orders... >::assembleProblem( - const DomainEquationKernel< Kernel, params >& kernel, + const Kernel& kernel, const util::ArrayOwner< d_id_t >& domain_ids, const SolutionManager::FieldValueGetter< n_fields >& fval_getter, util::ConstexprValue< field_inds > field_inds_ctwrpr, util::ConstexprValue< asm_opts > assembly_options, val_t time) - requires(params.n_rhs == n_rhs) + requires(Kernel::parameters.n_rhs == n_rhs) { L3STER_PROFILE_FUNCTION; assertState(State::OpenForAssembly, "`assembleProblem()` was called before `beginAssembly()`"); @@ -408,40 +327,6 @@ void AlgebraicSystem< max_dofs_per_node, CP, n_rhs, orders... >::assembleProblem #endif } -template < size_t max_dofs_per_node, CondensationPolicy CP, size_t n_rhs, el_o_t... orders > -template < typename Kernel, - KernelParams params, - ArrayOf_c< size_t > auto field_inds, - size_t n_fields, - AssemblyOptions asm_opts > -void AlgebraicSystem< max_dofs_per_node, CP, n_rhs, orders... >::assembleProblem( - const BoundaryEquationKernel< Kernel, params >& kernel, - const util::ArrayOwner< d_id_t >& boundary_ids, - const SolutionManager::FieldValueGetter< n_fields >& fval_getter, - util::ConstexprValue< field_inds > field_inds_ctwrpr, - util::ConstexprValue< asm_opts > assembly_options, - val_t time) - requires(params.n_rhs == n_rhs) -{ - L3STER_PROFILE_FUNCTION; - assertState(State::OpenForAssembly, "`assembleProblem()` was called before `beginAssembly()`"); - assembleGlobalSystem(kernel, - *m_mesh, - boundary_ids, - fval_getter, - *m_matrix, - m_rhs_views, - m_node_dof_map, - m_condensation_manager, - field_inds_ctwrpr, - assembly_options, - time); - -#ifdef L3STER_PROFILE_EXECUTION - m_rhs->getMap()->getComm()->barrier(); -#endif -} - template < size_t max_dofs_per_node, CondensationPolicy CP, size_t n_rhs, el_o_t... orders > void AlgebraicSystem< max_dofs_per_node, CP, n_rhs, orders... >::applyDirichletBCs() {