const bool use_coloring;
- Number *inv_jac;
- // We would like to use
+ // FIXME We would like to use
+ // Kokkos::Subview<Kokkos::View<Number **[dim][dim],
+ // MemorySpace::Default::kokkos_space>, int, decltype(Kokkos::ALL),
+ // decltype(Kokkos::ALL), decltype(Kokkos::ALL)> but we get error:
+ // incomplete type is not allowed. I cannot reproduce outside of deal.II.
+ // Need to investigate more.
+ Kokkos::Subview<
+ Kokkos::View<Number **[dim][dim], MemorySpace::Default::kokkos_space>,
+ int,
+ Kokkos::pair<int, int>,
+ Kokkos::pair<int, int>,
+ Kokkos::pair<int, int>>
+ inv_jac;
+ // FIXME We would like to use
// Kokkos::Subview<Kokkos::View<Number **,
// MemorySpace::Default::kokkos_space>, int, decltype(Kokkos::ALL)>
// but we get error: incomplete type is not allowed. I cannot reproduce
, mf_object_id(data->id)
, constraint_mask(data->constraint_mask[cell_id])
, use_coloring(data->use_coloring)
+ , inv_jac(Kokkos::subview(
+ data->inv_jacobian,
+ cell_id,
+ Kokkos::pair<int, int>(0, Utilities::pow(n_q_points_1d, dim)),
+ Kokkos::pair<int, int>(0, dim),
+ Kokkos::pair<int, int>(0, dim)))
, JxW(Kokkos::subview(
data->JxW,
cell_id,
, values(shdata->values)
{
local_to_global = data->local_to_global + padding_length * cell_id;
- inv_jac = data->inv_jacobian + padding_length * cell_id;
for (unsigned int i = 0; i < dim; ++i)
gradients[i] = shdata->gradients[i];
// TODO optimize if the mesh is uniform
const unsigned int q_point = internal::compute_index<dim, n_q_points_1d>();
- const Number * inv_jacobian = &inv_jac[q_point];
gradient_type grad;
for (unsigned int d_1 = 0; d_1 < dim; ++d_1)
{
Number tmp = 0.;
for (unsigned int d_2 = 0; d_2 < dim; ++d_2)
- tmp += inv_jacobian[padding_length * n_cells * (dim * d_2 + d_1)] *
- gradients[d_2][q_point];
+ tmp += inv_jac(q_point, d_2, d_1) * gradients[d_2][q_point];
grad[d_1] = tmp;
}
{
// TODO optimize if the mesh is uniform
const unsigned int q_point = internal::compute_index<dim, n_q_points_1d>();
- const Number * inv_jacobian = &inv_jac[q_point];
for (unsigned int d_1 = 0; d_1 < dim; ++d_1)
{
Number tmp = 0.;
for (unsigned int d_2 = 0; d_2 < dim; ++d_2)
- tmp += inv_jacobian[n_cells * padding_length * (dim * d_1 + d_2)] *
- grad_in[d_2];
+ tmp += inv_jac(q_point, d_1, d_2) * grad_in[d_2];
gradients[d_1][q_point] = tmp * JxW[q_point];
}
}
types::global_dof_index *local_to_global;
/**
- * Pointer to the inverse Jacobian.
+ * Kokkos::View of the inverse Jacobian.
*/
- Number *inv_jacobian;
+ Kokkos::View<Number **[dim][dim], MemorySpace::Default::kokkos_space>
+ inv_jacobian;
/**
* Kokkos::View of the Jacobian times the weights.
std::vector<types::global_dof_index *> local_to_global;
/**
- * Vector of pointer to the inverse Jacobian associated to the cells of each
- * color.
+ * Vector of Kokkos::View of the inverse Jacobian associated to the cells of
+ * each color.
*/
- std::vector<Number *> inv_jacobian;
+ std::vector<
+ Kokkos::View<Number **[dim][dim], MemorySpace::Default::kokkos_space>>
+ inv_jacobian;
/**
* Vector of Kokkos::View to the Jacobian times the weights associated to
if (update_flags & update_gradients)
{
data_host.inv_jacobian.resize(n_elements * dim * dim);
- Utilities::CUDA::copy_to_host(data.inv_jacobian,
+ Utilities::CUDA::copy_to_host(data.inv_jacobian.data(),
data_host.inv_jacobian);
}
std::vector<types::global_dof_index> local_to_global_host;
std::vector<Point<dim, Number>> q_points_host;
Kokkos::View<Number **, MemorySpace::Default::kokkos_space> JxW;
- std::vector<Number> inv_jacobian_host;
+ Kokkos::View<Number **[dim][dim], MemorySpace::Default::kokkos_space>
+ inv_jacobian;
std::vector<dealii::internal::MatrixFreeFunctions::ConstraintKinds>
constraint_mask_host;
// Local buffer
dofs_per_cell);
if (update_flags & update_gradients)
- inv_jacobian_host.resize(n_cells * padding_length * dim * dim);
+ inv_jacobian =
+ Kokkos::View<Number **[dim][dim], MemorySpace::Default::kokkos_space>(
+ Kokkos::view_alloc("inv_jacobian_" + std::to_string(color),
+ Kokkos::WithoutInitializing),
+ n_cells,
+ dofs_per_cell);
constraint_mask_host.resize(n_cells);
}
if (update_flags & update_gradients)
{
+ // FIXME too many deep_copy
+ auto inv_jacobian_host = Kokkos::create_mirror_view_and_copy(
+ MemorySpace::Host::kokkos_space{}, inv_jacobian);
const std::vector<DerivativeForm<1, dim, dim>> &inv_jacobians =
fe_values.get_inverse_jacobians();
- std::copy(&inv_jacobians[0][0][0],
- &inv_jacobians[0][0][0] +
- q_points_per_cell * sizeof(DerivativeForm<1, dim, dim>) /
- sizeof(double),
- &inv_jacobian_host[cell_id * padding_length * dim * dim]);
+ for (unsigned int i = 0; i < q_points_per_cell; ++i)
+ for (unsigned int j = 0; j < dim; ++j)
+ for (unsigned int k = 0; k < dim; ++k)
+ inv_jacobian_host(cell_id, i, j, k) = inv_jacobians[i][j][k];
+ Kokkos::deep_copy(inv_jacobian, inv_jacobian_host);
}
}
// Inverse jacobians
if (update_flags & update_gradients)
{
- // Reorder so that all J_11 elements are together, all J_12 elements
- // are together, etc., i.e., reorder indices from
- // cell_id*q_points_per_cell*dim*dim + q*dim*dim +i to
- // i*q_points_per_cell*n_cells + cell_id*q_points_per_cell+q
- transpose_in_place(inv_jacobian_host,
- padding_length * n_cells,
- dim * dim);
-
- alloc_and_copy(&data->inv_jacobian[color],
- ArrayView<const Number>(inv_jacobian_host.data(),
- inv_jacobian_host.size()),
- n_cells * dim * dim * padding_length);
+ data->inv_jacobian[color] = inv_jacobian;
}
alloc_and_copy(
Utilities::CUDA::free(local_to_global_color_ptr);
local_to_global.clear();
- for (auto &inv_jacobian_color_ptr : inv_jacobian)
- Utilities::CUDA::free(inv_jacobian_color_ptr);
- inv_jacobian.clear();
-
for (auto &constraint_mask_color_ptr : constraint_mask)
Utilities::CUDA::free(constraint_mask_color_ptr);
constraint_mask.clear();