Kokkos build and numerics support

include/base/libmesh_common.h

@@ -30,6 +30,10 @@
 // The library configuration options
 #include "libmesh/libmesh_config.h"
 
+// Device compilation support — must be included before assert macros
+// so that LIBMESH_DEVICE_ASSERT is available for the Kokkos path.
+#include "libmesh/libmesh_device.h"
+
 // Use actual timestamps or constant dummies (to aid ccache)
 #ifdef LIBMESH_ENABLE_TIMESTAMPS
 #  define  LIBMESH_TIME __TIME__
@@ -183,33 +187,33 @@ typedef std::complex<Real> COMPLEX;
 
 // Helper functions for complex/real numbers
 // to clean up #ifdef LIBMESH_USE_COMPLEX_NUMBERS elsewhere
-template<typename T> inline T libmesh_real(T a) { return a; }
-template<typename T> inline T libmesh_imag(T /*a*/) { return 0; }
-template<typename T> inline T libmesh_conj(T a) { return a; }
+template<typename T> LIBMESH_DEVICE_INLINE T libmesh_real(T a) { return a; }
+template<typename T> LIBMESH_DEVICE_INLINE T libmesh_imag(T /*a*/) { return 0; }
+template<typename T> LIBMESH_DEVICE_INLINE T libmesh_conj(T a) { return a; }
 
 template<typename T>
-inline T libmesh_real(std::complex<T> a) { return std::real(a); }
+LIBMESH_DEVICE_INLINE T libmesh_real(std::complex<T> a) { return std::real(a); }
 
 template<typename T>
-inline T libmesh_imag(std::complex<T> a) { return std::imag(a); }
+LIBMESH_DEVICE_INLINE T libmesh_imag(std::complex<T> a) { return std::imag(a); }
 
 template<typename T>
-inline std::complex<T> libmesh_conj(std::complex<T> a) { return std::conj(a); }
+LIBMESH_DEVICE_INLINE std::complex<T> libmesh_conj(std::complex<T> a) { return std::conj(a); }
 
 // std::isnan() is in <cmath> as of C++11.
 template <typename T>
-inline bool libmesh_isnan(T x) { return std::isnan(x); }
+LIBMESH_DEVICE_INLINE bool libmesh_isnan(T x) { return std::isnan(x); }
 
 template <typename T>
-inline bool libmesh_isnan(std::complex<T> a)
+LIBMESH_DEVICE_INLINE bool libmesh_isnan(std::complex<T> a)
 { return (std::isnan(std::real(a)) || std::isnan(std::imag(a))); }
 
 // std::isinf() is in <cmath> as of C++11.
 template <typename T>
-inline bool libmesh_isinf(T x) { return std::isinf(x); }
+LIBMESH_DEVICE_INLINE bool libmesh_isinf(T x) { return std::isinf(x); }
 
 template <typename T>
-inline bool libmesh_isinf(std::complex<T> a)
+LIBMESH_DEVICE_INLINE bool libmesh_isinf(std::complex<T> a)
 { return (std::isinf(std::real(a)) || std::isinf(std::imag(a))); }
 
 // Define the value type for unknowns in simulations.
@@ -287,7 +291,13 @@ extern bool warned_about_auto_ptr;
 #endif
 
 // The libmesh_assert() macro acts like C's assert(), but throws a
-// libmesh_error() (including stack trace, etc) instead of just exiting
+// libmesh_error() (including stack trace, etc) instead of just exiting.
+//
+// In .K translation units (LIBMESH_KOKKOS_COMPILATION defined),
+// LIBMESH_DEVICE_ASSERT is provided by libmesh_device.h using
+// printf + Kokkos::abort() — device-safe across CUDA/HIP/SYCL.
+// The assert macros delegate to it so that both host and device
+// code in the same file get assertion checking.
 #ifdef NDEBUG
 
 #define libmesh_assert_msg(asserted, msg)  ((void) 0)
@@ -299,6 +309,18 @@ extern bool warned_about_auto_ptr;
 #define libmesh_assert_less_equal_msg(expr1,expr2, msg)  ((void) 0)
 #define libmesh_assert_greater_equal_msg(expr1,expr2, msg)  ((void) 0)
 
