diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 97056378a0b89cd5abb172027cdaffcff92e0b16..069870bdc85a79725b8a4d8a70c8eaa9d8adfb25 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -93,6 +93,26 @@ IF (NEKTAR_USE_CUDA)
         -DTEST_PATH="${CMAKE_SOURCE_DIR}")
 ENDIF()
 
+add_executable(test_helmholtz test_helmholtz.cpp)
+target_link_libraries(test_helmholtz PRIVATE Operators)
+target_link_libraries(test_helmholtz PRIVATE ${NEKTAR++_LIBRARIES})
+target_link_libraries(test_helmholtz PRIVATE Boost::unit_test_framework)
+target_include_directories(test_helmholtz PRIVATE "${CMAKE_SOURCE_DIR}")
+target_include_directories(test_helmholtz PRIVATE ${NEKTAR++_INCLUDE_DIRS})
+target_compile_definitions(test_helmholtz PRIVATE
+    -DTEST_PATH="${CMAKE_SOURCE_DIR}")
+
+IF (NEKTAR_USE_CUDA)
+    add_executable(test_helmholtzcuda test_physderivcuda.cpp)
+    target_link_libraries(test_helmholtzcuda PRIVATE Operators)
+    target_link_libraries(test_helmholtzcuda PRIVATE ${NEKTAR++_LIBRARIES})
+    target_link_libraries(test_helmholtzcuda PRIVATE Boost::unit_test_framework)
+    target_include_directories(test_helmholtzcuda PRIVATE "${CMAKE_SOURCE_DIR}")
+    target_include_directories(test_helmholtzcuda PRIVATE ${NEKTAR++_INCLUDE_DIRS})
+    target_compile_definitions(test_helmholtzcuda PRIVATE
+        -DTEST_PATH="${CMAKE_SOURCE_DIR}")
+ENDIF()
+
 add_test(
   NAME BwdTrans
   COMMAND test_bwdtrans
@@ -154,3 +174,17 @@ IF (NEKTAR_USE_CUDA)
       WORKING_DIRECTORY "${CMAKE_SOURCE_DIR}"
     )
 ENDIF()
+
+add_test(
+  NAME Helmholtz
+  COMMAND test_helmholtz
+  WORKING_DIRECTORY "${CMAKE_SOURCE_DIR}"
+)
+
+IF (NEKTAR_USE_CUDA)
+    add_test(
+      NAME HelmholtzCUDA
+      COMMAND test_helmholtzcuda
+      WORKING_DIRECTORY "${CMAKE_SOURCE_DIR}"
+    )
+ENDIF()
diff --git a/tests/init_fields.hpp b/tests/init_fields.hpp
index 1fd31225e5224385ede7bc3bee4d136257b8be7b..a40c622082e22b6c58f7be50fd77f0b2e6e551b0 100644
--- a/tests/init_fields.hpp
+++ b/tests/init_fields.hpp
@@ -128,11 +128,13 @@ public:
             {
                 for (size_t phys = 0; phys < block.num_pts; ++phys)
                 {
-                    if (fabs(*(outptr++) - *(expptr++)) > tol)
+                    if (fabs(*outptr - *expptr) > tol)
                     {
                         printf("%04zu %04zu %20.16f %20.16f %20.16f\n", el,
                                phys, *outptr, *expptr, fabs(*outptr - *expptr));
                     }
+                    expptr++;
+                    outptr++;
                 }
             }
         }
