Skip to content

Debug meshFields backend to fit 1.0.0 API #346

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

Open
Sichao25 wants to merge 1 commit into
base: develop
Choose a base branch
from
Open

Debug meshFields backend to fit 1.0.0 API #346

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
60 changes: 56 additions & 4 deletions src/pcms/field/evaluator/mesh_fields_backend.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,7 @@

#include <Kokkos_Core.hpp>
#include <MeshField.hpp>
#include <MeshField_Config.hpp>
#include <memory>

#include "pcms/field/layout/mesh_fields.h"
Expand Down Expand Up @@ -52,8 +53,13 @@ class MeshFieldBackendImpl : public MeshFieldBackend<T>
Kokkos::View<LO*> offsets) const override
{
auto self = const_cast<MeshFieldBackendImpl<T, Dim, Order>*>(this);
#if MeshFields_VERSION < 10000
return self->mesh_field_.triangleLocalPointEval(localCoords, offsets,
shape_field_);
#else
return self->mesh_field_.triangleLocalPointEval(localCoords, offsets,
shape_field_.field);
#endif
}

void SetData(Rank1View<const T, DeviceMemorySpace> data, size_t num_nodes,
Expand All @@ -65,7 +71,11 @@ class MeshFieldBackendImpl : public MeshFieldBackend<T>
mesh_.nents(dim), KOKKOS_CLASS_LAMBDA(size_t ent) {
for (size_t n = 0; n < num_nodes; ++n) {
for (size_t c = 0; c < num_components; ++c) {
#if MeshFields_VERSION < 10000
shape_field_(ent, n, c, topo) =
#else
shape_field_.field(ent, n, c, topo) =
#endif
data[ent * stride + n * num_components + c];
}
}
Expand All @@ -82,7 +92,11 @@ class MeshFieldBackendImpl : public MeshFieldBackend<T>
for (size_t n = 0; n < num_nodes; ++n) {
for (size_t c = 0; c < num_components; ++c) {
data[ent * stride + n * num_components + c] =
#if MeshFields_VERSION < 10000
shape_field_(ent, n, c, topo);
#else
shape_field_.field(ent, n, c, topo);
#endif
}
}
});
Expand All @@ -91,9 +105,8 @@ class MeshFieldBackendImpl : public MeshFieldBackend<T>
private:
Omega_h::Mesh& mesh_;
MeshField::OmegahMeshField<DefaultExecutionSpace, Dim> mesh_field_;
using ShapeField =
decltype(mesh_field_.template CreateLagrangeField<T, Order, 1>());
ShapeField shape_field_;
using FWC = decltype(mesh_field_.template CreateLagrangeField<T, Order, 1>());
FWC shape_field_; // FWC = FieldWithController; keeps ctrlr alive
};

// ---------------------------------------------------------------------------
Expand Down Expand Up @@ -240,7 +253,11 @@ struct FillCoordinatesAndIndicesFunctor
const auto owner_idx = owning_elem_ids_(i);
LO count = Kokkos::atomic_sub_fetch(&elem_counts_(owner_idx), 1);
LO index = offsets_(owner_idx) + count - 1;
#if MeshFields_VERSION < 10000
for (int j = 0; j < (dim_ + 1); ++j) {
#else
for (int j = 0; j < dim_; ++j) {
#endif
coordinates_(index, j) = coord[j];
}
indices_(index) = i;
Expand Down Expand Up @@ -404,8 +421,18 @@ struct FillCoordinatesDeviceFunctor
global_coords_(orig_idx, 1)};
const auto local =
Omega_h::barycentric_from_global<2, 2>(point, vertex_coords);
#if MeshFields_VERSION < 10000
for (int j = 0; j < (dim_ + 1); ++j)
coordinates_(index, j) = local[j];
#else
// Computation in barycentric_from_global sometimes produce values
// slightly outside [0, 1] for points on or near element boundaries.
for (int j = 0; j < dim_; ++j) {
coordinates_(index, j) = local[j] < 0.0 ? 0.0
: local[j] > 1.0 ? 1.0
: local[j];
}
#endif
indices_(index) = orig_idx;
}
};
Expand Down Expand Up @@ -494,8 +521,13 @@ struct MeshFieldsAdapter2LocalizationHint
}

offsets_ = Kokkos::View<LO*, HostMemorySpace>("offsets", mesh.nelems() + 1);
#if MeshFields_VERSION < 10000
coordinates_ = Kokkos::View<Real**, HostMemorySpace>(
"coordinates", num_valid_, mesh.dim() + 1);
#else
coordinates_ = Kokkos::View<Real**, HostMemorySpace>(
"coordinates", num_valid_, mesh.dim());
#endif
indices_ = Kokkos::View<LO*, HostMemorySpace>("indices", num_valid_);