+#elif defined(LIBMESH_DEVICE_ASSERT)
+
+// Kokkos compilation: use the device-safe assert from libmesh_device.h.
+#define libmesh_assert_msg(asserted, msg)  LIBMESH_DEVICE_ASSERT(asserted)
+#define libmesh_exceptionless_assert_msg(asserted, msg) LIBMESH_DEVICE_ASSERT(asserted)
+#define libmesh_assert_equal_to_msg(expr1,expr2, msg)  LIBMESH_DEVICE_ASSERT((expr1) == (expr2))
+#define libmesh_assert_not_equal_to_msg(expr1,expr2, msg)  LIBMESH_DEVICE_ASSERT((expr1) != (expr2))
+#define libmesh_assert_less_msg(expr1,expr2, msg)  LIBMESH_DEVICE_ASSERT((expr1) < (expr2))
+#define libmesh_assert_greater_msg(expr1,expr2, msg)  LIBMESH_DEVICE_ASSERT((expr1) > (expr2))
+#define libmesh_assert_less_equal_msg(expr1,expr2, msg)  LIBMESH_DEVICE_ASSERT((expr1) <= (expr2))
+#define libmesh_assert_greater_equal_msg(expr1,expr2, msg)  LIBMESH_DEVICE_ASSERT((expr1) >= (expr2))
+
 #else
 
 #define libmesh_assertion_types(expr1,expr2)                            \
@@ -674,7 +696,6 @@ inline Tnew restrict_int (Told oldvar)
   return oldvar;
 }
 
-
 /**
  * This is a helper variable template for cases when we want to use a default compile-time
  * error with constexpr-based if conditions. The templating delays the triggering

include/base/libmesh_device.h

@@ -0,0 +1,74 @@
+// The libMesh Finite Element Library.
+// Copyright (C) 2002-2026 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
+
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+// Lesser General Public License for more details.
+
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+
+#ifndef LIBMESH_LIBMESH_DEVICE_H
+#define LIBMESH_LIBMESH_DEVICE_H
+
+// Defines LIBMESH_DEVICE_INLINE, mirroring MetaPhysicL's METAPHYSICL_INLINE
+// pattern (metaphysicl_device.h / METAPHYSICL_KOKKOS_COMPILATION).
+//
+// When compiling a .K translation unit (LIBMESH_KOKKOS_COMPILATION is defined
+// by kokkos.mk), this expands to KOKKOS_INLINE_FUNCTION so that annotated
+// methods are callable from both host and device code.  In all other
+// translation units it expands to plain `inline`.
+#ifdef LIBMESH_KOKKOS_COMPILATION
+#  include <Kokkos_Macros.hpp>
+#  include <Kokkos_Abort.hpp>
+#  define LIBMESH_DEVICE_INLINE KOKKOS_INLINE_FUNCTION
+
+// Backend-neutral device-code detection for Kokkos .K translation units.
+// This lets error/exception plumbing share a single predicate instead of
+// hardcoding per-backend checks in multiple headers.
+#  if defined(__CUDA_ARCH__) || defined(__HIP_DEVICE_COMPILE__) || defined(__SYCL_DEVICE_ONLY__)
+#    define LIBMESH_IN_DEVICE_CODE 1
+#  else
+#    define LIBMESH_IN_DEVICE_CODE 0
+#  endif
+
+// Device-safe assert: uses printf (supported on CUDA/HIP) and
+// Kokkos::abort() for backend-portable device termination.
+// Defined here (not in libmesh_common.h) because Kokkos headers
+// are only available in .K translation units.
+#  ifndef NDEBUG
+#    define LIBMESH_DEVICE_ASSERT(asserted)                             \
+       do { if (!(asserted)) {                                          \
+         printf("libMesh assert failed: %s, file %s, line %d\n",       \
+                #asserted, __FILE__, __LINE__);                         \
+         ::Kokkos::abort("libmesh_assert failed");                     \
+       } } while (0)
+#  else
+#    define LIBMESH_DEVICE_ASSERT(asserted) ((void) 0)
+#  endif
+
+#  define LIBMESH_DEVICE_ERROR_MSG(msg)                                 \
+     do {                                                               \
+       printf("libMesh error: %s, file %s, line %d\n",                 \
+              msg, __FILE__, __LINE__);                                 \
+       ::Kokkos::abort(msg);                                            \
+     } while (0)
+
+#  define LIBMESH_DEVICE_ERROR_MSG_IF(cond, msg)                        \
+     do { if (cond) { LIBMESH_DEVICE_ERROR_MSG(msg); } } while (0)
+
+#else
+#  define LIBMESH_DEVICE_INLINE inline
+#  define LIBMESH_IN_DEVICE_CODE 0
+#  define LIBMESH_DEVICE_ERROR_MSG(msg) libmesh_error_msg(msg)
+#  define LIBMESH_DEVICE_ERROR_MSG_IF(cond, msg) libmesh_error_msg_if(cond, msg)
+#endif
+
+#endif // LIBMESH_LIBMESH_DEVICE_H

include/base/libmesh_exceptions.h

@@ -23,6 +23,7 @@
 #include "libmesh/libmesh_config.h"
 
 #include "libmesh/libmesh_abort.h"
+#include "libmesh/libmesh_device.h"
 
 #include <stdexcept>
 #include <string>
@@ -212,7 +213,13 @@ class TerminationException
 #ifdef LIBMESH_ENABLE_EXCEPTIONS
 #define libmesh_noexcept noexcept
 
+#if LIBMESH_IN_DEVICE_CODE
+// Kokkos device code does not support C++ exceptions.
+#define LIBMESH_THROW(e) do { LIBMESH_DEVICE_ERROR_MSG((e).what()); } while (0)
+#else
 #define LIBMESH_THROW(e) do { throw e; } while (0)
+#endif
+
 #define libmesh_rethrow throw
 #define libmesh_try try
 #define libmesh_catch(e) catch(e)

include/libmesh_config.h.in

@@ -434,6 +434,9 @@
 /* Define to 1 if you have the <inttypes.h> header file. */
 #undef HAVE_INTTYPES_H
 
