diff --git a/.gitmodules b/.gitmodules
index 1291021..326f2bd 100644
--- a/.gitmodules
+++ b/.gitmodules
@@ -1,9 +1,6 @@
 [submodule "utils/spiner"]
 	path = utils/spiner
 	url = https://github.com/lanl/spiner.git
-[submodule "utils/variant"]
-	path = utils/variant
-	url = https://github.com/mpark/variant.git
 [submodule "utils/herumi-fmath"]
 	path = utils/herumi-fmath
 	url = https://github.com/herumi/fmath.git
diff --git a/CMakeLists.txt b/CMakeLists.txt
index b34b014..79ae778 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -80,13 +80,6 @@ endif()
 
 set(CMAKE_EXPORT_COMPILE_COMMANDS On)
 
-# Patches variant to be compatible with cuda
-# Assumes "patch" is present on system
-message(STATUS "Patching mpark::variant to support GPUs")
-execute_process(COMMAND patch -N -s -V never
-                ${CMAKE_CURRENT_SOURCE_DIR}/utils/variant/include/mpark/variant.hpp
-                ${CMAKE_CURRENT_SOURCE_DIR}/utils/cuda_compatibility.patch)
-
 # xl fix
 target_compile_options(singularity-opac::flags INTERFACE
                        $<$<COMPILE_LANG_AND_ID:CXX,XL>:-std=c++1y;-qxflag=disable__cplusplusOverride>)
diff --git a/singularity-opac/base/indexers.hpp b/singularity-opac/base/indexers.hpp
index 4f4657c..c642668 100644
--- a/singularity-opac/base/indexers.hpp
+++ b/singularity-opac/base/indexers.hpp
@@ -18,7 +18,7 @@
 
 #include <fast-math/logs.hpp>
 #include <spiner/databox.hpp>
-#include <variant/include/mpark/variant.hpp>
+#include <ports-of-call/variant.hpp>
 
 #include <singularity-opac/chebyshev/chebyshev.hpp>
 
diff --git a/singularity-opac/neutrinos/mean_neutrino_s_variant.hpp b/singularity-opac/neutrinos/mean_neutrino_s_variant.hpp
index 58ad41b..baeb39a 100644
--- a/singularity-opac/neutrinos/mean_neutrino_s_variant.hpp
+++ b/singularity-opac/neutrinos/mean_neutrino_s_variant.hpp
@@ -19,10 +19,10 @@
 #include <utility>
 
 #include <ports-of-call/portability.hpp>
+#include <ports-of-call/variant.hpp>
 #include <singularity-opac/base/opac_error.hpp>
 #include <singularity-opac/base/radiation_types.hpp>
 #include <singularity-opac/neutrinos/neutrino_s_variant.hpp>
