diff --git a/CMakeLists.txt b/CMakeLists.txt index 43526d8534cefb67b8ae7dfcfb41f571845b98a1..a628e8608e034cf7bbde313da690e3b5c27b6bfc 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -47,6 +47,7 @@ set(SRC Field.cpp Operators/Operator.cpp Operators/OperatorFwdTrans.cpp Operators/FwdTrans/FwdTransImpl.cpp Operators/OperatorHelmSolve.cpp Operators/HelmSolve/HelmSolveImpl.cpp Operators/OperatorMatrix.cpp Operators/Matrix/MatrixImpl.cpp + Operators/OperatorAddTraceIntegral.cpp Operators/AddTraceIntegral/AddTraceIntegralImpl.cpp ) if (NEKTAR_USE_CUDA AND NEKTAR_USE_SIMD) diff --git a/Operators/AddTraceIntegral/AddTraceIntegralImpl.cpp b/Operators/AddTraceIntegral/AddTraceIntegralImpl.cpp new file mode 100644 index 0000000000000000000000000000000000000000..fd19d79fbd3588e8f0f77a450c54ad48053b61d1 --- /dev/null +++ b/Operators/AddTraceIntegral/AddTraceIntegralImpl.cpp @@ -0,0 +1,13 @@ +#include "AddTraceIntegralStdMat.hpp" + +namespace Nektar::Operators::detail +{ + +// Add different AddTraceIntegral implementations to the factory. +template <> +std::string OperatorAddTraceIntegralImpl<double, ImplStdMat>::className = + GetOperatorFactory<double>().RegisterCreatorFunction( + "AddTraceIntegralStdMat", + OperatorAddTraceIntegralImpl<double, ImplStdMat>::instantiate, "..."); + +} // namespace Nektar::Operators::detail diff --git a/Operators/AddTraceIntegral/AddTraceIntegralStdMat.hpp b/Operators/AddTraceIntegral/AddTraceIntegralStdMat.hpp new file mode 100644 index 0000000000000000000000000000000000000000..eddba30a09486ff926dadfb11c7a610763ab387c --- /dev/null +++ b/Operators/AddTraceIntegral/AddTraceIntegralStdMat.hpp @@ -0,0 +1,102 @@ +#pragma once + +#include "Operators/IProductWRTBase/IProductWRTBaseStdMat.hpp" +#include "Operators/OperatorAddTraceIntegral.hpp" + +using namespace Nektar::MultiRegions; + +namespace Nektar::Operators::detail +{ + +// standard matrix implementation +template <typename TData> +class OperatorAddTraceIntegralImpl<TData, ImplStdMat> + : public OperatorAddTraceIntegral<TData> +{ +public: + OperatorAddTraceIntegralImpl( + const MultiRegions::ExpListSharedPtr &expansionList) + : OperatorAddTraceIntegral<TData>(std::move(expansionList)), + m_wsp_trace(Field<TData, FieldState::Coeff>::create( + GetBlockAttributes(FieldState::Coeff, expansionList->GetTrace()))) + { + // Initialise Trace + const MultiRegions::ExpListSharedPtr &traceExpansionList = + expansionList->GetTrace(); + nTraceCoeffs = traceExpansionList->GetNcoeffs(); + nFieldCoeffs = expansionList->GetNcoeffs(); + + // Get Trace-to-Element Map + m_locTraceToTraceMap = expansionList->GetLocTraceToTraceMap(); + + // Initialise IProductWRTBase operator + m_IProductWRTBaseOp = IProductWRTBase<>::create( + this->m_expansionList->GetTrace(), "StdMat"); + } + + void apply(Field<TData, FieldState::Phys> &in, + Field<TData, FieldState::Coeff> &out) override + { + // Step 1: Inner product for trace integral + m_IProductWRTBaseOp->apply(in, m_wsp_trace); + + // Step 2: LEGACY Map Trace to Element + // TODO: Make this Field-only + Array<OneD, TData> inarray(nTraceCoeffs); + Array<OneD, TData> outarray(nFieldCoeffs, 0.0); + + // Copy data from input field + // Field -> Array + auto *inarrptr = inarray.data(); + auto *inptr = m_wsp_trace.GetStorage().GetCPUPtr(); + for (auto const &block : m_wsp_trace.GetBlocks()) + { + auto nSize = block.block_size; + auto nElmts = block.num_elements; + auto nmTot = block.num_pts; + + std::copy(inptr, inptr + nElmts * nmTot, inarrptr); + + inarrptr += nElmts * nmTot; + inptr += nSize; + } + + // The legacy map (that needs to be vectorised) + m_locTraceToTraceMap->AddTraceCoeffsToFieldCoeffs(inarray, outarray); + + // Copy data to output field + // Array -> Field + auto *outarrptr = outarray.data(); + auto *outptr = out.GetStorage().GetCPUPtr(); + for (auto const &block : out.GetBlocks()) + { + auto nSize = block.block_size; + auto nElmts = block.num_elements; + auto nmTot = block.num_pts; + + std::copy(outarrptr, outarrptr + nElmts * nmTot, outptr); + + outarrptr += nElmts * nmTot; + outptr += nSize; + } + } + + static std::unique_ptr<Operator<TData>> instantiate( + const MultiRegions::ExpListSharedPtr &expansionList) + { + return std::make_unique< + OperatorAddTraceIntegralImpl<TData, ImplStdMat>>(expansionList); + } + + static std::string className; + +private: + Field<TData, FieldState::Coeff> m_wsp_trace; + std::shared_ptr<OperatorIProductWRTBase<TData>> m_IProductWRTBaseOp; + Nektar::MultiRegions::LocTraceToTraceMapSharedPtr m_locTraceToTraceMap; + Array<OneD, Array<OneD, int>> m_traceCoeffsToElmtMap; + Array<OneD, Array<OneD, int>> m_traceCoeffsToElmtSign; + int nTraceCoeffs; + int nFieldCoeffs; +}; +} // namespace Nektar::Operators::detail diff --git a/Operators/OperatorAddTraceIntegral.cpp b/Operators/OperatorAddTraceIntegral.cpp new file mode 100644 index 0000000000000000000000000000000000000000..d32b2a9bf00834c8cc9ecdc763ddaf4cfeb0c50c --- /dev/null +++ b/Operators/OperatorAddTraceIntegral.cpp @@ -0,0 +1,9 @@ +#include "OperatorAddTraceIntegral.hpp" + +namespace Nektar::Operators +{ + +template <> const std::string AddTraceIntegral<>::key = "AddTraceIntegral"; +template <> const std::string AddTraceIntegral<>::default_impl = "StdMat"; + +} // namespace Nektar::Operators diff --git a/Operators/OperatorAddTraceIntegral.hpp b/Operators/OperatorAddTraceIntegral.hpp new file mode 100644 index 0000000000000000000000000000000000000000..37130a722869c14f35c5f55989ccc4738c027125 --- /dev/null +++ b/Operators/OperatorAddTraceIntegral.hpp @@ -0,0 +1,56 @@ +#pragma once + +#include "Field.hpp" +#include "Operator.hpp" + +namespace Nektar::Operators +{ + +// AddTraceIntegral base class +// Defines the apply operator to enforce apply parameter types +template <typename TData> +class OperatorAddTraceIntegral : public Operator<TData> +{ +public: + virtual ~OperatorAddTraceIntegral() = default; + + OperatorAddTraceIntegral( + const MultiRegions::ExpListSharedPtr &expansionList) + : Operator<TData>(expansionList) + { + } + + virtual void apply(Field<TData, FieldState::Phys> &in, + Field<TData, FieldState::Coeff> &out) = 0; + virtual void operator()(Field<TData, FieldState::Phys> &in, + Field<TData, FieldState::Coeff> &out) + { + apply(in, out); + } +}; + +// Descriptor / traits class for AddTraceIntegral +template <typename TData = default_fp_type> struct AddTraceIntegral +{ + using class_name = OperatorAddTraceIntegral<TData>; + using FieldIn = Field<TData, FieldState::Phys>; + using FieldOut = Field<TData, FieldState::Coeff>; + static const std::string key; + static const std::string default_impl; + + static std::shared_ptr<class_name> create( + const MultiRegions::ExpListSharedPtr &expansionList, + std::string pKey = "") + { + return Operator<TData>::template create<AddTraceIntegral>(expansionList, + pKey); + } +}; + +namespace detail +{ +// Template for AddTraceIntegral implementations +template <typename TData, typename Op> class OperatorAddTraceIntegralImpl; +} // namespace detail + +} // namespace Nektar::Operators diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index be911b26ad232d10c0abd982747ecf4467dd42f9..a37c18c314c20fe7740070e5321314bfd78cd464 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -340,6 +340,15 @@ IF (NEKTAR_USE_CUDA) -DTEST_PATH="${CMAKE_SOURCE_DIR}") ENDIF() +add_executable(test_addtraceintegral test_addtraceintegral.cpp) +target_link_libraries(test_addtraceintegral PRIVATE Operators) +target_link_libraries(test_addtraceintegral PRIVATE ${NEKTAR++_LIBRARIES}) +target_link_libraries(test_addtraceintegral PRIVATE Boost::unit_test_framework) +target_include_directories(test_addtraceintegral PRIVATE "${CMAKE_SOURCE_DIR}") +target_include_directories(test_addtraceintegral PRIVATE ${NEKTAR++_INCLUDE_DIRS}) +target_compile_definitions(test_addtraceintegral PRIVATE + -DTEST_PATH="${CMAKE_SOURCE_DIR}") + IF (NEKTAR_USE_CUDA) add_test( NAME CUDAKernels @@ -573,3 +582,9 @@ IF (NEKTAR_USE_CUDA) WORKING_DIRECTORY "${CMAKE_SOURCE_DIR}" ) ENDIF() + +add_test( + NAME AddTraceIntegral + COMMAND test_addtraceintegral + WORKING_DIRECTORY "${CMAKE_SOURCE_DIR}" +) diff --git a/tests/init_addtraceintegralfields.hpp b/tests/init_addtraceintegralfields.hpp new file mode 100644 index 0000000000000000000000000000000000000000..fc8943bbe2e67fbc92b55bf0516599e783f3bcb5 --- /dev/null +++ b/tests/init_addtraceintegralfields.hpp @@ -0,0 +1,198 @@ +#include "init_fields.hpp" + +using namespace std; +using namespace Nektar::Operators; +using namespace Nektar::LibUtilities; +using namespace Nektar; + +class AddTraceIntegralField + : public InitFields<double, FieldState::Phys, FieldState::Coeff, + MultiRegions::DisContField> +{ +public: + AddTraceIntegralField() + : InitFields<double, FieldState::Phys, FieldState::Coeff, + MultiRegions::DisContField>() + { + } + + /* + * Re-Initialise the input blocks based on the Trace-ExpList for this + * operator Delete previouisly defined fixt_in (also for CUDA) and re-define + * input based on TraceExpList + */ + void ReConfigure(size_t nin = 1, size_t nout = 1) + { + if (fixt_in) + { + delete fixt_in; + } + const FieldState stateIn = FieldState::Phys; + auto blocks_in = + GetBlockAttributes(stateIn, fixt_explist->GetTrace(), vec_t::width); + auto f_in = + Field<double, stateIn>::create(blocks_in, nin, vec_t::alignment); + fixt_in = new Field<double, stateIn>(std::move(f_in)); + +#ifdef NEKTAR_USE_CUDA + if (fixtcuda_in) + { + delete fixtcuda_in; + } + auto fcuda_in = + Field<double, stateIn>::template create<MemoryRegionCUDA>(blocks_in, + nin); + fixtcuda_in = new Field<double, stateIn>(std::move(fcuda_in)); +#endif + } + + void SetTestCase(const std::vector<BlockAttributes> &blocks, double *inptr, + bool padding = true) + { + for (auto const &block : fixt_in->GetBlocks()) + { + for (size_t el = 0; el < block.num_elements; ++el) + { + for (size_t phys = 0; phys < block.num_pts; ++phys) + { + *(inptr++) = phys; + } + } + if (padding) + { + for (size_t el = 0; el < block.num_padding_elements; ++el) + { + for (size_t phys = 0; phys < block.num_pts; ++phys) + { + inptr++; + } + } + } + } + } + + void ExpectedSolution(const std::vector<BlockAttributes> &blocks, + double *inptr) + { + auto fixt_explist_trace = fixt_explist->GetTrace(); + Array<OneD, NekDouble> inTracephys(fixt_explist_trace->GetNpoints(), + 0.0); + Array<OneD, NekDouble> outFieldcoeffs(fixt_explist->GetNcoeffs(), 0.0); + + // Set test case + SetTestCase(blocks, inTracephys.get(), false); + + // Calculate expected result from Nektar++ + fixt_explist->AddTraceIntegral(inTracephys, outFieldcoeffs); + + // Copy expected result from Array to fixt_expected + double *ptr = outFieldcoeffs.get(); + for (auto const &block : blocks) + { + for (size_t el = 0; el < block.num_elements; ++el) + { + for (size_t coeff = 0; coeff < block.num_pts; ++coeff) + { + (*inptr++) = (*ptr++); + } + } + for (size_t el = 0; el < block.num_padding_elements; ++el) + { + for (size_t coeff = 0; coeff < block.num_pts; ++coeff) + { + inptr++; + } + } + } + } +}; + +class Seg : public AddTraceIntegralField +{ +public: + Seg() + { + meshName = "run/line.xml"; + } +}; + +class Quad : public AddTraceIntegralField +{ +public: + Quad() + { + meshName = "run/square.xml"; + } +}; + +class Tri : public AddTraceIntegralField +{ +public: + Tri() + { + meshName = "run/tri.xml"; + } +}; + +class SquareAllElements : public AddTraceIntegralField +{ +public: + SquareAllElements() + { + meshName = "run/square_all_elements.xml"; + } +}; + +class Hex : public AddTraceIntegralField +{ +public: + Hex() + { + meshName = "run/hex.xml"; + } +}; + +class Prism : public AddTraceIntegralField +{ +public: + Prism() + { + meshName = "run/prism.xml"; + } +}; + +class Pyr : public AddTraceIntegralField +{ +public: + Pyr() + { + meshName = "run/pyr.xml"; + } +}; + +class Tet : public AddTraceIntegralField +{ +public: + Tet() + { + meshName = "run/tet.xml"; + } +}; + +class CubePrismHex : public AddTraceIntegralField +{ +public: + CubePrismHex() + { + meshName = "run/cube_prismhex.xml"; + } +}; + +class CubeAllElements : public AddTraceIntegralField +{ +public: + CubeAllElements() + { + meshName = "run/cube_all_elements.xml"; + } +}; diff --git a/tests/init_fields.hpp b/tests/init_fields.hpp index 381c33380a7acecca65e61061a376b82cbf8f729..97a752a98263b7dbecf0ec5f8d9516abb320dd6b 100644 --- a/tests/init_fields.hpp +++ b/tests/init_fields.hpp @@ -6,6 +6,7 @@ #include "Field.hpp" #include <MultiRegions/ContField.h> +#include <MultiRegions/DisContField.h> #include <MultiRegions/ExpList.h> #ifdef NEKTAR_USE_CUDA @@ -96,6 +97,14 @@ public: Collections::eNoCollection); } + else if constexpr (std::is_same_v<TExpList, MultiRegions::DisContField>) + { + fixt_explist = + MemoryManager<MultiRegions::DisContField>::AllocateSharedPtr( + session, graph, "u", true, true, + Collections::eNoCollection); + } + if constexpr (std::is_same_v<TExpList, MultiRegions::ExpList>) { fixt_explist = diff --git a/tests/test_addtraceintegral.cpp b/tests/test_addtraceintegral.cpp new file mode 100644 index 0000000000000000000000000000000000000000..8363487ddaea921733029dae6b36540d5a6e0f9d --- /dev/null +++ b/tests/test_addtraceintegral.cpp @@ -0,0 +1,193 @@ +#define BOOST_TEST_DYN_LINK +#define BOOST_TEST_MODULE TestAddTraceIntegral +#include <boost/test/tools/output_test_stream.hpp> +#include <boost/test/unit_test.hpp> + +#include <iostream> +#include <memory> + +#include "Operators/OperatorAddTraceIntegral.hpp" +#include "init_addtraceintegralfields.hpp" + +BOOST_AUTO_TEST_SUITE(TestAddTraceIntegral) + +using namespace std; +using namespace Nektar::Operators; +using namespace Nektar::LibUtilities; +using namespace Nektar; + +/* + * Currently fails in GetBlockAttributes + * GEometry is not initialised for Expansion(0) + * Possibly, because trace is not working for 1D expansions +BOOST_FIXTURE_TEST_CASE(addtraceintegral_seg, Seg) +{ + Configure(); + ReConfigure(); + SetTestCase(fixt_in->GetBlocks(), fixt_in->GetStorage().GetCPUPtr()); + AddTraceIntegral<>::create(fixt_explist, "StdMat")->apply(*fixt_in, +*fixt_out); ExpectedSolution(fixt_expected->GetBlocks(), + fixt_expected->GetStorage().GetCPUPtr()); + BOOST_TEST(fixt_out->compare(*fixt_expected, 1.0E-12)); + boost::test_tools::output_test_stream output; + { + OutputIfNotMatch(fixt_out->GetStorage().GetCPUPtr(), + fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12); + } +} +*/ + +BOOST_FIXTURE_TEST_CASE(addtraceintegral_quad, Quad) +{ + Configure(); + ReConfigure(); + SetTestCase(fixt_in->GetBlocks(), fixt_in->GetStorage().GetCPUPtr()); + AddTraceIntegral<>::create(fixt_explist, "StdMat") + ->apply(*fixt_in, *fixt_out); + ExpectedSolution(fixt_expected->GetBlocks(), + fixt_expected->GetStorage().GetCPUPtr()); + BOOST_TEST(fixt_out->compare(*fixt_expected, 1.0E-12)); + boost::test_tools::output_test_stream output; + { + OutputIfNotMatch(fixt_out->GetStorage().GetCPUPtr(), + fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12); + } +} + +BOOST_FIXTURE_TEST_CASE(addtraceintegral_tri, Tri) +{ + Configure(); + ReConfigure(); + SetTestCase(fixt_in->GetBlocks(), fixt_in->GetStorage().GetCPUPtr()); + AddTraceIntegral<>::create(fixt_explist, "StdMat") + ->apply(*fixt_in, *fixt_out); + ExpectedSolution(fixt_expected->GetBlocks(), + fixt_expected->GetStorage().GetCPUPtr()); + BOOST_TEST(fixt_out->compare(*fixt_expected, 1.0E-12)); + boost::test_tools::output_test_stream output; + { + OutputIfNotMatch(fixt_out->GetStorage().GetCPUPtr(), + fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12); + } +} + +BOOST_FIXTURE_TEST_CASE(addtraceintegral_square_all_elements, SquareAllElements) +{ + Configure(); + ReConfigure(); + SetTestCase(fixt_in->GetBlocks(), fixt_in->GetStorage().GetCPUPtr()); + AddTraceIntegral<>::create(fixt_explist, "StdMat") + ->apply(*fixt_in, *fixt_out); + ExpectedSolution(fixt_expected->GetBlocks(), + fixt_expected->GetStorage().GetCPUPtr()); + BOOST_TEST(fixt_out->compare(*fixt_expected, 1.0E-12)); + boost::test_tools::output_test_stream output; + { + OutputIfNotMatch(fixt_out->GetStorage().GetCPUPtr(), + fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12); + } +} + +BOOST_FIXTURE_TEST_CASE(addtraceintegral_hex, Hex) +{ + Configure(); + ReConfigure(); + SetTestCase(fixt_in->GetBlocks(), fixt_in->GetStorage().GetCPUPtr()); + AddTraceIntegral<>::create(fixt_explist, "StdMat") + ->apply(*fixt_in, *fixt_out); + ExpectedSolution(fixt_expected->GetBlocks(), + fixt_expected->GetStorage().GetCPUPtr()); + BOOST_TEST(fixt_out->compare(*fixt_expected, 1.0E-12)); + boost::test_tools::output_test_stream output; + { + OutputIfNotMatch(fixt_out->GetStorage().GetCPUPtr(), + fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12); + } +} + +BOOST_FIXTURE_TEST_CASE(addtraceintegral_prism, Prism) +{ + Configure(); + ReConfigure(); + SetTestCase(fixt_in->GetBlocks(), fixt_in->GetStorage().GetCPUPtr()); + AddTraceIntegral<>::create(fixt_explist, "StdMat") + ->apply(*fixt_in, *fixt_out); + ExpectedSolution(fixt_expected->GetBlocks(), + fixt_expected->GetStorage().GetCPUPtr()); + BOOST_TEST(fixt_out->compare(*fixt_expected, 1.0E-12)); + boost::test_tools::output_test_stream output; + { + OutputIfNotMatch(fixt_out->GetStorage().GetCPUPtr(), + fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12); + } +} + +BOOST_FIXTURE_TEST_CASE(addtraceintegral_pyr, Pyr) +{ + Configure(); + ReConfigure(); + SetTestCase(fixt_in->GetBlocks(), fixt_in->GetStorage().GetCPUPtr()); + AddTraceIntegral<>::create(fixt_explist, "StdMat") + ->apply(*fixt_in, *fixt_out); + ExpectedSolution(fixt_expected->GetBlocks(), + fixt_expected->GetStorage().GetCPUPtr()); + BOOST_TEST(fixt_out->compare(*fixt_expected, 1.0E-12)); + boost::test_tools::output_test_stream output; + { + OutputIfNotMatch(fixt_out->GetStorage().GetCPUPtr(), + fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12); + } +} + +BOOST_FIXTURE_TEST_CASE(addtraceintegral_tet, Tet) +{ + Configure(); + ReConfigure(); + SetTestCase(fixt_in->GetBlocks(), fixt_in->GetStorage().GetCPUPtr()); + AddTraceIntegral<>::create(fixt_explist, "StdMat") + ->apply(*fixt_in, *fixt_out); + ExpectedSolution(fixt_expected->GetBlocks(), + fixt_expected->GetStorage().GetCPUPtr()); + BOOST_TEST(fixt_out->compare(*fixt_expected, 1.0E-12)); + boost::test_tools::output_test_stream output; + { + OutputIfNotMatch(fixt_out->GetStorage().GetCPUPtr(), + fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12); + } +} + +BOOST_FIXTURE_TEST_CASE(addtraceintegral_cube_prism_hex, CubePrismHex) +{ + Configure(); + ReConfigure(); + SetTestCase(fixt_in->GetBlocks(), fixt_in->GetStorage().GetCPUPtr()); + AddTraceIntegral<>::create(fixt_explist, "StdMat") + ->apply(*fixt_in, *fixt_out); + ExpectedSolution(fixt_expected->GetBlocks(), + fixt_expected->GetStorage().GetCPUPtr()); + BOOST_TEST(fixt_out->compare(*fixt_expected, 1.0E-12)); + boost::test_tools::output_test_stream output; + { + OutputIfNotMatch(fixt_out->GetStorage().GetCPUPtr(), + fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12); + } +} + +BOOST_FIXTURE_TEST_CASE(addtraceintegral_cube_all_elements, CubeAllElements) +{ + Configure(); + ReConfigure(); + SetTestCase(fixt_in->GetBlocks(), fixt_in->GetStorage().GetCPUPtr()); + AddTraceIntegral<>::create(fixt_explist, "StdMat") + ->apply(*fixt_in, *fixt_out); + ExpectedSolution(fixt_expected->GetBlocks(), + fixt_expected->GetStorage().GetCPUPtr()); + BOOST_TEST(fixt_out->compare(*fixt_expected, 1.0E-12)); + boost::test_tools::output_test_stream output; + { + OutputIfNotMatch(fixt_out->GetStorage().GetCPUPtr(), + fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12); + } +} + +BOOST_AUTO_TEST_SUITE_END()