+/* Define if Kokkos support is enabled in libMesh */
+#undef HAVE_KOKKOS
+
 /* Flag indicating whether the library will be compiled with LASPACK support
    */
 #undef HAVE_LASPACK

include/numerics/tensor_tools.h

@@ -45,92 +45,92 @@ namespace TensorTools
 // Vector specializations will follow.
 
 template <typename T, typename T2>
-inline
+LIBMESH_DEVICE_INLINE
 typename std::enable_if<ScalarTraits<T>::value && ScalarTraits<T2>::value,
                         typename CompareTypes<T, T2>::supertype>::type
 inner_product(const T & a, const T2& b)
 { return a * b; }
 
 template <typename T, typename T2>
-inline
+LIBMESH_DEVICE_INLINE
 typename CompareTypes<T, T2>::supertype
 inner_product(const TypeVector<T> & a, const TypeVector<T2> & b)
 { return a * b; }
 
 template <typename T, typename T2>
-inline
+LIBMESH_DEVICE_INLINE
 typename CompareTypes<T, T2>::supertype
 inner_product(const TypeTensor<T> & a, const TypeTensor<T2> & b)
 { return a.contract(b); }
 
 template <unsigned int N, typename T, typename T2>
-inline
+LIBMESH_DEVICE_INLINE
 typename CompareTypes<T, T2>::supertype
 inner_product(const TypeNTensor<N,T> & a, const TypeNTensor<N,T2> & b)
 { return a.contract(b); }
 
 template<typename T>
-inline
+LIBMESH_DEVICE_INLINE
 auto norm(const T & a)
 { using std::abs; return abs(a); }
 
 template<typename T>
-inline
+LIBMESH_DEVICE_INLINE
 T norm(std::complex<T> a) { using std::abs; return abs(a); }
 
 template <typename T>
-inline
+LIBMESH_DEVICE_INLINE
 auto norm(const TypeVector<T> & a) -> decltype(TensorTools::norm(T()))
 {using std::sqrt; return sqrt(a.norm_sq());}
 
 template <typename T>
-inline
+LIBMESH_DEVICE_INLINE
 auto norm(const VectorValue<T> & a) -> decltype(TensorTools::norm(T()))
 {using std::sqrt; return sqrt(a.norm_sq());}
 
 template <typename T>
-inline
+LIBMESH_DEVICE_INLINE
 auto norm(const TypeTensor<T> & a) -> decltype(TensorTools::norm(T()))
 {using std::sqrt; return sqrt(a.norm_sq());}
 
 template <typename T>
-inline
+LIBMESH_DEVICE_INLINE
 auto norm(const TensorValue<T> & a) -> decltype(TensorTools::norm(T()))
 {using std::sqrt; return sqrt(a.norm_sq());}
 
 
 template<typename T>
-inline
+LIBMESH_DEVICE_INLINE
 auto norm_sq(const T & a)
-{ using std::norm; return norm(a); }
+{ return a * libmesh_conj(a); }
 
 template<typename T>
-inline
+LIBMESH_DEVICE_INLINE
 T norm_sq(std::complex<T> a) { using std::norm; return norm(a); }
 
 template <typename T>
-inline
+LIBMESH_DEVICE_INLINE
 auto norm_sq(const TypeVector<T> & a)
 {return a.norm_sq();}
 
 template <typename T>
-inline
+LIBMESH_DEVICE_INLINE
 auto norm_sq(const VectorValue<T> & a)
 {return a.norm_sq();}
 
 template <typename T>
-inline
+LIBMESH_DEVICE_INLINE
 auto norm_sq(const TypeTensor<T> & a)
 {return a.norm_sq();}
 
 template <typename T>
-inline
+LIBMESH_DEVICE_INLINE
 auto norm_sq(const TensorValue<T> & a)
 {return a.norm_sq();}
 
 
 template<typename T>
-inline
+LIBMESH_DEVICE_INLINE
 bool is_zero(const T & a){ return a.is_zero();}
 
 // Any tensor-rank-independent code will need to include