diff --git a/CHANGELOG b/CHANGELOG index 875e02e989..1b393bb1c9 100644 --- a/CHANGELOG +++ b/CHANGELOG @@ -16,24 +16,27 @@ GPU support - Added initial support for hardware devices, such as GPUs, and programming models, such as CUDA, OCCA, RAJA and OpenMP. -- The GPU/device support is based on MFEM's new "okina" kernel device interface - that supports various backends with dynamic execution, and MFEM's new - lightweight host/device memory manager "mm", designed to work seamlessly with - these device kernels. See the mem_manager.hpp and okina.hpp files in the - general/ directory. +- The GPU/device support is based on MFEM's new backends and kernels plus a new + lightweight device/host memory manager. The kernels can be implemented either + in OCCA, or as a simple wrapper around for loops, which can then be dispatched + to RAJA and native backends. See the forall.hpp and mem_manager.hpp files in + the general/ directory. -- Several of the MFEM example codes (ex1, ex1p, ex6, ex6p) can now take - advantage of GPU acceleration with backend selectable at runtime, see the +- Several of the MFEM example codes (ex1, ex1p, ex6, and ex6p) can now take + advantage of GPU acceleration with the backend selectable at runtime, see the Device::Configure method in general/device.hpp. - Many of the linear algebra and finite element operations (e.g. partially - assembled bilinear forms) have been extended to take advantage of the device - kernels by simply replacing loops with the MFEM_FORALL() macro. + assembled bilinear forms) have been extended to take advantage of kernel + acceleration by simply replacing loops with the MFEM_FORALL() macro. - In addition to pure CUDA, we support OCCA, RAJA and OpenMP kernels, which could be mixed and matched in different parts of the same application. Adding - support for additional programming models and devices should be easy without - modifying the user code. + support for additional programming models and devices in the future will not + require significant modifications to the user code. + +- The full list of currently supported backends is: "occa-cuda", "raja-cuda", + "cuda", "occa-omp", "raja-omp", "omp", "occa-cpu", "raja-cpu", and "cpu". Support for wedge elements and meshes with mixed element types -------------------------------------------------------------- @@ -138,9 +141,9 @@ Miscellaneous - Various other simplifications, extensions, and bugfixes in the code. - Renamed the option MFEM_USE_OPENMP to MFEM_USE_LEGACY_OPENMP. This legacy - option is now deprecated and planned for removal in a future release. The - original option name, MFEM_USE_OPENMP, is now used to enable the new OpenMP - backends in the okina device kernels. + option is deprecated and planned for removal in a future release. The original + option name, MFEM_USE_OPENMP, is now used to enable the new OpenMP backends in + the new kernels. API changes ----------- diff --git a/config/config.hpp b/config/config.hpp index 48cb6b1486..5f21e69f6f 100644 --- a/config/config.hpp +++ b/config/config.hpp @@ -15,6 +15,9 @@ // // Otherwise, use the local file: _config.hpp. +#ifndef MFEM_CONFIG_HPP +#define MFEM_CONFIG_HPP + #ifdef MFEM_BUILD_DIR #define MFEM_QUOTE(a) #a #define MFEM_MAKE_PATH(x,y) MFEM_QUOTE(x/y) @@ -52,3 +55,5 @@ #error Building with PUMI (MFEM_USE_PUMI=YES) requires MPI (MFEM_USE_MPI=YES) #endif #endif // MFEM_USE_MPI not defined + +#endif // MFEM_CONFIG_HPP diff --git a/examples/ex1.cpp b/examples/ex1.cpp index 44da30bb4c..1b5d6c2bfe 100644 --- a/examples/ex1.cpp +++ b/examples/ex1.cpp @@ -25,12 +25,13 @@ // ex1 -m ../data/mobius-strip.mesh // ex1 -m ../data/mobius-strip.mesh -o -1 -sc // -// Device runs: ex1 -p -d cuda -// ex1 -p -d 'cuda raja' -// ex1 -p -d 'cuda occa' -// ex1 -p -d 'omp raja' -// ex1 -p -d 'omp occa' -// ex1 -m ../data/beam-hex.mesh -p -d cuda +// Device sample runs: +// > ex1 -p -d cuda +// > ex1 -p -d raja-cuda +// > ex1 -p -d occa-cuda +// > ex1 -p -d raja-omp +// > ex1 -p -d occa-omp +// > ex1 -m ../data/beam-hex.mesh -p -d cuda // // Description: This example code demonstrates the use of MFEM to define a // simple finite element discretization of the Laplace problem @@ -61,7 +62,7 @@ int main(int argc, char *argv[]) int order = 1; bool static_cond = false; bool pa = false; - const char *device = ""; + const char *device = "cpu"; bool visualization = true; OptionsParser args(argc, argv); @@ -75,7 +76,7 @@ int main(int argc, char *argv[]) args.AddOption(&pa, "-p", "--pa", "-no-p", "--no-pa", "Enable Partial Assembly."); args.AddOption(&device, "-d", "--device", - "Device configuration, e.g. 'cuda', 'omp', 'raja', 'occa'."); + "Device configuration string, see Device::Configure()."); args.AddOption(&visualization, "-vis", "--visualization", "-no-vis", "--no-visualization", "Enable or disable GLVis visualization."); diff --git a/examples/ex1p.cpp b/examples/ex1p.cpp index a0f403cc3a..f3bd4ee362 100644 --- a/examples/ex1p.cpp +++ b/examples/ex1p.cpp @@ -25,9 +25,10 @@ // mpirun -np 4 ex1p -m ../data/mobius-strip.mesh // mpirun -np 4 ex1p -m ../data/mobius-strip.mesh -o -1 -sc // -// Device runs: mpirun -np 4 ex1p -p -d cuda -// mpirun -np 4 ex1p -p -d occa -// mpirun -np 4 ex1p -p -d 'raja omp' +// Device sample runs: +// > mpirun -np 4 ex1p -p -d cuda +// > mpirun -np 4 ex1p -p -d occa +// > mpirun -np 4 ex1p -p -d raja-omp // // Description: This example code demonstrates the use of MFEM to define a // simple finite element discretization of the Laplace problem @@ -64,7 +65,7 @@ int main(int argc, char *argv[]) int order = 1; bool static_cond = false; bool pa = false; - const char *device = ""; + const char *device = "cpu"; bool visualization = true; OptionsParser args(argc, argv); @@ -78,7 +79,7 @@ int main(int argc, char *argv[]) args.AddOption(&pa, "-p", "--pa", "-no-p", "--no-pa", "Enable Partial Assembly."); args.AddOption(&device, "-d", "--device", - "Device configuration, e.g. 'cuda', 'omp', 'raja', 'occa'."); + "Device configuration string, see Device::Configure()."); args.AddOption(&visualization, "-vis", "--visualization", "-no-vis", "--no-visualization", "Enable or disable GLVis visualization."); @@ -209,11 +210,6 @@ int main(int argc, char *argv[]) a->FormLinearSystem(ess_tdof_list, x, *b, A, X, B); - if (myid == 0) - { - cout << "Size of linear system: " << A->Height() << endl; - } - // 13. Solve the linear system A X = B. if (!pa) { diff --git a/examples/ex6.cpp b/examples/ex6.cpp index 596099a890..a76914364b 100644 --- a/examples/ex6.cpp +++ b/examples/ex6.cpp @@ -15,9 +15,10 @@ // ex6 -m ../data/square-disc-surf.mesh -o 2 // ex6 -m ../data/amr-quad.mesh // -// Device runs: ex6 -p -d cuda -// ex6 -p -d occa -// ex6 -p -d 'raja omp' +// Device sample runs: +// > ex6 -p -d cuda +// > ex6 -p -d occa +// > ex6 -p -d raja-omp // // Description: This is a version of Example 1 with a simple adaptive mesh // refinement loop. The problem being solved is again the Laplace @@ -48,7 +49,7 @@ int main(int argc, char *argv[]) const char *mesh_file = "../data/star.mesh"; int order = 1; bool pa = false; - const char *device = ""; + const char *device = "cpu"; bool visualization = true; OptionsParser args(argc, argv); @@ -59,7 +60,7 @@ int main(int argc, char *argv[]) args.AddOption(&pa, "-p", "--pa", "-no-p", "--no-pa", "Enable Partial Assembly."); args.AddOption(&device, "-d", "--device", - "Device configuration, e.g. 'cuda', 'omp', 'raja', 'occa'."); + "Device configuration string, see Device::Configure()."); args.AddOption(&visualization, "-vis", "--visualization", "-no-vis", "--no-visualization", "Enable or disable GLVis visualization."); diff --git a/examples/ex6p.cpp b/examples/ex6p.cpp index 0e36c86766..4dd272d2b6 100644 --- a/examples/ex6p.cpp +++ b/examples/ex6p.cpp @@ -15,9 +15,10 @@ // mpirun -np 4 ex6p -m ../data/square-disc-surf.mesh -o 2 // mpirun -np 4 ex6p -m ../data/amr-quad.mesh // -// Device runs: mpirun -np 4 ex6p -p -d cuda -// mpirun -np 4 ex6p -p -d occa -// mpirun -np 4 ex6p -p -d 'raja omp' +// Device sample runs: +// > mpirun -np 4 ex6p -p -d cuda +// > mpirun -np 4 ex6p -p -d occa +// > mpirun -np 4 ex6p -p -d raja-omp // // Description: This is a version of Example 1 with a simple adaptive mesh // refinement loop. The problem being solved is again the Laplace @@ -54,7 +55,7 @@ int main(int argc, char *argv[]) const char *mesh_file = "../data/star.mesh"; int order = 1; bool pa = false; - const char *device = ""; + const char *device = "cpu"; bool visualization = true; OptionsParser args(argc, argv); @@ -65,7 +66,7 @@ int main(int argc, char *argv[]) args.AddOption(&pa, "-p", "--pa", "-no-p", "--no-pa", "Enable Partial Assembly."); args.AddOption(&device, "-d", "--device", - "Device configuration, e.g. 'cuda', 'omp', 'raja', 'occa'."); + "Device configuration string, see Device::Configure()."); args.AddOption(&visualization, "-vis", "--visualization", "-no-vis", "--no-visualization", "Enable or disable GLVis visualization."); diff --git a/fem/bilinearform.cpp b/fem/bilinearform.cpp index e24acc4477..28e9fda304 100644 --- a/fem/bilinearform.cpp +++ b/fem/bilinearform.cpp @@ -85,7 +85,7 @@ BilinearForm::BilinearForm (FiniteElementSpace * f, switch (assembly) { case AssemblyLevel::FULL: - if (Device::UsingOkina()) + if (Device::IsEnabled()) { mfem_error("Full assembly not supported yet in device mode!"); } diff --git a/fem/bilinearform_ext.cpp b/fem/bilinearform_ext.cpp index 9a85de10d0..37c13fc700 100644 --- a/fem/bilinearform_ext.cpp +++ b/fem/bilinearform_ext.cpp @@ -12,10 +12,8 @@ // Implementations of classes FABilinearFormExtension, EABilinearFormExtension, // PABilinearFormExtension and MFBilinearFormExtension. -#include "fem.hpp" -#include "bilininteg.hpp" -#include "../general/okina.hpp" -#include "../linalg/dtensor.hpp" +#include "../general/forall.hpp" +#include "bilinearform.hpp" namespace mfem { diff --git a/fem/bilinearform_ext.hpp b/fem/bilinearform_ext.hpp index 415a998e3a..f6242d8fc7 100644 --- a/fem/bilinearform_ext.hpp +++ b/fem/bilinearform_ext.hpp @@ -13,13 +13,7 @@ #define MFEM_BILINEARFORM_EXT #include "../config/config.hpp" -#include "../linalg/linalg.hpp" #include "fespace.hpp" -#include "gridfunc.hpp" -#include "linearform.hpp" -#include "bilininteg.hpp" -#include "staticcond.hpp" -#include "hybridization.hpp" namespace mfem { diff --git a/fem/bilininteg.hpp b/fem/bilininteg.hpp index 43432fa0b3..1b66943872 100644 --- a/fem/bilininteg.hpp +++ b/fem/bilininteg.hpp @@ -1651,13 +1651,13 @@ private: int dim, ne, dofs1D, quad1D; public: /// Construct a diffusion integrator with coefficient Q = 1 - DiffusionIntegrator() { Q = NULL; MQ = NULL; maps=NULL; geom = NULL; } + DiffusionIntegrator() { Q = NULL; MQ = NULL; maps = NULL; geom = NULL; } /// Construct a diffusion integrator with a scalar coefficient q - DiffusionIntegrator (Coefficient &q) : Q(&q) { MQ = NULL; maps=NULL; geom = NULL; } + DiffusionIntegrator (Coefficient &q) : Q(&q) { MQ = NULL; maps = NULL; geom = NULL; } /// Construct a diffusion integrator with a matrix coefficient q - DiffusionIntegrator (MatrixCoefficient &q) : MQ(&q) { Q = NULL; maps=NULL; geom = NULL; } + DiffusionIntegrator (MatrixCoefficient &q) : MQ(&q) { Q = NULL; maps = NULL; geom = NULL; } /** Given a particular Finite Element computes the element stiffness matrix elmat. */ diff --git a/fem/bilininteg_ext.cpp b/fem/bilininteg_ext.cpp index 8a2ec769e0..c9a0a6913f 100644 --- a/fem/bilininteg_ext.cpp +++ b/fem/bilininteg_ext.cpp @@ -9,23 +9,25 @@ // terms of the GNU Lesser General Public License (as published by the Free // Software Foundation) version 2.1 dated February 1999. -#include "../config/config.hpp" -#include "../general/okina.hpp" -#include "../linalg/vector.hpp" +#include "../general/forall.hpp" +#include "bilininteg.hpp" +#include "gridfunc.hpp" -#include "fem.hpp" #include #include #include #include -#include "bilininteg.hpp" -#include "bilininteg_ext.hpp" using namespace std; namespace mfem { +#ifdef MFEM_USE_OCCA +typedef std::pair id_t; +typedef std::map occa_kernel_t; +#endif // MFEM_USE_OCCA + static const IntegrationRule &DefaultGetRule(const FiniteElement &trial_fe, const FiniteElement &test_fe) { @@ -58,16 +60,22 @@ static void OccaPADiffusionAssemble2D(const int D1D, const double COEFF, double *op) { - const occa::memory o_W = mfem::OccaPtr(W); - const occa::memory o_J = mfem::OccaPtr(J); - occa::memory o_op = mfem::OccaPtr(op); - occa::properties props; props["defines/D1D"] = D1D; props["defines/Q1D"] = Q1D; - - MFEM_NEW_OCCA_KERNEL(DiffusionSetup2D, "fem/occa.okl", props); - DiffusionSetup2D(NE, o_W, o_J, COEFF, o_op); + const occa::memory o_W = mfem::OccaPtr(W); + const occa::memory o_J = mfem::OccaPtr(J); + occa::memory o_op = mfem::OccaPtr(op); + const id_t id = std::make_pair(D1D,Q1D); + static occa_kernel_t OccaDiffSetup2D_ker; + if (OccaDiffSetup2D_ker.find(id) == OccaDiffSetup2D_ker.end()) + { + const occa::kernel DiffusionSetup2D = + mfem::OccaDev().buildKernel("occa://mfem/fem/occa.okl", + "DiffusionSetup2D", props); + OccaDiffSetup2D_ker.emplace(id, DiffusionSetup2D); + } + OccaDiffSetup2D_ker.at(id)(NE, o_W, o_J, COEFF, o_op); } static void OccaPADiffusionAssemble3D(const int D1D, @@ -78,16 +86,22 @@ static void OccaPADiffusionAssemble3D(const int D1D, const double COEFF, double *op) { - const occa::memory o_W = mfem::OccaPtr(W); - const occa::memory o_J = mfem::OccaPtr(J); - occa::memory o_op = mfem::OccaPtr(op); - occa::properties props; props["defines/D1D"] = D1D; props["defines/Q1D"] = Q1D; - - MFEM_NEW_OCCA_KERNEL(DiffusionSetup3D, "fem/occa.okl", props); - DiffusionSetup3D(NE, o_W, o_J, COEFF, o_op); + const occa::memory o_W = mfem::OccaPtr(W); + const occa::memory o_J = mfem::OccaPtr(J); + occa::memory o_op = mfem::OccaPtr(op); + const id_t id = std::make_pair(D1D,Q1D); + static occa_kernel_t OccaDiffSetup3D_ker; + if (OccaDiffSetup3D_ker.find(id) == OccaDiffSetup3D_ker.end()) + { + const occa::kernel DiffusionSetup3D = + mfem::OccaDev().buildKernel("occa://mfem/fem/occa.okl", + "DiffusionSetup3D", props); + OccaDiffSetup3D_ker.emplace(id, DiffusionSetup3D); + } + OccaDiffSetup3D_ker.at(id)(NE, o_W, o_J, COEFF, o_op); } #endif // MFEM_USE_OCCA @@ -170,6 +184,24 @@ static void PADiffusionAssemble3D(const int Q1D, }); } +namespace internal +{ + +#ifdef MFEM_USE_OCCA +// This function is currently used to determine if an OCCA kernel should be +// used. +static bool DeviceUseOcca() +{ + return Device::Allows(Backend::OCCA_CUDA) || + (Device::Allows(Backend::OCCA_OMP) && + !Device::Allows(Backend::DEVICE_MASK)) || + (Device::Allows(Backend::OCCA_CPU) && + !Device::Allows(Backend::DEVICE_MASK|Backend::OMP_MASK)); +} +#endif + +} + static void PADiffusionAssemble(const int dim, const int D1D, const int Q1D, @@ -183,7 +215,7 @@ static void PADiffusionAssemble(const int dim, if (dim == 2) { #ifdef MFEM_USE_OCCA - if (Device::UsingOcca()) + if (internal::DeviceUseOcca()) { OccaPADiffusionAssemble2D(D1D, Q1D, NE, W, J, COEFF, op); return; @@ -194,7 +226,7 @@ static void PADiffusionAssemble(const int dim, if (dim == 3) { #ifdef MFEM_USE_OCCA - if (Device::UsingOcca()) + if (internal::DeviceUseOcca()) { OccaPADiffusionAssemble3D(D1D, Q1D, NE, W, J, COEFF, op); return; @@ -237,6 +269,9 @@ static void OccaPADiffusionMultAdd2D(const int D1D, const double* x, double* y) { + occa::properties props; + props["defines/D1D"] = D1D; + props["defines/Q1D"] = Q1D; const occa::memory o_B = mfem::OccaPtr(B); const occa::memory o_G = mfem::OccaPtr(G); const occa::memory o_Bt = mfem::OccaPtr(Bt); @@ -244,20 +279,30 @@ static void OccaPADiffusionMultAdd2D(const int D1D, const occa::memory o_op = mfem::OccaPtr(op); const occa::memory o_x = mfem::OccaPtr(x); occa::memory o_y = mfem::OccaPtr(y); - - occa::properties props; - props["defines/D1D"] = D1D; - props["defines/Q1D"] = Q1D; - - if (!Device::UsingDevice()) + const id_t id = std::make_pair(D1D,Q1D); + if (!Device::Allows(Backend::OCCA_CUDA)) { - MFEM_NEW_OCCA_KERNEL(DiffusionApply2D_CPU, "fem/occa.okl", props); - DiffusionApply2D_CPU(NE, o_B, o_G, o_Bt, o_Gt, o_op, o_x, o_y); + static occa_kernel_t OccaDiffApply2D_cpu; + if (OccaDiffApply2D_cpu.find(id) == OccaDiffApply2D_cpu.end()) + { + const occa::kernel DiffusionApply2D_CPU = + mfem::OccaDev().buildKernel("occa://mfem/fem/occa.okl", + "DiffusionApply2D_CPU", props); + OccaDiffApply2D_cpu.emplace(id, DiffusionApply2D_CPU); + } + OccaDiffApply2D_cpu.at(id)(NE, o_B, o_G, o_Bt, o_Gt, o_op, o_x, o_y); } else { - MFEM_NEW_OCCA_KERNEL(DiffusionApply2D_GPU, "fem/occa.okl", props); - DiffusionApply2D_GPU(NE, o_B, o_G, o_Bt, o_Gt, o_op, o_x, o_y); + static occa_kernel_t OccaDiffApply2D_gpu; + if (OccaDiffApply2D_gpu.find(id) == OccaDiffApply2D_gpu.end()) + { + const occa::kernel DiffusionApply2D_GPU = + mfem::OccaDev().buildKernel("occa://mfem/fem/occa.okl", + "DiffusionApply2D_GPU", props); + OccaDiffApply2D_gpu.emplace(id, DiffusionApply2D_GPU); + } + OccaDiffApply2D_gpu.at(id)(NE, o_B, o_G, o_Bt, o_Gt, o_op, o_x, o_y); } } @@ -273,6 +318,9 @@ static void OccaPADiffusionMultAdd3D(const int D1D, const double* x, double* y) { + occa::properties props; + props["defines/D1D"] = D1D; + props["defines/Q1D"] = Q1D; const occa::memory o_B = mfem::OccaPtr(B); const occa::memory o_G = mfem::OccaPtr(G); const occa::memory o_Bt = mfem::OccaPtr(Bt); @@ -280,20 +328,30 @@ static void OccaPADiffusionMultAdd3D(const int D1D, const occa::memory o_op = mfem::OccaPtr(op); const occa::memory o_x = mfem::OccaPtr(x); occa::memory o_y = mfem::OccaPtr(y); - - occa::properties props; - props["defines/D1D"] = D1D; - props["defines/Q1D"] = Q1D; - - if (!Device::UsingDevice()) + const id_t id = std::make_pair(D1D,Q1D); + if (!Device::Allows(Backend::OCCA_CUDA)) { - MFEM_NEW_OCCA_KERNEL(DiffusionApply3D_CPU, "fem/occa.okl", props); - DiffusionApply3D_CPU(NE, o_B, o_G, o_Bt, o_Gt, o_op, o_x, o_y); + static occa_kernel_t OccaDiffApply3D_cpu; + if (OccaDiffApply3D_cpu.find(id) == OccaDiffApply3D_cpu.end()) + { + const occa::kernel DiffusionApply3D_CPU = + mfem::OccaDev().buildKernel("occa://mfem/fem/occa.okl", + "DiffusionApply3D_CPU", props); + OccaDiffApply3D_cpu.emplace(id, DiffusionApply3D_CPU); + } + OccaDiffApply3D_cpu.at(id)(NE, o_B, o_G, o_Bt, o_Gt, o_op, o_x, o_y); } else { - MFEM_NEW_OCCA_KERNEL(DiffusionApply3D_GPU, "fem/occa.okl", props); - DiffusionApply3D_GPU(NE, o_B, o_G, o_Bt, o_Gt, o_op, o_x, o_y); + static occa_kernel_t OccaDiffApply3D_gpu; + if (OccaDiffApply3D_gpu.find(id) == OccaDiffApply3D_gpu.end()) + { + const occa::kernel DiffusionApply3D_GPU = + mfem::OccaDev().buildKernel("occa://mfem/fem/occa.okl", + "DiffusionApply3D_GPU", props); + OccaDiffApply3D_gpu.emplace(id, DiffusionApply3D_GPU); + } + OccaDiffApply3D_gpu.at(id)(NE, o_B, o_G, o_Bt, o_Gt, o_op, o_x, o_y); } } #endif // MFEM_USE_OCCA @@ -616,7 +674,7 @@ static void PADiffusionMultAssembled(const int dim, double* y) { #ifdef MFEM_USE_OCCA - if (Device::UsingOcca()) + if (internal::DeviceUseOcca()) { if (dim == 2) { @@ -783,25 +841,38 @@ static void OccaPAMassMultAdd2D(const int D1D, const double* x, double* y) { + occa::properties props; + props["defines/D1D"] = D1D; + props["defines/Q1D"] = Q1D; const occa::memory o_B = mfem::OccaPtr(B); const occa::memory o_Bt = mfem::OccaPtr(Bt); const occa::memory o_op = mfem::OccaPtr(op); const occa::memory o_x = mfem::OccaPtr(x); occa::memory o_y = mfem::OccaPtr(y); - - occa::properties props; - props["defines/D1D"] = D1D; - props["defines/Q1D"] = Q1D; - - if (!Device::UsingDevice()) + const id_t id = std::make_pair(D1D,Q1D); + if (!Device::Allows(Backend::OCCA_CUDA)) { - MFEM_NEW_OCCA_KERNEL(MassApply2D_CPU, "fem/occa.okl", props); - MassApply2D_CPU(NE, o_B, o_Bt, o_op, o_x, o_y); + static occa_kernel_t OccaMassApply2D_cpu; + if (OccaMassApply2D_cpu.find(id) == OccaMassApply2D_cpu.end()) + { + const occa::kernel MassApply2D_CPU = + mfem::OccaDev().buildKernel("occa://mfem/fem/occa.okl", + "MassApply2D_CPU", props); + OccaMassApply2D_cpu.emplace(id, MassApply2D_CPU); + } + OccaMassApply2D_cpu.at(id)(NE, o_B, o_Bt, o_op, o_x, o_y); } else { - MFEM_NEW_OCCA_KERNEL(MassApply2D_GPU, "fem/occa.okl", props); - MassApply2D_GPU(NE, o_B, o_Bt, o_op, o_x, o_y); + static occa_kernel_t OccaMassApply2D_gpu; + if (OccaMassApply2D_gpu.find(id) == OccaMassApply2D_gpu.end()) + { + const occa::kernel MassApply2D_GPU = + mfem::OccaDev().buildKernel("occa://mfem/fem/occa.okl", + "MassApply2D_GPU", props); + OccaMassApply2D_gpu.emplace(id, MassApply2D_GPU); + } + OccaMassApply2D_gpu.at(id)(NE, o_B, o_Bt, o_op, o_x, o_y); } } @@ -815,25 +886,38 @@ static void OccaPAMassMultAdd3D(const int D1D, const double* x, double* y) { + occa::properties props; + props["defines/D1D"] = D1D; + props["defines/Q1D"] = Q1D; const occa::memory o_B = mfem::OccaPtr(B); const occa::memory o_Bt = mfem::OccaPtr(Bt); const occa::memory o_op = mfem::OccaPtr(op); const occa::memory o_x = mfem::OccaPtr(x); occa::memory o_y = mfem::OccaPtr(y); - - occa::properties props; - props["defines/D1D"] = D1D; - props["defines/Q1D"] = Q1D; - - if (!Device::UsingDevice()) + const id_t id = std::make_pair(D1D,Q1D); + if (!Device::Allows(Backend::OCCA_CUDA)) { - MFEM_NEW_OCCA_KERNEL(MassApply3D_CPU, "fem/occa.okl", props); - MassApply3D_CPU(NE, o_B, o_Bt, o_op, o_x, o_y); + static occa_kernel_t OccaMassApply3D_cpu; + if (OccaMassApply3D_cpu.find(id) == OccaMassApply3D_cpu.end()) + { + const occa::kernel MassApply3D_CPU = + mfem::OccaDev().buildKernel("occa://mfem/fem/occa.okl", + "MassApply3D_CPU", props); + OccaMassApply3D_cpu.emplace(id, MassApply3D_CPU); + } + OccaMassApply3D_cpu.at(id)(NE, o_B, o_Bt, o_op, o_x, o_y); } else { - MFEM_NEW_OCCA_KERNEL(MassApply3D_GPU, "fem/occa.okl", props); - MassApply3D_GPU(NE, o_B, o_Bt, o_op, o_x, o_y); + static occa_kernel_t OccaMassApply3D_gpu; + if (OccaMassApply3D_gpu.find(id) == OccaMassApply3D_gpu.end()) + { + const occa::kernel MassApply3D_GPU = + mfem::OccaDev().buildKernel("occa://mfem/fem/occa.okl", + "MassApply3D_GPU", props); + OccaMassApply3D_gpu.emplace(id, MassApply3D_GPU); + } + OccaMassApply3D_gpu.at(id)(NE, o_B, o_Bt, o_op, o_x, o_y); } } #endif // MFEM_USE_OCCA @@ -1071,7 +1155,7 @@ static void PAMassMultAssembled(const int dim, double* y) { #ifdef MFEM_USE_OCCA - if (Device::UsingOcca()) + if (internal::DeviceUseOcca()) { if (dim == 2) { diff --git a/fem/bilininteg_ext.hpp b/fem/bilininteg_ext.hpp index 4ad420a366..6f6b114311 100644 --- a/fem/bilininteg_ext.hpp +++ b/fem/bilininteg_ext.hpp @@ -12,6 +12,8 @@ #ifndef MFEM_BILININTEG_EXT #define MFEM_BILININTEG_EXT +#include "fespace.hpp" + namespace mfem { diff --git a/general/CMakeLists.txt b/general/CMakeLists.txt index 5ff7b4e7bc..ffabe783e1 100644 --- a/general/CMakeLists.txt +++ b/general/CMakeLists.txt @@ -42,7 +42,7 @@ list(APPEND HDRS mem_alloc.hpp mem_manager.hpp occa.hpp - okina.hpp + forall.hpp optparser.hpp osockstream.hpp sets.hpp diff --git a/general/array.cpp b/general/array.cpp index 9cd75dc572..f1280edae5 100644 --- a/general/array.cpp +++ b/general/array.cpp @@ -13,7 +13,7 @@ // Abstract array data type #include "array.hpp" -#include "../general/okina.hpp" +#include "../general/forall.hpp" #include namespace mfem diff --git a/general/cuda.cpp b/general/cuda.cpp index 0c2b576498..f56ff33ae2 100644 --- a/general/cuda.cpp +++ b/general/cuda.cpp @@ -9,7 +9,7 @@ // terms of the GNU Lesser General Public License (as published by the Free // Software Foundation) version 2.1 dated February 1999. -#include "okina.hpp" +#include "cuda.hpp" namespace mfem { diff --git a/general/cuda.hpp b/general/cuda.hpp index ce4d0f8e07..aaa104fef4 100644 --- a/general/cuda.hpp +++ b/general/cuda.hpp @@ -12,30 +12,67 @@ #ifndef MFEM_CUDA_HPP #define MFEM_CUDA_HPP -#include +#include "../config/config.hpp" +#include "error.hpp" #ifdef MFEM_USE_CUDA +#include #include #endif -namespace mfem -{ +// CUDA block size used by MFEM. +#define MFEM_CUDA_BLOCKS 256 #ifdef MFEM_USE_CUDA -#define MFEM_DEVICE __device__ -#define MFEM_HOST_DEVICE __host__ __device__ -inline void CuCheck(const unsigned int c) -{ - MFEM_ASSERT(c == cudaSuccess, cudaGetErrorString(cudaGetLastError())); -} +#define MFEM_ATTR_DEVICE __device__ +#define MFEM_ATTR_HOST_DEVICE __host__ __device__ +// Define the CUDA debug macros: +// - MFEM_CUDA_CHECK_DRV(x) where 'x' returns/is type 'CUresult' +// - MFEM_CUDA_CHECK_RT(x) where 'x' returns/is type 'cudaError_t' +#ifdef MFEM_DEBUG +#define MFEM_CUDA_CHECK_DRV(x) \ + do \ + { \ + CUresult err = (x); \ + if (err != CUDA_SUCCESS) \ + { \ + const char *error_string; \ + cuGetErrorString(err, &error_string); \ + _MFEM_MESSAGE("CUDA error: (" << #x \ + << ") failed with error:\n --> " \ + << error_string, 0); \ + } \ + } \ + while (0) +#define MFEM_CUDA_CHECK_RT(x) \ + do \ + { \ + cudaError_t err = (x); \ + if (err != cudaSuccess) \ + { \ + _MFEM_MESSAGE("CUDA error: (" << #x \ + << ") failed with error:\n --> " \ + << cudaGetErrorString(err), 0); \ + } \ + } \ + while (0) +#else +#define MFEM_CUDA_CHECK_DRV(x) x +#define MFEM_CUDA_CHECK_RT(x) x +#endif #else // MFEM_USE_CUDA -#define MFEM_DEVICE -#define MFEM_HOST_DEVICE +#define MFEM_ATTR_DEVICE +#define MFEM_ATTR_HOST_DEVICE typedef int CUdevice; typedef int CUcontext; typedef void* CUstream; #endif // MFEM_USE_CUDA + +namespace mfem +{ + +// Define 'atomicAdd' function. #ifdef __CUDA_ARCH__ #if __CUDA_ARCH__ < 600 static __device__ inline double atomicAdd(double* address, double val) @@ -56,7 +93,7 @@ static __device__ inline double atomicAdd(double* address, double val) return __longlong_as_double(old); } #endif // __CUDA_ARCH__ < 600 -template MFEM_DEVICE +template MFEM_ATTR_DEVICE inline T AtomicAdd(T volatile *address, T val) { return atomicAdd((T *)address, val); @@ -72,6 +109,7 @@ template inline T AtomicAdd(T volatile *address, T val) } #endif // __CUDA_ARCH__ + /// Allocates device memory void* CuMemAlloc(void **d_ptr, size_t bytes); diff --git a/general/device.cpp b/general/device.cpp index 31cedc0d0b..3f9a00a232 100644 --- a/general/device.cpp +++ b/general/device.cpp @@ -9,29 +9,105 @@ // terms of the GNU Lesser General Public License (as published by the Free // Software Foundation) version 2.1 dated February 1999. -#include "okina.hpp" +#include "forall.hpp" #include "cuda.hpp" #include "occa.hpp" +#include +#include + namespace mfem { -CUstream *cuStream; +// Place the following variables in the mfem::internal namespace, so that they +// will not be included in the doxygen documentation. +namespace internal +{ + +CUstream *cuStream = NULL; static CUdevice cuDevice; static CUcontext cuContext; OccaDevice occaDevice; +// Backends listed by priority, high to low: +static const Backend::Id backend_list[Backend::NUM_BACKENDS] = +{ + Backend::OCCA_CUDA, Backend::RAJA_CUDA, Backend::CUDA, + Backend::OCCA_OMP, Backend::RAJA_OMP, Backend::OMP, + Backend::OCCA_CPU, Backend::RAJA_CPU, Backend::CPU +}; + +// Backend names listed by priority, high to low: +static const char *backend_name[Backend::NUM_BACKENDS] = +{ + "occa-cuda", "raja-cuda", "cuda", "occa-omp", "raja-omp", "omp", + "occa-cpu", "raja-cpu", "cpu" +}; + +} // namespace mfem::internal + +void Device::Configure(const std::string &device, const int dev) +{ + std::map bmap; + for (int i = 0; i < Backend::NUM_BACKENDS; i++) + { + bmap[internal::backend_name[i]] = internal::backend_list[i]; + } + std::string::size_type beg = 0, end; + while (1) + { + end = device.find(',', beg); + end = (end != std::string::npos) ? end : device.size(); + const std::string bname = device.substr(beg, end - beg); + std::map::iterator it = bmap.find(bname); + MFEM_VERIFY(it != bmap.end(), "invalid backend name: '" << bname << '\''); + Get().MarkBackend(it->second); + if (end == device.size()) { break; } + beg = end + 1; + } + + // OCCA_CUDA needs CUDA or RAJA_CUDA: + Get().allowed_backends = Get().backends; + if (Allows(Backend::OCCA_CUDA) && !Allows(Backend::RAJA_CUDA)) + { + Get().MarkBackend(Backend::CUDA); + } + + // Activate all backends for Setup(). + Get().allowed_backends = Get().backends; + Get().Setup(dev); + + // Enable only the default host CPU backend. + Get().allowed_backends = Backend::CPU; +} + +void Device::Print(std::ostream &out) +{ + out << "Device configuration: "; + bool add_comma = false; + for (int i = 0; i < Backend::NUM_BACKENDS; i++) + { + if (Get().backends & internal::backend_list[i]) + { + if (add_comma) { out << ','; } + add_comma = true; + out << internal::backend_name[i]; + } + } + out << '\n'; +} + #ifdef MFEM_USE_CUDA static void DeviceSetup(const int dev, int &ngpu) { cudaGetDeviceCount(&ngpu); MFEM_VERIFY(ngpu>0, "No CUDA device found!"); cuInit(0); - cuDeviceGet(&cuDevice,dev); - cuCtxCreate(&cuContext, CU_CTX_SCHED_AUTO, cuDevice); - cuStream = new CUstream; - MFEM_VERIFY(cuStream, "CUDA stream could not be created!"); - cuStreamCreate(cuStream, CU_STREAM_DEFAULT); + cuDeviceGet(&internal::cuDevice, dev); + cuCtxCreate(&internal::cuContext, CU_CTX_SCHED_AUTO, internal::cuDevice); + internal::cuStream = new CUstream; + MFEM_VERIFY(internal::cuStream, "CUDA stream could not be created!"); + cuStreamCreate(internal::cuStream, CU_STREAM_DEFAULT); } #endif @@ -39,36 +115,45 @@ static void CudaDeviceSetup(const int dev, int &ngpu) { #ifdef MFEM_USE_CUDA DeviceSetup(dev, ngpu); -#else - MFEM_ABORT("CUDA requested but MFEM was not built with MFEM_USE_CUDA=YES"); #endif } static void RajaDeviceSetup(const int dev, int &ngpu) { -#if defined(MFEM_USE_CUDA) && defined(MFEM_USE_RAJA) - DeviceSetup(dev, ngpu); -#elif !defined(MFEM_USE_RAJA) - MFEM_ABORT("RAJA requested but MFEM was not built with MFEM_USE_RAJA=YES"); +#ifdef MFEM_USE_CUDA + if (ngpu <= 0) { DeviceSetup(dev, ngpu); } #endif } static void OccaDeviceSetup(CUdevice cu_dev, CUcontext cu_ctx) { #ifdef MFEM_USE_OCCA - const bool omp = Device::UsingOmp(); - const bool cuda = Device::UsingCuda(); + const int cpu = Device::Allows(Backend::OCCA_CPU); + const int omp = Device::Allows(Backend::OCCA_OMP); + const int cuda = Device::Allows(Backend::OCCA_CUDA); + if (cpu + omp + cuda > 1) + { + MFEM_ABORT("Only one OCCA backend can be configured at a time!"); + } if (cuda) { - occaDevice = OccaWrapDevice(cu_dev, cu_ctx); +#if OCCA_CUDA_ENABLED + internal::occaDevice = occa::cuda::wrapDevice(cu_dev, cu_ctx); +#else + MFEM_ABORT("the OCCA CUDA backend requires OCCA built with CUDA!"); +#endif } else if (omp) { - occaDevice.setup("mode: 'OpenMP'"); +#if OCCA_OPENMP_ENABLED + internal::occaDevice.setup("mode: 'OpenMP'"); +#else + MFEM_ABORT("the OCCA OpenMP backend requires OCCA built with OpenMP!"); +#endif } else { - occaDevice.setup("mode: 'Serial'"); + internal::occaDevice.setup("mode: 'Serial'"); } std::string mfemDir; @@ -88,32 +173,46 @@ static void OccaDeviceSetup(CUdevice cu_dev, CUcontext cu_ctx) occa::io::addLibraryPath("mfem", mfemDir); occa::loadKernels("mfem"); #else - MFEM_ABORT("OCCA requested but MFEM was not built with MFEM_USE_OCCA=YES"); + MFEM_ABORT("the OCCA backends require MFEM built with MFEM_USE_OCCA=YES"); #endif } void Device::Setup(const int device) { + MFEM_VERIFY(ngpu == -1, "the mfem::Device is already configured!"); + + ngpu = 0; dev = device; - MFEM_ASSERT(ngpu==-1, "Only one MFEMDeviceSetup allowed"); - ngpu = 0; +#ifndef MFEM_USE_CUDA + MFEM_VERIFY(!Allows(Backend::CUDA_MASK), + "the CUDA backends require MFEM built with MFEM_USE_CUDA=YES"); +#endif +#ifndef MFEM_USE_RAJA + MFEM_VERIFY(!Allows(Backend::RAJA_MASK), + "the RAJA backends require MFEM built with MFEM_USE_RAJA=YES"); +#endif +#ifndef MFEM_USE_OPENMP + MFEM_VERIFY(!Allows(Backend::OMP|Backend::RAJA_OMP), + "the OpenMP and RAJA OpenMP backends require MFEM built with" + " MFEM_USE_OPENMP=YES"); +#endif + // The check for MFEM_USE_OCCA is in the function OccaDeviceSetup(). - // We initialize CUDA first so OccaDeviceSetup() can reuse the same - // initialized cuDevice and cuContext objects - if (cuda) { CudaDeviceSetup(dev, ngpu); } - if (raja) { RajaDeviceSetup(dev, ngpu); } - if (occa) { OccaDeviceSetup(cuDevice, cuContext); } - - if (cuda && ngpu==0) + // We initialize CUDA and/or RAJA_CUDA first so OccaDeviceSetup() can reuse + // the same initialized cuDevice and cuContext objects when OCCA_CUDA is + // enabled. + if (Allows(Backend::CUDA)) { CudaDeviceSetup(dev, ngpu); } + if (Allows(Backend::RAJA_CUDA)) { RajaDeviceSetup(dev, ngpu); } + if (Allows(Backend::OCCA_MASK)) { - MFEM_ABORT("CUDA requested but MFEM was not built with MFEM_USE_CUDA=YES"); + OccaDeviceSetup(internal::cuDevice, internal::cuContext); } } Device::~Device() { - if (cuda) { delete cuStream; } + delete internal::cuStream; } } // mfem diff --git a/general/device.hpp b/general/device.hpp index 47191b9ded..5bbd8462e6 100644 --- a/general/device.hpp +++ b/general/device.hpp @@ -17,23 +17,95 @@ namespace mfem { -/// The MFEM Device class that abstracts hardware devices, such as GPUs, and -/// programming models, such as CUDA, OCCA, RAJA and OpenMP. +/// MFEM backends. +/** Individual backends will generally implement only a subset of the kernels + implemented by the default CPU backend. The goal of the backends is to + accelerate data-parallel portions of the code and they can use a device + memory space (e.g. GPUs) or share the memory space of the host (OpenMP). */ +struct Backend +{ + /** @brief In the documentation below, we use square brackets to indicate the + type of the backend: host or device. */ + enum Id + { + /// [host] Default CPU backend: sequential execution on each MPI rank. + CPU = 1 << 0, + /// [host] OpenMP backend. Enabled when MFEM_USE_OPENMP = YES. + OMP = 1 << 1, + /// [device] CUDA backend. Enabled when MFEM_USE_CUDA = YES. + CUDA = 1 << 2, + /** @brief [host] RAJA CPU backend: sequential execution on each MPI rank. + Enabled when MFEM_USE_RAJA = YES. */ + RAJA_CPU = 1 << 3, + /** @brief [host] RAJA OpenMP backend. Enabled when MFEM_USE_RAJA = YES + and MFEM_USE_OPENMP = YES. */ + RAJA_OMP = 1 << 4, + /** @brief [device] RAJA CUDA backend. Enabled when MFEM_USE_RAJA = YES + and MFEM_USE_CUDA = YES. */ + RAJA_CUDA = 1 << 5, + /** @brief [host] OCCA CPU backend: sequential execution on each MPI rank. + Enabled when MFEM_USE_OCCA = YES. */ + OCCA_CPU = 1 << 6, + /// [host] OCCA OpenMP backend. Enabled when MFEM_USE_OCCA = YES. + OCCA_OMP = 1 << 7, + /** @brief [device] OCCA CUDA backend. Enabled when MFEM_USE_OCCA = YES + and MFEM_USE_CUDA = YES. */ + OCCA_CUDA = 1 << 8 + }; + + /** @brief Additional useful constants. For example, the *_MASK constants can + be used with Device::Allows(). */ + enum + { + /// Number of backends: from (1 << 0) to (1 << (NUM_BACKENDS-1)). + NUM_BACKENDS = 9, + /// Biwise-OR of all CUDA backends + CUDA_MASK = CUDA | RAJA_CUDA | OCCA_CUDA, + /// Biwise-OR of all RAJA backends + RAJA_MASK = RAJA_CPU | RAJA_OMP | RAJA_CUDA, + /// Biwise-OR of all OCCA backends + OCCA_MASK = OCCA_CPU | OCCA_OMP | OCCA_CUDA, + /// Biwise-OR of all OpenMP backends + OMP_MASK = OMP | RAJA_OMP | OCCA_OMP, + /// Biwise-OR of all device backends + DEVICE_MASK = CUDA_MASK + }; +}; + + +/** @brief The MFEM Device class abstracts hardware devices, such as GPUs, as + well as programming models, such as CUDA, OCCA, RAJA and OpenMP. */ +/** This class represents a "virtual device" with the following properties: + - There a single object of this class which is controlled by its static + methods. + - Once configured, the object cannot be re-configured during the program + lifetime. + - MFEM classes use this object to determine where (host or device) to + perform an operation and which backend implementation to use. + - Multiple backends can be configured at the same time; currently, a fixed + priority order is used to select a specific backend from the list of + configured backends. See the Backend class and the Configure() method in + this class for details. + - The device can be disabled to restrict the backend selection to only the + default host CPU backend, see the methods Enable() and Disable(). */ class Device { private: - enum MODES {HOST, DEVICE}; + enum MODES {SEQUENTIAL, ACCELERATED}; MODES mode; - int dev = 0; - int ngpu = -1; - bool cuda = false; - bool raja = false; - bool occa = false; - bool omp = false; - bool isTracking = true; + int dev = 0; ///< Device ID of the configured device. + int ngpu = -1; ///< Number of detected devices; -1: not initialized. + unsigned long backends; ///< Bitwise-OR of all configured backends. + /** Bitwise-OR mask of all allowed backends. All backends are active when the + Device is enabled. When the Device is disabled, only the host CPU backend + is allowed. */ + unsigned long allowed_backends; - Device(): mode{Device::HOST} {} + Device() + : mode(Device::SEQUENTIAL), + backends(Backend::CPU), + allowed_backends(backends) { } Device(Device const&); void operator=(Device const&); static Device& Get() { static Device singleton; return singleton; } @@ -41,90 +113,74 @@ private: /// Setup switcher based on configuration settings void Setup(const int dev = 0); + void MarkBackend(Backend::Id b) { backends |= b; } + public: + /// Configure the Device backends. + /** The string parameter @a device must be a comma-separated list of backend + string names (see below). The @a dev argument specifies the ID of the + actual devices (e.g. GPU) to use. + * The available backends are described by the Backend class. + * The string name of a backend is the lowercase version of the + Backend::Id enumeration constant with '_' replaced by '-', e.g. the + string name of 'RAJA_CPU' is 'raja-cpu'. + * The 'cpu' backend is always enabled with lowest priority. + * The current backend priority from highest to lowest is: 'occa-cuda', + 'raja-cuda', 'cuda', 'occa-omp', 'raja-omp', 'omp', 'occa-cpu', + 'raja-cpu', 'cpu'. + * Multiple backends can be configured at the same time. + * Only one 'occa-*' backend can be configured at a time. + * The backend 'occa-cuda' enables the 'cuda' backend unless 'raja-cuda' + is already enabled. + * After this call, the Device will be disabled. */ + static void Configure(const std::string &device, const int dev = 0); - /// Configure and enable the device. - /** The string parameter will enable a backend (cuda, omp, raja, occa) if the - corresponding substring is present (for now, the order is ignored). The - dev argument specifies which of the devices (e.g. GPUs) to enable. */ - static inline void Configure(std::string device, const int dev = 0) - { - if (device.find("cuda") != std::string::npos) { Device::UseCuda(); } - if (device.find("omp") != std::string::npos) { Device::UseOmp(); } - if (device.find("raja") != std::string::npos) { Device::UseRaja(); } - if (device.find("occa") != std::string::npos) { Device::UseOcca(); } - EnableDevice(dev); - } + /// Print the configuration of the MFEM virtual device object. + static void Print(std::ostream &out = mfem::out); - /// Print the configured device + programming models in order of priority - static inline void Print(std::ostream &out = mfem::out) - { - const bool omp = Device::UsingOmp(); - const bool cuda = Device::UsingCuda(); - const bool occa = Device::UsingOcca(); - const bool raja = Device::UsingRaja(); - out << "Device configuration: "; - if (cuda && occa) { out << "OCCA:CUDA\n"; return; } - if (omp && occa) { out << "OCCA:OpenMP\n"; return; } - if (occa) { out << "OCCA:CPU\n"; return; } - if (cuda && raja) { out << "RAJA:CUDA\n"; return; } - if (cuda) { out << "CUDA\n"; return; } - if (omp && raja) { out << "RAJA:OpenMP\n"; return; } - if (raja) { out << "RAJA:CPU\n"; return; } - if (omp) { out << "OpenMP\n"; return; } - out << "CPU\n"; - } + /// Return true if Configure() has been called previously. + static inline bool IsConfigured() { return Get().ngpu >= 0; } + + /// Return true if an actual device (e.g. GPU) has been configured. + static inline bool IsAvailable() { return Get().ngpu > 0; } /// Enable the use of the configured device in the code that follows. - /** In particular, use the device version of the okina kernels encountered, - with the device versions of the data registered in the memory manager - (copying host-to-device if necessary). */ - static inline void Enable() { Get().mode = Device::DEVICE; } + /** After this call MFEM classes will use the backend kernels whenever + possible, transferring data automatically to the device, if necessary. + + If the only configured backend is the default host CPU one, the device + will remain disabled. */ + static inline void Enable() + { + if (Get().backends & ~Backend::CPU) + { + Get().mode = Device::ACCELERATED; + Get().allowed_backends = Get().backends; + } + } /// Disable the use of the configured device in the code that follows. - /** In particular, use the host version of the okina kernels encountered, - with the host versions of the data registered in the memory manager - (copying device-to-host if necessary). */ - static inline void Disable() { Get().mode = Device::HOST; } - - constexpr static inline bool UsingMM() + /** After this call MFEM classes will only use default CPU kernels, + transferring data automatically from the device, if necessary. */ + static inline void Disable() { -#ifdef MFEM_USE_MM - return true; -#else - return false; -#endif + Get().mode = Device::SEQUENTIAL; + Get().allowed_backends = Backend::CPU; } - static inline void EnableDevice(const int dev = 0) { Get().Setup(dev); } - static inline bool DeviceEnabled() { return Get().ngpu > 0; } - static inline bool DeviceDisabled() { return Get().ngpu == 0; } - static inline bool DeviceHasBeenEnabled() { return Get().ngpu >= 0; } + /// Return true if the Device is enabled. + static inline bool IsEnabled() { return Get().mode == ACCELERATED; } - static inline bool UsingDevice() { return DeviceEnabled() && Get().mode == DEVICE; } - static inline bool UsingHost() { return !UsingDevice(); } + /// The opposite of IsEnabled(). + static inline bool IsDisabled() { return !IsEnabled(); } - static inline void DisableTracking() { Get().isTracking = false; }; - static inline void EnableTracking() { Get().isTracking = true; }; - static inline bool IsTracking() { return Get().isTracking; }; - - static inline bool UsingCuda() { return Get().cuda; } - static inline void UseCuda() { Get().cuda = true; } - - static inline bool UsingOmp() { return Get().omp; } - static inline void UseOmp() { Get().omp = true; } - - static inline bool UsingRaja() { return Get().raja; } - static inline void UseRaja() { Get().raja = true; } - - static inline bool UsingOcca() { return Get().occa; } - static inline void UseOcca() { Get().occa = true; } - - static inline bool UsingOkina() - { - return DeviceEnabled() && Get().mode == DEVICE && - (UsingCuda() || UsingOmp() || UsingRaja() || UsingOcca()); - } + /** @brief Return true if any of the backends in the backend mask, @a b_mask, + are allowed. The allowed backends are all configured backends minus the + device backends when the Device is disabled. */ + /** This method can be used with any of the Backend::Id constants, the + Backend::*_MASK, or combinations of those. */ + static inline bool Allows(unsigned long b_mask) + { return Get().allowed_backends & b_mask; } ~Device(); }; diff --git a/general/okina.hpp b/general/forall.hpp similarity index 71% rename from general/okina.hpp rename to general/forall.hpp index 993a5db337..24064edbe4 100644 --- a/general/okina.hpp +++ b/general/forall.hpp @@ -9,16 +9,11 @@ // terms of the GNU Lesser General Public License (as published by the Free // Software Foundation) version 2.1 dated February 1999. -#ifndef MFEM_OKINA_HPP -#define MFEM_OKINA_HPP +#ifndef MFEM_FORALL_HPP +#define MFEM_FORALL_HPP #include "../config/config.hpp" #include "error.hpp" - -#include -#include -#include - #include "cuda.hpp" #include "occa.hpp" #include "device.hpp" @@ -32,16 +27,15 @@ namespace mfem { -// OKINA = Okina Kernel Interface for Numerical Analysis - -// Implementation of MFEM's okina device kernel interface and its CUDA, OpenMP, -// RAJA, and sequential backends. +// Implementation of MFEM's "parallel for" (forall) device/host kernel +// interfaces supporting RAJA, CUDA, OpenMP, and sequential backends. // The MFEM_FORALL wrapper -#define MFEM_FORALL(i,N,...) \ - OkinaWrap(N, \ - [=] MFEM_DEVICE (int i) {__VA_ARGS__}, \ - [&] (int i) {__VA_ARGS__}) +#define MFEM_FORALL(i,N,...) \ + ForallWrap(N, \ + [=] MFEM_ATTR_DEVICE (int i) {__VA_ARGS__}, \ + [&] (int i) {__VA_ARGS__}) + /// OpenMP backend template @@ -49,7 +43,7 @@ void OmpWrap(const int N, HBODY &&h_body) { #ifdef MFEM_USE_OPENMP #pragma omp parallel for - for (int k=0; k void RajaCudaWrap(const int N, DBODY &&d_body) @@ -69,6 +64,7 @@ void RajaCudaWrap(const int N, DBODY &&d_body) #endif } + /// RAJA OpenMP backend template void RajaOmpWrap(const int N, HBODY &&h_body) @@ -80,6 +76,7 @@ void RajaOmpWrap(const int N, HBODY &&h_body) #endif } + /// RAJA sequential loop backend template void RajaSeqWrap(const int N, HBODY &&h_body) @@ -91,8 +88,10 @@ void RajaSeqWrap(const int N, HBODY &&h_body) #endif } + /// CUDA backend #ifdef MFEM_USE_CUDA + template __global__ static void CuKernel(const int N, BODY body) { @@ -100,6 +99,7 @@ void CuKernel(const int N, BODY body) if (k >= N) { return; } body(k); } + template void CuWrap(const int N, DBODY &&d_body) { @@ -109,28 +109,34 @@ void CuWrap(const int N, DBODY &&d_body) const cudaError_t last = cudaGetLastError(); MFEM_VERIFY(last == cudaSuccess, cudaGetErrorString(last)); } -#else + +#else // MFEM_USE_CUDA + template void CuWrap(const int N, DBODY &&d_body) {} + #endif -#define MFEM_CUDA_BLOCKS 256 -/// The okina kernel body wrapper +/// The forall kernel body wrapper template -void OkinaWrap(const int N, DBODY &&d_body, HBODY &&h_body) +void ForallWrap(const int N, DBODY &&d_body, HBODY &&h_body) { - const bool omp = Device::UsingOmp(); - const bool gpu = Device::UsingDevice(); - const bool raja = Device::UsingRaja(); - if (gpu && raja) { return RajaCudaWrap(N, d_body); } - if (gpu) { return CuWrap(N, d_body); } - if (omp && raja) { return RajaOmpWrap(N, h_body); } - if (raja) { return RajaSeqWrap(N, h_body); } - if (omp) { return OmpWrap(N, h_body); } - for (int k=0; k(N, d_body); } + + if (Device::Allows(Backend::CUDA)) + { return CuWrap(N, d_body); } + + if (Device::Allows(Backend::RAJA_OMP)) { return RajaOmpWrap(N, h_body); } + + if (Device::Allows(Backend::OMP)) { return OmpWrap(N, h_body); } + + if (Device::Allows(Backend::RAJA_CPU)) { return RajaSeqWrap(N, h_body); } + + for (int k = 0; k < N; k++) { h_body(k); } } } // namespace mfem -#endif // MFEM_OKINA_HPP +#endif // MFEM_FORALL_HPP diff --git a/general/mem_manager.cpp b/general/mem_manager.cpp index 885eb00689..fab9b128e7 100644 --- a/general/mem_manager.cpp +++ b/general/mem_manager.cpp @@ -9,7 +9,9 @@ // terms of the GNU Lesser General Public License (as published by the Free // Software Foundation) version 2.1 dated February 1999. -#include "../general/okina.hpp" +#include "../general/forall.hpp" + +#include // std::memcpy #include #include @@ -74,7 +76,7 @@ MemoryManager::~MemoryManager() void* MemoryManager::Insert(void *ptr, const std::size_t bytes) { - if (!Device::UsingMM()) { return ptr; } + if (!UsingMM()) { return ptr; } const bool known = IsKnown(ptr); if (known) { @@ -86,7 +88,7 @@ void* MemoryManager::Insert(void *ptr, const std::size_t bytes) void *MemoryManager::Erase(void *ptr) { - if (!Device::UsingMM()) { return ptr; } + if (!UsingMM()) { return ptr; } if (!ptr) { return ptr; } const bool known = IsKnown(ptr); if (!known) @@ -169,10 +171,10 @@ bool MemoryManager::IsAlias(const void *ptr) static inline bool MmDeviceIniFilter(void) { - if (!Device::UsingMM()) { return true; } - if (Device::DeviceDisabled()) { return true; } - if (Device::IsTracking() == false) { return true; } - if (!Device::DeviceHasBeenEnabled()) { return true; } + if (!mm.UsingMM()) { return true; } + if (!mm.IsEnabled()) { return true; } + if (!Device::IsAvailable()) { return true; } + if (!Device::IsConfigured()) { return true; } return false; } @@ -184,7 +186,7 @@ static void *PtrKnown(internal::Ledger *maps, void *ptr) const bool host = base.host; const bool device = !host; const std::size_t bytes = base.bytes; - const bool gpu = Device::UsingDevice(); + const bool gpu = Device::Allows(Backend::DEVICE_MASK); if (host && !gpu) { return ptr; } if (bytes==0) { mfem_error("PtrKnown bytes==0"); } if (!base.d_ptr) { CuMemAlloc(&base.d_ptr, bytes); } @@ -208,7 +210,7 @@ static void *PtrKnown(internal::Ledger *maps, void *ptr) // if necessary. static void *PtrAlias(internal::Ledger *maps, void *ptr) { - const bool gpu = Device::UsingDevice(); + const bool gpu = Device::Allows(Backend::DEVICE_MASK); const internal::Alias *alias = maps->aliases.at(ptr); const internal::Memory *base = alias->mem; const bool host = base->host; @@ -240,7 +242,7 @@ void *MemoryManager::Ptr(void *ptr) if (ptr==NULL) { return NULL; }; if (IsKnown(ptr)) { return PtrKnown(maps, ptr); } if (IsAlias(ptr)) { return PtrAlias(maps, ptr); } - if (Device::UsingDevice()) + if (Device::Allows(Backend::DEVICE_MASK)) { mfem_error("Trying to use unknown pointer on the DEVICE!"); } @@ -273,7 +275,8 @@ void MemoryManager::Push(const void *ptr, const std::size_t bytes) if (MmDeviceIniFilter()) { return; } if (IsKnown(ptr)) { return PushKnown(maps, ptr, bytes); } if (IsAlias(ptr)) { return PushAlias(maps, ptr, bytes); } - if (Device::UsingDevice()) { mfem_error("Unknown pointer to push to!"); } + if (Device::Allows(Backend::DEVICE_MASK)) + { mfem_error("Unknown pointer to push to!"); } } static void PullKnown(const internal::Ledger *maps, @@ -303,25 +306,26 @@ void MemoryManager::Pull(const void *ptr, const std::size_t bytes) if (MmDeviceIniFilter()) { return; } if (IsKnown(ptr)) { return PullKnown(maps, ptr, bytes); } if (IsAlias(ptr)) { return PullAlias(maps, ptr, bytes); } - if (Device::UsingDevice()) { mfem_error("Unknown pointer to pull from!"); } + if (Device::Allows(Backend::DEVICE_MASK)) + { mfem_error("Unknown pointer to pull from!"); } } -extern CUstream *cuStream; +namespace internal { extern CUstream *cuStream; } void* MemoryManager::Memcpy(void *dst, const void *src, const std::size_t bytes, const bool async) { void *d_dst = Ptr(dst); void *d_src = const_cast(Ptr(src)); - const bool host = Device::UsingHost(); if (bytes == 0) { return dst; } + const bool host = !Device::Allows(Backend::DEVICE_MASK); if (host) { return std::memcpy(dst, src, bytes); } if (!async) { return CuMemcpyDtoD(d_dst, d_src, bytes); } - return CuMemcpyDtoDAsync(d_dst, d_src, bytes, cuStream); + return CuMemcpyDtoDAsync(d_dst, d_src, bytes, internal::cuStream); } void MemoryManager::RegisterCheck(void *ptr) { - if (ptr != NULL && Device::UsingMM()) + if (ptr != NULL && UsingMM()) { if (!IsKnown(ptr)) { diff --git a/general/mem_manager.hpp b/general/mem_manager.hpp index 9379a63bab..806708efa8 100644 --- a/general/mem_manager.hpp +++ b/general/mem_manager.hpp @@ -9,21 +9,23 @@ // terms of the GNU Lesser General Public License (as published by the Free // Software Foundation) version 2.1 dated February 1999. -#ifndef MFEM_MEM_MANAGER -#define MFEM_MEM_MANAGER +#ifndef MFEM_MEM_MANAGER_HPP +#define MFEM_MEM_MANAGER_HPP #include "globals.hpp" - namespace mfem { -// Implementation of MFEM's lightweight host/device memory manager designed -// to work seamlessly with the okina device kernel interface. +// Implementation of MFEM's lightweight device/host memory manager designed to +// work seamlessly with the OCCA, RAJA, and other kernels supported by MFEM. /// The memory manager class class MemoryManager { +private: + bool enabled; + public: MemoryManager(); ~MemoryManager(); @@ -34,7 +36,32 @@ public: /// Remove the address from the map, as well as all the address' aliases void *Erase(void *ptr); - /// Return a host or device address, corresponding to the Device::mode + /// Return true if the memory manager is used: pointers seen by mm::New and + /// mm::Delete will be inserted in the ledger and erased from it + static inline bool UsingMM() + { +#ifdef MFEM_USE_MM + return true; +#else + return false; +#endif + } + + /// Disable the memory manager: mm::ptr, mm::push and mm::pull will be no-op + void Disable() { enabled = false; } + + /// Enable the memory manager: mm::ptr, mm::push and mm::pull wont be no-op + void Enable() { enabled = true; } + + /// Return true if the memory manager is used and enabled + bool IsEnabled() { return UsingMM() && enabled; } + + /// The opposite of IsEnabled(). + bool IsDisabled() { return !IsEnabled(); } + + /** @brief Translates ptr to host or device address, depending on what + backends are currently allowed by the Device class and on the ptr + state. */ void *Ptr(void *ptr); const void *Ptr(const void *ptr); @@ -148,4 +175,4 @@ inline void Pull(const void *ptr, const std::size_t bytes = 0) } // namespace mfem -#endif // MFEM_MEM_MANAGER +#endif // MFEM_MEM_MANAGER_HPP diff --git a/general/occa.cpp b/general/occa.cpp index 080e4b43f5..6ddd37d745 100644 --- a/general/occa.cpp +++ b/general/occa.cpp @@ -9,56 +9,53 @@ // terms of the GNU Lesser General Public License (as published by the Free // Software Foundation) version 2.1 dated February 1999. -#include "okina.hpp" +#include "forall.hpp" namespace mfem { -extern OccaDevice occaDevice; +// This variable is defined in device.cpp: +namespace internal { extern OccaDevice occaDevice; } static OccaMemory OccaWrapMemory(const OccaDevice dev, const void *d_adrs, const size_t bytes) { -#if defined(MFEM_USE_OCCA) && defined(MFEM_USE_CUDA) + // This function is called when an OCCA kernel is going to be used. +#ifdef MFEM_USE_OCCA void *adrs = const_cast(d_adrs); - // OCCA & UsingCuda => occa::cuda - if (Device::UsingCuda()) +#if defined(MFEM_USE_CUDA) && OCCA_CUDA_ENABLED + // If OCCA_CUDA is allowed, it will be used since it has the highest priority + if (Device::Allows(Backend::OCCA_CUDA)) { return occa::cuda::wrapMemory(dev, adrs, bytes); } +#endif // MFEM_USE_CUDA && OCCA_CUDA_ENABLED // otherwise, fallback to occa::cpu address space return occa::cpu::wrapMemory(dev, adrs, bytes); -#else // MFEM_USE_OCCA && MFEM_USE_CUDA -#ifdef MFEM_USE_OCCA - return occa::cpu::wrapMemory(dev, const_cast(d_adrs), bytes); -#else +#else // MFEM_USE_OCCA return (void*)NULL; #endif -#endif } OccaMemory OccaPtr(const void *ptr) { - OccaDevice dev = occaDevice; - if (!Device::UsingMM()) { return OccaWrapMemory(dev, ptr, 0); } + // This function is called when 'ptr' needs to be passed to an OCCA kernel. + OccaDevice dev = internal::occaDevice; + if (!mm.UsingMM()) { return OccaWrapMemory(dev, ptr, 0); } const bool known = mm.IsKnown(ptr); if (!known) { mfem_error("OccaPtr: Unknown address!"); } - const bool host = mm.IsOnHost(ptr); + const bool ptr_on_host = mm.IsOnHost(ptr); const size_t bytes = mm.Bytes(ptr); - const bool gpu = Device::UsingDevice(); - if (host && !gpu) { return OccaWrapMemory(dev, ptr, bytes); } - if (!gpu) { mfem_error("OccaPtr: !gpu"); } + const bool run_on_host = !Device::Allows(Backend::DEVICE_MASK); + // If the priority of a host OCCA backend is higher than all device OCCA + // backends, then we will need to run-on-host even if the Device allows a + // device backend. + if (ptr_on_host && run_on_host) { return OccaWrapMemory(dev, ptr, bytes); } + if (run_on_host) { mfem_error("OccaPtr: !ptr_on_host && run_on_host"); } void *d_ptr = mm.GetDevicePtr(ptr); return OccaWrapMemory(dev, d_ptr, bytes); } -OccaDevice OccaWrapDevice(CUdevice dev, CUcontext ctx) -{ -#if defined(MFEM_USE_OCCA) && defined(MFEM_USE_CUDA) - return occa::cuda::wrapDevice(dev, ctx); -#else - return 0; -#endif -} +OccaDevice OccaDev() { return internal::occaDevice; } } // namespace mfem diff --git a/general/occa.hpp b/general/occa.hpp index 3804f37852..4ef13ae718 100644 --- a/general/occa.hpp +++ b/general/occa.hpp @@ -12,20 +12,19 @@ #ifndef MFEM_OCCA_HPP #define MFEM_OCCA_HPP +#include "../config/config.hpp" +#include "cuda.hpp" // for CUdevice, CUcontext + #ifdef MFEM_USE_OCCA #include -#ifdef MFEM_USE_CUDA +#if defined(MFEM_USE_CUDA) && OCCA_CUDA_ENABLED #include #endif typedef occa::device OccaDevice; typedef occa::memory OccaMemory; -#define MFEM_NEW_OCCA_KERNEL(ker, filepath, prop) \ - static occa::kernel ker = NULL; \ - if (ker==NULL) { \ - ker = occaDevice.buildKernel("occa://mfem/" filepath, #ker, prop); \ - } + #else // MFEM_USE_OCCA typedef void* OccaDevice; @@ -36,10 +35,10 @@ typedef void* OccaMemory; namespace mfem { -extern OccaDevice occaDevice; +// Function called when the pointer 'a' needs to be passed to an OCCA kernel. OccaMemory OccaPtr(const void *a); -OccaDevice OccaWrapDevice(CUdevice dev, CUcontext ctx); +OccaDevice OccaDev(); -} // mfem +} // namespace mfem #endif // MFEM_OCCA_HPP diff --git a/general/table.cpp b/general/table.cpp index b9dff3a54b..a2df7d7be0 100644 --- a/general/table.cpp +++ b/general/table.cpp @@ -15,7 +15,7 @@ #include "table.hpp" #include "error.hpp" -#include "../general/okina.hpp" +#include "../general/forall.hpp" #include #include diff --git a/linalg/densemat.cpp b/linalg/densemat.cpp index 71878aef87..8ba013b954 100644 --- a/linalg/densemat.cpp +++ b/linalg/densemat.cpp @@ -17,7 +17,7 @@ #include "matrix.hpp" #include "densemat.hpp" #include "dtensor.hpp" -#include "../general/okina.hpp" +#include "../general/forall.hpp" #include "../general/table.hpp" #include "../general/globals.hpp" @@ -133,7 +133,7 @@ static void Transpose(const int height, const int width, const DeviceVector d_mdata(mdata); MFEM_FORALL(i, height, { - for (int j=0; j class TensorInd { public: - MFEM_HOST_DEVICE static inline int result(const int* sizes, T first, - Args... args) + MFEM_ATTR_HOST_DEVICE + static inline int result(const int* sizes, T first, Args... args) { #ifndef MFEM_USE_CUDA MFEM_ASSERT(first::result(sizes, - args...); + return first + sizes[N - 1] * TensorInd < N + 1, Dim, Args... > + ::result(sizes, args...); } }; + // Terminal case template class TensorInd { public: - MFEM_HOST_DEVICE static inline int result(const int* sizes, T first, - Args... args) + MFEM_ATTR_HOST_DEVICE + static inline int result(const int* sizes, T first, Args... args) { #ifndef MFEM_USE_CUDA MFEM_ASSERT(first class Init @@ -58,6 +61,7 @@ public: return first * Init < N + 1, Dim, Args... >::result(sizes, args...); } }; + // Terminal case template class Init @@ -70,11 +74,11 @@ public: } }; + /// A basic generic Tensor class, appropriate for use on the GPU template class DeviceTensor { - protected: int capacity; Scalar *data; @@ -82,8 +86,7 @@ protected: public: /// Default constructor - explicit DeviceTensor() : capacity(0), data(NULL) - { mfem_error("No default constructor."); } + DeviceTensor() = delete; /// Constructor to initialize a tensor from the Scalar array _data template @@ -117,7 +120,7 @@ public: } /// Copy constructor - MFEM_HOST_DEVICE DeviceTensor(const DeviceTensor& t) + MFEM_ATTR_HOST_DEVICE DeviceTensor(const DeviceTensor& t) { for (int i = 0; i < Dim; ++i) { @@ -130,14 +133,14 @@ public: inline operator Scalar *() const { return data; } /// Const accessor for the data - template MFEM_HOST_DEVICE inline + template MFEM_ATTR_HOST_DEVICE inline Scalar& operator()(Args... args) const { static_assert(sizeof...(args) == Dim, "Wrong number of arguments"); return data[ TensorInd<1, Dim, Args...>::result(sizes, args...) ]; } - MFEM_HOST_DEVICE inline Scalar& operator[](int i) const + MFEM_ATTR_HOST_DEVICE inline Scalar& operator[](int i) const { return data[i]; } diff --git a/linalg/operator.cpp b/linalg/operator.cpp index d1047693ac..e541eaad06 100644 --- a/linalg/operator.cpp +++ b/linalg/operator.cpp @@ -12,7 +12,7 @@ #include "vector.hpp" #include "dtensor.hpp" #include "operator.hpp" -#include "../general/okina.hpp" +#include "../general/forall.hpp" #include #include diff --git a/linalg/sparsemat.cpp b/linalg/sparsemat.cpp index b6ef629061..9d9ba82542 100644 --- a/linalg/sparsemat.cpp +++ b/linalg/sparsemat.cpp @@ -13,7 +13,7 @@ #include "linalg.hpp" #include "dtensor.hpp" -#include "../general/okina.hpp" +#include "../general/forall.hpp" #include "../general/table.hpp" #include "../general/sort_pairs.hpp" @@ -594,7 +594,7 @@ void SparseMatrix::AddMult(const Vector &x, Vector &y, const double a) const { double d = 0.0; const int end = d_I[i+1]; - for (int j=d_I[i]; j < end; j+=1) + for (int j=d_I[i]; j < end; j++) { d += d_A[j] * d_x[d_J[j]]; } diff --git a/linalg/vector.cpp b/linalg/vector.cpp index 241520df41..0c34ff4391 100644 --- a/linalg/vector.cpp +++ b/linalg/vector.cpp @@ -13,7 +13,7 @@ #include "vector.hpp" #include "dtensor.hpp" -#include "../general/okina.hpp" +#include "../general/forall.hpp" #if defined(MFEM_USE_SUNDIALS) && defined(MFEM_USE_MPI) #include @@ -860,10 +860,10 @@ static double cuVectorMin(const int N, const double *X) static CUdeviceptr gdsr = (CUdeviceptr) NULL; if (!gdsr) { ::cuMemAlloc(&gdsr,bytes); } cuKernelMin<<>>(N, (double*)gdsr, x); - CuCheck(cudaGetLastError()); + MFEM_CUDA_CHECK_RT(cudaGetLastError()); ::cuMemcpy((CUdeviceptr)h_min,(CUdeviceptr)gdsr,bytes); double min = std::numeric_limits::infinity(); - for (int i=0; i>>(N, (double*)gdsr, x, y); - CuCheck(cudaGetLastError()); - CuCheck(::cuMemcpy((CUdeviceptr)h_dot,(CUdeviceptr)gdsr,bytes)); + MFEM_CUDA_CHECK_RT(cudaGetLastError()); + MFEM_CUDA_CHECK_DRV(::cuMemcpy((CUdeviceptr)h_dot,(CUdeviceptr)gdsr,bytes)); double dot = 0.0; - for (int i=0; i::infinity(); - for (int i=0; i