if (num_missing_ > 0) {
Expand Down Expand Up @@ -537,8 +569,18 @@ struct MeshFieldsAdapter2LocalizationHint
global_coords(orig_idx, 1)};
const auto local =
Omega_h::barycentric_from_global<2, 2>(point, vertex_coords);
#if MeshFields_VERSION < 10000
for (int j = 0; j < (mesh.dim() + 1); ++j)
coordinates_(index, j) = local[j];
#else
// Computation in barycentric_from_global sometimes produce values
// slightly outside [0, 1] for points on or near element boundaries.
for (int j = 0; j < mesh.dim(); ++j) {
coordinates_(index, j) = local[j] < 0.0 ? 0.0
: local[j] > 1.0 ? 1.0
: local[j];
}
#endif
indices_(index) = static_cast<LO>(orig_idx);
}

Expand Down Expand Up @@ -631,8 +673,13 @@ struct MeshFieldsAdapter2LocalizationHint
set_final_functor);

// Step 6: Fill coordinates and indices on device
#if MeshFields_VERSION < 10000
coordinates_d_ =
Kokkos::View<Real**>("coordinates_d", num_valid_, mesh.dim() + 1);
#else
coordinates_d_ =
Kokkos::View<Real**>("coordinates_d", num_valid_, mesh.dim());
#endif
indices_d_ = Kokkos::View<LO*>("indices_d", num_valid_);
const auto tris2verts = mesh.ask_elem_verts();
const auto mesh_coords = mesh.coords();
Expand All @@ -653,8 +700,13 @@ struct MeshFieldsAdapter2LocalizationHint

// Create host mirrors for compatibility (lazy copy - only if needed)
offsets_ = Kokkos::create_mirror_view(offsets_d_);
#if MeshFields_VERSION < 10000
coordinates_ = Kokkos::View<Real**, HostMemorySpace>(
"coordinates_", num_valid_, mesh.dim() + 1);
#else
coordinates_ = Kokkos::View<Real**, HostMemorySpace>(
"coordinates_", num_valid_, mesh.dim());
#endif
DeepCopyMismatchLayouts(coordinates_, coordinates_d_);
indices_ = Kokkos::create_mirror_view(indices_d_);
if (num_missing_ > 0) {
Expand All @@ -681,4 +733,4 @@ struct MeshFieldsAdapter2LocalizationHint

} // namespace pcms

#endif // PCMS_ADAPTER_MESHFIELDS_MESH_FIELDS_BACKEND_H
#endif // PCMS_ADAPTER_MESHFIELDS_MESH_FIELDS_BACKEND_H
11 changes: 10 additions & 1 deletion src/pcms/transfer/mass_matrix_integrator.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,7 @@

