diff --git a/CMakeLists.txt b/CMakeLists.txt
index acca84370..e412eb785 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,4 +1,5 @@
 cmake_minimum_required(VERSION 3.15)
+set(CMAKE_CXX_COMPILER_WORKS 1)
 project(simsoptpp)
 #set(CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/cmake ${CMAKE_MODULE_PATH})
 
@@ -41,7 +42,7 @@ include_directories(${CMAKE_CURRENT_BINARY_DIR})
 include(CheckCXXCompilerFlag)
 IF(DEFINED ENV{CI})
     message(STATUS "CI environment detected. Set compilation flags targetting Westmere microarch.")
-    set(CMAKE_CXX_FLAGS "-O3 -march=westmere")
+    set(CMAKE_CXX_FLAGS "-O3 -march=znver3")
 elseif(DEFINED ENV{CONDA_BUILD})
     message(STATUS "conda build environment detected. Let conda set compilation flags accordingly.")
     # set(CMAKE_CXX_FLAGS "-O3 -march=ivybridge -mfma -ffp-contract=fast")
diff --git a/src/simsoptpp/python_tracing.cpp b/src/simsoptpp/python_tracing.cpp
index c862bbcd7..17c468846 100644
--- a/src/simsoptpp/python_tracing.cpp
+++ b/src/simsoptpp/python_tracing.cpp
@@ -46,15 +46,16 @@ void init_tracing(py::module_ &m){
         py::arg("vacuum"),
         py::arg("noK"),
         py::arg("solveSympl")=false,
+        py::arg("GPU")=false,
         py::arg("zetas")=vector<double>{},
         py::arg("omegas")=vector<double>{},
         py::arg("stopping_criteria")=vector<shared_ptr<StoppingCriterion>>{},
-        py::arg("vpars")=vector<double>{},
         py::arg("phis_stop")=false,
         py::arg("vpars_stop")=false,
         py::arg("forget_exact_path")=false,
         py::arg("axis")=0,
-        py::arg("predictor_step")=true
+        py::arg("predictor_step")=true,
+        py::arg("vpars")=vector<double>{}
         );
 
     m.def("particle_guiding_center_boozer_perturbed_tracing", &particle_guiding_center_boozer_perturbed_tracing<xt::pytensor>,
@@ -73,7 +74,6 @@ void init_tracing(py::module_ &m){
         py::arg("zetas")=vector<double>{},
         py::arg("omegas")=vector<double>{},
         py::arg("stopping_criteria")=vector<shared_ptr<StoppingCriterion>>{},
-        py::arg("vpars")=vector<double>{},
         py::arg("phis_stop")=false,
         py::arg("vpars_stop")=false,
         py::arg("Phihat")=0,
@@ -82,8 +82,9 @@ void init_tracing(py::module_ &m){
         py::arg("Phin")=0,
         py::arg("phase")=0,
         py::arg("forget_exact_path")=false,
-        py::arg("axis")=0
+        py::arg("axis")=0,
+        py::arg("vpars")=vector<double>{}
     );
-    
+
     m.def("get_phi", &get_phi);
 }
diff --git a/src/simsoptpp/tracing.cpp b/src/simsoptpp/tracing.cpp
index b0dd43562..287f4e54b 100644
--- a/src/simsoptpp/tracing.cpp
+++ b/src/simsoptpp/tracing.cpp
@@ -1196,9 +1196,9 @@ particle_guiding_center_boozer_perturbed_tracing(
         shared_ptr<BoozerMagneticField<T>> field, array<double, 3> stz_init,
         double m, double q, double vtotal, double vtang, double mu, double tmax, double abstol, double reltol,
         bool vacuum, bool noK, vector<double> zetas, vector<double> omegas,
-        vector<shared_ptr<StoppingCriterion>> stopping_criteria, vector<double> vpars,
+        vector<shared_ptr<StoppingCriterion>> stopping_criteria,
         bool phis_stop, bool vpars_stop, double Phihat, double omega, int Phim,
-        int Phin, double phase, bool forget_exact_path, int axis)
+        int Phin, double phase, bool forget_exact_path, int axis, vector<double> vpars)
 {
     typename BoozerMagneticField<T>::Tensor2 stz({{stz_init[0], stz_init[1], stz_init[2]}});
     field->set_points(stz);
@@ -1234,9 +1234,9 @@ tuple<vector<array<double, 5>>, vector<array<double, 6>>>
 particle_guiding_center_boozer_tracing(
         shared_ptr<BoozerMagneticField<T>> field, array<double, 3> stz_init,
         double m, double q, double vtotal, double vtang, double tmax, double dt, double abstol, double reltol, double roottol,
-        bool vacuum, bool noK, bool solveSympl, vector<double> zetas, vector<double> omegas,
-        vector<shared_ptr<StoppingCriterion>> stopping_criteria, vector<double> vpars,
-        bool phis_stop, bool vpars_stop, bool forget_exact_path, int axis, bool predictor_step)
+        bool vacuum, bool noK, bool GPU, bool solveSympl, vector<double> zetas, vector<double> omegas,
+        vector<shared_ptr<StoppingCriterion>> stopping_criteria,
+        bool phis_stop, bool vpars_stop, bool forget_exact_path, int axis, bool predictor_step, vector<double> vpars)
 {
     typename BoozerMagneticField<T>::Tensor2 stz({{stz_init[0], stz_init[1], stz_init[2]}});
     field->set_points(stz);
@@ -1282,6 +1282,8 @@ particle_guiding_center_boozer_tracing(
         auto f = SymplField<T>(field, m, q, mu);
         return solve_sympl(f, y, tmax, dt, roottol, zetas, omegas, stopping_criteria,
         vpars, phis_stop, vpars_stop, forget_exact_path, predictor_step);
+    } else if (GPU) {
+        return;
     } else {
         if (vacuum) {
           auto rhs_class = GuidingCenterVacuumBoozerRHS<T>(field, m, q, mu, axis);
@@ -1299,19 +1301,21 @@ particle_guiding_center_boozer_tracing(
     }
 }
 
-template
-tuple<vector<array<double, 6>>, vector<array<double, 7>>> particle_guiding_center_boozer_perturbed_tracing<xt::pytensor>(
+template<template<class, std::size_t, xt::layout_type> class T>
+tuple<vector<array<double, 6>>, vector<array<double, 7>>>
+particle_guiding_center_boozer_perturbed_tracing<xt::pytensor>(
         shared_ptr<BoozerMagneticField<xt::pytensor>> field, array<double, 3> stz_init,
         double m, double q, double vtotal, double vtang, double mu, double tmax, double abstol, double reltol,
         bool vacuum, bool noK, vector<double> zetas, vector<double> omegas,
         vector<shared_ptr<StoppingCriterion>> stopping_criteria,
-        vector<double> vpars={}, bool phis_stop, bool vpars_stop, double Phihat,
-        double omega, int Phim, int Phin, double phase, bool forget_exact_path, int axis);
+        bool phis_stop, bool vpars_stop, double Phihat,
+        double omega, int Phim, int Phin, double phase, bool forget_exact_path, int axis, vector<double> vpars);
 
-template
-tuple<vector<array<double, 5>>, vector<array<double, 6>>> particle_guiding_center_boozer_tracing<xt::pytensor>(
+template<template<class, std::size_t, xt::layout_type> class T>
+tuple<vector<array<double, 5>>, vector<array<double, 6>>>
+particle_guiding_center_boozer_tracing<xt::pytensor>(
         shared_ptr<BoozerMagneticField<xt::pytensor>> field, array<double, 3> stz_init,
         double m, double q, double vtotal, double vtang, double tmax, double dt, double abstol, double reltol, double roottol,
-        bool vacuum, bool noK, bool solveSympl, vector<double> zetas, vector<double> omegas,
+        bool vacuum, bool noK, bool GPU, bool solveSympl, vector<double> zetas, vector<double> omegas,
         vector<shared_ptr<StoppingCriterion>> stopping_criteria,
-        vector<double> vpars={}, bool phis_stop, bool vpars_stop, bool forget_exact_path, int axis, bool predictor_step);
\ No newline at end of file
+        bool forget_exact_path, int axis, bool predictor_step, bool phis_stop, bool vpars_stop, vector<double> vpars);
diff --git a/src/simsoptpp/tracing.h b/src/simsoptpp/tracing.h
index 677536157..546ab04df 100644
--- a/src/simsoptpp/tracing.h
+++ b/src/simsoptpp/tracing.h
@@ -117,51 +117,52 @@ class StepSizeStoppingCriterion : public StoppingCriterion {
     private:
         int min_dt;
     public:
-        StepSizeStoppingCriterion(int min_dt) : min_dt(min_dt) {};
-        bool operator()(int iter, double dt, double t, double x, double y, double z, double vpar=0) override {
-        // bool operator()(int iter, double t, Array& y) override {
-            return dt<min_dt;
-        };
-};
+			  StepSizeStoppingCriterion(int min_dt) : min_dt(min_dt) {};
+			  bool operator()(int iter, double dt, double t, double x, double y, double z, double vpar=0) override {
+			  // bool operator()(int iter, double t, Array& y) override {
+					return dt<min_dt;
+			  };
+	};
 
-template<class Array>
-class LevelsetStoppingCriterion : public StoppingCriterion {
-    private:
-        shared_ptr<RegularGridInterpolant3D<Array>> levelset;
-    public:
-        LevelsetStoppingCriterion(shared_ptr<RegularGridInterpolant3D<Array>> levelset) : levelset(levelset) { };
-        // bool operator()(int iter, double t, Array2& y) override {
-        bool operator()(int iter, double dt, double t, double x, double y, double z, double vpar=0) override {
-            double r = std::sqrt(x*x + y*y);
-            // double r = std::sqrt(y[0]*y[0] + y[1]*y[1]);
-            double phi = std::atan2(y, x);
-            // double phi = std::atan2(y[1],y[0]);
-            if(phi < 0)
-                phi += 2*M_PI;
-            double f = levelset->evaluate(r, phi, z)[0];
-            // double f = levelset->evaluate(r, phi, y[2])[0];
-            //fmt::print("Levelset at xyz=({}, {}, {}), rphiz=({}, {}, {}), f={}\n", x, y, z, r, phi, z, f);
-            return f<0;
-        };
-};
+	template<class Array>
+	class LevelsetStoppingCriterion : public StoppingCriterion {
+		 private:
+			  shared_ptr<RegularGridInterpolant3D<Array>> levelset;
+		 public:
+			  LevelsetStoppingCriterion(shared_ptr<RegularGridInterpolant3D<Array>> levelset) : levelset(levelset) { };
+			  // bool operator()(int iter, double t, Array2& y) override {
+			  bool operator()(int iter, double dt, double t, double x, double y, double z, double vpar=0) override {
+					double r = std::sqrt(x*x + y*y);
+					// double r = std::sqrt(y[0]*y[0] + y[1]*y[1]);
+					double phi = std::atan2(y, x);
+					// double phi = std::atan2(y[1],y[0]);
+					if(phi < 0)
+						 phi += 2*M_PI;
+					double f = levelset->evaluate(r, phi, z)[0];
+					// double f = levelset->evaluate(r, phi, y[2])[0];
+					//fmt::print("Levelset at xyz=({}, {}, {}), rphiz=({}, {}, {}), f={}\n", x, y, z, r, phi, z, f);
+					return f<0;
+			  };
+	};
 
-template<template<class, std::size_t, xt::layout_type> class T>
-tuple<vector<array<double, 6>>, vector<array<double, 7>>>
-particle_guiding_center_boozer_perturbed_tracing(
-        shared_ptr<BoozerMagneticField<T>> field, array<double, 3> stz_init,
-        double m, double q, double vtotal, double vtang, double mu, double tmax, double abstol, double reltol,
-        bool vacuum, bool noK, vector<double> zetas, vector<double> omegas,
-        vector<shared_ptr<StoppingCriterion>> stopping_criteria, vector<double> vpars,
-        bool zetas_stop=false, bool vpars_stop=false,
-        double alphahat=0, double omega=0, int alpham=0, int alphan=0, double phase=0,
-        bool forget_exact_path=false, int axis=0);
+	template<template<class, std::size_t, xt::layout_type> class T>
+	tuple<vector<array<double, 6>>, vector<array<double, 7>>>
+	particle_guiding_center_boozer_perturbed_tracing(
+			  shared_ptr<BoozerMagneticField<T>> field, array<double, 3> stz_init,
+			  double m, double q, double vtotal, double vtang, double mu, double tmax, double abstol, double reltol,
+			  bool vacuum, bool noK, vector<double> zetas, vector<double> omegas,
+			  vector<shared_ptr<StoppingCriterion>> stopping_criteria,
+			  bool zetas_stop=false, bool vpars_stop=false,
+			  double alphahat=0, double omega=0, int alpham=0, int alphan=0, double phase=0,
+			  bool forget_exact_path=false, int axis=0,  vector<double> vpars={});
 
-template<template<class, std::size_t, xt::layout_type> class T>
-tuple<vector<array<double, 5>>, vector<array<double, 6>>>
-particle_guiding_center_boozer_tracing(
-        shared_ptr<BoozerMagneticField<T>> field, array<double, 3> stz_init,
+	template<template<class, std::size_t, xt::layout_type> class T>
+	tuple<vector<array<double, 5>>, vector<array<double, 6>>>
+	particle_guiding_center_boozer_tracing(
+			  shared_ptr<BoozerMagneticField<T>> field, array<double, 3> stz_init,
         double m, double q, double vtotal, double vtang, double tmax, double dt, double abstol, double reltol, double roottol,
-        bool vacuum, bool noK, bool solveSympl, vector<double> zetas, vector<double> omegas,
+        bool vacuum, bool noK, bool GPU, bool solveSympl, vector<double> zetas, vector<double> omegas,
         vector<shared_ptr<StoppingCriterion>> stopping_criteria,
-        vector<double> vpars, bool zetas_stop=false, bool vpars_stop=false,
-        bool forget_exact_path=false, int axis=0, bool predictor_step=true);
+        bool forget_exact_path=false, int axis=0, bool predictor_step=true,
+        bool zetas_stop=false, bool vpars_stop=false, vector<double> vpars={});
+