Skip to content

Commit

Permalink
Unified Domain and Boundary kernel interface
Browse files Browse the repository at this point in the history
  • Loading branch information
@kubagalecki
kubagalecki committed Jan 3, 2024
1 parent 82b5803 commit 88477c7
Show file tree
Hide file tree
Showing 2 changed files with 34 additions and 149 deletions.
2 changes: 1 addition & 1 deletion include/l3ster/common/KernelInterface.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
181 changes: 33 additions & 148 deletions include/l3ster/glob_asm/AlgebraicSystem.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;

Expand Down Expand Up @@ -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");
Expand All @@ -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");
Expand All @@ -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 >
Expand Down Expand Up @@ -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()`");
Expand All @@ -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()
{
Expand Down

0 comments on commit 88477c7

Please sign in to comment.