diff --git a/src/main.cpp b/src/main.cpp
index b6f9446..d41cf81 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -47,7 +47,11 @@ int main(int argc, char *argv[]) {
   ReducedModelOptimizer::xRange beamDimensionsRatio{"bOverh", 0.7, 1.3};
   ReducedModelOptimizer::xRange beamE{"E", 0.1, 1.9};
   // Test set of full patterns
-  std::string fullPatternsTestSetDirectory = "../TestSet";
+  std::string fullPatternsTestSetDirectory = "TestSet";
+    if (!std::filesystem::exists(std::filesystem::path(fullPatternsTestSetDirectory))) {
+        std::cerr << "Full pattern directory does not exist:" << fullPatternsTestSetDirectory << std::endl;
+        return 1;
+    }
   //      "/home/iason/Documents/PhD/Research/Approximating shapes with flat "
   //      "patterns/Pattern_enumerator/Results/1v_0v_2e_1e_1c_6fan/3/Valid";
   std::vector<std::pair<FlatPattern *, FlatPattern *>> patternPairs;
@@ -74,7 +78,7 @@ int main(int argc, char *argv[]) {
   //  for (double rangeOffset = 0.15; rangeOffset <= 0.95; rangeOffset += 0.05)
   //  {
   ReducedModelOptimizer::Settings settings;
-  for (settings.maxSimulations = 100; settings.maxSimulations < 3500;
+  for (settings.maxSimulations = 100; settings.maxSimulations < 5000;
        settings.maxSimulations += 100) {
     std::string xRangesString = beamWidth.toString() + " " +
                                 beamDimensionsRatio.toString() + " " +
@@ -96,7 +100,7 @@ int main(int argc, char *argv[]) {
         resultsPerPattern(patternPairs.size());
     auto start = std::chrono::high_resolution_clock::now();
 
-    //#pragma omp parallel for
+    #pragma omp parallel for
     for (int patternPairIndex = 0; patternPairIndex < patternPairs.size();
          patternPairIndex++) {
       // const auto filepathString = filepath.string();
@@ -153,7 +157,11 @@ int main(int argc, char *argv[]) {
     //      resultsPerPattern[patternPairIndex].second.save(saveToPath);
     //    }
 
-    csvfile statistics(std::filesystem::path("../OptimizationResults")
+    if (!std::filesystem::exists(std::filesystem::path("OptimizationResults/"))) {
+        std::filesystem::create_directory(std::filesystem::path("OptimizationResults"));
+    }
+
+    csvfile statistics(std::filesystem::path("OptimizationResults")
                            .append("statistics.csv")
                            .string(),
                        false);
diff --git a/src/reducedmodeloptimizer.cpp b/src/reducedmodeloptimizer.cpp
index ab5025f..ba914f2 100644
--- a/src/reducedmodeloptimizer.cpp
+++ b/src/reducedmodeloptimizer.cpp
@@ -32,8 +32,22 @@ struct GlobalOptimizationVariables {
   std::vector<std::vector<double>> failedSimulationsXRatio;
   int numOfSimulationCrashes{false};
   int numberOfFunctionCalls{0};
-} global;
+} ;
+
+//static GlobalOptimizationVariables global;
+
+const static int MAX_THREAD = 64;
+
+struct MY_TLS_ITEM
+{
+    std::map<int, int> theMap;
+    char padding[64 - sizeof(theMap)];
+};
+
+__declspec(align(64)) GlobalOptimizationVariables tls[MAX_THREAD];
+
 //#pragma omp threadprivate(global)
+
 // struct OptimizationCallback {
 //  double operator()(const size_t &iterations, const Eigen::VectorXd &x,
 //                    const double &fval, Eigen::VectorXd &gradient) const {
@@ -115,6 +129,7 @@ double ReducedModelOptimizer::computeError(
     const SimulationResults &reducedPatternResults,
     const Eigen::MatrixX3d &optimalReducedPatternDisplacements) {
   double error = 0;
+  auto& global = tls[omp_get_thread_num()];
   for (const auto reducedFullViPair : global.g_reducedToFullViMap) {
     VertexIndex reducedModelVi = reducedFullViPair.first;
     //    const auto pos =
@@ -141,6 +156,7 @@ double ReducedModelOptimizer::computeError(
 }
 
 void updateMesh(long n, const double *x) {
+  auto& global = tls[omp_get_thread_num()];
   std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh =
       global
           .g_reducedPatternSimulationJob[global.g_simulationScenarioIndices[0]]
@@ -214,6 +230,7 @@ double ReducedModelOptimizer::objective(double x0, double x1, double x2,
 }
 
 double ReducedModelOptimizer::objective(long n, const double *x) {
+  auto& global = tls[omp_get_thread_num()];
   //  std::cout.precision(17);
 
   //  for (size_t parameterIndex = 0; parameterIndex < n; parameterIndex++) {
@@ -356,6 +373,7 @@ void ReducedModelOptimizer::computeMaps(
   //  std::endl;
 
   // Save excluded edges
+  auto& global = tls[omp_get_thread_num()];
   global.g_reducedPatternExludedEdges.clear();
   const size_t fanSize = 6;
   const size_t reducedBaseTriangleNumberOfEdges = reducedPattern.EN();
@@ -491,6 +509,7 @@ void ReducedModelOptimizer::initializePatterns(
   FlatPattern copyReducedPattern;
   copyFullPattern.copy(fullPattern);
   copyReducedPattern.copy(reducedPattern);
+  auto& global = tls[omp_get_thread_num()];
   global.g_optimizeInnerHexagonSize = copyReducedPattern.EN() == 2;
   if (global.g_optimizeInnerHexagonSize) {
     const double h = copyReducedPattern.getBaseTriangleHeight();
@@ -521,6 +540,7 @@ void ReducedModelOptimizer::initializePatterns(
 
 void ReducedModelOptimizer::initializeOptimizationParameters(
     const std::shared_ptr<SimulationMesh> &mesh) {
+  auto& global = tls[omp_get_thread_num()];
   const int numberOfOptimizationParameters = 3;
   global.g_initialParameters.resize(global.g_optimizeInnerHexagonSize
                                         ? numberOfOptimizationParameters + 1
@@ -613,6 +633,7 @@ void ReducedModelOptimizer::computeDesiredReducedModelDisplacements(
 ReducedModelOptimizer::Results ReducedModelOptimizer::runOptimization(
     const Settings &settings,
     double (*pObjectiveFunction)(long, const double *)) {
+  auto& global = tls[omp_get_thread_num()];
 
   global.gObjectiveValueHistory.clear();
 
@@ -981,6 +1002,7 @@ void ReducedModelOptimizer::visualizeResults(
         simulator.executeSimulation(pFullPatternSimulationJob);
     fullModelResults.registerForDrawing();
     fullModelResults.saveDeformedModel();
+  auto& global = tls[omp_get_thread_num()];
     const std::shared_ptr<SimulationJob> &pReducedPatternSimulationJob =
         global.g_reducedPatternSimulationJob[simulationScenarioIndex];
     SimulationResults reducedModelResults =
@@ -1011,6 +1033,7 @@ void ReducedModelOptimizer::visualizeResults(
 ReducedModelOptimizer::Results ReducedModelOptimizer::optimize(
     const Settings &xRanges,
     const std::vector<SimulationScenario> &simulationScenarios) {
+  auto& global = tls[omp_get_thread_num()];
 
   global.g_simulationScenarioIndices = simulationScenarios;
   if (global.g_simulationScenarioIndices.empty()) {