#include <KokkosController.hpp>
#include <MeshField.hpp>
#include <MeshField_Config.hpp>
#include <MeshField_Element.hpp>
#include <MeshField_Fail.hpp>
#include <MeshField_For.hpp>
Expand Down Expand Up @@ -38,7 +39,11 @@ class MassMatrixIntegrator : public MeshField::Integrator
// std::cerr << " Number points per Elem : " << numPtsPerElem << "\n";
assert(numPtsPerElem >= 1);
const size_t ptDim = p.extent(1);
#if MeshFields_VERSION < 10000
assert(ptDim == fe.MeshEntDim + 1);
#else
assert(ptDim == fe.MeshEntDim);
#endif
// Copy values needed in the kernel to avoid capturing host references
// (mesh and fe are host objects and cannot be dereferenced on the device)
const auto numElems = mesh.nelems();
Expand All @@ -50,8 +55,12 @@ class MassMatrixIntegrator : public MeshField::Integrator
const auto first = elm * numPtsPerElem;
const auto last = first + numPtsPerElem;
for (auto pt = first; pt < last; pt++) {
// FIXME better way to fill? pass kokkos::subview to getValues?
// FIXME better way to fill? pass kokkos::subview to getValues?
#if MeshFields_VERSION < 10000
Kokkos::Array<MeshField::Real, FieldElement::MeshEntDim + 1>
#else
Kokkos::Array<MeshField::Real, FieldElement::MeshEntDim>
#endif
localCoord;
for (auto i = 0; i < localCoord.size(); i++) {
localCoord[i] = p(pt, i);
Expand Down
24 changes: 23 additions & 1 deletion src/pcms/transfer/omega_h_form_integrator_utils.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,7 @@
#include "pcms/utility/assert.h"
#include <MeshField_Integrate.hpp>
#include <MeshField_Shape.hpp>
#include <MeshField_Config.hpp>
#include <Omega_h_shape.hpp>
#include <Kokkos_MathematicalFunctions.hpp>

Expand Down Expand Up @@ -72,14 +73,27 @@ inline void CheckOmegaHScalarLagrangeLayout(
}

[[nodiscard]] OMEGA_H_INLINE Omega_h::Vector<2> GlobalFromBarycentric(
#if MeshFields_VERSION < 10000
const MeshField::Vector3& barycentric_coord,
#else
const MeshField::Vector2& barycentric_coord,
#endif
const Omega_h::Few<Omega_h::Vector<2>, 3>& verts_coord)
{
Omega_h::Vector<2> real_coords = {0.0, 0.0};
#if MeshFields_VERSION < 10000
for (int i = 0; i < 3; ++i) {
real_coords[0] += barycentric_coord[i] * verts_coord[i][0];
real_coords[1] += barycentric_coord[i] * verts_coord[i][1];
}
#else
const Omega_h::Real xi3 = 1.0 - barycentric_coord[0] - barycentric_coord[1];
const Omega_h::Real xi[3] = {barycentric_coord[0], barycentric_coord[1], xi3};
for (int i = 0; i < 3; ++i) {
real_coords[0] += xi[i] * verts_coord[i][0];
real_coords[1] += xi[i] * verts_coord[i][1];
}
#endif
return real_coords;
}

Expand Down Expand Up @@ -250,15 +264,23 @@ OMEGA_H_INLINE void ForEachIntersectionSubtriangle(
// source_order + target_order), which are only known at construction.
struct IntegrationData
{
#if MeshFields_VERSION < 10000
Kokkos::View<MeshField::Vector3*> bary_coords; // barycentric coordinates
Kokkos::View<Omega_h::Real*> weights; // quadrature weights
#else
Kokkos::View<MeshField::Vector2*> bary_coords; // barycentric coordinates
#endif
Kokkos::View<Omega_h::Real*> weights; // quadrature weights

explicit IntegrationData(int order)
{
auto ip_vec = MeshField::getIntegrationPoints<MeshField::Triangle>(order);
const std::size_t num_ip = ip_vec.size();

#if MeshFields_VERSION < 10000
bary_coords = Kokkos::View<MeshField::Vector3*>("bary_coords", num_ip);
#else
bary_coords = Kokkos::View<MeshField::Vector2*>("bary_coords", num_ip);
#endif
weights = Kokkos::View<Omega_h::Real*>("weights", num_ip);

auto bary_coords_host = Kokkos::create_mirror_view(bary_coords);
Expand Down
4 changes: 4 additions & 0 deletions src/pcms/transfer/omega_h_mass_integrator.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -123,7 +123,11 @@ OmegaHMassIntegrator::OmegaHMassIntegrator(
omf(mesh);
auto coordField = omf.getCoordField();
const auto [shp, map] = MeshField::Omegah::getTriangleElement<1>(mesh);
#if MeshFields_VERSION < 10000
MeshField::FieldElement coordFe(mesh.nelems(), coordField, shp, map);
#else
MeshField::FieldElement coordFe(mesh.nelems(), coordField.field, shp, map);
#endif
auto elm_mass_dev = buildElementMassMatrix(mesh, coordFe);

// Build COO sparsity pattern on device: each element contributes a 3x3 block.
Expand Down
5 changes: 5 additions & 0 deletions test/test_omega_h_mass_integrator.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,7 @@
#include <pcms/transfer/omega_h_mass_integrator.hpp>
#include <KokkosController.hpp>
#include <MeshField.hpp>
#include <MeshField_Config.hpp>
#include <MeshField_Element.hpp>
#include <petscksp.h>
#include <map>
Expand Down Expand Up @@ -60,7 +61,11 @@ std::map<std::pair<pcms::GO, pcms::GO>, pcms::Real> BuildReferenceMassMap(
omf(mesh);
auto coordField = omf.getCoordField();
const auto [shp, map] = MeshField::Omegah::getTriangleElement<1>(mesh);
#if MeshFields_VERSION < 10000
MeshField::FieldElement coordFe(mesh.nelems(), coordField, shp, map);
#else
MeshField::FieldElement coordFe(mesh.nelems(), coordField.field, shp, map);
#endif
auto elm_mass_dev = buildElementMassMatrix(mesh, coordFe);
auto elm_mass_host =
Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace{}, elm_mass_dev);
Expand Down
Loading