Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

feature: Prepare to introduce RAJA for all loops #1018

Merged
merged 14 commits into from
Nov 3, 2023
Merged

feature: Prepare to introduce RAJA for all loops #1018

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
14 commits
Select commit Hold shift + click to select a range
65f8fbe
Refactor extended lambda and enable RAJA usage inside domain integral…
johnbowen42 Sep 5, 2023
3aef667
Add static assert for generic lambdas
johnbowen42 Sep 7, 2023
6fa35e5
Add docs string and limit visiblity of class used by static assert.
johnbowen42 Sep 7, 2023
b30e904
Exercise basic CUDA enabled functionality with RAJA
johnbowen42 Oct 18, 2023
3a204ea
Clang format
johnbowen42 Oct 18, 2023
cf6125e
Cleanup print statements and commented code
johnbowen42 Oct 18, 2023
830a7df
Separate out host device declarations and template parameter deductio…
johnbowen42 Oct 25, 2023
55ad49f
Fix build issues
johnbowen42 Oct 25, 2023
d1fcd28
Add maybe unused flag to empty trial element
johnbowen42 Oct 26, 2023
37eb4e1
Add comments explaining trial_elements_tuple and refactor for clarity…
johnbowen42 Nov 1, 2023
5931dc5
Rename some variables and types
johnbowen42 Nov 1, 2023
cf090f5
Move variables back to constexpr
johnbowen42 Nov 2, 2023
a1d55a8
style
johnbowen42 Nov 2, 2023
98a5462
Merge branch 'develop' into feature/bowen/refactor-for-extended-host-…
jamiebramwell Nov 3, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
20 changes: 9 additions & 11 deletions src/serac/numerics/functional/boundary_integral_kernels.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -153,18 +153,14 @@ auto batch_apply_qf(lambda qf, const tensor<double, 3, n>& positions, const tens
return outputs;
}

template <uint32_t differentiation_index, int Q, mfem::Geometry::Type geom, typename test, typename... trials,
typename lambda_type, typename derivative_type, int... indices>
void evaluation_kernel_impl(FunctionSignature<test(trials...)>, const std::vector<const double*>& inputs,
/// @trial_elements the element type for each trial space
template <uint32_t differentiation_index, int Q, mfem::Geometry::Type geom, typename test_element,
typename trial_element_type, typename lambda_type, typename derivative_type, int... indices>
void evaluation_kernel_impl(trial_element_type trial_elements, test_element, const std::vector<const double*>& inputs,
double* outputs, const double* positions, const double* jacobians, lambda_type qf,
[[maybe_unused]] derivative_type* qf_derivatives, uint32_t num_elements,
std::integer_sequence<int, indices...>)
camp::int_seq<int, indices...>)
johnbowen42 marked this conversation as resolved.
Show resolved Hide resolved
{
using test_element = finite_element<geom, test>;

/// @brief the element type for each trial space
static constexpr tuple<finite_element<geom, trials>...> trial_elements{};

// mfem provides this information as opaque arrays of doubles,
// so we reinterpret the pointer with
constexpr int dim = dimension_of(geom) + 1;
Expand Down Expand Up @@ -331,9 +327,11 @@ std::function<void(const std::vector<const double*>&, double*, bool)> evaluation
signature s, lambda_type qf, const double* positions, const double* jacobians,
std::shared_ptr<derivative_type> qf_derivatives, uint32_t num_elements)
{
auto trial_elements = trial_elements_tuple<geom>(s);
auto test = get_test_element<geom>(s);
return [=](const std::vector<const double*>& inputs, double* outputs, bool /* update state */) {
evaluation_kernel_impl<wrt, Q, geom>(s, inputs, outputs, positions, jacobians, qf, qf_derivatives.get(),
num_elements, s.index_seq);
evaluation_kernel_impl<wrt, Q, geom>(trial_elements, test, inputs, outputs, positions, jacobians, qf,
qf_derivatives.get(), num_elements, s.index_seq);
};
}

Expand Down
7 changes: 4 additions & 3 deletions src/serac/numerics/functional/detail/hexahedron_H1.inl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -167,7 +167,7 @@ struct finite_element<mfem::Geometry::CUBE, H1<p, c> > {
}