-#include <variant/include/mpark/variant.hpp>
 
 namespace singularity {
 namespace neutrinos {
@@ -60,11 +60,11 @@ class MeanSVariant {
           !std::is_same<MeanSVariant, typename std::decay<Choice>::type>::value,
           bool>::type = true>
   Choice get() {
-    return mpark::get<Choice>(s_opac_);
+    return PortsOfCall::get<Choice>(s_opac_);
   }
 
   MeanSVariant GetOnDevice() {
-    return mpark::visit(
+    return PortsOfCall::visit(
         [](auto &s_opac) {
           return s_opac_variant<SOpacs...>(s_opac.GetOnDevice());
         },
@@ -74,7 +74,7 @@ class MeanSVariant {
   PORTABLE_INLINE_FUNCTION Real PlanckMeanTotalScatteringCoefficient(
       const Real rho, const Real temp, const Real Ye,
       const RadiationType type) const {
-    return mpark::visit(
+    return PortsOfCall::visit(
         [=](const auto &s_opac) {
           return s_opac.PlanckMeanTotalScatteringCoefficient(rho, temp, Ye,
                                                              type);
@@ -84,7 +84,7 @@ class MeanSVariant {
   PORTABLE_INLINE_FUNCTION Real RosselandMeanTotalScatteringCoefficient(
       const Real rho, const Real temp, const Real Ye,
       const RadiationType type) const {
-    return mpark::visit(
+    return PortsOfCall::visit(
         [=](const auto &s_opac) {
           return s_opac.RosselandMeanTotalScatteringCoefficient(rho, temp, Ye,
                                                                 type);
@@ -93,8 +93,8 @@ class MeanSVariant {
   }
 
   inline void Finalize() noexcept {
-    return mpark::visit([](auto &s_opac) { return s_opac.Finalize(); },
-                        s_opac_);
+    return PortsOfCall::visit([](auto &s_opac) { return s_opac.Finalize(); },
+                              s_opac_);
   }
 };
 
diff --git a/singularity-opac/neutrinos/mean_neutrino_variant.hpp b/singularity-opac/neutrinos/mean_neutrino_variant.hpp
index 0675fed..4e9f5a0 100644
--- a/singularity-opac/neutrinos/mean_neutrino_variant.hpp
+++ b/singularity-opac/neutrinos/mean_neutrino_variant.hpp
@@ -19,10 +19,10 @@
 #include <utility>
 
 #include <ports-of-call/portability.hpp>
+#include <ports-of-call/variant.hpp>
 #include <singularity-opac/base/opac_error.hpp>
 #include <singularity-opac/base/radiation_types.hpp>
 #include <singularity-opac/neutrinos/neutrino_variant.hpp>
-#include <variant/include/mpark/variant.hpp>
 
 namespace singularity {
 namespace neutrinos {
@@ -60,25 +60,25 @@ class MeanVariant {
           !std::is_same<MeanVariant, typename std::decay<Choice>::type>::value,
           bool>::type = true>
   Choice get() {
-    return mpark::get<Choice>(opac_);
+    return PortsOfCall::get<Choice>(opac_);
   }
 
   MeanVariant GetOnDevice() {
-    return mpark::visit(
+    return PortsOfCall::visit(
         [](auto &opac) { return opac_variant<Opacs...>(opac.GetOnDevice()); },
         opac_);
   }
 
   PORTABLE_INLINE_FUNCTION RuntimePhysicalConstants
   GetRuntimePhysicalConstants() const {
-    return mpark::visit(
+    return PortsOfCall::visit(
         [](auto &opac) { return opac.GetRuntimePhysicalConstants(); }, opac_);
   }
 
   PORTABLE_INLINE_FUNCTION Real PlanckMeanAbsorptionCoefficient(
       const Real rho, const Real temp, const Real Ye,
       const RadiationType type) const {
-    return mpark::visit(
+    return PortsOfCall::visit(
         [=](const auto &opac) {
           return opac.PlanckMeanAbsorptionCoefficient(rho, temp, Ye, type);
         },
@@ -87,7 +87,7 @@ class MeanVariant {
   PORTABLE_INLINE_FUNCTION Real RosselandMeanAbsorptionCoefficient(
       const Real rho, const Real temp, const Real Ye,
       const RadiationType type) const {
-    return mpark::visit(
+    return PortsOfCall::visit(
         [=](const auto &opac) {
           return opac.RosselandMeanAbsorptionCoefficient(rho, temp, Ye, type);
         },
@@ -95,12 +95,14 @@ class MeanVariant {
   }
 
   inline void Finalize() noexcept {
-    return mpark::visit([](auto &opac) { return opac.Finalize(); }, opac_);
+    return PortsOfCall::visit([](auto &opac) { return opac.Finalize(); },
+                              opac_);
   }
 
 #ifdef SPINER_USE_HDF
   void Save(const std::string &filename) const {
-    return mpark::visit([=](auto &opac) { return opac.Save(filename); }, opac_);
+    return PortsOfCall::visit([=](auto &opac) { return opac.Save(filename); },
+                              opac_);
   }
 #endif
 };
diff --git a/singularity-opac/neutrinos/neutrino_s_variant.hpp b/singularity-opac/neutrinos/neutrino_s_variant.hpp
index feb09f6..b3829de 100644
--- a/singularity-opac/neutrinos/neutrino_s_variant.hpp
+++ b/singularity-opac/neutrinos/neutrino_s_variant.hpp
@@ -19,16 +19,16 @@
 #include <utility>
 
 #include <ports-of-call/portability.hpp>
+#include <ports-of-call/variant.hpp>
 #include <singularity-opac/base/opac_error.hpp>
 #include <singularity-opac/base/radiation_types.hpp>
-#include <variant/include/mpark/variant.hpp>
 
 namespace singularity {
 namespace neutrinos {
 namespace impl {
 
 template <typename... Ts>
-using s_opac_variant = mpark::variant<Ts...>;
+using s_opac_variant = PortsOfCall::variant<Ts...>;
 
 template <typename... S_Opacs>
 class S_Variant {
@@ -62,11 +62,11 @@ class S_Variant {
           !std::is_same<S_Variant, typename std::decay<Choice>::type>::value,
           bool>::type = true>
   Choice get() {
-    return mpark::get<Choice>(s_opac_);
+    return PortsOfCall::get<Choice>(s_opac_);
   }
 
   S_Variant GetOnDevice() {
-    return mpark::visit(
+    return PortsOfCall::visit(
         [](auto &s_opac) {
           return s_opac_variant<S_Opacs...>(s_opac.GetOnDevice());
         },
@@ -79,7 +79,7 @@ class S_Variant {
   PORTABLE_INLINE_FUNCTION Real TotalCrossSection(
       const Real rho, const Real temp, const Real Ye, const RadiationType type,
       const Real nu, Real *lambda = nullptr) const {
-    return mpark::visit(
+    return PortsOfCall::visit(
         [=](const auto &s_opac) {
           return s_opac.TotalCrossSection(rho, temp, Ye, type, nu, lambda);
         },
@@ -92,7 +92,7 @@ class S_Variant {
   PORTABLE_INLINE_FUNCTION Real DifferentialCrossSection(
       const Real rho, const Real temp, const Real Ye, const RadiationType type,
       const Real nu, const Real mu, Real *lambda = nullptr) const {
-    return mpark::visit(
+    return PortsOfCall::visit(
         [=](const auto &s_opac) {
           return s_opac.DifferentialCrossSection(rho, temp, Ye, type, nu, mu,
                                                  lambda);
@@ -106,7 +106,7 @@ class S_Variant {
   PORTABLE_INLINE_FUNCTION Real TotalScatteringCoefficient(
       const Real rho, const Real temp, const Real Ye, const RadiationType type,
       const Real nu, Real *lambda = nullptr) const {
-    return mpark::visit(
+    return PortsOfCall::visit(
         [=](const auto &s_opac) {
           return s_opac.TotalScatteringCoefficient(rho, temp, Ye, type, nu,
                                                    lambda);
@@ -116,24 +116,24 @@ class S_Variant {
 
   PORTABLE_INLINE_FUNCTION
   int nlambda() const noexcept {
-    return mpark::visit([](const auto &s_opac) { return s_opac.nlambda(); },
-                        s_opac_);
+    return PortsOfCall::visit(
+        [](const auto &s_opac) { return s_opac.nlambda(); }, s_opac_);
   }
 
   template <typename T>
   PORTABLE_INLINE_FUNCTION bool IsType() const noexcept {
-    return mpark::holds_alternative<T>(s_opac_);
+    return PortsOfCall::holds_alternative<T>(s_opac_);
   }
 
   PORTABLE_INLINE_FUNCTION
   void PrintParams() const noexcept {
-    return mpark::visit([](const auto &s_opac) { return s_opac.PrintParams(); },
-                        s_opac_);
+    return PortsOfCall::visit(
+        [](const auto &s_opac) { return s_opac.PrintParams(); }, s_opac_);
   }
 
   inline void Finalize() noexcept {
-    return mpark::visit([](auto &s_opac) { return s_opac.Finalize(); },
-                        s_opac_);
+    return PortsOfCall::visit([](auto &s_opac) { return s_opac.Finalize(); },
+                              s_opac_);
   }
 };
 
diff --git a/singularity-opac/neutrinos/neutrino_variant.hpp b/singularity-opac/neutrinos/neutrino_variant.hpp
index b53551a..2fc255a 100644
--- a/singularity-opac/neutrinos/neutrino_variant.hpp
+++ b/singularity-opac/neutrinos/neutrino_variant.hpp
@@ -19,16 +19,16 @@
 #include <utility>
 
 #include <ports-of-call/portability.hpp>
+#include <ports-of-call/variant.hpp>
 #include <singularity-opac/base/opac_error.hpp>
 #include <singularity-opac/base/radiation_types.hpp>
-#include <variant/include/mpark/variant.hpp>
 
 namespace singularity {
 namespace neutrinos {
 namespace impl {
 
 template <typename... Ts>
-using opac_variant = mpark::variant<Ts...>;
+using opac_variant = PortsOfCall::variant<Ts...>;
 
 template <typename... Opacs>
 class Variant {
@@ -62,18 +62,18 @@ class Variant {
           !std::is_same<Variant, typename std::decay<Choice>::type>::value,
           bool>::type = true>
   Choice get() {
-    return mpark::get<Choice>(opac_);
+    return PortsOfCall::get<Choice>(opac_);
   }
 
   Variant GetOnDevice() {
-    return mpark::visit(
+    return PortsOfCall::visit(
         [](auto &opac) { return opac_variant<Opacs...>(opac.GetOnDevice()); },
         opac_);
   }
 
   PORTABLE_INLINE_FUNCTION RuntimePhysicalConstants
   GetRuntimePhysicalConstants() const {
-    return mpark::visit(
+    return PortsOfCall::visit(
         [](auto &opac) { return opac.GetRuntimePhysicalConstants(); }, opac_);
   }
 
@@ -83,7 +83,7 @@ class Variant {
   PORTABLE_INLINE_FUNCTION Real AbsorptionCoefficient(
       const Real rho, const Real temp, const Real Ye, const RadiationType type,
       const Real nu, Real *lambda = nullptr) const {
-    return mpark::visit(
+    return PortsOfCall::visit(
         [=](const auto &opac) {
           return opac.AbsorptionCoefficient(rho, temp, Ye, type, nu, lambda);
         },
@@ -96,7 +96,7 @@ class Variant {
                         const RadiationType type, FrequencyIndexer &nu_bins,
                         DataIndexer &coeffs, const int nbins,
                         Real *lambda = nullptr) const {
-    mpark::visit(
+    PortsOfCall::visit(
         [&](const auto &opac) {
           opac.AbsorptionCoefficient(rho, temp, Ye, type, nu_bins, coeffs,
                                      nbins, lambda);
@@ -110,7 +110,7 @@ class Variant {
   PORTABLE_INLINE_FUNCTION Real AngleAveragedAbsorptionCoefficient(
       const Real rho, const Real temp, const Real Ye, const RadiationType type,
       const Real nu, Real *lambda = nullptr) const {
-    return mpark::visit(
+    return PortsOfCall::visit(
         [=](const auto &opac) {
           return opac.AngleAveragedAbsorptionCoefficient(rho, temp, Ye, type,
                                                          nu, lambda);
@@ -123,7 +123,7 @@ class Variant {
       const Real rho, const Real temp, const Real Ye, const RadiationType type,
       FrequencyIndexer &nu_bins, DataIndexer &coeffs, const int nbins,
       Real *lambda = nullptr) const {
-    mpark::visit(
+    PortsOfCall::visit(
         [&](const auto &opac) {
           opac.AngleAveragedAbsorptionCoefficient(rho, temp, Ye, type, nu_bins,
                                                   coeffs, nbins, lambda);
@@ -136,7 +136,7 @@ class Variant {
   PORTABLE_INLINE_FUNCTION Real EmissivityPerNuOmega(
       const Real rho, const Real temp, const Real Ye, const RadiationType type,
       const Real nu, Real *lambda = nullptr) const {
-    return mpark::visit(
+    return PortsOfCall::visit(
         [=](const auto &opac) {
           return opac.EmissivityPerNuOmega(rho, temp, Ye, type, nu, lambda);
         },
@@ -149,7 +149,7 @@ class Variant {
                        const RadiationType type, FrequencyIndexer &nu_bins,
                        DataIndexer &coeffs, const int nbins,
                        Real *lambda = nullptr) const {
-    mpark::visit(
+    PortsOfCall::visit(
         [&](const auto &opac) {
           return opac.EmissivityPerNuOmega(rho, temp, Ye, type, nu_bins, coeffs,
                                            nbins, lambda);
@@ -163,7 +163,7 @@ class Variant {
                                                 const RadiationType type,
                                                 const Real nu,
                                                 Real *lambda = nullptr) const {
-    return mpark::visit(
+    return PortsOfCall::visit(
         [=](const auto &opac) {
           return opac.EmissivityPerNu(rho, temp, Ye, type, nu, lambda);
         },
@@ -176,7 +176,7 @@ class Variant {
                   const RadiationType type, FrequencyIndexer &nu_bins,
                   DataIndexer &coeffs, const int nbins,
                   Real *lambda = nullptr) const {
-    mpark::visit(
+    PortsOfCall::visit(
         [&](const auto &opac) {
           return opac.EmissivityPerNu(rho, temp, Ye, type, nu_bins, coeffs,
                                       nbins, lambda);
@@ -189,7 +189,7 @@ class Variant {
                                            const Real Ye,
                                            const RadiationType type,
                                            Real *lambda = nullptr) const {
-    return mpark::visit(
+    return PortsOfCall::visit(
         [=](const auto &opac) {
           return opac.Emissivity(rho, temp, Ye, type, lambda);
         },
@@ -201,7 +201,7 @@ class Variant {
                                                  const Real temp, Real Ye,
                                                  const RadiationType type,
                                                  Real *lambda = nullptr) const {
-    return mpark::visit(
+    return PortsOfCall::visit(
         [=](const auto &opac) {
           return opac.NumberEmissivity(rho, temp, Ye, type, lambda);
         },
@@ -212,7 +212,7 @@ class Variant {
   PORTABLE_INLINE_FUNCTION Real
   ThermalDistributionOfTNu(const Real temp, const RadiationType type,
                            const Real nu, Real *lambda = nullptr) const {
-    return mpark::visit(
+    return PortsOfCall::visit(
         [=](const auto &opac) {
           return opac.ThermalDistributionOfTNu(temp, type, nu, lambda);
         },
@@ -223,7 +223,7 @@ class Variant {
   PORTABLE_INLINE_FUNCTION Real
   DThermalDistributionOfTNuDT(const Real temp, const RadiationType type,
                               const Real nu, Real *lambda = nullptr) const {
-    return mpark::visit(
+    return PortsOfCall::visit(
         [=](const auto &opac) {
           return opac.DThermalDistributionOfTNuDT(temp, type, nu, lambda);
         },
@@ -233,7 +233,7 @@ class Variant {
   // Integral of thermal distribution over frequency and angle
   PORTABLE_INLINE_FUNCTION Real ThermalDistributionOfT(
       const Real temp, const RadiationType type, Real *lambda = nullptr) const {
-    return mpark::visit(
+    return PortsOfCall::visit(
         [=](const auto &opac) {
           return opac.ThermalDistributionOfT(temp, type, lambda);
         },
@@ -243,7 +243,7 @@ class Variant {
   // Integral of thermal distribution/energy over frequency and angle
   PORTABLE_INLINE_FUNCTION Real ThermalNumberDistributionOfT(
       const Real temp, const RadiationType type, Real *lambda = nullptr) const {
-    return mpark::visit(
+    return PortsOfCall::visit(
         [=](const auto &opac) {
           return opac.ThermalNumberDistributionOfT(temp, type, lambda);
         },
@@ -253,7 +253,7 @@ class Variant {
   // Energy density of thermal distribution
   PORTABLE_INLINE_FUNCTION Real EnergyDensityFromTemperature(
       const Real temp, const RadiationType type, Real *lambda = nullptr) const {
-    return mpark::visit(
+    return PortsOfCall::visit(
         [=](const auto &opac) {
           return opac.EnergyDensityFromTemperature(temp, type, lambda);
         },
@@ -263,7 +263,7 @@ class Variant {
   // Temperature of thermal distribution
   PORTABLE_INLINE_FUNCTION Real TemperatureFromEnergyDensity(
       const Real er, const RadiationType type, Real *lambda = nullptr) const {
-    return mpark::visit(
+    return PortsOfCall::visit(
         [=](const auto &opac) {
           return opac.TemperatureFromEnergyDensity(er, type, lambda);
         },
@@ -273,7 +273,7 @@ class Variant {
   // Number density of thermal distribution
   PORTABLE_INLINE_FUNCTION Real NumberDensityFromTemperature(
       const Real temp, const RadiationType type, Real *lambda = nullptr) const {
-    return mpark::visit(
+    return PortsOfCall::visit(
         [=](const auto &opac) {
           return opac.NumberDensityFromTemperature(temp, type, lambda);
         },
@@ -282,22 +282,24 @@ class Variant {
 
   PORTABLE_INLINE_FUNCTION
   int nlambda() const noexcept {
-    return mpark::visit([](const auto &opac) { return opac.nlambda(); }, opac_);
+    return PortsOfCall::visit([](const auto &opac) { return opac.nlambda(); },
+                              opac_);
   }
 
   template <typename T>
   PORTABLE_INLINE_FUNCTION bool IsType() const noexcept {
-    return mpark::holds_alternative<T>(opac_);
+    return PortsOfCall::holds_alternative<T>(opac_);
   }
 
   PORTABLE_INLINE_FUNCTION
   void PrintParams() const noexcept {
-    return mpark::visit([](const auto &opac) { return opac.PrintParams(); },
-                        opac_);
+    return PortsOfCall::visit(
+        [](const auto &opac) { return opac.PrintParams(); }, opac_);
   }
 
   inline void Finalize() noexcept {
-    return mpark::visit([](auto &opac) { return opac.Finalize(); }, opac_);
+    return PortsOfCall::visit([](auto &opac) { return opac.Finalize(); },
+                              opac_);
   }
 };
 
diff --git a/singularity-opac/neutrinos/opac_neutrinos.hpp b/singularity-opac/neutrinos/opac_neutrinos.hpp
index 290190f..9967168 100644
--- a/singularity-opac/neutrinos/opac_neutrinos.hpp
+++ b/singularity-opac/neutrinos/opac_neutrinos.hpp
@@ -16,7 +16,7 @@
 #ifndef SINGULARITY_OPAC_NEUTRINOS_OPAC_NEUTRINOS_
 #define SINGULARITY_OPAC_NEUTRINOS_OPAC_NEUTRINOS_
 
-#include <variant/include/mpark/variant.hpp>
+#include <ports-of-call/variant.hpp>
 
 #include <singularity-opac/neutrinos/brt_neutrinos.hpp>
 #include <singularity-opac/neutrinos/gray_opacity_neutrinos.hpp>
diff --git a/singularity-opac/neutrinos/s_opac_neutrinos.hpp b/singularity-opac/neutrinos/s_opac_neutrinos.hpp
index b81982c..7e326fb 100644
--- a/singularity-opac/neutrinos/s_opac_neutrinos.hpp
+++ b/singularity-opac/neutrinos/s_opac_neutrinos.hpp
@@ -16,7 +16,7 @@
 #ifndef SINGULARITY_OPAC_NEUTRINOS_S_OPAC_NEUTRINOS_
 #define SINGULARITY_OPAC_NEUTRINOS_S_OPAC_NEUTRINOS_
 
-#include <variant/include/mpark/variant.hpp>
+#include <ports-of-call/variant.hpp>
 
 #include <singularity-opac/neutrinos/gray_s_opacity_neutrinos.hpp>
 #include <singularity-opac/neutrinos/neutrino_s_variant.hpp>
diff --git a/singularity-opac/photons/mean_photon_s_variant.hpp b/singularity-opac/photons/mean_photon_s_variant.hpp
index 515f1f9..699e5e1 100644
--- a/singularity-opac/photons/mean_photon_s_variant.hpp
+++ b/singularity-opac/photons/mean_photon_s_variant.hpp
@@ -19,10 +19,10 @@
 #include <utility>
 
 #include <ports-of-call/portability.hpp>
+#include <ports-of-call/variant.hpp>
 #include <singularity-opac/base/opac_error.hpp>
 #include <singularity-opac/base/radiation_types.hpp>
 #include <singularity-opac/photons/photon_s_variant.hpp>
-#include <variant/include/mpark/variant.hpp>
 
 namespace singularity {
 namespace photons {
@@ -60,11 +60,11 @@ class MeanSVariant {
           !std::is_same<MeanSVariant, typename std::decay<Choice>::type>::value,
           bool>::type = true>
   Choice get() {
-    return mpark::get<Choice>(s_opac_);
+    return PortsOfCall::get<Choice>(s_opac_);
   }
 
   MeanSVariant GetOnDevice() {
-    return mpark::visit(
+    return PortsOfCall::visit(
         [](auto &s_opac) {
           return s_opac_variant<SOpacs...>(s_opac.GetOnDevice());
         },
@@ -73,7 +73,7 @@ class MeanSVariant {
 
   PORTABLE_INLINE_FUNCTION Real
   PlanckMeanTotalScatteringCoefficient(const Real rho, const Real temp) const {
-    return mpark::visit(
+    return PortsOfCall::visit(
         [=](const auto &s_opac) {
           return s_opac.PlanckMeanTotalScatteringCoefficient(rho, temp);
         },
@@ -81,7 +81,7 @@ class MeanSVariant {
   }
   PORTABLE_INLINE_FUNCTION Real RosselandMeanTotalScatteringCoefficient(
       const Real rho, const Real temp) const {
-    return mpark::visit(
+    return PortsOfCall::visit(
         [=](const auto &s_opac) {
           return s_opac.RosselandMeanTotalScatteringCoefficient(rho, temp);
         },
@@ -89,8 +89,8 @@ class MeanSVariant {
   }
 
   inline void Finalize() noexcept {
-    return mpark::visit([](auto &s_opac) { return s_opac.Finalize(); },
-                        s_opac_);
+    return PortsOfCall::visit([](auto &s_opac) { return s_opac.Finalize(); },
+                              s_opac_);
   }
 };
 
diff --git a/singularity-opac/photons/mean_photon_variant.hpp b/singularity-opac/photons/mean_photon_variant.hpp
index a307d0f..7762a85 100644
--- a/singularity-opac/photons/mean_photon_variant.hpp
+++ b/singularity-opac/photons/mean_photon_variant.hpp
@@ -19,11 +19,11 @@
 #include <utility>
 
 #include <ports-of-call/portability.hpp>
+#include <ports-of-call/variant.hpp>
 #include <singularity-opac/base/opac_error.hpp>
 #include <singularity-opac/base/radiation_types.hpp>
 #include <singularity-opac/photons/mean_photon_types.hpp>
 #include <singularity-opac/photons/photon_variant.hpp>
-#include <variant/include/mpark/variant.hpp>
 
 namespace singularity {
 namespace photons {
@@ -61,24 +61,24 @@ class MeanVariant {
           !std::is_same<MeanVariant, typename std::decay<Choice>::type>::value,
           bool>::type = true>
   Choice get() {
-    return mpark::get<Choice>(opac_);
+    return PortsOfCall::get<Choice>(opac_);
   }
 
   MeanVariant GetOnDevice() {
-    return mpark::visit(
+    return PortsOfCall::visit(
         [](auto &opac) { return opac_variant<Opacs...>(opac.GetOnDevice()); },
         opac_);
   }
 
   PORTABLE_INLINE_FUNCTION RuntimePhysicalConstants
   GetRuntimePhysicalConstants() const {
-    return mpark::visit(
+    return PortsOfCall::visit(
         [](auto &opac) { return opac.GetRuntimePhysicalConstants(); }, opac_);
   }
 
   PORTABLE_INLINE_FUNCTION Real
   PlanckMeanAbsorptionCoefficient(const Real rho, const Real temp) const {
-    return mpark::visit(
+    return PortsOfCall::visit(
         [=](const auto &opac) {
           return opac.PlanckMeanAbsorptionCoefficient(rho, temp);
         },
@@ -86,7 +86,7 @@ class MeanVariant {
   }
   PORTABLE_INLINE_FUNCTION Real
   RosselandMeanAbsorptionCoefficient(const Real rho, const Real temp) const {
-    return mpark::visit(
+    return PortsOfCall::visit(
         [=](const auto &opac) {
           return opac.RosselandMeanAbsorptionCoefficient(rho, temp);
         },
@@ -95,8 +95,8 @@ class MeanVariant {
 
   PORTABLE_INLINE_FUNCTION
   Real AbsorptionCoefficient(const Real rho, const Real temp,
-                            const int gmode = Rosseland) const {
-    return mpark::visit(
+                             const int gmode = Rosseland) const {
+    return PortsOfCall::visit(
         [=](const auto &opac) {
           return opac.AbsorptionCoefficient(rho, temp, gmode);
         },
@@ -104,10 +104,9 @@ class MeanVariant {
   }
 
   PORTABLE_INLINE_FUNCTION
-  Real Emissivity(const Real rho, const Real temp,
-                  const int gmode = Rosseland,
+  Real Emissivity(const Real rho, const Real temp, const int gmode = Rosseland,
                   Real *lambda = nullptr) const {
-    return mpark::visit(
+    return PortsOfCall::visit(
         [=](const auto &opac) {
           return opac.Emissivity(rho, temp, gmode, lambda);
         },
@@ -115,12 +114,14 @@ class MeanVariant {
   }
 
   inline void Finalize() noexcept {
-    return mpark::visit([](auto &opac) { return opac.Finalize(); }, opac_);
+    return PortsOfCall::visit([](auto &opac) { return opac.Finalize(); },
+                              opac_);
   }
 
 #ifdef SPINER_USE_HDF
   void Save(const std::string &filename) const {
-    return mpark::visit([=](auto &opac) { return opac.Save(filename); }, opac_);
+    return PortsOfCall::visit([=](auto &opac) { return opac.Save(filename); },
+                              opac_);
   }
 #endif
 };
diff --git a/singularity-opac/photons/photon_s_variant.hpp b/singularity-opac/photons/photon_s_variant.hpp
index 4ae20ce..a7abc2a 100644
--- a/singularity-opac/photons/photon_s_variant.hpp
+++ b/singularity-opac/photons/photon_s_variant.hpp
@@ -19,16 +19,16 @@
 #include <utility>
 
 #include <ports-of-call/portability.hpp>
+#include <ports-of-call/variant.hpp>
 #include <singularity-opac/base/opac_error.hpp>
 #include <singularity-opac/base/radiation_types.hpp>
-#include <variant/include/mpark/variant.hpp>
 
 namespace singularity {
 namespace photons {
 namespace impl {
 
 template <typename... Ts>
-using s_opac_variant = mpark::variant<Ts...>;
+using s_opac_variant = PortsOfCall::variant<Ts...>;
 
 template <typename... S_Opacs>
 class S_Variant {
@@ -62,11 +62,11 @@ class S_Variant {
           !std::is_same<S_Variant, typename std::decay<Choice>::type>::value,
           bool>::type = true>
   Choice get() {
-    return mpark::get<Choice>(s_opac_);
+    return PortsOfCall::get<Choice>(s_opac_);
   }
 
   S_Variant GetOnDevice() {
-    return mpark::visit(
+    return PortsOfCall::visit(
         [](auto &s_opac) {
           return s_opac_variant<S_Opacs...>(s_opac.GetOnDevice());
         },
@@ -79,7 +79,7 @@ class S_Variant {
   PORTABLE_INLINE_FUNCTION Real
   TotalCrossSection(const Real rho, const Real temp, const Real nu,
                     Real *lambda = nullptr) const {
-    return mpark::visit(
+    return PortsOfCall::visit(
         [=](const auto &s_opac) {
           return s_opac.TotalCrossSection(rho, temp, nu, lambda);
         },
@@ -92,7 +92,7 @@ class S_Variant {
   PORTABLE_INLINE_FUNCTION Real
   DifferentialCrossSection(const Real rho, const Real temp, const Real nu,
                            const Real mu, Real *lambda = nullptr) const {
-    return mpark::visit(
+    return PortsOfCall::visit(
         [=](const auto &s_opac) {
           return s_opac.DifferentialCrossSection(rho, temp, nu, mu, lambda);
         },
@@ -105,7 +105,7 @@ class S_Variant {
   PORTABLE_INLINE_FUNCTION Real
   TotalScatteringCoefficient(const Real rho, const Real temp, const Real nu,
                              Real *lambda = nullptr) const {
-    return mpark::visit(
+    return PortsOfCall::visit(
         [=](const auto &s_opac) {
           return s_opac.TotalScatteringCoefficient(rho, temp, nu, lambda);
         },
@@ -114,24 +114,24 @@ class S_Variant {
 
   PORTABLE_INLINE_FUNCTION
   int nlambda() const noexcept {
-    return mpark::visit([](const auto &s_opac) { return s_opac.nlambda(); },
-                        s_opac_);
+    return PortsOfCall::visit(
+        [](const auto &s_opac) { return s_opac.nlambda(); }, s_opac_);
   }
 
   template <typename T>
   PORTABLE_INLINE_FUNCTION bool IsType() const noexcept {
-    return mpark::holds_alternative<T>(s_opac_);
+    return PortsOfCall::holds_alternative<T>(s_opac_);
   }
 
   PORTABLE_INLINE_FUNCTION
   void PrintParams() const noexcept {
-    return mpark::visit([](const auto &s_opac) { return s_opac.PrintParams(); },
-                        s_opac_);
+    return PortsOfCall::visit(
+        [](const auto &s_opac) { return s_opac.PrintParams(); }, s_opac_);
   }
 
   inline void Finalize() noexcept {
-    return mpark::visit([](auto &s_opac) { return s_opac.Finalize(); },
-                        s_opac_);
+    return PortsOfCall::visit([](auto &s_opac) { return s_opac.Finalize(); },
+                              s_opac_);
   }
 };
 
diff --git a/singularity-opac/photons/photon_variant.hpp b/singularity-opac/photons/photon_variant.hpp
index 0b2a1e7..33df636 100644
--- a/singularity-opac/photons/photon_variant.hpp
+++ b/singularity-opac/photons/photon_variant.hpp
@@ -19,16 +19,16 @@
 #include <utility>
 
 #include <ports-of-call/portability.hpp>
+#include <ports-of-call/variant.hpp>
 #include <singularity-opac/base/opac_error.hpp>
 #include <singularity-opac/base/radiation_types.hpp>
-#include <variant/include/mpark/variant.hpp>
 
 namespace singularity {
 namespace photons {
 namespace impl {
 
 template <typename... Ts>
-using opac_variant = mpark::variant<Ts...>;
+using opac_variant = PortsOfCall::variant<Ts...>;
 
 template <typename... Opacs>
 class Variant {
@@ -62,18 +62,18 @@ class Variant {
           !std::is_same<Variant, typename std::decay<Choice>::type>::value,
           bool>::type = true>
   Choice get() {
-    return mpark::get<Choice>(opac_);
+    return PortsOfCall::get<Choice>(opac_);
   }
 
   Variant GetOnDevice() {
-    return mpark::visit(
+    return PortsOfCall::visit(
         [](auto &opac) { return opac_variant<Opacs...>(opac.GetOnDevice()); },
         opac_);
   }
 
   PORTABLE_INLINE_FUNCTION RuntimePhysicalConstants
   GetRuntimePhysicalConstants() const {
-    return mpark::visit(
+    return PortsOfCall::visit(
         [](auto &opac) { return opac.GetRuntimePhysicalConstants(); }, opac_);
   }
 
@@ -83,7 +83,7 @@ class Variant {
   PORTABLE_INLINE_FUNCTION Real
   AbsorptionCoefficient(const Real rho, const Real temp, const Real nu,
                         Real *lambda = nullptr) const {
-    return mpark::visit(
+    return PortsOfCall::visit(
         [=](const auto &opac) {
           return opac.AbsorptionCoefficient(rho, temp, nu, lambda);
         },
@@ -95,7 +95,7 @@ class Variant {
   AbsorptionCoefficient(const Real rho, const Real temp,
                         FrequencyIndexer &nu_bins, DataIndexer &coeffs,
                         const int nbins, Real *lambda = nullptr) const {
-    mpark::visit(
+    PortsOfCall::visit(
         [&](const auto &opac) {
           opac.AbsorptionCoefficient(rho, temp, nu_bins, coeffs, nbins, lambda);
         },
@@ -108,7 +108,7 @@ class Variant {
   PORTABLE_INLINE_FUNCTION Real AngleAveragedAbsorptionCoefficient(
       const Real rho, const Real temp, const Real nu,
       Real *lambda = nullptr) const {
-    return mpark::visit(
+    return PortsOfCall::visit(
         [=](const auto &opac) {
           return opac.AngleAveragedAbsorptionCoefficient(rho, temp, nu, lambda);
         },
@@ -119,7 +119,7 @@ class Variant {
   PORTABLE_INLINE_FUNCTION void AngleAveragedAbsorptionCoefficient(
       const Real rho, const Real temp, FrequencyIndexer &nu_bins,
       DataIndexer &coeffs, const int nbins, Real *lambda = nullptr) const {
-    mpark::visit(
+    PortsOfCall::visit(
         [&](const auto &opac) {
           opac.AngleAveragedAbsorptionCoefficient(rho, temp, nu_bins, coeffs,
                                                   nbins, lambda);
@@ -132,7 +132,7 @@ class Variant {
   PORTABLE_INLINE_FUNCTION Real
   EmissivityPerNuOmega(const Real rho, const Real temp, const Real nu,
                        Real *lambda = nullptr) const {
-    return mpark::visit(
+    return PortsOfCall::visit(
         [=](const auto &opac) {
           return opac.EmissivityPerNuOmega(rho, temp, nu, lambda);
         },
@@ -144,7 +144,7 @@ class Variant {
   EmissivityPerNuOmega(const Real rho, const Real temp,
                        FrequencyIndexer &nu_bins, DataIndexer &coeffs,
                        const int nbins, Real *lambda = nullptr) const {
-    mpark::visit(
+    PortsOfCall::visit(
         [&](const auto &opac) {
           return opac.EmissivityPerNuOmega(rho, temp, nu_bins, coeffs, nbins,
                                            lambda);
@@ -156,7 +156,7 @@ class Variant {
   PORTABLE_INLINE_FUNCTION Real EmissivityPerNu(const Real rho, const Real temp,
                                                 const Real nu,
                                                 Real *lambda = nullptr) const {
-    return mpark::visit(
+    return PortsOfCall::visit(
         [=](const auto &opac) {
           return opac.EmissivityPerNu(rho, temp, nu, lambda);
         },
@@ -168,7 +168,7 @@ class Variant {
   EmissivityPerNu(const Real rho, const Real temp, FrequencyIndexer &nu_bins,
                   DataIndexer &coeffs, const int nbins,
                   Real *lambda = nullptr) const {
-    mpark::visit(
+    PortsOfCall::visit(
         [&](const auto &opac) {
           return opac.EmissivityPerNu(rho, temp, nu_bins, coeffs, nbins,
                                       lambda);
@@ -179,7 +179,7 @@ class Variant {
   // emissivity integrated over angle and frequency
   PORTABLE_INLINE_FUNCTION Real Emissivity(const Real rho, const Real temp,
                                            Real *lambda = nullptr) const {
-    return mpark::visit(
+    return PortsOfCall::visit(
         [=](const auto &opac) { return opac.Emissivity(rho, temp, lambda); },
         opac_);
   }
@@ -188,7 +188,7 @@ class Variant {
   PORTABLE_INLINE_FUNCTION auto NumberEmissivity(const Real rho,
                                                  const Real temp,
                                                  Real *lambda = nullptr) const {
-    return mpark::visit(
+    return PortsOfCall::visit(
         [=](const auto &opac) {
           return opac.NumberEmissivity(rho, temp, lambda);
         },
@@ -198,7 +198,7 @@ class Variant {
   // Specific intensity of thermal distribution
   PORTABLE_INLINE_FUNCTION Real ThermalDistributionOfTNu(
       const Real temp, const Real nu, Real *lambda = nullptr) const {
-    return mpark::visit(
+    return PortsOfCall::visit(
         [=](const auto &opac) {
           return opac.ThermalDistributionOfTNu(temp, nu, lambda);
         },
@@ -208,7 +208,7 @@ class Variant {
   // Temperature derivative of specific intensity of thermal distribution
   PORTABLE_INLINE_FUNCTION Real DThermalDistributionOfTNuDT(
       const Real temp, const Real nu, Real *lambda = nullptr) const {
-    return mpark::visit(
+    return PortsOfCall::visit(
         [=](const auto &opac) {
           return opac.DThermalDistributionOfTNuDT(temp, nu, lambda);
         },
@@ -218,7 +218,7 @@ class Variant {
   // Integral of thermal distribution over frequency and angle
   PORTABLE_INLINE_FUNCTION Real
   ThermalDistributionOfT(const Real temp, Real *lambda = nullptr) const {
-    return mpark::visit(
+    return PortsOfCall::visit(
         [=](const auto &opac) {
           return opac.ThermalDistributionOfT(temp, lambda);
         },
@@ -228,7 +228,7 @@ class Variant {
   // Integral of thermal distribution over energy per frequency and angle
   PORTABLE_INLINE_FUNCTION Real
   ThermalNumberDistributionOfT(const Real temp, Real *lambda = nullptr) const {
-    return mpark::visit(
+    return PortsOfCall::visit(
         [=](const auto &opac) {
           return opac.ThermalNumberDistributionOfT(temp, lambda);
         },
@@ -238,7 +238,7 @@ class Variant {
   // Energy density of thermal distribution
   PORTABLE_INLINE_FUNCTION Real
   EnergyDensityFromTemperature(const Real temp, Real *lambda = nullptr) const {
-    return mpark::visit(
+    return PortsOfCall::visit(
         [=](const auto &opac) {
           return opac.EnergyDensityFromTemperature(temp, lambda);
         },
@@ -248,7 +248,7 @@ class Variant {
   // Temperature of thermal distribution
   PORTABLE_INLINE_FUNCTION Real
   TemperatureFromEnergyDensity(const Real er, Real *lambda = nullptr) const {
-    return mpark::visit(
+    return PortsOfCall::visit(
         [=](const auto &opac) {
           return opac.TemperatureFromEnergyDensity(er, lambda);
         },
@@ -258,7 +258,7 @@ class Variant {
   // Number density of thermal distribution
   PORTABLE_INLINE_FUNCTION Real
   NumberDensityFromTemperature(const Real temp, Real *lambda = nullptr) const {
-    return mpark::visit(
+    return PortsOfCall::visit(
         [=](const auto &opac) {
           return opac.NumberDensityFromTemperature(temp, lambda);
         },
@@ -267,22 +267,22 @@ class Variant {
 
   PORTABLE_INLINE_FUNCTION
   int nlambda() const noexcept {
-    return mpark::visit([](const auto &opac) { return opac.nlambda(); }, opac_);
+    return PortsOfCall::visit([](const auto &opac) { return opac.nlambda(); }, opac_);
   }
 
   template <typename T>
   PORTABLE_INLINE_FUNCTION bool IsType() const noexcept {
-    return mpark::holds_alternative<T>(opac_);
+    return PortsOfCall::holds_alternative<T>(opac_);
   }
 
   PORTABLE_INLINE_FUNCTION
   void PrintParams() const noexcept {
-    return mpark::visit([](const auto &opac) { return opac.PrintParams(); },
+    return PortsOfCall::visit([](const auto &opac) { return opac.PrintParams(); },
                         opac_);
   }
 
   inline void Finalize() noexcept {
-    return mpark::visit([](auto &opac) { return opac.Finalize(); }, opac_);
+    return PortsOfCall::visit([](auto &opac) { return opac.Finalize(); }, opac_);
   }
 };
 
diff --git a/singularity-opac/photons/s_opac_photons.hpp b/singularity-opac/photons/s_opac_photons.hpp
index 45e5562..472ffa8 100644
--- a/singularity-opac/photons/s_opac_photons.hpp
+++ b/singularity-opac/photons/s_opac_photons.hpp
@@ -16,7 +16,7 @@
 #ifndef SINGULARITY_OPAC_PHOTONS_S_OPAC_PHOTONS_
 #define SINGULARITY_OPAC_PHOTONS_S_OPAC_PHOTONS_
 
-#include <variant/include/mpark/variant.hpp>
+#include <ports-of-call/variant.hpp>
 
 #include <singularity-opac/photons/gray_s_opacity_photons.hpp>
 #include <singularity-opac/photons/non_cgs_s_photons.hpp>
diff --git a/utils/cuda_compatibility.patch b/utils/cuda_compatibility.patch
deleted file mode 100644
index d9f10e6..0000000
--- a/utils/cuda_compatibility.patch
+++ /dev/null
@@ -1,142 +0,0 @@
-diff --git a/include/mpark/variant.hpp b/include/mpark/variant.hpp
-index 2fb2ac549..aeef9dcdd 100644
---- a/include/mpark/variant.hpp
-+++ b/include/mpark/variant.hpp
-@@ -1021,9 +1021,15 @@ namespace mpark {
-     friend struct access::recursive_union;                                 \
-   }
- 
-+#ifdef __CUDACC__
-+#define V_GPU_FUNCTION __host__ __device__
-+#else
-+#define V_GPU_FUNCTION
-+#endif
-+
-     MPARK_VARIANT_RECURSIVE_UNION(Trait::TriviallyAvailable,
-                                   ~recursive_union() = default;);
--    MPARK_VARIANT_RECURSIVE_UNION(Trait::Available,
-+    MPARK_VARIANT_RECURSIVE_UNION(Trait::Available,V_GPU_FUNCTION
-                                   ~recursive_union() {});
-     MPARK_VARIANT_RECURSIVE_UNION(Trait::Unavailable,
-                                   ~recursive_union() = delete;);
-@@ -1133,9 +1139,9 @@ namespace mpark {
-         });
- 
-     MPARK_VARIANT_DESTRUCTOR(
--        Trait::Available,
-+        Trait::Available,V_GPU_FUNCTION
-         ~destructor() { destroy(); },
--        inline void destroy() noexcept {
-+        inline constexpr void destroy() noexcept {
-           if (!this->valueless_by_exception()) {
-             visitation::alt::visit_alt(dtor{}, *this);
-           }
-@@ -1169,14 +1175,14 @@ namespace mpark {
- #endif
- 
-       template <std::size_t I, typename T, typename... Args>
--      inline static T &construct_alt(alt<I, T> &a, Args &&... args) {
-+      inline constexpr static T &construct_alt(alt<I, T> &a, Args &&... args) {
-         auto *result = ::new (static_cast<void *>(lib::addressof(a)))
-             alt<I, T>(in_place_t{}, lib::forward<Args>(args)...);
-         return result->value;
-       }
- 
-       template <typename Rhs>
--      inline static void generic_construct(constructor &lhs, Rhs &&rhs) {
-+      inline constexpr static void generic_construct(constructor &lhs, Rhs &&rhs) {
-         lhs.destroy();
-         if (!rhs.valueless_by_exception()) {
-           visitation::alt::visit_alt_at(
-@@ -1223,7 +1229,7 @@ namespace mpark {
- 
-     MPARK_VARIANT_MOVE_CONSTRUCTOR(
-         Trait::Available,
--        move_constructor(move_constructor &&that) noexcept(
-+        constexpr move_constructor(move_constructor &&that) noexcept(
-             lib::all<std::is_nothrow_move_constructible<Ts>::value...>::value)
-             : move_constructor(valueless_t{}) {
-           this->generic_construct(*this, lib::move(that));
-@@ -1256,11 +1262,11 @@ namespace mpark {
-   }
- 
-     MPARK_VARIANT_COPY_CONSTRUCTOR(
--        Trait::TriviallyAvailable,
-+        Trait::TriviallyAvailable,V_GPU_FUNCTION
-         copy_constructor(const copy_constructor &that) = default;);
- 
-     MPARK_VARIANT_COPY_CONSTRUCTOR(
--        Trait::Available,
-+        Trait::Available,V_GPU_FUNCTION
-         copy_constructor(const copy_constructor &that)
-             : copy_constructor(valueless_t{}) {
-           this->generic_construct(*this, that);
-@@ -1281,7 +1287,7 @@ namespace mpark {
-       using super::operator=;
- 
-       template <std::size_t I, typename... Args>
--      inline /* auto & */ auto emplace(Args &&... args)
-+      inline constexpr /* auto & */ auto emplace(Args &&... args)
-           -> decltype(this->construct_alt(access::base::get_alt<I>(*this),
-                                           lib::forward<Args>(args)...)) {
-         this->destroy();
-@@ -1304,7 +1310,7 @@ namespace mpark {
- #endif
- 
-       template <std::size_t I, typename T, typename Arg>
--      inline void assign_alt(alt<I, T> &a, Arg &&arg) {
-+      inline constexpr void assign_alt(alt<I, T> &a, Arg &&arg) {
-         if (this->index() == I) {
- #ifdef _MSC_VER
- #pragma warning(push)
-@@ -1332,7 +1338,7 @@ namespace mpark {
-       }
- 
-       template <typename That>
--      inline void generic_assign(That &&that) {
-+      inline constexpr void generic_assign(That &&that) {
-         if (this->valueless_by_exception() && that.valueless_by_exception()) {
-           // do nothing.
-         } else if (that.valueless_by_exception()) {
-@@ -1376,11 +1382,11 @@ namespace mpark {
-   }
- 
-     MPARK_VARIANT_MOVE_ASSIGNMENT(
--        Trait::TriviallyAvailable,
-+        Trait::TriviallyAvailable,V_GPU_FUNCTION
-         move_assignment &operator=(move_assignment &&that) = default;);
- 
-     MPARK_VARIANT_MOVE_ASSIGNMENT(
--        Trait::Available,
-+        Trait::Available,V_GPU_FUNCTION
-         move_assignment &
-         operator=(move_assignment &&that) noexcept(
-             lib::all<(std::is_nothrow_move_constructible<Ts>::value &&
-@@ -1394,6 +1400,7 @@ namespace mpark {
-         move_assignment &operator=(move_assignment &&) = delete;);
- 
- #undef MPARK_VARIANT_MOVE_ASSIGNMENT
-+#undef V_GPU_FUNCTION
- 
-     template <typename Traits, Trait = Traits::copy_assignable_trait>
-     class copy_assignment;
-@@ -1421,7 +1428,7 @@ namespace mpark {
- 
-     MPARK_VARIANT_COPY_ASSIGNMENT(
-         Trait::Available,
--        copy_assignment &operator=(const copy_assignment &that) {
-+        constexpr copy_assignment &operator=(const copy_assignment &that) {
-           this->generic_assign(that);
-           return *this;
-         });
-@@ -1694,8 +1701,8 @@ namespace mpark {
- 
-     ~variant() = default;
- 
--    variant &operator=(const variant &) = default;
--    variant &operator=(variant &&) = default;
-+    /*V_GPU_FUNCTION*/constexpr variant &operator=(const variant &) = default;
-+    /*V_GPU_FUNCTION*/constexpr variant &operator=(variant &&) = default;
- 
-     template <typename Arg,
-               lib::enable_if_t<!std::is_same<lib::decay_t<Arg>, variant>::value,
diff --git a/utils/ports-of-call b/utils/ports-of-call
index 58ce118..a284b6b 160000
--- a/utils/ports-of-call
+++ b/utils/ports-of-call
@@ -1 +1 @@
-Subproject commit 58ce1181b2d835bd32673ad70550c9130381f91b
+Subproject commit a284b6b2d42e70afeb99babb23522c869200d6ea
diff --git a/utils/spiner b/utils/spiner
index aa2e15d..e40c75d 160000
--- a/utils/spiner
+++ b/utils/spiner
@@ -1 +1 @@
-Subproject commit aa2e15dd1eaae1d5c69655c72df41d9d36111107
+Subproject commit e40c75df722aee786a6bf18ea88650630b0aed35
diff --git a/utils/variant b/utils/variant
deleted file mode 160000
index 23cb94f..0000000
--- a/utils/variant
+++ /dev/null
@@ -1 +0,0 @@
-Subproject commit 23cb94f027d4ef33bf48133acc2695c7e5c6f1e7