diff --git a/tests/init_helmholtzfields.hpp b/tests/init_helmholtzfields.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..ed2c007b31ec61f7eda6a69cfbc677611b93ad93
--- /dev/null
+++ b/tests/init_helmholtzfields.hpp
@@ -0,0 +1,188 @@
+#include "init_fields.hpp"
+
+using namespace std;
+using namespace Nektar::Operators;
+using namespace Nektar::LibUtilities;
+using namespace Nektar;
+
+class HelmholtzField
+    : public InitFields<double, FieldState::Coeff, FieldState::Coeff>
+{
+public:
+    HelmholtzField() : InitFields<double, FieldState::Coeff, FieldState::Coeff>()
+    {
+    }
+
+    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 coeff = 0; coeff < block.num_pts; ++coeff)
+                {
+                    *(inptr++) = coeff;
+                }
+            }
+            if (padding)
+            {
+                for (size_t el = 0; el < block.num_padding_elements; ++el)
+                {
+                    for (size_t coeff = 0; coeff < block.num_pts; ++coeff)
+                    {
+                        inptr++;
+                    }
+                }
+            }
+        }
+    }
+
+    void ExpectedSolution(const std::vector<BlockAttributes> &blocks,
+                          double *inptr)
+    {
+        Array<OneD, NekDouble> incoeffs(fixt_explist->GetNcoeffs());
+        Array<OneD, NekDouble> bwdtrans(fixt_explist->GetTotPoints());
+
+        Array<OneD, NekDouble> deriv(fixt_explist->GetCoordim(0) *
+                                       fixt_explist->GetTotPoints());
+        Array<OneD, NekDouble> deriv0 = deriv;
+        Array<OneD, NekDouble> deriv1 = deriv0 + fixt_explist->GetTotPoints();
+        Array<OneD, NekDouble> deriv2 = deriv1 + fixt_explist->GetTotPoints();
+        Array<OneD, Array<OneD, NekDouble>> derivarray(fixt_explist->GetCoordim(0));
+        if (fixt_explist->GetCoordim(0) > 0)
+        {
+            derivarray[0] = deriv0;
+        }
+        if (fixt_explist->GetCoordim(0) > 1)
+        {
+            derivarray[1] = deriv1;
+        }
+        if (fixt_explist->GetCoordim(0) > 2)
+        {
+            derivarray[2] = deriv2;
+        }
+        Array<OneD, NekDouble> mass(fixt_explist->GetNcoeffs());
+        Array<OneD, NekDouble> outcoeffs(fixt_explist->GetNcoeffs());
+
+        // Set test case
+        SetTestCase(blocks, incoeffs.get(), false);
+
+        // Calculate expected result from Nektar++
+        fixt_explist->BwdTrans(incoeffs, bwdtrans);
+        fixt_explist->PhysDeriv(bwdtrans, deriv0, deriv1, deriv2);
+        fixt_explist->IProductWRTDerivBase(derivarray, outcoeffs);
+        fixt_explist->IProductWRTBase(bwdtrans, mass);
+
+        // Copy expected result from Array to fixt_expected
+        double *coeffptr = outcoeffs.get();
+        double *massptr  = mass.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++) = (*coeffptr++) + (*massptr++);
+                }
+            }
+            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 HelmholtzField
+{
+public:
+    Seg()
+    {
+        meshName = "line.xml";
+    }
+};
+
+class Quad : public HelmholtzField
+{
+public:
+    Quad()
+    {
+        meshName = "square.xml";
+    }
+};
+
+class Tri : public HelmholtzField
+{
+public:
+    Tri()
+    {
+        meshName = "run/tri.xml";
+    }
+};
+
+class SquareAllElements : public HelmholtzField
+{
+public:
+    SquareAllElements()
+    {
+        meshName = "run/square_all_elements.xml";
+    }
+};
+
+class Hex : public HelmholtzField
+{
+public:
+    Hex()
+    {
+        meshName = "run/hex.xml";
+    }
+};
+
+class Prism : public HelmholtzField
+{
+public:
+    Prism()
+    {
+        meshName = "run/prism.xml";
+    }
+};
+
+class Pyr : public HelmholtzField
+{
+public:
+    Pyr()
+    {
+        meshName = "run/pyr.xml";
+    }
+};
+
+class Tet : public HelmholtzField
+{
+public:
+    Tet()
+    {
+        meshName = "run/tet.xml";
+    }
+};
+
+class CubePrismHex : public HelmholtzField
+{
+public:
+    CubePrismHex()
+    {
+        meshName = "run/cube_prismhex.xml";
+    }
+};
+
+class CubeAllElements : public HelmholtzField
+{
+public:
+    CubeAllElements()
+    {
+        meshName = "run/cube_all_elements.xml";
+    }
+};
diff --git a/tests/test_helmholtz.cpp b/tests/test_helmholtz.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..a9af5137b5dc76bb4d08055d8bbd4bd25c8fae3f
--- /dev/null
+++ b/tests/test_helmholtz.cpp
@@ -0,0 +1,168 @@
+#define BOOST_TEST_MODULE TestHelmholtz
+#include <boost/test/tools/output_test_stream.hpp>
+#include <boost/test/unit_test.hpp>
+
+#include <iostream>
+#include <memory>
+
+#include "Operators/OperatorHelmholtz.hpp"
+#include "init_helmholtzfields.hpp"
+
+BOOST_AUTO_TEST_SUITE(TestHelmholtz)
+
+using namespace std;
+using namespace Nektar::Operators;
+using namespace Nektar::LibUtilities;
+using namespace Nektar;
+
+BOOST_FIXTURE_TEST_CASE(helmholtz_seg, Seg)
+{
+    Configure();
+    SetTestCase(fixt_in->GetBlocks(), fixt_in->GetStorage().GetCPUPtr());
+    Helmholtz<>::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(helmholtz_quad, Quad)
+{
+    Configure();
+    SetTestCase(fixt_in->GetBlocks(), fixt_in->GetStorage().GetCPUPtr());
+    Helmholtz<>::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(helmholtz_tri, Tri)
+{
+    Configure();
+    SetTestCase(fixt_in->GetBlocks(), fixt_in->GetStorage().GetCPUPtr());
+    Helmholtz<>::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(helmholtz_square_all_elements, SquareAllElements)
+{
+    Configure();
+    SetTestCase(fixt_in->GetBlocks(), fixt_in->GetStorage().GetCPUPtr());
+    Helmholtz<>::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(helmholtz_hex, Hex)
+{
+    Configure();
+    SetTestCase(fixt_in->GetBlocks(), fixt_in->GetStorage().GetCPUPtr());
+    Helmholtz<>::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(helmholtz_prism, Prism)
+{
+    Configure();
+    SetTestCase(fixt_in->GetBlocks(), fixt_in->GetStorage().GetCPUPtr());
+    Helmholtz<>::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(helmholtz_pyr, Pyr)
+{
+    Configure();
+    SetTestCase(fixt_in->GetBlocks(), fixt_in->GetStorage().GetCPUPtr());
+    Helmholtz<>::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(helmholtz_tet, Tet)
+{
+    Configure();
+    SetTestCase(fixt_in->GetBlocks(), fixt_in->GetStorage().GetCPUPtr());
+    Helmholtz<>::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(helmholtz_cube_prism_hex, CubePrismHex)
+{
+    Configure();
+    SetTestCase(fixt_in->GetBlocks(), fixt_in->GetStorage().GetCPUPtr());
+    Helmholtz<>::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(helmholtz_cube_all_elements, CubeAllElements)
+{
+    Configure();
+    SetTestCase(fixt_in->GetBlocks(), fixt_in->GetStorage().GetCPUPtr());
+    Helmholtz<>::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()
diff --git a/tests/test_helmholtzcuda.cpp b/tests/test_helmholtzcuda.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..5a8de80086dde6b7315bd6531db5bf828e94278a
--- /dev/null
+++ b/tests/test_helmholtzcuda.cpp
@@ -0,0 +1,198 @@
+#define BOOST_TEST_MODULE TestHelmholtzCUDA
+#include <boost/test/tools/output_test_stream.hpp>
+#include <boost/test/unit_test.hpp>
+
+#include <iostream>
+#include <memory>
+
+#include "Operators/OperatorHelmholtz.hpp"
+#include "init_helmholtzfields.hpp"
+
+BOOST_AUTO_TEST_SUITE(TestHelmholtzCUDA)
+
+using namespace std;
+using namespace Nektar::Operators;
+using namespace Nektar::LibUtilities;
+using namespace Nektar;
+
+BOOST_FIXTURE_TEST_CASE(helmholtzcuda_seg, Seg)
+{
+    Configure();
+    SetTestCase(fixt_in->GetBlocks(), fixt_in->GetStorage().GetCPUPtr());
+    SetTestCase(
+        fixtcuda_in->GetBlocks(),
+        fixtcuda_in->template GetStorage<MemoryRegionCUDA>().GetCPUPtr());
+    Helmholtz<>::create(fixt_explist, "StdMat")->apply(*fixt_in, *fixt_expected);
+    Helmholtz<>::create(fixt_explist, "CUDA")
+        ->apply(*fixtcuda_in, *fixtcuda_out);
+    BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
+    boost::test_tools::output_test_stream output;
+    {
+        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
+                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+    }
+}
+
+BOOST_FIXTURE_TEST_CASE(helmholtzcuda_quad, Quad)
+{
+    Configure();
+    SetTestCase(fixt_in->GetBlocks(), fixt_in->GetStorage().GetCPUPtr());
+    SetTestCase(
+        fixtcuda_in->GetBlocks(),
+        fixtcuda_in->template GetStorage<MemoryRegionCUDA>().GetCPUPtr());
+    Helmholtz<>::create(fixt_explist, "StdMat")->apply(*fixt_in, *fixt_expected);
+    Helmholtz<>::create(fixt_explist, "CUDA")
+        ->apply(*fixtcuda_in, *fixtcuda_out);
+    BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
+    boost::test_tools::output_test_stream output;
+    {
+        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
+                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+    }
+}
+
+BOOST_FIXTURE_TEST_CASE(helmholtzcuda_tri, Tri)
+{
+    Configure();
+    SetTestCase(fixt_in->GetBlocks(), fixt_in->GetStorage().GetCPUPtr());
+    SetTestCase(
+        fixtcuda_in->GetBlocks(),
+        fixtcuda_in->template GetStorage<MemoryRegionCUDA>().GetCPUPtr());
+    Helmholtz<>::create(fixt_explist, "StdMat")->apply(*fixt_in, *fixt_expected);
+    Helmholtz<>::create(fixt_explist, "CUDA")
+        ->apply(*fixtcuda_in, *fixtcuda_out);
+    BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
+    boost::test_tools::output_test_stream output;
+    {
+        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
+                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+    }
+}
+
+BOOST_FIXTURE_TEST_CASE(helmholtzcuda_square_all_elements, SquareAllElements)
+{
+    Configure();
+    SetTestCase(fixt_in->GetBlocks(), fixt_in->GetStorage().GetCPUPtr());
+    SetTestCase(
+        fixtcuda_in->GetBlocks(),
+        fixtcuda_in->template GetStorage<MemoryRegionCUDA>().GetCPUPtr());
+    Helmholtz<>::create(fixt_explist, "StdMat")->apply(*fixt_in, *fixt_expected);
+    Helmholtz<>::create(fixt_explist, "CUDA")
+        ->apply(*fixtcuda_in, *fixtcuda_out);
+    BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
+    boost::test_tools::output_test_stream output;
+    {
+        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
+                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+    }
+}
+
+BOOST_FIXTURE_TEST_CASE(helmholtzcuda_hex, Hex)
+{
+    Configure();
+    SetTestCase(fixt_in->GetBlocks(), fixt_in->GetStorage().GetCPUPtr());
+    SetTestCase(
+        fixtcuda_in->GetBlocks(),
+        fixtcuda_in->template GetStorage<MemoryRegionCUDA>().GetCPUPtr());
+    Helmholtz<>::create(fixt_explist, "StdMat")->apply(*fixt_in, *fixt_expected);
+    Helmholtz<>::create(fixt_explist, "CUDA")
+        ->apply(*fixtcuda_in, *fixtcuda_out);
+    BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
+    boost::test_tools::output_test_stream output;
+    {
+        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
+                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+    }
+}
+
+BOOST_FIXTURE_TEST_CASE(helmholtzcuda_prism, Prism)
+{
+    Configure();
+    SetTestCase(fixt_in->GetBlocks(), fixt_in->GetStorage().GetCPUPtr());
+    SetTestCase(
+        fixtcuda_in->GetBlocks(),
+        fixtcuda_in->template GetStorage<MemoryRegionCUDA>().GetCPUPtr());
+    Helmholtz<>::create(fixt_explist, "StdMat")->apply(*fixt_in, *fixt_expected);
+    Helmholtz<>::create(fixt_explist, "CUDA")
+        ->apply(*fixtcuda_in, *fixtcuda_out);
+    BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
+    boost::test_tools::output_test_stream output;
+    {
+        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
+                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+    }
+}
+
+BOOST_FIXTURE_TEST_CASE(helmholtzcuda_pyr, Pyr)
+{
+    Configure();
+    SetTestCase(fixt_in->GetBlocks(), fixt_in->GetStorage().GetCPUPtr());
+    SetTestCase(
+        fixtcuda_in->GetBlocks(),
+        fixtcuda_in->template GetStorage<MemoryRegionCUDA>().GetCPUPtr());
+    Helmholtz<>::create(fixt_explist, "StdMat")->apply(*fixt_in, *fixt_expected);
+    Helmholtz<>::create(fixt_explist, "CUDA")
+        ->apply(*fixtcuda_in, *fixtcuda_out);
+    BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
+    boost::test_tools::output_test_stream output;
+    {
+        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
+                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+    }
+}
+
+BOOST_FIXTURE_TEST_CASE(helmholtzcuda_tet, Tet)
+{
+    Configure();
+    SetTestCase(fixt_in->GetBlocks(), fixt_in->GetStorage().GetCPUPtr());
+    SetTestCase(
+        fixtcuda_in->GetBlocks(),
+        fixtcuda_in->template GetStorage<MemoryRegionCUDA>().GetCPUPtr());
+    Helmholtz<>::create(fixt_explist, "StdMat")->apply(*fixt_in, *fixt_expected);
+    Helmholtz<>::create(fixt_explist, "CUDA")
+        ->apply(*fixtcuda_in, *fixtcuda_out);
+    BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
+    boost::test_tools::output_test_stream output;
+    {
+        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
+                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+    }
+}
+
+BOOST_FIXTURE_TEST_CASE(helmholtzcuda_cube_prism_hex, CubePrismHex)
+{
+    Configure();
+    SetTestCase(fixt_in->GetBlocks(), fixt_in->GetStorage().GetCPUPtr());
+    SetTestCase(
+        fixtcuda_in->GetBlocks(),
+        fixtcuda_in->template GetStorage<MemoryRegionCUDA>().GetCPUPtr());
+    Helmholtz<>::create(fixt_explist, "StdMat")->apply(*fixt_in, *fixt_expected);
+    Helmholtz<>::create(fixt_explist, "CUDA")
+        ->apply(*fixtcuda_in, *fixtcuda_out);
+    BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
+    boost::test_tools::output_test_stream output;
+    {
+        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
+                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+    }
+}
+
+BOOST_FIXTURE_TEST_CASE(helmholtzcuda_cube_all_elements, CubeAllElements)
+{
+    Configure();
+    SetTestCase(fixt_in->GetBlocks(), fixt_in->GetStorage().GetCPUPtr());
+    SetTestCase(
+        fixtcuda_in->GetBlocks(),
+        fixtcuda_in->template GetStorage<MemoryRegionCUDA>().GetCPUPtr());
+    Helmholtz<>::create(fixt_explist, "StdMat")->apply(*fixt_in, *fixt_expected);
+    Helmholtz<>::create(fixt_explist, "CUDA")
+        ->apply(*fixtcuda_in, *fixtcuda_out);
+    BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
+    boost::test_tools::output_test_stream output;
+    {
+        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
+                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+    }
+}
+
+BOOST_AUTO_TEST_SUITE_END()
diff --git a/tests/test_physderiv.cpp b/tests/test_physderiv.cpp
index ee582ca8c979d5cb5903f154116e2ab055678fd9..356f1f761bcef425a4c2ed4120b8f173ab3c4e2b 100644
--- a/tests/test_physderiv.cpp
+++ b/tests/test_physderiv.cpp
@@ -8,6 +8,8 @@
 #include "Operators/OperatorPhysDeriv.hpp"
 #include "init_physderivfields.hpp"
 
+BOOST_AUTO_TEST_SUITE(TestPhysDeriv)
+
 BOOST_FIXTURE_TEST_CASE(physderiv_seg, Seg)
 {
     Configure(1, 1);
@@ -157,3 +159,5 @@ BOOST_FIXTURE_TEST_CASE(physderiv_cube_all_elements, CubeAllElements)
                          fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
+
+BOOST_AUTO_TEST_SUITE_END()
diff --git a/tests/test_physderivcuda.cpp b/tests/test_physderivcuda.cpp
index 260c7557c16cbe282fd6e50cc787b42f6e4c5d7d..7bf5fd2ef2040eeb9c147f08350e5a8d144a1cb7 100644
--- a/tests/test_physderivcuda.cpp
+++ b/tests/test_physderivcuda.cpp
@@ -8,6 +8,8 @@
 #include "Operators/OperatorPhysDeriv.hpp"
 #include "init_physderivfields.hpp"
 
+BOOST_AUTO_TEST_SUITE(TestPhysDerivCUDA)
+
 BOOST_FIXTURE_TEST_CASE(physderivcuda_seg, Seg)
 {
     Configure(1, 1);
@@ -187,3 +189,5 @@ BOOST_FIXTURE_TEST_CASE(physderivcuda_cube_all_elements, CubeAllElements)
                          fixt_expected->GetStorage().GetCPUPtr(), 1.0E-10);
     }
 }
+
+BOOST_AUTO_TEST_SUITE_END()