template <int q>
static auto interpolate(const dof_type& X, const TensorProductQuadratureRule<q>&)
SERAC_HOST_DEVICE static auto interpolate(const dof_type& X, const TensorProductQuadratureRule<q>&)
{
// we want to compute the following:
//
Expand Down Expand Up @@ -228,8 +228,9 @@ struct finite_element<mfem::Geometry::CUBE, H1<p, c> > {
}

template <typename source_type, typename flux_type, int q>
static void integrate(const tensor<tuple<source_type, flux_type>, q * q * q>& qf_output,
const TensorProductQuadratureRule<q>&, dof_type* element_residual, int step = 1)
SERAC_HOST_DEVICE static void integrate(const tensor<tuple<source_type, flux_type>, q * q * q>& qf_output,
const TensorProductQuadratureRule<q>&, dof_type* element_residual,
int step = 1)
{
if constexpr (is_zero<source_type>{} && is_zero<flux_type>{}) {
return;
Expand Down
8 changes: 4 additions & 4 deletions src/serac/numerics/functional/detail/hexahedron_Hcurl.inl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -334,7 +334,7 @@ struct finite_element<mfem::Geometry::CUBE, Hcurl<p>> {
}

template <int q>
static auto interpolate(const dof_type& element_values, const TensorProductQuadratureRule<q>&)
SERAC_HOST_DEVICE static auto interpolate(const dof_type& element_values, const TensorProductQuadratureRule<q>&)
{
constexpr bool apply_weights = false;
constexpr tensor<double, q, p> B1 = calculate_B1<apply_weights, q>();
Expand Down Expand Up @@ -395,9 +395,9 @@ struct finite_element<mfem::Geometry::CUBE, Hcurl<p>> {
}

template <typename source_type, typename flux_type, int q>
static void integrate(const tensor<tuple<source_type, flux_type>, q * q * q>& qf_output,
const TensorProductQuadratureRule<q>&, dof_type* element_residual,
[[maybe_unused]] int step = 1)
SERAC_HOST_DEVICE static void integrate(const tensor<tuple<source_type, flux_type>, q * q * q>& qf_output,
const TensorProductQuadratureRule<q>&, dof_type* element_residual,
[[maybe_unused]] int step = 1)
{
constexpr bool apply_weights = true;
constexpr tensor<double, q, p> B1 = calculate_B1<apply_weights, q>();
Expand Down
7 changes: 4 additions & 3 deletions src/serac/numerics/functional/detail/hexahedron_L2.inl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -171,7 +171,7 @@ struct finite_element<mfem::Geometry::CUBE, L2<p, c> > {
}

template <int q>
static auto interpolate(const dof_type& X, const TensorProductQuadratureRule<q>&)
SERAC_HOST_DEVICE static auto interpolate(const dof_type& X, const TensorProductQuadratureRule<q>&)
{
// we want to compute the following:
//
Expand Down Expand Up @@ -232,8 +232,9 @@ struct finite_element<mfem::Geometry::CUBE, L2<p, c> > {
}

template <typename source_type, typename flux_type, int q>
static void integrate(const tensor<tuple<source_type, flux_type>, q * q * q>& qf_output,
const TensorProductQuadratureRule<q>&, dof_type* element_residual, int step = 1)
SERAC_HOST_DEVICE static void integrate(const tensor<tuple<source_type, flux_type>, q * q * q>& qf_output,
const TensorProductQuadratureRule<q>&, dof_type* element_residual,
int step = 1)
{
if constexpr (is_zero<source_type>{} && is_zero<flux_type>{}) {
return;
Expand Down
13 changes: 7 additions & 6 deletions src/serac/numerics/functional/detail/qoi.inl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -24,15 +24,15 @@ struct finite_element<g, QOI> {
SERAC_HOST_DEVICE static constexpr double shape_functions(double /* xi */) { return 1.0; }

template <int Q, int q>
static void integrate(const tensor<zero, Q>&, const TensorProductQuadratureRule<q>&, dof_type*,
[[maybe_unused]] int step = 1)
SERAC_HOST_DEVICE static void integrate(const tensor<zero, Q>&, const TensorProductQuadratureRule<q>&, dof_type*,
[[maybe_unused]] int step = 1)
{
return; // integrating zeros is a no-op
}

template <int Q, int q>
static void integrate(const tensor<double, Q>& qf_output, const TensorProductQuadratureRule<q>&,
dof_type* element_total, [[maybe_unused]] int step = 1)
SERAC_HOST_DEVICE static void integrate(const tensor<double, Q>& qf_output, const TensorProductQuadratureRule<q>&,
dof_type* element_total, [[maybe_unused]] int step = 1)
{
if constexpr (geometry == mfem::Geometry::SEGMENT) {
static_assert(Q == q);
Expand Down Expand Up @@ -77,8 +77,9 @@ struct finite_element<g, QOI> {
// this overload is used for boundary integrals, since they pad the
// output to be a tuple with a hardcoded `zero` flux term
template <typename source_type, int Q, int q>
static void integrate(const tensor<serac::tuple<source_type, zero>, Q>& qf_output,
const TensorProductQuadratureRule<q>&, dof_type* element_total, [[maybe_unused]] int step = 1)
SERAC_HOST_DEVICE static void integrate(const tensor<serac::tuple<source_type, zero>, Q>& qf_output,
const TensorProductQuadratureRule<q>&, dof_type* element_total,
[[maybe_unused]] int step = 1)
{
if constexpr (is_zero<source_type>{}) {
return;
Expand Down
7 changes: 4 additions & 3 deletions src/serac/numerics/functional/detail/quadrilateral_H1.inl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -202,7 +202,7 @@ struct finite_element<mfem::Geometry::SQUARE, H1<p, c> > {
// A(dy, qx) := B(qx, dx) * X_e(dy, dx)
// X_q(qy, qx) := B(qy, dy) * A(dy, qx)
template <int q>
static auto interpolate(const dof_type& X, const TensorProductQuadratureRule<q>&)
SERAC_HOST_DEVICE static auto interpolate(const dof_type& X, const TensorProductQuadratureRule<q>&)
{
static constexpr bool apply_weights = false;
static constexpr auto B = calculate_B<apply_weights, q>();
Expand Down Expand Up @@ -245,8 +245,9 @@ struct finite_element<mfem::Geometry::SQUARE, H1<p, c> > {
// flux can be one of: {zero, tensor<double,dim>, tensor<double,dim,dim>, tensor<double,dim,dim,dim>,
// tensor<double,dim,dim,dim>}
template <typename source_type, typename flux_type, int q>
static void integrate(const tensor<tuple<source_type, flux_type>, q * q>& qf_output,
const TensorProductQuadratureRule<q>&, dof_type* element_residual, int step = 1)
SERAC_HOST_DEVICE static void integrate(const tensor<tuple<source_type, flux_type>, q * q>& qf_output,
const TensorProductQuadratureRule<q>&, dof_type* element_residual,
int step = 1)
{
if constexpr (is_zero<source_type>{} && is_zero<flux_type>{}) {
return;
Expand Down
8 changes: 4 additions & 4 deletions src/serac/numerics/functional/detail/quadrilateral_Hcurl.inl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -285,7 +285,7 @@ struct finite_element<mfem::Geometry::SQUARE, Hcurl<p> > {
}

template <int q>
static auto interpolate(const dof_type& element_values, const TensorProductQuadratureRule<q>&)
SERAC_HOST_DEVICE static auto interpolate(const dof_type& element_values, const TensorProductQuadratureRule<q>&)
{
constexpr bool apply_weights = false;
constexpr tensor<double, q, p> B1 = calculate_B1<apply_weights, q>();
Expand Down Expand Up @@ -323,9 +323,9 @@ struct finite_element<mfem::Geometry::SQUARE, Hcurl<p> > {
}

template <typename source_type, typename flux_type, int q>
static void integrate(const tensor<tuple<source_type, flux_type>, q * q>& qf_output,
const TensorProductQuadratureRule<q>&, dof_type* element_residual,
[[maybe_unused]] int step = 1)
SERAC_HOST_DEVICE static void integrate(const tensor<tuple<source_type, flux_type>, q * q>& qf_output,
const TensorProductQuadratureRule<q>&, dof_type* element_residual,
[[maybe_unused]] int step = 1)
{
constexpr bool apply_weights = true;
constexpr tensor<double, q, p> B1 = calculate_B1<apply_weights, q>();
Expand Down
7 changes: 4 additions & 3 deletions src/serac/numerics/functional/detail/quadrilateral_L2.inl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -201,7 +201,7 @@ struct finite_element<mfem::Geometry::SQUARE, L2<p, c> > {
// A(dy, qx) := B(qx, dx) * X_e(dy, dx)
// X_q(qy, qx) := B(qy, dy) * A(dy, qx)
template <int q>
static auto interpolate(const dof_type& X, const TensorProductQuadratureRule<q>&)
SERAC_HOST_DEVICE static auto interpolate(const dof_type& X, const TensorProductQuadratureRule<q>&)
{
static constexpr bool apply_weights = false;
static constexpr auto B = calculate_B<apply_weights, q>();
Expand Down Expand Up @@ -244,8 +244,9 @@ struct finite_element<mfem::Geometry::SQUARE, L2<p, c> > {
// flux can be one of: {zero, tensor<double,dim>, tensor<double,dim,dim>, tensor<double,dim,dim,dim>,
// tensor<double,dim,dim,dim>}
template <typename source_type, typename flux_type, int q>
static void integrate(const tensor<tuple<source_type, flux_type>, q * q>& qf_output,
const TensorProductQuadratureRule<q>&, dof_type* element_residual, int step = 1)
SERAC_HOST_DEVICE static void integrate(const tensor<tuple<source_type, flux_type>, q * q>& qf_output,
const TensorProductQuadratureRule<q>&, dof_type* element_residual,
int step = 1)
{
if constexpr (is_zero<source_type>{} && is_zero<flux_type>{}) {
return;
Expand Down
10 changes: 5 additions & 5 deletions src/serac/numerics/functional/detail/segment_H1.inl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -94,7 +94,7 @@ struct finite_element<mfem::Geometry::SEGMENT, H1<p, c> > {
}

template <typename T, int q>
static auto batch_apply_shape_fn(int jx, tensor<T, q> input, const TensorProductQuadratureRule<q>&)
SERAC_HOST_DEVICE static auto batch_apply_shape_fn(int jx, tensor<T, q> input, const TensorProductQuadratureRule<q>&)
{
static constexpr bool apply_weights = false;
static constexpr auto B = calculate_B<apply_weights, q>();
Expand All @@ -121,7 +121,7 @@ struct finite_element<mfem::Geometry::SEGMENT, H1<p, c> > {
}

template <int q>
static auto interpolate(const dof_type& X, const TensorProductQuadratureRule<q>&)
SERAC_HOST_DEVICE static auto interpolate(const dof_type& X, const TensorProductQuadratureRule<q>&)
{
static constexpr bool apply_weights = false;
static constexpr auto B = calculate_B<apply_weights, q>();
Expand Down Expand Up @@ -155,9 +155,9 @@ struct finite_element<mfem::Geometry::SEGMENT, H1<p, c> > {
}

template <typename source_type, typename flux_type, int q>
static void integrate(const tensor<tuple<source_type, flux_type>, q>& qf_output,
const TensorProductQuadratureRule<q>&, dof_type* element_residual,
[[maybe_unused]] int step = 1)
SERAC_HOST_DEVICE static void integrate(const tensor<tuple<source_type, flux_type>, q>& qf_output,
const TensorProductQuadratureRule<q>&, dof_type* element_residual,
[[maybe_unused]] int step = 1)
{
if constexpr (is_zero<source_type>{} && is_zero<flux_type>{}) {
return;
Expand Down
8 changes: 4 additions & 4 deletions src/serac/numerics/functional/detail/segment_L2.inl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -121,7 +121,7 @@ struct finite_element<mfem::Geometry::SEGMENT, L2<p, c> > {
}

template <int q>
static auto interpolate(const dof_type& X, const TensorProductQuadratureRule<q>&)
SERAC_HOST_DEVICE static auto interpolate(const dof_type& X, const TensorProductQuadratureRule<q>&)
{
static constexpr bool apply_weights = false;
static constexpr auto B = calculate_B<apply_weights, q>();
Expand Down Expand Up @@ -155,9 +155,9 @@ struct finite_element<mfem::Geometry::SEGMENT, L2<p, c> > {
}

template <typename source_type, typename flux_type, int q>
static void integrate(const tensor<tuple<source_type, flux_type>, q>& qf_output,
const TensorProductQuadratureRule<q>&, dof_type* element_residual,
[[maybe_unused]] int step = 1)
SERAC_HOST_DEVICE static void integrate(const tensor<tuple<source_type, flux_type>, q>& qf_output,
const TensorProductQuadratureRule<q>&, dof_type* element_residual,
[[maybe_unused]] int step = 1)
{
if constexpr (is_zero<source_type>{} && is_zero<flux_type>{}) {
return;
Expand Down
7 changes: 4 additions & 3 deletions src/serac/numerics/functional/detail/tetrahedron_H1.inl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -361,7 +361,7 @@ struct finite_element<mfem::Geometry::TETRAHEDRON, H1<p, c> > {
}

template <int q>
static auto interpolate(const tensor<double, c, ndof>& X, const TensorProductQuadratureRule<q>&)
SERAC_HOST_DEVICE static auto interpolate(const tensor<double, c, ndof>& X, const TensorProductQuadratureRule<q>&)
{
constexpr auto xi = GaussLegendreNodes<q, mfem::Geometry::TETRAHEDRON>();

Expand All @@ -384,8 +384,9 @@ struct finite_element<mfem::Geometry::TETRAHEDRON, H1<p, c> > {
}

template <typename source_type, typename flux_type, int q>
static void integrate(const tensor<tuple<source_type, flux_type>, nqpts(q)>& qf_output,
const TensorProductQuadratureRule<q>&, tensor<double, c, ndof>* element_residual, int step = 1)
SERAC_HOST_DEVICE static void integrate(const tensor<tuple<source_type, flux_type>, nqpts(q)>& qf_output,
const TensorProductQuadratureRule<q>&,
tensor<double, c, ndof>* element_residual, int step = 1)
{
if constexpr (is_zero<source_type>{} && is_zero<flux_type>{}) {
return;
Expand Down
7 changes: 4 additions & 3 deletions src/serac/numerics/functional/detail/tetrahedron_L2.inl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -366,7 +366,7 @@ struct finite_element<mfem::Geometry::TETRAHEDRON, L2<p, c> > {
}

template <int q>
static auto interpolate(const tensor<double, c, ndof>& X, const TensorProductQuadratureRule<q>&)
SERAC_HOST_DEVICE static auto interpolate(const tensor<double, c, ndof>& X, const TensorProductQuadratureRule<q>&)
{
constexpr auto xi = GaussLegendreNodes<q, mfem::Geometry::TETRAHEDRON>();

Expand All @@ -389,8 +389,9 @@ struct finite_element<mfem::Geometry::TETRAHEDRON, L2<p, c> > {
}

template <typename source_type, typename flux_type, int q>
static void integrate(const tensor<tuple<source_type, flux_type>, nqpts(q)>& qf_output,
const TensorProductQuadratureRule<q>&, tensor<double, c, ndof>* element_residual, int step = 1)
SERAC_HOST_DEVICE static void integrate(const tensor<tuple<source_type, flux_type>, nqpts(q)>& qf_output,
const TensorProductQuadratureRule<q>&,
tensor<double, c, ndof>* element_residual, int step = 1)
{
if constexpr (is_zero<source_type>{} && is_zero<flux_type>{}) {
return;
Expand Down
7 changes: 4 additions & 3 deletions src/serac/numerics/functional/detail/triangle_H1.inl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -266,7 +266,7 @@ struct finite_element<mfem::Geometry::TRIANGLE, H1<p, c> > {
}

template <int q>
static auto interpolate(const tensor<double, c, ndof>& X, const TensorProductQuadratureRule<q>&)
SERAC_HOST_DEVICE static auto interpolate(const tensor<double, c, ndof>& X, const TensorProductQuadratureRule<q>&)
{
constexpr auto xi = GaussLegendreNodes<q, mfem::Geometry::TRIANGLE>();
static constexpr int num_quadrature_points = q * (q + 1) / 2;
Expand All @@ -290,8 +290,9 @@ struct finite_element<mfem::Geometry::TRIANGLE, H1<p, c> > {
}

template <typename source_type, typename flux_type, int q>
static void integrate(const tensor<tuple<source_type, flux_type>, q*(q + 1) / 2>& qf_output,
const TensorProductQuadratureRule<q>&, tensor<double, c, ndof>* element_residual, int step = 1)
SERAC_HOST_DEVICE static void integrate(const tensor<tuple<source_type, flux_type>, q*(q + 1) / 2>& qf_output,
const TensorProductQuadratureRule<q>&,
tensor<double, c, ndof>* element_residual, int step = 1)
{
if constexpr (is_zero<source_type>{} && is_zero<flux_type>{}) {
return;
Expand Down
7 changes: 4 additions & 3 deletions src/serac/numerics/functional/detail/triangle_L2.inl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -275,7 +275,7 @@ struct finite_element<mfem::Geometry::TRIANGLE, L2<p, c> > {
}

template <int q>
static auto interpolate(const tensor<double, c, ndof>& X, const TensorProductQuadratureRule<q>&)
SERAC_HOST_DEVICE static auto interpolate(const tensor<double, c, ndof>& X, const TensorProductQuadratureRule<q>&)
{
constexpr auto xi = GaussLegendreNodes<q, mfem::Geometry::TRIANGLE>();
static constexpr int num_quadrature_points = q * (q + 1) / 2;
Expand All @@ -299,8 +299,9 @@ struct finite_element<mfem::Geometry::TRIANGLE, L2<p, c> > {
}

template <typename source_type, typename flux_type, int q>
static void integrate(const tensor<tuple<source_type, flux_type>, q*(q + 1) / 2>& qf_output,
const TensorProductQuadratureRule<q>&, tensor<double, c, ndof>* element_residual, int step = 1)
SERAC_HOST_DEVICE static void integrate(const tensor<tuple<source_type, flux_type>, q*(q + 1) / 2>& qf_output,
const TensorProductQuadratureRule<q>&,
tensor<double, c, ndof>* element_residual, int step = 1)
{
if constexpr (is_zero<source_type>{} && is_zero<flux_type>{}) {
return;
Expand Down
Loading