From fc073ac178cc2d7a7a4e2353b385038d9ed35b07 Mon Sep 17 00:00:00 2001 From: kloepj2 Date: Wed, 29 Jan 2025 13:01:22 -0500 Subject: [PATCH 01/21] put cabana controller OmegahElement test --- src/MeshField.hpp | 3 +-- src/MeshField_ShapeField.hpp | 7 +++---- test/testOmegahElement.cpp | 6 +++--- 3 files changed, 7 insertions(+), 9 deletions(-) diff --git a/src/MeshField.hpp b/src/MeshField.hpp index 4c00df8..2f16335 100644 --- a/src/MeshField.hpp +++ b/src/MeshField.hpp @@ -163,8 +163,7 @@ template auto getTriangleElement(Omega_h::Mesh &mesh) { namespace MeshField { -template typename Controller = - MeshField::KokkosController> +template typename Controller> class OmegahMeshField { private: Omega_h::Mesh &mesh; diff --git a/src/MeshField_ShapeField.hpp b/src/MeshField_ShapeField.hpp index ca3a448..aebeffe 100644 --- a/src/MeshField_ShapeField.hpp +++ b/src/MeshField_ShapeField.hpp @@ -175,7 +175,7 @@ auto CreateLagrangeField(const MeshInfo &meshInfo) { } using Ctrlr = Controller; // 1 dof with 1 component per vtx - Ctrlr kk_ctrl({/*field 0*/ 1, 1, meshInfo.numVtx}); + Ctrlr kk_ctrl; auto vtxField = MeshField::makeField(kk_ctrl); using LA = LinearAccessor; using LinearLagrangeShapeField = ShapeField; @@ -191,8 +191,7 @@ auto CreateLagrangeField(const MeshInfo &meshInfo) { using Ctrlr = Controller; // 1 dof with 1 comp per vtx/edge - Ctrlr kk_ctrl({/*field 0*/ 1, 1, meshInfo.numVtx, - /*field 1*/ 1, 1, meshInfo.numEdge}); + Ctrlr kk_ctrl; auto vtxField = MeshField::makeField(kk_ctrl); auto edgeField = MeshField::makeField(kk_ctrl); using QA = QuadraticAccessor; @@ -234,7 +233,7 @@ auto CreateCoordinateField(const MeshInfo &meshInfo) { using MemorySpace = typename ExecutionSpace::memory_space; using Ctrlr = Controller; const int numComp = meshInfo.dim; - Ctrlr kk_ctrl({/*field 0*/ 1, numComp, meshInfo.numVtx}); + Ctrlr kk_ctrl; auto vtxField = MeshField::makeField(kk_ctrl); using LA = LinearAccessor; using LinearLagrangeShapeField = ShapeField; diff --git a/test/testOmegahElement.cpp b/test/testOmegahElement.cpp index b86bae5..4ff07e1 100644 --- a/test/testOmegahElement.cpp +++ b/test/testOmegahElement.cpp @@ -144,9 +144,9 @@ int main(int argc, char **argv) { MeshField::Debug = true; { auto mesh = createMeshTri18(lib); - MeshField::OmegahMeshField omf( + MeshField::OmegahMeshField omf( mesh); - +/* // setup field with values from the analytic function static const size_t OnePtPerElem = 1; static const size_t ThreePtsPerElem = 3; @@ -223,7 +223,7 @@ int main(int argc, char **argv) { if (failed) doFail("quadratic", "linear", testCase.name); } - } + }*/ } Kokkos::finalize(); return 0; From 666b17b0c1d2d391d213ecab446c65a9233e0098 Mon Sep 17 00:00:00 2001 From: kloepj2 Date: Sun, 2 Feb 2025 22:48:23 -0500 Subject: [PATCH 02/21] [skip ci] Added cabana controller to shape field --- src/MeshField_ShapeField.hpp | 18 +++++++++++++++--- test/testOmegahElement.cpp | 9 ++++++++- 2 files changed, 23 insertions(+), 4 deletions(-) diff --git a/src/MeshField_ShapeField.hpp b/src/MeshField_ShapeField.hpp index bf8d617..19f64c4 100644 --- a/src/MeshField_ShapeField.hpp +++ b/src/MeshField_ShapeField.hpp @@ -2,6 +2,9 @@ #define MESHFIELD_SHAPEFIELD_HPP #include "KokkosController.hpp" +#ifdef MESHFIELDS_ENABLE_CABANA +#include "CabanaController.hpp" +#endif #include "MeshField_Field.hpp" #include "MeshField_Shape.hpp" #include //decltype @@ -221,6 +224,8 @@ auto CreateLagrangeField(const MeshInfo &meshInfo) { * @param meshInfo defines on-process mesh metadata * @return a linear ShapeField */ + + template typename Controller = MeshField::KokkosController> auto CreateCoordinateField(const MeshInfo &meshInfo) { @@ -229,15 +234,22 @@ auto CreateCoordinateField(const MeshInfo &meshInfo) { } using DataType = Real; using MemorySpace = typename ExecutionSpace::memory_space; - using Ctrlr = Controller; const int numComp = meshInfo.dim; - Ctrlr kk_ctrl({/*field 0*/ 1, numComp, meshInfo.numVtx}); + using Ctrlr = std::conditional_t, MeshField::CabanaController>, + Controller, + Controller>; + auto createController = [](const int numComp, auto numVtx) + {if constexpr (std::is_same_v, MeshField::CabanaController>) + {return Ctrlr(numVtx);} + else + {return Ctrlr({/*field 0*/ 1, numComp, numVtx});}}; + Ctrlr kk_ctrl = createController(numComp, meshInfo.numVtx); auto vtxField = MeshField::makeField(kk_ctrl); using LA = LinearAccessor; using LinearLagrangeShapeField = ShapeField; LinearLagrangeShapeField llsf(kk_ctrl, meshInfo, {vtxField}); return llsf; -}; + }; } // namespace MeshField diff --git a/test/testOmegahElement.cpp b/test/testOmegahElement.cpp index b86bae5..217f37a 100644 --- a/test/testOmegahElement.cpp +++ b/test/testOmegahElement.cpp @@ -4,6 +4,9 @@ #include "MeshField_Fail.hpp" //remove? #include "MeshField_For.hpp" //remove? #include "MeshField_ShapeField.hpp" //remove? +#ifdef MESHFIELDS_ENABLE_CABANA +#include "CabanaController.hpp" +#endif #include "Omega_h_build.hpp" #include "Omega_h_file.hpp" #include "Omega_h_simplex.hpp" @@ -140,8 +143,12 @@ void doFail(std::string_view order, std::string_view function, int main(int argc, char **argv) { Kokkos::initialize(argc, argv); - auto lib = Omega_h::Library(&argc, &argv); + auto lib = Omega_h::Library(&argc, &argv); MeshField::Debug = true; + { + auto mesh = createMeshTri18(lib); + MeshField::OmegahMeshField omf(mesh); + } { auto mesh = createMeshTri18(lib); MeshField::OmegahMeshField omf( From 289148d819ba3e2a9b11cd495ff27d69578a07bf Mon Sep 17 00:00:00 2001 From: kloepj2 Date: Tue, 4 Feb 2025 18:27:47 -0500 Subject: [PATCH 03/21] added cabana controller initalizations where needed in shape field. Worked on copying over kokkos tests to test cabana --- src/MeshField_ShapeField.hpp | 26 ++++++++++++++++++------ test/testOmegahElement.cpp | 38 ++++++++++++++++++++++++++++++++++++ 2 files changed, 58 insertions(+), 6 deletions(-) diff --git a/src/MeshField_ShapeField.hpp b/src/MeshField_ShapeField.hpp index 3081477..059a30e 100644 --- a/src/MeshField_ShapeField.hpp +++ b/src/MeshField_ShapeField.hpp @@ -173,9 +173,16 @@ auto CreateLagrangeField(const MeshInfo &meshInfo) { if (meshInfo.numVtx <= 0) { fail("mesh has no vertices\n"); } - using Ctrlr = Controller; + using Ctrlr = std::conditional_t, MeshField::CabanaController>, + Controller, + Controller>; // 1 dof with 1 component per vtx - Ctrlr kk_ctrl; + auto createController = [](const int numComp, auto numVtx) + {if constexpr (std::is_same_v, MeshField::CabanaController>) + {return Ctrlr(numVtx);} + else + {return Ctrlr({/*field 0*/ 1, numComp, numVtx});}}; + Ctrlr kk_ctrl = createController(1, meshInfo.numVtx); auto vtxField = MeshField::makeField(kk_ctrl); using LA = LinearAccessor; using LinearLagrangeShapeField = ShapeField; @@ -188,10 +195,17 @@ auto CreateLagrangeField(const MeshInfo &meshInfo) { if (meshInfo.numEdge <= 0) { fail("mesh has no edges\n"); } - using Ctrlr = - Controller; + using Ctrlr = std::conditional_t, MeshField::CabanaController>, + Controller, + Controller>; // 1 dof with 1 comp per vtx/edge - Ctrlr kk_ctrl; + auto createController = [](const int numComp, auto numVtx, auto numEdge) + {if constexpr (std::is_same_v, MeshField::CabanaController>) + {return Ctrlr(numVtx * numEdge);} + else + {return Ctrlr({/*field 0*/ 1, numComp, numVtx, + /*field 1*/ 1, numComp, numEdge});}}; + Ctrlr kk_ctrl = createController(1, meshInfo.numVtx, meshInfo.numEdge); auto vtxField = MeshField::makeField(kk_ctrl); auto edgeField = MeshField::makeField(kk_ctrl); using QA = QuadraticAccessor; @@ -235,7 +249,7 @@ auto CreateCoordinateField(const MeshInfo &meshInfo) { using MemorySpace = typename ExecutionSpace::memory_space; const int numComp = meshInfo.dim; using Ctrlr = std::conditional_t, MeshField::CabanaController>, - Controller, + Controller, Controller>; auto createController = [](const int numComp, auto numVtx) {if constexpr (std::is_same_v, MeshField::CabanaController>) diff --git a/test/testOmegahElement.cpp b/test/testOmegahElement.cpp index 217f37a..a2fc175 100644 --- a/test/testOmegahElement.cpp +++ b/test/testOmegahElement.cpp @@ -148,6 +148,44 @@ int main(int argc, char **argv) { { auto mesh = createMeshTri18(lib); MeshField::OmegahMeshField omf(mesh); + + // setup field with values from the analytic function + static const size_t OnePtPerElem = 1; + static const size_t ThreePtsPerElem = 3; + auto centroids = createElmAreaCoords( + mesh.nfaces(), {1 / 3.0, 1 / 3.0, 1 / 3.0}); + auto interior = + createElmAreaCoords(mesh.nfaces(), {0.1, 0.4, 0.5}); + auto vertex = + createElmAreaCoords(mesh.nfaces(), {0.0, 0.0, 1.0}); + // clang-format off + auto allVertices = createElmAreaCoords(mesh.nfaces(), + {1.0, 0.0, 0.0, + 0.0, 1.0, 0.0, + 0.0, 0.0, 1.0}); + const auto cases = {TestCoords{centroids, OnePtPerElem, "centroids"}, + TestCoords{interior, OnePtPerElem, "interior"}, + TestCoords{vertex, OnePtPerElem, "vertex"}, + TestCoords{allVertices, ThreePtsPerElem, "allVertices"}}; + // clang-format on + + auto coords = mesh.coords(); + const auto MeshDim = 2; + for (auto testCase : cases) { + using ViewType = decltype(testCase.coords); + { + const auto ShapeOrder = 1; + auto field = omf.CreateLagrangeField(); + auto func = LinearFunction(); + setVertices(mesh, func, field); + using FieldType = decltype(field); + //auto result = omf.triangleLocalPointEval( + // testCase.coords, testCase.NumPtsPerElem, LinearFunction{}, + // field); + + + } + } } { auto mesh = createMeshTri18(lib); From ddd3b79b761d99c16966dd36cbb768cd95daa6b9 Mon Sep 17 00:00:00 2001 From: kloepj2 Date: Tue, 18 Feb 2025 23:28:12 -0500 Subject: [PATCH 04/21] completed adding cabana tests. Adjusted ordering. --- src/MeshField.hpp | 14 ++++----- src/MeshField_Element.hpp | 13 +++++++-- src/MeshField_ShapeField.hpp | 16 +++++----- test/testOmegahCoordField.cpp | 4 +-- test/testOmegahElement.cpp | 55 ++++++++++++++++++++++++++++++----- 5 files changed, 74 insertions(+), 28 deletions(-) diff --git a/src/MeshField.hpp b/src/MeshField.hpp index 2f16335..38ff887 100644 --- a/src/MeshField.hpp +++ b/src/MeshField.hpp @@ -40,13 +40,13 @@ createCoordinateField(MeshField::MeshInfo mesh_info, Omega_h::Reals coords) { const auto meshDim = mesh_info.dim; auto coordField = MeshField::CreateCoordinateField(mesh_info); - auto setCoordField = KOKKOS_LAMBDA(const int &i) { - coordField(0, 0, i, MeshField::Vertex) = coords[i * meshDim]; - coordField(0, 1, i, MeshField::Vertex) = coords[i * meshDim + 1]; - }; - MeshField::parallel_for(ExecutionSpace(), {0}, {mesh_info.numVtx}, - setCoordField, "setCoordField"); - return coordField; + auto setCoordField = KOKKOS_LAMBDA(const int &i) { + coordField(i, 0, 0, MeshField::Vertex) = coords[i * meshDim]; + coordField(i, 0, 1, MeshField::Vertex) = coords[i * meshDim + 1]; + }; + MeshField::parallel_for(ExecutionSpace(), {0}, {mesh_info.numVtx}, + setCoordField, "setCoordField"); + return coordField; } struct LinearTriangleToVertexField { diff --git a/src/MeshField_Element.hpp b/src/MeshField_Element.hpp index bde61f6..3ca9a64 100644 --- a/src/MeshField_Element.hpp +++ b/src/MeshField_Element.hpp @@ -59,8 +59,15 @@ struct FieldElement { const ElementDofHolderAccessor elm2dofIn) : numMeshEnts(numMeshEntsIn), field(fieldIn), shapeFn(shapeFnIn), elm2dof(elm2dofIn) {} - - using ValArray = Kokkos::Array + struct baseType { + using type = T; + }; + template + struct baseType { + using type = typename baseType::type; + }; + using ValArray = Kokkos::Array::type, ShapeType::numComponentsPerDof>; static const size_t NumComponents = ShapeType::numComponentsPerDof; @@ -91,7 +98,7 @@ struct FieldElement { for (int ci = 0; ci < shapeFn.numComponentsPerDof; ++ci) { auto map = elm2dof(ni, ci, ent, topo); const auto fval = - field(map.node, map.component, map.entity, map.topo); + field(map.entity, map.node, map.component, map.topo); c[ci] += fval * shapeValues[ni]; } } diff --git a/src/MeshField_ShapeField.hpp b/src/MeshField_ShapeField.hpp index 059a30e..cc85a12 100644 --- a/src/MeshField_ShapeField.hpp +++ b/src/MeshField_ShapeField.hpp @@ -174,14 +174,14 @@ auto CreateLagrangeField(const MeshInfo &meshInfo) { fail("mesh has no vertices\n"); } using Ctrlr = std::conditional_t, MeshField::CabanaController>, - Controller, + Controller, Controller>; // 1 dof with 1 component per vtx auto createController = [](const int numComp, auto numVtx) {if constexpr (std::is_same_v, MeshField::CabanaController>) {return Ctrlr(numVtx);} else - {return Ctrlr({/*field 0*/ 1, numComp, numVtx});}}; + {return Ctrlr({/*field 0*/ numVtx, 1, numComp});}}; Ctrlr kk_ctrl = createController(1, meshInfo.numVtx); auto vtxField = MeshField::makeField(kk_ctrl); using LA = LinearAccessor; @@ -196,15 +196,15 @@ auto CreateLagrangeField(const MeshInfo &meshInfo) { fail("mesh has no edges\n"); } using Ctrlr = std::conditional_t, MeshField::CabanaController>, - Controller, + Controller, Controller>; // 1 dof with 1 comp per vtx/edge auto createController = [](const int numComp, auto numVtx, auto numEdge) {if constexpr (std::is_same_v, MeshField::CabanaController>) - {return Ctrlr(numVtx * numEdge);} + {return Ctrlr(std::max(numVtx, numEdge));} else - {return Ctrlr({/*field 0*/ 1, numComp, numVtx, - /*field 1*/ 1, numComp, numEdge});}}; + {return Ctrlr({/*field 0*/ numVtx, 1, numComp, + /*field 1*/ numVtx, 1, numComp});}}; Ctrlr kk_ctrl = createController(1, meshInfo.numVtx, meshInfo.numEdge); auto vtxField = MeshField::makeField(kk_ctrl); auto edgeField = MeshField::makeField(kk_ctrl); @@ -249,13 +249,13 @@ auto CreateCoordinateField(const MeshInfo &meshInfo) { using MemorySpace = typename ExecutionSpace::memory_space; const int numComp = meshInfo.dim; using Ctrlr = std::conditional_t, MeshField::CabanaController>, - Controller, + Controller, Controller>; auto createController = [](const int numComp, auto numVtx) {if constexpr (std::is_same_v, MeshField::CabanaController>) {return Ctrlr(numVtx);} else - {return Ctrlr({/*field 0*/ 1, numComp, numVtx});}}; + {return Ctrlr({/*field 0*/ numVtx, 1, numComp});}}; Ctrlr kk_ctrl = createController(numComp, meshInfo.numVtx); auto vtxField = MeshField::makeField(kk_ctrl); using LA = LinearAccessor; diff --git a/test/testOmegahCoordField.cpp b/test/testOmegahCoordField.cpp index 4a8fb6d..79997bf 100644 --- a/test/testOmegahCoordField.cpp +++ b/test/testOmegahCoordField.cpp @@ -79,8 +79,8 @@ bool triangleLocalToGlobal(Omega_h::Mesh mesh) { MeshField::CreateCoordinateField(meshInfo); auto setCoordField = KOKKOS_LAMBDA(const int &i) { - coordField(0, 0, i, MeshField::Vertex) = coords[i * MeshDim]; - coordField(0, 1, i, MeshField::Vertex) = coords[i * MeshDim + 1]; + coordField(i, 0, 0, MeshField::Vertex) = coords[i * MeshDim]; + coordField(i, 0, 1, MeshField::Vertex) = coords[i * MeshDim + 1]; }; MeshField::parallel_for(ExecutionSpace(), {0}, {meshInfo.numVtx}, setCoordField, "setCoordField"); diff --git a/test/testOmegahElement.cpp b/test/testOmegahElement.cpp index a2fc175..82448bc 100644 --- a/test/testOmegahElement.cpp +++ b/test/testOmegahElement.cpp @@ -88,7 +88,7 @@ void setVertices(Omega_h::Mesh &mesh, AnalyticFunction func, ShapeField field) { // - TODO should be encoded in the field? const auto x = coords[vtx * MeshDim]; const auto y = coords[vtx * MeshDim + 1]; - field(0, 0, vtx, MeshField::Vertex) = func(x, y); + field(vtx, 0, 0, MeshField::Vertex) = func(x, y); }; MeshField::parallel_for(ExecutionSpace(), {0}, {mesh.nverts()}, setFieldAtVertices, "setFieldAtVertices"); @@ -109,7 +109,7 @@ void setEdges(Omega_h::Mesh &mesh, AnalyticFunction func, ShapeField field) { const auto x = (coords[left * MeshDim] + coords[right * MeshDim]) / 2.0; const auto y = (coords[left * MeshDim + 1] + coords[right * MeshDim + 1]) / 2.0; - field(0, 0, edge, MeshField::Edge) = func(x, y); + field(edge, 0, 0, MeshField::Edge) = func(x, y); }; MeshField::parallel_for(ExecutionSpace(), {0}, {mesh.nedges()}, setFieldAtEdges, "setFieldAtEdges"); @@ -177,14 +177,53 @@ int main(int argc, char **argv) { const auto ShapeOrder = 1; auto field = omf.CreateLagrangeField(); auto func = LinearFunction(); - setVertices(mesh, func, field); - using FieldType = decltype(field); - //auto result = omf.triangleLocalPointEval( - // testCase.coords, testCase.NumPtsPerElem, LinearFunction{}, - // field); + setVertices(mesh, func, field); + using FieldType = decltype(field); + auto result = omf.triangleLocalPointEval( + testCase.coords, testCase.NumPtsPerElem, LinearFunction{}, + field); + auto failed = checkResult(mesh, result, omf.getCoordField(), testCase, + LinearFunction{}); + if (failed) + doFail("linear", "linear", testCase.name); + } + { + const auto ShapeOrder = 2; + auto field = + omf.CreateLagrangeField(); + auto func = QuadraticFunction(); + setVertices(mesh, func, field); + setEdges(mesh, func, field); + using FieldType = decltype(field); + auto result = + omf.triangleLocalPointEval( + testCase.coords, testCase.NumPtsPerElem, QuadraticFunction{}, + field); + auto failed = checkResult(mesh, result, omf.getCoordField(), testCase, + QuadraticFunction{}); + if (failed) + doFail("quadratic", "quadratic", testCase.name); + } +{ + const auto ShapeOrder = 2; + auto field = + omf.CreateLagrangeField(); + auto func = LinearFunction(); + auto coords = mesh.coords(); + setVertices(mesh, func, field); + setEdges(mesh, func, field); + using FieldType = decltype(field); + auto result = + omf.triangleLocalPointEval( + testCase.coords, testCase.NumPtsPerElem, LinearFunction{}, + field); + auto failed = checkResult(mesh, result, omf.getCoordField(), testCase, + LinearFunction{}); + if (failed) + doFail("quadratic", "linear", testCase.name); + } - } } } { From a243d046cf70e04ed21ecc0a218d314258677581 Mon Sep 17 00:00:00 2001 From: kloepj2 Date: Tue, 18 Feb 2025 23:33:42 -0500 Subject: [PATCH 05/21] reverted previous change --- src/MeshField.hpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/MeshField.hpp b/src/MeshField.hpp index 38ff887..892601c 100644 --- a/src/MeshField.hpp +++ b/src/MeshField.hpp @@ -163,7 +163,8 @@ template auto getTriangleElement(Omega_h::Mesh &mesh) { namespace MeshField { -template typename Controller> +template typename Controller = + MeshField::KokkosController> class OmegahMeshField { private: Omega_h::Mesh &mesh; From 483f39b22abb01a37450d59f28d27910d6090ca9 Mon Sep 17 00:00:00 2001 From: kloepj2 Date: Wed, 19 Feb 2025 12:37:41 -0500 Subject: [PATCH 06/21] [skip ci] Returning correct instantiation using the Cabana Controller --- src/MeshField_ShapeField.hpp | 54 ++++++++++++++++++++++++++---------- 1 file changed, 39 insertions(+), 15 deletions(-) diff --git a/src/MeshField_ShapeField.hpp b/src/MeshField_ShapeField.hpp index cc85a12..5d302ce 100644 --- a/src/MeshField_ShapeField.hpp +++ b/src/MeshField_ShapeField.hpp @@ -248,21 +248,45 @@ auto CreateCoordinateField(const MeshInfo &meshInfo) { using DataType = Real; using MemorySpace = typename ExecutionSpace::memory_space; const int numComp = meshInfo.dim; - using Ctrlr = std::conditional_t, MeshField::CabanaController>, - Controller, - Controller>; - auto createController = [](const int numComp, auto numVtx) - {if constexpr (std::is_same_v, MeshField::CabanaController>) - {return Ctrlr(numVtx);} - else - {return Ctrlr({/*field 0*/ numVtx, 1, numComp});}}; - Ctrlr kk_ctrl = createController(numComp, meshInfo.numVtx); - auto vtxField = MeshField::makeField(kk_ctrl); - using LA = LinearAccessor; - using LinearLagrangeShapeField = ShapeField; - LinearLagrangeShapeField llsf(kk_ctrl, meshInfo, {vtxField}); - return llsf; - }; + if constexpr (std::is_same_v, MeshField::CabanaController>) { + if (numComp == 1) { + using Ctrlr = Controller; + Ctrlr kk_ctrl(meshInfo.numVtx); + auto vtxField = MeshField::makeField(kk_ctrl); + using LA = LinearAccessor; + using LinearLagrangeShapeField = ShapeField; + LinearLagrangeShapeField llsf(kk_ctrl, meshInfo, {vtxField}); + return llsf; + } + else if (numComp == 2) { + using Ctrlr = Controller; + Ctrlr kk_ctrl(meshInfo.numVtx); + auto vtxField = MeshField::makeField(kk_ctrl); + using LA = LinearAccessor; + using LinearLagrangeShapeField = ShapeField; + LinearLagrangeShapeField llsf(kk_ctrl, meshInfo, {vtxField}); + return llsf; + } + else { + using Ctrlr = Controller; + Ctrlr kk_ctrl(meshInfo.numVtx); + auto vtxField = MeshField::makeField(kk_ctrl); + using LA = LinearAccessor; + using LinearLagrangeShapeField = ShapeField; + LinearLagrangeShapeField llsf(kk_ctrl, meshInfo, {vtxField}); + return llsf; + } + } + else { + using Ctrlr = Controller; + Ctrlr kk_ctrl({/*field 0*/ meshInfo.numVtx, 1, numComp}); + auto vtxField = MeshField::makeField(kk_ctrl); + using LA = LinearAccessor; + using LinearLagrangeShapeField = ShapeField; + LinearLagrangeShapeField llsf(kk_ctrl, meshInfo, {vtxField}); + return llsf; + } +}; } // namespace MeshField From 6741769af4aefeae050f993100dfaaa3c0ce2a38 Mon Sep 17 00:00:00 2001 From: kloepj2 Date: Thu, 20 Feb 2025 00:34:08 -0500 Subject: [PATCH 07/21] completed final formatting and revisions --- src/MeshField.hpp | 14 ++--- src/MeshField_Element.hpp | 15 +++-- src/MeshField_ShapeField.hpp | 110 +++++++++++++++++------------------ test/testOmegahElement.cpp | 39 +++++++------ 4 files changed, 89 insertions(+), 89 deletions(-) diff --git a/src/MeshField.hpp b/src/MeshField.hpp index 892601c..562522f 100644 --- a/src/MeshField.hpp +++ b/src/MeshField.hpp @@ -40,13 +40,13 @@ createCoordinateField(MeshField::MeshInfo mesh_info, Omega_h::Reals coords) { const auto meshDim = mesh_info.dim; auto coordField = MeshField::CreateCoordinateField(mesh_info); - auto setCoordField = KOKKOS_LAMBDA(const int &i) { - coordField(i, 0, 0, MeshField::Vertex) = coords[i * meshDim]; - coordField(i, 0, 1, MeshField::Vertex) = coords[i * meshDim + 1]; - }; - MeshField::parallel_for(ExecutionSpace(), {0}, {mesh_info.numVtx}, - setCoordField, "setCoordField"); - return coordField; + auto setCoordField = KOKKOS_LAMBDA(const int &i) { + coordField(i, 0, 0, MeshField::Vertex) = coords[i * meshDim]; + coordField(i, 0, 1, MeshField::Vertex) = coords[i * meshDim + 1]; + }; + MeshField::parallel_for(ExecutionSpace(), {0}, {mesh_info.numVtx}, + setCoordField, "setCoordField"); + return coordField; } struct LinearTriangleToVertexField { diff --git a/src/MeshField_Element.hpp b/src/MeshField_Element.hpp index 3ca9a64..c1c2bff 100644 --- a/src/MeshField_Element.hpp +++ b/src/MeshField_Element.hpp @@ -59,16 +59,15 @@ struct FieldElement { const ElementDofHolderAccessor elm2dofIn) : numMeshEnts(numMeshEntsIn), field(fieldIn), shapeFn(shapeFnIn), elm2dof(elm2dofIn) {} - template - struct baseType { - using type = T; + template struct baseType { + using type = T; }; - template - struct baseType { - using type = typename baseType::type; + template struct baseType { + using type = typename baseType::type; }; - using ValArray = Kokkos::Array::type, - ShapeType::numComponentsPerDof>; + using ValArray = + Kokkos::Array::type, + ShapeType::numComponentsPerDof>; static const size_t NumComponents = ShapeType::numComponentsPerDof; /** diff --git a/src/MeshField_ShapeField.hpp b/src/MeshField_ShapeField.hpp index 5d302ce..a1d4786 100644 --- a/src/MeshField_ShapeField.hpp +++ b/src/MeshField_ShapeField.hpp @@ -173,15 +173,23 @@ auto CreateLagrangeField(const MeshInfo &meshInfo) { if (meshInfo.numVtx <= 0) { fail("mesh has no vertices\n"); } - using Ctrlr = std::conditional_t, MeshField::CabanaController>, - Controller, - Controller>; + using Ctrlr = std::conditional_t< + std::is_same_v< + Controller, + MeshField::CabanaController>, + Controller, + Controller>; // 1 dof with 1 component per vtx - auto createController = [](const int numComp, auto numVtx) - {if constexpr (std::is_same_v, MeshField::CabanaController>) - {return Ctrlr(numVtx);} - else - {return Ctrlr({/*field 0*/ numVtx, 1, numComp});}}; + auto createController = [](const int numComp, auto numVtx) { + if constexpr (std::is_same_v< + Controller, + MeshField::CabanaController>) { + return Ctrlr(numVtx); + } else { + return Ctrlr({/*field 0*/ numVtx, 1, numComp}); + } + }; Ctrlr kk_ctrl = createController(1, meshInfo.numVtx); auto vtxField = MeshField::makeField(kk_ctrl); using LA = LinearAccessor; @@ -195,16 +203,24 @@ auto CreateLagrangeField(const MeshInfo &meshInfo) { if (meshInfo.numEdge <= 0) { fail("mesh has no edges\n"); } - using Ctrlr = std::conditional_t, MeshField::CabanaController>, - Controller, - Controller>; + using Ctrlr = std::conditional_t< + std::is_same_v< + Controller, + MeshField::CabanaController>, + Controller, + Controller>; // 1 dof with 1 comp per vtx/edge - auto createController = [](const int numComp, auto numVtx, auto numEdge) - {if constexpr (std::is_same_v, MeshField::CabanaController>) - {return Ctrlr(std::max(numVtx, numEdge));} - else - {return Ctrlr({/*field 0*/ numVtx, 1, numComp, - /*field 1*/ numVtx, 1, numComp});}}; + auto createController = [](const int numComp, auto numVtx, auto numEdge) { + if constexpr (std::is_same_v< + Controller, + MeshField::CabanaController>) { + return Ctrlr(std::max(numVtx, numEdge)); + } else { + return Ctrlr({/*field 0*/ numVtx, 1, numComp, + /*field 1*/ numEdge, 1, numComp}); + } + }; Ctrlr kk_ctrl = createController(1, meshInfo.numVtx, meshInfo.numEdge); auto vtxField = MeshField::makeField(kk_ctrl); auto edgeField = MeshField::makeField(kk_ctrl); @@ -238,7 +254,6 @@ auto CreateLagrangeField(const MeshInfo &meshInfo) { * @return a linear ShapeField */ - template typename Controller = MeshField::KokkosController> auto CreateCoordinateField(const MeshInfo &meshInfo) { @@ -248,44 +263,29 @@ auto CreateCoordinateField(const MeshInfo &meshInfo) { using DataType = Real; using MemorySpace = typename ExecutionSpace::memory_space; const int numComp = meshInfo.dim; - if constexpr (std::is_same_v, MeshField::CabanaController>) { - if (numComp == 1) { - using Ctrlr = Controller; - Ctrlr kk_ctrl(meshInfo.numVtx); - auto vtxField = MeshField::makeField(kk_ctrl); - using LA = LinearAccessor; - using LinearLagrangeShapeField = ShapeField; - LinearLagrangeShapeField llsf(kk_ctrl, meshInfo, {vtxField}); - return llsf; - } - else if (numComp == 2) { - using Ctrlr = Controller; - Ctrlr kk_ctrl(meshInfo.numVtx); - auto vtxField = MeshField::makeField(kk_ctrl); - using LA = LinearAccessor; - using LinearLagrangeShapeField = ShapeField; - LinearLagrangeShapeField llsf(kk_ctrl, meshInfo, {vtxField}); - return llsf; - } - else { - using Ctrlr = Controller; - Ctrlr kk_ctrl(meshInfo.numVtx); - auto vtxField = MeshField::makeField(kk_ctrl); - using LA = LinearAccessor; - using LinearLagrangeShapeField = ShapeField; - LinearLagrangeShapeField llsf(kk_ctrl, meshInfo, {vtxField}); - return llsf; + // FIXME Oversized cabana datatypes when numComp = 1|2 + using Ctrlr = std::conditional_t< + std::is_same_v< + Controller, + MeshField::CabanaController>, + Controller, + Controller>; + auto createController = [](const int numComp, auto numVtx) { + if constexpr (std::is_same_v< + Controller, + MeshField::CabanaController>) { + return Ctrlr(numVtx); + } else { + return Ctrlr({/*field 0*/ numVtx, 1, numComp}); } - } - else { - using Ctrlr = Controller; - Ctrlr kk_ctrl({/*field 0*/ meshInfo.numVtx, 1, numComp}); - auto vtxField = MeshField::makeField(kk_ctrl); - using LA = LinearAccessor; - using LinearLagrangeShapeField = ShapeField; - LinearLagrangeShapeField llsf(kk_ctrl, meshInfo, {vtxField}); - return llsf; - } + }; + Ctrlr kk_ctrl = createController(numComp, meshInfo.numVtx); + auto vtxField = MeshField::makeField(kk_ctrl); + using LA = LinearAccessor; + using LinearLagrangeShapeField = ShapeField; + LinearLagrangeShapeField llsf(kk_ctrl, meshInfo, {vtxField}); + return llsf; }; } // namespace MeshField diff --git a/test/testOmegahElement.cpp b/test/testOmegahElement.cpp index 82448bc..d1309f7 100644 --- a/test/testOmegahElement.cpp +++ b/test/testOmegahElement.cpp @@ -143,11 +143,12 @@ void doFail(std::string_view order, std::string_view function, int main(int argc, char **argv) { Kokkos::initialize(argc, argv); - auto lib = Omega_h::Library(&argc, &argv); + auto lib = Omega_h::Library(&argc, &argv); MeshField::Debug = true; { auto mesh = createMeshTri18(lib); - MeshField::OmegahMeshField omf(mesh); + MeshField::OmegahMeshField omf( + mesh); // setup field with values from the analytic function static const size_t OnePtPerElem = 1; @@ -172,22 +173,24 @@ int main(int argc, char **argv) { auto coords = mesh.coords(); const auto MeshDim = 2; for (auto testCase : cases) { - using ViewType = decltype(testCase.coords); - { - const auto ShapeOrder = 1; - auto field = omf.CreateLagrangeField(); - auto func = LinearFunction(); - setVertices(mesh, func, field); - using FieldType = decltype(field); - auto result = omf.triangleLocalPointEval( - testCase.coords, testCase.NumPtsPerElem, LinearFunction{}, - field); - auto failed = checkResult(mesh, result, omf.getCoordField(), testCase, + using ViewType = decltype(testCase.coords); + { + const auto ShapeOrder = 1; + auto field = + omf.CreateLagrangeField(); + auto func = LinearFunction(); + setVertices(mesh, func, field); + using FieldType = decltype(field); + auto result = + omf.triangleLocalPointEval( + testCase.coords, testCase.NumPtsPerElem, LinearFunction{}, + field); + auto failed = checkResult(mesh, result, omf.getCoordField(), testCase, LinearFunction{}); if (failed) doFail("linear", "linear", testCase.name); - } - { + } + { const auto ShapeOrder = 2; auto field = omf.CreateLagrangeField(); @@ -204,14 +207,14 @@ int main(int argc, char **argv) { if (failed) doFail("quadratic", "quadratic", testCase.name); } -{ + { const auto ShapeOrder = 2; auto field = omf.CreateLagrangeField(); auto func = LinearFunction(); auto coords = mesh.coords(); setVertices(mesh, func, field); - setEdges(mesh, func, field); + setEdges(mesh, func, field); using FieldType = decltype(field); auto result = omf.triangleLocalPointEval( @@ -222,8 +225,6 @@ int main(int argc, char **argv) { if (failed) doFail("quadratic", "linear", testCase.name); } - - } } { From 5b4a046aac1bc2621584857ed0c6fad9d63a0b24 Mon Sep 17 00:00:00 2001 From: kloepj2 Date: Wed, 26 Feb 2025 02:47:30 -0500 Subject: [PATCH 08/21] added explanation for baseType struct --- src/MeshField_Element.hpp | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/MeshField_Element.hpp b/src/MeshField_Element.hpp index c1c2bff..ed3c231 100644 --- a/src/MeshField_Element.hpp +++ b/src/MeshField_Element.hpp @@ -59,9 +59,14 @@ struct FieldElement { const ElementDofHolderAccessor elm2dofIn) : numMeshEnts(numMeshEntsIn), field(fieldIn), shapeFn(shapeFnIn), elm2dof(elm2dofIn) {} + /* general template for baseType which simply sets type + */ template struct baseType { using type = T; }; + /* template specialization to recursively strip type to get base type + * Example: int[5][6] => int[6] => int + */ template struct baseType { using type = typename baseType::type; }; From fd21bd1ecd28d54ad9834c58340458a9fbd89373 Mon Sep 17 00:00:00 2001 From: kloepj2 Date: Wed, 26 Feb 2025 11:14:17 -0500 Subject: [PATCH 09/21] Documented change to indexing --- src/MeshField.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/MeshField.hpp b/src/MeshField.hpp index 562522f..595457e 100644 --- a/src/MeshField.hpp +++ b/src/MeshField.hpp @@ -180,6 +180,7 @@ class OmegahMeshField { getMeshInfo(mesh_in), mesh_in.coords())) {} template + // Ordering of field indexing changed to 'entity, node, component' auto CreateLagrangeField() { return MeshField::CreateLagrangeField(meshInfo); From 526a287a97b4535e307fdda511af3a7586cf3f8d Mon Sep 17 00:00:00 2001 From: kloepj2 Date: Tue, 18 Mar 2025 23:11:39 -0400 Subject: [PATCH 10/21] Changed Cabana controller to create an aosoa for each datatype --- src/CabanaController.hpp | 41 +++++++++++++++++++++++++----------- src/MeshField_ShapeField.hpp | 6 +++--- test/testCabana.cpp | 22 +++++++++---------- 3 files changed, 43 insertions(+), 26 deletions(-) diff --git a/src/CabanaController.hpp b/src/CabanaController.hpp index 7fd3716..83008de 100644 --- a/src/CabanaController.hpp +++ b/src/CabanaController.hpp @@ -94,8 +94,6 @@ class CabanaController { private: // all the type defenitions that are needed us to get the type of the slice // returned by the underlying AoSoA - using soa_t = Cabana::SoA; - template using member_data_t = typename Cabana::MemberTypeAtIndex::type; @@ -112,6 +110,15 @@ class CabanaController { template using wrapper_slice_t = CabanaSliceWrapper, T>; + template + static auto construct(std::integer_sequence, + std::vector &obj) { + return std::tuple< + Cabana::AoSoA, MemorySpace, vecLen>...>{ + Cabana::AoSoA, MemorySpace, vecLen>( + "sliceAoSoA", obj[I])...}; + } + template void construct_sizes() { // There are no dynamic ranks w/ cabana controller. // So we only need to know extent. @@ -134,26 +141,35 @@ class CabanaController { } break; } - extent_sizes[theta][0] = num_tuples; + extent_sizes[theta][0] = num_tuples[this->theta]; this->theta += 1; if constexpr (sizeof...(Tx) != 0) construct_sizes(); } // member vaiables - Cabana::AoSoA aosoa; - const int num_tuples; + std::tuple, MemorySpace, vecLen>...> + aosoa; + // Cabana::AoSoA aosoa; + int num_tuples[sizeof...(Ts)]; unsigned short theta = 0; int extent_sizes[sizeof...(Ts)][MAX_RANK]; public: - CabanaController() : num_tuples(0) { + CabanaController() { static_assert(sizeof...(Ts) != 0); + for (int i = 0; i < sizeof...(Ts); ++i) { + num_tuples[i] = 0; + } construct_sizes(); } - - CabanaController(int n) : aosoa("sliceAoSoA", n), num_tuples(n) { + CabanaController(const std::initializer_list items) { static_assert(sizeof...(Ts) != 0); + std::vector obj(items); + for (std::size_t i = 0; i < obj.size(); ++i) { + num_tuples[i] = obj[i]; + } + aosoa = construct(std::make_integer_sequence{}, obj); construct_sizes(); } @@ -162,18 +178,19 @@ class CabanaController { assert(j >= 0); assert(j <= MAX_RANK); if (j == 0) - return this->tuples(); + return this->tuples(i); else return extent_sizes[i][j]; } - int tuples() const { return num_tuples; } + int tuples(int i) const { return num_tuples[i]; } template auto makeSlice() { // Creates wrapper object w/ meta data about the slice. using type = std::tuple_element_t; - const int stride = sizeof(soa_t) / sizeof(member_value_t); - auto slice = Cabana::slice(aosoa); + const int stride = sizeof(Cabana::SoA, vecLen>) / + sizeof(member_value_t); + auto slice = Cabana::slice<0>(std::get(aosoa)); int sizes[MAX_RANK]; for (int i = 0; i < MAX_RANK; i++) sizes[i] = this->size(index, i); diff --git a/src/MeshField_ShapeField.hpp b/src/MeshField_ShapeField.hpp index a1d4786..a88dc63 100644 --- a/src/MeshField_ShapeField.hpp +++ b/src/MeshField_ShapeField.hpp @@ -185,7 +185,7 @@ auto CreateLagrangeField(const MeshInfo &meshInfo) { Controller, MeshField::CabanaController>) { - return Ctrlr(numVtx); + return Ctrlr({numVtx}); } else { return Ctrlr({/*field 0*/ numVtx, 1, numComp}); } @@ -215,7 +215,7 @@ auto CreateLagrangeField(const MeshInfo &meshInfo) { Controller, MeshField::CabanaController>) { - return Ctrlr(std::max(numVtx, numEdge)); + return Ctrlr({numVtx, numEdge}); } else { return Ctrlr({/*field 0*/ numVtx, 1, numComp, /*field 1*/ numEdge, 1, numComp}); @@ -275,7 +275,7 @@ auto CreateCoordinateField(const MeshInfo &meshInfo) { Controller, MeshField::CabanaController>) { - return Ctrlr(numVtx); + return Ctrlr({numVtx}); } else { return Ctrlr({/*field 0*/ numVtx, 1, numComp}); } diff --git a/test/testCabana.cpp b/test/testCabana.cpp index 02988ba..7c4d035 100644 --- a/test/testCabana.cpp +++ b/test/testCabana.cpp @@ -36,7 +36,7 @@ using MemorySpace = Kokkos::DefaultExecutionSpace::memory_space; void testMakeSliceCabana(int num_tuples) { printf("== START testMakeSliceCabana ==\n"); using Ctrl = MeshField::CabanaController; - Ctrl c(num_tuples); + Ctrl c({num_tuples}); auto field0 = MeshField::makeField(c); @@ -54,7 +54,7 @@ void testParallelReduceCabana() { printf("== START testParallelReduceCabana ==\n"); using Ctrl = MeshField::CabanaController; const int N = 9; - Ctrl c(N); + Ctrl c({N}); { double result = 0, verify = 0; @@ -132,7 +132,7 @@ void testCabanaControllerSize() { using simple = MeshField::CabanaController; - simple c1(a); + simple c1({a}); for (int i = 0; i < 2; i++) { assert(c1.size(0, i) == psi[i]); } @@ -140,7 +140,7 @@ void testCabanaControllerSize() { using multi = MeshField::CabanaController; - multi c2(a); + multi c2({a, a, a}); for (int i = 0; i < 3; i++) { for (int j = 0; j < 4; j++) { assert(c2.size(i, j) == psi[j]); @@ -150,7 +150,7 @@ void testCabanaControllerSize() { using varied = MeshField::CabanaController; - varied c3(a); + varied c3({a, a, a, a}); for (int i = 0; i < 3; i++) assert(c3.size(0, i) == psi[i]); assert(c3.size(1, 0) == psi[0]); @@ -182,9 +182,9 @@ void testCabanaFieldSize() { MeshField::CabanaController; using empty = MeshField::CabanaController; - simple c1(a); - multi c2(a); - varied c3(a); + simple c1({a}); + multi c2({a, a, a}); + varied c3({a, a, a, a}); empty c4; const int MAX_RANK = 4; @@ -241,7 +241,7 @@ void testCabanaParallelFor() { using simd_ctrlr = MeshField::CabanaController; - simd_ctrlr c1(x); + simd_ctrlr c1({x, x, x, x}); auto field0 = MeshField::makeField(c1); auto field1 = MeshField::makeField(c1); auto field2 = MeshField::makeField(c1); @@ -285,7 +285,7 @@ void testParallelScan() { printf("== START testParallelScan ==\n"); using s_cab = MeshField::CabanaController; - s_cab c1(N); + s_cab c1({N, N}); auto pre = MeshField::makeField(c1); auto post = MeshField::makeField(c1); @@ -312,7 +312,7 @@ void testSetField() { const int N = 10; using cab1 = MeshField::CabanaController; - cab1 c1(N); + cab1 c1({N, N, N, N}); auto f1 = MeshField::makeField(c1); auto f2 = MeshField::makeField(c1); auto f3 = MeshField::makeField(c1); From 32a84465053345dcee0f77815cd7228a5cb895ce Mon Sep 17 00:00:00 2001 From: kloepj2 Date: Wed, 19 Mar 2025 12:22:11 -0400 Subject: [PATCH 11/21] adjusted size testing to also test different tuple sizes for each datatype --- test/testCabana.cpp | 37 +++++++++++++++++++++++++++++++++++++ 1 file changed, 37 insertions(+) diff --git a/test/testCabana.cpp b/test/testCabana.cpp index 7c4d035..dbd5227 100644 --- a/test/testCabana.cpp +++ b/test/testCabana.cpp @@ -231,6 +231,43 @@ void testCabanaFieldSize() { assert(field0.size(i) == 0); } } + multi diffc2({a, a + 1, a + 2}); + varied diffc3({a, a + 3, a + 4, a + 5}); + { + auto field0 = MeshField::makeField(diffc2); + auto field1 = MeshField::makeField(diffc2); + auto field2 = MeshField::makeField(diffc2); + assert(field0.size(0) == a); + assert(field1.size(0) == a + 1); + assert(field2.size(0) == a + 2); + for (int i = 1; i < MAX_RANK; ++i) { + assert(field0.size(i) == psi[i]); + assert(field1.size(i) == psi[i]); + assert(field2.size(i) == psi[i]); + } + } + { + auto field0 = MeshField::makeField(diffc3); + auto field1 = MeshField::makeField(diffc3); + auto field2 = MeshField::makeField(diffc3); + auto field3 = MeshField::makeField(diffc3); + assert(field0.size(0) == a); + assert(field1.size(0) == a + 3); + assert(field2.size(0) == a + 4); + assert(field3.size(0) == a + 5); + for (int i = 1; i < 3; i++) { + assert(field0.size(i) == psi[i]); + } + for (int i = 1; i < 1; i++) { + assert(field1.size(i) == psi[i]); + } + for (int i = 1; i < 4; i++) { + assert(field2.size(i) == psi[i]); + } + for (int i = 1; i < 2; i++) { + assert(field3.size(i) == psi[i]); + } + } printf("== END testCabanaFieldSize ==\n"); } From 91f973646961e84df18ce3214c3a4aae482f711f Mon Sep 17 00:00:00 2001 From: kloepj2 Date: Wed, 19 Mar 2025 12:48:54 -0400 Subject: [PATCH 12/21] adjusted other tests to also deal with different tuple sizes --- test/testCabana.cpp | 184 +++++++++++++++++++++++++++++++++++--------- 1 file changed, 147 insertions(+), 37 deletions(-) diff --git a/test/testCabana.cpp b/test/testCabana.cpp index dbd5227..815a1b4 100644 --- a/test/testCabana.cpp +++ b/test/testCabana.cpp @@ -164,6 +164,29 @@ void testCabanaControllerSize() { for (int i = 0; i < 4; i++) { assert(c4.size(0, i) == 0); } + + multi diffc2({a, b, c}); + assert(diffc2.size(0, 0) == psi[0]); + assert(diffc2.size(1, 0) == psi[1]); + assert(diffc2.size(2, 0) == psi[2]); + for (int i = 0; i < 3; i++) { + for (int j = 1; j < 4; j++) { + assert(c2.size(i, j) == psi[j]); + } + } + + varied diffc3({a, b, c, d}); + assert(diffc3.size(0, 0) == psi[0]); + assert(diffc3.size(1, 0) == psi[1]); + assert(diffc3.size(2, 0) == psi[2]); + assert(diffc3.size(3, 0) == psi[3]); + for (int i = 1; i < 3; i++) + assert(diffc3.size(0, i) == psi[i]); + for (int i = 1; i < 4; i++) + assert(diffc3.size(2, i) == psi[i]); + for (int i = 1; i < 2; i++) + assert(diffc3.size(3, i) == psi[i]); + printf("== END testCabanaControllerSize ==\n"); } @@ -313,6 +336,45 @@ void testCabanaParallelFor() { MeshField::simd_parallel_for(c1, {0, 0, 0, 0}, {x, y, z, a}, vectorKernel4, "simple_loop"); } + { + using simd_ctrlr = + MeshField::CabanaController; + simd_ctrlr c1({x, x + 1, x + 2, x + 3}); + auto field0 = MeshField::makeField(c1); + auto field1 = MeshField::makeField(c1); + auto field2 = MeshField::makeField(c1); + auto field3 = MeshField::makeField(c1); + + auto vectorKernel = KOKKOS_LAMBDA(const int &i) { + field0(i) = i; + assert(field0(i) == i); + }; + MeshField::simd_parallel_for(c1, {0}, {x}, vectorKernel, "simple_loop"); + + auto vectorKernel2 = KOKKOS_LAMBDA(const int &i, const int &j) { + field1(i, j) = i + j; + assert(field1(i, j) == i + j); + }; + MeshField::simd_parallel_for(c1, {0, 0}, {x + 1, y}, vectorKernel2, + "simple_loop"); + + auto vectorKernel3 = + KOKKOS_LAMBDA(const int &i, const int &j, const int &k) { + field2(i, j, k) = i + j + k; + assert(field2(i, j, k) == i + j + k); + }; + MeshField::simd_parallel_for(c1, {0, 0, 0}, {x + 2, y, z}, vectorKernel3, + "simple_loop"); + + auto vectorKernel4 = + KOKKOS_LAMBDA(const int &i, const int &j, const int &k, const int &l) { + field3(i, j, k, l) = i + j + k + l; + assert(field3(i, j, k, l) == i + j + k + l); + }; + MeshField::simd_parallel_for(c1, {0, 0, 0, 0}, {x + 3, y, z, a}, + vectorKernel4, "simple_loop"); + } printf("== END testCabanaParallelFor() ==\n"); } @@ -349,43 +411,91 @@ void testSetField() { const int N = 10; using cab1 = MeshField::CabanaController; - cab1 c1({N, N, N, N}); - auto f1 = MeshField::makeField(c1); - auto f2 = MeshField::makeField(c1); - auto f3 = MeshField::makeField(c1); - auto f4 = MeshField::makeField(c1); - - Kokkos::View v1("1", N); - Kokkos::View v2("2", N, N); - Kokkos::View v3("3", N, N, N); - Kokkos::View v4("4", N, N, N, N); - - Kokkos::Array start = MeshFieldUtil::to_kokkos_array<4>({0, 0, 0, 0}); - Kokkos::Array end = MeshFieldUtil::to_kokkos_array<4>({N, N, N, N}); - Kokkos::MDRangePolicy> p(start, end); - - Kokkos::parallel_for( - "", p, - KOKKOS_LAMBDA(const int &i, const int &j, const int &k, const int &l) { - v1(i) += i; - v2(i, j) += i + j; - v3(i, j, k) += i + j + k; - v4(i, j, k, l) += i + j + k + l; - }); - - f1.set(v1); - f2.set(v2); - f3.set(v3); - f4.set(v4); - - Kokkos::parallel_for( - "", p, - KOKKOS_LAMBDA(const int &i, const int &j, const int &k, const int &l) { - assert(f1(i) == v1(i)); - assert(f2(i, j) == v2(i, j)); - assert(f3(i, j, k) == v3(i, j, k)); - assert(f4(i, j, k, l) == v4(i, j, k, l)); - }); + { + cab1 c1({N, N, N, N}); + auto f1 = MeshField::makeField(c1); + auto f2 = MeshField::makeField(c1); + auto f3 = MeshField::makeField(c1); + auto f4 = MeshField::makeField(c1); + + Kokkos::View v1("1", N); + Kokkos::View v2("2", N, N); + Kokkos::View v3("3", N, N, N); + Kokkos::View v4("4", N, N, N, N); + + Kokkos::Array start = MeshFieldUtil::to_kokkos_array<4>({0, 0, 0, 0}); + Kokkos::Array end = MeshFieldUtil::to_kokkos_array<4>({N, N, N, N}); + Kokkos::MDRangePolicy> p(start, end); + + Kokkos::parallel_for( + "", p, + KOKKOS_LAMBDA(const int &i, const int &j, const int &k, const int &l) { + v1(i) += i; + v2(i, j) += i + j; + v3(i, j, k) += i + j + k; + v4(i, j, k, l) += i + j + k + l; + }); + + f1.set(v1); + f2.set(v2); + f3.set(v3); + f4.set(v4); + + Kokkos::parallel_for( + "", p, + KOKKOS_LAMBDA(const int &i, const int &j, const int &k, const int &l) { + assert(f1(i) == v1(i)); + assert(f2(i, j) == v2(i, j)); + assert(f3(i, j, k) == v3(i, j, k)); + assert(f4(i, j, k, l) == v4(i, j, k, l)); + }); + } + { + cab1 diffc1({N, N + 1, N + 2, N + 3}); + auto f1 = MeshField::makeField(diffc1); + auto f2 = MeshField::makeField(diffc1); + auto f3 = MeshField::makeField(diffc1); + auto f4 = MeshField::makeField(diffc1); + + Kokkos::View v1("1", N); + Kokkos::View v2("2", N + 1, N); + Kokkos::View v3("3", N + 2, N, N); + Kokkos::View v4("4", N + 3, N, N, N); + + Kokkos::Array start = MeshFieldUtil::to_kokkos_array<4>({0, 0, 0, 0}); + Kokkos::Array end = MeshFieldUtil::to_kokkos_array<4>({N + 3, N, N, N}); + Kokkos::MDRangePolicy> p(start, end); + Kokkos::parallel_for( + "", p, + KOKKOS_LAMBDA(const int &i, const int &j, const int &k, const int &l) { + if (i < N) + v1(i) += i; + if (i < N + 1) + v2(i, j) += i + j; + if (i < N + 2) + v3(i, j, k) += i + j + k; + if (i < N + 3) + v4(i, j, k, l) += i + j + k + l; + }); + + f1.set(v1); + f2.set(v2); + f3.set(v3); + f4.set(v4); + + Kokkos::parallel_for( + "", p, + KOKKOS_LAMBDA(const int &i, const int &j, const int &k, const int &l) { + if (i < N) + assert(f1(i) == v1(i)); + if (i < N + 1) + assert(f2(i, j) == v2(i, j)); + if (i < N + 2) + assert(f3(i, j, k) == v3(i, j, k)); + if (i < N + 3) + assert(f4(i, j, k, l) == v4(i, j, k, l)); + }); + } printf("== END testSetField ==\n"); } From 58742ae1158d3f058187a99a265ae590b3f385df Mon Sep 17 00:00:00 2001 From: kloepj2 Date: Wed, 19 Mar 2025 16:01:45 -0400 Subject: [PATCH 13/21] attempt CI fix --- src/MeshField_ShapeField.hpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/MeshField_ShapeField.hpp b/src/MeshField_ShapeField.hpp index a88dc63..faa6ba0 100644 --- a/src/MeshField_ShapeField.hpp +++ b/src/MeshField_ShapeField.hpp @@ -2,9 +2,7 @@ #define MESHFIELD_SHAPEFIELD_HPP #include "KokkosController.hpp" -#ifdef MESHFIELDS_ENABLE_CABANA #include "CabanaController.hpp" -#endif #include "MeshField_Field.hpp" #include "MeshField_Shape.hpp" #include //decltype From 1f1f6421f8b03791e0525520ffb3bbeafe888f36 Mon Sep 17 00:00:00 2001 From: kloepj2 Date: Wed, 19 Mar 2025 16:05:08 -0400 Subject: [PATCH 14/21] clang format --- src/MeshField_ShapeField.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/MeshField_ShapeField.hpp b/src/MeshField_ShapeField.hpp index faa6ba0..39c7f1d 100644 --- a/src/MeshField_ShapeField.hpp +++ b/src/MeshField_ShapeField.hpp @@ -1,8 +1,8 @@ #ifndef MESHFIELD_SHAPEFIELD_HPP #define MESHFIELD_SHAPEFIELD_HPP -#include "KokkosController.hpp" #include "CabanaController.hpp" +#include "KokkosController.hpp" #include "MeshField_Field.hpp" #include "MeshField_Shape.hpp" #include //decltype From fa6da3dfd0265ff6a6ba6ba6dcbe062147b631d2 Mon Sep 17 00:00:00 2001 From: kloepj2 Date: Wed, 19 Mar 2025 16:21:24 -0400 Subject: [PATCH 15/21] ci fix --- src/MeshField_ShapeField.hpp | 21 ++++++++++++++++++++- 1 file changed, 20 insertions(+), 1 deletion(-) diff --git a/src/MeshField_ShapeField.hpp b/src/MeshField_ShapeField.hpp index 39c7f1d..0e6541f 100644 --- a/src/MeshField_ShapeField.hpp +++ b/src/MeshField_ShapeField.hpp @@ -1,7 +1,9 @@ #ifndef MESHFIELD_SHAPEFIELD_HPP #define MESHFIELD_SHAPEFIELD_HPP +#ifdef MESHFIELDS_ENABLE_CABANA #include "CabanaController.hpp" +#endif #include "KokkosController.hpp" #include "MeshField_Field.hpp" #include "MeshField_Shape.hpp" @@ -171,6 +173,7 @@ auto CreateLagrangeField(const MeshInfo &meshInfo) { if (meshInfo.numVtx <= 0) { fail("mesh has no vertices\n"); } +#ifdef MESHFIELDS_ENABLE_CABANA using Ctrlr = std::conditional_t< std::is_same_v< Controller, @@ -189,6 +192,10 @@ auto CreateLagrangeField(const MeshInfo &meshInfo) { } }; Ctrlr kk_ctrl = createController(1, meshInfo.numVtx); +#else + using Ctrlr = Controller; + Ctrlr kk_ctrl({/*field 0*/ meshInfo.numVtx, 1, 1}); +#endif auto vtxField = MeshField::makeField(kk_ctrl); using LA = LinearAccessor; using LinearLagrangeShapeField = ShapeField; @@ -201,6 +208,7 @@ auto CreateLagrangeField(const MeshInfo &meshInfo) { if (meshInfo.numEdge <= 0) { fail("mesh has no edges\n"); } +#ifdef MESHFIELDS_ENABLE_CABANA using Ctrlr = std::conditional_t< std::is_same_v< Controller, @@ -220,6 +228,12 @@ auto CreateLagrangeField(const MeshInfo &meshInfo) { } }; Ctrlr kk_ctrl = createController(1, meshInfo.numVtx, meshInfo.numEdge); +#else + using Ctrlr = + Controller; + Ctrlr kk_ctrl({/*field 0*/ meshInfo.numVtx, 1, 1, + /*field 1*/ meshInfo.numEdge, 1, 1}); +#endif auto vtxField = MeshField::makeField(kk_ctrl); auto edgeField = MeshField::makeField(kk_ctrl); using QA = QuadraticAccessor; @@ -261,7 +275,8 @@ auto CreateCoordinateField(const MeshInfo &meshInfo) { using DataType = Real; using MemorySpace = typename ExecutionSpace::memory_space; const int numComp = meshInfo.dim; - // FIXME Oversized cabana datatypes when numComp = 1|2 +// FIXME Oversized cabana datatypes when numComp = 1|2 +#ifdef MESHFIELDS_ENABLE_CABANA using Ctrlr = std::conditional_t< std::is_same_v< Controller, @@ -279,6 +294,10 @@ auto CreateCoordinateField(const MeshInfo &meshInfo) { } }; Ctrlr kk_ctrl = createController(numComp, meshInfo.numVtx); +#else + using Ctrlr = Controller; + Ctrlr kk_ctrl({/*field 0*/ meshInfo.numVtx, 1, numComp}); +#endif auto vtxField = MeshField::makeField(kk_ctrl); using LA = LinearAccessor; using LinearLagrangeShapeField = ShapeField; From 3a4d17bbf3ac32c0ed2c8c1d84536ec51944af7a Mon Sep 17 00:00:00 2001 From: kloepj2 Date: Wed, 19 Mar 2025 16:28:26 -0400 Subject: [PATCH 16/21] fixing test for ci --- test/testOmegahElement.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/test/testOmegahElement.cpp b/test/testOmegahElement.cpp index d1309f7..8a26426 100644 --- a/test/testOmegahElement.cpp +++ b/test/testOmegahElement.cpp @@ -145,6 +145,7 @@ int main(int argc, char **argv) { Kokkos::initialize(argc, argv); auto lib = Omega_h::Library(&argc, &argv); MeshField::Debug = true; +#ifdef MESHFIELDS_ENABLE_CABANA { auto mesh = createMeshTri18(lib); MeshField::OmegahMeshField omf( @@ -227,6 +228,7 @@ int main(int argc, char **argv) { } } } +#endif { auto mesh = createMeshTri18(lib); MeshField::OmegahMeshField omf( From ef6b98d4388e79ac6e8ebc13b2dd8f76855c98fe Mon Sep 17 00:00:00 2001 From: kloepj2 Date: Mon, 24 Mar 2025 16:46:21 -0400 Subject: [PATCH 17/21] removed commented code --- src/CabanaController.hpp | 1 - 1 file changed, 1 deletion(-) diff --git a/src/CabanaController.hpp b/src/CabanaController.hpp index 83008de..2f0b84a 100644 --- a/src/CabanaController.hpp +++ b/src/CabanaController.hpp @@ -150,7 +150,6 @@ class CabanaController { // member vaiables std::tuple, MemorySpace, vecLen>...> aosoa; - // Cabana::AoSoA aosoa; int num_tuples[sizeof...(Ts)]; unsigned short theta = 0; int extent_sizes[sizeof...(Ts)][MAX_RANK]; From b291111ae31c5da6f5727fded6439414da6151b1 Mon Sep 17 00:00:00 2001 From: kloepj2 Date: Fri, 28 Mar 2025 03:29:13 -0400 Subject: [PATCH 18/21] adjusted tests to be a function for both controllers --- test/testOmegahElement.cpp | 276 +++++++++++++++++-------------------- 1 file changed, 124 insertions(+), 152 deletions(-) diff --git a/test/testOmegahElement.cpp b/test/testOmegahElement.cpp index 8a26426..bd1bc20 100644 --- a/test/testOmegahElement.cpp +++ b/test/testOmegahElement.cpp @@ -141,26 +141,71 @@ void doFail(std::string_view order, std::string_view function, MeshField::fail(msg); } -int main(int argc, char **argv) { - Kokkos::initialize(argc, argv); - auto lib = Omega_h::Library(&argc, &argv); - MeshField::Debug = true; #ifdef MESHFIELDS_ENABLE_CABANA - { - auto mesh = createMeshTri18(lib); - MeshField::OmegahMeshField omf( - mesh); +template +auto makeLagrangeField( + MeshField::OmegahMeshField + omf) { + return omf.CreateLagrangeField(); +} +#endif +template +auto makeLagrangeField( + MeshField::OmegahMeshField + omf) { + return omf.CreateLagrangeField(); +} +#ifdef MESHFIELDS_ENABLE_CABANA +template +auto doTrianglePoint( + MeshField::OmegahMeshField omf, + Kokkos::View coords, int NumPtsPerElem, + LinearFunction LinFunc, FieldType field) { + return omf.triangleLocalPointEval( + coords, NumPtsPerElem, LinFunc, field); +} +#endif +template +auto doTrianglePoint( + MeshField::OmegahMeshField omf, + Kokkos::View coords, int NumPtsPerElem, + LinearFunction LinFunc, FieldType field) { + return omf.triangleLocalPointEval( + coords, NumPtsPerElem, LinFunc, field); +} +#ifdef MESHFIELDS_ENABLE_CABANA +template +auto doTrianglePoint( + MeshField::OmegahMeshField omf, + Kokkos::View coords, int NumPtsPerElem, + QuadraticFunction QuadFunc, FieldType field) { + return omf.triangleLocalPointEval( + coords, NumPtsPerElem, QuadFunc, field); +} +#endif +template +auto doTrianglePoint( + MeshField::OmegahMeshField omf, + Kokkos::View coords, int NumPtsPerElem, + QuadraticFunction QuadFunc, FieldType field) { + return omf.triangleLocalPointEval( + coords, NumPtsPerElem, QuadFunc, field); +} - // setup field with values from the analytic function - static const size_t OnePtPerElem = 1; - static const size_t ThreePtsPerElem = 3; - auto centroids = createElmAreaCoords( - mesh.nfaces(), {1 / 3.0, 1 / 3.0, 1 / 3.0}); - auto interior = - createElmAreaCoords(mesh.nfaces(), {0.1, 0.4, 0.5}); - auto vertex = - createElmAreaCoords(mesh.nfaces(), {0.0, 0.0, 1.0}); - // clang-format off +template